0 引　言

由于水下爆炸在军事和国防领域有着极大的重要性，因而一直受到大家的关注。多数时候，炸药都以球形装药的形式出现。对于球形炸药来说，得到其冲击波压力并不容易。早期对冲击波压力的研究以实验为主^{[1 – 2]}。其后，虽然出现了许多数值模拟，但水下爆炸的实验仍未中断^{[3 – 4]}。但考虑保密性的要求，很多结果未能发表。

然而，由于实验条件成本大，实施难，水下爆炸已经越来越多的被数值模拟所代替。但球形炸药的三维特性使得它在普通的数值模拟中很难得到广泛应用。其中很大一个原因就是普通的数值模拟中设置三维计算网格将消耗更多的计算机资源。因此，关于水下爆炸的研究工作多基于商业软件，如Autodyn^{[5 – 6]}，LS-DYNA^[7]，Abaqus^[8]，MSC-Dytran^[9]等。而目前的水下爆炸数值模拟仍然多停留在二维平面上，如张阿曼的无网格SPH方法^[10]、Daramizadeh^[11]的五方程方法。除此以外，还有一个原因也对球形炸药的数值模拟起着限制作用，那就是涉及到复杂状态方程下的Riemann问题求解。目前对于水下爆炸问题，采用的状态方程多为形式复杂Mie-Grüneisen状态方程^{[5, 7, 10]}，给冲击波的模拟造成很大的难度。

本文将求解坐标系由传统的二维网格体系转到球坐标体系下。由于药球爆炸后冲击波传播的球面对称形式，冲击波的参数仅和传播距离有关^[12]。这样，可将冲击波物理参数表示成球坐标系半径的函数，无需建立三维的直角坐标网格来求解，为计算带来很大的方便。对于Riemann问题的处理，本文在Mie-Grüneisen状态方程下对能量守恒Euler方程做变量分离，通过近似处理把参考压力、参考单位内能等复杂的物理量分离出来，使Euler方程组简单易于计算。利用这种算法，对水下爆炸的近场情况进行数值计算。在激波的传播过程中，整个流场被视为一个理想气体和水的混合体，当中的气体和水用质量分数加以区分。由于药球爆炸后，水和水蒸气共存的状态大量存在于流场周围，该方法能很好地反映这一状态，同时也能较好地计算出水下爆炸周围流场的情况。

1 控制方程

水下爆炸问题中，无论是炸药还是介质，其物理参数密度ρ、压力p和单位体积内能ρe都可写成Mie-Grüneisen状态方程的形式^[12]，其表达式如下：

$p - {p_{ref}}(\rho ) = \rho \varGamma (\rho ,e)\left( {e - {e_{ref}}(\rho )} \right),$

(1)

其中，Г为Mie-Grüneisen系数，可以表示成

$\varGamma (\rho ,e) = {\left. {\frac{1}{\rho }\frac{{\partial p}}{{\partial e}}} \right|_\rho }{\text{。}}$

(2)

因此，Г可以写成密度ρ的函数。p_ref和e_ref为参考点的压力和单位体积内能，与密度ρ有关。假设流场无粘，那么根据质量、动量和能量守恒关系，可以得出欧拉方程：

$\frac{{\partial U}}{{\partial t}} + \frac{{\partial F(U)}}{{\partial x}} = S(U),$

(3)

$U \!=\!\! {r^2}\left[ \begin{array}{l}\rho \\\rho u\\\rho E\end{array} \right],\;F( U) \!=\! {r^2}\left[ \begin{array}{l}\rho u\\\rho {u^2} + p\\(\rho E + p)u\end{array} \right],\;S( U) \!=\! r\left[ \begin{array}{l}0\\2p\\0\end{array} \right]{\text{。}}$

式中：ρ，p，u，E分别为密度，压力，速度和单位体积能量，并且对于单位体积能量E，可以写成内能和动能的总和，于是有

$E = e + \frac{1}{2}{u^2}{\text{。}}$

在炸药和水的接触间断面上，大多数物理参数为间断状态但流体的压力和速度仍然保持连续，间断的物理参数导数可以视为无穷大，因此u和p的导数相对其他导数来说可以忽略不计，那么式（3）中质量守恒和能量守恒可以近似表示成：

$\frac{{\partial \rho }}{{\partial t}} + u\frac{{\partial \rho }}{{\partial x}} = 0,$

(4a)

$\frac{{\partial (\rho e)}}{{\partial t}} + u\frac{{\partial (\rho e)}}{{\partial x}} = 0{\text{。}}$

(4b)

将方程（1）代入式（3）中的能量守恒方程，有

$\frac{\partial }{{\partial t}}\left( {\frac{{p - {p_{ref}}}}{\varGamma } + \rho {e_{ref}}} \right) + u\frac{\partial }{{\partial x}}\left( {\frac{{p - {p_{ref}}}}{\varGamma } + \rho {e_{ref}}} \right) = 0{\text{。}}$

(5)

由于Г，p_ref和e_ref都可以写成只和密度ρ有关的数，再忽略p的变化，因此，可以对式（5）作变量分离：

$\frac{\partial }{{\partial t}}\left( {\frac{1}{\varGamma }} \right) + u\frac{\partial }{{\partial x}}\left( {\frac{1}{\varGamma }} \right) + \rho \left[ {\frac{\partial }{{\partial \rho }}\left( {\frac{1}{\varGamma }} \right)} \right]\frac{{\partial u}}{{\partial x}} = 0,$

(6a)

$\frac{\partial }{{\partial t}}\left( {\frac{{{p_{ref}}}}{\varGamma }} \right) + u\frac{\partial }{{\partial x}}\left( {\frac{{{p_{ref}}}}{\varGamma }} \right) + \rho \left[ {\frac{\partial }{{\partial \rho }}\left( {\frac{{{p_{ref}}}}{\varGamma }} \right)} \right]\frac{{\partial u}}{{\partial x}} = 0,$

(6b)

$\frac{\partial }{{\partial t}}\left( {\rho {e_{ref}}} \right) + u\frac{\partial }{{\partial x}}\left( {\rho {e_{ref}}} \right) + \rho \left[ {\frac{\partial }{{\partial \rho }}\left( {\rho {e_{ref}}} \right)} \right]\frac{{\partial u}}{{\partial x}} = 0{\text{。}}$

(6c)

这里将 $\frac{\partial }{{\partial \rho }}\left( {\frac{1}{\varGamma }} \right),\frac{\partial }{{\partial \rho }}\left( {\frac{{{p_{ref}}}}{\varGamma }} \right),\frac{\partial }{{\partial \rho }}\left( {\rho {e_{ref}}} \right)$ 分别用 $\varphi ,\chi ,\psi $ 来表示。对于多介质混合体，假设单位体积内，第i种介质的质量分数为α_i（i=1.2，…，m），那么该介质在流场内依然满足质量守恒，所以有

$\frac{{\partial ({\alpha _i}\rho )}}{{\partial t}} + \frac{{\partial ({\alpha _i}\rho )u}}{{\partial x}} = 0,\;\;\;\;i = 1,2, \cdot \cdot \cdot {\text{。}}$

(7)

这样Mie-Gruneisen状态方程下的m种介质的控制方程可以表示成包含ρ，p，u，E和Г，p_ref，e_ref以及质量分数α_i的如下的欧拉关系式。方程（3）中的变量U可以写成：

$U = {r^2}{\left[ {\rho ,\rho u,\rho E,\frac{1}{\varGamma },\frac{{{p_{ref}}}}{\varGamma },\rho {e_{ref}},\rho {\alpha _i}} \right]^{\rm T}}{\text{。}}$

通量F则表示成：

$\begin{array}{l}F = {r^2}\left[ {\rho u,\;\rho {u^2} + p,\;\left( {\rho E + p} \right)u,\;\displaystyle\frac{1}{\varGamma }u + \rho \varphi u,} \right.\\\left. {\displaystyle\frac{{{p_{ref}}}}{\varGamma }u + \rho \chi u,\;\rho {e_{ref}}u + \rho \psi u,\;\rho {\alpha _i}u} \right]^{\rm T}{\text{。}}\end{array}$

这样， $\frac{1}{\varGamma },\frac{{{p_{ref}}}}{\varGamma },\rho {e_{ref}}$ 成为3个与ρ，p，u，E线性无关的变量。当Г，p_ref和e_ref解析式比较复杂时，这样处理能使计算变得简单。

2 计算方法

本文采用MUSCL-TVD格式的有限体积法进行计算。假设问题为一维Euler方程，其有限体积形式有：

$\frac{{\overline U _j^{n + 1} - \overline U _j^n}}{{Dt}} + \frac{{{{\overline F }_{j + 1/2}} - {{\overline F }_{j - 1/2}}}}{{Dx}} = 0{\text{。}}$

式中：Dt，Dx分别为时间和空间步长， $\overline U _j^n$ 为第n时间步节点j处的守恒变量， $\overline F _{j + 1/2}^n$ 为节点j和节点j+1之间的数值通量。对于通量 $\overline F _{j + 1/2}^n$ 的计算，采取2阶精度的TVD格式：

$\begin{split}\bar F_{j + 1/2}^n = \displaystyle\frac{1}{2}\left[ {F\left( {U_{L,j + 1/2}^n} \right) + F\left( {U_{R,j + 1/2}^n} \right) - } \right.\\[8pt]\quad\quad\quad\quad\!\!\! \left. {\hat R\left| {\hat \Lambda } \right|\hat L\left( {U_{R,j + 1/2}^n - U_{L,j + 1/2}^n} \right)} \right]{\text{。}}\end{split}$

$\widehat \Lambda ,\widehat L,\widehat R$ 为Roe解下的特征值和对应的左右特征向量构成矩阵。对于守恒变量 $U_{L,j + 1/2}^n$ 和 $U_{R,j + 1/2}^n$ ，有

$\left\{ \begin{array}{l}U_{L,j + 1/2}^n = U_j^n + \frac{1}{2}{R_{j + 1/2}}{W_{L,j}},\\U_{R,j + 1/2}^n = U_{j + 1}^n - \frac{1}{2}{R_{j + 1/2}}{W_{R,j}}{\text{。}}\end{array} \right.$

其中： ${R_{j + 1/2}}$ 为特征阵的右特征矩阵，列向量的 ${W_{L,j}},{W_{R,j}}$ 子元素 ${w_{L,j}},{w_{R,j}}$ 有：

$\left\{ \begin{array}{l}{w_{L,j}} = D_j^n\left( {1 - {\lambda _{j + 1/2}}\frac{{Dt}}{{Dx}}} \right),\\{w_{R,j}} = D_{j + 1}^n\left( {1 + {\lambda _{j + 1/2}}\frac{{Dt}}{{Dx}}} \right){\text{。}}\end{array} \right.$

式中： $D_j^n \!=\! \min od \left[ {{l_{j \!+\! 1/2}}\left( {{U_{j \!+\! 1}} - {U_j}} \right),{l_{j + 1/2}}\left( {{U_j} \!-\! {U_{j \!-\! 1}}} \right)} \right]$ ， ${l_{j + 1/2}}$ 和 ${\lambda _{j + 1/2}}$ 分别为左特征向量和特征值。

3 计算实例

半径为R₀裸药球在无限水域中爆炸，假设流体无粘且药球质量分布均匀，t=0时刻，炸药爆炸并形成高压气体，气体密度ρ₀等于炸药球的密度，即1 630 kg/m³，爆压pH=2.1×10¹⁰ Pa，气团半径等于药球半径^[1]。假设TNT炸药爆炸后形成的气体为理想气体，状态方程如下：

$p = A\left( {1 - \frac{{\omega \theta }}{{{R_1}}}} \right){e^{ - {R_1}/\theta }} + B\left( {1 - \frac{{\omega \theta }}{{{R_2}}}} \right){e^{ - {R_2}/\theta }} + \omega \theta {\rho _0}e{\text{。}}$

这里，θ=ρ/ρ₀。水的状态方程用多项式形式来表达：

$\left\{ \begin{array}{l}p = {a_1}\mu + {a_2}{\mu ^2} + {a_3}{\mu ^3} + \left( {{b_0} + {b_1}\mu + {b_2}\mu + {b_3}\mu } \right){\rho _0}E,\;\mu > 0;\\p = {a_1}\mu + \left( {{b_0} + {b_1}\mu } \right){\rho _0}E,\;\mu < 0{\text{。}}\end{array} \right.$

式中：μ=ρ/ρ₀–1，其系数如表1所示。

这样，转化成如（1）所示的Mie-Grüneisen的形式以后，相关参数为：

$\begin{array}{l}{p_{ref}} = \left\{ \begin{array}{l}{a_1}\mu + {a_2}{\mu ^2} + {a_3}{\mu ^3},\;\mu > 0\\{a_1}\mu ,\;\mu < 0\end{array} \right.,\;{e_{ref}} = 0,\\\varGamma = \left\{ \begin{array}{l}\left( {{b_0} + {b_1}\mu + {b_2}{\mu ^2} + {b_3}{\mu ^3}} \right)\rho /{\rho _0},\;\mu > 0\\\left( {{b_0} + {b_1}\mu } \right)\rho /{\rho _0},\;\mu < 0\end{array} \right.\text{。}\end{array}$

图1为压力峰值与文献[2]的对比结果。从图中可以看出，2条压力峰值曲线与大致保持一致，且在距爆点约12R₀处吻合得比较好。当冲击波过了12R₀以后，冲击波的压力衰减速度降低，曲线将有所缓和，与本文的计算存在少量差距。而在距爆点很近的地方（6R₀以内），误差会有所增加。由于这部分区域的爆炸冲击波的作用很强，压力变化非常剧烈，有些实验结果甚至未能在6R₀以内给出合理的经验公式^[1]。

图2为距爆心不同距离处的压力峰值曲线，横坐标是测点与爆心的距离R和初始半径R₀的比值。本文计算结果与文献[2]的经验公式相比，近场压力值总体有些偏大。但基本反映了压力在冲击波前后的变化情况。在冲击波传播到某点时，该点压力瞬间上升到峰值；随后冲击波继续前进，此时由于冲击波作用减弱该点压力则迅速下降。随后，向内收缩的稀疏波聚集到爆心处又重新扩散开来，造成压力的脉动变化，形成第2次峰值。

图3为不同时刻压力的分布图，从图中可以看出，随着时间的推移，压力曲线逐渐变得平稳，并在最大峰值过后一定距离出现了压力脉动，形成二次脉动压力。图4给出了不同位置的二次压力峰值，可以看出其数值比较小。值得注意的是，由于本文的计算是基于能量守恒的假设，未考虑能量的损耗，因此本文的计算结果比实际情况要大，由脉动压力形成的第2次的压力峰值也比较明显。实际情况下，第2次压力峰值很难在曲线上清楚显现出来^[1]。

图5理想气体的质量分数α和质量ρα在不同时刻的变化情况，从图5中可以看出，冲击波传播开来以后，中央的气团逐渐由单一的气体状态变为气体和液态水并存的一种状态。这情况与实际的气水模型存在一定差别，实际的气水模型中并未有爆轰气体溶解在水中，因此本计算模型仍然存在着一些数值上的不准确，但并不影响冲击波压力的计算。

4 结　语

利用本文的分离变量方法能很好地解决Mie-Grüneisen状态方程下的多介质可压缩流动的激波问题，并且能有效地简化计算过程。该方法在处理水下爆炸近场问题的时候能取得很好的效果。按本文的计算方法处理得到的计算结果和文献[1]对比，取得了不错的效果，并对计算第2次的脉动压力有着很好的帮助。

炸药在水中爆炸以后，水的状态包括液态水、水蒸气和液气混合态。再加上爆炸产生的气体、液态水、气体和气水混合体共存于流场中。利用质量分数可以区分开流场中的各种成分。炸药爆炸后，早期形成高密度气团，随后气团扩散很快，在气团扩散过程中，气体的质量分数和密度都急剧下降，最后形成一个中心大部分区域稀薄，两边小部分区域稠密的状态。

由于本文仅考虑了动量和能量守恒的理想情况，而实际情况比较复杂，能量损耗不可避免，因此，本文所计算的压力值比实际数值大，冲击波的传播速度也比实际要快。