Ex 2 - Finding the Von Mises Stress

[federicodalinger]

Hi @kinnala and community, I hope you are well :)

First, I would like to say "thanks a lot" for this library you have shared with the world. Besides, it's really nice and easier to install/understand (in comparison with Fenics, for example). I started to use your library to understand and learn more about FEA.

I am writing a Python script that calculates the deflections of a plate, with certain dimensions, material, and hydrostatic load (non-uniform load). I have found the correct deflections and validated them on Ansys (it's very accurate, with less than 1% error).

However, I am facing some issues when calculating the Von Mises Stress...

In example 4 this is calculated like this:

for i in [0, 1]:
for j in [0, 1]:
s[i, j] = dg.project(C(sym_grad(x))[i, j])

vonmises = np.sqrt(.5 * ((s[0, 0] - s[1, 1]) ** 2 +
(s[1, 1] - s[2, 2]) ** 2 +
(s[2, 2] - s[0, 0]) ** 2 +
6. * s[0, 1] ** 2))

But my first problem starts with "sym_grad(x)" giving me this error:

Traceback (most recent call last):
File "G:\kirchoff_plate_solver\kirchoff_hydro_plate.py", line 178, in
s2[i, j] = dg2.project(C(sym_grad(u2))[i, j])
File "C:\Users\fede_\AppData\Local\Programs\Python\Python39\lib\site-packages\skfem\helpers.py", line 74, in sym_grad
return .5 * (u.grad + transpose(u.grad))
ValueError: operands could not be broadcast together with shapes (2,4096,6) (4096,2,6)

Later on, if I manually do the ".5 * (u.grad + transpose(u.grad))", using NumPy for the transpose of u.grad (otherwise it doesn't work, I don't know why), I have the shape of (2, 4096, 6), for sym_grad(x).

Finally, I have the error in "dg.project(C(sym_grad(x))[i, j])" as:
File "G:\kirchoff_plate_solver\kirchoff_hydro_plate.py", line 170, in
C(result)
File "C:\Users\fede_\AppData\Local\Programs\Python\Python39\lib\site-packages\skfem\models\elasticity.py", line 39, in C
return 2. * Mu * T + Lambda * eye(trace(T), T.shape[0])
File "C:\Users\fede_\AppData\Local\Programs\Python\Python39\lib\site-packages\skfem\helpers.py", line 140, in trace
return np.einsum('ii...', T)
File "C:\Users\fede_\AppData\Local\Programs\Python\Python39\lib\site-packages\numpy\core\einsumfunc.py", line 1371, in einsum
return c_einsum(*operands, **kwargs)
ValueError: dimensions in operand 0 for collapsing index 'i' don't match (2 != 4096).

Really, I don't know what to do to solve this and plot the stresses. Please, Could you guide me on how to calculate it using the displacement of the plate?

Thanks in advance,
Federico

[kinnala]

I believe the issue is the following. The formulation in ex04 is based on plane strain solid and not Kirchhoff plate. For this reason, the equation to calculate strain and stress are different. I don't right now have time to derive the corresponding equations for Kirchhoff plate, but I will try to find some time in the near future.

[federicodalinger]

Thanks! I look forward to seeing it :)
I didn't realize it was based on plane strain solid. I will try to search more about on how to do it using a Kirchoff plate.

Thanks again,
Federico

[kinnala]

Here is a variant of ex02 which, I believe, calculates the different components of stress:

from skfem import *
from skfem.models.poisson import unit_load
from skfem.helpers import dd, ddot, trace, eye
import numpy as np

m = (
    MeshTri.init_symmetric()
    .refined(3)
    .with_boundaries(
        {
            "left": lambda x: x[0] == 0,
            "right": lambda x: x[0] == 1,
            "top": lambda x: x[1] == 1,
        }
    )
)

e = ElementTriMorley()
ib = Basis(m, e)

d = 0.1
E = 200e9
nu = 0.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))


K = asm(bilinf, ib)
f = 1e6 * asm(unit_load, ib)

D = np.hstack([ib.get_dofs("left"), ib.get_dofs({"right", "top"}).all("u")])

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


from skfem.visuals.matplotlib import draw, plot
ax = draw(m)
plot(ib,
     x,
     ax=ax,
     shading='gouraud',
     colorbar=True,
     nrefs=2).show()


u = ib.interpolate(x)
stress = C(dd(u))
p1_ib = ib.with_element(ElementTriP1())

for i in [0, 1]:
    for j in [0, 1]:
        plot(p1_ib,
             p1_ib.project(stress[i, j]),
             colorbar=True,
             nrefs=2).show()

[federicodalinger]

Thanks so much mate. I just returned to this little project. I only needed to multiply it by the 0.5*thickness to measure well the tensions :)

Kind regards