Static, free vibration, and buckling analysis of plates using strain-based Reissner–Mindlin elements

Abstract

A quadrilateral and a triangular element based on the strain approach are developed for static, free vibration and buckling analyses of Reissner–Mindlin plates. The four-node triangular element SBTP4 has the three essential external degrees of freedom at each of the three corner nodes and at a mid-side node; whereas the quadrilateral element SBQP has the same degrees of freedom at each of the four corner nodes. Both elements use the same assumed strain functions which are in the linear variation where bending and transverse shear strains are independent and satisfy the compatibility equations. The use of the strain approach allows obtaining elements with higher-order terms for the displacements field. The formulated elements have been proposed to improve the strain-based rectangular plate element SBRP previously published. Several numerical examples demonstrate that the present elements are free of shear locking and provide high-accuracy results compared to the available published numerical and analytical solutions.

Introduction

Analyses of static, buckling and free vibration of plate structures play a large role in structural engineering applications. Considerable research works on analysis of plates are still being conducted (Mackerle 1997, 2002; Leissa 1969, 1987; Liew et al. 1995, 2004).

Designers prefer low-order Reissner–Mindlin plate elements due to their simplicity and efficiency. However, for thin plates, these elements often suffer from the shear locking phenomenon when dealing with thin plates. To overcome shear locking, many research works have been undertaken where the use of the selective reduced integration was first intervened (Zienkiewicz et al. 1971; Hughes et al. 1978; Malkus and Hughes 1978). The formulation procedure used is to divide the strain energy into two parts, one of bending and the other of shear. Then, two different integration rules for these two parts are used. For low-order polynomial elements based on displacement model, such as the four-node classical bilinear element, an exact integration (two Gauss points in each direction) is taken for the bending strain energy; whereas a reduced integration (one Gauss point) is used for the shear strain energy. This selective integration can be provided with a more efficient element but often leads to numerical instability. Considerable investigations have been oriented to develop robust elements using different improved formulations and numerical techniques to avoid shear locking such as mixed formulation, enhanced assumed strain methods, assumed natural strain methods, discrete shear gap method and smoothed finite element method (Lee and Wong 1982; Ayad et al. 1998; Lovadina 1998; César de Sá and Natal Jorge 1999; César de Sá et al. 2002; Cardoso et al. 2008; MacNeal 1982; Bathe and Dvorkin 1985, 1986; Zienkiewicz et al. 1990; Batoz and Katili 1992; Bletzinger et al. 2000; Nguyen-Xuan et al. 2008; Liu and Nguyen-Thoi 2010).

The strain approach has been employed as an alternative to formulate robust plate elements (Belarbi and Charif 1999; Belounar and Guenfoud 2005; Belounar and Guerraiche 2014; Guerraiche et al. 2018; Belounar et al. 2018) to increase the accuracy and stability of the numerical solutions as well as to eliminate shear locking phenomena. The use of the strain approach (Belarbi and Charif 1999; Belounar and Guenfoud 2005; Belounar and Guerraiche 2014; Guerraiche et al. 2018; Belounar et al. 2018; Djoudi and Bahai 2004a, b; Rebiai and Belounar 2013; 2014) has several advantages where it enables to obtain efficient elements with high-order polynomial terms for the displacement functions without the need of including internal nodes. The first developed strain-based Mindlin plate element SBRP (Belounar and Guenfoud 2005) has been adopted for the linear analysis of plates having only rectangular shapes. However, this element suffers from shear locking for very thin plates (Belounar et al. 2018). Then, the formulation of a new three-node strain-based triangular Mindlin plate element SBTMP (Belounar et al. 2018) has been developed for static and free vibration of plate bending. The assumed curvatures and transverse shear strains for the SBRP element (Belounar and Guenfoud 2005) are coupled and contain quadratic terms. The key idea used in this paper is to formulate new elements to overcome shear locking for very thin plates and to improve the accuracy for plates with regular and distorted shapes.

In this paper, a quadrilateral and a triangular strain-based plate element have been formulated for static, free vibration and buckling analyses of plates using Reissner–Mindlin theory. The opportunity is taken to explore the displacements field obtained from the strain-based quadrilateral plate element (SBQP) by applying it to a four-node triangular element strain-based triangular plate with four nodes (SBTP4) having the same degrees of freedom (W, β_x, and β_y) at each of the three corner nodes and a mid-side node. In the process of formulation, these elements are based on linear variation for the five strain components where bending and transverse shear strains are independent and satisfying the compatibility equations. The numerical study shows that the SBQP and SBTP4 elements pass the patch test, are free of shear locking, and can be found numerically more efficient than the SBRP element (Belounar and Guenfoud 2005).

Formulation of the proposed elements

Derivation of the displacements field

For Reissner–Mindlin plate elements (Fig. 1), the strains in terms of the displacements are given as:

$$\kappa_{x} = \frac{{\partial \beta_{x} }}{\partial x},\quad \kappa_{y} = \frac{{\partial \beta_{y} }}{\partial y},\quad \kappa_{xy} = \left( {\frac{{\partial \beta_{x} }}{\partial y} + \frac{{\partial \beta_{y} }}{\partial x}} \right),\quad \gamma_{xz} = \beta_{x} + \frac{\partial W}{\partial x},\quad \gamma_{yz} = \beta_{y} + \frac{\partial W}{\partial y}.$$

(1a)

Fig. 1

Quadrilateral and triangular Reissner–Mindlin plate elements

In matrix form, it can be given as

$$\left\{ {\begin{array}{*{20}l} {\kappa_{x} } \hfill \\ {\kappa_{y} } \hfill \\ {\kappa_{xy} } \hfill \\ {\gamma_{xz} } \hfill \\ {\gamma_{yz} } \hfill \\ \end{array} } \right\} = \left[ {\begin{array}{*{20}l} 0 \hfill & {{\partial \mathord{\left/ {\vphantom {\partial {\partial x}}} \right. \kern-0pt} {\partial x}}} \hfill & 0 \hfill \\ 0 \hfill & 0 \hfill & {{\partial \mathord{\left/ {\vphantom {\partial {\partial y}}} \right. \kern-0pt} {\partial y}}} \hfill \\ 0 \hfill & {{\partial \mathord{\left/ {\vphantom {\partial {\partial y}}} \right. \kern-0pt} {\partial y}}} \hfill & {{\partial \mathord{\left/ {\vphantom {\partial {\partial x}}} \right. \kern-0pt} {\partial x}}} \hfill \\ {{\partial \mathord{\left/ {\vphantom {\partial {\partial x}}} \right. \kern-0pt} {\partial x}}} \hfill & 1 \hfill & 0 \hfill \\ {{\partial \mathord{\left/ {\vphantom {\partial {\partial y}}} \right. \kern-0pt} {\partial y}}} \hfill & 0 \hfill & 1 \hfill \\ \end{array} } \right]\left\{ {\begin{array}{*{20}c} W \\ {\beta_{x} } \\ {\beta_{y} } \\ \end{array} } \right\}.$$

(1b)

The five strains, bending (κ_x, κ_y and κ_xy) and transverse shear (γ_xz and γ_yz), given in Eq. (1a) cannot be considered independent, for they are in terms of the displacements W, β_x and β_y and therefore, they must satisfy the compatibility equations (Belounar and Guenfoud 2005) given as:

$$\frac{{\partial^{2} \kappa_{x} }}{{\partial y^{2} }} + \frac{{\partial^{2} \kappa_{y} }}{{\partial x^{2} }} = \frac{{\partial^{2} \kappa_{xy} }}{\partial x\partial y},\quad \frac{{\partial^{2} \gamma_{xz} }}{\partial x\partial y} - \frac{{\partial^{2} \gamma_{yz} }}{{\partial x^{2} }} + \frac{{\partial \kappa_{xy} }}{\partial x} = 2\frac{{\partial \kappa_{x} }}{\partial y},\quad \frac{{\partial^{2} \gamma_{yz} }}{\partial x\partial y} - \frac{{\partial^{2} \gamma_{xz} }}{{\partial y^{2} }} + \frac{{\partial \kappa_{xy} }}{\partial y} = 2\frac{{\partial \kappa_{y} }}{\partial x}.$$

(2)

The field of displacements due to the three rigid body modes is obtained by having Eq. (1a) equal to zero and the following results are obtained:

$$W = \alpha_{1} - \alpha_{2} x - \alpha_{3} y,\quad \beta_{x} = \alpha_{2} ,\quad \beta_{y} = \alpha_{3} .$$

(3)

The proposed quadrilateral and triangular elements (SBQP and SBTP4) have three degrees of freedom (W, β_x and β_y) at each of the four nodes. Therefore, the displacements field should contain twelve independent constants and having used three (α₁, α₂, α₃) for the representation of the rigid body modes, the remaining nine constants (α₄, α₅, …, α₁₂) are to be apportioned among the five assumed strains of the two elements.

The interpolation of the assumed strains field for the present elements (SBQP and SBTP4) is given as:

$$\begin{aligned} & \kappa_{x} = \alpha_{4} + \alpha_{5} y,\,\kappa_{y} = \alpha_{6} + \alpha_{7} x,\quad \kappa_{xy} = \alpha_{8} + (2\alpha_{5} x) + (2\alpha_{7} y), \\ & \gamma_{xz} = \alpha_{9} + \alpha_{10} y,\quad \gamma_{yz} = \alpha_{11} + \alpha_{12} x. \\ \end{aligned}$$

(4)

Assumed bending (κ_x, κ_y and κ_xy) and transverse shear (γ_xz and γ_yz) strains given in Eq. (4) of the proposed elements are independent and have only linear terms contrarily for the SBRP element (Belounar and Guenfoud 2005), where bending and transverse shear strains are coupled and quadratic terms are included in the assumed shear strain components.

The bracketed terms of the assumed strains (Eq. 4) are added to have the compatibility equations (Eq. 2) to be satisfied. The strain functions (κ_x, κ_y, κ_xy, γ_xz, γ_yz) given by Eq. (4) are substituted into Eq. (1a) and after integration, we obtain:

$$\begin{aligned} & W = - \,\alpha_{4} \frac{{x^{2} }}{2} - \alpha_{5} \frac{{x^{2} y}}{2} - \alpha_{6} \frac{{y^{2} }}{2} - \alpha_{7} \frac{{xy^{2} }}{2} - \alpha_{8} \frac{xy}{2} + \alpha_{9} \frac{x}{2} + \alpha_{10} \frac{xy}{2} + \alpha_{11} \frac{y}{2} + \alpha_{12} \frac{xy}{2} \\ & \beta_{x} = \alpha_{4} x + \alpha_{5} xy + \alpha_{7} \frac{{y^{2} }}{2} + \alpha_{8} \frac{y}{2} + \alpha_{9} \frac{1}{2} + \alpha_{10} \frac{y}{2} - \alpha_{12} \frac{y}{2} \\ & \beta_{y} = \alpha_{5} \frac{{x^{2} }}{2} + \alpha_{6} y + \alpha_{7} xy + \alpha_{8} \frac{x}{2} - \alpha_{10} \frac{x}{2} + \alpha_{11} \frac{1}{2} + \alpha_{12} \frac{x}{2}. \\ \end{aligned}$$

(5a)

The displacement functions obtained from Eq. (5a) are summed to the displacements of rigid body modes given by Eq. (3) to obtain the final displacement shape functions:

$$\begin{aligned} & W = \alpha_{1} - \alpha_{2} x - \alpha_{3} y - \alpha_{4} \frac{{x^{2} }}{2} - \alpha_{5} \frac{{x^{2} y}}{2} - \alpha_{6} \frac{{y^{2} }}{2} - \alpha_{7} \frac{{xy^{2} }}{2} - \alpha_{8} \frac{xy}{2} + \alpha_{9} \frac{x}{2} + \alpha_{10} \frac{xy}{2} + \alpha_{11} \frac{y}{2} + \alpha_{12} \frac{xy}{2} \\ & \beta_{x} = \alpha_{2} + \alpha_{4} x + \alpha_{5} xy + \alpha_{7} \frac{{y^{2} }}{2} + \alpha_{8} \frac{y}{2} + \alpha_{9} \frac{1}{2} + \alpha_{10} \frac{y}{2} - \alpha_{12} \frac{y}{2} \\ & \beta_{y} = \alpha_{3} + \alpha_{5} \frac{{x^{2} }}{2} + \alpha_{6} y + \alpha_{7} xy + \alpha_{8} \frac{x}{2} - \alpha_{10} \frac{x}{2} + \alpha_{11} \frac{1}{2} + \alpha_{12} \frac{x}{2}. \\ \end{aligned}$$

(5b)

The displacement functions of Eq. (5b) and the strain functions of Eq. (4) can be given in matrix form, respectively, as:

$$\{ U\} = [P]\{ \alpha \} = [N]\{ q_{\text{e}} \} ,$$

(6)

$$\{ \varepsilon \} = [Q]\{ \alpha \} = [B]\{ q_{\text{e}} \}$$

(7)

$${\text{Where}}\,[N] = [P][C]^{ - 1}, \quad [B] = [Q][C]^{ - 1}.$$

(8)

And the matrices [P] and [Q] are given as:

$$[P] = \left[ {\begin{array}{*{20}c} 1 & { - x} & { - y} & { - \frac{{x^{2} }}{2}} & { - \frac{{x^{2} y}}{2}} & { - \frac{{y^{2} }}{2}} & { - \frac{{xy^{2} }}{2}} & { - \frac{xy}{2}} & {\frac{x}{2}} & {\frac{xy}{2}} & {\frac{y}{2}} & {\frac{xy}{2}} \\ 0 & 1 & 0 & x & {xy} & 0 & {\frac{{y^{2} }}{2}} & {\frac{y}{2}} & {\frac{1}{2}} & {\frac{y}{2}} & 0 & { - \frac{y}{2}} \\ 0 & 0 & 1 & 0 & {\frac{{x^{2} }}{2}} & y & {xy} & {\frac{x}{2}} & 0 & { - \frac{x}{2}} & {\frac{1}{2}} & {\frac{x}{2}} \\ \end{array} } \right],$$

(9)

$$[Q] = \left[ {\begin{array}{*{20}c} 0 & 0 & 0 & 1 & y & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & x & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & {(2x)} & 0 & {(2y)} & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & y & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & x \\ \end{array} } \right].$$

(10)

And the displacements field, the strains field, and constant parameters vectors are:

$$\{ U\} = \{ W,\beta_{x} ,\beta_{y} \}^{T} ,\quad \{ \varepsilon \} = \{ \kappa_{x} ,\kappa_{y} ,\kappa_{xy} ,\gamma_{xz} ,\gamma_{yz} \}^{T} ,\quad \{ \alpha \} = \{ \alpha_{1} ,\alpha_{2} , \ldots ,\alpha_{12} \}^{T} .$$

(11)

The geometrical strains can be expressed as:

$$\{ \varepsilon^{\text{g}} \} = \left[ {\begin{array}{*{20}c} {\frac{\partial }{\partial x}} & 0 & 0 \\ {\frac{\partial }{\partial y}} & 0 & 0 \\ 0 & {\frac{\partial }{\partial x}} & 0 \\ 0 & {\frac{\partial }{\partial y}} & 0 \\ 0 & 0 & {\frac{\partial }{\partial x}} \\ 0 & 0 & {\frac{\partial }{\partial y}} \\ \end{array} } \right]\left\{ {\begin{array}{*{20}c} W \\ {\beta_{x} } \\ {\beta_{y} } \\ \end{array} } \right\} = \left[ {\begin{array}{*{20}c} {\frac{\partial }{\partial x}} & 0 & 0 \\ {\frac{\partial }{\partial y}} & 0 & 0 \\ 0 & {\frac{\partial }{\partial x}} & 0 \\ 0 & {\frac{\partial }{\partial y}} & 0 \\ 0 & 0 & {\frac{\partial }{\partial x}} \\ 0 & 0 & {\frac{\partial }{\partial y}} \\ \end{array} } \right][P]\{ \alpha \} ,$$

(12)

$${\text{where}}\,[G] = \left[ {\begin{array}{*{20}c} {\frac{\partial }{\partial x}} & 0 & 0 \\ {\frac{\partial }{\partial y}} & 0 & 0 \\ 0 & {\frac{\partial }{\partial x}} & 0 \\ 0 & {\frac{\partial }{\partial y}} & 0 \\ 0 & 0 & {\frac{\partial }{\partial x}} \\ 0 & 0 & {\frac{\partial }{\partial y}} \\ \end{array} } \right][P].$$

We substitute Eq. (6) into Eq. (12), we obtain:

$$\{ \varepsilon^{\text{g}} \} = [G]\{ \alpha \} = [B^{\text{g}} ]\{ q_{\text{e}} \} ,$$

(13)

$${\text{where}}\,[B^{\text{g}} ] = [G][C]^{ - 1} .$$

(14)

And the matrix [G] is given as:

$$[G] = \left[ {\begin{array}{*{20}c} 0 & { - 1} & 0 & { - x} & { - xy} & 0 & { - \frac{{y^{2} }}{2}} & { - \frac{y}{2}} & {\frac{1}{2}} & {\frac{y}{2}} & 0 & {\frac{y}{2}} \\ 0 & 0 & { - 1} & 0 & { - \frac{{x^{2} }}{2}} & { - y} & { - xy} & { - \frac{x}{2}} & 0 & {\frac{x}{2}} & {\frac{1}{2}} & {\frac{x}{2}} \\ 0 & 0 & 0 & 1 & y & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & x & 0 & y & {\frac{1}{2}} & 0 & {\frac{1}{2}} & 0 & { - \frac{1}{2}} \\ 0 & 0 & 0 & 0 & x & 0 & y & {\frac{1}{2}} & 0 & { - \frac{1}{2}} & 0 & {\frac{1}{2}} \\ 0 & 0 & 0 & 0 & 0 & 1 & x & 0 & 0 & 0 & 0 & 0 \\ \end{array} } \right].$$

(15)

The transformation matrix [C] which relates the element nodal displacements ({q_e}^T = (W₁, β_x1, β_y1, …, W₄, β_x4, β_y4)) to the 12 constants ({α}^T = (α₁, …, α₁₂)) can be given as:

$$\left\{ {q_{e} } \right\} = \left[ C \right]\left\{ \alpha \right\}$$

(16a)

The constant parameters vector {α} can be derived from Eq. (16a) as follows:

$$\{ \alpha \} = [C]^{ - 1} \{ q_{\text{e}} \} .$$

(16b)

The matrices [N] (Eq. 8), [B] (Eq. 8) and [B_g] (Eq. 14) are obtained, respectively, by substituting Eq. (16b) into Eqs. (6), (7) and (13):

$${\text{Where}}\,\,[C] = \left[ {\begin{array}{*{20}c} {[P_{1} ]} & {[P_{2} ]} & {[P_{3} ]} & {[P_{4} ]} \\ \end{array} } \right]^{T} .$$

(17)

And the matrix [P_i] calculated from Eq. (9) for each of the four element nodes coordinates (x_i, y_i), (i = 1, 2, 3, 4) to obtain:

$$[P_{i} ] = \left[ {\begin{array}{*{20}c} 1 & { - x_{i} } & { - y_{i} } & { - \frac{{x_{i}^{2} }}{2}} & { - \frac{{x_{i}^{2} y_{i} }}{2}} & { - \frac{{y_{i}^{2} }}{2}} & { - \frac{{x_{i} y_{i}^{2} }}{2}} & { - \frac{{x_{i} y_{i} }}{2}} & {\frac{{x_{i} }}{2}} & {\frac{{x_{i} y_{i} }}{2}} & {\frac{{y_{i} }}{2}} & {\frac{{x_{i} y_{i} }}{2}} \\ 0 & 1 & 0 & {x_{i} } & {x_{i} y_{i} } & 0 & {\frac{{y_{i}^{2} }}{2}} & {\frac{{y_{i} }}{2}} & {\frac{1}{2}} & {\frac{{y_{i} }}{2}} & 0 & { - \frac{{y_{i} }}{2}} \\ 0 & 0 & 1 & 0 & {\frac{{x_{i}^{2} }}{2}} & {y_{i} } & {x_{i} y_{i} } & {\frac{{x_{i} }}{2}} & 0 & { - \frac{{x_{i} }}{2}} & {\frac{1}{2}} & {\frac{{x_{i} }}{2}} \\ \end{array} } \right].$$

(18)

Element matrices

The standard weak form for free vibration and buckling can, respectively, be expressed as:

$$\int\limits_{{S_{\text{e}} }} {\delta \{ \varepsilon \} }^{T} [D]\{ \varepsilon \} {\text{d}}S + \int\limits_{{S_{\text{e}} }} {\delta \{ U\} }^{T} [T]\{ \ddot{U}\} {\text{d}}S = 0,$$

(19)

$$\int\limits_{{S_{\text{e}} }} {\delta \{ \varepsilon \} }^{T} [D]\{ \varepsilon \} {\text{d}}S + \int\limits_{{S_{e} }} {\delta \{ \varepsilon^{\text{g}} \} }^{T} [\tau ]\{ \varepsilon^{\text{g}} \} {\text{d}}S = 0.$$

(20)

By substituting Eqs. (6), (7) and (13) into Eqs. (19) and (20), we obtain:

$$\delta \{ q_{\text{e}} \}^{T} \left( {\int\limits_{{S_{\text{e}} }} {[B]^{T} [D][B]} dS} \right)\{ q_{\text{e}} \} + \delta \{ q_{\text{e}} \}^{T} \left( {\int\limits_{{S_{\text{e}} }} {[N]^{T} [T][N]} {\text{d}}S} \right)\{ \ddot{q}_{\text{e}} \} = 0,$$

(21)

$$\delta \{ q_{\text{e}} \}^{T} \left( {\int\limits_{{S_{\text{e}} }} {[B]^{T} [D][B]} {\text{d}}S} \right)\{ q_{\text{e}} \} + \delta \{ q_{\text{e}} \}^{T} \left( {\int\limits_{{S_{\text{e}} }} {[B^{\text{g}} ]^{T} [\tau ][B^{\text{g}} ]} {\text{d}}S} \right)\{ q_{\text{e}} \} = 0.$$

(22)

Where the element stiffness, mass and geometrical stiffness matrices ([K^e], [M^e], $[K_{\text{g}}^{\text{e}} ]$), are, respectively, as:

$$\begin{aligned} & [K^{\text{e}} ] = \int_{{S_{\text{e}} }} {\left[ B \right]^{T} [D][B]} {\text{d}}S \\ & [K^{\text{e}} ] = [C]^{ - T} \underbrace {{\left( {\int {[Q]^{T} [D][Q]} \det (J){\text{d}}\xi {\text{d}}\eta } \right)}}_{{[K_{0} ]}}[C]^{ - 1} = [C]^{ - T} [K_{0} ][C]^{ - 1} , \\ \end{aligned}$$

(23)

$$\begin{aligned} \,[M^{\text{e}} ] = \int_{{S_{\text{e}} }} {[N]^{T} [T][N]} {\text{d}}S \\ & [M^{\text{e}} ] = \left[ C \right]^{ - T} \underbrace {{\left( {\int {[P]^{T} [T][P]} \det (J){\text{d}}\xi {\text{d}}\eta } \right)}}_{{[M_{0} ]}}[C]^{ - 1} = [C]^{ - T} [M_{0} ][C]^{ - 1} , \\ \end{aligned}$$

(24)

$$\begin{aligned} & \left[ {K_{\text{g}}^{\text{e}} } \right] = \int_{{S_{e} }} {\left[ {B^{\text{g}} } \right]^{T} [\tau ][B^{\text{g}} ]} {\text{d}}S \\ & \left[ {K_{\text{g}}^{\text{e}} } \right] = [C]^{ - T} \underbrace {{\left( {\int {[G]^{T} [\tau ][G]} \det (J){\text{d}}\xi {\text{d}}\eta } \right)}}_{{[K_{g0} ]}}[C]^{ - 1} = [C]^{ - T} [K_{g0} ][C]^{ - 1} . \\ \end{aligned}$$

(25)

The stress–strain relationship is given by:

$$\{ \sigma \} = [D]\{ \varepsilon \} ,$$

(26)

where $\{ \sigma \} = \{ M_{x} ,M_{y} ,M_{xy} ,T_{x} ,T_{y} \}^{T} ,$ $\{ \varepsilon \} = \{ \kappa_{x} ,\kappa_{y} ,\kappa_{xy} ,\gamma_{xz} ,\gamma_{yz} \}^{T} .$

where [D], [D]_b, [D]_s are, respectively, rigidity, bending rigidity, shear rigidity matrices and [T] is the matrix containing the mass material density:

$$\begin{aligned} & [D] = \left[ {\begin{array}{*{20}c} {[D]_{\text{b}} } & 0 \\ 0 & {[D]_{\text{s}} } \\ \end{array} } \right],\quad [D]_{\text{b}} = \frac{{Eh^{3} }}{{12(1 - \upsilon^{2} )}}\left[ {\begin{array}{*{20}c} 1 & \upsilon & 0 \\ \upsilon & 1 & 0 \\ 0 & 0 & {\frac{(1 - \upsilon )}{2}} \\ \end{array} } \right] \\ & [D]_{\text{s}} = khG\left[ {\begin{array}{*{20}c} 1 & 0 \\ 0 & 1 \\ \end{array} } \right], \\ \end{aligned}$$

(27)

$$[T] = \rho \left[ {\begin{array}{*{20}c} h & 0 & 0 \\ 0 & {\frac{{h^{3} }}{12}} & 0 \\ 0 & 0 & {\frac{{h^{3} }}{12}} \\ \end{array} } \right],$$

(28)

$$[\sigma_{0} ] = \left[ {\begin{array}{*{20}c} {\sigma_{x}^{0} } & {\sigma_{xy}^{0} } \\ {\sigma_{xy}^{0} } & {\sigma_{y}^{0} } \\ \end{array} } \right],\quad [\tau ] = \left[ {\begin{array}{*{20}c} {h[\sigma_{0} ]} & 0 & 0 \\ 0 & {\frac{{h^{3} }}{12}[\sigma_{0} ]} & 0 \\ 0 & 0 & {\frac{{h^{3} }}{12}[\sigma_{0} ]} \\ \end{array} } \right],$$

(29)

where $\sigma_{x}^{0}$, $\sigma_{y}^{0}$ and $\sigma_{xy}^{0}$ are the in-plane stresses.

The matrices [K₀], [M₀] and [K_g0] given in Eqs. (23), (24) and (25) are numerically computed with exact Gauss and Hamer rule integration, respectively, for quadrilateral and triangular elements (SBQP and SBTP4). The element stiffness, mass and geometrical matrices ([K^e], [M^e] and $[K_{\text{g}}^{\text{e}} ]$) can then be obtained. These are assembled to obtain the structural stiffness, mass and geometrical matrices ([K], [M] and [K_g]).

For static analysis, we use

$$[K]\{ q\} = \{ F\} .$$

(30)

For free vibration, we use

$$([K] - \omega^{2} [M])\{ q\} = 0.$$

(31)

For the buckling analysis, we use

$$([K] - \lambda_{\text{cr}} [K_{\text{g}} ])\{ q\} = 0.$$

(32)

Numerical validation

To validate the accuracy and efficiency of the formulated quadrilateral and triangular elements (SBQP and SBTP4), several numerical examples have been investigated for static, free vibration and buckling analysis of isotropic plates where the patch test of rigid body modes and the mechanic patch test are first carried out. The obtained results of the SBQP and SBTP4 elements are compared with other numerical and analytical solutions available in the literature.

Patch test of rigid body modes

To verify that both SBQP and SBTP4 elements pass the patch test of rigid body modes, the eigenvalues of the stiffness matrix for a single element are computed for various shapes and different aspect ratio. The only three zero eigenvalues obtained correspond to the three rigid displacement modes for a plate.

Mechanic patch test

In this patch test, a rectangular plate of (L = 2a = 40) length and (2b = 20) width simply supported at the three corner 1,2 and 3 (W₁ = W₂ = W₃ = 0) is considered where the plate is modeled by several elements as shown in Fig. 2 (Batoz and Dhatt 1990) for various side–thickness L/h ratio (10,100 and 1000). The plate boundaries are subjected to solicitations that produce the state of constant moments (or stresses). For the case of M_n = 1 applied on all sides (Fig. 2), the obtained results are M_x = M_y = 1 everywhere in the plate (Table 1). Whereas for the case of M_ns = 1 applied on all sides (Fig. 2), the obtained results at any points of the plate are M_xy = 1 (Table 1). The results given in Table 1 confirm that both SBQP and SBTP4 elements fulfill the mechanic patch test.

Fig. 2

Quadrilateral and triangular meshes for the patch test (E = 1000, ν = 0.3)

Table 1 Results of mechanic patch test

Square plates

A classical benchmark is first studied of square plate bending problem (Fig. 3) with different boundary conditions and various thickness–side (h/L) ratios subjected to a uniform load (q = 1), where the shear locking free test and convergence investigation of central deflection are considered in this study.

Fig. 3

Square plate with a mesh of N × N elements (L = 10, E = 10.92, ν = 0.3, k = 5/6)

Shear locking free test is considered for a clamped square plate with several values of ratios (L/h = 10–1,000,000) using a mesh of 12 × 12. The central deflection results of the plate illustrated in Table 2 and Fig. 4, confirm that the new formulated elements (SBQP and SBTP4) are able to solve the shear locking problem when the plate thickness becomes gradually small. However, it is observed that the SBRP element (Belounar and Guenfoud 2005) exhibits from shear locking phenomena for (L/h > 100).

Table 2 Deflections at the center [(WD/qL⁴)100] of a clamped square plate with different aspect ratios

Fig. 4

Shear locking test (W_c/W_Ref) of a clamped square plate

Now, convergence tests of a square plate are investigated with three cases of boundary conditions [clamped, soft simply supported SS1 (W = 0), and hard simply supported SS2 (W = β_s = 0)]. Various values of h/L ratios of 0.1, 0.01, and 0.001 are considered for thick, thin and very thin plates, respectively. The obtained results of the vertical displacement at the center of the plate are presented in Tables 3, 4 and 5 and Figs. 5, 6 and 7, which show that:

Table 3 Central deflection [(WD/qL⁴)100] for clamped square plates with uniform load

Table 4 Central deflection [(WD/qL⁴)100] for SS1 square plates with a uniform load

Table 5 Central deflection [(WD/qL⁴)100] for SS2 square plates with a uniform load

Fig. 5

Central deflection [(WD/qL⁴)100] for clamped square plates

Fig. 6

Central deflection [(WD/qL⁴)100] for SS1 square plates

Fig. 7

Central deflection [(WD/qL⁴)100] for SS2 square plates

• Faster convergence towards analytical solutions (Taylor and Auricchio 1993) is obtained using only a small number of elements for all cases of ratios (h/L = 0.1, 0.01, and 0.001) and boundary conditions.

• The SBQP and SBTP4 elements have similar behaviors for thin and very thin plates (h/L = 0.01, 0.001); whereas, for thick plates (h/L = 0.1), the SBTP4 element is a little better than the SBQP element.

• Both proposed elements are free from shear locking phenomena where they are able to provide excellent results for thin and very thin plates (h/L = 0.01, 0.001).

• Slow convergence to analytical solutions (Taylor and Auricchio 1993) is obtained using the SBRP element (Belounar and Guenfoud 2005) for thick and thin plates (h/L = 0.1, 0.01) and suffers from shear locking for very thin plates (h/L = 0.001).

Skew plates

To show the performance of the present elements to the sensitivity of mesh distortion, two examples of thin skew plates subjected to a uniform load (q = 1) are considered which are known in the literature as severe tests and studied by many researchers (Razzaque 1973; Morley 1963). The first is concerned with Razzaque’s skew plate (Razzaque 1973) (β = 60°) with simply supported on two sides and free on the other sides (Fig. 8). The results of the vertical displacement at the center of the plate using uniform meshes N = 2, 4, 8, 12 and 16 are given in Table 6 and Fig. 9 for (h/L = 0.001). The obtained results for both elements (SBQP and SBTP4) are in quite a good agreement with the reference solution given by Razzaque (1973). But it can be seen that the SBTP4 element is a little better than the SBQP and MITC4 (Nguyen-Xuan et al. 2008) elements.

Fig. 8

Skew plates (a Razzaque, b Morley) with N × N meshes (L = 100, E = 10.92, ν = 0.3, k = 5/6)

Table 6 Convergence of central displacement (W_c) for the Razzaque’s skew plate

Fig. 9

Central displacement (W_c/W_Ref) for the Razzaque’s skew plate

The second example treated by Morley (β = 30°) (Morley 1963) is simply supported (W = 0) on all sides (Fig. 8). Using meshes of N = 4, 8, 16 and 32, the obtained vertical displacement at the center of the plate are presented in Table 7 and Fig. 10 for h/L = 0.01 and 0.001. It can be observed that for h/L = 0.01, the results of the SBTP4 and SBQP elements are in good agreement with the reference solution (Morley 1963); whereas, for h/L = 0.001, the SBTP4 element is more efficient than the SBQP and MITC4 (Chen and Cheung 2000) elements.

Table 7 Convergence of central displacement (W_c) for the Morley’s skew plate

Fig. 10

Central displacement (W_c/W_Ref) for the Morley’s skew plate

Free vibration of square plates

Convergence tests of the formulated quadrilateral and triangular elements are first undertaken for simply supported (W = β_s = 0) and clamped plates with two thickness–side ratios (h/L = 0.005 and 0.1) (Fig. 3). The results of the first six non-dimensional frequencies (λ = (ω²ρL⁴h/D)^1/4) using the SBQP and SBTP4 elements with four regular meshes (N = 4, 8, 16 and 22) are presented in Tables 8, 9, 10 and 11 and Figs. 11 and 12 together with the four-node mixed interpolation of tensorial component MITC4 (Nguyen-Thoi et al. 2012), the discrete shear gap triangle DSG3 (Nguyen-Thoi et al. 2012) and the edge-based smoothed discrete shear gap triangular ES-DSG (Nguyen-Thoi et al. 2012) elements. It can be demonstrated that:

Table 8 Six first nondimensional frequency parameters (λ) of a SSSS thin square plate (h/L = 0.005)

Table 9 Six first nondimensional frequency parameters (λ) of a SSSS thick square plate (h/L = 0.1)

Table 10 Six first nondimensional frequency parameters (λ) of a CCCC thin square plate (h/L = 0.005)

Table 11 Six first nondimensional frequency parameters (λ) of a CCCC thick square plate (h/L = 0.1)

Fig. 11

Six first frequencies of a simply supported square plate with a 4 × 4 mesh

Fig. 12

Six first frequencies of a clamped square plate with a 22 × 22 mesh

• Both elements (SBQP and SBTP4) agree well with analytical solutions (Abbassian et al. 1987) and other elements (MITC4, DSG3, and ES-DSG) (Nguyen-Thoi et al. 2012).

• Figures 11 and 12 show that the SBQP and SBTP4 elements produce more accurate results than those given by other elements (MITC4, DSG3, and ES-DSG) (Nguyen-Thoi et al. 2012) when few elements are employed (4 × 4 mesh).

Having verified the convergence rate of the formulated elements, thin square plates (h/L = 0.005) with five different kinds of boundary conditions (SSSF, SFSF, CCCF, CFCF, and CFSF) for a 22 × 22 mesh are considered. The results of the four non-dimensional frequencies (λ = ωL²(ρh/D)^1/2) are presented in Table 12 and the first four mode shapes of SSSF and CFCF plates are plotted in Figs. 13 and 14. For all cases of boundary condition, the following can be concluded:

Table 12 Four first nondimensional frequency parameters (λ) of a thin square plate (h/L = 0.005)

Fig. 13

First four mode shapes of SSSF square plate using the SBQP element

Fig. 14

First four mode shapes of CFCF square plate using the SBTP4 element

• The present results are very close to analytical solutions (Leissa 1969) and are more accurate than those of the MITC4, DSG3 and ES-DSG elements (Nguyen-Thoi et al. 2012).

• The two elements (SBQP and SBTP4) have similar behavior, are shear locking free and their accuracy is insensitive to boundary conditions.

Free vibration of parallelogram plates

A cantilever parallelogram plate of skew angle = 60° with two h/L ratios (0.001 and 0.2) is studied (Fig. 15) using 22 × 22 mesh. The computed six non-dimensional frequencies (λ = ωL²/π²(ρh/D)^1/2) and the mode shapes are illustrated in Table 13 and Fig. 16, respectively. These results are compared with other numerical (DSG3, ES-DSG3, and MITC4) (Nguyen-Thoi et al. 2012) and analytical solutions (Karunasena et al. 1996). It can be seen that the SBQP and SBTP4 elements have a good accuracy compared to exact solutions (Karunasena et al. 1996) and are good competitors to ES-DSG3 and MITC4 (Nguyen-Thoi et al. 2012) and better than DSG3 (Nguyen-Thoi et al. 2012).

Fig. 15

Cantilever skew plate with a mesh of N × N elements

Table 13 Frequency parameters (λ) of cantilever skew plates (CFFF)

Fig. 16

Mode shapes of a cantilever skew plate with h/L = 0.2

Buckling of square plates subjected to uniaxial compression

Square plates subjected to uniaxial compression (Fig. 17) with h/L of 0.01 is analyzed for both simply supported (SSSS) and clamped (CCCC). The buckling load factor is defined as K^h = λ_crL²/(π²/D). The results of the buckling load factor for the SBQP and SBTP4 elements using 4 × 4, 8 × 8, 12 × 12, 16 × 16 and 20 × 20 meshes are presented in Table 14 and Fig. 18. For all cases of boundary condition, the two elements (SBQP and SBTP4) have similar results and converge to analytical solutions (Timoshenko and Gere 1970). In addition, these elements have excellent accuracy compared to other elements (DSG3 and ES-DSG3) (Nguyen-Xuan et al. 2010a, b).

Fig. 17

Square plate subjected to axial compression

Table 14 Convergence of uniaxial buckling load factor (K^h) of square plates with (h/L = 0.01)

Fig. 18

Convergence of uniaxial buckling load factor (K^h/K_exact) of square plates with h/L = 0.01

The results of the buckling load factor (K^h) and the relative error using 20 × 20 mesh are presented in Table 15. Numerical results of the SBQP and SBTP4 elements are in good agreement with analytical solutions (Timoshenko and Gere 1970) and other numerical solutions (Nguyen-Xuan et al. 2010a, b; Tham and Szeto 1990; Vrcelj and Bradford 2008; Liew and Chen 2004).

Table 15 Uniaxial buckling load factor (K^h) of square plates with (h/L = 0.01)

Buckling of square plates subjected to biaxial compression

Square plate subjected to biaxial compression (Fig. 19) with three essential boundary conditions (SSSS, CCCC, SCSC) is considered for h/L = 0.01 using a mesh of 16 × 16. The buckling load factor results (K^h = λ_crL²/(π²/D)) of the proposed elements are presented in Table 16 with analytical (Timoshenko and Gere 1970) and other numerical solutions (Nguyen-Xuan et al. 2010a, b; Tham and Szeto 1990; Vrcelj and Bradford 2008). It can be seen that both elements (SBQP and SBTP4) provide results which agree well with analytical solutions (Timoshenko and Gere 1970) and other solutions (Nguyen-Xuan et al. 2010a, b; Tham and Szeto 1990; Vrcelj and Bradford 2008) for all cases of boundary condition.

Fig. 19

Square plate subjected to biaxial compression

Table 16 Biaxial buckling load factor (K^h) of square plates with (h/L = 0.01)

Conclusion

A simple and efficient quadrilateral and triangular strain-based finite elements have been presented for static, free vibration and buckling analyses of Reissner–Mindlin plates. The four-node strain-based triangular element SBTP4 has the three engineering external degrees of freedom at each of the three corner nodes and one mid-edge point, while the quadrilateral element SBQP has the same engineering degrees of freedom at each of the four corner nodes. These developed elements passed successfully both patch and benchmark tests for plate bending problems. Numerical results show that the SBQP and SBTP4 elements are shear locking free, stable and superior to the original strain-based rectangular plate element (SBRP) (Belounar and Guenfoud 2005) which suffers from shear locking when the plate thickness becomes progressively very thin and has less rate of convergence to analytical solutions for thick and thin plates. The obtained results using both strain-based elements (SBQP and SBTP4) show that a rapid convergence to analytical solutions can be achieved with relatively coarse meshes compared with other robust elements based on different methods. In perspective, these elements can be superposed with membrane robust elements to construct shell elements for the analysis of complex shell structures.