Numerical results for mimetic discretization of Reissner–Mindlin plate problems

Abstract

A low-order mimetic finite difference method for Reissner–Mindlin plate problems is considered. Together with the source problem, the free vibration and the buckling problems are investigated. Details about the scheme implementation are provided, and the numerical results on several different types of meshes are reported.

Appendix

In this section we describe how to build the local bilinear forms appearing in Sect. 3.

In what follows $m=m(E) \in {\mathbb N }$ will indicate the number of vertices of the polygon $E$. We number the vertices in counterclockwise sense as ${\mathsf v }_1,\ldots ,{\mathsf v }_m$ and analogously for the edges ${\mathsf e }_1,\ldots ,{\mathsf e }_m$, so that ${\mathsf v }_j$ and ${\mathsf v }_{j+1}$ are the endpoints of edge ${\mathsf e }_j$, $j=1,2,\ldots ,m$. Note that here and in the sequel all such indexes are considered modulus $m$, so that the index $m+1$ is identified with the index $1$. There are a total of $3m$ local degrees of freedom associated to each element of the mesh, three for each vertex. We order such local degrees of freedom first with all rotations and then all deflections, ordered as the vertices

$$\begin{aligned} \left\{ \, \varvec{\eta }_E^{{\mathsf v }_1}, \varvec{\eta }_E^{{\mathsf v }_2},\ldots ,\varvec{\eta }_E^{{\mathsf v }_m}, v_E^{{\mathsf v }_1}, v_E^{{\mathsf v }_2},\ldots , v_E^{{\mathsf v }_m} \, \right\} , \end{aligned}$$

where $( \varvec{\eta }_E, v_E) \in H_h|_E \times W_h |_E$.

The final local bilinear forms ${\mathbb M }={\mathbb M }(E) \in \mathbb R ^{3m \times 3m}$ associated to each element $E$ will be the sum of two parts

$$\begin{aligned} {\mathbb M }\ = \ {\mathbb M }_1 \ + \ \kappa t^{-2} {\mathbb M }_2 , \end{aligned}$$

(20)

the first one being associated to the $a_h(\cdot ,\cdot )$ term and the second one to the shear energy term. Once the elemental matrices ${\mathbb M }$ are built, the global stiffness matrix is implemented with a standard assembly procedure as in classical finite elements.

1.1 Matrix for the bilinear form $a_h(\cdot ,\cdot )$

We start from the bilinear form $a_h(\cdot ,\cdot )$, which is the sum of local bilinear forms that we express as matrices ${\mathsf M }={\mathsf M }(E) \in \mathbb R ^{2 m \times 2m}$

$$\begin{aligned} a_{h}^E(\varvec{\beta }_E,\varvec{\eta }_E) = \varvec{\beta }_E^T {\mathsf M }\varvec{\eta }_E {\quad \forall }E \in {\mathcal T }_h, \ \forall \varvec{\beta }_E,\varvec{\eta }_E \in H_h |_E. \end{aligned}$$

The first and main step is to build the matrix ${\mathsf M }$. With this purpose we introduce the matrices ${\mathsf N }={\mathsf N }(E)$ and ${\mathsf R }={\mathsf R }(E)$ in $\mathbb R ^{2 m \times 6}$. Note again that for ease of notation we do not make explicit the dependence on the involved matrices from $E$. Let $\mathbf q _1,\ldots ,\mathbf q _6$ be the following basis for the first order vector polynomials (with $2$ components) defined on $E$:

$$\begin{aligned} \mathbf q _1 = \begin{pmatrix} 1 \\ 0 \end{pmatrix} , \quad \mathbf q _2 = \begin{pmatrix} 0 \\ 1 \end{pmatrix} , \quad \mathbf q _3 = \begin{pmatrix} \bar{y}\\ -\bar{x}\end{pmatrix} , \quad \mathbf q _4 = \begin{pmatrix} \bar{y}\\ \bar{x}\end{pmatrix} , \ \mathbf q _5 = \begin{pmatrix} \bar{x}\\ \bar{y}\end{pmatrix} , \quad \mathbf q _6 = \begin{pmatrix} \bar{x}\\ -\bar{y}\end{pmatrix}. \end{aligned}$$

where we recall that $\bar{x},\bar{y}$ represent cartesian coordinates with the origin in the barycenter of the element. Then, the six columns ${\mathsf N }_1,\ldots ,{\mathsf N }_6$ of ${\mathsf N }$ are vectors in $\mathbb R ^{2m}$ defined by the interpolation of the polynomials $\mathbf q _1,\ldots ,\mathbf q _6$ into the space $H_h|_E$ (see (4))

