The TDNNS method for Reissner-Mindlin plates

TL;DR

This paper introduces a new family of locking-free finite elements for shear deformable Reissner-Mindlin plates, based on a novel elasticity formulation that treats bending moments as separate unknowns, achieving optimal convergence without shear treatment.

Contribution

It presents a new finite element formulation for Reissner-Mindlin plates that avoids shear locking and does not require reduced integration, with optimal convergence properties.

Findings

Elements are locking-free and achieve optimal convergence.

No special shear treatment like reduced integration is needed.

The formulation effectively models shear deformable plates.

Abstract

A new family of locking-free finite elements for shear deformable Reissner-Mindlin plates is presented. The elements are based on the "tangential-displacement normal-normal-stress" formulation of elasticity. In this formulation, the bending moments are treated as separate unknowns. The degrees of freedom for the plate element are the nodal values of the deflection, tangential components of the rotations and normal-normal components of the bending strain. Contrary to other plate bending elements, no special treatment for the shear term such as reduced integration is necessary. The elements attain an optimal order of convergence.

Tables

Table 1. Table 1. Different types of essential and corresponding natural boundary conditions and the incorporation of the latter into the right hand side of ( 35 ) - ( 36 ).

essential bc	natural bc	surface term
$w$	$\mu t^{-2}\left({\partial }_{n}w-{\theta }_n\right)=g_0$	${\int }_{\Gamma }{g}_{0}v𝑑s$
${𝜽}_{𝝉}$	${𝐦}_{𝐧𝝉}={𝐠}_1$	${\int }_{\Gamma }{𝐠}_1\cdot {𝜼}_{𝝉}𝑑s$
$m_{nn}$	${\theta }_n=g_2$	${\int }_{\Gamma }{g}_2{\tau }_{nn}𝑑s$

Equations

$${\boldsymbol{\tau}}=\mathbf{n}^{\bot}=\left(\begin{array}[]{c}-n_{y}\\ n_{x}\end{array}\right).$$

$$v_n:=\mathbf{v}\cdot \mathbf{n},{\mathbf{v}}_{\mathbf{\tau }}=\mathbf{v}-v_n\mathbf{n}.$$

$${\sigma }_{nn}={\mathbf{\sigma }}_{\mathbf{n}}\cdot \mathbf{n},{\mathbf{\sigma }}_{\mathbf{n}\mathbf{\tau }}={\mathbf{\sigma }}_{\mathbf{n}}-{\sigma }_{nn}\mathbf{n}.$$

$$⟨f,g⟩=f(g)\forall f\in V^{\ast },g\in V.$$

$${\mathbf{H}}^{-1}(\mathrm{div},\Omega )=\{\mathbf{f}\in {\mathbf{H}}^{-1}(\Omega ):\mathrm{div}\mathbf{f}\in H^{-1}(\Omega )\}=({\mathbf{H}}_0(\mathrm{curl},\Omega ){)}^{\ast }.$$

$$⟨\mathrm{div}\mathbf{f},w⟩=-⟨\mathbf{f},\nabla w⟩\forall w\in C_0^{\infty }(\Omega ).$$

$$\Vert \mathbf{f}{\Vert }_{{\mathbf{H}}^{-1}(\mathrm{div},\Omega )}^2$$

$$=$$

$$\underline{\mathbf{A}}\underline{\mathbf{\sigma }}$$

$$\mathbf{div}\underline{\mathbf{\sigma }}$$

$$\mathbf{u}=0\text{mbox}on\partial \Omega .$$

$${\int }_{\Omega }\underline{\mathbf{C}}\underline{\mathbf{\epsilon }}(\mathbf{u}):\underline{\mathbf{\epsilon }}(\mathbf{v})dx={\int }_{\Omega }\mathbf{f}\cdot \mathbf{v}dx\forall \mathbf{v}\in {\mathbf{H}}_0^1.$$

$${\int }_{\Omega }\underline{\mathbf{A}}\underline{\mathbf{\sigma }}:\underline{\mathbf{\tau }}dx+{\int }_{\Omega }\mathbf{div}\underline{\mathbf{\tau }}\cdot \mathbf{u}dx$$

$${\int }_{\Omega }\mathbf{div}\underline{\mathbf{\sigma }}\cdot \mathbf{v}dx$$

$${\mathbf{H}}_0(\mathrm{curl})=\{\mathbf{v}\in \mathbf{H}(\mathrm{curl}):{\mathbf{v}}_{\mathbf{\tau }}=0\text{mbox}on\partial \Omega \}.$$

$$\underline{\mathbf{H}}(\mathrm{div}\mathbf{div})=\{\underline{\mathbf{\tau }}\in {\underline{\mathbf{L}}}_{sym}^2:\mathrm{div}\mathbf{div}\underline{\mathbf{\tau }}\in H^{-1}\}.$$

$$⟨\mathbf{div}\underline{\mathbf{\tau }},\mathbf{v}⟩=-{\int }_{\Omega }\underline{\mathbf{\tau }}:\underline{\mathbf{\epsilon }}(\mathbf{v})dx+{\int }_{\partial \Omega }{\tau }_{nn}{v}_{n}ds.$$

$$\Vert \underline{\mathbf{\tau }}{\Vert }_{\underline{\mathbf{H}}(\mathrm{div}\mathbf{div})}^2=\Vert \underline{\mathbf{\tau }}{\Vert }_{{\underline{\mathbf{L}}}^2}^2+\Vert \mathrm{div}\mathbf{div}\underline{\mathbf{\tau }}{\Vert }_{H^{-1}}^2=\Vert \underline{\mathbf{\tau }}{\Vert }_{{\underline{\mathbf{L}}}^2}^2+\underset{w\in H_0^1\cap H^2}{\sup }\frac{⟨\mathbf{div}\underline{\mathbf{\tau }},\nabla w{⟩}^2}{\Vert \nabla w{\Vert }_{{\mathbf{L}}^2}^2}.$$

$${\int }_{\Omega }\underline{\mathbf{A}}\underline{\mathbf{\sigma }}:\underline{\mathbf{\tau }}dx+⟨\mathbf{div}\underline{\mathbf{\tau }},\mathbf{u}⟩$$

$$⟨\mathbf{div}\underline{\mathbf{\sigma }},\mathbf{v}⟩$$

$$⟨\mathbf{div}\underline{\mathbf{\tau }},\mathbf{v}⟩$$

$$=$$

$$\mathbf{u}=\left(\begin{array}[]{c}-z\theta_{1}\\ -z\theta_{2}\\ w\end{array}\right),$$

$$-\mathbf{div}({\underline{\mathbf{C}}}_b\underline{\mathbf{\epsilon }}(\mathbf{\theta }))-\mu t^{-2}(\nabla w-\mathbf{\theta })$$

$$-\mu t^{-2}\mathrm{div}(\nabla w-\mathbf{\theta })$$

$$\mathbf{\theta }$$

$$w$$

$$\underline{\mathbf{C}}_{b}=\frac{E}{12(1-\nu^{2})}\left(\begin{array}[]{ccc}1&\nu&0\\ \nu&1&0\\ 0&0&\frac{1-\nu}{2}\end{array}\right),\qquad\mu=\frac{k_{s}E}{2(1+\nu)}.$$

$$\underline{\mathbf{A}}_{b}=\underline{\mathbf{C}}_{b}^{-1}=\frac{12}{E}\left(\begin{array}[]{ccc}1&-\nu&0\\ -\nu&1&0\\ 0&0&2(1+\nu)\end{array}\right).$$

$${\int }_{\Omega }{\underline{\mathbf{C}}}_b\underline{\mathbf{\epsilon }}(\mathbf{\theta }):\underline{\mathbf{\epsilon }}(\mathbf{\eta })dx+\mu t^{-2}{\int }_{\Omega }(\nabla w-\mathbf{\theta })\cdot (\nabla v-\mathbf{\eta })dx$$

$$\nabla w-\mathbf{\theta }=0$$

$${\underline{\mathbf{A}}}_b\underline{\mathbf{m}}-\underline{\mathbf{\epsilon }}(\mathbf{\theta })$$

$$-\mathbf{div}(\underline{\mathbf{m}})-\mu t^{-2}(\nabla w-\mathbf{\theta })$$

$$-\mu t^{-2}\mathrm{div}(\nabla w-\mathbf{\theta })$$

$${\int }_{\Omega }{\underline{\mathbf{A}}}_b\underline{\mathbf{m}}:\underline{\mathbf{\tau }}dx+⟨\mathbf{div}\underline{\mathbf{\tau }},\mathbf{\theta }⟩$$

Full text

1. Introduction

In this paper we are concerned with finite elements for shear deformable plates based on the Reissner-Mindlin model [42, 34]. A direct discretization of the equations leads to shear locking phenomena as the plate thickness becomes small. In the limit of zero thickness, the Kirchhoff assumption is enforced, where the shear strain vanishes and the deflection gradient equals the rotations. Over the last decades, a vast amount of different elements overcoming shear locking by different kinds of remedies has been proposed. In most standard, conforming finite element methods, the Kirchhoff constraint of vanishing shear stress is alleviated or modified in some way. An alternative are discontinuous Galerkin (DG) methods, mixed, or hybrid methods.

An example for the alleviation of the Kirchhoff constraint is the “assumed shear strain method” introduced by MacNeal [33]. A special operator for the displacement-strain relation relying on local averaging is used, this approach was further developed by [31, 10, 21]. In a further method referred to as “linked interpolation”, the displacement gradient in the Kirchhoff constraint is augmented by a “kinematic linking operator”. Pioneers in this field were Zienkiewicz and co-authors [49] and Taylor and Auricchio [47]. In [9], Auricchio and Lovadina provide an analysis of general linked-interpolation elements.

Other methods employ additional unknowns for the shear stress quantity, which allows to pose the Kirchhoff condition of vanishing shear in weak sense. Examples are the Falk-Tu element [25] or the MITC element [16]. In implementations, the further unknown can be eliminated element-wise, which leads to a projection of the shear stresses in the penalty term. This projection is referred to as “reduction” in [16]. In some cases, it may be achieved by reduced integration of the shear term. This approach is analyzed minutely for the one-dimensional case of a thick beam in [2]. Reduced integration is also used without the background of an additional shear stress unknown [48, 26], where care has to be taken to avoid spurious modes.

For MITC elements, error analysis has been provided; see e.g. the works by Brezzi, Fortin and Stenberg [17], where additionally a postprocessing step for the deflections is proposed, or by Stenberg and Suri [46] for an $hp$ error analysis.

Nonconforming elements have been constructed by Arnold and Falk [6], where Crouzeix-Raviart elements are used for the deflection. More recently, Brezzi and Marini [18] developed a nonconforming element in the framework of discontinuous Galerkin methods. In both works, the shear strain is projected into a lower-order finite element space to alleviate the Kirchhoff constraint. Other DG approaches allow for a direct enforcement of the Kirchhoff constraint, as the rotation space can be chosen such that it contains the deflection gradient. Deflections, rotations and the shear are approximated using different continuity assumptions in [5]. In [14, 27] DG methods for deflection and rotation without reduced integration techniques are presented. We also mention [19] for a discontinuous Petrov-Galerkin method, where optimal test functions of higher polynomial order are chosen to suit the trial functions. A quadrilateral hybrid finite element method was introduced in [20], where shear stress and bending moment are discontinuous and the corresponding finite element basis is constructed to satisfy a local equilibrium condition.

Hughes and Franca [30] added an additional stabilization term to the variational equations. Also, Chapelle and Stenberg [22] augmented the equations by a stabilization term, which then allows for an analysis ensuring an optimal order of convergence. An entirely different approach by Pontaza and Reddy [41] is to use a least squares method instead of the standard Galerkin equations.

Mixed method with weak symmetry for the tensor of bending moments are proposed in [12, 11]. In both works, the bending moments are approximated in the normal-continuous space ${\underline{\mathbf{H}}}(\mathbf{{div}})$. Since symmetric ${\underline{\mathbf{H}}}(\mathbf{{div}})$-conforming elements are hard to construct and of high polynomial order [8], imposing symmetry weakly using a Lagrange multiplier has often been proposed in the context of continuum mechanics, see e.g. [44, 45, 4, 7]. The continuum mechanics formulation used in the current paper overcomes this problem.

The plate elements proposed in the current paper are based on the “tangential-displacement normal-normal-stress” (TDNNS) formulation of elasticity introduced by the authors in [38]. This leads to a formulation containing deflection, rotations and bending moments as separate unknowns. While the deflection is sought in the standard Sobolev space $H^{1}$, the rotation is assumed to be in the less regular space ${\mathbf{H}(\operatorname{curl})}$. Additionally the TDNNS stress space ${\underline{\mathbf{H}}(\operatorname{div}\mathbf{{div}})}$ is chosen for the bending moments. Accordingly, we use standard continuous finite elements for the deflection, tangential continuous Nédélec elements for the rotations, and normal-normal continuous tensor valued elements from [38] for the bending moments. The main benefit of this choice is that the gradient of the deflection space is a subset of the rotation space both for the infinite dimensional and for the finite element problem, and the Kirchhoff constraint of vanishing shear strain does not lead to locking. Thus, the proposed formulation seems to be a very natural alternative to the established ones based on $H^{1}$ continuity and continuous finite elements. In [39] the authors have shown that three-dimensional anisotropic TDNNS elements are suitable for the discretization of slim domains. The proposed Reissner-Mindlin plate elements show an optimal order of convergence, which is confirmed in our numerical results.

The proposed elements are closely related to the Hellan-Herrmann-Johnson element for the bending problem of a Kirchhoff plate [28, 29, 32]. Also in the HHJ formulation, the normal-normal component of the bending moment is continuous across interfaces. As the Hellan-Herrmann-Johnson element is restricted to the biharmonic problem, the rotations are not treated independently as done in the current work. On the other hand, for vanishing thickness one can eliminate the rotations in the current formulation, and arrives at the HHJ plate formulation. Thus, the HHJ formulation may be seen as the limiting case of the proposed method. In [24], Lagrangian multipliers for the normal-normal component of the bending moment are introduced. Postprocessing then leads to a faster convergence of the deflection gradient.

This work is organized as follows: in Section 2, the TDNNS method is shortly introduced and applied to the Reissner-Mindlin problem. An analysis of the infinite dimensional problem using the TDNNS spaces for positive as well as vanishing thickness is provided in Section 3. Finite elements are introduced in Section 4, and a-priori error estimates are provided. We mention hybridization of the bending moments by Lagrangian multipliers resembling the normal component of the rotation, which results after static condensation in a symmetric positive system matrix. Finally, Section 5 contains numerical examples verifying the claimed convergence orders.

2. Problem formulation

2.1. Notation

In the following, all vectors are denoted by boldface letters, tensors are boldface and underlined. Let $\Omega\subset\mathbb{R}^{2}$ be a bounded, connected, polygonal Lipschitz domain. Its unit outward normal $\mathbf{n}$ is defined almost everywhere on the boundary $\partial\Omega$. The unit tangential vector in counter-clockwise direction ${\boldsymbol{\tau}}$ is given as the rotation of $\mathbf{n}$,

[TABLE]

For a vector field $\mathbf{v}$ on $\Omega$, normal and tangential component on the boundary are denoted by

[TABLE]

For a tensor field of second order $\underline{\boldsymbol{\sigma}}$, its normal component is given by $\boldsymbol{\sigma}_{\mathbf{n}}:=\underline{\boldsymbol{\sigma}}\mathbf{n}$. The normal component $\boldsymbol{\sigma}_{\mathbf{n}}$ can further be split into a (scalar-valued) normal-normal component $\sigma_{nn}$ and a (vector-valued) normal-tangential component $\boldsymbol{\sigma}_{\mathbf{n}{\boldsymbol{\tau}}}$,

[TABLE]

Rotation and divergence of a two-dimensional vector field shall be denoted by $\operatorname{curl}$, $\operatorname{div}$, respectively. The operator $\mathbf{{div}}$ is the row-wise divergence operator, mapping tensor to vector fields.

For a general Hilbert space $V$, its inner product and norm are denoted by $(\cdot,\cdot)_{V}$ and $\|\cdot\|_{V}$, respectively. The duality product between $V$ and its dual $V^{*}$ is denoted by angles,

[TABLE]

Let $L^{2}(\Omega)$ denote the usual Lebesgue space. Moreover, let $H^{1}(\Omega)$ be the usual Sobolev space of weakly differentiable functions, and let $H^{1}_{0}(\Omega)$ be the space of $H^{1}$ functions satisfying zero boundary conditions. We also use the Sobolev spaces $H^{s}(\Omega)$ for integer $s$. Moreover, ${\mathbf{H}(\operatorname{curl},\Omega)}$ and ${\mathbf{H}_{0}(\operatorname{curl},\Omega)}$ shall be the spaces of vector-valued functions with weak rotation, the latter satisfying zero boundary conditions for the tangential component of the vector fields, see [35].

The dual space of $H^{1}_{0}(\Omega)$ shall be denoted by $H^{-1}(\Omega)$. It is well established [13, Equation (10.4.52)] that the dual space of ${\mathbf{H}_{0}(\operatorname{curl},\Omega)}$ is ${\mathbf{H}^{-1}(\operatorname{div},\Omega)}$ being the space of $\mathbf{H}^{-1}$ vector fields with distributional divergence in $H^{-1}$,

[TABLE]

The distributional divergence operator is defined by the relationship

[TABLE]

Using this definition, a natural norm of ${\mathbf{H}^{-1}(\operatorname{div},\Omega)}$ is

[TABLE]

In all further occurrences of Sobolev spaces on $\Omega$, the domain can also be omitted. We thus write $H^{1}$ for $H^{1}(\Omega)$ or ${\mathbf{H}(\operatorname{curl})}$ for ${\mathbf{H}(\operatorname{curl},\Omega)}$. When defined on domains other than $\Omega$, the domain must be indicated as above.

2.2. The TDNNS method

In [38, 39] we introduced the tangential-displacement normal-normal-stress (TDNNS) method for elasticity, and refined the analysis in [40]. In the current section, we briefly cover the main idea of the TDNNS method for elasticity problems in the two-dimensional continuum, as the proposed Reissner-Mindlin elements will be based on this method.

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded, connected, polygonal domain with Lipschitz boundary $\partial\Omega$. The displacement vector $\mathbf{u}=(u_{1},u_{2})$ and symmetric stress tensor $\underline{\boldsymbol{\sigma}}\in\mathbb{R}^{2\times 2}_{sym}$ are connected via Hooke’s law (9) and the equilibrium equation (10)

[TABLE]

Here, we use the linearized strain tensor $\underline{\boldsymbol{\varepsilon}}(\mathbf{u})=\frac{1}{2}(\nabla\mathbf{u}+(\nabla\mathbf{u})^{T})$. Stress and strain are connected via the compliance tensor $\underline{\mathbf{A}}$, which is the inverse of the elasticity tensor $\underline{\mathbf{C}}$ depending on Young’s modulus $E$ and the Poisson ratio $\nu$ in the well known way. For simplicity, we assume that homogeneous displacement boundary conditions are posed on $\partial\Omega$,

[TABLE]

Most standard methods for the elasticity problem rely on a primal formulation, which is obtained eliminating the stress tensor $\underline{\boldsymbol{\sigma}}$ from equations (9) and (10). Then one searches for the displacement $\mathbf{u}$ in $\mathbf{H}^{1}_{0}$ such that

[TABLE]

In a conforming finite element method, the displacement vector $\mathbf{u}$ is approximated by a continuous finite element function.

On the other hand, the dual Hellinger-Reissner formulation is obtained directly from system (9), (10) when multiplying with test functions and using integration by parts,

[TABLE]

Here, the solution spaces are $\underline{\mathbf{H}}_{sym}(\mathbf{{div}})$ for the stress and $\mathbf{L}^{2}$ for the displacement. In a conforming finite element method, the displacement elements can be totally discontinuous, while the stress elements need to be tensor-valued symmetric and normal-continuous. These requests lead to computationally expensive finite elements of at least 24 degrees of freedom in two dimensions (see [8]) and 162 degrees of freedom in three dimensions (see [1, 3]).

The TDNNS formulation is “in between” those two concepts, where the displacements are not assumed totally continuous or discontinuous, but where the tangential component is assumed to be continuous across element borders. In a mathematical setting, the displacement space is chosen as

[TABLE]

In [38, 40] we have shown that the corresponding stress space is the space of symmetric ${\underline{\mathbf{L}}}^{2}$ tensors with weak divergence in the dual space of ${\mathbf{H}_{0}(\operatorname{curl})}$. Due to (5), the stress space is given by

[TABLE]

The duality product $\langle\mathbf{{div}}\underline{\boldsymbol{\tau}},\mathbf{v}\rangle$, where the divergence of a stress tensor is applied to a displacement field, plays an important role in the TDNNS method. In [40] we elaborated on the meaning of the duality product $\langle\mathbf{{div}}\underline{\boldsymbol{\tau}},\mathbf{v}\rangle$ for $\underline{\boldsymbol{\tau}}\in{\underline{\mathbf{H}}(\operatorname{div}\mathbf{{div}})}$ and $\mathbf{v}\in{\mathbf{H}_{0}(\operatorname{curl})}$. We state that, for smooth $\underline{\boldsymbol{\tau}}\in{\underline{\mathbf{H}}(\operatorname{div}\mathbf{{div}})}$ and $\mathbf{v}\in{\mathbf{H}_{0}(\operatorname{curl})}$, i.e. for $\mathbf{v}$ smooth with vanishing tangential component $\mathbf{v}_{\boldsymbol{\tau}}=0$ on $\partial\Omega$, the duality product can be evaluated by

[TABLE]

A natural norm of the stress space ${\underline{\mathbf{H}}(\operatorname{div}\mathbf{{div}})}$ uses this duality product and is given by, see [40]

[TABLE]

It is well known that finite elements for ${\mathbf{H}(\operatorname{curl})}$ have to be tangential continuous, such as Nédélec elements introduced in [36, 37]. In [38] it was shown that finite elements for the stress space ${\underline{\mathbf{H}}(\operatorname{div}\mathbf{{div}})}$ are normal-normal continuous, meaning that the normal component $\sigma_{nn}$ of the normal stress vector $\boldsymbol{\sigma}_{\mathbf{n}}$ is continuous across element borders.

