[논문 리뷰] 이류와 압력을 완전히 갈라놓는다 — Baer–Nunziato flux의 TV 분리

Tokareva와 Toro의 2016년 논문은 재구현이 30분이면 끝날 것처럼 보였다. flux를 둘로 쪼개고, 각각 따로 계산해서 더하라는 것이 전부였으니까. 실제로 구현은 40줄로 끝났다. 강한 충격파 두 개가 부딪히는 Riemann 문제에서 Rusanov flux보다 1.5배 정확했다. 그리고 Noh 문제로 넘어가자 세 번째 스텝에서 압력이 −0.12가 됐다. 오늘은 그 분리가 정확히 무엇을 하는지, 무엇을 공짜로 얻고 어디서 값을 치르는지를 코드로 따라간다.

논문이 쪼개는 방식#

• 저자: S. A. Tokareva, E. F. Toro
• 제목: A flux splitting method for the Baer–Nunziato equations of compressible two-phase flow
• 출처: Journal of Computational Physics 323 (2016) 45–74
• DOI: 10.1016/j.jcp.2016.07.019

한 줄 요약: Toro–Vázquez(TV) flux 분리를 Baer–Nunziato(BN·두 상에 각자의 속도와 압력을 주는 7-방정식 모델) 방정식으로 확장하고, 쪼개진 부분계 중 압력 쪽만 값싸게 풀어 전체 flux를 만든다.

BN 방정식은 보존형으로 쓸 수 없다. 1D x-분할 형태는 비보존항을 따로 안고 있다.

$\mathbf{Q}$는 7개 보존량, $\mathbf{F}$는 보존 flux, $\mathbf{T}\,\partial_x\bar\alpha$는 체적분율 구배에 붙는 비보존항이다. 바(¯)는 고상, 맨 기호는 기상을 뜻한다.

논문의 출발점은 $\mathbf{F}$의 보존 부분을 두 조각으로 자르는 것이다. 기상 3성분만 적으면 이렇다.

$\mathbf{A}$에는 압력이 한 항도 없다. $\mathbf{P}$에는 $\rho u$ 이류가 한 항도 없다. 둘을 더하면 원래 Euler flux가 정확히 복원된다. 여기서 이류계(A-system)와 압력계(P-system)가 나온다. 비보존항은 압력계 쪽에 붙는다.

쪼개고 나면 음속이 한쪽에만 남는다#

쪼갠 보람은 각 부분계의 특성속도를 보면 드러난다. 이류계의 Jacobian은 세 번째 열이 0이라 고윳값이 $\{0, u, u\}$다. 음속이 없다. 압력계는 원시변수 $(\rho, u, p)$에 대해

$h$는 비엔탈피, $e_p = \partial e/\partial p$다. 이상기체를 넣어 보면 $\rho e_p = 1/(\gamma-1)$이므로 $h/(\rho e_p) = \gamma p/\rho = a^2$, 즉 $A = \sqrt{u^2 + 4a^2}$다. 고상의 stiffened EOS도 똑같이 $\bar A = \sqrt{\bar u^2 + 4\bar a^2}$로 정리된다.

$u = 0$을 넣으면 $\lambda_{1,3} = \mp a$다. 음속이 압력계에만 통째로 들어가 있다는 뜻이다. 시간 스텝을 묶는 것은 이류계가 아니다.

아래 시뮬레이션에서 네 파라미터를 직접 움직여보자.

a를 0.1에서 2.0으로 끌어올리면 가운데 패널(이류계)의 부채는 전혀 벌어지지 않고 오른쪽 패널(압력계)만 넓어진다. 각 패널 아래 |λ|max 막대를 같이 보면, 음향 CFL 제약이 어느 쪽에서 오는지가 한눈에 보인다.

BN이 더하는 것은 λ₇ = ū 하나#

기상·고상 각각 3개씩 6개, 거기에 일곱 번째 고윳값 $\lambda_7 = \bar u$가 붙는다. 논문이 고유벡터를 전개해 얻은 결론이 실무적으로 중요하다. 체적분율 $\bar\alpha$는 오직 $\lambda_7$ 장을 지나며 뛴다. 기상 밀도는 $\lambda_{1,3}$을 지나며 변하지 않고, 기상 압력은 $\lambda_2$에서 일정하다.

덕분에 Godunov 상태 표집이 싸진다. 아음속 배치 $S_L < \bar u < S_R$을 가정하면 $\lambda_1 < \lambda_2 = 0 < \lambda_3$이 항상 성립하므로, 계면 상태는 $\lambda_7 = \bar u$의 부호만으로 결정된다. 엔트로피 픽스도 필요 없다. 논문이 CPU 시간을 아꼈다고 말하는 지점이다. 위 시뮬레이션 아래쪽 초록/빨강 배지가 이 아음속 조건을 실시간으로 검사한다. ū를 ±2로 밀어보면 논문이 다루지 않는 영역으로 넘어간다.

