爆轟加載下彈塑性固體Richtmyer-Meshkov流動的擾動增長規律?

殷建偉 潘昊 吳子輝 郝鵬程 胡曉棉

1)(北京理工大學機電學院,北京 100081)

2)(北京應用物理與計算數學研究所,計算物理重點實驗室,北京 100094)

3)(中國工程物理研究院研究生部,北京 100088)

1 引 言

沖擊波加載不同密度的流體界面時,界面擾動將增長演化為渦結構并最終形成湍流混合層,這種流體動力學現象稱為Richtymer-Meshkov(RM)流動不穩定性[1?4],大量研究表明,流體的表面張力、黏性等耗散性質將削弱RM不穩定性的增長.對于固體介質,其RM流動過程受到固體物性的顯著影響,呈現出不同的物理現象:在沖擊波作用下,材料的強度致穩效應將導致擾動振幅增長一定程度后圍繞某均值周期性振蕩[5?9],若沖擊波壓力極高,如慣性約束聚變(ICF)內爆中加載靶丸的壓力達到100 GPa量級[10],擾動將劇烈增長,伴隨升溫、破碎、噴射、多物質混合等復雜物理現象,是武器內爆動力學、物質界面微噴射與混合動力學過程重點關注的問題.

目前對固體中RM流動的理論研究基于彈塑性材料均勻性假設,未考慮材料的相變、斷裂等復雜物理因素.Piriz[8]提出了剛塑性材料擾動面上尖釘振幅的極大值表達式,并在Buttler等[11,12]和Dimonte等[13]的實驗中得到初步驗證.Jensen等[14]在平面撞擊實驗中利用該表達式推算了鈰在極端加載條件下的強度,結果表明不同方法得到的強度結果偏差約20%—50%.因此,如何改進極值振幅表達式,更準確地描述彈塑性固體RM流動的擾動增長規律是目前亟待解決的問題,也是理解強度致穩機制的核心問題,有助于提高對界面失穩早期階段的物理認識.

本文采用AFE2D程序對炸藥爆轟加載下金屬材料的界面擾動增長過程進行數值模擬研究,定性分析了初始擾動的振幅與波長對擾動增長的影響以及材料強度的致穩效應,利用相關沖擊波關系式及流體RM不穩定性的研究成果,進一步推導得到了極值振幅的理論關系式,該式能夠描述擾動增長被抑制時振幅與初始擾動、材料性質之間的關系,通過分析關系式中各項的相關性,初步得到了線性關系式的適用范圍.

2 界面擾動增長問題的數值模擬

考慮如圖1所示的界面擾動問題,高能炸藥爆轟驅動金屬樣品,該問題是Buttler等[11]開展的銅、錫等材料RM實驗的簡化模型,其中高能炸藥厚度12 mm,金屬樣品厚6 mm,樣品中的沖擊波向介質/真空界面傳播時加載擾動自由面,擾動面為正弦形,其初始波長為λ0,初始振幅為ξ0.

圖1 沖擊波加載界面擾動問題示意圖Fig.1.Illustration of the surface perturbation problem loaded by the shock wave.

炸藥爆轟采用時間起爆方式,狀態方程為JWL類型狀態方程[15]

其中R1,R2和ω是JWL模型參數;A和B是JWL模型系數,由高能炸藥的CJ爆速DCJ、CJ爆壓PCJ及JWL參數確定;V=ρ/ρ0是相對比容,ρ0是初始密度,ρ是材料密度;E是單位初始體積的內能.固體材料的狀態方程為Mie-Grüneisen狀態方程[16]

其中μ=ρ/ρ0?1;γ0為Grüneisen系數;b為Grüneisen系數的一階體積修正量;c0,S1,S2和S3為材料沖擊波-粒子速度關系式的系數.固體材料本構采用Prandtl-Reuss關聯塑性流動律和von Mises屈服準則描述[7?9],即

