Large eddy simulation of premixed hydrogen-rich gas turbine combustion based on reduced reaction mechanisms

Abstract

This study numerically investigates the combustion characteristics in a swirl burner for CO-H₂ syngas by large eddy simulation. Firstly, experimental data are used to verify the applicability of the Boivin mechanism and the Sandia mechanism. The simulation results show that the Boivin mechanism predicts a more reasonable flame structure. Then, Syngas25, Syngas45, and Syngas100 (composed of 25%, 45%, and 100% H₂ by volume, respectively) are set to numerically investigate the effect of hydrogen addition ratio on combustion characteristics by using the dynamic Smagorinsky model and the partially stirred reactor model. The detailed simulation data shows that the addition of hydrogen greatly influences the flame structure and flow field; with an increase in hydrogen addition, the flame length shortens, the temperature rises, and flashback caused by combustion induced vortex breakdown occurs. Furthermore, the precessing vortex core cannot exist stably in Syngas100, though it can be clearly observed in other cases.

KEYWORDS:

• Large eddy simulation
• syngas
• premixed combustion
• reduced reaction mechanism
• swirl flame
• hydrogen addition

Nomenclature list

$\rho$	=	Density
$t$	=	Time
$u_i$	=	Component of velocity vector in the i direction, i= 1, 2, 3.
$x_i$	=	Component of spatial coordinate in the i direction, i= 1, 2, 3.
$p$	=	Pressure
${\tau }_{ij}$	=	Stress tensor
${\tau }_{sgs}$	=	Subgrid shear stress tensor
$Y_s$	=	Mass fraction of the species s
${\dot{Y}}_s$	=	Source term of the species s
${\dot{Q}}_s$	=	Energy source term
$D_s$	=	Mass diffusion coefficient of the species s
${\mu }_t$	=	Subgrid viscosity
$k_{sgs}$	=	Subgrid kinetic energy
${\tilde{S}}_{ij}$	=	Rate-of-strain tensor
${\mathrm{\Theta }}_{sgs}$	=	Sub-grid fluxes of species
$h_j$	=	Enthalpy
$J_j$	=	Mass fluxes
$T$	=	Temperature
$H_{sgs}$	=	Sub-grid fluxes of energy
$S$	=	Swirl number
$\phi$	=	Equivalence ratio
${\nu }_t$	=	Subgrid kinematic viscosity
$P{r}_t$	=	Turbulent Prandtl number
$S{c}_t$	=	Turbulent Schmidt number
${\omega }_s$	=	Global reaction rate in the PaSR model
${\omega }_a$	=	Unrevised global reaction rate
${\tau }_c$	=	Chemical time scale
${\tau }_t$	=	Turbulent time scale
$\nu$	=	Kinematic viscosity
$\Delta$	=	Mesh scale
$s_l$	=	Laminar flame speed
$E$	=	Total energy
$Le$	=	Lewis number
$Re$	=	Reynolds number
$Ka$	=	Karlovitz number
$CFL$	=	Courant–Friedrichs-Lewy number
$D$	=	Diameter of center body
$d$	=	Diameter of cylindrical combustion chamber
$t_{cycle}$	=	cycle time of swirling

1. Introduction

Despite the supply of energy by using different low-cost resources with lower carbon footprints is an on-going research interest (Eivaz et al., 2018), the major global energy supply still depend on fossil fuels (Jafari-Sejahrood et al., 2019). With the environmental protection regulations for power plants becoming increasingly stringent in recent years (Issakhov et al., 2020), integrated gasification combined cycle (IGCC) using coal gasification fuel and gas turbines has received much research attention due to its low emissions and high efficiency (E.O. Oluyede, 2007). The coal gasification fuel called syngas is a mixture of H₂, CO, N₂, and other minor constituents. Depending on the coal producer and gasification process, the relative proportion of H₂ and CO can vary widely. The changeable fraction ratio thus induces variable combustion characteristics for syngas.

H₂-CO syngas has a low volumetric heating value, about a third of methane (S. Taamallah et al., 2015). H₂, one of the main components of syngas, has high gravimetric heating value and little pollution released from its combustion (in the presence of sufficient oxygen) that produces only water and energy (Faizollahzadeh Ardabili et al., 2018). Most hydrogen is currently produced via steam reforming of natural gas. As a fuel, hydrogen leads to major operability concerns in premixed combustion, because it is different from standard fuels on chemical and transport properties (S. Taamallah et al., 2015).

Among available fuels, hydrogen stands out in terms of its transport properties. Heat and momentum diffusivity for hydrogen are over six times larger than that for CH₄, CO, and air, and mass diffusivity for hydrogen in air is about four times larger than CH₄ and CO in air. CH₄ and CO have similar transport properties, but different densities (S. Taamallah et al., 2015). Thus, CO-H₂ syngas should exhibit different combustion properties compared with CH₄ and CO.

