A Layer-Wise Coupled Thermo-Elastic Shell Model for Three-Dimensional Stress Analysis of Functionally Graded Material Structures

: In this work, a coupled 3D thermo-elastic shell model is presented. The primary variables are the scalar sovra-temperature and the displacement vector. This model allows for the thermal stress analysis of one-layered and sandwich plates and shells embedding Functionally Graded Material (FGM) layers. The 3D equilibrium equations and the 3D Fourier heat conduction equation for spherical shells are put together into a set of four coupled equations. They automatically degenerate in those for simpler geometries thanks to proper considerations about the radii of curvature and the use of orthogonal mixed curvilinear coordinates α , β , and z . The obtained partial differential governing the equations along the thickness direction are solved using the exponential matrix method. The closed form solution is possible assuming simply supported boundary conditions and proper harmonic forms for all the unknowns. The sovra-temperature amplitudes are directly imposed at the outer surfaces for each geometry in steady-state conditions. The effects of the thermal environment are related to the sovra-temperature proﬁles through the thickness. The static responses are evaluated in terms of displacements and stresses. After a proper and global preliminary validation, new cases are presented for different thickness ratios, geometries, and temperature values at the external surfaces. The considered FGM is metallic at the bottom and ceramic at the top. This FGM layer can be embedded in a sandwich conﬁguration or in a one-layered conﬁguration. This new fully coupled thermo-elastic model provides results that are coincident with the results proposed by the uncoupled thermo-elastic model that separately solves the 3D Fourier heat conduction equation. The differences are always less than 0.5% for each investigated displacement, temperature, and stress component. The differences between the present 3D full coupled model and the the advantages of this new model are clearly shown. Both the thickness layer and material layer effects are directly included in all the conducted coupled thermal stress analyses.