One obtains a variational problem of a form similar to the dual problem (13) - (14),

[TABLE]

In a finite element method, it is necessary to evaluate duality products of the form $\langle\mathbf{{div}}\underline{\boldsymbol{\tau}},\mathbf{v}\rangle$ for piecewise smooth functions with $\tau_{nn}$ and $\mathbf{v}_{\boldsymbol{\tau}}$ continuous on a finite element mesh $\mathcal{T}=\{T\}$. In this case, the definition from (17) can be extended to

[TABLE]

It was shown (see [38, 40]) that the infinite dimensional problem (19), (20) is well posed. Moreover, a stable family of mixed finite elements was constructed, using Nédélec’s elements for the displacement space and a new family of tensor-valued symmetric normal-normal continuous elements for the stress space. The two-dimensional mixed finite elements shall be used in the Reissner-Mindlin elements proposed in this work.

2.3. Reissner-Mindlin model

Let again $\Omega\subset\mathbb{R}^{2}$ be a bounded, connected domain with Lipschitz boundary $\partial\Omega$. We consider a plate of thickness $t$ corresponding to the three-dimensional domain $\Omega\times(-t/2,t/2)$. In the Reissner-Mindlin model, the displacement vector $\mathbf{u}$ is assumed to take the form

[TABLE]

where $\boldsymbol{\theta}=(\theta_{1},\theta_{2})$ are rotations and $w$ is the deflection in vertical $z$ direction. Both the rotations $\boldsymbol{\theta}$ and the deflection $w$ are assumed to depend on the in-plane coordinates $(x_{1},x_{2})$ only.

Assuming a vertical volume load $\mathbf{f}=(0,0,t^{2}g)^{T}\in{\mathbf{L}}^{2}$ to be given, the Reissner-Mindlin problem for a clamped plate is to find the deflection $w$ and rotations $\boldsymbol{\theta}$ such that

[TABLE]

Here $\underline{\mathbf{C}}_{b}$ is the tensor of bending moduli and $\mu$ is the shear modulus with shear correction factor $k_{s}$, which depend on Young’s modulus $E$ and Poisson’s ratio $\nu$ via

[TABLE]

Additionally, we provide the compliance tensor $\underline{\mathbf{A}}_{b}=\underline{\mathbf{C}}_{b}^{-1}$

[TABLE]

Of course, other boundary conditions such as simply supported or free boundaries, or boundary tractions and moments, may be prescribed. Although the analysis of the proposed finite element formulation is done for the clamped case for sake of simplicity, we shall comment shortly on the implementation of other boundary conditions in the end of the current section. We will see that all common types of boundary conditions can be treated in a very natural way.

As in the continuum problem in Section 2.2, a primal, displacement-based variational formulation of the Reissner-Mindlin problem (24)-(27) can be obtained directly. Both the rotations $\boldsymbol{\theta}$ and the deflection $w$ are assumed weakly differentiable, with $\boldsymbol{\theta}\in\mathbf{\Theta}_{prim}=\mathbf{H}^{1}_{0}$ and $w\in W_{prim}=H^{1}_{0}$. The primal variational formulation of the Reissner-Mindlin problem is to find $\boldsymbol{\theta}\in\mathbf{\Theta}_{prim}$ and $w\in W_{prim}$ such that for all test functions $\boldsymbol{\eta}\in\mathbf{\Theta}_{prim}$ and $v\in W_{prim}$

[TABLE]

A straightforward finite element discretization of this primal problem choosing continuous finite element spaces $\mathbf{\Theta}_{prim,h}\subset\mathbf{\Theta}_{prim}$ and $W_{prim,h}\subset W_{prim}$ leads to shear locking as the thickness $t$ tends to zero. In the limit case of a Kirchhoff plate with $t=0$, the condition of vanishing shear strain

[TABLE]

has to be satisfied. For conventional finite element discretizations one observes that $\nabla W_{prim,h}\not\subset\mathbf{\Theta}_{prim,h}$, thus the Kirchhoff constraint (31) cannot be satisfied by the discrete solution, the formulation locks. Different methods have been proposed to reduce this phenomenon by alleviating the Kirchhoff constraint (31), see Section 1. In this work, the rotation space will be chosen such that both $\nabla W\subset\mathbf{\Theta}$ and $\nabla W_{h}\subset\mathbf{\Theta}_{h}$. This ensures stability and an optimal order of convergence of the method.

The main idea of the current work is to use the TDNNS method presented in Section 2.2 for the discretization of rotations $\boldsymbol{\theta}$. To this end, additional unknowns $\underline{\mathbf{m}}=\underline{\mathbf{C}}_{b}\underline{\boldsymbol{\varepsilon}}(\boldsymbol{\theta})$ for the tensor of bending moments are introduced. This leads to the following system of partial differential equations

[TABLE]

Now, we obtain a variational formulation for finding $\boldsymbol{\theta}\in\mathbf{\Theta}={\mathbf{H}_{0}(\operatorname{curl})}$, $\underline{\mathbf{m}}\in\underline{\mathbf{M}}={\underline{\mathbf{H}}(\operatorname{div}\mathbf{{div}})}$ and $w\in W=H^{1}_{0}$ in the same manner as in Section 2.2,

[TABLE]

We shortly comment on the changes in the variational formulation (35) - (36) which are necessary for the incorporation of different types of boundary conditions. Essential boundary conditions, which have to be enforced by the finite element space, are the deflection $w$, the tangential component of the rotation $\boldsymbol{\theta}_{\boldsymbol{\tau}}$ and the normal component of the bending moment $m_{nn}$. Note that these expressions are also degrees of freedom of the corresponding finite element spaces. The corresponding natural conditions, in the same order, are the shear $\mu t^{-2}(\partial_{n}w-\theta_{n})$, the tangential component of the bending moment $\mathbf{m}_{nt}$, and the normal component of the rotation $\theta_{n}$. Natural homogeneous conditions are satisfied whenever the corresponding essential condition is dropped, inhomogeneous conditions result in additional surface integrals on the right hand side, see Table 1.

The analysis of system (35), (36) is subject of Section 3, while a finite element method is constructed and analyzed in Section 4.

3. Analysis of the TDNNS Reissner-Mindlin formulation

In the current section we show existence and uniqueness of the solution to the Reissner-Mindlin problem in the TDNNS setting. We show stability for decreasing thickness $t\to 0^{+}$ as well as the limit case $t=0$.

To this end, a further unknown $\boldsymbol{\gamma}=-\mu t^{-2}(\nabla w-\boldsymbol{\theta})$ related to the shear stresses is introduced, see e.g. [13, Chapter 10.4]. For positive thickness $t$, we choose the corresponding space $\mathbf{\Gamma}={\mathbf{L}^{2}(\Omega)}$. Problem (35), (36) transforms to

[TABLE]

We observe that for system (37) - (39) also the limit case of $t=0$ is well-defined, where the term $\mu^{-1}t^{2}\int_{\Omega}\boldsymbol{\gamma}\cdot\boldsymbol{\delta}\,dx$ vanishes. We reorder terms to obtain a mixed problem in the spirit of [13]. We introduce bilinear forms $a_{t}:(\underline{\mathbf{M}}\times\mathbf{\Gamma})\times(\underline{\mathbf{M}}\times\mathbf{\Gamma})\to\mathbb{R}$ depending on $t$ and $b:(\underline{\mathbf{M}}\times\mathbf{\Gamma})\times(\mathbf{\Theta}\times W)\to\mathbb{R}$ independent of $t$