Wang et al. (2009) found that using hydrogen, especially pure hydrogen, in fuel mixtures is prone to severe flashback. However, when CH₄ is added to H₂, flashback probability becomes lower. Daniele et al. (2009) concluded that in gas turbine, compared with natural gas, hydrogen-rich syngas are characterized by higher burning rates and, consequently, higher flashback propensity.

The Lewis number (Le) is defined as the ratio of thermal diffusivity to mass diffusivity (Abadi et al., 2020); for fuels with a small Le, the chemical energy contributed by the mass diffusion to the reaction zone is much higher than the energy lost from the reaction zone due to thermal diffusion, resulting in the formation of a local high reactant concentration zone and high temperature zone in the flow field (Fu et al., 2014). This phenomenon is caused by preferential diffusion. Park et al. (2009) analysed the effect of preferential diffusion on temperature for non-premixed flames of syngas and found that chemical reactions are closely related to the formation of a local temperature distribution in the flow field. Bazooyar and Darabkhani (2019) studied premixed hydrogen flame stabilization in a micro-combustor, and evaluated effect of preferential transport of species. Daniele et al. (2011) successfully investigated the influence of the hydrogen concentration of syngas on turbulent burning velocity on the basis of OH-PLIF data. The results confirmed that hydrogen addition raises the turbulent-to-laminar flame speed ratio at the same turbulent level, for preferential diffusive-thermal effects.

In burning syngas in a swirling premixed chamber, the flow and combustion characteristics are important for chamber performance and should thus be investigated systemically. Using large eddy simulation (LES) to study combustion chambers has become an important auxiliary design method for the combustion chambers of gas turbines. In the LES method, the large-scale eddies are resolved directly while the small-scale eddies are modeled, so that the dominant vortex structures in the turbulence flow can be computed with limited computational expense (XSun et al., 2020). LES can reduce the number of experimental conditions, reduce the time and resource cost of experiments, and obtain data that are otherwise difficult to acquire in experiments. In addition, it can capture the dynamic process of the interaction between combustion and turbulence. Huang used LES in studying the dynamic combustion processes in a conventional swirling chamber that burns methane (Huang & Yang, 2009). The results revealed that the precessing vortex core (PVC) strongly affects the flow and flame evolution in combustion systems. Zheng et al. (2013) performed an LES of the premixed swirling combustion of syngas by using reduced chemistry based on the Sandia mechanism. Their study revealed that the hydrogen concentrations of syngas does not affect the characteristic frequencies of the hydrodynamic instabilities. De and Acharya (2012b) found that an LES-based artificially thickened flame model, even with simplified chemistry (3 reactions/6 species), is capable of predicting the behavior of 30% hydrogen-addition methane flame in a swirl burner. Another research by De and Acharya (2012a) showed that high hydrogen content in CH₄-H₂ syngas makes the central toroidal recirculation zone (CTRZ) compact, leading to flashback for hydrogen has the high diffusivity and reactivity. Combustion-induced vortex breakdown (CIVB) was the main reason for the flashback observed numerically.

The behavior of CIVB depends on several factors, including the location of vortex breakdown, the inlet temperature, and axial velocity (Tangermann et al., 2010). Dam et al. (2011) experimentally investigated the flashback mechanism caused by CIVB. Syngas with various compositions and different swirl numbers were analyzed. The results indicated that flashback is more sensitive to fuel composition. Under the same air flow and swirl number, syngas with higher H₂ fraction leads to flashback at considerably lean conditions.

Although extensive investigations into premixed hydrogen-rich syngas using LES have been performed, most studies focused on CH₄-H₂ syngas. The current research on the effect of hydrogen addition on the combustion characteristics of CO-H₂ syngas has not fully taken advantages of LES, that is, obtaining instantaneous results. As CIVB occurs primarily in swirl-stabilized combustors without a center body (Tangermann et al., 2010), the number of studies using LES on flashback caused by CIVB in combustors with center body is limited. The objective of this study is thus to numerically investigate the effect of hydrogen addition ratio in CO-H₂ syngas on combustion characteristics via LES.

2. Flow configuration

2.1. Experimental device

To study the syngas flame in a swirl burner on the basis of different reaction mechanisms, this study used a baseline case with a steady inflow and experimental data for simulation validation. A new swirl burner designed by (Wang et al., 2016; Yu et al., 2013) was used to study the combustion characteristics. The burner called arrayed-vanes combustor is composed of an arrayed-vanes premixer and a chamber.

The arrayed-vanes combustor was experimentally studied by Yu et al. (2013). In the experiment, the arrayed-vanes premixer was assembled on a cylindrical chamber made of quartz glass, which facilitated the optical measurement. A typical operating condition of the combustor was set using syngas fuel, and the operating condition under this base case is summarized in Table 1. Under the current condition, the global equivalence ratio φ = 0.29, and the swirl number S = 0.54.

Table 1. Operating condition of the experiment.

	Volume flow rate m^3/s	Temperature K	Component Vol%
Air	0.4202	694	O2: 21, N2: 79
Syngas Exp	0.0306	293	H2:36.7, CO:44.4, CO2: 18.9

In the experiment, the combustion proceeded stably without flashback or blow out. The time-averaged temperature in the center of the outlet of the chamber was 1,393 K, which was measured by a thermocouple. The image of the flame in the chamber was captured by a digital camera (Figure 1). The mean flame length in the chamber was about 82 mm, which was measured by scaling the luminous area with the chamber size in the picture. However, given the difficulty of obtaining optical measurements in such a complicated chamber, the details of the distribution of the flow variables were not obtained. The dynamic combustion processes and the interaction between the flame and turbulent vortexes could not be obtained either. Thus, the LES of the swirl burner was conducted in the current work to reveal the complex process in the chamber. The accuracy of the results of the simulation was validated by the experimental data.

Figure 1. Flame image in experiment.

2.2. Computational domain

Since no reaction occurs in the complex arrayed-vanes premixer, to accelerate the computation process, we simulated the arrayed-vanes premixer and the chamber individually (Figure 2). LES of the premixer was firstly conducted, the simulation collected timing sequence data after achieving full development, including velocity, temperature, and species concentration on the premixer outlet plane. These data, which contain details of the outflow of the premixer, were used as the inlet boundary condition for the chamber. A similar approach was used successfully by (Tangermann et al., 2010; Tangermann & Pfitzner, 2008), who adopted a steady profile of velocity as the inlet boundary condition. By contrast, the dynamic inlet boundary condition adopted in the present research is more reasonable compared with the steady inlet boundary condition, because it can obtain the turbulence structure.

Figure 2. Computational domains of (a) premixer and (b) combustion chamber.

On the other hand, for simulating the combustion characteristics in the chamber, the computation domain for the LES encompasses an inflow passage with a diameter of 100 mm and a length of 50 mm and a cylindrical combustion chamber with a diameter of 210 mm and a length of 450 mm.

2.3. Mesh and boundary

An unstructured mesh system with 3.4 million grids was employed in the premixer. The mean grid scale was 0.5 mm, which is sufficient for LES in previous research. The chamber adopted a hexahedral mesh with optimization in the area where reactions occur. To check the grid independence, we employed three different strategies with 0.5, 1.5, and 2.5 million cells. Figure 3 illustrates the profiles of the mean axial velocity and temperature obtained by the three meshes. The coarse mesh yielded an obviously different result from the others, whereas the data from the medium mesh were similar to those from the fine mesh. Thus, the mesh with 1.5 million grids for original model was used in this study.

Figure 3. Profiles of (a) mean axial velocity and (b) temperature profiles predicted by different mesh schemes.

In the simulation of the combustor, the syngas and air inlet were specified according to Table 1, and the walls were regarded as the nonslip adiabatic boundary. In the simulation of the combustion chamber, the inflow data from the outlet of the premixer was adopted, and the walls were regarded as nonslip walls with constant temperature.

Running on an Intel cluster, an in-house source code package developed by our workgroup was applied to carry out the simulation. The code employs a semi-implicit Crank–Nicolson difference method with second-order accuracy as time–space coupled discrete method. For pressure term, velocity convection term and velocity viscous term, the code employs implicit difference methods with second-order accuracy. For scalar convection term, the code employs an implicit WENO difference method with third-order accuracy. Given the huge amount of computation involved in the LES of the chamber, a parallel strategy was adopted herein.

The simulation of each case was run until the flow time reached 2.0 s to obtain full development. Then, the statistical data were obtained in another 0.5 s. The time steps used in all cases were 2×10⁻⁶ s, and the max value of CFL is 0.6, guaranteed the computational stability and captured the dynamic processes of combustion. The simulations were completed on a massively parallel cluster, with one whole case costing about 2,400 CPU hours to obtain the final results.

2.4. Simulation cases

The original model of combustion chamber was simplified to further improve the calculation efficiency. The simplified model ignores two coaxial air inlets, which account for approximately 5% of the total air flow. As shown in Figure 7, given the low air velocity and no reaction occurs, their impact on the mainstream area can be neglected. Similar grid independence check has also been performed for the simplified model. Here, the mesh with 1.06 million grids for simplified model was used in this study.

Figure 4 shows Sketch of the mesh for combustion chamber for original model and simplified model. Herein, D = 0.12 m, and the profiles of four axial locations are shown in Figure 4.Three simulation cases were set to investigate the effect of hydrogen addition on combustion characteristics. The details are listed in Table 2. By analyzing the results of the three cases, we can study the combustion characteristics and flame dynamics caused by the addition of hydrogen.