Introduction
The thermal stress analysis, for innovative and modern aerospace industries, is fundamental for new projects about space vehicles and aircraft where high temperature gradients and fast changes in temperature are classical operational conditions [1][2][3][4].The study of the temperature gradient influence on strains and stresses is mandatory for the correct analysis of the performances of avant-garde structures.The main topic is related to the heat flux loads acting on the external structures of the airplanes and launch vehicles.These thermal loads and the mechanical loads can cause serious damage to the structure.For this reason, a proper mathematical formulation to analyse all these load effects must be developed.In addition, thanks to the technological improvement, the Functionally Graded Materials (FGMs), capable of continuously varying their mechanical and thermal properties along a given direction thanks to two or more constituent phases that vary over a defined volume, permit enhancing performances of the classical composite materials.The FGMs can vary their properties in each direction, but, in the aerospace sector, the thickness direction is the most common one.FGMs are even used for the thermal applications because they can better stand high temperatures at the same strength-to-weight ratio and there are not any layer interfaces [5].In this paper, a full coupled thermo-elastic model is presented where the 3D equilibrium equations and the 3D Fourier heat conduction equation for spherical shells are employed.In this way, the sovra-temperature and the three displacement components are unknown variables of the problem.The most important characteristics of full coupled thermo-elastic models were shown in [6][7][8][9][10][11][12] where divergence and gradient equations, constitutive equations, boundary conditions, variational principles, field equations, proportionality equations between the heat flux and gradient of the thermal variable, energy balance equations, and the initial conditions were deeply discussed.
Several papers about the thermal stress analysis of FGM structures, both numerical and analytical models, are present in the open literature.In order to better remark the new characteristics of this new formulation, thermo-mechanical models for plates and shells embedding FGM layers are grouped as those related to 1D exact solutions, 1D numerical models, 2D exact solutions, 2D numerical models, 3D exact solutions, and 3D numerical models.
For what concerns the literature about the 1D exact solutions, Kapuria et al.
[13] presented a third-order zigzag theory, together with the modified rule of mixtures for the effective modulus of elasticity for layered FGM beams.This model was evaluated for static and free vibration analyses.Ghiasian et al. [14] presented the buckling of FGM beams under different thermal loads.The beam was resting over a three-parameter elastic foundation with hardening/softening cubic nonlinearity, which acted in tension as well as in compression.Kiani and Eslami [15] showed the static and dynamic buckling of an FGM beam under uniform temperature rise and uniform compression.The material properties in [14,15] varied along the thickness direction.The Timoshenko beam model resting on a two-parameter non-linear elastic foundation was implemented in [16] in the case of the thermal buckling of FGM beams subjected to a temperature rise.Ma and Lee [17] presented a closed-form solution for the nonlinear static responses of FGM beams under uniform in-plane thermal loads.The three governing equations for the axial and transverse deformations of beams were reduced to a single nonlinear fourth-order integral-differential equation.The geometrically non-linear post-buckling load-deflection behaviour of FGM Timoshenko beams under in-plane thermal loadings was discussed in [18].The thermal loads were applied by providing a non-uniform temperature rise across the beam thickness in steady-state conditions.Zhang et al. [19] proposed a thermal buckling of ceramic-metal FGM beams subjected to a transversely non-uniform temperature rise.The investigation method was the symplectic theory in the Hamiltonian system.
In the framework of 1D numerical models, Chakraborty et al. [20] presented a beam element for thermoelastic analysis based on the first-order shear deformation theory where elastic and thermal properties changed along the thickness direction.The interpolating polynomials were created thanks to the exact solution of the same beam model.In [21], a thermo-elastic vibration analysis of FGM beams with general boundary conditions was presented.A higher-order shear beam deformation theory with material properties dependent on the temperature was used.An improved Finite Element Method (FEM) was shown in [22].Thanks to this new FE model, the transverse and axial vibrations of FGM beams under thermal fields and exposed to a moving mass were analysed.Esfahani et al. [23] proposed a thermal buckling and post-buckling analysis of FGM structures using the Timoshenko beam model resting on a non-linear elastic foundation.Both thermal and mechanical properties were functions of the temperature and position.Tang and Li [24] proposed FGM slender beam analyses.The numerical model was created via the principle of the minimum for the total potential energy in order to derive the non-linear governing equations of the structures.The buckling of FGM beams with different boundary conditions was presented in [25].The numerical model was based on the Timoshenko beam theory and the critical buckling load of the structures was evaluated using the Ritz method.Ziane et al. [26] presented a numerical model, using the Galerkin method, to analyse the thermal buckling of simply supported and clamped-clamped FGM box beams.
In the case of 2D exact solutions, Jahaveri and Eslami [27] presented stability and equilibrium equations for rectangular plates constituted of FGM layers under thermal loads.The equations were derived from the high-order shear deformation theory.Akbaş [28] proposed free vibration and static analyses of simply supported FGM plates; the porosity effect was included.Saad and Hadji [29] showed a thermal buckling problem for porous thick rectangular plates constituted of FGM layers using a high-order shear deformation theory.Sangeetha et al. [30] presented a closed form solution for FGM plates under thermal loads using a refined model based on the first-order shear deformation theory.The effects of thermal stresses were studied for several temperature variations across the thickness direction of the plate.In [31], Zenkour and Mashat presented thermal buckling responses using a Sinusoidal shear deformation Plate Theory (SPT).The governing equations were derived using SPT and solved in a closed form.Yaghoobi and Ghannad [32] proposed a thermal analysis for FGM cylinders under non-uniform heat fluxes.The governing equations were based on the first-order temperature theory and the energy method.Zeighami and Jafari [33] proposed a solution for the thermo-mechanical analysis of functionally graded carbon nanotube-reinforced composite plates with a central hole.This solution was possible thanks to the Lekhnitskii complex potential approach and the proper conformal mapping functions.
In the area of 2D numerical models, Praveen and Reddy [34] proposed a finite element able to take into account transverse shear strains, rotary inertia, and large rotations for FGM ceramic-metal plates.Static and dynamic analyses were conducted.Thai et al. [35] presented a four-unknown shear and normal deformation theory for static, dynamic and buckling analyses of FGM plates where the 3D material matrix was used.The system of equations was derived using the Galerkin weak form and the isogeometric analysis.In [36], a non-linear finite element model was presented to study the dynamic response of FGM structures under the thermal and mechanical harmonic loads.In this case, the FGM properties depended on the temperature and they varied in the thickness direction via a power law distribution.The effects of the material variation through the thickness and the size of the FGM were studied using the finite element method in [37] using the Crank-Nicolson-Galerkin scheme.Alibeigloo [38] showed the bending analysis of the FGM sandwich circular plates under thermo-mechanical loads.The used method was the Generalized Differential Quadrature method (GDQ) and the temperature distribution in the 3D form was computed by solving the heat conduction governing the equation in closed form.Hong [39] proposed a GDQ third-order shear deformation plate theory for FGM structures under thermal vibrations.Karakoti et al. [40] showed an eight-nodes isoparametric finite element in order to obtain a nonlinear transient response of porous FGM sandwich plates and shells.The FGM structure was subjected to blast loadings and thermal loads.Jooybar et al. [41] presented a numerical model where the equations of motion (and related boundary conditions), derived thanks to the Hamilton principle, were solved with the use of the differential quadrature method.The embedded FGM was temperature dependent.In [42], a thermal buckling analysis of the FGM sandwich plates, using an improved mesh-free Radial Point Interpolation Method (RPIM), was presented.This buckling formulation for plates was derived from an improvement of RPIM employing a new radial basis function.Therefore, the shape functions were built without any supporting fixing parameters based on the higher-order shear deformation plate theory.Qi et al. [43] showed a dynamic analysis for stiffened doubly curved sandwich composite panels with an FGM core and two isotropic layers under thermal loads.This analysis was based on von Kármán non-linear strain-displacement relationships and classical plate theory.The mathematical problem was solved by adopting the finite difference model and the Newmark method.Taj et al. [44] presented the static analysis of FGM plates using the HSDT (High Shear Deformation Theory) where the transverse shear stress was represented as quadratic along the thickness direction; the material properties of the FGM varied along the same direction.
For what concerns the 3D analytical solutions, they can be used to validate all the previously discussed models because they can consider all the peculiarities for geometry and lamination.Reddy and Cheng [45] proposed a theory for the bending of rectangular plates embedding FGM layers and piezoelectric actuators.In this way, when the FGM surface was subjected to thermal loads, displacements and stresses can be controlled.In [46], an analytical solution was presented for 3D steady and transient heat conduction problems of double-layer plates (including a coating layer and an FGM layer) with a local heat source.For this solution, the Poisson method and the layerwise approach were employed.Chen et al. [47] proposed a method based on state-space formulations with laminate approximations.The employed FGM was temperature dependent.Ootao e Tanigawa [48] showed a theoretical method for transient thermoelastic problems involving orthotropic FGM rectangular plates and non-uniform heat supply.The transient 3D temperature was analysed with an exponential law in the thickness direction.The same authors also [49] proposed a theoretical analysis of a 3D thermal stress problem for FGM structures subjected to a partial heat supply in a transient state.Jabbari et al. [50] discussed an exact solution for the steady state thermo-elastic problem of 3D simply supported circular FGM plates.Thermal and mechanical loads were axisymmetrically applied at the outer surfaces.Vel and Batra [51] presented an exact solution for 3D deformations of simply supported FGM rectangular plates subjected to mechanical and thermal loads at the external surfaces.Proper temperature and displacement functions were used in order to satisfy the boundary conditions at the edges and to reduce the system of partial differential equations for the thermo-elastic problem.Liu [52] discussed a 3D axisymmetric FGM circular plate under thermal loads at outer surfaces.A proper temperature function for the thermal boundary conditions at the edges was used and the variable separation method was employed to reduce the order of the set of governing equations in steady-state heat conduction.Alibeigloo [53] showed a 3D thermo-elastic model for FGM rectangular plates with simply supported edges under thermo-mechanical loads.The analytical solutions for temperature, stress, and displacement fields were proposed by using the Fourier series and the statespace method.
In the framework of 3D numerical (FEM, meshless methods, and GDQ) models, a study about the thermal elastic residual stresses occurring in Ni-Al2O3, Ni-TiO2, and Ti-SiC FG plates, due to different temperature fields through the plate thickness, was presented in [54].A 3D eight-nodes isoparametric-layered finite element with three degrees of freedom per node was implemented here.In [55], Hajlaoui et al. proposed a modified first-order enhanced solid-shell element formulation for the thermal buckling of functionally graded shells.The material properties varied in the thickness direction via a power law.In [56], a 3D free vibration analysis of shells constituted of laminated FGMs was shown thanks to the use of the quadrature element method.The shell geometry was analysed both in thermal and non-thermal configurations.Burlayenko et al. [57] presented a 3D analysis for free vibrations of thermally FGM sandwich plates.The material properties varied along the thickness direction and the analysis was conducted using the ABAQUS FE code.A 3D analysis of FGM cylinders containing semi-elliptical circumferential surface crack and thermo-mechanical loading was presented in [58].The variation law for the Young modulus was exponential through the thickness.Naghdabadi and Kordkheili [59] proposed an FE formulation for the thermoelastic analysis of FGM plates and shells.The power law distribution for the composition of the constituent phase varied in the thickness direction.Qian and Batra [60] proposed a transient heat conduction analysis for thick FGM plates by using a higher-order plate theory and a meshless local Petrov-Galerkin method.Mian and Spencer [61] proposed models for both laminated and FGM structures in a Cartesian coordinate system.The comparisons between the two materials were discussed.
This new 3D coupled thermo-elastic shell model is valid for several geometries: spherical shells and cylindrical panels, cylinders, and plates.A closed-form solution is proposed using the simply supported hypotheses for the edges and the harmonic forms for the primary variables.The materials embedded in the proposed structures could be classical or functionally graded along the thickness direction.The 3D equilibrium equations coupled with the 3D heat conduction equation for shells are solved thanks to the exponential matrix method along the thickness direction.A layer-wise approach is employed.The present 3D exact coupled thermo-elastic shell model can be seen as the general case of the pure mechanical model already developed by Brischetto in [62,63] in the case of free vibration and bending analyses, respectively.The addition of the thermal-related equation in the orthogonal mixed curvilinear coordinates (see [64][65][66][67]) to 3D equilibrium relations creates a homogeneous differential equation set that can be solved via the procedure shown in [68,69].The new results are presented in terms of displacements, stresses, and temperature profiles.They can be used as benchmarks for the development and testing of new 3D, 2D, and 1D numerical models for thermal stress analyses of FGM structures.
The paper is organized as follows.Section 2 is about the coupled thermo-elastic governing equations for spherical shells and the related solution methodology.Section 3 covers the results and is split in a first subsection for preliminary assessments and a second subsection for new benchmarks.Section 4 provides the main conclusions.

