Abstract
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.
1. 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β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β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β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 H2 combustion [20], CH4 combustion [21], and 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 H2-air turbulent nonpremixed flame [23], a CO-H2-air turbulent nonpremixed flame [24], and a CH4-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 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-H2-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.
2. 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 user-defined 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; , the velocity in the direction; , the enthalpy; , the coordinate in the direction; , the turbulent viscosity; , the effective Prandtl number; , the thermal conductivity. The symbols β and ~ on the variables represent the time and density average, respectively. is the temperature computed as where . denotes the standard enthalpy of formation of chemical species ; , the mass fraction of chemical species ; , the mean specific heat at constant pressure of chemical species .
was computed as where denotes the specific heat at constant pressure of chemical species .
2.2. 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 is the mass fraction and and are the mass fractions in the fine structure and surrounding gas, respectively. is the time, and is the reaction rates in the fine structure. denotes the length of the fine structures and is expressed as where , denotes the viscosity; is the turbulent dissipation rate, and denotes the turbulent kinetic energy. denotes the fine structure region residence time, and is given as where .
The mean reaction rate of chemical species is computed as
2.3. 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 CO2 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 nonuniformly divided since the time scale ranges from 1.0β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 two-dimensional lookup table is illustrated in Figure 2. If the variables, that is, the mixture fraction and progress variable, are at the sampled point , then the reaction rates of the chemical species are stored at the points , , , and 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 and denote the new and old data of the reaction rates of the chemical species at the stored points , , and , respectively. denotes the reaction rates at the sampled point . ββis calculated as where is the mass fractions of chemical species, and denotes the mass fractions of chemical species in the fine structure. Thus, the mean reaction rates are computed as where denotes the reaction rates of chemical species obtained by the interpolation.
and ββare the th weight function and sum of the weight function, which are computed as where ββis converted to an integer and varies from 1 to 8. In the above equations, , is the distance between the sampled point and the stored points , , , and , is the number of the writing reaction rates.
The weight function is shown in Figure 3. indicates that coordinates of the stored and sampled points are equal, which means the data obtained at is favorable for the lookup table. expresses the coefficient of the old data of the reaction rates, while implies the coefficient of the new data of the reaction rates. Thus, the data is well updated when . Once the reaction rates are written in the stored points , , , and , 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.0β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 denotes the norm; is the error value; is the reaction rates of chemical species at the sampled point ; is the reaction rates obtained by the interpolation; and ββis the scale factor chosen from the maximum .
The mixture fraction is computed as where denotes the mixture fraction, and is the effective Schmidt number.
In this study, the progress variable is expressed as where denotes the mass fraction of CO2.
3. Results
3.1. Numerical Conditions
A CO-H2-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, H2, and N2 in the fuel were 0.393, 0.0332, and 0.574, respectively. In the ambient air, the mass fractions of O2 and N2 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.0βK. A second-order upwind difference was applied to the governing equations of the momentum, mass fraction of chemical species , 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 and the Schmidt number was . 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] was used to compute the reaction rates. It involved 15 species, including H, O2, O, OH, H2, H2O, HO2, H2O2, CH, CH2O, CO, CO2, HCO, N, and N2.
3.2. 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.5β5βm and 1β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.6β5βm and 6.7β6βm, respectively.
Figures 5 and 6 show plots of the radial mixture fractions predicted with the two different grids, where denotes the radial coordinate and 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.
3.3. 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 . In all the figures, 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.
3.4. 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 . It can be noted from Figures 14, 15, and 16 that the mixture fraction and mass fractions of CO and H2 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 H2O is predicted well in spite of the fact that the radial profile of the mass fraction of H2O is slightly different when and when . The peak of the mass fraction of H2O 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 H2 obtained by using our technique are slightly underestimated relative to the experimental data when but are close when . However, the mass fraction of H2 obtained by the laminar flamelet model is better than that by our technique. The radial mass fraction of H2O computed by our technique is similar to the experimental data, while the laminar flamelet model underestimates the peak of the mass fraction of H2O.
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.
3.5. 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 our technique because the mass fractions equations of the chemical species are not required in the laminar flamelet model. The proposed technique is faster than the EDC model with the ISAT technique in FLUENT.
3.6. 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 (15βMbyte), (1.2βMbyte), and ββ(960βMbyte), respectively. The table with 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 and that with are in good agreement, while that obtained using the lookup table with is different. Note from Figure 25 that the discrepancy is minimal, although both peak values of OH exhibit differences at . However, the result obtained using the lookup table with is different from the mass fraction of OH obtained using the look-up table with and . Thus, if there is insufficient memory in the computer, the lookup table with can be used, but the reaction time scale cannot be eliminated.
4. Conclusions
In this 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-H2-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 (15βMbyte), (1.2βMByte), and (960βMbyte) were tested; it was found that if there was insufficient memory in the computer, the lookup table with could be used, but the reaction time scale cannot be eliminated.