金屬錫Rayleigh-Taylor不穩定性對模型參數敏感性的數值分析

王 濤，汪 兵，林健宇，柏勁松，李 平，鐘 敏，陶 鋼

（1.中國工程物理研究院流體物理研究所，四川 綿陽 621999；2.南京理工大學能源與動力工程學院，江蘇 南京 210094）

當低密度流體在慣性力作用下加速高密度流體時，界面上的預制初始擾動會增長，這種現象稱為Rayleigh-Taylor（RT）不穩定性[1-2]。同樣，在金屬材料界面也會發生RT不穩定性，產生與流體界面不穩定性類似的尖釘和氣泡結構，但是與流體界面不穩定性的最大不同之處在于金屬材料具有強度效應，而研究表明金屬材料強度能夠抑制擾動增長[3-6]。并且，金屬的強度效應也會使金屬界面不穩定性問題比流體界面不穩定性更復雜，需要考慮的因素更多，比如金屬材料的本構方程、損傷斷裂、相變、熔化、微結構特征等。金屬界面不穩定性的發生常常預示著材料經歷高溫、高壓、高應變率等極端環境，因而常存在于慣性約束聚變[7-9]、超新星爆炸[10]、小行星撞擊[11-12]、地球內核運動和板塊構造[13-14]以及爆炸焊接等領域，因此具有重要的研究意義。

目前，金屬材料界面不穩定性的理論研究主要利用能量平衡[15-17]或力平衡[18-20]原理，推導擾動增長所滿足的色散關系，即線性分析；由于這些線性分析方法均采用理想塑性本構關系，因此預測擾動增長受到較大的限制。Bai等[21]基于能量平衡原理，采用Steinberg-Guinan（SG）和Johnson-Cook（JC）本構模型及變壓力的準等熵加載歷程，針對有限厚度金屬平板問題，推導了相應的擾動增長方程，對線性分析方法進行了拓展，使其能夠用于平面幾何中金屬材料不穩定性線性至非線性段增長的預測，并對爆轟驅動和等離子體驅動的金屬界面不穩定性實驗[6,20,22-23]進行了分析，結果和實驗吻合較好，但是這種線性分析方法仍然存在較大的局限性，如目前還無法考慮加載過程等。

金屬材料RT不穩定性的實驗研究始于20世紀70年代，首次由Barnes等[22,24]通過炸藥爆轟準等熵驅動鋁平板實現，并利用高能X射線照相觀察界面擾動增長。之后的平面爆轟驅動金屬材料RT不穩定性實驗[25-27]均采用相似的裝置。美國和俄羅斯開展了大量有關金屬材料界面不穩定性的研究工作，包括界面擾動增長規律及影響因素等。Olson等[28]通過爆轟驅動實驗，研究了銅材料的微結構和機加工對RT不穩定性增長的影響，發現單晶晶向和材料加工導致的應變硬化會影響擾動增長，而多晶材料的晶粒尺寸和晶界強化在所研究的加載條件下對擾動增長無影響。Wang等[27]通過爆轟驅動實驗和數值模擬，研究了純鋁的RT不穩定性增長規律和強度因素（基于SG本構模型的初始剪切模量和初始屈服強度）對擾動增長的影響。何長江等[29]通過數值模擬研究了初始擾動波長對平面爆轟驅動鋁RT不穩定性增長的影響。郝鵬程等[30]和劉軍等[31]通過對內爆加載金屬鋼殼的數值模擬，研究了材料強度（剪切模量和基于理想彈塑性本構的屈服強度以及是否考慮彈塑性）和初始波長對界面擾動增長的影響。平面爆轟準等熵加載的壓力一般在30 GPa左右，為了提高加載壓力，Frahan等[26]采用兩級炸藥爆轟驅動技術，即第一級炸藥驅動鐵飛片撞擊起爆第二級炸藥，再加載樣品，使加載壓力提高到50 GPa，由此研究了該加載壓力下鈹的RT不穩定性問題，并結合數值模擬對材料不同本構模型的適用性進行了檢驗。為了研究極端加載壓力和應變率條件下的金屬界面不穩定性及材料混合行為，人們又發展了電磁加載[32-33]和激光加載技術[5-6,23]，并開展了大量研究工作，如Park等[5-6,34]利用Omega激光器，結合數值模擬，研究了釩和鉭分別在100 GPa和350 GPa的準等熵加載壓力下的擾動增長行為。目前，美國國家點火裝置NIF已經能夠實現超過1 TPa的準等熵加載[35]。

由于材料強度對金屬界面不穩定性發展有抑制作用，所以近年來人們又將其作為一種手段來研究金屬材料的強度—擾動增長法[4,25,36]。Lorenz等[23]利用基于激光準等熵加載技術的RT擾動增長法研究了峰值壓力為20 GPa時6061-T6鋁的屈服強度，結果顯示該條件下6061-T6鋁的屈服強度比正常條件下高出3.6倍。Park等[5-6]利用激光準等熵加載實驗和數值模擬研究發現了釩在100 GPa加載條件下具有一種新的RT不穩定性致穩機制—聲子拖曳導致有效晶格黏性增大，進而導致釩的強度增大，對RT不穩定性的致穩作用增強；而且研究顯示釩的強度增大了3.5倍，而350 GPa加載壓力下鉭的強度增大了13倍[34]。Belof等[37]通過爆轟驅動鐵的RT不穩定性實驗和數值模擬，確認了鐵在經歷α→ε→α′相變過程中強度大大提高。

本研究利用自研的爆轟與沖擊動力學歐拉計算程序和SG本構模型，針對爆轟驅動金屬錫RT不穩定性增長對兩類模型初始參數的敏感性進行數值模擬分析，一類是樣品的初始參數，包括初始擾動振幅、初始擾動波長和樣品初始厚度，另一類是SG本構模型的初始參數。

1 計算方法

對于考慮爆轟和材料彈塑性行為的多物質、大變形及強沖擊波物理問題，發展了高保真度的爆轟與沖擊動力學歐拉有限體積計算程序，其守恒型控制方程組為

式中：下標i、j分別代表x、y、z3個方向，遵循張量運算法則，V為控制體體積，S為控制體表面積，ni為外法向單位矢量，ρ、uk(k=i,j)、p、E分別為密度、速度、壓強和單位質量總能量，sij為偏應力張量分量。

采用維數分裂技術，將式(1)描述的物理問題分解為多個一維問題。對于每個一維問題，采用PPM（Piecewise parabolic method）方法對單元內物理量的分布進行插值和重構，然后通過歐拉型兩步算法進行計算，即物理量的Lagrange推進求解，再將拉氏網格上的物理量映射到靜止的歐拉網格上。材料強度效應、炸藥爆轟過程和人工黏性在Lagrange步中實現。多物質界面采用體積分數方法進行捕捉。

數值模擬中，炸藥和金屬材料分別采用JWL和Mie-Grüneisen狀態方程（EOS）描述

式中：v=ρ0/ρ為相對比容，T為溫度，α、σ、R1、R2、ω為常數，為單位體積內能。

金屬材料采用SG本構模型描述。SG本構模型是在彈塑性本構方程中引入壓力、溫度和應變率項，并且壓力與應變率對流動應力的耦合效應具有可分離變量特性。因為SG本構模型中流動應力依賴于壓力，所以材料的本構方程與狀態方程之間存在某種耦合關系。這一耦合關系反映了高壓下金屬材料的壓力硬化特性。SG本構模型給出的動態屈服強度YSG和剪切模量G分別為

式中：Y0和G0分別為初始屈服強度和剪切模量，β和n分別為材料應變硬化系數和硬化指數，A為壓力硬化系數，η=ρ/ρ0為材料壓縮比，B為熱軟化系數。

2 計算模型和條件

簡化的二維計算模型見圖1，炸藥為JO-9159，樣品為金屬錫。炸藥高度為30 mm；樣品直徑為32 mm，厚度3 mm，樣品面向炸藥一側的界面上預置了初始振幅a0= 0.1 mm、初始波長λ0= 3 mm的正弦形式初始擾動。炸藥和樣品之間的真空間隙尺寸為10 mm。該真空間隙保證了爆轟產物經過真空間隙加載樣品的過程中壓力連續增大，形成輕介質對重介質的準等熵加載，進而誘發RT不穩定性；它的另一個作用是確保在高壓下樣品溫度依然低于熔化溫度[23]。

圖1 二維計算模型Fig.1 Two dimensional computational model

JO-9159炸藥的JWL狀態方程參數見表1，錫的Mie-Grüneisen狀態方程參數和SG本構模型標準參數見表2和表3。其中：pCJ和DCJ分別為炸藥的爆壓和爆速，Ymax為材料最大屈服強度。

表1 JO-9159炸藥JWL狀態方程參數Table 1 EOS parameters of JO-9159 explosive

表2 錫的Mie-Grüneisen狀態方程參數Table 2 Mie-Grüneisen EOS parameters of Sn

表3 錫的SG本構模型參數Table 3 SG constitutive model parameters of Sn

3 結果分析

3.1 程序確認

在開展金屬錫RT不穩定性數值模擬研究之前，先通過對Lindquist等[38]的爆轟驅動T-6061鋁的RT不穩定性實驗進行模擬，確認本計算程序的可靠性。該實驗采用HMX炸藥爆轟加載1.5 mm厚的鋁樣品，加載壓力峰值為30 GPa。初始擾動波長λ0= 2 mm，初始擾動振幅a0= 0.15 mm, 0.11 mm。圖2給出了本研究通過自研歐拉程序計算得到的擾動振幅a和文獻中采用ARES程序計算的結果，及其與實驗結果的比較（橫軸y為自由面運動位移）。可以看出，吻合度較好，表明自研的歐拉程序模擬這類爆轟驅動金屬材料界面不穩定性是準確的。

圖2 Lindquist等爆轟驅動鋁實驗的擾動振幅比較Fig.2 Comparison of perturbation amplitudes of Lindquist et al.'s experiments driven by explosion

3.2 樣品初始參數敏感性分析

首先，考察模擬結果對計算網格的收斂性問題。圖3和圖4分別給出了網格尺寸為10.0、5.0和2.5μm時加載壓力剖面和擾動振幅增長曲線。可以看出，2.5 μ m時模擬結果已經達到收斂，所以后續的所有計算均采用2.5 μ m的計算網格。從圖3還可以看出：加載壓力p是連續光滑增大的，近似形成對界面的準等熵加載；之后由于爆轟產物膨脹壓力持續減小及由自由面多次反射的稀疏波抵達擾動界面形成持續的卸載作用，加載壓力逐漸減小。

圖3 不同網格尺寸時的加載壓力剖面Fig.3 Loading pressure profiles for different grid size

圖4 不同網格尺寸時的擾動振幅增長曲線Fig.4 Perturbation amplitude growth for different grid size

接下來，通過數值模擬研究金屬錫樣品初始參數，包括初始擾動振幅、初始擾動波長、樣品初始厚度對爆轟驅動金屬錫RT不穩定性增長的影響。在研究初始振幅的影響時，保持初始波長λ0為3 mm，樣品初始厚度h為3 mm，初始振幅a0分別取0.3、0.1、0.05、0.025和0.012 5 mm。圖5顯示了不同初始振幅時的擾動振幅增長曲線。可以看出：隨著初始振幅的減小，擾動振幅增長也減弱；當初始振幅減小到一定值后，擾動振幅不再增長，即爆轟驅動金屬錫RT不穩定性增長存在一個截止初始振幅。

圖5 不同初始振幅時的擾動振幅增長曲線Fig.5 Perturbation amplitude growth for different initial amplitude

在研究初始波長的影響時，保持初始振幅為0.1 mm，樣品初始厚度為3 mm，而初始波長λ0分別取9、6、3、1和0.5 mm。圖6給出了不同初始波長時的擾動振幅增長曲線。從初始波長為9 mm開始，隨著初始波長的減小，擾動振幅逐漸增大；但是，當初始波長減小到一定程度后，擾動振幅增長又隨著初始波長的減小而減小。所以爆轟驅動金屬錫RT不穩定性增長存在一個最不穩定的波長，或者增長速度最快的模態。

