Interfacial characteristics of a carbon nanotube-polyimide nanocomposite by molecular dynamics simulation

With molecular dynamics simulations, nanocomposites were characterized that are comprised of a polyimide (PI) polymer and carbon nanotubes (CNTs) with the same outer diameter but with one, two, or three walls. The simulations indicate that the PI/CNT interaction is strong, regardless of the number of CNT walls, and that there is some degree of alignment of the PI chains near the CNT interface. As the number of CNT walls increased, the density of PI chains near the CNT interface also increased and the average radius of gyration of the PI chains decreased, and these observations were attributed to changes due to the intertube van derWaals interactions. From simulations of the constant force pullout process of the CNT from the PI matrix, the limiting pullout force was calculated to be higher for the triple-walled CNT than for the single-walled one. The interfacial shear strength of the nanocomposites was also calculated from the pullout energy, and the results indicate that increasing the number of walls is a critical factor for enhancing the interfacial stress transfer during tension.


Introduction
Carbon nanotubes (CNTs), due to their high aspect ratios, have been regarded as one of the most promising reinforcement material in polymer-based nanocomposites [1][2][3][4][5][6]. Using pullout experiments with atomic force microscopy, the interfacial shear strength of CNTs within polymer matrices have been measured in the range of 35 to 380 MPa [7,8]. For systems without chemical bonding between CNTs and the matrix, the molecular-level physical interactions have been suggested to lead to nano-mechanical interlocking [9]. The existence of multiple walls in the CNT reinforcement can also impact the properties of the nanocomposite; for example, Raman spectroscopy was used to identify that only outer layers of MWCNTs are stressed in tension whereas all of the layers respond in compression [10]. Therefore, the wall number of CNTs plays an important role in stress transfer when interaction happens between CNTs and matrix under different types of loading.
Since the nanocomposite properties are closely related to the interfacial characteristics between CNT and the polymer matrix, researches have been done on the efficiency of interfacial stress transfer [11][12][13][14]. The stress transfer efficiency between a multi-wall CNT (MWCNT) and matrix was estimated through fragmentation method to be at least one order of magnitude larger in nanocomposite than in conventional fiber-reinforced composites [15]. MWCNTs were determined with SEM to provide a much longer pullout length and toughening as compared to carbon nanofibers [16]. Unfortunately, due to the difficulty in measuring nanoscale forces and strains in composites, large disparities have been observed in the experimental methods [17]. Therefore, computational simulations are drawing more attention in providing some crucial insights into the interfacial behavior of CNT/polymer nanocomposites.
Molecular dynamics (MD) simulations have been utilized to investigate CNT/polymer interactions, and are an effective approach because they can eludicate molecularlevel details that are challenging to do experimentally. In previous work [18][19][20], the interactions between CNTs and polymers with different CNT diameter, chirality, orientation, chain rigidity, chain length and functional groups were predicted with MD simulations [21][22][23]. MD simulations were also used to predict the degree of structural ordering within the polymer matrix, specifically the alignment of the chains parallel to the CNT axis, thus impacting the mechanical modulus of the material but not large scale ordering within the bulk [24].
Although the mechanical properties of polymer nanocomposites are dominated mainly by dispersion, alignment, and interface between CNT and the matrix, the interface is thought to be the dominant factor [25,26]. Recent simulations suggest that the small intertube spacing in MWCNTs can provide an effective means of load transfer between CNTs [27], although the pullout force for both double-walled CNTs (DWCNTs) and MWCNTs with more than two walls are found to be proportional to the diameter of the wall at the sliding surface [28], and independent of the length and chirality of the CNT. The interfacial shear strength (IFSS) between inner and outer CNT walls was estimated to be around 100 MPa from simulations [30]. In addition, binding energies and sliding frictional stresses between pristine CNTs and a range of polymer substrates have been estimated from molecular mechanics calculations [31]. Other significant factors appear to be the end cap of CNT, its diameter, and the strength of van der Waals interactions [32]. Recent MD simulations suggested that the interfacial energy and IFSS can be tuned through surface modifications of the CNT, especially hydroxyl groups at the ends [33][34][35][36].
Polymers that contain aromatic rings in the backbone tend to give better mechanical, thermal properties and good environmental resistance, especially for polyimides (PIs). However, CNTs/polyimide nanocomposites exhibit unsuitable mechanical strength [37], which is an issue for its desired applications. Thus, investigating the interaction mechanism between polyimide and CNTs is a way to figure out the phase response when nanocomposite is subjected to external force, and it is beneficial for property enhancement design. In this work, we present the results of MD simulations on the interfacial interactions of PI/CNT nanocomposite, and the effect of CNT wall number on interfacial characteristics has been analyzed. The specific PI presented in this report is 3,3',4,4'-biphenyltetracarboxylic dianhydride, p-phenylenediamine (BPDA/PDA). This PI was chosen because it has desirable mechanical and thermal properties, and thus the MD simulations will provide details about the interface between the CNTs and PI matrix, which can be used to guide the design of nanocomposites with desired characteristics through structural characterization and simulated pull-out studies.

