A higher order triangular plate finite element using Airy functions

Abstract

The present paper describes the formulation of a new moderately thick plate bending triangular finite element based on Mindlin–Reissner plate theory. It is called a Great Triangular Moderately Thick Plate Finite Element, or GTMTPFE. The formulation is based on the strain approach, on solution of Airy’s function and on the analytical integration in the construction of the stiffness matrix. The strengths associated with this approach consist of:

• automatic verification of equilibrium conditions and kinematic compatibility conditions,

• the enrichment of the degrees of the interpolation polynomials of displacements, strains and constraints (refinement p),

• the consideration distortions sections related to Poisson effects,

• the treatment of blocking phenomena related to transverse shear.

In general, this approach results in a competitive, robust and efficient new moderately thick plate finite element. This is visible, on the one hand, through its stability against patch tests (constant twists, state of constants moments, transverse shear locking phenomenon, isotropy test). This is visible, through its good response to the patch tests to which it is subjected (constant torsions, state of constant moments, phenomenon of blocking in transverse shears, isotropy test). As has excellent convergence to the reference solution. Thus, it exhibits better performance behavior than other existing plate elements in the literature, particularly for moderately thick plates and for thin plates (L/h ratio greater than 4).

Introduction

In the civil engineering field, it is often necessary to calculate shell-type structures that can be relatively complex (caisson decks, slabs, dams, vaults, etc.). For such structures, it is very unusual for classical theories of plates and shells to give analytical solutions to the problem directly. Therefore, most of the time, the use of numerical methods such as finite element methods, is necessary. Beam finite elements often provide a first approximation of the behavior of many simple structures, but they prove to be inadequate for the study of more complex works. In these latter cases, it is necessary to use more efficient elements such as shell finite elements or three-dimensional finite elements. A brief non-exhaustive review of the literature allows us to identify with the aim of identifying a particular work for solving these problems.

Ahmad et al. ¹ presented a three-dimensional approach articulated around the curved iso-parametric elements, but their application shows a multitude of constraints ² :

• Need for a dense mesh, leading to numerical problems,

• Difficulties related to the study of thin shells, giving systems of poorly conditioned equations.

These difficulties was been avoided by the concept of the three-dimensional degenerate model introduced by Ahmad et al. This concept is based on the generalization of the formalism of the iso-parametric elements used for the shell elements. Unfortunately, this simplicity leads to the inherent difficulties with thin shells: membrane and transverse shear locking.

Many shell elements in the literature was been constructed by using the flat faceted approach. We can refer to elements DKT18 and DST18, where the membrane part has described by the T3 triangular element and the bending part by the DKT and DST elements. This approach presents a better convergence results and better representation of rigid modes, but accompanied by a number of essentially related disadvantages:

• The need for a dense mesh to ensure good convergence,

• The absence of the sixth degree of freedom (DDL) at the time of the assembly of coplanar elements, causing singularities at the level of the overall stiffness matrix. The avoidance of the singularity problem is supported, in some works, by the introduction at the membrane element formulation of a fictitious stiffness.

The approach of curved shell elements is generally described in the framework of the classical shell theory. In general, two models can be distinguished:

• The Koïter–Sanders model based on Love–Kirchhoff hypotheses, used by several authors^3–5;

• The Reissner–Naghdi model where the influence of transverse shear strains is taken into account.^6,7

From a practical point of view, these elements are arduous (a very heavy specification of data). They often require an analytical description of the shell, and either do not lend themselves or lend themselves inadequately to the calculation of the intersections of the shells.

In recent years, the shell element with planar facets (membrane element superimposed on the bending element) has become an important technique of research, due to the simplicity in its formulation. The main disadvantage of these elements is the blocking problem related to transverse shear when modeling plates taking into account the transverse shear. This problem is associated with a stiffening of the stiffness matrix of the element due to an overestimation of the matrix Kc. To solve this problem, many methods was been introduced. We illustrated them without limitation:

• the hybrid mixed variational method ⁸ using shear strain interpolation obtained by Timoshenko’s beam functions, ⁹

• the smoothed finite element method, ¹⁰ – the combined hybrid method, ¹¹

• the discrete shear gap method, ¹² – the enhanced assumed strain method, ¹³

• polygonal finite element method, ¹⁴

• hybrid displacement function element method and Mindlin–Reissner plate theory, ¹⁵ and

• the hybrid-Trefftez method used by Rezaiee-Pajand and Karkon ¹⁶ developing finite elements based on two independent displacement fields (internal and at the limits).

Furthermore, other methods was been used by a multitude of researchers based, particularly, on selective integration, Hourglass mode, substitution strain method and the strain approach used in finite element formulations. In addition to the blocking problem related to transverse shear, bending elements are sensitive to their geometry, especially when the mesh is irregular. Based on displacement model, the work of Wilson et al. ¹⁷ proposed an iso-parametric method by developing a hexahedral element at 8 nodes. Using the cell-based strain smoothing technique and discrete hear gap method, Nguyen-Thoi et al., ¹² presented a three-node triangular plate bending element. Although these elements show good performance, they remained sensitive to geometric distortion. To get around this problem, Cen et al., ¹⁸ presented “the hybrid stress-function element method which was based on the minimum complementary energy”. This approach resulted in an efficient element for very complex meshes. However, it is the stress field, which governs the equations in most of the Hybrid-Trefftz bending elements. As these elements, have zero modes parasites. As a solution, Cen et al., ¹⁹ presented “the hybrid displacement function element method”, which was formulated based on the extended hybrid-Trefftz stress element method. In this method, an assumed displacement function was been employed as variational functional of complementary energy and the displacements along each element edge were determined by the Timoshenko’s beam theory. All these works and others have confirmed the effectiveness of methods based on hybrid models where it is the stress field governing equations. As the problem of sensitivity of the elements to geometric distortion is not asserted in most models (moving, hybrid, mixed models). Based on the work of, Ayad R, ²⁰ other models based on the concept called the SFR principle (rotation of the spatial fiber) have been implemented by other researchers: Ghomari et al. ²¹ and Zouari et al. ²² Very recently Ayad et al. ²³ published a stabilized version of the SFR8 element using a technique inspired by the work. ²⁴ The strain approach was been used to formulate robust elements. The first elements was been developed on curved elements. ²⁵ Then this model is extended to the plane elasticity and to plate finite elements^26–29 and three-dimensional elastic elements.^30,31

In this paper, the objective of this contribution is therefore the formulation of a moderately thick plate finite element based on the strain approach, having as an objective the avoidance of these difficulties on the one hand, and the construction of a simple and efficient finite shell element for the analysis of complex structures on the other. To perform this step, we have enhanced our approach with concepts and development techniques based on:

• The adoption of the strain approach,

• The use of analytic integration for the evaluation of the stiffness matrix,

• The construction of polynomial interpolations of the different kinematic fields from the bi-harmonic solution of the Airy function.

• The consideration distortions sections related to Poisson effects.

This last concept aims to take the physical phenomenon related to the Poisson effect directly into account at the level of the approximation functions of the displacement interpolation fields.

Our first work ³² resulted in the construction of membrane triangular finite elements that can be easily combined with flexural finite elements (plates, beams). The present work constitutes a further development of this research, focused on the development of a bending triangular finite element based on Mindlin–Reissner plate theory. We call this GTMTPFE (Great Triangular Moderately Thick Plate Finite Element).

This element was been formulated using the strain approach. These benefits are summarized as follows:

• The automatic verification of kinematic compatibility conditions,

• The enrichment of the degrees of the polynomials of interpolation of displacements, strain, and constraints (refinement p); which allows you to quickly approach the exact solution during convergence tests,

• To raise the degrees of the polynomials of interpolation of the strain and displacements (degree 2 for the fields of the strain, degree 3 for the fields of rotations and degree 4 for the fields of displacements for the case of this element). This also makes it possible to circumvent the blocking problems linked to transverse shearing,

• The introduction of the bi-harmonic solution of the Airy function designed by Teodorescu, ³³ to introduce the Poisson’s ratio as variables directly in the fields of interpolation of the kinematics. This allows taking into account the distortion of the sections following the Poisson effects.

We use the principle of virtual displacements as a variational principle to formulate this element. The analytical integration for the evaluation of the stiffness matrix is highly interesting to avoid the loss of convergence, a phenomenon observed in the iso-parametric elements (when we use the numerical integration), which are very sensitive in the presence of distorted meshes. In order to validate our proposed new element GTMTPFE, we have undertaken successive tests. For each test, our results are compared with the corresponding reference solution and with the results given by different plate finite elements found in the literature. In general, this approach has resulted in a competitive, robust and efficient plate finite element GTMTPFE. This is noticeable through its stability in the face of patch tests (constant twists, state of constants moments, transverse shear locking phenomenon, isotropy test) on the one hand and in view of its excellent degree of convergence toward the reference solution on the other hand. It exhibits superior performance behavior to other existing triangular plate finite elements in the literature, particularly for moderately thick plates and for thin plates (L/h ratio greater than 4).

Formulation

The element characteristics

The element considered here is a flat plate bending triangular finite element of elasticity, called GTMTPFE (Figure 1). Each node has three degrees of freedom: the transverse displacement w_i and the rotations θ_xi and θ_yi. It was been formulated by using the strain approach.

Figure 1.

Characteristics of the element and the local-global transformation ${{\mathit{d}}^{\prime }}^{\left(\mathrm{e}\right)}={\mathit{T}}^{\left(\mathrm{e}\right)}{\mathit{d}}^{\left(\mathrm{e}\right)}$

First, the polynomial interpolation was been established to describe the strain fields. Subsequently, the interpolation functions of the displacement fields were been obtained by integration of these polynomials. For zero volume strengths, the choice of the approximation functions is done through the introduction of the Airy function φ (x, y), which makes it possible to reduce the problem of the equilibrium conditions:

$$\left[\begin{array}{ccc}\frac{\partial }{\partial x}\hfill & 0\hfill & \frac{\partial }{\partial y}\hfill \\ \multicolumn{1}{c}{0}& \frac{\partial }{\partial y}\hfill & \frac{\partial }{\partial x}\hfill \end{array}\right].\left\{\begin{array}{c}{\sigma }_x\hfill \\ \multicolumn{1}{c}{{\sigma }_y}\\ \multicolumn{1}{c}{{\tau }_{\mathrm{xy}}}\end{array}\right\}=\left\{\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\end{array}\right\}$$

(1)

to the simple bi-harmonic equation:

$${\nabla }^4\varphi (x,y)=0$$

(2)

With ${\sigma }_{\mathrm{x}}=\frac{{\partial }^2\varphi (\mathrm{x},\mathrm{y})}{\partial {\mathrm{y}}^2};{\sigma }_{\mathrm{y}}=\frac{{\partial }^2\varphi (\mathrm{x},\mathrm{y})}{\partial {\mathrm{x}}^2};{\tau }_{\mathrm{xy}}=-\frac{{\partial }^2\varphi (\mathrm{x},\mathrm{y})}{\partial \mathrm{x}\partial \mathrm{y}}$

For the plane stress state and isotropic materials, the relations (Hooke’s law) that relate the strain to the stress are as follows:

$$\begin{array}{c}\multicolumn{1}{c}{\left\{\begin{array}{c}{\sigma }_x\hfill \\ \multicolumn{1}{c}{{\sigma }_y}\\ \multicolumn{1}{c}{{\tau }_{\mathrm{xy}}}\end{array}\right\}=\frac{E}{1-{\upsilon }^2}\left[\begin{array}{ccc}1\hfill & \upsilon \hfill & 0\hfill \\ \multicolumn{1}{c}{\upsilon }& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \frac{1-\upsilon }{2}\hfill \end{array}\right].\left\{\begin{array}{c}{\epsilon }_x\hfill \\ \multicolumn{1}{c}{{\epsilon }_y}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{xy}}}\end{array}\right\};\left\{\begin{array}{c}{\epsilon }_x\hfill \\ \multicolumn{1}{c}{{\epsilon }_y}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{xy}}}\end{array}\right\}}\\ =\frac{1}{E}\left[\begin{array}{ccc}1\hfill & -\upsilon \hfill & 0\hfill \\ \multicolumn{1}{c}{-\upsilon }& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 2(1+\nu )\hfill \end{array}\right].\left\{\begin{array}{c}{\sigma }_x\hfill \\ \multicolumn{1}{c}{{\sigma }_y}\\ \multicolumn{1}{c}{{\tau }_{\mathrm{xy}}}\end{array}\right\}\hfill \end{array}$$

(3)

By introducing the Airy function φ (x, y), the stains can be expressed as follows:

$$\{\begin{array}{c}{\epsilon }_x=\frac{1}{E}(\frac{{\partial }^2\varphi (x,y)}{\partial y^2}-\upsilon \frac{{\partial }^2\varphi (x,y)}{\partial x^2})\hfill \\ {\epsilon }_y=\frac{1}{E}(-\upsilon \frac{{\partial }^2\varphi (x,y)}{\partial y^2}+\frac{{\partial }^2\varphi (x,y)}{\partial x^2})\hfill \\ {\gamma }_{\mathrm{xy}}=-\frac{2(1+\nu )}{E}\left(\frac{{\partial }^2\varphi (x,y)}{\partial x\partial y}\right)\hfill \end{array}$$

(4)

The bi-harmonic solutions of the Airy function φ (x, y) considered for the construction of the interpolation polynomial of the strain field does ³³ establish those. Each node of the considered element has three degrees of freedom. Thus the displacement fields, formulated by the use of the strain approach, have nine (9) dependent constants (a₁……., a₉). As a result, we consider the first nine solutions (see Appendix I):

• 3 terms leading to rigid body movements,

• 3 terms leading to linear displacements,

• 3 terms leading to high order displacements.

Kinematic

Considering the hypotheses of the Mindlin–Reissner plate theory, the displacements at a point of the solid with reference to a Cartesian axis system are as follows:

$$\{\begin{array}{c}u(x,y,z)=z{\beta }_x(x,y)=z{\theta }_y(x,y)\hfill \\ \multicolumn{1}{c}{v(x,y,z)=z{\beta }_y(x,y)=-z{\theta }_x(x,y)}\\ \multicolumn{1}{c}{w(x,y,z)=w(x,y)}\end{array}$$

(5)

The components of the infinitesimal strain tensor are then:

$$\{\begin{array}{c}{\epsilon }_x=\frac{\partial u(x,y,z)}{\partial x}=z\frac{\partial {\beta }_x(x,y)}{\partial x}\hfill \\ {\epsilon }_y=\frac{\partial v(x,y,z)}{\partial y}=z\frac{\partial {\beta }_y(x,y)}{\partial y}\hfill \\ {\gamma }_{\mathrm{xy}}=\frac{\partial u(x,y,z)}{\partial y}+\frac{\partial v(x,y,z)}{\partial x}=z(\frac{\partial {\beta }_x(x,y)}{\partial y}+\frac{\partial {\beta }_y(x,y)}{\partial x})\hfill \\ {\gamma }_{\mathrm{xz}}={\beta }_x(x,y)+\frac{\partial w(x,y)}{\partial x}\hfill \\ {\gamma }_{\mathrm{yz}}={\beta }_y(x,y)+\frac{\partial w(x,y)}{\partial y}\hfill \end{array}$$

(6)

The moments–curvatures relation are as follows:

$$\begin{array}{c}{K}_x=\frac{\partial {\beta }_x(x,y)}{\partial x};K_y=\frac{\partial {\beta }_y(x,y)}{\partial y};\hfill \\ K_{\mathrm{xy}}=(\frac{\partial {\beta }_x(x,y)}{\partial y}+\frac{\partial {\beta }_y(x,y)}{\partial x})\hfill \end{array}$$

(7)

Constitutive law

For the plane stress state and for isotropic materials, the generalized behavior law according to the Mindlin–Reissner plate theory is as follow:

$$\begin{array}{c}\left\{\begin{array}{c}{M}_x\hfill \\ \multicolumn{1}{c}{M_y}\\ \multicolumn{1}{c}{M_{\mathrm{xy}}}\\ \multicolumn{1}{c}{T_x}\\ \multicolumn{1}{c}{T_y}\end{array}\right\}=\frac{E{h}^3}{12(1-{\nu }^2)}\hfill \\ \left[\begin{array}{ccccc}1\hfill & \nu \hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\nu }& 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \frac{1-\nu }{2}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & \frac{6k}{h^2}(1-\nu )\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & \frac{6k}{h^2}(1-\nu )\hfill \end{array}\right]\left\{\begin{array}{c}{K}_x\hfill \\ \multicolumn{1}{c}{K_y}\\ \multicolumn{1}{c}{K_{\mathrm{xy}}}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{xz}}}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{yz}}}\end{array}\right\}\hfill \end{array}$$

(8)

With:

• E: Young’s modulus, h: the thickness of the plate, ν: Poisson’s ratio, k: the reduction coefficient of the section, generally taken equal to 5/6 for the rectangular sections.

• M_x, M_y, M_xy, T_x, T_y respectively represent the bending moments, the torsion moment and the shear forces for per unit of length.

Strains and displacements fields

The first nine bi-harmonic solutions of the Airy function φ (x, y) constitute the functions associated with the nine parameters of the interpolation polynomials of the stress fields, the strain fields and finally by integration of the displacement fields. The interpolation polynomials of the bending strain fields (detailed in Appendix I) are obtained by introducing these functions in equation (4). These are as follows:

$$\begin{array}{c}\left\{\begin{array}{c}{K}_x(x,y)\hfill \\ \multicolumn{1}{c}{K_y(x,y)}\\ \multicolumn{1}{c}{K_{\mathrm{xy}}(x,y)}\end{array}\right\}\hfill \\ =\frac{1}{\mathrm{zE}}\left[\begin{array}{ccccccccc}0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & -6\nu x\hfill & 6y\hfill & (2x^2-2\nu y^2)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 6x\hfill & -6\nu y\hfill & (-2\nu x^2+2y^2)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 2\hfill & 0\hfill & 0\hfill & -8\nu \mathrm{xy}\hfill \end{array}\right].\hfill \\ \left\{\begin{array}{c}{a}_1\hfill \\ \multicolumn{1}{c}{a_2}\\ \multicolumn{1}{c}{a_3}\\ \multicolumn{1}{c}{a_4}\\ \multicolumn{1}{c}{a_5}\\ \multicolumn{1}{c}{a_6}\\ \multicolumn{1}{c}{a_7}\\ \multicolumn{1}{c}{a_8}\\ \multicolumn{1}{c}{a_9}\end{array}\right\}\hfill \end{array}$$

(9)

These fields systematically check the equilibrium conditions (equation (1)) and the kinematic compatibility condition:

$$\frac{{\partial }^2{K}_x}{\partial y^2}+\frac{{\partial }^2{K}_y}{\partial x^2}=\frac{{\partial }^2{K}_{\mathrm{xy}}}{\partial x\partial y}$$

(10)

The interpolation polynomials of the shear strain fields (detailed in Appendix II) are obtained by satisfying the two other kinematic compatibility equations:

