碳纳米管(carbon nanotube,CNT)和石墨烯(graphene,Gr)因其特殊的结构和优异的机械、光学及电学性能而受到广泛关注，被认为是未来的革命性材料[1-5]。 Gr 是一种由碳原子以sp2 杂化轨道组成的六边形蜂窝晶格结构的二维纳米材料[6]。 CNT 的基本结构与Gr 相似，是由Gr 围绕中心轴按一定螺旋角卷曲而成的一维纳米材料[3]。 近年来，CNT 和Gr 两种纳米碳材料作为增强材料被广泛应用于金属基复合材料(metal matrix composites, MMCs)的制备[7-9]。特别是在铝基复合材料(aluminium matrix composites,AMCs)中，与传统陶瓷颗粒相比，CNT 和Gr 具有更优越的强化效果，可以在很大程度上提高材料的力学性能[10-11]。
大量实验和计算研究表明，CNT 和Gr 都能显著改善AMCs 的力学性能。 然而由于实验制备技术和纳米表征技术的限制，碳纳米材料在金属基体中的增强机制还未被深刻理解。 碳纳米材料的分布和制备方法等对复合材料力学性能有很大影响。 想通过实验来研究还需要较长的时间， 分子动力学(molecular dynamics,MD)模拟方法是研究复杂纳米材料系统最有效的方法之一[12-14]。 近年来，有大量报道利用MD模拟方法研究CNT 和Gr 增强金属复合材料的力学性能。 Silvestre 等[15]采用MD 方法研究了CNT/Al的压缩性能，与纯Al 相比，杨氏模量增加约50%。而屈服应力和屈服应变都没有增加，这是由于CNT的局部屈曲导致碳纳米管-铝复合材料过早失效。Zhu 等[16]采用MD 方法对纳米层状结构的Gr/Al 进行了压缩实验，研究发现Gr 会限制位错的滑移，导致界面上位错密度的增加。 Lei 等[17]研究发现Gr 涂层会使Al 基体的承载面积增加约5.36 倍， 并改变了Al 基体在纳米压痕过程中的变形行为。Peng 等[18]研究了具有Gr 涂层的铜基体纳米压痕塑性变形，与纯铜相比，单、双和3 层Gr 涂层的临界载荷分别提高了114%、204%和344%。 但是目前对于结构一致性的CNT 和Gr 对金属基复合材料性能影响的对比研究还未见报道。
本文采用分子动力学方法对纯铝、CNT/Al 和Gr/Al 进行了相同条件下的纳米压痕模拟，得到了材料的载荷-压痕深度曲线， 分析了CNT 和Gr 对铝基复合材料压痕性能的提升作用。 通过剪切应变分布、 原子位移矢量和压痕形貌3 个方面， 阐明了CNT 和Gr 对铝基复合材料纳米压痕塑性变形行为的影响， 最后结合材料中的位错演化过程揭示了位错与CNT 和Gr 的相互作用。
1 模拟方法
在本研究中，采用LAMMPS[19]软件进行纳米压痕模拟。 原子模型由ATOMSK[20]和LAMMPS 软件共同构建，模型如图1 所示。纯铝的纳米压痕模型如图1a 所示，CNT/Al 和Gr/Al 复合材料模型如图1b 和c 所示。Al 基体晶体取向为x-[112]，y-[110]，z-[111]，基体整体尺寸为158.7 Å×147.3 Å×150.3 Å，共包含213 696 个原子。 为了保证纳米压痕模拟过程中模型结构的稳定性，基体由3 部分组成，分别是固定层、恒温层和牛顿层[21]。固定层和恒温层厚度分别为7 和10 Å。 压头为金刚石压头， 半径为30 Å，距Al 基体表面15 Å。 CNT 和Gr 在Al 基体内部同一高度位置处，距Al 基体上表面30 Å，图1b 和c 中的CNT 和Gr 结构分别如图1d 和e 所示。为了保证模型所用CNT 和Gr 具有可对比性，两者手性、方向及原子数均相同。
图1 纳米压痕模拟原子模型：(a)Al;(b)CNT/Al;(c)Gr/Al;(d)CNT;(e)Gr
Fig.1 Atomic models used for nanoindentation simulations：(a)Al;(b)CNT/Al;(c)Gr/Al;(d)CNT;(e)Gr
原子间相互作用势函数的准确性是决定分子动力学模拟结果准确性的重要因素之一。在本研究中，金属Al 基体中Al-Al 原子间势函数为EAM[22-23]，CNT和Gr 结构中的C-C 原子间势函数为AIREBO[24]，CNT 和Gr 与Al 基体之间的Al-C 原子间势函数为Morse[25]。金刚石压头与Al 基体之间的C-Al 原子间作用势函数为LJ[26-27]。 金刚石压头原子间相互作用力被忽略， 这是由于金刚石比金属基体硬得多，这样处理可以节约计算资源。
在进行纳米压痕模拟时， 体系参数设置如下：x和y 方向上为周期性边界条件，z 方向为固定边界条件。 固定边界层原子， 恒温层原子温度预设为300 K， 牛顿层速度不预设， 控温方式为速度标定法，并不断进行温度修正。 步长为0.001 ps，系综为NVE，弛豫时间为500 ps。 在实现纳米压痕过程时，压头下降速度为0.1 Å/ps，压入最大深度为30 Å，利用压头质心受力计算载荷，保载时间为20 ps，随后压头上升回复至原位。 模拟结果采用了OVITO[28]软件实现可视化，利用原子应变[29]分析剪切应变，共近邻分析[30]识别原子晶体结构，以及位错提取算法[31]识别位错线和晶体缺陷。
2 实验结果及讨论
2.1 纳米压痕载荷-压痕深度曲线
纯Al、CNT/Al 和Gr/Al 复合材料纳米压痕模拟加载过程中载荷与压痕深度关系如图2 所示。 加载曲线分为2 部分， 第1 部分为弹性变形阶段，根据赫兹弹性理论[32]，球面接触的弹性响应载荷F 和压痕深度h 的关系为：
图2 Al、CNT/Al 和Gr/Al 复合材料的载荷-压痕深度曲线
Fig.2 Load-indentation depth curves of the pure Al,CNT/Al and Gr/Al composites
式中，E*为约化弹性模量；压头半径R=3 nm。E*可以通过下式计算得到：
式中，Es 为弹性模量；νs 为泊松比， 纯铝的泊松比νs=0.33。 计算得到Gr/Al、CNT/Al 和Al 的弹性模量分别为115.92、109.82 和97.61 GPa，Gr/Al 和CNT/Al相比与纯Al 分别提高了18.76%和12.51%。
随着压痕深度的增大，材料进入塑性变形阶段，通过载荷曲线可以看出材料的抗压能力。 在压痕深度为5.6~7.8 Å 时，Gr/Al 相比于CNT/Al 和纯Al 抗压强度有明显上升。 随后3 种材料的压痕载荷相差不大， 通过计算塑性变形阶段的载荷平均值，Gr/Al、CNT/Al 和Al 分别为0.1876、0.1795 和0.186 4 μN。由此可见，Gr 对Al 基体的抗压能力有提升作用，而CNT 对Al 基体的抗压能力反而有弱化作用。 这主要是由于Gr 的平面状结构能够在面内有效承载压力，使Gr/Al 复合材料的纳米压痕强度提升；而CNT的中空管状结构在径向抗压能力较差， 且在压头压入方向的垂直面上有效承载面积较小， 故而导致CNT/Al 复合材料的纳米压痕强度较低。
2.2 纳米压痕塑性变形行为
为了进一步分析CNT 和Gr 在铝基复合材料中的纳米压痕塑性变形行为，提取了纯Al(图3a1~a5)、CNT/Al(图3b1~b5)和Gr/Al(图3c1~c5)在不同压痕深度时原子剪切应变的截面快照，如图3 所示。在压痕深度为5 Å 时(图3a1~c1)，材料开始发生塑性变形，基体中出现滑移带，此时抗压强度会出现小幅下降。 在纯铝中(111)滑移面被激活，压头下发生塑性变形的原子呈倒梯状分布； 在CNT/Al 中同样是(111) 滑移面被激活， 但压头下塑性变形原子呈与(111)面平行向下扩展；在Gr/Al 中塑性变形原子沿(111)面平行向下扩展过程中再次出现(111)滑移面，最终在塑性变形区前端汇聚成尖凸状， 增大了进一步塑性变形的难度，从而提高材料的整体抗压强度。随着压痕深度的增大，当d=10 Å 时(图3a2~c2)，可以明显观察到3 种材料塑性变形扩展行为的差异。在纯铝中，塑性变形区以台阶状向下扩展，高应变原子所在位置即发生下一次变形处；在CNT/Al 中，由于CNT 的界面引导作用，导致滑移区并没有像纯铝以台阶状连续向下扩展，而且直接沿(111)滑移面原子扩展至CNT-Al 界面处，在CNT 上方包围出一个倒梯形低塑性变形区。 这种弧形界面的引导作用使得CNT/Al 复合材料的抗压强度在此处发生较大下降，这也是图2 中蓝色实线在10 Å 处迅速下降的原因； 在Gr/Al 中，Gr 的平面状结构将塑性变形区引导至呈现出与Gr 宽度一致的矩形区，Gr 的存在能够有效阻碍高塑性变形区向Gr 下方基体的应变传递，从而提升Gr/Al 复合材料的抗压强度。在压痕深度为20 和30 Å 时，纯Al 仍然按照台阶式进行塑性变形区扩展；CNT 在压头的持续下降过程中被压扁，CNT/Al 基体中变形区沿(111)滑移面进行多条滑移带的扩展；Gr 持续阻碍高应变向下传递， 基体高塑性变形始终在Gr 所在平面以上进行， 原子大量向上运动并堆积在表面处。 图3(a5~c5)是提取了材料在压痕深度为30 Å 时剪切应变值大于0.1 的原子区域， 通过对比3 种材料的高应变区深度，纯Al、CNT/Al 和Gr/Al 分别为109、62 和40 Å，可以直观得看出CNT 和Gr 对阻碍压痕在金属基体中产生更大塑性变形的有益效果。
图3 在5、10、20 和30 Å 压痕深度下原子剪切应变截面快照：(a1~a5)Al;(b1~b5)CNT/Al;(c1~c5)Gr/Al
Fig.3 Snapshots of atomic shear strain cross sections at indentation depths of 5,10,20 and 30 Å：(a1~a5)Al;(b1~b5)CNT/Al;(c1~c5)Gr/Al
图4为纯铝、CNT/Al 和Gr/Al 在不同压痕深度时原子位移矢量的截面快照，图中箭头所指方向为原子位移矢量方向。3 种材料在压痕作用力下，原子运动方向的规律表现出明显的差异。 在压痕深度为10 Å 时，纯铝(图4a1)的原子位移矢量方向整体为沿Z 轴负方向，原子位移矢量大小以压头与基体接触点为中心呈放射状渐变分布；CNT/Al 复合材料(图4b1)由于CNT 的存在使上方基体原子的运动方向汇聚于CNT 的上表面， 原子位移矢量方向呈锥形分布，与图3b2 中的高剪切应变分布一致。 这种CNT 的引导作用会使CNT 的上表面受较大应力分布， 导致CNT 上表面的凹陷；Gr/Al 复合材料(图4c1)由于Gr 的存在，使得上表面基体原子均以垂直向下方向进行原子运动，Gr 表面的载荷分布均匀，且片状结构承载力高，能够进一步阻碍原子位移的向下传递。 随着压痕深度进一步增大到20 Å 时，纯铝(图4a2)的原子位移矢量方向以被激活滑移系方向为主导， 整体仍呈放射状分布；CNT/Al 复合材料(图4b)的原子位移矢量方向一旦突破CNT 上表面的承载，就和纯铝的位移矢量方向一致，呈放射状并以滑移系方向为主导；Gr/Al 复合材料(图4c2)由于Gr的阻碍作用更强，原子位移无法通过Gr 界面进行有效传递，因此出现原子被向上反推现象。 当压痕深度到30 Å 时，材料中都出现原子被反推向上运动，纯铝以多个滑移系方向为主要运动方向进行扩展(图4a3)，CNT/Al 复合材料则以CNT 下部垂直方向为主要运动方向进行扩展(图4b3)，Gr/Al 复合材料中Gr 的阻碍作用阻止了大部分原子的向下运动，主要承力点仍在Gr 上，大部分基体原子被反推向上进行扩展(图4b3)。
图4 在10、20 和30 Å 压痕深度下原子位移矢量截面快照：(a1~a3)Al;(b1~b3)CNT/Al;(c1~c3)Gr/Al
Fig.4 Snapshot of atomic displacement vector cross sections at 10,20 and 30 Å indentation depths：(a1~a3)Al;(b1~b3)CNT/Al;(c1~c3)Gr/Al
图5a 和c 展示了压痕深度为30 Å 时CNT 和Gr 的俯视和正视原子形貌以及原子剪切应变。CNT在受压痕载荷时中间段发生凹陷，而Gr 仍然保持平面状， 没有由于压痕载荷产生大变形。 同时， 观察CNT 和Gr 的原子剪切应变分布，以Gr 受压中心点沿x 方向两侧分布的剪切应变值最高，说明Gr 在受压时中心线两侧原子提供了较大的支撑力， 避免了Gr 产生较大的中心凹陷； 而CNT 中能起到支撑的原子数少，导致CNT 产生了较大的中心凹陷。 值得强调的是，原始模型中CNT 在z 方向上的最高点和Gr 相同，都是距离基体表面30 Å。 当压痕深度为30 Å 时，将模型垂直于x 方向进行截面观察，可以发现在CNT/Al 复合材料模型(图5b)中压头底部距离被压的CNT 顶部10 Å， 而Gr/Al 复合材料模型(图5d)中压头底部距离被压的Gr 顶部3 Å。图5b 和d中蓝色实线框出的是压痕深度为30 Å 时，CNT/Al 和Gr/Al 基体受压向上运动后形成的表面堆积原子。CNT/Al 的表面堆积原子数和堆积高度均小于Gr/Al， 受压变形的CNT 为基体原子的移动提供了空间，降低了原子运动所需能量，导致材料抗压性能下降。
图5 压痕深度为30 Å 时原子形貌和剪切应变：(a)CNT 的俯视与正视原子形貌及剪切应变；(b)CNT/Al 的原子截面形貌；(c)Gr的俯视与正视原子形貌及剪切应变；(d)Gr/Al 的原子截面形貌
Fig.5 Atomic morphology and atomic shear strain at an indentation depth of 30 Å：(a)top-view and front-view atomic morphology and shear strain of CNT;(b)atomic cross-sectional morphology of CNT/Al;(c)top-view and front-view atomic morphology and shear strain of Gr;(d)atomic cross-sectional morphology of Gr/Al
通过统计CNT 中间段压缩量与压痕深度曲线，如图6 所示，CNT 凹陷过程可以分为4 个阶段，分别是平稳阶段、缓慢波动上升阶段、线性快速上升阶段和稳定阶段。 当压痕深度小于6 Å 时，CNT的压缩量变化较小，属于正常原子运动范围，还未发生凹陷现象。 当压痕深度为6~16 Å 时，此时压缩量开始缓慢波动增大，CNT 出现被压缩现象，CNT表面发生一定程度的弯曲。 当压痕深度为16~24 Å时，CNT 的压缩量与压痕深度呈线性关系并快速增加，CNT 发生较大压缩变形，几乎失去变形抗力。 直到压痕深度大于24 Å 时，CNT 中段完全被压扁，此时CNT 中间段压缩量不再变化， 进一步变形以CNT 中间段向两侧传递方式进行。
图6 CNT 中间段压缩量与压痕深度曲线
Fig.6 Compression amount of the middle part of the CNT-indentation depth curve
2.3 纳米压痕位错演化行为
为了进一步研究材料在纳米压痕过程中的滑移行为及位错与CNT 和Gr 的相互作用， 利用共近邻原子分析(common neighbor analysis,CNA)和位错分析(dislocation analysis, DXA)获得不同压痕深度下材料的位错与缺陷结构信息，如图7 所示。 首先，在压痕深度为5 Å 时，材料开始发生塑性变形，纯Al(图7a1)中的位错以棱柱位错环形核，所激活的滑移面数量较多；CNT/Al(图7b1)中的位错以单位错线被激活，滑移面为较难激活的(111)面；Gr/Al(图7c1)中的位错同样以单位错线被激活， 滑移面与纯铝中类似，没有出现(111)滑移面。 当压痕深度为10 Å 时，纯Al(图7a2)中不仅在压头正下方出现了位错，同时在未被压头压入的位置也出现了位错环。 这表明纯铝在受压痕载荷作用时，基体内部的塑性变形传播速度更快，传播范围更广；CNT/Al(图7b2)中由于CNT的阻碍作用， 位错从压头-Al 界面连接到CNT-Al界面，形成压头下方的刃型位错；Gr/Al(图7c2)中，同样由于Gr 的阻碍作用， 位错主要分布在两个界面之间，此时CNT 和Gr 与位错相互作用的差异性还未体现；随着压痕深度增大到20 Å 时，纯Al(图7a3)中已出现大量缠结的位错环，位错滑移难度低，滑移距离大，基体内部发生较大塑性变形；CNT/Al(图7b2)中由于CNT 的管状结构上表面呈弧形，且承载面积较小，对位错的阻碍作用有限，当位错突破CNT 的阻碍后， 出现与纯Al 中相似的位错环并向下滑移；Gr/Al(图7c3)中Gr 对于位错的阻碍作用显著高于CNT，且形成的位错中肖克莱位错数量较少，全位错和梯杆位错数量提升，不可动位错比例的提高能够进一步提高材料的压痕抗力。 除此之外，Gr/Al 中的位错分布也较为集中，大多位于压头底部，既无法向下滑移，也不向周围扩展；当压痕深度到30 Å 时，纯Al(图7a4)中已经出现脱落的位错环，压头正下方位错较少，沿滑移系方向大量滑动；CNT/Al(图7b4)中虽然在CNT 下方出现位错环， 并有向下滑动趋势，但由于CNT 下界面的存在， 对位错的钉扎作用明显，未出现脱落的位错环；Gr/Al(图7c4)中的位错分布情况与压痕深度在20 Å 时类似，位错并没有随着压痕深度的增加而突破Gr 的阻碍。
图7 在5、10、20 和30 Å 压痕深度下位错及缺陷结构代表性快照：(a1~a4)Al;(b1~b4)CNT/Al;(c1~c4)Gr/Al(原子根据计算的CNA 着色，图中只显示了HCP 原子，其他原子未显示)
Fig.7 Representative snapshots of dislocations and defect structures at indentation depths of 5,10,20 and 30 Å：(a1~a4)Al;(b1~b4)CNT/Al;(c1~c4)Gr/Al(atoms are colored according to the calculated CNA,only HCP atoms are shown in the figure;other atoms are not shown)
3 结论
(1)载荷-压痕深度曲线中塑性变形段Gr/Al、CNT/Al 和Al 的平均载荷值分别为0.187 6、0.179 5 和0.186 4 μN。 Gr 对Al 基体的抗压能力有提升作用，而CNT 对Al 基体的抗压能力反而有弱化作用。
(2)压痕过程中Gr/Al 复合材料的塑性变形区域被限制在Gr 以上部分，变形区明显小于CNT/Al 和纯Al，且表面原子堆积数要明显高于CNT/Al，具有更高的抗压强度和较低的塑性变形能力。CNT 在压痕深度大于16 Å 后发生较大压缩量变形，这一过程是导致材料抗压能力较弱的主要原因。
(3)CNT 对位错的阻碍作用有限，位错一旦越过CNT 基体中会出现与纯Al 基体中类似的位错环，但由于CNT-Al 下界面的存在， 对位错还有一定的钉扎作用， 使得位错脱落的难度增大；Gr 对位错的阻碍作用明显优于CNT， 位错被阻挡在Gr 以上部分难以滑过， 且Gr/Al 中的不可动位错比例随压痕深入不断增大，能够进一步提升位错强化作用。
[1] GEIM A K,NOVOSELOV K S.The rise of graphene[J].Nature Materials,2007,6：183-191.
[2] PARK S, RUOFF R S. Chemical methods for the production of graphenes[J].Nature Nanotechnology,2009,4：217-224.
[3] IIJIMAS.Helical microtubulesof graphitic carbon[J].Nature,1991,354：56-58.
[4] ESAWI A M K,FARAG M M. Carbon nanotube reinforced composites： Potential and current challenges[J]. Materials & Design,2007,28(9)：2394-2401.
[5] POPOV V N. Carbon nanotubes： Properties and application[J].MaterialsScienceandEngineering：R：Reports,2004,43(3)：61-102.
[6] YAZYEV O V, CHEN Y P. Polycrystalline graphene and other two-dimensional materials[J]. Nature Nanotechnology, 2014, 9：755-767.
[7] CHEN B, SHEN J, YE X, JIA L, LI S, UMEDA J, TAKAHASHI M, KONDOH K. Length effect of carbon nanotubes on the strengthening mechanisms in metal matrix composites[J].Acta Materialia,2017,140：317-325.
[8] CHEN B,KONDOH K,LI J S,QIAN M.Extraordinary reinforcing effect of carbon nanotubes in aluminium matrix composites assisted by in-situ alumina nanoparticles[J]. Composites Part B： Engineering,2020,183：107691.
[9] ZHANG X, ZHAO N Q, HE C N. The superior mechanical and physical properties of nanocarbon reinforced bulk composites achieved by architecture design - A review[J].Progress in Materials Science,2020,113：100672.
[10] KINLOCH I A, SUHR J, LOU J, YOUNG R J, AJAYAN P M.Composites with carbon nanotubes and graphene： An outlook[J].Science,2018,362(6414)：547-553.
[11] ZHAO Z Y,BAI P K,DU W B,LIU B,PAN D,DAS R,LIU C T,GUO Z H.An overview of graphene and its derivatives reinforced metal matrix composites： Preparation, properties and applications[J].Carbon,2020,170：302-326.
[12] YOO S C,LEE D J,RYU S W,KANG B,RYU H J,HONG S H.Recent progress in low-dimensional nanomaterials filled multifunctional metal matrix nanocomposites[J]. Progress in Materials Science,2023,132：101034.
[13] YANG Y J,LIU M B,DU J X,ZHANG W,ZHOU S Q,REN W J,ZHOU Q H,SHI L.Construction of graphene network in Ni matrix composites：A molecular dynamics study of densification process[J].Carbon,2022,191：55-66.
[14] CHEN S, AITKEN Z H, PATTAMATTA S, WU Z X, YU Z G,SROLOVITZ D J, LIAW P K, ZHANG Y W. Simultaneously enhancing the ultimate strength and ductility of high-entropy alloys via short-range ordering[J].Nature Communications,2021,12：4953.
[15] SILVESTRE N,FARIA B,CANONGIA LOPES J N.Compressive behavior of CNT-reinforced aluminum composites using molecular dynamics[J]. Composites Science and Technology, 2014, 90：16-24.
[16] ZHU J Q,LIU X,YANG Q S.Dislocation-blocking mechanism for the strengthening and toughening of laminated graphene/Al composites[J].Computational Materials Science,2019,160：72-81.
[17] LEI Y X,YAN Y P,LYU J J.Atomistic study of the strengthening mechanisms of graphene coated aluminum [J]. Nanotechnology,2020,31(5)：055708.
[18] PENG W X, SUN K, ZHANG M, SHI J Q, CHEN J. Effects of graphene coating on the plastic deformation of single crystal copper nanocuboid under different nanoindentation modes[J].Materials Chemistry and Physics,2019,225：1-7.
[19] THOMPSON A P,AKTULGA H M,BERGER R,BOLINTINEANU DS,BROWNWM,CROZIERPS,IN'TVELDPJ,KOHLMEYERA,MOORE S G,NGUYEN T D,SHAN R,STEVENS M J,TRANCHIDAJ,TROTT C,PLIMPTON S J.LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic,meso,and continuum scales[J]. Computer Physics Communications, 2022,271：108171.
[20] HIREL P.Atomsk：A tool for manipulating and converting atomic data files[J].Computer Physics Communications, 2015,197： 212-219.
[21] HUA D P,XIA Q S,WANG W,ZHOU Q,LI S,QIAN D,SHI J Q,WANG H.Atomistic insights into the deformation mechanism of a CoCrNi medium entropy alloy under nanoindentation[J].International Journal of Plasticity,2021,142：102997.
[22] LIU C L, PLIMPTON S J. Molecular dynamics simulations of grain boundary diffusion in Al using embedded atom method potentials[J].Journal of Materials Research,1995,10：1589-1592.
[23] MISHIN Y, FARKAS D, MEHL M J, PAPACONSTANTOPOULOS D A. Interatomic potentials for monoatomic metals from experimental data and ab initio calculations[J].Physical Review B,1999,59(5)：3393-3407.
[24] STUART S J,TUTEIN A B,HARRISON J A.A reactive potential for hydrocarbons with intermolecular interactions[J].The Journal of Chemical Physics,2000,112(14)：6472-6486.
[25] DANDEKAR C R,SHIN Y C.Molecular dynamics based cohesive zone law for describing Al-SiC interface mechanics[J].Composites Part A：Applied Science and Manufacturing,2011,42(4)：355-363.
[26] PIAιTEK A, NOWAK R, GBURSKI Z. A titanium-decorated fullerene cluster - a molecular dynamics simulation[J].Solid State Phenomena,2008,140：109-116.
[27] HALICIOGˇLU T, POUND G M. Calculation of potential energy parameters form crystalline state properties[J].Physica Status Solidi(a),1975,30(2)：619-623.
[28] STUKOWSKI A. Visualization and analysis of atomistic simulation data with OVITO - the open visualization tool[J].Modelling and Simulation in Materials Science and Engineering,2010,18(1)：015012.
[29] SHIMIZU F,OGATA S,LI J.Theory of shear banding in metallic glasses and molecular dynamics calculations[J].Materials Transactions,2007,48(11)：2923-2927.
[30] HONEYCUTT J D, ANDERSEN H C. Molecular dynamics study ofmeltingandfreezingofsmallLennard-Jonesclusters[J].The Journal of Physical Chemistry,1987,91(19)：4950-4963.
[31] STUKOWSKI A, BULATOV V V, ARSENLIS A. Automated identification and indexing of dislocations in crystal interfaces[J].Modelling and Simulation in Materials Science and Engineering,2012,20(8)：085007.
[32] HEINRICH H. Ueber die Berührung fester elastischer körper[M].CRELLE A L, BORCHARDT C W, et al. Band 92. Berlin： De Gruyter,1882：156-171.