Algorithmic Analysis of Chemical Dynamics of the Autoignition of NH3–H2O2/Air Mixtures

The dynamics of a homogeneous adiabatic autoignition of an ammonia/air mixture at constant volume was studied, using the algorithmic tools of Computational Singular Perturbation. Since ammonia combustion is characterized by both unrealistically long ignition delays and elevated NO x emissions, the time frame of action of the modes that are responsible for ignition was analyzed by calculating the developing time scales throughout the process and by studying their possible relation to NO x emissions. The reactions that support or oppose the explosive time scale were identified, along with the variables that are related the most to the dynamics that drive the system to an explosion. It is shown that reaction H 2 O 2 (+M) → OH + OH (+M) is the one contributing the most to the time scale that characterizes ignition and that its reactant H 2 O 2 is the species related the most to this time scale. These findings suggested that addition of H 2 O 2 in the initial mixture will influence strongly the evolution of the process. It was shown that ignition of pure ammonia advanced as a slow thermal explosion with very limited chemical runaway. The ignition delay could be reduced by more than two orders of magnitude through H 2 O 2 addition, which causes only a minor increase in NO x emissions.

Ammonia is considered a source of hydrogen that can be easily manufactured, liquified, and economically stored safely at relatively low pressure, providing energy per unit volume greater than that of pure hydrogen [1]. When compared to hydrocarbons and hydrogen, ammonia is the least expensive onboard fuel (in terms of cost per gigajoule) with the cheapest 100 km driving range, and it has narrower flammability limits and lower risk of explosion [2]. Another attractive feature is its high octane rating leading to higher compression ratios inside internal combustion (IC) engines and thus to higher efficiency [1]. Furthermore, ammonia can be easily made from hydrogen, or produced from water electrolysis with renewable energy and atmospheric nitrogen through the Haber process [3].
There are several challenges in realizing ammonia as a fuel, such as long ignition delay [4], low flame speed and temperature, narrow flammability range, and high ignition energy and autoignition temperature [1]. Application of ammonia to fuel engines goes back to World War II, when the shortage in fossil fuel led some European nations to use ammonia as a fuel in buses [5]. Ammonia regained its interest as an alternative fuel in the past decade, during which plenty of studies have been reported on ammonia usage in compression ignition (CI) and spark ignition (SI) engines [3,[5][6][7][8]. Several studies have identified strategies and additives that would assist ammonia in reaching the optimal combustion behavior for its application in combustion engines. Hydrogen, with its high flammability range and flame speed, can be considered as one of the potential carbon-free additives to be used. Frigo and Gentili [3] studied the effect of introducing hydrogen onboard by catalytic reforming of ammonia to a 4-stroke twin-cylinder SI engine. By producing hydrogen from ammonia and introducing it in the fuel mixture, the minimum ignition energy dropped from 8 mJ to 0.018 mJ, with an increase in the laminar flame velocity from 0.015 ms −1 to 3.51 ms −1 was reported [9]. The results showed that with a minimum content of hydrogen, which differed as a function of engine load, ammonia combustion could be accelerated to achieve proper engine behavior [3].
Li et al. [10] conducted an experimental study to identify the combustion characteristics and emissions of the hydrogen-ammonia-air mixture. The flame velocity was measured using a Bunsen burner. The authors observed that the flame velocity of the mixture increased as the hydrogen mole fraction increased. Controlling the concentration of ammonia in the mixture gave the flexibility to control the flame velocity, thus, widening the possible applications of such fuel. The maximum flame velocity of the mixture occured at φ = 1.01, which was similar to the ammonia-air mixture, from which it can be concluded that hydrogen accelerates the combustion while ammonia has the major effect on the maximum flame speed. NO x emission increased with increasing equivalence ratio.
In [4], an experiment was reported in a dual-fuel engine with diesel being injected along with ammonia to ignite the mixture. Since the stoichiometric air-fuel ratio of ammonia is more than half that of diesel, additional ammonia can be introduced in the combustion cylinder to compensate for its lower energy density [7]. The results indicated that a minimum amount of diesel was needed in order to run the engine with the same output. The amount of NO x released decreased with the increase in ammonia. This happened because the cylinder temperature decreased with higher ammonia content, due to the low flame temperature of ammonia [4,8]. For the same reason, the amount of unburnt hydrocarbons (UHC) increased due to incomplete combustion of diesel [8], except for low diesel ratios [4]. Carbon dioxide emission dropped drastically without any effect on the output power of the engine when 40% to 60% of the power delivered to the engine was replaced by ammonia. By substituting only 3% of intake air with the gaseous fuel, 10% less diesel was used, which significantly reduced the CO 2 emissions [8]. "Dissociated" ammonia, i.e., ammonia mixed with hydrogen derived from its catalytic processing, was shown to have reduced ammonia and NOx emissions compared to pure ammonia [8]. In [11], CI engine operation with a mixture of ammonia with dimethyl ether (DME) was studied and it was shown that ammonia addition limited the operation range. This was explained by the high latent heat of ammonia and lower heating value that caused the cylinder temperature to drop by nearly 100 o C.
Li et al. [12] numerically studied ammonia laminar burning velocity at various contents of hydrogen using CHEMKIN [13]. The study also looked into the ignition delay using a homogeneous charge compression ignition model of CHEMKIN with realistic engine conditions. The simulation results proved that ammonia burning velocity increased with the mole fraction of blended hydrogen. This enhancement was caused by chemical and transport effects [12]. Using CHEMKIN, the same author studied the effect of oxygen enrichment on the combustion characteristics of ammonia [14]. As of the short-lived H, O, and OH radicals, they tended to increase with oxygen-enriched combustion [12,14]. HNO, which is a significant predecessor for NO formation, also increased at lean conditions, which lead to the increase in NO. At rich conditions, NH 2 , NH, and N radicals reduced the reaction rate of NO decreasing its mole fraction [14].
All this rich phenomenology points to the need for a rigorous exploration of the pertinent chemical dynamics that seems to be mostly missing. A notable exception is Ref. [15], where the autoignition of NH 3 /H 2 /O 2 mixtures were studied and branching reactions were identified, to which the results were most sensitive. Several experimental studies point to the crucial importance of additives, which have so far been determined empirically, rather than algorithmically. Of particular importance is also the relatively slow dynamics that determines NO x formation. In this work we have performed a novel study of the autoignition dynamics of NH 3 /air mixture using the method of Computational Singular Perturbation (CSP). The focus is on the determination of the reactions that contribute significantly to the generation of the time scale that characterizes ignition at engine relevant conditions and on the algorithmic determination of an additive that would reduce ignition delay, possibly without increasing NO x emissions.

The CSP Methodology
The homogeneous adiabatic autoignition of ammonia/air mixtures at constant volume are considered here. A chemical kinetics mechanism consisting of N = 34 species, E = 5 elements (H, N, O, Ar, and He), and K = 229 elementary reactions was employed. This mechanism was constructed by removing all carbon chemistry from the Aramco 2.0 mechanism [16,17] and by adding NH 3 -chemistry [18]. In what follows, for the K elementary reactions, each direction will be considered as a distinct unidirectional reaction. The system of the species mass fractions and temperature governing equations may be considered as of the form: where y is the state column vector of the species mass fraction of dimension N, ρ is the density of the mixture, T is the temperature, c v is the heat capacity, R is the universal gas constant, W is an N × N diagonal matrix with the species molecular weights, h c is the N-dimensional vector of the species absolute enthalpies and U = [1, 1, . . . , 1] [19,20]. Furthermore, S k denotes the stoichiometric vector and R k the reaction rate of the k-th unidirectional reaction. The species and energy governing Equations (1) and (2) can be considered as a (N + 1)-dimensional system of the form: where z = [y, T] T is a (N + 1)-dim. vector that consists of the vector of the species mass fractions y and the temperature T, while g(z) is the (N + 1)-dim. column vector field. Since the chemical source term g(z) represents K reversible reactions, whose two directions are treated separately, then: whereŜ k is the generalized stoichiometric vector [21,22]. Equation (3) can be further cast into a new from by resolving the vector field g(z) along the CSP basis vectors, so that: where the (N + 1)-dimensional column vector a n is the CSP basis vector of the n-th mode and b n is the (N + 1)-dimensional row n-th dual basis vector (b i · a j = δ i j ) [23][24][25]; f n is the related amplitude and is always considered positive, by appropriate altering the sign of b n (and of a n in order to preserve orthogonality). Each CSP mode, say the n-th (a n f n ), in Equation (5) is characterized by • its time scale, τ n , which provides a measure of the time frame of its action • its amplitude, f n , which provides a measure of its contribution to the evolution of the system • the variables that pertain to this mode.
When the M fastest time scales in the dynamics of the system in Equation (5) are of dissipative character and much faster than the rest (τ 1 < · · · < τ M ), the system can be reduced to the following two relations: The first relation is an M-dim. system of algebraic equations, which defines an N + 1 − M-dim. surface in the (N + 1)-dim. tangent space. This surface approximates the slow invariant manifold (SIM). On this surface, the solution evolves. The second relation in Equation (6) is an (N + 1)-dim. system of ordinary differential equations that describes the way the slow system evolves along the SIM. This system does not contain fast time scales and its dynamics is characterized by the slow ones. These approximations are valid in the case when the trajectory is not close to the edges of the SIM [26,27].
The fastest time scales in chemical processes originate from chemical kinetics [21,[28][29][30][31][32][33][34]. These time scales are approximated by the relation τ n = |λ n | −1 (n = 1, . . . N + 1 − E), where λ n is the n-th nonzero eigenvalue of the Jacobian J of the vector field g(z) in Equation (3). A positive (negative) real part of λ n is generated by components of the vector field g(z) that tend to drive the system away from (towards) equilibrium and, thus, τ n is called explosive (dissipative). Denoting the n-th right (column) and left (row) eigenvectors of J as α n and β n , the respective eigenvalue is defined as λ n = β n · J · α n . If the n-th eigenvalue λ n is real and nonzero (for complex pairs see [35]), then it can be expressed as a summation of 2K terms: where grad = (∂/∂z 1 , . . . , ∂/∂z N+1 ) and according to Equation (4), [22,36]. From the expression in Equation (7) the Time scale Participation Index (TPI) can be introduced as: where n = 1, . . . , (N − E + 1), k = 1, . . . , 2K [36]. J n k provides a measure of the relative contribution of the k-reaction to the eigenvalue λ n and, therefore, to the time scale τ n . Positive J n k indicates that the k-th reaction reinforces the explosive character of τ n . On the other hand, negative J n k indicates that the k-th reaction reinforces the dissipative character of τ n . The use of TPI has been already used in a large number of reacting flow, biological and pharmacokinetics problems in order to identify the driving processes [22,[33][34][35][35][36][37][38][39][40][41][42][43][44][45][46][47].
Given the expression the chemical kinetics term g in Equation (4), the amplitude of the n-th CSP mode f i (i = 1, . . . , N + 1), defined in Equation (5), can be expressed as the sum of 2K additive terms: where d n k = (b n ·Ŝ k )R k (k=1,2K) is the relative contribution of the k-th reaction to the amplitude f n . The amplitudes f N−E+2 to f N+1 equal zero, due to the conservation of the E elements. The relative contribution of each reaction to the amplitude f n can be evaluated by the Amplitude Participation Index (API): where n = 1, . . . , N + 1 and k = 1, . . . , 2K [21,22,24,33,[48][49][50][51]. Since by definition f n > 0, positive (negative) values of P n k indicate a contribution of the k-th process towards strengthening (weakening) the impact of the n-th mode. In the case of an exhausted mode f m ≈ 0 (m=1,M) (see Equation (6)), Therefore, a relatively large value of P n k indicates a large contribution from the k-th process to (i) the cancellations between the various terms in the expression in Equation (9) and (ii) the m-th component of the SIM, the development of which is characterized by the m-th fast time scale, τ m .
The variables that relate the most to the m-th fast time scale are identified by the CSP Pointer (Po) for the m-th mode: where due to orthogonality a 1

Autoignition of Ammonia
The dynamics of homogeneous adiabatic autoignition of a lean NH 3 /air mixture at constant volume was analyzed for the initial conditions of T 0 = 800 K, p 0 = 20 atm, and a stoichiometric ratio of φ= 0.7. This is the reference case and will be analyzed first. For future reference, the reactions influencing the most the dynamics of the system in Equation (3) are summarized in Table 1. Unlike hydrocarbon fuels, where ignition takes place in few milliseconds, the ignition delay time, which is defined as the period from t = 0 to the point in time where temperature gradient is maximum, was calculated to be t ign ≈ 2.4 s, indicating that ammonia has high resistance to ignition and a long ignition delay. As a result, it is difficult to time ammonia ignition realistically, making the operation of the engine practically impossible.
The evolution of the developing dissipative (grey) and explosive (red) time scales is depicted in Figure 1, along with the profile of the temperature (blue). The time scales of the system extend over a wide range of values, varying from 10 −12 to 10 10 s. Among them, two explosive time scales exist throughout the process, one fast and one slow, which are denoted as as τ e, f and τ e,s , respectively. These two time scales are present until shortly before the end of the steep temperature rise. This point marks the completion of the ignition delay. The period extending from t = 0 to the point where the two explosive time scales meet is known as the explosive stage [22,41]. Throughout this period, the fast explosive timescale characterizes the prevailing dynamics [22], since it is much faster than the slow one.  Figure 1 indicates that τ e, f initially decelerates until ∼0.75 s and then accelerates slowly until ∼2.3 s, when the temperature starts to increase much faster than previously. Furthermore, it is illustrated that the period, in which τ e, f undergoes a slow acceleration and its magnitude does not vary much, occupies the largest portion of the explosive stage. After ∼2.3 s, τ e, f starts a rapid acceleration until it reaches its minimum value. From that point, τ e, f decelerates rapidly until it meets the slow time scale τ e,s , which is also rapidly accelerating; they both disappear when they meet. Note that the point where the two time scales meet is right before the point where the temperature completes its steep increase. This feature appears also in the autoignition of H 2 /air, CH 4 /air, n-hexane/air, n-heptane/air, DME/air, and ethanol/air mixtures [22,34,37,38,41,43,44,46,[52][53][54][55][56].
In order to get insight of the dynamics that influences the autoignition process and is responsible for the long ignition delay, a set of diagnostics was obtained in five indicative points throughout the process along the explosive stage. The points are depicted in Figure 2 and are indicative of significant features the τ e, f . Point P 1 was chosen to represent the early period, where τ e, f exhibits a deceleration. Point P 2 was chosen in the second period, where the magnitude of the slowly accelerating τ e, f does not change significantly. Finally, points P 3 to P 5 where chosen in the final period, where τ e, f exhibits a steep acceleration and, in particular, P 5 was chosen as the point where τ e, f obtains its minimum value. The reactions that either support or oppose (i) the explosive action and (ii) the impact of the explosive mode are identified by the CSP Timescale Participation Index (TPI) and Amplitude Participation Index (API) and are displayed in Table 2. This Table reports also values of the Pointer (Po), which identifies the variables that affect the explosive mode the most. Only the largest contributions (larger than 5%) are included in the Table. (a) (b)   Figure 2); numbers in parenthesis denote powers of ten. Only contributions larger than 5% are displayed. According to Table 2, at the start of the process, the chain initiating reaction 18f: H 2 O 2 (+M) → OH + OH (+M) is the largest contributor to the τ e, f and it is the one that sets the time frame of action of the explosive mode. Furthermore, reactions 110f: NH 3 + NH 2 → N 2 H 3 + H 2 and 71f: NH 2 + HO 2 → H 2 NO + OH are shown to contribute to a lesser degree to τ e, f and, along with reaction 18f support its explosive character. The largest opposition to τ e, f comes from the chain terminating reaction 64f: NH 2 + NO → N 2 + H 2 O. According to API, the same set of reactions are found to provide the largest contribution to the amplitude of the explosive mode. This implies that the reactions which support (oppose) the explosive character of the explosive mode tend to increase (decrease) its impact.
As the process evolves in the period where τ e, f exhibits a mild acceleration, reaction 156b: N 2 H 4 + OH ← NH 3 + H 2 NO, which initially had a smaller relative contribution, becomes the dominant contributor to both the timescale and amplitude of the mode, while the relative contribution of reaction 18f decreases slightly. These are the reactions that provide the system with OH radicals, thus contributing to the most reactive component of the radical pool, which is needed in the final part of the explosive stage [41]. The main opposition to τ e, f now originates from reaction 208f: H 2 NO + NH 2 → HNO + NH 3 , which converts NH 2 back to NH 3 by depleting the reactant of 156b that promotes ignition. Throughout the slow acceleration of τ e, f (up to P 3 ), the reactions that influence the explosive dynamics remain the same; i.e., 18f and 156b being the major supporters of the explosive mode and 208f being the major opponent. Furthermore, reaction 70b: NH 3 + HO 2 ← NH 2 + H 2 O 2 appears to become an important contributor opposing both the explosive character of τ e, f and the magnitude of its amplitude f e, f . This reaction is competing with reaction 18f that promotes ignition for the consumption of H 2 O 2 , which is shown by the CSP pointer to be the most important species in this period of the autoignition process. While 18f uses H 2 O 2 in order to produce OH that contributes to the radical pool, 70b converts it to the relatively long-lived HO 2 and NH 3 .
Towards the end of the explosive stage, during the steep acceleration of τ e, f , hydrogen-related reactions become important. The major relative contribution to τ e, f becomes the chain branching reaction 5f: O 2 + H → O + OH, which is strongly exothermic and continues the OH radical production. Reactions 63f and 64f are found to contribute substantially at P 4 , however they compete with each other. Reaction 63f produces the radicals NNH and OH (thus promoting ignition), while reaction 64f produces the stable molecules N 2 and H 2 O (thus opposing ignition). As a result of this competition, the relative contribution of reaction 3f: H 2 + OH → H + H 2 O becomes dominant. This reaction is strongly exothermic, thus promoting the explosive character of τ e, f and its impact by increasing f e, f .
The CSP Pointer (Po) identified H 2 O 2 , which is a reactant of the ignition promoting reactions 18f and 156b, to be most related species to the explosive mode throught the explosive stage, except close to its end, where the steep rise of temperature and steep acceleration of τ e, f are manifested. After t = ∼0.75 s , temperature is identified by Po to be the variable mostly related to the explosive dynamics. The stage where temperature is not strongly pointed by CSP is defined as the chemical runaway and emerges in the initial portion of the explosive stage [22,29,42]. On the other hand, when the temperature is pointed strongly, then the system undergoes the thermal runaway, which usually develops in the last part of the explosive stage [22,29,42]. In previous studies relating to the autoignition of H 2 /air, CH 4 /air, n-hexane/air, n-heptane/air, DME/air, and ethanol/air mixtures, it was shown that chemical runaway regime practically occupies all the explosive stage, while the thermal runaway regime only appears at the very final portion of the explosive stage [22,34,37,38,41,43,44,46,[52][53][54][55][56]. The results diaplayed in Table 2 suggest that autoignition of ammonia is qualitatively different than this of hydrocarbons in that the chemical runaway process is very short and the process develops as a thermal explosion.
The use of CSP Pointer, as an indicator of probable additives in the starting fuel/air mixture from the set of intermediate species, has been successfully employed in previous works [42,43,45,57]. This tool was shown most successful when employed during the chemical runaway regime, where the explosive activity is chemically driven. The results in Table 2 indicate that H 2 O 2 is pointed throughout the largest part of the explosive stage, which means that is strongly related to τ e, f . Furthermore, the reaction that affects the most τ e, f is reaction 18f: H 2 O 2 (+M) → OH + OH (+M), which consumes H 2 O 2 . Therefore, the addition of H 2 O 2 to the initial mixture, is expected to reinforce the action of the explosive mode. This issue will be explored in the sections that follow. Figure 3 depicts the temperature profiles for the reference case of pure ammonia and all cases of H 2 O 2 addition that were considered. H 2 O 2 addition has a drastic effect on the ignition delay time. In particular, with just 2% molar addition of H 2 O 2 to the fuel, the ignition delay drops significantly by a factor of about 30; i.e., from 2.395 s to 8.375 × 10 −2 s. When increasing H 2 O 2 further to 10%, there appears to be an even more significant ignition delay reduction (by a factor of 248) compared to the reference case. This is very important, since only tenth of H 2 O 2 in the initial mixture can cause a reduction of more than two orders of the time needed for ignition in CI engines. Furthermore, by increasing H 2 O 2 concentration further up to 30% the ignition delay time is decreased by a factor more than 600 compared to the reference case.  Table 3, where the values of the ignition delay times and final temperatures for all cases considered are presented. Table 3. Absolute and relative change of ignition delay times and final temperatures for all cases considered; numbers in parenthesis denote power of ten.

Case t ign (s) Reduced % Change T f inal (K) Increased % Change
Pure Although the difference in final temperature is not as pronounced as the difference in the ignition delay times, the effect of the increasing final temperature changes the dynamics of the process leading to increased NO x species concentration. Figures 4-6 display the evolution of selected species of nitrogen chemistry against time (scaled to ignition delay), for all cases considered.   Figure 4 shows that for all cases considered the addition of H 2 O 2 results in larger mass fractions of NO throughout the ignition process in comparison to the reference case. The final mass fraction of NO increases monotonically with H 2 O 2 content in the initial mixture. In particular, for the cases of 0%, 2%, 10%, 20%, and 30% of H 2 O 2 in the initial mixture the final NO mass fractions are 0.0103, 0.0105, 0.0115, 0.0128, and 0.0144, respectively. Considering the 2% case, it is noted that a 96.50% reduction of ignition delay is achieved with only 2% increase in NO mass fraction. In order to study the dynamics of H 2 O 2 addition, the case of 20% H 2 O 2 addition was further analyzed by examining the developing explosive time scales and by comparing them to those of the reference case. As shown in Figure 7, the fast explosive timescale, which is the one characterizing the evolution of the ignition process, is much faster in the case of 20% H 2 O 2 addition, almost by three orders of magnitude. This indicates that the process in the case of 20% H 2 O 2 addition is faster than the reference one by the same order, confirming the results of Table 3.  Figure 7a indicates that, compared to the reference case (red), the profile of τ e, f in the case of H 2 O 2 addition (purple), does not exhibit any significant deceleration, and remains approximately steady during most of the explosive stage. Close to the end, the process accelerates and a transition to hydrogen chemistry is expected, similarly to the reference case. However, at around 99.5% of the process, there appears to be a sudden and short deceleration of τ e, f , as shown in Figure 7b, of about one order in magnitude. Inspection of Figures 4-6 indicates that in this period there are significant changes in the species profiles of NOx. However, soon, τ e, f regains its fast acceleration towards the end of the explosive stage.

Conclusions
Pure ammonia demonstrates unrealistically long ignition delays that make its application as a fuel impractical. CSP analysis of the isochoric and adiabatic ignition of an ammonia-air mixture with an equivalence ratio of 0.7 showed that the reaction chain initiating reaction 18f: H 2 O 2 (+M) → OH + OH (+M) is related for the most of the process to the explosive mode, supporting its explosive action. Throughout the most part of the explosive stage, this reaction competes with the ammonia dissociation reaction 156b: N 2 H 4 + OH ← NH 3 + H 2 NO for the production of OH radical, building the radical pool needed in the final part of the explosive stage. The action of this reaction is reinforced by that of the ammonia dissociation reaction N 2 H 4 + OH ← NH 3 + H 2 N, which along with 18f produces OH radicals. The major opposition to the ignition arises from the chain terminating reaction 208f H 2 NO + NH 2 → HNO+NH 3 .
The variable related to the explosive mode for the largest portion of the explosive stage was the temperature. This indicates that, contrary to hydrocarbons, ammonia demonstrates a brief chemical runaway and a long thermal one. H 2 O 2 was the species pointed the largest portion of the explosive stage. Therefore, this species was selected as an additive in order to shorten ignition delay. It was demonstrated that even a small addition of H 2 O 2 (e.g., 2%) could decrease the ignition delay by as much as factor of 30. This came with a relatively small increase in NO emissions by just 2%. These results show that H 2 O 2 addition can drastically decrease ammonia ignition delay without increasing NO emissions substantially. Practical application of H 2 O 2 addition will only be possible if the issues are addressed that relate to the corrosive action and the stability of this molecule, which tends to dissociate to H 2 O and O 2 after long term storage.
The mechanism by which H 2 O 2 accelerates the process will be studied in a future work using CSP, by considering different NH 3 /H 2 O 2 mixtures at various initial temperatures and pressures.

Abbreviations
The following abbreviations are used in this manuscript: