Simple Program for Step-by-Step Time Integration in Chemical Kinetics, Applied to Simple Model for Hydrogen Combustion

Panagis G. Papadopoulos¹, Christopher G. Koutitas², Panos D. Kiousis³, Christos G. Karayannis⁴, Yannis N. Dimitropoulos⁵

¹Structural Engineering Laboratory, Civil Engineering Department, Aristotle University of Thessaloniki, Thessaloniki, Greece.
²Computational Fluid Engineering Laboratory, Civil Engineering Department, Aristotle University of Thessaloniki, Thessaloniki, Greece.
³Structural and Geotechnical Engineering Laboratory, Civil and Environmental Engineering Department, Colorado School of Mines, Golden, USA.
⁴Concrete Structures Laboratory, Civil Engineering Department, Democritus University of Thrace, Xanthi, Greece.
⁵Physical Chemistry Laboratory, Chemistry Department, University of Ioannina, Ioannina, Greece.

DOI: 10.4236/ojpc.2020.102006                

Abstract

A simple algorithm is proposed for step-by-step time integration of stiff ODEs in Chemical Kinetics. No predictor-corrector technique is used within each step of the algorithm. It is assumed that species concentrations less than 10^-6 mol·L^-1 do not activate any chemical reaction. So, within each step, the time steplength Δt of the algorithm is determined from the fastest reaction rate maxR by the formula Δt = 10^-6mol·L^-1/max R. All the reversible elementary reactions occur simultaneously; however, by a simple book-keeping technique, the updating of species concentrations, within each step of the algorithm, is performed within each elementary reaction separately. The above proposed simple algorithm for Chemical Kinetics is applied to a simple model for hydrogen combustion with only five reversible elementary reactions (Initiation, Propagation, First and Second Branching, Termination by wall destruction) with six species (H₂, O₂, H, O, HO, H₂O). These five reversible reactions are recommended in the literature as the most significant elementary reactions of hydrogen combustion [1] [2]. Based on the proposed here simple algorithm for Chemical Kinetics, applied to the global mechanism of proposed five reversible elementary reactions for hydrogen combustion, a simple and short computer program has been developed with only about 120 Fortran instructions. By this proposed program, the following are obtained: 1) The total species concentration of hydrogen combustion, starting from the sum of initial reactants concentrations [H₂] + [O₂], gradually diminishes, due to termination reaction by wall destruction, and tends to the final concentration of the product [H₂O], that is to the 2/3 of its initial value, in accordance to the established overall stoichiometric reaction of hydrogen combustion 2H₂ + O₂ → 2H₂O. 2) Time-histories for concentrations of main species H₂, O₂, H, H₂O of hydrogen combustion, in explosion and equilibrium regions, obtained by the proposed program, are compared to corresponding ones obtained by accurate computational studies of [3]. 3) In the first step of the algorithm, the only nonzero species concentrations are those of reactants [H₂], [O₂]. So, the maximum reaction rate is that of the forward initiation reaction max R = R_if = k_if[H₂] [O₂], where the rate constant k_if is very slow. Thus, the first time steplength Δt₁ = 10^-6mol·L^-1/max R results long in sec. After the first step, the sequences of all the following Δt’s are very short, in μsec. So, the first time steplength Δt₁ can be considered as ignition delay time. 4) It is assumed that explosion corresponds to ignition delay time Δt₁ < 10 sec. So, the curve on T-P plane, obtained by proposed program for Δt₁ = 10 sec., can be considered as explosion limit curve. This curve is compared to the corresponding one obtained by the accurate computational studies of [2].

Keywords

Chemical Kinetics, Step-by-Step Time Integration, Elementary Reactions, Hydrogen Combustion, Ignition Delay, Explosion

1. Introduction

Refined algorithms of high accuracy have been proposed for step-by-step time integration of stiff ODEs in Chemical Kinetics [2] - [8]. Here a simplified algorithm is proposed.

The above proposed simple algorithm for Chemical Kinetics is applied to hydrogen combustion. Accurate models for hydrogen combustion have been proposed in the literature with 29 reversible elementary reactions, in [2] and 25 reversible elementary reactions in [3], both with 9 species H₂, O₂, H, O, HO, HO₂, O₃, H₂O₂, H₂O. Here a simple model is proposed, with only five reversible elementary reactions (Initiation, Propagation, First and Second Branching, Termination by Wall Destruction) with six species (H₂, O₂, H, O, HO, H₂O). The above five reversible elementary reactions are recommended in the literature as the most significant elementary reactions in hydrogen combustion [1] [2].

Based on the proposed here simple algorithm for Chemical Kinetics, applied on the global mechanism of the proposed five reversible elementary reactions of hydrogen combustion, a simple and short computer program has been developed.

Results, obtained by the above program, are compared with corresponding ones of published computational studies, based on refined algorithms of high accuracy for Chemical Kinetics and accurate models of hydrogen combustion.

2. The Initial Value Problem in Chemical Kinetics

Usually a stoichiometric chemical reaction is the resultant of a number of elementary, bi-molecular, reversible reactions of the form

$\text{D}+\text{E}\rightleftarrows \text{F}+\text{G}$

Thanks to the bi-molecular nature of chemical collisions of above reactions, their rates R_f, R_b, in forward and backward direction, respectively, can be written in the simple form

$\begin{array}{l}{R}_f=\text{d}{c}_f/\text{d}t=k_f\left[\text{D}\right]\left[\text{E}\right],\\ R_b=\text{d}{c}_b/\text{d}t=k_b\left[\text{F}\right]\left[\text{G}\right]\end{array}$, (1)

where c_f, c_b total species concentrations of left and right part of above reaction, respectively, t time, k_f, k_b rate constants in forward and backward direction, respectively, and [D], [E], [F], [G] concentrations of species D, E, F, G, respectively.

The above first order differential Equation (1), with respect to time t, constitute a system of ODEs (ordinary differential equations). Because, usually, many rate constants are very fast, the resulting time scales are very short, so the ODEs are called stiff. Also, because of large variation in orders of magnitude of rate constants k, the time t scales largely vary, too, from sec, msec = 10⁻³ sec., μsec = 10⁻⁶ sec. up to nsec = 10⁻⁹ sec.

Here, only gases are considered obeying the ideal gas law:

$PV/RT={\displaystyle \sum n}$,

where R = 8.314 J·mol⁻¹·K⁻¹ gas constant, T temperature, P pressure, V volume and ${\displaystyle \sum n}$ sum of moles of all species of the considered gas mixture. If the initial number n_i of moles of any i-th gas of the mixture is known, then the initial volume of the mixture is

$V={\displaystyle \sum n_i}RT/P$

and the initial concentration of any i-th gas is

$c_i=n_i/V$ (2)

The system of ODEs of Equation (1) (first derivatives of reactions concentrations with respect to time t), where $c_f=\left[\text{D}\right]+\left[\text{E}\right]$ and $c_b=\left[\text{F}\right]+\left[\text{G}\right]$, together with the initial values of species concentrations of Equation (2), constitute the Initial Value Problem of Chemical Kinetics.

3. Simple Proposed Step-by-Step Time Integration Algorithm in Chemical Kinetics

For the step-by-step time integration of the Initial Value Problem of Chemical Kinetics, described in previous Section 2, a simple algorithm is here proposed. No P-C (Predictor-Corrector) technique is used, within each step of the algorithm. It is reasonably assumed that, for species concentrations less than 10⁻⁶ mol·L⁻¹, no chemical reaction is activated. So, within each step, the time-steplength Δt of the algorithm is determined from the fastest reaction rate maxR by the formula $\Delta t={10}^{-\text{6}}\text{mol}\cdot {\text{L}}^{-\text{1}}/\max R$. Although all elementary reactions occur simultaneously, a simple book-keeping technique is proposed, by which, within each step of the algorithm, the updating of species concentrations is performed for every elementary reaction separately. In order for the proposed algorithm, to become more clearly understood, it will be applied to a specific case, that of a simple proposed model for hydrogen combustion, which is presented in the following sections 4 up to 8.

4. Simple Proposed Model for Hydrogen Combustion

For the computational study of hydrogen combustion, a simplified model, with only five elementary reversible reactions and six species (reactants H₂, O₂, intermediate species H, O, HO, product H₂O), is proposed. The five proposed reactions also recommended in [1], page 607,are the following:

1) Initiation ${\text{H}}_{\text{2}}+{\text{O}}_{\text{2}}\rightleftarrows \text{HO}+\text{HO}$,

2) Propagation $\text{HO}+{\text{H}}_{\text{2}}\rightleftarrows \text{H}+{\text{H}}_{\text{2}}\text{O}$,

3) First Branching $\text{H}+\text{HO}\rightleftarrows \text{O}+{\text{H}}_{\text{2}}$,

4) Second Branching $\text{O}+\text{HO}\rightleftarrows \text{H}+{\text{O}}_{\text{2}}$,

5) Termination by Wall Destruction $\text{H}+\text{H}+\text{M}\rightleftarrows {\text{H}}_{\text{2}}+\text{M}$, where the species M is inert.

5. Rate Constants for the Proposed Five Reversible Elementary Reactions of Hydrogen Combustion

The coefficients A, B, C, for the modified Arrhenius expression of Rate Constants $k\left(T\right)=A{T}^B{\text{e}}^{C/T}$, for the proposed five reversible elementary reactions of hydrogen combustion, based on test data, have been obtained from the literature [2] and are presented in Table 1. Similar values are given by [3]. The units of Rate Constants k(T) of the nine bi-molecular reactions are cm³/(molecules × sec). Here, for species concentrations, units mol·L⁻¹ are used. So, a corrective factor (6.022 × 10²³ molecules/mol)/(10³ cm³/L) is needed in the nine bi-molecular reactions. Only for the last backward termination reaction, which is uni-molecular, the unit of Rate Constant k(T) is sec⁻¹ and no corrective factor is needed. Also, only for the coefficient A of forward termination reaction, a value different from that given by [2] is used here after trials, A = 1.80 × 10⁻¹¹, as it is considered that gives better results.

6. Step-by-Step Time Integration of Individual Reversible Elementary Reactions

From the five proposed reversible elementary reactions of hydrogen combustion, the three (propagation, first and second branching) are regular bi-molecular reactions of the form $\text{D}+\text{E}\rightleftarrows \text{F}+\text{G}$ and their step-by-step time integration, by the proposed here algorithm, is described by the flow-chart of Figure 1. Based on this flow-chart, a simple and short computer program has been developed with only about 30 Fortran instructions.

Table 1. Coefficients A, B, C, for the Rate Constants $k\left(T\right)=A{T}^B{\text{e}}^{C/T}$ of proposed five reversible elementary reactions for hydrogen combustion, obtained from [2].

In the flow-chart of Figure 1, first temperature T and pressure P are read, as well as the initial number n_i of moles of every i-th species and the time t_u, for which the algorithm has to be interrupted. Then, the coefficients A, B, C of rate constants KF, KB, for forward and backward reaction, respectively, are read and KF, KB are determined. The initial volume of air-mixture $V={\displaystyle \sum n}\times RT/P$ is found and the initial concentration $c_i=n_i/V$ of every i-th species is determined. The species concentrations are noted by D, E, F, G.

Within each step of the algorithm, first, the reaction rates $RF=KF\times D\times E$ and $RB=KB\times F\times G$, of forward and backward reaction, respectively, are determined. From the maximum MAXR of these two reaction rates, the present time steplength $Dt={10}^{-\text{6}}\text{mol}\cdot {\text{L}}^{-\text{1}}/\text{MAXR}$ is found. The previous time instant t and present time steplength Dt are written. The time t is increased by Dt, $t\leftarrow t+\text{Dt}$.

The concentration increment of forward reaction is $DC=RF\times Dt$ and the corresponding updating of species concentrations: $D\leftarrow D-DC/2$, $E\leftarrow E-DC/2$, $F\leftarrow F+DC/2$, $G\leftarrow G+DC/2$.

The concentration increment of backward reaction is $DC=RB\times Dt$ and the corresponding updating of species concentrations: $D\leftarrow D+DC/2$, $E\leftarrow E+DC/2$, $F\leftarrow F-DC/2$, $G\leftarrow G-DC/2$.

If any i-th species concentration results negative c_i < 0, it is put equal to zero c_i = 0.

The present step of the algorithm, time instant t and present species concentrations D, E, F, G are written as output.

Then, we go to the next step of the algorithm. When the time t exceeds the predetermined maximum time t_u, the algorithm is interrupted.

The above algorithm runs for temperature T = 1000 K and pressure P = 1.0 atm. First, for forward reaction, with initial number of moles: n = 1.0 mole for each one of the reactants D, E and n = 0.0 mole for each one of products, F, G.

Figure 1. Flow-chart for step-by-step time integration of an individual reversible bi-molecular elementary reaction of the form $\text{D}+\text{E}\rightleftarrows \text{F}+\text{G}$, by the proposed algorithm.

The initial volume of gas mixture is

$\begin{array}{c}V={\displaystyle \sum n}\times RT/P\\ =2.0\text{mol}\times 8.314\text{N}\cdot \text{m}\cdot {\text{mol}}^{-1}\cdot {\text{K}}^{-1}\times 1000\text{K}/\left({10}^5\text{N}/{\text{m}}^{\text{2}}\right)\\ =0.1663{\text{m}}^3=166.3\text{L}\end{array}$

So, the initial concentrations are $c=1.0\text{mol}/166.3\text{L}=6.014\times {10}^{-3}\text{mol}\cdot {\text{L}}^{-1}$ for each one of the reactants D, E and zero c = 0.0 for each one of the products F, G. Then, the program runs, for a second time, in the backward direction, with initial concentrations of new reactants $\text{F}=\text{G}=6.014\times {10}^{-3}\text{mol}\cdot {\text{L}}^{-\text{1}}$ and of new products $\text{D}=\text{E}=0.0$. That is, the second time, the initial concentrations of reactants and products have been interchanged.

The above two runs of the program, forward and backward, are applied for the specific rate constants (Table 1) of the three proposed bi-molecular reversible elementary reactions of hydrogen combustion (propagation, first and second branching) and the time-histories, of these three reversible elementary reactions, forward and backward, have been obtained and are presented in Figures 2(b)-(d), respectively.

It is observed, in these figures that, if, by going from forward to backward reaction, the initial concentrations of reactants and products are interchanged, then the same equilibrium state is reached, in both cases. The fastest of these reactions is the forward second branching, with time scale in nsec = 10⁻⁹ sec. Also, fast are the other reactions in μsec = 10⁻⁶ sec., except of the backward propagation, which is slow in msec = 10⁻³ sec.

The Initiation reversible elementary reaction exhibits collisions between molecules of the same species HO. So, some modifications are needed in the above presented program: The backward reaction rate is $\text{RB}=\text{KB}\times \text{HO}\times \text{HO}/\text{2}$. The updating of HO concentration is $\text{HO}\leftarrow \text{HO}+\text{DC}$ in the forward and $\text{HO}\leftarrow \text{HO}-\text{DC}$ in the backward direction. Also, the initial concentration in the backward reaction, for HO, is 2 × 6.014 × 10⁻³ molL⁻¹. With these modifications, in the program, the time histories of species, for the Initiation reaction, forward and backward, have been obtained and shown in Figure 2(a).

The termination reversible elementary reaction, in addition to the modifications needed because it also exhibits collisions between the same species H, it needs further modifications of the program, because its right part is uni-molecular. So, the backward reaction rate is $\text{RB}=\text{KB}\times \text{H2}$. The updating of H2 concentration is $\text{H2}\leftarrow \text{H2}+\text{DC}/2$ in the forward and $\text{H2}\leftarrow \text{H2}-\text{DC}$ in the backward reaction. And the initial concentrations are H = 12.028 × 10⁻³ molL⁻¹, H2 = 0. in forward direction and H2 = 6.014 × 10⁻³ mol·L⁻¹, H = 0 in backward direction.

Also, only in the termination reaction we have to determine the total concentration S = H + H2 within each step of the algorithm. In the other reactions, the total concentration S remains constant, whereas, in this termination reaction by wall destruction, the S gradually diminishes and tends to the 1/2 of its initial value.

Figure 2. Time-histories of species, for the specific examples, of the individual proposed five reversible elementary reactions for hydrogen combustion. (a) Initiation, (b) Propagation, (c) First Branching, (d) Second Branching (e) Termination by Wall Destruction.

The so modified program run and the time-histories of species of termination reaction have been obtained and are shown in Figure 2(e).

In the time histories of species of initiation and termination reaction, it is again observed that, by going from forward to backward reaction, if the initial concentrations of products and reactants are interchanged, the same equilibrium state is reached in both cases. However, in the termination reaction, the comment “by following the stoichiometry of the reaction” has to be added, because, in termination reaction, the left part is bi-molecular and the right part uni-molecular.

It is observed that the backward termination reaction is extremely slow, with hr time units. The Initiation reactions, forward and backward, are slow, with time in sec. On the contrary, the forward termination reaction is fast with time in μsec = 10⁻⁶ sec.

7. Global Mechanism of All Simultaneous Five Proposed Reversible Elementary Reactions. Simple Book-Keeping Technique for Updating of Species Concentrations, within Each Elementary Reaction Seperately

In the flow-chart of Figure 3, the proposed algorithm, for the global mechanism of the five proposed reversible elementary reactions of hydrogen combustion, is described. Based on this flow-chart, a simple and short computer program has been developed, with only about 120 Fortran instructions.