$$\begin{aligned} {\mathsf N }_j = (\mathbf q _j)_{{\mathbf I }, E}\quad j=1,\ldots ,6. \end{aligned}$$

The columns of the matrix ${\mathsf N }$ thus represent the linear polynomials $\mathbf q _j$ written in terms of the degrees of freedom of $H_h|_E$.

The columns ${\mathsf R }_j$ of the matrix ${\mathsf R }$ are instead defined as the vectors in $\mathbb R ^{2m}$ associated to the right hand side of the consistency condition (S2 $_a$), computed with respect to the polynomials $\mathbf q _j$, $j=1,\ldots ,6$. In other words ${\mathsf R }_j$ is the unique vector in $\mathbb R ^{2m}$ such that for all $\varvec{\eta }_E \in H_h|_E \equiv \mathbb R ^{2m}$

$$\begin{aligned} ({\mathsf R }_j)^T \varvec{\eta }_E = \sum _{i=1}^m \left(\mathbb C \varvec{\varepsilon }({\mathbf q _j})\mathbf n _E^{{\mathsf e }_i} \right) \cdot \left(\frac{|{\mathsf e }_i|}{2}\left[\varvec{\eta }_E^{{\mathsf v }_i}+\varvec{\eta }_E^{{\mathsf v }_{i+1}}\right]\right) \end{aligned}$$

see equation (S2 $_a$). Note that, since $\varvec{\varepsilon }({\mathbf q _j}) = {\mathbf 0 }$ for $j=1,2,3$, the first three columns ${\mathsf R }_1,{\mathsf R }_2,{\mathsf R }_3$ of ${\mathsf R }$ have all zero entries.

From the definition of the vectors ${\mathsf N }_j$ and ${\mathsf R }_j$, it is clear that the consistency condition (S2 $_a$) translates into the algebraic condition

$$\begin{aligned} {\mathsf M }{\mathsf N }_j = {\mathsf R }_j j=1,\ldots ,6 \ \Leftrightarrow \ {\mathsf M }{\mathsf N }= {\mathsf R }. \end{aligned}$$

(21)

We therefore introduce the matrix ${\mathsf K }\in \mathbb R ^{6 \times 6}$ defined by

$$\begin{aligned} {\mathsf K }= {\mathsf N }^T {\mathsf R }= {\mathsf R }^T {\mathsf N }. \end{aligned}$$

It is easy to check that such matrix is symmetric and semi-positive definite. Moreover, it is of the form

$$\begin{aligned} {\mathsf K }= \begin{pmatrix} {\mathbf 0 }_{3 \times 3}&\quad {\mathbf 0 }_{3 \times 3} \\ {\mathbf 0 }_{3 \times 3}&\quad {\mathsf K }_\star \end{pmatrix} \end{aligned}$$

with ${\mathsf K }_\star$ positive definite. Therefore is it immediate to compute the pseudo inverse of ${\mathsf K }$

$$\begin{aligned} {\mathsf K }^\dagger = \begin{pmatrix} {\mathbf 0 }_{3 \times 3}&\quad {\mathbf 0 }_{3 \times 3} \\ {\mathbf 0 }_{3 \times 3}&\quad {\mathsf K }_\star ^{-1} \end{pmatrix}. \end{aligned}$$

We are now ready to define the local matrix ${\mathsf M }$. Let ${\mathsf P }$ be a projection on the space orthogonal to the columns of ${\mathsf N }$

$$\begin{aligned} {\mathsf P }= {\mathbb I }_{2m\times 2m} - {\mathsf N }({\mathsf N }^T {\mathsf N })^{-1} {\mathsf N }^T \end{aligned}$$

with ${\mathbb I }_{2m \times 2m}$ the identity matrix. We then set

$$\begin{aligned} {\mathsf M }= {\mathsf R }{\mathsf K }^\dagger {\mathsf R }^T + \alpha {\mathsf P }\end{aligned}$$

with $\alpha \in \mathbb R$ any positive number, typically scaled as the trace of the first part of the matrix. Then, it is immediate to check that ${\mathsf M }$ satisfies the consistency condition (21).

