## Abstract

A new method is developed for the simulation of carbon nanotubes (CNTs) based on the molecular mechanics (MM) and cellular automata (CA) algorithms. The proposed atomic-based CA algorithm (ACAA) is as accurate as the MM method, but much faster for the simulation of CNTs. In the ACAA model, a CNT is treated as a discrete system in which the interaction between atoms is described by using Tersoff–Brenner's many-body potential, but the solution framework is based on the conventional continuum mechanics. The buckling and post-buckling behaviours of defective CNTs are studied by using the ACAA. The numerical results show that the buckling occurs at the vacancy defect site. This phenomenon indicates that vacancy defects lead to a geometrical imperfection of CNTs. The study demonstrates that vacancy defects can significantly reduce the buckling load of CNTs.

## 1. Introduction

Much research has focused on the mechanical properties of carbon nanotubes (CNTs) since their discovery by Iijima in 1991. However, most of these studies have considered defect-free CNTs. For example, continuum mechanics models have been developed for the analysis of CNTs (Qian *et al*. 2002), a finite deformation continuum theory has been established for the mechanics of CNTs based on the exponential Cauchy–Born rule (Arroyo & Belytschko 2004), and a molecular mechanics (MM) model has been proposed for the critical buckling strain of multi-walled CNTs under axial compression and bending (Chang *et al*. 2005). Relatively little is known about the mechanical behaviour of CNTs with defects, even though CNTs can contain a range of structural defects, such as Stone–Wales defects, bond-breaking defects and atomic vacancies, all playing a significant role in determining their mechanical properties.

The effects of vacancy defects have been studied using atomistic methods. For example, molecular dynamics (MD) simulation has shown that a vacancy has a weak effect on the Young's modulus of a CNT, but a strong effect on its tensile strength (Sammalkorpi *et al*. 2004), and that local buckling occurs easily at the locations of vacancy defects (Xin *et al*. 2007). The MM simulation has been used to investigate the mechanical behaviour of CNTs (Hou & Xiao 2007) and has indicated that vacancy defects decrease the strength of CNTs, but improve their bending stiffness. Many researchers have studied the role played by vacancy defects in CNT fractures by employing both MD and MM simulations (Yakobson 1998; Banhart 1999; Wagner 2002; Xiao & Hou 2006). One study (Mielke *et al*. 2004) has considered the effect of vacancy location on the rupture progress of defective CNTs, with the results indicating that even a one-atom vacancy defect can reduce the failure stress of a CNT by as much as 26 per cent and that the reduction in failure strength caused by a vacancy defect is much larger than that caused by a Stone–Wales defect. A comparison (Belytschko *et al*. 2002) between the MD simulation results and the experimental results (Yu *et al*. 2000) of the fracture strain of a zigzag CNT has shown reasonably good agreement.

Other researchers have discussed the use of continuum mechanics approaches in the analysis of defective CNTs. For example, Zhang *et al*. (2005) employed a coupled MM–continuum mechanics technique to investigate the effects of single- and multi-atom vacancy defects on the fracture strength of CNTs. Their results demonstrated that these vacancy defects reduce the fracture strength by 20–33%. Under the assumption of a high defect concentration, an analytical expression (Sammalkorpi *et al*. 2004) was obtained to calculate the Young's modulus of CNTs with vacancy defects. Wang *et al*. (2007) used an elastic shell model to study the fracture strain of defective CNTs. They introduced a coefficient *α* to represent the effect of a single vacancy defect. Xiao & Hou (2007) proposed a multi-scale model for resonant oscillators by combining MD simulation and the continuum mechanics method. They then investigated the effects of vacancy defects on the behaviour of CNT-based oscillators. Walgraef (2007) proposed a nonlinear Donnell's shallow-shell model for the dynamic analysis of defective CNTs, but assumed defect density.

