水气二相流光滑粒子流体动力学方法研究进展

摘要

从状态方程、连续方程、动量方程和其他数值技术等方面，对光滑粒子流体动力学（SPH）解决水气二相流问题的数值技术和最新进展进行了综述，并对各项技术进行了评价。最后，提出了未来的研究展望。

关键词

光滑粒子流体动力学（SPH）; 水气二相流; 数值技术

光滑粒子流体动力学（SPH）方法是一种强拉格朗日性质的粒子法，在处理大变形、二相流、自由液面等问题上有其独特优势。不同学者对SPH基本方法^［1-2］、SPH方法在海岸和海洋工程中的应用^［3］、复杂流体模拟技术^［4］等做了很好的综述。也有学者对SPH和其他粒子法（如MPS （moving particle semi-implicit）‍，ISPH（incompressible smoothed particle hydrodynamics）等）的表现做了比较^［5-6］。二相流是SPH方法的重要应用方向，又包括水气二相流和水沙二相流等。单相流SPH方法在自由液面处压强为零，不考虑空气的影响，但在存在气泡、风、柔性边界等对流体产生显著影响的问题中，空气的动力学影响不能忽略。因此，水气二相流作为SPH的一个重要研究分支被广泛讨论。Colagrossi等^［7］第一次将SPH方法应用于水气二相流中，随后国内外学者在水气二相流SPH方面的研究取得了一系列进展，如含气静水水箱^［8］、水中气泡上升^［9-12］、两相溃坝^［7，13-15］、两相振荡液滴^［16-18］、两相驻波^［17-18］、两相液体晃荡^［14-15，19］等。

水气二相流SPH方法所要解决的主要问题是界面间物理参数的不连续引起的计算困难和数值不稳定，各国学者通过对水气交界面处的控制方程进行修正、添加界面力（表面张力和界面斥力等）、添加额外耗散项、采用额外的数值技术等方式来解决水气界面稳定性问题。本研究对水气二相流SPH方法涉及的数值方法进行综述，首先简述了单相流SPH方法的控制方程，然后详述了水气二相流SPH方法对状态方程、连续方程、动量方程的修正处理，最后讨论了背景压力、压力梯度、表面张力、界面斥力以及其他数值技术。兼顾了对水气二相流SPH数值实现方法的系统介绍及应用效果评价，以期为相关研究者提供详尽的数值方法选择依据。在此基础上，对未来水气二相流SPH方法的发展方向进行了展望。

1 单相流SPH方法控制方程

单相流是水气二相流研究的基础，首先对单相流SPH方法的控制方程进行简述。单相流SPH方法的控制方程通常采用拉格朗日形式的纳维斯托克斯（Navier-Stokes）方程，形式^［1］如下所示：

$$\frac{\mathrm{d}\rho }{\mathrm{d}t}=-\rho \nabla \cdot \mathit{u}$$

（1）

$$\frac{\mathrm{d}\mathit{u}}{\mathrm{d}t}=-\frac{\mathrm{1}}{\rho }\nabla p+\mathit{\Theta }-\mathit{g}$$

（2）

$$\frac{\mathrm{d}\mathit{r}}{\mathrm{d}t}=\mathit{u}$$

（3）

$$p=f\left(\frac{\rho }{{\rho }_{\mathrm{0}}}\right)$$

（4）

式中：$\rho$为密度；${\rho }_{\mathrm{0}}$为基准密度（压力为零时的密度）；t为时间；u为速度；p为压强；$\mathit{\Theta }$为黏度项；g为重力加速度；r为位移。在流体弱可压缩情况下，状态方程（4）指压力为ρ/ρ₀的函数，通常由以下2种形式来关联压力和密度^［20-21］：

$$p=B\left({\left(\frac{\rho }{{\rho }_{\mathrm{0}}}\right)}^{\gamma }-\mathrm{1}\right)$$

（5）

$$p=c_{\mathrm{0}}^{\mathrm{2}}\left(\rho -{\rho }_{\mathrm{0}}\right)$$

（6）

式中：$\gamma$为系数，对于流体$\gamma$取7，对于气体$\gamma$取1.4；$c_{\mathrm{0}}$为数值声速，弱可压缩假设要求密度变化范围在1%内，因此$c_{\mathrm{0}}$需满足大于流体最大速度的10倍以上；B为参数，B=$\left({\rho }_{\mathrm{0}}{c}_{\mathrm{0}}^{\mathrm{2}}\right)/\gamma$。

2 水气二相流SPH方法

SPH方法常处理的二相流问题包括水气二相流和水沙二相流。处理水沙二相流时通过增加新的粒子类型来表示泥沙粒子或水沙混合物，也有研究在水粒子中增加含沙量属性的选项，但总体来说水沙之间出现掺混，没有明显的二相界面；处理水气二相流时除了考虑两相之间密度差异巨大之外，还要处理水气界面处由密度不连续导致的压力振荡问题。针对水气二相流SPH问题的特点，相关学者提出了施加背景压力、压力梯度修正、表面张力、界面斥力等技术，本节将分别从状态方程、连续方程、动量方程以及其他技术的角度分别进行阐述。

2.1　状态方程

2.1.1　弱可压SPH方法状态方程

绝大多数单相流SPH方法不施加背景压力（部分不涉及自由面的单相流SPH方法引入了背景压力，如Marrone等^［22］在圆柱绕流模拟中添加了背景压力），但在水气二相流SPH方法中，为了防止流场中出现负压，通常需要在状态方程中添加一个固定的背景压力来避免张力不稳定，而张力不稳定是导致粒子聚集的原因之一^［23］。多数多相流SPH模拟仍然依赖背景压力来防止张力不稳定性的发生。这个背景压力通常是一个适中的数值，过小会导致流体质点的压力产生负压，而过大将导致流场额外的数值黏性和数值噪声。添加背景压力后的弱可压SPH状态方程形式为

$$p=B\left({\left(\frac{\rho }{{\rho }_{\mathrm{0}}}\right)}^{\gamma }-\mathrm{1}\right)+\chi$$

（7）

$$p=c_{\mathrm{0}}^{\mathrm{2}}(\rho -{\rho }_{\mathrm{0}})+\chi$$

（8）

式中：$\chi$为背景压力。文献［1］、［7］、［12］、［24］中采用式（7），为了保证界面处初始压力的平滑，气体的参数B和流体的相同，依此计算得到气体中的数值声速约为流体中的14倍，气体声速远大于流体声速，可保证数值稳定。文献［13］、［15］、［19］、［25］中采用式（8），数值声速$c_{\mathrm{0}}$需满足流体的弱可压缩性，气体和流体使用相同的数值。式（7）和式（8）对数值计算结果影响不大，但需注意两者气相数值声速取值不同。

背景压力$\chi$的取值方法主要分为2类。第1类取一个定值，如100 Pa^［13］、1 000 Pa^［26］（溃坝模型）等；第2类取最大静水压力的倍数，如$\mathrm{0.1}{\rho }_{\mathrm{l}}gH$^［8］（静水水箱模型）、$\mathrm{0.5}{\rho }_{\mathrm{l}}gH$^［27］（溃坝模型）、${\rho }_{\mathrm{l}}gH$^［8］（溃坝模型）等，其中${\rho }_{\mathrm{l}}$为流体密度，H为初始水深。随着流体雷诺数的增大，背景压力取值也应相应增大来避免大范围负压的产生。可以看到，后两者均为溃坝模型，但倍数取值有所不同，这是因为文献［27］采用了张力不稳定性控制技术（详见本节下文论述），允许存在部分负压。除了上述2类背景压力取值方法外，Ming等^［28］针对表面张力不可忽略的气泡上升和破裂案例，提出了一种与表面张力相关的背景压力形式，如下所示：

$$\chi =\mathrm{60}\frac{\left|{\rho }_{\mathrm{l}}-{\rho }_{\mathrm{g}}\right|}{{\rho }_{\mathrm{l}}+{\rho }_{\mathrm{g}}}\frac{\sigma }{R}$$

（9）

式中：${\rho }_{\mathrm{g}}$为气体密度；R为气泡的初始曲率半径；$\sigma$为表面张力系数。Ming等^［28］的方法能在邦德数较小时获得稳定二相界面，适用于表面张力占主导作用的小尺度气泡上升问题。

在水气二相流状态方程中添加背景压力是解决由负压导致的张力不稳定问题的最简单直接的方法，适用于各种工况，其难点在于合适背景压力数值的确认。事实上，负压是真实存在的，最近也出现了负压条件下张力不稳定性控制技术^［29］，可以在不添加背景压力的情况下防止张力不稳定。该方法最初是为了解决高雷诺数下自由液面流场中的数值空洞问题而提出的^［29］，后来扩展到多相流问题并适用于任意雷诺数下的流动^［27，30-32］。背景压力和张力不稳定性控制技术在多相流中可以同时应用，前者能够减少后者造成的动量不守恒^［27］，而后者通过保留负压使得前者的应用更加灵活。

除背景压力外，气相状态方程中额外增加的内聚力项^［7，33-35］也是保证清晰锐利的二相界面形状的数值方法，该方法基于分子作用力提出，主要适用于水中气泡运动案例，状态方程可以写成

