Hygro-Elastic Coupling in a 3D Exact Shell Model for Bending Analysis of Layered Composite Structures

Abstract

In this work, a 3D fully coupled hygro-elastic model is proposed. The moisture content profile is a primary variable of the model’s displacements. This generic fully coupled 3D exact shell model allows the investigations into the consequences arising from moisture content and elastic fields in terms of stresses and deformations on different plate and shell configurations embedded in composite and laminated layers. Cylinders, plates, cylindrical and spherical shells are analyzed in the orthogonal mixed curvilinear reference system. The 3D equilibrium equations and the 3D Fick diffusion equation for spherical shells are fully coupled in a dedicated system. The main advantage of the orthogonal mixed curvilinear coordinates is related to the degeneration of the equations for spherical shells to simpler geometries thanks to basic considerations of the radii of curvature. The exponential matrix method is used to solve this fully coupled model based on partial differential equations in the thickness direction. The closed-form solution is related to simply supported sides and harmonic forms for displacements and the moisture content. The moisture content amplitudes are directly applied at the top and bottom outer faces through steady-state hypotheses. The final system is based on a set of coupled homogeneous second-order differential equations. A first-order differential equation system is obtained by redoubling the number of variables. The moisture field implications are evaluated for the static analysis of the plates and shells in terms of displacement and stress components. After preliminary validations, new benchmarks are proposed for several thickness ratios, geometrical and material data, lamination sequences and moisture values imposed at the external surfaces. In the proposed results, there is clearly accordance between the uncoupled hygro-elastic model (where the 3D Fick diffusion law is separately solved) and this new fully coupled hygro-elastic model: the differences between the investigated variables (displacements, moisture contents, stresses and strains) are always less than 0.3%. The main advantages of the 3D coupled hygro-elastic model are a more compact mathematical formulation and lower computational costs. Both effects connected with the thickness layer and the embedded materials are included in the conducted hygro-elastic analyses.

1. Introduction

Modern composite aircraft tend to absorb moisture from the humid air that can accelerate the degradation of structural performances. This effect is more remarkable for the matrix material than for the fibers. The hygrothermal effects can be analyzed as the degradation of material properties and/or as hygroscopic and thermal load applications on structures to evaluate their typical bending behaviors [1,2,3,4,5]. In the literature, several works have been related to the hygro-elastic plate and shell bending analyses of multilayered and composite configurations. Both numerical and analytical solutions are proposed. In order to understand the novelties implemented in this new formulation, two different sections are discussed. The first one is related to analytical models in 1D, 2D, and 3D form. The second one is about 1D, 2D and 3D numerical approaches.

For what concerns the exact models, Bouazza and Zenkour [6] proposed the analytical approach of buckling of cross-ply and angle-ply composite beams under hygro-thermo-mechanical loadings. The steady-state dynamic problem of a composite double-beam system under hygrothermal environments was proposed in [7] in a systematic way. The governing equations were derived using the Euler–Bernoulli beam hypotheses considering the expansion caused by temperature and moisture conditions. Li et al. [8] proposed a rotating composite thin-walled beam structure introduced in the hygrothermal environment that can influence the material properties of composites and modify the mechanical behavior. An analytical model was presented by Tsokanas and Loutas [9] to calculate the first two mode components of delamination between two sub-laminates, modeled as Timoshenko beams, considering the residual hygrothermal stress effects. Yu and Sun [10] presented the large deformation post-buckling analyses of a linear-elastic beam in an hygrothermal environment, considering analytical approximate solutions for nonlinear problems. Boukert et al. [11] proposed the investigation of thick composite laminates using a high-order method and hygrothermal stress calculations. The moisture distribution was calculated along the z direction of the laminate with the Fick equation. Brischetto [12] analyzed the hygrothermal loading effects in the bending of multilayered composite plates with the Fick moisture diffusion law and the Fourier heat conduction equation. Refined two-dimensional models were used to evaluate these effects. The extension to shell geometries was presented in [13], where the refined two-dimensional models for the static hygro-thermo-elastic analysis of multilayered composite and sandwich shells were discussed. Chien et al. [14] analyzed the stability equations for perfect/imperfect composite plates subjected to hygro-thermal loads. In [15], a hygro-magneto-electro wave propagation analysis was presented for axially moving circular cylindrical nanoshells. Panduro and Mantani [16] discussed the hygro-thermo-mechanical behavior of multilayered composite and sandwich plates based on a layer-wise approach. Hygro-thermo-mechanical loads were applied considering the various profiles obtained via Fick moisture diffusion and Fourier heat conduction laws. The same authors presented the introduction of five different nonpolynomial shear strain shape functions in [17] in the models already proposed in [16]. Peron et al. [18] proposed the hygro-mechanical behavior of a bio-sourced composite material through the characterization of its moisture diffusion and elastic properties. Zenkour and El Sharany [19] showed the vibration of a simply supported rectangular composite plate with four actuating magnetostrictive layers in an hygrothermal environment. The plate was supported by the two-parameter elastic Pasternak’s foundations. The same authors presented in [20] the vibration study of a simply supported smart sandwich plate embedded in an elastic substrate. The sandwich plate contained fiber-reinforced and magnetostrictive layers and a core made of a viscoelastic material. Brischetto and Torre showed in [21] the stress analysis of single-layered and multilayered composite and sandwich plates and shells under steady state moisture conditions. This 3D solution was based on an exact layer-wise approach where the formulation was unique for different geometries. The moisture conditions were externally defined in an uncoupled way at the outer surfaces and evaluated in the thickness direction following three different procedures. The same authors presented in [22] the same model for the hygrometric loading effects analysis of single-layered and multilayered composite plates and shells embedding FGM layers.

In the framework of the numerical models, Jena et al. [23] investigated the vibration characteristics of a nanobeam embedded in a Winkler–Pasternak elastic foundation using the Haar wavelet method. The nanobeam was subjected to a longitudinal magnetic field in hygroscopic and thermal environments. Motamedian and Kulachenko [24] presented a method for modeling the moisture or thermal expansion of interconnected fiber networks using only beam elements. Oceallaigh et al. [25] proposed a fully coupled moisture–displacement finite element model to predict the viscoelastic behavior of reinforced timber elements when under loads and immersed in a relative humidity field. Zhou et al. [26] presented a study of piezoelectric beams under multi-physical hygro-thermo-electro-mechanical fields. Akbas proposed in [27] a post-buckling analysis of laminated composite beams under hygrothermal effects using the finite element method and the first-order shear beam theory. Amoushahi and Goodarzian [28] investigated the effects of hygrothermal conditions on free vibration frequencies and the buckling loads of composite laminated plates. Bandyopadhyay et al. showed in [29] a finite element-based method to investigate the hygrothermal effects on the transient dynamic response of a delaminated composite pre-twisted conical shells using an eight-node shell element based on Mindlin theory. A finite element-based method was developed in [30] to study the influence of elevated temperature and moisture absorption on the free vibration behavior of rotating pre-twisted sandwich conical shells. Kolahchi et al. [31] presented a study on the post-buckling behaviors of three different types of defective quadrilateral single-layered graphene sheets subjected to temperature, moisture, and in-plane magnetic loads. Kolahchi and Kolahdouzan analyzed in [32] the dynamic stability of the viscoelastic defective single-layered graphene sheet already seen in [31]. Lal et al. [33] investigated the effect of random system properties on the transverse nonlinear central deflection of composite spherical shell panels subjected to hygro-thermo-mechanical loadings. The higher-order shear deformation theory and von Karman nonlinear kinematics were used for the formulation. A control method was presented in [34] for the vibration suppression of the piezoelectric composite cantilever rectangular plates subjected to the aerodynamic forces in a hygrothermal environment thanks to a classical composite plate theory and Hamilton principle. Ma et al. [35] experimentally and analytically investigated the hygrothermal effects on the shear behaviors of composite stiffened panels. A moisture absorption test and a new moisture absorption model based on the non-Fickian model were included. Rajanna et al. [36] presented a model for a panel and stiffener, under severe external conditions, adopting nine-node heterosis plate elements and three-node beam elements. Rajia et al. showed in [37] the active stiffening and compensation analyses on the dynamic behavior of piezo-hygro-thermo-elastic laminates using a coupled piezoelectric finite element formulation involving a hygrothermal strain field. Gholami et al. [38] investigated the effects of hygrothermal conditions on the elastic properties of polymeric composite materials using a micromechanical degradation model based on the Fick second law via an ABAQUS scripting micro modeling and on an ABAQUS parallel finite element analysis. A 3D micromechanical model studying stresses due to moisture diffusion through a unidirectional fiber-reinforced polymeric matrix composite was analyzed in [39]. Knarud and Geving [40] proposed several benchmarks for hygrothermal simulation models developed in the COMSOL Multiphysics software. Vinyas and Kattimani [41] showed the static response of magneto-electro-elastic plates subjected to hygrothermal loads using a numerical formulation derived from the principle of total potential energy, taking into account the thermal and hygroscopic field effects. The coupled finite element equilibrium equations were solved using the condensation procedure. Ye et al. [42] proposed the reduction study of the strength and stiffness of composites due to the rise in moisture concentration and temperature. The 3D hygrothermal vibration analysis of multilayered cylindrical shells under general boundary conditions was conducted. Yi and Hilton [43] investigated the temperature and moisture content effects on the viscoelastic responses of composites using a numerical approach. Time-dependent strains and stresses were presented for the involved composite plates. In [44], Li et al. proposed the investigation of the active control effect on the bending and vibration responses of the magnetorheological elastomer (MRE) multifunctional grid composite sandwich plates using the first-order shear deformation theory, the energy principle, the Ritz method, and the Duhamel integral approach. A similar study for the static and dynamic performances of sandwich plates with the MRE core was proposed in [45]. In ref. [46], a finite element model updating (FEMU) procedure is shown to discover the best creep parameters for filament-wound cylinders under radial compression under harsh environmental conditions. Three different winding angles are considered, each under three different hygrothermal conditions.

This new 3D coupled hygro-elastic shell model could analyze different geometries: spherical and cylindrical panels, plates and cylinders. It is formulated in closed-form using the simply supported hypotheses for the sides and the harmonic forms for displacements and moisture contents. The 3D equilibrium equations and the 3D diffusion equation for multilayered composite structures were solved using the exponential matrix method along the thickness direction. A layer-wise approach was employed. The present 3D exact coupled hygro-elastic shell model can be seen as the general case of the pure mechanical model already developed by Brischetto in [47,48,49] in the case of free vibration and bending analyses of multilayered composite structures. The addition of the moisture-related equation in orthogonal mixed curvilinear coordinates (see [50,51,52,53]) to 3D shell equilibrium relations generates a set of homogeneous differential equations that can be solved via the procedure shown in [54,55]. The new results are presented in terms of displacements, stresses, and moisture profiles. They can be used as benchmarks for the development and testing of new 3D, 2D, and 1D numerical models for the hygro-elastic stress analysis of composite and sandwich structures. The main novelty with respect to the 3D uncoupled hygro-elastic shell model proposed in [21,22] is the full coupling between the displacement field and the hygroscopic field. Therefore, the moisture content profile is now directly obtained from the system solution and not “a priori” separately calculated. A similar procedure was already used by the authors in [56,57] in the case of full coupling between the displacement field and the thermal field for the 3D bending analysis of composite and functionally graded plates and shells.

2. Geometrical and Constitutive Equations for Hygro-Elastic Problems

The 3D fully coupled hygro-elastic exact model considers plate and shell structures via a general formulation for different geometries. The innovative feature with respect to the similar 3D model proposed by Brischetto and Torre in [21] is the full coupling between the elastic and moisture fields. Therefore, the primary variables of the formulation are the three displacement components u, v, and w and the scalar moisture content $\mathcal{M}$. The appropriate constitutive and geometrical equations will be introduced in 3D equilibrium equations and 3D Fick moisture diffusion relation for the spherical shells. The multilayered shell is subjected to a moisture content $\mathcal{M}$($\alpha$, $\beta$, z) at the outer faces of the structure. For each k physical layer, the geometrical relations written in an orthogonal mixed curvilinear reference system ($\alpha$, $\beta$, z) are:

$${\mathit{ϵ}}^k=(\Delta \left(z\right)+\mathit{G}\left(z\right)){\mathit{u}}^k-{\mathit{\eta }}^k{\mathcal{M}}^k={\mathit{ϵ}}_u^k-{\mathit{ϵ}}_{\mathcal{M}}^k$$

(1)

where ${\mathit{ϵ}}^k$ is the $6\times 1$ strain vector, $\Delta \left(z\right)$ is a $6\times 3$ matrix containing the derivatives for the shell geometry, $\mathit{G}\left(z\right)$ is a $6\times 3$ matrix that includes pure geometrical terms (radii of curvature $R_{\alpha }$ and $R_{\beta }$ and parametric coefficients $H_{\alpha }$ and $H_{\beta }$). ${\mathit{u}}^k$ is the $3\times 1$ displacement vector, ${\mathit{\eta }}^k$ is the $6\times 1$ vector containing the moisture expansion coefficients evaluated in the mixed orthogonal reference system, and $\mathcal{M}$ is the scalar moisture content. The previous matrices and vectors are:

$$\begin{array}{cc}\hfill & {\mathit{ϵ}}^k=\left[\begin{array}{c}{ϵ}_{\alpha \alpha }^k\\ {ϵ}_{\beta \beta }^k\\ {ϵ}_{zz}^k\\ {\gamma }_{\beta z}^k\\ {\gamma }_{\alpha z}^k\\ {\gamma }_{\alpha \beta }^k\end{array}\right],\Delta \left(z\right)=\left[\begin{array}{ccc}\frac{\partial }{\partial \alpha }\frac{1}{H_{\alpha }\left(z\right)}& 0& 0\\ 0& \frac{\partial }{\partial \beta }\frac{1}{H_{\beta }\left(z\right)}& 0\\ 0& 0& \frac{\partial }{\partial z}\\ 0& \frac{\partial }{\partial z}& \frac{\partial }{\partial \beta }\frac{1}{H_{\beta }\left(z\right)}\\ \frac{\partial }{\partial z}& 0& \frac{\partial }{\partial \alpha }\frac{1}{H_{\alpha }\left(z\right)}\\ \frac{\partial }{\partial \beta }\frac{1}{H_{\beta }\left(z\right)}& \frac{\partial }{\partial \alpha }\frac{1}{H_{\alpha }\left(z\right)}& 0\end{array}\right],\mathit{G}\left(z\right)=\left[\begin{array}{ccc}0& 0& \frac{1}{H_{\alpha }\left(z\right)R_{\alpha }}\\ 0& 0& \frac{1}{H_{\beta }\left(z\right)R_{\beta }}\\ 0& 0& 0\\ 0& -\frac{1}{H_{\beta }\left(z\right)R_{\beta }}& 0\\ -\frac{1}{H_{\alpha }\left(z\right)R_{\alpha }}& 0& 0\\ 0& 0& 0\end{array}\right],\hfill \\ & {\mathit{u}}^k=\left[\begin{array}{c}{u}^k\\ v^k\\ w^k\end{array}\right],{\mathit{\eta }}^k=\left[\begin{array}{c}{\eta }_{\alpha }^k\\ {\eta }_{\beta }^k\\ {\eta }_z^k\\ 0\\ 0\\ 0\end{array}\right],\hfill \end{array}$$

(2)

where parametric coefficients $H_{\alpha }\left(z\right)$, $H_{\beta }\left(z\right)$ and $H_z$ are:

$$H_{\alpha }\left(z\right)=\left(1+\frac{z}{R_{\alpha }}\right)=\left(1+\frac{\tilde{z}-h/2}{R_{\alpha }}\right),H_{\beta }\left(z\right)=\left(1+\frac{z}{R_{\beta }}\right)=\left(1+\frac{\tilde{z}-h/2}{R_{\beta }}\right),H_z=1.$$

(3)

$H_{\alpha }$ and $H_{\beta }$ are functions of the thickness coordinate z (which varies from $-h$/2 to +h/2 with the origin located in the correspondence of the ${\mathrm{\Omega }}_0$ surface) or $\tilde{z}$ (which varies from 0 to h with the zero located at the bottom surface). They represent the curvature terms of the shell in the two in-plane directions $\alpha$ and $\beta$.

It is possible to define the moisture content $\mathcal{M}$ in a non-dimensional form (or in percentage by simply multiplying the non-dimensional form for 100) as:

$$\mathcal{M}=\frac{W-W_d}{W_d}=\frac{W_d+W_c-W_d}{W_d}=\frac{W_c}{W_d},$$

(4)

where W is the mass of the moist material that is composed by the mass of the dry material ($W_d$) plus the mass of the moisture ($W_c$) present in the considered material. The mass of the moisture present in the material can be computed by integrating the moisture concentration c (in kg/m${}^3$) in the volume of the material as follows:

$$W_c={\int }_{V}c·dV.$$

(5)

The same idea is valid for the mass of the dry material that is calculated whilst integrating the mass density of the dry material in the volume:

$$W_d={\int }_V{\rho }_d·dV,$$

(6)

where ${\rho }_d$ (kg/m${}^3$) is the mass density of the dry material. Introducing Equations (5) and (6) in Equation (4), the definition of $\mathcal{M}$ can be rewritten as:

$$\mathcal{M}=\frac{W_c}{W_d}={\int }_V\frac{c·dV}{{\rho }_d·dV}=\frac{cV}{{\rho }_{d}V}=\frac{c}{{\rho }_d}.$$

(7)

The moisture content is in non-dimensional form because it is computed as a ratio between two mass density values in kg/m${}^3$. It should be noted that, if $\mathcal{M}$ is in non-dimensional form, the related moisture expansion coefficients ${\eta }_i$ have to be written in non-dimensional form; if the moisture content is written in percentage, the coefficients ${\eta }_i$ must be divided by 100.

The constitutive equations for the proposed problem permit the linking between the six-strain components with the six stress components via the $6\times 6$ elastic coefficient matrix $\mathit{C}$. The constitutive equations (Hooke law) are valid for each k physical layer:

$${\mathit{\sigma }}^k={\mathit{C}}^k{\mathit{ϵ}}^k={\mathit{C}}^k({\mathit{ϵ}}_u^k-{\mathit{ϵ}}_{\mathcal{M}}^k),$$

(8)

where the stress vector ${\mathit{\sigma }}^{\mathit{k}}=$[${\sigma }_{\alpha \alpha }^k{\sigma }_{\beta \beta }^k{\sigma }_{zz}^k{\sigma }_{\beta z}^k{\sigma }_{\alpha z}^k{\sigma }_{\alpha \beta }^k{]}^T$ has dimensions of $6\times 1$ and the elastic coefficient matrix ${\mathit{C}}^k$ has dimensions of $6\times 6$. The matrix of elastic coefficients in the structural reference system for the orthotropic case, with lamination angles of $0^{\circ }$ or ${90}^{\circ }$, has the following form:

$${\mathit{C}}^k=\left[\begin{array}{cccccc}{C}_{11}^k& C_{12}^k& C_{13}^k& 0& 0& 0\\ C_{12}^k& C_{22}^k& C_{23}^k& 0& 0& 0\\ C_{13}^k& C_{23}^k& C_{33}^k& 0& 0& 0\\ 0& 0& 0& C_{44}^k& 0& 0\\ 0& 0& 0& 0& C_{55}^k& 0\\ 0& 0& 0& 0& 0& C_{66}^k\end{array}\right].$$

(9)

The solution of the problem is only developed for orthotropic angles $0^{\circ }$ or ${90}^{\circ }$ because the impositions $C_{16}^k=C_{26}^k=C_{36}^k=C_{45}^k=0$ permit to solve the mathematical problem in a closed form. The substitution of Equation (1) into Equation (8) allows to link the stresses with the primary variables of the problem ${\mathit{u}}^k$ and ${\mathcal{M}}^k$. The explicit form of Equation (8) defines the hygro-elastic coupling coefficients ${\xi }_{\alpha }^k$, ${\xi }_{\beta }^k$ and ${\xi }_z^k$ in the mixed orthogonal structural reference system ($\alpha$, $\beta$, z). These are defined as follows:

$$\begin{array}{c}\hfill {\xi }_{\alpha }^k=C_{11}^k{\eta }_{\alpha }^k+C_{12}^k{\eta }_{\beta }^k+C_{13}^k{\eta }_z^k,\end{array}$$

(10)

$$\begin{array}{c}\hfill {\xi }_{\beta }^k=C_{12}^k{\eta }_{\alpha }^k+C_{22}^k{\eta }_{\beta }^k+C_{23}^k{\eta }_z^k,\end{array}$$

(11)

$$\begin{array}{c}\hfill {\xi }_z^k=C_{13}^k{\eta }_{\alpha }^k+C_{23}^k{\eta }_{\beta }^k+C_{33}^k{\eta }_z^k.\end{array}$$

