[论文解读] 用界面速度而不是声速来定Δt，23,865步变成了201步

低马赫下让声速决定时间步长，计算量增加100倍，界面也被完全抹平。只把压力在参考态附近线性化并隐式求解，Δt就可以由界面速度决定。

搬运一个气泡为什么要两万步

水中一个气泡以1 m/s漂移。而水中声速是1,500 m/s。全显式可压缩求解器会按声速缩小时间步长。气泡连一个网格都没走完，求解器已经跑了几千步。更糟的是，每一步都会给界面叠加声速量级的数值扩散。

本文解读Lukáčová-Medvid'ová、Peshkov和Thomann于2024年发表在JCP上的论文 An implicit-explicit solver for a two-fluid single-temperature model。核心只有一点：只把压力项在常数参考态附近线性化并隐式处理。这样时间步长就由界面移动的速度决定。读完本文，你会看到这种分裂为何只需解一个线性系统，以及在1D实验中步数和界面厚度如何变化。

两个时钟：声学CFL与物质CFL#

用马赫数 $M$ （流速除以声速）无量纲化后，压力项前面出现 $1/M^2$ 。它产生的波速为 $|u| \pm c/M$ 。全显式格式的稳定条件受最快的波支配。

\Delta t_{\text{ex}} \le \nu \frac{\Delta x}{\max\left(|u| + c/M\right)}

$\nu$ 为CFL数， $c$ 为无量纲声速， $\Delta x$ 为网格间距。

论文采用的条件不同，即式(59)的物质CFL（由物质运动速度确定的时间步长）。

\Delta t_{\text{IMEX}} \le \nu \frac{\Delta x}{\max\left(|\lambda^{v}_1|, |\lambda^{v}_2|\right)}, \qquad \lambda^{v}_1 = \boldsymbol v\cdot\boldsymbol n,\ \ \lambda^{v}_2 = (\boldsymbol v + (1-2c_1)\boldsymbol w)\cdot \boldsymbol n

$\boldsymbol v$ 为混合物速度， $\boldsymbol w$ 为两相相对速度， $c_1$ 为第1相质量分数。声速完全从式中消失了。两个时钟之比约为 $1 + c/M$ ， $M=0.01$ 时约为120倍。

问题不只是时间步长。Rusanov通量（按局部最大波速施加扩散的最简单近似黎曼求解器）的数值扩散为 $\lambda \Delta x/2$ 。一旦 $\lambda$ 里含有 $c/M$ ，体积分数 $\alpha$ 也会以声速量级被抹平。

下面的模拟可以直接动手操作。

Mach M = 0.010 cells N = 100 slow-motion ×6

explicit: Δt = νΔx/(|u|+c/M) = 4.19e-5   steps so far 0   steps to t=1 ≈ 23860
IMEX:     Δt = νΔx/|u|       = 5.00e-3   steps so far 0   steps to t=1 ≈ 200
ratio of step counts ≈ 1 + c/M = 119.3

把Mach滑块往左拖，橙色压力脉冲会糊满整根管道，蓝色刻度变得密密麻麻。绿色气泡则保持原来的速度前进。下方面板中， $M$ 越小，蓝色 $\alpha$ 曲线越快变平。

隐式处理什么：压力与化学势

论文的模型是单温度双流体模型，在对称双曲热力学相容（SHTC）框架下推导。未知量为 $\alpha_1$ 、各相质量 $\alpha_l\rho_l$ 、混合动量 $\rho\boldsymbol v$ 、相对速度 $\boldsymbol w$ 和总能量 $\rho E$ 。低马赫极限下的刚性（时间尺度相差悬殊、难以显式求解的性质）来自两处：动量方程中的压力 $p$ ，以及相对速度方程中的化学势差 $\mu$ 。

好在压力对 $\rho E$ 是线性的（式42）。

p = (\phi_p - 1)\left(\rho E - \rho E_{\text{kin}}\right), \qquad \phi_p = \frac{\gamma_1\alpha_1\rho_1 c_{v,1} + \gamma_2\alpha_2\rho_2 c_{v,2}}{\alpha_1\rho_1 c_{v,1} + \alpha_2\rho_2 c_{v,2}}

$\gamma_l$ 为比热比， $c_{v,l}$ 为定容比热， $E_{\text{kin}}$ 为动能项。单相时 $\phi_p=\gamma-1$ ，退化为欧拉方程。

麻烦在于 $\mu$ 。它通过相熵非线性地依赖于 $\rho E$ 。若直接隐式求解，每一步都要解非线性方程组。

