A coupled hygro-elastic 3D model for steady-state analysis of functionally graded plates and shells

: This 3D coupled hygro-elastic model proposes the three-dimensional (3D) equilibrium equations associated with the 3D Fick di ﬀ usion equation for spherical shells. The primary unknowns of the problem are the displacements and the moisture content. This coupled 3D exact shell model allows to understand the e ﬀ ects of the moisture ﬁ eld in relation with the elastic ﬁ eld on stresses and deformations in di ﬀ erent plates and shells. This model is speci ﬁ cally developed for con ﬁ gurations including functionally graded material (FGM) layers. Four di ﬀ erent geometries are analyzed using an orthogonal mixed curvilinear reference system. The main advantage of this reference system for spherical shells is the degeneration of the equations to those for simpler geometries. The solving method is the exponential matrix method in the thickness direction. The closed-form solution is possible because of simply supported sides and harmonic forms for displacements and moisture content. The moisture content amplitudes are directly applied at the top and bottom outer faces through steady-state hypotheses. The ﬁ nal system is based on a set of coupled homogeneous second-order differential equations. The moisture ﬁ eld e ﬀ ects are evaluated for the static analysis in terms of displacement, strain, and stress components. After preliminary validations, used to better understand how to properly de ﬁ ne the calculation of the curvature-related terms and FGM properties, four new benchmarks are proposed for several thickness ratios, geometrical data, FGM con ﬁ gurations, and moisture values imposed at the external surfaces. From the results, it is clear the accordance between the uncoupled hygro-elastic model and this new coupled hygro-elastic model when the 3D Fick di ﬀ usion law is employed. Both e ﬀ ects connected with the thickness layer and the embedded material are included in the 3D hygro-elastic analyses proposed. The 3D coupled hygro-elastic model is simpler than the uncoupled one because the 3D Fick di ﬀ usion law does not have to be separately solved.


Introduction
In the framework of aerospace structural analyses, the external environmental conditions become extremely important to be analyzed.Severe temperature and humidity gradients are the classic conditions where these aircraft and spacecraft usually operate [1][2][3][4].In particular, modern aerospace structures can have problems concerning the absorption of humidity from the external air.This phenomenon can create irreversible damages related to the debonding and the mechanical degradation of composite materials, with severe problems related to the safety of aircraft structures.For these reasons, the use of FGMs, whose peculiarity is the continuous variation of its own mechanical, thermal and hygroscopic properties along a given direction, can be relevant.These peculiarities of the FGMs can strongly mitigate all the aforementioned possible damages typical of composite materials.
The scientific and technical literature fulfills various features related to hygro-elastic models for plates and shells embedding FGM layers.Numerical and analytical models (1D, 2D, and 3D approaches) are deeply discussed in this section to better understand the originality of the new proposed model.In all the discussed papers, the moisture field is involved.Two different sections are proposed: in the first one, all the numerical models are deeply described; meanwhile, in the second one, analytical models are shown.
In the context of numerical models, Aria and Friswell [5] showed a hygro-thermal investigation of functionally graded (FG) sandwich microbeams based on Eringen nonlocal elasticity theory where motion equations are derived using Hamilton principle and the first order shear deformation theory (FSDT).The finite element (FE) procedure is based on a five-noded beam element.In the study by Garg et al. [6], the bending analysis of sandwich FGM beams under combined hygrothermomechanical loadings was evaluated.The study was carried out using a FE based on a high order zigzag theory (HOZT).Li and Tang [7] proposed three-directional FGMs to carry out the composition of slender beams.The nonlinear governing equations were developed using the principle of minimum for potential energy.Liu and Mohammadi [8] proposed the bending response of FG nanobeams under hygrothermomechanical loadings.The nonlocal strain gradient theory (NSGT) was implemented to consider both hardening and softening impacts on the system response.Nguyen et al. [9] discussed the hygro-thermal effects on vibration and buckling analyses of FG beams.This work was based on equations of motion obtained from Lagrange equations using a higherorder shear deformation theory.Tang and Ding [10] showed the nonlinear hygro-thermal dynamics of a bidirectionally FG beam where material and hygro-thermal properties gradually change along both thickness and length directions.Dastjerdi et al. [11] proposed a quasi three-dimensional approach for the bending analysis of plates embedding FGMs with porosities.The bending equations were developed using Hamilton principle, and the solution was obtained using a semi-analytical polynomial method.Lee and Kim [12] showed the postbuckling behaviors of FGM plates in hygrothermal environments.The model was based on the FSDT and the von Karman strain-displacement relations, and it was developed using the FE method and the Newton-Raphson technique.Sobhy [13] gave out the bending analysis of sandwich-curved beams with graphene platelets/aluminum (GPLs/Al) nanocomposite face sheets and aluminum honeycomb.The differential quadrature method (DQM), based on a new shear and normal deformation curved beam theory, was employed.Sobhy also [14] presented a three-dimensional asymmetric bending analysis of solid circular and annular plates, made of aluminum matrix reinforced with uniformly distributed GPLs, lying on an elastic medium with different boundary conditions.In-plane magnetic field, transverse external loads, thermal loads, and humid conditions were included.Zhao et al. [15] presented the dynamic analysis of GPL sandwich shallow shells reinforced with FG porous and embedding shape memory alloy (SMA) wires under hygrothermal loadings with general boundary conditions.Dai et al. [16] showed a mechanical model to investigate the coupling between temperature and moisture for a porous FG-carbon nanotube reinforced composite rotating annular plate with variable thickness.Numerical results were achieved by combining the DQM, the Runge-Kutta method, and the Newmark method.Saadaftar and Aghaie-Khafri [17] proposed a rotating FGM cylindrical shell with imperfect surfaces and piezoelectric layers subjected to an axisymmetric hygro-thermo-electromechanical loading.The Fourier series expansion method through the longitudinal direction and the DQM along the thickness direction were used to solve governing differential equations.Nie et al. [18] showed the hygro-thermo-electromechanical coupling problem of FG piezoelectric structures using the edge-based smoothed point interpolation method (ES-PIM).The basic equations for these structures were derived in the multi-physical field.The resolution of the coupled problem was obtained via ES-PIM and FE method under different hygrothermal conditions.
In the case of exact analytical models, Akbarzadeh and Chen [19] showed analytical solutions for hygrothermal stresses in 1D FG piezoelectric media where material properties varied through the thickness direction.The media were subjected to an external constant magnetic field.Governing equations were written in a unified form in cylindrical and spherical coordinates.Ebrahimi and Barati [20] presented an hygro-thermo-mechanical vibration analysis of FG size-dependent nanobeams exposed to various hygrothermal loadings.The analysis was performed via a semi-analytical differential transform method.Three hygro-thermal loadings types (uniform, linear, and sinusoidal) were investigated.Ebrahimi and Barati also [21] proposed a study where the combined effects of moisture and temperature on free vibration characteristics of FG nanobeams resting on elastic foundation were investigated.Different refined beam theories, including shear deformations, were involved.Jouneghani et al. [22] investigated the bending behavior of FG nanobeams with internal porosity and subjected to a hygro-thermomechanical loading.Eringen nonlocal theory was applied for the numerical study considering a uniform porosity within the nanobeam.The bending response of porous FG Bernoulli-Euler nanobeams under hygrothermomechanical loadings was analyzed by Penna et al. [23].The governing equations of the problem associated with the local/nonlocal stress and strain-driven gradient models of elasticity were derived using the principle of virtual work.Wang et al. [24] investigated the influence of hygrothermomechanical coupling loadings on buckling behaviors of porous bidirectional FG Timoshenko nanobeams.The governing equations and boundary conditions were derived using a two-phase local strain gradient theory/NSGT employing the principle of virtual work.A hygro-thermal stress analysis of rotating FG graphene/metal sandwich cylindrical shells embedding an auxetic honeycomb core was proposed the study by Allam et al. [25].The simply supported sandwich shell was spinning about its axial axis with a constant angular speed.Arshid et al. proposed in [26] a study about the vibrational behavior of three-layered cylindrical shells embedding graphene nanoplatelets reinforced composite core and piezoelectro-magnetic face sheets.The whole structure was rested on the visco-Pasternak foundation, and it was exposed to hygro-magneto-electro-thermal loads.Karimiasl et al. [27] investigated the post-buckling behaviors of doubly curved composite shells in a hygro-thermal environment employing multiple scales perturbation methods.Three-phase composite shells with polymer/carbon nanotube/fiber and polymer/GPL/fiber and SMA/matrix were taken into consideration according to the Halpin-Tsai model.Zidi et al. [28] investigated the bending response of an FGM plate resting on elastic foundation and subjected to hygrothermomechanical loadings.Zenkour and Radwan [29] investigated the effects of moisture and temperature on the bending analysis of FG porous plates resting on two-parameter elastic foundation.The effects of transverse normal and shear strains were taken into account.Tang et al. [30] proposed a study where the influence of the hygrothermal effects on center deflections, fundamental natural frequencies, and vibration deflection amplitudes for shells were discussed to understand the service life of polymer matrix composites during operating conditions.In the study by Mudhaffar et al. [31], the bending behavior of an advanced FG ceramic-metal plate subjected to a hygrothermomechanical load and resting on a viscoelastic foundation was studied using a simple higher-order integral shear deformation theory.
All the proposed exact models are related to 1D or 2D theories.For this reason, a 3D hygro-elastic model was developed by Brischetto and Torre in [32] to understand the behavior of shells embedding FGM layers when external moisture conditions were imposed.The present work can be seen as a general extension of the works presented by Brischetto in [33][34][35] where the mechanical and free vibration analyses of FGM shells were proposed.In the present article, the important novelty with respect to the 3D hygroelastic shell model proposed by Brischetto and Torre in [32,36] is the full coupling between the displacement field and the hygroscopic field.Therefore, the moisture content profile is directly obtained from the solution of the system, and it is not "a priori" defined.A similar procedure was already used by the authors in [37,38] in the case of full coupling between the displacement field and the thermal field for the 3D bending analysis of composite and FG plates and shells, and in the study by Brischetto and Cesare [39] for the case of full coupling between the displacement field and the hygroscopic field for the 3D bending analysis of classical multilayered composite plates and shells.The Fick diffusion equation [40,41] coupled with 3D equilibrium equations gives a system in mixed curvilinear orthogonal coordinates [42,43].The system is solved via the exponential matrix method [44,45].
This article is organized in a section containing the theoretical formulation of the 3D coupled hygro-elastic model, a section about results (divided into preliminary validations and new benchmarks) and a section about conclusions where the achieved results are summarized.

