Quadratic Solid–Shell Finite Elements for Geometrically Nonlinear Analysis of Functionally Graded Material Plates

In the current contribution, prismatic and hexahedral quadratic solid–shell (SHB) finite elements are proposed for the geometrically nonlinear analysis of thin structures made of functionally graded material (FGM). The proposed SHB finite elements are developed within a purely 3D framework, with displacements as the only degrees of freedom. Also, the in-plane reduced-integration technique is combined with the assumed-strain method to alleviate various locking phenomena. Furthermore, an arbitrary number of integration points are placed along a special direction, which represents the thickness. The developed elements are coupled with functionally graded behavior for the modeling of thin FGM plates. To this end, the Young modulus of the FGM plate is assumed to vary gradually in the thickness direction, according to a volume fraction distribution. The resulting formulations are implemented into the quasi-static ABAQUS/Standard finite element software in the framework of large displacements and rotations. Popular nonlinear benchmark problems are considered to assess the performance and accuracy of the proposed SHB elements. Comparisons with reference solutions from the literature demonstrate the good capabilities of the developed SHB elements for the 3D simulation of thin FGM plates.


Introduction
Over the last decades, the concept of functionally graded materials (FGMs) has emerged, and FGMs were introduced in the industrial environment due to their excellent performance compared to conventional materials. This new class of materials was first introduced in 1984 by a Japanese research group, who made a new class of composite materials (i.e., FGMs) for aerospace applications dealing with very high temperature gradients [1,2]. These heterogeneous materials are made from several isotropic material constituents, which are usually ceramic and metal. Among the many advantages of FGMs, their mechanical and thermal properties change gradually and continuously from one surface to the other, which allows for overcoming delamination between interfaces as compared to conventional composite materials. In addition, FGMs can resist severe environment conditions (e.g., very high temperatures), while maintaining structural integrity.
Thin structures are widely used in the automotive industry, especially through sheet metal forming into automotive components. In this context, the finite element (FE) method is considered nowadays as a practical numerical tool for the simulation of thin structures. Traditionally, shell and solid elements are used in the simulation of linear and nonlinear problems. However, the simulation results require very fine meshes to obtain accurate solutions due to the various locking phenomena

Quadratic Approximation for the SHB Elements
Conventional quadratic interpolation functions for traditional continuum prismatic and hexahedral elements are used in the formulation of the SHB elements. According to this formulation, the spatial coordinates i x and the displacement field i u are approximated using the following interpolation functions: where iI d are the nodal displacements, 1, 2, 3 i = correspond to the spatial coordinate directions, and I varies from 1 to K, with K being the number of nodes per element, which is equal to 15 for the SHB15 element and 20 for the SHB20 element.

Strain Field and Gradient Operator
Using the above approximation for the displacement within the element, the linearized strain tensor ε can be derived as: By combining Equations (1) and (2) with the help of the interpolation functions, the nodal displacement vectors i d write: where ( ) 1 2 3 , , , ,  are the nodal coordinate vectors. In Equation (4), index α goes from 1 to 11 for the SHB15 element, and from 1 to 16 for the SHB20 element. In addition, vector ( ) 1, 1, , 1 T = s  has fifteen constant components in the case of the SHB15 element, and twenty constant components vector for the SHB20 element. Vectors α h are composed of h α functions, which are evaluated at the element nodes, and the full details of their expressions can be found in [23].
With the help of some well-known orthogonality conditions and of the Hallquist [24] vectors

Quadratic Approximation for the SHB Elements
Conventional quadratic interpolation functions for traditional continuum prismatic and hexahedral elements are used in the formulation of the SHB elements. According to this formulation, the spatial coordinates x i and the displacement field u i are approximated using the following interpolation functions: where d iI are the nodal displacements, i = 1, 2, 3 correspond to the spatial coordinate directions, and I varies from 1 to K, with K being the number of nodes per element, which is equal to 15 for the SHB15 element and 20 for the SHB20 element.

Strain Field and Gradient Operator
Using the above approximation for the displacement within the element, the linearized strain tensor ε can be derived as: By combining Equations (1) and (2) with the help of the interpolation functions, the nodal displacement vectors d i write: x iK ) are the nodal coordinate vectors. In Equation (4), index α goes from 1 to 11 for the SHB15 element, and from 1 to 16 for the SHB20 element. In addition, vector s T = (1, 1, · · · , 1) has fifteen constant components in the case of the SHB15 element, and twenty constant components vector for the SHB20 element. Vectors h α are composed of h α functions, which are evaluated at the element nodes, and the full details of their expressions can be found in [23].
With the help of some well-known orthogonality conditions and of the Hallquist [24] vectors , where vector N contains the expressions of the interpolation functions N I , the unknown constants a ji and c αi in Equation (4) can be derived as: where the complete details on the expressions of vectors γ α can be found in [22]. By introducing the discrete gradient operator B, the strain field in Equation (3) writes: where the expression of the discrete gradient operator B is:

Hu-Washizu Variational Principle
The SHB solid-shell elements are based on the assumed-strain method, which is derived from the simplified form of the Hu-Washizu variational principle [25]. In terms of assumed-strain rate . ε, interpolated stress σ, nodal velocities . d, and external nodal forces f ext , this principle writes The assumed-strain rate is expressed in terms of the discrete gradient operator B as: Substituting the expression of the assumed-strain rate given by Equation (9) into the variational principle (Equation (8)), the expressions of the stiffness matrix K e and the internal forces f int for the SHB elements are where K GEOM is the geometric stiffness matrix. As to the fourth-order tensor C e (ζ), it describes the functionally graded elastic behavior of the FGM material. Its expression is given hereafter.

Description of Functionally Graded Elastic Behavior
In the framework of large displacements and rotations, the formulation of the SHB elements requires the definition of a local frame with respect to the global coordinate system, as illustrated in Figure 2. The local frame, which is designated as the "element frame" in Figure 2, is defined for each element using the associated nodal coordinates. In such an element frame, where the ζ-coordinate represents the thickness direction, the fourth-order elasticity tensor C e (ζ) for the FGM is specified. In this work, a two-phase FGM is considered, which consists of two constituent mixtures of ceramic and metal. The ceramic phase of the FGM can sustain very high temperature gradients, while the ductility of the metal phase prevents the onset of fracture due to the cyclic thermal loading. In such FGMs, the material at the bottom surface of the plate is fully metal and at the top surface of the plate is fully ceramic, as illustrated in Figure 3. Between these bottom and top surfaces, the elastic properties vary continuously through the thickness from metal to ceramic properties, respectively, according to a power-law volume fraction. The corresponding volume fractions for the ceramic phase c f and the metal phase m f are expressed as (see, e.g., [26,27]): where n is the power-law exponent, which is greater than or equal to zero, and , with t the thickness of the plate. For 0 n = , the material is fully ceramic, while when n → ∞ the material is fully metal (see Figure 4).   In this work, a two-phase FGM is considered, which consists of two constituent mixtures of ceramic and metal. The ceramic phase of the FGM can sustain very high temperature gradients, while the ductility of the metal phase prevents the onset of fracture due to the cyclic thermal loading. In such FGMs, the material at the bottom surface of the plate is fully metal and at the top surface of the plate is fully ceramic, as illustrated in Figure 3. Between these bottom and top surfaces, the elastic properties vary continuously through the thickness from metal to ceramic properties, respectively, according to a power-law volume fraction. The corresponding volume fractions for the ceramic phase f c and the metal phase f m are expressed as (see, e.g., [26,27]): where n is the power-law exponent, which is greater than or equal to zero, and z ∈ [−t/2, t/2], with t the thickness of the plate. For n = 0, the material is fully ceramic, while when n → ∞ the material is fully metal (see Figure 4). In this work, a two-phase FGM is considered, which consists of two constituent mixtures of ceramic and metal. The ceramic phase of the FGM can sustain very high temperature gradients, while the ductility of the metal phase prevents the onset of fracture due to the cyclic thermal loading. In such FGMs, the material at the bottom surface of the plate is fully metal and at the top surface of the plate is fully ceramic, as illustrated in Figure 3. Between these bottom and top surfaces, the elastic properties vary continuously through the thickness from metal to ceramic properties, respectively, according to a power-law volume fraction. The corresponding volume fractions for the ceramic phase c f and the metal phase m f are expressed as (see, e.g., [26,27]): where n is the power-law exponent, which is greater than or equal to zero, and , with t the thickness of the plate. For 0 n = , the material is fully ceramic, while when n → ∞ the material is fully metal (see Figure 4).    For an isotropic elastic behavior, the constitutive equations are governed by the Hooke elasticity law, which is expressed by the following relationship: where 1 denotes the second-order unit tensor, λ and μ are the Lamé constants given by:  For an isotropic elastic behavior, the constitutive equations are governed by the Hooke elasticity law, which is expressed by the following relationship: where 1 denotes the second-order unit tensor, λ and µ are the Lamé constants given by: with E and ν the Young modulus and the Poisson ratio, respectively. For FGMs that are made of ceramic and metal constituents, it is commonly assumed that only the Young modulus E varies in the thickness direction, while the Poisson ratio ν is kept constant. Therefore, the constant Young modulus in Equation (14) is replaced by E(z), whose value evolves according to the following power-law distribution: where E c and E m are the Young modulus of the ceramic and metal, respectively.

Nonlinear Benchmark Problems
In this section, the performance of the proposed elements is assessed through the simulation of several popular nonlinear benchmark problems. The static ABAQUS/Standard solver has been used to solve the following static benchmark problems. More specifically, the classical Newton method is considered for most benchmark problems, aside from limit-point buckling problems for which the Riks arc-length method is used.
To accurately describe the variation of the Young modulus through the thickness of the FGM plates, only five integration points within a single element layer is used in the simulations. For each benchmark problem, the simulation results given by the proposed elements are compared to the reference solutions taken from the literature. In the subsequent simulations, it is worth noting that the elastic properties of the metal and ceramic constituents of the FGM plates do not reflect a real metallic or ceramic material. Indeed, the terms metal and ceramic are commonly used in the literature to emphasize the difference between the properties of the FGM constituents (see, e.g., [15,26,27]).
Regarding the meshes used in the simulations, the following mesh strategy is adopted: (N 1 × N 2 ) × N 3 for the hexahedral SHB20 element, where N 1 is the number of elements along the length, N 2 is the number of elements along the width, and N 3 is the number of elements along the thickness direction. As to the prismatic SHB15 element, the mesh strategy consists of (N 1 × N 2 × 2) × N 3 , due to the in-plane subdivision of a rectangular element into two triangles. Figure 5a shows a simple cantilever FGM beam with a bending load at its free end. This is a classical popular benchmark problem, which has been widely considered in many works for the analysis of cantilever beams with isotropic material (see, e.g., [28,29]). The Poisson ratio of the FGM beam is assumed to be ν = 0.3, while the Young modulus of the metal and ceramic constituents are E m = 2.1 × 10 5 Mpa and E c = 3.8 × 10 5 Mpa, respectively. Figure 5b illustrates the final deformed shape of the cantilever beam with respect to its undeformed shape, as discretized with SHB20 elements, in the case of fully metallic material. Figure 6 shows the load-deflection curves obtained with the quadratic SHB elements, along with the reference solutions taken from [15], for various values of the power-law exponent n. One recalls that fully metallic material is obtained when n → ∞ , and fully ceramic material for n = 0. Overall, the SHB elements show excellent agreement with the reference solutions corresponding to the various values of exponent n. More specifically, it can be observed that the bending behavior of the FGM beam lies between that of the fully ceramic and fully metal beam, which is consistent with the power-law distribution of the Young modulus in the thickness direction. Another advantage of the proposed SHB elements is that, using the same in-plane mesh discretization as in reference [15], only five integration points through the thickness are sufficient for the SHB elements, while ten integration points have been considered in [15] to simulate this benchmark problem.  Figure 5b illustrates the final deformed shape of the cantilever beam with respect to its undeformed shape, as discretized with SHB20 elements, in the case of fully metallic material. Figure 6 shows the load-deflection curves obtained with the quadratic SHB elements, along with the reference solutions taken from [15], for various values of the power-law exponent n. One recalls that fully metallic material is obtained when n → ∞ , and fully ceramic material for 0 n = . Overall, the SHB elements show excellent agreement with the reference solutions corresponding to the various values of exponent n. More specifically, it can be observed that the bending behavior of the FGM beam lies between that of the fully ceramic and fully metal beam, which is consistent with the power-law distribution of the Young modulus in the thickness direction. Another advantage of the proposed SHB elements is that, using the same in-plane mesh discretization as in reference [15], only five integration points through the thickness are sufficient for the SHB elements, while ten integration points have been considered in [15] to simulate this benchmark problem.

Slit Annular Plate
In this section, the well-known slit annular plate problem is considered (see, e.g., [29][30][31]). The annular plate is clamped at one end and loaded by a line shear force P, as illustrated in Figure 7a Figure 7b illustrates the undeformed and deformed shapes of the annular plate, as discretized with SHB20 elements, in the case of fully metallic material. Figure 8 reports the load-out-of-plane vertical deflection curves at the outer point A of the annular plate as obtained with the SHB elements, along with the reference solutions taken from [15]. One can observe that the SHB elements perform very well with respect to the reference solutions for all considered values of exponent n. Similar to the previous benchmark problem, the same in-plane mesh discretization as in [15] with only five integration points through the thickness has been adopted by the proposed SHB elements for this nonlinear test, while ten integration points have been considered in [15].

Slit Annular Plate
In this section, the well-known slit annular plate problem is considered (see, e.g., [29][30][31]). The annular plate is clamped at one end and loaded by a line shear force P, as illustrated in Figure 7a. The inner and outer radius of the annular plate are R i = 6 m and R o = 10 m, respectively, while the thickness is t = 0.03 m. The Poisson ratio of the annular plate is ν = 0.3, while the Young modulus of the metal and ceramic constituents are E m = 21 Gpa and E c = 38 Gpa, respectively. Figure 7b illustrates the undeformed and deformed shapes of the annular plate, as discretized with SHB20 elements, in the case of fully metallic material. Figure 8 reports the load-out-of-plane vertical deflection curves at the outer point A of the annular plate as obtained with the SHB elements, along with the reference solutions taken from [15]. One can observe that the SHB elements perform very well with respect to the reference solutions for all considered values of exponent n. Similar to the previous benchmark problem, the same in-plane mesh discretization as in [15] with only five integration points through the thickness has been adopted by the proposed SHB elements for this nonlinear test, while ten integration points have been considered in [15].  Figure 7b illustrates the undeformed and deformed shapes of the annular plate, as discretized with SHB20 elements, in the case of fully metallic material. Figure 8 reports the load-out-of-plane vertical deflection curves at the outer point A of the annular plate as obtained with the SHB elements, along with the reference solutions taken from [15]. One can observe that the SHB elements perform very well with respect to the reference solutions for all considered values of exponent n. Similar to the previous benchmark problem, the same in-plane mesh discretization as in [15] with only five integration points through the thickness has been adopted by the proposed SHB elements for this nonlinear test, while ten integration points have been considered in [15].    Figure 9b illustrates the undeformed and deformed shapes of the square plate, as discretized with SHB20 elements, in the case of fully metallic material. The pressure-displacement curves for the SHB elements (where the displacement is computed at the center of the plate), along with the reference solutions taken from [15], are all depicted in Figure 10.

Clamped Square Plate under Pressure
The results obtained with the SHB elements, by adopting only five integration points in the thickness direction and the same in-plane mesh discretization as in [15], are in excellent agreement with the reference solutions that required ten through-thickness integration points.

Clamped Square Plate under Pressure
Figure 9a illustrates a fully clamped square plate, which is loaded by a uniformly distributed pressure. The length and thickness of the square plate are L = 1000 mm and t = 2 mm, respectively. The Poisson ratio is ν = 0.3, while the Young modulus of the metal and ceramic constituents are E m = 2 × 10 5 Mpa and E c = 3.8 × 10 5 MPa, respectively. Considering the problem symmetry, a quarter of the plate is discretized. Figure 9b illustrates the undeformed and deformed shapes of the square plate, as discretized with SHB20 elements, in the case of fully metallic material. The pressure-displacement curves for the SHB elements (where the displacement is computed at the center of the plate), along with the reference solutions taken from [15], are all depicted in Figure 10. The results obtained with the SHB elements, by adopting only five integration points in the thickness direction and the same in-plane mesh discretization as in [15], are in excellent agreement with the reference solutions that required ten through-thickness integration points. of the plate is discretized. Figure 9b illustrates the undeformed and deformed shapes of the square plate, as discretized with SHB20 elements, in the case of fully metallic material. The pressure-displacement curves for the SHB elements (where the displacement is computed at the center of the plate), along with the reference solutions taken from [15], are all depicted in Figure 10. The results obtained with the SHB elements, by adopting only five integration points in the thickness direction and the same in-plane mesh discretization as in [15], are in excellent agreement with the reference solutions that required ten through-thickness integration points.  Figure 11a shows a hinged cylindrical roof subjected to a concentrated force at its center. Two types of roofs are considered, thick and thin, with thicknesses t = 12.7 mm and t = 6.35 mm, respectively. Because this nonlinear benchmark test involves geometric-type instabilities (limit-point buckling), the Riks path-following method is used to follow the load-displacement curves beyond the limit points. The Poisson ratio of the cylindrical roof is  Figure 11b illustrates the undeformed and deformed shapes of the hinged cylindrical roof, as discretized with SHB20 elements, in the case of fully metallic material. The load-vertical displacement curves at the central point A of the thick and thin hinged cylindrical roofs are shown in Figures 12 and 13, and compared with the reference solutions taken from [30]. From these figures, it can be seen that the results obtained with the proposed quadratic SHB elements are in good agreement with the reference solutions for the different values of exponent n, corresponding to different volume fractions (from fully metal to fully ceramic). More specifically, the snap-through and snap-back phenomena, which are typically exhibited in such limit-point buckling problems, are very well reproduced by the proposed SHB elements. Note that, for the thick roof (i.e., t = 12.7 mm), the converged solutions in Figure 12 are obtained by using a mesh of (8 × 8 × 2) × 1 in the case of prismatic SHB15 elements, and a mesh of 8 × 8 × 1 with hexahedral SHB20 elements. As to the thin roof (i.e., t = 6.35 mm), finer meshes of (16 × 16 × 2) × 1 for the prismatic SHB15 elements, and 16 × 16 × 1 for the hexahedral SHB20 elements are required to obtain converged results (see Figure 13). These mesh refinements are similar to those used by Sze et al. [29] for the thick and thin roof in the case of an isotropic material as well as for multilayered composite materials.    Figure 11a shows a hinged cylindrical roof subjected to a concentrated force at its center. Two types of roofs are considered, thick and thin, with thicknesses t = 12.7 mm and t = 6.35 mm, respectively. Because this nonlinear benchmark test involves geometric-type instabilities (limit-point buckling), the Riks path-following method is used to follow the load-displacement curves beyond the limit points. The Poisson ratio of the cylindrical roof is ν = 0.3, while the Young modulus of the metal and ceramic constituents are E m = 70 × 10 3 Mpa and E c = 151 × 10 3 Mpa, respectively. Owing to the symmetry, only one quarter of the cylindrical roof is modeled. Figure 11b illustrates the undeformed and deformed shapes of the hinged cylindrical roof, as discretized with SHB20 elements, in the case of fully metallic material. The load-vertical displacement curves at the central point A of the thick and thin hinged cylindrical roofs are shown in Figures 12 and 13, and compared with the reference solutions taken from [30]. From these figures, it can be seen that the results obtained with the proposed quadratic SHB elements are in good agreement with the reference solutions for the different values of exponent n, corresponding to different volume fractions (from fully metal to fully ceramic). More specifically, the snap-through and snap-back phenomena, which are typically exhibited in such limit-point buckling problems, are very well reproduced by the proposed SHB elements. Note that, for the thick roof (i.e., t = 12.7 mm), the converged solutions in Figure 12 are obtained by using a mesh of (8 × 8 × 2) × 1 in the case of prismatic SHB15 elements, and a mesh of 8 × 8 × 1 with hexahedral SHB20 elements. As to the thin roof (i.e., t = 6.35 mm), finer meshes of (16 × 16 × 2) × 1 for the prismatic SHB15 elements, and 16 × 16 × 1 for the hexahedral SHB20 elements are required to obtain converged results (see Figure 13). These mesh refinements are similar to those used by Sze et al. [29] for the thick and thin roof in the case of an isotropic material as well as for multilayered composite materials.

