Kirchhoff-Love plates with the Hellan-Herrmann-Johnson method

10. Kirchhoff-Love plates with the Hellan-Herrmann-Johnson method#

The Kirchhoff-Love plate is the shear-free limit \(t\to0\) of Reissner-Mindlin plates. It has only one displacement unknown \(w\), and thus no shear locking can occur. However, it is a fourth-order problem making the discretization more challenging. The Hellan-Herrmann-Johnson (HHJ) method [Comodi. The Hellan-Herrmann-Johnson method: some new error estimates and postprocessing. Math. Comp. 52 (1989), 17-29] rewrites it as a mixed second-order saddle-point problem by introducing the bending moment as an additional unknown.

10.1. Kirchhoff-Love equation#

The Kirchhoff-Love problem is: Find \(w\in H^2_0(\Omega)\) such that

\[ \int_\Omega \mathbb{D}\nabla^2 w:\nabla^2 v\,dx = \int_\Omega f v\,dx \qquad \forall v\in H^2_0(\Omega). \]

The strong form with clamped boundary conditions is

\[ \operatorname{div}\operatorname{div}(\mathbb{D}\nabla^2w)=f \quad \text{in } \Omega,\qquad w=\partial_n w=0 \quad \text{on } \partial \Omega. \]

A conforming Galerkin method would require \(H^2\)-conforming finite elements. HHJ avoids this by introducing the bending moment as an additional unknown yielding a mixed method

\[ \sigma = \mathbb{D}\nabla^2 w. \]

10.2. HHJ mixed formulation#

The moment tensor lives in

\[ H(\operatorname{div}\operatorname{div}) = \{\tau\in L^2(\Omega;\mathbb{S}) : \operatorname{div}\operatorname{div}\tau\in H^{-1}(\Omega)\}. \]

The mixed formulation for the clamped case is: Find \((\sigma,w)\in H(\operatorname{div}\operatorname{div})\times H^1_0\) such that for all \((\tau,v)\in H(\operatorname{div}\operatorname{div})\times H^1_0\)

\[\begin{split} \begin{align*} &\int_\Omega \mathbb{D}^{-1}\sigma:\tau\,dx- \langle \tau,\nabla^2 w\rangle &&= 0,\\ &-\langle \sigma,\nabla^2 v\rangle &&= -\int_\Omega f v\,dx. \end{align*} \end{split}\]

The duality pairing is the same as the TDNNS pairing setting \(\beta =\nabla w\) evaluated element-wise

\[ \langle\sigma,\nabla^2 w\rangle = \sum_T\left(\int_T \sigma:\nabla^2w\,dx - \int_{\partial T}\sigma_{nn}\frac{\partial w}{\partial n}\,ds\right). \]

By integrating once and twice by parts, the pairing reads

\[\begin{split} \begin{align*} \langle\sigma,\nabla^2 w\rangle &= -\sum_T\left(\int_T \mathrm{div}(\sigma)\cdot\nabla w\,dx - \int_{\partial T}\sigma_{nt}\frac{\partial w}{\partial t}\,ds\right)\\ &=\sum_T\left(\int_T \mathrm{div}\mathrm{div}(\sigma) w\,dx - \int_{\partial T}(\mathrm{div}(\sigma)\cdot n + \frac{\sigma_{nt}}{\partial t})w\,ds + \sum_{V\subset T}[\![\sigma_{nt}]\!]_V^T\, w(V)\right). \end{align*} \end{split}\]

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw

L = 0.5
wp = WorkPlane().MoveTo(-0.5 * L, -0.5 * L)
for i in range(4):
    wp.Rotate(-90).Line(L).Rotate(90).Line(L / 4).Rotate(90).Line(L)
shape = wp.Face()
shape.edges.name = "free"
shape.edges.Max(X).name = "force"
shape.edges.Max(Y).name = "clamped"
shape.edges.Min(Y).name = "clamped"

mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.05))
Draw(mesh);

E = 2e3
nu = 0.3
thickness = 0.1
force = 1


def DMatInv(mat):
    return 12 * (1 - nu**2) / (thickness**3 * E) * (
        1 / (1 - nu) * mat - nu / (1 - nu**2) * Trace(mat) * Id(2)
    )


def tang(u):
    n = specialcf.normal(2)
    return u - (u * n) * n

