基于Lammps仿真模拟下铁在拉伸过程中动力学研究

闫琴琴

山西江淮重工有限责任公司 山西省晋城市 048000

摘要：分子动力学模拟软件是一种利用分子动力学方法来模拟分子系统的计算机程序，具有广泛的应用和重要性，可以用于研究各种分子系统的结构和性质，对其进行可视化和分析。Lammps是一个广泛使用的分子动力学模拟软件，可以用于模拟各种材料的结构和性质，如金属、聚合物、纳米颗粒等。它具有高度可扩展性和灵活性，可以进行并行计算和自定义代码开发，需要一定的编程知识来使用和定制。下面是基于Lammps仿真模拟下Q235钢的分子动力学拉伸研究。

关键词：分子动力学；铁；拉伸；应变率；抗拉强度

一、初始模型的建立

典型的金属原子结构排布方式主要有体心立方(bcc)、面心立方(fcc)和密排六方(hcp)等结构。纯铁在常温下的排布方式是体心立方结构，这也是我们采取的晶体结构。根据Q235钢材中元素含量的占比，采用Atomsk软件建立初始原胞，并进行扩胞（a0=0.2855nm）产生20000个Fe原子，通过lammps软件按比例随机替换Fe原子的方式建立Q235钢模型，这种方法在其他人的工作中也被使用。C原子个数40、Mn原子个数280、Si原子个数69，最终得到如图4-2所示的含碳量为 0.2% 的钢板合金体系。其中X、Y、Z坐标轴分别对应面心立方的[100]、[010]、[001]晶向，横截面（XY面）为长方形。

二、模拟过程

使用分子动力学模拟进行拉伸试验时，通常通过改变体系的：应变率，时间步长，体系温度，模型尺寸，拉伸方向，边界条件等获得不同的测试结果，其中时间步长主要与计算机性能有关，时间步长越短，通常计算结果越精确，同时也会耗费更多的机时，在分子动力学计算中，选取合适的时间步长非常重要，我们通过前期研究确定时间步长为1fs，这也是近几年对金属材料做分子动力学模拟最常用的时间步长之一。

模拟过程中使用Plimpton等开发的大规模原子/分子并行模拟器LAMMPS 研究拉伸下Q235纳米钢板的应力状态。对模型施加三维周期性边界条件消除模型的边界效应，在z轴[001]方向上施加拉伸载荷。模拟初始阶段，采用共轭梯度法对构建的钢材合金体系进行能量最小化，获得能量最小结构； 利用Nose-Hoover热浴法，采用npt系综控制体系的温度保持在300K，控制X、Y、Z轴方向上的压力为0，设置时间步长为0.001ps，进行30ps的弛豫，得到体系最稳定结构。将此时的模型作为单轴拉伸的初始位形，采用Aslam,I开发的MEAM势，该势函数准确的描述了Fe-Fe、Fe-Mn、Fe-Si、Fe-C、Mn-Si、Mn-C、Si-C、Mn- Mn、Si-Si、C-C之间的相互作用关系并用其测试了Q235纳米钢薄膜的弹性常数如图所示

B是体模量、E杨氏模量、G剪切模量

弹性常数	C11	C12	C44	B	E	G
实验值(GPa)	226	140	116	169	-	-
模拟值(GPa)	226	129	111	162	203	78

模拟结果与实验得到结果对比，计算结果在合理范围内。拉伸阶段，沿z轴方向上对模型进行均匀拉伸，应变率为1×10^9/s，用npt系综控制拉伸过程中的温度，并确保x，y轴方向上的压力为0，直到拉伸到模型扩大30%。在整个模拟过程中，采用Lammps软件进行模拟计算，借助软件OVITO功能板块如共近邻分析(CNA)、中心对称参数(CSP)、位错提取法(DXA)、原子应力云图、原子位移云图和多面体模板匹配(PTM)等可视化手段来观察原子位置和内部结构的变化过程。

在弛豫过程中发现，整个模型中的原子在位置上并没有发生太大的变化，只是在表面的个别原子会向模型边缘聚集，并与边缘的BCC结构的原子形成新的BCC结构。由于结构尺寸很小，弛豫后的原子的总能量和弛豫步数的变化趋势如图所示。

图1给出了在控温控压操作下体系的能量和温度变化过程。在经历100万步的弛豫后，能量从一开始的以一个较高的值为基线上下波动，但在95万步左右产生了突变，最后体系能量逐渐稳定，表明体系能量最小化，达到平衡状态。由图可以看出，在弛豫过程中温度基本保持稳定状态，这是由于采用了Nose-Hoover方法进行了控温操作。

三、结果与讨论

1.应力应变关系曲线

图2为纳米钢板在300k时沿[001]方向,应变率为1*10^9的拉伸应力-应变曲线。根据图中的曲线变化特征可将变化过程分为，弹性变形阶段，塑性变形阶段，脆性变形阶段以及延性断裂阶段，红色箭头是弹性变形阶段，黑色箭头是塑性变形阶段，蓝色箭头是脆性变形阶段，绿色箭头是延性断裂阶段。

（１）弹性变形阶段是,0<ε<0.06。在此过程，拉伸应力与应变呈现出正比例关系，随着应变的增减，应力值急剧上升。Fe-C-Mn-Si钢材在应力达到5.25GPa时开始缓慢上升，相对于纯 Fe 的弹性屈服极限8.7GPa 低了约 3.45GPa。这种弹性屈服极限降低的现象是C、Mn、Si原子添加到Fe晶格中产生的结果。除在锰，碳，硅原子周围形成非晶结构以外，晶格排列完整，材料在拉伸过程中变得细长，原子分布均匀，未产生位错线等，为弹性变形阶段，也是可逆变形阶段，通过拟合弹性阶段的计算点数据得到，说明在纳米尺度下的胡克定律仍然成立。模型变化过程如图a和b所示。

（2）塑性变形阶段是0.06<ε<0.12。在此过程，应力先缓慢上升后保持不变得到抗拉强度5.51Gpa，C、Mn、Si原子过于的聚集程度也会对抗拉强度产生影响，随后应力下降。此阶段材料开始出现不可逆变形，部分晶格结构开始发生转变，材料多处都开始出现位错并且多处出现孪生变形，通过计算得到屈服应变为0.06.模型变化过程如图c所示。

（3）脆性变形阶段是0.12<ε<0.16。在此过程，随着应变的增加应力急剧下降，整个模型中出现了明显有序的孪生带，此阶段模型的原始结构被破70%，材料的承载能力基本下降到极限。如图

（4）延性断裂阶段0.16<ε。此过程，应力缓慢下降，且逐步趋于零。此时材料内，部会恢复原始BCC结构，材料开始出现断口，材料继续拉伸，最终发生断裂，断裂处的原子结构被彻底破坏。如图

2、单轴拉伸下的结构变化

材料的不同晶体结构力学性质相差很大。拉伸过程中晶体结构的转变对材料力学性能影响很明显。图中给出了在不断应变的过程中，BCC，FCC，Other，HCP，四种晶体结构在模型中所占比例。结合应力应变曲线以及图示对材料在拉伸过程中的晶体结构进行了如下分析：

1.当应变为0<ε<0.06，除Mn，C，Si，原子附近发生的晶格畸变以为，其余区域晶格排列完整，整体结构并未发生转变。因此在弹性变形阶段晶格未发生转变。

2.0.06<ε<0.12，此阶段，原本BCC结构开始快速转变为FCC结构和无序结构，由于铁中FCC结构和无序结构的能量高于BCC结构，因此在拉伸过程中新产生的无序结构和FCC快速消耗因形变而产生的系统晶格能量，当应变到达0.12时FCC结构不再增加，晶体的塑性变形阶段后期，此时还有极少的HCP结构。

3.0.12<ε<0.16，此阶段，部分FCC结构开始恢复为BCC结构，无序结构继续增加，此时的无序结构呈现短程有序，长程无序，材料开始出现明显的滑移区域，此时的结构相比第二阶段稳定些，滑移带边缘会有一些BCC结构，材料开始有断裂的趋势。

4.0.16>ε，此阶段，bcc结构，fcc结构和无序结构含量缓慢增长并趋于稳定，新增长的无序结构破坏了原有的FCC和HCP结构，最终导致材料出现裂纹，随着拉伸的进行，裂纹周围还出现孔洞，并不断扩大最终断裂。断裂后由于材料内部缺少应力支持，孪生带迅速消失，而断裂区域仍保持无序结构，系统其他区域则恢复为BCC结构，应力最终消失。

更进一步的微观结构分析发现：

