Converging nozzle from a UI-exported case

Load a network and its boundary conditions from a UI-produced YAML file and solve the mean flow.

Source notebook: examples/getting-started/converging_nozzle.ipynb. This page is executed from it at build time.

This notebook loads the network and boundary conditions from the YAML file produced by the user interface. We solve for the steady mean flow, and perform an acoustic scattering analysis on top of the converged mean flow.

The case is a converging nozzle with pipework: a high-pressure reservoir (a total-pressure inlet) discharges through a feed pipe, an isentropic area contraction, and a tailpipe to a fixed back pressure.

The lossless duct elements do not affect the mean flow, but their lengths model the wave propagation in the perturbation network that drives the phase.

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

import nefes
from nefes.plotting import use_nefes_theme, COLORWAY

# Use the bundled theme instead of plotly's default.
use_nefes_theme()

CASE = "converging_nozzle.yaml"

1. Load the case and solve

load_case parses the following YAML keys:

model.globalAttributes containing model level parameters
model.nodes containing node definitions
model.edges containing topological information and edge areas

The first run after a fresh installation may take some time because of the JIT compilation. The compiled objects are then cached, and unless they are manually erased won’t be compiled again.

handles and area) into a Network. Network.solve() compiles the immutable problem and runs the damped Newton solve from a quiescent start using default settings.

net = nefes.load_case(CASE)

# Network.solve() returns a Solution holding the converged state and metadata.
# Control the level of verbosity (default 0 - silent).
sol = net.solve(verbose=1)
assert sol.converged, (sol.residual_norm, sol.print_residuals())

print("converged :", sol.converged)
print("iterations:", sol.iterations)
print("||R_hat|| :", f"{sol.residual_norm:.2e}")

kappa=0.1   ->   8 iters, ||R_hat||=5.895e-14, converged=True
kappa=0.01  ->   2 iters, ||R_hat||=2.020e-15, converged=True
kappa=0     ->   1 iters, ||R_hat||=5.270e-12, converged=True
converged : True
iterations: 20
||R_hat|| : 3.05e-12

sol

Converged: True

Iterations: 20

Residual norm: 3.0455924980651284e-12

2. Inspect the converged edge states

Every edge carries the recovered state (mdot, p, h_t, rho, u, T, c, M, p_t). Total pressure is uniform across the lossless contraction; the throat and tailpipe (area 0.010) run faster than the feed pipe (area 0.020).

# Read every edge state by name; total pressure is uniform across the lossless
# contraction, and the throat (area 0.010) runs faster than the feed pipe (0.020).
sol.print_states()

# The choking edge is simply the fastest one.
throat = int(np.argmax(sol.field("M")))
print(f"\nthroat = edge {throat}  (M = {sol.edge(throat)['M']:.3f})")

edge	mdot [kg/s]	p [Pa]	h_t [J/kg]	rho [kg/m^3]	u [m/s]	T [K]	c [m/s]	M [-]	p_t [Pa]	area [m^2]
0	4.1247	1.9035e+05	3.0135e+05	2.2422	91.978	295.79	344.74	0.2668	2e+05	0.02
1	4.1247	1.9035e+05	3.0135e+05	2.2422	91.978	295.79	344.74	0.2668	2e+05	0.02
2	4.1247	1.5e+05	3.0135e+05	1.8914	218.08	276.33	333.21	0.65447	2e+05	0.01
3	4.1247	1.5e+05	3.0135e+05	1.8914	218.08	276.33	333.21	0.65447	2e+05	0.01


throat = edge 2  (M = 0.654)

3. Back-pressure sweep: emergent choking

Lower the outlet static pressure and re-solve. Above the critical ratio (p/p_t \approx 0.528 for air) the flow is subsonic and the mass flow rises as the back pressure drops. Below it the throat chokes: it reaches M=1 and the mass flow saturates at the analytic maximum

\dot m_{\max} = \frac{p_t}{\sqrt{T_t}}\,A^\ast\,\sqrt{\tfrac{\gamma}{R}}\left(\tfrac{2}{\gamma+1}\right)^{\frac{\gamma+1}{2(\gamma-1)}}.

No regime switch is coded: choking emerges from the smooth complementarity row. Each sweep point solves a with_params copy of the loaded case (the base stays pristine), addressing the outlet’s pressure by name ("back-pressure.p") and warm-starting from the previous converged state. nefes.parameter_study packages this loop; it is written out here to show the idiom.

pt, Tt, R, gamma = 200000.0, 300.0, 287.0, 1.4
A_throat = net.get("throat.area")
flux_star = np.sqrt(gamma / R) * (2.0 / (gamma + 1.0)) ** ((gamma + 1.0) / (2.0 * (gamma - 1.0)))
mdot_max = pt / np.sqrt(Tt) * flux_star * A_throat
p_crit = (2.0 / (gamma + 1.0)) ** (gamma / (gamma - 1.0)) * pt

ratios = np.linspace(0.30, 0.95, 26)
mdots, machs = [], []
prev = None
for r in ratios:
    s = net.with_params({"back-pressure.p": r * pt}).solve(x0=prev.x if prev else None)
    assert s.converged, (s.residual_norm, s.print_residuals())
    prev = s  # warm-start the next point (parameter writes never change the layout)
    est = s.edge(throat)
    mdots.append(est["mdot"])
    machs.append(est["M"])

mdots, machs = np.array(mdots), np.array(machs)
print(f"throat area         = {A_throat:.4f} m^2")
print(f"mdot_max (analytic) = {mdot_max:.4f} kg/s")
print(f"critical ratio      = {p_crit/pt:.4f}")

throat area         = 0.0100 m^2
mdot_max (analytic) = 4.6671 kg/s
critical ratio      = 0.5283

fig = make_subplots(
    rows=1,
    cols=2,
    subplot_titles=("Mass-flow saturation (choking)", r"$\text{Throat reaches } M = 1$"),
    horizontal_spacing=0.12,
)

fig.add_trace(go.Scatter(x=ratios, y=mdots, mode="lines+markers", name="solved"), row=1, col=1)
fig.add_hline(
    y=mdot_max,
    line=dict(dash="dash", color=COLORWAY[4]),
    annotation_text=r"$\dot{m}_{\max}$",
    annotation_position="bottom right",
    row=1,
    col=1,
)
fig.add_vline(
    x=p_crit / pt,
    line=dict(dash="dot", color="#52606d"),
    annotation_text="critical ratio",
    annotation_position="top left",
    row=1,
    col=1,
)

fig.add_trace(
    go.Scatter(x=ratios, y=machs, mode="lines+markers", name="throat M", line=dict(color=COLORWAY[2])), row=1, col=2
)
fig.add_hline(y=1.0, line=dict(dash="dash", color=COLORWAY[4]), row=1, col=2)
fig.add_vline(x=p_crit / pt, line=dict(dash="dot", color="#52606d"), row=1, col=2)

fig.update_xaxes(title_text=r"$\text{back-pressure ratio } p_\mathrm{out}/p_t$", row=1, col=1)
fig.update_xaxes(title_text=r"$\text{back-pressure ratio } p_\mathrm{out}/p_t$", row=1, col=2)
fig.update_yaxes(title_text=r"$\text{throat mass flow } \dot{m}\;[\mathrm{kg/s}]$", row=1, col=1)
fig.update_yaxes(title_text=r"$\text{throat Mach number } M$", row=1, col=2)
fig.update_layout(
    height=400, width=900, showlegend=False, hovermode="x", title_text="Emergent choking of the converging nozzle"
)
fig.show()

4. Perturbation network: transfer & scattering matrices

The same compiled network supports a linear perturbation analysis on top of the mean flow (theory §12). The 1-D Euler system carries three characteristics : two acoustic (f, g) and one entropy / convected wave (h): so the perturbation network has the same variable count as the mean flow, and with the entropy wave its transfer / scattering matrices are genuinely 3×3 (growing further with reacting scalars).

perturbation_response freezes the converged mean state, builds the perturbation operator A(\omega), and forces it once per frequency. Every incoming wave is prescribed: the two acoustic terminals on their boundary rows, the incoming entropy wave on the inlet edge’s total-enthalpy transport row (the edge view of nodal energy conservation): so nothing is ever left floating at a boundary. The excite argument chooses which families are driven: the default ("acoustic",) drives the acoustic waves and pins the incoming entropy to zero (a clean, well-conditioned 2×2 acoustic response); adding "entropy" drives the entropy wave too for the full 3×3. Any transfer matrix between an edge pair, or the scattering matrix between the terminals, is then reconstructed for the whole frequency array without re-solving.

We drive the full set, read the 3\times3 transfer matrix T(\omega) (in the characteristic flavor f,g,h) and the 3\times3 scattering matrix S(\omega), plus the clean acoustic 2×2: then plot them with the nefes.plotting complex-matrix viewers (magnitude over phase, one cell per entry).

from nefes.perturbation import find_terminals, verify_acoustic

verify_acoustic(sol)  # subsonic, flow-aligned, length > 0: preconditions OK

terms = find_terminals(sol)
inlet = next(t for t in terms if t.at_tail).edge  # 'feed'     (f enters; carries entropy in)
outlet = next(t for t in terms if not t.at_tail).edge  # 'tailpipe' (g enters)
print(f"inlet  = edge {inlet}")
print(f"outlet = edge {outlet}")

freq = np.linspace(20.0, 1500.0, 600)  # Hz
# Drive the entropy wave too for the full 3x3.  The incoming entropy is pinned to
# zero on the acoustic columns, so the acoustic sub-block stays clean either way.
resp = sol.perturbation_response(freq, excite=("acoustic", "entropy"))

T = resp.transfer_matrix(inlet, outlet)  # (n_freq, 3, 3) in (f, g, h)
S = resp.scattering_matrix(inlet, outlet)  # (n_freq, 3, 3) incoming -> outgoing
incoming, outgoing = resp.scattering_labels(inlet, outlet)
print("T, S shapes :", T.shape, S.shape)
print("S incoming  :", incoming, "\nS outgoing  :", outgoing)

# Read off the familiar acoustic coefficients AND the entropy-noise coupling.
R_in = S[:, 0, 0]  # reflection   f_in -> g_in
T_fwd = S[:, 1, 0]  # transmission f_in -> f_out
entropy_sound = S[:, 1, 1]  # entropy noise: incoming entropy -> downstream sound
print(f"|T_fwd|         : {np.abs(T_fwd).min():.3f} .. {np.abs(T_fwd).max():.3f}")
print(f"|entropy->f_out|: {np.abs(entropy_sound).min():.3f} .. {np.abs(entropy_sound).max():.3f}")

inlet  = edge 0
outlet = edge 3
T, S shapes : (600, 3, 3) (600, 3, 3)
S incoming  : [('a', 0), ('a', 2), ('b', 1)] 
S outgoing  : [('a', 1), ('b', 0), ('b', 2)]
|T_fwd|         : 1.103 .. 1.103
|entropy->f_out|: 11.010 .. 11.010

from nefes.plotting import plot_transfer_matrix, plot_complex_matrix

# (a) the full 3x3 transfer matrix: two acoustic waves AND the entropy wave.
# Grid layout: one (magnitude-over-phase) cell per entry, labelled f->f, f->g, ...
fig = plot_transfer_matrix(
    T,
    freq,
    edges=(inlet, outlet),
    x_title=r"$f\;[\mathrm{Hz}]$",
    width=920,
    height=720,
    title=f"Perturbation transfer matrix: edge {inlet} -> edge {outlet}  "
    f"(mean throat M ~ {sol.edge(throat)['M']:.2f})",
)
fig.show()

# (b) the clean acoustic 2x2 scattering matrix (incoming entropy pinned to zero): # flat preset (magnitude row over phase row, the familiar Bode array). Magnitudes
# are flat (lossless ducts are pure delays); only the phases wind down with frequency.
S2 = resp.acoustic_scattering_matrix(inlet, outlet)  # (n_freq, 2, 2)
fig2 = plot_complex_matrix(
    S2,
    freq,
    layout="flat",
    x_title=r"$f\;[\mathrm{Hz}]$",
    width=1000,
    height=420,
    title="Acoustic scattering matrix (rows: g_in, f_out;  cols: f_in, g_out)",
)
fig2.show()

The transfer-matrix magnitudes are flat in frequency: between these two stations the only frequency dependence is the duct propagation, and in wave variables a lossless duct is a pure delay: a unit-magnitude phase: so it winds each entry’s phase down linearly (the sawtooth ramps) without ever changing its amplitude. The entropy column (\,\cdot \leftarrow h) is large: an incoming temperature/density spot converts strongly to sound as it accelerates through the contraction: the entropy noise an acoustics-only (2×2) model cannot see: and its phase winds at the slow convective delay L/u, not the acoustic L/(c\pm u). The acoustic 2×2 block is the clean acoustic scattering matrix with the incoming entropy pinned to zero, so it is flat in magnitude; nothing is left floating at the boundary to fold the (large) entropy→sound coupling back into the acoustic waves. The whole picture comes from a single forced solve per frequency, re-used to extract any edge pair, all on the same converged mean flow and Jacobian from §1.

Export for Nemo

The full network is available as network (a Network) and its converged mean flow as sol (a Solution). Save either to a UI-readable YAML: sol.to_yaml(path) embeds the mean-flow result fields, network.to_yaml(path) writes the topology only: then open the file in Nemo.

import os, tempfile

_out = os.path.join(tempfile.mkdtemp(), "converging_nozzle.yaml")
sol.to_yaml(_out)  # embeds the mean-flow results; use net.to_yaml(_out) for topology only
print("saved case:", _out)

saved case: /tmp/tmp4ji13ewn/converging_nozzle.yaml