[论文综述] 用两个压力先解,再合成一个 — Saurel(2009) 可压缩多相流松弛法

用六方程压力非平衡 + 刚性松弛绕开五方程扩散界面模型陷阱的 Saurel·Petitpas·Berry(2009) 综述。

爱达荷国家实验室的一次核燃料安全模拟突然停下来。在液态钠之间被困住的小气泡那一格,混合声速骤降到约 25 m/s。这就是假设压力平衡的五方程多相流模型的极限。Richard Saurel 与两位合作者在 2009 年给出退一步的答案 — 让两相先各自携带不同的压力,再在每个时间步末尾强行把它们拉回一致。本文跟随这套六方程非平衡 + 松弛算法为何能跑通,界面压力 $\hat{p}_I$ 的阻抗加权平均从何而来,以及 Wood 与 frozen 两条音速曲线分别在告诉我们什么。

一句要点#

平衡的五方程在声速处崩溃。退一步用六方程求解,再用刚性松弛把两个压力合并。阻抗平均 $\hat{p}_I = (Z_2 p_1 + Z_1 p_2)/(Z_1+Z_2)$ 这一行支撑起整套体系。

Wood 声速是个沙袋#

Kapila(2001) 的五方程模型假设压力平衡:同一格内两相共享单一压力。乍看简洁,但混合声速会跑偏。混合声速 $c_{eq}$ 满足 Wood(1930) 关系

\frac{1}{\rho c_{eq}^2} = \frac{\alpha_1}{\rho_1 c_1^2} + \frac{\alpha_2}{\rho_2 c_2^2}

其中 $\alpha_k$ 是相 $k$ 的体积分数, $\rho_k$ 是密度, $c_k$ 是纯相声速。问题在于这条曲线不是单调的。水(1500 m/s)与空气(340 m/s)在同一格内的平衡声速可以掉到约 24 m/s — 远低于两端任何一个纯相值。

下面的模拟可以亲手操作。拉动滑条改变 $\rho_k$ 、 $c_k$ ,观察两条曲线如何分道扬镳。

ρ₁ = 1000 kg/m³c₁ = 1500 m/sρ₂ = 10.0 kg/m³c₂ = 340 m/scursor α₁ = 0.500log scale on y-axis

Wood (equilibrium) sound speed dips far below either pure-fluid value near intermediate α; frozen speed stays inside the [c₂, c₁] band.

数值扩散区会出现两个纯态之间的虚拟混合,声波在其中突然变慢,波到达时间错位,激波追踪失灵。再加上 $\alpha$ 在接近 0 或 1 时,体积分数的正性也变得难以维持。

绕路 — 让两相压力故意不等#

Saurel 选择退一步。放弃压力平衡假设,让两相各自带 $p_1$ 、 $p_2$ 。从七方程 Baer–Nunziato 模型出发,只先取速度松弛极限,得到的就是六方程模型。

\partial_t \alpha_1 + u\,\partial_x \alpha_1 = \mu(p_1 - p_2)

\partial_t (\alpha_k \rho_k) + \partial_x(\alpha_k \rho_k u) = 0

\partial_t (\rho u) + \partial_x(\rho u^2 + \alpha_1 p_1 + \alpha_2 p_2) = 0

\partial_t (\alpha_k \rho_k e_k) + \partial_x(\alpha_k \rho_k e_k u) + \alpha_k p_k\,\partial_x u = \pm p_I \mu (p_1 - p_2)

其中 $\mu$ 是压力松弛率(刚性), $p_I$ 是界面压力, $e_k$ 是各相内能。该模型的混合声速变为

c_f^2 = Y_1 c_1^2 + Y_2 c_2^2

其中 $Y_k = \alpha_k \rho_k / \rho$ 为质量分数。该曲线单调,Wood 留下的沙袋消失。

界面压力 $\hat{p}_I$ — 阻抗加权平均#

出现在内能方程右端的 $p_I$ 来自七方程模型不对称项在速度平衡极限下的整理结果,落在简洁形式

p_I = \frac{Z_2 p_1 + Z_1 p_2}{Z_1 + Z_2}, \quad Z_k = \rho_k c_k

其中 $Z_k$ 是声学阻抗。更"硬"的相获得更大权重。直观上,当两相通过界面交换功时,硬的一侧动得少,软的一侧动得多 — 比例正好就是阻抗比。这一行同时保证了能量守恒与熵不等式。

松弛步 — 单变量方程一行#

数值上不直接积分 $\mu \to \infty$ 极限。采用算子分裂:先解双曲部分,再单独积分松弛源。松弛过程中 $u$ 、 $\rho E$ 、 $\alpha_k \rho_k$ 保持不变,只有比体积 $v_k = 1/\rho_k$ 演化。配合 stiffened gas EOS

p_k = (\gamma_k - 1)\rho_k e_k - \gamma_k p_{\infty,k}

积分有闭式解。

v_k(p) = v_k^0\,\frac{p^0 + \gamma_k p_{\infty,k} + (\gamma_k - 1)\hat{p}_I}{p + \gamma_k p_{\infty,k} + (\gamma_k - 1)\hat{p}_I}

剩下只需解饱和约束 $\sum_k (\alpha\rho)_k v_k(p) = 1$ 求出唯一的松弛压力 $p$ 。一次二分法即可完成。这是算法核心。

重置步 — 用守恒能量重新对齐#

松弛得到的 $p$ 、 $\alpha_k$ 与各相 EOS 一致,但与守恒的混合能量可能存在偏差。这是因为非守恒内能方程在激波处会累积小误差。Saurel 再做一次修正:从守恒的混合总能量,经混合 EOS

p(\rho, e, \alpha_1, \alpha_2) = \frac{\rho e - \big(\frac{\alpha_1 \gamma_1 p_{\infty,1}}{\gamma_1 - 1} + \frac{\alpha_2 \gamma_2 p_{\infty,2}}{\gamma_2 - 1}\big)}{\frac{\alpha_1}{\gamma_1 - 1} + \frac{\alpha_2}{\gamma_2 - 1}}

重新算压力,再用该压力反求各相 $e_k$ 。模型起步是非守恒的,但在激波处表现得像真正的守恒模型。代价仅为每格一次额外 EOS 评估。

代码 — 一格之内两个压力相遇的瞬间#

仅松弛步,大约 30 行 NumPy。Stiffened gas + 二分法求解饱和方程的核心循环。

import numpy as np
 
def stiffened_gas_volume(p, p0, v0, gamma, p_inf, p_hat_I):
    num = p0 + gamma * p_inf + (gamma - 1.0) * p_hat_I
    den = p   + gamma * p_inf + (gamma - 1.0) * p_hat_I
    return v0 * num / den
 
def saturation_residual(p, alpha_rho, p0, v0, gamma, p_inf, p_hat_I):
    total = 0.0
    for k in range(2):
        v_k = stiffened_gas_volume(p, p0[k], v0[k], gamma[k], p_inf[k], p_hat_I)
        total += alpha_rho[k] * v_k
    return total - 1.0
 
def stiff_pressure_relax(p0, v0, alpha_rho, gamma, p_inf, Z, tol=1e-9):
    p_hat_I = (Z[1] * p0[0] + Z[0] * p0[1]) / (Z[0] + Z[1])
    lo, hi = 1.0, 1.0e10
    for _ in range(80):
        mid = 0.5 * (lo + hi)
        r = saturation_residual(mid, alpha_rho, p0, v0, gamma, p_inf, p_hat_I)
        if r > 0:
            lo = mid       # 压力过低 -> 体积之和大于 1
        else:
            hi = mid
        if hi - lo < tol * mid:
            break
    return 0.5 * (lo + hi), p_hat_I
 
# 水(alpha=0.7) + 空气(alpha=0.3), 两个压力不同
gamma  = np.array([4.4, 1.4])
p_inf  = np.array([6.0e8, 0.0])
rho0   = np.array([1000.0, 1.0])
alpha0 = np.array([0.7, 0.3])
p0     = np.array([5.0e5, 1.0e5])
Z      = np.array([1500.0 * 1000.0, 340.0 * 1.0])
 
alpha_rho = alpha0 * rho0
v0        = 1.0 / rho0
 
p_eq, p_hat_I = stiff_pressure_relax(p0, v0, alpha_rho, gamma, p_inf, Z)
print(f"阻抗平均 ^pI = {p_hat_I:.3e} Pa")
print(f"松弛后压力 p = {p_eq:.3e} Pa")

水的阻抗( $Z_1 = 1.5\times 10^6$ )约为空气( $Z_2 = 340$ )的 4400 倍,因此 $\hat{p}_I$ 实际上几乎等于水的压力。最终的松弛压力是两相相互微量压缩与膨胀的结果。

观察 — 两个压力相遇的时间#

只有一格还看不出动态。在 1D 网格上跑一个简化 ODE 模型,每一格的 $p_1, p_2$ 都向阻抗平均 $\hat{p}_I$ 指数收敛 — $\mu \cdot \Delta t$ 越大越快。

Z₁ = 1.50 MPa·s/mZ₂ = 3.40 kPa·s/mμ·Δt = 0.050p₁(left) = 5.00 barp₂(left) = 1.00 barp₁(right) = 1.00 barp₂(right) = 5.00 bar

Stiff relaxation drives both p₁ and p₂ toward the impedance-weighted average p̂ᴵ = (Z₂p₁ + Z₁p₂)/(Z₁+Z₂). Increase μ·Δt for faster pull-in.

增大 $Z_1$ , $\hat{p}_I$ 会更快被拉向相 1 的压力一侧。阻抗比就是权重。把 $\mu \cdot \Delta t$ 调到 1 附近,一步即可达到平衡 — 这就是刚性松弛的含义。Saurel 论文把这一过程无条件用于每一时间步末尾,从而恢复五方程模型的刚性极限。

局限 — 七方程的影子#

六方程并非万能。速度平衡仍是假设。多孔介质或粉末炸药这类两相速度差异极大的问题,需要回到七方程 Baer–Nunziato。在极强激波下,非守恒内能误差不可忽略,需要 Petitpas 的人工热交换修正(论文第 5 节)。而 $\mu \to \infty$ 假设界面响应无限快;一旦表面张力或质量传递参与,就需要单独建模。

尽管如此,六方程 + 松弛仍是当前压缩性多相模拟最广泛使用的基础。AP-IMEX 时间积分(2026-05-03 文章)与 MUSCL-THINC-BVD 界面重构(2026-05-02 文章)都建立在这一六方程骨架之上。

三个要带走的事实#

五方程的陷阱在声速里。 Wood 平衡音速的非单调性会在界面处误导声波。
六方程 + 松弛在刚性极限下恢复五方程。 不必畏惧非守恒系统 — 解一次,合一次。
阻抗平均 $\hat{p}_I = (Z_2 p_1 + Z_1 p_2)/(Z_1+Z_2)$ 是关键。 能量守恒、熵不等式、单调声速 — 这一行全包了。