2.1 Flow field equations

A well-validated LES code, developed at the Center of Turbulence Research (CTR) at Stanford University, was used in the current study. This code has been successfully applied for the solution of large scale turbulent combustion problems (see Ham et al. (Reference Ham, Apte, Iaccarino, Wu, Herrmann, Constantinescu, Mahesh and Moin2003) and You et al. (Reference You, Ham and Moin2008)).

The governing equations solved here are the unsteady, low Mach number filtered Navier–Stokes equations. The flow field variables are decomposed into resolved and subgrid-scale terms, by applying a grid-filtering operation. Small-scale fluctuations are removed by the filter, while their statistical effects remain in the residual terms of the filtered equations (Ferziger Reference Ferziger and Essers1983). For the subgrid-scale model, eddy viscosity and eddy diffusivities are evaluated using the dynamic Smagorinsky model, developed by Germano et al. (Reference Germano, Piomelli, Moin and Cabot1991). The filtered LES equations for continuity, momentum and scalar transport are,

where

where

state relation:

Here, $\unicode[STIX]{x1D70C}$ is the density, $\boldsymbol{u}_{j}$ is the velocity vector, $p$ is the pressure, $\unicode[STIX]{x1D707}$ is the viscosity, $\unicode[STIX]{x1D6FF}_{ij}$ is the Kronecker symbol, $\unicode[STIX]{x1D719}_{i}$ represents a scalar quantity, $\unicode[STIX]{x1D6FC}_{i}$ is the diffusivity of scalar $i$ and $\unicode[STIX]{x1D714}_{i}$ is the chemical source term of scalar $i$ . The filtering operation is denoted by an overbar, and Favre (density weighted) filtering by a tilde. The unclosed transport terms in the momentum and scalar equations are: the residual stress $t_{ij}$ , residual scalar flux $q_{ik}$ , chemical source term $\tilde{\unicode[STIX]{x1D714}_{i}}$ and the turbulent diffusivity $\unicode[STIX]{x1D6FC}_{t}$ . These terms require a closure model. The state relation is reduced to an expression for density in terms of the transported scalars $\unicode[STIX]{x1D719}_{i}$ .

2.2 Chemical model

The FPV approach developed by Pierce & Moin (Reference Pierce and Moin2004) is utilized in the present study. In the FPV approach, a tabulated chemical table is computed in advance. Parameterization of the chemical properties is performed by using the mixture fraction, mixture fraction variance and a reaction progress variable. The conservation equations for the mixture fraction and progress variable scalars are

All the thermo-chemical properties, denoted by the vector $\unicode[STIX]{x1D733}$ , are parametrized in terms of the mixture fraction and the reaction progress variable:

where $\unicode[STIX]{x1D74D}=(T,y_{i},\dot{\unicode[STIX]{x1D714}}_{C},\unicode[STIX]{x1D70C},\unicode[STIX]{x1D708},\unicode[STIX]{x1D6FC})^{T}$ . $y_{i}$ is the mass fraction of species $i$ . The filtered thermo-chemical variables are obtained by integrating equation (2.8) over a joint probability density function (PDF). That way, subgrid fluctuations over the filtered chemical variables are incorporated.

The joint PDF is modelled by

Here, $\widetilde{P}(Z)$ is given by the presumed beta PDF, $\unicode[STIX]{x1D6FD}(Z;\widetilde{Z},\widetilde{Z^{\prime \prime 2}})$ . $P(C|Z)$ is the conditional PDF which is represented by a delta function under the assumption that the scale over which a single flamelet is valid is greater than or equal to the size of the local computational cell, defined as follows:

where $C\mid Z$ is given by the flamelet solution. The resulting state-space parametrization is a function of the mixture fraction, the mixture fraction variance and the reaction progress variable $\widetilde{\unicode[STIX]{x1D74D}}=\widetilde{\unicode[STIX]{x1D74D}}(\widetilde{Z},\widetilde{Z^{\prime \prime 2}},\widetilde{C})$ . The residual scalar variance of the mixture fraction can be written as $\widetilde{Z^{\prime \prime 2}}=\widetilde{Z^{2}}-\tilde{Z}^{2}$ . An additional transport equation is solved here for $\widetilde{Z^{2}}$ as this equation does not include a source term (Shunn & Moin Reference Shunn and Moin2009).

The flamelet libraries used in this study were computed for methane using GRI-Mech 2.11, a detailed kinetic mechanism containing 277 elementary chemical reactions and 49 species (Bowman et al. Reference Bowman, Hanson, Davidson, Gardiner, Lissianski, Smith, Golden, Frenklach and Goldenberg1995). The application of PDF methods for the approximation of reactive scalar subgrid fluctuations is comprehensively discussed by Ihme et al. (Reference Ihme, Chong and Pitsch2005). Using a DNS database for validation, Ihme et al. (Reference Ihme, Chong and Pitsch2005) have shown that the steady flamelet assumption in the context of the FPV approach, which is used here, leads to good predictions of local extinction and reignition.

