Ctrl+K

Linear Kirchhoff-Love and Reissner-Mindlin shells

9.2. Linear Kirchhoff-Love and Reissner-Mindlin shells#

In this section we discretize the linear Kirchhoff-Love shell equations with the Hellan-Herrmann-Johnson (HHJ) method and the linear Reissner-Mindlin shell equation by adding shearing degrees of freedom on top of the HHJ method in terms of the tangential-displacement-normal-normal-stress (TDNNS) method. Further, we will observe and discuss the problem of membrane locking of shells and how to cure it by using an interpolation operator into the so-called Regge finite elements following [Neunteufel, Schöberl. Avoiding membrane locking with Regge interpolation. Comput. Methods Appl. Mech. Eng. 373, (2021)].

9.2.1. HHJ for Kirchhoff-Love shells#

We show how linear Koiter (Kirchhoff-Love) shells can be handled with the Hellan-Herrmann-Johnson method and discuss the effect of membrane locking. The Kirchhoff-Love shell energy consists of the membrane and bending term and reads

\[\begin{align*} \mathcal{W}_{\mathrm{KL}}(u)=\int_{\mathcal{S}}\Big(\frac{t}{2}\|\mathrm{sym}(\nabla_{\mathcal{S}}^{\mathrm{cov}}u)\|^2_{\mathbb{C}} + \frac{t^3}{24}\Big\|\sum_{i=1}^3\nabla_{\mathcal{S}}^2u_i \nu_i\Big\|^2_{\mathbb{C}}\Big)\,ds - \int_{\mathcal{S}}f\cdot u\,ds, \end{align*}\]

where \(\nabla_{\mathcal{S}}^{\mathrm{cov}}u=P_{\mathcal{S}}\nabla_{\mathcal{S}}u\) is the covariant surface derivative, \(t\) the shell thickness, \(\nabla_{\mathcal{S}}^2u_i\) the surface Hessian, \(\nu\) the unit normal vector on \(\mathcal{S}\), and \(\mathbb{C}\) the material tensor. We use the HHJ method to rewrite the bending energy term, which would lead to a fourth order problem, into a mixed second order saddle-point problem. To this end, we introduce the bending moment stress tensor

\[\begin{align*} \sigma =\mathbb{C}(\mathcal{H}_{\nu}),\qquad \mathcal{H}_{\nu}:=\sum_{i=1}^3\nabla_{\mathcal{S}}^2u_i \nu_i \end{align*}\]

yielding the Lagrangian [Neunteufel, Schöberl. The Hellan-Herrmann-Johnson and TDNNS method for linear and nonlinear shells. arXiv, (2023)].

\[\begin{align*} \mathcal{L}_{\mathrm{KL}}^{\mathrm{HHJ}}(u,\sigma)= \int_{\mathcal{S}}\Big(\frac{t}{2}\|\mathrm{sym}(\nabla_{\mathcal{S}}^{\mathrm{cov}}u)\|^2_{\mathbb{C}} -\frac{6}{t^3}\|\sigma\|_{\mathcal{C}^{-1}}\Big)\,ds + \sum_{T\in\mathcal{T}}\int_T \sigma:\mathcal{H}_{\nu}\,ds - \sum_{E\in\mathcal{E}}\int_E [\![(\nabla_{\mathcal{S}}u)_{\mu\nu}]\!]\sigma_{\mu\mu}\,dl- \int_{\mathcal{S}}f\cdot u\,ds, \end{align*}\]

where \(\mu\) is the co-normal vector at the element-interfaces. Note that the duality pairing

\[\begin{align*} \langle \sigma, \mathcal{H}_{\nu}\rangle := \sum_{T\in\mathcal{T}}\int_T \sigma:\mathcal{H}_{\nu}\,ds - \sum_{E\in\mathcal{E}}\int_E [\![(\nabla_{\mathcal{S}}u)_{\mu\nu}]\!]\sigma_{\mu\mu}\,dl \end{align*}\]

reduces for plates to the HHJ pairing (noting that \(\nu=e_3\), \(\mu=n\), and \(w=u_3\))

\[\begin{align*} \langle \sigma, \mathcal{H}_{\nu}\rangle = \sum_{T\in\mathcal{T}}\int_T \sigma:\nabla^2w\,ds - \sum_{E\in\mathcal{E}}\int_E [\![\partial_{n}w]\!]\sigma_{nn}\,dl. \end{align*}\]

We consider a hyperboloid with free ends as a benchmark example. The shell described by the equation

\[\begin{align*} x^2+y^2 = 1+z^2,\qquad z\in [-1,1] \end{align*}\]

