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

于 明

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

爆轰是复杂的流体力学与剧烈的化学动力学耦合的一种物理现象。固体炸药的爆轰能够提供强大的能量来驱动和压缩材料，在国防科技和国民经济中获得广泛应用，因此，固体炸药爆轰与惰性介质相互作用问题是工程应用中十分重要的问题。固体炸药爆轰与惰性介质相互作用问题本质上是一种可压缩多物质流动问题，除了可压缩单物质流动中出现的激波、稀疏波、接触间断外，一种新的间断，即物质界面出现了，它隔开热力学性质或状态方程不一样的两种物质，这种物质界面实质上也是一种接触间断，它满足物质界面压力相等与速度相等的条件。对可压缩多物质流动，传统的数值模拟方法多是拉格朗日方法[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表示；其中 ρ 表示密度，e表示内能，p表示压力，u表示速度矢量， α 表示体积分数，v表示比容，Q表示单位质量的固体炸药由固相反应物转化为气相生成物时所释放的热量。在力学平衡条件下，混合物质物理量与组分物质物理量的关系式为：

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

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

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

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

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

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

2 数值方法

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

这里状态变量q=[ρsαs, ρgαg, ρkαk, ρu, ρv, ρE, αs, αk,p]T（不失一般性，惰性介质只写一种物质），矩阵A(q) 与矩阵B(q) 是具有类似表达形式的Jacobi 矩阵，其中

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

3 典型算例

3.1 一维爆轰增长问题

图1 一维爆轰的压力增长过程Fig. 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 滑移爆轰约束构型图Fig. 2 Configuration of confinement effect

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

图3 铜约束爆轰波传播的密度及压力分布Fig. 3 Distribution of density and pressure in detonation flowfield under copper confinement

图4 铜约束下爆轰波阵面形态Fig. 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 爆轰波绕射构型图Fig. 5 The configuration for the diffraction of detonation wave

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

图6 爆轰波绕射流场图Fig. 6 The flowfield for the diffraction of detonation wave at various simulation times

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

4 结 论

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