Forced flow response analysis of a turbulent swirling annular jet flame

This paper describes an investigation into the response of a reacting annular swirling jet flow to external acoustic forcing. Even in the absence of forcing, natural interactions among axial and azimuthal shear layers, centrifugal forces, and inertial waves endow swirling jets with complex physical underpinnings, which often manifest a remarkable array of distinct vortical structures. When combustion is present, such hydrodynamic features disturb and interact with the flame, often leading to significant heat release and pressure oscillations. This phenomenon, known as combustion instability, bears special relevance in engineering applications such as gas turbines, where it detrimentally impacts emissions, performance, and hardware integrity. In such systems, resonant acoustic modes associated with the combustor geometry can filter pressure fluctuations from the turbulent broadband, leading to a significant narrowband forcing effect and opening pathways for positive hydrodynamic and thermoacoustic feedback.^1,2 This motivates our current effort to understand how reacting turbulent swirling flows respond to narrowband forcing, and to devise low-order models which can predict their response.

The influence of external forcing toward a given flowfield is a strong function of the flow's natural dynamics and the strength and structure of the excitation. Fundamentally, shear flows tend to behave either as self-excited oscillators with intrinsic dynamics or as passive amplifiers which filter extrinsic dynamics. When clear, spatiotemporally coherent structures appear without imposed forcing, the turbulent system is self-excited and may be interpreted to lie in the vicinity of a low-dimensional attractor in phase space. This “noisy” attractor must be somewhat robust to perturbations, as turbulent fluctuations would otherwise carry the system away from it.³ Under such conditions, external excitation with an improper frequency, insufficient amplitude, or incorrect spatial pattern is unlikely to have a significant effect on the flow dynamics.⁴ However, if forcing with a sufficiently high amplitude, an appropriate temporal distribution, and a suitable spatial arrangement is applied, the system may bifurcate to a new set of dynamics governed by a different attractor. This is known as the “lock-in” phenomenon.^5,6

On the other hand, amplifier flows, which are the main subject of this paper, generally do not show such clear bifurcations. Naturally, these flows tend to exhibit more spectrally distributed spatiotemporal behavior. However, forcing, even with small amplitudes and low authority, can fundamentally alter the flow response by directly introducing and controlling the narrowband, spatiotemporally coherent structures present in the flow.⁴ In amplifier flows, the scale of the response typically depends directly on both the amplitude of the imposed forcing and the flow's underlying receptivity to the forcing frequency and structure.^6,7 This sensitivity implies a linear relationship between the imposed acoustic field and the dominant flow response, provided the forcing amplitude is sufficiently small to permit a linear approach.

Over the past few decades, linear hydrodynamic stability analysis has emerged as a prominent technique for understanding and modeling the leading-order coherent dynamics in turbulent shear flows. In such models, the nonlinear dynamics of the turbulent coherent structures are approximated by the linear dynamics of infinitesimal instability waves superposed upon the turbulent mean flow. While this assumption is quite nuanced, several theoretical studies^8–10 have rigorously outlined its implications and limitations, and examples of successful applications abound in the fluid mechanics literature.^11–21 

In the particular context of swirling flows, Manoharan et al.²⁰ generalized the weakly nonlinear analyses of Sipp and Lebedev⁸ to a turbulent swirling jet, providing evidence that mean flow linear stability analysis should yield meaningful predictions for self-excited oscillations. This explains several studies that have successfully applied mean flow linear stability analysis to study the coherent self-excited hydrodynamic behavior. For example, Oberleithner et al.¹⁵ used mean-flow stability analysis to explain the self-excited dynamics of a swirling turbulent water jet, demonstrating excellent agreement between the stability model and experimental observations. More recently, Mukherjee et al.²¹ reported similar close agreement for the self-excited behavior in a swirling jet burner configuration nearly identical to the one analyzed in this paper.

The main thrust of this paper is directed toward a mean flow stability analysis, which is used to construct a linear hydrodynamic transfer function between an input forcing and an output response of an amplifier-type turbulent reacting swirling jet. Setting aside the important aspect of combustion, our study may be interpreted as the amplifier-flow analogue of the study by Oberleithner et al.¹⁵ toward oscillator-type swirling jet flows. The remainder of this paper is organized as follows. Section II describes the flow configuration and introduces the modeling formulation and numerical methods. Section III characterizes the flow measurements and presents the mean flow about which the model is linearized. Section IV provides an overview of the stability analysis calculations and describes the mode shapes, which serve as the basis for the predicted response. Finally, the predictions are presented and compared with measurements in Sec. V before the conclusion in Sec. VI.

Before moving on to the modeling approach, a brief overview of the experimental facility is given for context. The experiments which form the foundation of this work have been reported previously in the works of Smith et al.;^7,22,23 however, the essential details will be repeated here. The data were obtained in the optically accessible combustor facility shown in Fig. 1, which features an in-line array of three burners. The combustor box is equipped with speakers to deliver acoustic excitation and has a high aspect ratio cross section to accommodate transverse planar acoustic waves. This design emulates an “unwrapped” version of the environment within annular gas turbine combustors where transverse acoustic waves can have a significant influence on the system behavior.²⁴ For the cases considered here, which are listed in Table I, the flow consists of a lean mixture of 500 K preheated methane and air injected from swirled annular nozzles with outer diameters of $D⋆=36.83$ mm and geometric swirl numbers of 0.6 at a nominal velocity of $U⋆=25$ m/s. Considering the equilibrium properties of the unburned reactants, we here define the Reynolds number $Re=ρ⋆U⋆D⋆/μ⋆=23 813$ and Péclet number $Pe=U⋆D⋆/α⋆=19 041$⁠, where $ρ⋆, μ⋆$⁠, and $α⋆$ refer to the dimensional density, molecular viscosity, and thermal diffusivity, respectively. Henceforth, unless indicated otherwise, all quantities reported in this paper are made dimensionless using the starred quantities listed above.

The experiments utilize high-speed planar stereoscopic particle image velocimetry (sPIV) and synchronized OH-planar laser-induced fluorescence (PLIF) measurements for flow and flame characterization. The specific test cases investigated here include one unforced single-nozzle flow case which forms the basis of the model and two transversely forced triple-nozzle flow cases which are used to evaluate the model performance. In these transversely excited cases, the lateral speakers were set to oscillate either in-phase (IP) or out-of-phase (OP) to generate standing transverse acoustic waves in the combustor box which correspond either to pressure antinodes or nodes, respectively, at the center of the central nozzle. To obtain three-dimensional flow information, measurements were performed on the central jet in its meridional plane as well as four separate axial planes at $x=[0.15,0.67,1.2,1.7]$ in the forced test cases. Representative composite meridional plane images from the sPIV and OH-PLIF cameras are shown in Fig. 2. As shown in the figure, a cylindrical coordinate system $x=(x,r,θ)$ is adopted with the origin centered on the dump plane at the central nozzle. Further experimental details and descriptions of post-processing methods can be found in the original works.^7,22,23

The mean-flow stability analysis approach used in this work is derived following the arguments outlined by Tammisola and Juniper¹⁷ for nonreacting flows and applied to a reacting flow by Manoharan et al.¹⁸ It is based on the triple decomposition proposed by Hussain and Reynolds,²⁵ which separates the fluctuating part of the usual mean/fluctuating Reynolds decomposition into coherent and incoherent parts, i.e.,

for any spatiotemporally varying quantity $ψ(x,t)$⁠. The steady part is defined by the weighted time-average,

where $w(x,t)$ is a weighting function. The coherent part is based on the ensemble-average of zero time-mean periodic features with frequency f at the same point of phase (i.e., the phase-average less the time-mean). The weighted phase-average is defined formally as

leading to the following definition for the coherent oscillations:

However, in practice, the constant sampling frequency of the diagnostics permits the use of a discrete Fourier transform to identify the coherent oscillations. Finally, the incoherent part $ψ´(x,t)$ consists of the leftover fluctuations unexplained by the sum of the mean and coherent parts. Here, as in Manoharan et al.,¹⁸ we adopt a Favré averaging approach, where the decompositions for the density and hydrodynamic pressure (ρ and p, respectively) are based on unweighted averages (i.e., w = 1) and the decompositions for the velocity and temperature (u and T, respectively) are based on mass averages (i.e., $w=ρ/ρ¯$ and $w=ρ/⟨ρ⟩$ for the time- and phase-averages, respectively).

To formulate the mean-flow stability analysis, the turbulent flow is assumed to be governed by the low-Mach number Navier–Stokes equations and all variables are expanded in terms of (1). The equations governing the coherent oscillations are identified by taking the difference of the phase-averaged and time-averaged equations and invoking the assumption that the incoherent fluctuations represent a zero-mean stochastic process. Then, the time- and phase-averaged effect of nonlinear interactions involving these incoherent oscillations is modeled using a turbulence model described below, while the nonlinear interactions among only coherent oscillations are assumed to have a negligible effect on the coherent oscillations themselves.

The coherent transport of momentum by incoherent motions is modeled by invoking the Boussinesq hypothesis and using the mean eddy viscosity $μ¯T(x)$⁠. Practically, this assumption amounts to the requirements that

where $k=12ρ(u−u¯)·(u−u¯)$ is the turbulent kinetic energy (TKE) and $S=∇u+(∇u)T−13(∇·u)I$ is the deviatoric stress tensor. Following Tammisola and Juniper,¹⁷ the TKE and eddy viscosity are assumed to be unaffected by the coherent oscillation such that $k̃=μ̃T=0$⁠. The value of the mean eddy viscosity is then determined by projecting the measured Reynolds stress tensor from the experiments onto the scalar eddy viscosity model of (5). This process results in an overdetermined system of equations, which is best satisfied in the least squares sense via the relationship,

where $ρu´u´¯$⁠, $k¯$⁠, and $S¯$ are evaluated directly from the measurements. Altogether, the resulting model for the coherent transport of momentum by incoherent fluctuations is

Similarly, the coherent thermal energy transport by incoherent motions is treated with a gradient diffusion model based on a constant turbulent Prandtl number assumption as in Manoharan et al.¹⁸ In this approach, the turbulent thermal transport is captured by the same behavior, resulting in the dimensionless turbulent diffusion coefficient, $D¯T=μ¯T$⁠.

Synthesizing the triple-decomposed low-Mach number Navier–Stokes equations and the ideal gas equation of state with the assumptions and turbulence model described above, the resulting linear system for the coherent oscillation is

where $f̃(x,t)$ represents a coherent external body force. Finally, under the assumption of axisymmetry, the normal modes expansion