직선 두 개로 끝내는 P-system Riemann solver#

압력계의 비선형 파에 걸리는 일반화 Riemann 불변량은 이렇다.

첫 등식이 $\rho = \text{const}$를 준다. 남은 ODE $du/dp = 2/[\rho(u-A)]$를 논문은 특성곡선의 발에서 계수를 얼려 선형화한다. $C_L = \rho_L(u_L - A_L)$로 고정하면 적분이 직선 하나로 끝난다.

$C_R = \rho_R(u_R + A_R)$이다. 두 직선을 연립하면 반복 없이 닫힌 해가 나온다.

$A > |u|$이므로 $C_L < 0 < C_R$이 항상 성립한다. 분모가 0이 되는 경우가 구조적으로 없다. 나눗셈 보호 코드를 넣을 자리가 없다는 뜻이고, 실제로 넣지 않아도 됐다.

다만 이 선형화가 얼마나 거친지는 직접 봐야 한다. 아래에서 좌·우 상태를 움직여보자.

흰 곡선이 불변량을 RK4로 정직하게 적분한 것이고, 점선이 논문이 실제로 쓰는 접선이다. p_L/p_R을 1.1로 낮추면 두 별 상태가 2% 안에서 겹친다. 기본값인 3에서 이미 12%, p_L을 40까지 올리면 67%까지 벌어진다. 접선은 곡선이 간 곳에 도달하지 못한다.

코드 — 40줄로 끝나는 TV flux#

이상기체 1상으로 축약해 구현했다. 압력계 상태에서 $\rho^*$가 필요 없다는 점이 편하다. $\rho e = p/(\gamma-1)$이라 P-flux 세 번째 성분이 $\gamma/(\gamma-1)\,u^* p^*$로 정리되기 때문이다.

import numpy as np
 
def prim(q, g):
    rho = q[0]; u = q[1] / rho
    return rho, u, (g - 1.0) * (q[2] - 0.5 * rho * u * u)
 
def pressure_wave_speed(rho, u, p, g):
    """A = sqrt(u^2 + 4h/(rho e_p)) -> sqrt(u^2 + 4a^2)   [논문 eq. 9]"""
    return np.sqrt(u * u + 4.0 * g * p / rho)
 
def advection_flux(rho, u):
    """A(Q): 압력항이 하나도 없다   [논문 eq. 5]"""
    return np.array([rho * u, rho * u * u, 0.5 * rho * u ** 3])
 
def p_system_star(rhoL, uL, pL, rhoR, uR, pR, g):
    """압력계 선형화 Riemann solver   [논문 eq. 15, 19]"""
    CL = rhoL * (uL - pressure_wave_speed(rhoL, uL, pL, g))   # 항상 < 0
    CR = rhoR * (uR + pressure_wave_speed(rhoR, uR, pR, g))   # 항상 > 0
    us = (2.0 * (pR - pL) + CL * uL - CR * uR) / (CL - CR)
    return us, pL + 0.5 * CL * (us - uL)
 
def tv_face_flux(qL, qR, g):
    rhoL, uL, pL = prim(qL, g)
    rhoR, uR, pR = prim(qR, g)
    us, ps = p_system_star(rhoL, uL, pL, rhoR, uR, pR, g)
    # 이류계 고윳값은 {0, u, u} -> 속도 부호 하나로 업윈드가 끝난다
    fa = advection_flux(rhoL, uL) if us >= 0.0 else advection_flux(rhoR, uR)
    # 압력계: rho e = p/(g-1) 이므로 rho* 를 알 필요가 없다
    fp = np.array([0.0, ps, g / (g - 1.0) * us * ps])
    return fa + fp

Toro Test 4에서 Rusanov와 붙여봤다#

토이 문제는 강한 충격파 두 개가 마주 달리는 Toro Test 4를 골랐다. 좌 $(\rho,u,p) = (5.99924,\ 19.5975,\ 460.894)$, 우 $(5.99242,\ -6.19633,\ 46.0950)$, $\gamma = 1.4$, $t = 0.035$, CFL 0.9. 정지 영역이 없어 업윈드 방향이 흔들리지 않는다. 정확해는 Toro의 반복 해법으로 따로 구현해 $p^* = 1691.647$, $u^* = 8.6898$을 얻었다.

같은 1차 정확도, 같은 비용대에서 TV 분리가 일관되게 1.5배 정확하다. 흥미로운 건 앞에서 본 별 상태 오차다. 이 조건의 압력비는 10이고, 압력계 별 상태는 정확값에서 65% 벗어난다. 그런데도 스킴은 이긴다. 1차 FV에서 지배 오차는 $O(\Delta x)$ 뭉갬이고, P-flux는 정확할 필요 없이 일관되고 산일적이면 되기 때문이다.

