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. Today we look at what a boundary can actually compute, what it has to invent, and what a single coefficient controlling that invention costs you.

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$.

The whole game is in 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.

Push the Mach number in the diagram below. Which way each of the five characteristics crosses the boundary changes as you drag.

Raise $M$ from $+0.3$ to $+1.4$ and the incoming $\mathcal{L}_1$ flips direction, dropping the required condition count from one to zero. Flip the sign negative and the same face becomes an inlet, where the count jumps to four. If a code that runs happily at subsonic outflow diverges when the exit goes supersonic, it is usually imposing one condition more than this ledger permits.

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$.

Try it in the simulation below. The duct is closed on the left and open on the right. Hit fire pulse to launch a pressure wave, then move the sigma slider.

Watch two things. First, when the pulse reaches the outlet, does the red curve — the incoming amplitude $A^-$ — rise? That bump is the reflection. Second, turn up source q and watch the pressure history at the bottom: at $\sigma = 0$ the white curve settles above the $p_\infty$ line and stays there. You killed the reflection and lost the pressure.

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

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.

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.

One Line for Choosing an Outlet Condition#

$\sigma$ is not a tuning knob. It is a coordinate between two failures — an unanchored pressure at one end, a resonating computational box at the other. The reason 0.25 keeps getting recommended is that it grabs the low frequencies and passes everything else.

When an outlet starts oscillating for no visible reason, work in this order. Count the incoming characteristics at that face and check that the number of conditions you impose matches. Then check whether the oscillation period is a multiple of $2L/c$ — if it is, you are looking at boundary reflection, not physics. Finally, lower $\sigma$. If the oscillation shrinks, the boundary was the cause; if the mean pressure starts wandering off, you have crossed into the failure on the other side.