3D Exact and Coupled Thermo-Elastic Governing Equations for Spherical Shells
This section shows the development of the 3D coupled exact thermo-elastic shell model for FGM structures.Plates, cylinders, cylindrical, and spherical shells (see the differences between the four geometries in Figure 1) are analysed using the same model thanks to appropriate considerations about the radii of the curvature R α and R β .
Geometries for assessments and benchmarks.Spherical shell and related particular cases as plate, cylindrical shell, and cylinder.
In Figure 1, the orthogonal mixed curvilinear reference system (α, β, z) has its origin in a corner.The in-plane directions α and β are parallel to the lateral curved sides and they lie on the middle surface Ω 0 .The Ω 0 surface is the reference surface for the computation of all the geometrical parameters.The thickness direction z is normal to the Ω 0 surface and directed from the bottom to the top surface.This model has, as unknown variables, the three displacement components u, v, and w in the three directions of the reference system and the scalar sovra-temperature θ.The 3D elastic equilibrium equations and the 3D Fourier heat conduction equation are put together in the same system that is solved as a system of second-order partial differential equations where the unknowns are displacements, temperature, and related derivatives created with respect to z.

3D Equilibrium and Heat Conduction Equations for Spherical Shells
The starting point for the mathematical formulation of the present 3D full coupled exact thermo-elastic shell model is the writing of the 3D equilibrium equations and the 3D Fourier heat conduction relation for spherical shells: Equations ( 1)-(3) are the 3D equilibrium equations, and they are linked with the 3D Fourier heat conduction equation in steady-state condition (Equation ( 4)) in order to couple the mechanical field with the thermal one.The 3D Fourier heat conduction relation in Equation ( 4) is proposed in a steady-state form.Therefore, the dependence on the time is discarded.R α and R β are the radii of curvature in the α and β directions, respectively, and they are constant values.The H α (z) and H β (z) coefficients presented in Equations ( 1)-( 4) are linear functions of the thickness coordinate z or z.They introduced the curvature effects of the shells and they are defined, for each direction, as: The term h shown in Equations ( 5)-( 7) represents the total thickness of the structure.h is always considered as constant.The variable z is in the range between −h/2 and h/2 and z is in the range between 0 and h.In Equation ( 4), the conduction coefficients κ * where κ k 1 (z), κ k 2 (z) and κ k 3 (z) are the conduction coefficients in the three directions of the mixed orthogonal reference system.If the k layer is an FGM, they depend on z. κ * 1 k (z), κ * 2 k (z), and κ * 3 k (z) also depend on the curvature via the use of H α (z) and H β (z).It must be noted that, in the case of a plate, the κ * i k (z) conduction coefficients degenerate in κ k i (z) because there are no curvatures involved (H α = H β = 1).

3D Geometrical and Constitutive Relations
The geometrical equations used for this thermo-mechanical model consider both the strains linked to the displacements u k and the strains linked to the sovra-temperature θ k .The equations can be written in matrix form as: where k is the 6 × 1 strain vector, ∆(z) is a 6 × 3 matrix containing the differential terms for the shell configuration, G(z) is a 6 × 3 matrix that includes the pure geometrical curvature terms, u k is the 3 × 1 displacement vector, and µ(z) k is the 6 × 1 vector containing the thermal expansion coefficients evaluated in α, β, and z directions.The explicit form of these matrices and vectors are here presented: The superscript k indicates that the matrix and the vector are valid for each physical layer.In order to obtain these coefficients, a rotation from the material reference system (1,2,3) to the structural reference system (α, β, and z) has been employed.
Thanks to the constitutive equations, the strain components can be related to the six stress components The constitutive relation is the well-known Hooke law: where the vector of mechanical strains k is the algebraic summation of geometrical strains and thermal strains.The elastic coefficient matrix C k (z) has 6 × 6 dimension and it depends on z in the case of a k layer constituted of an FGM.The elastic coefficient matrix used in this formulation has angle) to solve, in closed form, the mathematical problem with respect to the unknown displacements and temperature amplitudes by means of the Navier methodology in the α and β plane directions and the exponential matrix method in the z direction.
Using the constitutive relation expressed in Equation ( 11) and substituting Equation ( 9) in it: where M k (z) has 6 × 3 dimension and it indicates the pure mechanical coefficients; λ k (z) denotes the presence of the thermo-mechanical coupling coefficients in the structural reference system.Vector λ k (z) is defined as: According to the previously described substituting steps, the stresses are linked with u, v, w, and θ, which are the unknown variables of the formulation.

Exponential Matrix Methodology and Layer Wise Approach
In order to solve the problem in an exact and closed form, there is the necessity of the harmonic form for the primary unknowns u, v, w, and θ: Thanks to these impositions, the boundary conditions for all the structures are the simply supported sides.The amplitudes of the unknowns are identified with the appropriate capital letter (e.g., the amplitude of the displacement u is written as U).The ᾱ and β coefficients depend on the in-plane dimensions a and b and on the half-wave numbers m and n in the in-plane directions as follows: Substituting the unknowns written in harmonic form (Equations ( 14)-( 17)), the geometrical relations (Equation ( 9)) and the constitutive relations (Equation ( 11)) into the equilibrium relations (Equations ( 1)-( 4)), four second-order differential equations are written as follows: where coefficients A k s (z) (s rises from 1 to 19) and J k r (z) (r varies from 1 to 7) are not constant (they depend on z).Note that Equations ( 19)-( 22) are written in a generic k-th physical layer.
The unknowns, from this moment forward, are the amplitudes of the displacements and sovra-temperature and the related first-order derivatives (only created along the z direction because the ones created in in-plane directions are exactly evaluated and they become constant values).All four equations have coefficients that depend on z because of the curvature terms H α (z) and H β (z) (see Equations ( 5)-( 7)) and some FGM layers.In order to transform this set of non-constant coefficient second-order differential equations into a constant coefficients set, each k physical layer is split up into several fictitious layers.So, a number of M fictitious layers are employed along the thickness direction of the structure to discretise curvature terms and FGM properties.To define these fictitious layers, a new index has to be introduced: j is the index that counts these layers and it varies from 1 to M (total number of fictitious layers employed in the thickness direction).M can be easily calculated with the relation M = s • k where s indicates the number of subdivisions of each k physical layer.In these fictitious layers, the coefficients H α (z) and H β (z) and FGM properties can be exactly calculated in their middle points and the obtained values are constant within each fictitious layer.From this point forward, all the equations proposed are written in the j fictitious layer and they have constant coefficients.After the introduction of the fictitious layers, the system expressed in Equations ( 19)-( 22) can be transformed into a first-order system by redoubling the number of unknowns (as stated in [68,69]).After the redoubling of the unknowns, the set of first-order differential equations can be reported in compact form as: where [U j , V j , W j , Θ j , U j , V j , W j , and Θ j ] T is the new unknown vector containing both the displacement and sovra-temperature amplitudes and the related first-order derivatives along the z direction (indicated with the apex ).Different from the model proposed by Brischetto and Torre in [70], the sovra-temperature profile along the z direction is computed without using an external tool.It is calculated as a primary variable by means of Equation ( 4) combined with the equilibrium equations.Equation ( 23) can be written in a compact form as: A further rewriting of Equation ( 24) is: with A * j = D j −1 A j .X j is the 8 × 1 vector containing the unknowns and X j includes the derivatives along z of these unknowns.To solve the problem written in Equation ( 25), it is possible to use the method of the exponential matrix proposed in [68,69].The solution of the problem in Equation ( 25) can be written as: Knowing the exponential matrix term e (A * j zj ) , it is possible to obtain the unknown vector X j t (corresponding to the unknown vector at the top of the j-th fictitious layer) related with the unknown vector X j b (bottom of the j-th fictitious layer).The e (A * j zj ) term must be calculated introducing the thickness value h j of each j layer.The exponential matrix can be expanded in a power series and it must be computed for each fictitious layer for the values h j : where I is the 8 × 8 identity matrix.Equation ( 26) links the top with the bottom within each j fictitious layer.To link the top of the j fictitious layer with the bottom of the j + 1 fictitious layer, proper interlaminar continuity conditions must be imposed.These interlaminar continuity conditions must be imposed to u, v, and w displacements, sovra-temperature θ, transverse shear, and transverse normal stresses σ αz , σ βz , σ zz , and transverse normal heat flux q z at each fictitious layers' interface.The continuity of all these terms can be written in matrix form as: where Equation ( 28) is related to the displacements and sovra-temperature, Equation ( 29) is related to stresses and Equation (30) to the heat flux along the z direction.The conditions expressed in Equation ( 28) can be easily imposed for the amplitudes U j , V j , W j , and Θ j .It is possible to derive an amplitude form of Equations ( 28)-( 30) using the constitutive Equation ( 11) and the harmonic forms in Equations ( 14)- (17).These conditions can be rewritten in matrix form as: where the diagonal submatrix including 1 is Equation (28) written in matrix form.The other coefficients T i and τ i are Equations ( 29) and ( 30) in matrix form.A compact form of Equation ( 31) is: where T j+1,j is the transfer matrix.All the analyses can be conducted using this mathematical formulation considering the simply supported boundary conditions.This constraint configuration is automatically assured by Equations ( 14)- (17).It is possible to write: The load boundary conditions must be enforced at the outer faces of the structure; it is possible to write them as: Thanks to Equation (35), the sovra-temperature is directly imposed at the outer faces of the structure.Equation ( 35) can be rewritten in a compact way as: Equations ( 36) and ( 37) can be further compacted as: where vectors P t and P b contain the impositions related to the mechanical load conditions and the sovra-temperature.For the present thermal stress analysis, the mechanical load conditions are enforced to zero.
In order to include Equations ( 38) and ( 39) into an algebraic system written in compact form, it is possible to write X M t in terms of X 1 b to link the top of the last fictitious layer with the bottom of the first fictitious layer.This operation can be achieved by recursively introducing Equation (32) into Equation ( 26) as follows: Equation ( 40) defines the 8 × 8 matrix H m for structures embedding FGM layers.This matrix has a different size with respect to the H m matrix proposed in [63] for the pure mechanical case developed by Brischetto.The matrix size difference is due to the coupling of the 3D Fourier heat conduction equation directly with the elastic equilibrium relations of the problem.Introducing Equation (40) in Equation ( 38), this last one can be rewritten in terms of X 1 b as: Equations ( 39) and ( 41) can be now compacted as: The main characteristic of the final system written in Equation ( 42) is the fact that matrix E independently has an 8 × 8 dimension by the number of fictitious layers employed, even if the method uses a layer-wise approach.This new formulation can be considered as the generalization of the pure mechanical model proposed in [63] by Brischetto.Vector P now contains all the load impositions, both mechanical and thermal ones.The system in Equation ( 42) is formally the same as shown in [63,70], and in [71], but the addition of the 3D Fourier heat conduction equation to the 3D equilibrium relations is now considered.
This formulation can be implemented in a Matlab code where stresses, strains, and displacements can be evaluated along the thickness direction z of several structures embedding different FGM configurations.Once the displacements at the bottom of the first fictitious layer have been calculated, Equations ( 26) and (32) can be progressively used to compute the displacements and sovra-temperature (and related derivatives with respect to z) through all points in the z direction of the structure.

