존재하지 않는 유체로 경계면을 메운다 — Ghost Fluid Method

물과 공기가 만나는 한 칸에, 실제로는 거기 없는 가짜 공기를 채워 넣는다. 이 유령 유체(ghost fluid, 경계 너머로 외삽한 가상의 유체)는 격자 어디에도 존재하지 않는다. 그런데 바로 그 유령이 충격파가 경계면을 지날 때 계산을 안정시킨다. 복사할 양과 외삽할 양이 갈린다. 오리지널 GFM이 떠는 자리도 그 갈림에서 나온다.

밀도 점프를 키우면 계단은 가팔라지고, 압력은 빨간 점선을 그대로 지난다. "show ghost"는 경계 너머 연장선을 켰다 끈다.

점선을 켜 두면 각 유체가 자기 엔트로피를 들고 반대편으로 넘어간다. 경계면 위치를 옮겨도 압력 곡선은 한 줄로 이어진다. 끊기는 것은 밀도다.

한 칸에 두 물질이 부딪힐 때#

기체-액체, 폭약-금속처럼 물성이 크게 다른 두 물질이 한 격자 위에서 만난다. 단순히 평균 물성을 쓰면 경계면이 수치적으로 번진다(smearing). 밀도비가 1000배인 물-공기 경계에서는 이 번짐이 곧 비물리적 진동으로 번진다.

GFM의 발상은 다르다. 경계면을 흐리게 섞지 않는다. 각 유체를 자기 영역 밖으로 가짜로 연장한 뒤, 한쪽 물질만 있는 것처럼 단일물질 솔버를 두 번 돌린다. 경계면은 레벨셋(level set, 부호로 영역을 나누는 거리함수) $\phi$의 영점으로 추적한다.

여기서 $\phi$는 부호 거리함수, 영점 $\phi=0$이 경계면이다.

경계면에서 끊기는 것과 이어지는 것#

경계면을 가로질러 무엇이 연속인지 정확히 알아야 유령을 옳게 채운다. 압력 $p$와 수직 속도 $u_n$은 경계면에서 연속이다. 반면 엔트로피 $s$와 접선 속도 $u_t$는 불연속이다.

대괄호 $[\cdot]$는 경계면을 가로지른 점프를 뜻한다. 이 네 줄이 GFM의 전부다. 연속량은 복사하고, 불연속량은 외삽한다. 앞의 그림에서 압력만 매끄럽고 밀도만 계단인 이유가 여기 있다.

경계를 한 셀 안에서 평균해 버리면 다른 고장이 난다. 6방정식 two-fluid 모델은 밀도비가 큰 경계에서 고유값이 복소로 넘어가고, 격자를 줄이면 발산이 더 빨라진다. GFM은 그 평균을 하지 않는다. 칸은 언제나 한쪽 물질만 품는다.

물을 기체처럼 — Stiffened EOS#

물은 거의 비압축성이라 이상기체 상태방정식으로는 다룰 수 없다. GFM은 두 물질을 같은 형식의 Mie–Grüneisen 계열로 묶기 위해 Stiffened Gas EOS를 쓴다.

$\rho$는 밀도, $e$는 단위질량당 내부에너지, $\gamma$는 비열비, $p_\infty$는 물질 강성을 나타내는 상수다. 이상기체는 $p_\infty = 0$, 물은 대략 $\gamma \approx 4.4$, $p_\infty \approx 6 \times 10^8\,\mathrm{Pa}$로 맞춘다. 음속은 다음과 같다.

$p_\infty$ 덕분에 같은 압력에서도 물의 음속이 공기보다 훨씬 크게 나온다. 한 줄로 EOS를 통일하면 두 유체를 같은 솔버에 넣을 수 있다.

유령 셀을 채우는 두 가지 손길#

이제 핵심 절차다. 유체 A의 영역을 풀 때, B가 차지한 칸을 A의 유령으로 메운다. 손길은 둘로 나뉜다.

1. 복사: 연속량인 압력과 수직 속도는 경계면 너머 실제 B 셀의 값을 그대로 가져온다.
2. 외삽: 불연속량인 엔트로피(즉 밀도)와 접선 속도는 가장 가까운 실제 A 셀에서 경계면을 넘어 등엔트로피로 외삽한다.

등엔트로피 외삽은 엔트로피 $s = p / \rho^\gamma$를 보존하며 밀도를 다시 푼다.

$s_A$는 실제 A 유체의 엔트로피, $p_\text{copy}$는 복사해 온 압력이다. 외삽은 보통 경계면 법선 방향으로 빠르게 퍼뜨리는 fast marching 방법으로 수행한다. 한 번의 스윕으로 A의 유령장이 완성되면, B의 존재를 잊고 A만의 단일물질 Euler 방정식을 푼다. B에 대해서도 같은 절차를 반복한다.

