1950년, von Neumann이 충격파를 일부러 흐려놓았다

밀도가 1과 2로 톱니를 타도 두 칸 차분 압력은 0이다. collocated 격자는 그 힘을 보지 못한다.

1950년 Los Alamos. John von Neumann과 Robert Richtmyer가 Journal of Applied Physics에 짧은 논문을 실었다. 제목은 "A Method for the Numerical Calculation of Hydrodynamic Shocks". 충격이 격자 위에 한 칸으로 안 앉을 거면, 먼저 흐려 놓자는 한 줄이었다.

왼쪽 $\rho = 1, p = 1$ , 오른쪽 $\rho = 0.125, p = 0.1$ 인 Sod 관이다. ξ 슬라이더를 0에서 5로 올리면 충격 뒤 sawtooth가 잠기고, 불연속은 셀 몇 칸으로 퍼진다.

ξ artificial viscosity0.0

Sod shock tube. ξ=0: post-shock oscillations. ξ≈2-3: shock smeared over a few cells, wiggles damped.

ξ = 0이면 충격 뒷면이 톱니이고, 2 근처에서 진동은 죽으며 충격은 두세 칸으로 남고, 5면 plateau까지 흐물거린다. von Neumann이 권장한 $\xi \approx 2$ ~ $3$ 이 그 사이 자리다.

인공점성이 왜 먹히는지, 두께를 얼마나 남기는지, collocated 격자의 odd-even까지가 이 노트의 범위다. 끝에서 NumPy로 Mach 2 관을 한 번 잡는다.

격자 위의 충격파는 진동한다#

이상화된 Euler 방정식에는 점성이 없다. 그래서 진짜 충격파는 분자 평균자유경로 정도의 폭으로 압축된다. CFD 격자 셀은 그 폭의 백만 배쯤 된다. 셀 한 칸 안에서 밀도와 압력이 두 배로 뛰어야 한다는 뜻이다.

이 불연속을 어떤 보존적 차분 스킴에 통과시키면 거의 반드시 발생하는 현상이 있다. 충격 뒷면에서 솟구치는 sawtooth 모양의 진동, 흔히 wiggle이라 부르는 것이다. von Neumann은 1944년 Manhattan 프로젝트 시절 이미 이 문제를 경험했다. 폭발 시뮬레이션이 충격파 뒷부분에서 들쭉날쭉해지면, TNT 위력 계산이 통째로 어긋났다.

진동의 원인은 단순하다. 1차 upwind 같은 강한 수치 점성은 충격을 잡지만 다른 모든 흐름까지 무디게 만든다. 2차 이상의 스킴은 정확하지만 불연속 근처에서는 Gibbs 현상처럼 진동한다. 자연은 분자 점성으로 충격 안에서 엔트로피를 생성하는데, 격자는 그 통로가 없다. 칸 하나의 충격관을 리만 문제로 푸는 길은 Godunov가 칸마다 충격관을 넣는 방법이다. 인공점성은 그 불연속을 세 칸으로 늘려 진동을 죽인다. 불연속을 대하는 방식이 다르다.

1950년의 한 줄#

von Neumann과 Richtmyer의 답은 자연을 흉내내는 것이었다. 분자 점성을 직접 모사할 수는 없으니, 셀 폭 정도로 부풀린 가짜 점성 항을 압력에 더한다.

Q_i = \begin{cases} \xi^2 \, (u_{i+1} - u_{i-1})^2 \, \rho_i & \text{if } u_{i+1} \le u_{i-1} \\ 0 & \text{otherwise} \end{cases}

$Q_i$ 는 셀 $i$ 에서의 인공 bulk pressure, $\xi$ 는 충격을 몇 셀로 펼칠지 정하는 튜닝 파라미터(보통 2~3), $u_{i+1} - u_{i-1}$ 은 속도의 1차 차분, $\rho_i$ 는 셀 밀도. 모든 모멘텀·에너지 식에서 $P$ 를 $P + Q$ 로 치환하기만 하면 된다.

$Q$ 는 $(\Delta u)^2$ 에 비례한다. 매끈한 흐름에서는 거의 0이고, 충격 근처에서만 커진다. 켜는 조건은 $u_{i+1} \le u_{i-1}$ , 속도가 한쪽으로 모이는 칸만이다. 팽창파와 음파는 그대로 둔다. $\rho_i$ 를 곱하면 무거운 매질에 압력 보강이 더 실린다.

본문 식의 ξ와 코드의 0.25는 같은 두께가 아니다#

줄여 쓰면 $Q = \xi^2 (\Delta u)^2 \rho$ 다. 켜지는 쪽은 $\Delta u < 0$ 일 때뿐이다. 아래 목록은 0.25 * xi**2 * du**2 * rho를 쓴다. du는 두 칸을 건너뛴 차분이다. 셀 하나의 기울기로 보면 분모에 $2\Delta x$ 가 들어가고, 제곱하면 $1/4$ 이 붙는다. 코드의 0.25는 그 중심차분이다. 위의 식은 그 계수를 빠뜨린 채 $\xi^2$ 만 적었다.

$\xi = 2.5$ 이면 $\xi^2 = 6.25$ 다. 코드가 실제로 곱하는 수는 $0.25 \times 6.25 = 1.5625$ 다. 0.25 없이 같은 $\xi$ 를 식에 넣으면 $Q$ 는 네 배다. $Q$ 는 $\xi^2$ 에 비례하니, 네 배는 $\xi$ 를 두 배 키운 것과 같다. 슬라이더에서 $\xi = 2$ 로 잠잠해진 충격이, 식을 그대로 베낀 구현에서는 $\xi = 4$ 에 해당한다. $\xi = 5$ 면 plateau까지 흐물거린다. $\xi = 2.5$ 를 0.25 없이 넣으면 그 둔한 쪽에 붙는다.

격자 200칸, 시간 0.2, CFL 0.45는 목록에 적힌 그대로다. 출력 문장은 피크 밀도가 Mach 2의 Rankine–Hugoniot 값 2.667에서 한 자릿수 % 안이라고만 말한다. 0.25를 뺀 식을 그 문장으로 확인하면, 재는 충격 두께가 달라진다.

NumPy로 따라가는 인공점성#

50줄 안에 핵심 알고리즘을 담는다. donor-cell 이류 + 인공점성 항만 뽑아낸 최소 구현이다.

import numpy as np
 
