TDNNS method for nonlinear elasticity

6.5. TDNNS method for nonlinear elasticity#

For nonlinear, invertible material laws the TDNNS method for linear elasticity can not be directly extended to the nonlinear regime. The main problem is that the gradient of the displacement field \(u\in H(\mathrm{curl})\) is a distribution rather than a function. Therefore, multiplication is in general not well defined.

We use the Hu-Washizu principle following [Neunteufel, Pechstein, Schöberl. Three-field mixed finite element methods for nonlinear elasticity. Computer Methods in Applied Mechanics and Engineering, (2021).] introducing a new independent strain field. We “lift” the distribution \(\nabla u +I\) to the deformation gradient \(F\) as a new unknown, which will be a regular function again.

For a hyperelastic potential \(\Psi(\cdot)\) we define the Lagrangian

\[\begin{align*} \mathcal{L}(u, F,P)=\int_{\Omega}W(F)\,dx - \int_\Omega (F-I-\nabla u):P\, dx-\int_{\Omega}f\cdot u\, dx, \end{align*}\]

with the first Piola-Kirchhoff stress tensor \(P=\frac{\partial \Psi}{\partial F}\) as Lagrange multiplier to force \(F=I+\nabla u\). As \(\nabla u\) is not in \(L^2\) we use the TDNNS pairing \(\langle\sigma,\nabla u\rangle=\sum_T\big(\int_T \sigma:\nabla u\,dx-\int_{\partial T} \sigma_{nn}u_n\,ds\big)=-\langle \mathrm{div}(\sigma),u\rangle\).

As \(F\) and \(P\) are in general non-symmetric we use the augmented \(H(\mathrm{div}\mathrm{div})\) space with skew-symmetric \(L^2\)-matrices

\[\begin{align*} \tilde{H}(\mathrm{div}\mathrm{div})=H(\mathrm{div}\mathrm{div})\times [L^2]^{d\times d}_{\mathrm{skew}},\qquad \tilde{\Sigma}_h=\Sigma_h\times [Q_h]^{d\times d}_{\mathrm{skew}}, \end{align*}\]

where \(\Sigma_h\) and \(Q_h\) are the Hellan-Herrmann-Johnson and discontinuous finite elements, respectively.

Therefore, the TDNNS method for nonlinear elasticity reads: Find \((u,F,P)\in H(\mathrm{curl})\times [L^2]^{d\times d}\times \tilde{H}(\mathrm{div}\mathrm{div})\) for the Lagrangian

\[\begin{align*} \mathcal{L}(u, F, P)=\int_{\Omega}W(F)\,dx - \langle F-I-\nabla u,P\rangle-\int_{\Omega}f\cdot u\, dx. \end{align*}\]

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *
import ipywidgets as widgets


def GenerateMesh(t=0.1):
    bar = MoveTo(0, -t / 2).Rectangle(1, t).Face()
    bar.edges.Min(X).name = "left"
    bar.edges.Max(X).name = "right"
    bar.edges.Min(Y).name = "bottom"
    bar.edges.Max(Y).name = "top"
    bar.edges.Min(X).maxh = t / 8
    mesh = Mesh(OCCGeometry(bar, dim=2).GenerateMesh(maxh=t / 2.5))
    return mesh


mesh = GenerateMesh()
Draw(mesh)

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


def NeoHooke(F):
    C = F.trans * F
    return mu / 2 * (Trace(C - Id(2)) - log(Det(C))) + lam / 2 * (Det(F) - 1) ** 2


def SkewT(v):
    if v.dim == 3:
        return 0.5 * CoefficientFunction(
            (0, -v[2], v[1], v[2], 0, -v[0], -v[1], v[0], 0), dims=(3, 3)
        )
    else:
        return 0.5 * CoefficientFunction((0, -v, v, 0), dims=(2, 2))

order = 3
par = Parameter(0)
force = CF((0, -5 * par))

# clamped boundary at left, free elsewhere
V = HCurl(mesh, order=order, dirichlet="left")
P_sym = HDivDiv(mesh, order=order, dirichlet="bottom|right|top")
P_skw = L2(mesh, order=order)
F = MatrixValued(L2(mesh, order=order))
fes = V * P_sym * P_skw * F
gf_solution = GridFunction(fes)

u, Psym, Pskw, F = fes.TrialFunction()

P = Psym + CoefficientFunction((0, Pskw, -Pskw, 0), dims=(2, 2))
n = specialcf.normal(2)

a = BilinearForm(fes, symmetric=True, condense=True)
a += Variation((NeoHooke(F)) * dx).Compile()
a += Variation(InnerProduct(F - Grad(u) - Id(2), P) * dx).Compile()
a += Variation(u * n * Psym[n, n] * dx(element_boundary=True)).Compile()
a += Variation(-(force * u.Trace()) * ds("right")).Compile()

gf_u, gf_Psym, gf_Pskw, gf_F = gf_solution.components

# F=I
gf_F.Set(Id(2))
# set boundary data
gf_Psym.Set((force * n) * OuterProduct(n, n), definedon=mesh.Boundaries("right"))


scene = Draw(
    gf_u,
    mesh,
    deformation=True,
    settings={
        "camera": {
            "transformations": [{"type": "move", "dir": (0, 0, 1), "dist": -2.5}]
        }
    },
)
tw = widgets.Text(value="step = 0")
display(tw)
gf_history = GridFunction(fes, multidim=0)

num_steps = 20
with TaskManager():
    for step in range(num_steps):
        par.Set((step + 1) / num_steps)
        solvers.Newton(
            a,
            gf_solution,
            maxerr=1e-8,
            printing=False,
            inverse="",
        )
        scene.Redraw()
        tw.value = f"step = {step+1}/{num_steps}"
        gf_history.AddMultiDimComponent(gf_solution.vec)