[TABLE]

From (37), (38), (39) we obtain a saddle point problem of finding $\underline{\mathbf{m}}\in\underline{\mathbf{M}}$, $\boldsymbol{\gamma}\in\mathbf{\Gamma}$, $\boldsymbol{\theta}\in\mathbf{\Theta}$ and $w\in W$ such that

[TABLE]

In the limiting case of an infinitely thin (Kirchhoff) plate with $t=0$, it is well known (see e.g. [13, Proposition 10.4.3]) that for $t\to 0$ the shear $\boldsymbol{\gamma}$ stays bounded in $\mathbf{\Gamma}_{0}:={\mathbf{H}^{-1}(\operatorname{div})}$. We will see that $\mathbf{\Gamma}_{0}$ is the natural space for the analysis of the case $t=0$. Note that $a_{0}(\cdot,\cdot)$ is well-defined on $\underline{\mathbf{M}}\times\mathbf{\Gamma}_{0}$, while for $t>0$ $a_{t}(\cdot,\cdot)$, cannot be evaluated on the whole space $\underline{\mathbf{M}}\times\mathbf{\Gamma}_{0}$. The bilinear form $b(\cdot,\cdot)$ is also well-defined in the limiting case, as we shall show below.

For the stability analysis of (42), (43) by Brezzi’s theory [13, Theorem 4.2.3] a characterization of the kernel space $\operatorname{Ker}(B)$ is needed, which is provided in the following lemma.

Lemma 1**.**

Define the kernel space

[TABLE]

Then any $(\underline{\mathbf{m}},\boldsymbol{\gamma})\in\operatorname{Ker}(B)$ satisfies

[TABLE]

and

[TABLE]

These equalities also hold when $\mathbf{\Gamma}={\mathbf{L}}^{2}$ is replaced by $\mathbf{\Gamma}_{0}={\mathbf{H}}^{-1}(\operatorname{div})$.

Proof.

The proof follows directly, setting either $\boldsymbol{\eta}=\nabla v$ or $v=0$ in (44). Note that all duality products are well-defined due to the choice of spaces. ∎

3.1. The limiting case $t=0$

We prove existence and uniqueness of the solution for the limiting case of an infinitely thin (Kirchhoff) plate with $t=0$. As mentioned in the introduction, this case is closely related to the Hellan-Herrmann-Johnson formulation, when setting $\boldsymbol{\theta}=\nabla w$ and eliminating thereby the unknowns $\boldsymbol{\theta}$ and $\boldsymbol{\gamma}$. However, we will present an analysis of the full mixed system, as it will help understand the case of small thickness $t>0$.

As already mentioned, in this case the natural choice for the shear space is $\mathbf{\Gamma}_{0}:={\mathbf{H}^{-1}(\operatorname{div})}$. We show boundedness and stability estimates for the bilinear forms, which implies existence, uniqueness and stability of the solution [13, Theorem 4.2.3]. We use the following natural norms

[TABLE]

Lemma 2**.**

The bilinear form $a_{0}:(\underline{\mathbf{M}}\times\mathbf{\Gamma}_{0})\times(\underline{\mathbf{M}}\times\mathbf{\Gamma}_{0})$ is bounded, for all $\underline{\mathbf{m}},\underline{\boldsymbol{\tau}}\in\underline{\mathbf{M}}$ and $\boldsymbol{\gamma},\boldsymbol{\delta}\in\mathbf{\Gamma}_{0}$

[TABLE]

Moreover, it is coercive on $\operatorname{Ker}(B)$, for all $(\underline{\mathbf{m}},\boldsymbol{\gamma})\in\operatorname{Ker}(B)$

[TABLE]

Proof.

Boundedness of $a_{0}$ is straightforward, since $\underline{\mathbf{M}}\subset{\underline{\mathbf{L}}}^{2}$, and

[TABLE]

We proceed to showing coercivity of $a_{0}$ on $\operatorname{Ker}(B)$. We use equations (45) and (46) of Lemma 1 to bound $\boldsymbol{\gamma}$ by $\underline{\mathbf{m}}$ in their respective norms. Note that for $\mathbf{v}\in{\mathbf{H}}^{1}_{0}$ and $w\in H^{1}_{0}$ we have $\mathbf{v}\in{\mathbf{H}}_{0}(\operatorname{curl})$ and $\nabla w\in{\mathbf{H}}_{0}(\operatorname{curl})$.

[TABLE]

Now coercivity of $a_{0}$ is ensured by the minimal eigenvalue $\lambda_{min}(\underline{\mathbf{A}}_{b})$ of the compliance tensor $\underline{\mathbf{A}}_{b}$,

[TABLE]

∎

Lemma 3**.**

The bilinear form $b:(\underline{\mathbf{M}}\times\mathbf{\Gamma}_{0})\times(\mathbf{\Theta}\times W)$ is bounded, for all $\underline{\boldsymbol{\tau}}\in\underline{\mathbf{M}}$, $\boldsymbol{\delta}\in\mathbf{\Gamma}_{0}$, $\boldsymbol{\eta}\in\mathbf{\Theta}$ and $v\in W$

[TABLE]

Moreover, it is satisfies an inf-sup condition, for all $\boldsymbol{\theta}\in\mathbf{\Theta}$, $w\in W$ there exist $\underline{\mathbf{m}}\in\underline{\mathbf{M}}$, $\boldsymbol{\gamma}\in\mathbf{\Gamma}$ such that

[TABLE]

Proof.

Boundedness follows directly by the choice of spaces. We proceed to show the inf-sup condition. Let $\boldsymbol{\theta}\in\mathbf{\Theta},w\in W$ be fixed. From [40] we know that $\langle\mathbf{{div}}\underline{\mathbf{m}},\boldsymbol{\theta}\rangle$ is inf-sup stable on ${\underline{\mathbf{H}}(\operatorname{div}\mathbf{{div}})}\times{\mathbf{H}_{0}(\operatorname{curl})}$, i.e. there exists some $\tilde{}\underline{\mathbf{m}}\in{\underline{\mathbf{H}}(\operatorname{div}\mathbf{{div}})}$ such that

[TABLE]

and we have shown

[TABLE]

with $c_{1}>0$ and $0<\underline{c}\leq 1\leq\overline{c}$. Moreover, since ${\mathbf{H}^{-1}(\operatorname{div})}=\mathbf{\Gamma}_{0}$ is the dual space of ${\mathbf{H}_{0}(\operatorname{curl})}=\mathbf{\Theta}$, and since $\nabla W\subset\mathbf{\Theta}$, by the Riesz Isomorphism there exists some $\tilde{\boldsymbol{\gamma}}\in\mathbf{\Gamma}_{0}$ such that

[TABLE]

and

[TABLE]

In the remainder of the proof we verify that the pair $(\underline{\mathbf{m}},\boldsymbol{\gamma}):=(\frac{2}{c_{1}\underline{c}}\tilde{}\underline{\mathbf{m}},\tilde{\boldsymbol{\gamma}})$ satisfies the inf-sup condition (60) with stability constant $\underline{\beta}_{0}=\frac{c_{1}\underline{c}}{\sqrt{2}\overline{c}}$. First, we observe

[TABLE]

Using the triangle inequality in the form that $\|\boldsymbol{\theta}\|_{{\mathbf{H}(\operatorname{curl})}}^{2}+\|\nabla w-\boldsymbol{\theta}\|_{{\mathbf{H}(\operatorname{curl})}}^{2}\geq\frac{1}{2}\|\nabla w\|_{{\mathbf{H}(\operatorname{curl})}}^{2}$, and the fact that $\|\nabla w\|_{{\mathbf{H}(\operatorname{curl})}}=\|\nabla w\|_{{\mathbf{L}}^{2}}$, we obtain

[TABLE]

Since $2\|\boldsymbol{\theta}\|_{{\mathbf{H}(\operatorname{curl})}}\geq\frac{2}{\overline{c}}\|\tilde{}\underline{\mathbf{m}}\|_{{\underline{\mathbf{H}}(\operatorname{div}\mathbf{{div}})}}=\frac{c_{1}\underline{c}}{\overline{c}}\|\underline{\mathbf{m}}\|_{{\underline{\mathbf{H}}(\operatorname{div}\mathbf{{div}})}}$ due to equation (62) and the definition of $\underline{\mathbf{m}}$, and $\|\boldsymbol{\gamma}\|_{H^{-1}(\operatorname{div})}=\|\nabla w-\boldsymbol{\theta}\|_{H(\operatorname{curl})}$ due to equation (63), we conclude

[TABLE]

∎

Theorem 4**.**

For $t=0$, problem (42), (43) has a unique solution $\underline{\mathbf{m}}\in\underline{\mathbf{M}}$, $w\in W$, $\boldsymbol{\theta}\in\mathbf{\Theta}$ and $\boldsymbol{\gamma}\in\mathbf{\Gamma}_{0}={\mathbf{H}^{-1}(\operatorname{div})}$. The solution is bounded by

[TABLE]

with $c$ a generic constant.

3.2. The case of positive thickness $t>0$

In this section, we prove existence and uniqueness of a solution to the Reissner-Mindlin problem (37) - (39) in the case of positive thickness $t>0$. To this end, a different set of norms is introduced, that includes the thickness $t$.

[TABLE]

Lemma 5**.**

A norm equivalent to $\|\boldsymbol{\theta},w\|_{\mathbf{\Theta}\times W}$ can be defined omitting the term $\|w\|_{H^{1}}$ in (75), where the non-trivial bound is characterized by the Friedrichs constant $c_{F}$,

[TABLE]

Proof.

The statement of the lemma is clear from the following consideration, which uses Friedrichs’ inequality, the triangle inequality and the fact that $t<1$.

[TABLE]

The next two lemmas provide stability estimates for the bilinear forms in the $t$-dependent norms $\|\underline{\mathbf{m}},\boldsymbol{\gamma}\|_{\underline{\mathbf{M}}\times\mathbf{\Gamma},t}$ and $\|\boldsymbol{\theta},w\|_{\mathbf{\Theta}\times W,t}$. These estimates are used to ensure existence and uniqueness of a solution and to obtain a stability estimate not deteriorating with $t\to 0$. The proof of Lemma 6 is straightforward in the manner of the proof of Lemma 2:

Lemma 6**.**

For $t>0$, the bilinear form $a_{t}:(\underline{\mathbf{M}}\times\mathbf{\Gamma})\times(\underline{\mathbf{M}}\times\mathbf{\Gamma})$ is bounded, for all $\underline{\mathbf{m}},\underline{\boldsymbol{\tau}}\in\underline{\mathbf{M}}$ and $\boldsymbol{\gamma},\boldsymbol{\delta}\in\mathbf{\Gamma}$

[TABLE]

It is coercive on $\operatorname{Ker}(B)$, for all $(\underline{\mathbf{m}},\boldsymbol{\gamma})\in\operatorname{Ker}(B)$

[TABLE]

The constants $\bar{\alpha}$, $\underline{\alpha}$ are independent of the thickness $t$.

Lemma 7**.**

The bilinear form $b:(\underline{\mathbf{M}}\times\mathbf{\Gamma})\times(\mathbf{\Theta}\times W)$ is continuous with respect to the $t$-dependent norms, for all $\underline{\boldsymbol{\tau}}\in\underline{\mathbf{M}}$, $\boldsymbol{\delta}\in\mathbf{\Gamma}$, $\boldsymbol{\eta}\in\mathbf{\Theta}$ and $v\in W$

[TABLE]

Moreover, it is satisfies an inf-sup condition, for all $\boldsymbol{\theta}\in\mathbf{\Theta}$, $w\in W$ there exist $\underline{\mathbf{m}}\in\underline{\mathbf{M}}$, $\boldsymbol{\gamma}\in\mathbf{\Gamma}$ such that

[TABLE]

The constants $\bar{\beta}$, $\underline{\beta}$ are independent of the thickness $t$.

Proof.

Obviously, the bilinear form is bounded, as the divergence term is bounded in ${\underline{\mathbf{H}}}(\operatorname{div}\mathbf{{div}})\times{\mathbf{H}}(\operatorname{curl})$, and the integral is bounded by the respective scaled ${\mathbf{L}}^{2}$-norms $t\|\boldsymbol{\delta}\|_{{\mathbf{L}}^{2}}$ and $t^{-1}\|\nabla w-\theta\|_{{\mathbf{L}}^{2}}$.

To prove the inf-sup condition, assume $\boldsymbol{\theta}\in\mathbf{\Theta}$ and $w\in W$ are given. Similar to the proof of Lemma 3, we use the theory provided in [40]. We choose $\underline{\mathbf{m}}=\tilde{}\underline{\mathbf{m}}\in\underline{\mathbf{M}}={\underline{\mathbf{H}}(\operatorname{div}\mathbf{{div}})}$ from equations (61) and (62). Additionally, we choose $\boldsymbol{\gamma}=t^{-2}(\nabla w-\boldsymbol{\theta})$. This is possible since $\nabla w-\boldsymbol{\theta}\in{\mathbf{H}(\operatorname{curl})}\subset\mathbf{L}^{2}$. Then we have

[TABLE]

Combining the definition of $b(\cdot,\cdot)$ (41), equation (83) and the bounds from (61) and (62) we obtain

[TABLE]

Basic algebra for real numbers, the definition of the shear $\boldsymbol{\gamma}$ and the upper bound in (62) lead to

[TABLE]

We see inf-sup stability for $b(\cdot,\cdot)$ using the bound from Lemma 5

[TABLE]

∎

Theorem 8**.**

For $t>0$, problem (42), (43) has a unique solution $\underline{\mathbf{m}}\in\underline{\mathbf{M}}$, $w\in W$, $\boldsymbol{\theta}\in\mathbf{\Theta}$ and $\boldsymbol{\gamma}\in\mathbf{\Gamma}$, which is bounded as below

[TABLE]

where $c$ is a generic constant independent of $t$.

Proof.

Again, we use the statement of [13, Theorem 4.2.3], where coercivity (Lemma 6) and inf-sup stability (Lemma 7) ensure the existence and stability of a unique solution. Note that, for the solution $\underline{\mathbf{m}}\in\underline{\mathbf{M}}$, $w\in W$, $\boldsymbol{\theta}\in\mathbf{\Theta}$ and $\boldsymbol{\gamma}\in\mathbf{\Gamma}$, there holds $\boldsymbol{\gamma}=\mu t^{-2}(\nabla w-\boldsymbol{\theta})$. Thus $t\|\boldsymbol{\gamma}\|_{{\mathbf{L}}^{2}}=t^{-1}\|\nabla w-\boldsymbol{\theta}\|_{{\mathbf{L}}^{2}}$, which ensures the estimate

[TABLE]

