Finite Elements for the One Variable Version of Mindlin-Reissner Plate

Abstract

To analyze thin and thick plates, the paper presents two rectangular finite elements with high accuracy. In these elements, the proposed formulations of the displacement field utilize the Bergan-Wang approach, which depends only on one variable: the plate lateral deflection. This approach ensures that shear-locking problem will not happen as thickness decreases. The degrees of freedom of the proposed elements are twenty-four for the first element and it is named BWRE24, while the second one has thirty-six degrees of freedom and is named BWRE36. To evidence the efficiency of the two elements, a series of numerical examples for an isotropic plate subjected to various loadings and with different boundary conditions have been analyzed. Very good results are obtained suffering no numerical difficulties in case of very thin plates.

Keywords:
Bergan-Wang approach; Thin and Thick plate; Displacement finite element formulation

Graphical Abstract

1. Introduction

Plates have wide applications in different engineering fields. There are many theories for the analysis of plate bending problems. The classical plate theory (CPT) based on the assumption that the perpendicular line to the middle plane of the un-deformed plate remains straight and perpendicular to the deformed mid-surface. The only variable in CPT is the lateral deflection of the plate middle plane. Because of neglecting the effect of shear deformation, the (CPT) applies mainly to thin plates [¹, ²].

For moderately thick and thick plates with thickness to length ratio over 0.05, it is necessary to take into account the effect of transverse shear deformations [³]. The mostly used theory in engineering applications is that of Reissner - Mindlin plate theory. They modified the Kirchhoff's normality assumption by allowing uniform rotation (constant shear- Strain and stress) for the normal line about the middle surface after deformation [⁴].

According to Reissner-Mindlin plate theory, There are three independent variables: the lateral deflection $w$ of the plate middle plane and the two rotations of the normal line to the plate ${\theta }_x$ and${\theta }_y$. These three variables are governed by a coupled system of partial differential equations of six order. Therefore, it is necessary to satisfy three boundary conditions along the edges of a plate while only two boundary conditions are required in the classical plate theory [⁵-⁷].

Despite the simple formulation of The Reissner-Mindlin theory, the discretization by means of finite elements based on standard low-order schemas show an intense lack of convergence when the thickness of the plate is too small compared with the other dimensions of the plate. This unwanted phenomenon, known as shear locking, can be overcome by reducing the effect of the shear energy [⁸, ⁹].

P.G. Bergan and X. Wang [¹⁰] H. Abdalla and K. Hassan [¹¹] have employed an approximation that have led to general expressions of rotations of the normal as functions of displacement. Therefore, the strain energy is function of only one variable $w.$ Second, third, and fourth derivatives of this variable appear in the strain energy and consequently, finite elements of displacement type require polynomials of order larger than or equal to four. We notice that P.G. Bergan and X. Wang [¹⁰] used quadrilateral finite element, which is only complete for a polynomial of degree three.

Liu et al. [¹²] have also followed the same approximation. They present the formulation and assessment of a quadrilateral finite element for thin and moderately thick plates (NCQS). The NCQS element appears to be a simple and consistent 16-degree-of-freedom (DOF) conforming quadrilateral plate-bending element, but also complete for a third-degree polynomial only.

In this paper, new rectangular finite elements of a complete, forth degree polynomial and fifth degree polynomial based on Bergan - Wang principal approximation. The elements have been named as BWRE24 and BWRE36 respectively. We have tested the performance of these elements on different boundary condition of thin and thick plate.

For the reader's convenience, we start by briefly summarizing the theoretical formulation of the strain energy expression introduced by Bergan-Wang approximation for the case of isotropic material. Followed by the introduced displacement finite elements with the formulation of corresponding stiffness matrices. Accuracy illustrated on several numerical examples. Work is in progress to extend the formulation to anisotropic plates.

2. The governing equations of the first order shear plate theory

In Mindlin-Reissner theory, the normal on the middle un-deformed plane of the plate may rotate without remaining normal as shown in Fig. 1-a. Therefore, the displacement fields [¹³] are.

$u\left(x.y.z\right)=z{\theta }_x=z\left(-\frac{\partial w}{\partial x}+{\gamma }_x\right)$ (1.a)

$v\left(x.y.z\right)=z{\theta }_y=z\left(-\frac{\partial w}{\partial y}+{\gamma }_y\right)$ (1.b)

$w\left(x.y.z\right)=w$ (1.c)

Where,${\theta }_x$ and${\theta }_y$ are the rotations of the transverse line about x and y-axes, respectively. Moreover, ${\gamma }_x$ and${\gamma }_y$ are the shear strain of the mid-plane of the plate.

The strain deformation vector of the plate due to bending and shear can be written as.

$\left\{ϵ\right\}={\left\{{ϵ}_x{ϵ}_y{ϵ}_{xy}{\gamma }_x{\gamma }_y\right\}}^{\text{T}}={\left\{z\frac{\partial {\theta }_x}{\partial x}z\frac{\partial {\theta }_y}{\partial y}z\left(\frac{\partial {\theta }_x}{\partial y}+\frac{\partial {\theta }_y}{\partial x}\right){\gamma }_x{\gamma }_y\right\}}^{\text{T}}$ (2)

The stress vector due to bending and shear effects depending on stress-strain relationship in the isotropic homogenous elastic plate can be written as.

$\left\{\sigma \right\}={\left\{{\sigma }_x{\sigma }_y{\sigma }_{xy}{\tau }_{xz}{\tau }_{yz}\right\}}^T=\left[D\right]\{\int \}$ (3)

Where, $\left[D\right]$ the flexure bending and shear rigidity of isotropic elastic plate is defined as.

$\left[D\right]=\left[\begin{array}{cc}\left[D_b\right]& 0\\ 0& \left[D_s\right]\end{array}\right];$ (4)

$\left[D_b\right]=\frac{E{h}^3}{12\left(1-{\nu }^2\right)}\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& \frac{1-\nu }{2}\end{array}\right];\left[D_s\right]=Gkh\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$ (5)

In Eq. (5), $k$ is the shear correction factor, which is taken as 5/6 [⁶], $\nu$ is the Poisson’s ratio, $E$ is the Young’s modulus and $G$ is the modulus of rigidity.

By integrating the planar stresses$\sigma$ in Eq. (3) over the thickness of the plate, bending moments and twisting moment per unit length in Fig.1-b can obtained as

$\left[\begin{array}{c}{M}_x\\ M_y\\ M_{xy}\end{array}\right]=\underset{-h/2}{\overset{h/2}{{\displaystyle \int }}}z\left[\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\sigma }_{xy}\end{array}\right]dz=\frac{E{h}^3}{12\left(1-{\nu }^2\right)}\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& \frac{1-\nu }{2}\end{array}\right]\left[\begin{array}{c}\frac{\partial {\theta }_x}{{\partial }_x}\\ \frac{\partial {\theta }_y}{{\partial }_y}\\ \frac{\partial {\theta }_x}{{\partial }_y}+\frac{\partial {\theta }_y}{{\partial }_x}\end{array}\right]$ (6)

Moreover, the foundational equations of the shear force per unit length in Fig. 1-b are

$\left[\begin{array}{c}{Q}_x\\ Q_y\end{array}\right]=\underset{-h/2}{\overset{h/2}{{\displaystyle \int }}}\left[\begin{array}{c}{\tau }_{xz}\\ {\tau }_{yz}\end{array}\right]dz=Gkh\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}{\gamma }_x\\ {\gamma }_y\end{array}\right]$ (7)

Moreover, according to the principle of virtual work [¹⁴], we get:

$\underset{\Omega e}{{\displaystyle \int }}\left[M_x\frac{\partial \delta {\theta }_x}{\partial x}+M_y\frac{\partial \delta {\theta }_y}{\partial y}+M_{xy}\left(\frac{\partial \delta {\theta }_x}{\partial y}+\frac{\partial \delta {\theta }_y}{\partial x}\right)+Q_x\left(\delta {\theta }_x+\frac{\partial \delta w}{\partial x}\right)+Q_y\left(\delta {\theta }_y+\frac{\partial \delta w}{\partial y}\right)-\delta wP\right]dxdy-\underset{\Gamma e}{{\displaystyle \oint }}(\delta {\theta }_n{M}_n+\delta {\theta }_s{M}_{ns}+\delta w{Q}_n)\cdot ds=0$ (8)

Figure 1
(a) Un-deformed and deformed geometries of an edge of the plate according to Reissner and Mindlin theories. (b) Bending moments, shear forces and edge forces acting on the plate.

The governing differential equations of plate can written as

$\frac{\partial Q_x}{\partial x}+\frac{\partial Q_y}{\partial y}=-P$ (9.a)

$\frac{\partial M_x}{\partial x}+\frac{\partial M_{xy}}{\partial y}=Q_x$ (9.b)

$\frac{\partial M_y}{\partial y}+\frac{\partial M_{xy}}{\partial x}=Q_y$ (9.c)

3. Bergan-Wang energy expression.

Bergan and Wang obtained a mathematical expression for the shear strains in terms of the transverse displacement only by substituting the shear forces$Q_x$ and $Q_y$ from Eq. (9.b and 9.c) into Eq. (7) the shear strain${\gamma }_x$ can written as [¹²]:

${\gamma }_x=\frac{h^2}{5\left(1-\nu \right)}\left[\frac{{\partial }^{3}w}{\partial x^3}+\frac{{\partial }^{3}w}{\partial x\partial y^2}+\frac{{\partial }^2{\gamma }_x}{\partial x^2}+\frac{\left(1-\nu \right)}{2}\frac{{\partial }^2{\gamma }_x}{\partial y^2}+\frac{\left(1+\nu \right)}{2}\frac{{\partial }^2{\gamma }_y}{\partial x\partial y}\right]$ (10)

A similar expression can be obtained for ${\gamma }_y.$

Bergan and Wang assumed that, the shear deformations of the plate${\gamma }_x$ and${\gamma }_y$ can be approximated by piecewise linear polynomials. This assumption is suitable for finite element technique and leads to the following expressions [¹⁰]:

${\gamma }_x\cong -h_0^2\left(\frac{{\partial }^{3}w}{\partial x^3}+\frac{{\partial }^{3}w}{\partial x\partial y^2}\right)$ (11.a)

${\gamma }_y\cong -h_0^2\left(\frac{{\partial }^{3}w}{\partial y^3}+\frac{{\partial }^{3}w}{\partial y\partial x^2}\right)$ (11.b)

Where,

$h_0^2=\frac{h^2}{6k\left(1-\nu \right)}$ (12)

Substituting Eq. (11.a and 11.b), into the expressions for the bending and twisting moments (Eq. (6) and Eq. (7)) leads to:

$\left[\begin{array}{c}{M}_x\\ M_y\\ M_{xy}\end{array}\right]=\left[D_b\right]\left[\begin{array}{c}{w}_{.xx}+h_0^2\left(w_{.xxxx}+w_{.xxyy}\right)\\ w_{.yy}+h_0^2\left(w_{.yyyy}+w_{.xxyy}\right)\\ 2w_{.xy}+2h_0^2\left(w_{.xyyy}+w_{.xxxy}\right)\end{array}\right]$ (13)

$\left[\begin{array}{c}{Q}_x\\ Q_y\end{array}\right]=\left[D_s\right]\left[\begin{array}{c}-h_0^2\left(w_{.xxx}+w_{.xyy}\right)\\ -h_0^2\left(w_{.yyy}+w_{.xxy}\right)\end{array}\right]$ (14)

Finally, the strain energy expression $U$ per unit area, which is the sum of bending energy and shear energy according to Bergan and Wang approach [¹⁰, ¹¹], can be written as:

$U=U_b+U_s$ (15)

Where,

$U_b=\frac{1}{2}{\left[\begin{array}{c}{w}_{.xx}+h_0^2\left(w_{.xxxx}+w_{.xxyy}\right)\\ w_{.yy}+h_0^2\left(w_{.yyyy}+w_{.xxyy}\right)\\ 2w_{.xy}+2h_0^2\left(w_{.xyyy}+w_{.xxxy}\right)\end{array}\right]}^T\left[D_b\right]\left[\begin{array}{c}{w}_{.xx}+h_0^2\left(w_{.xxxx}+w_{.xxyy}\right)\\ w_{.yy}+h_0^2\left(w_{.yyyy}+w_{.xxyy}\right)\\ 2w_{.xy}+2h_0^2\left(w_{.xyyy}+w_{.xxxy}\right)\end{array}\right]$ (16)

$U_s=\frac{1}{2}{\left[\begin{array}{c}-h_0^2\left(w_{.xxx}+w_{.xyy}\right)\\ -h_0^2\left(w_{.yyy}+w_{.xxy}\right)\end{array}\right]}^T\left[D_s\right]\left[\begin{array}{c}-h_0^2\left(w_{.xxx}+w_{.xyy}\right)\\ -h_0^2\left(w_{.yyy}+w_{.xxy}\right)\end{array}\right]$ (17)

4. Formulation of finite element method.

The basic thought of the present approach in Eq. (16) and Eq. (17) can achieved by using a higher order polynomial approximation (higher than the third degree) of the displacement field. Therefore, two new rectangular finite elements of a complete fourth and fifth degree polynomial [¹⁵, ¹⁶] have used in the following plate bending problems and named as BWRE24 and BWRE36 respectively.

4.1 Derivation of element transverse displacement function

In both above plate elements, the transverse displacement function $w$ can expressed as generalized displacement modes

$w=w\left(x.y\right)=\left[N\right]\left\{{\beta }_i\right\}$ (18)

The number of terms for polynomial $\left[N\right]$ must be equal to the total number of degrees of freedom that each element has. The transverse displacement coefficients vector $\left\{{\beta }_i\right\}$ (where, $i$ from 1 to number of degrees of freedom of one element) can be evaluated by substituting with the local coordinates of each node of the element into the degrees of freedom mathematical formula ${\left\{{\psi }_i\right\}}^T$ for each plate element BWRE24 and BWRE36. This leads to a square matrix $\left[C\right]$ of size $i\times i$, consisting of numbers only [¹⁷]. Finally, the unknown coefficients vector $\left\{{\beta }_i\right\}$ can obtained from

$\left\{{\beta }_i\right\}={\left[C\right]}^{-1}\left\{{\psi }_i\right\}$ (19)

Therefore, the transverse displacement function $w$ can written as

$w=w\left(x.y\right)=\left[N\right]\left\{{\beta }_i\right\}=\left[N\right]{\left[C\right]}^{-1}\left\{{\psi }_i\right\}$ (20)

4.2 Derivation of element stiffness matrices and load vector

According to the strain energy expression (Eq. (16), (17)) and the generalized strains in Eq. (2), the generalized stiffness matrix of a thick isotropic plate element given by a bending part plus a shear part [¹⁸].

$\left[k^e\right]=\left[k_b^e\right]+\left[k_s^e\right]$ (21)

Where,

$\left[k_{\text{b}}^e\right]={\int }_{{\Omega }_{\text{e}}}{\left[{ϵ}_{\text{b}}\right]}^{\text{T}}\left[{\text{D}}_{\text{b}}\right]\left[{ϵ}_{\text{b}}\right]d\text{A};\left[k_{\text{s}}^e\right]={\int }_{{\Omega }_{\text{e}}}{\left[{ϵ}_{\text{s}}\right]}^{\text{T}}\left[{\text{D}}_{\text{s}}\right]\left[{ϵ}_{\text{s}}\right]d\text{A}$ (22)

The generalized bending strain vector $\left[{ϵ}_{\text{b}}\right]$ and shear strain vector $\left[{ϵ}_{\text{s}}\right]$ can written according to the approach presented in Eq. (13, 14) and according to the finite element formulation of the transverse displacement function $w$ in Eq. (20) as

$\left[S_b\right]=\left[\begin{array}{c}{N}_{.xx}+h_0^2\left(N_{.xxxx}+N_{.xxyy}\right)\\ N_{.yy}+h_0^2\left(N_{.yyyy}+N_{.xxyy}\right)\\ 2N_{.xy}+2h_0^2\left(N_{.xyyy}+N_{.xxxy}\right)\end{array}\right]\left\{{\beta }_i\right\};\left[S_s\right]=\left[\begin{array}{c}-h_0^2\left(N_{.xxx}+N_{.xyy}\right)\\ -h_0^2\left(N_{.yyy}+N_{.xxy}\right)\end{array}\right]\left\{{\beta }_i\right\}$ (23)

Finally, the finite element formulation of the element stiffness matrix is expressed by:

$\left[k_b^e\right]={\left\{{\beta }_i\right\}}^T\left(\iint {\left[S_b\right]}^T\left[D_b\right]\left[S_b\right]dxdy\right)\left\{{\beta }_i\right\}$ (24)

$\left[k_s^e\right]={\left\{{\beta }_i\right\}}^T\left(\iint {\left[S_s\right]}^T\left[D_s\right]\left[S_s\right]dxdy\right)\left\{{\beta }_i\right\}$ (25)

Similarly, according to the principle of virtual work (Eq. (8)), and the finite element formulation of the transverse displacement function $w$ (Eq. (20)), the finite element formulation of the element load vector [¹⁹] can written as

$\left\{F^e\right\}=-{\left\{{\beta }_i\right\}}^T\underset{{\Omega }_e}{{\displaystyle \iint }}{\left[N\right]}^{T}P\left(x.y\right)dxdy$ (26)

After assembling the elements stiffness matrix and the element force vector to obtain the formation of the overall plate stiffness matrix $\left[\mathbb{K}\right]$ and the resultant systems force vector$\left\{\mathbb{F}\right\}$, the final equilibrium equation can written as

$\left\{\mathbb{F}\right\}=\left[\mathbb{K}\right]\left\{\psi \right\}$ (27)

5. The selection of the displacement field.

Consider a rectangular plate of dimensions $\left(a\times b\right)$ and thickness (h) have been meshed by two models of rectangular finite element which the origin of its local coordinate system is located at the lower-left node of the element.

5.1 Finite element formulation of BWRE24

The first finite element model, BWRE24 has four nodes and the degrees of freedom at each node are the lateral deflection$w$, two rotations ${\theta }_x$ and${\theta }_y$, the two degrees of freedom ${\delta }_{M_x}$ and ${\delta }_{M_y}$ which are proportional to the internal bending moments per unit length $M_x$ and $M_y$, and finally ${\delta }_{M_{xy}}$ which is proportional to the internal twisting moment per unit length$M_{xy}$. So, the generalized displacement field at each node $\left(i\right)$ of the element can expressed as:

${\left\{{\psi }_i\right\}}^T=\left\{\begin{array}{cccccc}w& {\theta }_x& {\theta }_y& {\delta }_{M_x}& {\delta }_{M_y}& {\delta }_{M_{xy}}\end{array}\right\}$ (28)

$\left\{{\psi }_i\right\}=\left\{\begin{array}{c}w\\ -w_{.x}-h_0^2\left(w_{.xxx}+w_{.xyy}\right)\\ -w_{.y}-h_0^2\left(\left(w_{.yyy}+w_{.yxx}\right)\right)\\ w_{.xx}+h_0^2\left(w_{.xxxx}+w_{.xxyy}\right)\\ w_{.yy}+h_0^2\left(w_{.yyyy}+w_{.xxyy}\right)\\ 2w_{.xy}+2h_0^2\left(w_{.xyyy}+w_{.xxxy}\right)\end{array}\right\}$ (29)

Because of the six degrees of freedom at each node, only twenty-four terms of Pascal's triangle are selected [²⁰], Such that as $h_0^2$ tends to zero BWRE24 will be as 24-dof displacement type finite element model [²⁰]. Therefore, the deflection function as a complete polynomial of degree four can expressed as:

$w\left(x.y\right)=\left[\begin{array}{c}1\\ xy\\ x^{2}xy{y}^2\\ x^3{x}^{2}yx{y}^2{y}^3\\ x^4{x}^{3}y{x}^2{y}^{2}x{y}^3{y}^4\\ x^5{x}^2{y}^3{x}^3{y}^2{y}^5\\ \left(x^4\left(\frac{yb}{2}-\frac{2}{3}\frac{y^3}{b}\right)\right)\left(y^4\left(\frac{xa}{2}-\frac{2}{3}\frac{x^3}{a}\right)\right)\\ x^3{y}^3\\ \left(x^5\left(\frac{yb}{2}-\frac{2}{3}\frac{y^3}{b}\right)\right)\left(y^5\left(\frac{xa}{2}-\frac{2}{3}\frac{x^3}{a}\right)\right)\end{array}\right]\left\{\begin{array}{c}{\beta }_1\\ {\beta }_2\\ {\beta }_3\\ ⋮\\ ⋮\\ {\beta }_{24}\end{array}\right\}$ (30)

5.2 Finite element formulation of BWRE36

For the second finite element model, BWRE36 has nine nodes as shown in Fig. 2 and the degrees of freedom at each node are the lateral deflection$w$, two rotations ${\theta }_x$ and${\theta }_y$, and ${\delta }_{M_{xy}}$ which is proportional to the internal twisting moment per unit length$M_{xy}$. So, the generalized displacement field at each node $\left(i\right)$ of the element can be expressed as:

${\left\{{\psi }_i\right\}}^T=\left\{\begin{array}{cccc}w& {\theta }_x& {\theta }_y& {\delta }_{M_{xy}}\end{array}\right\}$ (31)

$\left\{{\psi }_i\right\}=\left\{\begin{array}{c}w\\ -w_{.x}-h_0^2\left(w_{.xxx}+w_{.xyy}\right)\\ -w_{.y}-h_0^2\left(\left(w_{.yyy}+w_{.yxx}\right)\right)\\ 2w_{.xy}+2h_0^2\left(w_{.xyyy}+w_{.xxxy}\right)\end{array}\right\}$ (32)

Because of the four degrees of freedom at each node, only thirty-six terms of Pascal's triangle have selected, Such that as $h_0^2$ tends to zero, BWRE36 will be as the finite element model in [¹⁶]. Therefore, the deflection function as a complete polynomial of degree five can expressed as:

$w\left(x.y\right)=\left[\begin{array}{c}1\\ xy\\ x^{2}xy{y}^2\\ x^3{x}^{2}yx{y}^2{y}^3\\ x^4{x}^{3}y{x}^2{y}^{2}x{y}^3{y}^4\\ x^5{x}^{4}y{x}^3{y}^2{x}^2{y}^{3}x{y}^4{y}^5\\ x^{5}y{x}^4{y}^2{x}^3{y}^3{x}^2{y}^{4}x{y}^5\\ x^5{y}^2{x}^4{y}^3{x}^3{y}^4{x}^2{y}^5\\ x^5{y}^3{x}^4{y}^4{x}^3{y}^5\\ x^5{y}^4{x}^4{y}^5\\ x^5{y}^5\end{array}\right]\left\{\begin{array}{c}{\beta }_1\\ {\beta }_2\\ {\beta }_3\\ ⋮\\ ⋮\\ {\beta }_{36}\end{array}\right\}$ (33)

Figure 2
BWRE36 rectangular plate element

6. Numerical tests

6.1 Bending of a square plate under uniformly distributed load with various boundary conditions

In order to verify the validity of the two finite element models BWRE24 and BWRE36, four cases of a thin and moderately thick square plate of side length $a=1$,having $\nu =$0.3, subjected to uniform load of intensity $q$ have been studied.

6.1.1 The results shown in Table 1 for case of two opposite edges simply supported and the other two edges clamped SCSC. With an illustration of the geometric shape of deflection and slopes for the plate at a certain case when h/a=0.1 as shown in Fig.3. And, comparing the plate’s central deflection and slopes results from BWRE24 and BWRE36 with these results from the Reissner-Mindlin plate theory as shown in Fig.4 and Fig.5

Table 1
The non-dimensional central deflections and bending moments for a uniformly distributed load with two opposite edges simply supported and the others are clamped SCSC.

Non-dimensional deflection $\overline{W}=\left(w_{\left(\frac{a}{2}.\frac{b}{2}\right)}D/q{a}^4\right)$
h/a	0.001	0.01	0.05	0.1	0.15	0.2
Reissner theory [21]	0.00192	0.00192	0.00199	0.00220	0.00254	0.00298
Mindlin theory [4]	0.00192	0.00192	0.00199	0.00221	0.00256	0.00302
BWRE36 (16x16)	0.001917	0.001926	0.002038	0.002240	0.002545	0.002981
BWRE24 (32x32)	0.001917	0.001944	0.002113	0.002327	0.002645	0.003139
Thin Plate theory [22]:0.0019171 for uniform load
Non-dimensional bending moment ${\overline{M}}_x=\left(M_x{}_{\left(\frac{a}{2}.\frac{b}{2}\right)}/q{a}^2\right)$
h/a	0.001	0.01	0.05	0.1	0.15	0.2
Reissner theory [21]	0.02439	0.02441	0.02491	0.02638	0.02863	0.03144
Mindlin theory [4]	0.02439	0.02440	0.02476	0.02579	0.02733	0.02921
BWRE24 (32x32)	0.02439	0.02473	0.02610	0.02567	0.02479	0.2448
Thin Plate theory [22]:0.024387 for uniform load
Non-dimensional bending moment ${\overline{M}}_y=\left(M_y{}_{\left(\frac{a}{2}.\frac{b}{2}\right)}/q{a}^2\right)$
h/a	0.001	0.01	0.05	0.1	0.15	0.2
Reissner theory [21]	0.03324	0.03325	0.03343	0.03395	0.03475	0.03579
Mindlin theory [4]	0.03324	0.03325	0.03326	0.03327	0.03321	0.03305
BWRE24 (32x32)	0.03323	0.03374	0.03529	0.03471	0.3373	0.03334
Thin Plate theory [22]:0.03324 for uniform load

Figure 3
Deflection w and slopes ${\theta }_x$ and ${\theta }_y$ of SCSC plate with h/a=0.1and b=a

Figure 4
Comparison between R-M, BWRE24 and BWRE36 in Deflection w of SCSC plate at midpoint with different h/a ratio.

Figure 5
Comparison between R-M, BWRE24 and BWRE36 in rotations ${\theta }_x$ and ${\theta }_y$ of SCSC plate at centroidal axes.

6.1.2 The results shown in Table 2 for case of two opposite edges simply supported and the other two edges free SFSF. With an illustration of the geometric shape of deflection and slopes for the plate at a certain case when h/a=0.1 as shown in Fig.6. And, comparing the plate’s central deflection and slopes results from BWRE24 and BWRE36 with these results from the Reissner-Mindlin plate theory as shown in Fig.7 and Fig.8

Table 2
The non-dimensional central deflections and bending moments for a uniformly distributed load with two opposite edges simply supported and the others are free SFSF.

Non-dimensional deflection $\overline{W}=\left(w_{\left(\frac{a}{2}.\frac{b}{2}\right)}D/q{a}^4\right)$
h/a	0.001	0.01	0.05	0.1	0.15	0.2
Reissner theory [21]	-	-	-	0.01341	0.01379	0.01433
Mindlin theory [4]	-	-	-	0.01346	0.01391	0.01454
BWRE36 (16x16)	0.01309	0.01309	0.01317	0.01344	0.01388	0.01441
BWRE24 (32x32)	0.01309	0.01309	0.01314	0.01350	0.01420	0.01513
Thin Plate theory [22]:0.013093 for uniform load
Non-dimensional bending moment ${\overline{M}}_x=\left(M_x{}_{\left(\frac{a}{2}.\frac{b}{2}\right)}/q{a}^2\right)$
h/a	0.001	0.01	0.05	0.1	0.15	0.2
Reissner theory [21]	0.1225	0.1225	0.1225	0.1225	0.1227	0.1230
Mindlin theory [4]	0.1225	0.1225	0.1225	0.1225	0.1226	0.1229
BWRE24 (32x32)	0.1225	0.1223	0.1217	0.1224	0.1225	0.1226
Thin Plate theory [22]:0.1225 for uniform load
Non-dimensional bending moment ${\overline{M}}_y=\left(M_y{}_{\left(\frac{a}{2}.\frac{b}{2}\right)}/q{a}^2\right)$
h/a	0.001	0.01	0.05	0.1	0.15	0.2
Reissner theory [21]	0.02704	0.02695	0.02641	0.02561	0.02466	0.02353
Mindlin theory [4]	0.02704	0.02695	0.02641	0.02564	0.02474	0.02372
BWRE24 (32x32)	0.02700	0.02649	0.0250	0.02660	0.02701	0.02703
Thin plate theory [22]: 0.0270 for uniform load

Figure 6
Deflection w and slopes ${\theta }_x$ and ${\theta }_y$ of SFSF plate with h/a=0.1and b=a

Figure 7
Comparison between R-M, BWRE24 and BWRE36 in Deflection w of SFSF plate at midpoint with different h/a ratio.

Figure 8
Comparison between R-M, BWRE24 and BWRE36 in rotations ${\theta }_x$ and ${\theta }_y$ of SFSF plate at centroidal axes.

6.1.3 The results shown in Table 3 for case of three edges simply supported and the other edge clamped SCSS. With an illustration of the geometric shape of deflection and slopes for the plate at a certain case when h/a=0.1 as shown in Fig.9. And, comparing the plate’s central deflection and slopes results from BWRE24 and BWRE36 with these results from the Reissner-Mindlin plate theory as shown in Fig.10 and Fig.11

Table 3
The non-dimensional central deflections and bending moments for a uniformly distributed load with two opposite edges simply supported and the others are clamped SCSS.

Non-dimensional deflection $\overline{W}=\left(w_{\left(\frac{a}{2}.\frac{b}{2}\right)}D/q{a}^4\right)$
h/a	0.001	0.01	0.05	0.1	0.15	0.2
R-M theory [23]	0.27835	0.27863	0.28535	0.30574	0.33852	0.38272
BWRE36 (16x16)	0.27857	0.28013	0.29518	0.31496	0.34500	0.38440
BWRE24 (32x32)	0.27856	0.28355	0.30855	0.32561	0.35203	0.39960
Thin Plate theory [22]:0.2786 for uniform load
Non-dimensional bending moment ${\overline{M}}_x=\left(M_x{}_{\left(\frac{a}{2}.\frac{b}{2}\right)}/q{a}^2\right)$
h/a	0.001	0.01	0.05	0.1	0.15	0.2
BWRE24 (32x32)	0.03916	0.03450	0.04227	0.04123	0.03980	0.03931
Thin Plate theory [22]:0.039 for uniform load
Non-dimensional bending moment ${\overline{M}}_y=\left(M_y{}_{\left(\frac{a}{2}.\frac{b}{2}\right)}/q{a}^2\right)$
h/a	0.001	0.01	0.05	0.1	0.15	0.2
BWRE24 (32x32)	0.03385	0.03990	0.03674	0.03576	0.03440	0.03397
Thin Plate theory [22]:0.034 for uniform load

Figure 9
Deflection w and Slopes ${\theta }_x$ and ${\theta }_y$ of SCSS plate with h/a=0.1 and b=a

Figure 10
Comparison between R-M, BWRE24 and BWRE36 in Deflection w of SCSS plate at midpoint with different h/a ratio.

Figure 11
Comparison between R-M, BWRE24 and BWRE36 in rotations ${\theta }_x$ and ${\theta }_y$ of SCSS plate at centroidal axes.

6.1.4 The results shown in Table 4 for case of three edges simply supported and the other edge free SFSS. With an illustration of the geometric shape of deflection and slopes for the plate at a certain case when h/a=0.1 as shown in Fig.12. And, comparing the plate’s central deflection and slopes results from BWRE24 and BWRE36 with these results from the Reissner-Mindlin plate theory as shown in Fig.13 and Fig.14

Table 4
The non-dimensional central deflections and bending moments for a uniformly distributed load with two opposite edges simply supported and the others are clamped SFSS.

Non-dimensional deflection $\overline{W}=\left(w_{\left(a.\frac{b}{2}\right)}D/q{a}^4\right)$
h/a	0.001	0.01	0.05	0.1	0.15	0.2
R-M theory [23]	0.01284	0.01287	0.01308	0.01349	0.01407	0.01481
BWRE36 (16x16)	0.01285	0.01292	0.01357	0.01388	0.01432	0.01492
BWRE24 (32x32)	0.01285	0.01312	0.01415	0.01432	0.01448	0.01503
Thin Plate theory [22]: 0.01286 for uniform load
Non-dimensional bending moment ${\overline{M}}_x=\left(M_x{}_{\left(a.\frac{b}{2}\right)}/q{a}^2\right)$
h/a	0.001	0.01	0.05	0.1	0.15	0.2
BWRE24 (32x32)	0.00031	0.00076	0.00072	0.00037	0.00085	0.00021
Non-dimensional bending moment ${\overline{M}}_y=\left(M_y{}_{\left(a.\frac{b}{2}\right)}/q{a}^2\right)$
h/a	0.001	0.01	0.05	0.1	0.15	0.2
BWRE24 (32x32)	0.1129	0.1139	0.1209	0.1185	0.1146	0.1129
Thin Plate theory [22]: 0.112 for uniform load

Figure 12
Deflection w and Slopes ${\theta }_x$ and ${\theta }_y$ of SFSS plate with h/a=0.1and b=a

Figure 13
Comparison between R-M, BWRE24 and BWRE36 in Deflection w of SFSS plate at midpoint with different h/a ratio.

Figure 14
Comparison between R-M, BWRE24 and BWRE36 in rotations ${\theta }_x$ and ${\theta }_y$ of SFSS plate about x and y axes that passes through point ($a.\frac{b}{2}$).

6.2 C. Chinosi and C. Lovadina Problem

Both finite element models BWRE24 and BWRE36 have performed on a model problem for which the exact solution explicitly known [²⁴]. This test consists of a unitary square block [0, 1]² with clamped boundary conditions on all four sides and a distributed load $q\left(x.y\right)$given by the following function:

$q\left(x.y\right)=\frac{E.h^3}{12\left(1-{\upsilon }^2\right)}\left[\left(12y\left(y-1\right)\left(5x^2-5x+1\right)\left(2y^2{\left(y-1\right)}^2\right)\left(5y^2-5y+1\right)x\left(x-1\right)\right)+\left(12x\left(x-1\right)\left(5y^2-5y+1\right)\left(2x^2{\left(x-1\right)}^2\right)+(y\left(y-1\right)\left(5x^2-5x+1\right)\right)\right]$ (34)

The exact solution for the displacement, rotations and bending moments have given in [²⁴, ²⁵] as:

$w\left(x.y\right)=\frac{1}{3}{x}^3{\left(x-1\right)}^3{y}^3{\left(y-1\right)}^3-\frac{2h^2}{5\left(1-\upsilon \right)}\left[y^3{\left(y-1\right)}^{3}x\left(x-1\right)\left(5x^2-5x+1\right)+x^3{\left(x-1\right)}^{3}y\left(y-1\right)\left(5y^2-5y+1\right)\right]$ (35)

${\theta }_x\left(x.y\right)=x^2\left(2x-1\right){\left(x-1\right)}^2{y}^3{\left(y-1\right)}^3$ (36)

${\theta }_y\left(x.y\right)=-y^2\left(2y-1\right){\left(y-1\right)}^2{x}^3{\left(x-1\right)}^3$ (37)

$M_{xx}\left(x.y\right)=\left(\frac{-E{h}^3}{6\left(1-{\nu }^2\right)}\right)\left(\left(y^3{\left(y-1\right)}^3\left(x-x^2\right)\left(5x^2-5x+1\right)\right)+\left(\nu (x^3{\left(x-1\right)}^3\left(y-y^2\right)\left(5y^2-5y+1\right)\right)\right)$ (38)

$M_{yy}\left(x.y\right)=\left(\frac{-E{h}^3}{6\left(1-{\nu }^2\right)}\right)\left(\left(x^3{\left(x-1\right)}^3\left(y-y^2\right)\left(5y^2-5y+1\right)\right)+\left(\nu (y^3{\left(y-1\right)}^3\left(x-x^2\right)\left(5x^2-5x+1\right)\right)\right)$ (39)

$M_{xy}\left(x.y\right)=M_{yx}\left(x.y\right)=\left(\frac{-3E{h}^3}{12\left(1+\nu \right)}\right)\left(2x-1\right)x^2{\left(x-1\right)}^2{y}^2{\left(y-1\right)}^2\left(2y-1\right)$ (40)

In the following table 5, 6 this plate problem is solved according to Reissner-Mindlin theory and due to Bergan-Wang approach by the two elements BWRE24 and BWRE36, With an illustration of the geometric shape of deflection and slopes for the plate at a certain case when h/a=0.1 as shown in Fig. 15. For each study a thin plate with h = 0.001 and a thick plate with h = 0.1 are considered, and comparing the plate’s central deflection and slopes results from BWRE24 and BWRE36 with these results from the Reissner-Mindlin plate theory as shown in Fig. 16 and Fig. 17 In order to test the performance of each element and investigate the presence of shear locking.

Table 5
The central deflections, rotation and bending moments for a functional distributed load with clamped boundary conditions on all four sides (h/a = 0.001), $E=$10.92$*{10}^9$ and $\nu =$0.3

Exact deflection at the midpoint of the plate $w_{\left(0.5.0.5\right)}=$8.138$e^{-5}$
No. of elements	2 X 2	4 X 4	8 X 8	16 X 16	32 X 32	40 X 40
R-M theory [23]	-4.185$e^{-7}$	-2.430$e^{-5}$	4.943$e^{-5}$	7.306$e^{-5}$	7.936$e^{-5}$	8.0124$e^{-5}$
BWRE24	-1.883$e^{-4}$	-1.956$e^{-6}$	5.815$e^{-5}$	7.547$e^{-5}$	8.088$e^{-5}$	8.2211$e^{-5}$
BWRE36	---	-2.353$e^{-5}$	5.573$e^{-5}$	7.540$e^{-5}$	8.026$e^{-5}$	8.0991$e^{-5}$
Exact slope of the plate at the midpoint (0.5, 0.5) ${\theta }_x={\theta }_y=0$
No. of elements	2 X 2	4 X 4	8 X 8	16 X 16	32 X 32	40 X 40
R-M theory[23]	1.058$e^{-22}$	1.4010$e^{-4}$	1.067$e^{-5}$	5.907$e^{-5}$	7.936$e^{-5}$	8.012$e^{-5}$
BWRE24	-1.883$e^{-4}$	-1.956$e^{-6}$	5.815$e^{-5}$	7.547$e^{-5}$	3.029$e^{-5}$	2.431$e^{-5}$
BWRE36	---	4.164$e^{-13}$	1.796$e^{-14}$	2.307$e^{-13}$	3.837$e^{-13}$	3.572$e^{-13}$
Exact internal twisting moment at midpoint of the plate $M_{xy}\left(0\cdot 5.0\cdot 5\right)=0$
No. of elements	2 X 2	4 X 4	8 X 8	16 X 16	32 X 32	40 X 40
BWRE24	0.0041	0.0011	4.085$e^{-4}$	1.331$e^{-4}$	8.134$e^{-5}$	8.933$e^{-5}$
BWRE36	---	1.280$e^{-13}$	3.635$e^{-15}$	9.412$e^{-14}$	2.217$e^{-11}$	7.037$e^{-12}$
Exact internal bending moment at midpoint of the plate $M_x=M_y=0\cdot 0025$
No. of elements	2 X 2	4 X 4	8 X 8	16 X 16	32 X 32	40 X 40
BWRE24	0.0030	1.172e-04	0.00140	0.00180	0.00194	0.001998

Table 6
The central deflections, rotation and bending moments for a functional distributed load with clamped boundary conditions on all four sides (h/a = 0.1), $E=$ 10.92$*{10}^3$ and $\nu =$0.3

Exact deflection at the midpoint of the plate $w_{\left(0.5.0.5\right)}=$9.25$e^{-5}$
No. of elements	2 X 2	4 X 4	8 X 8	16 X 16	32 X 32	40 X 40
R-M theory[23]	-4.185$e^{-5}$	-2.552$e^{-5}$	5.717$e^{-5}$	8.328$e^{-5}$	9.020$e^{-5}$	9.104 $e^{-5}$
BWRE24	-2.457$e^{-4}$	-8.717$e^{-6}$	6.965$e^{-5}$	9.247$e^{-5}$	9.744$e^{-5}$	9.761$e^{-5}$
BWRE36	---	-2.723$e^{-5}$	6.285 $e^{-5}$	8.664$e^{-5}$	9.254$e^{-5}$	9.330$e^{-5}$
Exact slope of the plate at the midpoint (0.5, 0.5) ${\theta }_x={\theta }_y=0$
No. of elements	2 X 2	4 X 4	8 X 8	16 X 16	32 X 32	40 X 40
R-M theory[23]	2.504$e^{-4}$	1.435$e^{-4}$	1.161$e^{-4}$	6.415$e^{-5}$	3.287$e^{-5}$	2.638$e^{-5}$
BWRE24	2.032$e^{-4}$	1.932$e^{-4}$	1.257$e^{-4}$	6.562$e^{-5}$	3.131$e^{-5}$	2.439$e^{-5}$
BWRE36	---	2.238$e^{-14}$	1.604$e^{-15}$	2.856$e^{-16}$	4.627$e^{-14}$	4.429$e^{-14}$
Exact internal twisting moment at midpoint of the plate $M_{xy}\left(0\cdot 5.0\cdot 5\right)=0$
No. of elements	2 X 2	4 X 4	8 X 8	16 X 16	32 X 32	40 X 40
BWRE24	0.0041	0.0011	4.085$e^{-4}$	1.332$e^{-4}$	8.134 $e^{-5}$	8.9337$e^{-5}$
BWRE36	---	1.281$e^{-13}$	3.636$e^{-15}$	9.412$e^{-14}$	2.218$e^{-11}$	7.037$e^{-12}$
Exact internal bending moment at midpoint of the plate $M_x=M_y=0\cdot 0025$
No. of elements	2 X 2	4 X 4	8 X 8	16 X 16	32 X 32	40 X 40
BWRE24	0.0035	3.5242$e^{-4}$	0.00150	0.00190	0.001991	0.0020

Figure 15
Deflection w and slopes ${\theta }_x$ and ${\theta }_y$ of CCCC plate with h/a=0.1and b=a

Figure 16
Comparison between exact solution, R-M, BWRE24 and BWRE36 in Deflection w of CCCC plate at midpoint with h/a =0.001 and h/a = 0.1.

Figure 17
Comparison between R-M, BWRE24 and BWRE36 in rotations ${\theta }_x$ and ${\theta }_y$ of CCCC plate at centroidal axes.

7. Conclusion

In this paper, the isotropic thick plate-bending problem has been analyzed by using two new elements. The formulations of these elements is based on a general expression for the strain energy of the plate as a function of plate deflection only. Therefore, a complete polynomial of fourth degree for the first element (BWRE24) and a complete polynomial of fifth degree for the second element (BWRE36) have taken in the constructing shape functions for both elements. It emphasized that the shear-locking problem in the thin plates has been avoided.

According to this hypothesis, the internal displacement fields for both thin and thick plate have calculated under different loads by using the same three boundary conditions of Mindlin Reissner plate theory. Numerical results show that the two elements offer good accuracy and convergence rate for both thin and thick plate bending problem. And of course, using complete polynomial of degree five (BWRE36 element) has a better performance than using complete polynomial of only degree four (BWRE24).

Applying this approach to orthotropic material and composite material is a relatively easy task. For general domain, triangular elements such as Bell and Argyris elements [¹⁵] are good candidates, one needs to modify them using the same degrees of freedom as BWRE24.

Acknowledgment

The authors thank professors: H. Asks of University of Sheffield, F. F. Mahmoud of Zagazig University, H. Abd Allah of Military Technical College and M. El Shazly of British University in Egypt for useful discussions and comments.

References

• ^[1]
  Gou, Y., Y. Cai, and H. Zhu, A Simple High-Order Shear Deformation Triangular Plate Element with Incompatible Polynomial Approximation. Applied Sciences, 2018. 8(6): p. 975.
• ^[2]
  Noor, A.K.; Burton, W.S. Assessment of shear deformation theories for multilayered composite plates.Appl. Mech. Rev. 1989, 42, 1-13.
• ^[3]
  Fuh-Gwo, Y. and R.E. Miller, 1988, A rectangular finite element for moderately thick flat plates. Computers & structures. 30(6): p. 1375-1387.
• ^[4]
  Lee, K., G. Lim, and C. Wang, Thick Lévy plates re-visited. International Journal of Solids and Structures, 2002. 39(1): p. 127-144.
• ^[5]
  Soh, A.-K., et al., 2001, A new twelve DOF quadrilateral element for analysis of thick and thin plates. European Journal of Mechanics-A/Solids. 20(2): p. 299-326.
• ^[6]
  Reissner, E. The effect of transverse shear deformation on the bending of elastic plates. J. Appl. Mech. 1945, 12, 69-77.
• ^[7]
  Mindlin, R.D. Influence of rotatory inertia and shear on flexural motions of isotropic elastic plates. J. Appl. Mech. 1951, 18, 31-38.
• ^[8]
  Da Veiga, L.B., et al., An isogeometric method for the Reissner-Mindlin plate bending problem. Computer Methods in Applied Mechanics and Engineering, 2012. 209: p. 45-53.
• ^[9]
  F. Auricchio, C. Lovadina, Analysis of kinematic linked interpolation methods for Reissner-Mindlin plate problems, Comput. Methods Appl. Mech. Engrg. 190 (2001) 2465-2482.
• ^[10]
  Bergan, P.G. and X. Wang, Quadrilateral plate bending elements with shear deformations. Computers & Structures, 1984. 19(1-2): p. 25-34.
• ^[11]
  Abdalla, H. and K. Hassan, On the Bergan‐Wang approach for moderately thick plates. Communications in applied numerical methods, 1988. 4(1): p. 51-58.
• ^[12]
  Ine-Wei Liu, Tienfuan Kerh, Chien-Chang Lin, A conforming quadrilateral plate bending element with shear deformations, Computers & structures, 1995, 56(1), 93-100.
• ^[13]
  Ozkul, T.A. and U. Ture, The transition from thin plates to moderately thick plates by using finite element analysis and the shear locking problem. Thin-Walled Structures, 2004, 42(10): p. 1405-1430.
• ^[14]
  Reddy, J.N., Energy principles and variational methods in applied mechanics, John Wiley & Sons, 2001.
• ^[15]
  P.G. Ciarlet, The finite element method for elliptic problems, North Holland, 1975.
• ^[16]
  O.C. Zienkiewicz, R. Taylor, J.Z. Zhu, The Finite element method: its basis and Fundamentals, 7 edition, Butterworth-Heinemann, 2013.
• ^[17]
  AboElsooud M.T., Vibration control of plates using periodically distributed shunted piezoelectric patches. (Doctoral dissertation) University of Maryland, College Park. 2003.
• ^[18]
  Oñate, E., Structural analysis with the finite element method. Linear statics: volume 2: beams, plates and shells, Springer Science & Business Media, 2013.
• ^[19]
  Ugural, A.C., Stresses in plates and shells, 2nd revised edition,McGraw-London, 1998.
• ^[20]
  Popplewell, N. and D. McDonald, Conforming rectangular and triangular plate-bending elements. Journal of Sound and Vibration, 1971. 19(3): p. 333-347.
• ^[21]
  Wang, C.M., et al., Relationships between bending solutions of Reissner and Mindlin plate theories. Engineering structures, 2001. 23(7): p. 838-849.
• ^[22]
  Timoshenko, S.P. and S. Woinowsky-Krieger, Theory of plates and shells. McGraw-hill, 1959.
• ^[23]
  Siva Srinivas Kolukula, isoparametric Q4 elements have used to dicretize the plate according to Mindlin-Reissner plate theory, Indira Gandhi Center for Atomic Research, http://www.mathworks.de/matlabcentral/fileexchange/32029-plate-bending, 2011.
  » http://www.mathworks.de/matlabcentral/fileexchange/32029-plate-bending
• ^[24]
  C. Chinosi, C. Lovadina, Numerical analysis of some mixed finite element methods for Reissner-Mindlin plates, Comput. Mech. 16 (1995), 36-44.
• ^[25]
  J. Kiendla, F. Auricchio, L. Beir˜ao da Veigac, C. Lovadinad, A. Realia, Isogeometric collocation methods for the Reissner-Mindlin plate problem, Computer Methods in Applied Mechanics and Engineering, 284 (2015), 489-507.

Edited by

Editor:

Marco L. Bittencourt.

Publication Dates

• Publication in this collection
  18 Sept 2020
• Date of issue
  2020

History

• Received
  26 June 2020
• Reviewed
  13 July 2020
• Accepted
  04 Aug 2020
• Published
  08 Aug 2020