Solving nonlinear Elasticity

54. Solving nonlinear Elasticity#

from ngsolve import *
from ngsolve.webgui import Draw

a rectangle with refinement at corners:

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} \]

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 \]

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:

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

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.811175555530797
it 2 energy =  8.748116767901415
it 3 energy =  8.747829234576708
it 4 energy =  8.74782915372321
it 5 energy =  8.747829153708558
it 6 energy =  8.747829153708558
it 7 energy =  8.747829153708558
it 8 energy =  8.747829153708558
it 9 energy =  8.747829153708558

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.

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)

54.1. Stress tensor#

Compute \(2^{nd}\) Piola Kirchhoff stress tensor by symbolic differentiation:

\[ \Sigma_{i,j} = \frac{\partial W}{\partial C_{i,j}} (C) \]

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:

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:

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 

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 scalar-matrix multiply, real, dims = 2 x 2
      coef unary operation 'exp', real
        coef binary operation '*', real
          coef -0.333333, real
          coef unary operation 'log', 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 scalar-matrix multiply, real, dims = 2 x 2
        coef -0.333333, real
        coef scalar-matrix multiply, real, dims = 2 x 2
          coef binary operation '/', real
            coef 1, 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 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
)

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: unary operation 'log'
     input: 7 
Step 9: binary operation '*'
     input: 1 8 
Step 10: unary operation 'exp'
     input: 9 
Step 11: 1
Step 12: binary operation '/'
     input: 11 7 
Step 13: cofactor, dims = 2 x 2
     input: 6 
Step 14: scalar-matrix multiply, dims = 2 x 2
     input: 12 13 
Step 15: scalar-matrix multiply, dims = 2 x 2
     input: 1 14 
Step 16: scalar-matrix multiply, dims = 2 x 2
     input: 10 15 
Step 17: scale 3, dims = 2 x 2
     input: 16 
Step 18: binary operation '+', dims = 2 x 2
     input: 0 17 
Step 19: scale 43.75, dims = 2 x 2
     input: 18 

Solving nonlinear Elasticity

Contents

54. Solving nonlinear Elasticity#

54.1. Stress tensor#