Source notebook: examples/thermoacoustics/rijke_tube.ipynb. This page is executed from it at build time.
A Rijke tube is the fundamental thermoacoustic oscillator: a duct with a heat source that, under the right conditions, features a self-excited thermoacoustic instability. It is the simplest system in which an active element (a heat-releasing element that responds to the acoustic field) feeds energy into a duct mode and drives it unstable. This notebook demonstrates the capability of Nefes to model this fundamental case in thermoacoustics.
We model the dynamic response of the flame with a flame transfer function \mathcal{F}(\omega) relating unsteady heat release to upstream velocity fluctuations:
\frac{\dot{Q}'(\omega)}{\overline{ \dot{Q}}} \;=\; \mathcal{F}(\omega)\,\frac{u'_{\mathrm{ref}}(\omega)}{\overline u_{\mathrm{ref}}}.
Here we take the classic n-\tau model,
\mathcal{F}(\omega) = n\,e^{-i\omega\tau},
so the heat release responds to the velocity an instant \tau earlier, with a constant gain n. The reference position is the acoustic “state” immediately before the flame.
In Nefes this is modeled in the source matrix \mathbf S(\omega) of the perturbation operator \mathbf A(\omega) = \mathbf J_{\mathrm{alg}} + i\omega\mathbf M + \mathbf P(\omega) + \mathbf S(\omega). The mean flame is acoustically passive (does not respond to acoustic fluctuations); attaching a DynamicSource makes \mathbf A non-self-adjoint, and the stability eigenproblem \det\mathbf A(\omega)=0 acquires complex roots: modal frequency \mathrm{Re}(\omega)/2\pi and growth rate -\mathrm{Im}(\omega).
In this notebook, we will (1) build the mean flow, (2) attach an n-\tau flame, (3) find the unstable mode, (4) map the stability boundary against the analytical compact-flame dispersion relation, and (5) repeat with the equilibrium flame model driven by the same FTF.
import warnings
import numpy as np
import plotly.graph_objects as go
import nefes
from nefes.elements import catalog as cat
from nefes.elements import n_tau, n_tau_flame
from nefes.perturbation import PerturbationBC, forced_power_balance
from nefes.plotting import use_nefes_theme, plot_transfer_function
use_nefes_theme()'nefes'Cold air enters a rigid (\dot{m}'=0) mass flow inlet, flows down through a duct of length L_1, crosses an acoustically compact flame, and leaves through a duct of length L_2 to an open end (p'=0). We keep the mean Mach number low so the acoustics are well represented by the quiescent conditions. This is not due to a limitation of Nefes, but to aid the analytical comparison.
The flame is a perfect-gas heat-release flame raising the total temperature by \Delta T = \dot Q/(\dot m c_p). An n–\tau flame transfer function is attached to it as a dynamic source; the steady solve ignores it (a constant mean source is acoustically passive), so the mean flow is the same with or without it.
# Constants that define the perfect gas model
R, GAMMA = 287.0, 1.4
CP = GAMMA * R / (GAMMA - 1.0)
# Tube cross-section
AREA = 0.01
# Cold (upstream) and hot (downstream) duct lengths
L1, L2 = 0.6, 0.4
# Mean mass flow rate through the duct
MDOT = 0.006
# Cold inlet temperature
T_IN = 300.0
# Open-end (and reference) static pressure
P_OUT = 1.0e5
# Mean temperature rise across the flame, used here to size the flame heat release rate.
DT_RISE = 400.0
# Integrated heat release rate this compact flame will produce [W]
QDOT = MDOT * CP * DT_RISE
def rijke(
n,
tau,
*,
R=R,
gamma=GAMMA,
area=AREA,
L_cold=L1,
L_hot=L2,
mdot=MDOT,
qdot=QDOT,
T_inlet=T_IN,
p_outlet=P_OUT,
drive=False,
):
"""Build and solve the Rijke tube with a compact n-tau flame.
Every geometric, mean-flow, and gas parameter is an explicit argument (defaulting to the
module-level constants), and the network is assembled node by node so the
inlet -> duct -> flame -> duct -> outlet topology is spelled out. The flame's velocity
reference is the edge just upstream of it; rather than guess that index up front, we
take it from the value ``connect`` returns and attach the n-tau source afterwards with
``set_dynamic_source``.
Parameters
----------
n : float
Flame-transfer-function gain of the n-tau source.
tau : float
Flame time lag [s] of the n-tau source.
R : float
Specific gas constant [J/(kg K)] of the perfect gas.
gamma : float
Heat-capacity ratio of the perfect gas.
area : float
Tube cross-sectional area [m^2] (uniform along the duct).
L_cold : float
Cold (upstream) duct length [m].
L_hot : float
Hot (downstream) duct length [m].
mdot : float
Mean mass flow rate through the duct [kg/s].
qdot : float
Integrated mean heat release rate of the compact flame [W].
T_inlet : float
Cold inlet temperature [K].
p_outlet : float
Open-end static pressure [Pa]; also the reference pressure.
drive : bool
When ``True`` the open end is additionally driven by a unit incoming acoustic wave;
the convective open-end reflection -(1-M)/(1+M) is kept, so the tube's modes are unchanged
and only a forcing is added: this is what the forced-response frequency sweep excites.
The mean flow is identical to ``drive=False``. Default ``False``.
Returns
-------
Solution
The converged mean-flow solution.
"""
net = nefes.Network(nefes.perfect_gas(R, gamma))
# The rigid inlet pins u' = 0 and the open end pins p' = 0, so the tube's ends are
# energy-neutral reflectors. When ``drive`` is set the open end also injects a unit
# incoming wave on top of its reflection, exciting the tube without altering its modes.
outlet_bc = PerturbationBC.open_end(driven=("acoustic",) if drive else ())
inlet_bc = PerturbationBC.hard_wall()
idx_inlet = net.add(cat.mass_flow_inlet(mdot, T_inlet, name="inlet", perturbation_bc=inlet_bc))
idx_cold = net.add(cat.duct(L_cold, name="cold duct"))
idx_flame = net.add(cat.heat_release_flame(qdot, name="flame"))
idx_hot = net.add(cat.duct(L_hot, name="hot duct"))
idx_outlet = net.add(cat.pressure_outlet(p_outlet, name="open end", perturbation_bc=outlet_bc))
# Assemble the network: cold air -> duct -> n-tau heat-release flame -> duct -> open end.
# connect() returns the new edge's id; the cold -> flame edge is the flame's velocity reference.
net.connect(idx_inlet, idx_cold, area)
idx_ref = net.connect(idx_cold, idx_flame, area)
net.connect(idx_flame, idx_hot, area)
net.connect(idx_hot, idx_outlet, area)
# Attach the n-tau flame transfer function now that the reference edge id is known.
net.set_dynamic_source(idx_flame, n_tau_flame(n, tau, ref_edge=idx_ref))
sol = net.solve()
assert sol.converged
return sol
sol0 = rijke(0.0, 0.0)
# edge 1 is just upstream of the flame; edge 2 the hot products just downstream
sol0.print_states()| 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 | 0.006 | 1e+05 | 3.0135e+05 | 1.1614 | 0.5166 | 300 | 347.19 | 0.0014879 | 1e+05 | 0.01 |
| 1 | 0.006 | 1e+05 | 3.0135e+05 | 1.1614 | 0.5166 | 300 | 347.19 | 0.0014879 | 1e+05 | 0.01 |
| 2 | 0.006 | 1e+05 | 7.0315e+05 | 0.49776 | 1.2054 | 700 | 530.34 | 0.0022729 | 1e+05 | 0.01 |
| 3 | 0.006 | 1e+05 | 7.0315e+05 | 0.49776 | 1.2054 | 700 | 530.34 | 0.0022729 | 1e+05 | 0.01 |
The n-\tau model is a pure gain n and a delay \tau: a flat magnitude and a phase that winds linearly with frequency. The plot_transfer_function helper offers tools to plot a complex valued transfer function.
N, TAU = 0.8, 4.0e-3 # gain and time lag [s]
ftf = n_tau(N, TAU)
freqs = np.linspace(0.0, 400.0, 400)
plot_transfer_function(
ftf,
freqs,
names=[f"n-tau (n={N}, tau={TAU*1e3:.1f} ms)"],
title=r"Flame transfer function",
).show()A self-excited mode is a complex root of \det\mathbf A(\omega)=0. Before hunting those roots directly, we perform a forced response analysis: drive one end of the tube and sweep a real frequency, watching two things - the acoustic energy stored inside the duct, and the net acoustic power crossing the boundaries. Directly hunting for complex roots may become very challenging for difficult problems, and doing a frequency sweep provides an idea about the landscape. Resonances show up as energy peaks; the boundary-power trace shows whether the domain is absorbing the drive or feeding energy back out. This locates the modes (and already reveals the flame’s effect) without solving an eigenvalue problem.
We then pin each mode with eigenmodes, which solves \det\mathbf A(\omega)=0 by Beyn’s contour-integral method, returning the exact complex frequency and growth rate. For this example we use the isentropic operator, e.g. \rho' = p'/\overline{c}^2 (the standard acoustic network model) so the flame’s energy jump stays physical while entropy convection is dropped. The boundaries are energy-neutral, so a positive growth rate is a genuinely self-excited (unstable) mode.
We drive the open end with a unit incoming acoustic wave (keeping its convective open-end reflection R=-(1-M)/(1+M), so the tube’s own modes are unchanged, we only excite them) and solve the forced field at each real frequency with forced_response. forced_power_balance then returns the node-wise acoustic-energy budget, frequency by frequency.
Conservation of the acoustic energy flux at every node gives a real-valued ledger: the power produced inside the domain plus the net flux across the boundaries equals the storage rate, which is zero on average for a steady forced state. So the budget closes, G + F_\mathrm{refl} + F_\mathrm{src} \;=\; 0, with three buckets:
The reflectors stay near zero and the drive mirrors the generation (F_\mathrm{src}\approx -G), so in this driven run energy crosses the boundary only at the excitation port, never through the reflecting ends. We unpack what that does (and does not) mean right after the plot.
freqs_sweep = np.linspace(60.0, 320.0, 600)
# forced_power_balance returns the node-wise energy budget at each frequency: stored energy E, the
# interior generation G (the flame's flux jump), and the boundary flux split into its passive
# reflector part (~0) and the excitation source (the drive). Steady state: G + reflection + source ~ 0.
sol_p = rijke(0.0, 0.0, drive=True)
sol_a = rijke(N, TAU, drive=True)
bal_p = forced_power_balance(sol_p.forced_response(freqs_sweep, isentropic=True), sol_p.problem)
bal = forced_power_balance(sol_a.forced_response(freqs_sweep, isentropic=True), sol_a.problem)
g_scale = np.abs(bal.generation).max() # shared scale for the three (active) power terms
fig = go.Figure()
fig.add_scatter(
x=freqs_sweep, y=bal_p.energy / bal_p.energy.max(), name="passive: stored energy E", line=dict(color="#9aa5b1")
)
fig.add_scatter(
x=freqs_sweep, y=bal.energy / bal.energy.max(), name="active: stored energy E", line=dict(color="#d62728")
)
fig.add_scatter(
x=freqs_sweep, y=bal.generation / g_scale, name="active: flame generation G", line=dict(color="#1f77b4")
)
fig.add_scatter(
x=freqs_sweep,
y=bal.boundary_source / g_scale,
name="active: boundary source (drive)",
line=dict(color="#1f77b4", dash="dot"),
)
fig.add_scatter(
x=freqs_sweep,
y=bal.boundary_reflection / g_scale,
name="active: boundary reflectors",
line=dict(color="#10b981", dash="dash"),
)
fig.add_hline(y=0.0, line_width=1, line_color="#888")
fig.update_layout(
title="Forced sweep: stored energy locates the modes; the budget is G + reflection + source = 0",
xaxis_title=r"$f\;[\mathrm{Hz}]$",
yaxis_title=r"$E,\;\; G,\;\; \text{boundary terms}\quad(\text{normalized})$",
)
# the reflectors carry ~0 and the drive mirrors G; the budget closes to numerical noise
closure = np.abs(bal.residual).max() / g_scale
print(f"active tube budget closure: max|G + reflection + source| / max|G| = {closure:.1e}")
fig.show()active tube budget closure: max|G + reflection + source| / max|G| = 7.2e-10The energy peaks are the tube’s resonances: our mode candidates. With the flame inert they sit at the bare-duct frequencies (\approx 111 and \approx 306 Hz); switching the flame on slides the fundamental up to \approx 136 Hz and broadens it.
The flame generation G is the physically meaningful trace; the net boundary flux is only its mirror, F_\mathrm{in}=-G. The inert flame is a passive scatterer: its mean density jump lets it only absorb (G<0): while the active flame produces acoustic power near its fundamental (G>0).
Now the subtle part: F_\mathrm{in}<0 means acoustic energy really does leave through the outlet in this driven run: it must, since a constant-amplitude forced state cannot store the flame’s output, so the G produced each cycle has to exit (F_\mathrm{in}=-G). How, when |R|\approx 1? The wave returning into the domain is g = R f + d: the reflection R f\approx -f and the imposed drive d\approx +f nearly cancel, so the returning wave is far weaker than the outgoing one (|g|\approx 0.3\,|f| here) and the net flux is outward. R fixes the reflection; it does not fix the total return once an independent wave is injected.
Undriven (d=0) the return is the full reflection, g=R f with |g|\approx|f|: the net flux vanishes and the same G instead accumulates: that trapping is the growing mode. So the outflow is a property of the driven state, not of the tube; keep the excitation separate from the reflection, and read stability off the drive-independent G: or, definitively, off the undriven eigenproblem below.
Passive baseline first (n=0):
# Passive baseline: with the flame inert (n=0) the network is acoustically passive, so these are the
# bare duct modes: the frequencies the active flame will act on.
# The growth rates are near-neutral; the small negative value is the O(M) mean-flow damping the
# energy-neutral ends leave behind.
sol_p = rijke(0.0, 0.0)
res_p = sol_p.eigenmodes(freq_band=(40.0, 320.0), growth_band=(-200.0, 200.0), isentropic=True)
for m in res_p.summary():
print(f"f = {m['freq_hz']:7.2f} Hz growth = {m['growth_rate']:+8.2f} 1/s")
res_p.plot_spectrum(title="Rijke spectrum: passive flame (n=0)").show()f = 111.16 Hz growth = -0.90 1/s
f = 306.11 Hz growth = -0.96 1/sTurning the flame on adds the source \mathbf S(\omega)=\overline{Q}\,\mathcal F(\omega) to the operator, with \mathcal F(\omega)=n\,e^{-i\omega\tau}. This is a complex, frequency-dependent term, so it does not simply add growth at the passive frequency: it moves the eigenvalue in the complex plane:
So the fundamental does not sit at the passive \approx 111\,\mathrm{Hz} and merely cross into instability: it shifts up to \approx 136\,\mathrm{Hz} as it goes unstable. This frequency pull is genuine thermoacoustics, not a numerical artefact: the analytic dispersion relation in §4 reproduces the same shifted, growing root at every \tau (and at n=0 collapses back onto the passive baseline above).
sol = rijke(N, TAU)
res = sol.eigenmodes(freq_band=(40.0, 320.0), growth_band=(-200.0, 200.0), isentropic=True)
for m in res.summary():
flag = " <-- UNSTABLE" if m["unstable"] else ""
print(f"f = {m['freq_hz']:7.2f} Hz growth = {m['growth_rate']:+8.2f} 1/s{flag}")
res.plot_spectrum(title=f"Rijke spectrum (n={N}, tau={TAU*1e3:.0f} ms)").show()f = 135.53 Hz growth = +57.18 1/s <-- UNSTABLE
f = 299.63 Hz growth = -37.13 1/sWe can show the trapping directly with the acoustic-power diagnostic. For a self-excited mode the energy balance is 2\sigma E = G + F_\mathrm{in}: the growth rate times the stored energy equals the flame generation plus the net boundary flux. boundary_power reads each terminal’s reflection magnitude |R| against the energy-neutral bound (1\pm M)/(1\mp M) and the net flux it passes into the domain. The rigid inlet (u'=0) and the open outlet (p'=0) sit at the energy-neutral magnitude, so F_\mathrm{in}\approx 0: only an O(M) convective leak crosses them. The flame’s generated energy is therefore trapped between the two reflectors and accumulates: that is the growth, 2\sigma E = G > 0, supplied by the flame alone.
# The boundaries are (energy-)neutral reflectors: for the unstable mode the net acoustic energy flux
# crossing the terminals is ~0, so the flame's generated energy is trapped between the ends and
# accumulates: that is the growth (each terminal is an acoustic node: the mass-flow inlet pins
# u'=0, the open outlet pins p'=0, so the Myers flux vanishes at the face but for an O(M) leak).
imode = int(np.argmax(res.growth_rates)) # the self-excited (unstable) mode
bp = res.boundary_power(imode)
print(f"unstable mode f = {bp.freq_hz:.2f} Hz growth = {bp.growth_rate:+.2f} 1/s\n")
print(f"{'terminal':<10}{'|R|':>8}{'energy-neutral |R|':>22}{'flux in':>14}")
for e in sorted(bp.entries, key=lambda d: d["edge"]):
print(f"{e['name']:<10}{e['reflection']:>8.3f}{e['passive_bound']:>22.3f}{e['power_in']:>14.1e}")
print(f"\nnet boundary flux = {bp.net:+.1e} (mode-scale units): ~0, the reflecting ends are")
print("energy-neutral: the +growth is the flame's generated energy trapped between them.")unstable mode f = 135.53 Hz growth = +57.18 1/s
terminal |R| energy-neutral |R| flux in
inlet 1.000 0.997 1.9e-15
open end 1.000 1.005 -2.4e-16
net boundary flux = +1.6e-15 (mode-scale units): ~0, the reflecting ends are
energy-neutral: the +growth is the flame's generated energy trapped between them.The same energy ledger yields the growth rate directly: no contour eigensolver. For a free mode the budget is 2\sigma E = G + F_\mathrm{in} (no excitation now), so \sigma \;=\; \frac{G + F_\mathrm{in}}{2E}, the generation plus the (near-zero) boundary flux, over twice the stored energy. res.energy_balance(mode) forms exactly this budget; its energy-derived growth rate must match the contour eigenvalue: an independent check on the solver, not a restatement.
# The same energy ledger, now read for the growth rate: sigma = (G + F_in) / (2 E).
eb = res.energy_balance(imode)
print(f"unstable mode f = {eb.freq_hz:.2f} Hz (mode-scale energy units)")
print(f" flame generation G = {eb.generation:+.3e}")
print(f" net boundary flux F_in = {eb.boundary_flux:+.3e} (~0: energy-neutral ends)")
print(f" stored energy E = {eb.stored_energy:.3e}")
print(f"\n growth rate (energy budget) sigma = (G + F_in)/(2E) = {eb.growth_rate_energy:+.3f} 1/s")
print(f" growth rate (contour eigenvalue) -Im(omega) = {eb.growth_rate:+.3f} 1/s")
print(f" agree: {eb.consistent}")unstable mode f = 135.53 Hz (mode-scale energy units)
flame generation G = +1.288e-13
net boundary flux F_in = +1.643e-15 (~0: energy-neutral ends)
stored energy E = 1.140e-15
growth rate (energy budget) sigma = (G + F_in)/(2E) = +57.184 1/s
growth rate (contour eigenvalue) -Im(omega) = +57.184 1/s
agree: TrueThe growth rate and frequency say when and how fast; the mode shape says where. Reconstructing the pressure field analytically inside each duct and animating it over one cycle shows the self-excited standing wave, with the compact flame marked at the cold/hot interface (x=L_1). The reconstruction is exact: the duct’s own phase relation evaluated at every interior station, not a discretization.
# Spatial shape of the self-excited mode: the pressure field reconstructed analytically inside each
# duct (f rides downstream at u+c, g upstream at c-u) and animated over a cycle. The flame sits at
# the cold/hot interface (x = L1); press play to watch the growing standing wave.
res.animate_mode(
imode,
variable="p",
title=f"Unstable Rijke mode: p' along the tube ({res.freqs[imode]:.0f} Hz)",
).show()Sweeping the time lag \tau traces the classic n–\tau stability band: the peak growth rate rises and falls, crossing zero into instability for a range of \tau.
We check Nefes against the analytical compact-flame dispersion relation: a two-region (cold/hot) acoustic tube with the zero-mean-Mach jump
p'_1 = p'_2, \qquad u'_2 - u'_1 = \frac{\gamma-1}{\gamma\,\overline{p}\,A}\,q',
and q' = \overline{Q}\,\mathcal{F}(\omega)\,u'_1/\overline{u}_1 with \mathcal{F}(\omega)=n\,e^{-i\omega\tau}.
Modal ansatz. Seek normal modes \hat\xi(x,t)=\Re\{\hat\xi(x)\,e^{\lambda t}\}. Nefes uses the harmonic convention e^{i\omega t}, so each eigenvalue is
\lambda = i\omega, \qquad \omega = \omega_r + i\omega_i \in \mathbb{C}.
The real and imaginary parts of \lambda are the quantities plotted in the spectrum:
In each uniform region j\in\{1,2\} (cold / hot) the spatial part is a pair of plane waves, k_j=\omega/c_j, Z_j=\rho_j c_j:
\hat p'_j(x) = Z_j\big(A_j^+ e^{ik_j x} + A_j^- e^{-ik_j x}\big),
with the flame at x=0. Matching the inlet hard wall (\hat u'_1=0 at x=-L_1), the open outlet (\hat p'_2=0 at x=L_2), pressure continuity, and the velocity jump above gives a 4\times4 system \mathbf M(\omega)\mathbf a=0 for the wave amplitudes \mathbf a. A nontrivial mode exists iff \det\mathbf M(\omega)=0 (equivalently \det\mathbf A(\omega)=0 in Nefes). Each root \omega_k yields an eigenvalue \lambda_k=i\omega_k; there is no closed form for \lambda_k in general, but Newton polish on \det\mathbf M recovers it.
At this low Mach the analytic and Nefes roots agree closely; the small residual is the O(M) mean-flow damping the zero-Mach analytic omits and the approximation \dot{m}' \approx u' used in the analytical evaluation.
def region_means(sol):
# The compact-flame analysis needs only the uniform mean state on each side of the flame.
# Edge 1 is the cold approach (region 1), edge 2 the hot products (region 2); the mean flow is the
# same with or without the FTF, so this is evaluated once on the inert (n=0) tube.
e1, e2 = sol.edge(1), sol.edge(2)
return dict(rho1=e1["rho"], c1=e1["c"], u1=e1["u"], p1=e1["p"], rho2=e2["rho"], c2=e2["c"])
def analytic_det(omega, m, Qdot, n, tau):
# Two-region (cold/hot) acoustic tube with a compact flame at x = 0. In each uniform region the
# field is a pair of plane waves; matching the four boundary/jump conditions gives a 4x4 system
# M(omega) a = 0 for the wave amplitudes a = [A1+, A1-, A2+, A2-]. A mode exists iff det M = 0.
k1, k2 = omega / m["c1"], omega / m["c2"] # wavenumbers (zero-Mach limit: k = omega / c)
Z1, Z2 = m["rho1"] * m["c1"], m["rho2"] * m["c2"] # characteristic impedances Z = rho c
# Plane-wave phase carried across each duct (cold length L1, hot length L2).
eP1, eM1 = np.exp(1j * k1 * L1), np.exp(-1j * k1 * L1)
eP2, eM2 = np.exp(-1j * k2 * L2), np.exp(1j * k2 * L2)
# Heat-release coupling. The velocity jump is driven by q' = Qdot * F(omega) * u'_1 / u_1 with the
# n-tau FTF F(omega) = n e^{-i omega tau}. Th gathers (gamma-1)/(gamma p A) * q', referred to the
# cold-side amplitudes through 1/(Z1 u1); it is what makes the row active (Th = 0 -> inert jump).
Th = (GAMMA - 1.0) * Qdot / (GAMMA * m["p1"] * AREA) * (n * np.exp(-1j * omega * tau)) / (Z1 * m["u1"])
# The four rows, top to bottom:
# (0) inlet hard wall u'_1 = 0 at the cold end,
# (1) pressure continuity p'_1 = p'_2 across the compact flame,
# (2) velocity jump u'_2 - u'_1 = (gamma-1)/(gamma p A) q' (the active row carrying Th),
# (3) open-end outlet p'_2 = 0 at the hot end.
M = np.array(
[[eP1, -eM1, 0, 0], [1, 1, -1, -1], [-(1 / Z1 + Th), (1 / Z1 + Th), 1 / Z2, -1 / Z2], [0, 0, eP2, eM2]],
dtype=complex,
)
return np.linalg.det(M)
def analytic_modes(m, Qdot, n, tau, fband=(40.0, 320.0)):
# Isolate the roots of det M in two stages: a coarse grid scan for seeds, then Newton polish.
# Stage 1: sample |det M| on a coarse omega = omega_r + i*omega_i grid over the search box.
wr = 2 * np.pi * np.linspace(*fband, 60) # real part: the angular-frequency band
wi = np.linspace(-220.0, 220.0, 60) # imag part: the growth-rate band (growth = -Im omega)
WR, WI = np.meshgrid(wr, wi)
D = np.abs(np.vectorize(lambda w: analytic_det(w, m, Qdot, n, tau))(WR + 1j * WI))
grid = WR + 1j * WI
# A grid node that is a local minimum of |det| (no smaller neighbour in its 3x3 stencil) sits near
# a root: take those as Newton seeds.
seeds = [
grid[i, j]
for i in range(1, D.shape[0] - 1)
for j in range(1, D.shape[1] - 1)
if D[i, j] == D[i - 1 : i + 2, j - 1 : j + 2].min()
]
out = []
for w0 in seeds:
# Stage 2: Newton on det M with a central-difference derivative (there is no closed form).
w = complex(w0)
for _ in range(60):
h = 1e-3 * (abs(w) + 1.0)
d = analytic_det(w, m, Qdot, n, tau)
dp = (analytic_det(w + h, m, Qdot, n, tau) - analytic_det(w - h, m, Qdot, n, tau)) / (2 * h)
if dp == 0:
break
s = d / dp
w -= s
if abs(s) < 1e-9: # converged on this root
break
# Keep converged roots inside the search box, dropping duplicates reached from adjacent seeds.
if 2 * np.pi * fband[0] < w.real < 2 * np.pi * fband[1] and -220.0 < w.imag < 220.0:
if all(abs(w - o) > 1e-1 for o in out):
out.append(w)
return out
# The mean flow ignores the FTF, so the two-region mean state is identical at every tau: build it once.
m_mean = region_means(rijke(0.0, 0.0))
taus = np.linspace(0.0, 6.0e-3, 25)
g_fns, g_an = [], []
for tau in taus:
# Nefes: the most unstable growth rate the contour eigensolver returns in the band.
sol = rijke(N, tau)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
r = sol.eigenmodes(freq_band=(40.0, 320.0), growth_band=(-260.0, 260.0), isentropic=True)
g_fns.append(max(r.growth_rates) if len(r) else np.nan)
# Analytic: the peak growth rate (-Im omega) over the roots of the compact-flame dispersion relation.
am = analytic_modes(m_mean, QDOT, N, tau)
g_an.append(max((-w.imag for w in am), default=np.nan))
fig = go.Figure()
fig.add_scatter(x=taus * 1e3, y=g_an, mode="lines", name="analytic")
fig.add_scatter(x=taus * 1e3, y=g_fns, mode="markers", name="Nefes eigenmodes", marker=dict(size=8))
fig.add_hline(y=0.0, line_dash="dash", annotation_text="stability boundary")
fig.update_layout(
title="Rijke stability map: peak growth rate vs. flame time lag",
xaxis_title=r"$\text{flame time lag } \tau\;[\mathrm{ms}]$",
yaxis_title=r"$\text{peak growth rate } \sigma\;[\mathrm{s^{-1}}]$",
)
fig.show()The Nefes markers sit on the analytic curve and cross the stability boundary at the same lags: the sign and magnitude of the heat-release coupling check out: the dynamic source genuinely drives the instability, not merely perturbs it.
The same FTF drives the reacting equilibrium flame: the source does not care which flame model supplies it. Here the mean flame temperature is not prescribed: stoichiometric H_2/air enters frozen and burns to chemical equilibrium across the flame (the reacting mean solver), and the mean heat release \overline{Q} for the FTF auto-derives from the converged temperature rise.
The reacting acoustics use the gas’s actual caloric coupling in the characteristic maps: the per-edge (\partial h/\partial\rho)_p,\ (\partial h/\partial p)_\rho from the equilibrium closure, not the perfect-gas c_p/R: so the reacting case is now quantitatively verified, not just demonstrated. Two physical differences from the fixed-power perfect-gas flame are worth noting:
We verify against the matching equilibrium-flame dispersion relation: five waves (the two acoustic + an entropy spot shed at the flame) with the exact equilibrium EOS derivatives and the same source coupling Nefes stamps.
from nefes.thermo import SpeciesDatabase, EQ_FROZEN, EQ_KERNEL
from nefes.chem import resolve_composition, enthalpy_mass
lib = SpeciesDatabase().select(["H2", "O2", "N2", "H2O", "OH", "H", "O", "HO2", "NO"]) # bundled CEA data
fuel_air = {"H2": 1.0, "O2": 0.5, "N2": 0.5 * 3.76} # stoichiometric, mole basis
Y, Z = resolve_composition(lib, fuel_air, basis="mole")
H_REACT = enthalpy_mass(lib, Y, 300.0)
MDOT_R = 0.02
def rijke_reacting(
n,
tau,
*,
area=AREA,
L_cold=L1,
L_hot=L2,
mdot=MDOT_R,
T_inlet=300.0,
p_outlet=1.0e5,
composition=fuel_air,
h_inlet=H_REACT,
):
"""Build and solve the reacting (equilibrium-flame) Rijke tube.
Same node-by-node topology as ``rijke``, but the flame temperature is not prescribed: the
reactant stream enters frozen (unburnt) and burns to chemical equilibrium across the flame.
The cold approach edges carry the frozen closure and the hot product edges the equilibrium
closure, set per edge via ``connect(..., edge_model=...)``. As in ``rijke``, the flame's
velocity reference is taken from the edge id ``connect`` returns and the n-tau source is
attached afterwards with ``set_dynamic_source``.
Parameters
----------
n : float
Flame-transfer-function gain of the n-tau source.
tau : float
Flame time lag [s] of the n-tau source.
area : float
Tube cross-sectional area [m^2] (uniform along the duct).
L_cold : float
Cold (upstream) duct length [m].
L_hot : float
Hot (downstream) duct length [m].
mdot : float
Mean mass flow rate through the duct [kg/s].
T_inlet : float
Reactant inlet temperature [K].
p_outlet : float
Open-end static pressure [Pa]; also the reference pressure.
composition : dict
Inlet (and backflow) reactant composition.
h_inlet : float
Absolute-enthalpy datum of the reactant stream [J/kg].
Returns
-------
Solution
The converged mean-flow solution.
"""
# unburnt reactant approach -> burnt products: frozen upstream, equilibrium downstream of the flame
net = nefes.Network(nefes.equilibrium(lib))
inlet = net.add(
cat.mass_flow_inlet(
mdot,
T_inlet,
composition=composition,
basis="mole",
name="inlet",
perturbation_bc=PerturbationBC.hard_wall(),
)
)
cold = net.add(cat.duct(L_cold, name="cold duct"))
flame = net.add(cat.equilibrium_flame(name="flame"))
hot = net.add(cat.duct(L_hot, name="hot duct"))
outlet = net.add(
cat.pressure_outlet(
p_outlet,
Tt_backflow=T_inlet,
composition=composition,
basis="mole",
name="open end",
perturbation_bc=PerturbationBC.open_end(),
)
)
# connect() returns the new edge's id; the cold -> flame edge is the flame's velocity reference
net.connect(inlet, cold, area, edge_model=EQ_FROZEN)
ref = net.connect(cold, flame, area, edge_model=EQ_FROZEN)
net.connect(flame, hot, area, edge_model=EQ_KERNEL)
net.connect(hot, outlet, area, edge_model=EQ_KERNEL)
# attach the n-tau flame transfer function now that the reference edge id is known
net.set_dynamic_source(flame, n_tau_flame(n, tau, ref_edge=ref))
sol = net.solve()
assert sol.converged
return sol
sol_r = rijke_reacting(0.0, 0.0)
# edge 1 is the unburnt reactant approach, edge 2 the burnt products across the (adiabatic) flame
sol_r.print_states()| 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 | 0.02 | 1.0003e+05 | 2572.6 | 0.83859 | 2.3849 | 300 | 408.7 | 0.0058355 | 1.0003e+05 | 0.01 |
| 1 | 0.02 | 1.0003e+05 | 2572.6 | 0.83859 | 2.3849 | 300 | 408.7 | 0.0058355 | 1.0003e+05 | 0.01 |
| 2 | 0.02 | 1e+05 | 2572.6 | 0.12271 | 16.299 | 2378.9 | 978.83 | 0.016651 | 1.0002e+05 | 0.01 |
| 3 | 0.02 | 1e+05 | 2572.6 | 0.12271 | 16.299 | 2378.9 | 978.83 | 0.016651 | 1.0002e+05 | 0.01 |
from nefes.thermo.api import thermo_state
def reacting_caloric(sol, e):
# (dh/drho)_p, (dh/dp)_rho at edge e via a complex step of the equilibrium closure
prob, x = sol.problem, sol.x
xi = np.ascontiguousarray(x[3 : 3 + prob.n_elem, e]).astype(np.complex128)
ht, p, d = complex(x[2, e]), float(sol.edge(e)["p"]), 1e-30
mid = int(prob.edge_model[e])
drho_dh = thermo_state(mid, prob.tf, prob.ti, xi, ht + 1j * d, complex(p))[1].imag / d
drho_dp = thermo_state(mid, prob.tf, prob.ti, xi, ht, p + 1j * d)[1].imag / d
return 1.0 / drho_dh, -drho_dp / drho_dh
def reacting_det(omega, mm, cal1, cal2, delta, n, tau):
# 5 waves [f1, g1, f2, g2, h2]: mass-flux + momentum + (h_t cont + source), equilibrium EOS
rho1, cc1, u1, p1, rho2, cc2, u2, p2 = mm
a1, b1 = cal1
a2, b2 = cal2
Z1, Z2 = rho1 * cc1, rho2 * cc2
F = n * np.exp(-1j * omega * tau)
k1p, k1m, k2p, k2m = omega / (u1 + cc1), omega / (cc1 - u1), omega / (u2 + cc2), omega / (cc2 - u2)
M = np.zeros((5, 5), complex)
M[0, 0], M[0, 1] = np.exp(1j * k1p * L1), -np.exp(-1j * k1m * L1) # inlet hard wall
M[1, 2], M[1, 3] = np.exp(-1j * k2p * L2), np.exp(1j * k2m * L2) # outlet open end
M[2] = [
AREA * (u1 * Z1 / cc1**2 + rho1),
AREA * (u1 * Z1 / cc1**2 - rho1),
-AREA * (u2 * Z2 / cc2**2 + rho2),
-AREA * (u2 * Z2 / cc2**2 - rho2),
-AREA * u2,
] # mass flux
# momentum
M[3] = [
Z1 + u1**2 * Z1 / cc1**2 + 2 * rho1 * u1,
Z1 + u1**2 * Z1 / cc1**2 - 2 * rho1 * u1,
-(Z2 + u2**2 * Z2 / cc2**2 + 2 * rho2 * u2),
-(Z2 + u2**2 * Z2 / cc2**2 - 2 * rho2 * u2),
-(u2**2),
]
ht1, ht2, src = a1 * Z1 / cc1**2 + b1 * Z1, a2 * Z2 / cc2**2 + b2 * Z2, delta * F / u1
M[4] = [-ht1 - src, -ht1 + src, ht2, ht2, a2] # energy + heat-release source
return np.linalg.det(M)
def reacting_mode(mm, cal1, cal2, delta, n, tau, seed):
w = complex(seed)
for _ in range(200):
h = 1e-3 * (abs(w) + 1.0)
dp = (
reacting_det(w + h, mm, cal1, cal2, delta, n, tau) - reacting_det(w - h, mm, cal1, cal2, delta, n, tau)
) / (2 * h)
if dp == 0:
break
s = reacting_det(w, mm, cal1, cal2, delta, n, tau) / dp
w -= s
if abs(s) < 1e-9:
break
return w.real / (2 * np.pi), -w.imag
def cp_eff(ed): # gamma R / (gamma - 1) from the mean state (sound-speed consistent)
g = ed["rho"] * ed["c"] ** 2 / ed["p"]
return g * (ed["p"] / (ed["rho"] * ed["T"])) / (g - 1.0)
# the same n-tau response, now on the reacting flame (a larger lag, since the equilibrium
# flame's intrinsic quasi-steady response must be overcome)
TAU_R = 6.0e-3
sol_r = rijke_reacting(N, TAU_R)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
res_r = sol_r.eigenmodes(freq_band=(50.0, 200.0), growth_band=(-250.0, 250.0), isentropic=True)
for mm in res_r.summary():
flag = " <-- UNSTABLE" if mm["unstable"] else ""
print(f"f = {mm['freq_hz']:7.2f} Hz growth = {mm['growth_rate']:+8.2f} 1/s{flag}")
# quantitative check against the equilibrium-flame dispersion relation
e1, e2 = sol_r.edge(1), sol_r.edge(2) # edge 1 unburnt approach, edge 2 burnt products
means = (e1["rho"], e1["c"], e1["u"], e1["p"], e2["rho"], e2["c"], e2["u"], e2["p"])
delta = 0.5 * (cp_eff(e1) + cp_eff(e2)) * (e2["T"] - e1["T"])
cal1, cal2 = reacting_caloric(sol_r, 1), reacting_caloric(sol_r, 2)
f_fns, g_fns = max(zip(res_r.freqs, res_r.growth_rates), key=lambda t: t[1])
f_an, g_an = reacting_mode(means, cal1, cal2, delta, N, TAU_R, 2 * np.pi * f_fns)
print(f"\nmost-unstable mode: Nefes f = {f_fns:6.2f} Hz, growth = {g_fns:+7.2f} 1/s")
print(f" analytic f = {f_an:6.2f} Hz, growth = {g_an:+7.2f} 1/s")
res_r.plot_spectrum(title=f"Reacting Rijke spectrum (H2/air, n={N}, tau={TAU_R*1e3:.0f} ms)").show()f = 110.99 Hz growth = +56.83 1/s <-- UNSTABLE
most-unstable mode: Nefes f = 110.99 Hz, growth = +56.83 1/s
analytic f = 110.99 Hz, growth = +56.83 1/seigenmodes then finds the self-excited modes.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(), "rijke_tube.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/tmp_nivlibu/rijke_tube.yaml