Arab J Sci Eng (2014) 39:8287–8305 DOI 10.1007/s13369-014-1398-7

RESEARCH ARTICLE - MECHANICAL ENGINEERING

Computational Analysis of the Extended Zeldovich Mechanism Lucky Anetor · Christopher Odetunde · Edward E. Osakue

Received: 21 February 2014 / Accepted: 30 June 2014 / Published online: 11 October 2014 © King Fahd University of Petroleum and Minerals 2014

Abstract A computational study was used to investigate the effects of temperature and pressure on the extended Zeldovich mechanism. Two temperatures, namely 2,600 and 1,900 K, were used while keeping the species concentration of all component products constant at 1.0 mole/m3 . These temperatures were selected because they represent typical operating conditions for internal combustion engines. Pressure was varied by decreasing the species concentration for all component species by 10 % while keeping the temperatures constant. The pressures were 1.513 × 105 and 1.362 × 105 N/m2 for T = 2,600 K and 1.106 × 105 and 9.95×104 N/m2 for T = 1,900 K. The estimate of the uncertainties in the global errors of [NO], [N], and [O] for the conditions investigated was found to vary between (±4.10 and ±4.34 %), (±9.33 and ±10.1 %), and (±5.73 and ±9.07 %), respectively. The numerical results show that high temperatures result in faster rate of production of both [NO] and [O] and faster rate of [N] depletion. At moderately low temperatures (≈1,900 K), the rate of depletion of both oxygen and nitrogen atoms is very high. It was observed that changes in pressure had a minor or at best marginal influence on the production of [NO], [N], and [O]. The work also shows that for all the equilibrium points investigated, only one in each case

was physical and the others were non-physically realizable because at least one of their equilibrium points was negative. Keywords Zeldovich mechanism · Chemical kinetics · Equilibrium point · ICE · NO X

L. Anetor (B) Department of Mechanical Engineering, Nigerian Defence Academy, Kaduna, Nigeria e-mail: [email protected] C. Odetunde Department of Aeronautics and Astronautics Engineering, Kwara State University, Ilorin, Nigeria

1 Introduction

E. E. Osakue Department of Industrial Technology, Texas Southern University, Houston, TX, USA

Fossil fuel reserves are presently been depleted at an unsustainable rate. It is estimated that in the United States of Amer-

123

8288

Arab J Sci Eng (2014) 39:8287–8305

ica alone, the automobile accounts for 50 % of the petroleum consumption, 20 % of the total energy consumed, and 50 % of the pollutants entering the atmosphere. Since combustion processes in the automotive, utility, and industrial applications are the main sources of oxides of nitrogen (NOx) production, it does not take much to realize that a relatively minor improvement in efficiency and emissions of these devices would yield significant economic, environmental, and health benefits. Modeling and simulation can play a crucial role in engine design, analyses, and optimization studies. Multidimensional, multizone simulations can highlight the effects of essential parameters such as fuel composition, engine geometry, and operating conditions on fluid-dynamic variables such as temperature, pressure, velocities, and species concentrations in a variety of engines. Such complex models can also use detailed chemical kinetics to describe fuel combustion and formation of soot, unburned hydrocarbon (UHC), and NOx. Conducting a transient multidimensional numerical simulation of the entire engine cycle consisting of intake, compression, power, and exhaust strokes is a formidable computational task [1,2]. Such simulations can take several hours to days, depending on the complexity of the flow field and/or chemistry model, the size of the computational domain, and the central processing unit (cpu) cycles available. Because of the immense computational efforts required, these methods are usually not used during the exploratory stages of engine design, development, and assessment of (NOx) reduction techniques. On the other hand, comprehensive, computationally efficient, physics-based, quasi-dimensional modeling tools require moderate computational resources, and the computational time is usually measured in seconds. Development of such modeling tools usually speed-up design, analyses, and optimization of a variety of NOx reduction techniques. These design tools can also be used for comparing the efficacy of various NOx reduction methods. It is pertinent to mention that experimental studies are useful in validating and calibrating these numerical models [3–5]. These simulation models, in turn, can be used to examine a range of operating engine conditions (cylinder temperature, pressure, revolutions per minute (rpm), compression ratio, equivalence ratio, etc.) and other emission-reducing technologies (catalytic converters, exhaust gas recirculation (EGR), etc.), thus complementing experimental efforts. The time evolution of NOx during an engine cycle can be computed by using finite-rate or equilibrium chemistry computations. These have been used in various quasi-dimensional models [6,7]. Both techniques have their advantages and disadvantages. Equilibrium chemistry computations depend only on the temperature, pressure, and initial elemental mixture composition (initial number of atoms of C, H, O, and N in most combustion systems). A set of about 20 species

123

is usually adequate to describe the combustion products in most fuel–air combustion systems. Furthermore, equilibrium NOx concentrations are not very sensitive to the initial concentration of trace species such as atomic nitrogen and oxygen. However, the equilibrium approach can be inaccurate in predicting the time history of NOx under normal engine operating conditions. This is because equilibrium computations assume that the elementary processes in NOx formation are infinitely fast compared to the residence time of the reacting mixture for a given set of conditions (temperature, pressure, and composition). During the compression and expansion cycle, the engine cylinder volume (and hence species concentration), pressure, and temperature are rapidly varying and are comparable to the time constants of the reactions producing NOx. This is especially true during the early phase of combustion when the cylinder temperature is not high or during the expansion stroke in the late phases of combustion. Therefore, the use of the equilibrium approach can lead to inaccurate prediction of the temporal evolution of thermal NOx concentration [8]. Finite-rate chemistry computations do not make any assumptions about the chemical kinetic timescales and therefore do not have the shortcomings of the chemical equilibrium models. However, finite-rate chemistry calculations do require the use of appropriate chemical kinetic mechanisms (for example, the extended Zeldovich mechanism) to describe the formation of NOx with reasonable accuracy. For stability and accuracy of time-marching process, it is necessary to ascertain that the time step, t, and initial conditions of trace species such as atomic nitrogen and oxygen are carefully selected. Finite-rate chemistry computations which uses an implicit variable time-stepping scheme and a reduced chemistry mechanism for NOx can be fast and robust [6]. While the knowledge of time history of NOx computed from finite-rate kinetics can provide useful information on its formation process from a kinetics standpoint, for most engine designs, the overall engine-out NOx is of greater interest.

2 Literature Review As stated above, the amount of NOx formed depends on several factors such as combustion temperature, equivalence ratio, turbulence and mixing of the fuel–air mixture, and the residence time of the combusting gases in the combustion chamber. A very good understanding of the role of these parameters in the overall process of NOx formation mechanism plays a crucial role in developing its prevention methods. Such methods can generally be grouped into two categories, namely combustion modification and aftercombustion processes. The principle of combustion modification involves changing the NOx formation mechanism either by varying the reaction rates of elementary processes

Arab J Sci Eng (2014) 39:8287–8305

through temperature or by altering the concentrations of nitrogen and oxygen in the reacting zone. The application of exhaust gas recirculation (EGR) method is a popular combustion modification technique, which is normally used to reduce NOx in automotive engines [9]. The use of exhaust gas recirculation lowers the peak combustion temperature and also reduces the overall oxygen content in the initial charge, thus leading to reduced NOx production levels. However, the use of EGR has some adverse effects on the overall performance of an engine, such as reduction in engine power output, increased soot, increased cycle-to-cycle variation, and probably unstable combustion at low loads. Nitrogen-enriched air is another technique that has shown promising results in NOx reduction in stationary natural gas engines [10]. A comparative study of the traditional method of NO computations (simplified extended Zeldovich finite-rate chemistry) in quasi-dimensional engine codes with full finite-rate chemistry was conducted by [11]. In this work, single-zone and two-zone models were used to compute the temporal variation in temperature in a spark ignition natural gas engine. Nitrogen oxide, NO, computations were undertaken for the simplified finite-rate chemistry and the full finite-rate chemistry using temperature and pressure as inputs. It was shown that the simplified finite-rate chemistry model could not accurately capture the physically expected behavior of NO formation, whereas the full finite-rate chemistry did so accurately. In [12], equilibrium chemistry computations using the equilibrium constant were used to investigate the engine-out NOx, and the results were compared with the experimental data. Because NOx formation is believed to freeze a few crank angle degrees after the combustion process is complete, the computed results, that is, the equilibrium NOx concentration from the end of combustion to exhaust valve open, were compared with the experimental data in order to assess the applicability of equilibrium methods for NOx predictions. Based on the accuracy of the equilibrium calculations, it was surmised that this method could be used for evaluating the relative merits of nitrogen-enriched air with other NOx reduction techniques such as exhaust gas recirculation. This paper is organized as follows: Sect. 1 gives a general introduction concerning the formation of NOx and various methods of reducing it. Section 2 presents a literature review of some work that has been done regarding the simulation and modeling of NOx, while Sect. 3 gives a comprehensive description of the objective of this work. In Sect. 4, we laid a solid foundation on which this study is based by presenting in great detail the rudiments of chemical kinetics. The four primary sources of NOx in combustion processes, the extended Zeldovich mechanism, and its analysis were presented in Sect. 5. In Sect. 6, the mathematical model of the extended Zeldovich mechanism was used to investigate the formation of NO, N, and O under the influence of temperature and pressure. The parameters that were used to investigate

8289

the dynamics and stability of the formation of these products were the eigenvalues and timescales of the Jacobian matrix. The concerns on how the issues of uncertainties (local and global errors) were addressed are highlighted in Sect. 7. The results of our findings and some pertinent comments were presented in Sect. 8. Lastly, the conclusions from this study were itemized in Sect. 9.

3 Objective of the Present Study From the foregoing, we see that detailed knowledge of the reaction mechanisms of the reacting products would go a long way in enhancing combustion efficiency, thereby leading to optimal fuel savings, reduced NOx formation, and environmental benefits. In view of these potential benefits, it was decided to undertake a computational study of the effects of temperature and pressure on the formation of NOx by using the extended Zeldovich mechanism. The two temperatures considered in this study were 2,600 and 1,900 K while keeping the species concentration of all component products constant at 1.0 mole/m3 . Furthermore, parametric studies were carried out with a view to understanding the effect of pressure on the extended Zeldovich mechanism. The pressures considered were 1.513 × 105 and 1.362 × 105 N/m2 for temperature, T = 2,600 K and 1.106×105 and 9.95×104 N/m2 for T = 1,900 K. These pressure variations were realized by decreasing the species concentration for all components by 10 % while keeping the relevant temperature constant.

