A Stress Approach Model of Functionally Graded Shells

-is paper presents a model called SAM-FG (stress approach model of functionally graded shells) for linear elastic, thin, and moderately thick shells made of functionally graded materials. -e model is an extension of the SAM-Hmodel, originally created for homogeneous shells. Assuming that the material is orthotropic and that one of its orthotropic directions is the thickness direction, the extension consists in considering that the 3D compliance tensormay depend on the thickness coordinate.-emodel starts with a tunable polynomial approximation of the 3D stress field that contains the same generalized forces as SAM-H. -is stress approximation verifies the 3D equilibrium equations and the stress boundary conditions at the faces of the shell. As in SAMH, 5 generalized displacements appear in SAM-FG. By applying the Hellinger–Reissner functional and Reissner’s variational method, the generalized forces, strains, and equations in SAM-FG turn out to be the same as in SAM-H, except for the generalized constitutive equations. To prove the accuracy of the model, SAM-FG is first applied to a simply supported, functionally graded plate and its results are compared to other models. To validate the model for shell-like structures, SAM-FG results are compared to those obtained with solid finite element calculations for three case studies of structures subjected to an internal pressure. -e first one deals with a hollow sphere made of an isotropic functionally graded material. -e second case considers a hollow cylinder made of an orthotropic functionally graded material. In the last case, a catenoid with an isotropic functionally graded material is studied. In all cases, the mean displacements are correctly predicted, even if the main purpose of the SAM-FG model is not to calculate these fields accurately.-e stress field approximations are very accurate, and since the implementation of the shell model in a finite element code would imply 5 degrees of freedom per node, SAM-FG is a good alternative to solid finite element calculations for the structural analysis of functionally graded shells with a reasonable computational cost.


Introduction
e aerospace industry has been looking for lighter materials that can withstand high stress levels and temperature gradients. For aerospace vessels, in certain conditions, a temperature difference of about 1000°C may be expected between the outer surface and the inner parts [1]. Homogeneous materials are not the best solution for these cases. An important improvement in material research came in 1987 by a group of Japanese scientists who created a functionally graded material (FGM) that served to increase thermal barriers between the inside and outside of a thin part of the plane [2]. In this heterogeneous material, a smooth transition from one homogeneous material to other helped to reduce thermal stresses and performed better than homogeneous materials with similar constituents. From that point forward, numerous scientific groups studied this concept to find more ways to build them [3][4][5] and exploit their benefits with new applications [6][7][8][9]. e design of structural parts made of FGMs has implied the modification and adaptation of beam, plate, and shell models originally developed for structures made of homogeneous materials. e classical Euler-Bernoulli and Timoshenko beam theories have been revisited in order to involve FGMs.
Chakraborty et al. [10] expanded the Timoshenko theory with FGM behavior and thermal deformation effects. More recently, Wang and Li [11] adapted the Levinson beam theory to analyze the natural frequencies of FGM beams and obtained accurate results. Xu and Meng [12] proposed a beam model that calculates bending, buckling, and free vibrations for several regular polygonal cross sections considering an FGM whose elastic properties vary from the center of the section to its boundary following a power law of the space coordinates.
Continuing with plate models, Chi and Chung [13,14] considered a moderately thick rectangular plate made of an FGM and adopted a Kirchhoff-Love-like theory: a displacement approach with 3 kinematic fields. Normal out-ofplane stresses were neglected. e model results (displacements and in-plane stresses) were compared with solid finite element (SFE) calculations at the middle surface of the plate. A good agreement with SFE was obtained. In [15], Altenbach and Eremeyev created a model based on Zhilin's theory and compared it with an equivalent of Kirchhoff and Mindlin models for case studies involving linear and power-law distributions for material properties. Even though the results were not validated in that paper, it was established that the different FGM approaches yield different results caused by the thickness coordinate behavior. Carrera et al. [16] created a model for FGM plates based on the classical principle of virtual displacements expanding taking into account both equivalent single-layer (ESL) and layerwise (LW) models proposing linear to fourth-grade polynomial approximations for the displacement field, as well as assuming that the stiffness matrix varies along the thickness coordinate. To validate its results, both closed-form and finite element solutions were compared to a 3D model, which yielded satisfactory results and an improvement on classical theories in two case studies. Nguyen et al. [17] proposed a model based on the first-order shear deformation theory (FSDT) applied to functionally graded material plates using shear correction factors (SCFs) to calculate accurately shear forces and strains. In another study [18], Nguyen et al. also collaborated in creating a formulation that uses three variables for shear deformation that eliminates the necessity of SCF. Other models equivalent to the FSDT are those developed by Zenkour [19] and Touratier [20] that include 5 generalized displacements and do not need SCF. Another interesting model was developed by SiddaRedddy et al. [21]; it is based on a higher-order shear deformation theory (HSDT) using 9 degrees of freedom and proved to yield accurate results for FGM plates. Other forms of calculation can be used to analyze FGM plates, like the sampling surfaces method (SaS) created by Kulikov in 2001 [22], which has been studied and improved in [23,24].
is method consists in creating virtual surfaces uniformly distributed through the thickness of the plate and considering the displacements at each surface as generalized displacements which help to build the 3D displacement field.
Shell-like structures made of FGM have also been the object of the development of several models. Numerous shell models have been published for specific geometries. Reddy and Chin [25] proposed a model based on an FSDTapproach to study the dynamic thermoelastic response of a functionally graded plate and a one-dimensional axisymmetric cylinder. Sarathchandra et al. [26] designed a model exclusively developed for cylindrical shells which evaluates stresses caused by mechanical and thermal loads. Also, Poorna and Ram [27] studied vibration modes for cylindrical and spherical shells using piezoelectric materials. Other interesting models for the analysis of vibrations in functionally graded shells are those developed by Kiani et al. [28,29] and Arefi et al. [30]. In [31], Kulikov and Plotnikova applied his SaS method to the study of static shells and validated the method for 7 reference surfaces by comparing the model results with those of SFE in the case of a homogeneous cylinder. Other shell models involving thermal, electric, or magnetic loads are worth to mention; for example, Zhu et al. [32] studied the effect of mechanical, thermal, and electric loads on a piezoelectric cylindrical shell reinforced with nanotubes distributed in various forms. Dini et al. [33] created a model that calculates stresses given by mechanical forces as well as thermal and magnetic fields in rotating rings made of an FGM with a layer of metal on the inside and ceramic on the outside.
Herein, the authors propose a stress approach model of functionally graded shells called SAM-FG. is model is inspired by the SAM-H model of homogeneous shells developed by Domínguez-Alvarado and Díaz-Díaz in [34], which has proved to approximate accurately displacements and stresses (particularly, out-of-plane stresses) in thin and moderately thick shells. e main features of SAM-H are the verification of the 3D equilibrium equations, the fulfillment of the stress boundary conditions at the faces of the shell, and the inclusion of Poisson's effect of out-of-plane normal stresses on in-plane strains. e present paper is structured in 4 sections. e first section details the notations adopted. e next section recalls briefly the assumptions adopted in SAM-H and how the equations of the model are obtained. e third part details the equations in the SAM-FG model. In the fourth section, the numerical results of SAM-FG are compared to those of SFE results in order to test the accuracy of the model.

Notations and Hellinger-Reissner Functional
First, let us define the notations adopted: (i) Boldface characters define vectors, matrices, and tensors; Greek subscripts take the values 1 and 2; if not otherwise specified, Latin subscripts take the values 1, 2 and 3; tensor product and simple and double contractions of tensors are denoted by: " ⊗ ," "·," and ":," respectively. (ii) Ω is the 3D domain in which lies the solid; h is the thickness of the shell; (ξ 1 , ξ 2 , ξ 3 ) are the orthogonal curvilinear coordinates considered and e 1 , e 2 , and e 3 are their related unit vectors; ω is the middle surface located at ξ 3 � 0; ξ 1 and ξ 2 are the coordinates along the lines of curvatures at ω; the Lamé coefficients at the middle surface are a 1 and a 2 ; the principal curvatures are κ 1 and κ 2 .
Additionally, the following assumptions are adopted: (1) In all calculations, the terms (ξ 3 /R α ) and (h/R α ) are considered, but higher-order terms (ξ 3 /R α ) n and (h/R α ) n (n ≥ 2) are neglected (2) e curvatures vary smoothly along the lines of curvature and their derivatives with respect to ξ α are neglected (3) For simplicity sake, the body loads are neglected (4) No 3D displacement constraints are applied on the inner and outer faces of the shell Assuming small displacements and strains, the Hellinger-Reissner (HR) functional [35] is expressed as follows: where (i) u * is a 3D displacement vector. e superindex * indicates that the field is not necessarily the solution of the problem. (ii) σ * is a second-order stress tensor. (iii) zΩ u and zΩ s are the boundaries where displacement and stress vectors are imposed, respectively. (iv) ε(u * ) is the Green-Cauchy 3D strain tensor.
(v) S is the 3D, 4th-order compliance tensor, and its components are S ijkl . (vi) n is the outward-pointing normal vector at the boundary of Ω. (vii) u g and s g are the imposed displacement and stress vectors at the boundaries.

