blastFoam水下爆炸数值仿真计算

目前，常用的水下爆炸数值仿真软件均为商业软件，例如AUTODYN、LS-DYNA和MSC.DYTRAN等。基于开源OpenFOAM框架，blastFoam提供了可用于进行爆炸数值仿真的模块，其源代码下载地址为https://github.com/synthetik-technologies/blastfoam，安装及其使用可参考其用户手册《blastFoam Theory and User Guide》。
本文的主要目的是建立blastFoam对应的水下爆炸数值仿真模型并进行计算。由于blastFoam已经给出装药空中爆炸对应的算例（对应路径为/home/zyh/ OpenFOAM/blastfoam/tutorials/blastFoam/axisymmeticCharge），因此本文将在此基础上对部分内容进行修改，主要集中于流体介质的状态方程，即将空气采用的理想气体状态方程修改为水介质采用的Tait状态方程或者多项式状态方程，此外，删除了原有算例中对自适应网格的设置，将C4装药改为了TNT。

1、JWL状态方程

用户手册中，Mie-Grüneisen状态方程的一般形式为
式（1）$$p=\left({\Gamma }_i-1\right){\rho }_{i}e-{\Pi }_i$$
式中，Γi为Grüneisen系数，Πi为与理想气体的压力偏差项，与具体状态方程形式有关。
用户手册给出的JWL状态方程形式为
式（2）$${\Pi }_i=A_i{\mathrm{e}}^{-R_{1,i}V}\left(\frac{{\omega }_i}{R_{1,i}V}-1\right)+B_i{\mathrm{e}}^{-R_{2,i}V}\left(\frac{{\omega }_i}{R_{2,i}V}-1\right)-{\rho }_i{\omega }_i{e}_{0,i}$$
式中，V=ρ0,i/ρi。
将式（2）代入式（1）可得
式（3）$$\begin{gathered} p=\left({\Gamma }_i-1\right){\rho }_{i}e-\left[A_i{\mathrm{e}}^{-R_{1,i}V}\left(\frac{{\omega }_i}{R_{1,i}V}-1\right)+B_i{\mathrm{e}}^{-R_{2,i}V}\left(\frac{{\omega }_i}{R_{2,i}V}-1\right)-{\rho }_i{\omega }_i{e}_{0,i}\right] \\ =A_i{\mathrm{e}}^{-R_{1,i}V}\left(1-\frac{{\omega }_i}{R_{1,i}V}\right)+B_i{\mathrm{e}}^{-R_{2,i}V}\left(1-\frac{{\omega }_i}{R_{2,i}V}\right)+{\rho }_i{\omega }_i\left(e_{0,i}+e\right) \\ =A\left(1-\frac{\omega }{R_{1}V}\right)\exp \left(-R_{1}V\right)+B\left(1-\frac{\omega }{R_{2}V}\right)\exp \left(-R_{2}V\right)+\rho \omega \left(e_0+e\right) \end{gathered}$$
式中，Γi=ω+1；A，B，R1，R2和ω均为状态方程系数；e0为补偿内能，即参考状态下的内能。
文献中标准的JWL状态方程为
式（4）
$$p=A\left(1-\frac{\omega }{R_{1}V}\right)\exp \left(-R_{1}V\right)+B\left(1-\frac{\omega }{R_{2}V}\right)\exp \left(-R_{2}V\right)+\frac{\omega e}{V}$$
通过对比式（3）和式（4）可以发现，除内能的物理含义不同外，二者相同。在式（3）中，内能是单位质量内能，式（4）中内能项为单位体积内能。因此，借助式（3）描述高能装药爆轰产物做工能力时，可直接借助式（4）给出的参数。
对于高能装药，还需设定其爆轰参数。特别需要注意的是，此时初始内能为单位体积内能，而非状态方程中的单位质量内能。表1给出了TNT对应的JWL状态方程和爆轰参数。
表1 TNT炸药JWL状态方程及爆轰参数
A/GPa B/GPa R1 R2 ω pCJ/GPa D/(m∙s-1) ρ/(kg/m3) eV0/GPa
371.2 3.231 4.15 0.95 0.3 21 6930 1630 7.0

2、Tait状态方程

