纯铁压缩的分子动力学模拟

《计算材料学导论》课程项目论文

纯铁压缩的分子动力学模拟?

摘要：本文介绍了利用lammps分子动力学模拟软件和VMD可视化处理软件模拟纯铁压缩过程的方法。得出了在压缩过程中，铁原子的排列和的运动规律。得到铁分子的运动细节，从而更加利于我们了解铁的性质。

关键词：纯铁压缩；lammps；VMD；分子动力学模拟；

New Roman：

Molecular dynamics simulation of the pure compression

WangQiang，LiuChang

(1. School of materials science and engineering, , Jiamusi University, Jiamusi 154007；)

Abstract：This paper introduces the simulation software and the method of VM visualization simulation software

of pure compression process using lammps molecular dynamics. In the compression process, the arrangement of the iron atoms and the movement of the iron atoms are obtained. The motion details of the iron molecules are obtained, which help us to understand the properties of iron.

Key words：Pure compression; lammps；VMD；molecular dynamics simulation;

1 引言（四号，黑体）

材料的力学性能随尺度发生变化尺度减小, 材料中缺陷存在的几率降低, 材料的强度提高同时尺度的变化可能导致材料内在变形竞争机制的改变, 例如多晶材料晶粒粒径在微米级以上时, 强度主要受位错强化机制控制, 而粒径进入纳米级后, 材料的变形主要来源于晶界滑移等机制原子尺度下, 微观效应占主导地位, 材料的理化、力学性能表现出与宏观不同、甚至相反的特性。Brenner发现金属单晶晶须拉伸强度与晶须直径呈反比,Fleck在微米级细铜丝的扭转试验中观察到尺寸效应纳米电机系统(NEMS)的出现同迫切要求了解纳米尺度下材料的力学行为, 当前从实验上较难获得详细的信息, 而分子动力学模拟可以提供相关细节。

分子动力学模拟，是指对于原子核和电子所构成的多体系统，用计算机模拟原子核的运动过程，并从而计算系统的结构和性质，其中每一原子核被视为在全部其它原子核和电子所提供的经验势场作用下按牛顿定律运动。我们用c语言编写程序，VMD动画演示得到原子在压缩过程中的变化。在控制温度不变的情况下，得到了金属铁分子的动力学模拟过程。通过不断压缩，观察每次压缩过程中以及压缩后铁原子的排列，得到金属铁的运动细节，从而更加利于我们了解铁的性质。

2 原理

材料的力学性能随尺度发生变化尺度减小, 材料中缺陷存在的几率降低, 材料的强度提高同时尺度的变化可能导致材料内在变形竞争机制的改变, 例如多晶材料晶粒粒径在微米级以上时, 强度主要受位错强化机制控制, 而粒径进入纳米级后, 材料的变形主要来源于晶界滑移等机制原子尺度下, 微观效应占主导地位, 材料的理化、力学性能表现出与宏观不

1

《计算材料学导论》课程项目论文

同、甚至相反的特性。Brenner发现金属单晶晶须拉伸强度与晶须直径呈反比,Fleck在微米级细铜丝的扭转试验中观察到尺寸效应纳米电机系统(NEMS)的出现同迫切要求了解纳米尺度下材料的力学行为, 当前从实验上较难获得详细的信息, 而分子动力学模拟可以提供相关细节。

分子动力学模拟，是指对于原子核和电子所构成的多体系统，用计算机模拟原子核的运动过程，并从而计算系统的结构和性质，其中每一原子核被视为在全部其它原子核和电子所提供的经验势场作用下按牛顿定律运动。我们用c语言编写程序，VMD动画演示得到原子在压缩过程中的变化。在控制温度不变的情况下，得到了金属铁分子的动力学模拟过程。通过不断压缩，观察每次压缩过程中以及压缩后铁原子的排列，得到金属铁的运动细节，从而更加利于我们了解铁的性质。

分子动力学通过直接模拟原子的运动过程, 使我们能够详细了解模拟对象的演化发展历史分子动力学模拟的一个关键在于原子势函数的选取原子势早期一般采用简单的对势, 但对势无法正确描述弹性常数, 其结果不理想世纪年代提出的镶嵌原子法、有效介质理论更客观地反映了原子间多体作用的本质, 可得到较合理的结果.认为体系总能量为：

Etol??rij??rij??2???Aexp[?p??r?1??]?????exp[?2q??r?1??]?ij?ii??0??0???j?i?1/2

有效地对分子动力学的数据进行后处理也是一个重要的研究方向, 因为对数据进行后处理, 获取有用的信息也是一个很烦琐的工作。现有的一些分子动力学软件功能不够强大, 只是偏重于某一行业, 通用性不高, 还有一些软件为自由软件, 可维护性不强。还没有出现集建模、求解、后处理于一体的分子动力学软件。因此, 针对特定的问题进行自行编程显得尤为重要。

通过分子动力学方法,模拟了铁分子锻造实验中的形变过程. 研究了晶体取向裂纹的形变特点和断裂机理,观察到各种形变现象,如位错形核和发射,位错运动,堆垛层错或孪晶的形成,纳米空洞的形成与连接等. 计算结果表明,裂纹扩展是塑性过程和弹性过程相结合的过程,其中塑性过程表现为由裂尖发射的位错导致的原子切变行为,而弹性过程的发生则是由无位错区中的原子断键所导致. 本研究采用VC++自行编程, VMD动画演示得到原子在拉伸过程中的变化,对所研究的问题进行求解。

其基本原理是使用一个含有有限个分子并且有周期性边界条件的立方盒子,从该体系某一设定的位能模型出发,通过计算机模拟求解微元中全部分子的牛顿运动方程,记录它们在各个不同时刻的位置、速度和受力等,然后统计得到体系的各种热力学、结构和性质,也就是由体系粒子的微观性质求算其宏观性质。

进行分子动力学模拟的首要问题是要得到准确的原子之间的相互作用势函数。知道原子间正确的相互作用势,从而必须知道相应的电子基态。

计算中根据以下基本假设:

（1）所有粒子的运动都遵循经典牛顿力学规律。

(2) 粒子之间的相互作用满足叠加原理。显然这两条忽略了量子效应和多体作用, 与

2

《计算材料学导论》课程项目论文

真实物理系统存在一定差别, 仍然属于近似计算。

II)VMD作图:VMD是一个强大的原子作图及动画演示软件,在运用C语言知识对上面求解计算后将得到一个文件，将此文件直接拖进VMD原子作图工具便可以得到原子运动动画。

3 方法

第一步: 即模型的设定，也就是势函数的选取。

势函数的研究和物理系统上对物质的描述研究息息相关。最早是硬球势，即小于临界值时无穷大，大于等于临界值时为零。常用的是LJ势函数，还有EAM势函数，不同的物质状态描述用不同的势函数。

第二步:给定初始条件，也就是给定原子的空间位置和速度。

运动方程的求解需要知道粒子的初始位置和速度，不同的算法要求不同的初始条件。如：verlet算法需要两组坐标来启动计算，一组零时刻的坐标，一组是前进一个时间步的坐标或者一组零时刻的速度值。

第三步：利用公式：

ri(t?dt)?ri(t)?dtvi(t)?dt2ai(t)2

计算在第n+1步时所有粒子所处的空间位置。计算在第n+1步时所有粒子的速度

vi(t?dt)?vi(t)?dt?ai(t?dt)?ai(t)?2，动能和速度标度因子：

??21(2N?2)kT?Ek??m(vin?1),????m(vn?1)2?2i?i???i?1/2

计算每个时间步系统的总势能U，考察系统达到平衡所需要的时间。

第四步：计算将速度乘以标度因子的值，并让该值作为下一次计算时，第n+1步粒子的速度。

第五步：返回第三步，开始第n+2步的模拟计算。重复几千次，一直到系统达到平衡。 第六步：对系统进行位移加载，每次加载后使系统平衡8000步，让系统经历一个准静态加载的过程。每次加载并让系统平衡后，计算宏观量。

4 In.shear文件修改

# 3d metal shear simulation

units metal boundary s s p atom_style atomic

lattice bcc 2.86 铁的点阵类型及点阵常数 region box block 0 16.0 0 10.0 0 2.828427

3

《计算材料学导论》课程项目论文

create_box 3 box

lattice bcc 2.86 orient x 1 0 0 orient y 0 1 1 orient z 0 -1 1 & origin 0.5 0 0 create_atoms 1 box pair_style eam/fs

pair_coeff * * Fe_mm.eam.fs Fe Fe Fe neighbor 0.3 bin neigh_modify delay 5

region lower block INF INF INF 0.9 INF INF region upper block INF INF 6.1 INF INF INF group lower region lower group upper region upper

group boundary union lower upper group mobile subtract all boundary set group lower type 2 set group upper type 3 # void

region void cylinder y 8 5 2.5 INF INF delete_atoms region void # temp controllers

compute new3d mobile temp

compute new2d mobile temp/partial 0 1 1 # equilibrate

velocity mobile create 300.0 5812775 temp new3d fix 1 all nve

fix 2 boundary setforce 0.0 0.0 0.0

fix 3 mobile temp/rescale 10 300.0 300.0 10.0 1.0 fix_modify 3 temp new3d thermo 25

thermo_modify temp new3d timestep 0.001 run 100

# shear 设置压缩速度

velocity upper set 0 -0.2 0

velocity mobile ramp vx 0.0 -0.2 y 1.4 8.6 sum yes unfix 3

fix 3 mobile temp/rescale 10 300.0 300.0 10.0 1.0 fix_modify 3 temp new2d

dump 1 all atom 100 dump.shear

#dump 2 all image 100 image.*.jpg type type &

# axes yes 0.8 0.02 view 0 0 zoom 1.5 up 0 1 0 adiam 2.0 #dump_modify 2 pad 4

#dump 3 all movie 100 movie.mpg type type &

# axes yes 0.8 0.02 view 0 0 zoom 1.5 up 0 1 0 adiam 2.0

4

《计算材料学导论》课程项目论文

#dump_modify 3 pad 4 thermo 100

thermo_modify temp new2d reset_timestep 0

run 8000 计算步数

5 结果

程序运行部分结果用VMD作图的结果：

图1 第2700步

由图一可见，锯齿状滑移带，滑移带之间成60度角，形成等边三角形。

图2 第5000步

由图二可见，随着压缩的进行，原有的三角形被破坏，形成更多小的等边三角形，滑移带进一步增加。

[参考文献] (References)

5