Recalling SAM-H Equations
e development of SAM-FG follows almost the same procedure as SAM-H. e difference between SAM-H and SAM-FG lies in the stress approximation and constitutive equations. A thorough explanation of the development of SAM-H equations can be found in [34]. Nevertheless, a summary of these equations is shown in this section.

Stress Approximation and Generalized Forces.
e stress components are approximated by ξ 3 -polynomials. e following basis of 4th-degree, ξ 3 -polynomials is considered: Mathematical Problems in Engineering 3 \openup5 P 0 ξ 3 � 1, e following stress approximation is obtained: where σ n * ij are the stress coefficients for the approximation, which are linear combinations of the generalized forces and components of the applied stress vector at the faces of the shell. In [34], the complete detail of the expressions of these stress coefficients is shown. e 2D tensors of generalized forces are membrane forces N * , shear forces Q * , and moments M * ; their components are defined as follows:

Generalized Displacements, Strains, and Equilibrium
Equations. Introducing the stress approximation in the HR functional helps to identify the 5 generalized displacements of the model: Also, the following generalized strains are obtained: where By applying the stationarity of the functional with respect to the generalized displacements, the term Ω (σ * : ε(u * ) − f · u * )dΩ in the functional yields five scalar equilibrium equations which can be written in a compact manner as follows: Moreover, this stationarity also yields the edge conditions on the generalized forces.

Generalized Constitutive Equations.
SAM-H assumes that the material is homogeneous and orthotropic and that one of its orthotropic directions is direction 3. e following compliance tensors and constants are defined: 4 Mathematical Problems in Engineering S αβcδ e α ⊗ e β ⊗ e c ⊗ e δ , S αβ33 e α ⊗ e β , All the abovementioned compliances are assumed uniform through the thickness of the shell. e following engineering vectors of generalized forces F and strains ε are defined: where superscript t denotes the transposition operation. By introducing the stress approximation in the HR functional and by applying the stationarity of the functional with respect to the generalized forces, one obtains the generalized constitutive equations of the model:

Stress
Approximation. SAM-FG uses the same generalized forces as SAM-H, but due to the assumed material heterogeneity, its stress approximation is different. In order to determine a correct stress approximation in SAM-FG, let us first deduce the form of the 3D in-plane strains ε 3 D αβ . As applied in the classical laminate theory (CLT), in a moderately thick FG shell, far enough from the edges, the 3D inplane displacements are virtually first degree polynomials of the thickness coordinate. For this reason, the vector of 3D in-plane strains take the following expression: where v ε,0 and v ε,1 are nondimensional vector fields that depend on the in-plane coordinates. Let us define the 3 × 3 plane-stress stiffness matrix R 3D (ξ 3 ); assuming a plane-stress state, the vector of in-plane stress components is given by In SAM-FG, the components of the 3D in-plane stiffness matrix are approximated across the thickness (for ξ 3 ∈ [− (h/2), (h/2)]) by a polynomial fit of order r. In this manner, where R 3D,i are constant matrices and P i (1 ≤ i ≤ n) is an orthogonal basis of polynomials ensuring the following condition: e first 5 polynomials of this basis are shown in equation (6); in the other polynomials, "1" is chosen as the leading coefficient, for example, . (24) Equations (21) and (22) yield Let us define the following vectors of generalized forces and moments: where superscript t denotes the transposition operator.
Using the stress approximation in equation (25) and owing to equations (8) and (10), one can obtain a system of linear equations that relates the 6 components of vector v NM with the 6 strains appearing in vectors v ε,0 and v ε,1 . After solving for the strain vectors, one obtains where C N,0 , C N,1 , C M,0 , and C M,1 are 3 × 4 compliance matrices. By introducing the equation above in equation (25), one obtains the approximation of in-plane stresses of the SAM-FG model: Mathematical Problems in Engineering ese in-plane stresses are ξ 3 -polynomials of degree r + 1: where σ n * αβ are stress coefficients that, for a given FGM, can be expressed as linear algebraic combinations of only 8 generalized fields: the generalized membrane forces and moments given in equations (8) and (10). e 3D equilibrium equation div σ * � 0 is used to deduce successively the degree of the polynomial approximation across the thickness for the shear stresses σ * α3 and normal stress σ * 33 . According to equations (2)-(4), the polynomial degrees of σ * α3 and σ * 33 turn out to be r + 3 and r + 4, respectively. e r + 4 coefficients of the σ * α3 polynomial are obtained by making use of a set of r + 4 linear equations. e first two equations are deduced from the 3D boundary condition at the faces of the shell: e third equation comes from the definition of the generalized shear forces in equation (9). e other r + 1 equations are obtained by following the next steps: is component is written as a series expansion of polynomials P i . For i > r + 2, the coefficients of these polynomials are zero or negligible because of being multiplied by second-order terms: (iii) e searched equations come from solving for the coefficients of polynomials P i (2 ≤ i ≤ r + 2) in the aforementioned series.
In this manner, the proposed expression for the out-ofplane shear stresses is given by where σ n * α3 are stress coefficients which can be written as linear combinations of τ + α , τ − α , Q α and the derivatives of N αβ and M αβ . For the determination of the coefficients of the polynomial approximation of the normal stress σ * 33 , a similar method is applied, but this time, the boundary conditions σ * 33 (ξ 3 � (h/2)) � σ + and σ * 33 (ξ 3 � − (h/2)) � − σ − are used and the third component of the vectorial equation div σ * � 0 is selected. e proposed expression for the out-of-plane normal stress is given by where the stress coefficients σ n * 33 can be written as linear combinations of the applied stresses at the faces of the shell and the derivatives of the generalized forces and moments defined in equations (8)-(10).

Generalized Displacements, Strains, and Equations of SAM-FG.
Following the same methodology of SAM-H, the stress approximation is introduced into the HR functional shown in equation (5) in order to identify the generalized displacements. It turns out that the same 5 generalized displacements shown in equation (11) are obtained for SAM-FG. Also, the same generalized strains as those in SAM-H (see equation (12)) are deduced. By applying the stationarity of the functional with respect to the generalized displacements, the term Ω (σ * : ε(u * ) − f · u * )dΩ yields the same equilibrium equations as those obtained for SAM-H and shown in equations (14)- (16).
Assuming that the material is orthotropic and that one of its principal directions is normal to the shell surface, the integral of the volumetric elastic energy appearing in the HR functional is expressed as follows: where w s * e , w n * e , w Q * e , and w c e are the surface densities of elastic energy defined by (summation over α, β, c, and δ): It is worth defining the total surface density of elastic energy w * e � (w s * e + w n * e + w Q * e + w c * e ). Let us recall that the energies w n * e and w c * e are usually neglected in plate and shell models but not in SAM-FG.
Owing to the polynomial expansions in equations (29)-(31), the surface densities of energy are expressed as follows: e in-plane stress coefficients σ m αβ are expressed as linear combinations of the generalized membrane forces and moments. e coefficients σ m α3 and σ n 33 originally involve the divergences of the generalized forces N , Q, and M . By making use of the equilibrium equations (14)- (16) and only for the approximation of the surface densities of elastic energy, these divergences may be approximated in terms of algebraic expressions of the generalized forces and applied stresses at the faces of the shell. us, the sum of surface densities of energy w * e � (w s * e + w n * e + w Q * e + w c * e ) appearing in equation (32) can be expressed as a linear combination of the following terms: where 1 ≤ p, q ≤ 8, 0 ≤ m, n ≤ 4, F p are the components of vector F defined in equation (18) and represent the generalized membrane forces and moments, Q α are the generalized shear forces, and A mn represents the integrals defined in equation (34). If possible, the integrals defined in the terms A mn are determined analytically; otherwise, they are calculated numerically. In this manner, w * e takes the form: where αβ , and D 2− αβ are compliance coefficients that are deduced from the integrals in equation (34) and the coefficients multiplying the generalized forces in the expressions of the stress coefficients σ m ij ; Z is a term that depends only on the applied stresses at the faces of the shell. After applying the stationarity of the HR functional, the term in equation (5) yields the following constitutive equations: where ε p are the components of vector ε defined in equation (18). ese constitutive equations can be written in a compact form as follows: e components of matrices C 1 , D 1 , D 2+ , and D 2− and those of vectors C 2+ , C 2− , and C 3 can easily be identified with the compliance coefficients used in equation (38).

