1 Introduction
Transitional regimes in fluid mechanics can be characterized as a function of time or spatial location by the passage from an initial state to a different state where each state has specific features. These phenomenon are widely observed in nature as well as in technological applications. One such example is the formation of shock waves, which follows a transient as a function of the Mach number where the acoustic waves propagating in the flow steepen to generate the shock: the initial state corresponds to a low Mach number while the transient ends when the shock is established. The boundary layer transition from laminar to turbulent also presents the characteristics of a transient between an initial and final state. Duct flows exhibit such behaviour as well as a function of the Reynolds number where the flow transitions from laminar to turbulent. In combustion dynamics, when the acoustic, the flow and the flame are coupled, transition from stable state to unstable state is a fundamental research topic. This article presents an approach to investigate this phenomenon with dynamic mode decomposition (DMD).
Reacting fluid dynamics has a central role for engineering and sciences, see the books of Poinsot & Veynante (Reference Poinsot and Veynante2011), Lieuwen (Reference Lieuwen2012), Kuo & Acharya (Reference Kuo and Acharya2012). In energy and aeronautics, most applications utilize combustion to generate power (gas turbines) or propulsion (jet engines). Swirling flames are widely studied in reacting fluid dynamics because of their significance in energy generation and aerospace propulsion applications. Gas turbine and jet engines combustors operate with swirling flames. The development of these applications induce many challenges that warrant scientific analysis. Among these analyses, combustion dynamics focuses on the interaction of waves, aerodynamic structures and the combustion processes inside the combustion chamber. Due to these challenges as well as fundamental interest, swirling flame dynamics has been extensively investigated, see the reviews by Huang & Yang (Reference Huang and Yang2009) and Candel et al. (Reference Candel, Durox, Schuller, Bourgouin and Moeck2014). Future advanced propulsion and power systems may operate in lean premixed mode, which occurs when fuel and air are homogeneously mixed to form a premixture to achieve lower pollutant emissions. The lean premixed mode is sensitive to combustion dynamics, see Candel (Reference Candel2002), Poinsot & Veynante (Reference Poinsot and Veynante2011) and Lieuwen (Reference Lieuwen2012). Combustion instability stems from the resonant closed loop involving the flow field, the flame region and the combustor acoustics. When an unsteady heat release (UHR) source is generated due to flow structures or waves, it induces an unsteady pressure field that feeds back to the initial flow disturbances through the combustor cavity acoustics and the mode conversion processes. The mechanisms by which flames respond to flow oscillations are of primary importance for unsteady heat release modelling. The present article is focused on understanding the physical mechanism driving combustion dynamics in a low swirl injector configuration during the transition to instability and at the limit cycle.
The key mechanisms for premixed swirling flames featuring a recirculating zone of hot combustion products, have been the subject of a series of articles by Palies (Reference Palies2010). Palies et al. (Reference Palies, Durox, Schuller and Candel2010) experimentally investigated the responses of these flames to incident acoustically generated flow perturbations by characterizing them with transfer functions. The gain of those transfer functions exhibited strong differences at specific frequencies and the analysis of the flow field and the flame dynamics was carried out. It was demonstrated analytically that the flame response is driven by two major mechanisms inducing the unsteady heat release: (i) the swirl number oscillation mechanism, which has been initially studied in Richards & Yip (Reference Richards and Yip1999), Straub & Richards (Reference Straub and Richards1999) and Komarek & Polifke (Reference Komarek and Polifke2010), due to the unsteady velocity field affecting the flame front by changing the flame base angle and (ii) the vortex shedding impacting the upper flame region. The mode conversion process inducing the swirl number oscillation was further analysed with numerical simulations, experimental measurements and analytical work in Palies et al. (Reference Palies, Durox, Schuller and Candel2011a ). It was shown that an acoustic wave impinging on an airfoil cascade induces a transverse in two dimensions (or azimuthal in three dimensions with a swirler) convective velocity oscillation on its downstream side. Numerical simulations were also performed to investigate the flow field of the experiment and included an upstream manifold, an injector and a flame tube, see Palies et al. (Reference Palies, Schuller, Durox, Gicquel and Candel2011d ). In this article, the mechanisms generating the unsteady heat release were retrieved with the simulations and further characterized with phase-locked flow visualizations. Based on the identification of the combustion dynamics mechanisms at work, an analytical model of the response of turbulent premixed swirling flames submitted to upstream flow modulation was derived with the G-equation. An explicit expression for this response has been published in Palies et al. (Reference Palies, Schuller, Durox and Candel2011c ). The modelling of combustion instability, i.e. the determination of frequencies and amplitudes was analytically achieved by coupling the flame describing function measurement with an acoustic model consisting of a set of three coupled ducts for multiple configurations, see Palies et al. (Reference Palies, Durox, Schuller and Candel2011b ). The dispersion relation describing the system was numerically solved. The results provided the frequency–growth rate trajectories and generally retrieved the self-sustained oscillations experimental measurements.
The results of this series of articles have been confirmed and explored by researchers and published in several articles. Among them, the article of Kim & Santavicca (Reference Kim and Santavicca2013) characterized in detail self-sustained oscillations in a swirl-stabilized combustor as a function of axial swirler location, oscillation frequency and mean nozzle velocity confirming previous findings. Stability maps with Rayleigh index on the
$y$
axis and phasing between the axial and azimuthal velocity fluctuations on the
$x$
axis were shown to well group the limit cycle and stable operating points. In the presence of equivalence ratio oscillation, Kim & Santavicca (Reference Kim and Santavicca2013) demonstrated that the oscillating flow field (azimuthal and acoustic velocity perturbation) controls the global flame dynamics in association with the flame–coherent structure interactions, depending on the phase difference between the two inlet perturbations. The article of Caux-Brisebois et al. (Reference Caux-Brisebois, Steinberg, Arndt and Meier2014) measured swirl number oscillation in a gas turbine model combustor, confirming previous findings, and demonstrated its impact on the helical vortex core and subsequent unsteady heat release oscillations. Acharya & Lieuwen (Reference Acharya and Lieuwen2015) investigated with numerical simulations premixed swirling flames and confirmed previous findings that the phasing between the azimuthal and axial velocity fluctuation was dependant on the frequency and control the flame response. Furthermore, the fluctuating vorticity field generated by the unsteady flow field and its impact on the flame response was analysed.
These analyses have focused on configurations where an inner recirculation zone stabilizes the flame. Low swirling flames have a different stabilization mechanism and it is a key aspect of the present article to investigate those flames. For the low swirl injector (LSI), phase-locked particle image velocimetry (PIV) and high-speed flame light imaging were carried out to identify combustion dynamics mechanisms at work on thermo-acoustically coupled
$\text{CH}_{4}$
(methane–air) and
$0.9\text{H}_{2}$
–
$0.1\text{CH}_{4}$
(hydrogen–methane–air) flames in Davis et al. (Reference Davis, Therkelsen, Littlejohn and Cheng2013) and Therkelsen et al. (Reference Therkelsen, Enrique Portillo, Littlejohn, Martin and Cheng2013). The results from the two studies showed that vortex shedding plays a significant role in the dynamics of both flames. Linear stability analysis was carried out based on a (
$n$
–
$\unicode[STIX]{x1D70F}$
) model to retrieve the frequency of oscillations. However, the effect of mode conversion at the swirler inducing swirl number fluctuation and the transition to instability were not investigated. The present paper fills this gap by using unsteady phase-resolved flow/flame dynamics analysis and analytical results. The susceptibility of the LSI to the mode conversion process is now evaluated by a non-dimensional analysis. The limit cycle flame dynamics is assumed to behave similarly to the dynamics of a flame submitted to acoustic waves at a frequency close to the self-sustained oscillation frequency to compare Strouhal numbers (
$St$
) of the present and previous work by Palies et al. (Reference Palies, Schuller, Durox and Candel2011c
). The interference length noted
$l_{i}$
between the LSI swirler outlet and the combustor inlet is used to form the Strouhal number based on a convective time of disturbances
$\unicode[STIX]{x1D70F}_{cv}=l_{i}/u_{cv}$
. For both flames being investigated in the LSI, this length is equal to 82 mm while the bulk velocity is
$u_{cv}=18~\text{m}~\text{s}^{-1}$
. Estimate of
$\mathit{St}=f\unicode[STIX]{x1D70F}_{cv}$
gives values of 1.4–1.6 for the LSI. The reference Strouhal number from Palies et al. (Reference Palies, Schuller, Durox and Candel2011c
) is 1.5 with
$l_{i}=43$
mm and
$u_{cv}=2.67~\text{m}~\text{s}^{-1}$
. One can see that the two numbers are close, indicating that the LSI may have a swirl number fluctuation mechanism at work. Following the work of Palies et al. (Reference Palies, Durox, Schuller and Candel2010) the vortex roll-up and swirl number fluctuation mechanisms of the two LSI flames are studied in the present paper. The approach undertaken in the present article is not based on a spatial splitting of the contributing mechanisms as initially carried out in Palies (Reference Palies2010) but on the global responses of the fluid mechanics mechanisms inducing unsteady heat release. A summary of the configurations investigated in the present paper along with previous works is given in figure 1.
An essential aspect of the present paper is the study of transition to instability, which has received little attention in the literature as pointed out by Lieuwen (Reference Lieuwen2002). In the experiment of Broda et al. (Reference Broda, Seo, Santoro, Shirhattikar and Yang1998) and Seo (Reference Seo1999), the inlet temperature of the premixture was shown to impact the transition from stability to instability. Numerical simulations of the experiment performed by Huang & Yang (Reference Huang and Yang2004) confirmed that the flame structure changed at an inlet temperature threshold with a transition to instability. The effect of an increase of the equivalence ratio on transition to instability has been investigated by Gotoda & Tachibana (Reference Gotoda, Nikimoto, Miyano and Tachibana2011), Therkelsen et al. (Reference Therkelsen, Enrique Portillo, Littlejohn, Martin and Cheng2013) and Taamallah et al. (Reference Taamallah, Labry, Shanbhogue and Ghoniem2015). The transition to instability has been investigated by Nagaraja, Kedia & Sujith (Reference Nagaraja, Kedia and Sujith2009) where numerical methods to capture transient growth rates were discussed. The present article attempts to analyse this transition by making use of dynamic mode decomposition directly on experimental flame data where different regimes of oscillation are observed, see Palies et al. (Reference Palies, Davis, Cheng and Ilak2015). DMD is a data-driven postprocessing algorithm for numerical simulation or experimental data. It provides a robust identification of key features when associated with other diagnostics to gain physical insight, as it decomposes the data into a set of spatial modes and corresponding frequencies and growth rates. Dynamic mode decomposition has been used by Bourgoin et al. (Reference Bourgoin, Durox, Moeck, Schuller and Candel2013) to study azimuthal acoustics modes in a system of 16 swirling injectors where DMD was applied directly to the light emission of the flames to identify the nature of the azimuthal mode: rotating or standing. Sayadi et al. (Reference Sayadi, Schmid, Richecoeur and Durox2015) investigated the transient regime of a flame exhibiting oscillation to capture the growth rate. Motheau, Nicoud & Poinsot (Reference Motheau, Nicoud and Poinsot2014) have applied DMD to filter data and capture oscillation frequencies in an aeroengine combustor presenting self-sustained entropy modes. DMD was also applied on the swirling flow field in the work of Labry et al. (Reference Labry, Taamallah, Kewlani, Shanbhogue and Ghoniem2014) where transition from stability to instability was experimentally investigated. The DMD algorithm of Schmid (Reference Schmid2010) is implemented in the present article. Three regimes of oscillations are characterized: stable, local transient to instability and limit cycle, each corresponding to a unique equivalence ratio.
Figure 1. Swirling flame stabilization configurations. The LSI
$\text{CH}_{4}$
flame (left) is stabilized at a distance
$h$
from the injector outlet, see Davis et al. (Reference Davis, Therkelsen, Littlejohn and Cheng2013) and Therkelsen et al. (Reference Therkelsen, Enrique Portillo, Littlejohn, Martin and Cheng2013). The LSI
$\text{H}_{2}$
(centre) flame is stabilized at the injector outlet rim. The
$\text{CH}_{4}$
swirling flame (right) is attached to the centre rod, see Palies et al. (Reference Palies, Durox, Schuller and Candel2010). The injector diameter
$d$
of the LSI is 58 mm and 22 mm for the EM2C configuration on the right.
The present article consists of six sections. The first section outlines the experimental set-up and diagnostics, including the dynamic mode decomposition algorithm. The second section investigates the regimes of oscillations. The third section focuses on flame spectral analysis and flame dynamics making use of ensemble averaging. The fourth and fifth sections contain the analysis of the mechanisms at work, respectively the vortex roll-up and the swirl number oscillation. The last section considers the relative contribution of the two mechanisms on the unsteady heat release.
2 Experimental set-up and diagnostics
2.1 Experimental set-up description
The low swirl injector is an ultra-low emissions technology for gas turbines and industrial heaters, see the reference article of Cheng et al. (Reference Cheng, Yegian, Miyasato, Samuelsen, Pellizzari, Loftus and Benson2000) and Cheng et al. (Reference Cheng, Littlejohn, Nazeer and Smith2008). It uses an aerodynamic flame stabilization method developed for fundamental research on turbulent premixed combustion, see Chan et al. (Reference Chan, Lau, Chin and Cheng2005) and Cheng (Reference Cheng2005). The specific flow field and self-similarity features as well as the flame characteristics have been investigated by laser diagnostics. The analyses showed that the low swirl flame stabilization mechanism is a divergent flow which enables the flame to freely propagate at the location where the local velocity is equal to the turbulent displacement flame speed. This mechanism is different than classical swirl-stabilized systems, where an inner recirculation zone anchors the flame.
The LSI flames being investigated are a
$\text{CH}_{4}$
–air flame of equivalence ratio
$\unicode[STIX]{x1D719}=0.7$
and a set of
$0.9\text{H}_{2}$
–
$0.1\text{CH}_{4}$
–air flames at equivalence ratio
$\unicode[STIX]{x1D719}=0.30,0.38$
and
$0.4$
with an injector bulk flow velocity
$u_{cv}$
equal to
$18~\text{m}~\text{s}^{-1}$
. The equivalence ratio for the mixture of hydrogen and methane with air is defined by:
$\unicode[STIX]{x1D719}=(m_{fuels}/m_{air})/(m_{fuels}/m_{air})_{stoechio}$
. The schematics of the experimental set-up for the LSI can be found in figure 2 of Therkelsen et al. (Reference Therkelsen, Enrique Portillo, Littlejohn, Martin and Cheng2013). The LSI swirler is located in the injector upstream the flame tube (quartz cylinder) of diameter
$d=18$
cm and length
$L=32$
cm. The pressure time trace measured synchronously with the high-speed video is obtained upstream of the swirler with levels of 151 dB at 305 Hz for the
$\text{CH}_{4}$
–air flame and 164 dB at 338 Hz for the
$0.9\text{H}_{2}$
–
$0.1\text{CH}_{4}$
–air flame. In Therkelsen et al. (Reference Therkelsen, Enrique Portillo, Littlejohn, Martin and Cheng2013), the mode shape of the LSI system was determined and is used here to map the acoustic pressure fluctuation (measured at one location) to the flame image data at each pixel location to form Rayleigh index maps.
The video camera equipped with an OH filter captured the flame images snapshots at a sampling frequency
$F_{s}$
of 4000 Hz with a
$800\times 600$
pixel sensor corresponding to
$123\times 92.3$
mm field of view. The relative fluctuation of chemiluminescence signals
$I^{\prime }/\bar{I}$
will be taken as the relative fluctuation of heat release
$\dot{Q}^{\prime }/\bar{\dot{Q}}$
. This is supported by previous experimental work from Price, Hurle & Sugden (Reference Price, Hurle and Sugden1969) and Palies et al. (Reference Palies, Durox, Schuller and Candel2009) and numerical simulations results from Palies et al. (Reference Palies, Schuller, Durox, Gicquel and Candel2011d
) for strictly perfectly premixed flames. High-speed videos and pressure time trace were analysed with 30 cycles of oscillations for the limit cycle study. Phase-locked data were formed with eight points per cycle corresponding to the phases of the PIV measurements documented in Therkelsen et al. (Reference Therkelsen, Enrique Portillo, Littlejohn, Martin and Cheng2013).
2.2 Dynamic mode decomposition
The dynamic mode decomposition algorithm implemented in the present article follows the work of Schmid (Reference Schmid2010). A set of 512 flame snapshots is selected for each operating condition and is organized into two matrices
$\unicode[STIX]{x1D651}_{1}$
and
$\unicode[STIX]{x1D651}_{2}$
such that
$\unicode[STIX]{x1D651}_{1}$
contains the data from instants
$t_{0}$
to
$t_{N-1}$
while
$\unicode[STIX]{x1D651}_{2}$
is made of data from instants
$t_{1}$
to
$t_{N}$
. The derivation of the algorithm is as follows. First, the singular value decomposition of the first matrix
$\unicode[STIX]{x1D651}_{1}$
is computed (
$\unicode[STIX]{x1D650}$
and
$\unicode[STIX]{x1D652}^{H}$
being unitary matrices):
$$\begin{eqnarray}\unicode[STIX]{x1D651}_{1}=\unicode[STIX]{x1D650}\unicode[STIX]{x1D72E}\unicode[STIX]{x1D652}^{H}.\end{eqnarray}$$
As described by Schmid (Reference Schmid2010), the underlying dynamics of the system, whose observables are the snapshot matrices
$\unicode[STIX]{x1D651}_{1}$
and
$\unicode[STIX]{x1D651}_{2}$
, is not necessarily linear. However, for linear dynamics, an operator
$\unicode[STIX]{x1D63C}$
can be defined for which the following mapping holds:
$\unicode[STIX]{x1D63C}\unicode[STIX]{x1D651}_{1}=\unicode[STIX]{x1D651}_{2}$
. An approximation to the spectrum of
$\unicode[STIX]{x1D63C}$
can be calculated using only the flame snapshots and the singular value decomposition of
$\unicode[STIX]{x1D651}_{1}$
, without knowing
$\unicode[STIX]{x1D63C}$
explicitly. This allows the determination of a spectrum associated with the dynamics in the nonlinear regime.
Making use of the previous singular value decomposition (2.1), the following expression is obtained:
$$\begin{eqnarray}\unicode[STIX]{x1D63C}\unicode[STIX]{x1D651}_{1}=\unicode[STIX]{x1D63C}\unicode[STIX]{x1D650}\unicode[STIX]{x1D72E}\unicode[STIX]{x1D652}^{H}=\unicode[STIX]{x1D651}_{2}.\end{eqnarray}$$
This can be rewritten as:
$$\begin{eqnarray}\unicode[STIX]{x1D650}^{H}\unicode[STIX]{x1D63C}\unicode[STIX]{x1D650}\unicode[STIX]{x1D72E}\unicode[STIX]{x1D652}^{H}=\unicode[STIX]{x1D650}^{H}\unicode[STIX]{x1D651}_{2}.\end{eqnarray}$$
The matrix
$\tilde{\unicode[STIX]{x1D64E}}$
is then obtained:
$$\begin{eqnarray}\tilde{\unicode[STIX]{x1D64E}}=\unicode[STIX]{x1D650}^{H}\unicode[STIX]{x1D63C}\unicode[STIX]{x1D650}=\unicode[STIX]{x1D650}^{H}\unicode[STIX]{x1D651}_{2}\unicode[STIX]{x1D652}\unicode[STIX]{x1D72E}^{-1}.\end{eqnarray}$$
The eigenvectors
$\boldsymbol{y}_{i}$
of the matrix
$\tilde{\unicode[STIX]{x1D64E}}$
are calculated and the DMD modes
$\unicode[STIX]{x1D719}_{i}$
are formed, with:
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{i}=\unicode[STIX]{x1D650}\boldsymbol{y}_{i}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D64E}}\boldsymbol{y}_{i}=\unicode[STIX]{x1D707}_{i}\boldsymbol{y}_{i}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D706}_{i}=\log (\unicode[STIX]{x1D707}_{i})/\unicode[STIX]{x1D6FF}t. & \displaystyle\end{eqnarray}$$
Finally, the eigenvalues
$\unicode[STIX]{x1D707}_{i}$
are used to compute the frequencies
$f$
and growth rates
$\unicode[STIX]{x1D70E}$
of the DMD modes:
$f=\text{Im}(\unicode[STIX]{x1D706}_{i}/2\unicode[STIX]{x03C0})$
and
$\unicode[STIX]{x1D70E}=\text{Re}(\unicode[STIX]{x1D706}_{i})$
with
$\unicode[STIX]{x1D6FF}t$
being the inverse of the data sampling frequency
$F_{s}$
. The snapshots used for the analysis in the next section are OH filtered flame images obtained as a function of time so that the captured growth rates correspond to the temporal growth rates.
The algorithm implementation was verified by applying it to a signal with known temporal and spatial characteristics, see Palies et al. (Reference Palies, Davis, Cheng and Ilak2015). The signal taken consisted of a standing wave pattern in a bi-dimensional closed–closed duct with a specific complex pulsation (oscillation frequency and growth rate) and with a superimposed random noise. The DMD algorithm captured the pulsation of the mode and its spatial distribution. The only positive growth rate obtained was the one prescribed in the standing wave signal and the highest amplitude mode that was captured corresponded to the prescribed frequency. Once this step was carried out, the flame image snapshots data were investigated by dynamic mode decomposition.
3 Regimes of oscillations as a function of the equivalence ratio
In this section, the analysis focuses on the hydrogen–methane–air flame at three equivalence ratios
$\unicode[STIX]{x1D719}$
: 0.30, 0.38 and 0.40. The pressure signal measured inside the injector is used to select specific windows within the entire time series. Each selected window corresponds to one equivalence ratio
$\unicode[STIX]{x1D719}$
and to one regime. For these three regimes, dynamic mode decomposition is applied to investigate the flames dynamics. The phase portrait and the fast Fourier transform for the selected windows are presented to characterize the regimes.
3.1 Stable regime
The stable regime is characterized by the time traces of the pressure signal and a selected window, corresponding to the equivalence ratio
$\unicode[STIX]{x1D719}=0.30$
, see figure 2(a). The corresponding integrated unsteady heat release signal is plotted in figure 2(c). The signals have low amplitudes. The pressure phase portrait is computed from the selected signal to characterize the oscillation regime, see figure 2(b). The spectral analysis of the selected pressure signal is achieved by computing the fast Fourier transform, see figure 2(d), and no distinct frequencies are identified.
Figure 2. Stable regime. (a) Pressure signal time traces: entire signal and selected signal (grey coloured). (b) Phase portrait
$f(p)=1/2\unicode[STIX]{x03C0}f\times \unicode[STIX]{x2202}p/\unicode[STIX]{x2202}t$
. (c) Unsteady heat release signal time trace corresponding to the selected window of the pressure signal. (d) Fast Fourier transform of the selected window of the pressure signal time trace.
Figure 3. Stable regime DMD spectra. (a) Amplitude spectrum. (b) Growth rate spectrum. (c) Close-up view of the amplitude spectrum. (d) Close-up view of the growth rate spectrum.
Figure 4. Stable regime DMD modes. (a) Non-oscillating (0 Hz) mode. (b) Mode associated with frequency
$f=335$
Hz in the spectrum.
The dynamic mode decomposition algorithm is applied to the flame images for the stable regime, corresponding to the grey pressure signal in figure 2(a). The frequency and the growth rate spectra are given in figure 3. Panels (a) and (b) correspond to spectra from 0 to 2000 Hz while (c) and (d) correspond to close-up views from 200 to 400 Hz. The results indicate that no dominant frequencies are highlighted in the spectrum. And for all frequencies, the growth rates are negative as no modes are amplified. The mode at 0 Hz has the highest amplitude near a value of unity. The corresponding spatial structure of the mode is given in figure 4(a). A characteristic mode of the spectrum taken at 335 Hz is documented in figure 4(b) for comparison with other regimes.
3.2 Transient regime: from stable to instability regime
Figure 5. Transient regime. (a) Pressure signal time traces: entire signal and selected signal (grey coloured). (b) Phase portrait
$f(p)=1/2\unicode[STIX]{x03C0}f\times \unicode[STIX]{x2202}p/\unicode[STIX]{x2202}t$
. (c) Unsteady heat release signal time trace corresponding to the selected window of the pressure signal. (d) Fast Fourier transform of the selected pressure signal time trace.
The transient regime is characterized by the time traces of the pressure signal and a selected window, corresponding to the equivalence ratio
$\unicode[STIX]{x1D719}=0.38$
, see figure 5(a) and the corresponding unsteady heat release signal is plotted in figure 5(c). The signal has a region of low amplitude, a region of growing amplitude and a region at high amplitude of oscillation. The pressure phase portrait is computed from the signal to characterize the oscillation regime, figure 5(b). The spectral analysis of the selected pressure signal is achieved by computing the fast Fourier transform (FFT) plotted in figure 5(d) and a peak frequency at
$f=328$
Hz is determined.
Figure 6. Transient regime DMD spectra. (a) Amplitude spectrum. (b) Growth rate spectrum. (c) Close-up view of the amplitude spectrum. (d) Close-up view of the growth rate spectrum.
The use of the dynamic mode decomposition on the flame images leads to the spectra plotted in figure 6. The amplitudes (a,c) and the growth rates (b,d) are plotted as a function of the frequency
$f$
. It is observed that the spectrum exhibits a peak frequency at
$f=321.8$
Hz with a growth rate of
$\unicode[STIX]{x1D70E}=5.15~\text{rad}~\text{s}^{-1}$
. This mode is associated with the transient oscillation inside the system. Indeed, this is the only mode close to the FFT pressure signal frequency with a positive growth rate. Other modes in the spectrum near the frequency of oscillations are damped, featuring a negative growth rate.
Figure 7. Transient regime DMD modes. (a) Non-oscillating (0 Hz) mode. (b) Mode associated with main peak in frequency spectrum
$f=321.8$
Hz and
$\unicode[STIX]{x1D70E}=5.15~\text{rad}~\text{s}^{-1}$
.
The spatial modes corresponding respectively to frequencies 0 Hz and 321.8 Hz are plotted in figure 7(a). The 0 Hz mode is slightly different than the stable case. Indeed, the flame is closer to the injector outlet than in the stable regime. This is attributed to the flame speed, which is higher than in the stable case due to the higher equivalence ratio. As a consequence, the flame location changes. The second mode plotted in figure 7(b) is the spatial mode obtained at
$f=321.8$
Hz. This mode is of critical interest as it is the main mode of oscillation for the heat release associated with the transient regime: it has a high amplitude and a positive growth rate. One can also observe the alternating transversal stripes (positive and negative regions) featured by this mode. The distance
$\unicode[STIX]{x1D706}/2$
between these stripes is equal to 20 mm and corresponds to the half of a wavelength. The characteristic velocity of this mode is then given by:
$u=\unicode[STIX]{x1D706}f=12.8~\text{m}~\text{s}^{-1}$
, this value being of the order of magnitude of the bulk velocity in the injector, respectively equal to
$18~\text{m}~\text{s}^{-1}$
. Therefore, the mode has a convective nature associated with a vortical wave.
The entire pressure time trace in figure 5(a) indicates an intermittency behaviour. Intermittency in combustion instability has been identified as a regime preceding limit cycle oscillation, see the articles of Gotoda & Tachibana (Reference Gotoda, Nikimoto, Miyano and Tachibana2011) and Nair & Sujith (Reference Nair, Thampi and Sujith2014). Such signals present local transient to instability characterized by strong pressure oscillation amplitude. The present work aims at isolating three distinct regimes, respectively stable, local transient to instability and limit cycle to determine the features inducing the transition and the mechanisms at work during limit cycle. The vortical wave is associated with the transition to instability and induces unsteady heat release. In addition, the non-oscillating 0 Hz mode provides significant indication on the role of the flame location that will affect the phasing between the heat release and the vortical wave as the latter has essentially a constant speed and constant frequency between the operating points. As the non-oscillating mode (flame location) shifts position between the stable, transient and limit cycle regimes, the unsteady heat release signal will also shift versus the pressure oscillation signal, which can lead to instability. Indeed, the pressure oscillation induced/emitted by the unsteady heat release propagates at a velocity near the sound speed and interacts with the acoustics of the system but weakly contribute to the phase difference between unsteady heat release and pressure. The delay between UHR and pressure in the feedback loop is mainly due to the upstream flame fluid dynamics mechanisms, specifically due to the acoustic–convective mode conversions occurring both at the swirler (swirl number oscillation) and at the dump plane of the combustor (vortical wave) that generate the UHR. This is supported by the following order of magnitude analysis: a leading edge flame location difference
$l=5$
mm will induce a corresponding phase difference equal to
$2\unicode[STIX]{x03C0}f\unicode[STIX]{x1D70F}$
where
$\unicode[STIX]{x1D70F}=l/u$
. The calculation results in a phase of 0.79 rad or
$45.2^{\circ }$
and demonstrates the strong impact of the non-oscillating DMD mode contribution on stability. The leading edge flame location difference for the stable and transient regime can be observe by comparing the non-oscillating DMD modes of figures 4(a) and 7(a). The mechanisms inducing unsteady heat release are further investigated in the next sections of this article.
3.3 Limit cycle regime
Figure 8. Limit cycle regime. (a) Pressure signal time traces: full signal and selected signal (grey coloured). (b) Phase portrait
$f(p)=1/2\unicode[STIX]{x03C0}f\times \unicode[STIX]{x2202}p/\unicode[STIX]{x2202}t$
. (c) Unsteady heat release signal time trace corresponding to the selected window of the pressure signal. (d) Fast Fourier transform of the selected pressure signal time trace.
The limit cycle regime is characterized by the time traces of the pressure signal for a selected window, corresponding to the equivalence ratio
$\unicode[STIX]{x1D719}=0.40$
, see figure 8(a). The corresponding unsteady heat release signal obtained is plotted in figure 8(c). The signal has a high amplitude of oscillation. The pressure phase portrait is computed from the signal to characterize the oscillation regime, figure 8(b). The spectral analysis of the selected pressure signal is carried out by computing the fast Fourier transform, figure 8(d) and shows a peak frequency at
$f=335$
Hz.
Figure 9. Limit cycle regime DMD spectra. (a) Amplitude spectrum. (b) Growth rate spectrum. (c) Close-up view of the amplitude spectrum. (d) Close-up view of the growth rate spectrum.
Figure 10. Limit cycle regime DMD modes. (a) Non-oscillating (0 Hz) mode. (b) Mode associated with the main peak in the frequency spectrum (
$f=340$
Hz).
The frequency and growth rate spectra of the dynamic mode decomposition are given in figure 9. Panels (a) and (b) correspond to the complete spectrum from 0 to 2000 Hz while (c) and (d) spectrum corresponds to close-up views from 200 to 400 Hz. The results indicate that multiple frequencies have peaks in the spectrum reflecting the nonlinear character of the flame. The main peaks in the spectrum correspond to the 0 Hz mode and to the frequency associated with the pressure oscillation at
$f=340$
Hz. One observes a set of higher harmonics of the fundamental frequency. It is interesting to note that the first harmonic has a positive growth rate while the main frequency has a near zero growth rate
$\unicode[STIX]{x1D70E}=-1.3~\text{rad}~\text{s}^{-1}$
as expected for a limit cycle oscillation.
The spatial modes corresponding respectively to frequencies 0 Hz and 340 Hz are plotted in figure 10(a). The 0 Hz mode is slightly different from previous regimes. Because the equivalence ratio is different, the flame location changes. The second mode plotted in figure 10(b) is the spatial mode obtained at
$f=340$
Hz. This mode is the main mode of oscillation for the heat release associated with the limit cycle regime: it has a high amplitude and a near zero growth rate. One observes as well that the mode consists of alternating transversal stripes. The distance
$\unicode[STIX]{x1D706}/2$
between these stripes is equal to 20 mm. The characteristic velocity of this mode is given by:
$u=\unicode[STIX]{x1D706}f=13.4~\text{m}~\text{s}^{-1}$
, this value being of the order of magnitude of the bulk velocity, respectively equal to
$18~\text{m}~\text{s}^{-1}$
. Therefore, the mode has a convective nature associated with a vortical wave.
Self-sustained oscillations which include coupling with entropy waves present a specific characteristic that the entropy wave propagates in a convective manner inside the flame tube up to its outlet where the entropy wave is converted into an upstream propagating acoustics wave feeding the resonant loop, see the article of Motheau et al. (Reference Motheau, Nicoud and Poinsot2014). In the present experimental study, the outlet of the flame tube does not have a boundary allowing such coupling.
3.4 Time average versus non-oscillating (0 Hz) DMD mode
In the previous sections, the stability of the present configuration is associated with the flame front location captured with the 0 Hz mode. The flame front location being determined by the effect of the equivalence ratio changing the flame speed. The present section demonstrates that the nonlinear character of the flame impacts weakly the non-oscillating (0 Hz) mode flame front location while it does impact the time averages. Indeed, the dynamic mode decomposition of a complex variable
$\unicode[STIX]{x1D6FC}(x,y,z,t)$
can take the following form:
$$\begin{eqnarray}\unicode[STIX]{x1D6FC}(x,y,z,t)=\mathop{\sum }_{k=1}^{N-1}\exp (\unicode[STIX]{x1D706}_{k}t)\unicode[STIX]{x1D719}_{k}(x,y,z),\end{eqnarray}$$
where
$N$
is the total number of snapshots,
$\unicode[STIX]{x1D706}_{k}=\unicode[STIX]{x1D70E}_{k}+\text{i}2\unicode[STIX]{x03C0}f_{k}$
its complex pulsation and
$\unicode[STIX]{x1D719}_{k}(x,y,z)$
its complex spatial distribution. One can now apply the time-average operator to this equation:
$$\begin{eqnarray}\overline{\unicode[STIX]{x1D6FC}(x,y,z,t)=\mathop{\sum }_{k=1}^{N-1}\exp (\unicode[STIX]{x1D706}_{k}t)\unicode[STIX]{x1D719}_{k}(x,y,z)},\end{eqnarray}$$
which leads to:
$$\begin{eqnarray}\overline{\unicode[STIX]{x1D6FC}(x,y,z,t)}=\overline{\exp (\unicode[STIX]{x1D70E}_{1}t)}\unicode[STIX]{x1D719}_{1}(x,y,z)+\overline{\mathop{\sum }_{k=2}^{N-1}\exp (\unicode[STIX]{x1D706}_{k}t)\unicode[STIX]{x1D719}_{k}(x,y,z)},\end{eqnarray}$$
which can be also written as:
$$\begin{eqnarray}\bar{\unicode[STIX]{x1D6FC}}(x,y,z)=\frac{\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x1D70E}_{1}T}(\exp (\unicode[STIX]{x1D70E}_{1}T)-1)+\frac{1}{T}\mathop{\sum }_{k=2}^{N-1}\frac{\unicode[STIX]{x1D719}_{k}}{\unicode[STIX]{x1D706}_{k}}(\exp (\unicode[STIX]{x1D706}_{k}T)-1).\end{eqnarray}$$
The time average
$\bar{\unicode[STIX]{x1D6FC}}$
over a range of time from 0 to
$T=N\unicode[STIX]{x1D6FF}t$
is the summation of two contributions: the first contribution is due to the non-oscillating (0 Hz) mode (
$k=1$
for which the frequency
$f_{1}$
is exactly zero and the growth rate
$\unicode[STIX]{x1D70E}_{1}$
is near zero); the second contribution contains the nonlinear modes interactions affecting the time average. This expression shows mathematically the explicit contribution of a static component and a dynamic component. As a consequence, the change in the equivalence ratio is the cause of the change in the 0 Hz mode across the stable, transient and limit cycle regime, as the nonlinear character of the flame is excluded.
4 Flame oscillation content and flame dynamics
The frequency content of the unsteady regions of the flame is characterized by performing fast Fourier transform on the chemiluminescence images to outline periodic motion. To extract frequency content, FFTs of time series of chemiluminescence signal at each pixel
$I(x,y,t)$
are computed to determine the peak frequency and its corresponding amplitude. Results are plotted in figure 11 with amplitudes at the top and frequencies at the bottom. For the methane–air (
$\text{CH}_{4}$
) flame (figure 11
a,c), most regions of the flame oscillate at a frequency close to 305 Hz while the root of the flame has a lower frequency that may be attributed to lift-off events. One can also see that the highest unsteadiness occurs mainly at the flame trailing region. For the
$0.9\text{H}_{2}$
–
$0.1\text{CH}_{4}$
-air flame (figure 11
b,d), the whole flame has a frequency peak at 345 Hz. For both flames, the first harmonic peak of the fundamental frequency (not shown here) is much weaker. The results demonstrate that the entire flame responds most strongly to the frequency of the pressure oscillation in the flame tube.
Figure 11.
$\text{CH}_{4}$
(a,c) and
$\text{H}_{2}$
(b,d) flame oscillation content. Time series
$I(x,y,t)$
of chemiluminescence signals are formed at each pixel before spectral analysis. Results indicate a main frequency of 305 Hz for
$\text{CH}_{4}$
flame and 345 Hz for
$\text{H}_{2}$
flame. (a,b) Map of the signals fast Fourier transform peak amplitudes. (c,d) Associated peak frequency maps where grey colour corresponds to the self-sustained oscillation frequencies. The images’ field of view is: 123
$\times$
92.3 mm.
The flame dynamics is now analysed by inverse Abel transformation to form phase-locked chemiluminescence images, see the articles by Herding et al. (Reference Herding, Snyder, Rolon and Candel1998) and Durox, Schuller & Candel (Reference Durox, Schuller and Candel2005). This mathematical transformation allows the construction of a slice for flame light images integrated over lines of sight. The inverse Abel transformed images at four phases, separated by
$90^{\circ }$
, at the fundamental frequency are shown in figure 12. For the
$\text{CH}_{4}$
flame, one observes that the bowl flame shape evolution is weak and most of the flame motion is confined to the flame trailing region. For the
$0.9\text{H}_{2}$
–
$0.1\text{CH}_{4}$
flame, the evolution is significantly stronger and the vortex flame roll-up within the cycle is quite obvious. At
$0^{\circ }$
, a vortex is formed at the LSI outlet injector while the preceding one is released at the flame trailing region. At
$90^{\circ }$
and
$180^{\circ }$
, the vortex is folding the flame front. At phase
$270^{\circ }$
, the vortex reaches the flame tip.
Figure 12. Inverse Abel transformed phase-locked images for
$\text{CH}_{4}$
(a,c,e,g) and
$\text{H}_{2}$
(b,d,f,h) flames. Phase instants of the cycle are (a,b)
$0^{\circ }$
, (c,d)
$90^{\circ }$
, (e,f)
$180^{\circ }$
and (g,h)
$270^{\circ }$
. The images’ field of view is:
$123\times 92.3$
mm.
5 Vortex roll-up mechanism
To characterize the vortex roll-up mechanism, the mode of propagation of the perturbation is determined by examining the phase of the flame transfer function (FTF) defined in the frequency domain as:
$$\begin{eqnarray}{\mathcal{F}}(\unicode[STIX]{x1D714})=\frac{I^{\prime }(x,y)/\bar{I}(x,y)}{p_{ref}^{\prime }(x_{ref},y_{ref})/\bar{p_{ref}}(x_{ref},y_{ref})}=G\exp (\text{i}\unicode[STIX]{x1D711}),\end{eqnarray}$$
where
$\unicode[STIX]{x1D714}=2\unicode[STIX]{x03C0}f$
,
$G$
is the gain and
$\unicode[STIX]{x1D711}$
is the phase difference between the chemiluminescence intensity
$I(x,y,t)$
and the reference pressure signal. Although the self-sustained oscillating flames have no external forcing,
${\mathcal{F}}(\unicode[STIX]{x1D714})$
is computed from the reference pressure signal. The FTF is a linear concept used to determine the response of flames to upstream weak disturbances leading to a gain and a phase as a function of the excitation frequency. In the present context, the phase of the transfer function at the dominant oscillation frequency is used to retrieve the convective speed of the vortical wave occurring at that frequency. The FTF has been extended to FDF (flame describing function) where the amplitude of the input forcing is not necessarily small but can have significant levels as in the present case. It has been shown that the phase of premixed swirling flames FDF have small variations for multiple input velocity forcing levels whereas the gain has strong variations, see Palies et al. (Reference Palies, Durox, Schuller and Candel2010). As a consequence, the propagation of the disturbance (vortical wave) along the flame generating unsteady heat release is weakly dependent on the oscillation amplitude level.
The chemiluminescence signal
$I(x,y,t)$
is used to identify the mode of propagation of the velocity perturbations. The phase plotted in figure 13 indicates that a convective mode of perturbation is at work for both flames. For the
$0.9\text{H}_{2}$
–
$0.1\text{CH}_{4}$
–air flame, this is specifically noticeable and it is possible to associate a convective time of perturbations along the flame front. Indeed, a phase
$\unicode[STIX]{x1D711}=\unicode[STIX]{x03C0}$
is obtained over a distance
$l=30~\text{mm}$
from figure 13 on the swirling jets. Since
$\unicode[STIX]{x1D711}=\unicode[STIX]{x1D714}\unicode[STIX]{x1D70F}=2\unicode[STIX]{x03C0}f\unicode[STIX]{x1D70F}$
, the time delay
$\unicode[STIX]{x1D70F}$
can be calculated as:
$\unicode[STIX]{x1D70F}=1/(2f)$
. The acoustics frequency being 338 Hz for the
$0.9\text{H}_{2}$
–
$0.1\text{CH}_{4}$
flame, the obtained convective velocity is
$20.6~\text{m}~\text{s}^{-1}$
. This confirms that velocity perturbations are convected along the flame front. While the convective velocity obtained from DMD (
$u_{cv}=13.4~\text{m}~\text{s}^{-1}$
) and from the transfer function phase lead to different values of the convective speed of the vortical mode, its convective nature is shown by the two methods. The difference is attributed to three sources. Firstly, the measurements from the first diagnostics are extracted from the flame transfer function phase result for only two points located on the swirling jet. Secondly, the determination of the velocity from the DMD limit cycle mode does not make use of an edge detection algorithm to detect the wavelength between the two stripes of the mode. Finally, the transfer function formalism decomposes signals into a set of periodic signals. As the DMD showed, the limit cycle presents a negative growth rate (
$\unicode[STIX]{x1D70E}=-1.3~\text{rad}~\text{s}^{-1}$
). This growth rate will have an impact on Fourier analysis results (flame transfer function). DMD and Fourier analysis should lead to the same results for perfectly oscillating signals (Rowley et al.
Reference Rowley, Mezic, Bagheri, Schlatter and Henningson2009) and not for a signal with a positive or negative growth rate as in the present case.
Figure 13.
$\text{CH}_{4}$
(a,c,e) and
$\text{H}_{2}$
(b,d,f) flames Rayleigh index and flame transfer function phase maps. (a,b) Time integrated Rayleigh index map. (c,d) and (e,f) Flame transfer function phase
$\unicode[STIX]{x1D711}$
and phase criterion
${\mathcal{C}}_{\unicode[STIX]{x1D711}}$
respectively. The phase criterion is equal to
$-1$
if
$|\unicode[STIX]{x1D711}|\leqslant 90^{\circ }$
and to 1 if
$|\unicode[STIX]{x1D711}|>90^{\circ }$
. The images’ field of view is:
$123\times 92.3$
mm.
To evaluate the regions where vortex roll-up induces driving regions for the oscillation, the two-dimensional Rayleigh index maps are formed and integrated in time over the whole time series. The Rayleigh index time integral is defined as:
$$\begin{eqnarray}RI(x,y)=\frac{1}{\bar{I}(x,y)\bar{p}(x,y)}\int I^{\prime }(x,y,t)\times p^{\prime }(x,y,t)\,\text{d}t.\end{eqnarray}$$
The quantity
$RI(x,y)$
for the
$\text{CH}_{4}$
flame (figure 13
a) shows that the driving region is at the flame trailing region and enclosed by two damping regions upstream and downstream. For the
$0.9\text{H}_{2}{-}\text{CH}_{4}$
flame, the driving regions correspond to the upper part of the flame, while the damping region is located on the sides of the flame front.
Figure 13(c,d) compares the phase of the flame transfer function with
$RI(x,y)$
maps. While the
$RI(x,y)$
presents the regions of driving and damping, the phase of the flame transfer function indicates the phasing between the reference pressure signal and the chemiluminescence signal and therefore the mode of propagation of perturbations. From these results, one observes that for driving regions, the absolute value of the phase
$\unicode[STIX]{x1D711}$
of the flame transfer function is such that
$|\unicode[STIX]{x1D711}|\leqslant 90^{\circ }$
, so that the pressure oscillation and the heat release oscillation are in phase, while for the damping region, the value of the absolute phase
$\unicode[STIX]{x1D711}$
is such that:
$|\unicode[STIX]{x1D711}|>90^{\circ }$
, meaning that the pressure and heat release oscillations are out of phase, as expected.
6 Swirl number fluctuation mechanism
The methodology for analysing swirl number fluctuation mechanism relies on examining the fluctuation of the swirl number within an ensemble-averaged oscillation cycle and its relationship to unsteady heat release.
The relative fluctuation of swirl number
$S^{\prime }/\bar{S}$
defined in Palies et al. (Reference Palies, Durox, Schuller and Candel2010) writes:
$$\begin{eqnarray}\frac{S^{\prime }}{\bar{S}}=\frac{u_{\unicode[STIX]{x1D703}}^{\prime }}{\bar{u}_{\unicode[STIX]{x1D703}}}-\frac{u_{z}^{\prime }}{\bar{u}_{z}},\end{eqnarray}$$
where
$u_{z}^{\prime }/\bar{u}_{z}$
and
$u_{\unicode[STIX]{x1D703}}^{\prime }/\bar{u}_{\unicode[STIX]{x1D703}}$
are the relative fluctuations of axial and azimuthal velocities. This definition of the swirl number is a linearized form of the classical swirl number definition based on momentum conservation equations. The hypothesis of linear regime is valid for the wave part of the flow (acoustic and vortical) as the main nonlinearities are due to the flame.
The axial and the radial velocity components are measured in Therkelsen et al. (Reference Therkelsen, Enrique Portillo, Littlejohn, Martin and Cheng2013) and the azimuthal and swirl number signals are reconstructed from those measurements based on analytical results. While the phase-locked PIV measured axial
$u_{z}$
and radial
$u_{r}$
velocities components, we use previous results indicating that
$u_{z}^{\prime }/\bar{u}_{z}$
and
$u_{\unicode[STIX]{x1D703}}^{\prime }/\bar{u}_{\unicode[STIX]{x1D703}}$
have similar amplitude, see Palies et al. (Reference Palies, Durox, Schuller and Candel2010). Accordingly, we estimated
$u_{\unicode[STIX]{x1D703}}^{\prime }/\bar{u}_{\unicode[STIX]{x1D703}}$
from the measured axial and radial velocity components by amplitude scaling and using the fact that
$u_{r}^{\prime }$
and
$u_{\unicode[STIX]{x1D703}}^{\prime }$
have the same phase as they propagate both perpendicularly to the axial direction of the oscillation, see appendix A. The scaling is performed by using averaged values obtained at the swirling jet (
$19~\text{mm}<r<31$
mm) at specific locations upstream of the flame front, respectively
$z=30~\text{mm}$
for the methane flame and
$z=3.5$
mm for the methane–hydrogen flame.
Results of this step are plotted on figures 14(b) and 15(b). For the
$\text{CH}_{4}$
flame,
$u_{z}^{\prime }/\bar{u}_{z}$
and
$u_{\unicode[STIX]{x1D703}}^{\prime }/\bar{u}_{\unicode[STIX]{x1D703}}$
amplitudes are weak and out of phase leading to low fluctuations in
$S^{\prime }/\bar{S}$
, see (6.1). For the
$0.9~\text{H}_{2}$
–
$0.1~\text{CH}_{4}$
flame, those velocity signals have higher amplitudes and are evolving with a smaller phase difference that leads to large swirl number fluctuation. Those results indicate that the swirl number is unsteady over the ensemble-averaged cycle. In addition, the phase-resolved divergence rate relative fluctuations
$a_{z}^{\prime }/\bar{a}_{z}$
taken from Therkelsen et al. (Reference Therkelsen, Enrique Portillo, Littlejohn, Martin and Cheng2013) are also given in figures 14 and 15 as an indication of the swirl number. It was shown in Cheng et al. (Reference Cheng, Littlejohn, Nazeer and Smith2008) that there is a dependency of the divergence rate
$a_{z}$
with the swirl number
$S$
. As presented in appendix A of Day et al. (Reference Day, Tachibana, Bell, Lijewski, V. and Cheng2012) the near field divergence rate
$a_{z}=(\text{d}u_{z}/\text{d}z)/U_{0}$
where
$u_{z}$
is the mean axial velocity and
$U_{0}$
the bulk flow velocity is a parameter to characterize the self-similar nature of the near field generated by the LSI. For the
$\text{CH}_{4}$
flame,
$a_{z}^{\prime }/\bar{a}_{z}$
evolves with
$S^{\prime }/\bar{S}$
for most of the phases instants while for the
$0.9\text{H}_{2}$
–
$0.1\text{CH}_{4}$
flame, the discrepancy is significant and is not used in the subsequent analysis.
Figure 14. (a) Relative fluctuations of pressure and chemiluminescence emission signals as a function of the phase instant for the
$\text{CH}_{4}$
flame. Chemiluminescence signals are obtained by spatially integrating the full OH emission or the full inverse Abel emission. (b) Relative fluctuations of axial and azimuthal velocities and swirl.
Figure 15. (a) Relative fluctuations of pressure and chemiluminescence emission signals as a function of the phase instant for the
$0.9\text{H}_{2}$
–
$0.1\text{CH}_{4}$
flame. Chemiluminescence signals are obtained by spatially integrating the full OH emission or the full inverse Abel emission. (b) Relative fluctuations of axial and azimuthal velocities and swirl.
It is now of interest to quantify the effect of the swirl number oscillation on the unsteady heat release. It can be investigated in a non-dimensional form by linking the swirl oscillations to the fluctuations of the turbulent flame speed. This is accomplished by using (6.1) and the definition of the fluctuating ensemble-averaged turbulent flame speed defined in Palies et al. (Reference Palies, Durox, Schuller and Candel2010). The assumption is made that this relationship holds for all type of swirling flames and is supported by the G-equation model presented in Palies et al. (Reference Palies, Schuller, Durox and Candel2011c ):
$$\begin{eqnarray}\frac{S_{T}^{\prime }}{\bar{S}_{T}}=\unicode[STIX]{x1D712}\frac{u_{\unicode[STIX]{x1D703}}^{\prime }}{\bar{u}_{\unicode[STIX]{x1D703}}}+\unicode[STIX]{x1D701}\frac{u_{z}^{\prime }}{\bar{u}_{z}},\end{eqnarray}$$
where
$\unicode[STIX]{x1D712}=-0.4$
and
$\unicode[STIX]{x1D701}=0.4$
. Introducing (6.1) into (6.2), one can show that the turbulent burning velocity is linked to the swirl number fluctuation:
$$\begin{eqnarray}\frac{S_{T}^{\prime }}{\bar{S}_{T}}=\unicode[STIX]{x1D712}\frac{S^{\prime }}{\bar{S}}.\end{eqnarray}$$
The turbulent flame speed is therefore directly dependent on the oscillation of the swirl number, which is important for modelling of unsteady heat release in swirling flames, and will be used in the next section.
7 Relative contributions of vortex roll-up and swirl number fluctuation
By utilizing hypotheses for the flow field and the flame structure, two mechanisms inducing unsteady heat release have been characterized in the previous sections: vortex roll-up and swirl number fluctuation. Within these hypotheses and combining the measured quantities and analytical methods, the contributions of the two mechanisms towards the global unsteady heat release can be calculated. The following equation expresses the heat released by a premixed flame front of surface area
$A$
propagating at a turbulent flame speed
$S_{T}$
:
$$\begin{eqnarray}\dot{Q}=I=\unicode[STIX]{x1D70C}S_{T}A\unicode[STIX]{x0394}hY_{f},\end{eqnarray}$$
where
$\unicode[STIX]{x1D70C}$
is the density of the premixed reactants,
$\unicode[STIX]{x0394}h$
is the heat released per unit mass of fuel during reaction and
$Y_{f}$
is the mass fraction of fuel in the premixture. The two main contributions to the heat release deduced from this equation are due to the flame surface area fluctuation and the turbulent flame velocity fluctuation.
The total unsteady heat release is determined by linearizing (7.1) around a mean state as defined by Palies et al. (Reference Palies, Durox, Schuller and Candel2010). One obtains:
$$\begin{eqnarray}\frac{I^{\prime }}{\bar{I}}\simeq \frac{\dot{Q}^{\prime }}{\bar{\dot{Q}}}=\left[\frac{A^{\prime }}{\bar{A}}+\frac{S_{T}^{\prime }}{\bar{S}_{T}}\right].\end{eqnarray}$$
This indicates that the unsteady heat release depends on a contribution that modulates the flame surface and on a contribution that modulates the turbulent burning velocity. It has been presented in previous sections that the flame surface oscillation is due to the vortex roll-up, while the swirl number oscillation mechanism is directly linked to the turbulent flame speed via (6.3).
Figure 16. Relative fluctuations of chemiluminescence emission signals as a function of the phase instant. Chemiluminescence signals are obtained by spatially integrating the full inverse Abel emission. (a)
$\text{CH}_{4}$
flame. (b)
$0.9\text{H}_{2}$
–
$0.1\text{CH}_{4}$
flame.
The first mechanism, the vortex roll-up, induces a flame surface area modification
$A^{\prime }/\bar{A}$
and is due to the vortical wave propagating along the flame. The evolution of the velocity signal
$u_{z}^{\prime }/\bar{u}_{z}$
with the total unsteady chemiluminescence signal is compared next. For the
$\text{CH}_{4}$
flame, results are plotted in figure 14. The axial velocity signal presents a small amplitude and is in phase opposition with the unsteady chemiluminescence signal. For the
$0.9\text{H}_{2}$
–
$0.1\text{CH}_{4}$
flame in figure 15, the chemiluminescence emission signal follows the axial velocity oscillation whose amplitude is significantly higher than for the
$\text{CH}_{4}$
flame. This confirms previous observations that a strong axial velocity fluctuation occurs upstream the flame front and induces this unsteady heat release. The flame area fluctuation
$A^{\prime }/\bar{A}$
is linked to the vortical wave along the flame front convecting from the base of the flame to its tip. This wave has a radial and an axial component (the vorticity vector being oriented on the azimuthal direction). In the case of the
$\text{CH}_{4}$
flame, the same mechanism of flame area fluctuation is at work but the amplitude is lower than in the case of the
$0.9\text{H}_{2}$
–
$0.1\text{CH}_{4}$
flame.
Figure 17. Schematic of the mechanisms at work during the self-sustained oscillation. The acoustic modulation induces a mode conversion process at the swirler that generate fluctuations of axial and azimuthal velocities. The swirl number is then modulated by this unsteady velocity field and give rise to an unsteady heat release (local time intensity variation) through (6.3). The axial velocity oscillation leads to an unsteady heat release rate through the vortex shedding (local spatial variation). These two contributions determine the total unsteady heat release. Schematic adapted from Palies (Reference Palies2010) and Candel et al. (Reference Candel, Durox, Schuller, Bourgouin and Moeck2014).
The second mechanism is due to the swirl number oscillation
$S^{\prime }/\bar{S}$
and induces an oscillation of the turbulent burning velocity
$S_{T}^{\prime }/\bar{S}_{T}$
as demonstrated by (6.3). One can also show that this mechanism is responsible for unsteady heat release through (7.1). The swirl number depends both on the axial and azimuthal velocities and the vortical disturbances will have an effect on this quantity. However, the swirl number fluctuation is determined upstream of the flame front where the axial velocity fluctuation is mainly due to the mode conversion process from the swirler (convective azimuthal component and acoustic axial component propagating downstream the swirler) where the vortical part is essentially small. The locations of the measurements for the swirl number oscillation are obtained at the swirling jet (
$19~\text{mm}<r<31~\text{mm}$
) at specific locations upstream of the flame front, respectively
$z=30~\text{mm}$
for
$\text{CH}_{4}$
and
$z=3.5~\text{mm}$
for
$0.9\text{H}_{2}$
–
$0.1\text{CH}_{4}$
. To support that point, these locations can be compared to the trajectory of the shed vortices in the article of Therkelsen et al. (Reference Therkelsen, Enrique Portillo, Littlejohn, Martin and Cheng2013). As a consequence, the two processes can be seen as independent. In Palies et al. (Reference Palies, Schuller, Durox and Candel2011c
), a model for premixed swirling flames was determined and it was shown that those two processes drive the flame dynamics.
It is now important to quantify the effect of each mechanism, respectively vortex shedding and swirl number oscillation, on the unsteady heat release. This is achieved through the use of equations (6.3) and (7.2). Indeed, the measured quantities
$I^{\prime }/\bar{I}$
(taken as
$\dot{Q}^{\prime }/\bar{\dot{Q}}$
) and
$S^{\prime }/\bar{S}$
are known and allow the determination of the unknowns, respectively
$A^{\prime }/\bar{A}$
and
$S_{T}^{\prime }/\bar{S}_{T}$
. The computed results are given in figure 16. The total unsteady heat release determined from the integrated heat release of the inverse Abel transform images is plotted as a function of each phase instants for both flames. The contribution UHR1 due to the turbulent burning velocity oscillation induced by the swirl number oscillation and the contribution UHR2 for the flame surface area modification due to the vortex roll-up are given. One can see that the two contributions are out of phase and of similar amplitudes. This result is in agreement with equations (18) and (22) of Palies et al. (Reference Palies, Schuller, Durox and Candel2011c
).
In figure 17, a schematic diagram is proposed to represent the effect of the two identified mechanisms on the unsteady heat release of the two LSI flames. The acoustic modulation induces a mode conversion process at the swirler that generate fluctuations of axial and azimuthal velocities. The swirl number is then modulated by this unsteady velocity field and give rise to an unsteady heat release through (6.3). The axial velocity oscillation leads to an unsteady heat release rate through the vortex shedding. Those two contributions determine the total unsteady heat release. The flame front location (captured with the non-oscillating DMD mode) is a consequence of the type of fuel (methane or hydrogen) and the equivalence ratio.
8 Conclusion
The first major result of this study is the determination of two key contributors for the transition to instability: the flame location and the vortical wave. Firstly, the flame location driven by the equivalence ratio impacts strongly the flame stability. The analysis of three different regimes: stable, local transient to instability and limit cycle respectively, demonstrates that the transition to instability is linked to the flame front location. The flame front location is captured with the non-oscillating DMD mode of the chemiluminescence signal. It is shown analytically that the time average has a non-oscillating component and a component encompassing all others modes. These components are the explicit contribution of a static component and a dynamic component This could be useful in fluid dynamics where the assumption is made that the flow field can be decomposed into a time average and a fluctuating component. Secondly, the vortical wave is responsible for the amplification during the transition. In a network modelling perspective, those two physics are respectively linked to a phase delay and a gain of a flame transfer function. Another key result is that low swirl injector and swirl-stabilized flames anchored with an inner recirculation zone exhibiting self-sustained oscillations are driven by the same mechanisms. In addition, the limit cycle DMD results highlight the nonlinear behaviour of the flame which consists of several peaks at different oscillating frequencies and different growth rates. This regime is studied at the main oscillating frequency to extract the contributing mechanisms. A scenario is developed for those mechanisms in which the swirl number oscillations induce turbulent burning velocity oscillation and the vortical wave rolls up the flame by modifying the flame surface. The unsteady heat release due to both mechanisms is determined and those contributions are in phase opposition with a similar amplitude confirming previous analytical work on premixed swirling flames. Those two mechanisms are respectively affecting the flame intensity through the turbulent flame speed and the flame surface through the vortex roll-up to generate the total unsteady heat release.
Acknowledgements
The authors are thankful to the EM2C Laboratory to allow the use of the Abel transform code, see references Herding et al. (Reference Herding, Snyder, Rolon and Candel1998) and Durox et al. (Reference Durox, Schuller and Candel2005). The authors are thankful to Dr Davis for initiating the cooperation among the authors. The work at LBNL was supported by the Assistant Secretary for Fossil Energy, Advanced Turbines Program, of the U.S. Department of Energy under Contract no. DE-AC02-05CH11231.
Appendix A. Azimuthal and radial waves amplitudes and wave vectors
In the present appendix, expressions for the azimuthal and radial fluctuating amplitudes and their wave vectors are determined. The starting point are the linearized Euler momentum equations in cylindrical coordinates along the radial and azimuthal directions with the axisymmetric assumption for the flow. These equations are the governing equations that can describe the non-reacting perturbed flow downstream of a swirler in linear regime. The fluctuating quantities are assumed to depend only on the longitudinal direction
$z$
. These equations are:
$$\begin{eqnarray}\displaystyle & \displaystyle \bar{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}u_{r}^{\prime }}{\unicode[STIX]{x2202}t}-\frac{\bar{\unicode[STIX]{x1D70C}}\bar{u}_{\unicode[STIX]{x1D703}}^{2}}{r}-\frac{\unicode[STIX]{x1D70C}^{\prime }\bar{u}_{\unicode[STIX]{x1D703}}^{2}}{r}-\frac{2\bar{\unicode[STIX]{x1D70C}}\bar{u}_{\unicode[STIX]{x1D703}}u_{\unicode[STIX]{x1D703}}^{\prime }}{r}+\bar{\unicode[STIX]{x1D70C}}\bar{u}_{z}\frac{\unicode[STIX]{x2202}u_{r}^{\prime }}{\unicode[STIX]{x2202}z}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \bar{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D703}}^{\prime }}{\unicode[STIX]{x2202}t}+\frac{\bar{\unicode[STIX]{x1D70C}}\bar{u}_{r}\bar{u}_{\unicode[STIX]{x1D703}}}{r}+\frac{\unicode[STIX]{x1D70C}^{\prime }\bar{u}_{r}\bar{u}_{\unicode[STIX]{x1D703}}}{r}+\frac{\bar{\unicode[STIX]{x1D70C}}\bar{u}_{r}u_{\unicode[STIX]{x1D703}}^{\prime }}{r}+\frac{\bar{\unicode[STIX]{x1D70C}}u_{r}^{\prime }\bar{u}_{\unicode[STIX]{x1D703}}}{r}+\bar{\unicode[STIX]{x1D70C}}\bar{u}_{z}\frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D703}}^{\prime }}{\unicode[STIX]{x2202}z}=0. & \displaystyle\end{eqnarray}$$
The fluctuating density has an upstream and downstream propagating component:
$$\begin{eqnarray}\unicode[STIX]{x1D70C}^{\prime }=(\hat{\unicode[STIX]{x1D70C}}^{+}\exp (\text{i}k_{\unicode[STIX]{x1D70C}}^{+}z)+\hat{\unicode[STIX]{x1D70C}}^{-}\exp (-\text{i}k_{\unicode[STIX]{x1D70C}}^{-}z))\exp (-\text{i}2\unicode[STIX]{x03C0}ft).\end{eqnarray}$$
To reflect the convective nature (downstream propagation) of the fluctuating velocities
$u_{r}^{\prime }$
and
$u_{\unicode[STIX]{x1D703}}^{\prime }$
along the axial direction, the followings canonical forms are used:
$$\begin{eqnarray}\displaystyle & \displaystyle u_{r}^{\prime }=\hat{u} _{r}^{\prime }\exp (\text{i}k_{u_{r}}z)\exp (-\text{i}2\unicode[STIX]{x03C0}ft), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u_{\unicode[STIX]{x1D703}}^{\prime }=\hat{u} _{\unicode[STIX]{x1D703}}^{\prime }\exp (\text{i}k_{u_{\unicode[STIX]{x1D703}}}z)\exp (-\text{i}2\unicode[STIX]{x03C0}ft). & \displaystyle\end{eqnarray}$$
Injecting (A 3)–(A 5) into (A 1) and (A 2) leads to two equations which are functions of the fluctuating quantities
$\hat{u} _{r}^{\prime }$
and
$\hat{u} _{\unicode[STIX]{x1D703}}^{\prime }$
and their respective wave vector components
$k_{u_{r}}$
and
$k_{u_{\unicode[STIX]{x1D703}}}$
. Coupling the two obtained equations leads to:
$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \hat{u} _{r}^{\prime }\left[-\bar{\unicode[STIX]{x1D70C}}\text{i}2\unicode[STIX]{x03C0}f+\bar{\unicode[STIX]{x1D70C}}\bar{u}_{z}\text{i}k_{u_{r}}+\frac{\bar{\unicode[STIX]{x1D70C}}\bar{u}_{\unicode[STIX]{x1D703}}^{2}}{r\bar{u}_{r}}\right]\exp (\text{i}k_{u_{r}}z)\nonumber\\ \displaystyle & & \displaystyle \displaystyle \quad +\hat{u} _{\unicode[STIX]{x1D703}}^{\prime }\left[\frac{-2\bar{\unicode[STIX]{x1D70C}}\bar{u}_{\unicode[STIX]{x1D703}}}{r}-\frac{\bar{\unicode[STIX]{x1D70C}}\text{i}2\unicode[STIX]{x03C0}f\bar{u}_{\unicode[STIX]{x1D703}}}{\bar{u}_{r}}+\frac{\bar{\unicode[STIX]{x1D70C}}\bar{u}_{\unicode[STIX]{x1D703}}}{r}+\text{i}k_{u_{\unicode[STIX]{x1D703}}}\frac{\bar{\unicode[STIX]{x1D70C}}\bar{u}_{z}\bar{u}_{\unicode[STIX]{x1D703}}}{\bar{u}_{r}}\right]\exp (\text{i}k_{u_{\unicode[STIX]{x1D703}}}z)=0.\end{eqnarray}$$
The next step consists in expressing the real and imaginary parts of this last equation (A 6) and expressed them for
$z=0$
, which corresponds to the swirler outlet. It leads for the real part to:
$$\begin{eqnarray}\hat{u} _{r}^{\prime }\left[\frac{\bar{\unicode[STIX]{x1D70C}}\bar{u}_{\unicode[STIX]{x1D703}}^{2}}{r\bar{u}_{r}}\right]+\hat{u} _{\unicode[STIX]{x1D703}}^{\prime }\left[\frac{-2\bar{\unicode[STIX]{x1D70C}}\bar{u}_{\unicode[STIX]{x1D703}}}{r}+\frac{\bar{\unicode[STIX]{x1D70C}}\bar{u}_{\unicode[STIX]{x1D703}}}{r}\right]=0.\end{eqnarray}$$
And for the imaginary part it leads to:
$$\begin{eqnarray}\hat{u} _{r}^{\prime }\left[-\bar{\unicode[STIX]{x1D70C}}2\unicode[STIX]{x03C0}f+\bar{\unicode[STIX]{x1D70C}}\bar{u}_{z}k_{u_{r}}\right]+\hat{u} _{\unicode[STIX]{x1D703}}^{\prime }\left[\frac{-\bar{\unicode[STIX]{x1D70C}}2\unicode[STIX]{x03C0}f\bar{u}_{\unicode[STIX]{x1D703}}}{\bar{u}_{r}}+\frac{k_{u_{\unicode[STIX]{x1D703}}}\bar{\unicode[STIX]{x1D70C}}\bar{u}_{z}\bar{u}_{\unicode[STIX]{x1D703}}}{\bar{u}_{r}}\right]=0.\end{eqnarray}$$
As a consequence, the ratio of radial and azimuthal convective wave amplitudes is obtained:
$$\begin{eqnarray}\frac{\hat{u} _{r}^{\prime }}{\hat{u} _{\unicode[STIX]{x1D703}}^{\prime }}=\frac{\bar{u}_{r}}{\bar{u}_{\unicode[STIX]{x1D703}}}.\end{eqnarray}$$
Making use of
$k_{u_{r}}=2\unicode[STIX]{x03C0}f/v_{r}$
and noting
$v_{r}$
the convective velocity of the radial fluctuations, the azimuthal component of the wave vector is:
$$\begin{eqnarray}k_{u_{\unicode[STIX]{x1D703}}}=\frac{2\unicode[STIX]{x03C0}f}{\bar{u}_{z}}\left[1+\left(\frac{\bar{u}_{r}}{\bar{u}_{\unicode[STIX]{x1D703}}}\right)^{2}\left(1-\frac{\bar{u}_{z}}{v_{r}}\right)\right].\end{eqnarray}$$
It has been shown in previous work (Palies et al.
Reference Palies, Durox, Schuller and Candel2011a
) that the azimuthal fluctuations propagate at the convective axial velocity of the flow implying that the term in brackets of (A 10) is unity. It implies that
$v_{r}=\bar{u}_{z}$
. This demonstrates that the radial fluctuations also propagate at the axial velocity of the flow. In addition, it has been shown that the phase between the axial and azimuthal components is zero (Palies et al.
Reference Palies, Durox, Schuller and Candel2011a
). As a consequence, as the waves propagate at the same speed, it is appropriate to indicate that the radial and azimuthal fluctuation have the same phase as their initial phasing at the swirler trailing edge is zero.