Finally, note that the matrix ${\mathsf M }\in \mathbb R ^{2m \times 2m}$ is defined only with respect to the rotation degrees of freedom, since the bilinear form $a_h(\cdot ,\cdot )$ is independent of the deflection variable. When it comes to build the local matrix ${\mathbb M }_1 \in \mathbb R ^{3m \times 3m}$ appearing in (20) one simply needs to introduce the restriction matrix ${\mathsf S } \in \mathbb R ^{3m \times 2m}$

$$\begin{aligned} {\mathsf S } = \begin{pmatrix} {\mathbb I }_{2m \times 2m} \\ {\mathbf 0 }_{m \times 2m} \end{pmatrix} \end{aligned}$$

and set

$$\begin{aligned} {\mathbb M }_1 = {\mathsf S } {\mathsf M }{\mathsf S }^T. \end{aligned}$$

 

1.2 Matrix for the shear term

The local matrices for the shear part are obtained as a product of matrices representing the operators and bilinear forms that appear in the second term of the left hand side of Method 1. We therefore start building a matrix $\overline{\mathsf M }=\overline{\mathsf M }(E) \in \mathbb R ^{m\times m}$ that represents the local scalar product

$$\begin{aligned}{}[\varvec{\gamma }_E , \varvec{\delta }_E ]_{\Gamma _h,E} = \varvec{\gamma }_E^T \overline{\mathsf M }\varvec{\delta }_E \quad \forall \varvec{\gamma }_E,\varvec{\delta }_E \in \Gamma _h|_E. \end{aligned}$$

We order the $m$ degrees of freedom of $\Gamma _h|_E$ as the edges of $E$. The construction follows the same philosophy as in the previous section and therefore is presented more briefly. Now, the two columns of the matrix $\overline{\mathsf N }\in \mathbb R ^{m \times 2}$ are defined by

$$\begin{aligned} \overline{\mathsf N }_j = ( \mathbf {curl} \,q_j)_{\text{ II},E} \quad j=1,2 , \end{aligned}$$

where the sub-index ${}_{II}$ represents the interpolation operator shown in (5) and $q_1,q_2$ denote the following basis of the (zero average) linear polynomials on $E$

$$\begin{aligned} q_1 = \bar{x}, \quad q_2 = \bar{y}. \end{aligned}$$

Analogously, the matrix $\overline{\mathsf R }\in \mathbb R ^{m \times 2}$ is defined through its columns as the right hand side of the consistency condition (16) in [11]

$$\begin{aligned} (\overline{\mathsf R }_j)^T \varvec{\delta }_E = - \sum _{i=1}^m \delta _E^{{\mathsf e }_i} \int _{{\mathsf e }_i} q_j \quad \forall j=1,2, \ \forall \varvec{\delta }_E \in \Gamma _h|_E \equiv \mathbb R ^m, \end{aligned}$$

where we neglected the $\mathrm{rot}_{\Gamma _{h}}$ part since $q_1$ and $q_2$ have zero average on $E$. Again, we need to introduce $\overline{\mathsf K }\in \mathbb R ^{2 \times 2}$ given by $\overline{\mathsf K }= \overline{\mathsf N }^T \overline{\mathsf R }= \overline{\mathsf R }^T \overline{\mathsf N }$ that is easily shown to be positive definite and symmetric. We can therefore finally set

$$\begin{aligned} \overline{\mathsf M }= \overline{\mathsf R }\, (\overline{\mathsf K })^{-1} \overline{\mathsf R }^T + \alpha \overline{\mathsf P }\end{aligned}$$

with $\alpha \in \mathbb R ^+$ and the projection matrix $\overline{\mathsf P }= {\mathbb I }_{m\times m} - \overline{\mathsf N }(\overline{\mathsf N }^T \overline{\mathsf N })^{-1} \overline{\mathsf N }^T.$ The consistency condition $\overline{\mathsf M }\overline{\mathsf N }=\overline{\mathsf R }$ follows by construction while the stability can be derived with the results in [20].

The local matrix ${\mathbb M }_2$ appearing in (20) can be built combining $\overline{\mathsf M }$ with a matrix ${\mathsf C }={\mathsf C }(E) \in \mathbb R ^{m \times 3m}$ representing the $\nabla _h$ and $\Pi _h$ operators that appear in Method 1. We therefore set ${\mathsf C }= \begin{pmatrix} - {\mathsf C }_1&{\mathsf C }_2 \end{pmatrix}$ with the matrix ${\mathsf C }_1 = {\mathsf C }_1(E) \in \mathbb R ^{m \times 2m}$ representing the $\Pi _h$ operator