Practical Use of SAM-FG.
e principal orthotropic directions of the material are L, T, and N; the unit vectors associated with these directions are as follows: e L � cos θ ξ 3 e 1 + sin θ ξ 3 e 2 , e T � − sin θ ξ 3 e 1 + cos θ ξ 3 e 2 , where θ is the angle between direction L and direction 1, and this angle may vary across the thickness. One may provide directly the values of 3D stiffness and compliance components in the orthotropic directions, but in practice, it is easier to identify Young's moduli, Poisson's ratios, and shear moduli in the orthotropic directions. ese properties depend on ξ 3 . Some authors, like Altenbach and Eremeyev in [15], analyze the particular case when these properties follow a power law of the thickness coordinate:

Mathematical Problems in Engineering
where ψ represents the elastic properties in the orthotropic directions; superscripts ext and int denote the value of the property at the outer and inner surfaces, respectively; n ≥ 0 is a constant which modifies the power law. In Figure 1, the example of Young's modulus is plotted against the normalized thickness coordinate (ξ 3 /h) for n � 0.2, 0.5, 1, 2, and 5. In this example, the chosen Young's moduli at the inner and outer surfaces were 100 GPa and 200 GPa, respectively. Once the functions representing the elastic properties E L , E T , E N , ] LT , ] LN , ] TN , G LT , G LN , G TN , and the angle θ being defined, by making use of rotation operations, one obtains the components of the 3 × 3 plane-stress stiffness matrix R 3D appearing in equation (21) and the components of compliances S , S c , S Q , and S σ in the basis (e 1 , e 2 , e 3 ).
After this, the degree r of the polynomial fit of the components of R 3D (ξ 3 ) is selected and the constant matrices R 3D,i (0 ≤ i ≤ r) appearing in equations (22) and (28) are determined. e components of R 3D,i are required for the calculation of the stress coefficients σ n ij . en, the integrals A mn ijkl in equation (34) can be determined numerically or if possible, analytically. e generalized compliances appearing in the constitutive equations (40) can then be obtained. In practice, it is more useful to express the generalized forces in terms of generalized strains; with a similar method to that used by Domínguez-Alvarado and Díaz-Díaz in [34], one can obtain the alternative form of the constitutive equations: where K and L are 6 × 6 and 2 × 2 stiffness matrices, respectively; k + and k − are 6-dimensional vectors; and L + and L − are 2-dimensional square matrices. e expressions of these matrices and vectors are given by In the previous definition, H 1 In this paper, for a given FGM shell, MATLAB software was first applied to calculate the integrals appearing in the generalized compliances and then to determine the coefficients of K, L, k + , k − , L + , and L − . All these coefficients were input manually in COMSOL Multiphysics (5.3a version), a commercial finite element software which was used to solve SAM-FG equations for the shell considered. e LiveLink module in COMSOL for MATLAB would help to make this task automatically.
Let us finally emphasize that, whatever the degree r of the polynomial fit of the components of the stiffness matrix R 3D (ξ 3 ), SAM-FG uses the same number of generalized displacements, strains, and forces. e polynomial fit is a preprocessing task similar to the determination of the 6 × 6 ABB D stiffness matrix in a laminated plate made up of a variable number of layers. SAM-FG equations (14)-(16) and (43) can be written in the form of a second-order partial differential equation system whose solution is the set of 5 generalized displacements. ese equations were implemented in the PDE (partial differential equation) solver module of COMSOL Multiphysics 5.3a. After solving for the generalized displacements, the generalized forces are determined in a postprocessing stage by making use of the constitutive equations (43). Finally, the polynomial approximation of stresses is performed and the results can be plotted as a color map on a reconstructed 3D solid geometry.

Results and Discussion
In this section, SAM-FG is applied to 4 mechanical problems. In all cases, a 5th-order polynomial fit (r � 5) was enough to accurately approximate the stiffness coefficients. Higher values of r yield virtually the same numerical results. In this sense, the results of SAM-FG shown in this section converge. SAM-FG results are compared to those of solid finite element (SFE), calculations are also performed, and these also converge.

Plate Subjected to a Sinusoidal
where q 0 � 1 MPa and λ � (π/a). is case was selected because several authors have used it to test the accuracy of their models: Reddy and Chin [25], Zenkour [19], Touratier [20], and SiddaRedddy et al. [21]. ese authors have used a displacement approach consisting of more than 5 generalized displacements (degrees of freedom). e model proposed by SiddaRedddy et al. [21] seems to be the most accurate since it uses 9 generalized displacements. In order to have an accuracy reference, SFE in COMSOL Multiphysics software was also applied. A functionally graded isotropic material is considered; Poisson's ratio is 0.3 and constant across the thickness, but Young's modulus follows the power law in equation (42): at the bottom face E int � 70 GPa and at the upper face E ext � 380 GPa. Several values of the exponent n in the power law are considered. Let us define the following dimensionless stresses evaluated at particular points of the structure: SAM-FG equations were solved numerically using the finite element software COMSOL Multiphysics. e selected maximum element size was 17.5 mm for both SFE and SAM-FG. e 3D SFE mesh was refined within the zones where stresses are evaluated.
In Table 1, the results obtained for the dimensionless stresses are summarized for various n values. Using the SFE results, the relative errors of SAM-FG and the model proposed in [21] are also shown in this table. It can be appreciated that, with the exception of a couple of values for σ x and σ y , the errors presented by SAM-FG are below 1%. For σ y , σ xz , and σ yz , nonetheless, SAM-FG shows greater accuracy than all references in almost all the cases. is effect is accentuated in σ xz for n � 5, where the results presented in [21] have an error of 7.82% and the one obtained by SAM-FG is only 0.29%. is proves the benefits of the present model.
Let us plot some of the dimensionless stresses evaluated by SAM-FG across the thickness at the (x, y) locations specified in equation (49). Figures 2 and 3 display the evolution of σ x , σ xy , and σ xz for n � 2 and n � 10, respectively. e results obtained with SFE are also displayed.
In both cases (n � 2 and n � 10), the stresses are accurately evaluated by SAM-FG.

Hollow Sphere Subjected to an Internal Pressure.
In this section, a hollow sphere subjected to an internal pressure is considered. e material is isotropic, but its Young's modulus and Poisson's ratio vary linearly through the thickness. At the inner surface of the sphere, Young's    which can be compared to the average value U 3D 3 of the 3D solution across the thickness of the sphere. e relative error in the calculation of the average value of the radial displacement is defined by In a similar manner, one may define the relative error δσ in the calculation of the maximum von Mises stress determined through the thickness of the sphere.
Let us first analyze the case when ] ext � ] int and η � 0.1. Figure 4 shows the graph of the relative error δU against the heterogeneity ratio (E int /E ext ).
ese figures consider the case of 5 values of Poisson's ratio: − 0.4, − 0.2, 0, 0.2, and 0.4. It is worth mentioning that, owing to material linearity, these curves do not depend on the magnitude of the applied pressure. e absolute value of the relative error is less than 1.5%. Let us now analyze the relative error δσ in the calculation of the maximum von Mises stress determined through the thickness of the sphere. In Figure 5, the relative error δσ is plotted against the (E int /E ext ) ratio. Once again, 5 values of Poisson's ratio were considered. In each case, the relative error δσ is less than 3.6%.
Let us now analyze the effect of the thickness-to-radius ratio η on the accuracy of the model. Four material types are considered: In Figure 6, the relative error δU is plotted against the thickness-to-radius ratio η. As η increases, the absolute value of the relative error increases. For moderately thick shells (η ≤ 0.1), |δU| is smaller than 0.8%, and for thicker shells (η > 0.1), this error remains small ( < 5.3%). erefore, SAM-FG is very accurate in the approximation of the mean 3D displacement of the sphere for all the materials considered. Figure 7 displays the chart of δσ vs. the thickness-to-radius ratio η. e error |δσ| is greater for materials 2 and 4 (] int � 0.2 and ] ext � 0.4) than for the other materials. For moderately thick shells, |δσ| is less than 1%, and for thick shells, |δσ| is less than 5.2%. Let us compare SAM-FG and 3D stress components results for material 1 for the following thickness-to-radius ratios η: 0.01, 0.1, and 0.3. In Figures 8  and 9, the normalized in-plane stress ((σ 11 h)/(pR)) and the normalized radial stress (σ 33 /p) are plotted against the    normalized position through the thickness (ξ 3 /h), respectively. It can be seen that SAM-FG predicts accurately both the in-plane stress and the radial stress. e maximum and minimum values of the radial stress are perfectly calculated by SAM-FG because in this model, the approximated stress σ 33 verifies the boundary conditions at the faces of the shell.

