高純銅界面Rayleigh-Taylor不穩定性擾動增長的數值模擬

李碧勇，彭建祥，谷巖，尹曉春，賀紅亮

(1.南京理工大學 理學院，江蘇 南京 210094;2.中國工程物理研究院 流體物理研究所，四川 綿陽 621999)

0 引言

當低密度流體(介質)加速高密度流體(介質)時，界面會出現Rayleigh-Taylor(RT)不穩定性現象[1-2]。RT不穩定性是超新星爆發和星系演變的重要過程，也是慣性約束聚變的關鍵因素[3]。殼體界面上一個微小的初始擾動在沖擊作用下將會隨著時間的推移而增長并導致后期的湍流混合，這種不穩定性會直接影響到內爆壓縮效率[4]。目前，對流體或者氣體界面RT不穩定性問題的理論研究[5-6]、實驗研究[7-8]和數值模擬[9-10]等方面都取得了許多重要成果，內容涉及流體黏性、表面張力以及不穩定性發展的線性和非線性等問題，但是對于固體界面RT不穩定問題的研究相對較少[11-12]。

在內爆動力學中，金屬界面的不穩定性發展機制與許多復雜物理機制(如微噴和熔化破碎等)共同存在并且相互耦合，還涉及材料本構、動態斷裂和沖擊相變(熔化)等多物理過程。因此，金屬界面的不穩定性狀態相對于流體更加復雜，在一定程度上導致了一些公開發表的文獻資料中出現明顯相反的結論。例如，針對半無限大彈性體，Miles[13]根據能量平衡方程給出了RT不穩定性擾動增長的截止波長判據，認為在該截止波長以下的界面擾動不增長。Drucker[14]采用針對半空間的二維剛性碰撞Prandtl模式解的方法給出了與材料強度有關的臨界幅值判據，認為界面擾動的增長與初始幅值有關，小于臨界幅值的界面擾動將不增長。Swegle等[15]通過拉式有限差分程序對金屬鋁板進行了數值模擬研究，認為界面擾動局部的應力梯度才是控制RT不穩定性擾動增長的主要因素。上述文獻中都近似將材料假設為理想彈性體，基于小擾動近似且考慮了其在加載條件為常壓下的不穩定性增長。但實際上在爆轟作用下，樣品表面的加載壓力始終是處于變化的，其加速度并非常數。此外，材料界面擾動處于大變形的運動狀態，不但有彈性變形，而且更多的應該是非線性的塑性變形。因此，上述擾動發展的相關判據均存在一定的局限性。由此可見，材料性質和界面擾動的初始特征對固體材料RT不穩定性的影響還沒有形成共識，需要進一步深入研究。

本文以工業中常用的高純銅(OFHC)材料[16]為例，采用二維有限差分程序開展數值模擬，研究平面幾何條件下金屬界面的RT不穩定性擾動增長行為，從界面的初始擾動特征、樣品厚度和材料強度特性等方面分析影響界面擾動發展的相關因素，促進對相關物理問題的認識和理解。

1 數值模擬

基于AUTODYN有限差分程序，建立爆轟加載條件下金屬界面RT不穩定性問題的二維數值模型，如圖1所示。圖1中，a0為擾動的初始幅值(波峰到波谷的距離)，λ為擾動的初始波長，t0=0 μs為初始時刻，a(x，t)為當前時刻t界面位置x處所對應的擾動函數，ρ為材料密度。數值模型采用平面應變單元，樣品后自由面為完整的平面(無擾動)，在樣品前界面(加載面)上預制特定波長和幅值的余弦型初始小擾動，在初始時刻界面位置x處對應的界面擾動函數形式可以表示為a(x，0)=a0cos (2π/λx)/2.在當前時刻t界面擾動幅值可以進一步寫為a(t)=|a(λ/2，t)-a(0，t)|.

圖1 RT不穩定性數值模型Fig.1 Numerical model of Rayleigh-Taylor instability

為了簡化計算，本文只選取初始波長的1/2進行計算(0≤x≤λ/2)，其中上下邊界分別對應初始擾動的波峰和波谷。對于RT不穩定性問題，界面擾動的波長在演化過程中基本保持不變[17-18]，因此根據對稱性可以在模型的上下邊界設置周期性邊界條件。

Barnes等[12]開展的實驗中，爆轟產物經過一定的間隙后作用于待研究的金屬樣品，間隙的作用使得在預制擾動的界面形成了無沖擊加載(近似于準等熵)，避免了樣品因為沖擊造成劇烈的溫升、層裂或其他形式的擾動發生。一方面，在對炸藥和金屬樣品的建模過程中，由于炸藥的厚度遠大于樣品的厚度，炸藥網格數量占總網格比例較大，對炸藥爆轟過程的計算將耗費大量時間，導致模擬計算的效率較低；另一方面，相對于金屬，炸藥爆轟產物的密度很小，但其聲速卻很大。因此，在數值模擬中可以將炸藥模型做進一步簡化，通過在加載面施加壓力邊界條件的方式替代炸藥的加載作用[18]。本文計算中所采用的壓力邊界條件數據取自于文獻[12]中對鋁材料開展的RT不穩定性實驗，如圖2所示。

圖2 樣品加載面處的壓力曲線[12]Fig.2 Applied pressure on loading surface[12]

對于大多數金屬材料而言，體內的沖擊波速度D與粒子速度u可以采用簡單的線性關系描述為D=C0+su，其中C0和s為實驗擬合得到的樣品沖擊雨貢紐參數。材料狀態方程選為常用的Mie-Gruneisen形式：

p-pH=Γρ(E-EH)，

(1)

式中：p、E和pH、EH分別表示當前狀態和Hugoniot狀態下的壓力和內能；Γ為Gruneisen系數。考慮到固體的RT不穩定性具有高壓和高應變率的特點，典型的壓力范圍從10～100 GPa、應變率更是高達105～109s-1.因此，本文采用目前對高壓區描述最好的Steinberg-Cochran-Guinan(SCG)模型來描述金屬材料在高應變率(≥105s-1)情況下的RT不穩定性彈塑性變形特性[19-20]，其具體形式表示為

(2)

式中：G為剪切模量；Y為屈服強度；G0和Y0為常溫常壓下的剪切模量和屈服強度；η為材料壓縮比；ε為等效塑性應變；εi為初始塑性應變(通常情況下為0)；β、n為應變硬化參數；T為材料內溫度；G′p、G′T和Y′p分別為在參考狀態下G和Y對壓力p或溫度T的1階偏導數。表1分別給出了無氧銅(OFHC)和工業純鋁(1100-Al)的SCG模型參數。

表1 金屬材料的SCG模型參數Tab.1 Material parameters of SCG model

注：Ymax為應變硬化允許的屈服強度最大值，Tm為與壓力相關的熔化溫度。

為了排除網格尺寸的影響并保證計算結果的收斂性，分別對不同網格尺寸下界面擾動幅值的增長行為進行數值計算，并將數值模擬的結果與文獻[12]中的實驗結果進行了對比。根據Barnes等[12]給出的實驗設計，計算模型中的樣品材料選為1100-Al，樣品厚度為2.54 mm，界面擾動的初始波長為5.08 mm，初始幅值為0.203 mm，模型參數見表1.不同網格尺寸下界面擾動幅值增長的計算和實驗結果如圖3所示。

圖3 樣品材料為1100-Al時不同網格尺寸下擾動幅值隨時間的演化關系Fig.3 Perturbation growth amplitudes of 1100-Al material under different mesh sizes

從圖3中可以看出，對于本文所建立的金屬界面RT不穩定性二維數值模型，隨著模型網格尺寸的減小，界面擾動幅值的計算結果逐漸趨于一致，體現了較好的收斂性，但是相對來講，網格越細、計算時間就越長，效率也越低。同時將本文計算結果與文獻[12]中的實驗結果進行對比，結果表明本文計算結果與實驗數據基本吻合。此外還將本文計算結果與文獻[12]中MAGEE二維拉格朗日程序所計算的結果進行了比較，發現擾動幅值隨時間增長的數據也相一致，進一步表明本文數值模擬方法和結果是可靠的。因此，在后文的數值計算中，OFHC數值模型均采用127×127規格的網格劃分，最小的網格尺寸約為20 μm.

