Skip to content
cfd-lab:~/ja/posts/2026-09-18-conservative-…● online
$_
Younglin's CFD Blog
毎日の記録 · 計算流体力学

この記事は著者が韓国語で書いた原文の翻訳です。 原文(한국어)を読む

NOTE #158DAY FRI CFD기법DATE 2026.09.18UPDATED 2026.10.19READ 7 min read#Allen-Cahn#Phase-Field#LBM#Multiphase#Interface-Capturing

移動度を10倍にしても界面の厚みは動かなかった — 保存型Allen–Cahnの反拡散項

界面の厚みを決めるのは移動度Mではなく W です。M が決めるのは厚みではなく、崩れた界面が戻る速さだけです。

拡散を10倍にしたのに界面は厚くならなかった#

相場(phase field、界面を0と1のあいだを滑らかに渡るスカラーで表す方法)のコードで、移動度 MM を0.02から0.20へ上げました。拡散係数を10倍にしたことになります。界面は10倍速くなまるはずです。測った厚みは4.44から4.24へ、わずかに減りました。

計算が壊れたわけではありません。方程式がそういう形をしています。格子ボルツマン法(LBM)で界面を保持する保存型Allen–Cahnは厚みを WW に縛り、反拡散は平衡分布の一行に入ります。∇ϕ\nabla\phi を格子でどう求めるかも、同じ数値の側に残ります。下の数値は、併記したコードを走らせて得たものです。

smear interface を押すと界面はいったん厚くなり、反拡散がオンの側だけ緑の破線の厚みへ戻ります。

W = 8.00mass drift = 0.0000%

Press smear interface and watch the orange width trace. With anti-diffusion ON the profile is pulled back onto the dashed tanh of width W whatever the mobility M is; with it OFF the same M only decides how fast the interface keeps swelling. The mass readout stays at zero either way, because both terms are written as a divergence of a flux.

smear interface を押してわざと界面を崩し、下段のオレンジ色の厚み履歴を見ます。反拡散がオンならどの MM でも緑の破線の高さへ戻り、オフなら MM は厚みが増える速さを決めるだけです。

界面を方程式ひとつで押さえるという発想

VOFやlevel setは移流を解いたあと、界面を再構成したり距離関数を再初期化したりします。別工程です。同じ保存型Allen–Cahnでも曲率項を落とすと、静止した液滴が溶けます。その場所は曲率項が消える保存型Allen–Cahnにあります。

保存型Allen–Cahnは、その後処理を同じ方程式の中に畳み込みます。

∂ϕ∂t+∇⋅(ϕ u)=∇⋅[ M(∇ϕ−4ϕ(1−ϕ)W n)]\frac{\partial \phi}{\partial t} + \nabla\cdot(\phi\,\mathbf{u}) = \nabla\cdot\left[\,M\left(\nabla\phi - \frac{4\phi(1-\phi)}{W}\,\mathbf{n}\right)\right]

ϕ\phi は相場変数(一方の相で0、他方で1)、MM は移動度、WW は目標とする界面厚み(格子セル数)、n=∇ϕ/∣∇ϕ∣\mathbf{n}=\nabla\phi/|\nabla\phi| は界面法線です。括弧内の第1項は界面を広げる通常の拡散、第2項が逆向きに押し戻す反拡散です。

右辺全体が発散形になっている点が効きます。周期境界では ∫ϕ dV\int\phi\,dV が機械精度まで保存されます。質量を保存しないことで知られるAllen–Cahnの一形式が「保存型」と呼ばれる理由がここにあります。

反拡散項がちょうどゼロになる場所

界面法線方向の座標を ξ\xi とすると、平衡プロファイルはこうなります。

ϕ(ξ)=12[1+tanh⁡2ξW]\phi(\xi) = \frac{1}{2}\left[1+\tanh\frac{2\xi}{W}\right]

微分すれば勾配が出ます。

dϕdξ=1W sech22ξW\frac{d\phi}{d\xi} = \frac{1}{W}\,\mathrm{sech}^2\frac{2\xi}{W}

一方、同じプロファイル上で 4ϕ(1−ϕ)=1−tanh⁡2(2ξ/W)=sech2(2ξ/W)4\phi(1-\phi) = 1-\tanh^2(2\xi/W) = \mathrm{sech}^2(2\xi/W) です。2つの式は同じ関数です。つまり平衡プロファイル上では

∇ϕ=4ϕ(1−ϕ)W n\nabla\phi = \frac{4\phi(1-\phi)}{W}\,\mathbf{n}

が恒等的に成り立ちます。括弧の中がゼロです。そして MM はそのゼロに掛かります。

ここで分かれることは二つです。平衡厚みは WW だけで閉じ、MM はその席に入りません。MM が決めるのは平衡から押し出されたあとの回復時間だけで、その時間はおおむね W2/MW^2/M で動きます。

格子の上ではそのゼロは自動では保たれない

連続系でぴたりと打ち消し合う2項も、離散化すると別々に動きます。拡散と反拡散を格子点上で中心差分により個別に評価すると、どちらのステンシルも i−1i-1 と i+1i+1 だけを見て ii を飛ばします。奇数点と偶数点が互いを知らない状態になり、その格子スケールのモードを反拡散が増幅します。実際その形で組んだところ、M=0.20M=0.20 で厚みが目標の4.0を通り越して2.88まで下がりました。

対処は、フラックスを格子点ではなく面で作ることです。面 i+1/2i+1/2 のフラックスを ϕi\phi_i と ϕi+1\phi_{i+1} の2値だけで組めば、隣り合う点が互いを見ます。分離モードの入る隙がありません。

同じ物差しで厚みを8回測った#

スラブを厚み8で初期化し、目標 W=4W=4 に向けて8,000ステップ回しました。反拡散のオン・オフを MM の4値それぞれで比べます。

import numpy as np
 
L, dx, dt = 384, 1.0, 0.2
W, u, NSTEP, R = 4.0, 0.05, 8000, 96.0
x = np.arange(L) * dx
 
def slab(width):
    d = (x - L / 2 + L / 2) % L - L / 2
    return 0.5 * (1.0 + np.tanh(2.0 * (R - np.abs(d)) / width))
 
def measured_width(phi):
    g = (np.roll(phi, -1) - np.roll(phi, 1)) / (2.0 * dx)
    return 1.0 / np.max(np.abs(g))
 
def face_flux(phi, M, sharp):
    # 面 i+1/2 の上で拡散と反拡散をまとめて作る
    g = (np.roll(phi, -1) - phi) / dx
    pf = 0.5 * (phi + np.roll(phi, -1))
    return M * (g - sharp * (4.0 * pf * (1.0 - pf) / W) * np.sign(g))
 
def step(phi, M, sharp):
    F = face_flux(phi, M, sharp)
    adv = u * (np.roll(phi, -1) - np.roll(phi, 1)) / (2.0 * dx)
    return phi + dt * (-adv + (F - np.roll(F, 1)) / dx)
 
def run(M, sharp, w_init=8.0):
    phi = slab(w_init)
    m0 = phi.sum() * dx
    for _ in range(NSTEP):
        phi = step(phi, M, sharp)
    return measured_width(phi), phi.sum() * dx - m0
 
print("ruler: an exact tanh of W=4.0 on this grid reads %.4f" % measured_width(slab(W)))
print("start W=8.0, %d steps at u=%.2f, target W=%.1f\n" % (NSTEP, u, W))
print("%-8s %-12s %-12s %-12s" % ("M", "anti-diff", "W_read", "mass drift"))
for M in (0.02, 0.05, 0.10, 0.20):
    for sharp in (1.0, 0.0):
        w, dm = run(M, sharp)
        print("%-8.2f %-12s %-12.4f %-+12.2e" % (M, "ON" if sharp else "OFF", w, dm))
w_lo, _ = run(0.05, 1.0, w_init=2.0)
print("\nfrom the thin side: W_init=2.0, M=0.05, ON -> %.4f" % w_lo)
ruler: an exact tanh of W=4.0 on this grid reads 4.3279
start W=8.0, 8000 steps at u=0.05, target W=4.0
 
