# A Stress Approach Model of Functionally Graded Shells

Authors: Carlos Humberto Rubio Rascón; Alberto Díaz Díaz; Axel Fernando Domínguez Alvarado
Journal: Mathematical Problems in Engineering
 (2020)

---

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


---

## Body

## 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; h is the thickness of the shell; ξ1,ξ2,ξ3 are the orthogonal curvilinear coordinates considered and e1, e2, and e3 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 a1 and a2; the principal curvatures are κ1 and κ2.(iii)
First- and second-order tensors of the 2D space defined by the basise1,e2 are noted with an underlined and double-underlined boldface character, respectively (for example, U¯ and M¯¯); κ¯¯=κ1e1⊗e1+κ2e2⊗e2 and κ′¯¯=κ2e1⊗e1+κ1e2⊗e2 are curvature tensors.(iv)
At the outer (ξ3=h/2) and inner (ξ3=−h/2) faces of the shell, the stress vectors s+ and s− are applied, respectively. These vectors are defined as follows:(1)s+ξ1,ξ2=τ¯++σ+e3,s−ξ1,ξ2=τ¯−+σ−e3,whereτ¯+=τ1+e1+τ2+e2,τ¯−=τ1−e1+τ2−e2.(v)
divσ denotes the 3D divergence operator applied to the second-order tensor σ:(2)divσ=1A1A2A3v1+v2,where vectorsv1 and v2 take the array form:(3)v1=∂A2A3σ11∂ξ1+∂A1A3σ12∂ξ2+∂A1A2σ13∂ξ3∂A2A3σ21∂ξ1+∂A1A3σ22∂ξ2+∂A1A2σ23∂ξ3∂A2A3σ31∂ξ1+∂A1A3σ32∂ξ2+∂A1A2σ33∂ξ3,and(4)v2=A2a2∂a1∂ξ2σ21−A1a1∂a2∂ξ1σ22+A2a1R1σ31A1a1∂a2∂ξ1σ12−A2a2∂a1∂ξ2σ11+A1a2R2σ32−A2a1R1σ11−A1a2R2σ22,whereA1=a11+ξ3κ1 and A2=a21+ξ3κ2.(v)
div¯ and div denote the divergence operators applied to second-order 2D tensors and first-order 2D tensors, respectively; grad¯¯ and grad¯ 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ξ3/Rα and h/Rα are considered, but higher-order terms ξ3/Rαn and h/Rαn (n≥2) 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 shellAssuming small displacements and strains, the Hellinger–Reissner (HR) functional [35] is expressed as follows:(5)H.R.u∗,σ∗=∫Ωσ∗:εu∗−12σ∗:S:σ∗dΩ−∫∂Ωuσ∗⋅n⋅u∗−ugdS−∫∂Ωssg⋅u∗dS,where(i)
u∗ 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)
∂Ωu and ∂Ω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 Sijkl.(vi)
n is the outward-pointing normal vector at the boundary of Ω.(vii)
ug and sg 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ξ3-polynomials. The following basis of 4th-degree, ξ3-polynomials is considered:(6)P0ξ3=1,P1ξ3=ξ3h,P2ξ3=−6ξ3h2+12,P3ξ3=−2ξ3h3+310ξ3h,P4ξ3=−143ξ3h4+ξ3h2−140.The following stress approximation is obtained:(7)σαβ∗ξ1,ξ2,ξ3=∑n=01σαβn∗ξ1,ξ2Pnξ3,σα3∗ξ1,ξ2,ξ3=∑n=03σα3n∗ξ1,ξ2Pnξ3,σ33∗ξ1,ξ2,ξ3=∑n=04σ33n∗ξ1,ξ2Pnξ3,where σijn∗ 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 N¯¯∗, shear forces Q¯∗, and moments M¯¯∗; their components are defined as follows:(8)Nαβ∗=∫−h/2h/2σαβ∗1+ξ3R3−βdξ3,(9)Qα∗=∫−h/2h/2σα3∗1+ξ3R3−αdξ3,(10)Mαβ∗=∫−h/2h/2ξ3σαβ∗1+ξ3R3−βdξ3.


### 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:(11)U1∗=∫−h/2h/2P0hu1∗dξ3,U2∗=∫−h/2h/2P0hu2∗dξ3,U3∗=∫−h/2h/2P0hu3∗dξ3,ϕ1∗=∫−h/2h/212P1h2u1∗dξ3,ϕ2∗=∫−h/2h/212P1h2u2∗dξ3.Also, the following generalized strains are obtained:(12)ε¯¯∗=grad¯¯U¯∗+U3∗κ¯¯,χ¯¯∗=grad¯¯Φ¯∗,d¯∗=grad¯U3∗+Φ¯∗−κ¯¯⋅U¯∗,where(13)U¯∗=U1∗e1+U2∗e2,Φ¯∗=ϕ1∗e1+ϕ2∗e2.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:(14)div¯N¯¯+κ¯¯⋅Q¯+h2κ1+κ2τ¯+−τ¯−+τ¯++τ¯−=0¯,(15)divQ¯+σ+−σ−2hκ1+κ2+σ++σ−−N11κ1−N22κ2=0,(16)div¯M¯¯−Q¯+h24κ1+κ2τ¯++τ¯−+h2τ¯+−τ¯−=0¯.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:(17)S¯¯¯¯=∑α,β,γ,δ=12Sαβγδeα⊗eβ⊗eγ⊗eδ,Sc¯¯=2∑α,β=12Sαβ33eα⊗eβ,SQ¯¯=4∑α,β=12Sα3β3eα⊗eβ,Sσ=S3333.All the abovementioned compliances are assumed uniform through the thickness of the shell.The following engineering vectors of generalized forcesF and strains ε are defined:(18)F=N11,N22,N12,N21,M11,M22,M12,M21t,ε=ε11,ε22,ε12,ε21,χ11,χ22,χ12,χ21t,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:(19)ε=CF+σ+c++σ−c−,d¯=DQ¯+D+τ¯++D−τ¯−,where(i)
C is an 8×8 compliance matrix(ii)
c+ and c− are 8-component compliance vectors(iii)
D, D+, and D− are 2×2 compliance matricesBesides, the aforementioned stationarity of the HR functional yields the edge conditions on the generalized displacements.


