液滴シミュレーションがΔtに足を取られる理由 — 表面張力を陰的に解く論文

格子を2倍細かくすると、時間刻みは約2.8倍小さくなります。物理が速くなったわけではありません。表面張力が効く界面では、毛細管時間刻み制約（capillary time-step constraint）がほぼ常に一番厳しい制約です。

$\Delta t/\Delta t_\sigma$ のスライダーを0.9より上げると、オレンジ（陽的）だけ振幅が飛びます。シアン（陰的）は同じ定在波を保ちます。

破線が解析的な減衰です。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年に初めて書いた形では、

$\rho_A, \rho_B$ は二流体の密度、$\sigma$ は表面張力係数、$\Delta x$ は格子間隔です。指数 $3/2$ が厄介です。界面を鮮明にしようと $\Delta x$ を半分にすると、時間刻みは $2^{1.5}\approx 2.8$ 倍縮みます。液滴が小さく、毛細管スケールの流れであるほど、ステップ数が爆発します。陽的な表面張力では逃げ道がありません。

なぜ Δx^1.5 なのか — 毛細管波の分散関係#

この制約は実はCFL条件です。界面上の小さな正弦波擾乱は、毛細管波（capillary wave、表面張力が復元力となる波）として振動します。非粘性・深水極限での分散関係は

$\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そのものです。

表面張力を陰的に — 一行の振動子で見る核心

ソルバ全体を持ち出す必要はありません。毛細管波の一モードは減衰調和振動子に帰着します。

$\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、1ステップ後に振幅が掛かる倍率）$|G|$ で境界が鮮明になります。numpy 30行で再現します。トイ問題は衝撃波管ではなく、定常毛細管波の減衰です。

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へ間隔だけが4分の1になるのに、同じ物理時間のステップ数は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$ 自体を動かします。

陽的（オレンジ）は $\Delta t/\Delta t_\sigma \approx 0.9$ まで $|G|=1$ を保ち、その先で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のような大きな密度比でも結合を保った点です。界面移流はCICSAMではなくTHINC/QQで離散化し、大きな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$ をさらに押したあとの窓は 5倍で開き、10倍で閉じた記録 にあります。