This page was generated from wta/elasticity3D.ipynb.

3D Solid Mechanics

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

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):

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

name faces for boundary conditions:

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

DrawGeo(geo);

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\).

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)

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( (1e-3,0,0) )
f = LinearForm(force*v*ds("force")).Assemble()

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

CG iteration 1, residual = 0.0001798904505629877
CG iteration 2, residual = 7.593089644801455e-05
CG iteration 3, residual = 9.553190911073941e-05
CG iteration 4, residual = 6.915155168636508e-05
CG iteration 5, residual = 6.408034527301741e-05
CG iteration 6, residual = 4.6957231935788e-05
CG iteration 7, residual = 4.0742471318053335e-05
CG iteration 8, residual = 2.96435104549506e-05
CG iteration 9, residual = 2.3189596705094882e-05
CG iteration 10, residual = 1.7242540984822337e-05
CG iteration 11, residual = 1.3840252847566043e-05
CG iteration 12, residual = 1.1371290859109155e-05
CG iteration 13, residual = 7.557236098508591e-06
CG iteration 14, residual = 5.772464004177619e-06
CG iteration 15, residual = 4.233341381713284e-06
CG iteration 16, residual = 3.4657844607239695e-06
CG iteration 17, residual = 2.362507411751025e-06
CG iteration 18, residual = 1.960334962442719e-06
CG iteration 19, residual = 1.468798970318026e-06
CG iteration 20, residual = 1.1178053345310865e-06
CG iteration 21, residual = 7.786087109406145e-07
CG iteration 22, residual = 5.912354552243875e-07
CG iteration 23, residual = 4.871534255435866e-07
CG iteration 24, residual = 3.1939698996758574e-07
CG iteration 25, residual = 2.6448994915477027e-07
CG iteration 26, residual = 1.7131785458603375e-07
CG iteration 27, residual = 1.5039339583039745e-07
CG iteration 28, residual = 1.0180139619024763e-07
CG iteration 29, residual = 7.508935450463193e-08
CG iteration 30, residual = 5.4855918275150506e-08
CG iteration 31, residual = 4.239828755821247e-08
CG iteration 32, residual = 2.7875622339812555e-08
CG iteration 33, residual = 2.1924668675301854e-08
CG iteration 34, residual = 1.6129431810821462e-08
CG iteration 35, residual = 1.11345061976361e-08
CG iteration 36, residual = 8.848793463601322e-09
CG iteration 37, residual = 5.85370483196082e-09
CG iteration 38, residual = 4.719941827518409e-09
CG iteration 39, residual = 3.069282116372148e-09
CG iteration 40, residual = 2.2827977160230614e-09
CG iteration 41, residual = 1.5635755233329909e-09
CG iteration 42, residual = 1.5013392923226421e-09
CG iteration 43, residual = 1.0839989468650738e-09
CG iteration 44, residual = 8.32444435970913e-10
CG iteration 45, residual = 5.034692397814603e-10
CG iteration 46, residual = 3.9895789292461415e-10
CG iteration 47, residual = 3.069599794676291e-10
CG iteration 48, residual = 2.3135237478378914e-10
CG iteration 49, residual = 2.135201982960793e-10
CG iteration 50, residual = 1.3488363337094907e-10
CG iteration 51, residual = 9.302020111971208e-11
CG iteration 52, residual = 6.915832693342542e-11
CG iteration 53, residual = 4.868535127548191e-11
CG iteration 54, residual = 5.104112901521503e-11
CG iteration 55, residual = 2.8151726040131484e-11
CG iteration 56, residual = 1.987090363402689e-11
CG iteration 57, residual = 1.4158216831531278e-11
CG iteration 58, residual = 1.0499583711922871e-11
CG iteration 59, residual = 7.383215719277148e-12
CG iteration 60, residual = 4.9350302332638065e-12
CG iteration 61, residual = 4.1698428291659885e-12
CG iteration 62, residual = 2.756398229507717e-12
CG iteration 63, residual = 1.950453609466774e-12
CG iteration 64, residual = 2.0065685812963534e-12
CG iteration 65, residual = 1.0498860864697147e-12

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

Draw (5e4*gfu, mesh);

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