2.1 Compressible stability equations

The viscous, compressible equations and the equation of state are formulated in cylindrical polar coordinates ( $r,\unicode[STIX]{x1D703},z$ ), non-dimensionalized for a CRHPF of mean density $\bar{\unicode[STIX]{x1D70C}}_{c}$ , mean temperature $\bar{T}_{c}$ and mean axial velocity $\bar{u}_{c}$ inside a pipe of radius $R$ , rotating at a constant angular speed of $\unicode[STIX]{x1D6FA}$ , with $(\ast )_{c}$ denoting the respective centreline quantities, to yield

where $\text{D}/\text{D}t$ , $\unicode[STIX]{x1D6E5}$ , $L$ and $\unicode[STIX]{x1D6F7}$ are respectively the material derivative, the Laplacian operator, the dilatation term and the scalar dissipation function, given by

with $\unicode[STIX]{x1D6FE}=1.4$ being the ratio of specific heats, $Re=\bar{u}_{c}R/\unicode[STIX]{x1D708}$ the Reynolds number, $Pr=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D6FC}$ the Prandtl number and $Ma=\bar{u}_{c}/\bar{a}_{c}$ the Mach number, where $\unicode[STIX]{x1D708}$ and $\unicode[STIX]{x1D6FC}$ are, respectively, the kinematic viscosity and thermal diffusivity of the mean flow, taken here as constants, and $\bar{a}_{c}$ is the mean acoustic speed.

The linear stability equations are obtained from (2.1) via separating the flow variable $\boldsymbol{q}=[\unicode[STIX]{x1D70C}\,u_{r}\,u_{\unicode[STIX]{x1D703}}\,u_{z}\,p\,T]^{\text{T}}$ into a mean $\bar{\boldsymbol{q}}$ and a fluctuating $\boldsymbol{q}^{\prime }$ , with the latter modelled to possess travelling-wave-like solutions along the axial $z$ and azimuthal $\unicode[STIX]{x1D703}$ directions with a periodic time $t$ , by defining

where $\hat{\boldsymbol{q}}(r)$ is the unknown complex eigenfunction, $\unicode[STIX]{x1D6FC}$ and $m$ are respectively the axial and azimuthal wavenumbers and $\unicode[STIX]{x1D714}$ is the frequency. Standard analytical results are used for the mean $\bar{\boldsymbol{q}}$ , described in § 2.2. Using (2.3) in (2.1a )–(2.1f ) and linearizing for small fluctuations yields the compressible stability equations, written here as

where $\hat{\boldsymbol{q}_{1}}=[\hat{\unicode[STIX]{x1D70C}}\,\hat{u} _{r}\,\hat{u} _{\unicode[STIX]{x1D703}}\,\hat{u} _{z}\,\hat{p}]^{\text{T}}$ and $\unicode[STIX]{x1D647}=\unicode[STIX]{x1D647}_{1}+\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D647}_{2}+\unicode[STIX]{x1D647}_{3}/Re+\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D647}_{4}/Re+\unicode[STIX]{x1D6FC}^{2}\unicode[STIX]{x1D647}_{5}/Re$ . Each of the operators $\unicode[STIX]{x1D647}_{i}$ , $i=1{-}5$ , is a $5\times 5$ matrix, whose details are given in appendix B.

2.2 Mean flow and boundary conditions

In general, both $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D714}$ in (2.3) can be complex numbers in a spatio-temporal analysis. In this work, we track the convective instability boundary via a spatial analysis which imposes $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FC}_{r}+\text{i}\unicode[STIX]{x1D6FC}_{i}$ , where $\unicode[STIX]{x1D6FC}_{r}$ is the streamwise wavenumber and $\unicode[STIX]{x1D6FC}_{i}$ is the growth rate, while $\unicode[STIX]{x1D714}$ and $m$ are now specified real numbers. The mean $\bar{\boldsymbol{q}}_{1}$ is specified as

where $\unicode[STIX]{x1D716}$ is the Rossby number defined as $\unicode[STIX]{x1D716}=\bar{u}_{c}/\unicode[STIX]{x1D6FA}R$ , the inverse of the swirl number, used to specify $\bar{u}_{\unicode[STIX]{x1D703}}$ , while $\bar{u}_{z}$ is from standard results of laminar pipe flow. The mean density $\bar{\unicode[STIX]{x1D70C}}(r)$ is specified, as discussed next, which yields via (2.1b ) a mean pressure

