Abstract
Lean premixed (LPM) combustion processes are of increased interest to the gas turbine industry due to their reduction in harmful emissions. These processes are susceptible to thermoacoustic instabilities, which are produced when energy added by an in-phase relationship between unsteady heat release and acoustic pressure is greater than energy dissipated by loss mechanisms. To better study these instabilities, quantitative experimental resolution of heat release is necessary, but it presents a significant challenge. Most combustion systems are partially premixed and therefore will have spatially varying equivalence ratios, resulting in spatially variant heat release rates. For laminar premixed flames, optical diagnostics, such as OH chemiluminescence, are proportionally related to heat release. This is not true for turbulent and partially premixed flames, which are common in commercial combustors. Turbulent eddies effect the strain on flame sheets which alter light emission, such that there is no longer a proportional relationship. In this study, phased, averaged, and spatially varying heat release measurements are performed during a self-excited thermoacoustic instability without and with porous inert media (PIM). Previous studies have shown that PIM can passively mitigate thermoacoustic instabilities, and to the best of the authors’ knowledge, this is the first-time that heat release rates have been quantified for investigating the mechanisms responsible for mitigating instabilities using PIM. Heat release is determined from high-speed PIV and Abel inverted chemiluminescence emission. OH chemiluminescence is used with a correction factor, computed from a chemical kinetics solver, to calculate heat release. The results and discussion show that along with significant acoustic damping, PIM eliminates the direct path in which heat release regions can be influenced by incoming perturbations, through disruption of the higher energy containing flow structures and improved mixing.
1. Introduction
There is an increased effort to utilize lean premixed (LPM) combustion processes in gas turbine combustors to reduce NOx emissions and fuel consumption. Unfortunately, these systems are also more susceptible to thermoacoustic instabilities [1]. These instabilities have been shown to damage hardware, produce unwanted noise, and limit the operational range of a combustor [1]. A system is defined to be experiencing thermoacoustic instability when the energy addition, due to an in-phase relationship between unsteady heat release and acoustic pressure, is greater than the energy dissipated by the natural damping of the system [2]. This phenomenon is defined as the Rayleigh criterion [3]. Gas turbine designers continue to look for insightful tools to predict instabilities and assess the effectiveness of mitigation strategies. To predict these instabilities, the models must capture multiple physical interactions including, but not limited to, turbulent flame vortex interaction, system damping, and acoustic wave propagation [4].
To combat these instabilities, different active and passive control methodologies have been implemented. Effective mitigation using active methods is defined by the ability to break the coupling between the unsteady heat release and acoustic pressure fluctuations [5], which will eliminate or reduce a thermoacoustic instability. Active mitigation techniques monitor system response and adjust air and/or fuel flowrates to prevent instability formation [6]. This adds additional complexity and cost to the system but provides a wider range of operating conditions at which the mitigation is possible. Passive mitigation methods, in general, are applicable to a narrower frequency range than active methods but do not require the additional control complexity. Richards et al. [7] present an in-depth and comprehensive discussion of passive mitigation techniques for thermoacoustic instabilities. For example, acoustic resonators can be added to target specific frequencies and extract acoustic energy. Common resonator techniques involve Helmholtz resonators or quarter wave tubes to damp a specific frequency. This was demonstrated by Gysling et al. [8] and Bourquard and Noiray [9] which studied the effectiveness of Helmholtz resonators for a thermoacoustically unstable combustion system. The resonators were shown to be effective at mitigating the instability but required tuning for the specific instability frequency. It has been shown that instabilities arise near a system resonant frequency [10], so modification of the physical geometry can also prevent instability formation. Lahiri and Bake [11] have studied the addition of acoustic liners along the combustion flow path as a mitigation technique. It was shown that the instabilities could be eliminated if the liner is sized properly.
Another area of research in passive mitigation is the addition of porous inert media (PIM). The usage of PIM for passive mitigation was initially proposed by Agrawal and Vijaykant [12]. Their work with PIM has been expanded on in works by Meadows et al. [13–16] and Kornegay et al. [17]. From these previous experimental investigations, it was concluded that PIM mitigates instabilities through its natural acoustic damping properties and the disruption to the recirculation zones in dump combustors. An optimization study was conducted by Sequera and Agrawal [18], where parameters of the porous structures were optimized for instability suppression. This work showed that ring-shaped PIM was effective while also minimizing the pressure drop. In addition, the study determined that it is advantageous to prevent the flame from stabilizing within the porous structure by reducing pore size below the quenching distance. The work by Johnson and Agrawal [19] compared a similarly sized solid ring insert to PIM in a combustion system to identify the independent effects of geometry change versus acoustic damping. The results show that the PIM consistently reduced the instability, whereas the solid insert was less effective at damping or could even amplify instabilities.
To study PIMs effect on thermoacoustic instabilities and to understand the flame response, heat release quantification is critical. Heat release measurements are made more challenging due to the need to study relevant flames. For example, operating under partially premixed combustion [20] instead of perfectly premixed combustion will result in spatially distributed equivalence ratios/heat release rates within the combustor. For laminar, perfectly premixed flames, OH, or CH radical emission are proportional to heat release [21]. Turbulent fluctuations and inhomogeneous fuel/oxidizer mixing in flames eliminate the applicability of direct measurement techniques relating optical chemiluminescence to heat release [21]. Turbulent flames add strain effects due to increased velocity gradients and create a nonlinear relationship between heat release and radical formation [22]. Lauer et al. [22] presented a method to quantify heat release from a lean premixed flame in an atmospheric burner. The method presented uses data from experimentally resolved equivalence ratio, OH intensity, and strain rate. For this method, numerical simulations of a counter flow flame using chemical kinetic solver are also needed to relate OH intensity to heat release as a function of strain rate.
The light emission of specific radicals is often used in combustion studies and these can be resolved using either chemiluminescence with a high-speed camera or planar laser–induced fluoresce (PLIF). Allen et al. [23] demonstrated PLIF measurements in a combustor, by capturing simultaneous acoustic and OH PLIF measurements during instability. PLIF offers the ability to measure spatial radical formation at a single plane, which can then be processed to draw conclusions regarding equivalence ratio and heat release. PLIF requires costly hardware and different laser wavelengths for different species under investigation, whereas chemiluminescence is less costly and complex, only requiring an optical filter change to capture new species. Chemiluminescence captures an integrated line of sight measurement of an excited radical species, which makes it more challenging to determine the spatial distribution of excited radicals. Spatially resolved radical formation is possible by deconvoluting an integrated line of sight measurements via analytical inversion or capture of multiple viewing angles for tomography [24]. With the assumption that the flow is symmetric around a central axis, an Abel inversion of an integrated line of sight intensities can spatially resolve radical formation, as shown by Dribinski et al. [25].
Particle Image Velocimetry (PIV) is another optical diagnostic technique for measuring velocity vectors, strain rates, and turbulence quantities. PIV utilizes a laser sheet and a single camera for 2D velocity components or 2 cameras with stereoscopic PIV to capture 3 components of velocity. Stereoscopic PIV has been used by Arndt et al. [26] to study thermoacoustic instabilities in a combustion chamber, while simultaneously measuring chemiluminescence and laser Ramen scattering. Meadows and Agrawal [14] also studied a thermoacoustically unstable system with PIV. The study investigated the effects of PIM on the flow field and how the change in flow field contributed to mitigation of the instability for a lean premixed system.
Decomposition techniques of resolved flow fields allow for the study of the underlying flow structures. Common decomposition methods are proper orthogonal decomposition, dynamic mode decomposition, and spectral proper orthogonal decomposition. The work by Towne et al. [27] provides a detailed overview of the different decomposition methods and the applicability to different data sets. The result of the decomposition is a series of energy containing modes or structures and the summation of these different modes results in the original flow field. Spectral proper orthogonal decomposition (SPOD) is used on data sets where there are temporally and spatially coherent structures evolving over time. SPOD is different from the other methods due to the resolution of frequencies at which the flow field varies over time and uses those to construct the spatial energy structures. The implementation of SPOD for practical use was presented by Schmidt and Colonius [28].
Previous studies of PIM have sought to elucidate the different mechanisms responsible for PIM’s effectiveness at mitigating instabilities in different combustion environments. From these works, the driving mechanism for instability mitigation was surmised to be natural acoustic damping and the alteration of flow field. These flow field effects likely disrupt the in-phase relationship between acoustic pressure and heat release fluctuations; however, quantification on the impact of heat release rates has not been performed. Past studies with ring-shaped PIM have focused on near perfectly premixed combustion instead of technically premixed systems. In contrast to perfectly premixed combustion, technically premixed combustion will have significant spatial variations in equivalence ratio which result in spatial heat release variations. Dowd and Meadows[29] investigated equivalence ratio fluctuations in a self-excited thermoacoustic instability with and without PIM. The work identified a coupling of the equivalence ratio oscillations with the instability frequency without PIM, and significant reduction in coupling was observed with PIM. In addition, the PIM redistributed the spatial variations in equivalence ratios and provided acoustic damping leading to a stable system. To the authors’ knowledge, the present study is the first to simultaneously investigate the role of equivalence ratio fluctuations, heat release fluctuations, and the velocity field in a thermoacoustically unstable system with PIM. The quantification of heat release provides the necessary insight to understand the instability formation and mechanisms for mitigating instabilities. The current study analyzes an atmospheric, swirl-stabilized combustor with ring-shaped PIM. Simultaneous chemiluminescence and PIV are used to quantify the spatially and phased averaged equivalence ratios and heat release rates. The quantification of heat release follows the methodology proposed by Lauer et al. [22]. The discussion of the paper is as follows: first the experimental setup and methodology are described, detailing the diagnostic tools and processing. Next, the strain rate model and counterflow flame simulations are presented. The results for the acoustic, velocity, and heat release fields are then presented. Finally, conclusions from the work are presented.
2. Methods
To describe the experimental methodology for the current investigation, first the experimental rig is introduced in Section 2.1. This includes details about the optical diagnostics used to quantify chemiluminescence intensities and velocity fields with PIV. Following the method outlined by Lauer et al. [22] to resolve heat release, equivalence ratio and strain rate need to be determined. Section 2.2 details the methodology used to phase average and process line of sight chemiluminescence intensities into spatially resolved, phase-averaged, equivalence ratio values using Abel deconvolution. Section 2.3 provides an overview of the pressure reconstruction methodology and Rayleigh index calculation. The Rayleigh index provides a spatial representation of the areas contributing or damping the instability. Turbulence quantities necessary to determine strain rate and macroscopic turbulence properties are defined in Section 2.4.
2.1. Experimental Setup
An atmospheric, swirl stabilized, dump combustor is utilized to study lean, partially premixed flames, using gaseous methane as a fuel. The dimensional breakdown and image of the physical system are given in Figures 1(a) and 1(b), respectively. The system consists of an air plenum, PIV seed supply, annular mixing tube, optically accessible combustor chamber, and steel combustor section. Porous media, when used in testing, rests in the combustion chamber and is wrapped in a graphite sheet. The graphite forms a tight fit between the quartz combustion chamber and prevents any leakage or burning between the PIM and quartz walls. The PIM is fabricated out of silicon carbide foam with 45 pores per inch and a porosity of 85%. The pore size was selected to prevent burning within the porous structure, as shown by measurements in Sequera and Agrawal [18]. The PIM is ring shaped with a 65 mm OD with a 22 mm hole in the center and is 50 mm long, as depicted in Figure 2. The location of PIM on the baseplate does not interfere with the measurements by the dynamic pressure transducer, PT1; however, the PIM will naturally reduce pressure fluctuations in the combustor due to the acoustic absorbing property of the material.
(a)
(b)
(a)
(b)
The annular mixing tube has an inner diameter of 15.8 mm, outer diameter of 19.0 mm, and is 510 mm long. Fuel is delivered 124 mm upstream of the baseplate of the combustor through the fuel lance which has a 9.52 mm outer diameter. A 45° flat vane, axial swirler is used to improve the mixing of reactants and produces a swirl number of 0.6. The swirler is flush mounted into the baseplate at the end of the annular mixing tube and fuel lance. The swirler has the same diameter as the mixing tube to maximize flow area and prevent any step into the flow path. The combustion chamber consists of a quartz cylinder with an inner diameter of 70.0 mm and a length of 410 mm, which sets concentric to a steel combustor section that has an inner diameter of 64.0 mm and is 102 mm long, used to hold the quartz cylinder and instrumentation.
The main portion of air is injected through a chocked orifice in the plenum and fuel is injected through 6 circumferentially spaced openings on the fuel lance that runs concentric to the mixing tube. Mass flow controllers (Alicat MCR) are used to measure and control the air and methane flowrates into the system. Air flow up to 10.34 g/s (500 SLPM) and fuel flow up to 0.645 g/s (50 SLPM) are possible with these two devices, with an uncertainty of ±0.8% of reading which correlates to a max uncertainty of 82.6 mg/s and 4.37 mg/s, respectively. A secondary air supply is brought in through the PIV seed supply, located between the plenum and mixing tube. This secondary air flow carries aluminum oxide (Al2O3) seed particulate used as tracers to visualize velocity vectors. The particulates have a mean diameter of 0.3 μm, resulting in Stokes number of 1e − 4. The seed injection is approximately 32 L/D upstream of the swirler. At nominal conditions, flow with a Reynolds number of approximately 21,000 is in the mixing tube, resulting in an approximate distance of 23 L/D for the flow to become fully developed. The seed flow will have sufficient time to mix after injection due to the upstream injection location and homogenous seeding is assumed. As the aluminum oxide particulate would quickly coat the quartz tube and distort optical diagnostics in the combustor, an alternate flowrate of air was supplied through the PIV seed supply that contained no aluminum oxide particulate. This alternate nonseed flow was used to ignite and get to the desired operating point with an instability. Since the flowrates between the seed and nonseed flow need to match to not alter the instability, a needle valve was installed to control the nonseeded flow. The flowrate of the secondary seed flow was measured with a 500 SLPM mass flow controller (Alicat MCR), with the same range and uncertainty as the controller used on the main air supply to the plenum. The flow controller was only used as a passive flowrate measurement device. A schematic detailing the seed flow paths is depicted in Figure 3.
Inlet air temperature is determined by a K-type thermocouple in the air plenum. Pressure perturbation measurements are obtained via microphones (PCB-Y378C10, ¼” ICP Microphone System) in the mixing tube and dynamic pressure transducers (Kistler 6025A) in the baseplate and steel combustor. The microphones are located 62.43 mm from each other and start 247.2 mm upstream of the baseplate, as detailed in Figure 1(a). The microphones and dynamic pressure transducers have a sensitivity of 2 mV/Pa and 103 pC/bar, respectively. The dynamic transducers can operate up to 700°C, and the signals are amplified with a charge amplifier with a 400 mV/pC gain. For this experiment, mass flow rate, temperature, and pressure transducers/microphones are sampled at 1 Hz, 1 Hz, and 100 kHz, respectively.
A high-speed camera (Photron FASTCAM SA5) and high-speed intensifier (SIL-25HG50D) are used to capture chemiluminescence emission intensity. A Nikon Rayfact UV lens (PF10545MF-UV) with 105 mm focal length was used and the imaging was focused near the baseplate through the quartz cylinder. Flame intensity is filtered with monochromatic filters, 310 nm with a FWHM of 10 nm and 430 nm with a FWHM of 10 nm for OH and CH, respectively. The images are captured at 5 kHz with an 1024 × 1024 pixel resolution. The spatial resolution of the images captured was 0.2 mm/pixel. A gain of 9, exposure time of 180 μs, and delay of 55 ns are used for the intensifier to properly amplify the images.
Velocity measurements of the flow field were obtained using time-resolved PIV. A Photonics dual head DM150-532DH Nd:YAG laser with a wavelength of 532 nm, pulse energy of 30 mJ, and 10 kHz repetition rate is used. The time between two laser pulses was 8 s. A LaVision VZ-Beam steering arm with a variable laser sheet optic with = 50 mm lens configured to generate a 1.5 mm thick laser sheet in the field of investigation was used. A Phantom v2512 high-speed camera with a Sigma 105 mm focal length lens and 1280 × 800 pixel CMOS sensor was used to capture data at 10 kHz, resulting in velocity measurements at 5 kHz. A 532 nm bandpass optical filter is used to filter out unwanted light. For this experiment, the spatial resolution was 52 μm/pixel, resulting in a minimum velocity resolution of 6.5 m/s.
Postprocessing of the PIV images was conducted in DaVis 10.0 from LaVision. The velocity data were preconditioned by using a time filter and geometric mask to better highlight the particulate in the flow. The calculation takes place using a multi-pass vector method, initially involving 2 passes with a window size of 48 × 48 pixels and a 50% overlap and then 4 final passes with a window of 12 × 12 pixels and 87.5% overlap. The computed vectors are then validated using a universal outlier detection method, where a local filter region of 5 × 5 pixels is used. Results are eliminated if residuals fall outside of the 2.5–3.5 range. Final validation is conducted, deleting any computed vectors that do not have a correlation value above 0.5. The resulting spatial spacing is 1.29 mm between vectors. Using this spacing with the minimum velocity resolution and properties of air at 1600 C, results in a minimum turbulent Reynolds number of 25. Using the Kolmogorov Reynolds number of 1 as reference, the velocity vectors calculated here are an order of magnitude larger than the Kolmogorov eddies in the system.
With the use of multiple diagnostics, proper alignment of laser and camera setups is critical so that there is no interference between measurements. To accomplish this, the VZ-Beam steering arm is positioned above the chemiluminescence camera, as shown in Figure 4. Precise timing was critical during testing, as large data collection and coating of PIV seed material on the quartz window limited data collection to a few seconds. Timing was specified in the system control logic to trigger instrumentation, control seed flow, and extinguish the flame after the test interval. System control and data collection were done in LabVIEW. The camera, PIV system, microphones, and pressure transducers are configured for synchronized data collection.
2.2. Phase-Averaged Abel Inverted Equivalence Ratios
In the presence of a strong thermoacoustic instability, there is a dominant frequency at which the system is fluctuating. This frequency can be used to phase average results at discrete phases along the period of the instability. By phase averaging, simultaneous imaging is not required to capture both OH and CH intensities. To obtain spatially defined intensities at the center plane, the Abel inversion method presented by Dribinski et al. [25] is used to reconstruct the line of sight OH and CH images. Equation (1), from Dribinski et al. [25], represents the governing equation of the inverse Abel transform. This equation relates a two-dimensional projected image on the detector plane, coordinates x and z, to the original three-dimensional light emission, coordinates r and z.where I is the intensity of the original 3D source, r is the radius, and P is the intensity on the detector plane. The reconstructions presented in the current work use 500 Gaussian basis set functions to reconstruct each image. The core assumption of the reconstruction is that the flow and fuel distribution is axisymmetric. For the current investigation, the flow can be assumed to be axisymmetric except for minor variations due to discrete fuel injection locations causing slight circumferential equivalence ratio asymmetry. The variations are considered negligible in comparison to the expected longitudinal fluctuation due to the instability.
The study by Hardapulas and Orain [30] showed a monotonic dependence of the equivalence ratio on the ratio of OH to CH radical light emission. Hardapulas and Orain [30] studied counter flow gaseous flames to understand the influence of strain rate and equivalence ratio on the production of OH, CH, and C2 radicals. Correlations relating to the ratio of OH/CH depend on the transmissivities of the various optical components. In the present study, a correlation developed by Hugger [31], which uses the exact same optical setup and settings in the present study, is used to relate the ratios of OH/CH to equivalence ratios as shown in the following equation:where is the equivalence ratio. Equation (2) from Hugger [31] is based on measurements of equivalence ratios ranging from 0.7 to 1.3 and extrapolation beyond this region will likely introduce error. For this reason, only ratios of OH to CH will be used that produce equivalence ratios from 0.7 to 1.3. Also, due to the way equation (2) is formulated, there are points at which the resulting equivalence ratio would be undefined. For this reason, some filtering of the data is required to avoid unrealistic values. Hardapulas and Orain [30] determined that the ratios of OH;/CH are independent of strain rate, thus, it is expected that the correlation developed by Hugger [31] is also independent of strain rate.
2.3. Rayleigh Index and Pressure Reconstruction
To better understand how the spatial distribution of heat release couples with the acoustic pressure in a combustor to form a thermoacoustic instability, the Rayleigh index is used [32, 33]. The Rayleigh index is defined in equation (3). The form of the equation is intended to study a set of heat release and pressure measurements that occur over the period of an instability for a spatially resolved cross section. The Rayleigh index quantifies the energy added from the in-phase relationship between pressure and heat release, where a positive value states that energy is being added to the system and negative values signify damping. The total number of instances, , have the same temporal spacing between measurements.where is the ratio of specific heats, is the period of the instability, is pressure, and is the heat release, with the overbar denoting mean pressure and apostrophe denoting fluctuating components. The acoustic pressure is used for in the above expression and is recorded via wall-mounted pressure transducers. To resolve the spatially distributed acoustic pressure, Hale et al. [34] presented an analysis technique to reconstruct the pressure mode shape using surface-mounted pressure transducers at different axial locations in the combustor. The methodology utilizes a reduced form of general acoustic pressure governing equation and a least square curve fit function to predict the variation in acoustic pressure between measurement locations. The reconstruction assumes that the mode is purely longitudinal with open acoustic inlet and exit boundary conditions that represent the system being analyzed.
2.4. Macroscopic Turbulence Properties
Turbulent velocity fluctuations, the rms values, can be calculated from the instantaneous velocity vectors for N instances [35] as shown in the following equation:where i represents the different components of velocity and represents an average over time. The integral length scale, , is defined as the integral of the two-point correlation function, , normalized by [36], as presented in the following equation:where and are the radius and unit vector for of each component of the velocity, respectively. Equation (5) presents a challenge as the velocity measurements have similar spatial resolution as the integral length scales, such that direct evaluation of equation (5) will result in large uncertainties. Instead, an analytical two-point correlation can be used as shown in [37]
The integral length scale can be determined by fitting experimental data with the exponential function, with as the only unknown. For additional details regarding the process, please refer to Lauer et al. [22].
3. Strain Rate Model
Due to the spatial resolution of the PIV employed in this study being larger than the Kolmogorov length scale and no out of plane velocity measurements, a method for modeling the strain rate given the known quantities is needed. The method presented here follows a discussion by Lauer et al. [22]; the method calculates a probability density function (PDF) to define the strain rate at each spatial location. The PDF for strain rates of different types of turbulence is given by Yeung et al. [38]. For the methodology to be applicable, it must be shown that the flame can be approximated as having constant property, homogeneity, and isotropic turbulence. Also, it must be determined whether the flame can be described as a material or randomly oriented surface. It has been concluded by Bradley et al. [39] and Pope [40] that the main contribution to flame straining of a premixed flame is due to the turbulent structures upstream of the burning. Due to this finding, the statistical strain rate distributions can be calculated from the turbulent properties of the unburnt mixture upstream of the flame, which are constant property and homogeneity. Bradley et al. [39] and Chen and Im [41] have compared predictions from the strain rate PDF to experimental and direct numerical simulation results. Due to the good agreement, it can be concluded that although Yeung et al. [38] does define the strain rate PDFs from simulations of nonreacting flow, the distributions reliably define the strain in isotropic turbulent flames.
For flow to be considered isotropic, the time averaged product of the orthogonal components of velocity fluctuations must equal zero, as shown in the following equation:
As expected for real flames, the turbulence defined in this study is not isotropic. The product of the orthogonal time averaged product does not equal zero, as seen in Figure 5. The axial coordinate distance, , and radial distance, , are normalized by the radius of the quartz, . The orthogonal product, , is compared to the axial component fluctuation product in Figure 5. Figure 5 only shows representative results for the first phase of the non-PIM case. All phases and PIM cases show similar results when comparing the fluctuating products. Comparing the two contours, the orthogonal component is significantly less than the axial component, justifying the assumption that the macroscopic turbulence can be approximated as isotropic. Furthermore, since the straining of the flame is not caused by macroscopic turbulence of the flow field, Kolmogorov’s hypothesis [42] of local isotropy is valid and small-scale turbulence can be considered isotropic. Therefore, it is acceptable to approximate the macroscopic turbulence as isotropic. Unfortunately, it is not possible to quantify the true error associated with this assumption given the current data; however, it is tenable that this approach will be able to predict regions of high strain rate and estimate the local impacts of strain rate on heat release rate fluctuations.
The flamelets observed in the current study must now be classified as either a material or randomly oriented surface. A material-oriented surface does not propagate by itself and is strongly influenced by the turbulent structures, whereas a randomly oriented surface moves independent of turbulence [38]. To define a flamelet as a material surface, the Kolmogorov velocity must be greater than the laminar flame speed. The laminar flame speed can be determined using the locally defined equivalence ratio and a reduced order chemical kinetics solver, Cantera [43], calculated from freely propagating flame of a premixed methane-air mixture. The Kolmogorov velocity calculation is reproduced in the following equation:where is the kinematic viscosity of the unburnt mixture, is the Kolmogorov velocity, and is the Taylor length scale. The Taylor length scale is approximated by which utilizes the spatially defined integral length scale. This definition of Taylor length scale has been used in works investigating turbulent length scales both experimentally and numerically [44–46]. The spatial variation in laminar flame speed and Kolmogorov velocity is compared in Figure 6. Since the equivalence ratio is only defined between 0.7 and 1.3, there are more undefined regions in the laminar flame speed contour than the Kolmogorov velocity. Again, only representative data for the first phase of the no PIM case is presented. As detailed in the contours, the laminar flame speed is clearly less than the Kolmogorov velocity. Therefore, the flamelets in the current study are approximated as material surfaces.
The flamelets have been defined as material surfaces and it has been shown that the turbulence can be approximated as isotropic, constant property, and homogenous. This leads to the statistical distribution of strain rate, , reproduced in equation (9), following the strain rate model proposed by Yeung et al. [38]. The turbulent time scale, , is reproduced in equation (10).
4. Counterflow Flame Simulation
To properly define a relationship between optically resolved radical emissions and heat release, numerical simulations of one-dimensional counterflow flames are used. Numerical simulations are used instead of experimental quantification, due to the difficulty and unreliability of experimental heat release. The 1D numerical simulation is conducted using the chemical kinetics solver, Cantera [43]. The following section briefly outlines how the relationship between OH chemiluminescence intensity, equivalence ratio, and volumetric heat release is obtained. The methodology presented here follows Lauer et al. [22], and a high-level overview is presented, for further detail the reader is referred to the paper.
The counterflow flame simulation is setup with opposed inlets creating a stagnation plane. A flame will anchor at or directly adjacent to the stagnation plane. The results of these simulations will produce a table of proportional correction values that can be applied to OH intensity values to obtain heat release. The correction value will be used at each spatial location and will be dependent on the local strain rate and equivalence ratio, as both factors strongly influence the amount of heat release. Therefore, multiple simulations will be conducted to independently vary strain rate and equivalence ratio to calculate the resulting ratio of heat release to OH intensity for each case. A fresh-to-burnt setup is used, meaning a fresh gas inlet brings in CH4, O2, and N2, corresponding to the equivalence ratio, and at the burnt gas inlet, a zero gradient boundary condition is specified for the volume fraction of all species and temperature. Fresh-to-burnt was shown to provide a more realistic description for premixed and partially premixed flames [39]. The inlet mass flow rates for the two boundaries are used to define the strain rate for each simulation. The reaction mechanism from Kathroita et al. [47] is used for the chemical kinetics. Strain rate, , is varied between 500 and 20,000 and equivalence ratios are varied between 0.7 and 1.3. The range for equivalence ratio is selected based on the range of applicability of equation (2) for the calculation of equivalence ratio.
The results to be extracted from each simulation will be the axially varying heat release profile, , and the volumetric OH-intensity, . The volumetric intensity of OH is calculated using the Einstein coefficient of spontaneous emission and the mole fraction of OH. The strain rate and equivalence ratio dependent proportionality factor, , can now be defined in equation (11), reproduced from Lauer et al. [22].
A few comments about using this expression and the resulting table are as follows: First, the table is discretely defined and to find points in between those, a bicubic interpolation is used. Second, these correction factors are specific to the current experimental setup; if the fuel composition, system pressure, or preheating is changed then a new set of simulations will need to be run.
4.1. Heat Release Calculation
The proportionality factor defined in equation (10) cannot directly relate heat release to OH∗ intensity for experimental data due to the PDF definition of strain rate from Section 3.2. To obtain the nonlinear relationship between OH∗ intensity, heat release, and strain rate, a new strain rate averaged proportionality factor, , must be defined. At each spatial location, the equivalence ratio and strain rate PDF, are used in equation (12), reproduced from Lauer et al. [22].where the quench strain rate, , is defined as the strain rate at which the integral heat release falls below two percent of its unstrained value and is equivalence ratio dependent. This results in a proportional term relating the time averaged OH∗ intensity to the heat release for each spatial location under investigation. As the light intensity measured for OH∗ chemiluminescence does not quantitatively define the photon count, the product of OH∗ chemiluminescence and results in relative heat release measurements. The above product can then be normalized by its mean value and used with the integral heat release rate, determined from an energy balance on the system, to find the actual heat release in system.
5. Results and Discussion
The following results section is divided into three subsections for discussion of the independent acoustic, velocity, and heat release results. A self-excited instability was found at an air flow rate of 4.48 g/s (225 SLPM) and global equivalence ratio of 0.85. The flow split was 55% of the total air flow being unseeded and the remainder from the secondary seed flow path.
5.1. Acoustic Results
The resulting self-excited instability without PIM resulted in a frequency of oscillation at 312 Hz and sound pressure level (SPL) of 157 dB. When PIM was added to the system, the prominent frequency shifted to 307 Hz at 119 dB. The addition of PIM results in a 38 dB reduction in peak frequency SPL. The SPL levels are calculated with reference pressure 20e − 6 Pa. To summarize the acoustic response of the system, the Fast Fourier Transform (FFT) results are presented for both PIM and No PIM in Figure 7. The results shown were collected from the dynamic pressure transducer mounted on the baseplate.
The no PIM spectrum in Figure 7 highlights a large peak at 312 Hz, representative of the instability frequency and its second harmonic at 624 Hz. The excited frequency is assumed to be the first longitudinal mode of the combustion chamber based on the geometric length of the combustor. Fundamental acoustic resonance prediction was conducted on a representative open-open series of ducts to determine resonance frequencies.
To sync and temporarily average the results from different measurement sources, phase averaging is used. The different phases are defined in Figure 8, where the period of instability is divided into 10 phases. The dynamic pressure sensor in the baseplate is used for phase synchronization. Figure 8 is representative of the acoustic pressure variations in the combustor during an instability period.
The different acoustic measurements at the various locations highlighted in Figure 1 are used to reconstruct the pressure mode shape in the combustor and mixing tube during the thermoacoustic instability. The plenum was not included in modeling due to the large area change entering the mixing tube, which responds similarly to an open boundary. Both inlet and exit boundaries are treated as acoustically open. The exit of the combustor is given artificial additional length to account for the effects of venting to ambient as detailed in Ref. [4]. One downfall of the reconstruction technique presented by Hale et al. [34] is its applicability to single ducts. This results in separate solutions for the mixing tube and combustor. The pressures are all normalized by the maximum pressure that occurs in the system, approximately 4 kPa. The reconstruction is done only for the no PIM case. The PIM produces a similar mode shape with a large reduction in magnitude. The resulting mode shape of the pressure reconstruction is presented in Figure 9. The locations of the microphones and dynamic pressure transducers are represented by squares and diamonds on the x-axis, respectively. The fuel injection location is represented by a circle for reference.
(a)
(b)
The mixing tube fluctuations are represented from X/L = 0 to 0.52 in Figure 9(a), and the combustor in Figure 9(b) fluctuations are from X/L = 0.52 to 1. The resulting mode shape shows larger amplitude pressure fluctuations in the mixing tube than in the combustor. The exit of the combustor shows that the mode shape prediction seems to be approaching zero but does not converge to zero as expected. This is due to the additional artificial entrainment length added to the combustor, resulting in the higher acoustic pressure at the physical exit of the combustor. Pressure from the mixing tube to the combustor should be continuous. The results in Figure 9 show a step in the pressure from the mixing tube to the combustor. This step is a result of the attenuation of the dynamic pressure transducer in the baseplate, which is in a recessed cavity to protect the device from overheating. The impact of attenuation due to the cavity will result in a reconstruction magnitude lower than the actual magnitude. The pressure reconstruction will be used in the heat release section for the calculation of Rayleigh Index.
The reconstruction in Figure 9 is similar to the mode shape prediction presented in Ref. [29] for a representative system, and this increases the confidence in the reconstruction. To validate the reconstruction presented in Figure 9, a FFT is performed on the raw pressure time histories for the three microphones and two differential pressure transducers (DPT) for the unstable no PIM case. The results are presented for the frequency of instability in Table 1. The gain, relative to the baseplate DPT, and the phase difference, between the measurement and the baseplate DPT, are reported. Phase difference is reported in degrees. The naming of measurements and their relative location corresponds to those called out in Figure 1.
The results in Table 1 confirm the mode shape reconstructed in Figure 9, with some explanation. The gain can be seen to be higher in the mixing tube than the combustor magnitudes, which supports the results in the reconstructed mode shapes. The gain also decreases when moving upstream in the mixing tube toward the plenum, X/L = 0. One discrepancy is the phase difference between combustor and mixing tube; the reconstruction is predicting an almost in-phase relationship, whereas the FFT results are higher than expected. This can be explained by the presence of the swirler, which cannot be modeled in the reconstruction. The swirler will partially reflect acoustic waves from both the mixing tube and combustor, which will lead to a discrepancy in the phase, and it will accelerate the flow leading to a convective wave as opposed to a purely standing wave. The microphones in the mixing tube are in phase with one another, and the two dynamic transducers in the combustor are also in phase. This phase difference is between the combustor and mixing tube and will not have an effect when phase relationships are made using the combustor measurements in upcoming sections.
The acoustic results have demonstrated the large reduction in acoustic SPL with the addition of PIM, which will significantly reduce acoustic fluctuations throughout the system. Acoustic perturbations generated by the flame or system resonance propagate at the speed of sound and are present at all axial locations in the system. This is critical at the fuel injection location, as the acoustic perturbations generate small flow field perturbation which will affect local mixing and equivalence ratio, which are then convected to the flame front. The work by Dowd and Meadows [29] further explain how equivalence ratio fluctuations are generated at the fuel injection point. Acoustic interactions were shown to alter the equivalence ratio, which are then convected to the baseplate, resulting in an in-phase relationship between equivalence ratio and acoustic pressure. The convected variations in equivalence ratio cause heat release fluctuations and grow the instability until a limit cycle is reached. The in-phase relationship between acoustic and heat release fluctuation has been stated as a primary feedback mechanism for thermoacoustic instabilities [4]. This highlights the importance in the decrease in the acoustic perturbations when adding PIM and a major factor for system stability at a nominally unstable operating point. To better understand what physical processes are contributing to the instability mitigation, the velocity fields are analyzed.
5.2. Velocity Results
Using PIV image capturing, and DaVis post processing software detailed in Section 2.1, the instantaneous velocity fields in the combustor can be computed. A Stokes number much less than one is desired for particles to closely follow streamlines and give an accurate representation of the flow field [48]. For the current investigation, seed with a Stokes number of 1e − 4 was used, resulting in the particulates closely following the flow streamlines. Analysis presented in Melling [49] show that the max frequency able to be resolved with this seeding is greater than 20 kHz. The seeding is sufficient to resolve the flow field investigated here. Using the phases defined in Figure 8, a phase averaged velocity field is calculated and presented in Figure 10 for the no PIM results. The results for the PIM case compared to a single phase of the no PIM velocity magnitude is shown in Figure 11. A single phase from PIM results is presented due to the similarity between all phases. The velocity was captured and reported from the centerline to the quartz wall, to maximize the spatial resolution captured by the PIV camera pixels. The radial and axial spatial coordinates in Figures 10 and 11 are normalized by the radius of the quartz. The term in Figure 11 is defined as the relative axial length divided by the quartz radius, due to the PIM results being 50 mm axially downstream of the No PIM results. A = 0 in the PIM and no PIM cases correspond to the PIM downstream surface and the swirler exit, respectively. The velocity results have regions with no flow data reported, due to reflections from the quartz cylinder interfering with PIV measurements at that location. The region outside = −0.9 is also omitted from the results for this reason.
The results of the phase averaged velocity for the no PIM results show three major flow structures. First, a high velocity jet exiting the swirler, a corner recirculation zone formed between that jet and the baseplate, and a central recirculation zone, traveling back to the baseplate along the centerline. A precessing vortex core (PVC) may be positioned in the central recirculation zone; however, out of plane velocity information is needed to confirm the presence of PVC. There has been significant discussion of the PVC and its effect on thermoacoustic instability formation [14], and both recirculation zones and PVC are feedback mechanisms that support coupling in combustion systems. From the no PIM results, the spatial distribution of velocity is unaffected by the instability but there is a strong effect on the magnitude. In the high velocity jet exiting the swirler, there is a variation in spatially phased averaged velocity fields of up to 20% from the time averaged value.
The PIM flow field presented in Figure 11 is dominated by a jet exiting the core of the PIM ring with a recirculation zone extending over the porous structure. It is interesting to note that the jet exiting the PIM opening does not extend across the entire opening, as seen by the lower flow near the centerline. The flow exiting the swirler encounters a radial PIM blockage, forcing the flow axially, with the higher velocity flow remaining in the near wall region at the edge of the opening. The PIM flow field shows lower velocities in comparison to those in the no PIM case, again due to the impingement, turning of the flow by the porous structure, and being a further distance downstream where mixing and dispersion will reduce the velocity. The effects of PIM on the flow energy breakdown will further be discussed in the SPOD results. Although recirculation zones are still present, their effect on system stability in the PIM case will be less. Instead of interacting directly with the baseplate and jet exiting the swirler, as in the no PIM case, the recirculation zone forms over the porous structure. This is not to say that recirculation zones do not form at the swirler exit in the PIM case. If these recirculation zones are present, they would not have the same interaction with the flame as in the no PIM case due to the lack of burning in the PIM. Sequera and Agrawal [18] have studied thermoacoustic mitigation using PIM and the effects of different sized PIM structures. The study monitored thermal radiation levels in the visible spectrum of the PIM to determine if burning was occurring in the porous structure. The PIM used in this investigation was sized similarly to a porous structure from Ref. [18] that had limited radiation, which is a result of minimal burning within the porous structure. The velocity results in Figure 11 exiting the central region of the PIM are also significantly higher than the turbulent flame speeds [50], making it unlikely for a flame to stabilize in that core flow region within the PIM. Turbulent flame speed between 0.4 and 2 m/s was found by Vargas et al. [50] for equivalence ratios between 0.8 and 1 for atmospheric methane combustion. As the majority of the burning is interacting with the recirculation zones at the exit of PIM, the pressure gradient induced by PIM and the physical distance to the baseplate significantly reduce burning interaction with the flow entering the combustor. Recirculation zones are pathways for fluctuations to be convected upstream and are both supportive of instability formation. Therefore, the elimination and reduction of these effects improve the system stability.
The velocity field is further investigated in Figure 12, where the spatial average velocity of each phase is plotted over the instability period. Figure 12 also includes results for the PIM velocity and the acoustic pressure fluctuations with phase. The results for velocity and pressure are presented as the variation from their respective mean value for each phase. The results are then normalized by dividing the max fluctuation for each set of results. The acoustic pressure is only a single curve for PIM and no PIM, since we are normalizing each curve by its respective max fluctuation. The absolute value of the two pressure traces is significantly different, as shown by the SPL levels in Figure 7. These results more quantitatively exemplify the conclusions drawn from Figures 10 and 11, involving the no PIM fluctuation or PIM lack of fluctuation in mean velocity.
The variation between results shows a phase shift of approximately 100 degrees between the acoustic pressure and mean velocity fluctuations of the no PIM case. This phase shift is physical and expected. The results in Figure 12 are for mean velocity fluctuation and not acoustic velocity, but the same rationale to explain why acoustic velocity and pressure are separated by 90 deg in a standing wave can be applied here. PIV is measuring the summation of turbulent and acoustic velocity fluctuations, but for the no PIM case, the acoustic effects are much more prominent than the random turbulent fluctuations. The velocity response is therefore expected to follow the acoustic relationship between velocity and pressure. At the acoustic pressure maximum, the particle displacement and resulting acoustic velocity will be zero, resulting in a 90 deg shift in phase. As the system is experiencing a thermoacoustic instability, it has been shown that an in-phase relationship between heat release and acoustic perturbations [4] is likely present. The actual phase relationship between acoustics and heat release will be further studied in Section 5.3.
To further study the flow field and the underlying energy containing structures, SPOD is performed on the velocity data. The SPOD methodology and tools provided in Schmidt and Colonius [28] are used here. The data set collected studies the product of the axial () and radial () velocity fluctuations. This term will be indicative of the turbulent kinetic energy in the flow field. The measurements in each data set are separated by 0.002 s, and a window size of 128 was used to divide the total set into different bins. A 50% overlap between data placed in each bin was specified. This window size was optimized by an independent study, where window size was varied until the “spectral leakage” of energy between frequencies was minimized and the further increase in window size does not improve the resolution of energy breakdown. A weighting factor was also applied to the data to account for the radial geometric dependence. The cells near the wall were given a larger weighting since they encompass a larger percentage of the total volume than cells near the centerline. To confirm that the flow field energy distribution was not varying with time and ensure that a sufficient number of time instances are used, various data sets throughout the PIV data collection period were analyzed. The results produced a similar energy breakdown and spatial distribution of energy at the instability frequency and highest mode energy. The final decomposition used 0.006 s worth of data to properly decompose the flow field. A breakdown of the amount of energy contained in each frequency at different modes is presented in Figure 13. The first 5 modes and results for frequencies up to 1000 Hz are shown for both the no PIM and PIM results.
(a)
(b)
The results in Figure 13 show the lack of distinct peaks, observed by comparing the 312 Hz frequency mode energy and significant reduction in kinetic energy in the modes, in the PIM case. It can be concluded that the energy in the system has been redistributed to higher order modes. In contrast, the no PIM results show that there is a significant flow structure producing energy at 312 Hz and its harmonic 625 Hz, which is at the instability frequency demonstrated in Figure 7. Similar results were presented in Meadows and Agrawal [14], where proper orthogonal decomposition (POD) of the flow field showed a similar redistribution of energy when PIM was added to an unstable system; however, POD does not determine which structures couple with the instability frequency. To better understand how PIM is disrupting the energy containing flow structures, the magnitude spatial distribution of energy is presented for PIM and no PIM cases in Figure 14. Only the spatial distribution of the highest energy containing mode is presented at the instability frequency. Again, represents the relative axial displacement since an axial distance of 50 mm separates no PIM and PIM results.
The spatial distribution of the PIM decomposition shows a small energy contributing structure near the centerline at the exit of the PIM ring. The remainder of the PIM flow field is composed of individual vortices being shed from the higher velocity core jet. The no PIM flow field shows a large coherent energy structure along the flow path exiting the swirler. The decomposition also shows that the shear layer along the lower portion of the core jet is contributing energy to the overall system as the baseplate recirculation zone is formed. The PIM addition is disrupting the energy producing flow field structures. The lack of these flow features in the PIM case is confirmed in Figure 11. The presence of these higher energy flow structures in the no PIM case will create larger velocity and mass flow fluctuations entering the combustor. Corner recirculation zones will generate flow oscillations entering the combustor due to the opposing forces generated by these structures. A higher energy flow in the recirculation zone will generate larger flow oscillations than a lower energy structure. The results in Figure 12 show that these flow structures affect the mean velocity with a 20% variation from the mean. The result of these fluctuations will be investigated further in Section 5.3.
This does not mean there are no mean flow feedback mechanisms to the incoming flow in the PIM case. The reduced energy and multiple structures reduce the likeliness that the energy produced by this interaction would be significant. Any coupling with PIM structures will be with a less energetic flow structure, which will have less influence on system response and will lead to a more stable system with the addition of PIM.
5.3. Heat Release Results
To highlight the oscillating flame behavior without PIM and the stable flame behavior with PIM, Figure 15 shows raw OH chemiluminescence images at three different time steps corresponding to 1/3 of the instability frequency. For determining heat release rates, please refer to Section 3 for a detailed description. Figure 16 presents the results of each major processing step. Starting with raw instantaneous OH and CH chemiluminescence imaging, a phase averaged Abel inverted processing takes place which can be used to compute spatial equivalence ratios. Finally, with PIV and the look-up tables generated from Cantera simulations, the actual heat release is calculated. It is important to note that due to the large spatial resolution and 2-D nature of the PIV data, strain rate is approximated assuming isotropy. Of course, the present flow field is not truly isotropic, and is a limitation of the current approach; however, since straining of the flame is dominated by small scale turbulence, and Kolmogorov [42] argued that the directional biases of the large scales are lost in the chaotic scale-reduction process as energy is transferred to successively smaller eddies, small-scale turbulent motions can be approximated as statistically isotropic. In addition, the orthogonal products, << suggest that approximation of isotropy for macroscopic turbulence is reasonable; however, the true error for this assumption could not be quantified given the present data set. The results are all presented for the 36 deg phase of the no PIM case, and results for chemiluminescence and Abel inverted values are normalized by their respective max intensities. Equivalence ratio and heat release, in Watts/m3, are presented as absolute units. The results for the PIM case are computed similarly. Equivalence ratio values are only presented between 0.7 and 1.3 as that is the applicable range of equation (2). Values outside of this range will not be included in the results. A weak Gaussian blending is performed on the equivalence ratio distribution to eliminate any local discontinuities between adjacent pixels. Heat release results are only presented where a valid equivalence ratios were calculated.
(a)
(b)
(c)
(d)
From the processing steps, it is important to highlight the necessity of the Abel inversion and resolved strain rate to calculate heat release. The effects of Abel inversion can clearly be seen moving from Figures 16(a) to 16(b), where incorrect conclusions would have been drawn from the line-of-sight integrated values about the flame shape. For example, in the region near the centerline from = 0.5 to 1, high intensity would have led to conclusions of a large flame presence but after inversion and processing, there is no significant heat release in that area. The strain rate importance will be further elucidated in upcoming figures, but the results here show that, moving from Figures 16(c) to 16(d), the strain rate effects the magnitude of heat release. If studying an unstrained laminar flame, similar equivalence ratios should yield similar heat release rates. Equivalence ratios around 1 are noted at and in Figure 16. This does not produce the same heat release as noted at those locations in the heat release in Figure 15. Higher heat release is noted for similar equivalence ratio further from the baseplate. The difference is caused by the straining of the flame. Using Figure 10, the location closer to the baseplate has a higher velocity magnitude/gradient in the 36 deg phase in comparison to the location further from the baseplate. The higher gradients result in a higher strain rate near the shear layer and less efficient burning, producing lower levels of heat release. It should be noted that direct strain rate results are not presented due to the method in which they are calculated, outlined in Section 3. A strain rate PDF is defined at each spatial location based on turbulence quantities and used to determine a strain rate averaged proportionality factor, which relates OH∗ intensity to heat release.
Figures 17 and 18 highlight the evolution of the heat release fluctuations, equivalence ratio fluctuations, and the SPOD mode at 312 Hz throughout the instability cycle for the no PIM case. The variation in heat release with phase shows that the majority of the fluctuating heat release is occurring near the impingement point of the flow on the quartz wall. As the flow impinges, it will slow down and the strain rate will be reduced, which results in more efficient burning. Furthermore, the SPOD mode is out of phase with the HRR rate. Pressure fluctuations inside the mixing tube will generate equivalence ratio oscillations, and due to the stiffness of the fuel injectors, there will be an inverse relationship with velocity and equivalence ratio. Therefore, when the SPOD mode increases, the strain rate will increase while simultaneously reducing equivalence ratio, both of which will reduce the local HRR. It is also shown that the flame shape is changing slightly as progressing through the instability. This conclusion is also supported by previous work by Dowd and Meadows [29], where equivalence ratio variation during a self-excited combustion instability was studied and showed strong fluctuation. The work noted an in-phase relationship between equivalence ratio and acoustic pressure for a similar combustion setup.
PIM heat release and equivalence results are presented in Figure 19. As the variations during the instability period are minimal, the results from a single phase are presented and representative of the distribution at each phase. Again, it should be noted that only equivalence ratios between 0.7 and 1.3 are presented. The PIM equivalence ratio and heat release figures have significantly different contours. This in contrast to what was seen in Figure 16, where the no PIM results had similar distribution of equivalence ratio and heat release. The difference in the contours highlights the importance of strain rate. A high shear region is created by the flow exiting the opening in the PIM and lower velocity flow above the PIM structure. This will cause less efficient burning and low levels of heat release in the jet. The mixture will burn more efficiently when the flow settles and disperses in the recirculation zone. There are also higher levels of heat release on the surface of the PIM, where a fuel-air mixture diffused through the PIM, at = 0.8. The burning is in the low velocity region next to the flame established in the recirculation zone.
(a)
(b)
The prominent burning locations demonstrate the relationship between the incoming reactants’ velocity and the heat release distribution. In the no PIM case, direct coupling, between the heat release and incoming reaction quantities, is observed. Also, the corner recirculation flow structure allows the acoustic or flow feedback from heat release to directly influence the incoming flow. In contrast, the PIM case shows that this direct coupling is less likely. Flow or equivalence perturbations must first traverse the additional axial distance, where diffusion and additional mixing take place, due to the turning of flow and porous structure. Then, any perturbations in the flow must follow the higher velocity jet exiting the PIM, where burning is less likely due to higher strain rates and mix further as the flow is shed from the main jet into the recirculation zone where it burns. The addition of PIM results in change in the flow field, improved mixing, and reduction in axial velocity fluctuations, and it is less likely that fluctuations introduced to the equivalence ratio or acoustic field can directly couple to the heat release due the redistribution of turbulent energy to lower order modes. This is supported by the dispersed equivalence ratio field for the PIM case, as seen in Figure 19(a). When compared to the equivalence ratio field in the no PIM case, the PIM equivalence ratio is less concentrated, meaning there was more time and flow structures supporting mixing.
The elimination of recirculation zones and redistribution of energy into lower energy flow structures with PIM addition supports the mitigation of instabilities, but as mentioned previously, an important effect of PIM is the natural damping of acoustic waves, which will also improve system stability. Acoustic waves have been shown to directly result in fluctuation of equivalence ratio coming into the combustor [29]. For partially premixed combustion systems, fluctuations in equivalence ratio are a major factor in heat release variations which can be a source mechanism for instabilities. The addition of PIM reduces the strength of the acoustic waves. The additional changes to the flow field improve the mixing and reduce magnitudes of the equivalence ratio fluctuations reaching the flame front. This reduces the magnitude of heat release fluctuations and thus reduces the energy addition to the system. The increased damping, which needs to be overcome to form an instability, makes it less likely for a system with PIM to become unstable.
To examine how the heat release is varying with acoustic pressure, Figure 20 presents the spatial sum of heat release at each phase. The heat release is then subtracted by the integral heat release and normalized by dividing by the integral heat release for both the PIM and no PIM results. Heat release is compared to the acoustic pressure and equivalence ratio fluctuations with phase for the no PIM case. Since the variations for PIM are small, the sum of the PIM spatial results will be similar to the integral heat release for all phases. The acoustic pressure and equivalence ratio are subtracted by their respective mean values and then normalized by the max acoustic pressure and equivalence ratio fluctuation, respectively. A single pressure curve is used to represent the pressure change throughout a period of instability for both PIM and no PIM results. As shown in Figure 7, the magnitude of pressure is significantly different between the no PIM and PIM pressure results, but the normalized acoustic pressure definition at each phase will produce a similar curve.
The results show an in-phase relationship between heat release and equivalence ratio fluctuations. A slight shift in phase of about 36 deg between the acoustic pressure and heat release is also noted. Since the heat release is in phase with the pressure fluctuation, the Rayleigh index is positive and adding energy to the system. This energy addition still needs to be larger than the damping in the system to form an instability. The significance of the difference in acoustic perturbation is again made clear. The effect of the large reduction in acoustic perturbation has been shown to reduce equivalence ratio fluctuation [29] and heat release fluctuations, as shown by the results in Figure 20. PIM also adds damping to the system, resulting in a larger amount of energy needing to be produced by an in-phase relationship to form an instability.
To further study the spatial contribution to the instability formation or mitigation in the no PIM or PIM cases, respectively, the Rayleigh index is presented in Figure 21. The results for Rayleigh index use the heat release results presented in this section, with the pressure reconstruction from Figure 9. The results for the pressure reconstruction of the PIM case were not presented above but follow similar methodology and produce a similar mode shape with a significant reduction in magnitude. A positive Rayleigh index denotes a location that contributes to the instability and a negative value signifies damping. The Rayleigh index is only presented in regions where heat release was nonzero in all phases. The results are averaged temporally over the phases, through equation (3). Overlaid on each of the contours is the time averaged velocity vectors for the respective case. The no PIM results have two locations denoted by a diamond (Location 1) and square (Location 2) used in the following discussion. The PIM results are multiplied by 2.5 to lie on a similar scale as the no PIM results.
(a)
(b)
The first thing to note is the positive Rayleigh index in the no PIM results and lack of any positive values in the PIM case. This means that in the PIM case, there are little to no spatial locations where the acoustic pressure and the heat release are in phase and contributing energy to an instability. The overlaid velocity vectors confirm that the prominent burning locations are only in areas of lower strain on the flame. In the no PIM case, the Rayleigh index shows that most of the heat release near the quartz is in phase with acoustic perturbations and contributing to the instability. In contrast, the burning that occurs in the jet is damping the instability. This is explained in the next paragraph by the mean velocity and its inverse relation to the convective time delay. The convective time delay is defined as the amount of time it takes the fuel to convect from the injection site to the flame front, which has been previously shown to be a critical parameter in the system stability [4, 50].
The convection time delay is calculated in two parts, first in the mixing tube and then in the combustor. The instability frequency is driving equivalence ratio fluctuations at the injection site, supported by in-phase trends in Figure 20. The equivalence ratio perturbations in the mixing tube convect to the combustor at a mean velocity determined by the mass flow in the mixing tube to approximate the convective time delay. The distance from injection to the swirler and the mean velocity in the mixing tube results in a time delay of 4.1 ms. To approximate the time delay in the combustor, the absolute distance traveled from the swirler exit to the two burning locations identified in Figure 21 must be calculated. The unstable high heat release region at = 0.9 and = 0.8 (Location 1) and the damping burning in the jet at = 0.4 and = 0.4 (Location 2) are used. The streamline passing through the two locations were used to determine the time delay from the dump plane to the physical location. The results are summarized in Table 2, which specify the phase difference for 5 different phases in the instability period at the two locations. The phase difference is calculated by dividing the convection time delay by the period of instability and multiplying by 360 deg, and is presented in degrees. This is only presented for the no PIM case to explain spatial damping and instability regions.
The difference between the two burning locations highlight why one is damping the instability while the other is contributing to the instability. The in-phase relationship at the higher heat release area with the injection site variations is confirmed through the convective time delay analysis here. In the jet at location 2, the convective time delay shows that the heat release is completely out of phase with the acoustic perturbations at the injection point, confirming the damping region presented in Figure 21. The mechanism for generating the instability in this configuration is shown in Figure 22. The pressure disturbances occurring in the mixing tube generate equivalence ratio fluctuations which then convect to the flame and generate a heat release oscillation. If the pressure and heat release are in phase, energy is added to the system and the amplitude of pressure grows creating a feedback mechanism. If the pressure and heat release are out of phase, damping occurs.
To explain the differences in Figure 21 and the improved stability with the PIM addition, the relationship between pressure and heat release is further investigated. Figure 23 shows the sum of the heat release fluctuation and acoustic pressure fluctuation for each phase of the PIM cases. The fluctuation at each phase is calculated by subtracting the mean of all phases by the value at each phase. As expected, there is little variation in PIM heat release. A sinusoidal function would be expected if any coupling were present. Instead, a random fluctuation around a mean value is observed. Comparing pressure and heat release at each phase, an equal number of phases produce a negative Rayleigh index as positive. The final negative Rayleigh index in Figure 21(b) is a result of the summation and the magnitude of the damping outweighing any coupling. As mentioned in Figures 15 and 18, the amount of heat release is dependent on equivalence ratio and strain rate or local velocity. The equivalence ratio and local velocity are more uniformly distributed due to the improved mixing by PIM addition. This causes a reduced coupling with the acoustic perturbations and produces less fluctuation around mean values, resulting in a more stable system. Positive Rayleigh index is still possible, as shown in Figure 21, but the large reduction in fluctuation reduces the level of any coupling between acoustic and heat release distribution in the domain.
6. Conclusions
Spatially resolved, phase averaged equivalence ratios and heat release measurements were performed on an atmospheric, partially premixed combustor. The system was operated at a self-excited thermoacoustic instability and porous inert media (PIM) were to passively mitigate the instability. The addition of PIM caused a significant reduction in peak frequency SPL of 38 dB. Time-resolved PIV measurements were used to investigate the flow field, and results show strong coupling of heat release, equivalence ratio, and SPOD modes without PIM, and a lack of coupling with PIM. The heat release rates were found to be in phase with the equivalence ratios and pressure fluctuations, whereas the SPOD dominant mode was out of phase. The fluctuations observed with PIM did not couple with the instability frequency, and the PIM redistributed the turbulent energy to higher order SPOD modes.
The spatial measurement of heat release elucidated the role of PIM in reducing the coupling between acoustic and heat release perturbations. The distributions of heat release and velocity in the domain were shown to be more uniform with PIM addition. It was shown that, while the no PIM case also had regions of damping, there were larger regions of the flow field that was adding energy to the system. The pressure fluctuations in the mixing tube cause equivalence ratio oscillations to convect downstream to the flame. The location of maximum heat release is impacted by the local strain rate, which causes local quenching if above a critical value, and it impacts the overall convective time delay. The overall stability of the system is determined by the convective time delay between fuel injection and flame location/heat release location. In conclusion, PIM improves thermoacoustic stability, dampens acoustic perturbations, alters the flow field by disrupting the recirculation zones, distributes the energy containing flow structures, improves mixing, and reduces the paths that heat release can couple with incoming fluctuations.
Data Availability
Data are available upon request to the corresponding author.
Conflicts of Interest
The authors declare no conflicts of interest.