M        anti-diff    W_read       mass drift  
0.02     ON           4.4403       -1.10e-11   
0.02     OFF          21.8974      -1.34e-11   
0.05     ON           4.2854       -2.95e-11   
0.05     OFF          32.9005      -6.76e-12   
0.10     ON           4.2442       -3.22e-11   
0.10     OFF          45.6971      -9.66e-13   
0.20     ON           4.2357       -3.60e-11   
0.20     OFF          64.0262      -2.27e-13   
 
from the thin side: W_init=2.0, M=0.05, ON -> 4.2854

最初の行が物差しです。厚みを 1/max⁡∣∇ϕ∣1/\max|\nabla\phi| で測っていますが、厳密な幅4のtanhを格子に載せて同じ方法で測ると4.0ではなく4.3279と出ます。物差し自身の離散化誤差です。

その物差しに当てれば、反拡散オンの4値4.24〜4.44はすべて ±3%\pm3\% の内側です。移動度が10倍動くあいだ、厚みはほぼ動いていません。反拡散を切ると同じ範囲で21.9から64.0まで開きます。質量のずれはどちらも 10−1110^{-11} 程度です。

W=4W=4、M=0.02M=0.02 なら W2/M=16/0.02=800W^2/M = 16/0.02 = 800 です。M=0.20M=0.20 なら80です。コードの8000ステップに dt=0.2dt=0.2 を掛けると積分時間は1600です。小さい移動度でも回復時間の2回分、大きい移動度では20回分を覆います。だからONの4行が4.4403、4.2854、4.2442、4.2357に集まります。物差し自身は幅4のtanhを4.3279と読みます。1/4.32791/4.3279 は約0.231で、連続の最大勾配 1/W=0.251/W = 0.25 より小さいです。4.24〜4.44が ±3%\pm3\% の内側という話は、4.0ではなくこの4.3279を基準にしています。

最後の行は、これが単なる「変化なし」ではないことを示します。厚み2から始めても4.2854まで上がり、上から下りてきた値と小数第4位まで一致します。平衡プロファイルは両側から引き寄せるアトラクターです。

LBMに移すと、この項はソース項ではない#

ϕ\phi に分布関数 hih_i を別に立てます。ここでよくある回り道が、反拡散を衝突後に足すソース項として組むことです。その必要はありません。

hi(x+ciδt, t+δt)=hi(x,t)−1τϕ[hi(x,t)−hieq]h_i(\mathbf{x}+\mathbf{c}_i\delta t,\ t+\delta t) = h_i(\mathbf{x},t) - \frac{1}{\tau_\phi}\left[h_i(\mathbf{x},t) - h_i^{\mathrm{eq}}\right]

平衡分布をこう取ります。

hieq=wi ϕ+wics2 ci⋅[ϕ u+M 4ϕ(1−ϕ)W n]h_i^{\mathrm{eq}} = w_i\,\phi + \frac{w_i}{c_s^2}\,\mathbf{c}_i\cdot\left[\phi\,\mathbf{u} + M\,\frac{4\phi(1-\phi)}{W}\,\mathbf{n}\right]

wiw_i は格子重み、cs2=1/3c_s^2=1/3 は格子音速の2乗、τϕ\tau_\phi は無次元緩和時間です。角括弧の中が平衡分布の1次モーメント、つまり ϕ\phi を運ぶ速度です。物理速度 u\mathbf{u} の隣に、反拡散が速度の顔をして並んでいます。

Chapman–Enskog展開をかけると、拡散係数は緩和時間から出てきます。

M=cs2(τϕ−12)δtM = c_s^2\left(\tau_\phi - \frac{1}{2}\right)\delta t

2つを合わせると、先の保存型Allen–Cahn方程式がそのまま復元されます。緩和が +M∇2ϕ+M\nabla^2\phi を作り、平衡分布に載せた1次モーメントが −∇⋅(M 4ϕ(1−ϕ)n/W)-\nabla\cdot(M\,4\phi(1-\phi)\mathbf{n}/W) を作ります。表の M=0.10M=0.10 を cs2=1/3c_s^2=1/3 で割ると (τϕ−1/2)δt=0.30(\tau_\phi - 1/2)\delta t = 0.30 です。M=0.02M=0.02 なら0.06、M=0.20M=0.20 なら0.60です。移動度を10倍にするとはこの積を伸ばすことであって、WW に触れることではありません。

実装上の利点ははっきりしています。streamingとcollisionの構造に手を入れません。界面捕捉の仕事が平衡分布の一行に収まります。運動量側の表面張力は別のつまみです。厚みのつまみと混ぜてはいけません。

代償は n\mathbf{n} を全格子点で求めること#

ただではありません。平衡分布に n=∇ϕ/∣∇ϕ∣\mathbf{n}=\nabla\phi/|\nabla\phi| が入っている以上、∇ϕ\nabla\phi をどこでも評価する必要があり、ステンシルの選択が結果を分けます。

素直な選択は軸方向の4近傍だけを使う中心差分です。格子ボルツマンにはもう一つの手があります。格子速度集合全体に重みを載せるやり方です。

∇ϕ≃1cs2∑iwi ci  ϕ(x+ciδt)\nabla\phi \simeq \frac{1}{c_s^2}\sum_i w_i\,\mathbf{c}_i\;\phi(\mathbf{x}+\mathbf{c}_i\delta t)

D2Q9では重みは軸方向が 1/91/9、対角方向が 1/361/36 です。対角近傍が入ることで4次の異方性誤差が打ち消されます。同じ液滴の上での二つのステンシルの誤差は、下の出力に数値で残っています。

平均誤差が大きいほうを選ぶ理由

半径24、厚み WW のtanh液滴を96×96格子に載せ、界面リング上の ∣∇ϕ∣|\nabla\phi| の相対誤差を角度24区間に分けて測りました。

import numpy as np
 
N, R = 96, 24.0
cx = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
cy = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1])
w = np.array([4/9] + [1/9]*4 + [1/36]*4)
 
def fields(W):
    i = np.arange(N)
    X, Y = np.meshgrid(i - N/2 + 0.5, i - N/2 + 0.5, indexing='ij')
    r = np.sqrt(X**2 + Y**2)
    phi = 0.5 * (1.0 + np.tanh(2.0 * (R - r) / W))
    exact = (1.0 / W) / np.cosh(2.0 * (R - r) / W)**2
    return X, Y, r, phi, exact
 
def grad_central(phi):
    gx = (np.roll(phi, -1, 0) - np.roll(phi, 1, 0)) / 2.0
    gy = (np.roll(phi, -1, 1) - np.roll(phi, 1, 1)) / 2.0
    return np.sqrt(gx**2 + gy**2)
 
def grad_lattice(phi):
    # 対角近傍まで格子重みでまとめる
    gx = np.zeros_like(phi); gy = np.zeros_like(phi)
    for k in range(9):
        s = np.roll(np.roll(phi, -cx[k], 0), -cy[k], 1)
        gx += w[k] * cx[k] * s
        gy += w[k] * cy[k] * s
    return np.sqrt((3.0 * gx)**2 + (3.0 * gy)**2)
 
print("%-6s %-12s %-14s %-14s" % ("W", "stencil", "mean err %", "anisotropy %"))
for W in (5.0, 4.0, 3.0):
    X, Y, r, phi, exact = fields(W)
    ring = np.abs(r - R) < 0.6
    theta = np.arctan2(Y, X)[ring]
    for name, g in (("central", grad_central(phi)), ("D2Q9", grad_lattice(phi))):
        err = (g[ring] - exact[ring]) / exact[ring] * 100.0
        b = ((theta + np.pi) / (np.pi / 12)).astype(int) % 24
        binned = np.array([err[b == j].mean() for j in range(24) if np.any(b == j)])
        print("%-6.1f %-12s %-14.3f %-14.3f" % (W, name, np.abs(err).mean(),
                                                binned.max() - binned.min()))