where $\bar{p}_{0}=1/\unicode[STIX]{x1D6FE}Ma^{2}$ , the mean pressure at the pipe centreline. Note here that the Rossby number may also be included to define an azimuthal Reynolds number $Re_{\unicode[STIX]{x1D703}}=\bar{u}_{c}R/\unicode[STIX]{x1D716}\unicode[STIX]{x1D708}$ , which we use extensively in this work.

Table 1.

Mean density cases considered.

Figure 1.

Mean profiles with radius-dependent density $\bar{\unicode[STIX]{x1D70C}}(r)$ , showing —— case B, – – – case C and – ⋅ – ⋅ – case C1 of table 1. The corresponding grey curves are the gradients $\bar{\unicode[STIX]{x1D70C}}^{\prime }(r)$ .

Studies of inviscid, incompressible swirling flows by Fung & Kurzweg (Reference Fung and Kurzweg1975) with monotonically increasing $\bar{\unicode[STIX]{x1D70C}}(r)$ and $\bar{\unicode[STIX]{x1D70C}}^{\prime }(r)$ have shown unconditional stability at all disturbance frequencies, while a density stratification that is radially decreasing is likely to induce instabilities. Physically speaking, this is unlikely to be straightforward for compressible viscous flows with more options for energy mechanisms available (see (2.11)), where the nature of density stratification along with the corresponding gradients are expected to be equally important in deciding the overall stability, besides the effects of core compressibility via the chosen Mach number. To assess the role of stratification, four types of mean density profile are considered in this work, summarized in table 1 and figure 1. The motivation behind such choices comes from the air-assisted and air blast atomizers of the modern gas turbine fuel injection systems, introduced in § 1, which have surprisingly wide design variations that in most cases amount to multiple swirling jets of fuel and air introduced at different locations (often at different temperatures) to better control the mixing between fuel and oxidizer ahead of the combustion zone. Case A in table 1, a constant mean density case, models uniform (or perfect) mixing between fuel and air. This case, in its compressible and incompressible forms (of Fernandez-Feria & del Pino Reference Fernandez-Feria and del Pino2002), also serves as a benchmark for the other compressible cases. These cases introduce different types of heterogeneities in $\bar{\unicode[STIX]{x1D70C}}(r)$ , likely when mixing remains incomplete, of which cases B and C of table 1, having higher densities near the centre, may be the result of configurations that introduce greater fractions of the usually higher-density air near the swirler core. Conversely, case C1 of table 1, similar to the monotonically increasing density case of Fung & Kurzweg (Reference Fung and Kurzweg1975), may be due to a fuel bias near the central region. Additionally, we use different mathematical models for the decreasing density cases, where case B has an exponential variation and case C models algebraic changes. The choice of such functions introduces differences in density gradients $\bar{\unicode[STIX]{x1D70C}}^{\prime }(r)$ , where the algebraic case has a monotonically increasing density gradient while the gradient in the exponential case decreases monotonically until $r=0.7$ , beyond which it rises a little (see figure 1). As noted before, in spite of similarly decreasing densities, differing density gradients are expected to yield dissimilar effects on the flow stability due to the presence of these gradient terms in (2.11). Finally, note that in practical configurations, the density stratification is perhaps a combination of all or some of these cases, but as we shall see in § 3, a better understanding of the associated energy mechanisms is obtained by considering them individually.

The boundary conditions are specified at the pipe axis $r=0$ (see e.g. Batchelor & Gill Reference Batchelor and Gill1962; Khorrami, Malik & Ash Reference Khorrami, Malik and Ash1989):

and at the pipe wall $r=1$ :

where $\unicode[STIX]{x1D712}_{1}$ and $\unicode[STIX]{x1D712}_{2}$ are constants, set here to zero, which has been found to yield no loss of accuracy while obtaining numerical solutions (see also Khorrami et al. Reference Khorrami, Malik and Ash1989).

2.3 Disturbance energy equation

