- Email: [email protected]

Evaluation of the laminar diffusion flamelet model in the calculation of an axisymmetric coflow laminar ethylene–air diffusion flame Fengshan Liu ∗ , Hongsheng Guo, Gregory J. Smallwood Institute for Chemical Process and Environmental Technology, National Research Council, Building M-9, 1200 Montreal Road, Ottawa, Ontario, Canada K1A 0R6 Received 22 September 2004; received in revised form 1 September 2005; accepted 10 September 2005 Available online 27 October 2005

Abstract A numerical study of an axisymmetric coflow laminar ethylene–air diffusion flame at atmospheric pressure was conducted using detailed chemistry and complex thermal and transport properties and two different methodologies: (1) the direct simulation method of solving the two-dimensional axisymmetric elliptic governing equations, and (2) the steady-state stretched diffusion flamelet model. Soot formation and radiative heat transfer were not taken into account in these calculations, both for simplicity and to avoid the complications associated with the issues of how to incorporate these chemical and physical processes into the flamelet model. The same reaction mechanism and thermal and transport properties were used in the 2D direct simulation and the generation of the flamelet library. The flamelet library was generated from the solutions of counterflow ethylene–air diffusion flames at a series of stretch rates. Results of the 2D direct simulation and the flamelet model are compared in physical space. Although the overall results of the flamelet model are qualitatively similar to those of the direct simulation, significant differences exist between the results of the two methods even for temperature and major species. The direct simulation method predicts that the peak concentrations of CO2 and H2 O occur in different regions in the flame, while the flamelet model results show that the peak concentrations of CO2 and H2 O occur in the same region. The flamelet model predicts an overly rapid approach to the equilibrium structure in the downstream region, leading to significantly higher flame temperatures. The main reason for the failure of the flamelet model in the downstream region is due to the neglect of the effects of multidimensional convection and diffusion and the fundamental difference in the chemical structure between a coflow diffusion flame and a counterflow diffusion flame. The findings of this paper are highly relevant to understanding the flamelet model results in the calculations of multidimensional turbulent diffusion flames. Crown Copyright 2005 Published by Elsevier Inc. on behalf of The Combustion Institute. All rights reserved. Keywords: Flamelet model; Coflow diffusion flame; Flow field effect

1. Introduction * Corresponding author. Fax: +1 613 957 7869.

E-mail address: [email protected] (F. Liu).

Since the comprehensive review of the steadystate stretched laminar diffusion flamelet model by

0010-2180/$ – see front matter Crown Copyright 2005 Published by Elsevier Inc. on behalf of The Combustion Institute. All rights reserved. doi:10.1016/j.combustflame.2005.09.005

606

F. Liu et al. / Combustion and Flame 144 (2006) 605–618

Peters [1], it has become a popular alternative modeling approach for laminar flames [2–4], for turbulent diffusion flames [5–10], and in many other applications, e.g., fires [11] and lifted turbulent flames [12]. The distinct advantage offered by this model, compared to conventional turbulent combustion models, is that flow field and chemical kinetics calculations are essentially decoupled. The influence of the flow field on the flamelet structure is assumed to be solely represented by the scalar dissipation rate, which can be interpreted as the inverse of a characteristic diffusion time [1]. As such, detailed chemistry and complex thermal and transport properties can easily be incorporated into multidimensional calculations of laminar and turbulent nonpremixed combustion through the concept of the stretched flamelet model. An underlying fundamental assumption made in the flamelet theory is that multidimensional effects (convection and diffusion in directions other than that normal to the local stoichiometric surface) are negligibly small. The consequence of this assumption is that there exists a universal state relationship of thermochemical scalars (temperature, density, and species) as functions of the mixture fraction parameterized at a characteristic scalar dissipation rate in diffusion flames of different configurations. Based on this assumption, computationally expensive calculations of multidimensional laminar diffusion flames can be conducted very efficiently through the flamelet model, as demonstrated in several studies [2–4,13]. Motivated by the promising approach offered by the flamelet model to simulating multidimensional laminar diffusion flames, several studies have made attempts to evaluate the diffusion flamelet model in the calculation of axisymmetric coflow laminar methane–air diffusion flames over the last decade or so [13–15]. The findings of these studies are inconsistent and somewhat controversial. Nishioka et al. [15] found that there is a remarkable similarity in the structures of the coflow and counterflow diffusion flames, both in the physical plane and in the mixture fraction plane. The agreement is even better for molar production rates than scalars themselves. However, it is worth noting that the computational domain of Nishioka et al. [15] covered less than half of the flame and the comparison was limited to the near burner region. On the other hand, Smooke et al. [14] found that there are substantial differences in the temperature and major species profiles between the coflow and the counterflow flames even within the first half of the coflow flame. The study of Heyl and Bockhorn [13] showed that there is a rather significant difference between the 2D calculation and the flamelet model in the temperature distribution along the flame centerline. Clearly, further study is required to evaluate the

