Linear analysis of a swirling jet with a realistic swirler model

Abstract

The dynamics of an axisymmetrical swirling jet is studied via global linear stability and resolvent analyses. The modeled flow represents a combustor-like swirling jet, that is turbulent, compressible, non-parallel, and enclosed. In particular, the computational domain embeds a realistic axisymmetrical swirler model to resolve the mode conversion process. Swirl fluctuations are non-negligible on this configuration representative of a swirl burner, and match the analytical mode shapes of inertial waves of an inviscid uniform flow as obtained from global stability analysis. The stability map presents two eigenvalues driving a modal amplification. These eigenmodes couple a standing acoustic wave sustained in the mixing duct and the combustion chamber with the Kelvin-Helmholtz mechanism at the mixing duct exit and the acoustic-vorticity mode conversion process at the swirler, and act as a frequency selection criterion. Finally, the most amplified forcing from the resolvent analysis is similar to an unsteady heat source in the combustion chamber, and the identified optimal amplification mechanism is likely to be triggered in reacting flow with unsteady heat release rate.

Introduction

Swirl-stabilized flames are widely used in combustion technologies to anchor flames. The formation of a central recirculation bubble downstream of the area expansion in a combustion chamber provides a stabilization region without contact with a solid part. This recirculation region is also beneficial for mixing reactive radicals with the unburnt mixture, thereby enhancing flame stabilization. ¹ Swirl stabilization is used for both premixed and non-premixed flames. An upstream swirler generates the swirling jet motion by converting axial or radial momentum into azimuthal momentum. Acoustic waves impinging on the swirler generate swirl fluctuations through momentum conservation, in a process called mode conversion. ² These swirl fluctuations are propagated downstream by an underlying inertial wave: in a column of fluid in rotation, velocity fluctuations in the azimuthal and radial directions couple through the Coriolis force and propagate in the axial direction.^3,4 Swirl fluctuations in swirl stabilized combustion systems have been evidenced experimentally.^2,5 and through Computational Fluid Dynamics (CFD)^6,7

In many industrial applications, flames are velocity-sensitive. This means that velocity fluctuations induce heat release rate oscillations which generate acoustic waves. These acoustic waves may, in turn, produce velocity fluctuations that provide feedback to the flame. This thermoacoustic feedback loop ⁸ can turn unstable and generate instability. Such thermoacoustic instabilities pose huge limitations to the safe and reliable operations of premixed combustion devices, ⁹ particularly for low-emission engines operated in the lean regime.

The mechanism by which inertial waves impact the flame is not yet fully established, despite undeniable progress based on experimental, analytical, and simulation studies. Indeed, the main characteristic of swirl waves is to carry azimuthal velocity fluctuations. However, if swirl stabilized flames exhibit a rotational symmetry, azimuthal velocity fluctuations happen to be tangential to the mean flame sheet and, therefore, cannot induce heat-release rate fluctuations by a change of the flame surface area However, azimuthal velocity fluctuations can induce heat-release rate fluctuations through flame speed fluctuation. ¹¹

If not a direct impact, azimuthal velocity fluctuations can have an indirect impact on the flame by triggering axial and radial velocity fluctuations upstream from the flame. Since axial and radial velocity fluctuations have a component normal to the flame sheet, they eventually induce heat release rate fluctuations. Some transfer mechanisms from azimuthal to radial and axial components, have been proposed, for example, by Acharya and Lieuwen, ⁷ who suggest that vortex tilting downstream the area expansion generates radial velocity fluctuations. However, Albayrak et al. ⁴ found that inertial waves also carry axial and azimuthal velocity fluctuations. These results were established on an analytical model based on an inviscid parallel flow under solid body rotation. The presence of axial and radial components brings an additional element to the feedback mechanism since it enables inertial waves to have a direct impact on the flame sheet. The importance of axial and radial components of swirl waves was established by Albayrak et al., ¹² where they dominate the flame front modulation. These results follow from applying the linearized reactive flow ¹³ model on a laminar swirl stabilized flame, with idealized swirl fluctuations imposed at the inlet and integrated in time.

Recently, Müller et al. ¹⁴ investigated the behavior of inertial waves in a turbulent flow by resolvent analysis and spectral proper orthogonal decomposition of experimental measurements for a radial swirler configuration. They confirm the existence of inertial waves as a driver for swirl fluctuations in turbulent duct flows, but of small amplitude compared to shear-driven effects. The latter comment has to be balanced because shear is stronger in radial than in axial swirler. Palies et al. ¹¹ showed that downstream of the duct, swirl fluctuations dominate over shear-layer contributions during the development of the instability and once the limit cycle is reached. In summary, although their existence is demonstrated, the mechanisms by which inertial waves impact the flame-flow dynamic still lack a comprehensive picture.

In addition to the type of velocity fluctuations carried by inertial waves, their propagation velocity is another crucial parameter for thermoacoustic instabilities. Indeed, thermoacoustic instability is at work when pressure and heat-release rate fluctuations are in-phase. Therefore, the time delay of the feedback through inertial waves, tightly linked to their propagation velocity, is a key parameter. Komarek and Polifke ⁵ actually found that the propagation speed of inertial waves exceeds the convection velocity, by up to 40%. In addition, the study of Albayrak et al. ⁴ demonstrated the dispersive behaviour of inertial waves. That study unveiled the existence of fast and slow traveling modes, respectively, which depart from the convective velocity depending on the swirl intensity. The amount by which the fast and slow modes depart from the convection velocity grows with the swirl intensity. The dominance of one of these two modes is relevant to controlling thermoacoustic instabilities in combustion devices since it impacts the time delay. Consequently, further knowledge of the generation of inertial waves is needed. In the present study, the swirler is materialized as source terms in the governing equations in lieu of an explicit geometrical treatment into the computational domain. The purpose is to retain a representative flowfield associated with mode conversion process yet without the swirler geometrical part itself, which is three dimensional by essence. Since three-dimensional linear analysis remains computationally expensive and is out of reach for the size of the configuration of interest, there exists, to date, no linear analysis of swirling flows that includes a swirler section. This one is usually modeled with ad-hoc boundary conditions, for the base- or meanflow and the fluctuations. In the present work, we propose to include the swirler section in the analysis with the versatile axisymmetrical swirler model proposed by Kiesewetter et al. ¹⁰ This model proved itself capable of describing axial and radial swirlers. Building on that model, we propose a method for an integrated operator-based analysis of a swirling jet representative of an axial swirl burner, that is, turbulent, compressible, enclosed, and non-parallel. As a result, the study reports on the generation and transport of inertial waves in such a flow, as well as their role in the mechanisms that optimally amplify small perturbations.