下標采用Einstein求和約定,Dij為應變率張量;Sij為偏應力張量;G為材料剪切模量,Y為屈服強度,采用彈性理想塑性模型,保持G和Y不變.幾種常見金屬材料的狀態方程參數與剪切模量見表1,炸藥的爆轟參數與JWL狀態方程參數見表2.

利用AFE2D程序數值模擬PBX9501炸藥爆轟驅動銅飛片的實驗(見文獻[11]Cu-pRad0426實驗),如圖1所示,樣品的自由面上引入kξ0(k=2π/λ0,初始擾動的波數)分別為0.12,0.35,0.75和1.5的初始擾動,初始波長均為550μm,銅樣品中的沖擊波壓力約36 GPa.圖2為實驗的數值模擬結果,可以看到沖擊波加載擾動面時,初始波谷位置首先啟動,形成向外的突起尖釘(Spike)結構;初始波峰位置稍后啟動,與尖釘等部位間的速度差導致其向材料內部凹陷,形成氣泡(Bubble)結構;反相后尖釘與氣泡結構進一步演化,初始擾動kξ0=0.12時尖釘基本不增長,kξ0=0.35時尖釘增長一段時間后被抑制,kξ0=0.75和1.5時增長不受抑制,尖釘演化為射流.

表1 金屬材料的狀態方程參數和剪切模量Table 1.Parameters of equation of state and shear modulus of several selected metals.

表2 PBX9501炸藥的爆轟參數與JWL狀態方程參數Table 2.Parameters of the detonation performance and JWL equation of state of the high explosives PBX9501.

定義尖釘的振幅ξSP為尖釘與無擾動區域(Bulk)之間在沖擊波傳播方向的相對位移,尖釘的增長速度vSP為尖釘與無擾動區域之間在沖擊波傳播方向的相對速度,即

氣泡的振幅ξBU和增長速度vBU的定義類似,

模擬結果顯示,kξ0=0.12時尖釘的最大振幅=max(ξSP)=10.55 μm,相比初始振幅=11μm幾乎沒有增長,但相位出現反轉;kξ0=0.35時=155.3 μm. 在Buttler對CupRad0426實驗數據的分析中[11],kξ0=0.12時=(7±4)μm,kξ0=0.35時=160 μm.表3對比了不同kξ0時尖釘增長速度的極大值=max(vSP),模擬結果與實驗數據基本一致,表明本文采用的程序、模型和參數適用于界面擾動問題,可用于分析界面擾動的增長規律.

表3 PBX9501炸藥驅動銅飛片實驗的模擬結果Table 3.Simulation results of the Cu-pRad0426 experiment where the copper fl yer plate was driven by the high explosive PBX9501.

圖2 Cu-pRad0426實驗的模擬結果 (a)初始擾動界面;(b)沖擊波首次加載界面后0.7μs時刻的模擬結果;(c)3.7μs時刻的模擬結果;(d)3.7μs時刻實驗質子照相圖像[11]Fig.2.The simulation results of the Cu-pRad0426 experiment:(a)Initial state;(b)plot of the simulation results at 0.7μs which relative to shockwave breakout at the free surface;(c)plot of the simulation results at 3.7μs;(d)proton radiographic picture of the experiment at 3.7μs[11].

3 界面擾動的增長規律

3.1 初始形態對擾動增長的影響

首先分析擾動初始形態對后期增長的影響,利用無量綱數kξ0表征初始擾動的振幅與波長之比,考慮圖1所示的界面擾動問題,選取初始波長λ0=450—600 μm,初始振幅ξ0=15—38 μm,則kξ0的變化范圍0.20—0.40,樣品材料為OFHC銅,屈服強度Y=0.47 GPa,炸藥為PBX9501.圖3比較了沖擊波首次加載界面后3.1μs時刻不同kξ0的擾動增長情況,尖釘的振幅隨著kξ0增加逐步變大,尖釘形狀在kξ0=0.20時仍接近正弦擾動形態,繼續提高kξ0尖釘將逐漸拉伸、變窄,最終演化為射流.上述結果表明,對于相同的加載條件和樣品材料,若初始擾動的振幅與波長比值越小,相應的初始界面越“平坦”,擾動增長越不明顯,界面保持穩定;振幅與波長比值越大,相應的初始界面越“溝壑”,擾動越易于增長,甚至形成局部的微射流,界面失穩.

