基于分子動力學方法的復雜合金Grüneisen狀態方程參數計算及驗證：以GH4169合金為例

張圣來，丁 軍，童 權，黃 霞，宋 鹍，路世青

（重慶理工大學 機械工程學院，重慶 400054）

1 研究進展

當前，航天航空等領域對性能優異的復雜合金的需求越來越大，對其使用環境要求也更加嚴格。例如，航空發動機在工作中一直處于高溫高壓狀態，因此選擇材料時必須重點考慮極端服役環境下的合金性能。

在超高壓下，應力偏量幾乎可以忽略，此時材料的本構模型退化為體積應力關系，即高壓物理領域中的狀態方程。Mie-Grüneisen狀態方程作為應用最廣泛的一種狀態方程，通常用于描述物質在極端條件下的狀態，如沖擊波傳播。在這些極端條件下，物質會處于高溫高壓狀態，性質發生變化。近幾十年來，等溫和絕熱狀態方程研究取得了很大進展。等溫狀態方程有兩種：一種是在0 K，也叫作冷曲線；另一種是在300 K。0 K等溫狀態方程一般根據300 K下的沖擊波數據通過溫度修正得到。300 K等溫狀態方程是靜態狀態方程，通過靜態壓縮保持溫度不變。Bridgman［1］開發了一種靜態壓縮裝置，用于探索材料在高壓下的壓縮特性。由于壓縮過程非常緩慢，材料內部產生的熱量能與外界完全交換，因此可將該壓縮過程認為是等溫壓縮。在此基礎上，許多學者對材料進行高壓等溫壓縮。張劍等［2］采用原位高壓同步輻射X射線衍射技術在室溫、21.5 GPa的壓力范圍測定鋁的等溫狀態方程，利用Murnaghan方程對實驗結果進行數值擬合，得到鋁的零壓體彈模量及其1級壓力導數。該結果同相關文獻資料研究結果在誤差范圍內符合得很好。Syassen等［3］采用12 GPa的高壓X射線衍射測定Ag和Al的等溫壓縮曲線，并選擇Birch方程作為等溫狀態方程來擬合實驗數據。Voˇcadlo等［4］和Shanker等［5］通過比較各種等溫狀態方程，研究了等溫體積模量及其1階壓力導數。Wang等［6］使用經典平均場方法計算了壓力高達1 TPa時5種金屬Al、Cu、Ta、Mo和W的293 K等溫線。Hrubiak等［7］采用雙參數Birch-Murnaghan方程描述了Hf在室溫下的等溫線，得到了67 GPa壓力下α相Hf的P（E，V）型狀態方程。Wang等［8］用平均場勢法對Al的靜態和動態壓縮及體積模量進行全面研究，檢驗了Hugoniot曲線和300 K等溫線。Akahama等［9］通過X射線衍射實驗，利用Birch-Murnaghan和Vinet狀態方程描繪了Bi、Pt和Au的300 K等溫線。Verma等［10］計算了Al、Ta、Mo和W的固態和液態Hugoniot曲線，討論了0 K等溫線和P－T曲線。

絕熱狀態方程通常基于Hugoniot數據得到，而Hugoniot數據則來源于通過輕氣炮或者爆炸載荷驅動的飛片撞擊實驗。Al’tshuler等［11］研究了高壓區Al、Cu和Pb的狀態方程。Zhang等［12］研究了93 W合金的Hugoniot曲線、沖擊溫度和多項式狀態方程。Dewale等［13］研究了Ta、Au、Pt、Al、Cu和W在94 GPa以上的靜態和動態狀態方程。Bringa等［14］和Agarwal等［15］對Mg、Cu單晶的絕熱沖擊線進行非平衡分子動力學（NEMD）研究，得到了不同晶向下沖擊波速度與粒子速度的線性關系。Mackenchery等［16］通過分子動力學（MD）討論了4種不同原子間勢對單晶Ti在0.5～2.0 km／s碰撞速度下的線性關系和Hugoniot曲線的影響。