The article is organized as follows: the configuration, models, and numerical choices are introduced in the “Modeling and numerical framework” section. The generation and propagation of inertial waves is presented in the “Propagation of inertial waves in non-parallel flows” section. In the “Amplification mechanisms” section, the principle of linear stability and resolvent analysis are introduced and the identified amplification mechanisms are then presented. These results are then discussed in the context of combustion dynamics before concluding.

Modeling and numerical framework

The configuration of interest depicted in Figure 1(a) and Table 1, is a cylindrical burner composed of an axial swirler with a central bluffbody and a flow expansion into a combustion chamber. This configuration is studied with an axial symmetry, in a cylindrical coordinate system $x$ , $r$ , and $\theta$ , for the axial, radial, and azimuthal directions, respectively. Since the azimuthal Fourier modes of order greater or equal to 1 do not affect the global heat release rate fluctuation of an axisymmetric swirling flame in the regime of small fluctuations, only the azimuthally symmetric mode is preserved, and the dependence in $\theta$ is dropped. The impact of azimuthal velocity fluctuations on the global heat-release rate through the flame speed, as evidenced by Palies et al., ¹¹ also cancels for Fourier modes greater than one. The flow is described by three state variables: the density $\rho$ , the velocity vector $\mathbf{u}$ , and the enthalpy $h$ , all gathered in a vector variable $q=[\rho ,\mathbf{u},h]$ .

Figure 1.

(a) Configuration of the axisymmetrical swirling jet. The dashed-circled portion represents the swirler section, where a the swirler model of Kiesewetter et al. ¹⁰ is applied. Dimensions are given in Table 1. (b) Baseflow velocity fields. Upper half: ${\overline{u}}_x$ with flow lines in the $(x,r)$ plan and lower half: ${\overline{u}}_{\theta }$ .

Table 1.

Dimensions of the configuration of Figure 1(a), inlet reflection coefficient R_in, swirl number Sw, Reynolds number Re, and constant $A_0$ from Kiesewetter’s et al. model ¹⁰ ((1) and (2)). The first row corresponds to the base case, and the second row corresponds to the independent parameter variations.

	$l$ chamber [mm]	$l$ tube [mm]	$l$ swirl [mm]	$r$ in [mm]	$r$ out [mm]	$r$ chamber [mm]	$R_{\mathrm{in}}$	$\mathrm{Sw}$	Re	$A_0$
Base config.	300	168	4	8	20	45	1	0.75	4.104	25
Param. variation	400	126					0

Bold values represent the base case.

Kiesewetter et al.^10,16 proposed to model the axial to azimuthal momentum conversion in a swirler by a coupling term added to the azimuthal momentum conservation equation in a restricted region of the mixing duct, depicted in Figures 1 and 2,

$$F_{\theta }=A\rho u_x^2$$

(1)

The prefactor $A$ is non-zero only in the swirler section. This volumetric force in the azimuthal direction leads to an unphysical volumetric power input $P_{\theta }=F_{\theta }{u}_{\theta }=A\rho u_x^2{u}_{\theta }$ . Therefore, this power addition in the azimuthal direction must be compensated by a volumetric force $F_x$ in the axial direction, such that $F_{\theta }{u}_{\theta }+F_x{u}_x=0$ . Adding a coupling term,

$$F_x=-A\rho u_x{u}_{\theta }$$

(2)

in the axial direction allows to preserve a zero net power addition to the flow. The prefactor $A$ is space-dependent and controls the intensity of the resulting swirl motion. Figures 3 and 4 in Kiesewetter et al. ¹⁰ and Fig. 61 in Kiesewetter ¹⁶ display the quantitative agreement of the axisymmetrical swirler model (equations (1) and (2)) with experimental and large eddy simulation (LES) velocity profiles.

Figure 2.

Radial profile of the base-flow downstream the swirler: time-averaged LES of a 3D resolved swirler ¹⁵ (dashed) and the present axisymmetrical model (solid). The region between the green dashed lines corresponds approximately to a solid body rotation. LES: large eddy simulation; 3D: three-dimensional.

Figure 3.

Zoom on the mixing tube. Azimuthal velocity fluctuation ${\hat{u}}_{\theta }$ for an axial velocity forcing at the inlet, at 300 Hz. The body forcing terms of the swirler model are applied in the region circled in black. The blue and black arrows indicate the two radii mentioned in Figure 7.

Figure 4.

Dispersion diagram of ${\hat{u}}_{\theta }$ : estimated wavenumber along the dashed center line of Figure 3) against the frequency of the inlet velocity forcing.

The state variables evolve according to Navier-Stokes equations, here in cylindrical coordinates,

$${\mathrm{\partial }}_t\rho +\frac{1}{r}{\mathrm{\partial }}_r(r\rho u_r)+{\mathrm{\partial }}_x(\rho u_x)=0$$

(3\rm a)