코드에 적힌 초기 조건만으로 유령 밀도를 계산하면, 외삽이 실제 밀도를 어떻게 바꾸는지가 보인다. 왼쪽은 밀도 1, 압력 1, 비열비 1.4다. 엔트로피 $s = p / \rho^\gamma$는 1이다. 오른쪽 압력 0.1을 복사해 왼쪽 유체를 연장하면 유령 밀도는 0.1의 $1/1.4$제곱이라 약 0.193이다. 실제 오른쪽 밀도 0.125가 아니다. 단일물질 솔버가 보는 밀도비는 8배가 아니라 약 5.2배다.

반대편도 같은 식이다. 오른쪽 비열비는 1.67이다. 엔트로피는 0.1을 $0.125^{1.67}$로 나눈 값, 약 3.22다. 왼쪽 압력 1을 이 엔트로피에 넣으면 유령 밀도는 약 0.496이다. 실제 왼쪽 밀도 1의 절반에 가깝다. 압력은 넘기고, 각 유체의 엔트로피는 그 유령에 남긴다.

음속도 상태방정식에서 바로 나온다. 강성 상수가 0이면 음속은 비열비, 압력, 밀도로 만든 제곱근이다. 왼쪽은 약 1.183, 오른쪽은 약 1.156이다. 점 400개에 길이 1이면 간격은 $1/399$, 약 0.00251이다. CFL 0.4를 왼쪽 음속에 대면 시간 간격은 약 0.00085다. 종료 시각 0.14를 나누면 165스텝 안팎이다. 음파가 길이 1을 건너는 시간은 약 0.85라, 0.14는 그 5분의 1이 안 된다. 물-공기의 밀도비 1000에 비하면 이 충격관의 8배는 온건하다. 그래도 유령 밀도는 0.193과 0.496으로 실제 값과 이미 어긋나 있다.

오리지널 GFM이 떨리는 자리#

오리지널 GFM은 경계면 한쪽이 매우 강성일 때(stiff, 물처럼 $p_\infty$가 클 때) 문제가 생긴다. 강한 충격파나 폭발파가 경계면을 지나면, 단순 복사·외삽으로 만든 유령 상태가 실제 Riemann 해와 어긋난다. 그 어긋남이 경계면 근처에서 압력 진동으로 나타나고, 진동은 종종 계산 발산으로 이어진다.

아래 시뮬레이션에서 압력 펄스를 강성 유체(분홍 영역)로 통과시켜 보자. stiffness ratio를 올릴수록 경계면을 지난 직후의 잔물결이 커진다.

stiffness ratio를 9까지 올리면 경계면 직후에 고주파 진동이 또렷하게 보인다. 이 진동은 격자를 세분화해도 사라지지 않는다. 원인이 격자가 아니라 경계조건이기 때문이다.

수정 GFM — 경계면 Riemann 문제로 값을 정하다#

해법은 유령 상태를 추측하지 말고 풀어서 구하는 것이다. 수정 GFM(Modified GFM, MGFM 또는 IGFM)은 경계면 양쪽 실제 상태를 좌·우로 삼아 Riemann 문제를 푼다.

$W_A, W_B$는 경계면 양쪽 실제 유체의 상태, $\mathcal{R}$은 Riemann 해법, $p^*, u^*$는 접촉 불연속의 압력·속도다. 이 $p^*, u^*$를 양쪽 유령에 똑같이 부여하면 압력·속도 연속 조건이 정확히 만족된다. 위 시뮬레이션에서 "use modified GFM"을 켜면 같은 강성에서도 진동이 사라진다.

Python — 기체-기체 GFM 충격관#

두 이상기체($\gamma_A = 1.4$, $\gamma_B = 1.67$)가 만나는 충격관을 등엔트로피 외삽 GFM으로 푼다. HLL 플럭스를 쓰고, 매 스텝 두 번의 유령 스윕 뒤 레벨셋으로 병합한다.

import numpy as np
 
GAMMA = (1.4, 1.67)   # 좌(공기), 우(단원자 기체)의 비열비
 
def primitive(U, g):
    rho = U[0]; u = U[1] / rho
    e = U[2] / rho - 0.5 * u * u
    p = (g - 1.0) * rho * e          # Stiffened EOS, p_inf = 0 (이상기체)
    return rho, u, p
 
def conserved(rho, u, p, g):
    E = p / (g - 1.0) + 0.5 * rho * u * u
    return np.array([rho, rho * u, E])
 
def hll(UL, UR, g):
    rL, uL, pL = primitive(UL, g); rR, uR, pR = primitive(UR, g)
    aL = np.sqrt(g * pL / rL); aR = np.sqrt(g * pR / rR)
    sL = min(uL - aL, uR - aR); sR = max(uL + aL, uR + aR)
    FL = np.array([rL * uL, rL * uL * uL + pL, uL * (UL[2] + pL)])
    FR = np.array([rR * uR, rR * uR * uR + pR, uR * (UR[2] + pR)])
    if sL >= 0: return FL
    if sR <= 0: return FR
    return (sR * FL - sL * FR + sL * sR * (UR - UL)) / (sR - sL)
 
