球坐标系下多介质混合物模型的数值模拟

吴宗铎; 赵勇; 严谨; 宗智; 高云

doi:10.11883/bzycj-2018-0075

球坐标系下多介质混合物模型的数值模拟

doi: 10.11883/bzycj-2018-0075

吴宗铎^1,,
赵勇^2, ,,
严谨¹,
宗智³,
高云⁴

1.
广东海洋大学海洋工程学院，广东湛江 524088
2.
大连海事大学交通运输工程博士后流动站，辽宁大连 116026
3.
大连理工大学船舶工程学院，辽宁大连 116024
4.
西南石油大学油气藏地质及开发工程国际重点实验室，四川成都 610500

基金项目: 国家自然科学基金（11702066，51679021，51609206，51679037）；广东省自然科学基金（2017A030313275）

详细信息

作者简介:
吴宗铎（1984－　），男，博士，讲师，wuzongduo0@aliyun.com

通讯作者:
赵　勇（1981－　），男，博士，副教授，fluid@126.com

中图分类号: O.382; O359+.1
计量
- 文章访问数: 5604
- HTML全文浏览量: 1513
- PDF下载量: 41
- 被引次数: 8
出版历程
- 收稿日期: 2018-03-12
- 修回日期: 2018-06-27
- 刊出日期: 2019-05-01

Numerical simulation about the multi-component mixture model under spherical coordinate system

WU Zongduo^1
,,
ZHAO Yong^{2
, ,},
YAN Jin¹,
ZONG Zhi³,
GAO Yun⁴

1.
College of Ocean Engineering, Guangdong Ocean University, Zhanjiang 524088, Guangdong, China
2.
Transportation Engineering postdoctoral research station, Dalian Maritime University, Dalian 116026, Liaoning, China
3.
Department of Naval Architecture, Dalian University of Technology, Dalian 116024, Liaoning, China
4.
State Key Laboratory of Oil and Gas Reservoir Gelogy and Exploration, Southwest Petroleum University,
Chengdu 610500, Sichuan, China

摘要

摘要: 利用多介质混合模型在求解球坐标系下的Riemann问题时，需要考虑界面处压力平衡性弱、奇点处理、状态方程复杂等多个难点。本文将原始基于体积分数的Mie-Grüneisen多介质混合模型扩展到球坐标系下，并对多个细节进行了修正和改进，包括：在界面处对热力学参数进行修正、采用质量分数导出新输送方程、利用质量分数加权计算偏导数、采用相邻网格点的物理量定义奇点等。经过改进后的计算模型，可以得到无振荡的数值解，而且可以准确捕捉到冲击波和界面的位置。另外，使用改进后的质量分数模型比原始的体积分数模型得到的计算结果更准确。
- 多介质 /
- 混合模型 /
- Mie-Grüneisen模型界面的位置方程 /
- 球坐标
Abstract: The aim of the paper is to extend the Mie-Grüneisen mixture model to spherical coordinate. As the multi-component mixture model is applied to the Riemann problem under spherical coordinate, many problems need to be taken into account: weak equilibrium, singular point treatment, complex equations of states and so on. In the article, the research work starts from the Mie-Grüneisen mixture model, then extend to the revision and modification about many details, include: revision of the thermo dynamical parameters at interface, deduction of new transport equation by mass fraction, weighting evaluation of partial derivatives by mass fraction, definition of physical parameters by the adjacent grid for singular point and other so on. The seriously modified numerical model, can not only obtain non-oscillation solutions, but also catch the positions of shock wave and interface clearly. In addition, the modified mass fraction model, can get more accurate results than the original model with mass fraction.
- multi-component /
- mixture model /
- Mie-Grüneisen equation /
- spherical coordinate

HTML全文

多介质的激波间断问题一直以来都是计算流体力学方面的一个备受关注的问题，目前广泛应用的一个处理方法是利用Level Set函数追踪气水界面，同时结合近似Riemann解下的虚拟流体质点完成对冲击波参数的计算^[1]。虽然得到的结果较理想，但是计算过程复杂，且计算依赖近似Riemann解，无法在JWL (Jones-Wilkins-Lee)及Mie-Grüneisen等复杂形式状态方程下展开更深入的研究^[2-3]。

相比之下，利用多介质混合模型，可以将流场视为一个整体并用体积分数或质量分数来区分流体，从而在不直接求近似Riemann解的前提下完成计算。多介质混合物模型早期由Abgrall提出^[4-5]，通过非守恒方程来计算由气体或低压缩度的固液体组成的流场^[5-6]。随后，Abgrall将该混合物模型扩展成一个七方程模型，能适应一般形式状态方程的流场^[7]。柏劲松等^[8]也利用此类多介质混合物方程来模拟两种以上介质的冲击过程，且配合等效方程可以适应一些复杂形式的状态方程^[9]。梁珊等^[10]和Liang等^[11]将Abgrall的七方程模型与多GPU的并行计算相结合来提高计算精确度。为了让算法更好的适应并行计算，刘娜等^[12]将谱体积方法运用到多介质混合物模型中，提升了计算效率。