def vnr_shock_solver(N=200, T=0.2, xi=2.5, gamma=1.4):
    dx = 1.0 / N
    rho = np.where(np.arange(N) < N // 2, 1.0, 0.125)
    mom = np.zeros(N)
    p   = np.where(np.arange(N) < N // 2, 1.0, 0.1)
    E   = p / (gamma - 1)
    t = 0.0
    while t < T:
        u = mom / rho
        c = np.sqrt(gamma * p / np.maximum(rho, 1e-9))
        dt = 0.45 * dx / np.max(np.abs(u) + c)
        # von Neumann-Richtmyer Q
        du = np.zeros(N)
        du[1:-1] = u[2:] - u[:-2]
        Q = np.where(du < 0, 0.25 * xi**2 * du**2 * rho, 0.0)
        # donor-cell advection of rho, mom, E
        uf = 0.5 * (u[:-1] + u[1:])
        for q_arr, q_new in [(rho, 'rho'), (mom, 'mom'), (E, 'E')]:
            flux = np.where(uf > 0, q_arr[:-1] * uf, q_arr[1:] * uf)
            q_arr[1:-1] -= dt / dx * (flux[1:] - flux[:-1])
        # pressure update + body force from P+Q
        p = (gamma - 1) * (E - 0.5 * rho * (mom / rho) ** 2)
        Ptot = p + Q
        mom[1:-1] -= dt / (2 * dx) * (Ptot[2:] - Ptot[:-2])
        E[1:-1]   -= dt / (2 * dx) * (Ptot[2:] * u[2:] - Ptot[:-2] * u[:-2])
        t += dt
    return rho, mom, E
 
rho, mom, E = vnr_shock_solver(xi=2.5)
print(f"shock peak rho ≈ {rho.max():.3f}   (Rankine-Hugoniot exact ≈ 2.667 for Mach 2)")

실행하면 충격 peak이 RH 해의 한 자릿수% 이내로 떨어진다. ξ = 0으로 같은 코드를 돌리면 진동 때문에 max·min이 비대칭으로 솟는다.

odd-even 디커플링이라는 또 다른 함정#

인공점성으로 충격을 잡았다고 끝이 아니다. 같은 collocated 격자(모든 변수를 셀 중심에 배치)에는 von Neumann이 알아챈 또 다른 병이 있다. 다음 정지 상태를 보자. 속도는 0, 비열에너지는 1, 그러나 밀도가 sawtooth로 진동한다.

\rho_i = \begin{cases} 1 & i = 1, 3, 5, \dots \\ 2 & i = 2, 4, 6, \dots \end{cases}

압력 $p_i = \rho_i$ (isothermal 가정)이면 셀 $i = 2$ 가 이웃보다 두 배 압력이다. 직관적으로는 즉시 풀려야 한다. 그런데 셀 $i$ 의 모멘텀 업데이트는

\frac{\partial (\rho u)_i}{\partial t} = -\frac{F_{i+1/2} - F_{i-1/2}}{\Delta x} = -\frac{p_{i+1} - p_{i-1}}{2\Delta x} = 0

가 된다. $p_{i+1}$ 과 $p_{i-1}$ 이 같은 격자 부분배열에서 나오므로 차이가 정확히 0이다. 격자가 두 개의 독립 부분격자(짝수 셀, 홀수 셀)로 갈라져 서로를 못 보는 셈이다. 이 sawtooth 모드는 자라지도, 줄지도 않는다.

Initial: ρ alternates 1↔2 (isothermal, p=ρ). Collocated: ∂p/∂x at cell center cancels — sawtooth survives. Staggered: face gradient sees the jump and damps it.

위에서 mode를 collocated → staggered로 바꿔 보자. staggered 격자는 모멘텀을 셀 면(face)에 둔다. 면 $i+1/2$ 의 압력 차이는 $p_{i+1} - p_i$ — 짝·홀 부분격자가 한 칸 이웃 사이를 본다. sawtooth는 곧바로 1차 차분의 사거리에 들어와 damping된다.

무엇을 잃고 무엇을 얻었나#

von Neumann의 1950년 트릭은 두 가지를 동시에 보여 줬다. 하나는 격자 한계를 모사로 메우는 길이 가능하다는 것. 다른 하나는 그 모사가 또 다른 격자 병을 만든다는 것. 두 사실 모두 현대 CFD에 그대로 살아 있다.

오늘날 ZEUS, Athena, PLUTO 같은 astrophysics 코드는 staggered 격자 + 인공점성을 그대로 쓴다. 한편 OpenFOAM이나 ANSYS Fluent는 collocated 격자에 Rhie-Chow 보간을 얹어 odd-even 디커플링을 우회한다. 둘은 같은 문제를 다른 각도에서 푼 결과다.

격자에서 충격을 본 적이 있다면, 그 충격 뒤의 sawtooth는 버그가 아니라 1944년부터 알려진 격자 한계의 다른 얼굴이다.

0.25를 뺀 ξ=2.5는 슬라이더의 ξ=5 쪽에 붙는다#

0.25를 빼고 $\xi = 2.5$ 를 본문 식에 넣으면, 코드가 $\xi = 5$ 로 켠 $Q$ 와 크기가 같아진다. 슬라이더에서 그 자리는 이미 충격이 둔하다. $\xi \approx 2$ ~ $3$ 은 코드의 0.25가 붙어 있을 때의 자리다. 계수를 뺀 채 그 숫자를 옮기면 두께 이야기가 바뀐다.

짝수 칸과 홀수 칸이 서로를 못 보는 문제는 인공점성과 별개다. 속도를 셀 면에 두면 한 칸 이웃의 압력을 본다. 그 배치가 충격 뒤 sawtooth와는 다른 오차로 어긋나는 기록은 컬 구속이 엇갈린 격자에서 깨진 경우에 있다.

1950년, von Neumann이 충격파를 일부러 흐려놓았다 — 인공점성과 odd-even 디커플링