# Large Eddy Simulation of Premixed Stratified Swirling Flame Using the Finite Rate Chemistry Approach

Authors: Yinli Xiao; Zhengxin Lai; Wenyan Song
Journal: International Journal of Aerospace Engineering
 (2019)

---

## Abstract

Large eddy simulations of a stratified swirling flow of a Cambridge swirl burner for both nonreacting and reacting cases are conducted using a finite rate chemistry approach represented by a partially stirred reactor model. The large eddy simulation predictions are compared with experimental measurements for velocity, temperature, and concentrations of major species. The agreement is found in overall trend of velocity prediction, but temperature and concentration of major species show slight discrepancies in the central region. Two reduced chemical mechanisms are examined in the present paper with the objective of assessing their capabilities in predicting swirling flame characteristics, and the distinct difference using two mechanisms is found in CO distribution profiles, which is considered the consequence of different kinetics of CO-CO2 equilibrium. Flow structures are qualitatively and quantitatively analyzed with numerical results. Large-scale vortex structures and precession motions are observed in both nonreacting and reacting cases. Frequency of vortex shedding is identified from the point data of instantaneous velocity in the discharging stream-induced shear layer. On this basis, the intensity and frequency of precession motion are shown to be enhanced in the presence of combustion. Large-scale wrinkling of the flame surface is resolved and characterized in the flame zone, and the effect of mixture stratification is then further discussed.


---

## Body

## 1. Introduction

Lean premixed combustion has been widely adopted in practical combustion systems for it can offer low NOx emission by decreasing peak reacting temperature. In practical combustion systems, however, constrained by physical space and fuel-air mixing time, inhomogeneous mixture and spatial gradient of mixture equivalence ratio often exhibit. This may lead to combustion taking place under stratified condition. Study shows that stratified mixture with nonuniform fuel concentration has wider flame stability and ignition range compared with homogeneous lean premixed condition [1], but complete understanding of stratified combustion is still out of reach till now. Complex turbulent flow and chemistry reactions make flow field measurements in practical combustion chamber a tough challenge. To investigate the stratification and swirl effects on turbulent combustion, the Cambridge swirl burner (SwB) was designed by Sweeney et al. [2–4] and detailed measurements of the flow fields were carried out by Zhou et al. [5]. The experimental measurements are valuable for combustion model assessment and numerical simulation validation.Large eddy simulation (LES) is under intense development in the last two decades and has become a promising method to study unsteady characteristics of flow and combustion because it allows for higher fidelity than RANS, at a lower computational cost than DNS [6–8]. The philosophy behind LES is to decompose turbulent motions into a grid-scale portion and a sub-grid-scale (SGS) portion. The grid-scale portion is solved directly whilst the sub-grid-scale portion needs to be modeled using the subgrid closure models [9]. LES has been proven to be competent to capture complex unsteady turbulent flow characteristics, e.g., large-scale vortex structure generation and evolution [10–12]. Nambully et al. [13] conducted LES of two nonswirling cases of the SwB associated with a filtered-laminar-flame model in combination with a presumed PDF SGS closure. LES well reproduced the equivalence ratio distribution in the vicinity of a bluff body recirculation zone, and the differential diffusion effect of mixture in stratified combustion was analyzed based on numerical results. Proch and Kempf [14] simulated three stratified nonswirling cases using the artificial thickened flame (ATF) approach in combination with flamelet-generated manifold (FGM) lookup tables, and the process of stratified combustion was characterized and investigated from simulation results. Brauner et al. [15] simulated six SwB cases with varying degrees of swirl and mixture stratification by using the LES-PDF method with the Eulerian stochastic field solution method. Velocity and major species profiles showed good agreement with experimental data, whilst local temperature profiles had slight discrepancies in the central regions for swirling cases. Zhang et al. [16] simulated both reacting and nonreacting SwB cases to assess different LES SGS models. Dynamic thickened flame combined with the FGM tabulation approach (DTF-FGM) and presumed PDF combined with the FGM approach (PPDF-FGM) were used to model premixed flame and stratified flame, where flame features and unsteady large-scale structure behaviors were investigated.Generally, combustion models are classified into two main categories: flamelet-like models [17] and finite rate chemistry models. Flamelet-like models regard turbulent flame as an ensemble of stretched laminar flamelets. Flamelets are very thin because the reactions that determine the fuel consumption rate are very fast and chemical time scales are very short. As a result, chemistry is most active within a thin layer, and flamelet structure is mainly affected by local turbulent eddies. Since the flamelet concept focuses on the location of the flame surface, which is defined as an isosurface of a nonreacting scalar, where the species composition and temperature are determined through a manifold library based on the flame coordinate, the reactive scalars are decoupled from actual chemistry reactions [18, 19]. In most situations, the flamelet-like models are computationally efficient. However, in stratified or partially premixed conditions such as the SwB configuration, where the local equivalence ratio distributes with considerable spatial gradient, local burning velocity and flame thickness vary. The flame regime could also vary throughout the flow field from wrinkled flamelets up to thin reaction zones. Actually, a sub-grid-scale combustion closure model that could be able to dynamically follow the change in the flame regime is necessary to model the entire flame region [13]. For the finite rate chemistry models, on the contrary, no assumptions are applied to a flow or flame, but the species equations and flame front propagation based on chemistry reactions and Arrhenius rates were directly solved. This feature makes it possible to take quenching/reignition into account dynamically and follow the change in the local flame regime. Nogenmyr et al. [20, 21] applied both the PaSR model and level-set G-equation flamelet model to simulate a stratified/partially premixed low swirl methane/air flame. Their research showed inhomogeneous equivalence ratio distribution may lead to local extinction in the shear layer, which is difficult to simulate. But since the finite rate chemistry method does not impose any restriction of the combustion regime, the PaSR model well captured the dynamic characteristics of flame such as local extinction and reignition phenomenon and reproduced the feature of stratified combustion.Actually, it is difficult to handle turbulent combustion due to the complicated physical process behind flow and combustion as well as the fact that there is a lack of generally acceptable and universally valid models. In the present work, Cambridge swirl flame SwB7 is investigated using the PaSR model with reduced chemical mechanisms. The objective is to evaluate the performance of the PaSR model and the predictive capabilities of reduced chemical mechanisms in LES of reacting flow and investigate the basic unsteady characteristics of the swirling flame. All the simulations are performed based on the open source CFD toolbox OpenFOAM [22].


## 2. Configuration and Case Description

The Cambridge swirl burner, referred to as SwB, is schematically presented in Figure1. The burner consists of two coannular tubes. A ceramic central bluff body is used to stabilize flame. Premixed methane/air mixtures flow through the coannular tubes, and coflow air is supplied around the outer annulus to prevent the entrainment of ambient air. In experimental measurements, the bulk velocity of the inner annulus was set at Uj=8.31m/s and the axial velocity in the outer annulus was Us=18.7m/s, which was more than twice the value of the velocity in the inner annulus in order to generate substantial levels of shear between the two flows. Coflow air was supplied with a bulk velocity Ue=0.4m/s. Reynolds numbers derived from bulk velocities and burner geometry are Rei=5960 and Reo=11500 for the inner and outer annulus, respectively. All flows were with room temperature of 295 K and atmospheric pressure. The swirl flow ratio, the ratio of the outer flow rate through the swirl plenum to the total outer annulus flow, was set at 0.33, which generated a corresponding geometric swirl number S=0.79. The swirl number was calculated as the ratio of measured mean tangential velocity to axial velocity S=Ws/Us. The equivalence ratio of the inner annulus, ϕi, and outer annulus, ϕo, are 1.0 and 0.5, respectively, indicating the stratification ratio sR=2.Figure 1
Schematic of the Cambridge swirl burner [16].


## 3. Computational Methods

### 3.1. LES SGS Closure Model