It should be mentioned that, in principle, the non-premixed flamelet library used in the present flamelet-progress variable approach is not expected to predict the propagation of premixed branches of the lifted flames. However, what one can see in the present study configuration are mostly non-premixed large-scale flames, which are partially detached from the inlet and are stabilized by the swirling motion of the recirculation zone. Nevertheless, in the present cases, air is preheated to 750 K, suggesting that the entrainment of hot combustion products may not be sufficient to fully stabilize the flame, which then indicates the importance of partially premixed flame propagation in those regions.

Recently, See & Ihme (Reference See and Ihme2015) have shown that both the FPV framework, assuming steady flamelets, and a premixed-based filtered chemistry model provide good prediction of the highly unsteady flame dynamics of the German aerospace centre (DLR) dual swirl gas turbine model combustor (GMTC) which operates in the partially premixed combustion regime. In their configuration, the flame is lifted from the inlet. See & Ihme (Reference See and Ihme2015) also report a flame lifted from the inlet (in simulation and experiment). These results support the use of the current FPV model for the combustion regime of our configuration, which is mainly non-premixed and may involve local regions of partially premixed combustion and lifted flames. Further details on the application of the FPV approach can be found in the paper of Pierce & Moin (Reference Pierce and Moin2004).

2.3 Numerical method

The solution scheme formulates the discretized problem over a staggered mesh, where velocity components are staggered both in space and time. The challenge of capturing the turbulent structures responsible for the rapid mixing between fuel and oxidizer, while maintaining numerical stability, is met by the use of a non-dissipative, energy conserving scheme (Ham et al. Reference Ham, Apte, Iaccarino, Wu, Herrmann, Constantinescu, Mahesh and Moin2003). The time integration method is of second-order accuracy and fully implicit.

Neumann conditions for pressure and scalars were set as wall boundary conditions. No-slip Dirichlet conditions were applied for the velocity. Outflow boundaries were implemented using convective conditions.

A possible approach to accurately generate the combustor’s inflow conditions that fit the experiment (Owen et al. Reference Owen, Spadaccini and Bowman1976) is to completely resolve the turbulent flow in the configuration located upstream. This includes the wind tunnel ducts and swirl vanes. However, this is a very demanding task. A more practical approach, and the one used in the present study, was suggested by Pierce & Moin (Reference Pierce and Moin2004) as follows: inflow conditions are generated using a separated code, which simulates a fully developed, turbulent coaxial pipe flow. For the pipe simulation, spatially periodic boundary conditions are applied at the inlet and outlet of the pipe. Recycling is then used to establish fully developed, turbulent inflow conditions. Velocity cross-section data are sampled every few time steps and saved into an inflow database. The main combustion simulation uses this inflow database to generate the time-dependent fuel and air inlet boundary conditions. This approach is also used to simulate swirling inflow conditions.

4.1 Non-swirling conditions

Solution-mesh sensitivity tests for time-averaged velocity and temperature were conducted for the non-swirling case. It should be noted that a mesh-independence test is considered a challenging task in LES. Still, significant convergence has been achieved in most of the solution region and particularly in the flame region. Further details on the sensitivity tests conducted here can be found in appendix A.

Radial profiles of time-averaged velocity, temperature, CO mass fraction and product mass fractions are compared with the experimental database (figure 2). Here, measurements that had been taken on the bottom half of the combustor were copied to the top half and presented in the radial range of: $r/R\in (0,1.3)$ .

Figure 2.

Validation of the present computations. Non-swirling inflow ( $Sn=0$ ) time-averaged radial profiles of: (a–c) non-dimensional axial velocity, (d–f) temperature, (g–i) CO mass fraction and (j–l) product mass fraction. Experiment (Owen et al. Reference Owen, Spadaccini and Bowman1976) ●, present computation ——. Velocity root mean square (RMS): experiment $\times$ (blue), computation — (blue).

The computed velocity profiles are in very good agreement with the measurements. As for the temperature, measured at the first two downstream stations, near the axis, our computation over-predicts the temperature by 300–400 K. This discrepancy in the upstream locations close to the inlet may also result from our combustion model which uses the non-premixed flamelet library and it may be that if one employs a partially premixed model, one would obtain a better fit of the computational results with the experimental data at upstream locations. However, the computed temperature profile is similar to the experimental profile. Furthermore, agreement at the external region ( $r/R$ close to unity) is significantly better than near the axis.

In addition, as one may note, the measurements presented in the experimental study are scattered and hence cause some limitations to our model validation. However, as can be noted, in most cases our computational results appear to fall between the scattered experimental values; i.e. between the results measured in the upper radial location ( $+r$ ) of the combustor and the lower radial location ( $-r$ ) of the combustor, in this axisymmetric configuration. As will be shown here in § 5, large-scale motions of local quenched zones dominate the mixing region, travelling at low frequencies (below 2 Hz). Downstream, however, the flow stabilizes, and consequently, the agreement between the computation and experiment improves. Product mass fraction profiles are predicted well at all stations. CO mole fraction computational value is generally lower than the measured values at all four measurement stations. This is also the case in the computation of Pierce & Moin (Reference Pierce and Moin2004). Shunn & Moin (Reference Shunn and Moin2009) reported that using a more detailed model for the chemistry, which included convective heat loss, the prediction of CO mole fraction was improved. However, the use of the improved model only had a minor effect on the temperature field. The purpose of the current effort is to characterize the slowly evolving turbulent structures, which is accomplished through rigorous tracking of the temperature field. Thus, the present chemical model can be adequately applied here.

