Geometric Non-Linear Analysis of Auxetic Hybrid Laminated Beams Containing CNT Reinforced Composite Materials

In the current work, a novel hybrid laminate with negative Poisson’s ratio (NPR) is developed by considering auxetic laminate which is composed of carbon nanotube-reinforced composite (CNTRC) and fiber-reinforced composite (FRC) materials. The maximum magnitude of out-of-plane NPR is identified in the case of (20 F/20 C/−20 C/20 C) S laminate as well. Meanwhile, a method for the geometric non-linear analysis of hybrid laminated beam with NPR including the non-linear bending, free, and forced vibrations is proposed. The beam deformation is modeled by combining higher-order shear-deformation theory (HSDT) and large deflection theory. Based on a two-step perturbation approach, the asymptotic solutions of the governing equations are obtained to capture the linear and non-linear frequencies and load-deflection curves. Moreover, a two-step perturbation methodology in conjunction with fourth-order Runge–Kutta method is employed to solve the forced-vibration problem. Several key factors, such as CNT distribution, variations in the elastic foundation, and thermal stress, are considered in the exhaustive analysis. Theoretical results for some particular cases are given to examine the geometric non-linearity behavior of hybrid beam with NPR as well as positive Poisson’s ratio (PPR).


Introduction
Materials and structures with negative Poisson's ratio (NPR) behave in a counter-intuitive manner: when compressed (stretched) in the axial direction, they contract (expand) transversely. The materials and structures that exhibit this feature are also termed as "auxetics" [1]. Lakes [2] first reported NPR behavior in polyurethane (PU) foam with re-entrant structures. Wojciechowski [3] presented the first thermodynamically stable molecular model to study the mechanisms for generating auxetic behavior of solid. NPR was analyzed in systems containing rigid rotating hexamers [4]. Thereafter, analysis on cellular auxetics [5,6], multi-material auxetics [7,8], and auxetic composites [9,10] was carried out. A more systematic review of the development and application for the investigation of auxetic material and structures was reported by Ren et al. [11] and Lakes [12].
In recent years, laminated beams or plates with NPR are developed with more applications as primary structural elements in many fields. Fascinating properties of composite materials with NPR have led researchers to find outstanding applications in the field of aerospace [13,14], automobile industry [15,16], and civil engineering [17,18]. Recently, studies on the design of auxetic structures have propelled to achieve improvements in the mechanical properties such as impact resistance [19,20] and energy absorption [21,22]. On the other hand, with the growing applications of auxetic materials in different industries, developing multi-scale models for design and capturing the non-linear static and dynamic responses of the auxetic laminates under various loading scenarios is critical, when related to the structural design. A series of work has been dedicated to the study of composite beams, in particular, investigations have been carried out on their static and dynamic behavior [23][24][25][26]. However, the value of NPRs (ranging from −0.2 to −1) is not for any particular composite material in real engineering, but is a hypothetical value in a model analysis.
Fiber-reinforced composite (FRC) materials have attracted the attention of researchers, primarily due to the strong anisotropy they offer when designing laminate with NPR. Zhang et al. [27] showed that both the particular stacking sequence and the individual ply material (strongly anisotropic) are essential for a laminate to exhibit NPR. The authors also presented an optimal angle for ply and particular stacking sequences were presented. Evans et al. [28] specially designed a software to predict the effective engineering constants. It was reported that the NPR property can be obtained by designing stacking sequences in the laminated plates. Lempriere [29] measured that the effective Poisson's ratio (EPR) in orthotropic materials is −0.4, occurring at θ = 45 • orientation. Clarke et al. [30] reported that the EPR of the laminates show negative values for lay angle in a range between 15 • and 30 • for a (±θ). It should be noted that when layers are oriented at angles equal in magnitude but with opposite signs, the appropriate + or -sign is used. Herakovich [31] investigated the auxetic characteristics of laminated structure made of graphite-epoxy to determine the value of v e 13 . Such laminates exhibited a very wide range of NPR ranging from a peak of 0.49 for a laminate with ply angle of 90 • to as low as −0.21 for laminate with ply angle of (±25) S . Hine et al. [32] reported that the out-of-plane Poisson's ratio reached −1/2 when a high modulus of elasticity carbon fiber is used in the laminates. Matsuda et al. [33] observed that peak values of NPR in carbon fiber-reinforced plastic laminates were around −0.7 when the axis is oriented at 25 o . The influence of Young's modulus ratio (E 1 /E 2 ), the type of resin and the volume of fraction on the ERP (v e 13 ) of an angle-ply [±θ] 2s plate was investigated by Harkati et al. [34,35]. It was shown that the NPR of Kevlar and carbon reinforced composite plate is −0.746 at θ = 20 • . Therefore, it can be concluded form the above discussion that the maximum value of NPR of the laminates greatly depends on the ply orientations and stacking sequences in laminates [36][37][38]. One of the practical applications of such composite structures with NPR is the low-velocity impact resistance [39,40] and static indentation resistance [41]. Alderson and Coenen [39] found that the auxetic FRC laminates showed increased load and energy absorption to failure at low levels of impact energy. Zhou et al. [40] carried out experimental studies on the low-velocity impact response of 3D auxetic composites. They found that 3D auxetic textile composite exhibited excellent impact protective performance compared to 3D non-auxetic textile composite. Coenen and Alderson [41] manufactured auxetic laminates and evaluated their static indentation resistance in comparison with two laminates having near zero and larger positive Poisson's ratio. The results showed enhancement in load sustained and energy absorption by auxetic laminates.
Carbon nanotubes (CNTs) are extensively used in different fields of industry and research with the development in fabrication technology [42][43][44]. A composite laminate structure made from CNTRCs is prone to exhibit an auxetic feature due to inherent special properties such as strong anisotropy. Furthermore, functionally graded (FG) materials are frequently employed in engineering applications and the laminae with various CNT volume fractions are used to achieve excellent mechanical properties [45]. Based on the studies mentioned above, the auxetic concept of FRC laminates is introduced for designing FRC/CNTRC hybrid laminate with NPR. This hybrid laminate is modeled by placing FRCs in the outermost layers and FG-CNTRC core. Shen was the first investigator to analyze the particular characteristics and behavior of the FG-CNTRC structures at different scales [45]. Following this landmark research, many studies on the forced vibration of Materials 2020, 13, 3718 3 of 23 FG-CNTRC beams [46,47] and plates [48] were carried out. The free vibration of FG-CNTRC plate were examined by Huu Quoc et al. [49] with a new refined plate theory. Based on Timoshenko beam theory (TBT), the non-linear vibration, forced vibration, and bending behavior of FG-CNTRC beams were performed by Mohammadimehr et al. [50] and Ansari et al. [51]. Considering the matrix cracks, Fan and Wang [52,53] presented an investigation on the non-linear static and dynamic responses of hybrid laminated structures made either from FRCs or CNTRCs. Based on an element-free numerical approach in conjunction with self-consistent model, the effect of matrix cracks on bending and vibration characteristics of hybrid laminated plates were reported by Lei et al. [54,55]. In addition to the applications of FRC/CNTRC hybrid beam or plate, the adoptions of the FRC/CNTRC hybrid concept were later extended to the design of composite blades by Zhang et al. [56,57]. They carried out several investigations exploring the non-linear vibration characteristics of hybrid composited blades. In addition, the design and analysis of the hybrid laminates comprising multi-materials with different strengths have been carried out by [58][59][60][61].
From the literature survey, it can be noted that no research has been carried out to analyze of hybrid beams using FRC/CNTRC materials with NPR at different external conditions. Therefore, the primary objective of the current work is to consider both hybrid configuration and NPR and analyze the static and dynamic characteristics of these auxetic structures. Furthermore, a non-linear model is developed based on the higher-order shear-deformation theory (HSDT) and the Von Kármán large deflection assumption. The motion equations due to thermal stress and reaction force due to foundation are derived. The effects of the CNT material distribution, environment condition, foundation type on the free and forced-vibration feature, and non-linear bending behavior of hybrid laminated beams are investigated in detail.

Design of Hybrid Laminate with NPR
As we know that both the ply orientations and stacking sequence have a significant effect on the mechanical response of laminated structures. Therefore, the effective engineering constants are usually used for the convenience of engineers in describing the mechanical behavior of the laminates. Sun et al. [36] and Chen et al. [62] presented a model for EPR for general thick laminates. However, only the extensional response was taken into consideration while the bending and bending-extension coupling characteristics of the laminates were neglected. Therefore, their model fails to provide accurate solutions for EPR of an asymmetric angle-ply laminated plate. Considering the effects of bending and bending-extension coupling, the general solutions of the EPRs for an arbitrary angle-ply laminate are derived in the previous works [63][64][65].
For a laminate where the material parameters of the layers are not distributed symmetrically along the section, the general formula of v e 13 can be expressed as follows: where the matrixes (A 11 , A 13 , B 11 , B 5-1 , B 5-3 , D) are defined in [63][64][65]. For symmetric laminates, the bending-extension coupling stiffnesses B ij (i,j = 1,2 . . . 6) are zero. The preceding expression simplifies as follows: where A ij represents the beam stiffnesses, which are defined in terms of the transformed elastic coefficients (C ij ) k as: Materials 2020, 13, 3718 4 of 23 where the compliance constants (S ij ) of a laminate whose fibers direction makes an angle θ with the direction of X-axis (see Figure 1) and c = cosθ, s = sinθ.
where the compliance constants of laminate S ij are given as follow: where the basic material parameters of each layer of laminate are introduced as follow by referring to Figure 1. 13,23) = Shear moduli in i-j planes, respectively v ij , ij = 12, 13, 23) = Poisson's ratios (the subscripts i and j represent the loading and strain directions, respectively). The theoretical solution presented above can predict the EPR of an arbitrary angle-ply laminates. A systematic investigation of the symmetric hybrid laminates has been carried out. The hybrid laminated beam is designed by placing FRC in the outermost layers and CNTRC in the rest of the layers which gives excellent performance. The lay-ups of the hybrid laminated beam are considered eight-layered and angle-ply (θ F /θ C /-θ C /θ C ) S . The subscript/index "S" indicates that the laminate is symmetric and θ is the angle of CNT or fiber orientation. Unless otherwise stated, the superscript F represents the layer of FRC while the superscript C represents the other layers with CNTRC. These configurations are for an identical material and have constant thickness of 0.125 mm and 0.25 mm for FRC and CNT layers, respectively. Several types of CNT distributions are taken into consideration and the type of CNT volume fraction (V CN ) is given. The temperature-related material properties of CNTRCs are predicted by the extended micromechanical model [45] as summarized in Table 1. It should be noted that the thermal expansion coefficients (α 11 , α 22 ) are calculated by Equation (15), which is obtained by the Schapery model [66].
The EPR (v e 13 ) of the hybrid FRC/CNTRC laminated beam is obtained by adjusting specific ply orientations and stacking sequence. The results in terms of NPR and PPR of hybrid laminated beam and the corresponding ply orientations are presented in Table 3. To assess how the distribution and volume fractions of CNT influence the static and dynamic behavior of the beam, we considered five configurations summarized in Table 2. The volume fraction of fiber is fixed as Vf = 0.6. It differ in the following characteristics: FG-Λ: 0.6/(0.11)2/(0.14)2/(0.17)2/0.6, FG-V: 0.6/(0.17)2/(0.14)2/(0.11)2/0.6, FG-X:(0.6/0.11/0.14/0.17)S, FG-O: (0.6/0.17/0.14/0.11)S, and uniform distribution (UD): (0.6/0.14/0.14/0.14)S. It can be assumed that G13 = G23 = G12. The EPR ( ) of the hybrid FRC/CNTRC laminated beam is obtained by adjusting specific ply orientations and stacking sequence. The results in terms of NPR and PPR of hybrid laminated beam and the corresponding ply orientations are presented in Table 3.

Theoretical Modeling of Hybrid Beams
Laminates consist of layers of composites reinforced with CNTRC and FRC. Consider a hybrid laminated beam composed of eight layers with lamination scheme (θ F /θ C /-θ C /θ C ) S resting on a continuous visco-Pasternak foundation as shown in Figure 1. Figure 1a defines the coordinate system used in development of the hybrid laminated beam analysis. The XYZ coordinate system is assumed to have its origin on the middle face of the beam so that the middle surface lies in the XY-plane. Q(X,t) is the out-of-plane static or dynamic load. The displacement at a point on the X, Y, and Z directions are U, V, and W, respectively.
continuous visco-Pasternak foundation as shown in Figure 1. Figure 1a defines the coordinate system used in development of the hybrid laminated beam analysis. The XYZ coordinate system is assumed to have its origin on the middle face of the beam so that the middle surface lies in the XY-plane. Q(X,t) is the out-of-plane static or dynamic load. The displacement at a point on the X, Y, and Z directions are U , V , and W , respectively. The simply supported beam is resting on a three-parameter foundation including the Winkler foundation ( ), shearing layer stiffness ( ), and damping parameter ̅ . The reaction force from the visco-Pasternak foundation P0 (X, ̅ ) is given by: The method of analysis is based on the HSDT [67] for the laminated beam undergoing large deflection. The effect of the elevated temperature is considered by introducing thermal stress resultants , , and as shown in Appendix A. The motion equations are given as follow: ( ) The simply supported beam is resting on a three-parameter foundation including the Winkler foundation (K 1 ), shearing layer stiffness (K 2 ), and damping parameter C d . The reaction force from the visco-Pasternak foundation P 0 (X,t) is given by: The method of analysis is based on the HSDT [67] for the laminated beam undergoing large deflection. The effect of the elevated temperature is considered by introducing thermal stress resultants where Ψ x denotes rotation about the longitudinal axes. The reduced stiffnesses of the beam (A 11 ,B 11 , E 11 ) and the coefficients S ij and inertias I i are defined and given in Appendix A.
The formulae for the forces owing to thermal stress are given as: where ∆T is the temperature increment from a reference state (T 0 = 300 K), ∆T = T − T 0 . T is set as 400 K or 500 K. The coefficient A x is known in terms of the thermal expansion coefficients (α 11 , α 22 ): where α ij (ij = 11,22) are the thermal expansion coefficients and subscripts 11 and 22 denote the longitudinal and transverse directions, respectively. The formulae for calculating the α 11 and α 22 are obtained from [66].
It should be noted that Equation (11) represent restricted boundary conditions. It is suitable for a beam with immovable boundary condition i.e., the longitudinal displacement is equal to zero at both ends of the beam. It should be noted that this immovable boundary condition is unacceptable for the compressive post-buckling analysis of beam.
The non-linear motion equations for the vibration can be solved by a two-step perturbation approach proposed by Shen [68,69]. By deriving the dimensionless forms of the motion Equations (9)-(11), we get: The non-dimensional parameters mentioned above can be written as follow: (19) in which E 0 = E m and ρ 0 = ρ m are the elastic modulus and density of the matrix. A T x , D T x , F T x are the functions of thickness and A x can be defined by:

Free Vibration Analysis
For non-linear vibration problem, the solutions for Equations (16) and (17) consist of an additional displacement term and initial displacement term as a result of the varying temperature. The initial deflection under the thermal loading for these structures are reported from Shen [63]. By considering, τ = εt, the solution equations can be expanded as a function with a small perturbation parameter mboxemphε j (j = 1,2,3, . . . .) as given below: The boundary conditions for τ = 0 can be expressed as: If we substitute τ = εt and Equation (21) into Equations (16) and (17) and collect all the terms of the same order of ε i (i = 1,2,3, . . . ), we will get a set of perturbation equations. The solution of Euler-Bernoulli beam with simply supported conditions can be used for perturbation equations with ε i (i = 1) Then, the equations with ε j (j = 1,2,3,...) can be solved step by step, we will get: Introducing expression τ = t and εA (1) 10 =W m in Equation (26) and applying Galerkin procedure yielded Equation (26) which can be re-written as: The numerical solution of Equation (27) corresponds to the forced-vibration response of the structure. It can be converted to a free vibration problem by setting the uniform load to zero. Therefore, the solution of the differential equation with g c = λ q (t) = 0 can be obtained as follow: where ω NL and ω L represent the dimensionless non-linear and linear frequency, respectively. A = W m = W m /L. W m is the amplitude of deflection. According to Equation (19), the corresponding linear frequency can be expressed as Ω = ω L (π/L)(E 0 /ρ 0 ) 1/2 . The details of g 30 , g 31 , g 32 , and g 33 are described in Appendix B.
For the purpose of verification, Table 4 and Figure 2 show the comparison of the solutions obtained by the present method and the method proposed by Fan & Wang [70] with the present solutions of a hybrid laminated beam. The reinforced material of the hybrid beam consists of two  The material properties for the same matrix are E m = 2.5 GPa, ρ m = 1150 kg/m 3 , and ν m =0.34. The value of geometric parameters used are: L = 20 h, h = 0.5mm. Table 4 compares the fundamental frequencies Ω = Ω(L 2 /h) 0 / 0 for (θ C /90 F )S hybrid beam with different values of VCN and the angle of the CNT orientation with θ C = 15°,30°, and 45°. It can be found that the predicted frequencies tally reasonably well with those results obtained from Fan & Wang [70]. Meanwhile, the non-linear vibration of a (0 C /90 F )S hybrid laminated beam is validated with Fan &Wang study, as plotted in Figure 2. It can be observed from Figure 2 that the results in terms of frequency ratio obtained from the presented model are very close to the predictions made by Fan & Wang [70]. Based on the above discussion, it can be concluded that the proposed model is accurate to analyze free vibration of the hybrid laminated beams.  In the following, a detailed investigation of the vibration of hybrid laminated beams has been carried out. The effect of temperature T, distribution pattern, and the foundation constants (k1, k2) on the free vibration of hybrid laminated beams are scrutinized. In the following study, the angle-ply (25 F /20 C /-20 C /20 C ) S and (25 F /90 C /-90 C /90 C ) S hybrid laminated beams are adopted.
The effect of temperature T on the first four frequencies Ω (i = 1,2,3,4) of hybrid laminated beams with five different distribution pattern is investigated in Table 5. It is observed that increase in temperature cause decrease in the frequencies of both the angle-ply hybrid laminated beams. In addition, the temperature effect of the first-order frequency is more significant than on the other order frequencies. Consequently, the influence of degradation caused by thermal stress should be taken into consideration in the design and service of hybrid laminated beams.  Table 4 compares the fundamental frequencies Ω = Ω(L 2 /h) ρ 0 /E 0 for (θ C /90 F ) S hybrid beam with different values of V CN and the angle of the CNT orientation with θ C = 15 • , 30 • , and 45 • . It can be found that the predicted frequencies tally reasonably well with those results obtained from Fan & Wang [70]. Meanwhile, the non-linear vibration of a (0 C /90 F ) S hybrid laminated beam is validated with Fan &Wang study, as plotted in Figure 2. It can be observed from Figure 2 that the results in terms of frequency ratio obtained from the presented model are very close to the predictions made by Fan & Wang [70]. Based on the above discussion, it can be concluded that the proposed model is accurate to analyze free vibration of the hybrid laminated beams.
In the following, a detailed investigation of the vibration of hybrid laminated beams has been carried out. The effect of temperature T, distribution pattern, and the foundation constants (k 1 , k 2 ) on the free vibration of hybrid laminated beams are scrutinized. In the following study, the angle-ply (25 F /20 C /-20 C /20 C ) S and (25 F /90 C /-90 C /90 C ) S hybrid laminated beams are adopted.
The effect of temperature T on the first four frequencies Ω i (i = 1, 2, 3, 4) of hybrid laminated beams with five different distribution pattern is investigated in Table 5. It is observed that increase in temperature cause decrease in the frequencies of both the angle-ply hybrid laminated beams. In addition, the temperature effect of the first-order frequency is more significant than on the other order frequencies. Consequently, the influence of degradation caused by thermal stress should be taken into consideration in the design and service of hybrid laminated beams. Table 6 shows the effects of foundation constants on the hybrid laminated beam with length/thickness ratio of 20. Three different foundation coefficients (k 1 , k 2 ) are (10 3 , 10), (10 3 , 0), and (0, 0), respectively. It should be noted that structures with (k 1 , k 2 ) = (0, 0) are selected for comparison. As can be seen from Table 6, the natural frequencies of the hybrid laminated beam for higher (k 1 , k 2 ) are more than the other cases.  The frequency ratio-deflection curves for all the five types of hybrid laminated beams are described in Figure 1. As expected, FG-X has the highest frequency ratio while FG-O has the lowest frequency ratio. At the room temperature, FG-Λ and FG-V show similar higher amplitude frequency (see Figure 3). (see Figure 3).
The influence of temperature field on the non-linear vibration behavior of UD and FG-X hybrid laminated beam with NPR is examined by considering three different temperature values (T = 300 K,400 K,500 K). Due to the presence of thermal stress, increase in the temperature leads to a higher frequency ratio, as shown in Figure 4. It is found that (25 F /20 C /−20 C /20 C ) S laminated beams are more sensitive to the temperature effect than those with (25 F /90 C −90 C /90 C ) S.
The effect of the foundation coefficients on the resulting lager amplitude vibration frequency of UD and FG-X are examined in Figure 5, for a member with L = 10 h. From Figure 5, it is confirmed that the beam becomes more rigid as the foundation coefficients of the beam increases. Hence, the frequency ratio of the UD and FG-X increases.  The influence of temperature field on the non-linear vibration behavior of UD and FG-X hybrid laminated beam with NPR is examined by considering three different temperature values (T = 300 K, 400 K, 500 K). Due to the presence of thermal stress, increase in the temperature leads to a higher frequency ratio, as shown in Figure 4. It is found that (25 F /20 C /−20 C /20 C ) S laminated beams are more sensitive to the temperature effect than those with (25 F /90 C −90 C /90 C ) S .  The effect of the foundation coefficients on the resulting lager amplitude vibration frequency of UD and FG-X are examined in Figure 5, for a member with L = 10 h. From Figure 5, it is confirmed that the beam becomes more rigid as the foundation coefficients of the beam increases. Hence, the frequency ratio of the UD and FG-X increases.

Forced-Vibration Analysis
In this section, the dynamic response of hybrid laminated beams is investigated by considering uniform load Q which is dependent on time as plotted in Figure 1. We need to determine the relationship between central deflection and time. Hence, a numerical procedure for solving the second order differential Equation (27) is employed here. Given the initial value Wm (t0) and (t0) at initial time t0 = 0, Equation (27) can be solved to obtain the central deflection-time relationship for the beam under applied load by employing the fourth-order Runge-Kutta method. Here, it should also be noted that the initial deflection for FG-Λ or FG-V will be triggered because of thermal stress.
In the following, the time response history of the verification analyses is presented. Figure

Forced-Vibration Analysis
In this section, the dynamic response of hybrid laminated beams is investigated by considering uniform load Q which is dependent on time as plotted in Figure 1. We need to determine the relationship between central deflection and time. Hence, a numerical procedure for solving the second order differential Equation (27) is employed here. Given the initial value W m (t 0 ) and . W m (t 0 ) at initial time t 0 = 0, Equation (27) can be solved to obtain the central deflection-time relationship for the beam under applied load by employing the fourth-order Runge-Kutta method. Here, it should also be noted that the initial deflection for FG-Λ or FG-V will be triggered because of thermal stress.
In the following, the time response history of the verification analyses is presented. Figure 6 shows the comparison of the solution of a cross-ply (0/90/0/90/0) S graphene-reinforced composite (GRC) laminated beam under a sudden load Q = 12 Mpa by the method presented by Fan et al. [71]. and the method presented in this study. The dimensions and material properties of the GRC layer are as follows: L = 20 h = 60 mm, E G 11 = 1.812 TPa, E G 22 = 1.87 TPa, G G 12 = 0.683 TPa, v G 12 = 0.177, ρ G = 4118 kg/m 3 , E m = 2.5 GPa, ρ m = 1150 kg/m 3 . v m = 0.34. The subscripts, G and m, represent the graphene and matrix, respectively. In the forced-vibration stage (time range 0 to 0.2 ms), the central deflection-time curve is described in Figure 6. It is found from Figure 6 that the forced response predicted from the proposed model follows the same trend and agrees well with that provided by Fan et al. [71].
Let us examine the forced-vibration behavior by applying the proposed method to the cases (25 F /20 C /−20 C /20 C ) S . Figure 7 exhibits the forced-vibration curves of the uniform distribution (UD) beam under the reference temperature (T = 300 K). Four different dynamic loads are taken into consideration. The values of the parameters used are L = 20 h, h = 1.75 mm. In the forced-vibration region (time range 0 to 0.8 ms), it can be observed that the step load leads to highest deflection among the four while lowest deflection occurs at the exponential load. At time range of 0.8 to 1.6 ms, the free vibration behavior of hybrid laminated beams with an initial deflection and velocity is triggered resulting of the dynamic load is released. Please note that a sudden uniform load is taken as the applied load in the parametric study. Figure 8 depicts central deflection versus time for UD and FG-CNTRC beam under the reference temperature (T = 300 K). Five different distribution patterns designed above are considered. The value of load amplitude is taken as Q = 0.4 MPa. Figure 8 shows deflection-time curves among the above designed distributing patterns. The curves of the UD pattern are used at the same temperature for reference. It can be observed that the lowest deflection is obtained by FG-Λ beam as the UD plate has the highest deflection. Thus, for the next parametric analysis, UD and FG-Λ are considered/studied. curve is described in Figure 6. It is found from Figure 6 that the forced response predicted from the proposed model follows the same trend and agrees well with that provided by Fan et al. [71]. Let us examine the forced-vibration behavior by applying the proposed method to the cases (25 F /20 C /−20 C /20 C ) S. Figure 7 exhibits the forced-vibration curves of the uniform distribution (UD) beam under the reference temperature (T = 300 K). Four different dynamic loads are taken into consideration. The values of the parameters used are L = 20 h, h = 1.75 mm. In the forced-vibration region (time range 0 to 0.8ms), it can be observed that the step load leads to highest deflection among the four while lowest deflection occurs at the exponential load. At time range of 0.8 to 1.6ms, the free vibration behavior of hybrid laminated beams with an initial deflection and velocity is triggered resulting of the dynamic load is released. Please note that a sudden uniform load is taken as the applied load in the parametric study.   Let us examine the forced-vibration behavior by applying the proposed method to the cases (25 F /20 C /−20 C /20 C ) S. Figure 7 exhibits the forced-vibration curves of the uniform distribution (UD) beam under the reference temperature (T = 300 K). Four different dynamic loads are taken into consideration. The values of the parameters used are L = 20 h, h = 1.75 mm. In the forced-vibration region (time range 0 to 0.8ms), it can be observed that the step load leads to highest deflection among the four while lowest deflection occurs at the exponential load. At time range of 0.8 to 1.6ms, the free vibration behavior of hybrid laminated beams with an initial deflection and velocity is triggered resulting of the dynamic load is released. Please note that a sudden uniform load is taken as the applied load in the parametric study.    Figure 8 shows deflection-time curves among the above designed distributing patterns. The curves of the UD pattern are used at the same temperature for reference. It can be observed that the lowest deflection is obtained by FG-Λ beam as the UD plate has the highest deflection. Thus, for the next parametric analysis, UD and FG-Λ are considered/studied. The above approach is also used to investigate the influence of changing the temperature field on UD and FG-Λ beams. Figure 9 shows the influence of the thermal stress on UD and FG-Λ beams. The applied load is 45 MPa for (25 F /20 C /−20 C /20 C ) S and (25 F /90 C /−90 C /90 C ) S with L = 5 h. Incremental temperatures are marked on these curves. The curves show that the amplitude and period of response increase with increase in temperature. This is again attributed to the fact that the development of thermal resultants reduces the overall plate stiffness. Figure 10 shows the maximum deflection with time variation of hybrid beams (25 F /20 C /−20 C /20 The above approach is also used to investigate the influence of changing the temperature field on UD and FG-Λ beams. Figure 9 shows the influence of the thermal stress on UD and FG-Λ beams. The applied load is 45 MPa for (25 F /20 C /−20 C /20 C ) S and (25 F /90 C /−90 C /90 C ) S with L = 5 h. Incremental temperatures are marked on these curves. The curves show that the amplitude and period of response increase with increase in temperature. This is again attributed to the fact that the development of thermal resultants reduces the overall plate stiffness.
The above approach is also used to investigate the influence of changing the temperature field on UD and FG-Λ beams. Figure 9 shows the influence of the thermal stress on UD and FG-Λ beams. The applied load is 45 MPa for (25 F /20 C /−20 C /20 C ) S and (25 F /90 C /−90 C /90 C ) S with L = 5 h. Incremental temperatures are marked on these curves. The curves show that the amplitude and period of response increase with increase in temperature. This is again attributed to the fact that the development of thermal resultants reduces the overall plate stiffness. Figure 10 shows the maximum deflection with time variation of hybrid beams (25 F /20 C /−20 C /20 C ) S and (25 F /90 C /−90 C /90 C )S for different foundations stiffnesses ((k1, k2 Cd)=(0, 0, 0), (10 2 ,0, 0), (10 2 , 10, 1), (10 2 , 10,2)). Beams without foundation (k1 = k2 = Cd = 0) are selected as a reference. The results confirm that the foundation stiffness decrease the transverse deflection. Moreover, the deflection of beams on the viscoelastic foundation decreases with increasing time.

Non-Linear Bending Analysis
In this section, non-linear bending response of hybrid laminated beams with NPR and PPR will be discussed in detail. For this purpose, we will determine relationship between applied pressure and deflection of the beam. For the static analysis, this relationship is independent of the change in time (t). Therefore, the transverse load is modified to be uniformly distributed, and Q (X, t) = Q (X) = q0. W is independent of time. Under these assumptions, Equations (16) and (17) can be simplified as:

Non-Linear Bending Analysis
In this section, non-linear bending response of hybrid laminated beams with NPR and PPR will be discussed in detail. For this purpose, we will determine relationship between applied pressure and deflection of the beam. For the static analysis, this relationship is independent of the change in time (t). Therefore, the transverse load is modified to be uniformly distributed, and Q (X, t) = Q (X) = q 0 . W is independent of time. Under these assumptions, Equations (16) and (17) can be simplified as: A perturbation approach is also used to derive the solutions of Equations (30) and (31). The solution equations can be expanded as a function with a small perturbation parameter ε j (j = 1,2,3, . . . ).
Following the perturbation solutions procedure, the solutions for the equations with the first order ε can be assumed as: The thermal load can be expanded as the Fourier modes as: For the free vibration analysis, the half-wavelength (m) can be assumed as m = 1,2,3 . . . , whereas, for the bending study m is adopted as 1. Substituting Equations (32) and (35) in the motion equations, we get the asymptotic solution as: and Equations (36) and (37)  10 ) is taken as the second perturbation parameter related to the dimensionless maximum deflection W m . From Equations (36a) and (37), the relationships between the load and central displacement are obtained as: The non-linear bending of a hybrid laminated beam with cross-ply (0 C /90 F /0 C /90 F /0 C ) under a uniformly distributed load is studied in the current analysis. Both the ends of the beam are taken as simply supported. The same analysis has been reported by Fan & Wang [72]. Constituent materials of the beam are similar to the cases of hybrid beam listed in Table 4. Present dimensionless deflections are compared with those of Fan & Wang [72] for the hybrid laminated beam with L/h = 10, 20, 40. For the sake of consistency, we have used q o L 3 /E 0 I, which is adopted by Fan & Wang [72]. Nevertheless, q o L 4 /E 0 I is non-dimensional load which is adopted in this paper. According to curves plotted in Figure 11, similar trend and good accordance may be observed between the results. The non-linear bending of a hybrid laminated beam with cross-ply (0 C /90 F /0 C /90 F /0 C ) under a uniformly distributed load is studied in the current analysis. Both the ends of the beam are taken as simply supported. The same analysis has been reported by Fan & Wang [72]. Constituent materials of the beam are similar to the cases of hybrid beam listed in Table 4. Present dimensionless deflections are compared with those of Fan & Wang [72] for the hybrid laminated beam with L/h = 10, 20, 40. For the sake of consistency, we have used qoL 3 /E0I, which is adopted by Fan & Wang [72]. Nevertheless, qoL 4 /E0I is non-dimensional load which is adopted in this paper. According to curves plotted in Figure  11, similar trend and good accordance may be observed between the results.  [72] 2&Present 2&Fan&Wang [72] 3&Present 3&Fan&Wang [72] (0 C /90 Figure 11. Comparisons of load-deflection relationships of a hybrid laminated beam. Figure 12 illustrates the dimensionless load-deflection curves for UD and FG beam under a uniform pressure. Five different distribution patterns designed above are considered. /h and qoL 4 /E0I represent the dimensionless defection and load, respectively. Figure 12 shows the deflectiontime curves among the above-mentioned distributing patterns. The curve of the UD pattern is used at the same temperature for reference. UD beam has the maximum central deflection as compared with the other distribution patterns while FG-Λ beam has the minimum central displacement. Thus, for the next comprehensive studies, UD and FG-Λ are taken as the case studies.  Figure 11. Comparisons of load-deflection relationships of a hybrid laminated beam. Figure 12 illustrates the dimensionless load-deflection curves for UD and FG beam under a uniform pressure. Five different distribution patterns designed above are considered. W m /h and q o L 4 /E 0 I represent the dimensionless defection and load, respectively. Figure 12 shows the deflection-time curves among the above-mentioned distributing patterns. The curve of the UD pattern is used at the same temperature for reference. UD beam has the maximum central deflection as compared with the other distribution patterns while FG-Λ beam has the minimum central displacement. Thus, for the next comprehensive studies, UD and FG-Λ are taken as the case studies. Figure 13 demonstrates the influence of the changing temperature field on the bending response of UD and FG-Λ with the above-mentioned conditions and L = 5 h. It may be concluded from Figure 13 that the central deflection increases with a rise in temperature. It is also observed that (25 F /20 C /−20 C /20 C ) S laminated beams are less sensitive to variations in the temperature than (25 F /90 C /−90 C /90 C ) S laminated beam. Figure 14 plots the effect of the foundation stiffnesses on the bending characteristics of UD and FG-Λ with reference environmental conditions and L = 20 h. The predicted load-deflection curves are given for (k 1 , k 2 ) = (0, 0), (10 3 , 0), (10 3 , 10 2 ). In Figure 14, it is evident that the dimensionless central deflection decreases by increasing the foundation stiffnesses. Meanwhile, the effect of reinforcement of foundation on (25 F /90 C /−90 C /90 C ) S laminated beam is more obvious than the other beams. qoL 4 /E0I represent the dimensionless defection and load, respectively. Figure 12 shows the deflectiontime curves among the above-mentioned distributing patterns. The curve of the UD pattern is used at the same temperature for reference. UD beam has the maximum central deflection as compared with the other distribution patterns while FG-Λ beam has the minimum central displacement. Thus, for the next comprehensive studies, UD and FG-Λ are taken as the case studies.  Materials 2020, 13, x FOR PEER REVIEW 19 of 25 Figure 13 demonstrates the influence of the changing temperature field on the bending response of UD and FG-Λ with the above-mentioned conditions and L = 5 h. It may be concluded from Figure  13 that the central deflection increases with a rise in temperature. It is also observed that (25 F /20 C /−20 C /20 C ) S laminated beams are less sensitive to variations in the temperature than (25 F /90 C /−90 C /90 C ) S laminated beam. Figure 14 plots the effect of the foundation stiffnesses on the bending characteristics of UD and FG-Λ with reference environmental conditions and L = 20 h. The predicted load-deflection curves are given for (k1, k2) = (0, 0), (10 3 , 0), (10 3 , 10 2 ). In Figure 14, it is evident that the dimensionless central deflection decreases by increasing the foundation stiffnesses. Meanwhile, the effect of reinforcement of foundation on (25 F /90 C /−90 C /90 C ) S laminated beam is more obvious than the other beams.

Conclusions
In this study, a model for free vibration, non-linear forced vibration, and bending analyses of hybrid laminated beams with NPR and PPR under various external conditions is proposed. The model developed accounts for the FG configurations and NPRs in hybrid laminated beams. Five   Figure 13 demonstrates the influence of the changing temperature field on the bending response of UD and FG-Λ with the above-mentioned conditions and L = 5 h. It may be concluded from Figure  13 that the central deflection increases with a rise in temperature. It is also observed that (25 F /20 C /−20 C /20 C ) S laminated beams are less sensitive to variations in the temperature than (25 F /90 C /−90 C /90 C ) S laminated beam. Figure 14 plots the effect of the foundation stiffnesses on the bending characteristics of UD and FG-Λ with reference environmental conditions and L = 20 h. The predicted load-deflection curves are given for (k1, k2) = (0, 0), (10 3 , 0), (10 3 , 10 2 ). In Figure 14, it is evident that the dimensionless central deflection decreases by increasing the foundation stiffnesses. Meanwhile, the effect of reinforcement of foundation on (25 F /90 C /−90 C /90 C ) S laminated beam is more obvious than the other beams.

Conclusions
In this study, a model for free vibration, non-linear forced vibration, and bending analyses of hybrid laminated beams with NPR and PPR under various external conditions is proposed. The model developed accounts for the FG configurations and NPRs in hybrid laminated beams. Five types of volume fractions of CNTs are considered which include: UD, FG-Λ, FG-V, and FG-X. For a

Conclusions
In this study, a model for free vibration, non-linear forced vibration, and bending analyses of hybrid laminated beams with NPR and PPR under various external conditions is proposed. The model developed accounts for the FG configurations and NPRs in hybrid laminated beams. Five types of volume fractions of CNTs are considered which include: UD, FG-Λ, FG-V, and FG-X. For a configuration of (25 F /20 C /−20 C /20 C ) S and (25 F /90 C /−90 C /90 C ) S , the out-of-plane Poisson's ratios exhibit NPR and PPR, respectively. The solutions for the static and dynamic responses are derived by using a two-step perturbation approach by using a fourth-order Runge-Kutta method. The theoretical model predicts reliable values of frequencies and deflections in comparison with the results available in the existing model in the literature. The influence of the temperature variations and elastic foundation on the static and dynamic behavior of hybrid laminated beams is investigated in detail.
Following comprehensive results are concluded: • The distribution type of CNTs show considerable effect on the static and dynamic behavior of FRC/CNTRC plate with NPR or PPR. It is concluded that the FG-X beams have a lower frequency ratio. Moreover, for the non-linear bending and dynamic response, the FG-Λ beam showed better performance than that with UD under different external conditions. • An increment in the temperature make the frequencies ratio and central defection considerably larger, whereas the increasing foundation stiffness will result in an opposite effect.

•
The dynamic responses considering viscosity of foundation show obvious differences from those obtained with Pasternak foundation.

•
The comparative studies reveal that the non-linear vibration of (25 F /20 C /−20 C /20 C ) S laminated beams are more sensitive to the temperature field, whereas, the non-linear bending in the (25 F /20 C /−20 C /20 C ) S laminated beams are less sensitive to the changes in temperature and foundation coefficient than those with (25 F /90 C /−90 C /90 C ) S .
The results of the current work demonstrate that the NPR has a significant effect on the non-linear bending, free vibration, and forced-vibration characteristics of FRC/CNTRC hybrid laminated beams. Furthermore, the results show that the thermal-mechanical characteristics of FRC/CNTRC hybrid laminated beam can be remarkably improved by the determination of proper distribution of CNTs. For example, to have satisfactory performance in high temperature environments from hybrid laminated beams, beams with FG-X or FG-Λ arrangement can be preferred. Meanwhile, it is hoped that the current work will give an insight into the non-linear behavior of hybrid laminated beams with NPR and will be helpful for the further investigations.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix B
In Equation (