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

*Corresponding author https://doi.org/10.1590/1679-78256170 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 element


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 [1,2].
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 [3]. 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 [4].
According to Reissner-Mindlin plate theory, There are three independent variables: the lateral deflection of the plate middle plane and the two rotations of the normal line to the plate and . 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 [5][6][7].
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 [8,9]. P.G. Bergan and X. Wang [10] H. Abdalla and K. Hassan [11] 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 . 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 [10] used quadrilateral finite element, which is only complete for a polynomial of degree three.
Liu et al. [12] 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.

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 [13] are.
Where, and are the rotations of the transverse line about x and y-axes, respectively. Moreover, γ and γ 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. The stress vector due to bending and shear effects depending on stress-strain relationship in the isotropic homogenous elastic plate can be written as.
Where, [D] the flexure bending and shear rigidity of isotropic elastic plate is defined as.
In Eq. (5), k is the shear correction factor, which is taken as 5/6 [6], ν is the Poisson's ratio, E is the Young's modulus and is the modulus of rigidity.
By integrating the planar stresses in Eq. (3) over the thickness of the plate, bending moments and twisting moment per unit length in Fig.1 Moreover, the foundational equations of the shear force per unit length in Fig. 1 Moreover, according to the principle of virtual work [14], we get: The governing differential equations of plate can written as

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 γ x can written as [12]: A similar expression can be obtained for .
Bergan and Wang assumed that, the shear deformations of the plate and can be approximated by piecewise linear polynomials. This assumption is suitable for finite element technique and leads to the following expressions [10]: Where, Substituting Eq. (11.a and 11.b)))))))))), into the expressions for the bending and twisting moments (Eq. (6) and Eq. (7)) leads to: Finally, the strain energy expression per unit area, which is the sum of bending energy and shear energy according to Bergan and Wang approach [10,11], can be written as: Where,

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 [15,16] have used in the following plate bending problems and named as BWRE24 and BWRE36 respectively.

Derivation of element transverse displacement function
In both above plate elements, the transverse displacement function can expressed as generalized displacement modes The number of terms for polynomial [ ] must be equal to the total number of degrees of freedom that each element has. The transverse displacement coefficients vector { } (where, 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 { } T for each plate element BWRE24 and BWRE36. This leads to a square matrix [ ] of size × , consisting of numbers only [17]. Finally, the unknown coefficients vector { } can obtained from Therefore, the transverse displacement function can written as

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 [18]. Where, The generalized bending strain vector [ϵ b ] and shear strain vector [ϵ s ] can written according to the approach presented in Eq. (13,14)))))))))) and according to the finite element formulation of the transverse displacement function in Eq. (20) as Finally, the finite element formulation of the element stiffness matrix is expressed by: Similarly, according to the principle of virtual work (Eq. (8)), and the finite element formulation of the transverse displacement function (Eq. (20)), the finite element formulation of the element load vector [19] can written as After assembling the elements stiffness matrix and the element force vector to obtain the formation of the overall plate stiffness matrix [ ] and the resultant systems force vector { }, the final equilibrium equation can written as

The selection of the displacement field.
Consider a rectangular plate of dimensions ( × b) 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.

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 , two rotations and , the two degrees of freedom and which are proportional to the internal bending moments per unit length and , and finally which is proportional to the internal twisting moment per unit length . So, the generalized displacement field at each node ( ) of the element can expressed as: Because of the six degrees of freedom at each node, only twenty-four terms of Pascal's triangle are selected [20], Such that as ℎ 0 2 tends to zero BWRE24 will be as 24-dof displacement type finite element model [20]. Therefore, the deflection function as a complete polynomial of degree four can expressed as:

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 , two rotations and , and which is proportional to the internal twisting moment per unit length . So, the generalized displacement field at each node ( ) of the element can be expressed as: Because of the four degrees of freedom at each node, only thirty-six terms of Pascal's triangle have selected, Such that as ℎ 0 2 tends to zero, BWRE36 will be as the finite element model in [16]. Therefore, the deflection function as a complete polynomial of degree five can expressed as:  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 = 1 ,having ν =0.3, subjected to uniform load of intensity have been studied. 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 2 for case of two opposite edges simply supported and the other two edges free SFSF.

The results shown in
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 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 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

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 [24]. This test consists of a unitary square block [0, 1] 2 with clamped boundary conditions on all four sides and a distributed load ( . ) given by the following function: The exact solution for the displacement, rotations and bending moments have given in [24,25] as: 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.

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 [15] are good candidates, one needs to modify them using the same degrees of freedom as BWRE24.