Unlike incompressible flows, where a linearized vorticity equation may provide important clues on instability mechanisms, a total energy-based formulation is required to paint a more complete picture for compressible flows like the CRHPF considered here. In this section, we develop the total disturbance energy equation following the lines of e.g. Chu (Reference Chu1965) and Lu & Lele (Reference Lu and Lele1999) to obtain the energy budget of the different instability mechanisms yielding the observed stability states. In particular, as the instability transitions from a convective to absolutely unstable state or from a stable to convectively unstable state, the focus would be on understanding the varying roles of the component energy mechanisms. It may be noted here that the incompressible version of this equation is the classical Reynolds–Orr equation (see e.g. Wang et al. Reference Wang, Rusak, Gong and Liu2016).

The kinetic energy $k$ represents the energy from the vortical modes via

which satisfies the kinetic energy equation that may be formed by combining (2.1b )–(2.1d ), while the total disturbance energy (e.g. Chu Reference Chu1965)

also includes contributions from the acoustic and entropy modes, governed directly by (2.1a ) and (2.1e ), respectively. Here, the entropy $s=p/\unicode[STIX]{x1D6FE}\bar{p}-\unicode[STIX]{x1D70C}/\bar{\unicode[STIX]{x1D70C}}$ is non-dimensionalized by $c_{p}$ , the specific heat capacity at constant pressure. The total disturbance energy equation for a given mean flow can now be formed from (2.1a )–(2.1f ), which on using (2.10) yields

where

with $()^{\ast }$ denoting the complex conjugate operation and c.c. representing the sum of all complex conjugate terms on the right of (2.11). Written in this way, (2.11) is valid for any compressible viscous flow with the mean $\bar{u}_{r}=0$ (see (2.5)), true for most swirling mean flow models. The different components of (2.11), as labelled, are as follows: $E_{u_{z}}$ is the production term due to shear of axial velocity $u_{z}$ , $E_{u_{\unicode[STIX]{x1D703}}}$ is due to shear of azimuthal velocity $u_{\unicode[STIX]{x1D703}}$ , $E_{p}$ is a term that redistributes energy from near the source of unstable energy to the interior of the pipe, $E_{\unicode[STIX]{x1D708}}$ accounts for the momentum dissipation from viscous terms, $E_{\unicode[STIX]{x1D6F7}}$ accounts for the dissipation from thermal conduction and $\unicode[STIX]{x1D6F7}$ of (2.1e ) and finally $E_{s}$ is the entropy fluctuation term, including the inviscid entropy fluctuation $E_{si}$ , a non-zero term for all swirling flows irrespective of the density stratification; $E_{u_{z}}$ is the predominant source of instability for incompressible swirling flows and at higher $\unicode[STIX]{x1D716}$ but is of diminishing importance, as we shall see, at higher $Re_{\unicode[STIX]{x1D703}}$ . A rigid-body rotation (see (2.5)) rules out $E_{u_{\unicode[STIX]{x1D703}}}$ , while $E_{p}$ usually promotes instability by moving energy from near the wall, if available, toward the centre of the pipe. The entropy fluctuations $E_{s}$ , which in the way defined in (2.11) directly include the density stratification and the pressure–density coupling via $E_{si}$ , are expected to be important, while the viscous part $E_{sv}$ may turn out to be a major stabilizing mechanism. The viscous dissipation term $E_{\unicode[STIX]{x1D708}}$ always acts as a sink of energy, resisting unstable growth and thus providing the traditional stabilizing role. The contribution from $E_{\unicode[STIX]{x1D6F7}}$ is usually very small for all the cases considered and will not be shown separately. Note that in (2.9)–(2.11), we have removed the primes from fluctuations for clarity. While computing (2.11), it is worth remembering that as the fluctuations are modelled via (2.3), the sign of the total energy $E_{total}(r)$ is automatically fixed via the sign of $-\unicode[STIX]{x1D6FC}_{i}$ , so that for a spatially unstable mode $(-\unicode[STIX]{x1D6FC}_{i}>0)$ , the total energy is positive at all radial locations and conversely is negative for a stable mode. Importantly, in § 3, the components of equation (2.11) are reported after normalizing each case with the corresponding maximum total energy $E_{total}(r)|_{max}$ at $r=r_{max}$ .

2.4 Numerical solution of stability equations

Equation (2.4) is numerically solved using a standard Chebyshev spectral collocation technique before transforming to a linear companion matrix (e.g. Bridges & Morris Reference Bridges and Morris1984) to yield

where $\hat{\boldsymbol{q}}_{2}=\unicode[STIX]{x1D6FC}\hat{\boldsymbol{q}}_{1}$ and $\unicode[STIX]{x1D644}$ is the $5\times 5$ identity matrix. The $-1\leqslant s\leqslant 1$ interval of the Chebyshev polynomials is mapped to a $0\leqslant r\leqslant 1$ interval, as required here, via the transformation $s=2r-1$ . Equation (2.13), discretized at the $N$ Gauss–Lobatto points, is solved along with the boundary conditions (2.7) and (2.8) using standard eigen solvers. Spurious eigenvalues are discarded by setting a convergence criterion for all the eigenvalues. A value of $N=50$ was found to be sufficient for the convergence of all physical eigenvalues up to at least seven decimal places.

Note that in (2.3), since $\unicode[STIX]{x1D6FC}(m,\unicode[STIX]{x1D714})=-\unicode[STIX]{x1D6FC}^{\ast }(-m,-\unicode[STIX]{x1D714})$ , we follow a convention where we set $m\leqslant 0$ and allow $\unicode[STIX]{x1D714}$ to take any sign, so that as $m<0$ and $\unicode[STIX]{x1D714}<0$ , (2.3) yields

where for a spatially unstable mode we still require $-\unicode[STIX]{x1D6FC}_{i}>0$ . The sign of the corresponding phase speed $c_{p}=\unicode[STIX]{x1D714}/\unicode[STIX]{x1D6FC}_{r}$ is decided by the relative signs of $\unicode[STIX]{x1D714}$ and $\unicode[STIX]{x1D6FC}_{r}$ . Note here that in this work we only consider the azimuthal modes $m=-1,-2$ , as the axisymmetric $m=0$ mode is known to be convectively stable for all rotations (see e.g. Pedley Reference Pedley1968, Reference Pedley1969), while the same is also checked to be true for our compressible cases with imposed density stratifications.

Figure 2.

The onset of absolute instability for the most unstable mode in case A of table 1 at $\unicode[STIX]{x1D716}=0.85$ and $Ma=0.8$ is tracked by monitoring (a) the spatial growth rate $\unicode[STIX]{x1D6FC}_{i}$ and (b) the real group velocity $c_{r}$ as a function of the disturbance frequency $\unicode[STIX]{x1D714}$ for successive values of $Re$ (– – –  $Re=217$ , – ⋅ – ⋅ –  $Re=222$ and ——  $Re=Re_{t}=227$ ) until a pinch point is found at $\unicode[STIX]{x1D714}=-0.971$ for the last Reynolds number.

2.5 Tracking the onset of absolute instability

In the Briggs–Bers formalism (e.g. Briggs Reference Briggs1964; Bers Reference Bers1983; Huerre & Rossi Reference Huerre and Rossi1998; Schmid & Henningson Reference Schmid and Henningson2001), absolute instability of a mode is tracked via its spatio-temporal behaviour, where the occurrence of pinch points in the complex $\unicode[STIX]{x1D6FC}$ -plane, together with branch points in a specific half of the complex $\unicode[STIX]{x1D714}$ -plane, renders the mode absolutely unstable. However, if the details of the absolute modal growth rate are not sought, a simpler spatial stability-based method, as described here (following e.g. Fernandez-Feria & del Pino Reference Fernandez-Feria and del Pino2002), may be used to predict the onset of absolute instability.

Here, the complex group velocity

is zero when an unstable pinch/branch point is reached (see Schmid & Henningson Reference Schmid and Henningson2001), which is an alternative to solving the corresponding complex dispersion relation. In this process, the real group velocity $c_{r}(\unicode[STIX]{x1D714})$ and $\unicode[STIX]{x1D6FC}_{i}(\unicode[STIX]{x1D714})$ are tracked by varying $Re$ so that at $Re=Re_{t}$ when the absolute instability is reached, $c=0$ , which via (2.15) yields $c_{r}=0$ and $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}_{i}/\unicode[STIX]{x2202}\unicode[STIX]{x1D714}\rightarrow \infty$ , the basic requirements for pinching. Figure 2, for example, demonstrates this for the most unstable mode in case A of table 1, where with $Re$ gradually increased for a given $\unicode[STIX]{x1D716}$ until $Re=Re_{t}=227$ , a pinching appears in figure 2(a) and simultaneously $c_{r}=0$ (figure 2 b) at $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{t}=-0.971$ . In addition, it is checked (not shown here) that the pinch point in the complex $\unicode[STIX]{x1D6FC}$ -plane is formed by two spatial branches $\unicode[STIX]{x1D6FC}^{+}$ and $\unicode[STIX]{x1D6FC}^{-}$ such that the negative branch always lies in the lower half-plane $(\unicode[STIX]{x1D6FC}_{i}<0)$ , while the positive branch moves to the upper half-plane as $\unicode[STIX]{x1D714}_{i}$ , the imaginary part of a complex $\unicode[STIX]{x1D714}$ , is increased from zero. This process is repeated for all $\unicode[STIX]{x1D716}$ , yielding the stability boundaries between the respective convective and absolutely unstable states.

