This page was generated from wta/elasticity3D.ipynb.
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:
Linear strain:
Stress by Hooke’s law:
Equilibrium of forces:
Displacement boundary conditions:
Traction boundary conditions:
Variational formulation:¶
Find:
holds for all
[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)
with TaskManager():
a = BilinearForm(InnerProduct(Stress(Sym(Grad(u))), Sym(Grad(v))).Compile()*dx)
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=True, tol=1e-8)
gfu.vec.data = inv * f.vec
CG iteration 1, residual = 0.00017982731508198196
CG iteration 2, residual = 7.527505823604999e-05
CG iteration 3, residual = 8.667360647586854e-05
CG iteration 4, residual = 6.611763953432786e-05
CG iteration 5, residual = 6.009530495083377e-05
CG iteration 6, residual = 4.240151181641339e-05
CG iteration 7, residual = 3.174659664109285e-05
CG iteration 8, residual = 2.7007176859764265e-05
CG iteration 9, residual = 2.0093499070754638e-05
CG iteration 10, residual = 1.620278465256741e-05
CG iteration 11, residual = 1.1164241035834938e-05
CG iteration 12, residual = 8.707361209404863e-06
CG iteration 13, residual = 6.263300218968363e-06
CG iteration 14, residual = 4.827662227893153e-06
CG iteration 15, residual = 3.5048301823912414e-06
CG iteration 16, residual = 2.6744872507307778e-06
CG iteration 17, residual = 1.9959626253864358e-06
CG iteration 18, residual = 1.4173939615000084e-06
CG iteration 19, residual = 1.0735929420645687e-06
CG iteration 20, residual = 7.342807879185329e-07
CG iteration 21, residual = 6.364131834892706e-07
CG iteration 22, residual = 4.1303078290730905e-07
CG iteration 23, residual = 3.667583721055049e-07
CG iteration 24, residual = 2.5942494636496794e-07
CG iteration 25, residual = 1.805629294458435e-07
CG iteration 26, residual = 1.3642541682195863e-07
CG iteration 27, residual = 9.965383340157361e-08
CG iteration 28, residual = 7.638450216780138e-08
CG iteration 29, residual = 5.282418030039626e-08
CG iteration 30, residual = 4.016696599502615e-08
CG iteration 31, residual = 2.853050965153251e-08
CG iteration 32, residual = 2.171386098777419e-08
CG iteration 33, residual = 1.6955931220410742e-08
CG iteration 34, residual = 1.1309617391619725e-08
CG iteration 35, residual = 8.114589692920307e-09
CG iteration 36, residual = 5.8518862397539e-09
CG iteration 37, residual = 4.02666935143064e-09
CG iteration 38, residual = 3.061722085798082e-09
CG iteration 39, residual = 2.524397820390125e-09
CG iteration 40, residual = 1.6419925053863179e-09
CG iteration 41, residual = 1.158470124444963e-09
CG iteration 42, residual = 9.165155677560578e-10
CG iteration 43, residual = 6.359342431270262e-10
CG iteration 44, residual = 4.120788555112914e-10
CG iteration 45, residual = 3.475360816124547e-10
CG iteration 46, residual = 2.913862570760715e-10
CG iteration 47, residual = 1.92728323632225e-10
CG iteration 48, residual = 1.4408886051090374e-10
CG iteration 49, residual = 1.122426507311188e-10
CG iteration 50, residual = 9.517641882334069e-11
CG iteration 51, residual = 6.691289015672579e-11
CG iteration 52, residual = 4.132646365973799e-11
CG iteration 53, residual = 3.0963083090874015e-11
CG iteration 54, residual = 2.025614648097733e-11
CG iteration 55, residual = 1.5070686043116477e-11
CG iteration 56, residual = 1.1829538190137948e-11
CG iteration 57, residual = 7.406493281312491e-12
CG iteration 58, residual = 5.294067613357582e-12
CG iteration 59, residual = 3.453776593803619e-12
CG iteration 60, residual = 2.390872304120967e-12
CG iteration 61, residual = 1.7659815586180077e-12
[10]:
with TaskManager():
fesstress = MatrixValued(H1(mesh,order=3), symmetric=True)
gfstress = GridFunction(fesstress)
gfstress.Interpolate (Stress(Sym(Grad(gfu))))
[11]:
Draw (5e4*gfu, mesh);
[12]:
Draw (Norm(gfstress), mesh, deformation=1e4*gfu, draw_vol=False, order=3);
[ ]: