Size dependency and potential field influence on deriving mechanical properties of carbon nanotubes using molecular dynamics

K.G.S.Dilrukshi，M.A.N.Dewpriy∗，U.G.A.Puswewl

aDepartment of Civil Engineering，University of Moratuwa，Moratuwa，Sri Lanka

bSchool of Engineering Science，Simon Fraser University，Burnaby，Canada

Size dependency and potential field influence on deriving mechanical properties of carbon nanotubes using molecular dynamics

K.G.S.Dilrukshia，M.A.N.Dewapriyab，∗，U.G.A.Puswewalaa

aDepartment of Civil Engineering，University of Moratuwa，Moratuwa，Sri Lanka

bSchool of Engineering Science，Simon Fraser University，Burnaby，Canada

A R T I C L E I N F O

Article history:

Accepted 11 March 2015

Available online 29 May 2015

Carbon nanotubes

Molecular dynamics

Potential field

Shear modulus

Poisson’s ratio

A thorough understanding on the mechanical properties of carbon nanotube（CNT）is essential in extendingtheadvancedapplicationsofCNTbasedsystems.However，conductingexperimentstoestimate mechanical properties at this scale is extremely challenging.Therefore，development of mechanistic models to estimate the mechanical properties of CNTs along with the integration of existing continuum mechanics concepts is critically important.This paper presents a comprehensive molecular dynamics simulationstudyonthesizedependencyandpotentialfunctioninfluenceofmechanicalpropertiesofCNT.Commonly used reactive bond order（REBO）and adaptive intermolecular reactive bond order（AIREBO）potential functions were considered in this regard.Young’s modulus and shear modulus of CNTs are derived by integrating classical continuum mechanics concepts with molecular dynamics simulations.The results indicate that the potential function has a significant influence on the estimated mechanical properties of CNTs，and the influence of potential field is much higher when studying the torsional behaviour of CNTs than the tensile behaviour.

Crown Copyright©2015 Published by Elsevier Ltd on behalf of The Chinese Society of Theoretical and Applied Mechanics.This is an open access article under the CC BY-NC-ND license（http://creativecommons.org/licenses/by-nc-nd/4.0/）.

Since the discovery by Iijima［1］，carbon nanotubes（CNTs）have attracted considerable attention due to their remarkable mechanical，electrical and thermal properties.In particular，CNTs have shown a broad range of applications in nanotechnology，especially in the development of nanosensors，nanomechanical resonators and as reinforcement in composite materials［2］.However，in order to further extend the CNT-based applications of nanotechnology，experimental demonstrations of versatility of CNTshavetobesupportedinparallelbytheoretical/computational models of CNT-based systems.In this regard，accurate knowledge on constitutive properties such as Young’s modulus（E），shear modulus（G），and Poisson’s ratio（ν）of CNT is essential to describe their mechanical behaviour under different loading，boundary and environmental conditions.However，the use of experiments at this scale to measure constitutive properties is extremely challenging and only few studies have been reported［2］.These limitations of experimental investigations can be overcome by developing mechanistic models to estimate the mechanical properties of CNTs along with the integration of CNT mechanics with existing continuum mechanics concepts.

Atomistic systems such as CNT can be accurately modelled by using the first principles methods.However，due to high computational cost，the first principles methods are only applicable to systems with several hundreds of atoms.Continuum modelling is the most computationally efficient method.Those are，however，unable to account for quantum effects of matter at nanoscale，and also the discrete nature of matter at nanoscale［3］.Combination of continuum concepts with properties derived using atomistic analysis provides an attractive method to analyse the systems at the nanoscale［4-6］.In this regard，molecular dynamics（MD）plays an important role in studying the nanoscale systems.

Many researchers have used atomistic modelling approaches to determine tensile strength，elastic modulus，shear modulus，and Poisson’s ratio by simulating tension and torsion tests of CNT［7-9］.The general approach involves the determination of potential energy of a CNT using MD for different strain levels and then computationofthestress-straincurveusingnumericalquadrature methods.However，the properties of CNT reported in the literature indicate different ranges of values depending（0.8-1.2 TPa）on the CNT diameter and methods of analysis［10-14］.The different trends of size dependency of CNT’s shear moduli reported in the literature are significant among these.For example，Yu et al.［7］have reported that shear modulus substantially increases with increasing diameter，whereas，experiments and several theoretical models show an opposite trend［15，16］.Some authors haveestimated the shear modulus of CNT from Young’s modulus and Poisson’s ratio obtained from MD simulation of tension tests and using the classical isotropic elasticity relationship between E，ν，and G.These results differ from the shear modulus estimated from MD simulation of a pure torsion test of CNT.Therefore，it is important to examine the reasons for the different trends of sizedependency of CNT moduli reported in the literature and differences in shear modulus values reported from different approaches（i.e.，tensionvs.torsion）.Inouropinion，advancedmodellingofCNT systems has limited value if the reasons for the above disagreements are not well understood.Furthermore，accurate estimation of elastic moduli is critical to the continuum analysis of CNT-based systems and devices.

Thispaper investigatesthesizedependencyand potentialfunction influence on the mechanical properties of CNT derived from MD simulations.Nanoscale MD simulations of CNTs subjected to direct tension and torsion loading were used in this regard，where thecommonlyusedreactivebondorder（REBO）andadaptiveintermolecular reactive bond order（AIREBO）potential functions representthemolecularinteraction.MDsimulationswereperformedfor different chiralities and diameters to understand the influence of potential field on the size dependency of the moduli.In addition，the validity of common isotropic elastic relationship between E，ν，and G was also examined.Also，the size and potential field dependence of torsional strengths was investigated.

REBO potential（also called Tersoff-Brenner potential）was originally developed to simulate chemical vapour deposition of diamond［17］，and later was modified to provide more accurate treatment of the energetic，elastic and vibration properties of solid carbonandsmallhydrocarbons（secondgenerationpotential）［18］.

In REBO potential，the energy stored in the bond between atom i and atom jis given as

wherebijisthebondordertermwhichmodifiesthestrengthofthe bond depending on the local bonding environment，andare repulsive and attractive potentials， respectively，rijis the distance between atoms i and j，and fis called the cut-off function.The purpose of the cut-off function is to limit the interatomic interactions to the nearest neighbours［18］.

Even though REBO potential is successful in describing the intermolecular interactions in carbon and hydrogen materials［17］，it is not appropriate for every hydrocarbon system due to its inability to explicitly capture non-bonded interactions and torsional interactions.By addressing these shortcomings of REBO potential，AIREBO potential has been developed［19］，which is an extension of the REBO potential

The original cut-off function for the REBO potential is given by

where R（1）and R（2）are the cut-off radii，which have the values of 1.7and 2 Å，respectively.The values of cut-off radii are defined based on the first and second nearest neighbouring distances of hydrocarbons.

It has been observed that the cut-off function could cause non-physical strain hardening in stress-strain curves of carbon nanostructures［20，21］whereas experiments［22］and ab-initio calculations［23］do not show any strain hardening.Therefore，researchers have modified the cut-off radii ranging from 1.9 Å to 2.2 Å［24，25］to eliminate this non-physical strain hardening.However，none of the previous studies have given much insight into the effect of cut-off function on the stress-strain relation，and this non-physical behaviour continues to prevail in recent simulation studies［21，26］.

To fully understand the influence of the cut-off function，and also to identify the best cut-off radii，a set of uniaxial tensile test was performed in this study on a 4.8 nm long（10，10）CNT using LAMMPS MD simulator［27］.The stress-strain curve obtained with defaultcut-offfunction（giveninEq.（3））indicatesstrainhardening at strain value around 0.16 as shown in Fig.1.In order to remove this non-physical behaviour，various modified cut-off functions were tested［6］.The results show that a modified cut-off function，fm（r）giveninEq.（4）eliminatesthestrainhardeningwhenthecutoff radius（R）is 2 Å.