def KirchhoffLovePlateHHJ(order=1):
    V = HDivDiv(mesh, order=order - 1, dirichlet="free|force")
    Q = H1(mesh, order=order, dirichlet="clamped")
    X = V * Q
    (sigma, w), (tau, v) = X.TnT()
    n = specialcf.normal(2)

    a = BilinearForm(X, symmetric=True)
    a += (
        InnerProduct(DMatInv(sigma), tau)
        + div(sigma) * Grad(v)
        + div(tau) * Grad(w)
    ) * dx
    a += -(
        sigma[n, :] * tang(Grad(v))
        + tau[n, :] * tang(Grad(w))
    ) * dx(element_boundary=True)
    a.Assemble()

    f = LinearForm(v * force * ds("force")).Assemble()
    gf_solution = GridFunction(X)
    gf_solution.vec.data = a.mat.Inverse(X.FreeDofs(), inverse="umfpack") * f.vec
    return gf_solution


gf_solution = KirchhoffLovePlateHHJ(order=2)
gf_sigma, gf_w = gf_solution.components
Draw(gf_w, mesh, name="displacement", deformation=True, euler_angles=[-60, 5, 30]);

10.3. Hybridization#

The direct HHJ formulation is a saddle-point problem. Hybridization breaks the normal-normal continuity of \(\sigma\) locally and restores it with a facet multiplier \(\alpha\)

\[ \alpha \approx \partial_n w. \]

This allows static condensation. After eliminating the local moment unknowns one obtains a symmetric positive definite system for the remaining global unknowns.

def KirchhoffLovePlateHHJ_Hybridized(order=1):
    V = Discontinuous(HDivDiv(mesh, order=order - 1))
    Q = H1(mesh, order=order, dirichlet="clamped")
    H = NormalFacetFESpace(mesh, order=order - 1, dirichlet="clamped")
    X = V * Q * H
    (sigma, w, alpha), (tau, v, beta) = X.TnT()
    n = specialcf.normal(2)

    a = BilinearForm(X, symmetric=True, condense=True)
    a += (
        InnerProduct(DMatInv(sigma), tau)
        + div(sigma) * Grad(v)
        + div(tau) * Grad(w)
    ) * dx
    a += -(
        sigma[n, :] * tang(Grad(v))
        + tau[n, :] * tang(Grad(w))
    ) * dx(element_boundary=True)
    a += (
        sigma[n, n] * beta * n
        + tau[n, n] * alpha * n
    ) * dx(element_boundary=True)
    a.Assemble()

    f = LinearForm(v * force * ds("force")).Assemble()
    gf_solution = GridFunction(X)

    f.vec.data += a.harmonic_extension_trans * f.vec
    gf_solution.vec.data += a.mat.Inverse(X.FreeDofs(True), inverse="sparsecholesky") * f.vec
    gf_solution.vec.data += a.harmonic_extension * gf_solution.vec
    gf_solution.vec.data += a.inner_solve * f.vec
    return gf_solution

order = 1
gf_hyb = KirchhoffLovePlateHHJ_Hybridized(order=order)
gf_sigma_hyb, gf_w_hyb, gf_alpha_hyb = gf_hyb.components
Draw(gf_w_hyb, mesh, name="hybridized displacement", deformation=True, euler_angles=[-60, 5, 30]);

10.4. Postprocessing property of HHJ method#

We can use lowest order elements, i.e. linear polynomials for the vertical deflection \(w\) and piece-wise constants for the moment tensor \(\sigma\). As \(\sigma= \mathbb{D}\nabla^2 w\) one might try to compute a new quadratic \(\tilde w\) as a postprocessing step element-wise with an improved order of convergence.

[Li. Recovery-based a posteriori error analysis for plate bending problems. J Sci Comput. 2021]

Correct displacement \(w \in V_h^k\) with bubble functions of one degree higher \(w_b\in V_h^{k+1}\) such that \(\tilde{w} = w+w_b\)

\[ \begin{align*} \widetilde w = \operatorname{arg}\min_{v_h \in V_h^{k+1} \atop v_h|_{V_h^k} = w_h} \| \mathbb{D}\nabla^2 v_h - \sigma \|_{L_2}^2. \end{align*} \]

We need to solve a global problem. However, it was proved that the condition number of a Jacobi preconditioner stays bounded. Therefore, the costs are comparable to solve local problems.

from ngsolve.krylovspace import CGSolver

X = H1(mesh, order=order + 1, dirichlet="clamped")
w, dw = X.TnT()

freedofs = BitArray([False] * X.ndof)
for el in mesh.edges:
    freedofs[X.GetDofNrs(el)[order - 1]] = True
freedofs &= ~X.GetDofs(mesh.Boundaries("clamped"))


a = BilinearForm(
    InnerProduct(w.Operator("hesse"), dw.Operator("hesse")) * dx, symmetric=True
).Assemble()
f = LinearForm(
    InnerProduct(
        DMatInv(gf_sigma) - gf_w.Operator("hesse"), dw.Operator("hesse")
    )
    * dx
).Assemble()

