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

별표 압력이 좌우 압력을 가르는 쪽으로 파를 고르고, 부채가 겹치는 스텝은 버린다.

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

슬라이더를 1.0 기준으로 넘겨 보자. 인접한 두 부채가 겹치면 빨간 경고가 켜지고, 1.0 아래에서는 경계 플럭스가 시간상 상수로 남는다.

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|λ|.

경고가 켜진 스텝은 자기 닮음이 깨진 스텝이다. $\Delta t$ 가 $\Delta x/\max|\lambda|$ 를 넘는 그림과 같다.

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$ 에서 좌·우 어느 쪽인지, 각 파가 충격인지 팽창인지에 따라 상태를 추출한다. 그 추출이 수치 플럭스로 굳는 경로는 Riemann 문제에서 Godunov로 이어지는 플럭스에 따로 있다.

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

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

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

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

앞의 빨간 경고가 이 부등식이 깨진 자리다.

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)이다.

인쇄된 두 상태만으로 별표값을 가른다. 좌측 압력 3.528, 우측 0.571, 별표 압력 3.40948. 왼쪽에서는 3.528에서 0.11852만 내려오고, 오른쪽에서는 0.571에서 2.83848이 올라온다. 별표 속도 0.62556은 좌측 0.698보다 0.07244 낮고, 우측 0보다는 위에 있다. 좌우 속도의 산술평균은 0.349다. 속도 식의 보정항이 0.349를 0.62556까지 밀어 올린다.

분기만 보면 충분하다. 별표 압력이 좌측 압력보다 작으므로 왼쪽 $f_L$ 은 팽창 가지다. 우측 압력보다 크므로 오른쪽 $f_R$ 은 충격 가지다. 허용오차는 $10^{-9}$ , 반복 상한은 80, 압력 바닥은 $10^{-6}$ 이다. 초기 추정은 $0.5(p_L+p_R)$ 에서 속도 점프를 더한다. $u_R-u_L=-0.698$ 이라 그 보정은 양수다. 두 압력의 평균은 2.0495. 추정은 그 평균보다 위에서 출발해 3.40948에 닿는다. 이 뺄셈은 코드에 적힌 식과 출력 주석만 쓴 것이다.

100칸에서 1차 업윈드가 충격을 다섯 칸으로 펴면, 그 폭은 격자 길이의 5%다. 슬라이더 1.0은 $\Delta t = \Delta x/\max|\lambda|$ 의 경계다. 1을 넘기면 인접 부채가 겹치고, 셀 경계 플럭스가 한 스텝 동안 상수라는 가정이 깨진다. $\gamma=1.4$ 는 $A_K$ , $B_K$ , $c_K$ 세 곳에 같이 들어간다.

부채를 직접 펴 보자#

아래 시뮬레이션은 위 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 발상의 골격은 그대로 살아남았다.

정확해를 셀 경계에 넣을 때#

$p^{\star}=3.40948$ 은 좌측 3.528보다 작고 우측 0.571보다 크다. 왼쪽은 팽창, 오른쪽은 충격이다. $u^{\star}=0.62556$ 은 좌측 0.698에서 0.07244 내려온 값이다. 이 한 쌍이 셀 경계마다 필요하면 정확해를 돌린다. 격자 100개에서 충격이 다섯 칸으로 퍼져도 되면 1차 업윈드로 충분하다.

정확해는 비싸다. 1981년 이후 Roe, HLL, HLLC가 그 자리를 차지한 이유다. 검증용 기준해로는 이 $p^{\star}$ , $u^{\star}$ 가 남아 있다. 충격을 일부러 두껍게 만드는 반대편은 1950년 인공점성이 충격을 흐린 방식이다. Godunov는 흐리지 않고 포착하고, 그 대가로 $\Delta t$ 가 부채에 묶인다.