4 Chemical Kinetics Combustion will largely be sustained through only those chemical mixtures that are capable of reacting quickly enough to be considered explosive. The expression explosive implies very rapid reaction. With respect to combustion, the interest in chemical kinetic mechanisms has generally focused on the conditions under which chemical systems undergo explosive reaction. Recently, however, great interest has developed in the rates and mechanisms of steady (non-explosive) chemical reactions, since most of the known complex pollutants form in zones of steady, usually lower temperature, reactions during, and even after, the combustion process. A more comprehensive coverage of these features of chemical kinetics can be found in [13–16]. All chemical reactions, whether of the hydrolysis, acid– base, or combustion type, take place at a definite rate and depend on the conditions of the system. The most important of these conditions are the concentration of the reactants, the temperature, radiation effects, and the presence of a catalyst or inhibitor. The rate of the reaction may be expressed in terms of the concentration of any of the reacting substances

123

8290

Arab J Sci Eng (2014) 39:8287–8305

or of any reaction product, that is, the rate may be expressed as the rate of decrease in the concentration of a reactant or the rate of increase in a reaction product. 4.1 Law of Mass Action The discussions in this section are closely adapted in part from [17,18]. A single (one-step), forward chemical reaction of arbitrary complexity can be represented by the following stoichiometric equation: N 

kf

ν1 [Mi ] →

i=1

N 

νi [Mi ]

(1)

i=1

In Eq. (1) above, [Mi ] is the molar concentration of the ith chemical species, (i = 1, . . ., N ), N is the total number of species involved, and νi and νi are the corresponding forward and reverse stoichiometric (molar concentration) coefficients of the species, and the subscript f in k f designates the (forward) direction of the reaction. If a species represented by [Mi ] does not occur as a reactant or product, its ν j or ν j equals zero. The net stoichiometric (molar concentration) coefficient, νi , for each species is defined as νi = νi − νi

(2)

The rate of change in the molar concentration [Mi ] (moles per unit volume) of species i is given by d[Mi ] dt and is uniquely related to  ω j of species j by

 ωi =

(3)

  d Mj  ωj 1 =ω = ν j − ν j ν j − ν j dt

(4)

Since ω (moles per unit volume per second) is speciesindependent, it can be defined as the reaction rate of Eq. (1). The law of mass action, which has been experimentally confirmed, states that the rate of disappearance of a chemical species i, defined as ω, is proportional to the product of the concentrations of the reacting chemical species, where each concentration is raised to a power equal to the corresponding stoichiometric coefficient, that is, ω = k f (T )

N 

ν

[M]i i

(5)

i=1

where the proportionality factor k f (T ) is called the specific reaction rate constant. In Eq. (5), νi is also given the symbol, n which is called the overall order of the reaction; νi itself would be the order of the reaction with respect to species i.

123

where AT β represents the collision frequency and the exponential term is referred to as the Boltzmann factor. This factor specifies the fraction of collisions that have energy levels greater than the activation energy, E a. The values of A, β and E a depend on the nature of the elementary reaction. For a given chemical reaction, these parameters are neither functions of the concentration nor functions of temperature. The exponent β lies between 0 and 1, that is, (0 < β < 1). In an actual reacting system, the rate of change in the concentration of a given species i is given by   ν  

 d[Mi ]  Mj j = νi − νi ω = νi − νi k f (T ) dt N

(7)

i=1

since ν j moles of [Mi ] are formed for every ν j moles of [Mi ] consumed. Equation (5) is based on the microscopic viewpoint that the reaction rate is proportional to the collision frequency, which in turn is proportional to the product of the species concentrations. Thus, implicit in Eq. (5) is the requirement that reaction (1) is an elementary one in which the reaction actually takes place at the molecular level. As an example of an elementary reaction, we have H + HO2 → OH + OH which could be an important step in the oxidation of hydrogen and hydrocarbons. With reference to the chemical reaction stated above vis-a-vis, Eq. (7)

 ωi 1 d[Mi ] =    νi − νi νi − νi dt =

For a given chemical reaction, k f (T ) is independent of concentrations [Mi ], and it is only a function of temperature. In general, k f (T ) is given as  Ea (6) k f (T ) = AT β exp − RT

N = 3, M1 = H, M2 = HO2 and M3 = OH ν1 = 1, ν2 = 1, and ν3 = 0

ν1 = 0, ν2 = 0, and ν3 = 2

Then, for this reaction, we can write d [H] d [HO2 ] 1 d [OH] =− = dt dt 2 dt while, from the law of mass action, the reaction rate is given by ω=−

ω = k f [H] [HO2 ]. 4.2 Reversible Reactions Associated with every forward reaction (1) is the corresponding backward reaction N  i=1

kr

νi [Mi ] →

N  i=1

ν  [Mi ]

(8)

Arab J Sci Eng (2014) 39:8287–8305

8291

It follows then that the net reaction rate in the presence of both forward and backward reactions is given by 



ωi, f +  ωi,b = νi − νi ω f − ωr  ωi = 

 = νi − νi ω = νi ω (9) with ω given by ω = kf

N 

ν

[M]i i − kr

i=1

N 

ν 

[M]i i

(10)

i=1

Only one of k f and kr needs to be determined. This is because at equilibrium, ω ≡ 0, the rate of forward reaction is balanced by that of backward reaction, and Eq. (10) yields  kf (ν  −ν  ) = [M]i i i kr N

(11)

i=1

By definition, Eq. (11) is just the equilibrium constant K c , which is given by  K p (T ) kf (ν  −ν  ) = K c (T ) = [M]i i i = N   kr (RT ) i=1 (νi −νi ) i=1

N 

k j, f

νi, j [Mi ] →

i=1 N 

νi, j [Mi ]

(16a)

νi, j [Mi ]

(16b)

i=1 k j,r

νi, j [Mi ] →

i=1

(12)

where Ru = 8.314416 J/(mol-K) is the universal gas constant and T is the temperature in Kelvin. If we substitute Eq. (12) into Eq. (10), we have N

N  1  νi νi [M]i − [M]i (13) ω = kf Kc

N  i=1

where ΔG 0 is the standard Gibbs free energy change in the reaction, R is the gas constant, and T is the absolute temperature. From Eqs. (2), (3), (4), (6), (9), (10), and (13), we have

N N   d[Mi ] νi νi  [M]i − kr [M]i = νi k f  ωi = dt i=1 i=1 N

N  1  νi νi  [M]i − [M]i = νi k f (15) Kc i=1

Equation (15) is called the law of mass action. 4.3 Multistep Reactions There are seldom reactions in which the original reactants interact with each other in a single step at the molecular level to produce the final products shown in Eq. (1). A large number

  d[M]i νi, j − νi, j ω j = dx N

 ωi =

(17)

i=1

such that ω j = k j, f

N N   ν  ν , j [M]i i − k j,r [M]i i, j i=1

i=1

Since it is easier to determine the equilibrium constant, K c to a much greater accuracy than kr, Eq. (13) is usually preferred to Eq. (10). It can be shown thermodynamically that the equilibrium constant, K c , is given by  ΔG 0 (14) K c = exp − RT

i=1

N 

and the generalized law of mass action is given by

N

i=1

of intermediate elementary steps with intermediate species are actually involved before the final products are formed. Now, consider a volume V in which a reactive mixture consisting of N molecular species composed of L atomic elements are interacting in J elementary reactions. Let us assume that the gas is an ideal mixture of ideal gases that satisfies Dalton’s law of partial pressures. The reaction will be considered to be driven by molecular collisions. We will not model individual collisions, but instead attempt to capture their collective effect. A compact notation of the reaction mechanism described above is given by

( j = 1, 2, . . ., J ) ,

i=1

(18) It should be noted that in the above, the subscript for  ωi is for species i, whereas that for ω j is for reaction j. The law of conservation of mass which was established in 1789 by French Chemist Antoine Lavoisier states that in any chemical reaction process (in the absence of nuclear reactions), the number of moles of each element in each reaction is conserved. This molar balance is enforced by the following stoichiometric constraint on element l in reaction j, thus N 

ϕl,i νi, j = 0, for (l = 1, 2, . . ., L), and ( j = 1, 2, . . ., J )

i=1

(19) where νi, j = νi, j − νi, j is the stoichiometric matrix, which gives the net stoichiometric coefficient for species i in reaction j. Equation (19) demands that νi, j lie in the right null space of ϕl,i . Physically, it means the mass and number of moles of each element are conserved in each reaction. In general, ϕl,i and νi, j are non-square matrices of dimensions L × N and N × J , respectively. However, ϕl,i is of full rank, while νi, j is of rank R ≤ (N − L). This implies that the chemical interaction between the N species can be described using R species.

123

8292

Arab J Sci Eng (2014) 39:8287–8305

5 Mechanisms of Oxides of Nitrogen Formation The combustion chemistry of nitrogen compounds has generated a lot of interest because of its role in the mechanism of formation of oxides of nitrogen, which are collectively referred to as NOx and its adverse effects on our environment and health. The increasing stringent environmental regulations have necessitated the introduction of highly sophisticated approaches for the control/amelioration of NOx emissions. A detailed understanding of chemical kinetics modeling plays a major role in the development and implementation of such approaches. Furthermore, another major motivation for the study of NOx formation mechanism stems from the fact that energetic materials, such as explosives, always contain nitrogen compounds and oxides of nitrogen are constantly present in their combustion products. In view of the foregoing, the knowledge of NOx formation mechanism is essential for understanding combustion of energetic materials. N O is known to be formed in a variety of ways, namely 1. Thermal NO which was first extensively studied by Zeldovich [19] is primarily a consequence of high-flame temperatures; 2. Fenimore [20] postulated that prompt NO is generated in fuel-rich parts of the flames; 3. Malte and Pratt [21] stated that N2 O mechanism can be important in high-pressure flames; 4. Fuel nitrogen NO [20] results from converting nitrogencontaining compounds in the fuel into NO, and; 5. Dean and Bozzelli [22] found that NNH mechanism is active in flame fronts where high atom concentrations appear. All the mechanisms stated above except number 4 convert atmospheric nitrogen into NO and relate to all fuels, while pathway (4) pertains to flames of petroleum, coal, and biomass. 5.1 Thermal or Extended Zeldovich NO Mechanism In combustion with clean fuels (which do not contain nitrogen compounds) in the presence of air which contains atmospheric N2 , N O is formed mainly by the Zeldovich mechanism. Detailed descriptions of all the mechanisms can be found in [20–23]. Thermal NOx refers to NOx formed through hightemperature oxidation of the diatomic nitrogen found in combustion air. The formation rate is primarily a function of temperature and the residence time of nitrogen at that temperature. At high temperatures, usually above 1,600◦ C (2,900◦ F), molecular nitrogen (N2 ) and oxygen (O2 ) in the

123

