2.1. Burner and operating conditions

The burner used in the current study is shown in figure 1; its geometry has been described in previous publications (Geigle et al. Reference Geigle, Köhler, O'Loughlin and Meier2015, Reference Geigle, Hadef, Stöhr and Meier2017). For completeness, it is briefly described again here. The injector consists of a pair of annular swirl nozzles separated by a ring of 60 fuel channels, each with an area of $0.5\times 0.4\ {\rm mm}^{2}$, for the introduction of the gaseous fuel, ethylene. The central air nozzle has a diameter of 12.3 mm and a nominal geometric swirl number of 0.82 (Lefebvre Reference Lefebvre1998). The annular outer air nozzle measures 19.8 mm in diameter and has a nominal geometric swirl number of 0.79. The combustion chamber has a square cross-section measuring $68\times 68\ {\rm mm}^{2}$. It is 120 mm long and enclosed by 3 mm-thick quartz windows for optical access. The water-cooled dome of the combustion chamber has a cylindrical exhaust hole – diameter 40 mm, length 24 mm – linked to the combustion chamber by a curved tube. The combustor was mounted in an optically accessible pressure vessel.

Figure 1. Experimental setup.

The flows of air (central, annular and seeding) and fuel were supplied by separate mass flow controllers (Bronkhorst). The flow rate of air is 841.7 standard litres per minute (SLM) and flow rate of fuel is 37.9 SLM, resulting in a thermal power of 38.4 kW. In the non-reacting studies, no fuel is injected into the system but the air flow rate is the same. The flow controllers were calibrated in house, resulting in an accuracy of better than 1.0 % of the controllers’ maximum flow. The majority of the ethylene flow was bubbled through a temperature-controlled reservoir ($T=25\,^{\circ }{\rm C}$, 298 K) filled with acetone, which provided an acetone-saturated ethylene flow of the desired composition. For the present study, the combustor was operated with ethylene/acetone fuel at a pressure of 5 bars and $\varPhi = 0.67$. The air split, defined as the fraction of air flow passing through the centre nozzle, was varied from 0 to 1 in the non-reacting cases and 0.2 to 0.5 in the reacting cases. The Reynolds number based on the diameter of the outer nozzle is $Re = 19\,490$.

2.2. Measurement techniques

Time-resolved stereo particle image velocimetry (sPIV), planar laser-induced fluorescence of OH (OH-PLIF) and acetone-PLIF were applied simultaneously with a repetition rate of 10 kHz to capture flame/vortex/fuel interaction. Data were obtained for 1 s in the technically premixed ethylene/air swirl flame stabilized in the model of the aero-engine combustor. This diagnostic system has been described previously in the scientific literature (Litvinov et al. Reference Litvinov, Yoon, Noren, Stöhr, Boxx and Geigle2021). For completeness, a brief description is included below.

2.3. Data analysis

2.3.2. Complex network analysis

Coherent self-excited flow oscillations are generated due to internal feedback of flow oscillations within critical regions in the flow field. This feedback is generated by nonlinearity to leading order by the interaction of the flow oscillation with itself. These flow regions are referred to as wavemakers. They can be identified by regions of high structural sensitivity or receptivity of self-excited linear modes as determined from hydrodynamic stability analysis, as shown in several prior studies of self-excited bluff-body wakes (Giannetti & Luchini Reference Giannetti and Luchini2007; Marquet, Sipp & Jacquin Reference Marquet, Sipp and Jacquin2008; Juniper, Tammisola & Lundell Reference Juniper, Tammisola and Lundell2011) and swirl nozzles (Tammisola & Juniper Reference Tammisola and Juniper2016; Mukherjee et al. Reference Mukherjee, Muthichur, More, Gupta and Hemchandra2021). However, applying these approaches to analyse the present flow problem is challenging. Firstly, the highly turbulent and density stratified nature of the present flow field would introduce large quantitative uncertainties in a hydrodynamic stability analysis of the current flow since unambiguous quantitatively accurate density fields are not directly measured. Secondly, it is known that regions of receptivity for hydrodynamic modes are localized inside nozzles and close to burner lips where shear layers separate (Wang et al. Reference Wang, Kaiser, Meindl, Oberleithner, Polifke and Lesshafft2022). Determining time-averaged flow fields in these regions for the present nozzle using PIV is challenging. Also, optical access to determine time-averaged flows within the nozzle is not possible in the current experimental setup. Therefore, we adopt an alternate data-driven approach from complex network theory to determine the position and shape of critical regions associated with various flow field oscillations from time-series PIV data directly. These in turn yield insight into the position and spatial extent of wavemaker regions as will be described in this section.

Our approach closely follows the methods described in Krishnan et al. (Reference Krishnan, Manikandan, Midhun, Reeja, Unni, Sujith, Marwan and Kurths2019). They apply complex network analysis to identify the critical region driving shear layer instability in an axisymmetric bluff-body-stabilized turbulent combustor. They show experimentally that modifying the flow in this region stabilizes the shear layer and the flow instability. This result is analogous to the study of Mukherjee et al. (Reference Mukherjee, Muthichur, More, Gupta and Hemchandra2021), wherein traditional structural sensitivity analysis was used to guide a change in nozzle geometry to induce a self-excited PVC oscillation. While comparing estimates of critical regions from traditional physics-based stability analysis and data-driven complex network theory methods still remains an open question to the best of our knowledge, we apply the latter method due to the lower uncertainty in the results that it affords for the data analysed in this paper.

Following Krishnan et al. (Reference Krishnan, Manikandan, Midhun, Reeja, Unni, Sujith, Marwan and Kurths2019), we construct an undirected weighted spatial correlation network using the velocity fields obtained from the sPIV measurements as follows. The points where sPIV measurements are available are considered as nodes of the network. The strength of the connection between two nodes $i$ and $j$ is established using the Pearson's correlation coefficient $R_{ij}$, defined using the in-plane transverse velocity component, $u_{r}$ as follows:

(2.1)\begin{equation} R_{ij} = \frac { \displaystyle\sum\limits_{n=1}^{N} (u^{n}_{r,i} - \langle u_{r,i} \rangle) (u^{n}_{r,j} - \langle u_{r,j} \rangle) } { \sqrt{ \displaystyle\sum\limits_{n=1}^{N} (u^{n}_{r,i} - \langle u_{r,i} \rangle)^{2} } \sqrt{ \displaystyle\sum\limits_{n=1}^{N} (u^{n}_{r,j} - \langle u_{r,j} \rangle)^{2} } }, \end{equation}

where the superscript ‘$n$’ indexes time and subscripts ‘$i$’ and ‘$j$’ index spatial location in the PIV field of view. Angle brackets represent a temporal average computed over $N = 9300$ samples for the non-reacting cases, and $N = 10\,000$ samples for the reacting cases. The connectivity between nodes in the network, i.e. spatial locations, is then represented as an adjacency matrix, $A_{ij}$ defined as follows:

(2.2)\begin{equation} A_{ij} = \begin{cases} |R_{ij}|, & |R_{ij}| \geqslant R_t\ \mbox{and}\ i \neq j \\ 0, & \text{otherwise} \end{cases} \end{equation}

where we choose the threshold $R_{t}=0.5$. Thus, the value of $A_{ij}$ for a given pair of nodes quantifies the magnitude of the correlation between flow motions at the two nodes and is taken as the strength of the link between these nodes. We have repeated the analysis using other velocity components (not shown). The resulting adjacency matrices are qualitatively and quantitatively similar in all cases.

We then compute the weighted closeness centrality of each node in the network as follows (Opsahl, Agneessens & Skvoretz Reference Opsahl, Agneessens and Skvoretz2010):

(2.3)\begin{equation} C_{i} = \sum\limits_{j=1,j \neq i}^{N} 2^{{-}d_{ij}}. \end{equation}

Here, $d_{ij}$ is a the least costly path between two nodes defined as follows:

(2.4)\begin{equation} d_{ij} = \text{min} \left( \frac{1}{A_{ih}} + \frac{1}{A_{hk}} + \cdots + \frac{1}{A_{rs}} + \frac{1}{A_{sj}} \right), \end{equation}

where $h,k,r,s,\ldots$ are the intermediary nodes on some path connecting nodes $i$ and $j$ in the network. The minimization in (2.2) is performed over all possible paths between nodes $i$ and $j$ as determined from $A_{ij}$ using Dijkstra's algorithm (Dijkstra Reference Dijkstra1959). From (2.3), it can be shown that $C_{i}$ is large when the node $i$ is connected to many nodes in the network by low cost paths. From (2.4) the cost of the path is determined by how closely velocity oscillations between two connected nodes are correlated.

Thus, the value of $C_{i}$ defined by (2.3), measures how closely flow oscillations at node $i$ are strongly correlated with flow oscillations at all other points in the flow. Thus, a spatial region corresponding to a cluster of nodes with large $C_{i}$ values, suggests the presence of strong feedback coupling between flow oscillations at points in this spatial region and that this coupling is essentially responsible for the flow oscillations. Therefore, clusters of nodes with large values of $C_{i}$ in coherently self-excited flows, identify critical regions in the flow that are analogous to ‘wavemaker’ regions determined from traditional stability analysis. Similarly, these high $C_{i}$ clusters in forced flows essentially identify regions of high receptivity where the presence of significant amplitude of imposed forcing drives coherent unsteady flow oscillations. For the sake of simplicity, in this paper, we refer to all critical regions identified by closeness centrality as wavemakers because broadly speaking, flow oscillations in these regions drive flow oscillations in other regions of the flow. Further discussion of the complex network analysis implementation and its relevance to wavemakers is provided in the supplementary material.

2.3.3. Flame dynamics analysis

The signal measured from the OH-PLIF is used to characterize the flame at every instant. In order to observe the impact that the flame structure has on the flow field, we need to define a flame edge. Here, we assume that the flame edge can be defined at the edge of the OH-PLIF signal. In a fully premixed flame, the flame edge can be identified using a maximum gradient definition from an OH-PLIF image (Fugger et al. Reference Fugger, Roy, Caswell, Rankin and Gord2019). Even though the flame in the present experiment is nominally partially premixed, instantaneous acetone-PLIF images obtained concurrently with the OH-PLIF images show good mixing of fuel and air ahead of the regions with significant OH-PLIF signal intensity. Figure 2 shows a typical result that superimposes an instantaneous snapshot of the OH-PLIF signal intensity (red contour) with the corresponding concurrently measured acetone PLIF intensity (green contour) on the left, along with a time-averaged representation of both fields for the same air split condition ($\text {air split}=0.25$). The drop off in intensity in the acetone-PLIF signal ahead of the regions of high OH-PLIF intensity justifies our use of a fully premixed flame edge definition based on the instantaneous OH-PLIF gradient. Further analysis of the two PLIF results and more detailed discussion of image processing can be found in Karmarkar, Boxx & O'Connor (Reference Karmarkar, Boxx and O'Connor2022).

