单晶Ce 冲击相变的分子动力学模拟

第伍旻杰 胡晓棉

1) (中国工程物理研究院研究生院， 北京 100088)2) (北京应用物理与计算数学研究所， 计算物理国家重点实验室， 北京 100088)(2020 年3 月2日收到; 2020 年3 月25日收到修改稿)

金属Ce在室温条件下当压力达到约0.7 GPa时会发生一阶相变， 体积突变减小14%—17%， 相变前后两相分别为γ-Ce和α-Ce， 均为面心立方结构. 实验中发现冲击波在Ce中传播， 其波形存在明显的多波结构，依次为 γ-Ce弹性前驱波、γ-Ce塑性波、γ-Ce → α-Ce相变波. 基于新发展的金属Ce的嵌入原子势， 对单晶Ce的冲击相变行为进行了分子动力学模拟. 模拟结果表明， 在一定强度下， 单晶Ce中的冲击波阵面分裂为多波结构， 波形结构与加载晶向明显相关: 在[001]和[011]晶向加载下表现为双波结构， 依次为前驱波和相变波; 在[111]晶向加载下波阵面分裂为弹性前驱波、γ-Ce塑性波、γ → α相变波， 与已有的实验观察相一致.冲击波速的Hugoniot关系在低强度加载下与实验符合得较好. 同时在此冲击相变过程中， 应力偏量对相变起促进作用， 相较于静水压加载， 冲击加载的相变压力条件更低一些.

1 引 言

金属Ce由于其特殊的电子结构， 具有特殊的相图和丰富的相变行为. 如在室温条件下， 当压力达到约0.7 GPa时， 会发生γ-Ce → α-Ce的一阶相变， 相变前后的两相均为面心立方(fcc)结构， 体积突变减小14%—17%[1-4].

在对金属Ce的动态压缩冲击实验中，Pavlovskii等[5]发现冲击波在Ce中传播， 其波形存在明显的多波结构， 按次序分别为: 1) γ-Ce弹性前驱波; 2) γ-Ce塑性波; 3) γ-Ce → α-Ce相变波.与静态加载下的情况相比， Ce在动态加载下同样会发生α-γ同构相变， 相变的时间间隔在0.1 μs以内， 同时室温条件下α-γ冲击相变的临界压强约为0.76 GPa， 与静态加载相变相近. Borisenok等[6]和Simakov等[7]利用PVDF(聚偏二氟乙烯)测压计对Ce在冲击加载下的α-γ相变进行了实验测量， 亦得到了Ce在冲击下载下的多波结构， 相变压力为1.0—1.2 GPa， 与之前Pavlovskii等[5]的实验结果相近.

Yelkin等[8]和El’kin等[9]曾基于Aptekar’和Ponyatovskii的赝二元固溶模型， 给出了动态加载下Ce的α-γ两相全物态方程. 当动态加载的压力达到α-γ相变的压力以上时， 会形成多波结构，其中包含等熵波和冲击波. 得到的动态加载下开始发生α-γ相变的压力为0.75 GPa， 与之前Pavlovskii等[5]的实验结果相近. 然而当加载强度进一步增高， 物态方程的计算得到的压力明显低于实验结果. 对此的解释是实验中相变是在远离热动平衡的状态下发生的， 而此物态方程基于准静态热力学. 之后El’kin等[10]在此模型基础上又加入了ε相和液相， 物态方程的预测同样与静态和动态加载的实验符合较好. 动态加载下的非平衡效应以及弱冲击加载下金属相变更加复杂， Pan等[2]和Hu等[11]考虑到过程中应力偏量、塑性流动、应变率等的效应， 基于多相Steinberg-Guinan本构模型对Ce的α-γ冲击相变进行了数值模拟研究. 相较于El’kin等的Ce物态方程能更好地描述相变的非平衡效应.

伴随着计算机水平的发展， 已经发展出多种计算材料学模拟方法， 其中分子动力学(MD)模拟已成为在微纳米尺度上理解材料动态响应的重要研究方法. 如Kadau等[12]在对单晶Fe的MD研究中发现， 冲击的加载晶向会显著影响冲击波阵面的结构和冲击相变过程. 而仅有的对Ce同构冲击相变MD研究中， Dupont等[13]利用修正的Voter-Chen原子嵌入法势模拟得出冲击波在Ce单晶中传播呈弹性波-塑性相变波双波结构， 未能得到实验中的三波结构， 认为这与该模拟选用的Ce原子间作用势有关.

在本文中， 利用最新发展的面心立方Ce的嵌入原子势[14]， 通过大规模MD模拟研究了面心立方Ce单晶在冲击加载下的同构相变行为.

2 模拟方法

本文分别考虑了沿面心立方Ce单晶[001]，[011]， [111]晶向冲击加载过程. 相应设置了三种单晶样品模型， 具体模型尺寸参数如表1所列. 加载之前先对各个单晶模型在300 K温度和零压条件下进行等温等压(NPT系综)弛豫， 消除模型内部的应力并使样品达到热力学平衡态. 弛豫完成后，对模型的x和y方向采用周期性边界条件而对z方向采取独立边界条件， 再以动量反射的方式沿模型的纵向(z方向)加载冲击波， 即选取纵向一端的若干层原子为活塞， 以一定的速度up向模型内部运动， 由此产生冲击波.

上述弛豫和冲击加载过程均采用LAMMPS开源程序[15]， 时间步长为1 fs. Ce的原子间相互作用势采用嵌入原子法(EAM)的形式， 相关介绍详见文献[14]. 该势函数可以描述金属Ce的同构相变行为， 适合本文模拟研究. 对于模拟结果， 在分析局部晶格结构时， 由于共同近邻法(CNA)[16，17]在局部畸变较大时不理想， 在本文中采用多面体模板比对法(PTM)[18]， 再利用位错抽取算法(DXA)[19]更清楚地显示位错、堆垛层错等缺陷.

表 1 单晶Ce计算模型详细参数Table 1. Parameters of single crystal Ce sample for MD simulation.

3 结果与讨论

3.1 波阵面结构

图1为分别沿[001]， [011]， [111]晶向不同冲击强度下的密度剖面图(时间为80 ps). 可以看出，在不同加载强度下， 冲击波形具有单波或多波结构的特征: 其中在加载强度较低(即活塞速度up较小)的情况下一般呈现为单波结构; 提高加载强度，冲击波阵呈现出多波结构.

图 1 不同加载晶向和强度(up)的密度剖面图(t = 80 ps)Fig. 1. Density profiles of different loading orientation and strength (up) for t = 80 ps.

由图1可以看出， 在较弱的冲击加载下， 冲击波阵面为单波结构， 仅有一个弹性压缩波. up=50 m·s—1或up= 100 m·s—1时， 在[001]晶向加载下， 波阵面处的粒子速度等物理量缓慢上升， 表现为斜波; [111]晶向加载下的波阵面与前者相比明显更加陡峭; 特别是[011]晶向加载下波阵面几乎为垂直上升可视为突变. 当up= 150 m·s—1时，[011]加载仍为单波结构， 仅有一个弹性波;[111]加载的弹性波后平台有明显的起伏， 再对结构进行分析发现除了fcc晶格之外还有位错、层错等缺陷(如图2)， 且密度不超过7.4 g·cm—3(见图1)，表明未发生γ → α相变， 此处的起伏为塑性变形的结果.

