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

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 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,16,17,18,19,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,26,27,28,29,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,32,33,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,58,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,61,62,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.

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 $M^k\left(z\right)$ indicates the moisture content amplitude. Referring to Figure 1 and Figure 2: m and n are the half-wave numbers in the two in-plane directions $\alpha$ and $\beta$, respectively. a and b are the shell dimensions referred to the mid-surface ${\Omega }_0$; they allow calculating the terms $\overline{\alpha }=\frac{m\pi }{a}$ and $\overline{\beta }=\frac{n\pi }{b}$. Note that the harmonic form of the moisture content reduces its assessment to its profile along with the thickness direction z. What happens in $\alpha$ and $\beta$ directions is already defined through the sinusoidal functions.

The hygrometric boundary conditions are the moisture content amplitudes at the bottom and the top of the shell, $M_b$, and $M_t$, 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:

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.

The 3D equilibrium equations for shells are written in the orthogonal mixed curvilinear coordinate system $(\alpha ,\beta ,z)$ shown in Figure 1 and Figure 2. These equations are modified using 3D constitutive equations for Functionally Graded Materials (FGMs) and general geometrical relations for shells in $(\alpha ,\beta ,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 $\alpha$ and $\beta$ 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 $\mathcal{M}(\alpha ,\beta ,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 ($\alpha$,$\beta$,z) are the algebraic summation of mechanical elastic parts (subscript m) and hygroscopic parts (subscript $\mathcal{M}$). The $6\times 1$ vector (${ϵ}_{\alpha \alpha }^k,{ϵ}_{\beta \beta }^k,{ϵ}_{zz}^k,{\gamma }_{\beta z}^k,{\gamma }_{\alpha z}^k,{\gamma }_{\alpha \beta }^k$) for each k layer is defined as

In the hygro-elastic strains, the three displacement components $u^k$, $v^k$, and $w^k$ and the scalar moisture content ${\mathcal{M}}^k$ are defined in the reference system ($\alpha$, $\beta$, z). Moisture expansion coefficients ${\eta }_{\alpha }^k\left(z\right)$, ${\eta }_{\beta }^k\left(z\right)$ and ${\eta }_z^k\left(z\right)$ 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 ($\alpha$, $\beta$, z) starting from the moisture expansion coefficients ${\eta }_1^k\left(z\right)$, ${\eta }_2^k\left(z\right)$ and ${\eta }_3^k\left(z\right)$ 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 displacement amplitudes are indicated as $(U^k\left(z\right),V^k\left(z\right),W^k\left(z\right))$. As already described, m and n are the half-wave numbers, a and b the mid-surface dimensions, and $\overline{\alpha }=\frac{m\pi }{a}$ and $\overline{\beta }=\frac{n\pi }{b}$.

The three-dimensional equilibrium equations for spherical shells having constant radii of curvature $R_{\alpha }=R_{\beta }$ and a total number $N_L$ of physical k (either classical or FGM) layers are

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$ is the stress vector having a $6\times 1$ dimension, ${\mathit{C}}^k\left(z\right)$ is the $6\times 6$ elastic coefficient matrix and it could depend on the z coordinate in the case of a kth 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^{\circ }$ or ${90}^{\circ }$ to obtain $C_{16}^k\left(z\right)=C_{26}^k\left(z\right)=C_{36}^k\left(z\right)=C_{45}^k\left(z\right)=0$ in the structural reference system ($\alpha$, $\beta$, z). Therefore,

The explicit form of the constitutive equations develops thanks to the introduction of the geometrical Equations (25)–(30) into the constitutive Equation (37): partial derivatives ($\frac{\partial }{\partial \alpha }$), ($\frac{\partial }{\partial \beta }$), and ($\frac{\partial }{\partial z}$) are indicated in Equations (39)–(44) through subscripts ($,\alpha$), ($,\beta$), and ($,z$), respectively. Terms ${\xi }_{\alpha }^k\left(z\right)$, ${\xi }_{\beta }^k\left(z\right)$, and ${\xi }_z^k\left(z\right)$ designate the hygro-mechanical coupling coefficients in the structural reference system ($\alpha$, $\beta$, z), and they could depend on z in the case of FGM layers:

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:

The above equations have no constant coefficients because the parametric coefficients $H_{\alpha }$ and $H_{\beta }$ 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_{\alpha }$ and $H_{\beta }$ 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_s^j$ (s from 1 to 19) and $L_r^j$ (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 $\alpha$ and $\beta$ 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^j$, $V^j$, $W^j$ and their derivatives in z. Redoubling the variables as proposed in [71,72] allows reducing the system into a first-order one in z. In each j layer, the $3\times 1$ vector of unknowns ($U^j$, $V^j$, $W^j$) becomes a $6\times 1$ unknown vector ($U^j$, $V^j$, $W^j$, $U^{j^{{}^{\prime }}}$, $V^{j^{{}^{\prime }}}$, $W^{j^{{}^{\prime }}}$); superscript ${}^{\prime }$ indicates derivatives performed to z (also written as $\frac{\partial }{\partial z}$). Terms $M^j$ and $M^{j^{\prime }}$ can be considered as known terms because they can be opportunely calculated, as will be demonstrated in the next section:

By defining vectors ${\mathit{U}}^j={\left[U^j{V}^j{W}^j{U}^{j^{{}^{\prime }}}{V}^{j^{{}^{\prime }}}{W}^{j^{{}^{\prime }}}\right]}^T$, ${{\mathit{U}}^j}^{{}^{\prime }}=\frac{\partial {\mathit{U}}^j}{\partial z}$ and ${\mathit{M}}^j$$={\left[M^j{M}^{j^{{}^{\prime }}}0000\right]}^T$ (where T means the transpose of a vector), the following compact form is possible: the above equation, thanks to the definitions ${\mathit{A}}^{{*}^j}={\mathit{D}}^{j^{-1}}{\mathit{A}}^j$ and ${\mathit{L}}^{{*}^j}={\mathit{D}}^{j^{-1}}{\mathit{L}}^j$, can be rewritten as

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 in a jth mathematical layer, $a_M^j$ and $b_M^j$ 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. ${\tilde{z}}^j$ 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^j$ of the same mathematical layer.

The system of first-order differential equations in $\tilde{z}$ or z shown in Equation (54) is not homogeneous because the hygroscopic term ${\mathit{L}}^{{*}^j}{\mathit{M}}^j$ depends on ${\tilde{z}}^j$ or $z^j$. In the case of a generic system of non-homogeneous first-order differential equations having an unknown $G\times 1$ vector $\mathit{x}$, a $G\times G$ |$\mathit{A}$ matrix containing constant coefficients, and a known function vector $\mathit{f}\left(t\right)={[f_1\left(t\right)...f_G\left(t\right)]}^T$,we can write a possible solution of the Equation (56) can be obtained via the exponential matrix method:

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

Therefore, Equation (54) can be rewritten as ${\mathit{M}}^{{*}^j}$ includes only linear and known functions in ${\tilde{z}}^j$ coordinate. The Equation (59) can be solved through the exponential matrix method:${\mathit{A}}^{{**}^j}={\mathrm{e}}^{\left({\mathit{A}}^{{*}^j}{h}^j\right)}$ and ${\mathit{L}}^{{**}^j}={\int }_0^{h^j}{\mathrm{e}}^{\left({\mathit{A}}^{{*}^j}(h^j-s)\right)}{\mathit{M}}^{{*}^j}\left(s\right)ds$ must be opportunely defined for each j layer having thickness $h^j$. 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^j$: where $\mathit{I}$ is the $6\times 6$ identity matrix, and the integral given in Equation (60) can be calculated in each j layer having thickness $h^j$ 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 where ${\mathit{U}}_t^j$ indicates ${\mathit{U}}^j\left(h^j\right)$ and it is defined at the top t of each j layer, and ${\mathit{U}}_b^j$ indicates ${\mathit{U}}^j\left(0\right)$ and it is defined at the bottom of each j layer. In this way, Equation (63) is defined as

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:

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

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 ${\mathit{U}}^j$ and ${\mathit{M}}^j$, 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:

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:

Equation (68) expresses the link between the displacements and their z derivatives, calculated at the bottom of the jth 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:

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 $\alpha$, $\beta$, and z. They are grouped into vectors $\mathit{P}={\left(P_{\alpha }{P}_{\beta }{P}_z\right)}^T$; the superscript G 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: 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 $\mathit{B}$ convert the mechanical and/or hygrometric loads into displacements. ${\mathit{B}}_t^G$, in Equation (71) handles the top (t) of the last layer (G); ${\mathit{B}}_b^1$, 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 ${\mathit{U}}_t^G={\mathit{U}}^G\left(h^G\right)$ as a function of ${\mathit{U}}_b^1={\mathit{U}}^1\left(0\right)$. 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: