[논문리뷰] 저마하에서 소리를 죽이지 않는 법 — 점근 보존 IMEX 파동 스킴

저마하 극한에서 풍상 점성이 해를 죽이는 이유와 AP-IMEX 처방

저마하(low Mach) 유동을 압축성 코드로 풀면, 음속이 빨라질수록 해가 더 또렷해질 것 같다. 실제로는 정반대다. 마하수가 작아질수록 표준 풍상(upwind) 스킴은 멀쩡한 파동을 격자에서 지워버린다. 이 글은 Arun·Das Gupta·Samantaray(2019)가 선형 파동 방정식계를 모델로 그 붕괴의 원인을 분석하고, IMEX-RK 시간 적분으로 마하수와 무관하게 안정한 점근 보존(AP) 스킴을 세우는 과정을 따라간다. 읽고 나면 "왜 음향 연산자를 음함수로 풀어야 저마하에서 정확한가"를 코드로 직접 확인할 수 있다.

ε가 0으로 갈 때 격자가 파동을 잡아먹는다#

스케일을 맞춘 등엔트로피 Euler 방정식에는 작은 매개변수 하나가 박혀 있다. 마하수 $\varepsilon$ 이다. 운동량 방정식의 압력항이 $\nabla p / \varepsilon^2$ 로 들어온다.

\partial_t \rho + \nabla\cdot(\rho \mathbf{u}) = 0, \qquad \partial_t (\rho\mathbf{u}) + \nabla\cdot(\rho\mathbf{u}\otimes\mathbf{u}) + \frac{\nabla p}{\varepsilon^2} = 0

여기서 $\rho$ 는 밀도, $\mathbf{u}$ 는 속도, $p=\rho^\gamma$ 는 압력, $\varepsilon$ 은 기준 속도/기준 음속의 비(마하수)다.

$\varepsilon \to 0$ 이면 해는 비압축성 극한으로 수렴한다. 수학적으로는 특이 섭동(singular perturbation) 문제다. 분석을 단순화하려고 저자들은 선형 파동 방정식계를 쓴다.

\partial_t \varrho + (\mathbf{u}\cdot\nabla)\varrho + \frac{a}{\varepsilon}\nabla\cdot\mathbf{u} = 0, \qquad \partial_t \mathbf{u} + (\mathbf{u}\cdot\nabla)\mathbf{u} + \frac{a}{\varepsilon}\nabla\varrho = 0

여기서 $\varrho$ 는 선형화된 밀도 섭동, $a$ 는 선형화 상태의 음속이다. $a/\varepsilon$ 이 실효 음속 $c$ 다. 마하수가 작아지면 이 음향 속도가 폭발한다.

well-prepared 공간 — 살아남아야 할 해#

극한 $\varepsilon\to 0$ 에서 해는 아무 데로나 가지 않는다. "잘 준비된(well-prepared)" 공간으로 수렴한다. 밀도가 상수이고 속도가 발산이 0인 장들의 집합이다.

\varrho^{(0)} = \text{const}, \qquad \nabla\cdot\mathbf{u}^{(0)} = 0

이 공간은 음향 연산자 $L(U)=(a\nabla\cdot\mathbf{u},\,a\nabla\varrho)$ 의 핵(kernel)이다. 좋은 스킴이라면 이산 해도 같은 공간으로 수렴해야 한다. Dellacherie(2010)는 이 불변성을 깨는 스킴이 저마하에서 정확도를 잃는다는 것을 보였다. 핵을 보존하면 점근 정확(AA), 안정 제약이 $\varepsilon$ 과 무관하면 점근 보존(AP)이다.

풍상 점성은 음속에 비례해 커진다#

문제의 정체는 풍상 스킴의 수치 점성이다. 파동계를 특성 변수로 대각화하면 속도 $\pm c$ 로 흐르는 두 스칼라 이류로 분리된다. 각 모드에 풍상(Rusanov) 스킴을 적용하면 한 스텝당 증폭 계수가 나온다.

|g|^2 = \big(1 - \nu(1-\cos\theta)\big)^2 + \nu^2\sin^2\theta, \qquad \nu = \frac{c\,\Delta t}{\Delta x},\ \ \theta = k\Delta x

여기서 $\nu$ 는 음향 CFL 수, $\theta$ 는 셀당 위상, $k$ 는 파수다.

핵심은 누적이다. 고정된 물리 시간 $T$ 에 도달하려면 음향 스텝 수가 $M = T/\Delta t = T c /(\nu\Delta x)$ 만큼 필요하다. 작은 $\theta$ 에서 한 스텝의 감쇠는 $\tfrac{1}{2}\nu(1-\nu)\theta^2$ 이고, $M$ 번 쌓으면 로그 감쇠가 이렇게 된다.

M\cdot\tfrac{1}{2}\nu(1-\nu)\theta^2 \;=\; \tfrac{1}{2}\,T\,c\,(1-\nu)\,k^2\Delta x \;\propto\; c \;\propto\; \frac{1}{\varepsilon}

마하수가 한 자릿수 작아질 때마다 실효 점성이 한 자릿수 커진다. 아래 그래프에서 직접 끌어보자.

Mach ε = 0.050

D_explicit ≈ 0.078

D_AP ≈ 7.8e-3

ratio ≈ 10.0×

$\varepsilon$ 을 왼쪽으로 옮기면 풍상 곡선(빨강)은 기울기 $-1$ 로 치솟고 AP 곡선(시안)은 평평하게 머문다. $\varepsilon=10^{-3}$ 에서 둘의 비가 수백 배로 벌어진다. 이것이 저마하 정확도 붕괴의 정량적 모습이다.

연산자 분리 — 대류는 양함수로, 음향은 음함수로#

처방의 출발점은 방정식을 시간 척도로 쪼개는 것이다. 진화형으로 다시 쓰면 두 연산자가 보인다.

\partial_t U + H(U) + \frac{1}{\varepsilon}L(U) = 0

여기서 $H$ 는 시간 척도 $O(1)$ 의 대류 연산자, $L/\varepsilon$ 은 시간 척도 $O(\varepsilon)$ 의 음향 연산자다. 뻣뻣함(stiffness)은 전부 $L/\varepsilon$ 에서 나온다. 그렇다면 이 부분만 음함수로 풀면 된다. 대류는 비싸지 않으니 양함수로 둔다. 이 양함수-음함수 혼합이 IMEX(Implicit-Explicit)다.

음향을 음함수로 처리하면 중심 차분을 써도 안정하다. 풍상 점성이 사라진다. 시간 간격은 빠른 음향이 아니라 느린 대류 척도가 정한다. $\varepsilon$ 이 작아도 $\Delta t$ 를 줄일 필요가 없다.

IMEX-RK가 뻣뻣함을 받아낸다#

저자들은 시간 적분으로 IMEX 룽게-쿠타(IMEX-RK)를 쓴다. 명시적 표 $(\tilde A,\tilde b)$ 가 $H$ 를, 암시적 표 $(A,b)$ 가 $L/\varepsilon$ 을 맡는다. 각 스테이지에서 음향 부분이 밀도에 대한 타원형 방정식(압력 포아송과 닮은 꼴)으로 응축된다.

\varrho - \beta^2 \Delta t^2\, a^2 \nabla^2 \varrho = (\text{명시적 항})

여기서 $\beta$ 는 IMEX 표의 대각 계수, $\nabla^2$ 은 라플라시안이다. 이 타원형 문제의 가역성은 변분 안장점(saddle point) 이론으로 보장된다. 이산 단계에서는 순환 행렬(circulant matrix) 이론으로 같은 결과를 얻는다. 결과는 세 가지다. 유일해의 존재, $\varepsilon$ 에 무관한 균일 안정, 그리고 well-prepared 공간의 불변성.

Python — 한 모드의 진폭 감쇠를 재현한다#

아래 코드는 단일 푸리에 모드의 풍상 증폭 계수를 계산하고, 고정 물리 시간 $T$ 뒤에 남는 진폭을 마하수별로 출력한다. 음함수 스킴은 모드를 보존하므로 진폭이 1로 남는다.

import numpy as np
 
a  = 1.0        # 선형화 음속
k  = 2 * np.pi  # [0,1] 위의 단일 푸리에 모드
T  = 1.0        # 고정 물리 종료 시간
nu = 0.5        # 음향 CFL 수
 
