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

균일 압력은 Γ를 질량으로 나눠 되찾을 때 깨지고, Γ 자체를 업윈드할 때 남는다.

방정식을 하나 지웠더니 코드가 더 안 터졌다. 압축성 다상유동(서로 다른 유체가 섞이는 압축성 흐름)에서는 보통 반대로 배운다. 상(phase, 물질의 서로 다른 형태)마다 압력·속도·온도를 따로 푸는 7-방정식이 가장 안전하고, 식을 줄이면 위태롭다고. Collis(2025)는 그 순서를 뒤집는다. 식이 가장 적은 4-방정식이 오히려 버틴다. 열쇠는 계면 평형이다.

Naive 토글을 켜고 right gas γ를 1.66이나 1.09로 밀면, 계면의 빨간 압력선이 올라가고 초록선은 평평하게 남는다.

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.

빨간 선의 높이는 γ 점프와 같이 커진다. 스텝을 900까지 올려도 초록선이 바닥에 붙어 있는지가, 이 조건이 살았는지 죽은 것인지다.

평형의 사다리 — 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$ 로 나눈다. 두 결과는 계면에서 갈린다. 균일 압력인데도 비보존 경로가 어긋나면 압력이 흔들린다는 자리가 Abgrall이 남긴 기준이다. 그 차이가 압력 진동이다.

앞에서 γ를 밀 때 빨간 선이 커진 것이 이 차이다. 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$ 를 숫자로 펼치면 점프의 키가 보인다. $\Gamma = 1/(\gamma-1)$ 이라 공기 $\gamma=1.4$ 는 $\Gamma=2.5$ 다. 헬륨 $\gamma=1.66$ 은 $\Gamma=1/0.66 \approx 1.515$ 다. SF6 $\gamma=1.09$ 는 $\Gamma=1/0.09 \approx 11.11$ 이다. 공기→헬륨 차이는 약 0.985다. 공기→SF6는 약 8.61이다. SF6 점프가 헬륨 점프의 약 8.74배다. γ를 1.09로 밀 때 빨간 선이 커지는 배율이다.

날카로운 한 칸을 코드의 $\nu=0.6$ 으로 한 번만 민 경우다. 오른쪽 셀은 $\rho=0.2$ , $\Gamma\approx 1.515$ 이고, 업윈드 이웃은 공기 $\rho=1$ , $\Gamma=2.5$ 다. $\rho\Gamma$ 를 보존형으로 옮긴 뒤 $\rho$ 로 나누면 $\Gamma_{\mathrm{naive}}\approx 2.384$ 다. $\Gamma$ 를 직접 밀면 $\Gamma_{\mathrm{cons}}\approx 2.106$ 이다. 위의 에너지 식을 압력으로 되돌리면 $P/P_0 = \Gamma_{\mathrm{cons}}/\Gamma_{\mathrm{naive}} \approx 0.883$ 이라, 한 스텝의 $|P-P_0|\approx 0.117$ 이다. 같은 한 칸을 SF6로 바꾸면 $|P-P_0|\approx 0.692$ 다. 출력의 $3.83\times 10^{-2}$ 는 이 날카로운 칸이 아니다. 램프 폭 $0.06$ 은 $\Delta x=1/200$ 에서 12셀이고, 300스텝 동안 계면은 $\nu \cdot 300 \cdot \Delta x = 0.9$ 만큼 이동한다. 번진 계면을 오래 밀면 최대 오차는 한 스텝의 점프보다 작아져 $0.0383$ 에 머문다. 날카로운 한 스텝의 $0.117$ 과 이 $0.0383$ 의 비는 약 $3.05$ 다. 12셀로 나뉜 뒤에도 오차는 한 자릿수 퍼센트에 남는다. consistent의 $4.44\times 10^{-15}$ 는 경로 차가 0이라 반올림만 남은 것이다.

논문이 실제로 한 일#

위 토이는 이상기체 $\gamma$ 하나만 밀었다. Collis(2025)의 무대는 더 넓다.

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

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

Γ를 ρ로 나눠 되찾으면 평평한 압력이 한 스텝에 꺾인다#

한 스텝의 날카로운 칸에서 헬륨은 이미 $|P-P_0|\approx 0.117$ 이다. SF6는 $0.692$ 다. 300스텝 출력의 $0.0383$ 은 12셀로 번진 램프가 그 점프를 나눠 가진 결과다. $\rho\Gamma$ 를 보존형으로 옮긴 뒤 $\rho$ 로 나누는 한 줄만 남겨도, 균일했던 $P_0$ 는 계면에서 바로 꺾인다. 식을 줄인 대가로 이 한 줄을 비보존으로 바꾸는 선택이 4-방정식의 전부다.

5-방정식은 속도·압력을 공유하면서 여분 식을 하나 남긴다. 그 여분이 계면의 음속을 무너뜨리는 쪽이라면, IEC를 맞추려고 식을 지운 것이 그 병도 지운 것인지가 남는다. 물과 공기를 반씩 섞었을 때 음속이 무너지는 5식이 4-방정식의 혼합물 상태식 안에 숨어 있는지, 식을 지울 때 같이 사라지는지가 다음 갈림이다.