(12)

3. Three-Dimensional Coupled Hygro-Elastic Governing Equations

Figure 1 shows the different proposed geometries in an orthogonal mixed curvilinear reference system $(\alpha ,\beta ,z)$. The global thickness of the structure h has a constant value for each configuration. The origin of the orthogonal mixed curvilinear reference system is placed in the corner. $\alpha$ and $\beta$ are parallel to the lateral curvilinear sides of the shell and located on the middle surface ${\mathrm{\Omega }}_0$. ${\mathrm{\Omega }}_0$ is the reference surface for the evaluation of all geometrical parameters. z is the coordinate in the thickness direction and it is orthogonal to ${\mathrm{\Omega }}_0$ and directed from the bottom to the top surface. Angle $\psi$ shows the curvature angle of the two in-plane dimensions.

For the present formulation, the governing equations consist of three 3D equilibrium equations and the 3D Fick diffusion equation for spherical shells. The explicit forms of the governing equations, valid for each k physical layer, are written as:

$$\begin{array}{cc}\hfill & H_{\beta }\left(z\right)\frac{\partial {\sigma }_{\alpha \alpha }^k}{\partial \alpha }+H_{\alpha }\left(z\right)\frac{\partial {\sigma }_{\alpha \beta }^k}{\partial \beta }+H_{\alpha }\left(z\right)H_{\beta }\left(z\right)\frac{\partial {\sigma }_{\alpha z}^k}{\partial z}+\left(\frac{2H_{\beta }\left(z\right)}{R_{\alpha }}+\frac{H_{\alpha }\left(z\right)}{R_{\beta }}\right){\sigma }_{\alpha z}^k=0,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill & H_{\beta }\left(z\right)\frac{\partial {\sigma }_{\alpha \beta }^k}{\partial \alpha }+H_{\alpha }\left(z\right)\frac{\partial {\sigma }_{\beta \beta }^k}{\partial \beta }+H_{\alpha }\left(z\right)H_{\beta }\left(z\right)\frac{\partial {\sigma }_{\beta z}^k}{\partial z}+\left(\frac{2H_{\alpha }\left(z\right)}{R_{\beta }}+\frac{H_{\beta }\left(z\right)}{R_{\alpha }}\right){\sigma }_{\beta z}^k=0,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill & H_{\beta }\left(z\right)\frac{\partial {\sigma }_{\alpha z}^k}{\partial \alpha }+H_{\alpha }\left(z\right)\frac{\partial {\sigma }_{\beta z}^k}{\partial \beta }+H_{\alpha }\left(z\right)H_{\beta }\left(z\right)\frac{\partial {\sigma }_{zz}^k}{\partial z}-\frac{H_{\beta }\left(z\right)}{R_{\alpha }}{\sigma }_{\alpha \alpha }^k-\frac{H_{\alpha }\left(z\right)}{R_{\beta }}{\sigma }_{\beta \beta }^k+\left(\frac{H_{\beta }\left(z\right)}{R_{\alpha }}+\frac{H_{\alpha }\left(z\right)}{R_{\beta }}\right){\sigma }_{zz}^k=0,\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill & \frac{{\mathcal{D}}_1^k}{H_{\alpha }^2\left(z\right)}\frac{{\partial }^2{\mathcal{M}}^k}{\partial {\alpha }^2}+\frac{{\mathcal{D}}_2^k}{H_{\beta }^2\left(z\right)}\frac{{\partial }^2{\mathcal{M}}^k}{\partial {\beta }^2}+{\mathcal{D}}_3^k\frac{{\partial }^2{\mathcal{M}}^k}{\partial z^2}=0.\hfill \end{array}$$

(16)

Equations (13)–(16) are valid for all the structures proposed in Figure 1. Thanks to the orthogonal mixed curvilinear reference system, the formulation is valid for plates, cylinders, and cylindrical panels by assuming proper values for the two radii of curvature $R_{\alpha }$ and $R_{\beta }$. Equation (16) is the 3D Fick diffusion equation under steady-state conditions specifically written for spherical shells and already presented by Brischetto and Torre in [21]. In this fully coupled hygro-elastic problem, the 3D Fick diffusion relation is directly coupled with the 3D equilibrium equations. Equation (16) can be rewritten as:

$${\mathcal{D}}_1^{*k}\left(z\right)\frac{{\partial }^2\mathcal{M}}{\partial {\alpha }^2}+{\mathcal{D}}_2^{*k}\left(z\right)\frac{{\partial }^2\mathcal{M}}{\partial {\beta }^2}+{\mathcal{D}}_3^{*k}\frac{{\partial }^2\mathcal{M}}{\partial z^2}=0$$

(17)

in order to highlight the terms ${\mathcal{D}}_i^{*k}\left(z\right)$ (for each direction $\alpha$, $\beta$ and z) that takes into account the moisture diffusion coefficients ${\mathcal{D}}_i^k$ and the curvature terms $H_{\alpha }\left(z\right)$ and $H_{\beta }\left(z\right)$. It must be noted that ${\mathcal{D}}_i^{*k}\left(z\right)$ are dependent on z because of the curvature terms.

The closed-form solution of the coupled governing Equations (13)–(16) can be obtained by imposing the harmonic form for displacements and moisture content. These impositions directly affect the boundary conditions, because only simply supported configurations are considered. The harmonic forms for displacements and moisture content can be written as:

$$\begin{array}{cc}\hfill u^k(\alpha ,\beta ,z)& =U^k\left(z\right)\cos \left(\overline{\alpha }\alpha \right)\sin \left(\overline{\beta }\beta \right),\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill v^k(\alpha ,\beta ,z)& =V^k\left(z\right)\sin \left(\overline{\alpha }\alpha \right)\cos \left(\overline{\beta }\beta \right),\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill w^k(\alpha ,\beta ,z)& =W^k\left(z\right)\sin \left(\overline{\alpha }\alpha \right)sin\left(\overline{\beta }\beta \right),\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill {\mathcal{M}}^k(\alpha ,\beta ,z)& =M^k\left(z\right)\sin \left(\overline{\alpha }\alpha \right)\sin \left(\overline{\beta }\beta \right),\hfill \end{array}$$

(21)

where the two coefficients $\overline{\alpha }$ and $\overline{\beta }$ are defined as $\overline{\alpha }=\frac{m\pi }{a}$ and $\overline{\beta }=\frac{n\pi }{b}$: m and n are the half-wave numbers in the $\alpha$ and $\beta$ directions, a and b are the shell in-plane dimensions. $U^k\left(z\right)$, $V^k\left(z\right)$, and $W^k\left(z\right)$ are the displacement amplitudes. $M^k\left(z\right)$ is the moisture content amplitude. These amplitudes are the unknowns of the problem.

The introduction of the harmonic form for displacements and moisture content (Equations (18)–(21)), and geometrical and constitutive relations (Equations (1)–(12)) into the coupled hygro-elastic Equations (13)–(16) gives four differential relations written as functions of amplitudes for displacements and moisture content and related derivatives in z. The derivatives in $\alpha$ and $\beta$ directions are exactly computed thanks to the harmonic forms in Equations (18)–(21). Therefore, they are calculated algebraic values. The final system of equations is composed of four second-order differential equations in the thickness direction z. They do not have constant coefficients because the curvature terms $H_{\alpha }$ and $H_{\beta }$ depend on z. In order to transform this set of four second-order differential equations having non-constant coefficients into a set of differential equations with constant coefficients, each k physical layer must be divided into an appropriate number of fictitious layers. As a consequence, a further index j is introduced to count the number of fictitious layers: this index j varies from 1 to the total number of fictitious layers F. The total number of fictitious layers con be easily defined as $F=s·k$ where s indicates the number of subdivisions for each physical layer k. Assuming a proper thickness for each fictitious layer, the parametric coefficients $H_{\alpha }$ and $H_{\beta }$ can be calculated in its middle point and considered as constant. As such, a set of four second-order differential equations with constant coefficients is obtained. The compact form of this second-order system of partial differential equations in z is:

$$\begin{array}{cc}\hfill & A_1^j{U}^j+A_2^j{V}^j+A_3^j{W}^j+J_1^j{M}^j+A_4^j{U}_{,z}^j+A_5^j{W}_{,z}^j+A_6^j{U}_{,zz}^j=0,\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & A_7^j{U}^j+A_8^j{V}^j+A_9^j{W}^j+J_2^j{M}^j+A_{10}^j{V}_{,z}^j+A_{11}^j{W}_{,z}^j+A_{12}^j{V}_{,zz}^j=0,\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill & A_{13}^j{U}^j+A_{14}^j{V}^j+A_{15}^j{W}^j+J_3^j{M}^j+A_{16}^j{U}_{,z}^j+A_{17}^j{V}_{,z}^j+A_{18}^j{W}_{,z}^j+J_4^j{M}_{,z}^j+A_{19}^j{W}_{,zz}^j=0,\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill & (J_5^j+J_6^j)M^j+J_7^j{M}_{,zz}^j=0.\hfill \end{array}$$

(25)

The system of Equations (22)–(25) can be reduced to a system of first-order partial differential equations via the methodology described in [54,55]. The order of derivatives along the thickness direction can be reduced by redoubling the number of unknowns from 4 ($U^j$, $V^j$, $W^j$, $M^j$) to 8 ($U^j$, $V^j$, $W^j$, $M^j$, $U^{j^{{}^{\prime }}}$, $V^{j^{{}^{\prime }}}$, $W^{j^{{}^{\prime }}}$, $M^{j^{{}^{\prime }}}$) by including in the unknown vector the first partial derivatives in z (indicated by the superscript ${}^{\prime }$). Terms $M^j$ and $M^{j^{{}^{\prime }}}$ are calculated in analogy with the displacement amplitudes and related derivatives by solving the system of Equations (22)–(25). For this fully coupled case, the moisture content in the thickness direction must not be a priori defined or separately calculated by solving the Fick diffusion equation. The matrix form of the problem can be written as:

