本文译自作者的韩语原文。 阅读原文(한국어)
声速1625 m/s的水在压力子系统里只有18 m/s — Toro–Vázquez分裂把刚性移到了哪里
分裂后压力子系统的波速A不是声速。水的A/c只有0.011，消失的刚性转移到了两步之间的状态方程重算中。
声速1625 m/s的水被算成了18 m/s#
水中声速约为1625 m/s。但在这篇论文的压力子系统里，1 bar水的波速是18.44 m/s。这不是算错了。它直接来自一个公式。
论文是Schropff、Petitpas和Daniel于2026年发表在JCP上的"Semi-implicit method for two-phase flows in pressure disequilibrium"。这是一种IMEX(隐式–显式混合)方法，用一个速度求解压力不同的两相。把材料切到Water或Spinel，橙线几乎竖直；压力和马赫数滑块会同时推动 和时间步长比。
灰色虚线是原系统的特征线 。橙色实线是压力子系统的 ，蓝色实线是输运速度 。右侧柱状图用对数刻度显示低马赫下时间步长增大了多少倍。水里橙线越接近竖直， 就越落到0.01以下。
1993年：把压力从对流中拆出来#
起点是Zha–Bilgen分裂(flux splitting，按物理作用拆分通量并分别处理的方法)。它把欧拉方程的通量分成两块。一块是随流速输运的对流部分。另一块是由压力推动的部分。
Toro和Vázquez在2012年把这种分裂改写成Godunov形式。这样每个子系统都能单独提出黎曼问题。Schropff等人把这一思路搬到压力非平衡的两相模型上。与压力效应有关的通量和非守恒项归入第一个子系统。剩下的输运项归入第二个子系统。
2026年：两相拆成两个子系统#
原始变量为 。压力子系统如下。
是相 的体积分数， 是相密度， 是比内能。 是混合压力。
输运子系统里没有压力。体积分数、质量、内能和动量都由 输运。因此7个特征值全部等于 。两个子系统相加即恢复原模型。分析中略去了松弛项。
压力步中密度是冻结的
再看压力子系统。 和 都不随时间变化，所以 也不变。改变压力的唯一途径是内能。写成原始变量形式：
是Grüneisen系数(密度固定时加入能量引起的压力升高率)。 是密度固定时的 。将此式与动量方程联立，得到波速 。
其中 是混合密度。 中不含 。因此论文说这种波系总是亚声速的。
对刚性气体(stiffened gas)状态方程，。于是得到论文式(14)。
是冻结声速(frozen sound speed，无相间松弛时的混合声速)， 是刚性气体的刚性压力。理想气体有 ，空气的 。水的情况不同。 项几乎占了 的全部，这一部分被整个去掉了。
把四种材料的A和c放进一张表#
用论文表1的刚性气体参数(空气、水、环氧树脂、尖晶石)计算 和 。同一段代码还单独取出空气的压力子系统做脉冲实验。在相同Courant数下比较显式Lax–Friedrichs和隐式后向欧拉。最后一块是低马赫下的时间步长比。
import math
# Stiffened-gas parameters (gamma, p_inf [Pa]) from the paper's Table 1
MATS = {"Air": (1.4, 0.0), "Water": (4.4, 6e8), "Epoxy": (2.43, 5.3e9), "Spinel": (1.62, 141e9)}
RHO = {"Air": 1.2, "Water": 1000.0, "Epoxy": 1185.0, "Spinel": 3622.0}
def split_speeds(gam, pinf, p, rho):
c = math.sqrt(gam * (p + pinf) / rho) # full sound speed
A = math.sqrt((gam - 1.0) * p / rho) # pressure-subsystem speed, Eq. (14)
return c, A
def thomas_solve(rhs, beta):
# (I - beta*L) x = rhs, compact Laplacian, zero-gradient ends
n = len(rhs)
cp, dp = [0.0] * n, [0.0] * n
b0 = 1 + beta
cp[0], dp[0] = -beta / b0, rhs[0] / b0
for i in range(1, n):
b = 1 + beta if i == n - 1 else 1 + 2 * beta
m = b + beta * cp[i - 1]
cp[i] = -beta / m
dp[i] = (rhs[i] + beta * dp[i - 1]) / m
x = [0.0] * n
x[-1] = dp[-1]
for i in range(n - 2, -1, -1):
x[i] = dp[i] - cp[i] * x[i + 1]
return x
def ddx(f, dx):
n = len(f)
return [(f[min(i + 1, n - 1)] - f[max(i - 1, 0)]) / (2 * dx) for i in range(n)]
def pulse_run(courant, implicit, n=200, travel=0.3):
# linearised pressure subsystem (air, 1 bar): p_t + K u_x = 0, u_t + p_x / rho = 0
rho, K = 1.2, 0.4 * 1e5
A = math.sqrt(K / rho)
dx = 1.0 / n
x = [(i + 0.5) * dx for i in range(n)]
p = [100.0 * math.exp(-((xi - 0.5) / 0.04) ** 2) for xi in x]
u = [0.0] * n
dt = courant * dx / A
steps = math.ceil(travel / A / dt)
for _ in range(steps):
if implicit:
ux = ddx(u, dx)
p = thomas_solve([p[i] - dt * K * ux[i] for i in range(n)], dt * dt * K / rho / dx ** 2)
px = ddx(p, dx)
u = [u[i] - dt / rho * px[i] for i in range(n)]
else: # Lax-Friedrichs
P = [p[0]] + p + [p[-1]]; U = [u[0]] + u + [u[-1]]
p = [0.5 * (P[i] + P[i + 2]) - dt * K * (U[i + 2] - U[i]) / (2 * dx) for i in range(n)]
u = [0.5 * (U[i] + U[i + 2]) - dt / rho * (P[i + 2] - P[i]) / (2 * dx) for i in range(n)]
if max(abs(v) for v in p) > 1e9:
return steps, float("inf"), float("nan")
right = max(range(n // 2, n), key=lambda i: p[i])
speed = (x[right] - 0.5) / (steps * dt)
return steps, max(abs(v) for v in p[n // 2:]), speed
print("material p[bar] c[m/s] A[m/s] A/c")
for name, (g, pinf) in MATS.items():
for pbar in (1.0, 100.0):
rho = RHO[name] * (pbar if name == "Air" else 1.0) # isothermal air
c, A = split_speeds(g, pinf, pbar * 1e5, rho)
print(f"{name:7s} {pbar:7.0f} {c:10.1f} {A:10.2f} {A / c:9.5f}")
print("\nwater/air mixture, 1 bar, Eq. (13)")
for aw in (0.01, 0.5, 0.99):
rho = aw * 1000.0 + (1 - aw) * 1.2
A = math.sqrt((aw * 1e5 * 3.4 + (1 - aw) * 1e5 * 0.4) / rho)
print(f"alpha_w={aw:4.2f} rho={rho:7.1f} A={A:7.2f} m/s")
A_air = math.sqrt(0.4e5 / 1.2)
print(f"\npressure pulse, air 1 bar, A = {A_air:.1f} m/s, c = {math.sqrt(1.4e5 / 1.2):.1f} m/s")
print("Courant_A steps exp_peak exp_speed imp_peak imp_speed")
for cr in (0.5, 0.9, 1.5, 5.0, 20.0):
s, pe, ve = pulse_run(cr, False)
_, pi, vi = pulse_run(cr, True)
print(f"{cr:7.1f} {s:6d} {pe:9.2f} {ve:9.1f} {pi:9.2f} {vi:9.1f}")
print("\nlow-Mach air: dt ratio at equal Courant number = (|u|+c)/|u|")
for M in (0.3, 0.1, 0.01, 0.001):
print(f"M={M:6.3f} ratio={(1 + M) / M:9.1f}")material p[bar] c[m/s] A[m/s] A/c
Air 1 341.6 182.57 0.53452
Air 100 341.6 182.57 0.53452
Water 1 1624.9 18.44 0.01135
Water 100 1638.3 184.39 0.11255
Epoxy 1 3296.7 10.99 0.00333
Epoxy 100 3299.8 109.85 0.03329
Spinel 1 7941.3 4.14 0.00052
Spinel 100 7941.6 41.37 0.00521
water/air mixture, 1 bar, Eq. (13)
alpha_w=0.01 rho= 11.2 A= 62.00 m/s
alpha_w=0.50 rho= 500.6 A= 19.48 m/s
alpha_w=0.99 rho= 990.0 A= 18.45 m/s
pressure pulse, air 1 bar, A = 182.6 m/s, c = 341.6 m/s
Courant_A steps exp_peak exp_speed imp_peak imp_speed
0.5 120 25.59 184.1 35.90 181.1
0.9 67 42.37 183.2 30.52 180.2
1.5 40 240.14 129.3 25.76 178.0
5.0 12 102438.55 13.7 16.02 168.9
20.0 3 3336.76 32.0 9.82 129.3
low-Mach air: dt ratio at equal Courant number = (|u|+c)/|u|
M= 0.300 ratio= 4.3
M= 0.100 ratio= 11.0
M= 0.010 ratio= 101.0
M= 0.001 ratio= 1001.0把 Courant 滑块从 0.5 推到 20，并在显式和隐式之间切换，表里的 exp_peak 和 imp_peak 会在画面上重算。
到 0.9 为止两种波速都在 180 m/s 附近。超过 1.5，显式峰值跳到 240 Pa。隐式在 20 也只是被压到 9.82 Pa，数字并不炸开。
1 bar水的 为0.01135。环氧树脂降到0.00333，尖晶石降到0.00052。把压力提高到100 bar，水的 正好变为10倍，因为 与 成正比。而 只从1624.9变到1638.3，变化0.8%。空气无论压力如何都固定在0.53452。
水–空气混合物也值得看一眼。水体积分数为0.5时， 为19.48 m/s。与0.99时的18.45 m/s几乎相同。因为混合密度和分子 都由水决定。只有空气占到99%时， 才跳到62.00 m/s。
压力从 1 bar 提到 100 bar，水的 从 18.44 变成 184.39。正好 10 倍。 与 成正比，压力变为 100 倍，波速就是 10 倍。同一段里 只从 1624.9 到 1638.3，也就是页上写的 0.8%。环氧树脂的 同样从 10.99 到 109.85，是 10 倍， 停在 3296.7 和 3299.8。尖晶石的 从 4.14 到 41.37， 是 7941.3 和 7941.6。 从 0.00052 增到 0.00521，仍低于 0.01。
空气在 1 bar 和 100 bar 是同一行。，，。理想气体把密度按压力成比例放大，温度固定， 也就固定。所以 写成 0.535。
只要有一点水， 就贴到水的数值上。 时 ，。 时 ，。 时 ，，和纯水的 18.44 只差 0.01 m/s。分母里的混合密度和分子里的 都由水带走。脉冲表里显式峰值在 Courant 0.9 是 42.37 Pa，1.5 是 240.14 Pa，5.0 是 102438.55 Pa。隐式在同样三格里反而下降：30.52、25.76、16.02 Pa。Courant 20 的隐式峰值 9.82 Pa 是右行波基准 50 Pa 的 1/5，测得速度掉到 129.3 m/s。
刚性没有消失，而是移到了两步之间
那么水的1625 m/s去哪了？两个子系统里都没有 。压力子系统只有 ，输运子系统只有 。
真实声速在两步交替推进时重新出现。输运步中 发生变化。状态方程以 的比率把这一变化转换为压力变化。这个压力成为下一个压力步的输入。水的刚性从子系统特征值中退出，转移到了步与步之间的状态方程重算中。把声学放进一步、输运放进下一步来搬走刚性，把声学和输运分成两步的IMEX也是这么做的。
以下是笔者的解读，并非论文原文。单看压力步，似乎可以按 显式处理。但这个不等式并不保证组合格式的稳定性。从这个角度看，论文即使在 很小的水中也一律隐式求解压力步，是顺理成章的。
隐式压力步让脉冲存活，但被压平
输出的第二块是脉冲实验。显式Lax–Friedrichs在Courant 0.9时还能撑住，到1.5就崩溃。1.5时最大压力达到240 Pa，超过了初始值。5.0时超过10万Pa。隐式格式在20时也不发散。代价是峰值被压到9.82 Pa，只有精确解50 Pa的1/5。
这种衰减正是IMEX的交换条件。低马赫问题中，声波振幅不是关注对象。只要压力场能立刻跟上速度场就够了。若问题需要精确追踪声波，这笔交易就不划算。
论文在隐式压力步中为每一相建立压力方程。由于 是压力的函数，方程是非线性的。用Picard迭代(把部分系数固定为上一次迭代值以实现线性化的迭代法)求解。空间离散采用紧致模板，因此矩阵是窄带矩阵。上面代码中的 thomas_solve 就相当于在1D中求解一个这样的带状矩阵。
低马赫下Δt变成1001倍的算术#
显式格式的时间步长受 限制。IMEX中只有输运步是显式的，因此放宽到 。论文使用 。假设Courant数相同，比值为：
是马赫数。 时为4.3倍， 时为1001倍。不过这是把每步成本视为相同的上限。隐式步每次Picard迭代都要解一个线性方程组。论文结论也这样说。低马赫问题收益很大，激波管问题收益很小甚至没有。
收回非守恒项损失能量的最后一步
分裂还附带一项善后工作。各相内能方程中有 这类非守恒项。离散后，这些项不能精确守恒总能量。
因此论文把修正算子作为第三步。质量、动量、混合总能量和体积分数直接采用预测值。只对各相内能加上 。这样把差值分配下去，使之与混合总能量一致。压力松弛在IMEX步完成后单独进行。
只看水的 就把压力步做成显式，组合格式撑不住#
尖晶石的 是 0.00052，环氧树脂是 0.00333，1 bar 的水是 0.01135。三个都在 0.01 附近或更低。只看压力步， 像是可以显式。这个不等式只讲子系统。组合格式的刚性在 ，不在 。1 bar 的水里 18.44 和 1624.9 是两个数。 时时间步比是 4.3 倍， 时是 1001 倍。在 1D 里，只有 超过 Picard 迭代的成本，这个比才算收益。
Courant 20 的后向欧拉把这份输出里的脉冲峰值削到 9.82 Pa，是 50 Pa 的 1/5。声波振幅进结果的话，这个数字得写下来。隐式压力步把 峰值压平的情形，和隐式压力步把 峰值压平是同一件事。振幅要拿来分析，时间步就按声学来定，或者离开后向欧拉，改用更高阶的隐式积分。
相关文章
如果对您有帮助,请分享。