[论文解读] 水和空气各占一半,声音只剩24 m/s — Murrone–Guillard 五方程模型中的K∇·u项

1930年教科书里的一行公式#

1930年,英国海军研究实验室的A. B. Wood出版了一本声学教科书。书中有一行公式,给出含气泡水的声速。它说混合物的压缩率等于各相压缩率的体积加权和。按此计算,水和空气各占一半的泡沫,声速约为24 m/s。这比水(1625 m/s)慢得多,也比空气(342 m/s)慢得多。

本文解读Murrone–Guillard的INRIA报告RR-4778(2003)。报告研究可压缩两相流的五方程约化模型(reduced model,即由七方程模型在瞬时松弛极限下得到的模型)。本文只问一个问题:体积分数方程中的$K\nabla\cdot u$项到底起什么作用?我们用Python一维计算来回答。同一根管、同一个脉冲,分别加入和去掉$K$,测量声音的速度。

先锁住速度、再靠压力松弛减少方程的路线,写在Saurel(2009)的压力松弛里。本文测量的是把那次松弛推到底之后,体积分数方程里留下的$K$。

α滑块一动,A管的Wood脉冲和B管的快脉冲就在同一画面上分开。

把α滑块放在0.5,A管的脉冲缓慢爬行,B管的脉冲却绕了几十圈。把α移到1%附近,蓝色曲线已经降到水声速的1/14。

2003年报告重新提出的问题#

Wood公式本身是老结论。真正的问题是数值模型能否自行再现它。七方程Baer–Nunziato模型为两相分别保留速度和压力。物理上很完整,但计算代价大。如果两相很快锁定到同一速度和同一压力,就可以减少变量。

Kapila等人(2001)在这一极限下得到了五方程模型。Murrone和Guillard用渐近展开(以小参数将解展开为级数的方法)重新推导了它。随后他们分析了双曲性(具有实数且有限波速的性质)和波结构,并测试了两种数值格式。再往下砍、只留界面平衡,是四方程模型走的那条路。本文保持五个方程,只开关$K$这一行。

报告中的一张图是本文的起点。图中画的是空气–水混合物声速随空气体积分数的变化。共有两条曲线。"5 Eq. Reduc"是约化模型,"5 Eq. Trans"是只输运体积分数的模型。两条曲线相差很大。

松弛无限快时剩下的五个方程

将七方程模型的速度松弛和压力松弛时间尺度取为0,剩下的是:

$\alpha_k$是相$k$的体积分数,$\rho_k$是相密度,$u$和$p$是共同的速度和压力,$E$是混合物总能量。前四个方程是守恒型。只有最后一个是非守恒型,右端的$K\nabla\cdot u$就是本文的主角。

K∇·u从何而来#

压力平衡时,两相经历相同的压力变化。在没有激波的区域,相$k$被等熵压缩,因此$D\ln\rho_k/Dt = (\rho_k c_k^2)^{-1} Dp/Dt$。再结合相质量守恒$D(\alpha_k\rho_k)/Dt = -\alpha_k\rho_k\nabla\cdot u$:

把两相的方程相加。由于$\alpha_1+\alpha_2=1$,左端为0,于是得到压力方程:

$c_W$就是Wood声速。代回第一个式子,$K$随之确定:

$\rho_k c_k^2$是相$k$的体积模量(bulk modulus,抵抗压缩的程度)。水比空气硬约1.9万倍。因此若以空气为相1,则$K\approx\alpha_1$。也就是说,混合物被压缩时,体积变化几乎全部由空气承担。

Wood声速:压缩率按体积,惯性按密度#

Wood声速之所以低,是因为两个性质来自不同的相。压缩率由柔软的空气决定。惯性(密度)由沉重的水决定。含气50%时,压缩率约为纯空气的一半。密度约为空气的400倍。声速按$\sqrt{1/(\text{密度}\times\text{压缩率})}$计算,因此降到空气声速的$1/\sqrt{200}$左右。