2 計算結果與分析

2.1 界面初始特征的影響

在內爆動力學所關注問題的研究中，金屬界面RT不穩定性與界面的初始擾動特征密切相關，但是目前沒有形成統一的認識，從彈性理論推導的各種穩定性判據都具有一定局限性，對加載壓力隨時間變化的RT不穩定性增長行為也沒有給出相應的結論。本文的數值模擬將通過改變樣品表面擾動的初始幅值和初始波長等，研究爆轟加載下界面的初始擾動特征對OFHC的RT不穩定性擾動增長的影響。為了便于比較，利用界面擾動幅值的無量綱增長因子描述擾動的增長，其定義為當前時刻材料界面擾動幅值a(t)與初始幅值a0之比。

2.1.1 初始幅值的影響

對于圖1所示的二維數值模型，OFHC樣品初始厚度H=2.54 mm，假設界面的初始擾動按照余弦函數形式分布，擾動的初始波長λ=5.08 mm，擾動初始幅值(波峰到波谷的距離)a0分別取值為0.10 mm、0.20 mm、0.25 mm、0.30 mm、0.40 mm、0.60 mm共6種情況，OFHC的模型參數見表1.通過模擬得到不同初始幅值下樣品內部的壓力分布云圖和界面擾動增長因子的計算結果，分別如表2和圖4所示。

表2 不同時刻界面擾動增長的壓力云圖Tab.2 Pressure contours of perturbation growth at different times

圖4 初始幅值對擾動增長的影響Fig.4 Effect of initial amplitude on perturbation growth factor

從表2中可以看出，擾動波峰和波谷位置處對應的界面壓力基本相同，然而在波峰位置處，樣品的等效厚度比波谷位置處的要大，因此，在相同初始波長和加載壓力作用下，界面擾動的初始幅值越大，波谷位置處擾動增長的加速度就越大，波峰和波谷位置由于加速度的差異所導致的界面運動速度差也就越大，初始幅值越大的界面擾動就更容易增長。

從圖4中可以看出，在界面擾動發展的初期，不同的初始幅值所對應的擾動增長因子的變化基本一致，但在大約1.5 μs之后，不同初始幅值對應的界面擾動增長趨勢呈現出明顯差異。從圖2中加載面的壓力曲線可以看到：1.5 μs正好對應于壓力達到峰值的時刻，這種變化與Swegle等[15]提出的關于RT不穩定性擾動增長的應力梯度理論一致；當擾動的初始幅值a0=0.10 mm時，在t≈5 μs時刻，擾動幅值的增長系數從1.0逐漸增大到2.3后就不再繼續增長，并呈現出上下振動的穩定發展形式；當a0=0.25 mm時，擾動幅值的增長因子隨時間基本呈現出線性增長，在t=8 μs時刻擾動增長因子等于6.4，擾動幅值大約為初始波長的0.3倍，按照Birkhoff[17]對RT不穩定性發展階段的劃分，此時界面擾動處于指數發展階段或小擾動階段，界面上的擾動按照指數形式發展，擾動的基本形狀不變；當a0>0.25 mm時，如a0分別取值為0.30 mm、0.40 mm、0.60 mm情況下，初始幅值越大，擾動增長速率也越大，且界面擾動的發展在后期呈現出非線性增長的趨勢，如果不考慮加載面卸載的情況，則界面擾動將一直增長，直到材料滿足相應的斷裂準則而發生斷裂為止。因此，可以認為a0=0.25 mm為界面擾動增長穩定與不穩定的臨界幅值，當初始幅值小于0.25 mm時，界面擾動的增長速率不斷減小，其狀態可以認為是穩定的；當初始幅值大于0.25 mm時，界面擾動的增長速率卻不斷增大，其狀態是不穩定的。根據Drucker的臨界幅值理論[14](ac=2Y0/(ρg)，ac為臨界幅值，g為重力加速度)，計算出擾動增長的臨界幅值為0.06 mm，明顯低于本文數值模擬的計算值，這也表明了臨界幅值理論的局限性。然而Drucker的理論雖然不適用于本文爆轟加載條件下的擾動增長分析，但是仍然可以借鑒并從定性角度判斷，界面擾動增長的臨界幅值與材料的屈服強度Y0、樣品的厚度H和界面加載壓力p的大小密切相關。材料的屈服強度和厚度越大或者加載壓力越小，材料就越不容易發生塑性變形，相應地擾動不穩定增長的臨界幅值ac也就越大，小于ac的界面擾動增長就處于穩定狀態。

