Godunov 기법 — 격자 한 칸 안에 충격관 하나

Riemann 부채와 정확해 Lax 충격관으로 본 Godunov 방법의 원리와 한계.

100개 격자를 박아도 1차 업윈드는 충격파를 다섯 칸 두께로 뭉갠다. 1959년 Sergei Godunov는 반대를 택했다. 격자 한 칸 한 칸의 경계에 작은 충격관(shock tube) 하나를 놓고, 그 정확한 해석해를 한 시간 스텝 동안 흘려보낸 다음 셀 평균을 다시 잡는다. 이 글은 Riemann 문제(좌·우 두 상태가 격막을 두고 마주 본 초기조건)의 다섯 영역 부채 구조를 정확해로 풀어 보고, Godunov 한 스텝의 시간 제약과 이 발상이 닿은 1차 정확도의 벽까지 따라간다.

1959년, 격자 한 칸 안에 충격관 하나#

지금까지의 격자식 해법은 셀 내부가 일정하다고 가정한다. 그러면 셀 경계 $x_{i+1/2}$ 에는 좌측 상태 $\mathbf{q}_i$ 와 우측 상태 $\mathbf{q}_{i+1}$ 사이의 도약이 생긴다. 이 도약을 잠시 우주 전체에 펼쳐 놓으면 — 좌측은 $-\infty$ 까지 $\mathbf{q}_i$ , 우측은 $+\infty$ 까지 $\mathbf{q}_{i+1}$ — 그것이 곧 Riemann 문제다.

Godunov의 아이디어는 단순했다. 모든 셀 경계에서 이 작은 Riemann 문제를 풀어, 그 정확한 자기 닮음(self-similar) 해를 한 시간 스텝 동안 흘려보낸다. 끝나면 셀 평균을 다시 잡는다. 압축성 Euler 방정식의 고유 구조(특성속도 $u, u \pm c_s$ )를 분리하지 않고 한꺼번에 다룬다는 점이 핵심이다.

Riemann 부채: 다섯 영역의 자기 닮음 해#

$\gamma$ -기체의 Euler 방정식

\partial_t \mathbf{q} + \partial_x \mathbf{f}(\mathbf{q}) = 0, \quad \mathbf{q} = (\rho, \rho u, \rho E)^\top

여기서 $\rho$ 는 밀도, $u$ 는 속도, $E$ 는 총 에너지, $\mathbf{f}$ 는 보존 플럭스.

Riemann 초기조건 $(\rho_L, u_L, p_L)$ 과 $(\rho_R, u_R, p_R)$ 이 주어지면 해는 $\xi = x/t$ 한 변수만의 함수가 된다. 격막에서 펼쳐지는 부채는 정확히 다섯 영역으로 나뉜다.

좌측 미교란 영역 — 좌 측 파(rarefaction/shock)의 머리(head)가 아직 닿지 않은 곳.
좌측 파의 내부 — 팽창파(rarefaction wave)면 연속 변화, 충격파면 점프.
별표 좌(star-L) 영역 — 속도와 압력은 별표값으로, 밀도는 좌측 정보로.
별표 우(star-R) 영역 — 같은 별표 압력·속도, 다른 별표 밀도.
우측 미교란 영역.

3번과 4번 사이의 경계가 접촉 불연속(contact discontinuity)이다. 좌우 별표 압력·속도는 같지만 밀도와 엔트로피만 점프한다.

$p^{\star}$ 한 변수의 함수방정식#

세 미지수 $(\rho^{\star}_L, \rho^{\star}_R, p^{\star}, u^{\star})$ 가 있지만, 별표 압력 $p^{\star}$ 하나만 풀면 나머지가 따라온다. Toro의 표준 정리는 다음 함수방정식에 닫힌다.

f_L(p^{\star}; \mathbf{q}_L) + f_R(p^{\star}; \mathbf{q}_R) + (u_R - u_L) = 0

여기서 좌·우 파별 $f_K$ 는 충격파면 Rankine–Hugoniot, 팽창파면 등엔트로피 관계로 갈린다.

f_K(p) = \begin{cases} (p - p_K)\sqrt{\dfrac{A_K}{p + B_K}} & p > p_K \;\;\text{(shock)} \\[6pt] \dfrac{2 c_K}{\gamma - 1}\left[\left(\dfrac{p}{p_K}\right)^{\!\frac{\gamma-1}{2\gamma}} - 1\right] & p \le p_K \;\;\text{(rarefaction)} \end{cases}

$A_K = 2/[(\gamma + 1)\rho_K]$ , $B_K = \frac{\gamma - 1}{\gamma + 1} p_K$ , $c_K = \sqrt{\gamma p_K / \rho_K}$ . 별표 속도는

u^{\star} = \tfrac{1}{2}(u_L + u_R) + \tfrac{1}{2}\bigl[f_R(p^{\star}) - f_L(p^{\star})\bigr].

Newton 반복으로 $p^{\star}$ 를 얻은 뒤 $\xi = x/t$ 에서 좌·우 어느 쪽인지, 각 파가 충격인지 팽창인지에 따라 상태를 추출한다.

시간 스텝 한계 — 부채가 닿는 순간#

Godunov의 본질적 제약은 인접한 두 부채가 한 시간 스텝 안에 겹치면 안 된다는 것이다.

\Delta t \le \min_i \frac{\Delta x}{\max_k |\lambda_k^{(i+1/2)}|}

$\lambda_k$ 는 Euler 방정식의 세 고유값 $u, u \pm c_s$ . 부채가 닿는 순간 자기 닮음 가정이 깨지고 보존도 부서진다.

time step Δt / Δt_CFL = 0.70max wave speed = 1.00

Each cell interface spawns a self-similar Riemann fan. Push the time step past 1 and the fans collide — Godunov requires Δt ≤ Δx / max|λ|.

슬라이더를 1.0 너머로 올리면 두 인접 부채가 충돌해 빨간 경고가 켜진다. 1.0 미만이면 셀 경계의 플럭스는 시간상 상수라는 가정이 그대로 통한다.

Python: Lax 충격관 정확해#

Toro 책의 표준 절차를 짧게 옮긴다. Lax 충격관 $(\rho_L, u_L, p_L) = (0.445, 0.698, 3.528)$ , $(\rho_R, u_R, p_R) = (0.5, 0, 0.571)$ .

import numpy as np
 
GAMMA = 1.4
 
def f_wave(p, rho_K, p_K, c_K, gamma=GAMMA):
    # 단일 파의 압력-속도 함수 (충격/팽창 분기)
    if p > p_K:
        A = 2.0 / ((gamma + 1.0) * rho_K)
        B = (gamma - 1.0) / (gamma + 1.0) * p_K
        return (p - p_K) * np.sqrt(A / (p + B))
    return (2.0 * c_K / (gamma - 1.0)) * ((p / p_K) ** ((gamma - 1.0) / (2.0 * gamma)) - 1.0)
 
def f_prime(p, rho_K, p_K, c_K, gamma=GAMMA):
    if p > p_K:
        A = 2.0 / ((gamma + 1.0) * rho_K)
        B = (gamma - 1.0) / (gamma + 1.0) * p_K
        return np.sqrt(A / (p + B)) * (1.0 - (p - p_K) / (2.0 * (p + B)))
    return (1.0 / (rho_K * c_K)) * (p / p_K) ** (-(gamma + 1.0) / (2.0 * gamma))
 
def solve_star(left, right, gamma=GAMMA, tol=1e-9, max_iter=80):
    rho_L, u_L, p_L = left
    rho_R, u_R, p_R = right
    c_L = np.sqrt(gamma * p_L / rho_L)
    c_R = np.sqrt(gamma * p_R / rho_R)
    # 1차 추정: 두 압력의 평균에서 속도 도약 보정
    p = max(1e-6, 0.5 * (p_L + p_R) - 0.125 * (u_R - u_L) * (rho_L + rho_R) * (c_L + c_R))
    for _ in range(max_iter):
        f  = f_wave(p, rho_L, p_L, c_L) + f_wave(p, rho_R, p_R, c_R) + (u_R - u_L)
        df = f_prime(p, rho_L, p_L, c_L) + f_prime(p, rho_R, p_R, c_R)
        dp = f / df
        p_new = max(1e-6, p - dp)
        if 2.0 * abs(p_new - p) / (p_new + p) < tol:
            p = p_new
            break
        p = p_new
    u = 0.5 * (u_L + u_R) + 0.5 * (f_wave(p, rho_R, p_R, c_R) - f_wave(p, rho_L, p_L, c_L))
    return p, u, c_L, c_R
 
# Lax 충격관
left  = (0.445, 0.698, 3.528)
right = (0.5,   0.0,   0.571)
p_star, u_star, c_L, c_R = solve_star(left, right)
print(f"p* = {p_star:.5f}, u* = {u_star:.5f}")
# → p* = 3.40948, u* = 0.62556

이 12줄 Newton 반복이 Godunov 한 스텝의 핵심이다. 정확해 자체는 비쌌기 때문에 1981년 Roe가 선형화로 대체했고, 그 이후 HLL·HLLC가 표준이 되었다. 그러나 정확해는 여전히 검증의 기준점(reference solution)이다.

부채를 직접 펴 보자#

아래 시뮬레이션은 위 Newton 반복을 브라우저에서 다시 풀고, 시간 슬라이더로 자기 닮음 부채를 펼친다. preset에서 123 problem(좌우로 빠지는 진공 가까운 흐름)으로 바꾸면 가운데에 깊은 압력 골이 생기고, Strong-blast(폭발파)는 거의 단일 충격파에 가까운 모양으로 풀린다.

presettime t = 0.140

Three waves emerge from the origin: a left-going rarefaction or shock, a contact discontinuity at u*, and a right-going wave. The vertical pink dashed line marks the initial diaphragm at x = 0.

t를 작게 두고 시작해 점차 늘려 보면, 좌·우에서 동시에 출발한 파가 일정 속도로 멀어진다. Lax에서는 좌측이 팽창파(rarefaction)로 부드럽게 풀리고 우측은 충격파(shock)로 점프한다. 가운데의 작은 밀도 단차가 접촉 불연속이다.

1차의 벽 — Godunov 정리#

Godunov 본인은 같은 1959년 논문에서 선형이며 1차 초과의 정확도를 가지는 어떤 단조 스킴도 존재하지 않는다는 정리를 증명한다. 자기 자신의 발상이 1차에 묶인다는 음울한 결과였다.

이 벽을 우회하는 길이 곧 1970–80년대 고차 충격 포착 스킴의 역사다. van Leer의 MUSCL, Harten의 TVD, Roe의 선형화, Colella·Woodward의 PPM, Shu·Osher의 WENO — 모두 셀 내부를 일정이 아닌 비선형 재구성(reconstruction)으로 표현해 단조성을 보존하면서 정확도를 끌어올린다. Godunov 발상의 골격은 그대로 살아남았다.

마지막에 남기고 싶은 것#

셀 경계마다 정확한 Riemann 부채를 한 번씩 풀면, 충격파를 별도 처리 없이 자연스럽게 포착한다.
시간 스텝은 부채가 닿지 않게 $\Delta t \le \Delta x / \max|\lambda|$ 로 묶인다.
1차 정확도의 한계는 Godunov 본인의 정리로 못 박혔다. 이후의 모든 고차 충격 포착 스킴은 이 골격 위의 재구성 전쟁이다.