'nefes'How the source stamps into the operator
1 Dynamic sources
Of all the blocks that build the acoustic operator, one carries physics that is prescribed rather than derived from the mean-flow residuals: the unsteady feedback an element injects into the acoustics, of which a flame’s heat-release response is the archetype. This is the source block \mathbf{S}(\omega) of the perturbation network, and it is what closes the loop of a thermoacoustic instability — the acoustics perturb the flame, and the flame’s modulated heat release perturbs the acoustics in return. Whether that loop grows or decays is the stability question of analyses; this document supplies the source models it feeds on.
A dynamic source is specified, not solved: the analyst provides a transfer function describing how a heat-release (or mass-injection) fluctuation follows a reference fluctuation in the flow.
1.1 The general form of a dynamic source
A dynamic source is the statement that a source quantity carried by one element — the heat release of a flame, or the injected mass flow of a fuel port — fluctuates in response to the unsteady flow elsewhere in the network. Neither the quantity that drives the response nor the location at which it is read is fixed by the construction: the response is a superposition of terms, each reading its own reference quantity on its own reference edge, and is given as:
\frac{\widehat{\omega}}{\overline{\omega}} \;=\; \sum_{k=1}^{K} \mathcal{F}_k(\omega)\;\frac{\widehat{\phi}_k\big|_{e_k}}{\overline{\phi}_k\big|_{e_k}},
where \widehat{q} is the complex amplitude of the modulated source quantity and \overline{q} its mean, K is the number of terms, \mathcal{F}_k is the complex transfer function of term k evaluated at the angular frequency \omega = 2\pi f, and \widehat{\phi}_k/\overline{\phi}_k is the fractional fluctuation of its reference quantity \phi_k, evaluated on its reference edge e_k and normalized by the mean of that quantity on that same edge. The modulated quantity \omega is the unsteady heat release \dot{Q} [W] for a flame and the unsteady injected mass flow \dot{m}_{\text{src}} [kg/s] for a mass source; the reference quantity \phi_k may be a velocity, a static pressure, a density, a mass flow, or any transported composition scalar. Each term can be interpreted as a sensor placed on edge e_k that reads the fractional fluctuation of \phi_k there, and \mathcal{F}_k(\omega) as the frequency response through which the source follows that sensor.
An important remark is that the reference edge e_k is an arbitrary edge of the network. It need not be adjacent to the element carrying the source, it need not lie upstream of it, and the K terms need not share an edge: a source may respond simultaneously to the velocity on one edge, the pressure on another, and a composition scalar on a third, each through a transfer function of its own (test: test_descriptor_analytic_and_max_delay, two terms on two distinct edges). This is what makes the source block \mathbf{S}(\omega) structurally a feedback across the graph rather than a local element property: term k places its coefficients in the columns of edge e_k, on the residual row of the element that owns the source, so the two may be arbitrarily far apart in the network.
A less obvious but equally important point concerns what the summation presumes about its terms. The terms are superposed as though each carried a distinct physical sensitivity, yet the construction places no constraint on the references: nothing prevents two terms from reading quantities that the acoustics themselves relate, and the assembled operator is well defined either way. The responsibility for keeping the terms non-redundant therefore rests with the analyst, and two consequences deserve to be stated plainly.
First, two terms that share both a reference edge and a reference quantity are exactly degenerate: their coefficients coincide, so the operator receives only the sum \mathcal{F}_1 + \mathcal{F}_2 of their transfer functions and cannot distinguish them from a single term carrying that sum. Second, and more generally, the individual \mathcal{F}_k acquire a separate physical meaning only insofar as their reference fluctuations are dynamically independent. Where two references are constrained to co-vary — the velocity and the mass flow on the same edge, say, which the linearized state relations tie together — the network response determines only their combined effect, and the split between the two terms is a modeling assumption rather than a consequence of the physics. The practical hazard is double counting: a gain calibrated from data in which a second reference quantity was free to co-vary already contains that second sensitivity, and adding an explicit term for it counts the same physics twice.
This distinction is exactly what makes the terms recoverable or not. In the forward problem the degeneracy is silent, since prescribed \mathcal{F}_k always assemble into a definite operator; in the inverse problem it is not, where separating several sensitivities from one measurement requires the excitation to render their reference fluctuations sufficiently independent, and the resulting condition number is the diagnostic that guards it (see identification).
Working with fractional rather than absolute fluctuations is what renders each \mathcal{F}_k dimensionless, and with it the model portable between operating points. The normalization carries one further restriction: the mean \overline{\phi}_k on the reference edge must be nonzero, since the fractional fluctuation is otherwise undefined, and a source referencing a quiescent edge is rejected at assembly rather than silently producing an infinite gain. In practice most flames are represented adequately by the single velocity term of the next section, and the form above is the general case that admits, for example, a separate pressure sensitivity or an equivalence-ratio term read at the injector.
1.2 The flame transfer function
The archetype of the general form, and by a wide margin its most common use, is the flame transfer function (Lieuwen 2012): a single term, K = 1, whose reference quantity is the velocity and whose reference edge is the one immediately upstream of the flame. Under this specialization the heat-release fluctuation is that upstream velocity fluctuation, scaled by a frequency response, given as:
\frac{\widehat{q}}{\overline{q}} \;=\; \mathcal{F}(\omega)\;\frac{\widehat{u}}{\overline{u}},
where \widehat{u}/\overline{u} is the fractional velocity fluctuation on the upstream edge, and \mathcal{F} is the single transfer function that remains once the term index is dropped. It should be noted that the restriction to a single upstream velocity term is a modeling choice about the flame, not a limitation of the framework; it is adopted here because the closures of the n–\tau section below are stated for it.
1.3 The mass-source response
The second quantity a dynamic source may modulate is the mass flow delivered by an inline injector, an element that adds a stream of prescribed total temperature and composition to the through-flow without reacting it. Setting q = \dot{m}_{\text{src}}, the injected mass flow, the general form reads:
\frac{\widehat{\dot{m}}_{\text{src}}}{\overline{\dot{m}}_{\text{src}}} \;=\; \sum_{k=1}^{K} \mathcal{F}_k(\omega)\;\frac{\widehat{\phi}_k\big|_{e_k}}{\overline{\phi}_k\big|_{e_k}},
where \widehat{\dot{m}}_{\text{src}} is the complex amplitude of the injected mass flow and \overline{\dot{m}}_{\text{src}} the mean the element delivers, the remaining symbols carrying the meanings assigned above. A fuel feed whose delivery responds to the pressure fluctuation across its injection holes is the physically natural case, taking \phi = p on the feed-side edge; a velocity-modulated feed is written just as readily, and is what the shipped builder assumes by default.
The consequences of a modulated injector are richer than those of a modulated flame, which perturbs the acoustics through a single energy balance. An injected mass pulse perturbs three balances at once: it adds mass, it adds the axial momentum \overline{\dot{m}}_{\text{src}}\,u_{\text{inj}} that the injected stream carries, and it modulates the mass-weighted mixing of the injected total enthalpy and composition into the outflow (the row placement is given in the stamping section; test: test_mass_source_feedback_on_node_rows). The momentum contribution is proportional to the injection velocity u_{\text{inj}} and therefore vanishes for transverse injection, u_{\text{inj}} = 0, which adds mass with no axial momentum. The mixing contribution is proportional to the contrast between the injected stream and the local mixture, \phi_{\text{src}} - \overline{\phi}_{\text{out}}, so an injector delivering a stream that already matches the outflow composition and enthalpy modulates neither, however strongly its mass flow fluctuates.
Intuitively, a fluctuating injector converts an acoustic fluctuation into convected waves: an equivalence-ratio wave and an enthalpy (entropy) spot, both carried downstream at the mean flow speed rather than at the speed of sound. These are the forward path to the fuel-flow-driven instabilities in which a downstream flame burns the modulated mixture and closes the loop through its own heat release, and, when such a spot reaches a compact nozzle, the origin of indirect combustion noise. We emphasize that the injector itself performs no reaction: it sets the mixture that a downstream flame element subsequently burns.
1.4 The n–\tau model and its roll-off
The canonical closure is the interaction-index / time-lag model, due to Crocco, who introduced it as the sensitive time-lag theory of combustion instability in liquid-propellant rocket motors (Crocco and Cheng 1956), and in which the flame responds with a fixed gain after a fixed delay, given as:
\mathcal{F}(\omega) \;=\; n\,e^{-\mathrm{i}\omega\tau}, \qquad \omega = 2\pi f,
where n is the interaction index (the gain) and \tau is the flame time lag, with f the frequency in hertz. Under the e^{+\mathrm{i}\omega t} convention of the acoustics, the factor e^{-\mathrm{i}\omega\tau} is the causal delay of the response behind the driving fluctuation, so the model has constant gain n and a phase that falls linearly with frequency at slope -\tau; both features are visible in Figure 1 (tests: test_ntau_value_and_phase, test_ntau_complex_frequency_is_analytic).
A pure n–\tau gain that never rolls off makes every high frequency equally excitable, which is unphysical and complicates the stability count; a first-order low-pass variant, the shape a kinematic model of a conical flame predicts (Fleifil et al. 1996; Dowling 1997), bounds the unstable band, given as:
\mathcal{F}(\omega) \;=\; \frac{n\,e^{-\mathrm{i}\omega\tau}}{1 + \mathrm{i} f/f_c},
where f_c is the cutoff frequency, and its pole sits at f = \mathrm{i} f_c in the upper half-plane — the stable side under the sign convention — so the roll-off adds damping at high frequency without introducing an instability of its own.
The first-order roll-off is monotone in gain, which suits a conical flame but not a V-shaped or swirl-stabilized one, whose measured response overshoots below the cutoff before falling off more steeply above it (Palies et al. 2011). A second-order roll-off reproduces that shape (Dowling 1997), given as:
\mathcal{F}(\omega) \;=\; \frac{n\,e^{-\mathrm{i}\omega\tau}}{1 - (f/f_c)^2 + 2\mathrm{i}\zeta\,f/f_c},
where \zeta is the damping ratio: below 1/\sqrt{2} the gain peaks near f_c, and the phase swings through a further half turn on top of the pure lag. Both poles again lie in the upper half-plane, at f = f_c\big(\mathrm{i}\zeta \pm \sqrt{1 - \zeta^2}\big), so the model stays analytically continuable for the eigenproblem as long as the search region does not reach up to the nearer of them (tests: test_lowpass2_matches_the_delayed_second_order_filter, test_lowpass2_is_analytic_below_its_poles, test_lowpass2_overshoots_where_the_first_order_form_cannot). This is the form OSCILOS prescribes for the flame of the EM2C combustor benchmark (Li et al. 2017).
Code
f = np.linspace(1.0, 500.0, 400)
F1 = np.array([complex(n_tau(1.0, 3.0e-3)(x)) for x in f])
F2 = np.array([complex(n_tau_lowpass(1.0, 3.0e-3, 150.0)(x)) for x in f])
fig = make_subplots(rows=2, cols=1, shared_xaxes=True,
subplot_titles=("gain |F|", "phase ∠F [deg]"))
fig.add_trace(go.Scatter(x=f, y=np.abs(F1), mode="lines", line=dict(color=COLORWAY[0]),
name="n–τ"), row=1, col=1)
fig.add_trace(go.Scatter(x=f, y=np.abs(F2), mode="lines", line=dict(color=COLORWAY[1], dash="dash"),
name="n–τ low-pass"), row=1, col=1)
fig.add_trace(go.Scatter(x=f, y=np.degrees(np.unwrap(np.angle(F1))), mode="lines",
line=dict(color=COLORWAY[0]), showlegend=False), row=2, col=1)
fig.add_trace(go.Scatter(x=f, y=np.degrees(np.unwrap(np.angle(F2))), mode="lines",
line=dict(color=COLORWAY[1], dash="dash"), showlegend=False), row=2, col=1)
fig.update_xaxes(title_text="frequency f [Hz]", row=2, col=1)
fig.update_yaxes(range=[0.0, 1.15], row=1, col=1)
figThe source enters the operator by adding to the rows the base Jacobian already populates, never overwriting them — so the element keeps its mean-flow relation and gains an unsteady term on top. For a flame the term lands on the downstream edge’s total-enthalpy transport row with a factor -\delta, given as the residual contribution:
\widehat{h}_{t} - \widehat{H}_{\text{donor}} - \frac{\widehat{q}}{\dot m} = 0, \qquad \delta = \frac{\overline{Q}}{\dot m},
where \widehat{H}_{\text{donor}} is the inherited donor enthalpy (see transport), \widehat{q} the heat-release fluctuation of the transfer-function form above, \dot m the through-flow, and \delta the mean specific enthalpy rise across the flame. An important consequence of the \omega = 0 value \mathcal{F}(0) = n \neq 0 is that the source block has a nonzero direct-current gain, so \mathbf{S}(0) \neq 0; this does not disturb the mean flow, because the steady solve ignores the source entirely and treats the flame as passive, but it means the acoustic operator at zero frequency is \overline{\mathbf{J}} + \mathbf{S}(0) rather than \overline{\mathbf{J}} alone (see perturbation network). A mass-injection source stamps analogously but onto the node’s mass and momentum rows and onto the scalar-mixing rows, modulating the enthalpy and composition waves a fuel pulse drags with it — the forward path to fuel-flow-driven instabilities (tests: test_heat_release_lands_on_downstream_energy_row_with_correct_sign, test_q_mean_override_scales_the_coupling).
1.5 Tabulated and continued transfer functions
A measured or simulated flame response arrives not as a formula but as a table of complex samples \mathcal{F}(f) on the real-frequency axis, and using such a table correctly requires care about where in the complex plane it will be evaluated. As a raw table it is defined only on the real axis, interpolating magnitude and phase between its samples, and it deliberately refuses a complex-frequency argument. This is sufficient for a forced response and for the real-axis Nyquist stability driver, both of which evaluate only on the real axis, but it is not sufficient for the contour eigensolver of analyses, which searches the complex plane for the roots of \det\mathbf{A}(\omega) and therefore needs an analytic function (tests: test_tabulated_recovers_samples_and_rejects_complex, test_stability_rejects_nonanalytic_transfer).
The bridge is analytic continuation, and its correct form follows from the physics of the response rather than from a fitting preference. A table of numbers does not determine its own extension off the real axis; one must first say what kind of function the response is, and the two physically meaningful answers give the two continuations below. If the response to a brief disturbance dies out after a finite time — true of a flame, which forgets a velocity disturbance after a finite convective time, and of any compact element without an internal resonator — the transfer function is a finite sum of pure delays, entire in the complex plane, and the impulse-response fit of the next subsection is the faithful continuation; this is the recommended default. If instead the response rings — a cavity damper, a resonant end plate — its transfer function genuinely has poles, and the rational fit further below is the honest description.
1.6 The impulse-response model
System identification from a broadband simulation or experiment delivers a flame response most naturally as the discrete impulse response h_j, the heat-release response to a unit velocity impulse j\,\Delta t seconds earlier (Polifke 2020). Its frequency response is the finite sum
\mathcal{F}(f) \;=\; \sum_{j} h_j\, e^{-\mathrm{i}\,2\pi f\, j\,\Delta t}, \tag{1}
whose zero-frequency gain is \sum_j h_j and whose longest lag is J\,\Delta t. A finite sum of exponentials is entire: it has no poles anywhere in the complex-frequency plane, so it continues into the eigensolver’s search region exactly, at any growth rate, with no fitted poles to place and none to avoid. The n–\tau model is its one-coefficient limit (tests: test_fir_single_spike_is_an_n_tau, test_fir_is_entire_and_matches_its_definition_off_the_real_axis).
When the data arrive as frequency samples rather than as an impulse response, the coefficients h_j are recovered by fit_impulse_response: a least squares on the samples with a penalty on the second difference of h_j, so the fit follows the trend of the data rather than every measurement wiggle, over a memory length the user states from the physics (a few times the largest transport delay, readable as the slope of the phase). By default the sample spacing is set from the top of the tabulated band, \Delta t = 1/(2 f_{\max}), so the fit carries no frequency content the data cannot constrain (tests: test_fit_impulse_response_recovers_a_known_response_exactly, test_fit_impulse_response_continues_off_the_real_axis_exactly, test_fit_impulse_response_smooths_noisy_data).
Because the response is sampled, Equation 1 repeats with period 1/\Delta t, and only frequencies below the Nyquist limit 1/(2\Delta t) are resolved (test: test_fir_is_periodic_at_the_sampling_rate). The BRS combustor benchmark uses this form: its flame response is published only as a figure, and reconstructing the impulse response behind that figure, rather than fitting a rational function to it, is what makes the intrinsic mode of the rig computable (Emmert et al. 2017).
The same reasoning extends to measured transfer and scattering matrices. The entries of a scattering matrix are causal responses of outgoing waves to incoming ones, so for a compact element they are finite-memory and continue entry by entry with the impulse-response fit. The entries of a transfer matrix are an algebraic rearrangement of those responses and do not share the property: they mix delays of both signs (a plain duct already has \cosh/\sinh entries that grow off the real axis in both directions), and the rearrangement divides by the transmission response, creating genuine poles. Measured matrix data should therefore be converted to scattering form, continued there, and converted back — the conversions preserve analyticity (test: test_impulse_continuation_of_a_finite_memory_scattering_matrix_is_exact).
1.7 The rational alternative
For a response with a genuine resonance, the barycentric rational fit given by the AAA algorithm (Nakatsukasa et al. 2018) represents the samples as a rational function — analytic everywhere except at its isolated poles — so that the same object serves both the real-axis sweep and the complex-plane search. Two refinements make the fit trustworthy near the real axis: the pure transport lag is estimated and peeled before fitting, since a delay e^{-\mathrm{i}\omega\tau} is entire and carries no poles, and re-applied analytically on evaluation — the descriptor form of an analytic delay times a low-order rational remainder — and spurious Froissart-doublet poles introduced by the fit are detected and removed. The resulting continuation exposes its poles and zeros for diagnosis, and warns when a pole falls inside the frequency–growth window the stability search will scan (tests: test_descriptor_analytic_and_max_delay, test_fit_is_analytic_and_feeds_dynamic_sources).
That warning is not a formality, and it is why the rational fit is the second choice rather than the default. A rational fit of samples taken on the real axis places its poles near the sampled interval, because that is where the data constrains it; a mode sitting ten or twenty hertz off the real axis is then read out of the region where the fit is least trustworthy, and a fit of noisy data can scatter artificial poles into the search window outright. A finite-memory response has no such hazard — its continuation has no poles to place — which is why the impulse-response route is preferred whenever the physics permits it.
1.8 The reacting-acoustics caveat
One modeling gap must be kept in view when the flow carries composition, and stating it precisely is the point of this section. The coupling by which a composition fluctuation generates sound — the compositional or indirect noise (Magri et al. 2016), captured by a coefficient R_\xi — is retained everywhere the acoustic linearization is inherited from the mean-flow kernel: at a flame, an area change, a resolved nozzle, and even the compact choked-nozzle element, whose critical-mass-flux row is complex-stepped through its composition dependence and so carries R_\xi automatically. It is dropped in exactly one place: the hand-written analytic terminal closures for a choked nozzle or a constant-mass-flow outlet, which overwrite the terminal row with a three-wave (f, g, h) relation that retains the entropy coupling R_s but has no composition column. Accordingly the solver raises a compositional-noise warning precisely when reacting scalars are present and the flow is terminated by one of those analytic closures, so the gap is surfaced rather than silent (test: test_inherited_nozzle_carries_compositional_noise_analytic_closure_drops_it). One further approximation of the same low-Mach character is worth recording, and it enters only when the mean heat release \overline{Q} is left to be inferred rather than supplied. The perfect-gas heat-release flame carries its power as an element parameter and is de-normalized exactly; every other flame falls back on the sensible-enthalpy rise \dot m\,\overline{c}_p\,\Delta T, which neglects the kinetic-energy difference across the flame and is therefore accurate to \mathcal{O}(M^2). Supplying the mean heat release explicitly removes the approximation altogether.
The network’s stability is a computable property of the operator — whether the flame–acoustic loop grows (see analyses). When the flame’s response is itself the unknown, it can be recovered from a measured network response by the identification procedure.
References
Crocco, Luigi, and Sin-I Cheng. 1956. Theory of Combustion Instability in Liquid Propellant Rocket Motors. AGARDograph No. 8. Butterworths Scientific Publications.
Dowling, Ann P. 1997. “Nonlinear Self-Excited Oscillations of a Ducted Flame.” Journal of Fluid Mechanics 346: 271–90. https://doi.org/10.1017/S0022112097006484.
Emmert, Thomas, Sebastian Bomberg, Stefan Jaensch, and Wolfgang Polifke. 2017. “Acoustic and Intrinsic Thermoacoustic Modes of a Premixed Combustor.” Proceedings of the Combustion Institute 36 (3): 3835–42. https://doi.org/10.1016/j.proci.2016.08.002.
Fleifil, M., A. M. Annaswamy, Z. A. Ghoneim, and A. F. Ghoniem. 1996. “Response of a Laminar Premixed Flame to Flow Oscillations: A Kinematic Model and Thermoacoustic Instability Results.” Combustion and Flame 106 (4): 487–510. https://doi.org/10.1016/0010-2180(96)00049-1.
Li, Jingxuan, Dong Yang, Charles Luzzato, and Aimee S. Morgans. 2017. Open Source Combustion Instability Low Order Simulator (OSCILOS–Long) Technical Report. Department of Mechanical Engineering, Imperial College London. https://www.oscilos.com.
Lieuwen, Tim C. 2012. Unsteady Combustor Physics. Cambridge University Press.
Magri, Luca, Jeffrey O’Brien, and Matthias Ihme. 2016. “Compositional Inhomogeneities as a Source of Indirect Combustion Noise.” Journal of Fluid Mechanics 799: R4.
Nakatsukasa, Yuji, Olivier Sète, and Lloyd N. Trefethen. 2018. “The AAA Algorithm for Rational Approximation.” SIAM Journal on Scientific Computing 40 (3): A1494–522.
Palies, Paul, Daniel Durox, Thierry Schuller, and Sébastien Candel. 2011. “Nonlinear Combustion Instability Analysis Based on the Flame Describing Function Applied to Turbulent Premixed Swirling Flames.” Combustion and Flame 158 (10): 1980–91. https://doi.org/10.1016/j.combustflame.2011.02.012.
Polifke, Wolfgang. 2020. “Modeling and Analysis of Premixed Flame Dynamics by Means of Distributed Time Delays.” Progress in Energy and Combustion Science 79: 100845. https://doi.org/10.1016/j.pecs.2020.100845.