爆轰产物组成和爆轰参数计算方法的理论研究

杜明燃; 汪旭光; 郭子如; 颜事龙

doi:10.11883/1001-1455(2015)04-0449-05

爆轰产物组成和爆轰参数计算方法的理论研究

doi: 10.11883/1001-1455(2015)04-0449-05

1.
安徽理工大学土木建筑学院，安徽淮南 232001
2.
北京矿冶研究总院，北京 100044
3.
安徽理工大学化学工程学院，安徽淮南 232001

基金项目: 国家自然科学基金项目(51134012)

详细信息

作者简介:
杜明燃(1987-), 男, 博士研究生, dumingranaust@163.com

中图分类号: O381
计量
- 文章访问数: 4784
- HTML全文浏览量: 674
- PDF下载量: 860
- 被引次数: 18
出版历程
- 收稿日期: 2013-04-12
- 修回日期: 2014-02-27
- 刊出日期: 2015-08-10

Theoretical studies for calculating the detonation products and properties of explosives

1.
School of Civil Engineering and Architecture, Anhui University of Science and Technology, Huainan 232001, Anhui, China
2.
Beijing General Research Institute of Mining and Metallurgy, Beijing 100044, China
3.
School of Chemical Engineering, Anhui University of Science and Technology, Huainan 232001, Anhui, China

摘要

摘要: 为实现爆轰产物组成和爆轰参数的计算，采用拉格朗日乘数法和牛顿迭代的方法预测爆轰产物组成，利用BKW状态方程预测爆轰参数，在0~600 GPa和300~15 000 K压力温度范围内选取金刚石作为碳的生成相；对爆轰产物系统采用最小自由能原理，结合牛顿迭代法求解爆轰产物的化学平衡方程组；对BKW状态方程参数提出修订，取α=0.5，β=0.298，θ=6 620，κ=9.50；采用自编程序实现计算过程。使用此方法和Hugoniot关系计算密度为1.77 g/cm³的PETN爆轰CJ点爆轰参数验证计算精度，结果显示计算与实验结果的误差均小于1%。利用此方法结合Hugoniot关系预测出爆轰CJ点的产物密度为2.43 g/cm³。
- 爆炸力学 /
- 爆轰参数 /
- 热力学性质 /
- 状态方程 /
- 吉布斯自由能
Abstract: In order to calculate the detonation products and parameters, Lagrange multiplier and Newton iterative method were used to predict detonation products. The state equation of BKW was used to predict detonation parameters. In a range of pressure from 0 to 600 GPa and temperature from 300 to 15 000 K, diamond was intended as the elemental carbon product. Based on the principle of minimum free energy, the equilibrium compositions of detonation products were calculated by using Newton iterative method, which need not calculate the free energy of each composition. The parameters of the state equation of BKW were modified.α=0.5; β=0.298; θ=6 620; κ=9.50. Using self-made program, the detonation properties at CJ point of PETN, whose density is 1.77 g/cm³, were calculated with the theory in this paper and the equation of Hugoniot. The results show satisfactory agreement with the experimental data, with the error less than 1%. The density of detonation products is also predicted easily. When the density of PETN is 1.77 g/cm³, the density of detonation products is 2.43 g/cm³.
- mechanics of explosion /
- detonation parameters /
- thermodynamic property /
- equation of state /
- Gibbs free energy

HTML全文

温压弹作为一种非常规武器, 具有耐高温高压、作用时间长、作用范围广等优点。在对固体云爆剂的研究中, 铝作为一种质轻含能高的金属被广泛采用, 因此对于炸药粉尘铝粉尘等多种粉尘形成的两相爆轰过程发展及传播的研究具有极其重要的意义。含铝复合装药具有爆温爆压高、作用范围广的优点^[1]；陈朗等^[2]研究了铝粉尘尺寸对于爆炸性能的影响, 发现了小尺寸粉尘作用时间短放能快；洪滔等^[3]数值模拟了炸药颗粒与铝颗粒混合粉尘的双波阵面爆轰, 分析对比了双颗粒与单颗粒爆轰的爆轰波参数。

爆轰模拟过程中最重要的是捕捉强间断问题, 近些年发展出了一系列的格式方法, 例如TVD格式、ENO格式、WENO格式、BGK格式和RKDG格式等。CE/SE(时-空守恒元解元)方法是一种解双曲守恒方程的新方法, 它在概念和方法上都不同于已有的方法^[4], 将时间和空间作为一个整体统一对待, 使得从整体和局部都满足守恒律。Y.Wu等^[5]将这种方法拓展应用到化学反应流。董贺飞等^[6]应用CE/SE方法模拟了RDX颗粒在空气中的爆轰问题。

本文中应用CE/SE方法并采用两相流模型, 分析颗粒不同密度对爆轰速度的影响, 并研究RDX颗粒和铝颗粒悬浮粉尘在复杂通道中的爆轰过程。

1. 数值方法及爆轰模型

时-空守恒元解元方法^{[4, 7-8]}与传统差分方法相比, 具有构造简单、格式精度高、分辨率较高等特点。对于二维守恒方程组：

$\frac{\partial U}{\partial t}+\frac{\partial E(U)}{\partial x}+\frac{\partial F(U)}{\partial y}=0$

(1)

CE/SE求解区域模型如图 1所示，图中A、B、C、D和P为记录物理量的网格点, P为在x-y平面的投影, 与A、B、C、D相差Δt。假设解元物理量连续, 采用泰勒展开, 对守恒方程进行积分, 得到CE/SE格式, 计算格式如下:

$U=\frac{1}{4}\left(\overline{U}+\frac{\Delta t}{\Delta x} \overline{E}+\frac{\Delta t}{\Delta y} \overline{F}\right)$

(2)

图 1 CE/SE求解区域模型

Figure 1. CE/SE area model

下载: 全尺寸图片幻灯片

式中：

$\overline{U}=U\left(\frac{\Delta x}{4}, \frac{\Delta y}{4}, 0\right)_{A}+U\left(-\frac{\Delta x}{4}, \frac{\Delta y}{4}, 0\right)_{B}+U\left(-\frac{\Delta x}{4},-\frac{\Delta y}{4}, 0\right)_{C}+ \\ \;\;\;\;\;\;U\left(\frac{\Delta x}{4},-\frac{\Delta y}{4}, 0\right)_{D}$

(3)

$\overline{E}=E\left(0, \frac{\Delta y}{4}, \frac{\Delta t}{4}\right)_{A}-E\left(0, \frac{\Delta y}{4}, \frac{\Delta t}{4}\right)_{B}-E\left(0,-\frac{\Delta y}{4}, \frac{\Delta t}{4}\right)_{C}- \\ \;\;\;\;\;\;E\left(0,-\frac{\Delta y}{4}, \frac{\Delta t}{4}\right)_{D}$

(4)

$\overline{F}=F\left(0, \frac{\Delta y}{4}, \frac{\Delta t}{4}\right)_{A}-F\left(0, \frac{\Delta y}{4}, \frac{\Delta t}{4}\right)_{B}-F\left(0,-\frac{\Delta y}{4}, \frac{\Delta t}{4}\right)_{C}- \\ \;\;\;\;\;\;F\left(0,-\frac{\Delta y}{4}, \frac{\Delta t}{4}\right)_{D}$

(5)

其中H(Δx, Δy, Δt)=H+H_XΔx+H_yΔy(H=U, E, F)为物理量的泰勒展开。

物理量偏导数为:

$U_{I}=\frac{\left|U_{I}^{-}\right|^{\alpha} U_{1}^{+}+\left|U_{I}^{+}\right|^{\alpha} U_{I}^{-}}{\left|U_{I}^{-}\right|^{\alpha}+\left|U_{I}^{+}\right|^{\alpha}}$

(6)

式中：I=x, y; α=1~2;U_I^-, U_I⁺为在相应方向上的左右导数, 对其采用中心差分。然后带入方程得到其他物理量导数。

$U_{x}^{-}=-\frac{U\left(0,0, \frac{\Delta t}{2}\right)_{A}+U\left(0,0, \frac{\Delta t}{2}\right)_{B}-2 U_{p}}{\Delta x}$

(7)

$U_{x}^{+}=\frac{U\left(0,0, \frac{\Delta t}{2}\right)_{C}+U\left(0,0, \frac{\Delta t}{2}\right)_{D}-2 U_{P}}{\Delta x}$

(8)

$U_{y}^{-}=-\frac{U\left(0,0, \frac{\Delta t}{2}\right)_{A}+U\left(0,0, \frac{\Delta t}{2}\right)_{C}-2 U_{p}}{\Delta y}$

(9)

$U_{y}^{+}=-\frac{U\left(0,0, \frac{\Delta t}{2}\right)_{B}+U\left(0,0, \frac{\Delta t}{2}\right)_{D}-2 U_{p}}{\Delta y}$

(10)

采用CE/SE方法求解不含源项欧拉方程, 对欧拉方程中的源项采用4阶龙格库塔方法求解。

模拟过程采用了两相流模型, 假设颗粒为球形, 初始直径都相同, 颗粒的温度都是均衡的, 忽略了粒子间的作用, 粒子化学反应产生的能量假定都被气体吸收, 气体的组分都是均匀的, 略了粒子和墙壁间的热传导, 忽略粒子与气体间的辐射作用, 固相和气相都满足方程。

气相方程：

$\frac{\partial \rho_{1} \varphi_{1}}{\partial t}+\frac{\partial \rho_{1} \varphi_{1} u_{1}}{\partial x}+\frac{\partial \rho_{1} \varphi_{1} v_{1}}{\partial y}=I_{\mathrm{a}}+I_{\mathrm{aa}}$

(11)

$\frac{\partial \rho_{1} \varphi_{1} u_{1}}{\partial t}+\frac{\partial\left(\rho_{1} \varphi_{1} u_{1}^{2}+p\right)}{\partial x}+\frac{\partial \rho_{1} \varphi_{1} u_{1} v_{1}}{\partial y}=I_{\mathrm{a}} u_{\mathrm{a}}-F_{\mathrm{a}x}+I_{\mathrm{aa}} u_{\mathrm{aa}}-F_{\mathrm{aa}x}$

(12)

$\frac{\partial \rho_{1} \varphi_{1} v_{1}}{\partial t}+\frac{\partial \rho_{1} \varphi_{1} u_{1} v_{1}}{\partial x}+\frac{\partial\left(\rho_{1} \varphi_{1} v_{1}^{2}+p\right)}{\partial y}=I_{\mathrm{a}} v_{\mathrm{a}}-F_{\mathrm{a}y}+I_{\mathrm{aa}} v_{\mathrm{aa}}-F_{\mathrm{aa}y}$

(13)

$\begin{array}{*{20}{c}} {\frac{{\partial {\rho _1}{\varphi _1}\left( {{e_1} + 0.5\left( {u_1^2 + v_1^2} \right)} \right)}}{{\partial t}} + \frac{{\partial {\varphi _1}{u_1}\left( {{\rho _1}{e_1} + 0.5{\rho _1}\left( {u_1^2 + v_1^2} \right) + p} \right)}}{{\partial x}} + }\\ {\frac{{\partial {\varphi _1}{v_1}\left( {{\rho _1}{e_1} + 0.5{\rho _1}\left( {u_1^2 + v_1^2} \right) + p} \right)}}{{\partial y}} = }\\ { - {Q_{\rm{a}}} + {I_{\rm{a}}}\left( {{e_{\rm{a}}} + 0.5\left( {u_{\rm{a}}^2 + v_{\rm{a}}^2} \right)} \right) + {I_{\rm{a}}}{q_{{\rm{RDX}}}} - {F_{{\rm{a}}x}}{u_{\rm{a}}} - {F_{{\rm{a}}y}}{v_{\rm{a}}} - {Q_{{\rm{aa}}}} + }\\ {{I_{{\rm{aa}}}}\left( {{e_{{\rm{aa}}}} + 0.5\left( {u_{{\rm{aa}}}^2 + v_{{\rm{aa}}}^2} \right)} \right) + {I_{{\rm{aa}}}}{q_{{\rm{Al}}}} - {F_{{\rm{aa}}x}}{u_{{\rm{aa}}}} - {F_{{\rm{aay}}}}{v_{{\rm{aa}}}}} \end{array}$

(14)

固相方程：

$\frac{\partial \rho_{2} \varphi_{2}}{\partial t}+\frac{\partial \rho_{2} \varphi_{2} u_{2}}{\partial x}+\frac{\partial \rho_{2} \varphi_{2} v_{2}}{\partial y}=-I_{2}$

(15)

$\frac{\partial \rho_{2} \varphi_{2} u_{2}}{\partial t}+\frac{\partial \rho_{2} \varphi_{2} u_{2}^{2}}{\partial x}+\frac{\partial \rho_{2} \varphi_{2} u_{2} v_{2}}{\partial y}=-I_{2} u_{2}+F_{2 x}$

(16)

$\frac{\partial \rho_{2} \varphi_{2} v_{2}}{\partial t}+\frac{\partial \rho_{2} \varphi_{2} u_{2} v_{2}}{\partial x}+\frac{\partial \rho_{2} \varphi_{2} v_{2}^{2}}{\partial y}=-I_{2} v_{2}+F_{2 y}$

(17)

$\begin{array}{c}{\frac{\partial \rho_{2} \varphi_{2}\left(e_{2}+0.5\left(u_{2}^{2}+v_{2}^{2}\right)\right)}{\partial t}+\frac{\partial \rho_{2} \varphi_{2} u_{2}\left(e_{2}+0.5\left(u_{2}^{2}+v_{2}^{2}\right)\right)}{\partial x}+} \\ {\frac{\partial \rho_{2} \varphi_{2} v_{2}\left(e_{2}+0.5\left(u_{2}^{2}+v_{2}^{2}\right)\right)}{\partial y}=Q_{2}-I_{2}\left(e_{2}+\left(u_{2}^{2}+v_{2}^{2}\right)\right)+F_{2 x} u_{2}+F_{2 y} v_{2}}\end{array}$

(18)

$\frac{\partial n_{2}}{\partial t}+\frac{\partial n_{2} u_{2}}{\partial x}+\frac{\partial n_{2} v_{2}}{\partial y}=0$

(19)

组分方程：

$\frac{\partial \rho_{1} \varphi_{1} y_{i}}{\partial t}+\frac{\partial \rho_{1} \varphi_{1} y_{i} u_{1}}{\partial x}+\frac{\partial \rho_{1} \varphi_{1} y_{i} v_{1}}{\partial y}=\omega_{i}$

(20)

气体状态方程：

$p=\rho R T \sum\limits_{1}^{m} \frac{y_{i}}{w_{i}}$

(21)

式中：角标1代表气体, 2=a, aa分别代表RDX和Al, 变量ρ为密度, u为横向速度, v为纵向速度, e为内能, p为压力, φ(φ₁+φ+φ_aa=1)为体积分数, I为单位体积内颗粒的质量变化率。气体中假设存在7种组分, 分别为O₂、N₂、CO、CO₂、H₂O、Al(gas)和Al₂O₃, 质量生成率 $\omega_{1}=-\frac{8}{9} I_{\rm{aa}}$ , $\omega_{2}=\frac{42}{111} I_{\rm{a}}$ , $\omega_{3}=\frac{42}{111} I_{\rm{a}}$ , ω₄=0, $\omega_{5}=\frac{27}{111} I_{\rm{a}}$ , ω₆=0, T₁≤T_Al₂O₃, $\omega_{7}=\frac{17}{9} I_{\rm{aa}}$ , 气体中考虑可逆反应CO+0.5O₂ $\rightleftharpoons$ CO₂。Y为组分浓度, W为分子量。

对于RDX炸药, 采用洪滔等^[9]提出的模型, 计算公式如下：

${I_{\rm{a}}} = \left\{ {\begin{array}{*{20}{l}} 0&{{T_{\rm{a}}} < {T_{\rm{m}}}}\\ {{Q_{\rm{a}}}/L}&{{T_{\rm{a}}} \ge {T_{\rm{m}}}} \end{array}} \right.$

(22)

式中：T_m为炸药颗粒的熔点, L为炸药潜热。炸药颗粒在激波作用下运动, 并在热传导作用下升温, 当到达炸药颗粒熔点后开始熔化, 发生剥离现象, 在高温气体环境中瞬时分解, 释放能量。

对于铝颗粒, 采用以下模型^[10]：

$I_{\mathrm{aa}}=-n_{\mathrm{aa}} \rho_{\mathrm{aa}} 4 \mathsf{ π} R_{\mathrm{aa}}^{2} \frac{\mathrm{d} R}{\mathrm{d} t}$

(23)

$\frac{1}{R} \frac{\mathrm{d} R}{\mathrm{d} t}=-\frac{1}{k d_{0}^{m} / \Psi^{0.9}}$

(24)

式中：n_aa为单位体积内粒子数, R_aa为粒子半径, ${R_{{\rm{aa}}}} = \sqrt[3]{{\frac{{3\varphi }}{{4{n_{{\rm{aa}}}}{\rm{\pi }}}}}}$ , d₀为粒子初始直径, ψ为气体中氧气的摩尔份额, m=1.75, F为气体对粒子的压力,

$F_{\mathrm{a}x}=\frac{n \mathsf{ π} R^{2} C_{\mathrm{a}} \rho_{1} \sqrt{\left(u_{1}-u_{2}\right)^{2}+\left(v_{1}-v_{2}\right)^{2}}\left(u_{1}-u_{2}\right)}{2}$

(25)

${C_{\rm{a}}} = \left\{ \begin{array}{l} \frac{{24 \times \left( {1 + \frac{{\mathit{R}{\mathit{e}^{\frac{2}{3}}}}}{6}} \right)}}{{\mathit{Re}}}\;\;\;\;\;\;\;\;\;\;\mathit{Re} < 1\;000\\ 0.44\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathit{Re} \ge 1\;000 \end{array} \right.$

(26)

式中：Re为雷诺数, Q_a为气体与粒子间的热传导,

${Q_{\rm{a}}} = 4n{\rm{ \mathsf{ π} }}{R^2}{\lambda _1}Nu\left( {{T_1} - {T_2}} \right)/\left( {2R} \right)$

(27)

式中：λ₁为气体导热系数, Nu为Nusselt数, q为单位质量的粒子的反应能。

当温度高于金属氧化物沸点T_Al₂O₃时, 金属氧化物会发生分解^[11], 生成气态Al和O₂, 因此温度将保持在金属氧化物沸点。

2. 模拟结果

为了验证程序的正确性, 分别对单独炸药颗粒和铝颗粒情况进行了模拟。炸药密度为750 g/m³, 半径为20 μm, 起爆条件为φ₁=1、ρ₁=3 kg/m³、u₁=1 000 m/s、T₁=3 600 K。粉尘从左端起爆, 边界条件为左端封闭, 右端开口。得到爆轰波传播速度约为1.899 km/s, 峰值压力约4.75 MPa。S.Eidelman等^[12]得到的峰值压力为4.84 MPa, D.L.Zhang等^[7]在计算中得到的爆轰波传播速度为1 802 m/s, 董贺飞等^[13]得到的爆轰波速度为1 916 m/s, 铝粉尘的密度为304 g/m³, 铝颗粒半径为1.7 μm, 爆轰管直径为15.2 cm。起爆条件为φ₁=1、ρ₁=3 kg/m³、u₁=1 000 m/s、T₁=3 600 K。从左端起爆, 边界条件为左端闭合, 右端开口, 上下为固壁。得到的爆轰波参数为D=1 630 m/s、ρ_CJ=2.43 kg/m³、u_CJ=673 m/s、P_CJ=2.04 MPa、T=3 800 K, 峰值压力为P=3.31 MPa, 文献[14]中结果为D=1 630 m/s、ρ_CJ=2.48 kg/m³、u_CJ=681 m/s、P_CJ=1.91 MPa、T=3 800 K, 与A.J.Tulis等^[15]由实验中得到的铝粉尘爆速1 650 m/s符合较好。

3. 变密度下双粉尘在空气中的爆轰模拟

由于炸药颗粒与铝颗粒混合爆轰问题的复杂性, 我们还研究了颗粒定直径下的爆轰波波速及峰值压力随密度的变化关系。炸药密度565 g/m³、半径为10 μm, 铝粉尘半径为3.5 μm, 起爆条件为φ₁=1、ρ₁=2.2 kg/m³、u₁=2 000 m/s、T₁=3 200 K。粉尘从左端起爆, 边界条件为左端封闭, 右端开口。如图 2~3所示, 通过拟合关系可以看出, 在不改变炸药颗粒密度的情况下，爆轰波波速和峰值压力基本随密度呈线性关系。

图 2 爆轰波速度随双粉尘中铝粉尘密度的变化

Figure 2. Wave speed varying with dust density

下载: 全尺寸图片幻灯片

图 3 爆轰波压力随双粉尘中铝粉尘密度的变化

Figure 3. Detonation pressure varying with dust density

下载: 全尺寸图片幻灯片

4. 复杂通道中的爆轰模型

模拟了在复杂通道中两种粉尘的起爆及爆轰波传播模型, 模型如图 4所示, 通道长3.2 m、宽1.5 m、障碍物长0.4 m、宽0.4 m, 分别位于0.5、1.4和2.3 m处, 网格数为600×240。炸药颗粒密度为450 g/m³、半径为20 μm, 铝颗粒密度为140 g/m³、半径为2 μm。从左端起爆, 起爆条件为φ₁=1、ρ₁=10 kg/m³、u₁=1 000 m/s、T₁=3 000 K, 边界条件为左端封闭、右端为开口, 其余为固壁。选取坐标分别为A(0.25, 0.75)、B(1.15, 0.75)、C(1.15, 0.25)、D(2.05, 0.75)、E(2.05, 0.25)、F(3.0, 0.75)、G(3.0, 0.25)的7个点进行研究。

图 4 模拟区域几何模型

Figure 4. Geometrical model of simulated area

下载: 全尺寸图片幻灯片

图 5为爆轰波在复杂通道中传播的压力图, p₁代表气体流场压力。在0.318 ms时, 爆轰波到达0.7 m处, 与0.5 m处障碍物发生反射, 形成反射冲击波, 在第1个障碍物左固壁处形成高压区域, 局部压力达到7 MPa；0.568 ms时, 爆轰波到达1.15 m处, 爆轰波从狭窄通道进入宽阔区域发生绕射, 在拐角区域形成低压区域, 在通道左端可见两道明显的左行反射波, 反射波经过0.5 m处拐点形成向右绕射波, 并与爆轰波形成叠加, 在爆轰波最右端可见明显高压区域；0.871 ms时爆轰波进入第2段狭窄通道, 与第2障碍物左固壁形成的反射波与初始冲击波叠加形成4个高压奇点, 并在通道左端可见反射形成第2道右传冲击波；1.178 ms时, 爆轰波最右端形成高压区域, 由于爆轰波绕射拐角处出现低压区域, 而第2段宽阔通道内由于各种反射波的相互作用, 流场变得复杂, 形成了2个高压区域；1.648 ms时在第3段宽阔通道内形成了类似0.871 ms时第2段宽阔通道内的高压区域, 而此时第2段宽阔通道内的2个高压区向上运动叠加形成了蘑菇状的高压区域。

图 5 不同时刻流场压力演化图

Figure 5. Detonation pressure varying with time

下载: 全尺寸图片幻灯片

图 6中从左至右依次为点A(0.25, 0.75)、B(1.15, 0.75)、D(2.05, 0.75)、F(3.0, 0.75)处的压力值随时间变化曲线。在爆轰波向右传播过程中逐渐趋于稳定, 可以看出当爆轰波到达D点时, 已经变成稳定爆轰, 此时爆轰波压力3.44 MPa, 爆轰波传播速度为1 656.78 m/s。从C、D2点的压力值曲线图中可以见到明显的2个峰值, 第1个峰值为初始爆轰波形成, 第2个峰值为爆轰波在与壁面作用形成的反射波。

图 6 不同点的压力随时间变化曲线

Figure 6. Pressure varying with time at different points

下载: 全尺寸图片幻灯片

图 7中曲线从左到右依次为A(0.25, 0.75)、B(1.15, 0.75)、D(2.05, 0.75)、F(3.0, 0.75)点处的温度随时间变化曲线。从B、D、F这3点的温度曲线可以看出, 当爆轰波到达时, 温度迅速提升到3 400 K左右, 而后温度经过缓慢下降, 在反射波到达后温度会再次小幅度升高, 流场温度始终保持在3 000 K以上。

图 7 不同点的温度随时间变化曲线

Figure 7. Temperature varying with time at different point

下载: 全尺寸图片幻灯片

对比模拟了铝颗粒单颗粒悬浮粉尘在同种情况下的爆轰过程, 图 8中曲线从左到右依次为A(0.25, 0.75)、B(1.15, 0.75)、D(2.05, 0.75)、F(3.0, 0.75)点处的压力随时间变化曲线。对比图 6可以看出, 到达稳定爆轰的距离也变长, 流场压力峰值比双颗粒情况低约1 MPa, 模拟得到的爆轰波传播速度约为1 468.25 m/s, 远小于双颗粒情况。对比图 6与图 8, 在F点处都为稳定爆轰状态, 而双粉尘爆轰时爆轰波经过F点处, 压力下降非常明显, 说明此时爆轰波宽度较窄, 而单颗粒情况下压力下降比双粉尘时明显变缓, 说明此时爆轰波的反应区明显加大。图 9中曲线从左到右依次为A(0.25, 0.75)、B(1.15, 0.75)、D(2.05, 0.75)、F(3.0, 0.75)点处的温度随时间变化曲线。流场温度变化规律与双颗粒情况基本相同, 但是流场温度明显低于双颗粒情况, 并且流场温度变化范围高于双颗粒情况。

图 8 铝单颗粒时不同点的压力随时间变化曲线

Figure 8. Pressure varying with time at different points of Al

下载: 全尺寸图片幻灯片

图 9 铝单颗粒时不同点的温度随时间变化曲线

Figure 9. Temperature varying with time at different points of Al

下载: 全尺寸图片幻灯片

5. 结论

用CE/SE方法模拟了空气中悬浮炸药-铝粉尘的二维爆轰模型, 通过对比文献验证了模型程序的正确性。研究了双粉尘爆轰过程中颗粒密度对爆轰波速度、爆轰波压力的影响, 得到了关于颗粒密度与波速压力的线性关系。模拟了炸药颗粒密度为450 g/m³、半径为20 μm, 铝颗粒密度为140 g/m³、半径为2 μm的悬浮粉尘在复杂通道中的传播过程及流场演化过程, 对比单种铝颗粒在同种情况下的爆轰波过程发现, 双颗粒爆轰明显提高了流场的峰值压力、温度及爆轰波传播速度, 并且爆轰波后的流场温度要比单颗粒的更稳定。

很好地模拟了双粉尘的两相爆轰过程, 为以后研究多粉尘复杂情况下的爆轰奠定了基础。本次模拟研究了微米尺度下颗粒的爆轰过程, 在纳米尺度下粉尘放能加快, 爆轰过程更复杂, 在以后的模拟中还要对纳米尺度下多种粉尘的爆轰过程进行细致的研究。

表 1 BKW状态方程参数取值

Table 1. The parameter values of BKW equition of state

炸药	α	β	θ	κ
RDX	0.5	0.16	400	10.91
工业炸药	0.5	0.16	2 890	10.91
本文中	0.5	0.298	6 620	9.50

下载: 导出CSV

表 2 每摩尔PETN炸药CJ点爆轰产物组成

Table 2. Compositions of detonation products per mole of PETN explosive at CJ point

方法	n(CO₂)/mol	n(CO)/mol	n(N₂)/mol	n(O₂)/mol	n(H₂)/mol
BKW^{[1, 2]}	3.950 9	0.096 2	1.999 2	1.359 1×10^-5	0.568 9×10^-4
LJD^{[1, 2]}	3.779 0	0.446 9	1.969 5	0	0.480 4×10^-1
WCA^[13]	3.934 9	0.154 6	1.991 7	4.494 3×10^-5	3.989 5×10^-4
本文	3.999 5	0.001 1	1.999 9	1.190 4×10^-5	2.077 7×10^-4
方法	n(H₂O)/mol	n(NH₃)/mol	n(NO)/mol	n(CH₄)/mol	n(C)/mol
BKW^{[1, 2]}	3.998 6	1.327 5×10^-4	2.244 1×10^-4	0.537 3×10^-5	0.951 6
LJD^{[1, 2]}	3.929 1	0	0.619 5×10^-1	1.197 9×10^-2	0.759 8
WCA^[13]	3.975 2	1.629 0×10^-2	3.141 7×10^-4	0	0.910 5
本文	3.999 7	7.354 8×10^-8	4.530 7×10^-6	5.482 0×10^-6	0.999 1

下载: 导出CSV

表 3 PETN(1.77 g/cm³)炸药爆轰CJ点的爆轰参数

Table 3. Detonation parameter of CJ for PETN explosive(1.77 g/cm³)

方法	p/GPa	D/(m·s^-1)	T/K	方法	p/GPa	D/(m·s^-1)	T/K
WCA^[2]	31.5	8 416	4 349	LJD^[6]	28.5	8 087	4 378
BKW^[5]	31.8	8 421	2 833	WCA^[7]	33.49	8 540	3 237
实验^[2]	33.5	8 300	3 400	本文	33.53	8 342	3 473

下载: 导出CSV

参考文献(13)

[1]	Chirat R, Pittion-Rossillon G. A new equation of state for detonation products[J]. Journal of Chemical Physics, 1981, 74(8): 4634-4645. doi: 10.1063/1.441653
[2]	Mader C L. Numerical modeling of explosives and propellants[M]. 3rd Edtion. Boca Raton, FL, USA: CRC Press, 2008: 31-63.
[3]	Thiel M V, Ree F H. Nonequilibrium effects of slow diffusion controlled reactions on the properties of explosives[C]∥Proceedings of 9th International Symposium on Detonation. Portland, USA: Oregon, 1991: 743-750.
[4]	Wu X. BKW equation of state for detonation products[C]∥Proceedings of 8th International Symposium on Detonation. Portland, USA: Oregon, 1985: 438.
[5]	李德华, 程新路, 杨向东, 等. PETN、RDX和HMX炸药爆轰参数的数值模拟[J].爆炸与冲击, 2005, 25(4): 325-328. Li De-hua, Cheng Xin-lu, Yang Xiang-dong, et al. Numerical simulation of detonation parameters for PETN, RDX and HMX explosives[J]. Explosion and Shock Waves, 2005, 25(4): 325-328.
[6]	赵艳红, 刘海风, 张弓木. PETN炸药爆轰产物状态方程的理论研究[J].高压物理学报, 2009, 23(2): 143-149. Zhao Yan-hong, Liu Hai-feng, Zhang Gong-mu. Equation of state of detonation products for PETN explosive[J]. Journal of High Pressure Physics, 2009, 23(2): 143-149.
[7]	Fried L E, Howard W M. Explicit Gibbs free energy equation of state applied to the carbon phase diagram[J]. Physical Review, 2000, B61(13): 8734-8743.
[8]	赵艳红, 刘海风, 张弓木.基于统计物理的爆轰产物物态方程研究[J].物理学报, 2007, 56(8): 4791-4797. Zhao Yan-hong, Liu Hai-feng, Zhang Gong-mu. Equation of state of detonation products based on statistical mechanical theory[J]. Acta Physica Sinica, 2007, 56(8): 4791-4797.
[9]	杨祖一.国外民爆器材法规和技术资料汇编[M].北京: 中国爆破器材行业协会, 2010: 290-298.
[10]	杨向东, 谢文, 武保剑.液氮的冲击压缩理论计算[J].高压物理学报, 1998, 12(1): 1-5. Yang Xiang-dong, Xie Wen, Wu Bao-jian. Theoretical calculation for the Hugoniot curves of liquid nitrogen[J]. Journal of High Pressure Physics, 1998, 12(1): 1-5.
[11]	刘福生, 陈先猛, 陈攀森, 等.液态CO₂高温高密度状态方程研究[J].高压物理学报, 1998, 12(1): 28-33. Liu Fu-sheng, Chen Xian-meng, Chen Pan-sen, et al. Equation of state of liquid CO₂ at high temperatures and high densities[J]. Journal of High Pressure Physics, 1998, 12(1): 28-33.
[12]	Ihsan B, Gregor P. Thermochemical Data of Pure Substances: Volumn 2[M]. 3rd Edtion. New York: VCH, 1995.
[13]	李德华, 程新路, 杨向东, 等. PETN炸药爆轰参数的数值模拟[J].四川师范大学学报, 2005, 28(4): 448-451. Li De-hua, Cheng Xin-lu, Yang Xiang-dong, et al. Numerical simulation of detonation parameters for PETN explosive[J]. Journal of Sichuan Normal University, 2005, 28(4): 448-451.

施引文献

期刊类型引用(8)

1.	聂根良，程晓红，周文海，胡才智，凌晓. 冲击载荷作用下炮孔近区节理岩体的动力响应. 矿业研究与开发. 2023(11): 72-80 . 百度学术
2.	刘小鸣，陈士海. 群孔微差爆破的地表振动波形预测及其效应分析. 岩土工程学报. 2020(03): 551-560 . 百度学术
3.	陈立军. 开采爆破作用下台阶岩体损伤区分布规律研究. 爆破. 2020(03): 85-89+114 . 百度学术
4.	程平，王林峰，郑志伟，曾韬睿，吴发友. 隐伏岩溶区小净距隧道爆破振动规律. 科学技术与工程. 2020(24): 10017-10024 . 百度学术
5.	梁瑞，周文海，余建平，李珍宝，杜超飞，王敦繁. 冲击载荷作用下岩体拉-压损伤破坏的边坡抛掷爆破模拟. 高压物理学报. 2019(01): 82-91 . 百度学术
6.	管晓明，聂庆科，李华伟，安建永. 隧道爆破振动下既有建筑结构动力响应及损伤研究综述. 土木工程学报. 2019(S1): 151-158 . 百度学术
7.	仪海豹，张西良，杨海涛，李明，高琪琪，金科. 基于重球触地实验的空区塌落振动分析及治理. 爆炸与冲击. 2019(07): 91-103 . 本站查看
8.	周文海，梁瑞，陈金林，朱冕，陈鹏辉，楼晓明，王敦繁. 时程稳定性系数确定的边坡逐孔起爆孔间微差降振时间. 爆炸与冲击. 2019(08): 183-191 . 本站查看