$$p=B\left({\left(\frac{\rho }{{\rho }_{\mathrm{0}}}\right)}^{\gamma }-\mathrm{1}\right)+\chi -\overline{a}{\rho }^{\mathrm{2}}$$

（10）

式中：$\overline{a}$为内聚力系数。这种处理方式在水中气泡上升和两相溃坝问题中均有用到。值得注意的是，若在气相粒子状态方程中添加了内聚力项，则对应的气相动量方程中也要添加一个内聚力项${\mathit{F}}_{\mathrm{c}}$，计算式如下所示：

$${\mathit{F}}_{\mathrm{c}}=-\mathrm{2}\overline{a}{\rho }_{}^{\mathrm{2}}\sum _j{V}_j{\nabla }_i{W}_{ij}$$

（11）

式中：V为粒子体积；W为核函数；下标i为当前粒子编号；下标j为i粒子支持域内的其他粒子编号。内聚力项的增加在小尺度气液二相流案例中能够有效防止气体粒子发生穿透^［35］，维持气液二相的清晰界面。

2.1.2　强可压SPH方法状态方程

对于水下爆炸、气泡脉动、高速砰击入水等具有大幅体积变化的问题，需要考虑流体介质的真实压缩性。因此，不是使用人工声速，而是使用真实物理声速。在状态方程中考虑内能变化，压强可以写为密度和内能的函数，统一形式^［36］如下所示：

$$p=f(\rho ,e)$$

（12）

式中：e为内能。

刚性气体状态方程^［37-39］可作为强可压SPH方法的状态方程使用，形式为

$$p=(\beta -\mathrm{1})\rho e-\beta b$$

（13）

式中：β为定压比热容和定容比热容之比，对于理想气体，β数值上等同于绝热指数；b为与流体性质相关的参数，理想气体b为零时，式（13）退化为理想气体状态方程。在计算水气二相流问题时，气相使用理想气体状态方程，液相的β和b可采用对应流体的取值^［38，30］。宗智等^［40］应用式（13）模拟了二维真空爆炸和水下爆炸，取得了合理的模拟结果。

涉及到爆炸初期的计算，理想气体状态方程不再适用，Liu等^［41］采用了凝聚炸药爆轰产物的状态方程，即Jones-Wilkins-Lee（LWJ）方程^［42］。Liu等^［41］使用的气体和流体的状态方程如下所示：

$$p=A\left(\mathrm{1}-\frac{\omega \eta }{R_{\mathrm{1}}}\right){\mathrm{e}}^{-\frac{R_{\mathrm{1}}}{\eta }}+B_{\mathrm{1}}\left(\mathrm{1}-\frac{\omega \eta }{R_{\mathrm{2}}}\right){\mathrm{e}}^{-\frac{R_{\mathrm{2}}}{\eta }}+\omega \eta {\rho }_{\mathrm{0}}E$$

（14）

$$p={\rho }_{\mathrm{0}}{c}_{\mathrm{0}}^{\mathrm{2}}\mu +({\gamma }_{\mathrm{0}}+a_{\mathrm{0}}\mu ){\rm E}$$

（15）

式中：$\eta$为爆炸产物密度与原炸药初始密度之比；E为单位质量比内能；A、B₁、R₁、R₂、$\omega$均为拟合系数；c₀为初始声速；$\mu$为扰动前后的水密度之比减1；${\gamma }_{\mathrm{0}}$为系数，取0.5；a₀为体积修正系数。Liu等^［41］采用该方程后，模型较好地模拟了一维三硝基甲苯（TNT）平板爆炸、自由空间的水下爆炸、密闭空间内的水下爆炸等案例。Zhang等^［43］也采用LWJ方程作为状态方程模拟了水下爆炸的气液过程，取得了良好的模拟效果。

除式（13）、（14）之外，对于一些流体能量较小但可压缩性不能忽略的案例，如船底斜度角小时形成射流的模拟^［44-45］，正压存在于全计算域，仍可采用弱可压SPH状态方程以简化计算，并通过声速的合理取值得到合理的强可压SPH模拟结果^［30，46］。

2.2　连续方程

在SPH方法中使用连续方程求解密度，针对水气二相流等大密度比多相流SPH问题，Hu等^［47］提出了一种密度表示方法，如下所示：

$${\rho }_i=m_i\sum _j{W}_{ij}$$

（16）

式中：$m$为粒子质量。这种形式符合质量守恒，类似于常见的SPH密度表示方法${\rho }_i=\sum _j{m}_{j}Wij$，但是在式（16）中，相邻粒子的质量对选定粒子的密度没有贡献，能较好地处理界面密度不连续问题。

水气二相流SPH方法对于数值稳定性有较高的要求，一些用于单相流SPH方法的数值稳定方法在水气二相流SPH方法中得到广泛应用。δ-SPH方法是连续方程中应用最广泛的数值稳定方法，在水气二相流中也得到普遍应用，因此本节对δ-SPH方法加以阐述。

Molteni等^［48］在2009年提出了带有密度耗散的δ-SPH方法，在弱可压缩SPH方法的基础上，在连续方程中加入额外的密度扩散项减少非物理的压力波动。在Antuono等^［49］、Marrone等^［50］的贡献下， δ-SPH方法逐步得到建立和完善。事实证明该方法能够显著减少压力场中的非物理波动^［51-52］。δ-SPH方法可以写为

$$\frac{D{\rho }_i}{Dt}=-m_i\sum _j({\mathit{u}}_j-{\mathit{u}}_i)\cdot {\nabla }_i{W}_{ij}+\delta h{c}_{\mathrm{0}}\sum _j{D}_i{\nabla }_i{W}_{ij}{V}_j$$

（17）

$$\frac{({\mathit{r}}_j-{\mathit{r}}_i)}{({\mathit{r}}_j-{\mathit{r}}_i{)}^{\mathrm{2}}+{(\mathrm{0.01}h)}^{\mathrm{2}}}$$

（18）

（19）

$$L_i={\left(\sum _j({\mathit{r}}_j-{\mathit{r}}_i)\otimes {\nabla }_i{W}_{ij}{V}_j\right)}^{-\mathrm{1}}$$

（20）

式中：D为微分符号；$\delta$为密度扩散系数，通常取固定值0.1；h为光滑长度；为重整的密度梯度。式（17）右侧第二项即为密度耗散项。δ-SPH方法不影响流体的整体演化，仅作为压力场平滑的一个参数作用于局部压力场。Sun等^［27，32］将经典δ-SPH方法扩展到多相问题，对水流冲击弹性板的溃坝模型和楔形体入水（见图1）等问题进行了模拟，表明SPH模拟结果与实验现象吻合良好。水气二相流SPH模拟中，δ-SPH方法成功用于振荡液滴、驻波等的水气二相模拟^［18］。Hammani等^［8］又采用2种方法对多相流中的δ-SPH方法进行推导，从理论上进一步论证了将δ-SPH方法拓展到多相流的正确性。值得注意的是，密度扩散项仅应用于同相粒子中，异相界面处不考虑密度扩散项。

图1  不同时刻楔形体入水实验拍摄现象^［53］与多相流δ-SPH模拟结果^［32］对比

Fig.1  Comparison between experimental snapshots of wedge entry^[53] and results of multi-phase δ-SPH model^[32] at different moments

2.3　动量方程

水气二相流SPH方法的动量方程一般包含压力梯度项、黏度项、重力项、表面张力项，以及防止异相粒子相互穿透的界面斥力项（有时还会添加内聚力项，在2.1节中已有介绍，这里不再赘述），其中表面张力项及界面斥力项是二相流特有的。动量方程如下所示：

$$\frac{D\mathit{u}}{Dt}=\frac{\mathrm{1}}{\rho }(-\nabla p+\mathit{\Theta }+\rho \mathit{g}+{\mathit{F}}_{\mathbf{s}}+{\mathit{F}}_{\mathbf{i}\mathbf{n}\mathbf{t}\mathbf{e}\mathbf{r}\mathbf{f}\mathbf{a}\mathbf{c}\mathbf{e}})$$

（21）

式中：${\mathit{F}}_{\mathbf{s}}$为表面张力；${\mathit{F}}_{\mathbf{i}\mathbf{n}\mathbf{t}\mathbf{e}\mathbf{r}\mathbf{f}\mathbf{a}\mathbf{c}\mathbf{e}}$为界面斥力。本节对压力梯度项、表面张力项和界面斥力项分别进行论述。黏度项（包括人工黏度和物理黏度）的处理与单相流SPH方法一致，此处不再赘述，需要注意的是黏度项仅应用于同相粒子相互作用中。

2.3.1　压力梯度项

水气的密度比高达1 000，界面处密度不连续，在计算压力时极易产生较大的误差。Hu等^［47，54］利用粒子平均空间导数近似处理密度，相邻粒子对选定粒子只有体积贡献而没有密度贡献，从而较好地处理界面密度不连续问题。Hu等^［47，54］提出的压力梯度公式如下所示：

$${\left(\frac{\mathrm{1}}{\rho }\nabla p\right)}_i=\frac{\mathrm{1}}{m_i}\sum _j(p_i{V}_i^{\mathrm{2}}+p_j{V}_j^{\mathrm{2}})\cdot {\nabla }_i{W}_{ij}$$

（22）

$${\left(\frac{\mathrm{1}}{\rho }\nabla p\right)}_i=\frac{\mathrm{1}}{m_i}\sum _j(V_i^{\mathrm{2}}+V_j^{\mathrm{2}})\frac{{\rho }_i{p}_j+{\rho }_j{p}_i}{{\rho }_i+{\rho }_j}\cdot {\nabla }_i{W}_{ij}$$

（23）

针对水气界面处压力是否连续的问题，不同学者展开了讨论。Chen等^［13］基于压力在界面处连续这一假设将界面处粒子支持域内异相粒子换成同相粒子，对于式（23）在大密度比情况下界面处轻粒子的压力占据主要作用。对式（23）右侧的压力项进行了以下修改：

$$P_{\mathrm{g}\mathrm{l}}=\frac{{\rho }_{\mathrm{g}}{p}_{\mathrm{l}}+{\rho }_{\mathrm{l}}{p}_{\mathrm{g}}}{{\rho }_{\mathrm{g}}+{\rho }_{\mathrm{l}}}=\frac{{\rho }_{\mathrm{l}}{p}_{\mathrm{g}}+o({\rho }_{\mathrm{l}}{p}_{\mathrm{g}})}{{\rho }_{\mathrm{g}}+{\rho }_{\mathrm{l}}}\approx$$

$$\frac{{\rho }_{\mathrm{l}}}{{\rho }_{\mathrm{g}}+{\rho }_{\mathrm{l}}}{p}_{\mathrm{g}}\approx p_{\mathrm{g}}$$

（24）

式中：$p_{\mathrm{l}}$为流体压强；$p_{\mathrm{g}}$为气体压强；$o({\rho }_{\mathrm{l}}{p}_{\mathrm{g}})$为${\rho }_{\mathrm{l}}{p}_{\mathrm{g}}$的高阶无穷小量；$P_{\mathrm{g}\mathrm{l}}$为流体粒子与气体粒子间的压力。然而，Zhou等^［19］提出式（24）没有考虑到界面压力梯度不连续问题，会导致气相粒子的垂直加速度出现误差。因此，对界面处气相粒子支持域内的流体粒子应做如下修改：

$$p_{\mathrm{g}\mathrm{n}\mathrm{e}\mathrm{w}}=p_{\mathrm{g}}+\frac{c_p}{{\rho }_{\mathrm{g}\mathrm{0}}}({\rho }_{\mathrm{l}\mathrm{0}}-{\rho }_{\mathrm{g}\mathrm{0}})(p_{\mathrm{g}}-p_{{}_{\mathrm{l}}})$$

（25）

式中：$p_{\mathrm{g}\mathrm{n}\mathrm{e}\mathrm{w}}$为修正后的气体粒子压力；${\rho }_{\mathrm{g}\mathrm{0}}$为气体粒子参考密度；${\rho }_{\mathrm{l}\mathrm{0}}$为流体粒子参考密度；$c_p$为压力修正系数，其取值在0到1之间。Zhou等^［19］采用式（25）模拟了两相溃坝和楔体入水案例，结果表明式（25）在处理界面小变形问题中优势明显，也具有处理复杂界面问题的能力。也有学者对压力采用了更简单直接的取均值的处理方式^［55］，在两相溃坝和水中气泡上升案例中起到了很好的模拟效果。具体处理方式如下所示：

$${\left(\frac{\mathrm{1}}{\rho }\nabla p\right)}_i=\frac{\mathrm{1}}{m_i}\sum _j(V_i^{\mathrm{2}}+V_j^{\mathrm{2}})\frac{p_i+p_j}{\mathrm{2}}\cdot {\nabla }_i{W}_{ij}$$

（26）

2.3.2　表面张力项

由于同相粒子间的引力大于异相粒子间的引力，因此在气液交界处存在引力差，产生表面张力。在很多小尺度水气二相流问题（如气泡运动等）中，表面张力是不可忽略的。表面张力通常用2种方法计算，第1种是从分子动力学的角度，在同相粒子和异相粒子间分别施加系数不同的分子间作用力项^［56］，从而产生表面张力，形式如下所示：