Draw(
    gf_history.components[0],
    mesh,
    animate=True,
    min=0,
    max=1,
    autoscale=True,
    deformation=True,
    settings={
        "camera": {
            "transformations": [{"type": "move", "dir": (0, 0, 1), "dist": -2.5}]
        }
    },
);

Noting that \(\mathrm{skew} \nabla u= \mathrm{skw}(\mathrm{curl}(u))\in L^2\) there is no need to lift the complete deformation gradient, only its symmetric part, \(F_{\mathrm{sym}}\). Then also only the symmetric part of \(P\) is needed. Further, we use hybridization to retain a minimization problem.

Find \(u\), \(F_{\mathrm{sym}}\), \(P_{\mathrm{sym}}\), and \(\alpha\) for the Lagrangian

\[\begin{align*} \mathcal{L}(u, F_{\mathrm{sym}},P_{\mathrm{sym}}, \alpha) = &\ \int_{\Omega}\Psi(F_{\mathrm{sym}} + \mathrm{skw}(\mathrm{curl}(u)))\,dx + \langle \varepsilon (u), P_{\mathrm{sym}}\rangle - \int_\Omega (F_{\mathrm{sym}}-I):P_{\mathrm{sym}}\, dx+ \sum_{T} \int_{\partial T} P_{\mathrm{sym},nn} \alpha_{n}\, ds. \end{align*}\]

from ngsolve import *
from ngsolve.webgui import Draw
import ipywidgets as widgets
from ngsolve.meshes import MakeQuadMesh

L_1 = 48  # 48m
L_2 = 44  # 44m
L_3 = 16  # 16m


def GenerateMesh(num_el=2):
    mesh = MakeQuadMesh(
        nx=num_el,
        ny=num_el,
        mapping=lambda x, y: (L_1 * x, L_2 * x + L_2 * y - (L_2 - L_3) * x * y),
    )
    return mesh


mesh = GenerateMesh()

mu = Parameter(80.194)  # N/mm^2
lam = Parameter(400889.8)  # N/mm^2

par = Parameter(0)
force = CoefficientFunction((0, par * 32))


def NeoHooke(F):
    C = F.trans * F
    return mu / 2 * (Trace(C - Id(2)) - log(Det(C))) + lam / 2 * (Det(F) - 1) ** 2


def SkewT(v):
    if v.dim == 3:
        return 0.5 * CoefficientFunction(
            (0, -v[2], v[1], v[2], 0, -v[0], -v[1], v[0], 0), dims=(3, 3)
        )
    else:
        return 0.5 * CoefficientFunction((0, -v, v, 0), dims=(2, 2))


order = 3

# clamped at left, free elsewhere
V = HCurl(mesh, order=order, dirichlet="left")
Sigma = Discontinuous(HDivDiv(mesh, order=order))
Hyb = NormalFacetFESpace(mesh, order=order, dirichlet="left")
Gamma = MatrixValued(L2(mesh, order=order + 1), symmetric=True)
fes = V * Sigma * Gamma * Hyb
gf_solution = GridFunction(fes)

u, Psym, Fsym, uh = fes.TrialFunction()
Fmat = Fsym + SkewT(curl(u))

n = specialcf.normal(2)

a = BilinearForm(fes, symmetric=True, condense=True)
a += Variation((NeoHooke(Fmat)) * dx).Compile()
a += Variation(InnerProduct(Fsym - Grad(u) - Id(2), Psym) * dx).Compile()
a += Variation((u - uh) * n * Psym[n, n] * dx(element_boundary=True)).Compile()
a += Variation(-(force * uh.Trace() + force * u.Trace()) * ds("right")).Compile()

gf_u, gf_Psym, gf_Fsym, gf_uh = gf_solution.components
gf_F = gf_Fsym + SkewT(curl(gf_u))

# F=I
gf_Fsym.Set(Id(2))

scene = Draw(
    gf_u,
    mesh,
    deformation=True,
    settings={
        "camera": {
            "transformations": [{"type": "move", "dir": (0, 0, 1), "dist": -1.5}]
        }
    },
)
tw = widgets.Text(value="step = 0")
display(tw)
gf_history = GridFunction(fes, multidim=0)

num_steps = 40
with TaskManager():
    for step in range(num_steps):
        par.Set((step + 1) / num_steps)
        solvers.Newton(
            a,
            gf_solution,
            maxerr=1e-8,
            printing=False,
            inverse="sparsecholesky",
        )
        scene.Redraw()
        tw.value = f"step = {step+1}/{num_steps}"
        gf_history.AddMultiDimComponent(gf_solution.vec)

Draw(
    gf_history.components[0],
    mesh,
    animate=True,
    min=0,
    max=25,
    autoscale=True,
    deformation=True,
    settings={
        "camera": {
            "transformations": [{"type": "move", "dir": (0, 0, 1), "dist": -1.5}]
        }
    },
);

Alternative TDNNS formulations for nonlinear elasticity:

Use Cauchy-Green strain tensor \(C=F^\top F\) instead of \(F\) as additional strain field. Advantage: symmetric field independent of rigid body motions. Challenge: Handling nonlinear distribution \(\nabla u^\top\nabla u\) [Neunteufel, Pechstein, Schöberl. Three-field mixed finite element methods for nonlinear elasticity. Computer Methods in Applied Mechanics and Engineering, (2021).]
Use Updated Lagrangian scheme for linear TDNNS formulation [Pechstein. Large deformation mixed finite elements for smart structures. Mechanics of Advanced Materials and Structures, (2020).]