Hollow Cylinder Subjected to an Internal Pressure.
In this section, a hollow cylinder subjected to a 1 MPa internal pressure is considered. e ends of the cylinder are axially blocked. e radius of the middle surface is R 0 � 1 m and the thickness-to-radius ratio is η � 0.1. e material is orthotropic and its principal directions 1, 2, and 3 coincide with the axial, azimuthal, and radial directions, respectively. e elastic moduli at the internal (superscript "int") and external (superscript "ext") surfaces of the shell are shown in Table 2.
Poisson's ratios ] 12 , ] 13 , and ] 23 are assumed to be uniform and equal to 0.3. ese elastic properties vary through the thickness of the cylinder following the power law shown in equation (42). Several values of the exponent n appearing in the power law will be tested in the following paragraph. Owing to symmetries, the 3D equations of the problem are reduced to a 1D problem which was solved, once again, in COMSOL Multiphysics 5.3a software; the solution obtained is called the 3D solution. For the SAM-FG model, after applying the symmetries, one obtains a set of algebraic equations. In Figures 10-12, the stress components σ 11 (axial stress), σ 22 (hoop stress), and σ 33 (radial normal stress) are plotted against the normalized position (ξ 3 /h) for n � (1/3), n � 1, and n � 3, respectively. In all cases, the normal stress σ 33 calculated by SAM-FG verifies the boundary conditions (σ 33 � − 1 MPa at (ξ 3 /h) � − 0.5 and σ 33 � 0 at (ξ 3 /h) � 0.5). For n � (1/3) and n � 3, the normal stresses σ 11 and σ 22 do not vary linearly across the thickness and SAM-FG predicts accurately this variation. In all cases, the stress approximation is very accurate.   Let us now compare the generalized displacement U SAM 3 provided by SAM-FG with the average value U 3D 3 of the 3D solution across the thickness of the cylinder. In Figure 13, the graphs of U SAM 3 and U 3D 3 are plotted against the exponent n appearing in the power law. e evolution of the curve of U 3D 3 is accurately approximated by SAM-FG. In this figure, the relative error δU in the calculation of the average value of the radial displacement is also plotted. e absolute value of the relative error is less than 0.4%; this proves the great accuracy of SAM-FG.

Catenoid Subjected to an Internal Pressure.
e last case considers a 0.1 m thick catenoid. It is a solid of revolution and its middle surface is given by the following parametric equations in the r − z plane (see Figure 14):  Internal  50  50  8  20  10  10  External  100  100  12  30  15   Since the principal curvatures κ 1 and κ 2 are not uniform, the thickness-to-radius ratio η is not constant. e maximum ratio η max and average ratio η for the catenoid are as follows: e material is isotropic and its Young's modulus E and Poisson's ratio ] vary linearly across the thickness in such a manner that these properties are E � 200 GPa and ] � 0.3 at the internal face and E � 550 GPa and ] � 0.31 at the external face. e boundary conditions of the 2D axisymmetric problem are shown in Figure 14: a roller is considered at the basis of the catenoid and an internal pressure p � 1 MPa is applied at the internal face.
For this case, SAM-FG equations cannot be reduced to a set of algebraic equations and a finite element resolution of the generalized equations is performed by making use of COMSOL Multiphysics 5.3a. SAM-FG results were compared with those obtained using solid finite elements for the axisymmetric problem. In Figure 15, the used meshings for the calculations are shown. A 32.5 mm maximum size of elements and cubic Lagrange shape functions were applied in both calculation methods. In the SFE model, the meshing in lines A and B shown in Figure 14 was refined so as to plot accurately the 3D stress components along these lines. e maximum element size along these lines was 1 mm.
In order to evaluate the accuracy of SAM-FG in the evaluation of stresses through the thickness of the structure, the stress approximation was plotted in the cross section (r − z plane) and also along lines A and B shown in Figure 14. In Figure 16, the stress components σ 11 and σ 22 are plotted across the section of the catenoid. SAM-FG stress color maps are very similar to those obtained by SFE. SAM-FG evaluates accurately both σ 11 and σ 22 stresses: the relative errors in the maximum evaluation of σ 11 and σ 22 are 0.47% and 0.27%, respectively. In Figure 17, the normal stress σ 33 and the equivalent stress σ eq (von Mises stress) are plotted across the section. Both stresses are evaluated with great accuracy by SAM-FG.
In Figure 18, the stress components σ 11 , σ 22 , σ 33 , and σ eq are plotted against the normalized position (ξ 3 /h) along line A for both SAM-FG and SFE. e internal and external faces of the shell are located at (ξ 3 /h) � − 0.5 and 0.5, respectively; as expected, the normal stress σ 33 evaluated by SAM-FG verifies the boundary conditions at these faces. One can observe that all the stress components are accurately evaluated along line A. In Figure 19, the evolution of the stress components σ 11 , σ 22 , σ 33 , and σ eq along line B is shown for SAM-FG and SFE. An outstanding accuracy of SAM-FG is observed for all stress components in both cases.  Let us finally compare the computational cost between the resolution of SAM-FG and 3D SFE. In this section, 2D SFE calculations were performed owing to symmetries; in contrast to SFE, the proposed SAM-FG finite element resolution did not make use of these symmetries to reduce the computational cost. Let us pose the hypothetical case of the same catenoid subjected to whatever loading condition. In this case, 3D SFE would be required, but the same SAM-FG model implemented in COMSOL Multiphysics would be used after assigning the correct boundary conditions. Let us assume that tetrahedral elements are used in the SFE case. If one considers third-order Lagrangian elements in both   models and uses one-fourth of the thickness as the maximum element size of elements, the degrees of freedom would be 5.6 × 10 6 for SFE and 4.7 × 10 5 for SAM-FG finite element model. is partially shows one of the advantages of SAM-FG. e main advantage of SAM-FG lies in the fact that larger elements can be considered without losing significant accuracy in the evaluation of stresses. In the SFE model, larger elements would imply a significant accuracy loss because the variation of the elastic properties of the material would not be captured whereas in SAM-FG the constitutive equations hold the necessary information of the material properties across the thickness without meshing along the thickness direction. Much larger elements would add geometric inaccuracies in the SFE but not in the shell finite element model: this is also the advantage of homogeneous shell models over SFE.

