界面追踪法求解流体流动的表面张力

本文基于OpenFOAM软件，采用移动网格界面追踪方法，模拟三维方向上不可压缩、不互溶的两相流界面流动问题。界面处网格被划分为两部分，通过施加一定的运动学和动力学条件，利用迭代方法（PISO）求解网格相应参数。

文章目录

1. 数学模型
2. 数值方法
- 2.1计算域离散
- 2.2 数学模型离散
- 2.3 表面追踪法
- - 2.3.1 Af 处的 p、vp、\mathbf vp、v 值
  - 2.3.2 Bf 处的 p、vp、\mathbf vp、v 值
  - 2.3.3 净质量通量修正
- 2.4 曲面导数的计算
- 2.5 表面张力
3. 解题步骤
4. 质量守恒和迪利克雷流体边界条件

1. 数学模型

前提假设：牛顿不可压缩流体、绝热、不互溶；网格体积 V，运动表面 S。

根据连续性方程以及动量定理，可知：

∮Sρn⋅vdS=0(1)\oint_S\rho\mathbf n\cdot \mathbf vdS=0\tag1∮Sρn⋅vdS=0(1)

ddt∫VρvdV+∮Sn⋅ρ(v−vS)dS=∮Sn⋅(μ∇v)dS−∫V∇pdV(2)\frac d{dt}\int_V\rho \mathbf vdV+\oint_S\mathbf n\cdot\rho(\mathbf v-\mathbf v_S)dS=\oint_S\mathbf n\cdot(\mu\nabla\mathbf v)dS-\int_V\nabla pdV\tag2 dtd∫VρvdV+∮Sn⋅ρ(v−vS)dS=∮Sn⋅(μ∇v)dS−∫V∇pdV(2)

其中，

n\bf nn：表面单位外法线；
v\bf vv：流体速度；
vS\mathbf v_SvS：表面流体速度；
ρ\rhoρ：流体密度；
μ\muμ：流体动力粘度；
ppp：流体动压（dynamic pressure），动压=总压-静水压（hydrostatic pressure）。

在该方程中，p=pρ(3)p=\frac p\rho\tag3p=ρp(3)

体积 V 的变化率与表面速度 vS\mathbf v_SvS 的关系为：

ddt∫VdV−∮Sn⋅vSdS=0(4)\frac d{dt}\int_VdV-\oint_S\mathbf n\cdot\mathbf v_SdS=0\tag4dtd∫VdV−∮Sn⋅vSdS=0(4)

为了保证界面的连续性，界面两侧速度必须相等（kinematic condition）：

vA=vB(5)\mathbf v_A=\mathbf v_B\tag5vA=vB(5)

同样地，必须满足动量守恒定律（dynamic condition），即流体边界上的力是连续的。其中，切向力与切向速度的方向导数有关，形式如下：

μB(n⋅∇vt)B−μA(n⋅∇vt)A=−∇Sσ−(μB−μA)(∇svn)(6)\mu_B(\mathbf n\cdot\nabla\mathbf v_t)_B-\mu_A(\mathbf n\cdot\nabla\mathbf v_t)_A=-\nabla_S\sigma-(\mu_B-\mu_A)(\nabla_sv_n)\tag6μB(n⋅∇vt)B−μA(n⋅∇vt)A=−∇Sσ−(μB−μA)(∇svn)(6)

n\bf nn：流体交界面单位法向量，由流体A指向流体B；
vt=(I−nn)⋅v\mathbf v_t=(\mathbf I-\mathbf n\mathbf n)\cdot \mathbf vvt=(I−nn)⋅v：切向速度；
∇S=∇−nn⋅∇\nabla_S=\nabla-\mathbf n\mathbf n\cdot\nabla∇S=∇−nn⋅∇：表面梯度算子；
σ\sigmaσ：表面张力系数；
vn=n⋅vv_n=\mathbf n\cdot\mathbf vvn=n⋅v：表面法向速度。

对于表面张力系数的切向梯度 (∇Sσ\nabla_S\sigma∇Sσ)，当表面活性剂分布不均匀或存在温度梯度时，该值不为零。

根据界面法向力平衡，压力跃变为：

pB−pA=σκ−2(μB−μA)∇S⋅v(7)p_B-p_A=\sigma\kappa-2(\mu_B-\mu_A)\nabla_S\cdot\mathbf v\tag7pB−pA=σκ−2(μB−μA)∇S⋅v(7)

其中，

κ=−∇S⋅n\kappa=-\nabla_S\cdot\bf nκ=−∇S⋅n：界面平均曲率的两倍；
RHS第二项：法向粘性力在界面处的跃变，表现为界面速度的表面离散。

2. 数值方法

这一部分可参考基于FVM的应力求解。

2.1计算域离散

在界面追踪过程中，需要根据界面随时间的变化，对有限体网格进行调整。本文采用变形网格方法，在网格拓扑结构不变的情况下，根据边界点的运动规则调整内部控制体顶点。这里，网格顶点的位移 u\bf uu 由拉普拉斯方程确定：
∇⋅(Γ∇u)=0(8)\nabla\cdot(\Gamma\nabla\mathbf u)=0\tag8∇⋅(Γ∇u)=0(8)

Γ\GammaΓ ：与到移动边界的距离的平方成反比；则越靠近移动界面，网格刚性越强，保证了边界处的网格质量；

2.2 数学模型离散

根据 FVM，积分守恒方程可以转化为面积分的加和形式；面积积分和体积积分采用中点法则逼近二阶精度；时间离散采用隐式的三层二阶格式(three-level second order scheme)，即向后差分。

对于移动的控制体 VPV_PVP，动量方程离散形式如下：
ρP3vPnVPn−4vPoVPo+vPooVPoo2Δt+∑f(m˙fn−ρfV˙fn)vfn=∑fμfnf⋅(∇v)fnSfn−(∇p)pnVPn(9)\begin{aligned}&\rho_P\frac{3\mathbf v_P^nV_P^n-4\mathbf v_P^oV_P^o+\mathbf v_P^{oo}V_P^{oo}}{2\Delta t}+\sum_f(\dot{\mathbf m}_f^n-\rho_f\dot{V}_f^n)\mathbf v_f^n\\&=\sum_f\mu_f\mathbf n_f\cdot(\nabla\mathbf v)_f^nS_f^n-(\nabla p)_p^nV_P^n\end{aligned}\tag9ρP2Δt3vPnVPn−4vPoVPo+vPooVPoo+f∑(m˙fn−ρfV˙fn)vfn=f∑μfnf⋅(∇v)fnSfn−(∇p)pnVPn(9)

下角标 P、f_P、_fP、f 分别代表网格体心和面心的值；
上角标 n,o,oon, o, oon,o,oo 分别代表新值（new time)，旧值（old time）、旧-旧值（old-old time），un=u(t+Δt)u^n=u(t+\Delta t)un=u(t+Δt)、uo=u(t)u^o=u(t)uo=u(t)、uoo=u(t−Δt)u^{oo}=u(t-\Delta t)uoo=u(t−Δt)；
面心质量流 m˙fn=ρfnfn⋅vfnSfn\color{red}{\dot{\mathbf m}_f^n=\rho_f\mathbf n_f^n\cdot\mathbf v_f^nS_f^n}m˙fn=ρfnfn⋅vfnSfn ，其必须满足质量守恒方程；
面心体积通量 VfnV_f^nVfn 必须满足离散化的 GCL 格式；
其他的通量做显式处理，用以线性化对流项。

面心处的值大多采用与相邻体心值线性插值求得：
ϕfn=(ϕPn)f‾=fxnϕPn+(1−fxn)ϕNn(10)\phi_f^n=\overline{(\phi_P^n)_f}=f_x^n\phi_P^n+(1-f_x^n)\phi_N^n\tag{10}ϕfn=(ϕPn)f=fxnϕPn+(1−fxn)ϕNn(10)

其中，

ϕ\phiϕ 为一般的因变量；
fx=fN‾/PN‾f_x=\overline{fN}/\overline{PN}fx=fN/PN 为插值系数；

这里假设线 PN‾\overline{PN}PN 与面 fff 相交于面心位置，否则，必须引入一个偏度（skewness）修正的线性插值。详见多材料线弹性固体的应力分析 2.2.1节。

为了保证方程二阶精度的同时具有有界性，式（9）对流项中的面心速度 vfn\mathbf v_f^nvfn 需要采用混合格式（又称 Gamma 差分格式）。

式（9）中的拉普拉斯项 nfn⋅(∇v)fn\mathbf n_f^n\cdot(\nabla \mathbf v)_f^nnfn⋅(∇v)fn 可离散为：

nfn⋅(∇v)fn=∣Δfn∣vNn−vPn∣dn∣⏟orthogonal+(nfn−Δfn)⋅(∇v)fn⏟non−orthogonal(11)\mathbf n_f^n \cdot(\nabla \mathbf v)_f^n= \underbrace{|\mathbf \Delta_f^n|\frac{\mathbf v_N^n-\mathbf v_P^n}{|\mathbf d^n|}}_{orthogonal} +\underbrace{(\mathbf n_f^n-\mathbf \Delta_f^n)\cdot(\nabla \mathbf v)_f^n}_{non-orthogonal} \tag{11} nfn⋅(∇v)fn=orthogonal∣Δfn∣∣dn∣vNn−vPn+non−orthogonal(nfn−Δfn)⋅(∇v)fn(11)
其中，Δf=dn/(dn⋅nf)\mathbf \Delta_f=\mathbf d^n/(\mathbf d^n\cdot \mathbf n_f)Δf=dn/(dn⋅nf)。

注意：正交部分采用隐式处理，非正交部分采用显式处理。

式（9）中的梯度项 ∇p\nabla p∇p 采用高斯积分理论进行计算。其他的体心梯度项则采用最小二乘法。这一方法不用考虑局部网格质量，而得到二阶精度梯度。

根据以上规则，动量方程离散形式可以写成线性代数方程的形式：

aPvPn+∑NaNvNn=rPn−(∇p)Pn(12)a_P\mathbf v_P^n+\sum_Na_N\mathbf v_N^n=\mathbf r_P^n-(\nabla p)_P^n\tag{12}aPvPn+N∑aNvNn=rPn−(∇p)Pn(12)

aPa_PaP：对角线系数
aNa_NaN：临近系数（非对角线元素）
rPn\mathbf r_P^nrPn：源项，由离散过程中一些显性处理的速度场构成（如：体心质量流、非正交修正项等）

原始的 Rhie-Chow 插值依赖于时间步长，对角线上的非稳态项及源项变得十分重要：

aP=aP∗+3ρP2Δt(13)a_P=a_P^*+\frac{3\rho_P}{2\Delta t}\tag{13}aP=aP∗+2Δt3ρP(13)

rP=rP∗+2ρPVPoVPnΔtvPo−ρPVPoo2VPnΔtvPoo(14)\mathbf r_P=\mathbf r_P^*+\frac{2\rho_PV_P^o}{V_P^n\Delta t}\mathbf v_P^o-\frac{\rho_PV_P^{oo}}{2V_P^n\Delta t}\mathbf v_P^{oo}\tag{14}rP=rP∗+VPnΔt2ρPVPovPo−2VPnΔtρPVPoovPoo(14)

上角标 ∗*∗ 代表扩散项和对流项等。

流体流动的数学模型采用分离程序求解，对速度方程与压力方程进行解耦。压力方程的构建基于连续性方程而来，其通过对动量方程做散度并相加获得。对于当前网格 P，连续性方程离散形式如下：
}
∑fm˙fn=∑fρfnfn⋅vfnSfn=0(15)\sum_f\dot{\mathbf m}_f^n=\sum_f\rho_f\mathbf n_f^n\cdot\mathbf v_f^nS_f^n=0\tag{15}f∑m˙fn=f∑ρfnfn⋅vfnSfn=0(15)

根据式（12)，单元体心速度如下：

vPn=HP(vn)aP−1aP(∇p)Pn(16)\mathbf v_P^n=\frac{\mathbf H_P(\mathbf v^n)}{a_P}-\frac 1{a_P}(\nabla p)_P^n\tag{16}vPn=aPHP(vn)−aP1(∇p)Pn(16)

其中，

HP(vn)=−∑faNvNn+rPn=HP∗(vn)+2ρPVPoVPnΔtvPo−ρPVPoo2VPnΔtvPoo(17)\mathbf H_P(\mathbf v_n)=-\sum_fa_N\mathbf v_N^n+\mathbf r_P^n=\mathbf H_P^*(\mathbf v^n)+\frac{2\rho_PV_P^o}{V_P^n\Delta t}\mathbf v_P^o-\frac{\rho_PV_P^{oo}}{2V_P^n\Delta t}\mathbf v_P^{oo}\tag{17}HP(vn)=−f∑aNvNn+rPn=HP∗(vn)+VPnΔt2ρPVPovPo−2VPnΔtρPVPoovPoo(17)

根据式（16）对单元面心速度进行改写：

vfn=(Ha)f−(1a)f(∇p)fn(18)\mathbf v_f^n=(\frac{\mathbf H}{a})_f-(\frac 1{a})_f(\nabla p)_f^n\tag{18}vfn=(aH)f−(a1)f(∇p)fn(18)

(Ha)f(\frac{\mathbf H}{a})_f(aH)f、(1a)f(\frac 1{a})_f(a1)f 是由式（16）中对应项对共享面 fff 的两个网格线性插值得到 (一般采用 Rhie-Chow 插值格式)。

将式（18）带入到式（15）中，可得：

∑f(1a)fnfn⋅(∇p)fnSf=∑fnfn⋅(Ha)fSfn(19)\sum_f(\frac 1a)_f\mathbf n_f^n\cdot(\nabla p)_f^nS_f=\sum_f\mathbf n_f^n\cdot(\frac {\mathbf H}a)_fS_f^n\tag{19}f∑(a1)fnfn⋅(∇p)fnSf=f∑nfn⋅(aH)fSfn(19)

法向压力拉普拉斯项 nfn⋅(∇p)fn\mathbf n_f^n\cdot(\nabla p)_f^nnfn⋅(∇p)fn 采用式（11）方式进行离散。对式（19）进行求解后，可以得到自由发散的面心质量通量：

m˙fn=ρf[nfn⋅(Ha)f−(1a)fnfn⋅(∇p)fn]Sf(20)\dot{\mathbf m}_f^n=\rho_f[\mathbf n_f^n\cdot(\frac{\mathbf H}{a})_f-(\frac 1{a})_f\mathbf n_f^n\cdot(\nabla p)_f^n]S_f\tag{20}m˙fn=ρf[nfn⋅(aH)f−(a1)fnfn⋅(∇p)fn]Sf(20)

根据原始的 Rhie-Chow 插值格式，有：

(Ha)fRC=[Hp(vn)aP]‾f，(1a)fRC=(1aP)‾f(21)(\frac{\mathbf H}{a})_f^{RC}=\overline{ \left[\frac{\mathbf H_p(\mathbf v^n)}{a_P}\right]}_f，(\frac 1a)_f^{RC}=\overline{(\frac 1{a_P})}_f \tag{21}(aH)fRC=[aPHp(vn)]f，(a1)fRC=(aP1)f(21)

()‾f\overline{()}_f()f 代表与相邻网格体心值进行线性插值。

注意：当时间步长非常小时，这种 Rhie-Chow 插值格式并不能保证压力场的准确性。

lim⁡Δt→0(Ha)fRC=43vPo‾f−13vPoo‾f，lim⁡Δt→0(1a)fRC=0(22)\lim_{\Delta t\to 0}(\frac{\mathbf H}{a})_f^{RC}=\frac 43\overline{\mathbf v_P^o}_f-\frac 13\overline{\mathbf v_P^{oo}}_f，\lim_{\Delta t\to 0}(\frac 1a)_f^{RC}=0 \tag{22}Δt→0lim(aH)fRC=34vPof−31vPoof，Δt→0lim(a1)fRC=0(22)

由于线性插值得到的速度并不满足连续性方程，在考虑极限的情况下会产生一个非物理的压力场。另外，对于原始的Rhie-Chow 插值格式，其体心质量流依赖于时间步长。

对此，提出一种改进的Rhie-Chow 插值格式（mRC）：

(Ha)fmRC=[Hp∗(vn)]‾f(aP)‾f+(2ρPVPoVPnΔt)‾f(aP)‾f(vPo)f−(ρPVPoo2VPnΔt)f‾(aP)‾f(vPoo)f，(1a)fmRC=1(aP)f‾(23)\begin{aligned}&(\frac{\mathbf H}{a})_f^{mRC}= \frac{\overline{[\mathbf H_p^*(\mathbf v^n)]}_f}{\overline{(a_P)}_f}+\frac{\overline{(\frac{2\rho_PV_P^o}{V_P^n\Delta t})}_f}{\overline{(a_P)}_f}(\mathbf v_P^o)_f-\frac{\overline{(\frac{\rho_PV_P^{oo}}{2V_P^n\Delta t})_f}}{\overline{(a_P)}_f}(\mathbf v_P^{oo})_f， \\[1.5ex]&(\frac 1a)_f^{mRC}=\frac 1{\overline{(a_P)_f }}\end{aligned}\tag{23}(aH)fmRC=(aP)f[Hp∗(vn)]f+(aP)f(VPnΔt2ρPVPo)f(vPo)f−(aP)f(2VPnΔtρPVPoo)f(vPoo)f，(a1)fmRC=(aP)f1(23)

在上述改进的 Rhie-Chow 插值格式中，(vpo)f(\mathbf v_p^o)_f(vpo)f、(vpoo)f(\mathbf v_p^{oo})_f(vpoo)f 均满足其所在时间步下的质量守恒方程：
(vpo)f=(vpo)f‾+nfo[m˙foρfSfo−nfo⋅(vpo)f‾],(vpoo)f=(vpoo)f‾+nfoo[m˙fooρfSfoo−nfoo⋅(vpoo)f‾](24)\begin{aligned}&(\mathbf v_p^o)_f=\overline{(\mathbf v_p^o)_f}+\mathbf n_f^o\left[\frac{\dot m_f^o}{\rho_fS_f^o}-\mathbf n_f^o\cdot\overline{(\mathbf v_p^o)_f} \right] ,\\[1.5ex] &(\mathbf v_p^{oo})_f=\overline{(\mathbf v_p^{oo})_f}+\mathbf n_f^{oo}\left[\frac{\dot m_f^{oo}}{\rho_fS_f^{oo}}-\mathbf n_f^{oo}\cdot\overline{(\mathbf v_p^{oo})_f} \right]\end{aligned}\tag{24}(vpo)f=(vpo)f+nfo[ρfSfom˙fo−nfo⋅(vpo)f],(vpoo)f=(vpoo)f+nfoo[ρfSfoom˙foo−nfoo⋅(vpoo)f](24)

mRC 插值格式具有以下优点：

即使时间步趋于零；
对于稳态问题，式（20）与时间步长无关。

瞬态情况下，为了保证动网格中的面心体积流 V˙fn\dot V_f^nV˙fn 满足DGCL（discretised geometric conservation law）, 有：

3VPn−4VPo+VPoo2Δt−∑fV˙fn=0(25)\frac{3V_P^n-4V_P^o+V_P^{oo}}{2\Delta t}-\sum_f\dot V_f^n=0\tag{25}2Δt3VPn−4VPo+VPoo−f∑V˙fn=0(25)

网格体积对时间的差分可写成以下形式：

VPn−VPo=∑fδVfn(26)V_P^n-V_P^o=\sum_f\delta V_f^n\tag{26}VPn−VPo=f∑δVfn(26)

其中，δVPn\delta V_P^nδVPn 为时间间隔内网格表面 fff 移动过程中扫过的体积。将式（26）带入到式（25）中，可得：

12Δt∑f(3δVfn−δVfo)=∑fV˙fn(27)\frac 1{2\Delta t}\sum_f(3\delta V_f^n-\delta V_f^o)=\sum_f\dot{V}_f^n\tag{27}2Δt1f∑(3δVfn−δVfo)=f∑V˙fn(27)

由此，网格面心体积流表达式如下：

V˙fn=32δVfnΔt−12δVfoΔt(28)\dot V_f^n=\frac 32\frac{\delta V_f^n}{\Delta t}-\frac 12\frac{\delta V_f^o}{\Delta t}\tag{28}V˙fn=23ΔtδVfn−21ΔtδVfo(28)

这种向后差分格式可以保证动网格具有二阶精度。

2.3 表面追踪法

对于带有尖锐表面的两相流，通常采用动网格表面追踪法（interface tracking procedure）进行数值模拟。将代求网格分为两部分，每一部分仅包含一种理想流体。

如上图所示，在相边界处，两部分网格的接触面分别为：SA、SBS_A、S_BSA、SB；接触面由一系列的交界面组成，接触面 SAS_ASA 上每一个表面 AfA_fAf 对应于 SBS_BSB 上的表面 BfB_fBf，两者在几何形状上完全相同。交界面处相匹配的网格采用界面追踪法；对于不参与匹配的网格，采用具有二阶精度的 inverse-distance 加权插值。

在相边界处施加合适的边界条件，以完成流动方程的耦合。若 ρA>ρB\rho_A>\rho_BρA>ρB，则：

相边界 Af：动压 pAp_ApA 与速度的法向导数 nA⋅(∇v)A\mathbf n_A\cdot\mathbf (\nabla \mathbf v)_AnA⋅(∇v)A 已知；
相边界 Bf：速度 vB\mathbf v_BvB 和动压的法向导数 nB⋅(∇p)B\mathbf n_B\cdot(\nabla p)_BnB⋅(∇p)B 已知。

外迭代过程中，采用以下形式更新边界条件：

注意：以下所有的值都是根据当前时间步求得的，因此上角标 n 省略

2.3.1 Af 处的 p、vp、\mathbf vp、v 值

界面处的切速度可根据式（6）与式（5）求得：

(vAf)t=μA(δBf)n(vAp)t+μB(δAf)n(vBp)tμA(δBf)n+μB(δAf)n+(δAf)n(δBf)n[(∇Sδ)Af+(μB−μA)(∇svn)Af]μA(δBf)n+μB(δAf)n(29)\begin{aligned}(\mathbf v_{Af})_t&=\frac{\mu_A (\delta_{Bf})_n(\mathbf v_{Ap})_t+\mu_B (\delta_{Af})_n(\mathbf v_{Bp})_t}{\mu_A (\delta_{Bf})_n+\mu_B (\delta_{Af})_n}\\[2ex]&+ \frac{(\delta_{Af})_n(\delta_{Bf})_n\left[(\nabla_S\delta)_{Af}+(\mu_B-\mu_A)(\nabla_sv_n)_{Af}\right]}{\mu_A (\delta_{Bf})_n+\mu_B (\delta_{Af})_n}\end{aligned}\tag{29}(vAf)t=μA(δBf)n+μB(δAf)nμA(δBf)n(vAp)t+μB(δAf)n(vBp)t+μA(δBf)n+μB(δAf)n(δAf)n(δBf)n[(∇Sδ)Af+(μB−μA)(∇svn)Af](29)

其中，

(vAp)t(\mathbf v_{Ap})_t(vAp)t、(vBp)t(\mathbf v_{Bp})_t(vBp)t：网格 AP 与 BP 体心处的切向速度；
(δAf)n(\delta_{Af})_n(δAf)n、(δBf)n(\delta_{Bf})_n(δBf)n：网格 AP 与 BP 体心到相应边界 Af 、Bf 的法向距离；
RHS 第一项采用了几何调和平均法（geometric harmonic mean）；
通过式（29)：式（5）中切向速度的法向离散采用了单边的一阶精度近似离散。

法向速度导数可根据更新后的边界速度求得：

nAfn⋅(∇v)Afn=μA(vAf)t−(vBf)t(δAf)n+μA(∇S⋅v)Afnf(30)\mathbf n_{Af}^n \cdot(\nabla \mathbf v)_{Af}^n=\mu_A\frac{(\mathbf v_{Af})_t-(\mathbf v_{Bf})_t}{(\delta_{Af})_n}+\mu_A(\nabla_S\cdot\mathbf v)_{Af}\mathbf n_f\tag{30}nAfn⋅(∇v)Afn=μA(δAf)n(vAf)t−(vBf)t+μA(∇S⋅v)Afnf(30)

nAf=nBf\mathbf n_{Af}=\mathbf n_{Bf}nAf=nBf ：交界面 Af 的单位法向量。

动压表达式为：

pAf=pBf−(ρA−ρB)g⋅rAf−(σκ)Af−2(μA−μB)(∇S⋅v)Af(31)p_{Af}=p_{Bf}-(\rho_A-\rho_B)\mathbf g\cdot \mathbf r_{Af}-(\sigma\kappa)_{Af}-2(\mu_A-\mu_B)(\nabla_S\cdot\mathbf v)_{Af}\tag{31}pAf=pBf−(ρA−ρB)g⋅rAf−(σκ)Af−2(μA−μB)(∇S⋅v)Af(31)

rAf\mathbf r_{Af}rAf：相边界 Af 的面心位置矢量;
(σκ)Af(\sigma\kappa)_{Af}(σκ)Af：表面张力的法向分量，参见式（55）。

2.3.2 Bf 处的 p、vp、\mathbf vp、v 值

根据式（4）可求得交界面 Bf 上的切速度：

(vBf)t=(vAf)t(32)(\mathbf v_{Bf})_t=(\mathbf v_{Af})_t\tag{32}(vBf)t=(vAf)t(32)

由于交界面 Bf 的净质量通量为零，即：
m˙Bf−ρBV˙Bf=0(33)\dot m_{Bf}-\rho_B\dot V_{Bf}=0\tag{33}m˙Bf−ρBV˙Bf=0(33)

可求得其法向速度分量：

(vBf)n=V˙BfSBfnBf(34)(\mathbf v_{Bf})_n=\frac{\dot{V}_{Bf}}{S_{Bf}}\mathbf n_{Bf}\tag{34}(vBf)n=SBfV˙BfnBf(34)

V˙Bf=V˙Af\dot{V}_{Bf}=\dot{V}_{Af}V˙Bf=V˙Af：交界面 Bf 的体积流量

界面处动压的法向导数计算结果如下：

nBf⋅(∇p)Bf=−nBf⋅(DvDt)Bf(35)\mathbf n_{Bf}\cdot(\nabla p)_{Bf}=-\mathbf n_{Bf}\cdot\left(\frac{D \mathbf v}{D t}\right)_{Bf}\tag{35}nBf⋅(∇p)Bf=−nBf⋅(DtDv)Bf(35)

2.3.3 净质量通量修正

通常情况下，外迭代结束时，通过相边界 Af 的净质量通量不为零，即：

m˙Af−ρAV˙Af≠0(36)\dot m_{Af}-\rho_A\dot V_{Af}\not=0\tag{36}m˙Af−ρAV˙Af=0(36)

为了修正该界面处的净质量通量，在单元面心上方引入一个控制点（control point）AcA_cAc，其位置见下图：

为了实现净质量通量的修正，必须移动界面点（interface points），以修正体积流量：

V˙Af′=m˙AfρA−V˙Af(37)\dot V_{Af}^{'}=\frac{\dot{m}_{Af}}{\rho_A}-\dot V_{Af}\tag{37}V˙Af′=ρAm˙Af−V˙Af(37)

界面点的位置需要根据控制点（control point）AcA_cAc 进行修正，具体方法如下：

a. 相界面 Af 从当前位置到修正位置所扫过的体积为：

δVAf′=23V˙Af′Δt(38)\delta V_{Af}^{'}=\frac 23 \dot V_{Af}^{'}\Delta t\tag{38} δVAf′=32V˙Af′Δt(38)

注：式（38）由式（28）衍生而来；

b. 控制点 AcA_cAc 沿 fAC\mathbf f_{A_C}fAC 方向的位移：

hAc′=δVAf′SAfnAf⋅fAC(39)h_{A_c}^{'}=\frac{\delta V_{Af}^{'}}{S_{Af}\mathbf n_{Af}\cdot\mathbf f_{A_C}}\tag{39}hAc′=SAfnAf⋅fACδVAf′(39)

则控制点的修正位置如下：

rAc=rAcp+hAc′fAc(40)\mathbf r_{A_c}=\mathbf r_{A_c}^p+h_{A_c}^{'}\mathbf f_{A_c}\tag{40}rAc=rAcp+hAc′fAc(40)

rAcp\mathbf r_{A_c}^prAcp：校正前 AcA_cAc 的位移矢量

c. 界面网格顶点 Ai 的位移

利用最小二乘法，将其放在相应控制点构成的平面上，其位置矢量表达式如下：

rAi=rAip+NAi⋅(p−rAip)NAi⋅fAifAi(41)\mathbf r_{Ai}=\mathbf r_{Ai}^p+\frac{\mathbf N_{Ai}\cdot(\mathbf p-\mathbf r_{Ai}^p)}{\mathbf N_{Ai}\cdot\mathbf f_{Ai}}\mathbf f_{Ai}\tag{41}rAi=rAip+NAi⋅fAiNAi⋅(p−rAip)fAi(41)

fAi\mathbf f_{Ai}fAi：边界网格顶点的位移方向；
NAi\mathbf N_{Ai}NAi：面的单位法向量

NAi=G−1⋅∑AcωAc2rAc∣G−1⋅∑AcωAc2rAc∣(42)\mathbf N_{Ai}=\frac{\mathbf G^{-1}\cdot\sum_{A_c}\omega_{A_c}^2\mathbf r_{A_c}}{|\mathbf G^{-1}\cdot\sum_{A_c}\omega_{A_c}^2\mathbf r_{A_c}|}\tag{42}NAi=∣G−1⋅∑AcωAc2rAc∣G−1⋅∑AcωAc2rAc(42)

p\mathbf pp：点在平面上的位置

p=∑AcωAc2rAc∑AcωAc2(43)\mathbf p=\frac{\sum_{A_c}\omega_{A_c}^2\mathbf r_{A_c}}{\sum_{A_c}\omega_{A_c}^2}\tag{43}p=∑AcωAc2∑AcωAc2rAc(43)

张量 G\mathbf GG：

G=∑AcωAc2rAc2(44)\mathbf G=\sum_{A_c}\omega_{A_c}^2\mathbf r_{A_c}^2\tag{44}G=Ac∑ωAc2rAc2(44)

ωAc\omega_{A_c}ωAc：权函数，为控制点 Ac 到顶点 Ai 距离的倒数。

2.4 曲面导数的计算

式（30)、(31）中的梯度项和散度项均采用高斯积分法进行计算。若已知张量 ϕ\bf \phiϕ，其定义在曲面 S上，以封闭曲线 ∂S\partial S∂S 为边界，则：

∫S∇SϕdS=∮∂SmϕdL−∫SκnϕdS(45)\int_S\nabla_S\phi\;dS=\oint_{\partial S}\mathbf m\phi\; dL-\int_S\kappa\mathbf n\phi \;dS\tag{45}∫S∇SϕdS=∮∂SmϕdL−∫SκnϕdS(45)

n\bf nn：面 S 的单位法向量；
m\bf mm：单位副法向量(bi-normal)，其垂直于线 ∂S\partial S∂S ，并与面 S 相切；

如图，对于对控制区域 SAfS_{Af}SAf，采用有限体积法对面心速度散度 (∇S⋅v)Af(\nabla_S\cdot\mathbf v)_{Af}(∇S⋅v)Af 和面心速度梯度 (∇Sv)Af(\nabla_S\mathbf v)_{Af}(∇Sv)Af 进行离散：

(∇S⋅v)Af=1SAf∑emeveLe−κAfnAf⋅vAf(46)(\nabla_S\cdot\mathbf v)_{Af}=\frac{1}{S_{Af}}\sum_e\mathbf m_e\mathbf v_e L_e-\kappa_{Af}\mathbf n_{Af}\cdot\mathbf v_{Af}\tag{46}(∇S⋅v)Af=SAf1e∑meveLe−κAfnAf⋅vAf(46)

(∇Sv)Af=1SAf∑emeveLe−κAfnAfvAf(47)(\nabla_S\mathbf v)_{Af}=\frac{1}{S_{Af}}\sum_e\mathbf m_e\mathbf v_e L_e-\kappa_{Af}\mathbf n_{Af}\mathbf v_{Af}\tag{47}(∇Sv)Af=SAf1e∑meveLe−κAfnAfvAf(47)

下角标 e ：边界 Le 的中点值；
对封闭曲面 Af 的所有边的中心值进行加和。

注：控制区域 SAfS_{Af}SAf 的面积分、长度为 LeL_eLe 的控制区域边缘 e 的线积分运用中点规则近似。

各边中心的单位副法向量 me\mathbf m_eme：

me=e^×ni+nj∣ni+nj∣(48)\mathbf m_e=\hat \mathbf e \times\frac{\mathbf n_i+\mathbf n_j}{|\mathbf n_i+\mathbf n_j|}\tag{48}me=e^×∣ni+nj∣ni+nj(48)

e^\hat\mathbf ee^：与边 e 平行的的单位向量；
ni、nj\mathbf n_i、\mathbf n_jni、nj：边 e 上点 i、j 的界面单位法向量；

上图采用了基于边的局部坐标系，坐标轴分别与单位向量 n、t、t′\bf n、\bf t、\bf t^{'}n、t、t′ 对齐。其中，t\bf tt 与 PeN‾\overline{PeN}PeN 相切。

由此，边 e 中心处的速度值线性插值形式如下：
ve=(Te)T⋅[exTP⋅vP+(1−ex)TN⋅vN](49)\mathbf v_e=(T_e)^T\cdot[e_x\mathbf T_P\cdot\mathbf v_P+(1-e_x)\mathbf T_N\cdot\mathbf v_N]\tag{49}ve=(Te)T⋅[exTP⋅vP+(1−ex)TN⋅vN](49)

exe_xex：插值因子，

ex=eN‾PeN‾(50)e_x=\frac{\overline{eN}}{\overline{PeN}}\tag{50}ex=PeNeN(50)

TP、TN、Te\mathbf T_P、\mathbf T_N、\mathbf T_eTP、TN、Te：全局变换的张量，由直角坐标系转化为局部正交坐标系得来。

2.5 表面张力

追踪多相流的相界面过程中，必须对表面张力加以计算。如果表面张力计算不准确，会在界面附近产生非物理的寄生电流(parasitic current)。这种寄生电流会干扰乃至破坏界面和计算过程。

根据第一原理，在一个封闭的曲面上，表面张力的总和一定等于零。这是推导计算表面张力的重要依据。

对于任意的多边形控制区，采用非结构性表面网格对相界面进行离散，则作用在控制区域 SAfS_{Af}SAf 上的表面张力为：
FAfσ=∮∂SAfmσdL=∑e∫LemσdL=∑e(σm)eLe(51)\mathbf F_{Af}^\sigma=\oint_{\partial S_{Af}}\mathbf m\sigma \;dL=\sum_e\int_{Le}\mathbf m\sigma\;dL=\sum_e(\sigma \mathbf m)_eL_e\tag{51}FAfσ=∮∂SAfmσdL=e∑∫LemσdL=e∑(σm)eLe(51)

LeL_eLe：边 e 的长度；
(σm)e(\sigma\mathbf m)_e(σm)e：单位边长的表面张力。

对于共享边界 e 的两个控制区，如果 (σm)e(\sigma\mathbf m)_e(σm)e 大小相等且方向相反，则表面张力的总和为零。

将表面张力分解为切向 (∇Sσ)Af(\nabla_S\sigma)_{Af}(∇Sσ)Af 和法向 (κσ)AfnAf(\kappa\sigma)_{Af}\mathbf n_{Af}(κσ)AfnAf 两部分。根据高斯积分理论，有：

FAfσ=∫SAf∇SσdS+∫SAfκσndS(52)\mathbf F_{Af}^\sigma=\int_{S_{Af}}\nabla_S\sigma\;dS+\int_{S_{Af}}\kappa\sigma\mathbf n\;dS\tag{52}FAfσ=∫SAf∇SσdS+∫SAfκσndS(52)

式（52）RHS 采用中点规则进行离散，式（51），可得：

(∇Sσ)Af+(κσ)AfnAf=1SAf∑e(σm)eLe(53)(\nabla_S\sigma)_{Af}+(\kappa\sigma)_{Af}\mathbf n_{Af}=\frac 1{S_{Af}}\sum_e(\sigma \mathbf m)_eL_e\tag{53}(∇Sσ)Af+(κσ)AfnAf=SAf1e∑(σm)eLe(53)

因此，表面张力的切向部分等于式（53）RHS 的切向分量：

(∇Sσ)Af=1SAf(1−nAfnAf)⋅∑e(σm)eLe(54)(\nabla_S\sigma)_{Af}=\frac 1{S_{Af}}(1-\mathbf n_{Af}\mathbf n_{Af})\cdot\sum_e(\sigma \mathbf m)_eL_e\tag{54}(∇Sσ)Af=SAf1(1−nAfnAf)⋅e∑(σm)eLe(54)

同理，表面张力的法向分量表达式为：
(κσ)AfnAf=1SAf(nAfnAf)⋅∑e(σm)eLe(55)(\kappa\sigma)_{Af}\mathbf n_{Af}=\frac 1{S_{Af}}(\mathbf n_{Af}\mathbf n_{Af})\cdot\sum_e(\sigma \mathbf m)_eL_e\tag{55}(κσ)AfnAf=SAf1(nAfnAf)⋅e∑(σm)eLe(55)

当表面张量系数为常数，即 σ=constant\sigma=constantσ=constant，若控制区域 SAfS_{Af}SAf 的单位法向量满足：

n=∑emeLe∣∑emeLe∣(56)\mathbf n=\frac{\sum_e\mathbf m_e L_e}{|\sum_e\mathbf m_e L_e|}\tag{56}n=∣∑emeLe∣∑emeLe(56)

则式（54）所代表的表面张力的切向分量为零。

目前，多边形的法线采用如下方式计算：将多边形分解成三角形，对每个三角形的法线进行平均。

为了保证表面张力的总和为零，即使在表面张力系数为常数的情况下，也要对式（54）加以分析。

式（54)、(55）并不能完全保证表面张力总和为零。尤其当表面张力局部计算不准确，界面附近会出现非物理流体流动，导致表面张力的计算误差进一步加大。

通过式（51）可以看出：表面张的计算主要取决于单位边长的表面张力 (σm)e(\sigma\mathbf m)_e(σm)e。使用梯形规则(trapezoidal rule)来代替中点规则进行积分逼近：

(σm)e=σie^×ni+σje^×nj2(57)(\sigma\mathbf m)_e= \frac{\sigma_i\hat\mathbf e\times\mathbf n_i+\sigma_j\hat\mathbf e\times\mathbf n_j}{2}\tag{57}(σm)e=2σie^×ni+σje^×nj(57)

σi、σj\sigma_i、\sigma_jσi、σj：点 i、j 处的表面张力系数；

结果表明：使用梯形规则进行积分逼近可有效地减少表面张力的震荡。

通过式（57）可以看出表面张力的计算精度取决于控制区域顶点单位法向量的精度。局部坐标系下，采用最小二乘双二次曲面(least squares biquadratic surface)拟合计算顶点法线。局部坐标系的原点与 i 点重合，z 轴与顶点 i 的发现平行。这一双曲面表达式如下：

z(x,y)=ax2+by2+cxy+dx+ey(58)z(x,y)=ax^2+by^2+cxy+dx+ey\tag{58}z(x,y)=ax2+by2+cxy+dx+ey(58)

该曲面通过点 i ，并能够拟合所有经过改点的曲面。

假定 σ=c\sigma=cσ=c，根据式（55）计算平均曲率 κ\kappaκ 。

3. 解题步骤

该解题步骤基于PISO算法

更新时间步，根据上一时间步长的相应值，初始化当前时间步（new time step）的速度场、压力场以及质量通量（根据上一时间步长界面网格点的位置确定网格的相界面体积流，进而初始化相界面的净质量流）。
定义界面处网格点和控制点的位移方向。
进行外迭代循环：
a. 计算界面处网格点的位移，以保证净质量通量为零。（参见 2.3.3）
b. 以界面处网格点的位移作为边界条件，用于求解网格运动问题。网格运动后，利用式（28）重新计算相界面的体积流。
c. 更新压力场与速度场的边界条件。
d. 根据当前相界面的形状，组装并求解离散后的动量方程（8)。压力场与面质量流采用外迭代求解。
e. 开始 PISO 迭代循环：
（i）根据上一步得到的速度场，组装并求解离散后的压力方程（19)。
（ii) 根据当前压力场，利用方程（20）求解通过网格表面的绝对质量通量。
（iii) 根据当前的压力场，运用式（16）求解网格体心速度。
（iv) 如果 PISO迭代次数少于指定值，则返回步骤(i)。
f. 计算通过相界面的净质量通量。
g. 检查是否收敛，如果残差或者表面净质量通量没有满足精度要求，则返回步骤 a。
如果没有达到结束时间（end time），则返回步骤 1。

在步骤 b 中，如果只有界面处网格发生移动，上述解题步骤的效率将大大提高。在这种情况下，需要在外迭代开始前，根据上一时间步长得到的界面总位移，对整体网格进行变形。基于此，每个时间步长只需进行一次网格运动。这种方法适用于相界面位移小于网格厚度的情况。

根据PISO算法以及时间精度的需求，时间步长受到库朗特数（Courant number）的影响。在表面张力很大的情况下，由于表面张力的显式处理，时间步长受到更严格的限制，其必须小于最短的毛细波（capillary waves）周期：

Δt<ρmin(LPeN)2πσ(59)\Delta t<\frac{\sqrt{\rho min(L_{PeN})}}{2\pi\sigma}\tag{59}Δt<2πσρmin(LPeN)(59)

LPeNL_{PeN}LPeN：大地距离（geodetic distance） PeN‾\overline{PeN}PeN。

值得注意的是，对于预先规定了边界运动的移动边界问题，上述的求解过程可以直接被修正。

4. 质量守恒和迪利克雷流体边界条件

不可压缩流体的质量守恒条件取决于边界处净质量通量的数值。经验表明，当界面处网格分辨率足够高、外迭代次数足够多。界面处的净质量通量可以达到一个令人满意的水平。

当不可压缩流体 B 被规定速度边界（Dirichlet boundary）完全包围时，随着迭代过程中净质量通量的逐渐减少，压力方程的求解会出现问题。对于隶属于相 B 的边界面，规定其净质量通量为零，即：