参考态线性化：剩下的非线性项是 $O(M^2)$ #

论文给出的方法是参考态（reference state，RS）法。在渐近分析给出的常数主阶状态 $\boldsymbol q_{\text{RS}}$ 附近对 $\mu$ 做泰勒展开（式45、48）。

\mu = \hat\mu + \bar\mu, \qquad \hat\mu = \hat\mu_{\text{RS}} + \left(\frac{\partial \mu}{\partial(\rho E)}\right)_{\text{RS}} \rho E

$\hat\mu$ 是对 $\rho E$ 线性的部分， $\bar\mu$ 是余项。论文证明，对于良好准备的初值（压力扰动为 $O(M^2)$ 的初始条件）， $\bar\mu = O(M^2)$ （式53）。低马赫下这个余项很小，因此显式处理也不会限制时间步长。

这样分裂的结果很干净。显式子系统(54)只剩输运，特征速度只有 $\lambda^v_1$ 和 $\lambda^v_2$ 。隐式子系统(57)成为常系数线性系统。参考态是常数，所以系数也是常数。无论直接法还是迭代法，解一次即可。

时间积分采用刚性精确（stiffly accurate）的IMEX Runge–Kutta格式，压力项用中心差分处理。论文证明这两点保证了渐近保持性（AP，即 $M\to0$ 极限下与不可压离散一致的性质）。若想把类似的分裂与双时间预处理对比，可以参考把稳态预处理放进双时间循环的实验。

代码：1D气泡输运，三种组合#

论文5.3节在周期区域中让气泡以 $\boldsymbol v=(1,1)$ 绕行一周。这里把它降到1D，模型也做了简化：等熵单相欧拉方程 $p=\rho^\gamma$ ，再把 $\alpha_1$ 作为标记函数带上。这是一个用 $\rho$ 代替论文中 $\rho E$ 的简化版。参考态取 $\rho_{\text{RS}}=1$ 。

\hat p = 1 + \gamma(\rho - 1), \qquad \bar p = \rho^\gamma - \hat p = O\left((\rho-1)^2\right)

$\hat p$ 隐式处理， $\bar p$ 显式处理。把面动量写成 $\rho^{n+1}$ 的函数并代入质量方程，就只剩下一个方程。

\left[1 + \beta\left(2 - S^{+} - S^{-}\right)\right]\rho^{n+1} = \text{rhs}, \qquad \beta = \frac{\gamma\,\Delta t^2}{M^2\Delta x^2}

$S^{\pm}$ 为单网格平移算子。这是常系数循环矩阵，一次FFT即可求解。比较三种组合：(A) 全显式＋声学CFL，(B) 全显式＋物质CFL，(C) RS-IMEX＋物质CFL。均取 $\nu=0.5$ 、 $N=100$ 、 $t=1$ （一周）。

import numpy as np
np.seterr(all='ignore')
 
N, GAM, NU, TEND = 100, 1.4, 0.5, 1.0
dx = 1.0 / N
x = (np.arange(N) + 0.5) * dx
 
def bubble_init(M):
    """1D analog of the paper's bubble test: u=1, p in equilibrium, sharp alpha."""
    r = np.abs(x - 0.5)
    alpha = 0.8 * np.arctan(-2000.0 * (r - 0.2)) / np.pi + 0.5
    rho = 1.0 + 0.1 * M**2 * np.cos(2 * np.pi * x)   # well-prepared O(M^2) perturbation
    return alpha, rho, rho * 1.0
 
def ddx_c(f):  # centered derivative, periodic
    return (np.roll(f, -1) - np.roll(f, 1)) / (2 * dx)
 
def lap(f):
    return (np.roll(f, -1) - 2 * f + np.roll(f, 1)) / dx**2
 
def alpha_update(alpha, u, lam, dt):
    """Rusanov form of alpha_t + u alpha_x = 0 with dissipation speed lam."""
    adv = u * ddx_c(alpha)
    diss = 0.5 * lam * lap(alpha) * dx
    return alpha - dt * (adv - diss)
 
def rusanov_full(rho, m, M, dt):
    u = m / rho
    p = rho**GAM
    lam = np.max(np.abs(u) + np.sqrt(GAM * p / rho) / M)
    U = np.stack([rho, m])
    F = np.stack([m, m * u + p / M**2])
    Up, Fp = np.roll(U, -1, axis=1), np.roll(F, -1, axis=1)
    H = 0.5 * (F + Fp) - 0.5 * lam * (Up - U)       # face i+1/2
    U = U - dt / dx * (H - np.roll(H, 1, axis=1))
    return U[0], U[1], lam
 
def rs_imex_step(rho, m, M, dt):
    """Explicit: advection + nonlinear pressure remainder. Implicit: linearized pressure."""
    u = m / rho
    lam_u = np.max(np.abs(u))
    pbar = rho**GAM - (1.0 + GAM * (rho - 1.0))      # p - p_hat, O((rho-rho_RS)^2)
    G = m * u + pbar / M**2
    H = 0.5 * (G + np.roll(G, -1)) - 0.5 * lam_u * (np.roll(m, -1) - m)
    m_star = m - dt / dx * (H - np.roll(H, 1))
    # face momentum: m_f = avg(m*) - dt*GAM/M^2 * drho/dx ;  rho^{n+1} = rho - dt*div(m_f)
    beta = dt**2 * GAM / M**2 / dx**2
    mf_star = 0.5 * (m_star + np.roll(m_star, -1))
    rhs = rho - dt / dx * (mf_star - np.roll(mf_star, 1))
    k = np.fft.fftfreq(N) * 2 * np.pi
    rho_new = np.real(np.fft.ifft(np.fft.fft(rhs) / (1 + beta * (2 - 2 * np.cos(k)))))
    m_new = m_star - dt * GAM / M**2 * ddx_c(rho_new)
    return rho_new, m_new, lam_u
 
def run(M, scheme, dt_mode):
    alpha, rho, m = bubble_init(M)
    a0 = alpha.copy()
    t, n = 0.0, 0
    while t < TEND - 1e-12:
        u = m / rho
        c = np.sqrt(GAM * rho**(GAM - 1))
        speed = np.max(np.abs(u) + c / M) if dt_mode == "acoustic" else np.max(np.abs(u))
        dt = min(NU * dx / speed, TEND - t)
        if scheme == "explicit":
            rho, m, lam = rusanov_full(rho, m, M, dt)
        else:
            rho, m, lam = rs_imex_step(rho, m, M, dt)
        alpha = alpha_update(alpha, u, lam, dt)
        t += dt; n += 1
        if not np.isfinite(rho).all() or np.max(np.abs(rho - 1.0)) > 0.5:
            return n, t
    err = np.sum(np.abs(alpha - a0)) * dx
    width = np.count_nonzero((alpha > 0.15) & (alpha < 0.85))
    du = np.max(np.abs(m / rho - 1.0))
    return n, (err, width, du, alpha.max())
 
print(f"N={N}, nu={NU}, T=1 (one period), initial interface width = 2 cells")
print(f"{'M':>6} {'case':<26} {'steps':>7} {'L1(alpha)':>10} {'width':>6} {'max a':>6} {'max|u-1|':>9}")
for M in [0.1, 0.03, 0.01]:
    for label, sch, mode in [("explicit, acoustic CFL", "explicit", "acoustic"),
                             ("explicit, material CFL", "explicit", "material"),
                             ("RS-IMEX, material CFL", "imex", "material")]:
        n, res = run(M, sch, mode)
        if not isinstance(res, tuple):
            print(f"{M:>6} {label:<26} {n:>7}  blew up at t={res:.3f}")
        else:
            err, w, du, am = res
            print(f"{M:>6} {label:<26} {n:>7} {err:>10.4f} {w:>6d} {am:>6.3f} {du:>9.1e}")

N=100, nu=0.5, T=1 (one period), initial interface width = 2 cells
     M case                         steps  L1(alpha)  width  max a  max|u-1|
   0.1 explicit, acoustic CFL        2568     0.3580    100  0.458   2.5e-03
   0.1 explicit, material CFL          16  blew up at t=0.074
   0.1 RS-IMEX, material CFL          201     0.0886     44  0.895   2.9e-08
  0.03 explicit, acoustic CFL        8089     0.3813    100  0.420   2.2e-05
  0.03 explicit, material CFL          11  blew up at t=0.048
  0.03 RS-IMEX, material CFL          201     0.0885     44  0.895   2.0e-13
  0.01 explicit, acoustic CFL       23865     0.3814    100  0.420   6.9e-09
  0.01 explicit, material CFL           9  blew up at t=0.031
  0.01 RS-IMEX, material CFL          201     0.0885     44  0.895   1.8e-12

下面在浏览器里直接运行上面的代码。

Mach M = 0.010 cells N = 100 CFL ν = 0.50

t = 0.000 / 1   (Step = advance physical time by 0.005)
A  explicit, acoustic CFL  steps      0  L1(α)=0.0000  width=0  max α=0.000  max|u−1|=0.0e+0
B  explicit, material CFL  steps      0  L1(α)=0.0000  width=0  max α=0.000  max|u−1|=0.0e+0
C  RS-IMEX,  material CFL  steps      0  L1(α)=0.0000  width=0  max α=0.000  max|u−1|=0.0e+0

把Mach从0.1降到0.003。A的步数随 $1/M$ 增长，C始终停在201步。B在任何 $M$ 下都会在几十步内发散。把ν调到0.9，还能看到C的界面宽度（width）变窄。

23,865与201，以及0.38与0.09#

从输出中可以读出三点。

第一是步数。A在 $M=0.1, 0.03, 0.01$ 时分别需要2,568 → 8,089 → 23,865步，比例恰好符合 $1+c/M$ 。C与 $M$ 无关，都是201步。 $M=0.01$ 时相差119倍。

第二是稳定性。B与C使用相同的时间步长，但密度在 $t=0.03\sim0.07$ 之间发散。声学项保持显式却用物质CFL是行不通的。能放大时间步长，功劳在于隐式压力求解。

第三是界面。A的 $\alpha$ 最大值从0.9降到0.42，而0.42正是区域平均值。气泡已完全混合消失，L1误差0.38就是结果。原因不在于步数多，而在于每一步都附加了 $(c/M)\Delta x/2$ 量级的扩散。总扩散量与 $\lambda\cdot t$ 成正比，所以 $M$ 越小越糟。C保住了0.895的最大值。不过由于是一阶格式，界面从2个网格扩展到44个。论文也写道RS-IMEX1"相当耗散"，二阶扩展（3.4节）能较好地保持初始形状。

再看速度扰动 $\max|u-1|$ 。C在 $3\times10^{-8}$ 以下。A在 $M=0.1$ 时为 $2.5\times10^{-3}$ ，声波仍然存在。低马赫下没有必要专门去解析这些声波。

这个1D实验看不到的东西#

简化版必然有缺失。真实模型中还有相对速度 $\boldsymbol w$ 和摩擦弛豫 $\tau^{(w)}$ 。论文在气泡算例中需要取 $\tau^{(w)}=10^{-8}\sim10^{-12}$ 才能压住 $\boldsymbol w$ 的扰动。两相马赫数也可能不同。论文算例中 $M_1=1.336$ 与 $M_2=2.67\times10^{-2}$ 共存于同一区域。此时混合马赫数在1.23到0.36之间，论文指出它并不能很好地反映各相的流动区域。

还有无扰动问题。使用混合状态方程时，等速界面上可能出现压力振荡。昨天讨论的基于压力的BN模型的无扰动条件正是这个话题。六方程双流体模型的双曲性在two-fluid模型复特征值一文中讨论过。论文的SHTC模型是对称双曲的，从而避开了这个问题。

在低马赫多相代码中加入IMEX之前要量的两个数#

先计算 $1 + c/M$ 。若超过10，显式格式会同时失去计算时间和界面锐度。再把压力在参考态附近展开，测量余项的大小，确认它是否按 $O(M^2)$ 减小。如果是，这一项就可以保留显式处理，而隐式部分只需求解一个常系数线性系统。

[论文解读] 用界面速度而不是声速来定Δt，23,865步变成了201步 — 参考态IMEX双流体求解器

搬运一个气泡为什么要两万步

两个时钟：声学CFL与物质CFL#

隐式处理什么：压力与化学势

参考态线性化：剩下的非线性项是 $O(M^2)$ #

代码：1D气泡输运，三种组合#

23,865与201，以及0.38与0.09#

这个1D实验看不到的东西#

在低马赫多相代码中加入IMEX之前要量的两个数#

评论

相关文章

[论文解读] 用界面速度而不是声速来定Δt，23,865步变成了201步 — 参考态IMEX双流体求解器

搬运一个气泡为什么要两万步

两个时钟：声学CFL与物质CFL#

隐式处理什么：压力与化学势

参考态线性化：剩下的非线性项是O(M2)O(M^2)O(M2)#

代码：1D气泡输运，三种组合#

23,865与201，以及0.38与0.09#

这个1D实验看不到的东西#

在低马赫多相代码中加入IMEX之前要量的两个数#

评论

相关文章

参考态线性化：剩下的非线性项是 $O(M^2)$ #