압력을 그대로 넣었더니 정지한 계면이 떨렸다

압력형과 자유에너지형 forcing은 깁스-듀엠으로 연속에서 같다. 하지만 격자 위에서는 압력형이 계면에 유령 힘을 남긴다.

정지한 방울이 스스로 흐르기 시작했다#

반데르발스 상태방정식을 격자 볼츠만(LBM, Lattice Boltzmann Method)에 얹어 2상 유체를 풀었다. 초기조건은 정지한 방울 하나. 밀도장은 매끈하고, 속도는 전부 0이다.

몇 스텝을 돌리자 계면 근처에서 작은 속도가 돋아났다. 아무도 밀지 않았는데 유체가 흐른다. 이런 인공 속도를 기생 전류(parasitic current, 물리적 원인 없이 계면에서 생기는 유령 흐름)라 부른다.

원인을 좁혀 들어가니 코드 한 곳에 닿았다. 힘 항에 무엇을 넣느냐다. 압력 $p$ 를 그대로 넣을 것인가, 화학퍼텐셜 $\mu$ 를 넣을 것인가. 교과서는 둘이 같다고 한다. 이 글은 그 "같다"가 어디까지 참인지를 장부로 적는다. 답은 한 줄이다 — 연속에서는 같고, 격자 위에서는 갈라진다.

힘은 어디로 들어오는가#

LBM은 분포함수를 이류·충돌시켜 거시 방정식을 복원한다. 이상기체라면 저절로 나오는 압력은 $c_s^2 \rho$ 하나뿐이다( $c_s$ 는 격자 음속). 반데르발스 같은 비이상 유체의 진짜 압력은 그것과 다르다. 그 차액을 메우는 것이 힘 항 $\mathbf{F}$ 의 역할이다.

목표는 다음 운동량 방정식을 되살리는 것이다.

\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u}\mathbf{u}) = -\nabla p + \nabla \cdot \boldsymbol{\tau} + \kappa \rho \, \nabla (\nabla^2 \rho)

$p$ 는 반데르발스 압력, $\boldsymbol{\tau}$ 는 점성 응력, 마지막 항은 계면을 세우는 코르테베그(Korteweg) 응력이며 $\kappa$ 는 그 세기다. 스트리밍이 $c_s^2\rho$ 를 주므로, 힘은 나머지 차액만 채우면 된다.

여기서 두 갈래가 생긴다. 압력을 그대로 쓰는 형태와, 자유에너지에서 온 화학퍼텐셜을 쓰는 형태다.

\mathbf{F}_p = -\nabla\!\left(p - c_s^2 \rho\right) + \kappa \rho \, \nabla(\nabla^2 \rho)

\mathbf{F}_\mu = -\rho \, \nabla \mu + c_s^2 \nabla \rho + \kappa \rho \, \nabla(\nabla^2 \rho)

앞이 압력형(오늘의 소재 " $p$ 를 그대로"), 뒤가 자유에너지형이다. 둘이 정말 같은지 보려면 먼저 반데르발스 상태방정식의 생김새를 알아야 한다. 아래에서 온도를 직접 내려 보자.

temperature0.900 Tc

Drag temperature down from 1.0: the isotherm folds into an S, and the amber tie line drops in so the two green lobes have equal area. The blue and pink dots are the vapor and liquid densities the LBM must hold apart — their gap is what the forcing term exists to sustain.

temperature를 1.0 아래로 내리면 등온선이 S자로 접힌다. 한 압력에 밀도가 셋 대응한다. 이때 물리적 2상 압력은 두 초록 로브의 넓이가 같아지는 자리(맥스웰 등면적 규칙)에서 정해진다. 파란 점과 분홍 점의 간격이 곧 forcing 항이 지탱해야 하는 밀도차다.

깁스-듀엠 — 압력과 화학퍼텐셜은 같은 것을 두 번 말한다#

두 형태가 같은지는 한 관계식으로 결판난다. 등온에서 압력과 화학퍼텐셜을 잇는 깁스-듀엠(Gibbs–Duhem) 관계다.

\mathrm{d}p = \rho \, \mathrm{d}\mu \qquad (T = \text{const})

이 식을 $\mathbf{F}_\mu$ 에 대입해 보자. $-\rho\nabla\mu = -\nabla p$ 이므로 화학퍼텐셜 항이 곧바로 압력 기울기로 바뀐다. 남는 $c_s^2\nabla\rho$ 는 압력형이 $-\nabla(-c_s^2\rho)$ 로 갖고 있던 항과 정확히 맞물린다. 결국 $\mathbf{F}_p = \mathbf{F}_\mu$ 다.

즉 " $p$ 를 그대로 넣어도 된다"는 주장의 유일한 근거는 깁스-듀엠이다. 반데르발스는 이 관계를 정확히 만족한다. 상태방정식과 자유에너지가 같은 열역학에서 나왔기 때문이다. Chapman–Enskog 전개로 확인한 LBM forcing 스킴들이 저마다 다른 얼굴을 해도 결국 같은 거시 방정식을 복원하는 것과 같은 종류의 등가다.

Python으로 확인한 공존 밀도와 깁스-듀엠#

말로만 두면 미덥지 않다. 반데르발스 상태방정식을 축약단위로 짜고, 온도별 공존 밀도를 뉴턴법으로 구한 뒤, 깁스-듀엠 결함 $\lvert \mathrm{d}p/\mathrm{d}\rho - \rho\,\mathrm{d}\mu/\mathrm{d}\rho \rvert$ 을 직접 재 보자.

import numpy as np
 
A, B, R = 9.0 / 8.0, 1.0 / 3.0, 1.0   # 반데르발스, 축약단위 (rho_c=1, T_c=1)
 
def p_eos(rho, T):                     # 반데르발스 압력
    return rho * R * T / (1.0 - B * rho) - A * rho * rho
 
def mu_eos(rho, T):                    # 화학퍼텐셜 mu = df/drho
    return R * T * (np.log(rho / (1.0 - B * rho)) + B * rho / (1.0 - B * rho)) - 2.0 * A * rho
 
def maxwell(T):                        # 등면적 규칙: p와 mu가 두 상에서 같아지는 밀도
    x = np.array([0.30, 1.90]); h = 1e-8
    for _ in range(80):
        f = np.array([p_eos(x[0], T) - p_eos(x[1], T), mu_eos(x[0], T) - mu_eos(x[1], T)])
        J = np.empty((2, 2))
        for k in range(2):
            y = x.copy(); y[k] += h
            g = np.array([p_eos(y[0], T) - p_eos(y[1], T), mu_eos(y[0], T) - mu_eos(y[1], T)])
            J[:, k] = (g - f) / h
        x -= np.linalg.solve(J, f)
    return x[0], x[1]
 
# 깁스-듀엠: dp = rho d(mu). 압력을 그대로 써도 되는 유일한 근거.
print("T/Tc  rho_vap  rho_liq |  max| dp/drho - rho*dmu/drho |")
for T in (0.95, 0.90, 0.85):
    rv, rl = maxwell(T)
    r = np.linspace(rv, rl, 400)
    dp = np.gradient(p_eos(r, T), r)
    dmu = np.gradient(mu_eos(r, T), r)
    err = np.abs(dp - r * dmu).max()
    print(f"{T:4.2f}  {rv:7.4f}  {rl:7.4f} |  {err:.2e}")

T/Tc  rho_vap  rho_liq |  max| dp/drho - rho*dmu/drho |
0.95   0.5790   1.4617 |  2.94e-04
0.90   0.4257   1.6573 |  9.44e-04
0.85   0.3197   1.8071 |  1.98e-03

결함은 $10^{-3}$ 수준이고, 그마저 유한차분으로 미분을 잰 탓이다. 온도를 내릴수록 공존 밀도의 간격이 벌어진다. 깁스-듀엠은 연속에서 정확히 성립한다. 여기까지는 압력형과 자유에너지형이 완전히 같다.