m˙Bfn=ρBV˙Bfn(60)\dot m_{Bf}^n=\rho_B\dot V_{Bf}^n\tag{60}m˙Bfn=ρBV˙Bfn(60)

在相界面 B 一侧，式（19）中的项 nfn⋅(H/a)f\mathbf n_f^n\cdot(H/a)_fnfn⋅(H/a)f 表达式如下：

nBfn⋅(Ha)Bf=m˙BfnρBfnSBfn+(1a)BfnBfn⋅(∇p)Bfn(61)\mathbf n_{Bf}^n\cdot(\frac Ha)_{Bf}=\frac{\dot m_{Bf}^n}{\rho_{Bf}^nS_{Bf}^n}+(\frac 1a)_{Bf}\mathbf n_{Bf}^n\cdot(\nabla p)_{Bf}^n\tag{61}nBfn⋅(aH)Bf=ρBfnSBfnm˙Bfn+(a1)BfnBfn⋅(∇p)Bfn(61)

动压的法向导数 nBfn⋅(∇p)Bfn\mathbf n_{Bf}^n\cdot(\nabla p)_{Bf}^nnBfn⋅(∇p)Bfn 可由式（35）求得。

对于不可压缩流 B，其相界面的总体质量通量 ∑Bfm˙Bf\sum_{Bf}\dot m_{Bf}∑Bfm˙Bf 应恒为零，然而外迭代过程中并不能满足这一条件。这里，对其进行修正：

m˙Bfc=m˙Bfp−ωBf∑Bfm˙Bfp(62)\dot m_{Bf}^c=\dot m_{Bf}^p-\omega_{Bf}\sum_{Bf}\dot m_{Bf}^p\tag{62}m˙Bfc=m˙Bfp−ωBfBf∑m˙Bfp(62)

m˙Bfp\dot m_{Bf}^pm˙Bfp：前一步外迭代得到的质量通量；
ωBf\omega_{Bf}ωBf：权重因子，用于保证质量通量修正与 B 侧的质量通量分布成正比。

ωBf=∣m˙Bfp∣∑Bf∣m˙Bfp∣(63)\omega_{Bf}=\frac{|\dot m_{Bf}^p|}{\sum_{Bf}|\dot m_{Bf}^p|}\tag{63}ωBf=∑Bf∣m˙Bfp∣∣m˙Bfp∣(63)

随着外迭代的进行，质量通量收敛到零。