Molecular Dynamics Simulation of Carbon Nanotubes

Molecular Dynamics Simulation of Carbon Nanotubes

1Sumit Sharma*, 2Rakesh Chandra, 3Pramod Kumar, 4Navin Kumar

1Research scholar, Department of Mechanical Engineering

Dr. B. R. Ambedkar National Institute of Technology, Jalandhar, India

2Professor, Department of Mechanical Engineering

Dr. B. R. Ambedkar National Institute of Technology, Jalandhar, India

3Associate Professor, Department of Mechanical Engineering

Dr. B. R. Ambedkar National Institute of Technology, Jalandhar, India

4Professor, School of Mechanical, Materials & Energy Engineering (SMMEE) Indian Institute of Technology, Ropar, India

Abstract

Elastic properties of single walled carbon nanotubes (SWCNTs) have been determined using molecular dynamics (MD) simulation. Mechanical properties of three types of SWCNTs viz., armchair, zigzag and chiral nanotubes have been evaluated. From computational results, it can be concluded that the Youngs moduli of SWCNTs decrease with increase in radius of SWCNT and increase with increase in CNT volume fractions (Vf) and aspect ratios (l/d).

Keywords: A. Carbon nanotubes; A. Nano composites; A. Short-fiber composites; C. Elastic properties.

1. Introduction

   Carbon nanotubes (CNTs) were first reported by Iijima [1] in 1991. Since then, CNTs have been attracting much attention to explore their exceptional electronic and material properties. Due to their large aspect ratios and small diameters, CNTs have emerged as potentially attractive materials as reinforcing elements in lightweight and high strength structural composites. As a one-dimensional structure, CNTs can be thought of as one sheet or multiple sheets of graphene rolled into a cylinder. There are single or multiple layers of carbon atoms in the tube thickness direction, called single-walled carbon nanotubes (SWCNTs) or multi-walled

   carbon nanotubes (MWCNTs), respectively. According to different chiral angles, SWCNTs can be classified into zigzag (=0), armchair (=30) and chiral tubule (0 < < 30).

2. Literature review

   The determination of Youngs modulus for CNTs has been a subject of considerable interest. The computation of Youngs modulus of CNTs may be classified into two categories. One is molecular dynamics (MD) simulation using a potential energy function obtained by empirical, tight-binding or ab initio methods. The other approach relies on the development of models based on molecular and continuum mechanics. In experimented approaches, the force displacement response of a nanotube is measured and the axial Youngs modulus is obtained by comparison to an equivalent elastic beam. Treacy et al., [2] showed an average value of 1.8 TPa (with large scatter) for the axial Youngs modulus from the direct measurements with a transmission electron microscope of a variety of multi-walled nanotubes (MWNTs) of different inner and outer diameters using the thermal vibration analysis of anchored tubes. The nanotubes with the smallest inner diameter were considerably stiffer, with a Youngs modulus of 3.70 TPa.

   Lourie and Wagner [3] obtained the axial Youngs modulus for a series of temperatures by micro-Raman spectroscopy from measurements of cooling-induced compressive deformation of nanotubes embedded in an epoxy matrix. At 81K, the experimental results gave 3 TPa for single-walled nanotubes (SWNTs) with an average radius of 0.7 nm, and 2.4 TPa for MWNTs with an average radius of 510 nm. Wong et al., [4] used an atomic force microscope (AFM) to measure forcedisplacement relations for anchored MWNTs on a substrate. They obtained the Youngs modulus by comparing their results with elastic beam theory. An average of 1.28 ± 0.5 TPa with little dependence of nanotube diameter was reported. Lu [5] used an empirical force- constant model to determine several elastic moduli of single- and multi-walled nanotubes and

   obtained the Youngs modulus of about 1 TPa and the rotational shear modulus of about 0.5 TPa. The analysis showed that the elastic properties were insensitive to the radius, helicity, and the number of walls. However, Yao and Lordi [6] used MD simulations and found that changes in structure such as radius and helicity of the SWNTs could affect the Youngs modulus because their results showed that the torsional potential energy, which is the dominant component of total potential energy, increased as the quadratic function of the decreasing tube radius. Zhou et al.,

   [7] claimed that both Youngs modulus and the wall thickness were independent of the radius and the helicity of SWNTs. They applied the strain energy of SWNTs directly from electronic band structure without introducing empirical potentials and continuum elasticity theory to describe the mechanical properties of SWNTs. The estimated value for the axial modulus was reported as 5.0 TPa, which is 5 times larger than the value of MWNTs.

   Liu et al., [8] reported the Youngs modulus of CNTs is 1-0.1 TPa with the diameter increasing from 8 to 40 nm by measuring resonance frequency of carbon nanotubes. Krishnan et al., [9] used TEM to observe the vibration of an SWCNT at room temperature and reported Youngs modulus of SWCNTs in the range from 0.90 to 1.70 TPa, with an average of 1.25 TPa. Tombler et al., [10] used AFM to bend an SWCNT and reported the Youngs modulus of SWCNTs around 1.2 TPa. Yu et al., [11] conducted nanoscale tensile tests of SWCNT ropes pulled by AFM tips under a scanning electron microscope and reported that the Youngs modulus of SWCNT ropes ranged from 0.32 to 1.47 TPa. Demczyk et al., [12] reported that the Youngs modulus of MWCNTs range from 0.8 to 0.9 TPa when TEM is used to bend an individual tube.

   Bao et al., [13] predicted the Youngs modulus of SWCNTs and graphite based on molecular dynamics (MD) simulation. The inter-atomic short-range interaction and long-range

   interaction of carbon nanotubes have been represented by a second generation reactive empirical bond order (REBO) potential and LennardJones (LJ) potential, respectively. The obtained potential expression has been used to calculate the total potential energies of carbon nanotubes. From the simulation, the Youngs moduli of SWCNTs are weakly affected by the tube chirality and tube radius. The numeric results are in good agreement with the existing experimental results. Youngs moduli of SWCNTs are in the range of 929.8 ± 11.5 GPa. The average of the Youngs modulus of graphite is 1026.176 GPa. Yang and Wei [14] investigated the mechanical properties of nano-single crystal gold and carbon nanotube-embedded gold (CNT/Au) composites under uni-axial tension and reported the Youngs modulus of the nano-single crystal gold as 66.22 GPa. Maximum yield stress has been reported to be 5.74 GPa at a strain of 0.092. The increase in Youngs modulus of long CNT-embedded gold composite over pure gold has been found to be very large.

   In spite of the variety of theoretical studies on the macroscopic elastic behavior of CNTs, there still remain controversial issues regarding the effect of geometric structure of CNTs on elastic moduli, as evidenced by the wide scatter among the elastic moduli reported in the literature. The objective of this paper is to reexamine the elastic behavior of CNTs in detail using MD simulation. The organization of this paper is as follows. In the next section, the morphological structure of carbon nanotubes will be briefly discussed. Section 4 provides the force fields and total potential energy that are related to the interatomic potentials for MD simulations. In addition, the bonding and nonbonding terms in the ttal potential energy are described. In Section 5, several elastic moduli are determined by applying different small-strain deformation modes. The elastic moduli are predicted using energy and force approaches. Numerical results and a summary are given in Sections 6 and 7.

3. SWCNT structure

   Single-walled nanotubes are formed by folding a grapheme sheet to form a hollow cylinder which is composed of hexagonal carbon ring units, which are also referred to as graphene units. The fundamental carbon nanotube structure can be classified into three categories: armchair, zigzag, and chiral, in terms of their helicity. Figure 1 shows a segment of single graphite plane that can be transformed into a carbon nanotube by rolling it up into a cylinder. To describe this structure, a chiral vector is defined as OA=na1+ma2, where a1 and a2 are unit vectors for the honeycomb lattice of the graphene sheet, n and m are two integers, along with a chiral angle which is the angle of the chiral vector with respect to the x direction shown in Figure 1. This chiral vector, OA, will be denoted as (n, m) which will also specify the structure of the carbon nanotube. Vector OB is perpendicular to the vector OA. To construct a CNT, we cut off the quadrangles OA B'B and roll it into a cylinder with OB and AB' overlapping each other. The relationship between the integers (n, m) and the nanotube radius, r, and chiral angle,

   is given by:

   = 3 (2 + + 2)1 2 /2 (1)

   = tan1 3 /( + 2) (2)

   Where, is the length of the C-C bond.

4. Molecular dynamics simulation methodology

   The concept of the MD method is rather straightforward and logical. The motion of molecules is generally governed by Newtons equations of motion in classical theory. In MD simulations, particle motion is simulated on a computer according to the equations of motion. If one molecule moves solely on a classical mechanics level, a computer is unnecessary because mathematical calculation with pencil and paper is sufficient to solve the motion of the molecule.

   However, since molecules in a real system are numerous and interact with each other, such mathematical analysis is impracticable. In this situation, therefore, computer simulations become a powerful tool for a microscopic analysis. If the mass of molecule i is denoted by mi, and the force acting on molecule i by the ambient molecules and an external field denoted by fi, then the motion of a particle is described by Newtons equation of motion:

   2

   2 = (3)

   If a system is composed of N molecules, there are N sets of similar equations, and the motion of N molecules interacts through forces acting among the molecules. Differential equations such as Eq. (3) are unsuitable for solving the set of N equations of motion on a computer. Computers readily solve simple equations, such as algebraic ones, but are quite poor at intuitive solving procedures such as a trial and error approach to find solutions. Hence, Eq. (3) will be transformed into an algebraic equation. To do so, the second-order differential term in Eq. (3) must be expressed as an algebraic expression, using the following Taylor series

   expansion:

   ()

   1 2 2()

   1 3 3 ()

   +

   =

   + + 2!

   2 + 3!

   3 + (4)

   Eq. (4) implies that x at time (t + h) can be expressed as the sum of x itself, the first- order differential, the second-order differential, and so on, multiplied by a constant for each term. If x does not significantly change with time, the higher order differential terms can be neglected for a sufficiently small value of the time interval h. In order to approximate the second-order differential term in Eq. (3) as an algebraic expression, another form of the Taylor series expansion is necessary:

   ()

   1 2 2()

   1 3 3()

   =

   + 2!

   2 3!

   3 + (5)

   If the first order differential term is eliminated from Eqs. (4) and (5), the second-order

   differential term can be solved as;

   2()

   + 2 + ( ) 2

   2 =

   2 + O

   (6)

   The last term on the right-hand side of this equation implies the accuracy of the approximation, and, in this case, terms higher than p are neglected. If the second order differential is

   approximated as;

   2()

   2 =

   + 2 + ( )

   2 (7)

   This expression is called the central difference approximation. With this approximation and the notation ri = (xi, yi, zi) for the molecular position and fi = (fxi, fyi, fzi) for the force acting on particle i, the equation of the x-component of Newtons equation of motion can be written as;

   +

   = 2

   2

   +

   (8)

   Similar equations are satisfied for the other components. Since Eq. (8) is a simple algebraic equation, the molecular position at the next time step can be evaluated using the present and previous positions and the present force. If a system is composed of N molecules, there are 3N algebraic equations for specifying the motion of molecules; these numerous equations are solved on a computer, where the motion of the molecules in a system can be pursued with the time variable. Eq. (8) does not require the velocity terms for determining the molecular position at the next time step. This scheme is called the Verlet method. A scheme using the positions and velocities simultaneously may be more desirable in order to keep the system temperature constant. Considering that the first and second-order differentials of the position are equal to the velocity and acceleration, respectively, and neglecting differential terms of higher order in Eq. (4);

   +

   =

   +

   2

   + 2

   (9)

   This equation determines the position of the molecules, but the velocity term arises on the right- hand side, so that another equation is necessary for specifying the velocity. The first-order

   differential of the velocity is equal to the acceleration;

   + = + (10)

   In order to improve accuracy, the force term in Eq. (10) is slightly modified and the following

   equation is obtained;

   + = + ( + + ) (11)

   2

   Scheme of using Eq. (9) and Eq. (11) for determining the motion of molecules is called the velocity Verlet method.

5. Total potential energies and inter-atomic forces

   The reliability of MD simulation technique depends on the use of appropriate inter- atomic energies and forces. In the context of molecular modeling force field refers to the form and parameters of mathematical functions used to describe the potential energy of a system of particles (typically molecules and atoms). Force field functions and parameter sets are derived from both experimental work and high-level quantum mechanical calculations. In this study, we have used the Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies (COMPASS) forcefield. This forcefield is a member of the consistent family of force fields (CFF91, PCFF, CFF and COMPASS), which are closely related second-generation force fields. They were parameterized against a wide range of experimental observables for organic compounds containing H, C, N, O, S, P, halogen atoms and ions, alkali metal cations, and several biochemically important divalent metal cations. COMPASS is the first force field that has been

   parameterized and validated using condensed phase properties in addition to empirical data for molecules in isolation. Consequently, this force field enables accurate and simultaneous prediction of structural, conformational, vibrational, and thermo-physical properties for a broad range of molecules in isolation and in condensed phases.

   The COMPASS force field consists of terms fr bonds (b), angles (), dihedrals (), out- of-plane angles () as well as cross-terms, and two non-bonded functions, a Coulombic function for electrostatic interactions and a 9-6 Lennard-Jones potential for van der Waals interactions.

   Etotal = Eb + E + E + E + Eb,b' + Eb, + Eb, + E, + E,' + E,', + Eq + EvdW (12) Where,

   = 2 0 2 + 3 0 3 + 4 0 4

   = 2 0 2 + 3 0 3 + 4 0 4

   (13)

   (14)

   = 1 1 + 2 1 2 + 3 1 3

   (15)

   = 22

   (16)

   0

   0

   , = 0 ( ) (17)

   , = 0 ( 0 )

   ,

   (18)

   , = 0 1 + 2 2 + 3 3

   ,

   (19)

   , = 0 1 + 22 + 3 3

   ,

   (20)

   0

   0

   , = 0 ( )

   (21)

   ,

   0

   0

   , , = 0

   (22)

   , ,

   =

   (23)

   0 9 0 6

   = 2

   3

   (24)

   Where,

   k, k1, k2, k3 and k4 = force constants determined experimentally

   b, = bond length and bond angle after stretching and bending respectively b0, 0 = equilibrium bond length and equilibrium bond angle respectively

   = bond torsion angle

   = out of plane inversion angle

   Eb,b' , E,' , Eb, , Eb, , E, , E,', = cross terms representing the energy due to interaction between bond stretch-bond stretch, bond bend-bond bend, bond stretch-bond bend, bond stretch- bond torsion, bond bend-bond torsion and bond bend-bond bend-bond torsion respectively.

   i, j= well depth or bond dissociation energy

   ij

   ij

   r0 = distance at which the interaction energy between the two atoms is zero rij= separation between the atoms/molecules

   qi, qj = atomic charges on the atoms/molecules 0= permittivity of free space