Figure 2. (a) Instantaneous and (b) time-averaged images of superimposed OH-PLIF (red contour) and acetone-PLIF (green contour) signals.

We obtain time-averaged progress variable fields by binarizing every frame based on the flame instantaneous flame edges extracted from the instantaneous OH-PLIF field snapshots. The raw OH-PLIF images are first filtered using a bilateral filter to remove noise and smooth sharp gradients. The images are then binarized using a multi-level thresholding function, multithresh in MATLAB, where the reactant regions are assigned a value of zero and the products are assigned a value of one. The thresholds are calculated using Otsu's method and a maximum of four levels are used for the binarization. Finally, the binarized images are averaged in time to obtain the time-averaged progress variable fields, $\bar {c}$, corresponding to every operating condition. We use $\bar {c}=0.5$ as a reference flame location at every condition to characterize the flame shape.

3.1. Time-averaged flow profiles

The main operational parameter varied in this study is the air split, which is the ratio of the mass flow through the central portion of the nozzle to the total mass flow. At the non-reacting condition, 11 air splits were tested, varying from 0 (no air flow through the centre) to 1 (all air flow through the centre) in steps of 0.1. At all conditions, the total mass flow rate is held constant. No fuel flow is present through the slots between the inner and outer swirlers.

Figure 3 shows the time-averaged streamwise velocity contours for all non-reacting air split conditions obtained from the sPIV measurements. The streamlines are depicted using white curves and the red dotted curves are the zero contours of time-averaged streamwise flow velocity, i.e. $\bar {u}_z=0$. As the air split increases, two features of the swirling jet change. First, as air split increases and more air flows through the central swirler, the annular jet becomes more compact. This change is because more flow is injected near the centreline and because the radial spread of the annular jet decreases as air split increases. Second, vortex breakdown results in the formation of a recirculation zone, which develops and becomes stronger and narrower as the air split increases.