$$\frac{\mathrm{1}}{\rho }{\mathit{F}}_{\mathbf{s}}=\left\{\begin{array}{l}\frac{\mathrm{2}{c}_{\mathrm{2}}}{m_i}\sum _j{\nabla }_i{W}_{ij},\mathrm{同相}\mathrm{粒子}\\ \frac{\mathrm{2}k{c}_{\mathrm{2}}}{m_i}\sum _j{\nabla }_i{W}_{ij},\mathrm{异相}\mathrm{粒子}\end{array}\right.$$

（27）

式中：c₂为范德华常数；k为小于1的系数，它决定了表面张力的大小。这种方法比较简单，但缺点是随着分辨率的增加，表面张力不会趋于一个定值。

第2种方法是从宏观粒子运动的角度，Brackbill等^［57］提出了连续表面张力（CSF）模型，通过计算界面的局部曲率来模拟表面张力，随着分辨率的增加表面张力会收敛到精确值，但CSF模型不满足动量守恒。此后CSF模型产生了一些变体，如Lafaurie等^［58］、Morris^［59］、Hu等^［47］采用动量法，避免了局部曲率的计算，但对于气液密度比高达1 000的水气二相流，需要对权函数进行插值处理来平滑界面。Adami等^［24］、Szewc等^［11］、Zainali等^［10］进一步改进了传统CSF模型，通过引入色标函数c来计算界面曲率，以Adami等^［24］的公式为例：

$${\mathit{F}}_{\mathbf{s}}=-\alpha \kappa \mathit{n}{\delta }_{\mathrm{s}}$$

（28）

$${\kappa }_i=-\nabla \cdot {\mathit{n}}_i$$

（29）

$${\mathit{n}}_i=\frac{\nabla c_i}{\left|\nabla c_i\right|}$$

（30）

$$c_i^j=\left\{\begin{array}{l}\mathrm{1},j\mathit{粒子}\mathit{和}i\mathit{粒子}\mathrm{异相}\\ \mathrm{0},j\mathit{粒子}\mathit{和}i\mathit{粒子}\mathit{同相}\end{array}\right.$$

（31）

$$\nabla c_i=\frac{\mathrm{1}}{V_i}\sum _j(V_i^{\mathrm{2}}+V_j^{\mathrm{2}})c_i^j{\nabla }_i{W}_{ij}$$

（32）

$$c_{ij}=\frac{{\rho }_j}{{\rho }_i+{\rho }_j}{c}_i^i+\frac{{\rho }_i}{{\rho }_i+{\rho }_j}{c}_j^i$$

（33）

式中：$\alpha$为表面张力系数；$\kappa$为曲率；n为界面处的单位法向量；${\delta }_{\mathrm{s}}$为表面狄拉克函数（表面取1，反之为0）；c为色标函数；c_ij为粒子间平均色标函数值。其他粒子法（如MPS等）也采用类似的做法求解表面张力^［60-61］。

2.3.3　界面斥力项

水气二相流界面易出现虚假破碎和不稳定问题，比较简单的方法是在不同粒子间的压力梯度中引入一个较小的斥力以保持界面清晰。Grenier等^［9］首次提出了这个概念并验证了可行性，界面斥力计算式为

$${\mathit{F}}_i^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}}=-\frac{{\epsilon }_{\mathrm{I}}}{V_i}\sum _j(\left|p_i\right|V_i^{\mathrm{2}}+\left|p_j\right|V_j^{\mathrm{2}}){\nabla }_i{W}_{ij}$$

（34）

式中：${\epsilon }_{\mathrm{I}}$为系数，其大小在0.01到0.1之间。注意，这里的下标i和j要求为异相粒子。

Monaghan等^［16］基于标准SPH算法提出了另外一种界面斥力的形式，这与Grenier等^［9］增加的斥力类似，但不需要其他复杂的算法，形式为

$$\frac{\mathrm{1}}{{\rho }_i}{\mathit{F}}_i^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}}=-\sum _j\mathrm{0.08}{m}_j\left|\frac{{\rho }_{\mathrm{0}i}-{\rho }_{\mathrm{0}j}}{{\rho }_{\mathrm{0}i}+{\rho }_{\mathrm{0}j}}\right|\left|\frac{p_i+p_j}{{\rho }_i{\rho }_j}\right|{\nabla }_i{W}_{ij}$$

（35）

式中：${\rho }_{\mathrm{0}i}$为i粒子的参考密度；${\rho }_{\mathrm{0}j}$为j粒子的参考密度。通过这种方式，当i和j粒子有相同的参考密度时界面力为零，即省去了界面判断环节。这种简单的处理方法可以有效地模拟两相密度比在1~1 000的情况^［16］。

后期有学者结合式（34）和式（35）的优点提出了新的界面斥力计算方式，如Ming等^［28］和徐丞君等^［55］。这里以Ming等^［28］提出的公式为例：

$${\mathit{F}}_i^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\mathrm{e}}=-\frac{\mathrm{0.08}}{V_i}\frac{\left|{\rho }_{\mathrm{l}}-{\rho }_{\mathrm{g}}\right|}{{\rho }_{\mathrm{l}}+{\rho }_{\mathrm{g}}}\sum _j(\left|p_i\right|V_i^{\mathrm{2}}+$$

$$\left|p_j\right|V_j^{\mathrm{2}}){\nabla }_i{W}_{ij}$$

（36）

