이동도를 10배 키웠는데 계면 두께는 그대로였다 — 보존형 Allen–Cahn의 반확산 항

확산을 10배로 키웠는데 계면이 두꺼워지지 않았다#

상장(phase field, 계면을 0과 1 사이의 부드러운 스칼라로 표현하는 방법) 코드에서 이동도 $M$을 0.02에서 0.20으로 올렸다. 확산 계수를 열 배 키운 셈이다. 계면은 열 배 빨리 뭉개져야 한다. 측정된 두께는 4.44에서 4.24로, 오히려 미세하게 줄었다.

계산이 틀린 것이 아니다. 방정식이 그렇게 생겼다. 격자 볼츠만(LBM)에서 계면을 붙잡는 보존형 Allen–Cahn은 두께를 $W$에 묶고, 반확산은 평형분포 한 줄로 들어간다. $\nabla\phi$를 격자에서 어떻게 구하는지도 같은 숫자에 남는다. 아래 숫자는 실은 코드를 돌려 얻은 것이다.

smear interface를 누르면 계면이 두꺼워졌다가, 반확산이 켜진 쪽만 초록 점선 두께로 돌아온다.

smear interface를 눌러 계면을 일부러 뭉갠 다음 아래쪽 주황색 두께 기록을 본다. 반확산이 켜져 있으면 어떤 $M$에서도 초록 점선 높이로 되돌아오고, 꺼져 있으면 $M$은 두께가 커지는 속도만 정한다.

계면을 방정식 하나로 붙잡는다는 발상#

VOF나 level set은 이송을 푼 뒤 계면을 재구성하거나 거리함수를 재초기화한다. 별도의 단계다. 같은 보존형 Allen–Cahn이라도 곡률 항을 빼면 정지한 물방울이 녹는다. 그 자리는 곡률 항이 지워지는 보존형 Allen–Cahn에 있다.

보존형 Allen–Cahn은 그 후처리를 같은 방정식 안으로 집어넣는다.

$\phi$는 상장 변수(한쪽 상에서 0, 다른 쪽에서 1), $M$은 이동도, $W$는 목표 계면 두께(격자 칸 수), $\mathbf{n}=\nabla\phi/|\nabla\phi|$는 계면 법선이다. 괄호 안 첫 항은 계면을 퍼뜨리는 보통의 확산이고, 둘째 항이 반대 방향으로 밀어 넣는 반확산이다.

우변 전체가 발산 형태라는 점이 중요하다. 주기 경계에서 $\int\phi\,dV$는 기계 정밀도까지 보존된다. Allen–Cahn이 원래 질량을 보존하지 않는 방정식인데도 "보존형"이라 부르는 이유가 여기 있다.

반확산 항이 정확히 0을 만드는 자리#

계면 법선 방향 좌표를 $\xi$라 하면 평형 프로파일은 이렇게 생겼다.

미분하면 기울기가 나온다.

한편 같은 프로파일에서 $4\phi(1-\phi) = 1-\tanh^2(2\xi/W) = \mathrm{sech}^2(2\xi/W)$이다. 두 식이 같은 함수다. 즉 평형 프로파일 위에서는

가 항등적으로 성립한다. 괄호 안이 0이다. 그리고 $M$은 그 0에 곱해진다.

여기서 갈라지는 것은 두 가지다. 평형 두께는 $W$ 하나로 닫히고 $M$은 그 자리에 끼지 못한다. $M$이 정하는 것은 평형에서 밀려났을 때의 회복 시간뿐이고, 그 시간은 $W^2/M$ 정도로 움직인다.

격자 위에서는 그 0이 저절로 유지되지 않는다#

연속식에서 정확히 상쇄되는 두 항도 이산화하면 따로 논다. 확산과 반확산을 각각 격자점에서 중심차분으로 계산하면 두 스텐실이 모두 $i-1$과 $i+1$만 보고 $i$를 건너뛴다. 홀수 점과 짝수 점이 서로 모르는 상태가 되고, 그 격자 규모 진동이 반확산에 증폭된다. 실제로 그렇게 짜서 돌렸더니 $M=0.20$에서 두께가 목표 4.0을 지나쳐 2.88까지 내려갔다.

해결은 플럭스를 격자점이 아니라 면에서 계산하는 것이다. 면 $i+1/2$의 플럭스를 $\phi_i$와 $\phi_{i+1}$ 두 값만으로 만들면 이웃한 점이 서로를 본다. 홀짝 분리가 생길 자리가 없다.

같은 자를 대고 두께를 여덟 번 쟀다#

슬래브 하나를 두께 8로 초기화한 뒤 목표 $W=4$로 8,000스텝 돌렸다. 반확산을 켠 경우와 끈 경우를 $M$ 네 값에 대해 비교한다.

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

첫 줄의 ruler가 눈금이다. 두께를 $1/\max|\nabla\phi|$로 재는데, 폭 4짜리 tanh를 격자에 얹어 같은 방법으로 재면 4.0이 아니라 4.3279가 나온다. 자 자체의 이산화 오차다.

그 눈금에 대면 반확산을 켠 네 값 4.24~4.44는 전부 $\pm3\%$ 안이다. $M$이 10배 움직이는 동안 두께는 사실상 고정되어 있다. 반확산을 끄면 같은 $M$ 범위에서 21.9에서 64.0으로 벌어진다. 질량 표류는 양쪽 다 $10^{-11}$ 수준이다.

$W=4$, $M=0.02$이면 $W^2/M = 16/0.02 = 800$이다. $M=0.20$이면 80이다. 코드의 8000스텝에 $dt=0.2$를 곱하면 적분 시간은 1600이다. 작은 이동도에서도 회복 시간 두 번, 큰 이동도에서는 스무 번을 덮는다. 그래서 ON 네 줄이 4.4403, 4.2854, 4.2442, 4.2357로 모인다. 자 자체는 폭 4의 tanh를 4.3279로 읽는다. $1/4.3279$는 약 0.231로, 연속식의 최대 기울기 $1/W = 0.25$보다 작다. 4.24~4.44가 $\pm3\%$ 안이라는 말은 4.0이 아니라 이 4.3279를 기준으로 한 것이다.

마지막 줄이 이것이 단순한 "변화 없음"이 아님을 보여준다. 두께 2로 시작해도 4.2854로 올라온다. 위에서 내려온 값과 소수점 넷째 자리까지 같다. 평형 프로파일은 양쪽에서 끌어당기는 끌개다.

이 항을 LBM으로 옮기면 소스항이 아니다#

$\phi$에 분포함수 $h_i$를 따로 세운다. 여기서 흔한 오해가 하나 있다. 반확산을 충돌 후에 더하는 소스항으로 짜는 것이다. 그럴 필요가 없다.

평형분포를 이렇게 잡는다.

$w_i$는 격자 가중치, $c_s^2=1/3$은 격자 음속 제곱, $\tau_\phi$는 무차원 완화시간이다. 대괄호 안이 평형분포의 1차 모멘트, 즉 $\phi$를 실어 나르는 유속이다. 물리적 속도 $\mathbf{u}$ 옆에 반확산이 속도처럼 나란히 앉아 있다.

Chapman–Enskog 전개를 돌리면 확산 계수가 완화시간에서 나온다.

두 개를 합치면 앞의 보존형 Allen–Cahn 방정식이 그대로 복원된다. 완화가 $+M\nabla^2\phi$를 만들고, 평형분포에 얹은 1차 모멘트가 $-\nabla\cdot(M\,4\phi(1-\phi)\mathbf{n}/W)$를 만든다. 표의 $M=0.10$을 $c_s^2=1/3$으로 나누면 $(\tau_\phi - 1/2)\delta t = 0.30$이다. $M=0.02$면 0.06, $M=0.20$이면 0.60이다. 이동도를 열 배 키운다는 것은 이 곱을 벌린다는 뜻이지, $W$를 건드린다는 뜻이 아니다.

구현 입장에서 이것이 주는 이득은 명확하다. streaming과 collision 구조를 건드리지 않는다. 계면을 잡는 일 전체가 평형분포 한 줄에 들어간다. 운동량 쪽 표면장력은 별개의 손잡이다. 두께 손잡이와 섞으면 안 된다.