6. Stiffness of SWCNTs

   The stiffness of SWCNTs having different chirality have been calculated using Materials

   Studio 5.5 MD software. The basic steps in calculation of stiffnesses of different types of CNTs have been explained in the following sub-sections.

   1. Modeling of SWCNTs

      The first step is to model the SWCNTs using Build tool in Materials Studio. We can construct SWCNTs having different chirality (n,m). In this study we have built three types of SWCNTs viz. zigzag (n,0), armchair (n,n) and chiral (n,m) nanotubes. Here, the integer n controls the overall size of the nanotube. The minimum value for n is 1. The integer m controls the chiral angle or twist of the graphite sheet used to construct the nanotube. The minimum value for m is 0. Some of the models constructed have been shown in Figures 2-4.

   2. Geometry optimization

      A frequent activity in molecular dynamics simulation is the optimization or minimization (with respect to potential energy) of the system being examined. For instance it is often desirable to optimize a structure after it has been sketched, since sketching often creates the molecule in a high energy configuration and starting a simulation from such an un-optimized structure can lead to erroneous results. There are a number of optimization techniques available in Materials Studio viz., steepest descent method, conjugate gradient and newton-raphson method.

      In the steepest descents method, the line search direction is defined along the direction of the local downhill gradient. Each line search produces a new direction that is perpendicular to the previous gradient; however, the directions oscillate along the way to the minimum. This inefficient behavior is characteristic of steepest descents, especially on energy surfaces having narrow valleys. Convergence is slow near the minimum because the gradient approaches zero, but the method is extremely robust, even for systems that are far from being harmonic. It is the method most likely to generate the true low-energy structure, regardless of what the function is

      or where the process begins. Therefore, the steepest descents method is often used when the gradients are large and the configurations are far from the minimum. This is commonly the case for initial relaxation of poorly refined crystallographic data or for graphically built models. The reason that the steepest descents method converges slowly near the minimum is that each segment of the path tends to reverse progress made in an earlier iteration. It would be preferable to prevent the next direction vector from undoing earlier progress. This means using an algorithm that produces a complete basis set of mutually conjugate directions such that each successive step continually refines the direction toward the minimum. If these conjugate directions truly span the space of the energy surface, then minimization along each direction in turn must, by definition, end in arrival at a minimum. The conjugate gradient algorithm

      constructs and follows such a set of directions. As a rule, N 2 independent data points are required

      to solve a harmonic function with N variables numerically.

      Since a gradient is a vector with N variables, the best we can hope for in a gradient-based minimizer is to converge in N steps. However, if we can exploit second-derivative information, an optimization could converge in one step, because each second derivative is an × matrix. This is the principle behind the variable metric optimization algorithms, of which Newton- Raphson is perhaps the most commonly used. Another way of looking at Newton-Raphson is that, in addition to using the gradient to identify a search direction, the curvature of the function (the second derivative) is also used to predict where the function passes through a minimum along that direction. Since the complete second-derivative matrix defines the curvature in each gradient direction, the inverse of the second-derivative matrix can be multiplied by the gradient to obtain a vector that translates directly to the nearest minimum.

      In this study we have used the smart algorithm which is a cascade of the above stated methods. The parameters used in the optimization of the nano-structures have been shown in Table 1.

   3. Dynamics

      Once an energy expression and, if necessary, an optimized structure have been defined for the system of interest, a dynamics simulation can be run. The basis of this simulation is the classical equations of motion which are modified, where appropriate, to deal with the effects of temperature and pressure on the system. The main product of a dynamics run is a trajectory file that records the atomic configuration, atomic velocities and other information at a sequence of time steps which can be analyzed subsequently. Different parameters used in dynamics run have been listed in the Table 2.

   4. Mechanical properties

      We have used the "Forcite" module to calculate the mechanical properties of SWCNTs. The Forcite mechanical properties task allows us to calculate mechanical properties for a single structure or a trajectory of structures. Forcite mechanical properties calculation may be performed on either a single structure or a series of structures generated, for example, by a dynamics run and stored in a trajectory file (.arc, .his, .trj, .xtd). The mechanical properties are then calculated using the classical simulation theory, averaged over all valid configurations, and reported in the output text document. Anybody or element thereof, which is acted on by external forces is in a state of stress. Moreover, if the body is in equilibrium, the external stress must be exactly balanced by internal forces. In general, stress is a second rank tensor with nine

      components as follows:

      11 12 13

      = 21 22 23 (25)

      31 32 33

      In an atomistic calculation, the internal strss tensor can be obtained using the so-called virial expression:

      = 1 ( ) +

      (26)

      0

      =1

      <

      where index i runs over all particles 1 through N; mi vi and fi denote the mass, velocity and force acting on particle i; and V0 denotes the (un-deformed) system volume. The application of stress to a body results in a change in the relative positions of particles within the body, expressed

      quantitatively via the strain tensor:

      11 12 13

      = 21 22 23 (27)

      31 32 33

      The elastic stiffness coefficients, relating the various components of stress and strain are defined

      by:

      1 2

      =

      ,

      = 0

      , ,

      (28)

      where A denotes the Helmholtz free energy. For small deformations, the relationship between the stresses and strains may be expressed in terms of a generalized Hooke's law:

      = (29)

      For calculating the mechanical properties of SWCNTs, the parameters shown in Table 3 have been used.

7. Results and discussion

   In this section, the results obtained for SWCNTs, have been discussed in detail. Several models of armchair, zigzag and chiral SWCNTs have been constructed. The arm-chair SWCNTs have been shown in Figure 5 and Figure 6. Variation of temperature with time for (10,10) SWCNT has been shown in Figure 7. This has been obtained using the Forcite module in

   Material Studio software. The dynamics run has been made for 5 ps. Two structures of zigzag SWCNTs have been shown in Figure 8 and Figure 9. Dynamics run for (10,0) SWCNT has been shown in Figure 10.

   Chiral SWCNTs have been shown in Figure 11 and Figure 12. Dynamics run showing the variation of temperature with time for (18,8) SWCNT has been shown in Figure 13. We have also calculated the moduli of clusters of SWCNTs. Simulation cell for cluster of seven CNTs has been shown in Figure 14. Similarly, Figure 15 and Figure 16 show the simulation cell for clusters of nine CNTs and nineteen CNTs respectively. Variation of Young's modulus (E11) with radius for different types of SWCNTs has been shown in Figure 17. Yao and Lordi [15] explored the dependence of Youngs modulus (Y) on diameter and helicity by using the universal force field developed by Rappe et al., [16]. It has been observed that Y decreases significantly with increasing tube diameter, e.g., Y for a (20,20) tube is 15% smaller than that for a (5,5) tube. They also noted the correlation between the variation of Y and the torsional energy with diameter and helicity. Yao and Lordi [15] also discussed the difference between such dependence of Y noted by them and near independence of Y on diameter and helicity as noted by Lu [17] on the basis of the absence of torsion energy or four-atom interaction energy term in the computations of Lu [17]. Thus, the results of Yao and Lordi [15] suggest that the diameter dependence observed in the present work is not surprising as Brenners potential [18] used here has four-atom interaction terms. All these observations suggest that either Y is constant with diameter or the variation of Y with diameter is small and it depends on some other factors such as four-atom interaction terms. Depending upon the strength of such terms at the configuration of the atoms of the tube, the computed results might show a decrease or increase or constancy of Y with diameter.

   The carbon atoms in a CNT are in sp2 configurations and connected to one another by three strong bonds. Due to the geometric orientation of the carbon-carbon bonds relative to the nanotube axis, armchair SWCNTs exhibit higher tensile strength and Youngs moduli values compared to the zigzag SWCNTs. In general, the Youngs moduli depend more on the radii than on the helicity. The Youngs modulus of chiral SWCNTs is slightly lower than that of armchair and zigzag SWCNTs as observed from Figure 17. The results obtained from MD simulations for SWCNTs have been tabulated in Table 4. The average of Youngs modulus, E11 of armchair SWCNTs is 602 GPa. The average transverse Youngs modulus, E22 of armchair SWCNTs is 190 GPa. Similarly, the average of Youngs modulus, E11 of zigzag SWCNTs is 599 GPa and the average transverse Youngs modulus, E22 of zigzag SWCNTs is 166 GPa. The average of Youngs modulus, E11 of chiral SWCNTs is 598 GPa and the average transverse Youngs modulus, E22 of chiral SWCNTs is 165 GPa.

   Figure 18 shows the variation of Young's modulus (22 ) with radius for different types of SWCNTs. The transverse Youngs modulus (22 ) also decreases with increase in radius of SWCNT. Also, the value of (22 ) is higher for armchair SWCNT in comparison to the zigzag and chiral SWCNTs. Table 5 shows the results obtained for Bulk modulus, shear modulus and poisson's ratio of SWCNTs for different tube radius. The average Bulk modulus of armchair SWCNTs is 94.2 GPa and of zigzag SWCNTs is 91.5 GPa. The average Bulk modulus of chiral SWCNTs is 134.82 GPa. The average shear modulus of armchair SWCNTs is 68.51 GPa and of zigzag SWCNTs is 72.36 GPa. The average shear modulus of chiral SWCNTs is 58.66 GPa, whereas the poissons ratio is the same (0.30) for both armchair and zigzag SWCNTs, the poissons ratio for chiral SWCNTs is 0.22. Figure 19 shows the variation of Bulk modulus (K) with radius for different types of SWCNTs. For diameters less than 12 Ã…, the average bulk

   modulus (K) of armchair SWCNTs is greater than that of zigzag nanotubes. For diameters greater than 12 Ã…, the average bulk modulus (K) of zigzag SWCNTs is greater than that of armchair nanotubes. Also, it can be observed from Figure 19 that for diameters greater than 12 Ã…, the average bulk modulus (K) of chiral SWCNTs is the greatest.

   Figure 20 shows that the shear modulus (G) decreases with increase in radius of SWCNTs. It can be inferred from Figure 20 that the shear modulus of armchair SWCNTs is greater than that of zigzag SWCNT. For diameters greater than 12 Ã…, the shear modulus is almost the same for all three types of SWCNTs. Figure 21 shows the variation of Poissons ratio with radius for different types of SWCNTs. Here also, the value of poissons ratio is the greatest for armchair SWCNT and also the poissons ratio decreases with increase in diameter of SWCNTs. Table 6 shows the simulation conditions for a cluster of (7,0) SWCNTs. The length of (7,0) SWCNT has been kept constant in each case. Conjugate gradient method has been used for geometry optimization and Universal force field has been used for dynamics run. Other parameters for MD simulation have been shown in the Table 6.

   Table 7 shows the MD simulation results for a cluster of (7,0) SWCNTs. Results show that the moduli decrease with increase in the number of SWCNTs. It confirms the fact that the clusters of SWCNTs have lower transverse properties in comparison to individual nanotubes. Figure 22 shows the variation of Youngs modulus (E11) with number of carbon nanotubes. In clusters of SWCNTs, the inter-tube force interactions are primarily due to non-bonding weak Van der Waals interactions. These weak cohesive properties of nanotube bundles cause the Youngs modulus (Y) to decrease with increase in the number of SWCNTs. Figure 23 shows the variation of transverse modulus (E22), shear modulus (G23) and bulk modulus (K23) with number of SWCNTs. The transverse modulus is greater than shear and bulk modulus.

8. Conclusions

Using MD simulations, we have evaluated the Youngs moduli of armchair, zigzag and chiral SWCNTs. Results show that the average of Youngs modulus, E11 of armchair SWCNTs is 602 GPa. The average transverse Youngs modulus, E22 of armchair SWCNTs is 190 GPa. Similarly, the average of Youngs modulus, E11 of zigzag SWCNTs is 599 GPa and the average transverse Youngs modulus, E22 of zigzag SWCNTs is 166 GPa. The average of Youngs modulus, E11 of chiral SWCNTs is 598 GPa and the average transverse Youngs odulus, E22 of chiral SWCNTs is 165 GPa. Main findings of the study have been listed below:

1. It can be concluded that the Youngs moduli of SWCNTs decrease with increase in radius and also with the number of SWCNTs. The shear moduli, bulk moduli and poissons ratio also decrease with increase in radius of SWCNTs.

2. For clusters of SWCNTs, as the number of SWCNTs is increased the moduli decrease. The SWCNTs exhibit poor shear properties and excellent longitudinal properties.

Our study illustrates that the simulation of Youngs modulus of SWCNTs with MD is quiet workable and reliable. This study will be further extended to study the effect of functionalization of CNTs on their mechanical properties.

References

1. Iijima S. Helical microtubules of graphitic carbon. Nature 1991; 354:56-8.

2. Treacy MMJ, Ebbesen TW, Gibson JM. Exceptionally high Youngs modulus observed for individual carbon nanotubes. Nature 1996; 381: 678-80.

3. Lourie O, Wagner HD. Evaluation of Youngs modulus of carbon nanotubes by micro- Raman spectroscopy. Materials Research 1998; 13: 2418-22.

4. Wong EW, Sheehan PE, Lieber CM. Nanobeam mechanics: elasticity, strength, and toughness of nanorods and nanotubes. Science 1997; 277:1971-5.

5. Lu JP. Elastic properties of carbon nanotubes and nanoropes. Physical Review Letters 1997; 79:1297-300.

6. Yao N, Lordi V. Youngs modulus of single-walled carbon nanotubes. Applied Physics 1998; 84:1939-43.

7. Zhou X, Zhou J, Ou-Yang Z. Strain energy and Youngs modulus of single-wall carbon nanotubes calculated from electronic energy-band theory. Physical Review B 2000; 62:13692-6.

8. Liu JZ, Zheng QS, Jiang Q. Effect of a rippling mode on resonances of carbon nanotubes.

   Physical Review Letters 2001; 86:4843-6.

9. Krishnan A, Dujardin E, Ebbesen TW, Yianilos PN, Treacy MMJ. Young's modulus of single-walled nanotubes. Physical Review B 1998; 58:14013-19.

10. Tombler TW, Zhou CW, Alexseyev L, Kong J, Dai HJ, Lei L. Reversible electromechanical characteristics of carbon nanotubes under local-probe manipulation. Nature 2000; 405:769- 72.

11. Yu MF, Files BS, Arepalli S, Ruoff RS. Tensile loading of ropes of single wall carbon nanotubes and their mechanical properties. Physical Review Letters 2000; 84:5552-5.

12. Demczyk BG, Wang YM, Cumings J, Hetman M, Han W, Zettl A, Ritchie RO. Direct mechanical measurement of the tensile strength and elastic modulus of multi-walled carbon nanotubes. Materials Science and Engineering A 2002; 334:173-8.

13. Bao WX, Zhu CC, Cui WZ. Simulation of Youngs modulus of single-walled carbon nanotubes by molecular dynamics. Physica B 2004; 352:156-63.

14. Yang SH, Wei ZX. Molecular dynamics study of mechanical properties of carbon nanotube- embedded gold composites. Physica B 2008: 403:559-63.

15. Yao N, Lordi V. Young's modulus of single-walled carbon nanotubes. Applied Physics 1998; 84:1939-43.

16. Rappe AK, Casewit CJ, Colwell KS, Goddard WA, Skiff WM. UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations. American Chemical Society 1992; 114:10024-39.

17. (a) Lu JP. Elastic Properties of Carbon Nanotubes and Nanoropes. Physical Review Letters 1997; 79:1297-1300; (b) Lu JP. The elastic properties of single and multilayered carbon nanotubes. Physics and Chemistry of Solids 1997; 58:1649-1652.

18. Brenner DW, Shenderova OA, Harrison JA, Stuart SJ, Ni B, Sinnott SB. A second generation reactive empirical bond order potential energy expression for hydrocarbons. Physics Condensed Matter 2002; 14:783-802.

19. Ostaz AA, Pal G, Mantena PR, Cheng A. Molecular dynamics simulation of SWCNT polymer nanocomposite and its constituents. Material Science 2008; 43:164-173.

Figure 1 Graphite plane of nanotube Figure 2 A zigzag (10,0) SWCNT. surface coordinates.

Figure 3 An armchair (20,20) SWCNT. Figure 4 A chiral (12,6) SWCNT.

Figure 5 An armchair (6,6) SWCNT. Figure 6 An armchair (10,10) SWCNT.

Figure 7 Dynamics run showing the variation of temperature with time for (10,10) SWCNT.

Figure 8 A zigzag (6,0) SWCNT. Figure 9 A zigzag (10,0) SWCNT.

Figure 10 Dynamics run showing the variation of temperature with time for (10,0) SWCNT.

Figure 11 A chiral (12,6) SWCNT. Figure 12 A chiral (18,8) SWCNT.

Figure 13 Dynamics run showing the variation of temperature with time for (18,8) SWCNT.

Figure 14 A simulation cell for seven Figure 15 A simulation cell for nine (7,0) SWCNTs. (7,0) SWCNTs.

Figure 16 A simulation cell for nineteen (7,0) SWCNTs.

Armchair

Zigzag

Chiral

Armchair

Zigzag

Chiral

1600

1400

Young's modulus ( E11), GPa

Young's modulus ( E11), GPa

1200

1000

800

600

400

200

0

0 2 4 6 8 10 12 14 16

Radius ()

Figure 17 Variation of Young's modulus (11 ) with radius for different types of SWCNTs.

400

350

300

250

200

150

100

50

0

0

5

10

Radius ()

15

20

Armchair

Zigzag

Chiral

400

350

300

250

200

150

100

50

0

0

5

10

Radius ()

15

20

Armchair

Zigzag

Chiral

Zigzag

Chiral

0

Zigzag

Chiral

0

Young's modulus ( E22), GPa

Young's modulus ( E22), GPa

Bulk modulus (K), GPa

Bulk modulus (K), GPa

Figure 18 Variation of Young's modulus (22 ) with radius for different types of SWCNTs.

Armchair

170

160

150

140

130

120

110

100

90

80

70

60

50

40

30

20

10

0

Armchair

170

160

150

140

130

120

110

100

90

80

70

60

50

40

30

20

10

0

Radius ()

Radius ()

2

2

4

4

6

6

8

8

10

10

12

12

14

14

16

16

Figure 19 Variation of Bulk modulus (K) with radius for different types of SWCNTs.

120

110

100

90

80

70

60

50

40

30

20

10

0

0 2 4 6 8 10 12 14 16

Radius ()

120

110

100

90

80

70

60

50

40

30

20

10

0

0 2 4 6 8 10 12 14 16/p>

Radius ()

Armchair

Zigzag

Chiral

Armchair

Zigzag

Chiral

Shear modulus (G), GPa

Shear modulus (G), GPa

Figure 20 Variation of Shear modulus (G) with radius for different types of SWCNTs.

Armchair

Zigzag

Chiral

Armchair

Zigzag

Chiral

0.4

0.35

Poisson's ratio (12)

Poisson's ratio (12)

0.3

0.25

0.2

0.15

0.1

0.05

0

0 2 4 6 8 10 12 14 16

Radius ()

Figure 21 Variation of poisson's ratio (12) with radius for different types of SWCNTs.

6

9

12

15

18

21

6

9

12

15

18

21

No. of nanotubes

No. of nanotubes

E11

E11

1200

1200

1000

1000

800

800

600

600

400

400

200

200

0

0

E22

G23

K23

E22

G23

K23

60

60

50

50

40

40

30

30

Young's modulus (E11), GPa

Young's modulus (E11), GPa

Modulus, GPa

Modulus, GPa

Figure 22 Variation of Young's modulus E11 with no. of SWCNTs.

6

9

No. of nanotubes

18

21

6

9

No. of nanotubes

18

21

20

20

10

10

0

0

12

12

15

15

Figure 23 Variation of modulus with no. of SWCNTs.

Table 1 Geometry optimization parameters for SWCNTs.

S.No.	Parameter	Value
1.	Algorithm	Smart
2.	Quality convergence tolerance	Fine
3.	Energy convergence tolerance	10-4 kcal/mol
4.	Force convergence tolerance	0.005 kcal/mol/
5.	Displacement convergence tolerance	5×10-5
6.	Maximum no. of iterations	500

Table 2 Dynamics run parameters for SWCNTs.

S.No	Parameter	Value
1.	Ensemble	NVT
2.	Initial velocity	Random
3.	Temperature	300 K
4.	Time step	1 fs
5.	Total simulation time	5 ps
6.	No. of steps	5000
7.	Frame output every	5000
8.	Thermostat	Andersen
9.	Collision ratio	1
10.	Energy deviation	5×1012 kcal/mol
11.	Repulsive cut-off	6

Table 3 Mechanical properties simulation parameters for SWCNTs.

S.No	Parameter	Value
1.	Number of strains	6
2.	Maximum strain	0.001
3.	Pre-optimize structure	Yes
4.	Algorithm	Smart
5.	Maximum number of iterations	10
6.	Forcefield	Compass
7.	Repulsive cut-off	6

Table 4 Young's modulus of SWCNTs for different tube radius.

(n,m)	Number of atoms	Radius ()	Length ()	Longitudinal Young's modulus (E11), GPa	Transverse Young's modulus (E22), GPa
Armchair
(2,2)	48	1.35	14.75	1455	367
(4,4)	192	2.71	29.51	970	302
(6,6)	408	4.07	41.81	652	220
(8,8)	768	5.42	59.08	587	194
(10,10)	1120	6.78	68.86	519	181
(12,12)	1680	8.13	86.08	440	178
(14,14)	2240	9.49	98.38	391	129
(16,16)	2880	10.85	110.67	360	118
(18,18)	3744	12.2	127.89	330	112
(20,20)	4560	13.56	140.19	317	98
Average				602	190
Zigzag
(2,0)	24	0.78	12.78	1050	275
(4,0)	80	1.56	21.30	850	263
(6,0)	168	2.35	29.82	650	214
(8,0)	288	3.13	38.34	560	193
(10,0)	400	3.92	42.60	550	167
(12,0)	624	4.7	55.38	500	153
(14,0)	840	5.48	63.9	498	108

