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

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

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

1 研究進展

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

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

絕熱狀態(tài)方程通?；贖ugoniot數(shù)據(jù)得到，而Hugoniot數(shù)據(jù)則來源于通過輕氣炮或者爆炸載荷驅動的飛片撞擊實驗。Al’tshuler等［11］研究了高壓區(qū)Al、Cu和Pb的狀態(tài)方程。Zhang等［12］研究了93 W合金的Hugoniot曲線、沖擊溫度和多項式狀態(tài)方程。Dewale等［13］研究了Ta、Au、Pt、Al、Cu和W在94 GPa以上的靜態(tài)和動態(tài)狀態(tài)方程。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狀態(tài)方程作為高壓下的固體模型得到了廣泛應用，許多仿真軟件可以通過直接輸入Mie-Grüneisen狀態(tài)方程的參數(shù)來定義材料。沈天樂等［19］通過實驗得到了3種不同密度的環(huán)氧樹脂－SiO2空心微球復合材料的Hugoniot關系，討論了基于Hugoniot關系得到的Mie-Grüneisen狀態(tài)方程，并與P-α模型進行了比較，顯示結果一致。Walsh等［20－21］通過爆炸系統(tǒng)測量了金屬的狀態(tài)方程，并估計了Grüneisen系數(shù)γ的體積依賴性。林華令等［22］給出了一種用沖擊壓縮數(shù)據(jù)計算Grüneisen物態(tài)方程參數(shù)的方法，計算了23種金屬的Grüneisen系數(shù)，使Grüneisen系數(shù)沿實驗沖擊絕熱線的計算值與晶格振動理論的計算值達到最優(yōu)擬合。Mc-Queen等［23］通過平面波爆炸系統(tǒng)在100～200 GPa壓力下獲得了19種金屬元素的狀態(tài)方程。通過Mie-Grüneisen理論，Hugoniot關系（P，V，E）被擴展為更完整的狀態(tài)方程（P，V，E，T）。吳紹曾等［24］利用剛性離子模型計算NaCl的Grüneisen參數(shù)及其隨壓力的變化，計算中對最近鄰排斥勢分別使用兩種修正Born-Mayer勢，并利用電子氣理論的計算結果處理次近鄰排斥勢和關聯(lián)勢的影響，最后得到的結果與實驗相符。Nellis等［25］通過80～440 GPa壓力范圍內的Hugoniot曲線測量Al、Cu和Ta的狀態(tài)方程，并討論了沖擊溫度和γ。Jin等［26］提出了基于絕熱沖擊線數(shù)據(jù)表征金屬的0 K等溫線的熱力學公式，并使用Born-Meyer勢預測了7種金屬的0 K等溫線和Grüneisen狀態(tài)方程。

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

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

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

2 理論框架

2.1 復雜合金的C0和λ

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

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

2.2 Hugoniot曲線及其內能

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

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

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

其中：C0是壓力為0時的體積聲速；λ是材料常數(shù)。這些量是通過擬合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范圍內進行溫度修正時，體積膨脹系數(shù)αv可以近似地看作一個常數(shù)。V0K、ρ0K、C′0和λ′可以寫成：

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

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

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

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

其中，A和B是材料常數(shù)。

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

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

2.4 Grüneisen狀態(tài)方程

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

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

Grüneisen系數(shù)量綱為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）中的參數(shù)，可以得到：

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

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

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

3 計算方法

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

所有模型尺寸均為20a×20a×20a，a為單個晶胞的晶格常數(shù)。在3個方向上均使用周期性邊界條件。在加載沖擊前，首先在NPT系統(tǒng)中對模型在300 K進行40 ps的弛豫，使其達到可以忽略系統(tǒng)中殘余應力的狀態(tài)。弛豫后通過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，F(xiàn)e的波速范圍為6～10 km／s，Mo的波速范圍為6～10 km／s。當模擬持續(xù)足夠長的時間，粒子速度、壓強、溫度和其他物理量均達到穩(wěn)定后，從穩(wěn)定狀態(tài)下的數(shù)據(jù)中提取需要的數(shù)據(jù)。

4 結果與討論

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

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

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

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

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

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

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

4.2 Hugoniot曲線及其內能

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

圖1 GH4169合金的Hugoniot曲線

圖2 GH4169合金的內能

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

4.3 0 K等溫線及其內能

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

圖3、4為GH4169合金的0 K等溫線Pc及其內能Ec。藍色方塊為使用文獻［41］方法得到的數(shù)據(jù)計算結果，黑線和紅線分別表示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狀態(tài)方程

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

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

圖5 P-V-E形式的Grüneisen狀態(tài)方程曲面，以式（40）作為Grüneisen系數(shù)

圖6 P-V-E形式的Grüneisen狀態(tài)方程曲面，以式（41）作為Grüneisen系數(shù)

5 結論

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

2）GH4169鎳基高溫合金的Born-Meyer勢材料常數(shù)Q、q和Morse勢材料常數(shù)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狀態(tài)方程在P-V-E空間中表現(xiàn)為一個凹面，Grüneisen系數(shù)的模型是影響狀態(tài)方程最大壓力值的重要因素。