Energy & Fuels 2005, 19, 1873-1878
1873
Space Discretization Error of Methane Combustion Simulations in Turbulent Flow J. Lindberg* Division of Energy Engineering, Lulea˚ , University of Technology, SE-97187 Lulea˚ , Sweden
M. J. Cervantes Division of Fluid Mechanics, Lulea˚ , University of Technology, SE-97187 Lulea˚ , Sweden Received March 16, 2005. Revised Manuscript Received June 15, 2005
Numerical investigation of methane combustion in a pipe with turbulent flow is studied. The space discretization error is investigated quantitatively and qualitatively, using the Richardson extrapolation and profiles comparisons. Comparison of the profiles indicates that the solution converges to a grid-independent solution. The Richardson method gives unsatisfactory results to determine the grid error, because of the rigidity of the method. A second-order polynomial is used as an alternative to the Richardson method. The results are more stable and have a better goodness of fit. The results of the simulations are compared with those of a similar experiment and the corresponding analytical solution.
Introduction The facility to obtain results in computational fluid dynamics (CFD) hides a complex process, which has many different sources of error. A careful analysis and quantification of the different sources of error is a necessity to assess the predictability of a model. Three different sources of errors are generally present in a simulation: modeling errors, iterative errors, and discretization errors. CFD software is used to model combustion processes in a variety of applications within the field of energy engineering. The large number of equations that need to be solved makes such simulations time-consuming and costly, in regard to computational resources. To keep the computational effort at a reasonable level, simple two equations (turbulence models and combustion models) with only 1-5 reaction steps are generally used with a relatively coarse grid. Furthermore, simulation of gas combustion involving chemical kinetics is often phased with numerical convergence problems, which are essentially due to the rigidity in the formulation of the chemistry.1 Therefore, today, the tolerances of a combustion simulation are relatively large and generally a qualitative assessment of the simulations is often sufficient. A systematic method, which quantifies the error, is necessary to increase trust and quality, especially for the grid dependency of the solution (i.e., the space discretization error). Today, the Richardson extrapolation method is widely used for error analysis within CFD (see, e.g., Bergstrom * Author to whom correspondence should be addressed. Phone: +46 70 677 62 73. Fax: +46 920 491074. E-mail:
[email protected]. (1) Norton, D. G.; Vlachos, D. G. Chem. Eng. Sci. 2003, 58, 48714882.
and Gebart2 for a detailed description of the method). In the method, simulations with at least three different grids but similar topology are performed. A suitable variable is chosen, and the results of which are fitted to the equation
y ) A + BhR
(1)
where A is the extrapolated solution for an infinite number of cells from which a grid error can be calculated, h the inverse averaged grid cell edge size, and R the average order of the discretization scheme that is used. Subtraction of the result obtained with the finest grid to A gives the grid error. The result of the simulations must be in the asymptotic region to obtain a reliable result. Unfortunately, the asymptotic region can only be guessed, which is the main drawback of the method. Such estimation may also be obtained analytically if the number of cells is doubled between each grid. The Richardson extrapolation method is based on the fact that the order of the discretization scheme is between 0 and 2. The results present often some noise, mostly due to the difficulty in obtaining topologically similar grids. The noise and the stiffness involved in the determination of the equation parameters in the Richardson method (three equations and three unknowns) often gives unsatisfying results. This is especially true with tetra mesh, where similar grids are difficult to obtain. Therefore, an extension of the Richardson’s method to an ordinary polynomial will allow some variation in the results but still allow an estimation of the error. (2) Bergstro¨m, J.; Gebart, R. Int. J. Numer. Methods Heat Fluid Flow 1999, 9 (4), 472-486.
10.1021/ef050066l CCC: $30.25 © 2005 American Chemical Society Published on Web 08/10/2005
1874
Energy & Fuels, Vol. 19, No. 5, 2005
Lindberg and Cervantes
The objective of the paper is to investigate the solution of a methane combustion simulation using Nicol’s3 reaction model with the standard energy-dissipation rate (k-) turbulence model in a pipe. Such a case is attractive for the simplicity of the geometry and the fact that the standard k- turbulence model is calibrated from pipe flows. Furthermore, experiments on a similar test case have been performed. Special attention is given to the space discretization error through qualitative and quantitative methods such as profile comparison, the Richardson extrapolation, and polynomial extrapolation. The results then are discussed in detail and compared to the experiments and the analytical solution. Simulation of Methane Combustion in a Turbulent Flow The simulated flow consists of a mixture of methane, air, carbon dioxide (CO2), and water vapor (H2O), which is supplied through the inlet of a pipe at a temperature of 1073 K. The pipe has a diameter of 0.035 m and a length of 1 m. The calculations are initiated with a domain temperature at 1500 K to ensure combustion. The temperature of the pipe wall is kept constant at 323 K, to simulate a cooled combustion chamber. Simulations were performed using the segregated solver of Fluent 6 and double precision. Models. Obtaining numerical simulations on combustion is demanding, in regard to the computational resources, because of the large number of differential equations to be solved. The axis symmetry of the geometry allows two-dimensional simulations, which substantially reduces the number of differential equations to be solved and the number of cells. The flow is modeled using the standard turbulence k- model, the radiation heat transfer is modeled using the discrete ordinates radiation model, and the chemical process of combustion is modeled with the Arrhenius4 theory. The reaction model described by Nicol et al.,3 consisting of the three reaction equations
3 CH4 + O2 f CO + 2H2O 2
(I)
1 CO + O2 f CO2 2
(II)
1 CO2 f CO + O2 2
(III)
is used. Because the chemical reaction rate is closely coupled to turbulence, the turbulence chemistry coupling is of great importance. For turbulent flames, the Eddy Dissipation Concept model (EDC),5 which is an expansion of the Eddy Dissipation models,6 is used in this study. The EDC model presumes that chemical reaction occurs in small turbulent structuressso-called “fine-scales”, where the mixing is rapid. These fine(3) Nicol, D. G.; Malte, P. C.; Hamer, A. J.; Roby, R. J.; Steele, R. C. J. Eng. Gas Turbines Power 1999, 121, 272-280. (4) Turns, S. R. An Introduction to CombustionsConcepts and Applications, 2nd Edition; McGraw-Hill: Boston, 2000; Chapter 4, pp 111-147. (5) Byggstøl, S.; Magnussen, B. F. Turbulent Shear Flows 4; Springer: Berlin, 1985; pp 381-395. (6) Magnussen, B. F.; Hjertager, B. H. In Proceedings of the 16th Symposium on Combustion; The Combustion Institute: Pittsburgh, PA, 1976; pp 719-729.
scales can be viewed as plug-flow reactors (PFRs), where reactants are mixed on a molecular level. The residence time of the components in the fine-scales, together with the time of the chemical reactions, describe the reaction rate. Boundary Conditions. At the inlet boundary, a gas mixture with a mass flow rate of 0.00685 kg/s is applied. The assortment contains the following gases, with their respective mass fraction: methane (CH4 ) 0.023), carbon dioxide (CO2 ) 0.047), water vapor (H2O ) 0.0386), oxygen (O2 ) 0.15), and nitrogen (N2 ) 0.7414). Turbulence is specified using the hydraulic diameter and turbulence intensity. Outlet boundary is defined as the pressure outlet; the wall temperature is fixed to 323 K, and the internal emissivity is set to 0.7. Grid and Schemes. The domain, a two-dimensional plane with a symmetry axis along the center line of the cylinder, was discretized using a two-dimensional structured mesh. Four grids have been investigated: 156 050, 247 455, 397 800, and 643 134 quadrilateral cells. Each grid is divided in two regions of different grid density. The finest region of the grid is found between the inlet and a distance x ) 0.4 m, with a density of 2.7, compared to the coarser part. The finer region is homogeneous and its purpose is to capture the large gradient generated near the inlet due to the rapid changes in properties of the different variables describing the process. The coarse region of the grid is located from a distance x ) 0.4 m to the end of the pipe, where the combustion process is supposed to be completed and the gradients are, therefore, negligible. The cells of the coarser grid are elongated in the x-direction with a factor of 1.002, to minimize the number of cells. Pressure is interpolated using the second-order interpolation scheme, and the SIMPLE algorithm is used for the pressure velocity coupling. The second-order scheme reconstructs the face pressure in a manner that gives accurate second-order convection terms. A secondorder upwind discretization scheme is applied to the transport and conservation equations. The underrelaxation factors were equal to 0.5 for the pressure, the velocity, and the concentration of the different species. For the other variable, the under-relaxation factor was equal to 0.8. Assessment of the Space Discretization Error. The grid has a fundamental importance in the result of a simulation. Therefore, careful grid error analysis is necessary to assess the quality of a simulation. Several methods exist, both qualitative and quantitative. Comparison of a profile obtained with different grids is one method. If the profiles are similar in shape, a gridindependent solution is assumed. This method is qualitative and does not allow a quantification of the grid error of special interest to assess the quality of a combustion model. To quantify the grid error, the Richardson method is an alternative. In this method, at least three simulations with different space discretizations are necessary. The result of the simulations must be in the asymptotic region to obtain a reliable result. Unfortunately, the asymptotic region can only be guessed, which is the main drawback of the method. The results of the simulations are fitted to eq 1 (presented earlier in this paper). The stiffness involved in the determination of
Methane Combustion Simulations in Turbulent Flow
Energy & Fuels, Vol. 19, No. 5, 2005 1875
Table 1. Scaled Residuala grid 1 continuity x-velocity y-velocity energy turbulent kinetic energy, k turbulent dissipation rate, D-O intensity CH4 O2 CO2 CO H2O a
grid 2
10-6
10-6
7.33 × 3.06 × 10-6 4.98 × 10-7 9.99 × 10-7 6.94 × 10-6 1.61 × 10-5 7.59 × 10-7 2.36 × 10-5 4.24 × 10-6 5.99 × 10-6 1.91 × 10-4 2.69 × 10-7
7.63 × 2.07 × 10-6 2.00 × 10-7 7.62 × 10-7 5.69 × 10-6 9.97 × 10-6 9.96 × 10-7 1.13 × 10-5 3.10 × 10-6 4.58 × 10-6 1.37 × 10-4 1.36 × 10-7
grid 3 10-5
2.39 × 2.64 × 10-6 1.69 × 10-7 9.99 × 10-7 7.35 × 10-6 1.43 × 10-5 9.10 × 10-7 1.04 × 10-5 2.05 × 10-6 3.04 × 10-6 8.94 × 10-5 5.58 × 10-8
grid 4 2.51 × 10-6 4.51 × 10-7 6.05 × 10-8 2.01 × 10-7 1.47 × 10-6 4.90 × 10-6 4.24 × 10-7 8.86 × 10-7 5.66 × 10-6 8.05 × 10-6 2.35 × 10-5 3.29 × 10-6
Grid 1 is the coarsest and grid 4 is the finest.
the equation parameters often gives unsatisfying results, such as a scheme order of 10 and unrealistic values of the extrapolated solution. However, even in the asymptotic region, the Richardson method has the tendency to give large errors in the extrapolated value when small variations of the parameter under study appear. This variation appears often when it is difficult to have a global refinement of the mesh for untrivial geometries. The order of the schemes involved in CFD simulations (equal to 2 or below) and the use of nonuniform grids gives an average scheme order of less than 2. Therefore, a curve fit to a polynomial of second order should be sufficient to obtain a value of the extrapolated value in cases where the Richardson method fails. Results The convergence criteria are first presented. An assessment of numerical error due to space discretization follows, with discussion of the simulations results, which are compared to the experiments and analytical solutions. Convergence Criteria. The scaled residual was used as a convergence criteria for the simulations. The scaled residual of a general variable φ is defined by
Rφ )
∑cells P|∑nb anbφnb + b - aPφP.| ∑cells P |aPφP|
(2)
where aP is the center coefficient, anb the influence coefficients for neighboring cells, and b the contribution of the constant portion of the source term. The convergence criteria for the scaled residual was set to 10-6 for the temperature and the D-O intensity. The other parameters (i.e., mass, velocity, turbulent kinetic energy, turbulent dissipation, and species concentration) had a limit of 10-3. The reached scaled residuals are presented in Table 1, where the highest value is observed for the carbon monoxide: R[CO] ) 1.37 × 10-4. Space Discretization Error. Carbon dioxide (CO2) is a product of the complete combustion of hydrocarbons and is present in a large concentration in most parts of the combustion chamber. The temperature is strongly associated with the combustion process and heat transfer to the wall. This makes these two parameters suited to study the simulation and, more particularly, the space discretization error. The temperature and CO2 concentration along the centerline of the pipe and at x
) 0.5 m are presented in Figure 1a-d. The scale in Figure 1d is enhanced for better resolution. The profiles for the different grids are similar, which indicates a grid-independent solution. A regression analysis indicates a relatively good fit using the Richardson extrapolation. The result of the extrapolations gives poor prediction of the theoretical discretization order. The values of R in Table 2 have a relevance for 3 of the extrapolations (Tw,0.5, T1,out, T2,out). The value of A represents the solution for an edge length of h ) 0 (i.e., the exact solution). A second-order polynomial has also been fitted to the data (see Table 3). The cross correlation r2 indicates a better fit for the second-order polynomial than with the Richardson extrapolation. This indicates a better local description of the data with a second-order polynomial. Figure 2 illustrates the Richardson extrapolation and the second-order polynomial fit of the area-weighted average of the temperature at the outlet and the mass fraction of carbon monoxide (CO) at x ) 0.5 m. If some variations in the data seem to be fitted (r2 * 1), the Richardson extrapolation can become unstable as h approaches zero, thereby explaining the sudden drop in Figure 2. Results from Simulations. Figure 3 presents the variation of the different parameters describing the physics and chemistry of the combustion process. The gas mixture is supplied through the pipe inlet at x ) 0 m. Combustion starts and the concentration of methane decreases as the concentration of CO increases through the process of combustion, according to reaction I. The concentration of CO2 is stable until x ) 0.01 m, where the increase in CO concentration initiates reaction II and the production of CO2 (see Figure 3). The majority of the energy is released during the methane combustion, which results in a large temperature increase up to x ) 0.1 m, where all methane is consumed. The rest of the energy released is generated through the second reaction, where CO is oxidized to CO2 (reaction II), resulting in a temperature maximum at x ) 0.25 m that is due to the completion of these two reactions. The maximum temperature of the simulation is 1895 K, compared to an estimation of the adiabatic flame temperature at 1996 K; thus, the result of the simulation is reasonable. With the cold wall of the pipe, a deviation from the adiabatic flame temperature is expected. The pressure decreases from the inlet of the pipe up to x ) 0.1 m and then increases again toward the outlet.
1876
Energy & Fuels, Vol. 19, No. 5, 2005
Lindberg and Cervantes
Figure 1. Profiles of the temperature and mass fraction of carbon dioxide (CO2): (a) temperature profiles along the center line, (b) temperature profiles of the cross section at x ) 0.5 m, (c) [CO2] profiles along the center line, and (d) [CO2] profiles of the cross section at x ) 0.5 m. Note the change in scale. Table 2. Results from the Richardson Extrapolationa A B r2 R
A B r2 R
T1,0.5
T2,0.5
Tw,0.5
T1,out
T2,out
Tw,out
1251 375 0.966 -0.13
2899 -1384 0.975 0.04
1407 -73.3 0.9712 0.625
1243 -35.3 0.9421 1.315
1154 -16.1 0.9237 1.965
1024 -10 0.9038 3.355
CO21,0.5
CO22,0.5
CO2w,0.5
CO21,out
CO22,out
CO2w,out
0.1061 0.00037 0.002 -0.105
0.1123 -0.0056 0.1705 0.1
0.1201 -0.0131 0.4346 0.1
0.0939 0.0143 0.963 0.1
0.08829 0.01975 0.9533 0.1
0.08800 0.01991 0.9533 0.1
a Data for various temperatures T and CO concentrations are 2 given. A, B, and R are values from the Richardson function (eq 1); r2 is the goodness of fit. Subscript “1” represents y ) 0.001 m, subscript “2” represents y ) 0.00875 m, and subscript “w” represents the area-weighted average. Subscript “0.5” represents x ) 0.5 m, and subscript “out” represents x ) 1 m (the outlet).
The pressure minimum coincides with the temperature and velocity maximum, which gives the largest losses. The velocity increases as the temperature increases, in accordance with the ideal gas law until the temperature is at its maximum at x ) 0.25 m. Thereafter, the temperature decreases, as a result of heat loss to the cooled wall, and the velocity decreases accordingly. CO2 is a product of complete combustion and its concentration increases as the methane is oxidized to CO, which, in turn, is oxidized to CO2, through reaction II. Assuming complete combustion, an analytical solution of the theoretical concentration of CO2 is calculated to a mass fraction of 0.110. The maximum CO concentration observed in the simulations is a mass fraction of 0.109. Taking into account the presence of a small concentration of CO of a mass fraction of ∼0.001 in the
Table 3. Results from Curve Fitting of a Second-Order Polynomiala A B C r2
A B C r2
T1,0.5
T2,0.5
Tw,0.5
T1,out
T2,out
Tw,out
1691 -80.2 15.9 0.9719
1586 -86 15.7 0.9788
1388 -60.9 7 0.9723
1251 -38.1 -4.8 0.9415
1150 4.9 -17.1 0.9238
1010 26.3 -22.4 0.9059
CO21,0.5
CO22,0.5
CO2w,0.5
CO21,out
CO22,out
CO2w,out
0.1107 -0.006 0.002 0.9881
0.1126 -0.0081 0.0026 0.777
0.1035 -0.1 0.003 0.7826
0.1059 0.0029 -0.0006 0.9736
0.1042 0.0049 -0.0011 0.9791
0.1035 0.0057 -0.0014 0.9991
a Data for various temperatures T and CO concentrations are 2 given. A, B, and C are coefficients of the second-order polynomial; 2 r is the goodness of fit. Subscript “1” represents y ) 0.001 m, subscript “2” represents y ) 0.00875 m, and subscript “w” represents the area-weighted average. Subscript “0.5” represents x ) 0.5 m, and subscript “out” represents x ) 1 m (the outlet).
simulations not included in the theoretic calculation, the results from the simulations agree well with theory. Measurements from a similar experiment have been performed. The only difference with the simulation was the wall boundary conditions, which were cooled with the surrounding air through natural convection. The experimental setup consisted of a vertical circular duct with a diameter of 0.035 m and a length of 3 m. The fuel and oxidizer were premixed and supplied through the inlet at the bottom at a velocity of ∼6 m/s at a temperature of 293 K. The combustion gases exited through the top, where an open flame incinerated any products of incomplete combustion. The gas was continuously ignited through spark ignition before it passed through a stainless-steel net that acted as a flame
Methane Combustion Simulations in Turbulent Flow
Figure 2. Results of the simulation, Richardson extrapolation, and second-order polynomial fit: (a) fit Tw,out and (b) fit to CO2w,0.5.
holder. The concentrations of the combustion gases were measured 0.35 m above the flame holder, where a watercooled probe was inserted horizontally into the combustion chamber. The probe was fitted with a ceramic tip that was 0.003 m in diameter and closed at the end, and it had a hole for sampling (0.0015 m in diameter) that was facing the flow. The temperature was measured by three thermocouples centered in the cylindrical duct and placed 0.9 m apart. The temperature of the wall was measured by three thermocouples that were placed in contact with the outside surface of the wall of the combustion chamber (Figure 4). The measurements indicate a mass fraction of CO2 that is equal to 0.122, which corresponds to the extrapolated value using the Richardson or second-order polynomial method. At high temperatures and a high concentration of CO2, the oxidized CO can be reformed through dissociation (reaction III). When the CO2 concentration is at its maximum at x ) 0.3 m and the temperature is equal to 1700 K, an increase in CO can be observed due to reaction III. As the number of molecules in the system increases, accompanied by CO2 dissociation, a slight increase in velocity, corresponding to the change in volume, can be observed. The change in CO concentration in Figure 3 should not be overestimated. Because values are normalized with their maximum value, a change in the CO concentration will seem larger than the same change in the CO2 concentration. Discussion A study of the temperature and CO2 profiles indicates a result converging to a grid-independent solution.
Energy & Fuels, Vol. 19, No. 5, 2005 1877
Figure 3. Variation of the pressure, velocity, temperature, [CH4], [CO2], and [CO] along the centerline. The parameters are normalized to their maximum values; maximum pressure ) 115.9 Pa, maximum velocity ) 38.9 m/s, and maximum temperature ) 1895.3 K, and the species mass fractions are CH4,max ) 0.023, CO2,max ) 0.109, and COmax ) 0.0106.
Figure 4. Experimental setup in Paper II.
The differences observed between the different grids are only a few percent. However, using a quantitative method such as the Richardson extrapolation to evaluate the grid error, the values obtained for the discretization order R indicate a solution far from the asymptotic range. This is certainly not the fact, and the poor results from the method are attributed to its rigidity. Results from extrapolation using the Richardson theory also vary, depending on the parameter studied, indicating that several parameters must be studied when performing a Richardson extrapolation. The engineers evaluating combustion simulations are often limited to a qualitative assessment of the results.
1878
Energy & Fuels, Vol. 19, No. 5, 2005
Here, a second-order polynomial fit can become a more useful tool in giving a qualitative assessment of the space discretization error. The fit of a second-order polynomial always resulted in a better goodness of fit than that with the Richardson extrapolation. The extrapolated solution for h ) 0 is more consequent than with the Richardson method. Further investigation of the method is necessary. The extrapolated values for the mass fraction of CO2 at x ) 0.5 m are comparable with the analytical solution and the experimental measurements. Comparison between the maximum temperature and the estimated adiabatic flame temperature also indicates that the results from the simulations are trustworthy.
Lindberg and Cervantes
Conclusion Although the profiles of temperature and carbon dioxide (CO2) agree well, indicating a result that converges to a grid-independent solution, the Richardson extrapolation method does not produce results in the asymptotic region. This is a result of the stiffness of the method, which does not allow variations in the results. A second-order polynomial fit presented a better goodness of fit to the data set than that obtained with the Richardson method. Further investigation of the second-order polynomial method is necessary. EF050066L
