Multiscale simulation of temperature- and pressure-dependent nonlinear dynamics of PMMA/CNT composite plates

Owing to the excellent mechanical properties, polymethyl methacrylate (PMMA)/carbon nanotube (CNT) composite has been increasingly adopted in the aerospace field, which is usually subjected to various temperature conditions and supersonic aerodynamic loads. Using a molecular dynamics (MD)-based multiscale simulation, the nonlinear forced vibration of PMMA/CNT composite plate is investigated under coupled temperature and aerodynamic load conditions. The longitudinal, transverse and shear moduli of PMMA/single-walled CNT (SWCNT) nanocomposite obtained from MD simulation at different temperature and pressure levels are substituted into the extended rule of mixtures to establish the constitutive equations of PMMA/CNT composite plate. Meanwhile, Poisson’s ratio and thermal expansion coefficient of nanocomposite are used in the constitutive equations. Based on the third-order shear deformation theory, von Karman nonlinear strain–displacement relation and Hamilton’s principle, the partial differential equation of composite plate is derived, which is reduced into a set of coupled ordinary differential equations by applying Galerkin method and is solved using the fourth-order Runge–Kutta method. The bifurcation diagrams, phase portraits, time histories and Poincare maps of composite plate are obtained under the complex loads including transverse harmonic excitation and aerodynamic pressure, which are applied to analyze the dynamic characteristics of system. This study reveals the nonlinear dynamic characteristics of PMMA/CNT composite plate, which contributes to the prediction of long-term performance of composite materials in aerospace field.