用户手册中Tait状态方程形式如下：
式（5）
$$p_i=\left({\gamma }_i-1\right){\rho }_i{e}_i-{\gamma }_i\left(b_i-a_i\right)$$
在水下爆炸理论中，常用的Tait等熵状态方程为
式（6）
$${\left(\frac{\rho }{{\rho }_0}\right)}^n=\frac{p+\bar{B}}{p_0+\bar{B}}=\frac{p+B-A}{B}$$
文献《A Non-oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (the Ghost Fluid Method)》中给出了式（6）对应的内能为
式（7）
$$e=\frac{B{\rho }^{n-1}}{\left(n-1\right){\rho }_0^n}+\frac{B-A}{\rho }$$
结合式（6），上式可简化为
式（8）
$$\begin{gathered} e=\frac{B{\rho }^{n-1}}{\left(n-1\right){\rho }_0^n}+\frac{B-A}{\rho }=\frac{B}{\left(n-1\right)\rho }\frac{{\rho }^n}{{\rho }_0^n}+\frac{B-A}{\rho } \\ =\frac{p+B-A}{\left(n-1\right)\rho }+\frac{B-A}{\rho } \end{gathered}$$
即
式（9）
$$\begin{gathered} \left(n-1\right)\rho e=p+B-A+\left(n-1\right)\left(B-A\right)=p+n\left(B-A\right) \\ p=\left(n-1\right)\rho e-n\left(B-A\right) \end{gathered}$$
通过对比式（5）和式（9）可以发现，二者完全一致，只是符号表示有所不同罢了。文献中同样给出了水介质对应的参数值，即γi=n=7.15，bi=B=3.31×108Pa，ai=A=101325Pa，ρ0=1000kg/m3。对Tait状态方程更多的分析可参考贾曦雨博士的学位论文《水下爆炸近场冲击波研究》。

3、 Linear Tilloston状态方程

用户手册中与多项式状态方程类似的状态方程为Linear Tilloston状态方程，形式如下
式（10）
$$\begin{gathered} p_i=p_{0,i}+{\omega }_i{\rho }_i\left(e_T-e_{0,i}\right)+A_i\mu +B_i{\mu }^2+C_i{\mu }^3 \\ =p_0+\omega \rho \left(e_T-e_0\right)+A\mu +B{\mu }^2+C{\mu }^3 \end{gathered}$$
式中，p0和e0分别为参考压力和内能。
目前，数值仿真软件中常用的多项式状态方程形式为
式（11）
$$p=\{\begin{array}{l}{A}_1\mu +A_2{\mu }^2+A_3{\mu }^3+\left(B_0+B_1\mu \right){\rho }_0{e}_M,\mu ⩾0\\ T_1\mu +T_2{\mu }^2+B_0{\rho }_0{e}_M,\mu <0\end{array}$$
式（12）
$$p=\{\begin{array}{l}{a}_1\mu +a_2{\mu }^2+a_3{\mu }^3+\left(b_0+b_1\mu +b_2{\mu }^2+b_3{\mu }^3\right){\rho }_0{e}_M,\mu ⩾0\\ a_1\mu +\left(b_0+b_1\mu \right){\rho }_0{e}_M,\mu <0\end{array}$$
式（13）
$$p=\{\begin{array}{l}{C}_0+C_1\mu +C_2{\mu }^2+C_3{\mu }^3+\left(C_4+C_5\mu +C_6{\mu }^2\right)e_V,\mu ⩾0\\ C_0+C_1\mu +C_3{\mu }^3+\left(C_4+C_5\mu \right)e_V,\mu <0\end{array}$$
通过对比可以发现，式（10）与式（13）最为接近，不过式（10）中与内能项有关的参数较小。水介质对应式（13）的参数如表2所示，将该组参数应用于式（10）时令C5=0，该操作将引入误差，未来应寻找式（10）对应的参数或者直接扩展blastFoam对应的状态方程，即将式（13）引入blastFoam。
表2 水介质状态方程参数
C0/Pa C1/GPa C2/GPa C3/GPa C4 C5 C6
101325 2.2 9.54 14.57 0.28 0.28 0

4、仿真结果

状态方程设置位于文件case/constant/phaseProperties。除状态方程外，还需在该文件内修改热力学模型及其参数，该部分主要参考自用户手册及相关教程案例，其可靠性本文并未考证（在商业软件中，根据已有经验，热力学模型很少使用）。对于网格、离散格式、线性方程求解算法、时间控制、观测点设置等，主要参考自原案例，只对部分参数进行了修改，例如将网格尺寸由1m改为0.1m，取消了自适应网格等，不再赘述。
t=0.01s时的压力场如图1所示。可以看出，采用不同水介质状态方程时，压力仿真结果有一定的差别。

(a)Tait

(b)Linear Tilloston
图1 压力场仿真结果
图2给出了5倍装药半径处的压力时间曲线。可以看出，不同状态方程对应的曲线在衰减阶段的大部分时间内重合，但波峰处差距较大。

图2 5倍装药半径处压力时间曲线

5、 小结

本文只是对blastFoam在自由场水下爆炸问题中应用的初步探索，并未对仿真结果的精度进行验证，因此可能存在一些问题。未来将对blastFoam给出的仿真结果与商业软件进行对比，并进行二次开发，进一步拓展blastFoam的使用范围。

---

blastFoam输入文件可在CSDN搜索下载，资源名为“blastFoam水下爆炸数值仿真文件”。