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=(Γi−1)ρie−Πip=\left( \varGamma _i-1 \right) \rho _ie-\varPi _ip=(Γi−1)ρie−Πi
式中，Γi为Grüneisen系数，Πi为与理想气体的压力偏差项，与具体状态方程形式有关。
用户手册给出的JWL状态方程形式为
式（2）Πi=Aie−R1,iV(ωiR1,iV−1)+Bie−R2,iV(ωiR2,iV−1)−ρiωie0,i\varPi _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 _ie_{0,i}Πi=Aie−R1,iV(R1,iVωi−1)+Bie−R2,iV(R2,iVωi−1)−ρiωie0,i
式中，V=ρ0,i/ρi。
将式（2）代入式（1）可得
式（3）p=(Γi−1)ρie−[Aie−R1,iV(ωiR1,iV−1)+Bie−R2,iV(ωiR2,iV−1)−ρiωie0,i]=Aie−R1,iV(1−ωiR1,iV)+Bie−R2,iV(1−ωiR2,iV)+ρiωi(e0,i+e)=A(1−ωR1V)exp⁡(−R1V)+B(1−ωR2V)exp⁡(−R2V)+ρω(e0+e)p=\left( \varGamma _i-1 \right) \rho _ie-\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 _ie_{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_1V} \right) \exp \left( -R_1V \right) +B\left( 1-\frac{\omega}{R_2V} \right) \exp \left( -R_2V \right) +\rho \omega \left( e_0+e \right) p=(Γi−1)ρie−[Aie−R1,iV(R1,iVωi−1)+Bie−R2,iV(R2,iVωi−1)−ρiωie0,i]=Aie−R1,iV(1−R1,iVωi)+Bie−R2,iV(1−R2,iVωi)+ρiωi(e0,i+e)=A(1−R1Vω)exp(−R1V)+B(1−R2Vω)exp(−R2V)+ρω(e0+e)
式中，Γi=ω+1；A，B，R1，R2和ω均为状态方程系数；e0为补偿内能，即参考状态下的内能。
文献中标准的JWL状态方程为
式（4）
p=A(1−ωR1V)exp⁡(−R1V)+B(1−ωR2V)exp⁡(−R2V)+ωeVp=A\left( 1-\frac{\omega}{R_1V} \right) \exp \left( -R_1V \right) +B\left( 1-\frac{\omega}{R_2V} \right) \exp \left( -R_2V \right) +\frac{\omega e}{V}p=A(1−R1Vω)exp(−R1V)+B(1−R2Vω)exp(−R2V)+Vωe
通过对比式（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）
pi=(γi−1)ρiei−γi(bi−ai)p_i=\left( \gamma _i-1 \right) \rho _ie_i-\gamma _i\left( b_i-a_i \right) pi=(γi−1)ρiei−γi(bi−ai)
在水下爆炸理论中，常用的Tait等熵状态方程为
式（6）
(ρρ0)n=p+Bˉp0+Bˉ=p+B−AB\left( \frac{\rho}{\rho _0} \right) ^n=\frac{p+\bar{B}}{p_0+\bar{B}}=\frac{p+B-A}{B}(ρ0ρ)n=p0+Bˉp+Bˉ=Bp+B−A
文献《A Non-oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (the Ghost Fluid Method)》中给出了式（6）对应的内能为
式（7）
e=Bρn−1(n−1)ρ0n+B−Aρe=\frac{B\rho ^{n-1}}{\left( n-1 \right) \rho _{0}^{n}}+\frac{B-A}{\rho}e=(n−1)ρ0nBρn−1+ρB−A
结合式（6），上式可简化为
式（8）
e=Bρn−1(n−1)ρ0n+B−Aρ=B(n−1)ρρnρ0n+B−Aρ=p+B−A(n−1)ρ+B−Aρ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}e=(n−1)ρ0nBρn−1+ρB−A=(n−1)ρBρ0nρn+ρB−A=(n−1)ρp+B−A+ρB−A
即
式（9）
(n−1)ρe=p+B−A+(n−1)(B−A)=p+n(B−A)p=(n−1)ρe−n(B−A)\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) (n−1)ρe=p+B−A+(n−1)(B−A)=p+n(B−A)p=(n−1)ρe−n(B−A)
通过对比式（5）和式（9）可以发现，二者完全一致，只是符号表示有所不同罢了。文献中同样给出了水介质对应的参数值，即γi=n=7.15，bi=B=3.31×108Pa，ai=A=101325Pa，ρ0=1000kg/m3。对Tait状态方程更多的分析可参考贾曦雨博士的学位论文《水下爆炸近场冲击波研究》。

3、 Linear Tilloston状态方程

用户手册中与多项式状态方程类似的状态方程为Linear Tilloston状态方程，形式如下
式（10）
pi=p0,i+ωiρi(eT−e0,i)+Aiμ+Biμ2+Ciμ3=p0+ωρ(eT−e0)+Aμ+Bμ2+Cμ3p_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 ^3pi=p0,i+ωiρi(eT−e0,i)+Aiμ+Biμ2+Ciμ3=p0+ωρ(eT−e0)+Aμ+Bμ2+Cμ3
式中，p0和e0分别为参考压力和内能。
目前，数值仿真软件中常用的多项式状态方程形式为
式（11）
p={A1μ+A2μ2+A3μ3+(B0+B1μ)ρ0eM,μ⩾0T1μ+T2μ2+B0ρ0eM,μ<0p=\begin{cases} A_1\mu +A_2\mu ^2+A_3\mu ^3+\left( B_0+B_1\mu \right) \rho _0e_M,\mu \geqslant 0\\ T_1\mu +T_2\mu ^2+B_0\rho _0e_M,\mu <0\\\end{cases}p={A1μ+A2μ2+A3μ3+(B0+B1μ)ρ0eM,μ⩾0T1μ+T2μ2+B0ρ0eM,μ<0
式（12）
p={a1μ+a2μ2+a3μ3+(b0+b1μ+b2μ2+b3μ3)ρ0eM,μ⩾0a1μ+(b0+b1μ)ρ0eM,μ<0p=\left\{ \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 _0e_M,\mu \geqslant 0\\ a_1\mu +\left( b_0+b_1\mu \right) \rho _0e_M,\mu <0\\\end{array} \right. p={a1μ+a2μ2+a3μ3+(b0+b1μ+b2μ2+b3μ3)ρ0eM,μ⩾0a1μ+(b0+b1μ)ρ0eM,μ<0
式（13）
p={C0+C1μ+C2μ2+C3μ3+(C4+C5μ+C6μ2)eV,μ⩾0C0+C1μ+C3μ3+(C4+C5μ)eV,μ<0p=\left\{ \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 \geqslant 0\\ C_0+C_1\mu +C_3\mu ^3+\left( C_4+C_5\mu \right) e_V,\mu <0\\\end{array} \right. p={C0+C1μ+C2μ2+C3μ3+(C4+C5μ+C6μ2)eV,μ⩾0C0+C1μ+C3μ3+(C4+C5μ)eV,μ<0
通过对比可以发现，式（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水下爆炸数值仿真文件”。