*Article* **3D Stress Analysis of Multilayered Functionally Graded Plates and Shells under Moisture Conditions**

**Salvatore Brischetto \* and Roberto Torre**

Department of Mechanical and Aerospace Engineering, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy; roberto.torre@polito.it

**\*** Correspondence: salvatore.brischetto@polito.it; Tel.: +39-011-090-6813; Fax: +39-011-090-6899

**Abstract:** This paper presents the steady-state stress analysis of single-layered and multilayered plates and shells embedding Functionally Graded Material (FGM) layers under moisture conditions. This solution relies on an exact layer-wise approach; the formulation is unique despite the geometry. It studies spherical and cylindrical shells, cylinders, and plates in an orthogonal mixed curvilinear coordinate system (*α*, *β*, *z*). The moisture conditions are defined at the external surfaces and evaluated in the thickness direction under steady-state conditions following three procedures. This solution handles the 3D Fick diffusion equation, the 1D Fick diffusion equation, and the a priori assumed linear profile. The paper discusses their assumptions and the different results they deliver. Once defined, the moisture content acts as an external load; this leads to a system of three non-homogeneous secondorder differential equilibrium equations. The 3D problem is reduced to a system of partial differential equations in the thickness coordinate, solved via the exponential matrix method. It returns the displacements and their z-derivatives as a direct result. The paper validates the model by comparing the results with 3D analytical models proposed in the literature and numerical models. Then, new results are presented for one-layered and multilayered FGM plates, cylinders, and cylindrical and spherical shells, considering different moisture contents, thickness ratios, and material laws.

**Keywords:** functionally graded materials; 3D shell model; steady-state hygro-elastic analysis; Fick moisture diffusion equation; moisture content profile; layer-wise approach

### **1. Introduction**

The environmental conditions characterizing the service life of structural components can be adverse in many applications. The aerospace field gives several examples of changing environments, which results in temperature gradients and moisture concentration variability. It is crucial to consider such two factors and to include them in a proper structural analysis: the thermal and hygrometric fields induce an internal stress distribution, which changes as soon as the environmental conditions change. As with any stress field, it might induce failures on the structure [1], either on its own or because it sums up to that caused by classical mechanical loads. Composite, multilayered, and FGM-embedding structures require a special focus. Composite materials can also be degraded by moisture absorption and the following diffusion through the matrix [2]; multilayered structures highlight a strong heterogeneity in the hygro/thermal/mechanical properties; FGMs induce a further complication due to non-constant terms in governing equations. However, their boosted implementation in critical structural applications recently increased the attention of the researchers on these effects.

Laminated structures and composite materials suffer from a clear variation of the properties at the interfaces, which is the critical source of the delamination process [3]. Eliminating this discontinuity and substituting it with a smooth trend is the key achievement of Functionally Graded Materials (FGMs). They are advanced composite materials made by two or more different phases mixed with a continuous graded distribution. As a result, they are heterogeneous materials, delivering optimized responses for each application: the

**Citation:** Brischetto, S.; Torre, R. 3D Stress Analysis of Multilayered Functionally Graded Plates and Shells under Moisture Conditions. *Appl. Sci.* **2022**, *12*, 512. https:// doi.org/10.3390/app12010512

Academic Editor: Victor Franco Correia

Received: 26 November 2021 Accepted: 25 December 2021 Published: 5 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

advancements in processing technologies made it possible to control a unidirectional and even multidirectional variation. Combinations of a metallic and a ceramic phase are classic examples of FGMs finding application in severe thermal environments: they overcome the differences in thermal properties of the two constituents and, at the same time, deliver a reduced thermal stress distribution. Stiffness coefficients, hardness, thermal conductivity, moisture diffusivity, and corrosion resistance are just some examples of the performance characteristics that can be combined and enhanced [4].

Several recent and relevant works focused on structures embedding FG layers underline their importance in practical applications. From an analytical perspective, Sobhani et al. recently proposed several analytical results in the frame of the vibration analysis of FG shells. First, the governing equations followed the First-order Shear Deformation Theory (FSDT); then, they used the Generalized Differential Quadrature Method (GDQM), a semi-analytical solution method, to solve the system of partial differential equations. They studied conical shells embedding hybrid matrix/fiber nanocomposites in [5]; the approach was then extended to paraboloidal and hyperboloidal shells embedding polymer matrix, carbon fiber, and graphene nanoplatelets in [6]. Using the same approach, they studied sandwich conical-cylindrical-conical shells [7]. The layers are reinforced with functionally graded carbon nanotubes and graphene nanoplatelets; they described the elastic coefficients following an Equivalent Single Layer approach. They also studied five different patterns of CNT fibers distribution inside of the matrix while defining the vibrational behavior of coupled conical-conical shells [8]. They used the five-parameter shell theory and solved the differential equations through GDQM. Three-phase nanocomposites were also studied in hemispherical-cylindrical shells exploiting a similar approach [9] and defining the governing equations through the first-order shear deformation theory. From a numerical perspective, Rezaiee-Pajand et al. studied sandwich beams embedding FGM in two ways: they applied the Ritz method and the principle of minimum total potential energy within the framework of Timoshenko and Reddy beam theories in [10] to assess the bending of beams with different cross sections; they also developed a four-node isoparametric beam element to study porous beams with FGM [11]. Concerning two-dimensional geometries, the same authors proposed the nonlinear analysis of FG shells in [12,13] by improving an isoparametric six-node TRIA element with strain interpolation functions; the mechanical properties grading followed a power low. They also proposed a three-node TRIA element [14] using a mixed strain finite element approach and demonstrated that it is possible to get the exact response of the beams with a low number of elements under large deformations.

The solutions available in the literature seem to be promising. However, the results of a hygrometric stress analysis depend not only on the capabilities of the elastic model implemented but also on how the moisture content field has been evaluated, given the external boundary conditions. The hygrometric field acts as a field load; its quantification necessarily influences the stress analysis results. Aerospace applications usually involve thin components. Consequently, evaluating the moisture content field translates into determining its profile along with the thickness direction, generally coinciding with the grading direction of the mechanical/hygrometric properties. Developing a mechanical/elastic model for a thin component can be done at different levels of detail: 3D or 2D approaches, coupled with an analytical or numerical method. However, this might not be sufficient in defining the refinement of a model when hygrometric stresses are concerned. Even an analytical, exact 3D model would give wrong results when fed from an inaccurate moisture content field. The molecular diffusion depends on the gradient of the concentration and is described by Fick's law. A solid mathematical analogy between Fick's law and the Fourier heat conduction equation strongly simplifies the analysis and the dissertation. An exact solution comes from resolving three-dimensionally the diffusion/conduction equations; however, simplified solutions might benefit from their unidirectional simplification or an assumed-linear profile. Those three situations require an external evaluation of the moisture content field before the mechanical analysis. An alternative consists in defining a

coupled hygro-mechanical model, in which the moisture content field is a primary variable of the problem in analogy with the thermal field [15–20].

This paper discusses an hygro-elastic shell model, which relies on the exponential matrix method. This model does not limit to a specific geometry or lamination scheme; on the contrary, it handles plates, cylindrical, and spherical shells. Furthermore, it accepts one or more FGM layers on their own or coupled with homogeneous layers (sandwiches give an example). It extends the authors' 3D exact thermo-elastic shell model discussed in [21] to hygro-elastic stress analysis. The problem is defined under steady-state conditions only; the solution requires the moisture content amplitudes at the top and bottom external surfaces to be specified. As a first step, the solution evaluates the moisture content profile through the thickness direction. The authors considered three possible options: the 3D Fick diffusion law, its 1D simplified version, and an a priori assumed linear trend. Despite this, the model handles the solution via a layer-wise approach and through the exponential matrix method. The present authors are not the only ones using the exponential matrix method for solving the differential equilibrium equations. Soldatos and Ye [22] already used it in analyzing the free vibrations of cylinders; Messina [23] applied it to study multilayered plates. The orthogonal mixed curvilinear coordinate system helped study spherical shells in [24] using three transverse stress and three displacement components as primary variables of the problem. The analytical procedure is in analogy with the free vibration analysis and the static analysis under mechanical load discussed by the first author in [25–30]. Those previous formulations already handled different materials and geometries but lacked loads other than the mechanical ones. In those simpler cases, a set of three homogeneous differential equations are at the basis of the problem. However, the hygrometric load adds an additional term that makes them not homogeneous. The exponential matrix method also handles this feature, as discussed in [21] for the thermoelastic stress analysis of shells with FGMs. The closed-form solution of the problem is possible given the harmonic forms imposed on the variables (displacements and hygrometric field) and the simply supported boundary conditions. Moreover, the whole formulation benefits an orthogonal mixed curvilinear coordinate system, following the suggestions in [31–34]. This strategy introduces a set of curvature terms, the elements through which the equations automatically adapt to the different geometries. Such terms are a function of the thickness coordinate, together with the elastic and hygrometric coefficients in FGM layers. Introducing a set of fictitious (mathematical) layers allows obtaining constant coefficients and using the method discussed in [35] to solve the problem.

The literature overview offers different analytical and numerical 2D models handling the hygrometric stress analysis; they specifically focus on multilayered structures. Far fewer discuss the problem of structures embedding FGMs. The literature overview will demonstrate that, to the authors' knowledge, no analytical 3D shell model exists, in which the structures are different, provided that they have constant radii of curvature, and the moisture content evaluation follows three different approaches. Laoufi [36] studied rectangular plates embedding FGMs when subjected to different boundary conditions, including moisture content and temperature field. They developed an analytical method through the hyperbolic shear deformation plate model, which satisfied the stress boundary conditions and required no shear corrections. The volume fractions of the ceramic-metal constituents were used to define the materials grading in the thickness direction following a power law. The same power law was included in Inala's work [37], devoted to studying how the hygrothermal environment affects the vibration characteristics of plates embedding FGMs. To do this, the author developed a finite element model of the structures under investigation. Dai [38] and his colleagues focused on circular plates with a variable thickness along with the radial direction. They derived nonlinear governing equations for temperature, moisture content, and displacement fields and solved them through the differential quadrature method. Zenkour also studied how the thermal and hygrometric fields affect the bending and the buckling of plates embedding FGMs [39]. In this case, also, the authors considered a material grading in the thickness direction of the structure. Their formulation relies on an

exponential shear deformation theory and applies to plates resting on elastic foundations. Hamilton's principle was used to derive the equilibrium equations, and Navier's method to obtain the results. Boukelf also studied plates resting on elastic foundations, embedding FGMs [40]. The author developed a novel higher-order shear deformation theory model, deducing the problem equations through the virtual work principle. Hygro, thermal, and mechanical loads are all properties that can be handled. Analogously, Zidi [41] developed a further model for the same issue. In this case, the author used a four-variable refined plate theory. Analogous research has been proposed in [42] concerning functionally graded beams. The author deepened the influence of moisture content field and temperature on the bost-buckling response of beams. A nonlinear finite element solution was considered, in which the authors handled the kinematics of the bost-buckling through the Lagrangian approach.

More substantial is the literature discussing the moisture content effects on isotropic, orthotropic, and laminated structures. Chiba and Sugano [43] studied multilayered plates and proposed a two-dimensional analytical solution for hygrothermal effects on multilayered plates. The Classical Lamination Plate Theory and an analytical 3D plate model were compared in [44,45] while performing the hygro-thermal stress analysis. The CLT also received the attention of Kalil [46] when he investigated composite plates. Kollar [47] made an analytical investigation of composite cylindrical shells, while Shen [48] focused on buckling and post-buckling behaviors. The CLT was found to be inadequate in hygrothermal mechanical analysis by Lee and his colleagues [49]. There is no shortage of numerical models in this field; an example is given in [50], where plates with a central hole were studied. Multilayered structures under hygrometric fields have been the target of the finite element model by Khoshbakht et al. [51]; Kundu and Han studied buckling and vibration in multilayered shells with double curvature [52] through a finite element model grasping the hygrothermal effects. They extended this research to dynamic instability of doubly-curved shells embedding composite materials through a nonlinear finite element under orthogonal curvilinear coordinates [53]. Marques and Creus included the Fick diffusion law [54] into a FE shell model devoted to isotropic and multilayered structures under a hygrothermal environment. With the same aim, Naidu and Sinha [55] investigated cylindrical shell panels under large deflections through a higher-order shear deformation theory. Patel also proposed a higher-order FE for laminated parts [56]. Sereir et al. [57–59] considered the elastic properties variation in composite plates with the temperature and the moisture content and discussed a transient hygroscopic stress analysis. In [60–63], simplified solutions for hygro-thermal stresses analysis of composite plates considering the variation in elastic properties due to the hygrothermal environment are also considered. An analogous evaluation devoted to dynamic behavior has been proposed in [64,65]. Ghosh [66] used a FE model to investigate how a severe hygrothermal environment can affect the initiation and progress of damages in composites.