Introduction
Polymethyl methacrylate (PMMA) has been widely applied as a replacement for glass in the applications of porthole and windshield of airplane owing to the excellent transparency and easy processing [1,2]. However, as the PMMA is fragile in stress concentration zone and possesses low melting temperature, it is considered as one of the vulnerable parts of the airplane. Because of the superior mechanical properties and the high resistance to thermal stress, PMMA/carbon nanotube (CNT) composite has been increasingly adopted in the aerospace application [3,4]. During the working condition, the airplane is inevitably subjected to various temperatures and aerodynamic pressures, which cause a series of complex nonlinear dynamic behaviors [5][6][7][8][9][10]. Specifically, the increase in temperature level leads to the reduced fundamental frequencies and the increased nonlinearto-linear frequency ratios of composite plate, which affect the structural dynamic response of the airplane. Meanwhile, the sequence of modes of the composite plate changes with the increased aerodynamic pressure [11]. The long-term vibration phenomenon at different temperature and aerodynamic pressure levels usually leads to the unpredictable dynamic behavior and the consequent failure of composite structure. In order to understand the long-term performance of PMMA/CNT composite in the airplane, it is a key issue to investigate the nonlinear dynamic response of PMMA/CNT composite in the complex temperature and aerodynamic pressure conditions.
Many studies have been devoted to the synthesis and preparation of PMMA/CNT composites using in situ polymerization, melt mixing, solvent casting and other methods, and the homogeneous dispersion of CNT in PMMA matrix could be achieved [3,[12][13][14][15][16]. Meanwhile, several techniques have been developed to achieve the orientation and high concentration of nano-filler in polymer matrix, including layer-by-layer assembly technique, solution phase mixing and spin coating technique [17]. Nevertheless, there are few reports on the manufacturing of oriented CNT-reinforced PMMA composite plate. Meanwhile, previous experimental works have focused on the vibration response of composite plate. With the help of microcontroller-type thermal chamber and the impacttype vibration analyzer with virtual programming, the experimental frequency of CNT-reinforced composite (CNTRC) plate was recorded under the elevated temperature conditions [18]. It was shown that the vibration frequency of CNTRC plate decreased with the increase in temperature, as the stiffness of plate decreased [18]. Moreover, the oscillation amplitudes of cantilevered flexible plate under the effect of axial flow were measured in a low-turbulence wind tunnel with the increased flow velocity by a laser sensor [19]. It was shown that when the flow velocity was relatively low but still above the flutter point, the oscillations of plates were periodic and almost harmonic according to the measured time history of displacement. However, to the author's knowledge, the simultaneous changes of material properties and aerodynamic pressure caused by altitude conditions are difficult to be implemented in the experiment. Therefore, the use of numerical methods is warranted to obtain preliminary insights on the nonlinear dynamics of the investigated composite plates.
The numerical simulation has been considered as an effective approach in analyzing the nonlinear behaviors of composite in the complex environment and load conditions [20][21][22][23]. With the rapid development of computer technology, a variety of commercial finite element software have been used in the investigation of vibration behavior for plate structures, and the simulation results are in good agreement with experimental and theoretical analysis [24][25][26]. It is noticed that the viscous damping is usually adopted in the theoretical analysis, which is proportional to the speed of system. Comparatively, for the nonlinear forced vibration analysis in the commercial finite element software, the input of damping parameter is Rayleigh damping, which is related to the scale factors, mass and stiffness matrices of the system. The scale factors are derived from the damping ratio (generally between 0 and 1) and obtained from the experimental measurement of the corresponding material structure. It is noted that the damping ratio of investigated PMMA/CNT composite structures has not been reported in the existing literature, and hence, there is still difficulty in performing the nonlinear dynamic simulations of PMMA/CNT composite structures using commercially available finite element software. Comparatively, the theoretical modeling and numerical solution methods are usually adopted to investigate the dynamic behaviors of the PMMA/CNT composite structures. It was found that with the temperature increased from 300 to 700 K, the natural frequency of PMMA/CNT composite plate decreased and the nonlinear-to-linear frequency ratio increased, which indicated that the increased temperature enhanced the nonlinear behavior of the plate [10,27,28]. Meanwhile, it was reported that the maximum deflection of functionally graded (FG)-based PMMA/CNT composite plate occurred at the center of the flat panel without the aerodynamic pressure, while for the flat panel subjected to the aerodynamic pressure, the maximum deflection of the panel moved to the opposite direction of the airflow and the vibration mode shapes were changed accordingly [29,30]. From these previous investigations, it is found that both temperature and aerodynamic pressure conditions cause the complex nonlinear behavior of composite system. However, it is not clear about the nonlinear dynamic response of PMMA/CNT composite plate in the coupled effect of temperature and pressure conditions. Meanwhile, the mechanism of nonlinear dynamic response under the complex environment and aerodynamic pressure conditions is still lacking for composite structures.
In order to fully understand the complex dynamic response of structure, several multiscale methods have been developed to study the vibration behavior of structure. Specifically, the scale effects (also called as nonlocal effects or size effects) and the long-range atomic interactions are considered in nonlocal theory [31,32], which is mainly applicable to the elastic static and dynamic analysis for nano-and micro-structure [33][34][35]. Meanwhile, fractional calculus has been used to model a variety of nonlocal, multiscale and anomalous problems [36][37][38]. In the fractional calculus approach, fractional derivatives are the differintegral type of operator, which are intrinsically multiscale and provide a way to consider nonlocal and memory effects. In the context of elasticity, spatial fractional derivatives are often used to describe nonlocal constitutive relations, which could be an additional way to model composite plates. Furthermore, based on the bottom-up sequential multiscale conception, the molecular dynamics (MD)-based multiscale approach was developed to investigate the nonlinear vibration of PMMA/CNT composite plate [17,[39][40][41]. The nanoscale effects including the size of reinforcement, polymer chain length and long-range atomic interactions are considered in MD simulation, which are derived as the efficiency parameters into the extended rule of mixtures into the simulation at macroscale. The extended rule of mixtures model has been validated to give the effective elastic properties of macroscopic PMMA/CNT composite plates [42]. In our recent simulation study, it was reported that the distributed CNT close to top and bottom surfaces caused high natural frequency and low nonlinear frequency of PMMA/CNT composite plate, as the CNT distribution is efficient to increase the stiffness of composite plate [43]. Meanwhile, the effect of CNT volume fraction was studied, where the natural frequency of FG-based PMMA/CNT composite plate increased with the increased CNT volume fraction, and the nonlinear-to-linear frequency ratio decreased correspondingly due to the enhanced stiffness [44]. Furthermore, the nonlinear dynamic responses of FGbased PMMA/CNT composite beam were investigated under the axial airflow in thermal conditions [45]. It was demonstrated that the thermal stress induced by the temperature difference significantly reduced the flutter boundary. In addition, it was shown that with the elevated aerodynamic pressure, the first-order frequency of FG-CNTRC panel increased, while the second-order frequency decreased gradually, and the frequency coalescence between the first and second modes leads to the flutter instability [46]. Among these studies, the temperature-independent efficiency parameters are introduced into the extended rule of mixtures (EROM), and the prediction of CNT-reinforced composite plate may not be strictly accurate as the effect of temperature is not considered in the calculated mechanical parameters. Meanwhile, there are few investigations relevant to the nonlinear dynamic responses in the coupled temperature and aerodynamic conditions.
In this paper, the nonlinear dynamic response of PMMA/CNT composite plate is investigated in the complex temperature and aerodynamic pressure conditions. The longitudinal, transverse, shear moduli, Poisson's ratio and thermal expansion coefficient (TEC) of PMMA/SWCNT nanocomposite are determined using MD simulation at different temperature and pressure levels. Subsequently, the obtained elastic properties from MD simulation are substituted into the EROM to calculate the effective longitudinal, transverse and shear moduli of PMMA/CNT composite plate. Both TEC and elastic properties are applied to establish the constitutive equations of composite plate in the partial differential equations, which are then reduced into a set of coupled nonlinear ordinary differential equations using Galerkin method. The ordinary differential equations are solved using the fourth-order Runge-Kutta method to obtain the bifurcation diagrams, time histories, phase portraits and Poincare maps of PMMA/CNT composite plate. This study contributes to the nonlinear dynamic behavior of PMMA/CNT composite plate in the complex temperature and aerodynamic pressure conditions.

Molecular simulations
The MD simulation adopted in this study includes the model construction of SWCNT, PMMA and PMMA/ SWCNT nanocomposite, and the structural relaxation and dynamic deformation. After that, the longitudinal, transverse, shear moduli, Poisson's ratio and TEC of molecular models are obtained, which are applied as the inputs in the constitutive equation for the subsequent macromechanical simulations. The details of modeling and simulation are provided in the following sections.

Force field and molecular model
In the molecular simulation, the Tersoff potential is applied to characterize the interactions of SWCNT, which consists of two-body interactions, three-body interactions and multibody term [47]. Meanwhile, the polymer consistent force field (PCFF) is adopted to characterize the interactions of PMMA polymer and between SWCNT and PMMA matrix in nanocomposite [48]. The potential functions in the PCFF include the bond stretching, angle bending, dihedral angle torsion, improper out-of-plane and cross-coupling terms of the bonded interactions, and the van der Waals (vdW) and Coulombic terms of the non-bonded interactions with a cutoff distance of 1.25 nm [49,50]. The potential functions of Tersoff and PCFF potentials are described in our recent study [43]. In the construction of the molecular models, a periodic (5, 5) SWCNT model with a diameter of 0.678 nm and a length of 6.0 nm is investigated, as shown in Fig. 1a, which has been extensively adopted as the reinforcing filler in polymer/SWCNT nanocomposites [49,51]. For the PMMA polymer, isotactic PMMA (i-PMMA) is investigated here, which is widely used in various engineering applications due to the excellent mechanical properties [52,53]. The molecular model of i-PMMA polymer chain with 50 repeating units is shown in Fig. 1b. The PMMA polymer developed in our previous study is used here, where 31 i-PMMA chains are packed into a periodic box with an initial density of 1.08 g cm -3 , which is close to the experimental range of 1.15-1.19 g cm -3 [41,[52][53][54][55][56]. For the nanocomposite, the configuration is constructed that only a small number of SWCNT segments are embedded into the matrix [49]. Hence, the single (5, 5) SWCNT segment is placed in the center of a periodic simulation box, and 31 i-PMMA chains are then packed randomly around the tube in non-overlapping positions. The molecular model of nanocomposite is considered as the representative model at the molecular level, as shown in Fig. 1c.

Equilibration and deformation
After the model construction, the equilibration and deformation of the molecular models are performed using the open-source code LAMMPS with a time step of 1 fs [57]. During the energy minimization, the SWCNT model is equilibrated in the microcanonical (NVE) ensemble for 50 ps followed by another equilibration in the isothermal and isochoric (NVT) ensemble at a constant temperature of 100 K for 50 ps. After that, the system is further equilibrated in the isobaric ensemble (NPT) ensemble. During the NPT simulation, the temperature ramps up from the 100 to 300 K with a ramp rate of 100 K per 100 ps, and the pressure is set as 1 atm in each case. For PMMA polymer and PMMA/SWCNT nanocomposite, the constructed models are equilibrated in the NVE ensemble for 50 ps, then in the NVT ensemble at 300 K for another 50 ps and finally in the NPT Fig. 1 Schematic diagram of MD-based multiscale approach for CNT-reinforced PMMA composite: molecular models of a SWCNT, b PMMA chain and c nanocomposite; d equivalent nanocomposite; and e macroscale PMMA/CNT composite plate ensemble at 300 K and 1 atm for 100 ps. After that, the PMMA polymer and PMMA/SWCNT nanocomposite models are further equilibrated using a validated equilibration procedure with seven relaxation cycles to achieve the fully relaxed state [58][59][60][61]. In every cycle, the models are firstly equilibrated in an NVT ensemble at 800 K temperature for 50 ps, then at 300 K for another 50 ps and finally in NPT ensemble at 300 K with a changing pressure. Specifically, the imposed pressure in the first three cycles increases steadily from 1 to 50,000 atm, during which the structure is compressed efficiently. In the subsequent four equilibration cycles, the imposed pressure decreases gradually to 1 atm, during which the model is decompressed efficiently. When the pressure reaches 1 atm, the model is subjected to a final relaxation process of 800 ps at the temperature of 300 K and pressure of 1 atm.
To investigate the coupled effect of temperature and pressure conditions, the equilibrated systems are conditioned at different temperature and pressure levels. Notably, the flight altitude of aircraft is generally less than 13 km above sea level [62]. Accordingly, the four flight conditions are investigated in this study, including 0.0, 4.0, 8.0 and 12.0 km. The temperature and pressure parameters corresponding to these four altitude conditions are derived from the 1976 U.S. Standard Atmosphere, as listed in Table 1 [63]. In order to reach these conditions, a further equilibration process is performed on the SWCNT, PMMA polymer and nanocomposite models, with the temperature decreasing from 288.15 to 216.65 K and the pressure decreasing from 1 to 0.1915 atm, respectively. In each investigated condition, the model is equilibrated in an NVT simulation for 0.5 ns followed by an NPT equilibration for 1.0 ns, and the detailed process is shown in Table 1. It is noted that the cooling rate is close to those used in previous simulation studies [43,64]. The root-mean-squared displacement, center of mass, temperature and potential energy of the molecular PMMA polymer and PMMA/SWCNT nanocomposite models are recorded during the 1.0-ns NPT equilibration in each investigated case, as shown in Fig. 2, and the constant oscillation demonstrates that the models have reached the equilibrium state.
After the equilibration procedure, the extended equilibration is performed to measure the TEC of PMMA/SWCNT nanocomposite in each investigated condition as shown in Table 1. Specifically, the model is subjected to a stepwise cooling process starting from 288.15, 262.15, 236.15 and 216.65 K with a ramp rate of 5 K per 100 ps until the temperature drops 25 K, respectively. In each investigated condition, the model is relaxed in the NVT ensemble for 100 ps followed by an NPT relaxation for 100 ps, and the pressure is constant in the entire simulation. The volume and temperature are recorded in the cooling procedure. Additionally, the SWCNT is subjected to a tensile deformation in an NVT ensemble to quantify the Poisson's ratio, in which the strain rate is set as 1.6 9 10 10 s -1 , and Poisson's ratio is determined by the ratio of longitudinal strain to transverse strain. Based on the obtained models, the longitudinal, transverse and shear moduli of SWCNT, PMMA polymer and PMMA/SWCNT nanocomposite at different temperature and pressure levels are quantified using a constant-strain method, where the strain rates are set as 2 9 10 10 s -1 in the longitudinal and transverse deformations, and the shear deformation along the xy plane [65,66]. The deformation process repeats until the strain reaches 3%. During the entire MD simulation, the constant temperature and pressure are controlled using the Nose-Hoover thermostat and the Andersen barostat, respectively [67,68].

Macroscale simulation
By using the extended rule of mixtures, the molecular model of PMMA/CNT nanocomposite is homogenized as a representative volume element, as shown in Fig. 1d. Meanwhile, the material parameters of SWCNT, PMMA and nanocomposite obtained from MD simulations are used as inputs to investigate the nonlinear dynamic response of FG-based composite plates. In the third-order shear deformation theory (TSDT), the higher-order terms of the in-plane displacement fields are considered without shear correction coefficient, and the condition of zero shear stress at the top and bottom surfaces is satisfied, which appropriately increases the computational complexity while widely applicable for plates with different length-to-thickness ratios [69].
Comparatively, the first-order shear deformation theory considers the effect of the transverse shear deformation by introducing a shear correction factor and provides the reliable results for the thin and moderately thick plates [70,71]. In order to avoid the selection and verification of shear correction coefficient, TSDT is adopted in this study. According to the TSDT, von Karman nonlinear strain-displacement relation and Hamilton's principle, the partial differential equation of composite plate is derived. The ordinary differential governing equation of motion is obtained by applying Galerkin method and solved by the fourth-order Runge-Kutta method. The phase portraits and time histories of FG-based composite plates at different altitudes are obtained under the transverse harmonic excitation and aerodynamic pressure to analyze the dynamic characteristics of system.

Effective material properties
The mechanical parameters obtained from MD simulation are transformed to the inputs of nonlinear dynamic equation of PMMA/CNT composite plate at different altitudes, as shown in Fig. 1e. In this study, the FG-based composite plate is constructed with the E composite where g x , g y and g xy are the corresponding efficiency parameters with respect to the longitudinal, transverse and shear deformation; E x , E y and G xy are the longitudinal, transverse and shear modulus of SWCNT, PMMA polymer and PMMA/SWCNT nanocomposite, respectively; and V CNT and V PMMA are the volume fraction of SWCNT and PMMA matrix, respectively. In this study, V CNT is set as 1.0 vol%. The CNT distribution of the composite plate is considered with uniformly distribution (UD) and FG-X distribution, where X denotes the shape of CNT distribution, as shown in Fig. 3b-c. The mathematical expressions of the two distributions are given as follows [73]: where V CNT Ã is the CNT volume fraction of composite plate, h denotes the thickness of the plate and z is the coordinate along the thickness direction. By substituting V CNT into Eq. (1), the longitudinal, transverse and shear moduli of composite plate are obtained based on the derived efficiency parameters. Poisson's ratio of PMMA/CNT composite is derived based on the following equations: where m composite xy , m CNT and m PMMA are Poisson's ratio of PMMA/SWCNT nanocomposite, SWCNT and PMMA polymer, respectively. Meanwhile, the linear TEC of composite is derived using the following equation [17]: where a is the linear TEC of composite, L 0 is the initial length, L i is the final length and DT is the temperature rise. For the plate subjected to the high supersonic flow, the quasi-steady first-order piston theory [74][75][76] is commonly used as where P a , M a and t are the aerodynamic pressure, Mach number and time, respectively, w is the transverse displacement of the plate in the middle plane and V a and q a denote the velocity of incoming flow and air mass density, respectively, which are obtained via the ideal gas law [77] and Bernoulli equation [78]  where P i and T i are the pressure and temperature in Case #i (i = 1,2,3,4), respectively, and M, R * , g, h and P 0 are the molar mass of air, universal gas constant, earth-surface gravitational acceleration, altitude above sea level and sea level standard atmospheric pressure of 1 atm, respectively.

Governing equations
Reddy's third-order shear deformation theory [79] is adopted to establish the governing equations as follows: where c 1 = 4/(3h 2 ), u 0 , v 0 and w 0 are displacement components of middle plane corresponding to the coordinates (x, y, z) and / x and / y are the rotations of the transverse normal in middle plane with respect to the y-and x-axes, respectively. The von Karman strain-displacement relation is obtained as where The constitutive relations considering the thermal effects are written as where In addition, we assume that G composite [27]. By applying Hamilton's principle where dT, dU and dW are the kinetic energy, potential energy and external work, respectively, raised from the virtual displacement with The strain can be written in terms of the generalized displacements using Eq. (10). Substituting Eq. (15) dT ¼ À N xy;y du 0 À N xy;x dv 0 À M y;y d/ y À N y;y dv 0 þ c 1 P y;y d/ y À c 1 P y;yy dw 0 À N y;y w 0;y dw 0 À N y w 0;yy dw 0 into Eq. (14), the equations of motion are obtained as follows: where The stress resultants are related to the strains by the relations where and the thermal stress resultants are represented as All terms of Eqs. (16) are shown in Appendix. By combining ''Appendix'' and Eqs. (16)(17)(18)(19)(20), the nonlinear equations of motion in terms of generalized displacements (u 0 , v 0 , u x , u y and w 0 ) are It is worth mentioning that structure under different boundary conditions presents various nonlinear dynamical behaviors or nonlinear properties. For instance, it was found that for the double-layered nanoplates with clamped conditions, the minimum points of nonlinearity corresponding to the aspect ratios for each nonlocal parameter are larger than those of the simply supported conditions [80]. Comparatively, the simply supported double-layered nanoplates with immovable edges more likely appear chaotic vibration compared with movable edges [81].
In this work, the common boundary conditions for the simply supported plate with movable edges are considered as [82] v For the boundary condition with immovable edge, it affects the cubic nonlinear term, which may increase the structural stability due to the increased displacement constraints in plate edge. In future research, we will consider the effects of different boundary conditions on structural nonlinear dynamic behavior in detail.
The modal expansions for the generalized displacements u 0 , v 0 , w 0 , / x and / y are described as  Fig. 6 Bifurcation diagrams, phase planes, time histories and Poincare maps of PMMA/CNT composite plate with UD distribution pattern under nonlinear forced vibration in Case #1: the (a, b) bifurcation diagrams of F versus W 1 and W 2 , (c, bifurcation diagrams, phase planes, time histories and Poincare maps of PMMA/CNT composite plate with UD distribution pattern under nonlinear forced vibration in Case #1: the (a, b) bifurcation diagrams of F versus W 1 and W 2 , (c, e) phase portraits on plane (W 1 , _ W 1 ) and plane (W 2 , _ W 2 ), (d, f) time history diagrams on plane (t, W1) and plane (t, W 2 ), (g, h) Poincare maps on plane (W 1 , _ W 1 ) and plane (W 2 , _ W 2 )) phase portraits on plane (W 1 , _ W 1 ) and plane (W 2 , _ W 2 ), (d, f) time history diagrams on plane (t, W1) and plane (t, W 2 ), (g, h) Poincare maps on plane (W 1 , _ W 1 ) and plane (W 2 , _ W 2 ) The transverse excitation is described as where F ij denotes the amplitude of the transverse excitation corresponding to the modes of transverse deflection. By substituting Eqs. (22)(23)(24) into Eqs. (21a)-21e) and applying the Galerkin procedure to discretize the equations, the four coupled algebraic equations are obtained, and the expressions of u, v, / x and / y only contain w. The nonlinear ordinary differential equation containing only W is obtained as where W = [W 1 , W 2 ] T is the transverse displacements vectors, l is the damping terms including the aerodynamic damping, M is the mass matrix, NL is the nonlinear stiffness vector and F is the amplitude excitation matrix. Because the transverse vibration of plate is mainly considered in this work, in the ordinary differential equations derived by Galerkin method, the in-plane inertia terms and the transverse inertia term in the expressions of u, v, / x and / y derived from the kinetic energy can be ignored to reduce the computational complexity and the accuracy of the results is guaranteed, as reported in Ref. [83]. Based on Eq. (25), the time histories and phase portraits of PMMA/CNT composite plate are obtained to analyze the dynamic characteristics.

Results and discussion
The molecular simulation results of SWCNT, PMMA polymer and PMMA/SWCNT nanocomposite in various temperature and pressure conditions are presented herein, which are used in the multiscale simulation of b Fig. 7   FG-based composite plates, and the obtained nonlinear dynamic responses are discussed.

Nano-mechanical analysis
During the last 1.0-ns NPT relaxation process, the densities of PMMA polymer and PMMA/SWCNT nanocomposites are recorded, and the data from the last 50 ps are averaged as shown in Table 2. The density of PMMA polymer in Case #1 is calculated as 1.184 ± 0.001 g cm -3 , which is close to the experiment result of 1.190 g cm -3 and the simulation result of 1.180 g cm -3 at room temperature [56,84]. With the decrease in temperature and pressure, the density shows a slight increase, with the value of 1.190 ± 0.002, 1.193 ± 0.001 and 1.195 ± 0.002 g cm -3 in Case #2, #3 and #4, respectively. The observed increasing trend of density agrees very well with the experimental measurement of other polymer, where the density continues to increase as the temperature decreases from 343.15 to 283.15 K, and the pressure decreases from 35 to 0.1 MPa [85]. With the addition of 1 vol% SWCNT, the density of PMMA/SWCNT nanocomposites shows a consistent increase compared with pure PMMA polymer in all the investigated cases, with the value of 1.202, 1.206, 1.210 and 1.214 g cm -3 in Case #1, #2, #3 and #4, respectively, as shown in Table 2. The increment of about 1.5% due to the added 1 vol% SWCNT has been observed in previous experimental and simulation studies of other polymer/CNT nanocomposites at the temperature of 300 K and pressure of 1 atm, with the increases of 1.4% in experiment and of 1.9% in MD simulation [86,87]. It is shown that a closely packed PMMA structure can be formed in the vicinity of the PMMA-SWCNT interface, which could contribute to the increase in the PMMA/SWCNT nanocomposite density [41,86]. After characterizing the density, the TEC of PMMA/SWCNT nanocomposites in the various temperature and pressure conditions is evaluated. Accordingly, the longitudinal and transverse lengths of simulation cell are recorded, as shown in Fig. 4(a-d) corresponding to Cases #1-#4. It is observed that the length of PMMA/SWCNT nanocomposite model increases linearly with the increased temperature. The longitudinal and transverse TEC of PMMA/ SWCNT nanocomposite is calculated according to Eq. (5), as shown in Table 2. In Case #1, the longitudinal and transverse TEC is calculated as 4.39 9 10 -5 and 4.41 9 10 -5 /K, which are close to the calculated results of the Schapery model based on the molecular simulation results of 4.46 9 10 -5 and 4.47 9 10 -5 /K, respectively [88]. With the decreased temperature and pressure in Case #2 to #4, the TEC decreases continuously, with the values of 4.28 9 10 -5 , 4.09 9 10 -5 and 3.69 9 10 -5 /K for longitudinal direction and 4.31 9 10 -5 , 4.17 9 10 -5 and 3.89 9 10 -5 /K for transverse direction, respectively. It is probably because that the decreased temperature leads to the shrinkage of simulation box and the decreased deformation of molecular cell [89]. The decrease trend of TEC with the decrease in temperature for PMMA/SWCNT nanocomposite agrees very well with other polymer/CNT nanocomposite [90][91][92]. It is noted that the TEC in the longitudinal direction is smaller than that in the transverse direction, as the SWCNT restricts the expansion of matrix in the longitudinal direction and enforces the matrix expansion in the transverse direction [93].
After determining the density and TEC, the Poisson's ratio of SWCNT, PMMA polymer and PMMA/ SWCNT nanocomposite in the various temperature b Fig. 8 Phase planes and time histories of PMMA/CNT composite plate under nonlinear forced vibration in Case #2: the a, c phase portraits with UD distribution pattern on planes (W 1 , _ W 1 ) and (W 2 , _ W 2 ), b, d time history diagrams with UD distribution pattern on planes (t, W 1 ) and (t, W 2 ), e, g phase portraits with FG-X distribution pattern on planes (W 1 , _ W 1 ) and (W 2 , _ W 2 ), f, h time history diagrams with FG-X distribution pattern on planes (t, W 1 ) and (t, W 2 ) and pressure conditions is calculated, as listed in Table 2. In Case #1, Poisson's ratios of SWCNT, PMMA polymer and PMMA/SWCNT nanocomposite are calculated as 0.2400, 0.3100 and 0.3093, respectively, which are close to the corresponding experimental and simulation results of 0.25, 0.33 and 0.33 [94][95][96]. It is seen that with the decrease in temperature and pressure in the investigated cases, Poisson's ratio of SWCNT, PMMA polymer and PMMA/ SWCNT nanocomposite shows a monotonous decrease, which agrees with finite element simulation, MD simulation and rule of mixtures [88,95,97], where the values show a similar decreased trend with the decreased temperature. It is probably because that the transverse moduli of SWCNT, PMMA polymer and PMMA/SWCNT nanocomposite increase with the decrease in temperature and pressure. As the deformation process is determined by the thermal energy and strain energy, the thermal energy is reduced with the decreased temperature, and thus, the strain energy required for the deformation increased [89], which leads to the decreased transverse deformation during the longitudinal tensile deformation, causing the decrease in Poisson's ratio [98].
Apart from the density, TEC and Poisson's ratio, the moduli of SWCNT, PMMA polymer and PMMA/ SWCNT nanocomposite are quantified using the MD deformation. The stress-strain curves of SWCNT, PMMA polymer and PMMA/SWCNT nanocomposite are shown in Fig. 5, where L, T and S denote the longitudinal, transverse and shear deformation, respectively. From the stress-strain curves depicted in Fig. 5(a), it is shown that the longitudinal stress of SWCNT is higher than those along the transverse and shear directions. In Case #1, the longitudinal, transverse and shear moduli of SWCNT are estimated as 2099.48 ± 0.02, 151.02 ± 0.01 and 720.25 ± 0.06 GPa, respectively, as shown in Table 2. The longitudinal and transverse moduli are close to the simulation results of 2000 and 142 GPa, respectively, and the shear modulus is in the range of 550-860 GPa [99][100][101]. With the decrease in temperature and pressure in other cases, there is no obvious variation of the longitudinal, transverse and shear stresses.
For PMMA polymer, it is observed from Fig. 5b that the longitudinal and shear stress level increases gradually with the decreased temperature and pressure in the investigated cases, and the longitudinal stress is higher than the transverse result. The longitudinal, transverse and shear moduli of PMMA polymer are summarized in Table 2. In Case #1, it is evaluated that the polymer has a longitudinal modulus of 2.72 ± 0.01 GPa, which is in the range between 2.24 and 3.80 GPa of experiments [41]. The shear modulus of the polymer is estimated as 1.04 ± 0.04 GPa, which is close to MD simulation result of 1.19 GPa [95]. With the decrease in temperature and pressure, the longitudinal modulus increases continuously, with the value of 3.82 ± 0.01, 4.88 ± 0.05 and 4.96 ± 0.03 Gpa in Case #2, #3 and #4, respectively. The increasing trend is consistent with the MD results of 3.85, 4.92 and 5.01 GPa obtained at the temperature level of 262.15, 236.15 and 216.65 K, respectively [102].
With the addition of 1.0 vol% SWCNT, the nanocomposite shows the increased longitudinal stress level compared with PMMA polymer case, while no obvious increase is observed for the transverse and shear deformations, as shown in Fig. 5c. In Case #1, the longitudinal, transverse and shear modulus of nanocomposite is calculated as 21.80 ± 0.05, 3.20 ± 0.01 and 1.60 ± 0.02 GPa, which are close to the simulation results of 22.30, 3.45 and 1.62 GPa, respectively [72,103]. Compared with the PMMA polymer, the increase in longitudinal modulus of nanocomposite is because that the SWCNT carries b Fig. 9 Phase planes and time histories of PMMA/CNT composite plate under nonlinear forced vibration in Case #3: the a, c phase portraits with UD distribution pattern on planes (W 1 , _ W 1 ) and (W 2 , _ W 2 ), b, d time history diagrams with UD distribution pattern on planes (t, W 1 ) and (t, W 2 ), e, g phase portraits with FG-X distribution pattern on planes (W 1 , _ W 1 ) and (W 2 , _ W 2 ), f, h time history diagrams with FG-X distribution pattern on planes (t, W 1 ) and (t, W 2 ) most of the load in the longitudinal deformation, and the high longitudinal modulus of SWCNT results in a significant increase. Comparatively, the effects of SWCNT on the transverse and shear moduli are not obvious, as that SWCNT has a comparatively low transverse and shear modulus than the longitudinal one. With the decrease in temperature and pressure, the longitudinal, transverse and shear moduli of PMMA/SWCNT nanocomposite show similar increase as the PMMA polymer case, with the value of 22.53 ± 0.01, 23.40 ± 0.08 and 24.84 ± 0.02 Gpa in longitudinal direction. The increasing trend is consistent with the results obtained from the extended rule of mixtures for PMMA/SWCNT nanocomposite at the decreased temperature level [88].

Nonlinear dynamic analysis
The efficiency parameters of PMMA/SWCNT composite are derived through matching the longitudinal, transverse and shear moduli calculated from the EROM with those from MD simulations. The results with various CNT volume fractions are shown in Table 3. Based on the obtained efficiency parameters, the longitudinal, transverse and shear moduli of UD and FG-X plates are calculated using Eqs. (1)(2)(3). The input parameters q a and V a of the quasi-steady firstorder piston theory are calculated via Eqs. (7)(8) and listed in Table 3. The Mach number used in this study is determined as 1.5 [104].
Lyapunov exponent is usually adopted to assess the chaos characteristics of system, and the mutual mapping of bifurcation diagram and Poincare map could also verify the nonlinear dynamic regimes. Accordingly, the bifurcation diagrams and Poincare maps in Case #1 and #4 are presented to reveal the change in the nonlinear dynamic behaviors for the minimum and maximum altitude. Specifically, the bifurcation diagrams of transverse excitation F versus displacements of composite plate with the UD and FG-X distributions in Case #1 are shown in Figs. 6a, b and 7a, b, and the bifurcation diagrams in Case #4 are shown in Figs. 10a, b and 11a, b. It is clearly seen that the motion state of the plate for both Cases #1 and #4 starts from periodic motion, followed by the complex nonlinear dynamic behavior of alternating periodic motion and chaotic motion. By comparing Figs. 6 a, b and 10 a, b with the UD distributions, it is observed that in Case #4 more periodic motion and less chaotic motion are occurred than that in Case #1 with the same distribution pattern. By comparing Figs. 10 a, b and 11 a, b in Case #4, it is seen that the FG-X plate shows periodic motion in the transverse excitation range of 0-10 9 10 5 N m -2 , while the UD plate emerges chaotic motion as the transverse excitation is greater than 7.6 9 10 5 N m -2 . The simulation results of bifurcation diagrams and Poincare maps demonstrate that from Case #1 to #4, the composite structure possesses the improved motion stability due to the increased elastic modulus, and the effect of elastic modulus on the motion stability of plate structure is consistent with Ref. [105]. According to the trend of structural stability and amplitude response observed from Case #1 to #4, the effect of altitude condition on nonlinear dynamic behaviors in Case #2 and #3 could be inferred. Meanwhile, Figs Fig. 10  Figs. 6(g, h), 7(g, h), 10(g, h) and 11(g, h). For different CNT distribution types, the chaotic motion in Case #1 is demonstrated in Figs. 6(g, h) and 7(g, h), and the monoperiodic motions in Case #4 are demonstrated in Figs. 10(g, h) and 11(g, h). For the results of the two CNT distributions in each investigated case, W 1 and W 2 represent the amplitudes of the first-and second-order vibrations, respectively, and _ W 1 and _ W 2 represent the velocities of the first-and second-order vibrations, respectively. It is clearly shown that the amplitude of first-order vibration is significantly larger than the second-order vibration, and such trend is consistent with Ref. [106]. Meanwhile, the maximum amplitude of nonlinear forced vibration for FG-X type is smaller than that of UD distribution, which is due to the fact that FG-X types of CNT distribution effectively enhances the stiffness of plates, as demonstrated in Table 3, where the longitudinal, transverse and shear moduli of the composite plate with the FG-X distribution are improved significantly compared with the UD distribution. From four investigated cases, it is indicated that the amplitudes of PMMA/CNT composite plate with UD and FG-X types both decrease with the increased altitude conditions. Meanwhile, it is observed from Figs. 6-11 that the amplitude of the composite plate decreases gradually from Case#1 to #4 under the same external excitation and distribution mode. Moreover, it can be seen from Table 3 that the elastic moduli of composites increase with the increase in altitude, resulting in the increase in stiffness and the decrease in the plate amplitude under the same external excitation condition. Additionally, it could be known from Eqs. 6-8 that the aerodynamic pressure increases with the increase in altitude, which leads to the increase in the plate amplitude. To sum up, it is illustrated that with the increase in altitude, the influence of aerodynamic pressure caused by air flow parallel to the plate surface is less than that of material stiffness variation caused by the change in temperature and pressure on the transverse displacement of plate.

Conclusions
In this study, the nonlinear dynamic response of PMMA/CNT composite plate in the different temperature and pressure conditions is investigated using a MD-based multiscale approach. Firstly, MD simulation is applied to calculate the longitudinal, transverse, shear modulus, the TEC and Poisson's ratio of SWCNT, PMMA polymer and PMMA/SWCNT nanocomposite in the coupled temperature and pressure conditions. Secondly, PMMA/SWCNT nanocomposites are utilized as micro-inclusions in a micro-mechanical approach to estimate the macroscopic response of composite plate. Finally, the nonlinear dynamic response of composite plate is analyzed using macroscale simulation. The nonlinear dynamical behaviors of the PMMA/CNT composite plate become more stable with the increase in structural stiffness. Moreover, it is illustrated that with the increase in altitude, the influence of aerodynamic pressure caused by air flow parallel to the plate surface is less than that of material stiffness variation caused by the change in temperature and pressure on the transverse displacement of plate. Meanwhile, the second-order maximum vibration amplitude is significantly lower than the first-order. The nonlinear dynamic response of PMMA/CNT composite plate is studied, which contributes to the understanding of the long-term dynamic behavior of composite materials during the service life in the aerospace field. Furthermore, the quantitative transfer of physical and mechanical properties between molecular and c Fig. 11 Bifurcation diagrams, phase planes, time histories and Poincare maps of PMMA/CNT composite plate with FG-X distribution pattern under nonlinear forced vibration in Case #4: the a, b bifurcation diagrams of F versus W 1 and W 2 , c, e phase portraits on plane (W 1 , _ W 1 ) and plane (W 2 , _ W 2 ), d, f time history diagrams on plane (t, W1) and plane (t, W 2 ), (g, h) Poincare maps on plane (W 1 , _ W 1 ) and plane (W 2 , _ W 2 ) macroscale scales is established to provide a theoretical basis for the design and optimization of FG-CNTRC.
Funding This work was supprored by National Natural Science Foundation of China (Grant numbers 12072003 and 51808020).
Data availability The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Declarations
Conflicts of interest The authors declare that they have no conflict of interest.