图 2 [111]晶向冲击加载后晶体中的微结构 (a) 全部原子; (b) DXA分析显示位错并隐去了fcc结构原子. 绿色原子为局部fcc， 红色为hcp， 蓝色为bcc. (b) 中用管表示位错: 绿色为Shockley偏位错， 深蓝色为全位错， 浅蓝色为梯杆位错. up = 150 m·s—1， 时间t = 80 psFig. 2. Microstructure of the sample shocked along [111]:(a) All atoms are shown; (b) only non-fcc atoms are shown.Color coding: Green for local fcc atoms; red for hcp; blue for bcc. Dislocations are illustrated with tubes in (b): Green for Shockley partials; deep blue for perfect fcc dislocations;light blue for stair-rod dislocations. up = 150 m·s—1， t =80 ps.

如图1的密度剖面曲线， 在较强的冲击加载下各晶向加载的冲击波阵面开始表现为多波结构. 对于[001]晶向加载， 当up[001]= 150 m·s—1， up[001]=200 m·s—1， 以及up[001]= 250 m·s—1时， 冲击波阵面呈现为双波结构: 前驱斜波， 以及一个相变波. 相变波后密度显著增大， 主要由局部fcc和bcc组成的混合物(组分如表2所列). 对于其中的前驱波与相变波之间的平台还可以发现， 此处的粒子速度、密度、压强等物理量随冲击强度的增大而降低. 对于[011]晶向加载， 当up[011]= 200 m·s—1或up[011]=250 m·s—1时， 波阵面表现为双波结构， 二者之间存在一个较稳定的平台. 同时在前驱波后可发现明显的回跳特征. 如图3 (up[011]= 200 m·s—1)， 前驱波后的晶格结构仍为fcc， 并生成了位错、层错等晶格缺陷. 在第二个冲击波到达之前压缩比不超过8%(对应的fcc晶格常数约5.03 Å)， 未发生相变; 第二个冲击波后体积明显减小， 压缩比超过20%， 晶格结构仍保持为fcc， 表明此处发生了面心立方Ce的γ → α同构相变. 另外如图3(b)， 在相变波到达前， 已生成Shockley偏位错和hcp结构， 表明在相变之前就已经进入塑性阶段. 但是此处的塑性对压力、密度等物理量未产生显著的影响. 而在相变波后， 位错明显更加密集形成位错林， 其中主要为Shockley偏位错， 同时还生成了全位错和梯杆位错. 对于[111]晶向加载， 当up[111]= 200 m·s—1或up[111]= 250 m·s—1时， 波阵面表现为三波结构，依次为: 弹性前驱波、一个上升较缓满的斜波、一个较陡峭的压缩波. 如图3的结构分析显示， 前驱波前后的晶格均为fcc结构， 而前驱波峰值处的密度为7.32 g·cm—3， 仍为γ-Ce未发生相变， 前驱波峰值后同样发生了回跳; 而在第二冲击波的起始处发现有位错和堆垛层错生成， 表明此阶段是偏位错主导的γ-Ce中的塑性波; 第三冲击波后的部分压缩比超过了17%， 除了包含有位错、层错等缺陷， 仍主要由fcc晶格构成， 表明此处发生了Ce的γ → α同构相变， 第三冲击波为γ → α相变波. 同时相对于[011]晶向加载， 塑性在此对密度、压强等物理量存在显著影响， 图1中塑性斜波的上升虽缓慢但仍十分明显.

图 3 冲击加载后晶体中的微结构. 冲击晶向为 (a) [001];(b)， (c) [011]; (d)， (e) [111]; (c)， (e)中隐去了fcc结构原子.up = 200 m·s—1， 时间t = 80 psFig. 3. Microstructure of the shocked samples. The shock orientation is along (a) [001]， (b) and (c) [011]， (d) and(e) [111]， respectively. Atoms in fcc structure are hidden in(c) and (e). up = 200 m·s—1， t = 80 ps.

冲 击 强 度 继 续 增 大(up= 300 m·s—1， up=400 m·s—1， up= 500 m·s—1)， 如沿[001]晶向加载，虽波阵面仍为双波结构(前驱波以及一个相变波)，然而其中的前驱波变为稀疏波: 两波之间平台部分的密度、压强等物理量低于γ-Ce受冲击前的数值，粒子向冲击传播的反方向运动， 并且冲击越强该效应越明显. 这与之前提到的趋势([001]晶向加载下前驱波与相变波之间平台处的粒子速度、密度、压强等物理量随冲击强度的增大而降低)是一致的，并且在此进一步延伸为稀疏波和负压状态. 对于沿[011]晶向加载， 波阵面结构除了包含弹性前驱波和γ → α相变波， 并且相变波是在一个突起信号之后紧跟着出现的. PTM分析表明， 此突起峰值处有少量的bcc结构产生. 在[111]晶向加载下，波阵面仍为弹性前驱波、塑性波、γ → α相变波的三波结构. 弹性前驱波和相变波间距减小， γ-Ce塑性阶段进一步被压缩.

3.2 冲击Hugoniot关系

图4展示了由模拟得出的不同晶向加载冲击Hugoniot关系， 其中图4(a)为Us-up关系， 并将其与实验结果[20]进行了对比. 在不同晶向加载下， 前驱波由于单晶弹性的各向异性， 差别比较明显 (其中[001]晶向较低强度加载下(up[001]= 50—250 m·s—1)， 前驱波的上升阶段几乎重合， 前驱波的扰动从波阵面的前沿开始， 故以此作为前驱波的位置). 相变波虽各有差异， 但与实验的结果[20]差别不大. 特别是[111]晶向， 当加载强度较低(up[111]=200—300 m·s—1)并且冲击波形表现出前驱波、塑性波、相变波三波结构时， 与实验符合得较好.

图4(b)为MD模拟得出的各晶向加载的P-u关系. 可以看出模拟所得的P-u关系与Jensen等[20]的多晶Ce的实验数据符合得较好. 同时单晶中相同粒子速度对应的压强略高于多晶实验中的数值，这与李俊等[21]对单晶铁的实验观察是一致的. 另一方面， 加载强度较高时(up= 400 m·s—1或up=500 m·s—1)， 不同加载晶向P-u点分布较为一致; 反之加载强度较弱时(up= 200 —300 m·s—1)， P-u关系呈现出加载晶向相关性， 其中[001]与[111]晶向加载的数值比较接近， 而[011]晶向加载的数值与之有明显的差异(压力偏高).

3.3 冲击固固相变

图 4 单晶Ce的冲击Hugoniot关系 (a) Us-up关系; (b) P-u关系. 实验数据来自文献[20]. (b)中的“ ”表示统计标准差Fig. 4. Shock Hugoniot for single crystal Ce: (a) Shock speed vs. piston velocity; (b) pressure vs. particle velocity.Experimental data is cited from Ref.[20]. The symbol in(b) represents the statistical standard error.

图 5 不同晶向加载下的压力剖面图， up = 200 m·s—1， 时间t = 80 psFig. 5. Pressure profile for each loading orientation at up =200 m·s—1 and t = 80 ps.

图5 为t = 80 ps时的压力剖面图(up=200 m·s—1)， 可以看出相变冲击波前的压强更接近平衡相变过程中的数值[14]， 以此认为波前的状态即为冲击相变的条件. 由此得到的图6即为不同晶向和强度加载下的相变温度-压力状态. 相比于静水压加载MD模拟(采用相同的Ce原子间作用势)中的相变条件[14]， 可以发现同样温度下[011]和[111]晶向冲击加载的γ → α相变压力偏低. 表明在冲击加载下， 应力偏量对Ce的γ → α相变起促进作用[22]， 这与相关的实验结果[23]是一致的.

图 6 冲击相变的相变温度-压力状态与静水压相变条件Fig. 6. Temperature-pressure condition of shock-induced and hydrostatic phase transition.

