Molecular modelling of the effect of loading rate on elastic properties of CNT-polyethylene nanocomposite and its interface

Polymer composites reinforced by carbon nanotubes (CNTs) are so attractive for many mechanical applications due to their higher strength and elastic modulus than those of polymer itself. Besides the mechanical performance of CNTs being superior, another main reason is related to the CNT-polymer interface. The interface at nanoscale level plays a key role of load transfer. These kind of polymeric materials are highly sensitive to temperature and loading rates (or strain rates), and further as the matrix with nanoparticle additives, the mechanical properties of the composite have also been found to be strongly dependent on temperature and loading rates. In this paper, specifically for the CNT-reinforced polyethylene nanocomposite, we studied the effect of loading rates on the elastic properties of the nanocomposite and its CNT-polyethylene interface by molecular dynamics (MD) modelling. Our simulated results are well agreed with the ones that are calculated on the rule-of-mixtures (ROM) and from other publications. The results show that higher loading rates lead to larger stress and Young’s modulus of the nanocomposite. And on the atomistic model with taking account of the CNT-polyethylene interface, we found that the interface thickness is independent on loading rates, while its interaction energy and Young’s modulus are related to selected CNTs, and dependent on loading rates. Both of them reach the largest value when the CNTs with the highest aspect ratios (correspondingly the highest specific surface areas) are selected and the nanocomposite being under the fastest loading rate. These results are very valuable for deeply understanding the elastic properties of CNT-polyethylene nanocomposites and their microstructural design.