대신 $\mathbf{n}$을 매 노드에서 구해야 한다#

공짜는 없다. 평형분포에 $\mathbf{n}=\nabla\phi/|\nabla\phi|$가 들어 있으니 $\nabla\phi$를 모든 격자점에서 계산해야 한다. 여기서 스텐실 선택이 결과를 가른다.

가장 단순한 선택은 축 방향 이웃 넷만 쓰는 중심차분이다. 격자 볼츠만에는 다른 선택지가 있다. 격자 속도 집합 전체에 가중치를 얹는 것이다.

D2Q9에서 $w_i$는 축 방향 $1/9$, 대각 방향 $1/36$이다. 대각선 이웃이 들어오면서 4차 이방성 오차가 상쇄된다. 같은 액적 위에서 두 스텐실의 오차는 아래 출력에 숫자로 남아 있다.

평균 오차가 더 큰 쪽을 골라야 하는 이유#

반지름 24, 두께 $W$인 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로 바꾸면 오른쪽 장미가 네 잎에서 원으로 접힌다.

$W$를 4에 두면 이방성 2.575%와 0.155%가 그 차이다. 2 쪽으로 내리면 둘 다 나빠지지만, 찌그러지는 쪽은 중심차분뿐이다.

$W=4$에서 D2Q9의 평균 오차는 6.803%로 중심차분의 5.338%보다 크다. 정확도만 보면 지는 쪽이다. 그런데 각도에 따른 편차는 2.575%에서 0.155%로, 17배 가까이 줄었다. $2.575/0.155 = 16.6$이다. 평균 오차 차이는 $6.803 - 5.338 = 1.465$, 곧 1.5%포인트다. $W=5$에서 이방성은 1.812에서 0.075로 줄고, $W=3$에서는 3.749에서 0.402로 준다. 평균 오차는 같은 방향으로 커진다. 중심차분 3.641, 5.338, 8.369. D2Q9는 4.675, 6.803, 10.539다. 얇게 만들수록 등방성 이득은 남고, 크기 오차는 양쪽 다 $W=5$ 대비 두 배 넘게 커진다.

이 둘은 값이 다른 것이 아니라 성격이 다르다. 계면 전체에 균일하게 얹힌 오차는 $\nabla\phi$의 크기를 일정 비율로 틀리게 만든다. 표면장력 계수를 그만큼 잘못 잡은 것과 같고, 정지한 액적은 여전히 정지해 있다. 각도에 따라 달라지는 오차는 다르다. 같은 원 위의 지점마다 다른 크기의 힘이 걸린다. 그 불균형이 정지해 있어야 할 액적 표면에 소용돌이를 만든다. 흔히 말하는 기생 전류(spurious current)다.

그래서 여기서는 평균 오차 1.5%포인트를 내주고 이방성 17배를 산다. 남는 장사다.

$W$를 정할 때 실제로 흥정하는 것#

표의 세 번째 열과 네 번째 열이 반대 방향으로 움직인다. $W$를 5에서 3으로 줄이면 계면은 물리적으로 날카로워지지만 $\nabla\phi$의 평균 오차는 3.6%에서 8.4%로, 중심차분 이방성은 1.8%에서 3.7%로 커진다. $W$를 키우면 오차는 줄지만 계면이 뭉툭해져 작은 액적이나 얇은 막을 해상하지 못한다.

실무에서 $W=4\sim5$가 굳어진 이유가 이 표에 있다. 계면을 얇게 만들고 싶을 때 $M$을 줄이는 것은 아무 소용이 없다. 첫머리의 4.44와 4.24가 그 이야기다. $M$을 건드리면 회복 속도만 바뀌고, 너무 줄이면 이송이 만든 수치 확산을 반확산이 따라잡지 못한다. 두께 손잡이는 $W$ 하나다. 세 상의 장력을 서로 다르게 줄 때 계면이 새는 자리는 쌍별 장력이 어긋나는 3상 Allen–Cahn에 따로 있다.