郝鵬程,馮其京,胡曉棉,2

(1.北京應(yīng)用物理與計算數(shù)學(xué)研究所，北京 100094;2.北京應(yīng)用物理與計算數(shù)學(xué)研究所計算物理國家重點(diǎn)實(shí)驗(yàn)室，北京 100088;3.中國工程物理研究院研究生部，北京 100088)

內(nèi)爆加載金屬界面不穩(wěn)定性的數(shù)值分析*

郝鵬程1,2,3,馮其京1,胡曉棉1,2

(1.北京應(yīng)用物理與計算數(shù)學(xué)研究所，北京 100094;2.北京應(yīng)用物理與計算數(shù)學(xué)研究所計算物理國家重點(diǎn)實(shí)驗(yàn)室，北京 100088;3.中國工程物理研究院研究生部，北京 100088)

采用自行研制的多介質(zhì)彈塑性流體力學(xué)歐拉程序，對柱形內(nèi)爆加載金屬界面不穩(wěn)定性進(jìn)行了數(shù)值研究，數(shù)值模擬結(jié)果與文獻(xiàn)實(shí)驗(yàn)數(shù)據(jù)吻合較好。數(shù)值結(jié)果表明：材料強(qiáng)度對界面不穩(wěn)定性發(fā)展有不可忽略的抑制作用；材料屈服強(qiáng)度對較高模數(shù)不穩(wěn)定性增長的抑制較強(qiáng)，而剪切模量對不穩(wěn)定性發(fā)展的影響相似但敏感性相對較弱；金屬界面不穩(wěn)定性增長存在最不穩(wěn)定模數(shù)，最不穩(wěn)定模數(shù)隨屈服強(qiáng)度增加而減小，并近似與屈服強(qiáng)度的對數(shù)呈線性關(guān)系；隨著殼的厚度減小，擾動增長加快。

流體力學(xué)；界面不穩(wěn)定性；彈塑性；內(nèi)爆；材料參數(shù)

輕介質(zhì)支撐或者加速重介質(zhì)，介質(zhì)界面會發(fā)生Rayleigh-Taylor不穩(wěn)定性(RTI)，不穩(wěn)定性誘發(fā)界面擾動增長，重介質(zhì)形成尖釘侵入輕介質(zhì)，輕介質(zhì)氣泡進(jìn)入重介質(zhì)；而當(dāng)激波掠過介質(zhì)界面，會發(fā)生Richtmyer-Meshkov不穩(wěn)定性(RMI)，RMI的發(fā)展結(jié)果往往成為RTI發(fā)展的初始條件[1-2]。不穩(wěn)定性現(xiàn)象廣泛存在于天體物理、地球物理、慣性約束聚變(ICF)等領(lǐng)域。無黏流體的不穩(wěn)定性增長最快，存在多種制穩(wěn)機(jī)制，比如黏性、表面張力、密度梯度、材料強(qiáng)度、燒蝕等。對于彈塑性介質(zhì)，J.W.Miles[3]給出了無限厚彈性介質(zhì)RT不穩(wěn)定性增長的截止波長；D.C.Drucker[4]給出了彈性介質(zhì)RT不穩(wěn)定性增長的截止振幅；A.I.Lebedev等[5]給出了有限厚彈性介質(zhì)RT不穩(wěn)定性發(fā)展截止擾動振幅與波長、厚度之間的關(guān)系。J.F.Barnes等[6]通過爆轟加載鋁片實(shí)驗(yàn)，確認(rèn)了金屬界面RT不穩(wěn)定性的發(fā)展存在截止波長和振幅。J.W.Swegle等[7]對彈塑性介質(zhì)RTI進(jìn)行了充分的理論分析以及數(shù)值模擬工作，認(rèn)為彈塑性介質(zhì)界面不穩(wěn)定性發(fā)展不僅依賴于擾動波長、振幅、介質(zhì)厚度，還與壓力加載方式等因素相關(guān)。K.Mikaelian[8]提出黏性流體RMI的指數(shù)增長模型，J.D.Colvin等[9]將塑性介質(zhì)應(yīng)力與黏性流體類比，分析彈塑性介質(zhì)的界面不穩(wěn)定性，A.R.Piriz等[1-2]對彈塑性介質(zhì)RTI以及RMI進(jìn)行了有意義的理論解析探索。界面不穩(wěn)定性是實(shí)驗(yàn)研究材料彈塑性參數(shù)及其動態(tài)響應(yīng)特性的重要手段之一[10-11]。內(nèi)爆過程中界面不穩(wěn)定性是重要的研究方向，是影響ICF靶丸壓縮、點(diǎn)火狀態(tài)的重要因素之一。V.Frachet等[12]通過套筒實(shí)驗(yàn)，研究了柱形內(nèi)爆加載下流體以及金屬介質(zhì)界面不穩(wěn)定性，H.S.Park等[13]在更高壓力、應(yīng)變率下研究金了屬界面不穩(wěn)定性。劉軍等[14]對文獻(xiàn)[12]中的實(shí)驗(yàn)進(jìn)行了數(shù)值模擬。

