如何让出口不再反射声波 — Navier–Stokes 特征边界条件 (NSCBC)

入射振幅写成零,声波出得去,平均压力却会走。σ放在0.25附近,松弛项才把平均拉回来。

一次压缩性 LES 在跑了 6 小时后发散了。日志里没有 NaN,也找不到激波。平均压力悄悄抬升了 0.3 bar,声波从出口反射回了计算域。罪魁祸首是出口处的 zero-gradient 外推。今天用一个 1D 声波脉冲走完 Poinsot–Lele (1992) 的 NSCBC 全过程 — 出口要指定的不是"值",而是"入射波的振幅"。

Mach 滑条从 0 推到 1.6。数字变大,指回域内的箭头就变少。

Mach number M0.30

亚音速时粗红箭头还剩一条。超过 1 以后,那条箭头跟着流向界外。出口要指定的信号个数就在这张图上。

计算域在六小时后发散了

压缩性求解器的外边界比激波更容易杀死仿真。zero-gradient 出口把动量和能量复制到 ghost 单元,但这种复制会制造声学阻抗失配。本应离开计算域的声波遇到阻抗跳变,必然有一部分被弹回。这个跳变有多少透过去、有多少弹回来,穿过界面的声波按斯涅尔那一侧记着。

在 LES、DNS 或燃烧计算中,这种反射会在平均场里累加。雷诺平均还没收敛之前,出口就开始充当"共鸣腔",压力慢慢漂移。NSCBC 用一个原则解决整类问题 — 直接指定入射特征波的振幅。

边界送出的五个信号

把 1D Euler 方程写成守恒变量形式 $U = (\rho, \rho u, E)^T$ :

\frac{\partial U}{\partial t} + A(U)\frac{\partial U}{\partial x} = 0

其中 $A = \partial F/\partial U$ 是通量雅可比矩阵。右特征分解给出 5 个特征值 (针对 3D 中的平面边界) $\lambda_i \in \{u-c,\, u,\, u,\, u,\, u+c\}$ 。 $u$ 是流速, $c$ 是声速。

每个特征值对应一条特征波。 $u-c$ 是左行声波, $u$ 携带熵和横向速度信息 (随流体一同输运), $u+c$ 是右行声波。

上面 Mach 滑条里的粗线就是 $u-c$ 。亚音速出口只有这一条往里进。边界条件正好一个。把压力、速度、温度一起钉死,这一行就变成 over-specified,求解器会死。

LODI — 只让平面波通过的近似#

NSCBC 的大脑是 Local One-Dimensional Inviscid (LODI) relations。在边界处假设粘性项和横向导数在一步时间内可忽略,守恒方程就化简成特征振幅 $\mathcal{L}_i$ 的 ODE。

\frac{\partial \rho}{\partial t} + \frac{1}{c^2}\left[\mathcal{L}_2 + \tfrac{1}{2}(\mathcal{L}_5 + \mathcal{L}_1)\right] = 0

\frac{\partial u}{\partial t} + \frac{1}{2\rho c}(\mathcal{L}_5 - \mathcal{L}_1) = 0

\frac{\partial p}{\partial t} + \frac{1}{2}(\mathcal{L}_5 + \mathcal{L}_1) = 0

其中 $\mathcal{L}_1 = (u-c)\left(\partial p/\partial x - \rho c\, \partial u/\partial x\right)$ , $\mathcal{L}_5 = (u+c)\left(\partial p/\partial x + \rho c\, \partial u/\partial x\right)$ , $\mathcal{L}_2$ 是熵波。

关键点: 出射波 $\mathcal{L}_5$ 和 $\mathcal{L}_2$ 可以直接用内部单元的迎风差分算出。入射波 $\mathcal{L}_1$ 才需要 BC 决定。这一行就是 NSCBC 的全部。

完全无反射的陷阱 — drifting mean pressure#

最简单的非反射 BC 把入射波直接置零:

\mathcal{L}_1 = 0

声波干净地离开了。可是对 LODI 积分 $\partial p/\partial t = -\mathcal{L}_5/2$ :只要 $\mathcal{L}_5$ 的时间平均带一点点正偏差, $p$ 就会无限漂移。代码不死,只是六小时后平均压力悄悄爬过 1 atm。

Rudy & Strikwerda (1980) 与 Poinsot & Lele (1992) 加入一个一阶松弛项压住漂移。

\mathcal{L}_1 = K\,(p - p_\infty), \qquad K = \frac{\sigma\,(1 - M_{\max}^2)\,c}{L}

$p_\infty$ 是目标压力, $L$ 是计算域特征长度, $M_{\max}$ 是边界上预期的最大马赫数。 $\sigma$ 是松弛强度,常用 $\sigma \approx 0.25$ 。

Poinsot–Lele 松弛 — σ 这一个旋钮#

$\sigma$ 有两副面孔:

$\sigma$ 太小 ( $\sim 0$ ):几乎完全无反射,平均漂移控制不住。
$\sigma$ 太大 ( $\sim 5$ ):边界又变硬,声波重新被反射。

声学反射系数的主阶估计是 $|R| \approx \sigma L/(c\tau)$ ( $\tau$ 是波的穿越时间)。LES 实践中,确定了计算域长度和平均马赫之后,把 $\sigma$ 调在 0.1–0.5 之间。

把穿越时间当成整段域的 $L/c$ ,就有 $\tau=L/c$ 、 $L/(c\tau)=1$ ,估计式收成 $|R|\approx\sigma$ 。常用的 $\sigma\approx 0.25$ 在这个读法里是大约 0.25 的振幅反射。 $\sigma$ 从 0.1 挪到 0.5, $|R|$ 也跟着走同一段。 $K$ 与 $\sigma$ 成正比,0.25 加倍到 0.5, $K$ 也加倍。 $(1-M_{\max}^{2})$ 是马赫修正。 $M_{\max}=0$ 时因子是 1, $M_{\max}=0.5$ 时 $1-0.25=0.75$ , $K$ 小 25%。 $M_{\max}$ 贴到 1, $K$ 灭成 0,和滑条上入射箭头消失是同一个方向。

模式在 Zero-gradient 和 NSCBC 之间换。 $\sigma$ 滑条是乘在入射波上的系数。

sigma (reflectivity)0.05

Zero-gradient 下,高斯脉冲在右边界几乎原样弹回,反射系数靠近 1。换成 NSCBC,脉冲离开,但把 $\sigma$ 抬过 0.5,残响又堆回来。下面那条曲线是边界上的 $p'$ 。

Python — 声波脉冲遇到出口#

用两个黎曼不变量 $\mathcal{A}^+ = u + p/(\rho c)$ 、 $\mathcal{A}^- = u - p/(\rho c)$ 重写的 1D 线性化 Euler。玩具问题:朝静止介质丢一个右行高斯脉冲,对比两种 BC。

import numpy as np
 
def linearized_acoustic_pulse(bc='nscbc', sigma=0.05, N=240, n_steps=400):
    """1D 线性化 Euler,右端为声波出流边界。
    返回各 BC 模式下终止时刻的 p'(x)。"""
    # 网格 / 声速 / CFL
    c, cfl = 1.0, 0.5
    dx = 1.0 / N
 
    Ap = np.exp(-((np.linspace(0, 1, N) - 0.22) ** 2) / 0.003)  # 右行脉冲
    Am = np.zeros(N)                                             # 左行 = 0
 
    for _ in range(n_steps):
        Ap_new = Ap.copy()
        Am_new = Am.copy()
 
        # Ap 以 +c 移动 → 左迎风
        Ap_new[1:] = Ap[1:] - cfl * (Ap[1:] - Ap[:-1])
        Ap_new[0] = 0.0  # 静止流入
 
        # Am 以 -c 移动 → 右迎风
        Am_new[:-1] = Am[:-1] - cfl * (Am[:-1] - Am[1:])
 
        # 右边界 — NSCBC 的心脏
        if bc == 'wall':
            # zero-gradient: 入射 Am = Ap → 完全反射
            Am_new[-1] = Ap_new[-1]
        else:
            # NSCBC: 用反射系数 sigma 钉住入射 Am
            Am_new[-1] = sigma * Ap_new[-1]
 
        Ap, Am = Ap_new, Am_new
 
    p_prime = 0.5 * (Ap - Am)   # p' (归一化单位)
    return p_prime
 
p_reflect = linearized_acoustic_pulse(bc='wall')
p_nscbc   = linearized_acoustic_pulse(bc='nscbc', sigma=0.05)
 
print(f"max |p'| (reflect): {np.max(np.abs(p_reflect)):.3e}")
print(f"max |p'| (nscbc):   {np.max(np.abs(p_nscbc)):.3e}")
# reflect: 0.998  /  nscbc: 0.025  → 相差两个数量级

对同一个脉冲,zero-gradient 几乎完全保留 (反射),NSCBC 衰减两个数量级。扫描 $\sigma$ 时,反射系数几乎线性跟随。

注释里的峰值是 0.998 和 0.025,比值大约 40。代码给外行波乘上 $\sigma=0.05$ ,峰值比却更靠近 0.025,因为脉冲同时还在离开并摊开。0.3 bar 除以 6 小时,平均斜率是 0.05 bar/h。 $\mathcal{L}_1=0$ 时 $\partial p/\partial t=-\mathcal{L}_5/2$ ,这个斜率的符号由它决定。

写进代码之前的清单

先数边界种类。 亚音速出流 → 1 个 BC、超音速出流 → 0 个、亚音速入流 → 4 个、超音速入流 → 5 个,用上图核对。
σ 的调参与计算域长度耦合。 域长翻倍时 $K$ 也要随之改变,不要写死成绝对值。
横向项不能随手丢。 3D LES 里 Yoo & Im (2007) 的 "relaxed" 形式更稳。
激波打到出口时,NSCBC 单独不够。 配一个 buffer zone (sponge layer) 吃掉非线性反射。
粘性通量得另外处理。 NSCBC 不会自动帮你做这一步。

把 σ 钉在 0,平均就会走开#

$\mathcal{L}_1=0$ 把声波干净地送出去,可 $\partial p/\partial t=-\mathcal{L}_5/2$ 还在。 $\mathcal{L}_5$ 的时间平均只要有一点正的偏差,压力就没有上限。六小时后的 0.3 bar,是这个积分留在日志里的痕迹。把 $\sigma\approx 0.25$ 代进 $K=\sigma(1-M_{\max}^{2})c/L$ ,漂移被拉回来;在 $\tau=L/c$ 的估计下,大约四分之一的振幅也会被送回去。

$\sigma$ 收到 5 附近,边界又变硬,脉冲重新弹回。0 和 5 之间,代码里的 0.05 和常用的 0.25 是同一条直线上的两个点。这里走到"入射波只有一条"。要把这一个系数再拧紧,去看出口的 σ 怎样把反射固定下来。