LEAN-BLOW-OUT SIMULATION OF NATURAL GAS FUELED, PREMIXED TURBULENT JET FLAME ARRAYS WITH LES AND FGM-MODELING

To ensure compliance with stricter regulations on exhaust gas emissions, new industrial burner concepts are being investigated. One of these concepts is the matrix burner, consisting of an array of premixed, non-swirling jet flames. For the design of such burners, the prediction of fundamental burner properties is mandatory. One of these essential quantities is the lean blowout limit (LBO), which has already been investigated experimentally. This study investigates the possibility of numerical LBO prediction using a tabulated chemistry approach in combination with Large-Eddy-Simulation turbulence modeling. In contrast to conventional swirl burners, the numerical description of blowout events of multi jet flames has not yet been studied in detail. Lean blowout simulations have therefore been conducted for multiple nozzle variants, varying in their diameter and global dump ratio for a variety of operating conditions, showing their general applicability. A procedure to induce LBO is introduced where a stepwise increase in total mass flow is applied. LBO is determined based on the temporal progress of the mean reaction rate. A comparison with measurements shows good agreement and demonstrates that the procedure developed here is an efficient way to predict LBO values. Further investigations focused on the flame behavior when approaching LBO. The flame shape shows a drastic change from single jet flames (stable conditions) to a joint conical flame approaching LBO, which increases in length for increasing inlet velocity, showing the importance of jet interaction at LBO. ∗Contact: alexander.schwagerus@kit.edu NOMENCLATURE

general applicability. A procedure to induce LBO is introduced where a stepwise increase in total mass flow is applied. LBO is determined based on the temporal progress of the mean reaction rate. A comparison with measurements shows good agreement and demonstrates that the procedure developed here is an efficient way to predict LBO values. Further investigations focused on the flame behavior when approaching LBO. The flame shape shows a drastic change from single jet flames (stable conditions) to a joint conical flame approaching LBO, which increases in length for increasing inlet velocity, showing the importance of jet interaction at LBO

MOTIVATION AND INTRODUCTION
A primary focus of current gas turbine research is the reduction of toxic exhaust emissions. This is usually achieved by operating combustion systems at very lean conditions. There, however, the flames become more susceptible to combustion instabilities, which can lead to lean blowout (LBO) of the flame, a key issue in combustion chamber design [1]. To extend the operating range and load flexibility of gas turbines, which are currently based on premixed swirl burners, new burner concepts are being researched and investigated. One of these new burner concepts is the so-called matrix burner, which consists of a matrix of non-swirling jet flames. While the blowout of jet flames has already been investigated in detail, especially for bluff-body and swirl stabilized flames [2][3][4][5], the effect of interacting nonswirling jets on the LBO-limits has only been the focus of recent studies.
One experimental study about the matrix burner has been conducted by Bhagwan et al. [6], who investigated the influence of the geometric variation of the hole diameter of the matrix on the LBO limits and the non-reactive flow field. They discovered that the nozzles at a constant global dump ratio behave geometrically similar and thus the LBO limits can be described by a correlation based on the Peclet number. Further LBO experiments on the matrix burner were performed by Weis et al. [7], who examined the correlation of the LBO limits at an additional variation of the dump ratio. They observed that nozzles of higher dump ratios have an increased operating range and using a correlation based on the Damköhler number they were able to predict the geometric and thermodynamic effects on the LBO limits. Due to the high costs of experiments and the limited access to physical quantities inside the combustion chamber, this paper explores the possibility of numerically investigating the blowout of the matrix burner.
Numerical simulation of the lean blowout represents a big challenge, since both a detailed description of the turbulence field due to its influence on the flame velocity, as well as combustion models, which take the flame extinguishing mechanisms into account, are needed. Of particular interest for industry-scale simulations are flame-generated-manifold (FGM) models, which provide a strong reduction in the number of reactive equations required and therefore in reduction of needed computing power. The blowout phenomenon in general has already been numerically investigated by multiple research groups in recent years using various turbulence and combustion modeling approaches. Akhtar et al. [8] investigated a turbulent premixed single jet flame at an increased preheating temperature and pressure using a combination of an FGM model and Reynolds-Averaged-Navier-Stokes (RANS) turbulence modeling. They found a dependence of the flame position on the inlet velocity or turbulence and were able to calculate LBO limits with an accuracy of 20% compared to the experimental blowout velocity. A comparison of the flamelet/progress variable (FPV) and thickened flame model in conjunction with LES were performed for a swirl stabilized burner by Ma et al. [9]. In their work they suggest using the integrated heat release as an early warning signal for detecting flame blowout. Both models underpredicted the blowout limits by 25% and 20% respectively. Nassini et al. [10] used an FGM model, based on an extended TFC approach in combination with LES modeling, to investigate the flame behavior during LBO and the fragmentation of the flame. Deviations of the two investigated operating conditions from the experimental LBO point were 5% and 10% of the calculated equivalence ratio.
These investigations show that FGM models can be successfully used for the investigation of the blowout limits for different burner types. However, they have been usually only applied for a few operating conditions and were mostly limited to single burner systems. Especially for multi-burner systems without active flow stabilization a concrete proof of applicability is missing. This study aims to close this gap by numerically calculating the blowout limits for the matrix burner using an FGM model in conjunction with a JPDF-and LES turbulence model. The numerical calculations will provide greater insight about the flame behavior towards LBO.