$$\begin{array}{cc}\hfill & \left[\begin{array}{cccccccc}{A}_6^j& 0& 0& 0& 0& 0& 0& 0\\ 0& A_{12}^j& 0& 0& 0& 0& 0& 0\\ 0& 0& A_{19}^j& 0& 0& 0& 0& 0\\ 0& 0& 0& J_7^j& 0& 0& 0& 0\\ 0& 0& 0& 0& A_6^j& 0& 0& 0\\ 0& 0& 0& 0& 0& A_{12}^j& 0& 0\\ 0& 0& 0& 0& 0& 0& A_{19}^j& 0\\ 0& 0& 0& 0& 0& 0& 0& J_7^j\end{array}\right]{\left[\begin{array}{c}{U}^j\\ V^j\\ W^j\\ M^j\\ U^{j^{{}^{\prime }}}\\ V^{j^{{}^{\prime }}}\\ W^{j^{{}^{\prime }}}\\ M^{j^{{}^{\prime }}}\end{array}\right]}^{{}^{\prime }}=\hfill \\ \hfill & \left[\begin{array}{cccccccc}0& 0& 0& 0& A_6^j& 0& 0& 0\\ 0& 0& 0& 0& 0& A_{12}^j& 0& 0\\ 0& 0& 0& 0& 0& 0& A_{19}^j& 0\\ 0& 0& 0& 0& 0& 0& 0& J_7^j\\ -A_1^j& -A_2^j& -A_3^j& -J_1^j& -A_4^j& 0& -A_5^j& 0\\ -A_7^j& -A_8^j& -A_9^j& -J_2^j& 0& -A_{10}^j& -A_{11}^j& 0\\ -A_{13}^j& -A_{14}^j& -A_{15}^j& -J_3^j& -A_{16}^j& -A_{17}^j& -A_{18}^j& -J_4^j\\ 0& 0& 0& -(J_5^j+J_6^j)& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{U}^j\\ V^j\\ W^j\\ M^j\\ U^{j^{{}^{\prime }}}\\ V^{j^{{}^{\prime }}}\\ W^{j^{{}^{\prime }}}\\ M^{j^{{}^{\prime }}}\end{array}\right].\hfill \end{array}$$

(26)

Equation (26) can be compacted as:

$${\mathit{D}}^j{{\mathit{X}}^j}^{{}^{\prime }}={\mathit{A}}^j{\mathit{X}}^j,$$

(27)

where ${\mathit{X}}^j={\left[U^j{V}^j{W}^j{M}^j{U}^{j^{{}^{\prime }}}{V}^{j^{{}^{\prime }}}{W}^{j^{{}^{\prime }}}{M}^{j^{{}^{\prime }}}\right]}^T$ and ${{\mathit{X}}^j}^{{}^{\prime }}=\frac{\partial {\mathit{X}}^j}{\partial z}$. Superscript T indicates a transposed vector. Equation (27) can be further written as:

$${{\mathit{X}}^j}^{{}^{\prime }}={\mathit{D}}^{j^{-1}}{\mathit{A}}^j{\mathit{X}}^j\Rightarrow {{\mathit{X}}^j}^{{}^{\prime }}={\mathit{A}}^{{*}^j}{\mathit{X}}^j,$$

(28)

where ${\mathit{A}}^{{*}^j}={\mathit{D}}^{j^{-1}}{\mathit{A}}^j$. The solution of the problem presented in Equation (28), which is a general homogeneous system of first-order differential equations, can be written as:

$${\mathit{X}}^j\left({\tilde{z}}^j\right)={\mathrm{e}}^{\left({\mathit{A}}^{{*}^j}{\tilde{z}}^j\right)}{\mathit{X}}^j\left(0\right),$$

(29)

$${\mathit{X}}^j\left(h^j\right)={\mathit{A}}^{{**}^j}{\mathit{X}}^j\left(0\right)\Rightarrow {\mathit{X}}_t^j={\mathit{A}}^{{**}^j}{\mathit{X}}_b^j,$$

(30)

where ${\mathit{A}}^{{**}^j}$ is the exponential matrix. Equation (30) is employed to define the unknown vector ${\mathit{X}}^j\left(h^j\right)$ that represents the unknown vector at the top (t) of the j-th layer, when the exponential matrix ${\mathit{A}}^{{**}^j}={\mathrm{e}}^{({\mathit{A}}^{{*}^j}{h}^j)}$ has been computed in each j fictitious layer considering $h^j$ as thickness. ${\mathit{X}}^j\left(h^j\right)$ is ${\mathit{X}}_t^j$ and ${\mathit{X}}^j\left(0\right)$ is ${\mathit{X}}_b^j$. t and b indicate the top and the bottom of the j-th fictitious layer, respectively. The exponential matrix can be expanded as a power series and calculated for each j fictitious layer having a thickness $h^j$ as:

$${\mathit{A}}^{{**}^j}={\mathrm{e}}^{({\mathit{A}}^{{*}^j}{h}^j)}=\mathit{I}+{\mathit{A}}^{{*}^j}{h}^j+\frac{{{\mathit{A}}^{{*}^j}}^2}{2!}{h^j}^2+\frac{{{\mathit{A}}^{{*}^j}}^3}{3!}{h^j}^3+\cdots +\frac{{{\mathit{A}}^{{*}^j}}^N}{N!}{h^j}^N,$$

(31)

where the identity matrix $\mathit{I}$ has dimensions of $8\times 8$.

Equation (30) links the top and bottom unknown vectors (primary variables and the related derivatives in z) within the same mathematical layer j. To move from the j-th fictitious layer to the next one, a set of interlaminar continuity conditions must be imposed. These conditions involve the displacements, the transverse shear/normal stresses, the moisture content, and the transverse normal moisture flux at each fictitious layer interface. The continuity of the displacement components and moisture content at each interface can be imposed as:

$$u_b^{j+1}=u_t^j,v_b^{j+1}=v_t^j,w_b^{j+1}=w_t^j,{\mathcal{M}}_b^{j+1}={\mathcal{M}}_t^j.$$

(32)

The conditions given in Equation (32) can be easily imposed for the related amplitudes $U^j$, $V^j$, $W^j$, and $M^j$. In the same way, the continuity of the transverse shear stresses, transverse normal stress and transverse normal moisture flux in the z direction can be written as:

$${\sigma }_{z{z}_b}^{j+1}={\sigma }_{z{z}_t}^j,{\sigma }_{\alpha z_b}^{j+1}={\sigma }_{\alpha z_t}^j,{\sigma }_{\beta z_b}^{j+1}={\sigma }_{\beta z_t}^j,g_{z_b}^{j+1}=g_{z_t}^j.$$

(33)

In Equations (32) and (33), for each considered variable, the continuity is imposed between the value at the top (t) of the general j-th layer and the value at the bottom (b) of the ($j+1$)-th layer. The use of constitutive equations and harmonic forms allows us to obtain an amplitude form of Equations (32) and (33). The methodology has already been discussed in [47,48,49], however, the only difference lies in the dimension of the involved matrices and vectors. The amplitude form of Equations (32) and (33) may be rewritten in a matrix form introducing the transfer matrix:

$${\left[\begin{array}{c}U\\ V\\ W\\ M\\ U^{{}^{\prime }}\\ V^{{}^{\prime }}\\ W^{{}^{\prime }}\\ M^{{}^{\prime }}\end{array}\right]}_b^{j+1}={\left[\begin{array}{cccccccc}1& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\\ T_1& 0& T_2& 0& T_3& 0& 0& 0\\ 0& T_4& T_5& 0& 0& T_6& 0& 0\\ T_7& T_8& T_9& {\tau }_1& 0& 0& T_{10}& 0\\ 0& 0& 0& 0& 0& 0& 0& {\tau }_2\end{array}\right]}^{j+1,j}{\left[\begin{array}{c}U\\ V\\ W\\ M\\ U^{{}^{\prime }}\\ V^{{}^{\prime }}\\ W^{{}^{\prime }}\\ M{}^{\prime }\end{array}\right]}_t^j\Rightarrow {\mathit{X}}_b^{j+1}={\mathit{T}}^{j+1,j}{\mathit{X}}_t^j.$$

(34)

where ${\mathit{T}}^{j+1,j}$ is the transfer matrix mentioned previously. The diagonal part composed of 1 indicates the congruence conditions in terms of displacement components and the continuity of moisture content written in Equation (32). The other coefficients $T_i$ and ${\tau }_i$ indicate the stress and moisture flux continuity conditions of Equation (33) in terms of displacements and moisture content (and related derivatives in z), respectively. Equation (34) allows to link the unknown vector defined at the top of the j-th layer with the unknown vector defined at the bottom of the ($j+1$)-th layer.

All configurations consider simply supported sides: these boundary conditions are automatically satisfied via the proper harmonic forms in Equations (18)–(21) used for all primary variables. Thus, it is possible to write:

$$\begin{array}{cc}\hfill & M=0,w=v=0\mathrm{for}\alpha =0, a,\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill & M=0,w=u=0\mathrm{for}\beta =0, b.\hfill \end{array}$$

(36)

In addition, load boundary conditions must be imposed at the outer faces of the structures. These conditions can be written as:

$$\begin{array}{cc}\hfill & {\sigma }_{zz}=0,{\sigma }_{\alpha z}=0,{\sigma }_{\beta z}=0,M=M_{ext}\mathrm{for}z=\pm h/2,\hfill \end{array}$$

(37)

Equation (37) can be rewritten, in matrix form, as:

$$\left[\begin{array}{cccccccc}-\frac{C_{13}^M}{H_{{\alpha }_t}^M}\overline{\alpha }& -\frac{C_{23}^M}{H_{{\beta }_t}^M}\overline{\beta }& \frac{C_{13}^M}{H_{{\alpha }_t}^M{R}_{\alpha }}+\frac{C_{23}^M}{H_{{\beta }_t}^M{R}_{\beta }}& -{\xi }_z^M& 0& 0& C_{33}^M& 0\\ 0& -\frac{C_{44}^M}{H_{{\beta }_t}^M{R}_{\beta }}& \frac{C_{44}^M}{H_{{\beta }_t}^M}\overline{\beta }& 0& 0& C_{44}^M& 0& 0\\ -\frac{C_{55}^M}{H_{{\alpha }_t}^M{R}_{\alpha }}& 0& \frac{C_{55}^M}{H_{{\alpha }_t}^M}\overline{\alpha }& 0& C_{55}^M& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\end{array}\right]{\left[\begin{array}{c}U\\ V\\ W\\ M\\ U^{{}^{\prime }}\\ V^{{}^{\prime }}\\ W^{{}^{\prime }}\\ M^{{}^{\prime }}\end{array}\right]}_t^F=\left[\begin{array}{c}0\\ 0\\ 0\\ M_t\end{array}\right],$$

(38)

$$\left[\begin{array}{cccccccc}-\frac{C_{13}^1}{H_{{\alpha }_b}^1}\overline{\alpha }& -\frac{C_{23}^1}{H_{{\beta }_b}^1}\overline{\beta }& \frac{C_{13}^1}{H_{{\alpha }_b}^1{R}_{\alpha }}+\frac{C_{23}^1}{H_{{\beta }_b}^1{R}_{\beta }}& -{\xi }_z^1& 0& 0& C_{33}^1& 0\\ 0& -\frac{C_{44}^1}{H_{{\beta }_b}^1{R}_{\beta }}& \frac{C_{44}^1}{H_{{\beta }_b}^1}\overline{\beta }& 0& 0& C_{44}^1& 0& 0\\ -\frac{C_{55}^1}{H_{{\alpha }_b}^1{R}_{\alpha }}& 0& \frac{C_{55}^1}{H_{{\alpha }_b}^1}\overline{\alpha }& 0& C_{55}^1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\end{array}\right]{\left[\begin{array}{c}U\\ V\\ W\\ M\\ U^{{}^{\prime }}\\ V^{{}^{\prime }}\\ W^{{}^{\prime }}\\ M^{{}^{\prime }}\end{array}\right]}_b^1=\left[\begin{array}{c}0\\ 0\\ 0\\ M_b\end{array}\right].$$

(39)

Equations (38) and (39) can be compacted as:

$$\begin{array}{cc}\hfill & {\mathit{B}}_t^F{\mathit{X}}_t^F={\mathcal{P}}_{\mathit{t}},\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill & {\mathit{B}}_b^1{\mathit{X}}_b^1={\mathcal{P}}_{\mathit{b}}.\hfill \end{array}$$

(41)

Superscript F indicates the last fictitious layer and superscript 1 shows the first fictitious layer. Vectors ${\mathcal{P}}_{\mathit{b}}$ and ${\mathcal{P}}_{\mathit{t}}$ include the impositions related to the load conditions in the three directions $\alpha$, $\beta$, and z, as well as the moisture content. Assuming a classical hygro-elastic stress analysis, the mechanical loads in the $\alpha$, $\beta$, and z directions are set equal to zero ($P_{{\alpha }_t}=P_{{\beta }_t}=P_{z_t}=P_{{\alpha }_b}=P_{{\beta }_b}=P_{z_b}=0$).

Considering Equations (40) and (41), it is possible to link the bottom of the first fictitious layer to the top of the F-th fictitious layer by recursively introducing Equation (34) into Equation (30). It can be obtained as:

$$\begin{array}{cc}\hfill {\mathit{X}}_t^F=& \left({{\mathit{A}}^{**}}^F{\mathit{T}}^{F,F-1}{{\mathit{A}}^{**}}^{F-1}{\mathit{T}}^{F-1,F-2}\dots \dots {{\mathit{A}}^{**}}^2{\mathit{T}}^{2,1}{{\mathit{A}}^{**}}^1\right){\mathit{X}}_b^1.\hfill \end{array}$$

(42)

Equation (42) defines the $8\times 8$ matrix ${\mathit{H}}_m$ for multilayered structures. This matrix is different from the ${\mathit{H}}_m$ matrix presented in [47,48,49] for the pure elastic analysis; the difference is due to the coupling of the 3D Fick diffusion relation with the 3D equilibrium equations. In a compact form, it can be rewritten as:

$${\mathit{X}}_t^F={\mathit{H}}_m{\mathit{X}}_b^1.$$

(43)

Using Equation (43), Equation (40) can be rewritten in terms of ${\mathit{X}}_b^1$:

$${\mathit{B}}_t^F{\mathit{H}}_m{\mathit{X}}_b^1={\mathcal{P}}_{\mathit{t}}.$$

(44)

Equations (44) and (41) can be now compacted as:

$$\mathit{E}{\mathit{X}}_b^1=\mathcal{P},$$

(45)

where

$$\mathit{E}=\left[\begin{array}{c}{\mathit{B}}_t^F{\mathit{H}}_m\\ {\mathit{B}}_b^1\end{array}\right],\mathcal{P}=\left[\begin{array}{c}{\mathcal{P}}_{\mathit{t}}\\ {\mathcal{P}}_{\mathit{b}}\end{array}\right].$$

(46)

The matrix $\mathit{E}$ is always $8\times 8$ even if the number F of fictitious layers changes and a layer-wise approach is employed. $\mathcal{P}$ is the vector containing all the hygro-elastic loading impositions. The matrix $\mathit{E}$ changes in size with respect to the pure elastic analysis in [47]. Considering the addition of the 3D Fick diffusion equation in the 3D equilibrium equations for shells, this new model can be defined as the generalization of the pure mechanical model proposed in [47] by Brischetto. The system in Equation (45) is formally the same as that seen in [47,48,49] for the pure elastic investigation, but the coupling of the 3D Fick diffusion equation with the 3D equilibrium equations is now considered. The vector of unknowns has dimensions of ${\mathit{X}}_b^1$ has $8\times 1$.

After the calculation of the unknowns at the bottom, Equations (30) and (34) can be subsequently substituted to define the displacements and moisture content (and the related derivatives with respect to z) through the global thickness of the shell or the plate. Therefore, the unknown trends along the z direction are obtained.

The present formulation can be implemented in a Matlab code permitting the evaluation of stresses, strains, and displacements along the thickness direction z for all the proposed geometries with different materials, thickness ratios, lamination schemes, and load conditions.

4. Results

The first part of this section proposes two validation results used to certify the correctness of the developed general 3D exact coupled hygro-elastic shell model. Comparisons with [21] are presented to define the proper order of expansion N for the exponential matrix in Equation (31) and the appropriate number of fictitious layers F for the computation of the constant curvature terms of the shell. After these two validation cases, four new benchmarks are proposed. N and F values are opportunely identified and the new results aim to discuss the effects of thickness ratios, geometry configurations, lamination sequences, embedded materials, and moisture content impositions when a 3D hygro-elastic analysis of plates and shells is executed. Geometrical data for the validation results and benchmarks are collected in

Table 1. Moisture contents, geometrical data, and lamination sequences, as already proposed in [21], for the preliminary validation results.

	Assessment 1	Assessment 2
	Rectangular Plate	Cylindrical Shell
a (m)	4	$\frac{\pi }{3}$R${}_{\alpha }$
b (m)	12	50
h (m)	2, 1, 0.4, 0.2, 0.08, 0.04	5, 2.5, 1, 0.2, 0.1, 0.02
R${}_{\alpha }$ (m)	∞	10
R${}_{\beta }$ (m)	∞	∞
h${}_1$	h/3	h
h${}_2$	h/3	-
h${}_3$	h/3	-
Lamination	0°/90°/0°	0°
${\mathcal{M}}_t$ (%)	+0.5	+0.5
${\mathcal{M}}_b$ (%)	+0.1	+0.1
m	1	1
n	1	1

and

Table 2. Moisture contents, geometrical data, and lamination sequences for the four benchmarking analyses. t can assume the values 0.5, 0.25, 0.1, 0.02, and 0.01. s can assume the values 2.5, 1, 0.2, 0.1, and 0.02.

	B1	B2	B3	B4
	Square Plate	Cylindrical Shell	Cylinder	Spherical Shell
a (m)	1	$\frac{\pi }{3}{R}_{\alpha }$	2$\pi$R${}_{\alpha }$	$\frac{\pi }{3}{R}_{\alpha }$
b (m)	1	30	30	$\frac{\pi }{3}{R}_{\alpha }$
h (m)	t	s	s	s
R${}_{\alpha }$ (m)	∞	10	10	10
R${}_{\beta }$ (m)	∞	∞	∞	10
h${}_1$ (m)	h/10	h/3	h/4	h/3
h${}_2$ (m)	h/10	h/3	h/4	h/3
h${}_3$ (m)	6h/10	h/3	h/4	h/3
h${}_4$ (m)	h/10	-	h/4	-
h${}_5$ (m)	h/10	-	-	-
Lamination	0°/90°/PVC/90°/0°	0°/90°/0°	0°/90°/0°/90°	0°/90°/0°
${\mathcal{M}}_t$ (%)	+1.0	+1.0	+1.0	+1.0
${\mathcal{M}}_b$ (%)	0.0	+0.5	+0.5	+0.5
$P_{z_t}$ (Pa)	10,000	0	0	0
m	1	1	2	1
n	1	1	1	1

. The material data shown in

Table 3. Mechanical and hygrometrical properties of the materials used in preliminary results and new benchmarks [21].

	Composite One	Composite Two	PVC
$E_1$ (GPa)	172	138	3
$E_2$ (GPa)	6.9	8.5	3
$E_3$ (GPa)	6.9	8.5	3
${\nu }_{12}$	0.25	0.29	0.4
${\nu }_{13}$	0.25	0.29	0.4
${\nu }_{23}$	0.25	0.36	0.4
$G_{12}$ (GPa)	3.4	4.5	$\frac{E}{2(1+\nu )}$
$G_{13}$ (GPa)	3.4	4.5	$\frac{E}{2(1+\nu )}$
$G_{23}$ (GPa)	1.4	3.2	$\frac{E}{2(1+\nu )}$
${\eta }_1$$(\frac{1}{\%})$	0	0	$0.28\times {10}^{-2}$
${\eta }_2$$(\frac{1}{\%})$	$0.4\times {10}^{-2}$	$0.4\times {10}^{-2}$	$0.28\times {10}^{-2}$
${\eta }_3$$(\frac{1}{\%})$	$0.4\times {10}^{-2}$	$0.4\times {10}^{-2}$	$0.28\times {10}^{-2}$
${\mathcal{D}}_1$ (kg/ms)	7.04	7.04	$9.324\times {10}^{-8}$
${\mathcal{D}}_2$ (kg/ms)	4.96	4.96	$9.324\times {10}^{-8}$
${\mathcal{D}}_3$ (kg/ms)	4.96	4.96	$9.324\times {10}^{-8}$

are the Young moduli, Poisson ratios, shear moduli, moisture expansion coefficients, and moisture diffusion coefficients in the three directions of the material reference system.

The acronyms employed represent the different models involved in the performed analyses: 3D-u-$\mathcal{M}$ means the 3D fully coupled hygro-elastic model where the displacements (indicated as u) and the moisture content profile (indicated as $\mathcal{M}$) are fully coupled. The displacements u, v, and w along the three directions and the moisture content profile $\mathcal{M}$ are the primary variables of the hygro-elastic problem. The general acronym 3D() is used for 3D uncoupled hygro-elastic models presented by Brischetto and Torre in [21]. The variable separately solved and the version of the Fick diffusion equation involved are indicated inside the parentheses. Specifically, the 3D(${\mathcal{M}}_c$, 3D) indicates the 3D uncoupled model where the 3D Fick diffusion equation is separately solved, the 3D(${\mathcal{M}}_c$, 1D) denotes the 3D uncoupled model where the 1D Fick diffusion equation is separately solved and 3D(${\mathcal{M}}_a$) denotes the 3D uncoupled model including the a priori assumed linear moisture content profile.

4.1. Comparisons for Validations

Two validation cases are proposed for this new three-dimensional general exact coupled hygro-elastic shell model (denoted by 3D-u-$\mathcal{M}$). A rectangular plate and a cylindrical shell panel are analyzed for different thickness ratios, moisture content profiles, and lamination schemes. When $N=3$ and $M=300$, the 3D-u-$\mathcal{M}$ model gives the same results as the 3D(${\mathcal{M}}_c$, 3D) model presented by Brischetto and Torre in [21]. The convergence of the 3D-u-$\mathcal{M}$ model (obtained thanks to the coupled formulation) occurs for fewer mathematical layers than those employed in the 3D(${\mathcal{M}}_c$, 3D) model.

The first preliminary validation results take into account a three-layered composite rectangular plate having simply supported sides. The geometrical data can be seen in the first column of Table 1. Material data employed for this case are collected in Table 3, in correspondence with the Composite One column. The results used as a comparison are those presented by Brischetto and Torre in [21]. The model used as a reference is the 3D(${\mathcal{M}}_c$, 3D) model. The new 3D-u-$\mathcal{M}$ model couples the 3D Fick diffusion relation with the 3D equilibrium relations. The moisture content profile is an unknown due to the displacements. The proposed new solution employs $N=3$ and $F=300$. This preliminary analysis uses these two values in order to understand whether the 3D-u-$\mathcal{M}$ model gives the same results obtained with the 3D(${\mathcal{M}}_c$, 3D) model.

Table 4. First preliminary validation results, $0^{\circ }/{90}^{\circ }/0^{\circ }$ rectangular composite plate with moisture content ${\mathcal{M}}_t=+0.5\%$, and ${\mathcal{M}}_b=+0.1\%$ for m = n = 1. The solution for comparisons is the 3D (${\mathcal{M}}_c$, 3D) model by Brischetto and Torre [21] employing the 3D calculated moisture content profile.

$\mathit{a}/\mathit{h}$	2	4	10	20	50	100
w (mm) in $\overline{z}=+h/2$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_c$, 3D) [21]	−0.0710	0.0133	0.0977	0.204	0.513	1.03
3D-u-$\mathcal{M}$	−0.0710	0.0133	0.0977	0.204	0.513	1.03
${\sigma }_{\alpha \alpha }$ (MPa) in $\overline{z}=0$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_c$, 3D) [21]	11.9	7.05	3.67	3.04	2.86	2.83
3D-u-$\mathcal{M}$	11.9	7.05	3.67	3.04	2.86	2.83
${\sigma }_{\beta \beta }$ (MPa) in $\overline{z}=+h$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_c$, 3D) [21]	−12.5	−13.0	−13.2	−13.3	−13.3	−13.3
3D-u-$\mathcal{M}$	−12.5	−13.0	−13.2	−13.3	−13.3	−13.3

proposes the transverse displacement and in-plane stress components for different thickness ratios $a/h$. The present 3D-u-$\mathcal{M}$ model is perfectly in accordance with the reference solution for each thickness ratio because both the material and thickness layer effects are evaluated thanks the use of the 3D Fick moisture diffusion Equation (fully coupled or separately solved).

The second preliminary validation results show a simply supported single-layered composite cylindrical shell panel. The geometrical data can be seen in the proper column of Table 1 and the material data can be seen in Table 3 (Composite One column). The employed reference results derive from the 3D (${\mathcal{M}}_c$, 3D) model presented in [21]. The 3D-u-$\mathcal{M}$ model adds the 3D Fick diffusion equation to the 3D equilibrium equations and it has all the peculiarities previously discussed in the first validation results.

Table 5. Second preliminary validation results, one-layered cylindrical composite shell with moisture content ${\mathcal{M}}_t=+0.5\%$, and ${\mathcal{M}}_b=+0.1\%$ for m = n = 1. The solution for comparisons is the 3D (${\mathcal{M}}_c$, 3D) model by Brischetto and Torre [21] employing the 3D calculated moisture content profile.

${\mathit{R}}_{\mathit{\alpha }}/\mathit{h}$	2	4	10	50	100	500
w (m) in $\overline{z}=0$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_c$, 3D) [21]	−0.00271	−0.00245	−0.00152	0.00823	0.0326	0.212
3D-u-$\mathcal{M}$	−0.00271	−0.00245	−0.00152	0.00823	0.0326	0.212
${\sigma }_{\alpha \alpha }$ (MPa) in $\overline{z}=0$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_c$, 3D) [21]	16.7	9.32	0.654	−12.3	−22.4	−28.2
3D-u-$\mathcal{M}$	16.7	9.32	0.654	−12.3	−22.4	−28.2
${\sigma }_{\alpha \alpha }$ (MPa) in $\overline{z}=+h$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_c$, 3D) [21]	23.8	14.3	7.03	12.6	22.2	28.1
3D-u-$\mathcal{M}$	23.8	14.3	7.03	12.6	22.2	28.1

indicates the same quantities evaluated in the previous case for different thickness ratios $R_{\alpha }/h$. The 3D-u-$\mathcal{M}$ model is coincident for each thickness ratio $R_{\alpha }/h$ with the reference solution [21] because both use the 3D Fick diffusion relation: coupled with the 3D equilibrium equations or separately solved in an external tool.

In the next section, four new benchmarking analyses are proposed where different geometrical data, thickness ratios, lamination stacking sequences, materials, and moisture content profiles are presented. For all the four benchmarks, a number of $F=300$ fictitious layers and $N=3$ order of expansion for the exponential matrix are employed.

4.2. New Benchmarks

Four benchmarking analyses are presented in this subsection. The benchmarks consider plates, cylinders, cylindrical panels, and spherical panels (see Figure 1). Outer moisture content impositions are applied directly to the structures; different amplitudes and half-wave numbers are taken into account. All geometrical and material data for the four benchmarks are presented in Table 2 and Table 3, respectively. In these new results, different materials, lamination schemes, and load configurations are considered. The 3D-u-$\mathcal{M}$ model will be used for these new analyses. The $N=3$ order of expansion for the exponential matrix and $F=300$ fictitious layers will be used for each benchmark. The 3D coupled hygro-elastic shell model will be extensively discussed and compared with past uncoupled 3D results where 3D calculated moisture content profiles (${\mathcal{M}}_c$, 3D), 1D calculated moisture content profiles (${\mathcal{M}}_c$, 1D), and assumed moisture content profiles (${\mathcal{M}}_a$) were externally defined. These new results may be considered as test cases for those scientists involved in the development of two-dimensional and/or three-dimensional analytical or numerical shell models for hygro-elastic stress analysis of multilayered structures.

The first benchmark shows a square sandwich composite plate with simply supported sides (see Figure 1). The geometrical and lamination data used for this benchmark are included in the Table 2 (column B1). The considered thickness ratios $a/h$ are also highlighted in Table 2. The five-layered configuration is made of two skins (the bottom one $0^{\circ }/{90}^{\circ }$ and the top one ${90}^{\circ }/0^{\circ }$) and a PVC core. The elastic and hygroscopic properties of the skins and the core are collected in Table 3 in the columns Composite Two and PVC, respectively. Figure 2 shows the moisture content profile through the z direction of a thick ($a/h=2$) and a thin ($a/h=100$) square plate. For the thin case, the moisture content profiles (calculated via (${\mathcal{M}}_c$, 3D), via (${\mathcal{M}}_c$, 1D), and obtained from (u, $\mathcal{M}$) as a primary variable) are the same. The linear assumed profile (${\mathcal{M}}_a$) is wrong because it did not consider the different hygroscopic properties of the skins and core in the sandwich configuration. For thick configurations, the u-$\mathcal{M}$ profile is in accordance with the (${\mathcal{M}}_c$, 3D) profile; these two models display behaviors through the z direction with respect to the (${\mathcal{M}}_c$, 1D) profile and the (${\mathcal{M}}_a$) profile because the two previous profiles take into account the thickness and material layer effects.

Table 6. Benchmark number one, sandwich $0^{\circ }/{90}^{\circ }/PVC/{90}^{\circ }/0^{\circ }$ square plate with moisture content ${\mathcal{M}}_t=+1.0\%$, and ${\mathcal{M}}_b=0\%$ and a mechanical load ${P_z}_t=$ 10,000 Pa for $m=n=1$. Three-dimensional uncoupled hygro-elastic models from [21] with the moisture content profile were assumed to be linear (3D (${\mathcal{M}}_a$)), defined via the one-dimensional Fick diffusion relation (3D (${\mathcal{M}}_c$, 1D)) and via the three-dimensional Fick diffusion relation (3D (${\mathcal{M}}_c$, 3D)).

