Simulation of COH 2-Air Turbulent Nonpremixed Flame Using the Eddy Dissipation ConceptModel with Lookup Table Approach

We present a new combustion simulation technique based on a lookup table approach. In the proposed technique, a flow solver extracts the reaction rates from the look-up table using the mixture fraction, progress variable, and reaction time. Look-up table building and combustion simulation are carried out simultaneously. The reaction rates of the chemical species are recorded in the look-up table according to the mixture fraction, progress variable, and time scale of the reaction. Once the reaction rates are recorded, a direct integration to solve the chemical equations becomes unnecessary; thus, the time for computing the reaction rates is shortened. The proposed technique is applied to an eddy dissipation concept (EDC) model and it is validated through a simulation of a CO-H2-air nonpremixed flame. The results obtained by using the proposed technique are compared with experimental and computational data obtained by using the EDC model with direct integration. Good agreement between our method and the EDC model and the experimental data was found. Moreover, the computation time for the proposed technique is approximately 99.2% lower than that of the EDC model with direct integration.


Introduction
The growing availability of computational resources allows intensive use of numerical method to predict the chemical structure of flames and to design combustion equipment.However, simulation of turbulent combustion that includes detailed chemical mechanisms still remains elusive, although several techniques have been proposed to overcome this problem.For example, the intrinsic low dimensional manifolds (ILDM) [1,2] is a method to automatically reduce the detailed chemical mechanisms.This method is based on a direct mathematical analysis of the dynamic behavior of the responses of nonlinear chemical equations.However, the ILDM is not well agreed in the low-temperature domain because the dimension and complexity of a lookup table increase tremendously in this domain.A flamelet-prolonged ILDM (FPI) [3] is developed to overcome this problem.
Also, in the ILDM, diffusion is neglected in the construction of the lookup table, and therefore the results in the region where both chemistry and diffusion are significant are less accurate.To improve the accuracy in the region, the phase space ILDM [4] is proposed.
Another alternative is the in-situ adaptive tabulation (ISAT) [5][6][7] method, which is based on in-situ generation of a lookup table that is constructed by solving for the time evolution of species concentrations directly.This method has been adopted in ANSYS FLUENT [8] by using the eddy dissipation concept (EDC) model [9][10][11], and the probability density function (PDF) model [12,13].
Combustion simulation by an artificial neural network (ANN) [14,15] is a recent alternative.This method consists in training the ANN off-line using a thermo-chemical database that is representative of the actual composition and turbulence.Then, a stiff ordinary differential equation (ODE) solver for computing the reaction rates of the chemical species is replaced using the database of the ANN.
A laminar flamelet model [16][17][18] is widely used for combustion simulation.Although this is a numerical model for computing the mass fractions of the chemical species in turbulent combustion, the computation time for performing a simulation by the laminar flamelet model is very short.This is because, in this model, the temperature and mass fractions of the chemical species are obtained by a lookup table that is built before the combustion simulation.
A flamelet-generated manifolds (FGM) technique [19] shares the idea of the laminar flamelet model and can take into account both the chemical and transport processes as a function of the mixture fraction and progress variable using the lookup table.This technique has been applied to the simulations of the H 2 combustion [20], CH 4 combustion [21], and NO X emission [22].
A combustion simulation that considers detailed chemical mechanisms requires significant computation time because the reaction calculation involves n-dimensional ordinary differential equations that must be solved according to the number of chemical species present.If the computation time would be easily reduced without losing accuracy, we could be able to obtain results more quickly.
To reduce the computation time, a combustion simulation technique based on combination of the chemical equilibrium method and the EDC model (CE-EDC) has been developed in a previous study.The technique was validated by simulating an H 2 -air turbulent nonpremixed flame [23], a CO-H 2 -air turbulent nonpremixed flame [24], and a CH 4 -air partial turbulent nonpremixed flame [25].An advantage of the aforementioned technique is the ease with which a reduced chemical mechanism can be built according to the accuracy requirements for the chemical species.Thus, the technique can predict intermediate species with high accuracy even when a reduced mechanism for the minor species is built.In this technique, a lower number of species, modified by the chemical equations, leads to lower computation time.However, the prediction accuracy also depends on the number of modified species included.
In this work, we present a new simulation technique based on the lookup table approach.The table construction and combustion simulation are performed simultaneously.The reaction rates are estimated by a direct integration solution to the chemical equations and then storing them in the lookup table.Once the reaction rates are recorded in the lookup table, further direct integration of the chemical equations becomes unnecessary, and hence, the computation time is shortened without compromising on accuracy.The advantage of our technique is that the construction of the lookup table is not needed prior to the combustion simulation.Moreover, this technique is easily applied to other combustion models if the proper variables are chosen to the indexes of the lookup table.Although, reduced mechanisms are sufficient to obtain the temperature and the mass fractions of the major species, our technique can reduce the time for computing the ordinary differential equations (ODEs) of the reduced mechanisms.However, the mass fractions of the intermediate and radical species are needed according to the purpose of the simulation.For example, the mass fractions of O, OH, and H are needed to estimate NO X emission.In such case, the computation time increases with an increase of the number of the chemical equations.
The proposed technique was applied to the EDC model and was validated by computing a CO-H 2 -air turbulent nonpremixed flame.The results were compared with experimental data reported by Correa and Gulati [26] as well as with computational data obtained by using the EDC model with direct integration.

