冲击压缩下金属钯的结构相变*

刘泽涛 陈博 令伟栋 包南云 康冬冬 戴佳钰

(国防科技大学物理系,长沙 410073)

钯作为典型高压标定材料,研究其在极端条件下的结构变化以及热力学性质具有广泛需求并充满了挑战,特别是冲击加载下钯的固-固相变过程研究仍然匮乏.本文基于嵌入原子势,使用经典分子动力学方法从原子角度揭示了冲击载荷加载下钯的结构相变路径,在0—375 GPa 的压力区间观察到一系列复杂的结构转变特征,从初始的面心立方(FCC)结构,至带密排六方(HCP)结构的层错体心立方(BCC)结构,直至完全熔化.在沿 〈100〉 晶向冲击下,在70.0 GPa 发现了FCC-BCC 相变过程,远低于之前研究中静高压的结果.此外,还发现了冲击方向依赖的相变点,在沿着 〈110〉 及〈111〉 晶向冲击时FCC-BCC 相变压力分别增加至135.8 和165.4 GPa,同时相比完美晶体,引入缺陷会使FCC-BCC 相变压强值有20—30 GPa 的增幅,并通过势能分布的分析予以验证.本文发现冲击加载下钯的FCC-BCC 相变压力大大降低的特殊现象,为钯在高压实验等极端条件下的应用提供了新的理论认识.

1 引言

极端条件下物质的结构及其相变动力学是行星物理、冲击波物理、凝聚态物理以及核物理中具有重要研究价值的课题之一.近年来随着高压实验设备与技术的快速发展,常压下无法观测的新物理、新现象和新规律在高压研究中逐渐被发现.快速变化的外部载荷可以产生超高的温度压力环境[1,2],从而形成与常压下不同的热力学状态[3-5]、独特的物理性能和亚稳态结构[6-10],并且材料在高温高压下引发相变后的物理、力学性质相比初相有着本质的不同[11],因此研究极端条件下材料的热力学性质及结构变化动力学过程对于科学研究和工业应用都具有重要的意义[12,13].钯是铂族金属中熔点最低,密度最小的金属元素,因其卓越的吸氢能力,钯在电子信息制造业及医疗器械制备等领域[14-16]均有重要应用,并且成为20 世纪80 年代末极具关注的冷聚变实验的关键成分[17].得益于优异的热力学性能与独特的物理特性,钯吸引了大量的关注[18-21],然而目前大多数研究集中于常态下,对于极端条件下钯的性质则鲜有所见.常压下钯为典型的面心立方(face-centered cubic,FCC)结构金属,在高压下仍能维持稳定的FCC 结构,因而在高压实验中常常被用作压力校准器.然而高温高压下钯的热力学性质以及结构变化过程仍不明晰,只有少数实验和理论研究在温度-压强同时作用下的相变演化,这使得研究极端条件下钯的相变过程具有重要意义.

