σ₁₂に1.00を入れて0.78が返ってきた — 3成分Allen–Cahnで存在しない相が界面に漏れる場所

σ₁₂に1.00を入れて0.78が返ってきた#

3相のphase-fieldコードに相1、相2、相3を入れました。1-2の表面張力の欄に1.00と書きました。平衡に達した界面を横切って自由エネルギーを積分しました。出てきたのは0.78でした。

抜けていたのは相3です。計算領域のどこにも相3はありません。初期値として置いた量は$10^{-3}$だけです。それでも平衡界面の中央で相3の体積率が0.21まで立ち上がっていました。

この記事はその21%をどこまでも追いかけます。対ごとに書いた$\sigma_{ij}$が成分ごとの勾配係数$\kappa_i$へ戻る計算、$\kappa_i$の符号がNeumann三角形の閉じる条件と同じ不等式であること、そして存在しない相が界面に漏れるのを止める項を見ます。最後に同じ格子へ閉じ方を4通り掛けて数字で比べます。

体積率N個は自由度N個ではない#

$N$個の非混和流体を扱うとき、体積率$c_i$は拘束を一つ抱えています。

自由度は$N-1$個です。ところがコードを書くと配列は$N$個でき、各配列にAllen–Cahn方程式を一つずつ付けるのがいちばん自然に見えます。2相コードをコピーして成分数だけ増やすと、まさにその形になります。

すると総和が壊れます。各相の化学ポテンシャル$\mu_i$を独立に書いて$\partial_t c_i = -M \mu_i$で進めると

となり、右辺がゼロになる理由はありません。実験では$\left|\sum c_i - 1\right|$が0.859まで開き、界面が溶けて表面張力は0.026まで落ちました。

教科書的な処方はLagrange乗数です。$\bar\mu = \frac{1}{N}\sum_j \mu_j$を引くと

となり、総和の時間変化が厳密にゼロになります。下のシミュレーションで二つの閉じ方を自分で入れ替えてみてください。

左のボタンで3 equations, no closureを選ぶと、$\sum c_i - 1$の帯が下へ崩れながら三本の曲線が一斉に潰れます。mean projectionに変えると帯は$10^{-13}$に張り付きますが、赤い$c_3$の山が界面の中央に立ち上がり、sigma_effが1.00を下回ります。Lambdaスライダーを上げるとその山が押さえられます。

対ごとのσᵢⱼを成分ごとのκᵢへ戻す

自由エネルギーは普通、対ごとに書きます。

$\beta_{ij}$は$i$-$j$対の二重井戸の深さ、$\kappa_i$は成分$i$の勾配エネルギー係数です。罰則は対ごと、勾配係数は成分ごと。この非対称が後ろのすべてを作ります。

$\beta_{ij}$も$\kappa_i$も自由に選べる値ではありません。2相極限で物理量に縛られます。$c_2 = 1 - c_1$を入れると上の式は

になり、見慣れた2相形$\beta c^2(1-c)^2 + \frac{\lambda}{2}|\nabla c|^2$と比べると$\lambda = \kappa_1 + \kappa_2$です。1次元平衡$\lambda c'' = dE_0/dc$を一度積分すると$\frac{\lambda}{2}c'^2 = E_0$、つまり$c' = \sqrt{2\beta/\lambda}\, c(1-c)$。解はtanhで

ここで$D$は界面厚さ(格子セル数)、$\sigma$は表面張力です。二つの関係を反転すると

$\beta_{ij}$は対ごとに即座に決まります。$\kappa_i$は三つの式を連立しないと出ません。3相なら

成分ごとの係数一つが、対ごとの表面張力三つの組み合わせです。これがこの記事の中心の式です。2相の保存型Allen–Cahnで移動度と反拡散項を扱った記事では$\lambda$は一つだけでした。相が三つになると$\lambda$は三つに分かれ、その分かれ方が引き算です。

$D=4$なら係数$\frac{3D}{4}$は$3$です。表の2行目、$\sigma$が$1.00$、$0.80$、$0.60$なら$\kappa_1$は$3$掛ける$1.20$、つまり$3.60$です。$\kappa_2$は$3$掛ける$0.80$で$2.40$、$\kappa_3$は$3$掛ける$0.40$で$1.20$です。三つを足すと$7.20$、1-2の対だけ足すと$6.00$です。2相の関係$\kappa_i+\kappa_j=\frac{3}{2}\sigma_{ij}D$に$1.00$と$4$を入れても$6.00$が出ます。$\lambda=6$を$\sigma=\frac{2\lambda}{3D}$に入れると$\frac{12}{12}=1.00$です。連続式は入力と一致します。格子でAが返した$0.9665$との$3.3\%$は、界面4セルが作った差です。同じ$3$を$\sigma=(1.00,\,0.40,\,0.40)$に掛けると$\kappa_3=-0.60$です。表で三角形が閉じない行が、この引き算です。$(1.00,\,0.55,\,0.55)$は$\kappa_3=0.30$で符号が生きていて、くさびは$155.4^\circ$、$155.4^\circ$、$49.2^\circ$です。負へ渡る前の余裕が$0.30$で、相3の角はすでに$49^\circ$まで狭まっています。$\kappa_3=1.20$は$\sigma_{12}$の帳簿に入りません。界面に相3が現れてはじめて勾配エネルギーになります。等しい$\sigma=(1.00,\,1.00,\,1.00)$なら三つの$\kappa$はすべて$3.00$、くさびは$120^\circ$三つで、和は$360^\circ$です。$0.55$の行の$49.2^\circ$は、その$120^\circ$より$70.8^\circ$狭いです。符号が死ぬ前に角が先に狭まります。

κ₃が負へ渡る場所

引き算ということは符号を失えるという意味です。$\sigma_{13} + \sigma_{23} < \sigma_{12}$なら$\kappa_3 < 0$になります。勾配エネルギーが下に開きます。界面を鋭くするほどエネルギーが下がるので、解は格子間隔まで崩れます。モデルが劣化するのではなく、積分がはじめから下に有界でないのです。

同じ条件は別の名前でも知られています。三重点で三つの表面張力が釣り合うと、各相のくさび角$\theta_i$が決まります。

$\theta_1$は相1が占める角で、残り二つは添字を回して得ます。この$\cos$が$[-1, 1]$を外れる条件が、まさに$\kappa_i < 0$の条件です。つまりNeumann三角形が閉じるかとモデルがwell-posedかは同じ不等式です。$\sigma_{12}=1.00$を固定して$\sigma_{13}+\sigma_{23}$が$\sigma_{12}$を下回ると、三角形が切れて相3のくさびが180°まで開きます。相3が1-2界面を完全に濡らす状態です。$\sigma=(1.00, 0.80, 0.60)$ではくさび角が143.1° / 126.9° / 90.0°です。3-4-5の直角三角形なので、相3の角がちょうど直角です。

同じ格子に閉じ方を4通り掛けた#

1次元格子200セル、$D=4$、$\sigma=(1.00, 0.80, 0.60)$。左は相1、右は相2、相3は界面に$10^{-3}$だけ置きます。8000ステップ緩和させた後、$\left|\sum c_i - 1\right|$、$\max c_3$、そして界面を横切る自由エネルギー積分(=実現された表面張力)を測ります。

import math
import numpy as np
 
D, MOB, DT, NX = 4.0, 1.0 / 3.0, 0.2, 200
 
 
def pair_to_component(s12, s13, s23):
    """対ごとの表面張力から成分ごとの勾配係数 kappa_i と対の罰則 beta_ij へ。"""
    f = 0.75 * D
    kap = np.array([f * (s12 + s13 - s23),
                    f * (s12 + s23 - s13),
                    f * (s13 + s23 - s12)])
    bet = 12.0 / D * np.array([[0, s12, s13], [s12, 0, s23], [s13, s23, 0]])
    return kap, bet
 
 
def neumann_wedges(s12, s13, s23):
    """Neumann 三角形が閉じれば三相のくさび角(度)、閉じなければ None。"""
    args = [(s23 ** 2 - s12 ** 2 - s13 ** 2) / (2 * s12 * s13),
            (s13 ** 2 - s12 ** 2 - s23 ** 2) / (2 * s12 * s23),
            (s12 ** 2 - s13 ** 2 - s23 ** 2) / (2 * s13 * s23)]
    if any(abs(a) > 1 for a in args):
        return None
    return [math.degrees(math.acos(a)) for a in args]
 
 
def _ext(a):
    return np.concatenate(([a[0]], a, [a[-1]]))
 
 
def lap1d(a):
    e = _ext(a)
    return e[:-2] - 2 * a + e[2:]
 
 
def grad1d(a):
    e = _ext(a)
    return (e[2:] - e[:-2]) / 2.0
 
 
def chem_potential(c, kap, bet, lam3):
    c1, c2, c3 = c
    dbulk = np.array([
        2 * bet[0, 1] * c1 * c2 ** 2 + 2 * bet[0, 2] * c1 * c3 ** 2 + 2 * lam3 * c1 * (c2 * c3) ** 2,
        2 * bet[0, 1] * c2 * c1 ** 2 + 2 * bet[1, 2] * c2 * c3 ** 2 + 2 * lam3 * c2 * (c1 * c3) ** 2,
        2 * bet[0, 2] * c3 * c1 ** 2 + 2 * bet[1, 2] * c3 * c2 ** 2 + 2 * lam3 * c3 * (c1 * c2) ** 2,
    ])
    return np.array([dbulk[i] - kap[i] * lap1d(c[i]) for i in range(3)])
 
 
def line_tension(c, kap, bet, lam3):
    """界面を横切る自由エネルギー積分 = 実現された表面張力。"""
    c1, c2, c3 = c
    e = (bet[0, 1] * (c1 * c2) ** 2 + bet[0, 2] * (c1 * c3) ** 2
         + bet[1, 2] * (c2 * c3) ** 2 + lam3 * (c1 * c2 * c3) ** 2)
    for i in range(3):
        e = e + 0.5 * kap[i] * grad1d(c[i]) ** 2
    return float(e.sum())
 
 
def initial_profile(seed):
    x = np.arange(NX) - NX / 2
    c3 = seed * np.exp(-(x / D) ** 2)
    c1 = 0.5 * (1 - np.tanh(2 * x / D)) * (1 - c3)
    c2 = 0.5 * (1 + np.tanh(2 * x / D)) * (1 - c3)
    return np.array([c1, c2, c3])
 
 
def relax(kap, bet, closure, lam3=0.0, steps=8000, seed=1e-3):
    c = initial_profile(0.0 if closure == 'two-phase' else seed)
    for _ in range(steps):
        mu = chem_potential(c, kap, bet, lam3)
        if closure == 'two-phase':
            m = 0.5 * (mu[0] + mu[1])
            rhs = np.array([mu[0] - m, mu[1] - m, np.zeros(NX)])
        elif closure == 'independent':
            rhs = mu
        else:
            rhs = mu - mu.mean(axis=0)
        c = c - MOB * DT * rhs
        if not np.isfinite(c).all() or np.abs(c).max() > 1e3:
            return None, None, None
    return c, float(np.abs(c.sum(axis=0) - 1.0).max()), line_tension(c, kap, bet, lam3)
 
 
print("== pairwise sigma_ij -> per-component kappa_i (D = %.1f) ==" % D)
print("%-20s %-26s %-8s %s" % ("sigma_12,13,23", "kappa_1,2,3", "closes", "wedge angles [deg]"))
for s in [(1.0, 1.0, 1.0), (1.0, 0.8, 0.6), (1.0, 0.55, 0.55), (1.0, 0.4, 0.4)]:
    kap, _ = pair_to_component(*s)
    th = neumann_wedges(*s)
    print("%-20s %-26s %-8s %s" % (
        "%.2f %.2f %.2f" % s,
        "%+.2f %+.2f %+.2f" % tuple(kap),
        "yes" if kap.min() > 0 else "NO",
        "-" if th is None else "%.1f %.1f %.1f (sum %.1f)" % (th[0], th[1], th[2], sum(th))))
 
kap, bet = pair_to_component(1.0, 0.8, 0.6)
print()
print("== 1-2 interface at equilibrium, phase 3 seeded at 1e-3 (input sigma_12 = 1.00) ==")
print("%-30s %-13s %-9s %s" % ("closure", "max|sum c-1|", "max c_3", "sigma_eff"))
cases = [("A two-phase (c_3 frozen)", 'two-phase', 0.0),
         ("B three equations, no closure", 'independent', 0.0),
         ("C mean projection", 'projected', 0.0),
         ("D mean projection + Lambda=400", 'projected', 400.0)]
for tag, closure, lam3 in cases:
    c, drift, sig = relax(kap, bet, closure, lam3)
    if c is None:
        print("%-30s diverged" % tag)
        continue
    print("%-30s %-13.2e %-9.4f %.4f" % (tag, drift, c[2].max(), sig))
 
print()
print("== triple-penalty Lambda scan (mean projection) ==")
print("%-10s %-10s %-10s %s" % ("Lambda", "max c_3", "sigma_eff", "vs input"))
for lam3 in (0.0, 10.0, 30.0, 60.0, 150.0, 400.0):
    c, drift, sig = relax(kap, bet, 'projected', lam3)
    print("%-10.0f %-10.4f %-10.4f %+.1f%%" % (lam3, c[2].max(), sig, (sig - 1.0) * 100))

== pairwise sigma_ij -> per-component kappa_i (D = 4.0) ==
sigma_12,13,23       kappa_1,2,3                closes   wedge angles [deg]
1.00 1.00 1.00       +3.00 +3.00 +3.00          yes      120.0 120.0 120.0 (sum 360.0)
1.00 0.80 0.60       +3.60 +2.40 +1.20          yes      143.1 126.9 90.0 (sum 360.0)
1.00 0.55 0.55       +3.00 +3.00 +0.30          yes      155.4 155.4 49.2 (sum 360.0)
1.00 0.40 0.40       +3.00 +3.00 -0.60          NO       -
 
== 1-2 interface at equilibrium, phase 3 seeded at 1e-3 (input sigma_12 = 1.00) ==
closure                        max|sum c-1|  max c_3   sigma_eff
A two-phase (c_3 frozen)       1.78e-15      0.0000    0.9665
B three equations, no closure  8.59e-01      0.0000    0.0257
C mean projection              2.59e-13      0.2114    0.7809
D mean projection + Lambda=400 8.77e-15      0.0512    0.9715
 
== triple-penalty Lambda scan (mean projection) ==
Lambda     max c_3    sigma_eff  vs input
0          0.2114     0.7809     -21.9%
10         0.1887     0.8031     -19.7%
30         0.1448     0.8377     -16.2%
60         0.0963     0.8704     -13.0%
150        0.0515     0.9111     -8.9%
400        0.0512     0.9715     -2.9%

三つの$\sigma$スライダーを動かすと、$\kappa$の棒の符号とくさび角が一緒に変わります。$(1.00,\,0.40,\,0.40)$まで下げると、表のcloses列がひっくり返ります。

$\sigma_{12}$を$1.00$に固定して$\sigma_{13}$と$\sigma_{23}$を下げると、赤い$\kappa_3$が$0$を越えた瞬間に三角形が切れます。$(1.00,\,0.40,\,0.40)$の行の$\kappa_3=-0.60$がその場所です。$(1.00,\,0.80,\,0.60)$では角が$143.1^\circ$ / $126.9^\circ$ / $90.0^\circ$のままです。

Σcは1e-13なのに界面が21.9%軽くなった#

Aは相3を完全に凍結した2相計算です。$\sigma_{\rm eff} = 0.9665$ — 入力1.00に対して3.3%低い。これはモデルの欠陥ではなく、界面4セルの離散化誤差です。$\kappa_1 + \kappa_2 = \frac{3}{2}\sigma_{12}D$の関係が正しく成り立っている確認になります。

Bは三つの式を独立に解いた場合です。総和が0.859ずれ、表面張力は0.0257まで消えます。界面が溶けました。

Cが問題の行です。平均投影で総和は$2.6\times10^{-13}$まで押さえました。ところが$\max c_3 = 0.2114$、$\sigma_{\rm eff} = 0.7809$です。入力に対して21.9%低い。拘束は守られて物理が間違っています。