Figure 4. Sketch of the mesh for combustion chamber for original model (above) and simplified model (below).

Table 2. Conditions of simulation cases.

	Volume flow rate m^3/s	Temperature K	Component Vol%
Air	0.3946	694	O2: 21, N2: 79
Syngas25	0.0306	293	H2:25, CO:75, CO2: 0
Syngas45	0.0306	293	H2:45, CO:55, CO2: 0
Syngas100	0.0306	293	H2:100, CO:0, CO2: 0

3. Numerical methods

3.1. Subgrid model and governing equations

In this study, the combustion process was simulated using LES, which could capture most of the turbulent energy in large eddies. To describe the transportation of mass, momentum, species, and energy in three dimensions, we used the following Favre-filtered conservation equations: (1) $\begin{array}{rl}& \frac{\mathrm{\partial }\overline{\rho }}{\mathrm{\partial }t}+\frac{\mathrm{\partial }\overline{\rho }{\tilde{u}}_i}{\mathrm{\partial }{x}_i}=0\end{array}$(1) (2) $\begin{array}{rl}& \frac{\mathrm{\partial }\overline{\rho }{\tilde{u}}_i}{\mathrm{\partial }t}+\frac{\mathrm{\partial }\overline{\rho }{\tilde{u}}_i{\tilde{u}}_j}{\mathrm{\partial }{x}_j}=-\frac{\mathrm{\partial }\overline{p}}{\mathrm{\partial }{x}_i}+\frac{\mathrm{\partial }({\tau }_{ij}+{\tau }_{sgs})}{\mathrm{\partial }{x}_j}\end{array}$(2) (3) $\begin{array}{rl}& \frac{\mathrm{\partial }\overline{\rho }{\tilde{Y}}_s}{\mathrm{\partial }t}+\frac{\mathrm{\partial }\overline{\rho }{\tilde{u}}_j{\tilde{Y}}_s}{\mathrm{\partial }{x}_j}=\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left(\overline{\rho }{D}_s\frac{\mathrm{\partial }{\tilde{Y}}_s}{\mathrm{\partial }{x}_j}+{\mathrm{\Theta }}_{sgs}\right)+{\dot{Y}}_s\end{array}$(3) (4) $\begin{array}{rl}& \frac{\mathrm{\partial }\overline{\rho }\overline{E}}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}[(\overline{\rho }\overline{E}+\overline{p}){\tilde{u}}_j]\\ & =\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left(\lambda \frac{\mathrm{\partial }\tilde{T}}{\mathrm{\partial }{x}_j}-\sum h_j{J}_j+H_{sgs}\right)+{\dot{Q}}_s\end{array}$(4)

In these equations, ‘∼’ and ‘-’ denote the Favre-filtered and Reynolds-filtered variables, respectively. The subgrid stress tensor in Equation (2), ${\tau }_{sgs}$, is calculated using a dynamic Smagorinsky model shown in Equation (5), where the subgrid viscosity, ${\mu }_t$, and the subgrid kinetic energy, $k_{sgs}$, are calculated dynamically (Germano et al., 1991). Subgrid species and energy fluxes are approximated by the gradient model, as shown in Equations (6) and (7): (5) $\begin{array}{rl}& \{\begin{array}{l}{\tau }_{sgs}=2{\mu }_t{\tilde{S}}_{ij}-\frac{2}{3}{k}_{sgs}{\delta }_{ij}\\ {\tilde{S}}_{ij}=\frac{1}{2}\left({\displaystyle \frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{x}_j}}+{\displaystyle \frac{\mathrm{\partial }{u}_j}{\mathrm{\partial }{x}_i}}\right)\end{array}\end{array}$(5) (6) $\begin{array}{rl}& {\mathrm{\Theta }}_{sgs}=\frac{v_t}{S{c}_t}\frac{\mathrm{\partial }{\tilde{Y}}_s}{\mathrm{\partial }{x}_i}\end{array}$(6) (7) $\begin{array}{rl}& H_{sgs}=\frac{v_t}{P{r}_t}\frac{\mathrm{\partial }\tilde{T}}{\mathrm{\partial }{x}_i}\end{array}$(7)

The chemical source species term, ${\dot{Y}}_s$, and chemical source energy term, ${\dot{Q}}_s$, in Equations (3) and (4) could be obtained from reaction rate, ${\omega }_s$, which is calculated by the PaSR model (Berglund et al., 2010). This model has been proved to be accurate in the LES of swirling premixed combustion, and its expression is defined in the following section. The source term of the energy introduced by the radiation in Equation (4), is neglected.

3.2. PaSR model

In the LES of combustion, various models are proposed to solve the interaction between the turbulence and the chemical reaction. One intuitive approach is the finite chemistry model, in which the reaction rate is calculated using the well-known Arrhenius equation. Every variable in the Arrhenius equation can be decomposed into filtered and subgrid components after filtering.

The PaSR model is considered as an accurate finite rate model. It regards each computing grid as a small reactor, where fuel and oxygen are mixed and chemical reactions occur. As the thickness of the area where combustion occurs is smaller than the grid scale, only a part of the components in the grid participates in the chemical reaction. In the PaSR model, each LES cell can be divided into fine structures, where mixing and reactions are assumed to take place, and surroundings, dominated by large-scale coherent flow structures. The PaSR model expresses the influence of turbulence in the form of superimposed correction on the calculated value of the chemical reaction rate, instead of simply using the flow field parameters to calculate the chemical reaction rate. The model calculates the chemical reaction rate of each reaction step separately, and it can use multistep reaction mechanisms.

The PaSR model ignores scalar pulsation in fine structures, thus the effect of subgrid pulsations on chemical reaction rates was not considered in the model. Researchers have widely applied the PaSR model to the combustion field, and the most prediction of PaSR-LES is consistent with the experimental data (Berglund et al., 2010).

Global reaction rate in the PaSR model, ${\omega }_s$, is defined as Equation (8). (8) ${\omega }_s\text{ = }\frac{{\tau }_c}{{\tau }_c+{\tau }_t}\cdot {\omega }_a$(8)

In Equation (8), ${\omega }_a$ is calculated by the Arrhenius equation, and the chemical time scale, ${\tau }_c$, and turbulent time scale, ${\tau }_t$, are defined as follows: (9) $\begin{array}{rl}{\tau }_c& =\frac{\upsilon }{s_l^2}\end{array}$(9) (10) $\begin{array}{rl}{\tau }_t& ={\upsilon }^{0.25}{\mathrm{\Delta }}^{0.75}{k_{sgs}}^{0.625}\end{array}$(10)

3.3. Reaction mechanism

Extensive studies have explored premixed syngas combustion, especially in syngas combustion mechanisms. Marzouk and David Huckaby (2010) compare the performance and computational cost of 8 syngas combustion mechanisms by using the PaSR model and a modified k-epsilon model. Carsten Olm et al. (2015) tested 16 syngas combustion mechanisms against experimental data, and investigated the dependence of their predictions on the types of experiment and the experimental conditions.

Comparison of syngas combustion mechanisms are listed in Table 3. Syngas combustion mechanisms could be divided into detailed mechanisms and reduced mechanisms. Using detailed mechanisms can predict reactions more accurately, whereas using reduced mechanisms can reduce more computing resources and time. The time cost is at least a cubic function of the number of species, even a minor reduction in the number of species can result in substantial reduction in the computation time (Lu & Law, 2009).

Table 3. Comparison of syngas combustion mechanisms in recent years.

Mechanism	Ref.	Species no	Reactions no	Type
Kéromnès-2013	(Keromne et al., 2013)	15	49	Detailed
NUIG-NGM-2010	(Kathrotia et al., 2010)	15	41	Detailed
Li-2007	(Li et al., 2007)	15	45	Detailed
Sun-2007	(H. Sun et al., 2007)	15	48	Detailed
Davis-2005	(Davis et al., 2005)	14	38	Detailed
Sandia-2002	(Chen, 2002)	10	6	Reduced
Boivin-2011	(Boivin et al., 2011)	8	4	Reduced
Watanabe-2006	(Watanabe & Otaka, 2006)	5	3	Reduced

In the PaSR model, all species are solved. Hence, accurate reduced mechanisms should be adopted in the simulation to reduce the time cost. In the present study, Boivin mechanism and Sandia mechanism were adopted after measuring computing resources cost and prediction accuracy. The results were then compared with the experimental outcomes to verify their accuracy.

Boivin mechanism involves 4 reaction steps and 8 species (Boivin et al., 2011). A previous study showed that the mechanism can obtain reasonable flame speed and ignition time under various operating conditions. This mechanism is written as follows: $\{\begin{array}{l}3H_2+O_2=2H_{2}O+2H\\ 2H+M=H_2+M\\ H_2+O_2=H{O}_2+H\\ CO+H_{2}O=C{O}_2+H_2\end{array}$

The Sandia mechanism is presented on the website of Sandia National Lab with 6 reaction steps and 10 species (Chen, 2002). Its accuracy was verified through a comparison with the results of a detailed mechanism. This mechanism is written as follows: $\{\begin{array}{l}{O}_2+H=OH+O\\ H_2+O=H+OH\\ H_2+OH=H+H_{2}O\\ 2H=H_2\\ CO+OH=H+C{O}_2\\ 2NO=O_2+N_2\end{array}$

The same working condition as that presented in Table 1 was set to verify the applicability of the two reaction mechanisms. The results are detailed in Section 4.1.

4. Results and discussion

4.1. Applicability of the two reaction mechanisms

To verify the applicability of the two reaction mechanisms, we conduct LES with the same working conditions as those in Table 1.

As the H radical is the most important intermediate product in the Boivin and Sandia mechanisms, the mass fraction of the H radical can reflect the reaction intensity in the chamber and is usually used to indicate flame structure. The isosurfaces of the H radical in the chamber obtained in both cases are shown in Figure 5. The Boivin mechanism could predict a similar flame structure to that predicted in the experiment (Figure 1). The flame exhibits a double-layer structure and is stabilized near the wall. However, the flame structure predicted by the Sandia mechanism is blown off, which indicates a low reaction intensity. The simulation result from the Sandia mechanism is obviously different from the stable combustion state shown in the experiment.

Figure 5. Isosurfaces of H radical obtained by (a) Boivin mechanism and (b) Sandia mechanism.

Figure 6 shows the time-averaged temperatures obtained by the two mechanisms. The temperatures on the outlet predicted by both mechanisms are approximately 1,420 K, which is close to the experimental data (1,393 K). However, the Sandia mechanism obtains a longer and wider low-temperature region than the Boivin mechanism. This result coincides with the different prediction of flame structure in Figure 5. The mid-temperature, $T_{\text{mid}}=(1/2)(T_{\text{inlet}}+T_{\text{outlet}})$, defined as average value of outlet temperature and inlet temperature (the black thick line) denotes the flame length. As the temperature rises sharply near this value, the mixture combusts rapidly. The result obtained by the Boivin mechanism is about 88 mm, and that obtained by the Sandia mechanism is about 135 mm. Obviously, the Boivin mechanism can predict a more reasonable flame length compared with that in the experiment. Measured by a flue gas analyzer, the CO concentration at the outlet in the experiment is about 8 ppm (corrected to 15% oxygen by volume). The simulation adopting the Boivin mechanism predicts 11 ppm, whereas that based on the Sandia mechanism predicts 1.5 ppm. The former obtains better results than the latter.

Figure 6. Time-averaged temperatures obtained by (a) Boivin mechanism and (b) Sandia mechanism.

Figure 7 shows radial profiles of mean axial velocity and mean heat release at four axial locations obtained by two mechanisms. As shown in Figure 7, for X= −0.38D, the mean axial velocity predicated by two mechanisms are exactly the same, since there is almost no reaction. For the other three cases, when R < 0.02 m, the axial velocity predicated by Boivin mechanism is lower compared with Sandia mechanism. When R > 0.02 m, the axial velocity predicated by Boivin mechanism is higher compared with Sandia mechanism, since heat release obtained by Boivin mechanism is higher.

Figure 7. Radial profiles of (a) mean axial velocity and (b) mean heat release at four axial locations obtained by two mechanisms.

The results imply that the Boivin mechanism can obtain more accurate simulation data compared with the experiment. Therefore, the following discussion is focused on the simulation results based on the Boivin mechanism, as well as the effect of hydrogen addition on the combustion characteristics in the chamber.

4.2. Flame structure

To investigate the effect of hydrogen addition ratio on combustion characteristics, we perform LES on three cases on the basis of the Boivin mechanism (Table 2).

Evolution of instantaneous concentration of H radical for the three cases based on the Boivin mechanism are shown in Figure 8. The thick black line in the figure indicates the mid-temperature. In Figure 8, the concentration of the H radical is the highest near the isotherm of the mid-temperature; this result indicates that the flame burns violently near the black line. As the isotherm can indicate the flame structure, the isotherm can be approximated as the flame front to study the differences across the three cases.

Figure 8. Evolution of instantaneous concentration of H radical for (a) Syngas25, (b) Syngas45, (c) Syngas100.

Syngas25 has the lowest hydrogen content, but its flame length is the longest, and its flame front extends radially. In addition, the high concentration of H radical downstream of the flame front covers a large area, and thus, the reaction rate is low. When the H₂ content increases from 20% to 45%, the flame shortens radially, and the flame front shrinks toward upstream. For Syngas100, the flame front continues to shrink toward upstream, and H radical are highly concentrated near the flame front. The area of high-concentration H radical downstream of the flame front decreases because the Le of hydrogen is lower than that of carbon monoxide and regions with high H radical concentrations caused by preferential diffusion appear.

Karlovitz number (Ka) also affects the distribution of H radical. Ka is defined as the ratio of the chemical reaction time scale to the turbulent Kolmogorov time scale. As the hydrogen combustion requires little time, the chemical reaction time scale becomes small, thus Ka decreases with hydrogen addition ratio increasing. With hydrogen decreasing, i.e. Ka increasing, the turbulent vortex structure could easily penetrate the flame surface and affect the internal flame structure, which will promote or suppress chemical reactions. According to the previous research performed by Wang et al. (2018), with increasing Ka, i.e. increasing turbulence intensity, the reactions will be promoted at weak reaction regions and suppressed at intense reaction regions, leading to the increase of H radical at weak reaction regions. As we can see form Figure 8, with hydrogen content increasing, affected by decreasing Ka and Le, H radical distributed narrowly, which indicates reaction zone is thinned.

With the increase in hydrogen addition ratio, the flame surface transforms smooth from wrinkly. The observation is consistent with the previous research performed by Wang and Abraham (2018); they concluded that with Ka increasing, the probability of forming highly curved flame surface becomes higher. The influence of the turbulent vortex structure on the flame surface decreases with hydrogen content increasing, eventually leading to the flame surface becoming smooth. The other characteristics of the flow field in the chamber will be detailed in the next section.

Figure 9 shows evolution of instantaneous temperatures for the three cases. For Syngas25, a high-temperature zone appears near the outlet of the combustion chamber, and the maximum mean temperature is 1,490 K. For Syngas45, the high-temperature zone moves upstream, and the maximum mean temperature rises to 1,510 K. For Syngas100, a high-temperature zone appears downstream of the flame front, and the maximum mean temperature is 1,570 K. The change in the location of the high-temperature zone is consistent with the prediction of flame structure in Figure 8. The observation is consistent with the main conclusion of Park et al. (2009). As the hydrogen fraction increases, the Le of the fuel decreases, leading to the appearance of a high-temperature zone caused by preferential diffusion. The result is similar to the previous research performed by Fu et al. (2014), they observed that preferential diffusion leads to the formation of a high temperature zone in the flow field.

Figure 9. Evolution of instantaneous temperature for (a) Syngas25, (b) Syngas45, (c) Syngas100.

4.3. Flow field and flashback

Figure 10(a) shows the radial profiles of the mean axial velocity at four axial locations across different cases. As shown in Figure 10(a), CTRZ appears in the downstream of the flame front near the centerline since the axial velocity in the region is negative. The CTRZ provides a dominant flame stabilization mechanism. Shear layer and turbulence affect the flow in CTRZ (Huang & Yang, 2009). Meanwhile, corner recirculation zone (CRZ) appears at the expansion corner of the chamber.

Figure 10. Radial profiles of mean and fluctuating axial velocity at four axial locations.

Although the hydrogen fraction in syngas varies greatly, the average axial velocity at each location presents the same trend. Compared with those in Syngas25 and Syngas45, the peaks of the mean velocity in Syngas100 occur at a larger radial location. The maximum velocity in Syngas100 is higher than that in the other cases due to the high combustion heat release shown in Figure 11(a). Figure 11(a) shows the radial profiles of the mean heat release rates at four locations for different cases.

Figure 11. Radial profiles of mean and fluctuating heat release rates at four axial locations.

In the combustion chamber, the change of axial velocity is affected by the two factors below: swirling motion and combustion heat release. The swirling motion diverges the flow and reduces axial velocity, whereas combustion heat release expands the gas volume and increases the axial flow velocity.

Figure 10(b) shows the radial profiles of the fluctuating axial velocity at four axial locations across different cases. The simulation results accurately capture the position of the inner and outer shear layers, which are two peaks of velocity pulsations in the radial direction. When X = 0.4D, 0.03 m < R < 0.04 m is the inner shear layer, and 0.06 m < R < 0.07 m is the outer shear layer. Large-scale coherent flow structures are generated in the shear layer regions and shed downstream sequentially due to Kelvin-Helmholtz instabilities (Huang & Yang, 2009). These vortex structures significantly influence the combustion process.

The combustion reaction can cause velocity pulsations and provide energy for the existence of the shear layer. As hydrogen has a high thermal conductivity, its reaction position is concentrated upstream of the flow field, and the flame length is short. Therefore, the effect of combustion on the shear layer downstream of the flow field is minimal. Relative to those of Syngas25 and Syngas45, the peak of the fluctuating velocity in Syngas100 is decreased.

The radial profiles of the fluctuating heat release rates at four locations for different cases are shown in Figure 11(b). With the increase in added hydrogen, the value change of the RMS heat release becomes sharp, and the range of the RMS heat release becomes narrow; hence, the flame surface becomes smooth. The result is consistent with the discussion on Figure 6.

Figure 12 shows the distribution of the mean axial velocity and streamline for every case. Several salient structures are observed in the flow field: vortex breakdown-induced CTRZ, PVC surrounding the CTRZ, shear layers (including inner shear layers and outer shear layers), and CRZ. Vortex breakdown usually manifests as an abrupt change in the core of a slender vortex. It provides the dominant flame stabilization mechanism in this region, which exists internal stagnation points and reversed flows.

Figure 12. Distribution of mean axial velocity and streamline for (a) Syngas25, (b) Syngas45, (c) Syngas100.

As we can see from Figure 8, flame kernel detaches from the shear layers intermittently, and the detached flame kernel becomes small with hydrogen content increasing. After detaching from flame, the flame kernel reaches the CRZ, and flame ignites intermittently, which could make the flame extend radially. Since Syngas25 has the larger flame kernel, the flame could extend radially, whereas Syngas100 hardly extends radially. Hence, CRZ shrinks and recirculation weakens with hydrogen content increasing as shown in Figure 12. The sudden and intermittent ignition in the CRZ leads to large fluctuations in heat release and temperature (as shown in Figure 9), which is to be considered as a potential mechanism leading to thermo-acoustic instability (Soufien Taamallah et al., 2016).

As hydrogen has a high thermal conductivity, a high hydrogen addition ratio leads to a low heat release rate downstream and a high heat release rate upstream. As the hydrogen addition ratio increases from 25 to 100, the size of the CTRZ decreases, the stagnation point of the vortex breakdown moves upstream and eventually approaches the center body, then flashback caused by CIVB occurs. Similar observations in CH₄-H₂ syngas have also been reported by De and Acharya (2012a), they found high hydrogen content in CH₄-H₂ syngas makes the CTRZ compact, leading to flashback.

CIVB is promoted by the high gravimetric heating value of hydrogen. As the flame moves upstream as the hydrogen fraction increases, it will encounter a fuel-rich zone, thereby resulting in hydrogen with a high flame propagation velocity. The basic mechanism of CIVB is that the flame contributes to vortex breakdown, leading to CTRZ and the stagnation point. As the flame moves upstream, the location of vortex breakdown advances further upstream. This process will continue as the flame moves further upstream.

The PVC is an unsteady flow structure that central vortex core precesses around the symmetry axis at a well-defined frequency. This phenomenon has often been reported in turbulent swirl burners (Froud et al., 1995; N Syred et al., 1997). The PVC is usually linked to vortex breakdown and the associated recirculation zone in a high Re flow.

Figure 13 shows evolution of instantaneous pressure over one cycle for three cases. The thick black line indicates the contour of zero axial velocity, and environmental pressure, p_e= 101,325 Pa. Here, time cycle of swirling, $t_{\text{cycle}}=(\pi d/w_{\text{ave}})\approx 10\text{ms}$, the tangential average velocity, $w_{ave}\approx$65 m/s.

Figure 13. Evolution of instantaneous pressure over one cycle for (a) Syngas25, (b) Syngas45, (c) Syngas100.

For Syngas25, we can clearly observe the existence of the PVC (low pressure core), the PVC is situated near the zero velocity line. This phenomenon is consistent with the observation reported by Syred and Beer (Nicholas Syred, 2006) that the PVC is usually situated on the boundary of the reverse flow zone between the zero velocity and the zero streamline. As we can see from Figure 13, the PVC occurs near the center body and precesses downstream.

For Syngas45, we can also observe the existence of the PVC (low pressure core), though the pressure of the PVC becomes higher relative to Syngas25. Since higher heat released by hydrogen makes pressure becomes higher. Similar to Syngas 25, the PVC starts occurs near the center body and precess downstream.

For Syngas100, we can observe the PVC form t = 2100 ms to t = 2106 ms, and the pressure of the PVC is higher relative to Syngas25. However, the PVC disappears form t = 2108 ms to t = 2110 ms, which means the PVC cannot exist stably in a time cycle.

As shown in Figure 9, when the hydrogen addition ratio increases, the flame temperature rises greatly. As the temperature rises, the viscosity of the gas improves greatly, and the turbulence, in particular, the turbulent stress increases significantly. These factors generally explain why the PVC cannot exist stably in the combustion (Huang & Yang, 2009).

5. Conclusions

In this work, the effects of hydrogen addition ratio for CO-H₂ syngas on combustion characteristics are numerically investigated in a swirl burner via LES. The main conclusions can be drawn as follows:

1. Compared with experimental data, the Boivin mechanism predicts a more reasonable flame structure.

2. As the H₂ addition increases, the flame temperature rises, and the flame length decreases due to preferential diffusion. Besides, H radical distributed narrowly, and the flame surface becomes smooth.

3. As hydrogen has a high thermal conductivity, the reaction location moves upstream as the hydrogen content increases. This condition leads to a high heat release and RMS of heat release. Meanwhile, the effect of combustion on the shear layer downstream decreases.

4. As the H₂ addition increases, the size of the CTRZ decreases, and flashback caused by CIVB occurs.

5. The PVC cannot exist stably in Syngas100, though it can be clearly observed in other cases. And the pressure of the PVC is higher in Syngas100 and Syngas45 relative to Syngas25.

Due to the limitation of computing resource, we just set three cases to investigate the effect of hydrogen addition ratio on combustion characteristics. In the future, more cases with various compositions should be set for further research.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work is Supported by Foundation of State Key Laboratory of High-efficiency Utilization of Coal and Green Chemical Engineering [grant number 2019KF19].