纳米管拉伸的分子动力学模拟

纳米管拉伸的分子动力学模拟
孙欢
【摘要】致力于对纳米单晶铝杆的拉伸过程进行仿真模拟研究.应用了分子动力学方法,所用势函数为紧束缚TB势函数,时间积分为Verlet速度算法.通过矢量叠加,得到系统的总能量及受力情况,进而通过系统达平衡后进行位移加载的方法,定性分析纳米尺度金属的拉伸特性,通过精确度分析,最终得出其相应的应力应变曲线与形状在宏观上大体一致的结论.
【期刊名称】《甘肃科学学报》
【年(卷),期】2013(025)002
【总页数】3页(P25-27)
【关键词】铝纳米杆;分子动力学;Verlet速度算法;应力一应变
【作者】孙欢
【作者单位】陕西师范大学数学与信息科学学院,陕西西安710062
【正文语种】中文
【中图分类】O795
纳米技术是在纳米尺度［1，2］0．1～100nm 之间上研究物质（原子、分子）的特性和相互作用的技术．某些普通材料成为纳米量级后就会具有防垢、韧性好、耐高温、耐腐蚀等特殊的功能．这就可能改变相关产品的设计和制造方式，实现生产方式的飞跃［3］．
分子动力学方法的模拟结果是否精确很大一部分在于原子间相互作用势函数的选取．我们所选用的势函数确定为紧束缚TB势并假设满足以下条件［4］：
（1）粒子的运动与牛顿力学运动定律一致；
（2）粒子之间的相互作用满足叠加原理．
以上条件没有考虑量子效应与多体作用的因素，并不完全符合真实的物理系统，所以这是一种近似计算．
假设N为系统的原子数，第i个原子的质量是mi，位置坐标向量为ri，速度为vi ＝，加速度为受到的作用力为Fi，原子i与原子j之间的距离为原子j对原子i的作用力为fij，原子i与原子j相互作用势能为φ（rij）．系统的所有的物理量都是随时间变化的，即A＝A（t），控制方程为
由此建立一个线性的微分方程组，给定初始位置和速度，可以得到任意时刻系统中所有原子的位置ri（t）和速度vi（t），进而可得到系统的总动能［5］，由于势函数已知，故可得到系统的总能量．进行位移加载可导致各原子受力的变化，也可起到拉伸的效果．施加应力求应变与施加位移求应力的道理是一样的．
1 原子模型分析
1．1 模型的建立
分子动力学模拟分为6个步骤：
（1）势函数的选取纳米铝杆初始结构模型按几何方法生成，原子按理想面心结构分布，X、Y坐标轴分别对应面心立方的［001］、［001］晶向．尺寸为
6a0×16a0，其中a0 为铝的晶格常数（4．050nm），X端为自由端，Y向采用周期性边界条件．原子总数为132个．我们选用函数表达式［6］：
对于铝而言，以上参数中
（2）初始条件的确定运动方程的求解要求给出粒子的初始位移和初始速度，而在一般情况下很难知道系统的初始条件，尽管给出的初始条件可以使得系统更快地达到平衡．即便如此，在选定的模拟时间足够长的情况下，就可以忽略初始条件［7］了．常用的初始条件［8］可以选择为：①令初始位置在差分划分网格的格子上，初始速度则从玻尔兹曼分布随机抽样得到；②令初始位置随机的偏离差分划分网格的格子上，初始速度为0；③令初始位置随机的偏离差分划分网格的格子上，初始速度也是从玻尔兹曼分布随机抽样得到．实验中由于铝金属的原子分布结构是面心立方结构，原子间距已知，则可以利用几何关系确定各个原子的初始位置，然后将各原子初始速度清零．
（3）确定步长
可计算得到在第n＋1步时所有的粒子所处的空间位移．计算在第n＋1步［9］时所有粒子的速度为
动能和速度标度因子为
计算每个时间步系统的总势能U，考察系统达到平衡所需要的时间．
在分子动力学模拟过程中，如何确定所用的时间步长是个十分关键的问题，它直接决定了模拟过程所需要的时间．当然，步长越小，最后的误差也就越小，但是对系统进行模拟的驰豫时间相应也就越长了．因此要根据实际选择步长，我们选用时间步长为0．5fs．
（4）计算n＋1步粒子的速度计算将速度乘以标度因子的值，并让该值作为下一次计算时第n＋1步粒子的速度．
（5）重复计算直至系统达到平衡返回第三步，开始第n＋2步的模拟计算．重复200次，一直到系统达到平衡．
（6）对系统进行位移加载我们选取的是对端面直接施加位移的方法，由此可以得到关于纳米单晶铝杆拉伸模拟的应力—应变关系曲线，即进行位移加载时，固定最下面两排原子，对最上面的两排原子进行位移加载，每次加载长度为第二排至倒数第二排原子间距的0．000 5倍．
每次加载后使系统平衡100步，让系统经历一个准静态加载的过程．每次加载并让系统平衡后，计算系统所受的力．应力—应变关系是按如下方法得出的：对应于每个载荷步，我们用所施加的总载荷除以横截面积，得到对应于这个载荷步的应力．实际上它只是横截面上的平均应力，对应于这个载荷步的应变也是平均应变．1．2 计算结果与分析
（1）原子的分布见图1．
图1 原子分布
该模型图与铝原子的面心立方结构相一致，说明模型的正确性．由此可得各原子的初始位移．
（2）循环次数与系统总能量的关系曲线见图2．
图2 循环次数与系统总能量的关系曲线
从图2中容易看出势能变化曲线随着循环次数的增加而减小，最终趋于平稳，即在200步左右，系统已基本上达到平衡．
（3）进行位移加载时，加载次数与系统受力的关系曲线见图3［10］．
图3 加载次数与系统受力的关系曲线
由于加载是位移加载，而原子量级的长度单位为纳米，即10－9 m．又每次加载长度为第二排至倒数第二排原子间距的0．000 5倍，由图1加载长度为
20×0．000 5＝0．01纳米级，即10－11 m 单位量级．图中力与加载次数的关
系曲线大致为一斜率为正的直线，表明系统受力与加载次数呈正相关关系．由于每次进行的是位移加载，此图也可以说明应力—应变的变化曲线也为这样的直线，
斜率为杨氏模量．这与非纳米状态下的铝的拉伸性质相符合．
1．3 精确度分析
图中纵坐标的单位取自参考文献［10］中图6．易知铝金属的杨氏模量大致为
72GPa（不考虑温度影响）．我们简单对图3的斜率数值计算检测模拟的精确度：取图中（5，1．985），（15，2．062），（25，2．14）三点计算直线斜率（因为两点确定一条直线，这里取3个间距均衡的点）：
则斜率的平均值为
相对误差为
由于在精确度分析中具体的运算没有按照求杨氏模量的公式严格导出，例如截面积大小．只在一维尺度数量级上进行了拉伸程度分析，可以看到误差较大．这里的问题只与计算方法相关，与建立的动力学模型没有太大关联，可以看到此模型简单易求，是有一定的优点的．
2 结论
实验采用了分子动力学的方法模拟了纳米单晶铝杆的拉伸实验．在模拟中选择用紧束缚TB势函数作为原子间的势函数．经过模拟，得出了铝纳米杆相应的应力—应变曲线，可以验证其曲线与其实际宏观条件下的形状基本一致［11］，可见此方
法得到了良好的实际效果．
【相关文献】
［1］张媛媛，张文飞，王鹏．不同温度下纳米单晶铜杆拉伸的分子动力学模拟［J］．兵器材料
科学与工程，2007，30（4）：52－57．
［2］文玉华，周富信，刘曰武，等．纳米晶铜单向拉伸变形的分子动力学模拟［J］．力学学报，2002，34（1）：29－36．
［3］文玉华，朱如曾，周富信，等．分子动力学模拟的主要技术［J］．力学进展，2003，33（1）：65．
［4］李世影．Ni纳米团簇热力学性质和纳米碳管的动力学性质研究［D］．扬州：扬州大学，2008．
［5］郭晓光．单晶硅纳米级磨削过程的分子动力学仿真研究［D］．大连：大连理工大学，2008．［6］宋海洋，查新未．碳纳米管／金复合材料拉伸力学性能的分子动力学模拟［J］．稀有金属
材料与工程，2008，37（7）：1 217－1 220．
［7］李作祥．煤中甲烷吸附与解吸的分子动力学模拟可行性分析［J］．矿业快报，2008，15（6）：51－53．
［8］曹鸣．硫酸钡晶体压电常数的分子动力学模拟［D］．武汉：华中科技大学，2007．
［9］魏育新．金原子纳米团簇反常尺寸效应的分子动力学模拟［D］．重庆：重庆大学，2002．［10］宋海洋，查新未．碳纳米管轴向压缩行为的分子动力学模拟［J］．微纳电子技术，2007，37：246－253．
［11］付称心，陈云飞．硅纳米杆拉伸的分子动力学模拟［J］．传感技术学报，2006，19（5）：1 697－1 699．。