이산 격자에서 둘이 갈라진다#

문제는 격자다. 연속의 $\mathrm{d}p = \rho\,\mathrm{d}\mu$ 가 이산 미분에서도 성립하리란 보장은 없다. 중심차분한 $\nabla p$ 와 $\rho\,\nabla\mu$ 는 계면처럼 밀도가 급히 꺾이는 곳에서 서로 어긋난다.

정지한 평면계면을 오일러-라그랑주 조건으로 풀고, 그 위에서 두 힘 형태를 직접 계산했다.

import numpy as np
 
A, B, R = 9.0 / 8.0, 1.0 / 3.0, 1.0
CS2, KAPPA, T = 1.0 / 3.0, 0.02, 0.90
 
def p_eos(rho):  return rho * R * T / (1.0 - B * rho) - A * rho * rho
def mu_eos(rho): return R * T * (np.log(rho / (1.0 - B * rho)) + B * rho / (1.0 - B * rho)) - 2.0 * A * rho
def dmu(rho):    return R * T * (1.0 / (rho * (1.0 - B * rho)) + B / (1.0 - B * rho) ** 2) - 2.0 * A
def diff1(a, dx): return (np.roll(a, -1) - np.roll(a, 1)) / (2.0 * dx)
def lap(a, dx):   return (np.roll(a, -1) - 2.0 * a + np.roll(a, 1)) / dx ** 2
 
def maxwell():
    x = np.array([0.30, 1.90]); h = 1e-8
    for _ in range(80):
        f = np.array([p_eos(x[0]) - p_eos(x[1]), mu_eos(x[0]) - mu_eos(x[1])])
        J = np.empty((2, 2))
        for k in range(2):
            y = x.copy(); y[k] += h
            g = np.array([p_eos(y[0]) - p_eos(y[1]), mu_eos(y[0]) - mu_eos(y[1])])
            J[:, k] = (g - f) / h
        x -= np.linalg.solve(J, f)
    return x[0], x[1]
 
rv, rl = maxwell()
mu_co = mu_eos(np.array([rv]))[0]
 
# (1) 정지 평면계면을 풀고, 두 forcing 형태를 비교한다
NX = 240
xs = np.arange(NX)
rho = 0.5 * (rl + rv) + 0.5 * (rl - rv) * (np.tanh((xs - NX / 4) / 6.0) - np.tanh((xs - 3 * NX / 4) / 6.0) - 1.0)
for _ in range(6000):                                    # 오일러-라그랑주 잔차 완화
    rho -= 0.15 * (mu_eos(rho) - KAPPA * lap(rho, 1.0) - mu_co) / dmu(rho)
 
Gp  = -diff1(p_eos(rho) - CS2 * rho, 1.0) + KAPPA * rho * diff1(lap(rho, 1.0), 1.0) - diff1(CS2 * rho, 1.0)
Gmu = -rho * diff1(mu_eos(rho) - KAPPA * lap(rho, 1.0), 1.0) + CS2 * diff1(rho, 1.0) - diff1(CS2 * rho, 1.0)
gd  = np.abs(diff1(p_eos(rho), 1.0) - rho * diff1(mu_eos(rho), 1.0)).max()
 
print(f"coexistence         rho_vap = {rv:.4f}   rho_liq = {rl:.4f}")
print(f"free-energy form    max|G_mu|        = {np.abs(Gmu).max():.2e}   (well-balanced)")
print(f"pressure form       max|G_p|         = {np.abs(Gp).max():.2e}   (spurious force)")
print(f"gap between forms   max|G_p - G_mu|  = {np.abs(Gp - Gmu).max():.2e}")
print(f"discrete Gibbs-Duhem defect          = {gd:.2e}   <- the gap, exactly")
 
