让声音与物质分开流动 — 声学-对流分裂（Lagrange–Projection）#

声速是每秒340米。一个人以每小时5公里散步，他身旁的空气仍以同样的340米传递声音。可压缩求解器必须在一个时间步里同时处理这两种速度。麻烦在于两者之比。在极慢的流动中，声波比物质本身快上数百倍。直接的Godunov类格式，其时间步被这些快速声波拴住。而真正关心的慢速对流，却被数值扩散抹平。

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

雅可比矩阵一分为二

把一维Euler方程写成原始变量 $\mathbf{W}=(\rho,u,p)^T$ ，便得到拟线性形式。

\partial_t \mathbf{W} + \mathbf{B}(\mathbf{W})\,\partial_x \mathbf{W} = \mathbf{0}

其中 $\rho$ 是密度， $u$ 是速度， $p$ 是压力。关键的观察是：系数矩阵 $\mathbf{B}$ 恰好裂成两块。

\mathbf{B}(\mathbf{W}) = \underbrace{\mathbf{A}(\mathbf{W})}_{\text{声学}} + \underbrace{\mathbf{C}(\mathbf{W})}_{\text{对流}}, \qquad \mathbf{C}(\mathbf{W}) = u\,\mathbf{I}

$\mathbf{A}$ 包含全部压力项，即声学（acoustic）部分。 $\mathbf{C}=u\mathbf{I}$ 只是原样搬运物质，即对流（convective）部分。这种分裂实质上就是把 $u\partial_x$ 项从拉格朗日导数 $\mathrm{D}/\mathrm{D}t = \partial_t + u\partial_x$ 中剥离出来。

特征值的加法分裂

这种分裂为何强大，看特征值便一目了然。全系统的波速是 $u-c,\ u,\ u+c$ 。它们恰好由声学与对流两部分相加而成。

\lambda_{1,2,3} = \underbrace{(-c,\ 0,\ +c)}_{\lambda^a\ \text{声学}} + \underbrace{(u,\ u,\ u)}_{\lambda^c\ \text{对流}}

$c=\sqrt{\gamma p/\rho}$ 是声速。声学波速为 $\pm c$ ，与流速无关。对流波速一律为 $u$ 。在低马赫极限 $M=u/c\to 0$ 下，声学带保持宽度 $2c$ ，而对流速度缩向零。这一落差正是刚性（stiffness）的根源。

在下面的模拟中动手试试。

Mach M = u/c = 0.150

stiffness ratio (u+c)/u = 7.7×

amber = 압력 섭동(음향파, u±c) · teal = 엔트로피 이상(대류/접촉, u). M을 0.02까지 내리면 음향파가 접촉파를 압도적으로 앞질러 달아난다.

把马赫数降到0.02，压力脉冲（amber）裂成两支声波向两侧飞奔，而熵扰动（teal）乘着接触波几乎不动。三个三角标记的速度 $u\!-\!c,\ u,\ u\!+\!c$ 正是上面的加法分裂。

声学步 — 拉格朗日坐标下的HLLC#

声学部分就是压力驱动的膨胀与压缩。论文把它移到质量（拉格朗日）坐标，用HLLC型黎曼求解器求解。界面 $j+1/2$ 的星号状态归结为两个公式。

u^*_{j+1/2} = \frac{u_j+u_{j+1}}{2} + \frac{p_j-p_{j+1}}{2\,a_{j+1/2}}, \qquad p^*_{j+1/2} = \frac{p_j+p_{j+1}}{2} + \frac{a_{j+1/2}}{2}\,(u_j-u_{j+1})

$a_{j+1/2}=\max(\rho_j c_j,\ \rho_{j+1}c_{j+1})$ 是声阻抗（密度乘声速）。仅凭这些星号速度与压力，就能推进Euler变量的更新。每个网格的膨胀率压缩成一个系数。

R_j = 1 + \frac{\Delta t}{\Delta x}\big(u^*_{j+1/2}-u^*_{j-1/2}\big)

$R_j$ 表示声学步把网格体积拉伸或压缩了多少。密度随即得到 $\rho^{n+1-}_j = \rho^n_j / R_j$ 。

对流步 — 迎风投影

现在用声学步产生的星号速度 $u^*$ 来搬运物质。用简单的迎风格式更新守恒量 $\varphi\in\{\rho,\rho u,\rho E\}$ 。

\varphi^{n+1}_j = R_j\,\varphi^{n+1-}_j - \frac{\Delta t}{\Delta x}\big(u^*_{j+1/2}\varphi^{\text{up}}_{j+1/2} - u^*_{j-1/2}\varphi^{\text{up}}_{j-1/2}\big)

当 $u^*\!\ge\!0$ 时 $\varphi^{\text{up}}_{j+1/2}$ 取 $\varphi_j$ ，否则取 $\varphi_{j+1}$ 。上一步的 $R_j$ 在此回归，质量、动量与能量得以精确守恒。依次走完这两步，一个时间步便完成。

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$ 又送了回去。

低马赫下真正分道扬镳的东西

分裂真正发挥价值的地方在低马赫。直接法的稳定时间步总是被 $|u|+c$ 拴住，哪怕关心的物理其实活在 $u$ 上。把这个速度比画成图看看。

M = 0.200

(u+c)/u = 6.0×  ← 직접법 시간전진이 견디는 속도비

음향 밴드(amber) 폭은 M과 무관하게 항상 2c. 대류 속도 u(teal)만 0으로 줄어든다. 이 격차가 저마하 강성의 정체.

声学带的宽度（amber）与 $M$ 无关，恒为 $2c$ 。只有对流速度 $u$ （teal）收敛到零。在 $M=0.05$ 时，比值 $(u+c)/u$ 已达21。分裂格式可以让对流步以不同于声学步的时间步推进，于是慢速物理不再被快速波挟持。低马赫下直接Godunov遭受的过度数值扩散，也能靠这种分裂避开。

批判性地看

有三点值得留意。其一，这里用的分裂在时间上是一阶的。要达到二阶，需用Strang分裂把声学–对流–声学包起来，代价也随之增加。其二，原论文真正的舞台是Kapila五方程两相流。体积分数方程中的非守恒项 $K\partial_x u$ 与正性保证才是真正的难关，本文向单相的化简恰好越过了它们。其三，阻抗 $a=\max(\rho c)$ 稳健却耗散。在强激波上，这种耗散会抹糊接触面，所以论文另行把精度与效率同直接法作对比。

从工程角度看，这个想法并不新奇。OpenFOAM基于压力的 rhoPimpleFoam 以及all-Mach系列求解器，用隐式与显式分别处理声学与对流，根子是一样的。Lagrange–Projection只是用黎曼求解器的语言，把这种分裂写得更清楚的那个版本。

这篇论文改变了什么

分裂的正当性：波速 $u\pm c,\ u$ 恰由声学 $(\pm c,0)$ 与对流 $(u,u,u)$ 相加而成。这个加法是整个分裂的根据。
低马赫的处方：声学带固定为 $2c$ ，对流缩向零。分开推进二者即可绕开刚性。
守恒不是免费的：必须把声学步的 $R_j$ 送回对流步，质量、动量与能量才得以存活。