凝聚相炸藥爆轟波沖擊空穴塌陷過程的高精度數(shù)值模擬

劉新橋， 王成

(北京理工大學(xué) 機(jī)電學(xué)院， 爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，北京 100081)

?

凝聚相炸藥爆轟波沖擊空穴塌陷過程的高精度數(shù)值模擬

劉新橋， 王成

(北京理工大學(xué) 機(jī)電學(xué)院， 爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，北京 100081)

利用高精度數(shù)值格式對凝聚相炸藥爆轟波沖擊空穴塌陷過程進(jìn)行了研究. 采用RGFM和Level-set相結(jié)合的物質(zhì)界面處理方法，克服了爆炸流場中高密度比、高壓強(qiáng)比的物質(zhì)界面容易引起數(shù)值振蕩這一問題. 采用5階WENO高精度有限差分格式進(jìn)行空間離散，采用3階TVD Runger-Kutta法進(jìn)行時(shí)間離散，開發(fā)了高精度動態(tài)并行程序，并運(yùn)用該程序?qū)δ巯嗾ㄋ幈Z波沖擊空穴后的塌陷過程進(jìn)行了數(shù)值模擬，給出了PBX-9404炸藥在不同入射壓力下的起爆距離及沖擊起爆過程中的壓力歷史曲線圖，并分別與實(shí)驗(yàn)結(jié)果對比后進(jìn)行了驗(yàn)證. 對凝聚相炸藥中的空穴塌陷過程進(jìn)行了高精度數(shù)值計(jì)算，給出了不同時(shí)刻的密度、壓力和粒子速度云圖，研究結(jié)果表明由于空穴的存在，爆轟波的傳播過程中會在空穴處形成射流現(xiàn)象以及密度和能量的匯聚，進(jìn)一步產(chǎn)生高溫高壓的熱點(diǎn).

凝聚相炸藥；高精度；數(shù)值模擬；空穴塌陷；熱點(diǎn)

凝聚相炸藥中的空穴塌陷不僅會導(dǎo)致熱點(diǎn)的形成和隨后的局部點(diǎn)火點(diǎn)的生成，在炸藥爆轟過程中，由于空穴的存著還會影響到炸藥的爆轟波傳播規(guī)律及爆轟性能. 因此，詳細(xì)了解空穴塌陷和熱點(diǎn)生成的背后物理機(jī)制有利于加深對爆轟現(xiàn)象的認(rèn)識，對于完善爆轟理論及炸藥裝藥的優(yōu)化設(shè)計(jì)起著十分重要的作用. 溫麗晶等[1]基于Kim的彈粘塑性單球殼塌縮模型，得到了新的沖擊起爆三項(xiàng)式細(xì)觀反應(yīng)速率模型，并把該模型加入DYNA2D中，模擬了PBX-9501炸藥的一維沖擊起爆過程. 劉源等[2]驗(yàn)證了Particle Level-set方法界面處理的高精度方法，將Particle Level-set方法與流體動力學(xué)控制方程進(jìn)行了耦合對激波與氣泡相互作用問題進(jìn)行了數(shù)值模擬. Udaykumar等[3]對HMX炸藥中的空穴塌陷過程進(jìn)行了數(shù)值模擬，研究了加載強(qiáng)度和孔隙大小的影響，結(jié)果表明，在微米尺度的空隙中，氣相化學(xué)反應(yīng)是一個(gè)重要因素，但對空隙崩潰本身的影響作用很小. 空隙完全塌陷時(shí)將會有大量反應(yīng)產(chǎn)物生成，由于空穴崩潰的時(shí)間太短，無法在崩潰之前建立能夠引發(fā)炸藥起爆的熱點(diǎn).

本文利用高精度數(shù)值格式對凝聚相炸藥沖擊空穴塌陷過程進(jìn)行了研究. 采用RGFM和Level-set相結(jié)合的物質(zhì)界面處理方法，在數(shù)值上采用5階WENO格式[4]和3階TVD Runge-Kutta法對描述炸藥起爆和爆轟波傳播的Euler方程組進(jìn)行離散，并開發(fā)了高精度動態(tài)并行程序. 用所開發(fā)的程序?qū)φㄋ幍目昭ㄋ葸^程進(jìn)行了模擬，給出了由氣泡存在而引起的局部射流現(xiàn)象的過程以及氣泡存在對爆轟波傳播規(guī)律的影響.

1 化學(xué)反應(yīng)流體動力學(xué)控制方程組

由于爆轟過程的時(shí)間極短，因而可以忽略爆轟過程中的黏性、熱傳導(dǎo)和擴(kuò)散等輸運(yùn)效應(yīng)，則非定常可壓縮Euler方程組可以作為爆轟過程中流體動力學(xué)控制方程組. 其形式為

(1)

(2)

(3)

(4)

(5)

(6)

式中：ρ為密度；p為壓強(qiáng)；u，v分別為x，y方向的速度；E為單位質(zhì)量總能量；e為比內(nèi)能.

狀態(tài)方程采用JWL狀態(tài)方程，反應(yīng)物和產(chǎn)物的狀態(tài)方程參數(shù)值的表達(dá)式如式(7)所示，對于PBX-9404炸藥，狀態(tài)方程中參數(shù)見表1[5].

(7)

式中：ν為比容；ν0為初始時(shí)刻炸藥的比容；T為溫度；A，B，R1，R2，Cυ和ω為常數(shù).

表1 PBX-9404炸藥JWL狀態(tài)方程參數(shù)

對于凝聚相爆轟的化學(xué)反應(yīng)率模型，使用較多的是Lee和Tarver提出的點(diǎn)火和成長反應(yīng)率模型， 對于炸藥PBX-9404，反應(yīng)速率模型中的參數(shù)見表2[5].

表2 PBX-9404炸藥Lee-Tarver模型參數(shù)

(8)

式中：I，G，Z均為常數(shù)；ν1為受沖擊但未反應(yīng)炸藥的比容；η為未反應(yīng)炸藥的相對壓縮度.

2 界面處理方法

2.1 Level-set界面追蹤

本文采用Level-set物質(zhì)界面跟蹤方法，通過構(gòu)造界面函數(shù)φ(x，y，t)，使得在任意時(shí)刻，運(yùn)動的物質(zhì)界面Γ(t)正好是φ(x，y，t)的零等值面，即

(9)

φ函數(shù)選取的初值應(yīng)滿足在Γ(0)附近為法向單調(diào)，而且在Γ(t)上為0. 一般取φ(x，y， 0)為點(diǎn)(x，y)到界面Γ(0)的符號距離，即

(10)

式中d((x，y)，Γ(0))為(x，y)到Γ(0)的距離. 因此，Level-set方法是一種隱式的界面追蹤方法. 隱式是指在任意時(shí)刻，僅僅給定計(jì)算節(jié)點(diǎn)到界面的法向符號距離，通過各個(gè)節(jié)點(diǎn)到界面的距離來隱式地刻畫出界面位置所在. 通常情況下，在任意時(shí)刻，所有計(jì)算節(jié)點(diǎn)的φ函數(shù)值都不會恰好為0，也即不存在任意一個(gè)節(jié)點(diǎn)正好處于運(yùn)動物質(zhì)界面上.

2.2 RGFM界面處理

RGFM方法首先在物質(zhì)界面處構(gòu)建局部Riemann問題，并將局部Riemann問題的求解結(jié)果賦予真實(shí)流體一側(cè)緊鄰物質(zhì)界面的帶狀區(qū)域，然后通過常值法向外推，在該介質(zhì)以外的流場賦予該介質(zhì)的Ghost物理量，將多介質(zhì)問題轉(zhuǎn)化為若干單介質(zhì)問題進(jìn)行處理. 通過求解局部Riemann問題，可以確保物質(zhì)界面兩側(cè)的壓強(qiáng)和法向速度始終保持一致，而密度和切向速度等物理量則允許存在間斷，這與實(shí)際流場中物質(zhì)界面兩側(cè)的間斷條件一致. 求解Riemann問題就是求解物質(zhì)界面附近的流場相互作用之后的流動狀態(tài)，因此，RGFM邊界條件定義方法可以真實(shí)的反映介質(zhì)之間的接觸作用過程.

3 數(shù)值方法

采用有限差分WENO格式構(gòu)造方法對控制方程的空間進(jìn)行離散，需求出通量F(U)、G(U)相應(yīng)的特征矩陣及其左右特征向量，由于在化學(xué)反應(yīng)區(qū)由于同時(shí)存在爆轟產(chǎn)物和未反應(yīng)炸藥，壓力沒有具體的表達(dá)形式，需要將問題歸結(jié)為一個(gè)線性方程組的求解，從而得到在化學(xué)反應(yīng)區(qū)內(nèi)包含各守恒量的壓力表達(dá)形式[6]. 并通過左特征向量將WENO格式中用到的值投影到局部特征空間，實(shí)現(xiàn)WENO重構(gòu)后，將所得到的結(jié)果通過右特征向量轉(zhuǎn)化到物理空間. 有限差分格式為