# (2) 같은 물리 계면을 dx만 조밀하게 하면 결함은 빠르게 0으로 수렴한다
print("\ncells/interface |  Gibbs-Duhem defect  order")
prev = None
for n in (10, 20, 40, 80):
    L = 40.0; N = int(L * n / 10)
    z = np.linspace(-L / 2, L / 2, N, endpoint=False); dx = z[1] - z[0]
    r = 0.5 * (rl + rv) - 0.5 * (rl - rv) * np.tanh(z / (0.1 * n))
    d = np.abs(diff1(p_eos(r), dx) - r * diff1(mu_eos(r), dx))[N // 4:3 * N // 4].max()
    order = "" if prev is None else f"{np.log(prev / d) / np.log(2.0):5.2f}"
    print(f"{n:9d}       |  {d:.3e}          {order}")
    prev = d

coexistence         rho_vap = 0.4257   rho_liq = 1.6573
free-energy form    max|G_mu|        = 2.02e-16   (well-balanced)
pressure form       max|G_p|         = 1.60e-02   (spurious force)
gap between forms   max|G_p - G_mu|  = 1.60e-02
discrete Gibbs-Duhem defect          = 1.60e-02   <- the gap, exactly
 
cells/interface |  Gibbs-Duhem defect  order
       10       |  1.140e-02          
       20       |  1.428e-03           3.00
       40       |  5.115e-05           4.80
       80       |  1.611e-06           4.99

세 줄이 핵심이다. 자유에너지형은 계면에서 힘이 $10^{-16}$ , 사실상 0이다. 압력형은 $1.6\times10^{-2}$ 의 힘을 남긴다. 그리고 두 형태의 차이가 이산 깁스-듀엠 결함과 소수점까지 정확히 같다. 압력형이 흘린 유령 힘의 정체가 바로 이 결함이다. 이 힘이 정지한 계면을 밀어 기생 전류를 만든다.

아래에서 계면을 얼마나 넓게 퍼뜨렸는지(격자 해상도)를 바꿔 보자.

resolution (cells / interface)2.5

The red hump is the force the pressure form adds over the free-energy form — the discrete Gibbs–Duhem defect. It lives exactly on the interface and pushes a fluid that should be at rest. Drag resolution up: the hump collapses onto the green zero line. It was never physics, only the price of writing p on too few cells.

빨간 봉우리는 압력형이 자유에너지형 위에 더 얹는 힘, 곧 $G_p - G_\mu$ (이산 깁스-듀엠 결함)다. 초록 점선은 자유에너지형의 0 기준선이다. resolution을 올려 계면을 더 많은 셀에 걸치게 하면 빨간 봉우리가 0으로 주저앉는다. 유령 힘은 물리가 아니라 이산화의 부산물이었다.

그래서 무엇을 넣을 것인가#

정리하면 선택은 두 가지다. 첫째, 자유에너지형( $-\rho\nabla\mu$ )을 쓴다. 이 형태는 정의상 계면에서 균형이 맞아, 거친 격자에서도 유령 힘이 없다. 둘째, 정말 압력을 그대로 쓰고 싶다면 $\nabla p$ 를 자유롭게 이산화하지 말고 $\rho\,\nabla\mu$ 와 일치하도록 짠다. 그러면 이산에서도 깁스-듀엠이 성립해 균형이 산다.

계면을 4~5셀 이상으로 충분히 풀 수 있다면 압력형의 오차는 위 표처럼 빠르게 사라진다. 하지만 실무의 계면은 대개 3셀 안팎으로 얇다. 그 영역에서 압력형은 밀도차의 제곱에 비례하는 유령 힘을 낸다. 같은 증상을 기생 전류와 well-balanced 계면 장력에서 표면장력 쪽으로 본 적이 있다. 뿌리는 하나다 — 연속에서 균형인 항을 이산에서도 균형으로 옮겼는가.

" $p$ 를 그대로 쓴다"는 편의는 공짜가 아니다. 그 대가는 계면 두께라는 화폐로 청구된다.

압력을 그대로 넣었더니 정지한 계면이 떨렸다 — 비이상 LBM forcing의 두 형태