Mie［17］和Grüneisen［18］提出了描述壓力和體積之間關系的固體理論。目前，Mie-Grüneisen狀態方程作為高壓下的固體模型得到了廣泛應用，許多仿真軟件可以通過直接輸入Mie-Grüneisen狀態方程的參數來定義材料。沈天樂等［19］通過實驗得到了3種不同密度的環氧樹脂－SiO2空心微球復合材料的Hugoniot關系，討論了基于Hugoniot關系得到的Mie-Grüneisen狀態方程，并與P-α模型進行了比較，顯示結果一致。Walsh等［20－21］通過爆炸系統測量了金屬的狀態方程，并估計了Grüneisen系數γ的體積依賴性。林華令等［22］給出了一種用沖擊壓縮數據計算Grüneisen物態方程參數的方法，計算了23種金屬的Grüneisen系數，使Grüneisen系數沿實驗沖擊絕熱線的計算值與晶格振動理論的計算值達到最優擬合。Mc-Queen等［23］通過平面波爆炸系統在100～200 GPa壓力下獲得了19種金屬元素的狀態方程。通過Mie-Grüneisen理論，Hugoniot關系（P，V，E）被擴展為更完整的狀態方程（P，V，E，T）。吳紹曾等［24］利用剛性離子模型計算NaCl的Grüneisen參數及其隨壓力的變化，計算中對最近鄰排斥勢分別使用兩種修正Born-Mayer勢，并利用電子氣理論的計算結果處理次近鄰排斥勢和關聯勢的影響，最后得到的結果與實驗相符。Nellis等［25］通過80～440 GPa壓力范圍內的Hugoniot曲線測量Al、Cu和Ta的狀態方程，并討論了沖擊溫度和γ。Jin等［26］提出了基于絕熱沖擊線數據表征金屬的0 K等溫線的熱力學公式，并使用Born-Meyer勢預測了7種金屬的0 K等溫線和Grüneisen狀態方程。

從上述研究中不難發現，Mie-Grüneisen狀態方程可由等溫線或Hugoniot曲線推導得到，因而兩者均可作為參考線。雖然已推導出了Mie-Grüneisen狀態方程的公式，但很難通過實驗獲得狀態方程的所有參數。確定這些參數的手段除實驗外還必須借助理論分析和分子動力學模擬方法。Yang等［27］使用分子動力學方法全面研究了Hugoniot曲線、等溫線、等熵線及其內能，通過分子動力學模擬獲得了單晶Al和Pb的Mie-Grüneisen狀態方程，計算并討論了沿Hugoniot曲線和等熵線的溫度、沿Hugoniot曲線的熵增量和Grüneisen系數γ，研究結果與已有文獻吻合較好。Chen等［28］通過分子動力學模擬研究了TiAl合金的Hugoniot曲線和Mie-Grüneisen狀態方程，發現材料成分嚴重影響沖擊波速與質點速度的線性關系、Hugoniot曲線和內能。Zhang等［29］使用分子動力學研究單晶鎳的沖擊行為。模擬得到的P－V曲線略高于實驗數據，實驗樣品為多晶Ni且內含缺陷，因此導致了數據差異。Zhang等［30］提出了一種新的與溫度相關的能熱力學模型，用以研究材料的溫度相關Grüneisen。通過分子動力學對Ni、Cu、Au進行分子動力學模擬，驗證了新建立的熱力學模型。

目前對于等溫線、Hugoniot曲線和Mie-Grüneisen狀態方程的研究較多，材料類型主要集中在單一元素和二元合金，對多元素的復雜合金的研究相對較少。眾所周知，復雜合金的物理性質與其組成元素的物理性質有很大不同，且服役環境通常較復雜。因此，有必要研究復雜合金在極端條件下的動態力學性能。GH4169作為一種鎳基高溫合金，是航空發動機葉片的主要材料。該材料在高溫下仍具有很好的力學性能，因此一直被研究人員所關注。

然而，由于缺少相應的勢函數，分子動力學目前無法通過模擬直接得到復雜合金Mie-Grüneisen狀態方程的各個參數，因此在Ji等［31］和Millett等［32］的工作基礎上，提出一種用于計算復雜合金的Mie-Grüneisen狀態方程參數，并以GH4169鎳基合金為例計算了Mie-Grüneisen狀態方程。

2 理論框架

2.1 復雜合金的C0和λ

由于目前已有的勢函數中，包含4種元素以上的勢函數較少，因此使用單一的勢函數無法直接計算出復雜合金us和up的關系。Ji等［31］在對累計滾焊Al／Ni復合材料的沖擊誘導反應特性進行研究時，通過純Al、純Ni的參數C0和λ，以質量占比為權重，擬合Al／Ni復合材料的C0和λ，從而得到Al／Ni復合材料的Grüneisen狀態方程。結果證明理論結果與實驗結果符合良好，驗證了所提方法的正確性。Millett等［32］通過平板撞擊實驗得到TiAl合金Ti-46.5Al-2Nb-2Cr的Hugoniot數據，并通過合金中各元素的參數C0和λ，以質量占比為權重擬合了Ti-46.5Al-2Nb-2Cr合金的C0和λ。結果表明，通過這種計算方式獲得的Hugoniot參數與實驗測得的參數誤差較大。鑒于結論存在矛盾，本文中提出一種新的思路：以質量分數為權重，通過合金中各穩定相的Hugoniot參數，對合金的Hugoniot參數進行擬合：

其中：mPi代表第i種相的質量分數；下標“A”和“Pi”分別代表合金和第i種相；n是合金中存在的相的總數。

2.2 Hugoniot曲線及其內能

Hugoniot曲線由質量、動量和能量守恒方程決定，其表達式為：

其中：ρ為密度；V＝1／ρ為受沖擊材料的比體積；up為粒子速度；us為沖擊波波速；P為壓力；E為比內能。下標“0”和“H”分別表示這些量處于未受沖擊的初始態和Hugoniot狀態。

常溫（300 K）下，在某個很寬泛的壓力范圍內，us和up的關系是線性的：

其中：C0是壓力為0時的體積聲速；λ是材料常數。這些量是通過擬合us和up得到的。

由式（3）～（6）可以推導出Hugoniot曲線PH（V）和內能EH（V）：

其中：cV是定容比熱；T是沖擊溫度。當EH遠大于E0時，可以忽略E0。

2.3 冷壓和冷能

0 K下，us0K與up0K之間線性關系和Hugoniot曲線PH0K（V）的表達式為：

通常，式（6）是在常溫（300 K）下通過實驗得到，因此式（10）中的C′0和λ′需要基于C0和λ進行溫度修正。

在0～300 K范圍內進行溫度修正時，體積膨脹系數αv可以近似地看作一個常數。V0K、ρ0K、C′0和λ′可以寫成：

其中：γ（V0）是300 K時的Gruneisen系數。從式（11）～（15）可以算出C′0和λ′，此時可以求得式（10）的Hugoniot曲線PH0K。

Born-Meyer勢可以精準描述離子晶體和金屬的可壓縮性，因此常用于計算冷壓Pc－BM和冷能Ec－BM，其表達式為：

其中：Q和q是材料常數；ρ0K是0K時的材料密度；δ＝V0K／V。

Morse勢同樣可以描述高壓0 K條件下原子間的相互作用，因此Morse勢的冷壓Pc－M和冷能Ec－M可以表達為：

其中，A和B是材料常數。

一般，假設PH0K、等熵線PS（V）和Pc（V）在初始點（P＝0，V＝V0K）的1階導和2階導近似相等：

由式（8）～（21）可得材料常數Q、q、A和B的表達式為：

2.4 Grüneisen狀態方程

Grüneisen狀態方程是高壓物理和爆炸力學領域描述固體行為最常用的狀態方程之一。作為特征方程，它描述了晶格熱振動的貢獻。Grüneisen狀態方程最常用的參考曲線是Hugoniot曲線，其經典形式為：

經典形式以冷壓Pc和冷能Ec作為參考線，但常溫下獲得的Hugoniot曲線PH（V）及其內能EH（V）同樣可以作為參考曲線。替換參考曲線后的狀態方程為：

Grüneisen系數量綱為1，表達式為［33］：

其中：t代表表征γ的不同模型，當t＝0、1、2，式（28）分別表示Slate模型γS，Dugdale-MacDonald模型γDM和自由體積模型γf。

PS在初始點（P＝0，V＝V0）近似滿足式（28），因此式（28）也可以寫成：

當t＝0、1和2時，有：

假設PS與PH的1階導和2階導在初始點（P＝0，V＝V0）處恒等，而P′H（V0）和P″H（V0）又可直接由實驗或者分子動力學模擬得到的PH求導得到，那么當獲得式（2）中的參數，可以得到：

將式（33）～（35）代入式（29）可得：

Grüneisen系數γ僅是V的函數，且γ＞0。γ的表達式一般是通過實驗數據擬合得到的經驗公式，而由于本文中研究在無實驗條件下僅通過數值模擬得到復雜合金的Grüneisen狀態方程，因此可用的經驗公式為：

基于式（1）～（2）、（6）～（8）、（17）和（40）～（41），最終可以得到復雜合金的Grüneisen狀態方程。

3 計算方法

采用的勢函數有Mishin等［34］開發的Ni的嵌入原子法（EAM）勢，Zhang等［35］開發的Ni-Nb的嵌入原子法（EAM）勢，Ko等［36］開發的Ni-Ti的修正嵌入原子法（MEAM）勢，Purja等［37］開發的Ni-Al的嵌入原子法（EAM）勢，Howells等［38］開發的Cr的角相關勢（ADP），Chamati等［39］開發的Fe的嵌入原子法勢和Zhou等［40］開發的嵌入法原子勢。

所有模型尺寸均為20a×20a×20a，a為單個晶胞的晶格常數。在3個方向上均使用周期性邊界條件。在加載沖擊前，首先在NPT系統中對模型在300 K進行40 ps的弛豫，使其達到可以忽略系統中殘余應力的狀態。弛豫后通過MSST在X軸方向對模型施加沖擊波，Ni的波速范圍為11～15 km／s，Ni3Nb的波速范圍為10～18 km／s，Ni3Ti的波速范圍為6～10 km／s，Ni3Al的波速范圍為6～10 km／s，Cr的波速范圍為11～15 km／s，Fe的波速范圍為6～10 km／s，Mo的波速范圍為6～10 km／s。當模擬持續足夠長的時間，粒子速度、壓強、溫度和其他物理量均達到穩定后，從穩定狀態下的數據中提取需要的數據。

4 結果與討論

4.1 GH4169鎳基高溫合金的C0和λ

GH4169鎳基高溫合金的化學成分復雜，且許多微量金屬元素的質量分數小于1%，因此對GH4169的成分進行簡化：首先，保留質量分數大于1%的金屬元素，滿足此條件的元素有Ni、Cr、Fe、Nb、Mo和Ti；其次，由于GH4169鎳基合金中具有γ′相，對增強GH4169的高溫性能具有重要意義，因此所有γ′相包含的金屬元素也要保留，滿足此條件的元素有Ni、Nb、Ti和Al。經簡化后的GH4169鎳基高溫合金各主要元素的質量分數如表1所示。該成分來源由上海康晟航材的GH4169高溫合金成分表簡化而來，對所選元素中Ni以外的元素，基于成分的含量范圍進行質量分數取整后，余量即為Ni元素的質量分數。

表1 GH4169鎳基合金各元素的質量分數 和原子分數 %

GH4169高溫鎳基合金中，γ″相為Ni3Nb，是主要強化相；γ′相包括Ni3Al和Ni3Ti，為輔助強化相；Ni是γ相，為基體。所涉及的相中無δ相，根據上海康晟航材給出的熱處理制度，無δ相GH4169的熱處理方式為（1 010～1 065）℃±10℃，1 h，油冷、空冷或水冷＋720℃±5℃，8 h，以50℃／h爐冷至620℃±5℃，8 h，空冷。以原子占比為依據進行計算，可得各相的理想質量分數如表2所示。此外，通過分子動力學的多尺度模擬技術（MSST）獲得的各相的C0和λ如表2中所示。

表2 GH4169鎳基合金中各相的質量分數 和Hugoniot參數

表3為通過擬合得到的GH4169鎳基高溫合金的C0和λ。通過與文獻［41］的研究成果進行對比發現，所提出的理論計算得到的C0和λ與通過一級輕氣炮進行平板撞擊實驗得到的C0和λ非常接近，證明所提出的計算理論是正確的。

表3 GH4169鎳基高溫合金的Hugoniot參數 和相對誤差

4.2 Hugoniot曲線及其內能

得到C0和λ后代入式（7）（8）可以得到GH4169合金的Hugoniot曲線及其內能曲線。

圖1 GH4169合金的Hugoniot曲線

圖2 GH4169合金的內能

紅色方塊為使用文獻［41］方法得到的數據計算結果，可以發現該結果與本文方法計算結果一致，表明本文中所提出方法的計算結果可靠。

4.3 0 K等溫線及其內能

通過分子動力學模擬升溫過程，將獲得的體積數據代入式（11）（12），得到GH4169合金的體積膨脹系數αv為251×10－6K－1。通過式（13）～（15），（22）～（25）可以計算得到Born-Meyer勢和Morse勢的0 K等溫線及其內能曲線。材料常數Q、q、A和B分別為70.55 GPa、11.55、151.52 GPa和4.45。

圖3、4為GH4169合金的0 K等溫線Pc及其內能Ec。藍色方塊為使用文獻［41］方法得到的數據計算結果，黑線和紅線分別表示Born-Meyer勢和Morse勢，這兩種勢均可以用來描述0 K等溫線。

圖3 GH4169合金的冷壓曲線

圖4 GH4169合金的冷能曲線

可以看到，0 K等溫線Pc及其內能Ec都隨著V0K／V的增大而增大，且V0K／V越大，壓力Pc和內能Ec的增長速率越大，其冷曲線特征與Cu和Pb一致［42］。當V0K／V小于1.5時，Born-Meyer勢與Morse勢的等溫線、內能曲線幾乎是重合在一起的，且與劉焦等［41］的結果吻合良好；當V0K／V大于1.5時，Born-Meyer勢與Morse勢的等溫線、內能曲線逐漸分離，此時不論是等溫線還是內能，Morse勢都更接近于文獻［41］的結果。這表明當V0K／V小于1.5時，Born-Meyer勢和Morse勢都是比較準確的，而當V0K／V大于1.5時，Morse勢更加準確。

4.4 Grüneisen狀態方程

通過式（1）～（2），（6）～（8），（27）和（40）～（41），最終可以得到Grüneisen狀態方程的3D曲面圖，如圖5、6所示。

圖5、6分別是以式（40）、（41）作為Grüneisen系數獲得的狀態方程。雖然從正面看，2個Grüneisen狀態方程的形狀均為一個凹面，且形狀相似，但可以發現最大壓力是不相同的。這是由Grüneisen系數方程不同造成的，因此可以推斷，Grüneisen系數是影響Grüneisen狀態方程壓力值的重要因素。

圖5 P-V-E形式的Grüneisen狀態方程曲面，以式（40）作為Grüneisen系數

圖6 P-V-E形式的Grüneisen狀態方程曲面，以式（41）作為Grüneisen系數

5 結論

1）使用所提出的計算復雜合金Hugoniot參數的理論對GH4169鎳基高溫合金的C0和λ進行計算，結果與實驗結果吻合度良好，證明該理論可行。

2）GH4169鎳基高溫合金的Born-Meyer勢材料常數Q、q和Morse勢材料常數A、B分別為70.55、11.55、151.52 GPa和4.45。當V0K／V小于1.5時，Born-Meyer勢和Morse勢都可以準確地描述0 K等溫線，而當V0K／V大于1.5時，Morse勢對于0 K冷曲線的描述更加準確。

3）Grüneisen狀態方程在P-V-E空間中表現為一個凹面，Grüneisen系數的模型是影響狀態方程最大壓力值的重要因素。