Gasdynamic equations with the number of gasdynamic variables exceeding the number of summational invariants

. An asymptotic new method of solving the generalized Boltzmann equation previously developed by the authors is presented, and its difference from the traditional Chapman-Enskog method is discussed. The application of this new method is considered by an example of an O 2 /O dissociating mixture flow with separation around a double cone. A comparison with available literature and experimental data is presented. The computed and measured surface heat fluxes are found to be in reasonable agreement. The computed surface heat flux is analyzed through comparisons to numerical results of other researchers and codes.


Introduction
In modeling complex flows with physical and chemical transformations, the gasdynamic equations are most often used for an extended list of variables, not only for those corresponding to summational invariants (SI) of the original Boltzmann equation (BE), but also for a number of other variables.For the case of chemically reacting gases, these are, first of all, chemical components' densities and, in a more detailed description, temperatures of the internal degrees of freedom, translational temperatures of various components (for mixtures with significantly different components' masses), etc. Extension of the list of gasdynamic variables is relevant due to the high computational cost of solving the kinetic equations, even in their simplified form.If the gasdynamic variable does not correspond to the SI, source terms appear in the corresponding gasdynamic equation (reaction rates, momentum and energy exchange rates, etc.), which are either determined from phenomenological considerations or calculated from the BE using equilibrium distribution functions.The traditional method of deriving gasdynamic equations from kinetic equations does not yield correct results from the asymptotic point of view [1].It may even lead to divergency of higher-order expansion terms with respect to the asymptotic expansion parameter in the corresponding source terms.It was shown [1,2] that the use of equilibrium distributions, which are traditionally employed for the source terms' calculation, can lead to significant inaccuracies.This makes the development of a new approach for derivation of gasdynamic equations relevant.
A new asymptotic method for solving the Boltzmann equation was proposed in [2].Initially it was developed for reaction rates' calculation [3][4][5].It was also applied for deriving a homogeneous nucleation rate [6] and rates of momentum and energy exchange in a mixture of heavy and light gases [7].In all cases, the resultant reaction and relaxation rates demonstrate a complex dependence of the source terms on the components' densities (partial pressures).This leads to the need to abandon the concept of rate constants for chemically reacting mixtures, for instance.Here, it should be mentioned that, in studying the nucleation process and deriving the expression for the nucleation rate within the quasi-chemical nucleation model, the resulting corrections to its classical (Frenkel-Zeldovich) expression were small.Now it is clear that this is a consequence not of the method [2] itself, but of the method of solving the corresponding equations for corrections to the equilibrium distribution function, connected with the clusters' internal vibrations.Those corrections were calculated by expansion into series over the complete set of polynomials.As shown in [8] by an example of an exactly solvable model, the series expansion method significantly underestimates the impact of vibrational non-equilibrium; therefore, we expect significant corrections to the nucleation rate obtained by building a new model for clusters' vibrations and applying our method to the quasichemical nucleation model, as it is done for the reaction rates in [5].
Here we present one more flow example, with nonequilibrium diatomic molecule dissociation behind a shock wave.A two-temperature model with quasistationary reaction rates is employed.
The paper is organized as follows.In the next section we briefly present the method leading to new expressions for reaction and relaxation rates, describe the problem statement, and discuss the computation results; the conclusions are drawn in the last section.

Asymptotic method for solving the Boltzmann equation
The starting point for building the asymptotic procedure is the generalized Boltzmann equation where t is the time, r is the spatial coordinate, v is the particle velocity, F is the distribution function, I is the collisional operator describing the chemical reactions,  =     ⁄ , and   and   are the characteristic time scales of the inelastic processes and the gas dynamics, respectively.The convective characteristic time is assumed to be of the order of the gasdynamic time, and therefore, it is sufficient to introduce only one dimensionless parameter  [2].
In our method, the collisional integral part responsible for the chemical reactions is not considered as a small correction (perturbation) of the order of .Additionally to the traditional exact summational invariants (ESI), for which 〈  ,   〉 = 0 (i numerates corresponding invariants, and  numerates the mixture species), approximate summational invariants (ASI),   , are introduced, which are invariants of the collisional operator I within the small parameter : 〈  ,   〉 ≤ , where 〈, ψ〉 = ∑ ψ ,    is the scalar product.Gasdynamic variables are introduced by corresponding ASI averaging: Γ  = 〈  ,   〉.The distribution function is represented as a sum of two terms,   =    + Φ  .The first term is a function of time and spatial coordinates only via the gasdynamic variables    =    (Γ  ), for which the function obtained from the maximum entropy principle is used.
The generalized Boltzmann equation is written as a singularly perturbed system of equations for the gasdynamic, Γ  , and fast, Φ  , variables.The formulated asymptotic method sets an algorithm to determine sequentially the corrections of different orders in  to the quasi-equilibrium distribution    .
As    is not a solution of the equation I(F) = 0, the correction to    occurs at the zeroth order.At the same time, I(   ) is assumed to be of the order of .As a result, the reaction and relaxation rates in the transport equations (called generalized relaxation rates (GRR)) contain a description of non-equilibrium effects in the zeroth-order (Euler) approximation.The resulting expressions for the reaction and relaxation rates can be referred to as quasi-stationary.