a/h	2	4	10	50	100
$\mathcal{M}$(%) in $\tilde{z}=+4h/5$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_a$) [21]	0.8000	0.8000	0.8000	0.8000	0.8000
3D (${\mathcal{M}}_c$, 1D) [21]	1.0000	1.0000	1.0000	1.0000	1.0000
3D (${\mathcal{M}}_c$, 3D) [21]	0.8914	0.9709	0.9952	0.9998	1.0000
3D-u-$\mathcal{M}$	0.8914	0.9709	0.9952	0.9998	1.0000
v (${10}^{-2}$ mm) in $\tilde{z}=+4h/5$, $\alpha =a/2$ and $\beta =0$
3D (${\mathcal{M}}_a$) [21]	−7.165	−9.489	−10.12	−12.11	−18.49
3D (${\mathcal{M}}_c$, 1D) [21]	−7.759	−10.18	−10.90	−12.94	−19.31
3D (${\mathcal{M}}_c$, 3D) [21]	−7.017	−9.879	−10.84	−12.93	−19.31
3D-u-$\mathcal{M}$	−7.017	−9.879	−10.84	−12.93	−19.31
w (${10}^{-1}$ mm) in $\tilde{z}=+h/2$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_a$) [21]	−0.3691	0.7748	2.976	26.55	120.8
3D (${\mathcal{M}}_c$, 1D) [21]	−1.043	0.6929	3.678	30.89	129.5
3D (${\mathcal{M}}_c$, 3D) [21]	−0.8156	0.6990	3.667	30.88	129.5
3D-u-$\mathcal{M}$	−0.8156	0.6990	3.667	30.88	129.5
${\sigma }_{zz}$ (${10}^4$ Pa) in $\tilde{z}=+h/5$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_a$) [21]	−95.33	−19.54	−2.631	0.02410	0.1033
3D (${\mathcal{M}}_c$, 1D) [21]	−76.76	−15.07	−1.941	0.05121	0.1101
3D (${\mathcal{M}}_c$, 3D) [21]	−55.91	−13.73	−1.908	0.05127	0.1101
3D-u-$\mathcal{M}$	−55.91	−13.73	−1.908	0.05127	0.1101
${\sigma }_{\alpha \alpha }$(MPa) in $\tilde{z}=+h$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_a$) [21]	61.71	47.41	43.69	58.53	105.7
3D (${\mathcal{M}}_c$, 1D) [21]	69.10	54.11	49.97	64.64	111.8
3D (${\mathcal{M}}_c$, 3D) [21]	58.79	51.83	49.64	64.62	111.8
3D-u-$\mathcal{M}$	58.79	51.83	49.64	64.62	111.8
${\gamma }_{\beta z}$(${10}^{-4}$) in $\tilde{z}=+h/3$, $\alpha =a/2$ and $\beta =0$
3D (${\mathcal{M}}_a$) [21]	−15.18	−7.892	−3.402	0.06178	1.279
3D(${\mathcal{M}}_c$, 1D) [21]	−16.93	−8.657	−3.628	0.02663	1.262
3D (${\mathcal{M}}_c$, 3D) [21]	−12.67	−8.024	−3.586	0.02696	1.262
3D-u-$\mathcal{M}$	−12.67	−8.024	−3.586	0.02696	1.262

shows the moisture content, the in-plane and transverse displacements, the transverse normal and in-plane normal stresses, and the transverse shear strains. Different thickness coordinate positions and different $a/h$ thickness ratios are considered. For thin structures, the 3D (${\mathcal{M}}_a$) model gives the wrong results because it discards both material and thickness layer effects. For thick or moderately thick configuration, the 3D (${\mathcal{M}}_a$) model remains inadequate. In this last case, the 3D (${\mathcal{M}}_c$, 1D), 3D (${\mathcal{M}}_c$, 3D), and 3D-u-$\mathcal{M}$ models provide the same results because they are able to include the material layer effect while the thickness layer effect is negligible. Figure 3 gives the moisture content profile, the in-plane and transverse displacement components, and three stress components through the thickness direction of a moderately thick structure ($a/h=4$). The two displacement components are nonlinear and not constant because of the thickness value considered and the sandwich configuration involved. The transverse normal stress component is continuous and it fulfills the boundary load conditions imposed at the outer surfaces (${\sigma }_{{zz}_t}=P_{z_t}=$10,000 Pa and ${\sigma }_{{zz}_b}=P_{z_b}=0$). The in-plane normal and in-plane shear stress components are discontinuous at each physical layer interface of the sandwich. The different values depend on the mechanical and hygrometric properties of the materials embedded in the sandwich configuration.

The second benchmark shows a composite cylindrical shell with simply supported sides, and equally divided into three layers ($0^{\circ }/{90}^{\circ }/0^{\circ }$ lamination sequence); the structure is shown in Figure 1. Geometrical data are available in Table 2 and material data are available in the column Composite One of Table 3. The moisture content profiles in the z direction (Figure 4) show two different $R_{\alpha }/h$ ratios: in the case of the thick cylindrical shell, the (${\mathcal{M}}_c$, 3D), and the (u-$\mathcal{M}$) solutions display the proper moisture content profile because they consider the thickness layer effect (the profile is not linear inside each single thick layer). The (${\mathcal{M}}_a$) profile is always rectilinear for the entire multilayered cylindrical shell for both thick and thin configurations. In the case of a thin cylindrical shell panel, all the moisture content profiles are correct because neither material nor thickness layer effects are involved. In

Table 7. Benchmark number two, $0^{\circ }/{90}^{\circ }/0^{\circ }$ composite cylindrical shell with moisture content ${\mathcal{M}}_t=+1\%$, and ${\mathcal{M}}_b=+0.5\%$ for $m=n=1$. Three-dimensional uncoupled hygro-elastic models from [21] with the moisture content profile assumed to be linear (3D (${\mathcal{M}}_a$)), defined via the one-dimensional Fick diffusion relation (3D (${\mathcal{M}}_c$, 1D)) and via the three-dimensional Fick diffusion relation (3D (${\mathcal{M}}_c$, 3D)).

${\mathit{R}}_{\mathit{\alpha }}/\mathit{h}$	4	10	50	100	500
$\mathcal{M}$(%) in $\tilde{z}=+h/2$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_a$) [21]	0.7500	0.7500	0.7500	0.7500	0.7500
3D (${\mathcal{M}}_c$, 1D) [21]	0.7500	0.7500	0.7500	0.7500	0.7500
3D (${\mathcal{M}}_c$, 3D) [21]	0.6850	0.7389	0.7495	0.7499	0.7500
3D-u-$\mathcal{M}$	0.6850	0.7389	0.7495	0.7499	0.7500
u (mm) in $\tilde{z}=+h/3$, $\alpha =0$ and $\beta =b/2$
3D (${\mathcal{M}}_a$) [21]	−1.643	−1.243	2.555	5.167	6.794
3D (${\mathcal{M}}_c$, 1D) [21]	−1.643	−1.243	2.555	5.167	6.794
3D (${\mathcal{M}}_c$, 3D) [21]	−1.514	−1.224	2.554	5.167	6.794
3D-u-$\mathcal{M}$	−1.514	−1.224	2.554	5.167	6.794
w (mm) in $\tilde{z}=+h/2$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_a$) [21]	−3.061	−2.621	8.252	16.11	21.17
3D (${\mathcal{M}}_c$, 1D) [21]	−3.061	−2.621	8.252	16.11	21.17
3D (${\mathcal{M}}_c$, 3D) [21]	−2.838	−2.582	8.249	16.11	21.17
3D-u-$\mathcal{M}$	−2.838	−2.582	8.249	16.11	21.17
${\sigma }_{zz}$ (MPa) in $\tilde{z}=+2h/3$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_a$) [21]	−1.314	−0.5024	−0.1344	−0.06058	−0.007403
3D (${\mathcal{M}}_c$, 1D) [21]	−1.314	−0.5024	−0.1344	−0.06058	−0.007403
3D (${\mathcal{M}}_c$, 3D) [21]	−1.207	−0.4957	−0.1343	−0.06057	−0.007403
3D-u-$\mathcal{M}$	−1.207	−0.4957	−0.1343	−0.06057	−0.007403
${\sigma }_{\beta \beta }$(MPa) in $\tilde{z}=+h$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_a$) [21]	−25.79	−26.27	−25.90	−25.73	−25.75
3D (${\mathcal{M}}_c$, 1D) [21]	−25.79	−26.27	−25.90	−25.73	−25.75
3D (${\mathcal{M}}_c$, 3D) [21]	−25.86	−26.28	−25.90	−25.73	−25.75
3D-u-$\mathcal{M}$	−25.86	−26.28	−25.90	−25.73	−25.75
${\gamma }_{\alpha \beta }$(${10}^{-3}$) in $\tilde{z}=0$, $\alpha =0$ and $\beta =0$
3D (${\mathcal{M}}_a$) [21]	−1.026	−0.7063	−0.2706	−0.07889	0.01048
3D (${\mathcal{M}}_c$, 1D) [21]	−1.026	−0.7063	−0.2706	−0.07889	0.01048
3D (${\mathcal{M}}_c$, 3D) [21]	−0.9645	−0.6992	−0.2705	−0.07888	0.01048
3D-u-$\mathcal{M}$	−0.9645	−0.6992	−0.2705	−0.07888	0.01048

, the results in terms of displacements and stresses are presented. For thin structures, the four three-dimensional models always provide the same values. For thicker plates, 3D (${\mathcal{M}}_c$, 3D) and 3D-u-$\mathcal{M}$ models are coincident and show bigger differences with respect to the other two models because they are able to consider the thickness layer effects. Figure 5 gives the displacement and stress components for the moderately thick ($R_{\alpha }/h=10$) cylindrical shell panel. The typical peculiarities of a multilayered composite cylindrical panel are shown: a zigzag form of displacements, a moisture content profile, and stresses. It is possible to notice the proper inclusion in the model of the interlaminar continuity conditions for displacements and transverse normal/shear stresses. In-plane stresses and strains are discontinuous due to the change of mechanical material properties. Moreover, the free mechanical load conditions at the external surfaces are satisfied in terms of the ${\sigma }_{zz}$ stress component.

The third benchmark considers a four-layered composite $0^{\circ }/{90}^{\circ }/0^{\circ }/{90}^{\circ }$ cylinder with simply supported sides (see Figure 1). All the geometrical characteristics are given in Table 2 and the employed material is shown in Table 3 (Composite One column). The moisture content profiles are shown in Figure 6 for thick ($R_{\alpha }/h=4$) and thin ($R_{\alpha }/h=500$) structures. In this benchmark, the thickness layer effect for thick cylinders is negligible due to the symmetry and rigidity of the closed cylindrical geometry: the difference between (${\mathcal{M}}_c$, 3D) or (u-$\mathcal{M}$) profiles and (${\mathcal{M}}_c$, 1D) or (${\mathcal{M}}_a$) profiles for the thick case is not so evident. These conclusions are still valid for the displacement, stress, and strain components discussed in

Table 8. Benchmark number three, $0^{\circ }/{90}^{\circ }/0^{\circ }/{90}^{\circ }$ composite cylinder with moisture content ${\mathcal{M}}_t=+1\%$ and ${\mathcal{M}}_b=+0.5\%$ for $m=2$ and $n=1$. Three-dimensional uncoupled hygro-elastic models from [21] with the moisture content profile assumed to be linear (3D(${\mathcal{M}}_a$)), defined via the one-dimensional Fick diffusion relation (3D (${\mathcal{M}}_c$, 1D)) and via the three-dimensional Fick diffusion relation (3D (${\mathcal{M}}_c$, 3D)).

