Chain Model for Carbon Nanotube Bundle under Plane Strain Conditions

Carbon nanotubes (CNTs) have record high tensile strength and Young’s modulus, which makes them ideal for making super strong yarns, ropes, fillers for composites, solid lubricants, etc. The mechanical properties of CNT bundles have been addressed in a number of experimental and theoretical studies. The development of efficient computational methods for solving this problem is an important step in the design of new CNT-based materials. In the present study, an atomistic chain model is proposed to analyze the mechanical response of CNT bundles under plane strain conditions. The model takes into account the tensile and bending rigidity of the CNT wall, as well as the van der Waals interactions between walls. Due to the discrete character of the model, it is able to describe large curvature of the CNT wall and the fracture of the walls at very high pressures, where both of these problems are difficult to address in frame of continuum mechanics models. As an example, equilibrium structures of CNT crystal under biaxial, strain controlled loading are obtained and their thermal stability is analyzed. The obtained results agree well with previously reported data. In addition, a new equilibrium structure with four SNTs in a translational cell is reported. The model offered here can be applied with great efficiency to the analysis of the mechanical properties of CNT bundles composed of single-walled or multi-walled CNTs under plane strain conditions due to considerable reduction in the number of degrees of freedom.


Introduction
There exist a huge number of carbon polymorphs, including a wide class of sp 2 structures such as fullerenes, carbon nanotubes (CNT), and graphene. Due to the action of relatively weak van der Waals forces, a great variety of secondary structures can be formed, and some of them can have a long-range order, for example, fullerite crystal composed of fullerenes [1][2][3], graphite made of graphene layers [4,5], and crystals made of CNTs [6][7][8]. Such crystalline structures are of great interest of simulation methods. Continuum mechanics is a powerful and effective tool to successfully describe macroscopic parameters of CNT bundles, but it has some limitations. For example, thermal fluctuations of CNTs and their fracture can be more adequately modelled in frame of atomistic models. On the other hand, full atomic models, as mentioned above, are very demanding on computational resources. One possible compromise is to use atomistic models for particular deformation modes, when the number of degrees of freedom can be substantially reduced.
In this study, in order to reduce the number of degrees of freedom, a full atomic model of CNT bundles under plane strain conditions is substituted by the chain model developed in the work [60] and successfully used to study structure and properties of secondary structures such as folds and scrolls of carbon nanoribbons [60][61][62][63][64] and dynamics of surface ripplocations on a graphite substrate [65]. For simplicity, here we will only consider the case of a bundle composed of single-walled CNTs of equal diameter, but the model can be applied to the cases of CNTs of different diameter, multi-walled CNTs, and include graphene scrolls and cylindrically crumpled graphene.

Materials and Methods
The computational model employed in this study is schematically shown in Figure 1. The nanotube bundle is aligned along the z-axis and CNTs of equal diameter create in cross-section a triangular lattice; they are numbered by the indices i = 1, ..., I and j = 1, ..., J (the case of I = J = 2 is shown). Only zigzag CNTs are considered for simplicity. The carbon atoms move on the (x, y) plane and each atom represents a rigid row of atoms oriented normal to the (x, y) plane. Within each CNT, carbon atoms are numbered by the index n = 1, ..., N, anti-clockwise, starting from the atom with maximal x-coordinate. Thus, total number of atoms in the computational cell is I × J × N. Atomic positions are defined by the radius-vectors r ijn = (x ijn , y ijn ). Periodic boundary conditions are imposed. Schematic of the computational cell that includes I × J carbon nanotubes (CNTs) (I = J = 2 in this case) numbered by the indices i = 1, ..., I and j = 1, ..., J. Carbon nanotubes in cross-section create a triangular lattice. Within each CNT, carbon atoms are numbered by the index n = 1, ..., N anti-clockwise, starting from the atom with maximal x-coordinate. Atoms move on the (x, y) plane. Each atom represents a row of atoms oriented normal to the (x, y) plane, which moves as a rigid body. The computational cell has the shape of a parallelogram with the sides I × A and J × A, where A is the distance between centers of neighboring CNTs. Periodic boundary conditions are used.
Let us describe the model geometry. The interatomic distance in graphene is equal to ρ = 1.418 Å. The distance between neighboring atomic rows oriented along the armchair direction in graphene is then a = ρ √ 3/2 = 1.228 Å, and this is the distance between atoms in the chain model (see Figure 1). Carbon nanotube diameter is D = a/ sin(π/N). Let d be the shortest distance between CNT walls, then the distance between centers of neighboring CNTs is A = D + d. The sides of the computational cell in the form of parallelogram are I × A and J × A. In our simulations, we consider CNTs with N = 30 having diameter D = 11.75 Å and equilibrium value of d = 3.088 Å, which can be compared to the interplanar distance of graphite equal to 3.3 Å.
Carbon nanotube bundle under uniform lateral compression can be efficiently described by the Hamiltonian of the chain model [60] which includes kinetic energy energy of valence bonds energy of valence angles and energy of van der Waals interactions In Equation (2) M is the carbon atom mass, which is 12 amu. In our simulations, time is measured in picoseconds, energy in eV, and distance in angstrom. In these units M = 12 × 1.0364 × 10 −4 . As it can be seen from Equation (3), harmonic potential with stiffness k is used to model deformation of the valence bonds. In order to reproduce the longitudinal stiffness of graphene sheet one should take k = 405 N/m [60], which in the units adopted here gives k = 25.279.
In Equation (4), cosine of the angle between two valence bonds, r ijn − r ijn−1 and r ijn+1 − r ijn , is calculated as Bending rigidity of graphene sheet is well reproduced with the value of the potential parameter = 3.50 eV [60]. The van der Waals interactions in Equation (5) are given by the Lennard-Jones potential (5,11) [65] with the interaction energy ε = 0.00166 eV and the equilibrium bond length σ = 3.61 Å. Further information on the chain model and on the procedure of fitting its parameters can be found in [60,65].
As it has been mentioned, a CNT of sufficiently large diameter can have either cylindrical or collapsed equilibrium configuration. In the present study we consider CNTs of relatively small diameter (N = 30, D = 11.75 Å) with only circular stable state when unloaded.
The aim of this study is to evaluate mechanical response of CNT bundle to lateral biaxial compression under plane strain condition with ε xx = ε yy ≤ 0 and ε xy = 0. Firstly, equilibrium configurations are found at zero temperature and then, their stability at room temperature (T = 300 K) is analyzed.
Perturbation-relaxation molecular dynamics simulations are done at zero temperature in order to find equilibrium structures at different values of applied strain and T = 0 K. The simulation protocol is as follows. The compressive strain is applied by increments ∆ε xx = ∆ε yy = −0.0025 starting from zero strain. After each increment, the positions of atoms are perturbed by adding small random displacements to their xand y-coordinates. The displacements are uniformly distributed in the range from −10 −6 to 10 −6 Å. Then the equilibrium structure is obtained by minimizing potential energy of the system with the help of the gradient method. Energy minimization stops when the absolute value of the maximal force acting on atoms becomes smaller than 10 −10 eV/Å. Different computational cell sizes were considered. Calculations with 24 × 24 = 576 CNTs have revealed that structures with period doubling are formed as a result of instability at particular value of compressive strain. Such structures can be analyzed with smaller cell size and most of the results reported here are for the cell that includes 6 × 6 = 36 CNTs.
Classical molecular dynamics was used to assess stability of equilibrium structures with respect to thermal fluctuations at T = 300 K. Temperature in our simulations is defined as where k B = 8.617 × 10 −5 eV·K −1 is the Boltzmann constant and the averaging time is t 2 − t 1 =10 ps. For a given temperature T, the initial velocities of atoms are assigned according to the Maxwellian distribution. Random initial displacements of atoms are assigned in a way to increase the potential energy of the system by the amount equal to the kinetic energy.
Equations of atomic motion that stem from the Hamiltonian Equation (1) are integrated numerically with the help of the Stormer method of order six with the time step of 0.1 fs. The structure is considered to be stable if no structure transformations are observed within 100 ps. Structure transformations can be very well seen on the time dependencies of kinetic and potential energies during our simulations with NVE ensemble (constant number of particles, volume, and total energy). When structure changes, kinetic energy increases in expense of potential energy.

