Thermal Buckling and Postbuckling Behaviors of Couple Stress and Surface Energy-Enriched FG-CNTR Nanobeams

: Small-sized structural elements such as beams, plates, and shells are usually used as nanomechanical resonators, nanoscale mass sensors, nanoelectromechanical actuators, and nanoen-ergy harvesters. At the nanoscale, the structures usually possess a high surface area-to-bulk volume ratio, leading to the free energy related to surface atoms becoming considerable compared to that of the bulk part. Earlier reports indicated several physical reasons for size-dependent phenomena, e.g., nonlocal stress, surface energy, and couple stress. To provide an in-depth insight into the mechanical behavior of small-scale structures, size-dependent continuum models including two or more physical factors have attracted the attention of the academic community. This research analyzes the thermal buckling and postbuckling characteristics of functionally graded carbon nanotube-reinforced (FG-CNTR) nanobeams with a tri-parameter, nonlinear elastic foundation and subjected to a uniform temperature rise. Chen-Yao’s surface energy theory and Yang’s symmetrical couple stress theory are combined to capture two types of size effects in nanobeams. The postbuckling model is formulated based on the Euler–Bernoulli deformation hypothesis and Euler–Lagrange equation. Using a two-step perturbation technique, the related postbuckling equilibrium path is determined. In numerical analysis, the impacts of surface energy, couple stress, elastic foundation, boundary conditions, geometric factor, layout type, and volume fraction of CNTs on the thermal buckling and postbuckling behaviors of nanobeams are revealed. It is indicated that considering couple stress or surface energy can lead to a signiﬁcant increase in the postbuckling stability of nanobeams compared to the case in which it is not considered. In addition, there is a reverse competition between couple stress or surface energy effects on the thermal buckling responses of nanobeams. As the temperature rise will cause the material elastic moduli softening, the thermal buckling load–deﬂection curves of nanobeams with the temperature-independent case are much higher than those with the temperature-dependent cases.


Introduction
Nanostructures possibly offer great potential designs and applications in a wide range of nanoprobes, optoelectronics, and biomedical implants [1,2]. There is also growing interest in utilizing nanostructures as reinforcements, e.g., carbon nanotubes (CNT), thanks to their excellent physical and chemical properties [3]. For instance, functionally graded carbon nanotube-reinforced (FG-CNTR) nanocomposites show many excellent performances (such as ultra-high elastic modulus, large, and thermal resistances), and hence have a number of mechanical applications in spaceflight, biomedicine, and sensors. The largest difference between CNTR-reinforced composites and traditional carbon fiber-reinforced composites is that the fiber volume content of CNTR composites can reach more than 60%, while the weight percentage of CNTs in CNTR composites is only 2~5%. At the micro/nanoscale, several experiments on the bending of beams [4][5][6][7][8] and torsion of copper wires [9][10][11][12] have observed size-dependent phenomena in elastic modulus and yield strength. To achieve a better design performance of nanodevices, it is crucial to deeply understand the size-dependent static and dynamic behaviors of nanostructures. However, classical continuum mechanics excludes size effects into constitutive formulation and is thus not adequate for the mechanical analysis of micro/nanostructures. Given such shortcomings, scholars have formulated different non-classical continuum theories, such as the nonlocal elasticity theory [13,14], couple stress theory [15][16][17][18], strain gradient theory [4,[19][20][21][22][23], and surface energy theory [24][25][26]. Among them, the symmetrical couple stress theory (CST) of Yang et al. [17] and the classical surface energy theory (SET) [24] of Gurtin and Murdoch have gained much attention in the field of micro/nanostructural mechanics.
The classical SET stems from the following physical fact: the atoms on the surface of the material lack the action of some atoms, and many suspended bonds would be generated; as a result, the forces on the atoms on the surface are different from those in the bulk materials, resulting in the reduction of symmetry and surface relaxation. The symmetrical CST starts from the perspective of introducing higher-order equilibrium relations. Due to the arising of the additional moment of the couple equilibrium equation, the couple stress tensor is forced to a be symmetric tensor; thus, the related constitutive equation contains one material length-scale parameter (MLSP) only. This feature reduces the difficulty of determining the MLSP in non-classical constitutive relation greatly. Recently, Thai et al. [27] and Roudbari et al. [28] provided two detailed overviews on the continuum mechanics modeling for micro/nanostructures within the framework of various sizedependent continuum theories (including the symmetrical CST and classical SET). Different geometries such as rods, beams, plates, and shells were considered.
To more comprehensively understand the vibration and buckling behaviors of micro/nanostructures, some scholars have carried out studies on the size-dependent continuum mechanics modeling of beams and plates with multiple size effects. For example, Gao and Mahmoud [29][30][31] developed couple stress and surface energy-enriched modes for the Euler-Bernoulli beam, Timoshenko beam, and Reddy beam, and gave the related Navier solutions for the simply supported case. Gao and Zhang [32,33] proposed the size-dependent Kirchhoff plate and Mindlin plate models resting on a two-parameter elastic foundation, combining the use of the symmetrical CST and classical SET. Based on the high-order continuity basis functions of non-uniform rational B-splines, Yin et al. [34] presented a new isogeometric Timoshenko beam model with couple stress and surface energy effects, fulfilling the higher-order continuity condition. Zhang et al. [35] proposed a size-dependent nonlinear beam model with structure-foundation interaction along with strain gradient and nonlocal effects, and to solve the deep postbuckling and nonlinear bending problems using a two-step perturbation method. Allahkarami et al. [36] carried out a vibration analysis of symmetrical CST-based agglomerated CNT curved beams. Thanh et al. [37] combined the isogeometric analysis and symmetrical CST to analyze the static and free vibration characteristics of an FG-CNTR nanoplate. Zhang et al. [38] developed a size-dependent Kirchhoff plate model that considered the effects of surface energy, the strain gradient, and inertia gradient on static bending and free vibration behaviors of thin microplates and constructed a C 2 -type differential quadrature plate element. In another work, Zhang et al. [39] presented surface energy-enriched gradient elastic Euler-Bernoulli, Timoshenko, and Reddy beam models and gave a novel numerical solution method. Khabaz et al. [40] studied the combined effects of strain gradient and surface energy on the dynamic behavior of an advanced sandwich composite microbeam including piezoelectric layers. Dangi et al. [41] proposed a theoretical model for bidirectional FG Euler-Bernoulli nanobeams with nonlocal stress, strain gradient, and surface energy effects and applied the model developed to evaluate three types of size effects on the natural frequencies of nanobeams. Attia and Shanab [42] made a combination of the symmetrical CST and classical SET to study the size-dependent geometric nonlinear behavior of FG Euler-Bernoulli and Timoshenko nanobeams. Shaat et al. [43] employed Newton's second law to establish the governing equation for surface energy and couple stress enriched Kirchhoff nanoplates. Lu et al. [44] developed three isotropic plate modes with nonlocal stress, strain gradient, and surface energy effects according to the Kirchhoff, Mindlin, and Reddy deformation assumptions, respectively. Duong et al. [45] explored the static responses and stress concentration phenomenon of FG-CNTR composite cylindrical shells with various boundary conditions using the higher-order, shear-normal deformation theory and Laplace transform. Doan et al. [46] examined the vibration response and static buckling of variable flexoelectric nanoplates using the FEM and novel hyperbolic sine shear deformation theory, in which the thickness is adjusted by linear and nonlinear rules. Thom et al. [47] proposed a phase-field model to investigate the thermal buckling of fractured FG plates based on the third-order shear deformation plate theory and demonstrated the difference between the plate's static stability response to temperature-dependent and temperature-independent cases. Do et al. [48] used the phase field model and Mindlin plate theory to predict the thermal buckling of cracked FG plates and considered two cases with and without the difference between the neutral surface and middle surface.
To provide an in-depth insight into the mechanical behavior of small-scale structures, size-dependent continuum models including two or more physical factors have attracted the attention of the academic community. However, through a thorough literature survey, one can see that little attention has been paid to the thermal buckling and postbuckling of FG-CNTR nanobeams considering two or more types of size effects and placed on the nonlinear elastic foundation. Thus, this paper intends to fill such a blank. This study is outlined as follows. Section 2 derives the theoretical formulation for the Euler-Bernoulli nanobeams embedded on a tri-parameter Winkler-Pasternak elastic foundation within the combined framework of symmetrical CST and Chen-Yao's SET. Section 3 adopts a two-step perturbation method to solve the related postbuckling equilibrium path of a nanobeam with two ends simply supported (SS) or two ends clamped-clamped (CC). Finally, selective numerical examples are presented to exhibit the influences of various factors on the thermal buckling and postbuckling responses of FG CNT-reinforced nanobeams.

Geometrical and Material Description
Consider an FG-CNTR nanobeam made of CNT agents and silicon (Si) matrix with length l, total thickness H b , and width b, as shown in Figure 1. The nanobeam is under the action of uniform temperature rise ∆T and is placed on tri-parameter elastic substrates. Moreover, the structure is referred to a Cartesian coordinate system (X, Y, Z). To capture the surface effect, the nanobeam is abstracted as a composite structure composed of an extremely thin surface layer and a bulk volume. Figure    For the UD pattern: For the X pattern: For the O pattern: where The subscripts "CNT" and "m" indicate the CNT reinforcements and matrix, respectively. ω and ρ stand for the total mass fraction and mass density, respectively.
The displacement field of the Euler-Bernoulli theory is written as where U indicates the axial displacement due to in-plane stretching, and W denotes the deflection in the Z direction. The nonlinear strain is

Governing Equations
According to the symmetrical CST [17] and Equation (5), the strain energy of the nanobeam is written as in which In Equations (6) and (7), the internal material length l is utilized to reflect the MCS effects; α (11) , E (11) , and G (12) , respectively, denote the thermal expansion coefficient, Young's modulus, and shear modulus of the resultant nanocomposites; ∆T = T − 300K denotes the temperature rise. The thermophysical parameters can be evaluated by a modified version of the mixture rule [50].
Based on [26], the surface energy of the nanobeam is given by with µ (1) and µ (2) where C nb denotes the perimeter of the cross section. The potential energy of a tri-parameter substrate is defined as where K 1 , K 2 and K 3 refer to the linear spring, shear spring, and nonlinear spring parameters, respectively. Thus, the total energy functional for a size-dependent FG-CNTR nanobeam is where L d is the following Lagrangian density: where n W is the unit normal vector along the Z-axis. In the prevailing circumstance, we have I S0 = 2b and I S2 = bH 2 b /2 + H 3 b /6. Using the Lagrange-Euler equation, the equilibrium equations of an FG-CNTR nanobeam are Substituting Equation (13) into Equations (15) and (16) and applying the immovable boundary conditions yield where , , Due to the introduction of couple stress and surface energy effects, Equation (17) is different from its classical counterpart.

Solution Methodology
To solve the nonlinear governing equation, the two-step perturbation technique [51] is used to determine the thermal buckling equilibrium path. This method has a wide application in nonlinear bending, buckling, and vibration of nanostructures [35,52,53]. According to [51], w and λ T can be assumed as where ε is a small perturbation parameter. Inserting Equation (23) in Equation (20) leads to the following equations: O ε 1 : O ε 2 : O ε 3 : It is clearly seen that the perturbation solutions of Equations (24), (25) and (26) depend on the related boundary conditions (see Equations (21) and (22)).

Immovable SS Ends
To satisfy the SS ends shown in Equation (21), the solutions of Equations (24), (25) and (26) are represented by where w i is the assumed ith-order perturbation solution. By solving these perturbation equations step by step with the help of the Galerkin procedure, asymptotic solutions for nonlinear buckling behaviors of an SS nanobeam can be determined as where In the present case, m and x are usually set to 1 and π/2, respectively. According to Equation (28), the maximum deflection w max becomes Substituting Equation (31) into Equation (29), the analytical postbuckling equilibrium path of the SS nanobeam is given as

Immovable CC Ends
The asymptotic solutions corresponding to the CC ends can be assumed as Similarly, substituting Equation (33) into Equations (24), (25) and (26), and upon simplification, yields with By introducing the quantities m = 1 and x = π/2, and by considering Equation (34), Applying Equation (37) to Equation (35) gives Currently, λ (i) CC in Equations (32) and (38) are related to the effective material properties of CNT-based reinforced nanocomposites.

Results and Discussion
Several numerical examples are presented to show the thermal buckling and postbuckling behaviors of FG-CNTR nanobeams. An FG nanobeam made of the silicon (Si) matrix and CNT agents is considered, where material properties are dependent on temperature.
Si matrix: CNT reinforcements: In what follows, the temperature-dependent (TD) case, belongs to the conditions when composite properties are calculated at current temperature, while the temperature independent (TID), on the other hand, refers to the case when properties are evaluated at 300 K. As presented in Equations (39), (40) and (41), the E m and G m of the Si matrix decrease as the temperature rises, while α m is just the opposite. For the CNT reinforcements in Equations (42), (43) and (44), the E (11) CNT and α (11) CNT decrease as the temperature rises. Then, the critical buckling load of the TID case is higher than that of the TD one. For the TD case, the thermal buckling load should be obtained for each temperature load step by the incremental method. In the first step, we calculate the material property parameters at the reference temperature T 0 = 300K and obtain the trial solution of the thermal buckling load ∆T 0 according to Equations (32) or (38). In the second step, we calculate the material property parameters at T 1 = T 0 + ∆T 0 and obtain the trial solution of the thermal buckling load ∆T 1 according to Equations (32) or (38). Then, follow the above steps until |(∆T i − ∆T i−1 )/∆T i−1 | < 0.1%.

Verification Study
The current model can degenerate into the size-independent one by setting S = 0 N/m and l = 0 nm. Here, the isotropic macrobeam made of metal/ceramic phases is selected, and effective material properties can be found in Ref. [54]. Tables 1 and 2 compare the results of the critical buckling temperature quantity ∆T cr α 0 (L/H b ) 2 for a macrobeam under various boundary conditions with results in previous works. Table 1 shows that the present results are larger than the available ones in [54]. The differences between the two are attributed to the fact that the existing theoretical model is based on the Reddy's higher-order shear deformation theory, while our model is based on the Euler-Bernoulli beam theory. The introduction of shear deformation effect will weaken the structural bending rigidity. Moreover, the bending moment is assumed as a function of coordinates in Ref. [54], which is actually unreasonable. In Table 2, one can observe that the present results and the reported ones in Ref. [55] have a good consistency. To provide a consistent comparison with previous works, Figure 3 presents the thermal postbuckling equilibrium path of a CC macrobeam with TD and TID cases, where ρ = A Z 2 dA /A denotes the radius of gyration. As expected, the results in this paper are consistent with those in the literature. (α 0 = 12.33 × 10 −6 (1/K)).

Effects of Surface Energy and Couple Stress
In this section, parametric studies are presented to the show surface energy and couple stress effects on the postbuckling behaviors of FG-CNTR nanobeams under different boundary conditions (BC). For simplicity, the following analysis uses the abbreviations "SE" and "CS", respectively, to refer to the case of only including surface energy (l = 0) and the case of only including couple stress (µ  Tables 3 and 4 tabulate the critical thermal buckling temperature ∆T cr for six groups of SE and CS parameters, respectively. Here, only the UD nanobeam with SS ends is selected. These results confirm that both SE and CS have pronounced effects on the critical buckling temperature ∆T cr . For smaller SE and CS parameters, 0.01µ (1) S and 0.01µ (2) S (see Table 3) or l/H b = 0.01 (see Table 4), the buckling load approaches the classical one (i.e., S = 0 and l/H b = 0). Generally, the stronger the size effects, the larger the buckling temperature ∆T cr . Furthermore, the buckling temperature ∆T cr of a nanobeam with a TID case is much higher than the nanobeam with TD case; hence, the TID condition may be overestimated with regards to the critical thermal buckling responses.    Figure 4 shows the SE and CE effects (i.e., µ S , and l/H b ) on the postbuckling behaviors of nanobeams. As the SE and CS parameters tend to zero, e.g., 0.01µ (1) S , 0.01µ (2) S , and l/H b = 0.01, the postbuckling response approaches the classical one. By increasing the SE and CS parameters, the critical buckling temperature load and the stability of the UD nanobeam increase, which implies the larger values of SE and CE parameters generate larger additional bending rigidity. According to Tables 3 and 4 and Figure 4, it can be concluded that an increase in µ (i) S (i = 1, 2) or l/H b postpones the branching point of the nanobeam. Additionally, for the uniform temperature field, the thermal postbuckling load-deflection curve of UD nanobeams is of the bifurcation buckling type. Meanwhile, the difference between the TID and TD cases becomes higher at the larger SE and CS parameters. The reason is that increasing the SE and CS parameters would lead to an increase of structural bending rigidity, especially for the TID case.  Figure 5 exhibits the SE and CE effects on the critical buckling load of a UD nanobeam with l/H b = 50. Clearly, the SE and CS make positive contributions to increasing the buckling load ∆T cr , which means small-scale effects may exert a considerable rigidity and make nanobeams harder to deform. These contributions are more prominent, especially when the l/H b is large. However, the TD case exhibits a negative impact on the buckling temperature when compared with those of the TID case, and this phenomenon is more remarkable when considering the SE effect. As predicted from Figure 6 for the SE, CS, and SE-CS analyses, the postbuckling loaddeflection curves are sensitive for the CS effects. Obviously, the biggest effect of the scale parameter is associated with the SE-CS response. Moreover, an interesting result observed from Figures 5 and 6 is that the SE is more dominant than CS for the present nanostructures. Thus, although the surface of a nanobeam system is modeled as an infinitely thin thickness (zero-thickness) layer, the characteristics of this layer could also significantly influence the buckling responses.

Effects of Geometrical Property and Substrate Stiffness
In this section, the impacts of the aspect ratio (l/H b ) and substrate stiffness constants (k 1 , k 2 , and k 3 ) on the buckling analysis of FG-O nanobeams with V * CNT = 0.12 are discussed in detail. The geometrical properties adopted are: L = 1000 nm, b = 40 nm. The combined effects of SE and CS and effective TD material properties are also involved in the present size-dependent model, and the same below.
In Figure 7, the effects of aspect ratio l/H b on the buckling load ∆T cr are depicted. It is again observed that the maximum buckling temperature is associated with SE-CS analysis, while the minimum is associated with classical elasticity analysis. Additionally, for classical results, the buckling temperature decrease with increasing the ratio l/H b , which is due to the reduced total stiffness of the structure. However, as predicted from Figure 7a for the size-dependent nanobeam with SS ends, the critical buckling temperature enhances as the ratio l/H b increases, and as the ratio l/H b increases, the SE and CE effects gradually become significant.  Figure 8 shows that, in the beginning, the paths of the load-deflection stage predicated by the size-dependent model are higher than that of the corresponding classical model, especially for large aspect ratios. Moreover, in the deep postbuckling region, the load-deflection curves calculated by both models have similar distribution patterns. For instance, the postbuckling response is the greatest for the nanobeam at l/H b = 60, followed by l/H b = 80 and l/H b = 100 third. The postbuckling load-deflection curves with different substrate models are plotted in Figure 9. An FG-O curved nanobeam with l/H b = 100 is selected. To capture the reactions acting on nanobeams from the foundation, the concepts of three substrate models are introduced. The elastic substrates are characterized by (k 1 , k 2 , k 3 ) = (100, 0, 0) for the Winkler substrate (W-S), (k 1 , k 2 , k 3 ) = (100, 10, 0) for the Pasternak substrate (P-S), and (k 1 , k 2 , k 3 ) = (100, 10, 1000) for a tri-parameter nonlinear substrate (T-S). Evidently, (k 1 , k 2 , k 3 ) = (0, 0, 0) stands for no substrate (N-S). It is noticed that the nonlinear stiffness (e.g., k 3 = 1000) has no contribution on the buckling load ∆T cr ; however, its effect is a dominant factor for the postbuckling analysis. Generally, the increasing substrate stiffness constants k 1 , k 2 , and k 3 lead to a lower deflection in the postbuckling behavior due to enhancement of the reaction force of the nanobeam from substrates.  Figure 10 also confirm that the structure with the FG-X pattern has the biggest postbuckling temperature, followed by the UD and FG-O types. Briefly, the FG-X nanobeam is the optimum pattern in this case.  Figure 11 shows the postbuckling behavior of UD and FG-X nanobeams with different CNT volume fractions. Obviously, the CNT volume percent V * CNT (= 0.12, 0.17, and 0.28) plays a key role in the buckling and postbuckling responses. The FG-X nanobeam with V * CNT = 0.12 has the highest critical buckling temperature and thermal postbuckling strength. Furthermore, from Figures 10 and 11, it can be drawn that the buckling load as well as postbuckling strength of FG-X nanobeam are higher than those of the UD nanobeam, while they slightly affect the responses when considering scale effects or decreasing V * CNT .

Concluding Remarks
According to Chen-Yao's SET and Yang's symmetrical CST, the thermal buckling and postbuckling analyses of FG-CNTR nanobeams under uniform temperature rise were con-ducted. All the thermomechanical material characteristics for CNT-based nanobeams are position-and temperature-dependent. The impacts of different parameters, particularly two types of size effects, on thermal buckling and post buckling were investigated via a two-step perturbation method. Major conclusions can be summarized as follows: (1) As the two types of size effects increase, the thermal buckling and postbuckling responses of nanobeams decrease. This indicates that the SE and CS effects can enhance the structural bending rigidity. Moreover, there is a reverse competition between the SE and CS effects on the thermal buckling responses of nanobeams. (2) As the slenderness ratio increases, the thermal buckling responses of the nanobeam increase. The larger the slenderness ratio is, the more significant the two types of size effects, which are quite different from the classical ones. The CNT distribution patterns and the volume percent play a key role in thermal buckling responses. (3) The TID's thermal buckling load-deflection curves of nanobeams are much higher than the TD's, especially for larger surface energy effects. The reason is that the increase of temperature will cause the material to soften. Therefore, further efforts should be devoted to developing nonlinear, finite element models of size-dependent theories, especially the symmetrical CST and SET-based models for practical applications.