${\mathit{R}}_{\mathit{\alpha }}/\mathit{h}$	4	10	50	100	500
$\mathcal{M}$(%) in $\tilde{z}=+h/2$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_a$) [21]	0.7500	0.7500	0.7500	0.7500	0.7500
3D (${\mathcal{M}}_c$, 1D) [21]	0.7500	0.7500	0.7500	0.7500	0.7500
3D (${\mathcal{M}}_c$, 3D) [21]	0.7354	0.7476	0.7499	0.7500	0.7500
3D-u-$\mathcal{M}$	0.7354	0.7476	0.7499	0.7500	0.7500
v (mm) in $\tilde{z}=+3h/4$, $\alpha =a/2$ and $\beta =0$
3D (${\mathcal{M}}_a$) [21]	−1.133	−1.252	−1.281	−1.283	−1.284
3D (${\mathcal{M}}_c$, 1D) [21]	−1.133	−1.252	−1.281	−1.283	−1.284
3D (${\mathcal{M}}_c$, 3D) [21]	−1.117	−1.249	−1.281	−1.283	−1.284
3D-u-$\mathcal{M}$	−1.117	−1.249	−1.281	−1.283	−1.284
w (mm) in $\tilde{z}=+h/2$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_a$) [21]	2.109	2.479	2.672	2.696	2.715
3D (${\mathcal{M}}_c$, 1D) [21]	2.109	2.479	2.672	2.696	2.715
3D (${\mathcal{M}}_c$, 3D) [21]	2.080	2.473	2.672	2.696	2.715
3D-u-$\mathcal{M}$	2.080	2.473	2.672	2.696	2.715
${\sigma }_{zz}$ (${10}^3$ Pa) in $\tilde{z}=+h/4$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_a$) [21]	−680.9	172.3	95.63	51.54	10.90
3D (${\mathcal{M}}_c$, 1D) [21]	−680.9	172.3	95.63	51.54	10.90
3D (${\mathcal{M}}_c$, 3D) [21]	−664.4	172.2	95.62	51.54	10.90
3D-u-$\mathcal{M}$	−664.4	172.2	95.62	51.54	10.90
${\sigma }_{\alpha \beta }$(${10}^3$ Pa) in $\tilde{z}=0$, $\alpha =0$ and $\beta =0$
3D (${\mathcal{M}}_a$) [21]	−703.8	−164.1	−24.90	−12.15	−2.389
3D (${\mathcal{M}}_c$, 1D) [21]	−703.8	−164.1	−24.90	−12.15	−2.389
3D (${\mathcal{M}}_c$, 3D) [21]	−692.9	−163.7	−24.90	−12.15	−2.389
3D-u-$\mathcal{M}$	−692.9	−163.7	−24.90	−12.15	−2.389
${\gamma }_{\beta z}$(${10}^{-5}$) in $\tilde{z}=+h/3$, $\alpha =a/2$ and $\beta =0$
3D (${\mathcal{M}}_a$) [21]	15.59	6.534	1.258	0.6243	0.1241
3D (${\mathcal{M}}_c$, 1D) [21]	15.59	6.534	1.258	0.6243	0.1241
3D (${\mathcal{M}}_c$, 3D) [21]	15.46	6.525	1.257	0.6243	0.1241
3D-u-$\mathcal{M}$	15.46	6.525	1.257	0.6243	0.1241

, where 3D (${\mathcal{M}}_a$), 3D (${\mathcal{M}}_c$, 1D), 3D (${\mathcal{M}}_c$, 3D), and 3D-u-$\mathcal{M}$ models provide the same results for thin cylinders and small distinctions for very thick structures: the thickness layer effect in the thick configuration is not so important thanks to the great rigidity of the closed and symmetric cylinder configuration. Figure 7 shows displacements, stresses, and strains through the thickness direction of the four-layered composite cylinder. The displacements and transverse normal stress shown in Figure 7 are continuous, even if the material properties change in the thickness direction because of the proper congruence and equilibrium impositions. The same is valid for the moisture content profile. The transverse normal stress satisfies the imposed external load conditions. The in-plane normal stress and transverse shear strain show discontinuities in correspondence with the layer interfaces.

The last benchmark investigates a composite three-layered spherical shell panel with simply supported sides (see Figure 1). Geometrical data are available in column B4 of Table 2 and the material properties are written in Table 3 (Composite One column). The results in Figure 8 demonstrate that the moisture content profile is linear when the spherical shell panel is thin. The material layer effect is not visible because the material of each physical layer is always the same and the change in the fiber orientation angle from $0^{\circ }$ to ${90}^{\circ }$ has no effects on material properties. This feature means that only the thickness layer effect is visible. For the thick configurations, the (${\mathcal{M}}_c$, 3D) and (u-$\mathcal{M}$) profiles give different trends with respect to the (${\mathcal{M}}_c$, 1D) and (${\mathcal{M}}_a$) cases because the thickness layer effect is evident. The displacement, stress, and strain components for different thickness ratios $R_{\alpha }/h$ are given in

Table 9. Benchmark number four, $0^{\circ }/{90}^{\circ }/0^{\circ }$ composite spherical shell with moisture content ${\mathcal{M}}_t=+1\%$ and ${\mathcal{M}}_b=+0.5\%$ for $m=n=1$. Three-dimensional uncoupled hygro-elastic models from [21] with the moisture content profile assumed to be linear (3D (${\mathcal{M}}_a$)), defined via the one-dimensional Fick diffusion relation (3D (${\mathcal{M}}_c$, 1D)) and via the three-dimensional Fick diffusion relation (3D (${\mathcal{M}}_c$, 3D)).

${\mathit{R}}_{\mathit{\alpha }}/\mathit{h}$	4	10	50	100	500
$\mathcal{M}$(%) in $\tilde{z}=+h/2$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_a$) [21]	0.7500	0.7500	0.7500	0.7500	0.7500
3D (${\mathcal{M}}_c$, 1D) [21]	0.7500	0.7500	0.7500	0.7500	0.7500
3D (${\mathcal{M}}_c$, 3D) [21]	0.6390	0.7301	0.7492	0.7498	0.7500
3D-u-$\mathcal{M}$	0.6390	0.7301	0.7492	0.7498	0.7500
v (mm) in $\tilde{z}=+h/3$, $\alpha =a/2$ and $\beta =0$
3D (${\mathcal{M}}_a$) [21]	−0.8904	−0.6471	−0.2607	−0.2528	−0.2663
3D (${\mathcal{M}}_c$, 1D) [21]	−0.8904	−0.6471	−0.2607	−0.2528	−0.2663
3D (${\mathcal{M}}_c$, 3D) [21]	−0.7641	−0.6328	−0.2606	−0.2528	−0.2663
3D-u-$\mathcal{M}$	−0.7641	−0.6328	−0.2606	−0.2528	−0.2663
w (mm) in $\tilde{z}=+h/2$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_a$) [21]	−0.1617	0.3608	1.651	1.703	1.681
3D (${\mathcal{M}}_c$, 1D) [21]	−0.1617	0.3608	1.651	1.703	1.681
3D (${\mathcal{M}}_c$, 3D) [21]	−0.0540	0.3680	1.650	1.703	1.681
3D-u-$\mathcal{M}$	−0.0540	0.3680	1.650	1.703	1.681
${\sigma }_{zz}$ (${10}^3$ Pa) in $\tilde{z}=+h/3$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_a$) [21]	−1080	−553.0	−59.78	−20.61	−2.559
3D (${\mathcal{M}}_c$, 1D) [21]	−1080	−553.0	−59.78	−20.61	−2.559
3D (${\mathcal{M}}_c$, 3D) [21]	−958.4	−544.2	−59.76	−20.60	−2.559
3D-u-$\mathcal{M}$	−958.4	−544.2	−59.76	−20.60	−2.559
${\sigma }_{\beta z}$(${10}^3$ Pa) in $\tilde{z}=+2h/3$, $\alpha =a/2$ and $\beta =0$
3D (${\mathcal{M}}_a$) [21]	−4584	−2144	−463.4	−233.6	−46.95
3D (${\mathcal{M}}_c$, 1D) [21]	−4584	−2144	−463.4	−233.6	−46.95
3D (${\mathcal{M}}_c$, 3D) [21]	−4312	−2121	−463.1	−233.6	−46.95
3D-u-$\mathcal{M}$	−4312	−2121	−463.1	−233.6	−46.95
${\gamma }_{\alpha \alpha }$(${10}^{-4}$) in $\tilde{z}=+h$, $\alpha =a/2$ and $\beta =b/2$
3D (${\mathcal{M}}_a$) [21]	2.998	2.315	1.370	1.118	0.9025
3D (${\mathcal{M}}_c$, 1D) [21]	2.998	2.315	1.370	1.118	0.9025
3D (${\mathcal{M}}_c$, 3D) [21]	2.801	2.280	1.369	1.117	0.9025
3D-u-$\mathcal{M}$	2.801	2.280	1.369	1.117	0.9025

, where it is clear that for thick cases, the thickness layer effect is the only effect acting due to the particular material configuration. Displacement and stress component evaluations through the thickness direction shown in Figure 9 remark these effects. Moisture content, in-plane displacement, transverse displacement, and transverse shear stress components are not discontinuous because the congruence and equilibrium conditions have been correctly imposed in the mathematical formulation. In-plane strain shows two discontinuities in correspondence with the two physical interfaces where the composite material changes the fiber orientation. Typical zigzag effects due to the transverse anisotropy are shown in Figure 9 and load boundary conditions are fully satisfied in terms of the ${\sigma }_{zz}$ stress component.

5. Conclusions

A fully coupled hygro-elastic 3D exact shell model for the static analysis of one-layered and multilayered isotropic, sandwich, and composite plates and shells has been discussed. The moisture content environment influences the external faces of the structures under steady-state conditions. The moisture content profile is evaluated through the thickness direction. The moisture content profile is a primary variable of the problem as the displacements: the 3D Fick equation and the 3D equilibrium equations for the shells are coupled. The 3D Fick equation allows a moisture content profile able to consider both the thickness layer and material layer effects for each possible configuration. The coupled system is solved in a closed form thanks to the exponential matrix method and harmonic forms for each variable. Several results, in terms of displacements, in-plane and out-of-plane stresses, and moisture content profiles have been discussed for different thickness ratios, geometrical properties, lamination sequences, moisture impositions, and materials. These analyses showed a complete match between the 3D uncoupled model that uses the 3D Fick diffusion relation and the 3D fully coupled model developed in this paper. This new method takes into account both the thickness and the material layer effects using a mathematical formulation that is simpler and more elegant because the 3D Fick diffusion relation is not separately solved via an external method. Moreover, a reduced number of mathematical layers F is requested in comparison with the uncoupled 3D model.