$$\begin{aligned} {\mathsf C }_1 = \frac{1}{2} \begin{pmatrix} (\mathbf t _E^{{\mathsf e }_1})^T&\quad (\mathbf t _E^{{\mathsf e }_1})^T&\quad {\mathbf 0 }_{1\times 2}&\quad {\mathbf 0 }_{1\times 2}&\quad \ldots&\quad {\mathbf 0 }_{1\times 2} \\ {\mathbf 0 }_{1\times 2}&\quad (\mathbf t _E^{{\mathsf e }_2})^T&\quad (\mathbf t _E^{{\mathsf e }_2})^T&\quad {\mathbf 0 }_{1\times 2}&\quad \ldots&\quad {\mathbf 0 }_{1\times 2} \\ {\mathbf 0 }_{1\times 2}&\quad {\mathbf 0 }_{1\times 2}&\quad (\mathbf t _E^{{\mathsf e }_3})^T&\quad (\mathbf t _E^{{\mathsf e }_3})^T&\quad \ldots&\quad {\mathbf 0 }_{1\times 2} \\ \vdots&\quad \vdots&\quad \vdots&\quad \vdots&\quad \vdots&\quad \vdots \\ (\mathbf t _E^{{\mathsf e }_m})^T&\quad {\mathbf 0 }_{1\times 2}&\quad \ldots&\quad {\mathbf 0 }_{1\times 2}&\quad {\mathbf 0 }_{1\times 2}&\quad (\mathbf t _E^{{\mathsf e }_m})^T \end{pmatrix} , \end{aligned}$$

and the matrix ${\mathsf C }_2 = {\mathsf C }_2(E) \in \mathbb R ^{m \times m}$ representing the $\nabla _h$ operator

$$\begin{aligned} {\mathsf C }_2 = \begin{pmatrix} -|{\mathsf e }_1|^{-1}&\quad |{\mathsf e }_1|^{-1}&\quad 0&\quad 0&\quad \ldots&\quad 0 \\ 0&\quad -|{\mathsf e }_2|^{-1}&\quad |{\mathsf e }_2|^{-1}&\quad 0&\quad \ldots&\quad 0 \\ 0&\quad 0&\quad -|{\mathsf e }_3|^{-1}&\quad |{\mathsf e }_3|^{-1}&\quad \ldots&\quad 0 \\ \vdots&\quad \vdots&\quad \vdots&\quad \vdots&\quad \vdots&\quad \vdots \\ -|{\mathsf e }_m|^{-1}&\quad 0&\quad \ldots&\quad 0&\quad 0&\quad |{\mathsf e }_m|^{-1} \end{pmatrix}. \end{aligned}$$

Finally, the local matrices for the shear part are given by

$$\begin{aligned} {\mathbb M }_2 = {\mathsf C }^T \overline{\mathsf M }{\mathsf C }. \end{aligned}$$

1.3 Right hand sides

The loading term for the source problem in Method 1 follows immediately from (13). One gets the local right hand vectors ${\mathbf b }={\mathbf b }(E) \in \mathbb R ^{3m}$ defined by

$$\begin{aligned} {\mathbf b }_j = \left\{ \begin{array}{l@{\quad }l} 0&\text{ if} j=1,2,\ldots ,2m \\ \bar{g}|_E \, \omega _E^{(j-2m)}&\text{ if} j=2m+1,2m+2,\ldots ,3m , \end{array} \right. \end{aligned}$$

that are then assembled as usual into the global load vector.

The mass matrix for the free vibration problem in Method 2, associated to the bilinear form (14) is built again by a standard assembly procedure. The local mass matrices ${\mathsf D }= {\mathsf D }(E) \in \mathbb R ^{3m \times 3m}$ associated to the elemental mass bilinear forms

$$\begin{aligned} m_h^E(\varvec{\beta }_E, w_E ; \varvec{\eta }_E , v_E) = (\varvec{\beta }_E , w_E)^T {\mathsf D }\, (\varvec{\eta }_E , v_E) \end{aligned}$$

$\forall E \in {\mathcal T }_h, \forall \varvec{\beta }_E,\varvec{\eta }_E \in H_h |_E, \forall w_E, v_E \in W_h|_E$ are diagonal and defined by

$$\begin{aligned} {\mathsf D }_{ii} = \left\{ \begin{array}{l@{\quad } l} t^2 \omega _E^{\lceil i/2 \rceil } / 12&\text{ if} i=1,2,\ldots ,2m \\ \omega _E^{(i-2m)}&\text{ if} i=2m+1,2m+2,\ldots ,3m \end{array} \right. \end{aligned}$$

where the symbol $\lceil \ \rceil$ stands for a round up to the nearest integer.

The stress matrix for the buckling problem in Method 3, associated to the bilinear form (15) is also built as the sum of local matrices $\widehat{\mathsf B }=\widehat{\mathsf B }(E) \in \mathbb R ^{m \times m}$

$$\begin{aligned} b_h^E( w_E, v_E) = w_E^T \widehat{\mathsf B }\, v_E \qquad \forall w_E, v_E \in W_h|_E. \end{aligned}$$

Note that the symmetric tensor $\varvec{\sigma }\in \mathbb R ^{2 \times 2}$ that appears in (S2 $_b$) can have either rank 2 or rank 1. In order to build the matrix $\widehat{\mathsf B }$, we start introducing $\{\hat{q}_1,\hat{q}_2,\hat{q}_3 \}$ a basis for the linear polynomials on $E$, such that $\hat{q}_1=1$ and $\hat{q}_2,\hat{q}_3$ have zero integral on $E$. Moreover, if $\mathrm rank (\varvec{\sigma }) = 1$, we also require that $\nabla \hat{q}_2 \in \mathrm Ker (\varvec{\sigma })$. We then define as usual the auxiliary matrices $\widehat{\mathsf N }=\widehat{\mathsf N }(E) \in \mathbb R ^{m \times 3}$ and $\widehat{\mathsf R }=\widehat{\mathsf R }(E) \in \mathbb R ^{m \times 3}$ through its columns. We set

$$\begin{aligned} \widehat{\mathsf N }_j = (\hat{q}_j)_{\text{ I},E} ,\quad j=1,2,3 , \end{aligned}$$

where the sub-index ${}_{\text{ I}}$ denotes the interpolation operator in (3), and define $\widehat{\mathsf R }_j$ as the unique vector in $\mathbb R ^m$ such that

$$\begin{aligned} \widehat{\mathsf R }_j^T v_E = \sum _{i=1}^m \left(\varvec{\sigma }\nabla \hat{q}_j \cdot \mathbf n _E^{{\mathsf e }_i} \right) \frac{|{\mathsf e }_i|}{2}[v_E^{{\mathsf v }_i} + v_E^{{\mathsf v }_{i+1}}] \quad \forall j=1,2,3,\;\forall v_E \in W_h|_E \equiv \mathbb R ^m \end{aligned}$$

in accordance with (S2 $_b$). Note that clearly $\widehat{\mathsf R }_1$ is null, and that, if $\text{ rank}(\varvec{\sigma })=1$ also $\widehat{\mathsf R }_2$ is null. One then defines as usual the semi-positive definite and symmetric matrix $\widehat{\mathsf K }=\widehat{\mathsf K }(E) \in \mathbb R ^{3 \times 3}$ given by $\widehat{\mathsf K }= \widehat{\mathsf R }^T \widehat{\mathsf N }= \widehat{\mathsf N }^T \widehat{\mathsf R }.$ Since $\widehat{\mathsf K }$ is block diagonal, with the first block of zeros and the second invertible, it is immediate to compute the pseudo inverse matrix $\widehat{\mathsf K }^\dagger$, in a way similar to the one used for ${\mathsf K }$ in Sect. 6.1. Then, we introduce $\widehat{\mathsf B }=\widehat{\mathsf B }(E)\in \mathbb R ^{m \times m}$

$$\begin{aligned} \widehat{\mathsf B }= \widehat{\mathsf R }\, (\widehat{\mathsf K })^\dagger \widehat{\mathsf R }^T + \alpha \widehat{\mathsf P }\end{aligned}$$

with $\alpha \in \mathbb R$ non negative and the projection matrix $\widehat{\mathsf P }= {\mathbb I }_{m\times m} - \widehat{\mathsf N }(\widehat{\mathsf N }^T \widehat{\mathsf N })^{-1} \widehat{\mathsf N }^T$. Note that, since no global coercivity conditions are required, differently from the previous matrices also the choice $\alpha =0$ can be taken.

Finally, note that the matrix $\widehat{\mathsf B }\in \mathbb R ^{m \times m}$ is defined only with respect to the deflection degrees of freedom, since the bilinear form $b_h(\cdot ,\cdot )$ is independent of the rotation variable. The remaining entries in the assembled (right hand side) stress matrix associated to Method 3 can be simply filled with zeros.