First, the temperature T and pressure P are read, as well as the initial number n_i of moles of every i-th species and the time t_u, for which the algorithm has to be interrupted.

Then, the coefficients A, B, C of rate constants $k\left(T\right)=A{T}^B{\text{e}}^{C/T}$ are read and the rate constants of five proposed elementary reactions, forward and backward, are determined:

The initial volume $V={\displaystyle \sum n}\times RT/P$ of the air-mixture is found, where ${\displaystyle \sum n}$ initial total number of moles of all gas species. Then the initial concentration $c_i=n_i/V$ of every i-th species is determined. The symbol of six species concentrations are H2, O2, H, O, HO, H2O.

Within each step of the algorithm, the reaction rates are determined:

Figure 3. Flow-chart for the global mechanism of all simultaneous five proposed reversible elementary reactions of hydrogen combustion. Simple book-keeping technique, within each step, for updating of species concentrations, within each elementary reaction seperately.

The maximum reaction rate MAXR of the above 5 × 2 = 10 reaction rates is determined. The present time steplength is found $\text{Dt}={10}^{-\text{6}}\text{mol}\cdot {\text{L}}^{-\text{1}}/\text{MAXR}$.

The previous time instant t and present time steplength Dt are written. The time t is increased by Dt, $t\leftarrow t+\text{Dt}$.

By a simple book-keeping technique, within each step, the species concentrations are updated, within each reversible elementary reaction separately, where DC is the total concentration increment of an elementary reaction, forward or backward:

If any i-th species concentration c_i results negative, c_i < 0.0, it is put equal to zero, c_i = 0.0.

The total sum of concentrations of all species is formed for the present step,

$\text{S}=\text{H2}+\text{O2}+\text{H}+\text{O}+\text{HO}+\text{H2O}$,

which is gradually reduced because of termination reaction by wall destruction.

The present step of algorithm, time t, concentrations of species H2, O2, H, O, HO, H2O and total concentration S are written as output of present step of algorithm.

Then, we continue with the next step of the algorithm. When the time t exceeds the pre-determined maximum time t_u, the algorithm is interrupted.

8. Application

The proposed program, for the global mechanism of all the five proposed reversible elementary reactions of hydrogen combustion, described by the flow-chart of Figure 3 of previous Section 7, is applied to the specific case, also studied in [3], of an air mixture H₂:O₂:N₂ with initial amounts of 2:1:3.76 moles, where nitrogen N₂ is inert to the examined mechanism. So, the initial total amount is ${\displaystyle \sum n}=2.0+1.0+3.76=6.76\text{moles}$. Temperature T = 1000 K and pressure P = 1.0 atm are considered. So, the initial volume of the gas mixture is

$\begin{array}{c}V={\displaystyle \sum n}\times RT/P\\ =6.76\text{moles}\times 8.314\text{N}\cdot \text{m}\cdot {\text{mol}}^{-1}\cdot {\text{K}}^{-1}\times 1000\text{K}/\left({10}^5\text{N}/{\text{m}}^{\text{2}}\right)\\ =0.5620{\text{m}}^3=562.0\text{L}\end{array}$

and the initial concentrations of reactants are $\left[{\text{H}}_{\text{2}}\right]=2.0\text{mol}/562.0\text{L}=3.559\times {10}^{-3}\text{mol}\cdot {\text{L}}^{-1}$, $\left[{\text{O}}_{\text{2}}\right]=1.0\text{mol}/352.0\text{L}=1.779\times {10}^{-3}\text{mol}\cdot {\text{L}}^{-1}$ and those of remaining species H, O, HO, H₂O are zero. The proposed program run, with these initial conditions, and produced, as output, the time-histories of species concentrations of H₂, O₂, H, O, HO, H₂O.

First, a macroscopic view of the above time-histories of species concentrations is shown in ignition delay, explosion and chemical equilibrium regions, with time scale in sec, in Figure 4(a). The duration of ignition delay is equal to the first time steplength Δt₁, as will be discussed below with help of Figure 4(c).

In Figure 4(b), a microscopic view of the same time-histories of species concentrations is shown, in the explosion and equilibrium regions, with time scale in msec = 10⁻³ sec. The same microscopic view of this application is presented by [3] and will be compared with the corresponding results of present work in following Section 9.2.

In the first step of the algorithm, the only nonzero species concentrations are those of reactants [H₂], [O₂]. So, the fastest reaction rate is that of forward initiation elementary reaction $\max R=R_{if}=k_{if}\left[{\text{H}}_{\text{2}}\right]\left[{\text{O}}_{\text{2}}\right]$, where the rate constant k_if is very slow. Thus, the first time steplength $\Delta t_{\text{1}}={10}^{-\text{6}}\text{mol}\cdot {\text{L}}^{-\text{1}}/\max R$ results long, Δt₁ = 0.3009 sec in the present application. In Figure 4(c), the sequence of the following time-steplenths Δt are shown. Initially they are in μsec = 10⁻⁶ sec time units. Then, they gradually diminish, first by oscillating, then by smooth variation and tend to nsec = 10⁻⁹ sec. After entering equilibrium region, they begin to increase gradually again. The maximum of all time-steplengths Δt, after the first one Δt₁, in the explosion region, is Δt₄ = 11.90 μsec. So, the first time steplength Δt₁ = 0.3009 sec is 0.3009 sec/11.90 μsec = 25,286 times longer than the maximum Δt₄ = 11.90 μsec of all the following time-steplengths Δt, in the explosion region; thus, the first time steplength Δt₁ = 0.3009 sec can be considered as the Ignition Delay Time.

Figure 4. Output of application (a) Macroscopic view of species concentrations of H₂, O₂, H, O, HO, H₂O, in sec time unit, in ignition delay, explosion and equilibrium regions (b) Microscopic view of species concentrations of H₂, O₂, H, O, HO, H₂O in msec = 10⁻³ sec time unit, in explosion and equilibrium regions c. Variation of the steplenth Δt of the algorithm in the first 70 μsec of explosion region.

9. Comparisons

9.1. Gradual Reduction of Total Concentration to 2/3 of Its Initial Value in Accordance to Overall Stoichiometric Reaction

By the proposed algorithm, it is obtained that the total concentration of all species ${\displaystyle \sum c_i}$ starting from the sum of initial concentrations of reactants $\left[{\text{H}}_{\text{2}}\right]+\left[{\text{O}}_{\text{2}}\right]=\text{5}.\text{339}\times {10}^{-\text{3}}\text{mol}\cdot {\text{L}}^{-\text{1}}$ in present application, gradually diminishes, due to the termination elementary reaction, by wall destruction, and tends to the final concentration of the product $\left[{\text{H}}_{\text{2}}\text{O}\right]=\text{3}.559\times {10}^{-3}\text{mol}\cdot {\text{L}}^{-\text{1}}$, in the present application, which is the 2/3 of its initial value, as shown in Figure 5. For time t = 2.0 sec., ${\displaystyle \sum n_i}$ is very close to the value 3.559 × 10⁻³ mol·L⁻¹. This is in accordance to the global stoichiometric chemical reaction of hydrogen combustion ${\text{2H}}_{\text{2}}+{\text{O}}_{\text{2}}\to {\text{2H}}_{\text{2}}\text{O}$, which implies 2.0 mol + 1.0 mol = 3.0 mol→2.0 mol.

9.2. Comparison of Time-Histories of Main Species Concentrations of H₂, O₂, H, H₂O, in Explosion and Equilibrium Regions, to Published Computational Data by [3]

Time-histories of main species concentrations [H₂], [O₂], [H], [H₂O] of hydrogen combustion, in explosion and equilibrium regions, obtained by the proposed algorithm and model, are compared to corresponding ones, obtained by the published accurate computational studies of [3], page 64, as shown in Figure 6 of present work. A temporal shift, about 0.2 msec, is observed between corresponding species concentrations time-history curves. Particularly, for the vertex of critical H curve, the temporal shift is only 0.135 msec. So, the approximation, between corresponding time-history curves, can be considered satisfactory.

9.3. Comparison of Explosion Limit Curve, on T-P Plane, to Published Computational Results of [2]

For an initial mixture H₂:O₂, of 2.0:1.0 moles, temperature T = 1000 K, pressure P = 1.0 atm, also studied in [2], curves, of equal ignition delay time Δt₁, have been drawn, by the proposed program, on the T-P (temperature-pressure) plane, for values of Δt₁: 0.1, 1, 10, 100, 1000 sec., as shown in Figure 7. It is reasonably assumed that explosion corresponds to ignition delay times Δt₁ less than 10 sec. So, the curve for Δt₁ = 10 sec. is considered as the explosion limit curve and is compared to the corresponding curve, obtained by the published accurate computational studies of [2], page 34, after transforming their logarithmic scale for

Figure 5. Time-history of total species concentration ${\displaystyle \sum n}$, of hydrogen combustion, obtained by the proposed program, in the present application, tending to 2/3 of its initial value, in accordance to the overall stoichiometric reaction ${\text{2H}}_{\text{2}}+{\text{O}}_{\text{2}}\to {\text{2H}}_{\text{2}}\text{O}$.

Figure 6. Comparison between time-histories of main species concentrations [H₂], [O₂], [H], [H₂O] of hydrogen combustion, in explosion and equilibrium regions, obtained by the proposed program, to corresponding ones obtained by the accurate published computational studies of [3], page 64.

Figure 7. Curves of equal Δt₁ = ignition delay time, drawn on T-P (temperature-pressure) plane, by proposed program, for Δt₁: 0.1, 1, 10, 100, 1000 sec. Explosion limit curve, for ignition delay time Δt₁ = 10 sec., is compared to corresponding one of accurate published computational studies by [2], page 34, for high pressures P: 0.4 up to 2.0 atm, after transforming their logarithmic scale of pressure P to regular scale.

pressure P to regular scale. And a satisfactory approximation is observed, between corresponding explosion limit curves, for high pressures P: 0.4 up to 2.0 atm, as shown in Figure 7.

10. Conclusions

1) The proposed here simple program for Chemical Kinetics, applied to the simple proposed model for hydrogen combustion, run, first, for the individual five proposed reversible elementary reactions. And the same equilibrium states reached from forward and backward reactions were observed, as well as a large variation of time scales ranging from sec up to nsec = 10⁻⁹ sec.

2) By the proposed program, it is obtained that the total species concentration of hydrogen combustion, starting from the sum of initial concentrations of reactants [H₂] + [O₂], gradually diminishes, due to termination reaction by wall destruction, and tends to the final concentration of product [H₂O], that is to the 2/3 of its initial value, in accordance to the overall stoichiometric reaction ${\text{2H}}_{\text{2}}+{\text{O}}_{\text{2}}\to {\text{2H}}_{\text{2}}\text{O}$.

3) Time-histories of concentrations of the main species H₂, O₂, H, H₂O of hydrogen combustion, in explosion and equilibrium regions, obtained by the proposed program, are compared to corresponding ones, obtained by the accurate published computational data of [3]. The approximation between corresponding time-history curves is satisfactory.

4) In the first step of the proposed algorithm, the only nonzero species concentrations are those of reactants [H₂], [O₂], so the maximum reaction rate is that of forward initiation reaction $\max R=R_{if}=k_{if}\left[{\text{H}}_{\text{2}}\right]\left[{\text{O}}_{\text{2}}\right]$, where the rate constant k_if is very slow. Thus, the first time steplength $\Delta t_{\text{1}}={10}^{-\text{6}}\text{mol}\cdot {\text{L}}^{-\text{1}}/\max R$ results long in sec time unit. After the first step, the sequences of all the following Δt’s are very short in μsec = 10⁻⁶ sec. So, the first time steplength Δt₁ can be considered as ignition delay time.

5) It is assumed that Explosion corresponds to ignition delay time Δt₁ < 10 sec. So, the curve on T-P (temperature – pressure) plane, obtained by the proposed program, for Δt₁ = 10 sec. can be considered as Explosion limit curve. This curve is compared to corresponding one obtained by the accurate published computational studies of [2], and a satisfactory approximation is observed between the corresponding Explosion limit curves.

6) As regards future research on hydrogen combustion, “hydrogen as fuel” is recently a very active research field, with significant practical applications, as shown in literature [9] - [19].

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Brady, J.E., Russell, J.W. and Holum, J.R. (2000) Chemistry: Matter and Its Changes. 3rd Edition, John Wiley & Sons Inc., New York. https://www.amazon.com/Chemistry-Study-Matter-Changes-2000-01-19/dp/B01F81WBQG
[2]	Burks, T.L. and Oran, E.S. (1981) A Computational Study of the Chemical Kinetics of Hydrogen Combustion. Naval Research Laboratory, NASA, Washington DC. https://apps.dtic.mil/dtic/tr/fulltext/u2/a094348.pdf
[3]	Mott, D.R. (1999) New Quasi-Steady-State and Partial Equilibrium Methods for Integrating Chemically Reacting Systems. Ph.D. Thesis Dissertation, Aerospace Engineering and Scientific Computing, University of Michigan, Ann Arbor, MI. https://pdfs.semanticscholar.org/f26c/513c99e2fe330bd78f209d262a3502c982f2.pdf
[4]	Young, T.R. and Boris, J.P. (1977) A Numerical Technique for Solving Stiff Ordinary Differential Equations Associated with the Chemical Kinetics of Reactive-Flow Problems. The Journal of Physical Chemistry, 81, 2424-2427. https://doi.org/10.1021/j100540a018
[5]	Young Jr., T.R. (1980) CHEMEQ: A Subroutine for Solving Stiff Ordinary Differential Equations. Laboratory for Computational Physics. Naval Research Laboratory, Washington DC. https://apps.dtic.mil/dtic/tr/fulltext/u2/a083545.pdf
[6]	Radhakrishnan, K. (1984) Comparison of Numerical Techniques for Integration of Stiff Ordinary Differential Equations Arising in Combustion Chemistry. Lewis Research Center, Cleveland, OH. https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19850001758.pdf
[7]	Mott, D.R. and Oran, E.S. (2001) CHEMEQ2 A Solvent for the Stiff Ordinary Differential Equations of Chemical Kinetics. Naval Research Laboratory, Washington DC. https://apps.dtic.mil/dtic/tr/fulltext/u2/a392490.pdf https://doi.org/10.21236/ADA385560
[8]	Sichel, M., Tonello, N., Oran, E.S. and Jones, D.A. (2002) A Two-Step Kinetics Model for Numerical Simulation of Explosions and Detonations in H2-O2 Mixtures. Proceedings of the Royal Society London A, 458, 49-82. https://doi.org/10.1098/rspa.2001.0853
[9]	Westbrook, C.K. and Dryer, F.L. (1984) Chemical Kinetics Modeling of Hydro-Carbon Combustion. Progress in Energy and Combustion Science, 10, 1-57. https://doi.org/10.1016/0360-1285(84)90118-7
[10]	Westbrook, C.K., Mizobushi, Y., Poinsot, T., Smith, P.J. and Warnatz, J. (2005) Computational Combustion. Proceedings of the Combustion Institute, 30, 125-157. https://doi.org/10.1016/j.proci.2004.08.275
[11]	Law, C.K. (2006) Combustion Physics. Cambridge University Press, Cambridge.
[12]	Wang, X. and Law, C.K. (2013) An Analysis of the Explosion Limits of Hydrogen-Oxygen Mixtures. Journal of Chemical Physics, 138, Article ID: 134305. https://doi.org/10.1063/1.4798459
[13]	Liang, W. and Law, C.K. (2018) An Analysis of the Explosion Limits of Hydrogen-Oxygen Mixtures with Non-Linear Chain Reactions. Physical Chemistry, Chemical Physics, 20, 742-751. https://doi.org/10.1039/C7CP05639G
[14]	Zhao, P. (2018) Detailed Kinetics in Combustion. Simulation: Manifestation, Model Reduction and Computational Diagnostics. In: De, S., Agarwal, A., Chaudhuri, S. and Sen, S., Eds., Modeling and Simulation of Turbulent Combustion, Springer, Berlin. https://doi.org/10.1007/978-981-10-7410-3_2
[15]	Friedrich, A., Grune, J., Sempert, K., Kuznetsov, M. and Jordan, T. (2019) Hydrogen Combustion Experiments in a Vertical Semi-Confined Channel. International Journal of Hydrogen Energy, 44, 9041-9049. https://doi.org/10.1016/j.ijhydene.2018.06.098
[16]	Kuznetsov, M. and Grune, J. (2019) Experiments on Combustion Regimes for Hydrogen/Air Mixtures in a Thin Layer Geometry. International Journal of Hydrogen Energy, 44, 8727-8742. https://doi.org/10.1016/j.ijhydene.2018.11.144
[17]	Faizal, M., Chuah, L.S., Lee, C., Hameed, A., Lee, J. and Shankar, M. (2019) Review of Hydrogen Fuel for Internal Combustion Engines. Journal of Mechanical Engineering Research and Developments, 42, 35-46. https://doi.org/10.26480/jmerd.03.2019.35.46
[18]	Wikipedia (2020) Hydrogen Fuel. https://en.wikipedia.org/wiki/Hydrogen_fuel
[19]	Enviroment and Energy Conferece. Hydrogen Days (2020) Organizer: Czech Hydrogen Technology Platform (HYTEP). To Be Held 25-27 March 2020, Prague, Czech Republic. https://www.hydrogendays.cz/2020/