劉超，石藝娜，秦承森，梁仙紅

(北京應(yīng)用物理與計算數(shù)學(xué)研究所，北京100088)

0 引言

沖擊相變是當(dāng)前沖擊波與爆轟物理研究領(lǐng)域的熱點(diǎn)問題之一，由于其對于高溫高壓極端條件下材料物性研究的科學(xué)意義及在高新技術(shù)領(lǐng)域的應(yīng)用背景，沖擊相變研究受到廣泛關(guān)注。同時，沖擊相變還是一個非線性、非平衡的復(fù)雜物理過程，具有時間相關(guān)效應(yīng);并且沖擊相變研究既是一個涉及物理、力學(xué)、材料科學(xué)及冶金學(xué)的交叉學(xué)科問題，又是一個典型的多個尺度問題。因此，沖擊相變也是沖擊波物理研究的難點(diǎn)問題之一。自1956 年Bancroft 等采用炸藥加載首次觀察到鐵的沖擊相變［1］以來，眾多研究者在沖擊相變領(lǐng)域積累了豐富的文獻(xiàn)資料，同時也發(fā)現(xiàn)了許多有趣的實(shí)驗(yàn)現(xiàn)象。

時至今日，沖擊相變研究仍集中于實(shí)驗(yàn)和唯象理論模型方面，對于其細(xì)觀機(jī)理研究相對較少，其中一個很重要的原因是缺乏細(xì)觀尺度數(shù)值模擬手段。盡管多尺度算法研究發(fā)展十分迅猛，但多集中于宏觀尺度計算方法與分子動力學(xué)等微觀尺度算法耦合的搭橋類算法，細(xì)觀尺度的計算方法很少。目前，研究多晶材料沖擊響應(yīng)的細(xì)觀尺度數(shù)值模擬方法，原則上可分為兩類:一類是以離散元為代表的粒子類方法，另一類是以連續(xù)介質(zhì)力學(xué)為基礎(chǔ)的傳統(tǒng)數(shù)值模擬方法。與傳統(tǒng)數(shù)值模擬方法相比，離散元具有模型構(gòu)建方法易行、晶粒間各相異性特性表征便捷、算法實(shí)現(xiàn)簡單等優(yōu)點(diǎn)。

離散元法是20 世紀(jì)70 年代初由美國學(xué)者cundall 首先提出的［2］，最初主要應(yīng)用于巖石力學(xué)、顆粒態(tài)群體及土壤力學(xué)問題分析中。離散元方法允許單元間的相對運(yùn)動，而且不象連續(xù)介質(zhì)模型那樣的依賴于高度簡化、規(guī)定性的本構(gòu)關(guān)系，并具有算法簡單、易于實(shí)現(xiàn)的優(yōu)點(diǎn)。20 世紀(jì)90 年代初，Sawamoto等［3］首先將離散元方法成功地用于混凝土動態(tài)沖擊破壞等非線性大變形問題的數(shù)值模擬研究。劉凱欣等在這一領(lǐng)域做了大量的研究工作［4-5］。1999年以來，Yano 等［6］、Case 等7］利用離散元方法研究了銅、鐵等多晶金屬的沖擊響應(yīng)。2000 年以后，王文強(qiáng)［8］、于繼東等［9］將離散元方法應(yīng)用到非均質(zhì)材料炸藥在沖擊作用下細(xì)觀損傷的研究。這些工作展示了離散元方法模擬細(xì)觀非均質(zhì)材料動力學(xué)問題的能力。由于可以方便地表征晶粒間取向的分布特征，離散元方法在模擬細(xì)觀非均質(zhì)材料的沖擊響應(yīng)方面具有獨(dú)特的優(yōu)勢。

本文利用離散元方法結(jié)合基于無擴(kuò)散相變的兩相模型，對于α 鐵的沖擊相變過程進(jìn)行了數(shù)值模擬研究，給出了鐵的相邊界及沖擊Hugoniot 關(guān)系。并在此基礎(chǔ)上對于多晶鐵的沖擊相變過程進(jìn)行了模擬，研究了晶粒大小及傳播距離對于沖擊波前沿不規(guī)則度的影響，并對沖擊相變過程中相變特征量隨局部沖擊壓力的變化規(guī)律進(jìn)行了統(tǒng)計。

1 計算方法與相變模型

離散元方法通過求解多體運(yùn)動的牛頓力學(xué)方程組，跟蹤全部單元的運(yùn)動軌跡，來揭示系統(tǒng)與外界的相互作用和自身響應(yīng)、演化規(guī)律。與傳統(tǒng)的數(shù)值模擬方法相比，離散元具有處理沖擊載荷作用下單元間常見的大變形、斷裂等問題方便，算法實(shí)現(xiàn)簡單等優(yōu)點(diǎn)。

1.1 單元間相互作用力模型

一般單元間可能包括以下相互作用力［10］:1)中心勢力;2)中心阻尼;3)彈塑性剪切力;4)切向阻尼;5)干摩擦力。如圖1 所示，圖中v 為線速度，ω為角速度。

圖1 單元間相互作用力的示意圖Fig.1 Interaction model of an element-pair

利用圖1 描述的單元間相互作用力模型，可將單元i 與j 間的作用力合力表示為

根據(jù)所研究的問題選取合適的單元間作用力模型是離散元方法研究的核心內(nèi)容，本文采用適合于描述材料沖擊響應(yīng)的單元間作用力模型，模型中單元間相互作用力主要包括中心勢力與中心阻尼，兩相單元間相互作用力模型參數(shù)取值見文獻(xiàn)［11］。

1.2 溫度與相變模型

影響單元溫度的力學(xué)過程可以分為可逆過程與不可逆過程，可逆的力學(xué)過程如中心勢力的作用過程，不可逆的力學(xué)過程為耗散力的作用過程［11］:

采取與Forbes［12］類似的推導(dǎo)方法，可以得到溫度的可逆部分，即等熵過程中溫度可表示為

式中:T0為初始溫度;p 為靜水壓力;B0與等溫體模量KT相關(guān)，B0= KT/β，β 為狀態(tài)方程參數(shù);γ0為Gruneisen 參數(shù)。

不可逆部分僅考慮了由熱傳導(dǎo)和粘性力帶來的能量耗散過程，不可逆過程帶來的溫升是上述兩類過程的累計效應(yīng):

連接與接觸單元間考慮了基于傅里葉定律的熱傳導(dǎo)過程:

式中:ΔQ 為在Δt 時間內(nèi)由單元j 向單元i 傳導(dǎo)的熱;κ 為熱傳導(dǎo)系數(shù);Ti和Tj分別為單元i 與單元j的溫度;dij為單元i 與單元j 間的距離;Aij為單元i與單元j 間的接觸面積。

由熱傳導(dǎo)導(dǎo)致的單元溫度變化，可表示為

由粘性力帶來的溫升:

式中:m 為單元質(zhì)量;cv為等容比熱;Cn為粘性系數(shù)為單元i 與單元j 間的相對速度在其中心連線方向的投影。

本文采用基于無擴(kuò)散相變的兩相模型，該模型基于如下假設(shè):1)每個單元中的相變均勻產(chǎn)生，并由局部能量條件控制;2)單元中的壓力和溫度滿足局部平衡條件。相變模型包括3 個重要的組成部分:1)平衡的相邊界;2)局部閾值條件;3)相變動力學(xué)方程。

首先，平衡的相邊界，即相邊界處兩相的Gibbs自由能相等Gα(p，T)= Gε(p，T)). 此處i 相的Gibbs 自由能的表達(dá)式為

式中:Hi和Si分別為比焓與熵，下標(biāo)i 分別表示α與ε 相。

在壓力不太高的情況下，沖擊壓縮產(chǎn)生的熵增不大，等熵過程與沖擊壓縮過程差別不大。因此，對于本文所研究的壓力范圍內(nèi)，可采用Murnaghan 狀態(tài)方程，并采取與Forbes［12］類似的推導(dǎo)方法，得到以下比焓與熵的表達(dá)式:

式中:下標(biāo)0 為初始狀態(tài);V0i為零壓下的比體積;具體溫度及相變模型參數(shù)取值見表1.

表1 溫度及相變模型參數(shù)Tab.1 Values of temperature and phase transition model parameters

其次，局部閾值條件，即當(dāng)自由能差額ΔG=Gα-Gε超過某一閾值時相變立即被觸發(fā)。計算中相變和逆相變的激活能分別為ΔGf和ΔGb.

再次，相變動力學(xué)方程，用于描述相變份額λE的變化率，本文分別采用1 階與2 階動力學(xué)方程，具體形式如下:

式中:λE為單元中ε 相的質(zhì)量份額;1 階相變動力學(xué)方程中λE∈［0，1］，初始時刻取λE=0;2 階相變動力學(xué)方程中λE∈［0.001，0.999］，初始取λE=0.001;τE為單元相變松弛時間，由單元直徑與馬氏體相變速率之比進(jìn)行估算，文中相變速率近似取常數(shù)1 km/s.

α-ε 相變伴隨著約5%的體積變化，數(shù)值模擬中通過依據(jù)單元中ε 相的質(zhì)量份額改變單元半徑來模擬這一效應(yīng):

式中:r0為單元初始半徑;V0α、V0ε分別為零壓狀況下a 和ε 相的比體積。

本文在計算中僅考慮了靜水壓力而未考慮偏應(yīng)力對于相變過程的影響。

2 α 鐵的數(shù)值模擬結(jié)果與分析

計算模型飛片與靶板均為6 400 μm×27 μm 的α 鐵，單元直徑為9 μm，計算單元總數(shù)約5 000 個。如圖2 所示，計算模型的上、下邊界采用周期性邊界條件，左、右邊界采用自由邊界條件。初始溫度300 K，飛片不同的速度從左端撞擊靶板。

圖3 給出了準(zhǔn)靜態(tài)加載條件下鐵的相圖。從圖3 可以看出，本文計算結(jié)果與文獻(xiàn)［16 -19］實(shí)驗(yàn)結(jié)果符合較好，證明本文采用的無擴(kuò)散兩相模型能夠正確反應(yīng)α 鐵的相變特性。

圖2 計算模型示意圖Fig.2 Initial model

圖3 準(zhǔn)靜態(tài)加載條件下鐵的相圖Fig.3 Phase diagram of iron under quasi-static loading

數(shù)值模擬結(jié)果顯示室溫300 K 下平衡態(tài)相邊界通過11.15 GPa 壓力點(diǎn)。Barker 等［20］通過實(shí)驗(yàn)發(fā)現(xiàn)鐵的沖擊相變壓力閾值在(12.9 GPa，13.7 GPa)范圍內(nèi)，逆相變壓力閾值在(9.4 GPa，10.2 GPa)范圍內(nèi)。因此，文中取相變的激活能為13.40 J/g，對應(yīng)室溫條件下相變的臨界壓力13.37 GPa;逆相變的激活能取為-8.24 J/g，對應(yīng)室溫條件下逆相變的臨界壓力為9.80 GPa.

圖4(a)給出了Fe 的p-V/V0Hugonoit 曲線;圖4(b)給出了Fe 的沖擊波波速D 與波后粒子速度up的關(guān)系圖。由圖4 可以看出，本文計算結(jié)果與文獻(xiàn)［1，20 -21］實(shí)驗(yàn)結(jié)果符合較好，證明計算所用單元間作用力模型及相變模型能夠正確反應(yīng)鐵的沖擊壓縮及相變特性。

此外，圖4(a)的模擬結(jié)果顯示，沖擊相變壓力閾值約12.95 GPa，比室溫下的結(jié)果13.37 GPa 低了約0.42 GPa，這說明沖擊導(dǎo)致靶板內(nèi)溫度升高，從而致使相變的壓力閾值降低。圖4(b)結(jié)果表明:當(dāng)沖擊壓力小于相變壓力閾值時，樣品內(nèi)波形為單波結(jié)構(gòu);當(dāng)沖擊壓力在相變壓力閾值至約36 GPa 之間，波形為雙波結(jié)構(gòu)(沖擊波與相變波);當(dāng)沖擊壓力繼續(xù)增大時，相變波與沖擊波匯合成穩(wěn)定的沖擊波。當(dāng)沖擊壓力低于相變壓力閾值時，沖擊波速度隨波后粒子速度近似呈線性增長;當(dāng)沖擊壓力高于相變壓力閾值時沖擊波近似以定常速度傳播，而相變波波速隨波后粒子速度增加較快。

圖4 鐵的沖擊Hugoniot 關(guān)系Fig.4 Hugoniot data of Iron

3 多晶鐵的數(shù)值模擬結(jié)果與分析

沖擊壓縮過程中多晶金屬，在晶粒尺度上是一個各向異性、非平衡過程，此時的沖擊波結(jié)構(gòu)不能為傳統(tǒng)的連續(xù)介質(zhì)力學(xué)描述［22］。Meyers 和Carvalho應(yīng)用晶粒取向的概然論模型研究了多晶鎳的沖擊響應(yīng)，結(jié)果表明傳入多晶中的沖擊波前沿變得不規(guī)則，其波面不規(guī)則度與晶粒尺寸同量級［23］。本文構(gòu)建了兩種不同晶粒尺度的多晶鐵模型(如圖5)，采用離散元方法研究了晶粒尺度對于沖擊波波面不規(guī)則程度的影響。

圖5 兩種不同晶粒尺度多晶鐵模型局部放大示意圖Fig.5 Enlarged view of initial model of polycrystalline iron

多晶鐵模型中的晶粒取向呈隨機(jī)分布，以此表征晶粒間取向的分布差異。各晶粒的取向?yàn)?0° ～60°)之間的隨機(jī)數(shù)，晶粒內(nèi)部單元為規(guī)則的密排六邊形，晶界處的間隙由小尺度單元填滿，圖5 中不同灰度代表不同晶粒取向，兩種模型的具體參數(shù)參見表2 所示。為避免邊界效應(yīng)的影響，計算模型的上、下邊界為周期型邊界條件，左側(cè)為固壁邊界，右側(cè)為自由邊界，初始時刻模型以一定的速度撞擊固壁。

表2 兩種晶粒尺度多晶鐵模型參數(shù)Tab.2 Parameters of polycrystalline iron models

圖6 為兩種不同晶粒大小的多晶鐵以700 m/s的速度撞擊固璧后，某時刻模型中的沖擊波前沿位置圖。此處，定義波后粒子速度等于0.5 倍撞擊速度的位置為沖擊波前沿位置，圖中深色部分為沖擊壓縮區(qū)，淺色部分為未受沖擊區(qū)域。數(shù)值模擬結(jié)果表明，沖擊波前沿分布并不均勻，并且大晶粒模型內(nèi)沖擊波前沿的不規(guī)則程度更高。

圖7 為沖擊波前沿位置中線定義示意圖，圖中的黑色粗實(shí)線代表沖擊波前沿位置，陰影部分為波后的沖擊壓縮區(qū)。如圖7 所示，沿y 向在沖擊波輪廓線的波峰(y 向最大值)與波谷(y 向最小值)之間選取一條直線(平行于x 軸)，使得直線上方的陰影面積之和等于直線下方的空白面積之和，并定義該直線位置為沖擊波前沿位置中線。

圖6 兩種不同晶粒尺度多晶鐵中的沖擊波前沿位置Fig.6 The shock wave front in polycrystalline iron

圖7 沖擊波前沿位置中線的定義Fig.7 The definition of shock wave front midline

統(tǒng)計不同時刻兩種不同晶粒大小的多晶鐵模型中，沖擊波前沿位置距離中線位置的標(biāo)準(zhǔn)偏差。圖8 給出了沖擊波前沿位置的標(biāo)準(zhǔn)偏差隨沖擊波傳播距離的變化關(guān)系。從統(tǒng)計結(jié)果看，沖擊波的波面不規(guī)則程度隨傳播距離的增加而增加;晶粒大小對于沖擊波前沿的不規(guī)則程度有一定影響，大晶粒模型內(nèi)沖擊波前沿的不規(guī)則程度更高。

圖8 沖擊波前沿位置的標(biāo)準(zhǔn)偏差隨沖擊波傳播距離的變化Fig.8 Variation of standard deviation in the shock wave front position

圖9為模型1 以600 m/s 速度撞擊固壁后，沖擊波前沿附近的質(zhì)量份額與壓力圖，圖中Ve 代表相變質(zhì)量份額，p 為壓力。從圖可以看到:由于晶粒間取向分布差異及晶粒邊界效應(yīng)的影響，壓力及質(zhì)量份額場的分布很不均勻，應(yīng)力遠(yuǎn)高于或低于波后平均應(yīng)力的區(qū)域集中于晶界附近;從沖擊波前沿位置可見，由于應(yīng)力集中效應(yīng)的影響，相變首先發(fā)生在晶界處，然后穿透到晶粒內(nèi)部;在沖擊波前沿處的一些晶粒中，可觀察到指狀傳播的相邊界。

圖9 沖擊波前沿附近的質(zhì)量份額與壓力Fig.9 Transformed mass fraction and pressure fields

室溫下靜壓實(shí)驗(yàn)的測量結(jié)果表明，即使在準(zhǔn)靜態(tài)條件下，鐵的α→ε 相變也并非沿平衡面進(jìn)行［24］。Boettger 等的理論計算結(jié)果表明，沖擊相變過程中得到的p-λ 曲線與靜壓實(shí)驗(yàn)的測量結(jié)果非常相近［25］。本文采用寬度27 μm 的采樣窗，沿著波的傳播方向統(tǒng)計相變質(zhì)量份額與局部平均壓力。

圖10 為采用以上統(tǒng)計方法得到的局部平均壓力與相變質(zhì)量份額曲線。圖中的離散點(diǎn)為靜壓實(shí)驗(yàn)的測量結(jié)果［24］，兩條曲線分別為采用1 階與2 階動力學(xué)方程得到的統(tǒng)計結(jié)果。從圖10 中可見，沖擊相變的局部壓力-相變質(zhì)量份額曲線均逐漸逼近準(zhǔn)靜態(tài)壓縮實(shí)驗(yàn)曲線上的點(diǎn)。本文的統(tǒng)計結(jié)果與Boettger 等［25］的結(jié)果可以互相佐證。采用兩種動力學(xué)方程得到結(jié)果的重要不同點(diǎn)在于相變臨界壓力，1 階模型的相變臨界壓力約為8 GPa，2 階模型的相變臨界壓力為10 GPa;與1 階動力學(xué)方程的統(tǒng)計結(jié)果相比，2 階動力學(xué)方程得到的局部平均壓力-相變質(zhì)量分額曲線在低壓段(小于13 GPa)，與實(shí)驗(yàn)結(jié)果符合更好。

4 結(jié)論

本文利用離散元方法，結(jié)合基于無擴(kuò)散相變的兩相模型，模擬了α 鐵的沖擊相變過程，數(shù)值模擬結(jié)果表明:

1)計算得到鐵的相邊界、沖擊Hugoniot 關(guān)系與實(shí)驗(yàn)結(jié)果符合較好，證實(shí)采用離散元方法結(jié)合兩相模型模擬α 鐵沖擊相變過程是可行的。

圖10 局部平均壓力與相變質(zhì)量份額Fig.10 Local average pressure and transformed mass fraction

2)對多晶鐵的沖擊響應(yīng)過程進(jìn)行了模擬，初步的數(shù)值模擬結(jié)果顯示沖擊波的波面不規(guī)則程度隨傳播距離的增加而增加，大晶粒模型內(nèi)沖擊波前沿的不規(guī)則程度更高。

3)對多晶鐵沖擊相變過程進(jìn)行了統(tǒng)計，結(jié)果表明沖擊相變的局部壓力-相變質(zhì)量份額曲線均逐漸逼近準(zhǔn)靜態(tài)壓縮實(shí)驗(yàn)曲線上的點(diǎn);與1 階動力學(xué)方程的相比，2 階動力學(xué)方程得到的局部平均壓力-相變質(zhì)量分額曲線在低壓段與實(shí)驗(yàn)結(jié)果符合更好。

References)

［1］Bancroft D，Peterson E L，Minshall S. Polymorphism of iron at high pressure［J］. J Appl Phys，1956，27(3):291 -298.

［2］Cundall P A. A computer model for simulating progressive large scale movement in block rock system［C］∥Symposium ISRM.Nancy，F(xiàn)rance:International Society of Rock Mechnics，1971:129 -136.

［3］Sawamoto Y，Tsubota H，Kasai Y，et al. Analytical studies on local damage to reinforced concrete structures under impact loading by discrete element method ［J］. Nucl Eng Des，1998，179:157 -177.

［4］劉凱欣，高凌天，鄭文剛. 混凝土動態(tài)破壞過程的數(shù)值模擬［J］. 工程力學(xué)，2000(增刊):470 -474.LIU Kai-xin，GAO Ling-tian，ZHENG Wen-gang. Numerical simulation for the concrete failure process under shock loading［J］.Engineering Mechanics，2000(Suppl)470 -474. (in Chinese)

