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

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 . 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 (In-itiation, Propagation, First and Second Branching, Termination by wall destruction) with six species (H 2 , O 2 , H, O, HO, H 2 O). These five reversible reactions are recommended in sum of initial reactants concentrations [H 2 ] + [O 2 ], gradually diminishes, due to termination reaction by wall destruction, and tends to the final concentration of the product [H 2 O], that is to the 2/3 of its initial value, in accordance to the established overall stoichiometric reaction of hydrogen combustion . 2) Time-histories for concentrations of main species H 2 , O 2 , H, H 2 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 2 ], [O 2 ]. So, the maximum reaction rate is that of the forward initiation reaction , where the rate constant k if is very slow. Thus, the first time steplength 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 1 can be considered as ignition delay time. 4) It is assumed that explosion corresponds to ignition delay time Δt 1 < 10 sec. So, the curve on T-P plane, obtained by proposed program for Δt 1 = 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].


sum of initial reactants concentrations [H 2 ] + [O 2 ]
, gradually diminishes, due to termination reaction by wall destruction, and tends to the final concentration of the product [H 2 O], that is to the 2/3 of its initial value, in accordance to the established overall stoichiometric reaction of hydrogen combustion 2  , where the rate constant k if is very slow. Thus, the first time steplength 6 1 1 10 mol L max t 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 1 can be considered as ignition delay time. 4) It is assumed that explosion corresponds to ignition delay time Δt 1 < 10 sec. So, the curve on T-P plane, obtained by proposed program for Δt 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 2 , O 2 [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.

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 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 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, re- 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 −3 sec., μsec = 10 −6 sec. up to nsec = 10 −9 sec.
Here, only gases are considered obeying the ideal gas law: where R = 8.314 J·mol −1 ·K −1 gas constant, T temperature, P pressure, V volume and 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 i V n RT P = ∑ and the initial concentration of any i-th gas is The system of ODEs of Equation (1)

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 −6 mol·L −1 , 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 Open Journal of Physical Chemistry . 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.

Rate Constants for the Proposed Five Reversible Elementary Reactions of Hydrogen Combustion
The coefficients A, B, C, for the modified Arrhenius expression of Rate Con- , 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 3 /(molecules × sec). Here, for species concentrations, units mol·L −1 are used. So, a corrective factor (6.022 × 10 23 molecules/mol)/(10 3 cm 3 /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 −1 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 −11 , as it is considered that gives better results.  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

Step-by-Step Time Integration of Individual Reversible Elementary Reactions
and the corresponding updating of species concentrations: 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. Open Journal of Physical Chemistry 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)  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. 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 equili-

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.  The maximum reaction rate MAXR of the above 5 × 2 = 10 reaction rates is determined. The present time steplength is found 6 1 Dt 10 mol L MAXR

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 2 :O 2 :N 2 with initial amounts of 2:1:3.76 moles, where nitrogen N 2 is inert to the examined mechanism. So, the initial total amount is 2.0 1.0 3.76 6.76 moles n = + + = ∑ . Temperature T = 1000 K and pressure P = 1.0 atm are considered. So, the initial volume of the gas mixture is ( )  Figure 4(a). The duration of ignition delay is equal to the first time steplength Δt 1 , 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 −3 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.

Comparison of Explosion Limit Curve, on T-P Plane, to Published Computational Results of [2]
For an initial mixture H 2 :O 2 , 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 1 , have been drawn, by the proposed program, on the T-P (temperature-pressure) plane, for values of Δt 1 : 0.1, 1, 10, 100, 1000 sec., as shown in Figure 7. It is reasonably assumed that explosion corresponds to ignition delay times Δt 1 less than 10 sec. So, the curve for Δt 1 = 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   plane, by proposed program, for Δt 1 : 0.1, 1, 10, 100, 1000 sec. Explosion limit curve, for ignition delay time Δt 1 = 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.

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 −9 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 2 ] + [O 2 ], gradually diminishes, due to termination reaction by wall destruction, and tends to the final concentration of product [H 2 O], that is to the 2/3 of its initial value, in accordance to the overall stoichiometric reaction 5) It is assumed that Explosion corresponds to ignition delay time Δt 1 < 10 sec. So, the curve on T-P (temperature -pressure) plane, obtained by the proposed program, for Δt 1 = 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].