[論文レビュー] 毛細管の足枷を外したら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$ が小さくなければなりません。論文はこの制約を避けるため、CICSAMではなくTHINC/QQを使いました。

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$ がスイープ途中の実際の値と異なり、その差が体積誤差として残ります。Rider–Kotheの単一渦（時間反転あり、 $T=2$ ）にTHINC移流を載せて測りました。

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）の非分割移流で膨張補正項を消すか、THINCをBVDで選ぶ再構成のように面値を移流と別に選ぶかです。三つとも代数的VOFの安いコストを一部差し出します。

残るのは新しい上限ではなく新しいボトルネックです。毛細管制約が空けた席に界面移流のCFLが座りました。 $\Delta t_\sigma$ は $\Delta x^{3/2}$ で減りますが、このボトルネックは $\Delta x$ でしか減りません。格子を細かくするほど、相対的には息がつきます。 $C = 0.05$ の体積誤差0.294%を十倍で受けるか、1595ステップを保つかは、その次に選びます。

[論文レビュー] 毛細管の足枷を外したらCFL 0.05が残った — VOF界面移流の本当の上限