may be introduced as an ansatz for the coherent oscillations, where $q=(ρ,u,p,T)T$ is the state vector, m is the integer azimuthal wavenumber, and ω is the complex frequency. Note that even though the experimental configuration is not strictly axisymmetric, prior work has demonstrated that the mean flow is highly axisymmetric within the region of interest,⁷ suggesting that the nonaxisymmetric Fourier components of the mean flow are negligible and that the above expansion is valid.

The inhomogeneous linear system (9) represents a boundary value problem, which describes the system's forced response, provided that the spatial distribution of the imposed forcing is known throughout the domain and that suitable boundary conditions can be identified. While such information is often easily accessible when analyzing computational datasets, it is not available in many experimental situations due to the practical limitations of optical diagnostics systems. For the experiments featured in this work, the measurement window is focused just above the nozzle to minimize laser reflections from the polished combustor walls and ensure that all of the essential mean flow gradients in the near field are well-resolved. However, this positioning also leaves near-wall velocity gradients unresolved and only provides velocity data in the immediate vicinity of the nozzle, impeding the application of no-slip and far-field boundary conditions necessary to determine the response modes in a “global” spatial framework.

Rather than pursuing a strictly global approach, the analysis performed in this paper leverages a so-called “weakly global” approach to analyze the response of the flow to acoustic forcing. This technique, which is more accommodating of the challenges associated with experimental datasets, is based on the WKBJ theory of Chomaz et al.²⁶ and constructs weakly global solutions for the natural hydrodynamic eigenmodes of the flow using a sequence of local parallel-flow stability analyses. In the unforced context of weakly global theory, these local analyses can be used to determine the natural “wavemaker” location and frequency in the flow, as shown by Juniper and Pier.²⁷ In the forced situation considered here, however, the imposed excitation represents the sole wavemaker for the hydrodynamic stability problem, provided the flow is not self-excited.⁴ At the small Mach numbers characterizing the flows considered in this paper, the dominant hydrodynamic effect from the acoustic forcing is strongly localized to the nozzle boundaries, where acoustic reflections induce small, coherent shear layer vortices that may be amplified or diminished as they convect along the flow. Thus, a leading-order prediction of the forced response is obtained by interpreting the input forcing as the wavemaker and projecting this forcing onto the natural weakly global modes at the forcing frequency along the most upstream axial station resolved by the measurements. Since it only requires local stability analyses, this method has the additional benefit of not requiring streamwise boundary conditions, which could represent a significant source of error in a fully global (i.e., non-local) framework due to the limited spatial extent of the measurements.

Before moving on, it is worth briefly mentioning that the framework outlined above hinges on the assumption that the flow is not self-excited. If the self-excited (i.e., globally unstable) behavior were to be detected, the weakly global analysis would become significantly more delicate as the upstream forcing would only represent one of the possible wavemaker regions.^4,27 However, as will be shown in the results, no self-excited behavior was detected for the flow considered in this paper.

To proceed with formulating the local parallel-flow stability problems at each axial station, the WKBJ ansatz is introduced such that the global modes of (10) can be expressed in terms of local modes as

where k is the slowly varying axial wavenumber, x₀ is the position of the wavemaker, and $q̂Lm$ are the local modes where x appears as a parameter rather than an argument. As pointed out by Wang and Rusak,²⁸ this ansatz is also associated with an assumption of axial homogeneity, which prevents inertial waves propagating along the central axis from accumulating at streamwise boundaries. However, in the present approach, this corollary assumption is desirable since it further minimizes the influence of the uncertain streamwise boundary conditions on the hydrodynamic behavior.

Next, the mean flow is assumed to vary slowly in the axial direction (in the same manner as k) such that mean flow streamwise derivatives and radial velocities are negligible, and x again appears as a parameter rather than an argument. With these assumptions incorporated, the local stability problems which appear at each x can be deduced from the linearized system of (9) as

Here, the linear operator contained in (9) after the application of (10) has been decomposed into $[A0(r;m)+kA1(r;m)+k2A2(r;m)]$⁠. In (12), $A0$ is an operator which represents the part of $Jm$ which is constant with respect to k, $A1$ represents the part which is linear, and so on. Each of these operators is defined explicitly in the first author's doctoral dissertation.²⁹ Equation (12) is a local dispersion relationship which governs the hydrodynamic flow behavior at each x. Therefore, in the classical stability analysis framework, the local temporal stability is deduced by solving (12) as a linear eigenvalue problem for complex ω given real k, the local spatiotemporal stability is determined by solving (12) for the stationary points where $∂ω/∂k=0$ with ω and k complex, and local spatial stability is determined by solving (12) as a quadratic eigenvalue problem for k with an imposed real ω.^30,31

Determination of each weakly global mode is achieved by solving (12) as a spatial eigenvalue problem for a single real frequency ω_f (the external excitation frequency) over a sequence of axial positions to obtain a sequence of $q̂Lm$ and k. Due to the quadratic nature of the spatial eigenvalue problem, each local spatial solution has two branches of solutions which propagate in opposite directions along x. Here, as the wavemaker is considered to represent the upstream forcing at the burner lips where the hydrodynamics and acoustics are coupled through boundary effects, the weakly global modes are constructed from the downstream-propagating k⁺ solution branch using (11). Hence, these modes physically describe the asymptotic response to a generic pointwise harmonic forcing located at x₀. To ensure consistency in amplitude and phase across all axial stations, the magnitude of each local mode is regularized to a unit norm such that

at each x and its phase is rotated such that the norm of the difference between the local mode at x = x_i and the mode at the previous station $x=xi−1$ is minimized. This ensures that the phase and amplitude of the local modes remain constant with x, and the overall phase and amplitude evolution arise only due to (11).

Finally, the coherent response is constructed as a linear combination of the spatially amplified weakly global modes,

where $âj$ is a constant coefficient). As will be further explained in Sec. V, the linear coefficients underlying the sum are determined by projecting the input forcing onto the basis of weakly global modes at x₀ using the inner product contained in (13). Therefore, the resulting $q̃$ yields a prediction of the response to a given forcing input based on weighted contributions from the various weakly global modes.

The weakly global hydrodynamic stability problem described above is solved using a Chebyshev spectral–Galerkin method,^32,33 a technique closely related to the pseudospectral collocation methods commonly used by earlier local hydrodynamic stability studies of swirling flows.^15,34–36 As will be discussed below, the problems are considered in an unbounded radial domain to avoid confinement effects.³⁷ To achieve this, points in the physical space r are mapped to a Chebyshev–Gauss quadrature in the computational space $ξ=[−1,1]$ using the algebraic mapping function,

where α controls grid stretching and β controls grid clustering. This mapping function is equivalent to that used by Fabre and Jacquin³⁸ when β = 0; the inclusion of $β>0$ allows grid points to be clustered at a finite r. This, in turn, allows the annular jet profile to be resolved with fewer grid points than would be required if the points were concentrated about only the centerline. In this work, $α=1/4$ and $β=3/4$ are selected such that at least 50% of the quadrature points are positioned on the interval $r=[0.1,1]$⁠, where the mean flow variation is strongest.

To pose the local stability problems, boundary conditions must be enforced along r. However, an approximation is needed to enforce these conditions as the experimental measurement window does not extend to the actual physical boundaries of the flow apparatus. Conveniently, near the outer edge of the measurement plane, it was observed that the mean flow is always comprised of nearly stagnant unburned fluid. Therefore, the mean flow velocities at radial locations beyond the measurement window are assumed to be zero and the associated fluid properties are chosen to be those of the unburned fluid. To avoid generating artificial boundary effects, this interval of stagnant fluid is then extended to infinity, where the Galerkin expansion ensures that the density and velocity perturbations vanish regularly in the $r→∞$ limit. Note that no such requirements are enforced on the hydrodynamic pressure in order to yield an inf–sup compatible formulation. To enforce the proper constraints on the axis, symmetry conditions are enforced at $r=ξ=0$⁠. These conditions are determined by the even or oddness of m and exploit the parity of the Chebyshev polynomials and each component of $q̂$⁠.³⁷ For odd (conversely, even) m, the density, axial velocity, and hydrodynamic pressure perturbations have odd (even) parity and the radial and azimuthal velocity perturbations have even (odd) parity, effectively halving the dimensions of the operator matrices. For this study, convergence tests showed that a Chebyshev polynomial expansion of order N = 201 was suitable to accurately and affordably compute all of the eigenvalues and eigenvectors of interest (see convergence criteria below).

All of the eigenvalue problems in this work are solved using the QZ and QR algorithm implementations within the Matlab eig and eigs functions, respectively. The temporal stability analysis involves the direct solution of (12) as a generalized eigenvalue problem for ω with k given. The spatiotemporal problem is solved using a variant of the algorithm proposed by Deissler,³⁹ which finds saddle points where $∂ω/∂k=0$ using a gradient-descent algorithm for the group velocity on the complex k–ω manifold. This process is facilitated by manipulating the dispersion relationship of (12) to yield an analytical expression for the group velocity as detailed in  Appendix. Finally, the spatial eigenvalue problem, which seeks the spectrum of complex k with ω given, is solved directly as a generalized eigenvalue problem via the companion matrix method.⁴⁰ 

Regardless of the problem type, even when the computational grid is very highly resolved, the raw eigenvalue spectrum may contain spurious and poorly resolved solutions due to the discretization and floating-point arithmetic.^15,36 These problematic eigenvalues are filtered from the eigenvalue spectrum by applying the following selection criteria: (i) eigenvalues must be associated with axial phase speeds smaller in magnitude than 10⁶, (ii) eigenvectors must vanish as $r→∞$ using the criterion of Parras and Fernandez-Feria³⁶ with a cutoff value of $10−6$⁠, and (iii) eigenvalues must have a difference in phase speed of less than 1% when computed with N = 201 and N = 181.

This investigation will focus solely on the flow through the center nozzle, for which the nominal flow conditions remain consistent across all of the test cases. To identify coherent structures within the measurements, the harmonic content is first analyzed on the basis of its energy spectra. The kinetic energy associated with each frequency is characterized by the spatial average of the squared magnitude of each Fourier component. This quantity represents a power spectral density (PSD) which is shown in Fig. 3. These results show that the unforced flow does not exhibit any clear energy peaks, indicating that the flow is not self-excited. It also demonstrates that transverse acoustic forcing has a controlling effect on the flow dynamics as it introduces strong, narrowband energy content at the forcing frequency. Finally, the spectra indicate that OP forcing yields a stronger response from a kinetic energy perspective than IP forcing. Note that the Fourier modes associated with the marked frequencies in Fig. 3 will henceforth be identified as the coherent part of their respective cases.