令$K=0$,情况就变了。体积分数只是随流体输运。此时混合物表现得像单一的stiffened gas(近似水这类坚硬液体的状态方程)。其声速由$\rho\xi c^2 = \sum_k \alpha_k\xi_k\rho_k c_k^2$给出,其中$\xi_k = 1/(\gamma_k-1)$。水巨大的体积模量通过能量加权混了进来。所以即使含气一半,声音仍然很快。

代码:同一根管里放两个模型

状态方程取空气$\gamma=1.4,\ \pi=0$,水$\gamma=4.4,\ \pi=6\times10^8$ Pa。在1 m长的管中充满含气50%的混合物。在中央放一个振幅0.1%的高斯压力脉冲。计算到右行脉冲应走过0.3 m为止,由峰值位置测速度。守恒变量用Rusanov通量,体积分数用中心差分加同等大小的人工扩散。

import math
 
# stiffened gas: 相 1 = 空气, 相 2 = 水
G = (1.4, 4.4)
PI = (0.0, 6.0e8)
P0 = 1.0e5
RHO0 = (1.2, 1000.0)
 
def wood_speed(a1):
    rho = a1 * RHO0[0] + (1 - a1) * RHO0[1]
    comp = a1 / (G[0] * (P0 + PI[0])) + (1 - a1) / (G[1] * (P0 + PI[1]))
    return 1.0 / math.sqrt(rho * comp)
 
def isobaric_speed(a1):
    """只输运 alpha 时 (K = 0): 混合物就是一种 stiffened gas"""
    rho = a1 * RHO0[0] + (1 - a1) * RHO0[1]
    xi = a1 / (G[0] - 1) + (1 - a1) / (G[1] - 1)
    gpi = a1 * G[0] * PI[0] / (G[0] - 1) + (1 - a1) * G[1] * PI[1] / (G[1] - 1)
    gam = 1 + 1 / xi
    return math.sqrt(gam * (P0 + gpi / xi / gam) / rho)
 
def mix_pressure(a1, rhoe):
    a2 = 1 - a1
    num = rhoe - a1 * G[0] * PI[0] / (G[0] - 1) - a2 * G[1] * PI[1] / (G[1] - 1)
    return num / (a1 / (G[0] - 1) + a2 / (G[1] - 1))
 
def kapila_K(a1, p):
    k1, k2 = G[0] * (p + PI[0]), G[1] * (p + PI[1])
    a2 = 1 - a1
    return a1 * a2 * (k2 - k1) / (a1 * k2 + a2 * k1)
 
def pulse_speed(a_bg, use_K, n=400, travel=0.3, amp=1e-3):
    dx = 1.0 / n
    x = [(i + 0.5) * dx for i in range(n)]
    a1 = [a_bg] * n
    p = [P0 * (1 + amp * math.exp(-((xi - 0.5) / 0.03) ** 2)) for xi in x]
    m1 = [a_bg * RHO0[0]] * n
    m2 = [(1 - a_bg) * RHO0[1]] * n
    mom = [0.0] * n
    E = [a_bg * (pi + G[0] * PI[0]) / (G[0] - 1) + (1 - a_bg) * (pi + G[1] * PI[1]) / (G[1] - 1) for pi in p]
    c = wood_speed(a_bg) if use_K else isobaric_speed(a_bg)
    s = 1.05 * c                        # Rusanov 速度 = 该模型自身的声速
    t_end = travel / c
    steps = math.ceil(t_end / (0.4 * dx / s))
    dt = t_end / steps
    g = lambda q, i: q[min(max(i, 0), n - 1)]   # 零梯度虚拟单元
    for _ in range(steps):
        u = [mom[i] / (m1[i] + m2[i]) for i in range(n)]
        pr = [mix_pressure(a1[i], E[i] - 0.5 * mom[i] * u[i]) for i in range(n)]
        U = (m1, m2, mom, E)
        F = ([m1[i] * u[i] for i in range(n)], [m2[i] * u[i] for i in range(n)],
             [mom[i] * u[i] + pr[i] for i in range(n)], [(E[i] + pr[i]) * u[i] for i in range(n)])
        new = []
        for q, f in zip(U, F):
            fl = [0.5 * (g(f, i - 1) + g(f, i)) - 0.5 * s * (g(q, i) - g(q, i - 1)) for i in range(n + 1)]
            new.append([q[i] - dt / dx * (fl[i + 1] - fl[i]) for i in range(n)])
        na = []
        for i in range(n):
            K = kapila_K(a1[i], pr[i]) if use_K else 0.0
            dudx = (g(u, i + 1) - g(u, i - 1)) / (2 * dx)
            dadx = (g(a1, i + 1) - g(a1, i - 1)) / (2 * dx)
            lap = g(a1, i + 1) - 2 * a1[i] + g(a1, i - 1)
            na.append(a1[i] + dt * (-u[i] * dadx + K * dudx + s * lap / (2 * dx)))
        m1, m2, mom, E = new
        a1 = na
    pr = [mix_pressure(a1[i], E[i] - 0.5 * mom[i] ** 2 / (m1[i] + m2[i])) for i in range(n)]
    ipk = max((i for i in range(n) if x[i] > 0.6), key=lambda i: pr[i])
    return c, (x[ipk] - 0.5) / t_end, steps
 
