음수 밀도가 새어나오는 자리 — 폭발파를 위한 Positivity-Preserving 플럭스 한정자

오전 2시 27분. Le Blanc 충격관 베리피케이션이 0.003초 만에 NaN을 토했다. CFL을 0.1까지 낮춰도, 격자를 두 배 더 박아도, 그래프 위쪽 두 곡선이 같이 무너졌다. 디버거를 띄우자 원인은 단순했다 — 강한 팽창파(rarefaction, 압력이 빠르게 떨어지는 영역)의 머리 쪽에서 한 셀의 밀도가 음수로 떨어진 것이다. 두 스텝 뒤 음속이 NaN이 되고, 다음 스텝에서 모든 셀이 같이 무너진다. 오늘은 그 음수 밀도가 어떤 자리에서 새어나오는지, 그리고 Hu·Adams·Shu(2013)가 제안한 "한 칸씩만 물러나는" θ-블렌드 플럭스 한정자가 그것을 어떻게 막는지를 따라간다. 끝에는 NumPy 50줄로 그 트릭을 직접 켜고 끄며 비교한다.

음수 밀도는 어디에서 새어나오는가#

고차 스킴은 안정성과 정확도 사이에서 위태롭게 균형을 잡는다. 1차 업윈드는 충격을 다섯 칸쯤 흐려놓는다. 2차 중심차분(central difference)은 충격 양쪽에 톱니 진동을 만든다. 진동이 가벼우면 압력만 살짝 흔들리고 끝난다. 두 셀 사이 밀도비가 100:1을 넘으면 — Le Blanc 충격관, Sedov 폭발, 초음속 노즐 후방 — 이야기가 달라진다.

밀도 보존식을 한 셀 $i$ 에 대해 정리한다.

\rho_i^{n+1} = \rho_i^n - \frac{\Delta t}{\Delta x}\bigl( F_{i+1/2}^\rho - F_{i-1/2}^\rho \bigr)

$F^\rho = \rho u$ 는 밀도 플럭스. $F_{i+1/2}^\rho \gg F_{i-1/2}^\rho$ 인 셀(한쪽으로 질량이 대량 빠져나가는 셀)에서 $\rho_i^{n+1}$ 은 음수가 될 수 있다. 1D Euler 시뮬에서 충격 직후 팽창 영역의 첫 셀이 정확히 이 상황이다.

음수가 되면 어떤 일이 생기는가. 음속 $a = \sqrt{\gamma p / \rho}$ 에서 $\rho \to 0^-$ 면 $a^2 < 0$ . 다음 플럭스 계산에서 $\sqrt{\cdot}$ 가 NaN을 반환한다. 한 셀의 NaN은 인접 두 셀로 번지고, 두 스텝 뒤에는 도메인 전체가 죽는다.

Lax–Friedrichs — 흐리지만 무너지지 않는다#

1954년 Lax가 제안한 LF 플럭스는 이 위험에서 자유롭다.

F_{i+1/2}^{LF} = \tfrac{1}{2}\bigl(F_L + F_R\bigr) - \tfrac{1}{2}\alpha\bigl(U_R - U_L\bigr)

$\alpha$ 는 좌·우 셀의 최대 신호 속도, $U = (\rho, \rho u, E)^T$ 는 보존 변수다. 둘째 항은 인공점성(artificial viscosity) — 이웃 셀끼리 평균화시켜 톱니 진동을 죽인다. Lax는 CFL $\leq$ 1 조건 아래 $\rho^{n+1} \geq 0$ 임을 증명했다 (Hu et al. 정리 1의 단순화 형태). 충격이 빠르건 약하건 LF는 무너지지 않는다.

대가는 정확도다. 두 칸짜리 충격은 다섯 칸이 되고, 접촉 불연속(contact discontinuity)은 한 시간 스텝마다 한 셀씩 더 흐려진다. 실무에서 LF만 쓰는 사람은 없다. 그래도 LF가 가진 안전 마진은 한정자의 출발점이 된다.

θ-블렌드 — 한 칸씩만 물러난다#

