松弛方法在計算凝聚炸藥爆轟問題中的應用*

陳秋陽，于 明

(1.中國工程物理研究院研究生部，北京 100088； 2.北京應用物理與計算數學研究所計算物理國家重點實驗室，北京 100094)

?

松弛方法在計算凝聚炸藥爆轟問題中的應用*

陳秋陽1，于 明2

(1.中國工程物理研究院研究生部，北京 100088； 2.北京應用物理與計算數學研究所計算物理國家重點實驗室，北京 100094)

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

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

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

凝聚炸藥爆轟的控制方程是帶強剛性化學反應源項和復雜形式狀態方程的非線性的雙曲型守恒系統。強剛性化學反應源項和復雜形式狀態方程這兩個特點給運用高精度、高分辨率格式進行數值計算帶來極大困難。

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

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

1 凝聚炸藥爆轟控制方程

凝聚炸藥爆轟在歐拉坐標下的一維化學反應流動方程為：

(1)

式中:

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

對凝聚炸藥，化學反應率源項中包含有很大的數(如PBX9404炸藥，I=7.43×1011、G1= 3.1、G2=400.0)，因此源項被認為是強剛性的。

本文中凝聚炸藥爆轟的未反應固體炸藥和爆轟氣體產物采用JWL形式狀態方程。在化學反應混合區采用壓力平衡和溫度平衡的假設下，混合區的狀態可以表示(下標s表示固相、g表示氣相)為：

式中:V=ρ0/ρ為相對體積，e為單位質量的內能，q為化學反應釋放的單位熱量。

2 松弛方程組的建立

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

(2)

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

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

下面簡單分析松弛率在數值離散格式中的角色[12]。

由于w趨近于f(u)，可對w使用Chapman-Enskog級數展開，有：

(3)

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

(4)

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

3 松弛方程組的離散求解

為方便求解，將式(2)進行對角化運算，有：

(5)

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

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

(6)

方程組(6)進行時間離散時，可寫成如下常微分方程組形式：

(7)

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

(8)

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

至此，獲得凝聚炸藥爆轟流動控制方程式(1)在空間和時間方向均具有五階離散精度的數值格式。可以看出來，爆轟控制方程式在離散過程中沒有進行Riemann問題求解。

4 典型算例分析

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

4.1 PBX9404炸藥平面一維定常爆轟波結構

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

數值模擬給出Δr=1/5 000 cm條件下炸藥爆轟起爆過程幾個典型時刻的壓力和速度變化趨勢，如圖3所示，圖中每條曲線對應的時刻為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后爆轟基本達到定常狀態。

圖1 PBX9404炸藥化學反應區內物理量分布Fig.1 Distrubitions of physical variables in chemical reaction zone of PBX9404

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

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

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

炸藥聚心爆轟波在傳播過程中將逐漸出現超壓爆轟現象，不正確的計算方法往往會獲得反常的雙波結構[1]。球狀炸藥外端起爆，計算獲得炸藥爆轟波傳播過程幾個典型時刻的壓力和速度變化趨勢，如圖4所示，圖中每條曲線對應的時刻為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，所用離散網格為Δr=1/5 000 cm。可以看出，隨著半徑減小爆轟波壓力和速度增大；爆轟波傳播約0.8 μs后出現超壓爆轟現象。

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

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

炸藥散心爆轟波在傳播過程中前導沖擊波陣面后面的物理量會急劇下降，差的計算方法不能正確處理發散幾何因素的影響會導致爆轟波熄火[1]。球狀炸藥在中心區域起爆，計算獲得炸藥爆轟波傳播過程幾個典型時刻的壓力和速度變化趨勢，如圖5所示，圖中每條曲線對應的時刻為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，所用離散網格為Δr=1/5 000 cm。從圖中可以看出：炸藥爆轟波的壓力和速度均隨波面半徑的增加而增加；爆轟波傳播2.8 μs后壓力值變化不明顯，散心爆轟的最大壓力尖峰值約為0.551、最大速度尖峰值約為0.342，均低于相應的平面一維爆轟的壓力和速度尖峰值。

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

5 結 論

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

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

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

[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] 孫承緯,衛玉章，周之奎.應用爆轟物理[M].北京：國防工業出版社，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.

(責任編輯 曾月蓉)

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

國家自然科學基金項目(11272064); 中國工程物理研究院科學技術發展基金重點項目(2010B0201030)

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

O381 國標學科代碼： 13035

A