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.81117555920626
it 2 energy = 8.7481167678828
it 3 energy = 8.747829234525453
it 4 energy = 8.747829153671926
it 5 energy = 8.747829153657277
it 6 energy = 8.747829153657273
it 7 energy = 8.747829153657273
it 8 energy = 8.747829153657271
it 9 energy = 8.747829153657273
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 unary operation ' ', real
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: unary operation ' '
input: 10
Step 12: binary operation 'pow'
input: 9 11
Step 13: scale 3
input: 12
Step 14: binary operation '+'
input: 7 13
Step 15: 1
Step 16: binary operation '-'
input: 14 15
Step 17: scale 43.75
input: 16
[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 unary 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 unary 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 unary 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 unary 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 unary 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 unary 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 unary 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 unary 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: unary operation ' '
input: 1
Step 3: Identity matrix, dims = 2 x 2
Step 4: N6ngcomp31GridFunctionCoefficientFunctionE, dims = 2 x 2
Step 5: binary operation '+', dims = 2 x 2
input: 3 4
Step 6: Matrix transpose, dims = 2 x 2
input: 5
Step 7: matrix-matrix multiply, dims = 2 x 2
input: 6 5
Step 8: Determinant
input: 7
Step 9: 1
Step 10: binary operation '-'
input: 2 9
Step 11: binary operation 'pow'
input: 8 10
Step 12: binary operation '*'
input: 2 11
Step 13: cofactor, dims = 2 x 2
input: 7
Step 14: UnitVectorCF 0, dim=4
Step 15: reshape, dims = 2 x 2
input: 14
Step 16: innerproduct, fix size = 4
input: 13 15
Step 17: binary operation '*'
input: 12 16
Step 18: 1
Step 19: binary operation '-'
input: 2 18
Step 20: binary operation 'pow'
input: 8 19
Step 21: binary operation '*'
input: 2 20
Step 22: cofactor, dims = 2 x 2
input: 7
Step 23: UnitVectorCF 1, dim=4
Step 24: reshape, dims = 2 x 2
input: 23
Step 25: innerproduct, fix size = 4
input: 22 24
Step 26: binary operation '*'
input: 21 25
Step 27: 1
Step 28: binary operation '-'
input: 2 27
Step 29: binary operation 'pow'
input: 8 28
Step 30: binary operation '*'
input: 2 29
Step 31: cofactor, dims = 2 x 2
input: 7
Step 32: UnitVectorCF 2, dim=4
Step 33: reshape, dims = 2 x 2
input: 32
Step 34: innerproduct, fix size = 4
input: 31 33
Step 35: binary operation '*'
input: 30 34
Step 36: 1
Step 37: binary operation '-'
input: 2 36
Step 38: binary operation 'pow'
input: 8 37
Step 39: binary operation '*'
input: 2 38
Step 40: cofactor, dims = 2 x 2
input: 7
Step 41: UnitVectorCF 3, dim=4
Step 42: reshape, dims = 2 x 2
input: 41
Step 43: innerproduct, fix size = 4
input: 40 42
Step 44: binary operation '*'
input: 39 43
Step 45: VectorialCoefficientFunction, dim=4
input: 17 26 35 44
Step 46: reshape, dims = 4 x 1
input: 45
Step 47: Matrix transpose, dims = 1 x 4
input: 46
Step 48: reshape, dims = 2 x 2
input: 47
Step 49: scale 3, dims = 2 x 2
input: 48
Step 50: binary operation '+', dims = 2 x 2
input: 0 49
Step 51: scale 43.75, dims = 2 x 2
input: 50
[ ]:
[ ]:
[ ]: