松弛方法在計(jì)算凝聚炸藥爆轟問(wèn)題中的應(yīng)用*

陳秋陽(yáng)，于 明

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

?

松弛方法在計(jì)算凝聚炸藥爆轟問(wèn)題中的應(yīng)用*

陳秋陽(yáng)1，于 明2

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

利用松弛近似，將非線性的凝聚炸藥爆轟控制方程轉(zhuǎn)化為線性的松弛方程組，并采用五階WENO格式和五階線性多步顯隱格式對(duì)線性松弛方程組進(jìn)行空間方向和時(shí)間方向的離散，由此建立具有高精度和高分辨率性質(zhì)的計(jì)算凝聚炸藥爆轟的松弛方法。建立的松弛方法可以避免求解Riemann問(wèn)題及計(jì)算非線性通量的Jacobi矩陣，同時(shí)無(wú)需分裂處理反應(yīng)源項(xiàng)。通過(guò)對(duì)凝聚炸藥的平面一維定常爆轟波結(jié)構(gòu)及球面一維聚心、散心爆轟起爆和傳播過(guò)程的數(shù)值模擬，驗(yàn)證了所建立的松弛方法能夠很好地計(jì)算凝聚炸藥爆轟問(wèn)題。

爆炸力學(xué)；松弛方法；高分辨格式；凝聚炸藥；爆轟波

目前工程應(yīng)用最廣泛的計(jì)算凝聚炸藥爆轟問(wèn)題的數(shù)值方法是Lagrange方法[1]，它存在如下主要問(wèn)題：(1)通常采用交錯(cuò)離散網(wǎng)格，總能量不守恒；(2)爆轟波在人為黏性捕捉過(guò)程中往往被抹平，且人為黏性系數(shù)靠經(jīng)驗(yàn)確定；(3)時(shí)間和空間方向均僅有一階精度，且難以推廣到高階精度；(4)使用細(xì)密離散網(wǎng)格時(shí)，網(wǎng)格容易扭曲從而導(dǎo)致計(jì)算失效。隨著爆轟物理精密化要求，對(duì)爆轟細(xì)致流動(dòng)圖像與規(guī)律需要更深入地了解。高精度、高分辨率的Euler方法，由于能夠較好地克服Lagrange方法的缺點(diǎn)，因此在計(jì)算凝聚炸藥爆轟問(wèn)題時(shí)得到越來(lái)越多的關(guān)注。國(guó)外應(yīng)用的一些典型Euler方法有如下主要特點(diǎn)[2-5]：使用二階Godunov型離散格式；采用簡(jiǎn)單形式狀態(tài)方程和化學(xué)反應(yīng)模型；運(yùn)用反應(yīng)率源項(xiàng)分裂格式等。

凝聚炸藥爆轟的控制方程是帶強(qiáng)剛性化學(xué)反應(yīng)源項(xiàng)和復(fù)雜形式狀態(tài)方程的非線性的雙曲型守恒系統(tǒng)。強(qiáng)剛性化學(xué)反應(yīng)源項(xiàng)和復(fù)雜形式狀態(tài)方程這兩個(gè)特點(diǎn)給運(yùn)用高精度、高分辨率格式進(jìn)行數(shù)值計(jì)算帶來(lái)極大困難。

對(duì)強(qiáng)剛性源項(xiàng)的處理，不管是否采用分裂格式，都可能計(jì)算出錯(cuò)誤的強(qiáng)間斷傳播速度。H.C.Yee等分析指出[6-7]，計(jì)算出錯(cuò)誤的強(qiáng)間斷傳播速度來(lái)源于流動(dòng)控制方程中對(duì)流項(xiàng)離散格式的數(shù)值耗散過(guò)大，從而引起虛假的化學(xué)反應(yīng)。因此，正確捕捉爆轟波傳播速度需采用數(shù)值耗散小的高分辨離散格式。目前，大多數(shù)高分辨離散格式都是基于完全氣體等簡(jiǎn)單形式狀態(tài)方程而采用Riemann求解器進(jìn)行的。凝聚炸藥爆轟過(guò)程中的未反應(yīng)固體炸藥和爆轟氣體產(chǎn)物經(jīng)常使用JWL形式、BKW形式、Davis形式、HOM形式等各種復(fù)雜形式狀態(tài)方程，甚至使用SESAME數(shù)據(jù)庫(kù)等形式，并且化學(xué)反應(yīng)混合區(qū)通常在壓力平衡和溫度平衡的假設(shè)下進(jìn)行溫度迭代運(yùn)算處理[8]，這樣，采用基于Riemann求解器的高分辨格式來(lái)離散凝聚炸藥爆轟控制方程是困難的。并且，為了更好地捕捉爆轟強(qiáng)間斷，除了空間項(xiàng)的離散需要高分辨率，時(shí)間項(xiàng)的離散也需要高分辨率。從目前典型成果看[9]，無(wú)分裂的顯隱格式是較好的處理方式：對(duì)流項(xiàng)離散采用高精度顯式格式，剛性源項(xiàng)離散采用高精度隱式格式。

近年來(lái)正在發(fā)展的松弛方法是數(shù)值求解雙曲型守恒系統(tǒng)的一種有力方法，其主要思想是利用松弛近似，將非線性的雙曲型守恒系統(tǒng)轉(zhuǎn)化為線性的雙曲型松弛方程組，當(dāng)松弛率趨近于零并滿足一定的限制條件后，松弛方程組的解趨近于原守恒系統(tǒng)的解[10-13]。松弛方法的優(yōu)點(diǎn)是可以直接采用一些簡(jiǎn)單、高效的高分辨率數(shù)值格式，勿需求解Riemann問(wèn)題及計(jì)算非線性通量的Jacobi矩陣，因此它特別適合Riemann問(wèn)題復(fù)雜或難以求解的雙曲型守恒律方程。本文中將松弛方法應(yīng)用于一維空間的凝聚炸藥爆轟問(wèn)題的數(shù)值計(jì)算。通過(guò)對(duì)凝聚炸藥PBX9404的平面一維定常爆轟波結(jié)構(gòu)及球面一維聚心、散心爆轟起爆和傳播過(guò)程的數(shù)值模擬，驗(yàn)證所建立的松弛方法能夠很好地計(jì)算凝聚炸藥爆轟問(wèn)題。

1 凝聚炸藥爆轟控制方程

凝聚炸藥爆轟在歐拉坐標(biāo)下的一維化學(xué)反應(yīng)流動(dòng)方程為：

(1)

式中:

其中N為幾何因子(N=0表示平面，N=1表示柱面，N=2表示球面)，R為化學(xué)反應(yīng)率，本文中取三項(xiàng)式Ignition-Growth化學(xué)反應(yīng)率[8]：

對(duì)凝聚炸藥，化學(xué)反應(yīng)率源項(xiàng)中包含有很大的數(shù)(如PBX9404炸藥，I=7.43×1011、G1= 3.1、G2=400.0)，因此源項(xiàng)被認(rèn)為是強(qiáng)剛性的。

本文中凝聚炸藥爆轟的未反應(yīng)固體炸藥和爆轟氣體產(chǎn)物采用JWL形式狀態(tài)方程。在化學(xué)反應(yīng)混合區(qū)采用壓力平衡和溫度平衡的假設(shè)下，混合區(qū)的狀態(tài)可以表示(下標(biāo)s表示固相、g表示氣相)為：

式中:V=ρ0/ρ為相對(duì)體積，e為單位質(zhì)量的內(nèi)能，q為化學(xué)反應(yīng)釋放的單位熱量。

2 松弛方程組的建立

利用松弛近似，將凝聚炸藥爆轟控制方程式(1)轉(zhuǎn)化為如下松弛方程組：

(2)

式中:w是中間變量，A=diag[a1,a2,a3,a4]是元素為正的對(duì)角矩陣，0<ε≤1是松弛率。

如此，將非線性的雙曲型守恒系統(tǒng)(1)轉(zhuǎn)化為線性的雙曲型松弛方程組(2)，這樣可以利用松弛方程組(2)的線性特征場(chǎng)構(gòu)造簡(jiǎn)單、高效的高分辨率數(shù)值格式。文獻(xiàn)[10-11]中分析指出，在滿足-ak≤?f(u)/?u≤ak(k=1,2,3,4)條件下，如果松弛率ε→0，則式(2)的解趨于式(1)的解：w→f(u)。

下面簡(jiǎn)單分析松弛率在數(shù)值離散格式中的角色[12]。

由于w趨近于f(u)，可對(duì)w使用Chapman-Enskog級(jí)數(shù)展開(kāi)，有：

(3)

將該展開(kāi)式代入到式(2)中，保留ε的一階項(xiàng)并消去w，可得：

(4)

從式(4)可以看出，引入松弛率相當(dāng)于引入了數(shù)值耗散。

3 松弛方程組的離散求解

為方便求解，將式(2)進(jìn)行對(duì)角化運(yùn)算，有：

(5)

可以看出，方程組(5)具有常系數(shù)線性性質(zhì)，其特征線為dr/dt=±A，特征線上的Riemann不變量為：w±Au。

考慮空間網(wǎng)格均勻的有限差分離散，式(5)的半離散表達(dá)式為：

(6)

方程組(6)進(jìn)行時(shí)間離散時(shí)，可寫(xiě)成如下常微分方程組形式：

(7)

采用保持總變差有界(TVB)和單調(diào)性質(zhì)的五階精度線性多步顯隱格式[9]來(lái)求解式(7)：

(8)

線性多步顯隱格式的起步值采用三階精度的顯隱格式Runge-Kutta方法確定。

至此，獲得凝聚炸藥爆轟流動(dòng)控制方程式(1)在空間和時(shí)間方向均具有五階離散精度的數(shù)值格式。可以看出來(lái)，爆轟控制方程式在離散過(guò)程中沒(méi)有進(jìn)行Riemann問(wèn)題求解。

4 典型算例分析

以凝聚炸藥PBX9404為例，考察其平面一維定常爆轟波結(jié)構(gòu)及球面一維聚心、散心爆轟起爆和傳播性能。選取炸藥尺度為4.0 cm，以CJ條件起爆。PBX9404炸藥的未反應(yīng)物JWL狀態(tài)方程主要參數(shù)為(g-cm-μs制)：As=69.69、Bs=-1.727、R1s=7.8、R2s=3.9、ws=0.857 8、ρ0=1.842；爆轟氣體產(chǎn)物JWL狀態(tài)方程主要參數(shù)為：Ag=8.524、Bg=0.180 2、R1g=4.55、R2g=1.30、wg=0.38、q=0.102；相應(yīng)的Ignition-Growth化學(xué)反應(yīng)率參數(shù)可見(jiàn)文獻(xiàn)[8]。PBX9404炸藥的Von Neumann尖峰值為：pN=0.563，VN=0.607 5，uN=0.347；相應(yīng)的CJ條件為：pCJ=0.370，VCJ=0.740 3，uCJ=0.229，DCJ=0.880 9。數(shù)值模擬時(shí)松弛率取為ε=10-7。

4.1 PBX9404炸藥平面一維定常爆轟波結(jié)構(gòu)

平面炸藥左端起爆，計(jì)算獲得爆轟達(dá)到定常狀態(tài)時(shí)化學(xué)反應(yīng)區(qū)內(nèi)的壓力p、速度v、相對(duì)比容V、反應(yīng)道λ分布，如圖1所示，所用離散網(wǎng)格為Δr=1/1 000、1/2 000、1/5 000、1/10 000 cm；同時(shí)獲得爆轟CJ速度和壓力Von Neumann尖峰值隨網(wǎng)格大小變化關(guān)系，如圖2所示。從計(jì)算結(jié)果看出，計(jì)算出的爆轟前導(dǎo)沖擊波間斷附近沒(méi)有出現(xiàn)非物理振蕩；隨網(wǎng)格變小時(shí)數(shù)值解逐漸趨近精確解；PBX9404炸藥反應(yīng)區(qū)寬度約為0.01 cm(相應(yīng)反應(yīng)區(qū)內(nèi)網(wǎng)格數(shù)目為10、20、50、100)，當(dāng)采用Δr=1/5 000 cm或更小的離散網(wǎng)格尺度(即反應(yīng)區(qū)內(nèi)的網(wǎng)格數(shù)目至少達(dá)到50個(gè))時(shí)，計(jì)算獲得的數(shù)值解與精確解符合良好。圖1中同時(shí)給出Δr=1/5 000 cm網(wǎng)格尺度下二階MUSCL[13]格式的計(jì)算結(jié)果，可以看出本文的五階精度格式結(jié)果更接近精確解。

數(shù)值模擬給出Δr=1/5 000 cm條件下炸藥爆轟起爆過(guò)程幾個(gè)典型時(shí)刻的壓力和速度變化趨勢(shì)，如圖3所示，圖中每條曲線對(duì)應(yīng)的時(shí)刻為t=0.05、0.1、0.2、0.4、0.8、1.2、1.6、2.0、2.4、2.8、3.2、3.6、4.0 μs?？梢钥闯?，起爆約0.2 μs后爆轟基本達(dá)到定常狀態(tài)。

圖1 PBX9404炸藥化學(xué)反應(yīng)區(qū)內(nèi)物理量分布Fig.1 Distrubitions of physical variables in chemical reaction zone of PBX9404

圖2 PBX9404炸藥CJ爆速和Von Neumann壓力與離散網(wǎng)格尺度的關(guān)系Fig.2 Relations of the CJ velocity and Von Neumann pressure to the mesh sizes in PBX9404

圖3 PBX9404炸藥平面爆轟波傳播過(guò)程中壓力和速度分布變化趨勢(shì)Fig.3 Distributions of pressure and velocity of planar detonation wave in PBX9404

4.2 PBX9404炸藥球面一維聚心爆轟的起爆和傳播

炸藥聚心爆轟波在傳播過(guò)程中將逐漸出現(xiàn)超壓爆轟現(xiàn)象，不正確的計(jì)算方法往往會(huì)獲得反常的雙波結(jié)構(gòu)[1]。球狀炸藥外端起爆，計(jì)算獲得炸藥爆轟波傳播過(guò)程幾個(gè)典型時(shí)刻的壓力和速度變化趨勢(shì)，如圖4所示，圖中每條曲線對(duì)應(yīng)的時(shí)刻為t=0.05、0.1、0.2、0.4、0.8、1.2、1.6、2.0、2.4、2.8、3.2、3.6、4.0 μs，所用離散網(wǎng)格為Δr=1/5 000 cm。可以看出，隨著半徑減小爆轟波壓力和速度增大；爆轟波傳播約0.8 μs后出現(xiàn)超壓爆轟現(xiàn)象。

4.3 PBX9404炸藥球面一維散心爆轟的起爆和傳播

圖4 PBX9404炸藥球面聚心爆轟波傳播過(guò)程中壓力和速度分布變化趨勢(shì)Fig.4 Distributions of pressure and velocity of spherically convergent detonation wave in PBX9404

炸藥散心爆轟波在傳播過(guò)程中前導(dǎo)沖擊波陣面后面的物理量會(huì)急劇下降，差的計(jì)算方法不能正確處理發(fā)散幾何因素的影響會(huì)導(dǎo)致爆轟波熄火[1]。球狀炸藥在中心區(qū)域起爆，計(jì)算獲得炸藥爆轟波傳播過(guò)程幾個(gè)典型時(shí)刻的壓力和速度變化趨勢(shì)，如圖5所示，圖中每條曲線對(duì)應(yīng)的時(shí)刻為t=0.05、0.1、0.2、0.4、0.8、1.2、1.6、2.0、2.4、2.8、3.2、3.6、4.0 μs，所用離散網(wǎng)格為Δr=1/5 000 cm。從圖中可以看出：炸藥爆轟波的壓力和速度均隨波面半徑的增加而增加；爆轟波傳播2.8 μs后壓力值變化不明顯，散心爆轟的最大壓力尖峰值約為0.551、最大速度尖峰值約為0.342，均低于相應(yīng)的平面一維爆轟的壓力和速度尖峰值。

圖5 PBX9404炸藥球面散心爆轟波傳播過(guò)程中壓力和速度分布變化趨勢(shì)Fig.5 Distributions of pressure and velocity of spherically divergent detonation wave in PBX9404

5 結(jié) 論

(1)將松弛方法應(yīng)用于計(jì)算凝聚炸藥的爆轟問(wèn)題，采用五階WENO格式和五階線性多步顯隱格式對(duì)爆轟松弛方程組進(jìn)行空間離散和時(shí)間離散，避免了由于復(fù)雜狀態(tài)方程引起的求解Riemann問(wèn)題及計(jì)算非線性通量的Jacobi矩陣的困難。

(2)通過(guò)對(duì)PBX9404炸藥的平面一維定常爆轟波結(jié)構(gòu)及球面一維聚心、散心爆轟波傳播過(guò)程的數(shù)值模擬，驗(yàn)證了所建立的松弛方法能夠很好地計(jì)算凝聚炸藥爆轟的主要特征，同時(shí)可以看出給出的離散格式具有高精度和高分辨率的性質(zhì)。