Methodology
To investigate the interfacial properties of a nanocomposite comprised of a CNT and BPDA/PDA polyimide polymers, MD simulations are carried out on the CNT-polymer interactions as well as a simulation of the CNT pullout process under constant externally applied force conditions. We varied the number of CNT walls as single-walled (SWCNT), double-walled (DWCNT), and triple-walled", each depicted in Figure 1 in order to observe the effect of the wall number on the properties of the system.

Molecular models
The BPDA/PDA PI and the CNT structure were built using Accelrys' Materials Studio version 8.0. The molecular models of different type of CNTs are shown in Figure 1(a-c). The basic repeat unit and single chain of the BPDA/PDA PI is illustrated in Figure 1(d). For the PI, we chose chains that are 20 repeat units long to control the number average molecular weight to 7300, which not only reflects the real value [38] but provides control over the simulation box size. The length of PI chains is 40.6 nm. In addition, the wall separation of the CNT models was set to 0.34 nm, which is close to the interlayer distance in graphite. The longitudinal axis of each CNT was set spatially along the Z-axis to simplify the computational analysis. Open-end CNTs were used in all simulations.
For the multi-chain simulations, which mimic an actual nanocomposite, the CNT chiralities for each wall are as follows: (13,13) SWCNT, (8,7)(13,13) DWCNT, and (2,2)(8,7)(13,13) TWCNT. All three of the CNTs had an outermost diameter of 17.6 Å. The CNT was put into the center of a cubic box and longitudinally aligned along the zdirection, then 36 PI chains were placed parallel to the CNT. The initial configurations of the polymer chains were chosen to be relatively linear and longitudinally aligned with the CNT since the BPDA/PDA PI is classified as semicrystalline and previous simulations indicated that structural ordering occurs at CNT interfaces [24,39]. All three of the CNT systems had the same initial configuration of the polymers prior to equilibration. For the equilibration process described below, we made the simulation box large enough initially so that the density at first is low, thus the minimization and equilibration process could take up a larger portion of the energy space. Different starting configurations were sampled, and we used the lowest energy structure that was produced after minimization and equilibration.
For the simulation of a single chain of the PI with a TWCNT, we did the same procedure as earlier work [40,41] where the initial configuration of the CNT was along the Z-axis with the relaxed polymer chain perpendicular to it; the parallel conformation was also tried and no significant difference was observed. The distance between the CNT and polymer chain initially was around 20 Å, implying that the polymer chain was outside of the cutoff radius (10.0 Å) of interaction at the initial stage, to give the polymer as much conformational freedom as possible prior to being influenced by the TWCNT. The TWCNT was of type (2,2)(8,7)(13,13) with a diameter of 17.6 Å and a total length of 245.9Å.

Molecular dynamics (MD) simulations
MD simulations were performed with the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS, at the "High Performance Computing Center at North Carolina State University") [42]. All simulations used the pcff force field [43] which has been previously validated for simulating π-π interactions and mechanical properties between aromatic portions of a polymer (including PIs) and a CNT surface [44][45][46][47]. For the nonbonded terms, pcff uses the Coulomb potential for electrostatic interactions, and the Lennard Jones 6-9 potential for the van der Waals interactions, which have been validated previously for adequately representing aromatic interactions and interlayer potentials for graphitic compounds [48,49].
The multi-chain simulations were performed using periodic boundary conditions and with an isobaric and isothermal (NPT) ensemble with a pressure of 3 atm and temperature of 300 K for 3.5 ns with a simulation time step of 1 fs. The purpose of this process was to gradually compress the structure with low residual stresses in order to obtain an equilibrated structure. The applied pressure was optimized by conducting simulations at different pressures, and 3 atm was chosen as it resulted in the most energetically favorable system with the least amount of voids. Equilibration was defined as when the coefficient of variance of the last 50000 steps was controlled within 0.5%. The final box dimensions are 8.7×6.4×60.7 nm 3 for the SWCNT, 8.7×6.4×60.8 nm 3 for the DWCNT, and 8.8×6.4×61.2 nm 3 for the TWCNT. For the single chain simulation, we used an ensemble with a constant number of molecules, constant volume, and constant temperature (NVT) at 300 K with cubic periodic boundary conditions at 1 atm pressure, using the Nose-Hoover temperature constraints. We used a time step of 1 fs and a total equilibration time of 1 ns.

Simulations of CNT pullout
The pullout simulations of CNT from PI matrix were carried out on each equilibrated structure by adapting the procedure developed previously by others [50,51]. We applied a force-controlled load (using forces ranging from 6.94 fN to 69.4 fN) at every time step to all CNT atoms during NVT (isochoric and isothermal) MD simulations at 300 K, where the boundary was set as periodic in the x and y directions, and no boundary was set in the z direction since the forcecontrolled load was applied along the z direction; the simulation box in the z direction was extended 24.5 nm to accommodate the pullout process. The total time was set to 0.5 ns to accommodate the complete pullout process for systems where the force was sufficient to do so.

Simulation analysis
In an MD simulation, two kinds of atomistic interactions are taken into consideration, the bonded interactions from covalent bonds (E bonded ), and the nonbonded (E nonbonded ) interactions, including van der Waals (E vdW ) and electrostatics (E elec ). Thus, the total potential energy of the system (E system ) is given by the following equation, The interaction energy (E inter ) can then be defined as, where E CNT is the potential energy of the CNT, and E PI is the potential energy of PI, which were each calculated by running a single step MD simulation of each isolated system at 300 K. Note that because the bonded contributions will be the same in both the total system and the isolated components, E inter only takes into account the nonbonded contributions (E vdw and E elec ).