The literature survey highlighted that the hygromechanical stress analysis of structures embedding FGM layers still misses general and exact solutions as benchmarks in new solutions. Instead, the only analytical models work on defined and specific boundary conditions, laminations, and geometry. This paper intends to fill this gap by extending the authors' previous work on multilayered structures to FGMs. The manuscript is organized as follows. Section 3 describes the hygro-elastic shell problem and its solution with the exponential matrix method. Section 2 explores the moisture diffusion problem with the three approaches introduced before. The Fick law of diffusion has the same mathematical formulation of the Fourier heat conduction equation; the moisture diffusion problem and the heat conduction problem are in analogy, indeed, as demonstrated in [67], and further confirmed by Tay and Goh [68,69]. Therefore, Section 2 quantifies the moisture content profiles. Solving the 3D Fick diffusion law is possible by exploiting the analogy with the heat conduction problem. Tungikar and Rao [35] give a methodology that can also be applied in this context; the unidimensional Fick diffusion law disregards the diffusion fluxes in directions other than the thickness coordinate and can be calculated in analogy. Section 4 gives two sets of results. The first set, labeled as Assessments, is introduced to validate the problem. In the absence of 3D exact solutions for hygro-elastic problems in shells with FGM layers, the section exploits existing results for thermal stress analysis. After validating the present model and an additional FE model, the section uses this last FE model to verify the results when the hygrometric load replaces the thermal one. The second set, labeled as Benchmarks, discusses new results, which introduce further comments on the effects of moisture content profiles, thickness ratio, material, and lamination scheme, together with the impact of the geometry. The main conclusions and the further development are then summarized in Section 5.

#### **2. Fick Moisture Diffusion Equation**

This model relies on a decoupled solution, which means that the moisture content is evaluated separately and enters the elastic part of the problem as a known term. The essential hypothesis to obtain exact closed-form solutions is that all the problem variables have a harmonic form. For what concerns the moisture content, this means

