# Large Eddy Simulation of a Swirl-Stabilized Pilot Combustor from Conventional to Flameless Mode

Authors: Ehsan Fooladgar; C. K. Chan
Journal: Journal of Combustion
 (2016)

---

## Abstract

This paper investigates flame and flow structure of a swirl-stabilized pilot combustor in conventional, high temperature, and flameless modes by means of a partially stirred reactor combustion model to provide a better insight into designing lean premixed combustion devices with preheating system. Finite rate chemistry combustion model with one step tuned mechanism and large eddy simulation is used to numerically simulate six cases in these modes. Results show that moving towards high temperature mode by increasing the preheating level, the combustor is prone to formation of thermalNOx with higher risks of flashback. In addition, the flame becomes shorter and thinner with higher turbulent kinetic energies. On the other hand, towards the flameless mode, leaning the preheated mixture leads to almost thermal NOx-free combustion with lower risk of flashback and thicker and longer flames. Simulations also show qualitative agreements with available experiments, indicating that the current combustion model with one step tuned mechanisms is capable of capturing main features of the turbulent flame in a wide range of mixture temperature and equivalence ratios.


---

## Body

## 1. Introduction

Ever increasing global energy consumption and environmental concerns combined with the lack of energy resources have put the designers of combustion devices to a difficult test in order to come up with new technologies that are more energy efficient and less polluting. One common method to increase energy efficiency in almost all combustion systems including lean premixed (LPM) combustion is increasing the mixture temperature using recovered exhaust heat directly. Employing this method in LPM combustion not only increases the efficiency of the system but also improves combustion stability and flammability limits. As outlined by Huang and Yang [1], LPM combustion is the most promising technology for environmentally friendly combustion systems since operating under fuel lean conditions can have low emissions and high efficiency. Thermal nitric oxide formation is reduced because flame temperature is generally low and, for hydrocarbon fuels which are leaned by excess air, hydrocarbon and carbon monoxide (CO) emissions are reduced due to complete burnout of fuel. Unfortunately, as explained by Dunn-Rankin [2], achieving these improvements and meeting the demands of practical combustion systems are complicated by low reaction rates, extinction, instabilities, mild heat release, and sensitivity to mixing.Panoutsos et al. [3] studied the effect of preheating of air/methane mixture up to 400°C on local equivalence ratio in a swirl-stabilized model gas turbine combustor using chemiluminescence sensor. They concluded that air preheating was beneficial for the premixing of fuel and air and as the temperature of the combustion air increased from 25°C to 400°C, the flame became shorter and moved upstream until it was stabilized at the boundaries of the inner recirculation zone. Seo [4] and Huang and Yang [5] presented experimental and LES results on flame structure in lean premixed combustors due to increasing inlet temperature. Huang and Yang reported that inlet temperature and equivalence ratio are key parameters determining stability characteristics of their combustor, where slight increase in inlet temperature above a critical value causes abrupt instability in the combustor. Foley et al. [6] investigated experimentally flame shapes for preheat temperatures ranging from 366 K to 533 K with equivalence ratios ranging from 0.40 to 0.70. They found that transition from one flame configuration to another is essentially due to the flame extinction phenomenon and that sensitivity of these transition points to fuel/air ratio and preheat temperature can be reasonably captured with extinction strain rate calculations.A nonpilot assisted combustion device works in conventional, high temperature combustion (HiTC) [7] and flameless or Moderate or Intense Low-Oxygen Dilution (MILD) combustion modes based on premixture temperature and composition [8]. Despite its benefits, increasing the reactants temperature above a certain level gives rise to thermal NOx formation in conventional and HiTC modes, thereby limiting the potential of this method to produce less NOx. The risk of flashback also grows in these regions. Low heat release ultralean combustion with significant preheat, usually referred to as flameless combustion, could thus be a solution to these problems [9]. However, no significant work has been done to compare the flame and flow evolution of a practical swirl flame in these different combustion modes.The objective of this paper is to investigate the flame and flow structure of a swirl-stabilized pilot combustor in conventional, HiTC, and flameless mode using finite rate chemistry combustion model and large eddy simulation (LES). In this paper, the flameless mode is achieved by preheating the ultralean premixture indirectly using an external preheater and/or a recoupaerator and is different from the common flameless concept in which the premixture is heated and diluted using exhaust gas recirculation. As experiments in premixed combustion with inlet temperature higher than 700 K is of safety concern, numerical simulation, especially LES, is the preferred choice for studying turbulent premixed combustion in this range of preheating. The present paper also aims to give a better insight to incorporate new technologies like MILD and HiTC in swirl-stabilized combustors.


## 2. Numerical Procedure

