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

English Abstract

全文HTML

--> --> -->

金属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单晶在冲击加载下的同构相变行为.

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

加载 晶向	x 轴晶向及 尺寸/nm	y 轴晶向及 尺寸/nm	z 轴晶向及 尺寸/nm	模型 原子数
[001]	[100]	[010]	[001]	5.00×106
	25.8	25.8	255.5
[011]	[100]	$ [0 1 \bar1] $	[011]	4.83×106
	25.8	25.5	251.9
[111]	$ [\bar1 \bar1 2] $	$ [1 \bar1 0] $	[111]	4.88×106
	25.9	25.6	253.7

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

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

2

-->

3.1.波阵面结构

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

由图1可以看出, 在较弱的冲击加载下, 冲击波阵面为单波结构, 仅有一个弹性压缩波. u_p = 50 m·s^–1或u_p = 100 m·s^–1时, 在[001]晶向加载下, 波阵面处的粒子速度等物理量缓慢上升, 表现为斜波; [111]晶向加载下的波阵面与前者相比明显更加陡峭; 特别是[011]晶向加载下波阵面几乎为垂直上升可视为突变. 当u_p = 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偏位错, 深蓝色为全位错, 浅蓝色为梯杆位错. u_p = 150 m·s^–1, 时间t = 80 ps
Figure2. 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. u_p = 150 m·s^–1, t = 80 ps.

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

up[001] /m·s–1	150	200	250	300	400	500
fcc	66.1	44.1	24.6	11.8	2.1	0.4
bcc	30.1	51.6	71.4	84.9	96.1	97.4
其他	3.8	4.8	4.0	3.3	1.8	2.2

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

图 3 冲击加载后晶体中的微结构. 冲击晶向为 (a) [001]; (b), (c) [011]; (d), (e) [111]; (c), (e)中隐去了fcc结构原子. u_p = 200 m·s^–1, 时间t = 80 ps
Figure3. 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). u_p = 200 m·s^–1, t = 80 ps.

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

3.2.冲击Hugoniot关系

图4展示了由模拟得出的不同晶向加载冲击Hugoniot关系, 其中图4(a)为U_s-u_p关系, 并将其与实验结果^[20]进行了对比. 在不同晶向加载下, 前驱波由于单晶弹性的各向异性, 差别比较明显 (其中[001]晶向较低强度加载下(u_p^[001] = 50—250 m·s^–1), 前驱波的上升阶段几乎重合, 前驱波的扰动从波阵面的前沿开始, 故以此作为前驱波的位置). 相变波虽各有差异, 但与实验的结果^[20]差别不大. 特别是[111]晶向, 当加载强度较低(u_p^[111] = 200—300 m·s^–1)并且冲击波形表现出前驱波、塑性波、相变波三波结构时, 与实验符合得较好.
图 4 单晶Ce的冲击Hugoniot关系　(a) U_s-u_p关系; (b) P-u关系. 实验数据来自文献[20]. (b)中的“”表示统计标准差
Figure4. 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.

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

3.3.冲击固固相变

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

图 6 冲击相变的相变温度-压力状态与静水压相变条件
Figure6. 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 冲击前后的径向分布函数
Figure7. 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. u_p = 200 m·s^-1, 时间t = 80 ps　(a) [001]; (b) [011]; (c) [111]
Figure8. 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. u_p = 200 m·s^-1, t = 80 ps.

在[001]晶向加载的模拟中, 利用PTM分析发现了γ-Ce向fcc和bcc混合物的转变. 如表2所列, 其中bcc占比与冲击强度呈正相关fcc占比反之, 特别是u_p^[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结构中比值c/a的分界点($c/a\!=\! (1+\sqrt{1/2})/{2} \!\approx 0.8536$), 此亚稳态是由采用的势函数^[14]决定的. 若原子处于此亚稳态, 将难以区分究竟是fcc还是bcc. 所谓的fcc和bcc混合物实际上是此亚稳态附近的结构, 是介于fcc和bcc之间的体心四方(bct). 此处[001]晶向加载的第二冲击波实际上是fcc(γ-Ce) → bct亚稳态结构的相变波.
图 9 [001]晶向不同加载强度下相变波后的微结构 (a) u_p^[001] = 150 m·s^–1; (b) u_p^[001] = 200 m·s^–1; (c) u_p^[001] = 250 m·s^–1; (d) u_p^[001] = 300 m·s^–1; (e) u_p^[001] = 400 m·s^–1; (f) u_p^[001] = 500 m·s^–1
Figure9. Microstructure of the sample after phase transition shock along [001] with listed piston velocity: (a) u_p^[001] = 150 m·s^–1; (b) u_p^[001] = 200 m·s^–1; (c) u_p^[001] = 250 m·s^–1; (d) u_p^[001] = 300 m·s^–1; (e) u_p^[001] = 400 m·s^–1; (f) u_p^[001] = 500 m·s^–1.

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

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

利用EAM模型势的分子动力学模拟表明, 单晶Ce在一定强度的冲击载荷下其冲击波阵面会呈现出多波结构, 并具有明显的晶向取向相关性: 在[001]和[011]晶向加载下表现为双波结构, 依次为前驱波和相变波; 而在[111]晶向加载下表现为三波结构, 在前驱波和相变波之间还存在一个γ-Ce中的塑性波, 该塑性波对应于偏位错和层错的生成. 相变波的U_s-u_p 以及P-u Hugoniot关系与实验符合得较好.
[011]和[111]晶向冲击加载下, Ce的γ → α相变压力值相较于静水压加载条件下偏低, 单轴冲击加载的剪切效应对Ce的γ → α相变起到了促进作用.
在[001]晶向冲击加载的模拟中发现了γ-Ce向fcc和bcc混合物的转变. 基于变形路径的分析表明, 此结果来自于转变路径中存在一个介于fcc和bcc之间的bct结构亚稳态, 该亚稳态是本文采用的Ce原子间作用势决定的.