The Outlet Is Open, Yet the Wave Comes Back — NSCBC and the Price of σ

The outlet is open. The wave comes back anyway. Handle the outflow boundary of a compressible code with a casual extrapolation, run a flame or a turbulent patch, and an oscillation with no physical cause grows in the middle of the domain. Measure its period and it usually matches the domain length divided by the speed of sound. Your computational box has quietly become a resonator. Split what that boundary can compute from what it has to invent. The coefficient that moves the invented part has a price, and the price is in the table.

Drag the Mach slider from a positive subsonic value through sonic. Incoming arrows leave. Flip the sign negative and the condition count on that same face jumps.

A subsonic outlet carries one condition. Past sonic the count is zero. Flip the face into an inlet and it becomes four. A wave is an invented number only while its arrow points into the box.

The Edge of the Box Is Not Physics#

A physical domain has no edge. Space continues past the combustor exit. A grid, though, has to stop somewhere, and the last cell has no neighbor. Without a neighbor there is no derivative, and without a derivative there is no governing equation to advance.

In RANS codes the problem stayed hidden for a long time. Turbulent and artificial viscosity are large, so a badly manufactured wave at the boundary dies within a few cells. LES and DNS change the accounting. Artificial viscosity is near zero and turbulent viscosity is at a minimum. Whatever error the boundary creates does not dissipate — it crosses the domain and comes back.

The recipe Poinsot and Lele assembled in 1992 inverts the approach. Instead of extrapolating variables at the boundary, you count the waves crossing it and fix each amplitude. It extends the Euler characteristic boundary conditions (ECBC) to Navier–Stokes with its viscous terms, hence NSCBC — Navier–Stokes Characteristic Boundary Conditions. Not a single line of extrapolation is used.

What a Boundary Can Count, and What It Must Invent#

Put the boundary at $x_1 = L$. Recast the $x_1$-direction terms as waves and the continuity equation reads:

Here $\rho$ is density, $u_i$ the velocity components, and $d_1$ collects the contributions normal to the boundary. That $d$ vector is the product of the characteristic analysis, and inside it live the wave amplitudes $\mathcal{L}_i$.

$c$ is the local speed of sound ($c^2 = \gamma p / \rho$) and $p$ the pressure. $\mathcal{L}_1$ is the amplitude variation of the acoustic wave running in the negative $x_1$ direction, $\mathcal{L}_5$ the one running positive. The remaining three ride along with the fluid: $\mathcal{L}_2$ carries entropy, $\mathcal{L}_3$ and $\mathcal{L}_4$ carry the transverse velocities $u_2$ and $u_3$, all three moving at $u_1$.

What splits the cases is the sign of those speeds. If a wave leaves the domain, its amplitude is computable from interior points. If it enters, that information is nowhere in the solution. You have to invent it. The number of entering waves is exactly the number of physical boundary conditions the problem allows.

How characteristics split sound across a whole domain is in why a supersonic flow does not know its upstream. Here the count is only the waves cut by one outlet face.

Raise $M$ from $+0.3$ to $+1.4$ and the arrow that flipped above is $\mathcal{L}_1$. The condition count goes from one to zero. A negative sign turns the face into an inlet and the count is four. A code that runs at subsonic outflow and diverges once the exit is supersonic is usually imposing one condition more than this ledger allows.

LODI — the Rule That Manufactures an Amplitude#

Inventing an incoming amplitude needs a justification. At each boundary point NSCBC builds a local one-dimensional inviscid system, dropping every transverse, viscous, and reaction term. Those are the LODI relations.

LODI relations are not physical conditions, and they are not the equations you actually solve. They exist only to estimate the incoming $\mathcal{L}_i$. The procedure is three steps: delete each conservation equation whose variable is physically imposed, use the matching LODI relation to write the unknown $\mathcal{L}_i$ in terms of the known ones, then advance the remaining equations for everything else.

A perfectly non-reflecting outlet makes the simplest possible choice here — it declares that no acoustic wave arrives from outside.

Zero incoming wave, zero reflection. It looks clean. But notice that the outer pressure $p_\infty$ appears nowhere in that statement.

What σ = 0 Costs: Pressure That Never Returns to p∞#

No $p_\infty$ means the boundary has no idea what the pressure is supposed to be. Heat released inside the domain raises the pressure, and nothing anywhere pulls that offset back. The problem stops being well posed.

The fix from Rudy and Strikwerda is to stop setting the incoming wave to zero and instead tie it to the pressure difference:

$L$ is a characteristic domain length, $M$ the maximum Mach number, and $\sigma$ is the only free parameter in the whole prescription. At $\sigma = 0$ you are back to perfectly non-reflecting. Raise $\sigma$ and the boundary starts dragging the pressure toward $p_\infty$.

The closed duct, with sigma and source q, sits immediately after the printed table.

The Narrow Window Between Two Failures#

$\sigma$ sits between two failures that pull in opposite directions.

Small $\sigma$ lets the mean pressure drift; large $\sigma$ stiffens the boundary until it throws acoustic energy back inside. The values Poinsot and Lele actually used were $\sigma \approx 0.25$, and $\sigma = 0.58$ for the equivalent coefficient in the extrapolation-based B1. Neither number comes from theory. Both were picked from between the two failures.

The frequency dependence shows why it is a compromise. For an acoustic wave of angular frequency $\omega$, this boundary reflects with

Lower frequencies reflect more. So $\sigma$ is really a filter that grabs low frequencies and lets high ones through. Mean pressure is the $\omega \to 0$ component, so it gets grabbed; the acoustics you want to evacuate pass. The window where that separation works is narrow.

Measuring Reflection and Offset in Code#

For one-dimensional linear acoustics the state splits exactly into two characteristic amplitudes: $A^\pm = p' \pm \rho c u'$, traveling at $\pm c$. Choose $\Delta t = \Delta x / c$ and each one shifts exactly one cell per step, so every wiggle on screen comes from the boundary condition rather than the scheme.

import numpy as np
 
N, C, L = 240, 1.0, 1.0
DX = L / N
DT = DX / C                       # exact one-cell shift: no scheme diffusion
 
 
def outlet_relax_k(sigma, mach=0.0):
    """NSCBC relaxation coefficient 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+ one cell right, A- one cell left. At the outlet A- must be invented."""
    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]                  # closed end at x = 0 (u = 0)
    ap += q * DT                   # weak uniform heat release
    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}   reflected={r * 100:5.1f}%   mean p - p_inf={p:+.4f}")

What it prints:

sigma= 0.00   reflected=  0.0%   mean p - p_inf=+0.3000
sigma= 0.25   reflected=  1.9%   mean p - p_inf=+0.0007
sigma= 1.00   reflected=  7.3%   mean p - p_inf=+0.0006
sigma= 4.00   reflected= 24.5%   mean p - p_inf=+0.0227
sigma=10.00   reflected= 46.8%   mean p - p_inf=-0.1142

When a fire pulse hits the outlet, a rising red curve $A^-$ is the reflection. With source q up and $\sigma = 0$, the white curve sits above $p_\infty$. The reflection is gone. The pressure is not.

At $\sigma = 0$ the reflection is gone completely, and the 0.3 offset built by the heat release sits there untouched. At $\sigma = 0.25$ reflection is under 2% and the pressure is held within 0.0007 of $p_\infty$. So far, as expected.

The surprise is the last two rows. Pushing $\sigma$ to 4 and 10 raises reflection to 25% and 47%, which is unsurprising — but the pressure offset gets worse again. A boundary that stiff starts ringing on its own, and the ringing moves the mean. Raising $\sigma$ is not a trade where you buy pressure control with reflection. Past a certain point you lose both.

The reflection column reads 0.0, 1.9, 7.3, 24.5, 46.8%. In the code $C = 1$, $L = 1$, and Mach is 0, so $K = \sigma \left( 1 - M^2 \right) \frac{c}{L}$ collapses to $K = \sigma$. Raising $\sigma$ from 0.25 to 10, a factor of 40, moves the reflection from 1.9 to 46.8. About 25 times. Not linear. $|R| = \frac{K}{\sqrt{K^2 + 4\omega^2}}$ saturates toward 1 as K grows. The pulse is not a single frequency. It is a Gaussian of width 0.09. Do not invert one $\omega$ out of these percentages. The direction is enough. Larger K reflects more, then the climb flattens.

The pressure column does not follow that monotone. In the run with $q = 0.3$, $\sigma = 0$ leaves a mean of +0.3000, the same size as the source. $\sigma = 0.25$ leaves +0.0007 and $\sigma = 1$ leaves +0.0006. Against 0.3000, 0.0007 is 0.23%. Up to there the relaxation nearly erases the offset. At $\sigma = 4$ the offset is +0.0227, about 32 times 0.0007. At $\sigma = 10$ it is −0.1142. The sign flips. The absolute value is about 160 times 0.0007. A boundary reflecting 46.8% does not hold the mean. It shakes it. The floor of the offset column sits at 0.25 and 1.0. The 0.25 Poinsot and Lele used is the left edge of that floor. 0.58 belongs to the B1 coefficient, so it does not land on this column as printed.

Waves the Grid Invented Travel the Wrong Way#

Physical acoustics are not the only thing reflecting off the boundary. Components with wavelengths shorter than about four mesh spacings are not solutions of Navier–Stokes at all; they are artifacts of the discretization. Poinsot and Lele called these "q waves" to separate them from the physical "p waves."

Their signature is group velocity. Even for the one-dimensional advection equation at speed $V$, the group velocity $u_g$ of short wavelengths has the opposite sign to $V$. The flow moves right, and the numerical error crawls upstream to the left. Worse, $|u_g / V|$ grows with the order of the spatial scheme, so high-order codes are more exposed, not less.

That is why a boundary treatment has to be judged with two reflection coefficients: $A_p / A_1$ for physical waves and $A_q / A_1$ for numerical ones. Any usable treatment needs $A_q / A_1 \ll 1$ under all circumstances; one that claims to be non-reflecting also needs $A_p / A_1 \ll 1$. Simply starting a run from an initial field with steep gradients generates q waves, and in a DNS nothing removes them afterward.

If the Period Is 2L/c, Suspect Outlet σ First#

$\sigma$ is not a tuning knob. One end is a pressure tied to nothing. The other end is a box that has become a resonator. Near 0.25 the boundary holds the low frequencies and lets the rest through.

An oscillation with no visible cause starts with a count: incoming characteristics versus imposed conditions. A period that is a multiple of $2L/c$ is the boundary, not the physics. Lower $\sigma$. If the oscillation shrinks, the boundary was the cause. If the mean pressure drifts, the run has crossed to the other side. Counting those waves before touching the outlet is the point of an earlier pass that sets NSCBC up as a non-reflecting outflow.