采用morse势模拟单晶铝预制裂纹的扩展过程

采用morse势模拟单晶铝预制裂纹的扩展过程

——作者水平有限，仅供学习交流使用，如有不足还请指出

        对于morse势与LJ势类似，都是无键相连的相互作用势，在lammps的语法为:

pair_style          morse cutoff  

pair_coeff          N M D0 α r0 (cutoff)  

        这里的N、M代表指定使用morse势的原子类型，D0为结合能，α为弹性模量，r0为平衡位置的原子间距，cutoff为截断半径；各量具体单位参见手册。

       对于大部分初学者来说，很纠结的问题是如何选定截断半径cutoff的值，本文将详细说明截断半径的选取。

      那首先明确何为截断半径cutoff？我们知道对于任何物质都是由很小的粒子组成，而粒子与粒子之间存在相互作用力。那我们想通过计算机来实现模拟，粒子相互作用的计算是不可避免需要解决的问题，而且粒子之间距离越近，作用力也就越大，也就是说任何一个粒子所受到的力都来自于周围粒子对他的合力。所以要计算一个粒子受到的力，首先要确定周围的粒子数量，那多大范围内的粒子才是他的周围?所以我们引入截断半径的概念，也就是以中央粒子为球心，以截断半径为半径画一个圆，只计算圆内的粒子对他的合力，圆外粒子对他的作用忽略不计。如图所示。

        明白了截断半径的含义，那就进入我们的主题如何选取截断半径？首先，很遗憾告诉大家的是对于截断半径没有一个明确的选取标准。可能一些朋友看到这里后暗骂一声垃圾便退出了。

       但是呢，我想通过一个例子去选取不同的截断半径对比得到的结果去说明截断半径的选取。

        该例子来自对2010年2月发表在《系统仿真学报》王晓娟[1] 的论文复现。

         in文件如下：

        只进行建模部分的代码书写，后续的裂纹扩展只需加上velocity设定原子速度即可。我们着重讨论势函数部分截断半径的选择。

1.        #复现Al的裂纹扩展  

2.        units               metal   

3.        atom_style          atomic  

4.        dimension           3  

5.        boundary            s s p  

6.        neighbor            2.0 bin   

7.        neigh_modify        every 1 delay 0 check yes   

8.        #建模  

9.        lattice             fcc 4.05  

10.     region              box block 0 40 0 16 0 10  

11.     create_box          1 box  

12.     create_atoms        1 box  

13.     mass                * 26.982  

14.     #势函数  

15.     pair_style          morse 7.0  

16.     pair_coeff          * * 0.2703 1.165 3.253  

17.     #划分区域  

18.     region              lower block INF INF INF 1 INF INF  

19.     region              upper block INF INF 15 INF INF INF   

20.     group               lower region lower  

21.     group               upper region upper   

22.     group               boundary union lower upper  

23.     group               mobile subtract all boundary   

24.     #预制裂纹  

25.     region              crack block INF 5 7 9 INF INF   

26.     delete_atoms        region crack compress yes   

27.     #能量最小化  

28.     min_style cg  

29.     minimize            1e-15 1e-15 5000 5000  

30.     #温度初始化、弛豫  

31.     velocity            mobile create 300 859646  

32.     fix                 1 all nve  

33.     thermo              100  

34.     thermo_style        custom step temp lx ly lz pe   

35.     run                 5000  

36.     write_data          AL-crack.xyz  

        可以看到这里的截断半径为7.0，那我们用OVITO将模型可视化后可以看到该模型基本符合原文献。

#势函数  

pair_style          morse 7.0  

pair_coeff          * * 0.2703 1.165 3.253  

        可能就有朋友问了，为什么7.0就行，5.0不行？8.3不行？10.5不行？那我们分别将截断半径改为5.0、8.0、10.5分别运行之后看一下得到的模型。

#将截断半径改为5.0  

pair_style          morse 5.0  

pair_coeff          * * 0.2703 1.165 3.253  

      截断半径选取太小，粒子受力计算太小都不能维持晶体结构。

#截断半径改为8.3

pair_style          morse 8.3 

pair_coeff          * * 0.2703 1.165 3.253  

        截断半径有些大，预制裂纹有愈合的倾向。

#截断半径改为10.5

pair_style          morse 10.5  

pair_coeff          * * 0.2703 1.165 3.253  

       可以看到截断半径太大，导致预制裂纹闭合。

       那说了这么多，截断半径到底怎么选？难道要一个一个的试吗？不，有方法，接下来就是激动人心的方法介绍。不知道大家有没有注意到该模型的预制裂纹的尺寸是多大？可能有些朋友想回到前面看一下in文件。不用再往前翻了，我已经将预制裂纹部分代码copy到这里啦。是不是很贴心。哈哈哈哈。

#预制裂纹  

region              crack block INF 5 7 9 INF INF   

delete_atoms        region crack compress yes   

         可以看到预制裂纹为5a*2a*10a(a为铝的晶格常数)的一个区域。换算成units box单位就是20.25*8.1*40.5（按a=4.05近似计算）也就是说该裂纹的三个尺寸最小的一个大约为8.1（这里的单位为埃）。是不是有种恍然大明白的感觉？对了，截断半径不能超过它预制裂纹的最小尺寸。那超过了会怎样？假如你的截断半径为8.5，在预制裂纹表面处有个原子受到的力，是不是以它为中心8.5半径的粒子给他的力都算进去，也就是裂纹的上下表面有相互作用，是不是会有裂纹愈合的倾向？

参考文献：

[1]王晓娟,朱宝全,王红梅.温度对单晶铝裂纹扩展影响的分子动力学模拟[J].系统仿真学报,2010(2):534-536

 

旺-旺啊     

2022年11月15日

北京