is loaded by a periodic volume force \(p=t^310^4\cos(2\zeta)\nu\), \(\zeta\in [0,2\pi)\). Due to symmetry of the problem, we only consider one eight of the geometry and prescribe symmetry boundary conditions on the new boundaries. The material parameters read

\[\begin{align*} \bar{E} = 2.85\times 10^4,\quad \bar{\nu}=0.3. \end{align*}\]

We consider as reference value the radial deflection at point \((0,0,1)\) which read for different thicknesses \(t=0.1\): \(-0.1856305\), \(t=0.01\): \(-0.1502913\), \(t=0.001\): \(-0.1498749\).

from ngsolve import *
import netgen.meshing as meshing
from ngsolve.webgui import Draw

ref_value = {0.1: -0.1856305, 0.01: -0.1502913, 0.001: -0.1498749}
thickness = 0.001
radius = 1
order = 2

E = 2.85e4
nu = 0.3
mu = E / (2 * (1 + nu))
lam = E * nu / (1 - nu**2)


free = "right"
symmetry = "left|top|bottom"

e_n = -specialcf.normal(3)
force = 1e4 * thickness * cos(2 * IfPos(y, atan(z / y), pi / 2)) * e_n * thickness**2

mapping = lambda x, y, z: (
    x,
    sqrt(1 + x**2) * cos(pi / 2 * y),
    sqrt(1 + x**2) * sin(pi / 2 * y),
)
geom = meshing.SurfaceGeometry(mapping)
mesh = Mesh(geom.GenerateMesh(quads=False, nx=5, ny=5)).Curve(order)
Draw(mesh)

# point of interest
P = (0, 0, radius)

For implementation we use a hybridization technique to eliminate \(\sigma\) on element-level. This is done exactly as for Kirchhoff-Love plates.

First we try out this method without further adaptions. We observe that for small thicknesses \(t\) the results are bad, the displacements \(u\) tend to be zero. This so-called membrane locking phenomena is caused by the limit constraint that the membrane energy has to be zero for \(t\to0\) (in the so-called bending dominated regime). This constraint cannot be fulfilled good enough with the finite element displacements. Thus, to fulfill the constraint, the displacements get close to zero.

To circumvent membrane locking we interpolate the membrane strains into a finite element space of lower polynomial degree. This will give the finite element solution more freedom to fulfill the “weakened” constraints without being zero. However, projecting too much allows for so-called spurious zero-energy modes making the method unstable in the membrane dominated regime (where no membrane locking occurs as membrane modes are dominant). A good balance is to interpolate into the tangential-tangential continuous Regge finite element space of one polynomial order less than the displacement field. We therefore replace the membrane energy

\[\begin{align*} \frac{t}{2}\|\mathrm{sym}(\nabla_{\mathcal{S}}^{\mathrm{cov}} u)\|_{\mathbb{C}}^2\qquad \Rightarrow\qquad\frac{t}{2}\|\mathcal{I}^{\mathrm{Reg}}\mathrm{sym}(\nabla_{\mathcal{S}}^{\mathrm{cov}} u)\|_{\mathbb{C}}^2. \end{align*}\]
# try out False and True
interpolateMembrane = False


def MaterialStress(mat):
    return E / (1 - nu**2) * ((1 - nu) * mat + nu * Trace(mat) * Id(3))


def MaterialStressInv(mat):
    return (1 + nu) / E * (mat - nu / (nu + 1) * Trace(mat) * Id(3))


fesU = VectorH1(
    mesh,
    order=order,
    dirichletx_bbnd="left",
    dirichlety_bbnd="top",
    dirichletz_bbnd="bottom",
)
fesS = HDivDivSurface(mesh, order=order - 1, discontinuous=True)
fesH = NormalFacetSurface(mesh, order=order - 1, dirichlet_bbnd="left|top|bottom")

fes = fesU * fesS * fesH
(u, sigma, uh), (du, dsigma, duh) = fes.TnT()
# Need to take traces as we are on the surface
sigma, dsigma, uh, duh = sigma.Trace(), dsigma.Trace(), uh.Trace(), duh.Trace()

fesRegge = HCurlCurl(mesh, order=order - 1)

# normal, tangential, and co-normal vectors
nv = specialcf.normal(3)
tv = specialcf.tangential(3)
cnv = Cross(nv, tv)

# Projection to tangent space and surface derivatives
Ptau = Id(3) - OuterProduct(nv, nv)
gradu = Grad(u).Trace()
graddu = Grad(du).Trace()

gf_solution = GridFunction(fes)

# sum_{i=1}^3 hesse(u_i) nv_i
Hn = (u.Operator("hesseboundary").trans * nv).Reshape((3, 3))
dHn = (du.Operator("hesseboundary").trans * nv).Reshape((3, 3))

