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

완화계로 리만 문제를 선형화하고 부분특성 조건으로 안정성을 사는 Thomann(2019) 리뷰

정확한 리만 해법기(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).

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

스칼라 보존법칙 $u_t + f(u)_x = 0$ 을 생각하자. $u$ 는 보존량, $f(u)$ 는 플럭스다. Godunov류 기법은 셀 경계마다 리만 문제를 풀어 파동 속도 $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$ )을 위 완화계로 푼 결과다.

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).

shock 초기조건에서 완화 속도 a를 1.0 아래로 내리면 프로필이 진동하며 발산한다. a를 충분히 크게 두면 매끈한 충격파가 잡힌다. 왜 그런지가 다음 절이다.

부분특성 조건 — 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$ 을 키우면 충격파가 두꺼워지는 이유다.

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$ 의 선형 이류만 반복했다.

비판적으로 보면#

완화는 공짜가 아니다. 인공 확산 $\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 \ge |f'(u)|$ —부분특성 조건이 곧 인공 점성을 음수로 만들지 않는 조건이다.
선형퇴화 구조가 실전을 연다. Suliciu형 완화는 저Mach·실기체·다상까지 확장되며, Thomann 외의 전속도 IMEX 기법이 그 대표 사례다.