前人对钯在高温高压条件下的相变研究主要集中在熔化等固-液相变过程.Folies 和Adams[22]使用由他们团队[23]早先开发的多体嵌入原子势函数(embedded atom method,EAM)计算金属钯的固相和液相的热力学性质,通过蒙特卡罗模拟和准简谐近似(quasi harmonic approximation,QHA)方法开展固相和液相的自由能计算,并利用固相和液相的吉布斯自由能曲线的交汇点来确定钯的熔点.2011 年刘中利等[24]通过基于EAM 势的分子动力学模拟计算了钯的熔化曲线,利用8000 个原子构造固-液共存系,计算得到了钯的熔点,同时还使用25 万个钯原子,沿着〈100〉晶向利用非平衡分子动力学方法(nonequilibrium molecular dynamics,NEMD)模拟冲击波诱导熔化,相比两相共存系的模拟结果在226 GPa 下产生了18.3 %的过热熔化.2013 年Errandonea[25]利用激光加热金刚石顶砧(diamond-anvil cell,DAC)的实验方法首次测量了钯的高压熔化曲线,与之前刘中利等[24]的理论计算结果有明显的分歧,实验发现在30 GPa 下与理论计算存在约为500 K 的差距,但与1999 年Jeong 和Chang[26]基于分子动力学的计算结果在20 GPa 以上的高压范围内拟合良好.2015 年刘中利等[27]提出一种冲击熔化的模拟方法并证明能够精确有效地确定金属的熔化曲线,仅用756 个钯原子测量了熔化曲线,且与两相法模拟结果符合,在理论研究方面显著提升了计算效率,但仍未解决理论与实验在熔化曲线上的差别.对于钯的结构性质方面刘中利等[24]利用第一性原理对钯在高压下的相变过程做了详细的分析,通过比较FCC 结构、密排六方结构(hexagonal close packed,HCP)和体心立方结构(body-centered cubic,BCC)在QHA下的吉布斯自由能,发现在0—500 GPa 的范围内不同温度条件下的FCC 和HCP 结构的自由能差值曲线没有交点,得出FCC 结构的钯原子在高达500 GPa 和5000 K 的条件下是结构稳定的.至今为止,关于钯在极端条件下的结构相变规律研究仍然匮乏,并且对于钯相变的理论研究大都停留在准静态压缩.而物质在冲击压缩这种快速载荷作用下,其内部动态行为表现与静态[28]或准静态载荷作用下截然不同.阐明冲击压缩下钯材料结构转变的微观动力学过程并探索其相变的物理机理,不论是从理论上进一步理解钯的物理性质,还是将其更好地应用于高温高压环境下,都具有重要意义.

实时观测冲击下材料结构的响应极具挑战性,伴随着计算机水平的快速发展,至今为止已经发展出多种计算材料学的模拟方法,其中分子动力学模拟已成为在微纳米尺度上理解材料动态响应的重要研究方法[29-32],分子动力学能够模拟材料中原子的实时运动轨迹并获得动态演变细节,在常规实验手段难以达到的时间和空间尺度上对材料的结构特征和相变现象进行模拟研究,对于理解材料在原子尺度上的冲击特性是极具优势的.本文利用分子动力学模拟研究了钯在高温高压条件下的热学性质,用数十万原子通过多尺度冲击压缩技术模拟了不同条件下的FCC 理想晶体钯在冲击速度范围为4.2—9.0 km/s 下的动态响应结构演化.研究发现,当沿着〈100〉晶向冲击时,压强达到35.6 GPa时FCC 结构的钯原子含量呈现下降趋势,同时HCP 结构在FCC 结构中以层错的形式出现;当压强达到约为57.5 GPa 时BCC 结构开始明显增加,HCP 结构消失;当压强约为70.0 GPa 时,BCC 的含量快速增加,并在94.9 GPa 时达到最大值,这意味着FCC 结构完成了向BCC 结构的转变过程.相较于前人静态压缩模拟的计算结果[23],本文采用更为贴近实验的动态压缩结果,揭示了理想晶体钯在1382.3 K 和94.9 GPa 时完成了FCC-BCC 的结构转变现象.此外,当沿着〈110〉晶向进行冲击压缩时,在冲击压强约为135.8 GPa 时发生了FCCBCC 的相变过程,而沿着〈111〉晶向时对应相变压强为165.4 GPa,显示出强烈的晶向依赖性.另外在初始的晶体结构中引入缺陷后,相较完美晶体的相变压强大约有20—30 GPa 的增幅.本文揭示了单晶钯在动态冲击压缩下FCC-BCC 相变规律,通过模拟原子尺度层面的动力学过程为实验研究提供了理论支撑,更加深入地理解钯的热力学性质及结构特征,从而达到拓展应用领域、挖掘其潜在特性的目的.

2 方法

2.1 多尺度冲击技术

冲击波的模拟采用LAMMPS (large-scale atomic/molecular massively parallel simulator)[33]程序包里的多尺度冲击技术(multi-scale shock technique,MSST)[34]模块.系统中所有原子都遵循修正的拉格朗日关系来更新位置和速度,以将整个体系约束在冲击雨贡纽条件下.当受到冲击速度为us的冲击加载时,物质遵循以下的雨贡纽关系:

这里up,ρ,p和e分别是粒子速度、密度、压强和单位质量的能量,下标为0 表示冲击波前的物理量.将守恒方程(1)—(3)作为额外的约束加入原子的拉格朗日量中,使整个体系中的原子在达到平衡时处于冲击速度为us的冲击波扫过后的准稳态.其拉格朗日量可以表示为

这里T,V,Q和u分别是单位质量的动能、势能、类质量参数以及比体积.通过调整盒子的单轴长度来实现应力的加载,并在雨贡纽方程的约束下进行演化模拟冲击加载的过程,直到系统达到平衡(=0),整个体系处于波后的准稳态.相较于动量镜等基于大体系的非平衡方法,MSST 方法通过在运动方程中额外加入雨贡纽方程作为约束,能将体系约束在冲击状态附近,因而本工作中所统计的热力学状态都是体系处于冲击波扫过后的准稳态,并且经过长时间的平衡弛豫过程以达到体系内部的应力与温度均在空间保持一致.在模拟中通过设置冲击波的速度来体现体系不同的冲击状态,每一个冲击速度对应于冲击波传播后达到准稳态的压力值和温度值,P-T-V的状态方程变化符合雨贡纽关系.具体的描述可参考Reed 等[34]的文章.

本文的MSST 计算中,使用Foiles 等[23]开发的多体嵌入原子势对钯原子间的相互作用进行描述.通过调整模拟盒子单一方向的长度变化来对体系施加单轴方向的应力,冲击压缩的方向确定为z轴方向,钯金属体系的初始结构为FCC 结构,计算晶胞尺寸为30×30×80,晶格常数为3.8907 Å,由288000 个原子组成,对3 个方向均采用周期性边界条件.冲击加载前先使用NPT 系综在300 K 和1 bar (1 bar=105Pa)下对整个体系进行50 ps 的弛豫.平衡后对体系加载冲击速度为4.2—9.0 km/s的冲击波,间隔为0.2 km/s.设置类质量参数为3600,人工黏度系数为0.0903,转换因子为0.01 以使得系统更快达到波后的准稳态.整个MSST 模拟持续500 ps,待体系充分平衡后统计轨迹及热力学数据.在实现冲击波沿不同晶向进行传递时,通过调整沿冲击波方向的晶格基矢来确定晶向,沿不同晶向体系加载冲击之后,待冲击波扫过之后的体系达到准稳态再获取波后结构信息;模拟构建含缺陷晶体类型时,通过随机删除指定占比的原子来设置初始晶格体系的缺陷率来模拟点缺陷的情况,用来进行本次工作中的对比计算.

2.2 结构分析方法

为对冲击过程中的结构进行分析,采取了两种结构分析方法,借助开源可视化工具OVITO (open visualization tool)[35]实现.其一为径向分布函数(radial distribution function,RDF)g(r),RDF 定义为以一个粒子为中心,在半径r—r+dr的空间范围内发现另一个粒子的概率,用以表征结构的无序程度.其定义为

其中V为系统的体积,N为系统中的原子数,rij为i原子到j原子的位矢.RDF 已被用作揭示材料结构相变的有效方法[36-38],曲线中的峰表示相邻原子相互靠近的概率.RDF 曲线中第二峰和第三峰出现融合趋势并形成“短波浪”形状,该特征的出现可被看作是FCC 结构向BCC 结构转变的一个判据[39].

另一种结构分析方法是由Honeycutt 和Andersen[40]提出的公共近邻分析(common neighbor analysis,CNA)方法,该方法已经广泛地运用在固-固及固-液相变等结构分析当中,并被证实是行之有效的方法[36,41].该方法是根据原子周围相邻原子键的拓扑结构进行分类,具体可以通过区分近邻原子的配位数来分辨结构类型,例如对于BCC 结构其第一近邻配位数为8,而FCC 结构和HCP 结构的第一近邻配位数为12.近邻配位数可以通过截断半径r处的RDF 积分得到:

