Non-linear simulation of a hyperelastic beam

We consider a simple beam with a hyperelastic material model. In particular, we will assume the beam has energy:

\[E(u) := \int_{\Omega} \frac{\mu}{2} tr(\mathbb{C}-I)+ \frac{2\mu}{\lambda} \det(\mathbb{C})^{-\frac{\lambda}{2\mu}-1}\, dx + \int_{\partial \Omega} \vec{f} \cdot \vec{u} \, ds\]

where \(\mathbb{C} = F^T F\) is the right Cauchy-Green tensor, \(F\) is the deformation gradient, \(\mu\) and \(\lambda\) are the Lamé parameters, and \(\vec{f}\) is the force applied to the top face of the beam. A simple discretisation of this energy leads to a non-linear problem that we solve using PETSc SNES.

from ngsolve import *
import netgen.gui
from netgen.occ import *
import netgen.meshing as ngm
from mpi4py.MPI import COMM_WORLD

if COMM_WORLD.rank == 0:
    d = 0.01
    box = Box ( (-d/2,-d/2,0), (d/2,d/2,0.1) ) + Box( (-d/2, -3*d/2,0.1), (d/2, 3*d/2, 0.1+d) )
    box.faces.Min(Z).name = "bottom"
    box.faces.Max(Z).name = "top"
    mesh = Mesh(OCCGeometry(box).GenerateMesh(maxh=0.05).Distribute(COMM_WORLD))
else:
    mesh = Mesh(ngm.Mesh.Receive(COMM_WORLD))

E, nu = 210, 0.2
mu  = E / 2 / (1+nu)
lam = E * nu / ((1+nu)*(1-2*nu))

def C(u):
    F = Id(u.dim) + Grad(u)
    return F.trans * F

def NeoHooke (C):
    # return 0.5*mu*InnerProduct(C-Id(3), C-Id(3))
    return 0.5*mu*(Trace(C-Id(3)) + 2*mu/lam*Det(C)**(-lam/2/mu)-1)

loadfactor = Parameter(1)
force = loadfactor * CF ( (-y, x, 0) )

fes = H1(mesh, order=3, dirichlet="bottom", dim=mesh.dim)
u,v = fes.TnT()

a = BilinearForm(fes)
a += Variation(NeoHooke(C(u)).Compile()*dx)
a += ((Id(3)+Grad(u.Trace()))*force)*v*ds("top")

Once we have defined the energy and the weak form, we can solve the non-linear problem using PETSc SNES. In particular, we will use a Newton method with line search, solving the arising linear systems with a direct solver.

from ngsPETSc import NonLinearSolver
gfu_petsc = GridFunction(fes)
gfu_ngs = GridFunction(fes)
gfu_ngs.vec[:] = 0; gfu_petsc.vec[:] = 0
gfu_history_ngs = GridFunction(fes, multidim=0)
gfu_history_petsc = GridFunction(fes, multidim=0)
numsteps = 30
for step in range(numsteps):
    print("step", step)
    loadfactor.Set(720*step/numsteps)
    solver = NonLinearSolver(fes, a=a, objective=False,
                             solverParameters={"snes_type": "newtonls",
                                               "snes_max_it": 200,
                                               "snes_monitor": "",
                                               "ksp_type": "preonly",
                                               "pc_type": "lu",
                                               "snes_linesearch_type": "basic",
                                               "snes_linesearch_damping": 1.0,
                                               "snes_linesearch_max_it": 100})
    gfu_petsc = solver.solve(gfu_petsc)
    solvers.Newton (a, gfu_ngs, printing=True, dampfactor=0.5)
    gfu_history_ngs.AddMultiDimComponent(gfu_ngs.vec)
    gfu_history_petsc.AddMultiDimComponent(gfu_petsc.vec)

We compare the performance of the PETSc SNES solvers and of NGSolve’s own implementation of Newton’s method:

Performance of different non-linear solvers

Solver

non-linear iteration

NGS Newton (damping=1.0)

Diverged

NGS Newton (damping=0.5)

Diverged

NGS Newton (damping=0.3)

12

PETSc SNES

10

This suggests that while NGSolve’s non-linear solver when finely tuned performs as well as PETSc SNES, it is more sensitive to the choice of the damping factor. In this case, a damping factor of 0.3 was found to be the best choice.