高純銅界面Rayleigh-Taylor不穩(wěn)定性擾動增長的數(shù)值模擬

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

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

0 引言

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

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

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

1 數(shù)值模擬

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

圖1 RT不穩(wěn)定性數(shù)值模型Fig.1 Numerical model of Rayleigh-Taylor instability

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

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

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

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

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

(1)

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

(2)

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

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

注：Ymax為應(yīng)變硬化允許的屈服強度最大值，Tm為與壓力相關(guān)的熔化溫度。

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

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

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

2 計算結(jié)果與分析

2.1 界面初始特征的影響

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

2.1.1 初始幅值的影響

對于圖1所示的二維數(shù)值模型，OFHC樣品初始厚度H=2.54 mm，假設(shè)界面的初始擾動按照余弦函數(shù)形式分布，擾動的初始波長λ=5.08 mm，擾動初始幅值(波峰到波谷的距離)a0分別取值為0.10 mm、0.20 mm、0.25 mm、0.30 mm、0.40 mm、0.60 mm共6種情況，OFHC的模型參數(shù)見表1.通過模擬得到不同初始幅值下樣品內(nèi)部的壓力分布云圖和界面擾動增長因子的計算結(jié)果，分別如表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中可以看出，擾動波峰和波谷位置處對應(yīng)的界面壓力基本相同，然而在波峰位置處，樣品的等效厚度比波谷位置處的要大，因此，在相同初始波長和加載壓力作用下，界面擾動的初始幅值越大，波谷位置處擾動增長的加速度就越大，波峰和波谷位置由于加速度的差異所導(dǎo)致的界面運動速度差也就越大，初始幅值越大的界面擾動就更容易增長。

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

2.1.2 初始波長的影響

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

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

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

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

2.1.3 樣品厚度的影響

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

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

2.2 材料力學(xué)性質(zhì)的影響

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

本文進(jìn)一步研究了表征材料彈塑性力學(xué)行為的兩個重要參數(shù)(材料的屈服強度Y0和剪切模量G0)對OFHC界面擾動增長的影響。數(shù)值計算中分別單獨調(diào)整SCG模型中的參數(shù)Y0和G0的值，與此同時其余參數(shù)均保持不變，當(dāng)材料參數(shù)Y0=G0=0 GPa時，SCG模型退化為純流體模型，即不考慮材料的本構(gòu)關(guān)系。經(jīng)過一系列數(shù)值計算，得到了界面擾動增長因子隨著屈服強度Y0和剪切模量G0變化的模擬結(jié)果，分別如圖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)，界面擾動的增長和純流體模型的計算結(jié)果相似，幾乎與材料的強度無關(guān)；當(dāng)t>1.5 μs時，擾動增長因子的變化就與材料的強度密切相關(guān)。從圖9中可見：屈服強度Y0越大，材料抵抗剪切變形的能力就越強，從而界面擾動幅值的增長也就越小，當(dāng)Y0=1.0 GPa時，界面擾動幅值的增長幾乎為0.圖9中剪切模量對擾動增長的影響與圖8中的結(jié)果不一致。在擾動增長到1.5 μs時刻之后，純流體模型(G0=0 GPa)情況下的擾動增長與考慮材料強度(G0>0 GPa)的擾動增長情況存在著顯著的區(qū)別?？傮w而言，隨著G0的增大，界面擾動增長因子也相應(yīng)地減小。在擾動增長中期(1.5 μs5 μs時，相互之間的差別隨著時間的推移呈現(xiàn)逐漸增大的現(xiàn)象，但是這種變化卻隨著G0的增大而減小，例如當(dāng)G0從47.7 GPa逐漸增大到95.4 GPa時，界面擾動的增長系數(shù)逐漸趨于一致。數(shù)值模擬結(jié)果表明：相對于材料的屈服強度Y0，剪切模量G0對擾動幅值增長過程的影響可以忽略。

3 結(jié)論

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

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

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

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

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