本文中，使用自行研制的多介質(zhì)彈塑性流體力學(xué)歐拉程序(MEPH2Y)[15]，對柱形內(nèi)爆加載金屬界面不穩(wěn)定性進(jìn)行數(shù)值模擬，分析影響金屬界面不穩(wěn)定性發(fā)展的因素，探索金屬界面不穩(wěn)定的發(fā)展規(guī)律。

1 計算方法

1.1 計算模型

圖1 柱形內(nèi)爆計算模型Fig.1 Computational model in cylindrical implosion

采用與文獻(xiàn)[12]實(shí)驗(yàn)相似的幾何構(gòu)型(如圖1所示)，外層為鋼殼，向中心分別為高能炸藥、薄鋼殼、硅橡膠。鋼外殼，厚度2 cm；炸藥，外壁半徑10 cm，厚5 cm，從外壁面起爆爆炸藥；內(nèi)層鋼殼，厚4 mm，在內(nèi)層鋼殼的內(nèi)壁刻振幅0.5 mm的正弦曲線。

1.2 數(shù)值方法

在界面不穩(wěn)定性問題中存在界面大變形現(xiàn)象，因此，在數(shù)值模擬中采用自行研制的歐拉程序。所使用的多介質(zhì)彈塑性流體力學(xué)程序(MEPH2Y)采用交替方向方法進(jìn)行時間推進(jìn)，在每個方向上進(jìn)行拉氏計算，而后重映到原始網(wǎng)格，界面采用VOF方法重構(gòu)界面。該流體力學(xué)程序能夠模擬流體和彈塑性介質(zhì)，能夠用于炸藥爆轟、射流射彈形成及其侵徹過程等多介質(zhì)復(fù)雜流場模擬，具有較強(qiáng)的解決實(shí)際工程計算問題的應(yīng)用能力。

1.3 狀態(tài)方程及本構(gòu)模型

金屬采用凝聚介質(zhì)實(shí)用狀態(tài)方程：

式中：p為壓力，ρ為密度，c0為初始聲速，ρ0為初始密度，γ為比熱比。材料參數(shù)分別為：鋼材料，ρ0=7.8 g/cm3，c0=3.8 km/s，γ=4.075；硅橡膠，ρ0=1.1 g/cm3，c0=0.5 km/s，γ=1.6。爆轟計算中，炸藥使用CJ爆轟結(jié)合Wilkins反應(yīng)率的時間起爆方式，爆轟產(chǎn)物采用理想氣體狀態(tài)方程。炸藥參數(shù)為：初始密度1.8 g/cm3，爆速7.5 km/s，比熱比3.0。