Fig.1 also shows a stress-strain curve based on the modified cut-off function，which qualitatively agrees with the behaviour observed in experiments［20］and ab-initio calculations［23］.In order to obtain further insight into the effect of cut-off function on fracture of an individual bond，the force-strain curve of the carbon-carbon（C-C）bond between atoms 1 and 2 in Fig.2 was studied by increasing the bond length r.The positions of other atoms，relative to atoms 1 and 2，were kept unchanged.As shown in Fig.2，the strain hardening of the force-strain curve disappears when both cut-off radii are equal to 2 Å，since bond breaking takes place before bond length reaches 2 Å.It should be noted that the modified cut-off function in Eq.（4）is not continuous as opposed to the original cut-off function in Eq.（3）.However，this discontinuity does not affect the simulation results since C-C bond length does not reach the cut-off radius of the modified function（i.e.，2 Å）.Therefore，the modified cut-off function given in Eq.（4）was used in the ensuing MD simulations.

Fig.1.Stress-strain curves of the（10，10）CNT.

Fig.2.Force-strain curves of the C-C bond between atoms 1 and 2 shown in inset（arrows indicate the strain direction）.

To study the effects of potential field on the mechanical properties obtained from MD simulations，armchair（10，10）and zigzag（17，0）CNTs were subjected to tension and torsion tests.The simulations were repeatedly done using REBO and AIREBO potential fields.Initial investigations of the effect of aspect ratio（length to diameter ratio）for the simulations revealed that the effect is negligible if the aspect ratio is greater than 12.Therefore，considering the computational efficiency，CNTs with aspect ratio of 15 were used for the ensuing simulations.Tensile and torsional loadings were applied on CNTs by changing the coordinates of atoms accordingly.Atoms at the two edges were fixed in position for both tensile and torsional simulations.For tensile loading，extension was applied at an interval of 1.25 Å（length of CNT is～200 Å），and CNTs are allowed to reach equilibrium over 40000 time steps（20 ps）.During torsion test simulations，a rotation angle of 5°was applied at a step，and then allowed to reach equilibrium over 80000 time steps（40 ps）.The equilibrium potential energy wasobtainedbyaveragingtheenergyover10000timesteps（5ps），after reaching the equilibrium configuration.Fig.3 shows（10，10）and（17，0）CNTs subjected to 50°rotation.The first 30 Å of CNTs from the fixed edge is shown.Only half of CNTs are shown for clarity.The axes of CNTs are along y-direction.

Variations of strain energy per atom of（10，10）and（17，0）CNTs subject to tension and torsion are shown in Figs.4 and 5， respectively.It can be seen that irrespective of the chirality，potential energy obtained using AIREBO potential is higher than that of REBO potential in tensile test.As opposed to the REBO，the AIREBO potential considers energy from non-bonded interactions and also energy from torsional interactions.Therefore，the energy given by AIREBO potential is higher compared to the energy given by REBO.Change in strain energy of CNTs due to torsional（shear）strain does not have considerable effect from potential field used for the simulation.However，a predominant effect was noted after the initiation of torsional buckling（which can be identified from the point of inflection of the graph），as shown in Fig.5，the rate of strain energy change reduced with the AIREBO potential，whereas，it was constant with the REBO potential.The torque and shear strain variation，which was derived from the energy-strain variation，is shown in Fig.6.

Fig.4.Variation of potential energy with axial strain for CNTs subjected to uniaxial tensile test.

Fig.3.（10，10）and（17，0）CNTs after been subjected to 50°rotation.

Thetorque-strain curvesshow thatirrespectiveof the chirality，torsional buckling of CNT occurred earlier with the AIREBO potential than REBO potential.This may be due to the ability of AIREBO potential to capture the change in energy due to torsional deformations，which gives ability to estimate energy corresponding to initiation of buckling accurately compared to REBO.Further，it was noticed that for the CNT of same diameter，armchair CNTs are subjected to torsional buckling earlier than the armchair CNT（see Fig.6）.Also，AIREBO potential indicated slightly higher torsional strength for zigzag CNTs than that of armchair，whereas REBO potential indicated the opposite trend.

Fig.5.Variation of potential energy with shear strain for CNTs subjected to torsion test.

Fig.6.Variation of torque with shear strain of CNTs.

The relationships between the strain energy and the modulus values，defined in structural mechanics，were used to estimate both elastic and shear moduli of CNTs.Since the wall thickness of CNTs is not accurately known（researchers have used various wall thicknesses ranging from 0.06 nm to 6 nm），the surface modulus values were considered.

The relationship between strain energy（U）and axial strain（ε）for a thin cylindrical shell of radius（a），length（L），and wall thickness（t）can be expressed as

whereE istheelasticmodulus.Itshouldbenotedthateventhough CNTs are nonlinear，the behaviour at small strains（up to～0.04）can be approximated as linear［20，21］.Therefore，U-εrelationship up to the axial strain of 0.04 was conceded in estimating elastic modulus values in the study.

Also，for a CNT subjected to tensile loading，Poisson’s ratio can also be defined as

Fig.7.Variation ofv2with strain.

where r0and rεare the initial radius and radius under axial strain ofε.This relation was used to study the strain dependency of Poisson’s ratio（ν）.The value ofνis calculated at the strains where potential energies are obtained in the uniaxial tensile test of CNTs.The value of r0is obtained at the middle length of CNTs（6 nm from edges）after allowing CNTs to relax over 15 ps.Initial radius（r0）of（10，10）and（17，0）CNTs are 6.68 nm and 6.58 nm，respectively.Also，the variation of r0with potential field was negligible.Similar method was used to calculate the values of rε.Fig.7 shows the variation ofνwith strain.It can be observed in the figure that νvalue of（10，10）CNT reduces with increasing strain，whereas the value of（17，0）CNT increases with strain.Also，Poisson’s ratio values obtained by AIREBO potential are always higher than the of REBO potential.

Calculated elastic modulus and Poisson’s ratio can be used to estimate the shear modulus by

Torsional test results were also used to estimate shear modulus values.In this case，the potential energy variation of CNTs before the initiation torsional buckling was used.The relationship between potential energy（U），twisting angle（φ），and torsional stiffness（K）can be expressed as

The torsional stiffness for a thin cylindrical shell of radius a，length L，wall thickness t，and shear modulus G is defined by K = 2πa3Gh/L.Hence，surface shear modulus（Gt）can be obtained as

Values calculated for surface elastic modulus，Poisson’s ratio and surface shear modulus for CNTs using two potential functions for（17，0）and（10，10）CNTs are given in Table 1.

Itcanbeseenthat，irrespectiveofchiralitysurfaceelasticmodulus（Et）aswellasthesurfaceshearmodulusvalues（G1t）estimated using REBO potential are higher than the values estimated using AIREBO potential.The effect is significant for（17，0）CNTs with respect to surface elastic modulus and for the armchair CNTs with respect to surface shear modulus.However，the effect is higher for surface shear modulus.This can be attributed to the differentgeometric variations related to the CNTs of different chirality under tensile and torsional strains and ability of capturing the energy variation with respect to them by the two potential functions.

Average Poisson’s ratio values obtained for CNTs subjected to axial strain up to 0.04 during the tensile test simulations are given in Table 1.It can be noted that the potential field also has significant effect on Poisson’s ratio.The effect for（17，0）CNT（represents zigzag）is higher than that of（10，10）CNT（represents armchair）.The values obtained for Poisson’s ratio were further usedtocalculatesurfaceshearmodulusvaluesusingthestructural mechanics relationship as per Eq.（7）（G2t）.As seen in Table 1，the difference between two values（G1t and G2t）is considerably less with AIREBO than the REBO.

Table 1Poisson’s ratio estimated from REBO and AIREBO potentials.

Table 2Variation of elastic and shear moduli with diameter.

Values for elastic modulus and shear modulus were estimated for CNT’s with different diameters and variations are shown in Fig.8.AIREBO potential is used for the simulations.The summary of the results is given in Table 2.

It can be seen that the value of the shear modulus decreases with the increase of the diameter，which is similar to the experimental observations of Salvetat et al.［15］.Also，elastic moduli of CNTs do not show a significant variation with the diameter if the diameter is greater than 1 nm.Irrespective of the potential field，the values of shear modulus of zigzag and armchair tubes are almost similar at particular diameter，whereas elastic moduli of zigzag CNTs are higher than that of armchair.

In conclusion，the influence of potential field on the estimated mechanical properties is much higher when studying the torsional behaviourofCNTsthanthetensilebehaviour.Whenestimatingthe torsional strength and also studying the torsional behaviour after the initiation of torsional buckling，the effect of potential field is considerable，and the AIREBO potential can produce more realistic results than the REBO potential.This should be attributed to the ability of AIREBO potential to compute the torsional interactions between atoms，which is included into the total energy.Geometry（or chirality）of the CNT also plays a key role in this regard.However，irrespectiveofchiralitysurfaceelasticmodulus（Et）aswellas the surface shear modulus values（G1t）estimated using REBO potentialarehigherthanthevaluesestimatedusingAIREBOpotential and effect seems higher with shear modulus.Also，if the structural mechanics relationships are used to estimate shear modulus using known values of shear modulus and Poisson’s ratio，AIREBO potentialcanproducemoreaccurateresults（closetothevaluesobtained by direct energy method）than the REBO potential.

Fig.8.Variation of elastic and shear moduli with CNT diameter.

Acknowledgements

Authors thank Dr.Nimal Rajapakse for fruitful discussions and useful comments.This work was financially supported by National ScienceFoundation（NSF）ofSriLankaandtheNaturalSciencesand Engineering Research Council（NSERC）of Canada.WestGrid and Compute/Calcul Canada provided computing resources.

References

［1］S.Iijima，Helical microtubules of graphitic carbon，Nature 354（1991）56-58.

［2］R.S.Ruoff，D.Qian，W.K.Liu，Mechanical properties of carbon nanotubes: theoretical prediction and experimental measures，C.R.Physique 4（2003）993-1008.

［3］L.Tapasztó，T.Dumitrică，S.J.Kim，et al.，Breakdown of continuum mechanics for nanometre-wavelength rippling of graphene，Nat.Phys.8（2012）739-742.

［4］M.A.N.Dewapriya，A.S.Phani，R.K.N.D.Rajapakse，Influenceoftemperatureand free edges on the mechanical properties of graphene，Modelling Simul.Mater.Sci.Eng.21（2013）065017.

［5］M.A.N.Dewapriya，R.K.N.D.Rajapakse，A.S.Phani，Atomistic and continuum modelling of temperature-dependent fracture of graphene，Int.J.Fract.187（2014）199-212.

［6］M.A.N.Dewapriya，Molecular dynamics study of effects of geometric defects onthemechanicalpropertiesofgraphene（MSthesis），TheUniversityofBritish Columbia，Canada，2012.

［7］W.Yu，W.X.Xi，N.Xianggui，Atomic simulation of the torsional deformation of carbon nanotubes，Modelling Simul.Mater.Sci.Eng.12（2004）1099.

［8］A.Montazeri，M.Sadeghi，R.Naghdabadi，et al.，Computational modelling of the transverse isotropic elastic properties of single-walled carbon nanotube，Comput.Mater.Sci.49（2010）544-551.

［9］A.Montazeri，J.Javadpour，A.Khavandi，et al.，Mechanical properties of multiwalled carbon nanotube/epoxy composites，Mater.Des.31（2010）4202-4208.

［10］O.Gülseren，T.Yildirim，S.Ciraci，Systematicabinitiostudyofcurvatureeffects in carbon nanotubes，Phys.Rev.B.65（2002）153405.

［11］D.Sánchez-Portal，E.Artacho，J.M.Soler，et al.，Ab initio structural，elastic，and vibrational properties of carbon nanotubes，Phys.Rev.B.59（1999）12678.

［12］S.Xiao，W.Hou，Studies of size effects on carbon nanotubes’mechanical properties by using different potential functions，Fullerenes，Nanotubes Carbon Nanostruct.14（2006）9-16.

［13］B.W.Xing，Z.C.Chun，C.W.Zhao，Simulation of Young’s modulus of singlewalled carbon nanotubes by molecular dynamics，Physica B 352（2004）156-163.

［14］J.Y.Hsieh，J.M.Lu，M.Y.Huang，et al.，Theoretical variations in the Young’s modulus of single-walled carbon nanotubes with tube radius and temperature:a molecular dynamics study，Nanotechnology 17（2006）3920.

［15］J.P.Salvetat，G.A.D.Briggs，J.M.Bonard，etal.，Elasticandshearmoduliofsinglewalled carbon nanotube ropes，Phys.Rev.Lett.82（1999）944-947.

［16］T.Natsuki，K.Tantrakarn，M.Endo，Prediction of elastic properties for singlewalled carbon nanotubes，Carbon 42（2004）39-45.

［17］D.W.Brenner，Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond，Phys.Rev.B 42（1990）9458.

［18］D.W.Brenner，O.A.Shenderova，J.A.Harrison，et al.，A second-generation reactive bond order（REBO）potential energy expression for hydrocarbons，J.Phys.:Condens.Matter.14（2002）783-802.

［19］S.J.Stuart，A.B.Tutein，J.A.Harrison，A reactive potential for hydrocarbons with intermolecular interactions，J.Appl.Phys.112（2000）6472.

［20］O.A.Shenderova，D.W.Brenner，A.Omeltchenko，et al.，Atomistic modeling of the fracture of polycrystalline diamond，Phys.Rev.B 61（2000）3877.

［21］A.J.Cao，J.M.Qu，Atomisticsimulationstudyofbrittlefailureinnanocrystalline graphene under uniaxial tension，Appl.Phys.Lett.102（2013）071902.

［22］B.Peng，M.Locascio，P.Zapol，et al.，Measurements of near-ultimate strength for multiwalled carbon nanotubes and irradiation-induced crosslinking improvements，Nature Nanotechnol.3（2008）626-631.

［23］F.Liu，P.Ming，J.Li，Abinitiocalculationofidealstrengthandphononinstability of graphene under tension，Phys.Rev.B 76（2007）064120.

［24］H.Zhao，N.R.Aluru，Temperature and strain-rate dependent fracture strength of graphene，J.Appl.Phys.108（2010）064321.

［25］M.A.N.Dewapriya，R.K.N.D.Rajapakse，Molecular dynamics simulations and continuum modeling of temperature and strain rate dependent fracture strength of graphene with vacancy defects，J.Appl.Mech.81（2014）081010.

［26］L.Wang，Q.Zhang，Elastic behavior of bilayer graphene under in-plane loadings，Curr.Appl.Phys.12（2012）1173-1177.

［27］S.J.Plimpton，Fast parallel algorithms for short-range molecular dynamics，J.Comput.Phys.117（1995）1-19.

28 October 2014

∗.

E-mail address:mandewapriya@sfu.ca（M.A.N.Dewapriya）.

http://dx.doi.org/10.1016/j.taml.2015.05.005

2095-0349/Crown Copyright©2015 Published by Elsevier Ltd on behalf of The Chinese Society of Theoretical and Applied Mechanics.This is an open access article under the CC BY-NC-ND license（http://creativecommons.org/licenses/by-nc-nd/4.0/）.

*This article belongs to Solid Mechanics