動態載荷下鐵相變的原子模擬研究進展

王 昆，肖時芳，祝文軍，陳 軍，胡望宇

（1. 湖南大學材料科學與工程學院，湖南 長沙 410082；2. 湖南大學物理與微電子科學學院，湖南 長沙 410082；3. 中國工程物理研究院流體物理研究所沖擊波物理與爆轟物理重點實驗室，四川 綿陽 621999；4. 北京應用物理與計算數學研究所，北京 100088）

絕大多數元素及其化合物在高壓下會發生多型相變，并且以馬氏體相變為主。在自然界和人類生產生活中的很多場景都存在高壓條件，如地心的極端高壓環境、隕石撞擊事件、車輛碰撞、金屬材料的沖壓塑性成型等。探測鐵的高壓多型相變有助于正確認識地心的多層結構[1–2]。金屬材料的多型相變/逆相變會改變材料的性質、塑性機制耦合、促進或抑制額外晶體缺陷的產生等，進而影響材料的斷裂及屈服行為[3–5]。鐵的 α?ε相變[6–7]是金屬高壓多型相變的典型范例，迄今為止，相關研究已持續近70 年，在實驗測試和機理分析方面取得了豐碩的成果。2000 年之前，相關研究主要以實驗與宏觀理論分析為主，在這方面已有不少經典綜述和專著報道[8–11]。隨著計算機技術的快速發展和原子尺度模擬方法的逐步推廣應用，鐵的沖擊相變研究逐步深入到微觀原子尺度，特別是原位在線診斷實驗技術的突破，極大地促進了該領域的蓬勃發展。近20 年鐵高壓相變的原子模擬研究可大致分為兩個階段：2011 年及以前是以Voter-Chen 勢為基礎的分子動力學（Molecular dynamics，MD）模擬研究，主要探討鐵的沖擊相變機制、相變產物類型以及相變形核與生長規律，這部分進展將在本文第1 節的開始部分進行簡要回顧；第2 階段是以改進的鐵原子間相互作用勢為基礎的MD 模擬，主要探討相變的變體選擇規則、相變與塑性、相變與層裂等問題，本文將重點梳理這一階段的研究進展。

在熱力學方面，高壓相變理論研究主要關注壓力對晶體結構穩定性的影響（如高壓相圖）；而在動力學方面，則主要關注壓力等因素對相轉變過程的影響（如相轉變路徑、相變成核條件、相變與塑性的相互作用等）[10]。前者從熱力學原理探索潛在新相存在的可能性，后者則針對實際相轉變過程，探討新相形核與長大的條件。兩者相輔相成，相變動力學研究不能脫離平衡相圖的約束，而平衡相圖的實驗研究也離不開相變動力學，例如：實驗中觀察到的相變-逆相變壓力回滯現象就是受相變動力學的影響[12]。目前，晶體高壓平衡相變的理論預測通常采用密度泛函理論結合適當的晶體結構搜索算法來完成，國際上已有大量的研究工作，尤其是AIRSS[13]、CALYPSO[14]、IM2ODE[15]等晶體結構預測軟件的出現，使晶體材料高壓結構相變的理論預測工作變得程式化。相變動力學研究主要依賴于實驗與大規模經典MD 模擬，尤其是激光加載結合原位飛秒X 射線衍射測量技術的出現[16]，使人們能夠直接觀測金屬在晶格尺度的動力學過程，結合MD 模擬結果，有望在晶格層次揭示整個動態響應過程及其物理本質，為高壓相變動力學物理建模提供可靠的理論支撐。

本文主要以鐵的 α?ε相變為例，綜述MD 模擬在相變動力學方面的研究進展。首先介紹經典MD 模擬研究所采用的鐵高壓相變的原子間相互作用勢函數，重點討論平面應變加載下鐵的相變機制及動力學的影響因素，包括加載強度、晶體各向異性、應變率、應變梯度、初始微結構等，然后簡要介紹 工程實際中的特殊加載對相變的影響，最后進行總結和展望。

1 鐵的原子間相互作用勢函數

隨著壓強的增大，鐵 α相的吉布斯自由能（Gα）逐漸接近甚至高于 ε 相的吉布斯自由能（Gε），在相變壓力閾值（pc）處兩者相等。由圖1 所示的鐵的相圖可知， α 與 ε 相邊界較陡，即pc受溫度的影響小。在零級近似下，令溫度為零，相變的熱力學條件實際上要求兩相的焓相等，即

圖1 鐵的高壓相圖[17]Fig. 1 High-pressure phase diagram of iron[17]

在MD 模擬中，能否正確描述鐵的 α?ε相變完全取決于所采用的原子間相互作用勢函數。由于大規模非平衡分子動力學（Non-equilibrium molecular dynamics，NEMD）模擬對計算效率的要求很高，金屬體系的勢函數通常采用嵌入原子模型（Embedded atom model，EAM）勢或Finnis-Sinclair（FS）勢，因此本文的討論僅限于這兩種類型的勢及相關研究工作。大約在2000 年以前，國際上絕大多數勢函數是針對常壓應用環境開發的，主要探討金屬在常壓下的熱力學行為，以擴散機制為主。將這些勢函數直接應用于高壓下金屬力學行為的模擬時，常常會遇到一些與實驗觀測不相符的現象。鐵是其中最典型的例子。2012 年以前國際上常用的鐵的勢函數主要有3 種：Voter-Chen 勢[18]、Meyer-Ente 勢[19]和Mendelev 勢[20]。由式(1)檢驗發現[21–22]，只有Voter-Chen 勢可以同時正確描述鐵的bcc（體心立方， α相）?hcp（密排六方， ε相）及其相變壓力閾值。然而，NEMD 模擬研究發現，Voter-Chen 勢只能正確描述[001]晶向的沖擊相變產物，其他晶向的沖擊相變產物主要為fcc（面心立方）相。盡管如此，基于Voter-Chen 勢的原子模擬研究在歷史上發揮了巨大的作用。例如，人們曾一度認為fcc 相是鐵 α?ε相變的中間產物[23–24]，即所謂的“三步轉變”機制，之后Hawreliak 等[25]結合沖擊波實驗和大規模NEMD 模擬方法排除了這種相變機制，基于密度泛函理論的計算分析同樣排除了fcc 相出現的可能性[26]。現在公認的鐵的相變機制[22,27–28]如圖2所示，即沿 α鐵的[001]晶向壓縮15%～18%，然后通過相鄰(110)面間的穿插，最終形成ε 相。受實際樣品尺寸的限制，在相變的第一步壓縮過程中可能出現一個額外的沿(110)[111]方向的微小剪切變形，使最終的變體繞c軸轉動了約5°[29]。對于大規模NEMD 模擬而言，由于在垂直沖擊波傳播方向一般采用周期性邊界條件，所以這種額外的小轉動不太明顯。

此外，無論是單晶還是多晶，Voter-Chen 勢的模擬結果都沒有觀察到位錯[22,30]，這與實驗不符。因此，要深入開展鐵的高壓相變動力學研究，必須發展新的原子間相互作用勢函數。如何構建能夠正確描述高壓相變現象的原子間相互作用勢函數是原子尺度模擬材料動態相變行為的一個亟待解決的重要應用基礎科學問題，與之密切相關的還有沖擊塑性（如孿晶、位錯等）的原子模擬問題。最近10 年，針對高壓應用環境，國內外在發展高壓原子間相互作用勢函數方面已經進行了一些重要的探索，除了本文將要介紹的鐵的勢函數以外，還有描述沖擊高壓下鉭形變孿晶的勢函數[31]、鎢位錯芯結構與形變孿晶的勢函數[32]、鉛塑性與相變的勢函數[33]等。本節將介紹2012 年以后鐵高壓勢函數研究的最新進展，主要包括兩種勢，即改進的Ackland 勢[21]和改進分析型嵌入原子模型勢（MAEAM）[34]，通過對比分析不同原子間相互作用勢的特點，總結出能夠模擬高壓相變與塑性行為的勢函數應具備的最基本特征，為將來進一步改進或發展其他元素的高壓勢函數提供解決思路。

圖2 單晶鐵沿[001]晶向沖擊的相變機制[27]Fig. 2 Phase transition mechanism of single crystalline iron under the shock along [001] direction[27]

1.1 改進的Ackland 勢

式中：f(rij)為原子電子密度。式（2）與式（3）中的求和遍及體系中所有原子，具體形式如下

式中：H(x)為階躍函數。模型中的勢參數見表1[36]。

表1 改進的Ackland 勢參數[36]Table 1 Parameters of modified Ackland potential[36]

1.2 MAEAM 勢

我國在鐵的高壓原子間相互作用勢函數研究方面與國外幾乎同時起步，湖南大學胡望宇課題組與中國工程物理研究院流體物理研究所祝文軍團隊合作，在MAEAM 模型勢框架下[37]發展了一個可以同時正確描述鐵沖擊相變與塑性的勢函數[34,38]。在該模型框架下，第i個原子的能量等于兩體相互作用能（V）、嵌入能（F）與能量修正項（M）之和

嵌入能F(ρi)表示為電子密度 ρi的函數

其中

式中：Pe為平衡態下的電子密度偏離函數值。針對鐵的高壓應用環境，優化的勢模型參數見表2[34]。

表2 MAEAM 勢參數[34]Table 2 Parameters of MAEAM potential[34]

1.3 高壓勢函數的驗證與評價

表3 改進的Ackland 勢與MAEAM 勢預測的鐵主要物理性質與有關第一原理計算及實驗結果的對比Table 3 Comparisons of key properties of iron predicted using modified Ackland potential and MAEAM potential with the first-principle-base calculations or experimental results

滿足表3 中的基本物性可以保證勢函數正確描述平衡態附近的基本熱力學行為和高壓結構相變，要正確模擬鐵的沖擊雨貢紐關系，還需要驗證高壓物態方程。圖3(a)為MAEAM 勢預測的鐵的3 種相的物態方程與有關實驗及第一原理計算結果的對比，可見MAEAM 勢可以很好地重現實驗高壓狀態方程直到200 GPa 以上。改進的Ackland 勢的預測結果見圖3(b)，其預測的高壓相的物態方程偏硬。

圖3 (a) MAEAM 勢預測的物態方程與實驗及第一原理計算結果的對比，(b) 采用MAEAM 勢與改進的Ackland 勢的計算結果對比[49,52–54]Fig. 3 (a) Comparison of equations of states predicted by MAEAM potential with experiments and the first-principle calculations,(b) Comparison of the results between MAEAM potential and modified Ackland potential[49,52–54]

圖4 MAEAM 勢(a)和改進Ackland 勢(b)預測的簡并<111>/2 螺位錯芯結構[21]Fig. 4 Degenerate core structure of <111>/2 screw dislocations predicted by (a) MAEAM and (b) modified Ackland potential[21]

圖5 彈性壓縮和相變的截面輪廓(a)，模擬的未沖擊樣品(b)和相變后樣品(c)的實驗衍射圖樣[36]Fig. 5 Line profile from elastically compressed and phase changed sections of samples (a), and the simulated diffraction patterns in an experimental geometry: (b) unshocked and (c) phase-changed sample[36]

2 平面加載下的結構相變

