음파와 물질을 따로 흘려보낸다 — 음향·대류 분리(Lagrange–Projection)#

음속은 초당 340미터다. 시속 5킬로로 걷는 사람 옆에서도 소리는 그 속도다. 압축성 솔버는 둘을 한 시간刻에 넣어야 한다. 비율이 문제다. 아주 느린 흐름이면 음파가 물질보다 수백 배 빠르다. 직접 Godunov 계열은 그 빠른 파에 시간刻이 묶인다. 느린 대류는 수치확산에 뭉개진다.

Mach 슬라이더를 내리면 압력 펄스가 갈라지는 속도와, 엔트로피 덩어리가 멈추는 자리가 같이 바뀐다.

0.02까지 내리면 amber 펄스가 좌우로 갈라져 달리고, teal 덩어리는 거의 제자리다. 삼각형 세 개의 속도가 $u\!-\!c$, $u$, $u\!+\!c$다. 느릴수록 가운데만 멈추고 양쪽 간격은 $2c$로 남는다.

ten Eikelder 등(2017)은 이 둘을 아예 떼어 놓는다. 지배방정식을 음향 부분과 대류 부분으로 쪼개, 각각 다른 솔버로 번갈아 푼다. 단상(single-phase) Euler에 그 Lagrange–Projection을 올려 "느린 흐름 속 음향 펄스"로 숫자를 맞춘다.

“

논문: M.F.P. ten Eikelder, F. Daude, B. Koren, A.S. Tijsseling, "An acoustic-convective splitting-based approach for the Kapila two-phase flow model", Journal of Computational Physics 331 (2017) 188–208. DOI: 10.1016/j.jcp.2016.11.031

야코비안이 둘로 쪼개진다#

1차원 Euler 방정식을 원시변수 $\mathbf{W}=(\rho,u,p)^T$로 쓰면 준선형 형태가 된다.

$\rho$는 밀도, $u$는 속도, $p$는 압력이다. 계수행렬 $\mathbf{B}$는 두 조각으로 정확히 갈라진다.

$\mathbf{A}$는 압력 항을 모두 담는 음향(acoustic) 부분이다. $\mathbf{C}=u\mathbf{I}$는 물질을 그대로 실어 나르는 대류(convective) 부분이다. 분리는 사실상 라그랑주 미분 $\mathrm{D}/\mathrm{D}t = \partial_t + u\partial_x$에서 $u\partial_x$ 항을 떼어내는 일이다.

고유값의 덧셈 분해#

이 분리가 왜 강력한지는 고유값을 보면 드러난다. 전체 계의 파속은 $u-c,\ u,\ u+c$다. 이것이 음향과 대류로 정확히 더해진다.

$c=\sqrt{\gamma p/\rho}$는 음속이다. 음향 파속은 $\pm c$로 흐름 속도와 무관하다. 대류 파속은 통째로 $u$다. 저마하 극한 $M=u/c\to 0$에서 음향 밴드는 폭 $2c$로 그대로인데 대류 속도만 0으로 줄어든다. 이 격차가 강성(stiffness)의 정체다. 앞 화면의 삼각형 간격이 이 덧셈과 같다. $M\to 0$이면 $\lambda^c$만 사라지고 $\lambda^a=\pm c$는 남는다.

음향 단계 — 라그랑주 좌표의 HLLC#

음향 부분은 압력이 만드는 팽창·압축이다. 논문은 이를 질량 좌표(라그랑주)로 옮겨 HLLC형 리만 솔버로 푼다. 면 $j+1/2$의 별표 상태는 두 식으로 끝난다.

$a_{j+1/2}=\max(\rho_j c_j,\ \rho_{j+1}c_{j+1})$는 음향 임피던스(impedance, 밀도×음속)다. 이 별표 속도·압력만으로 오일러 변수 갱신이 이뤄진다. 각 셀의 팽창률은 하나의 계수로 압축된다.

$R_j$는 음향 단계가 셀 부피를 얼마나 늘리거나 줄였는지를 나타낸다. 밀도는 $\rho^{n+1-}_j = \rho^n_j / R_j$로 바로 나온다.

배경은 $c_0=1$, $\rho_0=1$, $\gamma=1.4$라 $p_0=\rho_0 c_0^2/\gamma=1/1.4$다. 압력 펄스 진폭은 $10^{-3}$, 밀도 이상은 $5\times 10^{-2}$다. 중심이 $x=0.5$와 $x=0.25$로 떨어져 초기에 겹치지 않는다. 밀도 최댓값은 엔트로피 덩어리 쪽, 대략 $1.05$에서 시작한다. $R_j$가 1보다 크면 그 셀은 팽창하고 밀도는 줄어든다. 1보다 작으면 압축이다. $\Delta t/\Delta x$가 면의 $u^*$ 차이에 곱해지므로, CFL 0.8에서도 $R_j$는 1 근처를 크게 벗어나지 않는다. 격자 400, 영역 길이 1이면 $\Delta x=1/400$이다.

