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

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

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

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

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

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

temperature를 1 아래로 내리면 등온선이 접힌다. 한 압력에 밀도가 셋 붙는 구간이 열린다.

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.

가로선은 두 초록 로브의 넓이가 같아지는 맥스웰 압력이다. 파란 점과 분홍 점의 간격이, 힘 항이 나중에 버텨야 하는 밀도차다.

힘은 어디로 들어오는가#

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$ 를 그대로"), 뒤가 자유에너지형이다. 등온선이 접히는 온도는 앞 다이얼에서 이미 내려 봤다. 이제 두 힘이 연속에서 같은 이유만 장부에 적는다.

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

두 형태가 같은지는 한 관계식으로 결판난다. 등온에서 압력과 화학퍼텐셜을 잇는 깁스-듀엠(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}$ 수준이고, 그마저 유한차분으로 미분을 잰 탓이다. 온도를 내릴수록 공존 밀도의 간격이 벌어진다. 깁스-듀엠은 연속에서 정확히 성립한다. 여기까지는 압력형과 자유에너지형이 완전히 같다.

공존 밀도를 칸별로 빼 보자. $T/T_c=0.95$ 는 기체 $0.5790$ , 액체 $1.4617$ , 차 $0.8827$ 이다. $0.90$ 은 $1.6573-0.4257=1.2316$ 이고, $0.85$ 는 $1.8071-0.3197=1.4874$ 다. 온도를 $0.10$ 내리는 동안 밀도차가 $0.8827$ 에서 $1.4874$ 로 약 $1.7$ 배가 된다. 같은 구간에서 결함은 $2.94\times10^{-4}$ , $9.44\times10^{-4}$ , $1.98\times10^{-3}$ 으로 커진다. 항등식 $\mathrm{d}p=\rho\,\mathrm{d}\mu$ 가 깨진 게 아니다. 400점 np.gradient가 더 가파른 등온선에서 차분 오차를 더 낸 것이다. $T=0.90$ 의 액체/기체 비는 $1.6573/0.4257\approx 3.89$ 다. 이 비가 forcing이 버틸 밀도 점프의 크기다.

힘 두 줄을 겹치면 차이가 나는 자리는 하나다. $\mathbf{F}_p$ 의 $-\nabla(p-c_s^2\rho)$ 와 $\mathbf{F}_\mu$ 의 $-\rho\nabla\mu+c_s^2\nabla\rho$ 는 $\mathrm{d}p=\rho\,\mathrm{d}\mu$ 아래에서 같다. 코르테베그 항 $\kappa\rho\nabla(\nabla^2\rho)$ 는 양쪽에 같은 꼴이라 차이를 만들지 않는다. $c_s^2=1/3$ 은 격자 음속이 이미 넣어 둔 보정이지, 계면을 미는 물리 힘이 아니다. 연속 장부에서 두 형태가 같다는 문장은 여기서 끝난다. 격자에서 그 문장이 깨지는 크기는 뒤 표의 $1.60\times10^{-2}$ 다.

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

문제는 격자다. 연속의 $\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

resolution을 올리면 계면이 걸리는 셀 수가 늘고, 빨간 봉우리의 높이가 같이 죽는다.

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이다. 셀을 늘리면 봉우리가 주저앉는다. 유령 힘은 물리가 아니라 이 차분의 부산물이다.

자유에너지형은 계면에서 $2.02\times10^{-16}$ 이다. 압력형은 $1.60\times10^{-2}$ 를 남긴다. 두 값의 차와 이산 깁스-듀엠 결함이 소수점까지 같다. $1.60\times10^{-2}/2.02\times10^{-16}\approx 7.9\times10^{13}$ 이다. 기계 정밀도와 유령 힘의 자릿수가 한 줄에 있다.

격자 표도 같은 장부다. 계면당 10셀이면 결함 $1.140\times10^{-2}$ , 20셀이면 $1.428\times10^{-3}$ 이라 관측 차수 $3.00$ 이다. 40셀은 $5.115\times10^{-5}$ (차수 $4.80$ ), 80셀은 $1.611\times10^{-6}$ (차수 $4.99$ ). 10셀에서 80셀로 여덟 배만 조여도 $1.140\times10^{-2}/1.611\times10^{-6}\approx 7.1\times10^{3}$ 배 줄어든다. 중심차분이 $\mathrm{d}p=\rho\,\mathrm{d}\mu$ 를 계면 기울기에서 놓친 크기가, 두 힘의 차와 같다. 이 힘이 정지한 계면을 민다.

압력형을 남기면 얇은 계면이 값을 청구한다#

자유에너지형( $-\rho\nabla\mu$ )은 정의상 계면에서 균형이 맞는다. 거친 격자에서도 유령 힘이 $2.02\times10^{-16}$ 에 붙는다. 압력을 그대로 쓰려면 $\nabla p$ 를 마음대로 차분화하지 말고 $\rho\,\nabla\mu$ 와 같은 차분으로 짠다. 그래야 이산 깁스-듀엠이 살고, 표의 $1.60\times10^{-2}$ 가 내려간다.

4~5셀 이상으로 계면을 풀면 압력형 결함은 차수 $4.99$ 까지 떨어진다. 실무 계면은 3셀 안팎이다. 그 두께에서 유령 힘은 밀도차의 제곱을 따라 남는다. 힘 항이 연속 등식만 만족하고 이산 모멘트까지 비면 같은 잔차가 남는데, 그 장부는 힘 텀의 3차 모멘트 한계에 있다. " $p$ 를 그대로"의 대가는 계면 두께로 청구된다.

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