Results
In this section the equilibrium structures of CNT bundle under lateral compression are reported and their properties are analyzed. First, the potential energy and stress as functions of strain are given and then the change of CNT geometry with strain is presented. Finally, stability of equilibrium structures at T = 300 K is analyzed.

Energy and Stress in the System
We start with the analysis of potential energy per atom calculated for equilibrium structures at different values of compressive strain. In Figure 2a, total potential energy per atom is shown as a function of strain, while in Figure 2b-d this energy is decomposed into three parts: the energy of van der Waals interactions, the energy of valence bonds, and the energy of valence angles, respectively. Total potential energy in the range of strain below 3.75% increases with strain quadratically but for larger strain a linear increase of energy with strain can be observed. These two regimes are separated in Figure 2 by the vertical dashed line. This qualitative change in the behavior of potential energy is due to structural changes observed in the system for strain exceeding the critical value of 3.75%. Below this critical strain, all CNTs in the system have the same cross-section in the form of six-fold flattened cylinders (see Figure 3a), while above the critical value the CNTs become elliptic. Two different structures with elliptic CNTs can form, Structure I with the translational cell doubled in one direction (see Figure 3b) and Structure II with the translational cell doubled in two directions (see Figure 4). In Figure 2   From Figure 2a-d it is clear that below the critical strain the potential energy of all three kinds increase with strain. At the critical value of strain just before the structure transformation one has U VdW = 0.0041 eV, U B = 0.0056 eV, and U A = 0.0029 eV. The largest increase of energy is observed for the valence bonds because main mechanism of lattice deformation in this regime is contraction of valence bonds. The smallest contribution to the energy increase comes from the valence angles, since they do not change much during transformation of CNT cylindrical shape into six-fold flattened shape. This picture drastically changes for compressive strain larger than 3.75% when CNTs become elliptic. The deformation of structure in this regime is mainly due to change of valence angles and U A increases rapidly with strain while other two components of energy decrease with strain. The decrease of U VdW with strain in this regime is explained by formation of new van der Waals bonds with increasing ellipticity of CNTs. At strain of 8.75%, U A is already one order of magnitude larger than two other components of energy. Note that both Structure I and Structure II have very close energies, and only for strain above 7% is U VdW in Structure I slightly higher than in Structure II. Repetition of the simulations with different random atomic displacements has revealed that Structure I and Structure II are formed at the transition point with nearly equal probability, and this is because they have practically same energy.
Variation of stress components σ xx , σ yy , and σ xy with strain is shown in Figure 2e-g, respectively. For strain below the critical level, compressive stress increases so that σ xx = σ yy and σ xy = 0. This is because the structure with six-fold symmetry is elastically isotropic (see Figure 3a). For strain above the critical level, deformation occurs at nearly constant pressure p = −(σ xx + σ yy )/2. Both structures become anisotropic. In particular, Structure I is orthotropic with |σ xx | > |σ yy | and σ xy = 0. On the other hand, Structure II has general anisotropy with |σ xx | < |σ yy | and σ xy = 0.

Geometry of CNTs
To better understand the relation between structure and macroscopic parameters of the system, let us quantify the geometry of CNTs in different structures. For each CNT we calculate its minimal and maximal diameters, D min and D max , and the angle α between the x-axis and the maximal diameter, as shown in Figure 3b.
In Figure 5a,b the ratio D min /D max is shown for Structures I and II, respectively. In Figure 5c,d the angle α is shown for Structures I and II, respectively.
Since the translation cell of Structure I includes two CNTs (numbered as 1 and 2 in Figure 3b), in Figure 5a,c, two different values of D min /D max and α can be seen for strain above the critical value.
The translation cell of Structure II includes four CNTs (numbered as 1 to 4 in Figure 4) and hence, in Figure 5b,d, four different values of D min /D max and α can be seen for strain above the critical value. For strain below the critical value, all curves merge into one since translation cell includes single CNT, see Figure 3a. From Figure 5a, one can see that in Structure I the ellipticity of SNTs gradually increases with increasing compressive strain. At the value of strain 8.75%, D min /D max is smaller than 0.5. Recall that further increase of strain above 9% results in the formation of collapsed CNTs with non-convex cross-section, but we do not analyze such structures here. As Figure 5c suggests, the orientation of elliptic CNTs in Structure I is practically strain-independent. Carbon nanotubes 1 and 2 have orientation angles α = 35 and 145 deg., with the difference equal to 110 deg.
Structure II demonstrates more complicated evolution with strain. Two different regimes can be distinguished looking at Figure 5b. For strain between 4% and 6.5%, three CNTs (1, 2, and 4) in the translation cell have nearly the same ellipticity, while CNT 3 is less elliptic (see Figure 4a). For larger strain, CNTs 2 and 4 become more elliptic than CNTs 1 and 3 (see Figure 4b). Orientation angles of CNTs also change with strain (see Figure 5d). At 8.75% strain, CNT 4 is nearly aligned with the y-axis (α is close to 90 deg.), while CNT 2 with the x-axis (α is close to 180 deg.). Carbon nanotubes 1 and 3 have angles of about 50 and 130 deg., with a difference of about 80 deg.

Temperature Effect
The stability of all three types of equilibrium structures reported in Section 3.1 with respect to thermal oscillations is investigated here at temperature T = 300 K.
It was found that the structure with all identical CNTs (see Figure 3a) is stable at room temperature up to compressive strain of about 4.0%, which slightly exceeds the stability range of this structure at 0 K (3.75% strain). At a strain of 4.0% and temperature 300 K, the cross-sectional shapes of CNTs fluctuate in time but on average all CNTs remain the same, preserving the high symmetry of the structure. The pressure-induced phase transition from high-symmetry structure to the structures with period doubling in one or two directions is the second-order phase transition. This is justified by the absence of the jumps in macroscopic properties at the transition point (3.75% compressive strain), see Figure 2. Low-symmetry Structures I and II under heating transform into high-symmetry structures, if the compressive strain is not too high.
Within the range of compressive strain from 4.1% to 6%, Structures I and II are stable at 300 K; they are preserved within the simulation time of 100 ps and no jumps of macroscopic parameters are observed.
Both Structures I and II become unstable at 300 K for compressive strain exceeding 6%. The instability of Structure I is illustrated in Figure 6 for compressive strain of 7%. In (a,b), one can see the time evolution of kinetic energy per atom and components of compressive stress, respectively. Until t = 10 ps, kinetic energy oscillates near the value of 0.0259 eV, which corresponds to 300 K, but then it starts to increase in expense of the potential energy (total energy is conserved in the system). Pressure drops at the transition point from 175 to 140 MPa. Jumps in the macroscopic parameters of the system indicate that this phase transition is of the first order. In (c), a snapshot of the structure is presented at t = 20 ps. One can see that the long-range crystal order is lost and an irregular structure that includes collapsed CNTs with non-convex cross-section is formed.