其中ρN为配位原子数密度,即粒子个数与体系空间体积之比.对于传统CNA 方法,往往通过选取全局统一的截断半径对所有原子晶体结构进行表征.但对于本次工作中涉及到的多相系统,以及考虑到高温或者原子应力分布不均引起的位置波动,传统CNA 方法中全局固定的截断半径难以描述不同环境下的原子结构.因此使用自适应公共近邻分析(adaptive CNA,a-CNA)的方法[42]对处于不同原子环境中的粒子分别计算截断半径,该方法是传统CNA 的扩展,为多相系统的结构分析提供更加精准的判断方法.文献[43,44]详细介绍了这种方法.

3 结果与讨论

在分子动力学模拟中,势函数的准确性对于模拟结果的合理性至关重要[45].为了检验所选势函数的准确性和模拟方法的有效性,首先计算了钯的状态方程(equation of state,EOS).图1 给出了冲击波在沿〈100〉晶向传递过程中体积比随压强的变化关系,同时给出了以往的冲击压缩实验以及理论计算[24,46-49]中钯的EOS 关系,结果表明在0—375 GPa 的范围内使用EAM 势函数能得到与实验数据一致的EOS 预测结果.对不同大小体系的EOS 进行了收敛性测试,分别选用24000 (图中蓝色方形线条)和288000 (图中红色圆形线条)个原子的体系在相同条件下进行模拟,两支线条基本符合,因此证明在EOS 上24000 个原子体系已达到收敛标准.

图1 理论计算得到的雨贡纽P-V 曲线与冲击实验数据的比较Fig.1.Hugoniot P-V curves compared with the shock experimental data and theoretical results.

利用a-CNA 方法分析了不同冲击速度加载下体系的结构变化.不同结构原子的占有率随冲击压力的变化如图2(a)所示,随着冲击速度的逐渐增大,压强大于35 GPa 时,FCC 结构的钯原子在整个体系中所占比例开始下降,HCP 结构的原子开始显现,在图2(b)中观察到HCP 结构的钯原子在FCC 结构中以层错的形式存在.当冲击压强增加到57.5 GPa 时,BCC 结构开始增加,HCP 结构随之消失,同样从图2(b)可以看出,BCC 结构也是以层错的形式出现在FCC 结构中,没有出现BCC 的整体块状结构.在70.0—94.9 GPa 的压强范围中BCC 结构的含量快速增长直到其含量达到最大值,此时FCC 和HCP 的原子比例几乎降至0,意味着FCC 结构向BCC 结构转换完成.当压强超过232.0 GPa 时,晶体的结构特征全部消失,与图3(b)所示的P-T图中对应的间断点(228.3 GPa,5488.8 K)相对应,表明此刻出现了熔化,其微观表现为相对体积和原子平均能量都随着温度升高而线性增大,达到某一临界温度时体积骤降,原子平均能量骤升,之后两者又继续保持线性增大.该压强间断点与刘中利等[24]对〈100〉晶向的冲击模拟获得的数据(压强和温度于5855 K 和226 GPa时出现熔化)一致.总之,本文的结果准确地揭示了沿着〈100〉晶向冲击压缩过程中的相变过程及热力学响应.与小体系静态计算[24]中钯FCC-BCC 相变压力大于500 GPa 的结果不同,我们发现沿〈100〉晶向加载冲击波后,在较低的相变压力70.0 GPa时便能出现FCC-BCC 的结构转变.可见在冲击压缩这种载荷快速加载的过程中,温度与压力等热力学状态的耦合变化,以及原子剧烈位移的动态过程会降低钯的固-固相变压力.

图2 (a)沿 〈100〉 晶向冲击后不同结构特征的钯的原子分数与纵向应力的关系;(b)不同结构特征随纵向应力变化的可视化表示;(c)温度与纵向应力的关系,其中红色点线图为刘中利等[24]冲击熔化的分子动力学模拟数据Fig.2.(a)Fraction of atoms with different structural features versus longitudinal stress along the crystallographic direction 〈100〉 ;(b)visual representation of the variation of different structural features with longitudinal stress;(c)temperature versus longitudinal stress,the red dotted line plot shows the MD simulation data of Liu et al.[24] shock melting.

图3 沿不同晶向冲击后的结果 (a)〈100〉晶向;(b)〈110〉晶向;(c)〈100〉 晶向.其中实线表示冲击速度在5.0—7.8 km/s 范围内结构的RDF 图像,点划线表示与之对应的配位数曲线,黑色虚线位置指向初始状态下的第一峰的位置Fig.3.Results along different crystallographic directions:(a)〈100〉;(b)〈110〉;(c)〈111〉 .The solid lines indicate radial distribution functions of deformed micro-structure shocked at a range of shock velocity of 5.0—7.8 km/s.The dotted lines indicate the coordination number curves that correspond to the RDF,and the dashed lines points to the position of the first peak in the initial state.

为了探究钯在冲击压缩过程中影响其FCCBCC 相变压强的因素,首先模拟了冲击波沿不同晶向进行传播的过程.并通过RDF 和a-CNA 方法对冲击压缩后的结构进行了分析.图3 分别给出了沿〈100〉,〈110〉和〈111〉三个晶向加载冲击的RDF以及配位数曲线.对于〈100〉向,冲击速度为5.2—7.2 km/s 时可以观察到RDF 的第二峰与第三峰的融合趋势,呈现为“短波状”[39]的特征,在r=4.25 Å 处(对应于RDF 的波谷)附近的配位数呈现出“平台”状特征;而当相变发生后,该结构特征消失并且配位数的变化曲线趋于平滑,表明此径向处的原子配位发生改变,说明发生了相变,对应的温度和压强范围在1382—1699 K 和95—108 GPa,这与a-CNA 的判断基本一致.另外,同样的波状特征在〈110〉和〈111〉起始出现的冲击速度分别为6.6 及7.2 km/s,分别对应冲击压力135.8 和182.1 GPa,意味着相比〈111〉晶向,沿〈110〉晶向冲击会在较低压力下发生FCC-BCC 的相变,并且均比〈100〉晶向需要更高的冲击压力.除此之外,随着冲击速度的增加,在沿着各个晶向冲击下,第一峰出现左移(图3 中虚线位置),并有逐渐拓宽的趋势,峰值高度减小以及整体的展宽,在低速冲击(5.0 km/s＜us＜5.5 km/s)区域,结合图2(a)中的a-CNA 诊断方法及原子结构图发现,沿〈100〉晶向冲击会使体系出现带HCP 的层错结构,沿〈110〉晶向会出现FCC-BCC-HCP 多相共存的纳米结构,这些微观结构的出现会导致冲击后RDF 峰展宽;在高速冲击(us＞7.2 km/s)区域,RDF 峰的展宽与强度降低是由于熔化或接近熔化导致的非谐效应造成的,表明随压力升高体系的非谐效应增大,表现出向液相转变的趋势.为进一步分析沿不同晶向的冲击结构,利用a-CNA 方法对不同冲击方向的波后体系进行了对比.图4 给出了沿〈110〉和〈111〉晶向在冲击压缩后的结构随冲击压强的变化情况,将a-CNA 采样点中BCC 结构增长首次趋于稳定的压力值作为FCC-BCC 结构的相变点,综合图2(a)中沿〈100〉晶向冲击的结果,沿〈100〉,〈110〉和〈111〉这3 个晶向加载冲击对应的相变压力依次为70.0,135.8 和165.4 GPa,这与上述的RDF结果表现基本一致,相比其他两个晶向,沿〈111〉晶向冲击更不容易发生相的转变.

图4 利用a-CNA 方法得到的沿不同晶向冲击的不同结构比例随压强的变化 (a)〈110〉晶向;(b)〈111〉 晶向Fig.4.Fraction of atoms with different structural versus compressing stress along different crystallographic directions obtained by a-CNA:(a)〈110〉;(b)〈111〉.

