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

声速仍是每秒340米。旁边有人以每小时5公里散步，也一样。可压缩求解器把两种速度塞进同一个时间步。麻烦在比值。极慢的流动里，声波比物质快上数百倍。直接的Godunov类格式被快波拴住时间步。慢对流则被数值扩散抹平。

把Mach滑条往下拉，压力脉冲裂开的速度和熵团停住的位置一起变。

拉到0.02时，amber脉冲向两侧裂开飞走，teal团几乎不动。三个三角落在 $u\!-\!c$、$u$、$u\!+\!c$。越慢，只有中间那个停住，外侧间距仍是 $2c$。

ten Eikelder等（2017）把两者拆开。声学一份，对流一份，各用自己的格式轮流算。这里把这套Lagrange–Projection放到单相（single-phase）Euler上，用慢速流动里的声学脉冲对数字。

“

论文: 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$，便得到拟线性形式。

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

$\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$。它们恰好由声学与对流两部分相加而成。

$c=\sqrt{\gamma p/\rho}$ 是声速。声学波速为 $\pm c$，与流速无关。对流波速一律为 $u$。在低马赫极限 $M=u/c\to 0$ 下，声学带保持宽度 $2c$，而对流速度缩向零。这一落差正是刚性（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})$ 是声阻抗（密度乘声速）。仅凭这些星号速度与压力，就能推进Euler变量的更新。每个网格的膨胀率压缩成一个系数。

$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\}$。

当 $u^*\!\ge\!0$ 时 $\varphi^{\text{up}}_{j+1/2}$ 取 $\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$。取 $c=1$、$M=0.05$，则 $u=0.05$、$u+c=1.05$，比值 $(u+c)/u$ 为 21。分裂格式可以让对流步用不同于声学步的时间步。慢物理不再被快波拖着走。低马赫下直接 Godunov 的过度数值扩散，也能靠这种分裂躲开。

批判性地看

有三处别扭。时间精度是一阶。要二阶，就得用 Strang 分裂把声学–对流–声学包起来，多一次扫描，代价跟着上去。原论文的舞台是 Kapila 五方程两相流。体积分数方程里的非守恒项 $K\partial_x u$ 和正性保证才是难处，单相化简正好绕过它们。阻抗 $a=\max(\rho c)$ 稳，但耗散。强激波上，这份耗散会把接触面抹厚。论文因此另把精度和效率同直接法对照。

工程上这条分叉并不陌生。OpenFOAM 基于压力的 rhoPimpleFoam 以及 all-Mach 系列求解器，用隐式与显式分别处理声学与对流，根子一样。Lagrange–Projection 是用黎曼求解器的语言把这种分裂写下来的那一版。

漏掉 R_j，守恒先死#

声学步的 $R_j$ 不送回对流投影，质量、动量、能量的更新就会互相错开。上面那次运行误差落到 $2.2\times 10^{-16}$，只发生在这份回送还在的时候。用 Strang 包起来可以提高阶数。顺序拉成声学–对流–声学，扫描也多一次。

走到 Kapila 五方程，体积分数里的 $K\partial_x u$ 还在。缩成单相，就碰不到这个非守恒项。阻抗取 $a=\max(\rho c)$，激波扛得住，接触面变厚。从时间积分里把声学拿开的另一种切法，是 Peluchon的声-输运IMEX。$M=0.05$ 时 $(u+c)/u=21$。这 21 倍交给声学求解器，还是留在对流 CFL 上，是实现时的分叉。