bfa = BilinearForm(fes, symmetric=True, condense=True)
# membrane part
if interpolateMembrane:
    bfa += (
        thickness
        * InnerProduct(
            MaterialStress(Interpolate(Sym(Ptau * gradu), fesRegge)),
            Interpolate(Sym(Ptau * graddu), fesRegge),
        )
    ) * ds
else:
    bfa += (
        thickness * InnerProduct(MaterialStress(Sym(Ptau * gradu)), Sym(Ptau * graddu))
    ) * ds
# bending part
bfa += (-12 / thickness**3 * InnerProduct(MaterialStressInv(sigma), dsigma)) * ds
bfa += (InnerProduct(dsigma, Hn) + InnerProduct(sigma, dHn)) * ds
bfa += (
    sigma[cnv, cnv] * (duh - graddu[nv, :]) * cnv
    + dsigma[cnv, cnv] * (uh - gradu[nv, :]) * cnv
) * ds(element_boundary=True)


f = LinearForm(force * du * ds)

with TaskManager():
    bfa.Assemble()
    f.Assemble()

    if bfa.condense:
        f.vec.data += bfa.harmonic_extension_trans * f.vec
        inv = bfa.mat.Inverse(fes.FreeDofs(True), inverse="sparsecholesky")
        gf_solution.vec.data = inv * f.vec
        gf_solution.vec.data += bfa.harmonic_extension * gf_solution.vec
        gf_solution.vec.data += bfa.inner_solve * f.vec
    else:
        inv = bfa.mat.Inverse(fes.FreeDofs(), inverse="")
        gf_solution.vec.data = inv * f.vec

gf_u, gf_sigma, _ = gf_solution.components
print(f"value = {gf_u(mesh(*P, BND))[2]}, reference = {ref_value[thickness]}")
Draw(gf_u, mesh, deformation=True);

value = -0.0024043926800899175, reference = -0.1498749

9.2.2. TDNNS method for Reissner-Mindlin shells#

Next, we add shearing degrees of freedom to allow for shear deformations. The additional shearing field \(\gamma\) will be discretized by \(H\)(curl)-conforming Nedelec elements such that we obtain a TDNNS like discretization method.

The primal energy including shearing energy would read

\[\begin{align*} \mathcal{W}_{\mathrm{RM}}(u,\gamma)=\int_{\mathcal{S}}\Big(\frac{t}{2}\|\mathrm{sym}(\nabla_{\mathcal{S}}^{\mathrm{cov}}u)\|^2_{\mathbb{C}} + \frac{t^3}{24}\Big\|\sum_{i=1}^3\nabla_{\mathcal{S}}^2u_i \nu_i-\nabla_{\mathcal{S}}\gamma\Big\|^2_{\mathbb{C}} + \frac{t\kappa G}{2}\|\gamma\|_2\Big)\,ds - \int_{\mathcal{S}}f\cdot u\,ds, \end{align*}\]

where \(G=\frac{E}{2(1+\nu)}\) and \(\kappa=5/6\) denote the shearing modulus and shear correction factor, respectively.

We again add the bending moment tensor as unknown giving us the following Lagrangian [Neunteufel, Schöberl. The Hellan-Herrmann-Johnson and TDNNS method for linear and nonlinear shells. arXiv, (2023).]

\[\begin{align*} \mathcal{L}_{\mathrm{RM}}^{\mathrm{TDNNS}}(u,\sigma,\gamma)= \int_{\mathcal{S}}\Big(\frac{t}{2}\|\mathrm{sym}(\nabla_{\mathcal{S}}^{\mathrm{cov}}u)\|^2_{\mathbb{C}} -\frac{6}{t^3}\|\sigma\|_{\mathcal{C}^{-1}}\Big)\,ds + \sum_{T\in\mathcal{T}}\int_T \sigma:(\mathcal{H}_{\nu}-\nabla_{\mathcal{S}}\gamma)\,ds - \sum_{E\in\mathcal{E}}\int_E [\![(\nabla_{\mathcal{S}}u)_{\mu\nu}-\gamma_{\mu}]\!]\sigma_{\mu\mu}\,dl- \int_{\mathcal{S}}f\cdot u\,ds. \end{align*}\]

To avoid possible membrane locking we can again add the interpolant into Regge finite elements in the membrane energy term. Due to the hierarchical approach of adding the shearing field \(\gamma\) by construction no shear locking can occur, as the Kirchhoff-Love limit \(\gamma=0\) can always exactly be represented.

The reference value of the hyperboloid benchmark now change to \(t=0.1\): \(-0.18954566\), \(t=0.01\): \(-0.15046617\), \(t=0.001\): \(-0.1498902\).

