凝聚相炸药爆轰波冲击空穴塌陷过程的高精度数值模拟

刘新桥， 王成

(北京理工大学 机电学院， 爆炸科学与技术国家重点实验室，北京 100081)



凝聚相炸药爆轰波冲击空穴塌陷过程的高精度数值模拟

刘新桥， 王成

(北京理工大学 机电学院， 爆炸科学与技术国家重点实验室，北京 100081)

利用高精度数值格式对凝聚相炸药爆轰波冲击空穴塌陷过程进行了研究. 采用RGFM和Level-set相结合的物质界面处理方法，克服了爆炸流场中高密度比、高压强比的物质界面容易引起数值振荡这一问题. 采用5阶WENO高精度有限差分格式进行空间离散，采用3阶TVD Runger-Kutta法进行时间离散，开发了高精度动态并行程序，并运用该程序对凝聚相炸药爆轰波冲击空穴后的塌陷过程进行了数值模拟，给出了PBX-9404炸药在不同入射压力下的起爆距离及冲击起爆过程中的压力历史曲线图，并分别与实验结果对比后进行了验证. 对凝聚相炸药中的空穴塌陷过程进行了高精度数值计算，给出了不同时刻的密度、压力和粒子速度云图，研究结果表明由于空穴的存在，爆轰波的传播过程中会在空穴处形成射流现象以及密度和能量的汇聚，进一步产生高温高压的热点.

凝聚相炸药；高精度；数值模拟；空穴塌陷；热点

凝聚相炸药中的空穴塌陷不仅会导致热点的形成和随后的局部点火点的生成，在炸药爆轰过程中，由于空穴的存着还会影响到炸药的爆轰波传播规律及爆轰性能. 因此，详细了解空穴塌陷和热点生成的背后物理机制有利于加深对爆轰现象的认识，对于完善爆轰理论及炸药装药的优化设计起着十分重要的作用. 温丽晶等[1]基于Kim的弹粘塑性单球壳塌缩模型，得到了新的冲击起爆三项式细观反应速率模型，并把该模型加入DYNA2D中，模拟了PBX-9501炸药的一维冲击起爆过程. 刘源等[2]验证了Particle Level-set方法界面处理的高精度方法，将Particle Level-set方法与流体动力学控制方程进行了耦合对激波与气泡相互作用问题进行了数值模拟. Udaykumar等[3]对HMX炸药中的空穴塌陷过程进行了数值模拟，研究了加载强度和孔隙大小的影响，结果表明，在微米尺度的空隙中，气相化学反应是一个重要因素，但对空隙崩溃本身的影响作用很小. 空隙完全塌陷时将会有大量反应产物生成，由于空穴崩溃的时间太短，无法在崩溃之前建立能够引发炸药起爆的热点.

本文利用高精度数值格式对凝聚相炸药冲击空穴塌陷过程进行了研究. 采用RGFM和Level-set相结合的物质界面处理方法，在数值上采用5阶WENO格式[4]和3阶TVD Runge-Kutta法对描述炸药起爆和爆轰波传播的Euler方程组进行离散，并开发了高精度动态并行程序. 用所开发的程序对炸药的空穴塌陷过程进行了模拟，给出了由气泡存在而引起的局部射流现象的过程以及气泡存在对爆轰波传播规律的影响.

1 化学反应流体动力学控制方程组

由于爆轰过程的时间极短，因而可以忽略爆轰过程中的黏性、热传导和扩散等输运效应，则非定常可压缩Euler方程组可以作为爆轰过程中流体动力学控制方程组. 其形式为

(1)

(2)

(3)

(4)

(5)

(6)

式中：ρ为密度；p为压强；u，v分别为x，y方向的速度；E为单位质量总能量；e为比内能.

状态方程采用JWL状态方程，反应物和产物的状态方程参数值的表达式如式(7)所示，对于PBX-9404炸药，状态方程中参数见表1[5].

(7)

式中：ν为比容；ν0为初始时刻炸药的比容；T为温度；A，B，R1，R2，Cυ和ω为常数.

表1 PBX-9404炸药JWL状态方程参数

对于凝聚相爆轰的化学反应率模型，使用较多的是Lee和Tarver提出的点火和成长反应率模型， 对于炸药PBX-9404，反应速率模型中的参数见表2[5].

表2 PBX-9404炸药Lee-Tarver模型参数

(8)

式中：I，G，Z均为常数；ν1为受冲击但未反应炸药的比容；η为未反应炸药的相对压缩度.

2 界面处理方法

2.1 Level-set界面追踪

本文采用Level-set物质界面跟踪方法，通过构造界面函数φ(x，y，t)，使得在任意时刻，运动的物质界面Γ(t)正好是φ(x，y，t)的零等值面，即

(9)

φ函数选取的初值应满足在Γ(0)附近为法向单调，而且在Γ(t)上为0. 一般取φ(x，y， 0)为点(x，y)到界面Γ(0)的符号距离，即

(10)

式中d((x，y)，Γ(0))为(x，y)到Γ(0)的距离. 因此，Level-set方法是一种隐式的界面追踪方法. 隐式是指在任意时刻，仅仅给定计算节点到界面的法向符号距离，通过各个节点到界面的距离来隐式地刻画出界面位置所在. 通常情况下，在任意时刻，所有计算节点的φ函数值都不会恰好为0，也即不存在任意一个节点正好处于运动物质界面上.

2.2 RGFM界面处理

RGFM方法首先在物质界面处构建局部Riemann问题，并将局部Riemann问题的求解结果赋予真实流体一侧紧邻物质界面的带状区域，然后通过常值法向外推，在该介质以外的流场赋予该介质的Ghost物理量，将多介质问题转化为若干单介质问题进行处理. 通过求解局部Riemann问题，可以确保物质界面两侧的压强和法向速度始终保持一致，而密度和切向速度等物理量则允许存在间断，这与实际流场中物质界面两侧的间断条件一致. 求解Riemann问题就是求解物质界面附近的流场相互作用之后的流动状态，因此，RGFM边界条件定义方法可以真实的反映介质之间的接触作用过程.

3 数值方法

采用有限差分WENO格式构造方法对控制方程的空间进行离散，需求出通量F(U)、G(U)相应的特征矩阵及其左右特征向量，由于在化学反应区由于同时存在爆轰产物和未反应炸药，压力没有具体的表达形式，需要将问题归结为一个线性方程组的求解，从而得到在化学反应区内包含各守恒量的压力表达形式[6]. 并通过左特征向量将WENO格式中用到的值投影到局部特征空间，实现WENO重构后，将所得到的结果通过右特征向量转化到物理空间. 有限差分格式为

(11)

时间上采用3阶TVD Runge-Kutta方法进行离散：

(12)

4 数值模拟结果与分析

4.1 PBX-9404炸药冲击起爆数值模拟与实验对比

本算例对冲击起爆数值模拟结果与实验结果[7-8]进行了对比，初始参数如表3所示[5]. 图1(a)给出了不同炸药的入射压力与起爆距离之间的关系，曲线6为炸药PBX-9404的实验结果，本文计算结果在图中通过标记各点后连线得到. 图1(b)给出了炸药在2.5 GPa入射压力下压力随时间变化曲线图. 由图中可以看到，本文的数值模拟结果与实验吻合较好.

表3 PBX-9404炸药参数

4.2 PBX-9404炸药空穴塌陷过程数值模拟

取长宽分别为15和10 mm的PBX-9404炸药，空穴中心距离左端5 mm处，在炸药左端为2 mm范围内采用CJ条件起爆. 对物质界面采用Local Level Set方法进行追踪和Real Ghost Fluid Method(RGFM)的界面处理方法. 不同时刻的密度和压力云图如图2所示.

由图2可见，在爆轰波到达空穴的初始阶段，由于空穴对爆轰波的稀疏作用，将使得爆轰波阵面的密度和压力迅速降低并导致波速下降，但在爆轰波与空穴边界接触位置的密度和压力将会产生汇聚，略有增加，此时在空穴内部x方向的粒子速度将急剧增加. 随着爆轰的向前传播，在空穴内部将形成射流，并且空穴被进一步压缩，在空穴压缩塌陷过程中，会出现高温、高压和高密度的现象，导致热点出现，当爆轰波阵面完全扫过空穴时，空穴被压缩为一个极小的高密度区域，此时，压力和粒子速度最大的区域出现在爆轰波阵面绕过空穴后汇聚的位置.

图3～图5给出了3个不同时刻密度、压力及x方向粒子速度曲线. 由密度图可见，随着爆轰波的推进密度开始汇聚，并于0.78 μs在7.5 mm处达到最大值6.3 g/cm3，这已远大于稳定爆轰时的密度峰值. 由压力图可以看到，在0.25 μs时压力峰值约为54.5 GPa，而在0.46 μs时，气泡左端壁面的压力已小于CJ压力，而在气泡的右端由于汇聚出现了45.0 GPa的峰值. 在0.78 μs时出现两个压力峰值，气泡内部的压力进一步汇聚达到57.0 GPa，而此时在气泡右端内壁面处压力已达到78.0 GPa. 由粒子速度曲线图可见，随着时间的推进，粒子速度峰值由0.46 μs到0.78 μs急剧下降到4.5 km/s，同时在7.5 mm位置处粒子速度具有极小值，与压力曲线图对比可知此处的压力具有极大值，这是由于粒子在空穴壁面反射速度反向后与向前运动的粒子进行了碰撞而导致的.

5 结 论

采用5阶WENO高精度有限差分格式以及3阶TVD Runger-Kutta方法分别在空间和时间上进行离散，基于RGFM和Level-set相结合的物质界面处理方法，克服了爆炸流场中高密度比、高压强比物质界面引起的数值振荡，通过Fortran语言编程求解Euler方程，实现了凝聚相炸药爆轰波冲击空穴塌陷过程的高精度数值模. 模拟结果的起爆距离以及冲击起爆过程中的压力历史曲线图和试验结果吻合较好，验证了该方法的可靠性. 计算结果很好地反映了爆轰波在空穴处的传播以及压力、密度、能量等的汇聚过程，给出了在空穴处的高温高压热点形成以及空穴界面运动的演化.

[1] 温丽晶，段卓平，张震宇，等.弹粘塑性双球壳塌缩热点反应模型[J].高压物理学报，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] 刘源，王成，赵海涛.多物质流体动力学界面追踪的Particle Level Set方法[J].北京理工大学学报，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] 章冠人，陈大年.凝聚炸药起爆动力学[M].北京:国防工业出版社,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.

(责任编辑：刘雨)

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

国家自然科学基金资助项目(11325209)

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

王成(1972—)，男，教授，博士生导师，E-mail:wangcheng@bit.edu.cn.

O 38

A

1001-0645(2016)04-0354-05

10.15918/j.tbit1001-0645.2016.04.004