Numerical Modeling of NOx Formation in Turbulent Flames Using a

Article pubs.acs.org/EF

Numerical Modeling of NOx Formation in Turbulent Flames Using a Kinetic Post-processing Technique A. Cuoci,* A. Frassoldati, A. Stagni, T. Faravelli, E. Ranzi, and G. Buzzi-Ferraris Department of Chemistry, Materials, and Chemical Engineering “G. Natta”, Politecnico di Milano, Piazza Leonardo da Vinci 32, 20133 Milano, Italy ABSTRACT: The numerical prediction of formation of nitrogen oxides (NOx) in combustion devices (domestic and industrial burners and furnaces, internal combustion engines, etc.) is gaining rising importance because of strong limitations imposed by the legislation. Such predictions usually require the adoption of detailed kinetic mechanisms, because the numerical modeling of NOx can be accurate only considering all of the relevant chemical paths (thermal, prompt, N2O, fuel NOx, NNH, etc.). Despite the continuous increase of computational resources, the direct coupling of detailed kinetics and complex computational fluid dynamics (CFD) is very expensive in terms of central processing unit (CPU) time, and alternative approaches must be taken into account. In this work, the prediction of NOx in combustion systems is performed using a kinetic post-processing technique, which is based on the assumption that pollutants, such as NOx, only marginally affect the main combustion process (because they are usually present in very small amounts). In particular, the proposed numerical methodology, which can be applied to turbulent flames (premixed and diffusive) in arbitrarily complex geometries, transforms the original CFD simulation of the combustion device under investigation in a network of perfectly stirred reactors at a fixed temperature, which is solved by employing a very detailed kinetic scheme, accounting for the formation of NOx. The main novelty with respect to similar approaches documented in the literature is represented by the possibility to solve networks with a very large number of reactors (of the order of ∼100 000), also in the presence of very complex fluid dynamics (recirculations, swirled flows, multiphase flows, etc.). For this purpose, a robust and accurate numerical algorithm was developed, to reduce the CPU time, without sacrificing the reliability of the numerical predictions. The kinetic post-processing technique was validated on several turbulent jet flames (fed with syngas, hydrogen, and methane). The comparison to the experimental measurements is quite good, demonstrating the feasibility and reliability of the proposed approach. Then, the formation of NOx was successfully modeled in a small-scale combustor, operating in moderate to intensive low oxygen dilution (MILD) conditions.

1. INTRODUCTION Combustion devices, such as industrial and domestic burners, gas turbines, industrial furnaces, etc., need to respect always more stringent limitations concerning the emissions of pollutant species (especially NOx and CO). This explains the increasing demand for computational tools able to characterize the combustion systems, not only in terms of temperature and flow fields but also in terms of pollutant species, whose chemistry is in most cases quite complex. To be used for industrial applications, these tools have to be very accurate [pollutants are often present in parts per million (ppm)] and reliable. This goal requires the direct coupling of computational fluid dynamics (CFD) with a detailed description of the kinetics. However, even with the continuous increase of computer power and speed, the application of a detailed kinetic mechanism (DKM) in complex CFD is still too central processing unit (CPU)-consuming, especially when considering the typical dimensions of the computational grids used for industrial applications. In particular, the computational cost increases with the number of cells of the computational grid and more rapidly with the number of reacting species. Moreover, the turbulent flow of most practical combustion devices is characterized by strong interactions between fluid mixing and chemical reactions. Therefore, the prediction of pollutant species in complex, turbulent CFD is still a challenging, expensive task. In most cases, simplified approaches, specifically conceived for each class of pollutant © 2012 American Chemical Society

species, have to be taken into account. If we restrict our analysis to Reynolds averaged Navier−Stokes (RANS) models, two main approaches can be identified:1 (1) the turbulent flow is accurately described, but a reduced chemistry model is adopted; (2) the flow description is simplified, but a DKM is kept. The first modeling approach has been and is again widely used, according to several techniques. As an example, for premixed flames, the eddy break-up (EBU)2 and the Bray− Moss−Libby (BML)3 models are often adopted, while eddy dissipation (ED), eddy dissipation/finite rate (ED/FR), and eddy dissipation concept (EDC) models,4 infinitely fast chemistry, or equilibrium state are assumed in most cases for non-premixed flames. When transported probability density function (PDF) methods are applied, even if the chemical source term is closed, very simplified kinetic schemes have to be considered, to limit the CPU time.5 The second approach (on which this work is based) keeps an accurate description of the kinetics but simplifies the flow representation. This approach is very attractive, especially if the main interest is the prediction of pollutant species, whose characteristic chemical times are large, such as nitrogen oxides (NOx). Indeed, pollutant species only marginally affect the Received: October 18, 2012 Revised: December 27, 2012 Published: December 28, 2012 1104

dx.doi.org/10.1021/ef3016987 | Energy Fuels 2013, 27, 1104−1122

Energy & Fuels

Article

A similar approach was followed by Frassoldati et al.,15 who studied the formation of NOx in a swirl burner fed with natural gas and experimentally studied within the German TECFLAM cooperation.16 On the basis of the temperature and flow fields obtained from a conventional RANS simulation with reduced kinetics, a network of PSRs was built by grouping together the cells with similar temperature and composition (i.e., equivalent from a kinetic point of view). The network was then iteratively solved using a detailed kinetic scheme. Several networks were studied, with an increasing number of reactors, up to ∼300, and the effects of the grouping procedure on the predictions of NOx were discussed. The same technique was also applied for studying cases of industrial interest, such as ethylene furnaces operating in low-NOx regimes.17,18 Cuoci et al.19 applied the same technique to perform a kinetic analysis about the formation of NOx in turbulent jet flames fed with syngas and doped with various amounts of NH3. The effects of turbulent fluctuations of temperature were accounted for through the introduction of the so-called kinetic equivalent temperature, evaluated on the basis of the temperature variance, as predicted by the CFD simulation. Moreover, some improvements in the numerical procedure were introduced, to make the solution of the reactor network faster and more accurate. The same technique was applied with satisfactory results to model NOx formation in an advanced low-NOx combustor, experimentally studied in the frame of the recent European Program efficient and environmentally friendly aero-engine (EEFAE) for the component validator for environmentally friendly aero-engine (CLEAN)20 and new aero-engine core (NEWAC)21 demonstrators. In particular, it was possible to show the capabilities of the proposed procedure in very complex geometries and nonconventional operating conditions. Fichet et al.1 applied a very similar technique to post-process a gas turbine flame fed with natural gas, operating at the pressure of 15 bar, with the aim to predict the NOx emissions. Great attention was devoted to the splitting criteria to minimize the number of equivalent reactors needed for the construction of the network. Contrary to previous approaches described above, the reactor temperature was not assumed fixed from CFD results, but the energy balance was solved during the postprocessing phase. The reactor network was solved using an iterative approach around the PSR code developed by Kee et al.22 The NOx emissions predicted using a network with about 400 reactors were found in good agreement with the experimental data, confirming again the great potential of this technique to predict pollutant emissions. More recently, Monaghan et al.23 used a network of PSRs to study the pathways (thermal, prompt, N2O, and NO2) by which NOx are formed in the Sandia D piloted methane−air diffusion flame.24 The network was built by first dividing the CFD domain in several zones, on the basis of values of temperature, mixture fraction, and axial position. Within each zone, appropriate reactor criteria were applied to group neighboring computational cells. The resulting network, consisting of 1200 PSRs, was solved using a detailed kinetic scheme, containing 103 species and 582 reactions. Even if the effects of turbulent fluctuations were neglected during the postprocessing phase, the numerical results showed a satisfactory agreement with the experimental data. In this work, we describe and apply a new KPP framework, which can be applied to turbulent flows in complex geometries to predict the formation of pollutant species, such as NOx. On

main combustion process and, consequently, do not significantly influence the overall temperature and flow fields. Consequently, a kinetic post-processing (KPP) technique can be conveniently applied in three steps: (1) the temperature and flow fields are modeled using a standard RANS code and a simplified kinetic mechanism; (2) then, appropriate algorithms are applied on the CFD fields to create an ensemble of connected zones, which can be viewed as a network of perfectly stirred reactors (PSRs); and (3) finally, the reactor network is solved using a detailed kinetic scheme, to predict the pollutant species with a high degree of accuracy. The KPP technique, also known as reactor network analysis (RNA), has already been employed by various authors to postprocess CFD results and evaluate the formation of pollutants, using DKM for various applications using a different level of description and various numerical methodologies. Ehrhardt et al.6 were the first to propose and apply this postprocessing technique. In their work, on the basis of the computed three-dimensional (3D) results for flow, temperature, and stoichiometric fields, the volume of the combustor was reduced to a simplified network of ideal PSRs or plug flow reactors. Then, within each reactor, a detailed kinetic model was used to predict the concentrations of the minor pollutant species and especially NOx. It was found that it is possible to obtain fairly grid-independent results for NOx concentrations with few reactors. However, the proposed zone model was restricted only to flows without recycling zones and with downstream convection much larger than upstream diffusion. In 2001, Faravelli et al.7 extended this approach by applying the so-called simplified fluid dynamics by ideal reactor network (SFIRN) technique for the prediction of NOx in a 75 MWe furnace. The main difference with respect to the approach proposed in ref 6 was represented by the possibility to model flows with complex fluid dynamics with strong recirculations (i.e., strong interaction loops between the reactors), to apply more complex kinetic schemes and to extend the procedure also to liquid fuels. Falcitelli et al.8 proposed a general algorithm to construct the reactor network, without the limitations reported above, which was applied to practical industrial systems, such as glass melting furnaces and pilot- and full-scale boilers.9,10 The calculated emissions of NOx were found in very good agreement with the experimental data, demonstrating the feasibility of the RNA for the design of practical combustion systems. The same procedure was adopted by Mancini et al.11 to calculate NOx emissions from combustion of natural gas in moderate to intensive low oxygen dilution (MILD) conditions. Skjøt-Rasmussen et al.12 applied a computational tool to post-process the data extracted from RANS simulations performed with reduced chemistry. Instead of simplifying the flow field, in their post-processing technique, Skjøt-Rasmussen et al.12 retained all of the individual cells in the computational domain, which were treated as PSRs, with fixed temperature and mass flow rates, as predicted by the CFD simulation. The reactors were modeled using DKM and the effects of turbulent fluctuations were accounted for by splitting each reactor in reactive and inert parts, following the approach proposed by Gran and Magnussen.13,14 The proposed approach was successfully applied to a cylindrical furnace, which was postprocessed using a reactor network of 30 000 reactors with a very large kinetic scheme, involving 159 species and 773 reactions. 1105