3D coupled hygro-elastic model for FGM shells
The 3D formulation for full coupled exact hygro-elastic analyses of FGM plates and shells is here presented.This formulation allows hygro-elastic analyses of several geometries such as cylinders, plates, cylindrical panels, and spherical shells.The differences between these four geometries are highlighted in Figure 1. Figure 1 indicates also the mixed curvilinear orthogonal reference system with the origin point located in the left corner.The directions α and β are along to the in-plane curved sides, and they lie on the surface Ω 0 , defined as the middle and reference surface.Ω 0 surface will be used for the calculation of all the geometrical parameters.The thickness direction z is normal to the Ω 0 surface and directed upside and towards the top surface.A coupled formulation is possible because the three displacement components u, v, and w, and the scalar moisture content are the four unknown primary variables of the problem.

Geometrical and constitutive equations
Geometrical relations written in matrix form are here proposed for the coupled 3D hygro-elastic model.The matrix equation is: where ε k is the × 6 1 hygro-mechanical strain vector, ε u k is the × 6 1 strain vector related to the × 3 1 displacement vector u k , ε k is the × 6 1 strain vector related to the scalar moisture content k , z Δ k ( ) is the × 6 3 matrix containing differential geometrical terms, G z k ( ) is the × 6 3 matrix containing algebraic geometrical and curvature terms, and η z k ( ) is the × 6 1 vector containing hygrometric expan- sion coefficients.
Each k FGM layer has elastic and hygroscopic proper- ties that are functions of the z coordinate.Isotropic FGM layers are analyzed.Therefore, material reference system (1, 2, 3) is coincident with the structural reference system (α β z , , ).The extended forms of the previous matrices and vectors are as follows: R α and R β are the radii of curvature in the α and β in- plane directions, respectively; they are constant values.

H z
α ( ) and H z β ( ) coefficients introduce the curvature shell terms.They are defined, for each direction, as follows: The curvature shell terms are linear functions of the thickness coordinate z or z ˜. h is the total thickness of the struc- ture and it is always considered as constant.The z coordi- nate varies in the range between − ∕ h 2 and ∕ h 2; z ˜is in the range between 0 and h.
The moisture content is defined as follows: where W d is the mass of the dry material in kg, W c is the mass of the moisture present in the material in kg, W is the mass of the moist material in kg (that is the sum of W d and W c ), V is the volume of the structure in m 3 , c is the moisture concentration in kg/m 3 , and ρ d is the density of the dry material in kg/m 3 .Therefore, is in nondimensional or in percentage form.
Stresses are linked with strains via the 3D Hooke law written for a generic k FGM layer: where the stress vector is ) and the elastic coefficient matrix C z k ( ) has × 6 6 dimension.The substitution of Eq. (1) into Eq.( 5) gives: where z k M ( ) indicates the pure mechanical coefficients and it has dimension × 6 3; ξ z k ( ) includes the hygro-elastic coupling coefficients in the structural reference system: In Eq. ( 6), the six stress components are linked with displacements u, v, w, and the moisture content that are unknown variables of this formulation.

3D equilibrium equations for spherical shells
The equilibrium equations and the 3D Fick diffusion equation for spherical shells are the basic equations of the present coupled formulation.The set of three equilibrium equations for spherical shells written in the matrix form using the mixed curvilinear orthogonal coordinates (α, β, z) is as follows: Coupled hygro-elastic 3D model for steady-state analysis of FG plates and shells  5 where is the vector whose stress components are applied on the face perpendicular to α axis; is the vector whose stress components are applied on the face perpendicular to β axis; is the vector whose stress components are applied on the face perpendicular to z axis; is the vector for the stress components that is the vector for the stress components that considers the R β correction.
The final compact form is: where subscripts α, β, and, z indicate the partial derivatives Eq. ( 9) is coupled with the 3D Fick diffusion equation for spherical shells in steady-state conditions: where the diffusion coefficients z * k 1 ( ), z * k 2 ( ) and z * k 3 ( ) for each k layer are defined as follows: ), and z k 3 ( ) are the three diffusion coefficients in the mixed curvilinear reference system depending on z in the case of FGM layers.In Eq. ( 10), the diffusion coefficients take into account the curvature of the geometry along in-plane directions.For the plate cases, the z * i k ( ) diffusion coefficients degenerate in z i k ( ) because curvatures are not involved.

