A Computational Fluid Dynamics Study of Turbulence, Radiation, and Combustion Models for Natural Gas Combustion Burner

This paper presents a Computational Fluid Dynamics (CFD) study of a natural gas combustion burner focusing on the effect of combustion, thermal radiation and turbulence models on the temperature and chemical species concentration fields. The combustion was modelled using the finite rate/eddy dissipation (FR/EDM) and partially premixed flame models. Detailed chemistry kinetics CHEMKIN GRI-MECH 3.0 consisting of 325 reactions was employed to model the methane combustion. Discrete ordinates (DO) and spherical harmonics (P1) model were employed to predict the thermal radiation. The gas absorption coefficient dependence on the wavelength is resolved by the weighted-sum-of-gray-gases model (WSGGM). Turbulence flow was simulated using Reynolds-averaged Navier-Stokes (RANS) based models. The findings showed that a combination of partially premixed flame, P1 and standard k-ε (SKE) gave the most accurate prediction with an average deviation of around 7.8% of combustion temperature and 15.5% for reactant composition (methane and oxygen). The results show the multi-step chemistry in the partially premixed model is more accurate than the two-step FR/EDM. Meanwhile, inclusion of thermal radiation has a minor effect on the heat transfer and species concentration. SKE turbulence model yielded better prediction compared to the realizable k-ε (RKE) and renormalized k-ε (RNG). The CFD simulation presented in this work may serve as a useful tool to evaluate a performance of a natural gas combustor. Copyright © 2018 BCREC Group. All rights reserved


Introduction
The combustion system involving natural gas or methane is commonly used in the industry, especially in the power generation plant.Assessment of the combustor performance and pollution generation can be performed experimentally using a combination of equipment such as online gas chromatography, arrays of temperature sensors and advanced nonintrusive laser-based measurement.However, the cost to set such an experimental rig is often prohibitively too high for most researchers.In addition, the combustor chamber wall must be made of a special material such transparent quartz to enable laser-based measurement.In addition, measurement at high temperature (>1000 K) is potentially dangerous [1][2].Alternatively, CFD can perform a detailed evaluation of the combustion system, but it need to be validated before it can be routinely used.
A simple cylindrical burner is often used as a test bed to evaluate the accuracy of the CFD modelling strategy of a natural gas combustion system.One of those experimental works in a cylindrical combustion chamber with detailed measurement of gas composition and temperature profiles was presented by Garreton and Simonin [3].A CFD study on the cylindrical burner has been performed by several researchers [4][5][6][7][8][9][10].They emphasized on the chemistry reaction, turbulence mixing, thermal radiation, pollution formation and buoyancy effect inside the burner.Most of the early work evaluates the effect of models used to the prediction of heat transfer and chemical species profile.The findings obtained from prior works concluded that the modelling approach is vital to the prediction accuracy in cylindrical burner.Hence, the aim of this work is to develop an accurate modelling strategy by evaluating the effect of different turbulence, combustion and radiation models to the temperature and gas species profile in the burner.
It is important to predict the turbulent flow inside the combustion chamber accurately to enable better prediction of the detailed reaction chemistry involved in the combustion process.RANS-based turbulence model, such as: the SKE model is commonly used to resolve the turbulent flow due to its robustness and lesser computational demand [4,[7][8].Ronchetti et al. [9] performed a comparison of different turbulence models on the prediction of temperature and carbon monoxide mass fraction.They found that no substantial difference was obtained between the k-ε and k-w turbulence models.SKE is known to provide a reasonable prediction on the temperature and chemical species concentration for natural gas combustion, although it has a known issue to maintain a positive turbulence stresses besides giving a poor prediction of rotational and strained turbulence flow.The newer k-ε variant, i.e.RKE and RNG, are known to address the aforementioned issues.However, no previous work deals with various k-ε based models, such as: SKE, RKE, and RNG for the methane combustion in a cylindrical burner.Hence, the effect of SKE, RKE, and RNG on the prediction of temperature and the gas species profile was assessed in this work.
Combustion models, such as: the eddy dissi-pation concept (EDC) [4], eddy break-up (EBU) [10], presumed probability-distributionfunction (PPDF) [10], finite rate/eddy dissipation model (FR/EDM) [8][9], have been widely employed to predict the natural gas combustion.Among all the aforementioned combustion models, EDM is the most commonly used due to its reasonable predictions for methane combustion [11].Therefore, a two-step EDM was considered in the present work.Thermal radiation dominates the heat transfer process in most combustion systems like the natural gas combustion burner.Prior work has shown that thermal radiation accounts for 96% of the total heat transfer in the combustion system [12].Earlier works by da Silva et al. [6] on the similar case focused on the effect of thermal radiation on the temperature and chemical species concentration distribution.Their work indicated that the inclusion of thermal radiation gave a more uniform and accurate heat transfer prediction inside the combustion chamber.Wang et al. [13] reported that the combustion simulation without radiation model over-predicts the temperature field.Hence, it is vital to consider the radiation model to the heat transfer model for natural gas combustion.Most of the CFD studies dealing with the natural gas combustion employed discrete transfer radiation model (DTRM) for radiation [4][5][7][8].It has to be noted that DTRM does not include the effect of radiation scattering and can only be accurate when a large number of rays is modelled (CPUintensive).In addition the reflection of incident radiation at the surface is isotropic with respect to the solid angle, which is questionable, since the radiation should be a function of solid angle.
All the aforementioned issues are addressed in the DO and P1 models.However, no previous work used DO and P1 models for the case studied in this work, therefore, DO and P1 models with gas absorption coefficient WSGGM were assessed in this work.In this work, the modelling strategy was developed by evaluating the effect of different turbulence, radiation and combustion models in a successive way.The prediction was validated with the experimental data from Garreton and Simonin [3].

Geometry and computational grid
The geometry in this work is similar to the one measured by Garréton and Simonin [3].A two-dimensional axisymmetric cylindrical combustion burner was prepared by using GAM-BIT 2.4.6 as shown in Figure 1.The cylindrical burner has 170 cm of length and 25 cm in radius.The burner consists of two ducts of inlet which is air and fuel inlet.The natural gas (0.232 m 3 /s) is injected into the burner from the fuel inlet with a radius of 3 cm, while the air (0.728 m 3 /s) enters the chamber through a centered annular duct having a spacing of 2 cm.The outflow of the chamber is 12.5 cm in radius.The whole domain was prepared by quadrilateral mesh.Four different grids, (i.e.140k, 335k, 560k and 650k), were tested in this work.

Combustion modelling
In this work, combustion of natural gas was modelled by FR/EDM [14] and partial premixed flame model.The two models chosen is suitable for a fast reaction (Damkohler »1) like the one in this work i.e.Da ≈ 2.86 (volume averaged); although at the flame region the Da can reach as high as Da ≈ 63.In EDM a fast chemical reaction was assumed to be controlled by the turbulent mixing, while FR model abandoned the effect of turbulent mixing and computed the chemical reaction rate according to the Arrhenius equation.The FR/EDM model switches automatically between the two mode using the data obtained from the CFD simulation i.e. data on temperature and turbulent flow is automatically fed to the FR/EDM model during CFD simulation.The simplified two-steps combustion of natural gas is given by Equations (1-2). (1) (2) The species transport equation is given by: (3) where ūj is the mass-averaged velocity of mixture, is the mass fraction, Di,m is a diffusion coefficient for species i in mixture, μt is the turbulent viscosity, Sct is the turbulent Schmidt number, Ri is the net production rate and Si is the source term.The net production rate is given as: (4) where Mw,i is the molecular weight of species, NR is the total number of reactions, and is the Arrhenius reaction rate.In EDM, the production rate of species is modelled according to Magnussen and Hjertager (Eqs.(5-6)) [14]: (5) (6) where and are the stoichiometric coefficients of reactant and product, respectively, A and B are the Magnussen constant for reactant and product, respectively, YP and YR are the mass fractions of the species in product and reactant, respectively.
Only two simplified kinetic mechanisms were used in EDM to solve the natural gas combustion reaction rate.However, the combustion of natural gas is complex due to the multiple chemical reactions that occur simultaneously with the turbulence and heat transfer.Hence, an inclusion of a multi-step reaction model is vital in order to get an accurate prediction.A multi-step mechanism was introduced into the flamelet library to account for Figure 1.Two-dimensional geometry of combustion burner ,r j v the turbulence and non-equilibrium chemistry.A detail CHEMKIN GRI-MECH 3.0 reaction mechanism [15] which consisted of over 325 reactions and 53 species equipped with associated rate and thermodynamic data was used for partially premixed turbulent combustion model.In this work, the combustion occurs in both the non-premixed and premixed mode.Initially, natural gas and air are introduced separately into the burner (non-premixed).The natural gas and air are partially premixed at the base of the lifted diffusion flame.The inhomogeneous turbulent mixing separates the mixture into fuel-rich and fuel-lean regions.As the flame front propagates, the thin flame sheet separates the regions into unburnt and burnt mixture regions.Therefore, partially premixed combustion was considered by combining the both flamelet models from non-premixed combustion and premixed combustion, respectively [16][17].
The mixing of natural gas and air in a turbulent flow field is described using the mixture fraction model.The transport equation of Favre mean and variance of mixture fraction is modelled in Equations (7-8).( 7) (8) where is mean mixture fraction and Sct is turbulent Schmidt number.In premixed flame model, the reactive flows divided into burnt and unburnt region, which is separated by the flame sheet.In premixed combustion model, a progress variable is used to model the flame front propagation (Equation ( 9)): (9) where c is mean reaction progress variable, μt is turbulent viscosity, Sct is the turbulent Schmidt number, ρu is the density of unburnt mixture and Ut is the turbulent flame speed.The progress variable is computed as in Equation (10): (10) where n is the number of products, Yi is the mass fraction of product species, and Yi,eq is the equilibrium mass fraction of product species.It is given that c = 0 for unburnt mixture and c = 1 for burnt mixture.In partially premixed flame model, the reaction behind the flame front is modelled by mixture faction model and the flame front position is determined using progress variable.It is best suited for a fast reaction (i.e.Da » 1) especially for the case of chemical equilibrium or moderately nonequilibrium flamelet structure.

Radiation modelling
As discussed earlier radiative heat transfer account for about 96% of the total heat transfer in the combustion system and hence must be modelled accordingly.The P1 model is the simplest radiation model generate by the P-N model which is based on the expansion radiation intensity into an orthogonal spherical harmonic [18].Only zeroth and first order moments of the intensity are considered in the P1 model.P1 model solves isotropic radiative heat transfer and it requires low computational demand [19].Simulation via P1 model accounts for scattering effect and it is applicable for large optical thickness.The transport equation for P1 is modelled as [20] (Equation ( 11)): (11) where qr is radiative heat flux, a is the absorption coefficient, G is an incident radiation flux, σ is the Stefan-Bolzmann constant, and T is a temperature.The radiative heat flux at wall is given in Equation ( 12): (12) where σs is the scattering coefficient and C is a linear-anisotropic phase function coefficient (between -1 and 1).A backward scattering is given at a negative value, a forward scattering is in positive value and a zero value is denoted for an isotropic scattering.
The DO model utilizes a different approach to solve the radiation transfer equation (RTE) compared to P1 model.The solid angle at a certain point of domain is split up into a number of discrete directions and the radiative intensity is assumed to be constant within each division of the solid angle.DO model is more time consuming than the P1 model due to the solution is required for many different directions [19].The RTE is modelled in Equation ( 13): where Iλ is the spectral radiation intensity, λ is a wavelength, aλ is a spectral absorption coefficient, s is the path length, Ibλ is a black body intensity, n is a refractive index, ϕ is a phase function, is a scattering direction and dΩ' is a solid angle.

Turbulence modelling
Fluid flow in a combustion process is usually turbulent whereby the velocity and pressure fluctuate chaotically.Turbulent flows can affect the heat transfer and chemical reaction in the combustion process, and hence must be included in the CFD model.The accuracy and reliability of a CFD simulation is significantly depends on the model used.Miltner et al. [21] and Ilbas et al. [22] reported that no single turbulence model can be universally applied in all cases.Therefore, three RANS turbulence models, namely: SKE, RKE, and RNG, were compared in this work.The RANS transport equations are given in Equations (14-15).( 14) (15) where ūi is mean velocity, ρ is fluid density, is external forces, is mean pressure, μ is fluid viscosity, is mean strain tensor rate, and is Reynolds stresses tensor.The SKE is the most used k-ε turbulence model as it is easier to converge and requires relatively low computational demand.The turbulent kinetic energy equation for SKE is modelled in Equation ( 16).(16) where ρ is the fluid density, k is turbulent kinetic energy, μ is fluid viscosity, μt is turbulent viscosity, σk = 1.0 is Prandtl-Schmidtl number, Gk and Gb are the production rate due to mean velocity gradient and buoyancy, respectively, ε is dissipation rate, and YM is the dilatation dissipation term accounts for compressibility effect.The turbulent viscosity is computed by Equation (17).(17) where Cμ is given as 0.09.The production rate of SKE is given in Equations (18)(19).
(18) (19) where is normal stresses, Prt is turbulent Prandtl number with a constant value of 0.85, and gi is component of gravitational vector in ith direction.The destruction rate (turbulent dissipation rate) is given in Equation ( 20).(20) The transport equation of dissipation rate in SKE is modelled in Equations (21)(22).( 21) (22) where v and u are the component of velocity parallel and perpendicular to gravitational vector, respectively.The model constants are σε = 1.3, C1ε = 1.44, and C2ε = 1.92 [23].
The SKE model often gave a poor prediction of flow with a strong streamline curvature, gradient flow and rotation owing to its constant eddy viscosity formulation which can lead to a negative normal stresses computation under certain circumstances.It was known to give a poor prediction on the species concentration owing to the constant value of Cε in the transport equation of dissipation rate [24].Hence, other k-ε variant such as RKE [25] and RNG [26] were introduced to overcome the limitation of SKE.RKE differs from the SKE because its turbulent viscosity is no longer constant.The turbulent viscosity coefficient of RKE is computed as a function of local states of the flow to ensure a positive normal stresses ( ) under all flow conditions.Therefore, this model can provide a better prediction of the rotation, vortices and separation flows features [27].The Cμ for RKE is computed by Equation ( 23):

Modelling setup
The simulation of natural gas combustion in a two-dimensional cylindrical chamber was performed using ANSYS FLUENT 16.2 installed on the HP Compaq Pro 6300 MT workstation with a Quadcore i7-3770 processor (3.40 GHz) and 4 GB RAM.The simulation was firstly performed using first-order upwind scheme, steady-state SKE turbulence, DO radiation, and FR/EDM.The unsteady-state solver and higher-order discretization scheme was then enabled after a converged solution was achieved.The thermophysical properties (i.e., specific heat, dynamic viscosity and thermal conductivity) of each chemical species at temperature range from 300 to 2500 K were introduced as a piecewise linear function.In partially premixed flame model, the GRI-MECH 3.0 associated with 325 mechanisms was used for more detail prediction.NASA polynomials (Thermochemical Data for Combustion Calculations) were used to model the gas properties as a function of temperature.The data were recorded for over 1000 time steps after a pseudosteady solution was achieved and the value reported in this work is a statistical timeaveraged.The simulation setup used in this work is shown in Table 1.The CFD predictions from various model combinations (i.e., turbulence, radiation and combustion models) were compared with the experimentally measured temperature and chemical species concentration [3].

Grid density analysis
A two-dimensional axisymmetric burner was prepared by quadrilateral meshes.Four different grids of combustion burner (i.e., 140k, 335k, 560k and 650k) were tested in this work.The predictions for the four different grids were compared with the experimental data [3].
Figure 2 shows the temperature distribution at different position of combustion burner for the four different grids.In Figure 2(A), all grid types yielded similar predictions of temperature at 0 m < X < 0.9 m along the centre of the burner.However, the two coarser grids, namely 140k and 335k grids under-predicts the temperature at 0.9 m < X < 1.7 m, whereas, the two finer grids (560k and 650k) showed a fair prediction.Figures 2(B) to 2(D) show the temperature profiles along the radial position at three different axial positions.No substantial difference of predictions obtained from the four different grids in Figure 2(B).However, Figures 2(C) and 2(D) clearly showed that the finer grids (560k and 650k) yielded a better prediction compared to the two coarser grids.Figure 3 showed that the 560k and 650k grids provided more accurate predictions on the chemical species concentration, except for the Figure 3(D).Hence, the two finer grids are the suitable choice in this work.However, the higher grid densities need longer computational time, as shown in Table 2.The two coarser grids (i.e.140k and 335k) gave a poor prediction, although they require much lesser computational time.The 560k grid was selected, instead of the finest grid (650k) for the rest of this work to minimize the computational demand, since the predictions obtained by the 560k grid is comparable to the one by 650k grid.

Effect of radiation model
This work aims to investigate the importance of including radiation model in the natural gas combustion simulation.Therefore, the simulation with P1 radiation model was compared with the one without radiation model and without weighted-sum-of-gray-gases model (WSGGM).The predictions of temperature and   4 and 5, respectively.It was clearly shown that the simulation is in better agreement with the experimental data [3] when the radiation is included and WSGGM is enabled.Whereas, the simulation without the radiation model enabled produce a relatively poor prediction of temperature and chemical species concentration.This is because the inclusion of radiative heat transfer enhanced the homogenization of temperature inside the combustion burner by transferring the thermal heat from hot gas region to burner's wall and outlet [6].Therefore, the temperature becomes more uniform, unlike the one without a radiation model.For instance, Figures 4(A-C) show that the temperature is over predicted in the core region, but under-predicted near the outlet in Figure 4(D) due to lack of radiative heat transfer homogenization.
Inclusion of radiation model without the WSGGM to account for the gray gas absorption coefficient also produces a less accurate prediction on temperature and species profile.The absorption of gas species in combustion chamber is not constant, but is depends on the temperature.WSGGM was introduced to resolve the spectral gas absorption, and therefore it is important to be included for the radiative heat transfer like the simulation in this work.This work showed that the inclusion of radiation provided more accurate prediction of temperature and chemical species, and hence radiation model with WSGGM must be used.
The simulation using P1 radiation model was then compared with the one using the DO model.Figures 6 and 7 are the temperature and chemical species concentration profiles, respectively, using two different radiation models (i.e.DO and P1).Partially premixed flame and SKE models were employed.The CFD prediction in this work shows a reasonable agreement with the experimental data from Garreton and Simonin [3].Although both DO and P1 models over-predicts the temperature along the radial position at 0.912 m from the inlet of the burner as shown in Figure 6(C) and had a relatively poor prediction of carbon monoxide concentration as shown in Figure 7(D).The large deviation of the predicted temperature at radial position of X = 0.912 m is a follow through of poor gas fraction prediction in the same region (see Figure 7 at X ~ 0.9).The radiation through the gas inside the chamber is modelled using a WSGGM.The WSGGM uses a number of grey gases and weighting factor polynomials to model gas radiative properties, i.e. emissivity.Thus error in gas composition prediction may affect the radiation heat transfer rate, since radiation account for about 90% of heat   transfer in a combustion chamber.In addition, P1 is known to over-predict the radiative fluxes from localized heat sources, i.e. combustion flame.It was found that the simulation using DO and P1 models yielded a similar trend on the temperature and chemical species concentration profiles at various positions of the burner.
Theoretically, DO solves a finite number of discrete solid angles and hence DO requires higher computational demand than the P1 model.In the simplified two-dimensional case like the one presented in this work, the difference between the P1 and DO models is not pronounced.This is attributed by the limited solid angle available for radiation in twodimensional simulation.It is worth noting that the DO is more accurate than P1 for 3D simulation like the one presented in our previous work [2].Therefore, the P1 model was used for the remainder of this work to provide a quick estimation of heat transfer in the combustion burner.

Effect of turbulence model
The effect of turbulence model on the prediction of temperature and chemical species concentration inside the burner was evaluated using an unsteady-state, P1 and partially pre-mixed combustion models.The three different RANS turbulence models are SKE, RKE. and RNG.The predictions were taken along the symmetry line and the radial position of the burner.
Figures 8 to 9 shows the comparison of the predicted temperature and chemical species mass fraction inside the burner obtained using the three turbulence models with the experimental data [3].The results clearly showed that the simulation via SKE yielded the best agreement with the experimental data [3], whereas the RKE gave the poorest prediction.This may be attributed by the fact that the fluid flow inside the burner is mostly isotropic and homogenous turbulence, which favours SKE.The fluid mixing in the burner did not feature a strong swirling flow, which is suited the RKE and RNG as shown in Figures 10 and  11.Only a minor recirculating flow appeared in the region just above the inlet and outlet, respectively (refer Figure 11).Therefore, the SKE model is sufficient for the natural gas combustion modelling in a cylindrical burner.

Effect of combustion model
The SKE turbulence model was then used to evaluate the effect of combustion models on the temperature and chemical species concentra-  shows a relatively poor prediction of the temperature at various positions of the burner.In Figures 13(A) and (C), the mass fraction of methane and carbon dioxide are well resolved by both FR/EDM and partially premixed flame models.
Figure 13(B) shows the prediction using partially premixed flame model yielded a better agreement with the experimental data [3], whereas the FR/EDM over-predicts the oxygen mass fraction.This is attributed by the detailed turbulent flame modelling of partially premixed flame model by combining the modelling strategies for non-premixed and premixed flame models.The partially premixed flame model considers the turbulent mixing for both fuels-rich and fuel-lean regions and also determines the reactions for both burnt and unburn mixture regions.In addition, GRI-MECH 3.0 is optimized for the turbulence-chemistry in methane oxidation at the wide range in temperature, like the natural gas combustion in the present work.GRI-MECH 3.0 considers 325 reaction mechanisms, including the intermediate reactions and chemical species dissociation.However, FR/EDM only considers two-steps global mechanism.Therefore, the partially premixed flame with multiple mechanisms is more accurate compared to the FR/EDM.

Conclusions
The effect of modelling methods on the natural gas combustion in a two-dimensional cylindrical burner was performed by evaluating various turbulence, radiation and combustion models.The CFD simulation was successfully validated with the experimentally measured data.Application of the radiation model (i.e.DO and P1) in conjunction with the gas absorption coefficient model, WSGGM, improved markedly the prediction accuracy of radiation dominated the heat transfer in natural gas combustion burner.It was found that detailed mechanism GRI-MECH 3.0 provided more accurate prediction of the temperature and chemical species concentration (i.e.error of 7.8% and 15.5%) than a two-step FR/EDM (i.e.error of 17.5% and 31.4%).The findings obtained from this work showed that partially premixed combustion model coupled with the P1 radiation model and the SKE turbulence model is the best combination for the modelling of natural gas combustion in the burner.The CFD simulation presented in this work may serve as a useful tool to evaluate a performance of a natural gas combustor.

Figure 2 .
Figure 2. Grid density analysis on temperature distribution along the (A) symmetry line; (B) radial position at 0.312 m; (C) radial position at 0.912 m; (D) radial position at 1.312 m

Figure 3 .
Figure 3. Grid density analysis on chemical species concentrations along the symmetry line: (A) methane; (B) oxygen; (C) carbon dioxide; (D) carbon monoxide

Figure 4 .Figure 5 .
Figure 4. Comparison between with radiation model, without radiation and without WSGGM on temperature distribution along the (A) symmetry line; (B) radial position at 0.312 m; (C) radial position at 0.912 m; (D) radial position at 1.312 m

Figure 6 .Figure 7 .
Figure 6.Comparison of radiation model on temperature distribution along the (A) symmetry line; (B) radial position at 0.312 m; (C) radial position at 0.912 m; (D) radial position at 1.312 m

Figure 8 .
Figure 8.Comparison of turbulence model on temperature distribution along the (A) symmetry line; (B) radial position at 0.312 m; (C) radial position at 0.912 m; (D) radial position at 1.312 m

Table 1 .
CFD setup SKE turbulence model, DO radiation model, and partially premixed combustion model were employed for the grid densities comparison.

Table 2 .
Grid density analysis