Abstract

The rate-dependent interfacial behavior between a carbon nanotube (CNT) and a polyethylene (PE) matrix is investigated using molecular dynamics (MD) simulations. Various MD simulations were set up to determine the “size” effects on the interfacial properties, such as the molecular weight, or the length of the polymer, the diameter of the CNT, and the simulation model size. The interfacial rate-dependency was probed by applying various relative sliding velocities between the CNT and the polymer. Two quantities, directly obtained from the MD simulations, described the interfacial properties: the critical interfacial shear stress (CISS) and the steady interfacial shear stress (SISS). The simulations show that the SISS was not sensitive to the simulation size. In addition, the CISS was dependent upon the combined factors of the variation in PE stiffness, induced by simulation size changes and the effect of the fixed boundaries of the simulation models. The CISS increases almost linearly with the relative sliding velocity of CNTs. Also, a linear relationship between the SISS and the CNT-sliding velocity is observed when the SISS drops below a critical value. A clear size scaling is observed as the CISS and SISS decrease with increasing CNT radius and increase with the increasing polymer chain length.

1. 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 [26]. 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 pentagon-heptagon 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 temperature-dependent 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.

2. 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, 𝐸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: 𝐸tot=𝐸bs(𝑟)+𝐸be(𝜃)+𝐸to(𝜙)+𝐸vdW𝑟,𝐸bs(𝑟)=atom𝑘𝑟𝑟𝑟0,𝐸be(𝜃)=atom𝑘𝜃𝜃𝜃0,𝐸to(𝜙)=atom𝑉1cos𝜙+𝑉2cos2𝜙+𝑉3cos3𝜙+𝑉6,𝐸cos6𝜙vdW𝑟𝜎=4𝜀𝑟12𝜎𝑟6,(1) where 𝐸bs is the two-body potential for bond length 𝑟,𝑟0=0.1533 nm is the equilibrium bond length, and 𝑘𝑟 is the bond stretching stiffness. 𝐸be is the three-body potential for bond angle𝜃, 𝜃0=109.5 is the equilibrium bond angle, and 𝑘𝜃=56.25 kcal/mol is the bond bending stiffness. 𝐸to is the four-body potential for dihedral angle 𝜙. 𝐸vdW is another two-body potential for nonbonded distance 𝑟, 𝜀 is the potential well depth, and 𝜎 is distance where the van der Waals energy reaches zero. The van der Waals potential energy 𝐸vdW=𝐸vdW_PE+𝐸int, where 𝐸vdW_PE represents the van der Waals interaction between the PE chains, 𝐸int represents the interfacial interaction between the PE chains and the CNT carbon atoms. The potential functions and parameters for the PE are from [17], and, for the nonbonded interfacial interaction between the CNT and the polymer matrix, 𝜀=0.4492 kJ/mol and 𝜎=0.401nm [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 𝑥 and 𝑦 directions. Two edges of the unit cell in the 𝑧 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 106 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 𝐹(𝑟)=𝑘(𝑟𝑅)2on each united atom, where 𝑘 is the specified force constant, 𝑟 is the distance between the united atoms and the cylinder axis 𝑥, 𝑅=𝑅CNT+𝑅cuto,𝑅CNT is the radius of the CNT, and 𝑅cuto=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].

3. Simulation Results

3.1. 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=𝐹CNT/2𝜋𝑟𝐿, where 𝐹CNT is the CNT-axial component of the total force on CNT due to the pairwise interactions of PE molecules, 𝑟 is the CNT radius, and 𝐿 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 stick-slip 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), 5(a), 5(c), 5(e), and 6(b)) while SISS nonlinearly increases with the increase of the CNT’s sliding velocity (see Figures 4(b), 5(b), 5(d), 5(f), and 6(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, 𝐸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 𝐸int increases linearly and reaches a plateaus at the sliding displacement until 3 KCal/mol.

3.2. 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 46) 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 no-slip reinforcement at the PE-CNT interface. In contrast, the SISS is obtained from the steady sliding stage.

3.3. 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.

4. 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.

4.1. 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 108/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.

4.2. 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) 𝑔(𝑟). From the MD simulations, the expression for evaluating 𝑔(𝑟) from the output data is given as follows: 𝑔(𝑟)=𝑁𝑖=1𝑁𝑖(𝑟,Δ𝑟)𝑁𝜌𝑉(𝑟,Δ𝑟),(2) where 𝑁 is the total number of atoms in the simulation system, 𝜌=𝑁/𝑉the atom number density with 𝑉 being the representative volume, 𝑁𝑖 the number of atoms found in a spherical shell of radius 𝑟 and thickness Δ𝑟, with the cell centered on atom 𝐼, and 𝑉(𝑟,Δ𝑟) 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.

4.3. 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 108/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.

4.4. 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=𝐹CNT/𝐴int, where 𝐹CNT is the total force upon a CNT 𝐴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.

5. 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.

Acknowledgment

This paper was conducted collaboratively when the first three authors worked at the Center for Advanced Vehicular Systems (CAVS), Mississippi State University and thus would like to thank CAVS for the support.