Solution procedure
The resolution of this coupled hygro-elastic problem in a closed form needs the imposition of the Navier conditions for all the main unknowns u k , v k , w k , and k .These impositions can be written as follows: Using the Navier form for the unknowns, the boundary conditions are simply-supported sides for all the investigated structures.The capital letters proposed in Eq. ( 12) are the amplitudes of the unknowns.The α ¯and β ¯coefficients depend on the a and b dimensions and on the half-wave numbers m and n in the in-plane directions: All the primary variables expressed in Eq. ( 12) are 3D variables because they are functions of α, β, and z.Substituting the primary variables written in the harmonic forms as in Eq. ( 12), the geometrical relations as in Eq. ( 1) and the constitutive equations as in Eq. ( 5) into the governing Eqs. ( 9) and (10), a global set of four second-order differential equations is obtained.It can be properly transformed into a first-order differential equation system by redoubling the number of unknowns (as stated in [44] and [45]).After all these manipulations, the set of first-order differential equations can be written in a compact matrix form as follows: where Coupled hygro-elastic 3D model for steady-state analysis of FG plates and shells  7 X z k ( ) is the × 8 1 redoubled unknown vector containing displacements, moisture content amplitudes, and related first order derivatives along z direction (stated with the superscript ′ ).Eq. ( 14) is valid for each k physical layer.Coefficient matrices in Eq. ( 14) depend on z due to the curvature terms H z α ( ) and H z β ( ), and elastic and hygro- scopic terms for FGMs.To solve the system in Eq. ( 14), the coefficients must be constant.Introducing a number of L fictitious layers in each physical layer k along the thickness direction, H z α ( ), H z β ( ), D z k ( ), and A z k ( ) can be exactly calculated in the middle points of these layers, and the values become constant within each of it.A generic fictitious layer is indicated with j and goes from 1 to = × F L P (where P is the total number of physical layers).The system in Eq. ( 14) has now constant coefficients, and it can be rewritten as follows: where . The exponential matrix method, proposed in [44,45], is employed and the solution of the problem in Eq. ( 16) is: The unknown amplitude vector X t j at the top of the j ficti- tious layer can be obtained from the unknown amplitude vector X b j at the bottom of the same j fictitious layer.The exponential matrix is computed referring to the thickness value h j .The exponential matrix can be expanded and calculated in correspondence of h j for each fictitious layer: I is the × 8 8 identity matrix.Further relationships must be included for interlaminar continuity conditions.They are useful to link the top of the j fictitious layer and the bottom of the + j 1 fictitious layer.These interlaminar continuity conditions must be imposed for displacements u, v, and w, moisture content , and transverse shear/normal stresses σ αz , σ βz , and σ zz and transverse normal moisture flux g z in the z direction.The interlaminar continuity conditions in the compact form is expressed as follows: The interlaminar continuity conditions expressed in Eq. ( 19) can be easily imposed directly for the amplitudes U j , V j , W j , and j .It is possible to obtain an amplitude form of Eqs. ( 19)-( 21) taking into account the constitutive relations in Eq. ( 5) and the harmonic forms in Eq. ( 12).This amplitude form can be rewritten as follows: 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 where + T j j 1, is the transfer matrix that links the bottom of the + j 1 layer with the top of the j layer.The boundary conditions for simply supported sides are satisfied for harmonic forms: Load boundary conditions are imposed on the external faces as follows: The matrix forms of Eq. ( 25) are as follows:   (27) and they can be further compacted as follows: Superscript F indicates the last fictitious layer and superscript 1 is the first fictitious layer.Vectors t P and b P include the impositions related to the mechanical load in the three directions α, β, and z and the moisture content.Assuming a classical hygro-elastic stress analysis, the mechanical loads in α, β, and z directions are set equal to zero.
It is possible to write X t F in terms of X b 1 to directly link the top of the last fictitious layer with the bottom of the first one.This operation can be achieved introducing recursively Eq. ( 22) into Eq.( 17): Eq. ( 30) defines the × 8 8 matrix H m for structures embedding FGM layers.Introducing Eq. ( 30) into Eq. ( 28): Eqs. ( 29) and ( 31) can be now compacted in an unique equation as follows: Matrix E has × 8 8 dimension independently by the number of fictitious layers employed and layer-wise approach.Vector P contains all the load impositions, both mechanical and hygroscopic ones, at the top and at the bottom of the structure.The system in Eq. ( 32) is formally the same shown in the previous studies by Brischetto [32][33][34] even if the 3D Fick diffusion equation is added to equilibrium relations.This mathematical formulation is implemented in a Matlab code to evaluate stresses, strains, and displacements along the thickness direction for all the structures presented in the next section and including different FGM configurations (sandwich and single-layered ones).
Once the unknown vector at the bottom of the first fictitious layer ( = j 1) is computed, Eqs. ( 17) and ( 22) can be recursively applied to obtain the trends of all the variables through the z direction of the structures.
Coupled hygro-elastic 3D model for steady-state analysis of FG plates and shells  9  [GPa] 227.24G m [GPa]  65.55 K c [GPa]  125.83 G c [GPa]  58.077 Single FGM layer Single FGM layer -235.9 -236.8 -236.8 -236.9 -236.9 3D-u--235.9-236.9 -236.9 -236.9 -236.9 σ αα [Pa] 3D FEM [32] 3 Results This 3D coupled hygro-elastic shell model here developed is indicated with the acronym 3D-u-that clearly indicates the employed primary variables (displacement vector and scalar moisture content) .This new model is also compared with past 3D uncoupled hygro-elastic shell models where the moisture content profile is separately defined.These models are indicated as 3D( a ) when the moisture content profile is a priori assumed as linear through the entire thickness direction, 3D( c , 1D) when the moisture content profile is separately calculated by solving the 1D version of the Fick diffusion equation, and 3D( c , 3D) when the moisture content profile is separately calculated by solving the 3D version of the Fick diffusion equation.This section is splitted in two different parts: the first subsection is devoted to the validation of the 3D coupled hygro-elastic model for FGM shells by means of comparisons with the past 3D( c , 3D) model and the 3D solid finite element model by Patran and Nastran in the study by Brischetto and Torre [32].This assessment subsection is very important because it is useful to understand the number of fictitious layers F that must be used to prop- erly account for effects of the shell curvatures and FGM variable properties.In addition, the order of the expansion N for the exponential matrix must be also defined to obtain accurate results.In the second subsection, four new benchmarks are presented.They are useful to understand all the potentiality of the model because they consider different parameters such as the thickness ratio, the FGM configuration, the moisture content at the external surfaces, and the geometry.In both subsections, the FGM is composed by two phases: the metallic one is Monel (70Ni30Cu) and the ceramic one is Zirconia.The hygroelastic properties of these two materials are clearly shown in Table 1.
The volume fraction of the ceramic phase for all the cases proposed (both assessments and benchmarks) is: where z ˜FGM is the local thickness coordinate of the FGM layer, p is the exponential coefficient, and h FGM is the total dimension in the z direction of the FGM layer.The FGM layer is completely in Monel for = V 0 c , and it is completely in Zirconia for = V 1 c .Bulk and shear moduli of the FGM layer vary in the thickness direction following the Mori-Tanaka model: where ) .The moisture expansion coefficient η, written according to the Hatta and Taya model, is a function of the bulk modulus as follows: and the moisture diffusion coefficient can be computed using: where the dependence on V c is clear.The material data employed in this section were proposed by Reddy and Cheng in [46].

Comparisons for validations
Two preliminary validation results (PVRs) are proposed for this new three-dimensional general exact coupled hygro-   elastic shell model (3D-u-).A square plate and a cylindrical shell panel are analyzed for several thickness ratios, moisture content profiles and lamination schemes.For these validation cases, the 3D-u-model is compared with the 3D( c , 3D) model and the 3D FE model via Patran and Nastran (for the details about these models, see [32]).Table 2 shows all the data used for the assessment cases.When = N 3 and = F 300, the 3D-u-model obtains the same results as the 3D( c , 3D) presented by Brischetto and Torre in [32], but some differences occur with respect to the 3D FE analysis.The results obtained by the Patran and Nastran model are always the same for each number of fictitious layers F because the 3D solid FE model does not use fictitious layer discretization.
The first preliminary validation results (PVR1) propose a simply supported square plate composed by a single FGM layer.The hygro-elastic material properties are presented in Table 1, and the geometrical parameters are given in the first column of Table 2. Table 3 shows the results of the conducted analyses for different numbers of fictitious layers F and different thickness ratios ∕ a h.Transverse normal displacement w, in-plane displacement u, and in- plane normal stress σ αα are analyzed.There is a clear con- vergence of results for the increasing number of fictitious layers F .For the majority of the cases, a number of 200 fictitious layers is sufficient for both the 3D uncoupled model and for the 3D coupled model (both using the 3D Fick diffusion equation).
The second preliminary validation results (PVR2) show a simply supported cylindrical shell embedding a single FGM layer.The hygro-elastic material properties are collected in Table 1, and the geometrical parameters are given in the second column of Table 2. Table 4 proposes results  for different numbers of fictitious layers F and for two different thickness ratios: a thick configuration and a thin configuration.Transverse and in-plane displacements are analyzed and the reference solutions are the same seen in PVR1.The same conclusions seen for PVR1 are confirmed.For a number of fictitious layers = F 200 and order of expansion = N 3 for the exponential matrix, the 3D-u- model is coincident with the reference solution for both thickness ratios because these two different 3D analytical models take into account the same effects but in a different way.
After these PVRs, a number of = F 300 fictitious layers (that is overbundant with respect to = F 200) and an order of expansion for the exponential matrix = N 3 are fixed for all the future new benchmarks proposed in the next section.

