出口明明开着，波却折了回来 — NSCBC 与 σ 决定的反射

出口是开着的。波却折了回来。把可压缩程序的出口边界用外插草草处理，再去算火焰或湍流，域的正中央就会长出一段找不到来由的振荡。量一量它的周期，多半正好等于域长除以声速。计算域本身已经变成了一根共鸣管。今天要看的是：在制造这段共鸣的边界上，什么是能算出来的，什么是必须编造的，以及调节这个编造过程的一个系数，究竟索取什么代价。

计算域的尽头不是物理

物理域没有尽头。燃烧室出口之外，空间仍在继续。可网格必须在某处停下，最后一个单元没有邻居。没有邻居就求不出导数，求不出导数就推不动控制方程。

在 RANS 程序里，这个问题隐藏了很久。湍流粘性和人工粘性都很大，边界上错造出来的波走几个单元就消散了。到了 LES 和 DNS，账就不一样了。人工粘性接近零，湍流粘性也压到最低。边界造出的误差不会消失，它横穿整个域再折返回来。

Poinsot 和 Lele 在 1992 年整理的处方把思路翻了过来。不再在边界上外插变量，而是去数穿过边界的波，并逐个确定其振幅。它把 Euler 方程的特征边界条件（ECBC）推广到含粘性项的 Navier–Stokes，因此称为 NSCBC（Navier–Stokes Characteristic Boundary Conditions）。整个过程一行外插也不用。

边界上能数清的，和必须编造的

设边界位于 $x_1 = L$ 。把 $x_1$ 方向的项重新按波的形式归并，连续方程写成：

\frac{\partial \rho}{\partial t} + d_1 + \frac{\partial (\rho u_2)}{\partial x_2} + \frac{\partial (\rho u_3)}{\partial x_3} = 0

其中 $\rho$ 是密度， $u_i$ 是速度分量， $d_1$ 汇集了垂直于边界方向的贡献。这个 $d$ 向量正是特征分析的产物，波的振幅 $\mathcal{L}_i$ 就藏在里面。

\mathcal{L}_1 = (u_1 - c)\left( \frac{\partial p}{\partial x_1} - \rho c \frac{\partial u_1}{\partial x_1} \right)

\mathcal{L}_5 = (u_1 + c)\left( \frac{\partial p}{\partial x_1} + \rho c \frac{\partial u_1}{\partial x_1} \right)

$c$ 是当地声速（ $c^2 = \gamma p / \rho$ ）， $p$ 是压力。 $\mathcal{L}_1$ 是沿 $x_1$ 负方向传播的声波振幅变化率， $\mathcal{L}_5$ 则是沿正方向传播的。剩下三个随流体一起被输运： $\mathcal{L}_2$ 携带熵， $\mathcal{L}_3$ 和 $\mathcal{L}_4$ 携带切向速度 $u_2$ 、 $u_3$ ，三者都以速度 $u_1$ 移动。

关键在于速度的符号。波若离开计算域，其振幅可以由内部点算出来。波若进入计算域，这个信息根本不在解里面，只能编造。进入的波有几条，这个边界上允许给出的物理边界条件就有几个。

在下面的示意图里直接拖动 Mach 数。五条特征线穿过边界的方向会实时改变。

Mach M+0.30

boundary

把 $M$ 从 $+0.3$ 提到 $+1.4$ ，原本进入的 $\mathcal{L}_1$ 掉头，所需条件数从 1 降到 0。把符号翻成负数，同一个面变成入口，条件数跳到 4。如果一段在亚声速出口跑得好好的程序一到超声速就发散，通常是多强加了一个这张表不允许的条件。

LODI — 编造入射振幅的规则#

编造入射振幅需要依据。NSCBC 在边界的每一点上建立一个局部一维无粘系统，把切向项、粘性项和反应项全部抹掉。这就是 LODI（Local One Dimensional Inviscid）关系。

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

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

LODI 关系不是物理条件，也不是真正拿来求解的方程。它只用于估计入射的 $\mathcal{L}_i$ 。步骤分三步：把被物理条件约束的守恒方程从方程组中删去；用被删方程对应的 LODI 关系，把未知的 $\mathcal{L}_i$ 表示成已知的 $\mathcal{L}_i$ ；再用剩下的方程推进其余变量。

完全无反射出口在这里做了最简单的选择——宣布外面没有声波进来。

\mathcal{L}_1 = 0

入射为零，反射也为零。看上去很干净。可是这个式子里，哪儿都找不到外界压力 $p_\infty$ 。

σ = 0 的代价：压力回不到 p∞#

没有 $p_\infty$ ，意味着边界不知道"现在的压力该是多少"。域内放热把压力顶上去，没有任何地方提供把这个偏移拉回来的恢复力。问题不再适定。

Rudy 与 Strikwerda 的处方是：不把入射波置零，而是让它正比于压力差。

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

$L$ 是域的特征长度， $M$ 是最大 Mach 数， $\sigma$ 是整个处方里唯一的自由参数。 $\sigma = 0$ 就退回到完全无反射；把 $\sigma$ 调大，边界开始把压力往 $p_\infty$ 拽。

在下面的模拟里亲手操作一下。左端是封闭的管，右端是出口。按 fire pulse 打出一道压力波，再移动 sigma 滑块。

sigma0.00

source q0.30

fire a pulse and watch the red curve at the outlet — that is the reflection. then drop sigma to 0 and watch the white curve settle above p_inf instead of returning to it.

要看两件事。第一，脉冲抵达出口的瞬间，红色曲线（入射波 $A^-$ ）会不会隆起——那就是反射。第二，把 source q 调高后看下方的压力历史：在 $\sigma = 0$ 时白色曲线不会回到 $p_\infty$ 那条线，而是停在上方。反射消掉了，压力却丢了。

夹在两种失败之间的窄窗

$\sigma$ 处在两种方向相反的失败之间。

处方	出口上做的事	失败方式
B1（外插 + Riemann 不变量）	外插速度和密度，只松弛压力	外插造出的伪波
B2（NSCBC， $\sigma = 0$ ）	$\mathcal{L}_1 = 0$	平均压力不被 $p_\infty$ 锚定
B3（NSCBC， $\sigma > 0$ ）	$\mathcal{L}_1 = K(p - p_\infty)$	$\sigma$ 一大就反射
B4（反射出口）	压力固定， $\mathcal{L}_1 = \mathcal{L}_5$	完全反射——域变成共鸣管

$\sigma$ 小则平均压力漂走， $\sigma$ 大则边界变硬、把声能扔回域内。Poinsot 和 Lele 实际采用的值是 $\sigma \approx 0.25$ ，对应外插型 B1 的系数则取 $\sigma = 0.58$ 。这些数字都不是从理论推出来的，而是在两种失败之间挑出的折中。

看频率依赖性就明白为什么是折中。对角频率为 $\omega$ 的声波，这个边界的反射系数为

|R| = \frac{K}{\sqrt{K^2 + 4\omega^2}}

频率越低反射越强。也就是说 $\sigma$ 是一个"抓住低频、放走高频"的滤波器。平均压力属于 $\omega \to 0$ 的分量，会被抓住；我们想放出去的声波则通过。这种分离能奏效的区间，就是那扇窄窗。

用代码量一量反射率和压力偏移

对一维线性声学，状态恰好分裂成两个特征振幅： $A^\pm = p' \pm \rho c u'$ ，分别以 $\pm c$ 移动。取 $\Delta t = \Delta x / c$ ，每步正好平移一个单元，因此屏幕上的每一处抖动都来自边界条件，而非格式本身。

import numpy as np
 
N, C, L = 240, 1.0, 1.0
DX = L / N
DT = DX / C                       # 正好平移一个单元 — 无格式扩散
 
 
def outlet_relax_k(sigma, mach=0.0):
    """NSCBC 松弛系数 K = sigma (1 - M^2) c / L"""
    return sigma * (1.0 - mach ** 2) * C / L
 
 
def duct_step(ap, am, am_b, k_relax, q):
    """A+ 右移一格，A- 左移一格。出口处的 A- 必须编造。"""
    ap[1:] = ap[:-1].copy()
    ap[0] = 0.0
    am[:-1] = am[1:].copy()
 
    p_b = 0.5 * (ap[-1] + am[-1])
    am_b -= DT * k_relax * p_b     # L1 = K (p - p_inf)
    am[-1] = am_b
 
    ap[0] = am[0]                  # 左端封闭 (u = 0)
    ap += q * DT                   # 微弱的均匀放热
    am += q * DT
    return am_b
 
 
def measure_outlet(sigma, q, steps, pulse):
    ap, am, am_b = np.zeros(N), np.zeros(N), 0.0
    x = (np.arange(N) + 0.5) * DX
    if pulse:
        ap += np.exp(-((x - 0.30) / 0.09) ** 2)
    k = outlet_relax_k(sigma)
    refl = 0.0
    for n in range(steps):
        am_b = duct_step(ap, am, am_b, k, q)
        if pulse and n > 0.85 * N:
            refl = max(refl, np.abs(am[:-8]).max())
    return refl, float(np.mean(0.5 * (ap + am)))
 
 
for sigma in (0.0, 0.25, 1.0, 4.0, 10.0):
    r, _ = measure_outlet(sigma, q=0.0, steps=650, pulse=True)
    _, p = measure_outlet(sigma, q=0.3, steps=6000, pulse=False)
    print(f"sigma={sigma:5.2f}   反射={r * 100:5.1f}%   平均压力-p_inf={p:+.4f}")

运行结果：

sigma= 0.00   反射=  0.0%   平均压力-p_inf=+0.3000
sigma= 0.25   反射=  1.9%   平均压力-p_inf=+0.0007
sigma= 1.00   反射=  7.3%   平均压力-p_inf=+0.0006
sigma= 4.00   反射= 24.5%   平均压力-p_inf=+0.0227
sigma=10.00   反射= 46.8%   平均压力-p_inf=-0.1142

$\sigma = 0$ 把反射彻底清零，却把放热造出的 0.3 偏移原封不动留在那里。 $\sigma = 0.25$ 时反射不到 2%，压力被压在距 $p_\infty$ 0.0007 以内。到这里都在预料之中。

意外在最后两行。把 $\sigma$ 提到 4 和 10，反射跳到 25% 和 47% 并不奇怪，但压力偏移也重新变坏了。边界一旦太硬，就会自己产生振荡，而振荡会摇动平均值。调大 $\sigma$ 并不是"用反射换压力控制"的交易。过了某个点，两样都会输掉。

网格造出的波会朝反方向折回

在边界上反射的并不只有物理声波。波长短于约四倍网格间距的成分根本不是 Navier–Stokes 的解，而是离散化的产物。Poinsot 和 Lele 把它们称为"q 波"，与物理的"p 波"区分开。

q 波的标志是群速度。即便是速度为 $V$ 的一维对流方程，短波长的群速度 $u_g$ 与 $V$ 符号相反。流动往右走，数值误差却逆流向左爬。更麻烦的是 $|u_g / V|$ 随空间格式的阶数升高而增大，所以高阶格式暴露得更多，而不是更少。

因此评价一种边界处理必须用两个反射系数：物理波的 $A_p / A_1$ 和数值波的 $A_q / A_1$ 。可用的处理在任何情况下都要满足 $A_q / A_1 \ll 1$ ；若号称无反射，还要做到 $A_p / A_1 \ll 1$ 。仅仅是用带陡峭梯度的初场起算，就会生成 q 波，而在 DNS 里之后没有任何东西能把它们清掉。

选出口条件时的一句话

$\sigma$ 不是调参旋钮，而是两种失败之间的一个坐标。一端是无处锚定的压力，另一端是变成共鸣管的计算域。0.25 附近之所以反复被推荐，是因为它恰好抓住低频、放走其余。

出口出现原因不明的振荡时，按这个顺序排查。先数一数该面上入射特征波的条数，确认此刻强加的条件数与之相符。再看振荡周期是不是 $2L/c$ 的倍数——若是，那就是边界反射，不是物理。最后把 $\sigma$ 调低。振荡变小，说明病根在边界；平均压力开始漂走，说明已经跨到了另一侧的失败里。