대류 단계 — 업윈드 투영#

음향 단계가 만든 별표 속도 $u^*$로 이제 물질을 실어 나른다. 보존량 $\varphi\in\{\rho,\rho u,\rho E\}$를 단순 업윈드로 갱신한다.

$\varphi^{\text{up}}_{j+1/2}$는 $u^*\!\ge\!0$이면 $\varphi_j$, 아니면 $\varphi_{j+1}$이다. 앞 단계의 $R_j$가 여기서 되돌아오며 질량·운동량·에너지가 정확히 보존된다. 두 단계를 순서대로 밟으면 한 시간전진이 끝난다.

$u^*$가 양이면 면 값은 왼쪽 셀이다. 음수면 오른쪽 셀이다. $R_j\,\varphi^{n+1-}$가 먼저 부피 변화를 되돌리고, 그다음 면 플럭스가 물질을 옮긴다. 순서가 바뀌면 질량과 $\rho u$, $\rho E$의 합이 어긋난다. 주기 경계라 면 플럭스는 도메인 밖으로 새지 않는다. $R_j$ 되먹임이 맞으면 운동량 합에는 반올림만 남는다. 저마하에서 직접 Roe가 퍼지는 폭은 Hope–Collins의 인공확산 세 스케일링이 따로 잰다.

Python — 느린 흐름 속 음향 펄스#

전체 파이프라인을 numpy로 짠다. 주기 경계, 이상기체, 셀당 세 보존량이다. 초기조건은 배경 흐름 $u_0=Mc_0$ 위에 작은 압력 펄스와 밀도 이상을 얹는다.

import numpy as np
 
GAMMA = 1.4  # 이상기체 비열비
 
def primitives(rho, mom, Ene):
    """보존량 -> 원시변수 (속도, 압력, 음속)."""
    u = mom / rho
    e = Ene / rho - 0.5 * u * u            # 내부 비에너지
    p = (GAMMA - 1.0) * rho * e
    c = np.sqrt(GAMMA * p / rho)
    return u, p, c
 
def acoustic_faces(rho, u, p, c):
    """면 j+1/2의 HLLC 음향 상태 u*, p* (논문 식 34)."""
    rp, up, pp, cp = (np.roll(a, -1) for a in (rho, u, p, c))
    a = np.maximum(rho * c, rp * cp)       # 음향 임피던스 a = max(rho*c) (식 32)
    ustar = 0.5 * (u + up) + (p - pp) / (2 * a)
    pstar = 0.5 * (p + pp) + 0.5 * a * (u - up)
    return ustar, pstar
 
def lagrange_projection_step(rho, mom, Ene, dx, dt):
    """음향 단계 -> 대류 단계 (논문 식 38, 40)."""
    u, p, c = primitives(rho, mom, Ene)
    uf, pf = acoustic_faces(rho, u, p, c)          # j+1/2
    uf_m, pf_m = np.roll(uf, 1), np.roll(pf, 1)    # j-1/2
    lam = dt / dx
    # 1) 음향 단계: 압력이 만드는 팽창/압축
    R = 1.0 + lam * (uf - uf_m)                     # 식 (39)
    rho1 = rho / R
    mom1 = (mom - lam * (pf - pf_m)) / R
    Ene1 = (Ene - lam * (pf * uf - pf_m * uf_m)) / R
    # 2) 대류 단계: 물질을 u*로 실어 나른다 (업윈드 투영)
    def project(phi1):
        phi_f = np.where(uf >= 0, phi1, np.roll(phi1, -1))
        phi_f_m = np.roll(phi_f, 1)
        return R * phi1 - lam * (uf * phi_f - uf_m * phi_f_m)
    return project(rho1), project(mom1), project(Ene1)
 
def run_acoustic_pulse(mach, n=400, cfl=0.8, tmax=0.25):
    x = (np.arange(n) + 0.5) / n
    dx = 1.0 / n
    c0, rho0 = 1.0, 1.0
    p0 = rho0 * c0**2 / GAMMA
    u0 = mach * c0
    dp = 1e-3 * np.exp(-((x - 0.5) / 0.03)**2)     # 음향 압력 펄스
    ds = 5e-2 * np.exp(-((x - 0.25) / 0.03)**2)    # 엔트로피(밀도) 이상 -> 접촉파
    rho = rho0 + dp / c0**2 + ds
    u = np.full(n, u0)
    p = p0 + dp
    mom = rho * u
    Ene = p / (GAMMA - 1) + 0.5 * rho * u * u
    t, m0 = 0.0, mom.sum()
    while t < tmax:
        _, _, c = primitives(rho, mom, Ene)
        dt = min(cfl * dx / np.max(np.abs(mom / rho) + c), tmax - t)
        rho, mom, Ene = lagrange_projection_step(rho, mom, Ene, dx, dt)
        t += dt
    return rho, mom, m0, mom.sum()
 
for M in (0.02, 0.2):
    rho, mom, m0, m1 = run_acoustic_pulse(M)
    print(f"M={M}: 운동량 보존오차={abs(m1-m0)/abs(m0):.1e}, rho_max={rho.max():.4f}")
# M=0.02: 운동량 보존오차=2.2e-16, rho_max=1.0497
# M=0.2:  운동량 보존오차=0.0e+00, rho_max=1.0452

운동량이 기계정밀도로 보존된다. 압력은 유계, 밀도는 양수로 남는다. 분리했는데도 보존이 깨지지 않는 이유는 대류 단계에 $R_j$를 되돌려 넣었기 때문이다.

출력의 $\rho_{\max}$는 $M=0.02$에서 1.0497, $M=0.2$에서 1.0452다. 초기 봉우리 1.05에서 거의 안 깎였다. 압력 펄스가 밀도에 더하는 양은 $dp/c_0^2=10^{-3}$뿐이라, 최댓값을 만드는 쪽은 접촉파다. 운동량 오차 $2.2\times 10^{-16}$과 0은 분리가 보존을 깨지 않는다는 숫자다. 시간刻을 $|u|+c$로 잡았기 때문에 $M$을 0.02에서 0.2로 올려도 스텝 수는 거의 같다. $c_0=1$이면 $|u|+c$가 1.02와 1.2다. 대류가 열 배 빨라져도 음향이 시간 간격의 대부분을 가져간다.

저마하에서 갈리는 것#

분리가 사는 곳은 저마하다. 직접법의 안정 시간刻은 항상 $|u|+c$에 묶인다. 관심 물리가 $u$ 위에 있어도 그렇다.

Mach 슬라이더를 옮기면 amber 밴드 폭은 $2c$에 고정되고, teal인 $u$만 줄어든다.

음향 밴드의 폭은 $M$과 무관하게 항상 $2c$다. 대류 속도 $u$만 0으로 수렴한다. $c=1$, $M=0.05$이면 $u=0.05$이고 $u+c=1.05$라, 비 $(u+c)/u$는 21이다. 분리 기법은 대류 단계를 음향 단계와 다른 시간刻으로 밟을 수 있다. 느린 물리가 빠른 음파에 끌려가지 않는다. 저마하에서 직접 Godunov가 겪는 과도한 수치확산도 이 분리로 피할 수 있다.

비판적으로 보면#

걸린 게 세 개다. 시간 정확도는 1차다. 2차를 얻으려면 Strang 분리로 음향–대류–음향 순서를 감싸야 하고, 그만큼 비용이 는다. 원 논문의 무대는 Kapila 5방정식 이상류다. 부피분율 방정식의 비보존 항 $K\partial_x u$와 양수성 보장이 실제 난관인데, 단상 축약은 그 어려움을 건너뛴다. 임피던스 $a=\max(\rho c)$는 견고하지만 소산적이다. 강한 충격에서는 이 소산이 접촉면을 뭉갤 수 있다. 논문도 직접법 대비 정확도·효율을 별도로 견준다.

실무에서 이 갈림은 낯설지 않다. OpenFOAM의 압력기반 rhoPimpleFoam이나 all-Mach 계열 솔버가 음향과 대류를 암시적·명시적으로 나눠 다루는 것과 같은 뿌리다. Lagrange–Projection은 그 분리를 리만 솔버 언어로 적은 판본이다.

R_j를 빼 먹으면 보존이 먼저 죽는다#

음향 단계의 $R_j$를 대류 투영에 되돌리지 않으면 질량·운동량·에너지 식이 서로 어긋난다. 위 실행에서 오차가 $2.2\times 10^{-16}$까지 내려간 자리는 그 되먹임이 있을 때뿐이다. Strang으로 감싸 2차를 만들 수는 있다. 순서가 음향–대류–음향으로 늘고 패스도 하나 더 돈다.

Kapila 5방정식으로 가면 부피분율의 $K\partial_x u$가 남는다. 단상으로 줄이면 그 비보존 항을 안 만난다. 임피던스를 $a=\max(\rho c)$로 두면 충격은 견디고 접촉면은 두꺼워진다. 음향을 시간 적분에서 빼는 다른 절단은 Peluchon의 음향-수송 IMEX다. $M=0.05$에서 $(u+c)/u=21$이면, 그 21배를 음향 솔버에 맡길지 대류 CFL에 남길지가 구현의 갈림이다.