Figure 3. Time-averaged streamwise velocity for all non-reacting conditions with streamlines depicted as white curves and the limits of the recirculation zone ($\bar {u}_z=0$) depicted with dotted red contours.

The dynamic swirl number, $S_{d}$, is determined from the ratio of the axial fluxes of azimuthal and axial momentum, radially integrated from the $\bar {u}_z=0$ contour at the outer boundary of the recirculation zone ($r=r_c$) to the edge of the frame of data ($r=R$) at each downstream distance as follows:

(3.1)\begin{equation} S_d=\frac{\displaystyle\int_{r_c}^{R} \bar{u}_{\theta}(z,r) \bar{u}_{z}(z,r) r\,{\rm d}r} {\displaystyle\int_{r_c}^{R} \bar{u}_{z}(z,r)^{2} r \,{\rm d}r}. \end{equation}

We choose to integrate the momentum from the edge of the recirculation zone in order to capture only the swirling part of the flow. Figures 4(a) and 4(b) show the centreline variations of the time-averaged axial velocity ($\bar {u}_z$) and $S_{d}$, respectively, at every air split condition. As the air split increases, the magnitude of the reverse velocity component increases, indicating an increase in the strength of recirculation. Note that the lowest three air splits have very little recirculation along the centreline, indicating that the central recirculation zone structure is largely determined by the flow through the central swirler. Furthermore, the swirl number, $S_d$, does not vary significantly across different air split conditions, especially close to the inlet. This invariance is likely due to the fact that the swirl number of the inner and outer swirlers is very similar.

Figure 4. (a) Centreline time-averaged axial velocity, (b) swirl number and (c) backflow ratio and swirl number at the downstream location of minimum mean centreline axial velocity for all air split conditions.

These flow profiles have been condensed into two parameters – the backflow ratio and the swirl number – in figure 4(c). The backflow ratio is defined as follows (Yu & Monkewitz Reference Yu and Monkewitz1990):

(3.2)\begin{equation} \varLambda = \frac{|\bar{u}_{z,o}|-|\bar{u}_{z,max}|}{|\bar{u}_{z,o}|+|\bar{u}_{z,max}|}, \end{equation}

where $\bar {u}_{z,o}$ is the maximum streamwise reverse flow velocity magnitude on the axis, i.e. at $z\sim 8$ mm for all cases – see figure 4(a). The quantity $\bar {u}_{z,max}$ is the maximum off-axis streamwise velocity magnitude at this downstream distance. Therefore, values of $\varLambda$ close to zero or positive imply strong centreline reverse flow and thus strong recirculation.

Figure 4(c) shows the variation of $\varLambda$ with air split. Note that $\varLambda$ is only calculated for air splits of 0.3 and above, as cases with lower air splits do not show a minimum in centreline axial velocity in the field of view – see figure 4(a). The corresponding variation of swirl number, $S_{o}$, at the same axial location is also plotted for reference. Together, these results show that increasing air split results in an increase in recirculation strength and swirl up to an air split value of 0.5. Beyond this value, both swirl intensity and recirculation strength stay constant. Previous studies of swirl flow stability using similar parameters (Oberleithner et al. Reference Oberleithner, Sieber, Nayeri, Paschereit, Petz, Hege, Noack and Wygnanski2011; Manoharan et al. Reference Manoharan, Hansford, O'Connor and Hemchandra2015) show that increasing swirl and recirculation strength causes flow profiles in a swirl flow to become locally absolutely unstable. As such, these two parameters play a key role on the global stability of the PVC to leading order (Monkewitz, Huerre & Chomaz Reference Monkewitz, Huerre and Chomaz1993). The nearly constant values of $\varLambda$ and $S_{o}$ for air splits larger than 0.5 in figure 4(c) indicate that the flow through the central nozzle essentially governs the coherent unsteady dynamics of the swirling flow field. The impact of the reduction in outer nozzle mass flow rate is not significant. We will show further evidence of the importance of the centre-nozzle flow from SPOD analysis of the non-reacting velocity fields.

3.2. Spectral characterization of non-reacting flow

The time-averaged velocity fields for the non-reacting flow conditions shown in figure 3 suggest that the flow through the centre nozzle governs the structure of the swirling flow field and its transition to vortex breakdown. This transition results in the emergence of a coherent hydrodynamic oscillation in the flow, as has been observed in several prior studies (Escudier & Keller Reference Escudier and Keller1985; Billant, Chomaz & Huerre Reference Billant, Chomaz and Huerre1998; Oberleithner et al. Reference Oberleithner, Sieber, Nayeri, Paschereit, Petz, Hege, Noack and Wygnanski2011; Manoharan et al. Reference Manoharan, Frederick, Clees, O'Connor and Hemchandra2020). Figure 5 shows the variation of modal energy with oscillation frequency, referred to hereafter as modal energy spectra, for all modes obtained from SPOD for all air split conditions. The red curves show the energy spectra for the most energetic mode and the grey curves are for modes in decreasing order of energy at each frequency. As the air split increases, a dominant coherent mode can be seen, starting at 254 Hz at an air split of 0.2. Both the frequency and amplitude of this mode increases with increasing air split; harmonics of the mode are also present in all cases.

Figure 5. Modal energy spectra from SPOD for non-reacting conditions varying air splits.

Figure 6 shows the variation of the peak frequency ($\,f_{peak}$) with air split, along with the same data plotted as a non-dimensional Strouhal number $St=\,f_{peak} \bar {U}_b/D_{i}$, using the bulk velocity through the central nozzle as the reference velocity ($U_{b}$) and the diameter of the central nozzle ($D_{i}$) as the reference length scale. The frequency of the coherent mode increases linearly with increasing air split, while the Strouhal number remains constant. This Strouhal number behaviour is characteristic of hydrodynamic instability rather than an acoustic mode. The variation of the PVC frequency with flow split, and hence the velocity of the jet through the centre nozzle, is similar to results from several studies in a number of different experimental configurations (Syred & Beer Reference Syred and Beer1974; Martinelli, Olivani & Coghe Reference Martinelli, Olivani and Coghe2007; Oberleithner et al. Reference Oberleithner, Sieber, Nayeri, Paschereit, Petz, Hege, Noack and Wygnanski2011; Moeck et al. Reference Moeck, Bourgouin, Durox, Schuller and Candel2012; Manoharan et al. Reference Manoharan, Frederick, Clees, O'Connor and Hemchandra2020). In all these cases, the frequency scales with the nozzle flow velocity, resulting in a constant Strouhal number throughout. The length scale used in the Strouhal number is therefore the centre-nozzle diameter, rather than the effective diameter, since the flow through the centre nozzle determines the behaviour of the system. This trend was formally shown in the work of Manoharan et al. (Reference Manoharan, Frederick, Clees, O'Connor and Hemchandra2020), wherein the PVC frequency depends only on the time-averaged flow state at a critical swirl number $S_{c}$, where vortex breakdown is first initiated for a given fixed mass flow rate and the swirl number $S$. The PVC emerges as a result of a supercritical Hopf bifurcation, creating a stable limit cycle oscillation in the flow. For high $Re$ swirl flows, $S_{c}$ for a fixed inner nozzle mass flow rate is a weak function of $Re$ (Escudier Reference Escudier1988). Therefore, the nearly constant $St$ in figure 6 suggests that $S>S_{c}$ for the inner nozzle flow for all values of air split greater than 0.3. Further, the nearly constant values of $\varLambda$ and $S_{o}$ in figure 4(c) suggest that $S_{c}$ is the same for all values of air split beyond 0.3, i.e. where the onset of PVC is observed – see figure 5. Thus, analysis in Manoharan et al. (Reference Manoharan, Frederick, Clees, O'Connor and Hemchandra2020) shows that, for a fixed value $S_{c}$ and $S$, the PVC oscillation amplitude and frequency must scale linearly with inner nozzle bulk flow velocity yielding a constant $St$, consistent with the result in figure 6.

Figure 6. Peak frequency and Strouhal number of Mode 1 as a function of air split for non-reacting conditions.

To understand the structure of the spatial oscillation field of this mode, we reconstruct velocity SPOD modes corresponding to each frequency peak at every air split condition. Figure 7 shows the transverse velocity component from the SPOD mode at the first peak frequency for three air splits (0.3, 0.5 and 0.8), as well as the closeness centrality parameter, $C_{i}$ (see (2.3)), for four air splits (0.1, 0.3, 0.5 and 0.8). In the SPOD modes, the black dotted line represents $\bar {u}_{z}=0$ and demarcates the limits of the recirculation zone. The cases shown all have a strongly coherent transverse velocity oscillation at the base of the vortex breakdown bubble, showing that the flow precesses about the inner nozzle centreline. The mode also shows helical rollup along the inner shear layer. This spatial SPOD mode structure is characteristic of a PVC oscillation seen in past studies of this instability (Tammisola & Juniper Reference Tammisola and Juniper2016; Manoharan et al. Reference Manoharan, Frederick, Clees, O'Connor and Hemchandra2020; Datta et al. Reference Datta, Gupta, Chterev, Boxx and Hemchandra2021; Gupta et al. Reference Gupta, Shanbhogue, Shimura, Ghoniem and Hemchandra2021). The coherent rollup along the shear layer is due to the helical forcing imposed on the shear layers, as has been argued in recent studies (Oberleithner et al. Reference Oberleithner, Terhaar, Rukes and Paschereit2013; Manoharan et al. Reference Manoharan, Frederick, Clees, O'Connor and Hemchandra2020). Also, the fact that the coherent structures in the transverse velocity field are are symmetric about the central axis indicates that the shear layer rollup is radially asymmetric and helical. Thus, the mode shape and the helically paired vortical structures seen along the edge of the recirculation zone confirm that the coherent mode associated with the frequency peak in the modal spectrum in figure 5 is the manifestation of a global hydrodynamic limit cycle oscillation in the form of a PVC.

Figure 7. Spectral POD mode shapes corresponding to the PVC frequency (a–d) with the critical region identified by the weighted closeness centrality measure (e–h) from the complex network analysis computed for the transverse velocity component at various air splits.

Figure 7(e–h) shows regions of large $C_{i}$ in the flow at every air split value. The large region located at $z\sim 20$ mm at the lowest air split of 0.1 corresponds to turbulence with some residual level of correlation. We can conclude this because the SPOD results for this case show no evidence of narrowband tonal oscillations – see figure 5). We also verify this region of large $C_{i}$ is driven by turbulence by re-computing the closeness centrality parameter after removing the nodes from the upstream critical regions at the base of the bubble and along the shear layers. When the nodes from these critical regions are eliminated, the amplitude of the closeness centrality parameter in the central downstream region remains unchanged. This lack of change indicates that the nodes in this downstream region are not correlated with the nodes at the base of the bubble or along the shear layer, but rather well correlated with other nodes in the same downstream region, as would be expected for turbulence but not a coherent oscillation. A more detailed explanation of this method is provided in the supplementary material. The emergence of a PVC around air split value of 0.3 results in the emergence of a new compact large $C_{i}$ value region at the upstream end of the vortex breakdown region. As air split increases, this second large $C_{i}$ region becomes dominant with increasing air split. The prior studies of Tammisola & Juniper (Reference Tammisola and Juniper2016) and Mukherjee et al. (Reference Mukherjee, Muthichur, More, Gupta and Hemchandra2021) identify wavemaker regions in swirl nozzle flows from linear global hydrodynamic stability analysis. We note here that the cluster of large $C_{i}$ at air split values of 0.3–0.5 in figure 7 are located at a similar position relative to the central recirculation zone. Also, slightly smaller values of $C_{i}$ are observed along the shear layers that align closely with the shear layer oscillations in the corresponding SPOD result. This similarity leads us to conclude that the region upstream of high $C_{i}$ at the upstream end of the breakdown bubble corresponds to the wavemaker of the PVC oscillation in the flow, for reasons discussed earlier in § 2.3.2. Note also that figure 7 shows that the position and shape of the wavemaker remains largely unchanged with increasing air split. This result lends further support to the fact that the dynamics of the PVC in the present experiment is governed mainly by the flow through the centre nozzle.

The results for $C_{i}$ from the non-reacting result in figure 7 have several important general implications that we can learn from and apply to the reacting flow data next. Similarity between the spatial distributions of $C_{i}$ and SPOD mode oscillation amplitudes is an important indicator of a critical motion in the flow. Regions where $C_{i}$ and the SPOD spatial mode amplitude at a prominent frequency in the modal spectrum have similar spatial distributions are evidence of high-energy, highly correlated motions like vortex precession in the present case. This result is especially true because the method used to construct $C_{i}$ is applied on raw time-series data and is, prima facie, not spectrally selective. For the same reason, regions with large values of $C_{i}$ that do not align with coherent motions identified by SPOD mode shapes may have strong correlation but are not narrowband oscillations, as shown by the 0.1 air split result in figure 7. Put together, the information in the SPOD modes (narrow-band, high-energy fluctuations) and the closeness centrality parameter (high spatial and temporal correlation in a given region of space) provide important complementary insight into the motions in a flow and their origins.

4.1. Time-averaged flow and flame profiles

The time-averaged streamwise velocity contours with streamlines, obtained from the sPIV measurement, for all reacting air split conditions are shown in figure 8. In all cases, a region of negative axial velocity, the central recirculation zone (CRZ), can be seen close to the axis of symmetry, similar to the structure seen in the non-reacting cases. As air split increases, this recirculation becomes stronger and more well defined. As in the non-reacting data, the red dotted curve is the $\bar {u}_{z}=0$ contour showing the extent of the central recirculation zone. As the air split increases and the flow velocity through the central nozzle increases, the recirculation zone becomes more compact – see figure 8. The CRZ is quite wide for air splits of 0.2 and 0.25, whereas the CRZ width and overall shape for air splits of 0.35–0.5 are similar with the change in shape occurring at an air split of 0.3.

Figure 8. Time-averaged streamwise velocity component with streamlines (white solid lines) with the mean flame location ($\bar {c}=0.5$) depicted by the solid black lines and the limits of recirculation zone ($\bar {u}_z=0$) depicted as dotted red lines.

The black contour in figure 8 is the time-averaged progress variable $\bar {c}=0.5$, which is indicative of the most likely position of the flame. Thus, these results show that at all air splits, the flame is anchored at the base of the CRZ, near the centreline of the flow. For air splits of 0.2–0.3, the distinctive ‘M’-shape of the $\bar {c}=0.5$ contour suggests that the flame not only stabilizes in the shear layer around the CRZ, but also in the outer shear layer formed by the outer nozzle jet and the ambient fluid. In these cases, air flow through the outer nozzle is more substantial and the shear generated in that region possibly helps with flame stabilization. For air splits larger than 0.3, the $\bar {c}=0.5$ contour shape suggests that the flame is no longer attached at the outer lip of burner due to the decrease in the outer nozzle mass flow rate. Note that in all cases, however, the flame structure near the CRZ is somewhat insensitive to air split. The time-averaged flame near the centreline matches the base of the zero axial velocity contour in that region for all air splits above 0.2.

Figure 9 shows the variation of the time-averaged axial velocity along the centreline and swirl number $S_d$, as defined in (3.1), with downstream distance at every air split condition. Similar to the non-reacting case, figure 9(a) shows that as air split increases, the strength of recirculation increases because more air flows through the centre nozzle. Note that the streamwise position of the recirculation zone stays largely constant with varying air split. Additionally, figure 9(b) shows that the swirl number remains nominally unchanged by air split. These results mirror the non-reacting case for the air split range (0.2–0.5), although strong recirculation is observed at a lower air split in the reacting case (0.2) than in the non-reacting case (0.3). Enhancement of reverse flow velocity in the presence of a flame has been shown in both experimental (Syred & Beer Reference Syred and Beer1974) and theoretical (Rusak, Kapila & Choi Reference Rusak, Kapila and Choi2002) studies.

Figure 9. (a) Centreline time-averaged axial velocity, (b) swirl number and (c) backflow ratio and swirl number at the downstream location of minimum mean centreline axial velocity for all air split conditions.

Figure 9(c) shows the backflow ratio ($\varLambda$) (3.2) and swirl number ($S_{o}$) at the maximum axial backflow velocity location on the centreline, as was shown for the non-reacting conditions in figure 4. The trends for $\varLambda$ and $S_{o}$ between the non-reacting and reacting cases are very similar, recognizing that the air split range in the reacting case is half that of the non-reacting case. In both cases, $\varLambda$ increases from the onset of vortex breakdown until an air split of approximately 0.45 or 0.5, and then plateaus. Note that it is difficult to extrapolate a trend past an air split of 0.5 for the reacting case, but the trend ahead of that point is very similar to the non-reacting cases. The value of $\varLambda$ is less negative as a result of both the increasing axial velocity in the jet and the centreline recirculation velocity.

The values of $S_{o}$ are quite different for the non-reacting and reacting cases. In the non-reacting case, $S_{o}$ initially decreases with increasing air split because the location of minimum centreline axial velocity moves downstream. Figure 4(b) shows that the swirl numbers are decrease as a function of downstream distance for the first 15 mm downstream of the nozzle exit as a result of the shear-induced decay of the swirl velocity. In the reacting cases, however, flow expansion due to the presence of the flame makes the structure of the vortex breakdown bubble and the flow around it more uniform at air splits above 0.2.

4.2. Spectral characterization of reacting flow

The time-averaged flow profiles indicate that increasing air split leads to strengthening of recirculation in the vortex breakdown region. This increase in recirculation strength can lead to changes in the hydrodynamic instability characteristics of the flow. This may be attributed to several fundamental mechanisms that cause the onset of local absolute instability in the flow as shown by prior fundamental local hydrodynamic stability studies (Gallaire & Chomaz Reference Gallaire and Chomaz2003; Manoharan et al. Reference Manoharan, Hansford, O'Connor and Hemchandra2015; Oberleithner et al. Reference Oberleithner, Stöhr, Im, Arndt and Steinberg2015). A key difference between the reacting and non-reacting cases in the present study, however, is the presence of thermoacoustic instability in the reacting case.

Combustor acoustic pressure power spectral densities at each air split condition are shown in figure 10. Pressure is measured using a probe placed through the end cap of the combustion chamber, which is located 120 mm from the dump plane. The power spectral density shows multiple distinctive peaks at each condition, indicating the presence of multiple oscillations. The presence of multiple frequencies is characteristic of reacting systems exhibiting hydrodynamic and thermoacoustic instability. In the lower air split cases, a strong peak can be seen at around 600 Hz (marked by a circle), which is identified to be due to acoustic pressure oscillations resulting from thermoacoustic instability.

Figure 10. Acoustic pressure spectra at each air split condition.

Acoustic modelling of the combustor rig using COMSOL (not shown), shows two acoustic modes with strong coupling between the combustor and the inner passage of the injector at 570 and 960 Hz, similar to the frequencies measured in the data. The COMSOL model assumed only two temperatures in the system for reactants and products, where the air in the inlet plenum is at the reactant temperature and the air in the combustor and exhaust section is at the equilibrium flame temperature of products corresponding to the overall global equivalence ratio. Also, we have assumed the boundary conditions upstream of the injector plenum and downstream of the exhaust as hard wall. These simplifications were needed as no acoustic boundary conditions were characterized in the experiment. Still, the modes captured in the model are in reasonable agreement with the peak frequencies marked with circles and triangles in figure 10, showing that the latter are thermoacoustic in origin.

As air split increases, the peak at around 600 Hz persists, but decreases in amplitude. The second peak at around 850 Hz is identified as the second thermoacoustic mode, marked with a triangle. Furthermore, harmonics of the dominant peaks can also be seen. Note also that the two peaks do not change appreciably with air split, as would be expected of an acoustic, rather than hydrodynamic, oscillations. We refer to the coherent flow oscillation modes induced by thermoacoustic oscillations at $\sim$600 Hz as TA mode 1 and those at $\sim$900 Hz as TA mode 2 in the remainder of this paper. The flow oscillation mode associated with the PVC will be referred to simply as the PVC mode.

The spatial distribution of coherent flow oscillations at these acoustic frequencies are determined from corresponding velocity field measurements using SPOD. The modal spectra of the velocity data at every reacting air split condition are shown in figure 11. As before, the red curve in each plot depicts the highest-energy mode and the grey curves depict subsequent modes in order of decreasing energy. In all cases, the coherent content in the flow field measurements at the peak frequencies is captured mainly by the first SPOD mode. Therefore, we restrict our analysis to the mode shapes of the first SPOD mode alone.

Figure 11. Modal spectra of velocity fields obtained from spectral POD at all air splits for the reacting conditions.

Similar to the acoustic pressure spectra, multiple peaks can be seen at each condition, which indicates that there are multiple coherent oscillations in the flow. In order to understand the nature of the oscillation corresponding to every peak, the SPOD mode shapes of the first mode at the respective peak frequencies are reconstructed at every air split; these shapes will be discussed at length in this section. Using the mode shapes, we determined the source of each peak, where the PVC is indicated with a black star, TA mode 1 with a pink circle, and TA mode 2 with a blue triangle. Several modes at frequencies corresponding to the difference between the TA modes and the PVC are also seen and have been marked with green squares. We will refer to these modes hereafter as interaction modes.

The variation of the frequencies and amplitudes of the three basic flow oscillation modes with air split from SPOD are depicted in figure 12. From the plots in figure 12 and the modal spectra in figure 11, two important characteristics of the flow oscillations in the system can be observed. First, figure 12(a) shows that, while the frequency of the PVC mode increases linearly with increasing air split, the thermoacoustic frequencies remain largely the same at all air split conditions; this is further confirmation that TA modes 1 and 2 are indeed due to thermoacoustic oscillations. Second, figure 12(b) shows that increasing amplitude of the PVC corresponds to decreasing amplitude of TA mode 1; TA mode 2 is relatively weaker at all conditions. A key observation from figure 12 is that as the air split is increased, the PVC mode frequency sweeps through those of the two TA modes, the amplitude of the flow oscillation associated with the thermoacoustic mode closest to the PVC mode frequency is suppressed. For TA mode 1, this occurs at air splits of 0.35 and 0.4, and for TA mode 2, this occurs at an air split of 0.45 and 0.5. This result indicates that when present, the PVC not only dominates the velocity response of the flow but can also suppress the response of the flow to forcing imposed by thermoacoustic oscillations when their frequencies overlap.

Figure 12. Frequency (a) and amplitude (b) characterization of the three oscillation modes seen at each air split condition from SPOD.

The mechanism of the thermoacoustic mode suppression is investigated in two ways: through complex network analysis of the time-series velocity data coupled with the SPOD spatial mode shapes and through a weakly nonlinear theoretical analysis. First, we analyse the SPOD modes of the velocity field along with the results of the complex network analysis discussed in § 2. Figure 13 shows the SPOD spatial mode shapes based on the radial velocity component at three frequencies for each of the air split cases: the PVC frequency, the dominant TA frequency and the interaction, or difference, frequency of those two modes. In addition to these mode shapes, we show $C_{i}$ maps from complex network analysis. The first column of results shows the mode shapes for the PVC, whose frequency increases with air split and increasing flow velocity through the central nozzle. In these modes, a region of high-amplitude oscillation is located at the base of the recirculation zone, indicative of the precession of the vortex core around the nozzle centreline. Symmetric disturbances in the shear layers emanate from this region, where, because this mode shape depicts the radial velocity component, $u_r$, the symmetric disturbances indicate asymmetric oscillations. The structure of these oscillations is very similar to the structure of the PVC mode in the non-reacting case, shown in figure 7. The main difference between the two is that in the reacting case, the shear layers spread at a wider angle than in the non-reacting case, a result of the gas expansion from the flame. Similar differences can be seen in the structure of the recirculation zone from the time-averaged flow profiles in figures 3 and 8.

Figure 13. Spectral POD mode shapes at PVC, thermoacoustic and interaction frequencies, and weighted closeness centrality measures computed for the radial velocity component at all air split conditions.

The next two columns show the mode shapes at thermoacoustic mode frequencies, which stay relatively constant as a function of air split. These mode shapes have strong fluctuations in the shear layers and no fluctuation in the CRZ of the flow. The fluctuations in the shear layer look asymmetric, but because the transverse velocity component is pictured, this structure is indicative of symmetric oscillations in the shear layer. This symmetric motion is the expected hydrodynamic response of the flow when forced by a longitudinal acoustic field that is nominally axisymmetric near the burner lip, as shown in previous swirling flow studies by O'Connor & Lieuwen (Reference O'Connor and Lieuwen2012). We will refer to this response as an $m=0$ response, where $m$ is the azimuthal wavenumber of the flow disturbance. Analysis of the wavelength of the velocity disturbance as well as the bulk velocity confirm that this is the $m=0$ mode at both thermoacoustic frequencies, rather than a higher-order axisymmetric mode. The wavelength of the convective disturbances in the shear layers is shorter in TA mode 2 than TA mode 1 because TA mode 2 has a higher frequency. The impact of the symmetric oscillations on flame behaviour is discussed in detail in § 5.

The fourth column shows the mode shape for the interaction frequency between the PVC and the dominant thermoacoustic mode; no interaction modes are shown where the PVC is not present. Interaction frequencies, particularly sum and difference frequencies, commonly occur in the presence of nonlinear interaction between two periodic oscillations. In this case, the interaction frequency is the difference frequency between the PVC mode and the dominant TA mode. Spatially, the interaction mode shows up in the shear layers when the thermoacoustic mode is present (air splits of 0.3 and 0.35) and in the shear layers and the CRZ, particularly near the base, when the PVC mode is dominant (air splits of 0.4–0.5).

The final column shows the results of the complex network analysis for the reacting flow. We use these results in concert with the SPOD modes to propose a mechanism by which the thermoacoustic oscillations are suppressed by the PVC. In the non-reacting flow, the complex network analysis clearly identified the region of the PVC wavemaker, i.e. the region of the flow driving the self-excited oscillation. In a field where multiple oscillations are present, however, this time-series method provides valuable insight into the interaction between modes. Since the Pearson correlation coefficient is calculated in time, it naturally identifies close nonlinear coupling of oscillations with multiple frequencies, as is seen in the coupling between the PVC mode and the thermoacoustic mode. In the air split 0.2 case, the only oscillation present in the flow is TA mode 1. The $C_{i}$ field closely mimics the 664 Hz SPOD mode, indicating strong velocity oscillations in the shear layers. Both the SPOD mode and the centrality parameter identify the well-known velocity-coupled mechanism by which the longitudinal acoustic mode drives vortical oscillations in the shear layer, which in turn drive flame response; evidence of direct flame response is discussed next in § 5. The oscillations appear at the thermoacoustic oscillation frequency in the SPOD mode, as SPOD captures the location of high-energy fluctuations at a given frequency. The oscillations in the $C_{i}$ field now identify flow regions that drive the hydrodynamic response of the shear layers. Coherent, axisymmetric vortex shedding in the shear layer results in regions of high coherence where vortices are shed. These two regions of high SPOD amplitude and $C_{i}$ overlap spatially just as was the case in the non-reacting data with the PVC.

Figure 12(b) shows that, for the air split 0.3 case, both PVC and TA mode 1 modal energy amplitudes are comparable. Therefore, the spatial distribution of $C_i$, as shown in figure 13, is the result of both PVC and TA mode 1 oscillations. The somewhat different qualitative structure of $C_{i}$ for this case when compared with those for air splits greater than 0.3 is because of the presence of these two modes with similar amplitudes. The large region of high $C_{i}$ on the centreline for the air split 0.3 case is due to the PVC. We show this by comparing results for the spatial distribution of $C_{i}$ for reduced datasets that selectively exclude time-series contributions from the PVC and TA mode 1. The details of how this analysis is performed is discussed in § 1.2 of the supplementary material. This analysis shows that the spatial distribution of $C_{i}$ for the raw dataset and the reduced dataset that excludes the TA mode 1 contributions show the closest qualitative and quantitative match. Physically, this means that the PVC wavemaker region identified by complex network analysis is located further downstream of the centre nozzle for the air split 0.3 case than for the higher air splits. This location is the result of the lower mass flow rate through the centre nozzle as compared with the outer nozzle and the impact of that this fact has on flame attachment and flow field structure in the air split 0.3 case.

At an air split of 0.35, the PVC oscillation is strong, the TA mode 1 amplitude is small and the TA mode 2 amplitude is at its peak, according to the SPOD spectra in figure 11. Additionally, a strong interaction frequency appears in the spectrum at 293 Hz. The fourth row of figure 13 shows the SPOD modes at 508 Hz (PVC), 859 Hz (TA mode 2) and 352 Hz (PVC/TA2 difference frequency). The PVC mode shape shows asymmetric oscillations in the shear layer and at the base of the CRZ, whereas the thermoacoustic mode shape symmetric oscillations in the shear layers. The interaction mode shows evidence of both, with oscillations in the shear layers where the PVC and TA modes overlap, as well as motions at the base of the recirculation zone. The $C_{i}$ field at this condition shows oscillations, both at the base of the recirculation zone as well as in the shear layers. Oscillations in these regions are indicative of both thermoacoustic and PVC oscillations, as is the interaction SPOD mode shape at 352 Hz. The structure of the closeness centrality parameter for the higher air splits where there is significant interaction between the PVC and the TA modes are similar; closeness centrality is higher in shear layer regions where the PVC and the thermoacoustic mode could interact.

As the air split increases, the PVC oscillation increases in strength and the thermoacoustic oscillations decrease in strength, resulting in the emergence of a high $C_{i}$ region at the upstream end of the CRZ. This result is similar to what was observed with the emergence of PVC oscillations in the non-reacting flow. For example, at an air split 0.5, the TA mode 1 amplitude is quite low and the coherence of the SPOD mode at 586 Hz is relatively low; the oscillations diffuse quickly with downstream distance as compared with the low air splits where the axisymmetric oscillations from the TA mode 1 acoustic oscillations retain high coherence throughout the flow field. TA mode 2 at 918 Hz has a relatively weak amplitude, although the spatial coherence of the axisymmetric vortex shedding is visible. The interaction mode at 176 Hz between the PVC and TA mode 2 shows oscillations along both the centreline at the base of the recirculation zone and in the shear layers. The $C_{i}$ parameter again has a similar structure to the interaction mode, displaying a region of strong coherence at the base of the recirculation zone, as well as oscillations in the shear layer, indicative of interaction between the PVC and TA mode 2. These structures in both the interaction mode and $C_{i}$ are indicative of nonlinear interaction between the PVC and the thermoacoustic mode.

The second way in which we can understand the suppression of the thermoacoustic mode by the PVC is from a weakly nonlinear analysis of the flow governing equations. This analysis is presented in full detail in the supplementary material and is a generalization of the analysis of Manoharan et al. (Reference Manoharan, Frederick, Clees, O'Connor and Hemchandra2020) to include weak axisymmetric harmonic forcing. Here, the term ‘weak’ implies that the amplitude of the response of the flow, e.g. velocity amplitudes in shear layers, to forcing by thermoacoustic oscillations is much larger than the amplitude of the forcing itself. Note that in the present study SPOD modes in figure 13 for TA modes 1 and 2 show large-amplitude velocity oscillations along shear layers when compared with the velocity amplitudes elsewhere. This suggests that the flow in the present study is weakly forced by the thermoacoustic oscillation.

In this weakly forced limit, (2.19) in the supplementary material shows the solution to leading order for the oscillating flow field. This solution is comprised of several components as follows: $m=1$ PVC mode, the $m=0$ hydrodynamic response to forcing and $m=1$ interaction modes. This solution structure qualitatively reflects the nature of oscillations observed in the present SPOD results from reacting flow experiments. Further, (2.9) and (2.10) in the supplementary material show the governing equations for the amplitudes of the hydrodynamic response ($A_{o}$) to forcing by thermoacoustics and the PVC ($A_{1}$), respectively. From these equations, it is clear that $A_{1}$ is formally independent of the forcing, consistent with the present experiments because the PVC once established is always present. The amplitude of the hydrodynamic response, on the other hand, depends on the amplitude of the PVC and a coupling coefficient $\beta _{A_{o}A_{1}}$, determined by the interaction between the PVC and the axisymmetric hydrodynamic mode – see (2.11)–(2.15) in the supplementary material. The discussion corresponding to (2.15)–(2.18) in the supplementary material shows that the value of $\beta _{A_{o}A_{1}}$ has an inverse dependence on the difference between the PVC and thermoacoustic mode frequencies. Thus, (2.13) in the supplementary material can explain suppression of the TA modes by the PVC when their frequencies are close as being due to the suppression of the shear layer response and, thereby, heat release rate oscillations driving the thermoacoustic oscillations.

Put together, the complex network analysis and the weakly nonlinear analysis both show the direct interaction of the PVC $m=1$ mode with the acoustically excited $m=0$ mode, suppressing the action of the $m=0$ oscillations. These results also show that the thermoacoustic modes are not eliminated as a result of the flame moving away from a region of the acoustic mode where the Rayleigh criterion could be met. Because the $m=0$ oscillations are evident in both the experimental data and their suppression is evident through the impact of the $\beta _{A_{o}A_{1}}$ term in the theory, there is evidence to conclude that the PVC is indeed suppressing the acoustically induced symmetric motions.