#### Abstract

This 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. The model is an extension of the SAM-H model, 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 tensor may depend on the thickness coordinate. The model starts with a tunable polynomial approximation of the 3D stress field that contains the same generalized forces as SAM-H. This stress approximation verifies the 3D equilibrium equations and the stress boundary conditions at the faces of the shell. As in SAM-H, 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. The first one deals with a hollow sphere made of an isotropic functionally graded material. The 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. The 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.

#### 1. Introduction

The 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–5] and exploit their benefits with new applications [6–9]. The 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.

The 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-of-plane stresses were neglected. The 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]. This 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 FSDT approach 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. This 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. The 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. The present paper is structured in 4 sections. The first section details the notations adopted. The next section recalls briefly the assumptions adopted in SAM-H and how the equations of the model are obtained. The 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.

#### 2. 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; is the thickness of the shell; are the orthogonal curvilinear coordinates considered and , , and are their related unit vectors; is the middle surface located at ; and are the coordinates along the lines of curvatures at ; the Lamé coefficients at the middle surface are and ; the principal curvatures are and .(iii)First- and second-order tensors of the 2D space defined by the basis are noted with an underlined and double-underlined boldface character, respectively (for example, and ); and are curvature tensors.(iv)At the outer () and inner () faces of the shell, the stress vectors and are applied, respectively. These vectors are defined as follows:(v) denotes the 3D divergence operator applied to the second-order tensor : where vectors and take the array form: and where and .(v) and denote the divergence operators applied to second-order 2D tensors and first-order 2D tensors, respectively; and are the 2D gradient operators applied to first-order tensor and scalar fields, respectively.

Additionally, the following assumptions are adopted:(1)In all calculations, the terms and are considered, but higher-order terms and () are neglected(2)The 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) is a 3D displacement vector. The superindex indicates that the field is not necessarily the solution of the problem.(ii) is a second-order stress tensor.(iii) and are the boundaries where displacement and stress vectors are imposed, respectively.(iv) is the Green–Cauchy 3D strain tensor.(v) is the 3D, 4th-order compliance tensor, and its components are .(vi) is the outward-pointing normal vector at the boundary of .(vii) and are the imposed displacement and stress vectors at the boundaries.

#### 3. Recalling SAM-H Equations

The development of SAM-FG follows almost the same procedure as SAM-H. The 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.

##### 3.1. Stress Approximation and Generalized Forces

The stress components are approximated by -polynomials. The following basis of 4th-degree, -polynomials is considered:

The following stress approximation is obtained:where 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. The 2D tensors of generalized forces are membrane forces , shear forces , and moments ; their components are defined as follows:

##### 3.2. 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 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.

##### 3.3. Generalized Constitutive Equations

SAM-H assumes that the material is homogeneous and orthotropic and that one of its orthotropic directions is direction 3. The following compliance tensors and constants are defined:

All the abovementioned compliances are assumed uniform through the thickness of the shell.

The following engineering vectors of generalized forces and strains are defined:where superscript 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:where(i) is an compliance matrix(ii) and are 8-component compliance vectors(iii), , and are compliance matrices

Besides, the aforementioned stationarity of the HR functional yields the edge conditions on the generalized displacements.

#### 4. SAM-FG Model

##### 4.1. 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 . As applied in the classical laminate theory (CLT), in a moderately thick FG shell, far enough from the edges, the 3D in-plane 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 and are nondimensional vector fields that depend on the in-plane coordinates. Let us define the plane-stress stiffness matrix ; 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 ) by a polynomial fit of order . In this manner,where are constant matrices and () is an orthogonal basis of polynomials ensuring the following condition:

The first 5 polynomials of this basis are shown in equation (6); in the other polynomials, “1” is chosen as the leading coefficient, for example,

Let us define the following vectors of generalized forces and moments:where superscript 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 with the 6 strains appearing in vectors and . After solving for the strain vectors, one obtainswhere , , , and are compliance matrices. By introducing the equation above in equation (25), one obtains the approximation of in-plane stresses of the SAM-FG model:

These in-plane stresses are -polynomials of degree :where 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).

The 3D equilibrium equation is used to deduce successively the degree of the polynomial approximation across the thickness for the shear stresses and normal stress . According to equations (2)–(4), the polynomial degrees of and turn out to be and , respectively. The coefficients of the polynomial are obtained by making use of a set of linear equations. The first two equations are deduced from the 3D boundary condition at the faces of the shell: and . The third equation comes from the definition of the generalized shear forces in equation (9). The other equations are obtained by following the next steps:(i)The component of the vectorial equation is selected.(ii)This component is written as a series expansion of polynomials . For , the coefficients of these polynomials are zero or negligible because of being multiplied by second-order terms: , , or .(iii)The searched equations come from solving for the coefficients of polynomials () in the aforementioned series.

In this manner, the proposed expression for the out-of-plane shear stresses is given bywhere are stress coefficients which can be written as linear combinations of , , and the derivatives of and . For the determination of the coefficients of the polynomial approximation of the normal stress , a similar method is applied, but this time, the boundary conditions and are used and the third component of the vectorial equation is selected. The proposed expression for the out-of-plane normal stress is given bywhere the stress coefficients 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).

##### 4.2. 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 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 , , , and are the surface densities of elastic energy defined by (summation over , and ):

It is worth defining the total surface density of elastic energy . Let us recall that the energies and 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:

The in-plane stress coefficients are expressed as linear combinations of the generalized membrane forces and moments. The coefficients and originally involve the divergences of the generalized forces , , and . 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. Thus, the sum of surface densities of energy appearing in equation (32) can be expressed as a linear combination of the following terms:where , , are the components of vector defined in equation (18) and represent the generalized membrane forces and moments, are the generalized shear forces, and represents the integrals defined in equation (34). If possible, the integrals defined in the terms are determined analytically; otherwise, they are calculated numerically. In this manner, takes the form:where , , , , , and 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 ; 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 termin equation (5) yields the following constitutive equations:where are the components of vector defined in equation (18). These constitutive equations can be written in a compact form as follows:

The components of matrices , , , and and those of vectors , , and can easily be identified with the compliance coefficients used in equation (38).

##### 4.3. Practical Use of SAM-FG

The principal orthotropic directions of the material are , , and ; the unit vectors associated with these directions are as follows: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. These properties depend on . Some authors, like Altenbach and Eremeyev in [15], analyze the particular case when these properties follow a power law of the thickness coordinate:where represents the elastic properties in the orthotropic directions; superscripts and denote the value of the property at the outer and inner surfaces, respectively; is a constant which modifies the power law. In Figure 1, the example of Young’s modulus is plotted against the normalized thickness coordinate () for = 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 , , , , , , , , , and the angle being defined, by making use of rotation operations, one obtains the components of the plane-stress stiffness matrix appearing in equation (21) and the components of compliances , , , and in the basis .

After this, the degree of the polynomial fit of the components of is selected and the constant matrices () appearing in equations (22) and (28) are determined. The components of are required for the calculation of the stress coefficients . Then, the integrals in equation (34) can be determined numerically or if possible, analytically. The 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 and are and stiffness matrices, respectively; and are 6-dimensional vectors; and and are 2-dimensional square matrices. The expressions of these matrices and vectors are given by

In the previous definition, and are and matrices, respectively. All the components of and are zero, except for

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 , , , , , and . 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. The LiveLink module in COMSOL for MATLAB would help to make this task automatically.

Let us finally emphasize that, whatever the degree of the polynomial fit of the components of the stiffness matrix , SAM-FG uses the same number of generalized displacements, strains, and forces. The polynomial fit is a preprocessing task similar to the determination of the 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. These 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.

#### 5. Results and Discussion

In this section, SAM-FG is applied to 4 mechanical problems. In all cases, a 5th-order polynomial fit () was enough to accurately approximate the stiffness coefficients. Higher values of 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.

##### 5.1. Plate Subjected to a Sinusoidal Pressure

The first case study consists of a simply supported square plate lying in the volume defined by , , and ( = 1 m and = 0.1 m). The following pressure is applied at the upper face of the plate:where = 1 MPa and . This 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]. These authors have used a displacement approach consisting of more than 5 generalized displacements (degrees of freedom). The 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 = 70 GPa and at the upper face = 380 GPa. Several values of the exponent 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. The selected maximum element size was 17.5 mm for both SFE and SAM-FG. The 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 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 and , the errors presented by SAM-FG are below 1%. For , , and , nonetheless, SAM-FG shows greater accuracy than all references in almost all the cases. This effect is accentuated in for , where the results presented in [21] have an error of 7.82% and the one obtained by SAM-FG is only 0.29%. This proves the benefits of the present model.

Let us plot some of the dimensionless stresses evaluated by SAM-FG across the thickness at the locations specified in equation (49). Figures 2 and 3 display the evolution of , , and for and , respectively. The results obtained with SFE are also displayed. In both cases ( and ), the stresses are accurately evaluated by SAM-FG.

