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

post-shock 진동을 잠재운 ξ²(Δu)² 트릭과 그것이 남긴 또 하나의 문제

1950년 Los Alamos. John von Neumann과 Robert Richtmyer가 한 편의 짧은 논문을 Journal of Applied Physics에 실었다. 제목은 "A Method for the Numerical Calculation of Hydrodynamic Shocks". 본문의 트릭은 한 줄짜리였다 — 충격파가 격자 위에서 깔끔하게 떨어지지 않을 거라면, 우리가 먼저 그것을 의도적으로 흐려 놓자. 이 글은 그 한 줄의 인공점성(artificial viscosity)이 왜 작동하는지, 어떤 대가를 치르는지, 그리고 그것이 collocated 격자에 남긴 또 하나의 문제 — odd-even 디커플링 — 까지를 따라간다. 마지막에는 NumPy 50줄로 Mach 2 충격관을 직접 잡아 본다.

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

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

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

진동의 원인은 단순하다. 1차 upwind 같은 강한 수치 점성은 충격을 잡지만 다른 모든 흐름까지 무디게 만든다. 2차 이상의 스킴은 정확하지만 불연속 근처에서는 Gibbs 현상처럼 진동한다. 자연은 분자 점성으로 충격 안에서 엔트로피를 생성하는데, 격자는 그 통로가 없다.

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$ 를 곱해 무거운 매질에서 더 큰 압력 부스트를 준다.

직접 만져 보자#

아래 시뮬레이션은 고전적인 Sod 충격관이다. 왼쪽 $\rho = 1, p = 1$ , 오른쪽 $\rho = 0.125, p = 0.1$ 로 시작해 막을 떼는 순간을 흉내낸다. ξ 슬라이더로 인공점성 강도를 0부터 5까지 직접 조절한다.

ξ artificial viscosity0.0

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

ξ = 0에서 충격 뒷면에 작은 sawtooth가 보인다. ξ = 2 근처에서 진동은 잠잠해지지만 충격이 두세 셀로 펴진다. ξ = 5까지 올리면 충격 자체가 너무 둔해져서 plateau 근처가 흐물거린다. von Neumann이 권장한 $\xi \approx 2$ ~ $3$ 이 trade-off의 sweet spot이라는 것을 눈으로 확인할 수 있다.

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년부터 알려진 격자 한계의 다른 얼굴이다.

기억할 점#

**인공점성 $Q = \xi^2 (\Delta u)^2 \rho$ **는 압축( $\Delta u < 0$ )에서만 켜진다. 부드러운 흐름은 건드리지 않고 충격만 3셀 정도로 펴서 진동을 잠재운다.
trade-off은 명백하다. ξ가 작으면 진동이 남고, 크면 충격이 둔해진다. $\xi \approx 2$ ~ $3$ 이 1950년대부터 굳어진 황금비.
odd-even 디커플링은 collocated 격자의 구조적 결함이다. staggered 격자(ZEUS류) 또는 Rhie-Chow(OpenFOAM류)로 우회한다 — 둘 중 어떤 풀이를 보고 있는지 코드를 열어 확인하는 습관이 매번 유용하다.