圖3 不同kξ0值的界面擾動增長情況,沖擊波首次加載界面后3.1μs時刻Fig.3.Growth results of the perturbated surfaces with different values of kξ0,3.1 μs relative to shockwave breakout at the free surface.

3.2 材料強度對擾動增長的致穩效應

相關研究表明,金屬材料的動態屈服強度存在較強的應變硬化效應,并且在高應變率下表現出顯著的應變率硬化效應[17?21].本文所考慮的沖擊波加載擾動面問題中,擾動面上形成的尖釘和氣泡等特征結構均為局部高應變區,局部應變率達到106—107s?1.然而目前高應變率范圍內金屬材料的強度數據極少,對應變率的硬化規律認識亦有待進一步完善,本文所采用的彈性理想塑性模型雖不能反映材料真實的應變與應變率硬化效應,但在均勻性假設下通過調節屈服強度模擬不同的硬化效果也不失為一種可行的近似方法,能夠用于近似闡述材料強度的致穩效應.

選擇擾動的kξ0=0.251,取OFHC銅的屈服強度為Y=0.17 GPa,令其逐步增加至0.57 GPa以模擬不同的硬化效應.圖4為擾動界面的增長情況,強度較低時擾動面的尖釘演化為射流,提高材料強度后尖釘增長受到抑制,逐漸從尖銳的射流形態退化為類似正弦波形的擾動形態,定性地顯示出材料強度對擾動增長的致穩效應,表明在相同的加載和初始擾動條件下,強度越高的材料擾動增長越不明顯.

3.3 擾動增長規律的線性近似

由圖3和圖4可知,提高初始擾動的kξ0,擾動面的尖釘振幅將逐漸增大,在爆轟加載下擾動增長受到材料強度的顯著抑制作用,因此當擾動增長被抑制時,尖釘增長幅度的控制因素應包含初始擾動kξ0與強度Y的某種組合形式.

圖4 擾動增長情況,調節屈服強度以模擬不同程度的硬化效應,沖擊波首次加載界面后3.1μs時刻Fig.4. Growth results of the surfaces with varied yield strength for the different equivalent effects of hardening,3.1μs relative to shockwave breakout at the free surface.

利用(4)式定義的尖釘振幅、增長速度和Dimonte等[7,13]的分析方法,類似于尖釘最大振幅和最大增長速度的定義,得到氣泡最大振幅=max(ξBU)和氣泡最大增長速度=max(vBU),利用初始波數無量綱化尖釘和氣泡的最大振幅

利用材料初始密度ρ0、屈服強度Y無量綱化尖釘和氣泡的最大增長速度,分別標記為無量綱數和

圖5 (網刊彩色)尖釘和氣泡的最大振幅與最大增長速度的關系(銅強度分別取0.27,0.32和0.37 GPa)Fig.5. (color online)Relations between the nondimensional maximum amplitude and growth velocity of the spikes and bubbles(yield strength of copper varied as 0.27,0.32 and 0.37 GPa).

圖5和圖6分別顯示了銅、鋁樣品中尖釘和氣泡的最大振幅與最大增長速度無量綱數之間的關系,選擇初始擾動kξ0變化范圍0.15—0.3以避免尖釘演化為射流,屈服強度分別取三組值.模擬結果顯示,無量綱化的尖釘最大振幅與最大增長速度在一定范圍內表現出近似線性關系,斜率C約為0.3,即

圖6 (網刊彩色)尖釘和氣泡的最大振幅與最大增長速度關系(鋁強度分別取0.24,0.29和0.34 GPa)Fig.6. (color online)Relations between the nondimensional maximum amplitude and growth velocity of the spikes and bubbles(yield strength of aluminum varied as 0.24,0.29 and 0.34 GPa).

Piriz理論預測的斜率值為0.29[7],Dimonte利用PAGOSA程序擬合的數據點見圖5[13],忽略截距后斜率值為0.24,本文結果中線性擬合通過原點,斜率約為0.3,更接近Piriz的理論預測值.

隨著擾動kξ0增大,尖釘的增長速度逐步增加,最大振幅與增長速度之間表現出非線性關系,尖釘形態從擾動增長被抑制時的穩定態向持續增長的不穩定態過渡,直至尖釘頭部收縮變窄形成射流,最終脫離樣品主體形成噴射物.流體RM不穩定性分析研究表明[1,2,10,22],在線性近似假設下尖釘的最大增長速度可近似表示為初始擾動kξ0、沖擊波經過后的界面起跳速度Δv和界面兩側Atwood數的乘積,即

其中Atwood數A=(ρ2?ρ1)/(ρ2+ρ1)表示界面兩側物質的密度差,對于固體/真空界面,A=?1.沖擊波在彈塑性材料中傳播時,界面起跳速度Δv可由沖擊波強度和材料的Hugoniot沖擊線確定[11,17,23?25],若介質的初始壓力和速度為零,壓力為P的沖擊波傳播時波后粒子速度與壓力的關系為

其中η表示沖擊波前后介質的壓縮率,η=1?ρ0/ρ,與材料的聲阻抗和沖擊波速度有關;vp為波后粒子速度.對于介質/真空界面,自由面速度與波后粒子速度存在倍數關系,Δv≈2vp,將(9)和(10)式代入(8)式中得

因此尖釘振幅的增長因子為

(12)式表明,擾動增長被抑制時尖釘的最大振幅增長因子與參數組合kξ0/Y成線性關系,并且沖擊波壓力越高,尖釘增長越顯著.

3.4 線性近似的適用范圍

分析(12)式中制約尖釘增長因子的各因素:系數C是線性關系的斜率;壓縮率η與沖擊波壓力P共同表征了沖擊波的壓縮效應,為擾動增長提供啟動能量;kξ0描述了擾動界面的初始形態,為初始擾動的振幅與波長之比;屈服強度Y是材料的物性參數,擾動增長受其抑制.因此后三個因素之間是相互獨立的,本文分析系數C與其他因素之間的相關性,討論線性近似的適用范圍.

考慮三種材料,分別是6061鋁、OFHC銅和304不銹鋼,狀態方程為Mie-Grüneisen狀態方程,本構模型為彈性理想塑性模型,參數見表1,炸藥為PBX9501,相關模型參數見表2.每種材料的屈服強度取三組值以模擬不同程度的硬化效應,初始擾動kξ0的變化范圍0.15—0.35.圖7為不同材料、強度和kξ0時系數C隨參數組合kξ0/Y的變化情況,圖8為kξ0=0.223和0.279時鋁材料擾動面上尖釘、氣泡和無擾動區域的速度曲線,分別對應擾動增長被抑制和未被抑制的情況.

圖7 (網刊彩色)無量綱的尖釘最大振幅與增長速度之比隨參數組合kξ0/Y 的變化情況Fig.7.(color online)Plots of the ratios between the non-dimensional maximum amplitude and growth velocity of the spikes verse the combination of parameters kξ0/Y.

由圖7可知,系數分布存在三個階段,對應于界面擾動的穩定、穩定與不穩定之間的過渡和不穩定階段.

1)參數組合kξ0/Y≤0.8 GPa?1,系數C基本為常數,擾動處于被抑制的穩定態,尖釘的最大振幅與增長速度之間為線性關系,增長因子由(12)式描述,并且圖7中不同材料的系數一致,說明在穩定段系數C基本保持獨立性,不受其它因素的影響.

2)參數組合0.8 GPa?1≤kξ0/Y≤1.0 GPa?1,系數C不再為常數,擾動增長情況由穩定態向不穩定態過渡.此階段擾動可能會形成穩定的尖釘形態,也可能進一步發展成射流,取決于兩個特征時間的競爭,即沖擊波首次加載后在樣品內部波系多次反射形成的二次沖擊波加載界面所需時間twave,與沖擊波首次加載后尖釘振幅從開始增長至被抑制所需的穩定時間tstable.若twave＞tstable,表明在二次沖擊波到達前,強度的致穩效應已抑制了擾動增長,尖釘、氣泡和未擾動區域的速度趨于一致,如圖8所示,從被抑制到二次沖擊波到達,這段時間內尖釘、氣泡和無擾動區域之間的相對速度保持為零,意味著尖釘、氣泡的振幅不再變化,界面形狀保持不變;若twave≤tstable,表明在二次沖擊波到達自由面時,強度的致穩效應未能將尖釘充分減速,尖釘與其他位置相比存在速度差,沖擊波二次加載將破壞此時的尖釘和氣泡形態,界面形狀進一步演化.在過渡階段,(8)式的線性近似關系與(12)式的尖釘增長規律依然可能成立,取決于擾動增長在二次沖擊波達到前是否能被抑制.

圖8 (網刊彩色)擾動增長穩定與不穩定時尖釘、氣泡和未擾動區域的速度曲線Fig.8.(color online)Velocity pro files of the spikes,bubbles and the unperturbated areas in the stable and unstable growth of the perturbations.

3)參數組合kξ0/Y≥1.0 GPa?1,擾動增長顯著,尖釘頭部迅速演化為射流,界面失穩.

綜上所述,在擾動增長的穩定態,系數C保持獨立性,線性關系(8)式適用,描述尖釘增長規律的線性近似(12)式具備一定的材料普適性;在穩定態向不穩定態過渡階段以及不穩定態,系數C受其他因素影響失去獨立性,(8)式和(12)式不具備材料普適性;圖7中系數C保持為常數的平臺區域揭示了界面擾動從增長被抑制的線性關系區向形成射流的非線性區過渡時存在穩定性閾值,表明強度介質的RM擾動增長中也存在類似Rayleigh-Taylor擾動增長的穩定性閾值[26],初步研究結果表明閾值對應的物理狀態kξ0/Y取值約0.8 GPa?1.

4 結 論

本文研究了爆轟加載下不同材料飛片RM流動的擾動增長規律,分析結果表明擾動增長被抑制時,尖釘振幅增長因子與初始擾動的kξ0成正比,與材料的強度Y成反比,理論關系式符合相關實驗結果,與Piriz利用牛頓第二定律近似推導所得的結果基本一致,該理論關系式也為測量材料動態屈服強度提供了新的實驗技術途徑.對擾動增長規律影響因素的相關性分析表明,不同金屬材料的擾動增長均存在穩定區域,穩定區內描述尖釘振幅增長因子與加載條件、初始擾動和材料性質之間關系的(12)式具備一定的材料普適性,穩定區范圍kξ0/Y＜0.8 GPa?1,其中閾值0.8 GPa?1的物理意義有待進一步分析研究.

[1]Richtmyer R D 1960Commun.Pure Appl.Math.13297

[2]Meshkov E E 1969Sovit.Fluid Dyn.4151

[3]Mikaelian K O 2013Phys.Rev.E87031003

[4]Brouillette M 2002Annu.Rev.Fluid Mech.34445

[5]Piriz A R,Lopez Cela J J,Tahir N A,Hoff mann D H H 2006Phys.Rev.E74037301