New benchmarks
Four benchmark cases are here proposed.They take into account different characteristics for each case to have a global overview of all the peculiarities.The proposed four benchmarks consider the geometries given in Figure 1.The results obtained in this subsection consider different lamination schemes, FGM properties, and external moisture content impositions.= N 3 is used as the order of expansion for the exponential matrix, and = F 300 is the number of ficti- tious layers involved in each benchmark.This choice about N and F has been deeply discussed in the previous subsec- tion.The new 3D coupled hygro-elastic model (3D-u-) is compared with 3D uncoupled models presented by Table 7: Benchmark 2 (B2), closed cylinder embedding a single FGM layer.3D uncoupled hygro-elastic models [32] vs the 3D coupled hygroelastic model [32] 0.5000 0.5000 0.5000 0.5000 0.5000 3D( c , 1D) [32] 0.2627 0.2627 0.2627 0.2627 0.2627 3D( c , 3D) [32] 0 Brischetto and Torre in [32] (3D( c , 3D), 3D( c , 1D), and 3D ( a )).Now, some reference results have been corrected because of some residual errors in [32].All the geometrical features involved in the four benchmarks are collected in Table 5.The results discussed in this subsection can be used as test cases for scientists involved in the developing of 2D or 3D analytical/numerical models to evaluate the effects of the moisture content in these types of structures.The material data for the four benchmarks are given in Table 1.
The first benchmark (B1) considers a square plate with a single FGM layer and simply supported sides.The top of the structure is fully ceramic, and the bottom is fully metallic.All the data used for the hygroelastic analysis are presented in the first column of Table 5.Several thickness ratios are considered to evaluate the effects of geometrical parameters.The thickness ratios vary from ∕ = a h 2 to ∕ = a h 100.The moisture content is directly imposed at the external surfaces in harmonic form.In Figure 2, it is possible to see how the bulk modulus and the volume fraction of the ceramic phase vary across the thickness direction of the material.The trend adopted for the volume fraction of the ceramic phase is linear and p = 1 is the exponent used in Eq. (33).In Table 6, moisture content, displacements, stresses and strains are proposed.The results for the thick plate highlight the accordance between the 3D-u-and the 3D( c , 3D) because both consider the three-dimensional nature of the problem, taking into account both the material and the thickness layer effects.For this reason, these two solutions show the most accurate results.Due to the fact that the 3D( c , 1D) considers only the material effect and the 3D( a ) discards both material and thickness layer effects, the results proposed by these two models are not correct for the thick cases.For what concerns the thin configuration, as the thickness layer effect vanishes, even the 3D ( c , 1D) is reliable as the 3D-u-and the 3D( c , 3D); the only discrepancies are proposed by the 3D( a ) because it assumes a moisture content profile always linear.In Figure 3, it is possible to visualize the four moisture content profiles for a moderately thick case and a moderately thin one.It is possible to clearly see the differences involved for the moderately thick structure case among the ucoupled profile, the calculated ( c , 3D) profile and the calculated ( c , 1D) profile.For the moderately thin structure case, only the assumed linear profile ( a ) is wrong because it does not take into account both material and thickness layer effect.Figure 4 shows the trend along the z direction of two displacement components, the moisture content calculated with this new 3D coupled model, and stresses.All the variables are continuous along the thickness direction because this is a peculiarity of the FGMs, and it also demonstrates the correct implementation of the continuity conditions in the model.The transverse normal stress σ zz and the transverse shear stress σ βz satisfy the external mechanical loading boundary conditions: they equal zero at both the external faces because no external mechanical loads act on them.
The second benchmark (B2) is about a closed cylinder with a single FGM layer and simply supported sides.The metallic and ceramic phases and related mechanical and hygroelastic properties are the same seen in B1 and in Table 1.The geometrical data of the structure are shown  .The FGM law (Eq.( 33)) has exponent = p 2. In Figure 5, it is possible to see the bulk modulus and the volume fraction trend.In Table 7, results are shown to compare the 3D( c , 3D) model, the 3D( c , 1D) model, the 3D( a ) model, and the new 3D-u-θ model.The conclusions for both thin and thick structures are the same seen in the first benchmark.In Figure 6, the moisture content profile is given for a moderately thin and a moderately thick case.For the closed cylinder case, the differences among the ( c , 1D) profile, ( c , 3D) profile, and the uprofile are also small for the moderately thick case.Only the ( a ) profile is always incorrect because the hypothesis of linear profile is not true due to the material and thickness layer effects.Figure 7 shows some of interest to understand the main peculiarities of this case.All the variables presented are continuous: it means the correctness of the continuity conditions.The continuous trends of these variables are a specific characteristic of the FGMs.In addition, as no mechanical loads are applied, the stress σ zz is zero at the top and bottom surfaces.
The third benchmark (B3) proposes a sandwich cylindrical shell panel with simply supported sides.).For this sandwich configuration, the volume fraction of the ceramic phase V c of the core is 0 at the bottom (full metallic constituent) and 1 at the top (full ceramic constituent); it varies with continuity from 0 to 1 in the FGM core.The exponent p adopted for the FGM law is 0.5.The volume fraction V c and the bulk modulus trends   γ αz [10 8  -] at = + ∕ h z 3 ˜, = α 0, and = ∕ β b 2 3D( a ) [32] 0.000 0.000 0.000 0.000 0.0000 3D( c , 1D) [32] 0.000 0.000 0.000 0.000 0.0000 3D( c , 3D) [32] 0.000 0.000 0.000 0.000 0.0000 3D-u-0.0000.000 0.000 0.000 0.0000 be seen in Figure 8.The mechanical and hygroscopic properties are the same used for B1 and B2.The geometrical characteristics are shown in the third column of Table 5.Table 8 shows results for the four different models and for different ∕ R h α ratios (from thick to thin case).Even for the sandwich configuration, the thick cases ( ∕ = R h 2 α and ∕ = R h 4 α ) show the accordance between the 3D( c , 3D) model and the 3D-u-model.The reason is to be sought in the caption of the thickness effect that is crucial for thick structures.The 3D( c , 1D) cannot grasp the thickness layer effect but only the material layer effect that is not so much influent for these configurations.The 3D( a ) model is always inadequate for thick structures because the moisture content profile assumed as linear is not the actual one.For the thin structure case, as the thickness layer effect is less important, the 3D( c , 1D) model is similar to the 3D ( c , 3D) and the 3D-u-models.These conclusions can also be seen in Figure 9 where the moisture content profiles are computed for a moderately thick and a moderately thin configuration.In Figure 10, six variables are proposed: two displacements, the moisture content profile calculated with the 3D-u-model, and two stresses.Even for the sandwich case, all the variables are continuous as a consequence of the continuous variation in the thickness direction of the FGM properties.In addition, the transverse shear stress σ αz fulfills the external load conditions because no mechan- ical loads are applied at the external surfaces.
The last benchmark (B4) considers a sandwich spherical shell panel with simply supported sides.The sandwich lamination, the thickness values of the three layers, and the FGM law are the same seen in B3 (see also Figure 8).The hygro-elastic properties are given in Table 1, and the geometrical data are proposed in the fourth column of Table 5.Table 9 presents some results for different thickness ratios and for different variables.The considerations discussed for the B3 are the same.In Figure 11, it is possible to see the moisture content profile for a moderately thick and a moderately thin sandwich spherical shell configuration.It may be easily noted the fact that the ( c , 1D) profile has an important difference with respect to the ( c , 3D) and u-profiles because it does not take into account some important effects acting on this case.For what concerns the moderately thin case, the ( c , 1D) is sufficiently accurate.In the end, Figure 12 shows six variables in the thickness direction.The correct imposition of the continuity conditions are confirmed.The transverse shear stress σ βz and the trans- verse normal stress σ zz satisfy the load boundary conditions because no mechanical loads are applied at the external surfaces.