$$\mathcal{M}^k(\mathfrak{a}, \mathfrak{z}, z) = M^k(z) \sin(\hbar a) \sin(\mathfrak{z} \mathfrak{z}) \text{ , \tag{1}$$

*Mk*(*z*) indicates the moisture content amplitude. Referring to Figures 1 and 2: *m* and *n* are the half-wave numbers in the two in-plane directions *α* and *β*, respectively. *a* and *b* are the shell dimensions referred to the mid-surface Ω0; they allow calculating the terms *α*¯ = *<sup>m</sup><sup>π</sup> a* and *β*¯ = *<sup>n</sup><sup>π</sup> <sup>b</sup>* . Note that the harmonic form of the moisture content reduces its assessment to its profile along with the thickness direction *z*. What happens in *α* and *β* directions is already defined through the sinusoidal functions.

**Figure 1.** Geometrical data and coordinate system for plates and cylinders. The figure also shows the stacking sequence used in Benchmarks 1 and 2.

The hygrometric boundary conditions are the moisture content amplitudes at the bottom and the top of the shell, *Mb*, and *Mt*, respectively. Three approaches exist to evaluate the moisture content profile; this solution implements all three, allowing a homogeneous comparison among them and with other methods. The solution is accomplished in analogy with what the authors achieved in [21] concerning the temperature field. Indeed, the Fourier heat conduction equation, regulating the thermal phenomenon, features the same mathematical expression of the Fick diffusion law, which applies to the moisture diffusion problem. The three approaches are here listed:


**Figure 2.** Geometrical data and coordinate system for cylindrical and spherical shell panels. The figure also shows the stacking sequence used in Benchmarks 3 and 4.

All the acronyms report the term 3D: it states that the elastic part of the shell model relies on a three-dimensional solution. Instead, the part of the acronym inside the parentheses defines how the moisture content profile has been quantified. "c" means calculated, either via the 3D or the 1D Fick diffusion laws; "a" means linear assumed.

#### *2.1. 3D Fick Equation*

By analogy with the heat flux *q*, suppose to define a moisture flux *g*. Then, the differential equation of moisture diffusion into a homogeneous solid, in the absence of chemical reactions and under steady-state conditions ( *<sup>∂</sup>*<sup>M</sup> *<sup>∂</sup><sup>t</sup>* = 0) reads:

$$\nabla \cdot \mathbf{g}(\mu\_1, \mu\_2, \mu\_3) = 0 \,. \tag{2}$$

Equation (2) is written in a general orthogonal curvilinear coordinate system (*u*1, *u*2, *u*3); however, in these conditions, the divergence of the moisture flux takes the following expression:

$$\nabla \cdot \mathbf{g} = \frac{1}{a} \left[ \frac{\partial}{\partial u\_1} \left( \frac{a}{a\_1} \mathbf{g}\_1 \right) + \frac{\partial}{\partial u\_2} \left( \frac{a}{a\_2} \mathbf{g}\_2 \right) + \frac{\partial}{\partial u\_3} \left( \frac{a}{a\_3} \mathbf{g}\_3 \right) \right], \tag{3}$$

*g*1, *g*2, *g*<sup>3</sup> express the components of the 3D flux in directions 1, 2, and 3; their explicit form is

$$g\_i = -\mathcal{D}\_i \frac{1}{a\_i} \frac{\partial \mathcal{M}}{\partial u\_i} \, , \tag{4}$$

D*<sup>i</sup>* is the diffusion coefficient along with direction *i*; *a*1, *a*2, and *a*<sup>3</sup> are the so-called scale factors, with *a* the product of the three factors (*a* = *a*<sup>1</sup> *a*<sup>2</sup> *a*<sup>3</sup> ). The problem takes place in a mixed curvilinear orthogonal coordinate system (*α*, *β*, *z*); Povstenko [32] discussed that, in this context, Equation (3) might be rewritten as follows:

$$\frac{1}{H\_{\text{a}}H\_{\text{f}}} \left[ \frac{\partial}{\partial \mathfrak{a}} \left( \frac{H\_{\text{a}}H\_{\text{f}}}{H\_{\text{a}}} \mathcal{D}\_{1} \frac{1}{H\_{\text{a}}} \frac{\partial \mathcal{M}}{\partial \mathfrak{a}} \right) + \frac{\partial}{\partial \mathfrak{f}} \left( \frac{H\_{\text{a}}H\_{\text{f}}}{H\_{\text{f}}} \mathcal{D}\_{2} \frac{1}{H\_{\text{f}}} \frac{\partial \mathcal{M}}{\partial \mathfrak{f}} \right) \right] + \mathcal{D}\_{3} \frac{\partial^{2} \mathcal{M}}{\partial \mathfrak{z}^{2}} = 0 \,,\tag{5}$$

FGM layers distinguish from classical isotropic and orthotropic layers as the diffusion coefficients D1(*z*), D2(*z*), and D3(*z*) are a function of the thickness coordinate *z*. Despite the lamination scheme, *Hα*(*z*) and *Hβ*(*z*) are two parametric coefficients, defined as follows:

$$H\_{\mathfrak{a}}(z) = \left(1 + \frac{z}{R\_{\mathfrak{a}}}\right) = H\_{\mathfrak{a}}(\tilde{z}) = \left(1 + \frac{\tilde{z} - h/2}{R\_{\mathfrak{a}}}\right),\tag{6}$$

$$H\_{\mathbb{R}}(z) = \left(1 + \frac{z}{R\_{\beta}}\right) = H\_{\mathbb{R}}(\overline{z}) = \left(1 + \frac{\overline{z} - h/2}{R\_{\beta}}\right). \tag{7}$$

When considering shells with constant radii of curvature, *R<sup>α</sup>* and *Rβ*, they are a linear function of the thickness coordinate *z*, which varies from −*h*/2 to +*h*/2 or *z*˜, which varies from 0 to *h*. *h* is the global thickness. The thickness coordinate *z* (or *z*˜) is rectilinear; for consistency, a further coefficient might be defined: *Hz* = 1. Equation (5) points out two main blocks; the second block has a simpler formulation as the third coordinate *z* is rectilinear. Actually, this point is further confirmed by the work of Leissa [70], which showed that

$$a\_1 = H\_{\mathfrak{a}}, \ a\_2 = H\_{\mathfrak{f}}, \ a\_3 = H\_{\mathfrak{z}} = 1. \tag{8}$$

The coefficients of Equation (5) are not constant even inside a *k*-th physical layer. This issue is due to the parametric coefficients *H<sup>α</sup>* and *Hβ*, a function of the thickness coordinate *z*, and to D1, D2, and D<sup>3</sup> which are not constant in FGM layers. A possible solution consists in dividing each physical layer into sufficiently thin mathematical (fictitious) layers. There, calculating the parametric and the diffusion coefficients at its middle allows moving the differential operators to the moisture content only. Inside each mathematical layers it holds that

$$\mathcal{D}\_1^{\*j} \frac{\partial^2 \mathcal{M}}{\partial a^2} + \mathcal{D}\_2^{\*j} \frac{\partial^2 \mathcal{M}}{\partial \beta^2} + \mathcal{D}\_3^{\*j} \frac{\partial^2 \mathcal{M}}{\partial z^2} = 0 \,. \tag{9}$$

where

$$\mathcal{D}\_1^{\*\dot{j}} = \frac{\mathcal{D}\_1^{\prime}}{H\_a^{2\dot{j}}}, \quad \mathcal{D}\_2^{\*\dot{j}} = \frac{\mathcal{D}\_2^{\dot{j}}}{H\_\beta^{2\dot{j}}}, \quad \mathcal{D}\_3^{\*\dot{j}} = \mathcal{D}\_3^{\dot{j}}.\tag{10}$$

Equation (9) is automatically satisfied by the harmonic expression of the moisture content M(*α*, *β*, *z*) already discussed in Equation (1). Still, the form of the moisture content amplitude in the thickness direction *M*(*z*) needs to be clarified. A tentative function is

$$M^j(z) = S\_1^j \cosh(s\_1^j z) + S\_2^j \sinh(s\_1^j z) \tag{11}$$

Each mathematical layer feature a coefficients *s<sup>j</sup>* and a pair of coefficients *M<sup>j</sup>* <sup>0</sup>. The first can be computed, introducing the harmonic form of moisture content, provided with the assumption of Equation (11) on its amplitude, in Equation (9):

$$s\_{1,2}^j = \pm \sqrt{\frac{\mathcal{D}\_1^{\*j}\bar{\mathbf{a}}^2 + \mathcal{D}\_2^{\*j}\bar{\mathbf{b}}^2}{\mathcal{D}\_3^{\*j}}},\tag{12}$$

and choosing the positive solution, referred to as *s j* <sup>1</sup>. Equation (11) shows that a pair of coefficients per mathematical layer *j* is needed, in addition to *s j* <sup>1</sup>. Each layer features its own set of coefficients, which means that still 2 × *G* unknowns are involved. However, at each interface between the fictitious layers, the following continuity conditions must hold:

$$\mathcal{M}\_b^{(j+1)} = \mathcal{M}\_t^j \,. \tag{13}$$

$$\mathcal{D}\_3^{\*j+1} M^{(j+1)}\_{,z\_b} = \mathcal{D}\_3^{\*j} M^j\_{,z\_t} \,. \tag{14}$$

The physical meaning of both the equations is obvious; Equation (13) states that the moisture content at the bottom of any (*j* + 1)-th layer must equal that at the top of the *j*-th layer. Likewise, Equation (14) states that the moisture flux component in the thickness direction at the bottom of any (*j* + 1)-th layer must equal that at the top of the *j*-th layer. The matrix form of Equations (13) and (14) allows compacting the analysis:

$$
\begin{bmatrix} S\_1 \\ S\_2 \end{bmatrix}^{j+1} = \begin{bmatrix} V\_{M\_1}^{j+1,j} & V\_{M\_2}^{j+1,j} \\ V\_{M\_3}^{j+1,j} & V\_{M\_4}^{j+1,j} \end{bmatrix} \begin{bmatrix} S\_1 \\ S\_2 \end{bmatrix}^j. \tag{15}
$$

Both the continuity conditions hold at the interfaces between layers; consequently, 2 × (*G* − 1) conditions can be imposed. The recursive use of Equation (15) allows linking the coefficients of the bottom layer (*j* = 1) with those of the top (*j* = *G*); such an idea can be compacted by identifying the transfer matrices of Equation (15) with the name *V*(*j*+1,*j*) *<sup>M</sup>* .

$$
\begin{bmatrix} \mathbf{S}\_1 \\ \mathbf{S}\_2 \end{bmatrix}^G = \mathbf{V}\_M^{(G,G-1)} \mathbf{V}\_M^{(G-1,G-2)} \dots \dots \mathbf{V}\_M^{(3,2)} \mathbf{V}\_M^{(2,1)} \begin{bmatrix} \mathbf{S}\_1 \\ \mathbf{S}\_2 \end{bmatrix}^1 = \mathbf{V}\_M^{(G,1)} \begin{bmatrix} \mathbf{S}\_1 \\ \mathbf{S}\_2 \end{bmatrix}^1. \tag{16}
$$

Two conditions are still missing to quantify all the 2 × *G* coefficients. However, simply imposing the moisture content at the top and the bottom of the structure gives the missing information. The coefficients for the external layers derive from this input and Equation (16). Once they are known, Equation (15) allows calculating all the remaining ones. Then, the moisture content profile is determined along with the thickness direction. The 3D hygro-elastic model including this hygrometric profile is referred to as 3D(M*c*,3D).

#### *2.2. 1D Fick Equation*

The three-dimensional problem of defining the moisture content field can be simplified when the structure under investigation is thin enough; such condition is expressed through sufficiently high thickness ratios. By recalling Equation (4) and the harmonic form for the moisture content field, the three components of the moisture content flux inside a *k*-th layer are

$$\mathcal{R}\_1^k = \mathcal{D}\_1^{\*k}(z)\bar{a}\mathcal{M}^k(z)\cos(\bar{a}a)\sin(\bar{\beta}\beta) \text{ , \tag{17}$$

$$\log\_2^k = \mathcal{D}\_2^{\*k}(z)\bar{\beta}\mathcal{M}^k(z)\sin(\hbar\alpha)\cos(\bar{\beta}\beta) \text{ ,} \tag{18}$$

$$\log\_3^k = \mathcal{D}\_3^{\*k}(z) M\_{;z}^k(z) \sin(\bar{a}a) \sin(\bar{\beta}\beta) \,. \tag{19}$$

As already seen, the diffusion coefficients are, in general, a function of the thickness coordinate as they are not constant inside FGM layers. The relative weight of the moisture content fluxes *g*<sup>1</sup> and *g*<sup>2</sup> decreases if compared to *g*<sup>3</sup> the thinner is the shell. When this happens, the first two components can be disregarded; Equation (19) now becomes

$$\text{gg}\_3(z) = \left(\mathcal{D}\_3^\* \frac{\partial M}{\partial z}\right) = \text{const.}\tag{20}$$

As already seen in the previous section, there is no practical difference between D<sup>∗</sup> 3 and D<sup>3</sup> because the third coordinate 3 ≡ *z* is rectilinear, and its parametric coefficient *Hz* equals 1. Equation (20) implies that *g*3(*z*) is actually constant along with *z*; it can be specialized for a generic *j*-th mathematical layer. The coefficient D<sup>∗</sup> <sup>3</sup> does not change inside it; this implies that the derivative of the moisture content to *z* is constant: the moisture content is linear inside each mathematical layer. Consequently Equation (20) can be simplified and rewritten in an algebraic form:

$$\mathcal{G}^{j}\_{3} = -\mathcal{D}^{\*j}\_{3} \frac{\partial \mathcal{M}^{j}}{\partial z} = -\frac{\mathcal{D}^{\*j}\_{3}}{h^{\flat}} (\mathcal{M}^{j}\_{t} - \mathcal{M}^{j}\_{b}) \,. \tag{21}$$

For the term D<sup>∗</sup> 3 *j* /*h<sup>j</sup>* , it is possible to exploit the analogy with the electrical conductance, as already done in [21] where a sort of thermal conductance was defined: the layer diffuses more, becoming thinner, increasing D<sup>∗</sup> <sup>3</sup> . This analogy is beneficial, as it enables defining an equivalent moisture content diffusion resistance for the overall structure:

$$\mathcal{R}\_{z\_{eq}} = \sum\_{j=1}^{G} \frac{h^j}{\mathcal{D}\_3^{\*j}} \,. \tag{22}$$

The moisture flux in the thickness direction can now be easily quantified, as the dissertation demonstrated it is constant and quantified the equivalent diffusion resistance. Defined *Mt* the moisture content amplitude on the top external surface and *Mb* that on the bottom, *g*<sup>3</sup> is

$$\mathbf{g}\_3 = -\frac{1}{R\_{z\_{eq}}} (M\_\mathbf{f} - M\_\mathbf{b}) = \text{const.}\tag{23}$$

As *g*<sup>3</sup> is constant across all the layers, the moisture content at any interface (and even at any *z* coordinate) is easily obtained:

$$\mathbf{g}\_3^j = -\mathcal{D}\_3^{\*j} \frac{(M\_t^j - M\_b^j)}{h^j} = \mathbf{g}\_3^{j+1} = -\mathcal{D}\_3^{\*j+1} \frac{(M\_t^{j+1} - M\_b^{j+1})}{h^{j+1}} = \mathbf{g}\_3 = \text{const.} \,, \tag{24}$$

The coefficient D<sup>∗</sup> 3 *<sup>j</sup>* changes in the thickness direction when FGM layers are considered in the stacking sequence. As a consequence, the moisture flux can keep constant only if the slope of the moisture content profile modifies accordingly. In this perspective, the effect of the material is considered (which means the stacking sequence and the FGM law), but the impact of the thickness is disregarded. The 3D hygro-elastic model including this hygrometric profile is referred to as 3D(M*c*,1D).

#### *2.3. Assumed Linear Moisture Content through the Thickness Direction*

The analysis is further simplified if also the effect of the material is disregarded, in addition to that of the thickness. A common assumption in the literature is that the moisture content profile is linear throughout the thickness direction of the shell: it does not take into account any change in the hygroscopic properties of the layers, as well as the fluxes in *α* and *β* directions. It is close to reality only when the shell is really thin and embeds a single layer or even different, but with homogeneous hygroscopic properties. The profile is immediately determined once the top and bottom external sovra-temperatures are set. The 3D hygroelastic model including this assumed hygrometric profile is referred to as 3D(M*a*).

#### **3. 3D Exact Shell Model for Hygro-Elastic Stress Analysis**

The 3D equilibrium equations for shells are written in the orthogonal mixed curvilinear coordinate system (*α*, *β*, *z*) shown in Figures 1 and 2. These equations are modified using 3D constitutive equations for Functionally Graded Materials (FGMs) and general geometrical relations for shells in (*α*, *β*, *z*) coordinates. Therefore, the system includes three differential equations of second order in *z*, and the related coefficients are not constant because of the radii of curvature and elastic coefficients for FGMs. A reasonable number of mathematical layers is introduced to obtain constant-coefficient equations; redoubling the number of variables allows reducing the order of the differential equations. Simply-supported sides and harmonic variables allow the analytical calculation of the partial derivatives in *α* and *β* directions. The final system, including the moisture content profile, has only first-order partial derivatives in *z*, and the exponential matrix method allows determining both general and particular solutions.

The investigated multilayered shells and plates subjected to a moisture content M(*α*, *β*, *z*) at the external surfaces have *k* classical/composite layers and/or functionally graded material layers. Stains defined in an orthogonal mixed curvilinear reference

system (*α*,*β*,*z*) are the algebraic summation of mechanical elastic parts (subscript *m*) and hygroscopic parts (subscript M). The 6 × 1 vector (*<sup>k</sup> αα*, *<sup>k</sup> ββ*, *<sup>k</sup> zz*, *γ<sup>k</sup> <sup>β</sup>z*, *<sup>γ</sup><sup>k</sup> <sup>α</sup>z*, *γ<sup>k</sup> αβ*) for each *k* layer is defined as

$$
\varepsilon\_{aa}^{k} = \epsilon\_{aan}^{k} - \epsilon\_{aa,\mathcal{M}}^{k} = \frac{1}{H\_{\mathfrak{a}}(z)} \frac{\partial u^{k}}{\partial a} + \frac{w^{k}}{H\_{\mathfrak{a}}(z)R\_{\mathfrak{a}}} - \eta\_{a}^{k}(z)\mathcal{M}^{k},\tag{25}
$$

$$
\epsilon^{k}\_{\beta\beta} = \epsilon^{k}\_{\beta\beta w} - \epsilon^{k}\_{\beta\beta \mathcal{M}} = \frac{1}{H\_{\beta}(z)} \frac{\partial v^{k}}{\partial \beta} + \frac{w^{k}}{H\_{\beta}(z)R\_{\beta}} - \eta^{k}\_{\beta}(z)\mathcal{M}^{k} \,. \tag{26}
$$

$$
\epsilon\_{zz}^k = \epsilon\_{zzm}^k - \epsilon\_{zz,\mathcal{M}}^k = \frac{\partial w^k}{\partial z} - \eta\_z^k(z)\mathcal{M}^k,\tag{27}
$$

$$
\gamma^k\_{\beta z} = \gamma^k\_{\beta zw} = \frac{1}{H\_\beta(z)} \frac{\partial w^k}{\partial \beta} + \frac{\partial v^k}{\partial z} - \frac{v^k}{H\_\beta(z)R\_\beta} \,' \tag{28}
$$

$$
\gamma\_{az}^k = \gamma\_{azm}^k = \frac{1}{H\_\hbar(z)} \frac{\partial w^k}{\partial a} + \frac{\partial u^k}{\partial z} - \frac{u^k}{H\_\hbar(z)R\_\hbar} \tag{29}
$$

$$
\gamma\_{a\S}^k = \gamma\_{a\S m}^k = \frac{1}{H\_a(z)} \frac{\partial v^k}{\partial a} + \frac{1}{H\_\emptyset(z)} \frac{\partial u^k}{\partial \beta} \tag{30}
$$

In the hygro-elastic strains, the three displacement components *uk*, *vk*, and *w<sup>k</sup>* and the scalar moisture content M*<sup>k</sup>* are defined in the reference system (*α*, *<sup>β</sup>*, *<sup>z</sup>*). Moisture expansion coefficients *η<sup>k</sup> <sup>α</sup>*(*z*), *η<sup>k</sup> <sup>β</sup>*(*z*) and *<sup>η</sup><sup>k</sup> <sup>z</sup>* (*z*) in the *k* physical layer could depend on the *z* coordinate in the case of Functionally Graded Material (FGM) layers; they are defined in the structural reference system (*α*, *β*, *z*) starting from the moisture expansion coefficients *ηk* <sup>1</sup>(*z*), *<sup>η</sup><sup>k</sup>* <sup>2</sup>(*z*) and *<sup>η</sup><sup>k</sup>* <sup>3</sup>(*z*) in the material reference system (1, 2, 3). The partial derivatives are defined via the symbol *∂*.

As anticipated, an essential hypothesis for exact closed-form solutions is the harmonic forms for all the variables: displacement components and moisture content. While the latter has already been discussed, for the displacement components it holds that

$$u^k(a,\beta,z) = \mathcal{U}^k(z)\cos(\hbar a)\sin(\vec{\beta}\beta)\ ,\tag{31}$$

$$\upsilon^k(\alpha, \beta, z) = V^k(z) \sin(\hbar \alpha) \cos(\vec{\beta} \beta) \ , \tag{32}$$

$$w^k(\alpha, \beta, z) = \mathcal{W}^k(z) \sin(\mathbb{H}\alpha) \sin(\beta \beta) \text{ , \tag{33}$$

displacement amplitudes are indicated as (*Uk*(*z*), *Vk*(*z*), *Wk*(*z*)). As already described, *m* and *n* are the half-wave numbers, *a* and *b* the mid-surface dimensions, and *α*¯ = *<sup>m</sup><sup>π</sup> <sup>a</sup>* and *β*¯ = *<sup>n</sup><sup>π</sup> b* .

The three-dimensional equilibrium equations for spherical shells having constant radii of curvature *R<sup>α</sup>* = *R<sup>β</sup>* and a total number *NL* of physical *k* (either classical or FGM) layers are

$$H\_{\not\!\!\!} (z) \frac{\partial \sigma\_{aa}^{k}}{\partial a} + H\_{\not\!\!\!a} (z) \frac{\partial \sigma\_{a\beta}^{k}}{\partial \beta} + H\_{\not\!\!\!\!\/)} H\_{\not\!\!\!/} (z) \frac{\partial \sigma\_{az}^{k}}{\partial z} + (\frac{2H\_{\not\!\!\!\!} (z)}{R\_{a}} + \frac{H\_{\not\!\!\!\!\/)} \sigma\_{az}^{k} = 0,\tag{34}$$

$$H\_{\not\beta}(z)\frac{\partial \sigma^{k}\_{\text{a}\beta}}{\partial a} + H\_{\text{a}}(z)\frac{\partial \sigma^{k}\_{\beta\beta}}{\partial \beta} + H\_{\text{a}}(z)H\_{\beta}(z)\frac{\partial \sigma^{k}\_{\beta z}}{\partial z} + (\frac{2H\_{\text{a}}(z)}{R\_{\beta}} + \frac{H\_{\beta}(z)}{R\_{\text{a}}})\sigma^{k}\_{\beta z} = 0,\tag{35}$$

$$H\_{\not\beta}(z)\frac{\partial v\_{az}^{k}}{\partial a} + H\_{\not\mathfrak{k}}(z)\frac{\partial v\_{\not\beta z}^{k}}{\partial \beta} + H\_{\not\mathfrak{k}}(z)H\_{\not\beta}(z)\frac{\partial v\_{zz}^{k}}{\partial z} - \frac{H\_{\not\mathfrak{k}}(z)}{R\_{\not\mathfrak{k}}}v\_{aa}^{k} - \frac{H\_{\not\mathfrak{k}}(z)}{R\_{\not\mathfrak{k}}}v\_{\not\beta \ell}^{k} + (\frac{H\_{\not\mathfrak{k}}(z)}{R\_{\not\mathfrak{k}}} + \frac{H\_{\not\mathfrak{k}}(z)}{R\_{\not\mathfrak{k}}})\sigma\_{zz}^{k} = 0. \tag{36}$$

Three-dimensional equilibrium equations for cylinders/cylindrical panels and plates are obtained from Equations (34)–(36) when one of the two radii of curvature is infinite and when both the radii of curvature are infinite, respectively.

The constitutive equations are developed by considering the algebraic summation of mechanical elastic strains and hygroscopic strains:

$$
\sigma^k = \mathcal{C}^k(z)\,\epsilon^k = \mathcal{C}^k(z)(\epsilon^k\_{\mathfrak{m}} - \mathfrak{e}^k\_{\mathcal{M}})\,\,\prime \tag{37}
$$

*<sup>σ</sup><sup>k</sup>* is the stress vector having a 6 <sup>×</sup> 1 dimension, *Ck*(*z*)is the 6 <sup>×</sup> 6 elastic coefficient matrix and it could depend on the *z* coordinate in the case of a *k*th FGM layer. The strains are included in the constitutive equations by using the form shown in Equations (25)–(30). Closed-form solutions are possible when orthotropic angles are 0◦ or 90◦ to obtain *Ck* <sup>16</sup>(*z*) = *<sup>C</sup><sup>k</sup>* <sup>26</sup>(*z*) = *<sup>C</sup><sup>k</sup>* <sup>36</sup>(*z*) = *<sup>C</sup><sup>k</sup>* <sup>45</sup>(*z*) = 0 in the structural reference system (*α*, *β*, *z*). Therefore,

$$\mathbf{C}^{k}(z) = \begin{bmatrix} \mathbf{C}\_{11}^{k}(z) & \mathbf{C}\_{12}^{k}(z) & \mathbf{C}\_{13}^{k}(z) & 0 & 0 & 0\\ \mathbf{C}\_{12}^{k}(z) & \mathbf{C}\_{22}^{k}(z) & \mathbf{C}\_{23}^{k}(z) & 0 & 0 & 0\\ \mathbf{C}\_{13}^{k}(z) & \mathbf{C}\_{23}^{k}(z) & \mathbf{C}\_{33}^{k}(z) & 0 & 0 & 0\\ 0 & 0 & 0 & \mathbf{C}\_{44}^{k}(z) & 0 & 0\\ 0 & 0 & 0 & 0 & \mathbf{C}\_{55}^{k}(z) & 0\\ 0 & 0 & 0 & 0 & 0 & \mathbf{C}\_{66}^{k}(z) \end{bmatrix} . \tag{38}$$

The explicit form of the constitutive equations develops thanks to the introduction of the geometrical Equations (25)–(30) into the constitutive Equation (37):

$$\sigma^{k}\_{\text{aa}} = \frac{\mathbb{C}^{k}\_{11}(z)}{H\_{\text{a}}(z)} \boldsymbol{u}^{k}\_{\text{a}} + \frac{\mathbb{C}^{k}\_{11}(z)}{H\_{\text{a}}(z)R\_{\text{a}}} \boldsymbol{w}^{k} + \frac{\mathbb{C}^{k}\_{12}(z)}{H\_{\text{\beta}}(z)} \boldsymbol{v}^{k}\_{\text{\beta}} + \frac{\mathbb{C}^{k}\_{12}(z)}{H\_{\text{\beta}}(z)R\_{\text{\beta}}} \boldsymbol{w}^{k} + \mathbb{C}^{k}\_{13}(z)\boldsymbol{w}^{k}\_{z} - \mathbb{S}^{k}\_{\text{a}}(z)\boldsymbol{\mathcal{M}}^{k}, \tag{39}$$

$$\sigma^{k}\_{\beta\beta} = \frac{\mathbb{C}^{k}\_{12}(z)}{H\_{a}(z)}u^{k}\_{\mu} + \frac{\mathbb{C}^{k}\_{12}(z)}{H\_{a}(z)R\_{a}}w^{k} + \frac{\mathbb{C}^{k}\_{22}(z)}{H\_{\beta}(z)}v^{k}\_{\beta} + \frac{\mathbb{C}^{k}\_{22}(z)}{H\_{\beta}(z)R\_{\beta}}w^{k} + \mathbb{C}^{k}\_{23}(z)w^{k}\_{,z} - \mathfrak{z}^{k}\_{\beta}(z)\mathcal{M}^{k} \,, \tag{40}$$

$$\sigma\_{zz}^{k} = \frac{\mathbb{C}\_{13}^{k}(z)}{H\_{a}(z)}u\_{\mathcal{A}}^{k} + \frac{\mathbb{C}\_{13}^{k}(z)}{H\_{a}(z)R\_{a}}w^{k} + \frac{\mathbb{C}\_{23}^{k}(z)}{H\_{\beta}(z)}v\_{\mathcal{A}}^{k} + \frac{\mathbb{C}\_{23}^{k}(z)}{H\_{\beta}(z)R\_{\beta}}w^{k} + \mathbb{C}\_{33}^{k}(z)w\_{z}^{k} - \mathbb{f}\_{z}^{k}(z)\mathcal{M}^{k} \,, \tag{41}$$

$$
\sigma\_{\beta z}^{k} = \frac{\mathbb{C}\_{44}^{k}(z)}{H\_{\beta}(z)} w\_{\gamma \theta}^{k} + \mathbb{C}\_{44}^{k}(z) v\_{z}^{k} - \frac{\mathbb{C}\_{44}^{k}(z)}{H\_{\beta}(z) R\_{\beta}} v^{k},\tag{42}
$$

$$
\sigma\_{az}^{k} = \frac{\mathbb{C}\_{55}^{k}(z)}{H\_{a}(z)} \boldsymbol{w}\_{\boldsymbol{\mu}}^{k} + \mathbb{C}\_{55}^{k}(z) \boldsymbol{u}\_{\boldsymbol{z}}^{k} - \frac{\mathbb{C}\_{55}^{k}(z)}{H\_{a}(z)R\_{a}} \boldsymbol{u}^{k},\tag{43}
$$

$$
\sigma\_{a\beta}^{k} = \frac{\mathcal{C}\_{\theta\theta}^{k}(z)}{H\_{\theta}(z)} \upsilon\_{,a}^{k} + \frac{\mathcal{C}\_{\theta\theta}^{k}(z)}{H\_{\beta}(z)} \upsilon\_{,\beta}^{k} \tag{44}
$$

partial derivatives ( *<sup>∂</sup> ∂α* ), ( *<sup>∂</sup> ∂β* ), and ( *<sup>∂</sup> <sup>∂</sup><sup>z</sup>* ) are indicated in Equations (39)–(44) through subscripts (, *α*), (, *β*), and (, *z*), respectively. Terms *ξ<sup>k</sup> <sup>α</sup>*(*z*), *ξ<sup>k</sup> <sup>β</sup>*(*z*), and *<sup>ξ</sup><sup>k</sup> <sup>z</sup>* (*z*) designate the hygromechanical coupling coefficients in the structural reference system (*α*, *β*, *z*), and they could depend on *z* in the case of FGM layers:

$$\mathcal{G}\_{a}^{k}(z) = \mathbb{C}\_{11}^{k}(z)\eta\_{a}^{k}(z) + \mathbb{C}\_{12}^{k}(z)\eta\_{\beta}^{k}(z) + \mathbb{C}\_{13}^{k}(z)\eta\_{z}^{k}(z) \ , \tag{45}$$

$$\mathfrak{J}\_{\beta}^{k}(z) = \mathbb{C}\_{12}^{k}(z)\eta\_{a}^{k}(z) + \mathbb{C}\_{22}^{k}(z)\eta\_{\beta}^{k}(z) + \mathbb{C}\_{23}^{k}(z)\eta\_{z}^{k}(z) \, , \tag{46}$$

$$\mathcal{G}\_z^k(z) = \mathbb{C}\_{13}^k(z)\eta\_a^k(z) + \mathbb{C}\_{23}^k(z)\eta\_\beta^k(z) + \mathbb{C}\_{33}^k(z)\eta\_z^k(z) \,. \tag{47}$$

The final system is obtained thanks to the substitution of the harmonic form equations for displacements and moisture content Equations (1)–(31) and the modified constitutive relations Equations (39)–(44) into the three-dimensional equilibrium Equations (34)–(36) developed for spherical shells:

$$A\_1^\dagger \mathcal{U}^\dagger + A\_2^\dagger V^\dagger + A\_3^\dagger \mathcal{W}^\dagger + A\_4^\dagger \mathcal{U}\_{\neq}^\dagger + A\_5^\dagger \mathcal{W}\_{\neq}^\dagger + A\_6^\dagger \mathcal{U}\_{\neq}^\dagger + L\_1^\dagger \mathcal{M}^\dagger = \mathbf{0} \,, \tag{48}$$

$$A\_{\tau}^{\dagger} \mathcal{U}^{\dagger} + A\_{8}^{\dagger} V^{\dagger} + A\_{9}^{\dagger} \mathcal{W}^{\dagger} + A\_{10}^{\dagger} V\_{\neq}^{\dagger} + A\_{11}^{\dagger} \mathcal{W}\_{\neq}^{\dagger} + A\_{12}^{\dagger} V\_{\neq \ge}^{\dagger} + L\_{2}^{\dagger} \mathcal{M}^{\dagger} = 0 \,,\tag{49}$$

$$A\_{13}^{\dagger} \mathcal{U}^{\dagger} + A\_{14}^{\dagger} V^{\dagger} + A\_{15}^{\dagger} \mathcal{W}^{\dagger} + A\_{16}^{\dagger} \mathcal{U}\_{\overline{z}}^{\dagger} + A\_{17}^{\dagger} V\_{\overline{z}}^{\dagger} + A\_{18}^{\dagger} \mathcal{W}\_{\overline{z}}^{\dagger} + A\_{19}^{\dagger} \mathcal{W}\_{\overline{z}\overline{z}}^{\dagger} + L\_3^{\dagger} M\_{\overline{z}}^{\dagger} + L\_4^{\dagger} M^{\dagger} = 0. \tag{50}$$

The above equations have no constant coefficients because the parametric coefficients *H<sup>α</sup>* and *H<sup>β</sup>* depend on *z* in the case of shell geometries and/or elastic and hygrometric coefficients depend on *z* in FGM layers. For these reasons, each *k* physical layer is divided into a suitable number of mathematical layers indicated with a new index *j* which changes from 1 to the global number of mathematical layers *G*. In each *j* mathematical layer, the parametric coefficients *H<sup>α</sup>* and *H<sup>β</sup>* and the variable elastic and hygrometric coefficients for FGM layers can be exactly defined by using the *z* coordinate in the middle of each *j* layer. Therefore, coefficients *A<sup>j</sup> <sup>s</sup>* (*<sup>s</sup>* from 1 to 19) and *<sup>L</sup><sup>j</sup> <sup>r</sup>* (*r* from 1 to 4) become constant terms in the compact form of the system of differential equations in *z* defined in Equations (48)–(50).

Equations (48)–(50) indicate a system of three differential equations of second order in *z*. The unknowns are the displacement and moisture content amplitudes and the associated derivatives calculated to *z*. The derivatives in *α* and *β* have already been calculated via the harmonic forms for displacements and moisture content previously introduced. In this system, decoupling the variables is possible: this means a separate quantification of the moisture content profile through the thickness direction *z*, addressed in an appropriate section. Therefore, it becomes a known term.

The consequence of these choices is that the system contains second-order differential equations only, in the displacement amplitudes *U<sup>j</sup>* , *V<sup>j</sup>* , *W<sup>j</sup>* and their derivatives in *z*. Redoubling the variables as proposed in [71,72] allows reducing the system into a firstorder one in *<sup>z</sup>*. In each *<sup>j</sup>* layer, the 3 × 1 vector of unknowns (*U<sup>j</sup>* , *V<sup>j</sup>* , *W<sup>j</sup>* ) becomes a 6 × 1 unknown vector (*U<sup>j</sup>* , *V<sup>j</sup>* , *W<sup>j</sup>* , *U<sup>j</sup>* , *V<sup>j</sup>* , *W<sup>j</sup>* ); superscript indicates derivatives performed to *z* (also written as *<sup>∂</sup> <sup>∂</sup><sup>z</sup>* ). Terms *<sup>M</sup><sup>j</sup>* and *<sup>M</sup><sup>j</sup>* can be considered as known terms because they can be opportunely calculated, as will be demonstrated in the next section:

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *Aj* <sup>6</sup> 0 000 0 0 *A<sup>j</sup>* <sup>12</sup> 000 0 0 0 *A<sup>j</sup>* <sup>19</sup> 00 0 00 0 *A<sup>j</sup>* <sup>6</sup> 0 0 00 00 *A<sup>j</sup>* <sup>12</sup> 0 00 000 *A<sup>j</sup>* 19 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *Uj Vj W<sup>j</sup> Uj Vj W<sup>j</sup>* ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 000 *A<sup>j</sup>* <sup>6</sup> 0 0 0000 *A<sup>j</sup>* <sup>12</sup> 0 00000 *A<sup>j</sup>* 19 −*A<sup>j</sup>* <sup>1</sup> <sup>−</sup>*A<sup>j</sup>* <sup>2</sup> <sup>−</sup>*A<sup>j</sup>* <sup>3</sup> <sup>−</sup>*A<sup>j</sup>* <sup>4</sup> <sup>0</sup> <sup>−</sup>*A<sup>j</sup>* 5 −*A<sup>j</sup>* <sup>7</sup> <sup>−</sup>*A<sup>j</sup>* <sup>8</sup> <sup>−</sup>*A<sup>j</sup>* <sup>9</sup> <sup>0</sup> <sup>−</sup>*A<sup>j</sup>* <sup>10</sup> <sup>−</sup>*A<sup>j</sup>* 11 −*A<sup>j</sup>* <sup>13</sup> <sup>−</sup>*A<sup>j</sup>* <sup>14</sup> <sup>−</sup>*A<sup>j</sup>* <sup>15</sup> <sup>−</sup>*A<sup>j</sup>* <sup>16</sup> <sup>−</sup>*A<sup>j</sup>* <sup>17</sup> <sup>−</sup>*A<sup>j</sup>* 18 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *Uj Vj W<sup>j</sup> Uj Vj W<sup>j</sup>* ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ + ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0000 0 0 0000 0 0 0000 −*L<sup>j</sup>* <sup>1</sup> 0 0000 −*L<sup>j</sup>* <sup>2</sup> 0 0000 −*L<sup>j</sup>* <sup>4</sup> <sup>−</sup>*L<sup>j</sup>* <sup>3</sup> 0000 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *M<sup>j</sup> M<sup>j</sup>* 0 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (51)

By defining vectors *U<sup>j</sup>* = [*U<sup>j</sup> <sup>V</sup><sup>j</sup> <sup>W</sup><sup>j</sup> <sup>U</sup><sup>j</sup> Vj W<sup>j</sup>* ] *<sup>T</sup>*, *U<sup>j</sup>* = *<sup>∂</sup>U<sup>j</sup> <sup>∂</sup><sup>z</sup>* and *<sup>M</sup><sup>j</sup>* = [*M<sup>j</sup> M<sup>j</sup>* 0000] *<sup>T</sup>* (where *T* means the transpose of a vector), the following compact form is possible:

$$\mathbf{D}^{j}\mathbf{U}^{j'} = \mathbf{A}^{j}\mathbf{U}^{j} + \mathbf{L}^{j}\mathbf{M}^{j} \, , \tag{52}$$

the above equation, thanks to the definitions *A*∗*<sup>j</sup>* <sup>=</sup> *D<sup>j</sup>* −1 *A<sup>j</sup>* and *L*∗*<sup>j</sup>* <sup>=</sup> *D<sup>j</sup>* −1 *Lj* , can be rewritten as

$$\mathbf{U}^{\dot{j}'} = \mathbf{D}^{j-1} A^j \mathbf{U}^j + \mathbf{D}^{j-1} L^j \mathbf{M}^j,\tag{53}$$

$$\mathbf{U}^{\dot{j}} = \mathbf{A}^{\*\dot{j}} \mathbf{U}^{\dot{j}} + \mathbf{L}^{\*\dot{j}} \mathbf{M}^{\dot{j}}.\tag{54}$$

The moisture content profile through the thickness direction *z* can be calculated using one of the three methods proposed in the previous section. This profile can be reconstructed using a linear approximation of the moisture content in each *j* mathematical layer. This reconstruction can be defined as

$$M^j(\vec{z}^j) = a\_M^j \vec{z}^j + b\_{M'}^j \tag{55}$$

in a *j*th mathematical layer, *a j <sup>M</sup>* and *b j <sup>M</sup>* are two constant coefficients. The first represents the slope of the moisture content profile inside a mathematical layer; the second the moisture content at the bottom. *z*˜*<sup>j</sup>* is a local thickness coordinate for each *j* mathematical layer, and it changes from *0* at the bottom of the considered *j* mathematical layer to the top value *h<sup>j</sup>* of the same mathematical layer.

The system of first-order differential equations in *z*˜ or *z* shown in Equation (54) is not homogeneous because the hygroscopic term *L*∗*<sup>j</sup> M<sup>j</sup>* depends on *<sup>z</sup>*˜*<sup>j</sup>* or *<sup>z</sup><sup>j</sup>* . In the case of a generic system of non-homogeneous first-order differential equations having an unknown *<sup>G</sup>* <sup>×</sup> 1 vector *x*, a *<sup>G</sup>* <sup>×</sup> *<sup>G</sup>* <sup>|</sup>*A* matrix containing constant coefficients, and a known function vector *f*(*t*)=[ *<sup>f</sup>*1(*t*)... *fG*(*t*)]*T*,we can write

$$\frac{d\mathbf{x}}{dt} = A\mathbf{x} + f(t) \, \tag{56}$$

a possible solution of the Equation (56) can be obtained via the exponential matrix method:

$$\mathbf{x}(t) = \mathbf{e}^{\mathbf{A}(t-t\_0)}\mathbf{x}\_0 + \int\_{t\_0}^t \mathbf{e}^{\mathbf{A}(t-s)}f(s)ds\,. \tag{57}$$

The known term in Equation (54) can be given in the following complete form:

$$\mathbf{M}^{\circ^{j}} = \boldsymbol{L}^{\circ^{j}}\boldsymbol{M}^{\circ^{j}} = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ -\boldsymbol{L}\_{1}^{\circ^{j}} & 0 & 0 & 0 & 0 & 0 \\ -\boldsymbol{L}\_{2}^{\circ^{j}} & 0 & 0 & 0 & 0 & 0 \\ -\boldsymbol{L}\_{4}^{\circ^{j}} & -\boldsymbol{L}\_{3}^{\circ^{j}} & 0 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} \boldsymbol{a}\_{M}^{j}\boldsymbol{z}^{j} + \boldsymbol{b}\_{M}^{j} \\ \boldsymbol{a}\_{M}^{j} \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ -\boldsymbol{L}\_{1}^{\circ^{j}}(\boldsymbol{a}\_{M}^{j}\boldsymbol{z}^{j} + \boldsymbol{b}\_{M}^{j}) \\ -\boldsymbol{L}\_{2}^{\circ^{j}}(\boldsymbol{a}\_{M}^{j}\boldsymbol{z}^{j} + \boldsymbol{b}\_{M}^{j}) \\ -\boldsymbol{L}\_{2}^{\circ^{j}}(\boldsymbol{a}\_{M}^{j}\boldsymbol{z}^{j} + \boldsymbol{b}\_{M}^{j}) \\ -\boldsymbol{L}\_{4}^{\circ^{j}}(\boldsymbol{a}\_{M}^{j}\boldsymbol{z}^{j} + \boldsymbol{b}\_{M}^{j}) - \boldsymbol{L}\_{3}^{\circ^{j}}\boldsymbol{a}\_{M}^{j} \end{bmatrix}. \tag{58}$$

Therefore, Equation (54) can be rewritten as

$$\mathbf{U}^{\stackrel{\circ}{\prime}} = \mathbf{A}^{\*\stackrel{\circ}{\prime}} \mathbf{U}^{\stackrel{\circ}{\prime}} + \mathbf{M}^{\*\stackrel{\circ}{\prime}},\tag{59}$$

*M*∗*j* includes only linear and known functions in *z*˜*<sup>j</sup>* coordinate. The Equation (59) can be solved through the exponential matrix method:

$$\mathbf{U}^{j}(\tilde{z}^{j}) = \mathbf{e}^{(\mathbf{A}^{\*\tilde{j}}, \tilde{z}^{j})} \mathbf{U}^{j}(0) + \int\_{0}^{\mathcal{Z}^{j}} \mathbf{e}^{(\mathbf{A}^{\*\tilde{j}}(\tilde{z}^{j} - s))} \mathbf{M}^{\*\tilde{j}}(s) ds \,. \tag{60}$$

*A*∗∗*j* <sup>=</sup> <sup>e</sup>(*A*∗*<sup>j</sup> hj* ) and *L*∗∗*<sup>j</sup>* = *<sup>h</sup><sup>j</sup>* <sup>0</sup> <sup>e</sup>(*A*∗*<sup>j</sup>* (*hj* <sup>−</sup>*s*))*M*∗*<sup>j</sup>* (*s*)*ds* must be opportunely defined for each *j* layer having thickness *h<sup>j</sup>* . In this way, the displacement vector at the top of each *j* mathematical layer is calculated. Both terms are defined via the exponential matrix opportunely expanded and evaluated in each *j* mathematical layer with thickness *h<sup>j</sup>* :

$$A^{\ast \ast^j} = \mathbf{e}^{(A^{\ast^j}h^j)} = \mathbf{I} + A^{\ast^j}h^j + \frac{{\mathbf{A}^{\ast^j}}^2}{2!}h^{j^2} + \frac{{\mathbf{A}^{\ast^{j^3}}}^3}{3!}h^{j^3} + \dots + \frac{{\mathbf{A}^{\ast^{j^N}}}^N}{N!}h^{j^N},\tag{61}$$

$$\begin{split} \mathbf{L}^{\*,i^{j}} &= \int\_{0}^{\mathbf{i}^{j}} \mathbf{e}^{(\mathbf{A}^{i}(\mathbf{h}^{j}-\mathbf{s}))} \mathbf{M}^{\*,i}(\mathbf{s}) d\mathbf{s} = \int\_{0}^{\mathbf{i}^{j}} \left( \mathbf{I} + \mathbf{A}^{\*^{j}}(\mathbf{h}^{j}-\mathbf{s}) + \frac{\mathbf{A}^{\*^{j}}}{2!}(\mathbf{h}^{j}-\mathbf{s})^{2} + \frac{\mathbf{A}^{\*^{j}}}{3!}(\mathbf{h}^{j}-\mathbf{s})^{3} + \dots \right) \\ &\cdots + \frac{\mathbf{A}^{\*^{j}}}{N!}(\mathbf{h}^{j}-\mathbf{s})^{N} \right) \mathbf{M}^{\*^{j}}(\mathbf{s}) d\mathbf{s}, \end{split} \tag{62}$$

where *I* is the 6 <sup>×</sup> 6 identity matrix, and the integral given in Equation (60) can be calculated in each *j* layer having thickness *h<sup>j</sup>* by using the exponential matrix and the same order *N* already shown in Equation (61). By using Equations (61) and (62), Equation (60) is modified as

$$\mathbf{U}^{j}(\hbar^{j}) = \mathbf{A}^{\ast \ast^{j}} \mathbf{U}^{j}(0) + \mathbf{L}^{\ast \ast^{j}},\tag{63}$$

where *U<sup>j</sup> <sup>t</sup>* indicates *<sup>U</sup><sup>j</sup>* (*h<sup>j</sup>* ) and it is defined at the top *<sup>t</sup>* of each *<sup>j</sup>* layer, and *U<sup>j</sup> <sup>b</sup>* indicates *Uj* (0) and it is defined at the bottom of each *j* layer. In this way, Equation (63) is defined as

$$\mathbf{U}\_t^j = \mathbf{A}^{\ast \ast^j} \mathbf{U}\_b^j + \mathbf{L}^{\ast \ast^j}. \tag{64}$$

Using Equation (64) allows to connect displacements and their derivatives in *z* defined at the top of the *j* mathematical layer with the same variables defined at the bottom of the same *j* layer.

The general three-dimensional shell model is developed using a layer-wise approach. Inter-laminar continuity conditions in displacements and transverse stresses must be defined at each interface between the two adjacent mathematical layers. The inter-laminar continuity conditions for displacements come through congruence hypotheses:

$$u\_b^j = u\_t^{j-1}, \ v\_b^j = v\_t^{j-1}, \ w\_b^j = w\_t^{j-1} \,. \tag{65}$$

The inter-laminar continuity conditions for transverse shear and normal transverse stresses are defined employing equilibrium hypotheses:

$$
\sigma^j\_{zz\_b} = \sigma^{j-1}\_{zz\_t} \; , \; \sigma^j\_{az\_b} = \sigma^{j-1}\_{az\_t} \; , \; \sigma^j\_{\beta z\_b} = \sigma^{j-1}\_{\beta z\_t} \; . \tag{66}
$$

The formulation of this solution requires rewriting the inter-laminar continuity conditions, Equations (65) and (66), in a displacement form. Achieving this task is possible by recalling the constitutive Equations (39)–(44) and the harmonic form of the variables of the problem, Equations (1) and (31)–(33). This procedure allows writing the inter-laminar continuity condition in the amplitude displacements and their derivatives to the thickness coordinate. Then, compacting the notation is possible, recalling the vectors *U<sup>j</sup>* and *M<sup>j</sup>* , and introducing a pair of transfer matrices. The procedure is similar to that employed in [25]; exception made for an additional coefficient multiplying the moisture content amplitude:

$$
\begin{bmatrix} U \\ V \\ W \\ M' \\ V' \\ W' \end{bmatrix}\_b = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ T\_1 & 0 & T\_2 & T\_3 & 0 & 0 \\ 0 & T\_4 & T\_5 & 0 & T\_6 & 0 \\ T\_7 & T\_8 & T\_9 & 0 & 0 & T\_{10} \end{bmatrix}^{j-1,j} \begin{bmatrix} U \\ V \\ W \\ M' \\ V' \\ V' \\ W' \end{bmatrix}\_t + \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ T\_{11} & 0 & 0 & 0 & 0 & 0 \end{bmatrix}^{j-1,j} \begin{bmatrix} M \\ M \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}\_t \end{bmatrix} \tag{67}
$$

The first three rows of Equation (67) denote the displacement continuity equations; the last three are the continuity conditions for stresses. The compact form of Equation (67) follows:

$$\mathbf{U}\_{b}^{j} = T\_{\mathbf{U}}^{j-1,j} \mathbf{U}\_{t}^{j-1} + T\_{\mathbf{M}}^{j-1,j} \mathbf{M}\_{t}^{j-1}. \tag{68}$$

Equation (68) expresses the link between the displacements and their *z* derivatives, calculated at the bottom of the *j*th layer, with their corresponding values plus the moisture content (and its *z* derivative), at the top of the (*j* − 1)th layer.

The harmonic form implemented in all the variables of the problem automatically satisfies the simply supported boundary conditions, here reported for exhaustiveness:

$$w = \upsilon = 0, \ \sigma\_{aa} = 0 \quad \text{ for} \quad a = 0, a \,, \tag{69}$$

$$w = u = 0, \; \sigma\_{\mathbb{R}^6} = 0 \quad \text{for} \quad \beta = 0, b \tag{70}$$

This solution extends the model already seen in [25] for the static analysis of shells subjected to static loads. In that context, mechanical loads can act on the external top and bottom surfaces, with components defined in the three directions *α*, *β*, and *z*. They are grouped into vectors *<sup>P</sup>* = (*P<sup>α</sup> <sup>P</sup><sup>β</sup> Pz*)*T*; the superscript *<sup>G</sup>* means the one acting on the top surface; the one with superscript 1 identifies the one acting on the lower one. The effect of the external loads affects the displacements:

$$\mathbf{P}\_t^G \mathbf{U}\_t^G = \mathbf{P}\_t^G \, \!/ \tag{71}$$

$$B\_b^1 \mathcal{U}\_b^1 = P\_b^1. \tag{72}$$

more details are available in [25]. The hygro-elastic analysis does not involve anything different; once applied, the moisture content induces an equivalent load, which sums up the (possible) mechanical load. The matrices *B* convert the mechanical and/or hygrometric loads into displacements. *B<sup>G</sup> <sup>t</sup>* , in Equation (71) handles the top (t) of the last layer (G); *<sup>B</sup>*<sup>1</sup> *b*, in Equation (72), the bottom (b) of the first layer (1).

The algebraic system of Equations (71) and (72) can be reformulated into a matrix form, rewriting the displacements *U<sup>G</sup> <sup>t</sup>* <sup>=</sup> *<sup>U</sup>G*(*hG*) as a function of *<sup>U</sup>*<sup>1</sup> *<sup>b</sup>* <sup>=</sup> *<sup>U</sup>*1(0). This rewriting is possible through a recursive substitution of Equation (68) into Equation (64), linking the displacements at the top of the last layer (and their z derivatives) to those at the bottom of the first layer:

$$\begin{aligned} \mathcal{W}\_{i}^{G} &= \left( A^{\*,\*G} \mathcal{T}\_{\mathcal{U}}^{G-1,\mathcal{G}} A^{\*,\*G-1} \mathcal{T}\_{\mathcal{U}}^{G-2,\mathcal{G}-1} \dots A^{\*,\*2} \mathcal{T}\_{\mathcal{U}}^{1,2} A^{\*,\*} \right) \mathcal{U}\_{0}^{1} + \\ & \left( A^{\*,\*G} \mathcal{T}\_{\mathcal{U}}^{G-1,\mathcal{G}} A^{\*,\*G-1} \dots A^{\*,\*2} \mathcal{T}\_{\mathcal{U}}^{1,2} L^{\*,\*} + \\ & A^{\*,\*G} \mathcal{T}\_{\mathcal{U}}^{G-1,\mathcal{G}} A^{\*,\*G-1} \dots A^{\*,\*3} \mathcal{T}\_{\mathcal{U}}^{2,3} L^{\*,\*2} + \\ & \vdots \\ & \quad \vdots \\ & A^{\*,\*G} \mathcal{T}\_{\mathcal{U}}^{G-1,\mathcal{G}} L^{\*,\*G-1} + \\ & A^{\*,\*G} \mathcal{T}\_{\mathcal{U}}^{G-1,\mathcal{G}} A^{\*,\*G} \dots A^{\*,\*2} \mathcal{T}\_{\mathcal{U}}^{1,2} M\_{l}^{1} + \\ & A^{\*,\*G} \mathcal{T}\_{\mathcal{U}}^{G-1,\mathcal{G}} A^{\*,\*G} \dots A^{\*,\*3} \mathcal{T}\_{\mathcal{U}}^{1,3} M\_{l}^{1} + \\ & \vdots \\ & \quad \vdots \\ & A^{\*,\*G} \mathcal{T}\_{\mathcal{U}}^{G-1,\mathcal{G}} A^{\*,\*G-1} \mathcal{T}\_{\mathcal{M}}^{1,2} M\_{l}^{1} + \\ & \quad \vdots \\ & A^{\*,\*G} \mathcal{T}\_{\mathcal{M}}^{G-1,\*G} M\_{l}^{G-1} \end{aligned} \tag{7.2}$$

Equation (73) consists of two main blocks. The first has as a common multiplication factor the bottom displacements of the first layer *U*<sup>1</sup> *<sup>b</sup>*; the content in the brackets is a <sup>6</sup> <sup>×</sup> 6 matrix, which is identical to that defined *<sup>H</sup><sup>m</sup>* for classical elastic analysis in [25]. The hygrometric field brings with it additional terms included in the second block. The terms *<sup>L</sup>*∗∗*<sup>j</sup>* explicitly include the hygrometric profile; they are *<sup>G</sup>* as each mathematical layer features its own profile. The terms *M<sup>j</sup> <sup>t</sup>* set the moisture content at each interface; consequently, they are *G* − 1. This block is a known and constant term; it takes the form of a 6 <sup>×</sup> 1 vector, from now on referred to as *HM*. Therefore, the compact expression of Equation (73) takes the following form:

$$
\mathbf{U}\_t^G = \mathbf{H}\_m \mathbf{U}\_b^1 + \mathbf{H}\_M. \tag{74}
$$

Given this result, Equation (71) can be rewritten in terms of *U*<sup>1</sup> *b*:

$$B\_t^G H\_m \mathcal{U}\_b^1 = -B\_t^G H\_M \,. \tag{75}$$

Equations (72) and (75) now share the same unknown; they can be put to system as follows:

"

*E* = ! *BG <sup>t</sup> H<sup>m</sup>*

$$E\mathbf{U}\_b^1 = \mathbf{P}\_{M,b} \tag{76}$$

(77)

where

and

$$
\begin{array}{c}
\mathbf{L} = \begin{bmatrix} & \mathbf{B}\_b^1 & \\ \end{bmatrix} \\\\ \mathbf{P}\_M = \begin{bmatrix} -\mathbf{B}\_t^G \mathbf{H}\_M \\ \mathbf{0} \end{bmatrix}.
\end{array}
\tag{78}
$$

One of the main advantages of this solution is that the dimensions of the matrix *E* keep low, 6 × 6, despite the number *G* of mathematical layers and the layer-wise approach to the problem. Furthermore, *E* is the same as that needed for the classical elastic analysis (see in [25]). The vector *<sup>P</sup><sup>M</sup>* adds the hygrometric load as an equivalent mechanical action and sums up the (possible) mechanical load *P*. Solving the system allows getting the bottom displacement components (and their *z* derivatives); Equations (64) and (68) allow then calculating their values at any coordinate in the thickness direction.

**0**

### **4. Results**

This section is of fundamental importance. First of all, it defines the properties of the Functionally Graded (FM) layers, presenting the mechanical and hygrometric properties of their constituents and the law defining their variation in composition. Then, it features two subsections. The first one validates this exact 3D solution for shells embedding layers made of FGM. Validations of new solutions are often conducted comparing the new outputs with those already available in the literature. Besides verifying the accuracy of this solution, this phase helps define how many mathematical layers should be used to consider with confidence the effects of the curvature and those of the FGM law and the order of expansion to use in the exponential matrix calculation. Strengthened by those results, the second subsection presents a set of new results. The effect of the moisture content field is studied on different geometries, featuring different stacking sequences, thickness ratios, and moisture content boundary conditions.

In all the assessments and benchmarks, the FGM layers rely on two constituents: a metal and ceramic. The metal is Monel, 70Ni30Cu, a nickel-based alloy; the ceramic is Zirconia. An estimate of the mechanical properties is given in terms of the bulk modulus *K* and the shear modulus *μ* of the two materials; the moisture expansion coefficients *η* and the moisture diffusion coefficients D are given explicitly:

*Km* <sup>=</sup> 227.24 GPa, *<sup>μ</sup><sup>m</sup>* <sup>=</sup> 65.55 GPa, *<sup>η</sup><sup>m</sup>* <sup>=</sup> <sup>2</sup> <sup>×</sup> <sup>10</sup>−<sup>3</sup> <sup>1</sup> %, <sup>D</sup>*<sup>m</sup>* <sup>=</sup> <sup>10</sup>−<sup>9</sup> kg ms , (79)

*Kc* <sup>=</sup> 125.83 GPa, *<sup>μ</sup><sup>c</sup>* <sup>=</sup> 58.077 GPa, *<sup>η</sup><sup>c</sup>* <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>3</sup> <sup>1</sup> %, <sup>D</sup>*<sup>c</sup>* <sup>=</sup> <sup>10</sup>−<sup>10</sup> kg ms , (80)

The metal and the ceramic phases are denoted by the subscripts *m* and *c*. The estimate of the thermal properties is also given, as they will be necessary in the assessment phase:

$$
\gamma\_m = 15 \times 10^{-6} \frac{1}{K}, \quad \kappa\_m = 25 \frac{\text{W}}{\text{mK}} \tag{81}
$$

$$
\gamma\_{\rm c} = 10 \times 10^{-6} \frac{1}{K'} \quad \kappa\_{\rm c} = 2.09 \frac{\text{W}}{\text{mK}} \tag{82}
$$

*γ* denotes the thermal expansion coefficient, while *κ* the conductivity coefficient. This work assumes that the volume fraction of the ceramic phase follows a power law of order *p*. Introducing the thickness of the FG layer *hFG* and a local thickness coordinate *z*˜*FG* inside it (0 at its bottom, *h* at its top), the FG law takes the following expression:

$$\mathbf{V\_{\mathcal{E}}} = (\tilde{z}\_{\mathrm{FG}} / \mathbf{h\_{FG}})^p \tag{83}$$

At the bottom of the FG layer, where *z*˜*FG* = 0, Equation (83) implies *Vc* = 0, meaning that it is made of metal only. At the top of the FG layer, where *z*˜*FG* = *hFG* , Equation (83) implies *Vc* = 1, meaning that it is made of ceramic only. The bulk and the shear moduli of the FG layer evolve along with the thickness direction following the Mori–Tanaka estimates:

$$\frac{K - K\_m}{K\_c - K\_m} = \frac{V\_c}{1 + (1 - V\_c)\frac{K\_c - K\_m}{K\_m + \frac{4}{3}\mu\_m}}, \quad \frac{\mu - \mu\_m}{\mu\_c - \mu\_m} = \frac{V\_c}{1 + (1 - V\_c)\frac{\mu\_c - \mu\_m}{\mu\_m + f\_m}}, \quad f\_m = \frac{\mu\_m (9K\_m + 8\mu\_m)}{6(K\_m + 2\mu\_m)}\tag{84}$$

The same applies to the effective moisture expansion and moisture diffusion coefficients, by analogy with the estimates of Hatta and Taya for the corresponding thermal properties:

$$\frac{\mathcal{D} - \mathcal{D}\_m}{\mathcal{D}\_c - \mathcal{D}\_m} = \frac{V\_c}{1 + (1 - V\_c)\frac{\mathcal{D}\_c - \mathcal{D}\_m}{3\mathcal{D}\_m}}, \qquad \frac{\eta - \eta\_m}{\eta\_c - \eta\_m} = \frac{\frac{1}{\mathcal{K}} - \frac{1}{\mathcal{K}\_m}}{\frac{1}{\mathcal{K}\_c} - \frac{1}{\mathcal{K}\_m}}\tag{85}$$

The estimates reported in Equation (85) are also valid for the thermal properties, following the parallel of the moisture diffusion coefficient D with the thermal conductivity coefficient *κ*, and that of the moisture expansion coefficient *η* with the the thermal expansion coefficient *γ*.

#### *4.1. Assessments*

The present solution handles several geometries and different load cases. As discussed, it can study plates, cylinders, cylindrical and spherical shells under mechanical, thermal, and hygrometric load. Furthermore, it is not limited to isotropic monolayer structures, but it also handles orthotropic and multi-layered lamination schemes, and it can go up to layers embedding FGM. However, confident use of the model is possible only after validated against established solutions already offered in the literature. The authors did not find hygro-elastic results from exact 3D solutions in the literature which were applied to FGMs. For this reason, the validation process of the model is built by separately validating its different sections against the results provided by other researchers, exploiting the parallel of the moisture content field and the thermal field, and complementing this process with the assistance of 3D FE (Finite Element) models. A static 3D FE model is solved through the Nastran solver SOL101. IsoMesh meshed the geometry of each plate with 3D HEX8 elements; the mesh did not change with the thickness ratio and included 25 elements in the

thickness direction and 30 in both the in-plane ones. Solid elements were necessary to define the mechanical properties evolution in the thickness direction and build significant thermal, hygrometric, and mechanical variables profiles. Depending on the geometry, the coordinate system of the model has one, two, or none curvilinear coordinate. For consistency with the analytical model, the displacement-related boundary conditions are defined on the lateral surfaces of the structure following the harmonic form hypotheses: a pair of displacement coordinates is set 0 on each surface. First, the harmonic thermal/hygrometric field is introduced in an equation, following Equation (1). It is applied to the top and bottom surfaces of the structure; the preprocessor automatically calculates the field values at each node of the external surfaces through their coordinates. The first run solves the thermal part of the problem and returns the temperature (or moisture content) at each node. This field is then applied to a further (and identical) FE model, which solves the elastic part. The thermal, hygrometric, and mechanical properties of the FGM layer are given imagining to split the structures into a number of fictitious layers coinciding with the number of elements in the thickness direction. The properties of each fictitious layer are calculated at its midpoint, following Equations (83)–(85). Needless to say, they are not exact solutions, but they can guide in benchmarking the proposed model.

The first assessment considers a simply-supported one-layered FGM square plate. It investigates different thickness ratios (*a*/*h* = 4, 10, 50) in plates with in-plane dimensions of *a* = *b* = 100 m. The FGM layer relies on a metallic constituent, and a ceramic one, whose mechanical, thermal, and hygrometric properties are defined at the beginning of Section 4. The volume fraction power law considers *p* = 2 as the exponent. An external sovratemperature field acts on the top (*θ<sup>t</sup>* = +1 K) and bottom (*θ<sup>b</sup>* = 0 K) surfaces. The thermal field has a harmonic form, with half-wave numbers *m* = *n* = 1. The reference solution is the asymptotic method of Reddy and Cheng [73], which considers a 3D temperature profile along with thickness direction. Table 1 proposes a pair of results for each thickness ratio at different coordinates along with direction *z* in terms of a displacement, *w* or *u*, and an in-plane shear component, *σαα*. The results show that the 3D shell model always coincides with Reddy and Cheng's asymptotic method, despite the thickness ratio and the considered variable, when the number of mathematical layers is sufficiently high. *NL* = 300, coupled with an order of expansion *N* = 3 for the exponential matrix, always delivers the correct results. Therefore, this assessment verified that the 3D shell model correctly handles the thermomechanical analysis of FGM plates. Furthermore, it simultaneously assessed the 3D FE model, which will be helpful in the following assessment to validate the hygromechanical analysis.

The second assessment is meant to validate the hygroelastic part of the 3D solution for plates embedding an FGM layer. To this end, it considers the previous test case as a reference and removes the thermal field in favor of a hygrometric one. Consequently, it focuses on a simply-supported one-layered FGM square plate; the in-plane dimensions are *a* = *b* = 100 m, the thickness varies to obtain different thickness ratios (*a*/*h* = 4, 10, 50). The FGM layer relies on the same metallic and ceramic constituents, whose volume fraction follows the same power-law with *p* = 2. An external hygrometric field acts on the top (*Mt* = 1.0%) and bottom (*Mb* = 0.5%) surfaces; it has a harmonic form, with half-wave numbers *m* = *n* = 1. The reference results are obtained through the same 3D FE model of the previous assessment, in which the hygrometric field replaces the thermal one. Its previous validation allows considering it as a reliable source for reference results. Table 2 proposes a pair of results for each thickness ratio at different coordinates along with direction *z* in terms of a displacement, *w* or *u*, and an in-plane shear component, *σαα*. Consistent with the previous test case, the results show that the 3D shell model always gives comparable results with the 3D FE model, despite the thickness ratio and the considered variable, when the number of mathematical layers is sufficiently high. *NL* = 300, coupled with an order of expansion *N* = 3 for the exponential matrix, always delivers the correct results. Therefore, this assessment confirmed the capabilities of the 3D shell model in handling the hygromechanical analysis of FGM plates.

**Table 1.** First assessment. One-layered FGM square plate (*a*/*b* = 1), featuring different thickness ratios. The volume fraction power law considers *p* = 2 as the exponent. An external sovra-temperature field acts on the top (*θ<sup>t</sup>* = +1 K) and bottom (*θ<sup>b</sup>* = 0 K) surfaces; *m* = *n* = 1. The reference solution is the asymptotic method of Reddy and Cheng [73], considering a 3D temperature profile along with thickness direction. A 3D FE model is also assessed. The results of the present solution are obtained with *N* = 3 and for a varying number of mathematical layers NL.


The third assessment considers a simply-supported one-layered FGM cylindrical shell. Again, the dimensions of the reference surface are fixed, as in the previous cases; they are *a* = 1 m and *b* = *<sup>π</sup>* <sup>3</sup> *Rβ*, with *R<sup>β</sup>* = 10 m. This test case also considers different thickness ratios *R<sup>β</sup>* = 50, 1000) to evaluate their influence on the performance of the solution. The constituents of the FGM layer are the same as those defined at the beginning of Section 4 and considered in the previous assessments. The power law is also the same, with *p* = 2. The top and bottom external surfaces are subjected to an external sovra-temperature field with amplitudes *θ<sup>t</sup>* = +1 K and *θ<sup>b</sup>* = 0 K. The half-wave numbers of the thermal field are *m* = *n* = 1. The reference solution is a refined 2D layer-wise solution based on the Unified Formulation [74], which considers a 3D temperature profile along with thickness direction. Table 3 proposes six results for each thickness ratio: the transverse displacement w and an in-plane displacement, evaluated at three different coordinates along with direction *z*. The table also assesses a 3D FE model, solved through the Nastran solver, which helps evaluate the hygroelastic analysis of the following assessment. IsoMesh meshed the geometry of each shell with 3D HEX8 elements; the mesh did not change with the thickness ratio and included 25 elements in the thickness direction and 30 in both the in-plane ones. Solid elements were necessary to define the mechanical properties evolution in the thickness direction and build significant thermal and mechanical variables profiles. The results show that the 3D shell model always coincides with the quasi-3D method [74], despite the thickness ratio and the considered variable, when the number of mathematical layers is sufficiently high. *NL* = 300, coupled with an order of expansion *N* = 3 for the exponential matrix, always delivers the correct results. Therefore, this assessment verified that the 3D