[6]Piriz A R,Lopez Cela J J,Cortazar O D,Tahir N A,Hoff mann D H H 2005Phys.Rev.E72056313

[7]Piriz A R,Lopez Cela J J,Tahir N A,Hoff mann D H H 2008Phys.Rev.E78056401

[8]Piriz A R,Lopez Cela J J,Tahir N A 2009Phys.Rev.E80046305

[9]Lopez Ortega A,Lombardini M,Pullin D I,Meiron D I 2014Phys.Rev.E89033018

[10]Remington B A,Rudd E R,Wark J S 2015Phys.Plasmas22090501

[11]Buttler W T,Oro D M,Preston D L,Mikaelian K O,Cherne F J,Hixson R S,Mariam F G,Morris C,Stone J B,Terrones G,Tupa D 2012J.Fluid Mech.70360

[12]Buttler W T,Oro D M,Olsen R T,Cheren F J,Hammerberg J E,Hixson R S,Monfared S K,Pack C L,Rigg P A,Stone J B,Terrones G 2014J.Appl.Phys.116103519

[13]Dimonte G,Terrones G,Cheren F J,Germann T C,Dunpont V,Kadau K,Buttler W T,Oro D M,Morris C,Preston D L 2011Phys.Rev.Lett.107264502

[14]Jensen B J,Cheren F J,Prime M B,Fezzaa K,Iverson A J,Carlson C A,Yeager J D,Ramos K J,Hooks D E,Cooley J C,Dimonte G 2015J.Appl.Phys.118195903

[15]Sun Z F,Xu H,Li Q Z,Zhang C Y 2010Chin.J.High Pressure Phys.2455(in Chinese)[孫占峰,徐輝,李慶忠,張崇玉2010高壓物理學報2455]

[16]Robinson A C,Swegle J W,1989J.Appl.Phys.662838

[17]Zhu J S,Hu X M,Wang P,Chen J,Xu A G 2010Adv.Mech.40400(in Chinese)[朱建士,胡曉棉,王裴,陳軍,許愛國2010力學進展40400]

[18]Vogler T J,Chhabildas L C 2006Int.J.Impact Engng.33812

[19]Barton N R,Bernier J V,Becker R,Arsenlis A,Cavallo R,Marian J,Rhee M,Park H S,Remington B A,Olson R T 2011J.Appl.Phys.109073501

[20]Smith R F,Eggert J H,Rudd R E,Swift D C,Blome C A,Collins G W 2011J.Appl.Phys.110123515

[21]Park H S,Rudd R E,Cavallo R M,Barton N R,Arsenlis A,Belof J L,Blobaum K J M,El-dasher B S,Florando J N,Huntington C M,Maddox B R,May M J,Plechaty C,Prisbrey S T,Remington B A,Wallace R J,Wehrenberg C E,Wilson M J,Comley A J,Giraldex E,Nikroo A,Farrell M,Randall G,Gray III G T 2015Phys.Rev.Lett.114065502

[22]Wouchuk J G 2001Phys.Rev.E63056303

[23]Pan H,Wu Z H,Hu X M,Yang K 2013Chin.J.High Pressure Phys.27778(in Chinese)[潘昊,吳子輝,胡曉棉,楊堃2013高壓物理學報27778]

[24]Pan H,Hu X M,Wu Z H,Dai C D,Wu Q 2012Acta Phys.Sin.61206401(in Chinese)[潘昊,胡曉棉,吳子輝,戴誠達,吳強2012物理學報61206401]

[25]Yu Y Y,Tan H,Hu J B,Dai C D,Chen D N,Wang H R 2008Acta Phys.Sin.572352(in Chinese)[俞宇穎,譚華,胡建波,戴誠達,陳大年,王煥然 2008物理學報 572352]

[26]Colvin J D,Legrand M,Remington B A,Schurtz G,Weber S V 2003J.Appl.Phys.935287