A gas passage cannot carry an arbitrarily large mass flow. As the back pressure is lowered, the flow in the narrowest section accelerates until it reaches the local speed of sound, M = 1, and the mass-flow rate then stops increasing however much further the outlet is drawn down. This saturation is choking, the most consequential compressible-flow phenomenon in network analysis — the operating principle of orifice flow meters, control valves, and rocket nozzles alike.
In this framework choking is not switched on by a regime flag; it is emergent. The either/or of “subsonic and lossless” versus “sonic and capped” is encoded once, as a complementarity condition inside the area-change and outlet rows already written in elements, and the network then discovers by itself which passage chokes and caps the mass flow there, with no regime declaration and no branching logic. This document derives that mechanism and exhibits the operating map it produces, computed directly by the solver. It restricts itself, in keeping with the overview scope, to subsonic operation up to a sonic throat; the supersonic and shock-bearing regimes are deferred, and the closing section states precisely which part of the physics that defers and why.
1.1 The choking limit
For a given upstream stagnation state, the mass flux through a passage is bounded, and the bound is reached exactly when the passage becomes sonic. The maximum dimensionless mass flux and the pressure ratio at which it occurs are given as:
where \gamma is the ratio of specific heats, R the specific gas constant, A the passage area, and p_t and T_t the upstream total pressure and temperature. The second relation is the critical pressure ratio: while the ratio of back pressure to total pressure stays above it the passage is subsonic and the mass flow rises as the back pressure falls, and once it drops to the critical value the passage is sonic and the mass flow can rise no further.
1.2 The throat is always an edge
A structural convention of the model makes the sonic state easy to place: every (atomic) element is taken to be internally monotone in area, interpolating between its port areas and hiding no interior extremum. The minimum area along any flow path is therefore an edge area, so any throat coincides with an edge, and the natural home of the sonic state is that edge — a state the recovery of state and recovery returns without any singularity, since the density root remains single and regular through M = 1. Consequently the choking constraint is attached to the existing elements adjacent to their small-port edges, rather than to a special nozzle element with hidden internal geometry; the physics is expressed on the geometry the network already carries.
1.3 Emergent choking via complementarity
Two facts conspire to make choking invisible to the plain isentropic equation, and understanding them is what motivates the complementarity form. First, the choking inequality — that the dimensionless mass flux cannot exceed its sonic value — holds identically on the edge-state manifold, so it can never be violated and therefore cannot itself act as a limiter. Second, at the choked state the total-pressure coupling of the isentropic element becomes tangent to zero, a fold: a double root at which the Jacobian is singular, observed as a smallest singular value collapsing to \sigma_{\min} \sim 5\times10^{-5} with over-driven cases stalling pinned at M = 1. The missing ingredient is the explicit either/or, supplied as a complementarity between the subsonic margin and the total-pressure drop, given as:
where M^{\text{in}}_{\text{small}} is the Mach number at the small port oriented into the element and the symbol \perp requires that at least one of the two non-negative quantities vanish. The condition reads directly: either the small port is subsonic, M < 1, and the total-pressure drop is zero — the classical lossless element, valid in either flow direction — or the small port is exactly sonic, M = 1, and a total-pressure drop becomes admissible, the lumped normal shock standing in the diverging continuation.
An important and easily missed point is that the distance from choking must be measured as 1 - M, a regular zero that is linear in the approach to sonic, and not as the mass-flow deficit \dot m^\ast - \dot m formed from the same edge. The latter vanishes quadratically at the sonic point — it is the fold in disguise — and using it would re-create the very singularity the complementarity is meant to remove; with the 1 - M measure the choked state is a regular equation and the smallest singular value at a choked solution improves by more than an order of magnitude (tests: test_choked_nozzle_saturates_mass_flow, test_critical_pressure_ratio_is_the_knee).
The either/or is realized as a single smooth residual by the Fischer–Burmeister complementarity function, given as:
\varphi_\varepsilon(a, b) = a + b - \sqrt{a^2 + b^2 + \varepsilon^2},
where a and b are the two non-negative quantities above and \varepsilon is a small, fixed smoothing. On the smoothed root manifold 2ab = \varepsilon^2, so the inactive branch is pinned to a quadratically small bias, and \varepsilon must be kept small and fixed — the bias \tfrac{1}{2}\varepsilon^2 p_t is a fictitious total-pressure loss that must stay far below the smallest driving pressure difference in the network, which is why \varepsilon is not relaxed along the artificial-resistance continuation of well-posedness. The same complementarity governs a converging pressure outlet, whose exit chokes at exactly M = 1 and whose static pressure then detaches upward from the specification — the underexpanded discharge of a choked orifice, whose expansion to the ambient occurs outside the network (test: test_choked_nozzle_outlet_critical_mass_flux).
1.4 The validated operating map
The whole of the subsonic-to-choked behaviour follows from these rows without any regime input, and it is best seen in the operating map of a converging nozzle as the back pressure is swept down, shown in Figure 1. Every point is an independent converged solve of the same two-element network — a reservoir feeding an isentropic area change that discharges through a pressure outlet — so the curves are the solver’s own output, not an analytic overlay.
Code
prs = np.linspace(0.95, 0.30, 28)M, mdot, pexit = np.array([solve_nozzle(pr) for pr in prs]).Tfig = go.Figure()fig.add_trace(go.Scatter(x=prs, y=M, mode="lines+markers", line=dict(color=COLORWAY[0]), marker=dict(size=5), name="exit Mach M"))fig.add_trace(go.Scatter(x=prs, y=mdot / mdot_max, mode="lines+markers", line=dict(color=COLORWAY[1]), marker=dict(size=5), name="mass flow ṁ / ṁ*"))fig.add_trace(go.Scatter(x=prs, y=pexit / PT, mode="lines+markers", line=dict(color=COLORWAY[2]), marker=dict(size=5), name="exit pressure p / p_t"))fig.add_vline(x=crit_pr, line=dict(color="#888888", dash="dash", width=1))fig.add_annotation(x=crit_pr, y=1.18, text="critical ratio ≈ 0.528", showarrow=False, font=dict(size=12, color="#888888"), xanchor="left")fig.update_layout(xaxis_title="back-pressure ratio p_b / p_t", yaxis_title="")fig.update_xaxes(autorange="reversed") # back pressure falls to the rightfig.update_yaxes(range=[0.0, 1.25])fig
Figure 1: Operating map of a converging nozzle as the back-pressure ratio p_b/p_t is lowered, every point an independent solve of the same network. The exit Mach number (top) rises to unity and then pins; the mass flow (middle), normalized by its sonic value \dot m^\ast, rises and then saturates at 1; the exit static pressure (bottom) follows the back pressure while subsonic and then detaches, holding at the critical value p^\ast/p_t \approx 0.528 as the back pressure keeps falling — the underexpanded choked discharge. The knee of all three curves is the critical pressure ratio (dashed), and nothing in the network was told where it is.
The three regimes are all emergent from the one set of rows, and each matches the classical construction:
Above the critical ratio — the passage runs as a subsonic venturi: lossless, with mass flow below the sonic value and the exit pressure equal to the back pressure (validated against the exact isentropic relation, test_subsonic_nozzle_matches_isentropic).
At and below the critical ratio — the throat edge reaches M = 1, the mass flow saturates at \dot m^\ast and is never exceeded, and the saturation is continuous and monotone through the knee (test_choked_nozzle_saturates_mass_flow, test_critical_pressure_ratio_is_the_knee).
Converging exit below critical — the exit behaves as a choked orifice: sonic at M = 1, mass flow at \dot m^\ast, and the exit static pressure detached upward to p^\ast, independent of the specified back pressure (test_choked_nozzle_outlet_critical_mass_flux).
Beyond these lies the shock regime — a normal shock standing in a diverging passage downstream of a supersonic patch — which requires supersonic flow inside the domain and is outside the present subsonic scope. It is worth noting how it would be recovered: a normal shock conserves mass and total enthalpy, so the converged total-pressure loss determines it in post-processing, and the shock’s Mach number and location are deduced from the Rankine–Hugoniot relations rather than imposed as unknowns — the location is an output, not an input.
One subtlety follows from the complementarity being written on the total-pressure loss rather than on the Mach number directly: the lossless branch is satisfied by the supersonic isentropic state (M > 1, no loss) as well as the subsonic one, so a purely isentropic supersonic branch is also a root of the steady rows. That branch is non-physical within the subsonic scope — a lossless duct cannot carry a supersonic edge into a subsonic one at equal total pressure — but a cold start can occasionally land on it, since the two branches share the same stagnation state. The solver therefore carries a subsonic-scope guard: after a converged solve it checks for a genuinely supersonic edge (beyond the sonic-throat band) and, if one is present, re-solves once from a near-stagnation seed that reliably reaches the physical subsonic branch, so a bare solve() returns the subsonic root without any hand-seeding. The guard is on by default and can be turned off with nefes.config.enforce_subsonic = False (or per solve, solve(enforce_subsonic=False)) for the deferred supersonic work; the choking mechanism itself is untouched — a real throat still pins at M = 1.
Within the subsonic scope, the choking mechanism of this document is complete: the network finds the choke point, caps the flow, and matches the classical operating map, all from the element rows written for the mean flow.
The same rows, linearized at the operating point, form the acoustic layer, where the compact choked nozzle contributes its Marble–Candel reflection, in perturbation network.