고차 플럭스 $F_{i+1/2}^{HO}$ (예: WENO, MUSCL, 중심차분)가 LF보다 정확하지만 음수 밀도를 만들 수 있다고 하자. Hu·Adams·Shu의 아이디어는 단순하다. 두 플럭스를 한 줄로 묶는다.

F_{i+1/2} = F_{i+1/2}^{LF} + \theta_{i+1/2}\bigl( F_{i+1/2}^{HO} - F_{i+1/2}^{LF} \bigr)

$\theta \in [0, 1]$ 은 면(face)마다 따로 정해지는 블렌딩 가중치다. $\theta=1$ 이면 순수 고차, $\theta=0$ 이면 순수 LF. 한정자의 역할은 단 하나 — 모든 셀이 다음 스텝에서 양수 밀도를 유지하도록 가장 큰 $\theta$ 를 고른다. 너무 멀리 가지 말고, 한 칸씩만 물러난다.

LF로 미리 계산한 결과 $\rho^{LF}_i \geq \rho_{\min}$ 임을 알고 있으므로, 보정 플럭스 $\Delta F = F^{HO} - F^{LF}$ 가 더해진 뒤에도 양수가 유지되는 $\theta$ 의 상한을 면별로 푼다. 한 면은 양쪽 두 셀에 영향을 주고, 두 셀 모두 floor 위에 있어야 하므로 더 작은 쪽이 한계다.

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

Per-face θ search · single rarefaction face · curve ρ_iⁿ⁺¹(θ) at fixed left face

log₁₀(ρ_R) = -2.00 (ρ_R = 1.00e-2)outflow velocity u_R = 4.00

The cyan curve is ρⁿ⁺¹ after one Euler step as a function of how much central-flux we mix into the right face. The red zone is where ρ has crossed the floor 10⁻³. The yellow dot is the largest θ that still lands above the floor — that is the θ the per-face limiter picks. Drop ρ_R below 10⁻³ or push u_R past 6 and θ collapses toward zero: the limiter falls back to first-order Lax–Friedrichs.

ρ_R을 10⁻³ 아래로 떨어뜨리면 곡선이 빠르게 빨간 영역으로 들어가고, 한정자가 고르는 θ가 0.2 아래로 주저앉는다. 그 자리에서는 사실상 1차 LF로 풀린다.

두 단계 한정 — 밀도 먼저, 압력은 그 다음#

밀도가 양수여도 압력이 음수면 솔버는 똑같이 죽는다 (음속이 허수가 된다). 그래서 Hu et al.은 두 패스로 나눈다.

1단계 (밀도 한정). 위 식으로 $\theta^\rho$ 를 면별로 결정. 모든 셀의 다음 스텝 밀도가 $\rho_{\min}$ 이상이 되도록.

2단계 (압력 한정). $\theta^\rho$ 로 밀도-flux를 고정한 뒤, 운동량·에너지 flux에 대해 또 다른 $\theta^p \leq \theta^\rho$ 를 구한다. 압력 양수 조건

p = (\gamma - 1)\bigl(E - \tfrac{1}{2}\rho u^2\bigr) \geq p_{\min}

은 비선형이라 일반적으로 이차방정식이지만, 보수적 근사(예: $\theta^p$ 를 작은 쪽으로 안전하게 잡는다)로 한 줄에 풀린다.

최종. $\theta = \min(\theta^\rho, \theta^p)$ 로 모든 보존 변수의 flux를 같은 가중치로 블렌딩. 한 면, 한 $\theta$ , 모든 보존량 일관성 유지.

50줄 NumPy — Le Blanc 충격관에서 θ 자동 결정#

Le Blanc 문제는 $\rho_L=1, p_L=100$ 대 $\rho_R=0.02, p_R=0.05$ 의 강한 충격관이다. 강한 팽창파가 좌측으로 진행하면서 첫 셀의 밀도가 음수로 떨어지는 고전적 시험 케이스다.

import numpy as np
 
GAMMA = 1.4
RHO_FLOOR = 1e-3
N, NSTEP = 200, 220
dx = 1.0 / N
dt = 9.0e-5
 
def cons_to_prim(U):
    rho, mu, E = U[0], U[1], U[2]   # 보존 → 원시 변수
    u = mu / rho
    p = (GAMMA - 1.0) * (E - 0.5 * rho * u * u)
    return u, p
 
def euler_flux(U):
    u, p = cons_to_prim(U)
    return np.array([U[1], U[1] * u + p, u * (U[2] + p)])
 
def alpha_face(UL, UR):
    uL, pL = cons_to_prim(UL); uR, pR = cons_to_prim(UR)
    aL = np.sqrt(GAMMA * max(pL, 0.0) / UL[0])
    aR = np.sqrt(GAMMA * max(pR, 0.0) / UR[0])
    return max(abs(uL) + aL, abs(uR) + aR)
 
def density_theta(U, dF, dt_dx):
    """면별 θ 결정 — ρ_LF 위에서 양수 유지 최대값."""
    theta = np.ones(N + 1)
    for f in range(1, N):
        d = dF[0, f]                          # 밀도 anti-diffusion
        if abs(d) < 1e-14: continue
        rL_room = U[0, f - 1] - RHO_FLOOR     # 좌 셀 여유
        rR_room = U[0, f]     - RHO_FLOOR     # 우 셀 여유
        cap = (rL_room if d > 0 else rR_room) / (dt_dx * abs(d))
        theta[f] = max(0.0, min(theta[f], cap))
    return theta

(전체 코드는 1차 LF 한 패스 + 중심차분 anti-diffusion 한 패스 + θ 결정 한 패스. 220 스텝, 200 셀. 한정자를 끄면 step 30 부근에서 NaN; 켜면 끝까지 살아남는다.)

위에서 만든 솔버를 시간축으로 직접 펴 보자.

1D strong shock tube · ρ_L=1, p_L=100 / ρ_R=0.02, p_R=0.05 · central + θ·anti-diffusion · 200 cells

t-index = 80 / 220show limited

Red is the raw centered scheme — its rarefaction tail dives below zero around t ≈ 0.006 and the solver blows up. Cyan is the same scheme blended with Lax–Friedrichs by a per-face θ ≤ θ_max so that ρ stays above the floor 10⁻³. Drag the time slider to watch the two solutions separate.

빨간 곡선(한정자 OFF)이 step 30 부근에서 ρ를 0 아래로 끌고 내려가는 순간을 잡을 수 있다. 청록은 같은 시점에 ρ ≈ 0.02 부근에서 평평하게 멈춰선다. 정확도 차이는 5% 미만이지만 한 쪽은 살고 한 쪽은 죽는다.

한정자가 가져가는 것#

θ가 0 으로 내려가는 면에서는 그 면이 1차 LF로 풀린다. 즉 충격 앞면이나 강한 팽창의 머리는 정확도가 떨어진다. 다행히 그 영역은 원래 전체 도메인의 5% 미만이라 전역 정확도는 거의 영향이 없다. 평균적으로 θ는 0.95~1.0 사이를 오간다.

함정 두 가지. 첫째, $\rho_{\min}$ 을 너무 작게 잡으면(예: $10^{-12}$ ) 한정자가 너무 늦게 발동해서 같은 NaN이 재발한다. 실무 기준은 시뮬레이션 평균 밀도의 $10^{-4}$ 정도. 둘째, RK 다단계 시간적분에서 매 단계마다 한정자를 재실행해야 한다 — 첫 단계만 통과시키면 두 번째 단계에서 그대로 발산한다.

내일도 기억할 것#

음수 밀도는 어디서나 새어나올 수 있다 — 강한 충격·팽창·접촉 불연속에서.
LF는 정답은 아니지만 안전 마진이다 — 그 마진을 기준으로 고차 보정을 한 칸씩만 허용한다.
밀도 먼저, 압력은 그 다음 — 두 패스로 끝난다. 더 복잡한 한정자는 더 많은 미세 조정을 요구할 뿐이다.

참고: X. Y. Hu, N. A. Adams, C.-W. Shu. Positivity-preserving method for high-order conservative schemes solving compressible Euler equations. J. Comput. Phys. 242 (2013) 169–180.