(3)下一步的工作是將建立的松弛方法推廣到二維空間下的凝聚炸藥爆轟問(wèn)題。

[1] Mader C L. Numerical modeling of explosives and propellants[M]. 2nd ed. New York: CRC Press, 1998.

[2] Henshaw W D, Schwedeman D W. An adaptive numerical scheme for high-speed reactive flow on overlapping grids[J]. Journal of Computational Physics, 2003,19(2):420-447.

[3] Kapila A K, Schwedeman D W, Bdzil J B. A study of detonation diffraction in the ignition-and-growth model[J]. Combustion Theory and Modeling, 2007(11):781-822.

[4] Banks J W, Schwedeman D W, Kapila A K. A study of detonation propagation and diffraction with compliant confinement[R]. UCRL-JRNL-233735, 2007.

[5] Schwedeman D W, Kapila A K, Henshaw W D. A study of detonation diffraction and failure for a model of compressible reactive flow[R]. UCRL-JRNL-M43735, 2010.

[6] Yee H C, Kotov D V, Sjogreen B. Numerical dissipation and wrong propagation speed of discontinuties for stiff source terms[C]∥Proceedings of ICCFD. Hawaii, 2011.

[7] Yee H C, Kotov D V, Shu C W. Spurious behavior of shock-capturing methods: Problems containing stiff source terms and discontuities[C]∥Proceedings of ICCFD7, 2012.

[8] 孫承緯,衛(wèi)玉章，周之奎.應(yīng)用爆轟物理[M].北京：國(guó)防工業(yè)出版社，2000.

[9] Hundsdorfer W, Ruuth S J. IMEX extensions of linear multistep methods with general monotonicity and boundedness properties[J]. Journal of Computational Physics, 2007(225): 2016-2042.

[10] Liu T P. Hyperbolic conservation laws with relaxation[J]. Communications in Mathematical Physics, 1987(108):153-175.

[11] Jin S, Xin Z. The relaxation schemes for systems of conservation laws in arbitrary space dimensions[J]. Communications of Pure and Applied Mathematics, 1995(48):235-276.

[12] Chalabi A. Convergence of relaxation schemes for hyperbolic conservation laws with stiff source[J]. Mathematics of Computation, 1999(68):955-970.

[13] Tang T. Convergence of MUSCL relaxing scheme to the relaxed scheme for conservation laws with stiff source terms[J]. Journal of Scientific Computing, 2000,15(2):173-195.

[14] Henrick A K, Aslam T D, Powers J M. Mapped weighted essentially non-oscillatory schemes: Achieving optimal order near critical points[J]. Journal of Computational Physics, 2005(207):542-567.

(責(zé)任編輯 曾月蓉)

Application of relaxation method for calculating detonation in condensed explosives

Chen Qiu-yang1, Yu Ming2

(1.GraduateSchool,ChineseAcademyofEngineeringPhysics,Beijing100088,China; 2.KeyLaboratoryforComputationalPhysics,InstituteofAppliedPhysicsandComputationalMathematics,Beijing100094,China)

Based on the theory of relaxation approximation, the nonlinear governing equations of detonation in condensed explosives are transformed into linear relaxation systems, which can easily adopt simple and effective high resolution scheme. A fifth-order WENO reconstruction scheme in spatial discretization and a fifth-order IMEX scheme of linear multistep methods with monotonicity and TVB in time discrtiezation are utilized, where it can leave out solving Riemann problem and computing the Jacobian matrix of the nonlinear flux, and make it unnecessary to split the stiff source term resulting from the chemical reaction. The developed method is applied to simulating the steady structure of one-dimensional planar detonation wave and unsteady initiation and propagation of one-dimensional spherically convergent and divergent detonation wave in condensed explosives PBX9404, with the results demonstrating the excellent performance of the present method.

mechanics of explosion; relaxation method; high resolution scheme; condensed explosives; detonation wave

10.11883/1001-1455(2015)06-0785-07

2014-05-07；

2014-10-07

國(guó)家自然科學(xué)基金項(xiàng)目(11272064); 中國(guó)工程物理研究院科學(xué)技術(shù)發(fā)展基金重點(diǎn)項(xiàng)目(2010B0201030)

陳秋陽(yáng)(1990— )，男，碩士研究生;通訊作者： 于 明,yu_ming@iapcm.ac.cn。

O381 國(guó)標(biāo)學(xué)科代碼： 13035

A