Rhie-Chow 插值：为什么需要它，以及在多相流中为什么棘手#

Rhie-Chow 插值（动量加权插值，MWI）改的是同置网格上的面速度，用来压掉棋盘格压力振荡。单相只有一项。到了多相，压力、重力、表面张力、非定常修正在同一个面上叠在一起。错一项，界面上就留下假速度。

从 Rhie & Chow (1983) 的稳态 SIMPLE 出发，经过 Bartholomew 等人 (2018) 的统一式，再到 interFoam 和 compressibleInterFoam。

勾上 Rhie–Chow，交错的压力花纹熄掉；取消勾选，花纹马上回来。

不修正时，模板跳过 $p_i$，只剩 $p_{i-1}$ 和 $p_{i+1}$。修正打开后，紧致梯度把被跳过的那一格重新接上。

1. 从半离散动量方程推导面速度#

所有 Rhie-Chow 变体公式的起点都是单元中心 $P$ 点处的半离散（semi-discretized）动量方程：

• $a_P$：对角项系数（包含时间导数、对流、扩散的贡献）
• $H(\mathbf{U}) = \sum a_{Nb}\mathbf{U}_{Nb} + \text{sources}$：相邻单元的贡献 + 显式源项

求解单元中心速度：

问题：棋盘格振荡

如果直接通过简单线性插值将单元中心速度投影到面上，压力梯度只会连接 $p_{i-1}$ 和 $p_{i+1}$ 而跳过 $p_i$，从而导致奇偶解耦（棋盘格现象）。

Rhie-Chow 的解决方案#

在面上重构动量方程本身，用紧致（compact）的面中心压力梯度代替宽模板的插值压力梯度：

其中 $\nabla p\big|_f^{\text{compact}} = (p_Q - p_P)/|\mathbf{d}_{PQ}|$ 是直接利用相邻单元值计算的面梯度，上划线表示单元中心值的线性插值。

压力修正项 $-(1/a_P)_f[\nabla p|_f - \overline{\nabla p}_f]$ 起到了类似于 $\Delta x^2 \cdot \partial^3 p/\partial x^3$ 的低通滤波器作用，从而抑制棋盘格。

2. 核心论文的贡献#

Rhie & Chow (1983)：原型#

发表在 AIAA Journal 上的原文是专为稳态 SIMPLE 算法设计的：

$\hat{u} = u + d \cdot \Delta p$ 是去掉压力梯度的伪速度，$d = V/a_P$。它仅包含压力修正项，没有处理非定常项或松弛（relaxation）项。后来 Majumdar (1988) 指出了收敛解对松弛因子的依赖性，Choi (1999) 指出了非定常流中对时间步长的依赖性。压力打成格子纹时 SIMPLE 用 Rhie–Chow 接上的那一项，到了多相还要再分给重力和表面张力。

Jasak (1996)：OpenFOAM 的基石 - H/A 方法#

在其伦敦帝国理工学院的博士论文中，系统化了非结构多面体网格的有限体积法。其特点是隐式实现 Rhie-Chow，而没有显式的修正项：

HbyA = H/A 是去除了压力影响的速度预测值，压力方程推导为 $\nabla\cdot[(1/A)_f \nabla p] = \nabla\cdot(H/A)_f$。Rhie-Chow 效应自然产生，因为拉普拉斯离散（面 snGrad）和速度修正的梯度（高斯定理）使用了不同的模板。

非定常修正通过 fvc::ddtCorr 处理，但并非完全一致，因此仍存在时间步长依赖性。

Cubero & Fueyo (2007)：紧致动量插值 (CMI)#

提出了一个独立于时间步长和松弛因子的面速度公式。核心思想是分别且一致地插值动量方程的各个系数：

• $a_t = \rho V/\Delta t$：时间系数
• $a_\tau = \rho V/\Delta\tau$：松弛系数
• $a_s$：空间系数

收敛时，$u_f = u^{\text{prev}}_f$ 且 $u_f = u^O_f$，因此非定常/松弛修正项会自动消失。这保证了稳态解与 $\Delta t$ 及松弛因子无关。

Bartholomew et al. (2018)：统一框架#

在 JCP 375 中提出了适用于结构化和非结构化网格的统一 MWI 公式。三个核心贡献：

1) MWI 系数应仅由空间系数组成。 如果将时间导数系数 $\rho V/\Delta t$ 放在分母中，当 $\Delta t \to 0$ 时，压力滤波器会变弱，导致棋盘格再次出现。

2) “驱动压力梯度 (Driving pressure gradient)”概念： 当存在体心力时，MWI 滤波器应仅作用于 $\nabla P = \nabla p - \mathbf{S}_{\text{total}}$（减去所有源项后的驱动压力梯度）。这是平衡力（balanced-force）离散的核心。

3) 密度加权： 单元中心压力梯度按 $\nabla P_P/\rho_P$ 加权，面密度使用调和平均插值，从而在密度比高达 $10^{24}$ 时仍能获得稳定的结果。

Denner & van Wachem (2014)：平衡力 VOF#

通过在同一位置（面）使用相同的离散梯度算子评估 CSF 表面张力 ($\sigma\kappa\nabla\alpha$) 和压力梯度，在静态液滴中实现了达到求解器容差（tolerance）级的精确力平衡。该方法在密度比超过 $10^6$ 时依然稳定，并大幅减少了寄生速度（spurious velocity）。

3. 多相流中四种修正项的处理#

多相流面速度的四种修正项各有其物理意义和处理原则。

压力修正项（必选）

这是 Rhie-Chow 的本质。它以产生与 $\Delta x^2 \partial^4 p/\partial x^4$ 成正比的质量守恒误差为代价来抑制棋盘格。在多相流中，由于界面附近压力剧烈变化，误差可能增大，因此需要平衡力方法。

非定常修正项（非定常流必选）

如果没有这一项，当 $\Delta t \to 0$ 时 MWI 滤波器强度会发生变化，导致出现锯齿状压力振荡。Cubero & Fueyo (2007) 和 Bartholomew 等人 (2018) 建议将时间系数从 MWI 系数中完全分离出来，作为独立的修正项处理。

重力修正项（多相流必选）

在存在密度不连续的多相流中必须包含此项。Mencinger & Zun (2007) 证明了如果没有此项，界面处会产生非物理速度尖峰。在静压平衡 ($\mathbf{U}=0$) 下，$\nabla p = \rho\mathbf{g}$ 在离散意义上必须成立。

在 OpenFOAM 中，使用修正压力 $p_{rgh} = p - \rho\mathbf{g}\cdot\mathbf{r}$，并在面上直接以 $-g_h \cdot \text{snGrad}(\rho)$ 的形式进行评估。

表面张力修正项（最具争议）

在平衡力方法中，不应对表面张力应用 Rhie-Chow 修正。原因显而易见：

“

如果压力梯度和表面张力使用相同的面中心梯度算子进行离散化， 那么在静态平衡时 $\nabla p = \sigma\kappa\nabla\alpha$ 精确成立。 此时添加 Rhie-Chow 修正会破坏平衡，导致寄生电流（parasitic current）。

因此，面上的表面张力通过 $(1/a_P)_f \cdot (\sigma\kappa)_f \cdot \text{snGrad}(\alpha)$ 直接评估，而与插值的单元中心值之间的差异修正则设为 0。

修正项总结

4. 抑制寄生速度的设计原则#

寄生电流（寄生速度，spurious velocity）是在具有界面曲率的多相流中，由于表面张力和压力梯度的离散不平衡而产生的非物理速度场。其大小与 $\sigma$ 成正比，与 $\mu$ 成反比，在低粘度流动中尤为严重。

平衡力算法

Francois 等人 (2006) 提出的核心原则：

