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

出口是开着的。波却折了回来。把可压缩程序的出口边界用外插草草处理，再去算火焰或湍流，域的正中央就会长出一段找不到来由的振荡。量一量它的周期，多半正好等于域长除以声速。计算域本身已经变成了一根共鸣管。在制造这段共鸣的边界上，把能算出来的值和必须编造的值分开。拨动那个编造值的系数，代价写在同一张表里。

把 Mach 滑块从正的亚声速推过声速，进来的箭头会出去。符号改成负，同一个面的条件数就跳。

亚声速出口的条件是一个。过了声速变成 0，翻成入口变成 4。箭头指向域内的时候，那条波才是编出来的数。

计算域的尽头不是物理

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

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

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

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

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

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

$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$ 移动。

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

特征线怎样在整个区域里把声波切开，写在超声速流动为什么不知道上游。这里只数这些线在出口一个面上被切断的条数。

把 $M$ 从 $+0.3$ 提到 $+1.4$，前面掉头的箭头就是 $\mathcal{L}_1$。条件数从 1 到 0。符号改成负，面变成入口，条件是 4。亚声速出口跑得好好的程序一到超声速就发散，多半是多挂了一个超出这个数目的条件。

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

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

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

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

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

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

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

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

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

封闭管上的 sigma 和 source q，放在输出表的正后方。

夹在两种失败之间的窄窗

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

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

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

频率越低反射越强。也就是说 $\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

fire pulse 打出的脉冲碰到出口时，红色曲线 $A^-$ 隆起就是反射。把 source q 调高且 $\sigma = 0$，白线会停在 $p_\infty$ 上方。反射没有了，压力却没抓住。

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

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

反射那一列是 0.0、1.9、7.3、24.5、46.8%。代码里 $C = 1$、$L = 1$、Mach 为 0，于是 $K = \sigma \left( 1 - M^2 \right) \frac{c}{L}$ 收成 $K = \sigma$。把 σ 从 0.25 加到 10，是 40 倍，反射从 1.9 到 46.8，大约 25 倍。不是线性。$|R| = \frac{K}{\sqrt{K^2 + 4\omega^2}}$ 会随 K 变大而饱和到 1。脉冲不是单一频率，而是宽度 0.09 的高斯。不要用一个 ω 把这些百分数倒推回来。方向对就够了。K 变大，反射增加，然后变钝。

压力那一列不跟这个单调走。$q = 0.3$ 的那次运行里，$\sigma = 0$ 的平均是 +0.3000，和热源一样大。$\sigma = 0.25$ 是 +0.0007，$\sigma = 1$ 是 +0.0006。相对 0.3000，0.0007 是 0.23%。到这里，松弛几乎把偏移擦掉了。$\sigma = 4$ 时是 +0.0227，大约是 0.0007 的 32 倍。$\sigma = 10$ 时是 −0.1142，符号反过来。绝对值大约是 0.0007 的 160 倍。反射 46.8% 的边界抓不住平均，只会晃它。表上偏移的底部在 0.25 和 1.0。Poinsot 和 Lele 用的 0.25 是这个底部的左端。0.58 是 B1 的对应系数，不会原样落在这一列上。

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

在边界上反射的并不只有物理声波。波长短于约四倍网格间距的成分根本不是 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 里之后没有任何东西能把它们清掉。

周期若是 2L/c，先怀疑出口的 σ#

$\sigma$ 不是调参旋钮。一端是无处锚定的压力，另一端是变成共鸣管的计算域。0.25 附近只抓住低频，其余放走。

没有来由的振荡，先数入射特征波的条数，看它和强加的条件数是否相同。周期是 $2L/c$ 的倍数，那就是边界，不是物理。把 $\sigma$ 调低，振荡变小，病根在边界；平均压力开始漂，就已经跨到另一侧。从波的条数数起、再去动出口，写在把 NSCBC 做成无反射出口的上一篇整理里。