Next, the flow and flame motions are considered from a structural perspective. Each case is visualized in Fig. 4 by an instantaneous snapshot of the flow field and the mean and coherent components of the triple decomposition from (1). For the decomposition, instantaneous three-dimensional velocity vectors come directly from the processed sPIV image pairs, and the instantaneous density is deduced from product/reactant fields inferred from the binarized OH-PLIF measurements and a burned/unburned density ratio determined from methane–air chemical equilibrium calculations. All of the data considered in this paper come from the interval $x∈[0.08,1.4]$⁠, which corresponds to the dimension of the window for which reliable simultaneous PLIF and sPIV measurements are available. Several aspects of the data presented in Fig. 4 are worthy of comment. To begin, the figure shows that the mean flow field is essentially unaffected by forcing as the mean flow is similar across all cases. The coherent parts, however, differ greatly. While the unforced flow does not indicate any significant coherent oscillations, the forced cases do present a clear response at the forcing frequency $ωf=3.61$⁠. Specifically, the OP case shows a clear asymmetric response which is indicative of helical shear layer oscillations. These flow oscillations then manifest oscillations of the flame, which is kinematically coupled to the flow, as shown by the staggering of the phase-mean progress variable field (PVF) contours. Conversely, the IP case exhibits a largely symmetric response which is indicative of symmetric helical or ring-like shear layer oscillations. The stronger response of the flow to asymmetric OP-forcing, in comparison to symmetric IP-forcing, suggests an underlying preference of the flow toward asymmetric helical disturbances. This will be discussed further in Sec. IV and V.

The mean flow used for the stability analysis is derived entirely from the unforced single-nozzle experiment. Previous experimental investigations have shown the mean flow from this single-nozzle experiment approximates that of the three-nozzle experiment, and that both cases have similar disturbance profiles.²² The single-nozzle dataset consists of 2727 synchronized sPIV and OH-PLIF measurements, which are used to identify the Favré-averaged mean flow state as described above and shown in Fig. 4. After the mean velocity and flame progress variable fields are determined from the images, they are spatially averaged by reflection across the centerline to yield an azimuthally homogeneous mean flow state (not shown in Fig. 4).

As indicated in Sec. II, the Boussinesq hypothesis is used to model turbulent transport in the stability calculations. The eddy viscosity is determined directly from the measured Reynolds stress tensor using the projection of (7). For the present data set, the overdetermined system yielded a spatially varying eddy viscosity which was slightly negative in some small regions. Though certain conditions do exist where this may be physically significant,^41,42 in this work, the negative eddy viscosity values are set to zero such that $μ¯T=D¯T≥0$⁠, consistent with the approach taken by previous related studies.^17,18,43

A continuous representation of the velocity, density, and eddy viscosity data are obtained by fitting a two-dimensional piecewise-cubic spline curve to the azimuthally averaged and Favré-averaged basic flow fields. A two-dimensional Gaussian smoothing filter was also applied to the $μ¯T$ field. The resulting smooth, axisymmetric basic flow profiles are therefore directly based on the measurements rather than an assumed model equation.

Profiles of the axisymmetric, Favré-averaged mean flow used for the stability analysis are shown in Fig. 5. It shows first that the axial component displays a significant amount of reverse flow in the central recirculation region, as well as very clear inner and outer shear layers. It is well known from inviscid theory that such features heavily influence a flow's stability properties.⁴⁴ Second, the presence of the central recirculation region caused by the combined effects of vortex breakdown and the center body wake causes the average rotating component to rapidly transition from a well-defined azimuthal jet at the dump plane to a more radially distributed vortex profile, even as the axial velocity profile retains its annular jet-like character. Third, the time-average of the binarized OH-PLIF signals forms a mean flame PVF which shows that the flame exists primarily within the inner shear layer. The presence of the flame here has a drastic effect on turbulent transport as measured through the Reynolds stresses, and this effect is clearly reflected in the resultant eddy viscosity.

In the local temporal stability analysis framework, complex-valued frequencies are sought for imposed real wavenumbers. In more physical terms, this form of stability analysis corresponds to interrogating the time evolution of perturbations to the mean flow from the perspective of a moving reference frame which travels with each perturbation. In that sense, the temporal stability analysis indicates a “worst case” scenario for instability, since a perturbation that does not grow in its own reference frame cannot appear to grow in any other reference frame.

The temporal analysis is performed over a range of axial and azimuthal wavenumbers at all axial locations to identify all of the possibly relevant instability modes, which are then set aside for further analysis. Each perturbation has an axial phase velocity $ω/k$⁠, an azimuthal phase velocity $ω/m$⁠, and a helical winding orientation $−k/m$⁠. Henceforth, we consider $m≤0$ without loss of generality, such that $ω>0$ indicates counter-rotating perturbations, k > 0 indicates co-winding perturbations, and vice versa. Complementary to their azimuthal and axial phase velocities and winding orientation, each mode may be further classified according to its radial velocity distribution. However, unlike the phase velocities and helix orientation classifications which have precise definitions, these distinctions are somewhat subjective. A given mode may be a center mode that displays its peak oscillatory amplitude in the central recirculation region, a shear layer mode that peaks within the annular jet's inner and outer shear layers, a ring mode which peaks in the fluid surrounding the jet, or, more generally, a combination of these.^35,38,45

