[论文评述] 毛细枷锁解开后，CFL 0.05 还在

压缩格式用来立起界面的那个值，是被 Courant 数除过的商。时间步长一大，最先消失的就是这个商。

Janodet 等人 2025 年论文的结论里有这样一句："所提出的算法可以用大于毛细时间步长约束的时间步长计算真实的气液流动 — 只要其他时间步长约束得到满足。"着重是这里加的。表面张力改成隐式，解开了一道枷锁。论文实际跑到的最大 CFL 仍是 0.05。输运色函数的格式把时间步长接了过去。这个上限在 NVD 方框的哪一条边上折下来，涡对流又向体积收多少账，用数字对一下。

毛细约束的公式后面会原样出现。这里只处理把那一项划掉之后剩下的 $C_{\max}$ 。

立起界面的那个值来自下风

代数 VOF（不重构界面、直接对流色函数的做法）里界面变厚只有一个原因：迎风面值总会把界面抹平。于是压缩格式把面值朝下风单元拉。极端地直接取下风值，阶跃就被关进一个单元里。

麻烦在于下风值不保证有界性（boundedness）。色函数一旦跌破 0 或越过 1，密度就变成负数，计算到此为止。所以需要一条规则来规定能拉到多远。那条规则以 Courant 数为自变量。

把 Courant C 从 0.05 推到 0.9，可行区域的天花板就朝迎风对角线折下去。

Courant C0.05

blend gamma_f0.60transition cells0slab width0

Start at C = 0.05: the amber window fills almost the whole box, the blue HYPER-C curve pins itself to phi~_f = 1, and the slab keeps a two-cell edge forever. Drag C toward 0.9 and the ceiling min(1, phi~_D/C) folds down onto the dashed upwind diagonal — the window is a sliver, the pink face states have nowhere to sit, and the slab bleeds out over a dozen cells.

把 Courant C 放在 0.05，左侧琥珀色的可行区域几乎填满整个方框，右侧的方块保持两个单元宽的边缘。把 C 推到 0.9，天花板就塌到虚线对角线（即迎风）上。粉色的点无处落脚，方块随时间不断洇开。格式没变，网格没变，变大的只有时间步长。

在 NVD 方框上，Courant 数把天花板压低#

归一化变量图（NVD）以一个面为基准，把迎风单元 $U$ 、施主单元 $D$ 、受主单元 $A$ 的值这样规格化。

\tilde{\phi}_D = \frac{\phi_D - \phi_U}{\phi_A - \phi_U}, \qquad \tilde{\phi}_f = \frac{\phi_f - \phi_U}{\phi_A - \phi_U}

$\tilde{\phi}_D$ 表示施主单元位于迎风与受主之间的何处， $\tilde{\phi}_f$ 表示面值落在何处。 $\tilde{\phi}_f = \tilde{\phi}_D$ 是迎风， $\tilde{\phi}_f = 1$ 是下风。

Leonard 的对流有界性判据（CBC, Convection Boundedness Criterion）钉死了显式时间推进下面值必须待的位置。

\tilde{\phi}_D \le \tilde{\phi}_f \le \min\!\left(1,\ \frac{\tilde{\phi}_D}{C}\right), \qquad 0 \le \tilde{\phi}_D \le 1

其中 $C = u\,\Delta t/\Delta x$ 是该面的 Courant 数。天花板 $\tilde{\phi}_D/C$ 的含义很直接：一步之内从施主单元抽走体积 $C$ ，若这份体积里携带的色函数量超过单元原有的量，单元就变成负值。这个条件就是 $C \cdot \tilde{\phi}_f \le \tilde{\phi}_D$ ，整理后即上面的不等式。

$C$ 取 0.05 时天花板是 $20\,\tilde{\phi}_D$ 。 $\tilde{\phi}_D$ 刚过 0.05 就能把面值抬到 1。 $C$ 取 0.8 时天花板是 $1.25\,\tilde{\phi}_D$ ，只剩紧贴对角线的一条窄带。可用的压缩量正比于 $1/C$ 。

CICSAM 停在 0.01、THINC/QQ 停在 0.05 的原因#

CICSAM 在这个方框里混合两条曲线：一条是紧贴天花板的 HYPER-C，另一条是较平缓的 ULTIMATE-QUICKEST。混合权重 $\gamma_f$ 由界面法向与面向量的夹角决定。界面垂直于面就偏向 HYPER-C，倾斜就偏向 UQ。因为压缩倾斜界面会产生阶梯状的人为褶皱。

上面模拟中的 CICSAM blend 按钮和 blend gamma_f 滑块就是这个混合。调低 $\gamma_f$ ，曲线从天花板下滑，方块立刻变厚。同时也能看到，一维对齐界面上 $\gamma_f = 1$ ，CICSAM 实际上退化成了 HYPER-C。

CICSAM 的实用 CFL 上限落在 0.01 附近，原因正在于此。在 $\gamma_f$ 于 0 和 1 之间来回摆动的真实三维界面上，只有 HYPER-C 分量自身满足 CBC，UQ 分量必须另行再限制一次。要让混合结果不越过天花板， $C$ 就得足够小。论文为了避开这个约束，用 THINC/QQ 取代了 CICSAM。

一个 tanh 在单元内部重画界面#

THINC（Tangent of Hyperbola for INterface Capturing）不去挑选面值，而是干脆把单元内部的分布画出来。在把单元规格化到 $[0,1]$ 的坐标 $\tilde{x}$ 上，令

\Phi(\tilde{x}) = \frac{1}{2}\left[1 + \gamma \tanh\!\big(\beta(\tilde{x} - \tilde{x}_c)\big)\right]

其中 $\beta$ 是界面锐度（通常取 2 附近）， $\gamma = \pm 1$ 是从相邻单元读出的界面朝向， $\tilde{x}_c$ 是 tanh 跳变所在的位置。 $\tilde{x}_c$ 由精确复现单元平均 $\bar{\phi}$ 的要求确定，并且有闭式解。

\tilde{x}_c = \frac{1}{\beta}\,\mathrm{artanh}\!\left(\frac{\cosh\beta - e^{\beta\gamma(2\bar{\phi}-1)}}{\sinh\beta}\right)

穿过面的量由这条曲线在出发区域上的积分给出。

F_{i+1/2} = \int_{1-C}^{1} \Phi(\tilde{x})\,\mathrm{d}\tilde{x}

THINC/QQ 在此基础上再加一层二次曲面（quadratic surface）重构，对有曲率的界面抓得更准。在激波一侧使用同一个 tanh 的做法，见 TENO-THINC 重构。

这个积分把 $C$ 直接写进了积分限。 $C$ 变大，积分区间就逼近整个单元宽度，最终等同于只输运一个单元平均值。tanh 不使用下风值，所以 CBC 自动满足，但压缩力随 $C$ 衰减这一性质完全相同。

时间步长预算里有三个条目

现在看整份预算。显式求解毛细波的约束是

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

旁边还并列着流动的 CFL 约束 $\Delta t_{\text{CFL}} = C_{\max}\Delta x / U$ 。论文做的事，是把第一个条目从预算里划掉。剩下的是界面对流格式允许的 $C_{\max}$ 。

两条赛道的差距，在输出正后方的赛跑里再看一次。U 和 interface CFL cap 是那里的旋钮。

用一个涡量出的体积误差

$C$ 变大时究竟什么在变坏？在方向分裂（directional splitting）对流中，每个方向的扫掠看不到无散的速度场，因此必须加入膨胀修正项。

\phi^{*}_{i} = \phi^{n}_{i} - \left(F_{i+1/2} - F_{i-1/2}\right) + \phi^{n}_{i}\left(C_{i+1/2} - C_{i-1/2}\right)

这一项保证均匀的 $\phi = 1$ 区域不会被一次扫掠破坏。但在界面单元里， $\phi^n_i$ 与扫掠过程中的实际值并不相同，这个差额便以体积误差的形式留下。我把 THINC 对流放到 Rider–Kothe 单涡（含时间反转， $T=2$ ）上做了测量。

from math import atanh, cos, cosh, exp, log, log1p, pi, sin, sinh
 
N, BETA, EPS = 40, 2.0, 1e-6
H = 1.0 / N
 
 
def lncosh(z):
    a = abs(z)
    return a + log1p(exp(-2.0 * a)) - log(2.0)
 
 
def thinc_slab(pbar, g, a, b):
    """把施主单元的 tanh 重构在 [a, b] 区间上积分。"""
    s = g * (2.0 * pbar - 1.0)
    r = max(-0.999999, min(0.999999, (cosh(BETA) - exp(BETA * s)) / sinh(BETA)))
    xc = atanh(r) / BETA
    return 0.5 * ((b - a) + (g / BETA) * (lncosh(BETA * (b - xc)) - lncosh(BETA * (a - xc))))
 
 
