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 Δ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 (H2, O2, H, O, HO, H2O). 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 [H2] + [O2], gradually diminishes, due to termination reaction by wall destruction, and tends to the final concentration of the product [H2O], that is to the 2/3 of its initial value, in accordance to the established overall stoichiometric reaction of hydrogen combustion 2H2 + O2 → 2H2O. 2) Time-histories for concentrations of main species H2, O2, H, H2O 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 [H2], [O2]. So, the maximum reaction rate is that of the forward initiation reaction max R = Rif = kif[H2] [O2], where the rate constant kif is very slow. Thus, the first time steplength Δt1 = 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 Δt1 can be considered as ignition delay time. 4) It is assumed that explosion corresponds to ignition delay time Δt1 < 10 sec. So, the curve on T-P plane, obtained by proposed program for Δt1 = 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].