Figure 3.

Solutions of (2.13) for $Re=100$ , $\unicode[STIX]{x1D716}=0.5$ , $Ma=0.8$ , $m=-1$ , $\unicode[STIX]{x1D714}=-1$ and $\overline{\unicode[STIX]{x1D70C}}=1$ showing (a) the eigenvalue spectrum, with (●) being the marginally unstable mode, and (b) absolute magnitude of the eigenfunctions with ——  $\tilde{u} _{z}$ , – – –  $\tilde{u} _{\unicode[STIX]{x1D703}}$ , –  $\cdot$  –  $\cdot$   $\tilde{u} _{r}$ and –  $\cdot \cdot$  –  $\cdot \cdot$   $\tilde{p}$ shown.

3.1 Eigenspectra and mode shapes

Figure 3(a) shows the eigenspectrum for parameters satisfying case A of table 1 with $Re=100$ , $\unicode[STIX]{x1D716}=0.5$ and $Ma=0.8$ . Only the physical eigenvalues with $\unicode[STIX]{x1D6FC}_{r}>0$ are shown, while unphysical solutions that appear naturally in Chebyshev spectral collocation methods may be arbitrarily shifted away from the physical ones by a capacitance matrix method (e.g. Hagan & Priede Reference Hagan and Priede2013), if needed. The flow of figure 3(a) is only marginally unstable, since it has a mode with a small negative $\unicode[STIX]{x1D6FC}_{i}$ . The corresponding mode shapes are in figure 3(b), where the observed variation in $\tilde{u} _{r}$ may couple with a favourable mean density gradient to generate strong entropy waves (see inviscid $E_{s}$ in (2.11)) or it may act in tandem with $\tilde{p}$ to cause redistribution of shear energy within the pipe flow (see $E_{p}$ in (2.11)); $\tilde{u} _{\unicode[STIX]{x1D703}}$ rises close to the wall, generating a fluctuating centrifugal force characteristic of swirling flows. The fluctuating $\tilde{u} _{z}$ , $\tilde{\unicode[STIX]{x1D70C}}$ (not shown) and $\tilde{p}$ all attain their respective peaks inside the core region of the pipe flow, away from the wall.

Figure 4.

Neutral curves showing regions of convective stability, convective instability and absolute instability for ——  $m=-1$ and – – –  $m=-2$ perturbations for (a,b) constant mean density (case A), (c) exponentially varying mean density (case B) and (d) algebraically varying mean density (case C) (for cases see table 1), all at $Ma=0.8$ with the grey curves showing the incompressible case $(Ma\rightarrow 0)$ . In (a,b), the set points described in table 2 are marked with black dots, while labels against the dotted lines in (b–d) indicate the respective limiting values of $Re$ and $Re_{\unicode[STIX]{x1D703}}$ . In all figures, the thicker curves denote the convective–absolute boundary.

A detailed understanding of the flow instability mechanisms is incomplete unless the fluctuating quantities are considered in conjunction with the mean flow parameters, which is precisely what the disturbance energy formulation (2.11) aims at, and this will be our focus in § 3.3.

3.2 Effect of density stratification and compressibility on neutral curves

Neutral curves marking the onset of convective and absolute instabilities for the first and second asymmetric azimuthal modes $m=-1,-2$ at a fixed Mach number of $Ma=0.8$ are shown in figure 4 for the mean density cases A, B and C of table 1, with the incompressible case of Fernandez-Feria & del Pino (Reference Fernandez-Feria and del Pino2002) plotted alongside for reference. Note here that the unconditionally stable case C1 does not possess a neutral curve. Further, for swirling flows, the neutral curves may be plotted in at least two different ways: as functions of $Re$ versus $\unicode[STIX]{x1D716}$ or $Re$ versus $Re_{\unicode[STIX]{x1D703}}$ , both of which include $\bar{u}_{c}$ and $\unicode[STIX]{x1D6FA}$ (the latter via $\unicode[STIX]{x1D716}$ ), the two mean quantities known to act as potential sources of instability at the incompressible limit. However, as we shall see, once compressibility is introduced, mean pressure (and density) dominates at certain parametric configurations, especially at lower $\unicode[STIX]{x1D716}$ , when it is more prudent to use $Re_{\unicode[STIX]{x1D703}}$ as the independent variable since it better describes the different limiting configurations.