In a large eddy simulation, the grid-scale motions are solved directly based on the spatial filtered governing equations by governing equations whilst the sub-grid-scale motions, which contain the influence on grid-scale motions, are taken into account with closure models. In the present work, the sub-grid-scale modeling is achieved by a dynamic one-equation eddy viscosity model, which solves a transport equation for the sub-grid-scale kinetic energy rather than calculating the sub-grid-scale kinetic energy from known resolved quantities. Due to no assumption of local balance between the sub-grid-scale energy production and dissipation rate, the one-equation model is expected to be more appropriate than the algebraic model such as the Smagorinsky model in the regions where local balance is violated [23]. The transport equation for sub-grid-scale kinetic energy is formulated as follows:
(1)∂ksgs∂t+ui¯∂ksgs∂xi=−τij∂ui¯∂xj−Cϵksgs3/2Δ+∂∂xiνkσk∂ksgs∂xi,where ksgs is the subgrid kinetic energy, σk=1 is usually adopted, and subgrid stresses τij are described by
(2)τij=−2μsgsSij~+23ksgsδij,where Sij is the filtered strain rate tensor and μsgs is the SGS eddy viscosity, which is calculated as follows:
(3)μsgs=ρ¯CkΔksgs,where Δ=ΔxΔyΔz1/3 is the filtered width. The values of model coefficients Cϵ and Ck are obtained by the dynamical method described by Germano et al. [24]


### 3.2. Turbulent Combustion Model

From the mathematical point of view, reactive flow equations are the balance equations of mass, momentum, and energy, describing convection, diffusion, and reactions [25, 26]. LES solves the low-pass filtered equations on the grids of manageable size. The general Favre filtered three-dimensional Navier-Stokes equations are given as follows [27]:
(4)∂ρ¯∂t+∂∂xiρ¯u~i=0,(5)∂ρ¯u~j∂t+∂∂xiρ¯u~iu~j=−∂ρ¯∂xj+∂∂xiτ¯ij−τijsgs,(6)∂ρ¯Y~i∂t+∂∂xiρ¯u~iY~i=∂∂xiμsc∂Y~i∂xi−ρ¯uiYi¯−u~iY~i,(7)∂ρ¯h~s∂t+∂∂xiρ¯u~ih~s=∂P¯∂t+u~i∂P¯∂xi+∂∂xiμPr∂h~s∂xi−ρ¯uihs¯−u~ih~s+τij∂¯ui∂xj+ω̇T¯,where ρ is the density, P is the pressure, Yi is the mass fraction of the ith species, hs is the sensible enthalpy, sc is the Schmidt number, and Pr is the Prandtl number. The heat release ω̇T¯ due to combustion is modeled by heat of formation ∆hf,i0 and reaction rate ω̇i¯ as follows:
(8)ω̇T¯=∑i=1NΔhf,i0ω̇i¯.With the understanding of turbulence-combustion interaction, the different finite rate chemistry models are developed to model the low-pass filtered reaction rates. The PaSR model has been widely tested and proven to have a good performance in turbulent-combustion simulation [28, 29]. The development of the PaSR model is based on the assumption that the combustion takes place in the fine structures where the turbulence kinetic energy dissipates. These well-mixed regions are not evenly distributed in time and space but concentrate in particular regions that occupy a certain fraction of the flow [30]. The PaSR model incorporates the sequential processes of micromixing and chemical reactions. Microscale processes responsible for turbulence energy dissipation and molecular mixing are highly intermittent and just occupy a small fraction of fluid volume, but have a significant impact on chemical reactions. The characteristic dimensions of fine structures are small compared with the LES filter width. Thus, for a computational cell, it can be considered as a partially stirred reactor which consists of two parts, one for well-mixed reacting volume (R) and the other for nonreacting volume (N). Nonreacting volume and reacting volume are assumed as locally uniform elements, and the interaction of all elements occurs by “exchange with the mean” [31]. The cell mean value such as species concentration and chemical reaction rates can be calculated as follows [31]:
(9)Yi=1−κYiN+κYiR,(10)ω̇i¯=κω̇i,where YiN and YiR denote the species properties in nonreacting and reacting volume, respectively, and Yi is the cell mean value. Similarly, ω̇i and ω̇i¯ denote the change rate of species concentration in reacting volume and the cell mean. Microscale processes and chemical reactions and their interactions are incorporated by the reacting volume fraction κ, which can be described as follows:
(11)κ=τchτm+τch.The chemical reaction time scaleτch, which is determined by chemical kinetics, is estimated as follows:
(12)1τch=max−ω̇fuelYfuel,−ω̇o2Yo2,and the micromixing time scale τm, which is determined by turbulent microprocesses, is estimated as follows:
(13)τm=Cmixνeffϵ.The effective viscosityνeff is the sum of laminar viscosity ν and SGS viscosity νsgs, and the model coefficient Cmix, which relates to turbulent mixing time scale, needs to be estimated a priori by
(14)Cmix=11+CμRet,where Cμ=0.09 and Ret=κ2/ν⋅ε denotes the turbulence Reynolds number.


### 3.3. Chemical Kinetics

The finite rate chemistry approach has been proven to be capable of handling turbulent combustion. However, a drawback with the finite rate method is that it requires high spatial resolution and it is sensitive to the selection of the reaction mechanism. It is necessary for the finite rate chemistry models to solve transport equations for each species. As a consequence, using a detailed mechanism can drastically increase the computational cost because a detailed mechanism usually contains tens of species and hundreds of intermediate reactions. The expensive computational cost caused by a detailed mechanism makes the choice of reaction mechanism a trade-off between computational capability and accuracy; also, it spurs the use of reduced mechanisms or an adaptable flamelet assumption approach [32–34] in complex situations to capture the flame characteristics. Several studies have been conducted to assess the predictive capabilities of the skeletal mechanisms and global mechanisms. Benim et al. [35] investigated a swirling flame in a model combustor by using the EDC model with a 2-step global mechanism. It is found that global mechanism provided a good overall accuracy compared with the experimental data. Fedina et al. [36] assessed several skeletal and global mechanisms in simulating premixed combustion in the SGT-100 dry low emission burner. The studies suggested that a reaction mechanism plays a significant role in numerical combustion, and the reduced reaction mechanisms are feasible in LES.In the present work, the 2-step global mechanism for methane oxidation is selected to handle the combustion and the chemical reactions are formulated as equations (15) and (16). Two sets of Arrhenius parameters are adopted to describe the chemical kinetics, one for the 2sCM2 scheme [37] (referred to as Mech I) and the other for the Arrhenius parameter-optimized 2-step scheme (referred to as Mech II) reduced based on the GRI 3.0.
(15)CH4+32o2→Co+2H2o(16)Co+12o2↔Co2The first reaction (15) is irreversible, the reaction rate is given by equation (17), and the kinetic parameters are reported in Table 1. The second reaction (16) is reversible and the reaction rates described by the two schemes are quite different, which were listed as equation (18) for Mech I and equation (19) for Mech II, where Ke is the equilibrium constant whereas kinetic parameters are reported in Table 2.
(17)ω̇1=A1Tβ1ρYCH4WCH4α1CH4ρYo2Wo2α1o2exp−Ea1RT,(18)ω̇2=A2fTβ2ρYCoWCoα2CoρYo2Wo2α2o2−1KeρYCo2WCo2α2Co2exp−Ea2RT,(19)ω̇2=A2fTβ2fρYCoWCoα2CoρYo2Wo2α2o2exp−Ea2fRT−A2rTβ2rρYCo2WCo2α2Co2exp−Ea2rRT.Table 1
Arrhenius parameters for the first reaction (units: kmol, m, s, kJ, and K).
A1
β1
Ea1
α
Mech I
2.0×1012
0
1.466×105
αCH4=0.9, αo2=1.1
Mech II
4.56×109
-0.0711
1.184×105
αCH4=0.5819, αo2=0.9223Table 2
Arrhenius parameters for the second reaction (units: kmol, m, s, kJ, and K).
A2
β2
Ea2
α
Mech I
6.32×107
0
51000
αCo=1.0, αo2=0.5, αCo2=1.0
Mech II
3.78×1011
0.221
77966
αCo=1.822, αo2=0.928
6.07×105
-0.108
132620
αCo2=1.073


### 3.4. Numerical Setup and Boundary Conditions

The three-dimensional computational domain with a diameter of 200 mm for the current simulation extends from 20 mm upstream of the bluff body to 300 mm downstream of the burner exit. The computational domain is discretized using nonuniform structured mesh. A section view of mesh is shown in Figure2. The mesh is refined near the bluff body and shear layers, which generates a corresponding minimal mesh space of 0.4 mm and overall cell number of about 1.6 million. The minimum grid spacing in the vicinity of the burner exit is considered a reasonable size to resolve the turbulence in the reacting zone [15]. The OpenFOAM reactingFoam solver is employed in the large eddy simulation of reaction cases, which is an unsteady reactive flow solver, and the thermophysical model for reacting mixture is based on compressibility ψ=RT−1. The transport property and thermodynamic property such as μ and Cp are calculated using Sutherland’s law and from JANAF tables of thermodynamics, respectively, and the energy variable used in the solution is based on sensible enthalpy hs.Figure 2
Grid resolution of the computational domain.At the inlets, Dirichlet conditions are adopted for all variables except for pressure. For pressure, zero Neumann condition is used. At the outlet, Dirichlet condition is used for pressure and zero Neumann conditions are adopted for other variables. The experimental data for velocity fluctuation is available only at the near exit, e.g., at 2 mm or 5 mm downstream of the burner exit [4, 5]. Uniform velocity boundary condition is adopted at the inlet, and turbulent intensity is controlled by κ and ε, which are estimated based on the near-exit velocity fluctuation data. All the walls are treated as no slip conditions with wall functions, and the lateral boundary of coflow is the freestream condition. Pressure-velocity coupling is performed using the PIMPLE algorithm, which is a combination of SIMPLE and PISO algorithms. The diffusion terms are discretized through the second-order central differencing scheme, and the convection terms are discretized through the second-order limited central differencing scheme. Temporal discretization is the first-order implicit Euler scheme with an adjustable time step to ensure that the maximal Courant number is below 0.5. The LES was initialized by three-dimensional steady-state RANS results. After the flow field fully developed, statistical sampling has been carried out for at least 5τ. τ is a characteristic domain flow through time and is defined as τ=L/Uj, where L is the streamwise domain length and Uj is the inner annulus jet velocity.


## 3.1. LES SGS Closure Model

In a large eddy simulation, the grid-scale motions are solved directly based on the spatial filtered governing equations by governing equations whilst the sub-grid-scale motions, which contain the influence on grid-scale motions, are taken into account with closure models. In the present work, the sub-grid-scale modeling is achieved by a dynamic one-equation eddy viscosity model, which solves a transport equation for the sub-grid-scale kinetic energy rather than calculating the sub-grid-scale kinetic energy from known resolved quantities. Due to no assumption of local balance between the sub-grid-scale energy production and dissipation rate, the one-equation model is expected to be more appropriate than the algebraic model such as the Smagorinsky model in the regions where local balance is violated [23]. The transport equation for sub-grid-scale kinetic energy is formulated as follows:
(1)∂ksgs∂t+ui¯∂ksgs∂xi=−τij∂ui¯∂xj−Cϵksgs3/2Δ+∂∂xiνkσk∂ksgs∂xi,where ksgs is the subgrid kinetic energy, σk=1 is usually adopted, and subgrid stresses τij are described by
(2)τij=−2μsgsSij~+23ksgsδij,where Sij is the filtered strain rate tensor and μsgs is the SGS eddy viscosity, which is calculated as follows:
(3)μsgs=ρ¯CkΔksgs,where Δ=ΔxΔyΔz1/3 is the filtered width. The values of model coefficients Cϵ and Ck are obtained by the dynamical method described by Germano et al. [24]


## 3.2. Turbulent Combustion Model

From the mathematical point of view, reactive flow equations are the balance equations of mass, momentum, and energy, describing convection, diffusion, and reactions [25, 26]. LES solves the low-pass filtered equations on the grids of manageable size. The general Favre filtered three-dimensional Navier-Stokes equations are given as follows [27]:
(4)∂ρ¯∂t+∂∂xiρ¯u~i=0,(5)∂ρ¯u~j∂t+∂∂xiρ¯u~iu~j=−∂ρ¯∂xj+∂∂xiτ¯ij−τijsgs,(6)∂ρ¯Y~i∂t+∂∂xiρ¯u~iY~i=∂∂xiμsc∂Y~i∂xi−ρ¯uiYi¯−u~iY~i,(7)∂ρ¯h~s∂t+∂∂xiρ¯u~ih~s=∂P¯∂t+u~i∂P¯∂xi+∂∂xiμPr∂h~s∂xi−ρ¯uihs¯−u~ih~s+τij∂¯ui∂xj+ω̇T¯,where ρ is the density, P is the pressure, Yi is the mass fraction of the ith species, hs is the sensible enthalpy, sc is the Schmidt number, and Pr is the Prandtl number. The heat release ω̇T¯ due to combustion is modeled by heat of formation ∆hf,i0 and reaction rate ω̇i¯ as follows:
(8)ω̇T¯=∑i=1NΔhf,i0ω̇i¯.With the understanding of turbulence-combustion interaction, the different finite rate chemistry models are developed to model the low-pass filtered reaction rates. The PaSR model has been widely tested and proven to have a good performance in turbulent-combustion simulation [28, 29]. The development of the PaSR model is based on the assumption that the combustion takes place in the fine structures where the turbulence kinetic energy dissipates. These well-mixed regions are not evenly distributed in time and space but concentrate in particular regions that occupy a certain fraction of the flow [30]. The PaSR model incorporates the sequential processes of micromixing and chemical reactions. Microscale processes responsible for turbulence energy dissipation and molecular mixing are highly intermittent and just occupy a small fraction of fluid volume, but have a significant impact on chemical reactions. The characteristic dimensions of fine structures are small compared with the LES filter width. Thus, for a computational cell, it can be considered as a partially stirred reactor which consists of two parts, one for well-mixed reacting volume (R) and the other for nonreacting volume (N). Nonreacting volume and reacting volume are assumed as locally uniform elements, and the interaction of all elements occurs by “exchange with the mean” [31]. The cell mean value such as species concentration and chemical reaction rates can be calculated as follows [31]:
(9)Yi=1−κYiN+κYiR,(10)ω̇i¯=κω̇i,where YiN and YiR denote the species properties in nonreacting and reacting volume, respectively, and Yi is the cell mean value. Similarly, ω̇i and ω̇i¯ denote the change rate of species concentration in reacting volume and the cell mean. Microscale processes and chemical reactions and their interactions are incorporated by the reacting volume fraction κ, which can be described as follows:
(11)κ=τchτm+τch.The chemical reaction time scaleτch, which is determined by chemical kinetics, is estimated as follows:
(12)1τch=max−ω̇fuelYfuel,−ω̇o2Yo2,and the micromixing time scale τm, which is determined by turbulent microprocesses, is estimated as follows:
(13)τm=Cmixνeffϵ.The effective viscosityνeff is the sum of laminar viscosity ν and SGS viscosity νsgs, and the model coefficient Cmix, which relates to turbulent mixing time scale, needs to be estimated a priori by
(14)Cmix=11+CμRet,where Cμ=0.09 and Ret=κ2/ν⋅ε denotes the turbulence Reynolds number.


## 3.3. Chemical Kinetics

The finite rate chemistry approach has been proven to be capable of handling turbulent combustion. However, a drawback with the finite rate method is that it requires high spatial resolution and it is sensitive to the selection of the reaction mechanism. It is necessary for the finite rate chemistry models to solve transport equations for each species. As a consequence, using a detailed mechanism can drastically increase the computational cost because a detailed mechanism usually contains tens of species and hundreds of intermediate reactions. The expensive computational cost caused by a detailed mechanism makes the choice of reaction mechanism a trade-off between computational capability and accuracy; also, it spurs the use of reduced mechanisms or an adaptable flamelet assumption approach [32–34] in complex situations to capture the flame characteristics. Several studies have been conducted to assess the predictive capabilities of the skeletal mechanisms and global mechanisms. Benim et al. [35] investigated a swirling flame in a model combustor by using the EDC model with a 2-step global mechanism. It is found that global mechanism provided a good overall accuracy compared with the experimental data. Fedina et al. [36] assessed several skeletal and global mechanisms in simulating premixed combustion in the SGT-100 dry low emission burner. The studies suggested that a reaction mechanism plays a significant role in numerical combustion, and the reduced reaction mechanisms are feasible in LES.In the present work, the 2-step global mechanism for methane oxidation is selected to handle the combustion and the chemical reactions are formulated as equations (15) and (16). Two sets of Arrhenius parameters are adopted to describe the chemical kinetics, one for the 2sCM2 scheme [37] (referred to as Mech I) and the other for the Arrhenius parameter-optimized 2-step scheme (referred to as Mech II) reduced based on the GRI 3.0.
(15)CH4+32o2→Co+2H2o(16)Co+12o2↔Co2The first reaction (15) is irreversible, the reaction rate is given by equation (17), and the kinetic parameters are reported in Table 1. The second reaction (16) is reversible and the reaction rates described by the two schemes are quite different, which were listed as equation (18) for Mech I and equation (19) for Mech II, where Ke is the equilibrium constant whereas kinetic parameters are reported in Table 2.
(17)ω̇1=A1Tβ1ρYCH4WCH4α1CH4ρYo2Wo2α1o2exp−Ea1RT,(18)ω̇2=A2fTβ2ρYCoWCoα2CoρYo2Wo2α2o2−1KeρYCo2WCo2α2Co2exp−Ea2RT,(19)ω̇2=A2fTβ2fρYCoWCoα2CoρYo2Wo2α2o2exp−Ea2fRT−A2rTβ2rρYCo2WCo2α2Co2exp−Ea2rRT.Table 1
Arrhenius parameters for the first reaction (units: kmol, m, s, kJ, and K).
A1
β1
Ea1
α
Mech I
2.0×1012
0
1.466×105
αCH4=0.9, αo2=1.1
Mech II
4.56×109
-0.0711
1.184×105
αCH4=0.5819, αo2=0.9223Table 2
Arrhenius parameters for the second reaction (units: kmol, m, s, kJ, and K).
A2
β2
Ea2
α
Mech I
6.32×107
0
51000
αCo=1.0, αo2=0.5, αCo2=1.0
Mech II
3.78×1011
0.221
77966
αCo=1.822, αo2=0.928
6.07×105
-0.108
132620
αCo2=1.073


## 3.4. Numerical Setup and Boundary Conditions

The three-dimensional computational domain with a diameter of 200 mm for the current simulation extends from 20 mm upstream of the bluff body to 300 mm downstream of the burner exit. The computational domain is discretized using nonuniform structured mesh. A section view of mesh is shown in Figure2. The mesh is refined near the bluff body and shear layers, which generates a corresponding minimal mesh space of 0.4 mm and overall cell number of about 1.6 million. The minimum grid spacing in the vicinity of the burner exit is considered a reasonable size to resolve the turbulence in the reacting zone [15]. The OpenFOAM reactingFoam solver is employed in the large eddy simulation of reaction cases, which is an unsteady reactive flow solver, and the thermophysical model for reacting mixture is based on compressibility ψ=RT−1. The transport property and thermodynamic property such as μ and Cp are calculated using Sutherland’s law and from JANAF tables of thermodynamics, respectively, and the energy variable used in the solution is based on sensible enthalpy hs.Figure 2
Grid resolution of the computational domain.At the inlets, Dirichlet conditions are adopted for all variables except for pressure. For pressure, zero Neumann condition is used. At the outlet, Dirichlet condition is used for pressure and zero Neumann conditions are adopted for other variables. The experimental data for velocity fluctuation is available only at the near exit, e.g., at 2 mm or 5 mm downstream of the burner exit [4, 5]. Uniform velocity boundary condition is adopted at the inlet, and turbulent intensity is controlled by κ and ε, which are estimated based on the near-exit velocity fluctuation data. All the walls are treated as no slip conditions with wall functions, and the lateral boundary of coflow is the freestream condition. Pressure-velocity coupling is performed using the PIMPLE algorithm, which is a combination of SIMPLE and PISO algorithms. The diffusion terms are discretized through the second-order central differencing scheme, and the convection terms are discretized through the second-order limited central differencing scheme. Temporal discretization is the first-order implicit Euler scheme with an adjustable time step to ensure that the maximal Courant number is below 0.5. The LES was initialized by three-dimensional steady-state RANS results. After the flow field fully developed, statistical sampling has been carried out for at least 5τ. τ is a characteristic domain flow through time and is defined as τ=L/Uj, where L is the streamwise domain length and Uj is the inner annulus jet velocity.


## 4. Results and Discussions

### 4.1. Nonreacting Flow Statistic Results

To assess the grid resolution, turbulent kinetic energy is collected at 30 mm downstream of the burner exit in the shear layer between the annular swirl and axial jets. The power spectrum is constructed by applying Fast Fourier Transform (FFT) to the instantaneous turbulent kinetic energy, as presented in Figure3. The comparison with the Kolmogorov -5/3 power law indicates the simulation is able to predict the characteristics of the inertial subrange. Mean velocity magnitude and its fluctuation in the region behind the bluff body are shown in Figure 4. Velocities at several representative locations downstream of the burner exit (∆X) are extracted to compare with experimental data, as presented in Figure 5.Figure 3
Power spectrum of turbulent kinetic energy.Figure 4
Mean (r<0mm) and rms (r>0mm) velocity behind the bluff body.Figure 5
Comparison of mean (a) and rms (b) axial velocity at∆X=2mm, 10 mm, 30 mm, and 50 mm (circles, experimental results; solid line, LES results).
(a)
(b)From the numerical results, the flow field is divided into three regions along the central axis. From∆X=0mm to 28 mm behind the bluff body, a closed primary recirculation zone (PRZ) is shown up along the axis, which indicates the occurrence of vortex breakdown of the inner annulus stream. The location of PRZ is far from the outer annulus swirl stream, and it is also found that the outer annulus swirl stream has minor effect on the formation of PRZ. From ∆X=28mm to 50 mm, the inner annular axial jet merges behind the PRZ and forms a positive velocity region. Above 50 mm, the annular swirl jet induces a centrifugal force, which causes the flow to expand and diverge outward from the central region and results in a secondary recirculation zone (SRZ). As shown in Figure 4, the fluctuation of velocity in shear layers downstream of PRZ shows slightly wavy patterns, which is also observed in experiment [5], and this phenomenon may be due to the processing vortex with swirl.In Figure5, at ∆X=2mm and 10 mm, mean velocities agree with experimental results well but in the central region, the velocities are slightly over predicted. There are two apparent steps, corresponding to the inner axial jet and outer annular swirl jet, respectively. Three velocity fluctuating peaks are captured from the center to the outside, corresponding to three shear layers, where intense flow mixing and sharp velocity gradient exist. With the mixing of flows, peaks of velocity fluctuation are smeared out gradually. At ∆X=30mm and 50 mm, underpredictions of velocity are observed in the central region. The underprediction of the axial velocity at downstream should contribute to the discrepancy in predicting the inner annulus stream. As shown in both experiments and simulation, for the nonreacting case, the inner annulus stream would converge downstream of the exit, which plays an important role in forming the closed primary recirculation zone. The experiment shows that the length of the recirculation zone behind the bluff body is about 15 mm. However, in this simulation, the length of the recirculation zone is about 27 mm, which is slightly longer than the measurement. As a result, the flow field downstream recirculation is affected.


### 4.2. Reacting Flow Statistic Results

For the reacting case, the LES predictions are compared with experimental data for velocity, temperature, and concentrations of major species such as CH4, CO, and CO2 at three typical locations (∆X=10mm, 30 mm, and 50 mm). Figure 6 shows the contour of mean velocity magnitude and its fluctuation. The mean velocity field of the reacting case is different from that of the nonreacting case as both axial and swirl jets diffuse outward the centerline due to the expansion of hot products. The behavior of recirculation zone is substantially affected by the presence of combustion. The nonreacting case shows a narrow and short conical recirculation zone. But in the reacting case, the shear zone moves rapidly outwards, which leads to the merging of closed PRZ with SRZ and generates big oval recirculation zone as denoted by the white solid line. For the fluctuation, combustion generates a big velocity fluctuating region within the flame zone, but the peak fluctuation falls slightly compared with the nonreacting case, which is believed to be caused by the relaminarization effect in the flame zone, where the fluctuation presented in the nonreacting case is suppressed in the reacting case [13]. Figure 7 shows the profiles of mean and fluctuating axial velocity. At ∆X=10mm, from a qualitative point of view, the trend of the mean and fluctuating velocity is reproduced, but the maximum recirculating velocity is underpredicted, and the predicted fluctuation is lower than the measurements. For the inlet inside the burner, the temporal velocity data is lacking. Uniform velocity boundary condition is adopted at the inlet, and turbulent intensity is controlled by κ and ε, which are estimated based on the near-exit velocity fluctuation data e.g., at 2 mm or 5 mm downstream of the burner exit. The simulation result shows that the flow at the exit may be not enough to establish a fully developed turbulence through 20 mm channel. This deviation on rms velocity prediction especially at the near exit (∆X=10mm) is caused by the velocity inlet boundary condition. At ∆X=30mm and 50 mm, the simulated velocities are too low compared with the measurements, which suggests that the reversed flow in the recirculation zone is not well reproduced. This is to a large extent in relation to the turbulence model. At ∆X=50mm, the velocity fluctuation predicted by both reacting simulations show slightly high values compared with that of the nonreacting case, but actually the fluctuation near the recirculation should have been suppressed.Figure 6
Mean (r<0mm) and rms (r>0mm) velocity behind the bluff body.Figure 7
Comparison of mean (a) and rms (b) axial velocity at∆X=10mm, 30 mm, and 50 mm (circles, experimental results; dash line, Mech I prediction; solid line, Mech II prediction).
(a)
(b)Figure8 shows the qualitative comparison of mean temperature distribution predicted by different mechanisms, and Figure 9 quantitatively compares numerical predictions with experimental data for temperature. From Figure 8, both chemical mechanisms predict similar flame shape and spreading flame angle at the base. The peak temperature is found in the leading edge of the flame near the bluff body, where fuel is consumed rapidly and heat release is intense. From Figure 9, within the flame zone, the mean and fluctuating temperature are underpredicted, but at the edge, where temperature fluctuation is high, the profiles are well reproduced. It is believed that the PaSR model can precisely capture the flame surface, but for the discrepancies in the central region at ∆X=10mm, a possible cause is the lack of many important intermediate radicals and reactions to depict the complex physical processes of combustion in the reacting zone. The two simulations exhibit a similar temperature profile at ∆X=30mm and 50 mm, which may indicate that deviation within the recirculation zone is caused by the turbulent combustion model.Figure 8
Comparison of simulated mean temperature distribution contour.Figure 9
Comparison of mean (a) and rms (b) temperature at∆X=10mm, 30 mm, and 50 mm (circles, experimental results; dash line, Mech I simulation; solid line, Mech II simulation).
(a)
(b)Figures10–12 compare simulation results with the experimental data for mean and fluctuation of major species. CH4 concentration is well reproduced with slight discrepancy in fluctuation. CH4 concentration rises sharply from zero in coflow to a maximal value with a step in a mixing layer (r=16mm), which is due to equivalence ratio stratification between the inner and outer annulus streams, whereas in the recirculation zone, CH4 concentration decreases towards the center with large gradient near the flame brush. At downstream, the peak of CH4 shifts outwards and the stratification step is smeared out. The highest fluctuation of CH4 is found near the flame brush. As expected, the same observations made for the temperature are evident for CO2 because the formation of CO2 is closely related to heat release. CO2 profiles are substantially affected by the flame brush behaviors. CO is a significant intermediate radical to determine chemical kinetics as well as to affect flame behaviors. It is necessary to take CO into account to evaluate the performance of chemical mechanisms. In Figure 12, CO as well as its fluctuation decreases sharply from the peak within the flame brush to zero by both sides of the flame brush, which indicates that CO is substantially exhibiting in a reacting layer with a pretty short resident time. Therefore, it is usually hard to be precisely captured. For the Mech I simulation, the mean and rms of CO are too low whilst Mech II simulation gives a reasonable prediction of CO concentration although the fluctuation is still too low. It indicates that the CO-CO2 equilibrium reaction in Mech II is more reasonable in this configuration. It also suggests that the flexibility of chemical mechanisms is deeply dependent on practical situations.Figure 10
Comparison of mean (a) and rms (b) CH4 mass fraction at ∆X=10mm, 30 mm, and 50 mm (circles, experimental results; dash line, Mech I simulation; solid line, Mech II simulation).
(a)
(b)Figure 11
Comparison of mean (a) and rms (b) CO2 mass fraction at ∆X=10mm, 30 mm, and 50 mm (circles, experimental results; dash line, Mech I simulation; solid line, Mech II simulation).
(a)
(b)Figure 12
Comparison of mean (a) and rms (b) CO mass fraction at∆X=10mm, 30 mm, and 50 mm (circles, experimental results; dash line, Mech I simulation; solid line, Mech II simulation).
(a)
(b)


### 4.3. Large-Scale Structures and Vortex Shedding

In premixed combustion, swirl-induced large-scale structures play a significant role in flame stabilization as well as flame surface wrinkling. To visualize the large-scale vortex structures,λ2-criterion [38] is applied to the velocity field to identify vortices. λ2-criterion identifies the vortices in the regions where the symmetric tensor S2+Ω2 has two negative Eigenvalues. S and Ω represent the symmetric and antisymmetric portions of the velocity gradient tensor matrix, which are defined as follows:
(20)S=12∇U+∇UT,(21)Ω=12∇U−∇UT.Figures13 and 14 show the vortex structures for the nonreacting and reacting cases illustrated using the isosurface of a single positive value of λ2. From the nonreacting flow field, there are two helical vortex structures which spiral in the same direction. The inner vortex originates from the shear layer between the annular axial and swirl jets, and the outer vortex structure originates from the shear layer between the annular swirl jet and coflow. The structures precess around the central axis with the same direction of swirl whilst they spin around their own axis. Near the burner exit, vortex structures are pretty regular, then at about ∆X=30~40 mm, large vortex structures break down into smaller ones, which distribute rather randomly and irregularly. In the inner shear layer, several streamlines are shown. From the figure, the streamlines and the vortex axis are orthogonal, which indicates that the vortex structures are generated because of Kelvin–Helmholtz instability of the shear layer [39]. And for the reacting flow field shown in Figure 14, the evolution of large vortex structures exhibit similar feature with that in the nonreacting flow field but have a larger spreading angle due to hot product expansion. An evident central vortex core is observed far downstream of the burner exit, whilst at the base, a PVC is generated caused by the swirl flow oscillation. The low-pressure vortex core moves circumferentially along the axis with a certain frequency.Figure 13
Instantaneous vortex structures (a) and PVC (b) for the nonreacting case.
(a)
(b)Figure 14
Instantaneous vortex structures (a) and PVC (b) for the reacting case.
(a)
(b)Vortex shedding is another important feature of flow instability, which may cause precession motion and flow field oscillation. To quantitatively investigate the vortex shedding, axial velocity point data at shear layer, namely, (∆X=25mm, r=13mm) for nonreacting flow and (∆X=25mm, r=17.5mm) for reacting flow, are sampled and analyzed through the FFT transform. The frequency is usually normalized by the Strouhal number (St) for its vital role played in detecting precession motion such as precessing vortex core (PVC) [40]. St is defined as st=2fRi/Ws, where f is the frequency, Ri is the inner radius of the annular swirl jet, and Ws is the tangential velocity of the swirl jet. The representative velocity power spectra of the monitoring points (∆X=25mm, r=13mm) and (∆X=25mm, r=17.5mm) are shown in Figure 15. The nonreacting flow exhibits two apparent characteristic frequency peaks at st=0.23 and st=1.25, respectively. From a qualitative point of view, the highest frequency peak at st=1.25 represents the occurrence of vortex shedding, and the lower frequency peak at st=0.23 represents the influence by the downstream precession motion, which is qualitatively depicted in Figure 13. The reacting flow shows the characteristic vortex shedding frequency at the st=2.45, and the representative frequency of precession motion at st=0.53, as presented in Figure 14.Figure 15
Power spectrum of instantaneous axial velocity and characteristic St for the nonreacting and reacting cases.
(a)
Nonreacting case
(b)
Reacting caseBy comparing the power spectra of the reacting case with the nonreacting case, both the amplitude and frequency increase in the reacting case, which indicates that the vortex shedding and precession motion are enhanced due to heat release. Namely, the unsteady characteristics and instability of flow field are intensified in presence of combustion.


### 4.4. Flame Behavior Analysis

The representative instantaneous temperature distribution predicted by the two mechanisms is illustrated in Figure16, which shows different flame behaviors between the mechanisms. In general, Mech I produces a compact high temperature zone within PRZ, and the flame surface tends to converge towards the centerline downstream of the PRZ, but in Mech II simulation, the flame surface extends farther. As expected, the species distribution in the downstream region is distinguishing, which may be attributed to the different chemical time scale described in two kinetics. In both simulations, flame surface wrinkling and high temperature pockets shedding from the flame surface are observed, which is caused by the unsteady characteristics of flow-flame interaction. It suggests that the PaSR model can precisely capture the unsteady features.Figure 16
Comparison of instantaneous temperature simulated using different mechanisms.Fuel stratification and its effect on flame behavior is analyzed based on Mech II simulation. Figure17 shows the sample simulation results comprising snapshots of instantaneous and mean CH4 mass fraction distributions. The mixtures are homogeneous at the two inlets; also, the mixing of inner and outer annular streams and the diluting of the outer annular stream with the coflow stream can be seen. Fuel stratification and gradient exhibit in the mixing layers.Figure 17
Snapshot of instantaneous and mean CH4 mass fraction of Mech II simulation.The equivalence ratio (ER) distribution at the same instance is shown in Figure18, in which the flame surface is derived from the heat release rate by solid line. The ER used in the experiment by Sweeney et al. [3] is derived from an atomic balance of carbon, hydrogen, and oxygen as follows:
(22)ϕ=XCo2+2XCH4+XCo+0.5XH2o+XH2XCo2+Xo2+0.5XCo+XH2o.Figure 18
Snapshot of instantaneous ER distribution of Mech II simulation (solid line, flame surface).It is based on the ratio of the demand for oxygen by local hydrogen and carbon to the locally available oxygen and is conserved from reactants to products across the premixed flame surface. From the CH4 concentration and ER distribution, three characteristic regions are distinguished within the flame zone. Below 25 mm, combustion first takes place at the inner annulus stream, where fuel is rich and mixtures are homogeneous, and ER is almost uniform. The combustion is stabilized by premixed reactants from the inner annulus stream. Flame interacts with the mixing layer and stratified mixtures from 25 mm to 50 mm, where the vortex breaking down enhances stream mixing and fuel diffusion, and CH4 gradient begins to appear in the reaction zone. As a result, the wrinkling of flame surface increases due to stratification effect. By comparing with experimental result [4], there is a noticeable expansion in the stratified region with the increase of swirl. Above 50 mm, the reaction zone starts to approach the flammability limit (ϕ≈0.5) region, where the flame angle decreases and flame tends to be unstable. Quenching/reignition happens and continuous flame surface break in the region. The characterized flame behavior in the swirling case is shown to be similar to that in the nonswirling case [14].It is worth noting that even in the so-called premixed combustion, the nonpremixed or partially premixed characteristics still exist in combustion processes. So any single combustion mode assumption in the complex flow field may induce certain error especially in flamelet-like model simulation. Combustion mode investigation could be helpful to improve the prediction accuracy of the flow field. Namely, by analyzing the combustion mode in different locations of flow field, reasonable flame regimes could be automatically introduced to describe turbulent combustion. Fortunately, some works have been reported on the combustion mode analysis method such as Takeno Flame Index (TFI) [41, 42], multiregime flamelet (MRF) [43], and Chemical Explosive Mode Analysis (CEMA) [44, 45].In this work, the TFI is used to describe the mixing mode between CH4 and O2, which is evaluated as G=∇YCH4∙∇Yo2. In the premixed reaction zone, G takes positive values, which means that CH4 and O2 are consumed in the same physical direction. Whilst in the nonpremixed reaction zone and out of the reaction zone, G takes negative values and zero value, respectively. A snapshot of G values contour plot from Mech II simulation is shown in Figure 19. As shown, near the bluff body, G takes positive values in the inner shear layer, which is the primary reaction zone according to Figure 18. The flame in this region exhibits evident premixed combustion characteristics. G takes negative values in small regions at the shear layer between the inner and outer stream, where equivalence gradient and flow mixing take place. These regions have minor impact on the flame. In the stratification region above 25 mm, nonpremixed characteristics are found in the reaction zone, which may be effected by equivalence stratification. In summary, it demonstrates that the combustion takes place in premixed mode dominated environment in the reaction zone, but nonpremixed characteristics still exhibit due to the stratification effect. On the basis of flame index analysis, it suggests that partially premixed model should be developed to improve the prediction accuracy.Figure 19
Snapshot of instantaneous and meanG value contour behind the bluff body (green, positive; purple, negative; and gray, zero).


## 4.1. Nonreacting Flow Statistic Results

To assess the grid resolution, turbulent kinetic energy is collected at 30 mm downstream of the burner exit in the shear layer between the annular swirl and axial jets. The power spectrum is constructed by applying Fast Fourier Transform (FFT) to the instantaneous turbulent kinetic energy, as presented in Figure3. The comparison with the Kolmogorov -5/3 power law indicates the simulation is able to predict the characteristics of the inertial subrange. Mean velocity magnitude and its fluctuation in the region behind the bluff body are shown in Figure 4. Velocities at several representative locations downstream of the burner exit (∆X) are extracted to compare with experimental data, as presented in Figure 5.Figure 3
Power spectrum of turbulent kinetic energy.Figure 4
Mean (r<0mm) and rms (r>0mm) velocity behind the bluff body.Figure 5
Comparison of mean (a) and rms (b) axial velocity at∆X=2mm, 10 mm, 30 mm, and 50 mm (circles, experimental results; solid line, LES results).
(a)
(b)From the numerical results, the flow field is divided into three regions along the central axis. From∆X=0mm to 28 mm behind the bluff body, a closed primary recirculation zone (PRZ) is shown up along the axis, which indicates the occurrence of vortex breakdown of the inner annulus stream. The location of PRZ is far from the outer annulus swirl stream, and it is also found that the outer annulus swirl stream has minor effect on the formation of PRZ. From ∆X=28mm to 50 mm, the inner annular axial jet merges behind the PRZ and forms a positive velocity region. Above 50 mm, the annular swirl jet induces a centrifugal force, which causes the flow to expand and diverge outward from the central region and results in a secondary recirculation zone (SRZ). As shown in Figure 4, the fluctuation of velocity in shear layers downstream of PRZ shows slightly wavy patterns, which is also observed in experiment [5], and this phenomenon may be due to the processing vortex with swirl.In Figure5, at ∆X=2mm and 10 mm, mean velocities agree with experimental results well but in the central region, the velocities are slightly over predicted. There are two apparent steps, corresponding to the inner axial jet and outer annular swirl jet, respectively. Three velocity fluctuating peaks are captured from the center to the outside, corresponding to three shear layers, where intense flow mixing and sharp velocity gradient exist. With the mixing of flows, peaks of velocity fluctuation are smeared out gradually. At ∆X=30mm and 50 mm, underpredictions of velocity are observed in the central region. The underprediction of the axial velocity at downstream should contribute to the discrepancy in predicting the inner annulus stream. As shown in both experiments and simulation, for the nonreacting case, the inner annulus stream would converge downstream of the exit, which plays an important role in forming the closed primary recirculation zone. The experiment shows that the length of the recirculation zone behind the bluff body is about 15 mm. However, in this simulation, the length of the recirculation zone is about 27 mm, which is slightly longer than the measurement. As a result, the flow field downstream recirculation is affected.


## 4.2. Reacting Flow Statistic Results

For the reacting case, the LES predictions are compared with experimental data for velocity, temperature, and concentrations of major species such as CH4, CO, and CO2 at three typical locations (∆X=10mm, 30 mm, and 50 mm). Figure 6 shows the contour of mean velocity magnitude and its fluctuation. The mean velocity field of the reacting case is different from that of the nonreacting case as both axial and swirl jets diffuse outward the centerline due to the expansion of hot products. The behavior of recirculation zone is substantially affected by the presence of combustion. The nonreacting case shows a narrow and short conical recirculation zone. But in the reacting case, the shear zone moves rapidly outwards, which leads to the merging of closed PRZ with SRZ and generates big oval recirculation zone as denoted by the white solid line. For the fluctuation, combustion generates a big velocity fluctuating region within the flame zone, but the peak fluctuation falls slightly compared with the nonreacting case, which is believed to be caused by the relaminarization effect in the flame zone, where the fluctuation presented in the nonreacting case is suppressed in the reacting case [13]. Figure 7 shows the profiles of mean and fluctuating axial velocity. At ∆X=10mm, from a qualitative point of view, the trend of the mean and fluctuating velocity is reproduced, but the maximum recirculating velocity is underpredicted, and the predicted fluctuation is lower than the measurements. For the inlet inside the burner, the temporal velocity data is lacking. Uniform velocity boundary condition is adopted at the inlet, and turbulent intensity is controlled by κ and ε, which are estimated based on the near-exit velocity fluctuation data e.g., at 2 mm or 5 mm downstream of the burner exit. The simulation result shows that the flow at the exit may be not enough to establish a fully developed turbulence through 20 mm channel. This deviation on rms velocity prediction especially at the near exit (∆X=10mm) is caused by the velocity inlet boundary condition. At ∆X=30mm and 50 mm, the simulated velocities are too low compared with the measurements, which suggests that the reversed flow in the recirculation zone is not well reproduced. This is to a large extent in relation to the turbulence model. At ∆X=50mm, the velocity fluctuation predicted by both reacting simulations show slightly high values compared with that of the nonreacting case, but actually the fluctuation near the recirculation should have been suppressed.Figure 6
Mean (r<0mm) and rms (r>0mm) velocity behind the bluff body.Figure 7
Comparison of mean (a) and rms (b) axial velocity at∆X=10mm, 30 mm, and 50 mm (circles, experimental results; dash line, Mech I prediction; solid line, Mech II prediction).
(a)
(b)Figure8 shows the qualitative comparison of mean temperature distribution predicted by different mechanisms, and Figure 9 quantitatively compares numerical predictions with experimental data for temperature. From Figure 8, both chemical mechanisms predict similar flame shape and spreading flame angle at the base. The peak temperature is found in the leading edge of the flame near the bluff body, where fuel is consumed rapidly and heat release is intense. From Figure 9, within the flame zone, the mean and fluctuating temperature are underpredicted, but at the edge, where temperature fluctuation is high, the profiles are well reproduced. It is believed that the PaSR model can precisely capture the flame surface, but for the discrepancies in the central region at ∆X=10mm, a possible cause is the lack of many important intermediate radicals and reactions to depict the complex physical processes of combustion in the reacting zone. The two simulations exhibit a similar temperature profile at ∆X=30mm and 50 mm, which may indicate that deviation within the recirculation zone is caused by the turbulent combustion model.Figure 8
Comparison of simulated mean temperature distribution contour.Figure 9
Comparison of mean (a) and rms (b) temperature at∆X=10mm, 30 mm, and 50 mm (circles, experimental results; dash line, Mech I simulation; solid line, Mech II simulation).
(a)
(b)Figures10–12 compare simulation results with the experimental data for mean and fluctuation of major species. CH4 concentration is well reproduced with slight discrepancy in fluctuation. CH4 concentration rises sharply from zero in coflow to a maximal value with a step in a mixing layer (r=16mm), which is due to equivalence ratio stratification between the inner and outer annulus streams, whereas in the recirculation zone, CH4 concentration decreases towards the center with large gradient near the flame brush. At downstream, the peak of CH4 shifts outwards and the stratification step is smeared out. The highest fluctuation of CH4 is found near the flame brush. As expected, the same observations made for the temperature are evident for CO2 because the formation of CO2 is closely related to heat release. CO2 profiles are substantially affected by the flame brush behaviors. CO is a significant intermediate radical to determine chemical kinetics as well as to affect flame behaviors. It is necessary to take CO into account to evaluate the performance of chemical mechanisms. In Figure 12, CO as well as its fluctuation decreases sharply from the peak within the flame brush to zero by both sides of the flame brush, which indicates that CO is substantially exhibiting in a reacting layer with a pretty short resident time. Therefore, it is usually hard to be precisely captured. For the Mech I simulation, the mean and rms of CO are too low whilst Mech II simulation gives a reasonable prediction of CO concentration although the fluctuation is still too low. It indicates that the CO-CO2 equilibrium reaction in Mech II is more reasonable in this configuration. It also suggests that the flexibility of chemical mechanisms is deeply dependent on practical situations.Figure 10
Comparison of mean (a) and rms (b) CH4 mass fraction at ∆X=10mm, 30 mm, and 50 mm (circles, experimental results; dash line, Mech I simulation; solid line, Mech II simulation).
(a)
(b)Figure 11
Comparison of mean (a) and rms (b) CO2 mass fraction at ∆X=10mm, 30 mm, and 50 mm (circles, experimental results; dash line, Mech I simulation; solid line, Mech II simulation).
(a)
(b)Figure 12
Comparison of mean (a) and rms (b) CO mass fraction at∆X=10mm, 30 mm, and 50 mm (circles, experimental results; dash line, Mech I simulation; solid line, Mech II simulation).
(a)
(b)


## 4.3. Large-Scale Structures and Vortex Shedding

In premixed combustion, swirl-induced large-scale structures play a significant role in flame stabilization as well as flame surface wrinkling. To visualize the large-scale vortex structures,λ2-criterion [38] is applied to the velocity field to identify vortices. λ2-criterion identifies the vortices in the regions where the symmetric tensor S2+Ω2 has two negative Eigenvalues. S and Ω represent the symmetric and antisymmetric portions of the velocity gradient tensor matrix, which are defined as follows:
(20)S=12∇U+∇UT,(21)Ω=12∇U−∇UT.Figures13 and 14 show the vortex structures for the nonreacting and reacting cases illustrated using the isosurface of a single positive value of λ2. From the nonreacting flow field, there are two helical vortex structures which spiral in the same direction. The inner vortex originates from the shear layer between the annular axial and swirl jets, and the outer vortex structure originates from the shear layer between the annular swirl jet and coflow. The structures precess around the central axis with the same direction of swirl whilst they spin around their own axis. Near the burner exit, vortex structures are pretty regular, then at about ∆X=30~40 mm, large vortex structures break down into smaller ones, which distribute rather randomly and irregularly. In the inner shear layer, several streamlines are shown. From the figure, the streamlines and the vortex axis are orthogonal, which indicates that the vortex structures are generated because of Kelvin–Helmholtz instability of the shear layer [39]. And for the reacting flow field shown in Figure 14, the evolution of large vortex structures exhibit similar feature with that in the nonreacting flow field but have a larger spreading angle due to hot product expansion. An evident central vortex core is observed far downstream of the burner exit, whilst at the base, a PVC is generated caused by the swirl flow oscillation. The low-pressure vortex core moves circumferentially along the axis with a certain frequency.Figure 13
Instantaneous vortex structures (a) and PVC (b) for the nonreacting case.
(a)
(b)Figure 14
Instantaneous vortex structures (a) and PVC (b) for the reacting case.
(a)
(b)Vortex shedding is another important feature of flow instability, which may cause precession motion and flow field oscillation. To quantitatively investigate the vortex shedding, axial velocity point data at shear layer, namely, (∆X=25mm, r=13mm) for nonreacting flow and (∆X=25mm, r=17.5mm) for reacting flow, are sampled and analyzed through the FFT transform. The frequency is usually normalized by the Strouhal number (St) for its vital role played in detecting precession motion such as precessing vortex core (PVC) [40]. St is defined as st=2fRi/Ws, where f is the frequency, Ri is the inner radius of the annular swirl jet, and Ws is the tangential velocity of the swirl jet. The representative velocity power spectra of the monitoring points (∆X=25mm, r=13mm) and (∆X=25mm, r=17.5mm) are shown in Figure 15. The nonreacting flow exhibits two apparent characteristic frequency peaks at st=0.23 and st=1.25, respectively. From a qualitative point of view, the highest frequency peak at st=1.25 represents the occurrence of vortex shedding, and the lower frequency peak at st=0.23 represents the influence by the downstream precession motion, which is qualitatively depicted in Figure 13. The reacting flow shows the characteristic vortex shedding frequency at the st=2.45, and the representative frequency of precession motion at st=0.53, as presented in Figure 14.Figure 15
Power spectrum of instantaneous axial velocity and characteristic St for the nonreacting and reacting cases.
(a)
Nonreacting case
(b)
Reacting caseBy comparing the power spectra of the reacting case with the nonreacting case, both the amplitude and frequency increase in the reacting case, which indicates that the vortex shedding and precession motion are enhanced due to heat release. Namely, the unsteady characteristics and instability of flow field are intensified in presence of combustion.


## 4.4. Flame Behavior Analysis

The representative instantaneous temperature distribution predicted by the two mechanisms is illustrated in Figure16, which shows different flame behaviors between the mechanisms. In general, Mech I produces a compact high temperature zone within PRZ, and the flame surface tends to converge towards the centerline downstream of the PRZ, but in Mech II simulation, the flame surface extends farther. As expected, the species distribution in the downstream region is distinguishing, which may be attributed to the different chemical time scale described in two kinetics. In both simulations, flame surface wrinkling and high temperature pockets shedding from the flame surface are observed, which is caused by the unsteady characteristics of flow-flame interaction. It suggests that the PaSR model can precisely capture the unsteady features.Figure 16
Comparison of instantaneous temperature simulated using different mechanisms.Fuel stratification and its effect on flame behavior is analyzed based on Mech II simulation. Figure17 shows the sample simulation results comprising snapshots of instantaneous and mean CH4 mass fraction distributions. The mixtures are homogeneous at the two inlets; also, the mixing of inner and outer annular streams and the diluting of the outer annular stream with the coflow stream can be seen. Fuel stratification and gradient exhibit in the mixing layers.Figure 17
Snapshot of instantaneous and mean CH4 mass fraction of Mech II simulation.The equivalence ratio (ER) distribution at the same instance is shown in Figure18, in which the flame surface is derived from the heat release rate by solid line. The ER used in the experiment by Sweeney et al. [3] is derived from an atomic balance of carbon, hydrogen, and oxygen as follows:
(22)ϕ=XCo2+2XCH4+XCo+0.5XH2o+XH2XCo2+Xo2+0.5XCo+XH2o.Figure 18
Snapshot of instantaneous ER distribution of Mech II simulation (solid line, flame surface).It is based on the ratio of the demand for oxygen by local hydrogen and carbon to the locally available oxygen and is conserved from reactants to products across the premixed flame surface. From the CH4 concentration and ER distribution, three characteristic regions are distinguished within the flame zone. Below 25 mm, combustion first takes place at the inner annulus stream, where fuel is rich and mixtures are homogeneous, and ER is almost uniform. The combustion is stabilized by premixed reactants from the inner annulus stream. Flame interacts with the mixing layer and stratified mixtures from 25 mm to 50 mm, where the vortex breaking down enhances stream mixing and fuel diffusion, and CH4 gradient begins to appear in the reaction zone. As a result, the wrinkling of flame surface increases due to stratification effect. By comparing with experimental result [4], there is a noticeable expansion in the stratified region with the increase of swirl. Above 50 mm, the reaction zone starts to approach the flammability limit (ϕ≈0.5) region, where the flame angle decreases and flame tends to be unstable. Quenching/reignition happens and continuous flame surface break in the region. The characterized flame behavior in the swirling case is shown to be similar to that in the nonswirling case [14].It is worth noting that even in the so-called premixed combustion, the nonpremixed or partially premixed characteristics still exist in combustion processes. So any single combustion mode assumption in the complex flow field may induce certain error especially in flamelet-like model simulation. Combustion mode investigation could be helpful to improve the prediction accuracy of the flow field. Namely, by analyzing the combustion mode in different locations of flow field, reasonable flame regimes could be automatically introduced to describe turbulent combustion. Fortunately, some works have been reported on the combustion mode analysis method such as Takeno Flame Index (TFI) [41, 42], multiregime flamelet (MRF) [43], and Chemical Explosive Mode Analysis (CEMA) [44, 45].In this work, the TFI is used to describe the mixing mode between CH4 and O2, which is evaluated as G=∇YCH4∙∇Yo2. In the premixed reaction zone, G takes positive values, which means that CH4 and O2 are consumed in the same physical direction. Whilst in the nonpremixed reaction zone and out of the reaction zone, G takes negative values and zero value, respectively. A snapshot of G values contour plot from Mech II simulation is shown in Figure 19. As shown, near the bluff body, G takes positive values in the inner shear layer, which is the primary reaction zone according to Figure 18. The flame in this region exhibits evident premixed combustion characteristics. G takes negative values in small regions at the shear layer between the inner and outer stream, where equivalence gradient and flow mixing take place. These regions have minor impact on the flame. In the stratification region above 25 mm, nonpremixed characteristics are found in the reaction zone, which may be effected by equivalence stratification. In summary, it demonstrates that the combustion takes place in premixed mode dominated environment in the reaction zone, but nonpremixed characteristics still exhibit due to the stratification effect. On the basis of flame index analysis, it suggests that partially premixed model should be developed to improve the prediction accuracy.Figure 19
Snapshot of instantaneous and meanG value contour behind the bluff body (green, positive; purple, negative; and gray, zero).


## 5. Conclusion

Large eddy simulations are carried out on Cambridge swirl burner moderately stratified configuration with both the nonreacting and reacting cases. Combustion is handled using a partially stirred reactor model associated with reduced chemical mechanisms. Based on flow and flame investigations, the main conclusions are as follows:(1)
LES predictions are compared with measurements for mean and fluctuation of velocity, temperature, and major species concentrations. From a qualitatively point of view, the velocity profiles are reproduced. The shape of the recirculation zone changes obviously with the combustion from cone to oval. Deviations are found in temperature and major species concentrations occurring near the bluff body region
(2)
The location of the flame surface and flame surface wrinkling are well captured by the PaSR model. Two mechanisms showed evident difference in CO prediction and flame behavior, which indicates that predictive capabilities of reduced mechanisms are still limited due to the essentially complex physical processes in combustion and the flexibility of chemical mechanisms is deeply dependent on practical situations
(3)
Large-scale vortex structures are observed in nonreacting and reacting cases, and the evolution of large-scale structure, which is related to PVC formation, is illustrated in detail. Precession motions are found occurring at a certain frequency in both nonreacting and reacting cases. Based on the axial velocity point data at the shear layer, both the intensity and frequency of precession motion and vortex shedding are shown to be enhanced in the presence of combustion
(4)
Three characteristic reaction regions are distinguished in the flame zone, and the stratification effects on flame surface wrinkling are investigated. Near the bluff body, the combustion is stabilized in premixed reactants from the inner annulus stream and shows premixed combustion characteristics. Whilst in the equivalence gradient region, the flame surface is sensitive to stratification and flame surface wrinkling increases due to the stratification effect



---

*Source: 7042692-2019-03-31.xml*