performance of the flamelet model in the calculation of multidimensional laminar diffusion flames. The laminar flamelet model has been popularly used in the literature to conduct numerical simulation of turbulent nonpremixed combustion in various applications, from turbulent jet flames to fire spread in compartments. However, relatively little attention has been paid to the fundamental aspects of the flamelet model in applying this model to multidimensional flame modeling. A careful and direct evaluation of the laminar flamelet model, such as the one presented in this paper, has not been published in the literature. The two previous studies conducted by Smooke et al. [14] and Nishioka et al. [15] constructed the relationships of scalar versus the mixture fraction based on 2D numerical results and compared them with the flamelet model results. In addition, these studies were limited to the near-burner field. Although the findings of Smooke et al. and Nishioka et al. are relevant to the validity of the flamelet model in 2D flame calculations, they do not reveal the direct consequences of applying the laminar flamelet model to calculate a 2D laminar coflow diffusion flame. This study reveals for the first time the quantitative differences between the results of the laminar flamelet model and the 2D direct numerical simulation in the physical space in an atmospheric coflow laminar diffusion flame. While it is possible to evaluate the flamelet model in the calculation of a turbulent jet flame and compare the results with experimental data, such an evaluation is subject to uncertainties from both numerical (namely turbulence model) and experimental (measurement errors) sides. The present study does not suffer such uncertainties and can best reveal the fundamental drawback of the flamelet model. In the present study, the steady-state laminar diffusion flamelet model was evaluated in the calculation of an axisymmetric laminar coflow ethylene diffusion flame by comparing its results against those from direct numerical simulation in the physical space. Soot formation and radiation heat transfer were not taken into account in these calculations, both for simplicity and to avoid the complications associated with the issues of how to incorporate these chemical and physical processes into the flamelet model. The same reaction mechanism and thermal and transport properties were used in both calculations.

2. Theoretical and numerical models 2.1. Laminar flamelet model The laminar flamelet model has been extensively reviewed by Peters [1] and Bray and Peters [16]. The development of the flamelet model was based on

F. Liu et al. / Combustion and Flame 144 (2006) 605–618

rewriting the energy and species equations in a local coordinate system and subsequent asymptotic analysis of the order of magnitude of individual terms [1]. In diffusion flames, combustion takes place in a thin layer in the vicinity of the surface of the stoichiometric mixture if the local mixture fraction gradient is assumed to be sufficiently high. This thin layer and the surrounding nonreacting mixing region are defined as a laminar diffusion flamelet by Peters [1]. Based on the conservation equations written in Eulerian form, assuming unity Lewis number, and using the Croccotype coordinate transformation, Peters [1] has derived the unsteady governing equations of species and temperature in the mixture fraction space, ρ

1 ∂ 2 Yi ∂Yi = χ + ω˙ − R(Yi ), ∂t 2 ∂Z 2

(1)

ρ

1 ∂ 2 T hi qr ∂T = χ − ω˙ i + − R(Yi ), 2 ∂t 2 ∂Z cp cp

(2)

n i

where χ is the scalar dissipation rate, defined as χ = 2D(∂Z/∂xk )2 , and D is a common diffusion coefficient. The new coordinates Z, Z2 , and Z3 are respectively the mixture fraction, which is locally normal to the surface of stoichiometric mixture, and the other two independent variables (Z2 = x2 , Z3 = x3 ). The operator R is defined as ∂ ∂ + v3 R = ρ v2 ∂Z2 ∂Z3 ∂(ρD) ∂ ∂(ρD) ∂ − − ∂x2 ∂Z2 ∂x3 ∂Z3 3 ∂Z ∂ 2 ∂2 2 − ρD + (3) ∂xk ∂Z∂Zk ∂Zk2 k=2 and it physically represents effects of multidimensional convection, diffusion, and curvature in the two dimensions other than that normal to the surface of the stoichiometric mixture. In the limit of the thin reaction zone, where the mixture fraction gradient is high and the scalar dissipation rate is large, all terms involving derivatives with respect to Z2 and Z3 are of lower order than other terms in Eqs. (1) and (2). Under these conditions, the flamelet formulation reduces to a one-dimensional structure normal to the surface of the stoichiometric mixture. The unity Lewis number assumption in the derivation of the species equation, Eq. (1), can be removed, as has done by Mauß et al. [17]. The steady-state form of Eqs. (1) and (2) (with the neglect of the operator R) can be solved for the flamelet structure in the mixture fraction coordinate with properly defined boundary conditions at Z = 0 (oxidant side) and Z = 1 (fuel side) and a specified

607

characteristic scalar dissipation rate. The flamelet solutions for a series of scalar dissipation rates ranging from a very small value (the equilibrium limit) to a large value (the quenching limit) form what is widely known as the flamelet library. An alternative way to establish the flamelet library is to solve the governing equations of counterflow diffusion flames along the stagnation streamline (physical space) at a series of specified stretch rates. The link between the physical space and the mixture fraction coordinate is a properly defined mixture fraction. Once the solutions of these counterflow diffusion flames in the physical space are obtained, the mixture fraction distribution can be obtained in postprocessing, either by solving a transport equation for the mixture fraction or by using the mixture fraction definition of Bilger et al. [18]. It has been shown by Chou et al. [3] that the solutions (the structures of the flamelet in the mixture fraction space) obtained by these two different methods are not identical, but the differences remain relatively small. A detailed comparison of these two different methods for generating the flamelet library is beyond the scope of the present study. In this study, the latter approach was used, since it allows identical reaction mechanisms and thermal and transport properties to be consistently used in both the counterflow flame and the 2D coflow flame calculations. It is worth pointing out that preferential diffusion (nonunity Lewis number effect) was retained in the calculations of both the counterflow diffusion flame (for the generation of the flamelet library) and the 2D coflow diffusion flame. The governing equations and numerical methods for the counterflow diffusion flame have been described previously in [19]. The numerical results of thermochemical properties in the physical domain were then converted to the mixture fraction space using two different approaches. In the first approach, the following transport equation for the mixture fraction Z, suggested by Pitsch and Peters [20], was solved, d dZ dZ = ρDZ , ρu (4) dx dx dx where DZ = λ/ρcp with λ, ρ, and cp being the mixture thermal conductivity, density, and specific heat at constant pressure, respectively. In the second approach, the mixture fraction was defined based on the elemental mass fractions of carbon, hydrogen, and oxygen [18], Z=

2ξC /WC + 0.5ξH /WH + (ξO,2 − ξO )/WO , (5) 2ξC,1 /WC + 0.5ξH,1 /WH + ξO,2 /WO

where ξj and Wj are the elemental mass fractions and atomic masses, respectively. Subscripts 1 and 2 represent fuel and oxidizer stream, respectively. The scalar dissipation rate χ as a function of the mixture fraction is then calculated according to its definition (using

608

F. Liu et al. / Combustion and Flame 144 (2006) 605–618

DZ ). The objective of solving the diffusion flame in the counterflow configuration at a series of specified stretch rates is to establish the flamelet library, which consists of not only the thermochemical scalars such as temperature, density, and mass fractions of species, but also the thermal and transport properties of the mixture, i.e., the viscosity µ and the diffusion coefficient DZ , which are required in the calculation of the axisymmetric coflow diffusion flame using information in the flamelet library. The correspondence between the solution in the flamelet library to the coflow flame is established by a characteristic scalar dissipation rate χc and the mixture fraction. There exists a certain ambiguity when the flamelet library is compiled in this way, since the resultant scalar dissipation rate at a given stretch rate varies across the flamelet, i.e., a function of the mixture fraction. The characteristic scalar dissipation rate χc may be defined as the scalar dissipation rate either at the stoichiometric mixture fraction Zst or at Zmax , where the maximum flame temperature occurs. The elliptic governing equations of mass and momentum, and the mixture fraction (the counterpart of Eq. (4)), written in axisymmetric cylindrical coordinates were solved. The gravitational acceleration term was maintained in the momentum equation in the x direction (vertically upward). The scalar dissipation rate at each point in the coflow flame is calculated as ∂Z 2 ∂Z 2 . + χ = 2DZ (6) ∂x ∂r The local values of the mixture viscosity µ, density ρ, and diffusion coefficient DZ are extracted interactively from the flamelet library based on the local values of the mixture fraction Z and the scalar dissipation rate χ . Once the converged solution of the flow field is obtained, other thermochemical properties such as temperature and species concentrations are calculated based on the flamelet library in the same way as the density and the viscosity. 2.2. 2D direct numerical simulation In the direct numerical simulation of axisymmetric coflow laminar diffusion flames, besides the governing equations of mass and momentum, the governing equations of energy and species, Eqs. (7) and (8) below, in elliptic form were also solved simultaneously, ∂T ∂T + ρu cp ρv ∂r ∂x ∂T ∂ ∂T 1 ∂ rλ + λ = r ∂r ∂r ∂x ∂x KK ∂T ∂T + Vkx ρcpk Yk Vkr − ∂r ∂x k=1

−

KK

hk Wk ωk ,

(7)

k=1

ρv

∂Yk 1 ∂ ∂Yk + ρu =− (rρYk Vkr ) ∂r ∂x r ∂r ∂ − (ρYk Vkx ) + Wk ωk , k = 1, 2, . . . , KK, ∂x

(8) where T is the mixture temperature, and Wk , cpk , and ωk are respectively the molecular weight, the specific heat at constant pressure, and the mole production rate per unit volume of the kth species. Symbols hk and Yk represent the specific enthalpy and the mass fraction of the kth species. Vkr and Vkx are the diffusion velocities of the kth species in the r and x directions, respectively and KK is the total number of species considered. The method of correction diffusion velocity [21] was employed to ensure that the net diffusion flux of all species sums to zero in both r and x directions. For example, the diffusion velocity in the x-direction is calculated as Vkx = −

1 ∂Yk + VTkx + Vcx , Dk Yk ∂x

(9)

where VTkx is the thermal diffusion velocity in the xdirection for the kth species and Vcx is the correction diffusion velocity. In this study, only the thermal diffusion velocities of H2 and H are accounted for using the expression given in [21]; i.e., the thermal diffusion velocity of all other species is set to zero. The mixture-averaged diffusion coefficient Dk is related to the binary diffusion coefficient Dkj through 1 − Xk . Dk = KK j =k Xj /Dkj

(10)

The diffusion velocity in the r-direction is calculated in a way similar to Eq. (9). Note that in the 2D direct simulation the diffusion velocities are calculated with a method identical to that used for counterflow diffusion flames. The system of governing equations is closed by the ideal gas state equation. Although the mixture fraction Z is not required in the direct simulation method, the mixture fraction field was nevertheless obtained by the two methods discussed above, i.e., either by solving its transport equation or from the definition of Bilger et al. [18], expressed in Eq. (5). In the former case, however, the diffusion coefficient DZ was now calculated directly from its definition and the associated thermal properties extracted from CHEMKIN subroutines based on the local conditions, rather than extracted from the flamelet library. The scalar dissipation rate field was also obtained in the direct simulation. The GRI-Mech 3.0 chemical reaction mechanism [22], along with its thermal and transport database,

F. Liu et al. / Combustion and Flame 144 (2006) 605–618

was used in both the counterflow flame calculations and the direct numerical simulation. Reactions and species related to NOx formation (except N2 ) were removed in the present calculations. Although this reaction mechanism is optimized for methane combustion as opposed to ethylene diffusion flames, for the purpose of the present study it is adequate to use this mechanism. In the 2D direct simulation, reaction rates and thermal and transport properties were calculated using the CHEMKIN subroutines. 2.3. Numerical method The elliptic governing equations were discretized using the control volume method. The SIMPLE numerical algorithm [23] was used to deal with pressure and velocity coupling. The diffusion and convection terms were respectively discretized by the central and upwind schemes. In the 2D direct simulation, the discretized equations of species were solved in a fully coupled fashion at each control volume to handle the stiffness of the chemistry using the method given in [24]. The discretized equations of momentum, pressure correction, and energy were solved sequentially using the standard tridiagonal matrix algorithm (TDMA).

3. Results and discussion 3.1. Counterflow diffusion flame Calculations of the counterflow diffusion flame were carried out at atmospheric pressure for a series of stretch rates a ranging from 0.5 to 1200 s−1 using the CHEMKIN codes. The air stream issued from the oxidant jet on the left and pure ethylene issued from the fuel jet on the right. The stretch rate was specified at the fuel stream. The separation distance of the two jets varied with the stretch rate. In general a larger separation distance is required to ensure that the solution has negligible gradients near both the air and fuel jets at lower stretch, since the reaction zone becomes broader. The boundary conditions at the air jet are T = 300 K, Z = 0, YO2 = 0.232, YN2 = 0.768. The inlet conditions at the fuel stream are T = 300 K, YC2 H4 = 1, Z = 1. Variations of the calculated peak flame temperature, mixture fraction at the peak flame temperature, ZTmax , and scalar dissipation rates at the peak flame temperature and at the stoichiometric mixture fraction with the stretch rate are shown in Fig. 1. At the lowest stretch rate considered in this study, the peak flame temperature is about 2336.5 K, which is only about 50 K lower than the adiabatic flame temperature of ethylene/air combustion of 2387 K (at 1 atm and initial temperature

609

300 K). At the highest stretch rate of 1200 s−1 , the peak flame temperature drops to about 1696 K. Therefore, the solutions in this flamelet library cover a wide range from near-equilibrium to near-quenching and the flamelet library is considered adequate to compute the axisymmetric coflow laminar ethylene diffusion flame. Extension of the flamelet library to an even lower stretch rate not only requires excessive computing time to obtain the counterflow diffusion flame structure, but also makes only a slight quantitative change to the results presented below. Fig. 1a shows that the mixture fraction at the peak flame temperature calculated from both methods increases significantly with the stretch rate, indicating that the peak flame temperature shifts toward the fuel-rich side as the stretch rate increases. When the mixture fraction is obtained by solving its transport equation, the resultant scalar dissipation rates at the peak flame temperature are consistently larger than those at the stoichiometric mixture fraction (Zst = 0.0634); compare the curve with open circles with that with triangles in Fig. 1b. When the definition of Bilger et al. is used for the mixture fraction, the dissipation rates at the peak temperature are even lower. At the lowest stretch rate considered, the scalar dissipation rate at ZTmax is respectively about 0.0627 s−1 when Z is solved from its transport equation and 0.0309 s−1 when Z is calculated from the definition of Bilger et al. To assist the discussion of the results presented below, the flamelet model using the library compiled based on the transport equation for Z is called Flamelet Model 1 and the one using the library based on the definition of Bilger et al. for Z is termed Flamelet Model 2. The effect of using a different choice of the characteristic dissipation rate to represent a flamelet, i.e., either χ at ZTmax or χ at Zst , was investigated using Flamelet Model 1. The difference in the results of these two different choices of χ is very similar to that between the results of Flamelet Model 1 and Flamelet Model 2 when χ at ZTmax is used as the characteristic scalar dissipation rate. Therefore, only results based on χ at ZTmax are presented below. 3.2. Structure of the coflow diffusion flame In the two-dimensional calculations using both the flamelet model and the 2D direct simulation method, the same solution domain and computational grid were used. The simulated atmospheric, axisymmetric, coflow laminar ethylene–air diffusion flame was generated with a burner in which pure ethylene flows vertically upward through a 10.9-mm inner diameter tube and the air flows from the annular region between the fuel tube and a 90-mm inner diameter concentric tube. The wall thickness of the fuel tube is 0.55 mm. The computational domain used in the calculations is

610

F. Liu et al. / Combustion and Flame 144 (2006) 605–618

Fig. 1. Variations of the peak flame temperature, the mixture fraction at the peak flame temperature, and scalar dissipation rates at the peak flame temperature and the stoichiometric mixture with the stretch rate in the counterflow flame.

Fig. 2. Schematic of the computational domain.

schematically shown in Fig. 2. The inlet conditions at the air stream are u¯ = 77.62 cm/s, v = 0 cm/s, T = 300 K, Z = 0, YO2 = 0.232, YN2 = 0.768. The inlet conditions at the fuel stream are u¯ = 3.53 cm/s, v = 0 cm/s, T = 300 K, YC2 H4 = 1, Z = 1. The inlet velocity profiles of the fuel and air streams are respectively parabolic and uniform. Along the centerline of

r = 0 and the outer boundary of r = 4.5 cm, zero gradient is applied to all variables except the radial velocity component v, which is set to 0. At the top boundary of x = 13.12 cm, zero gradient condition is assumed for all the variables. The computational domain of 13.12 cm × 4.5 cm was divided into 228 × 95 control volumes. Nonuniform grid division was used with finer grids placed in the near burner exit region in the x-direction and between 0 and 1.2 cm in the r-direction. The smallest grid sizes in the x- and rdirections were respectively 0.25 and 0.165 mm. The computational mesh was confirmed to be sufficiently fine to ensure that the numerical results presented below are grid independent. Fig. 3 displays the mixture fraction distributions from the direct simulation method and the two flamelet models. As mentioned earlier, the mixture fraction is not required in the 2D direct simulation calculation. In this case, the mixture fraction was calculated after the iteration converged by either solving the transport equation or simply from the definition of Bilger et al. [18]. Also plotted in Fig. 3 are the contours of the stoichiometric value. The black solid lines are based on the transport equation for the mixture fraction. The black dashed line, Fig. 3a, is from the definition of Bilger et al. [18]. Overall the mixture fraction fields are qualitatively similar. The two definitions of the mixture fraction for bridging the

F. Liu et al. / Combustion and Flame 144 (2006) 605–618

counterflow diffusion flame solutions in the physical space and the mixture fraction space do not cause significant differences in the mixture fraction field, Figs. 3b and 3c. On the other hand, the flame height (defined as the axial location where the mixture fraction is stoichiometric) is significantly dependent on the way the mixture fraction is calculated in the 2D direct simulation, Fig. 3a. Clearly the mixture fraction definition of Bilger et al. leads to better agreement in the flame height between the flamelet model and the 2D direct simulation. Fig. 4 compares the temperature distributions in the coflow laminar ethylene diffusion flame computed using the direct simulation method and the two flamelet models. Although the two flamelet models predict rather significantly higher temperatures in the centerline regions beyond about x = 2 cm, the overall shapes of the temperature distribution are qualitatively similar with the peak temperatures occur in the centerline region in all three results. Compared to the temperature distribution of the direct simulation method in the centerline region between x = 3 and 4 cm, results of the Flamelet Model 2 are slightly better than those from the Flamelet Model 1. The centerline variations of temperature, mixture fraction, and scalar dissipation rate are compared in Fig. 5. The mixture fraction and the scalar dissipation rate in the 2D simulation are obtained by solving the transport equation of Z. The centerline distributions of the mixture fraction and the scalar dissipation rate calculated by the two flamelet models are almost the same. The two flamelet models consistently predict a higher flame temperature along the flame centerline than the direct simulation, especially Flamelet Model 1. The peak temperature in the two flamelet models reaches the maximum value of 2336.5 K in the flamelet library at about x = 4.4 cm, where the temperature from the direct simulation method is only about 2160 K. The higher temperatures predicted by Flamelet Model 1 in the centerline region between x = 2 and 4 cm is due to the higher dissipation rates associated with this model, Fig. 1b. It is seen from Fig. 1b that for a given characteristic scalar dissipation rate, Flamelet Model 1 corresponds to a flamelet of lower stretch rate that has a higher peak flame temperature, Fig. 1a. Distributions of the scalar dissipation rate along the flame centerline from the 2D simulation and the flamelet models are very similar: it first rises rapidly to a maximum at about x = 0.5 cm and then decays monotonically. Note that the scalar dissipation rate decays to 0.0627 s−1 (the lowest value in the flamelet library) at about x = 2.3 cm in Flamelet Model 1 and 0.0309 s−1 at about x = 2.7 cm in Flamelet Model 2, which implies that from these locations downstream the flame approaches its “equilibrium” structure (the flamelet solution of a = 0.5 s−1 )

611

in the centerline region and the thermochemical properties of the flame are solely defined by the mixture fraction in the flamelet model. As a consequence of the rapid approach to equilibrium in the flamelet methodology, the peak temperature predicted by the two flamelet models is essentially the adiabatic flame temperature of ethylene. On the other hand, the substantially lower peak flame temperature than the adiabatic value in the 2D direct simulation implies that the equilibrium state cannot be reached in the coflow diffusion flame investigated here due to finite residence time and multidimensional effects. Distributions of mass fractions of two major species H2 O and CO2 are respectively compared in Figs. 6 and 7. It is evident that the flamelet model predicts somewhat higher H2 O mass fractions and significantly lower CO2 mass fractions than the direct numerical simulation. This trend is in agreement with that observed by Nishioka et al. [15], who found that the 2D direct simulation yields lower H2 O concentration and higher CO2 concentration than in results of counterflow flames, even at locations close to the burner exit. Besides the rather significant differences in their peak values, especially for CO2 , there exists a qualitative difference in the distributions of these two major species. The 2D direct simulation predicts that H2 O mass fraction peaks low in the flame in a region below x = 2 cm and around a narrow radial region around r = 0.5 cm, Fig. 6a, while the mass fraction of CO2 peaks in the centerline region between about x = 5 and 7 cm, Fig. 7a. In other words, the mass fractions of CO2 and H2 O peak in very different regions in the flame in the 2D direct simulation. These are the essential features of CO2 and H2 O concentration distributions in coflow hydrocarbon flames calculated using detailed chemistry as discussed in detail by Smooke et al. [14] in their numerical study of an axisymmetric coflow laminar methane–air diffusion flame. In contrast, the flamelet model predicts a quite different flame structure in terms of the distributions of CO2 and H2 O; i.e., they peak in the same regions in the flame (centerline) and the distributions correlate quite well with the temperature. Although the results of the two flamelet models are somewhat quantitatively different, they are qualitatively similar. The coflow flame structure predicted by the two flamelet models resembles that expected from a fast chemistry model such as the equilibrium model or the one-step irreversible reaction, i.e., the Burke–Schumann flame sheet model, in the sense that the concentrations of CO2 and H2 O and the temperature peak in the same region in the flame. To illustrate why the flamelet model predicts such a flame structure, the distributions of CO2 , H2 O, CO, and temperature in three representative counterflow flames of a = 0.5, 100, and 1000 s−1 are shown in

612

F. Liu et al. / Combustion and Flame 144 (2006) 605–618

Fig. 3. Comparison of the mixture fraction fields from the 2D direct simulation and the two flamelet models. The black solid lines are the stoichiometric mixture fraction contour from solving the transport equation of Z. The black dashed line shows the stoichiometric value of Z from the definition of Bilger et al. [18].

Fig. 4. Comparison of the temperature distributions in the coflow flame.

F. Liu et al. / Combustion and Flame 144 (2006) 605–618

613

Fig. 5. Distributions of temperature, mixture fraction, and scalar dissipation rate along the centerline of the coflow flame.

Fig. 8 in physical space. It is obvious that the concentrations of CO2 and H2 O and temperature always peak in the same location in counterflow flames either close to equilibrium at low stretches or close to extinction at high stretches. Given the methodology of the flamelet model, this method should always predict a flame structure similar to the counterflow one when applied to multidimensional flames; i.e., it is unable to predict the coflow flame structure in which the concentrations of CO2 and H2 O peak in different regions. It is also noticed that the peak values of CO2 and H2 O mass fraction predicted by the two flamelet models shown in Figs. 6 and 7 are those in the near equilibrium counterflow flame, Fig. 8a. Nevertheless, results shown in Figs. 6–8 once again demonstrate that multidimensional diffusion flames have features different from those observed in counterflow flames. Such a different structure between multidimensional and counterflow diffusion flames is not totally unexpected based on the following two considerations. First, the structure of a diffusion flame not only is controlled by chemical kinetics but also is strongly influenced by the flow field. Second, multidimensional effects (convection and diffusion in directions other than that normal to the flame) are absent in counterflow flames. It is also very important to realize that in the downstream region of a coflow diffusion flame the boundary conditions associated with a flamelet normal to the stoichiometric surface are not those of the pure oxidizer or fuel streams specified in the genera-

tion of the flamelet library, since in this region a local flamelet is formed by some fuel-lean and fuel-rich combustion products. These observations highlight the fundamental drawbacks of the flamelet methodology when applied to multidimensional diffusion flames. As mentioned earlier, preferential diffusion of species was accounted for in both the 2D direct simulation and the flamelet model (in the generation of flamelet library in counterflow configuration) calculations. Therefore, the preferential diffusion is not expected to be the cause for the significant discrepancies between the results of the 2D direct simulation and the flamelet model. However, to quantify the role of preferential diffusion in the calculated coflow flame structure, an extra run of the 2D direct simulation was conducted with the assumption of unity Lewis for all species. To our surprise, the preferential diffusion indeed has a rather significant impact on the calculated coflow flame structure. Under the unity Lewis number assumption, the peak flame temperature in the 2D direct simulation increases from 2188.2 to 2249.8 K, which is still significantly below the adiabatic flame temperature. While the peak mass fraction of H2 O is almost unaffected (increases from 0.08 to 0.087), the centerline region of H2 O concentration is significantly enhanced. When preferential diffusion is neglected, the predicted peak mass fraction of CO2 is lowered from 0.206 to 0.175. Even so, the discrepancies between the direct simulation results and those of the flamelet model still remain substantial.

614

F. Liu et al. / Combustion and Flame 144 (2006) 605–618

Fig. 6. Comparison of the distributions of mass fraction of H2 O in the coflow flame.

Fig. 7. Comparison of the distributions of mass fraction of CO2 in the coflow flame.

F. Liu et al. / Combustion and Flame 144 (2006) 605–618

615

Fig. 8. Distributions of mass fractions of CO2 , H2 O, and CO and temperature in three representative counterflow diffusion flames at a = 0.5, 100, and 1000 s−1 , respectively.

Distributions of mass fractions of CO are compared in Fig. 9. There are again rather significant quantitative differences between the results of the 2D direction simulation and the two flamelet models. The mass fraction of CO peaks lower on the flame centerline region in the flamelet model prediction than in the direct simulation, i.e., about x = 2.5 cm compared with about x = 3.75 cm. The distributions of H2 calculated by the direct simulation method and the two flamelet models (not shown) are qualitatively similar to those of CO shown in Fig. 9. The last species compared in this paper is the mass fraction of OH shown in Fig. 10. Again, the flamelet model predicts the distribution of OH concentration in qualitative agreement with the direction simulation. Results of the two flamelet models are again quantitatively similar. However, the OH concentrations in the coflow flame predicted by the flamelet model are consistently higher and spread much further downstream than those by the direct simulation. The comparison between the results of the 2D direct simulation and the flamelet model shows that there are substantial differences between them. While the effect of flow is partially built in the flamelet model to describe the degree of nonequilibrium, the present results indicate that the scalar dissipation rate

is insufficient to quantify the nonequilibrium structure of a diffusion flame in an axisymmetric coflow configuration. 3.3. Views on the failure of the flamelet model The present study demonstrates that there are substantial differences between the results of the flamelet model and the 2D direct simulation, especially in the downstream region beyond about x = 2 cm. These differences have also been previously observed by Smooke et al. [14] and Haworth et al. [25]. In fact, it has been made clear that the stretched laminar steadystate diffusion flamelet model predicts an overly rapid approach to equilibrium in the far downstream field in the modeling of turbulent diffusion jet flames [25, 26] due to the rapid decay of the scalar dissipation rate. The cause for this severe drawback of the steadystate flamelet model has been believed to be the neglect of the transient term in the flamelet equation based on a Lagrangian point of view in which the flamelet is convected from the burner rim downstream [25,26]. Through a reexamination of the flamelet theory, it is suggested in this study that the failure of the flamelet model in the calculation of 2D coflow diffusion flames (laminar or turbulent) in the down-

616

F. Liu et al. / Combustion and Flame 144 (2006) 605–618

Fig. 9. Comparison of the distributions of mass fraction of CO in the coflow flame.

Fig. 10. Comparison of the distributions of mass fraction of OH in the coflow flame.

F. Liu et al. / Combustion and Flame 144 (2006) 605–618

stream region is caused by the following two fundamental drawbacks: (i) the neglect of multidimensional convection and diffusion in directions other than that normal to the local stoichiometric surface, and (ii) inadequate use of pure oxidizer and pure fuel boundary conditions to characterize the local flamelet structure. It is important to reiterate that the flamelet formulation reduces to a one-dimensional structure only when the scalar dissipation rate is large. As the scalar dissipation rate χ becomes sufficiently small, which occurs in the far field of jet diffusion flames due to small mixture fraction gradients, the neglect of the R terms in Eqs. (1) and (2), which represent multidimensional convection and diffusion in directions other than that normal to the surface of stoichiometric mixture, is no longer justified. In other words, effects of multidimensional convection and diffusion become significant and have to be acknowledged in this situation. However, as noted by Bray and Peters [16], consideration of these higher order effects reintroduce the flow field effects into the description of the flamelet structure and obviate its supposed advantage. As revealed in the present results and in the study of Smooke et al. [14], some aspects of the structure of a coflow diffusion flame predicted by the flamelet model are qualitatively different from that predicted by the 2D direct simulation. The improvement made by recently proposed transient stretched laminar flamelet model [25,26] over the steady-state flamelet model in jet diffusion flame calculations is entirely through an ad hoc perturbation to slow down the approach to equilibrium. This remedy, however, seems unable to overcome the fundamental drawback of the flamelet methodology in predicting the correct structure of a coflow diffusion flame revealed in Figs. 6 and 7.

4. Conclusions A numerical study of an axisymmetric coflow laminar ethylene–air diffusion flame was conducted using detailed chemistry and thermal and transport properties and two approaches: the two-dimensional direct simulation method and the flamelet model. Although ambiguity exists when the flamelet library is generated from counterflow diffusion flame solutions computed in physical space, the flamelet model results of the coflow flame investigated in this study depend only somewhat quantitatively, but not qualitatively, on how the mixture fraction is defined or how the characteristic scalar dissipation rate is chosen. The structure of the coflow flame computed by the flamelet model is in overall qualitative agreement with that from the direct simulation. However, there are substantial differences between the results of the flamelet model

617

and the direct simulation method even for temperature and major species CO2 and H2 O. In addition, the flamelet model fails to predict a distinct feature of the coflow flame structure: the peak concentrations of CO2 and H2 O occur in different regions in the flame. The flamelet model predicts a rapid approach to equilibrium in the downstream region while the direct simulation predicts significant departure from equilibrium even at far downstream locations. An analysis of the flamelet model theory indicated that the failure of the flamelet model at far downstream locations is due to neglect of multidimensional convection and diffusion and inadequate use of boundary conditions for local flamelets in the downstream region. The present study demonstrates that the flamelet model is incapable of adequately incorporating the flow effect through the scalar dissipation rate. The flamelet model should be used cautiously when applied to calculate multidimensional diffusion flames.

References [1] N. Peters, Prog. Energy Combust. Sci. 10 (1984) 319– 339. [2] K. Seshadri, F. Mauß, N. Peters, J. Warnatz, Proc. Combust. Inst. 23 (1990) 559–566. [3] C.-P. Chou, J.-Y. Chen, C.G. Yam, K.D. Marx, Combust. Flame 114 (1998) 420–435. [4] M. Balthasar, A. Heyl, F. Mauß, F. Schmitt, H. Bockhorn, Proc. Combust. Inst. 26 (1996) 2369–2377. [5] B. Rogg, F. Behrendt, J. Warnatz, Proc. Combust. Inst. 21 (1986) 1533–1541. [6] D.C. Haworth, M.C. Drake, R.J. Blint, Combust. Sci. Technol. 60 (1988) 287–318. [7] D. Lentini, Combust. Sci. Technol. 100 (1994) 95–122. [8] J.-Y. Chen, W.-C. Chang, Proc. Combust. Inst. 26 (1996) 2207–2214. [9] J.P.H. Sanders, J.-Y. Chen, I. Gökalp, Combust. Flame 111 (1997) 1–15. [10] X.S. Bai, M. Balthasar, F. Mauss, L. Fuchs, Proc. Combust. Inst. 27 (1998) 1623–1630. [11] J.X. Wen, L.Y. Huang, Fire Safety J. 34 (2000) 1–24. [12] M. Chen, M. Herrmann, N. Peters, Proc. Combust. Inst. 28 (2000) 167–174. [13] A. Heyl, H. Bockhorn, Chemosphere 42 (2001) 449– 462. [14] M.D. Smooke, Y. Xu, R.M. Zurn, P. Lin, J.H. Frank, M.B. Long, Proc. Combust. Inst. 24 (1992) 813–821. [15] M. Nishioka, Y. Takemoto, H. Yamashita, T. Takeno, Proc. Combust. Inst. 26 (1996) 1071–1077. [16] K.N.C. Bray, N. Peters, in: P.A. Libby, F.A. Williams (Eds.), Turbulent Reacting Flows, Academic Press, New York, 1994, p. 63. [17] F. Mauß, D. Keller, N. Peters, Proc. Combust. Inst. 23 (1990) 693–698. [18] R.W. Bilger, S.H. Stårner, R.J. Kee, Combust. Flame 80 (1990) 135–149. [19] M.D. Smooke, I.K. Puri, K. Seshadri, Proc. Combust. Inst. 21 (1986) 1783–1792.

618

F. Liu et al. / Combustion and Flame 144 (2006) 605–618

[20] H. Pitsch, N. Peters, Combust. Flame 114 (1998) 26– 40. [21] R.J. Kee, G. Dixon-Lewis, J. Warnatz, M.E. Coltrin, J.A. Miller, A Fortran Computer Code Package for the Evaluation of Gas-Phase, Multicomponent Transport Properties, Report SAND 86-8246, Sandia National Laboratories, 1986. [22] G.P. Smith, D.M. Golden, M. Frenklach, N.W. Moriarty, B. Eiteneer, M. Goldenberg, C.T. Bowman, R.K. Hanson, S. Song, W.C. Gardiner Jr., V.V. Lissian-

[23] [24] [25] [26]

ski, Z. Qin, available at: http://www.me.berkeley.edu/ gri_mech/. S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, New York, 1980. Z. Liu, C. Liao, C. Liu, S. McCormick, AIAA Paper 95-0205, 1995. D.C. Haworth, M.C. Drake, S.B. Pope, R.J. Blint, Proc. Combust. Inst. 22 (1988) 589–597. H. Pitsch, M. Chen, N. Peters, Proc. Combust. Inst. 27 (1998) 1057–1064.