Pull-Out of an Open-Ended Cylinder
The well-known pull-out test for an open-ended cylinder is considered in this section. As illustrated in Figure 14a, the cylinder is pulled by two opposite radial forces, which results in the deformed shape shown in Figure 14b. The isotropic material case as well as the laminated composite material case have been considered by many authors in the literature (see, e.g., [29,30,32]). The Poisson ratio of the cylinder is ν = 0.3, while the Young modulus of the metal and ceramic constituents are E m = 0.7 × 10 9 Pa and E c = 1.51 × 10 9 Pa, respectively. Owing to the symmetry of the problem, only one eighth of the cylinder is modeled. The force-radial displacement curves at points A, B and C (as depicted in Figure 14a), obtained with the SHB elements, are shown in Figures 15-17, respectively, along with the reference solutions taken from [13]. It can be observed that the developed SHB elements successfully pass this benchmark test as compared to the reference solutions. More specifically, the transition zone in the load-radial displacement curves, which is marked by the snap-through point, is well reproduced by both prismatic and hexahedral SHB elements for the various values of the power-law exponent n. Note that the converged solutions in Figures 15-17 are obtained with the proposed elements by using only five integration points in the thickness direction, and meshes of (24 × 36 × 2) × 1 and 12 × 18 × 1 in the case of the prismatic SHB15 element and hexahedral SHB20 element, respectively. Hence, the required meshes for convergence are coarser than those used by Sze et al. [29] in the case of an isotropic material.

Pull-Out of an Open-Ended Cylinder
The well-known pull-out test for an open-ended cylinder is considered in this section. As illustrated in Figure 14a, the cylinder is pulled by two opposite radial forces, which results in the deformed shape shown in Figure 14b. The isotropic material case as well as the laminated composite material case have been considered by many authors in the literature (see, e.g., [29,30,32] Figure 14a), obtained with the SHB elements, are shown in Figures  15-17, respectively, along with the reference solutions taken from [13]. It can be observed that the developed SHB elements successfully pass this benchmark test as compared to the reference solutions. More specifically, the transition zone in the load-radial displacement curves, which is marked by the snap-through point, is well reproduced by both prismatic and hexahedral SHB elements for the various values of the power-law exponent n. Note that the converged solutions in Figures 15-17 are obtained with the proposed elements by using only five integration points in the thickness direction, and meshes of (24 × 36 × 2) × 1 and 12 × 18 × 1 in the case of the prismatic SHB15 element and hexahedral SHB20 element, respectively. Hence, the required meshes for convergence are coarser than those used by Sze et al. [29] in the case of an isotropic material.

Pull-Out of an Open-Ended Cylinder
The well-known pull-out test for an open-ended cylinder is considered in this section. As illustrated in Figure 14a, the cylinder is pulled by two opposite radial forces, which results in the deformed shape shown in Figure 14b. The isotropic material case as well as the laminated composite material case have been considered by many authors in the literature (see, e.g., [29,30,32] Figure 14a), obtained with the SHB elements, are shown in Figures  15-17, respectively, along with the reference solutions taken from [13]. It can be observed that the developed SHB elements successfully pass this benchmark test as compared to the reference solutions. More specifically, the transition zone in the load-radial displacement curves, which is marked by the snap-through point, is well reproduced by both prismatic and hexahedral SHB elements for the various values of the power-law exponent n. Note that the converged solutions in Figures 15-17 are obtained with the proposed elements by using only five integration points in the thickness direction, and meshes of (24 × 36 × 2) × 1 and 12 × 18 × 1 in the case of the prismatic SHB15 element and hexahedral SHB20 element, respectively. Hence, the required meshes for convergence are coarser than those used by Sze et al. [29] in the case of an isotropic material.