(16,0)	1088	6.26	72.42	490	102
(18,0)	1368	7.04	80.94	430	94
(20,0)	1680	7.83	89.46	410	89
Average 599					166
Chiral
(12,6)	1176	6.21	78.89	750	265
(14,6)	1896	6.96	113.59	711	230
(16,6)	2328	7.71	125.86	690	219
(18,8)	3192	9.03	147.38	639	204
(20,8)	3744	9.78	159.62	610	180
(20,10)	3360	10.36	135.25	598	142
(24,11)	3844	12.13	132.06	556	130
(30,8)	4816	13.58	147.81	512	102
(30,10)	6240	14.11	184.31	498	93
(30,12)	6552	14.66	186.22	424	85
Average 598					165

Table 5 Bulk modulus, shear modulus and poisson's ratio of SWCNTs for different tube radius.

(n,m)	Number of atoms	Radius ()	Length ()	Bulk modulus (K), GPa	Shear modulus (GVoight), GPa	Poisson's ratio (12)
Armchair
(2,2)	48	1.35	14.75	160.46	113.62	0.34
(4,4)	192	2.71	29.51	106.32	89.01	0.36
(6,6)	408	4.07	41.81	84.56	84.70	0.34
(8,8)	768	5.42	59.08	75.59	61.20	0.29
(10,10)	1120	6.78	68.86	83.78	62.22	0.31
(12,12)	1680	8.13	86.08	89.9	58.07	0.32
(14,14)	2240	9.49	98.38	94.16	56.56	0.3
(16,16)	2880	10.85	110.67	89.93	56.94	0.29
(18,18)	3744	12.2	127.89	83.88	54.07	0.26
(20,20)	4560	13.56	140.19	73.43	48.78	0.27
Average 94.20					68.51	0.30
Zigzag
(2,0)	24	0.78	12.78	126.82	110.51	0.34
(4,0)	80	1.56	21.30	119.92	97.5	0.38
(6,0)	168	2.35	29.82	90.34	73.31	0.35
(8,0)	288	3.13	38.34	80.46	70.21	0.33
(10,0)	400	3.92	42.60	75.48	67.84	0.31
(12,0)	624	4.7	55.38	74.22	69.02	0.28
(14,0)	840	5.48	63.9	74.96	60.28	0.29
(16,0)	1088	6.26	72.42	84.59	60.2	0.3
(18,0)	1368	7.04	80.94	90.87	58.63	0.26
(20,0)	1680	7.83	89.46	97.37	56.19	0.27

Average 91.50					72.36	0.30
Chiral
(12,6)	1176	6.21	78.89	138.41	80.49	0.26
(14,6)	1896	6.96	113.59	128.98	75.12	0.28
(16,6)	2328	7.71	125.86	132.58	58.06	0.27
(18,8)	3192	9.03	147.38	138.68	56.33	0.22
(20,8)	3744	9.78	159.62	141.05	60.87	0.23
(20,10)	3360	10.36	135.25	147.49	54.78	0.19
(24,11)	3844	12.13	132.06	137.19	52.25	0.18
(30,8)	4816	13.58	147.81	130.66	51.43	0.19
(30,10)	6240	14.11	184.31	125.8	47.62	0.17
(30,12)	6552	14.66	186.22	127.44	49.7	0.173
Average 134.82					58.66	0.22

Table 6 Simulation conditions for cluster of (7,0) SWCNTs.

Simulation cell size (3)	Number of atoms	Geometry optimization parameters	Dynamics run parameters	Mechanical properties parameters
Seven nanotubes
26.48×26.48×105	4900	Algorithm: Conjugate gradient	Ensemble: NPT	No. of strains: 06
		Maximum no. of iterations: 1000	Pressure: 1MPa	Maximum strain: 0.0005
		Force field: Universal	Time step: 1 fs	Algorithm: Conjugate gradient
			No. of steps: 100000	Force field: Universal
			Temperature: 298 K	Maximum no. of iterations: 1000
Nine nanotubes
26.48×26.48×105	6300	Algorithm: Conjugate gradient	Ensemble: NPT	No. of strains: 06
		Maximum no. of iterations: 1000	Pressure: 1MPa	Maximum strain: 0.0005
		Force field: Universal	Time step: 1 fs	Algorithm: Conjugate gradient
			No. of steps: 100000	Force field: Universal

			Temperature: 298 K	Maximum no. of iterations: 1000
Nineteen nanotubes
43.86×43.86×105	13300	Algorithm: Conjugate gradient	Ensemble: NPT	No. of strains: 06
		Maximum no. of iterations: 1000	Pressure: 1MPa	Maximum strain: 0.0005
		Fore field: Universal	Time step: 1 fs	Algorithm: Conjugate gradient
			No. of steps: 100000	Force field: Universal
			Temperature: 298 K	Maximum no. of iterations: 1000

Table 7 MD simulation results of cluster of (7,0) SWCNTs.

No. of SWCNTs	Longitudinal Young's modulus (E11), GPa	Transverse Young's modulus (E22), GPa	Poisson's ratio (12)	Shear modulus (G23), GPa	Bulk modulus (K23), GPa
7	1055	41.25	0.31	10.75	54.07
9	1047	39.06	0.29	8.50	52.97
19	797	12.25	0.29	6.41	41.80