shell model correctly handles the thermomechanical analysis of FGM shells. Furthermore, it simultaneously assessed the 3D FE model, which will be helpful in the following assessment to validate the hygromechanical analysis.

**Table 2.** Second assessment. One-layered FGM square plate (*a*/*b* = 1), featuring different thickness ratios. The geometry, materials, and FGM power law are the same as the first assessment. The thermal load is substituted by an external moisture content acting on the top (*Mt* = 1.0%) and bottom (*Mb* = 0.5%) surfaces; *m* = *n* = 1. The reference solution is the 3D FE model, already validated in the previous assessment. The results of the present solution are obtained with *N* = 3 and for a varying number of mathematical layers NL.


The fourth assessment is meant to validate the hygro-elastic part of the 3D solution for shells embedding an FGM layer. To this end, it considers the previous test case as a reference and removes the thermal field in favor of a hygrometric one. Consequently, it focuses on a simply-supported one-layered FGM square shell. The dimensions of the reference surface are *a* = 1 m and *b* = *<sup>π</sup>* <sup>3</sup> *Rβ*, with *R<sup>β</sup>* = 10 m, the thickness varies to obtain different thickness ratios (*R<sup>β</sup>* = 50, 100). The FGM layer relies on the same metallic and ceramic constituents of the previous test cases, whose volume fraction follows the same powerlaw with *p* = 2. An external hygrometric field acts on the top (*Mt* = 1.0%) and bottom (*Mb* = 0.5%) surfaces; it has a harmonic form, with half-wave numbers *m* = *n* = 1. The reference results are obtained through the same 3D FE model of the previous assessment, in which the hygrometric field replaces the thermal one. Its previous validation allows considering it as a reliable source for reference results. Table 4 proposes six results for each thickness ratio: the transverse displacement w and an in-plane displacement, evaluated at three different coordinates along with direction *z*. Consistent with the previous test case, the results show that the 3D shell model always gives comparable results with the 3D FE model, despite the thickness ratio and the considered variable, when the number of mathematical layers is sufficiently high. *NL* = 300, coupled with an order of expansion *N* = 3 for the exponential matrix, always delivers the correct results. Therefore, this assessment

confirmed the capabilities of the 3D shell model in handling the hygromechanical analysis of FGM shells.

**Table 3.** Third assessment. One-layered FGM cylindrical shell (*R<sup>β</sup>* = 10 m), featuring different thickness ratios. The volume fraction power law considers *p* = 2 as the exponent. An external sovra-temperature field acts on the top (*θ<sup>t</sup>* = +1 K) and bottom (*θ<sup>b</sup>* = 0 K) surfaces; *m* = *n* = 1. The reference solution is a refined 2D layer-wise solution based on the Unified Formulation [74], considering a 3D temperature profile along with thickness direction. A 3D FE model is also assessed. The results of the present solution are obtained with *N* = 3 and for a varying number of mathematical layers NL.


**Table 4.** Fourth assessment. One-layered FGM cylindrical shell (*R<sup>β</sup>* = 10 m) featuring different thickness ratios. The geometry, materials, and FGM power law are the same as the third assessment. The thermal load is substituted by an external moisture content acting on the top (*Mt* = 1.0%) and bottom (*Mb* = 0.5%) surfaces; *m* = *n* = 1. The reference solution is the 3D FE model, already validated in the previous assessment. The results of the present solution are obtained with *N* = 3 and for a varying number of mathematical layers NL.


#### *4.2. New Benchmarks*

This section proposes a set of four benchmarks; those new results examine simply supported structures that undergo different moisture content profiles in steady-state conditions. They follow the harmonic form previously defined, precondition to get an exact solution to the problem. The assessments of the previous subsection validated the results of this new

3D shell model when applied to FGM layers: the results converge and are exact when the order of expansion *N* = 3 for the exponential matrix is coupled with a minimum number of *M* = 300 mathematical layer for the through the thickness mechanical properties and curvature approximation. The results of all the following benchmarks consider *N* = 3 and *M* = 300 as an a priori prerequisite for results accuracy.

The first benchmark studies a square plate with a single FGM layer and simplysupported sides. The plate has *a* = *b* = 10 m as in-plane dimensions but comes with several and different thicknesses, which allow the effect of this geometrical parameter to be evaluated. In fact, the thickness ratio goes from *a*/*h* = 2 to *a*/*h* = 100, thus ranging from very thick to very thin plates. In this benchmark, the volume fraction *Vc* of the ceramic phase evolves linearly: *p* = 1 is set in the material law defined through Equation (83). This relation also implies a fully ceramic top surface and a fully metallic bottom one. The moisture content is imposed on the top and the bottom external surfaces; it has a harmonic form on both with amplitudes M*<sup>t</sup>* = 1.0% and M*<sup>b</sup>* = 0.0%, on top and bottom, respectively. The harmonic form of the moisture content has *m* = *n* = 1 as half-wave numbers in *α* and *β* directions, respectively. The elastic and hygroscopic properties of the metallic and ceramic phases are the same as those introduced at the beginning of this section. The FGM nature of the layer makes its mechanical properties evolve through the thickness direction; Figure 3a,b, respectively, shows how the volume fraction *Vc* and the bulk modulus *K* evolve with respect to non-dimensional thickness coordinate *z*˜/*h*. Note that *K* is not linear as *Vc* due to Equation (85). Table 5 and Figures 4 report an extract of the main results. The results in tabular form give the amplitude of some variables of the problem; they reflect the three different ways of evaluating the moisture content. 3D(M*a*) implies the assumed linear moisture content profile; 3D(M*c*, 1*D*) relies on the 1D version of the Fick moisture diffusion equation; finally, 3D(M*c*, 3*D*) relies on a 3D solution of the moisture diffusion problem. The prefix 3D underlines that the elastic part of the solution is three-dimensional in all three models. Such an analysis allows grasping the differences between the three approaches. The 3D(M*c*, 3*D*) shell model considers both the mechanical/hygrometric properties evolution through z and the three-dimensional nature of the problem. As discussed in the previous section, it always delivers an accurate result. Table 5 underlines that the 3D(M*c*, 1*D*) model results get closer to those of the 3D(M*c*, 3*D*) model as the thickness decreases. It considers how the mechanical/hygrometric properties evolve through z but disregards the moisture diffusion through alpha and beta direction, which have a negligible weight in thin structures. The results of the 3D(M*a*) model are always unreliable, as they are built on a moisture content profile that is far from the actual scenario in a layer embedding an FGM. Figure 3c,d further facilitates understanding these concepts; it compares the three moisture content profiles for a thick and a thin plate. In thick structures, the difference between the three profiles is very pronounced: the threedimensionality of the problem and the mechanical/hygrometric properties variability in the thickness direction make the 3D profile differ from the linear assumption. Even the 1D profile differs from the linearity: the hygrometric properties vary through *z*, reflecting on the moisture content at different thickness coordinates. These concepts also apply to thin structures; however, the three-dimensionality of the problem is insignificant, and the 3D and 1D profiles coincide. Figure 4 shows the complete profile of the tree displacement components, two stresses, and a strain. Note that all the quantities evolve with continuity, which is essential as it demonstrates both the graded elastic/hygrometric properties and the correct introduction of the continuity conditions. The transverse stress *σzz* and the transverse shear strain *γβ<sup>z</sup>* satisfy the external mechanical boundary conditions: it equals 0 at both the top and the bottom surfaces as no external load acts on them.

**Figure 3.** First benchmark: one-layered FGM (*p* = 1) square plate with an imposed moisture content on the top and bottom surfaces. The figures show the volume fraction of the ceramic phase, the bulk modulus, and the moisture content profiles for a thick and a thin structure through their thickness. (**a**) Volume fraction of the ceramic phase *Vc*. (**b**) Bulk modulus K. (**c**) Moisture content profile of the *a*/*h* = 2 plate. (**d**) Moisture content profile of the *a*/*h* = 100 plate.

**Figure 4.** First benchmark: one-layered FGM square plate with an imposed moisture content on the top and bottom surfaces. The results are calculated for a moderately thick (*a*/*h* = 4) structure via the 3D(M*c*, 3*D*) model. (**a**) Amplitude of *u* displacement component. (**b**) Amplitude of *v* displacement component. (**c**) Amplitude of *w* displacement component. (**d**) Amplitude of *σzz* stress component. (**e**) Amplitude of *σαα* stress component. (**f**) Amplitude of *γβ<sup>z</sup>* strain component.