Ce的γ → α相变为一阶同构相变. 相较于异构相变， 同构相变的局部晶格对称性不发生改变[24]，很难通过分析晶格结构类型来判断是否发生了同构相变. 对于一阶同构相变只能定量地分析是否发生了体积坍缩. 图7为不同晶向加载下冲击前和依次各冲击波后径向分布函数(RDF). 未受冲击结构RDF的前三个峰值对应的间距r均满足fcc结构中最近三阶近邻间距之间的关系， 且第二个峰值对应的距离r = 5.16 Å， 即γ-Ce在零压下的晶格常数. 前驱波后， Ce晶体被压缩同时保持fcc结构.[001]晶向加载下第二冲击波后的RDF发生明显变化， 不再具有fcc结构RDF的特征; [011]晶向加载下第二冲击波后的RDF仍符合fcc结构， 各峰值对应的距离r明显减小， 其中第二个峰值处的距离r = 4.83 Å， 说明此处的Ce晶体即α-Ce， 可以确认第二个冲击波为Ce的γ → α相变波;[111]晶向加载下的第二冲击波后RDF相比于此波前(前驱波后)几乎不变， 此冲击波是γ-Ce中的塑性波， 而第三冲击波后的RDF则类似于[011]晶向加载下第二冲击波后的情况， 仍然保持fcc结构中RDF的特征且第二个峰值对应的距离r = 4.83 Å， 这表明[111]晶向加载下的第三冲击波为γ → α相变波.

图 7 冲击前后的径向分布函数Fig. 7. Radial distribution function of the sample before and after the shocks.

图8 中筛选并隐去了原子体积较大的fcc结构， 这样就可以显示α， γ两相的交界面. 图8(a)中[001]晶向加载的两相交界较平整， 几乎是一个垂直于加载方向的平面， 并且也与相变波阵面重合(z ≈ 90 nm); 图8(b)中[011]晶向加载的两相交界(z = 50—61 nm)有明显的起伏， α， γ两相相互交错; 图8(c)中的加载[111]晶向， 从开始有α-Ce大量生成(z ≈ 100 nm)到几乎全部转变为α相(z ≈ 70 nm)存在厚度约30 nm的过渡区， 无法在MD尺度上确定分明的α， γ两相边界. 另外对于[011]和[111]晶向加载， 相比于γ → α相变发生前， 相变后区域的位错和层错的密度更高(如图3)，相变导致生成了更多的位错、层错等结构缺陷.

图 8 Ce冲击相变分界面， 隐去了原子体积较大的fcc. up =200 m·s-1， 时间t = 80 ps (a) [001]; (b) [011]; (c) [111]Fig. 8. Phase boundary of shock induced transition. Shock orientation: (a) [001]; (b) [011]; (c) [111]. The atoms of fcc structure with larger atomic volume are hidden. up =200 m·s-1， t = 80 ps.

在[001]晶向加载的模拟中， 利用PTM分析发现了γ-Ce向fcc和bcc混合物的转变. 如表2所列， 其中bcc占比与冲击强度呈正相关fcc占比反之， 特别是up[001]= 400和500 m·s—1几乎全部为bcc. 生成的fcc结构符合局部的γ-Ce → α-Ce的转变， 但是在此处生成bcc结构的温度、压力条件与金属Ce的实验相图[1]有所差异， 我们认为这结果与晶格结构的稳定性有关. 在此考虑两种fcc与bcc结构间的变形路径[25]: 四方变形路径和三方变形路径. 对于四方变形路径， 可以用比值c/a来衡量. 晶胞的初始结构为fcc， 如果用c来表示晶胞[001]晶向的晶格常数， 用a表示其[100]和[010]晶向的晶格常数， 此时c = a以及c/a = 1.变形过程中改变c/a同时保持晶胞的体积不变. 如果将晶格沿[001]晶向进行压缩， 比值c/a将不再等于1， 当c/a = 1/√2时结构变为bcc. 三方变形则可以视为原胞夹角的改变， 同时原胞的体积不变. fcc原胞的夹角α = β = γ = 60°; 如果夹角变大到90°， 则变形为sc (简单立方)结构; 夹角继续变大到109.47°结构变为bcc. 图9中沿[001]相变波后的微结构图表明， 晶胞的[100]， [010]， [001]晶向在冲击波后均未发生改变， 而是仅仅沿[001]晶向压缩， 比值c/a变小， 类似于四方变形路径的结构转变. 同时冲击压缩过程中[100]和[010]晶向的晶格常数不变都为a， [001]晶向被压缩晶格常数c变小， 同时原子体积也变小， 这与四方变形路径有所差异. 图10中比较了这两种变形路径的能量差异， 此变形路径(a不变， a = 5.16 Å)当c/a ≈0.86时存在一个亚稳态， 此数值十分接近fcc和bcc结构中比值的分界点， 此亚稳态是由采用的势函数[14]决定的.若原子处于此亚稳态， 将难以区分究竟是fcc还是bcc. 所谓的fcc和bcc混合物实际上是此亚稳态附近的结构， 是介于fcc和bcc之间的体心四方(bct). 此处[001]晶向加载的第二冲击波实际上是fcc(γ-Ce) → bct亚稳态结构的相变波.

