A Simple High-Order Shear Deformation Triangular Plate Element with Incompatible Polynomial Approximation

Due to the mathematical complexity raised by a high continuity requirement, developing simple/efﬁcient standard ﬁnite elements with general polynomial approximations applicable for arbitrary HSDTs seems to be a difﬁcult task at the present theoretical level. In this article, a series of High-order Shear Deformation Triangular Plate Elements (HSDTPEs) are developed using polynomial approximation for the analysis of isotropic thick-thin plates, through-thickness functionally graded plates, and cracked plates. The HSDTPEs have the advantage of simplicity in formulation, are free from shear locking, avoid using a shear correction factor and reduced integration, and provide stable solutions for thick and thin plates. The work can be further applied to plates and shells analysis with arbitrary shapes of elements, as well as more general problems related to the shear deformable effect, such as fracture and functionally graded plates. Abstract: The High-order Shear Deformation Theories (HSDTs) which can avoid the use of a shear correction factor and better predict the shear behavior of plates have gained extensive recognition and made quite great progress in recent years, but the general requirement of C 1 continuity in approximation ﬁelds in HSDTs brings difﬁculties for the numerical implementation of the standard ﬁnite element method which is similar to that of the classic Kirchhoff-Love plate theory. As a strong complement to HSDTs, in this work, a series of simple High-order Shear Deformation Triangular Plate Elements (HSDTPEs) using incompatible polynomial approximation are developed for the analysis of isotropic thick-thin plates, cracked plates, and through-thickness functionally graded plates. The elements employ incompatible polynomials to deﬁne the element approximation functions u / v / w , and a ﬁctitious thin layer to enforce the displacement continuity among the adjacent plate elements. The HSDTPEs are free from shear-locking, avoid the use of a shear correction factor, and provide stable solutions for thick and thin plates. A variety of numerical examples are solved to demonstrate the convergence, accuracy, and robustness of the present HSDTPEs.


Introduction
The classic Kirchhoff-Love plate theory based on the assumption that a plane section perpendicular to the mid-plane of the plate before deformation remains plane and perpendicular to the deformed mid-plane after deformation is the simplest plate theory in engineering analysis.However, the Kirchhoff-Love plate theory is only applicable for thin plates due to the neglecting of the shear deformation effects.The most well-known and earliest plate theories that take into account the shear deformation effects were proposed by Reissner [1] and Mindlin [2], in which the Mindlin plate theory was based on an assumption of a linear variation of in-plane displacements through the thickness of the plate, referred to as the First-order Shear Deformation Theory (FSDT).The plate elements derived from the FSDT only require C 0 continuity in approximation fields, have the advantages of physical clarity and simplicity of application [3], and hence were widely accepted and used to model thick-thin plates by scientists and engineers.Unfortunately, the FSDT elements suffer from the shear-locking problem when the thickness to length ration of the plate becomes very small, due to inadequate dependence among transverse deflection and rotations using an ordinary low-order finite element [4].Quite a large number of techniques have been developed to overcome this problem, such as the assumed shear strain approach, the discrete Kirchhoff/Mindlin representation, the mixed/hybrid formulation, and the reduced/selected integration [5][6][7][8][9][10][11][12][13][14][15].These formulations are free from shear locking and are applicable to a wide range of practical engineering problems, but in general, it is rather complex and time consuming to include the transverse shear effects for thick plates, which would also lead to complexity and difficulty in the programming.Moreover, the assumption of FSDT causes constant transverse shear strains and stresses across the thickness, which violates the conditions of zero transverse shear stresses on the top and bottom surfaces of plates.A shear correction factor is therefore required to properly compute the transverse shear stiffness.The finding of such a shear correction factor in FSDT is difficult since it depends on geometric parameters, material, loading and boundary conditions, etc. [16].
In recent years, High-order Shear Deformation Theories (HSDTs) have gained extensive recognition and made quite great progress [4,.Based on polynomial or non-polynomial transverse shear functions, various HSDTs have been proposed to avoid the use of a shear correction factor, and to better predict the shear behavior of the plate, for instance, the third-order shear deformation theory [17,18], the fifth-order shear deformation theory [19], the exponential shear deformation theory [20], the hyperbolic shear deformation theory [21], and the combined or mixed HSDT [22,23].Please see Thai and Kim [16] and Caliri et al. [24] for a comprehensive review of HSDTs.In HSDTs, the bending angles of rotation and shear angles can be treated as independent variables, and the shear-locking problem encountered in FSDT can be well-solved [4].In [30][31][32][33][34][35][36][37][38][39][40][41][42][43], two well-know HSDTs named as equivalent single layer (ESL) and layer-wise (LW) models are developed to evaluate the effective mechanical behavior of composite structures correctly.The accuracy and reliability of HSDTs have been illustrated by numerous examples in the literature [4,17,[25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43].However, the general requirement of C 1 continuity in approximation fields in HSDTs brings difficulties for the numerical implementation of the standard Finite Element Method (FEM), which is similar to that of the classic Kirchhoff-Love plate theory.Most examples in the literature are focused on the analytical/numerical solutions of simple Navier-type or Levy-type square plates.The numerical examples reported for the C 1 rectangular finite element using Lagrange interpolation and Hermite interpolation proposed by Reddy [25] and the C 0 continuous isoparametric Lagrangian finite element with 63 Degrees Of Freedom (DOFs) per element proposed by Gulshan et al. [44] are also limited to the rectangular plate or skew plate.Owing to the striking feature of capturing the high-order continuity well, the Meshless Methods (MM) and IsoGeometric Analysis methods (IGA) appear to be suitable potential methods to construct the numerical formulations for the plate based on HSDTs.The successful implementation of MM [45][46][47][48] and IGA [19,23,[49][50][51][52][53][54] in a number of thick-thin plates with arbitrary geometries can be found in the literature.
From the above literature review, it is observed that, due to the mathematical complexity raised by the high continuity requirement, developing simple/efficient standard finite elements with general polynomial approximations applicable for arbitrary HSDTs seems to be a difficult and unreachable task at the present theoretical level.In Cai and Zhu [55], a locking-free MTP9 (Mindlin type Triangular Plate element with nine degrees of freedom) using incompatible polynomial approximation is proposed.It also provides a new way and methodology to develop simple and efficient plate/shell elements based on HSDTs.In this work, with a similar procedure as the MTP9, a series of simple High-order Shear Deformation Triangular Plate Elements (HSDTPEs) using incompatible polynomial approximation are developed for the analysis of isotropic thick-thin plates and through-thickness functionally graded plates.In the HSDTPEs, different orders of general polynomials can be easily employed as element approximation functions, the displacement continuity among the adjacent plate elements can be equivalently enforced by a fictitious thin layer which has a definite physical meaning, and consequently, there are no continuity requirements under the theoretical framework of the present HSDTPEs.The HSDTPEs avoid the shear-locking problem and the use of a shear correction factor, and have a good convergence rate and high accuracy for both thick and thin plates.Several representative numerical examples are solved and compared to validate the performance of the present HSDTPEs.

Incompatible Polynomial Approximation over Each Triangular Element
Consider a linear elastic plate with a length a, width b, and thickness h undergoing infinitesimal deformation, as illustrated in Figure 1.The mid-plane of the plate is divided into arbitrary triangular elements, as shown in Figure 2. The displacement function of the most well-known HSDTs [17] for each triangular element e i is generally defined by: where u 0 , v 0 and w 0 are the in-plane and transverse displacements at the mid-plane, respectively; u = [u, v, w] T denotes the displacements of a point x on the plate; θ x and θ y are the rotations of the normal to the cross section; z is the coordinate in the transverse direction; and g(z) describes the distribution of shear effect in the thickness direction.A review of transverse shear functions g(z) can be found in Nguyen et al. [56].For isotropic plates with infinitesimal strains, the in-plane displacements u 0 (x, y) and v 0 (x, y) can be neglected because, the thickness h is much smaller than the characteristic length a and b, and the transverse displacement is much smaller than the thickness h, which leads to u 0 (x, y) ≈ 0 and v 0 (x, y) ≈ 0 at the mid-plane.The transverse normal displacement w can also be assumed as w = w(x, y, z), which is not a constant along the z axis, and can be defined by the ESL or LW models [30][31][32][33][34][35][36][37][38][39][40][41][42][43] to capture the effective mechanical behavior along the thickness of composite structures well.To demonstrate the performance of the present theory for various transverse shear functions, the third-order shear function g(z) = − 4z 3 3h 2 [17] and the fifth-order shear function g(z) = − z 8 − 2z 3 h 2 + 2z 5 h 4 [19] are used to develop the Third-order Shear Deformation Triangular Plate Element (TrSDTPE) and Fifth-order Shear Deformation Triangular Plate Element (FfSDTPE), respectively.For the special case g(z) = 0, Equation ( 1) is actually the expression of Mindlin plate theory (or FSDT).The corresponding plate element is referred to as FiSDTPE (First-order Shear Deformation Triangular Plate Element) for comparison in the paper.
We assume that: where , and a w = a 25 a 26 • • • a 34 T are the vector of generalized approximation DOFs (degrees of freedom) of the triangular element e i ; P 2 is the second-order polynomial basis function; and P 3 is the third-order polynomial basis function in which: where x 0 = x − x i , y 0 = y − y i , (x i , y i ) are the coordinates of the central point of element e i .It should be noted that only triangular elements, as well as second-order and third-order polynomial functions, are implemented in the paper, but actually, the arbitrary shape of elements and arbitrary orders of polynomials can also be easily employed to derive high-order shear deformation plate elements in the present work.
Substituting Equation (2) into Equation (1), the displacement approximation over element e i can be further expressed as: where α z = −z − g(z), P 3 ,x = ∂P 3 ∂x , P 3 ,y = ∂P 3 ∂y , ) The strain-displacement relations of the linear elastic problem are given by: Substituting Equation ( 5) into Equation ( 8), we have: where ε = ε x , ε y , γ xy , γ yz , γ xz T is the strain vector and B is the strain matrix, where: L is a differential operator where: For an isotropic linear elastic material, the stress-strain relations in element e i are given by: where σ = σ x , σ y , τ xy , τ yz , τ xz T , the transverse stress σ z is assumed to be ignored for plate structures, and the elasticity matrix is: where D 0 = E 1−v 2 , E is the elastic modulus and v is the Poisson ratio.As mentioned above, a shear correction factor k is required to properly compute the transverse shear stiffness in the FiSDTPE with the assumption of FSDT, which causes constant transverse shear strains and stresses across the thickness and violates the conditions of zero transverse shear stresses on the top and bottom surfaces of plates [16].Usually, k is taken as k = 1.2 for the special case of FiSDTPE in Equation ( 13) according to the principle of the equivalence of strain energy.However, the high-order shear deformation theory gives a parabolic distribution of the transverse stresses/strains directly and avoids the use of a shear correction factor, and thus k is taken as k = 1.0 in Equation ( 13) for the rest of the HSDTPEs.
Therefore, the strain energy of element e i can be derived as: For the plate made of Functionally Graded (FG) materials which is created by mixing two distinct material phases, the composition of the FG materials is in general assumed to be varied continuously from the top to the bottom surface.There are many kinds of FG materials made from all classes of solids.But for the sake of simplicity and convenience, only a ceramic-metal composite is considered and implemented to test the performance of the HSDTPEs in the present study, and the power-law [25,32,45] is used to describe the through-the-thickness distribution of FG materials, which is expressed as: where n is the volume fraction exponent, V c is the volume fraction of the ceramic, P m represents the material property of the metal, P c represents the material property of the ceramic, and P denotes the effective material property.In this work, the Young's modulus E in Equation ( 13) varies according to Equation ( 16) and the Poisson ratio v is assumed to be constant for the analysis of functionally graded plates.

Fictitious Thin Layer between Adjacent Triangular Elements
According to the definition of the displacement approximation in Equation ( 5), the deformation along the share boundary of the adjacent elements e i and e j is discontinuous, which means that u e i x p , y p , z p = u e j x p , y p , z p for an arbitrary point p along the share boundary shown in Figure 2, where point p has local coordinates s p , n p , z p and global coordinates x p , y p , z p .Here, we introduce a fictitious thin layer e l shown in Figure 3 to enforce the continuous condition over the share boundary of the elements.The geometry dimensions of e l are the length l, width d, and height h, where d l, d h, and h is the thickness of the plate in the transverse direction.Because d l and d h, the strain-displacement relations ε l = [γ ns , ε n , γ nz ] T in thin layer e l can be simplified as: where is the displacement of point p in the local coordinate (s, n, z) computed by the approximation of triangular element e i , and u p j = u p j , v p j , w p j T is the displacement of point p in local coordinate (s, n, z) computed by the approximation of triangular element e j .u p i and u p j can be calculated using Equation (5).Substitution of Equation ( 5) into Equation ( 17) yields: where where N e i x p , y p , z p is the shape function of point p in triangular element e i and N e j x p , y p , z p is the shape function of point p in triangular element e j .λ l is the transformation matrix of point p from the global coordinate x p , y p , z p to the local coordinate s p , n p , z p , where: and where a e i is the DOFs of element e i and a e j is the DOFs of element e j .
The stress-strain relations in fictitious thin layer e l are then given by: where σ l = [τ ns , σ n , τ nz ] T and where and k = 1.2.Similar to Equation ( 13), the shear correction factor k is taken as k = 1.2 for the special case of FiSDTPE, and k = 1.0 for the rest of the high-order shear deformation plate elements, including the TrSDTPE and FfSDTPE, without the need to use a shear correction factor.
The width d of fictitious thin layer e l is an important artificial parameter for the present HSDTPEs, but it is easy to select a reasonable d to satisfy d l and d h for simplifying the strain-displacement relations of thin layer e l in Equation (17).Numerical studies show that the variation of d in a large range has little effect on the accuracy of the calculation results.In this paper, width d is taken as d = 0.0001l.
Thus, the strain energy of thin layer e l can be derived as:

Imposing Displacement Boundary Condition
As illustrated in Figure 4, along boundary 1-2 of element e k , rotations θ s , θ n in the local coordinate (s, n, z) or the displacements (u 0 , v 0 , w 0 ) in the local coordinate (s, n, z) are fixed, where (u 0 , v 0 , w 0 ) represents the in-plane and transverse displacements at the mid-plane.A fictitious thin layer e b over boundary 1-2 shown in Figure 4 is also introduced to enforce the displacement boundary condition.We divide the displacement approximation in Equation ( 5) into two parts to the follow to separately enforce the rotation and mid-plane displacement boundaries.
where u p r represents the displacement function of point p for rotations in triangular element e k , u p t represents the displacement function of point p for mid-plane displacements in triangular element e k , and a r and a t are the corresponding DOFs of element e k .
where "0" in bold in Equations ( 27) and ( 28) represents zero matrix , and λ b is the transformation matrix similar to Equation (20) where Using the same derivation process as in Equation ( 24), the strain energy of the thin layer e b is derived as: where D b is calculated using Equation ( 23), N r = λ b N r , and Please refer to Cai and Zhu [55] for the detailed derivation of the displacement boundary condition fixed at a point or the given displacement boundary condition.

Load Boundary Condition
A distributed force f 0 = 0 0 f z (x, y) T along the transverse direction z is applied at element e d , as illustrated in Figure 5.By using Equation ( 5), the external force potential energy of element e d is written as: where N e d is the shape function of element e d calculated by Equation ( 5), and a e d is the DOFs of element e d .Similarly, a distributed resultant moment M 0 = M s M sn 0 T is applied to the edge of element e m , as illustrated in Figure 6.By using Equation ( 5), the external force potential energy of element e m is written as: where a e m is the DOFs of element e m , and N e m is the shape function of element e m corresponding to the moment, where: and the transformation matrix Figure 6.Moment boundary condition.

Equilibrium Equation
From Equations ( 14), ( 24), ( 30), (31), and (32), the total potential energy of a plate is obtained as: The variation of total potential energy Π results in the following discrete equation: where Assembling the above stiffness matrix and force vector, the equilibrium equation for a plate is then obtained as: where K is the global stiffness matrix, F is the force vector, and U is the vector of DOFs to be solved.
As described in Senjanović et al. [4], the shear-locking problem could be well and naturally solved because the bending angles of rotation and shear angles are treated as independent variables in HSDTs.The regular full integration can be applied to make HSDTPEs valid for the thick-thin plates for the computation of Equation ( 42), for instance, seven quadrature points for each triangular element [57], fou Gauss quadrature points for transverse direction z (where the analytical integration can also be applied for the direction z), and four Gauss quadrature points for the local direction s of each fictitious layer are used for the integration of the TrSDTPE using the third-order shear function g(z) [58][59][60].

Analysis of Cracked Plates
The present HSDTPEs are also applied to the calculation of Stress Intensity Factors (SIFs) of cracked thick-thin plates.As illustrated in Figure 7, the mid plane of a cracked plate is taken as the x-y plane and is divided into arbitrary triangular elements.Accurate computation of SIFs remains challenging in the field of fracture mechanics.For plates loaded by a combination of bending and tension, the SIFs can also be computed by the Virtual Crack Closure Technique (VCCT) [61][62][63], the path-independent J-integral technique or interaction integral [64,65], and the stiffness derivative method [66].In this paper, the Virtual Crack Closure Technique (VCCT) [61][62][63] is employed to calculate the SIFs of the cracked plate.For the convenience of implementing the VCCT, point T 2 shown in Figures 7 and 8 is temporarily moved to T 3 along the extended line direction of T 1 − T. In the local coordinates (s, n, z) shown in Figure 9, the relative displacements [∆u(s, z), ∆v(s, z), ∆w(s, z)] of T 1 − T and the stresses [τ ns (s, z), σ n (s, z), τ nz (s, z)] of T − T 1 can be easily calculated using Equations ( 18) and ( 22) for the fictitious thin layer e l .For example, assuming that the T 1 − T is simulated by a fictitious thin layer with width d shown in Figure 3, the relative displacements can be evaluated by ∆u(s, z) = γ ns d, ∆v(s, z) = ε n d and ∆w(s, z) = γ nz d.Then, the energy release rate at crack tip T is obtained by the VCCT as: where G I is the energy release rate of crack mode I, G II is that of crack mode II, and G III is that of crack mode III.Then, the SIFs of the crack tip can be computed by means of the relations between the energy release rate and SIFs for the plate theory, for instance, K 1 = √ 3EG I [63].

Simply Supported Square Plate Subjected to Uniform Load
A simply supported square plate subjected to a uniform load q is tested to show the reliability and convergence of the present elements.The side length of the plate is L, and the thickness of the plate is h.A quarter of the plate is modeled as a result of symmetry, as illustrated in Figure 10.For the isotropic plates, the in-plane displacements (u 0 , v 0 ) and their DOFs in Equations ( 1), (2), and (5) are neglected in the following analyses.The displacement boundary conditions of the present theory along the simply supported edges in local coordinates are θ s = 0 and w = 0.The n × n regular mesh and irregular mesh illustrated in Figures 11 and 12 are employed for convergence studies.The elements DST-BL (Discrete Mindlin triangular plate element) [7] and RDKTM (Re-constituting discrete Kirchhoff triangular plate element) [14] have been selected for comparison with the present elements based on HSDTs.The reference solutions in the following Tables 1-6 are taken from Long et al. [67], which are also labeled as analytical solutions in [67].Table 1 lists   12(1−v 2 ) .Table 2 reports the normalized bending moment M 0 = M c / qL 2 10 of the simply supported square plate, where M c is the central bending moment of the plate.The convergence of the deflection for the simply supported square plate using different elements when the aspect ratio h/L = 0.1 is shown in Figure 13.It is observed that all the present TrSDTPE, FfSDTPE, and FiSDTPE shows a good convergence rate and high accuracy, and avoids the shear-locking problem.The results also indicate that the present elements are insensitive to element distortions of the irregular mesh shown in Figure 12.The width-to-length ratio d/l of the fictitious thin layer e l plays an important role in the present formulations.Table 3 reports the effect of the ratio d/l on the normalized defection W 0 for the simply supported square plate, where 16 × 16 regular mesh and an aspect ratio of h/L = 0.1 are employed.The results in Table 3 indicate that the artificial parameter d/l has little effect on the solution accuracy when d/l ≤ 0.001, and it is easy to select a reasonable d in the current formulation.In this work, width d is taken as d = 0.0001l.The condition numbers of the global stiffness matrices of the simply supported square plate using the present elements are also computed and reported in Figure 14.As seen, the variation of the condition number in Figure 14 reflects that the present elements show a good conditioning and stability in the case of mesh refinement.A clamped square plate subjected to a uniformly distributed load q is further investigated to test the performance of the present elements for clamp boundary conditions.The geometry and material parameters of the clamped plate are the same as those of the above simply supported plate.The displacement boundary conditions of the present theory along the clamped edges in Figure 10 in local coordinates are θ s = 0, θ n = 0 and w = 0.
The results for the normalized central deflection W 0 of the clamped square plate are compared in Table 4.It is seen that, for the plate with the clamped boundary conditions, the predictions of FiSDTPE, DST-BL, and RDKTM based on FSDT agree well with the reference solutions [67] for plates, but the TrSDTPE and FfSDTPE based on HSDTs seem to underestimate the deflections compared with the reference solutions [67] for thick plates of h/L ≥ 0.2.To further illustrate the accuracy of the present shear elements, the comparisons of the predictions by different elements and the solutions by 3D elasticity FEM software ANSYS using 20-nodes hexahedron isoparametric element and an element side length of 0.05 are listed in Table 5.By taking the 3D FEM solutions as the benchmark, Table 5 indicates that the present TrSDTPE and FfSDTPE show better solution accuracy than the elements DST-BL, RDKTM, and FiSDTPE based on FSDT for the plates involving clamp boundaries.Moreover, the TrSDTPE with a third-order shear function g(z) shows the best solution accuracy among all the elements for the clamped plate.
The DOFs of the different methods employed in Table 5 are compared by taking the case of h/L = 0.3 for the clamped plate.Assuming that a 16 × 16 regular mesh is employed for the plate element discretization and a 16 × 16 × 9 regular mesh for the 20-nodes hexahedron 3D element discretization, we can see that the total DOFs of the DST-BL/RDKTM element, HSDTs, and 3D 20-nodes hexahedron element are 867, 17,408, and 41,337, respectively.It is seen that the total DOFs and the efficiency of HSDTs are between the DST-BL/RDKTM and the 3D FEM.Although the number of DOFs only decreases to 42% of the 3D FEM method, the present 2D HSDTs have the advantage of simplicity and flexibility in the mesh generation compared with 3D FEM for the plates with different thicknesses.Compared with other 2D plate elements such as DST-BL and RDKTM, the computational DOFs of the present 2D HSDTs seem to be relatively higher, but the formulation and the numerical implementation of the high-order shear deformation theory in the present HSDTs are much simpler than those of the DST-BL/RDKTM.From the point of view of the 2D analysis, the total computational cost of the present elements is bearable and worthy in terms of its advantages in formulation and implementation.HSDTs which have almost the same computational efficiency of DST-BL/RDKTM could also be constructed using the reduced integral method similar to our previous work [55], but the present HSDTs avoiding the reduced integration by paying a certain computational cost are more practicable in engineering analysis.

Clamped Circular Plate Subjected to Uniform Load
A clamped circular plate subjected to a uniformly distributed load q is taken into consideration in this section.The thickness of the plate is h.The radius of the plate is r = 100.A quarter of the plate with symmetry conditions on axes x and y is modeled in Figure 15 6, where W c is the central deflection of the circular plate.Again, an excellent agreement between the present solutions and the reference solutions is observed for this problem.7, along with the reference solutions by Tanaka et al. [68] and Boduroglu et al. [69] based on FSDT for comparison.In Tanaka et al. [68], a cracked plate is analyzed by employing the mesh-free reproducing kernel approximation formulated by Mindlin-Reissner plate theory, and the moment intensity factor is evaluated by the J-integral with the aid of nodal integration.In Boduroglu et al. [69], the crack problem is solved by the dual boundary element method based on Reissner plate formulation, and the stress resultant intensity factor is calculated by employing the J-integral techniques.The SIFs in Table 7 are It is observed that the present elements show a high solution accuracy for the calculation of the SIFs.

Symmetric Edge Cracks in a Rectangular Plate
As illustrated in Figure 18, a rectangular plate with symmetric double edge cracks is analyzed in this example.The geometry dimensions, material properties, and boundary conditions are the same as those of the center cracked problem described in Section 4.4.The crack length is a and the plate thickness is h.Division of 2728 triangular elements shown in Figure 17 is also employed for this analysis.The normalized SIFs obtained by the present elements for different values of d/b and b/h are presented in Tables 8 and 9, along with reference solutions [68,69].The numerical methods and plate theories for solving the problem are the same as the above rectangular plate problem involving a center crack.As expected, the present results are in good agreement with the reference solutions.In this section, a simply supported FG plate subjected to a uniformly distributed load q is analyzed and compared.The FG plate is comprised of aluminum (E m = 70 GPa, v m = 0.3) and ceramic (E c = 151 GPa, v c = 0.3).The side length of the plate is L = 1m, and the thickness of the plate is h.A quarter of the FG plate is modeled as a result of symmetry, as shown in Figure 19.The displacement boundary conditions for the symmetric and simply supported sides in local coordinates are also illustrated in Figure 19.The 16 × 16 regular mesh similar to Section 4.1 is employed for computation.Tables 10 and 11 list the normalized defection W 0 of the FG plate for different aspect ratios h/L and different exponents n in Equation (15), where W 0 = W c / qL 4 E m h 3 and W c is the central deflection of the plate.In Ferreira et al. [45], the FG plate is solved by the meshless collocation method with multiquadric radial basis functions and a third-order shear deformation theory.The problem is also solved by Talha and Singh [70] using the C 0 isoparametric finite element with 13 degrees of freedom per node, and the power-law similar to Equations ( 15) and ( 16) is used to describe the through-the-thickness distribution of FG materials in the HSDT model.The results obtained by the present TrSDTPE, FfSDTPE, and FiSDTPE are in good agreement with the meshless solutions of Ferreira et al. [45], which compute the effective elastic moduli by the rule of mixture.

Conclusions
In this work, a series of novel HSDTPEs using incompatible polynomial approximation are developed for the analysis of isotropic thick-thin plates and through-thickness functionally graded plates.The HSDTPEs are free from shear-locking, avoid the use of a shear correction factor, and provide stable solutions for thick and thin plates.The present formulation, which defines the element approach with incompatible polynomials and avoids the need to satisfy the requirement of high-order continuity in approximation fields in HSDTs, also provides a new way and methodology to develop simple plate/shell elements based on HSDTs.The accuracy and robustness of the present elements are well demonstrated through various numerical examples.
Only two types of HSDTPEs including TrSDTPE and FfSDTPE, and one special type of first-order shear deformation triangular plate element FiSDTPE, have been studied and discussed in the paper.The present formulation can be further extended to plates and shells with arbitrary shapes of elements, and further applied to more general problems related to the shear deformable effect such as the thermomechanical, vibration, and buckling analysis of functionally graded plates and laminated/sandwich structures.

Figure 2 .
Figure 2. Triangular elements for the mid-plane of the plate.

Figure 3 .
Figure 3.A fictitious thin layer between adjacent triangular elements.

Figure 7 .
Figure 7. Triangular elements for a cracked plate.

Figure 8 .
Figure 8. Minor movement for the implementation of VCCT.

Figure 10 .
Figure 10.A quarter model of the square plate.
the normalized defection W 0 = W c / qL 4 100D b of the simply supported square plate, where W c is the central deflection of the plate and D b = Eh

Figure 12 .
Figure 12.Irregular mesh for the square plate.

Figure 13 .
Figure 13.Convergence of normalized deflection for the simply supported square plate.

Figure 14 .
Figure 14.Variation of condition number versus number of elements.
. The displacement boundary conditions of the present theory along the clamped edges in local coordinates are θ s = 0, θ n = 0, and w = 0. Divisions of 32, 128, 512, and 2048 triangular elements are employed for the convergence studies.Typical meshes of 512 and 2048 triangular elements for the circular plate are shown in Figure 16.The results for normalized deflection W 0 = W c / qr 4 100D b of the clamped circular plate are listed in Table

Figure 15 .
Figure 15.Model of the circular plate.

Figure 16 .
Figure 16.Typical meshes for the circular plate.

4. 4 .
Rectangular Plate Involving a Center Crack Consider a rectangular plate involving a center crack as shown in Figure 17.The material properties are E = 1.0 × 10 6 Pa and v = 0.3.The width and length are 2b = 1m and 2c = 2m, respectively.The crack length is 2a and the plate thickness is h.Divisions of 2728 and 7338 triangular elements are employed for the calculation of SIFs of the center crack plate.The displacement and moment boundary conditions are also illustrated in Figure 17.The numerical results obtained by TrSDTPE, FfSDTPE, and FiSDTPE for different a/h values are reported in Table

Figure 19 .
Figure 19.Modal of the FG plate.

Table 1 .
Normalized deflection for the simply supported square plate.

Table 2 .
Normalized bending moment for the simply supported square plate.

Table 3 .
Convergence of normalized deflection with different width-to-length ratios.

Table 4 .
Normalized deflection for the clamped square plate.

Table 5 .
Comparisons of the normalized deflections with 3D FEM solutions for thick plates.
Note: Value in parentheses is the relative error with respect to 3D FEM.

Table 6 .
Normalized deflection for the clamped circular plate.

Table 7 .
Normalized SIFs F 1 for the center cracked plate.
Figure 17.Model of the center cracked plate.
Figure 18.Model of symmetric edge cracks.