def face_flux(pm, p0, pp, c):
    """穿过面的色函数量。p0 是施主单元，c 是该面的 Courant 数。"""
    if abs(c) < 1e-14:
        return 0.0
    g = 1.0 if pp > pm else (-1.0 if pp < pm else 0.0)
    if g == 0.0 or p0 < EPS or p0 > 1.0 - EPS:
        return c * p0
    return thinc_slab(p0, g, 1.0 - c, 1.0) if c > 0 else -thinc_slab(p0, g, 0.0, -c)
 
 
def line(col, vel, k):
    """周期边界上的一次一维扫掠，含膨胀修正项。"""
    n, out = len(col), [0.0] * len(col)
    for i in range(n):
        cw, ce = vel[i] * k, vel[i + 1] * k
        fw = face_flux(col[(i - 2) % n], col[(i - 1) % n], col[i], cw) if cw > 0 else \
            face_flux(col[(i - 1) % n], col[i], col[(i + 1) % n], cw)
        fe = face_flux(col[(i - 1) % n], col[i], col[(i + 1) % n], ce) if ce > 0 else \
            face_flux(col[i], col[(i + 1) % n], col[(i + 2) % n], ce)
        out[i] = col[i] - (fe - fw) + col[i] * (ce - cw)
    return out
 
 
def run(courant, tend=2.0):
    uf = [[-sin(pi * i * H) ** 2 * sin(2 * pi * (j + .5) * H) for i in range(N + 1)] for j in range(N)]
    vf = [[sin(pi * j * H) ** 2 * sin(2 * pi * (i + .5) * H) for i in range(N)] for j in range(N + 1)]
    nstep = max(1, int(tend * max(abs(x) for r in uf for x in r) / (courant * H)))
    dt = tend / nstep
    f = [[1.0 if ((i + .5) * H - .5) ** 2 + ((j + .5) * H - .75) ** 2 < .15 ** 2 else 0.0
          for i in range(N)] for j in range(N)]
    f0, m0 = [r[:] for r in f], sum(sum(r) for r in f)
    for n in range(nstep):
        w = cos(pi * (n + .5) * dt / tend)          # Rider-Kothe 时间反转
        for ax in ((0, 1) if n % 2 == 0 else (1, 0)):
            if ax == 0:
                f = [line(f[j], [x * w for x in uf[j]], dt / H) for j in range(N)]
            else:
                cols = [[f[j][i] for j in range(N)] for i in range(N)]
                vv = [[vf[j][i] * w for j in range(N + 1)] for i in range(N)]
                cols = [line(cols[i], vv[i], dt / H) for i in range(N)]
                f = [[cols[i][j] for i in range(N)] for j in range(N)]
    lo = min(min(r) for r in f)
    hi = max(max(r) for r in f)
    dm = (sum(sum(r) for r in f) - m0) / m0
    err = sum(abs(f[j][i] - f0[j][i]) for j in range(N) for i in range(N)) / m0
    return nstep, lo, hi, dm, err
 
 
print("   C   steps    min(f)     max(f)-1     dM/M    (dM/M)/C   shape err")
for c in (0.05, 0.1, 0.2, 0.4, 0.8):
    ns, lo, hi, dm, err = run(c)
    print(f"{c:5.2f} {ns:6d}  {lo:9.2e}  {hi - 1.0:9.2e}  {dm:8.2e}  {dm / c:8.4f}   {err:8.3e}")

   C   steps    min(f)     max(f)-1     dM/M    (dM/M)/C   shape err
 0.05   1595   2.23e-29  -6.03e-07  2.94e-03    0.0588   1.988e-01
 0.10    797  -6.51e-07  -6.33e-07  5.87e-03    0.0587   2.119e-01
 0.20    398  -1.62e-06  -5.27e-07  1.16e-02    0.0580   1.850e-01
 0.40    199  -3.87e-06   1.38e-07  2.33e-02    0.0583   1.771e-01
 0.80     99  -3.20e-06   1.43e-06  4.62e-02    0.0577   2.214e-01

把 U 调低，毛细主导的一侧 A 道被 $\Delta t_\sigma$ 拴住，B 道超到前面。把 interface CFL cap 拉到 0.01，B 退回 A 旁边；推到 0.5，下方的体积误差变大。

D / dx120U (flow speed)2.0e-2

interface CFL cap0.05breach factor S8.0

Leave U small — capillary-driven flow — and lane A is bound by dt_sigma while lane B runs away: that is the paper’s selling point. Now drag the interface CFL cap down to 0.01, the CICSAM value: lane B collapses back onto lane A even though surface tension is still implicit. Push the cap up to 0.5 instead and the volume-error readout is what pays for it.

帽在 0.05 时，表里是 1595 步，体积误差是 $2.94\times 10^{-3}$ 。把帽往 0.5 推，步数下降，第五列的 0.058 乘上 $C$ ，账就按这个比例变大。

有界性完好。下冲停在 $10^{-6}$ ，THINC 守住了 0 和 1。体积误差与 $C$ 成正比。第四列除以 $C$ 得到的第五列钉在 0.058 附近。 $C$ 从 0.05 放到 0.80，放大 16 倍，系数依次是 0.0588、0.0587、0.0580、0.0583、0.0577，散开不到 2%。

用积分窗看，这个正比并不奇怪。穿过面的量，积分下限是 $1 - C$ ，窗宽就是 $C$ 。 $C = 0.05$ 只搬走单元最后的 5%， $C = 0.8$ 则一次搬走 80%。步数一起往下掉。 $T = 2$ 、 $N = 40$ 时， $C = 0.05$ 是 1595 步，0.1 是 797，0.2 是 398，0.4 是 199，0.8 是 99。差不多每次折半。 $C$ 加倍，窗加倍，步数减半，单位物理时间上的体积误差照样加倍。

$C = 0.05$ 的 $2.94\times 10^{-3}$ 是 0.294%。 $C = 0.8$ 的 $4.62\times 10^{-2}$ 是 4.62%。两者之比等于 $C$ 之比 16。这是二维面积，折成液滴直径大约是一半，2.3%。曲率是半径的倒数，Laplace 压力跳跃也偏这 2.3%。 $\beta = 2$ 的 tanh 跳跃位置只由单元平均闭合，窗变宽，下冲仍停在 $10^{-6}$ 以下。锐度还在，跟着 $C$ 走的只有体积账本。

CBC 的天花板读的是同一个 $C$ 。 $C = 0.05$ 时天花板是 $20\,\tilde{\phi}_D$ ，施主刚过 0.05 就能把面值抬到 1。 $C = 0.8$ 时天花板是 $1.25\,\tilde{\phi}_D$ 。20 除以 1.25 等于 16，和体积误差放大的倍数相同。压缩余地按 $1/C$ 变小的时候，账上的误差按 $C$ 变大。

最后一列的形状误差与 $C$ 无关，依次是 0.199、0.212、0.185、0.177、0.221。那是 $40^2$ 网格定下的地板。

把网格加密就够了吗

为了弄清系数 0.058 从何而来，我把 $40^2$ 换成 $60^2$ 重跑了 $C = 0.2$ 。系数从 0.0580 降到 0.0384。比值 0.66 几乎就是网格间距之比 $40/60$ 。也就是说，

\frac{\Delta M}{M} \approx 2.3\,C\,\Delta x \propto U\,\Delta t

体积误差在时间上是一阶的。保持 $C$ 不变去加密网格，误差按 $\Delta x$ 成比例下降。但固定 $\Delta t$ 只加密网格， $C$ 会同比例变大，误差原地不动。在界面对流里加大时间步长不是免费的，账单严格按 $\Delta t$ 比例开出。

这张账单若叠上曲率误差，账就错两次。体积偏 0.5%，半径也偏到那个附近，曲率是半径的倒数，就按同样的比例歪。对不上的表面张力会再次弄脏速度场，脏掉的速度接着晃下一步扫掠的 $C$ 。把 $\Delta t$ 放大、窗口仍会关上的隐式表面张力空出来的位子，界面对流的 $C$ 若停在 0.05，就还是关着的。

要把 0.05 提到 0.5，需要改变什么#

论文在结论里先抓住高度函数（height function）的稳健性。欠分辨的界面上高度函数一垮，曲率就整个死掉。接下来才是对流格式。论文的说法是，能允许更大 CFL 的对流，才是把这个方法的性能再抬上去的余地。把对流本身做成隐式，走出 CBC 的 $1/C$ 天花板；改用几何 VOF（PLIC）的非分裂对流，去掉膨胀修正项；或者像用 BVD 挑选 THINC 的重构那样，把面值和输运分开来选。三条路都要交回代数 VOF 一部分便宜的成本。

剩下的不是新上限，而是一个新瓶颈。毛细约束腾出的位子，由界面对流的 CFL 坐了上去。 $\Delta t_\sigma$ 按 $\Delta x^{3/2}$ 缩小，这个瓶颈只按 $\Delta x$ 缩小。网格越细，它相对越喘得过气。是接受 $C = 0.05$ 时 0.294% 体积误差的十倍，还是保住 1595 步，是下一步要做的选择。

[论文评述] 毛细枷锁解开后，CFL 0.05 还在 — VOF 界面对流的真正上限