静止した液滴が半分に溶けた — 保存形 Allen

保存形 Allen–Cahn のシャープニング項は界面厚みを保つ仕掛けではなく、曲率駆動の消滅をちょうど打ち消す項です。

流速はゼロなのに液滴が半分になった

曲率だけがあって流れのない界面は、自分から動くべきでしょうか。phase-field モデル(界面を有限の厚みをもつ遷移層として扱うモデル)をコードに入れ始めるとき、この問いは抜け落ちやすいところです。静止状態で初期化した2次元液滴を Allen–Cahn 方程式で4000ステップ回すと答えが出ます。半径は20.083から14.286まで縮み、体積率の総和は49.4%消えます。速度場は最後までちょうどゼロのままです。

この記事ではその消滅がどの項から来るのかを項単位で突き止めます。古典 Allen–Cahn と保存形(conservative)Allen–Cahn を並べると、違いは曲率項ひとつです。続いて同じ項が格子ボルツマン(LB)コードのどのモーメントに座るべきかまで追いかけます。読み終えると、コードの W とソース項がそれぞれ何を約束している数なのかが分かれます。

二つの Allen–Cahn の項の対照表#

項目	古典 Allen–Cahn	保存形 Allen–Cahn
形	反応拡散型	フラックス発散型
界面の維持	二重井戸 $\phi(1-\phi)(1-2\phi)$	逆フラックス $\frac{4}{W}\phi(1-\phi)\hat{\mathbf n}$
$\int\phi\,dV$	保存されない	離散レベルで保存
1次元平衡解	tanh	同じ tanh
曲率への応答	$v_n=-M\kappa$	ゼロ

二つの式は次のように書けます。

\frac{\partial \phi}{\partial t} = M\left(\nabla^2\phi - \frac{16}{W^2}\phi(1-\phi)(1-2\phi)\right)

\frac{\partial \phi}{\partial t} + \nabla\cdot(\phi\mathbf u) = \nabla\cdot\left[M\left(\nabla\phi - \frac{4\phi(1-\phi)}{W}\hat{\mathbf n}\right)\right]

$\phi$ は体積率、 $M$ は移動度(界面の緩和の速さを決める係数)、 $W$ は界面厚み、 $\hat{\mathbf n}=\nabla\phi/|\nabla\phi|$ は界面法線、 $\kappa=\nabla\cdot\hat{\mathbf n}$ は曲率です。M と W のスライダーが、左の半円(古典形)と右の半円(保存形)の半径を別々に動かします。

R0 = 20 luW = 4.0 luM = 0.050

t = 0  | classical R = 20.000, m/m0 = 1.0000 | conservative R = 20.000, m/m0 = 1.000000000

左の半円が古典形、右の半円が保存形です。M を上げると左半分だけが速く溶け、点線で描いた $\sqrt{R_0^2-2Mt}$ の上をたどります。W を変えても左側の消滅速度はほとんど動きません。原因は厚みではなく曲率だという意味です。

tanh 界面がただでくれる恒等式#

$W$ が二つの式に別の形で入る理由は、1次元平衡解から出てきます。平面界面の解は次の tanh です。

\phi(x) = \frac{1}{2}\left(1 + \tanh\frac{2x}{W}\right)

これを微分すると恒等式が二つただで手に入ります。

\frac{d\phi}{dx} = \frac{4}{W}\phi(1-\phi), \qquad \frac{d^2\phi}{dx^2} = \frac{16}{W^2}\phi(1-\phi)(1-2\phi)

一つ目は勾配を $\phi$ だけで書けるという意味です。界面座標は要らず、局所値だけでシャープニングフラックスが作れます。二つ目が古典形の二重井戸係数 $16/W^2$ の出どころです。両方をそれぞれの方程式に入れると右辺はゼロになります。1次元では二つのモデルは同じ平衡解を持ちます。

曲率項は残るのか消えるのか

違いは界面が曲がった瞬間に現れます。法線方向の座標 $s$ を取るとラプラシアンは二つに割れます。

\nabla^2\phi = \frac{\partial^2\phi}{\partial s^2} + \kappa\frac{\partial\phi}{\partial s}

ここに上の恒等式を代入します。古典形の右辺では $\partial^2\phi/\partial s^2$ と二重井戸項がちょうど打ち消し合い、曲率の断片だけが残ります。

\left.\frac{\partial\phi}{\partial t}\right|_{\text{classical}} = M\kappa\frac{\partial\phi}{\partial s} = M\kappa\frac{4\phi(1-\phi)}{W}

