Detailed Multi-dimensional Study of Pollutant Formation in a Methane

Feb 21, 2012 - *Telephone: +353-85-150-9336. ... Benefits of the method described here include its applicability to a wide range of combustion systems...
4 downloads 6 Views 5MB Size
Article pubs.acs.org/EF

Detailed Multi-dimensional Study of Pollutant Formation in a Methane Diffusion Flame Rory F. D. Monaghan,*,† Råbi Tahir,‡ Alberto Cuoci,§ Gilles Bourque,‡ Marc Füri,‡ Robert L. Gordon,‡ Tiziano Faravelli,§ Alessio Frassoldati,§ and Henry J. Curran† †

Combustion Chemistry Centre, National University of Ireland, Galway, Ireland Rolls-Royce Canada Limited, Montreal, Quebec H8T 1A2, Canada § Departmento di Chimica, Materiali e Ingegneria Chimica “Guilio Natta”, Politecnico di Milano, 20133 Milan, Italy ‡

ABSTRACT: This paper describes a method to produce chemical reactor networks (CRNs) consisting of large numbers of perfectly stirred reactors (PSRs) from computational fluid dynamics (CFD) simulations to predict pollutant emissions from combustion systems accurately, flexibly, and efficiently using detailed kinetic schemes and the kinetic post-processor (KPP) developed at Politecnico di Milano. Benefits of the method described here include its applicability to a wide range of combustion systems, its ability to predict emissions of a variety of pollutant species, and its speed. CFD and CFD−CRN simulation results of the Sandia D piloted methane−air diffusion round-jet flame are successfully validated against experimental data for axial velocity, mixture fraction, temperature, and speciation, including CO and NO mass fractions. A CRN consisting of a large number of PSRs is found to be required to simulate the system accurately, while ensuring independence of the solution from CRN size. The results of CFD−CRN analysis for a 1114 PSR network are used to study the pathways (thermal, prompt, N2O, and NO2) by which NO and NO2, the constituents of NOx, are formed in the flame. Results of CFD−CRN analysis show that NO is produced in the high-temperature (T > 1850 K) flame brush by a combination of the prompt, N2O, and thermal pathways and in the intermediate-temperature (1000 < T < 1600 K) post-flame region by a reversal of the NO2 pathway. NO is consumed in the fuelrich (mixture fraction, f > 0.43) region, where a low O atom concentration encourages a reversal of the prompt pathway (i.e., NO reburning), and in low-temperature (T < 1000 K) regions by the NO2 pathway, which oxidizes NO to NO2. Rate of production analysis, performed using CHEMKIN PRO at specified locations throughout the flame, shows that the trends of NO production and consumption observed in these simulations agree with expected and published results. Finally, the study predicts that, of the total NOx produced by the Sandia D flame, 47% is due to the prompt pathway, 32% is due to the N2O pathway, and 21% is due to the thermal pathway. As future steps in this work, the CFD−CRN method will be adapted and used to predict and study emissions from a range of more complex combustion systems.



with combustion systems can range from 10−7 m for reaction zone length scales to 100 m for combustor geometry, while timescales can range from 10−6 s for ignition and extinction phenomena to 100 s for pollutant formation.5 To describe pollutant formation accurately for all but the simplest fuels, 102 species and 103 reactions must also be considered. For these physical reasons, it is currently commercially impractical to include all phenomena of importance for pollutant formation in a single accurate computational fluid dynamics (CFD) simulation for a combustion system. This has led to the development of two distinct approaches to emission modeling in industry: (i) detailed treatment of fluid dynamics with reduced chemistry and (ii) reduced treatment of fluid dynamics with detailed chemistry. The former approach involves the use of CFD simulations that include simplified combustion models, such as infinitely fast chemistry or highly reduced kinetic mechanisms, sometimes consisting of as few as 1−3 reaction steps. These combustion models are linked to the turbulent flow field through turbulence− chemistry interaction models, such as eddy break-up (EBU)6 or eddy dissipation concept (EDC).7 Alternatively, probability

INTRODUCTION Gas turbines fueled by natural gas generate more than 20% of global electricity demand.1 This value is likely to increase because of (i) recent technological advances allowing for the cost-effective extraction of abundant alternative gas resources, (ii) decarbonization of coal-dependent energy systems for greenhouse gas emission reduction, (iii) uncertainty regarding the future role of nuclear energy, and (iv) the practical need to retain a fast-responding backup to increasing penetrations of intermittent renewables in the form of open-cycle gas turbines (OCGTs). Increasingly stringent legislative limits on carbon monoxide (CO) and oxides of nitrogen (NOx) emissions have driven development of low-emission gas turbine technologies, such as wet non-premixed systems and dry, lean premixed systems. European Union legislation, for example, limits NOx emissions to 37 ppmv at 3% O2 for high-efficiency combinedcycle gas turbine (CCGT) and combined heat and power (CHP) plants and to 25 ppmv for OCGTs.2 Refer to the comprehensive reviews by Correa3 and Correa et al.4 for overviews of gas turbine technologies with respect to NOx and CO emissions, respectively. Computer-based emission modeling techniques are used by gas turbine manufacturers and others in the field to evaluate new combustor designs. Bilger states that length scales associated © 2012 American Chemical Society

Received: November 27, 2011 Revised: February 7, 2012 Published: February 21, 2012 1598

