16. TDNNS method for Reissner-Mindlin shells

First, we considered linear Kirchhoff-Love shells: membrane plus bending, but no transverse shear. We used the Hellan-Herrmann-Johnson (HHJ) method to handle the bending term through a bending moment tensor.

We now move to linear Reissner-Mindlin shells by adding shear degrees of freedom on top of the HHJ formulation. This leads to the tangential-displacement-normal-normal-stress (TDNNS) method.

The additional shear field \(\gamma\) is discretized by \(H(\operatorname{curl})\)-conforming Nédélec elements, yielding a TDNNS discretization.

The primal energy including shear 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_{\hat{\mathcal{S}}}^2u_i \hat\nu_i-\nabla_{{\mathcal{S}}}\gamma\Big\|^2_{\mathbb{C}} + \frac{t\kappa G}{2}\|\gamma\|_2^2\Big)\,d s - \int_{{\mathcal{S}}} f\cdot u\,d s, \end{align*} \]

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

Rewriting with the rotation \(\beta = \nabla_{{\mathcal{S}}} u^T\hat\nu+\gamma\), one obtains the classical formulas for the bending and shear energies

\[ \begin{align*} \frac{t^3}{24}\|\mathrm{sym}(\nabla_{{\mathcal{S}}}^{\mathrm{cov}}\beta)\|^2_{\mathbb{C}},\qquad \frac{t\kappa G}{2}\|\beta - \nabla_{{\mathcal{S}}} u^T\hat\nu\|_2^2. \end{align*} \]

However, the hierarchical approach with the shear degrees of freedom is better suited. In the implementation below we keep the hierarchical shear variable \(\gamma\) directly; the rotation \(\beta\) is only used here to connect with the classical Reissner-Mindlin notation.

We again add the bending moment tensor as unknown giving us the following Lagrangian [Neunteufel, Schöberl. The Hellan-Herrmann-Johnson and TDNNS methods for linear and nonlinear shells. Computers & Structures, 305 (2024).]

\[\begin{split} \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\|^2_{\mathbb{C}^{-1}} + \frac{t\kappa G}{2}\|\gamma\|_2^2\Big)\,d s + \sum_{ T\in{\mathcal{T}}}\int_{ T} \sigma:(\mathcal{H}_{\hat\nu}(u)-\nabla_{{\mathcal{S}}}\gamma)\,d s \\ &\quad- \sum_{ E\in{\mathcal{E}}}\int_{ E} [\![(\nabla_{{\mathcal{S}}}u)_{\hat\mu\hat\nu}-\gamma_{\hat\mu}]\!]\sigma_{\hat\mu\hat\mu}\,d l- \int_{{\mathcal{S}}} f\cdot u\,d s. \end{align*} \end{split}\]

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 shear field \(\gamma\), no shear locking can occur by construction, as the Kirchhoff-Love limit \(\gamma=0\) can always be represented exactly.

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

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

It is loaded by a periodic normal surface load \(\hat p=t^3 10^4\cos(2\zeta)\hat\nu\), \(\zeta\in [0,2\pi)\). Due to symmetry of the problem, we only consider one eighth 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*} \]

The reference values of the hyperboloid benchmark are: \(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

young_modulus = 2.85e4
poisson = 0.3
G = young_modulus / (2 * (1 + poisson))
kappa = 5 / 6


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 young_modulus / (1 - poisson**2) * ((1 - poisson) * mat + poisson * Trace(mat) * Id(3))


def MaterialStressInv(mat):
    return (1 + poisson) / young_modulus * (mat - poisson / (poisson + 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_gamma, 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 = 454044010.1417629, reference = -0.1498902