炸药装药在发射环境下的安全性问题，一直制约着高能炸药在大口径榴弹中的应用。炸药装药在膛内早炸实质是多种因素共同作用的一个复杂过程。大量实验和研究结果表明，装药质量是影响炸药发射安全性的最主要原因，即装药的底隙、裂纹、气泡、缩孔等是导致装药膛炸的主要因素^[1]。因此，避免炸药装药发生早炸，研究炸药装药在发射环境下的安全性问题具有十分重要的意义。

炸药装药发射安全性研究一直处于发展完善阶段^[2-3]。学者们主要是从实验和数值仿真2个方面来研究。美国的学者们最早建立了用火药驱动的膛压模拟装置，系统研究了装药对底隙要求的安全性限定条件^[4-6]。建立了膛压式试验模拟装置和试验方法；分析炸药装药缺陷对炸药性能的影响^[7-8]；研究炸药装药在发射过程中的本构关系和结构响应等^[9]。由于现有实验条件有限，从微观层次上理解炸药内部的变形机制还比较困难，只能通过建立各种近似的本构模型来解释炸药的变形规律。由Bennett等^[10]建立的粘弹性裂纹本构模型是研究含能材料损伤的代表性细观损伤模型，该模型基于裂纹的演化、扩展机制，其特点是在裂纹体基础上耦合了广义粘弹性体，描述了炸药的粘弹性力学响应，并采用微裂纹体描述其损伤的演化过程。由于该模型能更贴近实际地模拟炸药的力学性能，且Dienes^[11]等已经验证了统计裂纹模型计算摩擦生热引起的温升能够引起炸药的点火反应。本文采用粘弹性裂纹本构模型，基于有限体积法自编程序模拟了炸药在弹底压力作用下的动力学响应，研究了炸药的点火温度及起爆的临界弹底压力，得到的结论可为炸药的安全性评估提供参考。

1 粘弹性统计裂纹模型

粘弹性统计裂纹模型(viscoelastic statistical cracking model)是在Addessio等研究工作的基础上发展而来的^[12]。该模型是由2个模型组成，一个是Maxwell研究的由n个粘弹性体元组成的粘弹性体模型；另一个是由Dienes^[10]研究的统计裂纹力学模型(statistical crack mechanics，SCRAM)。SCRAM模型本质上是应用细观力学研究易碎材料的大变形和裂纹损伤。在变形的过程中，假设裂纹是随机分布的。

采用张量的形式，由运动方程给出固体的应变率：

$ {{\mathit{\boldsymbol{\dot \varepsilon }}}_{ij}} = \frac{1}{2}\left( {\frac{{\partial {{\mathit{\boldsymbol{\dot u}}}_j}}}{{\partial {\mathit{\boldsymbol{x}}_i}}} + \frac{{\partial {{\mathit{\boldsymbol{\dot u}}}_i}}}{{\partial {\mathit{\boldsymbol{x}}_j}}}} \right) $

(1)

式中：${\mathit{\boldsymbol{\dot \varepsilon }}_{ij}} $为应变率；${\mathit{\boldsymbol{\dot u}}_j} $是材料的变形率。总的应变率可以分解成以下形式：

$ {{\mathit{\boldsymbol{\dot \varepsilon }}}_{ij}} = {{\mathit{\boldsymbol{\dot e}}}_{ij}} + {\mathit{\boldsymbol{\delta }}_{ij}}{{\mathit{\boldsymbol{\dot \varepsilon }}}_m} $

(2)

式中：${\mathit{\boldsymbol{\dot e}}_{ij}} $是偏应变率；δ_ij是Kronecker符号；${\mathit{\boldsymbol{\dot \varepsilon }}_m} $是平均应变率，${\mathit{\boldsymbol{\dot \varepsilon }}_m} = \frac{1}{3}{\mathit{\boldsymbol{\dot \varepsilon }}_{ii}} $。

在Maxwell的粘弹性体模型中，弹性元的偏应力与偏应变的关系为S_ij=2Ge_ij，粘性元的偏应力与偏应变率的关系为S_ij=2η ${\mathit{\boldsymbol{\dot e}}_{ij}} $。粘弹性体的偏应力率与偏应变率的关系为：

$ {{\mathit{\boldsymbol{\dot S}}}_{ij}} = 2\mathit{\boldsymbol{G\dot e}}_{ij}^{ev} - \frac{{{\mathit{\boldsymbol{S}}_{ij}}}}{\mathit{\boldsymbol{\tau }}} $

(3)

式中：${\mathit{\boldsymbol{\dot S}}_{ij}} $为偏应力率；τ=η/G，为松弛时间；G为弹性元的剪切模量；η为粘性元的粘性系数；$\mathit{\boldsymbol{\dot e}}_{ij}^{ev} $为粘弹性体的偏应变率。粘弹性体一般由N个粘弹性体元组成，粘弹性体的偏应力为N个粘弹性体元的偏应力之和，为${\mathit{\boldsymbol{S}}_{ij}} = \sum\limits_{n = 1}^N {\mathit{\boldsymbol{S}}_{ij}^{\left( n \right)}} $，则粘弹性体的偏应力率与偏应变率之间的关系为：

$ {{\mathit{\boldsymbol{\dot S}}}_{ij}} = \sum\limits_{n = 1}^N {\left( {2{\mathit{\boldsymbol{G}}^{\left( n \right)}}\mathit{\boldsymbol{\dot e}}_{ij}^{ev} - \frac{{\mathit{\boldsymbol{S}}_{ij}^n}}{\mathit{\boldsymbol{\tau }}}} \right)} $

(4)

式中：G⁽ⁿ⁾为第n个弹性元的剪切模量；S_ij⁽ⁿ⁾为第n个粘弹性体的偏应力。

在粘弹性裂纹本构模型中，偏应变率由粘弹性体的偏应变率和微裂纹体的偏应变率组成，即：

$ {{\mathit{\boldsymbol{\dot e}}}_{ij}} = \mathit{\boldsymbol{\dot e}}_{ij}^{ev} + \mathit{\boldsymbol{\dot e}}_{ij}^c $

(5)

Addessio等在1990年的研究工作中得到了微裂纹体的偏应变率和偏应力之间的关系:

$ \mathit{\boldsymbol{\dot e}}_{ij}^c = {\mathit{\boldsymbol{\beta }}^e}{\mathit{\boldsymbol{c}}^3}{\mathit{\boldsymbol{S}}_{ij}} $

(6)

式中：c为平均裂纹半径; β^e是与剪切模量和裂纹的初始分布N₀有关的参数，即:

$ 2\mathit{\boldsymbol{G}}{\mathit{\boldsymbol{\beta }}^e} = \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{N}}_0} \equiv \frac{1}{{{a^3}}} $

(7)

式中：a为初始裂纹尺寸。偏应变与裂纹半径的关系为:

$ 2\mathit{\boldsymbol{Ge}}_{ij}^c = {\left( {\frac{c}{a}} \right)^3}{\mathit{\boldsymbol{S}}_{ij}} $

(8)

而偏应变率与裂纹半径的关系为:

$ 2\mathit{\boldsymbol{G\dot e}}_{ij}^c = 3{\left( {\frac{c}{a}} \right)^2}{\mathit{\boldsymbol{S}}_{ij}} + {\left( {\frac{c}{a}} \right)^3}{{\mathit{\boldsymbol{\dot S}}}_{ij}} $

(9)

由式(4)、(5)、(9)可得到粘弹性裂纹本构模型的偏应力率的表达式为：

$ {{\mathit{\boldsymbol{\dot S}}}_{ij}} = \frac{{2\mathit{\boldsymbol{G}}{{\mathit{\boldsymbol{\dot e}}}_{ij}} - \sum\limits_{n = 1}^N {\frac{{\mathit{\boldsymbol{S}}_{ij}^{\left( n \right)}}}{{{\mathit{\boldsymbol{\tau }}^{\left( n \right)}}}}} - 3{{\left( {\frac{c}{a}} \right)}^2}{\mathit{\boldsymbol{S}}_{ij}}}}{{1 + {{\left( {\frac{c}{a}} \right)}^3}}} $

(10)

式中：$\mathit{\boldsymbol{\dot S}}_{ij}^{\left( n \right)} $为每个粘弹性体元的偏应力率；$\mathit{\boldsymbol{\dot S}}_{ij}^{\left( n \right)} $的表达式为：

$ \begin{array}{l} \mathit{\boldsymbol{\dot S}}_{ij}^{\left( n \right)} = 2{\mathit{\boldsymbol{G}}^{\left( n \right)}}{{\mathit{\boldsymbol{\dot e}}}_{ij}} - \frac{{\mathit{\boldsymbol{S}}_{ij}^{\left( n \right)}}}{{{\mathit{\boldsymbol{\tau }}^{\left( n \right)}}}} - \frac{{{\mathit{\boldsymbol{G}}^{\left( n \right)}}}}{\mathit{\boldsymbol{G}}}\left[ {3{{\left( {\frac{c}{a}} \right)}^2}\frac{{\dot c}}{a}{\mathit{\boldsymbol{S}}_{ij}} + } \right.\\ \;\;\;\;\;\;\;\;\left. {{{\left( {\frac{c}{a}} \right)}^3}{{\mathit{\boldsymbol{\dot S}}}_{ij}}} \right] \end{array} $

(11)

式中：$\mathit{\boldsymbol{G = }}\sum\limits_{n = 1}^N {{\mathit{\boldsymbol{G}}^{\left( n \right)}}} $。

Dienes在研究工作中描述了裂纹的增长速率，假设裂纹的增长速率与应力强度有关，即

$ \left\{ \begin{array}{l} \dot c = {\mathit{\boldsymbol{v}}_{\max }}{\left( {\frac{K}{{{K_1}}}} \right)^m}, \;\;\;\;\;K < K'\\ \dot c = {\mathit{\boldsymbol{v}}_{\max }}\left( {1 - {{\left( {\frac{{{K_0}}}{K}} \right)}^2}} \right), \;\;\;\;\;K \ge K' \end{array} \right. $

(12)

式中: $K\mathit{' = }{\mathit{K}_0}\sqrt {1 + \left( {\frac{2}{m}} \right)} , $ ${K_1}\mathit{ = }{\mathit{K}_0}\sqrt {1 + \left( {\frac{2}{m}} \right)} \cdot {\left[ {1 + \frac{m}{2}} \right]^{\frac{1}{m}}} $。

应力张量可分解为静水压力和偏应力张量，即:

$ {\mathit{\boldsymbol{\sigma }}_{ij}} = {\mathit{\boldsymbol{S}}_{ij}} + \mathit{\boldsymbol{p}}{\mathit{\boldsymbol{\delta }}_{ij}} $

(13)

式中：σ_ij为应力张量, 偏应力张量通过粘弹性裂纹本构模型更新。

2 状态方程

研究的过程中，炸药的压力通过材料的状态方程更新。炸药的固态状态方程一般采用Shock状态方程或者JWL状态方程。本文采用固态的JWL状态方程:

$ p = A\left( {1 - \frac{w}{{{R_1}V}}} \right){{\rm{e}}^{ - {R_1}V}} + B\left( {1 - \frac{w}{{{R_2}V}}} \right){{\rm{e}}^{ - {R_2}V}} + \frac{{wE}}{V} $

(14)

式中：p为压力；E为单位初始体积内能；V为相对体积；A、B、R₁、R₂、w为材料参数，具体取值见表 1^[14]。

3 宏观体积加热模型

为了预测炸药在外界刺激下反应过程温度场即描述整体热温度场的行为，文中引入了粘弹性统计裂纹模型中的宏观连续尺度内的体积加热模型。在粘弹性裂纹本构模型中，宏观控制体积单元的加热包括了描述绝热体积压缩，粘性，裂纹损伤及化学分解热等，基于Arrhenius一级化学反应动力学方程，忽略整体热传导项，整体温度的变化为：

$ \mathit{\boldsymbol{\dot T}} = - \gamma T{{\mathit{\boldsymbol{\dot \varepsilon }}}_{ii}} + \frac{I}{{\rho {C_V}}}\left( {{{\mathit{\boldsymbol{\dot W}}}_{ve}} + {{\mathit{\boldsymbol{\dot W}}}_{cr}}} \right) + \frac{{\Delta H}}{{{C_V}}}Z\exp \left( { - \frac{E}{{RT}}} \right) $

(15)

式中：T为温度；γ代表Gruneisen系数；C_V是定容比热；${\dot W_{ve}} $为粘性功率；${\dot W_{cr}} $为裂纹损伤功率；ΔH、E、Z为含能材料热分解反应单位质量分解热、活化能和指前因子。

4 炸药在弹底压力作用下的温升计算

炸药装药模型如图 1所示，炸药装药模型总长56.7 cm，最大半径6.3 cm。

炸药模型采用PBX9501炸药，它是一种固体的复合炸药，它的主要成份是HMX炸药，其性能参数，密度为1 810 kg/m³, C_v为1 200 J·kg^-1·K^-1，ΔH为5.5×10⁶，Z为5×10⁹，E/R的值为2.652×10⁴。由于粘弹性统计裂纹模型涉及的参数较多，目前现有文献中只给出了PBX9501炸药的相关参数，并做了相应的研究工作，说明了PBX9501炸药适合采用粘弹性统计裂纹本构模型进行分析。采用粘弹性裂纹本构模型进行计算时，选用5个粘弹性体元并联的粘弹性体，再与统计裂纹模型结合。PBX9501炸药的粘弹性裂纹本构模型参数见文献[10-13]，如表 2、3、4所示。

在炸药模型的底部施加模拟的弹底压力曲线，其分布曲线如图 2所示。在计算时，假设炸药模型初始温度为300 K。

基于粘弹性裂纹本构模型，采用有限体积法自编程序计算了炸药装药模型在弹底压力作用下的动力学响应。计算结果如图 3所示。

从计算结果可以看出，炸药模型在弹底压力作用下，炸药内部的整体温度不断升高，应力也不断变化，且应力的变化趋势与施加的弹底压力的变化曲线基本是一致的，说明计算结果是可信的。在计算中考虑了炸药的裂纹瑕疵，由于微裂纹的存在，炸药在弹底压力的作用下会产生裂纹损伤功，根据宏观体积加热模型可知，裂纹损伤功会使炸药的整体温度升高。结果显示炸药模型前部的最高温度为300.6 K，最大应力为2.1 MPa；炸药中部的最高温度在342 K，最大应力为113 MPa；炸药底部的最高温度为379 K，最大应力为202 MPa，这说明在弹底压力的作用下，整个炸药内部的温度不断的在变化，由于弹底压力施加在炸药的底部，所以炸药的底部温度升高的较快，然后温升慢慢从炸药的底部传到炸药的中部及前部，因此炸药的前部温度升高较为缓慢。

如果炸药模型不变，不断的增大弹底压力，炸药内部的整体温度也会不断的升高。当弹底压力增大到一定程度时，炸药底部的温度也会增加到一定值，当炸药底部温度达到了炸药的点火温度时，炸药最终发生了点火反应。采用Arrhenius一级化学反应动力学方程后，点火温度的判据就是dT/dt→∞。为了确定炸药起爆的临界弹底压力，在计算时仍采用图 2的弹底压力曲线，然后成倍的增加弹底压力的值，计算了炸药的起爆情况，计算结果如表 6所示。

图 9为炸药在初始压力为0.505 MPa，峰值压力为1.635 GPa的弹底压力作用下，炸药底部的温度变化曲线。如图可知，计算时间在16 ms之前，温度处于较为平稳的增长状态，计算时间在16 ms之后时，温度突然呈直线上升，且升高到无穷大，说明炸药的温度达到了点火温度，炸药发生了点火反应，从而确定了炸药起爆的临界弹底压力。从图中还可以看出，温度曲线突然上升的阶段，炸药温度大致在700~800 K，说明炸药的点火温度应该在700~800 K。

图 10为炸药在初始压力为0.505 MPa，峰值压力为1.635 GPa的弹底压力作用下，炸药温度的变化云图。如图所示，炸药底部的温度已经超过了800 K，也说明了炸药发生了点火反应。

5 结论

1) 在计算中采用粘弹性统计裂纹模型，由于微裂纹的存在，炸药在弹底压力的作用下会产生裂纹损伤功，根据宏观体积加热模型可知，裂纹损伤功会使炸药的整体温度升高，能够更准确地描述炸药的力学行为。

2) 炸药模型在弹底压力的作用下，炸药内部的应力变化与弹底压力的变化趋势大体是一致的，说明计算结果是正确的。由于弹底压力施加在炸药的底部，炸药底部的温度升高较为明显，而炸药中部和前部的温度相对增长较为缓慢。

3) 计算了在不断增长的弹底压力作用下炸药的起爆情况，初步估算了炸药的点火温度在700~800 K，确定了炸药起爆的临界弹底压力。研究结果能够为炸药的安全性评估提供依据。