[논문 리뷰] 음향을 음해로 풀었는데 Δt는 100배 줄었다 — 상대속도에 숨은 두 번째 1/M²

음향을 음해로 풀었는데 Δt는 줄어들기만 했다#

저마하 2상 코드를 짜다 겪는 전형적인 실패가 있다. 압력 항이 stiff하다는 것을 알고 음해(implicit)로 옮긴다. 음속에서 풀려났다고 생각하고 $\Delta t$를 물질속도 기준으로 키운다. 그런데 터진다. $M=0.01$에서는 $M=1$보다 100배 작은 $\Delta t$를 써야 겨우 돈다.

Lukáčová-Medviďová 등의 2022년 논문은 이 실패의 정확한 위치를 알려준다. 등엔트로피 2상 모델에는 $1/M^2$이 두 군데 있다. 하나는 익숙한 압력이고, 다른 하나는 상대속도 방정식에 걸린 엔탈피 차다. 이 글은 두 번째 것을 찾아내고, 그것을 음해로 옮기면 $\Delta t$ 상한이 마하수에서 실제로 풀리는지를 64셀 격자에서 직접 재본다.

log10 M을 왼쪽으로 밀면 빠른 부채가 같이 벌어진다. f^c implicit을 켜야 허용 $\Delta t$ 막대가 마하수에서 풀린다.

파란 부채는 압력, 분홍 부채는 엔탈피 차다. f^a implicit만 켠 채로는 아래 막대가 여전히 $M$을 따라 줄어든다.

등엔트로피 2상 모델이 들고 있는 두 개의 1/M²#

논문이 다루는 상태벡터는 $W = (\rho,\ \alpha\rho,\ \alpha_1\rho_1,\ \rho u,\ u_1-u_2)^T$다. 혼합 밀도, 부피분율(volume fraction, 셀 안에서 상 1이 차지하는 부피 비율), 부분밀도, 혼합 운동량, 그리고 두 상의 속도차다. 마지막 성분이 이 모델의 특징이다. 상마다 속도를 따로 들고 가지 않고, 혼합 속도 하나와 속도차 하나로 쓴다.

무차원화하면 플럭스에 $M_k$가 붙는다. 운동량 방정식에는 압력이

로 들어온다. 여기서 $\hat{p}_k$는 기준 상태 주변에서 잘라낸(truncated) 압력, $M_k$는 상 $k$의 마하수, $\varrho_k$는 기준 밀도비다. $M\to 0$에서 발산하는 이 항이 첫 번째 stiff 항이고, 단상 IMEX와 똑같이 음해로 처리한다.

문제는 다섯 번째 방정식이다. 속도차는 엔탈피 차로 구동된다.

$\hat{h}_k$는 상 $k$의 잘라낸 엔탈피다. 분모에 같은 $M_k^2$가 있다. 압력을 음해로 옮겨도 이 항은 그대로 남는다.

세 부분계의 특성곡선을 한 장에 겹쳐본다#

논문의 분할은 플럭스를 $f^a$(빠른 음향), $f^b$(물질 이동), $f^c$(stiff 혼합 항) 셋으로 쪼갠 뒤 완화 소스항 $r(W)$를 마지막에 붙인다. 각 조각이 자기 고유값을 갖는다. 앞 다이얼에서 파란 부채와 분홍 부채가 같이 열린 이유가 이 분할이다. f^c까지 음해로 묶는 일이 논문 3.2.3절의 전부다.

상대속도 부분계의 고유값을 직접 적어보기#

$f^c$에 남은 방정식은 부분밀도, 운동량, 속도차 셋이다. 이 $3\times 3$ 계의 고유값은

이고 판별식은

다. $\chi$는 질량분율, 상첨자 $RS$는 기준해(reference solution) 값이다. $\mathcal{A}$의 둘째 항에 $1/M_k^2$이 그대로 살아 있다. 즉 $\lambda^{\pm} \sim \mathcal{O}(1/M)$이다.

$\mathcal{A} > 0$이므로 고유값은 실수다. 부분계는 쌍곡형이고, 그래서 따로 풀 수 있다. 그런데 특성속도가 $1/M$이다. 양해로 풀면 CFL이 $\Delta t \le \Delta x \cdot M / \text{const}$가 된다. 압력을 음해로 옮긴 노력이 여기서 모두 반납된다.

운동량 방정식은 $\lambda^0=0$에 대응하므로 시간스텝을 제약하지 않는다. 그래서 논문은 부분밀도와 속도차, 두 방정식만 음해로 묶는다. 음향과 수송을 나눈 IMEX에서는 빠른 음향 하나만 음해로 넘기면 됐다. 상이 둘이면 stiff 항도 둘이다. 엔탈피 차의 $1/M^2$가 그 두 번째다.

분할을 고를 때의 진짜 제약#

논문이 명시하는 분할 조건은 두 가지다. 각 부분계가 그 자체로 쌍곡형이어야 하고, 쉽게 풀려야 한다. 앞 조건이 까다롭다. 항을 아무렇게나 모으면 고유값이 복소수가 되고, 부분계 단계에서 불안정해진다.

음해 단계를 어떻게 푸는지도 조건에 걸린다. 부분밀도 업데이트에 속도차 업데이트를 대입하면 $(\alpha\rho_1)^{\star\star\star}$에 대한 비선형 타원형 문제가 나온다. 논문은 그것을 선형화 반복으로 푼다. 계수행렬이 대칭 양정치가 되는 근거는 두 양

이 둘 다 양수라는 것이다. $s$는 $\varrho_1(\alpha\rho_1)/\rho = \chi \in (0,1)$이므로 양수다. 이 부호 하나가 "표준 선형 솔버로 효율적으로 풀린다"는 문장을 지탱한다. 연속에서 음속이 빠져도 시간 간격이 사는 성질은 이산이 자동으로 물려받지 않는다. 비이상기체 AP-IMEX가 그 간극을 음속 쪽에서 적었다.

Python으로 같은 격자에 양해·음해를 나란히 걸었다#

$u_1-u_2$를 0 근처에서 선형화하면 $f^c$에 남는 것은 2×2 파동계다.

$m = \alpha\rho_1$, $w = u_1-u_2$다. 파속은 $\sqrt{s\tilde{s}}$이고 $\tilde{s} \propto 1/M^2$이므로 $1/M$에 비례한다. 같은 주기 격자(64셀)에 두 방식을 걸었다. 양해는 정확한 $\lambda_c$로 수치점성을 준 Rusanov, 음해는 $w$를 소거해 얻은 타원형 문제를 주기 삼대각으로 푼다. $\Delta t$는 물질속도 CFL 0.5로 고정한다.

import math
 
N, L = 64, 1.0
dx = L / N
RHO, ALPHA, CHI, U_MAT = 1.0, 0.5, 0.5, 0.2
M0 = CHI * RHO                 # alpha*rho_1, 배경 상태
S = M0 * (1.0 - M0 / RHO)      # 식 (3.25)의 s(W)
 
def stiff_coeff(mach):
    """s-tilde: 엔탈피 계수, O(1/M^2)"""
    return 1.0 / (mach**2 * ALPHA) + 1.0 / (mach**2 * (1.0 - ALPHA))
 
def relative_wave_speed(mach):
    return math.sqrt(S * stiff_coeff(mach))
 
def initial_state():
    m = [M0 + 1e-3 * math.sin(2*math.pi*(i+0.5)*dx) for i in range(N)]
    return m, [0.0] * N
 
def explicit_enthalpy(mach, dt, nstep):
    """(3.19a)+(3.19c)에 Rusanov, 엔탈피 항 양해"""
    k, lam = stiff_coeff(mach), relative_wave_speed(mach)
    m, w = initial_state()
    for _ in range(nstep):
        fm = [0.0]*N; fw = [0.0]*N
        for i in range(N):
            j = (i+1) % N
            fm[i] = 0.5*S*(w[i]+w[j]) - 0.5*lam*(m[j]-m[i])
            fw[i] = 0.5*k*(m[i]+m[j]) - 0.5*lam*(w[j]-w[i])
        m = [m[i] - dt/dx*(fm[i]-fm[i-1]) for i in range(N)]
        w = [w[i] - dt/dx*(fw[i]-fw[i-1]) for i in range(N)]
        if not all(abs(v) < 1e12 for v in m):
            return float('inf')
    return max(abs(v-M0) for v in m)
 
def periodic_tridiag(a, b, c, d):
    """주기 삼대각계를 Sherman-Morrison으로 푼다"""
    n = len(d)
    gam = -b[0]
    bb = b[:]; bb[0] -= gam; bb[n-1] -= c[n-1]*a[0]/gam
    def thomas(rhs):
        cp = [0.0]*n; dp = [0.0]*n
        cp[0] = c[0]/bb[0]; dp[0] = rhs[0]/bb[0]
        for i in range(1, n):
            den = bb[i] - a[i]*cp[i-1]
            cp[i] = c[i]/den
            dp[i] = (rhs[i] - a[i]*dp[i-1])/den
        x = [0.0]*n; x[n-1] = dp[n-1]
        for i in range(n-2, -1, -1):
            x[i] = dp[i] - cp[i]*x[i+1]
        return x
    u = [0.0]*n; u[0] = gam; u[n-1] = c[n-1]
    y = thomas(d); z = thomas(u)
    fac = (y[0] + a[0]*y[n-1]/gam) / (1.0 + z[0] + a[0]*z[n-1]/gam)
    return [y[i] - fac*z[i] for i in range(n)]
 
def implicit_enthalpy(mach, dt, nstep):
    """식 (3.23): w를 소거하고 m에 대한 타원형 문제를 푼다"""
    k = stiff_coeff(mach)
    beta = dt*dt*S*k/(dx*dx)
    a = [-beta]*N; b = [1.0 + 2*beta]*N; c = [-beta]*N
    m, w = initial_state()
    for _ in range(nstep):
        rhs = [m[i] - dt/(2*dx)*S*(w[(i+1) % N]-w[i-1]) for i in range(N)]
        m = periodic_tridiag(a, b, c, rhs)
        w = [w[i] - dt/(2*dx)*k*(m[(i+1) % N]-m[i-1]) for i in range(N)]
    return max(abs(v-M0) for v in m)
 
def max_stable_dt(mach, nstep=200):
    lo, hi = 1e-9, 1.0
    for _ in range(40):
        mid = math.sqrt(lo*hi)
        if explicit_enthalpy(mach, mid, nstep) == float('inf'):
            hi = mid
        else:
            lo = mid
    return lo
 
dt_mat = 0.5*dx/U_MAT
print(f"dx={dx:.5f}  material CFL dt={dt_mat:.5f}  (u={U_MAT})")
print(f"{'M':>8} {'lambda_c':>10} {'dt_exp_max':>11} {'ratio':>7} {'|dm|_exp':>10} {'|dm|_imp':>10}")
for mach in (1.0, 0.3, 0.1, 0.03, 0.01):
    lam = relative_wave_speed(mach)
    dte = max_stable_dt(mach)
    e = explicit_enthalpy(mach, dt_mat, 200)
    i = implicit_enthalpy(mach, dt_mat, 200)
    es = "blow-up" if e == float('inf') else f"{e:.3e}"
    print(f"{mach:8.3g} {lam:10.2f} {dte:11.2e} {dte/dt_mat:7.3f} {es:>10} {i:10.3e}")

dx=0.01562  material CFL dt=0.03906  (u=0.2)
       M   lambda_c  dt_exp_max   ratio   |dm|_exp   |dm|_imp
       1       1.00    1.87e-02   0.479    blow-up  1.718e-06
     0.3       3.33    5.61e-03   0.144    blow-up  1.771e-14
     0.1      10.00    1.87e-03   0.048    blow-up  4.731e-13
    0.03      33.33    5.60e-04   0.014    blow-up  6.318e-12
    0.01     100.00    1.87e-04   0.005    blow-up  5.195e-11

enthalpy implicit을 끄고 log10 M을 $-1$ 아래로 내리면, 고정된 물질 $\Delta t$에서 로그 진폭이 천장에 닿는다.

다시 켜면 같은 $\Delta t$에서 진폭이 가라앉는다. material CFL을 줄여 양해가 살아나는 지점을 찾으면, 그 값이 위 표의 ratio다.

Δt 상한이 M에 정확히 비례한 이유#

dt_exp_max 열을 보면 $M$이 1에서 0.01로 100배 줄 때 상한도 $1.87\times 10^{-2}$에서 $1.87\times 10^{-4}$로 정확히 100배 줄었다. 우연이 아니다. $\lambda_c = \sqrt{s\tilde{s}}$에서 $s$는 $M$과 무관하고 $\tilde{s} \propto M^{-2}$이므로 $\lambda_c \propto M^{-1}$이다. Rusanov의 안정 조건 $\lambda_c \Delta t/\Delta x \le 1$을 그대로 대입하면 $\Delta t_{\max} \propto M$이다.

ratio 열이 실무에서 아픈 숫자다. $M=0.01$에서 물질 CFL로 쓰고 싶은 $\Delta t$의 0.5%만 쓸 수 있다. 같은 물리 시간을 보려면 스텝을 200배 더 밟는다. 반면 음해 열은 고정 $\Delta t$에서 200스텝을 다 돌았고 발산하지 않았다. $M=0.01$에서 음향 CFL은 256이었다.

표의 $\lambda_c$를 거꾸로 읽어 보자. $M=1$에서 $1.00$, $M=0.1$에서 $10.00$, $M=0.01$에서 $100.00$이다. 코드의 $s=0.25$와 $\tilde{s}=4/M^2$를 넣으면 $\sqrt{s\tilde{s}}=1/M$이 되어 인쇄된 파속과 칸이 같다. 물질 간격 $\Delta t=0.03906$을 $M=1$에 대입하면 $\lambda_c\Delta t/\Delta x=0.03906/0.01562=2.50$이다. Rusanov 한계 1을 넘으므로 ratio $0.479$가 나온다. $1.87\times10^{-2}/0.03906=0.479$. 엔탈피만 양해로 두면 $M=1$에서도 물질 CFL 0.5를 못 버틴다. $M=0.01$의 비는 $0.005$다. $0.479/0.005\approx 96$이라, 서두의 100배에 붙는다.

|dm|_imp는 다른 칸이다. 초기 진폭 $10^{-3}$ 대비 $M=1$의 잔차는 $1.718\times10^{-6}$, $M=0.01$은 $5.195\times10^{-11}$이다. 양해 열은 다섯 줄 모두 blow-up이다. 고정 $\Delta t$에서 200스텝을 마친 것은 음해뿐이다. 표가 적은 음향 CFL 256은 그 양해 제한이 어디로 가는지 보여 준다. $\tilde{s}$의 $1/M^2$를 음해로 넘기지 않으면 이 수가 양해 CFL로 되돌아온다.

무조건 안정은 무조건 정확이 아니다#

|dm|_imp 열을 정직하게 읽자. 초기 진폭은 $10^{-3}$이었는데 음해는 $10^{-11}$ 이하로 남겼다. 후진 Euler가 빠른 모드를 사실상 소멸시킨 것이다. 이것은 버그가 아니라 AP(asymptotic preserving, 점근 보존) 스킴이 의도한 동작이다. $\Delta t$가 음향 시간보다 훨씬 크면 그 스케일의 파동은 해상되지 않고 걸러진다. 남기고 싶은 것은 물질속도 스케일의 동역학이다.

다만 이 감쇠를 "정확도"로 착각하면 안 된다. 압축성 영역($M \sim 1$)에서도 같은 스킴을 쓰려면 잘라낸 항 $\bar{p}_k,\ \bar{h}_k$를 되살려 양해로 더해야 한다. 논문 3.2.4절이 그 일을 한다. 잘린 고차 항은 well-prepared 초기자료의 저마하 극한에서 사라지므로, 양해로 넣어도 CFL에 $M$ 의존성을 다시 들여오지 않는다. 그 발산과 감쇠는 출력 직후의 로그 그래프에서 이미 갈라져 있다.

완화항 두 개가 서로 다른 대접을 받는 까닭#

균질 계를 끝내고 남는 것은 완화 소스항이다. 부피분율에는 압력 완화가, 속도차에는 마찰이 걸린다.

$\tau$는 압력 완화 시간, $\zeta$는 마찰 계수다. 둘 다 stiff할 수 있어 음해로 적분하는데, 비용이 다르다.

마찰항은 $u_1-u_2$에 선형이다. 그래서 후진 Euler가 나눗셈 한 번으로 닫힌다.

압력 완화는 다르다. 상 압력이 $\alpha$에 비선형이다. 논문은 $g(\alpha)=0$에 Newton법을 쓰고, 야코비안에 음속이 들어온다.

$\alpha \to 0$ 또는 $1$에서 이 미분이 발산한다. 순수 상 근처에서 Newton이 힘들어지는 지점이 여기다. 코드에서 $\alpha$에 하한·상한을 두는 관행이 이 식 하나로 설명된다. $\tau \to 0$은 순간 완화, $\tau \to \infty$는 균질 모델이다.

엔탈피 차를 음해에서 빼면 Δt가 다시 M에 묶인다#

$1/M^2$이 압력에만 있는 줄 알고 음해로 옮기면, 상대속도 방정식의 $\hat{h}_1/M_1^2-\hat{h}_2/M_2^2$가 양해로 남는다. 그 항의 특성속도가 $1/M$이라, $M=0.01$에서 물질 CFL의 $0.5\%$만 쓸 수 있다. 표의 ratio $0.005$와 음향 CFL 256이 그 자리다.

압력 완화는 또 다른 음해다. $\alpha\to 0$ 또는 $1$에서 $\partial g/\partial\alpha$가 발산하고, Newton이 순수 상 근처에서 버거워진다. 두 압력을 풀었다가 한 압력으로 합치는 이완은 Saurel의 압력 이완이 이미 따로 장부를 갖고 있다. 엔탈피 차를 그 목록에서 빼면, 음향을 음해로 둔 보람이 $\Delta t$에 돌아오지 않는다.