Energy & Fuels

Article

built from the computational mesh adopted for the CFD simulation (pre-processing step); the number of reactors can be arbitrarily large, up to the original number of computational cells. (3) Turbulence− chemistry corrections: to take into account the effects of the turbulent fluctuations of the temperature, the reaction rates are corrected following the procedure described in section 2.3. (4) Reactor network solution: the reactor network is solved using a detailed gas phase mechanism. (5) Visualization: at the end of the calculations, the results are written into a new post-processing file to be visualized using the appropriate graphical post-processor (in most cases, the same code used for the CFD simulation). Phases 1 and 5 are specific to the CFD and visualization software and must be adapted accordingly to the specific codes. However, the core of the post-processing procedure (i.e., phases 2−4) is completely independent of the CFD code, which means that the KPP can be used with various CFD codes without modifications (see also refs 20 and 21, where proprietary CFD codes were employed). 2.1. Data Extraction. In this phase, all of the relevant data are extracted from the CFD simulation. (1) geometry: mesh topology and effective volumes of computational cells; (2) mass flow rates (convection and diffusion) between the computational cells; and (3) pressure, mean temperature and its variance, and composition of each computational cell. No information about the geometry of each computational cell is required by the KPP. The only relevant data are the connections between the cells. This means that the KPP can be applied without any difference to one-dimensional (1D), two-dimensional (2D), and 3D cases, with arbitrarily complex structured and unstructured meshes. The composition extracted in this phase is used only as a first guess in the post-processing procedure. 2.2. Grouping of Cells. The temperature and composition fields obtained through the CFD code allow for the identification of the critical zones in the computational domain, i.e., the specific regions where large temperature and/or composition gradients are present. In such zones, it is more convenient to retain the original detail of the CFD grid. However, there are usually large volumes of the system that are less critical from a kinetic point of view (especially in terms of NOx), in particular cold and/or nonreactive zones. As a consequence, the detail of the grid can be locally reduced by clustering and combining several cells into a single equivalent reactor. Of course, the lumped volume is simply the sum of the volumes of the grouped cells. The original CFD mesh is thus transformed into a network of PSRs, where the links between the different reactors simply combine and reflect the original flow field, as evaluated by the CFD code. According to the desired level of accuracy, the total number of equivalent reactors can be modified; as an example, the original ∼105−106 CFD cells can be conveniently grouped into ∼103−104 equivalent reactors, maintaining a reasonable description of the flame structure. This technique makes feasible to couple complex, large CFD meshes to very DKM. The mesh-coarsening algorithm was designed to prevent possible dangerous situations, such as the creation of geometrical irregularities and/or non-smooth transition between zones with very different volumes. The interlinking flows are evaluated on the basis of the convection/diffusion mass flow rates exchanged between the cells belonging to the different cells. Temperature and initial compositions in the equivalent reactors are the volume-averaged values of the combined cells. Different clustering levels are sequentially adopted, and calculations are iteratively performed by increasing the number of reactors up to the final convergence (i.e., up to the point where a further increase in the dimensions of the reactor network makes no significant difference in the final solution). The accuracy and convergence of the solution, together with the effect of the coarsening of the mesh, need to be carefully monitored, as better discussed in section 3. 2.3. Turbulence−Chemistry Corrections. As already reported, the KPP uses the temperature field as obtained by the CFD computations. A fixed average temperature is assumed in each equivalent reactor, and the rates of all of the reactions involved in the kinetic scheme are evaluated accordingly. However, in turbulent combustion conditions, these reaction rates cannot simply be

the basis of a proper CFD simulation of the combustion system under investigation, the KPP framework builds and solves a network of PSRs with detailed kinetics. The main novelty with respect to the previous works summarized above is represented by the possibility to solve a network with a very large number of equivalent reactors (on the order of ∼100 000) and very detailed kinetic schemes (∼100 species), also in presence of very complex fluid dynamics (recirculations, swirled flows, etc.). For this purpose, a very robust and accurate numerical algorithm was developed, to reduce the CPU time, without sacrificing the reliability of the numerical predictions. The paper is organized as reported in the following. In section 2, we present the KPP methodology, by describing all of the different phases required for the prediction of pollutant species. Then, in section 3, the numerical algorithm adopted to solve the reactor network is discussed. In section 4, the whole methodology is validated by performing numerical simulations of several laboratory turbulent flames and comparing the numerical predictions to the experimental measurements. Then, NOx formation in a small-scale combustor is analyzed and discussed in section 5. Finally, some conclusions and perspectives are drawn for the improvement of the KPP framework and its extension to even more complex combustion systems.

2. KPP FRAMEWORK The KPP framework works on the basis of the CFD simulation of the combustion system, where a simplified kinetic approach was adopted. From the simplified kinetic approach, we intend either the use of global mechanisms (in the context of ED, ED/FR, EDC, EBU, and BML models), the application of the equilibrium, or the steady-state laminar flamelet approaches.25 Then, an equivalent network of ideal PSRs is derived on the basis of the flow, thermal, and composition fields resulting from the CFD. The temperature in each reactor and the mass flow rates (convection and diffusion) between the reactors are kept constant during the post-processing calculations, according to the CFD simulation. On the contrary, the concentration of the species in the different reactors is calculated by solving the corresponding mass balance equations using a detailed kinetic scheme. The presence of possible recycle streams, recirculation zones, etc. is automatically handled during the construction and the solution of the reactor network. It is thus evident that the accurate CFD solution becomes a fundamental prerequisite for a correct estimation of pollutant formation. The KPP applies two major simplifications, making this technique feasible and attractive over the direct coupling of a detailed kinetic scheme inside the CFD code.19 (1) Reduction of the number of equations: the solution of the CFD code provides the detailed flow, composition, and temperature fields, and this information allows for critical and noncritical zones in the overall reacting system to be identified. Therefore, several adjacent and very similar cells are lumped (or grouped) into single equivalent reactors. This means that the description detail can be reduced in several regions without significantly affecting the results. This leads to a strong reduction of the dimensions of the overall system. (2) Reduction of the nonlinearities: a second way of making the numerical computations more stable and viable is to define a fixed temperature in each reactor. This reduces the large nonlinearity of the original system, which is mainly associated with the reaction rates (through the kinetic constants, depending upon the exponential of the temperature because of the Arrhenius law). The post-processing procedure developed in this work is based on five distinct phases, which are summarized in the following. (1) Data extraction: all of the relevant data (topology of the mesh, volumes of the computational cells, temperature, mass flow rates, mass fractions, etc.) are extracted from the CFD simulations and written into a file to be used by the KPP code. (2) Grouping of cells: a network of PSRs is 1106


Energy & Fuels

Article

⎛μ ⎞ ∂ ε (ρσ 2) + ∇(ρvσ 2) = ∇⎜ t ∇σ 2⎟ + Cgμt (∇T̅ )2 − Cdρ σ 2 ∂t κ ⎝ σt ⎠

calculated as a function of the mean temperature and composition, mainly because of the highly nonlinear dependence of reaction rates upon temperature. Temperature dependence of rate constants is usually described via the modified Arrhenius equation

⎛ Ej ⎞ kj(T ) = A j T βj exp⎜− ⎟ ⎝ RT ⎠

(3) where the constants σt, Cg, and Cd take the values 0.85, 2.86, and 2.0, respectively. The major simplification in the proposed approach consists of neglecting the effects of composition fluctuations on the average reaction rates. This choice is necessary to avoid recalculating the correction coefficients at every iteration during the post-processing phase. However, in most cases, the weight of composition fluctuations on the mean reaction rate is very small if compared to the temperature fluctuations. In addition, to account for nonperfect mixing, each equivalent reactor of the network can be separated in two parts: one part is inert with respect to chemical reactions, while the other part is modeled as a PSR, and the chemical reactions are allowed to progress according to the DKM.12 The fractional split between the reactive and inert part is determined from γ, the volume fraction of the fine structures in the reactor, which is calculated as13,14

(1)

where the subscript j refers to the jth generic reaction. By inspection of eq 1, it is evident that especially temperature fluctuations have a significant effect on the average rates of reactions with high activation energy. This effect is very important for the reactions involved in NOx formation and needs to be taken into account. The average rate value (which accounts for temperature fluctuations over time) is very different from the reaction rate calculated at the mean temperature. This difference obviously increases for large temperature fluctuations and for reactions with high activation energies. Hence, if time-averaged composition and temperature are employed to predict the mean NOx formation rate, significant errors will result. To tackle this problem with reasonable computational efforts, temperature fluctuations are taken into account by considering the PDF describing the time variation. The PDF method has proven very useful in the theoretical description of turbulent flow.26 In particular, a single PDF in terms of a normalized temperature is used to correct the reaction rates during the post-processing technique kj̅ =

∫T

Tmax

min

kj(T )P(T )dT = Cjkj(T̅ )

⎛ νε ⎞1/4 γ = 2.13⎜ 2 ⎟ ⎝κ ⎠

(4)

