Hyperelasticity + Follower Loads

NGSolve
1

Hi all,

I started to use NGSolve for solving some problems in modeling hyperelastic materials. I followed the examples here closely.

The example I want to calculate is the inflation of an ellipsoidal structure.

I generated my mesh as follows (units are in mm)

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw
from ngsolve.solvers import NewtonMinimization
ngsglobals.msg_level = 1

import matplotlib.pyplot as plt
from scipy.stats import qmc
import scipy as sp
import numpy as np
np.random.seed(137)


rs_endo = 2.5e-2 * 10 ** 3
rl_endo = 9.0e-2 * 10 ** 3
rs_epi  = 3.5e-2 * 10 ** 3
rl_epi  = 9.7e-2 * 10 ** 3

ax = Axes((1,0,0), n=X, h=Y)
ax2 = Axes((1,0,0), n=X, h=Y)
maxh=3.0

ell_endo = Ellipsoid (Axis ((0,0,0), X), rl_endo, rs_endo)
ell_endo.faces.name = "endocardium"
ell_endo.faces.col = (0.5,0.5,0)

ell_epi = Ellipsoid (Axis ((0,0,0), X), rl_epi, rs_epi)
ell_epi.faces.name = "epicardium"
ell_epi.faces.col = (1,0,0)

cutplane=HalfSpace((-25,0,0), X)
cutplane.faces.name="top"
cutplane.faces.col=(0,0,1)
lv = (ell_epi-ell_endo)-cutplane

geo = OCCGeometry(lv)
mesh=Mesh(geo.GenerateMesh(maxh=maxh))

This works fine. Now I can replicate some of the features laid out in the elasticity examples of the NGS24 workshop. I adapted things for almost incompressible Elasticity

fes = VectorH1(mesh, order=1, dirichlet="top")
u, deltau = fes.TnT()
gfu = GridFunction(fes)
loadpar = Parameter(0)
normal = specialcf.normal(mesh.dim)

kPa2MPa = 1e-3
mu = 33.3 * kPa2MPa
kappa = 1000 * kPa2MPa
appl_pressure = Parameter(0.0) # in kPa

def psi(F):
    C = F.trans * F
    J = Det(F)
    Cbar = J**(-2/3) * C
    Wiso = mu/2*(Trace(Cbar)-3)
    Wvol = kappa/2 * log(J)**2
    return Wiso + Wvol

a = BilinearForm(fes, symmetric=False)
defgrad = Id(mesh.dim) + Grad(u)

# Energy volumetric
a += Variation(psi(defgrad)*dx)

# dead load (pressure acts against normal direction this swallows the minus sign
a += appl_pressure *kPa2MPa*InnerProduct(normal, deltau)*ds("endocardium")

def SolveNonlinearElast(fes, pressure, numloadsteps=10):
    print("Solve Nonlinear Elast for p=", pressure)
    # Solution
    gfu = ng.GridFunction(fes)
    gfu.vec[:] = 0
    # Set params
    incpr = pressure / numloadsteps
    newval = 0.0
    conv = 0
    converged  = True
    with TaskManager():
        for loadstep in range(numloadsteps):
            print ("loadstep", loadstep+1, end='')
            newval += incpr
            print(" | val", newval)
            appl_pressure .Set(newval)
            i1, _ = NewtonMinimization(a, gfu, printing=True, linesearch=False, inverse="pardiso", maxerr=1e-6, maxit=10)
            conv = conv + i1
    if conv < 0:
        converged = False
    return gfu, converged

However, what I would actually like to model is a followerload which is defined as

P\boldsymbol{N} = -p_\mathrm{appl}JF^{-\top}\boldsymbol{N} \quad\quad\mathrm{on}~\Gamma_\mathrm{follow}

with
F = I + \nabla u,~J=\mathrm{det}(F),~P=\frac{\partial \Psi(F)}{\partial F} and N being the normal vector in the reference configuration (aka specialcf.normal).

The above example actually models
P\boldsymbol{N} = -p_\mathrm{appl}\boldsymbol{N} \quad\quad\mathrm{on}~\Gamma_\mathrm{follow}

Since follower loads don’t allow for a potential formulation I wanted to ask how I could introduce this here in my example?

2

Hi Elias,

as you say, you don’t have a potential, so you cannot use the NewtonMinimization -solver (which does a line-search), but you can use the built-in Newton - solver.

You want the normal on the deformed configuration:

n = J F^{-T} N = Cof(F) N

On the boundary, we don’t have access to the whole \nabla u, only to the surface gradient \nabla_S u which can be computed from u on the boundary. If my quick guess is correct, then

Cof(F) N = Cof(I + \nabla_S u) N

best, Joachim

3

Hi Joachim,
thanks for the quick reply

so if I do something like

fes = ng.VectorH1(mesh, order=1, dirichlet="top")
#fes = ng.H1(mesh, order=1,dim=mesh.dim)
u_elast = fes.TrialFunction()
v_elast = fes.TestFunction()
ubd = u_elast.Trace()
Fbd = F(ubd)
FbdT = Fbd.trans
J = ng.Det(Fbd)
cofactor = J * ng.Inv(FbdT)
normal = ng.specialcf.normal(3)
a_elast = ng.BilinearForm(fes, symmetric=False)
a_elast += ng.Variation( Psi(F(u_elast)).Compile(realcompile=False)*ng.dx
a_elast += inn_val * (ng.InnerProduct(cofactor * normal, v_elast)).Compile(realcompile=False) * ng.ds(definedon="endocardium")

this should work? I’m still a bit confused that in some examples you set realcompile=False and in some you don’t, what’s the difference here?

Cheers
Elias

4

here is a notebook explaining compile and real-compile:

https://jschoeberl.github.io/talk-pdesoft/wta/elasticity.html

Obviously, for real-compile you need a compiler installed. Thus, we avoid it in the tutorials.
Joachim