4.2 Swirl effect

Time-averaged radial profiles of the thermo-chemical variables computed are shown in figure 3. The effect of swirl on the axial velocity (mean and RMS) is reconstructed very well at all downstream locations. Still, a local disagreement exists at $x/R=0.14$ at the peak velocity. Downstream, at locations $x/R=0.76,x/R=1.44$ , the swirl effect is twofold: higher axial velocities in the regions close to the combustor axis and a shift in the location of maximum axial velocity toward the combustor walls. This phenomenon is captured well by the present computation.

Figure 3.

Swirling inflow ( $Sn=0.3$ ) time-averaged radial profiles of: (a–c) non-dimensional axial velocity, (d–f) temperature, (g–i) CO mass fraction and (j–l) product mass fraction. Experiment (Owen et al. Reference Owen, Spadaccini and Bowman1976) ●, present computation ——. Velocity RMS: experiment $\times$ (blue), computation — (blue).

The level of agreement for CO mole fractions appears to be significantly better than the one observed in the non-swirling case at all cross-sections. Furthermore, better prediction of temperature profiles at the first two cross-sections is obtained for the flow with the inflow swirl. Product mass fractions are captured very well, except for regions close to the combustor’s axis. This result can be attributed to the fact that the swirl dominates the mixing in the flow field and its effect is well captured by the present model.

6.1 Tracking functions formulation

Due to the complexity of turbulent flame evolution, the instantaneous three-dimensional flame structure is commonly presented graphically (e.g. Schefer Reference Schefer1997; Pierce & Moin Reference Pierce and Moin2004), in the same manner as displayed here in § 5. In the current study, we try to derive a quantitative characterization of the flame large-scale fluctuations. First, two-dimensional cross-sections of instantaneous snapshots of three-dimensional temperature distributions are extracted at each time step. Then, a set of one-dimensional functions is derived. This new approach provides a quantitative insight into the cross-section’s temporal flame evolution. The reduced functional representation lends itself to comparisons between flame evolution at different conditions.

Flame extinction and reignition regions were identified during the whole duration of the simulation. The quenched zone in the mixing region fluctuates in all directions. However, even in the case where no swirl was imposed on the inlet, large-scale motions of the quenched zone are identified in the azimuthal direction, which rotates about $180^{\circ }$ in less than one second (figure 5).

Figure 7.

Isosurfaces of temperature: 2000 K (yellow) and 2300 K (red) at different swirl conditions: (a)  $Sn=0$ , (b)  $Sn=0.15$ and (c)  $Sn=0.3$ . In addition, temperature distributions are shown at two axial downstream cross-sections: $x/R=0.1$ and $x/R=1.0$ .

It is suggested to track the quenched zone dynamics using the following formulation: the cross-sections are divided into sectors of a central angle $\unicode[STIX]{x1D719}$ , with variable circumferential angle $\unicode[STIX]{x1D703}$ , measured counterclockwise from the $x$ -axis to the bisector of each sector. The temperature field is scanned by rotating a sector counterclockwise, using small intervals of $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}=0.01~\text{rad}$ , so that the whole circular cross-section is covered, allowing an overlap during the shifting process. A one-dimensional function is derived,

which tracks the minimum temperature sector. The averaged sector temperature $\bar{T}_{\unicode[STIX]{x1D719}}$ is calculated for each sector of $\unicode[STIX]{x1D719}$ , where

and $\bar{T}=\bar{T}_{\unicode[STIX]{x1D719}}(0,2\unicode[STIX]{x03C0})$ is the average temperature over the entire cross-section.

$\unicode[STIX]{x1D702}$ can also be used to evaluate the cross-section circumferential flame homogeneity as previously used by Dagan, Arad & Tambour (Reference Dagan, Arad and Tambour2011). As the sector rotates from $0$ to $2\unicode[STIX]{x03C0}$ , the computations of $\unicode[STIX]{x1D702}$ are employed to monitor the location of the minimum temperature sector: $\unicode[STIX]{x1D703}(\unicode[STIX]{x1D719},t)$ . Evaluation of cross-sections at various downstream locations, extracted from the current computational data base at multiple time steps, led to the selection of $\unicode[STIX]{x1D719}=60^{\circ }$ , as the most efficient section size for capturing the quenched zone. Furthermore, the dependency of $\unicode[STIX]{x1D703}(\unicode[STIX]{x1D719},t)$ on the size of the sector ( $\unicode[STIX]{x1D719}$ ) was found to be very weak. The parameter $\unicode[STIX]{x1D703}$ exhibits discontinuity for every $2\unicode[STIX]{x03C0}$ cycle. For presentation purposes, this undesirable behaviour can be avoided by the use of a functional representation: $\sin ^{2}(0.5\unicode[STIX]{x1D703})$ . As flame dynamics is the focus of the present study, a clear criterion to differentiate between self-sustained burning zones and extinguished zones is needed. A dichotomic-selection method was used for flame location. Using this criterion, for each computational cell, there is a dichotomic classification: unity (flame region) and zero (no flame). The threshold between burning and quenched zone was set at a temperature of 1800 K as follows.