Results
This section is related to the comparison of results between the presented coupled thermo-elastic model and the past uncoupled thermo-elastic model proposed by Brischetto and Torre [70][71][72][73].The section is divided into two different parts: in the first one, two preliminary assessments are presented to validate the proposed general 3D exact coupled thermo-elastic shell model and to clearly understand the proper choice of N (order of expansion for the exponential matrix in Equation ( 27)) and the appropriate number of M (mathematical layers) for the calculation of FGM properties and constant curvature terms related to shell geometries.In the second part, four new benchmarks are proposed; after the opportune choice of N and M parameters, the present 3D coupled thermo-elastic model is considered validated and it allows for the discussion of effects due to the geometry of structures, thickness ratios, FGM laws, and temperature impositions for several cases.All the presented assessments and benchmarks consider an FGM layer with two constituent phases: a metallic one (Monel 70Ni-30Cu) and a ceramic one (Zirconia).These two constituent phases have the following elastic and thermal properties: where K m and K c indicate the bulk modulus of the metallic and ceramic phase, G m and G c indicate the shear modulus of the metallic and ceramic phase, µ m and µ c indicate the thermal expansion coefficient of the metallic and ceramic phase, and k m and k c indicate the conductivity coefficient of the metallic and ceramic phase.For each presented case (both assessments and benchmarks), the volume fraction of the ceramic phase V c follows the indicated law: where zFGM is the local thickness coordinate for the FGM layer (it goes from 0 at the bottom to h FGM at the top of the FGM layer), h FGM is the thickness of the FGM layer, and p is the related exponential coefficient.The FGM layer used for the proposed results is a full metallic at the bottom and full ceramic at the top.This consideration is the same for both sandwich and one-layered configurations.The bulk and shear moduli along the thickness direction depend on the volume fraction of the ceramic phase V c .These two material properties are estimated using the Mori-Tanaka model as follows: where . The heat conduction coefficient k is a function of the volume fraction of the ceramic phase V c as: and the thermal expansion coefficient µ can be computed using: In Equation ( 48), the dependence on V c is delegated to the bulk modulus K.The material data employed in this section were proposed by Reddy and Cheng in [72].