利用原子模擬方法在材料中實現高應變率變形主要有3 種途徑（見圖6）：第1 種是用一個“剛性”的活塞（飛片）從左側撞擊靶材料，在靶材料中產生一個向右傳播的沖擊波；第2 種是兩個靶碰撞，在兩個靶中分別產生向相反方向運動的沖擊波；第3 種是沿加載方向以一定的應變率均勻單軸壓縮材料，該加載方式不會在材料中形成沖擊波，便于分析材料力學行為機理的控制因素。通過第3 種途徑模擬的樣品相當于斜波壓縮樣品中的一個拉格朗日單元，在一些文獻[55]中也將該途徑等同于斜波（或“準等熵”）加載。采用這些模擬技術結合第1 節提到的鐵的高壓勢函數，近10 年在鐵的高壓相變機制、相變對織構的影響、相變形核閾值、生長規律以及相變與塑性的微觀耦合機制等諸多方面取得了大量突破性認識，呈現出包括材料物理、冶金、固體力學、沖擊動力學等多學科交叉融合的特點。

圖6 實現高應變率加卸載的3 種方法：(a)活塞加載，(b)對稱撞擊，(c)均勻單軸壓縮（“M”表示材料，箭頭表示活塞或材料的運動方向，灰色梯度表示材料受到碰撞時形成的沖擊波）Fig. 6 Three simulation approaches of de-/compression under high strain rates: (a) piston loadings, (b) symmetric impingements and(c) uniform uniaxial compressions, where “M” represents materials and the arrows show the motion directions of the piston and materials, and shock waves initiated around the impinging moment are illustrated by the grayscale

2.1 加載強度與晶體各向異性

圖7 單晶鐵的沖擊雨貢紐關系[34]Fig. 7 Shock Hugoniot relation of single crystalline iron[34]

用MAEAM 勢模擬的單晶鐵沿不同晶向的沖擊雨貢紐關系如圖7 所示（其中 σzz為應力的zz分量，v為粒子速度），結果表明沿不同晶向沖擊時會出現不同程度的過壓，沿[001]、[110]和[111]沖擊的平均相變壓力閾值分別為18.0、22.3 和23.8 GPa。為了與第1 節的平衡相變相關概念進行區分，本文用相變壓力閾值表示實際相變成核的起始壓力，用平衡相變壓力閾值表示平衡態的轉變壓力。根據第1 節的討論可知，bcc→fcc 的平衡相變壓力閾值高于[001]晶向的過壓均值，略低于[110]和[111]晶向的過壓均值，因此用該勢模擬兩個晶向的沖擊相變時，相變產物中會出現少量的fcc 相（如圖8(a)所示），但是仍以hcp 相為主，尤其在沖擊強度接近相變閾值附近。盡管單晶鐵的相變機制在晶體學上已經清楚（見圖2），但是由于沖擊加載使晶體發生了單軸變形，破壞了bcc 鐵原有的立方晶系對稱性，即原來等價的晶面或晶向相對于沖擊方向變得不等價。由此可知，在相同的相變晶體學機制下，不同的沖擊晶向將導致不同的相變路徑，產生不同的相變變體，得到不同的產物織構，最終影響晶體塑性行為。圖8(b)給出了沿單晶鐵3 個低指數晶向沖擊的相變產物類型，用hcp 相的基面（慣習面）表征。沿[001]和[111]沖擊時，最終的變體類型分別為2 種和3 種，相變后的樣品為多晶，與實驗觀測的單晶沖擊結果[27]一致。然而，沿[110]晶向沖擊時，變體類型只有一種，相變后的樣品仍為單晶，只是其中含有小角度晶界，與相變前的塑性有關[34]。最近的激光加載實驗也觀察到了這種單一變體的結果[56]。采用改進的Ackland 勢和斜波壓縮技術，Amadou 等[57]觀察到了單晶鐵相變前的形變孿晶，并且隨著沖擊強度的增大或沖擊波的形成，形變孿晶的體積分數逐漸減小直至消失，這種現象被認為與高沖擊強度下滑移主導的塑性有關。

圖8 (a)不同沖擊晶向與沖擊強度下單晶鐵的相變（紅色、黃色、淺藍色和深藍色分別表示hcp 原子、fcc 原子、bcc 原子和缺陷原子），(b)沖擊晶向與慣習面的關系[34]Fig. 8 (a) Phase transition of single crystalline iron under different shock directions and shock strength(The meaning of the colors is: red (hcp atoms), yellow (fcc atoms), light blue (bcc atoms) and dark blue (defects).); (b) relationship between the shock direction and habit plane[34]

2.2 應變率與應變梯度效應

2.2.1 應變率效應

激光驅動的斜波壓縮實驗表明[58]，鐵的相變壓力閾值（σα→ε）隨應變率（ ε˙）的變化規律用如下分段函數擬合

圖9 鐵的相變壓力閾值隨加載應變率的變化[61]Fig. 9 Transition pressure of iron as a function of uniaxial strain rate[61]

這種規律類似于屈服強度與應變率的關系[59–60]，因此σα→ε的變化趨勢在106s?1附近發生改變的機制也被認為與塑性屈服機制類似，即與熱激活的位錯運動到聲子拖拽機制的轉變相關。最近的斜波壓縮實驗中，通過相變平臺弛豫時間隨應變率的升高而縮短的現象，推測出相變成核與生長速率均增加[61]。采用NEMD 方法模擬單晶鐵在斜波壓縮下的相變時發現[62]，相變壓力閾值在應變率低于1010s?1時表現出的規律（見圖9）與式(13)描述的第2 段規律（即冪次律）類似，然而當應變率高于1010s?1時，相變壓力閾值隨應變率的變化規律逐漸偏離冪次律。這是因為在極端應變率的斜波壓縮下，波前存在很大的應變梯度，例如沿[001]晶向斜波壓縮單晶鐵的相變，波前應變梯度高達106m?1，使得耦合的應變梯度效應無法忽略。

沿[001]晶向斜波壓縮單晶鐵的研究結果表明[62]，相變的發生分為彈性單軸壓縮導致的晶格失穩、形核與生長直至完全相變幾個階段。隨著沖擊波的傳播，可以發現失穩點與形核點的位置逐漸接近，直到重合，這與形核過程的超聲速傳播有關[63]。在相變點附近，新相形核初期呈橢球狀，其中短軸沿某個(110)面法向，如圖10 所示。新相的生長分別通過沿(110)面內的<110>晶向滑移和沿短軸方向相變進行，其滑移方向的生長速度比短軸方向的相變速度快[64]。

盡管斜波加載可以避免沖擊加載引起的巨大溫升效應，但它仍屬于動態加載手段，除了在固體中產生一定的應變率外，還會形成一定的宏觀應變梯度。因此，檢驗應變梯度的影響顯得尤為重要。在原子模擬中，采用均勻單軸壓縮技術（見圖6(c)）不僅可以模擬斜波加載過程中產生的應變率效應，還能排除宏觀應變梯度效應的干擾。利用這種技術，邵建立等[65]研究了應變率分別為1.25 × 109和9.1 × 108s?1下相變的形核與生長，觀察到了柱狀晶粒與層狀結構的形成（見圖11），與Pang 等[63]在沖擊條件下觀察到的橢球狀晶粒的形核與生長不同。考慮到兩種結果都基于Voter-Chen 勢，為此推測出現這種差異的原因與應變率相關。Pang 等[63]發現晶粒的形核與生長主要發生在沖擊波前沿掃過之后的區域，該區域的宏觀應力和應變的變化很小，對應的應變率非常低。

基于MAEAM 勢模擬準等熵壓縮下單晶鐵相變的研究結果表明[62]：當應變率小于或等于109s?1時，臨界應變幾乎不變；而當應變率高于109s?1時，臨界應變快速增大。這主要是因為相變形核需要一定的孕育時間，當應變率較小時，孕育時間可以忽略，當應變率非常大，以至于在極短的孕育時間內材料仍會被進一步壓縮至更高的應力狀態才發生相變。然而，在斜波壓縮下，臨界應變隨應變率增加到一定程度后反而減小，說明極端應變率下應變梯度對相變形核有重要影響。

圖10 沖擊下鐵相變的微結構演化：(a)～(e)顯示相變機制，(f)～(g)顯示新相生長[63]Fig. 10 Microstructure evolutions during the process of shock-induced phase transition of iron:(a)–(e) phase transition mechanism, (f)–(g) growth of the new phase[63]

圖11 在應變率為1.25 × 109 s–1 的單軸均勻壓縮下hcp 相的形核與生長[65]Fig. 11 Nucleation and growth of hcp phase under unform uniaxial compressions with a strain rate of 1.25 × 109 s–1[65]

2.2.2 基于連續介質假設的應變-應變梯度耦合模型

位移的二階梯度與應變梯度分別定義為

其中，位移滿足uk=xk?Xk。由于應變梯度的出現，晶體中的質量密度（或比體積）將變得不均勻，因此與材料物性相關的密度除了依賴空間位置之外，還依賴計算密度所采用的體積，在梯度相關的高階連續理論中，該體積與所研究問題的特征空間尺度相對應，為了方便描述，本文稱之為特征體積。原子模擬結果表明[62,67]，晶格失穩首先從漲落最大的晶格原子開始，然后迅速波及周圍原子導致初始塑性或相變的出現。因此，解釋晶格層次的失穩時，特征體積應該定義在原子尺度。現考慮一個簡單的無限晶格，以X為中心的一個小體積（VX）的自由能的二階展開為

根據式(14)和式(15)，可以得到如下線性近似關系

接下來，考慮該小體積在虛的均勻小應變場擾動下自由能的變化，即式(18)對 δE取一次變分，截斷到二階小量，并注意應變梯度也會隨著應變場擾動的出現而發生一個小的變化，可得

將式(23)和式(24)代入式(25)，并化簡保留至應變或應變梯度擾動的二階小量，得到

式中： Γ為小體積的外表面，其單位法向量為n，B和 Λ的定義如下

式(26)中等號右邊的第2 個積分代表體積外部材料通過表面對該體積內材料的能量貢獻，若僅考慮材料內部（遠離邊界部分）的穩定性，表面的影響可以忽略，因此體自由能密度的變化為

式(29)中等號右邊第1 項等于平衡態下單位體積內外力所做的功，即

式(29)中等號右邊第2 項和第3 項為晶體中存儲的應變能密度。形變的晶體要在小的應變場擾動下保持穩定，其應變能密度在任意位置（對應于任意應變）處必須為正，即

則晶體穩定性條件（式(31)）變為

2.2.3 應變-應變梯度耦合模型的微觀基礎

要利用2.2.2 節中的應變-應變梯度耦合模型來判斷晶體的彈性穩定性，需要知道晶體在當前形變狀態下的各階彈性系數張量。然而，梯度相關的彈性系數在以前的高階連續理論中僅作為一個調和量，用于擬合一些已知的結果。王昆等[62,66]從原子尺度系統地提出了一種計算各階彈性系數張量的靜力學方法，并指出2.2.2 節所采用的現象性理論模型的微觀理解及其適用范圍。求解高階彈性系數的基本思路：將體系受到應變及應變梯度擾動的能量表示成應變和應變梯度的函數，為此首先需要將ΔU用位移及其各階梯度表示，然后再利用關系式(14)～式(17)。在假設近鄰原子數不變的前提下，經過一系列繁雜的代數運算，最終可得

根據2.2.2 節的分析，微觀勢能密度對應于連續尺度的應變能密度，因此，在等溫條件下，式(37)與式(18)中自由能密度的各項相對應，由此可得晶體高階穩定性條件（式(34)）所需的各階彈性系數張量。由式(40)、式(51)、式(45)及式(46)可知，對于中心對稱類型的晶體， τ和W恒為零。此外，值得注意的是，上述晶體各階彈性系數僅僅依賴于晶體的當前應變，與應變梯度無關，因為這里選取的特征體積是空間某格點處的單原子體積，由于微觀原子排列的離散性，應變梯度的概念僅僅在特征體積包含多個晶胞的情況下才存在，因此，此處提供的彈性系數計算方法僅僅適用于受到小應變梯度擾動的均勻形變晶體，或者用于不均勻形變晶體的局部彈性穩定性分析。前一種情況在進行靜態理論模型分析中經常遇到，如2.2.2 節的應變梯度穩定性分析；后一種情況在原子模擬結果分析中會遇到，本節的結果相當于提供了一種局域晶格彈性穩定性的分析方法。結合本節中各階彈性系數張量的計算方法與2.2.2 節建立的應變-應變梯度耦合模型，預測了斜波壓縮單晶鐵在極端應變率下的臨界失穩應變，如圖12 所示。從圖12 中可以發現，改進的波恩準則預測的應變失穩點在B處（見圖12(b)），基于應變-應變梯度耦合模型的高階失穩準則預測的失穩點在A處。NEMD 模擬結果（見圖12(a)）表明，當應變梯度低于10?3??1時，晶體失穩首先發生在A處，與高階失穩準則的預測結果一致。在更高的應變梯度下，失穩點逐漸偏離模型預測結果，這是由于 ～B和 ～TL是在小應變梯度擾動假設下得到的，用于高應變梯度情況時誤差較大。

2.3 初始微結構對相變的影響

晶體缺陷對相變的影響是多方面的。首先，從熱力學角度看，晶體缺陷處的原子排列混亂，并且缺陷原子具有更高的遷移率，因此在缺陷處通過原子重排形成低能結構的概率比完整晶體要高。其次，從動力學角度看，缺陷引起的晶格畸變能有助于相變的形核。此外，在沖擊壓縮下，缺陷處可能優先出現塑性，如發射位錯、晶界滑移與遷移等，然后發生相變。為了深入認識晶體缺陷對高壓相變的影響規律，需要從原子尺度上系統研究缺陷對相變形核孕育時間、相變機制以及塑性與相變微觀耦合機制等的影響，MD 模擬已在相關方面取得了大量進展。

2.3.1 初始孔洞缺陷

邵建立等[68]研究了有納米孔洞（約2 nm）存在時不同沖擊強度對相變形核時間的影響（見圖13），發現形核時間與過壓程度之間滿足指數關系，與沖擊波實驗研究結果[69]一致。此外，他們還采用均勻壓縮技術進一步比較了一維和三維應變下溫度對相變形核的影響，發現一維應變條件下的形核壓力閾值明顯低于三維應變情況，并且兩者隨溫度升高的變化趨勢截然不同[70]。崔新林等[71]研究了含有納米孔洞的單晶鐵的沖擊相變，觀察到了孔洞附近沿[100]壓縮與(011)[0 11]剪切的相變機制，導致相變首先在孔洞附近形核并形成蝴蝶狀圖案（見圖14）。這些研究都是基于Voter-Chen 勢的結果，缺乏相變前的塑性行為。Wu 等[72]使用MAEAM 勢研究了含圓柱形孔洞的單晶鐵的沖擊取向效應，發現沿[110]與[111]晶向沖擊時，孔洞附近首先出現位錯環發射現象，然后發生相變，而沿[001]晶向沖擊時沒有觀察到位錯。

圖13 相變形核時間與過壓程度（ Δp）的關系：(a)有納米孔洞存在時MD 模擬結果[68]，（b）實驗結果[69]Fig. 13 Nucleation time of phase transition versus overpressurization observed in(a) MD simulations at presence of a nonvoid[68] and (b) experiments[69]

圖14 相變在納米孔洞附近的形核（紅色表示hcp 原子，黃色表示bcc 原子，其他顏色為非晶原子）[71]Fig. 14 Nucleation of phase transition around a nonvoid, where red denotes hcp atoms,yellow denotes bcc atoms and the other color denotes amorphous atoms[71]

2.3.2 晶 界

Gunkelmann 等[73]分別采用Mendelev 勢（無相變）和改進的Ackland 勢（有相變）模擬了多晶鐵的沖擊響應行為，觀察到隨沖擊強度的增大沖擊波由彈塑性兩波結構到彈性-塑性-相變三波結構的轉變，其中塑性通過形成位錯環-穿過晶界-留下螺位錯的機制進行，如圖15 所示。王昆等[67]采用MAEAM 勢模擬了不同沖擊強度下納米多晶鐵的沖擊相變行為，發現納米多晶鐵在沖擊下具有兩種截然不同的相變耦合模式，并依賴于相變-塑性/形變耦合機制。第1 種為應力援助類型的相變耦合模式，該模式在單晶鐵中也觀察到了；第2 種為應變誘導相變耦合模式。其中，應力援助相變耦合模式分兩步進行，首先沿相變不變平面的<001>晶向壓縮形成六角密排結構，然后通過密排面間的穿插形成新相，是機械能貢獻相變驅動力的過程，在MD 模擬的應變率范圍內，這種相變模式通常以晶格過壓失穩引起的無勢壘形核的形式出現，相變所需的宏觀壓力閾值較高；而應變誘導相變耦合模式中的兩個步驟在時間上是重疊的，相變形核源于晶體的局域變形和能量弛豫，即相變可以在形變過程中產生的局域應力集中處形核，形成低能結構，降低體系的局域應力，以響應外部的施加載荷，這種相變模式所需的宏觀壓力閾值較低。

圖15 在0.7 km/s 的沖擊速度下加載15 ps 時多晶鐵樣品的局部結構（圖中僅顯示缺陷原子，顏色代表原子離觀察者的距離）[73]Fig. 15 Local structure of polycrystalline iron at 15 ps under the shock with a impacting velocity of 0.7 km/s,where only defect atoms are shown (The color represents the distance from the observer.)[73]

王昆等[74]研究了初始孿晶與相變的相互作用，發現孿晶激活的相變比單晶鐵更容易發生（見圖16），但是兩者的相變變體選擇規則都滿足最大轉變功準則，這與共格孿晶界缺乏界面位錯使得界面附近的應力較均勻有關。在此基礎上，張學陽等[75]進一步研究了 Σ5和兩種類型的 Σ3（扭轉與傾斜）晶界的沖擊響應，發現相變變體選擇規則均滿足最大應變功準則，然而只有沿兩種 Σ3晶界沖擊時相變前有明顯的位錯活動（見圖17）。由于 Σ5樣品的沖擊方向為[001]，而兩種 Σ3樣品的沖擊方向分別為[110] 和[111]，晶界附近出現位錯活動的沖擊晶向與含孔洞的單晶鐵一致[72]，因此推斷在沖擊下缺陷附近位錯的活動主要與單軸應變加載的晶向有關。

圖16 含有孿晶的雙晶鐵樣品在最大粒子速度為0.4 km/s 的斜波壓縮下31 ps 時的局部構型（紅色表示hcp 原子，藍色表示bcc 原子，右圖hcp 晶胞中黃色的原子面對應于左圖中黃色方框標記的原子面）[74]Fig. 16 Local configuration of bicrystal iron sample containing a twin at 31 ps under ramp compressions with a maximum particle velocity of 0.4 km/s (Red denotes hcp atom and blue denotes bcc atom. The right figure shows a unit cell of hcp phase,where the yellow atomic face corresponds to the atomic face marked with yellow square in the left figure.)[74]

圖17 在0.5 km/s 的粒子速度下 Σ3扭轉晶界(a)與 Σ 3傾斜晶界(b)的沖擊塑性與相變（沖擊波沿z 方向傳播，在兩種樣品中分別對應于<110>和<111>晶向）[75]Fig. 17 Shock-induced plasticity and phase transition of Σ3 twist grain boundary (a) and Σ3 tilt grain boundary (b)under the shock with a particle velocity of 0.5 km/s, where shock wave propagates along z direction,corresponding to <110> and <111> crystallographic direction, respectively[75]

2.3.3 位 錯

位錯是實際晶體中一類非常重要的晶體缺陷。高壓相變在預存位錯處的形核問題非常值得關注然而相關的原子尺度研究非常少。最近，Huang 等[76]研究了含刃型位錯單晶鐵的沖擊塑性與相變行為，觀察到位錯的活躍滑移系從{110}<111>到{112}<111>的轉變。鐵的相變發生在位錯滑移系啟動之后。圖18 顯示了受塑性控制的相變機制，其中與單晶鐵的相變機制的差別在于沿[001]晶向壓縮過程中(110)面繞[111]晶向發生了一定角度的轉動，造成最終相變變體的取向與單晶鐵的相變變體存在一定差異。事實上，對于一些更高尺度的方法，如多場耦合的相場模型[77–78]，位錯與相變相互作用的微觀機制，尤其是晶體形變學上的機制，對于發展基于物理的微力學模型至關重要。目前，仍缺乏該方面的相 關基礎研究工作。，

2.4 逆相變與層裂

鐵的高壓相變是可逆的，研究相變-逆相變的規律對于理解沖擊后回收樣品中的微結構特征具有重要意義[79]，同時對金屬材料加工及后處理技術也具有重要的指導意義[10]。層裂發生在沖擊波從后表面反射卸載階段，源于反射卸載波與入射稀疏波相互作用而產生的巨大拉伸應力，當拉伸應力超過材料的抗拉強度時，材料將出現損傷和斷裂。顯然，鐵的逆相變發生在層裂行為之前，因此逆相變留下的微結構將直接影響材料的后續拉伸斷裂行為。

形變孿晶是bcc 金屬中常見的一種塑性形變機制，動態相變前產生的形變孿晶與相變的相互作用對材料微結構及其形貌有重要影響，關系著材料的宏觀動態力學行為，這也是理解含相變材料動力學行為的物理基礎。王昆等[74]從原子尺度研究了雙晶鐵中垂直和平行孿晶與相變的相互作用，發現垂直孿晶在經歷一次正相變和逆相變之后可以完全“恢復”，即存在記憶效應，然而恢復孿晶的尺寸被顯著細化了，這為加工制造超細納米孿晶材料提供了一種可能的途徑；當加載強度足夠高時，相變產物中出現滑移帶，導致恢復孿晶中出現一些二次孿晶；這些細化孿晶的尺寸分布取決于孿晶-相變的過渡勢壘。對于平行孿晶界，在相變過程中形成大量滑移帶，導致逆相變孿晶出現幾乎所有可能的變體，從而在恢復樣品中留下類似于正六邊形的特征孿晶結構。

圖18 受塑性控制的相變（bcc→hcp）過程示意圖（紅色與藍色原子分別位于相鄰的兩個原子面上）[76]Fig. 18 Schematic diagram showing the plasticity-controlled bcc→hcp phase transition,where the red and yellow atoms locate in two adjacent atom planes[76]

Gunkelmann 等[4]采用準等熵壓縮技術研究了多晶鐵的壓縮與卸載行為，觀察到了類似于實驗恢復樣品中的bcc 孿晶[79]，通過跟蹤原子運動軌跡，發現殘留孿晶的出現與逆相變過程中的剪切作用有關。在二次壓縮下，孿晶會消失，然后再發生相變。鐵相變對層裂行為的影響研究也有一些進展。通過對比無相變（Machová勢）與有相變（改進的Ackalnd 勢）的原子間相互作用勢的模擬結果，發現相變降低了多次層裂和裂紋形成的可能性[80]（見圖19）。關于鐵層裂的原子模擬工作還非常少，實驗方面卻取得了不少進展，并且非常有必要結合原子模擬方法進行解析。例如，激光沖擊實驗發現[81]，低溫下鐵的層裂具有脆性斷裂特征，而高溫下層裂面上出現了很多小面結構，表明層裂過程中有很多直裂紋的傳播，既有晶間斷裂，又有穿晶斷裂，并且分布范圍可達幾十微米（見圖20）。目前，這種斷裂模式發生轉變的原子尺度機制還不清楚。此外，沖擊波實驗還觀察到純鐵的反常“二次”層裂，與相變對沖擊波傳播的影響有關[82]，晶格層次的機理仍需進一步探索。

圖19 采用Machová勢(a)和改進的Ackland 勢(b)得到的多晶鐵層裂的MD 模擬結果[4]Fig. 19 Spallation of polycrystalline iron by MD simulations using (a) Machová potential and (b) modified Ackland potential[4]

圖20 激光沖擊實驗回收的鐵樣品：(a)樣品初始溫度293 K，層裂面處的最高加載狀態處于α→ε相變邊界上；(b) 樣品初始溫度673 K，層裂面處的最高加載狀態處于兩相共存區[81]Fig. 20 Iron samples recovered after laser shocks (a) at 293 K, where the maximum loading state of the spall plane is on the α→ε phase transition boundary, and (b) at 673 K, where the maximum loading state of the spall plane locates at the two-phase-coexistence region[81]

3 非平面沖擊載荷下的結構相變

在沖擊波加載技術中，除了易于實現、應用最廣泛的平面沖擊外，還有柱面沖擊、球面或準球面沖擊等非平面動態加載[83]。例如：武器中的彈體結構多采用圓柱殼結構，彈體爆炸過程中外殼受到的就是典型的柱面動態載荷；柱面內爆磁通量壓縮實驗裝置中的金屬套筒在工作過程中受到的柱面會聚沖擊載荷[84]；在慣性約束聚變中，球形燃料靶丸受到嚴格控制的球面會聚沖擊[85–86]。與平面沖擊相比，非平面沖擊載荷（柱面沖擊、球面沖擊或準球面沖擊）的應變條件更復雜，材料的微觀響應也明顯不同。例如：在球面和準球面會聚沖擊波下，純鐵的塑性變形主要以滑移形式發生，而球面沖擊條件下試樣內部形成的局部變形帶和微孔結構明顯多于準球面沖擊[87]。再如，多晶銅在球面會聚沖擊波作用下，在高強度沖擊中出現空位位錯環，而在低強度沖擊[88]中卻未見空位位錯環[89]。單晶銅在柱面會聚沖擊波作用下，位錯形核位點與沖擊波速度有關[90]。此外，對304 不銹鋼應變誘導馬氏體相變的研究發現，馬氏體變形體的形成強烈依賴于應力條件[91]。在相同的應變下，多軸（雙軸或三軸）加載比單軸加載產生更多的馬氏體變體。受實驗加載和檢測手段的限制，到目前為止，有關非平面沖擊載荷下材料結構、塑性和相變等響應行為的研究還比較少，Asay 等[92]在2017 年出版的專著《Impactful Times: Memories of 60 Years of Shock Wave Research at Sandia National Laboratories》的展望部分，明確提出要發展研究復雜沖擊條件下材料響應規律的實驗技術和理論方法。

對于非平面的柱面和球面沖擊來說，由于沖擊面的大小隨著時間的延長而變化，因此適用于平面沖擊的活塞加載、對稱撞擊等方法均難以直接應用于柱、球面加載中。類比平面沖擊加載模擬中采用剛性活塞以一定速度撞擊靶材產生沖擊波的基本思想，胡望宇課題組通過引入鈍化勢能面與原子間的相互作用來實現加載。在柱面沖擊模擬中，柱狀勢能面的半徑以一定速率收縮（或膨脹），實現柱面內爆（或外爆）沖擊加載（見圖21）。柱狀勢能面與原子之間的相互作用描述為[90]

圖21 柱面內爆(a)和外爆(b)沖擊加載示意圖（藍色表示試樣，紫色圓環表示柱形勢能面，紅色箭頭表示加載方向）Fig. 21 Schematic diagram of the shock loading by cylindrical(a) implosion and (b) explosion (The blue part represents the sample, the purple ring represents the energy surface of the column, and the red arrow represents the loading direction.)

式中：K為力常數，K= 1 nN/?2；θ(r?R)為標準階躍函數；R為圓柱形勢能面的半徑（內爆加載時R隨時間線性減小，外爆加載時R隨時間線性增加）；r為原子到柱狀勢能面中心軸的垂直距離。沖擊模擬采用微正則系綜（NVE）。

單晶Fe 的柱面內爆沖擊（柱面軸沿[001]晶向）模擬研究發現，在低沖擊強度下，由于晶體的各向異性特征，盡管施加的是圓柱面形會聚沖擊，沖擊波前的形狀由開始的圓柱面逐漸轉變為方柱面（見圖22）。同平面沖擊相比，材料在柱面沖擊下受到的應變約束完全不同，材料內部單元之間存在相對運動，沖擊壓縮區的應力和溫度分布要復雜得多。在沿<110>方向的塑性變形區域出現了典型的1/2<111>位錯環，相變在位錯環區域形核。沖擊相變過程中形成了大量的相變孿晶（見圖23）。