2.1.2 初始波長的影響

為了研究初始波長對界面擾動增長的影響，分別計算初始幅值a0為0.20 mm和0.40 mm兩種情況下，初始波長λ從1.27 mm逐漸增大到10.16 mm不同情況下的擾動增長因子，結果如圖5所示。由圖5可見：在界面擾動增長早期，雖然初始波長不同，但擾動增長的趨勢卻基本一致，且增長幅度也很小，直到1.5 μs左右樣品加載面處的輸入壓力達到最大值，擾動幅值才開始迅速地增長；當初始幅值為0.20 mm時，各種初始波長情況下界面擾動幅值的增長因子均增加到一定值后都不再隨時間繼續增長，而是處于一種上下振動的穩定狀態。

圖5 初始波長對擾動增長的影響Fig.5 Effect of initial wavelength on perturbation growth factor

圖6所示為初始波長對擾動界面加速度的影響。由圖6(a)可知：當初始幅值小于臨界幅值時，初始波長較小的界面擾動的波峰和波谷位置處加速度呈現交替變化規律，二者差值隨著時間演化逐漸變小，因此后期波峰和波谷位置的擾動增長速度基本相同，從而界面擾動幅值趨于穩定，而初始波長較大的界面擾動，在擾動增長早期，波谷位置處的加速度始終要大于波峰位置處的加速度，直到后期二者才基本一致，從而界面擾動具有較大的增長趨勢并在后期趨于穩定，初始波長越大，最終的擾動增長因子也越大；當初始幅值為0.40 mm時，不同初始波長下的界面擾動幅值都呈現出不穩定發展趨勢，隨著初始波長的增加，在相同的時刻界面擾動幅值增長因子卻反而減小。從圖6(b)中可以看出，當初始幅值大于臨界幅值時：對于初始波長較小的擾動，波谷位置的加速度始終要大于波峰位置處的加速度，因此界面擾動增長越來越大；對于初始波長較大的擾動，二者呈現交替變化的規律，導致波峰和波谷位置處的速度差變小，從而初始波長越大，在相同時刻界面擾動增長也就越小。

圖6 初始波長對擾動界面加速度的影響Fig.6 Effect of initial wavelength on perturbation growth acceleration

2.1.3 樣品厚度的影響

RT不穩定性的本質是輕介質加速重介質，因此介質中的加速度也是一個非常重要的參數。本文假設樣品材料為不可壓縮的，在外部載荷不變的情況下，樣品的質量(厚度)為決定其加速度的唯一參數。通過數值模擬得到不同樣品厚度H對OFHC界面擾動增長行為的影響如圖7所示。從圖7中可以清楚地看到：樣品表面RT不穩定性擾動的增長存在著一個臨界厚度Hc，此臨界厚度Hc又與界面擾動初始幅值密切相關，當初始幅值a0=0.20 mm時，1.91 mm

圖7 樣品厚度對擾動增長的影響Fig.7 The effect of specimen thickness on the perturbation growth factor

2.2 材料力學性質的影響