**Table 5.** First benchmark: square plate with a single FGM layer subject to an external moisture content applied to the top and bottom surfaces. The results of all the 3D models consider *N* = 3 and *NL* = 300.

The second benchmark focuses on a closed cylinder, featuring a single FGM layer and simply-supported sides. The dimensions of the reference mid-surface, *a* = 2*πR<sup>α</sup>* and *b* = 30 m, are a function of the radii of curvature of the shell, one of which is infinite: *R<sup>α</sup>* = 10 m and *R<sup>β</sup>* = ∞. Different thicknesses have been considered; the thickness ratio *Rα*/*h* is expressed with respect to *R<sup>α</sup>* and ranges from 2 to 100 also in this second case study. The material volume fraction of the ceramic phase is a quadratic function of the thickness coordinate; the material law defined through Equation (83) consider *p* = 2. Given that a single layer is considered, the cylinder is metallic in the inner surface and ceramic in the outer. The moisture content is imposed on the outer external surface, M*<sup>t</sup>* = 1.0%, and on the inner one, M*<sup>b</sup>* = 0.0%. The half wave numbers of both the harmonic forms are the same, *m* = 2 and *n* = 1. The elastic and hygroscopic properties of both the phases introduced previously also apply here. Figure 5a,b, respectively, shows the volume fraction *Vc* and moisture diffusion coefficient D vs. the non-dimensional thickness coordinate *z*˜/*h*. *Vc* follows a power-law of order *p* = 2, D follows Equation (85). Table 6 and Figure 6 summarize an extract of the main results. This second benchmark also reports three different sets of results: the elastic model is the same (prefix 3D), but the moisture content profile follows the different approaches. This leads to models 3D(M*a*), in which the moisture content is a priori assumed, and 3D(M*c*, 1*D*)–3D(M*c*, 3*D*), in which the profile is calculated following a monodimensional or three-dimensional approach. This analysis allows highlighting the distinctions between the three methods. The last one is the only model in which no assumptions are made concerning the three-dimensionality of the problem as the moisture content amplitude derives from Fick's law of diffusion. The results coming from 3D(M*a*) are always wrong because the moisture content evaluation is inaccurate. The differences between 3D(M*c*, 1*D*) and 3D(M*c*, 3*D*) are less pronounced if compared with the previous benchmark and decrease with the thickness. The differences

in the three moisture content profiles are shown in Figure 5c,d for two different cylinders: a thick and a thin one. The discrepancies between the calculated and assumed fields are really pronounced; the 1D and 3D computed profiles do not significantly differ, which is even more true as the thickness ratio increases. Figure 6 shows the profiles of the tree displacement components: two stresses and a strain. There is continuity in all the plots: this qualifies the correct introduction of the continuity conditions and elastic/hygrometric properties grading. No external mechanical loads are applied, and this is coherent with the transverse stress values at the bottom and top surfaces, 0.

**Figure 5.** Second benchmark: one-layered FGM (*p* = 2) cylinder with an imposed moisture content on the top and bottom surfaces. The figures show the volume fraction of the ceramic phase, the diffusion coefficient, and the moisture content profiles for a thick and a thin structure through their thickness. (**a**) Volume fraction of the ceramic phase *Vc*. (**b**) Diffusion coefficient D. (**c**) Moisture content profile of the *Rα*/*h* = 2 cylinder. (**d**) Moisture content profile of the *Rα*/*h* = 100 cylinder.

**Figure 6.** Second benchmark: one-layered FGM closed cylinder with an imposed moisture content on the top and bottom surfaces. The results are calculated for a moderately thin (*Rα*/*h* = 5) structure via the 3D(M*c*, 3*D*) model. (**a**) Amplitude of *u* displacement component. (**b**) Amplitude of *v* displacement component. (**c**) Amplitude of *w* displacement component. (**d**) Amplitude of *σzz* stress component. (**e**) Amplitude of *σββ* stress component. (**f**) Amplitude of *γαβ* strain component.