def upwind_amp_factor(eps, N):
    """파동계 Rusanov 한 특성 모드의 증폭 |g|."""
    c  = a / eps              # 음향 속도 c = a/ε
    dx = 1.0 / N
    th = k * dx               # 셀당 위상 θ = kΔx
    re = 1.0 - nu * (1.0 - np.cos(th))
    im = nu * np.sin(th)
    return np.hypot(re, im)   # |g| <= 1
 
def retained_amplitude(eps, N):
    """물리 시간 T 뒤 모드 진폭 (풍상 vs 음함수)."""
    c  = a / eps
    dx = 1.0 / N
    dt = nu * dx / c          # 음향 CFL 스텝
    M  = T / dt               # 명시적 스텝 수
    explicit = upwind_amp_factor(eps, N) ** M
    implicit = 1.0            # 중심-음함수는 모드를 보존
    return M, explicit, implicit
 
if __name__ == "__main__":
    N = 64
    print(f"{'eps':>8} {'steps M':>10} {'explicit':>12} {'AP':>6}")
    for eps in [0.4, 0.2, 0.1, 0.05, 0.02, 0.01]:
        M, ex, ap = retained_amplitude(eps, N)
        print(f"{eps:8.3f} {M:10.0f} {ex:12.3e} {ap:6.2f}")

출력은 이렇게 나온다.

     eps    steps M     explicit     AP
   0.400        320    6.804e-01   1.00
   0.200        640    4.629e-01   1.00
   0.100       1280    2.143e-01   1.00
   0.050       2560    4.591e-02   1.00
   0.020       6400    4.430e-04   1.00
   0.010      12800    1.980e-07   1.00

$\varepsilon$ 이 10배 작아질 때마다 스텝 수가 2배가 되고, 풍상 진폭은 지수적으로 0으로 떨어진다. $\varepsilon=0.01$ 에서는 모드가 사실상 사라진다. 같은 격자에서 음함수 스킴은 100%를 지킨다.

아래 시뮬레이션에서 직접 조작해보자. Mach $\varepsilon$ 슬라이더를 내리면 빨간 곡선(풍상)이 물리 시간이 흐르기도 전에 바닥으로 꺼진다.

Mach ε: 0.050acoustic CFL ν: 0.50grid cells N: 64

CFL $\nu$ 를 1에 가깝게 올리면 $\nu(1-\nu)$ 가 작아져 감쇠가 잠깐 누그러진다. 하지만 $\varepsilon$ 을 더 줄이면 그 이득은 곧 $1/\varepsilon$ 누적에 묻힌다. 격자 $N$ 을 키워도 $\theta^2\propto N^{-2}$ 와 $\Delta x \propto N^{-1}$ 가 부분적으로만 상쇄돼 근본 문제는 남는다.

점근 보존(AP)과 점근 정확(AA)#

두 성질을 헷갈리기 쉽다. AP는 "이산 스킴의 $\varepsilon\to 0$ 극한이 극한 방정식의 올바른 이산화가 되는가"를 묻는다. 안정 제약이 $\varepsilon$ 과 무관해야 한다. AA는 한 발 더 나간다. 유한한 격자에서, 작지만 0이 아닌 $\varepsilon$ 에 대해 해가 well-prepared 공간 근처에 머무는가를 본다. 풍상 스킴은 둘 다 실패한다. IMEX-음함수 스킴은 둘 다 만족한다. 핵심은 시간 이산화의 선택이지 공간 차분의 정밀도가 아니다.

핵심 3줄#

저마하에서 풍상 점성은 음속 $c=a/\varepsilon$ 에 비례해 $1/\varepsilon$ 로 커지고, 고정 물리 시간에 누적돼 파동을 지운다.
well-prepared 공간(상수 밀도·발산 0 속도)을 보존하는지가 정확도의 갈림길이며, 그것을 좌우하는 건 시간 이산화다.
음향 연산자 $L/\varepsilon$ 만 IMEX로 음함수 처리하면 $\Delta t$ 가 대류 척도에서 풀려 마하수와 무관하게 안정·정확해진다(AP+AA).