虽然该多介质混合物模型在采用直角坐标系处理平面波的问题时，取得不错的进展，但是对于常用于处理球面波的球坐标系下的多介质问题，仍然具备一定的局限性。第一，多介质混合物模型建立在速度和压力在界面处维持平衡的基础，仅对平面波适用，但是球面波的压力和速度都会随着传播半径的增加而衰减；第二，球坐标系下，原点位置处为冲击波运动的奇点位置，必须作出合适的处理；第三，为了增加工程实用性，计算模型需要能适应复杂形式的状态方程，而不仅仅是理想气体方程等形式简单的状态方程。

考虑到以上因素，本文将直角坐标系下的基于Mie-Grüneisen状态方程下的多介质混合模型延伸到球坐标系下。为了压力和速度不平衡问题以及坐标原点的奇点问题，在原计算模型下作了参数修正并用质量分数替换体积分数，并对奇点处数值差分作了特殊处理—利用相邻网格点参数赋值给奇点参数但奇点处速度设定为0。最后，将通过数值算例验证其数值稳定性和计算准确性。

1. 直角坐标系下的计算模型

考虑多介质的可压缩流体组成的混合流场，其运动形式可以用欧拉方程来描述：

$\left\{ \begin{aligned} &\frac{{\partial \rho }}{{\partial t}} + \frac{\partial }{{\partial x}}\left( {\rho u} \right) = 0\\ &\frac{{\partial \left( {\rho u} \right)}}{{\partial t}} + \frac{\partial }{{\partial x}}\left( {\rho {u^2} + p} \right) = 0\\ &\frac{{\partial \left( {\rho E} \right)}}{{\partial t}} + \frac{\partial }{{\partial x}}\left( {\rho Eu + pu} \right) = 0 \end{aligned} \right.$

(1)

式中：ρ、p、u、E分别为密度、压力、速度和单位质量能量^[13]。单位质量的总能量为内能和动能的总和，且E=e+u²/2。

流场中介质的物理特性，可以用Mie-Grüneisen状态方程来表述：

$p - {p_{{\rm{ref}}}}(\rho ) = \rho \varGamma (\rho ,e)[e - {e_{{\rm{ref}}}}(\rho )]$

(2)

式中：Γ为Mie-Grüneisen系数，p_ref和e_ref为参考点的压力和单位质量内能，与密度ρ有关。

1.1 界面处热力学关系

界面处，流体的压力p和速度u均保持平衡^[4]，此时p与u在界面处不发生变化，而密度ρ与内能e则出现数值间断。因此∂p/∂x和∂u/∂x约为0，式(1)中的质量守恒与能量守恒可以分别表示为：

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

(3)

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

(4)

将式(2)代入式(4)，有：

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

(5)

由于Γ、p_ref、e_ref都只和密度ρ有关的函数，且式(5)对于任意的表达式p_ref(ρ)和e_ref(ρ)均成立。这里首先假设p_ref、e_ref都为0，那么有^[2,9]

$\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$

(6)

同理，可依次类推得到关于p_ref和e_ref的非守恒形式的方程：

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

(7)

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

(8)

1.2 波面处热力学关系

由于式(6)～(8)的推导过程需要基于p和u在界面处连续这一前提条件，但是在左右波面处，p和u均出现间断面，式(3)和(4)不再成立。此时，对(6)作如下考虑：

$\begin{split} \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}} =& \left[ {\frac{\partial }{{\partial \rho }}\left( {\frac{1}{\varGamma }} \right)} \right]\frac{{\partial \rho }}{{\partial t}} + u\left[ {\frac{\partial }{{\partial \rho }}\left( {\frac{1}{\varGamma }} \right)} \right]\frac{{\partial \rho }}{{\partial x}} + \rho \left[ {\frac{\partial }{{\partial \rho }}\left( {\frac{1}{\varGamma }} \right)} \right]\frac{{\partial u}}{{\partial x}} \\ {\rm{ = }}&\left[ {\frac{\partial }{{\partial \rho }}\left( {\frac{1}{\varGamma }} \right)} \right]\left(\frac{{\partial \rho }}{{\partial t}} + \frac{{\partial \rho u}}{{\partial x}}\right) {\rm{ = }}0 \end{split}$

(9)

因此，在左右波面处虽然p和u均间断，但由于在同类介质中1/Γ等参数的导数保持连续，因此，式(6)依然成立。同理，式(7)和(8)在左右波面处也成立。

1.3 输送方程

对于多介质流场，必须利用输送方程来捕捉界面的运动。流场中，某种介质的界面随时间变化的运动方程，称为输送方程。一般利用体积分数z_i来建立介质的输送方程：

$\frac{{\partial {{\textit z}_i}}}{{\partial t}} + u\frac{{\partial {{\textit z}_i}}}{{\partial x}} = 0\;\;\;\;\;\;\;i = 1,2, \cdots$

(10)

式中：i为介质的种类序号。非守恒变量1/Γ、p_ref/Γ、ρe_ref均为独立变量，但它们对的ρ偏导数可表达成关于体积分数z_i的加权叠加。这样，综合式(1)、(6)～(8)和(10)，可完成直角坐标系下的多介质流场的求解。

2. 球坐标系下的计算模型

球坐标系下，则需要考虑更多的问题。

第一，在采用直角坐标系下计算平面波问题时，p与u在界面处维持不变且偏导数∂p/∂x与∂u/∂x基本为0。该数值特性一直从界面处向左右延伸，直到波面处才形成间断，如图1所示。但是采用球坐标系处理球面波问题时，当冲击波沿着半径r向外扩散时，传播的形式则不再是平面冲击波，而是球面冲击波。此时，虽然p与u依然在界面处连续，但是球面波的扩散过程中p与u会随着半径增加而发生衰减，并非保持不变。因此，球面冲击波在界面处的偏导数∂p/∂r与∂u/∂r不再是0，而是一个稳定的负数。这种情况下，式(3)～(4)非守恒性方程的建立的理论基础并不十分准确。

图 1 平面冲击波界面处p的数值特征

Figure 1. The numerical property of p at the interfacefor plane shock wave

下载: 全尺寸图片幻灯片

第二，球坐标系下的内外介质产生相互作用时，内部的稀疏波(或激波)汇聚到一起后又反向朝着外部扩散，原点处便形成了一个奇点。奇点处的稀疏波形成了独特的反射机制。因此奇点的物理参数变化需仔细考虑。

2.1 球坐标系下控制方程

这里将球坐标系下的控制方程表示成守恒变量U与通量F的关系为^[13-14]:

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

(11)

其中： ${U} = \left[ \begin{array}{l}\rho \\\rho u\\\rho E\end{array} \right],\;{F} = \left[ \begin{array}{l}\rho u\\\rho {u^2} + p\\\rho Eu + pu\end{array} \right],\;{S} = - \displaystyle\frac{\alpha }{r}\left[ \begin{array}{l}\rho u\\\rho {u^2} + p\\\rho Eu + pu\end{array} \right]$

式中：系数α=2。式(11)可以改写成：

$\frac{{\partial \tilde{ U}}}{{\partial t}} + \frac{{\partial \tilde{ F}}}{{\partial r}} = {S}(\tilde{ U})$

(12)

其中： $\tilde{ U} = {r^2}{U},\;\tilde{ F} = {r^2}{F},\;{S}(\tilde{ U}) = {[0,\;2p/r, \;0]^{\rm T}}$

此时，将非守恒形式的方程组添加到式(11)中，即可得到球坐标系下的控制方程。

2.2 热力学参数修正

需要注意的是，非守恒形式的方程式(6)～(8)的推导，是以界面处于平衡状态为前提的。在球坐标系下，由于冲击波是沿着球面向四周扩散。随着扩散范围的增加，几乎所有的冲击波参数都会随着半径增加而降低。此时，p与u在界面处将呈现出持续的减小趋势。这样，式(6)～(8)无法直接使用到式(12)中。在直角坐标系和球坐标系下，界面间断的关系分别为

直角坐标系 $\displaystyle\frac{{\partial p}}{{\partial x}},\;\displaystyle\frac{{\partial u}}{{\partial x}} \to {\rm{0,}}\quad\;\;\;\displaystyle\frac{{\partial \rho }}{{\partial x}}, \;\displaystyle\frac{{\partial (\rho e)}}{{\partial x}},\;\displaystyle\frac{\partial }{{\partial x}}\left( {\displaystyle\frac{1}{\varGamma }} \right),\;\displaystyle\frac{\partial }{{\partial x}}\left( {\displaystyle\frac{{{p_{{\rm{ref}}}}}}{\varGamma }} \right),\;\displaystyle\frac{\partial }{{\partial x}}\left( {\rho {e_{{\rm{ref}}}}} \right) \to \infty$

球坐标系 $\displaystyle\frac{{\partial p}}{{\partial r}}{\rm{\text{＜} 0}},\; \displaystyle\frac{{\partial u}}{{\partial r}}{\rm{ \text{＜} 0,}} \quad\;\;\;\displaystyle\frac{{\partial \rho }}{{\partial r}},\displaystyle\frac{{\partial (\rho e)}}{{\partial r}},\displaystyle\frac{\partial }{{\partial r}}\left( {\displaystyle\frac{1}{\varGamma }} \right),\displaystyle\frac{\partial }{{\partial r}}\left( {\displaystyle\frac{{{p_{{\rm{ref}}}}}}{\varGamma }} \right),\displaystyle\frac{\partial }{{\partial r}}\left( {\rho {e_{{\rm{ref}}}}} \right) \to \infty$

考虑到三个非守恒变量1/Γ、p_ref/Γ、ρe_ref均为密度ρ的函数，这里每隔一个时间步对它们做修正:

$\frac{1}{\varGamma }{\rm{ = }}\frac{1}{{\varGamma ({\rho _i})}},\;\;\; \frac{{{p_{{\rm{ref}}}}}}{\varGamma } = \frac{{{p_{{\rm{ref}}}}({\rho _i})}}{{\varGamma ({\rho _i})}},\;\;\;\rho {e_{{\rm{ref}}}} = {\rho _i}{e_{{\rm{ref}}}}({\rho _i})\;\;\;\;\;\;\frac{{{\rho _i}}}{\rho } \text{＞} 0.99\;\;\;\;\;i = 1,2, \cdots$

