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

水与空气同格时 Wood 声速掉到约 24 m/s，阻抗平均贴在空气压力一侧，而不是水。

爱达荷国家实验室的一次核燃料安全模拟突然停下来。液态钠之间困着小气泡的那一格,混合声速骤降到约 25 m/s。这是假设压力平衡的五方程模型的尽头。2009 年 Saurel 与合作者退了一步。两个压力分开解,只在步末贴到一起。

拨动 $\rho_k$ 和 $c_k$ 滑条,Wood 曲线与 frozen 曲线拉开的宽度会变。

ρ₁ = 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.

水 1500 m/s、空气 340 m/s 附近,只有 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 — 远低于两端任何一个纯相值。

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

把代码里那一格代进 Wood 公式,正文的 24 m/s 会再出现。 $\alpha_1 = 0.7$ , $\rho_1 = 1000$ , $c_1 = 1500$ 时, $\rho_1 c_1^2 = 2.25\times 10^9$ 。空气一侧 $\alpha_2 = 0.3$ , $\rho_2 = 1$ , $c_2 = 340$ ,故 $\rho_2 c_2^2 = 1.156\times 10^5$ 。右边空气项约 $2.6\times 10^{-6}$ ,水项约 $3\times 10^{-10}$ ,几乎看不见。混合密度是 $0.7\times 1000 + 0.3 = 700.3$ 。 $1/(\rho c_{eq}^2)$ 被空气项吃掉后, $c_{eq}$ 约 24 m/s。同一格的 frozen 式 $c_f^2 = Y_1 c_1^2 + Y_2 c_2^2$ 相反。质量几乎全在水里, $c_f$ 贴在 1500 m/s 旁边。一格之内两个声速差了大约六十倍。这处塌陷,水和空气各半的声速挖得更深。

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

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 = (Z_2 p_1 + Z_1 p_2)/(Z_1 + Z_2)$ 。大的权重落在小的压力上。代入 $p_1 = 5\times 10^5$ , $p_2 = 1\times 10^5$ ,分子里 $Z_1 p_2$ 为 $1.5\times 10^{11}$ , $Z_2 p_1$ 只有 $1.7\times 10^8$ 。所以 $\hat{p}_I$ 不站在水压旁边,而站在空气压力 $1\times 10^5$ 紧邻,约 $1.001\times 10^5$ Pa。硬的一相并不独占界面。松弛之后的共同压力,要等饱和条件上的二分法再去找。

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

只有一格还看不出动态。在 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 时间积分与 MUSCL-THINC-BVD 重构都跑在这副六方程骨架上。

那一格还留着五方程吗

Wood 声速掉到 24 m/s 的格子若交给五方程,声波到达会拖后,激波位置会晃。那一格退回六方程,两个压力分开带,只在步末用阻抗平均贴上。

速度也对不齐时,这个妥协就到头了。多孔介质和粉末炸药回到七方程。方程再减、界面怎么闭合,是四方程平衡另问的事。