gf_wb = GridFunction(X)

preJpoint = a.mat.CreateSmoother(freedofs)
gf_wb.vec.data = CGSolver(a.mat, pre=preJpoint, printrates="", maxiter=100) * f.vec

Draw(gf_wb, mesh, "bubble_correction", order=3)
Draw(
    gf_w + gf_wb,
    mesh,
    "corrected_displacement",
    order=3,
    deformation=True,
    euler_angles=[-60, 5, 30],
);

Alternative post-processing: [Stenberg. Postprocessing schemes for some mixed finite elements. ESAIM: M2AN 25 , 1 (1991), 151-167]

Solve the following problem element-wise

\[ \begin{align*} \widetilde w = \operatorname{arg}\min_{v_h \in P^{k+1} \atop v_h(V) = w_h(V)} \| \mathbb{D}\nabla^2 v_h - \sigma \|_{L_2}^2 \end{align*} \]

under the constraint that the vertex values are preserved. Therefore, a discontinuous Lagrange space can be used, where on each element the dofs corresponding to a vertex are fixed.

Advantage: Local problems
Disadvantage: Corrected displacement are in general discontinuous

10.5. Plate eigenvalue#

We use the HHJ method to compute plate vibration modes. The Kirchhoff-Love eigenvalue problem reads

\[ \operatorname{div}\operatorname{div}(\mathbb{D}\nabla^2 w) = \lambda w. \]

We use the hybridized HHJ method:

\(w_h\in H^1\) stores the transverse displacement,
a normal-facet variable stores the normal derivative trace,
a discontinuous \(H(\operatorname{div}\operatorname{div})\) variable stores the local bending moment,
static condensation eliminates local moment unknowns,
LOBPCG computes eigenvectors of the generalized problem

\[ A u = \lambda M u. \]

import ngsolve.solvers


def MakePlateEVPMesh(maxh=0.35):
    shape = MoveTo(0, 0).RectangleC(5, 5).Face()
    cvert = Vertex(Pnt((0, 0, 0)))
    cvert.name = "support"
    shape = Glue([shape, cvert])

    mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=maxh))
    return mesh


def SolvePlateEVP(order=2, maxh=0.35, num=12, maxit=30):
    mesh = MakePlateEVPMesh(maxh=maxh)

    Vw = H1(mesh, order=order, dirichlet_bbnd=".*")
    Vn = NormalFacetFESpace(mesh, order=order - 1)
    Vs = PrivateSpace(Discontinuous(HDivDiv(mesh, order=order - 1)))
    X = Vw * Vn * Vs

    (w, wn, sigma), (dw, dwn, dsigma) = X.TnT()

    n = specialcf.normal(2)

    def tang(u):
        return u - (u * n) * n

    def normal(u):
        return (u * n) * n

    dS = dx(element_boundary=True)

    a = BilinearForm(X, condense=True, symmetric=True)
    a += -InnerProduct(sigma, dsigma) * dx
    a += (grad(w) * div(dsigma) + grad(dw) * div(sigma)) * dx
    a += -((tang(grad(w)) - normal(wn)) * (dsigma * n)\
           +(tang(grad(dw)) - normal(dwn)) * (sigma * n))* dS
    
    a += 1e-3 * w * dw * dx
    a.Assemble()

    m = BilinearForm(w * dw * dx, symmetric=True).Assemble()

    pre = a.mat.Inverse(freedofs=X.FreeDofs(), inverse="sparsecholesky")
    evals, evecs = ngsolve.solvers.LOBPCG(
        a.mat,
        m.mat,
        pre=pre,
        num=num,
        maxit=maxit,
        printrates=False,
    )

    modes = GridFunction(X, multidim=0)
    for evec in evecs[3:]:
        modes.AddMultiDimComponent(evec)

    return mesh, X, evals, modes


mesh_evp, fes_evp, evals_evp, modes_evp = SolvePlateEVP(order=2, maxh=0.35, num=100, maxit=20)
print("first computed eigenvalues:")
for i, lam in enumerate(evals_evp[:8]):
    print(f"  {i:2d}: {lam:.6e}")

Draw(
    modes_evp.components[0],
    mesh_evp,
    name="plate eigenmodes",
    min=-0.1,
    max=0.1,
    deformation=True,
    scale=1,
    animate=True,
);

first computed eigenvalues:
4.769740e-01
4.769797e-01
8.019030e-01
9.035002e-01
2.477452e+00
4.617966e+00
4.618104e+00
9.542445e+00