$$\{\begin{array}{c}\multicolumn{1}{c}{\frac{{\partial }^2{\gamma }_{\mathrm{xz}}}{\partial x\partial y}-\frac{{\partial }^2{\gamma }_{\mathrm{yz}}}{\partial x^2}+\frac{\partial K_{\mathrm{xy}}}{\partial x}=2\frac{\partial K_x}{\partial y}}\\ \multicolumn{1}{c}{\frac{{\partial }^2{\gamma }_{\mathrm{yz}}}{\partial x\partial y}-\frac{{\partial }^2{\gamma }_{\mathrm{xz}}}{\partial y^2}+\frac{\partial K_{\mathrm{xy}}}{\partial y}=2\frac{\partial K_y}{\partial x}}\end{array}$$

(11)

These are as follows:

$$\begin{array}{c}\left\{\begin{array}{c}{\gamma }_{\mathrm{xz}}(x,y)\hfill \\ \multicolumn{1}{c}{{\gamma }_{\mathrm{yz}}(x,y)}\end{array}\right\}=\frac{1}{z.E}\hfill \end{array}$$

$$\begin{array}{c}\left[\begin{array}{ccccccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & -6y^2\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & -6x^2\hfill & 0\hfill \end{array}\right].\left\{\begin{array}{c}{a}_1\hfill \\ \multicolumn{1}{c}{a_2}\\ \multicolumn{1}{c}{a_3}\\ \multicolumn{1}{c}{a_4}\\ \multicolumn{1}{c}{a_5}\\ \multicolumn{1}{c}{a_6}\\ \multicolumn{1}{c}{a_7}\\ \multicolumn{1}{c}{a_8}\\ \multicolumn{1}{c}{a_9}\end{array}\right\}\hfill \end{array}$$

(12)

By introducing the strain fields given by equations (9)–(12) into equation (6) and after integration, we obtain the following displacement fields:

$$\left\{\begin{array}{c}q(x,y)\hfill \end{array}\right\}=\frac{1}{z.E}\left[\mathbf{f}(\mathbf{x},\mathbf{y})\right].\left\{a_i\right\}$$

(13)

With: $\left\{\begin{array}{c}q(x,y)\hfill \end{array}\right\}=\left\{\begin{array}{c}w(x,y)\hfill \\ \multicolumn{1}{c}{{\beta }_x(x,y)}\\ \multicolumn{1}{c}{{\beta }_y(x,y)}\end{array}\right\}$

$$\begin{array}{c}\left[\mathbf{f}(\mathbf{x},\mathbf{y})\right]\hfill \\ =\left[\begin{array}{ccccccccc}1\hfill & -x\hfill & -y\hfill & (-\frac{x^2}{2})\hfill & (-\frac{y^2}{2})\hfill & -\mathrm{xy}\hfill & (-3x{y}^2+\nu x^3)\hfill & (-3y{x}^2+\nu y^3)\hfill & (-\frac{1}{6}(x^4+y^4)+\nu x^2{y}^2)\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & x\hfill & 0\hfill & y\hfill & (-3\nu x^2-3y^2)\hfill & -6\mathrm{xy}\hfill & (2/3x^3-2\nu x{y}^2)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill & y\hfill & x\hfill & 6\mathrm{xy}\hfill & (-3\nu y^2-3x^2)\hfill & (2/3y^3-2\nu y{x}^2)\hfill \end{array}\right]\hfill \end{array}$$

$${\left\{a_i\right\}}^T=\{a_1;a_2;a_3;a_4;a_5;a_6;a_7;a_8;a_9\}$$

Knowing the nodal coordinates (x_i, y_i) corresponding to the nodes i (i  =  1,…,3) and in application of the relation (13), the nodal displacements vector, at the elementary level, is given as follows:

$$\left\{\begin{array}{c}{q}^e\hfill \end{array}\right\}=\frac{1}{z.E}\left[\begin{array}{c}\left[\mathbf{f}({\mathbf{x}}_1,{\mathbf{y}}_1)\right]\hfill \\ \multicolumn{1}{c}{\left[\mathbf{f}({\mathbf{x}}_2,{\mathbf{y}}_2)\right]}\\ \multicolumn{1}{c}{\left[\mathbf{f}({\mathbf{x}}_3,{\mathbf{y}}_3)\right]}\end{array}\right]\left\{a_i\right\}$$

(14)

With $\left\{q^e\right\}$ the nodal displacements vector, such as:

$${\left\{q^e\right\}}^T=\{w_1;{\beta }_{x1};{\beta }_{y1};w_2;{\beta }_{x2};{\beta }_{y2};w_3;{\beta }_{x3};{\beta }_{y3}\}$$

[A]: the nodal coordinate’s matrix such as:

$$\left[\mathbf{A}\right]=\frac{1}{z.E}\left[\begin{array}{c}\left[\mathbf{f}({\mathbf{x}}_1,{\mathbf{y}}_1)\right]\hfill \\ \multicolumn{1}{c}{\left[\mathbf{f}({\mathbf{x}}_2,{\mathbf{y}}_2)\right]}\\ \multicolumn{1}{c}{\left[\mathbf{f}({\mathbf{x}}_3,{\mathbf{y}}_3)\right]}\end{array}\right]$$

(15)

From equations (14)–(15), we deduce the value of the parameters a_i :

$$\left\{a_i\right\}={\left[\mathbf{A}\right]}^{-1}\left\{\begin{array}{c}{q}^e\hfill \end{array}\right\}$$

(16)

These transformations make it possible to rewrite the kinematics equations as follows:

• Displacement fields

$$\left\{\begin{array}{c}q(x,y)\hfill \end{array}\right\}=\left\{\begin{array}{c}w(x,y)\hfill \\ \multicolumn{1}{c}{{\beta }_x(x,y)}\\ \multicolumn{1}{c}{{\beta }_y(x,y)}\end{array}\right\}=\left[\mathbf{f}(\mathbf{x},\mathbf{y})\right].{\left[\mathbf{A}\right]}^{-1}\left\{\begin{array}{c}{q}^e\hfill \end{array}\right\}$$

(17)

• Curvature fields related to the moments

$$\left\{\begin{array}{c}{K}_x(x,y)\hfill \\ \multicolumn{1}{c}{K_y(x,y)}\\ \multicolumn{1}{c}{K_{\mathrm{xy}}(x,y)}\end{array}\right\}=\left[{\mathbf{B}}_{\mathbf{f}}\right].{\left[\mathbf{A}\right]}^{-1}\left\{\begin{array}{c}{q}^e\hfill \end{array}\right\}$$

(18)

With: $\left[{\mathbf{B}}_{\mathbf{f}}\right]=\frac{1}{\mathrm{zE}}\left[\begin{array}{ccccccccc}0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & -6\nu x\hfill & 6y\hfill & (2x^2-2\nu y^2)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 6x\hfill & -6\nu y\hfill & (-2\nu x^2+2y^2)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 2\hfill & 0\hfill & 0\hfill & -8\nu \mathrm{xy}\hfill \end{array}\right]$

• Shear strain fields

$$\left\{\begin{array}{c}{\gamma }_{\mathrm{xz}}(x,y)\hfill \\ \multicolumn{1}{c}{{\gamma }_{\mathrm{yz}}(x,y)}\end{array}\right\}=\left[{\mathbf{B}}_{\mathbf{C}}\right].{\left[\mathbf{A}\right]}^{-1}\left\{\begin{array}{c}{q}^e\hfill \end{array}\right\}$$

(19)

With: $\left[{\mathbf{B}}_{\mathbf{C}}\right]=\frac{1}{\mathrm{zE}}\left[\begin{array}{ccccccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & -6y^2\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & -6x^2\hfill & 0\hfill \end{array}\right]$

Elementary stiffness matrix

The elementary stiffness matrix are obtained by application of the virtual displacements principle as follows:

• Stiffness related to flexural effects:

$$\left[{\mathbf{K}}_{\mathbf{f}}^{\mathbf{e}}\right]={\left[{\mathbf{A}}^{-1}\right]}^{\mathbf{T}}[\int \int {\left[{\mathbf{B}}_{\mathbf{f}}\right]}^{\mathbf{T}}\left[{\mathbf{D}}_{\mathbf{f}}\right]\left[{\mathbf{B}}_{\mathbf{f}}\right]\mathrm{dx}.\mathrm{dy}]\left[\left[{\mathbf{A}}^{-1}\right]\right]$$

(20)

• Stiffness related to shear effects:

$$\left[{\mathbf{K}}_{\mathbf{C}}^{\mathbf{e}}\right]={\left[{\mathbf{A}}^{-1}\right]}^{\mathbf{T}}[\int \int {\left[{\mathbf{B}}_{\mathbf{C}}\right]}^{\mathbf{T}}\left[{\mathbf{D}}_{\mathbf{C}}\right]\left[{\mathbf{B}}_{\mathbf{C}}\right]\mathrm{dx}.\mathrm{dy}]\left[\left[{\mathbf{A}}^{-1}\right]\right]$$

(21)

With: $\left[{\mathbf{D}}_{\mathbf{f}}\right]=\frac{E{h}^3}{12(1-{\nu }^2)}\left[\begin{array}{ccc}1\hfill & \nu \hfill & 0\hfill \\ \multicolumn{1}{c}{\nu }& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \frac{1-\nu }{2}\hfill \end{array}\right]\mathrm{and}\left[{\mathbf{D}}_{\mathbf{C}}\right]=\frac{\mathrm{Ekh}}{2(1+\nu )}\left[\begin{array}{cc}1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right]$

The evaluation of the expressions $[\int \int {\left[{\mathbf{B}}_{\mathbf{f}}\right]}^{\mathbf{T}}\left[{\mathbf{D}}_{\mathbf{f}}\right]\left[{\mathbf{B}}_{\mathbf{f}}\right]\mathrm{dx}.\mathrm{dy}]$ and $[\int \int {\left[{\mathbf{B}}_{\mathbf{C}}\right]}^{\mathbf{T}}\left[{\mathbf{D}}_{\mathbf{C}}\right]\left[{\mathbf{B}}_{\mathbf{C}}\right]\mathrm{dx}.\mathrm{dy}]$ is established by analytical integration of the different components resulting from the matrix products ${\left[{\mathbf{B}}_{\mathbf{f}}\right]}^{\mathbf{T}}\left[{\mathbf{D}}_{\mathbf{f}}\right]\left[{\mathbf{B}}_{\mathbf{f}}\right]$ and ${\left[{\mathbf{B}}_{\mathbf{C}}\right]}^{\mathbf{T}}\left[{\mathbf{D}}_{\mathbf{C}}\right]\left[{\mathbf{B}}_{\mathbf{C}}\right],$ whose expressions take the form $H_{(\alpha +1)(\beta +1)}=C.x^{\alpha }{y}^{\beta }$ .

Validation

Mechanical test patch

Several researchers to evaluate finite elements at sate of constants moments^34–37 have studied this patch test. A rectangular domain of sides 2a = 40 units and 2b = 20 units, is considered. We impose on the contour of this domain a loading reflecting the state of constants moments. In this example, the domain is discretized by four bending finite elements. In the first case, the moments on the contour are in agreement with the field defined by ${\overline{M}}_n=1$ (see Figure 2(a)) in order to have all along ${\overline{M}}_x={\overline{M}}_y=1$ (see Figure 3(a)). In the second case, the moments on the contour are in agreement with the field defined by ${\overline{M}}_{\mathrm{ns}}=1$ (see Figure 2(b)) in order to have all along ${\overline{M}}_{\mathrm{xy}}=1$ (see Figure 3(b)). We used thicknesses varying from 0.04 to 4.0 units to make these tests. The other data are defined as follows: E = 1000; ν = 0.3; h = (0.04; 0.4; 0.8; 2.0; and 4.0); W₁ = W₂ = W₃ = 0. The results for the moments inside the element (node 05) are given in Table 1. They show that the element can pass the test patch.

Figure 2.

Loading on the contour reflecting the constants moments state: (a) moment defined by ${\overline{M}}_n=1$ and (b) moment defined by ${\overline{M}}_{\mathrm{ns}}=1$ .

Figure 3.

Equivalents node’s forces: (a) first case and (b) second case.

Table 1.

Results of “mechanical” patch test of element Reissner-Mindlin type.

Loading on the contour	Thickness		h = 0.04	h = 0.4	h = 0.8	h = 2.0	h = 4.0
	Calculation results
${\overline{M}}_n=1$	${\overline{M}}_x={\overline{M}}_y=1$ everywhere	$M_{x\left(\mathrm{node}5\right)}$	1	1	1	1	1
		$M_{y\left(\mathrm{node}5\right)}$	1	1	1	1	1
		$M_{\mathrm{xy}\left(\mathrm{node}5\right)}$	0	0	0	0	0
${\overline{M}}_{\mathrm{ns}}=1$	${\overline{M}}_{\mathrm{xy}}=1$ everywhere	$M_{x\left(\mathrm{node}5\right)}$	0	0	0	0	0
		$M_{y\left(\mathrm{node}5\right)}$	0	0	0	0	0
		$M_{\mathrm{xy}\left(\mathrm{node}5\right)}$	1	1	1	1	1

Constant torsion

This test is used to evaluate the ability of finite elements by representing the constant torsion state. It is a square plate (Figure 4) of side L = 8.0, thickness h = 1.0, simply supported (W_B = W_C = W_D = 0) and subjected to a transverse force P = 5 applied to the corner A. The parameters of the material are Young’s modulus E = 10,000 and Poisson’s ratio ν = 0.30. The results obtained by GTMTPFE element are presented in Table 2. The comparison of these results with the exact solution and with those obtained by other elements extracted from the literature, shows a good performance of this element in the representation of constant torsion state.

Figure 4.

Torsion of a square plate.

Table 2.

Simulation results of the torsion of a square plate.

Type of element	Displacement “W”
	Point A	Point F
GTMTPFE (2 × 2)	0.24966	0.06
ACM (8 × 8)	0.24972	0.06
HCT (8 × 8)	0.25002	0.06
DKT	0.24960	0.06
Exact solution (Thin theory)	0.24960	0.06

Isotropy test

This test is used by Gholami et al., ³⁸ Kashtalyan, ³⁹ and Plevako ⁴⁰ in analysis of anisotropic functionally graded plates. It is a square plate simply supported on sides a = 0.9 m and thickness h = 0.3 m. The Young’s modulus varies according to the exponential distribution $E\left(z\right)=E_1{e}^{n(\frac{z}{h}-1)}$ and the Poisson’s ratio is constant υ = 0.3. The mechanical data are as follows: E₁ = 200GPa, uniform transverse loading equal to 5000 Mpa. The results obtained (see Table 3) and compared with those given by Gholami et al., ³⁸ Kashtalyan, ³⁹ and Plevako. ⁴⁰ Show a performance of the GTMTPFE element:

• Agreement with the general Plevako solution and the DKM method, with regard to the isotropic plate (homogeneous, n = 0),

• Rapid convergence of the GTMTPFE element: The result given by the functional gradation plate for the absolute value of the smallest Inhomogeneity parameter (n = 10⁻⁶) and in line with that given by the isotropic plate (n = 0).

Table 3.

Performance of GTMTPFE in comparison with current methods and DKM (GTMTPFE discretized with 200 meshes).

Inhomogeneity parameter (n)		Deflection at (a/2, a/2, h/2)
		W (a/2, a/2, h/2) (m)
		Current method	DKM 38	GTMTPFE
Functionally graded plates 39	−10−1	−0.024842415	−0.024842414	−0.024900715
	−10−2	−0.026430300	−0.026431240	−0.026308419
	−10−3	−0.026168025	−0.026166178	−0.026163792
	−10−4	−0.026178360	−0.026178513	−0.026178186
	−10−5	−0.026179725	−0.026179747	−0.026176746
	−10−6	−0.026179725	−0.026179871	−0.026176899
	10−6	−0.026179920	−0.026179898	−0.026176899
	10−5	−0.026179920	−0.026180022	−0.026177008
	10−4	−0.026181285	−0.026181255	−0.026178186
	10−3	−0.026193570	−0.026193598	−0.026189972
	10−2	−0.026317200	−0.026317326	−0.026308419
	10−1	−0.027585480	−0.027585469	−0.027554607
Homogeneous plate 40	0	−0.0261885	−0.0261799	−0.0261769

Convergence tests

Cantilever beam case

We use the present example to evaluate the behavior of the finite elements in the appearance of the locking phenomenon related to transverse shear. It is a cantilever console beam (Figure 5) subjected to simple bending, clamped at one end and subjected to a point force P = 0.1 N at the other end. The geometric and mechanical data used are length L = 10 m; width b = 1.0 m; Young’s modulus E = 1.2 × 10⁶ N/m² and Poisson’s ratio ν = 0. This test makes it possible to check the behavior of the GTMTPFE element in simple bending according to the slenderness (ratio L/h). The behavior of the finite element evaluated through the displacement W of the point A in the direction Oz is seen for several reports L/h, h being the thickness of the beam.

Figure 5.

Beam-Console subject to a point load.

Based on the Timoshenko beam theory, the exact solution of tip displacement for a plate is available as follows:

$$w_t=\frac{4P{L}^3}{\mathrm{Eb}{h}^3}(1+\frac{1}{2k}{\left(\frac{h}{L}\right)}^2)$$

(22)

With k reduction parameter for rectangular sections equal to 5/6,:

The results obtained by the GTMTPFE element, are given in Figure 6 and Table 4. They show a total absence of the locking phenomenon related to transverse shear, since the displacements obtained converge well toward the theoretical solution when the slenderness is important. The relative error is less than 0.05% for L/h ratio between 4 and 10 and equal to 0.00% for L/h ratio greater than 10, indicating its robustness compared to some existing elements in the literature, particularly for moderately thick plates and for thin plates (L/h ratio greater than 4). It is noting that its modeling has established with twenty (20) elements. We can conclude that:

• The GTMTPFE element converges to the analytical solution faster than the DKTM and DSTM elements for thin and moderately thick structures (L/h ratio greater than 4),

• The GTMTPFE element converges to the analytical solution faster than the ANST6 element for slender structures.

Figure 6.

Displacement W_A according to slenderness L/h.

Table 4.

Displacement W_A according to slenderness L/h.

L/h	DSTM	ANST6	DKTM	GTMTPFE	Theoretical solution
1	5.100E–07	5.300E–07	3.100E–07	3.330E–07	5.333E–07
2	2.900E–06	3.000E–06	2.500E–06	2.664E–06	3.067E–06
4	2.000E–05	2.200E–05	2.000E–05	2.131E–05	2.213E–05
5	4.000E–05	4.200E–05	3.900E–05	4.162E–05	4.267E–05
10	/	/	/	3.330E–04	3.353E–04
50	/	/	/	4.162E–02	4.168E–02
100	3.130E–01	3.230E–01	3.130E–01	3.330E–01	3.334E–01

N.B: The details of the elements using for comparison in this test are given in Annex III.

Cantilever plate

Rezaiee-Pajand and Yaghoobi ⁴¹ have analyzed the behavior of inflectional finite elements regarding their convergence use for test. As it has used to evaluate their response to the blocking problem linked to transverse shear encountered with certain finite elements based on the Mindlin–Reissner theory. It is a cantilever plate (Figure 7) clamped at one end and subjected to distributed shear force q = 4 at the other end. The geometric and mechanical data used are length L = 10 m; width b = 1.0 m; Young’s modulus E = 1.2 × 10⁶ MPa and Poisson’s ratio ν = 0. Based on the Timoshenko beam theory, the exact solution of tip displacement for a plate is available as follows:

$$w_t=\frac{4P{L}^3}{\mathrm{Eb}{h}^3}(1+\frac{1}{2k}{\left(\frac{h}{L}\right)}^2)===\to w_t=13,33413\mathrm{mm}$$

(23)

Figure 7.

Cantilever plate with distributed shear force at the end.

Two steps are considered in this test. First, the GTMTPFE element is subjected to a convergence test with a slender structure (L/h = 100). Then we use the optimal mesh obtained to evaluate its behavior according to a variation in its thickness: from the thickest (L/h = 5) to the thinnest (L/h = 1000). The results obtained by the GTMTPFE element, are given in Figures 8 and 9. They show a total absence of the blocking phenomenon linked to transverse shearing. The displacements obtained converge well to the theoretical solution and are in line with the results of the study realized by Rezaiee-Pajand. ⁴² Optimal discretization is obtained with 72 meshes. It is noting that, the GTMTPFE element is appropriate for any type of structures (thin and moderately thick). The relative error is less than 0.07% for moderately thick plates (L/h between 5 and 20) and less than 0.04% for thin plates (L/h greater than 20).

Figure 8.

Convergence curve for a cantilever plate with load distributed at its end.

Figure 9.

Relative error (%) on the displacement of the free side of the structure according to the variation of the thickness.

Isotropic square plate case

Many authors in literature have studied this example. It is used to check the capacity of finite elements, with respect to their convergence according to the size of the mesh. For this reason, we consider an isotropic square plate with a side a and a thickness h (Figure 10).

Figure 10.

Isotropic square plate.

This verification covers a number of different situations depending on the boundary conditions, the type of loading and the slenderness of the plate (a/h ratios). The reference solutions are as follows:

• Plate loaded by uniformly distributed load q:

$$w_{\mathrm{ref}.}={\alpha }_{q}q{a}^4/100D$$

(24)

• Plate loaded by a force P applied at the center (point C):

$$w_{\mathrm{ref}.}={\alpha }_{P}P{a}^2/100D$$

(25)

With: $D=\frac{E{h}^3}{12(1-{\nu }^2)}.$

According to the boundary conditions and slenderness of the plate, the values of α_q and α_P are given in Table 5.

Table 5.

Values of ${\alpha }_q$ et ${\alpha }_P$ .

Boundary conditions	Plate simply supported on its four sides		Plate recessed onits four sides
Slenderness	a/h = 10	a/h = 100	a/h = 10	a/h = 100
${\alpha }_q$	0.46169	0.4062	0.15046	0.1260
${\alpha }_p$	1.44267	1.1600	0.77775	0.5600

The results obtained for the deflection at the center of the plate (point C) are compared with reference solutions and those of other existing elements in the literature. They are presented in Tables 6 to 10 and illustrated by the Figures 11 to 18. The analysis of these results shows that:

• The GTMTPFE element converges well to the reference solution for all simulated situations (boundary conditions, loading type, and slenderness).

• There is an excellent performance in approaching the solution. Table 6 show the relative errors observed for a mesh size N = 8.

• In view of these observations and comparisons of the results with those of other existing elements in the literature, we notice the robustness of the GTMTPFE element in the treatment of moderately thick and thin structures.

Table 6.

Relative error $100(1-w_{\mathrm{cal}.}/w_{\mathrm{ref}.})$ .

Slenderness	Point load P in the center		Uniformly distributed load q
	Plate simply supported on all sides	Recessed plate on all sides	Plate simplysupported on all sides	Recessed plate on all sides
a/h = 10	−0.90%	−0.19%	2.14%	−5.31%
a/h = 100	−0.35%	0.48%	1.20%	0.82%

Table 7.

Values of α for simply supported plate loaded by uniformly distributed loading q.

a/h = 10–αq = 0.46169
N	R4	SBH8	GTMTPFE
2	0.23169	0.35935	0.31434
4	0.36519	0.43161	0.43584
8	0.43142	0.45299	0.47161
a/h = 100–αq = 0.4062
N	R4	SBH8	GTMTPFE
2	0.00446	0.08817	0.1932
4	0.01727	0.31235	0.3090
8	0.06128	0.39628	0.4111

Table 8.

Values of α for recessed plate on all sides, loaded by uniformly distributed load q.

a/h = 10–αq = 0.15046
N	R4	SBH8	GTMTPFE
2	0.06989	0.09089	0.12186
4	0.11518	0.13871	0.13092
8	0.13954	0.14789	0.14247
a/h = 100–αq = 0.1260
N	R4	SBH8	GTMTPFE
2	0.00101	0.00273	0.05049
4	0.00367	0.05624	0.09875
8	0.01322	0.11797	0.12704

Table 9.

Values of α for simply supported plate loaded by point loading P applied in center.

a/h = 10–αP = 1.44267
N	R4	SBH8	GTMTPFE
2	0.7358	1.0634	0.9426
4	1.1295	1.3014	1.2157
8	1.3461	1.4020	1.4297
a/h = 100–αP = 1.160
N	R4	SBH8	GTMTPFE
2	0.0134	0.2557	0.5825
4	0.0485	0.8563	1.1086
8	0.1690	1.1167	1.1559

Table 10.

Values of α for recessed plate loaded by a point loading P applied in the center.

a/h = 10–αP = 0.77775
N	R4	SBH8	GTMTPFE
2	0.32949	0.39486	0.38548
4	0.57104	0.65482	0.56384
8	0.71499	0.74492	0.77629
a/h = 100–αP = 0.56
N	R4	SBH8	GTMTPFE
2	0.00458	0.01078	0.1597
4	0.0169	0.24204	0.3603
8	0.06043	0.51216	0.5627

Figure 11.

Plate loaded by uniformly distributed loading q – simply supported, normalized displacement for a/h = 10.

Figure 12.

Plate loaded by uniformly distributed loading q – simply supported, normalized displacement for a/h = 100.

Figure 13.

Plate loaded by uniformly distributed loading q– recessed, normalized displacement for a/h = 10.

Figure 14.

Plate loaded by uniformly distributed loading q– recessed, normalized displacement for a/h = 100.

Figure 15.

Plate loaded by point loading P applied in center – simply supported, normalized displacement for a/h = 10.

Figure 16.

Plate loaded by point loading P applied in center – simply supported, normalized displacement for a/h = 100.

Figure 17.

Plate loaded by point loading P applied in center – recessed, normalized displacement for a/h = 10.

Figure 18.

Plate loaded by point loading P applied in center – recessed, normalized displacement for a/h = 100.

Circular plate

Several authors in literature, including Batoz and Dhatt, ³⁴ have studied this example. This study is performed according to several nuanced situations based on the boundary conditions of the plate under uniform loading f and slenderness R/h. The plate is discretized info three mesh concentrations: 16, 28, and 112 (Figure 19). The geometric and mechanical data are (R = 5, ν = 0.36, f = 0.1, E = 10⁷).

Figure 19.

Circular plate – mesh densities considered.

The reference solutions are:

$$W_{\mathrm{ref}}=\frac{f{R}^4}{64D}(\frac{5+\nu }{1+\nu }+Ø)\mathrm{for}\mathrm{simply}\mathrm{supported}\mathrm{plate}$$

(26)

$$W_{\mathrm{ref}}=\frac{f{R}^4}{64D}(1+Ø)\mathrm{for}\mathrm{simply}\mathrm{recessed}\mathrm{plate}$$

(27)

With: $Ø=\frac{16}{5}{\left(\frac{h}{R}\right)}^2\frac{1}{1-\nu }$ and $D=\frac{E{h}^3}{12(1-{\nu }^2)}$

The results are show in Figures 20 and 21 and the Table 11. Therefore, we have to conclude that the GTMTPFE element exhibits a good behavior with respect to this test and a quick convergence for the different simulated cases. This performance is more significant for the simply supported plate.

Figure 20.

Circular plate simply supported, relative error of the displacement in the center.

Figure 21.

Circular plate recessed, relative error of the displacement in the center.

Table 11.

Circular plate – displacement in the center.

h/R	Mesh	Wref	Wcalc	Wcal/Wref
Simply supported
0.04	16	5.035E–04	3.439E–04	0.68
	28	5.035E–04	4.284E–04	0.85
	112	5.035E–04	4.992E–04	0.99
0.1	16	3.257E–05	2.247E–05	0.69
	28	3.257E–05	2.880E–05	0.88
	112	3.257E–05	3.324E–05	1.02
Recessed
0.04	16	1.285E–04	3.474E–05	0.27
	28	1.285E–04	9.203E–05	0.72
	112	1.285E–04	1.055E–04	0.82
0.1	16	8.568E–06	2.676E–06	0.31
	28	8.568E–06	5.156E–06	0.60
	112	8.568E–06	7.578E–06	0.88

Skew plates

In this test, the behavior of two diamond skew plates was been investigated. The first case is that of a 60° bias plate subjected to a uniformly distributed load assumed to be simply supported on two sides and free on the other two sides. Its behavior was been studied according to three mesh densities: 4, 8, and 16 elements per side (Figure 22). Three forms of slenderness (L/h = 50, L/h = 100, and L/h = 1000) were considered. The geometric and mechanical data are (L = 100, ν = 0.3, f = 1.0, E = 10⁹).

Figure 22.

60° skew plate – mesh densities considered.

The second test, proposed by Morley, ⁴³ with 30° skew plate simply supported on all four sides. It presents a case of very severe behavior since the bending moment has a singularity at the obtuse angle. Due to its symmetry with respect to the diagonals of the plate, a quarter of the plate was been considered by confusing the Ox and Oy axes with the diagonals of the plate (Figure 23). In the same way as in the previous test, three forms of slenderness (L/h = 50, L/h = 100, and L/h = 1000) were studied. The geometric and mechanical data are (L = 100, ν = 0.3, f = 1.0, E = 10⁹).

Figure 23.

30° skew plate – mesh densities considered.

The reference solutions obtained by finite differences are as follows:

$$w_{\mathrm{ref}.}={\alpha }_{\mathrm{ref}}q{a}^4/100D$$

(28)

With: ${\alpha }_{\mathrm{ref}}=\alpha \_\mathrm{ref}\_60°=0.7945.{10}^{-2};{\alpha }_{\mathrm{ref}}=\alpha \_\mathrm{ref}\_30°=0.408.{10}^{-3}$ ;

$$D=E{h}^3/12(1-{\nu }^2)$$

The results were obtained in Tables 12 and 13 and Figures 24 and 25. It can be seen from these results that the GTMTPFE element exhibits a good behavior with respect to these tests and a rapid convergence for the different simulated cases. This performance is more significant for the slender plate. For the 30° bias plate, the configuration of the mesh arranged around the axes of symmetry of the plate (diagonal) permitted significantly faster convergence. Indeed, for a mesh density of 28 elements, the relative error is less than 2% for the L/h ratios equal to 50 and 100 and near to 5% for L/h equal to 1000.

Table 12.

60° skew plate simply supported on two opposite sides – displacement to the center.

L/h	Mesh	α_ref_60°	α_cal_60°	Wref	Wcal	Ratio Wcal/Wref
1000	4	0.7945	0.4483	8.676E+06	4.896E+06	0.56
	8		0.6831	8.676E+06	7.459E+06	0.86
	16		0.7632	8.676E+06	8.334E+06	0.96
100	4	0.7945	0.4697	8.676E+03	5.130E+03	0.59
	8		0.7208	8.676E+03	7.871E+03	0.91
	16		0.7889	8.676E+03	8.615E+03	0.99
50	4	0.7945	0.5000	1.084E+03	6.825E+02	0.63
	8		0.7656	1.084E+03	1.045E+03	0.96
	16		0.7945	1.084E+03	1.079E+03	1.00

Table 13.

30° skew plate simply supported on all sides - displacement to the center.

L/h	Mesh	α_ref_30°	α_cal_30°	Wref	Wcal	Ratio Wcal/Wref
1000	12	0.4080	0.3691	4.455E–01	4.030E–01	0.90
	28		0.3874	4.455E–01	4.230E–01	0.95
100	12	0.4080	0.3744	4.455E–04	4.089E–04	0.92
	28		0.4020	4.455E–04	4.390E–04	0.99
50	12	0.4080	0.3895	5.569E–05	5.316E–05	0.95
	28		0.4154	5.569E–05	5.670E–05	1.02

Figure 24.

60° skew plate simply supported on two opposite sides, relative error (%) on center displacement.

Figure 25.

30° skew plate simply supported on all sides, relative error (%) on center displacement.

Conclusion

In this study, we present a triangular finite element for moderately thick plate. Its formulation was been based on the strain approach, on solutions of Airy’s functions and on the analytical integration in the construction of the stiffness matrix.

The strengths associated with this approach consist of:

• The automatic verification of equilibrium conditions and kinematic compatibility conditions,

• The enrichment of the degrees of the interpolation polynomials of displacements, strain, and constraints (refinement p),

• The consideration distortions sections related to Poisson effects,

• The treatment of blocking phenomena related to transverse shear.

This approach has resulted in a competitive, robust and efficient plate finite element GTMTPFE:

• The strain approach for its functional discretization has given the possibility of enhancing the displacement fields, hence a greater precision in the approximation of the solution. This result is very marked in convergence tests, where we observe a rapid tendency toward the exact solution.

• The use of the analytic integration in the evaluation of the stiffness matrix gave this element robust behavior with respect to the geometrical distortion of the mesh. This result was significant in tests using distorted and irregular meshes, while for these same tests, the iso-parametric elements (using the numerical integration) are very sensitive (their convergence is conditioned by a regular, undistorted mesh).

• The using of the bi-harmonic solutions of the Airy function (by introducing the Poisson’s ratio as variable directly in the fields of interpolation of the kinematics) has allowed taking into account the distortion of the sections following the Poisson effects.

Based on our tests, the first findings indicate:

• The stability of the GTMTPFE element with respect to the patch tests carried out, particularly in its behavior with respect to the rigid modes and the modes representing the constant moments.

• A behavior marked by the absence of the transverse shear locking.

• The suitability of the GTMTPFE element for the treatment of both moderately thick and thin structures.

This is visible through its stability in the face of patch tests (constant twists, state of constants moments, transverse shear locking phenomenon, isotropy test), on the one hand, and in view of its excellent degree of convergence toward the reference solution, on the other hand. It exhibits superior performance behavior to other existing triangular plate finite elements in the literature, particularly for moderately thick plates and for thin plates (L/h ratio greater than 4).

Thus, we present an element that has convincing advantages that call for its use. In perspective, its further development remains necessary. It will be interesting to study its behavior in stresses, in dynamics, in buckling and to reinforce the performances observed here through additional validation tests.

In addition, the GTMTPFE element which treat the bending behavior problems can be combined with the membrane element developed by Himeur et al. ³² to treat problems of complex structures shell.