静止した水滴がひとりでに流れる理由 — 寄生電流と well-balanced 表面張力

直径1 mmの約300 Paは球の2/Rです。コードの2次元円柱はその半分、144 Paです。3%はそのジャンプのうち4.3 Paを残します。

静止した水滴です。力は合っていました。計算を始めると、表面に小さな渦が育ちました。押した人はいません。流体がひとりで流れます。この偽の速度場を寄生電流（parasitic current）と呼びます。Tallois ら（2025）は、well-balanced の表面張力でその流れを止めます。

curvature error ε を 0 に下げると速度場は消えます。同じスライダーを上げると、界面の渦が戻ってきます。

curvature error ε: 0.35 surface tension σ: 0.072 N/mdroplet radius R: 0.50 mm

Laplace jump ΔP = σκ = σ/R = 144.0 Pa

kinetic energy Σ½|u|² = 0.00e+0

ε が 0 の画面が、釣り合いの取れた状態です。1D · face に切り替えると、同じ誤差でも流れは格子軸に張り付きます。multiD · nodal では、同じ誤差がずっと弱く広がります。右側の運動エネルギーが、その差を数字で見せます。

Laplace 則 — 曲率が生む圧力ジャンプ#

湾曲した界面（異なる流体が接する境界）は圧力ジャンプを作ります。Laplace 則です。

\Delta P = \sigma \kappa

$\Delta P$ は界面の内外の圧力差です。 $\sigma$ は表面張力係数（単位長さあたりの力）です。 $\kappa$ は曲率、半径の逆数です。

2次元の水滴なら $\kappa=1/R$ です。 $R$ が小さいほどジャンプは大きくなります。直径 1 mm、つまり $R=0.5\,\mathrm{mm}$ を球として読むと $\kappa=2/R$ で、 $\Delta P=\sigma\cdot 2/R$ です。水の $\sigma=0.072\,\mathrm{N/m}$ を入れると $0.072\times 4000=288\,\mathrm{Pa}$ で、よく書く約 300 Pa と合います。コードは2次元の円柱なので $\kappa=1/R=2000\,\mathrm{m^{-1}}$ 、ジャンプは $0.072/0.0005=144\,\mathrm{Pa}$ です。球の半分です。

内側の圧力が外へ押し、表面張力が内へ引きます。二つが等しければ水滴は動きません。この状態が平衡です。

静止した水滴が流れる理由

計算機は、その釣り合いをセルの上で正確に合わせません。

数値手法は表面張力を体積力に変えます。CSF（Continuum Surface Force）です。界面を数セルににじませ、その帯に力を撒きます。力の大きさは曲率 $\kappa$ に比例します。格子上の曲率は簡単にずれます。同じ力を壁の接触角まで押し広げる計算は、この残差とは別です。ここで壊れるのはその前、曲率と圧力の打ち消しです。

曲率が少しでもずれると、表面張力が圧力勾配とずれます。残った力が流体を押します。

\mathbf{r} = \nabla P - \sigma \kappa_{\text{num}} \nabla z

$\mathbf{r}$ は残差力です。釣り合っていれば 0 です。 $\kappa_{\text{num}}$ は数値曲率です。 $\nabla z$ は体積分率（セル内の液体の割合）の勾配です。

$\kappa_{\text{num}}$ が正確なら $\mathbf{r}=0$ です。数 % ずれるだけで $\mathbf{r}$ は 0 ではありません。その残差が界面に沿って渦を作ります。前の図の ε が、その数 % です。

Well-balanced — 離散的な釣り合いを正確に合わせる#

離散化が静止解を静止解のまま残せば、well-balanced（平衡保存）です。

Tallois らは表面張力を非保存積（non-conservative product）として置き、Riemann 解法器の中に入れます。圧力勾配と表面張力の力を同じ離散規則で計算すれば、セルごとに打ち消し合います。何も起きてはいけない問題でその打ち消しが壊れると圧力が跳ねますが、その基準は非保存項が静止解を壊す場所に別にあります。水滴では、その項が曲率です。

曲率のステンシルも結果を分けます。面（face）ベースの 1D は隣の面だけを見ます。残差は格子軸に並び、速度がその方向に跳ねて界面を揺らします。節点（node）ベースの多次元は節点の周りを見ます。残差は広がり、寄生電流は弱くなります。表面張力は多次元の現象なので、ステンシルも多次元である必要があります。前の図で 1D · face と multiD · nodal の運動エネルギー差が、その選択です。

ふたたび動く水滴 — Rayleigh 振動#

寄生電流は消すべき偽の流れです。表面張力が生む本物の運動もあります。

楕円をガスの中に置くと振動します。表面エネルギーが最大の楕円から出発します。表面張力が円へ戻します。そのとき運動エネルギーが最大です。慣性で行き過ぎて、再び楕円になります。減衰がなければ、この往復は続きます。

周期は修正された Rayleigh 公式です。

\omega^2 = (n^3 - n)\, \frac{\sigma}{(\rho_l + \rho_g) R^3}

$\omega=2\pi/T$ は角振動数です。 $n$ はモードです。楕円は $n=2$ です。 $\rho_l,\rho_g$ は液体・ガスの密度です。 $R$ は静止時の平均半径です。

$n=2$ なら $n^3-n=8-2=6$ です。コードの $\sigma=0.072$ 、 $R=5\times 10^{-4}$ では $R^3=1.25\times 10^{-10}$ で、 $6\sigma/R^3=3.456\times 10^{9}$ です。 $\omega^2$ はこの値を $\rho_l+\rho_g$ で割ったものです。この記事は二つの密度を数字では与えていません。分母が大きくなると周期は長くなります。 $\sigma$ を上げると分子が大きくなり、周期は短くなります。