Figure 4(a,b) show neutral curves for case A, when the mean density is constant and $Ma=0.8$ . As the rotational speed $\unicode[STIX]{x1D6FA}$ is gradually raised for this Mach number, the $m=-1,-2$ neutral curves for convective instability follow the corresponding incompressible curves until $\unicode[STIX]{x1D716}\approx 2$ , before breaking away (see figure 4 a). In figure 4(b), at lower $Re\lessapprox 45$ , the convectively unstable neutral curves asymptotically approach $Re_{\unicode[STIX]{x1D703}}=Re_{\unicode[STIX]{x1D703}c}\simeq 80.3$ for $m=-1$ and $Re_{\unicode[STIX]{x1D703}c}\simeq 75.0$ for $m=-2$ (labelled in figure 4 b), cutoff values below which the flow is convectively stable at these conditions. Meanwhile, the incompressible curves follow a completely different path as $\unicode[STIX]{x1D716}$ is lowered, instead asymptotically approaching $Re=Re_{c}\simeq 82.9$ for $m=-1$ and $Re_{c}\simeq 91.0$ for $m=-2$ (also labelled in figure 4(b), confirming previous findings of e.g. Mackrodt Reference Mackrodt1976; Cotton & Salwen Reference Cotton and Salwen1981; Fernandez-Feria & del Pino Reference Fernandez-Feria and del Pino2002), respectively, below which unconditionally stable states are reached. At the higher $Re\gtrapprox 200$ , both the constant-density compressible and incompressible curves approach $Re_{\unicode[STIX]{x1D703}}=Re_{\unicode[STIX]{x1D703}c}\simeq 27.0$ for $m=-1$ and $Re_{\unicode[STIX]{x1D703}c}\simeq 43.5$ for $m=-2$ , identical to the previous spatial and temporal analysis findings for incompressible conditions. The neutral curves marking the onset of absolute instability show similar trends. As discussed by Fernandez-Feria & del Pino (Reference Fernandez-Feria and del Pino2002), incompressible swirling pipe flow demonstrates a minimum $Re_{\unicode[STIX]{x1D703}}$ below which absolute instability is impossible for any $Re$ . Here, for $m=-1$ , this limit is found to be $Re_{\unicode[STIX]{x1D703}}=Re_{\unicode[STIX]{x1D703}t}\simeq 254.2$ , while for $m=-2$ it is at $Re_{\unicode[STIX]{x1D703}t}\simeq 326.9$ . In figure 4, the compressible curves are not seen to possess any such clear limits (although at lower Mach numbers there appears to be one; see especially figure 5 a), and instead at lower $Re\rightarrow 0$ these asymptotically approach the respective convective instability curves (as per their azimuthal orders), and thus a convectively stable compressible swirling pipe flow directly transitions to a state of absolute instability at this $Re_{\unicode[STIX]{x1D703}}$ ( $=Re_{\unicode[STIX]{x1D703}c}=Re_{\unicode[STIX]{x1D703}t}$ ) once $Re$ is sufficiently low, as described.

The neutral curves of the exponential mean density case B of table 1, shown in figure 4(c), are not too dissimilar to case A, except perhaps for an earlier break away from the incompressible curve for the $m=-1$ convective instability as $Re$ is lowered. The $m=-2$ curve for this case is shifted rightwards (see figure 4 c for the limiting $Re_{\unicode[STIX]{x1D703}}$ values), and this coupled with the shape of the $m=-1$ curve points to the little extra convective stability provided via such a density stratification at these higher $Re$ values as compared to case A. As $Re$ is lowered, the absolute instability curves along with the convective ones asymptotically approach $Re_{\unicode[STIX]{x1D703}}\simeq 90.3$ for $m=-1$ and $Re_{\unicode[STIX]{x1D703}}\simeq 83.6$ for $m=-2$ . Thus, as in case A, the higher azimuthal mode actually seems to be slightly less stable at these lower $Re$ numbers, and this effect increases dramatically for case C, as we discuss next.

