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

저마하(low Mach) 유동을 압축성 코드에 넣으면, 음속이 빠를수록 그림이 선명해질 것 같다. 반대다. 마하수가 작아질수록 표준 풍상(upwind)은 멀쩡한 파동을 격자에서 지운다. Arun·Das Gupta·Samantaray(2019)는 선형 파동 방정식계로 그 붕괴를 짚고, IMEX-RK로 마하수와 무관하게 안정한 점근 보존(AP) 스킴을 세웠다. 음향 연산자를 왜 음함수로 풀어야 하는지는 아래 숫자로 확인된다.

ε 슬라이더를 왼쪽으로 밀면 빨간 풍상 곡선은 기울기 −1로 치솟고, 시안 AP 곡선은 평평하게 남는다.

$\varepsilon=10^{-3}$에서 두 곡선의 비는 수백 배까지 벌어진다. 마하수를 줄일수록 커지는 것은 파동이 아니라 풍상의 수치 점성이다.

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

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

$\rho$는 밀도, $\mathbf{u}$는 속도, $p=\rho^\gamma$는 압력이다. $\varepsilon$은 기준 속도 나누기 기준 음속이다.

$\varepsilon \to 0$이면 해는 비압축 극한으로 간다. 특이 섭동(singular perturbation)이다. 저자들은 분석을 줄이려고 선형 파동계를 쓴다.

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

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

극한 $\varepsilon\to 0$에서 해는 아무 데로나 새지 않는다. 잘 준비된(well-prepared) 공간으로 모인다. 밀도는 상수, 속도의 발산은 0이다.

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

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

범인은 풍상의 수치 점성이다. 파동계를 특성 변수로 대각화하면 속도 $\pm c$의 스칼라 이류 둘이 된다. 각 모드에 Rusanov 풍상을 씌우면 한 스텝의 증폭이 나온다.

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

한 스텝의 손실은 작다. 쌓이는 쪽이 문제다. 물리 시간 $T$까지 가려면 음향 스텝이 $M = T/\Delta t = T c /(\nu\Delta x)$번 필요하다. $\theta$가 작으면 한 스텝 감쇠는 $\tfrac{1}{2}\nu(1-\nu)\theta^2$이고, $M$번 곱하면 로그 감쇠가 이렇게 된다.

마하수를 한 자릿수 내리면 실효 점성도 한 자릿수 커진다. 앞의 그림에서 빨간 곡선의 기울기 −1이 이 비례다.

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

방정식의 시간 척도를 갈라 보면 연산자가 둘이다.

$H$는 시간 척도 $O(1)$의 대류다. $L/\varepsilon$은 시간 척도 $O(\varepsilon)$의 음향이다. 뻣뻣함(stiffness)은 전부 $L/\varepsilon$에 있다. 그 부분만 음함수로 푼다. 대류는 양함수로 둔다. 이 혼합이 IMEX(Implicit-Explicit)다. 같은 분리를 마하 0까지 가져갔을 때 시간간격이 버티는지는 IMEX TVD 재현이 이어서 다룬다.

음향을 음함수로 두면 중심 차분도 안정하다. 풍상 점성이 빠진다. $\Delta t$는 빠른 음향이 아니라 느린 대류가 정한다. $\varepsilon$이 작아도 간격을 줄일 이유가 이 모델 안에서는 없다.

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

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

$\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

Mach $\varepsilon$를 내리면 빨간 진폭이 물리 시간 $T=1$에 닿기 전에 바닥으로 간다. CFL $\nu$를 1 쪽으로 밀면 $\nu(1-\nu)$가 줄어 감쇠가 한숨 꺾인다. $\varepsilon$을 더 내리면 그 이득은 $1/\varepsilon$ 누적에 묻힌다. $N$을 키워도 $\theta^2\propto N^{-2}$와 $\Delta x\propto N^{-1}$가 일부만 상쇄된다.

코드의 수를 식에 넣으면 스텝 수가 닫힌 형태로 나온다. $a=1$, $T=1$, $\nu=0.5$, $N=64$이면 $\Delta x=1/64$이고 $c=1/\varepsilon$이다. $\Delta t=\nu\Delta x/c$이므로 $M=T/\Delta t=128/\varepsilon$이다. 표의 320, 640, 1280, 2560, 6400, 12800은 $\varepsilon=0.4,0.2,0.1,0.05,0.02,0.01$을 이 식에 넣은 값이다. $\varepsilon$을 반으로 줄이면 $M$은 두 배다. 열 배로 줄이면 $M$도 열 배다.

로그 감쇠 $\tfrac12 T c(1-\nu)k^2\Delta x$에 같은 수를 넣으면 지수는 $0.154/\varepsilon$이다. $k=2\pi$라서 계수가 그렇게 나온다. $\varepsilon=0.4$에서 지수는 0.386이고 $\mathrm{e}^{-0.386}=0.680$이다. 출력의 $6.804\times 10^{-1}$과 소수점 세 자리까지 같다. $\varepsilon=0.1$에서는 $\mathrm{e}^{-1.542}=0.214$로, 출력 $2.143\times 10^{-1}$과 맞는다. $\varepsilon=0.01$에서 지수는 15.4다. $\mathrm{e}^{-15.4}$는 $2.0\times 10^{-7}$ 근처이고, 출력은 $1.980\times 10^{-7}$이다. $\theta=k\Delta x=2\pi/64\approx 0.098$이라 작은 각 가정이 맨 아래까지 정확하진 않다. 표의 끝이 근사보다 조금 더 죽는 이유다. $\nu=0.5$는 $\nu(1-\nu)=0.25$로, 이 곱이 가장 큰 CFL이다. 그래서 이 열은 한 스텝 감쇠가 제일 센 조건의 기록이다.

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

둘은 다르다. AP는 이산 스킴의 $\varepsilon\to 0$ 극한이 극한 방정식의 올바른 이산화인지 묻는다. 안정 제한이 $\varepsilon$과 무관해야 한다. AA는 한 걸음 더다. 격자는 유한하고 $\varepsilon$은 작지만 0이 아닐 때, 해가 well-prepared 공간 근처에 머무는가. 풍상은 둘 다 실패한다. IMEX 음함수는 둘 다 통과한다. 갈리는 자리는 공간 차분의 차수가 아니라 시간 이산화다.

ν=0.5의 최악 감쇠에서도 AP 열은 1.00이다#

이 표의 AP 열은 여섯 줄 모두 1.00이다. 음함수가 풍상 점성을 빼서, 지수 $0.154/\varepsilon$이 진폭에 붙지 않기 때문이다. $\nu=0.5$는 한 스텝 감쇠 인자가 가장 큰데도 $M=128/\varepsilon$번을 쌓은 뒤 1.00이다. 이 선형 파동계에서는 $L/\varepsilon$만 음함수로 두는 선택이 맞다.

그 선택이 저마하 방정식 전체로 그대로 옮지는 않는다. 음향을 음함수로 넣었는데 시간간격이 다시 크게 줄어드는 계가 있다. 상대속도에 두 번째 $1/M^2$가 남는 자리는 엔탈피 쪽에 숨은 강성에 적혀 있다.