2012.MD simulation of carbon nanotube pullout behavior and its use

MD simulation of carbon nanotube pullout behavior and its use in determining mode I delamination toughnessLin Yang a ,b ,Liyong Tong b ,Xiaodong He a ,⇑a Centre for Composite Materials and Structures,Harbin Institute of Technology,Harbin 150080,ChinabSchool of Aerospace,Mechanical and Mechatronic Engineering,The University of Sydney,NSW 2006,Australiaa r t i c l e i n f o Article history:Received 31August 2011Received in revised form 7November 2011Accepted 10December 2011Available online 15January 2012Keywords:Carbon nanotube Molecular dynamics Traction law Pullout behaviora b s t r a c tThis study presents a multi scale approach of understanding the influence of multi-walled carbon nanotube (MWNT)on delamination toughness.Molecular dynamic (MD)simulation was used to investigate the influence of diameter and length of a MWNT and the boundary and loading conditions on its pullout behav-ior from polymer matrix.Numerical results reveal that the MWNT experiences non-uniform stretching and necking prior to the movement of its embedded end and pullout,and the force–displacement curve features with four distinct regions.A quadruple-linear model was proposed via a parametric study and validated by comparing with existing experimental result,and then used as a traction law in deriving the formula for determining the enhancement in mode I delamination toughness due to MWNT presence.Ó2011Elsevier B.V.All rights reserved.1.IntroductionDelamination in fiber-reinforced laminated composites has been identified as an important limiting factor in practical struc-tural applications.To enhance the delamination toughness and other poor through-the-thickness properties,various techniques,such as weaving,pinning,and stitching,have been attempted to place additional fiber tows of microscale in the thickness direction [1].Recently,carbon nanotubes (CNTs)have been used as nano-scale reinforcement entities,nanopins or stitches,by joining car-bon nanotubes (CNTs)to a fiber or ply in order to use their super properties to enhance the interfacial strength or out-of-plane resin toughness between two adjacent fibers or plies [2–7].The chemical vapor deposition (CVD)method was used to grow well-aligned multi-walled carbon nanotubes (MWNTs)of 60l m long on the 2D woven SiC fabric,and two such fabrics were then laminated to-gether to form MWNT reinforced laminated adjacent plies [2].The three-point bending test and DCB test results show a remarkable improvement in mode I delamination toughness [2].An alternative approach is the chemical modification to connect CNTs onto a car-bon fiber surface.Amido link and denderimer (macromolecule)link that connects CNTs onto a carbon fiber graphite sheet have been successfully realized [3].Fig.1depicts a scanning electron microscopy (SEM)micrograph of a carbon fiber surface with CNTs connected via denderimer (macromolecule)[3,4].It was shown inthe microdroplet experiment that the addition of the chemically-connected CNTs can lead to a noticeable improvement in the fiber’s interfacial shearing strength [3,8].Similar to a pin or stitch in a DCB test,a CNT can play the role of bridging the two delaminated or cracked surfaces with one end being connected to the SiC fiber and the other embedded in resin [2].This CNT bridging mechanism can be schematically illustrated in Fig.2,and the axial tensile force in a CNT together with its end displacement play a role of retarding crack propagation and absorbing deformation energy.The relationship between the axial pull force and its associated end displacement for an embedded CNT as schematically depicted in Fig.2c is often referred to as a traction law,and determination of this relationship is important in assessing the role of CNTs in affecting the interfacial shear strength and delamination toughness in CNT reinforced 3D carbon fiber composite materials [5–7].The interface between CNT and polymer could be controlled by non-bond forces or cross-links or a combination of both.For the first case,only non-bond forces,comprising of Coulombic electrostatic force and van der Waals force,exist between a non-chemically-mod-ified CNT and its surrounding polymer in composites.The TEM graphs [10,18]show that some CNTs pulled out from polymer have smooth surface,which indicates that the cross-links number be-tween CNT and polymer could be very low.For the second case,a chemical modification of the CNTs can achieve an enhanced compat-ibility to composite material via cross-links between CNT and ma-trix.In this case,CNTs coated with a polymer sheath are protruded from the composite fracture surface,indicating the adherence of the polymer to the nanotube based on the TEM images [11,23–25].0927-0256/$-see front matter Ó2011Elsevier B.V.All rights reserved.doi:10.1016/matsci.2011.12.014Corresponding author.E-mail address:xdhe@ (X.D.He).In fact,the adherence between CNT and polymer is strong enough to transfer the sliding mechanism of a CNT from its polymeric matrix to case of the polymer sheath being broken from resin matrix.The TEM images show that the CNTs coated with polymer sheaths could be entirely pulled out relatively easily,and thus leading to that the CNT pulling force is less than the CNT’s strength.As the distributionCNT-polymer cross-links can be random,this research focuses non-bond interaction of a non-chemically-modified CNT sliding ranging from440nm to600nm.This range is nearly equal to the length of the empty cylindrical hole(160nm)left after the MWNT’s partial pullout which represents the MWNT’s embedded end movement distance.It should be noted that the force-distance curve ranges from530nm to600nm drops back to contacting force with Kaptonfilm level which means the MWNT pulling force is low and ignorable.Cooper et al.[9]observed that the bending of the WMNTs before pullout is typically<6°.For the TEM images case,the MWNT crossing the hole in epoxy matrix is about 240nm,which,by multiplying sine(6°),is equal to25nm repre-senting the MWNT’s stretch elongation when considering bending. This distance is nearly the same as the SPM tip’s movement range from440nm to475nm during which the pulling force is slightly larger than contacting force with Kaptonfilm andfluctuating.As there are limited works studying the interaction between a CNT and polymer molecules by employing MD simulation,in this study a MD simulation approach was adopted to investigate the influences of various parameters on the interaction between a CNT and its surrounding resin molecules in a pullout process as schematically shown in Fig.2d.In the limited works available in the literature[13–15],used were relatively smaller molecular units containing shorter CNT and without considering the influences of pulling velocity,MWNT diameter and length in a uniform extrac-tion process.This study focuses on investigating MWNT’s extrac-tion with different pulling speeds,MWNT diameters and lengths.A simplified traction law was then established based on the MD re-1.An SEM image showing denderimer treated CNTs connected onto carbonfiber surface[4].of modelling the CNT bridging role:(a)interlaminar crack bridged by CNTs;(b)the crack tip where CNTs connect the upper leg and theyields pull-out of a CNT;and(d)a MD model.Science55(2012)356–364357side view of the overall model;(b)cross-sectional view of the model;(c)three-walled CNT;(d)rein matrixE pot¼E with-bondþE non-bondð1Þwhere E with-bond can be further decomposed intofive items,E bond, E angle,E dihedral,E improper,E coupling,and E non-bond comprises of two items:the Coulumbic electrostatic interaction energy E elec and the van der Waals’energy E vdw,which are[16]E non-bond¼E elecþE vdwð2Þwhere detail expressions of E ,E,E,E,E,to polymerize a three-dimensional network epoxy resin.We use the same method in[15]to form cross-links between epoxy resin and curing agent,the cross-links expanded in all directions in resin matrix.To ensure a reasonable molecular configuration,atoms are assumed to achieve equilibrium in4000time steps(1fs per time step)to simulate the curing process.After this equilibrium process, the epoxy resin molecules distribute uniformly around the MWNT as shown in Fig.3b.There is obviously an annular vacancy aroundDeformation configurations of three-walled CNT pull-out from polymer matrix of the stretching stage(a–d)and the extraction stageL.Yang et al./Computational Materials Science55(2012)356–364359360L.Yang et al./Computational Materials Science55(2012)356–3643.Basic observationsThis section presents the MD simulation results for a three-walled carbon nanotube being pulled out from an EPON862epoxy resin.These results were obtained using the following baseline parameters:MWCT embedded length=100Å,MWCNT external ra-dius=10.764Å,loading rate=0.2Å/fs,radius of the boundary cylin-der=24Å.In12,000simulation steps(0.08fs per time step),the MWNT segment were completely pulled out from the epoxy resin.3.1.MWNT’s deformation and axial force distributionFig.4depicts the MWNT’s deformations at the selected eight time steps.It is noted that the entire MWNT extraction process can be divided into two stages,i.e.,stage one featured with pro-gressive development of stretching,necking,and debonding until reaching full MWNT’s length(see Fig.4a–d),and stage two com-menced with sliding movement and then diameter enlargement of the MWNT’s embedded end(see Fig.4e–h).In stage one,as shown in Fig.4a–d,starting from the loading end,a portion of the MWNT is stretched,yielding necking due toPoisson’s contraction,and a combined stretching and necking in turn causes its debonding from epoxy resin.This continues until the entire MWNT is debonded.The MWNT elongation,necking and debonding from the matrix were observed in the TEM images of crack propagation in the MWNT–polystyrene thinfilms with some gaps between a MWNT and matrix when the MWNT was ex-tracted by the fracture surface[10].The phenomenon of MWNT shrinking was observed in the MWNT tensile tests[19].As shown in Fig.5b–d,the pulling force at the loading end increases quickly to its peak prior to debonding and then decreases gradually as deb-onding propagates along the MWNT.The axial force in the MWNT is positive but not uniformly distributed,and its undulation is be-lieved to be caused by the progressive debonding of MWNT from its surrounding resin via overcoming the Columbic electrostatic interaction force and van der Waals force[13].This oscillation re-flects the occurrence of unevenly debonding which creates stress waves transmitting along the MWNT’s longitudinal direction in the extracted segment in Fig.4b–d.In stage two,as shown in Fig.4e–h,the diameter near the embedded end of the MWNT enlarges remarkably compared to the reduced value.This is due to the creation of a segment in com-pression because the embedded end moves to the loading end at a rate faster than speed of MWNT exiting the epoxy resin.Fig.5e–h shows that there exists a segment of the embedded MWNT that is in compression with a compressive force of7l N at6000time step. The enlarged diameter in the embedded segment subsequently creates additional resistance to the pullout,which yields a sharp increase in the axial tensile force at the loading end as shown in Fig.5f and g.This MWNT compression and diameter enlargement phenomenon was observed in a TEM image of the buckled and compressed region of an MWNT[11].As shown in Fig.5e–g,the length of MWNT segment in tension reduces from approximately 140Åat4500time steps to90Åat7000time steps,and the axial force at the loading end increases from0.65l N to2.76l N.Fig.4e–h also shows that the embedded end exhibits a bubble shape,which indicates the inner walls move slower than the out-most one.This observation illustrates viscosity between the move-ments of neighboring walls.This phenomenon of relative slippage between the walls in a MWNT was observed in a high-resolution image of a MWNT pullout[18].Fig.6depicts the variation of the distance between cross sec-tions A and B near the embedded end of the MWNT as marked in Fig.3a.Fig.6shows that the segment length between cross sec-tions A and B experiences a small change up to3700time step,then varies from13Åto19Åfrom3700to4600time steps.This segment is in compression4600–9100time steps leading to en-larged diameter,and intension after9100time steps.The segment experiences a nominal strain of46%,possible due to MWNT’s super plastic phenomenon[20].3.2.Force–displacement curveFig.7depicts two curves of the axial force versus displacement of the MWNT at the loading end plotted as original MD and smoothed results.The axial force is defined by multiplying the average stress with the MWNT cross-sectional area at cross section near the loading region show in Fig.3a.The original MD simulation curve exhibits a noticeable oscillation range,which is similar to the measured behaviors[9].In Fig.7,the smoothed curve was deter-mined by averaging10–20fore-and-aft points.In Fig.7,there are four characteristic points with force labeled as F1,F2,F3,and F4(the displacements related to F1,F2,F3,F4can be denoted as d1,d2,d3and d4).The eight hollow diamond markers reflect sequentially the forces and displacements at the selected eight time steps in Figs.4and5.As shown in Fig.7,the axial force jumps to a value(F1)then decreases gradually and settles to F2 (about three fourths of F1)at the end of stage one(see Fig.5d).L.Yang et al./Computational Materials Science55(2012)356–364361When entering stage two,the axial force increases steeply to a peak(F3)(see Fig.5f)and then decreases to zero(F4)gradually. The subsequent axial force is small and negligible and maybe due to stress balance along the MWNT’s axial direction.This is sim-ilar to the observation reported in the force–displacement curve [9].Fig.8depicts the influence of the radius of the boundary cylin-drical surface on the axial force–displacement curves with the ra-dius being20Å,23Åand26Å,respectively.The radius of the boundary surface together with the external radius of the MWNT can be used to calculate the area density of carbon nanotube, namely,22.7%,17.2%and13.5%for the selected boundary radii with the three-walled MWNT.For the simulated cases,it is noted that the radius of the boundary cylinder has a notable effect on F4only.It is noted that the axial force–displacement curve at the load-ing end can be described approximately by four straight line seg-ment of with seven characteristic factors(F1,d1),(F2,d2),(F3,d3) and(F4=0,d4).This simplified force-curve with seven parameters can be used as a traction law for characterizing the influence of MWNT on mode I delamination or fracture toughness[5,6].These seven parameters can be determined experimentally by using measured force–displacement curves or numerically by using MD simulations.Hence the next section investigates the influences of the length and external diameter of a MWNT as well as the pulling rate at the loading end on these seven parameters.4.Parametric study and a traction law4.1.Influence of pull rateTo understand the effect of the pull rate at the loading end,MD simulations were performed for a three-walled MWNT at different rates,namely,0.02,0.04,0.06,0.08,0.10,0.12,0.14,0.16,0.18, 0.20Å/fs,respectively.Fig.9depicts the pullout force versus dis-placement curves predicted for ten different pulling rates.It is evident that the pulling rate has a significant effect on the load–displacement curves.Although the topological shape of the curves are similar and can be described by the idealized four segment model with seven parameters,the actual values of all the parameters are quite different.It is noted that the higher the pulling rate,the larger the characteristic forces and displacement.This trend is similar to that identified in previous research,e.g.[21,22],which reported that the rupture force increases with the pulling rate ranging from0.0002 to0.05nm/ps in the MD simulations of pulling individual lipid mol-ecules from a phospholipid membrane in the liquid–crystalline phase.It was also noted that the local viscosity due to local friction is a contributing factor to the difference of the measured rupture forces and the slowest MD simulation results.The MD results show that F2is about two thirds of F1and F4is zero and the ratios d1/d4,d2/d4and d3/d4are approximately0.02,0.5 and0.75respectively.This indicates that the influences of pulling rate on the force–displacement curves can be characterized by its effect on the three parameters:F1,F3and d4.By using curvefitting, these three parameters can be expressed in terms of the non-dimensionalised pulling rates considered.4.2.Influence of MWNT lengthMD simulations were conducted for two-walled MWNTs with different lengths of10nm,28nm,42nm,56nm,72nm.The pulling rate at the loading end is0.2Å/fs.Fig.10depicts the influence of MWNT length on the force–displacement curves.As shown in Fig.10,F1and F2remain nearly the same for all different lengths of the MWNTs,whereas F3increases with the length.Through curvefit-ting,F3and d4can be expressed in terms of the MWNT’s length. 4.3.Influence of MWNT radiusTwo-,four-,five-and six-walled MWNTs of nearly100Ålong were simulated with a pulling rate of0.2Å/fs in the same way as that of the three-walled MWNT model.It should be noted that160180200362L.Yang et al./Computational Materials Science55(2012)356–364the radii boundary for the two-,four-,five-andare 16,26,30and 34Å,respectively.This 23.1%,26.6%and 29.4%MWNT area density.Fig.11depicts the force–displacement four-,five-and six-walled MWNTs.It is noted external radius has a significant influence on ment curves,and in particular F 1and F 3Fig.11also reveals that,except for the d 4is nearly the same for all cases,which means has minimal effect on d 4.Hence only two need to be curved fitted by using the smoothed 4.4.Traction law and mode I fracture toughnessFrom the above parametric study and the equations,we propose the following parameters (F 0=d 0=0),(F 1,d 1),(F 2,d 2),(F 3,d 3)an idealized quadruple-linear model for the epoxy system:F 1¼0:78ÂðR Â0:219À1:11ÞÂ6:56V F 2¼0:52ÂðR Â0:219À1:11ÞÂ6:56VF 3¼0:5ÂðÀ4:2þ0:65R ÞÂð4:1V þ107:3V 2À291:7V 3ÞÂð0:83þ0:00052L Þð3Þwhere the unit of length and radius of the MWNT is Å.As the ratios d 1/d 4,d 2/d 4,d 3/d 4can approximated as 0.02,0.5,0.75,respectively,and the d 4is only proportional to the loading rate and the MWNT length.The displacements of the MWNT,d 1,d 2,d 3and d 4,are given byd 1¼0:17VL d 2¼4:31VL d 3¼6:47VL d 4¼8:62VL ð4ÞThe area beneath the force–displacement curve represents the additional fracture toughness offered by a MWNT.When the MWNTs are distributed evenly along the length of the saturated bridging zone,the additional mode I fracture toughness D G con-tributed by MWNTs across the width of the DCB specimen can be determined by [7]:D G ¼2SZdF ðd Þd d ¼2SE MWNT%2ÂS Âf 2:8VL ½ðR Â0:219À1:11ÞÂ6:56Vþ1:08VL ðÀ4:2þ0:65R ÞÂð4:1V þ107:3V 2À291:7V 3ÞÂð0:83þ0:00052L Þgð5Þwhere S is the number of MWNTs per unit length.The above equa-tion is important as it relates the MD simulation determined trac-tion law to the increment of the mode I fracture toughness in a standard DCB specimen.The area beneath the full traction law curve plays an important role in capturing the contribution of bridging MWNT between two cracked parison with experimental resultIn this section,we compare the proposed traction law,i.e.,four-line force–displacement curve,with the experimental results re-ported in [9].Cooper et al.[9]conducted an experiment of detaching nanotubes from an epoxy by using a SPM tip and reported the mea-sured force-distance curve for a MWNT of 8.2nm in diameter and 484nm in length.The maximum diameter and length of the MWNT in the MD simulation is 4.16nm and 72nm,respectively,so Eqs.(1)and (2)were used in this prediction.From Fig.3in Ref.[9],it is noted that d 4measures about 60nm (440nm to 500nm)between the mo-ment of the SPM tip coming in contact with the CNT and that of CNT being pulled out.As for Eq.(4),the MD simulation velocity is chosen as 0.015Å/fs corresponding to the real experimental speed.The SPM tip speed was not given in Ref.[9],but is believed to be on the order of magnitude of l m/s.As the MD simulation velocity is a virtual velocity and its corresponding real velocity could be influenced by different factors,e.g.,such as test setup,crack propagation speed and data acquisition method.If we assume the SPM tip speed in [9]is 1l m/s,Eqs.(3)and (4)can be useful in MWNT reinforced com-posite delamination analysis.From Eqs.(3)and (4),we can deter-mine that F 1=0.6l N,F 2=0.4l N and F 3=3.2l N;and d 1=1.2nm,d 2=31.3nm,d 3=47.0nm and d 4=62.6nm.Fig.12compares the predicted idealized force–displacement curve with the measured one from Cooper et al.[9]in the range of 440nm to 530nm in which the friction force with Kapton film (about 2.6l N)was subtracted.A broad agreement exists between the two curves although there are differences in actual values.One reason could be due to difference in epoxy systems,i.e.,EPON 862epoxy resin reacted with DETDA cur-ing agent was used in the MD model and Araldite LY564with hard-ener HY560was used in the experiment [9].6.ConclusionThe present simulation reveals that the MWNT extraction con-sists of two stages.In the simulation,the MWNT’s progressive stretching,necking and debonding phenomena have been noted prior to the movement of the end embedded in resin.The pulling force versus displacement curve is characterized by four distinct pairs of forces and displacements.The influence of the radius andparisons of the pull force–displacement curves between present prediction and the experimental measurement [9].L.Yanglength of MWNT and pulling rate were studied and then curved-fit-ted to obtain two equations describing the proposed quadruple-linear traction law.The proposed quadruple-linear traction law is then validated by comparing favorably with existing experimental result.By using the MD based quad-linear traction law,a simple formula was developed to estimate the increase in mode I energy release rate due to presence of MWNTs.AcknowledgementsThe authors are grateful to thefinancial support of the China Scholarship Council to Lin Yang during his study at the University of Sydney.This work wasfinancially supported by the National Natural Science Foundation of China(NSFC,Grant2009-90916008,2011-11002042)and Key Laboratory Opening Funding of National Key Laboratory for Advanced Composites in Special Environment(HIT.KLOF.2009038).The authors also wish to thank Lei Mei of HIT for her work in preparing carbonfiber with con-nected carbon nanotubes.Appendix AE with-bond can be further decomposed intofive items(E bond,E angle, E hideral,E improper,E coupling)and they are given by[7]:E with-bond¼E bondþE angleþE dihedralþE improperþE couplingðA:1ÞE bond¼Xb½k2ðbÀb0Þ2þk3ðbÀb0Þ3þk4ðbÀb0Þ4 ðA:2ÞE angle¼Xh½k2ðhÀh0Þ2þk3ðhÀh0Þ3þk4ðhÀh0Þ4 ðA:3ÞE dihedral¼X/½k1ð1Àcos/Þ2þk2ð1Àcos2/Þ3þk3ð1Àcos3/Þ4 ðA:4ÞE improper¼Xvk v v2ðA:5ÞE coupling¼Xb Xb0kbb0ðbÀb0Þðb0Àb0ÞþXbXhk b hðbÀb0ÞðhÀh0ÞþXb X/k b/ðbÀb0ÞðF1b/ðcos/ÞþF2b/ðcos2/ÞþF3b/ðcos3/ÞþXh X/k b/ðhÀh0ÞðF1h/ðcos/ÞþF2h/ðcos2/ÞþF3h/ðcos3/ÞþXh Xh0k hh0ðhÀh0Þðh0Àh0ÞþXhX/XhðhÀh0ÞÂðh0Àh00Þcos/ðA:6Þwhere b is the bond length,h is the bond angle,u is the dihedral an-gle,v is the out-of-plane angle;b0,h0,u0and v0are the equilibriumbond length,bond angle,dihedral angle,and out-of-plane angle;k iand k a b are the quadratic,cubic and quartic forcefield and couplingparameters.The E non-bond is composed of two items,the Coulumbic electro-static interaction energy E elec and the var der Waals energy E vdw,E vdw¼Xe2dr ij9À3dr ij6"#for r<r c1ðA:7ÞE coul¼Xi;jqiqjr ijfor r<r c2ðA:8Þwhere e is the negative of minimum van der Waals energy,d is thedistance parameter.The r ij is the distance between atoms i and j,r c1and r c2are the cut-off distance for var der Waals energy E vdW andColumbic electrostatic interaction energy E elec,and q i and q j arethe charges on the two atoms.References[1]L.Tong,A.P.Mouritz,M.K.Bannister,3D Fiber Reinforced Polymer Composites,Elsevier Science Ltd.,Oxford,UK,2002.[2]V.P.Veedu,A.Cao,X.Li,et al.,Nat.Mater.5(2006)457–462.[3]X.D.He,F.H.Zhang,R.Wang,et al.,Carbon45(13)(2007)2559–2563.[4]Mei.L,X.D.He,Y.Li,et al.,Mater.Lett.64(2010)2505–2508.[5]E.J.Garcia,B.L.Wardle,A.J.Hart,Compos.Part A–Appl.Sci.Manuf.39(2008)1065–1070.[6]L.Tong,X.Sun,P.Tan,pos.Mater.42(1)(2008)5–23.[7]J.Blanco,E.J.Garcia,R.G.De Villoria,et al.,pos.Mater.43(8)(2009)825–841.[8]F.H.Zhang,Study on Carbon Nanotube/Carbon Fiber Multi-scale Reinforcementand Its Composite Interfacial Properties,PhD Dissertation,Harbin Institute ofTechnology,2008,pp.60–61.[9]C.A.Cooper,S.R.Cohen,A.H.Barber,et al.,Appl.Phys.Lett.81(2002)3873–3875.[10]D.Qian,E.C.Dickey,R.Andrews,et al.,App.Phys.Lett.76(20)(2000)2868–2870.[11]C.Bower,R.Rosen,L.Jin,et al.,Appl.Phys.Lett.74(1999)3317–3319.[12]A.H.Barber,S.R.Cohen,H.D.Wanger,Appl.Phys.Lett.82(2003)4140–4142.[13]K.Liao,S.Li,Appl.Phys.Lett.79(2001)4225–4227.[14]S.J.V.Frankland,V.M.Harik,Surf.Sci.525(2003)103–108.[15]J.Gou,B.Minaie,B.Wang,et al.,Comput.Mater.Sci.31(2004)325–336.[16]H.Sun,J.Phys.Chem.B102(1998)7338–7364.[17]R.J.Meier,J.R.Maple,M.-J.Hwang,et al.,J.Phys.Chem.99(5)(1995)5445–5456.[18]F.Deng,T.Ogasawara,N.Takeda,Compos.Sci.Technol.67(2007)2959–2964.[19]M.F Yu,O.Lourie,M.J.Dyer,et al.,Science287(2000)637–640.[20]J.Y Huang,S.Chen,Z.Q Wang,et al.,Nature439(2006)281.[21]S.-J.Marrink,O.Berger,P.Tieleman,et al.,Biophys.J.74(1998)931–943.[22]A.C.Lorenzo,P.G.Pascutti,P.M.Bisch,put.Chem.24(3)(2003)328–339.[23]G.L.Hwang,K.C.Hwang,Nano Lett.1(2001)435–438.[24]Z.L.Wang,P.Poncharal,W.A.de Heer,J.Phys.Chem.Solids61(2000)1025–1030.[25]W.Ding,A.Eitan,F.T.Fisher,X.Chen,D.A.Dikin,R.Andrews,L.C.Brinson,L.S.Schadler,R.S.Ruoff,Nano Lett.3(2003)1593–1597.[26]A.F.Fonseca,T.Borders,R.H.Baughman,K.J.Cho,Kyeongjae,Phys.Rev.B81(4)(2010).Article Number:045429.364L.Yang et al./Computational Materials Science55(2012)356–364。