def extrapolate_ghost(rho, u, p, mat, side):
    # side 유체로 전체 영역을 채운다: 압력·속도는 복사, 밀도는 등엔트로피 외삽
    rg, ug, pg = rho.copy(), u.copy(), p.copy()
    g = GAMMA[side]
    idx = np.where(mat == side)[0]
    lo, hi = idx[0], idx[-1]
    ref = lo if side == 1 else hi      # 경계면에 가장 가까운 실제 셀
    s_ref = p[ref] / rho[ref] ** g     # 보존할 엔트로피
    for i in range(len(rho)):
        if mat[i] != side:
            j = min(max(i, lo), hi)            # 가장 가까운 실제 셀
            ug[i] = u[j]; pg[i] = p[j]         # 복사
            rg[i] = (pg[i] / s_ref) ** (1.0 / g)   # 등엔트로피 외삽
    return rg, ug, pg
 
def run_ghost_fluid_tube(N=400, t_end=0.14, cfl=0.4):
    x = np.linspace(0, 1, N); dx = x[1] - x[0]
    x0 = 0.5                                   # 초기 경계면 위치
    mat = np.where(x < x0, 0, 1)
    rho = np.where(x < x0, 1.0, 0.125)
    u   = np.zeros(N)
    p   = np.where(x < x0, 1.0, 0.1)
    t = 0.0
    while t < t_end:
        a = np.sqrt(np.where(mat == 0, GAMMA[0], GAMMA[1]) * p / rho)
        dt = min(cfl * dx / np.max(np.abs(u) + a), t_end - t)
        new = {}
        for side in (0, 1):                    # 두 번의 유령 스윕
            rg, ug, pg = extrapolate_ghost(rho, u, p, mat, side)
            g = GAMMA[side]
            U = np.array([conserved(rg[i], ug[i], pg[i], g) for i in range(N)]).T
            F = np.zeros((3, N + 1))
            for i in range(1, N):
                F[:, i] = hll(U[:, i - 1], U[:, i], g)
            F[:, 0] = F[:, 1]; F[:, N] = F[:, N - 1]
            new[side] = U - dt / dx * (F[:, 1:] - F[:, :-1])
        U = np.where(mat[None, :] == 0, new[0], new[1])   # 레벨셋 병합
        for i in range(N):
            rho[i], u[i], p[i] = primitive(U[:, i], GAMMA[mat[i]])
        x0 += dt * u[np.argmin(np.abs(x - x0))]           # 경계면 이동
        mat = np.where(x < x0, 0, 1)
        t += dt
    return x, rho, u, p
 
if __name__ == "__main__":
    x, rho, u, p = run_ghost_fluid_tube()
    print(f"t=0.14:  max p = {p.max():.3f}, "
          f"접촉면 부근 밀도 {rho[180]:.3f} -> {rho[220]:.3f}")

출력은 경계면을 가로지른 밀도 점프가 살아 있고 압력은 매끄럽게 이어졌음을 보여준다. 두 번의 단일물질 스윕이 GFM의 뼈대다.

코드 짤 때 빠지지 말아야 할 함정#

외삽을 격자 축으로 흘리면 비스듬한 경계에서 엔트로피가 섞인다. 퍼뜨리는 방향은 경계면 법선뿐이다. 강성비가 큰 물-공기는 오리지널 복사와 외삽이 압력 잔물결을 거의 항상 만든다. 그 경우는 처음부터 경계면 Riemann 해로 유령을 채운다. 레벨셋은 부호 거리로 다시 맞춰 주지 않으면 영점은 남는데 거리 값이 흐트러지고, 경계면 속도가 어긋난다.

기체-기체에서 통한 외삽을 물에 그대로 쓸 때#

$p_\infty$가 0인 이상기체 둘($\gamma = 1.4$와 $1.67$)은 등엔트로피 외삽만으로 접촉면 밀도 계단이 남는다. 물은 $\gamma \approx 4.4$, $p_\infty \approx 6 \times 10^8\,\mathrm{Pa}$다. 같은 압력에서 음속 식의 분자에 이 상수가 더해지면, 데모에서 stiffness ratio를 9까지 올려 본 잔물결과는 다른 크기의 문제가 된다. 격자를 잘게 해도 그 잔물결은 경계조건 안에 남는다.

유령에 넣을 압력과 속도를 복사로 둘지, 경계면 Riemann 해의 $p^*, u^*$로 바꿀지는 밀도 점프를 세기 전에 정한다. 외삽 밀도가 음수로 무너지는 자리는 리미터 문제가 따로다. 양수 보존 flux 리미터가 질량은 맞는데 밀도가 음수가 되는 칸을 어떻게 막는지, 강성 경계 다음에 이어서 볼 만하다.