Description of the combustion system
The matrix burner was developed and investigated at the Engler-Bunte-Institute, Division of Combustion Technology at Karlsruhe Institute of Technology [6] [7]. The combustion system consists of a combustion chamber with a hexagonal base and a nozzle, which contains cylindrical boreholes parallel to the direction of the flow. The boreholes are arranged in a specially defined pattern to ensure equal distances between adjacent holes and to the wall. Different nozzles, varying in their number of holes, the hole diameter d Hole and the global dump ratio DR (ratio of the combustion chamber area divided by the sum of the area of the holes) have been used to investigate the influence of geometrical scaling. It is noteworthy, that for a constant dump ratio the distance between holes, normalized by the hole diameter, is constant. Details of the nozzles are listed in Table 1, while a schematic sketch is shown in Fig. 1. The nozzle D 5 employs the same hole diameter and number of holes as the nozzle D 2 ,

V03AT04A026-2
Copyright © 2021 by ASME; reuse license CC-BY 4.0 but the holes are spread over a bigger combustor area, increasing the jet flame distances. The length of the combustion chamber was varied in the experiments for each nozzle to ensure that the flame was burning fully inside the combustion chamber and ranged from 32 to 40 d Re f . The combustion system is operated with a premixed natural gas-air mixture, preheated to a temperature T 0 at atmospheric pressure conditions. The composition of natural gas used in the experiments was measured via gas chromatography and is listed in Table 2.

Turbulent combustion modeling
The numerical prediction of flame instability and flame blowout has always represented a very challenging task for turbulent combustion CFD modelling. Due to the transient nature of the blowout and the small incremental changes in operating conditions required to determine the LBO point accurately, a long simulation time is needed. Therefore, it is crucial to reduce the number of computed transport equations as much as possible. In case of reactive simulations with detailed reaction mechanisms, the computation of the large amount of transport equations for the many species is a major problem. One possibility to reduce the calculation effort is offered by the FGM methods [11], which allow a reduction of needed transport variables by reducing the possible number of reaction paths. In the model used in the current work, the number of reaction paths is reduced by using a single reaction progress variable c. The definition is shown in Eq. 1 and is based on a characteristic variable Y C that relates a sum of mass fractions to their values in unburnt (u) and burnt (b) state. In the current model the characteristic variable is defined as the sum of the mass fractions of CO and CO 2 , the main products of the carbon oxidation. In order to map the reaction progress to the thermodynamic and composition space a model reaction system is needed, which was chosen to be the one-dimensional, premixed flame.
A transport equation is solved for the characteristic variable Y C and is shown in Eq. 2. For the Schmidt number a constant value of Sc t = 0.7 has been used. In order to close the equation, a calculation method for the LES filtered source termω Y C is required. For this purpose, the joint presumed probability density function model (JPDF-model) is used, which has been already applied successfully to a range of different burners [12][13][14]. The basic idea behind the JPDF model is that the mean source term can be calculated using a known probability density function (PDF) of the reaction progress (see Eq. 3).
Assuming the principal form to be a beta distribution, the function of the PDF can be determined with the help of two variables, the average characteristic mass fractionỸ C and a subgridscale (SGS) variance Y C of latter. As an exact equation for the closure of the SGS variance cannot be provided, an analytical approach analogous to the well-known Smagorinsky model in LES turbulence modeling is used in the current work (Eq. 4). In this