$$\begin{array}{rl}& {\mathrm{\partial }}_t(\rho u_r)+\frac{1}{r}{\mathrm{\partial }}_r(r\rho u_r^2)+{\mathrm{\partial }}_x(\rho u_r{u}_x)-\frac{1}{r}\rho u_{\theta }^2+{\mathrm{\partial }}_{r}p\\ & ={\mathrm{\partial }}_z{\tau }_{rz}+\frac{1}{r}{\mathrm{\partial }}_r(r{\tau }_{rr})-\frac{{\tau }_{\theta \theta }}{r}\end{array}$$

(3\rm b)

$$\begin{array}{rl}& {\mathrm{\partial }}_t(\rho u_x)+\frac{1}{r}{\mathrm{\partial }}_r(r\rho u_x{u}_r)+{\mathrm{\partial }}_x(\rho u_x^2)+{\mathrm{\partial }}_{x}p\\ & ={\mathrm{\partial }}_x{\tau }_{xx}+\frac{1}{r}{\mathrm{\partial }}_r(r{\tau }_{rx})-A\rho u_x{u}_{\theta }\end{array}$$

(3\rm c)

$$\begin{array}{rl}& {\mathrm{\partial }}_t(\rho u_{\theta })+\frac{1}{r^2}{\mathrm{\partial }}_r(r^2\rho u_r{u}_{\theta })+{\mathrm{\partial }}_x(\rho u_x{u}_{\theta })\\ & ={\mathrm{\partial }}_z{\tau }_{\theta z}+{\mathrm{\partial }}_r{\tau }_{r\theta }+\frac{2{\tau }_{r\theta }}{r}+A\rho u_x^2\end{array}$$

(3\rm d)

$$\begin{array}{rl}& {\mathrm{\partial }}_t(\rho h-p)+\frac{1}{r}{\mathrm{\partial }}_r(r\rho u_{x}h)+{\mathrm{\partial }}_x(\rho u_{x}h)\\ & ={\mathrm{\partial }}_x{\varphi }_x+\frac{1}{r}{\mathrm{\partial }}_r(r{\varphi }_r)\end{array}$$

(3\rm e)

with $r$ is the radial coordinate, $p$ is the pressure, and $\underset{\mathbf{\_}}{\mathit{\tau }}$ is the viscous shear stress tensor,

$$\underset{\mathbf{\_}}{\mathit{\tau }}=\mu (\mathrm{\nabla }\otimes \mathit{v}+(\mathrm{\nabla }\otimes \mathit{v}{)}^t-\frac{2}{3}(\mathrm{\nabla }\cdot \mathit{v})\underset{\mathbf{\_}}{\mathit{I}})$$

(4)

The dynamic viscosity $\mu$ is obtained through Sutherland’s law. ¹⁷ The vorticity in the plan ( $x,r$ ) is defined as $\omega ={\mathrm{\partial }}_x{u}_r-{\mathrm{\partial }}_r{u}_x$ . The heat conduction flux follows from Fourier’s law, $\mathit{\varphi }=-\alpha \mathrm{\nabla }T$ , where the heat diffusivity $\alpha$ is derived from the viscosity with a Prandtl number $\text{Pr}=0.85$ . Pressure $p$ , temperature $T$ , and density $\rho$ are related to each other by the equation of state of an ideal gas. Equation (3) is recast in the compact form,

$${\mathrm{\partial }}_{t}q=\mathcal{L}(q)$$

(5)

The flow is turbulent and hence decomposed using the triple decomposition. ¹⁸ The flow variables are decomposed into three parts,

$$q=\overline{q}+\tilde{q}+q^{\prime }$$

(6)

The mean part $\overline{q}$ is solution to the steady Reynolds-averaged Navier-Stokes (RANS) equation,

$${\mathcal{L}}_{\mathrm{RANS}}(\overline{q})=0$$

(7)

where the Reynolds stress and turbulent diffusion fluxes are modeled via the k- $\omega$ shear-stress transport model. ¹⁹ The third term, $q^{\prime }$ , is a stochastic fluctuation of zero temporal and phase average, that is, $\overline{q^{\prime }}=⟨q^{\prime }⟩=0$ . The second term $\tilde{q}$ is the harmonic fluctuation defined as

$$\tilde{q}=⟨q-\overline{q}⟩$$

(8)

The harmonic fluctuation is assumed to evolve linearly around the mean state $\overline{q}$ , that is, the evolution equation for $\tilde{q}$ ,

$${\mathrm{\partial }}_t\tilde{q}={\mathsf{L}}_{\overline{q}}\cdot \tilde{q}$$

(9)

stems from the linearization of (5). The details of the linear operator ${\mathsf{L}}_{\overline{q}}$ are given by Meindl et al. ²⁰ All walls are no-slip and non-penetrating, the outlet is non-reflecting, and the inlet has an acoustic reflection coefficient $R_{\text{in}}$ given in Table 1. Equation (7) is solved with finite volumes on a ${34.10}^3$ cells mesh via the semi-implicit method for pressure linked equations (SIMPLE) algorithm in OpenFoam. ²¹ A pseudo-time stepping approach finds the solution. The resulting base velocity fields $\overline{\mathbf{u}}$ are shown in Figure 1(b). Equation (9) is discretized in space using the discontinuous Galerkin finite element method (DG-FEM) approach implemented in the in-house code FelicitaX. ²⁰

Assessment of the axisymmetrical swirler model

Linear fluctuations solution to (9) for an input forcing are verified against 3D LES results. Equation (9) is Fourier-transformed in time,

$$i\mathrm{\Omega }\hat{q}={\mathsf{L}}_{\overline{q}}\cdot \hat{q}$$

(10)

The solution $\overline{q}$ is computed at $Re={4.10}^4$ and swirl number $S=0.7$ , with $A=A_0\sqrt{r}$ ( $A_0$ is given in Table 1). This choice of parameters mimics the radial profile of azimuthal velocity downstream the 3D resolved swirler of the BRS burner ⁵ which has the same dimensions (Figure 2). In this flow, an inner region of approximate solid body rotation is identified, as recently observed by Müller et al. ¹⁴ The linearized perturbations $\hat{q}$ are solved on the baseflow of Figure 2 for the axial velocity fluctuations ${\hat{u}}_{x,\text{in}}$ applied at the inlet.

Azimuthal velocity fluctuations corresponding to swirl waves are indeed developing in the swirler section (Figure 3). The wavelength and wavenumber $k$ of ${\hat{u}}_{\theta }$ along the center line (dashed line in Figure 3) is measured for various frequencies and reported in a $k-\omega$ plot in Figure 4. The slope of this curve agrees with the slope of the phase-frequency plot obtained from numerical simulation by Acharya and Lieuwen ⁷ : both give the propagation speed of ${\hat{u}}_{\theta }$ under the assumption that the bulk flow velocity is much lower than the acoustic speed. We conclude that the axisymmetrical model of the swirler can reproduce both the static and dynamic behavior of a realistic swirler.

Propagation of inertial waves in non-parallel flows

In this section, the velocity fluctuations from (9) subject to an inlet axial velocity forcing are compared against wave solutions from Albayrak et al., ⁴ that reports on the dispersive behavior of swirl waves on a uniform inviscid flow. However, Albayrak et al. ⁴ considered that inertial waves propagate as planar waves in the axial direction but the present results are spatially resolved in the $(x,r)$ plan.

To compare the present results with the planar waves of Albayrak et al., ⁴ the perturbation fields $\hat{q}$ (i.e. Fourier-transformed in time) are Fourier-transformed in the axial direction along the center line of the mixing duct (dashed line in Figure 3). The resulting Fourier coefficients are plotted for the three velocity components in Figure 5.

Figure 5.

Absolute value of the spatial Fourier coefficients of ${\hat{u}}_{\theta }$ along the center line of the mixing tube at ${10}^3$  Hz. The vertical line denotes the convective wavenumber.

The three velocity components display a peak around the convective wavenumber ¹ $k$ _conv, indicating that the inertial wave impacts all components of the velocity and not only the azimuthal component. The broadness of these peaks suggests that these inertial waves are not monochromatic and, henceforth, have a dispersive behavior, as reported by Albayrak et al. ⁴ The maxima are all shifted to wavenumbers smaller than the convective one, showing that faster modes are preferred in this configuration. For the sake of the interpretation of Figure 5, the authors note that the first peak on ${\hat{u}}_x$ at $k<100{\mathrm{m}}^{-1}$ corresponds to the acoustic waves traveling in the mixing duct.

The radial profiles of velocity fluctuations $\hat{\mathbf{u}}$ are then compared to the mode shapes of inertial waves derived by Albayrak et al. ⁴ on a uniform inviscid flow. These mode shapes are defined as,

$${\hat{u}}_r^n(r)={\hat{u}}_{\theta }^n(r)=M_1({\lambda }_{n}r),\text{and}{\hat{u}}_x^n(r)=M_0({\lambda }_{n}r)$$

(11)

with $M_i$ being the Bessel function of $i$ -th order, $n$ is the mode number, and ${\lambda }_n$ is the eigenvalue of mode $n$ dictated by the non-penetrating condition at the walls. ⁴ To carry out the comparison, we define the projection operator of two fields— $a$ and $b$ —depending on the radial coordinate as,

$$⟨a(r),b(r)⟩={\int }_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}a(r)b(r)r\text{d}r$$

(12)

with $r_{\mathrm{in}}$ and $r_{\mathrm{out}}$ spanning the entire annulus radius, or restricted to the region of solid body rotation depicted in Figure 2. The perturbations ${\hat{u}}_x$ , ${\hat{u}}_r$ , and ${\hat{u}}_{\theta }$ are projected on ${\hat{u}}_x^n$ , ${\hat{u}}_r^n$ , and ${\hat{u}}_{\theta }^n$ , respectively,

$${\mathrm{proj}}_i^n=\frac{⟨{\hat{u}}_i,{\hat{u}}_i^n⟩}{\Vert {\hat{u}}_i\Vert \Vert {\hat{u}}_i^n\Vert },\text{with}\Vert a\Vert =\sqrt{⟨a,a⟩}$$

(13)

for $i=x$ , $r$ , and $\theta$ . Those projections displayed in Figure 6 show how much each of the velocity components ${\hat{u}}_i$ align with the successive mode shapes ${\hat{u}}_i^n$ .

Figure 6.

Projection of $\hat{\mathbf{u}}$ on radial mode shapes defined by Albayrak et al., ⁴ for an inlet forcing at 100 Hz. Legend “solidbody” corresponds to the case for which only the region of solid body rotation defined in Figure 2 is considered, “full” corresponds to the whole extent of the radius.

The projection of the azimuthal velocity fluctuation decays regularly, with the exception of an odd-even mode oscillation when the projection is made on the region of solid body rotation. This indicates that the first modes of equation (11) capture well the azimuthal velocity fluctuations. The steeper decay for the projection restricted to the inner region of approximate solid body rotation reveals that this agreement is even better in that region, which is expected since the mode shapes of equation (11) are derived on a uniform flow with solid body rotation. The same observations do not hold for the axial and radial components of the velocity fluctuations. In particular, the agreement of the axial velocity perturbation and the basis of mode shapes is worse in the region of approximate solid body rotation than when the entire radius is considered. The projection of the radial velocity fluctuation does not show smooth decay. Thereby, we conclude that the axial and radial velocity fluctuations associated with inertial waves in the present non-uniform and viscous flow do not correspond to axial and radial mode shapes of equation (11).

A possible reason for that observation comes from the non-parallel effects that affect the axial and radial components, that is, axial acoustic waves and shear in the boundary layer. Since the flow is non-uniform, these mechanisms can interplay in the linear regime. This is supported by the weaker agreement of ${\hat{u}}_{\theta }$ with the mode shapes of Albayrak et al. ⁴ when the entire radius is considered for the projection rather than a restriction to the region of approximate solid body rotation.

Fast and slow mode