Liu *et al*. (2004, 2005) recently have developed an atomic-scale finite element method (AFEM) that is as accurate as MD simulation. The AFEM calculates CNTs using the framework of the conventional finite element method and is thus much faster than MD simulation. Leung *et al*. (2006) and Guo *et al*. (2007) used the AFEM to investigate the buckling and post-buckling behaviours of CNTs by incorporating it into the Abaqus commercial software package. However, because the AFEM performs its calculations on the assembled stiffness matrix, it cannot ensure that the stiffness matrix is positive definite for a large CNT system.

Atomistic and continuum mechanics methods are widely applied in the analysis of CNTs. However, the computational time and effort needed for the former limit the size of the CNTs that can be studied, and the latter cannot be used to study the mechanical properties of defective CNTs because they assume that materials are defect free. Hence, to allow the study of large CNTs with defects, in this paper, we develop an atomic-based cellular automata algorithm (ACAA). This algorithm enables the large-scale calculation of CNTs to be performed. We use the Tersoff–Brenner potential (Brenner *et al*. 2002) to describe multi-body interaction and apply the cellular automata (CA) algorithm (Wolfram 1986) for the computation. In CNT simulation, the CNT is treated as a discrete system with *n* nodes. Each atom is placed at a node, and an elementary calculation element is formed using a central node surrounded by the nearest three nodes and the second-nearest six nodes. Instead of assembling the global stiffness matrix, the proposed ACAA calculates the solution on the elemental level. Thus, only three simultaneous equations need to be solved, meaning that the computer capacity required is very low. Simulation on the elemental level also means that the ACAA model is suitable for parallel simulation and hence is a very promising method for the simulation of large CNT systems.

## 2. Theoretical formulation

As stated, the Tersoff–Brenner potential (Brenner *et al*. 2002) is employed in our ACAA simulation. The potential *U*(** x**) stored in the atomic bonds is expressed as(2.1)where

*V*

_{R}and

*V*

_{A}are the repulsive and attractive pair interactions, respectively;

*b*

_{ij}is the reactive empirical bond order between atoms;

*N*is the number of total atoms in the system; refers to the positions of

*N*atoms; and

*r*

_{ij}is the distance between atoms

*i*and

*j*. Each carbon atom interacts with both its nearest and second-nearest neighbouring atoms. Total energy

*Π*(

**) is(2.2)where**

*u***={(**

*F**F*

_{x},

*F*

_{y},

*F*

_{z})

_{1}, (

*F*

_{x},

*F*

_{y},

*F*

_{z})

_{2}, …, (

*F*

_{x},

*F*

_{y},

*F*

_{z})

_{N}}

^{T}is the external force vector;

**={(**

*u**u*

_{x},

*u*

_{y},

*u*

_{z})

_{1}, (

*u*

_{x},

*u*

_{y},

*u*

_{z})

_{2}, …, (

*u*

_{x},

*u*

_{y},

*u*

_{z})

_{N}}

^{T}is the displacement vector of all atoms;

**=**

*x*

*x*^{0}+

**; and**

*u*

*x*^{0}={(

*x*

^{0},

*y*

^{0},

*z*

^{0})

_{1}, (

*x*

^{0},

*y*

^{0},

*z*

^{0})

_{2}, …, (

*x*

^{0},

*y*

^{0},

*z*

^{0})

_{N}}

^{T}is the initial equilibrium position of

*N*atoms.

The total energy request of the equilibrium configuration of the system is minimal, and thus we have(2.3)Expressing *Π*(** u**) using the Taylor expansion to the second order in term

**and then substituting it into equation (2.3) yields(2.4)where**

*u***is the stiffness matrix, which is expressed as(2.5)and**

*K***is a non-equilibrium force vector(2.6)Owing to nonlinearity, equation (2.4) is solved iteratively until**

*P***reaches zero. Liu**

*P**et al*. (2004, 2005) proposed an AFEM to solve equation (2.4). However, the AFEM needs to assemble the element stiffness matrix in the global stiffness matrix and solve the simultaneous equilibrium equations in each iteration. For a large system, with thousands or even millions of atoms, assembling the global stiffness matrix and solving its inverse requires a tremendous amount of memory. Thus, the AFEM's application in the solution of large CNT systems is limited. The ACAA, by contrast, directly calculates the maximum 30×30 matrix for the unknown equilibrium configuration without forming the global stiffness matrix. Hence, a computer that has the capacity to calculate a 30×30 matrix can be used for the simulation of large CNTs.

## 3. Atomic-based CA algorithm

In this section, an ACAA is established for the simulation of CNTs. In the present algorithm, a cellular automaton consists of a hexagonal grid of cells and the state of a cell is determined by the previous state of the surrounding neighbourhood of cells. The cellular space is shown in figure 1, in which atoms are taken as cells and the grid is a two-dimensional hexagon. The cellular state is described as , where *i* and *k* denote the *i*th cell and the *k*th iteration step, respectively, in a discrete system of space and time; *u*_{x}, *u*_{y} and *u*_{z} are the displacements of cell *i* in the *x*, *y* and *z* directions, respectively; and *P*_{x}, *P*_{y} and *P*_{z} are the non-equilibrium forces in the *x*, *y* and *z* directions, respectively.

The first core part of a CA model is the definition of the cellular neighbourhood. Figure 2 shows the definition of this neighbourhood for the simulation of a CNT: it consists of a central node surrounded by the nearest three nodes and the second-nearest six nodes. The state of the cellular neighbourhood indicates that the energy stored in bond not only depends on its length, but also on the angles between bond and bonds , , and . Similarly, the energy stored in bond depends on both its length and the angles between bond and bonds , , and , and the energy stored in the bond depends on both its length and the angles between bond and bonds , , and . Hence, changing the position of node *N*_{i} has an effect on the energy stored in the surrounding nine bonds. This basic CNT structure indicates the nature of the multi-body interaction, and thus it is adopted as a cellular neighbourhood for CNT simulation. The cellular neighbourhood at node *i* is denoted as *D*_{i} and defined as(3.1)where *N*_{i} is the number of node *i* and (*l*=1, 2, …, 9) are the numbers of neighbour nodes in the neighbourhood *D*_{i}.

The second core part of a CA model is the definition of updated rules for CNT simulation. These updated rules are proposed according to the equilibrium condition. The carbon atoms of a CNT are taken as cells; thus, the CNT is treated as a discrete spatial structure, as shown in figure 1. The state of a cellular neighbourhood is determined by node force *P*_{i} and the previous states of a surrounding neighbourhood of cells. First, all of the neighbouring nodes are fixed artificially, and then the displacement increment Δ*u*_{i} of node *N*_{i} and the reaction increments of the neighbouring nodes are solved from equation (2.4). Note that the constraints on the neighbouring nodes are the ‘dummy constraints’. Thus, the negative reactions are applied to the neighbouring nodes to transfer the non-equilibrium forces to the entire domain of the cellular space. The real solution of equation (2.4) can be obtained by repeating the solution procedure until the displacement increment Δ*u*_{i} approaches zero. At the *k*th iteration step, the node displacement increments and the force increments of node *N*_{i} can be calculated as follows.

First, to capture the multi-body atomistic interactions, the Tersoff–Brenner's potential is employed as the potential to form the stiffness matrix of cell *i* and the node force vector . According to equations. (2.5) and (2.6), the stiffness matrix and the node force vector are expressed, respectively, as(3.2)and(3.3)It should be noted that the stiffness matrix used in this study is different from the element stiffness matrix used in Liu *et al*. (2004, 2005). Therefore, the local equilibrium equation of cell *N*_{i} can be written as(3.4)or as the matrix form(3.5)where refers to the displacement increments of node *N*_{i}, to the force increments of the neighbouring nodes, and(3.6)(3.7)and(3.8)Note that because all of the neighbouring nodes are assumed to be fixed. Then, the local equilibrium equations can be simplified as(3.9)Hence, there are at most three equilibrium equations, which are easily solved. When the displacement increments are obtained through equation (3.9), the force increments exerted on the neighbouring nodes are calculated by(3.10)where is a 27×3 matrix. Once is obtained, the displacement of node *N*_{i} at the (*k*+1)th iterative step is updated as(3.11)Accordingly, the node forces exerted on the neighbouring nodes are updated as(3.12)The iteration ends when all of the cells *N*_{i}(*i*=1, 2, …, *N*) have been calculated. The maximum absolute error between two adjacent iterations is taken as the convergence condition(3.13)where er is the tolerance.

Equations (3.9) and (3.10) show that, rather than assemble the global stiffness matrix, the ACAA model calculates the solution on the elemental level, as shown in figure 2. Thus, only the equations with a 3×3 matrix need to be solved, and the computer capacity required is very low. As the calculation is only on the elemental level, the ACAA model is suitable for parallel simulation and hence is a very promising method for the simulation of large CNT systems.

## 4. Results and discussion

To ensure the accuracy of the proposed ACAA model, we calculate the average energy per atom for a perfect (8, 0) zigzag single-walled carbon nanotube (SWCNT) with a length of 4.07 nm under axial compression and compare it with the quantum-generalized tight-binding MD (GTBMD) method (Srivastava *et al*. 1999) and MD simulations (Liew *et al*. 2004; Xiao *et al*. 2004), as shown in figure 3. The results obtained using the ACAA are in good agreement with those obtained using the GTBMD method and MD simulations.

Having validated the proposed ACAA to a certain extent, we examine the role of vacancy defects in SWCNT buckling. In all the examples, the SWCNT vacancy defect is formed as a hole by removing six atoms and the associated bonds from the lattice. One end of the SWCNT is completely fixed. The other end is free in the axial direction, but fixed in the cross section. Axial compression is achieved by applying a compressive displacement on one end of the SWCNT, while the other end is fixed.

The first example is a cantilever (20, 20) SWCNT that consists of 80 crystal cells and has length *L*=19.7869 nm. The vacancy defect is, respectively, are located at *L*/4, 2*L*/5 and *L*/2 from the fixed end, as shown in figure 4*a*. The total force on a section is plotted as the function of the strain in figure 4*b* for a perfect SWCNT and a defective SWCNT. Compared with the simulation of the perfect SWCNT, where the critical buckling strain is *ϵ*_{cr}=0.032, the critical buckling strains of defective SWCNTs are smaller and almost the same as *ϵ*_{cr}=0.031 for various defect locations at *L*/4, 2*L*/5 and *L*/2, and therefore the latter collapses earlier than the former. The vacancy defect that locates at *L*/4, 2*L*/5 and *L*/2 decreases the buckling load of the defective SWCNT by approximately 7.3, 5.4 and 4.7 per cent, respectively, in comparison with the perfect SWCNT, and the critical buckling load decreases as the vacancy defect approaches the fixed end. Figure 4*a* shows the buckling deformations of an SWCNT at strain *ϵ*_{cr}=0.031. The buckling configurations show that all of the buckling deformations are consistent with shell-buckling modes. The walls are severely wrinkled at the defect locations, implying that a vacancy defect acts as a geometrical imperfection. In addition, it can be seen that the force–strain relationships are not linear for either the perfect or the defective SWCNTs. The slopes of the force–strain curves of the defective SWCNTs are slightly smaller than those of the perfect SWCNT. Thus, it is demonstrated that vacancy defects can reduce the strength of SWCNTs slightly.

Figure 5 depicts the buckling and post-buckling behaviours of defective (20, 20) SWCNTs with length *L*=19.7869 nm. The two vacancy defects are located at *L*/2 of the SWCNTs with relative angles of 45°, 90°, 135° and 180°. The buckling configurations are shown in figure 5*a*. Again, these indicate that the buckling of a defective SWCNT is local buckling, as in the buckling mode of cylindrical shells, and the buckling occurs at the defect sites due to the imperfection caused. Figure 5*b* shows that an increase in the number of defects lowers an SWCNT's buckling load. The critical buckling strain is *ϵ*_{cr}=0.029 for both relative angles 90° and 180°. Compared with the SWCNT with one vacancy defect at *L*/2, the reduction of the buckling load for the relative angles 90° and 180° is 5.8 and 6.6 per cent, respectively. However, for both relative angles 45° and 135°, the critical buckling strain is *ϵ*_{cr}=0.03, with a 9.2 per cent reduction in the buckling load. This indicates that the effect of the relative location of defects along the circumferential direction cannot be ignored.

