刘小娥, 马保吉
(西安工业大学机电工程学院, 西安 710021)



Fi=-riU(r1,r2,…,rN)=
(1)
采用模拟体系为金属和无机离子两种模拟体系相混合的复杂模拟体系,模拟所选取的力场为Compass力场及Reaxff力场。用Compass力场和用Reaxff力场计算时体系总势能、粒子总势能的计算公式分别为[10]

(2)
E=Ebond+ELP+Eover+Eunder+Eangle+Edihedral+
EvdW+Eq
(3)


由于镁的密排六方晶格结构,选取其密排面Mg(0001)面进行研究。在构建Mg(0001)晶胞结构时,为了验证所采用方法和模型的合理性,首先构建镁的单晶胞,如图2(a)所示。然后采用Castep模块对镁的单晶胞进行结构优化[13],优化后三条边长尺寸为a=b=0.322 nm,c=0.517 nm,a、b为单晶胞底边边长,c为单晶胞的高,实验值a=b=0.320 nm,c=0.521 nm[14],优化后的结果与实验值非常接近。取Mg(0001)面,如图2(b)所。然后构建Mg(0001)表面超晶胞结构,体系的尺寸为31.56 nm×31.56 nm×27.38 nm,如图2(c)所示。体系包含1 584个镁原子。在模拟过程中,由于金属1~3层原子容易发生明显的弛豫现象,因此固定Mg(0001)晶胞四层以下的原子作为相体原子,上面的三层原子设置为可自由移动[15]。

图1 各离子模型结构优化结果Fig.1 Structure optimization results of each ion model

表1 不同离子优化结构与实验值[15]

图2 Mg(0001)晶胞结构的构建Fig.2 Construction of cell structure of Mg(0001)

腐蚀离子在Mg(0001)表面的扩散行为模拟采用正则系综(NVT)[17],模拟温度为310 K,控温方式为Andersen[16],范德华和库仑相互作用的计算方法为Ewald[17]。选取的截断半径为0.5 nm。模拟步长为1 fs,总时长为200 ps,每500步记录一次轨迹信息(帧)[18-19]。

图3 4种腐蚀离子在Mg(0001)表面的扩散体系模型Fig.3 Diffusion system model of four corrosion ions on the surface of Mg(0001)
在分子动力学模拟中,用温度平衡条件和能量平衡条件来判断体系平衡[20]。以Cl-离子溶液作用于Mg(0001)表面的模拟模型为例,图4、图5分别为在310 K温度下,Cl-离子溶液作用于Mg(0001)表面的模拟体系进行NVT模拟的温度变化和能量变化。……