Conclusion
A full coupled hygro-elastic 3D exact shell model for the static analysis of FGM structures has been shown.The moisture content along the thickness direction is evaluated in steady-state conditions.The moisture content profile is a primary variable of the problem together with the displacement components.The 3D Fick diffusion equation for spherical shells allows a moisture content profile that takes into account both thickness layer and material layer effects.The coupled system is solved in a closed form using the exponential matrix method and a layer wise approach.Different results, in terms of displacements, stresses, and moisture content profiles, have been discussed for several thickness ratios, geometrical properties, FGM configurations (sandwich or single layer), and moisture content impositions.These analyses showed a very comforting match between the past 3D uncoupled model that separately solved the 3D Fick diffusion equation and the present 3D full coupled model.The main advantage of this new coupled formulation is the introduction of both thickness and material layer effects using a simpler and more elegant mathematical formulation having a faster convergence: in fact, a reduced number of fictitious layers F is requested in comparison with the uncoupled 3D model.

Figure 1 :
Figure 1: Geometries in the mixed curvilinear orthogonal reference system.

Figure 2 :
Figure 2: Benchmark 1 (B1), bulk modulus, and volume fraction of the ceramic phase for a square plate embedding a single FGM layer ( = p 1.0).

Figure 5 :
Figure 5: Benchmark 2 (B2), bulk modulus and volume fraction of the ceramic phase for a circular cylinder embedding a single FGM layer ( = p 2.0).

Figure 8 :
Figure 8: Benchmarks 3 and 4 (B3-B4), bulk modulus and volume fraction of the ceramic phase for a sandwich cylindrical shell panel and a sandwich spherical shell panel embedding a FGM core ( = p 0.5).

Table 1 :
Mechanical and hygroscopic properties for the two PVR and for benchmarks (B)

Table 2 :
Moisture contents, geometrical data, and lamination configurations for the two PVR

Table 3 :
[32]iminary validation results (PVR1), square plate embedding a single FGM layer.Moisture content applied at the external surfaces.The reference solutions are the 3D uncoupled analytical model and the 3D FE model via Patran and Nastran proposed by Brischetto and Torre[32]

Table 4 :
[32]iminary validation results (PVR2), cylindrical shell embedding a single FGM layer.Moisture content applied at the external surfaces.The reference solutions are the 3D uncoupled analytical model and the 3D FE model via Patran and Nastran proposed by Brischetto and Torre[32]Coupled hygro-elastic 3D model for steady-state analysis of FG plates and shells  11

content profile (a/h = 50)
Coupled hygro-elastic 3D model for steady-state analysis of FG plates and shells  17 in the second column of Table Different thickness ratios are considered to better understand the effects of geometrical parameters.The thickness ratio varies from ∕ = R h 2