Hygrothermal environment effect on the critical buckling load of FGP microbeams with initial curvature integrated by CNT-reinforced skins considering the influence of thickness stretching


 Due to the need for structures with refined properties to bear against different loading conditions, recently, carbon nanotubes (CNTs) have been used widely to reinforce them. The extremely high stiffness of CNTs makes them significant as one of the best reinforcements to improve the mechanical behaviors of structures. This work focuses on microbeam buckling response with an initial curvature that includes three layers. The mid-layer that is known as the core is constituted of functionally graded porous (FGP) materials and two CNT-reinforced composite skins are bonded to the core to integrate it. The whole structure is affected by the hygrothermal environment and springs and shear layers are put below it. For the first time, for such a structure, a refined shear deformation theory (RSDT) as a higher-order theory that considers thickness stretching effect in polar coordinates is used that presents more accurate results, especially for deeply curved beams. Modified couple stress theory (MCST) in combination with the virtual displacement principle is utilized to establish the governing equations. The obtained results demonstrate the significance of porosity percentage and CNTs’ addition to the skins on the critical nanotubes buckling load. Also, the different behaviors of the microstructure at various temperatures are analyzed and discussed in detail.


Introduction
Nowadays, scientists are searching for ways to improve the mechanical behaviors of different kinds of structures against different loading conditions to get the desired response. Nevertheless, with the traditional materials and structures, reaching this aim is not possible. Therefore, the researchers were motivated to try different combinations of both materials and structures to analyze their mechanical behaviors [1][2][3][4][5][6][7][8][9]. One of these ways is using composite materials. Composites due to their more proper specifications are used for decades [10], but they have limitations too. Using nanofillers as the second phase of the composite is one of the most appropriate ways to reinforce the structures. These kinds of composites that are reinforced by nanofillers are known as nanocomposites [11][12][13][14][15]. The nanofillers are distributed within the matrix by different methods [16]. Among the different kinds of nanofillers, carbon nanotubes (CNTs) attracted researchers more than others due to their unique properties as high stiffness, sound and energy-absorbing specifications, thermal conductivity, and are known for their electromagnetic characteristics [17][18][19]. It may be interesting to know that CNTs were discovered first by Iijima in 1991 [20]. After that, the rapid development of CNT studies has occurred and different aspects of their specifications are considered by many authors [21][22][23][24][25][26][27][28]. Also, another way to improve the behavior of the structures is to use several layers together. This kind of structure is usually called sandwich structure and includes a mid-layer or a core, and two skins that are bonded to the core. The use of sandwich structures leads to reach approximately the desired results. It is noteworthy that the core of the sandwich is selected in such a way that it is softer than the skins. Also, its thickness is much more than that of skins. In the field of these kinds of engineering structures, researchers are publishing their new findings [29][30][31][32][33][34][35][36]. As an instance, vibration analysis of a CNT-reinforced composite plate was performed by Zhu et al. [37]. They applied the finite element analysis by using the first-order shear deformation plate theory (FSDT). They considered different types of CNT distributions and boundary conditions. In another study, the effect of piezoelectric actuators that were fully bonded to a porous core on the vibrational behavior of the system was considered by Arshid and Khorshidvand [38]. Loghman and Cheraghbak [39] addressed agglomeration effects on the behavior of the nanocomposite of the piezoelectric cylinder. They assumed that the structure was under internal pressure and CNTs were selected as the reinforcement. A refined mixture rule was used by Mirzaei and Kiani [40] to analyze cylindrical panels made from functionally graded-CNT (FG-CNT)-reinforced composites. They investigated the different effects of CNTs on the results. Stability analysis of the three-phase sandwich panels that were nonlinearly imperfect was presented by Thu and Duc [41]. Using the Ritz numerical approach, Selim et al. [42] attained natural frequencies of FG-CNT-reinforced composites plates. They included the thermal effects on the results and investigated the uniform and FG distribution of CNTs and presented their results. The dynamic response of the piezoelectric FG sandwich cylindrical shells was investigated by Nguyen [43] considering the shear deformation effect. Also, Mehar et al. [44] applied a finite element method to solve the nonlinear equations of a doubly curved FGCNT-reinforced composites shell. The structure was modeled by using higher-order shear deformation theory with Green-Lagrange strains. Yas and Samadi [45] explored the CNT-reinforced composite Timoshenko beams by investigating the vibration response and buckling load using the extended rule of the mixture to achieve effective properties.
In recent years, many kinds of research have illustrated that when the characteristic length scale related to nonuniform plastic deformation is in the micro/nanoscale, materials show strong size effects [46][47][48][49][50]. The theories of classical plasticity cannot be used to calculate this size dependence of material behavior at small scales because their constitutive models do not possess an inner length scale. Since investigating the size-dependent mechanical behavior of nano/microscaled structural systems can deliver important statistics for designers of micro/nanoscaled mechanisms, some macroscopic continuum mechanics were effectively industrialized and selected to evaluate the size-dependent effect. Thus, to give details about the size effects, developing strain gradient theory for small-scale devices is required [51][52][53]. Utilizing strain gradient theory, Hosseini et al. [54] analyzed various stress distributions in single-wall carbon nanotubes (SWCNTs) under internal pressure to obtain the size-dependent behavior of SWCNTs. Malikan and Nguyen [55] obtained buckling of piezo-magneto-electric nanoplates in a hygro-thermal atmosphere by applying the higher-order nonlocal strain gradient theory and a new variable plate method. Imani Aria and Biglari [56] proposed a nonlocal strain gradient model to investigate vibrational responses and buckling of microtubules utilizing the finite element method to investigate different size effect parameters. Arshid et al. [57] used the modified couple stress theory (MCST) to analyze both vibrational response and buckling behavior of curvy microbeams. Their under consideration model was composed of a porous core that was embedded between two nanocomposite patches. MCST was used by Jouneghani et al. [58] to analyze nonuniform composite laminated beams based on the Ritz formulation. They used Euler-Bernoulli beams theory to describe the displacements. Also, hygrothermal modeling of the buckling behavior of sandwich plates with nanocomposite face sheets resting on a Pasternak foundation was provided by Kiarasi et al. [59]. Soltani et al. [60] presented their study about the nonlocal elasticity theory for lateral stability of tapered thinwalled nanobeams with axially varying materials. In another study, Karami et al. [61] analyzed nonlocal buckling of composite curved beams, which were reinforced with FGCNTs. Size-dependent free vibrations of the FG polymer composite curved nanobeams reinforced with nanocomposite resting on Pasternak foundations were investigated by Arefi et al. [62]. It should be noted that they used the Cartesian coordinate system to perform their study, which may lead to inaccurate results, especially for higher values of the central angle. To overcome this shortcoming of the use of the Cartesian coordinate system, recently, the authors have been trying to analyze curved beams via the polar coordinate system. Sobhy [63] in his published work has used the polar coordinate system, natural frequencies, and critical buckling loads of a multilayer curved beam. He proved that at higher values of the central opening angle, the obtained results using the Cartesian system are not accurate as of the polar system. The influence of thickness stretching was taken into account by Talebizadehsardari et al. [64] to analyze wave features of a curved graphene-reinforced beam. In another work, Arshid and Amir [65] used a similar system using higher-order shear deformation theory, considering the effect of graphene platelets reinforcement on the frequency response of a curvy microbeam.
The literature review has shown that the buckling response of microbeams is not well investigated with initial curvatures that are composed of three multilayers of functionally graded porous (FGP) core with two CNT-reinforced composite skins. The whole structure is affected by both temperature and moisture variations. To obtain reliable results of the central opening angle, a refined shear deformation theory (RSDT) as a higher-order theory capable of considering the thickness stretching effect is used in the polar coordinates to represent the displacements of the microstructure. Also, MCST that uses a length-scale parameter is used to investigate the scale influence. The governing equations are initiated using virtual displacements and variational approaches and then solved analytically. The results are presented based on various parameters' variations and discussed in detail. The novelty of the present work lies in considering the three-layered curved microbeam in polar coordinates with such a layer arrangement and, also, in the hygrothermal environment. The findings of this work will help in designing and creating more optimal engineering and structures, especially those that have curvature.

Basic relations 2.1 Model description
As shown in Figure 1, a three-layered microbeam with an initial curvature composed of two skins and a core is considered. The core of the microstructure is made of FGP materials that consist of pores; while the top and bottom layers are CNT-reinforced composites that have the role of integration of the whole structure. The whole structure is affected by both moisture and temperature variations and lay on an elastic base that is modeled by the Pasternak foundation type. The geometrical parameters of the microbeam, namely its central opening angle and curvature mid-radius, are given by β and R, respectively. Also, the thickness of each layer is h i (i = f for the skins; i = c for the core).

Constitutive law
As mentioned earlier, the core and the skins of the curved microbeam are composed of FGP and CNT-reinforced composite materials, respectively. The constitutive law for both the core and skins in the presence of hygrothermal effect can be written as follows [66]: The superscripts c and f denote the core and skins, respectively. α ij (r) and β ij (r) are the thermal and moisture expansion coefficients, respectively, which vary through the thickness of the layers, the temperature and moisture changes are given by ΔT and ΔH, while Q ijkl are the stiffness components that can be determined for the FGP core via the following relations [67]: ii c c c E c (r) is Young's elasticity modulus and ν c represents the core Poisson's ratio. As it is clear, Young's elasticity modulus depends on the thickness variation. Therefore, its variations obey different functions and depend on the patterns of pores' placement. If the pores are considered to be distributed symmetrically throughout the core's thickness (which is known as the "even" pattern in this article), the following relation may be used to express this distribution and show its effect on the elasticity modulus [68]: For "uneven" porosity pattern, the pores are considered as distributed asymmetrically concerning the core's midplane and in the "uniform" pattern, they are independent of thickness [52,69]:  ( ) , for uniform pattern, In the above relations, e 1 is the porosity and changes between 0 to 100%. This dimensionless parameter shows the void to the bulk volume and its increment shows the increase of the void volume.
In addition to thickness, the properties of the FGP core also depend on the temperature changes. It means that the core's properties vary with changes in temperature. To consider this effect, the following relation is suggested by the previous studies for the temperature dependency of the SUS304 FGP core [67]: where P c represents the core's properties and T is the temperature (K), and P i (i = 0, 1, 2) coefficients are dependent on the temperature. Furthermore, to describe stiffness components Q ij for the CNT-reinforced composite skins, the following relations can be used [70]: Similar to the FGP core, the skins' Young's modulus also depends on both the thickness and temperature. First, its dependency on thickness is considered. The following relations show the effective properties of CNTreinforced composite skins [57,67]: where the longitudinal and shear Young's moduli of CNTs and matrix are given by E 11 CNT , G 12 CNT , E m , and G m , respectively. The efficiency parameters of CNT that could be found via molecular dynamics simulations (MDS) are given by η i (i = 1, 3) and their values for different volume fractions of CNTsare presented in Table 1 [71].
Throughout the article, m refers to the matrix of the skins. The hygrothermal properties of the skins can be obtained via RM with the following relation [72]: Also, for Poisson's ratio, the following relation hold [73]: In the above relation, the volume fraction of CNTs is given by V CNT ⁎ , which can be determined by the following relation [74]: where w CNT represents the mass density of reinforcements and CNTs and the matrix densities are given by ρ CNT and ρ m . Five different patterns for the CNT distribution, namely UU, FG-VA, FG-XX, and FG-OO, are considered in this work. The following relations hold for each of the mentioned distribution patterns [66,75]: for FG X pattern, In the above relations, i = t, b where t and b denote top and bottom skins. Also, for the upper skin, a plus sign and for the lower skin, a minus sign is used. The mentioned patterns are shown in the cross-sectional view of the beam in Figure 2.
Next, the temperature dependency of the skins is considered. The skins are made from polymethyl methacrylate (PMMA), which is reinforced by CNTs. The polynomial function considered for CNTs reinforcement dependence to temperature is as follows [76]: where P 1 , P 2 , and P 3 are coefficients that represent the dependence of the CNT properties on the temperature.

Derivation of the equations
To analyze the aforementioned microstructure, for the first time, the polar coordinate system is used, especially for deeply curved beams. The origin of this polar coordinate system lies on the left corner side of the microstructure's midplane. Also, the displacement components of its arbitrary point are described based on a novel RSDT, which considers the thickness stretching effect that, again, results in more reliable values compared to the theories that neglect this effect. Therefore, the strain components of an arbitrary point using RSDT are demonstrated as [63] ( ) ( ) In the above relations, the displacements of the microplate's midlayer in axial and transverse directions are, respectively, given by u and w and ψ and φ are additional displacements. Also, f(r) is the shear deformation function that obeys the following relation based on RSDT: where h is the total thickness of the structure. The variational approach and the virtual displacement principle were applied [77]: In the above equation, U c and U m are, respectively, the classical and nonclassical strain energies and the work of externally applied loads is shown by W.
The total strain energy (i.e., U c + U m ) of the under consideration curved microbeam is given by U T and includes those of the core and skins; therefore [78], As noted previously, the MCST is used to demonstrate the scale influence on the critical buckling load of the microstructure. In equation (19), ( ) m ij s is the higher-order stress: In the above equations, λ (=E/2(1 + ν)) is Lame's parameter and δ ij is the Kronecker delta.
Also, μ m is the length-scale parameter relevant to the symmetric rotation gradients ( ) χ ij s [79]: This leads to obtaining the nonzero terms as To calculate variations of the total strain energy, the following relation is obtained: In this study, the external work includes three terms: the first term is related to the Pasternak foundation, the second one is attributed to the applied in-plane loads, and the third is due to the hygrothermal loads.
As expressed before, the microstructure was rested on an elastic foundation bearing both normal and shear loads using springs and shear layers, respectively. Consequently, the work of this foundation is applied to the structure's lower surface and can be determined by [80] where k s and k p are Pasternak foundation moduli that relate to the springs and shear layer, respectively. Also, h is defined as To obtain the work applied by the in-plane load, the following expression may be employed [81]: and the variation of the hygrothermal load work can be evaluated by equation (26) [82]: Differential governing equations can be achieved by inserting the derived terms for strain energy variations and the work applied by external loads in the virtual displacement relation, equating expressions for each coefficient to zero, which leads to the following equations: in which Rewriting the above classical and higher-order stress resultants as a function of displacements and substituting in the derived governing equations, yields where    It should be noted that for all of the mentioned stress resultants and coefficients, the integration interval for the bottom layer is from R − h c /2 − h f R − h c /2, for the core it is from R − h c /2 up to R + h c /2, and for the top layer it is from R + h c /2 up to R + h c /2 + h f .

Solution procedure
The displacements can be considered as functions satisfying the simply supported boundary conditions at the ends of the microbeam that are immovable. So, the following functions can be considered for displacement components [84,85]: where U, Ψ, W, and Φ are amplitudes of displacements of the microbeam, respectively. Substituting the proposed functions in the governing equations of (33)-(36) yields Here, [K e ] is stiffness, and [K g ] is the geometric stiffness matrixes. N cr is the critical buckling load that is attained by solving the eigenvalue problem (39).

Results and discussion
In the first step, to ensure the correctness of the results, they were compared with the literature. To achieve this aim, due to the lack of similar work, we had to compare the obtained results with a simpler structure to make them comparable. By neglecting some parameters such as the initial curvature, skins, pores, hygrothermal environment, and foundation, the results are achieved and compared to the published work in the literature. The achieved dimensionless critical buckling loads are listed in Table 2 for two different length to thickness ratios, which, according to this table, indicates that the current results and those of others [57,[86][87][88] coincide suggesting a good agreement. The dimensionless critical buckling load N ⁎ in Table 2 is calculated using / N L EI cr 2 in which L = 10 m, ν = 0.3, and E = 30 GPa. Noteworthily, the slight difference between the results can be due to differences in displacement field theories and, also, the different approaches to derive the differential equations. Therefore, the results and their accuracy are verified. In the following, the results regarding all of the parameters are presented. The hygro-thermomechanical properties of the FGP core and the nanocomposite skins are presented in Tables 3-5, which are dependent on the temperature changes, as stated previously [89,90]. Also, Poisson's ratio of CNTs is 0.175 and for their moisture coefficient β 11 = 3 × 10 −4 /K is assumed. Moreover, the length-scale parameter, i.e., μ m = 17.6 μm is based on previous studies, and the microbeam's height is ten times more than it, while that for the skins is 0.1h c . It must be noted that for all of the following graphs unless otherwise expressed, the beam's length and the central opening angle are 10 m and π/3. Also, the relation between the curvature radius and central angle is R = Lβ; 40% of porosity and 0.14 of CNTs' volume fraction are used to extract the results. Moreover, even patterns of pores' placement and uniform one for the CNT distribution are used. Also, the presented dimensionless critical buckling load ( ) N ⁎ obeys the following relation: where E NP represents Young's modulus of the core without any pores.
In the next phase, the results are presented and discussed in detail. Figure 3 depicts the influence of two factors. First, the porosity percentage, and the second the pores' placement patterns. From Figure 3 it is found that by increasing the porosity from zero up to 80%, the critical buckling load reduces between 27 and 45% and depends on the pores' placement pattern. As known, the critical buckling load relates to the stiffness of the structure and increasingporosity means a reduction in stiffness; accordingly, by reducing the stiffness, the critical buckling load also reduces. But, comparing different patterns of pores' placement shows that there is a meaningful difference between the even pattern and two others. It may arise from the distribution of the elasticity modulus that leads to stiffer surfaces in comparison to the midplane of the microstructure.
In Figure 4, in addition to reconsidering the porosity variations effect, the changes in the central opening angle are also taken into consideration. As it is clear, increasing this angle from 30 to 90 degrees has reduced the critical buckling load, but by enhancing it from 90 to 150 degrees, a sharp enhance in the critical buckling load is seen.  Property Value The reason for this effect backs to the geometrics of the structure that reduces or increases the stiffness of the whole structure. Different patterns of CNTs' distribution are compared in Figure 5. As can be understood, there is no significant difference between their three types including O-O, U-U, and X-X, but the critical buckling load for the V-A pattern is the highest values, and that for A-V patterns are the lowest. It should be noted although it is stated that there is no significant difference between the mentioned three types; however, with more focus on the results it is seen that among these types, X-X and O-O have minimum and maximum values of the critical buckling loads. Figure 5 also shows the effect of another geometrical parameter of the structure, i.e., thickness to length. As can be expected, increasing this parameter resulted in enhancing the critical buckling load because the stiffness of the microbeam enhances.     Next, the influence of the volume fraction of CNTs on the buckling response of the curved microbeam is considered. Regardless of the nanotubes distribution patterns, it is seen in Figure 6 that by enhancing the amount of CNTs reinforcement within the matrix, due to their extreme stiffness, the stiffness of the whole skins and, in the following, the stiffness of the whole microbeam enhances, which means an increase in the critical buckling load.
The effect of length-scale and temperature variations for two values of volume fraction (i.e., 0.11 and 0.17) of CNTs is illustrated in Figure 7a and b. An interesting outcome of this figure is that the influence of the length scale on the critical buckling load strongly depends on the temperature. As can be seen, at room temperature (i.e., at ΔT = 0), increasing the h c /μ m ratio leads the results to reduce but as the temperature increases, the structure behaves differently. It means increasing the h c /μ m results in enhancing the critical buckling load. It is true for both considered values for the volume fraction of CNTs. By comparing Figure 7a and b, it is found that the results of critical buckling load in Figure 7b, which relates to the CNTs' volume fraction of 0.17, are more than 0.11 (i.e., Figure 7a). Figure 8a and b depicts the various effects of moisture on the critical buckling load. Based on these figures, (a) 17     increasing the moisture results in a softer structure, which means the critical buckling load is reduced. Also, it is seen, unlike the temperature, that at all values of moisture difference (i.e., ΔH) increasing the h c /μ m ratio results in the decrease of the buckling load. A similar finding with previous figures about the effect of volume fraction of CNTs is seen in these figures too.
Tables 6 and 7 illustrate the effect of different parameters including thicknesses ratio, temperature, and moisture changes, as well as different pores' placement patterns on the corresponding critical buckling load. These tables confirm the previous findings that stated that at higher temperature and moisture values, the critical buckling load decreases. Also, the conclusions about the effect of pores' placement patterns on the results are confirmed again. Figure 9 compares three different cases, namely without foundation, just spring included foundation, and both springs and shear layer included foundation and investigates their effect on the critical buckling load of the under consideration model. As can be guessed, in the absence of the foundation, due to less structural rigidity, critical buckling load has the least value, but the addition of the foundation, either just springs or both springs and shear layer, results in increasing the critical buckling load. Since the presence of both springs and shear layer increases the rigidity of the foundation, therefore, its critical buckling load has more value too.
In Table 8, the exact values of dimensionless critical buckling load are listed for different values of the foundation's moduli. As can be found, enhancing both moduli (either springs or shear layer) produces a structure that is more rigid and it leads to enhancement in the critical buckling load.

Conclusion
This presented study focused on considering the influence of the hygrothermal environment on the buckling behavior of a deeply curvy microbeam that is composed of three layers, namely an FGP core and two CNT-reinforced skins. Using MCST for accounting scale influence and RSDT in polar coordinates, which involves the thickness stretching effect and produces more accurate results, the   governing differential equations are extracted. By ensuring the validity of the results, the effect of different parameters on the critical buckling load is investigated. Finally, some of the most prominent findings of this study are noted in the following. It is seen that the porosity percentage enhances, reduces the microbeam stiffness that leads to lower values of the critical buckling load. Also, the addition of CNTs to the matrix of the skins makes the skins and, accordingly, the whole structure stiffer, which means enhancing the critical buckling load. Among the three considered patterns of pores placement through the FGP core, the even pattern leads to the highest critical buckling loads. By comparing different types of distribution patterns of CNTs, it is seen that, respectively, V-A and A-V patterns lead to the highest and the lowest results of critical buckling load. Although increasing the central opening angle of the beam from 30 up to 90 degrees results in lower critical buckling load, but above 90 degrees, the critical buckling load tends to increase rapidly. Enhancing both temperature and moisture reduces the structure's stiffness, which finally leads to decreasing the critical buckling load. Increasing h c /μ m as a criterion of length-scale parameter generally leads to the decrease of the critical buckling load. But at a higher temperature, it behaves differently. The addition of springs and shear layer to the foundation increases the structure's rigidity, which causes increases in the critical buckling load. The author believes that considering such a structure in the polar coordinate, which results in more accurate results, can be utilized especially for deeply curved beams and the outcomes can be accounted as benchmarks for future works.
Funding information: The author states no funding involved.
Author contributions: The author has accepted responsibility for the entire content of this manuscript and approved its submission.