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

曲率誤差が生む寄生電流と、それを消し去る well-balanced 手法

静止した水滴は動かないはずです。力は完璧に釣り合っています。ところがシミュレーションを起動すると、水滴の表面に小さな渦がじわじわと育ちはじめます。誰も押していないのに流体が流れるのです。この幽霊のような流れを寄生電流(parasitic current・数値誤差から生じる偽の速度場)と呼びます。Tallois ら(2025)は、この問題を well-balanced 表面張力手法で抑え込みます。今日は、なぜ水滴がひとりでに流れるのか、そしてどうすれば止められるのかを追っていきます。

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

湾曲した界面(異なる流体どうしが接する境界)は圧力ジャンプを生みます。これが Laplace 則です。

\Delta P = \sigma \kappa

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

2次元の水滴なら $\kappa = 1/R$ です。半径 $R$ が小さいほど圧力ジャンプは大きくなります。直径 1mm の水滴は内側が約 300Pa 高くなります。この圧力差が、水滴を球形にまとめる力です。

ここで重要なのは、これが平衡だという点です。内側の高い圧力が外へ押し、表面張力が内へ引きます。二つの力が正確に等しくなります。だから水滴はじっとしています。

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

問題は、コンピュータがこの釣り合いを正確に合わせられないところにあります。

数値手法は、表面張力を体積力に変換して組み込みます。これを CSF(Continuum Surface Force・界面を数セルにわたってにじませ、その帯に力を撒く方式)と呼びます。力の大きさは曲率 $\kappa$ に比例します。ところが格子の上で曲率を正確に計算するのが難しいのです。

曲率をほんの少し間違って計算すると、表面張力の力が圧力勾配とずれます。残った力が流体を押します。

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

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

曲率が正確なら $\mathbf{r} = 0$ となり、水滴は静かなままです。しかし $\kappa_{\text{num}}$ に数 % でも誤差が生じると $\mathbf{r} \neq 0$ になります。この残差が界面に沿って渦を作ります。それが寄生電流です。

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

解決の鍵は well-balanced(平衡保存)という性質です。離散化された方程式が、静止解を正確に保存するように設計することです。

Tallois らは表面張力を非保存積(non-conservative product)として扱います。そしてこの項を Riemann 解法器(界面で波を解く数値ツール)の中に直接組み込みます。圧力勾配と表面張力の力を同じ離散規則で計算すれば、両者がセル単位で正確に打ち消し合います。

曲率の計算方式も重要です。論文では、曲率を面(face)ではなく節点(node)ベースのステンシルで計算します。なぜでしょうか。

1D(面ベース): 隣接する面だけを見ます。残差が格子軸に整列します。速度場が格子方向に跳ね、界面を不安定にします。
多次元(節点ベース): 節点の周囲をすべて見ます。残差が滑らかに広がります。寄生電流がはるかに弱くなります。

表面張力は本質的に多次元の現象です。だから真の多次元ステンシルが有利になります。

下のシミュレーションで実際に操作してみましょう。curvature error ε を 0 に下げると速度場が消えます。それが well-balanced の状態です。

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

ε を上げると界面に渦が育ちます。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$ は静止時の平均半径です。

ここで well-balanced が再び重要になります。理想的には、水滴は減衰なしで振動するはずです。しかし拡散の大きい手法は振幅を素早く削ります。水滴が数周期でつぶれて円に落ち着いてしまいます。論文では low-Mach 補正でこの数値減衰を抑え込み、複数の周期にわたって振動を維持します。

下で実際に操作してみましょう。

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

numerical damping を 0 にすると振幅が維持されます。これが拡散のない理想的な手法です。値を上げると水滴が素早く円へと沈みます。σ を大きくするか ρ_l を小さくすると周期 $T$ が短くなります。公式のとおりです。

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

正確な曲率では残差が機械誤差の水準です。運動エネルギーも事実上 0 です。曲率をたった 3% 間違えるだけで残差が爆発します。運動エネルギーが 24 桁も跳ね上がります。これが画面で見た寄生電流の正体です。

覚えておくべきこと

Laplace 則 $\Delta P = \sigma \kappa$ は静止水滴の平衡です。内側の圧力と表面張力が正確に打ち消し合います。
寄生電流は、曲率誤差がその打ち消しを壊すことで生じる偽の速度場です。曲率を 3% 間違えるだけで運動エネルギーが爆発します。
Well-balanced 手法は、圧力と表面張力を同じ離散規則で計算して平衡を正確に保存します。節点ベースの多次元ステンシルが格子整列の不安定を防ぎます。