式（36）能够较好地模拟气泡在水中的演化过程。可以发现，式（34）―（36）的本质是相同的，均是对2个不同相的粒子在动量方程中额外添加一个很小的压力梯度项从而产生两相间的排斥力，只是式（35）及式（36）不需要判断两相界面，一定程度上提高了计算效率。使用和不使用界面力模拟的液滴振荡过程如图2所示，可见界面力极大地提高了两相界面捕捉的精度。除了界面斥力以外，在模拟气泡的研究中有学者引入了界面阻力。Wan等^［34］指出，流体中气泡的速度是由阻力和浮力的平衡决定的，因此有必要考虑阻力的影响，这个力与气泡的运动方向相反，在流体动量方程中为正，在气体动量方程中为负，形式为

$$f_{\mathrm{g}\mathrm{l}}=\frac{\mathrm{3}}{\mathrm{4}}{C}_{\mathrm{D}}{\rho }_{\mathrm{l}}\frac{{\epsilon }_{\mathrm{g}}}{d_{\mathrm{p}}}\left|{\mathit{u}}_{\mathrm{l}}-{\mathit{u}}_{\mathrm{g}}\right|({\mathit{u}}_{\mathrm{l}}-{\mathit{u}}_{\mathrm{g}})$$

（37）

$$F_{i,\mathrm{g}\mathrm{l}}^{}={\displaystyle \sum _{j=\mathrm{1}}^{N_{\mathrm{r}\mathrm{w}}}}{f}_{j,\mathrm{g}\mathrm{l}}$$

（38）

式中：$f_{j,\mathrm{g}\mathrm{l}}$为作用于j粒子上的阻力；$C_{\mathrm{D}}$为阻力系数；$d_{\mathrm{p}}$为粒子间距；${\mathit{u}}_{\mathrm{l}}$为流体粒子速度；${\mathit{u}}_{\mathrm{g}}$为气体粒子速度；${\epsilon }_{\mathrm{g}}$为空间中某一点气体粒子的体积分数；$F_{i,\mathrm{g}\mathrm{l}}^{}$为i粒子的总阻力；$N_{\mathrm{r}\mathrm{w}}$为支持域内的流体粒子数。Wan等^［34］指出，考虑该界面阻力可以使气泡上浮的模拟结果与物理试验结果更加吻合。

图2  采用和不采用界面斥力项对振荡液滴的SPH模拟结果^［28］

Fig.2  Simulation results of oscillating liquid drop in SPH model with and without interface sharpness force^[28]

2.4　其他数值技术

除了针对控制方程的修正处理，其他一些提高计算精度和稳定性的数值技术在水气二相流SPH方法中也经常使用，由于水气二相流SPH方法对数值稳定性有较高要求，因此这些数值技术在水气二相流SPH方法中得到广泛应用。

2.4.1　密度过滤

在弱可压SPH方法中，压力直接由密度决定，因此密度的精确计算非常重要。水气二相流的界面会出现密度、质量间断的问题，相比单相流更易出现密度振荡和压力不稳定。因此，周期性地采用密度求和法过滤密度成为一种常用的数值稳定方法^［7，13，18］。以Colagrossi等^［7］提出的一种一阶精确插值方法为例，如下所示：

（39）

式中：为过滤后i粒子的密度；$W_j^{\mathrm{W}\mathrm{L}\mathrm{S}}$为j粒子的移动最小二乘核函数。式（39）于每20个时间步长更新一次。

2.4.2　XSPH修正

为了防止粒子发生非物理穿透，Monaghan^［20］提出了一种速度修正方法（称为XSPH方法），用来修正速度散度项中的粒子速度，从而改进粒子分布。XSPH方法采用支持域内其他粒子的速度加权平均对当前粒子进行修正，修正后的速度项用于式（3）的计算中，因此XSPH方法本质上是对粒子位置修正而不是对速度修正，形式为

（40）

$$\mathrm{\Delta }{\mathit{u}}_i=\frac{\epsilon }{\mathrm{2}}\sum _j\frac{m_j}{\overline{{\rho }_{ij}}}({\mathit{u}}_j-{\mathit{u}}_i)W_{ji}$$

（41）

$$\overline{{\rho }_{ij}}=\frac{{\rho }_i+{\rho }_j}{\mathrm{2}}$$

（42）

式中：为修正后的当前粒子速度；$\epsilon$为系数；$m$为粒子质量；$\rho$为粒子密度。在单相流SPH方法中，XSPH方法在堰流和孔洞出流^［62］、水跃湍流模型^［63］中均有应用。Colagrossi等^［7］将XSPH方法引入水气二相流SPH方法中，XSPH方法有助于获得更加光滑清晰的异相界面。此后，XSPH方法修正为水气二相流SPH的常用方法^［35，64］。值得注意的是，在流体雷诺数偏高时使用XSPH方法会导致粒子混乱^［65］，这个影响可以通过增加粒子分辨率来减轻。

2.4.3　位移修正

在SPH方法数值计算中，由于对粒子的局部扩散性处理不好，因此易出现粒子局部聚集问题，进而影响计算精度。Xu等^［66］提出了粒子位移修正方法，用来解决粒子聚集问题。Lind等^［67］基于Fick扩散定律对粒子位移修正算法进行了改进。Sun等^［68］将该算法引入δ-SPH模型中，用于模拟自由液面流动问题。Antuono等^［69］将粒子位移修正技术与任意拉格朗日-欧拉格式相结合，在确保粒子均匀分布的同时，进一步提高了粒子位移修正算法的精度。在水气二相流SPH模拟中，由于气体密度较小，更易产生非物理空洞和两相分离问题，导致数值计算的不稳定，因此在二相流中引入位移修正方法也是很有必要的。Mokos等^［14］将位置修正扩展到多相流问题中，在流体和气体中分别采用不同的位移修正公式，使得气相中的空洞问题得到显著改善。王平平等^［18］指出上述方法过于复杂，并基于δ-SPH方法改进了二相流的位移修正算法，两相流体的处理方法一致，通过取不同的经验系数来体现两相中位移修正的差异，获得了良好的模拟效果。

2.4.4　黎曼求解器

黎曼求解器是通过对激波问题的更好求解，使得SPH方法求解压力场时非物理波动显著减弱。黎曼求解器较早由Inutsuka^［70］引入SPH算法，可以不依赖人工黏度和δ-SPH方法即能取得较好的压力分布。Rezavand等^［15］将黎曼求解器扩展到水气二相流SPH问题，该方法消除了剧烈流动中的非物理空洞，获得了清晰准确的异相界面，并且两相采用相同的声速，提高了计算效率。Yang等^［26］在Rezavand等^［15］方法的基础上，对黎曼耗散项进行了改进，解决了Rezavand等^［15］方法不能应对黏度变化的问题，使得基于黎曼求解器的水气二相流方法能够稳定模拟大密度、不同黏度的二相流过程。

2.4.5　粒子体积自适应技术

传统SPH方法大多应用于不可压或弱可压水气二相流问题，针对强可压水气二相流问题，如水下爆炸、空化等，常常涉及到气体相（或称气泡）体积的大幅变化。气体强压缩性会导致气体粒子体积膨胀或收缩，而水作为一种几乎不可压介质，其粒子体积变化相对气体而言变化很小，因此导致水气界面两侧的粒子分布严重不均匀。为解决这个问题，Sun等^［30-31，71］提出一种粒子体积自适应技术（VAS）。该算法中，对过度膨胀的气体粒子进行撕裂，对过度收缩的气体粒子进行融合，有效确保计算全程每个气体粒子的体积变化控制在一定范围以内，从而提高强可压多相流SPH方法数值计算的精度。

3 结论与展望

为了保证计算稳定性，水气二相流SPH方法在状态方程、连续方程、动量方程中均采用了一些针对性的处理方式。状态方程中，水气二相流SPH方法普遍添加一个背景压力，背景压力常取定值或最大静水压力的倍数。连续方程中，水气二相流SPH方法没有区别于单相流SPH方法的独有处理方式，但常引入一些增加数值稳定性的技术。动量方程中，水气二相流SPH方法基于两相的粒子体积和密度差异对压力梯度进行修正，并通常会在方程中添加表面张力项和界面斥力项。还有其他一些数值技术常用于水气二相流SPH方法，如密度过滤、XSPH修正、位移修正、黎曼求解器、粒子体积自适应技术等，以获得更好的密度场和更清晰的两相界面。不同物理过程（如两相溃坝、气泡上升、水下爆炸等）的水气二相流SPH模拟需要采用的处理方式均有所区别，要根据具体问题来选用合适的控制方程与数值技术。

随着GPU加速技术的普遍应用，SPH方法在实际工程问题中有了更为广泛的应用前景，水气二相流SPH方法的模拟技术也日趋成熟，因此期待未来研究中水气二相流SPH方法在解决实际问题方面发挥更大的作用，如极端条件下的海工结构受力、水下爆破、空蚀等。具体技术方面，涉及水气二相流的开边界问题、湍流和空化模型的建立以及两相之间速度差较大时界面处的摩擦力等问题需要进一步研究解决。

作者贡献声明

潘 毅：框架设计，写作与修改。

刘玉钰：写作与修改。

汪明智：查阅资料，论文写作。

匡翠萍：论文写作。

倪兴也：论文审阅。

孙鹏楠：学术指导，图表绘制。