engine.stability_hifi.validation.v2_temperature_jump Namespace Reference

Classes
class	ModeComparison

Functions
float	_dispersion_g (float f_hz, float c1, float L1, float c2, float L2)
List[float]	analytic_frequencies (float c1, float L1, float c2, float L2, float f_max_hz, int samples_per_period=40)
List[ModeComparison]	run_case (*, float c1, float L1, float c2, float L2, float R, int nx1, int nx2, int nr, int n_compare=4)

Detailed Description

Verification case V2: 1-D duct with a temperature (sound-speed) jump, passive
(paper Section VII.A: "checks: nonuniform-c-bar handling").

The paper *names* this case but does not give its analytic reference solution, so the
derivation below is original (not copied from the paper) — done directly from the
governing equation the paper does give, Eq. (6):

    lambda^2 p_hat - d/dx( c_bar^2 dp_hat/dx ) = 0      (1-D, no heat release)

Setup: a duct of length L = L1 + L2, split at x = L1 into two uniform zones with sound
speeds c1 (zone 1, 0 <= x <= L1) and c2 (zone 2, L1 <= x <= L1+L2). Rigid ("closed") ends
at x=0 and x=L, matching Section III.E's rigid-wall Neumann condition. Only the
sound-speed field is nonuniform; nothing about the boundary condition or the flame
(there is none) is new relative to V1 — this isolates exactly the one new piece of
physics V2 is meant to check.

Derivation
----------
Looking for undamped natural frequencies (lambda = i*omega, omega real — same as V1;
a lossless medium can't grow or decay), Eq. (6) in each uniform zone reduces to the
textbook 1-D Helmholtz equation

    d^2 p_hat/dx^2 = -k_i^2 p_hat,      k_i = omega / c_i     (i = 1, 2)

with general solution p_hat(x) = A*cos(k_i*x) + B*sin(k_i*x). The rigid-end conditions
p_hat'(0) = 0 and p_hat'(L) = 0 fix the form in each zone up to one free amplitude:

    zone 1:  p_hat_1(x) = A1 * cos(k1*x)                    [only cos satisfies p'(0)=0]
    zone 2:  p_hat_2(x) = A2 * cos(k2*(x - L))               [satisfies p'(L)=0 by construction]

At the interface x=L1, we need two matching conditions. Pressure continuity, (a), is
uncontroversial (nothing sources or sinks p there). Condition (b) is the one that is
easy to get wrong, and an earlier draft of this derivation *did* get it wrong (kept
below as a documented correction, since the mistake and the fix are both instructive):

    WRONG first attempt: continuity of mass flux rho_bar*u_hat, on the (mistaken)
    reasoning that this is the standard duct-acoustics matching condition when gas
    properties change across an interface. Using u_hat = -grad(p_hat)/(lambda*rho_bar)
    (given right after Eq. 6), rho1*u_hat_1 = rho2*u_hat_2 reduces (both rho AND lambda
    cancel) to dp_hat_1/dx = dp_hat_2/dx, giving k1*tan(k1*L1) + k2*tan(k2*L2) = 0.

    This numerically MISMATCHED the FEM solution by several percent at every mesh
    resolution, with no improvement on refinement — a strong signal (per the mesh-
    convergence-order test in V1) of a wrong reference formula, not a discretization
    error. Tracking it down: the standard "mass-flux continuity" rule from general duct
    acoustics assumes a mean flow physically carrying mass across the interface. But
    Section III.C's Helmholtz reduction explicitly assumes a *quiescent* mean flow
    (u_bar ~ 0) — there is no throughflow here at all, just a spatially-varying sound
    speed in still gas. The right way to find the correct condition is to go back one
    step further than Eq. 6, to the *pair* of first-order equations it was combined
    from (Eq. 4): the energy equation "dp'/dt + gamma*p_bar*div(u') = (gamma-1)*q_dot'"
    uses the volumetric dilatation div(u') directly (a purely kinematic quantity, not a
    mass flux), and Appendix B's assumption (ii) is that gamma*p_bar is spatially
    UNIFORM (mean pressure barely drops across the chamber) even though rho_bar is not.
    Integrating that energy equation across a vanishingly thin control volume straddling
    the interface forces the u' term's jump to vanish for the equation to stay finite —
    i.e. **u_hat itself (not rho_bar*u_hat) must be continuous.** Equivalently, since
    c_bar^2 = gamma*p_bar/rho_bar and gamma*p_bar is the same constant on both sides,
    continuity of u_hat (using u_hat = -grad(p_hat)/(lambda*rho_bar)) is the same
    statement as continuity of c_bar^2 * dp_hat/dx — which is also exactly the *natural*
    (weak-form) interface condition that Eq. (10)'s FEM discretization enforces
    automatically at any element boundary where c jumps. That the corrected physics
    argument and the FEM's own automatic behavior agree is a good consistency check.

Redoing the algebra with the corrected condition (b): c1^2 * dp_hat_1/dx = c2^2 *
dp_hat_2/dx at x=L1. Substituting the two zone solutions into "p continuous" and
"c^2 p' continuous" and eliminating A1, A2 (using k_i = omega/c_i, so c_i^2 * k_i =
c_i * omega):

    A1 cos(k1*L1) = A2 cos(k2*L2)                                (p continuous)
   -c1*omega*A1 sin(k1*L1) = c2*omega*A2 sin(k2*L2)               (c^2 p' continuous)

gives the dispersion relation:

    c1 * tan(k1*L1) + c2 * tan(k2*L2) = 0.

Sanity check (uniform limit): if c1 = c2 = c, both this and the wrong first attempt
collapse to the same tan(k*L1) = -tan(k*(L-L1)), i.e. f = n*c/(2L) — V1's closed-closed
spectrum. **This is why the uniform-limit check alone did not catch the error**: any
relation of the schematic form (function of c1, k1*L1) + (function of c2, k2*L2) = 0
that is antisymmetric the same way trivially passes it regardless of which power of
c_i actually belongs there. The uniform limit is a necessary but not sufficient check;
comparing against the independently-built FEM solution (which does not share whatever
mistake was made in the by-hand derivation) is what actually caught this one.

Numerically, ``tan`` has poles wherever k_i*L_i passes through (2n+1)*pi/2, which are
not physical roots but *look* like sign changes to a naive root-finder. To avoid that
trap, the code below root-finds the equivalent POLE-FREE form obtained by multiplying
through by cos(k1*L1)*cos(k2*L2):

    g(omega) = c1*sin(k1*L1)*cos(k2*L2) + c2*cos(k1*L1)*sin(k2*L2) = 0

which is smooth (entire) in omega, so every sign change of g really is a root.

Function Documentation

_dispersion_g()

float engine.stability_hifi.validation.v2_temperature_jump._dispersion_g ( float f_hz, float c1, float L1, float c2, float L2 )

protected

Pole-free dispersion residual (see module docstring). Root at f_hz => a true mode.

Corrected form ``c1*sin(k1L1)cos(k2L2) + c2*cos(k1L1)sin(k2L2)`` (velocity-continuity
interface condition) — NOT ``k1*sin(...)+k2*sin(...)`` (mass-flux continuity), which
was tried first and shown to be wrong; see the module docstring for the full story.

Definition at line 118 of file v2_temperature_jump.py.

Referenced by analytic_frequencies().

analytic_frequencies()

List[float] engine.stability_hifi.validation.v2_temperature_jump.analytic_frequencies	(	float	c1,
		float	L1,
		float	c2,
		float	L2,
		float	f_max_hz,
		int	samples_per_period = 40 )

All positive roots of the dispersion relation below ``f_max_hz``.

Dense-sample ``g`` finely enough to resolve its fastest oscillation (set by the
shorter, faster zone) and Brent's method to refine each sign change found. Because
``g`` (unlike the raw tan-tan form) has no poles, every sign change IS a root —
no risk of mistaking an asymptote for a mode, which is exactly why this pole-free
form was used instead of the more directly-derived tan-tan equation.

Definition at line 131 of file v2_temperature_jump.py.

References _dispersion_g(), max, and min.

Referenced by run_case().

run_case()

List[ModeComparison] engine.stability_hifi.validation.v2_temperature_jump.run_case	(	*	,
		float	c1,
		float	L1,
		float	c2,
		float	L2,
		float	R,
		int	nx1,
		int	nx2,
		int	nr,
		int	n_compare = 4 )

Solve the FEM two-zone duct at m=0 and match the lowest ``n_compare`` modes
against the dispersion-relation roots (same sort-and-pair strategy as V1; see the
completeness pitfall noted in ``docs/stability/stability_hifi_p0_v1_notes.md`` --
the analytic root search above already returns every root up to ``f_max_hz``, with
no truncation-by-branch that could silently skip one, so that specific bug class
does not apply here, but the general "did I search far enough" question still does,
hence ``f_max_hz`` below is set generously past the highest mode we intend to compare).

Definition at line 162 of file v2_temperature_jump.py.

References analytic_frequencies().