另外,由于在实际材料生长及应用过程中,很有可能会引入一些缺陷,因此还探讨了包含不同缺陷率的初始结构对钯在冲击压缩下FCC-BCC 相变压力的影响,结果如图5(a)所示.可以看出,缺陷的引入给状态方程带来了影响,在同等压力下相比完美晶体,含缺陷的体系会获得更高温度,在3%缺陷率时达到最大值.同样,通过a-CNA 跟踪了不同缺陷率体系沿〈100〉晶向在冲击后结构特征随着压力的变化情况,如图5 所示.相比无缺陷晶体(图2(a))可以发现,随着缺陷的引入相变压力整体呈现增大趋势,但与缺陷的占比并不遵循正相关关系,相变压强在3%缺陷率时达到最大值99.1 GPa,在5%缺陷率时又下降至84.8 GPa.为进一步分析引入缺陷增大相变压力的原因,分别计算了完美晶体、3%缺陷率以及5%缺陷率的缺陷晶体在冲击前原始状态及冲击后发生相变时的原子势能,冲击速度均设置为6 km/s,图6 给出了冲击前后的势能差(ΔEp)在体系中的分布情况,用以表征体系中FCC 转变为BCC 结构所需要克服的能垒.通过图6(a)的颜色分布以及图6(b)的概率分布图可看出,无缺陷晶体在冲击后的势能相比初始状态是均匀增加的,引入3%和5%的点缺陷后势能增幅变大且分布不均匀,因此相比完美晶体,含有缺陷的体系需要克服更大的能垒用于完成结构相变,而在缺陷为5%时所需要克服的能量有明显减少,这与本文a-CNA 结果一致,因此引入点缺陷会增加钯在冲击压缩下FCC-BCC 的转变压力.

图5 (a)不同缺陷率的晶体结构在冲击压缩后温度与压力的关系;(b),(c),(d)分别表示冲击压缩下含不同缺陷率体系的原子类别占比随冲击压强的变化Fig.5.(a)Comparison of Hugoniot P-T curve among crystal structures with different porosities defects;(b),(c),(d)comparison of phase transition processes of atoms of system with variation porosities defects under shock compression.

图6 (a)完美晶体与含缺陷晶体在冲击前后的势能差值分布图;(b)不同体系中势能差值的概率分布曲线Fig.6.(a)Distribution diagram of potential energy difference between perfect crystal and defective crystal before and after compression;(b)probability distribution curves of potential energy differences in perfect crystal and defective crystal.

4 结论

基于分子动力学的方法,模拟了冲击压缩下钯的结构相变过程,并探讨了沿不同晶向冲击及不同缺陷率的初始构型对相变过程的影响.首先得到了与实验一致的雨贡纽状态方程,并结合RDF 图谱的特征峰和配位数曲线以及a-CNA 的结构分析手段,在0—375 GPa 的压力区间观察到了一系列结构变化特征,从初始的FCC 结构至出现HCP 层错,到FCC 转换成BCC 结构,最后直至熔化.在沿〈100〉晶向的冲击加载下,发现了低至70.0 GPa的FCC-BCC 相变压力,本文的状态方程与以往的动态实验及NEMD 的分子动力学模拟结果拟合良好,低于之前静态计算中所报道的结果,由此表明了冲击压缩过程会降低钯的固-固相变压力.此外还发现,沿其他晶向加载冲击波会增大FCC-BCC的相变压力,沿〈110〉与〈111〉方向体系的FCC-BCC相变压力能分别增大至135.8 和165.4 GPa.最后,在初始结构中引入不同缺陷率的点缺陷,发现对体系的状态方程产生了影响,相变压力会明显提升20—30 GPa,通过比较冲击前后体系的势能差值情况可以看出,含有缺陷的结构在发生相变的过程中需要克服更大的能垒,因此点缺陷的引入会增加FCC-BCC 的转变压力.通过此次工作系统地研究了钯在高压动态载荷加载下的热力学性质与相变过程,发现了冲击压缩加载下降低固-固相变压力的特殊现象,不仅为高速冲击下的相变过程提供参考方向,也能为钯的高压实验等极端条件下的应用提供参考和理论依据.