engine.stability_hifi.validation.v3_ntau_duct Namespace Reference

Classes
class	V3Result

Functions
complex	dispersion_G (complex lam, *, float c, float L, float x_f, float y_noz, float beta, float tau)
complex	dispersion_root (complex lam_seed, *, float c, float L, float x_f, float y_noz, float beta, float tau, int iters=60, float max_step=20.0)
V3Result	run_case (*, float c, float L, float R, float x_f, float gamma, float Mbar_e, float flame_gain, float tau, int nx1, int nx2, int nr, bool rigid_end=False)

Detailed Description

Verification case V3: duct with a compact n-tau flame, closed/choked ends
(paper Section VII.A: "active flame term, delay nonlinearity, complex lambda,
all three solvers agree").

This is the first case where the eigenproblem is genuinely NONLINEAR in lambda
(the flame delay exp(-lambda*tau), Eq. 7a/11), so it is the first exercise of the
paper's Section IV solver hierarchy: the frozen-delay fixed point (Algorithm 1) and
the bordered Newton polish (Section IV.C). (The third solver of the hierarchy, Beyn
contour integration, is a completeness AUDIT rather than a mode solver and is deferred
to P3 per the paper's phasing — so "all three solvers" is, at P0, these two. Flagged
rather than silently reinterpreted.)

Configuration
-------------
A duct of length L, uniform sound speed c, rigid ("closed") end at x=0, compact
choked-nozzle admittance y_noz = (gamma-1)*Mbar_e/2 (Appendix C) at x=L, and a compact
(delta-in-x, uniform-in-r) pressure-coupled flame at x = x_f:

    q_hat(x) = gain * exp(-lambda*tau) * p_hat(x_f) * delta(x - x_f)

i.e. Eq. (7a) lumped to a single reference point and a single flame sheet, with
``gain`` absorbing (gamma-1) * n_p * (qbar_dot/pbar) * (flame axial thickness)
[units m/s]. The flame is uniform across the cross-section, so for the r-independent
longitudinal modes compared here the 2-D axisymmetric FEM and the 1-D analytic model
below describe *identical* physics.

Discrete form solved (Eq. 11, rank-1 flame):

    N(lambda) p = [K + lambda*C + lambda^2*M2 - lambda*gain*exp(-lambda*tau)*g b^T] p = 0

with b = point sample of p at (x_f, 0) and g = the r-weighted disk load at x_f
(two DIFFERENT vectors; see acoustics/assembly.py for why conflating them is a bug).

Analytic reference (derived, not from the paper)
------------------------------------------------
Zone solutions satisfying the end conditions, with s = lambda/c:

    zone 1 (0..x_f):   P(x) = cosh(s*x)                          [P'(0)=0, rigid]
    zone 2 (x_f..L):   Q(x) = r2*exp(s*(x-L)) + exp(-s*(x-L))

The admittance end condition p'(L) = -s*y_noz*p(L) (from Eq. 8 with y constant) gives

    r2 = (1 - y_noz) / (1 + y_noz)      [lambda-independent for a compact nozzle].

CAUTION (documented wrong turn): an earlier draft had r2 = (y-s')/(y+s') -> the
NEGATIVE of the correct value. The passive limit exposes it: with the wrong sign the
flame-off roots do not reproduce the independently verified tanh(s*L) = -y_noz
spectrum (sigma ~ -y*c/L, f ~ n*c/2L). Always check the passive limit of an
active-flame dispersion relation first.

At the flame, p is continuous and the slope jumps. The jump follows from the FEM's own
strong form (integrate lambda^2*p - d/dx(c^2 p') = lambda*gain*e^{-lambda*tau}
* p(x_f) * delta(x-x_f) across the sheet):

    [p']_{x_f-}^{x_f+} = -lambda * beta * exp(-lambda*tau) * p(x_f),
    beta = gain / c^2.

CAUTION (the bug that cost the most time): beta is gain/c^2, NOT gain*(R^2/2)/c^2.
The r-weighted disk load g sums to R^2/2, but the r dOmega measure multiplies every
OTHER matrix (K, M2) by the same R^2/2, so it cancels exactly in the 1-D reduction.
Double-counting it made the analytic flame coupling 1/(R^2/2) = 5000x too weak, which
presented as "FEM growth rate 5000x larger than analytic" — with frequencies agreeing
to 0.003%, the signature of a coupling-scale error rather than a discretization or
solver error. (The FEM answer was mesh-converged and confirmed independently by
first-order eigenvalue perturbation theory of the discrete system; the closed-form
sensitivity for the fundamental mode of the rigid-rigid case,
d(lambda)/d(gain) = exp(-lambda0*tau) * cos^2(pi*x_f/L) / L, in which R cancels
entirely, is re-checked in the test suite.)

Eliminating the two amplitudes with the two interface conditions gives the dispersion
relation, root-found in pole-free product form (V2's lesson):

    G(lambda) = P'(x_f) Q(x_f) - Q'(x_f) P(x_f)
                - lambda*beta*exp(-lambda*tau) * P(x_f) Q(x_f) = 0.

Phase convention worth knowing (physics finding, verified numerically three ways):
with the pure-delay coupling (7a), the Rayleigh driving of a mode goes as
+cos(omega*tau) — in-phase heat release drives, anti-phase damps — because the
cycle-averaged p'q' is |p_hat|^2 * cos(omega*tau): the REAL part of the transfer. The
lumped model (engine/pipeline/stability/acoustic.py) uses the DIFFERENCE form of the
Crocco response, n*[p'(t) - p'(t-tau)], whose transfer n*(1-e^{-i*omega*tau}) has real
part n*(1 - cos(omega*tau)) -- driving everywhere, most at omega*tau = pi. (It took the
imaginary part, n*sin, until 2026-09-26; that is the quadrature component and does no
work.) The two forms differ by the instantaneous n*p'(t) term, so their tau-bands are
not the same; comparisons between this framework and the lumped tier must translate
conventions first.

Function Documentation

dispersion_G()

complex engine.stability_hifi.validation.v3_ntau_duct.dispersion_G	(	complex	lam,
		*	,
		float	c,
		float	L,
		float	x_f,
		float	y_noz,
		float	beta,
		float	tau )

Pole-free dispersion residual G(lambda) (module docstring). Root <=> eigenvalue.

Definition at line 125 of file v3_ntau_duct.py.

Referenced by dispersion_root().

dispersion_root()

complex engine.stability_hifi.validation.v3_ntau_duct.dispersion_root	(	complex	lam_seed,
		*	,
		float	c,
		float	L,
		float	x_f,
		float	y_noz,
		float	beta,
		float	tau,
		int	iters = 60,
		float	max_step = 20.0 )

Bounded-step complex Newton on G, seeded near the expected root.

Deliberately LOCAL: we want the dispersion root nearest the seed (the FEM
eigenvalue being cross-checked), and an unbounded/global root-finder can wander to
a distant root and report success — scipy.fsolve did exactly that during
development, converging to an unrelated real root thousands of rad/s away.
Guaranteed-complete root surveys are the job of the Beyn audit tier (P3), not of a
verification cross-check. The derivative is numerical (central difference, step
h=1e-2 — G is smooth on O(1) lambda scales, so this is far below its variation).

Definition at line 138 of file v3_ntau_duct.py.

References dispersion_G().

Referenced by run_case().

run_case()

V3Result engine.stability_hifi.validation.v3_ntau_duct.run_case	(	*	,
		float	c,
		float	L,
		float	R,
		float	x_f,
		float	gamma,
		float	Mbar_e,
		float	flame_gain,
		float	tau,
		int	nx1,
		int	nx2,
		int	nr,
		bool	rigid_end = False )

Full V3 run: FEM NLEVP (both solvers) vs the 1-D dispersion relation.

``rigid_end=True`` replaces the choked-nozzle admittance at x=L with a rigid wall
(y_noz = 0) — the pure-flame sub-case, useful for isolating flame handling from
boundary-damping handling when something disagrees.

The mesh is built with ``two_zone_duct_mesh`` even though there is no material
jump, purely to guarantee an exact node column at x = x_f for the disk load
(see ``disk_load_vector``'s docstring).

Definition at line 164 of file v3_ntau_duct.py.

References dispersion_root().