固体炸药爆轰与惰性介质相互作用的一种扩散界面模型

于明

doi:10.11883/bzycj-2019-0435

固体炸药爆轰与惰性介质相互作用的一种扩散界面模型

doi: 10.11883/bzycj-2019-0435

于明^{1, 2,}

1.
北京大学应用物理与技术中心，北京 100871
2.
北京应用物理与计算数学研究所计算物理重点实验室，北京 100083

基金项目: 国家自然科学基金（11772066，11272064）；国防基础科研核科学挑战专题（TZZT2016002）；中国工程物理研究院创新发展基金（CX2019026）

详细信息

作者简介:
于　明（1971－　），男，博士，研究员，yu_ming@iapcm.ac.cn

中图分类号: O381
计量
- 文章访问数: 4212
- HTML全文浏览量: 1432
- PDF下载量: 69
- 被引次数: 1
出版历程
- 收稿日期: 2019-11-18
- 修回日期: 2020-06-12
- 网络出版日期: 2020-09-25
- 刊出日期: 2020-10-05

An improved diffuse interface model for the numerical simulation of interaction between solid explosive detonation and inert media

YU Ming^{1, 2
,}

1.
Center for Applied Physics and Technology, Peking University, Beijing 100871, China
2.
Key Laboratory for Computational Physics, Institute of Applied Physics and Computational Mathematics, Beijing 100083, China

摘要

摘要: 提出一种保持热力学一致性的扩散界面模型，用来数值模拟固体炸药爆轰与惰性介质的相互作用问题。基于混合网格内各组分物质间可以达到力学平衡状态而不能达到热学平衡状态的假设，由混合网格能量守恒以及压力相等条件，推导出每种组分物质的体积分数演化方程。由此获得的扩散界面模型包括组分物质的质量守恒方程、混合物质的动量及总能量守恒方程，同时包括组分物质的体积分数演化方程和混合物质的压力演化方程。该扩散界面模型的主要特点是考虑了化学反应以及热学非平衡的影响。提出的扩散界面模型在物质界面附近不会出现物理量的非物理振荡现象、适用于任意表达形式的物质状态方程以及任意数目的惰性介质。
- 扩散界面模型 /
- 固体炸药爆轰 /
- 多物质流动 /
- 热学非平衡 /
- 多相流模型
Abstract: In the article a thermodynamically consistent diffuse interface model is proposed in order to numerically simulate the interaction between solid explosive detonation and compressible inert media. The chemical reaction of detonation course in solid explosive is simplified as the solid-phase reactant changing into the gas-phase product, thus the mixture within a control volume is regarded to be composed by three kinds of components: solid-phase reactant, gas-phase product and inert media, and all components are thought to be in mechanical equilibrium and thermal nonequilibrium because of their having distinct thermodynamic properties or equations of state. The starting point based on the energy conservation of the mixture and pressure equivalence among components within the control volume is adopted to derive the evolution equation for volume fraction of every component, in which the equation for energy conservation of the mixture is decomposed into a family of equations for energy conservation of the all components with the heat exchange resulting from thermal nonequilibrium. Thus, the governing equations of proposed diffuse interface model include: the conservation equation for mass of every component and the conservation equations for momentum and total energy of the mixture, and the evolution equations for volume fraction of every component and for pressure of the mixture. The important character of the present model is that the mass transfer due to chemical reaction and the heat exchange due to thermal nonequilibrium are involved. In this model, pressure is solved directly from the governing equations instead of computed next from the obtained conservative variables. The present model may apply to arbitrary expression of equation of state and allow for any number of inert media. The partially differential governing equations of the diffuse interface model are numerically solved by a temporal-spatial second-order Godunov-type finite volume scheme with wave propagation algorithm. From numerical examples, the proposed diffuse interface model can eliminate the unphysical oscillations near the material interfaces, and obtain the agreeable results with the physical mechanism.
- diffuse interface model /
- solid explosive detonation /
- multimaterial flow /
- thermal nonequilibrium /
- multiphase flow model

HTML全文

爆轰是复杂的流体力学与剧烈的化学动力学耦合的一种物理现象。固体炸药的爆轰能够提供强大的能量来驱动和压缩材料，在国防科技和国民经济中获得广泛应用，因此，固体炸药爆轰与惰性介质相互作用问题是工程应用中十分重要的问题。固体炸药爆轰与惰性介质相互作用问题本质上是一种可压缩多物质流动问题，除了可压缩单物质流动中出现的激波、稀疏波、接触间断外，一种新的间断，即物质界面出现了，它隔开热力学性质或状态方程不一样的两种物质，这种物质界面实质上也是一种接触间断，它满足物质界面压力相等与速度相等的条件。对可压缩多物质流动，传统的数值模拟方法多是拉格朗日方法^[1-3]，它不能有效地处理物质大变形运动，为了较好地处理物质大变形运动，欧拉方法成为一个很好的选择^[3-4]。根据对处理物质界面的不同方式，欧拉方法大致可以分为4类：阵面追踪（front tracking）法^[5-6]、流体体积（volume of fluid）法^[7-8]、水平集合（level set）法^[9-10]以及扩散界面（diffuse interface）法^[11-17]等，其中，扩散界面法由于具有逻辑结构简单、物理量守恒性良好、对界面形状没有几何拓扑限制、对多物理问题适应性较强的特点，近年来获得了越来越多的关注。

扩散界面法将离散网格视为包含有多种组分物质的混合网格，物质界面被认为是具有一定厚度的虚拟混合区，即把无厚度的物质界面当作有一定厚度的扩散界面，扩散界面内部采用“虚拟状态方程”或混合规则来描述。根据对多物质混合状态的不同处理方式，扩散界面模型可以分为两种^[18]：基于单物质流动欧拉方程组的扩展欧拉模型（augmented Euler model）以及基于多相流动Baer-Nunziato方程组的多相流模型（multiphase flow model）。

在扩展欧拉模型中，混合规则通常采用基于力学平衡和热学平衡的假设。事实上，热学平衡仅适用于均匀混合状况，而包含有物质界面的混合网格内各组分物质通常不是均匀混合的，由此将热学平衡的假设用于计算包含界面接触间断的混合物质的热力学性质时，往往使得物质界面附近的压力和速度等物理量出现非物理振荡现象^[19]。为了消除物质界面附近的非物理振荡现象，不得不采取另外的技术方案进行修正处理，如附加状态方程参数演化方程^[20]、总能量调整^[21]、守恒与原始变量转换^[22]等。

不同于扩展欧拉模型，多相流模型^[23]中的混合物质被认为由处于热学、力学、化学非平衡的多种组分物质组成，每种组分物质具有各自的物理状态并按照各自的动力学规律分别进行演化，演化方程含有用于表达由组分物质非平衡引起的质量、动量和能量相互转化的各种源项，同时采用组分物质体积分数方程来描述物质界面的运动过程。由于考虑了每种组分物质各自的变化规律，能够保证混合物质的热力学性质自动满足一致性，从物理建模上消除了物质界面附近的非物理振荡，对诸如爆轰这种带有化学非平衡的可压缩流动，消除物质界面附近的非物理压力振荡非常重要，因为爆轰化学反应的激发与压力等量密切相关，如果压力出现非物理振荡会激发虚假的化学反应，进一步会引起错误的爆轰特性。针对具体物理过程，热学、力学、化学非平衡这些状态可以全部满足也可以部分满足，在全部满足非平衡条件下即为著名的Baer-Nunziato模型，在部分满足非平衡条件下形成各种简化多相流模型，如在力学平衡、热学非平衡条件下有Kapila模型^[24]、Allaire模型^[25]、Massoni模型^[26]、Murrone模型^[27]、Grove模型^[28]等多种典型模型，这些不同模型的差异主要在于以不同方式处理组分物质混合状态，使得组分物质体积分数方程有不同的表达形式。多相流模型由于考虑了更精确的物理机制，更加符合物理意义，近年来在可压缩多物质流动数值模拟领域获得了极大的重视，本文的工作正是基于多相流模型。

在固体炸药爆轰与惰性介质的相互作用过程中，爆轰化学反应过程通常简化为固相反应物转化成气相生成物，这样组分混合物质通常包括固相反应物、气相生成物、惰性介质这3种成分，由于这3种组分物质的材料物态性质和热力学性质差异极大，因此它们之间的相互作用过程被认为满足力学平衡和热学非平衡状态^[13,24,29]，即在流场控制体中每种组分物质拥有相同的压力和速度、以及不同的温度和内能。基于多相流思想，固体炸药爆轰与惰性介质相互作用过程首先由各种组分物质的质量守恒方程、混合物质的动量与总能量守恒方程描述，组分物质的质量守恒方程还需考虑由化学反应引起的质量转化的影响，然后需要补充每种组分物质的体积分数的控制方程。由混合物质能量守恒方程可以分解获得每种组分物质的内能变化方程，再利用每种组分物质压力相等的条件并结合组分物质的质量守恒方程，可以推导出每种组分物质的体积分数控制方程。由于涉及到热学非平衡状态，组分物质之间存在热量交换，分解获得的组分物质内能变化方程还需考虑热量交换的影响。这样，推导出的组分物质体积分数控制方程同时包含了化学反应和热学非平衡的影响。并且，混合物质的压力演化方程也被纳入到模型方程组中，这样压力通过直接离散求解演化方程获得而不是由流动守恒变量计算获得，这种方案可以增强消除物质界面非物理振荡的效果，也可以避免由状态方程非线性形式引起的压力迭代求解^[28]或者由组分体积方程非守恒型引起的压力松弛求解^[13]，还可以避免计算压力对守恒变量的导数而简化高阶精度格式的运算过程。最终，获得的扩散界面模型方程组包括组分物质的质量守恒方程、混合物质的动量及总能量守恒方程、组分物质的体积分数演化方程以及混合物质的压力演化方程。获得的扩散界面模型的最主要特点是考虑了化学反应以及热学非平衡的影响，因此具有良好的热力学一致性，同时，该扩散界面模型能够适用于任意表达形式的状态方程以及任意数目的多种惰性介质。

对所获得的扩散界面模型方程组采用一个具有波传播性质的时空二阶精度的有限体积法进行数值求解，典型算例结果显示，数值模拟图像与物理规律符合，物质界面附近不会出现物理量的非物理振荡现象。

1. 扩散界面模型的推导

首先给出物理量的定义。考虑一个控制体包含有炸药固相反应物、炸药气相生成物以及惰性物质，固相反应物的物理量用下标s表示、气相生成物的物理量用下标g表示，惰性物质设有K种物质，每种物质的物理量用下标k表示；其中 $\rho$ 表示密度， $e$ 表示内能， $p$ 表示压力， ${{u}}$ 表示速度矢量， $\alpha$ 表示体积分数， $v$ 表示比容， $Q$ 表示单位质量的固体炸药由固相反应物转化为气相生成物时所释放的热量。在力学平衡条件下，混合物质物理量与组分物质物理量的关系式为：

$\rho = {\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}} + \sum\limits_{k = 1}^K {{\rho _k}{\alpha _k}}$

(1)