Noh 문제에서 갈라진 자리#

Noh 문제로 옮겼다. $\gamma = 5/3$, $\rho = 1$, $u = -1$이 벽으로 밀려들고, 정확해는 $\rho = 4$, $p = 4/3$, $u = 0$, 충격파가 $x = t/3$에 선다. 후방 영역 전체가 $u \approx 0$인 정체 영역이다.

TV 분리는 두 갈래로 죽었다. 이류 flux를 $u^*$의 부호로 업윈드하면 $N = 200$에서 3번째 스텝, 벽에서 두 번째 셀에 $p = -0.12$가 떴다. 대신 $\tfrac12(u_L + u_R)$의 부호로 업윈드하면 살아남기는 하는데 밀도가 격자를 잘게 할수록 더 심하게 톱니를 만든다.

$\mathrm{TV}(\rho)$는 후방 영역 밀도의 전변동이다. 격자를 8배 늘리면 Rusanov는 6배 줄어드는데 TV 분리는 3배 늘어난다. 수렴하지 않는다. 평균값 자체는 $\bar\rho = 4.11$, Rusanov는 3.99로 둘 다 맞지만, 프로파일이 다르다.

원인은 분리 그 자체에 있다. 질량 flux는 $\mathbf{A}$의 첫 성분 $\rho u$ 하나뿐이고, $\mathbf{P}$의 첫 성분은 0이다. 즉 압력계는 밀도를 건드리지 않는다. 정체 영역에서 $u \to 0$이면 질량 flux가 함께 0으로 가고, 밀도장에 남는 수치 산일이 없어진다. 어느 쪽 상태를 업윈드로 고르든 결과가 같아지므로 홀짝 셀이 분리된다. 같은 이류 flux 선택이 Toro Test 4에서는 $L_1$을 0.594에서 0.568로 바꿀 뿐이었다. $u \neq 0$인 곳에서는 거의 무관하고, $u \approx 0$인 곳에서는 결정적이다.

논문이 넘어간 것#

논문의 시험 문제 6개는 전부 Riemann 문제다. 초기 불연속을 가로질러 속도가 0이 아니다. 정체 영역을 만드는 문제, 벽 반사, 정상 유동은 한 개도 없다. 위 결과는 그 사각지대가 실제로 존재한다는 것을 보여준다. 논문이 이류계 flux를 Toro–Vázquez 원논문에 넘긴 것도 아쉽다. "straightforward"라고만 적혀 있는데, 실제로는 이 선택이 정체 영역의 생사를 가른다.

아음속 제약도 남는다. $S_L < \bar u < S_R$은 "많은 저자가 물리적으로 더 타당하다고 본다"는 근거로 도입되지만, 고상이 기상 음속을 넘는 배치는 폭연-폭발 전이(DDT)에서 실제로 나타난다. 고상 접촉면에서는 여전히 얇은층(thin-layer) 비선형 연립을 반복법으로 풀어야 한다는 점도 "반복 없음"이라는 인상과는 다르다. 다만 논문은 수렴 문제를 겪지 않았다고 명시한다.

효율 수치는 설득력이 있다. 선형화 solver가 가장 빠르고 가장 가까운 경쟁자가 8배 느리다. 반면 HLLEM-TV는 125배, 수치 Roe는 67배 비싸다. 고유벡터를 평가해야 하는 순간 TV 분리의 매력은 사라진다는 저자들의 결론이 그대로 숫자로 나온다. OpenFOAM 계열에서 이 구조를 흉내내려면 surfaceScalarField 두 개를 따로 만들어 더하는 형태가 되는데, 압력계 쪽만 별도 Riemann solver로 채우는 것은 가능하다.

재현 가능성 점수#

8/10. 분리 정의(식 3–5), 압력계 고유구조(식 9–10), 선형화 관계식(식 15·19·23·26)은 그대로 코드가 됐고 검산도 맞았다. $u=0$에서 $A/2 = a$가 기계 정밀도로 일치하고, 접선이 정확 곡선의 $O(\Delta p^2)$ 근사인 것도 확인됐다. 감점 2점은 이류계 flux가 이 논문 안에 없다는 것, 그리고 고상 접촉면의 얇은층 연립 전개가 압축돼 있어 완전한 BN 재현까지는 원논문 [1]과 [4]를 함께 봐야 한다는 것이다.

다음에 읽을 것은 Toro & Vázquez(2012), Computers & Fluids 70, 1–12. 이 분리가 처음 나온 곳이고, 오늘 문제가 된 이류계 flux가 거기 있다.