［5］劉凱欣，高凌天. 離散元法在求解三維沖擊動力學(xué)問題中的應(yīng)用［J］. 固體力學(xué)學(xué)報，2004，25(2):181 -185.LIU Kai-xin，GAO Ling-tian. The application of discrete element method in solving three-dimensional impact dynamics problems［J］. Acta Mechanica Solida Sinica，2004，25(2):181 -185.(in Chinese)

［6］Yano K，Horie Y. Discrete-element modeling of shock compression of polycrystalline copper ［J］. Phys Rev B，1999，59:13672 -13680.

［7］Case S，Horie Y. Mesomechanics of the α-δ transition in iron［J］. J Mech Phys Solids，2007，55:589 -614.

［8］王文強(qiáng). 離散元方法及其在材料和結(jié)構(gòu)力學(xué)響應(yīng)分析中的應(yīng)用［D］. 合肥. 中國科學(xué)技術(shù)大學(xué)，2000.WANG Wen-qiang. Discrete element method and its use in analysis of response of materials and structures［D］. Hefei:University of Science and Technology of China，2000. (in Chinese)

［9］于繼東，王文強(qiáng)，劉倉理，等. 炸藥沖擊響應(yīng)的二維細(xì)觀離散元模擬［J］. 爆炸與沖擊，2008，28(6):488 -493.YU Ji-dong，WANG Wen-qiang，LIU Cang-li，et al. Two-dimensional mesoscale discrete element simulation of shock response of explosives［J］. Explosion and Shock Waves，2008，28(6):488 -493. (in Chinese)

［10］Tang Z P，Horie Y，Psakhie S G. Discrete meso-element modeling of shock processes in powders［C］∥High Pressure Shock Compression of Solids Ⅳ，Response of Highly Porous Solid to Shock Loading. N Y:Springer，1997:143 -176.

［11］Yano K，Horie Y. Mesomechanics of the α-ε transition in iron［J］. International Journal of Plasticity，2002，18:1427 -1446.

［12］Forbes J W. Experimental investigation of the kinetics of the shock-induced alpha to epsilon phase transformation in armco iron，WSU-SDL 76-01［R］. Pullman，Washington:Washington State University，1976:112 -120.

［13］Anderson O L. Equations of state of solids for geophysics and ceramic science［M］. New York:Oxford University，1995:195.

［14］Lide D R，Kehiaian H V. CRC handbook of thermophysical and thermodynamical data［M］. US:CRC，1994:136.

［15］Jephcoat A P，Mao H K，Bell P M. The static compression of iron to 78 GPa with rare gas solids as pressure-transmitting media［J］. J Geophys Res，1986，91:4677 -4684.

［16］Giles P M，Longenbach M H，Marder A R. High-pressure α-ε martensitic transformation in iron ［J］. J Appl Phys，1971，42(11):4290 -4295.

［17］Bundy F P. Pressure-temperature phase diagram of iron to 200 kbar，900 ℃［J］. J Appl Phys，1965，36(2):616 -620.

［18］Kennedy G C，Newton R C. Solids under pressure［M］. New York，McGrawHill，1963:163.

［19］Clougherty E V，Kaufman L. High pressure measurement［M］.Washington，Giadini A A，Lloyd E C，1963:152.

［20］Barker L M，Hollenbach R E，Shock wave study of the αTMε phase transition in iron［J］. J Appl Phys，1974，45:4872 -4887.

［21］Brown J M，F(xiàn)ritz J N，Hixson R S. Hugoniot data for iron［J］.J Appl Phys，2000，88:5496 -5498.

［22］Wallace D C. Irreversible thermodynamics of flow in solids［J］.Phys Rev B，1980，22(4):1477 -1486.

［23］Meyers M A ，Carvalho M S. Shock-front irregularities in polycrystalline metals［J］. Mater Sci Eng，1976，24:131 -135.

［24］Taylor R D，Pasternak M P，Jeanloz R. Hysteresis in the high pressure transformation of bcc-to hcp-iron ［J］. J Appl Phys，1991，69(8):6126 -6128.

［25］Boettger J C，Wallace D C. Metastability and dynamics of the shock-induced phase transition in iron［J］. Phys Rev B，1997，55:2840 -2849.