Preliminary Assessments
The new proposed 3D coupled thermo-elastic exact solution for shells embedding FGM layers, here indicated as 3D-u-θ, is validated using two preliminary assessments.A square one-layered FGM plate and a one-layered FGM cylindrical shell are investigated considering different thickness ratios.After these assessments, the 3D full coupled shell model will be defined as validated: the validation occurs for N = 3 (order of expansion for the exponential matrix) and M = 300 mathematical layers.For the cases presented in this section, the convergence of the 3D-u-θ model is obtained even for fewer mathematical layers than those used for the 3D uncoupled thermo-elastic model presented in [70] and called 3D(θ c , 3D) in Tables 1 and 2. However, higher values for N and M are chosen in a conservative sense.The results for these preliminary assessments are provided in non-dimensional forms as: where C = 10 −6 is a constant value, a is the in-plane dimension of the structure in the α direction, and K * = 10 9 Pa.

3D(θ c ,3D) [70]
3D-u-θ The first preliminary assessment shows a simply supported square plate (a = b) with thickness h = 1 m and bi-sinusoidal (m = n = 1) sovra-temperature imposed at the top and bottom surfaces (Θ t = +1 K and Θ b = −1 K).The structure is composed of a single FGM layer with the top part constituted of a ceramic phase and the bottom a metallic phase, as it is defined in Equations ( 43)- (48).The volume fraction of the ceramic phase V c is a function of the thickness coordinate, with the exponential coefficient p = 2 in Equation (45).The reference results are based on the 3D uncoupled thermo-elastic solution proposed by Brischetto and Torre [70] where the temperature is calculated by separately solving the 3D Fourier heat conduction Equation (3D(θ c , 3D)) by means of hyperbolic functions (see the mathematical formulation proposed by the authors in [70]).The present new solution uses an order of expansion N = 3 for the calculation of the exponential matrix and M = 300 mathematical layers for the calculation of curvature terms and FGM properties.The novelty of the 3D-u-θ model is the inclusion of the 3D Fourier heat conduction equation directly in the system including the 3D elastic equilibrium equations.This feature permits to take into account both the material layer and the thickness layer effects without the calculating the temperature profile using an external tool.Table 1 shows no-dimensional transverse and in-plane displacements and no-dimensional in-plane normal and transverse shear/normal stresses for different thickness ratios a/h.The present 3D coupled model provides the same results as the reference solution [70] for each thickness ratio a/h because, in both models, the material and thickness layer effects are properly evaluated.
The second preliminary assessment shows a simply supported cylindrical shell with radii of curvature R α = ∞ and R β = 10 m, bi-sinusoidal (m = n = 1) sovra-temperature imposed at the top and bottom surfaces (Θ t = +1 K and Θ b = 0 K), and in-plane dimensions a = 1 m and b = π 3 R β .The cylindrical shell is composed of a single FGM layer whose top part is constituted of a ceramic phase and the bottom of a metallic phase, as it is defined in Equations ( 43)- (48).The power law for the volume fraction of the ceramic phase V c is quadratic (p = 2) and it depends on the thickness coordinate zFGM .The results used as a reference are based on the previously discussed uncoupled 3D thermo-elastic solution proposed by Brischetto and Torre in [70].This second assessment uses the same values of N and M seen in the first one.The new proposed 3D coupled shell model (3D-u-θ) includes the 3D Fourier heat conduction equation in the same way seen in the first assessment.Table 2 shows no-dimensional transverse and in-plane displacements and no-dimensional in-plane and transverse stresses for two different thickness ratios (R β /h = 50 and R β /h = 1000).The 3D-u-θ model is always coincident with the reference solution [70] because both use the 3D Fourier heat conduction equation: coupled with the 3D elastic equilibrium Equations (3D-u-θ model) or separately solved using an external tool (3D(θ c ,3D) model).The results confirm all the previous considerations.
The benchmarks proposed in the next section consider the cases where different geometries, thickness ratios, FGM configurations, and temperature profiles are deeply evaluated.For this purpose, M = 300 mathematical layers combined with an N = 3 order of expansion will be employed in these new benchmarks, as suggested by the validation here conducted using the preliminary assessments.

