把Δt放大5倍快了1.9倍，放大10倍却原地踏步 — 隐式表面张力打开的窗口

一个液滴算了三天

一个二维振荡液滴已经算到第三天。流速不快，网格也不大。可时间步是 $10^{-6}$ 秒。表面张力一旦按显式处理， 决定时间步的就不再是流速，而是毛细波。

这时常见的建议是把表面张力改成隐式(用新时刻的值一起求解)。突破约束后时间步据说可以放大5倍、10倍。那么三天 会变成一天吗？

答案是"到5倍为止"。Janodet、van Wachem和Denner在2025年发表的全耦合算法论文，在密度比1000下把这个窗口的 两端都实测了一遍。上面是稳定极限，下面是每步的代价。本文沿着这两堵墙各自站在哪里往下走，并说明为什么加密 网格之后误差会停住。

捆住时间步的是毛细波，不是流速

界面上存在表面张力时，网格能分辨的最短毛细波只有一条，波长是 $\lambda_\sigma = 2\Delta x$。如果时间步大于 这条波穿过一个网格的时间，显式表面张力项就会发散。Denner和van Wachem整理的形式如下。

$\rho_A, \rho_B$ 是两种流体的密度，$\sigma$ 是表面张力系数，$\Delta x$ 是网格尺度。问题出在指数 $3/2$ 上。 网格减半，时间步缩小2.8倍，比对流CFL条件的 $\Delta x^1$ 收得更快。扩散项可以隐式求解并移出约束，表面张力 却长期做不到。这个约束从何而来、如何隐式化，写在 讨论毛细时间步约束的那篇里。

在下面的模拟中亲手调一调。

这是密度比1000的两种流体之间的毛细波在黏性作用下衰减的过程。灰色虚线是Prosperetti的解析解，蓝线是离散 求解器给出的振幅。提高 lambda/dx，两条线会靠拢。在 curvature dx^0.5 打开的情况下把 dt/dt_sigma 从0.5推到8，会发现误差几乎不动，这正是隔一节要讲的主题。

第一堵墙倒下的地方还立着第二个上限

把表面张力写成隐式确实可以越过 $\Delta t_\sigma$，但并不意味着时间步可以任意大。沿着Galusinski和Vigneaux 的分析，Denner等人给出的上限来自两个时间尺度的竞争。

其中 $\tau_{vc} = \hat\mu \lambda_\sigma / \sigma$ 是黏性-毛细时间尺度， $\tau_\sigma = \sqrt{\hat\rho \lambda_\sigma^3 / \sigma}$ 是毛细时间尺度 ($\hat\rho = \rho_A + \rho_B$，$\hat\mu = \mu_A + \mu_B$)。$a_1, a_2$ 是随算例变化的常数；取 $a_1 = 1/(16\pi)$ 且 $a_2 = 0$ 就恰好还原出原来的 $\Delta t_\sigma$。

两个尺度之比就是网格奥内佐格数。

$\mathrm{Oh}_{\Delta x} \ll 1$ 时惯性主导，$\Delta t^{*} \propto \tau_\sigma$；反之黏性主导， $\Delta t^{*} \propto \tau_{vc}$。工程上真正有利的是后者：黏度大或毛细波短，上限就开得很大。

数字上刺眼的是密度比。在静止液滴(Laplace平衡)算例中，$\mathrm{Oh}_{\Delta x} \ll 1$ 区域的上限在密度比 1000时是 $1.5\,\Delta t_\sigma$。同一系列算法在密度比1时可以用到 $15\,\Delta t_\sigma$，也就是十倍不见了。 在大 $\mathrm{Oh}_{\Delta x}$ 区域同样相差一个量级。真实的气液密度比会把窗口压窄。

网格加密8倍，误差只减了一半#

论文的第二个验证算例是衰减的毛细波。密度比和黏度比都取1000，拉普拉斯数 $\mathrm{La} = \rho\lambda\sigma/\mu^2 = 300$，网格 $\lambda/\Delta x = \{25, 50, 100, 200\}$，时间步 $\Delta t/\Delta t_\sigma = \{0.5, 2, 8\}$。与解析解的差别用振幅的 $L_2$ 范数衡量。

结果表里扎眼的不是误差大小，而是收敛阶。绝大多数格子落在0.46到0.95之间。同样的问题在密度比1下是出二阶的 地方。网格加密8倍，误差只减到大约一半。

论文指出的原因不是时间离散，而是界面输运，逻辑只有两行。所用的界面捕捉格式最好也就是二阶精度。曲率是颜色 函数的二阶导数，精度要掉两阶。也就是说曲率最好也只是零阶。在足够细的网格上，振幅误差的收敛阶最终趋于 零：误差停在一个常数上，不再下降。

用Python确认收敛阶的主人#

这条逻辑可以缩成一个阻尼振子来验证。线性范围内毛细波振幅满足 $A'' + 2\nu k^2 A' + \omega_0^2 A = 0$。求解器看到的不是 $\omega_0$，而是掺进曲率误差的 $\omega_{num} = \omega_0\sqrt{1 + C(\Delta x/\lambda)^q}$。只改变 $q$，用梯形法推进，取出与解析解之间的 $L_2$ 范数和收敛阶。

import math
 
SIGMA, RHO_HAT, LAMBDA, K, LA = 1.0, 1.0, 2*math.pi, 1.0, 300.0
MU = math.sqrt(RHO_HAT * LAMBDA * SIGMA / LA)
NU = MU / RHO_HAT
A0, T_END = LAMBDA / 100.0, 25.0
 
 
def capillary_omega(dx, q, c_kappa=0.6):
    """离散求解器实际看到的频率 — 曲率误差为 O(dx^q)"""
    w0 = math.sqrt(SIGMA * K**3 / RHO_HAT)
    return w0 * math.sqrt(1.0 + c_kappa * (dx / LAMBDA) ** q)
 
 
def analytic_amplitude(t):
    """A'' + 2*nu*k^2*A' + w0^2*A = 0 的解析解"""
    w0 = math.sqrt(SIGMA * K**3 / RHO_HAT)
    g = NU * K**2
    wd = math.sqrt(w0**2 - g**2)
    return A0 * math.exp(-g*t) * (math.cos(wd*t) + g/wd * math.sin(wd*t))
 
 
def march_amplitude(dt, w, n_steps):
    """用梯形法(Crank-Nicolson)推进 [A, A']"""
    g = NU * K**2
    a, v, hist = A0, 0.0, [A0]
    for _ in range(n_steps):
        h = 0.5 * dt
        rhs_a, rhs_v = a + h*v, v + h*(-w**2 * a - 2*g*v)
        det = (1 + 2*g*h) + h*h*w**2
        a = ((1 + 2*g*h) * rhs_a + h * rhs_v) / det
        v = (-h * w**2 * rhs_a + rhs_v) / det
        hist.append(a)
    return hist
 
 
def l2_amplitude(hist, dt):
    """振幅的 L2 误差范数 (论文 Eq. 61)"""
    acc = 0.0
    for i, a in enumerate(hist):
        w = 0.5 if i in (0, len(hist)-1) else 1.0
        acc += w * (a - analytic_amplitude(i*dt))**2 * dt
    return math.sqrt(acc / (len(hist)-1) / dt) / A0
 
 
def order_of(e_coarse, e_fine):
    return math.log(e_coarse / e_fine) / math.log(2.0)
 
 
for label, q in [("curvature error ~ dx^2", 2.0), ("curvature error ~ dx^0.5", 0.5)]:
    print(f"\n{label}")
    print("lam/dx |  dt/dt_s=0.5        dt/dt_s=2          dt/dt_s=8")
    prev = {}
    for n in [25, 50, 100, 200]:
        dx = LAMBDA / n
        dt_sigma = math.sqrt(RHO_HAT * dx**3 / (2*math.pi*SIGMA))
        w = capillary_omega(dx, q)
        row = []
        for s in [0.5, 2.0, 8.0]:
            dt = s * dt_sigma
            e = l2_amplitude(march_amplitude(dt, w, int(T_END/dt)), dt)
            tag = "  (-- )" if s not in prev else f" ({order_of(prev[s], e):4.2f})"
            row.append(f"{e:.3e}{tag}")
            prev[s] = e
        print(f"{n:6d} | " + "  ".join(row))

