Assigning displacement on nodes (when using condense)

[federicodalinger]

Hi @kinnala, I hope you are well.

I'm continuing with my plate deflection project (using some stuff from your Ex 2), and am now adding some stiffeners (like braces).

When using the condense function, I know it's possible to replace nodes with displacement equal to zero, for others values.

Is there a way to do such a thing by assigning the value by his node?
I have the rectangle as a plate, and the 2 horizontal braces (all components are welded), and I have calculated the deflection for each brace, I mean, for each brace node position. So I already have a list of the global nodes (belonging to those braces) with their respective displacements.

Thanks in advance and I wish you a lovely weekend.

Cheers,
Federico

[kinnala]

Something like:

x = basis.zeros()
D1 = basis.get_dofs('brace1').all('u')
x[D1] = [...]  # your values go here in the correct order
# the order is the same as in basis.doflocs[:, D1]
x = solve(*condense(K, f, x=x, D=D1))

If you have multiple braces then you can do D1, D2, etc. and concatenate them: D=np.concatenate((D1, D2)).

[federicodalinger]

Hi @kinnala , thanks so much for your answer and help. I figured out how to solve the displacement field with your last response (which are the normal deflections along the plate / brace).
However, when plotting the tension fields (sigma_xx, sigma_yy, tau_xy, Vm) I have excessive values on my results, but, when picking the values, from the plot, the results are accurate.

Please, see my Von Mises Stresses plot (with the selected picked point value at the bottom) versus the FEA software I am using to check it (SW simulation).

What do you think about this discrepancy between the plots? Any suggestions?

Thanks in advance mate, and I wish you a good week.

Kind regards,
Federico

[kinnala]

Argyris element (is that what you are using?) is somewhat challenging to use correctly. It has 21 DOFs per element including the derivatives along the x and y axes. If you have fixed the deflection along the x axis, it means that you must also know and, therefore, be able to fix the derivative along the x axis. I believe that you have not fixed all the DOFs you know on the braces and this is visible in the stress.

[kinnala]

Perhaps using something like ElementTriMorley is easier in your case?

[federicodalinger]

Hi @kinnala , thanks a lot for your quick answers.

I am using ElementTriMorley:

Also, just to add some results, please see the displacement plots (normal direction to the plate) versus the one obtained from the code:

I can assume that the normal displacements ("z" axis) are okay, which are very close between both plots.

The boundary conditions are defined as:

As I am coding the script with the idea of including 2 settings (without or with braces on the plate), I generate a new mesh variable that includes the braces, over the mesh for the braceless setting (called "M"):

(The function list_of_faces_and_nodes_indexes.... gives you the close faces and belonging node indexes (for certain x and y positions in the space). For this case I only pick the faces between the start and end point of that brace)

Later, they ("left", "top", "right", "bottom", "horiz_1") are restricted as:

Then, I calculate the displacement solutions as:

The next step I do is to define a variable "u" that stores the interpolation results of "w" into the basis I'm using, and then I use that "u" to obtain the stress fields (sigma_xx, sigma_yy, tau_xy, sigma_vm), using the constitutive relationships.

I hope this explanation of what I've done helps.

Really, thanks a lot @kinnala for your effort with all my queries.

Kind regards,
Federico

[kinnala]

I do observe some small variations also in the displacement plot near y=0.06. Can you double check that your ordering in for index_node_brace in D1 is correct? I do not see you imposing any particular order, unless it is done inside the function calculate_list_displacement_brace.

[federicodalinger]

Hi @kinnala.
Please see the print(D1) to see the order of the node indexes:

Besides, with the function "calculate_list_displacement_brace" I obtain a list of displacements, using the list from D1. So, I obtain a list of displacements and the index position of D1 corresponds to the index position of "calculate_list_displacement_brace". Please tell me if this is what you mean with having the same order.

Thanks in advance and wish you a good day,

Cheers,
Federico

[kinnala]

It's hard to say for sure without debugging the entire code. Here I tried making something similar in a very minimal fashion, by setting the brace displacement from sin(pi * x):

from skfem import *
from skfem.helpers import *
from skfem.models import unit_load
import numpy as np


m = (MeshTri.init_tensor(np.linspace(0, 1, 2 * 20),
                         np.linspace(0, .5, 2 * 10 + 1))
     .with_boundaries({
         'horiz_1': lambda x: np.isclose(x[1], 3 * .5 / 10.),
     }, boundaries_only=False))


@BilinearForm
def bilinf(u, v, w):
    d = .1
    E = 200e9
    nu = .3

    def C(T):
        return E / (1 + nu) * (T + nu / (1 - nu) * eye(trace(T), 2))

    return d**3 / 12.0 * ddot(C(dd(u)), dd(v))


basis = Basis(m, ElementTriMorley())

K = bilinf.assemble(basis)
f = unit_load.assemble(basis)

x = basis.zeros()
D1 = basis.get_dofs('horiz_1').all('u')
x[D1] = np.sin(np.pi * basis.doflocs[0, D1])

D = np.concatenate((
    D1,
    basis.get_dofs({'left', 'right', 'bottom'}).all(),
))

y = solve(*condense(K, f, x=x, D=D))

m.draw(boundaries=True).show()
basis.plot(y, shading='gouraud', nrefs=2).show()

I do not observe anything strange in the displacement field in this case.

[kinnala]

Here I made a small modification to draw the yy-stress as well:

from skfem import *
from skfem.helpers import *
from skfem.models import unit_load
import numpy as np


m = (MeshTri.init_tensor(np.linspace(0, 1, 2 * 20),
                         np.linspace(0, .5, 2 * 10 + 1))
     .with_boundaries({
         'horiz_1': lambda x: np.isclose(x[1], 3 * .5 / 10.),
     }, boundaries_only=False))


d = .1
E = 200e9
nu = .3
def C(T):
    return E / (1 + nu) * (T + nu / (1 - nu) * eye(trace(T), 2))

@BilinearForm
def bilinf(u, v, w):

    return d**3 / 12.0 * ddot(C(dd(u)), dd(v))


basis = Basis(m, ElementTriMorley())

K = bilinf.assemble(basis)
f = unit_load.assemble(basis)

x = basis.zeros()
D1 = basis.get_dofs('horiz_1').all('u')
x[D1] = np.sin(np.pi * basis.doflocs[0, D1])

D = np.concatenate((
    D1,
    basis.get_dofs({'left', 'right', 'bottom'}).all(),
))

y = solve(*condense(K, f, x=x, D=D))

m.draw(boundaries=True).show()
basis.plot(y, shading='gouraud', nrefs=2).show()

basis0 = basis.with_element(ElementTriP0())
syy = basis0.project(C(dd(basis.interpolate(y)))[1, 1])
basis0.plot(syy, shading='gouraud').show()

It looks fine.