We add the bound on $\|\boldsymbol{\gamma}\|_{{\mathbf{H}}^{-1}(\operatorname{div})}$: since ${\mathbf{H}}^{-1}(\operatorname{div})$ is dual to $\mathbf{\Theta}={\mathbf{H}}(\operatorname{curl})$ and the solution $\boldsymbol{\gamma}$ satisfies the variational equation (38) with $v=0$ we have

[TABLE]

Thus, the statement of the theorem is shown.

∎

4. Finite Elements

Throughout this section, let $\Omega\subset\mathbb{R}^{2}$ be a polygonal Lipschitz domain, and let $(\mathcal{T}_{h})$ be a family of decompositions into triangular elements $T$. We assume that the family $(\mathcal{T}_{h})$ is regular, shape-regular and quasi-uniform with mesh size $h$ (see e.g. [15]). Moreover, $(\mathcal{F}_{h})$ shall denote the set of edges in the mesh.

We propose a finite element method without using the additional unknown shear $\boldsymbol{\gamma}$, having the form of (35)-(36). We use the TDNNS finite element spaces for the bending moments and rotations, and a fully continuous Lagrange space for the deflection. For integer $k\geq 1$, $\underline{\mathbf{m}}$ is approximated in the normal-normal continuous space of order $k$ introduced in [38], while $\boldsymbol{\theta}$ is discretized by order $k$ Nédélec elements of the second kind [37]. The deflection elements are order $k+1$ continuous elements. In detail, we choose

[TABLE]

The discrete system used in implementations reads

[TABLE]

Note that the finite element space $\underline{\mathbf{M}}_{h}$ is (slightly) non-conforming, $\underline{\mathbf{M}}_{h}\not\subset\underline{\mathbf{M}}={\underline{\mathbf{H}}}(\operatorname{div}\mathbf{{div}})$. This is due to lacking continuity of $\mathbf{m}_{\mathbf{n}{\boldsymbol{\tau}}}$ at the corner points (in the interior) of each element, see [40, page 13] for a detailed discussion. However, the duality product $\langle\mathbf{{div}}\underline{\boldsymbol{\tau}}_{h},\boldsymbol{\theta}_{h}\rangle$ can now be understood as the $\mathbf{\Theta}_{h}^{*}\times\mathbf{\Theta}_{h}$ duality product, and can be evaluated by the relations (21) - (22). Moreover, the norm $\|\cdot\|_{\underline{\mathbf{M}}}$ from (18) is not well-defined for $\underline{\mathbf{m}}_{h}\in\underline{\mathbf{M}}_{h}$. In [40] we provided a discrete norm and a corresponding stability analysis for the TDNNS continuum mechanics elements. We will use this discrete norm in the current paper, defining

[TABLE]

and the parameter-dependent norm

[TABLE]

For finite element tensors $\underline{\mathbf{m}}_{h}\in\underline{\mathbf{M}}_{h}$, the edge $L^{2}$ terms in the norm above can also be omitted, as they are bounded by the domain ${\underline{\mathbf{L}}}^{2}$ norm. The divergence operator $\mathbf{{div}}:\underline{\mathbf{M}}_{h}\to\mathbf{\Theta}_{h}^{*}$ is bounded and LBB-stable, see [40]:

Theorem 9**.**

There exist positive constants $\beta_{1},\beta_{2}>0$ such that for any $\underline{\mathbf{m}}_{h}\in\underline{\mathbf{M}}_{h}$ and $\boldsymbol{\theta}_{h}\in\mathbf{\Theta}_{h}$

[TABLE]

and

[TABLE]

4.1. Discrete stability

For the analysis, it is convenient to introduce a finite element discretization for the shear $\boldsymbol{\gamma}$, which leads to a discrete system equivalent to (96), (97), but which is of the standard saddle point form (42), (43). The equivalence of the discrete systems is due to the inclusion $\nabla W_{h}\subset\mathbf{\Theta}_{h}$, and to our choice $\mathbf{\Gamma}_{h}=\mathbf{\Theta}_{h}$. Thus, for $w_{h}\in W_{h}$, $\boldsymbol{\theta}_{h}\in\mathbf{\Theta}_{h}$ and $\boldsymbol{\gamma}_{h}\in\mathbf{\Gamma}_{h}$ the discrete variational equation

[TABLE]

is equivalent to $\boldsymbol{\gamma}_{h}=\mu t^{-2}(\nabla w_{h}-\boldsymbol{\theta}_{h})$. This implies that $\boldsymbol{\gamma}_{h}$ can be eliminated, and the smaller system (96), (97) may be used in implementations.

The stability analysis is similar to the analysis of the infinite dimensional problem for positive thickness presented in Section 3.2.

Lemma 10**.**

The bilinear form $a_{t}:(\underline{\mathbf{M}}_{h}\times\mathbf{\Gamma}_{h})\times(\underline{\mathbf{M}}_{h}\times\mathbf{\Gamma}_{h})$ is coercive on $\operatorname{Ker}(B_{h}):=\{(\underline{\mathbf{m}}_{h},\boldsymbol{\gamma}_{h})\in\underline{\mathbf{M}}_{h}\times\mathbf{\Gamma}_{h}:b(\underline{\mathbf{m}}_{h},\boldsymbol{\gamma}_{h};\boldsymbol{\theta}_{h},w_{h})=0\ \forall\boldsymbol{\theta}_{h}\in\mathbf{\Theta}_{h},w_{h}\in W_{h}\}$. There exists a constant $\alpha_{1}>0$ independent of $t,h$ such that

[TABLE]

Proof.

For any $(\underline{\mathbf{m}}_{h},\boldsymbol{\gamma}_{h})\in\operatorname{Ker}(B_{h})$ we have by definition, setting $\boldsymbol{\theta}_{h}=\nabla w_{h}$,

[TABLE]

Thus, it follows

[TABLE]

∎

Lemma 11**.**

The bilinear form $b:(\underline{\mathbf{M}}_{h}\times\mathbf{\Gamma}_{h})\times(\mathbf{\Theta}_{h}\times W_{h})$ is bounded and inf-sup stable, for any $\boldsymbol{\theta}_{h}\in\mathbf{\Theta}_{h},w_{h}\in W_{h}$ there exist $\underline{\mathbf{m}}_{h}\in\underline{\mathbf{M}}_{h},\boldsymbol{\gamma}_{h}\in\mathbf{\Gamma}_{h}$ such that

[TABLE]

Boundedness is clear from the discrete boundedness of the divergence operator, see Theorem 9.

The discrete inf-sup condition is shown in the same manner as in the infinite dimensional case in Theorem 8. The arguments shall not be repeated here, but only shortly commented on.

We use the discrete inf-sup stability of the divergence operator from Theorem 9. As in Theorem 8, we set $\boldsymbol{\gamma}_{h}=t^{-2}(\nabla w_{h}-\boldsymbol{\theta}_{h})$, which is possible since $\nabla W_{h}\subset\mathbf{\Theta}_{h}=\mathbf{\Gamma}_{h}$. The remainder of the proof involves the same steps as shown in eq. (84)-(89), only replacing the infinite-dimensional norms by the discrete ones.

4.2. A-priori error estimates

To get a-priori error estimates, it is necessary to have interpolation error estimates. We use the standard nodal interpolation operator $\mathcal{I}_{W}$ for $H^{1}$ and the standard interpolator $\mathcal{I}_{\mathbf{\Theta}}$ of the Nédélec space defined using its degrees of freedom, see e.g. [35] for their definition. The following approximation properties for sufficiently smooth functions are provided there for $1\leq m\leq k$,

[TABLE]

An important property of the interpolation operators is that they commute with the gradient operator, see e.g. [35, Theorem 5.49],

[TABLE]

For the bending moments $\underline{\mathbf{m}}$ we use the nodal interpolation operator $\mathcal{I}_{\underline{\mathbf{M}}}$, which is provided and analyzed in [40]. An error estimate in the discrete ${\underline{\mathbf{H}}}(\operatorname{div}\mathbf{{div}})$ norm was found for $\mathcal{I}_{\underline{\mathbf{M}}}$ for $0\leq l\leq k$

[TABLE]

We shall not provide the degrees of freedom of the stress space in detail here, only note that in two space dimensions there are edge-based degrees of freedom coupling the normal-normal component $m_{nn}$ of $\underline{\mathbf{m}}$, and inner degrees of freedom which can be eliminated by static condensation. Corresponding polynomial basis functions can be found in [38].

We can now venture to show a convergence result of the proposed finite element method.

Theorem 12**.**

Let $\underline{\mathbf{m}}\in\underline{\mathbf{M}}$, $w\in W$ and $\boldsymbol{\theta}\in\mathbf{\Theta}$ be the exact solution to the Reissner-Mindlin problem (35), (36), and let $\underline{\mathbf{m}}_{h}\in\underline{\mathbf{M}}_{h}$, $w_{h}\in W_{h}$ and $\boldsymbol{\theta}_{h}\in\mathbf{\Theta}_{h}$ be the corresponding finite element solution. Then we have the a-priori error estimate for $1\leq m\leq k$

[TABLE]

Proof.

Since the finite element method is slightly nonconforming, $\underline{\mathbf{M}}_{h}\not\subset\underline{\mathbf{M}}$, see [40], we use techniques from Strang’s second lemma. We bound the total error (115) by interpolation error (118) and consistency error (121).

[TABLE]

Clearly, the interpolation error (118) can be bounded as stated above (114). We concentrate on the consistency error.

As stated in [13, Theorem 5.2.1], discrete stability ensures

[TABLE]

The first term, (130), was treated in [40], and is bounded by

[TABLE]

For the second term (130), we used the commuting diagram property of the interpolation operators $\mathcal{I}_{W}$ and $\mathcal{I}_{\mathbf{\Theta}}$ (111) and the linearity of $\mathcal{I}_{\mathbf{\Theta}}$,

[TABLE]

The discrete solution $\boldsymbol{\gamma}_{h}$ satisfies $\boldsymbol{\gamma}_{h}=\mu t^{-2}(\boldsymbol{\theta}_{h}-\nabla w_{h})$ (see (102)). As the solution $\boldsymbol{\gamma}$ satisfies $\boldsymbol{\gamma}=\mu t^{-2}(\boldsymbol{\theta}-\nabla w)$, we obtain

[TABLE]

We proceed to estimating the last term (130). Since $(\underline{\mathbf{m}}_{h},\boldsymbol{\gamma}_{h})$ and $(\underline{\mathbf{m}},\boldsymbol{\gamma})$ are solutions to the discrete and infinite-dimensional variational equations, and since $(\boldsymbol{\eta}_{h},w_{h})\in\mathbf{\Theta}_{h}\times W_{h}\subset\mathbf{\Theta}\times W$, we have

[TABLE]

The first term above and the divergence of similar differences of $\underline{\mathbf{m}}$ and discrete tensors in $\underline{\mathbf{M}}_{h}$ is well-defined in the sense

[TABLE]

We may rewrite

[TABLE]

In [40] we have shown that

[TABLE]

Thus, we deduce

[TABLE]

Consequently, we have arrived at the desired result. ∎

4.3. Hybridization

To avoid the implementation of normal-normal continuous finite elements and an indefinite system matrix, a hybridization technique in the spirit of [13, Chapter 7.1] was mentioned in [38] and analyzed in [43]. Here, the normal-normal continuity of the tensor of bending moments is broken and imposed by Lagrangian multipliers defined on element edges. The Lagrangian multipliers resemble the normal component of the rotation $\gamma_{n}$. As the Lagrangian multipliers are chosen of the same polynomial order as the normal-normal component of the bending moment, the discrete systems are equivalent. However, now the bending moment $\underline{\mathbf{m}}$ is completely local and can be eliminated element-wise (static condensation). The remaining system contains only displacement-based unknowns. It is symmetric positive definite, which makes it easier to be solved by sparse direct solver or standard iterative solvers.

5. Numerical example

5.1. Clamped square plate

The first example is taken from [23], where the solution is known analytically. We consider a clamped square plate $\Omega=(0,1)^{2}$, i.e. at the boundary deflection $w=0$ and rotation $\boldsymbol{\theta}=0$ vanish. The plate thickness varies from $t=0.1$ to $t=10^{-5}$. Young’s modulus and Poisson ratio are chosen as $E=12$ and $\nu=0$. The shear correction factor is set to $k_{s}=5/6$. The vertical component of the volume load is chosen as

[TABLE]

The solution $(\boldsymbol{\theta},w)$ is given by

[TABLE]

Two discretization methods are compared: the MITC7 element [16] and the TDNNS element for $k=1$ and $k=2$. In case of the MITC7 element and the TDNNS element with $k=1$, the deflections are discretized by polynomials of order two. For the higher-order TDNNS element, the deflections are of order three.

First, we compare the different methods for a thickness of $t=10^{-3}$. In Figure 1, the convergence of $\|w-w_{h}\|_{L^{2}(\Omega)}$ is shown, Figure 2 displays the convergence of $\|\boldsymbol{\theta}-\boldsymbol{\theta}_{h}\|_{{\mathbf{L}}^{2}(\Omega)}$. It shows that for the deflection $w$, both the MITC7 element and the lowest-order TDNNS element with $k=1$ show convergence order three, while the TDNNS element with $k=2$ converges, as expected, at order four. However, for the rotations $\theta$, the MITC7 element and the TDNNS element with $k=2$ converge at the same rate of order three, while the lowest-order TDNNS element shows a convergence rate of order two. Thus, from the point of view of convergence, the MITC7 element lies in between the TDNNS elements with $k=1$ and $k=2$.

Next, we plot the convergence of the lowest order TDNNS method for different thicknesses. Figure 3 shows the convergence of the method for thicknesses varying between $t=0.1$ and $t=10^{-5}$. The error curves are very close, as the method does not suffer from the degrading thickness.

5.2. Square plate with hole

In the second example, we consider a square plate of dimensions $100\times 100$ mm, in which a circular hole of diameter $d=30$ mm is cut. The Young’s modulus $E=2.1\times 10^{5}$ N/mm2 and Poisson ratio $\nu=0.3$ are chosen as those of steel. The shear correction factor is set to $k_{s}=5/6$. The plate is clamped at the left hand side ($x=0$), and a surface traction $\sigma_{zz}=0.1(y-50)$ N/mm2 acts on the right hand side. All other boundaries are free. See Figure 4 for a sketch of the setup.

An initial mesh consisting of 56 elements of mesh size approximately $h=30$ mm is used. A two-level geometric refinement towards the corners and the free boundary at the center hole is applied to catch singularities, leading to a total number of 95 elements. The TDNNS method with $k=4$ is applied, which leads to 2593 coupling degrees of freedom.

The bending moments $\underline{\mathbf{m}}_{xy}$ and $\underline{\mathbf{m}}_{yy}$ are depicted in Figure 5 and Figure 6, respectively. Note that in Figure 5, different scales are used for the original plate and the zoom to the interior hole, such that the steep gradient of the bending moment becomes visible.