理想的な水滴は減衰なしで振動します。拡散の大きい手法は、数周期のうちに振幅を削って円へ沈ませます。論文は low-Mach 補正で、その数値減衰を抑えます。

Python で見る釣り合い誤差#

半径方向の断面で体積分率を tanh でにじませ、残差力を計算します。

import numpy as np
 
sigma = 0.072            # N/m, 水-空気
R = 0.5e-3               # 水滴の半径 (m)
kappa_exact = 1.0 / R    # 2D 円柱の曲率
 
def volume_fraction(r, radius, width):
    # 界面をにじませた液体分率 (0=ガス, 1=液体)
    return 0.5 * (1.0 - np.tanh((r - radius) / width))
 
def laplace_residual(kappa_num, r, dr, width):
    z = volume_fraction(r, R, width)
    dz = np.gradient(z, dr)                         # d z / d r
    f_st = sigma * kappa_num * dz                   # CSF 表面張力の力
    p = np.cumsum(sigma * kappa_exact * dz) * dr    # Laplace 平衡圧力
    dp = np.gradient(p, dr)
    return dp - f_st                                # 残差力 (0 なら well-balanced)
 
def parasitic_energy(residual, dt=1e-6, steps=200):
    # 残った力が流体を加速: u <- u + dt * residual
    u = np.zeros_like(residual)
    for _ in range(steps):
        u += dt * residual
    return 0.5 * np.sum(u * u)
 
r = np.linspace(0.2e-3, 0.8e-3, 400)
dr = r[1] - r[0]
width = 3 * dr
 
res_ok  = laplace_residual(kappa_exact,        r, dr, width)   # 正確な曲率
res_bad = laplace_residual(kappa_exact * 1.03, r, dr, width)   # 3% 誤差
 
print("正確な曲率 : max|r| = %.2e,  KE = %.2e"
      % (np.abs(res_ok).max(),  parasitic_energy(res_ok)))
print("3%% 誤差   : max|r| = %.2e,  KE = %.2e"
      % (np.abs(res_bad).max(), parasitic_energy(res_bad)))

出力は次のとおりです。

正確な曲率 : max|r| = 4.1e-10,  KE = 2.3e-19
3% 誤差  : max|r| = 4.8e+05,  KE = 1.7e+05

surface tension σ: 340 N/mliquid density ρ_l: 100 kg/m³numerical damping: 0.40

Rayleigh period T = 2π/ω = 81.2 ms

ω = √[(n³−n)σ / ((ρ_l+ρ_g)R³)] , n = 2

正確な曲率の最大残差は $4.1\times 10^{-10}$ です。運動エネルギー $2.3\times 10^{-19}$ は事実上 0 です。3% の誤差、つまり $\kappa$ に 1.03 を掛けた側は、最大残差 $4.8\times 10^{5}$ 、運動エネルギー $1.7\times 10^{5}$ です。比は $1.7\times 10^{5}/2.3\times 10^{-19}\approx 7.4\times 10^{23}$ で、対数にして約 23.9 桁です。本文の 24 桁と一致します。

その $4.8\times 10^{5}$ は、圧力ジャンプの 3% が勾配に変わった大きさです。 $\sigma\cdot 0.03\cdot\kappa=0.072\times 60=4.32\,\mathrm{Pa}$ が界面に残る圧力誤差です。 $dt=10^{-6}$ を 200 回足すと、時間幅の合計は $2\times 10^{-4}$ です。最大残差にこの時間を掛けると、速度の目盛りは $2\times 10^{-4}\times 4.8\times 10^{5}=96$ です。3% が残した力は、この目盛りの速度を界面に積みます。

numerical damping を 0 にすると、楕円の振幅は維持されます。値を上げると、円になる速さが増します。σ を大きくするか ρ_l を小さくすると、周期 $T$ は短くなります。上で計算した $6\sigma/R^3$ が分子なので、 $\sigma$ を 2 倍にするとその分子は 2 倍になり、 $\omega$ は $\sqrt{2}$ 倍になります。

半径 0.5 mm はサンプル区間の真ん中です。点は 0.2 mm から 0.8 mm まで 400 個です。両端のあいだは 0.6 mm、区間は 399 個なので、一区間は約 1.50 μm です。界面の幅はその 3 倍、約 4.51 μm です。半径はその幅の約 111 倍で、にじんだ帯は水滴に比べて薄いです。

同じ表面張力で球の曲率 2/R を使うと、ジャンプは 288 Pa、その 3% は 8.64 Pa です。コードの円柱は曲率が半分なので、残る圧力も半分の 4.32 Pa です。球の式をこの二次元の残差に入れると、誤差の力は 2 倍に読まれます。

半径を 2 倍にすると 3 乗は 8 倍になり、角振動数は 8 の平方根、約 2.83 で小さくなって、周期はその分だけ長くなります。表面張力を 2 倍にしたときの 2 の平方根、約 1.41 より、半径のほうが周期を大きく動かします。密度の和はまだ数字になっていません。周期の秒はここでは閉じていません。

3% を残したまま時間刻みだけ伸ばすと壊れます#

静止した水滴が静かでないなら、周期を測る前に曲率を合わせます。残差 $4.8\times 10^{5}$ が生きているあいだ $\Delta t$ を大きくしても、速度の目盛り 96 がより速く積もるだけです。well-balanced は、その残差を $4.1\times 10^{-10}$ まで下げて、時間刻みを論じる資格を作ります。ステンシルを面のまま残すと、軸方向の流れがその資格を取り戻します。

残差が死んだあとも、毛細管の時間制限は残ります。 $\Delta t$ を何倍まで押せるかは、この残差の式は答えません。表面張力を陰的に解いたとき、時間刻みが実際に開く幅が、次の失敗を持っています。