Diatomic gas dissociation model
The gasdynamic equations for a dissociating diatomic gas were derived in [5].The main model positions are as follows.A mixture of molecules A2 and atoms A is considered.The corresponding reactions include: A2 + A ↔ 3A. ( A set of ASIs that leads to the traditional set of gasdynamic variables is chosen: the Kronecker symbol as an ASI corresponds to the species number density Γ   ≡   ,  =  2 , .
The momentum p corresponds to Γ  = , where  and  are the mixture mass density and velocity, respectively.
The total energy corresponds to the total energy density of the mixture: Γ  ≡  =  ̃+  2  finalizes the list of gasdynamic variables.
For our specific case, where only one diatomic species is considered, the subscript in the expressions concerning the internal energy characteristics can be omitted.Denoting the molecule number density as   , where  ≡  2 , its vibrational energy density can be written as   =      .Further on, the subscripts indicating a diatomic molecule are omitted where this does not lead to confusion.
The cutoff harmonic oscillator model is employed for the diatomic molecules' vibrational energy    = ℏ  ( + 1 2 ⁄ ), where   is the oscillator frequency and  is the vibrational quantum number.Though it leads to some inaccuracies, it is widely used to describe molecule internal vibrations when complex chemical reactions are studied [9,10].On the other hand, this model yields analytical expressions for final quantities (reaction and relaxation rates), which makes the parametric analysis clear.In this case, instead of the energy density,   = ℏ  ( ̅ + 1 2 ⁄ )  , the mean quantum value density, Γ  =    ̅, can be used as a gasdynamic variable, where  ̅ is the traditional mean quantum value.
For that set of gasdynamic variables, the quasiequilibrium distribution reads Here   (  ) is the statistical weight of the rotational quantum levels,   (,   ) is the partition function,   =   −  is the peculiar species velocity, and    is the rotational energy of the species .
The equations for these variables,   , , , and Γ  , were derived in [5]: Here , where   is the linearized form of the collisional operator , Φ (1) is the first-order correction to the quasi-equilibrium distribution function (further the superscript (1) will be omitted because only the first-order correction will be considered),   and   are the quasi-equilibrium and non-equilibrium parts of the reaction and relaxation rates.Equations derived within the traditional Chapman-Enskog method contain only quasi-equilibrium parts of the corresponding source terms.

Quasi-stationary reaction and relaxation rates
The quasi-stationary GRRs are determined not only by the equilibrium distribution, but by the first-order correction determined by the equation where    () = () + ∑ is its linearized form [2], while   are the ASIs.The summation is performed over all gasdynamic variables (M is their total number), and   are the GRRs.As the corresponding source terms vanish for ESIs, only ASIs contribute to this sum.The presence of additional terms in Kolesnichenko's collisional operator, , reflects the impact of variation of the gasdynamic variables through the terms containing   during collisional processes.
Assuming translational-rotational equilibrium, after integration over the velocity and summation over the rotational quantum number, Eq. ( 9) can be written as where  =   ,   = ∑ ∫ F    is the vibrational distribution function,    is its quasi-equilibrium value, and Ψ  is the corresponding correction which determines the quasi-stationary GRRs in their nonequilibrium parts.The exact solution of Eq. ( 10) was obtained in [5] for the cut-off oscillator model.Within that model, only one-quantum vibrational transitions are allowed with probabilities  ,±1 for VT-transitions and ,∓1 for VV-transitions.Dissociation is assumed to occur only from the highest vibrational level.As a result, the following quasi-stationary expressions for the GRRs were derived: where    and    are the quasi-equilibrium parts of the corresponding GRRs obtained within the traditional Chapman-Enskog method: Here   is the dissociation rate constant of the considered molecule from the level ,   is the corresponding recombination rate constant with the formation of the molecule at the vibrational state ,   is the highest vibrational level,   is the number density of the atoms A, and  10 is the probability of the VTtransition from level 1 to level 0.
There are no spatially inhomogeneous nonequilibrium corrections (proportional to hydrodynamic velocity divergency) within the harmonic oscillator model [5].The spatially homogeneous corrections are ). ( The renormalizing coefficients are . The coefficients   and   are proportional to the ratio of the vibrational and chemical time scales, which increase with temperature and may reach rather large values.That ratio is independent of the Knudsen number; therefore, the renormalizing coefficients cannot be neglected even if the Knudsen number approaches zero.On the other hand, it means that expressions (11) and ( 12) cannot be derived by summing over the corresponding terms of a different order of the traditional Chapman-Enskog expansion.
Expressions (11) and ( 12) transform into each other by means of permutation of the subscripts  and .This symmetry emphasizes the mathematical similarity of these quantities.From the physical point of view, it means the mutual impact of the reaction and relaxation rates.While the reactions' impact on the relaxation rates is well known [11], the opposite effect was obtained within our approach for the first time.Nevertheless, it should be noted that this opposite effect is negligibly small for the considered regimes due to the small values of the   coefficients.

Problem formulation
The governing equations ( 5) -( 8) are solved numerically with Ansys Fluent: system (5) -( 8) is adopted to a form consistent with the Ansys Fluent's form.Equation ( 8) is added as a user defined scalar.GRRs (11) and ( 12) are computed as source terms in user defined functions.A density-based solver is used for computations with the second-order upwind scheme and AUSM for the inviscid flux.Other details can be found in [5].

Computational domain and boundary conditions
An O2/O binary mixture flow around a double cone is simulated numerically under conditions of Run 87 considered in [12].The computational domain is highlighted in Fig. 1 with the grey color.The following boundary conditions are used.Segment AB is the axis.Segment BC is a supersonic flow with u∞ = 4019 m/s, p∞ = 165 Pa, T∞ = 625 K, Tv,∞ = 712 K, and YO2 = 0.9245 (mass fraction).Segment CD is the outlet where all variables are extrapolated from the interior.Segment DA is an isothermal wall at Tw = 300 K with no-slip and zero diffusion flux conditions (non-catalytic wall).The free-stream parameters correspond to the Mach number M∞ = 8.06, and the total enthalpy is h0 = 9.85 MJ/kg.Numerical simulation is performed on a hexa-mesh with 500×300 cells (along and normal to the double cone surface, respectively).The mesh is refined towards the double cone surface in order to resolve the boundary layer.The height of the first cell from wall is 14 µm.Further mesh refinement by two times in both directions showed no influence on the surface heat flux distribution.Therefore, the 500×300 mesh was chosen for further analysis.11) - (12), and quasi-equilibrium, qe, ( 13) -(14), GRRs.