## 3.1. Stress Approximation and Generalized Forces

The stress components are approximated byξ3-polynomials. The following basis of 4th-degree, ξ3-polynomials is considered:(6)P0ξ3=1,P1ξ3=ξ3h,P2ξ3=−6ξ3h2+12,P3ξ3=−2ξ3h3+310ξ3h,P4ξ3=−143ξ3h4+ξ3h2−140.The following stress approximation is obtained:(7)σαβ∗ξ1,ξ2,ξ3=∑n=01σαβn∗ξ1,ξ2Pnξ3,σα3∗ξ1,ξ2,ξ3=∑n=03σα3n∗ξ1,ξ2Pnξ3,σ33∗ξ1,ξ2,ξ3=∑n=04σ33n∗ξ1,ξ2Pnξ3,where σijn∗ 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 N¯¯∗, shear forces Q¯∗, and moments M¯¯∗; their components are defined as follows:(8)Nαβ∗=∫−h/2h/2σαβ∗1+ξ3R3−βdξ3,(9)Qα∗=∫−h/2h/2σα3∗1+ξ3R3−αdξ3,(10)Mαβ∗=∫−h/2h/2ξ3σαβ∗1+ξ3R3−βdξ3.


## 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:(11)U1∗=∫−h/2h/2P0hu1∗dξ3,U2∗=∫−h/2h/2P0hu2∗dξ3,U3∗=∫−h/2h/2P0hu3∗dξ3,ϕ1∗=∫−h/2h/212P1h2u1∗dξ3,ϕ2∗=∫−h/2h/212P1h2u2∗dξ3.Also, the following generalized strains are obtained:(12)ε¯¯∗=grad¯¯U¯∗+U3∗κ¯¯,χ¯¯∗=grad¯¯Φ¯∗,d¯∗=grad¯U3∗+Φ¯∗−κ¯¯⋅U¯∗,where(13)U¯∗=U1∗e1+U2∗e2,Φ¯∗=ϕ1∗e1+ϕ2∗e2.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:(14)div¯N¯¯+κ¯¯⋅Q¯+h2κ1+κ2τ¯+−τ¯−+τ¯++τ¯−=0¯,(15)divQ¯+σ+−σ−2hκ1+κ2+σ++σ−−N11κ1−N22κ2=0,(16)div¯M¯¯−Q¯+h24κ1+κ2τ¯++τ¯−+h2τ¯+−τ¯−=0¯.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:(17)S¯¯¯¯=∑α,β,γ,δ=12Sαβγδeα⊗eβ⊗eγ⊗eδ,Sc¯¯=2∑α,β=12Sαβ33eα⊗eβ,SQ¯¯=4∑α,β=12Sα3β3eα⊗eβ,Sσ=S3333.All the abovementioned compliances are assumed uniform through the thickness of the shell.The following engineering vectors of generalized forcesF and strains ε are defined:(18)F=N11,N22,N12,N21,M11,M22,M12,M21t,ε=ε11,ε22,ε12,ε21,χ11,χ22,χ12,χ21t,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:(19)ε=CF+σ+c++σ−c−,d¯=DQ¯+D+τ¯++D−τ¯−,where(i)
C is an 8×8 compliance matrix(ii)
c+ and c− are 8-component compliance vectors(iii)
D, D+, and D− are 2×2 compliance matricesBesides, 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εαβ3D. 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:(20)vεξ1,ξ2,ξ3=ε113Dε223D2ε123D≃vε,0ξ1,ξ2P0ξ3+vε,1ξ1,ξ2P1ξ3,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 R3Dξ3; assuming a plane-stress state, the vector of in-plane stress components is given by(21)vσξ1,ξ2,ξ3=σ113Dσ223Dσ123D=R3D⋅vε≃R3D⋅vε,0P0ξ3+vε,1P1ξ3.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,(22)R3Dξ3≃∑i=0rPiξ3R3D,i,where R3D,i are constant matrices and Pi (1≤i≤n) is an orthogonal basis of polynomials ensuring the following condition:(23)∫−h/2h/2Pmξ3Pnξ3dξ3=0,ifm≠n.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,(24)P6ξ3=ξ3h6−1544ξ3h4+5176ξ3h2−514784.Equations (21) and (22) yield(25)vσ≃∑i=0rPiR3D,ivε,0+∑i=0rPiR3D,iP1vε,1.Let us define the following vectors of generalized forces and moments:(26)vNM=N11,N22,N12+N21,M11,M22,M12+M21t,vN=N11,N22,N12,N21t,vM=M11,M22,M12,M21t,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 vNM with the 6 strains appearing in vectors vε,0 and vε,1. After solving for the strain vectors, one obtains(27)vε,0=1hCN,0⋅vN+1h2CM,0⋅vM,vε,1=1hCN,1⋅vN+1h2CM,1⋅vM,where CN,0, CN,1, CM,0, and CM,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:(28)σ11∗σ22∗σ12∗=1h2∑i=0rPiR3D,i⋅hCN,0+P1CN,1⋅N11∗N22∗N12∗N21∗+CM,0+P1CM,1⋅M11∗M22∗M12∗M21∗.These in-plane stresses areξ3-polynomials of degree r+1:(29)σαβ∗ξ1,ξ2,ξ3=∑n=0r+1σαβn∗ξ1,ξ2Pnξ3,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).The 3D equilibrium equationdivσ∗=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. The r+4 coefficients of the σα3∗ polynomial are obtained by making use of a set of r+4 linear equations. The first two equations are deduced from the 3D boundary condition at the faces of the shell: σα3∗ξ3=h/2=τα+ and σα3∗ξ3=−h/2=−τα−. The third equation comes from the definition of the generalized shear forces in equation (9). The other r+1 equations are obtained by following the next steps:(i)
The componentα of the vectorial equation divσ∗=0 is selected.(ii)
This component is written as a series expansion of polynomialsPi. For i>r+2, the coefficients of these polynomials are zero or negligible because of being multiplied by second-order terms: h2κ12, h2κ22, or h2κ1κ2.(iii)
The searched equations come from solving for the coefficients of polynomialsPi (2≤i≤r+2) in the aforementioned series.In this manner, the proposed expression for the out-of-plane shear stresses is given by(30)σα3∗ξ1,ξ2,ξ3=∑n=0r+3σα3n∗ξ1,ξ2Pnξ3,where σα3n∗ 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. The proposed expression for the out-of-plane normal stress is given by(31)σ33∗ξ1,ξ2,ξ3=∑n=0r+4σ33n∗ξ1,ξ2Pnξ3,where the stress coefficients σ33n∗ 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 ∫Ωσ∗:ε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:(32)12∫Ωσ∗:S:σ∗dΩ=∫ωwes∗+wen∗+weQ∗+wec∗dω,where wes∗, wen∗, weQ∗, and wec are the surface densities of elastic energy defined by (summation over α,β,γ, and δ):(33)wes∗=12∫−h/2h/2σαβ∗Sβαδγσγδ∗1+ξ3κ11+ξ3κ2dξ3,wen∗=12∫−h/2h/2σ33∗Sσσ33∗1+ξ3κ11+ξ3κ2dξ3,weQ∗=12∫−h/2h/2σα3∗SαβQσβ3∗1+ξ3κ11+ξ3κ2dξ3,wec∗=12∫−h/2h/2σαβ∗Sβαcσ33∗1+ξ3κ11+ξ3κ2dξ3.It is worth defining the total surface density of elastic energywe∗=wes∗+wen∗+weQ∗+wec∗. Let us recall that the energies wen∗ and wec∗ 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:(34)wes∗≃∑1≤α,β,γ,δ≤2∑m=0r+1∑n=0r+1σαβmσγδnAαβγδmn,wen∗≃∑m=0r+4∑n=0r+4σ33mσ33nA3333mn,wec∗≃∑1≤α,β≤2∑m=0r+1∑n=0r+4σαβmσ33nAαβ33mn,weQ∗≃∑1≤α,β≤2∑m=0r+3∑n=0r+3σα3mσβ3nAα3β3mn,whereAαβγδmn=12∫−h/2h/2Sαβγδξ3Pmξ3Pnξ31+ξ3κ1+κ2dξ3,A3333mn=12∫−h/2h/2Sσξ3Pmξ3Pnξ31+ξ3κ1+κ2dξ3,Aαβ33mn=12∫−h/2h/2Sαβcξ3Pmξ3Pnξ31+ξ3κ1+κ2dξ3,Aα3β3mn=12∫−h/2h/2SαβQξ3Pmξ3Pnξ31+ξ3κ1+κ2dξ3.The in-plane stress coefficientsσαβm are expressed as linear combinations of the generalized membrane forces and moments. The coefficients σα3m and σ33n 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. Thus, the sum of surface densities of energy we∗=wes∗+wen∗+weQ∗+wec∗ appearing in equation (32) can be expressed as a linear combination of the following terms:(35)FpFqAmn,Fpσ+Amn,Fpσ−Amn,σ−2Amn,σ+2Amn,σ−σ+Amn,τα+τβ+Amn,τα−τβ−Amn,τα+τβ−Amn,QαQβAmn,Qατβ+Amn,Qατβ−Amn,where 1≤p,q≤8, 0≤m,n≤4, Fp 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 Amn represents the integrals defined in equation (34). If possible, the integrals defined in the terms Amn are determined analytically; otherwise, they are calculated numerically. In this manner, we∗ takes the form:(36)we∗=∑1≤p,q≤8Cpq1FpFq+∑1≤p≤8FpCp2+σ++Cp2−σ−+∑1≤α,β≤2Dα,β1QαQβ+∑1≤α,β≤2QαDαβ2+τβ++Dαβ2−τβ−+Z,where Cpq1, Cp2+, Cp2−, Dαβ1, Dαβ2+, 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 σijm; 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(37)∫Ωσ∗:εu∗−12σ∗:S:σ∗dΩ,in equation (5) yields the following constitutive equations:(38)εp=∑1≤q≤8Cpq1+Cqp1Fq+Cp2+σ++Cp2−σ−,for1≤p≤8,dα=∑1≤β≤2Dαβ1+Dβα1Qβ+∑1≤β≤2Dαβ2+τβ++Dαβ2−τβ−,for1≤α≤2,where εp are the components of vector ε defined in equation (18). These constitutive equations can be written in a compact form as follows:(39)ε=C1+C1tF+σ+C2++σ−C2−,d¯=D1+D1tQ¯+D2+τ¯++D2−τ¯−.The components of matricesC1, D1, D2+, and D2− and those of vectors C2+, C2−, and C3 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 areL, T, and N; the unit vectors associated with these directions are as follows:(40)eL=cosθξ3e1+sinθξ3e2,eT=−sinθξ3e1+cosθξ3e2,eN=e3,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 ξ3. Some authors, like Altenbach and Eremeyev in [15], analyze the particular case when these properties follow a power law of the thickness coordinate:(41)ψξ3=ψint+ψext−ψint×ξ3h+12n,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.Figure 1
Young’s modulusE vs. normalized position ξ3/h assuming power laws with different exponent coefficients n.Once the functions representing the elastic propertiesEL, ET, EN, νLT, νLN, νTN, GLT, GLN, GTN, and the angle θ being defined, by making use of rotation operations, one obtains the components of the 3×3 plane-stress stiffness matrix R3D appearing in equation (21) and the components of compliances S¯¯¯¯, Sc¯¯, SQ¯¯, and Sσ in the basis e1,e2,e3.After this, the degreer of the polynomial fit of the components of R3Dξ3 is selected and the constant matrices R3D,i (0≤i≤r) appearing in equations (22) and (28) are determined. The components of R3D,i are required for the calculation of the stress coefficients σijn. Then, the integrals Aijklmn 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:(42)F=Kε+σ+k++σ−k−,Q¯=Ld¯+L+τ¯++L−τ¯−,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. The expressions of these matrices and vectors are given by(43)K=H2H1C1+C1tH2−1H1,k+=−KC2+,k−=−KC2−,L=D1+D1t−1,L+=−LD2+,L−=−LD2−.In the previous definition,H1 and H2 are 6×8 and 8×6 matrices, respectively. All the components of H1 and H2 are zero, except for(44)H111=H221=H331=H341=H451=H561=H671=H681=1,H112=H222=H332=H432=H542=H652=H762=H862=1,H462=κ2−κ1,H832=h212κ2−κ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 ofK, 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. The LiveLink module in COMSOL for MATLAB would help to make this task automatically.Let us finally emphasize that, whatever the degreer of the polynomial fit of the components of the stiffness matrix R3Dξ3, 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 6×6ABBD 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.