Introduction
Polymer nanocomposites filled with carbon nanotubes (CNTs) have an exciting prospect among nanotechnology applications due to their remarkable enhancements in mechanical, electrical, and thermal properties in comparison with neat polymers [1][2][3]. Even carbon nanotubes with small fractions can significantly affect the mechanical properties of composites because CNTs have extremely high elastic modulus. For example, additions of 12% and 17% volume fractions of CNTs (Young's modulus: 1-5 TPa) in a polyethylene matrix (Young's modulus: several GPas) can significantly improve the Young's modulus of CNTpolyethylene nanocomposites up to 94.6 GPa and 138.9 GPa, respectively [4].
To deeply understand these attractive characteristics of CNT-polyethylene nanocomposites, we can start from their component and composite structure. In fact, for the CNT-polyethylene nanocomposite, its elastic properties, obviously, depend on the properties of the matrix of polyethylene and additive CNTs. And further, previous studies found that, to a great extent, the elastic properties are closely related to loading rates, i.e., strain rates [2]. Besides, the CNT-polyethylene interface, as a critical role in the nanocomposite of controlling stress transfer, is closely related to the adhesion and interfacial strength [1]. Of these interfacial problems, the range size and elastic property of the interface have to be firstly significant. To intensively study above the problems of CNT-polyethylene nanocomposites, besides experimental works on synthesis and preparation, microstructure analysis, characterization and properties [2,[5][6][7], the investigations through atomic and multiscale modelling of predicting the properties of CNTs and other nano-filled nanocomposites have been absolutely necessary [8][9][10][11]. In all above approaches, molecular dynamics has been as the most effective method because it bases Newton's second law of motion equations and accurately shows the nature of particle interactions [10]. Thus, for various nanocomposites of polymer matrix embedded with nano additive (1D particles, 2D sheets or/and 3D CNTs), through molecular dynamics, many studies on the mechanical behaviors [10,[12][13][14][15][16], thermal properties [17], temperature and strain rate dependences [5,[18][19][20][21], and interfacial issues [1,3,11,[22][23][24][25] have been extensively carried out nearly two decades. However, for the polymeric materials with CNTs additives, improvements of the nanocomposite in properties are by no means guaranteed, and the results are usually sensitive to the particular polymer chosen, in addition to the CNTs used in the composite.
In this paper, for the CNT-reinforced polyethylene nanocomposite, on the atomistic model of the unit cell of carbon nanotube (CNT)-reinforced polyethylene (PE) nanocomposites and by molecular dynamics simulation, we investigated the elastic properties of the nanocomposite and the CNT-polyethylene interface, emphatically the effect of loading rates. 2. Atomistic model and simulation procedure 2.1. Atomistic model A CNT-zigzag type with a chiral vector (9, 0)-having diameters of 7.05 Å, was selected as reinforcing fillers, and amorphous PE was used as the polymer matrix. In order to simplify and generically represent the polymer matrix, we selected two single chains of PE with 250 repeat units-CH 2 -in each chain as building blocks to model the actual polymer. In addition, the periodic atomistic model of CNT-reinforced polyethylene was established to avoid the termination effect of CNTs. All the atomistic models, such as serrated carbon nanotubes reinforced PE nanocomposites and pure PE matrix, were similar and modeled under the same external conditions by commercial software: Materials studio 8.0. The polymer uniform force field (PCFF) was used throughout the simulation, which was considered in the force field (section 2.2).
The atomistic model was constructed as shown in figure 1. Firstly a pristine zigzag CNT is placed in the center of a periodic unit cell (30.48 Å×30.48 Å×30.48 Å) with its length aligned alone the x-direction, as shown in figure 1(a), and, subsequently, PE chains were packed into the periodic unit cell with a density of 0.8 g cm −3 , as shown in figure 1(b). Finally as shown in figure 1(c), at the atomic level, the model of a representative unit cell of carbon nanotube (CNT)-reinforced polyethylene (PE) nanocomposites was constructed. It is noted that the steepest descent method is here used to minimize the total potential energy of the nanocomposites with energy convergence tolerance of 0.001 kcal mol −1 .

Force field
Molecular simulations are based on the accurate calculation of interactions between atoms, including bonding interactions between atoms that make up the same molecule, van der Waals interactions between different molecules, and hydrogen bond interactions among some molecules. Generally, there are two ways to describe these interactions between atoms: one by quantum chemical calculations and the other by molecular force fields. A molecular force field is sometimes called a potential function. Usually, a molecular force field potential function includes some necessary parts as follows. Bond stretching energy (E bond ): the change in energy caused by the stretching motion of the various chemical bonds that make up a molecule in the direction of the bond axis. Bond angle bending energy (E angle ): the change in molecular energy caused by a change in bond angle. Dihedral angular distortion energy (E dihedral ): change in energy caused by the distortion of the molecular skeleton caused by the rotation of a single bond. Cross energy term (E improper ): the change in energy caused by the coupling of the above actions. Non-bond interactions (E non-bond ): include non-bond interactions related to energy, such as van der Waals forces and electrostatic interactions.
In this study, to describe the covalent interactions between CNT and polymer atoms and the non-covalent interactions at the CNT-PE matrix interface, PCFF force field was adopted (as a member of the consistent family of force fields CFF91, PCFF, CFF and COMPASS [26]), which are closely related second-generation force fields). Thus, the total potential energy (E pot ) of the atomic system can be given as, where E bond , E angle , E dihedral and E improper are the bonded interaction energies corresponding to different deformations indicated before, and the E non-bond which represents the non-bonded interaction energy includes van der Waals (vdW) and electrostatic interactions in PCFF.
Since van der Waals interaction is defined just by Lennard-Jones (LJ) 9-6 atomic potential, only van der Waals non-bond interaction energy is considered. For the 9-6 function form in simulation, potential cutoff radius (i.e., r cutoff ) of 10 Å with the inter-atomic potential parameters is given in table 1.

Simulation procedure
Based on the initial structure of the unit cell, MD simulation of LAMMPS was carried out by using PCFF force field. The potential energy of the CNT-PE model (see figure 1(c)) was the minimized by the conjugate gradient method in simulation. It was annealed to 600 K for 100 ps and then quenching to 100 K in 200 ps using the Nose-Hoover thermostat by NVT ensemble. Finally the material system has balanced at 200 ps under the NVE ensemble. With a temperature of 100 K and a pressure of 1 atm in all directions, the volume and potential energy  Table 1. Lennard-Jones (LJ) 9-6 inter atomic potential parameters.
0.054 4.010 Note: C-C and H-H represent the interactions between carbon and hydrogen atoms, respectively, and subscripts sp 2 and sp 3 specify hybridization types.
of the material become stable and relaxed. At this status, the equilibrium density of the polymer group was around 0.93 g cm −3 , which was slightly lower than the experimental value of 0.95 g cm −3 [27]. The main reason lies the non-defective material and only the non-bonding energy at the interface were taken account in the simulation of force field. After the equilibration procedure, uniaxial tension was simulated by the NVT ensemble, as shown in figure 2. The atoms at the bottom side were fixed and the ones at the top side stretched with a constant velocity. In simulation we presupposed three velocities of 3×10 −4 Åps −1 , 3×10 −5 Åps −1 and 3×10 −6 Åps −1 , respectively. On this atomistic model through MD simulation, we obtain the complete the stress-strain curve. And on the data of this variation curve, we can further estimate the accurate elastic modulus of the nanocomposite. Usually in MD simulation, the stress was given by the formula as [28] å å

Results and discussion
On the atomistic model of CNT-reinforced polyethylene nanocomposites and through MD simulation, the aim of this work is to predict the elastic properties of the CNT-polyethylene nanocomposite and the CNTpolyethylene interface under different loading rates. Specifically, in this section, the elastic properties of the nanocomposite, the interface characterization and the comparative verification are respectively given as follows.  the same conditions. Besides, data in the table 2 shows that the mass density increases with loading rates (every 10 times increment), and the elastic modulus increases accordingly. The results indicate a reduction of loading rates actually lead to the PE matrix material to soften, or in other words, lowering loading rates is equivalent to raising the temperature. Meanwhile, a CNT-zigzag typed with chiral vector (9, 0)-having diameters of 7.05 Å was also determined through MD simulations. The Young's modulus of CNTs at 300 K was obtained as 1056.95Gpa, consistent with the results obtained by Singh et al. [21]. All above these results partially support the validation of our modelling and simulation.

Elastic properties of the nanocomposite
To investigate the elastic properties of CNT-reinforced polyethylene nanocomposites, the atomistic model of a CNT with chiral vector (9, 0) embedded into pure PE polymer matrix is considered (see figure 2). On the model, the simulated stress-strain curves under different loading rates were obtained as shown in figure 3. In each case of loading rates the calculated data points of the stress are not equal at all but stable almost within a small range, especially for larger loading rates (for example, the rate 3×10 −4 Åps −1 ). For these data points we made a curve fitting through the least square method. Then, based on the stress-strain curves, usually the Young's modulus of the nanocomposite can be obtained by estimating the slope of the linear elastic zone (3% of the strain) of the curve. Thus, Young's modulus of the nanocomposite from the stress-strain curves under three different loading rates were obtained as listed in table 2 (the last column). Similarly with the results of the mass density and Young's modulus of the PE matrix in table 2, and obviously in figure 3, the faster the loading rate is, the bigger the stress is ( figure 3) and also the higher the Young's modulus (table 2). In other words, increasing of loading rates leads to the matrix material to harden. These results are the same conclusion reached by Bao et al. [13]. In fact, these Polymeric materials are found to be highly sensitive to loading rates (or strain rates) [30]. But for specific different polymeric material, for example polyethylene and polypropylene, they have different dependence of loading rates on stress, strength and Young's modulus, and the deformation process and mechanism are quite different. Therefore, for CNT-reinforced polyethylene nanocomposites, to interpret these interesting results, deformation processes involved microstructure evolution has to be emphatically considered. Besides, In figure 3 we can also find that the stress is not 0 when the strain just appear. The reason is that the  constant motion and interaction of atoms lead to the stress inside the nanocomposite during the initial relaxation process.

Characterization of the CNT-polyethylene interface and its properties
As shown in figure 4, for the CNT-polyethylene nanocomposite, the interface denotes the region between the CNT and polyethylene matrix. In this section, in order to completely characterize the interface, we will investigate the interfacial interaction energy, elastic modulus, and the thickness of the interfacial region, respectively. Firstly, to featuring the interface, we can start at predicting the interfacial interaction energy between CNT and PE matrix. In general, the interfacial interaction energy (E int ) is mainly determined by the non-bonded interaction energy (i.e., here E vdW interaction energy only). It can be calculated by subtracting the sum of the potential energy of CNT (E CNT ) and polyethylene matrix (E pol ) from the total potential energy of nanocomposite (E nc ), as given as the formula, where ε is the applied strain. Usually the interfacial interaction energy is closely related to the feature of CNTs. Before calculation, we need to give a description of CNTs. In general, for CNTs, besides the aspect ratio, the specific surface area (SSA), as another characteristic parameter, was also used. It can be calculated by where L and r are the length and diameter of CNT, respectively. In calculation, we choose three CNTs with the aspect ratios of 4.83, 3.62 and 2.41, and accordingly the SSAs are 0.628, 0.611 and 0.601. Then, the interface interaction energy E int can be calculated under applied strain after each step of the respective applied constant velocity, and the result of E int with strain ε was obtained as shown in figure 5. On these data in figure 5, the correlation error estimation of the fitted curves was estimated by regression analysis, and the determination coefficients (R 2 ) of the different curves were obtained as follows: R 2 >0.99 for 3×10 −4 Å ps −1 , R 2 >0.97 for 3×10 −5 Å ps −1 , and R 2 >0.96 for 3×10 −6 Å ps −1 . It is found that, the interfacial interaction energy (E int ) gets larger with applied strain, and it is much higher instead when the loading rate is lower.
Next, the Young's modulus of the interface region can be calculated by the following equation as, where E int is the interaction energy value at the interface (see figure 5). V int is the volume of the interface region between CNT and PE matrix, as shown in figure 4, which can be calculated by the following formula as, where L is the length of CNT placed parallel to the x-axis (see figure 1(a)) in the CNT-PE nanocomposite. r CNT is the radius of CNT, and t CNT the vdW separation distance which is the thickness of the interface region. Usually the thickness of the interface region can be obtained from a radial distribution function (RDF) curve, and the RDF values of PE particles around CNT can be observed through VMD in periodic unit cell [27]. As shown in figure 6, the RDF value of PE particles keep unchanged and always at 0 within the range of the radial distance from 0 to 2.1 Å. In other words, the probability of PE particles appearing within the distance away from CNT in 2.1 Å is 0. Therefore, under this set of parameters for the atomistic nanocomposite model of CNT-polyethylene nanocomposites, we can take the characteristic value 2.1 Å as the thickness of the interface region. According to equations (5) and (7), Young's modulus of the interface, to a large extent, depends on the selected CNTs. On above the atomistic model ( figure 1(c)), we kept the same mass fraction PE matrix and just replaced the CNT with different aspect ratios, and obtained the Young's modulus of the CNT-PE interface under different loading rates for three CNTs with different aspect ratios as listed in table 3. For the CNT-polyethylene nanocomposite model with a given CNT (aspect ratio is fixed), the higher the loading rate is, the larger Young's modulus of the interface is (corresponding to the data in each column); under the same loading rates, Young's modulus of the interface increases with the aspect ratio of CNTs.
And for CNTs, different aspect ratio corresponds to different specific surface area (SSA). For the CNTpolyethylene nanocomposite with a CNT which has different aspect ratio and correspondingly different SSA, as a result, the interfacial interaction energy (E int ) has to be different. Similarly on above the atomistic model, the interfacial interaction energy for three CNTs with different SSAs were obtained as listed in table 4. It is noted that, for the characteristic parameters of CNTs, SSAs of 0.628, 0.611 and 0.601 are just corresponding to the  aspect ratios of 4.83, 3.62 and 2.14. Obviously, if ignore the minus of the calculated values of the interfacial interaction energy, we can find similar variation with loading rates and SSA of CNTs. Combining the data in tables 3 and 4, for the CNT-polyethylene nanocomposite including the CNT with the maximum aspect ratio, correspondingly the maximum SSA, its Young's modulus and interfacial interaction energy reach the largest under the fastest loading rate.
The simulated results show that the external loading rate has a high impact on the elastic modulus of the CNT-polyethylene nanocomposite and its interface. In fact, the loading rate itself does not affect the elastic modulus of the overall nanocomposite. However, due to different loading rates, energy changes in the interface region are caused, and energy changes lead to temperature changes in the interface and composite. And polymer-based nanocomposites and their interfaces are very sensitive to temperature changes, and the fast stretching is equivalent to temperature falling, while the slow stretching is equivalent to temperature rising.

Comparative verification
To verify the model and simulation procedure in this work, we calculated the stress-strain relation and Young's modulus of the CNT polymer nanocomposite, and compared them with experimental and theoretically predicted ones in previously published works.
For the stress-strain relation of the CNT polymer nanocomposite, since there was no relevant works of CNT-polyethylene nanocomposites in which sufficient parameters are available, instead we choose the measured data of CNT polypropylene (PP) nanocomposites. Here simulated the stress-strain relation based on the MD model and made a comparison with the experimental data by Bao et al [13]. As given in figure 7, under the two different loading rates, the MD simulated results shows a very good agreement with the experimental data.
And, for Young's moduli of CNTs and the Integral nanocomposite, we still hope make a comparison of our calculated results on MD simulation with the ones on other approach and other theoretical predicted ones, which is also very necessary to verify the model and simulation procedure. Firstly, we compared our simulated Young's modulus of the Integral CNT-polyethylene nanocomposite with the ones by using the continuumbased rule-of-mixtures (ROM). ROM is a method based on the theory of continuum mechanics, which assumes perfect bonding of three phases in composite materials. The mathematical expressions of Young's modulus of three phase nanocomposites are given by strain compatibility equation and profit equilibrium equation as where V pol , V CNT and V int represent the volume fraction of the polymer matrix, CNT and interface, respectively, and Y pol , Y CNT , Y int represent the elastic moduli of the polymer matrix, CNT and interface in the nanocomposite. In addition, the three phases volume fraction of the nanocomposite must satisfy the following relation Based on the micromechanical hypothesis, we modeled the hollow CNT as a solid cylinder and calculated its volume fraction. As reported by Rafiee and Madhavi [31], considering that the cross-sectional area of hollow CNT is equal to the one of solid CNT, the equivalent radius (r sCNT ) and volume fraction of solid CNT can be calculated by where t CNT and r CNT are respectively the wall-thickness and radius of the hollow CNT, and A nc is the crosssectional area of the CNT-PE model (30.48 Å×30.48 Å), and in fact the cross-section area of the final model after NVE equilibrium. In simulation, the sample was stretched under the loading rate of 3×10 −4 Å ps −1 for verification. As listed in in table 5, Young's modulus of the nanocomposite obtained on the continuum-based ROM at the loading rate of 3×10 −4 Å ps −1 are very close to the one simulated by MD method (the relative error ∼ just 8%). The calculation value on the ROM is slightly larger because perfect bonding of composite materials is assumed among three phases of the CNT, PE and interface in the nanocomposite. Therefore, this comparison of Young's modulus on MD simulation with that on the ROM method supports the atomistic modeling and simulation procedure developed in this work.
Further, we compared our simulated Young's modulus of the nanocomposite with other theoretical predicted one. In Ref. [15], for the CNT reinforced methyl methacrylate (PMMA) composite, on taking account of interface effect and in the cases of different the aspect ratio of CNTs, the Young's modulus varying with the aspect ratio of CNTs was given. Thus, based on the molecular unit cell model of PMMA matrix with a size of 8.0 Å×3.7 Å×3.7 Å nm 3 and a mass density of 1.17 g cm −3 , the whole material was stretched and the young's modulus of the material was statistically calculated under the condition of changing the aspect ratio of CNTs. Finally, as shown in figure 8, we obtained the Young's modulus vs. the aspect ratio of CNTs, and compared the result with the one in [15]. Obviously, varying with the aspect ratio of CNTs, our simulated Young's modulus  keeps almost consistent with the one obtained by Arash et al [15]. And, the result indicates that carbon nanotubes with high aspect ratio will greatly enhance the strength of the composite.

Summary and conclusions
The elastic modulus of the CNT-polyethylene nanocomposite and its interface region were simulated by MD method under PCFF field. We selected the CNT-zigzag type with the chiral vector (9, 0) as reinforcing fillers for the nanocomposite. The external loading rates are between 3×10 −6 Å ps −1 and 3×10 −4 Å ps −1 , and the elastic modulus of the nanocomposite and its interface were calculated. Under different loading rates, the interfacial interaction energies of the CNT-polyethylene nanocomposite are accordingly different, and further, with different aspect ratios of CNTs, the interface energies are also different, which indicate that the interface energy is related to the contact area of the two phases. And, through comparing with the ones that are on ROM method and from other published works, the atomistic model and simulation procedure are verified to be accurate and reliable. Based on the results and discussion, the remarkable points are concluded as, (1) For the CNT-polyethylene nanocomposite with CNTs of higher aspect ratios, it has larger CNTpolyethylene interface area-to-volume ratios, and correspondingly the role of the interface is getting much more prominent in dominating the mechanical properties. And consequently, the elastic moduli of the interface region and the whole nanocomposite material are getting greater.
(2) In the case of uniaxial tension, if the full potential of CNTs enhancers in the nanocomposite is wanted to be utilized, the selected CNTs with higher aspect ratios and volume fractions have to be considered.