一种新爆轰产物状态方程及其在炸药爆轰性能预测上的应用

彭钺; 张蕾; 谢明伟; 麻宏亮; 袁晓霞; PENG Yue; ZHANG Lei; XIE Ming-wei; MA Hong-liang; YUAN Xiao-xia

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

一种新爆轰产物状态方程及其在炸药爆轰性能预测上的应用 PDF

- ORCID：
彭钺
✉
- ORCID：
张蕾
✉
- ORCID：
谢明伟
- ORCID：
麻宏亮
- ORCID：
袁晓霞

陕西应用物理化学研究所瞬态化学效应与控制全国重点实验室，陕西西安 710061

中图分类号： TJ55； O64

最近更新：2024-10-11

DOI：10.11943/CJEM2024021

摘要

为准确描述爆轰产物在高温高压下的热力学关系，拓展应用范围，实现CHNO单质/混合炸药与含能金属盐爆轰性能的可靠预测，基于Exp‑6势的维里系数理论值建立了新气态爆轰产物状态方程Virial‑Peng‑Long（VPL），并基于更准确的“冷压”项与考虑金属高压热运动变化的“晶格振动”项建立了新凝聚态金属产物三项式状态方程Wu‑Chen‑Peng（WCP）。应用VPL计算了几种典型CHNO单质/混合炸药的爆轰参数，应用VPL与WCP评价了几种典型含能金属盐的爆轰CJ参数与驱动做功能力。结果表明，相比VLW和BKW，基于VPL能够更准确评价CHNO单质/混合炸药的爆轰性能，对太安（PETN）爆速预测偏差在±2.1%以内，最高不超过2.5%；对RHT‑901驱动圆筒稳定速度预测误差绝对值<1%。而VPL与WCP能够准确预测含能金属盐爆轰性能，其中对叠氮化铅爆轰CJ参数计算相对误差不超过±4%；对铜叠氮化物驱动飞片速度预测相对误差绝对值小于1%。

图文摘要

关键词

爆炸化学; 爆轰产物状态方程; 含能金属盐; 爆轰性能; 维里系数

0 引言

基于爆轰产物状态方程进行炸药爆轰参数的热力学计算是研究炸药爆轰过程的基础，也是预测未知炸药爆轰性能、指导炸药装药配方设计与新炸药合成的重要方法。近年来，高能量密度含能材料分子设计得到大量研究，爆轰性能热力学计算有助于筛选具有合适爆轰性能的炸药分子^［

1］。因此，如何构建更加准确的爆轰产物状态方程成为国内外学者致力研究的课题，不同种类爆轰产物状态方程被提出，如Becker‑Kistiakosky‑Wilson（BKW）状态方程^{［参考文献 2

百度学术}2］、Jacobs‑Cowperthwaite‑Zwisler（JCZ）状态方程^{［参考文献 3

百度学术}3］、Virial‑Wu（VLW）状态方程^{［参考文献 1

百度学术}1］等。其中，BKW和VLW状态方程基于真实气体的维里状态方程而建立，具有完美的理论背景和简洁的形式，在实际研究中得到广泛应用^{［参考文献 1-2}1-2］。

目前常用的BKW与VLW状态方程仍存在一些不足。其中BKW状态方程提出的较早，也应用最为广泛。但一方面，BKW使用多套参数计算不同类型炸药，缺乏通用性^［

4］；另一方面，BKW状态方程对炸药做功能力的预测明显偏高，韩勇^{［参考文献 4

百度学术}4］认为这是其3阶维里项过于简化所致。与BKW相比，VLW将最高维里阶数拓展到5阶并使用更准确的3阶维里项形式，对炸药爆轰性能预测的准确性显著提高^{［参考文献 4

百度学术}4］，但仍存在不足：一方面，VLW3阶及以上阶数的维里项通过与2阶维里项间的“相似”关系建立，在描述3阶维里系数时与理论值符合较好，但在描述4阶和5阶维里系数时与理论值存在较大偏差^{［参考文献 4

百度学术}4］，使VLW不能有效描述气态产物高温高压热力学关系^{［参考文献 5

百度学术}5］；另一方面，赵雅琦^{［参考文献 6

百度学术}6］指出VLW状态方程使用的L‑J 6‑12分子势并不适于描述爆轰气体产物在高密度下的相互作用。以上问题势必降低VLW对炸药爆轰性能预测的准确性。此外，目前对炸药爆轰性能的热力学计算主要针对CHNO炸药，对含能金属盐的计算研究相对较少，缺乏能准确描述含能金属盐凝聚态金属产物的状态方程。爆轰热力学计算软件EXPLO5对金属产物使用不含温度项的Murnaghan状态方程^{［参考文献 2

百度学术}2］，无法完整描述其pVT关系。Fortran VLW使用Cowan‑Fickett状态方程描述含能金属盐中金属产物^{［参考文献 7

百度学术}7］，但其“冷压”项参数未参考金属零温压缩p‑V关系标定，“热压”项参数标定未考虑金属热振动频率在高压下随体积的变化，存在形式与参数上的不足。

针对以上问题，为提高爆轰热力学计算的准确性，同时拓展计算的材料范围实现对含能金属盐爆轰性能的准确预测。本研究通过理论分析与计算，建立新的气态产物状态方程Virial‑Peng‑Long（VPL）与凝聚态金属产物状态方程Wu‑Chen‑Peng（WCP）。将VPL与WCP应用在CHNO炸药与含能金属盐爆轰参数计算中，通过与实验值及其他状态方程计算值对比分析，并评价新状态方程在预测这两类含能材料的爆轰CJ参数及做功能力方面的准确性与适用性。

1 新爆轰产物状态方程

1.1 气态产物维里状态方程Virial‑Peng‑Long

凝聚炸药的气态爆轰产物在CJ状态下可达到几至几十吉帕的压力与几千开尔文的温度，密度超过炸药自身装药密度，其热力学关系已经不能用理想气体状态方程来描述。维里状态方程考虑了分子体积与分子间相互作用，理论上适合描述气体在高温高压下的热力学关系^［

8］。

\frac{p V_{m}}{R T} = 1 + \frac{B (T)}{V_{m}} + \frac{C (T)}{V_{m}^{2}} + \frac{D (T)}{V_{m}^{3}} + \frac{E (T)}{V_{m}^{4}} + \dots

（1）

式中，V_m是气体的摩尔体积，cm³·mol^-1；B（T）、C（T）、D（T）、E（T）…分别是气体的2阶、3阶、4阶、5阶等维里系数，描述不同温度下气体中2分子、3分子、4分子、5分子等相互作用。低压下维里状态方程可简化为理想气体状态方程，随着压力升高高阶维里项贡献连续增加。这一特性使维里状态方程适合描述炸药爆轰产物从爆轰波CJ点到波后爆轰产物膨胀做功过程发生的热力学关系变化。

维里状态方程中的n阶维里系数可根据Mayer展开理论^［

8］计算：

B_{n} = \frac{1 - n}{n!} N_{A} \int \dots \int \sum_{G} [\sum_{i j \in G} f_{i j}] d ω_{1} \dots d ω_{n - 1}

（2）

f_{i j} = e^{- u (r_{i j}) / k_{B} T} - 1

（3）

式中，k_B是Boltzmann常数，1.3806488×10^-23 J·K^-1；N_A是Avogadro常数，6.022×10²³ mol^-1；u（r_ij）是描述两分子间相互作用的势函数。如前所述，VLW状态方程中使用的L‑J 6‑12势函数并不适于描述气态爆轰产物。对分子间相互作用的理论研究表明，分子间排斥作用使用指数函数描述更为合适^［

6］，因此本研究使用Exp‑6势函数描述气态产物分子：