dx.doi.org/10.1021/ef201853k | Energy Fuels 2012, 26, 1598−1611

Energy & Fuels

Article

highly studied flame, (ii) its simplicity compared to industrial combustion systems allows us to isolate and account for fundamental phenomena of importance, and (iii) for method development, it is important to have as accurate of a CFD solution as possible. The CFD−CRN method is then used to analyze the principal pathways of NOx formation in the flame and to visualize the locations in the flame where these pathways are most important. Finally, conclusions are drawn for the application of the CFD−CRN method to more complex combustion systems.

distribution function (PDF) methods can be used to close the chemical source term.8 Turbulent flow can be modeled using steady-state Reynolds averaged Navier−Stokes (RANS) methods, which include k−ε,9 k−ω,10 and Reynolds stress method (RSM)11 techniques. Alternatively, turbulent flow can be simulated at large scales and modeled at small scales using unsteady large eddy simulation (LES).12 The computational expense required to predict pollutant emissions to required levels of accuracy means that the approach of using detailed fluid dynamics and reduced chemistry is still very much under development for this application. The latter approach, reduced fluid dynamics with detailed chemistry, involves the use of chemical reactor networks (CRNs), which consist of combinations of idealized chemical reactors, such as perfectly stirred reactors (PSRs), plug-flow reactors (PFRs), and partially stirred reactors (PaSRs). When applied correctly, CRNs can adequately approximate the complex fluid flow fields present in combustion systems, allowing computational resources to be devoted to detailed chemistry and/or other physical and chemical phenomena of importance.13−22 Because of its empirical or semiempirical nature, a CRN is only reliably applicable to the exact geometry and set of operating conditions for which it was created. To increase the flexibility, accuracy, and computational efficiency of CRNs, recent work has focused on the coupling of CFD and CRN techniques. Such an approach, to which we refer as CFD− CRN, typically involves the following steps: (1) CFD simulation of a combustion system using a simplified combustion model that adequately captures the fluid flow field and temperature distribution but not necessarily the distributions of minor species, (2) application of a set of criteria to the CFD solution that divides the CFD domain into a CRN and calculates the mass flows between reactors, and (3) simulation of the combustion system using the newly created CRN with detailed chemical kinetic mechanisms. The advantage of the CFD−CRN approach as a design tool is that accurate fluid flow and chemical kinetic simulations can be combined in a cost- and time-efficient manner because (1) the number of equations solved using the CFD−CRN approach is small relative to that required for full CFD simulation using detailed chemistry and (2) the fact that CFD−CRN methods usually fix reactor temperatures when solving for detailed chemistry eliminates much of the numerical stiffness associated with chemical kinetic problems. CRN analysis can also use PSR energy equations to calculate reactor temperatures but at greater computational expense. Most of the CFD−CRN work presented in the literature use reactor-creation criteria based on the temperature and stoichiometry (i.e., local equivalence ratio or mixture fraction) and generally show good accuracy in predicting NOx concentrations at the combustor exit when the CFD solution accurately describes the temperature field.23−29 This paper presents a CFD−CRN method for the prediction of both CO and NOx for combustion systems. CFD simulation configuration, reactor-creation criteria, reactor properties calculation, and inter-reactor mass flow evaluation are discussed. The benefits of the CFD−CRN method include (i) the applicability of the method to a wide range of combustion systems, (ii) its ability to predict emissions of both CO and NOx species, as opposed to being tuned to predict one of them, and (iii) the speed with which accurate results can be obtained. Rigorous validation for the laboratory flame “Sandia flame D” is described. This involves a comparison of CFD−CRN-predicted distributions of CO and NO to experimental values at numerous locations within the flame, as well as at the domain exit. Sandia flame D is chosen for this work for the following reasons: (i) it is a stable, well-understood, and



CFD−CRN METHOD

Sandia Flame D. Sandia flame D or Sandia D is an atmospheric pressure piloted methane−air round-jet flame that is used to develop computational techniques for combustion.30,31 It has been the subject of numerous extensive RANS32 and LES33,34 simulation campaigns. Figure 1 shows the configuration of Sandia D along with the

Figure 1. Configuration of Sandia D showing the computational mesh for CFD simulation. computational mesh for CFD analysis, which is discussed in the next section. The size of the CFD domain shown is 0.72 × 0.72 m. The inlet velocities of the jet, pilot, and co-flow streams are ujet = 49.6 m s−1, upilot = 11.4 m s−1, and uco‑flow = 0.9 m s−1, respectively. Inlet temperatures are Tjet = 294 K, Tpilot = 1880 K (burned temperature), and Tco‑flow = 291 K. The inlet diameter of the jet, djet, is 7.2 mm, while the outer diameter of the annular pilot, dpilot, is 18.2 mm. The jet composition is 25% CH4 and 75% air by volume. The pilot burns a mixture of C2H2, H2, air, CO2, and N2 having the same enthalpy and equilibrium composition as CH4−air at an equivalence ratio (ϕ) of 0.77. The Reynolds number of the jet (Rejet) is roughly 22 400. The co-flow stream consists of dry air. Experimental mean and root-mean-square (rms) values of the temperature (T), mixture fraction (f), and mass fractions (Y) of CH4, O2, N2, H2O, CO2, H2, CO, OH, and NO at locations throughout the flame are available from ref 30. Experimental values of the mean axial velocity on the centerline and radially at an axial position of 45 jet diameters (z/djet = 45) are available from ref 35. The mixture fraction is defined in eq 1,31 where YC, YH, wC, and wH are elemental mass fractions and atomic weights of carbon and hydrogen, respectively. Subscripts 1 and 2 refer to values in the main jet and co-flow streams, respectively. Using this definition, the stoichiometric value of the mixture fraction for this system, fst, is 0.351.

f = (2(YC − YC,2)/wC + (YH − YH,2)/2w H) /(2(YC,1 − YC,2)/wC + (YH,1 − YH,2)/2w H)

(1)

CFD Simulation. A three-dimensional steady-state CFD simulation of Sandia D has been created in ANSYS FLUENT 12.1.36 The computational mesh consists of 193 968 cells, a cross-section of which is shown in Figure 1. The reason for using a three-dimensional mesh, as opposed to a two-dimension mesh, is to ensure that this approach is 1599

dx.doi.org/10.1021/ef201853k | Energy Fuels 2012, 26, 1598−1611

Energy & Fuels

Article

Table 1. Zone Limits and Reactor Criteria Used in CFD−CRN Analysis zone

zone limits

upstream zone (USZ) fuel-rich zone (FRZ) low-temperature flame zone (LFZ) high-temperature flame zone (HFZ)

z z z z

jet expansion zone (JEZ) co-flow zone (CFZ)

z > 0 m, 0.01 < f < 0.1 z > 0 m, 0 < f < 0.01

< > > >

0 0 0 0

reactor criteria

m m, 0.9 < f < 1.0 m, 0.1 < f < 0.9, T < 1400 K m, 0.1 < f < 0.9, T > 1400 K

n/a Δf = 0.01, Δz = 0.01 m ΔT = 100 K ΔT = 100 K (T < 1800 K), ΔT = 50 K (1800 K < T < 2000 K), ΔT = 2 K (T > 2000 K) ΔT = 100 K Δz = 0.2 m

Figure 2. (a) Mass-averaged zone temperatures and (b) mass-averaged reactor temperatures in kelvin. developed in as general of a way as possible and can be adapted for more complex combustion systems that may not exhibit axisymmetric behavior. The mesh has been highly refined in regions of importance for this flame. Turbulence boundary conditions are defined as follows for the inlet streams: turbulent intensity of 5% for all inlets, with hydraulic diameters of 7.2 and 10 mm for the jet and pilot inlets, respectively, and turbulent viscosity ratio of 10% for the co-flow inlet. The RSM with the quadratic pressure-strain model37 is used to model viscosity. This level of detail was found to be required to correctly predict the axial extent of the flame. Thermal radiation is not modeled because its inclusion was found by this work to have a negligible effect on the predicted temperature field. Combustion is modeled using the DRM22 skeletal mechanism, which consists of 22 reacting species, 2 nonreacting species, and 104 reactions.38 EDC is used to model the turbulence−chemistry interaction. Note that DRM22 does not include a NOx submodel; therefore, the FLUENT NOx post-processor is used to produce CFD predictions of NOx for comparison to CFD−CRN predictions. The NOx post-processor has the following configuration: thermal, prompt, and N2O pathways activated, instantaneous thermal [O] and [OH] models, quasi-steady N2O model, temperature-only β-PDF turbulence interaction with 10 PDF points, and algebraic temperature variance with a global maximum temperature. SIMPLE pressure−velocity coupling is used, and all transport equations are solved using second-order discretization. Validation and results of CFD simulation are presented in later sections. CRN Creation. Reactor Criteria and Properties. The CFD domain is discretized into a network of PSRs using a FLUENT user-defined function (UDF) coded in C.39 The domain is first divided into zones that generally define important regions in nonpremixed flames. The limits of each zone are defined primarily by the mixture fraction, temperature, and axial position and are shown in Table 1. Zones are defined as follows: the upstream zone (USZ) is the region between the furthest upstream part of the domain and the outlets of the jet and pilot at z = 0 m (i.e., before they mix and react); the fuel-rich zone (FRZ) is where initial mixing and fuel decomposition occurs; the low- and hightemperature flame zones (LFZ and HFZ) are where heat release rates are highest (i.e., the flame brush); the jet expansion zone

(JEZ) is the region of combustion products at elevated temperatures; and the co-flow zone (CFZ) is the region of lowtemperature, slow-moving air, whose purpose is to prevent recirculation in the domain. The locations of the zones within the CFD domain are shown by their mass-averaged temperatures in Figure 2a. The extents of the domain shown are from −0.036 to 0.72 m (from −5djet to 100djet) axially and 0.36 m (50djet) radially. The vertical lines indicate divisions of 0.072 m (10djet). Within each zone, reactor criteria, which are also shown in Table 1, are then applied on the basis of the general physical and chemical phenomena that are expected to be found in them. Reactor criteria are expressed as ΔZ, where Z is a cell property. Consider neighboring computational cells, whose z coordinates, temperatures, and mixture fractions place them within the limits of the same zone. If the differences in properties between the cells are less than those defined for their assigned zone (Δf, ΔT, and Δz), they are considered to be the same reactor in the CRN. Massaveraged temperatures of the PSRs created by this method are shown as a contour plot in Figure 2b. Applying these reactor criteria produces a CRN consisting of 1114 PSRs. Note that validation of the CFD−CRN method is not presented here but in later sections. The necessity of such a large number of PSRs is also discussed in later sections. Mass-averaged properties, which include the temperature, mixture fraction, and axial velocity for each PSR are calculated using eq 2, in which ρi, Vi, and Zi are the density, volume, and property of interest of cell i within a PSR, which consists of Ncell,PSR cells. Temporal temperature variance for each cell (σT2 ) is calculated, for reasons explained later, using the approximate (i.e., algebraic) version of the temperature variance transport equation (eq 13.1-113 in ref 36) shown in eq 3, where μt is the turbulent viscosity, given by eq 4 for RANSbased viscosity models, k is the turbulent kinetic energy, ε is the turbulent dissipation rate, Cg, Cd, and Cμ are constants of values 2.86, 2.00, and 0.09, respectively, and ∇T is the magnitude of the temperature gradient. The EDC model, which is used in the CFD simulation, assumes that the chemical reaction occurs only in small turbulent structures, which have volume fraction ξ* as defined by eq 5, where ν is the kinematic viscosity and Cξ is a constant of value 2.1377. 1600

dx.doi.org/10.1021/ef201853k | Energy Fuels 2012, 26, 1598−1611

Energy & Fuels

Article

Mass-averaged values for the temperature variance and EDC volume fraction are calculated for each PSR from individual cell values using eq 2. The volume of each PSR available for reaction (VPSR) is the product of the total volume of cells in the PSR and the mass-averaged EDC volume fraction for the PSR. Ncell,PSR



= ZPSR ̅

Ncell,PSR

ρiVZ i i/

i=1 2 σT =

μ t k Cg (∇T )2 ρ ε Cd

μt = Cμρ

the detailed kinetic scheme for the CRN. During this phase, each reactor is solved using a local Newton method with the possible use of a time-stepping approach to improve robustness. When the residuals of all of the equations are sufficiently low, a modified global Newton method is applied to the whole system. If the global Newton method does not converge, a “false transient” method is applied to ensure a better approach to the solution of the whole system. After the time derivatives are added to the steady-state equations of each reactor, the full system becomes a system of ordinary differential equations (ODEs) rather than a nonlinear system. In this case, the solution is obtained via the backward-Euler method. The solution methods employed by KPP are described in detail in refs 44 and 45. During evaluation for this work, KPP proved to be robust and computationally efficient for a range of CRN sizes and three different kinetic mechanisms when compared to 64-bit CHEMKIN PRO.46 Table 2 shows the relative performances during software evaluation of CHEMKIN PRO and KPP in solving CRNs of 108−1114 PSRs using

k2 ε

⎛ νε ⎞1/4 ξ* = Cξ⎜ 2 ⎟ ⎝k ⎠

∑ i=1

ρiVi

(2)

(3)

(4)

(5)

Table 2. Relative Performances of CHEMKIN PRO and KPP in Solving CRNs of Different Sizes and Mechanisms

Inter-reactor Mass Flows. Two types of mass exchange are allowed between reactors: advective and turbulent diffusive. Advective mass flow between two PSRs is calculated in the UDF by summing all of the cell−cell mass flow rates at the boundary between the two PSRs. Turbulent diffusive mass flow between two PSRs is calculated through the indirect use of a Peclet number for turbulent mass transfer (Pem,t), defined in eq 6, where u is the velocity, lcell is the cell length scale, defined in eq 7, and Dm,t is the turbulent mass diffusivity, defined in eq 8. The turbulent Schmidt number (Sct) is an input parameter for the CFD simulation and is set to its default value of 0.7. For each cell−cell advective mass flow (ṁ adv), eq 9 can be used to calculate an equivalent turbulent diffusive mass flow (ṁ diff,t). Because diffusion does not involve a bulk movement of mass but rather an exchange of species, there are equal and opposite turbulent diffusive mass flows between cells at PSR boundaries. In a similar manner to advective mass flow, total turbulent diffusive mass flows between two PSRs are calculated by summing all of the cell−cell flow rates at the boundary.

Pe m,t =

ulcell Dm,t

1/3 lcell = Vcell

Dm,t = ṁ diff,t =

time to solution (h:min:s)

(6) (7)

ρADm,t ṁ adv = Pe m,t lcell

mechanism

CHEMKIN PRO

KPP

108 PSRs 108 PSRs 108 PSRs 211 PSRs 211 PSRs 255 PSRs 255 PSRs 372 PSRs 411 PSRs 523 PSRs 617 PSRs 722 PSRs 809 PSRs 1114 PSRs

GRI 1.2 GRI 3.0 NUIG C2 GRI 1.2 NUIG C2 GRI 1.2 NUIG C2 NUIG C2 NUIG C2 NUIG C2 NUIG C2 NUIG C2 NUIG C2 NUIG C2

0:14:20 0:25:11 failed 0:23:45 failed 1:57:48 failed failed failed failed failed failed failed failed

0:00:06 0:00:30 0:00:39 0:00:16 0:00:26 0:00:16 0:01:28 0:01:48 0:02:04 0:02:26 0:02:38 0:03:20 0:03:50 0:09:30

improvement factor 143 50 89 442

three different kinetic mechanisms: GRI 1.2 (32 species and 175 reactions),47 GRI 3.0 (53 species and 325 reactions),48 and NUIG C2 (103 species and 582 reactions). All evaluations were performed on the same 2.80 GHz Intel central processing unit (CPU) with a 64-bit operating system. As a result of the evaluation, KPP was selected to solve the 1114 PSR CRN using the NUIG C2 mechanism. Upon solution of the CRN in KPP, the graphical post-processor in ANSYS FLUENT is used to visualize results. The necessity of such large CRNs is discussed later. KPP uses mass-averaged reactor temperatures (T̅ ) obtained from eq 2. In principle, for a given reactor, the rates of all of the reactions involved in the detailed kinetic scheme should be evaluated at this temperature. In reality, however, turbulent combustion conditions result in temporal fluctuations of the temperature and species concentration that will strongly influence the characteristics of the flame and especially the formation of NOx. The relationships between the NOx formation rate, temperature, and species concentration are highly nonlinear. Hence, if time-averaged composition and temperature are employed to predict the mean NOx formation rate, significant errors will result. As a first approximation, fluctuations of species concentration can be neglected. On the contrary, fluctuations of the temperature can be expected to have a large effect, because of the highly nonlinear dependence of reaction rates upon the temperature in the Arrhenius law (eq 10). As a consequence, KPP accounts for the temperature fluctuations by considering a PDF, which describes the corresponding time variation.36 The mean turbulent reaction rate (r̃) is then calculated using eq 11, where r is the instantaneous reaction rate and p(T) is the temperature PDF. The limits of integration are determined from the minimum and maximum values of

μt ρSc t

CRN size

(8)

(9)

CRN Solution. Detailed Kinetic Mechanism. The detailed kinetic mechanism used for this work is the C0−C2 portion of the C5 mechanism developed in the Combustion Chemistry Centre (C3) at the National University of Ireland, Galway (NUIG).40−42 The NOx submechanism developed at the Centre National de la Recherche Scientifique (CNRS) in Orléans, France, has been added for NOx emission prediction.43 The C0−C2 portion of the NUIG mechanism combined with the NOx submechanism contains 103 species and 582 reactions in total and is henceforth referred to in this paper as “NUIG C2”. Software. The software used to solve the CRN with detailed chemistry is the kinetic post-processor (KPP), developed at Politecnico di Milano (Polimi). KPP solves the CRN, which corresponds to a large system of nonlinear equations, using a hybrid approach, combining the application of (i) successive substitutions, (ii) a false transient method, and (iii) a global Newton method. This procedure is briefly summarized as follows. A first guess solution is found by iteratively solving the sequence of individual reactors with successive substitutions. This is needed to have physically meaningful first estimations of chemical species, which are not considered in the CFD simulation but which are available in 1601

dx.doi.org/10.1021/ef201853k | Energy Fuels 2012, 26, 1598−1611

Energy & Fuels

Article

Figure 3. Experimental (□) and CFD−CRN-predicted profiles of (a) CO and (b) NO mass fractions in the radial direction near the exit of the CFD domain at z/djet = 75 for CRNs consisting of 57 PSRs (black), 108 PSRs (red), 211 PSRs (blue), and 532 and 1114 PSRs (both green).

Figure 4. Experimental (□) and CFD (- - -) values of the axial velocity on the (a) centerline and (b) z/djet = 45 (red).

Figure 5. Experimental (□), CFD (- - -), and CFD−CRN () values of the temperature on the (a) centerline and (b) z/djet = 30 (black), 45 (red), and 75 (blue). the temperature in the CFD solution. The PDF of the temperature can be assumed to be a β-PDF,36 a clipped-Gaussian PDF,49 or a double delta PDF.50 In most cases, the three approaches give very similar results. To build the PDF of the temperature, T̅ and σT2 of each reactor are required and are calculated by the UDF and supplied to KPP as inputs.

k = AT β exp(− Ea /RT ) r̃ =

primarily due to the fact that larger CRNs have been achieved by refining reactor criteria mostly at high temperatures, which are generally important for NO production. Predicted exit mass fractions of CO and NO given by 532 PSR and 1114 PSR networks are identical and are thus shown using the same green lines in Figure 3. The solution produced by the 532 PSR and 1114 PSR networks are therefore determined to be independent of the CRN size. The larger of these two networks, 1114 PSRs, is used to obtain the results that follow.

(10)

T

∫T max r(T , P , Yk)p(T )dT min



(11)

CFD−CRN METHOD VALIDATION Results of CFD and CFD−CRN simulation of Sandia D are validated with experimental data described in previous sections.30,35 Predictions of the temperature, mixture fraction, and mass fractions of O2, CO2, CO, and NO are compared to experimental values axially along the flame centerline and radially at z/djet = 30, 45, and 75 in Figures 5−10. CFD predictions and experimental data for axial velocity on the centerline and at z/djet = 45 are shown in Figure 4. For all profiles, open symbols, dashed lines, and solid lines refer to experimental results, CFD results, and CFD−CRN results,

CRN Size Independence. Experimental and CFD−CRN-predicted profiles of CO and NO mass fractions in the radial direction near the exit of the CFD domain are shown in Figure 3. CFD−CRN results are shown for CRNs consisting of 57 PSRs (black), 108 PSRs (red), 211 PSRs (blue), and 532 and 1114 PSRs (both green). Note that full method validation is described in the next section. This section focuses solely on ensuring CRN size independence. The 57 PSR network greatly overestimates the CO concentration and underestimates the NO concentration. Increasing the number of PSRs to 108 shows considerable improvement for both species. Further increases in the size of the CRN to 211 and 532 PSRs have a limited effect on CO accuracy but result in steady improvement of NO prediction. This is 1602

dx.doi.org/10.1021/ef201853k | Energy Fuels 2012, 26, 1598−1611

Energy & Fuels

Article

Figure 6. Experimental (□), CFD (- - -), and CFD−CRN () values of the mixture fraction on the (a) centerline and (b) z/djet = 30 (black), 45 (red), and 75 (blue).

Figure 7. Experimental (□), CFD (- - -), and CFD−CRN () values of the O2 mass fraction on the (a) centerline and (b) z/djet = 30 (black), 45 (red), and 75 (blue).

Figure 8. Experimental (□), CFD (- - -), and CFD−CRN () values of the CO2 mass fraction on the (a) centerline and (b) z/djet = 30 (black), 45 (red), and 75 (blue).

Figure 9. Experimental (□), CFD (- - -), and CFD−CRN () values of the CO mass fraction on the (a) centerline and (b) z/djet = 30 (black), 45 (red), and 75 (blue).

respectively. For radial profiles, black, red, and blue refer to z/djet = 30, 45, and 75, respectively. The locations of PSRs within the CRN are seen in the discrete nature of the solid lines, particularly in low-temperature regions, where reactor criteria

are less exacting, i.e., where reactor criteria are coarser. Note the different x axis scale in Figure 9b. CFD predictions of the axial velocity, mixture fraction, and major species (O2 and CO2) show satisfactory agreement with 1603

dx.doi.org/10.1021/ef201853k | Energy Fuels 2012, 26, 1598−1611

Energy & Fuels

Article

Figure 10. Experimental (□), CFD NOx post-processor (- - -), and CFD−CRN () values of the NO mass fraction on the (a) centerline and (b) z/djet = 30 (black), 45 (red), and 75 (blue).

Figure 11. Plots of experimental (● ), CFD ( □), and CFD−CRN (■) mass fractions of (a) OH and (b) NO versus mixture fraction at z/d jet = 45.

experimental results. Radial temperature profiles are slightly overpredicted by the CFD simulations (Figure 5b), but this does not appear to have a large impact on major species predictions because of the fact that Sandia D is a diffusion flame, in which the overall rate of fuel consumption is at least partly controlled by turbulent diffusion. For major species and CO, CFD−CRN results are of generally similar accuracy to those of CFD. As previously stated, to provide a comparison to CFD−CRN NO predictions, the FLUENT NOx postprocessor is used. As illustrated by the dashed lines in Figure 10, CFD NOx post-processors frequently overestimate NO mass fractions for two reasons: (1) global chemistry is used for reaction pathways, which, for the prompt pathway, in particular, lead to inaccuracies, and (2) mass fractions of radicals are frozen, meaning NOx-forming reactions, which, in reality, would deplete certain species, are modeled to proceed at too high of a rate. CFD−CRN-predicted NO profiles compare favorably to previous CFD-based NO predictions for Sandia D.32,34 Raw experimental mass fractions (scattered points) of OH and NO as well as CFD-predicted (open symbols) and CFD− CRN-predicted (filled symbols) values at axial position z/djet = 45 are plotted against the mixture fraction, as defined in eq 1, in Figure 11. This shows that the values predicted by the CFD− CRN method (filled symbols) are more accurate than those predicted by CFD only (open symbols) and generally fall within experimental scatter. The data presented in Figure 11 agree well with published results for H2−CO diffusion flames.51

The CFD−CRN method shows overall satisfactory accuracy in predicting CO and NO mass fractions throughout the flame and is considered to be adequately validated for the purposes of this work.



RESULTS Overall Rates of NOx Production. Contour plots of NO, NO2, NOx (NO + NO2), and N2O mass fractions predicted using the CFD−CRN method are shown in Figure 12. All plots in the figure use the same scale to allow for comparison. Peak NO and NOx mass fractions occur in the region 45 < z/djet < 55, which is slightly on the lean side of the region of the peak temperature seen in Figure 2b. NO2 is formed in lower temperature regions away from the flame. This agrees with the current understanding and the findings of previous work for hydrocarbon diffusion flames.51,52 Contour plots of calculated NO, NO2, NOx, and N2O rates of production in units of kg m−3 s−1 are shown in Figure 13. In this figure, red represents high rates of species production, blue represents high rates of species consumption, and green represents no production or consumption. Again, all plots in the figure use the same scale for comparison. Figure 13a shows high rates of NO production in the highest temperature regions of the flame. Immediately inside the high NO production zone is a band of NO consumption (“rich NO consumption zone”), where radially inward diffusing NO is consumed. There is also a region of NO consumption and consequent NO2 production that roughly corresponds to a mixture fraction contour of 0.1 1604

dx.doi.org/10.1021/ef201853k | Energy Fuels 2012, 26, 1598−1611

Energy & Fuels

Article

Figure 12. Contour plots of CFD−CRN predictions of (a) NO, (b) NO2, (c) NOx (NO + NO2), and (d) N2O mass fractions.

procedure involves indentifying the key reactions for each of the formation pathways and isolating, in turn, each of the pathways to determine their contribution to the NO rate of production at each point in the flame. Key reactions for each pathway are identified as follows52 and shown in Table 3: (i) thermal pathway, reactions 1 and 2, which convert N2 and O2 in the oxidant stream to NO; (ii) prompt pathway, any reactions involving HCN, which is a product of the rate-limiting step of this pathway, reaction 3; (iii) N2O pathway, any reactions involving N2O, including reactions 4−8, which Kuo52 identifies as key; and (iv) NO2 pathway, any reactions involving NO2, of which reactions 9−11 are important. The high activation energy required for reaction 1 to break the triple bond of N2 (Ea = 316 kJ mol−1, 75.5 kcal mol−1) means that the thermal pathway is important at high temperatures, hence, its name. The prompt pathway relies on the presence of hydrocarbon radicals to react with N2 and, thus, is generally assumed to be important for rich systems or moderately lean systems where highly non-equilibrium chemistry takes place, e.g., in high-intensity flame zones. Note that the N2O pathway as defined in this study covers two subpathways as defined by Kuo:52 (i) N2, O, and M reacting to produce N2O and (ii) NO reacting with N-containing radicals to form N2O. In both cases, N2O is subsequently attacked by O and H atoms. The first subpathway, especially reaction 6, is thought to resemble a low-temperature equivalent of the

(“lean NO consumption zone”). NO2 production is mostly confined to this region, as seen in Figure 13b. When the production rates of NO and NO2 are combined into a single NOx production rate, Figure 13c shows two distinct regions of NOx formation: (i) production of NO in the highest temperature region and (ii) production of NO2 in cooler regions. The next section analyses the chemical kinetic pathways by which NO and NO2 are produced in these zones. NOx Production Pathways. Oxides of nitrogen are generally thought to form via five production pathways:3,52−54 (i) thermal (or Zel’dovich), (ii) prompt (or Fenimore), (iii) N2O (nitrous oxide) intermediate, (iv) oxidation of NO to NO2, and (v) release of fuel-bound nitrogen. An additional route, via the formation of NNH, has been postulated by Bozzelli and Dean55 and is currently included in some although not all detailed kinetic schemes. The NNH pathway is not represented in the CNRS NOx submechanism43 used in this analysis and, therefore, is not considered further in this work. The fuel in the case of Sandia D is pure CH4; thus, the fuelbound pathway is also omitted from further consideration. This leaves thermal, prompt, N2O, and NO2 pathways. The overall interaction of the pathways is illustrated in Figure 14, with thermal in red, prompt in blue, N2O in magenta, and NO2 in green. To quantify the relative importance of the NOx formation and consumption pathways, the procedure described by Fackler et al.56 in their NOx pathway analysis of landfill gas combustion in a jet-stirred reactor (JSR) is employed. This 1605

dx.doi.org/10.1021/ef201853k | Energy Fuels 2012, 26, 1598−1611

Energy & Fuels

Article

Figure 13. Contour plots of CFD−CRN predictions of (a) NO, (b) NO2, (c) NOx (NO + NO2), and (d) N2O rates of production in kg m−3 s−1.

flame regions, where HO2 oxidizes NO to NO2 via reaction 9. It is reversible in regions near the flame via reactions 10 and 11. Figure 15 shows the results of NOx pathway analysis using the CFD−CRN method. Predicted contributions of thermal (red), prompt (blue), N2O (magenta), and NO2 (green) pathways to the NO rate of production (in kg m−3 s−1) are found by isolating each pathway by deactivating the key reactions mentioned above. Also shown in Figure 15 are experimental and CFD−CRN-predicted NO mass fractions on the centerline and at z/djet = 30, 45, and 75, which are identical to the data displayed in Figure 10. NO mass fraction scales are shown on the right of each plot. The following regimes, which are labeled in Figure 15a, are apparent: (i) NO production in the high-temperature flame brush via a combination of the prompt (blue), thermal (red), and N2O (magenta) pathways, (ii) NO consumption on the rich side of the flame brush (“rich NO consumption zone” from Figure 13) through a reversal of the prompt pathway (or reburning), and (iii) NO production on both the rich and lean sides of the flame brush, followed by (iv) NO consumption via the NO2 pathway. Although these regimes are highlighted in Figure 15a only, they are also visible in the radial profiles. The regimes of NOx production and consumption are discussed in the next section.

Figure 14. Schematic of reaction channels for NOx pathways: thermal (red), prompt (blue), N2O (magenta), and NO2 (green).

thermal pathway, because they involve the reaction of N2 from the oxidant stream. The presence of a third body in reaction 6 means that the N2O pathway is thought to be particularly important in lean (low radical concentration), premixed (low temperature), and high-pressure (high third body concentration) combustion systems. The NO2 pathway is generally thought to be of importance in lower temperature



DISCUSSION To determine the most important specific reactions within each NOx formation or consumption pathway and to ensure agreement 1606

dx.doi.org/10.1021/ef201853k | Energy Fuels 2012, 26, 1598−1611

Energy & Fuels

Article

pathway. The other species shown in Figure 16, NH and HNO, are intermediates in the prompt pathway.

Table 3. Key Reactions for NOx Production Pathways pathway

reaction

k12f

NCO + O XooooY NO + CO

k1f

O + N2 XoooY NO + N thermal

k1r

(reaction 1)

k12r

(reaction 2)

HCN + O XooooY NCO + H

k2f

N + O2 XoooY NO + O k2r

prompt

k13f

k13r

k3f

CH + N2 XoooY HCN + N k3r

k4r

(reaction 4)

k5f

NCO + NO XoooY N2O + CO k5r

N2O

(reaction 5)

k 6f

N2 + O + M XoooY N2O + M k 6r

(reaction 6)

k 7f

N2O + H XoooY N2 + OH k 7r

(reaction 7)

k8f

N2O + O XoooY NO + NO k8r

(reaction 8)

k 9f

NO + HO2 XoooY NO2 + OH k 9r

NO2

(reaction 9)

k10f

NO2 + H XooooY NO + OH k10r

(reaction 10)

k11f

NO2 + O XooooY NO + O2 k11r

(reaction 13)

Regime II: NO Consumption via Reburning in the Fuel-Rich Region. The CFD−CRN-predicted axial profile of the NO mass fraction shows a dip at around z/djet = 35, which in Figure 15a is seen to correspond with the reversal of the prompt pathway, labeled as regime ii. This is referred to as reburning (CxHy + NO = HCN + products). Reburning is also observed near the flame centerline in Figure 15b. The predicted decrease of NO in this region is explained primarily by the relatively low concentration of O atoms under locally rich conditions (ϕ > 1.2 for a premixed system or f > 0.43 for this system). The main effects of this are (1) to slow the conversion of HCN formed by the primary prompt reaction (reaction 3) to NO, (2) to reduce the rate of conversion of CH3 to CH2O via O addition, which, in turn, leads to increased concentrations of C, CH, CH2, and CH3, leading to increased rates of NO−HCN recycling via reactions 14−17, and (3) to reverse reaction 1, which consumes NO.52 ROP analysis for regime ii shows reactions 14−17 are responsible for the vast majority of NO consumption, as shown in Figure 17, which indicates agreement with the expected behavior.

(reaction 3)

k4f

NH + NO XoooY N2O + H

(reaction 12)

k14f

(reaction 11)

CH3 + NO XooooY HCN + H2O k14r

(reaction 14)

k15f

with observed and expected trends, a further rate of production (ROP) analysis is performed for PSRs, in which the regimes identified in the previous section are most apparent. The temperatures, residence times, and inlet compositions shown in Table 4 are exported from the CFD−CRN solution to CHEMKIN PRO for ROP analysis for the four regimes below. Recall the regimes identified in the table are the (i) dominance of the prompt pathway in the flame brush, (ii) reburning in the fuel-rich region, (iii) intermediate-temperature NO production via NO2 pathway, and (iv) low-temperature NO consumption via NO2 pathway. Regime I: NO Production via Prompt, Thermal, and N2O Pathways in the High-Temperature Flame Brush. The highest rate of NOx production is predicted by the CFD− CRN method to occur in the turbulent flame brush (see Figure 13) via a combination of the prompt, thermal, and N2O pathways (see regime i in Figure 15a). ROP analysis reveals the complex interactions between these pathways in the hightemperature regime, which are illustrated in Figure 16. The thickness of each line denotes its contribution to overall the NO rate of production. The role of the thermal pathway is seen in the N2 + O channel, i.e., reaction 1, as is the role of the N2O pathway in the N2−N2O−NO channels, i.e., reactions 6−8. The rest of the channels seen in Figure 16 represent the prompt pathway. The largest single contributor to NO production is seen to be the NCO + O channel, i.e., reaction 12, with NCO being supplied primarily by the oxidation of HCN (reaction 13), which agrees with the expected behavior.52 HCN is a product of reaction 3, the key reaction in the prompt

CH2 + NO XooooY HCN + OH k15r

(reaction 15)

k16f

CH2 + NO XooooY HCNO + H k16r

(reaction 16)

k17f

CH2(s) + NO XooooY HCN + OH k17r

(reaction 17)

Regime III: NO Production via NO2 Pathway in Intermediate-Temperature Regions. The region near the end of the centerline in Figure 15a shows dominance of the NO2 pathway and is labeled as regime iii. ROP analysis shows that the most important reactions for NO production are reactions 11 and 10 (shown above), in which NO2 is reduced to NO by O and H atoms, respectively. The production of NO is slowed somewhat by reaction 18. This again indicates agreement with the observed behavior. ROP analysis also shows that reaction 9 (NO + HO2 = NO2 + OH) proceeds in both directions at various stages in the regime. k18f

NO + OH + M XooooY HONO + M k18r

(reaction 18)

Regime IV: NO Consumption via NO2 Pathway in LowTemperature Regions. Low-temperature (1850 K) flame brush by a combination of the prompt (CHx + N2), N2O, and thermal (N2 + O) pathways and (2) in the intermediatetemperature (1000−1600 K) post-flame region by the reversed NO2 pathway, in which NO2 is reduced to NO by O and H atoms (reactions 10 and 11). NO is consumed in two regions: (1) in the fuel-rich ( f > 0.43) region, where a low O atom concentration encourages recycling of NO to HCN via CHx + NO reburning reactions (reactions 14−17), and (2) in lowtemperature (