Table 2.

Disturbance energy set points considered. Stability characteristic refers to cases A–C of table 1, except case C1, which is CS at all the set points.

Figure 5.

Same as in figure 4, the flow Mach number $Ma$ is varied for (a) the constant density case A ( $m=-1$ ) and (b) the algebraic density case C ( $m=-2$ ) of table 1, with –  $\cdot \,\cdot$  –  $\cdot \,\cdot$  –  $Ma=0.01$ , — — —  $Ma=0.2$ , –  $\cdot$  –  $\cdot$  –  $Ma=0.4$ , – – –  $Ma=0.6$ and ——  $Ma=0.8$ .

The stability of the algebraic mean density case C of table 1 is significantly different at higher $Re$ (see figure 4 d), which directly follows from the nature of this density stratification (see figure 1), where stronger density gradients at $r\rightarrow 0$ are unlike those of cases A and B. As we shall see in § 3.3.3, this brings about a strong change in the viscous entropy fluctuations even as $\unicode[STIX]{x1D716}\rightarrow \infty$ , yielding greater convective stability at higher $Re$ (check the limiting $Re_{\unicode[STIX]{x1D703}}$ values in figure 4 d), but especially at the higher azimuthal order of $m=-2$ . In contrast, as $Re\rightarrow 0$ (or at sufficiently low $Re$ ), the lower azimuthal-order perturbations (here, $m=-1$ ) are more convectively stable (i.e.  $Re_{\unicode[STIX]{x1D703}c}$ higher), as is clear from the figure.

Although as $Re\rightarrow \infty$ the convective stability is improved via stratification, especially in case C, the boundary of absolute instability remains unaltered and simply follows that of the incompressible case. Hence, density stratification plays no role at the convective–absolute boundary for sufficiently high $Re$ and, as we shall see next, neither does flow compressibility.

The surprising evolution of the neutral curves at lower $Re$ requires further investigation. Accordingly, in figure 5, we plot the same neutral curves but with compressibility now varied via changing the mean flow Mach number. Two cases from table 1 are shown: the constant density case A for $m=-1$ in figure 5(a) and the algebraic density case C for $m=-2$ in figure 5(b), since these are the furthest apart in terms of their respective variations with Reynolds number, as seen in figure 4. The results for the other cases, not shown, are essentially in between the two shown. As compressibility is gradually raised from the incompressible $Ma\rightarrow 0$ limit, neutral curves are seen to dramatically shift in such a way as to settle in a state increasingly governed by the azimuthal Reynolds number $Re_{\unicode[STIX]{x1D703}}$ over any constant $Re$ . For a given $Ma$ and $Re$ , as $Re_{\unicode[STIX]{x1D703}}$ is gradually raised, or as swirl rates go up, the mean pressure starts to dominate via (2.6). If the flow Mach number is low, say $Ma=0.01$ , the background (centreline) pressure $\bar{p}_{0}=1/\unicode[STIX]{x1D6FE}Ma^{2}$ can be a relatively large quantity, which then requires the Rossby number $\unicode[STIX]{x1D716}$ to be significantly low (i.e. very high $Re_{\unicode[STIX]{x1D703}}$ ) for the integral term in (2.6) to exceed $\bar{p}_{0}$ . In figure 5 this translates to convective instability curves following the incompressible curve up to a higher $Re_{\unicode[STIX]{x1D703}}$ for lower $Ma$ , but breaking away at much lower $Re_{\unicode[STIX]{x1D703}}$ as the compressible effects increase at higher $Ma$ . Now, if there is density stratification as in figure 5(b), the convective stability curves show the stability to be radically altered even if the flow is close to incompressible, where the corresponding $Ma=0.01$ curve evolves very differently to the $m=-2$ incompressible curve (without any stratification). In this case, the effect of the Mach number seems to start a little higher at $Re\simeq 205$ , and although increased compressibility via increasing $Ma$ still shrinks the region of convective stability, this effect is much less than in figure 5(a).

To summarize, the overall stability of compressible swirling pipe flow turns out to be quite interesting. Whereas for incompressible pipe flows, if the Reynolds number is below a critical value ( $Re_{c}\simeq 82.9$ for case A at $m=-1$ and $Re_{c}\simeq 91.0$ at $m=-2$ ) the flow is convectively stable for all rotation rates, increased compressibility seems to make it unstable, especially more absolutely unstable with fast disappearance of the convectively stable zones. The swirling pipe flow now almost immediately transitions to an absolutely unstable state from an initially convectively stable state as $Re_{\unicode[STIX]{x1D703}}$ is raised beyond a critical value, which is increasingly lower for higher $M$ . In this case, compressibility actually promotes the onset of instability, which is in direct contrast to prior findings on vortex flows, which were seen to be stabilized with increased compressibility (see e.g. Herrada et al. Reference Herrada, Pérez-Saborid and Barrero2003; Rusak et al. Reference Rusak, Choi and Lee2007, Reference Rusak, Choi, Bourquard and Wang2015). Further, if $Re$ is above this critical value, there seems to be hardly any effect of compressibility. Density stratification alters this, where a radially increasing mean density (case C1) can make the flow convectively stable at all Mach numbers, including when the flow is treated as incompressible. If this stratification nature is flipped, the stabilizing role of compressibility is now recovered, but only above a critical flow Reynolds number (e.g.  $Re\gtrapprox 205$ for case C at $m=-2$ ), and even then it depends on the nature of mean density stratification (e.g. the gradient of density) and the order of the azimuthal modes, with the higher modes becoming more stable. Finally, at these higher Reynolds numbers, for a convectively unstable flow, the onset of absolute instability is independent of compressibility and stratification levels.

In the next section, we discuss the reasons behind such observations, but before introducing the energy mechanisms it may be useful to reflect upon the $Re$ and $Ma$ ranges, as plotted in figures 4 and 5, for the purpose of trying to relate these to practical configurations. Flow Reynolds numbers in standard swirler designs are of the order of $10^{4}$ (see e.g. Wang et al. Reference Wang, Yang, Hsiao, Hsieh and Mongia2007), which could significantly alter depending upon the viscosities of the fuel–air mixture or in the case of compact combustors where swirler diameters are of the order of millimetres. Nevertheless, in figures 4 and 5(b,d) the neutral curves for convective instability are clearly independent of $Re$ at $Re>10^{3}$ , while the absolute instability curves follow a linear trend (in log scale) that is easily extrapolated. As far as the Mach numbers go, as discussed previously, for certain atomizer designs with air-assisted fuel breakup, this can be potentially high, perhaps reaching near-choking conditions near the nozzle exit (see e.g. Lefebvre & Ballal Reference Lefebvre and Ballal2010).

Figure 6.

Disturbance energies for case A (constant mean density) of table 1 for (a,c,e) $m=-1$ and (b,d,f) $m=-2$ , at the set points (a,b) SP1 (c,d) SP2 and (e,f) SP3 of table 2, showing the following components of (2.11): – ⋅ – ⋅ –  $E_{u_{z}}$ , – – –  $E_{\unicode[STIX]{x1D708}}$ , — — —  $E_{p}$ , ——  $E_{s}$ . The components of $E_{s}$ , –  $\bullet$  –  $\bullet$  –  $E_{si}$ and –  $\circ$  –  $\circ$  –  $E_{sv}$ , are only shown in some important cases. The thick solid grey line in each panel represents the total energy.

3.3 Energy mechanisms

In this section, instability mechanisms for the $m=-1$ and $m=-2$ modes are discussed with reference to the set points defined in table 2 and marked in figure 4(a,b), which represent all three stability states for cases A, B and C of table 1. The set point SP1 with a high Rossby number (or low $Re_{\unicode[STIX]{x1D703}}$ ) ensures all the cases of table 1 are convectively stable (CS). At SP2, all the cases except case C1 are convectively unstable (CU) and at SP3 these are absolutely unstable (AU). As we shall see in the following subsections, the way in which the energy budget due to (2.11) reaches identical stability states is quite distinct for the different stratified density cases. In what follows, for the CS cases (at SP1) the least stable modal solution is shown while unstable solutions refer to the most unstable mode. Table 3 summarizes the disturbance waves for cases A, B and C of table 1, shown next in figures 6–9.

Table 3.

Frequencies and wavenumbers of the perturbation waves calculated at the set points of table 2. Negative $\unicode[STIX]{x1D6FC}_{i}$ (growth rates) are unstable. Cases in parenthesis refer to those of table 1.