print("alpha_air  wood[m/s]  K=0[m/s]")
for a in [0.0, 1e-4, 1e-3, 0.01, 0.1, 0.5, 0.9, 0.99, 1.0]:
    print(f"{a:9.4f}  {wood_speed(a):9.1f}  {isobaric_speed(a):8.1f}")
amin = min((k / 1000 for k in range(1, 1000)), key=wood_speed)
print(f"Wood minimum: {wood_speed(amin):.2f} m/s at alpha_air = {amin:.3f}")
print(f"K at alpha_air=0.5: {kapila_K(0.5, P0):.6f}")
print("\n1D pulse, alpha_air = 0.5, N = 400")
for useK in (True, False):
    c, meas, st = pulse_speed(0.5, useK)
    tag = "with K" if useK else "K = 0 "
    print(f"{tag}  theory {c:7.2f}  measured {meas:7.2f}  ratio {meas / c:.3f}  steps {st}")

alpha_air  wood[m/s]  K=0[m/s]
   0.0000     1624.9    1624.9
   0.0001      956.6    1624.3
   0.0010      364.8    1618.9
   0.0100      118.6    1567.2
   0.1000       39.4    1228.3
   0.5000       23.6     745.3
   0.9000       39.2     581.7
   0.9900      112.4     540.9
   1.0000      341.6     341.6
Wood minimum: 23.65 m/s at alpha_air = 0.501
K at alpha_air=0.5: 0.499947
 
1D pulse, alpha_air = 0.5, N = 400
with K  theory   23.65  measured   23.55  ratio 0.996  steps 316
K = 0   theory  745.30  measured  742.19  ratio 0.996  steps 315

在浏览器里直接运行上面的代码。同样的Rusanov计算在两根管中逐步推进。

默认参数(N = 400, α = 0.5)下按Play,表中measured最终停在23.55和742.19。把N降到100,两个比值一起下降,但两模型之间31倍的差距不变。移动α,就能亲自测出输出表里没有的混合比下的脉冲速度。

23.55与742.19说明了什么#

两次计算使用相同的网格、相同的通量、相同的初始脉冲。唯一的区别是$K$这一行。但脉冲速度却相差31.5倍:23.55 m/s对742.19 m/s。两者都与各自模型的理论声速相差不到0.4%。剩下的0.4%是一阶Rusanov扩散把峰值往后拖的结果。

也就是说,数值格式没有错。两个模型求解的是不同的物理。$K=0$模型准确地传播了这个脉冲,只是速度不对。无论怎样加密网格,742都不会降到23。删掉$K$后界面附近的水被过度压缩,是同一项的另一张脸。今天的数字,是这一项在均匀混合区留下的痕迹。

把表中两列在 $\alpha = 0.5$ 相除,$745.3/23.6 \approx 31.6$。脉冲测量是 $742.19/23.55 \approx 31.5$,几乎一样。每次测量都是理论值的 0.996 倍。23.65 减去 23.55 是 0.10 m/s,745.30 减去 742.19 是 3.11 m/s。绝对误差在快的一侧更大,比值却相同。一阶 Rusanov 把波峰按同一比例往后拖。步数是 316 和 315。每次计算用自己的声速取 CFL,所以慢声和快声走同样的步数。变的是花费的时间,比值和声速比一致。

$K$ 也印在输出里。$\alpha = 0.5$ 时 $K = 0.499947$。两相体积模量差到无穷,这个式子里的 $K$ 就是 0.5。印出的值比 0.5 小 $5.3\times10^{-5}$。这和水比空气硬大约 1.9 万倍的说法相符。$\alpha = 0.01$ 时 Wood 声速是 118.6 m/s,约为纯水 1624.9 m/s 的 1/14。$\alpha = 0.0001$ 时已经是 956.6 m/s,从 1624.9 掉了 668 m/s。同一位置的 $K=0$ 列是 1624.3 m/s。气泡把声速削下去,只发生在有 $K$ 的时候。

最小值 23.65 m/s 在 $\alpha = 0.501$。和 0.5 处的 23.6 m/s 实质上相同。体积模量比大,最小值就贴在 0.5。$\alpha = 0.99$ 时 $K=0$ 是 540.9 m/s,比纯空气 341.6 m/s 快 199 m/s。那是 1% 的水把 $\pi$ 漏到了空气一侧。同一 $\alpha$ 的 Wood 值是 112.4 m/s,反而比纯空气慢。

从含气1%到99%,两条曲线分开的地方#

左列对一点点空气都很敏感。只要 0.01%,1625 m/s 就变成 957 m/s,1% 时是 119 m/s。最小值贴在 $\alpha\approx0.5$,是因为体积模量差得大。

右列在含气 99% 时比纯空气还快。那是 1% 的水把 $\pi$ 漏到空气一侧,报告里 "5 Eq. Trans" 偏离约化模型的也是这里。水滴把声波几乎原路送回,是声阻抗和 sharp 界面里的另一件事。今天的两条曲线问的是界面反射之前,均匀混合物自己有怎样的声速。

报告与实验对照了什么

报告并没有停在声速曲线上。它在空气–水激波管中并排绘出两个模型的相密度和内能。接着绘出环氧树脂–尖晶石粉末混合物的冲击Hugoniot,即激波速度$U_s$与粒子速度$U_p$的关系。图中同时给出五方程约化模型、七方程模型和实验点。

这个对比之所以重要,有其原因。$K\nabla\cdot u$是非守恒项,因此仅凭模型无法唯一确定激波跳跃条件。数值黏性的形式可能改变激波后的状态。所以需要通过图逐一确认约化模型与七方程模型及实验有多接近。对于今天这样的弱波,差别会消失。对于强激波,则需要单独验证。

脉冲不在23 m/s附近，说明代码丢掉了K#

在均匀混合物 $\alpha=0.5$ 里放一个小压力脉冲,测它的速度。靠近 23 m/s,$K$ 还在。靠近 745 m/s,也就是表里 $K=0$ 的理论值,体积分数只是被运着走。CFL 里用的声速要一起看。按 Wood 声速传播的代码若塞进冻结声速,时间步会保守几十倍;$K=0$ 模型若用上 23.6 m/s,就会越过 CFL。这股刚性分给声学还是输运,和Peluchon 的声-输运 IMEX选的位置叠在一起。

含气只要 0.1%,Wood 声速就是 365 m/s,不到水 1625 m/s 的四分之一。空化刚开始,或者泵吸进气泡时,这个数决定压力波何时到达。那个格子若仍印出 1600 m/s,模型已经走到没有 $K$ 的一侧。0.01% 时也已经掉到 957 m/s,还用纯水声速,到达时间就是错的。