${\kern 1pt} e = \frac{{{\rho _{\rm{s}}}{\alpha _{\rm{s}}}}}{\rho }{e_{\rm{s}}} + \frac{{{\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{\rho }{e_{\rm{g}}} - \frac{{{\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{\rho }Q + \sum\limits_{k = 1}^K {\frac{{{\rho _k}{\alpha _k}}}{\rho }{e_k}}$

(2)

${\kern 1pt} v = \frac{{{\alpha _{\rm{s}}}}}{\rho } + \frac{{{\alpha _{\rm{g}}}}}{\rho } + \sum\limits_{k = 1}^K {\frac{{{\alpha _k}}}{\rho }}$

(3)

$p = {p_{\rm{s}}} = {p_{\rm{g}}} = {p_k}\quad\quad\quad k = 1,2, \cdots ,K$

(4)

${{u}} = {\kern 1pt} {{{u}}_{\rm{s}}} = {\kern 1pt} {{{u}}_{\rm{g}}} = {\kern 1pt} {{{u}}_k}\quad\quad\quad k = 1,2, \cdots ,K$

(5)

${\kern 1pt} {\alpha _{\rm{s}}} + {\alpha _{\rm{g}}} + \sum\limits_{k = 1}^K {{\alpha _k}} = 1$

(6)

不考虑各种耗散因素及外力做功情况，则在力学平衡状态条件下固体炸药爆轰与惰性介质的可压缩流动方程组可以表达为：

$\frac{{\partial {\rho _{\rm{s}}}{\alpha _{\rm{s}}}}}{{\partial {\kern 1pt} t}} + \nabla \cdot ({\rho _{\rm{s}}}{\alpha _{\rm{s}}}{\kern 1pt} {{u}}) = - ({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}})\,r$

(7)

$\frac{{\partial {\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{{\partial {\kern 1pt} t}} + \nabla \cdot ({\rho _{\rm{g}}}{\alpha _{\rm{g}}}{\kern 1pt} {{u}}) = ({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}})\,r$

(8)

$\frac{{\partial {\rho _k}{\alpha _k}}}{{\partial {\kern 1pt} t}} + \nabla \cdot ({\rho _k}{\alpha _k}{\kern 1pt} {{u}})= 0\quad\quad\quad k = 1,2, \cdots ,K$

(9)

$\frac{{\partial \rho {\kern 1pt} {{u}}}}{{\partial {\kern 1pt} t}} + \nabla \cdot (\rho {\kern 1pt} {{u}} {\kern 1pt} {{u}} + p{{I}}) = {\bf 0}$

(10)

$\frac{{\partial \rho {\kern 1pt} E}}{{\partial {\kern 1pt} t}} + \nabla \cdot (\rho {\kern 1pt} E + p){{u}} = 0$

(11)

式中： $E = e + \dfrac{1}{2}{{{u}}^2}$ 为总能量， $r$ 为固体炸药的化学反应率， ${{I}}$ 为单位张量。

式(7)为固体炸药固相反应物的质量守恒方程，式(8)为固体炸药气相生成物的质量守恒方程，式(9)为惰性介质的质量守恒方程，式(10)为混合物质的动量守恒方程，式(11)为混合物质的总能量守恒方程。

在给定每种组分物质的状态方程条件下，并确定了每种组分物质的体积分数的控制方程之后，上述流动方程组成为封闭系统进而可以获得定解。组分物质体积分数的控制方程可由混合物质的能量守恒方程并结合混合网格内每种组分物质压力相等的条件导出，混合物质的能量守恒方程可以分解为每种组分物质的内能变化方程。

上述流动方程组式(7)～(11)可以经过变换得到混合物质的内能守恒方程：

$\frac{{{\rm d}e}}{{{\rm d}t}} + p\frac{{{\rm d}v}}{{{\rm d}t}} = 0$

(12)

将式(1)～(3)代入式(12)，有：

$\frac{{\rm d}}{{{\rm d}t}}\left( {\frac{{{\rho _{\rm{s}}}{\alpha _{\rm{s}}}}}{\rho }{e_{\rm{s}}} + \frac{{{\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{\rho }{e_{\rm{g}}} - \frac{{{\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{\rho }Q + \sum\limits_{k = 1}^K {\frac{{{\rho _k}{\alpha _k}}}{\rho }{e_k}} } \right) + p\frac{\rm d}{{{\rm d}t}}\left( {\frac{{{\alpha _{\rm{s}}}}}{\rho } + \frac{{{\alpha _{\rm{g}}}}}{\rho } + \sum\limits_{k = 1}^K {\frac{{{\alpha _k}}}{\rho }} } \right) = 0$

(13)

式(13)也可化为：

$\left( {\frac{\rm d}{{{\rm d}t}}\frac{{{\rho _{\rm{s}}}{\alpha _{\rm{s}}}}}{\rho }{e_{\rm{s}}} + p\frac{\rm d}{{{\rm d}t}}\frac{{{\alpha _{\rm{s}}}}}{\rho }} \right) + \left( {\frac{\rm d}{{{\rm d}t}}\frac{{{\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{\rho }{e_{\rm{g}}} - \frac{\rm d}{{{\rm d}t}}\frac{{{\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{\rho }Q + p\frac{\rm d}{{{\rm d}t}}\frac{{{\alpha _{\rm{g}}}}}{\rho }} \right) + \sum\limits_{k = 1}^K {\left( {\frac{\rm d}{{{\rm d}t}}\frac{{{\rho _k}{\alpha _k}}}{\rho }{e_k} + p\frac{\rm d}{{{\rm d}t}}\frac{{{\alpha _k}}}{\rho }} \right)} = 0$

(14)

式(14)中3个括号中分别表示炸药固相反应物、炸药气相生成物、惰性介质的内能变化，再考虑到这多种物质由于处于热学非平衡状态而存在温度差异，物质之间会产生热量交换，则可以将式(14)分裂为各种组分物质的内能变化方程：

$\frac{\rm d}{{{\rm d}t}}\frac{{{\rho _{\rm{s}}}{\alpha _{\rm{s}}}}}{\rho }{e_{\rm{s}}} + p\frac{\rm d}{{{\rm d}t}}\frac{{{\alpha _{\rm{s}}}}}{\rho } = {q_{\rm{s}}}$

(15)

$\frac{\rm d}{{{\rm d}t}}\frac{{{\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{\rho }{e_{\rm{g}}} - Q\frac{\rm d}{{{\rm d}t}}\frac{{{\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{\rho } + p\frac{\rm d}{{{\rm d}t}}\frac{{{\alpha _{\rm{g}}}}}{\rho } = {q_{\rm{g}}}$

(16)

$\frac{\rm d}{{{\rm d}t}}\frac{{{\rho _k}{\alpha _k}}}{\rho }{e_k} + p\frac{\rm d}{{{\rm d}t}}\frac{{{\alpha _k}}}{\rho } = {q_k}\quad\quad\quad k = 1,2, \cdots ,K$

(17)

式中： ${q_{\rm{s}}}$ 、 ${q_{\rm{g}}}$ 、 ${q_k}$ ( ${\kern 1pt} k = 1,2, \cdots ,K$ )表示单位质量的物质由于温度差产生的热流率，并有约束条件 ${q_{\rm{s}}} + {q_{\rm{g}}} + \displaystyle\sum\limits_{k = 1}^K {{q_k}} = 0$ 。

针对炸药固相反应物流动过程，式(15)借助混合物质质量守恒方程可以化成如下形式：

$\frac{{{\rm d}{\rho _{\rm{s}}}{e_{\rm{s}}}{\alpha _{\rm{s}}}}}{{{\rm d}t}} + {\rho _{\rm{s}}}{e_{\rm{s}}}{\alpha _{\rm{s}}}\nabla \cdot {\kern 1pt} {\kern 1pt} {{u}} + {\alpha _{\rm{s}}}p\nabla \cdot {{u}} + p\frac{{{\rm d}{\alpha _{\rm{s}}}}}{{{\rm d}t}} = \rho {\kern 1pt} {\kern 1pt} {q_{\rm{s}}}$

(18)

由式(7)、(18)可得固相反应物的内能变化方程：

$\frac{{{\rm d}{e_{\rm{s}}}}}{{{\rm d}t}} = \frac{{({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}})\,{e_{\rm{s}}}r + \rho {\kern 1pt} {\kern 1pt} {q_{\rm{s}}}}}{{{\rho _{\rm{s}}}{\alpha _{\rm{s}}}}} - \frac{p}{{{\rho _{\rm{s}}}}}\nabla \cdot {\kern 1pt} {{u}} - \frac{p}{{{\rho _{\rm{s}}}{\alpha _{\rm{s}}}}}\frac{{{\rm d}{\alpha _{\rm{s}}}}}{{{\rm d}t}}$

(19)

为简化推导过程可不失一般性地认为，每种物质具有 ${p_i}({\rho _i},{e_i})$ （ $i = s,g,1,2, \cdots ,K$ ）表达形式的状态方程，则在混合网格内各种组分物质压力相等： $p = {p_{\rm{s}}}({\rho _{\rm{s}}},{e_{\rm{s}}}) = {p_{\rm{g}}}({\rho _{\rm{g}}},{e_{\rm{g}}}) = {p_k}({\rho _k},{e_k})$ ( $k = 1,2, \cdots ,K$ )的条件下，可以得到固相反应物的压力控制方程：

$\frac{{{\rm d}p}}{{{\rm d}t}} = \frac{{\partial {p_{\rm{s}}}({\rho _{\rm{s}}},{e_{\rm{s}}})}}{{\partial {\rho _{\rm{s}}}}}\frac{{{\rm d}{\rho _{\rm{s}}}}}{{{\rm d}t}} + \frac{{\partial {p_{\rm{s}}}({\rho _{\rm{s}}},{e_{\rm{s}}})}}{{\partial {e_{\rm{s}}}}}\frac{{{\rm d}{e_{\rm{s}}}}}{{{\rm d}t}}$

(20)

为求出固相反应物的密度变化，由式(7)可得：

$\frac{{{\rm d}{\rho _{\rm{s}}}}}{{{\rm d}t}} = - \frac{{{\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{{{\alpha _{\rm{s}}}}}r - {\rho _{\rm{s}}}\nabla \cdot {\kern 1pt} {{u}} - \frac{{{\rho _{\rm{s}}}}}{{{\alpha _{\rm{s}}}}}\frac{{{\rm d}{\alpha _{\rm{s}}}}}{{{\rm d}{\kern 1pt} t}}$

(21)

将式(19)、(21)代入式(20)中，可得：

$\frac{{{\rm d}p}}{{{\rm d}t}} = - \frac{{\partial {p_{\rm{s}}}({\rho _{\rm{s}}},{e_{\rm{s}}})}}{{\partial {\rho _{\rm{s}}}}}\left( {\frac{{({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}})\,r}}{{{\alpha _{\rm{s}}}}} + {\rho _{\rm{s}}}\nabla \cdot {\kern 1pt} {{u}} + \frac{{{\rho _{\rm{s}}}}}{{{\alpha _{\rm{s}}}}}\frac{{{\rm d}{\alpha _{\rm{s}}}}}{{{\rm d}{\kern 1pt} t}}} \right) + \frac{{\partial {p_{\rm{s}}}({\rho _{\rm{s}}},{e_{\rm{s}}})}}{{\partial {e_{\rm{s}}}}}\left( {\frac{{({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}})\,{e_{\rm{s}}}r + \rho {q_{\rm{s}}}}}{{{\rho _{\rm{s}}}{\alpha _{\rm{s}}}}} - \frac{p}{{{\rho _{\rm{s}}}}}\nabla \cdot {\kern 1pt} {{u}} - \frac{p}{{{\rho _{\rm{s}}}{\alpha _{\rm{s}}}}}\frac{{{\rm d}{\alpha _{\rm{s}}}}}{{{\rm d}t}}} \right)$

(22)

式(22)也可以转化成如下固相反应物的体积分数变化方程：

$\frac{{{\rm d}{\alpha _{\rm{s}}}}}{{{\rm d}t}} = \frac{{ - ({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}})\,r\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {\rho _{\rm{s}}}}} + \dfrac{{({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}})\,{e_{\rm{s}}}r + \rho {q_{\rm{s}}}}}{{{\rho _{\rm{s}}}}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}}}}{{{\rho _{\rm{s}}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {\rho _{\rm{s}}}}} + \dfrac{p}{{{\rho _{\rm{s}}}}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}}}} - {\alpha _{\rm{s}}}\nabla \cdot {\kern 1pt} {{u}} - \dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {\rho _{\rm{s}}}}} + \dfrac{p}{{{\rho _{\rm{s}}}}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}}}}\frac{{{\rm d}p}}{{{\rm d}t}}$

(23)

同理可得到气相生成物、惰性物质的体积分数变化方程，这样各组分物质的体积分数变化方程有：

$\frac{{{\rm d}{\alpha _{\rm{s}}}}}{{{\rm d}t}} = - \frac{{\partial {p_{\rm{s}}}}}{{\partial {\rho _{\rm{s}}}}}\frac{{({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}})r}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}}{\kern 1pt} {\kern 1pt} {\kern 1pt} + \frac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}}\frac{{({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}}){e_{\rm{s}}}r + \rho {q_{\rm{s}}}}}{{\rho _{\rm{s}}^2c_{\rm{s}}^2}} - {\alpha _{\rm{s}}}\nabla \cdot {\kern 1pt} {{u}} - \frac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}}\frac{{{\rm d}p}}{{{\rm d}t}}$

(24)

$\frac{{{\rm d}{\alpha _{\rm{g}}}}}{{{\rm d}t}} = \frac{{\partial {p_{\rm{g}}}}}{{\partial {\rho _{\rm{g}}}}}\frac{{({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}})r}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}}{\kern 1pt} {\kern 1pt} {\kern 1pt} + \frac{{\partial {p_{\rm{g}}}}}{{\partial {e_{\rm{g}}}}}\frac{{({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}})(Q - {e_{\rm{g}}})r + \rho {q_{\rm{g}}}}}{{\rho _{\rm{g}}^2c_{\rm{g}}^2}} - {\alpha _{\rm{g}}}\nabla \cdot {\kern 1pt} {{u}} - \frac{{{\alpha _{\rm{g}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}}\frac{{{\rm d}p}}{{{\rm d}t}}$

(25)

${\kern 1pt} \frac{{{\rm d}{\alpha _k}}}{{{\rm d}t}} = \frac{{\partial {p_k}}}{{\partial {\rho _k}}}\frac{{\rho {q_k}}}{{\rho _k^2c_k^2}} - {\alpha _k}\nabla \cdot {\kern 1pt} {{u}} - \frac{{{\alpha _k}}}{{{\rho _k}c_k^2}}\frac{{{\rm d}p}}{{{\rm d}t}}\quad \quad \quad k = 1, \cdots ,K,$

(26)

式中： $c_{\rm{s}}^2 = \dfrac{{\partial {p_{\rm{s}}}}}{{\partial {\rho _{\rm{s}}}}} + \dfrac{p}{{\rho _{\rm{s}}^2}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}}$ ， $c_{\rm{g}}^2 = \dfrac{{\partial {p_{\rm{g}}}}}{{\partial {\rho _{\rm{g}}}}} + \dfrac{p}{{\rho _{\rm{g}}^2}}\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {e_{\rm{g}}}}}$ ， $c_k^2 = \dfrac{{\partial {p_k}}}{{\partial {\rho _k}}} + \dfrac{p}{{\rho _k^2}}\dfrac{{\partial {p_k}}}{{\partial {e_k}}}$ 。

在固相反应物、气相生成物、惰性物质满足体积分数的总和为一的条件下，即有恒等式 ${\alpha _{\rm{s}}} + {\alpha _{\rm{g}}} + \displaystyle\sum\limits_{k = 1}^K {{\alpha _k}} = 1$ ，则由式(24)～(26)可以获得混合物质的压力演化方程：

$\begin{split} \frac{{{\rm d}p}}{{{\rm d}t}} + \frac{1}{{\dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} + \dfrac{{{\alpha _{\rm{g}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}} + \displaystyle\sum\limits_{k = 1}^K {\dfrac{{{\alpha _k}}}{{{\rho _k}c_k^2}}} }}\nabla \cdot {\kern 1pt} {{u}}\, =& - \left( {\frac{{\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {\rho _{\rm{s}}}}} - \dfrac{{{e_{\rm{s}}}}}{{{\rho _{\rm{s}}}}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} - \dfrac{{\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {\rho _{\rm{g}}}}} - \dfrac{{{e_{\rm{g}}} - Q}}{{{\rho _{\rm{g}}}}}\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {e_{\rm{g}}}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}}} \right)\dfrac{{{\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{{\dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} + \dfrac{{{\alpha _{\rm{g}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}} + \displaystyle\sum\limits_{k = 1}^K {\dfrac{{{\alpha _k}}}{{{\rho _k}c_k^2}}} }}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} r + \\&\frac{{\dfrac{{\rho {q_{\rm{s}}}}}{{\rho _{\rm{s}}^2c_{\rm{s}}^2}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}} + \dfrac{{\rho {q_{\rm{g}}}}}{{\rho _{\rm{g}}^2c_{\rm{g}}^2}}\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {e_{\rm{g}}}}} + \displaystyle\sum\limits_{k = 1}^K {\dfrac{{\rho {q_k}}}{{\rho _k^2c_k^2}}\dfrac{{\partial {p_k}}}{{\partial {e_k}}}} }}{{\dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} + \dfrac{{{\alpha _{\rm{g}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}} + \displaystyle\sum\limits_{k = 1}^K {\dfrac{{{\alpha _k}}}{{{\rho _k}c_k^2}}} }} \end{split}$

(27)

将式(27)代入式(24)，可以获得固相反应物的体积分数演化方程：

$\begin{split} &\dfrac{{{\rm d}{\alpha _{\rm{s}}}}}{{{\rm d}t}} \!+\! {\alpha _{\rm{s}}}\left( {1 - \dfrac{1}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}}\frac{1}{{\dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} + \dfrac{{{\alpha _{\rm{g}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}}\! +\! \displaystyle\sum\limits_{k = 1}^K {\dfrac{{{\alpha _k}}}{{{\rho _k}c_k^2}}} }}} \right)\nabla \cdot {\kern 1pt} {{u}} = \dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}}\left( {\frac{{\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {\rho _{\rm{s}}}}} - \dfrac{{{e_{\rm{s}}}}}{{{\rho _{\rm{s}}}}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} - \frac{{\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {\rho _{\rm{g}}}}} - \dfrac{{{e_{\rm{g}}} - Q}}{{{\rho _{\rm{g}}}}}\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {e_{\rm{g}}}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}}} \right)\frac{{{\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{{\dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} \!+\! \dfrac{{{\alpha _{\rm{g}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}} \!+ \!\displaystyle\sum\limits_{k = 1}^K {\dfrac{{{\alpha _k}}}{{{\rho _k}c_k^2}}} }}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} r -\\ & \quad\quad\quad\quad\quad\quad\quad \frac{{\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {\rho _{\rm{s}}}}} - \dfrac{{{e_{\rm{s}}}}}{{{\rho _{\rm{s}}}}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}}({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}})r- \frac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}}\frac{{\dfrac{{\rho {q_{\rm{s}}}}}{{\rho _{\rm{s}}^2c_{\rm{s}}^2}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}} + \dfrac{{\rho {q_{\rm{g}}}}}{{\rho _{\rm{g}}^2c_{\rm{g}}^2}}\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {e_{\rm{g}}}}} + \displaystyle\sum\limits_{k = 1}^K {\dfrac{{\rho {q_k}}}{{\rho _k^2c_k^2}}\dfrac{{\partial {p_k}}}{{\partial {e_k}}}} }}{{\dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} + \dfrac{{{\alpha _{\rm{g}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}} + \displaystyle\sum\limits_{k = 1}^K {\dfrac{{{\alpha _k}}}{{{\rho _k}c_k^2}}} }} + \dfrac{{\rho {q_{\rm{s}}}}}{{\rho _{\rm{s}}^2c_{\rm{s}}^2}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}}\\[-26pt] \end{split}$

(28)

同理可以得到气相生成物与惰性介质的体积分数控制方程，这样各组分物质的体积分数控制方程均可获得，同时，混合物质的压力演化方程也成为一个流动控制方程，这样固体炸药爆轰与惰性介质相互作用的扩散界面模型被推导出：

$\begin{aligned} {\frac{{\partial {\rho _{\rm{s}}}{\alpha _{\rm{s}}}}}{{\partial t}} + \nabla \cdot {\rho _{\rm{s}}}{\alpha _{\rm{s}}} {{u}} = - ({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}})\,r} \end{aligned}$

(29)

$\begin{aligned} {\frac{{\partial {\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{{\partial t}} + \nabla \cdot {\rho _{\rm{g}}}{\alpha _{\rm{g}}} {{u}} = ({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}})\,r} \end{aligned}$

(30)

$\begin{aligned} {\frac{{\partial {\rho _k}{\alpha _k}}}{{\partial t}} + \nabla \cdot {\rho _k}{\alpha _k} {{u}} = 0\quad\quad\quad k = 1,2, \cdots ,K} \end{aligned}$

(31)

$\begin{aligned} {\frac{{\partial \rho {{u}}}}{{\partial t}} + \nabla \cdot (\rho {{u}}{{u}} + p{{I}}) = {\bf{0}}} \end{aligned}$

(32)

$\begin{aligned} {\frac{{\partial \rho E}}{{\partial t}} + \nabla \cdot (\rho E + p){{u}} = 0} \end{aligned}$

(33)

$\begin{split} &\frac{{\partial {\alpha _{\rm{s}}}}}{{\partial t}} + {{u}} \cdot \nabla {\alpha _{\rm{s}}} + {\alpha _{\rm{s}}}\left( {1 - \frac{1}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}}\frac{1}{{\dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} + \dfrac{{{\alpha _{\rm{g}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}} + \displaystyle\sum\limits_{k = 1}^K {\frac{{{\alpha _k}}}{{{\rho _k}c_k^2}}} }}} \right)\nabla \cdot {{u}} = \frac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}}\left( {\frac{{\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {\rho _{\rm{s}}}}} - \dfrac{{{e_{\rm{s}}}}}{{{\rho _{\rm{s}}}}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} - \frac{{\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {\rho _{\rm{g}}}}} - \dfrac{{{e_{\rm{g}}} - Q}}{{{\rho _{\rm{g}}}}}\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {e_{\rm{g}}}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}}} \right)\times\\ &\quad \frac{{{\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{{\dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} \!\!+\!\! \dfrac{{{\alpha _{\rm{g}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}}\!\!+\!\! \displaystyle\sum\limits_{k = 1}^K {\dfrac{{{\alpha _k}}}{{{\rho _k}c_k^2}}} }} r -\frac{{\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {\rho _{\rm{s}}}}} - \dfrac{{{e_{\rm{s}}}}}{{{\rho _{\rm{s}}}}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}}({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}})r - \dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}}\dfrac{{\dfrac{{\rho {q_{\rm{s}}}}}{{\rho _{\rm{s}}^2c_{\rm{s}}^2}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}} + \dfrac{{\rho {q_{\rm{g}}}}}{{\rho _{\rm{g}}^2c_{\rm{g}}^2}}\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {e_{\rm{g}}}}} + \displaystyle\sum\limits_{k = 1}^K {\dfrac{{\rho {q_k}}}{{\rho _k^2c_k^2}}\dfrac{{\partial {p_k}}}{{\partial {e_k}}}} }}{{\dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} + \dfrac{{{\alpha _{\rm{g}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}} + \displaystyle\sum\limits_{k = 1}^K {\dfrac{{{\alpha _k}}}{{{\rho _k}c_k^2}}} }} + \frac{{\rho {q_{\rm{s}}}}}{{\rho _{\rm{s}}^2c_{\rm{s}}^2}}\frac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}} \end{split}$

(34)

$\begin{split} &\frac{{\partial {\alpha _k}}}{{\partial t}} + {{u}} \cdot \nabla {\alpha _k} + {\alpha _k}\left( {1 - \frac{1}{{{\rho _k}c_k^2}}\frac{1}{{\dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} + \dfrac{{{\alpha _{\rm{g}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}} + \displaystyle\sum\limits_{k = 1}^K {\dfrac{{{\alpha _k}}}{{{\rho _k}c_k^2}}} }}} \right)\nabla \cdot {{u}} = \frac{{{\alpha _k}}}{{{\rho _k}c_k^2}}\left( {\frac{{\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {\rho _{\rm{s}}}}} - \dfrac{{{e_{\rm{s}}}}}{{{\rho _{\rm{s}}}}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} - \frac{{\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {\rho _{\rm{g}}}}} - \dfrac{{{e_{\rm{g}}} - Q}}{{{\rho _{\rm{g}}}}}\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {e_{\rm{g}}}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}}} \right)\times\\ &\quad \quad\frac{{{\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{{\dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} + \dfrac{{{\alpha _{\rm{g}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}} + \displaystyle\sum\limits_{k = 1}^K {\dfrac{{{\alpha _k}}}{{{\rho _k}c_k^2}}} }} r - \frac{{{\alpha _k}}}{{{\rho _k}c_k^2}}\frac{{\dfrac{{\rho {q_{\rm{s}}}}}{{\rho _{\rm{s}}^2c_{\rm{s}}^2}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}} + \dfrac{{\rho {q_{\rm{g}}}}}{{\rho _{\rm{g}}^2c_{\rm{g}}^2}}\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {e_{\rm{g}}}}} + \displaystyle\sum\limits_{k = 1}^K {\dfrac{{\rho {q_k}}}{{\rho _k^2c_k^2}}\dfrac{{\partial {p_k}}}{{\partial {e_k}}}} }}{{\dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} + \dfrac{{{\alpha _{\rm{g}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}} + \displaystyle\sum\limits_{k = 1}^K {\dfrac{{{\alpha _k}}}{{{\rho _k}c_k^2}}} }} + \frac{{\rho {q_k}}}{{\rho _k^2c_k^2}}\frac{{\partial {p_k}}}{{\partial {e_k}}}\quad\quad k = 1,2, \cdots ,K \end{split}$

(35)

$\begin{split} &\frac{{\partial p}}{{\partial t}} + {{u}} \cdot \nabla p + \frac{1}{{\dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} + \dfrac{{{\alpha _{\rm{g}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}} + \displaystyle\sum\limits_{k = 1}^K {\dfrac{{{\alpha _k}}}{{{\rho _k}c_k^2}}} }}\nabla \cdot {{u}} = - \left( {\frac{{\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {\rho _{\rm{s}}}}} - \dfrac{{{e_{\rm{s}}}}}{{{\rho _{\rm{s}}}}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} - \frac{{\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {\rho _{\rm{g}}}}} - \dfrac{{{e_{\rm{g}}} - Q}}{{{\rho _{\rm{g}}}}}\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {e_{\rm{g}}}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}}} \right)\frac{{{\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{{\dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} + \dfrac{{{\alpha _{\rm{g}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}} + \displaystyle\sum\limits_{k = 1}^K {\dfrac{{{\alpha _k}}}{{{\rho _k}c_k^2}}} }} r +\\ &\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \frac{{\dfrac{{\rho {q_{\rm{s}}}}}{{\rho _{\rm{s}}^2c_{\rm{s}}^2}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}} + \dfrac{{\rho {q_{\rm{g}}}}}{{\rho _{\rm{g}}^2c_{\rm{g}}^2}}\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {e_{\rm{g}}}}} + \displaystyle\sum\limits_{k = 1}^K {\dfrac{{\rho {q_k}}}{{\rho _k^2c_k^2}}\dfrac{{\partial {p_k}}}{{\partial {e_k}}}} }}{{\dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} + \dfrac{{{\alpha _{\rm{g}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}} + \displaystyle\sum\limits_{k = 1}^K {\dfrac{{{\alpha _k}}}{{{\rho _k}c_k^2}}} }} \end{split}$

(36)

2. 数值方法

为了便于数值求解，获得的界面扩散模型方程组(29)～(36)在二维直角坐标系下可以写成如下矢量-矩阵形式的偏微分方程组：

$\frac{{\partial {\kern 1pt} {{q}}}}{{\partial {\kern 1pt} t}} + {{A}}({{q}})\frac{{\partial {\kern 1pt} {{q}}}}{{\partial {\kern 1pt} x}} + {{B}}({{q}})\frac{{\partial {\kern 1pt} {{q}}}}{{\partial {\kern 1pt} y}} = {{s}}({{q}})$

(37)

这里状态变量 ${{q}} = {[{\rho _{\rm{s}}}{\alpha _{\rm{s}}},\;{\kern 1pt} \,{\rho _{\rm{g}}}{\alpha _{\rm{g}}},\; {\rho _k}{\alpha _k},\;\rho {\kern 1pt} u,\;\rho {\kern 1pt} v,\;\rho E,{\kern 1pt} \;{\alpha _{\rm{s}}},\;{\alpha _k},{\kern 1pt} \,p]^{\rm T}}$ （不失一般性，惰性介质只写一种物质），矩阵 ${{A}}({{q}})$ 与矩阵 ${{B}}({{q}})$ 是具有类似表达形式的Jacobi矩阵，其中

${{A}}({{q}}) = \left[ {\begin{array}{*{20}{c}} {\left(1 - \dfrac{{{\rho _{\rm{s}}}{\alpha _{\rm{s}}}}}{\rho }\right)u}&{ - \dfrac{{{\rho _{\rm{s}}}{\alpha _{\rm{s}}}}}{\rho }u}&{ - \dfrac{{{\rho _{\rm{s}}}{\alpha _{\rm{s}}}}}{\rho }u}&{\dfrac{{{\rho _{\rm{s}}}{\alpha _{\rm{s}}}}}{\rho }}&0&0&0&0&0 \\ { - \dfrac{{{\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{\rho }u}&{\left(1 - \dfrac{{{\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{\rho }\right)u}&{ - \dfrac{{{\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{\rho }u}&{\dfrac{{{\rho _{\rm{g}}}{\alpha _{\rm{g}}}}}{\rho }}&0&0&0&0&0 \\ { - \dfrac{{{\rho _k}{\alpha _k}}}{\rho }u}&{ - \dfrac{{{\rho _k}{\alpha _k}}}{\rho }u}&{\left(1 - \dfrac{{{\rho _k}{\alpha _k}}}{\rho }\right)u}&{\dfrac{{{\rho _k}{\alpha _k}}}{\rho }}&0&0&0&0&0 \\ { - {u^2}}&{ - {u^2}}&{ - {u^2}}&{2u}&0&0&0&0&1 \\ { - uv}&{ - uv}&{ - uv}&v&u&0&0&0&0 \\ { - Hu}&{ - Hu}&{ - Hu}&H&0&u&0&0&u \\ { - \dfrac{{{\alpha _{\rm{s}}}}}{\rho }\left(1 - \dfrac{{\rho {c^2}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}}\right)u}&{ - \dfrac{{{\alpha _{\rm{s}}}}}{\rho }\left(1 - \dfrac{{\rho {c^2}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}}\right)u}&{ - \dfrac{{{\alpha _{\rm{s}}}}}{\rho }\left(1 - \dfrac{{\rho {c^2}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}}\right)u}&{\dfrac{{{\alpha _{\rm{s}}}}}{\rho }\left(1 - \dfrac{{\rho {c^2}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}}\right)}&0&0&u&0&0 \\ { - \dfrac{{{\alpha _k}}}{\rho }\left(1 - \dfrac{{\rho {c^2}}}{{{\rho _k}c_k^2}}\right)u}&{ - \dfrac{{{\alpha _k}}}{\rho }\left(1 - \dfrac{{\rho {c^2}}}{{{\rho _k}c_k^2}}\right)u}&{ - \dfrac{{{\alpha _k}}}{\rho }\left(1 - \dfrac{{\rho {c^2}}}{{{\rho _k}c_k^2}}\right)u}&{\dfrac{{{\alpha _k}}}{\rho }\left(1 - \dfrac{{\rho {c^2}}}{{{\rho _k}c_k^2}}\right)}&0&0&0&u&0 \\ { - {c^2}u}&{ - {c^2}u}&{ - {c^2}u}&{{c^2}}&0&0&0&0&u \end{array}} \right],$

源项 ${{s}}({{q}}) = {[ - ({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}}){\kern 1pt} r,\;({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}}){\kern 1pt} r,\;0,\,0,\,0,\,0,\,\varPhi ,\,\varPsi ,\,\varPi \,]^{\rm T}},$

式中：u、v为速度分量，c为声速， $H = E + \dfrac{p}{\rho },$ $\dfrac{1}{{\rho {c^2}}} = \dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} + \dfrac{{{\alpha _{\rm{g}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}} + \displaystyle\sum\limits_{k = 1}^K {\dfrac{{{\alpha _k}}}{{{\rho _k}c_k^2}}},$

$\varPhi = \rho {c^2}\dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} \left( {\dfrac{{\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {\rho _{\rm{s}}}}} - \dfrac{{{e_{\rm{s}}}}}{{{\rho _{\rm{s}}}}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} - \dfrac{{\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {\rho _{\rm{g}}}}} - \dfrac{{{e_{\rm{g}}} - Q}}{{{\rho _{\rm{g}}}}}\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {e_{\rm{g}}}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}}} \right) ({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}})r - \dfrac{{\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {\rho _{\rm{s}}}}} - \dfrac{{{e_{\rm{s}}}}}{{{\rho _{\rm{s}}}}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} ({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}})r -$ 　　　　　　　　　　　　　 $\rho {c^2}\dfrac{{{\alpha _{\rm{s}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}}\bigg(\dfrac{{\rho {q_{\rm{s}}}}}{{\rho _{\rm{s}}^2c_{\rm{s}}^2}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}} + \dfrac{{\rho {q_{\rm{g}}}}}{{\rho _{\rm{g}}^2c_{\rm{g}}^2}}\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {e_{\rm{g}}}}} + \displaystyle\sum\limits_{k=1}^{K} {\dfrac{{\rho {q_k}}}{{\rho _k^2c_k^2}}\dfrac{{\partial {p_k}}}{{\partial {e_k}}}}\bigg) + \dfrac{{\rho {q_{\rm{s}}}}}{{\rho _{\rm{s}}^2c_{\rm{s}}^2}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}}$ ，

$\begin{array}{l} \varPsi = \rho {c^2}\dfrac{{{\alpha _k}}}{{{\rho _k}c_k^2}}\left( {\dfrac{{\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {\rho _{\rm{s}}}}} - \dfrac{{{e_{\rm{s}}}}}{{{\rho _{\rm{s}}}}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} - \dfrac{{\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {\rho _{\rm{g}}}}} - \dfrac{{{e_{\rm{g}}} - Q}}{{{\rho _{\rm{g}}}}}\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {e_{\rm{g}}}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}}} \right)({\rho _{\rm{s}}}{\alpha _{\rm{s}}} + {\rho _{\rm{g}}}{\alpha _{\rm{g}}}){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} r - \end{array}$

　　　　　　　　　　　　　 $\rho {c^2}\dfrac{{{\alpha _k}}}{{{\rho _k}c_k^2}}\bigg(\dfrac{{\rho {q_{\rm{s}}}}}{{\rho _{\rm{s}}^2c_{\rm{s}}^2}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}} + \dfrac{{\rho {q_{\rm{g}}}}}{{\rho _{\rm{g}}^2c_{\rm{g}}^2}}\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {e_{\rm{g}}}}} + \displaystyle\sum\limits_{k=1}^{K} {\dfrac{{\rho {q_k}}}{{\rho _k^2c_k^2}}\dfrac{{\partial {p_k}}}{{\partial {e_k}}}} \bigg) + \dfrac{{\rho {q_k}}}{{\rho _k^2c_k^2}}\dfrac{{\partial {p_k}}}{{\partial {e_k}}} ,$

$\varPi \!=\! - \rho {c^2}\left( {\dfrac{{\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {\rho _{\rm{s}}}}} - \dfrac{{{e_{\rm{s}}}}}{{{\rho _{\rm{s}}}}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}}}}{{{\rho _{\rm{s}}}c_{\rm{s}}^2}} - \dfrac{{\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {\rho _{\rm{g}}}}} - \dfrac{{{e_{\rm{g}}} - Q}}{{{\rho _{\rm{g}}}}}\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {e_{\rm{g}}}}}}}{{{\rho _{\rm{g}}}c_{\rm{g}}^2}}} \right)({\rho _{\rm{s}}}{\alpha _{\rm{s}}} \!+\! {\rho _{\rm{g}}}{\alpha _{\rm{g}}}){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} r + \rho {c^2}\bigg(\dfrac{{\rho {q_{\rm{s}}}}}{{\rho _{\rm{s}}^2c_{\rm{s}}^2}}\dfrac{{\partial {p_{\rm{s}}}}}{{\partial {e_{\rm{s}}}}} + \dfrac{{\rho {q_{\rm{g}}}}}{{\rho _{\rm{g}}^2c_{\rm{g}}^2}}\dfrac{{\partial {p_{\rm{g}}}}}{{\partial {e_{\rm{g}}}}} \!+\! \displaystyle\sum\limits_{k=1}^K {\dfrac{{\rho {q_k}}}{{\rho _k^2c_k^2}}\dfrac{{\partial {p_k}}}{{\partial {e_k}}}} \bigg)\!\!\!\!$ 。

非齐次非守恒型双曲律方程组(30)可用Strang分裂方法^[30]来求解，求解过程依次为双曲律方程组求解步与常微分方程组求解步。在双曲律方程组求解步中，一个具有波传播特征的二阶精度Godunov型格式^[31]被采用，该格式能够很好地处理双曲律方程组中的非守恒型对流项。在常微分方程组求解步中，首先根据固体炸药爆轰通常具有快反应过程和慢反应过程的特点^[29]，将爆轰化学反应率分解成快反应率和慢反应率两部分，然后一个具有显-隐离散特性的二阶精度Runge-Kutta格式^[32]被采用，显式离散用于处理慢反应率源项、隐式离散用于处理快反应率源项。

3. 典型算例

3.1 一维爆轰增长问题

考察固体炸药的起爆及爆轰增长过程^[33]，固相反应物和气相生成物均采用刚性气体状态方程 $p = (\gamma - 1)\rho e - \lambda Q$ 形式，爆轰化学反应率采用 $r = {H}\sqrt {1 - \lambda }$ 形式，这里 $\lambda$ 为固相反应物的质量分数，且 $\gamma = 3$ ， $Q = 3.68\;{\rm{GPa}}$ ，H=2.0 μs⁻¹。计算获得爆轰增长过程中压力变化如图1所示，爆轰波头压力单调增加，同时获得定常状态下爆压和爆速的计算值分别约为44.8 GPa、8 460 m/s，与爆压和爆速的理论值45.0 GPa、8 500 m/s相比，计算值与理论值一致。

图 1 一维爆轰的压力增长过程

Figure 1. Growth of pressure in one-dimensional detonation

下载: 全尺寸图片幻灯片

3.2 滑移爆轰约束效应问题

考察一平面爆轰在金属铜的约束下向前传播情况，计算域如图2所示，固体炸药固相反应物和气相生成物均采用Jones-Wilkins-Lee形式状态方程^[34]；爆轰化学反应率采用Ignition-Growth模型^[34]；铜采用Mie-Grüneisen形式状态方程^[18]。起爆区域设置有压力为36.5 GPa的高压静止气相生成物，其余区域均设置为标准状态；计算域入口设置为定值，其余边界设置为无反射边界条件。

图 2 滑移爆轰约束构型图

Figure 2. Configuration of confinement effect

下载: 全尺寸图片幻灯片

计算获得爆轰波达到定常状态时的密度及压力的分布，如图3所示，可以看出，炸药与铜界面附近的爆轰波阵面呈弯曲状态，爆轰波向铜介质折射一激波。数值计算获得定常爆轰状态下炸药与铜界面附近的爆轰波阵面的形态，如图4所示，爆轰定常传播速度计算值为7 670 m/s（理论值为7655 m/s）；爆轰波阵面边缘角 $\alpha$ 的计算值为81.3°（理论值为78.8°），偏转角 $\theta$ 的计算值为4.86°（理论值为4.75°），计算结果与理论结果吻合较好。同时也可以看出，爆轰波向金属进行的折射为正规折射，界面处没有出现反射波。

图 3 铜约束爆轰波传播的密度及压力分布

Figure 3. Distribution of density and pressure in detonation flowfield under copper confinement

下载: 全尺寸图片幻灯片

图 4 铜约束下爆轰波阵面形态

Figure 4. Detonation flowfield nearby explosives under copper confinement

下载: 全尺寸图片幻灯片

3.3 平面爆轰波绕射问题

定常的平面爆轰波绕一个由金属铜构成的90°扩张通道进行传播，计算域如图5所示。固体炸药固相反应物和气相生成物均采用Jones-Wilkins-Lee形式状态方程^[34]；爆轰化学反应率采用Ignition-Growth模型^[34]；铜采用Mie-Grüneisen形式状态方程^[18]。起爆区域设置有压力为36.5 GPa的高压静止气相生成物，其余均为标准状态；计算域入口设置为定值，上边界设置为固壁条件，其余边界设置为无反射边界条件。

图 5 爆轰波绕射构型图

Figure 5. The configuration for the diffraction of detonation wave

下载: 全尺寸图片幻灯片

计算获得爆轰波在t = 4.25，5.00，5.75，6.50 μs这4个时刻的密度和压力分布，如图6所示。

图 6 爆轰波绕射流场图

Figure 6. The flowfield for the diffraction of detonation wave at various simulation times

下载: 全尺寸图片幻灯片

从密度图看，爆轰波在铜介质中产生正规折射激波，该折射激波跨过台阶时会预压炸药；从压力图看，铜介质中折射激波在向炸药折射过程中产生稀疏波引起铜介质出现一个低压区。

4. 结　论

本文中提出一种具有热力学一致性的数值模拟固体炸药爆轰与惰性介质相互作用的扩散界面模型，该模型基于混合网格内各组分物质处于力学平衡与热学非平衡假设，推导获得的模型流动控制系统包括组分物质的质量守恒方程、混合物质的动量及总能量守恒方程，以及组分物质的体积分数演化方程和混合物质的压力演化方程。典型算例表明，数值计算结果与物理规律符合，并且在物质界面附近不会出现物理量的非物理振荡现象、适用于任意表达形式的物质状态方程以及任意数目的惰性介质。

图 1 一维爆轰的压力增长过程

Figure 1. Growth of pressure in one-dimensional detonation

下载: 全尺寸图片幻灯片

图 2 滑移爆轰约束构型图

Figure 2. Configuration of confinement effect

下载: 全尺寸图片幻灯片

图 3 铜约束爆轰波传播的密度及压力分布

Figure 3. Distribution of density and pressure in detonation flowfield under copper confinement

下载: 全尺寸图片幻灯片

图 4 铜约束下爆轰波阵面形态

Figure 4. Detonation flowfield nearby explosives under copper confinement

下载: 全尺寸图片幻灯片

图 5 爆轰波绕射构型图

Figure 5. The configuration for the diffraction of detonation wave

下载: 全尺寸图片幻灯片

图 6 爆轰波绕射流场图

Figure 6. The flowfield for the diffraction of detonation wave at various simulation times

下载: 全尺寸图片幻灯片

参考文献(34)

