Developments of Mindlin-Reissner Plate Elements

Since 1960s, how to develop high-performance plate bending finite elements based on different plate theories has attracted a great deal of attention from finite element researchers, and numerousmodels have been successfully constructed. Among these elements, themost popularmodels are usually formulated by two theoretical bases: theKirchhoff plate theory and theMindlin-Reissener plate theory. Due to the advantages that only C0 continuity is required and the effect of transverse shear strain can be included, the latter one seems more rational and has obtained more attention. Through abundant works, different types of Mindlin-Reissener plate models emerged in many literatures and have been applied to solve various engineering problems. However, it also brings FEM users a puzzle of how to choose a “right” one. The main purpose of this paper is to present an overview of the development history of theMindlin-Reissner plate elements, exhibiting the state-of-art in this research field. At the end of the paper, a promisingmethod for developing “shape-free” plate elements is recommended.


Introduction
Plate structure is a kind of important engineering components and has been widely used in many industries, such as aerospace engineering, shipbuilding industry, and civil engineering.To date, the finite element method is usually considered as the most convenient and effective tool in analyses and designs of various plate structures .The history of plate bending elements can be traced back to the 1960s [1,2].Since then, how to develop high-performance models has attracted a great deal of attentions from finite element researchers .Among these elements, the most popular models are usually formulated by two theories: the Kirchhoff plate theory and the Mindlin-Reissner plate theory.Most early researches focused on the models following the classical Kirchhoff thin plate theory, which can be found in related review paper [5].However, because transverse shear deformation is neglected, the Kirchhoff plate theory is only appropriate for modeling thin structures.Furthermore, since  1 continuity between adjacent elements is required, the construction procedures of Kirchhoff plate elements are quite difficult and more complicated.Actually, construction of  1 approximation for arbitrary finite elements is still a very challenging problem nowadays.
On the other hand, in Mindlin-Reissner (thick) plate theory [6,7], the deflection  and the rotations   ,   are independently defined, so that the influences of the transverse shear strain can be included, and only  0 continuity is required when formulating the interpolation functions of kinematics variables.Hence, the construction of thick plate elements is indeed easier than that of the thin plate elements [8], and more and more attention was paid to develop plate bending elements based on the Mindlin-Reissner plate theory .However, the low-order standard isoparametric displacement-based plate elements without special treatments only produce poor results in the thin plate case due to the false shear strains which lead to the shear locking phenomenon.This shear locking problem is primarily caused by the inability to model the Kirchhoff-type constraints in thin plate limits.On the discrete level, only if the trial function space is rich enough, the element can possess good properties.Unfortunately, this is not an easy case to standard low-order elements [9].To overcome this difficulty, different methods have been proposed, such as the reduced integration 2 Mathematical Problems in Engineering [10,11] and the selected reduced integration [12,13].In addition to the shear locking problem, there are two other great numerical challenges in Mindlin-Reissner plate theory, that is, the singularity caused by the plate obtuse corners and the edge effect caused by certain boundary conditions [14,15].Although these two do not pose great effects on the entire structure, they do complicate the numerical analysis.But, no enough attention has been paid to them.
Through substantial attempts from scholars, various Mindlin-Reissner plate bending elements  are available in numerous literatures.Furthermore, by using the higherorder elasticity theories, some models based on the nonclassical Mindlin-Reissner plate theories are proposed [189][190][191], such as the modified couple stress theory [192,193].However, no one has been treated as the so-called "best" one.Thus, how to choose a "right" one is quite puzzling to the FEM users.At the present stage, the major drawback in Mindlin-Reissner plate elements is that the element's performance is highly sensitive to mesh distortion.How to develop a "shapefree" element is still an open topic.The expression "shapefree" means that the element can always perform well regardless of the mesh quality.In this paper, a review on the development of the Mindlin-Reissner plate elements is firstly proposed.Then, a promising "shape-free" Mindlin-Reissner plate element method is briefly recommended.

Developments of Mindlin-Reissner Plate Elements
The first attempt to employ the Mindlin-Reissner plate theory in finite element method was the 8-node isoparametric element derived from the degenerated shell approach [17].In early time, most Mindlin-Reissner plate elements were the displacement-based models.A common problem encountered was that these elements often gave poor results in thin plate limit if their element stiffness matrices were evaluated using full integration schemes.This phenomenon is the wellknown "shear locking" problem.The first explanation on the shear locking was presented by Zienkiewicz and Hinton [18].They found that, with the use of a penalty function formulation, the product of the shear stiffness matrix and the displacement vector should be zero when the thickness-span ratio was sufficiently small.This discovery led to the development of so-called "locking indicator" [11], which measured the singularity of an element stiffness matrix.Also based on this work [18], Hughes et al. [12,19] gave an alternative concept, "constraint index, " to predict locking behavior in different finite element meshes.However, this index often deduced an overly pessimistic assessment.Afterwards, Spilker and Munir [20,21] proposed a modified version, called "rotational constraint index, " which can achieve better correlations.In order to gain more complete understanding and better estimation of the shear locking, Averill and Reddy [23] proposed another analytical technique verifying the roles that shear constraints play in thin plate.These constraints were expressed in terms of the nodal DOFs and then were judged to be the proper Kirchhoff constraints or spurious locking constraints.
In 1971, in order to eliminate the shear locking problem existing in the Mindlin-Reissner plate elements, Zienkiewicz et al. [10] introduced the famous reduced integration scheme which can drastically improve the element performance by reducing the integration order for all stiffness matrix components.This work [10] was based on the element given by [17].Almost at the same time, an analogous work was proposed by Pawsey and Clough [24].Although the reduced integration method did not have theoretical explanations, it indeed was a significant improvement for plate bending elements [25].Another effective approach was the selective reduced integration, whose validity had been subsequently demonstrated by Hughes and other scholars [12,13,19,26].Its main idea was to split the strain energy into two parts: bending part and shear part.Lower order integration was used for the troublesome shear strain energy, while full integration was employed for the bending energy.The applications of selective reduced integration may need some modifications in code procedures [27], so that the computational cost will be increased.Since then, these two reduced integration schemes have been widely applied, helping eliminate or alleviate the shear locking phenomenon.Prathap and Somashekar [28] proposed a vagarious strategy with (2 × 1) and (1 × 2) Gauss integration, respectively, for the terms   and   .The derived element QUAD4-CC was free of locking and had a proper rank in rectangular forms.Unfortunately, this advantage deteriorated when the element was not rectangular.
Although improved results can been obtained by using reduced or selective reduced integration in most cases, these methods cannot always ensure the absence of shear locking.Belytschko et al. [31] proposed a so-called mode decomposition approach combined with one-point quadrature, which can show improved performances compared to a straightforward application of reduced integration.However, it was still found to shear lock with certain mesh patterns.Furthermore, there were even some cases that elements cannot be improved by them at all [29].Moreover, the most noteworthy and fatal problem was that reduced or selective reduced integration often brought about spurious zero energy modes [32], which will yield unreliable results under certain boundary conditions.An exception that possessed a correct rank and did not lock was the so-called "heterosis" element [26] which utilized an 8-node interpolation for rotations and 9-node interpolation for deflections.In order to eliminate the rank deficiency, the stabilization method was proposed by Belytschko et al. [33][34][35][36].The essential ingredient of this method was to superimpose a "stabilizing" stiffness matrix to the stiffness matrix obtained by reduced integration.The most well-known one was the -method initiated by Flanagan and Belytschko [37], in which a vector  orthogonal to displacement fields was introduced.From the view of computational cost, this procedure seemed attractive and economical for constructing simple elements.However, additional researches were still needed for understanding the physical meanings of the stabilizing parameters more clearly [27].
Another alternative approach to alleviate shear locking, named by "assumed natural strain (ANS)" method, was proposed by Hughes and Tezduyar [38] and MacNeal [39].In [38], Hughes and Tezduyar were the first to find that shear locking can be avoided without losing important energy modes.The main idea of ANS method was to define the troublesome shear strains independent of the approximation of displacements.Then, these independent assumed strains were connected with the nodal displacements through the strain-displacement relation calculated at certain discrete points.Finally, they were interpolated over the element by using specific shape functions.This method had been proved to be mathematically valid by Brezzi et al. [40], and the key to the success was related to the choices of these preselected points.After then, based on this concept, many models were successfully developed [41][42][43][44][45].For instance, Tessler and Hughes [42,43] proposed the quadrilateral and triangular plate elements, named as MIN4 and MIN3, in which "anisoparametric" representations were used together with continuous transverse shear edge constraints.In their work, a concept of "finite-element-appropriate" shear correction factor, derived through multiplying the classical shear factor with a factor relating to the element's geometrical and material properties, was also proposed.Sze et al. [44] developed a triangular element AST6, which had a favorable constraint index of shear locking and optimized strains with respect to a linear pure bending field.Afterwards, they gave an improved model [45] in where the matrix inversion in original formulations can be avoided.Some variations of the ANS method were then presented, including the well-known mixed interpolated tensorial components (MITC) method proposed by Bathe et al. [40,[46][47][48][49]. Their method consisted in interpolating the covariant shear strain components   and   in a particular manner.A series of MITC plate/shell elements were constructed for linear/nonlinear analysis [50][51][52][53][54][55][56], among which the 4-node element MITC4 [47] may be the most widely used one in commercial codes [51], although it experienced some slight locking and shear oscillation in distorted meshes.Based on MITC4, Kebari [57] proposed the one point integrated version, in which the Taylor series expansion approach [52] was employed to accommodate the linear variation of stress within the element.
The discrete shear gap (DSG) method [58] was another scheme analogous to the ANS approach.The difference between these two methods was that the lack of collocation points makes the application of the DSG method independent of the element's order [59].Thus, the DSG method was employed by many elements to eliminate the shear locking [60][61][62][63].
In 1993, Zeinkiewicz et al. proposed the linked interpolation method [75,76], which has also been applied by other researchers later.The main thought of this concept was to improve the approximated deflection fields by the rotational DOFs.In this method, one order higher polynomials can be used for the deflection than rotations without adding additional element parameters.Thus, it can help meet the requirements in thin plate limit that the rotations were equal to derivatives of the deflections.It is known that the linked interpolation on its own cannot ensure free of shear locking, unless some other techniques were simultaneously employed, such as bubble modes, selective/reduced integration, or others.Even so, many elements were still developed based on the linked interpolation techniques [25,[75][76][77][78][79][80][81][82][83][84].Besides, some researchers started to study the relation between the linked interpolation elements and other models.For example, T3BL was proved to be essentially the same as the MIN3 [87].Recently, inspired by the higher-order linked interpolation functions developed for the Timoshenko beam, Ribarić and Jelenić [80] found a distinct basis for developing plate elements.Later, they [81] proposed two higher-order 9-node quadrilateral plate finite elements in which the interpolation functions were expended with three bubble parameters for keeping the polynomial completeness.Although the higherorder linked interpolation can make elements perform well in distorted mesh, it indeed led to vast computational cost owing to additional bubble degrees of freedom.
The hybrid/mixed approach was another way which can effectively and thoroughly eliminate the shear locking problem (see [20,21,32,[75][76][77] and the references herein).Lee and Pian [88] explained why shear locking appears in the thin limit case and gave some solutions in hybrid/mixed models.To develop hybrid/mixed elements, various variational principles were available, including the principle of minimum complementary energy, the Hellinger-Reissner variational principle, the Hu-Washizu variational principle, and their modified forms.In the hybrid/mixed method, displacements or boundary displacement and stresses/strains within the element were independently assumed.Thus, performances of these hybrid/mixed elements were largely based on the appropriate choice of parameters.However, the lack of specific rules in a priori choosing the "best" stress combination may be a major drawback.A significant work was proposed by Malkus and Hughes [13].They showed that there was equivalence between mixed models and displacement models with reduced integration.After this work, a series of effective hybrid/mixed elements were proposed, including some models mentioned above.Spilker and Munir [20] developed a high-order hybrid-stress element with use of the 12node cubic serendipity shape functions and 36 stress parameters.This element had a favorable "rotational constraint index" [21].Saleeb and Chang [27] proposed a HMPL5 element with three bubble functions represented by a fifth internal node.However, the use of internal node brought inconvenience and seemed unpopular.Simo and Rifai [117] presented a three-field mixed formulation in terms of displacements, stresses, and an enhanced strain field.Since the stress interpolation was assumed L2-orthogonal to the enhanced strain interpolation, the stress field can be eliminated from the final formulations.Thus, it collapsed into a two-field mixed method.This method can be used to develop low-order elements with enhanced accuracy in coarse meshes.Ayad and his cooperators proposed the mixed shear projected (MiSP) elements, including the stand bilinear models [99] and the high-order ones [113,114].This approach represented the transversal shear forces and shear strains in a particular manner: the shear forces were derived from the bending moments using the equilibrium equations and the shear strains were defined in terms of the edge tangential strains.Brasile [25] constructed element TIP3, whose stress fields were assumed from a moment-shear uncoupled polynomial approximation and then forced to satisfy some equilibrium conditions so that the minimum number of the stress parameters can be obtained.de Miranda and Ubertini [100] also proposed a self-equilibrium hybrid element.Particularly, its stresses were expressed by uncoupled polynomial expansions in terms of skew coordinates [119].Pereira and Freitas [115,116] employed Legendre polynomials to express the resultants within the element domain and the deflections along boundaries.Because of the orthogonality of the Legendre approximation, a highly sparse finite element solving system was derived.By combination of one-point integration and a stabilization matrix, Gruttmann and Wagner [111] proposed a model whose stiffness matrix can be analytically integrated, so that some computational costs can be saved.Even though there were some inherent drawbacks existing in hybrid/mixed method, such as the need of matrices inversion, it still inspired researchers to develop new types of novel hybrid/mixed models.
Recently, some scholars, mostly among mathematicians, began to employ the discontinuous Galerkin (DG) finite element methods [17] to design plate elements [120][121][122][123][124][125][126].This DG method admitted the discontinuities in the element discrete space, leading to new types of conforming or nonconforming elements.Brezzi and Marini [121] used nonconforming piecewise linear functions for both rotations and transversal displacements and added internal degrees of freedom to make the element perform better.Later, Lovadina [123] proposed a simplified version by using piecewise constant shear stress approximation functions.Arnold et al. [120] developed two DG Mindlin-Reissner plate element families: one was the fully discontinuous case and the other was the case of continuous rotations and nonconforming deflections.Hansbo et al. [126,127] also proposed a quadrilateral element model with continuous displacements and discontinuous rotations.In this model, the rotations had the same order as the parametric derivatives of the displacements in the parametric reference plane, and they obeyed the same transformation law.Boesing et al. [125] developed an interior penalty DG model.In their work, the a priori error analysis of DG methods was also presented.Although the discontinuous Galerkin method seems somewhat unfamiliar to traditional FEM, it was clear that the DG technique can offer a different way to these problems difficult for a classical approach [120].
Unlike that most quadrilateral elements were constructed in Cartesian or isoparametric coordinates, Cen et al. [128] extended the quadrilateral area coordinate method [129,130] to formulate a new element AC-MQ4 by employing the general conforming theory [148,149].A distinct feature of this method was that the transformation relations between the area coordinates and Cartesian coordinates were always linear, so that the order of the displacement fields expressed in form of area coordinates will not deteriorate in distorted meshes.Hence, the resulting element AC-MQ4 was quite insensitive to mesh distortion.When constructing the element, the locking-free Timoshenko's beam formulas were employed to help eliminate the shear-locking problem.Recently, this technique has also been widely used in developing Mindlin-Reissner plate elements [131][132][133][134][135][136].
The smoothed finite element method (S-FEM), proposed by Liu et al. [137,138], was also extended to Mindlin-Reissner plate elements.This method utilized the strain smoothing technique of mesh-free methods [139] to improve the accuracy of standard FEM.When developing smoothed Mindlin-Reissner plate models, the locking-free curvature smoothing formulation is mainly used.In S-FEM models, smoothing domains can be located inside the elements or cover parts of adjacent elements.Thus, based on different smoothing procedures, different types of S-FEM elements were derived, including CS-FEM [51,140], NS-FEM [61], and ES-FEM [62,[141][142][143].Nguyen-Thoi et al. [140] applied the cellbased strain smoothing technique to modify the well-known element MIN3 [87] to formulate a new element CS-MIN3.In CS-FEM, the number of supporting nodes was the same as those in elements.Therefore, the bandwidth of stiffness matrix in the CS-FEM was similar to the standard FEM.Nguyen-Xuan et al. [61] proposed a node-based smoothed model (NS-FEM), in which the strain smoothing was associated with the nodes of the elements.This NS-FEM can provide an upper bound solution in energy norm for the elasticity problem [144].By introducing the curvatures smoothing technique into the well-known MITC4 [47] and incorporating stabilized conforming nodal integration, Nguyen-Xuan et al. [141] proposed the MISC models.These elements can provide more accurate results, especially in distorted meshes.Furthermore, also based on the concept of S-FEM, Liu et al. [60,[145][146][147] proposed the alpha finite element method (FEM), in which the element was enhanced by the additional averaged nodal strain terms with an adjustable parameter , resulting in a properly softer stiffness formulation.
Based on the mechanism of shear locking phenomenon and the potential functional of Mindlin-Reissner plate, Cai et al. [150] extended the generalized mixed variational principle (GMVP) [151], which was a type of parameterized variational principles, to Mindlin-Reissner plate.In their models, splitting factors were introduced to adjust the shear potential and complementary energy components.Through this way, the element stiffness matrix can be artificially changed, which can make the element free of shear locking and possess higher accuracy.
The so-called free formulation method, introduced in 1975 by Bergan and Hanssen [152], was a radically different approach.Instead of various variational principles, this method [152][153][154][155][156][157] was directly derived from the conditions that the elements satisfy the requirements of the patch test and the rigid body motions.Such requirements can be greatly relaxed by slight modification of the coupling stiffness between fundamental and higher-order displacement modes [152].Thus, the shear deformations were only associated with the higher modes and can be tackled in a special way [153,154].
In addition to the methods mentioned above, there are still many other models and methods proposed by numerous literatures .But among them only a fraction of typical models can be briefly introduced here.Through an adaptive mesh-refinement method, Holzer et al. [160] proposed a model employing the enough higher-order interpolation for shape functions.Similar to the model given by [151], Arnold and Brezzi [177] also developed a model containing a shear energy splitting parameter.What is different was that this parameter was used to determine how much shear energy will be exactly integrated.Suggestions for appropriately choosing such parameter were also presented by [178].Similarly, Wang et al. [171,172] applied the combined hybrid method to design quadrilateral elements, which can make the system's energy optimal by adjusting the combined parameter.Polit et al. [179] proposed an 8-node quadrilateral element, in which each monomial term of the interpolation functions for the normal rotations was matched by the derivatives of its counterpart deflection.Zhou [176] took the linear combination of the mixed and displacement-based formulations to present a new variational formulation.Cheung and Chen extended the refined nonconforming element method [161] to develop Mindlin-Reissner plate elements [162,163].Duan and Liang [174] reformulated the classical mixed variational principle by relaxing the continuity of the transverse displacement while strengthening the normal continuity of the shear stresses.Pontaza and Reddy [175] proposed a formulation based on the least-squares variational principle, in which the highorder nodal expansions and the full integration were used.Maunder and de Almeida [170] proposed equilibrium models in which the stress fields were divided into hyperstatic part and spurious kinematic part.Ming and Shi [173] provided a unified variational formulation for three different elements, MIN4 [42], Q9 [183], and FMIN4 [184], and proved their optimal error bounds.Castellazzi and Krysl [164] proposed a new assumed-strain technique with use of the nodal integration, in which the weighted residual method was employed to weakly enforce the equilibrium equation and the kinematic equation.Later, they [167] verified the insensitivity of this method to geometrical distortion in a variety of distorted meshes.
Other than these conventional models mentioned above, there were some transition elements being proposed for the compatible connection of nonmatching meshes [185][186][187][188]. Sohn and Im [185] developed a quadrilateral model which can include an arbitrary number of additional nodes on each edge.Choi and Park [186,187] proposed different types of transition plate bending elements for adaptive ℎ-refinement.Sofuoglu and Gedikli [188] proposed a 5-node transition element by using the refined nonconforming method [161].
Generally, in order to develop a good element, more than one technique mentioned above can be employed in a combination.After decades of development, the Mindlin-Reissner plate element has been in a rather mature status.However, there are still two major deficiencies: (i) the performances of low-order elements will deteriorate as the meshes distort and (ii) the precision of stress is lower than that of displacement.The former is the more important one, especially in the nonlinear analysis.At present, this problem was generally alleviated at the expense of mathematical complexity or high computational cost and no effective and succinct solutions are available.For this problem, a promising solution was proposed recently, that is, the hybrid displacement function (HDF) method [194] for Mindlin-Reissner plate element.The main characteristic of this method is to assume the stress fields within element through the assumption of the displacement function  [195], rather than directly assuming themselves.All variables of Mindlin-Reissner plate theory can be expressed by such displacement function , and the derived resultants can satisfy all the governing equations.Furthermore, these solutions can degenerate into the ones of thin plate when the thickness is sufficiently small.The HDF elements are free of shear locking and can perform well in the severely distorted mesh, even when the mesh contains degenerated triangle or concave quadrilateral.Here, the outline of this method is introduced.
For a Mindlin-Reissner plate, the deflection  and rotations   ,   can be expressed by the displacement function  [195] as follows: in which  and  are the bending and shear stiffness and in (2) denotes the transversely distributed load.Then, by substituting (1) into the related governing equations of Mindlin-Reissner plate, the resultants R expressed as the functional of displacement function  are derived as follows: where R 0  ( = 1 ∼ ) are the general solutions of resultants and   are the corresponding coefficients; R * is the particular part.The fundamental analytical solutions of displacement function  and the corresponding solutions of the resultant forces have been listed in [194].Then, substitution of these solutions into the complementary energy functional of the Mindlin-Reissner plate element yields where Π  *  is the complementary energy within the element and   *  is the complementary energy along the boundary.
Finally, by applying the principle of minimum complementary energy, the element stiffness matrix K  and the element nodal equivalent load vector P   can be obtained as follows: in which q  is the nodal displacement vector.The elements derived from the HDF method can produce high-precision results for all field variables.An outstanding feature is that they are quite insensitive to mesh distortions, even when a severely distorted mesh containing concave quadrilateral or degenerated triangular elements is employed.Hence, this HDF method can be treated as a kind of "shapefree" finite element method for Mindlin-Reissner plate.In addition, the HDF method can also be applied to construct special elements for dealing with the edge effect problems of the Mindlin-Reissner plate [196], and the displacement function  can also be used for developing displacementbased element models [197].

Conclusions
Due to the fact that only  0 continuity is required and the effect of transverse shear strain can be included, the Mindlin-Reissner plate theory has become dominant in developing plate elements.In early history of Mindlin-Reissner plate elements, the major difficulty is the shear locking problem.
To overcome this problem, some effective and novel methods have been presented.Through plenty of attempts from scholars, different kinds of techniques and concepts are proposed, leading to different types of models.However, no one has emerged as the so-called "best" one.Thus, how to choose a right one is quite puzzling to the FEM users.So, the main purpose of this paper is to present a review on the development of the Mindlin-Reissner plate elements, briefly illustrating features of different models.
Furthermore, even though the development of Mindlin-Reissner plate elements has been in a quite good status, there are still some important difficulties that have not been solved yet, such as that the sensitivity problem to mesh distortion, which is more significant in the nonlinear analysis.Thus, how to develop a "shape-free" model ("shape-free" meaning the performance does not depend on the mesh's quality) is still an open topic.Therefore, at the end of this work, a promising "shape-free" Mindlin plate element model is recommended.