[논문리뷰] 단 한 스테이지로 2차 정확도 — Kapila 다상류를 위한 GRP 솔버

2차 정확도를 얻으려면 Runge-Kutta를 여러 단 돌려야 한다고 외우고 있었다. 이 논문은 단 한 스테이지로 그 차수를 만든다. 계면에서 값과 함께 시간 미분을 돌려주는 GRP(Generalized Riemann Problem) 솔버다.

논문: Tuowei Chen, Zhifang Du, Generalized Riemann Problem Method for the Kapila Model of Compressible Multiphase Flows, arXiv:2310.08241 (2023).

GRP 체크를 끄면 기울기 슬라이더가 잠기고 계면 값이 시간에 얼어붙는다. CFL을 올리면 플럭스에 쓰이는 반스텝 지점만 앞으로 나간다.

기울기를 0으로 두면 시안 점이 멈추고 Godunov와 포개진다. 이류 속도 $a$의 부호를 뒤집으면 값을 가져오는 셀이 바뀐다. 주황 원은 $\Delta t/2$ 시각이다.

한 스테이지로 2차 — GRP가 푸는 문제#

유한체적법의 한 스텝은 계면 플럭스 $\mathbf{F}_{i+1/2}$로 결정된다. Godunov는 계면에서 Riemann 문제를 풀어 값 하나를 얻는다. 그 값은 시간에 대해 상수다. 시간 2차를 원하면 같은 공간 연산을 Runge-Kutta처럼 여러 번 반복하게 된다.

GRP는 플럭스를 중간 시각 $t^{n+1/2}$에서 평가한다.

$\mathbf{U}^{n,*}$는 Riemann 해(계면 값)이고, $(\partial \mathbf{U}/\partial t)^{n,*}$는 그 계면의 순간 시간 미분이다. 두 항을 Taylor 전개로 합치면 계면 값이 $\Delta t/2$만큼 미래로 간다. 평가는 한 번이다.

조각상수에서 조각선형으로 — Riemann과 GRP의 차이#

Riemann 문제는 계면 양쪽이 조각상수다. GRP 문제는 양쪽이 조각선형이다. 초기 불연속에서 뻗는 파동은 더 이상 직선이 아니고, 솔버는 Riemann 해와 시간 미분을 같이 돌려준다.

스칼라 이류 $u_t + a u_x = 0$에서 셀 $i$의 기울기를 $\sigma_i$라 하자. $a>0$이면 계면의 Riemann 해는 왼쪽 셀의 외삽이다.

$u_i$는 셀 평균, $\sigma_i$는 셀 안 기울기, $\Delta x$는 셀 크기다. 앞의 시안 점이 이 식이고, 기울기를 0으로 둔 환원도 이 식에서 나온다.

Lax-Wendroff 절차 — 시간 미분을 공간 변화로#

시간 미분은 지배방정식에서 꺼낸다.

왼쪽의 시간 미분을 오른쪽의 공간 미분으로 바꾼다. $\sigma_i$는 재구성에서 이미 알고 있다. 반스텝 계면 값은

이고, 플럭스는 $F = a\,U^{n+1/2}$다. Euler나 Kapila에서는 $\partial u/\partial t = -A\,\partial u/\partial x$로 일반화되며 $A$는 자코비안이다. 논문은 이 절차로 음향파를 선형화해 GRP를 근사한다.

Python — GRP 스칼라 솔버로 수렴 차수를 잰다#

매끄러운 초기조건 $u_0=\sin(2\pi x)$를 한 주기 이류시켜 Godunov와 GRP의 $L_1$ 오차 차수를 비교한다. 영역은 주기 경계의 선형 이류다.

import numpy as np
 
# 1D 선형 이류  u_t + a u_x = 0,  주기 영역 [0,1]
# 1차 Godunov vs 단일 스테이지 2차 GRP(Lax-Wendroff) 비교
 