Conclusions
A new model called SAM-FG was developed for elastostatic problems of functionally graded shells whose thickness may range from thin to moderately thick. is model is inspired by the SAM-H model conceived by Domínguez-Alvarado and Díaz-Díaz in [34] and limited to homogeneous shells. In contrast to the stress approximation in SAM-H, the degree of the polynomial expansion used in SAM-FG to approximate the stress field across the shell thickness is variable. is degree is given by a polynomial fit of order r of the components of the in-plane stiffness matrix. e polynomial degrees of the approximation of in-plane stresses, out-ofplane shear stresses, and out-of-plane normal stresses are r + 1, r + 3, and r + 4, respectively. e stress field verifies the boundary conditions at the faces of the shell and the 3D equilibrium equations. e use of Hellinger-Reissner functional allowed to identify 5 generalized displacements that are identical to those in SAM-H. e application of Reissner's variational theorem yielded the generalized equilibrium equations, boundary conditions, and constitutive equations. Except for the last ones mentioned, these equations turned out to be the same as in SAM-H. e generalized constitutive equations of SAM-FG have a similar form as those in SAM-H, but the generalized compliances appearing in these equations are radically different because of the material heterogeneity allowed in SAM-FG. e model equations were implemented and solved in the commercial finite element software COMSOL Multiphysics (5.3a version). To test the accuracy of SAM-FG, the model was first applied to the case of bending of a simply supported FGM plate and its results were first compared to those obtained by solid finite elements (SFE) and other models found in the literature that use more generalized displacements than SAM-FG. en, for the case of curved shells, SAM-FG results were compared to those obtained by solid finite elements (SFE) for the following problems of internally pressurized structures made of functionally graded materials: a hollow sphere, a hollow cylinder, and a catenoid. SAM-FG proved to yield accurate predictions of displacements and stresses for different thickness-to-radius ratios and through-the-thickness laws followed by the elastic properties of the functionally graded material. Owing to its accuracy and to the simplicity of the model (only 5 kinematic fields are considered), SAM-FG is now an interesting alternative to SFE calculations for the study of stresses and displacements in functionally graded shells subjected to mechanical loads.
A priority for future work will be the extension to dynamic problems and the inclusion of inelastic strains such as thermal expansion strains for enabling the analysis of the effects of a temperature variation on stresses. For this last case, if the values of the temperature field are not known, the determination of the temperature field by means of a solid finite element technique would imply a high computational cost. us, it would be worth developing a simplified model which, for example, approximates the heat flux vector by polynomials in a similar manner that SAM-FG approximates stresses. Finally, the development of an optimization algorithm based on artificial intelligence to predict the best material distribution for a functionally graded shell to obtain a particular response of the structure is envisioned as well.

Data Availability
e equations used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.