给每一相各自的方程 — Baer–Nunziato 七方程模型

按代码里的μ和时间步相乘,第一步的体积分数大约挪动0.00126。

1986年,Baer 和 Nunziato 想解开固体推进剂里的火焰怎样变成爆轰。火药是一堆细颗粒。热气体从缝里渗进去,把颗粒烧掉。同一个点上,固体和气体各有各的压力,各有各的速度。假定单一压力、单一速度的混合模型装不下这一刻。于是两个相各自拿到自己的守恒律。

把松弛率 μ 和 λ 放到 0,两条压力曲线碰不到一起;把值抬高,它们被拉向同一个点。

압력 완화율 μ: 18속도 완화율 λ: 14초기 압력차 Δp: 1.6초기 속도차 Δu: 0.8

종료 시 압력차 0.000 bar · 속도차 0.000 m/s

μ·λ를 0으로 내리면 두 상은 영영 평형에 도달하지 못한다 — 이것이 7-방정식의 “비평형” 상태다.

看压力曲线和速度曲线是不是分开贴合。只抬一边的滑块,另一条曲线是否还停在原地。

从火药颗粒中诞生的七个方程

扩散界面,也就是把两相边界看成有限厚度的看法,有分叉。最简单的四方程认为压力、速度、温度都相等。Baer–Nunziato 站在另一头。不假设任何平衡。每个相保留自己的密度、速度和压力。再加一行输运方程,运送一个相所占的比例,也就是体积分数 $\alpha$ 。

这份一般性守住了双曲性。允许非平衡,模型就是无条件双曲的。把压力绑成一个的六方程,在某些条件下会失去双曲性,甚至给出声速平方为负的解。BN 避开那个陷阱。代价是方程变成七个,还夹进写不成标准守恒形式的项。

把七个方程铺开

两个相写成 $k=1,2$ 。先是体积分数。

\frac{\partial \alpha_1}{\partial t} + u_I \frac{\partial \alpha_1}{\partial x} = \mu\,(p_1 - p_2)

左边是体积分数被界面速度 $u_I$ 带走。右边是压力松弛。 $\mu$ 越大,平衡越快。其余六个是各相的质量、动量和能量。

\frac{\partial (\alpha_k \rho_k)}{\partial t} + \frac{\partial (\alpha_k \rho_k u_k)}{\partial x} = 0

\frac{\partial (\alpha_k \rho_k u_k)}{\partial t} + \frac{\partial (\alpha_k \rho_k u_k^2 + \alpha_k p_k)}{\partial x} = p_I \frac{\partial \alpha_k}{\partial x} + \lambda\,(u_{k^*} - u_k)

\frac{\partial (\alpha_k E_k)}{\partial t} + \frac{\partial \big(\alpha_k (E_k + p_k)\,u_k\big)}{\partial x} = p_I u_I \frac{\partial \alpha_k}{\partial x} + \lambda u_I (u_{k^*} - u_k)

$\rho_k$ 是相密度, $u_k$ 是相速度, $p_k$ 是相压力。 $E_k=\rho_k e_k+\tfrac12\rho_k u_k^2$ 是相的总能密度。 $k^*$ 是另一个相。 $\lambda$ 是速度松弛率。每个相用自己的状态方程封闭。常用的是刚性气体。

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

$\gamma_k$ 是比热比。 $p_{\infty,k}$ 是表示分子之间排斥的压力常数。像水这样几乎压不动的液体, $p_\infty$ 有好几个 GPa。

非守恒项与界面变量 — 真正困难的地方

动量和能量的右边有 $\partial\alpha_k/\partial x$ 。这一项捆不进散度。写不成 $\partial_t U+\partial_x F(U)=0$ 。Godunov 一类的有限体积法在这种非守恒项面前会晃。 $\alpha$ 在间断上跳跃时,该乘哪个值并不清楚。离散错了,压力和速度本来均匀,界面上也会蹦出虚假振荡。

试金石是 Abgrall 准则。压力和速度一开始均匀,无论体积分数怎么铺,这份均匀都得保住。两相以同一个压力、同一个速度流过,界面只是被运走,不该造出新波。本来均匀、压力却跳了21%的记录就是非守恒项乘错时,这条准则立刻破裂的例子。

界面量 $u_I$ 、 $p_I$ 也是模型选定的值。常见的 Saurel–Abgrall 选择是质量加权速度和体积加权压力。

u_I = \frac{\alpha_1 \rho_1 u_1 + \alpha_2 \rho_2 u_2}{\alpha_1 \rho_1 + \alpha_2 \rho_2}, \qquad p_I = \alpha_1 p_1 + \alpha_2 p_2

这个选择让模型对称。每个单元都能用同一套方程、同一套数值方法。

松弛:走向平衡之路

右边的 $\mu(p_1-p_2)$ 和 $\lambda(u_{k^*}-u_k)$ 是松弛算子。把 $\mu$ 和 $\lambda$ 加大,两相几乎立刻被拉到同一个压力、同一个速度。把松弛率送到无穷大,得到的就是更简单的平衡模型。

前面的画面里把松弛率降到 0,两条曲线碰不到。这是七方程允许的非平衡。 $\mu$ 越大,压力曲线越陡地收到一个点上。压力和速度是分开松弛的。

层级:7 → 6 → 5 → 4#

把松弛当成瞬时的,也就是强制平衡,就有一个方程变成代数约束,然后消失。绑住速度,是六方程。再绑住压力,是 Kapila 的五方程。再绑住温度,是四方程。BN 站在这层阶梯的顶端。

五方程在界面捕捉里用得最多。混合声速对体积分数是非单调的 Wood 曲线,所以体积分数源项很硬。BN 把这份硬度拆成松弛源。数值上更灵活。

Python — 积分压力松弛算子#

只留下松弛的骨架。两相的部分质量和内能保持不变。只有体积分数 $\alpha$ 按 $d\alpha/dt=\mu(p_1-p_2)$ 移动。积分刚性气体压力怎样随 $\alpha$ 到达平衡。

import numpy as np
 
# 刚性气体(stiffened gas)两相的“压力松弛”算子。
# 每一相的部分质量 m_k 与内能 E_k 守恒,
# 只有体积分数 alpha 按松弛 ODE  da/dt = mu (p1 - p2) 变化。
 
def phase_pressures(alpha, mass, energy, gamma, p_inf):
    frac = np.array([alpha, 1.0 - alpha])
    rho  = mass / frac          # 相密度 = 部分质量 / 体积分数
    e    = energy / mass        # 比内能
    return (gamma - 1.0) * rho * e - gamma * p_inf
 
def relax_pressure(alpha0, mass, energy, gamma, p_inf,
                   mu=4.0e-8, dt=1.0e-3, nmax=200000, tol=1.0):
    alpha, hist = alpha0, []
    for n in range(nmax):
        p = phase_pressures(alpha, mass, energy, gamma, p_inf)
        hist.append((n * dt, alpha, p[0], p[1]))
        gap = p[0] - p[1]
        if abs(gap) < tol:                       # 压力差消失即终止
            break
        alpha += dt * mu * gap                    # da/dt = mu (p1 - p2)
        alpha = min(max(alpha, 1e-4), 1.0 - 1e-4)
    return alpha, np.array(hist)
 
# 类似水的相(1, 刚性气体) + 类似空气的相(2)
gamma  = np.array([4.4, 1.4])
p_inf  = np.array([6.0e8, 0.0])      # Pa
mass   = np.array([480.0, 0.55])     # 部分质量  (kg/m^3)
energy = np.array([3.93e8, 1.10e6])  # 部分内能 (J/m^3)
 
alpha_eq, hist = relax_pressure(0.50, mass, energy, gamma, p_inf)
p0  = phase_pressures(0.50, mass, energy, gamma, p_inf)
peq = phase_pressures(alpha_eq, mass, energy, gamma, p_inf)
print(f"初始: alpha=0.500  p1={p0[0]/1e5:8.3f} bar  p2={p0[1]/1e5:7.3f} bar")
print(f"平衡: alpha={alpha_eq:.3f}  p1={peq[0]/1e5:7.3f} bar  p2={peq[1]/1e5:7.3f} bar")
print(f"{len(hist)} 步,  最终压力差 {abs(peq[0]-peq[1]):.3e} Pa")

运行结果是这样。

初始: alpha=0.500  p1= 324.000 bar  p2=  8.800 bar
平衡: alpha=0.506  p1=  8.906 bar  p2=  8.906 bar
75 步,  最终压力差 9.012e-01 Pa

→

Baer–Nunziato (full)

없음 — 압력·속도·온도 모두 독립

미지수: α, (ρ, u, p)₁, (ρ, u, p)₂

장점
가장 일반적. 무조건 쌍곡형(hyperbolic). 상마다 독립 EOS.

대가
비보존 항 + 두 속도장. 계면 변수 PI·uI 모델링이 까다롭다.

평형을 하나씩 더 부과할수록 방정식은 줄지만 — 일반성도 함께 줄어든다.

按下某个模型,被强制的平衡是哪一种、声速曲线在哪里折成非单调,都会换掉。刚降到五方程的那一刻,体积分数源变硬的位置就在这个画面上。

印出的 324 bar 是两项大数的差。 $\alpha=0.5$ 时液体密度是 480 除以 0.5,等于 960。比内能是 $3.93\times10^{8}/480=8.1875\times10^{5}$ 。 $(\gamma-1)\rho e$ 是 $2.6724\times10^{9}$ Pa。 $\gamma p_\infty$ 是 $4.4\times6\times10^{8}=2.64\times10^{9}$ Pa。相减得到 $3.24\times10^{7}$ Pa,也就是 324 bar。每一项只剩下 1.2%。

斜率也从同一个式子里出来。 $\alpha=0.5$ 时液体压力对 $\alpha$ 的导数约为 $-5.34\times10^{9}$ Pa。体积分数挪 0.001,液体压力大约是 53 bar。印出来的移动是从 0.500 到 0.506,差 0.006。气体压力是 $4.4\times10^{5}/(1-\alpha)$ ,0.5 时 8.8 bar,0.506 附近约 8.9 bar。平衡是液体降到气体已经站着的地方。

第一步的大小是代码里写明的系数。 $\mu=4\times10^{-8}$ ,时间步是 $10^{-3}$ ,初始差是 324 bar 减去 8.8 bar,即 $3.152\times10^{7}$ Pa。体积分数增量是这三项相乘,约为 $1.26\times10^{-3}$ 。这大约是印出的总移动 0.006 的五分之一。75 步是此后差值变小、增量被切细的结果。1 Pa 的容差要到最后才有意义。

体积分数只从 0.500 挪到 0.506,动了 0.6%,液体压力却从 324 bar 塌到 8.9 bar。用刚性气体写的水几乎压不动,微小的体积变化会叫出巨大的压力变化。松弛平衡坐在气体压力附近,原因相同。硬的液体几乎不用膨胀,也能把压力对齐。

第一步松弛就把体积分数挪了0.001#

第一步已经把体积分数挪了 $1.26\times10^{-3}$ 。液体压力光凭这一下就掉几十 bar。非守恒项里的 $p_I$ 把这次崩溃乘到界面梯度上。松弛和通量放在同一步,假波的大小就在这第一步定下来。

如果把平流和压力在通量里彻底拆开,这个乘积去了哪里。把 Baer–Nunziato 通量拆成平流和压力的 TV 分裂还留着一个问题:这种拆分盖不盖得住非守恒项。324 bar 在一步里塌掉的时候,拆开的压力通量会让界面保持安静吗。