V03AT04A026-3
Copyright © 2021 by ASME; reuse license CC-BY 4.0 equation, the value of the constant CỸ C is set to 0.1 and ∆ denotes the filter width of the LES-model. Using the variablesỸ C and Y C , the PDF can be calculated and the source term can be integrated using an algorithm based on the works of Liu et al. [15]. To speed up the calculations, the mean reaction source term is precalculated and tabulated in a multidimensional table prior to the main simulation as a function ofỸ C and Y C for each investigated equivalence ratio.

Peclet correlation
Due to the multitude of geometric and thermodynamic parameters influencing LBO, dimensionless parameters are often used to correlate LBO limits. The Peclet number, which describes the ratio of convective to diffusive heat transport, is a parameter known to be applicable to this problem in other systems. With the help of the two Peclet numbers Pe U and Pe S L , a characteristic critical blowout velocity u char,LBO can be correlated to the thermal diffusivity a and the laminar flame speed S L according to Eq. 5 and 6. The quantities C and n describe burner specific constants which, as evaluated by Bhagwan et al. [6], can also be used for the current matrix burner for all nozzles of a constant dump ratio if the volumetric inlet velocity u in is used as the characteristic velocity. As measured values show some scatter due to the probabilistic character of the underlying turbulent flow, it is useful to compare numerically calculated LBO limits not with the individual measured data but with the experimentally determined Peclet correlation, which provides an average of the measured values. Furthermore, the two Peclet numbers allow for a non-dimensional representation of the phenomenon and thus the visualization of many different operating conditions and nozzles in one diagram.

Numerical Setup and methods
All numerical investigations that were performed for this study used the open source C++ toolbox OpenFOAM. Large Eddy Simulation (LES) turbulence modelling, which can be considered state of the art for academic investigation is used to investigate the transient behavior towards LBO. The closure of the subgrid stress tensor has been achieved through the WALE subgrid model [16], which was often found to better model the turbulence, especially in near wall regions, compared to classical models like the standard Smagorinsky model [17,18]. For both, the discretization of the temporal as well as the spatial derivatives, second order discretization schemes were used.
Computational domains for each investigated nozzle have been defined and one is exemplary shown for the nozzle D 2 in Fig. 2. The domains have been discretized with around 1.2 million cells each and can be divided into three subdomains: The short part of the multiple inlets on the left, which lead flow into the combustion chamber and lastly into an outflow domain. The inlets and the combustion chamber have been discretized using a fine grid with a cell edge size around 1 mm, while the outer domain is coarsely discretized with a cell edge size around 4 mm. The hole spacing, a critical geometrical quantity in the domain, has been resolved with a minimum of five cells. The outflow domain has been added to reduce the risk of numerical pressure wave reflections at the outlet boundary back into the combustion chamber, which could induce a false LBO. Additional grid coarsening could not be applied, as the position of the flame strongly changes toward LBO, requiring a fine resolution everywhere in the combustion chamber. The cell resolution has been selected based on a grid study in non-reactive simulations in order to resolve the mean velocity profiles sufficiently. The dimensionless wall distance y+ has been evaluated in reacting simulations and was found to be less than 10, so no wall functions have been employed. The Pope criterion [19] is fullfiled for over 80% of the cells in the combustion chamber.
As an LES simulation requires the time-resolved description of the incoming turbulent structures at the domain inlet a turbulence generator proposed by Klein et al. [20] has been implemented [21]. This method is based on digital filtering of a series of uncorrelated random data to generate correlated velocity fields according to user-defined turbulence properties. The needed turbulence values have been fitted in several nonreactive simulations of the combustor to match the measurements of mean velocity and turbulent kinetic energy reported by Bhagwan et al. [6].
The matrix burner was, like in the previous experiments, examined under atmospheric pressure conditions fueled by a preheated natural gas -air mixture. The required chemical tables were generated with the detailed GRI 3.0 mechanism [22] using the open-source kinetics software Cantera assuming the model reaction system of a perfectly premixed one dimensional flame with a detailed description of the diffusive fluxes. The system

V03AT04A026-4
Copyright © 2021 by ASME; reuse license CC-BY 4.0 includes both the mutual influence of chemical reaction, energy and mass transport on a molecular level, which are essential for the stability of a flame. Since the one-dimensional flames were calculated without enthalpy sinks, the resulting chemical tables can be considered adiabatic.

Numerical procedure and LBO criterion
Starting from a stable flame configuration, there are several possibilities to induce lean blowout. The fuel mass flow can be reduced at a constant air flow rate, as it was also done in the previous investigations of the matrix burner. With this method the fuel-air ratio is decreased, while keeping the total flow rate nearly constant. This procedure can be carried out numerically, too. However, an additional transport equation for the mixture fraction must be calculated in the transient simulation. Furthermore, to ensure that a uniform mixture is achieved after each change of the fuel-air ratio, a time-consuming simulation period is needed. A numerically more advantageous alternative is to keep the fuel-air ratio constant in each simulation and only to increase the total mass flow rate of the mixture, which can be specified by the volumetric inlet velocity u in . This method was employed in the current study. In particular, the volumetric mean velocity is increased from stable flame conditions until LBO is determined. In literature LBO in numerical calculations is often recognized via a "sudden" drop or increase in global variables, like mean temperature or simply just optically. The current study proposes a more precise criterion, based on the volume averaged reaction rateω Y C inside the combustion chamber.
For a stable flame, which converts the entire fuel within the combustion chamber, the integral of the fuel reaction rate is equal to the incoming mass flow of fuel (Eq. 7). Furthermore the averaged fuel conversion rate is proportional to the mean reaction rate of the characteristic mass fraction Y C , while the mass flow of fuel is dependent on the inlet velocity u in (Eq. 8). In case of the LBO simulations, where the incoming mass fraction of the fuel Y Fuel is constant, this implies that as long as the flame burns stably inside the combustor the mean reaction rateω Y C must increase with increasing inlet velocity u in (Eq. 9). This holds true until the flame begins to get blown out, where the mean reaction rate begins to drop. So, for this investigation the first velocity increase step which shows a decrease in mean reaction rate is considered as LBO. The used blow out velocity for later plots is then the average velocity between the last stable operating point and the operating point, where LBO is detected.
General flame shape First simulations investigated on the representation of the flame shape at fixed operating conditions. In order to give an example Fig. 3 shows slices of the time-averaged reaction progress field for three different fuel-air equivalence ratios φ at a fixed inlet velocity u in using the nozzle D 2 . At the top of the figure a near stoichiometric flame is shown, while the flames were operated leaner towards the bottom. It is noticeable that the flame shape shows a strong change depending on the equivalence ratio: While the upper, near stoichiometric flame clearly shows distinguishable single flame jets, a decrease of the fuel-air ratio to φ = 0.68 leads to a thickening of the flame front, which causes a partial merging of the flame fronts while the flames are slightly elongated. Under these conditions, the individual jet flames can still be identified, but their mutual influence increases. This effect becomes even stronger with an increase to φ = 0.56, where finally a uniform flame front is formed in which the individual

V03AT04A026-5
Copyright © 2021 by ASME; reuse license CC-BY 4.0 jet flames are no longer distinguishable. Under these conditions, the flame is significantly longer and forms a conical flame, with the tip reaching almost the end of the combustion chamber. For even leaner conditions no flames could be stabilized. A qualitatively similar behavior was also found for the other nozzles and is also observed when the equivalence ratio is kept constant while the inlet velocity u in is increased. It is obvious that the flame-to-flame interaction becomes more and more important for the matrix burner when approaching LBO. The merging of the flame fronts with a decrease of the equivalence ratio is a phenomenon that has also been observed in other multi-burner systems, such as Ciardiello et al. [23], who have experimentally studied a combustion chamber with several swirled bluff body stabilized flames.

Lean blow out simulations
This section discusses the results of the LBO simulations. A stepwise increase of the inlet velocity u in has been used to induce LBO. The time interval of constant velocity was adjusted in each simulation to ensure that a stationary flame could form again. Snapshots of the reaction progress field in the combustion chamber of an exemplary LBO simulation of the nozzle D 3 are shown in Fig. 4 . Four points in time of the LBO simulation are shown, characterized by the current inlet velocity u in , increasing from top to bottom. The uppermost image shows a flame shape similar to Fig. 3 for near-stoichiometric conditions, which are characterized by single isolated jet flames. This type of flame shape is found when the operating conditions are far from LBO. When the inlet velocity is increased (second image) a short unified conical flame is formed. With a further increase of the velocity to 32 m/s (third image) a strong elongation of the cone of the flame is observed and already reaches the end of the combustion chamber. This is the last stable operating point in this simulation. After the next velocity increase, as shown in the last picture, a sudden blowout of the flame takes place, and the combustion chamber is almost completely filled with unburnt gas. Two remaining types of reaction zones can still be seen: On one hand, between the inlets, where complete extinction cannot be detected in this short time due to the long residence time there and the associated slow response. The fact that the flame is already blown out while reaction zones are still present near the nozzle indicates that the flame is not stabilized by the flame root, but by the leading edge of the flame. On the other hand, there is also still ongoing reaction at the edges of the outlet of the combustion chamber, which can be identified by red zones near the right side of the picture. This can be explained by the inclusion of the outer domain (which is not shown in the pictures): This outflow zone causes a slow-down of the gas mixture due to the cross-sectional expansion. In reality, a dilution of the fluid at these positions by entrained air would occur, preventing such flame stabilization. Since the solver assumes a perfectly premixed composition, the flame always stabilizes in these low velocity regions and therefore is able to occasionally move upstream the edges of the combustion chamber. Despite this fact, a clear blowout of the flame can be identified, indicating that the flame stabilization in the outer domain is no problem for the simulation of LBO.
The simulation shows a strong change in flame shape from single flames far from LBO to increasingly longer conical flames and a sudden blowout when a critical velocity is exceeded. Although the LBO point is optically visible, a more accurate and general detection method is desirable. Since the flame blowout is not only a local phenomenon, a globally defined tracking variable is required. For this purpose Fig. 5 shows the current inlet velocity as well as various possible tracking variables in form of volumetrically averaged values over the simulation time. The simulation started at an inlet velocity of 12 m/s and was increased stepwise by 5 m/s every 0.04s. In these diagrams, the black lines represent the transient evolution of the values, while the red curve values are additionally averaged over the time interval of constant velocity. An often-used parameter to detect LBO is the mean combustion chamber temperature, which is shown in the second diagram. The temperature is starting to decrease gradually due to the elongation of the flame, resulting in an increase of the volume of unburnt gases. Then at 0.2s a steeper drop is detected due to LBO, which results in a strong reduction of the average temperature compared to the previous operating point. Although the drop at the LBO point is more significant compared to the drop in the velocity steps before, it is difficult to clearly identify the point of blowoff. It should be noted that

V03AT04A026-6
Copyright © 2021 by ASME; reuse license CC-BY 4.0 The mass fraction of OH is almost constant up to 0.16s, then dropping slightly. At 0.2s also the mean mass fraction drops significantly, indicating the lean blowout. A short spike can be detected at 0.21s. The reason is a pressure wave reflection due to the extinction of the flame, which results in a temporal thickening of the flame front. A more reliable way to determine the LBO point is to use the volume averaged reaction rateω Y C , for which the transient values are shown in the bottom diagram. As it has already been derived, for a stable flame an increase in the reaction rate is expected with increasing inlet velocity. This can also be seen in the diagram, which shows a stepwise increase up to 0.2s. It should be noted that the reaction rate diagram shows higher fluctuations compared to thermodynamic quantities such as temperature and mass fractions, so that it must be ensured that a sufficiently long averaging time is provided. However, if this pre-requisite is fulfilled, the method allows for a clear and general identification of LBO by examining the change of the gradient from an increase to a decrease, which in this example is clearly observed at 0.2s. This criterion for the determination of LBO has been successfully applied to all simulations. The examination of the flame shape shows that close to LBO no longer single jets, but a joint flame cone with a lifted flame root defines the flame shape. It can therefore be assumed that not the hole velocity u in is the characteristic velocity of the matrix burner, as it would be for single jet burners, but the mean combustion chamber velocity u comb = u in /DR. With this a modified Peclet number Pe U,DR can be calculated according to Eq. 10.
Twelve LBO simulations have been performed for different equivalence ratios and nozzles to show the general applicability of the numerical setup. These simulations always showed a similar flame behavior as already described above. Additionally, to investigate the cell size sensitivity on the calculation of LBO, one simulation has been repeated with a numerical grid with half the

V03AT04A026-7
Copyright © 2021 by ASME; reuse license CC-BY 4.0 . Furthermore vertical error bars has been added to each numerically calculated value to indicate the uncertainty due to the stepwise velocity increase. As already reported by Bhagwan et al. [6], the experimental values of one dump ratio can be correlated using the Peclet correlation (Eq. 6). In addition, Peclet fitting parameters based on the measurements of the nozzle D 5 have been calculated. The resulting fits of the experimental values are plotted as solid black lines. The used Peclet constants are listed in Table 3.
In the double logarithmic plot, the Peclet curve forms a straight line for each DR, each of which is located in a very similar range. However, they differ in their slope, with a higher slope at DR = 6. The different gradient can probably be explained by the presence of the turbulent conical flame near LBO. The flame velocity of this is significantly influenced by the turbulent diffusivity inside the combustion chamber. Turbulence is generated by the expansion of the flow cross-section and the resulting shear zones. This strength of the shear zone is influenced by the dump ratio and might be the explanation for the different slopes. A comparison of the simulated LBO limits shows very good agreement with the experiments for the whole investigated Pe S L range. The values of the nozzles D 2 and D 3 reside almost perfectly on the Peclet correlation, while there are only slight deviations for the nozzle D 1 . This can be also seen in the mean deviation of the numerical LBO values compared to the experimentally determined Peclet curve, which amounts to 11% and is comparable to the 6% mean deviation of the experimental values itself. Furthermore, two additional simulations have been conducted, that also predict the change in flame stability of the nozzle D 5 due to the higher dump ratio. While one operating point predicts a slightly lower flame stability limit compared to

V03AT04A026-8
Copyright © 2021 by ASME; reuse license CC-BY 4.0 the experiments, there is almost perfect agreement for the other point. This shows that the numerical determination of LBO is a valuable alternative to experimentally determined correlations like the Peclet criterion due to its applicability for different geometrical set-ups of the matrix burner. It should be noted that all results are obtained without the need of tuning specific model parameters, showing the general applicability. Furthermore, the simulations have been conducted assuming adiabatic walls, indicating that heat losses might not be strong enough to impact flame stability in the current multi jet burner system.

Conclusion
This study presents a comprehensive numerical study of the blowout limits of the matrix burner, which includes a variety of operating conditions and nozzles. A numerical setup utilizing LES turbulence modeling and a FGM-model has been examined. The simulations show a strong change in flame shape, depending on the operating condition. For stable flames (low velocity or high fuel-air ratio) single isolated jets are observed, while forming a conical flame for more unstable conditions, which is elongating, while approaching LBO. After exceeding a critical velocity, a sudden blowout of the flame is observed. For the determination of the LBO a criterion based on the transient course of the mean reaction velocity was chosen, which shows a sign change in the slope at LBO. The numerically determined LBO limits for all nozzles over the whole investigated range correspond very well with the experimental values. The success of LBO prediction with the assumption of no heat losses indicates that they do not significantly impact flame stability in the current setup. Future investigations will target the limitations of the used numerical setup. For gas turbine applications it will be of particular interest if there is still good agreement at higher pressure levels.