Method
2.1.Numerical Implementation.ANSYS FLUENT 13.0 was used in this study.The temperature, enthalpy, and reaction rate of the chemical species were calculated by a userdefined function (UDF).The UDF is a function that can be programmed to load the solver dynamically to enhance the standard features of ANSYS FLUENT.Moreover, the momentum equation, turbulence model, and mass fractions equations were solved by using ANSYS FLUENT.To solve the momentum equation, the velocity and pressure were coupled by using the SIMPLE method [27], and the basic Reynolds stress model in ANSYS FLUENT [28] was chosen as the turbulence model.The simplified transport model proposed by Smooke [29] was used to compute the viscosity, thermal conductivity, and mass diffusivity.Unfortunately, the energy equation in ANSYS FLUENT could not be used in this study because the reaction rates were computed by the user-defined function.The temperature was obtained by solving the following energy equation: where ρ denotes the density; u i , the velocity in the i direction; h, the enthalpy; x i , the coordinate in the i direction; μ t , the turbulent viscosity; σ h , the effective Prandtl number; λ, the thermal conductivity.The symbols -and ∼ on the variables represent the time and density average, respectively.T is the temperature computed as where T 0 = 298.15.Δ f h 0 J denotes the standard enthalpy of formation of chemical species J; Y J , the mass fraction of chemical species J; Cp mJ , the mean specific heat at constant pressure of chemical species J.
Cp mJ was computed as where Cp J denotes the specific heat at constant pressure of chemical species J.

Theory of the EDC Model.
In this study, the EDC model [30,31] was used as a combustion model.In this model, it is assumed that combustion occurs in the fine structure of the turbulent flow.The fine structure is considered a perfectly stirred reactor and the steady state of the following equation is solved under the initial conditions obtained from the current mass fractions, density, and temperature in each computational cell: where Y J is the mass fraction and Y * J and Y 0 J are the mass fractions in the fine structure and surrounding gas, respectively.t is the time, and w * J is the reaction rates in the fine structure.ξ * denotes the length of the fine structures and is expressed as where C ξ = 2.1377, μ denotes the viscosity; is the turbulent dissipation rate, and k denotes the turbulent kinetic energy.τ * denotes the fine structure region residence time, and is given as where C τ = 0.408.The mean reaction rate of chemical species J is computed as

Proposed Technique.
A lookup table was used to reduce the computation time.The flow solver (ANSYS FLUENT as mentioned above) extracted the reaction rates from the lookup table by using the mixture fraction, progress variable, and time scale of the reaction.The proposed technique was applied to the EDC model in this study.However, our technique is easily applied to other combustion models if the proper variables were chosen to the indexes of the lookup table .The progress variable and time scale of the reaction are defined by the mass fraction of CO 2 and the fine structure region residence time of the EDC model, respectively.
To use the proposed technique, the characteristics of the flame must be organized.The mixture fraction is generally used to model the turbulent diffusion flame.The EDC model is not based on the assumption of fast chemistry, and therefore the progress variable is needed to express how the reaction proceeds.Although, the mixture fraction and progress variable can almost express the characteristics of the turbulent diffusion flame, the fine structure region residence time of the EDC model was taken into account to improve the prediction accuracy in this study.The mixture fraction and progress variable were uniformly divided in the lookup table.However, the fine structure region residence time was  Stored point (white dots) Sampled point (black dot) nonuniformly divided since the time scale ranges from 1.0e-6 to 1.0 in this study.The lookup table construction and combustion simulation were carried out simultaneously.The reaction rates of the chemical species were recorded in the lookup table according to the mixture fraction, progress variable, and time scale of the reaction.The computation time was reduced with an increase in the number of the data of the reaction rates in the lookup table because solving the ODEs is unnecessary.
Figures 1 and 2 show the algorithm of the proposed technique and schematic representation of the lookup table construction, respectively.For simplicity, only a twodimensional lookup table is illustrated in Figure 2. If the variables, that is, the mixture fraction and progress variable, are at the sampled point P, then the reaction rates of the chemical species are stored at the points Q n,m , Q n+1,m , Q n,m+1 , and Q n+1,m+1 in the lookup table at one time.These points are used for the interpolation of the reaction rates.However, the stored reaction rates include an error because the coordinates of the stored and sampled points are slightly different.Therefore, to reduce this difference, the reaction rates can be averaged by using the following equation at each stored point: where W Q,new and W Q denote the new and old data of the reaction rates of the chemical species at the stored points , and Q n+1,m+1 , respectively.W P denotes the reaction rates at the sampled point P. W P is calculated as where Y is the mass fractions of chemical species, and Y * denotes the mass fractions of chemical species in the fine structure.Thus, the mean reaction rates are computed as where w J,interp denotes the reaction rates of chemical species J obtained by the interpolation.f 1,l and f 0 are the lth weight function and sum of the weight function, which are computed as where f 1,l is converted to an integer and varies from 1 to 8. In the above equations, σ = 0.05, r is the distance between the sampled point P and the stored points , and Q n+1,m+1 , N is the number of the writing reaction rates.
The weight function f 1,l is shown in Figure 3. r = 0 indicates that coordinates of the stored and sampled points are equal, which means the data obtained at r = 0 is favorable for the lookup table.f 0 expresses the coefficient of the old data of the reaction rates, while f 1,l implies the coefficient of the new data of the reaction rates.Thus, the data is well updated when r = 0. Once the reaction rates are written in the stored points Q n,m , Q n+1,m , Q n,m+1 , and Q n+1,m+1 , a linear interpolation can be performed to extract the reaction rates if the mixture fraction and progress variable are in the sampling range.
The stored data are sometimes updated randomly in the points that have larger error than a reference value (1.0e-4 in this study) to take the latest data into the lookup table.
Thus, the lookup table can be automatically updated if a new condition such as an artificial ignition is given.The error determines where the reaction rates should be updated.Generally, the large error is appeared in the outer of the flame, and the small error is estimated in the inside of the flame.The data in the outer of the flame is not updated because the changes of the reaction rates are so small in this area.A possibility to update the reaction rates at the stored points was computed by the pseudo random number and was set as 1/1000 in this study.This value can determine the theoretical maximum speed of the proposed technique.The frequency of taking the latest data into the lookup table increases if this value is set to the bigger value.The error is evaluated as where L 1 denotes the L 1 norm; err is the error value; W J,P is the reaction rates of chemical species J at the sampled point P; w J,interp is the reaction rates obtained by the interpolation; and F is the scale factor chosen from the maximum L 1 .
The mixture fraction is computed as where Z denotes the mixture fraction, and σ m is the effective Schmidt number.In this study, the progress variable c is expressed as where Y CO2 denotes the mass fraction of CO 2 .

Numerical Conditions.
A CO-H 2 -air turbulent nonpremixed flame was simulated to verify the accuracy of the proposed technique.The numerical conditions were based on Correa's experiment [26].Figure 4 shows a sketch of the computational domain used in this study.The mass fractions of CO, H 2 , and N 2 in the fuel were 0.393, 0.0332, and 0.574, respectively.In the ambient air, the mass fractions of O 2 and N 2 were 0.233 and 0.767, respectively.The fuel was injected from a nozzle with a radius of 3.18 mm.The velocities of the fuel and ambient air were 80.0 and 6.5 m/s, respectively.The temperature of the fuel and ambient air was considered to be 300.0K.A second-order upwind difference was applied to the governing equations of the momentum, mass fraction of chemical species J, and the energy equation.The 2D axisymmetric solver in ANSYS FLUENT was chosen using the axis indicated by the centerline of Figure 4.The Reynolds stress model was chosen as the turbulence model; the effective Prandtl number was σ h = 1.0 and the Schmidt number was σ m = 1.0.The model parameters were set to adjust the radial profile of the mixture fraction of the experimental data.In this study, the detailed chemical mechanism suggested by Maas and Warnatz [32]

Influence of the Total Number of Cells.
First, we verified the influence of the total number of cells in the computational domain on the numerical error, since it was needed to decrease this error as much as possible to increase the prediction accuracy of the proposed technique.However, increasing the number of cells to decrease the numerical error caused an increase in the computation time.Two different grids with approximately 140,000 and 330,000 cells were prepared to examine the influence of the number of cells.The cells were concentrated near the outlet of the nozzle and outside of the nozzle because of high gradient of the velocity.In the computational grid with 140,000 cells, the dimensions of first layer at the outlet and outside of the nozzle are 2.5e-5 m and 1e-5 m, respectively.In the computational grid with 330,000 cells, the dimensions of first layer at the outlet and outside of the nozzle are set to 1.6e-5 m and 6.7e-6 m, respectively.Figures 5 and 6 show plots of the radial mixture fractions predicted with the two different grids, where r denotes the radial coordinate and d denotes the diameter of the nozzle.The predicted radial mixture fractions were in good agreement between two grids.The discrepancy among the two computational grids is minimal and can be neglected.Therefore, the grid with 140,000 cells was chosen for the simulations discussed in the following sections.

Comparison with Computational Data Obtained by Direct
Integration.Figures 7,8,9,10,11,and 12 show the mass fractions computed by the EDC model (DI, direct integration) and that with the lookup table (proposed).The size (mixture fraction × progress variable × time scale) of the lookup table was 100 × 100 × 100.In all the figures, x denotes the axial distance.As shown in Figures 7-12, the results obtained by the proposed technique are in good agreement with those obtained by using the EDC model with direct integration.It is noteworthy that the mass fractions of OH, which is the radical species, were similar in both cases.

Comparison with Experimental Data and Computational Data by Laminar Flamelet Model.
In this section, to validate our technique, we compare the results obtained using our technique with the experimental data and the results obtained by the laminar flamelet since the laminar flamelet model is well used in the simulation of turbulent diffusion flames.In Figures 13-21, EXP and LF indicate the experimental data and computational data obtained by the laminar flamelet model in ANSYS FLUENT.
As can be seen in Figure 13, the axial mixture fractions obtained by the proposed technique and laminar flamelet model are slightly lower than the experimental data when 20 < x/d.It can be noted from Figures 14, 15, and 16 that the mixture fraction and mass fractions of CO and H 2 obtained by using our technique are close to the experimental data and computational data obtained by the laminar flamelet model.Note from Figure 17 that the peak of the mass fraction of H 2 O is predicted well in spite of the fact that the radial profile of the mass fraction of H 2 O is slightly different when r/d < 0.5 and when 1.5 < r/d < 2.5.The peak of the mass fraction of H 2 O obtained by the laminar flamelet model is smaller than the experimental data.In Figure 18, the mixture fractions obtained by using our technique and the laminar flamelet model are slightly different from the experimental data but the trend of the radial profile is computed well.Figures 19 and 20 show that the radial profiles of the mass fractions of CO and H 2 obtained by using our technique are slightly underestimated relative to the experimental data when r/d < 0.5 but are close when 0.5 < r/d.However, the mass fraction of H 2 obtained by the laminar flamelet model is better than that by our technique.The radial mass fraction of H 2 O computed by our technique is similar to the experimental data, while the laminar flamelet model underestimates the peak of the mass fraction of H 2 O.
We can see from Figures 13-21 that the mixture fraction and mass fraction obtained with our technique exhibited slight differences in the experimental data.The reason is that our technique only reduces the computation time for simulating the reaction rates of the chemical species.The prediction accuracy depends on the combustion and turbulence models.In general, the combustion simulation using a large eddy simulation (LES) model shows good agreement with experimental data [33,34].To minimize the difference between the prediction results and experimental data, the LES model should be used.

Comparison of Computation Time.
Figure 22 shows a comparison of the computation time for the EDC model with direct integration, proposed technique, and CE-EDC [23,24].The CPU used in these simulations was an Intel Xeon processor X5492 with 12 M cache, 3.4 GHz, and 1600 MHz FSB.GFortran 4.5.1 [35] was chosen as a compiler and CVODE 2.6 [36] was adopted to solve the ODEs of the reaction calculation.As can be seen in Figure 22, the computation times for the EDC model with direct integration, for the CE-EDC, and for the proposed technique were 50.0 s, 7.3 s, and 0.34 s, respectively.Thus, the computation time of our technique is approximately 99.3% lower than that of the EDC model, and the proposed technique is faster than the CE-EDC.In addition, the prediction accuracy of the proposed technique is similar to that of the EDC model with direct integration.Figure 23 shows the computation time of flow solver per iteration but not the time for solving the reaction equations since we cannot estimate time for solving the reaction equations in the flamelet model and EDC model using ISAT technique using ANSYS FLUENT.
As can be seen in Figure 23, the computation times for the flow solver using the EDC model with direct integration, using the laminar flamelet model, using the EDC model with ISAT technique in ANSYS FLUENT, using the CE-EDC, and using our technique are 53.2 s, 1.24 s, 6.83 s, 12.16 s, and 4.88 s, respectively.The computation time of the flow solver using the laminar flamelet model is lower than that using

Table and Memory
Size.In our technique, the amount memory required depends on the size of the lookup table, since the reaction rates are stored in the latter.The size of the memory used in our technique is roughly estimated as the size of memory for a variable (here, 8 bytes for double precision) times the number of chemical species times the size of the lookup table.In this study, 15 species were included; three different sizes of the lookup table were prepared to examine the influence of the size of the lookup table.The sizes of the lookup tables were 50 × 50 × 50 (15 Mbyte), 100 × 100 × 1 (1.2 Mbyte), and 200 × 200 × 200 (960 Mbyte), respectively.The table with 100 × 100 × 1 was prepared to validate the prediction accuracy using the two-dimensional lookup table.As it can be seen in Figure 24, the mass fraction of OH obtained using the lookup table with 50 × 50 × 50 and that with 200 × 200 × 200 are in good agreement, while that obtained using the lookup table with 100 × 100 × 1 is different.Note from Figure 25 that the discrepancy is minimal, although both peak values of OH exhibit differences at r/d = 2.0.However, the result obtained using the lookup table with 100×100×1 is different from the mass fraction of OH obtained using the look-up table with 50 × 50 × 50 and 200 × 200 × 200.Thus, if there is insufficient memory in the computer, the lookup table with 50 × 50 × 50 can be used, but the reaction time scale cannot be eliminated.

Conclusions
In   and 200 × 200 × 200 (960 Mbyte) were tested; it was found that if there was insufficient memory in the computer, the lookup table with 50×50×50 could be used, but the reaction time scale cannot be eliminated.

Figure 1 :
Figure 1: Algorithm of our technique.

Figure 2 :
Figure 2: Schematic representation of the building lookup table.

Figure 3 :Figure 4 :
Figure 3: Plot of the weight function.

Figure 23 :Figure 24 :
Figure 23: Comparison of computation times for flow solver per iteration.

Figure 25 :
Figure 25: Radial mass fraction of OH at x/d = 20 predicted by different table sizes.
[26] paper, we have proposed a new simulation technique based on the lookup table, where the lookup table building and combustion simulation were performed simultaneously.The reaction rates were estimated by direct integration by solving ODEs and storing them in the lookup table.The proposed technique was applied to the EDC model.Furthermore, a CO-H 2 -air turbulent nonpremixed flame was computed.The results obtained using our model were compared with experimental data reported by Correa and Gulati[26]as well as with those obtained by using the EDC model with direct integration.It was found that the prediction accuracy of our technique is very close to that of the EDC model.Thus, the prediction accuracy depends on the combustion model that is accelerated by our technique.Further, the computation time for the technique using the lookup table was approximately 99.3% lower than that using the EDC model.The memory sizes of the lookup table with 50 × 50 × 50 (15 Mbyte), 100 × 100 × 1 (1.2 MByte),