在研究樣品初始厚度的影響時，保持初始振幅為0.1 mm，初始波長為3 mm，而樣品初始厚度h分別取3、5、7和9 mm。圖7顯示了樣品初始厚度不同時的擾動振幅增長曲線。隨著樣品初始厚度的增大，擾動增長減弱，表明樣品厚度增大可以抑制擾動增長；而且當樣品厚度增大到一定程度時，擾動增長停止，即爆轟驅動金屬錫RT不穩定性增長存在一個樣品截止厚度。

圖6 不同初始波長時的擾動振幅增長曲線Fig.6 Perturbation amplitude growth for different initial wavelength

圖7 不同樣品初始厚度時的擾動振幅增長曲線Fig.7 Perturbation amplitude growth for different initial thickness of sample

3.3 SG本構模型參數敏感性分析

研究金屬錫RT不穩定性增長對SG本構模型參數的敏感性時，計算模型的初始厚度取3 mm，初始振幅取0.1 mm，初始波長取3 mm。SG本構模型的可調參數包括應變硬化系數β、應變硬化指數n、壓力硬化系數A和熱軟化系數B，對于錫材料，它們的標準值分別為2 000、0.06、0.086 6 GPa-1、2 120 K-1（見表3）。在研究其中一個參數的影響時，其余參數均取標準值。材料應變硬化系數β依次取1 000、1 500、2 000、2 500、3 000；應變硬化指數n依次取 0.01、0.03、0.06、0.12、0.24；壓力硬化系數A依次取 0.021 7、0.043 3、0.086 6、0.173 2、0.346 4 GPa-1；熱軟化系數B依次取 0.53×10?3、1.06×10?3、2.12×10?3、4.24×10?3、8.48×10?3K-1。

圖8～圖11分別給出了改變錫SG本構模型的應變硬化系數β、應變硬化指數n、壓力硬化系數A和熱軟化系數B時的擾動振幅增長曲線。可以看出：在爆轟準等熵驅動下金屬材料發生復雜流動，此時材料應變硬化系數和應變硬化指數的變化對擾動振幅增長的影響很小，即這兩個SG本構模型參數的改變對材料屈服強度的影響較小；而壓力硬化系數和熱軟化系數的變化對擾動增長的影響顯著，隨著壓力硬化系數的增大，擾動增長受到的抑制作用越顯著，而隨著熱軟化系數的減小，擾動增長受到抑制效果的變化卻逐漸減小，即熱軟化系數小到一定程度時其變化將很難改變材料屈服強度。根據文獻分析可知，高壓高應變率加載條件下，金屬材料的強度一般都會增大，而標準參數的本構模型一般都低估了材料強度。因此基于SG本構模型，采用擾動增長法預估高壓高應變率下錫的材料強度，修正壓力硬化系數以獲得合理的強度是一條合理的途徑。

圖8 應變硬化系數不同時的擾動振幅增長曲線Fig.8 Perturbation amplitude growth for different strain hardening coefficient

圖9 應變硬化指數不同時的擾動振幅增長曲線Fig.9 Perturbation amplitude growth for different strain hardening exponent

圖10 壓力硬化系數不同時的擾動振幅增長曲線Fig.10 Perturbation amplitude growth for different pressure hardening coefficient

圖11 熱軟化系數不同時的擾動振幅增長曲線Fig.11 Perturbation amplitude growth for different thermal softening coefficient

4 結 論

開展了爆轟驅動金屬錫RT不穩定性增長對樣品初始參數和SG本構模型初始參數敏感性的數值模擬分析，得出以下結論。

（1）樣品初始參數（初始振幅、初始波長、樣品初始厚度）對爆轟驅動錫RT不穩定性增長有重要影響。RT不穩定性增長隨著初始振幅的減小而減小，且存在一個截止初始振幅；存在一個最不穩定的模態（波長），當初始波長大于該波長時，RT不穩定性增長隨著初始波長的減小而增大，反之，RT不穩定性增長隨著初始波長的減小而減小；樣品厚度的增大可以抑制RT不穩定性增長，而且存在一個樣品截止厚度。

（2）金屬錫的RT不穩定性增長對其SG本構模型應變硬化系數和應變硬化指數的變化不敏感，而對壓力硬化系數和熱軟化系數比較敏感，但是從采用擾動增長法預估材料強度的角度來說，修正壓力硬化系數以獲得合理的錫材料強度是合理的途徑。