[논문 리뷰] 저Mach에서 Roe가 망가지는 이유 — Hope–Collins(2022) 인공확산 세 스케일링

음파를 남길 계산이면 (1,1) 인공확산은 M의 -2승이 아니라 -1승에 둔다.

1993년 봄, Volpe는 Mach 0.01에서 NACA 0012 날개를 풀었다. Roe 플럭스가 토해낸 압력장은 인비스시드 답과 닮은 구석이 없었다. 그 한 장의 그래프가 30년에 걸친 저Mach 인공확산 연구의 출발선이었다. Hope–Collins(2022)는 그 30년의 흩어진 처방을 세 갈래로 묶어 정리한다. 오늘은 그 세 갈래를 점근 분석과 60줄 NumPy로 직접 통과시킨다.

논문: J. Hope-Collins, L. di Mare. Artificial diffusion for convective and acoustic low Mach number flows I: Analysis of the modified equations, and application to Roe-type schemes. J. Comput. Phys. (2022). 단일상 압축성 Euler, Roe 계열.

$M$ 슬라이더를 $10^{-2}$ 와 $10^{-3}$ 으로 내려 보자. 세 행렬의 (1,1) 칸이 같은 Mach에서 얼마나 벌어지는지 칸 크기로 보인다.

How each scaling explodes as M → 0

Convective (Aᶜ)

M-210000

M01.00

M-210000

M01.00

Acoustic (Aᵃ)

M-1100

M01.00

M-210000

M-1100

Mixed (Aᵐ)

M-1100

M01.00

M-210000

M01.00

log₁₀ M = -2.00 (M = 1.00e-2)

The (1,1) and (2,1) cells of the convective matrix grow as 10⁴ when M = 10⁻². That is the catastrophe Hope–Collins traces to a parabolic limit equation for acoustic pressure.

$M=10^{-2}$ 에서 convective의 (1,1)은 $10^4$ , mixed는 $10^2$ 다. 100배다. acoustic의 (2,2)는 $M^{-1}$ 이라 대류 속도를 먼저 지운다.

1993년, 충격이었던 Volpe 그래프#

Volpe와 Godfrey가 같은 해에 보고한 것은 한 줄로 요약된다 — 고속 흐름용으로 설계된 표준 Roe 스킴을 그대로 Mach 0.01에 끌고 가면 압력장이 비물리적으로 깨진다. Godfrey는 "전제조건 행렬에 맞춘 인공확산"으로 그래프를 살려냈다. 이후 30년 동안 AUSM, Zha–Bilgen, Roe–LM, Thornber 등 수많은 변형이 등장했다.

왜 깨질까. Euler 방정식의 속도 벡터에는 $M^{-2}$ 계수의 압력 구배 항이 있다. Mach가 작아질수록 이 항이 폭주하고, 균형을 맞추는 다른 항이 압력의 위계(asymptotic order)를 강제한다. 위계가 어긋난 인공확산 한 줄은 곧바로 정확도를 부순다.

이 논문의 미덕은 이산화 무관이다. 유한체적이든 유한차분이든, central이든 upwind든, "어떤 종류의 인공확산을 어디에 넣었는가"만으로 그 스킴이 저Mach에서 어떻게 행동할지가 결정된다.

세 가지 저Mach 극한 — convective, acoustic, mixed#

비차원화한 1D 엔트로피 변수 Euler 방정식은 다음과 같다.

\partial_t p + u\,\partial_x p + \gamma p\,\partial_x u = A_{11}\partial_{xx}p + A_{12}\partial_{xx}u

\rho\,\partial_t u + M^{-2}\,\partial_x p + \rho u\,\partial_x u = A_{21}\partial_{xx}p + A_{22}\partial_{xx}u

여기서 $p,u,\rho$ 는 압력·속도·밀도, $M$ 은 기준 Mach 수, $A_{ij}$ 는 인공확산 계수 행렬이다. 우변이 0이면 원래 Euler. 우변의 위계를 어떻게 잡느냐가 저Mach 거동을 결정한다.

각 변수를 $\psi = \psi^{(0)} + M\psi^{(1)} + M^2\psi^{(2)} + \dots$ 로 전개해 $M$ 의 거듭제곱별로 모으면, 단일 시간척도 분석은 세 극한을 토해낸다.

Convective 극한: $\partial_x p^{(0)} = 0$ , $\partial_x p^{(1)} = 0$ . 압력은 $O(M^2)$ 변동만 허용. 발산은 0. — 비압축성에 점근.
Acoustic 극한: $\partial_\tau u^{(0)} + \partial_x p^{(1)}/\rho^{(0)} = 0$ . 음파만 살아남는 선형 음향.
Mixed 극한: 짧은 공간 스케일의 음파가 긴 공간 스케일의 대류 위에 얹혀 공존하는 다중척도 흐름.

세 극한은 서로 다른 압력 위계를 강제한다. 인공확산 행렬은 그 위계와 자기 자신을 동시에 맞춰야 한다.

인공확산 행렬을 세 갈래로 설계한다#

Turkel(1969 → 1994)의 절차는 단순하다. ① 좌변( $\mathcal{L}$ )에서 $M\to 0$ 일 때 가장 큰 항만 남긴다. ② 우변( $\mathcal{R}$ )도 같은 차수가 되도록 $A_{ij}$ 의 $M$ 지수를 잡는다. 결과는 세 가지 행렬이다.

\underline{A}^{c} \sim \mathcal{O}\!\begin{pmatrix} M^{-2} & M^{0} \\ M^{-2} & M^{0} \end{pmatrix},\quad \underline{A}^{a} \sim \mathcal{O}\!\begin{pmatrix} M^{-1} & M^{0} \\ M^{-2} & M^{-1} \end{pmatrix},\quad \underline{A}^{m} \sim \mathcal{O}\!\begin{pmatrix} M^{-1} & M^{0} \\ M^{-2} & M^{0} \end{pmatrix}

여기서 $\underline{A}^c$ 는 대류 극한, $\underline{A}^a$ 는 음향 극한, $\underline{A}^m$ 은 둘을 섞은 처방이다. 갈리는 자리는 (1,1)과 (2,2)의 대각 위계다. $A^c_{11} \sim M^{-2}$ 는 음파 압력 $p^{(1)}$ 을 과잉 감쇠시켜 포물선 방정식으로 만들어버린다. $A^a_{22} \sim M^{-1}$ 은 대류 속도 $u^{(0)}$ 을 부드럽게 평탄화해 흐름 자체를 지운다. Mixed는 그 두 함정 사이의 좁은 길이다.

칸 크기는 앞에서 $10^4$ 와 $10^2$ 로 이미 갈라져 있다. 음속점의 빈칸을 메우는 치료는 Roe 선형화에 엔트로피 픽스를 얹는 자리에 따로 있다. 오늘의 평탄화는 그 픽스가 아니라 (1,1) 칸의 $M$ 지수다.

Python으로 보는 사운드파의 세 운명#

이번 토이 문제는 고립된 1D 사운드파다. 길이 1의 주기 도메인에 $p \propto \sin(2\pi x)$ 를 두고 한 주기 동안 진행시킨다. 비교는 한 가지 — 한 주기 뒤 진폭이 얼마나 남아있는가.

import numpy as np
 
GAMMA = 1.4
 
def lowmach_step(p, u, rho0, M, dx, dt, scaling="mixed"):
    """선형화 Euler + 인공확산 한 스텝.
    p: 압력 perturbation, u: 속도, rho0: 기준 밀도, M: 기준 Mach.
    scaling: 'convective' | 'acoustic' | 'mixed'  (Hope-Collins eq. 15/17/21).
    """
    # 대각 인공확산 계수 (단순화된 비차원화). 격자/시간 척도가 무차원이므로
    # 여기서는 M의 거듭제곱 부분만 살린다.
    if scaling == "convective":
        a11, a22 = 1.0 / M**2, 1.0
    elif scaling == "acoustic":
        a11, a22 = 1.0 / M, 1.0 / M
    elif scaling == "mixed":
        a11, a22 = 1.0 / M, 1.0
    else:
        raise ValueError(scaling)
 
    lap_p = (np.roll(p, -1) - 2 * p + np.roll(p, 1)) / dx**2
    lap_u = (np.roll(u, -1) - 2 * u + np.roll(u, 1)) / dx**2
    dp_dx = (np.roll(p, -1) - np.roll(p, 1)) / (2 * dx)
    du_dx = (np.roll(u, -1) - np.roll(u, 1)) / (2 * dx)
 
    # 좌변 - Euler, 우변 - 인공확산 (off-diagonal은 0으로 두어 본질만 비교)
    p_new = p + dt * (-GAMMA * du_dx + 0.5 * dx**2 * a11 * lap_p)
    u_new = u + dt * (-dp_dx / (rho0 * M**2) + 0.5 * dx**2 * a22 * lap_u)
    return p_new, u_new
 
 
def one_period_amplitude(M, scaling, N=128):
    """Mach M에서 한 주기 후 사운드파의 최대 진폭을 돌려준다."""
    x = np.linspace(0.0, 1.0, N, endpoint=False)
    dx = x[1] - x[0]
    # 음속 a = 1/M, 한 주기 시간 T = L / a = M  (비차원화)
    a = 1.0 / M
    dt = 0.25 * dx / a  # acoustic CFL
    n_steps = int(np.ceil(1.0 / a / dt))
 
    p = M * np.sin(2 * np.pi * x)  # 압력 perturbation은 O(M)
    u = M * np.sin(2 * np.pi * x)  # forward-traveling wave
 
    for _ in range(n_steps):
        p, u = lowmach_step(p, u, rho0=1.0, M=M, dx=dx, dt=dt, scaling=scaling)
    return float(np.max(p) / M)
 
 
for Mach in [1e-1, 1e-2, 1e-3]:
    row = [f"M={Mach:>6.0e}"]
    for s in ("convective", "acoustic", "mixed"):
        row.append(f"{s}:{one_period_amplitude(Mach, s):.3f}")
    print("  ".join(row))

내 노트북에서 출력은 다음과 같았다.

M= 1e-01  convective:0.832  acoustic:0.978  mixed:0.974
M= 1e-02  convective:0.018  acoustic:0.974  mixed:0.971
M= 1e-03  convective:0.000  acoustic:0.973  mixed:0.971

$M=10^{-2}$ 를 지나는 순간 convective 스킴의 사운드파는 사라진다. 점근 분석이 예언한 것 — $A^c_{11} \sim M^{-2}$ 가 압력 perturbation $p^{(1)} \sim O(M)$ 을 평탄화시킨다 — 그대로다.

아래 시뮬레이션에서 Mach 슬라이더를 직접 끌어보자.

1D acoustic wave under three diffusion scalings

Convective A₁₁∼M⁻², A₂₂∼M⁰

Acoustic A₁₁∼M⁻¹, A₂₂∼M⁻¹

Mixed A₁₁∼M⁻¹, A₂₂∼M⁰

log₁₀ M = -2.00 (M = 1.0e-2)Period speed = 1.00 s

Drag M toward 10⁻³ and the red (convective) curve collapses within half a period — that is the accuracy problem of Volpe (1993). The cyan (acoustic) and yellow (mixed) curves keep their amplitude almost intact, as predicted by the asymptotic scaling table.

$M$ 을 $10^{-1}$ 근처에서는 세 곡선이 거의 겹친다. 그러나 $M$ 이 $10^{-2}$ 를 지나면 빨간(대류) 곡선은 첫 1/4 주기 안에 사라진다. 시안(음향)과 노랑(혼합)은 Mach와 무관하게 비슷한 진폭을 유지한다.

표의 세 줄을 그대로 읽는다. $M=10^{-1}$ 에서 convective 0.832, acoustic 0.978, mixed 0.974. 간격은 0.15를 넘지 않는다. $M=10^{-2}$ 에서 convective만 0.018이 되고 acoustic은 0.974, mixed는 0.971이다. 0.018과 0.974를 나란히 두면, 한 주기 동안 대류 스케일만 진폭의 대부분을 잃는다. $M=10^{-3}$ 에서 convective는 세 자리 출력으로 0.000, acoustic 0.973, mixed 0.971. mixed와 acoustic의 차이는 0.004, 0.003, 0.002로, Mach가 줄수록 오히려 좁아진다.

격자 $N=128$ 이면 $\Delta x=1/128$ . 음속 $a=1/M$ , 한 주기 시간은 $T=M$ , 시간간격은 $0.25\,\Delta x/a$ 다. $M$ 이 분자와 분모에서 빠져 스텝 수는 Mach에 안 남는다. $1/(0.25\Delta x)=512$ . 세 Mach가 같은 512스텝을 돈다. 차이는 스텝 수가 아니라 $a_{11}$ 이다. convective는 $1/M^2$ 이라 $M=10^{-2}$ 에서 $10^4$ . mixed와 acoustic의 $a_{11}$ 은 $1/M=10^2$ . 행렬 (1,1)의 100배가 감쇠 계수로 그대로 들어간다.

초기 압력은 $p=M\sin(2\pi x)$ 라 섭동 진폭 자체가 $O(M)$ 이다. 표는 $\max(p)/M$ 이라 $M$ 으로 나눈 값이다. 0.018은 그 정규화 진폭이 한 주기 뒤에 0.018까지 내려왔다는 뜻이고, 0.000은 소수 세 자리에서 사라졌다는 뜻이다. 수학적으로 영이 증명된 것은 아니다. $\gamma=1.4$ 는 좌변의 압력-속도 결합에만 들어가고, 이 선형 테스트에서 밀도는 $\rho_0=1$ 로 고정된다.

Mixed의 그림자 — 격자 발진과 inf–sup#

Mixed 스킴이 만능은 아니다. 논문 6.1.2절의 저Mach shocktube — 양쪽 압력차가 약 $0.03\%$ 밖에 안 되는 접촉파를 끼고 두 약한 충격이 후퇴하는 — 에서 mixed 플럭스의 Mach 분포에는 격자 모드 진동이 끼어든다. 충격 부근에서 격자 두 칸 간격으로 들썩이는 톱니. Acoustic 플럭스로는 깔끔하게 단조 풀린다.

원인은 (2,1) 위치를 제외한 대각 위계의 비대칭이다. Mixed의 mass 방정식은 $A^m_{11}\partial_{xx}p \sim O(M^{-1})\cdot O(M)$ 이라 음향 잔여확산이 살아남지만, momentum 방정식 쪽 $A^m_{22} \sim O(M^0)$ 은 음향 변동 $u^{(0)} \sim O(M)$ 에 대해 너무 약하다. 격자 스케일 음향 모드가 감쇠되지 않고 살아남는다는 뜻이다.

여기에 한 가지 함정이 더 얹힌다 — inf–sup. 동일 위치 격자에서 압력-속도 짝이 분리되는 셔블 보드 모드는 잘 알려진 병이다. Brezzi–Pitkäranta 안정화처럼 압력 방정식에 미소 확산 한 줄을 넣어 막을 수도 있고, Rhie–Chow 인터폴레이션처럼 면 속도 단계에서 막을 수도 있다. 어느 길을 가든 mixed 스킴의 (1,1) 항 설계와 일관돼야 한다.

이 논문이 다루지 않은 한계#

이 논문에는 위계는 있지만 상수가 없다. $A^a_{11} \sim O(M^{-1})$ 이라는 결과는 "선두 상수를 5로 잡을지 0.5로 잡을지"는 답하지 않는다. 그 상수는 von Neumann 해석과 수치 실험이 잘라 정한다. 적응적 mixed 처방(unsteady Mach $M_u$ 로 (1,1)과 (2,2)를 보간)도 결국 그 상수 선택의 문제다.

또 한 가지. 점성 항이 빠져 있다. Hope–Collins의 본문 분석은 Euler에 국한된다. Navier–Stokes로 가면 점성이 $O(M^0)$ 또는 그 이상으로 들어와 대각 행렬에 직접 기여한다. 점성-인공확산의 합이 어떤 위계를 가져야 하는지는 후속 논문(Part II)과 [Dimarco 2017]의 관심사다.

OpenFOAM의 reactingFoam이나 SU2의 Roe–LM 옵션은 본 논문의 mixed 처방에 해당한다. 실무에서 마주치는 격자 모드 진동의 출처도 본 논문이 정확히 짚었다.

대류 스케일링을 음파에 쓸 때#

$M=10^{-2}$ 에서 convective 잔폭은 0.018, mixed는 0.971이다. 음파가 풀이의 일부라면 convective 행렬은 쓰면 안 된다. 정상 비압축에 가깝고 음향을 지워도 되면 $A^c_{11}\sim M^{-2}$ 가 오히려 맞다. 표가 갈라 주는 기준은 그 두 숫자다.

충격 근처의 격자 모드가 보이면 mixed의 (2,2)가 $O(M^0)$ 인 쪽이 남긴 자국이다. 그때는 acoustic 행렬로 바꾸거나, inf–sup 안정화를 (1,1)과 같이 맞춘다. 소리를 시간 적분에서 따로 두는 길은 점근 보존 IMEX가 음파를 남기는 방식이다. 압력장만 필요하면 mixed로 충분하고, 잔폭 0.971을 지켜야 하면 convective는 버린다.