Interfacial Characteristics of Carbon Nanotube-Polyethylene Composites Using Molecular Dynamics Simulations

1 Center of Micro/Nano Science and Technology, Jiangsu University, Zhenjiang 212013, China 2 Department of Engineering Mechanics and State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Liaoning, Dalian 116024, China 3 Radiation and Nucleation Detection Materials and Analysis Department, Sandia National Laboratories, Livermore, CA 94550, USA 4 Department of Mechanical and Aerospace Engineering, Utah State University, Logan, UT 84322, USA 5 Center for Advanced Vehicular Systems, Mississippi State University, Mississippi State, MS 39762, USA


Introduction
Carbon nanotubes (CNTs) have excellent properties, such as very high elastic modulus and tensile strength, which make them suitable for making advanced polymer composites [1].Recently, several research groups have fabricated and mechanically evaluated CNT-polymer nanocomposites, to explore their potential as advanced lightweight materials [2,3].These composites exhibited enhanced moduli, indicating that CNTs carry a portion of the applied load [2][3][4][5][6].The interfaces between polymer matrix and CNTs transfer load from one to the other and thus play important roles in the nanocomposite's stiffness and strength.At this point, the interfacial properties between CNTs and polymers have not yet been clearly defined for optimization of the composite's mechanical behavior.This paper will examine the CNT/polymer interfacial interactions to determine the dominant mechanisms responsible for load transfer using molecular dynamics (MD) simulations.
The interfacial shear strength, a key property responsible for maximizing the load transfer in short-fiber composites, has been experimentally evaluated using fiber pullout or push-out tests [7,8].Experiments have demonstrated that the load transfer mechanism in nanocomposites is weak due to the difficulties in making strong bonds between CNTs and polymer matrices [2,6,9,10].Thus, CNTs often interact primarily with the polymer matrix through van der Waals forces.The van der Waals forces especially dominate the load transfer capability of the CNT/polymer interface when cross-links between CNTs and polymer matrices are not present.Due to difficulties in devising experiments to study the CNT/polymer interface, molecular mechanics and molecular dynamics simulations have become increasingly popular in investigating the reinforcement mechanisms in nanocomposites [11,12].
MD simulations offer a method, similar to the CNT pullout from polymer matrix of the previously mentioned experiments, to probe the mechanism of interfacial failure and strengthening.Previous studies, using molecular mechanics simulations of CNT/polymer composite's interfacial binding, predicted that the critical interfacial shear stress (CISS) ranged from 18 to 135 MPa for (10,10) CNTs embedded in a polymer [13].The results showed that the binding energy and frictional force played only a minor role in determining the strength of the interface, but the helical polymer conformations in which chains can wrap around nanotubes may enhance nonbonded CNT-polymer interactions.Liao and Li [14] evaluated the interface of a CNT-reinforced polystyrene (PS) composite through molecular mechanics and elasticity calculations.Again, no chemical bonding was considered in the simulation; thus, the interfacial adhesion primarily resulted from van der Waals interaction, mismatch in the coefficients of thermal expansion, and radial deformation induced by atomic interactions.Using MD simulations, Frankland et al. [15] studied the influence of chemical cross-links between a single-walled CNT and a polymer matrix on the interfaces tensile and shear strength.The results suggested that load transfer and modulus of CNT/polymer composites can be effectively enhanced by deliberately increasing the cross-links.Based on this result, they suggested that the chemical bonding between CNTs and polymer matrices is partially responsible for the load transfer at the interface.Frankland and Harik analyzed the entire CNT pullout process by introducing the critical pullout force concept, with the CNT/polymer interfacial sliding described based upon the principle of the Newton's friction law [16].The results showed that the interfacial interactions during CNT sliding process could be approximated by a linear force-velocity relationship with the slope as the effective viscosity coefficient of the interface.Chowdhury and Okabe [17] also simulated a CNT pullout from a polymer using MD simulations, considering the influence of the polymer matrix density, chemical cross-links at the interface, and geometrical defects in CNTs.They showed that the case without cross-links between the CNT and matrix had an ISS of 1320 MPa which was lower than the case with cross links.Chowdhury and Okabe also showed that the pentagonheptagon geometrical defect of CNTs has little effect on the interfacial shear stress (ISS).Based on MD simulations, Liu et al. [18] proposed a new boundary element method (BEM) using a cohesive nonbond interface model for CNT/polymer composites.Results demonstrated the usefulness, efficiency, and promise of the developed BEM as a fast numerical tool to characterize the CNT/polymer composites at larger scales.Gou et al. [19] predicted the interfacial bonding of single-walled CNT reinforced epoxy composites using MD simulation based on a cured epoxy resin model.The pullout simulations indicated an effective stress transfer from epoxy resin to the nanotube, which was comparable to their experimental results.Wei [20] studied the temperaturedependent adhesion behavior of CNT/polymer composites using MD simulations and discovered that the ISS resulting from van der Waals interactions increased linearly with the applied tensile strains along the nanotube axis direction and decreased as the temperature increased.
As mentioned above, the interface of CNT/polymer composites has been studied by many research groups, considering the morphology of polymer, the temperature of the simulation systems, the cross-links between the CNTs and polymer, and the defects of the CNTs, using MD simulations, molecular mechanics, and some continuum methods.The contribution of this paper extends the MD simulations on CNT/polymer composites to investigate loading rates and size effects, including the simulation size, the chain length of polymers, and the radius of an embedded CNT.

Computational Method and Model
This work used a classical MD model where the Newtonian equations of motion are solved using the velocity-Verlet algorithm at a time step of 10 −16 s.The MD simulations were conducted using the software package LAMMPS [21].The Brenner potential [22] was used for carbon-carbon interactions (or bond) for the CNTs.The modified AMBER potential functions [23] were used to describe the interactions between adjacent PE molecules, with each molecule represented using a united-atom.The total potential energy of the molecular system, E tot , includes bond stretch (bs), bending (be), torsion (to), and the nonbonding potential (van der Waals, vdW) between the united atoms as given by the following: where E bs is the two-body potential for bond length r, r 0 = 0.1533 nm is the equilibrium bond length, and k r is the bond stretching stiffness.E be is the three-body potential for bond angle θ, θ 0 = 109.5interfacial interaction between the CNT and the polymer matrix, ε = 0.4492 kJ/mol and σ = 0.401 nm [24].
A single-walled CNT embedded into an amorphous polyethylene matrix comprised the unit cell of the MD simulation.The CNT spanned the total length of the unit cell, as shown in Figure 1.Periodic boundaries were applied in x and y directions.Two edges of the unit cell in the z direction were fixed to prevent the rigid body motion of the PE, shown as the red boxes in Figure 1.The CNT/PE composites were created according to the following procedure.
(1) Initially, create the PE amorphous structure; the coordinates of chains united atoms are randomly generated on an FCC lattice, using a random walk method similar to that of Shepherd [25].
(2) Obtain the equilibrium configuration of amorphous PE; all the chains are equilibrated for 10 6 MD time steps in the isothermal-isobaric (NPT) ensemble at 300 K to allow for energy relaxation and sufficient entanglement.
(3) Generate a cylindrical pore; a cylinder indenter is placed in the center of the simulation box which exerts a force of magnitude F(r) = −k(r − R) 2 on each united atom, where k is the specified force constant, r is the distance between the united atoms and the cylinder axis x, R = R CNT +R cutoff , R CNT is the radius of the CNT, and R cutoff = 2.5 σ is the van der Waals cutoff distance.
(4) Embed an armchair CNT into the cylinder pore; equilibrate the CNT/PE composite for 500,000 MD time steps at 300 K using NPT dynamics to create a zero initial stress state.
Upon reaching an equilibrium configuration for the CNT/PE composites, the temperature of the composite is cooled down from 300 K to ∼0 K to avoid any thermal effects of MD simulations during the sliding process [17].Some small initial stresses, which are dependent on the initial distribution of the PE molecules, may be generated at the interface during this annealing process.After cooling an incremental displacement ranging from 10 −6 to 10 −5 nm, which determines the CNT sliding velocity ranging from 10 m/s to 100 m/s, is applied to the CNT atoms to simulate the pullout.In all simulations the CNTs are treated as rigid fibers due to their exceptionally high stiffness compared with PE matrices [18].