[1]	NEUMANN J V, RICHTMYER R D. A method for the numerical calculations of hydrodynamical shocks [J]. Journal of Applied Physics, 1950, 21: 232–238. DOI: 10.1063/1.1699639.
[2]	WILKINS M L. Calculation of elastic-plastic flow, methods in computational physics: Vol.3 [M]. New York: Academic Press, 1964: 211−263.
[3]	BENSON D J. Computational methods in Lagrangian and Eulerian hydrocodes [J]. Computer Methods in Applied Mechanics and Engineering, 1992, 99(2−3): 235–394. DOI: 10.1016/0045-7825(92)90042-I.
[4]	BENSON D J. A multi-material Eulerian formulation for the efficient solution of impact and penetration problems [J]. Computational Mechanics, 1995, 15(6): 558–571. DOI: 10.1007/BF00350268.
[5]	GLIMM J, ISAACSON E, MARCHESIN D, et al. Front tracking for hyperbolic systems [J]. Advances in Applied Mathematics, 1981, 2(1): 91–119. DOI: 10.1016/0196-8858(81)90040-3.
[6]	TRYGGVASON G, BUNNER B, ESMAEELI A, et al. A front-tracking method for the computations of multiphase flow [J]. Journal of Computational Physics, 2001, 169: 708–759. DOI: 10.1006/jcph.2001.6726.
[7]	HIRT C, NICHOLS B. Volume of fluid (VOF) method for the dynamics of free boundaries [J]. Journal of Computational Physics, 1981, 39: 201–225. DOI: 10.1016/0021-9991(81)90145-5.
[8]	SAUREL R, ABGRALL R. A multiphase Godunov method for compressible multifluid and multiphase flows [J]. Journal of Computational Physics, 1999, 150: 425–467. DOI: 10.1006/jcph.1999.6187.
[9]	OSHER S, SMEREKA P. A level set approach for computing solutions to incompressible two-phase flow [J]. Journal of Computational Physics, 1994, 114: 146–159. DOI: 10.1006/jcph.1994.1155.
[10]	ASLAM T, BDZIL J, STEWART D. Level set method applied to modeling detonation shock dynamics [J]. Journal of Computational Physics, 1996, 126: 390–409. DOI: 10.1006/jcph.1996.0145.
[11]	ANDERSON D M, MCFADDEN G B, WHEELER A A. Diffuse-interface methods in fluid mechanics [J]. Annual Review of Fluid Mechanics, 1998, 30: 139–165. DOI: 10.1146/annurev.fluid.30.1.139.
[12]	YUE P, FENG J J, LIU C, et al. , A diffuse-interface method for simulating two-phase flows of complex fluids [J]. Journal of Fluid Mechanics, 2004, 515: 293–317. DOI: 10.1017/S0022112004000370.
[13]	PETITPAS F, SAUREL R, FRANQUET E, et al. Modelling detonation waves in condensed energetic materials: multiphase CJ conditions and multidimensional computations [J]. Shock Wave, 2009, 19: 377–401. DOI: 10.1007/s00193-009-0217-7.
[14]	FAVRIS N, GAVRILYUK S, RAUREL R. , Solid-fluid diffuse interface model in cases of extreme deformations [J]. Journal of Computational Physics, 2009, 228: 6037–6077. DOI: 10.1016/j.jcp.2009.05.015.
[15]	SCHOCH S, NIKIFORAKIS N, LEE B, et al. Multi-phase simulation of ammonium nitrate emulsion detonation [J]. Combustion and Flame, 2013, 160: 1883–1899. DOI: 10.1016/j.combustflame.2013.03.033.
[16]	NDANOU S, FAVTIE N, GAVRILYUK S. Multi-solid and multi-fluid diffuse interface model: applications to dynamics fracture and fragmentation [J]. Journal of Computational Physics, 2015, 295: 523–555. DOI: 10.1016/j.jcp.2015.04.024.
[17]	SAUREL R, PANTANO C. Diffuse-interface capturing methods for compressible two-phase flows [J]. Annual Review of Fluid Mechanics, 2018, 50: 105–130. DOI: 10.1146/annurev-fluid-122316-050109.
[18]	MICHAEL L, NIKIFORAKIS N. A hybrid formulation for the numerical simulation of condensed phase explosives [J]. Journal of Computational Physics, 2016, 316: 193–217. DOI: 10.1016/j.jcp.2016.04.017.
[19]	TON V T. Improved shock-capturing methods for multicomponent and reactive flows [J]. Journal of Computational Physics, 1996, 128: 237–253. DOI: 10.1006/jcph.1996.0206.
[20]	SHYUE K. An efficient shock-capturing algorithm for compressible multicomponent problems [J]. Journal of Computational Physics, 1998, 142(1): 208–242. DOI: 10.1006/jcph.1998.5930.
[21]	BANKS J, SCHWENDEMAN D, KAPILA A. A high-resolution Godunov method for compressible multi-material flow on overlapping grids [J]. Journal of Computational Physics, 2007, 223(1): 262–297. DOI: 10.1016/j.jcp.2006.09.014.
[22]	LEE B J, TORO E F, CASTRO C E, et al. Adaptive Osher-type scheme for the Euler equations with highly nonlinear equations of state [J]. Journal of Computational Physics, 2013, 246: 165–183. DOI: 10.1016/j.jcp.2013.03.046.
[23]	BAER M R, NUNZIATO J W. A two-phase mixture theory for the deflagration to detonation transition (DDT) in reactive granular materials [J]. International Journal of Multiphase Flow, 1986, 12: 861–889. DOI: 10.1016/0301-9322(86)90033-9.
[24]	KAPILA A K, MENIKOFF R, BDZIL J B, et al. Two-phase modeling of deflagration-to-detonation transition in granular materials: reduced equations [J]. Physics of Fluids, 2001, 13(10): 3002–3024. DOI: 10.1063/1.1398042.
[25]	ALLAIRE G, CLERC S, KOKH S. A five-equation model for the simulation of interfaces between compressible fluids [J]. Journal of Computational Physics, 2002, 181: 577–616. DOI: 10.1006/jcph.2002.7143.
[26]	MASSONI J, SAUREL R, NKONGA B. Some models and Eulerian methods for interfaces between compressible fluids with heat transfer [J]. International Journal of Heat and Mass Transfer, 2002, 45(6): 1287–1307. DOI: 10.1016/S0017-9310(01)00238-1.
[27]	MURRONE A, GUILLARD H. A five equation reduced model for compressible two phase flow problems [J]. Journal of Computational Physics, 2005, 202(2): 664–698. DOI: 10.1016/j.jcp.2004.07.019.
[28]	GROVE J W. Pressure-velocity equilibrium hydrodynamics models [J]. Acta Mathematica Scientia, 2010, 30B(2): 563–594.
[29]	ZHANG F. Shock wave science and technology reference library: Vol.6 [M]. Berlin: Springer-Verlag Berlin Heidelberg, 2012.
[30]	STRANG G. On the construction and comparison of difference schemes [J]. SIAM Journal of Numerical and Analysis, 1968, 5: 506–517. DOI: 10.1137/0705041.
[31]	LEVEQUE R J. Wave propagation algorithms for multi-dimensional hyperbolic systems [J]. Journal of Computational Physics, 1997, 131(1): 327–353.
[32]	ZHONG X L. Additive semi-implicit Runge-Kutta schemes for computing high-speed nonequilibrium reactive flows [J]. Journal of Computational Physics, 1996, 128: 19–31. DOI: 10.1006/jcph.1996.0193.
[33]	FICKETT W, DAVIS W C. Detonation: theory and experiment [M]. New York: Dover, 1979.
[34]	TARVER C M, MCGUIRE E M. Reactive flow modeling of the interaction of TATB detonation waves with inert materials [C] // The 12th International Symposium on Detonation. San Diego, California, 2002: 641−649.

施引文献

期刊类型引用(1)
1. 吴宗铎，满庆洋，严谨，庞建华，孙一方. 三介质Mie-Grüneisen混合模型在水下爆炸防护层问题中的应用. 计算物理. 2023(05): 556-569 . 百度学术
其他类型引用(0)

资源附件(0)

访问统计

图(6)

计量

文章访问数: 4212
HTML全文浏览量: 1432
PDF下载量: 69
被引次数: 1

1. 扩散界面模型的推导
2. 数值方法
3. 典型算例
3.1 一维爆轰增长问题
3.2 滑移爆轰约束效应问题
3.3 平面爆轰波绕射问题
4. 结　论

固体炸药爆轰与惰性介质相互作用的一种扩散界面模型

doi: 10.11883/bzycj-2019-0435

作者简介: 于 明（1971－ ），男，博士，研究员，yu_ming@iapcm.ac.cn

计量

出版历程

An improved diffuse interface model for the numerical simulation of interaction between solid explosive detonation and inert media

1. 扩散界面模型的推导

2. 数值方法

3. 典型算例

3.1 一维爆轰增长问题

3.2 滑移爆轰约束效应问题

3.3 平面爆轰波绕射问题

4. 结 论

期刊类型引用(1)

其他类型引用(0)

计量

出版历程

目录

1. 扩散界面模型的推导

2. 数值方法

3. 典型算例

3.1 一维爆轰增长问题

3.2 滑移爆轰约束效应问题

3.3 平面爆轰波绕射问题

4. 结 论

作者简介:
于　明（1971－　），男，博士，研究员，yu_ming@iapcm.ac.cn

4. 结　论

4. 结　论