圖22 （a）柱面軸沿單晶Fe [001]晶向在內爆沖擊下垂直于柱面軸的橫截面內的原子速率分布，（b）10 ps 時垂直柱面軸橫截面內的局域溫度分布Fig. 22 (a) Atom velocity distribution in the cross section perpendicular to the cylindrical axis under implosion impact along [001] direction of Fe single crystal, and (b) the local temperature distribution in the cross section perpendicular to the cylindrical axis at 10 ps

圖23 柱面軸沿單晶Fe [001]晶向在內爆沖擊下（vp = 0.6 km/s）局域微結構隨時間的演變（紅色表示hcp 結構原子，綠色表示以層錯形式出現的fcc 結構原子，灰色表示位錯線周圍的其他結構類型的原子，為清晰起見，未相變的bcc 結構原子已被刪除）Fig. 23 Evolutions of local microstructures with time under the implosive impacting with the cylindrical axis along [001] direction of Fe single crystal (vp = 0.6 km/s) (Red indicates hcp atoms,green indicates fcc atoms in the forms of stacking faults, and gray indicates other atoms distributed around dislocation lines. For clarity, bcc atoms without phase transition have been removed.)

對比柱面內爆形成的會聚沖擊，外爆沖擊過程中沖擊波前隨時間的演化而增大（見圖24），波前的沖擊強度明顯降低。在<110>的4 個等價方向上，沒有觀察到位錯環的出現，但在相變波前有類似于沿[110]晶向平面沖擊時的位錯運動。柱面外爆時，在持續的沖擊載荷下，可以觀察到bcc-hcp-bcc-hcp 循環相變，hcp-bcc 的逆相變導致bcc 相相對于初始時出現了45°的重定向。在鐵的bcc-hcp 相變中，剪切應力驅動原子在{110}滑移面中移動。在平面沖擊載荷下，當{110}滑移面平行于沖擊方向時，在單軸應變形成的強側向約束作用下，會產生垂直于沖擊載荷的剪應力。相反，柱面發散沖擊載荷下，應變是多軸的，導致產生多個垂直于應變主軸的剪應力，使hcp 變體c軸取向多樣化，不再局限于典型的<110>方向（見圖25）。

圖24 柱面外爆沖擊后試樣內局域微結構、速度波剖面、應力波剖面和溫度分布Fig. 24 Local microstructure, velocity profile, stress profile and temperature distribution in the sample after cylindrical explosions

圖25 柱面外爆沖擊后試樣內部的局域微結構(a)以及hcp 相變變體c 軸取向極圖(b)Fig. 25 Local microstructure inside the sample (a) and polar figure of the c-axis showing different hcp variants (b)

4 結論與展望

本文綜述了動態加載下鐵相變的原子模擬研究現狀，尤其是近年來隨著新的原子間相互作用勢被廣泛應用于鐵高壓相變的原子模擬研究，鐵在各種不同加載條件（包括沖擊晶向、沖擊強度、應變率等）與初始微結構（如初始位錯、孿晶界及一般晶界等）下的相變機理不斷被發現，已初步形成加載條件-相變晶體學機制-變體類型選擇-產物織構間的定性認識。得益于鐵的勢函數的改進，高壓相變與塑性微觀耦合機制研究已成為該領域中新的研究方向，對于一些基于物理的微力學建模具有重要的指導意義，同時也是理解含相變金屬材料動態力學性能不可或缺的一部分。已有的研究表明，相變前的塑性與加載晶向及缺陷類型密切相關，對后續相變機制與形核有重要影響。然而，目前還未形成完整的認識，仍需針對不同類型的晶體缺陷進行系統的模擬研究。介觀尺度的方法在描述晶體的高壓相變與塑性方面取得了長足的進展，如基于反應路徑的相場模型[93–95]已具備較強的動態力學響應行為預測能力，基于經典密度泛函理論的晶體相場模型[96–97]可以直接獲取介觀時間尺度的相變與塑性的原子機制，這些模型與MD 模擬方法一起，能夠完全覆蓋材料在微介觀時空尺度上從原子分辨到連續尺度的動態力學響應。此外，相變與層裂也是材料動力學行為研究的一個重要方向。原子模擬在這方面的研究非常少，其中一種可能的原因是缺乏相關晶格尺度原位診斷實驗的支撐。

盡管原子模擬在鐵動態相變研究中已經取得了豐碩的成果，但是還存在至少如下幾方面的不足，甚至矛盾。其一，現有的兩種可以正確描述鐵 α?ε相變的勢函數在描述沖擊塑性方面存在一些明顯的差異。改進的Ackland 勢可以預測出單晶和多晶鐵的變形孿晶，但只在多晶鐵中觀察到位錯，而MAEAM勢可以預測單晶鐵的沖擊塑性，但沒有出現相變前的形變孿晶，澄清這些問題需要更多晶格層次的原位實驗觀測。其二，目前的原子模擬研究主要以基于晶體學的機理認識為主，對相變動力學的理解仍然以經典的各向同性的形核理論為依據，在理論上幾乎沒有實質性的發展。此外，需要系統地研究材料的動態力學行為與原子間相互作用勢的關系，從勢函數出發理解材料動態力學行為機制發生的原因無疑極其重要：一方面，有助于進一步改進原子間相互勢函數，指導開發其他材料的高壓原子勢函數，助力材料基因工程的發展；另一方面，有助于從熱力學上整體把握材料力學行為與體系基本物理性質之間的關系，這對于發展其他基于能量泛函的更高尺度模型（如相場模型等）具有直接的理論指導價值。