In the axial direction, however, the velocity fluctuation $\hat{\mathbf{u}}$ does not purely evolve as the planar wave described by Albayrak et al. ⁴ Instead, the wavy structures are tilted in the axial direction; see Figure 3. This is owed to the non-uniform axial mean flow due to viscous effects at walls. In this section, we strive to distinguish by a local analysis which of the fast or slow mode^4,12 dominates the propagation of inertial waves.

If the characteristic length of ${\hat{u}}_{\theta }$ —measured by its wavelength—is much smaller than the characteristic length of the base flow—measured by its axial gradient length $\overline{q}/{\mathrm{\partial }}_x\overline{q}$ —the local parallel flow assumption holds. Under that assumption, the phase shift $\varphi (x)$ of ${\hat{u}}_{\theta }$ locally behaves as

$$\varphi (x)\approx -kx+\text{cst.}$$

(14)

with its local wavenumber $k$ in the axial direction. This condition is fulfilled in the present case for forcing frequencies above 300 Hz. The local wavenumber is then estimated as the axial gradient of the phase shift equation (14), that is, $k$ _est $={\mathrm{\partial }}_x\varphi (x)$ .

This estimated wavenumber $k$ _est is plotted along the radius at the location drawn in Figure 3: at the swirler section exit and close to the mixing duct exit, respectively, and compared to the convective wavenumber $k$ _conv in Figure 7. For $k$ _est < $k$ _conv , the wavelength associated with the inertial wave is smaller than the convective wavelength, therefore, inertial waves propagate slower than the bulk flow at that location. For $k$ _conv, it is the opposite, and inertial waves propagate faster than the bulk flow. Applying this reasoning at the swirler exit (blue curves in Figure 7) shows that inertial waves propagate faster than the mean flow velocity, so the fast mode identified by Albayrak et al. ⁴ dominates. The opposite is observed close to the mixing duct exit (black curves in Figure 7), where the slow mode dominates.

Figure 7.

Solid: estimation of the axial wavenumber of ${\hat{u}}_{\theta }$ from equation (14). Dashed: convective wavenumber $k$ _conv corresponding to the convective speed. Both are plotted along two radii of the mixing tube depicted in Figure 3: swirler outlet (blue) and mixing tube exit (black).

This inversion of the fast and slow modes can be interpreted as a faster dispersion of the fast mode or amplification of the slow mode in the flow direction. However, larger wavenumbers tend to be damped faster by diffusion than smaller ones, so the effect of viscosity on the diffusion of small convective structures does not give a comprehensive explanation. This inversion of the dominant mode can result from the distortion of the azimuthal velocity profile ${\overline{u}}_{\theta }$ by boundary layers developing in the mixing duct. Since it is known that the propagation of inertial waves is affected by the swirl intensity, a local change in the azimuthal velocity of the base flow may disperse or enhance each of these two modes.

Amplification mechanisms

In this section, the role of inertial waves in amplification mechanisms is investigated through stability and resolvent analyses of linear fluctuations $\tilde{q}$ .

Linear stability analysis

Equation (9) is Laplace-transformed and takes the form of an eigenvalue problem,

$$s_o{q}_o={\mathsf{L}}_{\overline{q}}\cdot q_o$$

(15)

with $s_o$ is the complex eigenvalue and $q_o$ is the corresponding eigenvector. The stability analysis is global because it does not have recourse to the parallel flow assumption: the operator ${\mathsf{L}}_{\overline{q}}$ is spatially discretized over the 2D computational domain for $\overline{q}$ solution to equation (7), without further assumption. Temporal solutions developing from eigenmodes evolve as

$$\tilde{q}(x,t)=\exp (s_{o}t)q_o(x)$$

(16)

Therefore, $\mathrm{Re}(s_o)$ is the growth rate of the solution and $\mathrm{Im}(s_o)$ is its pseudofrequency. The eigenvalues $s_o$ are depicted in the stability map Figure 8 and show that the flow is globally stable. Two eigenvalues ( $s_{\text{A}}$ and $s_{\text{B}}$ ) dominate the spectrum at medium frequency. Unlike the rest of the spectrum, $s_{\text{A}}$ and $s_{\text{B}}$ are affected by dimension changes: a shortening of the mixing duct shifts them to higher frequencies, and an increase of the combustion chamber length shifts them to lower frequencies. Finally, considering a non-reflecting inlet ( $R_{\mathrm{in}}=0$ ) dampens them. This suggests that $s_{\text{A}}$ and $s_{\text{B}}$ involve acoustic standing modes developing in the mixing duct and/or combustion chamber.

The loci of $s_{\text{A}}$ and $s_{\text{B}}$ are retrieved by an acoustic network model (circles in Figure 8). The configuration is represented as a network model with interconnected state-space models of 1D acoustic elements ²² depicted in Figure 9. Such a network model does not involve any other feature than the acoustic characteristics of the system. In particular, most hydrodynamic-related effects are neglected. The good agreement between the eigenvalues obtained from the LNSE (15) and network model approaches confirms our interpretation based on the LNSE regarding $s_{\text{A}}$ and $s_{\text{B}}$ : these two eigenvalues are acoustically driven and do not depend on the hydrodynamics. In addition, the mode shapes of these two eigenvalues computed by the acoustic network model show that these two eigenmodes can be approximated by a quarter-wave mode in the mixing duct with a velocity antinode at the duct exit. The importance of these two eigenvalues in the amplification will inform us on the role of acoustics in the dynamic behaviour of the system.

Figure 8.

Stability map for the base case (filled stars) and the variations (empty stars), see Table 1. The two acoustic-like eigenvalues are denoted $s_{\text{A}}$ and $s_{\text{B}}$ . When a parameter is varied, the arrows denote the shift of $s_{\text{A}}$ and $s_{\text{B}}$ . The empty circles indicate the loci of $s_{\text{A}}$ and $s_{\text{B}}$ obtained from the interconnected state-space models.

Figure 9.

Acoustic network representation of the configuration of Figure 1(a).

Frequency selection in the mixing tube

Owing to the fact that all eigenvalues of equation (15) have a negative real part, none of the associated eigenmodes will be self-amplified. In other words, this flow behaves as an amplifier rather than an oscillator: it amplifies fluctuations when forced. As a consequence, the asymptotic stability analysis gives only a limited picture of the dynamic behavior of the system. This amplification behavior shows up in a resolvent analysis that unveils the interaction of the eigenmodes under a forcing. Due to non-normality, the stable eigenmodes interact together to build up the amplified response.^23,24

In that perspective, we now consider the flow under an a harmonic force $\hat{f}$ , in the axial ( $\hat{f}={\hat{f}}_x$ ) or azimuthal ( $\hat{f}={\hat{f}}_{\theta }$ ) direction,

$$i\mathrm{\Omega }\hat{q}-{\mathsf{L}}_{\overline{q}}\cdot \hat{q}=\hat{f}$$

(17)

This forcing is applied either in the combustion chamber or in the mixing duct.

The resolvent analysis solves the following optimal amplification problem at various frequencies $\mathrm{\Omega }$ ,

$$\begin{array}{rl}G(\mathrm{\Omega })& =max\{\Vert \hat{q}{\Vert }_q\}.\\ & \Vert \hat{f}{\Vert }_f=1\end{array}$$

(18)

with $\Vert \hat{f}{\Vert }_f$ is the L-2 norm of the forcing and $\Vert \hat{q}{\Vert }_q$ is the Chu norm of the response. ²⁵ It was first introduced by Trefethen and Trefethen ²⁴ and yielded amplification results in the field of hydrodynamics.^26–28 Equation (17) is discretized with the previously mentioned DG-FEM approach

$$R\hat{q}=B\hat{f}$$

(19)

with $R$ is the matrix resulting from the discretization of the contiuous operator $i\mathrm{\Omega }-{\mathsf{L}}_{\overline{q}}$ , and accounting for boundary conditions and $B$ an extensor built to match the forcing $\hat{f}$ onto the relevant degrees of freedom. The norm of the forcing and the response are equivalently expressed through weighting matrices: $\Vert \hat{f}{\Vert }_f^2={\hat{f}}^t{Q}_f\hat{f}$ and $\Vert \hat{q}{\Vert }_q^2={\hat{q}}^t{Q}_q\hat{q}$ . The optimization problem (18) is recast into the generalized Hermitian eigenvalue problem,

$$(R^{-1}B{)}^t{Q}_q{R}^{-1}B\hat{f}=G^2(\mathrm{\Omega })Q_f\hat{f}$$

(20)

The largest eigenvalue is the optimal gain, associated with the optimal forcing $\hat{f}$ _opt. The optimal response $\hat{q}$ _opt is found by solving the linear system equation (19) with $\hat{f}$ _opt on the right-hand-side.

The optimal gains for the four configurations of interest show a one- or two-peak resonance Figure 10(a) at frequencies corresponding to the eigenfrequencies—that is, the imaginary parts—of $s_{\text{A}}$ and/or $s_{\text{B}}$ (see Figure 8). The frequencies at which these peaks occur are also shifted by the same amount as $s_{\text{A}}$ and $s_{\text{B}}$ for changes in geometrical parameters. This indicates that these resonances are driven by the presence of the eigenvalues $s_{\text{A}}$ and $s_{\text{B}}$ close to the real axis: when the forcing frequency approaches Re( $s_{\text{A}}$ ) and Re( $s_{\text{B}}$ ), it excites the corresponding eigenmodes that drives the amplification. It is a modal amplification.

Figure 10.

Optimal gains: (a) $f=f_x$ in the combustion chamber for various geometrical parameters (see Table 1) and (b) forcing in the mixing tube for the base case with $f=f_x$ or $f_{\theta }$ and the inlet closed or non-reflecting.

This maximum amplification occurs when the forcing is applied on the axial momentum in the combustion chamber. The optimal forcing as shown in Figure 11(a) is distributed in the entire combustion chamber, with a spatial dependence mainly in the axial direction. This shape is similar to a longitudinal acoustic wave pattern. In addition, the same peak is obtained for a forcing on the axial momentum in the mixing duct as shown in Figure 10(b), but not for a forcing in the azimuthal momentum or when the inlet is non-reflecting. These elements indicate that the forcing that maximizes the amplification of fluctuations is the one that triggers the cavity modes of the mixing duct coupled with the combustion chamber. This means that the frequency of the optimal amplification mechanism is set by the acoustic mode that the coupled mixing tube and combustion chamber can sustain.

Figure 11.

Forcing on the axial momentum equation in the combustion chamber with (upper halves) closed inlet or (lower halves) non-reflecting inlet: real parts of (a) the forcing ${\hat{f}}_x$ , response (b) ${\hat{u}}_x$ , (c) ${\hat{u}}_{\theta }$ , and (d) $\hat{\omega }$ at 500 Hz. Figure (a) shows a zoom on the combustion chamber and Figure (b–d) on the mixing duct and the first half of the combustion chamber.

When forced in the mixing duct, the results can be compared to existing results on resolvent analysis of jets. The optimal forcings as shown in Figure 12(a) are dominated by tilted stripes in the direction of the shear, typical of the Orr mechanism, as shown by Schmidt et al. ²⁹ This mechanism is dominant in jets for Strouhal number $\mathrm{St}\in [0.3,1]$ , ³⁰ which translates to a frequency range of 190–630 Hz in the present case, based on the inlet velocity and the mixing duct hydraulic diameter. The optimal responses, as shown in Figures 11(b) and 12(b)to (d), are dominated by the Kelvin-Helmholtz mechanism, that is, wave packets shed at the mixing tube outlet with phase change cross the neutral line, ³¹ and the Orr mechanism, that is, tilted wave packets emerging from within the jet core and rotating as they propagate, until being perpendicular to the main stream. ²⁹ This rotation of the wave packets is known as an amplification mechanism during which the perturbations take energy from the mean flow. ³² Those qualitative observations are recovered in the various amplification regimes identified by resolvent analysis against spectral proper orthogonal decomposition (SPOD) ² of high-fidelity LES by Pickering et al. ³⁰

Figure 12.

Optimal amplification at 454 Hz, for forcing in the MT: forcing on (upper halves) the azimuthal momentum or (lower halves) the axial momentum: (a) optimal forcing, optimal responses, (b) $u_x$ , (c) $u_y$ , and (d) $u_{\theta }$ . Zoom on the mixing duct and the first half of the combustion chamber.

Although a modal resonance has been identified in the gain, the optimal response cannot simply be associated with the least stable eigenvector—that is, $s_{\text{A}}$ or $s_{\text{B}}$ —driving that resonance since the gain results from the interplay between several eigenmodes. However, inspecting the leading eigenmode can inform us about a part of the mechanism underlying the optimal amplification. The least stable eigenvalue, and, therefore, the peak gain frequency, changes from $s_{\text{A}}$ to $s_{\text{B}}$ when the mixing duct is shortened or the combustion chamber extended as shown in Figure 8.

In the three geometrical configurations investigated, the least stable eigenvalue is the one showing the strongest acoustic-like fluctuation in the mixing duct and the strongest vorticity shed from the edges of the mixing duct. For the base case, the vorticity of $s_{\text{A}}$ is mainly generated at the edges of the mixing duct and amplified downstream in the combustion chamber as shown in Figure 13(a) (upper half), whereas the vorticity of $s_{\text{B}}$ is primarily generated within the mixing duct as shown in Figure13(a) (lower half). The axial velocity perturbation of $s_{\text{A}}$ is dominated by the approximate quarter-wave acoustic mode as shown in Figure 14(a) (upper half), whereas it is dominated by convective-like patterns for $s_{\text{B}}$ as shown in Figure 14(a) (lower half). In the case of the shorter duct, for which $s_{\text{B}}$ dominates, these observations are inverted. The axial velocity response of $s_{\text{A}}$ is a superposition of acoustic and convective-like patterns, while $s_{\text{B}}$ is dominated by a quasi quarter-wave mode as shown in Figure 14(b). The vorticity of $s_{\text{A}}$ is generated within the mixing duct while it is shed from the edges of the duct for $s_{\text{B}}$ . The observation of the quasi quarter-wave mode is supported by the results of the network model, which show similar behavior of the velocity field. The above-mentioned observation indicates that the least stable eigenmode, hence driving the modal resonance, is the one showing the strongest Kelvin–Helmholtz mechanism and that this mechanism is coupled with a quasi quarter-wave mode sustained in the mixing duct. This mechanism does not describe the entire optimal response but is responsible for the observed frequency selection.

Figure 13.

Real part of the $\hat{\omega }$ component of the eigenvector $s_{\text{A}}$ (upper halves) and $s_{\text{B}}$ (lower halves) for (a) the base case and (b) $l$ _duct $=126$  mm (see Figure 10(b)). Zoom on the end of the mixing duct and first half of the combustion chamber.

Figure 14.

Real part of the ${\hat{u}}_x$ component of the eigenvector $s_{\text{A}}$ (upper halves) and $s_{\text{B}}$ (lower halves) for (a) the base case and (b) $l$ _duct $=126$  mm. Zoom on the mixing duct and first half of the combustion chamber.

Mode-conversion process in the optimal amplification

The optimal response to a forcing in the combustion chamber displayed in Figure 11 shows an acoustic-like wave pattern in the mixing duct and a vorticity pattern close to those of the leading eigenmode. Such an observation is consistent with the modal amplification identified on the gain curves in Figure 10. This vorticity is amplified while propagating downstream by the Kelvin–Helmholtz mechanism and eventually reaches the position where a flame can be stabilized. However, the optimal response is not similar to the least stable eigenmode since it includes azimuthal velocity fluctuations generated in the swirler section by velocity fluctuations of the upstream acoustic wave traveling in the mixing tube as shown in Figure 11(b). These azimuthal fluctuations are also generated in the case of a non-reflecting inlet in Figure 11 (lower halves), that is, in a purely anechoic setup. This proves that their generation is also possible without reflection at the inlet boundary as an actual conversion of upstream propagating acoustic waves originating from the combustion chamber.

In the purely anechoic case, only inertial waves can close the feedback loop by propagating back to the combustion chamber. The present results illustrate that the “mode conversion” process is at play in the optimal amplification mechanism. The mode conversion refers to the process that generates swirl fluctuations from acoustic waves impinging on a swirler. ⁶ Such a conversion is based on the interactions between sound and vorticity in the linear regime close to walls described by Chu and Kovasznay, ³⁴ and has often been modeled via the “actuator disk” theory. ³⁵ This process was already invoked in experiments and CFD simulations to explain the generation of swirl fluctuation, and it is now found to be part of optimal amplification. Additionally, it is not visible if not in the frequency range of the peak or for suboptimal gain, which makes the mode conversion process peculiar to the most efficient amplification mechanism.

The response of the swirler is non-negligible for all types of forcing. In addition, the optimal forcing in the swirler section of the mixing duct is non-negligible. It is even principally clustered in the swirler section for forcing on the azimuthal momentum, as shown in Figure 12(a). Firstly, this reveals that the swirler section is a region of interest with respect to the optimal forcing, and confirms the necessity of accounting for the knowledge of the flow close to the swirler suggested by Lückoff et al. ³⁶ and Müllet et al. ¹⁴ The reversibility of the mode conversion process is illustrated by the response to the forcing on the azimuthal momentum, as shown in Figure 12 (upper halves). In that case, the forced azimuthal momentum is preferentially converted into axial velocity fluctuations instead of generating azimuthal velocity fluctuations. Indeed, some tilted azimuthal velocity fluctuations are generated in the swirler but tend to vanish at the swirler exit. In contrast, axial velocity fluctuations are advected downstream and re-generate azimuthal velocity fluctuations further downstream, with wave packets perpendicular to the main stream. However, in the present case, these axial velocity fluctuations correspond to vorticity perturbations and not acoustic modes.

In the case of forcing on the axial momentum, as shown in Figure 12 (lower halves), the optimal forcing is a superposition of tilted convective structures with an acoustic structure corresponding to a quarter wavelength. The response is similar to the case of forcing on the azimuthal momentum, except that all three components of velocity fluctuations are produced within the swirler and amplified while propagating downstream. These two cases, that is, forcing in the mixing tube, illustrate that the mode conversion process is at play in all the optimal amplification mechanisms investigated.

Discussion of the results

The results of the resolvent analysis show that the optimal forcing in the combustion chamber corresponds to the pattern of an acoustic cavity mode, as shown in Figure 11(a). An unsteady heat source in the combustion chamber is very likely to produce such an acoustic field, and therefore to trigger the optimal amplification mechanism. This mechanism consists of a quarter-wave mode in the mixing duct coupled with the Kelvin–Helmholtz mechanism at the mixing duct edges and a mode conversion in the swirler producing swirl fluctuations.

The role of swirl fluctuations accounts here through their axial component. Axial velocity fluctuations stemming from both swirl and acoustic waves trigger vorticity shedding at the edges of the duct, drawing two paths for the feedback. This vorticity is then amplified by the Kelvin–Helmholtz mechanism. In the case of a flame stabilized in the combustion chamber, the vorticity is advected to the flame brush and increases or decreases the flame surface, which in turn modulates the heat release rate and produces some acoustic forcing on the axial momentum. The resulting feedback loop is depicted in Figure 15. The axial component of inertial waves was analytically evidenced by Albayrak et al. ⁴ and is observed in the present work on a realistic flow.

Figure 15.

Flame-flow feedback loop.

The existence of two paths for vorticity generation implies the existence of two time delays in the flame-flow interaction, and therefore a higher propensity for unstable feedback to arise. The time delay of the acoustic wave is negligible, whereas the time delay of the inertial wave depends on its propagation velocity and the position of the swirler. These two time delays also control whether the contributions of inertial and acoustic waves to the total axial velocity fluctuation are in phase, which is necessary to maximize the velocity at the mixing duct exit. This is a reminder that the propagation speed of inertial waves is a crucial parameter in the feedback loop. However, as explained in the second section, none of the fast or slow modes was clearly dominant in the present study. Nevertheless, the radial non-uniformity of the convection speed is expected to have a more substantial impact than the dominance of the fast or slow mode. Indeed, inertial wave structures are tilted by the radial non-uniformity of the mean flow. This change in the orientation of the convective structures impacts the constructive or destructive nature of the interplay between acoustic and inertial waves.

This optimal amplification mechanism has been shown to set off a frequency selection relying on acoustic resonance in the coupled mixing duct and combustion chamber. Since the ratio of the area jump in between these two elements is not large enough to fully decouple the acoustic modes existing in each of them, it is not possible to completely characterize these modes. However, the acoustic network model of the present configuration shows that both of the eigenmodes $s_{\text{A}}$ and $s_{\text{B}}$ could be approximated by a quarter-wave mode in the mixing duct, with a velocity antinode at the area jump. The presence of this velocity antinode at the exit of the mixing duct is a key feature since it produces a maximum velocity amplitude that triggers vortex shedding at the edges. This categorization depends on the acoustic behaviour at the inlet: for example, a constriction due to an injection lance or a narrow passage is likely to produce such a quarter-wave mode. On the other hand, an expansion into a plenum directly connected to the swirler inlet would result in an additional velocity antinode at the inlet and change the mode shape to a constant (0th order, or bulk mode) or any multiple of a half wavelength. However, the fully anechoic case as shown in Figure 11 (lower halves) presents the same feature as the closed inlet case, and in particular, the same dominance of the mode conversion process and vorticity shed at the mixing duct exit. Thereby, it is expected that the above-described mechanism remains active and dominant for various impedances of the inlet.

Conclusion and perspectives

The present work has shown the existence of inertial waves in a swirling flow representative of a swirl burner, that is, compressible, non-parallel, and turbulent. The mode conversion process and inertial waves are key features of the most efficient amplification mechanisms. In the context of combustion, unsteady heat release rate is likely to trigger such an optimal amplification mechanism involving inertial waves. The flow under study is computed via an axisymmetrical RANS model, including a modeled swirler. The capability of such an axisymmetrical modeled swirler to reproduce the static and dynamic behaviour of an actual swirler is verified against a 3D LES. The non-uniformity of the flow makes it difficult to conclude on the dominance of the “fast” or “slow” mode of the inertial waves, as discussed by Albayrak et al. ⁴ and Lückoff et al. ³⁶ However, the analytical mode shapes derived by Albayrak et al. ⁴ on an inviscid uniform flow agree with the azimuthal velocity fluctuations observed in the flow under investigation in the present work. Optimal amplification mechanisms are identified thanks to a resolvent analysis, including the modeled swirler. The results of the resolvent analysis first show that a significant part of the optimal amplification mechanism takes place in the swirler region, if not the optimal forcing itself. Hence, the inclusion of the swirler in the analysis is critical. The optimal amplification mechanism presents a frequency selection driven by a quarter-wave acoustic mode in the mixing duct. Two acoustic modes are competing in that case, and the dominant one is the mode enhancing the strongest streamwise Kelvin–Helmholtz mechanism at the mixing duct exit.

The conclusions of the present analysis can be further verified on a reactive LES, in particular, whether an unsteady heat source in the combustion chamber triggers the identified optimal amplification mechanism. The swirl intensity is also expected to significantly impact the results and should be investigated.