An approximate method for solving the problem of the establishment of chemical equilibrium in the products of explosion of gas mixture

Based on the assumption of the existence of the partial chemical equilibrium in the detonation products, an approximate method for calculating composition of the detonation products is developed. The method uses the assumption of the existence of extremum of Helmholtz free energy for a given density, temperature, and molecular weight of the detonation products mixture. Without significant loss of accuracy to the solution of stiff differential equations, detailed kinetic mechanism can be replaced by one differential equation and a system of algebraic equations. This method is always consistent with the detailed mechanism and can be used separately or in conjunction with the decision of a stiff system, replacing it when bimolecular reactions are near equilibrium.


Introduction
The chemical reactions that occur in the explosion products change their composition and, consequently, thermodynamic parameters. The time evolution of the composition is calculated using a stiff system of chemical kinetics equations. This is the most general approach. At a large number of individual species in the explosion products and a large number of chemical reactions involved, calculation of the mixture composition requires a long computation time. Typically, the time it takes for calculating the mixture composition is much longer than that spent on solving the gas dynamics equations.
It was demonstrated previously [1][2][3][4][5] that the process of establishment of chemical equilibrium in the products of explosion of gas mixtures occurs in several stages. First, partial chemical equilibrium is established, without changes in the total number of molecules in the system. This is because the rapid bimolecular reactions that ensure the establishment of such partial equilibrium occur much faster than the dissociation and recombination reactions, which are responsible for establishing full chemical equilibrium. That is, even in cases where there is no chemical equilibrium, it suffices to know the molar mass of the mixture and two required thermodynamic parameters (for example, temperature and density) to determine the composition of the mixture.
The existence of quasi-equilibrium in exchange reactions was discussed in [6]. The authors of [6] used this assumption to perform approximate calculations of the composition of a hydrogenoxygen mixture. This same assumption underlies an approximate model for calculating the molar mass of gas mixtures at high temperatures [7]. The author of [7,8] proposed an approximate kinetic equation for calculating the molar mass of the mixture and then the thermodynamic parameters.
In [9] we made an effort to speed up directly chemical integrator for some special case. Another general approach have been taken in [10]. This approach does not solve the chemical equations directly; rather, it relies on a special precalculated database capturing the chemical model's input-output relations.
In [9] we developed a method for calculating the changing composition of the explosion products in the case where there is no full chemical equilibrium, but the bimolecular reactions are in equilibrium. A detailed kinetic mechanism is used, but without the need to solve the stiff system of equations. The core of our approach is to replace the complete system of chemical kinetics equations by an equivalent system composed of a differential and a few algebraic equations. What we offer is not an approximate kinetic equation, but an approximate method for solving the complete system of equations of chemical kinetics. This method is always consistent with the detailed mechanism and can be used together with solution of the stiff system, replacing it at the stage when bimolecular reactions are in quasi-equilibrium state.
The proposed method is based on a physically reasonable assumption that equilibrium in the bimolecular reactions is established much faster than the full chemical equilibrium, being applicable when this assumption holds. According to [9] the method includes only one differential equation for calculating the change in the number of molecules per unit mass of the mixture due to the reactions of recombination and dissociation. This equation is used instead of dozens of equations of chemical kinetics. In the present work, we propose to supplement this equation with an algebraic equations for calculating the equilibrium composition at given values of the density, molar mass and internal energy or temperature. The equations are derived based on the characteristic function (entropy or Helmholtz free energy) extremum method.
This paper is based on the assumption that, under certain conditions, the systems of kinetic equations can be solved as follows. At each time step of numerical integration, the change in total number of molecules in a system is calculated using the single differential equation based on the chosen kinetic mechanism. Then, the mixture composition is calculated assuming that the entropy of the mixture reaches its maximum value at a given internal energy, density, and number of molecules per unit volume or that the Helmholtz free energy of the mixture reaches its extremum at a given temperature, density, and number of molecules per unit volume.
The purpose of the work described here was to investigate the possibility of using the method of "quasi-equilibrium" during the flow-filed of detonation after the end of the induction period. Zel'dovich-Neumann-Doering detonation wave model is selected for this study. To verify the proposed method, its predictions were compared with the results of integration of the complete system of chemical kinetics equations.

Equation for calculating the chemical equilibrium composition of the mixture at constant molar mass
According [9] we obtain equations to calculate the mixture composition assuming that the Helmholtz free energy of the mixture reaches its minimum value at a given temperature, density, and number of molecules per unit volume. The Helmholtz free energy of a mixture of ideal gases with known dependence of the heat capacity of the temperature is given by where k is the number of individual compounds in the mixture, M i is the number of moles of the i-th substance per 1 kg of mixture, R is the universal gas constant, T is the temperature of the mixture, v is the specific internal volume of the mixture, S i is the standard entropy, U i is the internal energy of one mole of mixture component i at a given temperature T and P 0 = 10 5 Pa. The amount of each chemical element in the mixture is known, since it is determined by the initial composition of the mixture, remaining unchanged during chemical reactions. For chemical element j we obtain equation where m is the number of different chemical elements in the mixture, n j is the number of kilogramatoms of the j -th element in 1 kg of mixture, n ij is the number of atoms of the j -th element in a molecule of the i -th substance.
If the number of molecules in the system is constant, then we can use equation where M 0 = c/ρ and c is the total number of moles of all mixture components per unit volume.
Following the Lagrange procedure, we write the function where λ j and λ µ are the Lagrange multipliers, which, along with M i , are the independent variables of the Lagrangian. At equilibrium, the derivatives of the Lagrangian with respect to the independent variables must be equal to zero. This condition allows us to obtain k + m + 1 algebraic equations for k + m + 1 unknowns: This system of algebraic equations has no analytical solution, so it is solved approximately by the Newton method described in [11]. The Newton method requires up to 8 iterations to achieve desired accuracy. Really the number of iterations is not important because for gas dynamic calculations we traditionally use precalculated database for "quasi-equilibrium" species concentrations at given internal energy, density and molar mass. To calculate U i , S i , and enthalpy H i at a given temperature, we used polynomials from [12].

Equation for calculating changes in the number of molecules in the system
When the chemical kinetics equations are applied to calculating the composition of the mixture, the following system of equations is used where ω i is the rate of production of species i per unit volume; c j are the molar concentrations of the mixture components; k j is the rate of reaction j calculated according to the law of mass action, χ i is the chemical formula of species i participating in the reaction j ; ν ij ′ and ν ij ′′ are the stoichiometric coefficients of the product or reactant i in reaction j. The rate of production of the total number of moles of all species per unit volume are given by Formally, in (7), the summation is over all reactions. However, bimolecular reactions do not change the number of molecules, so the summation can be restricted to the recombination and dissociation reactions. Now, to calculate how the composition of the explosion products changes, we can use a single differential equation (7) together with algebraic equations (1)-(3), instead of a stiff system of differential equations of chemical kinetics. Differential equation (7) makes it possible to determine the molar mass of the mixture, and then by solving the algebraic system of equations (1)-(3), to calculate the mixture composition. Calculations are performed at each step of numerical solution of differential equation (7).

ZND model and the model of "quasi-equilibrium"
The ZND model is used to describe the structure of the detonation wave. The governing equations are ρ 0 D = ρu, P 0 + ρ 0 D 2 = P + ρu 2 , H 0 + 0.5D 2 = H + 0.5u 2 , where D -shock wave velocity (D is slightly bigger than the CJ theoretical detonation velocity), u is the flow velocity, subscript "0" is used for variables of undisturbed initial mixture.
For the numerical solution of equations (8) and (9) we transform them to a system of three differential equations for the unknown u, T and µ: If the bimolecular reactions are in equilibrium then enthalpy H is the function of T and µ and derivatives in equations (11) and (12) are obtained as where derivatives of M i are the solutions of two systems of linear equations (systems are obtained from equations (1)- (3))

Results of calculations
In the present study we consider a hydrogen-oxygen stoichiometric mixture of initial density ρ 0 =0.496 kg/m 3 and initial temperature T = 298 K. The detonation velocity was D = 3087 m/s. This detonation velocity is slightly bigger than the CJ theoretical detonation velocity. We used detailed kinetic mechanism presented in [13]. This mechanism includes 12 exchange reactions and 7 recombination reactions for 8 species. At first, the calculations were performed by solving the full system of equations of chemical kinetics corresponding to the selected detailed mechanism. This system is There are 8 differential equations for concentrations c i in system (13). The calculations were performed by Gear based solver. The automatic selection of time step was used. Temperature profile is given in figure 1 by line 1 (thin line).
In addition, we made calculations within the framework of the proposed method. The second calculation was performed as follows. Until the moment when the temperature reached 2500 K, the calculation was carried out by the Gear method for system (13) the same as in first case. Proposed solver was not used during this first time interval. This first time interval includes the induction period. During the induction period, the bimolecular reactions are not in equilibrium, so the complete system of kinetics equations should be solved. After this point, the calculation is carried out for system (8), (9) using the method proposed in the present work, a method with a   figure 1 by curve 2 (thick line). As can be seen, at temperatures above 2500 K, curves 1 and 2 merge. Note that, only over a time constituting a few percent of the process duration, the calculation is performed by the Gear method. The rest of the calculation is performed without loss of accuracy by the proposed more efficient calculation method. Mole fractions profiles of H 2 O and H 2 are displayed in figure 2. The lines in figure 2 represent molar fractions profiles of OH and H. It can be seen that good agreement is obtained for the concentrations of H 2 O and H 2 and for intermediate products (OH and H). Moreover, it follows from results represented in figure 3 that the temperature selected for the "switching" on the approximate calculation method corresponds approximately to the point where the maximum of the intermediate product H is reached. This maximum is often used to evaluate the end of the induction period. Thus, the proposed method shows good accuracy just in the region, namely after the end of the induction period, where the method is supposed to be used.

Conclusions
We developed the approximate method to calculate composition and thermodynamic parameters of detonation products based on the assumption of the existence of a partial chemical equilibrium. Without significant loss in accuracy, the respective stiff system of detailed kinetics differential equations can be replaced by a single differential equation and a system of algebraic equations. This method is always consistent with the detailed kinetic mechanism and can be used separately or in conjunction with the stiff system, replacing it when the bimolecular reactions reach the quasi-equilibrium state. The constituent equations of the model were derived and the respective computer code written. The applicability of the model was demonstrated by solving a test problem. It was shown that the proposed model could be used to calculate the characteristics of the detonation products flow-field after the induction period.