1.当应变为0<ε<0.06时，Q235钢中Mn，Si，C原子周围会形成一些无序区域，约占总数的7%,这些无序结构破坏了晶体内部的均匀性，导致晶体局部出现应力差异，在Mn，Si，C原子周围会形成高应力区域，影响材料的稳定性。

2.当应变为0.06<ε<0.12时，材料内部的非晶区域开始扩展，部分BCC

结构转变为FCC结构和Other结构,FCC结构和无序结构稳步增长，原先的高应力区域逐渐转变为低应力区域，这些转变区域同时也是高应变区域，材料的稳定性进一步被破坏。

3.当应变为0.12<ε<0.16时，材料内部的非晶结构扩展到最大，约占总数的50%，部分FCC结构重新转变为BCC结构，导致FCC结构总数降低，同时在材料内部的应力积聚到最高，非晶聚集程度较高的区域，应变相比较其他区域更大，这也为后面断裂断裂埋下伏笔。

4.当应变为0.16<ε时，无序结构和FCC结构基本不再转变为BCC结构，材料内部的原子应力应变逐渐降低并恢复到正常水平。而无序结构在材料裂口区域不断聚集，破坏了材料原本的周期性，使得裂开不断增大，最终导致材料完全断裂。

由于Q235纳米钢板中，存在Mn，Si,C原子，这些原子导致体系内部存在许多无序区域，可以认为是材料内部的一些结构缺陷，在外力达到一定值后，这部分区域首先发生扩展，破坏原有的BCC结构并最终形成应力应变聚集区导致材料断裂。

四、应变率对纳米钢板拉伸破坏的影响

实验结果表明，同种材料的强度会随着应变率载荷的变化而改变、在单晶材料加载拉伸过程中屈服强度，拉伸强度，屈服应变等都会不同程度受到应变率的影响。当应变率超出某范围以后，材料出现了整体结构的非晶化现象，在结构设计中，考虑材料行为的应变速率和温度依赖性是非常重要的。了解金属在不同温度和应变速率范围下的变形行为，在金属成形、高速加工、高速冲击、穿透力学、爆炸-金属相互作用和其他类似的动态条件下具有重要意义。

应变率是宏观连续介质力学中衡量材料变形速率的指标，材料的强度随着应变率提高而增大，当应变率低于10^3~10^5/s时，材料强度随着应变率缓慢增长，当应变率高于10^3~10^5/s时，材料强度突然剧增并迅速趋向无穷。研究不同的应变率对于材料的力学行为有重要意义。

以10a×20a×50a的纳米钢板作为研究对象。对于边界条件、温度、势函数和系综同上一节保持一致，对弛豫后的纳米钢板进行均匀的拉伸，时间步长取为0.001ps，通过调整拉伸过程中的加载间隔以达到相应的应变率。选取五种应变率进行加载，分别为5×10^8/s、1×10^9/s、5×10^9/s、1×10^10/s、5×10^10/s，在拉伸过程中通过控制x，y方向的压力为零，以模拟单轴拉伸的效果。

图 是我们展示出的Q235钢纳米板在不同加载应变率下的应力应变曲线，通过分析可以发现，无论应变率大小如何，应力都随着应变速率的提高而增大，对于同一数量级的屈服强度相差不大，但当应变率增加数量级后，屈服极限应力从5.36GPa上升到了6.35GPa，材料的塑性变形能力也同步增强。

从图中可以看出，高应变率和低应变率下，材料变形趋势相同，但是高应变速率的杨氏模量相对会更高一些，

1.弹性变形阶段，应力应变曲线呈现线性关系。而较高的应变率，材料对应的屈服强度相对较高。但在这个阶段的应力应变曲线整体趋势相较来说比较一致，这是因为在这个阶段材料未发生大规模结构转变，主要还是BCC结构的铁原子在吸收拉伸过程中多余的能量，模型中由Mn，Si，C原子引起的无序结构未发生扩展。在单晶纯金属元素对应变率敏感性的研究中已经发现在弹性变形阶段，材料对应变率不敏感。

2.塑性变形阶段，不同应变率的塑性变形范围不同，随着应变率的提高，塑性变形阶段也在不断增大，并随着应变率提高材料会更早进入这一阶段，在这一阶段Mn，Si，C原子周围的无序结构开始扩展，这些结构对应变率的敏感性更高是导致这一现象的原因。有研究表明更高的应变率更容易使材料非晶化。

