Mechanical Response of Carbon Nanotube Bundle to Lateral Compression

: Structure evolution and mechanical response of the carbon nanotube (CNT) bundle under lateral biaxial compression is investigated in plane strain conditions using the chain model. In this model, tensile and bending rigidity of CTN walls, and the van der Waals interactions between them are taken into account. Initially the bundle in cross section is a triangular lattice of circular zigzag CNTs. Under increasing strain control compression, several structure transformations are observed. Firstly, the second-order phase transition leads to the crystalline structure with doubled translational cell. Then the ﬁrst-order phase transition takes place with the appearance of collapsed CNTs. Further compression results in increase of the fraction of collapsed CNTs at nearly constant compressive stress and eventually all CNTs collapse. It is found that the potential energy of the CNT bundle during deformation changes mainly due to bending of CNT walls, while the contribution from the walls tension-compression and from the van der Waals energies is considerably smaller.


Introduction
Carbon atoms can create a huge number of allotropes because four valence electrons can form the sp 1 , sp 2 and sp 3 bond configurations. The action of van der Waals forces produces a variety of secondary structures including structures with a translational symmetry such as fullerite crystal [1][2][3], graphite [4,5], and crystals made of carbon nanotubes (CNTs) [6][7][8]. Such van der Waals crystals have properties not exhibited by their structural elements, for example, chemical properties typical for molecular crystals and peculiar mechanical properties related to deformability of the structural units [1][2][3][4][5][6][7][8]. This work is devoted to the mechanical properties of CNT bundles under lateral compression.

Simulation Setup
Lateral biaxial compression of CNT bundle is considered in this study under plane strain condition. The CNTs are aligned along the z-axis creating a triangular lattice in cross-section. For definiteness, the zigzag CNTs are considered. In Figure 1, the computational cell is schematically shown and the geometrical parameters are specified: a is the distance between carbon atoms in the CNT wall, d is the distance between the walls of neighboring CNTs, D is CNT diameter and A = D + d is the distance between centers of neighboring CNTs. The atoms move on the (x, y) plane and each atom stands for a rigid row of carbon atoms oriented along the z-axis. The computational cell includes 20 × 24 CNTs (the case of 2 × 2 CNTs is presented in Figure 1). Each CNT has 30 atoms and total number of atoms in the computational cell is N = 20 × 24 × 30 = 14400. Much smaller computational cell can be used for simulation of structures with a long-range order, but for irregular structures a representative volume should be considered and, as it will be shown later, the cell with 20 × 24 CNTs is sufficiently large. . Each zigzag CNT is represented by 30 carbon atoms that move on the (x, y) plane. Each atom stands for a row of atoms oriented normal to the (x, y) plane, which moves as a rigid body. Geometry parameters are as follows: a is the distance between carbon atoms in the CNT wall, d is the distance between the walls of neighboring CNTs, D is CNT diameter and A = D + d is the distance between centers of neighboring CNTs. The computational cell in the form of a parallelogram has the sides I × A and J × A. The computational cell is subject to the periodic boundary conditions. Indices p, q, s, f , and f are used to describe the interatomic interactions.
Each atom has two degrees of freedom, the components of the displacement vector in the (x, y) plane, with the total number of degrees of freedom equal to 2N = 28800. Periodic boundary conditions are used in both directions.
Numerical values of the geometry parameters are as follows. The equilibrium interatomic distance in graphene is ρ = 1.418 Å. The equilibrium distance between rigid atomic rows oriented along the z-axis is a = ρ √ 3/2 = 1.228 Å, CNT diameter is D = a/ sin(π/30) = 11.75 Å, the equilibrium distance between CNT walls under zero pressure is d = 3.088 Å, and the distance between centers of neighboring CNTs is A = D + d = 14.838 Å. Note that CNT diameter is relatively small in our study and separate CNT in the absence of external forces can have only circular stable configuration.
In our simulations the units of time, energy and distance are picoseconds, eV, and angstrom, respectively. In these units the mass of carbon atom is M = 12 × 1.0364 × 10 −4 .
The CNT bundle is described by the Hamiltonian (total energy) [69,70] where the four terms in the right-hand side give the kinetic energy, the energy of valence bonds, valence angles, and van der Waals non-valence interactions, respectively. These energies are fully described in [70]. Here we reproduce the expressions for the energies for the sake of the reader. In Figure 1, indices p, q, and s number nearest atoms in the wall of a CNT and f is the atom belonging to a different CNT, while f belongs to the same CNT as the atom q and the distance between q and f is not less than 3a. Then the energies in the right-hand side of Equation (1) are specified as follows.
Kinetic energy is where r q is the radius-vector of q-th atom, overdot means differentiation with respect to time, and summation is over all atoms in the computational cell. We note that the kinetic energy term is not taken into account in this work, since only the relaxational dynamics of atoms is considered, but we present it for completeness of the description of the model. The energy of valence bonds is defined by the harmonic potential which is summation over the bonds connecting nearest atoms (see Figure 1) and the number of such bonds is equal to the number of atoms. The longitudinal stiffness of graphene is well reproduced when k = 405 N/m [69]. The energy of valence angles is given by the anharmonic potential which is summation over the valence angles formed by the atoms p, q, and s (see Figure 1) and the number of such angles is equal to the number of atoms. Bending stiffness of graphene is well reproduced when = 3.50 eV [69]. The energy of van der Waals non-valence interactions is given by the (5,11) Lennard-Jones potential with the parameters ε = 0.00166 eV and σ = 3.61 Å [69]. Here summation is over all interatomic bonds connecting atoms q and f within the cut-off radius of 6 Å, and for atoms belonging to the same CNT (denoted as f in Figure 1), as it was already mentioned, the bonds shorter than 3a are not taken into account. Note that the summands in Equation (5) give the interaction potentials between pairs of rigid atomic rows oriented along the z-axis. This potential was obtained by summation of the C-C interatomic energies described by the (6,12) Lennard-Jones potential. The net interaction energy between two atomic rows is best fitted by the (5,11) Lennard-Jones potential [69].
In the course of fitting the model parameters, strong sensitivity of the simulation results to small changes in the parameters was not observed.
In this study, mechanical properties of CNT bundle under lateral biaxial compression in plane strain conditions with ε xx = ε yy ≤ 0 and ε xy = 0 are evaluated. The strain state of the structure is characterized by the absolute value of the volumetric strain |θ| = |ε xx + ε yy |. Temperature effect is not taken into consideration and the perturbation-relaxation molecular dynamics simulations are carried out. Equilibrium structures are found after each increment of compressive volumetric strain |∆θ| = 0.02 at T = 0 K. The simulation algorithm is as follows. Incremental compressive strain is applied with the step ∆ε xx = ∆ε yy = −0.0025. Each strain increment is followed by the perturbation of coordinates of the atoms by adding random displacements uniformly distributed over the range from −10 −6 to 10 −6 Å. After that, the potential energy of the system is minimized using the gradient method and the equilibrium structure is obtained. It is assumed that the energy minimization is complete when the absolute value of the maximal atomic force does not exceed 10 −10 eV/Å.

Results and Their Discussion
Here we report the simulation results starting from the structure analysis and then moving to the analysis of the energy-strain and stress-strain relations.

Structure Evolution
During step-wise biaxial compression the structure of CNT bundle exhibits several qualitative transformations (the results for the small volumetric strain can be found in [70]). In Figure 2a  For simulation of the crystalline structures shown in Figure 2a,b a much smaller computational cell can be used and even DFT simulations are possible for these cases. For the irregular structures presented in Figure 2c,d, a computational cell with a few CNTs is insufficient for evaluation of the averaged mechanical properties.
In order to quantify the structural changes, we define the ellipticity of CNT cross section where a and b are the largest and smallest dimensions of the cross section. We also define the orientation angle α of the CNT cross section. These parameters are shown in the inset of Figure 3a, left panel.
Parameters η and α are calculated for each CNT and then the 20-bin histograms are constructed. In Figure 3a Figure 2b). In (c) the histogram for η has a sharp peak at η ≈ 0.2 and a broad hump centered at η ≈ 0.5. The sharp peak stands for the collapsed CNTs and the hump represents the non-collapsed CNTs with various ellipticity. The orientation angle histogram in (c) does not reveal preferable orientation angles (see Figure 2c). In (d), the broad hump disappears in the distribution of η since all CNTs collapse (see Figure 2d). On the other hand, the histogram of orientation angles in (d) features three distinct peaks at 0, 60 and 120 degrees.

Energy-Strain and Stress-Strain Curves
In Figure 4, as the functions of the volumetric compressive strain, we plot the potential energy of the relaxed CNT bundle, U, and its three components: the energy of valence bonds U B , valence angles U A , and van der Waals interactions U vdW , see Equation (1).
For compressive volumetric strain less then 0.09, before the collapsed CNTs appear in the system, the potential energy U increases with strain quadratically. When structure compression takes place due to the increase in the fraction of compressed CNTs, U grows with strain linearly. For |θ| > 0.37, when nearly all CNTs are collapsed, potential energy starts to increase with strain faster than linearly.
It can be seen from Figure 4 that dominant contribution to the potential energy U comes from the energy of the valence angles U A and the other two components are much smaller. In the range of volumetric strain 0.09 < |θ| < 0.37 a reduction of the energy of van der Waals interactions U vdW is observed. This is explained by the formation of new van der Waals bonds when CNTs collapse and their walls start to interact. Smallness of the energy of the valence bonds U B means that CNT walls in the considered type of loading can be assumed as inextensible. It is also instructive to analyze the stress-strain curves which are presented in Figure 5. This curves reflect the two phase transitions in the CNT bundle subject to biaxial compression. The second-order phase transition takes place at |θ| ≈ 0.07, when CNTs obtain elliptic form, the translational cell doubles, and the long-range order is still preserved. Stress components do not show discontinuity at this transition point but the slope of the stress-strain curves changes. The first-order phase transition occurs at |θ| ≈ 0.12 when sudden drop of the normal stress components is observed due to appearance of the collapsed CNTs. In the regime of compression with increasing fraction of collapsed CNTs normal stress components are almost constant and they start to grow when collapsed CNTs percolate.
We also note that during the whole studied range of compressive strain the equality of normal stress components nearly holds, σ xx ≈ σ yy , and shear stress remains nearly zero, σ xy ≈ 0, as it should be for biaxial compression. This means that the considered computational cell size is sufficiently large.

Conclusions
Relaxational dynamics simulations were performed in the frame of the chain model developed in Refs. [69,70] to analyze the mechanical response of CNT bundle to biaxial compression. It was found that the CNT bundle undergoes two phase transitions. At the value of the compressive volumetric strain of |θ| ≈ 0.07 the second-order phase transition was observed. As a result of this transition period doubling took place while crystallinity of the structure was preserved. Cross section of CNTs at the transition point changed from nearly circular (Figure 2a) to elliptic (Figure 2b). At |θ| ≈ 0.12 the first-order phase transition took place when collapsed CNTs appeared in the system (Figure 2c). Further compression was observed at nearly constant values of normal stress due to the increase in the portion of collapsed CNTs. At |θ| ≈ 0.4 all CNTs were found collapsed (Figure 2d).
It was found that the deformation mechanism of CNT bundle was bending of CNT walls with very little tensile deformation of the walls and the main contribution to the potential energy of the system comes from the energy of valence angles, U A , (see Figure 4). Van der Waals energy, U vdW , slightly decreased in the regime of growing portion of collapsed CNTs due to formation of new bonds.
The chain model used in the present study assumes that the atomic rows along the armchair direction are rigid. All structures analyzed here correspond to this assumption and, therefore, can be obtained in case of refusal of plane strain conditions. In the case of full atomic modeling CNTs can collapse in another fashion, by initiation of the local collapsed region and its propagation along the tube. If this scenario is energetically favorable, then the first-order phase transition should happen somewhat earlier than in the case of plane strain conditions. Another limitation of the chain model is related to the simulations at finite temperature. Phonon waves propagating along the CNTs are completely suppressed by the plane strain conditions. All these effects should be addressed in frame of a full atomic model. In the present study full atomic simulations were not conducted, but in the previous works with smaller number of degrees of freedom it was shown that in many cases the chain model gives physically meaningful results close to those obtained in full atomic simulations [69][70][71][72][73][74][75].
Modification of the chain model to the case of the armchair CNT bundle is straightforward but the model parameters will be somewhat different from that for the zigzag CNTs.
Overall, our results contribute to understanding of the mechanisms of deformation and mechanical response of CNT bundle to lateral compression.
Future research directions include consideration of bundles composed of CNTs of different diameter or multi-walled CNTs. It is possible to modify the chain model used here for modelling other graphene-analogous 2D nanomaterials [76]. Nonlinear dynamics and spatially localized excitations such as discrete breathers and rotobreathers [77] will also be analyzed in CNT bundles. Funding: This research was funded by the Russian Foundation for Basic Research, grant number 18-29-19135.