New Benchmarks
Here, four new benchmarks are proposed considering plates, cylinders, cylindrical shells, and spherical shells.The involved geometries can be seen in Figure 1.Several sovra-temperature impositions, different m and n half-wave numbers and different FGM laws are considered for each case.The variation of the FGM law occurs with the variation of the parameter p in the exponential law of Equation (45).For each possible geometry, a FGM layer is involved (both single-layer and sandwich configurations).In all the proposed 3D coupled results, the N = 3 order of expansion for the exponential matrix and M = 300 mathematical layers are used.The new 3D coupled exact thermo-elastic model (3D-u-θ) will be used for these benchmarks and it will be compared with previously uncoupled 3D results where the 3D Fourier heat conduction equation was separately solved (3D(θ c ,3D)), the 1D Fourier heat conduction equation was separately solved (3D(θ c ,1D)), and the assumed linear temperature profiles were considered a priori (3D(θ a )).The results presented in this subsection can be useful for those scientists involved in the development of 2D and 3D analytical and numerical models for the thermal stress analysis of plates and shells embedding FGM layers.
The first benchmark proposes a simply supported sandwich square plate (a = b = 10 m) with an FGM core (see Figure 1).The thickness of the entire structure is h and the analysed thickness ratios are a/h = 2, 5, 10, 20, 50, 100.The thickness of the top and bottom skins is h 1 = h 3 = 0.2 h and the FGM core has h 2 = 0.6 h.The sovra-temperature has an amplitude value at the top of Θ t = +0.5 K and at the bottom Θ b = −0.5 K with bi-sinusoidal form (half-wave numbers m = 1 and n = 1).The volume fraction law of the ceramic phase V c for the FGM core has exponential p = 2; all the mechanical and thermal property variations of the FGM core through the thickness are described by means of Equations ( 46)- (48).The bottom skin is full metallic and the top skin is full ceramic, while the FGM core has a continuous variation (see Figure 2).Figure 3 shows the temperature profile through the thickness of a moderately thick (a/h = 5) and a moderately thin (a/h = 20) plate.In the case of the moderately thin configuration (a/h = 20), the calculated temperature profiles using 3D Fourier heat conduction Equation (θ c ,3D), using 1D Fourier heat conduction Equation (θ c ,1D), and using the coupled model (u-θ) are perfectly coincident.The assumed linear temperature profile (θ a ) is inappropriate for this case.For the moderately thick plate (a/h = 5), the temperature profile calculated using the full coupled model (u-θ) is coincident with the uncoupled model that computes the temperature profile via the 3D Fourier heat conduction Equation (θ c ,3D).These two models introduce both the thickness and material layer effects.The temperature profile using the 1D Fourier heat conduction Equation (θ c ,1D) does not consider the thickness layer effect.The uncoupled model that considers the temperature profile using an a priori linear assumption (θ a ) does not take into account both the effects.Table 3 shows the in-plane and transverse displacement components and the in-plane normal, in-plane shear, transverse shear, and transverse normal stress components in different positions through the thickness for different a/h ratios.In the case of very thin plates (a/h = 100), the last 3D models presented in Table 1 are coincident.For thick or moderately thick plates, the 3D(θ c ,3D) and the 3D-u-θ model provide the same results because they properly consider the thickness layer and material layer effects thanks to the implementation of the 3D Fourier heat conduction equation.The 3D(θ a ) model provides different results because it did not properly consider any material layer and thickness layer effect.Figure 4 shows the in-plane and transverse displacement components and two stress components (σ zz and σ αz ) through the thickness of a moderately thick plate (a/h = 10).The in-plane displacement and transverse displacement are continuous through the thickness direction because the mechanical properties of the FGM layer vary with continuity along the thickness direction and the compatibility conditions are correctly introduced in the shell model.The transverse normal stress and the transverse shear stress satisfies the boundary load conditions imposed at the top and at the bottom external surfaces of the structures (σ t zz = σ b zz = P z = 0 and σ t βz = σ b βz = P β = 0) and they are continuous because of the correct impositions of equilibrium conditions.Table 3.First benchmark, simply supported sandwich square plate with FGM (p = 0.5) core.Sovratemperature imposed as Θ t = +0.5 K and Θ b = −0.5 K for m = 1 and n = 1.3D uncoupled thermoelastic models from [70].The new 3D coupled thermoelastic solution is 3D-u-θ.σ αβ [10 3 Pa] at (α = 0, β = 0, z = h/4) 3D(θ a ) [70] 292 For the case related to the second benchmark, a simply supported one-layered FGM cylinder (see Figure 1) is analysed.The radii of curvature of the structure are R α = 10 m and R β = ∞.The global thickness of the structure is h.The in-plane dimensions are a = 2πR α and b = 30 m.The applied sovra-temperatures at the external surfaces are the same as those already proposed for the first benchmark, but half-wave numbers m = 2 and n = 1 are now imposed.The volume fraction of the ceramic phase V c is linear (p = 1), the bottom of the structure is full metallic and the top is full ceramic as can be seen in Figure 5.The reference equations for all the material characteristics are proposed in Equations ( 46)- (48).The temperature profiles through the thickness proposed in Figure 6 do not change when the thickness ratio varies.This feature is due to the symmetry and rigidity of the cylinder.Therefore, the effect of the thickness is negligible even if the structure is really thick (R α /h = 5 case).Table 4 shows the results in terms of displacements and stresses for different R α /h ratios.For thicker cylinders, the 3D(θ c ,3D) and 3D-u-θ models are in accordance because they take into account both the material layer and the thickness layer effects.The 3D(θ c ,1D) only considers the material layer effect and 3D(θ a ) only considers a linear assumed temperature amplitude: for this reason, the results are not correct.For thin cylinders, only the 3D(θ a ) model shows important differences with respect to the other three models (3D(θ c ,1D), 3D(θ c ,3D), and 3D-u-θ.Figure 7 provides the displacements and stresses for a moderately thick (a/h = 10) cylinder with one FGM layer.In-plane and transverse displacements are continuous because the mechanical and thermal properties of the FGM layer vary continuously in the thickness direction and the compatibility conditions are correctly imposed for each mathematical interface.The same consideration is also valid for the stresses σ αα and σ βz shown in Figure 7.The transverse shear stress satisfies the free mechanical load conditions at the external surfaces (σ t βz = σ b βz = P β = 0).σ αα [10 3 Pa] at (α = a/2, β = b/2, z = 0) 3D(θ a ) [70] 2245.9 σ αβ [10 3 Pa] at (α = 0, β = 0, z = h/4) 3D(θ a ) [70] 165 σ βz [10 3 Pa] at (α = a/2, β = 0, z = 3h/4) 3D(θ a ) [70] −104.The third benchmark presents a simply supported sandwich cylindrical shell with an FGM core (see Figure 1).The radii of curvature are R α = 10 m and R β = ∞ and the in-plane dimensions are a = π 3 R α and b = 30 m.The FGM core has thickness h 2 = 0.6 h and the two skins have h 1 = h 3 = 0.2 h; h is the global thickness of the structure.The material properties of this configuration are based on Equations ( 46)-( 48) with a quadratic exponential coefficient (p = 2).The bottom skin is full metallic and the top skin is full ceramic with thermal and mechanical properties expressed as in Equations ( 43) and (44) (as is visible in Figure 8).The investigated thickness ratios are R α /h = 2, 5, 10, 20, 50, The sovra-temperature in harmonic form (half-wave m = 1 and n = 1) has amplitude values at the top Θ t = −1 K and at the bottom Θ b = 0 K.The temperature profiles are provided in Figure 9 for a moderately thick (R α /h = 5) and a moderately thin (R α /h = 20) cylindrical shell.For thicker configurations, the thickness layer effect is a little bit more evident than the second benchmark because of the rigidity of the previous closed and symmetric cylinder.For thinner configurations, 3D(θ c ,1D), 3D(θ c ,3D), and 3D-u-θ models have no differences in terms of their temperature profiles.The 3D(θ a ) model always proposes an inadequate temperature profile because it did not consider any material layer and thickness layer effect.The same considerations are still valid for the discussion of displacements and stresses in Table 5 where the 3D(θ c ,1D), 3D(θ c ,3D), and 3D-u-θ models are quite coincident for all the R α /h ratios proposed.The great difference involves only the 3D(θ a ) model that is always inadequate for this configuration.Figure 10 shows displacement and stress evaluations through the thickness of the simply supported sandwich cylindrical shell with an FGM core.The analysed cylindrical shell is moderately thick and both proposed displacements and stresses are continuous because the FGM properties continuously vary along the thickness direction.σ αβ [10 3 Pa] at (α = 0, β = 0, z = h/4) 3D(θ a ) [70] 211 σ βz [10 3 Pa] at (α = a/2, β = 0, z = 3h/4) 3D(θ a ) [70] − The last benchmark takes into account a simply supported spherical shell with one FGM layer (see Figure 1).The global thickness is h, the radii of curvature are R α = R β = 10 m, the imposed sovra-temperatures are Θ t = 0 K and Θ b = −1 K, the in-plane dimensions are a = π 3 R α = b = π 3 R β , and the half-wave numbers are m = 2 and n = 1.The configuration and material characteristics are the same as seen in the second benchmark (see Equations ( 46)-( 48) for the bulk modulus, shear modulus, thermal expansion coefficient, and conductivity coefficient and see Equation (45) and Figure 11 for the volume fraction of the ceramic phase V c with p = 2).In Figure 12, the temperature profile for the four models are shown.The 3D(θ c ,3D) and 3D-u-θ models are mandatory for the correct analysis of the thick shells.The 3D(θ c ,1D) model is a good approximation only for thin spherical shells, but it always denotes some differences with respect to the previous two models.The 3D(θ a ) model is always inappropriate for this benchmark for each R α /h ratio considered.For the spherical shell, the thick configuration shows an important thickness layer effect.As in the previous benchmark, the 3D-u-θ model provides the same results obtained with the 3D(θ c ,3D) model because they consider both the thickness layer and the material layer effects.The displacements and stresses for several R α /h are proposed in Table 6 where it is evident how the results obtained using the 3D(θ a ) model are inadequate for each thickness ratio because the actual temperature profile is never linear along the thickness direction.In Figure 13, the displacement and stress evaluations through the thickness a spherical shell with one FGM layer are shown.In-plane and transverse displacements, transverse shear, and transverse normal stress are continuous because the compatibility and equilibrium conditions have been correctly imposed and also because the FGM layer continuously varies its own mechanical and thermal properties along the thickness direction.Transverse normal stress σ zz satisfies the load boundary conditions (σ t zz = σ b zz = P z = 0) as can be seen in Figure 13.σ αβ [10 3 Pa] at (α = 0, β = 0, z = h/4) 3D(θ a ) [70] 827 σ βz [10 3 Pa] at (α = a/2, β = 0, z = 3h/4) 3D(θ a ) [70] −184.

Conclusions
A full coupled thermo-elastic 3D exact shell solution for thermal stress analysis of shells and plates embedding FGM layers has been proposed.The sovra-temperature amplitudes have been directly imposed at the external surfaces in steady-state conditions and the sovra-temperature profile was evaluated along the z direction.The sovra-temperature profile is a primary unknown variable similar to the displacements; this is possible because the 3D Fourier heat conduction equation and the 3D equilibrium equations for shells are put together into a set of four second-order differential equations.This temperature profile considers both the thickness layer and material layer effects for each possible geometry, without separately solving the related 3D Fourier heat conduction equation.The set of four second-order differential equations for shells is solved in a closed-form thanks to the Navier solution and the exponential matrix method.Different analyses, in terms of displacements, temperature profiles, in-plane, and out-of-plane stresses have been shown for several thickness ratios, geometries, FGM configurations, and temperature impositions.The proposed results showed a complete match between the model that separately solves the 3D Fourier heat conduction equation and the present 3D full coupled thermo-elastic model.For each investigated variable (temperatures, displacements and stresses) and for each thickness ratio, geometry, FGM configuration, and load imposition, the results proposed differences that were always less than 0.5%.This new coupled method permits taking into account both the material and the thickness layer effects using a simpler and more consistent mathematical formulation.Moreover, a reduced number of fictitious layers

5 Figure 2 .Figure 3 .Figure 4 .
Figure 2. First benchmark, volume fraction of the ceramic phase for a simply supported sandwich square plate with FGM (p = 0.5) core.

1 Figure 5 .Figure 6 .
Figure 5. Second benchmark, volume fraction of the ceramic phase for a simply supported cylinder with one FGM (p = 1) layer.

Figure 8 .Figure 9 .
Figure 8. Third benchmark, volume fraction of the ceramic phase for a simply supported sandwich cylindrical shell with FGM (p = 2) core.

2 Figure 11 .
Figure 11.Fourth benchmark, volume fraction of the ceramic phase for a simply supported spherical shells with one FGM (p = 2) layer.

Figure 12 .
Figure 12.Fourth benchmark, temperature profiles for thick and moderately thick simply supported spherical shell with one FGM (p = 2) layer.