材料采用理想彈塑性本構(gòu)模型，使用von Mises準(zhǔn)則確定介質(zhì)的彈塑性狀態(tài)，材料剪切模量G為常數(shù)，材料屈服強(qiáng)度Y滿足：Y=max(0,Y0(1-e/em))(Y0為初始屈服強(qiáng)度，em為融化能)。該本構(gòu)模型中，剪切模量與狀態(tài)無關(guān)，屈服強(qiáng)度僅與內(nèi)能相關(guān)，即僅有軟化的表現(xiàn)。材料本構(gòu)參數(shù)分別為：鋼殼，G=83 GPa，Y0=2 GPa，em=1.3 kJ/g；硅橡膠基本為流體模型，Y0=0.01 GPa，G=0.01 GPa，em=0.1 kJ/g；炸藥為流體。

2 結(jié)果及分析

2.1 一維數(shù)值模擬結(jié)果

圖2 金屬殼外壁壓力Fig.2 Pressure on outer interface of metal shell

圖3 金屬殼內(nèi)壁速度Fig.3 Velocity on inner interface of metal shell

通過對柱坐標(biāo)一維情況的數(shù)值分析，能夠認(rèn)識無擾動情況下界面運(yùn)動及受力情況，為進(jìn)一步分析界面不穩(wěn)定性發(fā)展提供參考。數(shù)值模擬中，采用0.1 mm均勻網(wǎng)格。圖2為鋼殼外壁的壓力，圖3為金屬殼內(nèi)的界面速度，時間原點(diǎn)為爆轟波到達(dá)鋼殼外壁的時刻。爆轟波在金屬外壁的壓力峰值約為45 GPa，波在金屬殼內(nèi)往返，能夠清晰地看到6 μs內(nèi)有3次沖擊波作用，殼層內(nèi)界面得到加速，界面擾動將發(fā)生RM不穩(wěn)定性。而后直到中心反射的沖擊波到達(dá)前，界面幾乎為勻減速運(yùn)動，加速度指向較重的鋼殼，界面擾動的不穩(wěn)定性增長將以RT不穩(wěn)定性為主。

2.2 二維數(shù)值模擬結(jié)果

在二維模擬中，金屬殼的初始擾動模數(shù)n取13和29。圖4為模擬得到的擾動模數(shù)為29的殼在某時刻的界面，從中可以定性地看到，材料強(qiáng)度對界面不穩(wěn)定性發(fā)展有較大影響。

圖5為不同初始擾動的擾動增長結(jié)果，包括純流體以及彈塑性計算。在純流體計算中，29模數(shù)的比13模數(shù)的擾動發(fā)展快，這與流體的擾動發(fā)展理論相符。對于相同模數(shù)擾動的界面，考慮彈塑性計算中的擾動比純流體計算中的擾動發(fā)展得慢，印證了材料強(qiáng)度是不穩(wěn)定性發(fā)展的致穩(wěn)因素。在彈塑性計算中，29模數(shù)的擾動未得到發(fā)展，而13模數(shù)的擾動發(fā)展雖然比純流體計算慢，但界面擾動仍然得到一定發(fā)展。這個彈塑性計算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)定性相符，即流體計算中發(fā)展較快的高模數(shù)在彈塑性計算中反而發(fā)展得較慢，表明金屬材料的不穩(wěn)定性發(fā)展與擾動模數(shù)及材料強(qiáng)度均相關(guān)。另外，彈塑性計算的結(jié)果與實(shí)驗(yàn)值有一定的差異，應(yīng)該與計算中所使用的狀態(tài)方程、本構(gòu)關(guān)系等材料特性的差異相關(guān)。

圖4 密度等值云圖(n=29)Fig.4 Density contour plot

圖5 金屬界面擾動增長Fig.5 Growth of metal interfacial perturbation amplitude

2.3 材料強(qiáng)度與不穩(wěn)定性的發(fā)展

圖6 材料屈服強(qiáng)度對界面擾動增長的影響Fig.6 Effect of yield strength on perturbation growth

圖7 材料剪切模量對界面擾動增長的影響Fig.7 Effect of shear strength on perturbation growth

固定材料剪切模量和屈服強(qiáng)度中的一個，單獨(dú)調(diào)節(jié)另一個。為了排除其他因素的干擾，材料本構(gòu)中屈服強(qiáng)度保持不變，即Y=Y0。圖6給出模數(shù)13、29在不同屈服強(qiáng)度下的擾動發(fā)展情況，可以看到，材料強(qiáng)度抑制不穩(wěn)定性發(fā)展，屈服強(qiáng)度越大擾動發(fā)展越慢。對不同模數(shù)的擾動，屈服強(qiáng)度的影響不同：對于13模擾動，材料強(qiáng)度對擾動發(fā)展影響相對較小；而對29模擾動發(fā)展影響顯著，當(dāng)Y0=2.0 GPa時，擾動甚至沒有發(fā)展。而圖7給出模數(shù)13當(dāng)Y0=1.5 GPa時在不同剪切模量下的擾動發(fā)展情況，可以看出，與屈服強(qiáng)度相比較，剪切模量在10～1 000 GPa這樣大范圍內(nèi)對不穩(wěn)定性發(fā)展影響甚微。實(shí)際上，從材料屈服強(qiáng)度以及剪切模量在控制方程中的角色，可以理解他們對于材料所表現(xiàn)出的強(qiáng)度的貢獻(xiàn)。屈服強(qiáng)度控制材料等效偏應(yīng)力張量，而剪切模量則僅控制材料偏應(yīng)力張量的增長速率，因此，在較高的應(yīng)變率加載下，材料的屈服強(qiáng)度決定了材料所體現(xiàn)出來的強(qiáng)度。另外，在擾動發(fā)展初期的RM不穩(wěn)定性發(fā)展階段，內(nèi)部壓強(qiáng)較高，材料強(qiáng)度相對較小，擾動增長受到材料強(qiáng)度的影響較弱，材料強(qiáng)度的作用主要體現(xiàn)在后期的RT不穩(wěn)定性主導(dǎo)階段。

從圖6可知，屈服強(qiáng)度對不同模數(shù)的擾動影響不同。這里，以擾動模數(shù)為變量分析不穩(wěn)定性發(fā)展的機(jī)制，并在給定的屈服強(qiáng)度下，將后期發(fā)展最快的擾動模數(shù)稱為最不穩(wěn)定模數(shù)。圖8展示了不同屈服強(qiáng)度下不同模數(shù)的擾動隨時間增長情況，所考察的材料屈服強(qiáng)度Y0=0～2.0 GPa、擾動模數(shù)n=8～40。計算結(jié)果顯示：純流體計算(見圖8(a))中，最不穩(wěn)性模數(shù)為最大擾動模數(shù)；中等程度的屈服強(qiáng)度(見圖8(b)～(c))下，隨著擾動模數(shù)增高，擾動發(fā)展速率先加快后減慢；最大的屈服強(qiáng)度(見圖8(d))下，最不穩(wěn)定模數(shù)為最小的模數(shù)。圖9為最不穩(wěn)定模數(shù)隨材料屈服強(qiáng)度的變化曲線，隨著屈服強(qiáng)度的增大，最不穩(wěn)定模數(shù)單調(diào)減小，如圖9中小圖所示，最不穩(wěn)定模數(shù)與屈服強(qiáng)度的對數(shù)近似呈線性關(guān)系：nmax≈12.2-17.0 lg(Y0/GPa)。實(shí)際上，如圖8所示，在擾動發(fā)展初期的RM不穩(wěn)定性發(fā)展階段，受材料強(qiáng)度影響較弱，隨擾動模數(shù)的增加擾動增長速率均加快。在后期RT不穩(wěn)定性發(fā)展階段，界面減速是失穩(wěn)因素，在該因素作用下高模數(shù)擾動增長較快；而材料強(qiáng)度是致穩(wěn)因素，且強(qiáng)度愈大對高模數(shù)擾動增長的抑制作用愈強(qiáng)；兩個相反因素作用的結(jié)果，使得在給定強(qiáng)度下擾動發(fā)展存在最不穩(wěn)定模數(shù)。這個最不穩(wěn)定模數(shù)下擾動增長最快，也可以理解為最危險的模數(shù)，或者說是該條件下的特征模數(shù)。