combustion air disassociate into their atomic states and participate in a series of chain reactions. The six principal elementary reactions of the extended Zeldovich mechanism that produce thermal NOx are k1 f

N + NO → N2 + O k1r

N2 + O → N+NO k2 f

N+O2 → NO+O k2r

NO+O → N+O2 k2 f

N+OH → NO+H k2r

NO+H → N+OH

(20a) (20b) (21a) (21b) (22a) (22b)

where k1 f = 1.8 × 1011 exp(−38,370/T )m3 /(kmol · s) k1r = 3.8 × 1010 exp(−425/T )m3 /(kmol · s) k2 f = 1.8 × 107 exp(−4,680/T )m3 /(kmol · s) k2r = 3.81 × 106 exp(−20,820/T )m3 /(kmol · s) k3 f = 7.1 × 1010 exp(−450/T )m3 /(kmol · s) k3r = 1.7 × 1011 exp(−24,560/T )m3 /(kmol · s) The forward global mechanism representing the system of kinetic reactions shown in Eqs. (20)–(22) is given by kg

N2 + O2 →2NO

(23)

where k g is the global reaction constant. Reactions (20), (21), and (22) are reversible, the letters f and r in the subscript of reaction constant k represent the forward and reverse directions, respectively. Zeldovich was the first to suggest the importance of the first two reactions. The last reaction of atomic nitrogen with the hydroxyl radical (HO− ) was added by Lavoie et al. [24] to the mechanism. This was a significant contribution to the knowledge of thermal NOx formation. Reactions (20) and (21) constitute a chain sequence, such that a small amount of atomic oxygen can produce large quantities of N O. The first reaction (20) is rate-limiting because of its high activation energy, namely 320 kJ/mol. This is due to the strong triple bond in the N2 , and so it is sufficiently fast only at very high temperatures. As a rule of thumb, the thermal mechanism is usually unimportant at temperatures below 1,800 K. Comment: The Zeldovich mechanism is highly sensitive to temperature not only because of the high activation energy of reaction (20) but also as a result of the rapidly increasing concentration of oxygen atoms in flames with increasing temperature. Since the rate coefficient of reaction (20) is well known, computation of the Zeldovich NO production rate is straightforward if the oxygen atom concentration

Arab J Sci Eng (2014) 39:8287–8305

is established. At regions far from the flame front, the Zeldovich NO production can be computed with little concern for details of the chemical kinetics so long as the equilibrium of the nitrogen atom concentrations is assumed. The situation is more complicated in the flame-front region where a preponderance of equilibrium concentrations of radicals is present. This does not have anything to do with any particular flame-front reactions but because flame-front chemistry is driven by fast bimolecular reactions while post-flame chemistry includes slow termolecular recombination reactions. It is absolutely necessary to properly describe the kinetics of flame-front chemistry in order to predict Zeldovich NO production accurately. Nitrogen and oxygen molecules are preferred at low temperature. As the temperature rises, N and O begin to appear. It is possible when they are mixed for NO to appear as a product. An examination of the elementary reactions in Eqs. (20)– (22) indicates the presence of seven intermediate species in the formation of nitric oxide NO, namely NO, N, N2 , O2 , O, H, and OH. Note that the global expression in Eq. (23) that defines the initial reactions and the final products fails to account for these intermediate species and is therefore not representative of the actual chemical process. 5.2 Mathematical Model The details of the mathematical model used to describe the version of the extended Zeldovich mechanism employed in this study are presented in this section. From the six elementary reactions in Eqs. (20)–(22) and the law of mass action, it is possible to obtain the elementary reaction rates for the extended Zeldovich mechanism thus: Let r = (r1r2 r3 )T where  1 r1 = k1 f [N][NO] − [N2 ][O] (24) K c1  1 [NO][O] (25) r2 = k2 f [N][O2 ] − K c2  1 [NO][H] (26) r3 = k3 f [N][OH] − K c3 and T is the matrix transpose operator and K c1 , K c2 and K c3 are as defined in Eq. (14). If we define the seven (N = 7) species present in the six (J = 6) reactions, of Eqs. (20)–(22) as species i = 1 ≡ [NO], species i = 2 ≡ [N], species i = 3 ≡ [N2 ], species i = 4 ≡ [O], species i = 5 ≡ [O2 ], species i = 6 ≡ [H], and species i = 7 ≡ [OH], then for reactions, j = 1, 2, . . ., J , we have the net stoichiometric coefficient (νi, j − νi, j ) for   = = 0, νi, j = ν1,1 species [NO], [i = 1], as νi, j = ν1,1       1, ⇒ ν1,1 − ν1,1 = 1, νi, j = ν1,2 = 0, νi, j = ν1,2 = 1, ⇒

8293  − ν      ν1,2 1,2 = 1 and νi, j = ν1,3 = 0, νi, j = ν1,3 = 1, ⇒   ν1,3 − ν1,3 = 1. We can use a similar procedure to determine the values of (νi, j − νi, j ) for the remaining species. The final result is given by the elements of the matrix ν, while the (N = 7) species concentration, [X ], is given by the vector form shown below

⎞ −1 1 1 ⎜ −1 −1 −1 ⎟ ⎟ ⎜ ⎜ 1 0 0 ⎟ ⎟ ⎜ 1 0 ⎟ ν=⎜ ⎟ ⎜ 1 ⎜ 0 −1 0 ⎟ ⎟ ⎜ ⎝ 0 0 1 ⎠ 0 0 −1 ⎛

⎛

⎞ [NO] ⎜ [N] ⎟ ⎜ ⎟ ⎜ [N2 ] ⎟ ⎜ ⎟ ⎟ [X] = ⎜ ⎜ [O] ⎟ ⎜ [O2 ] ⎟ ⎜ ⎟ ⎝ [H] ⎠ [OH]

(27)

The dimension of matrix ν is N × J = 7 × 3. Furthermore, the matrix ν has rank R = 3. In the extended Zeldovich mechanism adopted for this study, we have L = 3 elements, with N , (l = 1), O, (l = 2), and H (l = 3). The stoichiometric matrix, ϕ whose elements are ϕl,i , has a dimension of L × N = 3 × 7, and it is given by ⎛

ϕ1,1 ϕ = ⎝ ϕ2,1 ϕ3,1 ⎛ 1 1 = ⎝1 0 0 0

⎞ ϕ1,2 ϕ1,3 ϕ1,4 ϕ1,5 ϕ1,6 ϕ1,7 ϕ2,2 ϕ2,3 ϕ2,4 ϕ2,5 ϕ2,6 ϕ2,7 ⎠ ϕ3,2 ϕ3,3 ϕ3,4 ϕ3,5 ϕ3,6 ϕ3,7 ⎞ 2 0 0 0 0 0 1 2 0 1⎠ 0 0 0 1 1

(28)

It is instructive to note that the stoichiometric constraint on element conservation for each reaction as expressed in Eq. (19), that is, ϕ · ν = 0 holds ⎞ 1 1 1 ⎟ ⎜ ⎞ ⎜ 1 −1 −1 ⎟ ⎜ 0 ⎟ 2 0 0 0 0 ⎜ −1 0 ⎟ ⎜ ⎠ 0 ⎟ 0 1 2 0 1 ⎜ −1 1 ⎟ ⎟ 0 0 0 1 1 ⎜ ⎜ 0 −1 0 ⎟ ⎝ 0 0 1 ⎠ 0 0 −1 ⎞ 0 0⎠ 0 ⎛

⎛

1 1 ϕ · ν = ⎝1 0 0 0 ⎛

0 0 = ⎝0 0 0 0

(29)

The resulting matrix is a 3×3 matrix of zero elements because there are three reactions each with 3-element constraints. Mathematically, we say that each of the column vectors of ν lies in the right null space of ϕ. From the foregoing, the nonlinear differential equations modeling the extended Zeldovich reaction dynamics can be expressed as

123

8294

Arab J Sci Eng (2014) 39:8287–8305

d[X] =ν·r dt ⎞ ⎛ ⎞ ⎛ −1 1 1 [NO] ⎜ [N] ⎟ ⎜ −1 −1 −1 ⎟ ⎟ ⎜ ⎟ ⎜ ⎜ [O] ⎟ ⎜ 1 0 0 ⎟ ⎟ ⎜ ⎜ ⎟ d 1 1 0 ⎟ =⎜ [N2 ] ⎟ × ⎜ ⎟ ⎜ ⎜ ⎟ dt ⎜ ⎟ ⎟ ⎜ ⎜ [O2 ] ⎟ ⎜ 0 −1 0 ⎟ ⎝ [H] ⎠ ⎝ 0 0 1 ⎠ 0 0 −1 [OH] ⎛  ⎞ 1 k1 f [N][NO] − K c1 [N2 ][O] ⎜  ⎟ ⎜ ⎟ × ⎜ k2 f [N][O2 ] − K1c2 [NO][O] ⎟ ⎝  ⎠ k3 f [N][OH] − K1c3 [NO][H]

(30)

From the nonlinear differential equations in Eq. (30), it can be seen that the reaction rate for each intermediate species is a nonlinear ordinary differential equation (ODE) whose order is equal to the largest sum of the reaction orders for each term of the rate equation. For the extended Zeldovich mechanism, the largest sum of the reaction orders for each of the seven reaction rates is two; therefore, the overall order of reaction for each rate equation in this system is two. In order to find all the linear dependencies, we transform the model, namely Eq. (30), to the equivalent reduced row echelon form. The transformed reduced row echelon form is given by ⎛

⎞ ⎛ ⎞ [NO] −1/2 0 0 0 0 0 ⎜ ⎟ 1/2 1 0 0 0 0⎟ ⎟ ⎜ [N] ⎟ ⎜ [O] ⎟ −1 −1 0 0 0 0 ⎟ ⎟ d ⎜ ⎟ ⎜ ⎟ 1/2 0 1 0 0 0⎟ ⎟ dt ⎜ [N2 ] ⎟ ⎟ ⎜ 1/2 1 0 1 0 0 ⎟ ⎜ [O2 ] ⎟ ⎟ 1 1 0 0 1 0 ⎠ ⎝ [H] ⎠ −1 −1 0 0 0 1 [OH] ⎞ 1 0 0  ⎞ ⎜0 1 0⎟⎛ ⎟ k1 f [N][NO] − 1 [N2 ][O] ⎜ ⎜0 0 1⎟⎜ K c1  ⎟ ⎟⎜ ⎜ ⎟ ⎜ k2 f [N][O2 ] − 1 [NO][O] ⎟ 0 0 0 =⎜ ⎟ ⎟⎝ ⎜ K c2  ⎠ ⎜0 0 0⎟ ⎟ k3 f [N][OH] − 1 [NO][H] ⎜ K c3 ⎝0 0 0⎠ 0 0 0   

−1/2 ⎜ 1/2 ⎜ ⎜ 0 ⎜ ⎜ 1/2 ⎜ ⎜ 1/2 ⎜ ⎝ 0 0 ⎛

(32)

1 1 [NO] + [N] + [O] + [O2 ] = C2 2 2

(33)

[N] + [O] + [H] = C3

(34)

−[N] − [O] + [OH] = C4

(35)

The constants C1 , C2 , C3 , and C4 can be determined from the initial conditions on all seven state variables. When Eqs. (32)–(35) are expressed in matrix form, we have ⎞ ⎛ [NO] ⎛ ⎞ ⎜ [N] ⎟ ⎛ ⎞ ⎟ C1 1/2 1/2 0 1 0 0 0 ⎜ ⎟ ⎜ ⎜ 1/2 1/2 1 0 1 0 0 ⎟ ⎜ [O] ⎟ ⎜ C2 ⎟ ⎜ ⎟ ⎜ [N2 ] ⎟ = ⎜ ⎟ (36) ⎟ ⎝ C3 ⎠ ⎝ 0 1 1 0 0 1 0⎠⎜ ⎜ [O2 ] ⎟ ⎟ ⎜ C4 0 −1 −1 0 0 0 1 ⎝ [H] ⎠ [OH]

Assuming the three free variables, namely [NO], [N], and [O] are known, we can move them to the right-hand side of Eq.(36) to get ⎛ ⎞⎛ ⎞ 1 0 0 0 [N2 ] ⎜ 1 1 0 0 ⎟ ⎜ [O2 ] ⎟ ⎜ ⎟⎜ ⎟ ⎝ 0 0 1 0 ⎠ ⎝ [H2 ] ⎠ 0 0 0 1 [OH] ⎛ ⎞ C1 − 21 [NO] − 21 [N] ⎜ C2 − 1 [NO] − 1 [NO] − 1 [N] − [O] ⎟ ⎟ 2 2 2 (37) =⎜ ⎝ ⎠ C3 − [N] − [O] C4 + [N] + [O]

On solving Eq. (37) for the bound variables, we have ⎛ ⎞ ⎞ ⎛ [N2 ] C1 − 21 [NO] − 21 [N] ⎜ [O2 ] ⎟ ⎜ C2 − 1 [NO] − 1 [NO] − 1 [N] − [O] ⎟ ⎜ ⎟ ⎟ ⎜ 2 2 2 ⎝ [H2 ] ⎠ = ⎝ ⎠ C3 − [N] − [O] [OH] C4 + [N] + [O] (38)

=U

(31) Equation (31) defines three free variables, those associated with the nonzero pivots of matrix U in Eq. (31), namely [NO], [N], and [O]. The remaining four variables [N2 ], [O2 ], [H], and [OH] are bound variables, which can be expressed in terms of the free variables. The last four ordinary differential equations in Eq. (31) are homogeneous and can be integrated to get

123

1 1 [NO] + [N] + [N2 ] = C1 2 2

Since we considered an isothermal, isochoric, and molepreserving reaction, it will also be isobaric. We shall now analyze two representative scenarios. We shall consider an isothermal reaction at T = 2,600 K and T = 1,900 K. These values were considered because they represent typical operating temperatures for internal combustion engines (spark and compression ignition engines). We shall take the initial condition as [NO]i = [N2 ]i = [Ni = [2 ]i = [O]i = [H]i = [OH]i = 1.0 moles/m3 . Therefore, C1 = 2, C2 = 3, C3 = 3, and C4 = −1. For this temperature and concentrations, the pressure, which will remain constant throughout the reaction, is given by

Arab J Sci Eng (2014) 39:8287–8305

P = RT

n 

8295

[M j ] ⇒ P = RT ([NO]i + [N]i + [N2 ]i

j=1

+[O2 ]i + [H]i + [OH]i )   mol J × 2,600(K) × 7 × 1.0 × = 8.31446 K · mol m3 N = 1.513 × 105 2 (39) m

⎛ ⎞ ⎛ ⎞ [NO] −1 1 1 d ⎝ [N] ⎠ = ⎝ −1 −1 −1 ⎠ dt [O] 1 1 0  ⎛ ⎞ k1 f [N][NO] − K1c1 [N2 ] [O] ⎜  ⎟ ⎟ ⎜ × ⎜ k2 f [N] [O2 ] − K1c2 [NO][O] ⎟ ⎝  ⎠ k3 f [N] [OH] − K1c3 [NO] [H]

(43)

where [M j ] is the number of moles per unit volume (concentration) of species, j. Notice that the pressure is a little greater than atmospheric. The thermodynamic data required to calculate g N O, g N , g N2 , g O , g O2 , g H , and g O H ; and consequently, K c1 , K c2 , and K c3 were taken from [25]. For our system at T = 2,600 K, we find g N O = −564,722.2, g N = 4,859.4, g N2 = −602,244.2, g O = −240,061, g O2 = −642,730.8, g H = −149,482.2 and g O H = −539,708. All values are in kJ/kmol. For each reaction, the corresponding G cj is given by

On substituting values for C1 = 2, C2 = 3, C3 = 3, and C4 = −1 and K c1 , K c2 and K c3 , k1 f , k2 f and k3 f, we have

G = g · ν (40) ⎞ ⎛ G c1 

× ⎝ G c2 ⎠ = g N O g N g N2 g O g O2 g H g O H G c3 ⎛ ⎞ −1 1 1 ⎜ −1 −1 −1 ⎟ ⎜ ⎟ ⎜ 1 0 0 ⎟ ⎜ ⎟ 1 0 ⎟ ×⎜ (41) ⎜ 1 ⎟ ⎜ 0 −1 0 ⎟ ⎜ ⎟ ⎝ 0 0 1 ⎠ 0 0 −1

− 51980188.96[N][O] − 55848182.93[N]2 (45)

Therefore, for each reaction, in Eqs. (24)–(26), we can compute the corresponding G cj , where ( j = 1, 2, 3) as follows: G c1 = −g N O − g N + g O = −282,442.4 J/mol G c2 = g N O − g N + g O − g O2 = −166911.8 J/mol G c3 = g N O − g N + g H − g O H = −179355.8 J/mol At T = 2,600 K, the calculated equilibrium constants, K c1 , K c2 , and K c3 for the J = 3 reactions are   G c1 282442.4 K c1 = exp − = exp = 472355.83 RT (8.3144) · (2600)   G c2 166911.8 = 2255.58 K c2 = exp − = exp RT (8.3144) · (2600)   G c3 179355.8 = 4011.06 K c3 = exp − = exp RT (8.3144) · (2600)

(42) From Eqs. (31)–(42), we find that the three differential equations governing the evolution of the free variables are given by

d[NO] = −44607.99[NO] − 36508213.02[N] dt + 0.0003[O] − 3853194.83[NO][N] + 11442.29[NO][O] + 51980188.96[N][O] + 55848182.93[N]2 (44) d[N] = 44607.99[NO] + 36508213.02[N] + 0.0003[O] dt + 3853194.83[NO][N] − 11442.29[NO][O] d[O] = 23207963.91[N] − 0.0003[O] dt − 3867923.84[NO][N] − 3427.04[NO][O] − 7735987.97[N][O] − 3867994[N]2

(46)

5.3 Specification of Solution Parameters In order to solve for the equilibrium solutions to Eqs. (44)– (46) and (63)–(65), the modified Rosenbrock triple was used. The complete details concerning this algorithm can be found in [29]. Relative error tolerance, r tol of 1 × 10−6 , and absolute error tolerance, atol of 1 × 10−6 , for each component were found to yield satisfactory results as will be discussed in Sect. 8.3 below. The solution interval used was from time, t = 0 to 1 × 10−6 seconds. The initial conditions of [NO], [N ], and [O] were set to [1 1 1]T , respectively. Since actual (global) errors may exceed these local tolerances, the values of r tol and atol stated above were conservatively chosen (as stated above) so as to ensure that global errors do not exceed local tolerances. The computational cost statistics are shown in Tables 15 and 16 for two conditions investigated. On solving Eqs. (44)–(46) based on the specifications given above, we obtain the solution shown in Fig. 1. The computations show a relaxation to final concentration values of mole t→∞ m3 mole lim [N] = 3.333 × 10−3 3 t→∞ m 0 mole lim [O] = 1.203 × 10 t→∞ m3 lim [NO] = 1.997 × 100

(47) (48) (49)

123

8296

Arab J Sci Eng (2014) 39:8287–8305

In this study, we will use local linear analysis in the neighborhood of each equilibria to rigorously ascertain the stability of each root.

2

Concentration, moles/m

3

[NO] [N] [O]

1.5

5.4 Analysis of Physical Roots and Initial State Taylor series expansion of Eqs. (50)–(52) in the neighborhood of an equilibrium point yields

1

d ∂ f [NO] ([NO] − [NO]e ) = f ([NO]e ) + ([NO] − [NO]e )    ∂[NO]e dt

0.5

=0

0

0

0.2

0.4

0.6

0.8

time, (s)

= 1.513 × 105

/

2

1 x 10

,

-6

= 2600

Fig. 1 Temporal variation in the concentrations of [NO], [N], and [O]

∂ f [NO] ∂ f [NO] + ([N] − [N]e ) + ([O] − [O]e ) + · · · (56) ∂[N]e ∂[O]e ∂ f [N] d ([NO] − [NO]e ) ([N] − [N]e ) = f ([N]e )=0 +    ∂[NO]e dt ∂ f [N] ∂ f [N] + (57) ([N] − [N]e ) + ([O] − [O]e ) + · · · ∂[N]e ∂[O]e ∂ f [O] d ([NO] − [NO]e ) ([O] − [O]e ) = f ([O]e ) +    ∂[NO]e dt =0

∂ f [O] ∂ f [O] + ([N] − [N]e ) + ([O] − [O]e ) + · · · ∂[N]e ∂[O]e

Eqs. (44)–(46) are of the forms d[NO] = f [NO] ([NO], [N], [O]) dt d[N] = f [N] ([NO], [N], [O]) [dt d[O] = f [O] ([NO], [N], [O]) dt

(50) (51) (52)

At equilibrium, we have f [NO] ([NO], [N], [O]) = 0

(53)

f [N] ([NO], [N], [O]) = 0

(54)

f [O] ([NO], [N], [O]) = 0

(55)

(58)

It is worth noting that the subscript e in Eqs. (56)–(58) denotes equilibrium point. Evaluation of Eqs. (56)–(58) near the physical root, namely root 2, yields the system of differential equations given by ⎞ ⎛ [NO] − 1.1278 · 10−1 d ⎝ −5 ⎠ [N] − 6.1645 · 10 dt [O] − 1.6259 · 100 ⎞ ⎛ ⎞⎛ −2641.5 44596703.10 13350.00 [NO] − 1.1278 · 10−1 −5 ⎠ = ⎝ 2641.5 −44596703.10 −13350.00 ⎠ ⎝ [N] − 6.1645 · 10 [O] − 1.6259 · 100 −5810.46 7199925.60 −3515.58    ∂f =J= ∂x

(59) The critical points (that is, equilibrium Solutions) of the system described by Eqs. (53)–(55) are the constants and are time-independent for a given set of initial conditions. The critical points may vary with initial conditions. Because of negative concentrations, roots 3–6 are nonphysical and root 1 is the trivial case. However, root 2 is physically realizable. These results are tabulated in Table 1. Table 1 Equilibrium points, T = 2,600 K

123

Equilibrium points

[NO] mole/m3

where x = ([NO][N][O])T . Equation (59) is of the form d ∂f ([x] − [x]e ) = J · ([x] − [x]e ) ([x] − [x]e ) = dt ∂ [x]e (60)

[N] mole/m3

[O] mole/m3

Remark

1 2

0 1.1278 ×10−1

0 6.1645 ×10−5

0 1.6259 ×100

3

−3.1757 × 10−1

9.8153 ×10−5

−7.2898 × 100

Non-physical

4

−1.1663 × 10−1

−6.0457 × 10−5

×100

Non-physical

5

3.0558 ×10−1

−9.5532 × 10−5

−6.8273 × 100

Non-physical

6

−5.3449 × 10−6

×100

Non-physical

−4.0000

× 100

1.6490 5.0000

Trivial Physical

Arab J Sci Eng (2014) 39:8287–8305

8297

It is the eigenvalues of the Jacobian matrix Jm that give the timescales of evolution of the concentrations as well as determine the stability/dynamics of the local equilibrium point. Furthermore, it can be shown that the solution to component of Eq. (60) are given by [x]1 = C1 exp(λ1 t), [x]2 = C2 exp(λ2 t), [x]3 = C3 exp(λ3 t), . . .

Therefore, the initial stiffness, κ = 7.13 × 10−3 /8.08 × 10−9 = 8.82 × 105 , is less than that of the time scale of the physical equilibrium point. It can be seen from Fig. 1a, b that this initial timescale of the order (10−9 ∼ 10−7 s) well predicts where significant evolution of species concentrations commences.

In this study, the Jacobian matrix Jm of Eq. (59) has three real, negative eigenvalues in the neighborhood of the physical root 2. The eigenvalues of the Jacobian matrix Jm are

5.5 Analysis of Non-Physical Roots, P = 1.513 × 105 N/m2 and T = 2,600 K

λ1 = −4.4621 × 107 s−1 , λ2 = −9.2302 × 10−2 s−1 and

A review of the modes of stability of the non-physical roots, 3 to 6, are considered in this section. Analyses similar to those of Eq. (59) were performed to determine the eigenvalues of the Jacobian matrix at the relevant critical points. These eigenvalues can be utilized to classify the critical points and to determine their stability according to their local behavior. Analysis shows that all four non-physical equilibrium points have real eigenvalues. The eigenvalues, time constants, and the corresponding classification of the type of stability of each critical point are summarized in Tables 2, 3, and 4.

λ3 = −5.6716 × 103 s−1 The value of the eigenvalues shows that the physical equilibrium point is linearly stable. The local time constants, τ1 , τ2 , and τ3 near equilibrium are given by the reciprocal of the magnitude of the eigenvalues. These are given by 1 1 = 2.24 × 10−8 s, τ2 = = 1.08 × 101 s and |λ1 | |λ2 | 1 τ3 = = 1.76 × 10−4 s |λ3 |

τ1 =

The timescales show that this is in fact a multiscale problem. One of the major difficulties in the numerical simulation of combustion problems is the capturing of the effects at all relevant scales (time and length). The problem is made more difficult as the breadth of the scales expands. In this problem, the breadth of scales is highly challenging. Near equilibrium, the ratio of the slowest to the fastest timescale is called the stiffness κ, and it is given by κ=

τ2 1.08 × 101 τmax = = = 4.48 × 108 τmin τ1 2.24 × 10−8

Comments The six roots of Eqs. (53) to (55) define the equilibrium points of the extended Zeldovich mechanism. All of the six points are real and contain no imaginary components. One of the points, namely point 2, is a stable equilibrium point, while points 3, 4, 5, and 6 are all unstable, and point 1 is the trivial case since there can be no zero [NO], [N], and [O] at equilibrium, refer to Figs. 1 and 3. It is worth mentioning that, of the six equilibrium points, only one point, that is, (point 2), represents the physical dynamics of the extended

(61)

A lot of studies have shown that many combustion problems can have stiffness greater than 106 . This is more prevalent at lower temperatures. The cases we investigated in this work could be classified as been in a relatively low-temperature regime/range. Hence, our stiffness of 4.48 × 108 is quite consistent with the expected values. A similar linearization near the initial state gives the following local eigenvalues: λ1 = −4.4032 × 102 s−1 , λ2 = −1.1209 × 107 s−1 and

Table 2 Eigenvalues of non-physical equilibrium points, T = 2,600 K Equilibrium point

λ1 (1/s)

λ2 (1/s)

λ3 (1 /s)

3

4.14 ×108

−2.70 × 10−2

6.16 ×102

4

−4.97 × 107

1.67

×10−2

1.85 ×103

5

3.93

×108

2.74

×10−2

−5.88 × 102

6

2.81 ×102

3.10 ×107

2.39 ×108

λ3 = −1.2373 × 108 s−1

Table 3 Local time constants near non-physical equilibrium points, T = 2,600 K

Notice that the dynamics of the system is linearly stable near the initial state as well, and the corresponding local timescales are 1 1 = 7.13 × 10−3 s, τ2 = = 8.92 × 10−8 s and τ1 = |λ1 | |λ2 | 1 τ3 = = 8.08 × 10−9 s |λ3 |

Equilibrium point

τ1 (s)

τ2 (s)

τ3 (s)

3

2.42 ×10−9

3.70 ×101

1.63 ×10−3

4

2.01

×10−8

5

2.55

×10−9

6

3.56 ×10−3

6.60

×10−1

5.40 ×10−4

3.66

×101

1.70 ×10−3

3.23 ×10−8

4.19 ×10−9

123

8298

Arab J Sci Eng (2014) 39:8287–8305

Table 4 Classification of critical points, T = 2,600 K Equilibrium point

Type/sign of eigenvalues

Mode of stability

2

Sink node–negative values

Linearly stable

3

Saddle point–positive/negative

Unstable

4

Saddle point–positive/negative

Unstable

5

Saddle point–positive/negative

Unstable

6

Source node–positive

Unstable

Zeldovich mechanism. The other four points contain at least one negative component. Negative components imply negative concentrations, which is an absurdity. There can be no negative concentrations in the physical system at any time throughout the stages of the mechanism. One of the key challenges in computational chemical kinetics is accurately predicting species concentration evolution with time. The problem is made more intractable because of the prevalence of physical phenomena that evolve over a widely different set of temporal and length scales. Systems which evolve over a wide range of scales (time and length) are known as stiff, and studying them numerically has been found to be non-trivial.

6 Effects of Temperature and Pressure In this section, we examined the effect of temperature and pressure on timescales and stiffness of the extended Zeldovich mechanism.

matical model describing the extended Zeldovich mechanism at temperature, T = 1,900 K, is given by d[NO] = −1, 408.53[NO] − 47289238.46[N] dt +5.97 × 10−9 [O] − 1455862.25[NO][N] +338.15[NO][O] + 53114564.26[N][O] +54570895.71[N]2

−338.15[NO][O] − 53114564.26[N][O] −54570895.71[N]2

−131.36[NO][O] − 2912662.9[N][O] −1456331.45[N]2

n  [M j ] ⇒ P = RT ([NO]i + [N]i + [N2 ]i

(65)

Solving Eqs. (63) to (65) numerically for T = 1,900 K, we obtain the solution shown in Fig. 2. The computation shows a convergence to final species concentration values of mole m3 mole lim [N] = 5.502 × 10−1 3 t→∞ m −1 mole lim [O] = 3.647 × 10 t→∞ m3 lim [NO] = 1.450 × 100

(66)

t→∞

(67) (68)

2

[NO] [N] [O]

3

Concentration, moles/m

P = RT

(64)

d[O] = 8737988.7[N] − 5.97 × 10−9 [O] dt −1456331.15[NO][N]

6.1 Effect of Temperature: P = 1.106 × 105 N/m2 and T = 1,900 K In this section, we examined the behavior of the extended Zeldovich mechanism at a lower temperature of T = 1,900 K while keeping all other parameters, including the initial species concentrations constant as in the high-temperature case. When the temperature is reduced, there is a corresponding decrease in pressure. In this situation, the new pressure is given by

(63)

d[N ] = 1408.53[NO] + 47289238.46[N ] dt +5.97 × 10−9 [O] + 1455861.64[NO][N]

1.5

1

0.5

j=1

+ [O2 ]i + [H]i + [OH]i )   J mol = 8.31446 × 1, 900 (K ) × 7 × 1.0 × K · mol m3 N = 1.106 × 105 2 (62) m

which is close to atmospheric pressure. Following considerations similar to the derivation of Eqs. (44)–(46), the mathe-

123

0

0

0.2

0.4

0.6

time, (s) = 1.106 × 105

2

,

0.8

1 x 10

-6

= 1900

Fig. 2 Temporal variation in the concentrations of [NO], [N], and [O]

Arab J Sci Eng (2014) 39:8287–8305 Table 5 Equilibrium points, T = 1,900 K

8299

[NO] mole/m3

Equilibrium points

[N] mole/m3

[O] mole/m3

Remark

1

0

0

0

Trivial

2

2.6616 ×10−2

1.01 ×10−6

1.3626 ×100

Physical

3

−4.4500 × 10−2

6.6934 ×10−7

−1.5134 × 100

Non-physical

4

×10−2

×10−7

−1.4975 × 100

Non-physical

4.4183

6.6676

5

−2.6805 × 10−2

−1.0025 × 10−6

3,661 × 100

Non-physical

6

−2.4590 × 10−8

−4.0000 × 100

5.0000 ×100

Non-physical

Furthermore, the equivalent equilibrium point solutions for the lower temperature give the following six finite roots tabulated in Table 5. An eigenvalue analysis was performed for both the initial and at the equilibrium states in order to estimate the time scales of the reaction at the lower temperature. For this dynamical model which is a three-ordinary nonlinear differential equations in three unknowns, namely, [NO], [N], and [O], as expected, three eigenvalues were obtained, resulting in three timescales. Let us denote the eigenvalues and their corresponding timescales as (λ1 , τ1 ), (λ2 , τ2 ), and (λ3 , τ3 ). All of these timescales evolved with time, t. At the initial state, (i) , we have

6.2 Analysis of Non-Physical Roots, T = 1,900 K

λ1i = −6.10 × 10−1 s−1 , λ2i = −4.31 × 106 s−1 and

Equilibrium point

λ3i = −1.14 × 108 s−1

3

1.28 ×108

4

×108

Notice that the dynamics of the model is linearly stable near the initial state; and the corresponding local timescales are 1 1 = 1.64 × 100 s, τ2i = = 2.32 × 10−7 s and |λ1i | |λ2i | 1 = 8.77 × 10−9 s τ3i = |λ3i |

A review of the modes of stability of the non-physical roots, 3 to 6, are considered in this section. Analyses similar to those of Eq. (59) were performed to determine the eigenvalues of the Jacobian matrix at the relevant critical points. In the lower-temperature case, we found all eigenvalues to be real. The eigenvalues, time constants, and the corresponding classification of each critical point are summarized in Tables 6, 7, and 8.

Table 6 Eigenvalues of non-physical equilibrium points, T = 1,900 K λ1 (1/s)

1.27

5

−2.53 × 107

6

1.22 ×100

λ2 (1/s)

λ 3 (1/s)

−8.14 × 10−7 −3.22 1.22

6.02 ×100

× 10−10

−2.57 × 100

× 10−6

1.83 × 101

1.17 × 107

2.24 × 108

τ1i =

The stiffness, κ1i = 1.87 × 108 , which is somewhat stiff. At the physical equilibrium state, (e) we find 7 −1

λ1e = −2.51 × 10 s

, λ2e = −1.24 × 10

−6 −1

s

and

λ3e = −1.83 × 101 s−1 Notice that the dynamics of the system is linearly stable near the equilibrium state as well, and the corresponding local timescales are 1 1 = 3.99 × 10−8 s, τ2e = = 8.10 × 105 s and |λ1e | |λ2e | 1 = 5.47 × 10−2 s = |λ3e |

τ1e = τ3e

The stiffness κ2e = 2.03 × 1013 , which is very stiff.

Table 7 Local time constants near non-physical equilibrium points, T = 1,900 K Equilibrium point

τ1 (s)

τ2 (s)

τ3 (s)

3

7.81 ×10−9

1.23 ×106

1.66 ×10−1

4

7.87

×10−9

5

3.95

×10−8

6

8.20 ×10−1

3.11

×10−10

3.89 ×10−1

8.20

× 10−5

5.47 × 10−2

8.54 × 10−8

4.46 × 10−9

Table 8 Classification of critical points, T = 1,900 K Equilibrium point

Type/sign of eigenvalues

Mode of stability

2

Sink node–negative values

Linearly stable

3

Source node–positive

Unstable

4

Saddle point–positive/negative

Unstable

5

Saddle point–positive/negative

Unstable

6

Source node–positive

Unstable

123

8300

Arab J Sci Eng (2014) 39:8287–8305 2

[NO] [N] [O]

2

1.5

3

Concentration, moles/m

Concentration, moles/m

3

[NO] [N] [O]

1

0.5

0

0

0.2

0.4

time, (s)

= 1.362 × 105

0.6

0.8

,

1

0.5

1 x 10

2

1.5

0

-6

0

0.2

0.4

0.6

time, (s)

= 2600

= 9.995 × 104

2

,

0.8

1 x 10

-6

= 1900

Fig. 3 Temporal variation in the concentrations of [NO], [N], and [O]

6.3 Effect of Reduced Pressure of −P = 1.362 × 105 N/m2 and T = 2,600 K In order to investigate the effect of pressure on the time evolution of [NO], [N], and [O], we kept the initial temperature constant at T = 2,600 K and decreased the initial concentration of each species by ten percent. The initial species concentration thus became [NO]i = [N2 ]i = [N]i = [O2 ]i = [O]i = [H]i = [OH]i = 0.9 moles/m3 . Therefore, C1 = 1.8, C2 = 2.70, C3 = 2.70, and C4 = −0.90. On solving the resulting equations numerically for P = 1.362 × 105 N/m2 and T = 2,600 K, we obtain the solution shown in Fig. 3. The computation show a convergence to final species concentration values of mole t→∞ m3 mole lim [N] = 2.9374 × 10−3 3 t→∞ m 0 mole lim [O] = 1.1329 × 10 t→∞ m3 lim [NO] = 1.9971 × 100

(69) (70) (71)

6.4 Effect of Reduced Pressure of −P = 9.95 × 104 N/m2 and T = 1,900 K Similarly, solving the resulting equations numerically for P = 9.95 × 104 N/m2 and T = 1,900 K, we obtain the solution shown in Fig. 4. The computation shows a convergence to final species concentration values of lim [NO] = 1.4211 × 100

t→∞

123

mole m3

(72)

Fig. 4 Temporal variation in the concentrations of [NO], [N], and [O]—reduced pressure

mole m3 mole lim [O] = 2.4548 × 10−1 3 t→∞ m lim [N] = 5.7895 × 10−1

t→∞

(73) (74)

7 Uncertainty Considerations In this section, we give a review of analysis of the errors and uncertainties in the present work. The vast majority of initial value problems that are of practical importance and that are often encountered in science, engineering, and other disciplines do not have closed-form (exact) analytical solutions; therefore, we are often constrained to use some suitable numerical method(s) to solve them. Numerical algorithms introduce errors, not only in each single step that is taken but also across the entire solution interval. It is evident from the numerical solutions to the Zeldovich model, that is, Eqs. (44) to (46) and Eqs. (63) to (65) that the equations describing the model are quite stiff. In order to solve for the approximate solution(s) of these stiff problems, numerical schemes which use some form of implicit discretization formula are generally employed. As a result of this, a system of nonlinear equations would have to be solved and the most reliable approach is based on the Newton’s method, which is given by tn+1 = tn −

f (tn ) f  (tn )

(75)

However, this requires the evaluation of the user-specified Jacobian matrix at each iteration. In order to circumvent this

Arab J Sci Eng (2014) 39:8287–8305

8301

time-consuming and expensive numerics, Rosenbrock [28] implemented a scheme that merged the Jacobian matrix with the numerical integration formula. This scheme resulted in the integration formula, which is generally referred to as the Rosenbrock method. The Rosenbrock method is also called linear implicit or semi-implicit Runge-Kutta method because it was derived from the fully implicit Runge-Kutta method. It is worth mentioning, however, that the Rosenbrock scheme has an inherent disadvantage in that it creates some trouble for users to supply a function for the Jacobian. Moreover, it can be quite an inconvenience working out analytical expressions for the partial derivatives necessary to evaluate the Jacobian. In view of the shortcomings of the Rosenbrock method highlighted above, Shampine and Reichelt [29] developed an improved method called a modified Rosenbrock triple on which the code ode23s is based. The code ode23s was used for all the computations reported in the present work. The complete details concerning the modified Rosenbrock triple can be found in [29]. However, a skeletal sketch of the modified Rosenbrock algorithm for advancing from step (tn , yn ) to (tn+1 = tn + h) is stated below

F1 = (tn + 0.5h, yn + 0.5k1 )

where s is the stage number of the numerical scheme. Details concerning the derivation of Eq. (77) can be found in [29]. By inspection, we can see that Eq. (77) is a quadratic interpolant and it is a continuous function. In summary, the algorithm, its error, and interpolant constitute the modified implicit Rosenbrock triple. This algorithm is of orders 3 and 2 with error control and is highly suitable for solving stiff systems. It advances from yn to yn+1 with the second-order method (that is, without local extrapolation). Because the algorithm is a W -formula, it is less sensitive to the quality of the approximate Jacobian. Furthermore, it controls the local error by taking the difference between the third- and second-order numerical solutions. With this choice and a good approximation to the Jacobian, the lowerorder algorithm is L−stable and the higher-order is A-stable. With these attributes, the error estimate is highly conservative at infinity—a highly desirable characteristic. The truncation error of the modified implicit Rosenbrock triple is order 2,  that is, ∼ Oh 2 . In order to measure the accuracy (uncertainty) of the numerical method, we designate the single-step error as the local error, ei , and the accumulated error over an interval as the global error, eg . To calculate the local error at a given step, we simply ensured that the estimated error at each integration step satisfied

k2 = W −1 (F1 − k1 ) + k1

ei ≤ max (r toli · |yi | , atoli )

F0 = F(tn , yn ) k1 = W −1 (F0 + hdT )

yn+1 = yn + hk2 , → y1 = y0 + hk2 F2 = F(tn+1, yn+1 ) k3 = W −1 [F_2 − e_32(k_2 − F_1) − 2(k_1 − F_0) + hdT ] h err or ≈ [k1 − 2k2 + k3 ] 6

(76)

where W = I − hd J, J ≈ T ≈

∂f 1 (tn , yn ) , d =  √ , ∂y 2+ 2

√ ∂f (tn , yn ) and e32 = 6 + 2 ∂t

The method uses the result y1 for propagating the solution, and if the step is a success, the F2 of the current step will be the F0 of the next one. For this reason, we will not need any additional function evaluations. This property is called first-same-as-last (FSAL), which means that the first stage of a step is the same as the last one from the end of the previous step. Because the pair samples at both ends of the step, it can recognize a quasi-discontinuity. A second-order approximation (interpolant) to y(tn + sh) is given by   s (1 − s) s (s − 2d) (77) k1 + k2 yn+s = yn + h 1 − 2d 1 − 2d

(78)

The solution from Eq. (78) is acceptable only if the weighted error is no more than the tolerance, r tol. The estimate of h that will yield an error of r tol on the next step or the next try at taking this step, as the case may be, was multiplied by a “safety factor” of 0.8 in order to avoid failures. It is worth mentioning that the absolute error tolerance can either be specified as a positive scalar or vector. A scalar tolerance applies to all components of the solution vector, while the elements of a vector of tolerances apply to corresponding components of the solution vector. In the present work, a scalar tolerance was specified, that is, atol of 1 × 10−6 . Calculating the global error, eg , involves a bit more complication as to how to combine the effects of the single-step errors. A widely used method is to calculate the Euclidean norm, which involves estimating the accumulation of error. This method is feasible when there is an exact solution to the problem been solved. However, there are no exact solutions to the model been investigated in the present study. Therefore, we used the infinity norm which is comes with the modified implicit Rosenbrock triple to estimate the global error. The algorithm is given in Eq. (79) below.   h  (k1 − 2k2 + k3 )   (79) eg = max  · 6 (max (|y| , |ynew |) + thr eshold)  where threshold = atol/r tol and h is the time step which is internally determined by the code.

123

8302

Arab J Sci Eng (2014) 39:8287–8305

The root mean square, r msg , of the global error was used to estimate the uncertainty in the global error, thus eg r msg = √ (80) N where N = 188 are the total number of steps in this study. The values of the estimate of the uncertainties in the global error of [NO], [N], and [O] for conditions investigated are tabulated in Table 17.

Table 11 Effect of pressure on eigenvalues at equilibrium point Temperature at T = 2,600 K (N/m2 )

Eigenvalues, (1/s)

λ1

λ2

λ3

1.513 × 105

−4.46 × 107

−9.23 × 10−2

−5.67 × 103

1.362 × 105

−4.28 × 107

−1.46 × 104

−1.64 × 103

Table 12 Effect of pressure on timescales at equilibrium point

8 Discussion of Results The equilibrium points and the nature of their stability were investigated at 2,600 and 1,900 K. The results are tabulated in Tables 1, 2, 3, 4, 5, 6, 7, and 8. The results show that for all the equilibrium points investigated, only one in each case was physical and the others were non-physically realizable because at least one of their equilibrium points was negative. The physically realizable equilibrium points produced eigenvalues which show that they were linearly stable as well. The non-physical equilibrium points produced eigenvalues which depict them as unstable. In order to investigate the effect of temperature and pressure on the formation of [NO], [N], and [O], a combination of four scenarios was investigated. These consisted of performing numerical calculations at two temperatures namely T = 2,600 K and T = 1,900 K and their respective pressures. Pressure was varied by keeping the relevant temperature constant while reducing all species concentration by 10 %. Since the timescales of a reaction are important descriptors of its kinetics because this indicates the times when substantial species concentration changes occur, we chose these parameters and their corresponding eigenvalues as evaluation/comparison criteria. The computed results are summarized in Figs. 1 and 2 and Tables 9, 10, 11, 12, 13, and 14. Table 9 Effect of temperature on eigenvalues at equilibrium point Temperature (K)

Eigenvalues (1/s) λ1

λ2

λ3

2,600

−4.46 × 107

−9.23 × 10−2

−5.67 × 103

1,900

−2.51 × 107

−1.24 × 10−6

−1.83 × 101

Table 10 Effect of temperature on timescales at equilibrium point Temperature Time scales (s)

Stiffness, κ

2,600

τ1 = 1/ |λ1 | τ2 = 1/ |λ2 | τ3 = 1/ |λ3 | τmax /τmin 2.24 × 10−8 1.08 × 100 1.76 × 10−4 4.84 × 108

1,900

3.99 × 10−8 8.10 × 105

(K )

123

5.47 × 10−2 2.03 × 1013

Stiffness, κ

Temp at Time scales (s) T = 2,600 K (N/m2 )

τ1 = 1/ |λ1 | τ2 = 1/ |λ2 | τ3 = 1/ |λ3 | τmax /τmin 1.513 × 105 1.362

× 105

2.24 × 10−8 1.08 × 100

1.76 × 10−4 4.84 × 108

2.33 × 10−8

6.08 × 10−4 2.61 × 104

6.84 × 10−5

Table 13 Effect of pressure on eigenvalues at equilibrium point Temperature at T = 1,900 K (N/m2 )

Eigenvalues, (1/s)

λ1

λ2

λ3

1.106 × 105

−2.51 × 107

−1.24 × 10−6

−1.83 × 101

× 104

−2.26 × 107

−3.05 × 10−6

−5.96 × 100

9.952

Table 14 Effect of pressure on timescales at equilibrium point Temp at Time scales (s) T = 2,600 K (N/m2 )

Stiffness (κ)

τ1 = 1/ |λ1 | τ2 = 1/ |λ2 | τ3 = 1/ |λ3 | τmax /τmin 1.106 × 105

3.99 × 10−8 8.10 × 105

5.47 × 10−2 2.03 × 1013

9.952 × 104

4.42 × 10−8 3.28 × 105

1.68 × 10−1 7.42 × 1012

Figures 1 and 2 show the rate of production of [NO], [N], and [O], as a function of time at temperatures of T = 2,600 K and T = 1,900 K while maintaining the pressures constant at 1.513 × 105 and 1.106 × 105 N/m2 , respectively. Figure 1 shows a rapid rate of production of [NO] and very fast consumption of the [N] atoms as evidenced from sharp rate of descent of the [N] atoms. The rate of production of the oxygen atoms was moderate. The rate of these reactions was extremely very fast as it lasted only from about 10−9 to 10−7 s, thereafter all the species equilibrated to their respective equilibrium values. Figure 2 presents the time history of the rate of production of [NO], [N], and O, at the lower temperature of 1,900 K. The figure shows the rapid rate of production of [NO] and the sharp rate of decline of both the

Arab J Sci Eng (2014) 39:8287–8305

[N] and [O] atoms. The reason for the sharp decline in the calculated [O] atoms could be as a result of the fact that the temperature was a bit too low to produce enough [O] atoms to offset the rate at which it was reacting with the [N] atoms to produce [NO]. The maximum rate of production of [NO] and consumption of [N] occurred at about 10−7 s. It is interesting to note that the [N] atoms recovered somewhat at about 3 × 10−7 s after all the species relaxed to their steady-state values. The maximum value of [NO] produced was higher under the high-temperature regime. The calculated values of [NO] at these temperatures do in fact support experimental observation where it has been shown that the rate of [NO] production is highly temperature-dependent. Figures 3 and 4 show the rate of production of [N O], [N], and [O] as a function of time at pressures of 1.362 × 105 and 9.95 × 104 N/m2 while maintaining the temperatures constant at T = 2,600 K and T = 1,900 K, respectively. Figures 3 and 4 are both qualitatively and quantitatively similar to Figs. 1 and 2. The inference that could be drawn from this is that the effect of pressure on the rate of production of [NO] and consumption of [N] and [O] is at best marginal–Le Chatelier’s principle. As stated above, the results of the parametric studies are presented in Tables 9, 10, 11, 12, 13, and 14. The tables show that the reaction timescales were consistently faster at 2,600 K than at 1,900 K. This is expected since the rate of production of [NO] is highly dependent on temperature. The stiffness of the system was found to be lower at 2,600 than at 1,900 K. This could be attributed to the fact that at high temperatures, the dynamics of the system are uniform and homogeneous. The results also show that the values of the timescale and stiffness were superior at the higher pressure, see Tables 10, 12 and 14. 8.1 Phenomenological interpretation of Figs. 1 and 2 The Zeldovich reaction mechanism postulates that N O formation starts only if atomic oxygen, [O], is available. Only the oxygen atom is able to split up the strongly bound N2 molecule at temperatures below 3,000 K. The thermal decomposition of N2 to N begins at temperatures considerably above 3,000 K. With this preamble stated, we now describe the physical and chemical phenomenon that resulted in the shapes of Figs. 1, 2, 3, and 4. The sharp rise of [NO] production observed in Figs. 1, 2, 3, 4 is due to the fact that [NO] increases exponentially with temperature and is linearly dependent on residence time during the combustion phase. After the end of combustion, the temperature and pressure drop; consequently, the [N O] concentration predicted using the Zeldovich mechanism assumptions drops steeply, Figs. 2 and 4. However, the [N O] concentration soon reaches steady state

8303

from 3.0 × 10−7 s after the end of combustion reaction. In Figs. 1 and 3, the steady (frozen) state is attained at a much earlier time (1.0 × 10−7 s). This is because of the much higher temperature. The oxygen atoms were also found not to decrease instead there was a slight rise before leveling-off. In order to understand this behavior, a closer examination of k1 f in Eq. (20) is necessary. The drop in temperature and pressure at the end of combustion leads to a sharp drop in kif . This is expected since kif depends strongly on temperature Ea )). This leads to overall drop in N atom concen(∼ exp( RT tration after the end of combustion reaction. Moreover, from 3.0 × 10−7 s onwards, the formation of [N O] remains practically frozen and it then remains constant thereafter. The temporal variation in nitrogen concentration obtained from the solution of Eq. (45) is shown in Figs. 1, 2, 3, and 4. The reaction rate constants were computed using instantaneous temperatures specified. The extended Zeldovich mechanism predicts a sharp drop in the total N atoms due to recombination process with oxygen atoms. The net production of N atoms is determined by the source term based on the reaction mechanism showed in Eqs. (20)–(22). Numerical evaluation of the source terms for the formation of N atoms shows that the rate of depletion of N atoms is strongly controlled by the forward reaction rate constant, k2 f for reaction 2. Since reaction 2 has a relatively small activation energy (as compared to the forward reaction rate constant of reaction 1), there is a rapid depletion of N atoms at low temperatures. After the start of combustion reaction, the temperature of the unburned gas is still low, and hence, there is no appreciable increase in the formation of N atoms in this zone. The burnt zone, however, has a high temperature, which supports the spontaneous formation of N atoms as the burned volume increases. The rapid increase in the N atom specie is enhanced by the increasing temperature and the presence of higher amounts of heated N2 molecules. On completion of combustion, the temperature and pressure starts to drop and the overall rate of the formation of N atoms ceases and remains frozen thereafter. The temporal variation in the concentration of [O] at 1,900 K is shown in Figs. 2 and 4. It can be seen that there is a sharp decrease in [O] concentration during combustion, that is, from the start of combustion to the end of combustion followed by a gradual rise to a constant leveling-off (frozen state) value. The sharp decrease is due to the reaction between the O and N atoms, which leads to the formation of [NO]. The rise observed in O is due to the copious O atoms which were produced as a result of high combustion temperature. In Figs. 1 and 3, we see that there is a slight increase in the [O] concentration between t = 0 and t = 1.0 × 10−7 s. This is due to the excessive production of O as a result of the thermal decomposition of O2 to O during the combustion phase. Moreover, the dissociation equilibrium moves to higher O

123

8304

Arab J Sci Eng (2014) 39:8287–8305

concentrations at higher temperatures leading to more N O formation. In summary, because of the physical and chemical phenomenon described above, NOx reduction techniques are based on combustion modification modes, which have the effect of reducing the peak flame temperature, the residence time in the flame zone, and the excess oxygen concentration. From the foregoing, we see that the Zeldovich mechanism did in fact accurately describe the physically expected behavior of [N O], [N], and [O] species formation from the start of combustion to the end.

Table 16 Computation statistics: P = 1.513 × 105 N/m2 , T = 2,600 K—number of successful steps, failed attempts, and solutions of linear equations with relative and absolute error tolerances r tol of 1.0 × 10−6 and atol of 1.0 × 10−6 , respectively

8.2 Comparison with Referenced Related Work

Temperature (K)

A new reaction scheme for NOx production was developed by Miller et al. [30] and incorporated into the General Engine SIMulation (GESIM) engine combustion simulation model. The new reaction scheme is a super-extended Zeldovich mechanism (SEZM), which was found to predict NOx formation levels to within 10 % of engine test data for several engine test conditions, whereas the 3-reaction, extended Zeldovich mechanism (EZM) was shown to have errors of up to approximately 50 % or more for similar conditions. However, the global error result presented in this study for NOx production varies from (±4.10 to ±4.34 %). Moreover, in the present work, our results were both qualitatively and quantitatively comparable to those presented in Figs.4 to 9 of [6,11], Figs. 2 to 8 of [12], and Fig. 14 of [31]. 8.3 Uncertainties (Global Error) Since actual (global) errors may exceed the local tolerances, the values of the relative and absolute error tolerances, r tol and atol, respectively, were conservatively chosen so as to ensure that global errors do not exceed local tolerances. The results of the computation statistics are tabulated in Tables 15 and 16. It can be seen from the tables that for the conditions investigated, there were no failed attempts. This goes to show that the relative and absolute error tolerances of r tol of 1 × 10−6 and atol of 1 × 10−6 , respectively, were indeed chosen prudently. Table 17 shows the estimate of the uncertainties in the global errors of [NO], [N]n and [O] for the conditions invesTable 15 Computation statistics: P = 1.106 × 105 N/m2 , T = 1,900 K—number of successful steps, failed attempts, and solutions of linear equations with relative and absolute error tolerances r tol of 1.0 × 10−6 and atol of 1.0 × 10−6 , respectively P = 1.513 × 105 , T = 1,900 K

Computation statistics

Successful steps

Failed attempts

Solution of linear equations

187

0

561

123

P = 1.513 × 105 , T = 2,600 K

Computation statistics

Successful steps

Failed attempts

Solution of linear equations

127

0

381

Table 17 Estimate of the uncertainties in the global error of [N O], [N] and [O] Species concentrations mole m3 [NO] (%)

[N] (%)

[O] (%)

2,600

±4.1 0

±10.1 0

± 5.73

1,900

±4.34

± 9.33

±9.07

tigated in this study. Table 17 shows that the estimate of the uncertainties in the global errors of [NO], [N]n and [O] vary between (±4.10 and ±4.34 %), (±9.33 and ±10.1 %)n and (±5.73 and ±9.07 %)n respectively, for the conditions investigated. We can see from the table that at high temperature, the variation (uncertainty) in [NO] specie was less than those of [N] and [O]. This could be attributed to the fact that higher temperatures favor the uniform production of [NO]. Moreover, the vigorous breaking down and rearrangement of both the nitrogen and oxygen electronic bonds at high temperatures could be responsible for their higher uncertainties. Similar trend was observed at the lower temperature. However, the variability of [NO] specie was lower. This could be due to the fact that the reactions were slower and hence the increased non-uniformity (uncertainty) observable in both the [NO] and [O] species. 9 Conclusions A computational analysis of the extended Zeldovich mechanism was undertaken. The effect of temperature and pressure was investigated. The results of our investigation show that •

•

•

For all the equilibrium points investigated, only one in each case was physical and the others were nonphysically realizable because at least one of their equilibrium points was negative. The physically realizable equilibrium points produced eigenvalues, which show that they were linearly stable as well. The non-physical equilibrium points produced eigenvalues, which depict them as unstable. High temperatures result in faster rate of production of both [NO] and [O] and rapid rate of depletion [N].

Arab J Sci Eng (2014) 39:8287–8305

• • •

• •

•

At moderately low temperatures (∼1,900 K), the rate of depletion of both oxygen and nitrogen atoms was very high. It was observed that pressure has a minor or at best a marginal influence on the production of [NO], [N], and [O]. The results of our study show that the effect of lowering the initial species concentrations while leaving temperature constant lowers the pressure proportionately, with the concomitant slowing down of the collision time, as well as the fastest and slowest time scales. It was found that reducing or lowering pressure does not significantly affect the dynamics of the system, for example, stiffness. Since actual (global) errors may exceed the local tolerances, the values of r tol and atol stated above were conservatively chosen so as to ensure that global errors do not exceed local tolerances. The estimate of the uncertainties in the global errors of [NO], [N], and [O] for the conditions investigated were found to vary between (±4.10 and ±4.34 %), (±9.33 and ±10.1 %), and (±5.73 and ±9.07 %), respectively.

Acknowledgments The work reported herein was graciously supported by the Nigerian Defence Academy, Kaduna, Nigeria, the Texas Southern University, Houston, Texas, USA, and the Kwara State University, Ilorin, Nigeria. The authors wish to express their sincere thanks to the referees for their insightful and constructive remarks, which led to improving this paper considerably.

References 1. Anetor, L.: Experimental and Numerical Simulation of Charge Motion in Internal Combustion Engine, PhD Thesis, University of British Columbia, Vancouver, Canada, (1994). https:// circle.ubc.ca/bitstream/handle/2429/6995/ubc_1994-953047.pdf. (Accessed August 9, 2014) 2. Torres, D.J.; Trujillo, M.F.: KIVA-4: an unstructured ALE code for compressible gas flow with sprays. J. Comput. Phys. 219, 943– 975 (2006) 3. Singh, S.; Reitz, R.D.; Musculus, M.P.B.: Comparison of the characteristic time (CTC), representative interactive flamelet (RIF), and direct integration with detailed chemistry models against optical diagnostic data for multi-mode combustion in a heavy-duty DI diesel engine. SAE paper 2006-01-0055 (2006) 4. Hildenbrand, F.; Schulz, C.; Wolfrum, J.; Keller, F.; Wagner, E.: Laser diagnostic analysis of NO formation in a direct injection diesel engine with pump-line-nozzle and common rail injection systems. Proc. Combust. Inst. 28(1), 1137–1143 (2000) 5. De Zilwa, S.; Steeper, R.: Predicting NOX Emissions from HCCI Engines Using LIF Imaging. SAE Technical Paper 2006-01-0025 (2006) 6. Aithal, S.M.: Modeling of NOx formation in diesel engines using finite-rate chemical kinetics. Appl. Energy 87(7), 2256–2265 (2010) 7. Andersson, M.; Johansson, B.; Hultqvist, A.; Nhre, C.: A real time NOx model for conventional and partially premixed diesel combustion. SAE 2006-01-0195 (2006) 8. Eggels, R.L.G.M.: Modelling of combustion processes and NO formation with reduced reaction mechanisms. Ph.D. dissertation, Eindhoven. University of Technology (1995)

8305 9. Zheng, M.; Reader, G.T.; Hawley, J.G.: Diesel engine exhaust gas recirculation review on advanced and novel concepts. Energy Convers. Manag. 45, 883900 (2004) 10. Biruduganti, M.; Gupta, S.; Bihari, B.; McConnell, S.; Sekar, R.: Air separation membranes an alternative to EGR in large bore natural gas engines. ASME Internal Combustion Engine 2009 Spring technical conference paper ICES2009-76054 (2009) 11. Aithal, S.M.: A comparative study of NOx computation methods coupled to quasi-dimensional models in SI engines. In: Proceedings of the ASME 2012 Internal Combustion Engine Division Fall Technical Conference ICEF2012, Vancouver, BC, Canada, Sept 23–26 (2012) 12. Aithal, S.M.: Equilibrium chemistry calculations for assessment of NOx abatement strategies in internal combustion engines. In: Proceedings of the ASME 2011 Internal Combustion Engine Division Fall Technical Conference, ICEF2011-60030, Morgantown, West Virginia, USA Oct 2–5 (2011) 13. Benson, S.W.: The Foundations of Chemical Kinetics. McGrawHill, New York (1960) 14. Weston, R.E.; Schwartz, H.A.: Chemical Kinetics. Prentice-Hall, Englewood Cliffs, NJ (1972) 15. Smith, I.W.M.: Kinetics and Dynamics of Elementary Reactions. Butterworth, London (1980) 16. Laidler, K.J.: Theories of Chemical Reaction Rates. McGrawHill, New York (1969) 17. Law, C.K.: Combustion Physics. Cambridge University Press, New York (2006) 18. Glassman, I.; Yetter, R.A.: Combustion. Academic Press, Burlington, MA (2008) 19. Zeldovich, Y.B.: The oxidation of nitrogen in combustion and explosives. Acta Physicochimica USSR 21, 577 (1946) 20. Fenimore, C.P.: Reactions of fuel-nitrogen in rich flame gases. Combust. Flame 26, 249–256 (1976) 21. Malte, P.C.; Pratt, D.T.: The role of energy-releasing kinetics in NOx formation: fuel- lean, jet-stirred CO-air combustion. Combust. Sci. Technol. 9, 221–231 (1974) 22. Dean, A.M.; Bozzelli, J.W.: Combustion Chemistry of Nitrogen in Mechanical Engineering Handbook. pp. 125–343. CRC Press LLC, Frank Kreith, Boca Raton (1999) 23. Kuo, K.K.: Principles of Combustion. 2nd edn. WileyInterscience, (2005) 24. Lavoie, G.A.; Heywood, J.B.; Keck, J.C.: Experimental and theoretical study of nitric oxide formation in internal combustion engines. Combust. Sci. Orid Technol. 1, 313 (1970) 25. JANAF Thermochemical Tables 3rd Edn. Thermal Group, Dow Chemical USA, Midland, MI (1985) 26. Roache, P.: Verification and Validation of Computational Fluid Dynamics Simulations. Hermosa, Albuquerque, NM (1998) 27. Roache, P.: Quantification of uncertainty in computational fluid dynamics. Ann. Rev. Fluid Mech. 29, 123–160 (1997) 28. Rosenbrock, H.H.: Some general implicit processes for the numerical solution of differential equations. Comput. J. 5, 329–330 (1963) 29. Shampine, L.F.; Reichelt, M.W.: The Matlab ODE suite. SIAM J. Sci. Comput. 18(1), 122 (1997) 30. Miller, R.; Davis, G.; Lavoie, G.; Newman, C.; Gardner, T.; SAE Technical Paper 980781. A Super-Extended Zeldovich Mechanism for NOx Modeling and Engine Calibration. (1998) doi:10.4271/ 980781 31. Ericson, C.; Westerberg, B.; Andersson, M.; Egnell, R.: Modelling diesel engine combustion and NOx formation for model based control and simulation of engine and exhaust aftertreatment systems. SAE Technical Paper 2006-01-0687 (2006). doi:10.4271/ 2006-01-0687

123