Pinched Hemispherical Shell
It is worth noting that although the performance of the prismatic SHB15 element is similar to that of the hexahedral SHB20 element, as demonstrated in the above nonlinear benchmark problems, the main motivation in developing the prismatic solid-shell element is to use it for the mesh discretization of complex geometries. Indeed, it is well-known that complex geometries cannot be discretized with only hexahedral elements, and require either an irregular mesh with prismatic elements, or a mixture based on a combination of prismatic and hexahedral elements. In this section, a hemispherical shell is loaded by alternating radial forces as shown in Figure 18a. Note that this benchmark problem has been considered in the literature for an isotropic material as well as a laminated composite material (see, e.g., [30]), while the case of FGMs has not been considered yet. Consequently, only the simulation results obtained with the proposed SHB elements corresponding to the fully metallic shell can be compared to the reference solution taken from [30].

Pinched Hemispherical Shell
It is worth noting that although the performance of the prismatic SHB15 element is similar to that of the hexahedral SHB20 element, as demonstrated in the above nonlinear benchmark problems, the main motivation in developing the prismatic solid-shell element is to use it for the mesh discretization of complex geometries. Indeed, it is well-known that complex geometries cannot be discretized with only hexahedral elements, and require either an irregular mesh with prismatic elements, or a mixture based on a combination of prismatic and hexahedral elements. In this section, a hemispherical shell is loaded by alternating radial forces as shown in Figure 18a. Note that this benchmark problem has been considered in the literature for an isotropic material as well as a laminated composite material (see, e.g., [30]), while the case of FGMs has not been considered yet. Consequently, only the simulation results obtained with the proposed SHB elements corresponding to the fully metallic shell can be compared to the reference solution taken from [30].