##### 5.2. Hollow Sphere Subjected to an Internal Pressure

In this section, a hollow sphere subjected to an internal pressure is considered. The 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 modulus and Poisson’s ratio are and whilst those at the outer surface are and , respectively. The radius of the middle surface is = 1 m and the thickness-to-radius ratio is . The applied pressure at the internal face of the shell is . Owing to the spherical symmetry, the 3D equations of the problem are reduced to a 1D problem which was solved by making use of the Mathematics module of the commercial finite element software COMSOL Multiphysics 5.3a; the solution obtained is called the 3D solution. For the SAM-FG model, the spherical symmetry helps to reduce the 2D equations of the model into a set of algebraic equations. SAM-FG results are compared to those of the 3D solution to test the accuracy of the model. SAM-FG yields a generalized displacement which can be compared to the average value of the 3D solution across the thickness of the sphere. The 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 and . Figure 4 shows the graph of the relative error against the heterogeneity ratio . These 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. The 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 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:(1), , and (2), , and (3), , and (4), , and

In Figure 6, the relative error is plotted against the thickness-to-radius ratio . As increases, the absolute value of the relative error increases. For moderately thick shells (), is smaller than 0.8%, and for thicker shells (), this error remains small (). Therefore, 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 . The error is greater for materials 2 and 4 ( and ) 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 and the normalized radial stress are plotted against the normalized position through the thickness , respectively. It can be seen that SAM-FG predicts accurately both the in-plane stress and the radial stress. The maximum and minimum values of the radial stress are perfectly calculated by SAM-FG because in this model, the approximated stress verifies the boundary conditions at the faces of the shell.

##### 5.3. Hollow Cylinder Subjected to an Internal Pressure

In this section, a hollow cylinder subjected to a 1 MPa internal pressure is considered. The ends of the cylinder are axially blocked. The radius of the middle surface is = 1 m and the thickness-to-radius ratio is . The material is orthotropic and its principal directions 1, 2, and 3 coincide with the axial, azimuthal, and radial directions, respectively. The elastic moduli at the internal (superscript “”) and external (superscript “”) surfaces of the shell are shown in Table 2. Poisson’s ratios , , and are assumed to be uniform and equal to 0.3. These elastic properties vary through the thickness of the cylinder following the power law shown in equation (42). Several values of the exponent 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 (axial stress), (hoop stress), and (radial normal stress) are plotted against the normalized position for , , and , respectively. In all cases, the normal stress calculated by SAM-FG verifies the boundary conditions ( at and at ). For and , the normal stresses and 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 provided by SAM-FG with the average value of the 3D solution across the thickness of the cylinder. In Figure 13, the graphs of and are plotted against the exponent *n* appearing in the power law. The evolution of the curve of is accurately approximated by SAM-FG. In this figure, the relative error in the calculation of the average value of the radial displacement is also plotted. The absolute value of the relative error is less than 0.4%; this proves the great accuracy of SAM-FG.

##### 5.4. Catenoid Subjected to an Internal Pressure

The 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 plane (see Figure 14):

Since the principal curvatures and are not uniform, the thickness-to-radius ratio is not constant. The maximum ratio and average ratio for the catenoid are as follows:

The material is isotropic and its Young’s modulus and Poisson’s ratio vary linearly across the thickness in such a manner that these properties are GPa and at the internal face and GPa and at the external face. The 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 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. The 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 ( plane) and also along lines A and B shown in Figure 14. In Figure 16, the stress components and 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 and stresses: the relative errors in the maximum evaluation of and are 0.47% and 0.27%, respectively. In Figure 17, the normal stress and the equivalent stress (von Mises stress) are plotted across the section. Both stresses are evaluated with great accuracy by SAM-FG.

In Figure 18, the stress components , , , and are plotted against the normalized position along line A for both SAM-FG and SFE. The internal and external faces of the shell are located at and 0.5, respectively; as expected, the normal stress 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 , , , and 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 for SFE and for SAM-FG finite element model. This partially shows one of the advantages of SAM-FG. The 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.

#### 6. 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. This 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. This degree is given by a polynomial fit of order of the components of the in-plane stiffness matrix. The polynomial degrees of the approximation of in-plane stresses, out-of-plane shear stresses, and out-of-plane normal stresses are , , and , respectively. The stress field verifies the boundary conditions at the faces of the shell and the 3D equilibrium equations. The use of Hellinger–Reissner functional allowed to identify 5 generalized displacements that are identical to those in SAM-H. The 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. The 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. The 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. Then, 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. Thus, 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

The equations used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.