Problem with computing normal stress in 3d solid mechanics example

General
1

Hello

In this example 3D Solid Mechanics — NGS-Py 6.2.2404-2-ge6a630133 documentation (ngsolve.org) how can I draw the normal stress with the correct direction?
Function Norm does not reflect the direction of the normal stress (compressive or tensile). For example, when I use force = CF( (-1e-3,0,0) ) instead of force = CF( (1e-3,0,0) ), I still get the same stress but it should be negative.

I think I have to somehow implement n * stress * n, but don’t know how.

Kind regards
Hamed

2

Hi Hamed,

you can use the specialcf.normal(3) to get the 3D normal vector and compute n * stress * n via

n = specialcf.normal(3)
stress_nn = stress*n*n

Then the stress changes its sign.

Best,
Michael

3

Hi Michael

Thanks for your reply. I did what you mentioned above, but it returns 0 when I print instead of -5.
Here is my code:

from netgen.occ import *
import ngsolve
from ngsolve import *
from ngsolve.webgui import Draw
from ngsolve.krylovspace import CGSolver

box = Box((0,0,0), (3,0.6,1))
geo = box
geo.faces.Min(X).name = “fix”
geo.faces.Max(X).name = “force”
mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=0.1)).Curve(3)

E, nu = 210, 0.2
mu = E / 2 / (1+nu)
lam = E * nu / ((1+nu)(1-2nu))

def Stress(strain):
return 2mustrain + lam*Trace(strain)*Id(3)
fes = VectorH1(mesh, order=3, dirichlet=“fix”)
u,v = fes.TnT()
gfu = GridFunction(fes)

with TaskManager():
a = BilinearForm(InnerProduct(Stress(Sym(Grad(u))), Sym(Grad(v))).Compile()dx)
pre = Preconditioner(a, “bddc”)
a.Assemble()
force = CF((-5,0,0))
f = LinearForm(forcev*ds(“force”)).Assemble()
inv = CGSolver(a.mat, pre, printrates=True, tol=1e-6)
gfu.vec.data = inv * f.vec

with TaskManager():
fesstress = MatrixValued(H1(mesh,order=3), symmetric=True)
gfstress = GridFunction(fesstress)
gfstress.Interpolate (Stress(Sym(Grad(gfu))))

signed_fesstress = H1(mesh,order=3)
gf_signed_stress = GridFunction(signed_fesstress)
n = specialcf.normal(3)
gf_signed_stress.Set(Stress(Sym(Grad(gfu)))*n *n, BND)

gf_none_signed_stress = Norm(gfstress)

print(gf_signed_stress(mesh(3,0.3,0.5)))
print(gf_none_signed_stress(mesh(3,0.3,0.5)))

Can you please tell me what I am doing wrong?

Regards
Hamed

4

replace the line

gf_signed_stress.Set(Stress(Sym(Grad(gfu)))*n *n, BND)

by

gf_signed_stress.Set(gfstress*n *n, BND)

On the boundary you can only access the boundary values of gfu, and you don’t get the stress from that.

5

I replaced

gf_signed_stress.Set(Stress(Sym(Grad(gfu)))*n *n, BND)

with

gf_signed_stress.Set(gfstress*n *n, BND)

but still returns 0. So I decided to check if I could set the normal vector to GridFunction and used the below line:

signed_fesstress = VectorH1(mesh,order=3)
gf_signed_stress = GridFunction(signed_fesstress)
n = specialcf.normal(3)
gf_signed_stress.Set(n,BND)

But returned 0 vector. Now, my question is how I can use specialcf.normal in a GridFunction?

Regards
Hamed

6

the gf.Set(func,BND) sets values only at Dirichlet boundaries. Either define dirichlet everywhere (VectorH1(...,dirichlet=".*"), or use

gf_signed_stress.Set(n, mesh.Boundaries(".*"))
7

It did the trick. Much appreciated!