A flamelet solution which is employed by the current chemical model is presented in figure 8(a). The threshold of 1800 K was chosen to separate between a stable self-sustained branch and a non-stable or quenched flame. By doing so, we essentially define a flame thickness, which is on average rich in fuel, as can be seen in figure 8(b). Using the dichotomic representation, large locally quenched regions are automatically captured at three downstream cross-section stations: $x/R=0.1$ , $x/R=1.0$ and $x/R=2.0$ (figure 9).

Figure 8.

(a) Maximum temperature as a function of the scalar dissipation rate, extracted from the present study flamelet solutions. A typical, three-branch s-curve is shown. The region above the chosen threshold of 1800 K represent the self-sustained, stable regions of the flame. (b) A flamelet solution (of the stable branch) with a dissipation rate of $\unicode[STIX]{x1D705}=1\times 10^{-5}~\text{s}^{-1}$ . The vertical dashed line represents the stoichiometric mixture fraction $Z_{st}$ . The shaded area represents the flame thickness, which in average is rich in fuel with $Z>Z_{st}$ .

Figure 9.

Instantaneous flame distribution based on dichotomic selection at downstream stations: (a)  $x/R=0.1$ , (b)  $x/R=1.0$ and (c)  $x/R=2.0$ . Flame region is displayed in white, computed at $t=3.45~\text{s}$ after ignition. No swirl condition: $Sn=0$ . $\unicode[STIX]{x1D702}(\unicode[STIX]{x1D719},t)$ criterion was applied to identify the quenched zones using $\unicode[STIX]{x1D719}=60^{\circ }$ . The resulting sectors are displayed by red lines.

In addition to the circumferential oscillations, radially pulsating flames were observed. A one-dimensional function that tracks the radial oscillations of the flame was used to analyse it. The temperature is extracted from a ring of width $\unicode[STIX]{x0394}r$ at a radius $r$ indicating the median of the ring. The radial location of maximum temperature of a downstream cross-section is extracted using

where

Here, $\unicode[STIX]{x1D709}$ represents the instantaneous radial temperature inhomogeneity of the flame radial location. The maximum temperature is tracked using a dichotomic tracking function $\unicode[STIX]{x1D701}(t)$ (with a threshold of 1800 K).

6.2 Application of the tracking functions

In order to build a database for the inhomogeneity analysis, the non-swirling and swirling cases, discussed in the previous sections, were solved to simulate 17.5 seconds of self-sustained combustion process. In the simulation, 770 000 time steps of $\unicode[STIX]{x0394}t=0.01R/U$ (a total non-dimensional time of $7700$ ) were computed. For the flame tracking analysis, 77 000 samples of instantaneous temperature distributions were collected every tenth time step. Samples were extracted at cross-sections for three downstream locations along the combustion chamber’s axis: $x/R=0.1$ , $x/R=1.0$ , $x/R=2.0$ . The functions $\unicode[STIX]{x1D703}(\unicode[STIX]{x1D719},t)$ and $\unicode[STIX]{x1D701}(\unicode[STIX]{x1D719},t)$ were calculated over each cross-section. It was found that the circumferential level of homogeneity $\unicode[STIX]{x1D702}$ is kept constant throughout the simulation for the three cross-sections in the mixing region. The preservation of inhomogeneity level was reported by Dagan et al. (Reference Dagan, Arad and Tambour2011). Monitoring the dynamics of the circumferential flame location, the flame oscillates in a disordered manner into the non-burning zones, while additional burning zones quench. This phenomenon is captured in figure 10, which shows the temporal development of $\unicode[STIX]{x1D703}$ . The values of $\unicode[STIX]{x1D703}$ were low-pass filtered, by implementing a local regression using weighted linear least-squares for clarity. This method was chosen to preserve low frequencies as they reflect the spatial downstream evolution of the large-scale quenched zones which rotate downstream about the axis, and are the main focus of the current study.

Figure 10.

Temporal development of the quenched region azimuth $\unicode[STIX]{x1D703}$ at two downstream stations: ——, $x/R=0.1$ and - - - -, $x/R=1.0$ for (a) non-swirling conditions, (b) low swirl of $Sn=0.15$ . Due to the periodic nature of the swirling case only a short time interval is presented. Time is measured from ignition.

Figure 11.

Power spectral density of low-pass filtered $\unicode[STIX]{x1D703}(t)$ , extracted at cross-section station $x/R=1.0$ . Inflow swirl condition: $Sn=0.15$ . The most energetic frequency peak is denoted by $k_{p}$ . $k_{s}$ is the wavenumber representing swirl enforcement frequency. The red line represents the turbulent inertial slope of $k^{-5/3}$ .

Figure 12.

A comparison between $\unicode[STIX]{x1D703}(t)$ spectra at different swirl conditions. $Sn=0$ (a–c), $Sn=0.075$ (d–f) and $Sn=0.15$ (g–i) at three downstream stations. $x/R=0.1$ (a,d,g), $x/R=1.0$ (b,e,h) and $x/R=2.0$ (c,f,i). Dominant frequency is indicated in the figure by $f_{p}$ .

It should be mentioned that the low-pass filter does not affect the results of the low flame frequencies. The low-pass filter is used here only to provide a clearer picture of the low frequencies.

For the non-swirling case, irregular dynamics of the quenched zone is clearly visible at the two stations of $x/R=0.1,1.0$ , in the proximity of the mixing region. The difference between the two cross-sections is negligible. This implies that instabilities originating near the inlet create a locally quenched zone which rotates. That way, the quenched zone along the combustor axis exhibits a solid-body rotation, preserving the same angle at different axial locations with similar temporal variations.

Flame dynamics in the case of a swirling inflow is very different than the non-swirling case discussed above. When swirl is imposed on the inlet, the flame moves periodically in a stable manner. The locally quenched zones, which resemble in shape and size those observed in the non-swirling case, swirl at a nearly constant rate (figure 10 b).

A power spectral density plot of the low-pass filtered $\unicode[STIX]{x1D703}(t)$ is shown in figure 11 for the case of swirling inflow $(Sn=0.15)$ . Here, $k=\unicode[STIX]{x1D714}/U_{0}$ is the wavenumber, where $\unicode[STIX]{x1D714}$ is the frequency ( $\unicode[STIX]{x1D714}=2\unicode[STIX]{x03C0}f$ ) and $U_{0}$ is a characteristic velocity, which is taken as the air inlet velocity. At the region of high wavenumber values, the slope of this spectrum is in agreement with the $k^{-5/3}$ slope that characterizes the turbulent inertial subrange. This indicates that the turbulent nature of the combustion is well captured by $\unicode[STIX]{x1D703}(t)$ . Several dominant low frequencies are observed. The most energetic one in the case of $Sn=0.15$ is $k_{p}=1.2$ . The wavenumber that represents the forcing swirl for this case is higher: $k_{s}=5.7$ . Hence, it appears that low-frequency modes of the circumferential flame dynamics are excited in the swirling cases.

It was found that the dominant frequency ( $f_{p}$ ) is linearly dependent on the air velocity magnitude $v$ , diverted by the angle $\unicode[STIX]{x1D6FD}$ , where $v=u/\text{cos}(\unicode[STIX]{x1D6FD})$ . A high linear correlation was found between $f_{p}$ and $v$ (figure 13). This correlation suggests that the rotation mode of the quench zone is directly influenced by the imposed inflow velocity at swirling conditions. This conclusion is consistent with the experimental results reported by Cassidy & Falvey (Reference Cassidy and Falvey1970) and appears also in Gupta et al. (Reference Gupta, Lilley and Syred1984), whereby for high Reynolds numbers, a frequency parameter $fD^{3}/Q$ was found to be independent of the Reynolds number and to increase linearly with the swirl number for swirl numbers higher than $0.3$ (where $f,D$ and $Q$ are the frequency of a PVC developed in a swirled flow, a reference diameter and the volumetric flow rate, respectively). However, for swirl numbers lower than $0.3$ , $fD^{3}/Q$ appears to be constant, and that is the range of swirl numbers that we are dealing with. Thus, for this regime, the frequency $f$ is a linear function of the volumetric flow rate $Q$ , and since $Q$ is directly proportional to the value of the inflow velocity $v$ , then the frequency $f$ should change linearly with $v$ at low swirl numbers, as we found in our study (see figure 13).

The swirling flows have significantly different spectra than the spectra of non-swirling flows. Furthermore, it is possible to characterize the influence of swirl on flame evolution for each swirling case. Focusing on the dominant frequencies of each case, a comparison between the spectra of three different low inflow swirl conditions: $Sn=0,0.075$ and $0.15$ is shown in figure 12. In the non-swirling inflow case, some dominant low frequencies (below 0.5 Hz) are evident. These frequencies do not change with the cross-section station downstream.

In addition to the fact that each of the swirling cases’ spectrum is dominated by a single frequency, it is less oscillatory than the non-swirling case’s spectrum. This insight is in agreement with the premise that the non-swirling flame evolution is a more chaotic process, as discussed in the previous section on time domain analysis.

Figure 13.

Dominant computed frequency $f_{p}$ as a function of the diverted air velocity $v$ at swirl conditions: $Sn=0.075,0.15,0.3$ , respectively. The dashed line represents linear curve fitting between data points.

Figure 14.

Power spectral density of the radial location of maximum temperature: $\unicode[STIX]{x1D701}$ , extracted at three downstream stations $x/R=0.1$ (a,d), $x/R=1.0$ (b,e), $x/R=2.0$ (c,f). $Sn=0$ (a–c), $Sn=0.15$ (d–f) (note that the abscissa and the ordinate in (a–c) are of a different scale than those in (d–f)).

Figure 15.

Correlation between the dominant frequencies of circumferential and radial oscillations of the flame, calculated at cross-section station: $x/R=0.1$ . Two inflow swirl conditions were tested: $Sn=0$ (blue) and $Sn=0.15$ (red). The dashed black line represent identity between circumferential and radial frequencies.

By tracking $\unicode[STIX]{x1D701}(t)$ it was possible to find the flame radial location and extract the pulsating characteristics of the flame. The pulsating motion of the flame starts as a sudden expansion outwards, thus rapidly consuming more of the oxidizer. As the stoichiometric ratio in this region is lost, the flame shrinks and the contact area of the fuel with the oxidizer reduces. This process keeps repeating itself. This flame instability is assumed to be triggered by the turbulent motion of the shear layer between the gaseous fuel and oxidizer, which modifies significantly the contact area.

For all the swirl levels that were investigated, radial fluctuation amplitudes were higher downstream, while minor fluctuations occur at the cross-section station in the mixing region ( $x/R=0.1$ ). On the other hand, large-scale circumferential fluctuations appear to be damped downstream (figure 12). Hence, it is suggested that circumferential oscillations are mostly affected by upstream fluctuations originating at the mixing region close to the inlet while radial fluctuation are excited mostly downstream.

Using the same practice as applied for the circumferential analysis discussed in the previous section, the spectra of radial flame locations for different swirling inflow conditions were computed at different cross-section stations (figure 14). For non-swirling conditions, a concentration of energy at low frequencies was observed. However, for swirl flow, it seems that this energy is dispersed upon a set of higher, distinct frequencies.

Figure 16.

Spectrum of the POD temperature modes.

Figure 17.

First five most energetic POD modes of temperature fluctuations. Modes are extracted for the non-swirling case ( $Sn=0$ ) at downstream location: $x/R=2.0$ .

It appears that the broad band of frequencies in the circumferential direction leaves its signature on fluctuation in the radial direction.

Analysing one-dimensional functions facilitates unique comparisons between different aspects of three-dimensional flame evolution. Such a comparison has been conducted between the flame’s circumferential fluctuations spectrum and the flame’s radial fluctuations spectrum using $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D701}$ spectra, respectively. While the spectra in the different directions (circumferential and radial) are different, a closer look at the most dominant radial and circumferential frequencies reveals a significant correlation which is presented in figure 15. Here, the $x$ -axis represents the circumferential leading frequencies and the $y$ -axis represents the radial dominant frequencies. Seven pairs of frequencies were collected for the non-swirling case (in blue) and seven for the swirling case (in red). The pairs were selected according to the most energetic frequencies (high power spectral density). Pearson coefficients were calculated for each case. The coefficients were $0.952$ for the non-swirling case and $0.9763$ for the swirling case, which indicates a high correlation between the pairs compared. The very high correlation of dominant frequencies indicates an interaction and transfer of energy between fluctuations in different directions. In the following sections, it will be shown that simplified representations of the flow field via proper orthogonal decomposition and linear stability analyses can reveal the onset of fluctuations and energy transfer between the circumferential and radial directions.

6.3 Proper orthogonal decomposition

A POD was utilized to further investigate the coherent flame structures. The POD analysis is used to extract low-order representations of a complex flow field by decomposing temporal, or spatial snapshot of flow field data into a set of orthogonal eigenfunctions. Commonly, POD is applied over the velocity field (Oberleithner et al. Reference Oberleithner, Sieber, Nayeri, Paschereit, Petz, Hege, Noack and Wygnanski2011). However, for the analysis of turbulent combustion, See & Ihme (Reference See and Ihme2014) found the temperature to be more informative. This approach is followed in the present study.

The eigenfunctions $\unicode[STIX]{x1D6E9}_{i}$ are calculated using the minimization of the mean square error between the temperature fluctuations $T(t,\boldsymbol{x})-\bar{T}(\boldsymbol{x})$ and the POD modes as follows:

The snapshot method (Sirovich Reference Sirovich1987) is applied here which extracts the most energetic structures and represents them as linear combinations. Since the large-scale turbulent structures usually contain the major portion of the energy, the POD modes with high-energy content are expected to span the basis for the dominant coherent structures.

The POD is employed on the computed data of the non-swirling case at three downstream cross-section locations: $x/R=0.1,1.0$ and $2.0$ . A set of 770 instantaneous temperature snapshots with a time increment of $0.0227~\text{s}$ were used for the POD analysis. Since the picture was very similar at all the analysed stations, we focus here on the most downstream cross-section location: $x/R=2.0$ . The most energetic modes are presented in figure 16. Over 50 % of the energy is contained in the first 5 most energetic modes. These modes are shown in figure 17. The first two modes contain over $35\,\%$ of the energy in the system and are also similar in shape and have a phase shift of $\unicode[STIX]{x03C0}/2$ , which clearly indicates a travelling azimuthal wave with a wavenumber of $m=1$ . Modes 3, 4 and 5 have twice the azimuthal wavenumber, having a phase shift of $\unicode[STIX]{x03C0}/4$ between them, therefore describing a travelling wave with a wavenumber of $m=2$ . Very similar modes have been identified by Oberleithner et al. (Reference Oberleithner, Sieber, Nayeri, Paschereit, Petz, Hege, Noack and Wygnanski2011) for swirling jets. This finding supports the claim that was raised in § 5: swirling-like rotations characterize the flame behaviour even in cases with zero inflow swirl. The fact that most of the energy is contained only in five simple azimuthal modes suggests that linear stability methods may be applied to investigate the large-scale dynamics of the flame.

7.1 Model formulation

In order to formulate the stability equations, the velocity is assumed to consist of a mean part and an infinitesimally small perturbation part, denoted here by $()^{\prime }$ :

where $U$ , $V$ and $W$ are the axial, radial and tangential base-flow velocities, respectively. Following Olendraru & Sellier (Reference Olendraru and Sellier2002), the spectral collocation method for swirling base flows is described here, in which a ‘quasi-parallel’ base flow is assumed. The linear stability analysis is employed by superimposing on the base flow the disturbances in a modal decomposition form

where $H$ , $F$ , $G$ and $P$ are the eigenfunctions, $k=k_{r}+\text{i}k_{i}$ is the complex axial wavenumber, $m$ is an integer azimuthal wavenumber and $w=w_{r}+iw_{i}$ is the complex frequency. By substituting the modal decomposition into the linearized Navier–Stokes equations and eliminating second-order terms we obtain the following linear system of equations:

Here, $()^{\prime }$ denotes derivation by $r$ .

Assuming no-slip boundary conditions at the walls,

On the centreline, the requirement of single-valued perturbations imposes the following boundary conditions (Olendraru & Sellier Reference Olendraru and Sellier2002):

For the temporal stability analysis, the system of equations is written in as a generalized eigenvalue problem

for a discretized vector $\boldsymbol{X}=\boldsymbol{X}(\unicode[STIX]{x1D701})$ . For $N$ collocation points, matrices $\unicode[STIX]{x1D63C}$ and $\unicode[STIX]{x1D63D}$ are of the size $4\boldsymbol{N}\times 4\boldsymbol{N}$ . Using the Chebyshev derivative matrix $\unicode[STIX]{x1D63F}$ , we derive the expressions for matrices $\unicode[STIX]{x1D63C}$ and $\unicode[STIX]{x1D63D}$ , described in details in appendix B. A Chebyshev collocation approach was used in the present study, which simplifies the treatment of boundary conditions and distribute the numerical error in a homogeneous manner. For the case of concentric jets, $r$ varies between zero and $r_{max}=1$ . Thus, a transformation

is used, which maps the Chebyshev interval $-1\leqslant \unicode[STIX]{x1D701}\leqslant +1$ onto the physical domain $0\leqslant r\leqslant 1$ .

Figure 18.

Time-averaged axial and circumferential base-flow velocity analytical profiles derived from LES results of the present study for a swirling case of $Sn=0.15$ .

Figure 19.

A typical eigenvalue map showing an unstable mode ( $\unicode[STIX]{x1D714}_{i}>0$ ). Calculated with circumferential wavenumber $m=1$ , axial wavenumber $k=3.6$ and two swirl parameters: $q=0$ (empty circles) and $q=0.65$ (filled circles), corresponding to swirl number of $Sn=0$ and $0.15$ in the simulation, respectively.

Figure 20.

Growth rate represented as a function of the axial wavenumber $k$ for various helical modes $m$ for the non-swirling case ( $q=0$ ). The two most unstable modes are of wavenumbers $m=0$ and $m=1$ , denoted in $-\times -$ and dashed lines, respectively. Higher wavenumbers are presented in solid lines.

Figure 21.

(a) Power spectral density of the radial flame fluctuations obtained from simulation (evaluated at $x/R=1.27$ ). In red bullets, dimensionless frequencies obtained from the linear stability analysis $(k=3.6,m=0,1,2,3,4)$ ; (b) correlations between the dimensionless frequencies obtained from simulation $(\unicode[STIX]{x1D714}_{LES})$ and the dimensionless frequencies obtained from linear stability analysis $(\unicode[STIX]{x1D714}_{LS})$ for $m=0,1,2,3,4$ and $k=3.6$ .

Figure 22.

Variation of the least stable mode with swirl parameter $q$ . For this case, $k=3.6$ and $m=1$ , corresponding to maximum growth rate (figure 20).

Figure 23.

Variation of: (a) axial and (b) radial eigenfunctions profiles with the swirl parameter $q$ . Eigenfunctions were calculated for $Re_{s}=20\,000,m=5$ and $k=2$ . The case of $q=0.65$ correspond to swirl number $Sn=0.15$ in simulation. A linear connection between the eigenvalue peaks is shown in green line for $q>0$ .

7.2 Stability analysis of two concentric jets with back flow

The stability analysis is performed here for a case of two concentric jets shear flow with significant back flow, which represents the current study configuration. Unlike the full turbulent reacting case, a non-reacting, constant density, developed base flow is assumed. For the base flow, time-averaged axial and circumferential velocity profiles are extracted from the simulation, from which an analytical function is derived using least-squares functional fit as follows. For the axial velocity,

where $a_{0}=1.4,b=13,c=1.185,d=0.8$ and $a_{1}=0.2$ are constant coefficients. The circumferential velocity profile of the Batchelor vortex case § B.2,

is used here as it represents a similar velocity curve, compared with the computed data for swirl parameters $q<2$ . This is a common approach to evaluate the least stable modes of swirling as well as non-swirling flows from experimental data (e.g. Oberleithner et al. Reference Oberleithner, Sieber, Nayeri, Paschereit, Petz, Hege, Noack and Wygnanski2011). It should be noted that wavenumber growth rates might be weakly dependent on the location of the selected velocity profile. The purpose of this stability analysis is to gain physical insight into mechanism behind the onset of shear flow instabilities in the current configuration. Hence, a profile close to the inlet is considered. Profiles are shown in figure 18. The radial velocity profile is assumed to be zero. The velocity curves are extracted at axial downstream location: $x/R=1.27$ (see § 4, figure 2).

Since two concentric jets are considered, the characteristic velocity for the Reynolds number is defined as the velocity difference between the two jets $\unicode[STIX]{x0394}U$ . The Reynolds number for the stability analysis is therefore

where $R$ is the outer cylinder radius, $\unicode[STIX]{x1D708}$ is the kinematic viscosity.

The eigenvalue map for both the non-swirling and swirling jet is shown in figure 19. One unstable eigenvalue ( $\unicode[STIX]{x1D714}_{i}>0$ ) was found. The growth rate for the non-swirling base-flow $(q=0)$ is shown in figure 20 for several circumferential wavenumbers as a function of the axial wavenumber $k$ . The maximum growth rate is obtained with wavenumber $m=1$ . For wavenumber $m=0$ , a slightly lower growth rate is obtained. As expected, these two modes are the least stable for most axial wavenumbers. It is obvious that the most significant point in the growth rate function is at the peak of the circumferential wavelength $m=1$ , which corresponds to the axial wavenumber of approximately $k=3.6$ (see figure 20). So by choosing the value of $k=3.6$ , we extracted the dimensionless frequency from our analytical stability analysis at $m=0,1,2,3,4$ . Next, by drawing the power spectral density for the radial flame fluctuations, one obtains figure 21(a). Then, by adding to figure 21(a) the points of radial frequencies obtained by our analytical stability calculations, one can clearly see that they fall in a very close vicinity of the computational peaks. This information appears also in figure 21(b), where we added wavenumbers $k=3$ and $k=4$ to show that at $k=3.6$ (about the maximum growth rate), the best match is obtained (with a Pearson correlation coefficient of 0.9980). It should be noted that some of the most dominant modes were not captured by this analysis, and specifically the most dominant mode: $\unicode[STIX]{x1D714}~0.6$ . Though simplified, these analytical results support the transfer of energy from circumferential and axial fluctuations to the radial fluctuations, giving insight into the mechanism of the onset of radial flame fluctuations in the real combustion chamber.

When swirl is considered ( $q>0$ ), the real wave frequency $\unicode[STIX]{x1D714}_{r}$ increases linearly with increase of the swirl parameter $q$ (figure 22 a). This is shown in figure 22 for $k=3.6$ and $m=1$ , corresponding to the least stable spatial mode of the non-swirling case. This is in agreement with the linear correlation found here in simulation (figure 13). On the other hand, the imaginary part of the frequency $\unicode[STIX]{x1D714}$ exhibits a nonlinear behaviour, sharply increasing in lower swirl parameters and reaching a plateau around $q=1$ (figure 22 b). A similar behaviour was observed in our validation test case with a Batchelor vortex (see appendix B for more details).

For non-swirling conditions ( $q=0$ ), axial and radial eigenfunctions are excited. These excitations occur for each test case studied here. However, the most significant effect on the radial fluctuations was found for the following combination of flow conditions: $m=5$ , $k=2$ and $Re_{s}=20\,000$ , which is shown in figure 23. In this specific case, the differences between non-swirling and swirling eigenfunctions are clear: for the axial eigenfunction, the non-swirling curve is diffused, while the swirling cases’ curves are narrow. A double-peak shape exists for the radial non-swirling case which is not evident when swirl is imposed. On the other hand, for increasing swirl parameters ( $q>0$ ), the eigenfunction peak magnitude is linearly increased, both in the axial and radial directions (figure 23). Since we introduced circumferential perturbations in the stability analysis which resulted in radial excitation of modes, especially in the shear layer region, one may expect that these energetic modes will be correlated with the circumferential modes. This phenomenon can explain the correlation between circumferential and radial fluctuations that was found in our computation of the turbulent case, where radial fluctuations were excited for increasing swirl numbers. Moreover, the most energetic frequencies excited are low compared with the typical turbulent inertial subrange frequencies (non-dimensional frequencies of the order of one). We may therefore suggest that the driving mechanism behind the frequency excitation is essentially laminar, as in the stability analysis performed here.

The use of different circumferential perturbation modes have also affected the eigenfunctions shape (figure 23). A double-peak eigenfunction was observed with $m=5$ in the radial eigenfunction at $r=0.82$ . In this case, the introduced swirl significantly changes the radial eigenfunction shape, and moves the axial velocity peak from $0.8$ , where the shear region is, towards the centre at $r=0.45$ . Radial fluctuations therefore affect the mixing between the two jets, especially in non-premixed combustion, where the chemical reaction is governed by mixing between fuel and oxidizer. We therefore expect that it also affects the flame’s evolution. Nevertheless, it should be mentioned that in order to accurately assess the effect of such a phenomenon, a full reactive flow stability analysis is required which is beyond the scope of our study.