As a representative example, temporal stability analysis results from x = 0.25 are shown in Fig. 6 for $|k|≤20$ and $|m|≤8$⁠. Note that perturbations with higher azimuthal wavenumbers were also analyzed, but are not shown as their growth rates decayed monotonically with increasing $|m|$ until becoming stable for $|m|≥15$ at this location. These results indicate that the most temporally unstable axisymmetric modes have a downstream axial phase velocity (i.e., they are co-propagating), and that the most temporally unstable helical modes rotate and wind in opposite directions. In certain cases, two distinct instability modes exist simultaneously at each k. In all cases, the dominant axial wavenumbers associated with the largest temporal growth rates correspond to wavelengths on the order of one nozzle diameter [i.e., $2π/k∼O(1)$]. The unstable modes all feature significant amplitudes along the annular jet's inner and outer shear layers, but, in certain cases, also exhibit elements of center or ring modes. As such, no attempt to classify these modes in terms of their radial velocity distributions will be attempted.

In the local spatiotemporal stability analysis framework, perturbations are considered with complex-valued frequencies and wavenumbers in order to assess the linear system's impulse response. As discussed earlier, the validity of the chosen spatial stability analysis-based response framework depends on the coherent system being unable to manifest self-excited behavior, such that the external acoustic disturbances introduced at the burner lips represent the sole wavemaker in the domain. The temporal analysis showed that certain perturbations exhibit exponential amplification in time within their own reference frame. However, these perturbations represent propagating dispersive waves with group velocities which may or may not enable them to grow in time at a fixed point of impulse. Thus, the purpose of the spatiotemporal analysis performed here is to confirm that all perturbations decay to zero at any point of impulse in order to rule out the possibility of a self-excited global instability.

A spatiotemporal analysis was performed for all of the temporally unstable modes identified in the previous step at all axial locations. In this framework, it is essential to ensure that all of the identified stationary points $(k0,ω0)$ where $∂ω/∂k=0$ satisfy causality. This is ensured using the Briggs–Bers “pinch-point” criterion, which requires that a relevant saddle corresponds to a stationary wave “pinched” between upstream and downstream propagating waves.⁴ Exemplary results of the spatiotemporal analysis from x = 0.25 are shown in Fig. 7. The analysis reveals that all of the dominant stationary points correspond to situations where $Im{ω}<0$⁠, indicating that the impulse response decays in time at the point of perturbation. Thus, in the local stability analysis vernacular, the flow is—at most—convectively unstable, and cannot sustain the absolute instability necessary to manifest self-excited behavior. This indicates that the signaling problem is well-posed, allowing further analysis in the spatial stability framework considered below.

In the spatial analysis, perturbations are considered with real frequencies and complex wavenumbers. This situation physically represents a signaling problem, where harmonic perturbations are introduced at a point, and the analysis determines whether the perturbations grow or decay along the streamwise direction. A representative spatial stability analysis result from x = 0.25 is shown in Fig. 8 for $|ω|≤20$ and $|m|≤2$⁠. While this analysis has a different physical interpretation compared to the temporal analysis considered earlier, the overall characteristics of the spatially unstable coherent perturbations are qualitatively quite similar to their temporal analogues. In particular, the dominant instabilities still tend to represent helical waves which wind and rotate in opposite directions or axisymmetric waves which propagate downstream. These features are shown to exhibit spatial amplification for a wide range of frequencies and wavenumbers.

Building on the results of the local analysis methods described above, the WKBJ approach can now be used to determine the linear weakly global modes. Following the argument presented earlier, as no natural absolute instability exists to generate sustained intrinsic coherent oscillations, the acoustic coupling which occurs near the burner lips is interpreted to represent the flow's sole wavemaker. This external excitation locally excites small, pure-tone perturbations which may then be spatially amplified by hydrodynamic instability as they convect downstream.

Based on this physical model, the weakly global modes are constructed from a sequence of local analyses conducted at 67 axial locations evenly distributed over the measurement interval $x∈[0.08,1.40]$⁠. The eigenvalues of the spatial stability analysis at the forcing frequency of 390 Hz, which corresponds to a dimensionless angular frequency of $|ωf|=3.61$⁠, are shown in Fig. 9 for each x. Using these independent eigenvalues, the spatial stability results are assimilated through the WKBJ approximation of (11) with $x0=0.08$ to determine the streamwise evolution of each weakly global mode's amplitude and phase. This information is also presented in Fig. 9 for the $|m|≤5$ modes which experience significant spatial amplification. Note that modes with $|m|≥6$ were considered in the analysis, but are not shown as they were all found to be more stable than those presented.

The analysis reveals that the jet exhibits a wide range of spatially amplified weakly global modes associated with coherent oscillations at the forcing frequency, with the strongest spatial amplification exhibited by a co-rotating, counter-winding $|m|=1$ mode. Thus, if the forcing excited all modes equally, the linear response would consist primarily of this $|m|=1$ oscillation along with several distinct hydrodynamic structures with different periodicities, winding orientations, and rotation directions. This result also offers a convincing explanation for the stronger response associated with OP forcing compared to IP forcing as discussed in Sec. III and shown in Fig. 3. Despite this, as will be discussed in the next section, the transverse forcing used in the experiments is strongly selective and does not equally excite all of the spatially amplified modes identified here. As such, the co-rotating $|m|=1$ mode should not be expected to dominate the forced response in all situations.

In addition to the axial amplitude and phase variation discussed above, the weakly global analysis also provides insight into structural characteristics of the linear dynamics. Visualizations of the leading weakly global modes for $|ω|=3.61$ and $|m|≤3$ are shown in Fig. 10. The visualizations reveal a variety of axisymmetric (m = 0) and spiral (⁠$m≠0$⁠) vortical structures which are concentrated along the annular jet and particularly its shear layers. According to Fig. 9, the counter- and co-rotating structures tend to correspond to co- and counter-winding helical vortex structures, respectively. Since the flame stabilizes along the jet's inner shear layer, these mode shapes also suggest a hydrodynamic connection between the acoustic forcing and the flame response.

It should again be emphasized at this point that the weakly global modes considered in Figs. 9 and 10 are based only on the forcing frequency and not on any structural characteristics of the forcing itself. This is a crucial distinction as the transverse acoustic forcing used in the experiments features important structural characteristics which do not evenly excite all disturbances. In particular, the OP and IP cases are designed to preferentially excite modes with odd and even m, respectively. This will be discussed further in the following, where the generic linear modes described in this section are synthesized into response predictions and compared with the measured response under specific excitation conditions.

As mentioned previously, the structure of the IP and OP excitation patterns considered in the forced experiments were selected to elicit specific and distinct flow responses. In the IP case, the nozzle center coincides with a pressure anti-node of a standing planar acoustic wave in the combustor, leading to symmetric fluctuations in pressure with respect to the anti-nodal line. Expanding this planar symmetric fluctuation in an azimuthal Fourier series, it can be shown that IP forcing is also primarily axisymmetric at the central nozzle with only minor excitation of higher even-m disturbances (i.e., $|m|=2,4,…$⁠). Conversely, for the OP case, the nozzle center lies at a pressure node, causing an anti-symmetric pressure response with a 180° phase difference across the nodal line. In this case, the dominant effect of the planar OP forcing is the excitation of two superposed $|m|=1$ waves of equal strength which rotate in opposite directions with minor contributions toward higher odd-m disturbances (i.e., $|m|=3,5,…$⁠).

Hence, even though the general analysis above indicates a wide range of independent, spatially amplified modes, the nature of the acoustic forcing used in the experimental investigation introduces two key effects which both must be taken into consideration by the model. First, as mentioned above, transverse-wave excitation only excites structures with specific periodicities. As shown by previous studies,^2,46 the present jet may be considered to be acoustically compact as the wavelengths associated with the transverse acoustics at $|ωf|=3.61$ are O(100) times larger than the diameter of the nozzle. This assumption allows the higher-order terms associated with the azimuthal Fourier expansions of the IP and OP forcing patterns to be discarded when constructing the response predictions, such that the leading order excitation consists only of m = 0 and $|m|=1$ structures for the IP and OP case, respectively. Second, the various linear modes evolve independently in the domain where the acoustics are decoupled, but are excited by a single pressure wave which is hydrodynamically coupled near the burner lips. Thus, the amplitude and phase $âj$ associated with each weakly global mode $q̂j$ in the overall response mode $q̃$ is determined by projecting the acoustic axial velocity distributions for IP and OP transverse forcing reported by Blimbaum et al.⁴⁶ onto the axial velocity at the most upstream point of each weakly global mode. Finally, the predicted response is evaluated as the weighted sum of individual contributions from each spatially amplified weakly global mode using (14).

As remarked earlier, the experimental response to forcing has been thoroughly detailed.⁷ In part of that analysis, the authors identified and isolated azimuthally periodic structures in each axial plane and charted their spatial evolution as a function of x using helical mode coefficients, $B̂m(x,r)$⁠, which quantify the energy present within a given helical mode denoted by m. These results are recounted in Fig. 11 and compared against the response predictions using the area-weighted helical mode coefficients, $B¯m(x)=∫0∞B̂mr dr$ for the axial component of velocity. It should be noted that the $B¯x,m$ reported from Smith's measurements are not normalized by the input, whereas the predicted response modes are referenced to the amplitude at the most upstream x. Hence, the absolute values of the measurements and predictions are not quantitatively comparable even though the relative changes are.

Based on Fig. 11, the linear predictions match several important characteristics of the observed response. For example, the predicted and measured response under OP forcing both indicate a relative dominance of the co-rotating $|m|=1$ structure in comparison to its counter-rotating analogue as well as a slightly higher initial excitation. This suggests that the co-rotating oscillation is dominant both due to an underlying hydrodynamic instability of the mean flow profile, but also due to a slightly biased excitation of the co-rotating mode by the forcing. Additionally, the growth and decay trends of the predicted modes with respect to the axial coordinate are very similar to the measurements. In particular, the amplitude of the m = 0 response appears to grow with x in the IP case for both prediction and measurement. On the other hand, the $|m|=1$ response modes do not grow uniformly. Both the co- and counter-rotating $|m|=1$ modes exhibit an interval of comparable spatial amplification until $x∼0.6$⁠, beyond which the amplitudes diverge. For the co-rotating mode, the amplitude continues to increase while the counter-rotating mode experiences a decay.

To further aid model evaluation, the spatial distribution of the amplified linear response modes are also compared against the measured coherent components of the test cases with IP and OP forcing in Fig. 12. Under IP forcing, Fig. 11 showed that the measured flow exhibits a primarily axisymmetric response. The spatial distribution of the predicted m = 0 response shown in Fig. 12 is clearly very similar to the measured response. In particular, both consist of axisymmetric vortical structures of similar axial wavelengths (and therefore axial phase speeds $k/ω$⁠) which cascade downstream along the inner and outer shear layers of the mean flow. For ease of interpretation, these vortex ring structures are also visualized by three-dimensional isocontours of the azimuthal component of vorticity in the lower portion of Fig. 12.

In the case of OP forcing, the predicted response consists of both co- and counter-rotating $|m|=1$ modes. As shown in Fig. 11, the co-rotating $|m|=1$ mode dominates in both the experimental and predicted response beyond $x≳0.5$⁠, but the amplitudes of both modes remain significant throughout the domain. Once again, the predictions agree remarkably well with the measurements, showing a similar spatial distribution of velocity both radially and axially. The three-dimensional isocontour representation of the predicted response's azimuthal vorticity field also clearly shows the prevalence of the counter-winding, co-rotating disturbances as the jet evolves downstream. Overall, these comparisons suggest that the linear response predictions derived by projecting the excitation profiles onto weakly global mean-flow stability modes capture all of the major aspects which control the evolution of coherent disturbances in this flowfield.

Before closing, some brief comments are worthwhile to better frame these results. First, it is important to emphasize that the response predictions and measurements are completely independent of each other. The mean flow used in the model is based solely on measurements of the unforced flow, and no information from the forced experiments other than the frequency and type of forcing are provided to the model.

Second, the agreement between the model and experiments suggests that the many assumptions used to reduce the order of the model do not interfere with the physics at play—at least for the situations we tested. This raises several questions regarding the limitations of such a simplified approach, as well as what physical processes must be included and what can be neglected. For example, earlier work⁷ has shown that nonlinear interactions among the fundamental oscillation at the forcing frequency and its super- and sub-harmonics do not play a significant role in the coherent dynamics of swirling jets, a result which should be contrasted with non-swirling jets¹² where nonlinear interactions are more prominent. This supports the validity of the linearity assumption used in our analysis, which lumps the high-dimensional turbulent flow dynamics into only mean density, velocity, and eddy viscosity fields. As a consequence, this model neglects the nonlinear coherent–coherent interactions which arise through the triple decomposition and only permits incoherent–coherent interactions through the eddy viscosity.⁴⁷ Nonetheless, the success of this assumption is not particularly surprising as it has been investigated thoroughly in previous nonreacting flow studies.^10,17,43

Another important assumption made in this study was that of a weakly nonparallel, axisymmetric mean flow, which neglects streamwise boundary effects. The influences of these assumptions on stability and forced response of these flows are less well-explored. The streamwise boundary effects could be particularly relevant to swirling flows, where the effect of streamwise confinement of inertial waves has been linked to onset of vortex breakdown.^28,48,49 Nonetheless, the agreement between simplified model and measurements reported here suggests that these effects do not exhibit prominent effects on the coherent response to forced disturbances in our configuration.

In this paper, a low-order model was developed and assessed for its ability to predict the response of an amplifier-type turbulent reacting swirling annular jet to external acoustic forcing. The measured mean flow field associated with an unforced test case was used as the basis for linear hydrodynamic stability analysis. Consistent with experimental observations, the stability analysis identified the flow to behave as a convectively unstable system, and a sequence of local stability analyses at the forcing frequency were used to identify weakly global modes as a basis for the predicted response. Finally, linear predictions of the flow response under IP and OP transverse acoustic excitation patterns were constructed by projecting the acoustic forcing profiles described by⁴⁶ onto the basis of weakly global modes.

The results of the model explained two key observations reported in previous experimental investigation.^7,22,23 First, it indicates that the observed dominant response of OP forcing in comparison to IP forcing of equal amplitude can be understood based on the flow's natural preference toward spatial amplification of a co-rotating $|m|=1$ mode. Second, it suggests that the experimentally observed dominance of the co-rotating $|m|=1$ mode in comparison to its counter-rotating analogue under OP forcing should be attributed not only to this natural preference, but also to a slightly favorable excitation of the co-rotating mode by the forcing due to its structural characteristics near the nozzle.

This work was partially supported by the Air Force Office of Scientific Research under award number FA9550-20-1-0215 (contract monitor Dr. Chiping Li). The authors thank Travis E. Smith for sharing the experimental data, which form the basis of this work.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

To derive the algorithm for the spatiotemporal analysis, (12) is differentiated with respect to k to yield

where $J=A0+kA1+k2A2$⁠. Then, after left-multiplying by the adjoint (left) eigenvector $q̂Lm†$ and using the identity $q̂Lm†H[−iωM+J]q̂Lm=0$⁠, the expression (A1) can be rearranged to find

Using a similar approach, the following expression for $∂2ω/∂k2$ can be derived

where $∂q̂Lm/∂k$ is determined from (A1) as

and $[·]+$ denotes the Moore–Penrose pseudoinverse. In practice, (A4) is solved using the lsqminnorm function in Matlab.

Based on the above equations, an iterative algorithm to converge to a saddle point from a given initial wavenumber k and a guessed initial eigenvalue ω₀ proceeds as follows. First, (12), (A2), and (A3) are solved to obtain ω, $q̂Lm, q̂Lm†, ∂ω/∂k$⁠, and $∂2ω/∂k2$ at the current point k. Next, the modulus and argument of the gradient-descent step are deduced as

and

respectively, where $δkmax>0$ is a step length threshold, which aids numerical stability by avoiding excessively large steps. Finally, the step is completed by taking $k←k+δk$⁠, and the process is repeated until convergence $|∂ω/∂k|<ϵ$ is reached for some tolerance $ϵ>0$⁠.