**Table 6.** Second benchmark: closed cylinder with a single FGM layer subject to an external moisture content applied to the top and bottom surfaces. The results of all the 3D models consider *N* = 3 and *NL* = 300.

The third benchmark considers a cylindrical sandwich shell panel with an FGM core and simply-supported edges. The top and the bottom skin are in line with the FGM law: the top skin is ceramic as the top surface of the core is; at the same time, the bottom skin is metallic as the bottom surface of the core is. *p* = 0.5 is the coefficient for the volume fraction law across the FGM core, which defines how the elastic and hygrometric properties evolve in the thickness direction. The elastic and hygroscopic properties of both the phases already introduced in the assessments also apply here. The radii of curvature are coherent with those proposed in the previous benchmark, *R<sup>α</sup>* = 10 m and *R<sup>β</sup>* = ∞. The dimension of the reference mid-surface in *α* direction is a function of the radius of curvature *Rα* and equals *a* = *<sup>π</sup>* <sup>3</sup> *Rα*; the dimension in the remaining direction *β* is fixed and equals *b* = 30. *m* = 2 and *n* = 0 have been chosen as half-wave numbers for the harmonic form of the moisture content imposed at the bottom and the top of the shell. The amplitude of the external fields discussed so far are as follows: the moisture content amplitude is M*<sup>t</sup>* = 1.0% on the top and M*<sup>b</sup>* = 0.0% on the bottom. This third case study also considers different thickness ratios to evaluate the effects of this parameter; as in the previous case, it is expressed with respect to *R<sup>α</sup>* and ranges from 2 to 100. The volume fraction *Vc* of the ceramic phase runs from 0 to 1 inside the core; it equals 0 inside the bottom skin as it is fully metallic, 1 inside the top coat as it is fully ceramic. This is visible in Figure 7a,b, showing the volume fraction and the shear modulus along with the thickness coordinate z; the shear modulus of the top skin coincides with that of the ceramic; the shear modulus of the bottom skin coincides with that of the metal. The amplitudes of some variables are reported for all the thickness ratios and at different thickness coordinates in Table 7; Figure 8 explores six variables and shows their trend through *z*. The Figures rely on the 3D calculated moisture content profile; the table also reports the results obtained through M*<sup>a</sup>* and M*c*, 1*D*. The differences between the three models can be already seen at the moisture content level and directly