3.脆性变形阶段，应变率对材料的影响大致相同，但是在5*10^10应变率下，材料没有脆性变形阶段，我们进一步分析发现高应变率导致基本所有的bcc结构转变为其他结构，而bcc结构的缺失导致材料的脆性变形阶段消失。有研究表面BCC结构的纯铁在高应变率下有脆性变形阶段。而Q235在这一阶段基本全部转变未无序和Fcc结构，是导致这一现象的原因。

4.延性变形阶段，所有应变率都表现相同的趋势，不同应变率在这一阶段都会使材料发生断裂。在这一阶段，高应变率导致材料更迟断裂，主要是无序结构让材料的延性增强。

更近一步分微观结构分析发现，

应变率5*10^8

应变率1*10^9

应变率5*10^9

应变率1*10^10

应变率5*10^10

不同的应变率对比发现，在断裂前的材料内部位错线密度有所不同，在应变率不断提高的过程中，位错线密度会降低，当应变率到达5*10^10时，材料断裂前不会出现位错线，这些现象也是因为材料内部的结构转变比例随着应变率的增加而增加。如果断裂前剩余的BCC结构较少，就越不容易生成位错线。

图 是不同应变率加载条件下的对应钢材体系内bcc结构随应变不断增加的变化，从图中可以看出：

随着应变率的不断增加，弹性变形阶段和塑性变形阶段材料内结构转变基本保持相同的速率，但是到脆性变形阶段不同应变率表现出不同的现象，应变速率越大，最后剩余的BCC结构越少，同时断裂后BCC结构的恢复速率也会降低。在延性变形阶段BCC结构在总结构中数目基本保持不变，此时材料已经发生断裂。

我们的研究表明，体系内bcc结构最终占比数量与钢材的抗拉强度有直接的关系，bcc结构占比较少时对应得抗拉强度和断裂应变相对较小，bcc结构转化为其他类型结构越多，材料表现出更大的抗拉强度。分别将不同应变率下拉伸时结构内部变化的微观图截取出来如下图所示：在不同应变率下纳米钢板拉伸过程中达到屈服应力分别为7.973GPa、6.970GPa、6.880GPa、6.682GPa、7.709GPa 和 8.644GPa ；材料的局部结构由代表着 FCC 结构的绿色变为的代表着 HCP 结构的红色；有极少部分由代表着 FCC 的绿色变为代表着 BCC 的蓝色。结构在不断拉伸的过程中，逐渐出现滑移带，随着应变的增大，滑移带的数量逐渐增多，同时滑移带的位置逐渐发生移动，并且在无序结构周围出现了部分的交叉滑移，这些现象主要由Mn，Si，C原子引起。

五、尺寸对纳米钢板拉伸破坏的影响

大量试验表明对于同一材料同一加载状态的金属材料，改变材料的尺寸会对材料的力学行为有不同程度的影响。

为了研究尺寸对纳米钢板的力学性能的影响，在基本模型上模拟条件相同，截面尺寸分别设置为5a×5a、8a×8a、10a×10a、12a×12a，长度方向上分别设置为20a、60a、100a，如表1所示。

截面尺寸	杨氏模量(GPa)	体积模量(GPa)	泊松比	屈服强度(GPa)
	203.0	162.9	0.29	7.44
	201.7	159.8	0.28	8.02
	200.3	158.2	0.27	8.38

六、总结

本试验采用分子动力学模拟方法，详细分析了纳米钢板在拉伸作用下的微观层面的演化过程，并且详细的分析了纳米钢板在拉伸过程中，材料内部的结构转变，以及原子应力应变的分布。计算了不同应变率和不同截面尺寸下纳米钢板的力学性能，对比了应力应变曲线和弹性模量之间的关系，并提出屈服强度和应变率之间的关系式，并对其进行误差分析。并得出以下结论：

1.Mn，Si，C原子会在材料内部形成高应力场，并在其周围形成非晶结构，会导致材料力学性能降低。

2.应变率越高材料的抗拉强度和断裂应变越大。

3.材料的尺寸越大，抗拉强度越小。

作者简介：闫琴琴，女，出生1986年12月，籍贯山西晋城， 2009.1毕业于重庆邮电大学信息与计算科学专业本科毕业，工程系列工程师，就职于山西江淮重工有限责任公司，从事信息化工作。