def cell_gradient(u, dx):
    # 매끄러운 해의 차수 측정을 위해 중심(무제한) 기울기 사용
    return (np.roll(u, -1) - np.roll(u, 1)) / (2.0 * dx)
 
def godunov_step(u, a, dx, dt):
    # 조각상수: 업윈드 계면 값 (a > 0)
    flux = a * u                       # F_{i+1/2} = a u_i
    return u - dt / dx * (flux - np.roll(flux, 1))
 
def grp_step(u, a, dx, dt):
    g = cell_gradient(u, dx)           # 공간 기울기 sigma_i
    u_star = u + g * dx / 2.0          # 계면 Riemann 해 U*
    u_half = u_star - a * g * dt / 2.0 # Lax-Wendroff 반스텝 U^{n+1/2}
    flux = a * u_half                  # F_{i+1/2}^{n+1/2}
    return u - dt / dx * (flux - np.roll(flux, 1))
 
def l1_error(num, exact, dx):
    return np.sum(np.abs(num - exact)) * dx
 
def convergence_study(a=1.0, cfl=0.4, t_end=1.0):
    print(f"{'N':>5} {'Godunov':>16} {'GRP':>16}")
    prev = {}
    for N in [40, 80, 160, 320]:
        dx = 1.0 / N
        x = (np.arange(N) + 0.5) * dx
        u0 = np.sin(2 * np.pi * x)     # 매끄러운 초기조건
        dt = cfl * dx / abs(a)
        nsteps = int(round(t_end / dt))
        dt = t_end / nsteps            # t_end에 정확히 도달
        ug = u0.copy(); ur = u0.copy()
        for _ in range(nsteps):
            ug = godunov_step(ug, a, dx, dt)
            ur = grp_step(ur, a, dx, dt)
        exact = np.sin(2 * np.pi * (x - a * t_end))
        eg = l1_error(ug, exact, dx)
        er = l1_error(ur, exact, dx)
        og = f"{np.log2(prev['g']/eg):.2f}" if 'g' in prev else "  - "
        orr = f"{np.log2(prev['r']/er):.2f}" if 'r' in prev else "  - "
        print(f"{N:>5} {eg:.2e}({og}) {er:.2e}({orr})")
        prev = {'g': eg, 'r': er}
 
convergence_study()

출력은 다음과 같다.

    N          Godunov              GRP
   40 1.63e-01(  - ) 1.31e-03(  - )
   80 8.76e-02(0.90) 2.69e-04(2.28)
  160 4.54e-02(0.95) 6.32e-05(2.09)
  320 2.31e-02(0.97) 1.55e-05(2.03)

$\lambda\Delta t$ 슬라이더를 1.6 위로 올리면 명시적 Euler가 $[0,1]$ 밴드를 뚫는지가 바뀌고, 초기 부피분율과 평형값을 움직여도 Crank-Nicolson이 밴드 안에 남는지를 같이 본다.

$\lambda\Delta t=3$까지 시안은 평형으로 내려가고 빨강은 진동한다. 방금 표의 차수 2.03은 매끄러운 이류의 숫자이고, 이 슬라이더는 그 반스텝의 시간 미분이 소스 적분에 들어갈 때의 이야기다.

N=40에서 오차 비는 $1.63\times 10^{-1}/1.31\times 10^{-3}$, 약 124배다. 같은 격자에서 100배 작다는 말은 이 격자를 가리킨다. N=320에서는 $2.31\times 10^{-2}/1.55\times 10^{-5}$, 약 1490배다. 격자를 40에서 320으로 8배 잘게 하면 Godunov 오차는 7.1배만 줄고 GRP는 85배 준다. 8의 제곱은 64다. 인쇄된 차수가 2.28에서 2.09, 2.03으로 내려앉는 동안 64를 조금 넘긴 것이다. Godunov 쪽 0.90, 0.95, 0.97도 같은 방향으로 1에 붙는다.

코드의 CFL은 0.4다. 반스텝 보정의 크기는 $|a\,\sigma_i\,\Delta t/2|$이고, 계면 외삽의 크기는 $|\sigma_i\,\Delta x/2|$다. 둘의 비는 $|a|\Delta t/\Delta x$라서 CFL과 같다. 0.4는 기울기가 계면을 밀어 둔 거리의 40%를, 플럭스를 읽는 시각이 다시 당긴다는 뜻이다. CFL을 0.1로 내리면 그 보정은 외삽의 10%가 된다. 차수 증명 자체는 남는다. 한 스텝에서 시간 미분이 플럭스에 넣는 양만 작아진다.

N=80의 비는 $8.76\times 10^{-2}/2.69\times 10^{-4}$, 약 326배다. N=160은 $4.54\times 10^{-2}/6.32\times 10^{-5}$, 약 718배다. 격자를 두 배씩 줄일 때마다 오차 비는 124, 326, 718, 1490으로 늘어난다. 한 단계에 두 배를 조금 넘는다.

Kapila 모델의 뻣뻣한 부피분율 소스#

Kapila 5방정식 모델(BN 7방정식의 축소형)의 부피분율 방정식에는 뻣뻣한 소스가 붙는다.

$\alpha_1$은 1상의 부피분율, $K$는 Wood 음속에서 나오는 압축·팽창 계수다. 이 $K$가 부피분율에 대해 어떻게 꺾이는지는 Wood 음속이 혼합물에서 내려앉는 계산과 같은 자리다. 소스는 $\alpha_1$을 $[0,1]$ 밖으로 밀 수 있다. 명시적 적분은 강성이 커지면 경계를 넘어 발산한다.

GRP의 시간 미분은 계면 값을 새 시간 $t^{n+1}$에서 알게 한다. Gauss-Green으로 $\nabla\cdot\mathbf{u}$의 셀 평균을 그 시각에 추정할 수 있다. 소스를 Crank-Nicolson(반음함수)으로 적분하고, 부피분율 경계를 지키는 시간 간격 제한까지 세운다. 앞의 $\lambda\Delta t$ 곡선이 그 제한의 모형이다.

재현하며 의심스러웠던 것#

스칼라 GRP는 2차가 나왔지만 Kapila 시스템은 그보다 거칠다. 무제한 중심 기울기는 충격파에서 진동한다. 차수가 나온 것은 테스트가 매끄러워서다. 불연속에서는 minmod 같은 제한자가 필요하고, 그때 차수는 국소적으로 1차로 떨어진다.

Wood 음속은 부피분율에 대해 비단조다. 계면의 수치 확산대에서 마하수가 흔들리는 원인이다. GRP의 시간 미분은 압축과 팽창을 잘 따라가지만, 그 비단조성 자체를 지우지는 못한다. 논문도 그 부분은 정면으로 풀기보다 강건성으로 비킨다.

실무 코드는 대개 Strang 분할에 Runge-Kutta를 붙인 구조라, 단일 스테이지 GRP를 그대로 넣기 어렵다. 이점을 살리려면 플럭스 함수를 Lax-Wendroff 형태로 다시 써야 한다. 이식 비용이 남는다. 선형화 음향 근사의 출전은 Ben-Artzi & Falcovitz의 가스동역학 GRP다.

소스를 플럭스 밖에 둔 채 차수를 기대한 실수#

표의 2.03은 기울기가 매끄러울 때의 숫자다. 충격파에서 제한자를 켜면 그 자리는 국소적으로 1차로 돌아간다. 차수를 보고 스테이지를 지우기 전에, 플럭스 함수를 다시 쓸 예산이 있는지가 먼저다.

소스를 플럭스에 녹이지 않고 밖으로 빼면, 같은 강성을 다른 시간 적분이 떠안는다. 압력을 두 개로 풀고 하나로 합치는 이완은 그 분리 쪽에서 부피분율 경계를 지키는가, 아니면 $\Delta t$를 다시 잘게 쪼개는가. 이 출력의 2.03은 그 질문에 답이 아니다.