(13)

式中：ρ_i为介质i的密度。而的其他参数如ρ、p、u、E的计算，仍然借助当前时间步下的Γ、p_ref、e_ref。而式(6)～(8)则保留在系统方程中，并且每隔一个时间步，非守恒变量1/Γ、p_ref/Γ、ρe_ref都会得到修正。当时间步足够密时，非守恒形式在球坐标系下带来的误差可以控制在一定的范围内。

2.3 输送方程变换

由于每隔一个时间步，都会需要利用ρ_i进行修正。这里用质量分数代替体积分数来构造输送方程。利用式(12)中的质量守恒关系可得：

$\frac{{\partial ({r^2}\rho {y_i})}}{{\partial t}} + \frac{{\partial ({r^2}\rho {y_i}u)}}{{\partial r}} = 0\;\;\;\;\;\;i = 1,2, \cdots$

(14)

根据式(14)得到的质量分数y_i可以很容易得到各介质的密度ρ_i (i=1, 2, ···)。而非守恒方程(6)～(8)中的偏导数，则转换成与质量分数相关的加权叠加：

$\begin{aligned} &\phi {\rm{ = }}\frac{\partial }{{\partial \rho }}\left( {\frac{1}{\varGamma }} \right){\rm{ = }}{y_1}{\left. {\left( { - \frac{{{{\varGamma }_1'}}}{{\varGamma _1^2}}} \right)} \right|_{\rho {y_1}}} + {y_2}{\left. {\left( { - \frac{{{{\varGamma }_2'}}}{{\varGamma _2^2}}} \right)} \right|_{\rho {y_2}}} + \cdots \\ &\varphi {\rm{ = }}\frac{\partial }{{\partial \rho }}\left( {\frac{{{p_{{\rm{ref}}}}}}{\varGamma }} \right){\rm{ = }}{y_1}{\left. {\left( {\frac{{{\varGamma _1}{{p}_{{\rm{ref,}}1}'} - {p_{{\rm{ref,1}}}}{{\varGamma }_1'}}}{{\varGamma _1^2}}} \right)} \right|_{\rho {y_1}}} + {y_2}{\left. {\left( {\frac{{{\varGamma _2}{{p}_{{\rm{ref}},2}'} - {p_{{\rm{ref}},2}}{{\varGamma }_2}'}}{{\varGamma _2^2}}} \right)} \right|_{\rho {y_2}}} + \cdots \\ &\psi {\rm{ = }}\frac{\partial }{{\partial \rho }}\left( {\rho {e_{{\rm{ref}}}}} \right){\rm{ = }}{y_1}{\left. {\left( {{\rho _1}{{e}_{{\rm{ref,}}1}'} + {e_{{\rm{ref,}}1}}} \right)} \right|_{\rho {y_1}}} + {y_2}{\left. {\left( {{\rho _2}{{e}_{{\rm{ref,}}2}'} + {e_{{\rm{ref,}}2}}} \right)} \right|_{\rho {y_2}}} + \cdots \end{aligned}$

(15)

式中：ϕ、φ、ψ为非守恒变量1/Γ、p_ref/Γ、ρe_ref偏导数，下角标1、2代表流场中介质的序号。

相比式(10)给出的基于体积分数的输送方程，基于质量分数的输送方程能考虑到每一种介质的密度ρ_i，从而得到更为合理的偏导数ϕ、φ、ψ。

2.4 奇点处理

对于球坐标系来说，奇点是一个无法避免的问题。参考文献[15]中对奇点的处理方式，奇点处质点的运动速度为0，即u=0，但是奇点处的质量和能量守恒关系依然成立：

$\begin{aligned} \frac{{\partial \rho }}{{\partial t}} + 3\frac{{\partial (\rho u)}}{{\partial r}} = 0, \quad\;\quad\; \frac{{\partial \rho E}}{{\partial t}} + 3\frac{{\partial (\rho Eu + pu)}}{{\partial r}} = 0 \end{aligned}$

(16)

展开后，将u=0代入，有：

$\begin{aligned} \frac{{\partial \rho }}{{\partial t}} + 3\rho \frac{{\partial u}}{{\partial r}} = 0, \quad\;\quad\; \frac{{\partial \rho E}}{{\partial t}} + 3(\rho E + p)\frac{{\partial u}}{{\partial r}} = 0 \end{aligned}$

(17)

由于在奇点处主要涉及冲击波的汇聚和发散，不涉及界面的运动，因此，可不考虑输送方程。参数Γ、p_ref、e_ref由密度ρ直接计算得到。

文献[15]中，对偏导数∂u/∂r采用了如下的计算方式：

$\frac{{\partial u}}{{\partial r}} = 2\frac{{{u_1}}}{{\Delta r}} - \frac{{{u_2}}}{{2\Delta r}}$

(18)

式中：下角标1和2代表网格点序号。但考虑到本文的方程体系非守恒性较强，将u₁和u₂都替换为奇点处的速度0，这样偏导数∂u/∂r=0。再根据式(17)计算得到奇点处的其他参数。

2.5 差分格式

在有限体积差分下，第j个网格点的变量U_j，可按如下格式计算：

$\frac{{\tilde { U}_j^{(n + 1)} - \tilde { U}_j^{(n)}}}{{\Delta t}} + \frac{{{{\tilde { F}}_{j + 1/2}} - {{\tilde { F}}_{j - 1/2}}}}{{\Delta r}} = {\tilde { S}_j}$

(19)

式中：Δt和Δr分别为时间和空间步长，上标(n)代表第n步时间序号。对于物理量 ${\tilde { U}_j}$ ，第j和j+1个网格点之间的通量 ${\tilde { F}_{j + 1/2}}$ 和源项 ${\tilde { S}_j}$ 分别为：

${{\tilde{ U}}_j} = r_j^2\left[ \begin{array}{l} {\rho _j}\\ {\rho _j}{u_j}\\ {\rho _j}{E_j}\\ 1/{\varGamma _j}\\ {p_{ref\;j}}/{\varGamma _j}\\ {\rho _j}{e_{{\rm{ref}}\;j}}\\ \rho {y_i}_j \end{array} \right]\;\;\;\; {{\tilde{ F}}_{j + 1/2}} = r_j^2\left[ \begin{array}{l} {\rho _{j + 1/2}}{u_{j + 1/2}}\\ {\rho _{{\rm{j + }}1/2}}u_{{\rm{j + }}1/2}^2 + {p_{{\rm{j + }}1/2}}\\ ({\rho _{j + 1/2}}{E_{j + 1/2}} + {p_{j + 1/2}}){u_{j + 1/2}}\\ {u_j} \cdot {(1/\varGamma )_{j + 1/2}} + {\rho _j}{\phi _j} \cdot {u_{j + 1/2}}\\ {u_j} \cdot {({p_{{\rm{ref}}}}/\varGamma )_{j + 1/2}} + {\rho _j}{\varphi _j} \cdot {u_{j + 1/2}}\\ {u_j} \cdot {({\rho}{e_{{\rm{ref}}}})_{j + 1/2}} + {\rho _j}{\psi _j} \cdot {u_{j{\rm{ + }}1/2}}\\ {(\rho {y_i})_{j + 1/2}} \cdot {u_{{\rm{j + }}1/2}} \end{array} \right]\;\;\;\; {{\tilde{ S}}_j} = 2{r_j}\left[ \begin{array}{l} 0\\ {p_j}\\ 0\\ 0\\ 0\\ 0\\ 0 \end{array} \right]$

(20)

式中：角标j+1/2表示通量中的参数，角标j表示网格点的参数；在球面冲击波中，通量 ${{\tilde{ F}}_{j + 1/2}}$ 表示沿半径r方向通过j点的物理变量。

2.6 时间空间步长

时间步长Δt，满足收敛的条件为：

$\Delta t = {C_{{\rm{CFL}}}}\frac{{\min (\Delta r)}}{{\max (\left| {{u_j}} \right| + {c_j})}}$

(21)

式中：C_CFL是为了保证迎风型计算格式收敛而设置的一个0～1之间的系数，满足收敛的条件为C_CFL＜1, 并且当C_CFL取得越小时间步长越小，计算越精细, 这里取C_CFL=0.4； u_j和c_j表示第j个网格点的质点速度和声速。声速c_j的表达式为：

${{c}_{j}}=\sqrt{{\left.\left( {{H}_{j}}-\frac{1}{2}{{u}_{j}}{}^{2}-{{p}_{j}}{{\phi }_{j}}+{{\varphi }_{j}}-{{\psi }_{j}} \right)\right/}{\left( \frac{1}{{{\varGamma }_{j}}} \right)}\;}$

(22)

式中：H=E+p/ρ。

3. 计算结果

3.1 无因次参数的气水作用算例

这里考虑一个无因次的气水作用问题^[16]。初始时刻，中心气团的物理状态为：p=8.3×10³，ρ=1.27，u=0，外部水的物理状态为：p=1.0，ρ=1.27，u=0。气体和水的状态方程分别为

$\begin{array}{c} p = (\gamma - 1)\rho e,\quad\;\;\gamma = 1.4;\\ p = (\gamma - 1)\rho e - \gamma B,\quad\;\;\gamma = 7.0,\quad\;\;B = 3\;000 \end{array}$

式中：γ和B为描述水和气体在低压缩度时的热力学参数。气团初始半径为R₀。由于中心气团的压力和密度高于周边水的。因此，高压的气团会在水中形成冲击波并推动界面向外部扩散。对冲击波波面和介质界面的运动位置的变化情况进行记录（波面为水中压力的最大值位置，界面位置为从外向内气体的质量分数刚超过0.5的位置），并绘制了冲击波和界面随时间变化的曲线，如图2所示。从图2可以看出，本文所得到的冲击波和界面的运动位置与Liu等^[16]和Flores等^[17]的结果都比较接近。三个数据的计算结果，在气水作用的初期差别不大。但是随着冲击波和界面向外扩散，本文的计算结果与Flores等^[17]的计算结果相比，误差略偏大，但仍然在一个合理范围内。

图 2 冲击波与界面位置随时间的变化

Figure 2. Time evolutions of position of shock wave and interface

下载: 全尺寸图片幻灯片

图3为冲击波到达3倍气团初始半径时候的压力和速度曲线。此时，向内收缩的稀疏波已达到中心原点处。由于气团已扩散开，中心处的压力开始逐渐下降，但中心处的速度仍然为0。速度曲线中，由中心向外的第一个拐点为稀疏波的端点，第二个拐点则为界面所在位置。另外，从图3中可以看出，处于数值奇点的原点处，数值稳定性较好。证明本文所用的奇点处理方式能效果良好。

图 3 冲击波到达3倍R₀时，压力与速度的分布

Figure 3. The distributions of pressure and velocity when shock wave reaches 3R₀

下载: 全尺寸图片幻灯片

3.2 圆形药球的水下爆炸算例

假设一TNT的实心药球质量为50 kg，将空间步长取为炸药球半径的1/100。利用改进后的多介质混合模型计算球坐标系下的冲击波运动情况。

炸药的状态方程为JWL状态方程^[3]：

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

(23)

式中：A、B、R₁、R₂、ω均为常数，由圆筒试验标定得到^[3]；θ=ρ/ρ₀，ρ₀为物质的初始密度。水的多项式状态方程为^[18]：

$p = \left\{ {\begin{array}{*{20}{l}} {{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}&\quad{\mu \text{＞} 0}\\ {{a_1}\mu + \left( {{b_0} + {b_1}\mu } \right){\rho _0}E}&\quad{\mu \text{＜} 0} \end{array}} \right.$

(24)

式中：μ=θ−1；a₁～a₃, b₀～b₃为展开后的多项式系数。状态方程(23)～(24)的参数取值见表1^[19]。

表 1 状态方程参数

Table 1. Coefficients of EOS

材料	ρ₀/(kg·m⁻³)	A/GPa	B/GPa	R₁	R₂	ω
TNT炸药	1 630	371.2	3.21	4.15	0.95	0.35

下载: 导出CSV

| 显示表格

材料	ρ₀/(kg·m⁻³)	a₁/GPa	a₂/GPa	a₃/GPa	b₀	b₁	b₂	b₃
水	1 000	2.19	9.224	8.767	0.394	1.393 7	0	0

下载: 导出CSV

| 显示表格

本算例中，将计算范围扩大至20倍以上的药球半径。这样，可以得到冲击波运动到不同位置处的压力，如图4所示。从图4可以看出，由于在偏导数ϕ、φ、ψ的处理方式差异，两种计算模型得到的压力分布有差别。但都具备较好的数值稳定性。当冲击波运动到较远位置时，会出现多个压力峰值。

图 4 冲击波到达不同位置处的压力曲线

Figure 4. Pressure curves at the time instants when the shock wave reaches different positions

下载: 全尺寸图片幻灯片

在实际工程中，最受重视的是第一次的峰值。图5为冲击波峰值压力随运动距离变化的曲线。作为对比，将给初始的体积分数模型^[2]和改进后的质量分数模型同Zamyshlyayev的经验公式^[20]曲线作了对比。对比发现，改进后的质量分数模型可以得到更为准确的计算结果。在压力随运动距离发生急剧下降的过程中，利用体积分数得到的计算结果，得到的结果比经验公式的结果略微偏大。

图 5 冲击波在不同位置处的压力峰值变化

Figure 5. Peak value variation of pressure of the shock wave at different positions

下载: 全尺寸图片幻灯片

图6为利用体积分数与质量分数计算模型分别得到的压力时程曲线。从不同的计算结果中可以看出，利用改进后的质量分数模型，得到的计算结果好于体积分数模型。各个半径处的压力时程，质量分数模型均吻合得比较好。由于Zamyshlyayev的经验公式仅考虑冲击波压力最大峰值未考虑压力的二次峰值，本文也仅模拟了冲击波压力的最大峰值。

图 6 不同计算模型下压力随时间变化曲线

Figure 6. Time evolutions of pressure for different numerical models

下载: 全尺寸图片幻灯片

4. 结　语

针对球坐标下多介质界面捕捉困难的问题，本文对原始的Mie-Grüneisen多介质混合物模型作了热力学参数修正、输送方程变换、特殊奇点处理等多项数值改进。改进后的基于质量分数的模型，可以在常规状态方程下准确捕捉界面和冲击波的运动，且可以在适应JWL等形式复杂的状态方程。而适当的奇点处理也保证了数值计算在奇点处的稳定性。另外，数值算例证明，相比原始的体积分数模型，改进后的质量分数模型可以得到较为准确的计算结果。

图 1 平面冲击波界面处p的数值特征

Figure 1. The numerical property of p at the interfacefor plane shock wave

下载: 全尺寸图片幻灯片

图 2 冲击波与界面位置随时间的变化

Figure 2. Time evolutions of position of shock wave and interface

下载: 全尺寸图片幻灯片

图 3 冲击波到达3倍R₀时，压力与速度的分布

Figure 3. The distributions of pressure and velocity when shock wave reaches 3R₀

下载: 全尺寸图片幻灯片

图 4 冲击波到达不同位置处的压力曲线

Figure 4. Pressure curves at the time instants when the shock wave reaches different positions

下载: 全尺寸图片幻灯片

图 5 冲击波在不同位置处的压力峰值变化

Figure 5. Peak value variation of pressure of the shock wave at different positions

下载: 全尺寸图片幻灯片

图 6 不同计算模型下压力随时间变化曲线

Figure 6. Time evolutions of pressure for different numerical models

下载: 全尺寸图片幻灯片

表 1 状态方程参数

Table 1. Coefficients of EOS

材料	ρ₀/(kg·m⁻³)	A/GPa	B/GPa	R₁	R₂	ω
TNT炸药	1 630	371.2	3.21	4.15	0.95	0.35

下载: 导出CSV

材料	ρ₀/(kg·m⁻³)	a₁/GPa	a₂/GPa	a₃/GPa	b₀	b₁	b₂	b₃
水	1 000	2.19	9.224	8.767	0.394	1.393 7	0	0

下载: 导出CSV

参考文献(20)

[1]	LIUT G, KHOOB C, YEOK S. Ghost fluid method for strong shock impacting on material interface [J]. Journal of Computational Physics, 2003, 190(2): 651–681. DOI: 10.1016/S0021-9991(03)00301-2.
[2]	SHYUE K M. A fluid-mixture type algorithm for compressible multicomponent flow with Mie-Grüneisen equation of state [J]. Journal of Computational Physics, 2001, 171(2): 678–707. DOI: 10.1006/jcph.2001.6801.
[3]	张宝銔, 张庆民, 黄风雷. 爆轰物理学 [M]. 北京: 北京理工大学出版社, 2001: 160, 377-383.
[4]	ABGRALL R. How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach [J]. Journal of Computational Physics, 1996, 125(1): 150–160. DOI: 10.1006/jcph.1996.0085.
[5]	SAUREL R, ABGRALL R. A simple method for compressible multifuid flows [J]. SIAM Journal on Scientific Computing, 1999, 21(3): 1115–1145. DOI: 10.1137/S1064827597323749.
[6]	ABGRALL R, KARNI S. Computations of compressible multifluids [J]. Journal of Computational Physics, 2001, 169(2): 594–623. DOI: 10.1006/jcph.2000.6685.
[7]	SAUREL R, ABGRALL R. A multiphase Godunov method for compressible multifluid and multiphase flows [J]. Journal of Computational Physics, 1999, 150(2): 425–467. DOI: 10.1006/jcph.1999.6187.
[8]	柏劲松, 陈森华, 李平. 多介质流体非守恒律欧拉方程组的数值计算方法 [J]. 爆炸与冲击, 2001, 21(4): 265–271. BO Jinsong, CHEN Senhua, LI Ping. Numerical methods of multicomponent flows of non-conservative Euler equations [J]. Explosion and Shock Waves, 2001, 21(4): 265–271.
[9]	柏劲松, 陈森华, 李平, 等. 多介质可压缩流体动力学界面捕捉方法 [J]. 爆炸与冲击, 2004, 24(1): 37–43. BAI Jingsong, CHEN Senhua, LI Ping. Interface capturing method for compressible multi-fluid dynamics [J]. Explosion and Shock Waves, 2004, 24(1): 37–43.
[10]	梁姗, 刘伟, 袁礼. 七方程可压缩多相流模型的HLLC格式及应用 [J]. 力学学报, 2012, 44(5): 884–895. DOI: 10.6052/0459-1879-12-022. LIANG Shan, LIU Wei, YUAN Li. An HLLC scheme for the seven-equation multiphase model and its application to compressible multicomponent flow [J]. Chinese Journal of Theoretical and Applied Mechanics, 2012, 44(5): 884–895. DOI: 10.6052/0459-1879-12-022.
[11]	LIANG S, LIU W, YUAN L. Solving seven-equation model for compressible two-phase flow using multiple GPUs [J]. Computers and Fluids, 2014, 99(7): 156–171. DOI: 10.1016/j.compfluid.2014.04.021.
[12]	刘娜, 陈艺冰. 多介质流体力学计算的谱体积方法 [J]. 爆炸与冲击, 2017, 37(1): 114–119. DOI: 10.11883/1001-1455(2017)01-0114-06. LIU Na, CHEN Yibing. High order spectral volume method for multi-component flows [J]. Explosion and Shock Waves, 2017, 37(1): 114–119. DOI: 10.11883/1001-1455(2017)01-0114-06.
[13]	ELEUTERIO F T. Riemann solvers and numerical methods for fluid dynamics [M]. 3rd ed. Berlin Heidelberg: Springer, 2009: 1-40. DOI: 10.1007/978-3-540-49834-6.
[14]	师华强, 宗智, 贾敬蓓. 水下爆炸冲击波的近场特性 [J]. 爆炸与冲击, 2009, 29(2): 125–130. DOI: 10.11883/1001-1455(2009)02-0125-06. SHI Huaqiang, ZONG Zhi, JIA Jingbei. Short-range characters of underwater blast waves [J]. Explosion and Shock Waves, 2009, 29(2): 125–130. DOI: 10.11883/1001-1455(2009)02-0125-06.
[15]	LIU T G, KHOO B C, YEO K S. The numerical simulations of explosion and implosion in air: use of a modified Harten's TVD scheme [J]. International Journal for Numerical Methods in Fluids, 1999, 31(4): 661−680. DOI: 10.1002/(SICI)1097-0363(19991030)31:4<661::AID-FLD866>3.0.CO;2-G.
[16]	LIUT G, KHOOB C, YEOK S. The simulation of compressible multi-medium flow: I: A new methodology with test application to 1D gas-gas and gas-water cases [J]. Computers and Fluids, 2001, 30(3): 291–314. DOI: 10.1016/S0045-7930(00)00022-0.
[17]	FLORES J, HOLT M. Glimm’s method applied to underwater explosion [J]. Journal of Computational Physics, 1981, 44(2): 377–387. DOI: 10.1016/0021-9991(81)90058-9.
[18]	CHIESUM J E, SHIN Y S. Explosion gas bubbles near simple boundaries [J]. Shock and Vibration, 1997, 4(1): 11–25. doi: 10.1155/1997/615415
[19]	LIU G R, LIU M B. 光滑粒子动力学:一种无网格粒子法[M]. 韩旭, 等译. 长沙: 湖南大学出版社, 2005.
[20]	ZAMYSHLYAYEV B V, YAKOVLEV Y S. Dynamic loads in underwater explosion: AD-757183[R]. USA: Naval intelligence Support Center, 1973.

施引文献

期刊类型引用(4)

1.	汪泉，胡程，谢守冬，李孝臣，李志敏，涂唱畅，朱群龙，杨锐. 负压条件对T2/Q235爆炸焊接复合板界面的影响. 火炸药学报. 2024(01): 64-71 . 百度学术
2.	李晓杰，易生泰，闫鸿浩，王小红，李科斌. 两种介质中爆轰波和冲击波连续测量技术研究. 爆破器材. 2017(04): 55-58 . 百度学术
3.	周楠，蒋敬，李浩，唐松泽，陶秋辰. 爆炸焊接不锈钢/铝复合板界面形态研究. 重庆理工大学学报(自然科学). 2017(07): 91-95+101 . 百度学术
4.	王亚朋，李晓杰，张程娇，王小红，闫鸿浩. 柱状炸药水下冲击波连续测量研究. 工程爆破. 2015(04): 17-19+57 . 百度学术

其他类型引用(4)

资源附件(0)

访问统计

图(6) / 表(2)

计量

文章访问数: 5604
HTML全文浏览量: 1513
PDF下载量: 41
被引次数: 8

1. 直角坐标系下的计算模型
1.1 界面处热力学关系
1.2 波面处热力学关系
1.3 输送方程
2. 球坐标系下的计算模型
2.1 球坐标系下控制方程
2.2 热力学参数修正
2.3 输送方程变换
2.4 奇点处理
2.5 差分格式
2.6 时间空间步长
3. 计算结果
3.1 无因次参数的气水作用算例
3.2 圆形药球的水下爆炸算例
4. 结　语

1. 直角坐标系下的计算模型
1.1 界面处热力学关系
1.2 波面处热力学关系
1.3 输送方程
2. 球坐标系下的计算模型
2.1 球坐标系下控制方程
2.2 热力学参数修正
2.3 输送方程变换
2.4 奇点处理
2.5 差分格式
2.6 时间空间步长
3. 计算结果
3.1 无因次参数的气水作用算例
3.2 圆形药球的水下爆炸算例
4. 结　语

参考文献(20)

施引文献

资源附件(0)

访问统计

球坐标系下多介质混合物模型的数值模拟

doi: 10.11883/bzycj-2018-0075

作者简介: 吴宗铎（1984－ ），男，博士，讲师，wuzongduo0@aliyun.com

通讯作者: 赵 勇（1981－ ），男，博士，副教授，fluid@126.com

计量

出版历程

Numerical simulation about the multi-component mixture model under spherical coordinate system

1. 直角坐标系下的计算模型

1.1 界面处热力学关系

1.2 波面处热力学关系

1.3 输送方程

2. 球坐标系下的计算模型

2.1 球坐标系下控制方程

2.2 热力学参数修正

2.3 输送方程变换

2.4 奇点处理

2.5 差分格式

2.6 时间空间步长

3. 计算结果

3.1 无因次参数的气水作用算例

3.2 圆形药球的水下爆炸算例

4. 结 语

期刊类型引用(4)

其他类型引用(4)

计量

出版历程

目录

1. 直角坐标系下的计算模型

1.1 界面处热力学关系

1.2 波面处热力学关系

1.3 输送方程

2. 球坐标系下的计算模型

2.1 球坐标系下控制方程

2.2 热力学参数修正

2.3 输送方程变换

2.4 奇点处理

2.5 差分格式

2.6 时间空间步长

3. 计算结果

3.1 无因次参数的气水作用算例

3.2 圆形药球的水下爆炸算例

4. 结 语

作者简介:
吴宗铎（1984－　），男，博士，讲师，wuzongduo0@aliyun.com

通讯作者:
赵　勇（1981－　），男，博士，副教授，fluid@126.com

4. 结　语

4. 结　语