以往的MD模拟研究中显示面心立方铯(Cs)单晶同构冲击相变前后的晶向会发生改变[26]，而Dupont等[13]对面心立方Ce单晶同构冲击相变的MD研究中则未发现此现象. 对于本文的模拟工作， 可以从图7比较不同加载情况下相变波后区域与未受冲击部分的晶向. 如之前提到的， 沿[001]加载相变前后的晶向不变; 而沿[001]和[111]晶向加载， 相变波后区域相较于对未受冲击区域的晶向有所差异， 并且相变后区域的晶向也不完全一致， 同时存在多个取向， 类似于多晶中的晶粒.

表 2 [001]晶向加载相变波后区域微结构组分(依据PTM分析)(%)Table 2. Fraction for each type of microstructure(analyzed with PTM algorithm) in the part after phase transition shock along [001] (%).

图 9 [001]晶向不同加载强度下相变波后的微结构(a) up[001] = 150 m·s—1; (b) up[001] = 200 m·s—1; (c) up[001] =250 m·s—1; (d) up[001] = 300 m·s—1; (e) up[001] = 400 m·s—1;(f) up[001] = 500 m·s—1Fig. 9. Microstructure of the sample after phase transition shock along [001] with listed piston velocity: (a) up[001] =150 m·s—1; (b) up[001] = 200 m·s—1; (c) up[001] = 250 m·s—1;(d) up[001] = 300 m·s—1; (e) up[001] = 400 m·s—1; (f) up[001] =500 m·s—1.

图 10 沿四方变形路径(原子体积不变)以及单轴压缩(a不变仅改变c/a)路径变形的能量Fig. 10. Comparison of the energy along tetragonal deformation path (atomic volume preserved) and the path of constant a.

4 结 论

利用EAM模型势的分子动力学模拟表明， 单晶Ce在一定强度的冲击载荷下其冲击波阵面会呈现出多波结构， 并具有明显的晶向取向相关性: 在[001]和[011]晶向加载下表现为双波结构， 依次为前驱波和相变波; 而在[111]晶向加载下表现为三波结构， 在前驱波和相变波之间还存在一个γ-Ce中的塑性波， 该塑性波对应于偏位错和层错的生成. 相变波的Us-up以及P-u Hugoniot关系与实验符合得较好.

[011]和[111]晶向冲击加载下， Ce的γ → α相变压力值相较于静水压加载条件下偏低， 单轴冲击加载的剪切效应对Ce的γ → α相变起到了促进作用.

在[001]晶向冲击加载的模拟中发现了γ-Ce向fcc和bcc混合物的转变. 基于变形路径的分析表明， 此结果来自于转变路径中存在一个介于fcc和bcc之间的bct结构亚稳态， 该亚稳态是本文采用的Ce原子间作用势决定的.