單晶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原子間作用勢決定的.