reflect the mechanical quantities. Figure 7c,d compares the three moisture content profiles for a thick and a thin cylindrical shell panel and confirms that the differences are sharp not only at a specific thickness coordinate, but throughout all the thickness. The 3D(M*c*, 1*D*) model results get closer to those of the 3D(M*c*, 3*D*) model as the thickness decreases, and this is clear from the results of Table 7. Figure 8 gives the profile of the tree displacement components, two stresses, and a strain. As in the previous cases, all the quantities evolve with continuity: the mechanical properties are introduced into the model with continuity. The transverse stress *σzz* satisfies the external mechanical boundary conditions: it equals 0 at both the top and the bottom surfaces as no external load acts on them.

**Figure 7.** Third benchmark: cylindrical sandwich shell panel featuring an FGM (*p* = 0.5) core with an imposed moisture content on the top and bottom surfaces. The figures show the volume fraction of the ceramic phase, the shear modulus, and the moisture content profiles for a thick and a thin structure through their thickness. (**a**) Volume fraction of the ceramic phase *Vc*. (**b**) Shear modulus *μ*. (**c**) Moisture content profile of the *Rα*/*h* = 2 shell. (**d**) Moisture content profile of the *Rα*/*h* = 100 shell.

**Figure 8.** Third benchmark: cylindrical sandwich shell panel featuring an FGM (*p* = 0.5) core with an imposed moisture content on the top and bottom surfaces. The results are calculated for a thick (*Rα*/*h* = 4) structure via the 3D(M*c*, 3*D*) model. (**a**) Amplitude of *u* displacement component. (**b**) Amplitude of *v* displacement component. (**c**) Amplitude of *w* displacement component. (**d**) Amplitude of *σzz* stress component. (**e**) Amplitude of *σββ* stress component. (**f**) Amplitude of *γα<sup>z</sup>* strain component.


**Table 7.** Third benchmark: cylindrical sandwich shell panel with an FGM core subject to an external moisture content applied to the top and bottom surfaces. The results of all the 3D models consider *N* = 3 and *NL* = 300.

The fourth and last benchmark proposes a sandwich spherical shell panel, which embeds an FGM core and features simply supported edges. The lamination scheme is analogous to that discussed in the third benchmark: the bottom skin is metallic, and the top ceramic. Then, the volume fraction of the ceramic phase evolves inside the core through the thickness direction following an exponential law with *p* = 0.5 as chosen coefficients. The hygrometric and elastic properties of the sandwich skin are the same proposed in the previous benchmark and assessments for the metallic and ceramic phases, respectively; those of the core follow the volume fraction law. The exponential trend of the volume fraction *Vc* vs. the non-dimensional thickness coordinate *z*˜/*h* can be seen in Figure 9a; for completeness, the evolution of the moisture expansion coefficient *η* through the thickness direction is also given in Figure 10. The spherical shell panel is the only structure among those studied in which both the radii of curvature are non-infinite; furthermore, they take the same value, which equals *R<sup>α</sup>* = *R<sup>β</sup>* = 10 m. Furthermore, the dimensions of the reference mid-surface are the same in *α*, and *β* directions as both are a function of the radii of curvature; it holds *a* = *<sup>π</sup>* <sup>3</sup> *<sup>R</sup><sup>α</sup>* and *<sup>b</sup>* <sup>=</sup> *<sup>π</sup>* <sup>3</sup> *Rβ*. Those dimensions are fixed; however, a wide range of thinner/thicker shells is considered by choosing different thickness ratios: *Rα*/*h* ranges from 2 to 100. The amplitude of the moisture content is imposed on the top and the bottom surfaces; it equals M*<sup>t</sup>* = 1.0% and M*<sup>b</sup>* = 0.0%, respectively. As discussed, the external fields are required to have a harmonic form in order for the problem to be exactly solved; *m* = 2 and *n* = 2 are the half-wave numbers considered in this last case study. Table 8 and Figures 9 and 10 summarize an extract of the main results. This fourth benchmark also reports three different sets of results: the elastic model is the same (prefix 3D), but the moisture content profile follows the different approaches. This analysis highlights the distinctions between the three methods. The 3D one is the only model in which no assumptions are made concerning the three-dimensionality of

the problem as the moisture content amplitude derives from Fick's law of diffusion. The results coming from 3D(M*a*) are always wrong because the moisture content evaluation is inaccurate. Considerable differences are present between the calculated and assumed fields. The moisture content profiles of Figure 9c,d once again demonstrate that the 1D and 3D moisture fields get closer in thin structures; despite the thickness, they always differ from the assumed profile, which completely disregards the physics of the problem. This reflects on the results in terms of displacements, strains, and stresses: the differences are high, and 3D(M*a*) does not provide a reasonable estimate. 3D(M*c*, 1*D*) provides acceptable results, but only when the shell is sufficiently thin. As in the previous cases, three displacement components, two stresses, and a strain are shown in their entirety along with the thickness direction. Figure 10 further qualifies the correct introduction of the continuity conditions, elastic/hygrometric properties grading, and mechanical boundary conditions. The transverse stresses *σβ<sup>z</sup>* and *σzz* satisfy the external mechanical boundary conditions: they equal 0 at both the top and the bottom surfaces as no external load acts on them. All the quantities are continuous throughout the thickness; this qualifies the division into fictitious layers: they are thin enough to describe the mechanical properties evolution with continuity.

**Figure 9.** Fourth benchmark: spherical sandwich shell panel featuring a FGM (*p* = 0.5) core with an imposed moisture content on the top and bottom surfaces. The figures show the volume fraction of the ceramic phase, the moisture expansion coefficient, and the moisture content profiles for a thick and a thin structure through their thickness. (**a**) Volume fraction of the ceramic phase *Vc*. (**b**) Moisture expansion coefficient *η*. (**c**) Moisture content profile of the *Rα*/*h* = 2 shell. (**d**) Moisture content profile of the *Rα*/*h* = 100 shell.

**Figure 10.** Fourth benchmark: spherical sandwich shell panel featuring a FGM (*p* = 0.5) core with an imposed moisture content on the top and bottom surfaces. The results are calculated for a thin (*Rα*/*h* = 50) structure via the 3D(M*c*, 3*D*) model. (**a**) Amplitude of *u* displacement component. (**b**) Amplitude of *v* displacement component. (**c**) Amplitude of *w* displacement component. (**d**) Amplitude of *σzz* stress component. (**e**) Amplitude of *σβ<sup>z</sup>* stress component. (**f**) Amplitude of *αα* strain component.


**Table 8.** Fourth benchmark, sandwich spherical shell panel with an FGM core subject to an external moisture content applied to the top and bottom surfaces. The results of all the 3D models consider *N* = 3 and *NL* = 300.

#### **5. Conclusions**

The authors proposed a closed-form 3D shell solution that handles the hygro-elastic stress analysis of plates, cylinders, cylindrical shells, and spherical shells while embedding Functionally Graded Material (FGM) layers. First, the author imposed the external moisture content on the top and the bottom surfaces. The moisture conditions act in steady state as an external load; calculating the moisture content profile is a prerequisite. The authors showed that three approaches might be used to determine the moisture content profile along the thickness direction and coupled them with a consolidated elastic solution. The results demonstrated the importance of a correct moisture content profile evaluation in the thickness direction. The 3D Fick's law of diffusion is the only way to obtain exact results when the structures embed FG layers; it is also necessary when the structures are sufficiently thick. On the other hand, the 1D Fick's law of diffusion comes closer to it only when structures are thin; as a rule of thumb, the results of the two models are almost coinciding only from a thickness ratio of 50. The problem relies on a set of differential equations in the thickness direction. The authors demonstrated that the exponential matrix method is a reliable way to solve it, provided that the structures are divided into a sufficiently high number of mathematical layers. This layer-wise approach is critical to get a reliable description of the material properties grading; as a rule of thumb, 300 mathematical layers always deliver the correct results. This achievement is general and does not depend on the geometry, FGM law, and lamination sequence/scheme.

**Author Contributions:** Data curation, R.T.; Formal analysis, S.B. and R.T.; Methodology, S.B.; Software, S.B.; Supervision, S.B.; Validation, R.T.; Writing—original draft, R.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com