(11)

時(shí)間上采用3階TVD Runge-Kutta方法進(jìn)行離散：

(12)

4 數(shù)值模擬結(jié)果與分析

4.1 PBX-9404炸藥沖擊起爆數(shù)值模擬與實(shí)驗(yàn)對比

本算例對沖擊起爆數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果[7-8]進(jìn)行了對比，初始參數(shù)如表3所示[5]. 圖1(a)給出了不同炸藥的入射壓力與起爆距離之間的關(guān)系，曲線6為炸藥PBX-9404的實(shí)驗(yàn)結(jié)果，本文計(jì)算結(jié)果在圖中通過標(biāo)記各點(diǎn)后連線得到. 圖1(b)給出了炸藥在2.5 GPa入射壓力下壓力隨時(shí)間變化曲線圖. 由圖中可以看到，本文的數(shù)值模擬結(jié)果與實(shí)驗(yàn)吻合較好.

表3 PBX-9404炸藥參數(shù)

4.2 PBX-9404炸藥空穴塌陷過程數(shù)值模擬

取長寬分別為15和10 mm的PBX-9404炸藥，空穴中心距離左端5 mm處，在炸藥左端為2 mm范圍內(nèi)采用CJ條件起爆. 對物質(zhì)界面采用Local Level Set方法進(jìn)行追蹤和Real Ghost Fluid Method(RGFM)的界面處理方法. 不同時(shí)刻的密度和壓力云圖如圖2所示.

由圖2可見，在爆轟波到達(dá)空穴的初始階段，由于空穴對爆轟波的稀疏作用，將使得爆轟波陣面的密度和壓力迅速降低并導(dǎo)致波速下降，但在爆轟波與空穴邊界接觸位置的密度和壓力將會產(chǎn)生匯聚，略有增加，此時(shí)在空穴內(nèi)部x方向的粒子速度將急劇增加. 隨著爆轟的向前傳播，在空穴內(nèi)部將形成射流，并且空穴被進(jìn)一步壓縮，在空穴壓縮塌陷過程中，會出現(xiàn)高溫、高壓和高密度的現(xiàn)象，導(dǎo)致熱點(diǎn)出現(xiàn)，當(dāng)爆轟波陣面完全掃過空穴時(shí)，空穴被壓縮為一個(gè)極小的高密度區(qū)域，此時(shí)，壓力和粒子速度最大的區(qū)域出現(xiàn)在爆轟波陣面繞過空穴后匯聚的位置.

圖3～圖5給出了3個(gè)不同時(shí)刻密度、壓力及x方向粒子速度曲線. 由密度圖可見，隨著爆轟波的推進(jìn)密度開始匯聚，并于0.78 μs在7.5 mm處達(dá)到最大值6.3 g/cm3，這已遠(yuǎn)大于穩(wěn)定爆轟時(shí)的密度峰值. 由壓力圖可以看到，在0.25 μs時(shí)壓力峰值約為54.5 GPa，而在0.46 μs時(shí)，氣泡左端壁面的壓力已小于CJ壓力，而在氣泡的右端由于匯聚出現(xiàn)了45.0 GPa的峰值. 在0.78 μs時(shí)出現(xiàn)兩個(gè)壓力峰值，氣泡內(nèi)部的壓力進(jìn)一步匯聚達(dá)到57.0 GPa，而此時(shí)在氣泡右端內(nèi)壁面處壓力已達(dá)到78.0 GPa. 由粒子速度曲線圖可見，隨著時(shí)間的推進(jìn)，粒子速度峰值由0.46 μs到0.78 μs急劇下降到4.5 km/s，同時(shí)在7.5 mm位置處粒子速度具有極小值，與壓力曲線圖對比可知此處的壓力具有極大值，這是由于粒子在空穴壁面反射速度反向后與向前運(yùn)動的粒子進(jìn)行了碰撞而導(dǎo)致的.

5 結(jié) 論

采用5階WENO高精度有限差分格式以及3階TVD Runger-Kutta方法分別在空間和時(shí)間上進(jìn)行離散，基于RGFM和Level-set相結(jié)合的物質(zhì)界面處理方法，克服了爆炸流場中高密度比、高壓強(qiáng)比物質(zhì)界面引起的數(shù)值振蕩，通過Fortran語言編程求解Euler方程，實(shí)現(xiàn)了凝聚相炸藥爆轟波沖擊空穴塌陷過程的高精度數(shù)值模. 模擬結(jié)果的起爆距離以及沖擊起爆過程中的壓力歷史曲線圖和試驗(yàn)結(jié)果吻合較好，驗(yàn)證了該方法的可靠性. 計(jì)算結(jié)果很好地反映了爆轟波在空穴處的傳播以及壓力、密度、能量等的匯聚過程，給出了在空穴處的高溫高壓熱點(diǎn)形成以及空穴界面運(yùn)動的演化.

[1] 溫麗晶，段卓平，張震宇，等.彈粘塑性雙球殼塌縮熱點(diǎn)反應(yīng)模型[J].高壓物理學(xué)報(bào)，2011，25(6):493-500.

Wen Lijing， Duan Zhuoping， Zhang Zhenyu. An elastic/viscoplastic pore collapse model of double-layered hollow sphere for hot-spot ignition in shocked explosive[J]. Chinese Journal of High Pressure Physics， 2011，25(6):493-500. (in Chinese)

[2] 劉源，王成，趙海濤.多物質(zhì)流體動力學(xué)界面追蹤的Particle Level Set方法[J].北京理工大學(xué)學(xué)報(bào)，2013，33(增刊2):172-174.

Liu Yuan， Wang Cheng， Zhao Haotao. Application of particle level set method in the multi-material interface capturing[J]. Journal of Beijing Institute of Technology, 2013，33(suppl 2):172-174. (in Chinese)

[3] Tran L， Udaykumar H S. Simulation of void collapse in an energetic material， part 2: reactive case[J]. Journal of Propulsion & Power， 2006，22(5):959-974.

[4] Jiang G S， Shu C W. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics, 1996，126:202-228.

[5] 章冠人，陳大年.凝聚炸藥起爆動力學(xué)[M].北京:國防工業(yè)出版社,1991.

Zhang Guanren， Chen Danian. Dynamics of condensed explosives initiation detonation[M]. Beijing: National Defense Industry Press, 1991. (in Chinese)

[6] Wang C， Liu X Q. High resolution numerical simulation of detonation diffraction of condensed explosives[J]. International Journal of Computational Methods， 2015，12(2):68-80.

[7] Davison L， Shahinpoor M. High-pressure shock compression of solids III[M]. New York: Springer, 1998.

[8] Lee E L， Tarver C M. Phenomenological model of shock initiation in heterogeneous explosives[J]. Physics of Fluids， 1980，23(12):2362-2372.

(責(zé)任編輯：劉雨)

High Resolution Numerical Simulation of Detonation Wave Shock-Induced Void Collapse in Condensed Explosives

LIU Xin-qiao， WANG Cheng

(State Key Laboratory of Explosion Science and Technology, Beijing Institute of Technology， Beijing 100081， China)

The detonation wave shock-induced collapse of voids in condensed explosives was studied. Combining RGFM with Level-set method for the material interface, the nonphysical oscillation problem was overcome, which caused by high-density ratio and high-pressure ratio between material interfaces. Taking a third-order TVD Runge-Kutta method for time scattering and the fifth-order high-resolution weighted essentially non-oscillatory (WENO) finite difference scheme for space scattering, a parallel solver was developed to simulate the collapse of voids in condensed explosives. Firstly, the detonation distance with different incident pressure and the pressure history curves with the impact initiation process were presented, and then compared with the experimental results. Further the code was used to simulate the void collapse process, and the density, pressure and particle velocity contours in different moments were given. With the density and energy convergence, the jet phenomenon formed can be clearly seen in the hole, and then further form a hot spot in high temperature and high pressure.

condensed explosives; high resolution; numerical simulation; void collapse; hot spots

2015-10-08

國家自然科學(xué)基金資助項(xiàng)目(11325209)

劉新橋(1981—)，男，博士生，E-mail:3120100070@bit.edu.cn.

王成(1972—)，男，教授，博士生導(dǎo)師，E-mail:wangcheng@bit.edu.cn.

O 38

A

1001-0645(2016)04-0354-05

10.15918/j.tbit1001-0645.2016.04.004