방정식을 줄였더니 더 안 터졌다 — 4-방정식 모델과 계면 평형 조건

방정식을 하나 지웠더니 코드가 더 안 터졌다. 압축성 다상유동(서로 다른 유체가 섞이는 압축성 흐름)에서는 보통 반대로 배운다. 상(phase·물질의 서로 다른 형태)마다 압력·속도·온도를 따로 풀어주는 7-방정식 모델이 가장 안전하고, 방정식을 줄일수록 위태롭다고. Collis(2025)는 그 통념을 뒤집는다. 가장 방정식이 적은 4-방정식 모델을 오히려 가장 견고하게 만들었다. 오늘은 그 열쇠인 계면 평형 조건을 따라간다. 왜 '보존형' 스킴이 균일한 압력을 흔드는지 Python으로 직접 확인한다.

평형의 사다리 — 7·6·5·4 방정식#

확산 계면 기법(diffuse interface·경계면을 몇 셀에 걸쳐 번지게 푸는 방법)은 번진 셀 안에서 두 상의 상태를 정의할 규칙이 필요하다. 그 규칙이 곧 "무엇을 평형으로 볼 것인가"다.

7-방정식(Baer–Nunziato): 상마다 속도·압력·온도가 다 다르다. 가장 일반적이고 가장 비싸다.
6-방정식: 속도는 공유. 압력은 따로 두고 이완(relaxation)으로 맞춘다.
5-방정식: 속도·압력을 공유. 온도만 따로.
4-방정식: 속도·압력·온도를 모두 공유. 가장 간결하다.

물리적 계면의 두께는 나노미터 규모다. 수치 계면은 그걸 격자 몇 칸으로 부풀린 것뿐이다. 그 좁은 띠 안에서 열·역학 평형은 사실상 즉각 이뤄진다. 그러니 계면 근처 물리에는 4-방정식이 가장 가깝다.

T^{(1)} = T^{(2)}, \quad P^{(1)} = P^{(2)}, \quad \mathbf{u}^{(1)} = \mathbf{u}^{(2)}

위 첨자 $(1),(2)$ 는 두 상. 4-방정식은 이 세 등식을 셀 안에서 모두 강제한다.

방정식이 적을수록 왜 견고한가#

5-방정식 이상은 대개 부피분율(volume fraction·셀 안 한 상의 부피 비율)이나 비열비 관련 $1/(\gamma-1)$ 같은 여분의 이송 방정식을 더 푼다. PDE 수준에서 이들은 중복이다. 중복 방정식은 이산화하면 서로 어긋난다. 같은 양을 두 경로로 계산해 값이 갈리는 것이다.

4-방정식은 그 여분을 없앤다. 질량·운동량·에너지만 보존한다. 압력은 혼합물 상태방정식으로 닫는다. 방정식이 적으니 어긋날 자리도 적다. 비용도 준다.

문제는 하나뿐이다. 여분의 방정식이 주던 여유가 사라지면서, 보존형 이산화가 계면에서 압력을 어긋나게 만든다.

계면 평형 조건이 깨지는 순간#

계면 평형 조건(IEC·interface equilibrium condition)의 정의는 단순하다. 압력·온도·속도가 처음에 균일하면, 계면이 지나가도 그 균일함이 깨지면 안 된다.

이를 어긴 스킴은 계면 주위에 없던 압력 진동을 만든다. 진동은 시간이 갈수록 자란다. 밀도비가 크면 음의 밀도나 음의 내부에너지가 튀어나와 코드가 멈춘다.

왜 보존형이 이걸 어길까. 이상기체 혼합을 예로 들자. $\Gamma \equiv 1/(\gamma-1)$ 로 두면 총에너지는 이렇게 쓰인다.

E = P\,\Gamma + \tfrac{1}{2}\,\rho\,u^2

$E$ 총에너지, $P$ 압력, $\rho$ 밀도, $u$ 속도. $\Gamma$ 는 유체 종류에 따라 계면에서 점프한다.

균일한 $u,P$ 아래 1차 업윈드로 한 칸 밀어보자. $\nu = u\,\Delta t/\Delta x$ 는 CFL 수(한 스텝에 정보가 셀을 지나는 비율). 보존형 에너지 갱신은 $+P$ 항이 상쇄돼 깔끔하게 준다.

E_i^{n+1} = E_i - \nu\,(E_i - E_{i-1})

여기에 위 총에너지 정의를 대입하면 오른쪽이 이렇게 갈라진다.

E_i^{n+1} = P_0\big[\,\Gamma_i - \nu(\Gamma_i - \Gamma_{i-1})\,\big] + \tfrac{1}{2}\,u_0^2\,\rho_i^{n+1}

새 압력이 $P_0$ 로 남으려면 $\Gamma$ 가 딱 저 대괄호처럼 갱신돼야 한다. 즉 비보존 업윈드로. 그런데 보존형은 $\rho\Gamma$ 를 통째로 이송하고 다시 $\rho$ 로 나눈다. 두 결과는 계면에서 갈린다. 그 차이가 곧 압력 진동이다.

아래 시뮬레이션에서 직접 조작해보자.