Discussion
The chain model introduced in the works [60,65] was developed here to enable simulation of the mechanical properties of CNT bundles under plane strain conditions. The model was applied to the analysis of structure transformations and mechanical properties of CNT crystal subjected to biaxial lateral compression. Carbon nanotube diameter is relatively small, so that the collapsed shape is unstable in the absence of external forces.
Three different crystalline structures stable at zero temperature have been found. For compressive strain |ε xx | = |ε yy | < 3.75%, primitive translational cell of the crystal includes single CNT. In the range of compressive strain from 4% to 9%, two phases of nearly the same potential energy were found, one has two and the other one has four CNTs in a primitive translational cell. These structures are referred to as Structure I and II, respectively. Pressure-induced phase transition from the high-symmetry structure to the Structures I or II is of the second order because no jumps of macroscopic properties are seen in Figure 2 at the transition point at strain of 3.75%.
Thermal fluctuations increase the stability range of the high-symmetry structure with single CNT in the primitive translational cell. At 300 K, this structure is stable up to 4.0% compressive strain, while at 0 K the stability threshold is at 3.75% strain. The transformation of low-symmetry phase into high-symmetry phase under an increase in temperature is typical for the second-order phase transitions [66,67].
Thermal fluctuations reduce the stability range of Structures I and II. At zero temperature, they are stable from 3.75% to 8.75% of compressive strain, while at room temperature the stability range of compressive strain is from 4.1% to 6%. The transition above 6% strain is of the first order, it is accompanied by the jumps in macroscopic parameters (see Figure 6a,b), when crystalline structure with a long-range order transforms into an irregular structure (see Figure 6c).
As for the mechanical properties of CNT bundles under lateral compression, the transition to the structures with elliptic CNTs results in a considerable drop in rigidity of the bundle. Indeed, the deformation of the structure in the range from 3.75% to 9% compressive strain is at nearly constant pressure, see Figure 2e,f. Unloading of the system from any strain below 9% has shown that the structure is non-linear elastic with no hysteresis effect. Also note that the high-symmetry structure is elastically isotropic, while Structures I and II are anisotropic since for them σ xx = σ yy , see Figure 2e,f.
The results reported here are in agreement with the results of full atomic and continuum mechanics modelling [6,8,39,51,53,54,[56][57][58], but they were obtained at a very low computational cost. For example, full atomic modelling with the same accuracy would require the use of the periodic boundary conditions in the z direction with at least one zigzag carbon chain within the translational cell. The calculation of interatomic forces between atoms within the considered cell and its translation images would require additional summation, which is absent in the chain model since it has been done in derivation of the effective potentials between the rigid atomic chains oriented along the z axis. For this reason, the chain model gives at least one order of magnitude acceleration of computations with the same accuracy, as compared to full atomic modelling.
Recall that the harmonic, unbreakable potential is used in this work to describe the valence interactions between carbon atoms, which is sufficient for modelling structure transformations at relatively small pressure considered here. In order to model irreversible deformation of CNT bundles under very high pressure, see, e.g., the work [39], the harmonic potential Equation (3) should be substituted with the breakable anharmonic potential, such as Morse potential [68].
In future works, it is important to study the effect of CNT diameter since new effects can be expected for larger diameter when collapsed isolated CNT is stable. Crystals composed of such CNTs can demonstrate irreversible plastic deformation with very peculiar mechanisms of plasticity. As for the applications, CNT bundles under lateral compression can show hysteresis effect, when it acts as an elastic damper [69], in a similar way to the compressed, vertically aligned CNT brushes [31].
The chain model proposed here can be readily adjusted to a number of newly found graphene-analogous 2D nanomaterials [70] by fitting model parameters to the results of first-principle or molecular dynamics simulations.
Consideration of bundles composed of different CNTs or multi-walled CNTs is straightforward. Substitution of the harmonic valence bond potential in Equation (3) with a suitable breakable anharmonic potential will enable the simulation of structure transformations in the CNT bundle under very high pressure.

Conclusions
We thus conclude that the chain model can be applied with a high numerical efficiency and sufficient accuracy to the analysis of structural and mechanical properties of CNT bundles under plane strain conditions. The atomistic chain model proposed here, unlike continuum mechanics models, is able to describe high curvature of collapsed CNT wall and fracture of the walls under high pressure.
The chain model proposed here can be readily applied to the cases of CNTs of different diameter, collapsed CNTs, multi-walled CNTs, and even include graphene scrolls and cylindrically crumpled graphene.