curvature error ~ dx^2
lam/dx |  dt/dt_s=0.5        dt/dt_s=2          dt/dt_s=8
    25 | 5.763e-04  (-- )  5.692e-04  (-- )  1.680e-02  (-- )
    50 | 1.516e-04 (1.93)  6.731e-05 (3.08)  2.023e-03 (3.05)
   100 | 3.885e-05 (1.96)  2.548e-05 (1.40)  2.345e-04 (3.11)
   200 | 9.833e-06 (1.98)  8.085e-06 (1.66)  2.506e-05 (3.23)
 
curvature error ~ dx^0.5
lam/dx |  dt/dt_s=0.5        dt/dt_s=2          dt/dt_s=8
    25 | 7.565e-02  (-- )  7.483e-02  (-- )  6.051e-02  (-- )
    50 | 5.449e-02 (0.47)  5.439e-02 (0.46)  5.266e-02 (0.20)
   100 | 3.897e-02 (0.48)  3.896e-02 (0.48)  3.874e-02 (0.44)
   200 | 2.775e-02 (0.49)  2.775e-02 (0.49)  2.773e-02 (0.48)

上表在时间步较小时给出二阶。$\Delta t/\Delta t_\sigma = 8$ 那一列的三阶不是额外收益：因为 $\Delta t \propto \Delta x^{3/2}$，二阶时间误差按 $\Delta x^3$ 下降。

下表才是论文描述的情形。收敛阶锁在0.5附近。更值得注意的是三列的数值几乎相同：时间步缩小16倍，误差纹丝 不动。决定精度下限的是曲率，不是时间离散。论文实测的0.46～0.95正好落在这张图上。

这一点与讨论界面对流CFL上限的那篇互为一对： 那边放大时间步后剩下CFL 0.05，这边加密网格后剩下曲率。

放大5倍得到1.9倍，放大10倍一无所获#

第三个算例是二维椭圆液滴的衰减振荡。长半轴0.15 m、短半轴0.1 m，以 $n=2$ 模态振荡，被黏性逐渐止住。实际 采用的时间步是两个约束的较小者。

$\Sigma$ 决定毛细约束被突破多少倍。论文跑了 $\Sigma \in \{2, 5, 10\}$，CFL上限固定为0.05。

先看精度：$\Sigma = 2$ 和 $5$ 的振荡频率误差约为3%，小于同样分辨率下显式表面张力处理得到的约4.5%。而 $\Sigma = 10$ 没能跟上动能衰减。论文的解释是，表面张力驱动的界面运动在时间上没有被分辨，因此不能指望时间 离散的形式二阶精度。

代价这一侧的数字就是本文标题。把 $\Sigma$ 从2提到5，即2.5倍，整体墙钟时间缩短了1.9倍。可继续提到10时，这份 收益没有保住，每步的折算计算时间显著上升。原因只有一个：时间步越大，每一步的非线性迭代收敛越慢。步数变少， 单步代价同比例上升。

S 滑块就是 $\Sigma$。从1推到5，蓝色赛道明显更早冲线；推到10几乎原地不动。把 Oh_dx 调低，红墙 ($\Delta t^{*}$)会左移，快车道干脆直接崩掉。

这些墙站在不同的位置

这一个算例上挂着的约束和最优点排成一张表是五行，性质各不相同。

突破 $\Delta t_\sigma$ 的算法只擦掉这张表的第二行，其余照旧。静止液滴上的力平衡如何维持，写在 讨论寄生流的那篇里。

那么 $\Sigma$ 该怎么选#

论文的结论是 $\Sigma$ 存在随算例而异的最优值，本问题中是5。找到它有三步。

先算 $\mathrm{Oh}_{\Delta x}$。这个数远小于1，稳定窗口本身就窄；密度比越大越窄。没有理由一上来就把 $\Sigma$ 放到10。

再数物理时间尺度：关心的振荡模态一个周期被切成了几步。稳定不等于准确，$\Sigma = 10$ 是稳定的，却漏掉了能量 衰减。

最后从日志里读每步的非线性迭代次数。若提高 $\Sigma$ 时迭代次数成比例增加，那一点就是窗口的右缘，墙钟时间 已经过了最低点。