液滴模拟为何被 Δt 绊住 — 把表面张力隐式求解的论文

1.5倍步长下显式峰值是3.6e+07，隐式峰值停在1.000。解开枷锁的是后向欧拉。

网格加密一倍，时间步长就缩到约2.8倍。物理并没有变快。带表面张力的界面上，毛细时间步长约束（capillary time-step constraint）几乎总是最紧的那条。

把 $\Delta t/\Delta t_\sigma$ 滑块推过0.9，只有橙色（显式）的振幅炸开。青色（隐式）还按着同一条驻波。

Δt / Δt_σ = 0.60mesh Oh = 0.10

explicit surface tension: stableimplicit: stable—— analytic decay

虚线是解析衰减。mesh Oh 一抬，阻尼变强，两条都更快平息。Janodet、van Wachem 与 Denner 在2025年的耦合，是把这块橙色炸开的地方放进一个矩阵。

评述对象是 R. Janodet, B. van Wachem, F. Denner, "A fully-coupled algorithm with implicit surface tension treatment for interfacial flows with large density ratios", arXiv:2410.17757 (2024)。

毛细时间步长这道枷锁

不可压界面流动求解器通常把扩散项隐式处理，抹去粘性时间约束。于是作为最强约束剩下的，就是表面张力。以 Brackbill 在 1992 年最早写下的形式，

\Delta t_\sigma = \sqrt{\frac{(\rho_A + \rho_B)\,\Delta x^3}{2\pi\sigma}} \propto \Delta x^{3/2}

$\rho_A, \rho_B$ 是两种流体的密度， $\sigma$ 是表面张力系数， $\Delta x$ 是网格间距。指数 $3/2$ 才是问题所在。为了把界面看得更清而把 $\Delta x$ 减半，时间步长就缩小 $2^{1.5}\approx 2.8$ 倍。液滴越小、流动越接近毛细尺度，步数就越爆炸。用显式表面张力，无路可逃。

为什么是 Δx^1.5 — 毛细波的色散关系#

这道约束其实是 CFL 条件。界面上的小正弦扰动以毛细波（capillary wave，以表面张力为恢复力的波）的形式振荡。在无粘深水极限下，色散关系为

\omega^2 = \frac{\sigma k^3}{\rho_A + \rho_B}

$\omega$ 是角频率， $k$ 是波数。网格能分辨的最短波是 $k_{\max} = \pi/\Delta x$ 。若 $\Delta t$ 超过该波的振荡周期，显式积分跟不上振荡而发散。计算 $\omega_{\max}\Delta t_\sigma$ 得到 $\sqrt{(k_{\max}\Delta x)^3/(2\pi)} = \pi/\sqrt{2}\approx 2.22$ 。也就是说， $\Delta t_\sigma$ 正是套在最快毛细波上的 CFL。

把表面张力隐式化 — 用一行振子看清核心

不必搬出整个求解器。一个毛细波模态可归约为阻尼谐振子。

\ddot{\eta} + 2\gamma\dot{\eta} + \omega^2 \eta = 0

$\eta$ 是界面振幅， $\gamma$ 是粘性阻尼。显式表面张力等同于用辛欧拉法积分这个振子，在 $\omega\Delta t > 2$ 处发散 —— 正是上面的毛细 CFL。把表面张力项取在新时刻的值（后向欧拉），就变成无条件稳定。

上面的橙色在0.9附近炸开，是因为辛欧拉的极限 $\omega\Delta t = 2$ 除以输出里的 $\omega\Delta t_\sigma = 2.221$ 。 $2/2.221 \approx 0.900$ 。再往上，显式一侧就把振幅放大。青色到3倍仍沿着虚线衰减。mesh Oh（粘性尺度）一大，阻尼变强，两条都更快平息。

静止界面上留下来的假流动，是这股恢复力的离散失去平衡时出现的。那个位置是力平衡被打破的寄生电流。

透过放大因子看稳定边界

用振子的放大因子（amplification factor，每步振幅被乘的倍率） $|G|$ ，边界就清晰了。numpy 三十行即可复现。玩具问题不是激波管，而是驻立毛细波的衰减。

import numpy as np
 
def capillary_omega(sigma, k, rho_hat):
    # 色散关系: 无粘深水毛细波 omega^2 = sigma k^3 / (rhoA+rhoB)
    return np.sqrt(sigma * k**3 / rho_hat)
 
def dt_capillary(sigma, dx, rho_hat):
    # 式(4): Brackbill 毛细时间步长 ~ dx^{3/2}
    return np.sqrt(rho_hat * dx**3 / (2*np.pi*sigma))
 
def run_standing_wave(omega, gamma, dt, steps, implicit):
    eta, v = 1.0, 0.0
    peak = 0.0
    for _ in range(steps):
        if implicit:        # 后向欧拉: 表面张力取新时刻值 (无条件稳定)
            v = (v - dt*omega**2*eta) / (1 + dt**2*omega**2 + 2*gamma*dt)
            eta = eta + dt*v
        else:               # 辛欧拉: 显式表面张力 (omega*dt<2 约束)
            v = v - dt*omega**2*eta - dt*2*gamma*v
            eta = eta + dt*v
        peak = max(peak, abs(eta))
    return peak
 
sigma, rho_hat, L = 0.07, 1001.2, 1.0      # 水/空气
for N in (64, 128, 256):
    dx = L/N; k = np.pi/dx                  # 最短可分辨毛细波
    om = capillary_omega(sigma, k, rho_hat)
    dts = dt_capillary(sigma, dx, rho_hat)
    print(f"N={N:4d}  dt_sigma={dts:.2e}  omega*dt_sigma={om*dts:.3f}")
# 输出: 所有 N 下 omega*dt_sigma = 2.221 (= pi/sqrt(2))，CFL 即毛细约束本身
 
N = 128; dx = L/N; k = np.pi/dx
om = capillary_omega(sigma, k, rho_hat); dts = dt_capillary(sigma, dx, rho_hat)
dt = 1.5*dts; gamma = 0.1*om                # 强制取毛细约束的 1.5 倍
pe = run_standing_wave(om, gamma, dt, 200, implicit=False)
pi_ = run_standing_wave(om, gamma, dt, 200, implicit=True)
print(f"dt=1.5*dt_sigma  显式峰值={pe:.1e}  隐式峰值={pi_:.3f}")
# 输出示例: 显式峰值=3.6e+07  隐式峰值=1.000

代码放进去的是 $\sigma=0.07$ 、密度和 $1001.2$ 、长度 $1$ 。 $1001.2$ 是水加空气， $0.07$ 是表面张力。 $N=64$ 时 $\Delta x=1/64$ ，同一公式给出的 $\Delta t_\sigma$ 约 $9.32\times 10^{-2}$ 。 $N=128$ 再除以 $2^{3/2}\approx 2.828$ ，约 $3.29\times 10^{-2}$ 。 $N=256$ 再除一次，约 $1.16\times 10^{-2}$ 。从64到256，间距只变成四分之一，同一物理时间的步数却变成8倍，因为 $4^{3/2}=8$ 。omega*dt_sigma 与分辨率无关，总是2.221。注释里的例子取 $\Delta t=1.5\Delta t_\sigma$ 、200步、 $\gamma=0.1\omega$ ，显式峰值 $3.6\times 10^{7}$ ，隐式 $1.000$ 。那时 $\omega\Delta t=1.5\times 2.221=3.3315$ ，越过辛欧拉的2。

$|G|$ 曲线的横轴是 $\Delta t/\Delta t_\sigma$ 。mesh 和 $\sigma$ 滑块移动的是 $\Delta t_\sigma$ 本身。

Δt / Δt_σ = 0.90mesh L/Δx = 128σ = 0.07

explicit |G| = 1.000implicit |G| = 0.447Δt_σ ≈ 3.29e-2 (water/air, L=1)

显式（橙色）把 $|G|=1$ 守到 $\Delta t/\Delta t_\sigma \approx 0.9$ ，再往上就越过1。隐式（青色）在任何步长下都 $|G|<1$ 。把 $\sigma$ 提高4倍， $\Delta t_\sigma$ 变成 $\sqrt{1/4}=1/2$ 。这道枷锁的斜率并不只来自网格。

全耦合才是关键 — 把色函数、压力、速度装进一个矩阵

隐式化的难点在于，表面张力 $\mathbf{f}_\sigma = \sigma\kappa\nabla\psi$ 与色函数（colour function，用 0~1 标记流体种类的 VOF 变量） $\psi$ 非线性纠缠。曲率 $\kappa$ 也由 $\psi$ 的导数决定。论文对该项关于新时刻的 $\psi$ 做 Newton 线性化，并把界面输运方程与连续、动量方程一起作为单一线性系统 $\mathbf{A}\boldsymbol{\phi}=\mathbf{b}$ 求解，其中 $\boldsymbol{\phi}=(p, u, v, w, \psi)$ 被整体隐式耦合。Denner 等人此前的算法被锁死在密度比 1。Janodet 的贡献是把连续与动量写成守恒形式，并在瞬态项中把密度关于 $\psi$ 隐式处理，从而在水/空气 1000:1 这样的大密度比下仍保持耦合。界面输运用 THINC/QQ 而非 CICSAM 离散，在大 CFL 下也保持界面锐利。与 OpenFOAM 的分离式 PIMPLE 不同，这要在块耦合（block-coupled）求解器里才有意义。

复现时撞上的问题

“无条件稳定”是有附注的。即便隐式表面张力突破了毛细约束，可稳定的最大 $\Delta t$ 仍然有限。Galusinski–Vigneaux 的粘性-毛细上限给出了它。在密度比 1000 下，这个上限比密度比 1 时严苛一个数量级，论文自己也报告在无粘极限下只稳定到约 $1.5\,\Delta t_\sigma$ 。上面的振子模型简化了这种粘性-毛细耦合，因而让隐式一侧显得比现实更宽容。此外后向欧拉虽稳定，却会人为地衰减快毛细波。在需要能量守恒的问题里，这种数值衰减可能被误认为物理衰减，所以论文单独验证了力平衡与能量守恒。

停在 1.5 Δt_σ，还是丢掉密度比 1000#

隐式跨过毛细 CFL 之后，上限还在。无粘极限里论文说的稳定窗口大约是 $1.5\,\Delta t_\sigma$ 。注释里 1.5 倍、200 步的例子，显式峰值冲到 $3.6\times 10^{7}$ ，隐式停在 1.000，那是振子，不是求解器。密度比 1000 时，Galusinski–Vigneaux 上限比密度比 1 紧一个数量级。振子给的宽限，到求解器里留不下来。

后向欧拉用数值把快毛细波杀掉，换来稳定。能量必须守住的问题里，这种衰减很容易被当成物理。把 $\Delta t$ 再往前推之后的窗口，记在五倍处打开、十倍处关上的记录里。