Thermomechanical Buckling Analysis of the E&P-FGM Beams Integrated by Nanocomposite Supports Immersed in a Hygrothermal Environment

Due to the widespread use of sandwich structures in many industries and the importance of understanding their mechanical behavior, this paper studies the thermomechanical buckling behavior of sandwich beams with a functionally graded material (FGM) middle layer and two composite external layers. Both composite skins are made of Poly(methyl methacrylate) (PMMA) reinforced by carbon-nano-tubes (CNTs). The properties of the FGM core are predicted through an exponential-law and power-law theory (E&P), whereas an Eshelby–Mori–Tanaka (EMT) formulation is applied to capture the mechanical properties of the external layers. Moreover, different high-order displacement fields are combined with a virtual displacement approach to derive the governing equations of the problem, here solved analytically based on a Navier-type approximation. A parametric study is performed to check for the impact of different core materials and CNT concentrations inside the PMMA on the overall response of beams resting on a Pasternak substrate and subjected to a hygrothermal loading. This means that the sensitivity analysis accounts for different displacement fields, hygrothermal environments, and FGM theories, as a novel aspect of the present work. Our results could be replicated in a computational sense, and could be useful for design purposes in aerospace industries to increase the tolerance of target productions, such as aircraft bodies.


Introduction
FGM-based structures serve as bi-phase beams, plates and shells whose properties vary across their thickness or length continuously. Over the last century, FGMs have increasingly attracted the interest of the scientific community for their use in different high-tech industries due to their outstanding mechanical properties compared to singlephase materials and structures. Some natural types of FGMs are visible, such as bamboo trees, teeth, bones, and human skin, which have evolved to meet a specific requirement in humans and their environment. For the first time, Shen and Bever [1] considered graded material composites, despite their limited knowledge and a general lack of sophisticated fabrication equipment. Subsequently, Japanese scientists in 1984 [2] applied this technique to an aerospace project, due to their necessity to have a 10 mm thickness thermal barrier with an internal and external temperature of 1000 K and 2000 K, respectively. Fu et al. [3] also studied the thermoacoustic response of a porous FGM cylindrical plate with a random distribution of pores. Moreover, Duc and Cong [4] studied the vibrational behavior of FGM plates using a Runge-Kutta method, whose model was affected by thermomechanical coupled loading conditions. In another work, Shen [5] studied the thermomechanical post-buckling behavior of a simply-supported FGM plate equipped with piezoelectric fiber-reinforced composite patches for sensors, actuators, transducers and active damping devices. Zenkour [6] studied the static deflection response of an FGM plate in a hygrothermal environment, while discussing the response sensitivity to different environmental conditions. Transient analysis of porous FGM plates was also considered by Van et al. [7] in a nonlinear domain subjected to a coupled hygrothermal and mechanical loading condition. A post-buckling study of sandwich plates with FGM skins was performed by Kiani and Eslami [8] in a numerical sense, focusing on the effect of the power-law index, foundation parameters and imperfections of the overall response. In addition to the theoretical studies, many researchers have experimentally investigated the vibration and damping properties of nanoparticle-reinforced composite materials [9][10][11][12].
Besides FGMs, other composite nanomaterials have attracted the interest of engineering device production, with the rapid development of various analytical and computational models designed to simulate their behavior even in a nonlocal sense. In this context, Arshid et al. [13] studied the dynamic and static behavior of annular FG graphene nanoplatelets (FG-GNPs) reinforced nanocomposite with porosities, and applied modified strain gradient theory (MSGT) to account for size-dependent effects. In another work, Foroutan et al. [14] surveyed the nonlinear buckling and vibration behavior of imperfect FG carbon nanotubes reinforced composite (FG-CNTRC) cylindrical shells in a hygrothermal environment. Similarly, Arshid et al. [15] applied a 3D plate theory to study the vibrations of sandwich structures with a honeycomb core and GNP-reinforced epoxy face sheets, as commonly found in many electric devices. Safaei [16] developed a generalized vibration model including the possible presence of porosity within sandwich structures (both in core and skins) immersed in a Pasternak foundation. Among the recent scientific literature, several theoretical and numerical methods have been proposed to solve complicated structural problems, involving advanced composite materials. Moradi-Dastjerdi et al. [17] applied a higher-order theory and the Eshelby-Mori-Tanaka (EMT) approach to assess the buckling behavior of a sandwich plate made of a CNTRC porous core with external piezoelectric face sheets, while assuming different CNTs' agglomerations within the material. In a similar direction, the authors of Refs. [18,19] proposed a refined higher-order theory to evaluate the CNTs' agglomeration impact on the vibration control and stiffness of nanocomposite sandwich beams, respectively. In addition, further efforts in this direction have been performed by other scholars .
Motivated by the aforementioned studies, this paper aims to further contribute to the thermomechanical response of sandwich beams including an exponential-law/powerlaw-based functionally graded core (E&P-FGC) integrated by two CNTRC layers. A unified framework was originally proposed to account for different high-order kinematic assumptions for the systematic study of sandwich structures resting on an elastic Pasternak substrate, under coupled hygrothermal loading conditions. The governing equations of the problem are determined from the virtual displacement principle, accounting for various beam theories and CNTs' agglomeration effects. A Navier-type procedure was, thus, proposed to solve the problem analytically, whose results, based on different E&P-FGM relations and thermomechanical conditions could be very useful for further computational investigations on this topic, even from a practical design standpoint.

Analytical Model
Let us consider a three-layer sandwich beam with length a and thickness h. As shown in Figure 1, the total thickness is the assemblage of an FGC, and two external (top and bottom) CNTRC skins, with thickness h c , h t , and h b , respectively. The model is referred to as the cartesian coordinate system (x, y, z), and it is located at the midplane of the model, as depicted in Figure 1. The whole structure is embedded in an elastic Pasternak substrate, and it is subjected to a variable hygrothermal surrounding condition.
as the cartesian coordinate system (x, y, z), and it is located at the midplane of the model, as depicted in Figure 1. The whole structure is embedded in an elastic Pasternak substrate, and it is subjected to a variable hygrothermal surrounding condition. For this sandwich structure, we propose and compare a thermomechanical model based on four different kinematic assumptions, namely, the third-order shear deformation theory (TSDT), the hyperbolic shear deformation theory (HSDT), the sinusoidal shear deformation theory (SSDT) and the exponential shear deformation theory (ESDT). Based on these kinematic assumptions, the displacement field for an arbitrary point in the mid surface is defined as [20,21].
where u and w refer to the longitudinal and transverse displacement components of the mid surface. For simplicity reasons, Equation (1) can be rewritten in a more compact For this sandwich structure, we propose and compare a thermomechanical model based on four different kinematic assumptions, namely, the third-order shear deformation theory (TSDT), the hyperbolic shear deformation theory (HSDT), the sinusoidal shear deformation theory (SSDT) and the exponential shear deformation theory (ESDT). Based on these kinematic assumptions, the displacement field for an arbitrary point in the mid surface is defined as [20,21].
where u and w refer to the longitudinal and transverse displacement components of the mid surface. For simplicity reasons, Equation (1) can be rewritten in a more compact notation by introducing a general function Q(z) to define the terms in the square brackets, depending on the selected kinematic theory. Thus, by simply changing the definition of Q(z), we are able to switch among different displacement assumptions.
The kinematic relations between the strain and displacement fields are, thus, expressed as [22]: where superscripts c, t, and b refer to the core, top, and bottom layer. For every single layer of the sandwich structure, we refer to the following constitutive relations [23,24]: where A refers to the elastic constants, β and α denote the moisture and thermal expansion coefficients, respectively, ∆T and ∆H stand for the thermal and moisture variation, respectively, defined as ∆T = T(z) − T 0 and ∆H = H(z) − H 0 . Moreover, T 0 and H 0 denote the ambient temperature and moisture, here kept equal to 293 K and 0.1 wt%H 2 O, respectively. Furthermore, T(z) and H(z) represent a linear temperature and moisture variation across the thickness direction of the beam [25]: where subscripts t and b refer to the top and bottom surface of the model. At the same time, ∆T tb and ∆H tb stand for the thermal and moisture variation between the top and bottom sides, namely, ∆T tb = T t − T b and ∆H tb = H t − H b . The core layer of the selected sandwich beam is made up of FGMs, which means that the top and bottom surfaces of the FGC are made of pure ceramic and pure metal, whereas from the bottom to the top surfaces, the volume fraction of the metal and ceramic phases vary, keeping fixed V c + V m = 1. In what follows, we consider two different variations of material properties throughout the thickness of the FGC, evaluated here comparatively for their impact on the thermomechanical buckling behavior of the whole structure. Based on the P-FGC and E-FGC, the mechanical properties of the FGC are defined as: where subscripts c and m denote the ceramic and metallic properties, respectively. Moreover, J is a general notation that defines the arbitrary mechanical property, as the elasticity modulus E, Poisson's ratio ν, density ρ, β and α. In the case of P-FGC, S is the power-law index. Allocating a zero value to S, the FGC reverts to a fully ceramic layer; when S = ∞ the core layer is fully metallic. Based on the two-fold definition of J c (z) in Equation (5), the elasticity modulus varies in the FGC thickness as plotted in Figure 2. The sandwich beam system includes two CNTRC faces, with increased overall stiffness. Both external skins are assumed to be completely bonded to the FGC without any possible interlayer slip. The presence of CNTs within the PMMA matrix plays an efficient role in the general improvement of the mechanical properties. Due to the high aspect ratio of CNTs, their bending and bundling within the matrix phase and clusters formation seem to be almost predictable. In this paper, we consider the effect of the CNTs' agglomeration on the overall mechanical response, based on an EMT scheme [26]. Based on this approach, the bulk modulus K and shear modulus G inside and outside the clusters can be defined as [26,27]: where subscripts p and r refer to the polymeric matrix and CNTs. when = ∞ the core layer is fully metallic. Based on the two-fold definition of J c (z) in Equation (5), the elasticity modulus varies in the FGC thickness as plotted in Figure 2. The sandwich beam system includes two CNTRC faces, with increased overall stiffness. Both external skins are assumed to be completely bonded to the FGC without any possible interlayer slip. The presence of CNTs within the PMMA matrix plays an efficient role in the general improvement of the mechanical properties. Due to the high aspect ratio of CNTs, their bending and bundling within the matrix phase and clusters formation seem to be almost predictable. In this paper, we consider the effect of the CNTs' agglomeration on the overall mechanical response, based on an EMT scheme [26]. Based on this approach, the bulk modulus K and shear modulus G inside and outside the clusters can be defined as [26,27]: Furthermore, µ and η denote the CNTs' agglomeration region volume fraction and the cluster's CNTs' volume fraction. It should be mentioned that the value of µ can be lower or equal to η. Regarding the definition and value assumed by µ and η, some possible cases can occur, as detailed in Table 1. Table 1. Possible combinations of µ and η.

Case
Result The whole composite layer serves as a big agglomerated region.
There is no CNT inside the clusters. η = 1 All CNTs are agglomerated inside the clusters.
In order to capture the impact of different CNT distribution patterns inside the composite faces, the following equations can be used to define the CNT volume fraction inside the composite skins [28]: distributions. Also, it should be mentioned that V* is the CNTs' volume fraction in the case of their uniform dispersion. In addition, the terms α r , β r , η r and δ r can be written as follows [17]: where k r , l r , m r , p r and n r denote the CNTs' elastic Hill's constants which are varied for different types of CNTs with respect to their chiral index value. In the current study, such constants are defined for single-walled carbon nanofillers (SWCNTs) with a chirality index equal to 10 [29].
Molecules 2021, 26, 6594 7 of 19 where kr, lr, mr, pr and nr denote the CNTs' elastic Hill's constants which are varied for different types of CNTs with respect to their chiral index value. In the current study, such constants are defined for single-walled carbon nanofillers (SWCNTs) with a chirality index equal to 10 [29]. Thus, the bulk modulus and shear modulus of the external skins can be obtained as [30]: (1 ) CNTRC faces thickness Thus, the bulk modulus and shear modulus of the external skins can be obtained as [30]: in which subscript n denotes the nanocomposite faces, and: where, The mechanical properties of the composite faces including E, ν, ρ, α and β can be addressed as [17,31]: Finally, the elastic constants for the external layers can be defined as [32]:

Governing Equations
We now apply the virtual work principle to determine the governing equations of the problem [15]: where Ω is the strain energy and Ψ is the external work. The classical strain energy for a sandwich beam is defined as [33]: The external work includes three terms, namely, the thermal force, the external force related to humidity and the elastic substrate force, defined as follows [34]: where L 1 is the spring constant and L 2 stands for the shear layer constant. The total external work reads as follows [35]: By substitution of Equations (26) and (28) into Equation (25), after mathematical manipulation, we obtain the following governing equations:

Analytical Solution Procedure
Equations (29)-(31) are solved here analytically by means of a Navier-type procedure. Based on this technique, we use the following expressions to define the displacement components, which satisfy simply supported boundary conditions [36]: where, U, U 1, and W serve as unknown coefficients; m denotes the mode number along the ESB length. By placing Equation (33) into Equations (29)-(31), after mathematical manipulation, the governing equations are easily solved.

Numerical Results and Discussions
In this section, we perform various numerical examples to demonstrate the accuracy of the formulation proposed herein, to study the critical buckling load of an Euler-Bernoulli beam (EBB) and Timoshenko beam (TB), for different values of power-law exponent and geometrical characteristics. Our results are compared to predictions by Li and Batra [37], as summarized in Table 2, with an acceptable agreement for the geometrical length-tothickness ratios a/h = 5 and 10, and for all values of S. This confirms the reliability of the proposed formulation to handle such a topic.
A systematic analysis is repeated, by assuming different input material and geometrical properties in the model, as listed in Tables 3-5, in line with Refs. [17,38], in which the term N* stands for the nondimensional critical buckling load.  Table 4. CNTRC skin material properties [17].  In Figure 4, we plot the variation of the dimensionless critical buckling load N* vs. the CNTs' volume fraction inside clusters, η, while setting different combinations of the foundation parameters L 1 , L 2 .
Based on results in this figure, CNTs' inflation inside some concentrated regions (i.e., an increased value η) yields a reduced critical buckling load because of the stiffness reduction in the whole structure, in line with the experimental findings by Pan et al. [11]. On the other hand, by increasing the foundation parameters L 1 , L 2 , the critical buckling load is enhanced because of an increased stiffness configuration induced on the structure.
As also plotted in Figure 5, the model is very sensitive to the rational length a/h and to the CNTs' agglomeration region µ, with a monotonic reduction of N* for an increased a/h ratio, along with a gradual reduction of N* for a reduced level of µ, under a fixed value of a/h. This means that a decreased agglomeration region of CNTs causes a stiffness reduction of the structure, which can buckle more easily. In Figure 4, we plot the variation of the dimensionless critical buckling load N* vs. the CNTs' volume fraction inside clusters, η, while setting different combinations of the foundation parameters L1, L2. Based on results in this figure, CNTs' inflation inside some concentrated regions (i.e., an increased value η) yields a reduced critical buckling load because of the stiffness reduction in the whole structure, in line with the experimental findings by Pan et al. [11]. On the other hand, by increasing the foundation parameters L1, L2, the critical buckling load is enhanced because of an increased stiffness configuration induced on the structure.
As also plotted in Figure 5, the model is very sensitive to the rational length a/h and to the CNTs' agglomeration region μ, with a monotonic reduction of N* for an increased a/h ratio, along with a gradual reduction of N* for a reduced level of μ, under a fixed value of a/h. This means that a decreased agglomeration region of CNTs causes a stiffness reduction of the structure, which can buckle more easily.   The study considers the effect of different top-bottom moisture variations ΔHtb based on the E&P-FGC models, while varying the core thickness hc, as represented in Figure 6. Based on the plots in Figure 6, it is visible that the application of a P-FGC model obtains higher values of the buckling load and structural stiffness compared to an E-FGC model, for a fixed value of ΔHtb. As also expected, an increased level of humidity in the surrounding environment has a destructive effect on the mechanical properties of the structure with a consecutive reduction of N* under the same assumption of hc. This mechanical decay for different humidity conditions represents a key aspect for many aerospace and marine applications, as well as nanostructures like generators and sensors. At the same time, an increased core thickness makes the structure stiffer, thus leading to a monotonic increase in the buckling load for each selected model and moisture variation.
We now repeat a similar investigation based on the E&P-FGC models, while checking for the sensitivity of the response to the thermal variation ΔTtb among the top and bottom The study considers the effect of different top-bottom moisture variations ∆H tb based on the E&P-FGC models, while varying the core thickness h c , as represented in Figure 6. Based on the plots in Figure 6, it is visible that the application of a P-FGC model obtains higher values of the buckling load and structural stiffness compared to an E-FGC model, for a fixed value of ∆H tb . As also expected, an increased level of humidity in the surrounding environment has a destructive effect on the mechanical properties of the structure with a consecutive reduction of N* under the same assumption of h c . This mechanical decay for different humidity conditions represents a key aspect for many aerospace and marine applications, as well as nanostructures like generators and sensors. At the same time, an increased core thickness makes the structure stiffer, thus leading to a monotonic increase in the buckling load for each selected model and moisture variation.   We now repeat a similar investigation based on the E&P-FGC models, while checking for the sensitivity of the response to the thermal variation ∆T tb among the top and bottom sides of the sandwich structure together with a varying thickness of the external skins (see Figure 7). Based on the plots in this figure, an increased thickness of the top and bottom sides enhances the overall stiffness of the structure together with its buckling load, independently of the selected model and thermal variation. Moreover, for each selected model, the results are almost unaffected by the thermal variation ∆T tb , where a P-FGC model always obtains higher values of the buckling load with respect to an E-FGC model. Figure 8, instead, aims at evaluating the sensitivity of N* for different CNT volume fractions V* and different FGC material constituents. The monotonic increase in each curve for an increased level of V* demonstrates the beneficial effect of such a quantity on the overall buckling response of the structure. It was also found that the use of the Al in conjunction with Al 2 O 3 , Si 3 N 4 and ZrO 2 as FGC material constituents, obtains higher magnitudes of the buckling load.
As visible in Figure 9, an enhanced power-law index causes a general reduction of dimensionless critical buckling. Among different possibilities of the CNTs' distributions, the FGA-V and FGV-A distributions yield the lowest and highest magnitudes of the buckling load, respectively, where other types of distribution yield intermediate predictions. Figure 10 displays the critical buckling temperature against the volume fraction of CNTs for different CNTs' volume fractions, η. As visible in this figure, an increased value of V* improves the critical buckling temperature, especially for a reduced level of η.    As visible in Figure 9, an enhanced power-law index causes a general reduction of dimensionless critical buckling. Among different possibilities of the CNTs' distributions, the FGA-V and FGV-A distributions yield the lowest and highest magnitudes of the buckling load, respectively, where other types of distribution yield intermediate predictions. Figure 10 displays the critical buckling temperature against the volume fraction of CNTs for different CNTs' volume fractions, η. As visible in this figure, an increased value of V* improves the critical buckling temperature, especially for a reduced level of η.
At the same time, the critical buckling temperature seems to reduce for an increased length of the beam, as plotted in Figure 11, for different values of the power-law index S. Moreover, the overall thermal response of the structure, seems to be slightly affected by S, whose increase produces a meaningless reduction in Tcr under the same geometrical length.
The final 3D plot in Figure 12 shows the double variation of Tcr with the core thickness hc and substrate elastic constant L1. Note that the double increase in hc and L1 produce the      At the same time, the critical buckling temperature seems to reduce for an increased length of the beam, as plotted in Figure 11, for different values of the power-law index S. Moreover, the overall thermal response of the structure, seems to be slightly affected by S, whose increase produces a meaningless reduction in T cr under the same geometrical length.   The final 3D plot in Figure 12 shows the double variation of T cr with the core thickness h c and substrate elastic constant L 1 . Note that the double increase in h c and L 1 produce the highest value of T c . The opposite effect is obtained for a simultaneous decrease in both parameters, which corresponds to the most serious thermal buckling condition.   Finally, in Table 6, we summarize the impact of different beam theories on the dimensionless critical buckling load, while considering different FGC models and CNT volume fractions, in a numerical sense. It seems that TSDT and HSDT provide the highest and lowest predictions, respectively, whereas SSDT and ESDT always produce intermediate results. It is also found that P-FGC-based results are much higher than those stemming from an E-FGC model, for a fixed value of µ. Furthermore, a gradual increase in the buckling load is observed for an increased value of µ, independent of the selected model.

Conclusions
This paper presents a thermomechanical buckling analysis of sandwich beams with two identical CNTRC skins, under different hygrothermal environmental conditions and elastic foundation parameters, based on an E&P-FGC model. Different FGM theories are employed here in a unified framework for the first time, in conjunction with different CNTs' agglomeration assumptions. By using the virtual displacement approach and Navier-type solution, we determine and solve the governing equations of the problem. Based on a large systematic investigation aimed at determining the thermomechanical buckling response and its sensitivity under different input parameters, some considerable conclusions can be summarized, as follows: