σ₁₂에 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 삼각형의 닫힘 조건과 같은 식이라는 사실, 그리고 없는 상이 계면으로 새는 것을 막는 항 하나를 본다. 끝에는 같은 격자에 닫힘 방식 네 가지를 걸어 숫자로 비교한다.

체적율 N개는 자유도 N개가 아니다#

$N$개의 불혼합 유체를 다룰 때 체적율 $c_i$는 제약 하나를 들고 있다.

자유도는 $N-1$개다. 그런데 코드를 짜면 배열은 $N$개가 생기고, 각 배열에 Allen–Cahn 방정식 하나씩을 붙이는 것이 가장 자연스럽다. 두상 코드를 복사해서 성분 수만 늘리면 정확히 그 모양이 된다.

그러면 총합이 깨진다. 각 상의 화학포텐셜 $\mu_i$를 독립적으로 쓰고 $\partial_t c_i = -M \mu_i$를 적으면

이고, 오른쪽이 0일 이유가 없다. 실험에서 $\left|\sum c_i - 1\right|$이 0.859까지 벌어졌다. 계면이 녹아버리고 표면장력은 0.026으로 내려앉는다.

교과서적 처방은 Lagrange 승수다. $\bar\mu = \frac{1}{N}\sum_j \mu_j$를 빼면

합의 시간변화가 정확히 0이 된다. 아래 시뮬레이션에서 두 닫힘을 직접 갈아 끼워보자.

왼쪽 버튼으로 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$다. 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상이면

성분별 계수 하나가 쌍별 표면장력 셋의 조합이다. 이것이 이 글의 중심 식이다. 두상 보존형 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 삼각형이 닫히는가와 모델이 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의 각이 직각이다.

같은 격자에 닫힘 네 가지를 걸었다#

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을 아예 동결한 두상 계산이다. $\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$에서 0이다. $-\beta_{12}/4$라는 내리막 기울기는 그대로 남는다. $\Lambda$가 하는 일은 봉우리가 자란 뒤에 벽을 세우는 것이다. 그래서 $\Lambda$를 올려도 $c_3$이 0으로 가지 않고 0.05에서 멈춘다.

구멍 자체를 막으려면 자유도를 $N-1$개로 줄이고 $c_N = 1 - \sum_{i<N} c_i$로 닫거나, 국소적으로 없는 상의 방정식을 아예 풀지 않아야 한다. A가 바로 그 경우고, 그래서 $\max c_3$이 정확히 0이다.

적었던 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$을 기준으로 남기면 그 차이가 로그에 안 남는다.