この一項が $v_n=-M\kappa$ 、つまり曲率駆動運動です。液滴なら $\kappa=-1/R$ なので $dR/dt=-M/R$ 、積分すれば $R^2=R_0^2-2Mt$ になります。半径20、 $M=0.05$ なら寿命は4000時間単位です。

保存形は同じ項を発散形の中でもう一度作ります。逆フラックスの発散を展開すると

\nabla\cdot\left[\frac{4\phi(1-\phi)}{W}\hat{\mathbf n}\right] = \frac{16}{W^2}\phi(1-\phi)(1-2\phi) + \frac{4\phi(1-\phi)}{W}\kappa

となり、二重井戸項と曲率項を一緒に出します。古典形には前半の断片しかありませんでした。後半が加わることで $M\kappa\,\partial\phi/\partial s$ が消えます。3相 Allen–Cahn で存在しない相が界面に漏れていた場所も同じ構造でした。どの項を落とすとどの保存量が壊れるかは、あらかじめ決まっています。

同じ tanh、二つの更新式、4000ステップ#

軸対称1次元格子に二つの更新式をかけました。半径方向だけを解くと曲率が $1/r$ として自動的に入ります。ラプラシアンもシャープニングフラックスも面(face)フラックスで書き、発散部分は両モデルとも離散レベルで保存するようにしています。したがって古典形の質量損失はすべて二重井戸項だけから出ます。

import math
 
NR, DR, R0, W, M, DT, NSTEP = 256, 0.5, 20.0, 4.0, 0.05, 0.5, 4000
rc = [(i + 0.5) * DR for i in range(NR)]          # セル中心
rf = [(i + 1.0) * DR for i in range(NR - 1)]      # 内部面
 
def tanh_disk(w):
    return [0.5 * (1.0 - math.tanh(2.0 * (r - R0) / w)) for r in rc]
 
def disk_mass(f):
    return sum(r * v for r, v in zip(rc, f)) * DR
 
def disk_radius(f):
    return math.sqrt(2.0 * disk_mass(f))
 
def peak_slope(f):
    return max(abs(f[i + 1] - f[i]) / DR for i in range(NR - 1))
 
def well_drive(v, w):
    return (16.0 / w ** 2) * v * (1.0 - v) * (1.0 - 2.0 * v)
 
def flux_div(F):
    out = [0.0] * NR
    for i in range(NR):
        lo = rf[i - 1] * F[i - 1] if i > 0 else 0.0
        hi = rf[i] * F[i] if i < NR - 1 else 0.0
        out[i] = (hi - lo) / (rc[i] * DR)
    return out
 
def grad_faces(f):
    return [(f[i + 1] - f[i]) / DR for i in range(NR - 1)]
 
def sharp_faces(f, w, eps=1e-12):
    F = []
    for i in range(NR - 1):
        g = (f[i + 1] - f[i]) / DR
        pf = 0.5 * (f[i] + f[i + 1])
        F.append(g - (4.0 / w) * pf * (1.0 - pf) * g / (abs(g) + eps))
    return F
 
def step_classical(f, w, dt):
    lap = flux_div(grad_faces(f))
    return [v + dt * M * (lap[i] - well_drive(v, w)) for i, v in enumerate(f)]
 
def step_conservative(f, w, dt):
    div = flux_div(sharp_faces(f, w))
    return [v + dt * M * div[i] for i, v in enumerate(f)]
 
f0 = tanh_disk(W)
m0 = disk_mass(f0)
fk, fc = list(f0), list(f0)
print("    t   R_classic  R_conserv   mass_classic  mass_conserv")
for s in range(NSTEP + 1):
    if s % 1000 == 0:
        print("%5.0f  %9.3f  %9.3f  %13.9f  %13.9f"
              % (s * DT, disk_radius(fk), disk_radius(fc),
                 disk_mass(fk) / m0, disk_mass(fc) / m0))
    fk = step_classical(fk, W, DT)
    fc = step_conservative(fc, W, DT)
print("curvature law  sqrt(R0^2-2*M*t) = %.3f" % math.sqrt(R0 ** 2 - 2 * M * NSTEP * DT))
print("interface width  init %.3f  classic %.3f  conserv %.3f"
      % (1 / peak_slope(f0), 1 / peak_slope(fk), 1 / peak_slope(fc)))
print("mass drift  classic %+.3e  conserv %+.3e"
      % (disk_mass(fk) / m0 - 1, disk_mass(fc) / m0 - 1))

半径は20.083、質量誤差は1e-14#

    t   R_classic  R_conserv   mass_classic  mass_conserv
    0     20.083     20.083    1.000000000    1.000000000
  500     18.803     20.083    0.876607772    1.000000000
 1000     17.429     20.083    0.753192944    1.000000000
 1500     15.936     20.083    0.629689171    1.000000000
 2000     14.286     20.083    0.506056483    1.000000000
curvature law  sqrt(R0^2-2*M*t) = 14.142
interface width  init 4.021  classic 4.032  conserv 4.000
mass drift  classic -4.941e-01  conserv -2.631e-14

古典形の14.286は解析式の14.142から1.0%ずれています。格子幅0.5と面積換算半径の定義から来る差です。傾向は正確に $R^2=R_0^2-2Mt$ に従います。保存形は小数第3位まで20.083から動きません。

質量の数字のほうがはっきりしています。古典形は $-4.941\times10^{-1}$ 、保存形は $-2.631\times10^{-14}$ です。後者は倍精度の累積丸め誤差で、格子を細かくしても減りません。細かくしても減らない側ではなく、最初から存在しない側です。界面厚みは両モデルとも4近辺を保ちます。厚みの維持は二つのモデルの違いではありません。

t=0 では両側の半径が 20.083、質量比が 1.000000000 です。t=500 で動くのは古典形だけで、18.803 と 0.876607772 です。t=1000 は 17.429 と 0.753192944、t=1500 は 15.936 と 0.629689171、t=2000 は 14.286 と 0.506056483 です。保存形の列は 20.083 と 1.000000000 のままです。同じ時刻の解析半径 $\sqrt{R_0^2-2Mt}$ は 14.142 です。14.286 と 14.142 の差 0.144 が上の 1.0% です。質量ドリフト $-4.941\times10^{-1}$ は 1 から 0.506056483 を引いた値です。保存形の $-2.631\times10^{-14}$ とは、桁が 1e-1 と 1e-14 に分かれます。界面幅は初期 4.021、古典形 4.032、保存形 4.000 です。半径 20、 $M=0.05$ の寿命 4000 時間単位は $R^2=R_0^2-2Mt$ を $400=0.1\,t$ と読んだものです。出力の t=2000 はその寿命の半分で、14.142 はその時刻の解析半径です。

LB に移すとこの項はソース項に座る#

保存形を LB で解く標準的な方法は、 $\phi$ 用の分布関数 $h_i$ をもう一つ持たせることです。D2Q9 格子での更新式はこう書きます。

h_i(\mathbf x+\mathbf c_i\Delta t,\ t+\Delta t) = h_i(\mathbf x,t) - \frac{h_i-h_i^{eq}}{\tau_\phi} + F_i

h_i^{eq} = \phi\,w_i\left(1+\frac{\mathbf c_i\cdot\mathbf u}{c_s^2}\right), \qquad M = c_s^2\left(\tau_\phi-\frac{1}{2}\right)\Delta t

$h_i$ は $\phi$ を運ぶ分布関数、 $\tau_\phi$ は無次元緩和時間、 $w_i$ は格子重み、 $c_s^2=1/3$ です。平衡分布には移流項しか入っていません。シャープニング項は平衡分布ではなくソース項 $F_i$ が供給します。

F_i = \Delta t\,w_i\left(1-\frac{1}{2\tau_\phi}\right)\frac{4\phi(1-\phi)}{W}\frac{\mathbf c_i\cdot\hat{\mathbf n}}{c_s^2}

Chapman–Enskog 展開(平衡からの小さなずれを格子間隔比の次数ごとに整理する方法)が要求する条件はモーメント2行です。0次モーメントは $\sum_i F_i = 0$ でなければなりません。 $\phi=\sum_i h_i$ なので、0次モーメントがゼロでなければソース項が体積率を作ったり消したりします。1次モーメントは逆フラックスを正確に出さなければなりません。 $\sum_i w_i(\mathbf c_i\cdot\hat{\mathbf n})\mathbf c_i = c_s^2\hat{\mathbf n}$ なので、上の $F_i$ では $c_s^2$ が約分されて目標値になります。D2Q9 の平衡が多項式で閉じる場所と同じ種類の条件で、ここでは守るべき行が2本で終わります。

0次モーメントを触ると質量が漏れ出す#

よくある事故は $F_i$ を平衡分布のように書くことです。つまり括弧の中に1を一緒に入れて $w_i\,\mathrm{amp}\,(1+\mathbf c_i\cdot\hat{\mathbf n}/c_s^2)$ とする場合です。1次モーメントはそのまま合います。だから界面厚みも正常に保たれ、コードは静かに回ります。

CX = [0, 1, -1, 0, 0, 1, -1, -1, 1]
CY = [0, 0, 0, 1, -1, 1, 1, -1, -1]
WT = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]
CS2 = 1.0 / 3.0
W, TAU, PHI, NX, NY = 4.0, 0.8, 0.5, 0.6, 0.8
 
def sharp_source(phi, nx, ny, w, tau, spill=False):
    amp = (1.0 - 0.5 / tau) * 4.0 * phi * (1.0 - phi) / w
    out = []
    for k in range(9):
        proj = (CX[k] * nx + CY[k] * ny) / CS2        # 1次モーメント用の射影
        out.append(WT[k] * amp * (proj + 1.0 if spill else proj))
    return out
 
amp = (1.0 - 0.5 / TAU) * 4.0 * PHI * (1.0 - PHI) / W
print("target            sum c_i F_i = [%+.6f %+.6f]" % (amp * NX, amp * NY))
for tag, spill in (("moment-safe", False), ("0th-moment spill", True)):
    F = sharp_source(PHI, NX, NY, W, TAU, spill)
    sx = sum(CX[k] * F[k] for k in range(9))
    sy = sum(CY[k] * F[k] for k in range(9))
    print("%-17s sum F_i = %+.3e   sum c_i F_i = [%+.6f %+.6f]" % (tag, sum(F), sx, sy))

target            sum c_i F_i = [+0.056250 +0.075000]
moment-safe       sum F_i = +8.674e-19   sum c_i F_i = [+0.056250 +0.075000]
0th-moment spill  sum F_i = +9.375e-02   sum c_i F_i = [+0.056250 +0.075000]

W_init = 3.0 luW_model = 6.0 lu

t = 0  |  W_eff = 3.000 lu  |  max |net flux| = 0.00e+0  | mass / mass0 = 1.000000

source を入れると0次モーメントが壊れます。W_init と W_model を違う値にすると赤い net flux が立ち上がり、点線の tanh まで下がります。W_eff は W_model に張り付いたまま、質量曲線だけが1から離れます。厚みの数字だけではこの穴は見えません。

二つの変形の1次モーメントは小数第6位まで同じです。0次モーメントだけがゼロから $9.375\times10^{-2}$ に上がります。この値は $(1-1/2\tau_\phi)\cdot4\phi(1-\phi)/W$ そのものです。界面の格子点ごとに毎ステップその分の体積率を差し込んでいるという意味です。

ソース出力の1次モーメントは両モードとも +0.056250 と +0.075000 です。 $\phi=0.5$ 、 $\tau_\phi=0.8$ 、 $W=4$ なら $(1-1/(2\tau_\phi))\cdot 4\phi(1-\phi)/W = 0.09375$ です。 $n_x=0.6$ を掛けると 0.056250、 $n_y=0.8$ を掛けると 0.075000 です。0次だけが +8.674e-19 と +9.375e-02 に分かれます。9.375e-02 は 0.09375 を指数で書き直した値です。幅が 4.000 で合っていても、この和が 0 でなければ質量は漏れます。

二つのモデルを選ぶときの対照表

状況	古典形	保存形
平面界面の緩和	適する	適する
静止液滴の長時間保持	適さない	適する
結晶粒成長・曲率駆動運動の模擬	適する	適さない
体積率の総和の報告	補正が必要	そのまま
LB 実装コスト	ソース項なし	$\hat{\mathbf n}$ の計算が必要

表の3行目が重要です。保存形が常に正しいわけではありません。曲率駆動運動を模擬したい状況なら、古典形が正しいモデルです。流体界面を移流させたい状況でのみ、その項が誤差になります。どちらにしても $\hat{\mathbf n}$ を計算する必要があるので、保存形では $|\nabla\phi|$ がゼロに近いバルク領域での除算の保護が追加作業として残ります。

界面捕獲コードに W を書き込むとき#

コードに W = 4 と書いた瞬間に決まるのは厚みひとつではありません。二重井戸係数 $16/W^2$ と逆フラックス係数 $4/W$ が同じ数から出ています。どちらか一方だけを変えると平衡解が tanh から外れます。二つの係数を同じ $W$ から導いているかが最初の点検項目です。

次の点検は保存量です。体積率の総和を毎ステップ出力してみれば、1e-14と1e-1は一目で分かれます。厚みは両モデルで正常に見えるので、厚みだけを見ても何も濾し取れません。移動度を10倍にしても界面厚みがそのままだった保存形では、厚みと質量は別の約束でした。LB 実装なら点検は1行で終わります。 $\sum_i F_i$ を出力してゼロかどうか確認することです。

静止した液滴が半分に溶けた — 保存形 Allen–Cahn が消す曲率項