相比流體或氣體，固體材料中RT不穩定性的增長狀態相對復雜，其中一個重要的原因就是固體材料具有一定的強度[4，11]。材料強度在本質上反映的是材料在沖擊載荷下抵抗剪切變形的能力，然而在經典高壓物態方程的研究中，通常忽略剪切應力的影響，將材料受力狀態簡化為流體靜水壓力狀態，并將實驗測量中得到的正應力直接作為壓力代入狀態方程和熱力學計算。然而只要材料沒有發生熔化，其從沖擊壓縮狀態卸載時就開始表現出典型的彈- 塑性卸載特征[21]。Morris等[22]通過實驗研究發現，在沖擊加載壓力達到80 GPa左右時，2024-鋁和高導無氧銅仍然具有幾個兆帕的屈服強度。因此，材料的力學特性(如強度和模量等)會對界面的擾動發展產生重要影響。從美國、俄羅斯各大實驗室公開發表的文獻資料可以看出，金屬材料強度對RT不穩定性的影響也是目前研究的熱點問題。Kalantar等[23]對鋁材料開展了準等熵壓縮實驗，研究了峰值壓力高達180 GPa、應變率在107s-1狀態下的RT不穩定性，發現在擾動增長后期(>20 ns)，界面擾動的增長達到飽和或受到抑制，推測此時材料的強度效應發揮了重要作用。Park等[24]研究了金屬釩材料在106～108s-1超高應變率和10～100 GPa加載壓力下的RT不穩定性，結果發現如果數值模型中不考慮材料強度，則在t=60 ns時刻，擾動幅值的增長因子將高達實驗結果的5倍左右。

本文進一步研究了表征材料彈塑性力學行為的兩個重要參數(材料的屈服強度Y0和剪切模量G0)對OFHC界面擾動增長的影響。數值計算中分別單獨調整SCG模型中的參數Y0和G0的值，與此同時其余參數均保持不變，當材料參數Y0=G0=0 GPa時，SCG模型退化為純流體模型，即不考慮材料的本構關系。經過一系列數值計算，得到了界面擾動增長因子隨著屈服強度Y0和剪切模量G0變化的模擬結果，分別如圖8和圖9所示。

圖8 屈服強度對擾動增長的影響Fig.8 Effect of material yield strength on perturbation growth

圖9 剪切模量對擾動增長的影響Fig.9 Effect of material shear modulus on the perturbation growth

從圖8和圖9中可以看出：在擾動增長的初期(t<1.5 μs)，界面擾動的增長和純流體模型的計算結果相似，幾乎與材料的強度無關；當t>1.5 μs時，擾動增長因子的變化就與材料的強度密切相關。從圖9中可見：屈服強度Y0越大，材料抵抗剪切變形的能力就越強，從而界面擾動幅值的增長也就越小，當Y0=1.0 GPa時，界面擾動幅值的增長幾乎為0.圖9中剪切模量對擾動增長的影響與圖8中的結果不一致。在擾動增長到1.5 μs時刻之后，純流體模型(G0=0 GPa)情況下的擾動增長與考慮材料強度(G0>0 GPa)的擾動增長情況存在著顯著的區別。總體而言，隨著G0的增大，界面擾動增長因子也相應地減小。在擾動增長中期(1.5 μs5 μs時，相互之間的差別隨著時間的推移呈現逐漸增大的現象，但是這種變化卻隨著G0的增大而減小，例如當G0從47.7 GPa逐漸增大到95.4 GPa時，界面擾動的增長系數逐漸趨于一致。數值模擬結果表明：相對于材料的屈服強度Y0，剪切模量G0對擾動幅值增長過程的影響可以忽略。

3 結論

本文通過數值模擬，研究了爆轟加載下OFHC界面的初始擾動特征、樣品厚度和材料強度特性對RT不穩定性擾動增長的影響。得出主要結論如下：

1)界面擾動的增長與初始幅值和初始波長都密切相關，并存在一個臨界幅值(其與材料的屈服強度、樣品厚度和加載壓力密切相關)。

2)在初始幅值小于臨界幅值的情況下，界面擾動的增長隨波長的增加而增大；當初始幅值大于臨界幅值時，界面擾動的增長情況則相反，初始波長越大，擾動幅值的增長就越小。

3)樣品厚度對界面擾動增長也起到了很大作用，樣品厚度大于臨界厚度時，界面擾動的增長處于穩定狀態。

4)在擾動增長初期，擾動幅值的增長和純流體模型的計算結果相似，幾乎與材料的強度無關；在擾動增長后期，擾動增長因子的變化就與材料的強度密切相關。屈服強度越大，擾動幅值的增長也就越小，當屈服強度增長到一定程度后，界面擾動幅值的增長幾乎為0，相對于材料的屈服強度，剪切模量對擾動幅值增長過程的影響可以忽略。