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

非平衡七方程模型的结构、非守恒项,以及松弛生成的模型层级

1986 年,Mark Baer 与 John Nunziato 想弄清固体推进剂中的火焰如何转变为爆轰。火药是无数细小颗粒的堆积。高温气体从缝隙渗入,点燃颗粒。在同一个点上,固体与气体以不同的压力、不同的速度共存。当时假定单一压力、单一速度的混合模型,无法捕捉这一瞬间。于是两人建立了一个给每一相各自完整守恒律的模型。本文沿着这个七方程模型的结构展开:为什么需要七个方程,非守恒项为何棘手,以及松弛(relaxation)如何孕育出更简单的模型。最后我们直接积分压力松弛算子,重现两相走向平衡的过程。

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

扩散界面(diffuse interface,把两相边界看作有限厚度的过渡区)方法有不同分支。最简单的四方程模型假定两相的压力、速度、温度全部相等。Baer–Nunziato(BN)站在另一个极端。它不假定任何平衡。 每一相都拥有自己的密度、速度和压力。在此之上还多出一个体积分数(volume fraction) $\alpha$ 的输运方程, $\alpha$ 表示某一相所占体积的比例。

为什么这种一般性重要?允许非平衡使模型无条件双曲(hyperbolic)。强制单一压力的六方程模型,在某些条件下会丧失双曲性,产生诸如声速平方为负之类的非物理解。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$ 为速度松弛率。每一相由各自的状态方程封闭。无论是推进剂还是空化,常用的是刚性气体(stiffened gas)。

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 型有限体积法在这类非守恒(non-conservative)项面前会失稳,因为当 $\alpha$ 跨越间断跳变时,该乘以哪个值并不明确。离散化一旦出错,即使在均匀的压力、速度场中,也会在界面两侧迸出虚假振荡。

筛除这一陷阱的试金石是 Abgrall 判据。若压力与速度起初均匀,则无论体积分数如何分布,这种均匀性都必须保持。 当两相以相同压力、速度流动时,界面只应被输运,而不应产生新的波。非守恒项的离散化正是为精确满足这一条件而设计的。

界面变量 $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$ 调大,两相几乎瞬间被拉到相同的压力与速度。这正是 BN 的妙处:松弛率趋于无穷的极限,恰好是一个更简单的平衡模型。

在下面的模拟中亲手操作一下。两相从不同的压力与速度出发,松弛率 $\mu$ 、 $\lambda$ 把它们拉向平衡。

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

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

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

把松弛率降到 0,两条曲线便永不相交。这就是七方程模型所允许的非平衡状态。 $\mu$ 越大,压力曲线收敛到一点就越陡。也请注意,压力与速度彼此独立地松弛。

层级:7 → 6 → 5 → 4#

假定松弛为瞬时(instantaneous),即强制平衡,则有一个方程退化为代数约束而消失。把速度锁定到平衡就得到六方程,再加上压力就是 Kapila 的五方程,再加上温度就是四方程。BN 位于这一层级的顶端。

下面点击各模型,比较它们各自假定了哪些平衡、又失去了什么。

→

Baer–Nunziato (full)

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

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

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

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

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

五方程(Kapila)在界面捕捉中应用最广,但其混合声速关于体积分数是非单调的(Wood 曲线),使得体积分数源项变得刚硬(stiff)。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

引人注目的是,体积分数仅从 0.500 移动到 0.506,变化 0.6%,液体压力却从 324 bar 崩落到 8.9 bar。以刚性气体建模的水近乎不可压,微小的体积变化便引发巨大的压力变化。松弛平衡落在气体压力附近也是同样的道理:刚硬的液体几乎不必膨胀就能匹配压力。

值得记住的#

七方程 = 非平衡。 给每一相各自的压力、速度、能量,换来的是无条件双曲性。
非守恒项是核心难点。 $p_I \partial_x \alpha$ 无法归并为守恒形式。Abgrall 判据(保持均匀压力与速度)是离散化的试金石。
松弛构筑层级。 $\mu, \lambda \to \infty$ 的极限给出 6 → 5 → 4 方程模型。BN 居于顶端,给出最一般的图景。