right gas γ = 1.66CFL = 0.60

Uniform velocity, uniform pressure, only γ jumps across the shaded interface. The naive scheme grows a pressure spike at the interface; the consistent one keeps P flat to machine precision.

Naive를 켜면 계면에서 빨간 압력선이 튀어 오른다. right gas $\gamma$ 를 1.66(헬륨)이나 1.09(SF6)로 밀어 점프를 키우면 진동도 커진다. IEC-consistent로 바꾸면 초록선은 스텝이 900까지 가도 평평하다.

Python — 진동의 정체와 처방#

같은 계산을 코드로 옮긴다. 균일한 $u,P$ , $\gamma$ 만 점프하는 계면을 1차 업윈드로 민다. 두 스킴만 다르다.

import numpy as np
 
def advect_interface(scheme, gamma_R=1.66, nu=0.6, steps=300, N=200):
    x = np.linspace(0.0, 1.0, N, endpoint=False)
    ramp = np.clip((x - 0.28) / 0.06, 0.0, 1.0)   # 몇 셀에 번진 계면
    gamma = 1.4 + (gamma_R - 1.4) * ramp          # 공기(1.4) -> 선택 기체
    rho = 1.0 + (0.2 - 1.0) * ramp                # 밀도도 점프한다
    u0, P0 = 1.0, 1.0
    Gam = 1.0 / (gamma - 1.0)                      # 혼합물 강성 1/(gamma-1)
    mom = rho * u0
    E = P0 * Gam + 0.5 * rho * u0**2               # E = P*Gamma + 0.5*rho*u^2
 
    for _ in range(steps):
        rhoG = rho * Gam                           # 갱신 전 rho*Gamma
        r_up, m_up, e_up = np.roll(rho, 1), np.roll(mom, 1), np.roll(E, 1)
        rho = rho - nu * (rho - r_up)              # 질량: 보존형 업윈드
        mom = mom - nu * (mom - m_up)              # 운동량: 보존형
        E   = E   - nu * (E   - e_up)              # 에너지: 보존형
        if scheme == "naive":
            rhoG = rhoG - nu * (rhoG - np.roll(rhoG, 1))  # rho*Gamma 보존 이송
            Gam = rhoG / rho                              # 다시 rho로 나눈다
        else:  # consistent
            Gam = Gam - nu * (Gam - np.roll(Gam, 1))      # Gamma 직접 비보존 업윈드
 
    P = (E - 0.5 * mom**2 / rho) / Gam             # 압력 복원
    return float(np.max(np.abs(P - P0)))
 
for s in ("naive", "consistent"):
    print(f"{s:11s}  max|P-P0| = {advect_interface(s):.3e}")

출력은 두 자릿수 차이를 그대로 드러낸다.

naive        max|P-P0| = 3.83e-02
consistent   max|P-P0| = 4.44e-15

naive는 계면에서 압력이 약 4% 어긋난다. consistent는 기계 정밀도(부동소수점 반올림 한계)까지 평평하다. 차이를 만든 건 단 한 줄. $\Gamma$ 를 rho로 나눠 되찾느냐, 아니면 처음부터 $\Gamma$ 자체를 이송하느냐다.

논문이 실제로 한 일#

위 토이는 이상기체 하나의 $\gamma$ 만 다뤘다. Collis(2025)의 무대는 훨씬 넓다.

NASG 상태방정식(Noble–Abel Stiffened Gas·액체까지 다루는 상태식)으로 물·공기·헬륨·SF6을 함께 닫는다. 혼합 압력은 Amagat 법칙으로 닫힌 형태의 이차식으로 풀린다.
1차 업윈드 대신 ENO 계열(WENO/TENO) 고차 재구성을 쓴다. 핵심은 그 재구성을 4-방정식의 열·역학 평형 가정과 일관되게 짜서 IEC를 기계 정밀도로 만족시킨 것이다.
강한 충격파가 큰 밀도비 계면을 때릴 때를 위해 양수성 보존 한정자(positivity-preserving limiter)를 4-방정식으로 확장했다. 밀도와 음속제곱이 음수로 새지 않게 막는다.
이 모두를 여분의 부피분율·비열비 방정식 없이, 보존성도 스칼라 필터도 희생하지 않고 이뤘다.

한계도 있다. 액체 상은 단일 성분, 기체 상은 이상기체 성분으로 제한했다. 성분이 여럿이면 압력·온도 평형을 반복 해석기로 풀어야 한다. 아직 프리프린트(동료심사 전)라는 점도 기억할 것.

기억할 점#

4-방정식 모델은 속도·압력·온도를 상 사이에 모두 공유한다. 방정식이 가장 적고 물리에 가장 가깝다.
계면 평형 조건(IEC): 균일한 압력·속도가 계면을 지나도 흔들리지 않아야 한다. 보존형 스킴은 혼합물 EOS 파라미터를 잘못 이송해 이 조건을 깬다.
처방은 그 파라미터( $1/(\gamma-1)$ )를 평형 가정과 일관되게 이송하는 것. Collis(2025)는 이를 고차 ENO로 끌어올려 여분 방정식 없이 4-방정식을 견고하게 만들었다.