Equilibrated systems
Snapshots from the final equilibrated multi-chain systems are given in Figure 2 for each type of CNT, and reveal that the PI polymers both completely cover both the axial and the transverse directions of the CNTs. Figure 2 also provides density plots of the PI relative to the transverse direction of the CNTs within the simulation box, where the dips in the middle of the box (around x = 0) indicate the position of CNT. An interesting observation is the higher density closer to each CNT surface for all three types at around 1.2 nm, which is about 0.5 nm from the CNT surface, and that this density near the CNT interface increases in magnitude as the number of walls increases. In our previous work [47], we also observed density peaks for PI near the TWCNT interface, and its intensity varied depending on the volume fraction of TWCNTs within the system. Portions of the PI chains that are within 0.7 nm of the CNT surface are also indicated at the bottom of Figure 2, and reveal a few interesting observations at the interface: (1) portions of the PI chains are aligned with the CNT surface in various directions; (2) there is some degree of wrapping of PI chains, although not complete wrapping for most of the chains due to both the overall length of the polymer chains and due to confinement effects; and (3) some interchain alignment also occurs due to π-π interactions. In Figure 2(b), the lowest value showed around 2nm for DWCNT rather than TWCNT is attributed to the lower interaction energy caused by fewer walls, reflecting that the polyimide chains are less easily to be gathered or aligned around CNTs. Figure 4 presents histograms of the radius of gyration (Rg) for the PI chains, and distinct differences are observed depending on the type of CNT, where the average Rg is highest for the SWCNT system (7.2 nm) and lowest for the TWCNT system (5.8 nm). Since all three systems started with the same initial configuration prior to the equilibration process and the outer CNT wall for all three systems has the same diameter, this observation indicates that the differences in the van der Waals forces that result from multiple walls can have a significant impact on the structure of the PI matrix.  As a comparison, we also investigated the PI chain behavior relative to the CNTs if it were not confined by other PI polymers, and Figure 5 provides snapshots at specific timesteps along the trajectory for the single PI chain interacting with a TWCNT. The polymer evolves from first being linear in the initial timestep, and then once it "finds" the TWCNT, a portion of the chain begins to adsorb to the surface of the TWCNT, which causes almost instantaneous, helical-like wrapping of the rest of the chain. At later time steps, the overall polymer conformation remains helicallike, but its overall conformation is distorted. The entire length of the polymer chain is also observed to actively interact with the TWCNT molecule. This CNT-wrapping behavior was also observed with a SWCNT and a DWCNT as well as in our previous studies on stiff backbone polymer chains [41] and by others [52][53][54]. The average interaction energy between BPDA/PDA and the TWCNT and DWCNT from these MD simulations was calculated to be 0.38-0.40 kJ/kg, whereas it is 0.1 kJ/kg less for the SWCNT. However, rarely do polymers have this much conformational freedom unless at low concentration in solution, and thus the conformations are going to be influenced not only by the CNT but also by the presence of confinement by other polymer chains, especially in systems with significant entanglement and/or strong interchain π-π interactions.

Pullout simulations
During the process of pulling the CNT out from the polymer matrix, three different phases, depicted in Figure 6, can occur: Phase I: The CNT starts to be pulled out due to the externally applied constant force, and the polymer chains remain relatively intact.
Phase II: While the CNT is still being pulled out, polymer chains near the interface react by moving with the CNT, likely due to a combination of van der Waals forces and a change in density of the system due to the newly created void.
Phase III: The CNT is completely removed from the polymer matrix, and the polymer chains react by rearranging themselves in response to the void left by the CNT.
Note that some previous work [28] did not observe Phase II because the polymer positions were held fixed in order to reduce the computational cost. However, although the system is in the solid phase, some polymer motions would be expected at the molecular scale due to the creation of the void (and thus an increase in entropy), the mechanical action itself, and an increase in the local temperature from the pullout process. Note that the degree of structural reordering that is expected is also dependent on the degree of entanglements, and thus less entangled systems such as PI are expected to have more structural rearrangements than amorphous polymers like polypropylene.
In Figure 7, E inter (Eq. 2) from the pullout simulation between the PI polymers and each CNT is given as a function of time and the applied constant external force. For the SWCNT system in Figure 7(a), E inter remains relatively constant for low applied external forces (less than 20.82fN), indicating that the polymer chain displacement is on the same scale as the SWCNT displacement; thus, both Phases I and II are occurring in these systems. For the system with an applied force of 20.82fN, there is a marked decrease in E inter to 0.17 kJ/ kg as time progresses but then levels off. This trend is due to the stable interfacial binding interaction and the decrease of the contact area during the pullout simulation; in other words, Phase I dominates since the SWCNT displacement at this pullout force is too fast for the movement of the PI polymers to keep up. Consequently, 20.82fN can be regarded as a limiting force in this pullout system, below which the force is not   Table 1). Similar observations can also be made for the DWCNT in Figure 7(b) and the TWCNT in Figure 7(c). However, the limiting forces differ. For both the SWCNT and DWCNT, a force of 20.82 fN causes a relatively constant E inter after about 0.30 ns, while for the TWCNT, a larger force of 24.29 fN is required to reach the plateau. This difference could be attributed to more energy being needed to overcome the interactions among the walls, which is governed by the intertube van der Waals forces.

Variation of other properties with applied external force
Nanocomposite properties as a function of the external force are given in Figure 8; note, that just the external forces beyond the limiting force from Figure 7 are included. The change in displacement (dL) is defined as the difference between the initial CNT position and the CNT position at complete pullout. For the change in the displacements in Figure 8(a) and Table 1, the largest values are observed at lower forces due to the prevalence of coupled polymer motion (Phase II). The effect is comparable for the SWCNT and DWCNT, and even more prominent for the TWCNT. As Phases I and III begin to dominate, the change in the CNT displacement decreases significantly. Figure 8(a) also reflects that the displacements of the three CNT types are the same values when large forces are applied. In this work, constant external forces were applied on the CNT, where the work (W) done by the external force (F) can be described as the force-distance product, which can be calculated for each step. Thus, the pullout energy (E pullout ) can be defined as, [54] where E pot is the potential energy of the system. The state before pullout is labeled "1", and the state after pullout is labeled "2". The E pullout is plotted as a function of external force in Figure 8b and is also summarized in Table 1.
These calculations indicate that the difference between the SWCNT and DWCNT systems for E pullout is minor, but is slightly higher for the TWCNT, which is consistent with the trends for the interaction energies in Figure 7. It is also worth noting that E pullout changes dramatically between 27.76 fN and 34.70 fN for all three types of CNTs. By analyzing the values in Table 1, it appears that both the displacement and the potential energy of the system is significantly greater at this applied external force. The interfacial shear strength (IFSS) between the CNT and the matrix can strongly influence the mechanical properties of the CNT composite, and thus it is an important parameter to calculate. The IFSS (τ i ) can be expressed in terms of E pullout according to the following relation, where x is the coordinate along the longitudinal tube axis [51], and r and L are the outer radius and displacement, respectively, of the CNT during the pullout. Solving for IFSS gives the following equation: Figure 8(c) presents the trend of IFSS with the change of applied external forces. At low applied external forces, the IFSS is calculated to be around 2 to 3 MPa for all three CNT types. These observations indicate that at small external forces, the induced polymer movement plays a dominant negative contribution during pullout, which confines the IFSS to less than 3 MPa. However, at higher applied external forces, the IFSS for the TWCNT reaches values of over 18 MPa, which is significantly greater than for SWCNT (4.2 MPa) and DWCNT (7.4 MPa). According to Eq. 5, the IFSS is comprised of the ratio of E pullout to dL. Since dL is equivalent for all three CNT types at high applied external pressures (Figure 8a), this difference must arise due to E pullout . According to Eq. 4, E pullout is comprised of the work done by the external force, which has the same value for all three CNTs since dL is constant at higher applied external forces, as well as the ∆E pot and ∆E inter of the system. The high potential energy changes may be attributed to the inter-wall interactions and also active polymer configuration changes during the pullout simulations (Phase II and III in Figure 6).
Since the IFSS value of CNT/BPDA-PDA PI composites has not been tested experimentally yet, we have no accurate data to which to compare. However, it is still meaningful that our simulation value is in the reasonable range for other CNT/polymer composites [3][4][5]. Therefore, this work provides not only a prediction of the IFSS as a function of the CNT walls, but also molecular-level details of the interfacial interaction between BPDA/PDA and CNTs that impact those IFSS values.

Conclusions
MD simulations were conducted to investigate the interfacial properties between CNTs and a BPDA-PDA PI matrix. The E inter between CNTs and a single chain of BPDA-PDA was as high as 0.400 ± 0.004 kJ/kg, and significant ordering of the PI chains was observed at the CNT interface, especially as the number of walls increased, suggesting that intertube van der Waals interactions can impact the matrix structure, depending on how it is processed.
The corresponding E inter and E pullout values for the nanocomposite reflect that there is a limiting value for external force, below which the pullout process cannot complete. Also at low external forces, the polymer chains were observed to follow the movement of the CNT, which is very important aspect of the analysis of the interfacial shear strength. The wall number effect is distinct in the pullout energy and interfacial shear strength, which is attributed to the intertube van der Waals forces. The interfacial shear strength was calculated to be 18.7 MPa as the peak value, which is reasonable in CNT-reinforced composites. This simulation approach is applicable for any polymer/CNT reinforcement nanocomposite and thus can be utilized for optimizing the desired properties of the nanocomposite by considering the interfacial characteristics.