where ν is the laminar kinematic viscosity of the fluid, ε is the rate of dissipation of turbulent kinetic energy, and κ is the turbulent kinetic energy. This means that the effective volume V of each equivalent reactor is equal to γVg, where Vg is the geometrical volume. 2.4. Reactor Network Solution. In this step, the equations corresponding to the reactor network built in the previous step are solved to predict the concentrations for all of the species in the detailed kinetic scheme in every equivalent reactor. As better explained in section 3.1, the system of equations to be solved, corresponding to the steady-state mass balances for all of the species, is strongly nonlinear (because of the reaction terms) and can reach very large dimensions. Therefore, in this work, great attention was devoted to the formulation of an accurate and robust numerical algorithm. In particular, in contrast to previous works,1,12 the network is not solved using a sequential procedure, in which the reactors are solved in sequence, but a global Newton’s algorithm is applied for the following reasons:28 (1) from one side, it is necessary to speed-up the convergence rate, which is very low in the case of direct substitutions, especially when the number of reactors is large and the original flow field is characterized by strong recirculations (such as in industrial furnaces); (2) at the same time, a global method not only increases the computational efficiency but also ensures the complete convergence to the solution. It is important to notice that one of the aims of the KPP procedure is the evaluation of pollutant species, which are often present in very small amounts (on the order of a few ppm); in these conditions, it is very important to converge to the solution with great accuracy (i.e., with very small tolerances), and this goal cannot be reached using direct substitutions. Actually, the numerical procedure combines different techniques to obtain the final solution, because the global Newton’s method (or modified Newton’s methods) can be successfully applied only if the first-guess solution is close to the real solution. As better explained in section 3, the new numerical procedure described in the present work introduces some substantial improvements with respect to the methodolgies adopted in previous works: (1) the robustness was strongly improved, thanks to a new global strategy to solve the global nonlinear system (on the basis of Armijio’s rule); (2) a higher level of accuracy can be achieved through the adoption of the well-known, widely tested multifrontal massively parallel sparse direct solver (MUMPS) and library information system (LIS) libraries for the solution of large, sparse linear systems; and (3) larger reactor networks with more detailed kinetic schemes can now be solved thanks to a more efficient management of the memory. Indeed, the previous numerical methodology adopted by the same authors19,23 had some deficiencies in the management of the memory, which

(2)

where k̅j is the mean turbulent rate of reaction j, kj(T) is the corresponding instantaneous value, P(T) is the PDF of the temperature, and Cj is a correction coefficient, introduced for convenience. Equation 2 must be integrated in every reactor for all of the reactions. Because the temperature is kept constant during the post-processing phase, this operation must be performed only once, at the beginning of the calculations, by storing the correction coefficients Cj. For the PDF of the temperature, the limits of integration are determined from the minimum and maximum values of the temperature in the combustion solution. P(T) is assumed to be a two-moment β function, as appropriate for combustion calculations.25 Because the β function requires that the independent variable assume values between 0 and 1, the temperature field is normalized for the application of eq 2. Figure 1 shows the correction coefficient Cj as a function of temperature fluctuations for two different activation energies (30 000

Figure 1. Correction coefficient of the rate constants versus the intensity of temperature fluctuations at average temperatures of 1500 and 2000 K, with an activation energy of (left) 30 000 cal/mol and (right) 70 000 cal/mol.

and 70 000 cal/mol) at the mean temperatures of 1500 and 2000 K. As expected, the coefficient increases with an increase of the temperature variance and the activation energy. The evaluation of the correction coefficient requires the temperature variance σ2, which is calculated in the CFD code by solving the corresponding transport equation27 1107


Energy & Fuels

Article

limited the maximum number of the equations to be solved through the global procedure. 2.5. Visualization. Visualization is the final step in the postprocessing phase. In this step, the concentrations of all of the species in the detailed kinetic scheme are mapped from the reactor network solved by the KPP to the original CFD grid, to be used in the appropriate graphical post-processor (in most cases, the same CFD code used for the original RANS simulation). There is no feed of data from the KPP into the CFD code; the procedure works exclusively as a post-processing phase. Moreover, additional utilities for kinetic investigations are available in this step, with the possibility to plot the maps of reaction and formation rates and sensitivity coefficients for selected species and reactions.

species contained in the DKM. In particular, it is pretty evident that the nonlinear contribution in these equations is only in the reaction term, because the convection/diffusion terms are linear with respect to the mass fractions. To make the description of the numerical methodology adopted in this work easier, eq 8 can be written also in the following form, by grouping the NC equations for each reactor: g(k)(ω) = −C(k)ω(i) + f(k) + R (k)(ω(k))

(9)

where ω(i) is the (NR × 1) vector of all of the mass fractions of the ith species in all of the reactors, ω(k) is the (NC × 1) vector of mass fractions of all of the species in the kth reactor, and C(k) is the (NC × NR) matrix accounting for the convection/ diffusion terms of the kth reactor. Actually, all of the rows of this matrix are equal, because the convection/diffusion mass flow rates are the same for all of the species in the same reactor. g(k), f(k), and R(k) are the (NC × 1) subvectors of g, f, and R, respectively, in which only the values corresponding to the kth reactor are considered. The Jacobian matrix J associated with the system of eq 8 presents some important features that can be conveniently exploited by the numerical method used to obtain the solution. In particular, J is a very sparse matrix, which can be decomposed in (NR × NR) square blocks J(k,p) with dimensions equal to (NC × NC), with most of them being null. In particular, the NR block matrices along the main diagonal have the following form:

3. REACTOR NETWORK SOLUTION 3.1. Governing Equations. Each equivalent reactor is assumed as a PSR at a fixed temperature (as predicted by the CFD simulation), exchanging mass (by convection and diffusion) with the neighbor cells and with the external environment (if the reactors belong to a boundary of the computational domain). The mass balance for the ith species (i = 1, ..., NC, where NC is the number of species) in the kth reactor (k = 1, ..., NR, where NR is the number of reactors) can be easily written by imposing that the output mass flow rate of the ith species is equal to the sum of its input mass flow rate and the production term because of the chemical reactions. NR

Mkωki =

∑ Dkjωji + fki

+ ΩkiVk (5)

j=1

In the equation reported above, Mk is the output mass flow rate (convective and diffusive) from the kth reactor and Dkj is the mass flow rate (convective and diffusive) from the jth reactor to the kth reactor. f ki is the mass flow rate of species i entering the kth reactor from the external environment (i.e., through the boundaries of the computational domain) and/or from an evaporating liquid phase or from the devolatilization of a solidphase particle. Ωki is the formation rate of the ith species, and Vk is the effective volume (see section 2.3) of the kth reactor. The convection/diffusion terms (Mk and Dkj) are extracted from the CFD simulation and are kept constant during the post-processing phase. Because of conservation of mass, they are related by the following constraint: NR

Mk =

j=1

J(k , k) = −Mk I + J(Rk)

J(Rk) =

J(k , p) = Dkp I

j=1

+ ΩkiVk = 0 (7)

or in a more compact form g(ω) = −Cω + f + R(ω) = 0

i , j = 1, ..., NC

(11)

k , p = 1, ..., NR

k≠p

(12)

In Figure 2, an example of a typical Jacobian structure for a reactor network with six reactors and three species is reported. It is pretty evident that full blocks are present only on the main diagonal; off-diagonal blocks, representative of convection/ diffusion exchanges, are only diagonal. By inspection of eq 9, it is evident that the reaction term is local; i.e., it involves only the unknowns in the reactor under investigation. On the contrary, the convection/diffusion terms are linear with respect the unknowns. 3.2. Newton’s Global Method. As previously explained, the nonlinear system is solved using the global Newton’s method, to ensure the accuracy needed to correctly predict chemical species present in very small amounts. The first guess solution corresponds to the composition calculated by the CFD simulation. Newton’s method proceeds iteratively, by estimating at each iteration the Newton direction d, which is calculated by solving the following system of linear equations:

NR

∑ (Dkj − Mkδkj)ωji + fki

∂R ki ∂ωkj

The blocks outside the main diagonal are given by

Obviously most of the Dkj coefficients are equal to zero, because each reactor is connected only with a small number of different reactors. Equation 5 can be rewritten as gki =

(10)

where I is the identity matrix and the matrix is defined on the basis of the reaction term in the kth reactor.

(6)

i=1

k = 1, ..., NR

JR(k)

NC

∑ Dkj + ∑ fki

k = 1, ..., NR

(8)

where C is a (NE × NE) sparse matrix accounting for the convection and diffusion terms, f is the vector accounting for input mass flow rates from the external environment, and R is the vector corresponding to the reaction term. NE is the total number of equations; i.e., NE = NC × NR. δkj is Kronecker’s delta. Equation 8 corresponds to a large system of nonlinear equations, whose unknowns are the mass fractions of all of the

J(n) d = −g(n)

(13) (n)

in which the global Jacobian matrix J and the residuals g(n) are evaluated at the iteration n. Then, the current solution ω(n) is corrected by adding the Newton correction (or the Newton 1108


Energy & Fuels

Article

3.3. Direct Substitutions (Sequential Approach). The global Newton’s method is usually not robust enough to solve the system when the CFD composition is used as a first-guess solution.35 Therefore, it is more convenient to approach a better estimation of the solution by iteratively solving the sequence of individual reactors with successive substitutions (helped by suitable acceleration techniques), which have the advantage of being rather inexpensive in terms of CPU time and take the whole system closer to the solution. This means that each reactor is solved using a local Newton’s method (i.e., a conventional Newton’s method applied to the single reactor). To improve the robustness of this approach, especially in the first iteration, a false transient method is used to solve the single reactors. The nonlinear system (eq 9) is transformed into an ordinary differential equation (ODE) system of equations by adding the accumulation term mk Figure 2. Example of the sparsity pattern of the Jacobian matrix associated with eq 8.

dt

= −C(k)ω(i) + f(k) + R (k)(ω(k))

k = 1, ..., NR

(15)

where mk is the total mass in the kth reactor. The composition at the previous iteration is assumed as the initial condition of this ODE system. The integration is performed over a long time interval, to ensure steady-state conditions (i.e., the accumulation term is sufficiently low). The local ODE systems are very stiff and are solved using the BzzOde solver,36,37 which was specifically conceived for reacting systems with DKM. The BzzOde solver was successfully applied in the past for the numerical modeling of laminar flames,38−40 showing excellent performances in terms of accuracy and robustness. The sequential approach is applied for several iterations and stopped in one of these conditions: (1) when the residuals of all of the equations reach sufficiently low values (which means that we are close to the real solution of the global algebraic system) or (2) when the convergence rate is too slow or, even worse, if the residuals are increasing (this is typical of complex flows, with strong recirculations). In the first case, the global Newton’s method can then be applied with success in most cases. On the contrary, in the case of slow convergence (or bad convergence), the successive application of Newton’s method can fail, because the first guess solution is still too bad. In such situations, a global time-step method is applied to improve the first guess solution, as better explained in the next section. 3.4. Global Time Stepping. When complex flows, with strong recirculations, are investigated and the number of reactors in the network is quite large, the sequential approach (i.e., direct substitutions) is not enough to reduce the residuals of all of the equations to sufficiently small values to successfully apply the global Newton’s method. In such a case, a global time-stepping procedure must be taken into account. The original global nonlinear system is converted into a global ODE system. This conversion is accomplished by appending the accumulation term in each balance equation, obtaining

step), which is evaluated as a scalar multiple of the Newton direction d. ω(n + 1) = ω(n) + λd

dω(k)

(14)

The positive scalar λ is chosen to ensure a sufficient decrease of the norm of the new residuals ∥g(n+1)∥2 with respect to the old residuals ∥g(n)∥2. Its value is calculated using the so-called Armijo rule (additional details can be found in refs 29 and 30. Figure 3 summarizes the whole procedure.

Figure 3. Global Newton’s method.

The Jacobian matrix is evaluated as reported in the previous paragraph. In particular, the convection/diffusion terms are fixed and kept equal to those predicted by the CFD R simulations. Only the blocks along the main diagonal J(k) have to be re-evaluated at each iteration. To increase the computational efficiency, special attention is devoted to the calculation of JR(k). In particular, the derivatives of rate equations are evaluated analytically rather than numerically. The most expensive operation is the solution of the linear system (eq 13). In the current version of KPP, the solution can be obtained using either direct or iterative solvers, according to the memory available. In the first case, the solution is performed using the MUMPS 4.10 solver.31,32 In the second case, the linear system is solved through the LIS library.33 Several algorithms can be chosen by the user (as an example, the GMRES(m) algorithm restarted every 40 iterations or the BiCGSafe method), using different preconditioners.34

dω mk ki = gki = dt k = 1, ..., NR 1109

NR

∑ (Dkj − Mkδkj)ωji + fki

+ ΩkiVk

j=1

i = 1, ..., NC

(16)


Energy & Fuels

Article

with appropriate initial conditions. The time derivative is then replaced by a backward Euler approximation over the time step τ(n+1) = t(n+1) − t(n). mk

ωki(n + 1) − ωki(n) τ

(n + 1)

− gki(n + 1) = 0

i = 1, ..., NC

k = 1, ..., NR (17)

At each time step, the global Newton’s method is applied to solve the system. ⎛ mk ⎞ ⎜ (n + 1) I − J(n) ⎟Δω = g(n) ⎝τ ⎠

(18)

The equation reported above looks very similar to eq 13. The important difference is that the diagonal of the time-dependent Jacobian is approximately that of the steady-state Jacobian J(n) but with an additive term on the main diagonal. This increase in the size of the diagonal produces a better conditioned system, and the solution from the nth time step usually provides an excellent starting guess to the solution at the next n + 1 time level. If direct solvers are employed, the work per time step is similar to that of the steady-state global Newton method. On the contrary, because of the better properties of the system (eq 18) in terms of conditioning number, if iterative solvers are used, the CPU cost for solving the system (eq 18) is much smaller that the CPU time to solve the system (eq 13). In any case, the time-like continuation of the numerical solution produces an iteration strategy that is usually more robust and less sensitive to the initial starting estimate than if the system (eq 9) were solved directly.41,42 The time step is not kept constant during the time-stepping procedure. In particular, it is reduced if negative mass fractions are found after the solution of the system (eq 18); otherwise, it is gradually increased, to improve the convergence toward the steady-state solution. Because we are interested in solving neither the transient problem all of the way to steady state nor the accuracy of the transient path, the time-stepping procedure is not particularly sophisticated. The method is relatively inexpensive per each time step and employs the same Newton’s algorithm used in the steady-state solution. The transient solution is used only to estimate the first guess solution for Newton’s algorithm but does not provide an accurate solution of the transient process. Generally speaking, the time-stepping method is quite robust but computationally inefficient if compared to the global Newton’s method, which is more efficient but has less desirable convergence properties. Therefore, after solving for a specified number of time steps, the KPP attempts to solve the steady-state problem by Newton’s method. If the steady-state solution fails, the solvers goes back to the time-stepping procedure. Obviously, the better the initial estimate of the solution, the less likely the KPP will have to resort to time stepping. The whole numerical procedure is summarized in Figure 4. 3.5. Convergence. Great attention is devoted to check the convergence of the numerical algorithm. In the case of direct substitutions, the numerical convergence is generally controlled by the typical normalized error sum of squares ⎛ ω (n + 1) − ω (n) ⎞2 ki ⎟ e = ∑ ∑ ⎜⎜ ki ⎟ (n) ω ⎠ k=1 j=1 ⎝ ki NR

Figure 4. Schematic representation of the numerical procedure for the solution of a large nonlinear system (adapted from ref 28).

where the superscript n refers to the iteration number. The request for e to be less than a user-defined precision εg is a necessary but not sufficient condition. Actually, when working with direct substitutions, a small e value may just be the result of convergence failure rather than the approach to the numerical solution. Conversely, when a Newton method is adopted, if e < εg, then the method is converging to the solution within the required precision (also see ref 28). Our numerical procedure checks also for the residuals of all of the equations of the global nonlinear system to be less than a user-defined precision εl. By doing so, the effective convergence to the solution of the nonlinear system is extensively and consistently controlled. 3.6. Adjustments of Mass Flow Rates. To successfully solve the nonlinear system (eq 9) with a global Newton’s method, it is important to ensure that the total mass flow rates on each reactor are perfectly close to zero. This is not always verified in the CFD solution; small errors in the closure can be observed, because in most cases, the CFD codes work using segregated approaches. Before starting to solve the reactor network, the output mass flow rates from each reactor are then adjusted to ensure that, for each reactor, the sum of input flow rates perfectly matches with the sum of output mass flow rates j≠k

Mk −

∑ j = 1, NR

NC

αkjMj =

∑ fki i=1

(20)

where αkj = Ckj/Cjj defines how the output mass flow rate from the kth reactor is split between the remaining reactors of the network. If we keep the splitting coefficients αkj fixed (as predicted by the CFD calculations), eq 20 corresponds to a (NR × NR) linear system in the unknowns Mk, whose solution provides the adjusted output mass flow rates from all reactors. This adjustment is performed only once, at the beginning of the calculations, because as described above, the mass flow rates are kept constant during the post-processing procedure. In most

NC

(19) 1110


Energy & Fuels

Article

Figure 5. Comparison between experimental measurements (symbols)48,49 and numerical predictions (lines) of the axial velocity: (a) axial direction and (b) radial direction at several axial positions.

equation for ε, to take into account the effect of vortex stretching on the scalar dissipation and, therefore, on effective viscosity. In all of the flames, the steady laminar flamelet model (SLFM) was employed to model the turbulent combustion.25 This model is widely used for the CFD simulation of industrial combustion systems of large dimensions, because its requests in terms of CPU time are very low, especially if compared to alternative approaches [such as EDC, transported PDFs, conditional moment closure (CMC), etc.]. The radiative heat transfer was modeled according to the suggestions proposed in the TNF website (http://www.sandia. gov/TNF/). In particular, the optically thin approximation was adopted (although it is known that this approach slightly overestimates radiation because it neglects reabsorption). According to this model, the radiative loss rate per unit volume can be calculated as46

cases, these adjustments are very small (less than 0.1% of the CFD values).

4. VALIDATION Several turbulent flames were simulated and studied, to validate the proposed KPP procedure. In particular, the attention was focused on three axisymmetric turbulent jet flames studied in the context of the framework of the International Workshop on Measurements and Computation of Turbulent Nonpremixed Flames (http://www.sandia.gov/TNF/), for which a very detailed and accurate experimental characterization is available in terms of velocity, temperature, main species, and NOx. All of the flames were simulated with the commercial CFD code Ansys FLUENT 13.27 A 2D steady-state simulation of the physical domain was considered because of the cylindrical symmetry of the system. The simulations were performed with conformal grids and rectangular cells. The grid points were not evenly spaced, but they were more dense near the axis of the system (in the radial direction) and near the nozzle of the burner (in the axial direction), to improve the prediction of the spreading rate of the jet, which affects the shape of the flame. The width and length of the grid were chosen to avoid any effect of boundaries on the flame, which were assumed unconfined. For the spatial resolution, the first-order upwind scheme was adopted. The segregated implicit solver was used, with the semi-implicit method for pressure-linked equations (SIMPLE) procedure for pressure−velocity coupling.43 The pressure staggering options (PRESTO!) algorithm was adopted for pressure interpolation.27 Turbulence was modeled via the RANS approach.44 The flames under consideration are simple jet diffusion flames, in which the central jet spreads because of the turbulent diffusion, allowing the fuel to mix with the surrounding oxidizer and react. No swirl motion is introduced, and no bluff-body is employed to stabilize the flame. As a consequence, the fluid dynamics of the system is quite simple and no recirculation zones or other complex flow patterns arise. Therefore, the standard κ−ε model seems to be a good choice for modeling the turbulence, without resorting to more complex, expensive models. However, it is well-recognized that the standard κ−ε model poorly predicts the velocity field in axisymmetric round jets; in particular, it tends to overestimate the spreading rate (i.e., the decay rate).45 Also, for the flames under investigation, the standard κ−ε model overpredicts the diffusion of the central jet and predicts the axial velocity on the centerline lower than the velocity actually measured. To avoid such inaccuracies, the κ−ε model was applied with the correction proposed by Pope,45 which consists of adding a source term in the transport

Q = −4σaP(T 4 − Tenv 4)

(21)

where Tenv is the environment temperature and σ is the Stefan− Boltzmann constant. The aP coefficient accounts only for the contributions of H2O, CO2, CO, and CH4, according to the following expression: aP = pH O aP,H2O + pCO aP,CO2 + pCO aP,CO + pCH aP,CH4 2

2

4

(22)

where pi and aP,i are the partial pressure and the Planck mean absorption coefficients, respectively, of species i, which are derived from calculations performed by the RADCAL software.47 4.1. CO/H2 Flame. The first flame, experimentally investigated in refs 48 and 49 is fed with a fuel stream whose composition is 40% CO, 30% H2, and 30% N2 (by volume). The burner has a central duct constructed from straight tubing with squared-off ends, with an internal diameter of 4.58 mm. The thick wall of the tubing (∼0.88 mm) creates a small recirculation zone that helps the flame stabilization. The computational grid was refined in this zone to better resolve the details of the near-nozzle flow. The central fuel jet mixes with the co-flow air stream, resulting in a turbulent unconfined diffusion flame. The jet fuel velocity is ∼76 m/s; the co-flow air velocity is ∼0.70 m/s; and both of the streams are at 292 K. The resulting Reynolds number is ∼16 700. Experimental results include axial and radial profiles of mean and root-meansquare (rms) values of temperatures and major species concentrations as well as velocity statistics and Reynolds stresses. Radial profiles of NO and OH radical concentrations are also available at different locations. 1111


Energy & Fuels

Article

Figure 6. Comparison between experimental measurements (symbols)48,49 and numerical predictions (lines) of the mixture fraction: (a) axial direction and (b) radial direction at several axial positions.

Figure 7. Comparison between experimental measurements (symbols)48,49 and numerical predictions (lines) of the temperature: (a) axial direction and (b) radial direction at several axial positions.

which accounts for the chemistry of formation of NOx.52−54 Figure 8 reports the predicted maps of NO and NO2, with

The PolimiH2CO1201 kinetic scheme (freely available at http://creckmodeling.chem.polimi.it/ in CHEMKIN format) was used to build the flamelet look-up table.50,51 This mechanism was extensively discussed and validated in a wide range of conditions.52 Figure 5 reports a comparison between the numerical predictions and the experimental measurements of axial velocity along the centerline and in the radial direction at several axial locations (x/d = 20, 40, and 60, where x is the distance from the fuel nozzle and d is the fuel nozzle internal diameter). The velocity predictions can be considered satisfactory, especially if compared to the simulations performed with the standard κ−ε model (not reported here). However, it is still clear that the CFD simulation slightly overestimates the jet decay rate. The turbulent kinetic energy is well-captured along the axis. The calculated and measured profiles of the mixture fraction and temperature (together with their rms) are compared in Figures 6 and 7, respectively, along the axial and radial directions. The numerical predictions are quite satisfactory; in particular, the temperature peak along the centerline is wellpredicted, together with the flame length. Only a small overprediction of the temperature in the tail of the flame can be observed. We performed a sensitivity analysis, which confirmed that the predicted temperature and mixture fraction profiles are nearly insensitive to the number of computational cells, to the numerical schemes employed, and to the kinetic mechanism used to generate the flamelet library (also see ref 19). The overall quality of the CFD simulation was considered satisfactory to perform the KPP analysis. Of course, the correct prediction of the velocity and temperature fields is a necessary prerequisite for the application of the KPP. Thus, moving from the CFD simulation described above, the KPP was applied using the PolimiH2CO1201NOx detailed kinetic scheme,

Figure 8. Two-dimensional maps of calculated NOx, accounting (left) or neglecting (right) the effects of turbulent fluctuations of the temperature.

(left) and without (right) accounting for the turbulent fluctuations of the temperature. The significant formation of NO2 in the post-flame zone can be clearly observed. The effect of turbulent fluctuations on the formation of NOx is quite important; as an example, if they were neglected, the predicted amount of NO would be less than 50% with respect to the amount calculated by accounting for the fluctuations of the temperature. 1112


Energy & Fuels

Article

Figure 9. Comparison between experimental measurements (symbols)49 and numerical predictions (lines) of the NO mass fraction: (a) axial direction and (b) radial direction at several axial positions.

Figure 10. NO mass fraction at different distances from the fuel nozzle. Comparison between single shot measurements (symbols)49 and numerical predictions (lines).

Figure 9 reports a comparison between the experimental measurements and the predicted mass fractions of NO along the axial and radial directions, together with the effect of turbulent fluctuations on the formation of NO. The agreement is satisfactory, even if some discrepancies can be observed. The shape of the radial NO profile in the flame is correctly reproduced at various distances from the fuel nozzle and is in good agreement with the measurements along the axis, with the exception of the peak temperature region, where the model tends to slightly overestimate the NO mass fraction. Figure 10 shows the single-shot NO measurements as a function of the mixture fraction in the radial direction at different axial locations, together with a comparison to numerical predictions (with and without turbulent fluctuations). It is evident that the predicted NO mass fraction is in good agreement with the measured values, especially at large x/ d ratios, while close to the nozzle NO tends to be underestimated. It is interesting to note that the predicted NO profile obtained when suppressing the effect of turbulent fluctuations (i.e., Cj = 1 in eq 2) lies largely below the lower boundary of the scatter plot. Figure 11 analyses the different contribution to the formation of NO in the flame. In particular, the KPP procedure was

Figure 11. Predicted NO profiles along the centerline for the syngas flame. The profiles were obtained by suppressing the thermal, N2O, or NNH mechanisms.

applied not only with the NOx complete kinetic scheme but also by suppressing the thermal path, the N2O path, or the NNH path. As expected, the NNH path has a negligible role on the formation of NO in the investigated flame. On the contrary, the thermal and especially N2O paths play a more important role; in particular, the N2O mechanism seems to contribute to about 65% to the formation of NO, while the thermal 1113


Energy & Fuels

Article

Figure 12. Clustering procedure details: (a) predicted NO profiles along the axis for different numbers of reactors and (b) peak values of NO, NO2, and N2O as a function of the number of reactors in the network.

Figure 13. Comparison between experimental measurements (symbols)56,57 and numerical predictions (lines) of the mixture fraction and temperature: (a) mixture fraction along the axial direction, (b) mixture fraction along the radial direction, (c) temperature along the axial direction, and (d) temperature along the radial direction.

kinetic point of view). Figure 12 reports some results referring to the clustering procedure applied to the syngas flame under investigation. In particular, Figure 12a shows the predicted NO mass fraction profiles along the centerline calculated using several networks with a smaller number of reactors with respect to the original CFD mesh. Figure 12b reports the calculated peak values of NO, NO2, and N2O as a function of the number of reactors of the network. It is quite evident that, to obtain reliable and accurate predictions of NOx formation in this particular flame, a network with 2000 reactors is enough. This means that the dimensions of the post-processing problem are only 10% of the original CFD simulation, with an important gain in terms of the CPU time. 4.2. H2/He Flame. The second flame used for the validation of the KPP procedure was experimentally studied by Barlow and Carter.56,57 The experimental setup consists of a vertical jet of pure H2 burning in a coaxial air stream. The inner diameter of the tube is 3.75 mm, while the outer is 4.8 mm. The air velocity has been fixed at 1 m/s, while the fuel jet enters at 296 m/s. Both air and fuel streams have a temperature equal to 295 K. A fully developed turbulent pipe flow can be assumed at the nozzle exit. The Reynolds number (evaluated with respect to

mechanism seems to contribute to only the remaining 35%. The NO formation through N2O is initiated by the third-order reaction N2 + O + M = N2O + M, which is followed by several reactions with O, OH, and H radicals, ultimately leading to the formation of NO and N2. The selectivity of this process is ruled by the local temperature and composition of the flame. The significant role played by the N2O mechanism in syngas combustion is a consequence of the significantly enhanced production of O radicals.55 The thermal mechanism is initiated and controlled by the Zeldovich mechanism through the O + N2 = NO + N reaction, which is followed by the N + O2 = NO + O and N + OH = NO + H reactions. The thermal mechanism requires high temperatures to be activated, and this explains its minor role in the present flame, whose peak temperature is relatively low (∼1920 K). The results reported above were obtained by performing the post-processing calculations on a very large reactor network, corresponding to the original CFD mesh, without any simplification (20 440 cells). However, as reported in section 2.2, the KPP is able to build automatically simplified reactor networks from the CFD topology, by grouping the cells with similar temperature and composition (i.e., equivalent from a 1114


Energy & Fuels

Article

Figure 14. Comparison between experimental measurements (symbols)56,57 and numerical predictions (lines) of the NO mass fraction: (a) axial direction and (b) radial direction at several axial positions.

the fuel nozzle) is 10 000, and the visible flame length of the flame was 675 mm. Also in this case, the chemistry was described through the PolimiH2CO1201 kinetic scheme.52 Figure 13 reports the comparison between measured and calculated mixture fractions and temperature profiles, both along the centerline and in the radial direction at several axial locations. The agreement is quite satisfactory, particularly for the mixture fraction; the shape of axial and radial profiles is correctly captured, both close to the fuel nozzle and at larger distances. The temperature profile along the centerline is well-predicted, together with its rms. As already reported, the correct prediction of the thermal field is a fundamental prerequisite for the application of the KPP technique. The KPP was applied to the CFD results described above. Also in this case, the NOx formation was described using the PolimiH2CO1201NOx kinetic scheme.50,52 The axial and radial profiles of the calculated NO mass fraction are reported in Figure 14 and compared to the experimental measurements. The calculations were performed on several networks, with an increasing number of reactors, up to the convergence, which was reached with about 4000 reactors (while the original CFD simulation was performed on a mesh with 22 000 cells). The calculated profiles are not completely satisfactory if compared to experiments. Even if the peak value along the axis is wellcaptured, the simulations largely underpredict the formation of NO in the first part of the flame, in the region between the fuel nozzle and x/d = 120. A possible explanation could be related to the wrong prediction of O2 and OH mass fraction profiles of the original CFD simulation, reported in Figure 15. From these results, it seems quite evident that the experimental flame is characterized by a small lift-off, which allows O2 to penetrate in

the region close the fuel nozzle (x/d < 50) and the combustion reactions to occur (as evident from the small amount of OH at x/d = 25). The CFD simulation was not able to predict this flame lift-off. Therefore, because no O2 is available in the region x/d < 120, the predicted formation of NO is negligible here. Similar computational results were also obtained by Obieglo et al.58 using more advanced turbulent combustion models based on the scalar probability density function transport equation. The effect of turbulent fluctuations of the temperature is more significant than in the previous flame; if these fluctuations were neglected, the peak NO mass fraction along the axis would only be 20% of the experimental value, as evident from Figure 14. This result is quite simple to explain, especially if coupled with results reported in Figure 16, where the different

Figure 16. Predicted NO profiles along the centerline. The profiles were obtained by suppressing the thermal, N2O, or NNH mechanisms. The “complete” and “no NNH” curves overlap.

contributions to the formation of NO are reported. Because the temperatures in the H2/He flame are quite high, the thermal path is the dominant path leading to the formation of NO, while the N2O path has a marginal role (the NNH route is overwhelmed by the thermal mechanism). Therefore, the NO chemistry is very sensitive to the temperature fluctuations, because the reactions involved in the thermal path have a very high activation energy (also see Figure 1). 4.3. CH4/H2/N2 Flame. As a third example, a simple jet flame fed with a mixture of CH4, H2, and N2 in low-velocity coflow was numerically investigated. Original measurements at the German Aerospace Center (DLR),59 which include Raman/Rayleigh scattering and planar imaging, were complemented by laser Doppler velocimetry (LDV) measurements performed at the Darmstadt University of Technology60 and Raman/Raleigh/laser-induced fluorescence (LIF) measurements at Sandia National Laboratories.61 The fuel nozzle has an inner diameter of 8.0 mm and is tapered to the thin edge.

Figure 15. Comparison between experimental measurements (symbols)56,57 and numerical predictions (lines) of O2 and OH mass fractions along the centerline. 1115


Energy & Fuels

Article

Figure 17. Comparison between experimental measurements (symbols)61 and numerical predictions (lines) of the mixture fraction and temperature: (a) mixture fraction along the axial direction, (b) mixture fraction along the radial direction, (c) temperature along the axial direction, and (d) temperature along the radial direction.

Figure 18. Comparison between experimental measurements (symbols)61 and numerical predictions (lines) of the NO mass fraction: (a) axial direction and (b) radial direction at several axial positions.

Figure 19. NO mass fraction at different distances from the fuel nozzle. Comparison between single-shot measurements (symbols)61 and numerical predictions.

truncation of the original GRI-Mech,63 with the objective of developing a smaller, detailed mechanism, to reproduce closely the main combustion characteristics predicted by the full mechanism. The mechanism involves 22 chemical species and 111 reactions and was considered sufficiently accurate to model the methane flame under investigation.

The fuel stream (22.1% CH4, 33.2% H2, and 44.7% N2, by volume) is fed at 42.2 m/s for a resulting Reynolds number equal to 15 200. Experimental measurements, along the centerline and at several axial locations in the radial direction, include velocity, temperature, main species, and NO. The DRM22 kinetic scheme was adopted to perform the CFD simulations of the flame.62 This scheme was developed by 1116


Energy & Fuels

Article

Figure 20. Ansys FLUENT NOx post-processor predictions: (a) H2/CO flame, (b) H2/He flame, and (c) CH4/H2/N2 flame. For each flame, three lines are compared: results obtained without considering the effects of turbulence on the NOx chemistry (no fluctuations), considering only the fluctuations of the temperature (T fluctuations), or considering only the fluctuations of mixture fractions (z fluctuations).

N2O intermediate model is activated. NOx transport equations are solved in a post-processing phase based on a given flow field and combustion solution. The main difference with respect to the KPP procedure is that the mass balance equations are solved only for the NOx species reported above, while the remaining species are kept equal to the CFD values; i.e., they are not updated during the post-processing calculations. This approach is very fast, because only a small number of equations has to be solved, but it can lead to large inaccuracy, as demonstrated in the following. It is thus evident that an accurate combustion solution becomes a prerequisite of the NOx prediction. The simulations were performed using the following configuration: thermal, prompt, and N2O pathways activated, instantaneous thermal [O] and [OH] models, and quasi-steady NO2 model. The effects of turbulent fluctuations are taken into account by considering a β-PDF distribution with 10 PDF points for the temperature or the mixture fraction. The numerical results are reported in Figure 20, where the calculated NO mass fractions along the axis are compared to the experimental values for all three flames. In particular, for each flame, three curves are reported, depending upon the turbulence−chemistry interaction model, which was applied. It is evident that the predictions are not satisfactory; in particular, for the H2/CO flame, calculated NOx is largely overestimated, while for the H2/He flame, a reasonable agreement can be obtained only by disabling the turbulent correction, which is not easy to explain from a physical point of view. The NOx predictions in the CH4/H2/N2 flame are largely underestimated with respect to the experimental data.

Figure 17 reports the comparison between measured and calculated mixture fractions and temperature profiles, both along the centerline and in the radial direction at several axial locations. The agreement is satisfactory; the axial and radial profiles are predicted with reasonable accuracy, even if some discrepancies in terms of the shape can be observed, especially for the radial temperature profiles. The CFD results were post-processed using the GRI30 kinetic mechanism (53 species and 325 reactions).63 Again, several networks were calculated, from 200 to 3000 reactors, which were found enough for the correct description of NOx formation in the flame under investigation. NO mass fraction profiles were found in satisfactory agreement with the experimental measurements, in both axial and radial directions (Figure 18). However, the calculated peak value along the centerline (76 ppm) overestimates the measured value (60 ppm) by about 20%. Because the temperature characterizing the present flame are not so high (the peak value is 1850 K), the effect of turbulent fluctuations is less relevant than in the H2/He flame. However, if they were neglected, the numerical predictions would be affected by significant errors, as better evident from Figure 19. Here, single-shot NO measurements are reported as a function of the mixture fraction in the radial direction at x/d = 40 and 60, together with a comparison to numerical predictions (with and without turbulent fluctuations). Only if turbulent fluctuations of the temperature are taken into account are the experimental data well-captured. 4.4. Comparisons to the FLUENT Post-processor. To better show the reliability of the proposed post-processing approach and its advantages with respect to similar techniques used in the literature, the formation of NOx was calculated also using the kinetic post-processor available in Ansys FLUENT 13.27 To predict NOx emissions, Ansys FLUENT 13 solves a transport equation for the NO concentration, together with additional transport equations for intermediate species (HCN and/or NH3) if fuel NOx sources are present and for N2O if the

5. NOX FORMATION IN A SMALL-SCALE COMBUSTOR After the validation of the KPP procedure described in the previous section, a small-scale MILD combustor was taken under investigation. Here, we report only a short description of the experimental setup, which is described with more detail in 1117


Energy & Fuels

Article

condition at the outlet. The turbulent flow was modeled using the RANS approach and the standard κ−ε model. The EDC model was employed to account for the interactions between turbulence and chemistry.4 The DRM22 kinetic scheme was employed to describe the combustion of methane.62 The in situ adaptive tabulation (ISAT) method was used to reduce the computational cost of the chemical source term integration for the EDC model.66 Thermal radiation was modeled using the discrete ordinates method, and the radiative properties of the gas mixture were accounted for through the weighted sum of gray gases model.27 Figure 22 reports the maps of the velocity, temperature, and main species as predicted by the CFD code. As better explained in ref 65, the internal recirculation zone, originated by the large momentum of the air jet, is characterized by a large extension, up to 80% of the combustor length. As a consequence, the incoming fuel mixes together with the hot combustion products, which come back from the top of the burner. As required by MILD combustion, because the momentum of the fuel stream is much lower than the air stream, the fuel entrains the flow of recirculated combustion products and mixes with them before igniting in a diluted environment. Figure 23 shows the comparison between calculated and measured profiles of the temperature and O2, CO2, and CO mole fractions along the centerline of the combustion chamber and at several axial locations. The predicted mole fractions of O2 and CO2 along the axis and in the radial direction are quite satisfactory, but a significant overprediction of the CO peak mole fraction is quite evident, even if the location of this peak agrees well with the experimental measurements. Moreover, the predicted profiles of CO in the radial direction are closer to the centerline and narrower than experimental data, which spread more significantly. The temperature predictions are not completely satisfactory, especially along the axis of the combustor. In particular, the rise of the temperature predicted by the CFD simulations appears very delayed with respect to the measurements. However, the experimental profile of the temperature is not completely consistent with the measured profile of O2 and CO2, which are correctly captured by the model. Therefore, it seems to be reasonable to consider the possibility that the experimental temperature measurements were affected by some kind of error. The radial temperature

ref 64. The combustion chamber is a quartz glass cylinder, with an internal diameter of 100 mm and a length of 340 mm (Figure 21). The fuel is pure methane at ambient temperature.

Figure 21. Schematic of the combustor (adapted from ref 64).

The air is supplied through a central nozzle with an internal diameter of 10 mm, which is surrounded by 16 small orifices of 2 mm internal diameter each (positioned on a circle with a diameter of 30 mm) for the fuel supply. The combustion air is preheated at 400 °C before entering in the combustion chamber by an electrical heating system. Experimental measurements of the temperature and mole fractions of O2, CO2, CO, and NOx are available along the centerline and in the radial direction at several axial locations. The configuration simulated in the present work (called run 2 in ref 64) corresponds to an excess air coefficient of 1.3. The inlet velocity of the fuel stream is 6.2 m/s, while the air inlet velocity is 113.2 m/s. The calculations were performed using the Ansys FLUENT 13 CFD code. Because of symmetry of the geometry, the computational domain corresponds to a section of only 11.25° of the combustion chamber. The structured mesh has about 100 000 computational cells, with a higher density in the region close to the burner. The fuel and air ducts were included in the CFD simulation to reduce the uncertainty in the assignment of the turbulent kinetic energy and its dissipation rate at the inlets. Flat velocity profiles for the fuel and air streams were imposed. The intensity of turbulence at the inlets was fixed equal to 5%, and the wall temperature was assumed equal to 1400 K, according to the suggestions reported in ref 65. The outlet pressure with zero-gauge pressure was set for the boundary

Figure 22. Calculated maps of the temperature, axial velocity, and mole fractions of OH, CO, CO2, and H2O. 1118


Energy & Fuels

Article

Figure 23. Comparison between experimental measurements (symbols)64 and numerical predictions (lines) of the temperature and O2, CO2, and CO mole fractions.

profiles in the first half of the combustion chamber show a peak overpredicting the experimental data, while for x > 150, the agreement appears quite good. The predicted profiles are very similar to those obtained by Verissimo et al.,65 who simulated the same burner (with a different excess air coefficient) using several combustion models and skeletal mechanism or DKM. None of the employed models was able to give completely satisfactory results, especially in terms of the temperature. As explained by Verissimo et al.,65 the poor prediction of the temperature field is probably due to the boundary conditions on the chamber walls, where the temperature was simply kept equal to a fixed, uniform value. Because this value is quite high (1400 K), the radiation emitted by the wall is expected to be very significant if compared to the emissions from the medium, and this may influence the temperature of the gases. The CFD results were post-processed using the KPP code. In particular, the PolimiC1C3HT1201 kinetic mechanism was adopted to predict the NOx formation.51,67 Several reactor networks were analyzed, from 500 to 50 000 reactors. It was observed that about 15 000 reactors were enough to ensure the convergence of the clustering procedure with respect to the NOx predictions. Because of the MILD conditions, the temperature fluctuations are expected to have small values. This was confirmed by the numerical simulations (not reported here), which showed the marginal role that the formation of NOx played by the fluctuations of the temperature (less than 5% on the peak values). Figure 24 reports the calculated maps of NO and N2O with a reactor network of 15 000 reactors. The concentration of NO in the second half of the combustion chamber is practically uniform (about 20 ppm), because of the strong recirculations of hot gases. The peak values are not located on the axis but in the region where the fuel and air

Figure 24. Calculated maps of mole fractions of NO and N2O, as predicted by the KPP procedure.

streams mix together and where the temperature reaches its maximum value (about 1800 K). On the contrary, the concentration of N2O is less uniform, and the peak value is located along the axis of the burner. Figure 25 shows the comparison between the calculated and measured NO mole fractions along the centerline and in the radial direction at several axial locations. The calculated values at the burner exit are in good agreement with the experimental data. However, some discrepancies can be observed. In 1119


Energy & Fuels

Article

Figure 25. Comparison between experimental measurements (symbols)64 and numerical predictions (lines) of the NO mass fraction.

the different paths in the formation of NOx. The results clearly showed the feasibility and reliability of the proposed procedure. A small-scale combustor operating in MILD conditions was then investigated. The results were satisfactory, even if some discrepancies with respect to the experimental measurements are evident. This was attributed to the poor prediction of the temperature field in some regions of the combustion device, probably because of the simplified boundary conditions applied on the walls of the cylindrical chamber. Further investigations are needed in the future to improve the accuracy of the numerical predictions. The proposed numerical methodology can be considered a useful tool for the optimal design of new combustion devices or for the evaluation of the performances of existing devices. Future extensions of this work will move in the following two main directions. (1) Improvement of the clustering procedure: the current algorithm employed for the automatic construction of reactor networks from the CFD simulation is quite simple and collects the computational cells, which are equivalent from a kinetic point of view (small differences in the temperature and composition). As evidenced by the simulations reported in sections 4 and 5, using this procedure, the number of reactors needed to obtain accurate and network-independent NOx predictions can be quite large (thousands of reactors). Results reported by Fichet et al.1 and Monagham et al.23 seem to indicate that, using smarter procedures also based on the residence time and mixture fraction, the number of reactors can be strongly reduced, without sacrificing the accuracy of NOx predictions. (2) Parallelization of the code: the main limitation of the KPP procedure is represented by the computational time, which increases very rapidly with the number of reactors. The parallelization of the code on distributed memory architectures appears as a valid solution to overcome this limitation.

particular, the model is not able to account for the formation of NOx in the proximity of the air stream nozzle (0 < x < 50 mm), where about 5 ppm of NO was measured. The predicted formation of NO occurs only at x > 50, where the temperature start to increase. The predictions in the radial direction are not completely satisfactory, especially in the first half of the combustion chamber. This poor agreement reflects the temperature profiles, which were not correctly described for the reasons previously explained. Further downstream, in the second half, the predictions are in good agreement with the measurements, in terms of both shape and absolute values. To obtain more satisfactory predictions of NOx in the present burner, it is necessary to improve the CFD simulation, especially in terms of the temperature. As already reported, the correct prediction of the temperature field is a fundamental prerequisite for the reliable KPP. This goal deserves further investigations. In particular, the boundary conditions on the walls of the combustion chamber need to be better prescribed to predict the temperature field with more accuracy.

6. CONCLUSIONS AND FUTURE WORKS In this work, we described and applied a new KPP framework to model the formation of NOx in turbulent flames and arbitrarily complex geometries. The main novelty with respect to similar approaches reported in the literature is represented by the possibility to solve very large reactor networks (with thousands of reactors) with DKM (hundreds of species) using fully coupled numerical algorithms to achieve great accuracy in the prediction of NOx concentrations. The framework was validated on several turbulent jet flames fed with syngas, hydrogen, and methane. In particular, the CFD simulations were performed using relatively small kinetic schemes, without considering the formation of NOx. The KPP procedure was then applied with different degrees of accuracy, by performing calculations on networks with an increasing number of reactors. The agreement between numerical predictions and experimental measurements was quite good, for all of the cases taken into account. Kinetic analysis was also performed to show the relative importance of

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. URL: http://creckmodeling. chem.polimi.it/. 1120


Energy & Fuels

Article

Notes

and formation of nitrogen oxide. Symp. (Int.) Combust., [Proc.] 2000, 28 (1), 303−309. (17) Barendregt, S.; Risseeuw, I.; Waterreus, F.; Frassoldati, A.; Buzzi-Ferraris, G.; Faravelli, T.; Ranzi, E.; Li, X. J.; Patel, A. R. Combustion modeling of ethylene furnaces applying ultra low NOx LSV burners. Proceedings of the 2006 AIChE Spring National Meeting; San Francisco, CA, Nov 12−17, 2006. (18) Van Goethem, M.; Risseeuw, I.; Barendregt, S.; Frassoldati, A. The design of ultra-low NOx critical furnaces. Proceedings of the 2010 AIChE Spring Meeting and 6th Global Congress on Process Safety; San Antonio, TX, March 21−25, 2010. (19) Cuoci, A.; Frassoldati, A.; Buzzi Ferraris, G.; Faravelli, T.; Ranzi, E. The ignition, combustion and flame structure of carbon monoxide/ hydrogen mixtures. Note 2: Fluid dynamics and kinetic aspects of syngas combustion. Int. J. Hydrogen Energy 2007, 32 (15), 3486−3500. (20) Frassoldati, A.; Cuoci, A.; Faravelli, T.; Ranzi, E.; Colantuoni, S.; Martino, P.; Cinque, G. Experimental and modeling study of a low NOx combustor for aero-engine turbofan. Combust. Sci. Technol. 2009, 181 (3), 483−495. (21) Frassoldati, A.; Cuoci, A.; Faravelli, T.; Ranzi, E.; Colantuoni, S.; Di Martino, P.; Cinque, G.; Kern, M.; Marinov, S.; Zarzalis, N.; Da Costa, I.; Guin, C. Fluid dynamics and detailed kinetic modeling of pollutant emissions from lean combustion systems. Proc. ASME Turbo Expo: Power Land, Sea, Air 2010, 2, 451−459. (22) Kee, R.; Rupley, F.; Miller, J. Chemkin-II: A Fortran Chemical Kinetic Package for the Analysis of Gas Phase Chemical Kinetics; Sandia National Laboratories: Livermore, CA, 1989; Sandia Report SAND898009B UC-706. (23) Monaghan, R.; Tahir, R.; Cuoci, A.; Bourque, G.; Furi, M.; Gordon, R.; Faravelli, T.; Frassoldati, A.; Curran, H. Detailed multidimensional study of pollutant formation in a methane diffusion flame. Energy Fuels 2012, 26 (3), 1598−1611. (24) Barlow, R.; Frank, J. Effects of turbulence on species mass fractions in methane/air jet flames. Proc. Combust. Inst. 1998, 27, 1087−1095. (25) Peters, N. Laminar diffusion flamelet models in non-premixed turbulent combustion. Prog. Energy Combust. Sci. 1984, 10 (3), 319− 339. (26) Veynante, D.; Vervisch, L. Turbulent combustion modeling. Prog. Energy Combust. Sci. 2002, 28 (3), 193−266. (27) Ansys, Inc. ANSYS FLUENT 13.0 Documentation; Ansys, Inc.: Canonsburg, PA, 2011. (28) Manca, D.; Buzzi-Ferraris, G.; Cuoci, A.; Frassoldati, A. The solution of very large nonlinear algebraic systems. Comput. Chem. Eng. 2009, 33 (10), 1727−1734. (29) Kelley, C.; Keyes, D. Convergence analysis of pseudo-transient continuation. SIAM J. Numer. Anal. 1998, 35 (2), 508−523. (30) Coffey, T.; Kelley, C.; Keyes, D. Pseudo-transient continuation and differential algebraic equations. SIAM J. Sci. Comput. 2003, 25 (2), 553−569. (31) Amestoy, P.; Duff, I.; Koster, J.; L’Excellent, J. A fully asynchronous multifrontal solver using distributed dynamic scheduling. SIAM J. Matrix Anal. Appl. 2001, 23 (1), 15−41. (32) Amestoy, P.; Guermouche, A.; L’Excellent, J.; Pralet, S. Hybrid scheduling for the parallel solution of linear systems. Parallel Comput. 2006, 32 (2), 136−156. (33) Nishida, A. Experience in developing an open source scalable software infrastructure in Japan. Lect. Notes Comput. Sci. 2010, 6017, 448−462. (34) Saad, Y.; Schultz, M. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 1986, 7, 856−869. (35) Keyes, D.; Smooke, M. Flame sheet starting estimates for counterflow diffusion flame problems. J. Comput. Phys. 1987, 73 (2), 267−288. (36) Buzzi-Ferraris, G.; Manca, D. BzzOde: A new C++ class for the solution of stiff and non-stiff ordinary differential equation systems. Comput. Chem. Eng. 1998, 22 (11), 1595−1621.

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Prof. Mario Costa from Instituto Superior Tecnico (Lisboa, Portugal) is acknowledged for providing experimental details about the small-scale combustor studied in Section 5. The authors also acknowledge the support of Marco Van Goethem from Technip BV and the useful suggestions of Prof. Luca Formaggia and Prof. Carlo De Falco from Politecnico di Milano. This work was partially supported by the European Union (EU) as part of the EMICOPTER Project (CS-GA2009-251798).

■

REFERENCES

(1) Fichet, V.; Kanniche, M.; Plion, P.; Gicquel, O. A reactor network model for predicting NOx emissions in gas turbines. Fuel 2010, 89 (9), 2202−2210. (2) Spalding, D. Mixing and chemical reaction in steady confined turbulent flames. Symp. (Int.) Combust., [Proc.] 1971, 13 (1), 649−657. (3) Bray, K.; Moss, J. A unified statistical model of the premixed turbulent flame. Acta Astronaut. 1977, 4 (3−4), 291−319. (4) Magnussen, B.; Hjertager, B. On mathematical modeling of turbulent combustion with special emphasis on soot formation and combustion. Symp. (Int.) Combust., [Proc.] 1977, 16 (1), 719−729. (5) Pope, S. PDF methods for turbulent reactive flows. Prog. Energy Combust. Sci. 1985, 11 (2), 119−192. (6) Ehrhardt, K.; Toqan, M.; Jansohn, P.; Teare, J.; Beer, J.; Sybon, G.; Leuckel, W. Modeling of NOx reburning in a pilot scale furnace using detailed reaction kinetics. Combust. Sci. Technol. 1998, 131 (1− 6), 131−146. (7) Faravelli, T.; Bua, L.; Frassoldati, A.; Antifora, A.; Tognotti, L.; Ranzi, E. A new procedure for predicting NOx emissions from furnaces. Comput. Chem. Eng. 2001, 25 (4−6), 613−618. (8) Falcitelli, M.; Pasini, S.; Rossi, N.; Tognotti, L. CFD + reactor network analysis: An integrated methodology for the modeling and optimization of industrial systems for energy saving and pollution reduction. Appl. Therm. Eng. 2002, 22 (8), 971−979. (9) Falcitelli, M.; Tognotti, L.; Pasini, S. An algorithm for extracting chemical reactor network models from CFD simulation of industrial combustion systems. Combust. Sci. Technol. 2002, 174 (11−12), 27− 42. (10) Falcitelli, M.; Pasini, S.; Tognotti, L. Modelling practical combustion systems and predicting NOx emissions with an integrated CFD based approach. Comput. Chem. Eng. 2002, 26 (9), 1171−1183. (11) Mancini, M.; Schwoppe, P.; Weber, R.; Orsino, S. On mathematical modelling of flameless combustion. Combust. Flame 2007, 150 (1−2), 54−59. (12) Skjoth-Rasmussen, M.; Holm-Christensen, O.; Ostberg, M.; Christensen, T.; Johannessen, T.; Jensen, A.; Glarborg, P.; Livbjerg, H. Post-processing of detailed chemical kinetic mechanisms onto CFD simulations. Comput. Chem. Eng. 2004, 28 (11), 2351−2361. (13) Gran, I.; Magnussen, B. A numerical study of a bluff-body stabilized diffusion flame. Part 1. Influence of turbulence modeling and boundary conditions. Combust. Sci. Technol. 1996, 119 (1−6), 171− 190. (14) Gran, I.; Magnussen, B. A numerical study of a bluff-body stabilized diffusion flame. Part 2. Influence of combustion modeling and finite-rate chemistry. Combust. Sci. Technol. 1996, 119 (1−6), 191−217. (15) Frassoldati, A.; Frigerio, S.; Colombo, E.; Inzoli, F.; Faravelli, T. Determination of NOx emissions from strong swirling confined flames with an integrated CFD-based procedure. Chem. Eng. Sci. 2005, 60 (11), 2851−2869. (16) Schmittel, P.; Gunther, B.; Lenze, B.; Leuckel, W.; Bockhorn, H. Turbulent swirling flames: Experimental investigation of the flow field 1121


Energy & Fuels

Article

(37) Buzzi-Ferraris, G. Bzzmath 6.0 Library; http://homes.chem. polimi.it/gbuzzi/index.htm. (38) Cuoci, A.; Mehl, M.; Buzzi-Ferraris, G.; Faravelli, T.; Manca, D.; Ranzi, E. Autoignition and burning rates of fuel droplets under microgravity. Combust. Flame 2005, 143 (3), 211−226. (39) Cuoci, A.; Frassoldati, A.; Faravelli, T.; Ranzi, E. Frequency response of counter flow diffusion flames to strain rate harmonic oscillations. Combust. Sci. Technol. 2008, 180 (5), 767−784. (40) Dixon-Lewis, G.; Marshall, P.; Ruscic, B.; Burcat, A.; Goos, E.; Cuoci, A.; Frassoldati, A.; Faravelli, T.; Glarborg, P. Inhibition of hydrogen oxidation by HBr and Br2. Combust. Flame 2012, 159 (2), 528−540. (41) Smooke, M.; Miller, J.; Kee, R. Determination of adiabatic flame speeds by boundary value methods. Combust. Sci. Technol. 1983, 34 (1−6), 79−90. (42) Grcar, J.; Kee, R.; Smooke, M.; Miller, J. A hybrid newton/timeintegration procedure for the solution of steady, laminar, onedimensional, premixed flames. Symp. (Int.) Combust., [Proc.] 1988, 21 (1), 1773−1782. (43) Patankar, S. Numerical Heat Transfer and Fluid Flow; Taylor and Francis Group: London, U.K., 1980. (44) Poinsot, T.; Veynante, D. Theoretical and Numerical Combustion, 2nd ed.; R.T. Edwards, Inc.: Flourtown, PA, 2005. (45) Pope, S. An explanation of the round jet/plane jet anomaly. AIAA J. 1978, 16 (3), 279−281. (46) Barlow, R.; Karpetis, A.; Frank, J.; Chen, J.-Y. Scalar profiles and no formation in laminar opposed-flow partially premixed methane/air flames. Combust. Flame 2001, 127 (3), 2102−2118. (47) Grosshandler, W. RADCAL: A Narrow-Band Model for Radiation Calculations in a Combustion Environment; National Institute of Standards and Technology (NIST): Gaithersburg, MD, 1993; NIST Technical Note 1402. (48) Flury, M. Experimentelle analyse der mischungstruktur in turbulenten nicht vorgemischten flammen. Ph.D. Thesis, ETH Zurich, Zurich, Switzerland, 1998. (49) Barlow, R.; Fiechtner, G.; Carter, C.; Chen, J.-Y. Experiments on the scalar structure of turbulent CO/H2/N2 jet flames. Combust. Flame 2000, 120 (4), 549−569. (50) CRECK Modeling Group. Gas Phase Combustion Schemes, Version 1201; http://creckmodeling.chem.polimi.it/kinetic.html. (51) Ranzi, E.; Frassoldati, A.; Grana, R.; Cuoci, A.; Faravelli, T.; Kelley, A.; Law, C. Hierarchical and comparative kinetic modeling of laminar flame speeds of hydrocarbon and oxygenated fuels. Prog. Energy Combust. Sci. 2012, 38 (4), 468−501. (52) Frassoldati, A.; Faravelli, T.; Ranzi, E. The ignition, combustion and flame structure of carbon monoxide/hydrogen mixtures. Note 1: Detailed kinetic modeling of syngas combustion also in presence of nitrogen compounds. Int. J. Hydrogen Energy 2007, 32 (15), 3471− 3485. (53) Faravelli, T.; Frassoldati, A.; Ranzi, E. Kinetic modeling of the interactions between NO and hydrocarbons in the oxidation of hydrocarbons at low temperatures. Combust. Flame 2003, 132 (1−2), 188−207. (54) Frassoldati, A.; Faravelli, T.; Ranzi, E. Kinetic modeling of the interactions between NO and hydrocarbons at high temperature. Combust. Flame 2003, 135 (1−2), 97−112. (55) Steele, R.; Malte, P.; Nicol, D.; Kramlich, J.; Arai, N.; Dean, A.; Glassman, I.; Reisel, J. NOx and N2O in lean-premixed jet-stirred flames. Combust. Flame 1995, 100 (3), 440−449. (56) Barlow, R.; Carter, C. Raman/Rayleigh/LIF measurements of nitric oxide formation in turbulent hydrogen jet flames. Combust. Flame 1994, 97 (3−4), 261−280. (57) Barlow, R.; Carter, C. Relationships among nitric oxide, temperature, and mixture fraction in hydrogen jet flames. Combust. Flame 1996, 104 (3), 288−299. (58) Obieglo, A.; Gass, J.; Poulikakos, D. Comparative study of modeling a hydrogen nonpremixed turbulent flame. Combust. Flame 2000, 122, 176−194.

(59) Bergmann, V.; Meier, W.; Wolff, D.; Stricker, W. Application of spontaneous Raman and Rayleigh scattering and 2D LIF for the characterization of a turbulent CH4/H2/N2 jet diffusion flame. Appl. Phys. B: Lasers Opt. 1998, 66 (4), 489−502. (60) Schneider, C.; Dreizler, A.; Janicka, J.; Hassel, E. Flow field measurements of stable and locally extinguishing hydrocarbon-fuelled jet flames. Combust. Flame 2003, 135 (1−2), 185−190. (61) Meier, W.; Barlow, R.; Chen, Y.-L.; Chen, J.-Y. Raman/ Rayleigh/LIF measurements in a turbulent CH4/H2/N2 jet diffusion flame: Experimental techniques and turbulence−chemistry interaction. Combust. Flame 2000, 123 (3), 326−343. (62) Kazakov, A.; Frenklach, M. Reduced Reaction Sets Based on GRIMech 1.2; http://www.me.berkeley.edu/drm/. (63) Smith, G.; Golden, D.; Frenklach, M.; Moriarty, N.; Eiteneer, B.; Goldenberg, W.; Bowman, C.; Hanson, R.; Song, S.; Gardiner, W.; Lissianski, V.; Qin, Z. GRI-Mech 3.0; http://www.me.berkeley.edu/ gri_mech/. (64) Verissimo, A.; Rocha, A.; Costa, M. Operational, combustion, and emission characteristics of a small-scale combustor. Energy Fuels 2011, 25 (6), 2469−2480. (65) Verissimo, A.; Oliveira, R.; Coelho, P.; Costa, M. Numerical simulation of a small-scale mild combustor. Proceedings of the 6th European Thermal Sciences Conference; Poitiers, France, Sept 4−7, 2012. (66) Pope, S. Computationally efficient implementation of combustion chemistry using in situ adaptive tabulation. Combust. Theory Modell. 1997, 1 (1), 41−63. (67) Cuoci, A.; Frassoldati, A.; Faravelli, T.; Ranzi, E. Formation of soot and nitrogen oxides in unsteady counterflow diffusion flames. Combust. Flame 2009, 156 (10), 2010−2022.

1122


Numerical Modeling of NOx Formation in Turbulent Flames Using a

Recommend Documents