# 8.4. Reissner-Mindlin plate#

In this section we discretize and solve the Reissner-Mindlin plate (see Reissner-Mindlin and Kirchhoff-Love plates for a derivation) with two different methods and discuss the appearance of shear locking. The total energy consisting of bending and shearing reads with the vertical deflection $$w$$ and the linearized rotation vector $$\beta$$

\begin{align*} \mathcal{W}_{\mathrm{Reissner-Mindlin}}(w,\beta) = \frac{1}{2}\int_{\Omega}\| \varepsilon(\beta) \|_{\mathbb{D}}^2 + \frac{\kappa\,G}{t^2} \| \nabla w - \beta \|^2\,dx - \int_{\Omega}f\cdot w\,dx, \end{align*}

where $$t$$ denotes the thickness of the plate, $$\varepsilon(\cdot)$$ is the symmetric part of the gradient and the elasticity tensor is $$\mathbb{D}\varepsilon=\frac{E}{12(1-\nu^2)}((1-\nu)\varepsilon+\nu\,\mathrm{tr}(\varepsilon)I_{2\times 2})$$ (with $$E$$ Young’s modulus, $$\nu$$ Poisson ratio). Further, $$G= \frac{E}{2(1+\nu)}$$ and $$\kappa=5/6$$ denote the shearing modulus and shear correction factor, respectively.

We are especially interested in positive, but small thicknesses $$t>0$$. Therefore, the shear energy acts as a penalty enforcing the equality $$\beta = \nabla w$$. If there are not enough finite element functions such that the discrete so-called Kirchhoff constraint can be satisfied, we observe bad results known as shear locking.

We consider a benchmark example [Katili. A new discrete Kirchhoff-Mindlin element based on Mindlin-Reissner plate theory and assumed shear strain fields-Part I: An extended DKT element for thick-plate bending analysis. International Journal for Numerical Methods in Engineering 36, 11 (1993), 1859-1883] with exact solution. It consists of a circle domain, with clamped boundary and a uniform body force (e.g. gravity) is acting on it. Young modulus $$E$$, Poisson ratio $$\nu$$, and shear correction factor $$\kappa$$.

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

E, nu, k = 10.92, 0.3, 5 / 6
G = E / (2 * (1 + nu))
fz = -1

def DMat(mat, E, nu):
return E / (12 * (1 - nu**2)) * ((1 - nu) * mat + nu * Trace(mat) * Id(2))

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

R = 5

def GenerateMesh(order, maxh=1):
circ = Circle((0, 0), R).Face()
circ.edges.name = "circ"
return Mesh(OCCGeometry(circ, dim=2).GenerateMesh(maxh=maxh * R / 3)).Curve(order)

def ExactSolution(thickness=0.1):
r = sqrt(x**2 + y**2)
xi = r / R
Db = E * thickness**3 / (12 * (1 - nu**2))

w_ex = (
fz
* thickness**3
* R**4
/ (64 * Db)
* (1 - xi**2)
* ((1 - xi**2) + 8 * (thickness / R) ** 2 / (3 * k * (1 - nu)))
)
return w_ex

mesh = GenerateMesh(order=2, maxh=1)
w_ex = ExactSolution(thickness=0.1)
Draw(w_ex, mesh, deformation=True, euler_angles=[-60, 5, 30]);


## 8.4.1. Simple discretization with Lagrange elements#

First we consider a simple discretization by using Lagrange elements for the displacements $$w$$ and rotations $$\beta$$ of same polynomial degree $$k$$. Then we can directly discretize the Reissner-Mindlin plate equation: Find $$(w,\beta)\in \mathrm{Lag}^k\times [\mathrm{Lag}^k]^2$$ with $$w=\beta=0$$ on $$\partial\Omega$$ such that for all $$(v,\delta)\in \mathrm{Lag}^k\times [\mathrm{Lag}^k]^2$$

\begin{align*} \frac{1}{12}\int_{\Omega}\mathbb{D}\varepsilon(\beta):\varepsilon(\delta)\,dx + \frac{\kappa G}{t^2}\int_{\Omega}(\nabla w-\beta)\cdot(\nabla v-\delta)\,dx = \int_{\Omega}f\,v\,dx. \end{align*}
def SolveReissnerMindlinLagrange(mesh, order=1, thickness=0.1, draw=False):
fesB = VectorH1(mesh, order=order, dirichlet="circ")
fesW = H1(mesh, order=order, dirichlet="circ")
fes = fesB * fesW
(beta, w), (delta, v) = fes.TnT()