Reacting flows are governed by the balance equations of mass, momentum, species, and energy. The basic idea of LES is resolving the larger turbulent motions in a flow field and modelling only the effects of the small ones. The resolved contributionf- is obtained by applying the spatial LES filter to instantaneous variables f. Filtering the instantaneous governing equations and introducing the Favre filtered variable, f~=fρ¯/ρ-, where over-bars denote spatial filtering, leading to the following equations:(1)∂ρ-∂t+∂∂xiρ-u~i=0,(2)∂∂tρ-u~i+∂∂xjρ-u~iu~j+∂p-∂xi=∂∂xjτ-ij-ρ-uiuj~-u~iu~j,(3)∂∂tρ-Y~k+∂∂xiρ-u~iY~k=∂∂xiρD~k∂Y~k∂xi-ρ-uiYk~-u~iY~k+w˙¯k;k=1,N,(4)∂∂tρ-h~s+∂∂xiρ-u~ih~sk=∂p-∂t+u~i∂p-∂xi+∂∂xiλ∂T~∂xi-ρ-uihs~-u~ih~s-∑k=1Nw˙¯kΔhf,k0,where u is the velocity, ρ is the density, p  is the pressure, τ  is the viscous tensor, Yk is the species mass fraction, hs is the sensible enthalpy, λ is the thermal conductivity, T is the temperature, Dk is the species diffusivity, w˙k  is the species reaction rate, and Δhf,k0 is the species formation enthalpies. The viscous heating term and radiation sink term are neglected in (4) as they are negligible compared to the combustion source term. The species mass flux is described by Fick’s law in the first term on the right hand side of (3). Dk may therefore be expressed as the ratio of kinematic viscosity ν to Schmidt number Sck, which is assumed to be unity for all species in this paper. The dynamic mixture viscosity μ and thermal conductivity are modeled by Sutherland’s law and considering the mixture as a Newtonian fluid, τ-ij=2μS~Dij, where S~Dij  is the deviatoric part of filtered strain tensor defined as S~ij=(1/2)∂u~i/∂x~j+∂u~j/∂x~i.In this system of equations, the second terms in the brackets result from filtering the convective terms and contain subgrid flow physics. These unclosed quantities are modeled using one-equation eddy viscosity model in which the unresolved subgrid scale (SGS) stresses(uiuj~-u~iu~j) are modeled as 2νtS~Dij, unresolved species fluxes (uiYk~-u~iY~k) as (νt/Sct)∇Y~k and unresolved enthalpy fluxes (uihs~-u~ih~s) as (νt/Prt)∇h~s where Sct and Prt are turbulent Schmidt and Prandtl numbers, respectively. The subgrid viscosity is defined as μt=ρνt=ckρΔk, where Δ is the filter size calculated as the cube root of the local cell volume and is provided by an equation for the subgrid kinetic energy [10] k such that(5)∂∂tρ-k+∂∂xiρ-u~jk=∂∂xiμ+μt∂k∂xj-2μtS~Dij∂u~i∂xj-2ρk3∂u~j∂xj-ceρ-k3/2Δ,with dimensionless constants of ck=0.07 and ce=1.048 as proposed by Berglund et al. [11]. Finally, the filtered species reaction rate w˙- is determined via the Arrhenius expression and is modeled as described in Section 3.The CFD code used in this study is OpenFoam which was introduced by Weller et al. [12] and has been validated in a number of studies, including nonreacting flows [13] and reacting flows [14, 15], with generally good agreement with experiments. The code is compressible and employs an unstructured collocated Finite Volume method of Grinstein et al. [13] in which discretization is based on Gauss theorem together with a semi-implicit time-integration scheme. The filtered governing equations, (1)–(5), are discretised using 2nd-order spatial and first order temporal schemes. Time integration is performed implicitly using PISO method, in which pressure and velocity fields are decoupled and solved iteratively, with three PISO correctors. The equations are solved sequentially, with iteration over the nonlinear source terms to obtain rapid convergence, with a maximum CFL number of 0.2.The computational domain of this paper is based on the reduced-scale swirl-stabilized combustor developed by Nogenmyr et al. [16]. In order to create swirling flow, the burner is designed with a radial swirler with tangential inlets enabling a wide range of swirl numbers. Premixed fuel is injected through a premixing system just before both tangential and axial inlets. The benefit of a reduced-scale combustor from a numerical simulation point of view is that the filter size, which is the mesh size in implicit filtering, could be as small as the laminar flame thickness. Geometry of the combustor is given by Nogenmyr et al. [16] and an overview of the domain and mesh used in this paper is shown in Figure 1. Essentially, axial flow of reactants is supplied via a central tube of 16 mm diameter and 160 mm in length. Swirling flow of reactants is induced by four tangential ducts supplied via a plenum with two lateral entries of 9.53 mm diameter tubes to assure flow symmetry. Each tangential duct has a rectangular cross-section of 8 mm in height aligned in the axial direction and 4 mm in width to introduce angular momentum to the axial flow through the burner pipe. To assure a uniform axial velocity distribution as well as generating turbulence, a perforated plate with circular holes of 1.62 mm diameter is placed just upstream of the tangential ducts. Figure 1(b) also outlines the overall flame location in the combustor by showing the mean fuel mass fraction isosurface at 0.01 for the benchmark case as described in Section 4. Computationally, the mesh consists of 1.7 M hexahedral cells, which are body fitted inside the domain. The resolution inside the burner and combustor is 0.5 mm per cell, yielding 32 cells over the diameter of the axial inlet. The computational domain consists of the swirler, the combustor, and a dome shaped region downstream of the combustor contraction to take into account a part of the outside atmosphere. Detailed dimensions of the computational set up are given by Nogenmyr et al. [16].Figure 1
(a) Mesh and (b) numerical domain and the mean flame location for the benchmark case.
(a) (b)Swirl numberS is defined as the ratio of the axial flux of the tangential momentum to the product of the axial momentum flux and a characteristic radius. For a tangential jet-induced swirl flow, it can be estimated as(6)S=∫0RUzUθr2drR∫oRUz2rdr≈πD16whD-wm˙ax/m˙tan+12,where m˙ax/m˙tan is the ratio of axial and tangential mass fluxes, D=16 mm is the diameter of the axial inlet, and w=4 mm and h=8 mm are the width and height of the tangential jets, respectively. Overall, Reynolds number Re is defined as Re=(m˙ax/m˙tan)/μπD/4, where μ is the dynamic viscosity of the premixed reactants.


## 3. Combustion Modeling

The major difficulty in LES of reacting flow using finite rate chemistry (FRC) is modeling the strongly nonlinear filtered reaction rate without specifically considering the geometrical properties of the flame. The FRC models are capable of handling nonunity Lewis number and predicting intermediate species by solving one equation for each species in the chemistry mechanism and capturing extinction and reignition due to retaining Arrhenius chemistry in the modeled reacting rates. The FRC combustion model used in this paper is the partially stirred reactor (PaSR) model which has been widely used in recent years in LES of turbulent premixed swirling flame [15], diesel spray ignition and combustion modeling [17], axisymmetric dump premixed combustor, [18] and piloted lean premixed jet flame [14, 19].The PaSR model is based on the sequential processes of molecular mixing and chemical reactions. Although the PaSR model is suitable for complex chemistry consideration, only one step tuned mechanism is used in this paper as a preliminary study. Since the microscale processes responsible for the molecular mixing, as well as dissipation of turbulent kinetic energy, are concentrated in isolated regions and occupy only a small fraction of the fluid volume, combustion takes place in these well-mixed regions provided that the temperature is high enough. These structures are often organized as tubes or sheets, whose characteristic dimensions are small compared to the LES filter size [14]. Therefore, each LES cell may be divided into reacting fine structures, treated like a perfectly stirred reactor, exchanging mass and energy with its nonreacting surroundings, and dominated by large-scale coherent flow structures. LES PaSR may be mathematically presented as(7)w˙kρ,Yk,T¯=w˙kρ,Yk,T¯w˙kρ-,Y-k,T-w˙kρ~,Y~k,T~=κw˙kρ~,Y~k,T~,where the ratio of the filtered reaction rate to its resolvable counterpart κ represents the reacting volume fraction. To close the subgrid combustion model, κ may be estimated as the ratio between the chemical time scale τc and the total reaction time, such that(8)κ=τcτm+τc.In this paper, the mixing time scale is estimated by replacing the molecular viscosity in Kolmogorov time scaleτK=ν/ε1/2 with effective viscosity ν+νt such that(9)τm=Cmν+νtε=Cm2ck3ceτΔ2+269ceτK2,where ε is the dissipation rate of turbulent kinetic energy and τΔ=Δ/u′ is the time of the subgrid velocity stretch with u′=2k/3 and Cm is a tuning constant taken to be 0.15 in this paper. Finally, the chemical time scale τc is determined by solving the reaction system and finding its characteristic time as the ratio of total molar concentration in each LES cell to the average destruction rate of the reactant species.The turbulent flame speed and thickness of LES PaSR model may be defined assT=Ξsl and lT=Flf, respectively [20, 21], where Ξ=κ(1+νt/ν) is the wrinkling factor, F=(1+νt/ν)/κ is the thickening factor, sl is the laminar flame speed, and lf is the laminar flame thickness. By introducing subgrid Reynolds, Damkohler and Karlovitz numbers as(10)ReΔ=u′Δsllf,DaΔ=slΔu′lf,KaΔ=u′sl3lfΔ1/2,and estimating τc≈ν/sl2 and lf≈ν/sl, the reacting subgrid volume fraction, wrinkling, and thickening factors of this paper’s PaSR model could be dynamically reformulated as(11)κ=11+ΓDaΔ,Ξ=1+CrReΔ1+ΓDaΔ,F=1+CrReΔ1+ΓDaΔ,where(12)Γ=Cm231ce3/2ReΔ+ckce,Cr=ck32.


## 4. Boundary Conditions

In this paper, a premixed flame is simulated on the reduced-scale swirl-stabilized combustor developed by Nogenmyr et al. [16], with a Reynolds number Re of 11,000, swirl number of 0.55, equivalence ratio ϕ of 0.9, and an inlet mixture temperature of 300 K as the benchmark case. Fixed mass flow rate, calculated based on Re and swirl number, and uniform temperature are specified at all inlet boundaries. For all cases, turbulence at the axial inlet are imposed using a turbulence inflow generator with the same turbulence intensity and subgrid turbulent kinetic energy approximated based on previous PIV measurements [16]. No-slip boundary condition is used for velocity at the wall with zero gradients for scalars. Temperature of combustor wall is set with a fixed profile estimated based on the measurements. At the atmospheric boundary, a partially reflective outlet boundary [22] is used for pressure to avoid numerical waves at the boundary and flux conserving condition for velocity.As shown in Figure2, combustion systems may be categorized into conventional, pilot assisted, HiTC, and flameless combustion modes based on mixture temperature Tr, autoignition temperature of fuel, Tign, and the maximum temperature rise, ΔT=Tp-Tr, where Tp is the product temperature [8, 23]. To investigate the effect of preheating on the flame in a wide range of aforementioned combustion modes, the mixture temperature at all inlets is increased from the benchmark case at 300 K to a temperature slightly higher than self-ignition point at 900 K with an increment of 200 K as cases A1, A2, A3, and B1 (referred to as path I in this paper). Keeping all the conditions of the case B1, equivalence ratio of the mixture is then decreased from 0.9 to 0.5 and 0.3 in order to approach flameless zone in the cases B2 and C (referred to as path II in this paper). Since the mass fluxes of all the inlets are constant during preheating, global swirl number of the burner is kept constant. However, Reynolds number decreases because of the increase in dynamic viscosity of the mixture as temperature rises according to the Sutherland’s law for an ideal gas. Flow conditions of six cases simulated in this study are summarized in Table 1.Table 1
Case description.
Case ϕ Re S T r (K) Combustion mode A1 0.9 11,000 0.55 300 Conventional A2 0.9 73,00 0.55 500 Conventional A3 0.9 5,700 0.55 700 Conventional B1 0.9 4,800 0.55 900 HiTC B2 0.5 4,800 0.55 900 HiTC C 0.3 4,800 0.55 900 FlamelessFigure 2
Location of all cases in a combustion mode chart proposed by [8, 23].


## 5. Results and Discussion

### 5.1. Laminar Flame Structure

Incorporating detailed chemical mechanisms in LES is still computationally prohibitive due to the complexity of turbulence/combustion interaction. In addition, one additional equation is required to be solved for species in the reaction mechanism using very fine grids [14, 24]. Thus, the present study uses reduced reaction mechanisms, which is common in many simulations of turbulent flames using FRC models. In this paper, the following reduced one step tuned reaction mechanism (OTM) for lean methane/air combustion is used such that(13)CH4+2O2⟶CO2+2H2Owithw˙=9.49×1011e-Ta/TYCH40.7YO20.8,where Ta= 23,650 K is the activation temperature and the species exponents of methane and oxygen are individually tuned matching the GRI 3.0 mechanism.In order to investigate the accuracy of OTM, a series of 1D freely propagating flames with the same mixture condition as described in Table1 is simulated by CANTERA [25] using GRI 3.0 and one step tuned mechanism used in LES of this study. Figures 3(a) and 3(b) compare the velocity profile of 1D freely propagating flame for the cases in paths I and II, respectively, using both GRI and OTM chemistry. GRI results show that increasing the mixture temperature gives rise to the unstrained laminar velocity sl quadratically, whereas decreasing the equivalence ratio of a highly preheated mixture decreases this velocity linearly. As shown in Figure 3(c), temperature profiles of the 1D flame indicate the same dependency of adiabatic flame temperature on temperature and equivalence ratio of the mixture but both in a linear fashion. These profiles also show that the laminar flame thickness calculated based on lf=ν/Prsl [24] decreases slightly as the preheating level rises and increases significantly when the mixture becomes leaner for the cases in path II. These tendencies may be explained by sl∝μw˙/ρ2 and lf∝μ/ρsl [26]. Combining the variation of laminar flame speed and thickness by calculation of chemical time scale according to τc=lf/sl reveals that chemistry accelerates by increasing the temperature and decelerates by reducing the equivalence ratio in the lean mixture.Figure 3
(a), (b) Velocity, (c) temperature, and (d)NOx mole fraction profiles of 1D freely propagating flame.
(a) (b) (c) (d)Mole fraction profiles ofNO species for all cases are shown in Figure 3(d). It shows that the amount of NOx emission in the product is strongly dependent on product temperature and is important when Tp is more than 1900 K. This illustrates that the cases inside or near the flameless combustion mode, cases B2 and C, are almost NOx-free. Detailed results of the 1D laminar simulation are also presented in Table 2.Table 2
Detailed results of 1D laminar simulation.
Case code OTM GRI 3.0 l f (mm) s l (m/s) τ c (ms) T p (K) l f (mm) s l (m/s) τ c (ms) T p (K) A1 0.060 0.366 0.164 2165 0.067 0.337 0.199 2137 A2 0.055 0.991 0.055 2316 0.064 0.843 0.076 2255 A3 0.043 2.225 0.019 2472 0.056 1.691 0.033 2364 B1 0.032 4.483 0.007 2638 0.047 3.035 0.015 2465 B2 0.094 1.515 0.062 1978 0.087 1.640 0.053 1973 C 0.297 0.482 0.617 1588 0.245 0.583 0.419 1589Comparing the results of two different mechanisms indicates that although the one step tuned mechanism used in turbulent simulation of this study overpredicts the adiabatic temperature and laminar flame velocity and underpredicts the laminar flame thickness of highly preheated cases with equivalence ratio of 0.9 (cases A3 and B1), it predicts adequately the trends of laminar flame velocity, flame thickness, and adiabatic temperature in both paths. It is worth mentioning that the self-ignition temperature in each case is calculated using zero-dimension well stirred reactor simulation for a residence time of 1 second and atmospheric pressure, as suggested by [8], with the same reduced chemistry mechanism exploited in the present LES simulation.


### 5.2. Temperature Field and Turbulent Flame Structure

Figure4 shows the profiles of characteristic temperatures extracted from the 3D domain of each case in two different paths. It illustrates that increasing the preheat mixture temperature increases the maximum temperature of the product Tp in the combustor linearly, whereas decreasing the equivalence ratio has a reverse effect. Increase of product temperature along path I is slightly less than the increase in reactant temperature between corresponding case leading to moderate decline in the maximum temperature rise ΔT, but the temperature drop along path II is significant. It can be seen that the maximum temperature of case B1 is about 500 K higher than case A1 with the preheating level of 600 K while decreasing equivalence ratio from 0.9 to 0.5 results in a drop of 700 K in ΔT. This variation of ΔT leads to a large decrease in the temperature ratio Tp/Tr from 7 to 3 along path I and from 3 to 1.7 along path II, playing an important role in flame-turbulence interaction. Since wall temperature of the combustor is kept constant for all the simulations, increase in exit temperature of the combustor for the cases along path I is attributed to higher heat loss through the walls. Consequently, only increasing preheating temperature of the mixture enhances the overall heat loss in the combustor; therefore, exit temperature of the product does not change significantly from case A1 to B1.Figure 4
Characteristic temperatures.
(a)
Path I (b)
Path IIFigure5 shows the mean filtered temperature contours T~ on a vertical plane across the center of the combustor for all cases. The symbol f denotes time-averaged of variable f. As shown in Figures 5(a)–5(d), the zone with temperature higher than 1800 K expands along path I, increasing the area for formation of thermal NOx in the combustor. In contrast, comparing the isothermal lines of 1800 K in Figures 5(d)–5(f) reveals that leaning the preheated mixture along path II decreases the heat loss through the combustor walls as well as decreases the regions for the formation of thermal NOx. It is worth mentioning that without preheating the mixture, the minimum equivalence ratio that can be achieved with the current burner is around 0.7. Thus, preheating is necessary for expanding the lean flammability limits in order to approach the flameless zone.Figure 5
Mean temperature contours across the center of the combustor.
(a)
Without preheat (b)
200°K preheat (c)
400°K preheat (d)
600°K preheat (e)
600°K preheatϕ=0.5 (f)
600°K preheatϕ=0.3Figure5 also illustrates the variation of mean flame brush thickness δt, as an indication of the transition zone between burned and unburned states of a premixed flame [27]. In this paper, δt is taken to be equal to the distance between mean reduced filtered temperature θ~=0.1 and 0.9, where θ~=(T~-Tr)/(Tp-Tr). By increasing the level of preheating along path I, δt decreases significantly in the flame tip, from about 20 mm to 5 mm. On the other hand, leaning the highly preheated mixture increases the mean flame brush thickness in the flame tip, leading to the thickest flame in case C with δt≃27 mm at the tip, which is even thicker than the benchmark case.Figure6 shows 3D results of mean flame fronts for all cases using the isosurfaces of mean progress variable c~ and corresponding instantaneous flame fronts at c~=0.8, with temperature as indicated. Progress variable is defined as c~=(YF,u-YF)/(YF,u-YF,b) with YF,b and YF,u denoting the mass fraction of fuel species on the burned and unburned sides of the flame, respectively. For both paths I and II, by increasing Tr or ϕ, the flame surface becomes smaller and less wrinkled for the cases close or inside the HiTC zone (cases A3 and B1).Figure 6
Isosurface of mean (left) and instantaneous (right) progress variable at 0.8 colored by mean and instantaneous temperature, respectively, (a) case A1, (b) case A2, (c) case A3, (d) case B1, (e) case B2, and (f) case C.
(a) (b) (c) (d) (e) (f)Figure7 superimposes mean flame front surfaces in a scaled grid for both paths. It can be seen that the maximum mean flame length of case B1 is about one-third of the case A1 and one-quarter of the case C, indicating the same effect of temperature and equivalence ratio of the mixture on the flame length in lean premixed combustion. This illustrates that, by increasing the temperature or enriching the preheated lean premixture, the flame becomes shorter and stabilizes closer to the combustor’s main inlet, which agrees with experimental results of Panoutsos et al. [3]. These phenomena also increase the risk of flashback for the cases in or close to HiTC zones (cases A3 and B1).Figure 7
Superposition of flame surfaces in a scaled grid.
(a)
Path I (b)
Path II


### 5.3. Turbulent Flow Structure

Apart from the flame, flow structure also changes significantly for both paths due to chemistry and turbulence interaction. Figure8 shows streamlines of mean velocity coloured according to its downstream component, z, of mean filtered velocity u~z together with mean flame location. Mean progress variables c~ at 0.2, 0.5, and 0.8 on a vertical plane across the center of the combustor for all cases are shown by solid lines. Time series line charts show the temperature and u~z at 6 different locations or probes marked with the same color at the domain, with indicative dotted lines shown in Figure 8(a).Figure 8
Streamlines of mean velocity and velocity signals: (a) case A1, (b) case A2, (c) case A3, (d) case B1, (e) case B2, and (f) case C.
(a) (b) (c) (d) (e) (f)Figure8(a) illustrates the flow field of the benchmark case. It clearly shows three typical structures of a swirl combustor, namely, the center or inner recirculation zone (IRZ) downstream of the injector generated by vortex breakdown, shear layers originating from the outer edge of the chamber inlet, and corner recirculating zone (CRZ), created in confined configurations due to sudden expansion of the flow at the chamber inlet. As shown in Figures 8(a)–8(d), while the IRZ expands and moves downstream with enhanced preheating, the size of the CRZ increases significantly, with the small vortices located at the bottom corner of this zone growing rapidly into two strong counter rotating vortices for cases A3 and B1. This is supported by the negative close to zero u~z at corresponding locations, as illustrated by the time series line charts. It can be seen by comparing Figures 8(d)–8(f) that leaning the high temperature premixture almost results in reverse flow.Since mass flux of the flow is constant for all cases and gas density is inversely proportional to temperature, the bulk inlet velocity increases for the cases with high preheating level (cases A3 to C). On the other hand, the flames for cases A3 and B1 are more compact and are stabilized close to the inlet with an almost planar shape in the center due to higher laminar flame speed. Therefore, for the cases of A3 and B1, as illustrated byT~ values at locations 2 and 3, reactions take place at the sudden expansion of the combustor, where flow acceleration due to temperature ratio compensates the momentum loss of the expansion in the z-direction leading to modification of CRZ and IRZ as described above.Variation of resolved turbulent kinetic energy of the flow may also be verified by comparing velocity signals of each case across the flame in Figure8. For location 1, at the unburned side of the flame, axial root mean square fluctuation velocity uzrms increases as the inlet temperature increases showing the highest turbulence level at the chamber inlet for the cases in HiTC and flameless mode. This is because in these cases the deceleration of incoming flow with higher velocity happens at almost the same distance as other cases before reaching the IRZ leading to a steeper velocity gradient ∂u~i/∂xj which is a production term in the equation of turbulent kinetic energy. Thus, both resolved and unresolved turbulent kinetic energies of the incoming flow increase with increasing Tr, while decreasing equivalence ratio of highly preheated mixture does not have any significant effect, indicating that ReΔ remains constant for the cases along path II.At locations 2 to 4, similar dependency ofuzrms and T~ occur, although intermittent behavior can be seen. These intermittencies are less obvious in comparison with temperature signals due to other major sources that affect the velocity fluctuation in confined swirling reacting flows, such as vortex shedding and thermos-acoustic instabilities. Focusing on corresponding magnified windows at locations 5 and 6 in the burnt side, low amplitude high frequency fluctuations of u~z occur for the cases of A3 and B1, indicating thermos-acoustic instabilities for these cases.From the above, it can be seen that increasing the preheating temperature of the mixture increases the turbulent kinetic energy of the incoming flow, leading to smaller variations in local equivalence ratio for the cases with more preheating level when the mixture is partially premixed before combustion. This can also be verified more clearly by comparing the profiles of different components of the mean velocity and RMS velocity of all cases in Figure9, where the dump plane and chamber exit are located at Z = 0 and 76 mm, respectively.Figure 9
Profiles of mean velocity components for all cases at different axial positions: axial (blue solid lines), radial (green solid lines), tangential (red solid lines), and RMS (dashed lines).


### 5.4. Turbulence-Flame Interaction

In order to further explain the major trends discussed in the previous subsections, Figure10 shows the location of all cases in a regime diagram for LES and DNS of premixed turbulent combustion proposed by Pitsch [28]. In this regime diagram, η is the Kolmogorov length scale, lG is the Gibson scale which is the smallest scale of the subgrid flame front wrinkling, lm is the broadened flame thickness, δ is the reaction zone thickness and Kaδ is the Karlovitz number based on δ. Since u′ is a fluctuating quantity, ReΔ, DaΔ, and KaΔ are also fluctuating quantities that change in space and time. The markers and error bars show the mode of time-averaged and limits of u′ in the mean flame brush, respectively.Figure 10
Location of all cases on a regime diagram for LES and DNS of premixed turbulent combustion proposed by [28].As discussed in Sections5.1 and 5.3, by moving towards the HiTC zone along path I, u′ and sl increase. However, since the increase of u′ due to volume expansion of the mixture is linear and much less than sl∝Tr2 and sllf≈ν∝T1.7, ReΔ decreases slightly but DaΔ increases significantly along this path bringing the cases A3 and B1 to the wrinkled flamelet region where the flame is less sensitive to turbulence owing to KaΔ≪1 or lf≪η. For path II, ReΔ remains almost constant while DaΔ decreases notably due to significant decrease of sl with ϕ moving the cases B1 and C back to the thin reacting zone where turbulence increases the transport within the chemically inert preheat region owing to KaΔ>1 or lf>η. For case C, the flame is partly resolved because DaΔ>1 and for the other cases the flame is entirely on the subgrid scale.As suggested by Bradley et al. [29] the Karlovitz stretch factor defined as Kf=urmslf/λtsl is used in this paper to explain the evolution of the flame in both paths. With this definition, λt is the Taylor scale and urms is the RMS turbulent velocity. For Kf>1, local quenching by flame stretch becomes less significant [30]. Here, estimating λt≈2Δ, Kf is approximated by Kf≈αDaΔ-1 with α ranging from 9 to 11 in the flame region.For the case with atmospheric temperature (case A1), the flame is stabilized by hot products in both the CRZ and IRZ, as shown in Figure8(a). The flame base is anchored at the boundary of IRZ whereas the flame tip fluctuates between two recirculation zones resulting in high value of δt, particularly in the tip region, and highly wrinkled flame surface as described in Section 5.2. Since the value of Kf is around 2 at the inlet of the combustor, flame fluctuation at the base is limited to the borders of IRZ. By progressively increasing Tr along path I, Kf decreases from 2 to 0.3 in front of the chamber inlet due to considerable growth of DaΔ, as shown in Figure 10. Therefore, as shown in Figures 8(b)–8(d), the flame moves upstream anchoring at the IRZ weakens in case A2 and eventually is lost for the cases A3 and B1 where the flame is stabilized only at the borders of the CRZ. This leads to a remarkable reduction in resolved flame surface, length, and wrinkling for the cases close or inside the HiTC zone as shown in Figures 6 and 7.Along path II, an almost reverse phenomenon occurs since by decreasingϕ from 0.9 to 0.3, Kf increases from 0.3 to 10 in the region between IRZ and CRZ due to the substantial drop in DaΔ. Thus, the flame moves downstream and the role of IRZ becomes increasingly significant for the cases B2 and C. For case C, with highest strain, the flame is mostly stabilized at the boundary of IRZ and furthest from the combustor inlet, as shown in Figure 8(f). This leads to a marked increase in flame surface, length, and wrinkling for the cases close to or inside the flameless zone as shown in Figures 6 and 7.


## 5.1. Laminar Flame Structure

Incorporating detailed chemical mechanisms in LES is still computationally prohibitive due to the complexity of turbulence/combustion interaction. In addition, one additional equation is required to be solved for species in the reaction mechanism using very fine grids [14, 24]. Thus, the present study uses reduced reaction mechanisms, which is common in many simulations of turbulent flames using FRC models. In this paper, the following reduced one step tuned reaction mechanism (OTM) for lean methane/air combustion is used such that(13)CH4+2O2⟶CO2+2H2Owithw˙=9.49×1011e-Ta/TYCH40.7YO20.8,where Ta= 23,650 K is the activation temperature and the species exponents of methane and oxygen are individually tuned matching the GRI 3.0 mechanism.In order to investigate the accuracy of OTM, a series of 1D freely propagating flames with the same mixture condition as described in Table1 is simulated by CANTERA [25] using GRI 3.0 and one step tuned mechanism used in LES of this study. Figures 3(a) and 3(b) compare the velocity profile of 1D freely propagating flame for the cases in paths I and II, respectively, using both GRI and OTM chemistry. GRI results show that increasing the mixture temperature gives rise to the unstrained laminar velocity sl quadratically, whereas decreasing the equivalence ratio of a highly preheated mixture decreases this velocity linearly. As shown in Figure 3(c), temperature profiles of the 1D flame indicate the same dependency of adiabatic flame temperature on temperature and equivalence ratio of the mixture but both in a linear fashion. These profiles also show that the laminar flame thickness calculated based on lf=ν/Prsl [24] decreases slightly as the preheating level rises and increases significantly when the mixture becomes leaner for the cases in path II. These tendencies may be explained by sl∝μw˙/ρ2 and lf∝μ/ρsl [26]. Combining the variation of laminar flame speed and thickness by calculation of chemical time scale according to τc=lf/sl reveals that chemistry accelerates by increasing the temperature and decelerates by reducing the equivalence ratio in the lean mixture.Figure 3
(a), (b) Velocity, (c) temperature, and (d)NOx mole fraction profiles of 1D freely propagating flame.
(a) (b) (c) (d)Mole fraction profiles ofNO species for all cases are shown in Figure 3(d). It shows that the amount of NOx emission in the product is strongly dependent on product temperature and is important when Tp is more than 1900 K. This illustrates that the cases inside or near the flameless combustion mode, cases B2 and C, are almost NOx-free. Detailed results of the 1D laminar simulation are also presented in Table 2.Table 2
Detailed results of 1D laminar simulation.
Case code OTM GRI 3.0 l f (mm) s l (m/s) τ c (ms) T p (K) l f (mm) s l (m/s) τ c (ms) T p (K) A1 0.060 0.366 0.164 2165 0.067 0.337 0.199 2137 A2 0.055 0.991 0.055 2316 0.064 0.843 0.076 2255 A3 0.043 2.225 0.019 2472 0.056 1.691 0.033 2364 B1 0.032 4.483 0.007 2638 0.047 3.035 0.015 2465 B2 0.094 1.515 0.062 1978 0.087 1.640 0.053 1973 C 0.297 0.482 0.617 1588 0.245 0.583 0.419 1589Comparing the results of two different mechanisms indicates that although the one step tuned mechanism used in turbulent simulation of this study overpredicts the adiabatic temperature and laminar flame velocity and underpredicts the laminar flame thickness of highly preheated cases with equivalence ratio of 0.9 (cases A3 and B1), it predicts adequately the trends of laminar flame velocity, flame thickness, and adiabatic temperature in both paths. It is worth mentioning that the self-ignition temperature in each case is calculated using zero-dimension well stirred reactor simulation for a residence time of 1 second and atmospheric pressure, as suggested by [8], with the same reduced chemistry mechanism exploited in the present LES simulation.


## 5.2. Temperature Field and Turbulent Flame Structure

Figure4 shows the profiles of characteristic temperatures extracted from the 3D domain of each case in two different paths. It illustrates that increasing the preheat mixture temperature increases the maximum temperature of the product Tp in the combustor linearly, whereas decreasing the equivalence ratio has a reverse effect. Increase of product temperature along path I is slightly less than the increase in reactant temperature between corresponding case leading to moderate decline in the maximum temperature rise ΔT, but the temperature drop along path II is significant. It can be seen that the maximum temperature of case B1 is about 500 K higher than case A1 with the preheating level of 600 K while decreasing equivalence ratio from 0.9 to 0.5 results in a drop of 700 K in ΔT. This variation of ΔT leads to a large decrease in the temperature ratio Tp/Tr from 7 to 3 along path I and from 3 to 1.7 along path II, playing an important role in flame-turbulence interaction. Since wall temperature of the combustor is kept constant for all the simulations, increase in exit temperature of the combustor for the cases along path I is attributed to higher heat loss through the walls. Consequently, only increasing preheating temperature of the mixture enhances the overall heat loss in the combustor; therefore, exit temperature of the product does not change significantly from case A1 to B1.Figure 4
Characteristic temperatures.
(a)
Path I (b)
Path IIFigure5 shows the mean filtered temperature contours T~ on a vertical plane across the center of the combustor for all cases. The symbol f denotes time-averaged of variable f. As shown in Figures 5(a)–5(d), the zone with temperature higher than 1800 K expands along path I, increasing the area for formation of thermal NOx in the combustor. In contrast, comparing the isothermal lines of 1800 K in Figures 5(d)–5(f) reveals that leaning the preheated mixture along path II decreases the heat loss through the combustor walls as well as decreases the regions for the formation of thermal NOx. It is worth mentioning that without preheating the mixture, the minimum equivalence ratio that can be achieved with the current burner is around 0.7. Thus, preheating is necessary for expanding the lean flammability limits in order to approach the flameless zone.Figure 5
Mean temperature contours across the center of the combustor.
(a)
Without preheat (b)
200°K preheat (c)
400°K preheat (d)
600°K preheat (e)
600°K preheatϕ=0.5 (f)
600°K preheatϕ=0.3Figure5 also illustrates the variation of mean flame brush thickness δt, as an indication of the transition zone between burned and unburned states of a premixed flame [27]. In this paper, δt is taken to be equal to the distance between mean reduced filtered temperature θ~=0.1 and 0.9, where θ~=(T~-Tr)/(Tp-Tr). By increasing the level of preheating along path I, δt decreases significantly in the flame tip, from about 20 mm to 5 mm. On the other hand, leaning the highly preheated mixture increases the mean flame brush thickness in the flame tip, leading to the thickest flame in case C with δt≃27 mm at the tip, which is even thicker than the benchmark case.Figure6 shows 3D results of mean flame fronts for all cases using the isosurfaces of mean progress variable c~ and corresponding instantaneous flame fronts at c~=0.8, with temperature as indicated. Progress variable is defined as c~=(YF,u-YF)/(YF,u-YF,b) with YF,b and YF,u denoting the mass fraction of fuel species on the burned and unburned sides of the flame, respectively. For both paths I and II, by increasing Tr or ϕ, the flame surface becomes smaller and less wrinkled for the cases close or inside the HiTC zone (cases A3 and B1).Figure 6
Isosurface of mean (left) and instantaneous (right) progress variable at 0.8 colored by mean and instantaneous temperature, respectively, (a) case A1, (b) case A2, (c) case A3, (d) case B1, (e) case B2, and (f) case C.
(a) (b) (c) (d) (e) (f)Figure7 superimposes mean flame front surfaces in a scaled grid for both paths. It can be seen that the maximum mean flame length of case B1 is about one-third of the case A1 and one-quarter of the case C, indicating the same effect of temperature and equivalence ratio of the mixture on the flame length in lean premixed combustion. This illustrates that, by increasing the temperature or enriching the preheated lean premixture, the flame becomes shorter and stabilizes closer to the combustor’s main inlet, which agrees with experimental results of Panoutsos et al. [3]. These phenomena also increase the risk of flashback for the cases in or close to HiTC zones (cases A3 and B1).Figure 7
Superposition of flame surfaces in a scaled grid.
(a)
Path I (b)
Path II


## 5.3. Turbulent Flow Structure

Apart from the flame, flow structure also changes significantly for both paths due to chemistry and turbulence interaction. Figure8 shows streamlines of mean velocity coloured according to its downstream component, z, of mean filtered velocity u~z together with mean flame location. Mean progress variables c~ at 0.2, 0.5, and 0.8 on a vertical plane across the center of the combustor for all cases are shown by solid lines. Time series line charts show the temperature and u~z at 6 different locations or probes marked with the same color at the domain, with indicative dotted lines shown in Figure 8(a).Figure 8
Streamlines of mean velocity and velocity signals: (a) case A1, (b) case A2, (c) case A3, (d) case B1, (e) case B2, and (f) case C.
(a) (b) (c) (d) (e) (f)Figure8(a) illustrates the flow field of the benchmark case. It clearly shows three typical structures of a swirl combustor, namely, the center or inner recirculation zone (IRZ) downstream of the injector generated by vortex breakdown, shear layers originating from the outer edge of the chamber inlet, and corner recirculating zone (CRZ), created in confined configurations due to sudden expansion of the flow at the chamber inlet. As shown in Figures 8(a)–8(d), while the IRZ expands and moves downstream with enhanced preheating, the size of the CRZ increases significantly, with the small vortices located at the bottom corner of this zone growing rapidly into two strong counter rotating vortices for cases A3 and B1. This is supported by the negative close to zero u~z at corresponding locations, as illustrated by the time series line charts. It can be seen by comparing Figures 8(d)–8(f) that leaning the high temperature premixture almost results in reverse flow.Since mass flux of the flow is constant for all cases and gas density is inversely proportional to temperature, the bulk inlet velocity increases for the cases with high preheating level (cases A3 to C). On the other hand, the flames for cases A3 and B1 are more compact and are stabilized close to the inlet with an almost planar shape in the center due to higher laminar flame speed. Therefore, for the cases of A3 and B1, as illustrated byT~ values at locations 2 and 3, reactions take place at the sudden expansion of the combustor, where flow acceleration due to temperature ratio compensates the momentum loss of the expansion in the z-direction leading to modification of CRZ and IRZ as described above.Variation of resolved turbulent kinetic energy of the flow may also be verified by comparing velocity signals of each case across the flame in Figure8. For location 1, at the unburned side of the flame, axial root mean square fluctuation velocity uzrms increases as the inlet temperature increases showing the highest turbulence level at the chamber inlet for the cases in HiTC and flameless mode. This is because in these cases the deceleration of incoming flow with higher velocity happens at almost the same distance as other cases before reaching the IRZ leading to a steeper velocity gradient ∂u~i/∂xj which is a production term in the equation of turbulent kinetic energy. Thus, both resolved and unresolved turbulent kinetic energies of the incoming flow increase with increasing Tr, while decreasing equivalence ratio of highly preheated mixture does not have any significant effect, indicating that ReΔ remains constant for the cases along path II.At locations 2 to 4, similar dependency ofuzrms and T~ occur, although intermittent behavior can be seen. These intermittencies are less obvious in comparison with temperature signals due to other major sources that affect the velocity fluctuation in confined swirling reacting flows, such as vortex shedding and thermos-acoustic instabilities. Focusing on corresponding magnified windows at locations 5 and 6 in the burnt side, low amplitude high frequency fluctuations of u~z occur for the cases of A3 and B1, indicating thermos-acoustic instabilities for these cases.From the above, it can be seen that increasing the preheating temperature of the mixture increases the turbulent kinetic energy of the incoming flow, leading to smaller variations in local equivalence ratio for the cases with more preheating level when the mixture is partially premixed before combustion. This can also be verified more clearly by comparing the profiles of different components of the mean velocity and RMS velocity of all cases in Figure9, where the dump plane and chamber exit are located at Z = 0 and 76 mm, respectively.Figure 9
Profiles of mean velocity components for all cases at different axial positions: axial (blue solid lines), radial (green solid lines), tangential (red solid lines), and RMS (dashed lines).


## 5.4. Turbulence-Flame Interaction

In order to further explain the major trends discussed in the previous subsections, Figure10 shows the location of all cases in a regime diagram for LES and DNS of premixed turbulent combustion proposed by Pitsch [28]. In this regime diagram, η is the Kolmogorov length scale, lG is the Gibson scale which is the smallest scale of the subgrid flame front wrinkling, lm is the broadened flame thickness, δ is the reaction zone thickness and Kaδ is the Karlovitz number based on δ. Since u′ is a fluctuating quantity, ReΔ, DaΔ, and KaΔ are also fluctuating quantities that change in space and time. The markers and error bars show the mode of time-averaged and limits of u′ in the mean flame brush, respectively.Figure 10
Location of all cases on a regime diagram for LES and DNS of premixed turbulent combustion proposed by [28].As discussed in Sections5.1 and 5.3, by moving towards the HiTC zone along path I, u′ and sl increase. However, since the increase of u′ due to volume expansion of the mixture is linear and much less than sl∝Tr2 and sllf≈ν∝T1.7, ReΔ decreases slightly but DaΔ increases significantly along this path bringing the cases A3 and B1 to the wrinkled flamelet region where the flame is less sensitive to turbulence owing to KaΔ≪1 or lf≪η. For path II, ReΔ remains almost constant while DaΔ decreases notably due to significant decrease of sl with ϕ moving the cases B1 and C back to the thin reacting zone where turbulence increases the transport within the chemically inert preheat region owing to KaΔ>1 or lf>η. For case C, the flame is partly resolved because DaΔ>1 and for the other cases the flame is entirely on the subgrid scale.As suggested by Bradley et al. [29] the Karlovitz stretch factor defined as Kf=urmslf/λtsl is used in this paper to explain the evolution of the flame in both paths. With this definition, λt is the Taylor scale and urms is the RMS turbulent velocity. For Kf>1, local quenching by flame stretch becomes less significant [30]. Here, estimating λt≈2Δ, Kf is approximated by Kf≈αDaΔ-1 with α ranging from 9 to 11 in the flame region.For the case with atmospheric temperature (case A1), the flame is stabilized by hot products in both the CRZ and IRZ, as shown in Figure8(a). The flame base is anchored at the boundary of IRZ whereas the flame tip fluctuates between two recirculation zones resulting in high value of δt, particularly in the tip region, and highly wrinkled flame surface as described in Section 5.2. Since the value of Kf is around 2 at the inlet of the combustor, flame fluctuation at the base is limited to the borders of IRZ. By progressively increasing Tr along path I, Kf decreases from 2 to 0.3 in front of the chamber inlet due to considerable growth of DaΔ, as shown in Figure 10. Therefore, as shown in Figures 8(b)–8(d), the flame moves upstream anchoring at the IRZ weakens in case A2 and eventually is lost for the cases A3 and B1 where the flame is stabilized only at the borders of the CRZ. This leads to a remarkable reduction in resolved flame surface, length, and wrinkling for the cases close or inside the HiTC zone as shown in Figures 6 and 7.Along path II, an almost reverse phenomenon occurs since by decreasingϕ from 0.9 to 0.3, Kf increases from 0.3 to 10 in the region between IRZ and CRZ due to the substantial drop in DaΔ. Thus, the flame moves downstream and the role of IRZ becomes increasingly significant for the cases B2 and C. For case C, with highest strain, the flame is mostly stabilized at the boundary of IRZ and furthest from the combustor inlet, as shown in Figure 8(f). This leads to a marked increase in flame surface, length, and wrinkling for the cases close to or inside the flameless zone as shown in Figures 6 and 7.


## 6. Conclusions

In this paper, flame and flow structures of a swirl-stabilized pilot combustor are investigated in conventional, high HiTC, and flameless mode using large eddy simulation with a finite rate chemistry combustion model with single step tuned mechanism. Results show that moving towards HiTC mode solely by increasing the preheating level, the combustor is likely to have more thermal NOx due to higher temperature of the product. It also has a higher risk of flashback due to lower value of Karlovitz stretch rate. The flame becomes shorter and thinner due to modification of laminar flame speed and higher turbulent kinetic energy. On the other hand, approaching the flameless mode along path II, leaning the highly preheated mixture leads to an almost thermal NOx-free combustion with lower risk of flashback as well as thicker and longer flames without any significant modification of upstream turbulence. This indicates that designing a combustion system inside or close to the flameless mode could be an optimized solution for exploiting the mixture preheating method in lean premixed combustion. With qualitative agreements with experiments, this paper shows the capability of the LES PaSR model with one step tuned mechanism in capturing the main features of the turbulent flame in a wide range of mixture temperatures and equivalence ratios, although further development of the current method is needed to demonstrate its capability with complex chemistry, especially for prediction of NOx emission, which is the advantage of the PaSR model.



---

*Source: 8261560-2016-05-31.xml*