“

压力梯度和表面张力的梯度算子必须离散化方式相同，且评估位置相同。

具体来说，由于压力方程中的 $\nabla p$ 是用面 snGrad 计算的，因此 CSF 的 $\sigma\kappa\nabla\alpha$ 也必须在面上使用相同的 snGrad 计算。遵循这一原则，在给定精确曲率的静态液滴中，可以实现达到机器精度级的零速度。

Rhie-Chow 本身也可能是原因#

Mencinger & Zun (2007) 的重要发现：当体心力不连续（密度突变、表面张力）时，压力的四阶导数会变大，Rhie-Chow 修正项本身就会产生巨大误差。

因此，Bartholomew 等人 (2018) 的统一框架不对体心力应用 Rhie-Chow 修正，而仅对驱动压力梯度应用滤波器：

尽管如此，曲率估计的准确度仍决定了寄生速度的极限，高度函数法 (height-function, Popinet, 2009) 提供了二阶收敛性，最为有效。

5. 迁移率系数插值：为什么要用调和平均#

面上迁移率系数 $(1/a_P)_f$ 的插值对多相流的稳定性有很大影响。由于 $a_P$ 正比于密度 ($a_P \propto \rho$)，在大密度比流动中，界面两侧的 $1/a_P$ 值可能相差数个数量级。

算术平均的问题

在水-空气 ($\rho_w/\rho_a = 1000$) 界面处：

结果被轻相（空气）较大的 $1/a_P$ 值主导，从而低估了重相（水）的惯性。

调和平均的解决方案

对于水-空气，$\rho_f^* \approx 2.0$，这与算术平均的 $\approx 500.5$ 截然不同。通过妥善反映对加速度具有更大阻力的相（重相），Bartholomew 等人 (2018) 和 Denner 等人 (2022) 在密度比超过 $10^6$ 时仍获得了稳定的结果。

如果使用算术插值，压力方程的拉普拉斯系数与 MWI 修正的密度加权之间会出现不匹配，从而产生寄生速度。

500.5 从哪来，写清楚就行。水的密度 1000，空气的密度 1，$(1000+1)/2 = 500.5$。调和平均是 2 除以 $(1/1000 + 1)$。分母是 1.001，$2/1.001 \approx 1.998$，就是文中的 2.0。500.5 除以 2.0，大约 250。轻相把平均拉高了这么多倍。

$a_P$ 与密度成正比，$1/a_P$ 就按同样的比例拉开。空气侧的 $1/a_P$ 是水侧的 1000 倍。算术平均把那个大数的一半安到面上。水的惯性，也就是小的 $1/a_P$，几乎被抹掉。密度比 $10^6$ 时，这个差距是一百万倍。调和平均被写成能撑到 $10^{24}$。Denner 的平衡力在 $10^6$ 以上把静止液滴的力平衡收到求解器容差。

滤波本身是量级为 $\Delta x^2 \partial^3 p/\partial x^3$ 的低通。界面那种压力四阶导数很大的地方，这项修正会把质量误差 $\Delta x^2 \partial^4 p/\partial x^4$ 放大。所以体积力留在滤波外面。把时间系数 $\rho V/\Delta t$ 放进 MWI 分母，$\Delta t$ 走向 0 时滤波变弱，锯齿状的压力振荡会回来。分母里只留空间系数。

6. OpenFOAM 实现分析#

interFoam：不压缩双相流#

OpenFOAM 的 interFoam 基于 Jasak 的 H/A 方法。

动量方程组装（不含压力梯度）：

fvVectorMatrix UEqn(
    fvm::ddt(rho, U)
  + fvm::div(rhoPhi, U)
  + turbulence->divDevRhoReff(rho, U)
);

计算 rAU 和 HbyA：

volScalarField rAU("rAU", 1.0/UEqn.A());
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU)); // 线性插值（默认值）
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));

构建 phiHbyA（4 种修正项）：

// (H/A)_f · S_f + 非定常修正
surfaceScalarField phiHbyA("phiHbyA",
    fvc::flux(HbyA)
  + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi, Uf)
);
 
// 重力 + 表面张力：直接在面上评估
surfaceScalarField phig(
    (
        mixture.surfaceTensionForce()    // sigma*kappa*snGrad(alpha)
      - ghf*fvc::snGrad(rho)            // -g·r_f * snGrad(rho)
    ) * rAUf * mesh.magSf()
);
 
phiHbyA += phig;

核心要点：

• ghf = g & mesh.Cf()：直接在面中心评估的重力势
• mixture.surfaceTensionForce()：在面上计算 $\sigma\kappa\cdot\text{snGrad}(\alpha)$
• 两种体心力均直接在面上评估，近似遵循平衡力原则

压力方程与通量修正：

fvScalarMatrix p_rghEqn(
    fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
);
p_rghEqn.solve();
phi = phiHbyA - p_rghEqn.flux();

由于 pEqn.flux() 是 $-(1/A)_f \cdot \text{snGrad}(p_{rgh}) \cdot |S_f|$，最终通量是无压力速度 + 面评估体心力 + 面评估压力梯度的组合。拉普拉斯算子使用面 snGrad 而速度修正梯度使用单元面值，这种模板不对称性隐式地提供了 Rhie-Chow 效应。

ddtCorr 的局限性： fvc::ddtCorr(U, phi) 计算的是 $(1/\Delta t)\cdot[\phi_{\text{old}} - (\mathbf{U}_{\text{old}})_f\cdot\mathbf{S}_f]$。与 CMI 不同，它并非完全一致，因此仍存在时间步长依赖性，如 Komen 等人 (2020) 所证，这可能会影响时间精确度。

compressibleInterFoam：压缩双相流#

phiHbyA 的构建与 interFoam 相同，但压力方程由三部分组成：

在不压缩部分（拉普拉斯项）中加入了各相的压缩性项。通过 phid = fvc::interpolate(psi)*phi ($\psi = d\rho/dp$) 实现压力-密度耦合，并在 $\alpha$ 方程中加入微分压缩性项 dgdt 以反映两相密度变化率的差异。

7. 方法对比#

一致性方法（CMI、统一 MWI）虽然比原版数值扩散更大，但保留了时间精确度，并提供了与网格/时间步长无关的解。残余寄生速度的主要原因是曲率估计误差。

时间系数进了分母，棋盘格就回来

滤波只加在驱动压力上。重力和表面张力再被同一套式子滤一遍，静止液滴上对好的平衡就破了，寄生电流跟着出现。

$\rho\mathbf{g}$ 和 $\sigma\kappa\nabla\alpha$ 留在面的 snGrad 上，MWI 只放行 $\nabla P$。rAU 仍用算术平均时，密度比 1000 的面系数会倒向空气。调和平均 $\rho_f^* \approx 2.0$ 把水的惯性留下来。ddtCorr 不动，$\Delta t$ 依赖就还在。simpleFoam 里压力残差死不掉的情形和这个修正错位的地方叠在一起。曲率不换之前，寄生速度的地板就是曲率误差。

参考文献

• Rhie, C.M. & Chow, W.L. (1983). AIAA Journal 21(11)
• Jasak, H. (1996). PhD Thesis, Imperial College London
• Cubero, A. & Fueyo, N. (2007). Numerical Heat Transfer Part B 52(6)
• Francois, M.M. et al. (2006). JCP 213(1)
• Mencinger, J. & Zun, I. (2007). JCP 224(1)
• Denner, F. & van Wachem, B. (2014). Numerical Heat Transfer Part B 65(3)
• Bartholomew, P. et al. (2018). JCP 375
• Aguerre, H.J. et al. (2018). JCP 365
• Komen, E.M.J. et al. (2020). Nuclear Engineering and Design 370