Energy Change during the CNT Sliding Process.
The ISS is first studied by examining the energetic changes during the sliding of an armchair (10, 10) single-walled CNT (SWCNT) embedded in a block of amorphous PE consisting of 100 chains with 200 monomers/chain (mpc).ISS is calculated using τ I = F CNT /2πrL, where F CNT is the CNT-axial component of the total force on CNT due to the pairwise interactions of PE molecules, r is the CNT radius, and L is the CNT embedded length.In the calculation, the total force on CNT is the sum of the force interaction between the group of carbon atoms of CNT and the group of PE molecules.In MD, forces acting upon atoms are derived from the potential energy which depends on the particles' coordinates.Figure 2 shows ISS as a function of CNT displacement for different sliding velocities.Three distinct stages were observed during the sliding process: (1) a linear increase, which may represent an "elastic" response, before peaking at a possible "yield" or critical stress point, (2) slippage that begins with a transition region where the ISS decreases as the displacement increases, and finally (3) the steady sliding stage where the ISS remains the same as the displacement increases.In the "elastic" regime, the ISS increases almost linearly with displacement until reaching a maximum value.We define this maximum value as the critical interfacial shear stress (CISS).Upon reaching the CISS, the ISS drops sharply suggestive of sliding.After the transition region, the ISS eventually remains constant, which we refer to as a steady sliding stage.During steady sliding stage, ISS oscillates with a period approximately equal to the length of the CNT's unit cell along the CNT axis (2.5 Å), as shown in Figure 2. The results show that the periodic structure of the CNT surface results in an approximate stickslip sliding [26], with the bonding of the matrix jumping from one unit cell of the CNT to the next.The steady sliding interfacial shear stress (SISS) is obtained by averaging the ISS over the whole steady sliding time steps.The ISS increases as the CNT sliding velocity increases, as shown in Figure 2. Further, the CISS is linearly dependent on the sliding velocity of CNTs (see Figures 4(a Figure 3 shows each component of the potential energy (See the symbols of Figure 3, which are defined in (1) as a function of the CNT sliding displacement with a pulling velocity of 10 m/s.One obvious observation is that only the interfacial energy, E int , varies while tension, bending, torsional, and van der Waals energy do not change upon applying a load.This indicates that the interfacial energy dominates the composite's response, which is described strictly by nonbonded van der Waals interactions between  the PE chains and the CNT.The interfacial energy E int increases linearly and reaches a plateaus at the sliding displacement until 3 KCal/mol.

Effect of the PE Chain Length.
To investigate the effect of PE chain length on ISS, the PE matrices were built with chain lengths of 100, 200, or 300 molecules/chain (mpc), respectively.An armchair (10, 10) SWCNT was embedded and pulled at velocities ranging from 10 ∼ 100 m/s.For each case, in an attempt to remove any effects from the simulation size, the total number of monomers is held constant at 100,000.Also, to ensure the results are not dependent on initial PE velocity profiles, all reported CISS and SISS values are averages of three individual simulations with different initial PE velocity trajectories.Figure 4 shows that both the CISS and SISS increase as the chain length increases.One factor that could increase the CISS is the stiffness of the material.Longer chains tend to lead to a greater number of entanglements in the matrix that could result in a stiffer material.Moreover, the CISS increases linearly with the relative sliding velocity of CNTs, while the relationship between SISS and the CNT-sliding velocity approximately (see Figures 4-6) yields a linear relationship when the SISS is below a critical value.As for the CISS, the shear strain rate of the PE is proportional to the CNT velocity due to the noslip reinforcement at the PE-CNT interface.In contrast, the SISS is obtained from the steady sliding stage.

Effect of the Simulation Size.
The MD simulation cell size effects on the predicted material response are demonstrated from two aspects: the number of monomers and the monomer chain lengths.To capture the effect of the simulation size on the ISS, two simulation sizes, 20,000 and 100,000 monomers, are simulated for chain lengths of 100, 200, and 300 mpc, with the configurations shown in Table 1.
The sliding of an armchair (10, 10) SWCNT is examined for velocities ranging from 10 m/s to 100 m/s. Figure 5 shows the CISS and SISS for all cases as a function of the CNT sliding velocity.For the cases of 200 and 300 mpc, enlarging the simulation size from 100 chains to 500 chains clearly increases the CISS while the SISS changes only slightly (see Figures 5(c) and 5(d)).In contrast, for the simulations of chains with 100 mpc, the CISS decreased when the number of chains increased from 200 to 1000 (see Figure 5(a)).In short, for increasing simulation size, the CISS decreases for shorter chain lengths but increases for longer chains.In contrast, the SISS is not sensitive to the simulation size for these configurations (see Figures 5(b), 5(d), and 5(f)).

3.4.
Effect of the CNT Radius.Armchair (5, 5), (8,8), (10,10), and (15, 15) SWCNTs, with corresponding radii of 0.346, 0.554, 0.692, and 1.038 nm, respectively, were simulated to study the role of CNT radii on interfacial properties.Figure 6 shows that the SISS, as a function of the sliding velocity for different radii of CNTs, decreased as the CNT radius increased.For higher sliding velocities, the ISS plateaus and the differences in ISS due to CNT radius are evident.Interestingly, in the high sliding velocity region, the distance between the curves of (5,5) and (8,8) CNTs is larger than that between the curves of (8, 8) and (10, 10) CNTs, while there is only a little difference between the curves of (10, 10) and (15,15) CNTs.It implies that for an increasing CNT radius, there is a lower limit for the SISS for a curvature approaching zero (a flat graphite sheet).To show the lower limit the ISS is calculated for a graphite sheet embedded in a PE matrix and displaced at the same range of constant velocities as above.The results again show that the SISS increases as the sliding velocity increases.It indicates that the SISS in the case of graphite sheet is a lower bound value (1.1 MPa) in the high sliding velocity region (see Figure 6(a)).Furthermore, the effect of the radius of CNTs on CISS is investigated, as shown in Figure 6(b).The CISS also decreases as the CNT radius increases, and the lower bound value of CISS (2.4 MPa) is obtained from the graphite sheet sliding process.

Discussions
The interfacial behavior of CNT-polymer composites, such as the SISS and CISS, is affected significantly by the MD simulation model.The fundamental influential factors for the interfacial behavior are studied in the following.

Effect of Polymer Chain
Length.Supplemental tensile simulations of bulk polymer with chain lengths of 100, 200, and 300 mpc have been performed to determine the matrix elastic properties for a loading strain rate of 10 8 /s.The results show that the stiffness of the PE increases significantly as the chain length increases from 100 to 200 mpc while there is a little change when the chain length increases from 200 to 300 mpc, as shown in Figure 7.This is in agreement with previous work for similar MD simulations of PE [27].The increases in stiffness parallel the results of chain length effects on the CISS, which increases with increasing chain length from 100 to 200 mpc while there is a little increase from 200 to 300 mpc (see Figure 4).This suggests that the increase of the ISS with increasing chain length is a result of the increasing PE stiffness.

Effect of Polymer Molecular Density at the Interface.
One method for explaining the increase in shear stress is to examine the proximity of PE molecules to the surface of the CNT, which can be described using a radial distribution function (RDF) g(r).From the MD simulations, the expression for evaluating g(r) from the output data is given as follows: where N is the total number of atoms in the simulation system, ρ = N/V the atom number density with V being the representative volume, N i the number of atoms found in a spherical shell of radius r and thickness Δr, with the cell centered on atom I, and V (r, Δr) the volume of the spherical shell.In Figure 8, the slight increase in the RDF indicates the presence of more PE atoms at the interface for the critical point of CNT slippage relative to the relaxed state.This increased density at the critical point results in a higher interfacial energy thus, facilitating a larger critical ISS than the sliding ISS (see Figure 2) during the CNT-pulling process.

Effect of the Simulation Model Sizes.
There are two factors influencing the CISS as the simulation size changes.One factor is the stiffness of PE and the other is the fixed  boundaries used to prevent the PE chains from moving with CNTs in our simulation model (see Figure 1).A series of tension tests are again carried out at a strain rate of 10 8 /s to investigate the effect of the simulation size on the stiffness of PE.A set of stress-strain curves is obtained for each chain length and simulation size and shown in Figure 9.For small chain length, 100 mpc, enlarging the simulation size has almost no effect on the elastic modulus of PE (see Figure 9(a)), while there is a prominent increase for the chain length of >200 mpc (see Figures 9(b) and 9(c)).
The second effect due to the fixed boundary is caused by a localization of the polymer deformation near the nanotube surface.The displacements of PE molecules along the CNT sliding direction versus the distance from the tube wall to the fixed boundary at the point of CISS are plotted in Figure 10.For chains length of 100 mpc (See Figure 10(a)), the displacement of the PE molecules at the interface is smaller for the simulation size.This decrease in displacement implies that the fix boundary effect dominants for the shorter chain length cases due to the increase of the distance from the fixed boundary to the tube wall.Oppositely, for the chain length of 300 mpc (see Figure 10(c)), the displacement of the interfacial PE molecules increases with increasing the number of chains from 67 to 333, which reveals that the increase of the stiffness of PE dominates for the longer chain length cases, resulting in an increase in CISS with increasing simulation size.The chain length of 200 mpc (see Figure 10(b)) is similar to the 300 mpc case.

Effect of CNT Radii.
To study the various mechanisms of the SISS with an increasing CNT radius, we trace back to the calculation of ISS, τ ISS = F CNT /A int , where F CNT is the total force upon a CNT A int is the interface area.For one carbon atom of a CNT, the occupied area is a constant value, which is independent of the CNT radius.Thus, the ISS is directly proportional to the force upon a carbon atom of CNT.Moreover, the force upon a carbon atom is directly proportional to the interaction energy, when it experiences the same displacement.Thus, the ISS is directly proportional to the interaction energy for one carbon atom of CNT.As shown schematically in Figure 11, the PE region of interaction with a single appointed carbon atom decreases with enlarging the radius of CNTs.Hence, the ISS decreases as the CNT radius increases, or the ratio of volume of interaction versus CNT surface area decreases as the CNT radius increases.The reason the graphite sheet sets the lower limit for the ISS is that the PE region of interaction for a single carbon atom is smaller than any other radius.

Summary
MD simulations are used to study the interfacial characteristics of CNT/PE composites using a nonbonded interfacial model.The interfacial behaviors between a CNT and a PE matrix are described as a function of the MD simulation model size, sliding velocity of CNTs, radius of CNTs, and polymer chain length.The CISS exhibits a linear dependence with the sliding velocity of the CNTs.The SISS is not sensitive to the simulation size while the CISS is significantly affected by the simulation size, resulting from a trade-off between the proximity of the rigid boundary and the stiffness of PE.Both the CISS and SISS decrease as the CNT radius increases but increase as the polymer chain length increases.

Figure 1 :
Figure 1: Atomistic model of CNT/PE composite.The red boxes show the fixed edges.

Figure 2 :
Figure 2: The interfacial shear stress (ISS) as functions of the sliding displacement of (10, 10) SWCNT at different CNT velocities.

Figure 3 :
Figure 3: Change of the potential energy as a function of the displacement of SWCNT10, 10 showing that the interfacial energy, E int , changes the most.The sliding velocity of CNT is 10 m/s.

Figure 5 :Figure 6 :
Figure 5: The ISS as a function of the sliding velocity of CNT (10, 10) for the different simulation size.

Figure 7 :
Figure 7:  The comparison between the stress-strain behaviors of PE systems consisting of 1000 chains each with 100 monomers, 500 chains each with 200 monomers, and 333 chains each with 300 monomers.The strain loading rate is 10 8 /s.

Figure 8 :
Figure 8: The radial distribution function (RDF) in the case of armchair (10, 10) CNT sliding in the PE with 100 chains, 200 monomers per chain.The sliding velocity of CNT is 10 m/s.

Figure 9 :
Figure 9: Comparisons between the stress-strain behaviors of PE systems consisting of (a) 200 and 1000 chains each with 100 monomers, (b) 100 and 500 chains each with 200 monomers, and (c) 67 and 333 chains each with 300 monomers.The strain loading rate is 10 8 /s.

Figure 10 :Figure 11 :
Figure 10: Comparisons between the displacements of PE molecules each with (a) 100 monomers, (b) 200 monomers, and (c) 300 monomers in the CNT sliding direction at different simulation sizes.

Table 1 :
Simulation sizes for different configurations.l x , l y and l z are the length of the simulation box in x, y and z directions.