a = BilinearForm(fes, symmetric=True)
# bending part
# shearing part
a += k * G / thickness**2 * InnerProduct(Grad(w) - beta, Grad(v) - delta) * dx
a.Assemble()

f = LinearForm(fz * v * dx).Assemble()

gf_solution = GridFunction(fes)
_, gf_w = gf_solution.components

inv = a.mat.Inverse(fes.FreeDofs(), inverse="sparsecholesky")
gf_solution.vec.data = inv * f.vec

if draw:
Draw(gf_w, mesh, "w")

return gf_w, fes


Solve for different thicknesses for a sequence of refined meshes. Different polynomial order can be tested.

results = []
thicknesses = [1, 0.1, 0.01, 0.001]
# try out polynomial order 1 and 2
order = 2

for t in thicknesses:
results.append([])
w_ex = ExactSolution(thickness=t)

for i in range(5):
mesh = GenerateMesh(order=order, maxh=0.5**i)
gfw, fes = SolveReissnerMindlinLagrange(
mesh, order=order, thickness=t, draw=False
)

# relative L2-error of displacement
err = sqrt(Integrate((gfw - w_ex) * (gfw - w_ex), mesh)) / sqrt(
Integrate(w_ex * w_ex, mesh)
)
results[-1].append((fes.ndof, err))


Plot error. We observe immense shear locking for polynomial order $$k=1$$ as the pre-asymptotic regime increases rapidly. Increasing the polynomial order to $$k=2$$ (or $$k=3$$) mitigates the problem, but it is still present.

import matplotlib.pyplot as plt

plt.yscale("log")
plt.xscale("log")
plt.xlabel("ndof")
plt.ylabel("relative error")
for i, result in enumerate(results):
ndof, err = zip(*result)
plt.plot(ndof, err, "-*", label="thickness=" + str(thicknesses[i]))
plt.legend()
plt.show()


## 8.4.2. TDNNS method#

Next, we apply the TDNNS method to discretize the bending term $$\int_{\Omega}\| \varepsilon(\beta) \|_{\mathbb{D}}^2$$ as in [Pechstein and Schöberl. The TDNNS method for Reissner-Mindlin plates. Numerische Mathematik 137, 3 (2017), 713-740].

Instead of using Lagrangian elements, we discretize the rotations $$\beta$$ in the Nedelec space. Thus, all discrete gradients $$\nabla w$$ are in the space of rotations, and we don’t have the locking problem from before. By introducing the linearized bending moment tensor $$\sigma=\mathbb{D}\varepsilon(\beta)$$ we obtain the following formulation.

Find $$(\sigma,\beta, w) \in H(\mathrm{div} \mathrm{div})\times H(\mathrm{curl})\times H^1(\Omega)$$ such that for all $$(\tau,\delta, v) \in H(\mathrm{div} \mathrm{div})\times H(\mathrm{curl})\times H^1(\Omega)$$

\begin{align*} \begin{array}{ccccll} -\int_{\Omega} \mathbb{D}^{-1} \sigma : \tau\,dx & - & \left< \mathrm{div} \tau, \beta \right> & = & 0 \\ -\left< \mathrm{div} \sigma, \delta \right> & + & \frac{\kappa\, G}{t^2} \int_{\Omega} (\nabla w - \beta) (\nabla v - \delta)\,dx & = & \int_{\Omega} f v\,dx, \end{array} \end{align*}

with the TDNNS duality pairing

\begin{align*} \langle \mathrm{div} \boldsymbol{\sigma}, \beta\rangle &= \sum_{T\in\mathcal{T}}\int_T\mathrm{div}\boldsymbol{\sigma}\cdot\beta\,dx -\int_{\partial T}\boldsymbol{\sigma}_{nt}\beta_t\,ds\\ &=-\sum_{T\in\mathcal{T}}\int_T\boldsymbol{\sigma}:\nabla \beta\,dx +\int_{\partial T}\boldsymbol{\sigma}_{nn}\beta_n\,ds =-\langle \boldsymbol{\sigma}, \nabla \beta\rangle \end{align*}