Comparison to other numerical results
In order to assess the fidelity of the presented new model of the quasi-stationary expressions ( 11) -( 12) for the GRRs, we performed a comparison between the surface heat fluxes computed by different models and codes.The first numerical result is taken from [10] and is denoted as "Ninni, 2022."Additionally, we repeated the simulation with the HyCFS code being developed at ITAM.The HyCFS uses the model described in [13] (see model: Park, Cdis = 0.3Ed).
Figure 3 shows all three curves of the surface heat flux.The widely used models (HyCFS and Ninni, 2022) predict the separation point location further downstream compared to the present model.The widely used models also give an earlier rise of the surface heat flux.Strictly speaking, none of the models provides perfect agreement with the experimental data.

Free-stream uncertainty
The experimental facility usually includes a nozzle, which provides a free stream for the considered body.Numerically computed parameters of the free stream (being a jet from the nozzle) also depend on the reaction and relaxation rates.Therefore, the free-stream parameters [12] could be inconsistent with the quasistationary expressions ( 11) -( 12) for the GRRs.Accurate simulation of the nozzle flow under experimental conditions seems to be impossible due to the lack of the nozzle geometry.However, it is still possible to assess approximately how the surface heat flux distribution changes with the free-stream parameters.Here, we use alternatively the free-stream parameters computed in [10].These free-stream parameters are still inconsistent with the presented model; however, a comparison between the original and alternative free-stream parameters allow quantifying the uncertainty in this regard.The alternatively computed free-stream parameters [10]   Figure 4 shows the surface heat flux distributions for the original [12] and alternative [10] free-stream parameters.One can see that the locations of the separation point and heat flux peak differ noticeably.Also, the rise of the heat flux in the reattachment region slightly changed the location.

Conclusion
One more example of the new method for deriving gasdynamic equations is considered.The simulations performed using the quasi-stationary expressions (11) - (12) show that they provide an adequate flow structure around a double cone with a separation bubble.The surface heat flux computed using Eqs. ( 11) - (1)  agrees fairly well with the experimental data.The presented model predicts an earlier location of the separation point and a later surface heat flux rise as compared to the numerical data computed by other researchers and codes.
The work was performed with financial support from the Russian Science Foundation, project 22-11-00080.

Figure 1
Figure 1 also shows the temperature isolines computed using the quasi-stationary expressions (11) -(12) for the GRRs.It predicts a typical flow structure, which includes an oblique shock wave, bow shock, and separation bubble.Numerical simulation is performed on a hexa-mesh with 500×300 cells (along and normal to the double cone surface, respectively).The mesh is refined towards the double cone surface in order to resolve the boundary layer.The height of the first cell from wall is 14 µm.Further mesh refinement by two times in both directions showed no influence on the surface heat flux

Fig. 3 .
Fig. 3. Surface heat flux compared to other numerical results.