## 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εαβ3D. 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:(20)vεξ1,ξ2,ξ3=ε113Dε223D2ε123D≃vε,0ξ1,ξ2P0ξ3+vε,1ξ1,ξ2P1ξ3,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 R3Dξ3; assuming a plane-stress state, the vector of in-plane stress components is given by(21)vσξ1,ξ2,ξ3=σ113Dσ223Dσ123D=R3D⋅vε≃R3D⋅vε,0P0ξ3+vε,1P1ξ3.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,(22)R3Dξ3≃∑i=0rPiξ3R3D,i,where R3D,i are constant matrices and Pi (1≤i≤n) is an orthogonal basis of polynomials ensuring the following condition:(23)∫−h/2h/2Pmξ3Pnξ3dξ3=0,ifm≠n.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,(24)P6ξ3=ξ3h6−1544ξ3h4+5176ξ3h2−514784.Equations (21) and (22) yield(25)vσ≃∑i=0rPiR3D,ivε,0+∑i=0rPiR3D,iP1vε,1.Let us define the following vectors of generalized forces and moments:(26)vNM=N11,N22,N12+N21,M11,M22,M12+M21t,vN=N11,N22,N12,N21t,vM=M11,M22,M12,M21t,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 vNM with the 6 strains appearing in vectors vε,0 and vε,1. After solving for the strain vectors, one obtains(27)vε,0=1hCN,0⋅vN+1h2CM,0⋅vM,vε,1=1hCN,1⋅vN+1h2CM,1⋅vM,where CN,0, CN,1, CM,0, and CM,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:(28)σ11∗σ22∗σ12∗=1h2∑i=0rPiR3D,i⋅hCN,0+P1CN,1⋅N11∗N22∗N12∗N21∗+CM,0+P1CM,1⋅M11∗M22∗M12∗M21∗.These in-plane stresses areξ3-polynomials of degree r+1:(29)σαβ∗ξ1,ξ2,ξ3=∑n=0r+1σαβn∗ξ1,ξ2Pnξ3,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).The 3D equilibrium equationdivσ∗=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. The r+4 coefficients of the σα3∗ polynomial are obtained by making use of a set of r+4 linear equations. The first two equations are deduced from the 3D boundary condition at the faces of the shell: σα3∗ξ3=h/2=τα+ and σα3∗ξ3=−h/2=−τα−. The third equation comes from the definition of the generalized shear forces in equation (9). The other r+1 equations are obtained by following the next steps:(i)
The componentα of the vectorial equation divσ∗=0 is selected.(ii)
This component is written as a series expansion of polynomialsPi. For i>r+2, the coefficients of these polynomials are zero or negligible because of being multiplied by second-order terms: h2κ12, h2κ22, or h2κ1κ2.(iii)
The searched equations come from solving for the coefficients of polynomialsPi (2≤i≤r+2) in the aforementioned series.In this manner, the proposed expression for the out-of-plane shear stresses is given by(30)σα3∗ξ1,ξ2,ξ3=∑n=0r+3σα3n∗ξ1,ξ2Pnξ3,where σα3n∗ 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. The proposed expression for the out-of-plane normal stress is given by(31)σ33∗ξ1,ξ2,ξ3=∑n=0r+4σ33n∗ξ1,ξ2Pnξ3,where the stress coefficients σ33n∗ 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 ∫Ωσ∗:ε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:(32)12∫Ωσ∗:S:σ∗dΩ=∫ωwes∗+wen∗+weQ∗+wec∗dω,where wes∗, wen∗, weQ∗, and wec are the surface densities of elastic energy defined by (summation over α,β,γ, and δ):(33)wes∗=12∫−h/2h/2σαβ∗Sβαδγσγδ∗1+ξ3κ11+ξ3κ2dξ3,wen∗=12∫−h/2h/2σ33∗Sσσ33∗1+ξ3κ11+ξ3κ2dξ3,weQ∗=12∫−h/2h/2σα3∗SαβQσβ3∗1+ξ3κ11+ξ3κ2dξ3,wec∗=12∫−h/2h/2σαβ∗Sβαcσ33∗1+ξ3κ11+ξ3κ2dξ3.It is worth defining the total surface density of elastic energywe∗=wes∗+wen∗+weQ∗+wec∗. Let us recall that the energies wen∗ and wec∗ 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:(34)wes∗≃∑1≤α,β,γ,δ≤2∑m=0r+1∑n=0r+1σαβmσγδnAαβγδmn,wen∗≃∑m=0r+4∑n=0r+4σ33mσ33nA3333mn,wec∗≃∑1≤α,β≤2∑m=0r+1∑n=0r+4σαβmσ33nAαβ33mn,weQ∗≃∑1≤α,β≤2∑m=0r+3∑n=0r+3σα3mσβ3nAα3β3mn,whereAαβγδmn=12∫−h/2h/2Sαβγδξ3Pmξ3Pnξ31+ξ3κ1+κ2dξ3,A3333mn=12∫−h/2h/2Sσξ3Pmξ3Pnξ31+ξ3κ1+κ2dξ3,Aαβ33mn=12∫−h/2h/2Sαβcξ3Pmξ3Pnξ31+ξ3κ1+κ2dξ3,Aα3β3mn=12∫−h/2h/2SαβQξ3Pmξ3Pnξ31+ξ3κ1+κ2dξ3.The in-plane stress coefficientsσαβm are expressed as linear combinations of the generalized membrane forces and moments. The coefficients σα3m and σ33n 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. Thus, the sum of surface densities of energy we∗=wes∗+wen∗+weQ∗+wec∗ appearing in equation (32) can be expressed as a linear combination of the following terms:(35)FpFqAmn,Fpσ+Amn,Fpσ−Amn,σ−2Amn,σ+2Amn,σ−σ+Amn,τα+τβ+Amn,τα−τβ−Amn,τα+τβ−Amn,QαQβAmn,Qατβ+Amn,Qατβ−Amn,where 1≤p,q≤8, 0≤m,n≤4, Fp 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 Amn represents the integrals defined in equation (34). If possible, the integrals defined in the terms Amn are determined analytically; otherwise, they are calculated numerically. In this manner, we∗ takes the form:(36)we∗=∑1≤p,q≤8Cpq1FpFq+∑1≤p≤8FpCp2+σ++Cp2−σ−+∑1≤α,β≤2Dα,β1QαQβ+∑1≤α,β≤2QαDαβ2+τβ++Dαβ2−τβ−+Z,where Cpq1, Cp2+, Cp2−, Dαβ1, Dαβ2+, 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 σijm; 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(37)∫Ωσ∗:εu∗−12σ∗:S:σ∗dΩ,in equation (5) yields the following constitutive equations:(38)εp=∑1≤q≤8Cpq1+Cqp1Fq+Cp2+σ++Cp2−σ−,for1≤p≤8,dα=∑1≤β≤2Dαβ1+Dβα1Qβ+∑1≤β≤2Dαβ2+τβ++Dαβ2−τβ−,for1≤α≤2,where εp are the components of vector ε defined in equation (18). These constitutive equations can be written in a compact form as follows:(39)ε=C1+C1tF+σ+C2++σ−C2−,d¯=D1+D1tQ¯+D2+τ¯++D2−τ¯−.The components of matricesC1, D1, D2+, and D2− and those of vectors C2+, C2−, and C3 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 areL, T, and N; the unit vectors associated with these directions are as follows:(40)eL=cosθξ3e1+sinθξ3e2,eT=−sinθξ3e1+cosθξ3e2,eN=e3,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 ξ3. Some authors, like Altenbach and Eremeyev in [15], analyze the particular case when these properties follow a power law of the thickness coordinate:(41)ψξ3=ψint+ψext−ψint×ξ3h+12n,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.Figure 1
Young’s modulusE vs. normalized position ξ3/h assuming power laws with different exponent coefficients n.Once the functions representing the elastic propertiesEL, ET, EN, νLT, νLN, νTN, GLT, GLN, GTN, and the angle θ being defined, by making use of rotation operations, one obtains the components of the 3×3 plane-stress stiffness matrix R3D appearing in equation (21) and the components of compliances S¯¯¯¯, Sc¯¯, SQ¯¯, and Sσ in the basis e1,e2,e3.After this, the degreer of the polynomial fit of the components of R3Dξ3 is selected and the constant matrices R3D,i (0≤i≤r) appearing in equations (22) and (28) are determined. The components of R3D,i are required for the calculation of the stress coefficients σijn. Then, the integrals Aijklmn 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:(42)F=Kε+σ+k++σ−k−,Q¯=Ld¯+L+τ¯++L−τ¯−,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. The expressions of these matrices and vectors are given by(43)K=H2H1C1+C1tH2−1H1,k+=−KC2+,k−=−KC2−,L=D1+D1t−1,L+=−LD2+,L−=−LD2−.In the previous definition,H1 and H2 are 6×8 and 8×6 matrices, respectively. All the components of H1 and H2 are zero, except for(44)H111=H221=H331=H341=H451=H561=H671=H681=1,H112=H222=H332=H432=H542=H652=H762=H862=1,H462=κ2−κ1,H832=h212κ2−κ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 ofK, 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. The LiveLink module in COMSOL for MATLAB would help to make this task automatically.Let us finally emphasize that, whatever the degreer of the polynomial fit of the components of the stiffness matrix R3Dξ3, 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 6×6ABBD 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 (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.


### 5.1. Plate Subjected to a Sinusoidal Pressure

The first case study consists of a simply supported square plate lying in the volume defined byx∈0,a, y∈0,a, and z∈−h/2,h/2 (a = 1 m and h = 0.1 m). The following pressure is applied at the upper face of the plate:(45)p=q0sinλxsinλy,where q0 = 1 MPa and λ=π/a. 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 Eint = 70 GPa and at the upper face Eext = 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:(46)σ¯x=haq0σxx,ata2,a2,h2,σ¯y=haq0σyy,ata2,a2,h3,σ¯xy=haq0σxy,at0,0,−h3,σ¯xz=haq0σxz,at0,a2,0,σ¯yz=haq0σyz,ata2,0,h6.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 Table1, 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. This 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%. This proves the benefits of the present model.Table 1
Summary of the results obtained by SFE, SAM-FG, and some FG plate models (case of a simply supported plate).
nSourceσ¯xσ¯yσ¯xyσ¯yzσ¯xz1Ref. [19]3.0871.48940.6110.26220.2462Ref. [25]3.085011.48980.611110.190710.23817Ref. [20]3.0871.48940.6110.26220.2462Ref. [21]3.087821.490340.6107040.2547210.238405SAM-FG3.071741.50990.6135950.2507970.2387323D3.09711.50620.611180.250870.23829Err. [21]0.30%1.05%0.08%1.54%0.05%Err. SAM0.82%0.25%0.40%0.03%0.19%2Ref. [19]3.60941.39540.54410.27630.2265Ref. [25]3.606641.395750.544340.18070.22568Ref. [20]3.60941.39540.54410.27630.2265Ref. [21]3.616351.396380.5434210.2639080.222026SAM-FG3.585471.418620.5473840.2493450.2262143D3.62481.4130.543620.249530.22565Err. [21]0.23%1.18%0.04%5.76%1.61%Err. SAM1.08%0.40%0.69%0.07%0.25%5Ref. [19]4.24881.10290.57550.24290.2017Ref. [25]4.242931.105390.573680.173070.21609Ref. [20]4.24881.10290.57550.24290.2017Ref. [21]4.259831.10220.5749580.2377090.199975SAM-FG4.211411.129830.5807140.2325550.2175653D4.2581.11810.576510.232420.21693Err. [21]0.04%1.42%0.27%2.28%7.82%Err. SAM1.09%1.05%0.73%0.06%0.29%10Ref. [20]5.0890.87750.58940.20410.2198Ref. [21]5.089320.877430.589390.2054340.216235SAM-FG5.053290.900730.595230.2241530.2170473D5.06660.894720.592680.223620.2166Err. [21]0.45%1.93%0.56%8.13%0.17%Err. SAM0.26%0.67%0.43%0.24%0.21%Let us plot some of the dimensionless stresses evaluated by SAM-FG across the thickness at thex,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. The results obtained with SFE are also displayed. In both cases (n=2 and n=10), the stresses are accurately evaluated by SAM-FG.Figure 2
Dimensionless stresses for a FG plate (n=2).Figure 3
Dimensionless stresses for a FG plate (n=10).


### 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 areEint and νint whilst those at the outer surface are Eext and νext, respectively. The radius of the middle surface is R = 1 m and the thickness-to-radius ratio is η. The applied pressure at the internal face of the shell is p. 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 U3SAM which can be compared to the average value U33D 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(47)δU=U3SAM−U33DU33D.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 Eint/Eext. 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 Eint/Eext ratio. Once again, 5 values of Poisson’s ratio were considered. In each case, the relative error δσ is less than 3.6%.Figure 4
Relative displacement errors vs.Eint/Eext ratio for η = 0.1 (pressurized hollow sphere).Figure 5
Relative errors in the approximation of the maximum von Mises stress vs.Eint/Eext ratio for η = 0.1 (pressurized hollow sphere).Let us now analyze the effect of the thickness-to-radius ratioη on the accuracy of the model. Four material types are considered:(1)
Eint=2Eext, νint=0.4, and νext=0.2(2)
Eint=2Eext, νint=0.2, and νext=0.4(3)
Eint=0.5Eext, νint=0.4, and νext=0.2(4)
Eint=0.5Eext, νint=0.2, and νext=0.4In Figure6, 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%). 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 (ν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 σ11h/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. The 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.Figure 6
Relative displacement errors vs. thickness-to-radius ratioη (pressurized hollow sphere).Figure 7
Relative errors in the approximation of the maximum von Mises stress vs. thickness-to-radius ratioη (pressurized hollow sphere).Figure 8
Normalized in-plane stress vs.ξ3/h (case of material 1, pressurized hollow sphere).Figure 9
Normalized radial stress vs.ξ3/h (case of material 1, pressurized hollow sphere).


### 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 isR0 = 1 m and the thickness-to-radius ratio is η=0.1. 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 “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. These 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.Table 2
Elastic moduli (in GPa) at the internal and external surfaces of the cylinder.
SurfaceE1E2E3G12G13G23Internal50508201010External10010012301515Owing 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 Figures10–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=−1MPa 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.Figure 10
Stress componentsσ11, σ22, and σ33 vs. normalized position ξ3/h through the thickness for n=1/3 (pressurized hollow cylinder).Figure 11
Stress componentsσ11, σ22, and σ33 vs. normalized position ξ3/h through the thickness for n=1 (pressurized hollow cylinder).Figure 12
Stress componentsσ11, σ22, and σ33 vs. normalized position ξ3/h through the thickness for n=3 (pressurized hollow cylinder).Let us now compare the generalized displacementU3SAM provided by SAM-FG with the average value U33D of the 3D solution across the thickness of the cylinder. In Figure 13, the graphs of U3SAM and U33D are plotted against the exponent n appearing in the power law. The evolution of the curve of U33D 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. The absolute value of the relative error is less than 0.4%; this proves the great accuracy of SAM-FG.Figure 13
Average valueU3 of the radial displacement (3D and SAM-FG) and relative error δU vs. exponent n appearing in the power law (pressurized hollow cylinder).


### 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 ther−z plane (see Figure 14):(48)r=a∗coshz−aa,forz∈0,1m,a=1m.Figure 14
Boundary conditions for axisymmetric catenoid.Since the principal curvaturesκ1 and κ2 are not uniform, the thickness-to-radius ratio η is not constant. The maximum ratio ηmax and average ratio η¯ for the catenoid are as follows:(49)ηmax=0.1,η¯=∫ωηdω∫ωdω=0.071.The material is isotropic and its Young’s modulusE 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. 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 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 Figure15, 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.Figure 15
SFE and SAM-FG meshing.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.Figure 16
Stress componentsσ11 and σ22 (in MPa) obtained by SFE and SAM-FG; maximum and minimum values appear in red and black, respectively (pressurized catenoid).Figure 17
Stress componentsσ33 and σeq (in MPa) obtained by SFE and SAM-FG; maximum and minimum values appear in red and black, respectively (pressurized catenoid).In Figure18, 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. The 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.Figure 18
σ11, σ22, σ33, and σeq vs. normalized position ξ3/h along line A (pressurized catenoid).Figure 19
σ11, σ22, σ33, and σeq vs. normalized position ξ3/h along line B (pressurized catenoid).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 be5.6×106 for SFE and 4.7×105 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.


## 5.1. Plate Subjected to a Sinusoidal Pressure

The first case study consists of a simply supported square plate lying in the volume defined byx∈0,a, y∈0,a, and z∈−h/2,h/2 (a = 1 m and h = 0.1 m). The following pressure is applied at the upper face of the plate:(45)p=q0sinλxsinλy,where q0 = 1 MPa and λ=π/a. 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 Eint = 70 GPa and at the upper face Eext = 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:(46)σ¯x=haq0σxx,ata2,a2,h2,σ¯y=haq0σyy,ata2,a2,h3,σ¯xy=haq0σxy,at0,0,−h3,σ¯xz=haq0σxz,at0,a2,0,σ¯yz=haq0σyz,ata2,0,h6.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 Table1, 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. This 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%. This proves the benefits of the present model.Table 1
Summary of the results obtained by SFE, SAM-FG, and some FG plate models (case of a simply supported plate).
nSourceσ¯xσ¯yσ¯xyσ¯yzσ¯xz1Ref. [19]3.0871.48940.6110.26220.2462Ref. [25]3.085011.48980.611110.190710.23817Ref. [20]3.0871.48940.6110.26220.2462Ref. [21]3.087821.490340.6107040.2547210.238405SAM-FG3.071741.50990.6135950.2507970.2387323D3.09711.50620.611180.250870.23829Err. [21]0.30%1.05%0.08%1.54%0.05%Err. SAM0.82%0.25%0.40%0.03%0.19%2Ref. [19]3.60941.39540.54410.27630.2265Ref. [25]3.606641.395750.544340.18070.22568Ref. [20]3.60941.39540.54410.27630.2265Ref. [21]3.616351.396380.5434210.2639080.222026SAM-FG3.585471.418620.5473840.2493450.2262143D3.62481.4130.543620.249530.22565Err. [21]0.23%1.18%0.04%5.76%1.61%Err. SAM1.08%0.40%0.69%0.07%0.25%5Ref. [19]4.24881.10290.57550.24290.2017Ref. [25]4.242931.105390.573680.173070.21609Ref. [20]4.24881.10290.57550.24290.2017Ref. [21]4.259831.10220.5749580.2377090.199975SAM-FG4.211411.129830.5807140.2325550.2175653D4.2581.11810.576510.232420.21693Err. [21]0.04%1.42%0.27%2.28%7.82%Err. SAM1.09%1.05%0.73%0.06%0.29%10Ref. [20]5.0890.87750.58940.20410.2198Ref. [21]5.089320.877430.589390.2054340.216235SAM-FG5.053290.900730.595230.2241530.2170473D5.06660.894720.592680.223620.2166Err. [21]0.45%1.93%0.56%8.13%0.17%Err. SAM0.26%0.67%0.43%0.24%0.21%Let us plot some of the dimensionless stresses evaluated by SAM-FG across the thickness at thex,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. The results obtained with SFE are also displayed. In both cases (n=2 and n=10), the stresses are accurately evaluated by SAM-FG.Figure 2
Dimensionless stresses for a FG plate (n=2).Figure 3
Dimensionless stresses for a FG plate (n=10).


## 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 areEint and νint whilst those at the outer surface are Eext and νext, respectively. The radius of the middle surface is R = 1 m and the thickness-to-radius ratio is η. The applied pressure at the internal face of the shell is p. 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 U3SAM which can be compared to the average value U33D 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(47)δU=U3SAM−U33DU33D.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 Eint/Eext. 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 Eint/Eext ratio. Once again, 5 values of Poisson’s ratio were considered. In each case, the relative error δσ is less than 3.6%.Figure 4
Relative displacement errors vs.Eint/Eext ratio for η = 0.1 (pressurized hollow sphere).Figure 5
Relative errors in the approximation of the maximum von Mises stress vs.Eint/Eext ratio for η = 0.1 (pressurized hollow sphere).Let us now analyze the effect of the thickness-to-radius ratioη on the accuracy of the model. Four material types are considered:(1)
Eint=2Eext, νint=0.4, and νext=0.2(2)
Eint=2Eext, νint=0.2, and νext=0.4(3)
Eint=0.5Eext, νint=0.4, and νext=0.2(4)
Eint=0.5Eext, νint=0.2, and νext=0.4In Figure6, 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%). 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 (ν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 σ11h/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. The 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.Figure 6
Relative displacement errors vs. thickness-to-radius ratioη (pressurized hollow sphere).Figure 7
Relative errors in the approximation of the maximum von Mises stress vs. thickness-to-radius ratioη (pressurized hollow sphere).Figure 8
Normalized in-plane stress vs.ξ3/h (case of material 1, pressurized hollow sphere).Figure 9
Normalized radial stress vs.ξ3/h (case of material 1, pressurized hollow sphere).


## 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 isR0 = 1 m and the thickness-to-radius ratio is η=0.1. 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 “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. These 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.Table 2
Elastic moduli (in GPa) at the internal and external surfaces of the cylinder.
SurfaceE1E2E3G12G13G23Internal50508201010External10010012301515Owing 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 Figures10–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=−1MPa 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.Figure 10
Stress componentsσ11, σ22, and σ33 vs. normalized position ξ3/h through the thickness for n=1/3 (pressurized hollow cylinder).Figure 11
Stress componentsσ11, σ22, and σ33 vs. normalized position ξ3/h through the thickness for n=1 (pressurized hollow cylinder).Figure 12
Stress componentsσ11, σ22, and σ33 vs. normalized position ξ3/h through the thickness for n=3 (pressurized hollow cylinder).Let us now compare the generalized displacementU3SAM provided by SAM-FG with the average value U33D of the 3D solution across the thickness of the cylinder. In Figure 13, the graphs of U3SAM and U33D are plotted against the exponent n appearing in the power law. The evolution of the curve of U33D 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. The absolute value of the relative error is less than 0.4%; this proves the great accuracy of SAM-FG.Figure 13
Average valueU3 of the radial displacement (3D and SAM-FG) and relative error δU vs. exponent n appearing in the power law (pressurized hollow cylinder).


## 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 ther−z plane (see Figure 14):(48)r=a∗coshz−aa,forz∈0,1m,a=1m.Figure 14
Boundary conditions for axisymmetric catenoid.Since the principal curvaturesκ1 and κ2 are not uniform, the thickness-to-radius ratio η is not constant. The maximum ratio ηmax and average ratio η¯ for the catenoid are as follows:(49)ηmax=0.1,η¯=∫ωηdω∫ωdω=0.071.The material is isotropic and its Young’s modulusE 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. 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 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 Figure15, 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.Figure 15
SFE and SAM-FG meshing.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.Figure 16
Stress componentsσ11 and σ22 (in MPa) obtained by SFE and SAM-FG; maximum and minimum values appear in red and black, respectively (pressurized catenoid).Figure 17
Stress componentsσ33 and σeq (in MPa) obtained by SFE and SAM-FG; maximum and minimum values appear in red and black, respectively (pressurized catenoid).In Figure18, 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. The 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.Figure 18
σ11, σ22, σ33, and σeq vs. normalized position ξ3/h along line A (pressurized catenoid).Figure 19
σ11, σ22, σ33, and σeq vs. normalized position ξ3/h along line B (pressurized catenoid).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 be5.6×106 for SFE and 4.7×105 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 r 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 r+1, r+3, and r+4, 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.



---

*Source: 7385158-2020-08-13.xml*