This page was generated from wta/elasticity.ipynb.
Solving nonlinear Elasticity¶
[1]:
from ngsolve import *
from ngsolve.webgui import Draw
a rectangle with refinement at corners:
[2]:
from netgen.occ import *
shape = Rectangle(1,0.1).Face()
shape.edges.Max(X).name="right"
shape.edges.Min(X).name="left"
shape.edges.Max(Y).name="top"
shape.edges.Min(Y).name="bot"
shape.vertices.Min(X+Y).maxh=0.01
shape.vertices.Min(X-Y).maxh=0.01
mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.05))
Cauchy-Green tensor
\[C = F^T F \qquad \text{with} \qquad F = I + \nabla u\]
and hyperelastic energy density
\[W : {\mathbb R}^{d \times d, sym} \rightarrow {\mathbb R}\]
[3]:
E, nu = 210, 0.2
mu = E / 2 / (1+nu)
lam = E * nu / ((1+nu)*(1-2*nu))
def C(u):
F = Id(2) + Grad(u)
return F.trans * F
def NeoHooke (C):
return 0.5*mu*(Trace(C-Id(2)) + 2*mu/lam*Det(C)**(-lam/2/mu)-1)
stationary point of total energy:
\[\delta \int W(C(u)) - f u = 0\]
[4]:
factor = Parameter(0)
force = CoefficientFunction( (0,factor) )
fes = H1(mesh, order=4, dirichlet="left", dim=mesh.dim)
u = fes.TrialFunction()
a = BilinearForm(fes, symmetric=True)
a += Variation(NeoHooke(C(u)).Compile()*dx)
a += Variation((-InnerProduct(force,u)).Compile()*dx)
gfu = GridFunction(fes)
gfu.vec[:] = 0
The Variation function declares that the non-linear form is the derivative of the energy.
a simple Newton solver, using automatic differentiation for residual and tangential stiffness:
[5]:
def SolveNewton(printrates=False):
for it in range(10):
if (printrates):
print ("it", it, "energy = ", a.Energy(gfu.vec))
res = a.Apply(gfu.vec)
a.AssembleLinearization(gfu.vec)
inv = a.mat.Inverse(fes.FreeDofs() )
gfu.vec.data -= inv*res
[6]:
factor.Set(0.4)
SolveNewton(printrates=True)
scene = Draw (C(gfu)[0,0]-1, mesh, deformation=gfu, min=-0.1, max=0.1)
it 0 energy = 8.749999999999972
it 1 energy = 8.811175559205765
it 2 energy = 8.748116767882797
it 3 energy = 8.747829234525454
it 4 energy = 8.747829153671926
it 5 energy = 8.747829153657275
it 6 energy = 8.747829153657271
it 7 energy = 8.747829153657271
it 8 energy = 8.747829153657271
it 9 energy = 8.747829153657271
Often, we don’t have a good starting value for Newton’s method. This can be overcome by increasing the load step by step (assuming the solution depends continuously on the loading). The solution of the previous load-step is the initial guess for the next step.
[7]:
numsteps = 5
maxload = 2
for ls in range (numsteps):
factor.Set(maxload*(ls+1)/numsteps)
SolveNewton()
Draw (C(gfu)[0,0]-1, mesh, deformation=gfu, min=-0.2, max=0.2)
Stress tensor¶
Compute \(2^{nd}\) Piola Kirchhoff stress tensor by symbolic differentiation:
\[\Sigma_{i,j} = \frac{\partial W}{\partial C_{i,j}} (C)\]
[8]:
C_=C(gfu).MakeVariable()
sigma = NeoHooke(C_).Diff(C_)
Draw (sigma[0,0], mesh, "Sxx", deformation=gfu, min=-10.001, max=10.001);
The energy functional is represented as an expression tree:
[9]:
u = fes.TrialFunction()
print (NeoHooke(C(u)))
43.75*(coef binary operation '-', real
coef binary operation '+', real
coef trace, real
coef binary operation '-', real, dims = 2 x 2
coef matrix-matrix multiply, real, dims = 2 x 2
coef Matrix transpose, real, dims = 2 x 2
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef trial-function diffop = grad, real, dims = 2 x 2
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef trial-function diffop = grad, real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef scale 3, real
coef binary operation 'pow', real
coef binary operation '*', real
coef Determinant, real
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef trial-function diffop = grad, real, dims = 2 x 2
coef Determinant, real
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef trial-function diffop = grad, real, dims = 2 x 2
coef -0.333333, real
coef 1, real
)
With the Compile method, the tree is linearized, and common sub-expressions are merged:
[10]:
print (NeoHooke(C(u)).Compile())
Compiled CF:
Step 0: Identity matrix, dims = 2 x 2
Step 1: trial-function diffop = grad, dims = 2 x 2
Step 2: binary operation '+', dims = 2 x 2
input: 0 1
Step 3: Matrix transpose, dims = 2 x 2
input: 2
Step 4: matrix-matrix multiply, dims = 2 x 2
input: 3 2
Step 5: Identity matrix, dims = 2 x 2
Step 6: binary operation '-', dims = 2 x 2
input: 4 5
Step 7: trace
input: 6
Step 8: Determinant
input: 2
Step 9: binary operation '*'
input: 8 8
Step 10: -0.333333
Step 11: binary operation 'pow'
input: 9 10
Step 12: scale 3
input: 11
Step 13: binary operation '+'
input: 7 12
Step 14: 1
Step 15: binary operation '-'
input: 13 14
Step 16: scale 43.75
input: 15
[11]:
print (sigma)
43.75*(coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef scale 3, real, dims = 2 x 2
coef reshape, real, dims = 2 x 2
coef Matrix transpose, real, dims = 1 x 4
coef reshape, real, dims = 4 x 1
coef VectorialCoefficientFunction, real, dim=4
coef binary operation '*', real
coef binary operation '*', real
coef -0.333333, real
coef binary operation 'pow', real
coef Determinant, real
coef matrix-matrix multiply, real, dims = 2 x 2
coef Matrix transpose, real, dims = 2 x 2
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef N6ngcomp31GridFunctionCoefficientFunctionE, real, dims = 2 x 2
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef N6ngcomp31GridFunctionCoefficientFunctionE, real, dims = 2 x 2
coef binary operation '-', real
coef -0.333333, real
coef 1, real
coef innerproduct, fix size = 4, real
coef cofactor, real, dims = 2 x 2
coef matrix-matrix multiply, real, dims = 2 x 2
coef Matrix transpose, real, dims = 2 x 2
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef N6ngcomp31GridFunctionCoefficientFunctionE, real, dims = 2 x 2
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef N6ngcomp31GridFunctionCoefficientFunctionE, real, dims = 2 x 2
coef reshape, real, dims = 2 x 2
coef UnitVectorCF 0, real, dim=4
coef binary operation '*', real
coef binary operation '*', real
coef -0.333333, real
coef binary operation 'pow', real
coef Determinant, real
coef matrix-matrix multiply, real, dims = 2 x 2
coef Matrix transpose, real, dims = 2 x 2
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef N6ngcomp31GridFunctionCoefficientFunctionE, real, dims = 2 x 2
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef N6ngcomp31GridFunctionCoefficientFunctionE, real, dims = 2 x 2
coef binary operation '-', real
coef -0.333333, real
coef 1, real
coef innerproduct, fix size = 4, real
coef cofactor, real, dims = 2 x 2
coef matrix-matrix multiply, real, dims = 2 x 2
coef Matrix transpose, real, dims = 2 x 2
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef N6ngcomp31GridFunctionCoefficientFunctionE, real, dims = 2 x 2
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef N6ngcomp31GridFunctionCoefficientFunctionE, real, dims = 2 x 2
coef reshape, real, dims = 2 x 2
coef UnitVectorCF 1, real, dim=4
coef binary operation '*', real
coef binary operation '*', real
coef -0.333333, real
coef binary operation 'pow', real
coef Determinant, real
coef matrix-matrix multiply, real, dims = 2 x 2
coef Matrix transpose, real, dims = 2 x 2
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef N6ngcomp31GridFunctionCoefficientFunctionE, real, dims = 2 x 2
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef N6ngcomp31GridFunctionCoefficientFunctionE, real, dims = 2 x 2
coef binary operation '-', real
coef -0.333333, real
coef 1, real
coef innerproduct, fix size = 4, real
coef cofactor, real, dims = 2 x 2
coef matrix-matrix multiply, real, dims = 2 x 2
coef Matrix transpose, real, dims = 2 x 2
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef N6ngcomp31GridFunctionCoefficientFunctionE, real, dims = 2 x 2
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef N6ngcomp31GridFunctionCoefficientFunctionE, real, dims = 2 x 2
coef reshape, real, dims = 2 x 2
coef UnitVectorCF 2, real, dim=4
coef binary operation '*', real
coef binary operation '*', real
coef -0.333333, real
coef binary operation 'pow', real
coef Determinant, real
coef matrix-matrix multiply, real, dims = 2 x 2
coef Matrix transpose, real, dims = 2 x 2
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef N6ngcomp31GridFunctionCoefficientFunctionE, real, dims = 2 x 2
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef N6ngcomp31GridFunctionCoefficientFunctionE, real, dims = 2 x 2
coef binary operation '-', real
coef -0.333333, real
coef 1, real
coef innerproduct, fix size = 4, real
coef cofactor, real, dims = 2 x 2
coef matrix-matrix multiply, real, dims = 2 x 2
coef Matrix transpose, real, dims = 2 x 2
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef N6ngcomp31GridFunctionCoefficientFunctionE, real, dims = 2 x 2
coef binary operation '+', real, dims = 2 x 2
coef Identity matrix, real, dims = 2 x 2
coef N6ngcomp31GridFunctionCoefficientFunctionE, real, dims = 2 x 2
coef reshape, real, dims = 2 x 2
coef UnitVectorCF 3, real, dim=4
)
[12]:
print (sigma.Compile())
Compiled CF:
Step 0: Identity matrix, dims = 2 x 2
Step 1: -0.333333
Step 2: Identity matrix, dims = 2 x 2
Step 3: N6ngcomp31GridFunctionCoefficientFunctionE, dims = 2 x 2
Step 4: binary operation '+', dims = 2 x 2
input: 2 3
Step 5: Matrix transpose, dims = 2 x 2
input: 4
Step 6: matrix-matrix multiply, dims = 2 x 2
input: 5 4
Step 7: Determinant
input: 6
Step 8: 1
Step 9: binary operation '-'
input: 1 8
Step 10: binary operation 'pow'
input: 7 9
Step 11: binary operation '*'
input: 1 10
Step 12: cofactor, dims = 2 x 2
input: 6
Step 13: UnitVectorCF 0, dim=4
Step 14: reshape, dims = 2 x 2
input: 13
Step 15: innerproduct, fix size = 4
input: 12 14
Step 16: binary operation '*'
input: 11 15
Step 17: 1
Step 18: binary operation '-'
input: 1 17
Step 19: binary operation 'pow'
input: 7 18
Step 20: binary operation '*'
input: 1 19
Step 21: cofactor, dims = 2 x 2
input: 6
Step 22: UnitVectorCF 1, dim=4
Step 23: reshape, dims = 2 x 2
input: 22
Step 24: innerproduct, fix size = 4
input: 21 23
Step 25: binary operation '*'
input: 20 24
Step 26: 1
Step 27: binary operation '-'
input: 1 26
Step 28: binary operation 'pow'
input: 7 27
Step 29: binary operation '*'
input: 1 28
Step 30: cofactor, dims = 2 x 2
input: 6
Step 31: UnitVectorCF 2, dim=4
Step 32: reshape, dims = 2 x 2
input: 31
Step 33: innerproduct, fix size = 4
input: 30 32
Step 34: binary operation '*'
input: 29 33
Step 35: 1
Step 36: binary operation '-'
input: 1 35
Step 37: binary operation 'pow'
input: 7 36
Step 38: binary operation '*'
input: 1 37
Step 39: cofactor, dims = 2 x 2
input: 6
Step 40: UnitVectorCF 3, dim=4
Step 41: reshape, dims = 2 x 2
input: 40
Step 42: innerproduct, fix size = 4
input: 39 41
Step 43: binary operation '*'
input: 38 42
Step 44: VectorialCoefficientFunction, dim=4
input: 16 25 34 43
Step 45: reshape, dims = 4 x 1
input: 44
Step 46: Matrix transpose, dims = 1 x 4
input: 45
Step 47: reshape, dims = 2 x 2
input: 46
Step 48: scale 3, dims = 2 x 2
input: 47
Step 49: binary operation '+', dims = 2 x 2
input: 0 48
Step 50: scale 43.75, dims = 2 x 2
input: 49
[ ]:
[ ]:
[ ]: