σ₁₂填了1.00，拿回0.78 — 三组分Allen–Cahn中缺席相渗入界面的位置

σ₁₂填了1.00，拿回0.78#

我在一个三相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三角形闭合的条件是同一个不等式，并找出阻止缺席相渗入界面的那一项。最后在同一套网格上跑四种闭合方式，用数字对比。

N个体积率不是N个自由度#

处理$N$个互不相溶的流体时，体积率$c_i$带着一个约束。

自由度是$N-1$个。可是写程序时会开出$N$个数组，给每个数组挂一个Allen–Cahn方程看起来最自然。把两相程序复制过来只增加组分数，得到的正是这个形状。

于是总和就坏了。把各相的化学势$\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$都不是可以随意选的量。两相极限把它们绑到物理量上。代入$c_2 = 1 - c_1$，上式变为

与熟悉的两相形式$\beta c^2(1-c)^2 + \frac{\lambda}{2}|\nabla c|^2$对比，得$\lambda = \kappa_1 + \kappa_2$。把一维平衡$\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$必须把三个式子联立才出来。三相时

一个按组分的系数，是三个两两表面张力的组合。这是本文的中心式子。在两相守恒型Allen–Cahn中单独处理迁移率与反扩散项的那篇里只有一个$\lambda$。相变成三个之后，$\lambda$分成三个，而分法是减法。

$D=4$时系数$\frac{3D}{4}$等于$3$。表里第二行，$\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$。两相关系$\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三角形是否闭合与模型是否良定是同一个不等式。固定$\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的角正好是直角。

在同一套网格上跑四种闭合

一维网格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彻底冻结，是两相计算。$\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}$取什么值都无关。在两两四次井里，往二元界面加入第三相的方向永远是下坡。$\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$。

三重罚项Λ的位置

处方是直接惩罚三相在一点重叠的项。

从扫描看，$\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$，三角形不闭合。五相时这样的三元组有十个。只要有一个为负，那三个相相遇的单元里能量就向下张开。三个$\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$当作基准，这个差不会进日志。