from ngsolve import *
import netgen.meshing as meshing
from ngsolve.webgui import Draw

ref_value = {0.1: -0.18954566, 0.01: -0.15046617, 0.001: -0.1498902}
thickness = 0.001
radius = 1
order = 2

E = 2.85e4
nu = 0.3
mu = E / (2 * (1 + nu))
lam = E * nu / (1 - nu**2)
G = E / (2 * (1 + nu))
kappa = 5 / 6


free = "right"
symmetry = "left|top|bottom"

e_n = -specialcf.normal(3)
force = 1e4 * thickness * cos(2 * IfPos(y, atan(z / y), pi / 2)) * e_n * thickness**2

mapping = lambda x, y, z: (
    x,
    sqrt(1 + x**2) * cos(pi / 2 * y),
    sqrt(1 + x**2) * sin(pi / 2 * y),
)
geom = meshing.SurfaceGeometry(mapping)
mesh = Mesh(geom.GenerateMesh(quads=False, nx=5, ny=5)).Curve(order)

# point of interest
P = (0, 0, radius)


def MaterialStress(mat):
    return E / (1 - nu**2) * ((1 - nu) * mat + nu * Trace(mat) * Id(3))


def MaterialStressInv(mat):
    return (1 + nu) / E * (mat - nu / (nu + 1) * Trace(mat) * Id(3))


fesU = VectorH1(
    mesh,
    order=order,
    dirichletx_bbnd="left",
    dirichlety_bbnd="top",
    dirichletz_bbnd="bottom",
)
fesB = HCurl(mesh, order=order - 1)
fesS = HDivDivSurface(mesh, order=order - 1, discontinuous=True)
fesH = NormalFacetSurface(mesh, order=order - 1, dirichlet_bbnd="left|top|bottom")
fes = fesU * fesB * fesS * fesH

fesRegge = HCurlCurl(mesh, order=order - 1)

(u, gamma, sigma, uh), (du, dgamma, dsigma, duh) = fes.TnT()
# Need to take traces as we are on the surface
sigma, dsigma, uh, duh, gamma, dgamma = (
    sigma.Trace(),
    dsigma.Trace(),
    uh.Trace(),
    duh.Trace(),
    gamma.Trace(),
    dgamma.Trace(),
)

nv = specialcf.normal(3)
tv = specialcf.tangential(3)
cnv = Cross(nv, tv)

Ptau = Id(3) - OuterProduct(nv, nv)
gradu = Grad(u).Trace()
graddu = Grad(du).Trace()

gf_solution = GridFunction(fes)

Hn = (u.Operator("hesseboundary").trans * nv).Reshape((3, 3))
dHn = (du.Operator("hesseboundary").trans * nv).Reshape((3, 3))

bfa = BilinearForm(fes, symmetric=True, condense=True)
# membrane part
bfa += (
    thickness
    * InnerProduct(
        MaterialStress(Interpolate(Sym(Ptau * gradu), fesRegge)),
        Interpolate(Sym(Ptau * graddu), fesRegge),
    )
) * ds
# bending part
bfa += (-12 / (thickness**3) * InnerProduct(MaterialStressInv(sigma), dsigma)) * ds
bfa += (
    InnerProduct(dsigma, Hn - Grad(gamma)) + InnerProduct(sigma, dHn - Grad(dgamma))
) * ds
bfa += (
    sigma[cnv, cnv] * (duh - graddu[nv, :] + dgamma) * cnv
    + dsigma[cnv, cnv] * (uh - gradu[nv, :] + gamma) * cnv
) * ds(element_boundary=True)
# shear part
bfa += (thickness * kappa * G * InnerProduct(gamma, dgamma)) * ds


f = LinearForm(force * du * ds)

with TaskManager():
    bfa.Assemble()
    f.Assemble()

    if bfa.condense:
        f.vec.data += bfa.harmonic_extension_trans * f.vec
        inv = bfa.mat.Inverse(fes.FreeDofs(True), inverse="sparsecholesky")
        gf_solution.vec.data = inv * f.vec
        gf_solution.vec.data += bfa.harmonic_extension * gf_solution.vec
        gf_solution.vec.data += bfa.inner_solve * f.vec
    else:
        inv = bfa.mat.Inverse(fes.FreeDofs(), inverse="")
        gf_solution.vec.data = inv * f.vec

gf_u, gf_beta, gf_sigma, _ = gf_solution.components
print(f"value = {gf_u(mesh(*P, BND))[2]}, reference = {ref_value[thickness]}")
Draw(gf_u, mesh, deformation=True);

value = -0.14998434211853243, reference = -0.1498902
Contents