拘束違反は目に見えます。21.9%は見えません。だからこちらのほうが高くつきます。

存在しない相が界面に漏れる算数

なぜ漏れるかは、バルクエネルギーだけ切り出すとすぐ見えます。1-2界面の中央で対称を仮定し、$c_3 = p$、$c_1 = c_2 = (1-p)/2$を入れます。

$p$で微分して$p=0$で値を見ます。第二項は$p^2$から始まるので傾きには寄与しません。残るのは第一項だけです。

常に負です。$\sigma_{13}$や$\sigma_{23}$がいくつでも関係ありません。対ごとの4次井戸では、二元界面に第三の相を入れる方向が常に下り坂です。$\beta_{12} = 3.0$なら傾きは$-0.75$です。

実際の数字でも確認できます。$p=0$で$E_0 = \beta_{12}/16 = 0.1875$。測定された$p = 0.2114$を入れると$E_0 = 0.1018$で、45.7%下がります。バルクは相3を歓迎します。止めてくれるのは$\kappa_3 |\nabla c_3|^2$だけで、$\kappa_3 = 1.20$では0.21でようやく釣り合います。

止まった液滴が勝手に流れ出す寄生流と似たところがあります。張力側の離散化がずれると、総和は合っていても界面が静かに間違います。こちらは$\sum c_i$が$2.59\times10^{-13}$なのに$\sigma_{\rm eff}$だけが$0.7809$です。

3重罰則項Λが住む場所#

処方は、三相が一点で重なることを直接罰する項です。

スキャンを見ると$\Lambda$が上がるほど$\max c_3$が0.2114から0.0512へ押さえられ、$\sigma_{\rm eff}$は0.7809から0.9715へ戻ります。$\Lambda = 400$で入力比−2.9%、つまりAの離散化誤差−3.3%とほぼ同じ場所です。相3を凍結したものと区別がつきません。

ただしこの項は漏れ穴を塞ぎません。$\Lambda c_1^2c_2^2c_3^2$を$p$で微分しても$p=0$ではゼロです。$-\beta_{12}/4$という下り傾きはそのまま残ります。$\Lambda$がするのは山が育った後に壁を立てることです。だから$\Lambda$を上げても$c_3$はゼロへ行かず0.05で止まります。

穴そのものを塞ぐには自由度を$N-1$個に減らして$c_N = 1 - \sum_{i<N} c_i$で閉じるか、局所的に存在しない相の方程式をそもそも解かないことです。Aがまさにその場合で、だから$\max c_3$がちょうどゼロになります。

書いた1.00を信じるか、積分した0.7809を信じるか#

対の表を入れた瞬間に、$\kappa_i = \frac{3D}{4}(\sigma_{ij}+\sigma_{ik}-\sigma_{jk})$の符号とくさび角が一緒に固定されます。$\sigma=(1.00,\,0.80,\,0.60)$なら$\kappa$は$+3.60$、$+2.40$、$+1.20$で、角は$143.1^\circ$ / $126.9^\circ$ / $90.0^\circ$、和は$360^\circ$です。$\sigma=(1.00,\,0.40,\,0.40)$なら$\kappa_3=-0.60$で、三角形は閉じません。5相ならこのような三重組は10個です。一つでも負なら、その三つが出会うセルでエネルギーが下に開きます。三つの$\cos\theta_i$の和が$360^\circ$かどうかは、$\kappa_i$が正かどうかと同じ問いです。角は図と直接重ねられます。

検証に掛ける値は入力の$\sigma_{12}=1.00$ではありません。Cの積分値$0.7809$です。総和は$2.59\times10^{-13}$で合っているのに、張力は$21.9\%$小さいです。表面張力を陰的に扱っても毛管の時間刻みが残る場合と同じく、この$21.9\%$はLaplace圧と毛管波速度にそのまま乗ります。$\sum c_i=1$だけ見れば合格です。$\Lambda=400$まで上げても残差は$-2.9\%$の$0.9715$で、相3を切ったAの$0.9665$との差は$0.0050$です。入力欄の$1.00$を基準に残すと、その差はログに出ません。