정지한 물방울이 절반으로 녹았다 — 보존형 Allen

보존형 Allen–Cahn의 샤프닝 항은 계면 두께를 지키는 장치가 아니라, 곡률 구동 소멸을 정확히 상쇄하는 항이다.

유속이 0인데 물방울이 절반으로 줄었다#

곡률만 있고 유속은 0인 계면은 저절로 움직여야 하는가. phase-field 모델(계면을 두께가 있는 전이층으로 다루는 모델)을 코드에 처음 넣을 때 이 질문은 잘 빠진다. 정지 상태로 초기화한 2차원 물방울을 Allen–Cahn 방정식으로 4000 스텝 돌리면 답이 나온다. 반지름 20.083이 14.286으로 줄고, 체적율의 총합은 49.4% 사라진다. 속도장은 끝까지 정확히 0이다.

이 글은 그 소멸이 어느 항에서 오는지 항 단위로 짚는다. 고전 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$	보존되지 않음	이산 수준에서 보존
1D 평형해	tanh	같은 tanh
곡률 응답	$v_n=-M\kappa$	0

두 식은 이렇게 적힌다.

\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$ 은 이동도(mobility, 계면 이완 속도를 정하는 계수), $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$ 가 왜 저 두 식에 다른 모양으로 들어가는지는 1D 평형해에서 나온다. 평면 계면의 해는 다음 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$ 의 출처다. 두 식을 각 방정식에 넣으면 우변이 0이다. 1D에서는 두 모델이 같은 평형해를 가진다.

곡률 항이 남느냐 지워지느냐#

차이는 계면이 휘는 순간 생긴다. 법선 방향 좌표 $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$ 를 따른다. 보존형은 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 전개(평형에서의 작은 편차를 격자 간격 비로 차수별로 정리하는 방법)가 요구하는 조건은 모멘트 두 줄이다. 0차 모멘트는 $\sum_i F_i = 0$ 이어야 한다. $\phi=\sum_i h_i$ 이므로, 0차 모멘트가 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 평형이 다항식으로 닫히는 자리와 같은 종류의 조건이고, 여기서는 지켜야 할 줄이 두 개로 끝난다.

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차 모멘트는 소수 여섯째 자리까지 같다. 0차 모멘트만 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}$ 계산 필요

표의 세 번째 줄이 중요하다. 보존형이 항상 옳은 것이 아니다. 곡률 구동 운동을 모사하려는 상황이면 고전형이 맞는 모델이다. 유체 계면을 이송하려는 상황에서만 그 항이 오차가 된다. 어느 쪽이든 $\hat{\mathbf n}$ 을 계산해야 하므로, 보존형에서는 $|\nabla\phi|$ 가 0에 가까운 벌크 영역의 나눗셈 보호가 추가 작업으로 남는다.

계면 포착기 코드에 W를 적을 때#

코드에 W = 4를 적는 순간 정하는 것은 두께 하나가 아니다. 이중우물 계수 $16/W^2$ 와 역플럭스 계수 $4/W$ 가 같은 숫자에서 나온다. 둘 중 하나만 바꾸면 평형해가 tanh에서 벗어난다. 두 계수를 같은 $W$ 에서 유도하는지가 첫 점검 항목이다.

그 다음 점검은 보존량이다. 체적율 총합을 매 스텝 찍어 보면 1e-14와 1e-1은 한눈에 갈린다. 두께는 양쪽 모델에서 모두 정상으로 보이므로 두께만 보면 아무것도 걸러지지 않는다. 이동도를 10배 올려도 계면 두께가 그대로였던 보존형에서 이미 두께와 질량은 다른 약속이었다. LB 구현이라면 점검은 한 줄로 끝난다. $\sum_i F_i$ 를 출력해 0인지 확인하는 것이다.

정지한 물방울이 절반으로 녹았다 — 보존형 Allen–Cahn이 지우는 곡률 항