W      stencil      mean err %     anisotropy %  
5.0    central      3.641          1.812         
5.0    D2Q9         4.675          0.075         
4.0    central      5.338          2.575         
4.0    D2Q9         6.803          0.155         
3.0    central      8.369          3.749         
3.0    D2Q9         10.539         0.402

ステンシルのボタンを中心差分からD2Q9へ切り替えると、右のローズが四枚の花びらから円へ折れます。

anisotropy = 0.000%mean |err| = 0.000%

Switch the stencil and watch the rose, not the readout beside it. The central-difference rose is a four-lobed flower whose petals point along the lattice axes; the D2Q9 rose is almost a circle even though its mean error is the larger of the two. Shrink W towards 2 to see both get worse, but only one of them get lopsided.

WW を4に置くと、異方性2.575%と0.155%がその差です。2の方へ下げると両方とも悪化しますが、いびつになるのは中心差分だけです。

W=4W=4 でD2Q9の平均誤差は6.803%、中心差分の5.338%より大きい。精度だけ見れば負けています。ところが角度による広がりは2.575%から0.155%へ、17倍近く縮みます。2.575/0.155=16.62.575/0.155 = 16.6 です。平均誤差の差は 6.803−5.338=1.4656.803 - 5.338 = 1.465、約1.5ポイントです。W=5W=5 では異方性が1.812から0.075へ、W=3W=3 では3.749から0.402へ縮みます。平均誤差は逆向きに増えます。中心差分は3.641、5.338、8.369。D2Q9は4.675、6.803、10.539です。薄くするほど等方性の得は残り、大きさの誤差は W=5W=5 に対してどちらも2倍を超えます。

この2つは大きさが違うのではなく、性質が違います。界面全体に一様に乗った誤差は ∇ϕ\nabla\phi の大きさを一定比率で外します。表面張力係数をその分だけ取り違えたのと同じで、静止した液滴は静止したままです。角度によって変わる誤差は別物です。同じ円の上の点ごとに違う大きさの力がかかります。その不均衡が、止まっているはずの液滴表面に流れを生みます。よく言われる寄生流(spurious current)です。

つまりここでは平均誤差1.5ポイントを渡して等方性17倍を買っています。割に合う取引です。

WW を決めるとき実際に値切っているもの#

表の3列目と4列目は逆向きに動きます。WW を5から3へ下げると界面は物理的に鋭くなりますが、∇ϕ\nabla\phi の平均誤差は3.6%から8.4%へ、中心差分の異方性は1.8%から3.7%へ増えます。WW を上げれば誤差は減りますが、界面が鈍って小さな液滴や薄い膜を解像できなくなります。

W=4∼5W=4\sim5 が実務の定番になった理由がこの表にあります。界面を薄くしたいときに MM をいじっても何も起きません。冒頭の4.44と4.24がその話です。MM を回すと回復速度だけが変わり、下げすぎると移流が生む数値拡散に反拡散が追いつかなくなります。厚みのつまみは WW ひとつです。三相の張力を互いにずらしたとき界面が漏れる場所は、対ごとの張力がずれる三相Allen–Cahnに別にあります。

関連記事

静止した液滴が半分に溶けた — 保存形 Allen–Cahn が消す曲率項流速ゼロの液滴で古典 Allen–Cahn は質量を49.4%失い、保存形は1e-14に留まりました。その差を生む曲率項をLBのソース項まで追います。
Oct 9, 2026·7 min read
σ₁₂に1.00を入れて0.78が返ってきた — 3成分Allen–Cahnで存在しない相が界面に漏れる場所3成分Allen–CahnでΣc=1を平均投影で合わせると、存在しない相が界面に漏れて表面張力が21.9%削られる。3重罰則項がその穴を塞ぐ。
Oct 2, 2026·7 min read
[論文レビュー] 毛細管の足枷を外したらCFL 0.05が残った — VOF界面移流の本当の上限表面張力を陰的に解いて毛細管時間刻み制約を破っても、色関数を運ぶ圧縮スキームがCFL 0.05で時間刻みを掴み直します。
Aug 16, 2026·8 min read

役に立ったらシェアしてください。

GitHub でソースを見る ↗誤字・改善提案 ↗