u (r) = \{\begin{matrix} \infty & r < r_{m} \\ ε [\frac{6}{α - 6} e^{α (1 - r / r^{*})} - \frac{α}{α - 6} {(\frac{r^{*}}{r})}^{6}] & r \geq r_{m} \end{matrix}

（4）

式中，ε是分子间最大吸引能，J；r^*是达到最大吸引能时的分子间距，cm；α是常数，对于气态产物通常取α=13.0^［

3］。

对于Exp‑6势宽广温度范围内的2阶维里系数已有计算研究^［

9］，通过最小二乘法拟合文献值建立2阶维里项。Kihara^{［参考文献 10

百度学术}10］通过对三分子相互作用空间关系进行分析提出了一种简便的3阶维里系数计算公式：

C = - 4 π^{2} N_{A} \int_{0}^{\infty} \int_{0}^{0.75} \int_{1 - \sqrt[]{1 - y}}^{\sqrt[]{1 - y}} R^{5} f_{12} f_{23} f_{13} d R d y d x

（5）

式中，f_ij的表达式如式（3）所示，其下标i、j代表分子i和分子j；R是分子1和分子2之间相对距离；x是分子1与分子3之间相对距离；y是分子2与分子3之间相对距离。使用复合Simpson公式求解式（5）中的三重积分，编写了3阶维里系数计算程序Virial3th，使用该程序计算了k_BT/ε=1~300范围内Exp‑6势的3阶维里系数，并使用最小二乘法对参数进行了拟合建立了3阶维里项。

对于4阶和5阶维里项，四分子和五分子间相互作用的Mayer展开形式十分复杂，很难直接求解。

Barker等^［

11］提出利用Legendre正交积分求解4阶和5阶维里系数的方法，并使用该方法计算了L‑J 6‑12势在k_BT/ε=0.625~20下的4~5阶维里系数。本研究在其他基础上编写了4阶与5阶维里系数计算程序Virial4th与Virial5th。为检验准确性，首先计算了L‑J 6‑12势的4~5阶维里系数并与文献［11］记载的理论值进行了对比，如图1所示。

图1 L‑J 6‑12势4~5阶维里系数计算值与文献值对比

Fig.1 Comparison between the calculated values and the literature values of the 4-5th order Virial coefficients of L‑J 6‑12 potential

从图1中可以看到，使用Virial4th和Virial5th程序可以得到准确的4~5阶维里系数值。对Exp‑6势在k_BT/ε=1~300范围内的4~5阶维里系数进行了计算，对计算结果使用最小二乘法拟合建立了VPL状态方程的4~5阶维里项。

通过基于维里系数理论计算值构建的Exp‑6势下2~5阶维里项，建立了Virial‑Peng‑Long（VPL）状态方程，如下所示：

\frac{p V_{m}}{R T} = 1 + B^{*} (T^{*}) (b_{0} / V_{m}) + C^{*} (T^{*}) {(b_{0} / V_{m})}^{2} + D^{*} (T^{*}) {(b_{0} / V_{m})}^{3} + E^{*} (T^{*}) {(b_{0} / V_{m})}^{4}

（6）

B^{*} (T^{*}) = \{\begin{array}{l} \sum_{i = 1}^{5} b_{2 i - 1} T^{* b_{2 i}} & 1 \leq T^{*} \leq 50 \\ b_{11} + \sum_{i = 6}^{8} b_{2 i} e^{T^{*} / b_{2 i + 1}} & T^{*} > 50 \end{array}

（7）

C^{*} (T^{*}) = \{\begin{array}{l} \sum_{i = 1}^{5} c_{2 i - 1} T^{* c_{2 i}} & 1 \leq T^{*} \leq 7.5 \\ c_{11} + \sum_{i = 6}^{8} c_{2 i} e^{T^{*} / c_{2 i + 1}} & T^{*} > 7.5 \end{array}

（8）

D^{*} (T^{*}) = \{\begin{array}{l} \sum_{i = 1}^{6} d_{2 i - 1} T^{* d_{2 i}} & 1 \leq T^{*} \leq 10 \\ \sum_{i = 7}^{10} d_{2 i - 1} T^{* d_{2 i}} & T^{*} > 10 \end{array}

（9）

E^{*} (T^{*}) = \{\begin{array}{l} \sum_{i = 1}^{4} e_{2 i - 1} T^{* e_{2 i}} & 1 \leq T^{*} \leq 5 \\ e_{9} + \sum_{i = 5}^{7} e_{2 i} e^{T^{*} / e_{2 i + 1}} & T^{*} > 5 \end{array}

（10）

T^{*} = k_{B} T / ε

（11）

b_{0} = \frac{2}{3} π N_{A} r^{* 3}

（12）

式中，b_x、c_x、d_x、e_x分别代表2阶、3阶、4阶、5阶无量纲维里系数‑温度关系拟合式系数。VPL状态方程中的2~5阶维里项参数如表1所示。

表1 VPL状态方程的2~5阶维里项参数

Table 1 Parameters of 2-5th Virial term in VPL equation of state

type	value	type	value	type	value	type	value	type	value
b₁	-1.76358	b₁₅	-57.48673	c₁₂	0.04845	d₉	-0.71993	e₃	-17.56250
b₂	-2.22513	b₁₆	0.14818	c₁₃	-12.27360	d₁₀	-8.50898	e₄	-4.58881
b₃	-2.66841	b₁₇	-1305.43710	c₁₄	0.05856	d₁₁	-0.01212	e₅	18.27770
b₄	-0.75632	c₁	-1.14019	c₁₅	-43.70161	d₁₂	-15.88648	e₆	-8.71393
b₅	3.37443	c₂	-3.12627	c₁₆	0.04603	d₁₃	-7.38927	e₇	0.05914
b₆	-0.38402	c₃	-2.43735	c₁₇	-215.78612	d₁₄	-1.40048	e₈	-0.88679
b₇	1.08068	c₄	-0.83587	d₁	-2.35225	d₁₅	-1.31933	e₉	1.328E‑7
b₈	-2.22385	c₅	5.18488	d₂	-3.78631	d₁₆	-0.87972	e₁₀	0.00612
b₉	-0.81853	c₆	-0.66549	d₃	-2.49397	d₁₇	3.52562	e₁₁	-17.37990
b₁₀	-1.56972	c₇	1.80678	d₄	-1.23009	d₁₈	-0.99915	e₁₂	0.01819
b₁₁	0.05199	c₈	-2.44981	d₅	2.23154	d₁₉	5.25646	e₁₃	-5.24886
b₁₂	0.10393	c₉	-0.21886	d₆	-1.09693	d₂₀	-1.55510	e₁₄	0.00104
b₁₃	-191.46520	c₁₀	-6.32838	d₇	2.20097	e₁	-1.57831	e₁₅	-66.58590
b₁₄	0.06510	c₁₁	0.00992	d₈	-3.78632	e₂	-12.25167

Note： b_x， c_x， d_x， e_x are fitting coefficients of the dimensionless virial coefficient‑temperature relationship for second， third， fourth， and fifth orders， respectively.

在得到准确的维里项形式后，需要确定合适的产物分子势参数，使VPL状态方程能够准确描述气态产物热力学关系。考虑炸药几种主要气态产物，参考其在高温高压下pVT关系实验或分子动力学计算数据^［

12-17］使用MATLAB中的多目标优化遗传算法模块对气态产物分子势参数进行了优化。热力学数据参考文献、优化后的势参数及其与参考值之间的相对误差百分比（MAPE）如表2所示。

表2 VPL状态方程产物分子势参数优化结果

Table 2 Calibration results of molecular potential parameters and correction term parameters of VPL equation of state

species	ε/k_B / K	b₀ / cm³·mol^-1	MAPE / %	ref.
H₂	25.4	57.11	1.248	［12］
O₂	97.6	75.38	1.035	［12］
CO₂	164.4	105.33	1.157	［12-13］
CO	96.4	98.12	0.965	［12， 14］
NO	134.8	76.96	1.248	［15］
N₂	83.9	96.28	0.671	［16］
CH₄	145.3	84.25	1.184	［12-13］

Note： ε/k_B and b₀ are parameters of Exp‑6 potential.

可以看到，优化后的产物分子势参数使VPL能准确描述产物单组分气体的热力学关系，计算相对误差在1%左右。对于强极性产物水，考虑极性分子间相互作用与温度相关性，参考文献［

17］中记载的不同温度下高温高压pVT数据标定了水分子的势参数温度函数，如下所示：

ε / k_{B} = \{\begin{array}{l} \begin{array}{l} 528.25 - 1.28777 T \\ + 0.00171 T^{2} - 5.66667 \times 10^{- 7} T^{3} \end{array} & 773.15 K \leq T \leq 1373.15 K \\ \begin{array}{l} 950.978 - 0.43498 T \\ + 9.34123 \times 10^{- 5} T^{2} - 6.64705 \times 10^{- 9} T^{3} \end{array} & 1373.15 K < T \leq 4273.15 K \end{array}

（13）

b_{0} = \{\begin{array}{l} 46.91216 + \frac{8.47704}{1 + e^{(T - 1501.43796) / 122.12326}} & 773.15 K \leq T \leq 1873.15 K \\ 49.1 + \frac{- 2.60211}{1 + e^{(T - 3122.91058) / 351.45539}} & 1873.15 K < T \leq 4273.15 K \end{array}

（14）

1.2 凝聚态金属产物三项式状态方程Wu‑Chen‑Peng

理论上，对于金属，无论是固体还是液体，其高温高压热力学关系均可通过三项式状态方程描述^［

18］：

p = p_{C} (V) + p_{I o n} (V, T) + p_{E l e} (V, T)

（15）

式中，p_C（V）是描述金属在零温下p‑V关系的“冷压”项；p_Ion（V，T）是考虑金属微观晶格振动的“晶格振动”项；p_Ele（V，T）是考虑微观下金属自由电子运动贡献的“电子运动”项。

对于“冷压”项，本研究使用物理背景清晰，形式简单，普适性强，适合描述金属在较大压力范围内热力学关系的Vinet状态方程^［

18］来计算：

p_{C} (V) = \frac{3 B_{0 K} (1 - x)}{x^{2}} e x p [1.5 (B_{0 K}^{'} - 1) (1 - x)]

（16）

x = {(V / V_{0 K})}^{1 / 3}

（17）

V_{0 K}是金属在零温下的比容，可以根据金属热膨胀关系式近似给出^［

18］。

V_{0 K} \approx \frac{V_{0}}{{(1 + α_{L} Δ T_{0})}^{3}}

（18）

对于“冷压”项参数B_{0 K}和 $B_{0 K}^{'}$ ，根据金属在常温常压下的热力学Gruneisen系数Γ₀与室温（≈300 K）下金属冲击波速度u_s与粒子速度u_p关系式参数^［

19］可使用吴强^{［参考文献 18

百度学术}18］提出的方法计算得到：

B_{0 K} = ρ_{0 K} C_{0}^{2} \frac{{(1 - φ)}^{2}}{{(1 - S_{0} φ)}^{3}} [1 + (S_{0} - Γ_{0}) φ]

（19）

B_{0 K}^{'} = 4 S_{0} (\frac{1 - φ}{1 - S_{0} φ}) [\frac{1}{1 + (S_{0} - Γ_{0}) φ}] [1 + \frac{1}{2} (S_{0} - Γ_{0} - \frac{Γ_{0}^{2}}{2 S_{0}}) φ + \frac{φ^{2}}{4} (Γ_{0}^{2} - Γ_{0} S_{0})] - 1

（20）

φ = 1 - (V_{0 K} / V_{0})

（21）

对于金属三项式状态方程中的“晶格振动”项，本研究采用基于德拜模型得到的自由能表达式^［

20］：

\begin{array}{l} F_{I o n} (V, T) = \frac{9}{8} N k_{B} Θ_{D} (V) + 3 N k_{B} T l n [1 - e x p (- Θ_{D} (V) / T)] - \\ N k_{B} T d (Θ_{D} (V) / T) \end{array}

（22）

d (x) = \frac{3}{x^{3}} \int_{0}^{x} \frac{y^{3}}{e^{y} - 1} d y

（23）

式中，N=N_A/M，M是金属的摩尔质量，g·mol^-1。Θ_D（V）是德拜温度比容函数，可以从金属Gruneisen系数Γ随比容的变化规律出发，结合陈俊祥^［

21］通过理论分析建立的金属通用Gruneisen系数比容函数推导得到：

Θ_{D} (V) = Θ_{D} (V_{0}) e x p (0.2 Γ_{0} ρ_{0}^{3} V^{3} - 0.3 Γ_{0} ρ_{0}^{2} V^{2} - Γ_{0} ρ_{0} V + 1.1 Γ_{0})

（24）

式中，Θ_D（V₀）是常温常压下金属的德拜温度，K。将式（31）代入式（27），可以得到完善的“晶格振动”项。根据热力学基本关系可以得到“晶格振动”项的压力形式：

p_{I o n} (V, T) = N k_{B} \frac{\partial d [Θ_{D} (V) / T]}{\partial [Θ_{D} (V) / T]} \frac{\partial Θ_{D} (V)}{\partial V} - \frac{9}{8} N k_{B} \frac{\partial Θ_{D} (V)}{\partial V} - \frac{3 N k_{B} e^{- Θ_{D} (V) / T}}{1 - e^{- Θ_{D} (V) / T}} \frac{\partial Θ_{D} (V)}{\partial V}

（25）

对于“电子运动”项，使用基于电子系统费米统计理论得到的金属自由电子气状态方程^［

22］：

p_{E l e} (V, T) = \frac{1}{4} β_{0} ρ_{0} {(\frac{V_{0}}{V})}^{1 / 2} T^{2}

（26）

式中，β₀是金属的电子比热系数，J·kg^-1·K^-1。

式（16）、式（27）与式（28）结合，建立了三项式形式的凝聚态金属产物状态方程Wu‑Chen‑Peng（WCP），铜、银、铅、镍四种金属的WCP状态方程参数如表3所示，其中Γ₀和Θ_D（V₀）取自文献［

22］，β₀取自文献［23］。

表3 WCP状态方程参数值

Table 3 Parameters of WCP equation of state

species	ρ₀ / g·cm^-3	ρ_{0 K}	B_{0 K} / GPa	B′_{0 K}	Γ₀^{［参考文献 22 百度学术}22］	Θ_D（V₀）^［ 22］ / K	β₀^［ 23］ / J·g^-1·K^-2
Cu	8.930	9.068	142.1074	5.05	1.98	321	1.08×10^-5
Ag	10.490	10.675	116.0718	5.21	2.41	231	6.00×10^-6
Pb	11.340	11.642	49.8341	4.86	2.74	87	1.45×10^-5
Ni	8.875	8.979	192.6397	4.76	1.83	450	1.21×10^-4

Note： ρ₀， Γ₀， Θ_D（V₀）， β₀ are the density， Gruneisen coefficient， Debye temperature， and electron specific heat of metals at room temperature and atmospheric pressure； ρ_{0 K}， B_{0 K}， and B′_{0 K} are the density， bulk moduli， and the first‑order partial derivative of the B_{0 K} with respect to pressure of metals at zero temperature and atmospheric pressure.

将以上新建立的VPL和WCP状态方程植入基于最小自由能原理求解爆轰产物平衡组成的热力学程序^［

4］中，利用爆轰波CJ条件和爆轰Hugoniot关系求解炸药爆轰CJ参数，利用爆轰产物波后膨胀的等熵特点计算膨胀过程p‑V关系。同时，本研究直接引入Fried等^{［参考文献 24

百度学术}24］建立的多相碳状态方程描述计算中涉及的凝聚碳产物。

2 应用新状态方程计算炸药爆轰性能

2.1 CHNO单质/混合炸药爆轰性能计算

首先应用新建立的状态方程VPL计算CHNO单质/混合炸药的爆轰性能，为了更全面地评价新状态方程计算准确性，分别计算了不同类型炸药的不同种类爆轰性能参数，包括太安（PETN）的爆速，黑索今（RDX）的爆轰CJ压力，以及熔黑梯‑901（RHT‑901）的爆轰产物做功能力。

2.2.1 PETN爆速计算

PETN的爆轰过程接近理想爆轰^［

4］且产物中凝聚碳含量很少，适合作为气态爆轰产物状态方程准确性检验的“参考物”。由于PETN爆速实验值相对较为丰富，且分布在较宽广的装药密度范围。因此，应用VPL计算了PETN在0.25~1.76 g·cm^-3装药密度范围的爆速D。在将计算结果与实验值^{［参考文献 3

百度学术}3］对比以检验其准确性的同时，为评价VPL相比另两种维里型气态爆轰产物状态方程BKW和VLW在计算准确性上的提高，也将计算结果与文献［3］中使用BKWS，文献［25］中使用VLW以及文献［3］中使用国际著名爆轰热力学计算程序CHEETAH在同样装药密度条件下得到的计算值进行了对比，如图2所示。

a. comparison between calculated and experimental results

b. relative error between calculated and experimental results

图2 PETN不同装药密度下爆速计算值与实验值对比

Fig.2 Comparison between calculated and experimental detonation velocities of PETN under different loading densities

通过图2中的对比不难看出，VPL计算值与实验值之间符合得很好，不同密度下计算结果相对误差均控制在±2.1%以内，最大计算偏差不超过2.5%，计算准确性在不同密度下均明显优于BKWS和VLW。CHEETAH中的JCZ3状态方程参数参考产物热力学数据进行优化，因此能较准确预测PETN在不同密度下爆速，但高密度下VPL计算准确性明显好于CHEETAH。这表明相比CHEETAH，VPL状态方程能够更准确描述气态产物在较高密度下的热力学关系，JCZ3使用固体模型描述产物气体在高密度下的热力学关系，相比VPL使用的真实气体维里方程可能与实际情况存在一定的偏差。

2.2.2 RDX爆轰CJ压力计算

作为一种常用炸药，RDX在不同装药密度下的爆轰CJ压力实验数据较为丰富，使用VPL状态方程计算了RDX在0.56~1.80 g·cm^-3装药密度范围的爆轰CJ压力p_CJ，同样，在将计算结果与实验值^［

3］对比以检验其准确性的同时，为评价VPL相比另两种维里型气态爆轰产物状态方程BKW和VLW在计算准确性上的提高，也将计算结果与文献［3］中使用BKWS，文献［25］中使用VLW以及文献［3］中使用国际著名爆轰热力学计算程序CHEETAH在同样装药密度条件下得到的计算值进行了对比，如图3所示。

a. comparsion between calculated and experimental results

b. relative error between calculated and experimental results

图3 RDX不同装药密度下爆轰CJ压力计算值与实验值对比

Fig.3 Comparison between calculated and experimental detonation CJ pressure of RDX under different loading densities

通过图3中的对比可以看到，VPL计算值与实验值之间符合得较好，相对误差基本维持在±6%以内，计算准确性与VLW和CHEETAH相近且明显优于BKWS。不同于BKWS和VLW，VPL产物势参数未参考RDX爆压实验值标定，因此可以认为VPL能够准确预测炸药在不同密度下的爆轰CJ压力。

2.2.3 RHT‑901爆轰产物做功能力计算

为评价VPL对炸药爆轰产物做功能力计算准确性。使用VPL状态方程计算了典型CHNO混合炸药RHT‑901的爆轰产物等熵膨胀p‑V关系，如图4所示。

图4 应用VPL计算得到的RHT‑901 爆轰产物等熵膨胀p‑V关系

Fig.4 p‑V relationship of isentropic expansion of RHT‑901 detonation products calculated using VPL equation of state

根据所得数据使用基于非线性最小二乘法的SEQS程序拟合了RHT‑901的爆轰产物JWL状态方程参数，如表4所示。

表4 RHT‑901爆轰产物JWL状态方程参数标定结果

Table 4 Parameter calibration results of JWL equation of state of the detonation products of RHT‑901

parameters	ρ₀ / g·cm^-3	A / GPa	B / GPa	R₁	R₂	ω	E₀ / Mbar
value	1.705	297.3745	1.34495	3.3345895	0.6589378	0.3	0.0950146

Note： ρ₀ is loading density of RHT‑901. A， B， R₁， R₂， ω， E₀ are JWL equation of state parameters of RHT‑901.

圆筒试验是评价炸药做功能力的标准方法，使用基于VPL计算值拟合得到的爆轰产物Jones‑Wilkins‑Lee（JWL）状态方程参数，使用非线性动力学仿真软件LS‑DYNA对RHT‑901的标准圆筒试验进行动力学仿真，仿真模型如图5所示。采用二维1/4几何结构，其中炸药与圆筒结构尺寸与文献［

26］记载的结构尺寸相同。圆筒内径25 mm，外径30.12 mm，长度300 mm；炸药尺寸Φ25 mm×300 mm。圆筒材料为纯铜，使用Johnson‑Cook材料模型和Gruneisen状态方程，参数取自文献［27］。

图5 RHT‑901标准圆筒试验动力学仿真模型结构图

Fig.5 Simulation geometric structure model of the standard cylinder test of RHT‑901

使用LS‑PrePost软件对仿真结果进行后处理，输出圆筒壁速度曲线与文献［

26］记载的实验值以及文献［4］记载的基于BKW和VLW状态方程得到的仿真结果进行了对比，如图6所示，其中t=20 μs时刻圆筒壁速度u_t₌₂₀的实验值与仿真值对比如表5所示。

图6 RHT‑901标准圆筒试验圆筒壁速度历程仿真结果与实验值对比

Fig.6 Comparison between calculated and experimental values of cylindrical wall velocity history in standard cylinder test of RHT‑901

表5 RHT‑901标准圆筒试验t=20 μs时刻圆筒壁速度仿真值与实验值对比

Table 5 Comparison between calculated and experimental values of cylindrical wall velocity in t=20 μs in standard cylinder test of RHT‑901

type	reference	u_t₌₂₀ / m·s^-1	relative error / %
experimental results	［26］	1676	0
VPL	in this study	1661	-0.895
VLW	［4］	1711	2.088
BKW	［4］	1861	11.038

从图5和表5中的对比可以看到，基于VPL得到的仿真结果与实验值间符合得更好，t=20 μs时刻筒壁速度仿真值与实验值之间相差不到1%。基于BKW得到的筒壁速度明显高于实验值，基于VLW得到的仿真结果好于BKW，但依然不及VPL，这表明VPL能更准确计算炸药爆轰产物在波后的膨胀做功过程。

2.2 含能金属盐爆轰性能计算

2.2.1 含能金属盐的爆速计算

使用新爆轰产物状态方程VPL与WCP对叠氮化铅、斯蒂芬酸铅、硝酸肼镍、5‑硝基四唑亚铜、乙炔银‑硝酸银五种含能金属盐在不同密度下的爆速D进行了计算，并与相关文献记载的实验值^［

28-30］进行了对比，如图7所示，其中点为实验值，线为计算值。从图7可以看到，虽然相关文献记载的实验值有一定分散性，但爆速‑密度关系计算结果与实验值在趋势上符合得较好。除某些分散很大的实验值，计算结果与实验值之间相对误差基本维持在±4.5%以内，这表明VPL与WCP能准确预测不同种类含能金属盐在不同密度下的爆速。

图7 五种含能金属盐爆速计算结果与实验值对比

Fig.7 Comparison between calculated and experimental values of detonation velocity of five energetic metal salts

LA—lead azide， LSt—lead styphnate， NHN—nickel hydrazine nitrate， CuNT—copper（I）‑5‑nitrotetrazolate， AgCN—silver acetylide‑silver nitrate

2.2.2 叠氮化铅爆轰CJ参数计算

不难看到图7中的爆速实验值具有一定分散性，Mu^［

31］应用微尺度瞬态测试技术对叠氮化铅在装药密度3.571 g·cm^-3下的爆轰参数进行了精确测试，得到一组准确的爆速与爆轰CJ压力数据，将计算结果与其对比能够更准确而全面检验VPL与WCP对含能金属盐爆轰性能计算准确性。因此，使用方程VPL和WCP计算叠氮化铅在3.571 g·cm^-3下的爆速D与爆轰CJ压力p_CJ，并与文献［31］中记载的实验值进行了对比，如表6所示。可以看到，采用VPL与WCP能够准确预测叠氮化铅的爆轰CJ参数，计算结果与文献［31］中采用高精度测试得到的实验值间相差不到4%，这进一步表明VPL与WCP能够准确预测含能金属盐的爆轰CJ参数。

表6 叠氮化铅（3.571 g·cm^-3）爆速与爆轰CJ压力计算值与实验值对比

Table 6 Comparison between calculated and experimental values of detonation velocity and detonation CJ pressure of lead azide（3.571 g·cm^-3）

type	experimental results^{［参考文献 31 百度学术}31］	calculated results	relative error / %
D	4726 m·s^-1	4807 m·s^-1	1.714
p_CJ	17.12 GPa	17.78 GPa	3.855

2.2.3 铜叠氮化物爆轰性能计算

使用VPL与WCP计算了铜叠氮化物CuN_4.432₉^［

32］的爆轰CJ参数与爆轰产物等熵膨胀p‑V关系，并标定了爆轰产物JWL状态方程参数，如表7所示。

表7 CuN_4.4329的爆轰CJ参数与爆轰产物JWL状态方程参数计算结果

Table 7 Calculation results of detonation CJ parameters and JWL state equation parameters of detonation products for CuN_4.4329

parameters	ρ₀ / g·cm^-3	D / m·s^-1	p_CJ / GPa	A / GPa	B / GPa	R₁	R₂	ω	E₀ / Mbar
value	2.215	5640	15.36	179.11568	1.10224	3.319194	0.6540107	0.19	0.083062

Note： ρ₀， D， p_CJ are loading density， detonation velocity， and detonation CJ pressure of CuN_4.4329； A， B， R₁， R₂， ω， E₀ are JWL equation of state parameters of CuN_4.4329.

基于计算值使用LS‑DYNA对CuN_4.4329爆轰驱动飞片过程进行动力学仿真，仿真模型如图8所示。采用二维1/4几何结构，装药与飞片结构及尺寸与文献［

32］中相同，装药尺寸Φ1.0 mm×1.0 mm，飞片厚度0.05 mm。飞片材质为纯钛，使用Johnson‑Cook材料模型和Gruneisen状态方程，参数取自文献［27］。

图8 CuN_4.4329爆轰驱动飞片动力学仿真模型结构图

Fig.8 Simulation geometric structure model of the detonation driven flyer test for CuN_4.4329

使用LS‑PrePost软件对仿真结果进行后处理，输出飞片速度曲线并与文献［

32］记载的实验值进行对比，如图9所示，0.4~0.5 μs间飞片平均速度u_average的实验与仿真值对比如表8所示。可以看到，仿真结果与实验曲线之间在飞片运动各个阶段均符合得很好，飞片稳定速度仿真值与实验值相差不到1%，这表明新状态方程WCP与VPL能够准确评价含能金属盐的爆轰驱动能力。

图9 CuN_4.4329爆轰驱动飞片速度历程仿真结果与实验值对比

Fig.9 Comparison between calculated and experimental values of detonation driven flyer velocity history for CuN_4.4329

表8 飞片的材料模型及状态方程参数

Table 8 Material constitutive model and equation of state parameters of flyer

type	experimental results^{［参考文献 32 百度学术}32］	calculated results	relative error / %
u_average	2483 m·s^-1	2459 m·s^-1	-0.967
type	experimental results^{［参考文献 32 百度学术}32］	calculated results	relative error / %
u_average	2483 m·s^-1	2459 m·s^-1	-0.967

3 结论

研究通过理论计算得到了Exp‑6势在T^*=1~300无量纲温度范围内的3~5阶维里系数，对其拟合建立了新气态产物状态方程Virial‑Peng‑Long（VPL），并通过引入更准确的“冷压”项和考虑金属高压热运动变化的“晶格振动”项建立了新的凝聚态金属产物三项式状态方程Wu‑Chen‑Peng（WCP），并将其应用到典型CHNO单质/混合炸药的爆轰参数，和典型含能金属盐的爆轰CJ参数与驱动做功能力计算中。得到结果如下：

（1）使用VPL计算了PETN、RHT‑901等几种典型CHNO单质/混合炸药的爆轰参数，通过对比表明通过VPL能够准确评价CHNO单质/混合炸药爆轰性能。对PETN爆速计算偏差在±2.1%以内，最高不超过2.5%；对RDX爆轰CJ压力预测相对误差在±6%以内；对RHT‑901圆筒实验筒壁稳定速度预测偏差绝对值<1%。计算准确性相比BKW和VLW两种状态方程有了显著提高。

（2）使用VPL与WCP评价了叠氮化铅、硝酸肼镍等几种典型含能金属盐的爆轰CJ参数与驱动做功能力，通过对比表明基于VPL与WCP能够准确预测含能金属盐爆轰性能。其中对叠氮化铅爆轰CJ参数计算偏差不超过4%；对铜叠氮化物驱动飞片速度预测相对误差在1%以内。

本研究建立的新爆轰产物状态方程VPL与WCP可用于CHNO单质/混合炸药与含能金属盐爆轰参数的热力学计算。为进一步实现VPL与WCP在爆炸动力学仿真中的直接应用，提高仿真的准确性，下一步有待进行爆轰热力学计算与动力学仿真的耦合研究。

参考文献

吴雄. VLW爆轰产物状态方程的发展及应用［J］. 火炸药学报， 2021， 44（1）： 1-7. [百度学术]

WU Xiong. Development and application of VLW equation of state for detonation products［J］. Chinese Journal of Explosives & Propellants， 2021， 44（1）： 1-7. [百度学术]

SUCESKA M， ANG H G， STIMAC B， et al. BKW EOS： History of modifications and further improvement of accuracy with temperature‑dependent covolumes of polar molecules［J］. Propellants， Explosives， Pyrotechnics， 2023， 48（1）： e202100278. [百度学术]

HOBBS M L， SCHMITT R G， MOFFAT H K. JCZS3‑an improved database for EOS calculations［R］. SAND2018‑6389C， 2018. [百度学术]

韩勇. 炸药爆轰产物状态方程理论与应用研究［D］. 绵阳：中国工程物理研究院， 2015. [百度学术]

HAN Yong. Theoretical and applied research on the equation of state for detonation products［D］. Mianyang： China Academy of Engineering Physics， 2015. [百度学术]

韩勇，郭向利，龙新平. 高温高压CO₂状态方程研究［J］. 含能材料， 2016， 24（5）： 462-468. [百度学术]

HAN Yong， GUO Xiang‑li， LONG Xin‑ping. High tempertaure and high pressure equation of state of carbon dioxide［J］. Chinese Journal of Energetic Materials （Hanneng Cailiao）， 2016， 24（5）： 462-468. [百度学术]

赵雅琦. 短脉冲加载下小尺寸装药二维冲击起爆与爆轰特性［D］. 北京：北京理工大学， 2022. [百度学术]

ZHAO Ya‑qi. Two‑dimensional shock initiation characteristics of small‑size charges under short‑pulse loading［D］. Beijing： Beijing Institute of Technology， 2022. [百度学术]

覃文志，龙新平，蒋小华. 金属配合物类炸药的爆轰性能计算及数值模拟［J］. 含能材料， 2011， 19（5）： 540-543. [百度学术]

QIN Wen‑zhi， LONG Xin‑ping， JIANG Xiao‑hua. Calculation of detonation parameters of metal compounds explosives［J］. Chinese Journal of Energetic Materials （Hanneng Cailiao）， 2011， 19（5）： 540-543. [百度学术]

BAI J F， ZHANG P， ZHOU C W， et al. Theoretical studies of real‑fluid oxidation of hydrogen under supercritical conditions by using the virial equation of state［J］. Combustion and Flame， 2022， 243： 111945. [百度学术]

韩勇，龙新平，黄毅民，等. L‑J Exp‑6两种形式势能函数对计算无量纲第二维里系数的影响［J］. 含能材料， 2009， 17（5）： 574-577. [百度学术]

HAN Yong， LONG Xin‑ping， HUANG Yi‑min， et al. Effects of L‑J and Exp‑6 potential function on calculation of reduced second virial coefficient［J］. Chinese Journal of Energetic Materials （Hanneng Cailiao）， 2009， 17（5）： 574-577. [百度学术]

KIHARA B T. Determination of intermolecular forces from the equation of state of gases［J］. Journal of the Physical Society of Japan， 1948， 3（4）： 265-268. [百度学术]

BARKER J A， LEONARD P J， POMPE A. Fifth virial coefficients［J］. The Journal of Chemical Physics， 1966， 44（11）： 4206-4211. [百度学术]

BELONOSHKO A， SAXENA S K. A molecular dynamics study of the pressure‑volume‑temperature properties of supercritical fluids： II. CO₂， CH₄， CO， O₂， and H₂［J］. Geochimica et Cosmochimica Acta， 1991， 55（11）： 3191-3208. [百度学术]

DUAN Z H， MOLLER N， WEARE J H. Molecular dynamics simulation of PVT properties of geological fluids and a general equation of state of nonpolar and weakly polar gases up to 2000 K and 20000 bar［J］. Geochimica et Cosmochimica Acta， 1992， 56（10）： 3839-3845. [百度学术]

ZAUG J M， CARTER J A， BASTEA S， et al. Experimental measurement of speeds of sound in dense supercritical carbon monoxide and development of a high‑pressure， high‑temperature equation of state［J］. The Journal of Physical Chemistry B， 2013， 117（18）： 5675-5682. [百度学术]

MAZEVET S， BLOTTIAU P， KRESS J D， et al. Quantum molecular dynamics simulations of shocked nitrogen oxide［J］. Physical Review B， 2004， 69（22）： 224207. [百度学术]

SPAN R， LEMMON E W， JACOBSEN R T， et al. A reference quality equation of state for nitrogen［J］. Internation Journal of Thermophysics， 1998， 19（4）： 1121-1132. [百度学术]

ZHANG Z G， DUAN Z H. Prediction of the PVT properties of water over wide range of temperatures and pressures from molecular dynamics simulation［J］. Physics of the Earth and Plantary Interiors， 2005， 149（3‑4）： 335-354. [百度学术]

吴强. 金属材料高压物态方程及Gruneisen系数的研究［D］. 绵阳：中国工程物理研究院， 2004. [百度学术]

WU Qaing. Research on the high pressure equation of state and Gruneisen coefficient of metals［D］. Mianyang： China Academy of Engineering Physics， 2004. [百度学术]

MARSH S P. LASL shock Hugoniot data［M］. Berkeley： University of California Press， 1980. [百度学术]

GREEFF C W， BOETTGER J C， GRAF M J， et al. Theoretical investigation of the Cu EOS standard［J］. Journal of Physics and Chemistry of Solids， 2006， 67（9）： 2033-2040. [百度学术]

陈俊祥，于继东，李平，等. Gruneisen γ通用函数及完全物态方程［J］，物理学报， 2015， 64（8）： 086401. [百度学术]

CHEN Jun‑xiang， YU Ji‑dong， LI Ping， et al. Universal function of Gruneisen γ and the complete equation of state［J］. Acta Physica Sinica， 2015， 64（8）： 086401. [百度学术]

ALTSHULER L V， BRUSNIKIN S E， KUZMENKOV E A. Isotherms and Gruneisen functions for 25 metals［J］. Journal of Applied Mechanics and Technical Physics， 1987， 28（1）： 129-141. [百度学术]

经福谦. 实验物态方程导引（第二版）［M］. 北京：国防工业出版社， 1999： 122-128. [百度学术]

JING Fu‑qian. Introduction of experimental equation of state （Second Edition）［M］. Beijing： National Defense Industry Press， 1999： 44-46. [百度学术]

FRIED L E， HOWARD W M. Explicit Gibbs free energy equation of state applied to the carbon phase diagram［J］. Physical Review B. 2000， 61（13）： 8734-8743. [百度学术]

刘琴. 炸药爆轰过程中凝聚碳相态演化与能量释放的关联机制［D］. 合肥：中国科学技术大学， 2022. [百度学术]

LIU Qin. Correlation mechanism between phase evolution of condensed carbon and energy release during explosive detonation［D］. Mianyang： China Academy of Engineering Physics， 2011. [百度学术]

董海山，周芬芬. 高能炸药及相关物性能［M］. 北京：科学出版社， 1989： 150. [百度学术]

DONG Hai‑shan， ZHOU Fen‑Fen. Performance of high explosives and related substances［M］. Beijing： Science Press， 1989： 150. [百度学术]

辛春亮，薛再清，涂建，等. 有限元分析常用材料参数手册［M］. 北京：机械工业出版社， 2018： 124-138. [百度学术]

XIN Chun‑liang， XUE Zai‑qing， TU Jian， et al. Handbook of common material parameters for finite element analysis［M］. Beijing： China Machine Press， 2018： 124-138. [百度学术]

JAFARI M， KESHAVARZ M H， ZAMANI A， et al. A novel method for assessment of the velocity of detonation for primary explosives［J］. Propellants， Explosives， Pyrotechnics， 2018， 43（4）： 342-347. [百度学术]

TAPPAN A S， BALL J P， MILLER J C. DBX‑1 （copper（I）‑5‑nitrotetrazolate） reactions at sub‑millimeter diameters［R］. SAND2013‑4254C， 2013. [百度学术]

曹仕瑾.叠氮肼镍的结构与性能［D］.南京：南京理工大学，2007. [百度学术]

CAO Shi‑jin. Structure and properties of nickel azide hydrazine［D］. Nanjing： Nanjing University of Science and Technology， 2007. [百度学术]

MU Y F， ZHANG W， SHEN R Q， et al. Observations on detonation growth of lead azide at microscale［J］. Micromachines， 2022， 13（3）： 451. [百度学术]

PENG Y， YU Q X， YIN Q， et al. Performance investigation of a copper azide micro detonator using experiments and simulations［J］. Propellants， Explosives， Pyrotechnics， 2022， 47（4）： e202100343. [百度学术]

Home

About Journal

For Authors

For Reviewers

Ethics Policy

Contact Us

EMF Journal

一种新爆轰产物状态方程及其在炸药爆轰性能预测上的应用 PDF

摘要

图文摘要

关键词

0 引言

1 新爆轰产物状态方程

1.1 气态产物维里状态方程Virial‑Peng‑Long

1.2 凝聚态金属产物三项式状态方程Wu‑Chen‑Peng

2 应用新状态方程计算炸药爆轰性能

2.1 CHNO单质/混合炸药爆轰性能计算

2.2 含能金属盐爆轰性能计算

3 结论

参考文献

Home

About Journal

For Authors

For Reviewers

Ethics Policy

Contact Us

EMF Journal

一种新爆轰产物状态方程及其在炸药爆轰性能预测上的应用 PDF

摘要

图文摘要

关键词

0 引 言

1 新爆轰产物状态方程

1.1 气态产物维里状态方程Virial‑Peng‑Long

1.2 凝聚态金属产物三项式状态方程Wu‑Chen‑Peng

2 应用新状态方程计算炸药爆轰性能

2.1 CHNO单质/混合炸药爆轰性能计算

2.2 含能金属盐爆轰性能计算

3 结 论

参考文献

0 引言

3 结论