Pinched Hemispherical Shell
It is worth noting that although the performance of the prismatic SHB15 element is similar to that of the hexahedral SHB20 element, as demonstrated in the above nonlinear benchmark problems, the main motivation in developing the prismatic solid-shell element is to use it for the mesh discretization of complex geometries. Indeed, it is well-known that complex geometries cannot be discretized with only hexahedral elements, and require either an irregular mesh with prismatic elements, or a mixture based on a combination of prismatic and hexahedral elements. In this section, a hemispherical shell is loaded by alternating radial forces as shown in Figure 18a. Note that this benchmark problem has been considered in the literature for an isotropic material as well as a laminated composite material (see, e.g., [30]), while the case of FGMs has not been considered yet. Consequently, only the simulation results obtained with the proposed SHB elements corresponding to the fully metallic shell can be compared to the reference solution taken from [30]. symmetry, a quarter of the structure is discretized. The hemispherical shell is discretized with a mixture of prismatic and hexahedral elements, which consists of 90 SHB15 elements located at the top of the hemisphere (far from the load points, see Figure 18b) and 110 SHB20 elements for the remaining area. The simulation results in terms of force-radial deflection at point A, for various values of the power-law exponent n, are plotted in Figure 19. This figure shows that the results corresponding to a fully metallic shell (i.e., n → ∞ ), obtained by the combination of prismatic and hexahedral SHB elements, are in excellent agreement with those provided in [30] for an isotropic shell. Note that an equivalent in-plane mesh discretization has been used in [30], where a fully integrated shell element with several integration points has been considered. Moreover, Figure 19 reveals that for all values of the exponent n, the simulated load-radial deflection curves lie between that of the fully ceramic shell and that of the fully metal shell, which is consistent with the numerical results found in the previous benchmark problems. The radius and thickness of the hemispherical shell are R = 10 m and t = 0.04 m, respectively. The Poisson ratio of the hemispherical shell is ν = 0.3, while the Young modulus of the metal and ceramic constituents are E m = 6.825 × 10 7 Pa and E c = 1.46 × 10 8 Pa, respectively. Due to the symmetry, a quarter of the structure is discretized. The hemispherical shell is discretized with a mixture of prismatic and hexahedral elements, which consists of 90 SHB15 elements located at the top of the hemisphere (far from the load points, see Figure 18b) and 110 SHB20 elements for the remaining area.
The simulation results in terms of force-radial deflection at point A, for various values of the power-law exponent n, are plotted in Figure 19. This figure shows that the results corresponding to a fully metallic shell (i.e., n → ∞ ), obtained by the combination of prismatic and hexahedral SHB elements, are in excellent agreement with those provided in [30] for an isotropic shell. Note that an equivalent in-plane mesh discretization has been used in [30], where a fully integrated shell element with several integration points has been considered. Moreover, Figure 19 reveals that for all values of the exponent n, the simulated load-radial deflection curves lie between that of the fully ceramic shell and that of the fully metal shell, which is consistent with the numerical results found in the previous benchmark problems. symmetry, a quarter of the structure is discretized. The hemispherical shell is discretized with a mixture of prismatic and hexahedral elements, which consists of 90 SHB15 elements located at the top of the hemisphere (far from the load points, see Figure 18b) and 110 SHB20 elements for the remaining area. The simulation results in terms of force-radial deflection at point A, for various values of the power-law exponent n, are plotted in Figure 19. This figure shows that the results corresponding to a fully metallic shell (i.e., n → ∞ ), obtained by the combination of prismatic and hexahedral SHB elements, are in excellent agreement with those provided in [30] for an isotropic shell. Note that an equivalent in-plane mesh discretization has been used in [30], where a fully integrated shell element with several integration points has been considered. Moreover, Figure 19 reveals that for all values of the exponent n, the simulated load-radial deflection curves lie between that of the fully ceramic shell and that of the fully metal shell, which is consistent with the numerical results found in the previous benchmark problems. Metal from [30] Metal with SHB elements n = 5 with SHB elements n = 2 with SHB elements n = 1 with SHB elements n = 0.5 with SHB elements n = 0.2 with SHB elements Ceramic with SHB elements P Radial deflection at point A (m) Figure 19. Load-displacement curves at point A for the pinched hemispherical shell, obtained with a mixture of prismatic SHB15 and hexahedral SHB20 elements.

Conclusions
In this work, quadratic prismatic and hexahedral solid-shell SHB elements have been proposed for the 3D modeling of thin FGM structures. The formulation of the SHB elements adopts the in-plane reduced-integration technique along with the assumed-strain method to alleviate various locking phenomena. A local (element) frame has been defined for each element, in which the thickness direction is specified. In this local frame, elastic properties of the thin structure are assumed to vary gradually through the thickness according to a power-law volume fraction distribution. The resulting formulations are implemented into the finite element software ABAQUS/Standard in the framework of large displacements and rotations. A series of selective and representative benchmark problems, involving FGM thin structures, has been performed to evaluate the performance of the SHB elements in geometrically nonlinear analysis. The results obtained with the SHB elements have been compared with reference solutions. Note that the state-of-the-art ABAQUS solid and shell elements have not been considered in the simulations, because these elements do not allow modeling of FGM behavior with only a single layer of elements. Overall, the numerical results obtained with the SHB elements showed excellent agreement with the available reference solutions. More specifically, the load-displacement curves for each benchmark test lie between that of the fully ceramic and fully metal behavior, which is consistent with the power-law distribution of the Young modulus in the thickness direction of the FGM plates. This good performance of the proposed SHB elements only requires a few integration points in the thickness direction (i.e., only five integration points), as compared to the number of integration points used in the literature to model thin FGM structures. Furthermore, it has been shown that the prismatic SHB15 element can be naturally combined with the hexahedral SHB20 element, within the same simulation, to help discretize complex geometries. Overall, the proposed SHB elements showed good capabilities in 3D modeling of thin FGM structures with only a single layer of elements and few integration points, which is not the case of traditional solid and shell elements.