圖8 不同屈服強(qiáng)度下擾動模數(shù)與擾動增長的關(guān)系Fig.8 Effect of perturbation mode on perturbation growth under different yield strength

圖9 最不穩(wěn)定模數(shù)與屈服強(qiáng)度的關(guān)系Fig.9 Relationship between yield strength and most instable mode

2.4 金屬殼厚度對不穩(wěn)定性發(fā)展的影響

將金屬殼厚度4與2.5 mm的擾動發(fā)展情況進(jìn)行了對比，如圖10所示。相同擾動模數(shù)的薄殼擾動發(fā)展均快于厚殼， 在這里，以Y0=0.5 GPa為例，其他強(qiáng)度下擾動發(fā)展情況類似。另外，如圖11所示，2.5 mm厚的最不穩(wěn)定模數(shù)亦隨屈服強(qiáng)度的對數(shù)線性減小，nmax≈12.1-11.8 lg(Y0/GPa)，其斜率絕對值比4.0 mm厚的小，即較薄的殼其最不穩(wěn)定模數(shù)隨強(qiáng)度的增長變化較緩。

圖10 不同金屬殼厚度下的擾動增長Fig.10 Effect of shell thickness on perturbation growth

圖11 不同金屬殼厚度下最不穩(wěn)定模數(shù)與屈服強(qiáng)度的關(guān)系Fig.11 Relationship between yield strength and most instable mode with different shell thicknesses

3 結(jié) 論

利用自行研制的彈塑性流體力學(xué)的歐拉程序，對柱形內(nèi)爆加載下金屬界面不穩(wěn)定開展了數(shù)值研究，結(jié)果表明本文的計算結(jié)果與文獻(xiàn)實(shí)驗(yàn)結(jié)果吻合較好，進(jìn)而研究了材料參數(shù)對于內(nèi)爆加載金屬界面不穩(wěn)定發(fā)展的影響規(guī)律。

對于內(nèi)爆加載下金屬界面不穩(wěn)定問題，殼體擾動模數(shù)、材料強(qiáng)度對界面不穩(wěn)定性發(fā)展均有影響。對于純流體，擾動模數(shù)越高不穩(wěn)定發(fā)展越快；材料強(qiáng)度對界面不穩(wěn)定性發(fā)展有不可忽略的抑制作用，材料屈服強(qiáng)度對高模數(shù)不穩(wěn)定增長的抑制作用較大，而材料剪切模量對不穩(wěn)定性發(fā)展的影響相對較弱；對于相同擾動的金屬殼，存在最不穩(wěn)定模數(shù)，且該模數(shù)隨材料屈服強(qiáng)度的增加而減小，近似與屈服強(qiáng)度的對數(shù)呈線性關(guān)系；隨著殼的厚度減小，擾動增長加快。

感謝秦承森研究員的建設(shè)性討論。

[1] Piriz A R, Lopez Cela J J, Tahir N A, et al. Richtmyer-Meshkov instability in elastic-plastic media[J]. Physical Review E, 2008,78(5):056401.

[2] Piriz A R, Lopez Cela J J, Tahir N A. Linear analysis of incompressible Rayleigh-Taylor instability in solids[J]. Physical Review E, 2009,80(4):046305.

[3] Miles J W. Taylor instability of a flat plate[R]. San Diego: General Dynamics, GAMD-7335, 1966.

[4] Drucker D C. “Taylor instability” of the surface of an elastic-plastic plate[C]∥Nemat-Nasser C. Mechanics Today, Vol. 5. New York: Pergamon Press, 1980:37-47.

[5] Lebedev A I, Nizovtsev P N, Rayevsky V A. Rayleigh-Taylor instability in solids[C]∥4th International Workshop on the Physics of Compressible Turbulent Mixing (IWPCTM). Cambridge, England, 1993.

[6] Barnes J F, Blewett P J, McQueen R G, et al. Taylor instability in solids[J]. Journal of Applied Physics, 1974,45(2):727-732.

[7] Swegle J W, Robinson A C. Acceleration instability in elastic-plastic solids. I. Numerical simulations of plate acceleration[J]. Journal of Applied Physics, 1989,66(7):2838-2858.

[8] Mikaelian K. Effect of viscosity on Rayleigh-Taylor and Richtmyer-Meshkov instabilities[J]. Physical Review E, 1993,47(1):375-383.

[9] Colvin J D, Legrand M, Remington B A, et al. A model for instability growth in accelerated solid metals[J]. Journal of Applied Physics, 2003,93(9):5287-5301.

[10] Atchison W L, Zocher M A, Kaul A M. Studies of material constitutive behavior using perturbation growth in explosive and magnetically driven liner systems[J]. Russian Journal of Physical Chemistry B, 2008,2(3):387-401.

[11] Dimonte G, Terrones G, Cheme F J, et al. Use of the Richtmyer-Meshkov instability to infer yield stress at high-energy density[J]. Physical Review Letters, 2011,107(26):264502.

[12] Frachet V, Geleznikoff F, Guix R, et al. Rayleigh Taylor instability in cylindrical configuration[C]∥2nd International Workshop on the Physics of Compressible Turbulent Mixing (IWPCTM). Pleasanton, CA, 1989.

[13] Park H S, Lorenz K T, Cavallo R M, et al. Viscous Rayleigh-Taylor instability experiments at high pressure and strain rate[J]. Physical Review Letters, 2010,104(13):135504.

[14] 劉軍，馮其京，周海兵.柱面內(nèi)爆驅(qū)動金屬界面不穩(wěn)定性的數(shù)值模擬研究[J].物理學(xué)報,2014,63(15):155201. Liu Jun, Feng Qijing, Zhou Haibing. Simulation study of interface instability in metals driven by cylindrical implosion[J]. Acta Physica Sinica, 2014,63(15):155201.

[15] 馮其京，郝鵬程，杭義洪，等.聚能裝藥的歐拉數(shù)值模擬[J].爆炸與沖擊,2008,28(2):138-143. Feng Qijing, Hao Pengcheng, Hang Yihong, et al. Eulerian numerical simulation of a shaped charge[J]. Explosion and Shock Waves, 2008,28(2):138-143.

(責(zé)任編輯 丁 峰)

A numerical study of the instability of the metal shell in the implosion

Hao Pengcheng1,2,3, Feng Qijing1, Hu Xiaomian1,2

(1.InstituteofAppliedPhysicsandComputationalMathematics,Beijing100094,China;2.NationalKeyLaboratoryofComputationalPhysics,InstituteofAppliedPhysicsandComputationalMathematics,Beijing100088,China;3.GraduateSchool,ChinaAcademyofEngineeringPhysics,Beijing100088,China)

In this work the instability of the metal shell in the cylindrical implosion was studied numerically using a multi-component elastic-plastic hydrodynamic Eulerian code. Agreeing with those of the experiments, the numerical results show that the material strength restrains the growth of the interfacial perturbation with an effect not to be overlooked. The material yield strength has an obvious restraining effect on the higher mode of perturbation, while that of the shear module is similar but less sensitive. There exists a most instable mode number, which decreases as the yield strength increases and is approximately linear with the logarithm of the yield strength, and the perturbation grows faster as the shell grows thinner.

fluid mechanics; interfacial instability; elastic-plastic simulation; implosion; material parameters

10.11883/1001-1455(2016)06-0739-06

2015-04-30； < class="emphasis_bold">修回日期：2015-08-17

2015-08-17

國家自然科學(xué)基金項(xiàng)目(11402029,11372052,11371069,U1430235)； 中國工程物理研究院科學(xué)技術(shù)發(fā)展基金項(xiàng)目(2015B0101021)

郝鵬程(1980— )，男，博士研究生，副研究員,hao_pengcheng@iapcm.ac.cn。

O357.41 <國標(biāo)學(xué)科代碼：1302531 class="emphasis_bold"> 國標(biāo)學(xué)科代碼：1302531 文獻(xiàn)標(biāo)志碼：A國標(biāo)學(xué)科代碼：1302531

A