3D Solid Mechanics¶

[1]:

from netgen.occ import *
from netgen.webgui import Draw as DrawGeo
import ngsolve

[2]:

box = Box((0,0,0), (3,0.6,1))
box.faces.name="outer"
cyl = sum( [Cylinder((0.5+i,0,0.5), Y, 0.25,0.8) for i in range(3)] )
cyl.faces.name="cyl"
geo = box-cyl

DrawGeo(geo);


find edges between box and cylinder, and build chamfers (requires OCC 7.4 or newer):

[3]:

cylboxedges = geo.faces["outer"].edges * geo.faces["cyl"].edges
cylboxedges.name = "cylbox"
geo = geo.MakeChamfer(cylboxedges, 0.03)


name faces for boundary conditions:

[4]:

geo.faces.Min(X).name = "fix"
geo.faces.Max(X).name = "force"

DrawGeo(geo);

[5]:

from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=0.1)).Curve(3)
Draw (mesh);


Linear elasticity¶

Displacement: $$u : \Omega \rightarrow R^3$$

Linear strain:

$\varepsilon(u) := \tfrac{1}{2} ( \nabla u + (\nabla u)^T )$

Stress by Hooke’s law:

$\sigma = 2 \mu \varepsilon + \lambda \operatorname{tr} \varepsilon I$

Equilibrium of forces:

$\operatorname{div} \sigma = f$

Displacement boundary conditions:

$u = u_D \qquad \text{on} \, \Gamma_D$

Traction boundary conditions:

$\sigma n = g \qquad \text{on} \, \Gamma_N$

Variational formulation:¶

Find: $$u \in H^1(\Omega)^3$$ such that $$u = u_D$$ on $$\Gamma_D$$

$\int_\Omega \sigma(\varepsilon(u)) : \varepsilon(v) \, dx = \int_\Omega f v dx + \int_{\Gamma_N} g v ds$

holds for all $$v = 0$$ on $$\Gamma_D$$.

[6]:

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

def Stress(strain):
return 2*mu*strain + lam*Trace(strain)*Id(3)

[7]:

fes = VectorH1(mesh, order=3, dirichlet="fix")
u,v = fes.TnT()
gfu = GridFunction(fes)

pre = Preconditioner(a, "bddc")
a.Assemble()

[8]:

force = CF( (1e-3,0,0) )
f = LinearForm(force*v*ds("force")).Assemble()

[9]:

from ngsolve.krylovspace import CGSolver
inv = CGSolver(a.mat, pre, printrates='\r', tol=1e-8)
gfu.vec.data = inv * f.vec

CG converged in 59 iterations to residual 1.2332950480961994e-12

[10]:

with TaskManager():
fesstress = MatrixValued(H1(mesh,order=3), symmetric=True)
gfstress = GridFunction(fesstress)

[11]:

Draw (5e4*gfu, mesh);

[12]:

Draw (Norm(gfstress), mesh, deformation=1e4*gfu, draw_vol=False, order=3);

[ ]: