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

다음 스텝 밀도가 바닥 위에 남는 한에서만, 고차 플럭스와 LF 사이의 간격을 θ로 연다.

오전 2시 27분. Le Blanc 충격관이 0.003초 만에 NaN을 토했다. CFL을 0.1까지 낮춰도, 격자를 두 배로 해도, 곡선 두 개가 같이 무너졌다. 원인은 단순했다. 강한 팽창파(rarefaction, 압력이 빠르게 떨어지는 영역) 머리에서 한 셀의 밀도가 음수로 떨어졌다. 두 스텝 뒤 음속이 NaN이 되고, 그다음 스텝에 도메인 전체가 죽는다. Hu·Adams·Shu(2013)의 θ-블렌드는 그 자리에서 고차 플럭스를 한 칸만 뒤로 물린다.

오른쪽 밀도 ρ_R을 10⁻³ 아래로 내리면 허용 곡선이 빨간 구역으로 들어가고, 그 면의 θ는 0.2 밑으로 떨어진다.

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.

빨간 구역은 다음 스텝 밀도가 바닥을 뚫는 쪽이다. θ가 그 경계에서 멈추면 그 면은 사실상 1차 Lax–Friedrichs다.

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

고차 스킴은 안정성과 정확도 사이에서 흔들린다. 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)이다. 이웃을 평균해서 톱니를 죽인다. 발상은 충격에 일부러 점성을 바른 1950년의 처리와 같은 자리다. Lax는 CFL $\leq$ 1 아래에서 $\rho^{n+1} \geq 0$ 임을 보였다. Hu et al. 정리 1을 단순화한 형태다. 충격이 세든 약하든 LF는 안 무너진다.

대가는 정확도다. 두 칸짜리 충격이 다섯 칸이 되고, 접촉 불연속(contact discontinuity)은 스텝마다 한 셀씩 더 흐려진다. 실무에서 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$ 의 상한을 면마다 푼다. 한 면은 양쪽 셀에 닿는다. 둘 다 바닥 위여야 하니 더 작은 쪽이 한계다.

앞에서 ρ_R을 내렸을 때 θ가 0.2 아래로 주저앉던 것이 이 상한이다. 그 면은 사실상 1차 LF로 풀린다.

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

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

밀도 패스. 위 식으로 $\theta^\rho$ 를 면마다 정한다. 모든 셀의 다음 밀도가 $\rho_{\min}$ 이상이다.

압력 패스. $\theta^\rho$ 로 밀도 플럭스를 고정한 뒤, 운동량·에너지 플럭스에 $\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)$ 로 모든 보존 플럭스를 같은 무게로 섞는다. 한 면, 한 $\theta$ . 보존량 사이의 일관성이 유지된다.

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

Le Blanc은 $\rho_L=1$ , $p_L=100$ 대 $\rho_R=0.02$ , $p_R=0.05$ 인 강한 충격관이다. 강한 팽창이 왼쪽으로 가면서 첫 셀 밀도가 음수로 떨어지는 시험이다.

적힌 좌우를 나누면 밀도비는 $1/0.02 = 50$ , 압력비는 $100/0.05 = 2000$ 이다. 앞에서 100:1을 위험선의 예로 들었지만, 이 쌍의 밀도 점프는 그 선 아래다. 센 쪽은 압력이다. $\gamma = 1.4$ 이면 음속 $a = \sqrt{\gamma p / \rho}$ 라 왼쪽은 $\sqrt{140} \approx 11.83$ , 오른쪽은 $\sqrt{3.5} \approx 1.871$ 이다. 둘의 비는 $\sqrt{40} \approx 6.32$ 로 딱 떨어진다.

격자도 식에 넣으면 시간 간격 탓이 아님이 보인다. $N = 200$ 이라 $\Delta x = 1/200 = 0.005$ 이고, $\Delta t = 9.0 \times 10^{-5}$ 이다. $\Delta t/\Delta x = 0.018$ 이다. 왼쪽 음속만으로 잰 쿠란트 수는 $a_L \Delta t/\Delta x \approx 0.213$ 이다. 새벽에 CFL을 0.1로 내려도 죽었던 것과 맞다. $\rho^{n+1} = \rho^n - (\Delta t/\Delta x)\Delta F$ 에서 플럭스 차가 밀도 여유를 넘으면, 스텝을 줄여도 바닥을 뚫는다.

오른쪽 초기 밀도 $0.02$ 에서 바닥 $10^{-3}$ 까지 여유는 $0.019$ 다. $\theta$ 가 1에 머물려면 밀도 보정 $|d|$ 가 $0.019/0.018 \approx 1.06$ 을 넘으면 안 된다. $|d| = 2$ 이면 상한은 $0.019/(0.018 \times 2) \approx 0.53$ 이다. 여유가 $5 \times 10^{-4}$ 까지 줄면, 같은 $|d| = 1$ 에서 $\theta$ 상한은 $0.0005/0.018 \approx 0.028$ 이다. 그 면은 거의 LF다.

블렌드를 가중치로 다시 쓰면 $F = (1-\theta)F^{LF} + \theta F^{HO}$ 다. $\theta = 0.53$ 이면 고차 몫은 53%다. $\theta = 0.028$ 이면 고차 몫은 2.8%뿐이고 나머지는 LF다. 코드의 상한은 업데이트와 같은 부등식이다. 보정이 밀도를 움직이는 양은 $(\Delta t/\Delta x)\,\theta\,|d|$ 이고, 오른쪽 여유 $0.019$ 를 넘으면 안 된다. 그래서 $\theta \le 0.019 / (0.018\,|d|)$ 다. $|d| \approx 1.06$ 에서 등호가 $\theta = 1$ 에 닿는다.

압력 식도 같은 방식으로 읽힌다. $p = (\gamma - 1)(E - \tfrac{1}{2}\rho u^2)$ 이고 $\gamma - 1 = 0.4$ 다. $E$ 와 $\rho$ 를 고정하면 운동에너지 1이 압력을 0.4만큼 깎는다. 오른쪽 $p = 0.05$ 는 운동에너지 $0.05/0.4 = 0.125$ 에서 0에 닿는다. 왼쪽 $p = 100$ 은 같은 계산으로 $250$ 까지 버틴다. 밀도 바닥만 지키고 이 몫을 빼먹으면, 밀도는 양수인데 음속이 허수가 되는 쪽으로 샌다.

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셀. 한정자를 끄면 30스텝 부근에서 NaN이고, 켜면 끝까지 남는다.

한정자 스위치를 끄면 빨간 곡선이 30스텝 근처에서 밀도를 0 아래로 끌고 간다.

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.

그 순간이 붕괴다. 청록은 같은 시각에 $\rho \approx 0.02$ 근처에서 평평해진다. 정확도 차이는 5% 미만인데, 한쪽은 살고 한쪽은 죽는다.

한정자가 가져가는 것#

$\theta$ 가 0으로 내려가는 면은 1차 LF로 풀린다. 충격 앞면과 강한 팽창의 머리는 거기서 무뎌진다. 그 구역은 도메인의 5% 미만이라 전역 정확도는 거의 그대로다. 평균 $\theta$ 는 0.95에서 1.0 사이를 오간다.

$\rho_{\min}$ 을 $10^{-12}$ 처럼 너무 낮추면 한정자가 늦게 깨어나서 같은 NaN이 다시 난다. 실무에서는 평균 밀도의 $10^{-4}$ 근처를 바닥으로 둔다. RK 다단계에서는 단계마다 한정자를 다시 돌려야 한다. 첫 단계만 통과시키면 다음 단계에서 그대로 발산한다.

바닥을 1e-12로 내리면 한정자는 NaN 다음에 깨어난다#

$\rho_{\min} = 10^{-12}$ 는 바닥이 아니라 장식이다. 오른쪽 밀도 0.02에서 보면 여유의 거의 전부가 아직 남아 있고, 위의 식으로 $\theta$ 상한은 1 근처에 머문다. 밀도가 이미 $10^{-8}$ 아래로 꺼진 뒤에야 상한이 조여진다. 음속 제곱 $\gamma p / \rho$ 는 그 전에 부호가 뒤집힌다. 한정자는 죽은 셀을 보고 나서야 LF로 물러난다.

면의 $\theta$ 는 셀 평균만 바닥 위에 둔다. 평균은 양수인데 셀 안 재구성이 한 점에서 음수로 꺾이면 이 한 줄은 그 점을 보지 못한다. 그 구멍을 평균은 지키는데 점이 음수로 꺾이는 재구성이 따로 막는지가 다음 선택이다.

참고: 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.