and the inverted elasticity tensor $$\mathbb{D}^{-1}\varepsilon = \frac{12(1-\nu^2)}{Et^3}\big(\frac{1}{1-\nu}\varepsilon-\frac{\nu}{1-\nu^2}\mathrm{tr}(\varepsilon)I_{2\times2}\big)$$. Here $$\partial T$$ denotes the element-boundary (edges) of $$T$$, $$n$$ the normal vector and $$t$$ the tangential vector on $$\partial T$$. With e.g. $$\sigma_{nt}:= t^\top\sigma n$$ we denote the normal tangential component of $$\sigma$$.

We use Lagrangian elements of order $$k$$ for the displacement $$w$$, and one order less for the rotations $$\beta$$ and the bending moment tensor $$\sigma$$.

def SolveReissnerMindlinTDNNS(mesh, order=1, thickness=0.1, draw=False):
fesB = HCurl(mesh, order=order - 1, dirichlet="circ")
fesS = HDivDiv(mesh, order=order - 1, dirichlet="")
fesW = H1(mesh, order=order, dirichlet="circ")

fes = fesW * fesB * fesS
(w, beta, sigma), (v, delta, tau) = fes.TnT()

n = specialcf.normal(2)

a = BilinearForm(fes, symmetric=True)
# bending part
a += (
-InnerProduct(DMatInv(sigma, E, nu), tau)
) * dx
a += -(sigma[n, n] * delta * n + tau[n, n] * beta * n) * dx(element_boundary=True)
# shearing part
a += k * G / thickness**2 * InnerProduct(Grad(w) - beta, Grad(v) - delta) * dx
a.Assemble()

f = LinearForm(fz * v * dx).Assemble()

gf_solution = GridFunction(fes)
gf_w, _, _ = gf_solution.components

inv = a.mat.Inverse(fes.FreeDofs(), inverse="")
gf_solution.vec.data = inv * f.vec
if draw:
Draw(gf_w, mesh, "w")

return gf_w, fes

results = []

# try out polynomial order 1 and 2
order = 1

thicknesses = [1, 0.1, 0.01, 0.001]
for t in thicknesses:
results.append([])
w_ex = ExactSolution(thickness=t)

for i in range(5):
mesh = GenerateMesh(order=order, maxh=0.5**i)
gfw, fes = SolveReissnerMindlinTDNNS(
mesh, order=order, thickness=t, draw=False
)

# relative L2-error of displacement
err = sqrt(Integrate((gfw - w_ex) * (gfw - w_ex), mesh)) / sqrt(
Integrate(w_ex * w_ex, mesh)
)
results[-1].append((fes.ndof, err))


With the TDNNS method we obtain a locking free discretization already for polynomial degree $$k=1$$ for the displacements.

import matplotlib.pyplot as plt

plt.yscale("log")
plt.xscale("log")
plt.xlabel("ndof")
plt.ylabel("relative error")
for i, result in enumerate(results):
ndof, err = zip(*result)
plt.plot(ndof, err, "-*", label="t=" + str(thicknesses[i]))
plt.legend()
plt.show()


In the limit $$t \to 0$$ the shear energy term can be understood as a penalty formulation enforcing the Kirchhoff-Love assumption

\begin{align*} \beta = \nabla w. \end{align*}

Thus, in the limit case, the total energy simplifies by eliminating the rotation $$\beta$$ to

\begin{align*} \mathcal{W}_{\mathrm{Kirchhoff-Love}}(w)=\frac{1}{2}\int_{\Omega}\| \varepsilon(\nabla w) \|_{\mathbb{D}}^2-\int_{\Omega}f\,w, \end{align*}

which is the Kirchhoff-Love plate model.

Drawback of the current implementation of the TDNNS formulation for Reissner-Mindlin plates: We have a mixed saddle point problem leading to an indefinite system matrix. This can be easily overcome with hybridization techniques NGSolve docu - Static condensation by breaking the normal-normal continuity of $$\sigma$$ and reinforcing it by a Lagrange multiplier $$\alpha$$. After static condensation a symmetric positive definite system is regained.

Exercise:

• Implement hybridization for TDNNS for Reissner-Mindlin plates