In contrast to the previous examples of (20, 20) SWCNTs with length/diameter ratio of 7.3, this example considers a (5, 5) SWCNT with a large length/diameter ratio of 290 to examine the buckling and post-buckling behaviours. The defective (5, 5) SWCNT consists of 800 crystal cells (*L*=197.869 nm) with one vacancy defect locating at *L/*2. The plot of the strain energy per atom of the defective SWCNT against the strain *ϵ* is given in figure 6. It is observed that the strain energy increases with the increase of strain until the critical strain *ϵ*=0.052. Shortly after *ϵ*=0.052, there is an abrupt drop in the strain energy, and the SWCNT starts to buckle at the defect site, as shown in figure 7*b*. After the buckle, the SWCNT enters into post-buckling stage with a hinge being formed at the defect site. As the strain increases, the strain energy increases almost linearly until the second critical strain *ϵ*=0.06 is reached, after which the strain energy slightly decreases, and the second hinge is formed. Further increase of the deformation beyond *ϵ*=0.06 results in an increase of the energy until the SWCNT collapses suddenly at the strain *ϵ*=0.086.

Figure 7 presents the morphologies of buckling and post-buckling of the defective SWCNT at various strains from *ϵ*=0.052 to *ϵ*=0.086. Instead of a shell-like buckling mode, the SWCNT with a large length/diameter ratio exhibits beam-like buckling character, as shown in figure 7. Similarly, the first buckling occurs at the location of the vacancy defect, as shown in figure 7*b*, indicating again that the vacancy defect leads to an initial geometrical imperfection. After the buckle, a hinge starts to form at the location of the vacancy defect as the strain increases, as shown in figure 7*c*. Then, the hinge has been formed and another hinge starts to form at the critical strain *ϵ*=0.06, as shown in figure 7*d*. Finally, figure 7*e* shows that the second hinge has been formed at *ϵ*=0.07. It is observed that two hinges are formed at the post-buckling stage, and the first hinge is located at the defect site.

The corresponding force–strain relationship is presented in figure 8. It is evident that the force–strain curve behaves nonlinearly with an increasing gradient even at a small strain level. Similarly, tiny drops in the force occur at critical strains *ϵ*=0.052 and *ϵ*=0.06 due to the release of energy. At the maximum strain *ϵ*=0.086, the force decreases rapidly to almost zero, indicating that the SWCNT loses its capacity to withstand the force.

## 5. Conclusions

In this paper, we develop an ACAA that can be used for atomic-scale simulation. The ACAA employs the inter-atomic potential to account for multi-body interactions and is as accurate as the MM method, but much faster because it is based on the solution framework of the matrix methods of structural analysis. Specifically, the use of the ACAA effectively avoids the problem of a non-positive-definite matrix in the simulation of large CNT system because the maximum size of the stiffness matrix is only 30 rows×30 columns. Although we consider the simulation of CNTs, it should be noted that the proposed method covers general topics of other nanomaterials such as fullerenes C_{60}, grapheme sheets and even CNT-reinforced composites.

We use the ACAA to study the buckling behaviour of defective CNTs and show that vacancy defects lower the critical buckling load significantly. The location of a vacancy defect and the relative location of two vacancy defects are examined to determine their effects on buckling behaviour. As expected, the buckling load depends on the location of the vacancy defect. The ACAA simulations also show that the buckling occurs at the defect sites, indicating that a vacancy defect causes a significant geometrical imperfection in CNTs.

## Acknowledgments

The work described in this paper was fully supported by a research grant from the Research Grants Council of the Hong Kong Special Administrative Region, China [Project no. CityU 113406]. The authors are grateful for this financial support.

## Footnotes

- Received June 5, 2008.
- Accepted September 5, 2008.

- © 2008 The Royal Society