[논문 리뷰] 비선형 플럭스를 선형으로 속이는 법

여유 있는 a는 안정이지만, Chapman–Enskog 항 ε(a²−f'²)이 그 여유만큼 충격을 뭉갠다.

정확한 리만 해법기(Riemann solver, 셀 경계에서 파동 구조를 푸는 장치)를 짜다 보면 상태방정식에서 멈춘다. 비선형 플럭스 $f(u)$ 의 야코비안을 매번 대각화해야 한다. 실기체면 고유구조가 닫힌 꼴로 안 나온다. 1995년 Shi Jin과 Zhouping Xin은 방향을 뒤집었다. 비선형 방정식을 그대로 밀지 않는다. 플럭스 자체를 새 미지수로 올리고, 선형 계를 만든 뒤 원래 식으로 되돌린다. 완화(relaxation)다. Thomann 외(2019)는 이 뼈대를 전속도(all-speed) Euler까지 밀었다. 우회가 어디서 단단하고 어디서 무너지는지. 그 대가가 이 노트의 전부다.

논문 정보 — Andrea Thomann, Markus Zenk, Gabriella Puppo, Christian Klingenberg, An all speed second order IMEX relaxation scheme for the Euler equations, arXiv:1907.08398 (2019).

완화 속도 슬라이더 a를 $1$ 밑으로 내리면, 왼쪽 $1$ ·오른쪽 $-0.4$ 인 Burgers 충격의 파면이 흔들린다.

relaxation speed a = 1.30relaxation time ε = 1e-3

a = 1.30 ≥ max|u| = 1.00 — sub-characteristic condition satisfied. Smaller ε projects v onto f(u) faster (sharper shock, more relaxation diffusion trade-off).

a를 $1.3$ 근처에 두면 파면이 선다. eps만 키우면 속도는 그대로 두고 두께만 는다.

비선형 플럭스가 리만 문제를 어렵게 만든다#

스칼라 보존법칙 $u_t + f(u)_x = 0$ 을 생각하자. $u$ 는 보존량, $f(u)$ 는 플럭스다. 격자 한 칸의 충격관을 경계마다 풀려면 파동 속도 $f'(u)$ 가 먼저 필요하다.

문제는 실제 유동이다. 압축성 Euler에서 플럭스는 압력을 품고, 압력은 밀도·내부에너지의 비선형 함수다. 실기체(real gas) 상태방정식이 끼면 고유값·고유벡터가 해석적으로 안 풀린다. 매 스텝 수치적으로 야코비안을 대각화하는 비용은 무겁다. 저Mach 영역에서는 음속 파동이 $1/M$ 로 커져 명시적 시간 스텝이 질식한다.

그래서 나온 질문. 비선형성을 매 스텝 정면으로 상대하지 않고 미룰 수는 없을까?

Jin과 Xin의 승격 — 플럭스를 미지수로 올린다#

Jin–Xin의 답은 이렇다. 플럭스 $f(u)$ 를 값으로 계산하는 대신, 그 값을 담을 **새 변수 $v$ **를 도입한다. 그리고 $v$ 가 천천히 $f(u)$ 로 끌려가도록 소스항을 붙인다.

\begin{aligned} u_t + v_x &= 0, \\ v_t + a^2 u_x &= \frac{1}{\varepsilon}\bigl(f(u) - v\bigr). \end{aligned}

여기서 $u$ 는 원래 보존량, $v$ 는 플럭스를 대신하는 완화 변수, $a$ 는 상수 완화 속도(relaxation speed), $\varepsilon$ 은 완화 시간이다. 좌변은 완전히 선형이다. 계수 행렬의 고유값은 $\pm a$ 로 고정—비선형 대각화가 사라졌다.

$\varepsilon \to 0$ 이면 오른쪽 소스가 $v = f(u)$ 를 강제하고, 첫 식은 원래 보존법칙으로 돌아간다. 즉 이 $2\times 2$ 선형 계는 비선형 방정식의 점성 근사다. 파동 속도가 상수 $\pm a$ 이므로 리만 문제는 단 한 번, 손으로 풀 수 있다. 위 그림이 이 계로 민 Burgers, $f(u)=u^2/2$ 다.

부분특성 조건 — a가 작으면 계가 폭발한다#

완화 속도 $a$ 를 아무렇게나 고를 수는 없다. 안정적인 근사가 되려면 부분특성 조건(sub-characteristic condition, Whitham 조건)을 지켜야 한다.

a \ge \max_u |f'(u)|.

$f'(u)$ 는 원래 방정식의 진짜 파동 속도다. 완화계가 들고 다니는 얼린 속도 $\pm a$ 가 그 진짜 속도를 항상 감싸야 한다는 뜻이다. Burgers라면 $f'(u)=u$ 이므로 $a \ge \max|u|$ 면 된다.

아래 그림에서 완화 속도 $a$ 와 해의 범위 $\max|u|$ 를 바꿔보자. 주황 곡선 $f'(u)$ 가 $\pm a$ 띠를 벗어나는 순간이 불안정의 문턱이다.

relaxation speed a = 1.40solution range max|u| = 1.00

The orange curve f′(u) stays inside the ±a band over the whole solution range → sub-characteristic condition holds.

$a$ 를 해의 최대 파동 속도 아래로 내리면 곡선이 띠를 뚫고 나간다. 이때 완화 근사는 물리적 확산이 아니라 음의 확산을 만든다. 정보가 잘못된 방향으로 흐르고 계는 폭발한다. 이유는 다음 절의 전개에서 정확히 드러난다.

Chapman–Enskog로 본 숨은 확산#

$\varepsilon$ 이 작지만 0은 아닐 때 무슨 일이 벌어질까. $v = f(u) + \varepsilon v_1 + \cdots$ 로 전개(Chapman–Enskog)해 완화계에 대입하면, 유효 방정식이 나온다.

u_t + f(u)_x = \varepsilon\,\partial_x\Bigl[\bigl(a^2 - f'(u)^2\bigr)\,u_x\Bigr] + O(\varepsilon^2).

오른쪽이 완화가 몰래 집어넣는 확산항이다. 확산계수는 $a^2 - f'(u)^2$ . 이 값이 음수가 아니려면 정확히 $a \ge |f'(u)|$ 가 필요하다. 부분특성 조건의 정체는 **"완화가 만드는 인공 점성이 음수가 되지 않을 조건"**이었던 것이다.

여기서 트레이드오프가 보인다. $a$ 를 크게 잡으면 안정성은 넉넉하지만 $a^2-f'^2$ 이 커져 해가 뭉개진다. $a$ 를 부분특성 한계에 바짝 붙이면 날카롭지만 위험하다. $\varepsilon$ 은 확산의 전체 크기를 조절한다. 앞의 시뮬레이션에서 $\varepsilon$ 을 키우면 충격파가 두꺼워지는 이유다.

여유 폭이 점성의 크기다. 코드의 리만 데이터는 왼쪽 $u=1$ , 오른쪽 $u=-0.4$ 다. 출발할 때 $\max|u|=1$ 이라 부분특성 하한도 1이다. $a=1.3$ 이면 $a^2=1.69$ 다. 왼쪽 $f'=1$ 에서 확산계수는 $1.69-1=0.69$ . 오른쪽 $f'=-0.4$ 에서는 $1.69-0.16=1.53$ 으로 더 크다. 같은 $\varepsilon=10^{-4}$ 를 곱해도 오른쪽이 더 잘 퍼진다. $a=0.8$ 이면 $a^2=0.64$ 라 왼쪽 계수가 $0.64-1=-0.36$ 이다. 부호가 뒤집힌다. 파면 정보가 뒤로 밀리고, 진폭은 초기값 1에 머물지 못한다.

Suliciu — 플럭스 전부가 아니라 압력만 완화한다#

Jin–Xin은 플럭스의 모든 성분을 완화한다. 우아하지만 확산이 과하다. 실전에서 더 인기 있는 것은 Suliciu형 완화다. Euler 방정식에서 진짜 비선형은 압력 하나이므로, 압력만 새 변수 $\pi$ 로 완화한다.

(\rho\pi)_t + \nabla\!\cdot(\rho\pi\,u) + a^2\,\nabla\!\cdot u = \frac{\rho}{\varepsilon}\,(p - \pi).

$\rho$ 는 밀도, $u$ 는 속도, $p$ 는 실제 압력, $\pi$ 는 완화 압력이다. 부분특성 조건은 $a > \rho\sqrt{\partial_\rho p}$ 로, 음속을 감싸는 형태가 된다. 완화된 계의 고유값은 $u,\ u\pm a/\rho$ 로 모두 선형퇴화(linearly degenerate)—접촉 불연속처럼 다루기 쉬운 파동만 남는다.

Thomann 외(2019)의 기여가 여기서 시작된다. 이들은 저Mach 영역을 겨냥해 압력을 느린 성분과 빠른 음향 성분으로 쪼갠다. 느린 부분은 명시적으로, 빠른 음향 부분은 완화계 위에서 음함수(implicit)로 푼다. 새 속도 변수 $\hat u$ 를 추가해 Mach 수에 무관한 확산을 확보하고, 그 결과 저Mach 극한에서 비압축 Euler로 수렴하는 점근보존(asymptotic-preserving) 2차 IMEX 기법을 얻는다. 선형퇴화 구조 덕에 이 모든 것이 비선형 대각화 없이 굴러간다.

Python — Burgers를 완화계로 푼다#

완화 아이디어의 뼈대만 numpy로 재현해보자. Burgers 방정식을 Jin–Xin 완화계로 푼다. 특성 변수 $r,s$ 로 분해하면 좌변은 두 개의 단순 이류가 된다.

import numpy as np
 
def burgers_flux(u):
    return 0.5 * u * u                     # f(u) = u^2 / 2
 
def relaxed_step(u, v, a, dx, dt, eps):
    # 특성 변수 분해: r은 +a, s는 -a 속도로 이동
    r = 0.5 * (u + v / a)
    s = 0.5 * (u - v / a)
    nu = a * dt / dx                       # CFL 수 a*dt/dx
    # 1차 업윈드 (주기 경계, np.roll)
    r_new = r - nu * (r - np.roll(r, 1))   # 오른쪽 이동파
    s_new = s + nu * (np.roll(s, -1) - s)  # 왼쪽 이동파
    u_new = r_new + s_new
    v_new = a * (r_new - s_new)
    # 완화 소스: v를 f(u)로 끌어당김 (지수 적분)
    kappa = 1.0 - np.exp(-dt / eps)
    v_new += (burgers_flux(u_new) - v_new) * kappa
    return u_new, v_new
 
def run_relaxation(u0, a, eps, cfl=0.9, t_end=0.3):
    n = u0.size
    dx = 2.0 / n
    u = u0.copy()
    v = burgers_flux(u)                    # 평형 상태에서 출발
    dt = cfl * dx / a                      # 부분특성이 CFL을 정한다
    t = 0.0
    while t < t_end:
        u, v = relaxed_step(u, v, a, dx, dt, eps)
        t += dt
    return u
 
# Riemann 초기조건: 왼쪽 1.0, 오른쪽 -0.4 (충격파)
n = 400
x = np.linspace(-1, 1, n, endpoint=False) + 1.0 / n
u0 = np.where(x < 0.0, 1.0, -0.4)
 
u_ok = run_relaxation(u0, a=1.3, eps=1e-4)   # a >= max|u|=1.0  -> 안정
u_bad = run_relaxation(u0, a=0.8, eps=1e-4)  # a <  max|u|      -> 위반
 
print("a=1.3  peak |u| =", round(float(np.max(np.abs(u_ok))), 3))   # ~1.0
print("a=0.8  peak |u| =", round(float(np.max(np.abs(u_bad))), 3))  # 폭발

출력은 부분특성 조건을 그대로 증언한다. $a=1.3$ 이면 최대 진폭이 초기값 근처에 머물며 매끈한 충격파가 잡힌다. $a=0.8$ 이면 진폭이 몇 배로 튀어 오른다. 비선형 야코비안을 단 한 번도 대각화하지 않았다. 상수 속도 $\pm a$ 의 선형 이류만 반복했다.

걸음 수도 $a$ 가 정한다. $n=400$ , 구간 길이 2라 $\Delta x=0.005$ 다. CFL 0.9이면 $\Delta t=0.9\times 0.005/a$ 다. $a=1.3$ 에서 $\Delta t\approx 3.46\times 10^{-3}$ 이라 $t=0.3$ 까지의 비율은 약 87이다. $a=0.8$ 이면 $\Delta t\approx 5.63\times 10^{-3}$ 으로 성겨진다. $a$ 를 낮추면 시간 간격이 오히려 늘다. $\varepsilon=10^{-4}$ 이면 $\Delta t/\varepsilon$ 가 30을 넘어, $\kappa=1-\exp(-\Delta t/\varepsilon)$ 의 지수항이 사라진다. 매 걸음 $v$ 는 $f(u)$ 에 붙고, 남는 번짐은 $\varepsilon(a^2-f'^2)$ 뿐이다.

비판적으로 보면#

완화는 공짜가 아니다. 인공 확산 $\varepsilon(a^2-f'^2)$ 이 반드시 따라붙고, $a$ 를 안전하게 크게 잡을수록 해가 뭉개진다. 부분특성 조건을 지키려면 $\max|f'(u)|$ 의 전역 상한이 필요한데, 강한 충격이나 진공 근처에서 이 상한을 보수적으로 잡으면 확산이 과해진다. 또 스칼라에서는 깔끔하던 논리가 실기체 Euler에서는 $\rho\sqrt{\partial_\rho p}$ 의 국소 추정과 양의 압력·밀도 보존을 위한 추가 장치를 요구한다. Thomann 외의 IMEX 결합은 음향 부분에서 타원형 방정식을 풀어야 하므로, 저Mach 이득이 선형 솔버 비용과 상쇄되는 지점도 존재한다. 재현해보면 $\varepsilon$ 과 $a$ 의 동시 튜닝이 생각보다 예민하다.

OpenFOAM이나 Fluent의 밀도기반 솔버에는 이 완화 리만 해법기가 직접 들어 있지 않지만, HLLC·AUSM 계열의 근사 리만 해법기가 같은 철학—파동 구조를 단순화해 대각화를 피한다—을 공유한다. Suliciu 해법기는 SU2 같은 오픈소스 코드에 접촉 불연속을 정확히 잡는 옵션으로 구현되어 있다.

하한 바로 위에 고정한 a가 깨지는 자리#

이 데이터의 하한은 1이고, 코드는 1.3을 골랐다. 차이 0.3이 왼쪽 확산계수 0.69를 만든다. 여유를 0까지 깎아 $a=\max|u|$ 로 두면 파면 위 계수가 0이 된다. 포착은 날카롭다. 반올림이 $a$ 를 한 끝만 밑으로 밀면 부호가 뒤집힌다.

상한을 처음에만 재면 다음 충격에서 깨진다. $|u|$ 가 1을 넘는 순간, 안전해 보이던 $a=1.3$ 도 부분특성 밖이다. Suliciu로 옮기면 부등식이 $a>\rho\sqrt{\partial_\rho p}$ 로 바뀐다. 음속이 셀마다 튀는데 상한이 한 박자 늦으면, 음향을 음함수로 풀어도 완화 확산은 이미 음수다. 압력만 건드리겠다는 선택은 이류와 압력을 갈라 놓은 플럭스 분리와 같은 자리다.

[논문 리뷰] 비선형 플럭스를 선형으로 속이는 법 — Jin–Xin·Suliciu 완화 기법