Modeling of nitrogen and oxygen gas mixture with a novel diatomic kinetic model

Practical applications involve flows that often have more than one constituent. Therefore, the capability to model a gas mixture flow is important. Extending kinetic model equations of the Bhatnagar–Gross–Krook type from a single-species gas to multi-species gas mixtures presents a number of important challenges. This challenge is further pronounced when diatomic gas mixtures are considered due to the addition of internal energy modes. In this paper, a novel diatomic binary mixture model with separate translational, rotational, and vibrational temperatures is derived. The species drift-velocity and diffusion are considered by introducing separate species velocities and accounting for their relationship. The derivation is detailed as a logical build-up with a multi-step approach from a diatomic model for a single gas, known in the literature. Transport properties are obtained through the Chapman–Enskog type expansion. The diatomic mixture model is numerically evaluated for a gas mixture of nitrogen and oxygen. The model is validated against Monte Carlo results for normal shocks, showing good agreement for species density and temperature profiles. A parametric study demonstrates the variation in flow properties for different Mach numbers, vibrational collision numbers, and concentrations. Interesting results for the mixture properties are shown when the physics of the flow is discussed in greater detail. The effect of the different levels of vibrational excitation in the different species emphasizes the importance of modeling the flow as a mixture. The newly introduced diatomic gas mixture model demonstrates promising computational results for


I. INTRODUCTION
Most practical applications include flows involving gas mixtures, and in the vast majority of cases, these are not noble gases. Monatomic mixture kinetic models represent an active research area. [1][2][3][4][5][6][7][8][9] Diatomic molecules present a further modeling challenge due to the internal degrees of freedom (DoFs). Molecules, unlike atoms, are known to have rotational and vibrational DoFs. These extra degrees of freedom and the corresponding flow-field variables associated with the internal DoFs differ for each species in a diatomic gas mixture, which makes the modeling of a mixture of gases with internal degrees of freedom a complex task. At the same time, diatomic gas mixtures are particularly important for aerospace applications such that in cases without chemistry, a binary mixture of nitrogen and oxygen needs to be considered (when ignoring a small concentration of argon). Relevant engineering problems includ hypersonic rarefied gas flows, chemical reacting flows, and plasmas. 10,11 Other practical flows where mixture modeling techniques are important include gas separation in order to create multiple products or purify a single product. 12,13 Another prominent application is a microchannel flow, 14 involving a variety of industries: medical, food, manufacturing, etc. Clearly, a flowfield comprised of a mixture of gases with internal degrees of freedom is of practical interest for numerous domains.
The Navier-Stokes equations lose their validity in the rarefied regime, and other modeling techniques are required. The Boltzmann equation 15 is a fundamental equation in kinetic modeling, and variety flows can be modeled accurately with it. The crucial properties of the Boltzmann equation include the conservation of properties such as mass, momentum, and energy during collisions, entropy dissipation, and collision equilibria. However, as defined, the Boltzmann equation is valid for monoatomic gases and includes only elastic collisions. Extensions of the equations have been derived to include molecules with internal degrees of freedom. 16 The Wang-Chang-Uhlenbeck (WCU) equations include both mechanical and quantized flow modeling. However, the computational cost associated with a numerical technique using a direct, deterministic discretization approach makes the Boltzmann equation unfeasible for practical applications. Note that the diatomic extension (Wang-Chang-Uhlenbeck equations) is two orders of magnitude more expensive to simulate numerically than the monoatomic version (Boltzmann equation). 17 Therefore, for engineering applications, alternative methods are to be found. The direct simulation Monte Carlo (DSMC) is a statistical method to solve the Boltzmann equation. It is a reliable and widely used modeling approach, particularly well-suited to high Knudsen number flows at high speeds. However, in cases with only moderate rarefaction, flows with a low signal-tonoise ratio, i.e., low velocity flows or mixture flows with trace species gas, DSMC is not an efficient and effective modeling approach. In these cases, deterministic methods for the solution of kinetic models may be more effective and efficient. In the majority of applications, this computational efficiency of deterministic solvers can only be achieved in case the full collision term in the Boltzmann equation is replaced by simpler relaxation-type models. Due to the computational memory overhead of deterministic discretization methods, applications are typically limited to two-dimensional flows. Coupling deterministic methods for kinetic equations to computationally efficient continuum solvers creates an effective approach to the simulation of three-dimensional large-scale flows with the localized areas of rarefaction. Relaxation model equations have been widely used and are constantly being improved, with the complex collision term significantly simplified and reduced relative to the Boltzmann or WCU equations.
The simplest and most widely used simplification of the Boltzmann equation is the Bhatnagar-Gross-Krook (BGK) model 18 for single-species monoatomic gases. We will focus on the extensions of the BGK model to a model for diatomic molecules, which is to include the internal energy of the gas, and we will consider the different levels of details of each model. Then, we will introduce a new binary mixture model for diatomic gases. Most initial models for polyatomic gases include the internal DoFs, where the rotational and vibrational components are not separated but represent a part of the internal energy. This kind of formalization is demonstrated in the work of Morse, 19 where the internal energy is quantized, and the work of Andries et al., 20 where an alternative approach is used and the internal energy is continuously distributed. Note that the relationship between the first polyatomic relaxation model shown by Morse 19 to the more detailed and complex WCU 16 is similar to the BGK model 18 and the Boltzmann equation 15 for the monoatomic gas. The monoatomic BGK model is modified by the Shakhov 21 and Ellipsoidal Statistical (ES) 22 models to have a correct Prandtl number in the continuum limit. Similar extensions for a polyatomic gas with a correct Prandtl number were introduced by the Rykov 23 (defined with a continuously distributed internal energy) and ES 22 polyatomic models (as introduced by Holway with quantum energy levels), which are extensions of the Shakhov model and the monoatomic ES model, respectively. The four diatomic models, with the exception of the Andries et al. model, 20 have a multi-relaxation approach. Their collision term of the right-hand side (RHS) of the governing equation consists of two stages, splitting the process into elastic and inelastic relaxation, while the model of Andries et al. has a single-relaxation approach. The single-relaxation allows for a more simple formulation, but the multi-step approach better represents the slower excitation of the internal modes. All of the models introduce translational and internal temperatures as separate flow properties, since it is well-known that in flows with strong non-equilibrium effects, the single-temperature model is insufficient. 10,11,24 Note that from the four models, only the Morse kinetic model 19 does not have a Prandtl correction. A different approach to derive the model of Andries et al. 20 was detailed by Brull and Schneider,25 leading to the same expressions. More recently, the model by Andries et al. 20 was further generalized with a model by Kosuge et al., 26 extending it from a calorically perfect gas to a thermally perfect gas by introducing temperature-dependent specific heat. The model inherits an ellipsoidal statistical (ES) correction for the Prandtl number and the continuous distribution of the internal energy. A kinetic model, introduced by Wu et al., 27 presents a model with continuous rotational energy and quantum vibrational energy. The authors build upon a model previously presented by two of them 28 that has continuously distributed energy by coupling it with the model of Morse,19 which contributes to the quantization of the vibrational energy. The model has a three-step relaxation process and separate translational, rotational, and vibrational temperatures. This paper also features an initial numerical evaluation of the model and shows good comparison with the results obtained with DSMC for the selected initial conditions. Titarev and Frolova 29 generalized the Rykov 23 and polyatomic ES models 20 to take vibrational energy into account. The two extensions have a three-step relaxation process for the translational, rotational, and vibrational relaxation. The Prandtl number correction is present, but the vibrational heat flux in the continuum limit is not accurately modeled, according to the authors. Nevertheless, the model is comprehensive, has good mathematical properties, and is validated with the DSMC results. Another recent model is presented by Mathiaud and Mieussens. 30 The model has a single relaxation time to a Maxwellian distribution function, which is based on a single equilibrium temperature. The internal DoFs include continuous rotational distribution and discrete vibrational energy. It allows for calorically imperfect gases to be considered. The continuum limit and H-theorem are shown. The simplicity of this model and the detailed description allow us to extend it to a multi-relaxation model for gas mixtures in this article. As mentioned previously, the main novelty and contribution of the present work are the derivation of this new binarymixture kinetic model for diatomic gases, with a particular emphasis on the detailed modeling of vibrational energy relaxation in the different species. The original model 30 will be further discussed in Sec. II.
Modeling the flow as a mixture with separate speciesdependent macroscopic properties in the flow-field is important for practical applications. Even though there is a lot of work in the literature on creating models for monoatomic gas mixtures, such models for polyatomic gases (and more specifically diatomic gases) are a less explored area. The model by Mathiaud and Mieussens 30 is said to be extended to a mixture model in the work of Baranger et al. 31 However, this mixture model has a single Maxwellian distribution for both species (in particular, the nitrogen and AIP Advances ARTICLE scitation.org/journal/adv oxygen mixture is discussed due to their small mass ratio), and no species-specific variables or relationships are modeled. As such, the model is lacking detail, e.g., in describing different levels of vibrational excitation in the gases. A recent model, introduced by Klingenberg, Pirner, and Puppo, 32 is capable of modeling a binary gas mixture with internal degrees of freedom. It demonstrates an extension from a single-species diatomic model, 33 which has a separate rotational relaxation applied to the macroscopic variables rather than the RHS of the governing equation. The two gases in the mixture in Klingenberg et al. 32 have individual species-specific selfand cross-collisions, which is consistent with previous monoatomic models from the same authors. 6,7 An ellipsoidal correction (ES) is also demonstrated in this paper, which corrects the Prandtl number in the continuum limit. The other ES-based models for polyatomic gas mixtures have been detailed with a particular focus on chemically reactive gases. [34][35][36] However, these models do not model the internal energy in detail, e.g., with a multi-relaxation collision term or separate translational, rotational, and vibrational variables. Therefore, even though these models have extended capabilities beyond inert mixtures, they do not possess some of the features of recent models for non-reactive gas mixtures. A model focusing on a mixture of polyatomic and monoatomic gases or two polyatomic gases with different discrete internal energy levels was presented by Bisi and Travaglini. 37 The model is yet to separate the effects of the rotational and vibrational DoFs and to introduce the multi-relaxation approach.  39,40 The model by Rykov has also been extensively tested by Larina and Rykov, 41 Titarev et al., 42 the unified gas-kinetic scheme (UGKS), 43 Colonia et al., 44 etc. These papers presenting more recent kinetic models 26,27,31 have shown the results for the flow past a cylinder or the profile of a normal shock within the same paper. Validation of these models is limited. In contrast to the models created through the extension of previous monoatomic models, a more practical approach is taken by Xu and co-workers [45][46][47] in order to simulate a two or three step multi-relaxation process. A more phenomenological approach is taken rather than a detailed mathematical derivation of a new kinetic model. The benefit of this approach is the simplicity of the modeling procedure and the broad variety of test cases.
There are not many results for the binary mixture of nonreacting diatomic gases in the literature based on kinetic modeling. Tcheremissine et al. 17 showed the detailed results for the profile of a normal shock wave for a binary mixture of oxygen and nitrogen and also for a mixture of nitrogen and argon. The model used for the simulations is the Generalized Boltzmann Equation (GBE), which is identical to the Wang-Chang-Uhlenbeck (WCU) equations.
It is clear from the discussion that there is a need for a detailed kinetic model for diatomic gas mixtures and numerical testing to show the model's applicability. In this paper, a new kinetic model for a binary mixture of diatomic gases is presented. The model is for a non-reactive mixture, while the rotation and the vibration degrees of freedom are considered. The focus on this study is on diatomic molecules at high enough temperatures to consider the rotational DoFs fully excited and the vibrational DoFs partially excited while ignoring the electron DoFs. This means that the number of rotational DoFs is fixed at 2, while the vibrational DoFs are temperature dependent. Note that at the specified temperature domain, the rotational energy is treated as a continuum-like translational energy, while the vibrational energy is quantized. The model features species-specific governing equations and distribution functions. Furthermore, the different rates of rotational and vibrational relaxation are accounted for with a three-step relaxation process of the non-equilibrium distribution function. Each species has a separate translational, rotational, and vibrational temperatures. The model also benefits from having separate species mean-velocities, allowing for accurate representation of the species diffusion effect and the velocity drift. The emphasis in this paper is on the procedure of deriving the model. This derivation is demonstrated as a logical build-up with a multi-step approach from a well-described diatomic model known in the literature. 30 This model 30 is for one species that has rotational and vibrational excitation and no chemical reactions. It models the internal degrees of freedom with one equilibrium temperature and a single-relaxation process. In this paper, the governing equations of the model, which we are using as a starting point, are detailed with emphasis on the collision term on the right-hand size (RHS) of the governing equation and the equilibrium distribution function. Starting from the existing singlespecies, single-temperature model, a two-temperature model with common translational and rotational temperatures is introduced. The vibrational temperature is allowed to deviate from the common translation-rotational temperature. Next, this model is extended to a three-temperature model with separate translational, rotational, and vibrational temperatures. All three models' continuum limit is obtained through the Chapman-Enskog expansion, and the effect of each additional temperature is studied. Finally, we extend the three-temperature model with the capability to describe a mixture of two diatomic gases. Species velocity and diffusion effects are accounted for with a new relaxation parameter, similar to the monoatomic mixtures, derived by Groppi et al. 4 and more recently by Todorova et al. 5 The transport properties in the continuum limit are obtained through the Chapman-Enskog type expansion and compared against the diatomic single-species model that we started from Ref. 30 and the monoatomic mixture models. 4,5 The new diatomic model for the binary mixture is shown to have good mathematical properties and to be consistent with the models in the literature.
Further to the theoretical derivation of the kinetic model, a numerical evaluation and testing are conducted. The model is implemented in the parallel multi-block discrete-velocity solver within the MΦC (Multi-Physics Code) framework. 44,48,49 The profile of a normal shock wave is the simplest form of a flow with high levels of non-equilibrium and is a classic benchmark test case. Starting with the evaluation of the binary mixture model with two gases of equal mass and identical initial conditions, it is numerically verified that ARTICLE scitation.org/journal/adv the two gases are indifferentiable. The normal shock profile is shown for a gas mixture consisting of nitrogen and oxygen. A parametric study with cases at different Mach numbers, vibrational collision coefficients, and species concentrations is presented. A detailed comparison with results, acquired with a direct simulation Monte Carlo (DSMC) solver 50,51 under the same conditions, is presented and used for the evaluation of the mixture model. The physics of the flow and more specifically the effect of the different levels of vibrational excitation and relaxation are studied further. The equilibrium distribution function change through the shock is also shown. In summary, the main contribution of this paper is the theoretical derivation and numerical application and evaluation of a novel diatomic mixture model. Unlike other models, a three-step relaxation process is used for the separate evaluation of the translation, rotation, and vibration. Furthermore, the model has species-specific macroscopic variables, including three temperatures per species, and accounts for the velocity diffusion effect. Another particular feature is the numerical application and validation of the model, which is not demonstrated for other diatomic mixture models in the literature.
The rest of this paper is organized as follows: In Sec. II, the underlying kinetic model for a single species with single temperature 31 is summarized. Furthermore, the extension to two-and three-temperature single-species models is presented. In Sec. III, the single-species three-temperature model is expanded to a threetemperature model for a binary mixture of diatomic gases using an approach previously implemented for the successful development of a monoatomic model to a mixture monoatomic model. 4,5 The compatibility of the model and the additional macroscopic variables required for the complete closure of the model are also shown. Section IV focuses on the transport properties of the diatomic mixture model, while these are shown in Appendixes A-C for the models required for the build-up in Sec. II. In Sec. V, the details for the numerical method, problems setup, and non-dimensionalization are provided. In Sec. VI, all computational results for the profile of a normal shock under varied conditions are shown. This includes the parametric study, DSMC validation, different flow properties, and the velocity distribution function plots. Section VII includes concluding remarks and future work.

II. SINGLE-SPECIES DIATOMIC MODELS
In this section, we will build up a three-temperature model with translational, rotational, and vibrational degrees of freedom and separate temperatures. Note that the emphasis here is on three different single-species models: with one common temperature, with two separate temperatures (common translational and rotational and separate vibrational temperature), and with three temperatures. The goal is to observe the changes each additional temperature provides and examine the transport properties in the continuum limit for all three models.

A. Single-temperature model
This section summarizes the work by Mathiaud, Mieussens, and co-workers on a BGK model for diatomic gases, presented in Refs. 30 and 31. The model is well-defined and proved theoretically to be mathematically consistent. Furthermore, numerical tests for moderately rarefied gas flows (Knudsen Kn = 0.01) demonstrate good agreement in comparison to the results obtained by a compressible-flow Navier-Stokes solver. 31 The BGK model includes rotational and vibrational degrees of freedom, where rotation is considered fully excited and continuous, while vibrational energy is partially excited and has quantized levels. The distribution functions and the governing equation are presented for each quantized level i, but for brevity, the subscript i is omitted, e.g., the distribution function f i is written as f. The model contains a single temperature T and gas velocity u 0 . The governing equation is in the form ∂f ∂t where the collision term on the RHS of the equation is denoted as Q with Q = 1 τ (m eq −f ) and τ being the relaxation time. The equilibrium distribution function is in the form where T is the equilibrium temperature, I is the rotational energy, θv is the characteristic temperature for vibrational excitation, m is the molecular mass, and i is the considered quantized level. The distribution function is defined at each quantized level, which is a very detailed modeling approach. However, the complexity of the model is reduced by summing over all quantized levels i up to infinity, where the Boltzmann distribution of the vibrational states is assumed, leading to the removal of discrete vibrational levels and a reduced form of the distribution function m eq . 30 In practical applications, only a finite number of vibrational quantum states need to be considered, and this number may be different for different diatomic species.
The first order Chapman-Enskog expansion provides an expression for the non-equilibrium distribution function, where the collision parameter τ is equal to τ = ετ. Note that the internal energy for the single-species singletemperature model is in the form where R = kB/m and δ(T) is the number of excited vibrational degrees of freedom, and for the harmonic oscillator, δ(T) is expressed as The harmonic oscillator is one of the simplest and most widely used methods for modeling intermolecular forces for the considered nonreactive flows. Note that if the kinetic models were to be extended to involve chemistry, then the more advanced and realistic anharmonic oscillator should be considered. The number of vibrational degrees of freedom δ(T), as defined in Eq. (5), for low temperatures with ARTICLE scitation.org/journal/adv respect to θv approaches 0, while the high-temperature limit is 2 (the fully excited state). The expression for the internal energy e is always true, as defined with the equilibrium temperature T, as shown in Eq. (4). The continuum limit and the H-theorem are shown in detail by Mathiaud and Mieussens. 30 This is a well-defined single temperature model including rotational and vibrational degrees of freedom. We will build on it to extend the level of accuracy, physical details, and applications.
B. Two temperature model The described single temperature model accounts for variable degrees of freedom but has a distinct disadvantage: the different types of internal energies are all based on one temperature, which is the equilibrium temperature, while many models, e.g., Refs. 20, 23, 26, 27, and 29, introduce different types of temperatures, according to the associated DoFs. We use a similar approach in the extension of the single-temperature kinetic model. Here, we present a two-temperature kinetic model for a diatomic gas. The model contains a single translational and rotational temperature Ttr and a separate vibrational temperature Tv and gas velocity u 0 . The typical BGK relaxation from non-equilibrium to equilibrium state, i.e., f → m eq , now has an extra step. The relaxation is in the form f → m 2 → m eq , where m 2 is an intermediate equilibrium distribution function in which the temperatures Ttr and Tv are separate, while the final equilibrium function meq stays the same and has a common equilibrium temperature T. The two-step relaxation process is not uniform, and the time for the non-equilibrium function f to relax to the first (intermediate) equilibrium function m 2 is τ as before, while it takes longer for the vibration to relax to the translation and effectively to equilibrium, represented by a longer relaxation time Zvτ between m 2 and m eq (Zv > 1). Typical values of the Zv coefficient are in the tens to hundreds and in some cases in the thousands. The governing equation is in the form where the collision integral is ) and the non-equilibrium function f and the equilibrium distribution functions are for each quantized level i, as in the previous model. The intermediate distribution function is in the form while the equilibrium distribution function is in the form Following the expansion of the non-equilibrium distribution function with the CE procedure in powers of ε (with ε being a small parameter of the order of the Knudsen number Kn) and utilizing the assumption that the difference between m eq and m 2 is of order O(ε) and therefore ετ( leads to the expression for the non-equilibrium distribution function up to and including first order terms in ε for the two-temperature model, In comparison to the single-temperature model [Eq. (3)], the zeroth-order expansion is the intermediate equilibrium function m 2 , while also the difference to the Maxwellian-like equilibrium m eq is accounted for with the relaxation coefficient Zv involved.

C. Three-temperature model
Building upon the two-temperature model, a three-temperature model for a diatomic gas is introduced. The next step of this process is to split the common Ttr temperature into a separate translational and rotational temperature while preserving the separate vibrational temperature. This allows for each different type of energy to be represented by the corresponding temperature. It will also assist in solving problems with high levels of non-equilibrium where all three internal degrees of freedom are excited at different levels and differ significantly between each other due to the difference in characteristic relaxation rate for translation, rotation, and vibration.
The new relaxation process is a three-step approach and follows the process f → m 1 → m 2 → m eq , where m 1 is the first intermediate equilibrium function and is a function of all three temperatures, m 2 is the second intermediate equilibrium function and is in the same form as before (two temperatures), and m eq is the final equilibrium function with one temperature. The three-temperature model contains a translationalT, rotational Tr, and vibrational temperature Tv. The gas velocity is u 0 . The governing equation is in the form where the collision term is The three-step relaxation process is not uniform, and the delay of rotational to translational equilibration is via the rotational collision number Zr, while to vibration is via the vibrational collision number Zv (Zv > Zr > 1). These coefficients are temperature dependent and the ratio of the two coefficients of rotation and vibration with respect to τ is typically 1:5:100 17 , with Zv at least one order of magnitude higher than the other two (with values of hundreds and thousands). Here, for simplicity, Zr and Zv are taken as constants, i.e., assumed to be independent of temperatures, and therefore fixed throughout the considered flow domain.
The first intermediate distribution function is in the form

ARTICLE scitation.org/journal/adv
The second intermediate distribution function is in the form while the equilibrium distribution function is in the form Therefore, the expansion of the non-equilibrium distribution function in powers of ε where m 1 is the zeroth-order expansion and the same property as in Eq. (9) for the difference between two equilibrium functions leads to For the three-temperature model, both collision numbers are involved and we note that they are present in the first order expansion, which means that they will directly influence the transport properties [which are of order O(ε)].
To visualize the main difference in the three models, Fig. 1 demonstrates the different temperatures for each model and their relaxation process. The continuum limit for the three models is summarized in Appendixes A-C. Splitting the collision term into three different relaxation processes is a widely used approach, and the above kinetic model is very similar to the UGKS three-temperature model 46 in terms of the definition of the equilibrium functions and the governing equation formulation. Furthermore, if vibrational excitation is not considered, the model will have only two relaxation processes (without the Zv), which is similar to the UGKS with translational and rotational non-equilibrium, presented in Ref. 43. Another two-step approach is demonstrated by the Rykov model, 23 but the formulation is classified as a heat-flux correction for diatomic gas as an alternative to the Shakhov correction for monoatomic gas. This shows that this multistep build-up from the single temperature model to the three temperature model is consistent with models previously published in the literature.

III. DIATOMIC MIXTURE WITH THREE TEMPERATURES
After building up a diatomic model from a single-temperature model, available in the literature, 30,31 to a three-temperature model, which can model separately translation, rotation, and vibration, we focus our attention to expand this model to a diatomic mixture model. Note that we have defined the model for a binary mixture, similar to the monoatomic mixture models. 49 In this section, a diatomic mixture model with separate species velocities is presented. A species target velocity is also introduced similar to the monoatomic mixture models in previous work. 4,49 Furthermore, the model features three temperatures per species and has a common translational target temperatureT (used in the equilibrium distribution function), consistent with the monoatomic mixture model.
The governing equation for a diatomic mixture is per species s, while the rotational and vibrational collision numbers Zr and Zv are kept constant for all species. The characteristic temperature θv s is also per species s, which allows for specific gases with different levels of vibrational excitation to be considered. This is one of the key benefits of this model. The governing equation per species for the diatomic kinetic model is where f s is the non-equilibrium function, ms 1 is the initial (or first) intermediate equilibrium function, ms 2 is the second intermediate equilibrium, and m eq s is the final equilibrium function, resembling a Maxwellian distribution. All distribution functions follow the same notation as before but are per species s. Similarly, the governing equation and the distribution functions describing it are quantized (for the vibrational energy levels) as in the single-species multitemperature approach. Note that, in general, for different diatomic gases, different numbers of quantized levels need to be considered.
The first equilibrium function ms 1 that the non-equilibrium function f s relaxes to is in the form where the target translational temperature is common for both species and isT. Furthermore, the target species velocity is u s . Both of these variables are consistent with the monoatomic binary mixture model by Groppi et al. 4 and the two models presented in Ref. 5. The relationships between the target temperatureT and the species translational temperatures Ttran s and between the target species velocities u (g) s and the species velocities us and gas velocity u 0 will be discussed in detail. They are important for the compatibility of the model, but also to model the diffusion behavior of the two gases in the mixture. More terms specific for this equilibrium function are the species rotational temperatures Tr s and the species vibrational temperatures Tv s . Note that since this is the first or initial equilibrium function, the gas is furthest away from equilibrium in AIP Advances ARTICLE scitation.org/journal/adv the considered relaxation process. This is the reason all temperatures differ from each other.
In the next relaxation step, the first intermediate equilibrium function ms 1 relaxes to the second intermediate equilibrium function ms 2 with relaxation time Zr times longer than the initial relaxation. At the stage of the relaxation process represented by this equilibrium function, the species translational and species rotational temperatures should be in equilibrium leading to a species translationalrotational temperatures Ttr s . However, in the time of having translational and rotational temperatures equalize, the species translational temperature will be common for both species and similarly the species rotational temperature will be common as well or Tr 1 = Tr 2 . Therefore, having species translation-rotation temperature Ttr s is not necessary and a common temperature Ttr is used, whereT = Tr 1 = Tr 2 = Ttr. The second intermediate equilibrium function becomes where two species vibrational temperatures Tv s represent the fact that the vibrational temperatures have equilibrated yet with the other temperatures. Important to notice in Eq. (18) is the velocity, which has relaxed from species target velocity u (g) s to the common mixture velocity u 0 of the gas. There are two possible relaxation routes for the velocity, which is not directly related to the Zr and Zv collision numbers. It is considered that the velocity equilibration occurs with the quicker initial relaxation rather than the slower relaxation. For this reason, the ms 2 distribution function is defined with the gas velocity u 0 .
Finally, the last equilibrium distribution function is defined with the equilibrium temperature T and the gas mixture velocity u 0 as At this last stage, all temperatures are in equilibrium, which makes this distribution function to resemble a species Maxwellian distribution function for a diatomic gas. The relaxation process of the distribution function is as fs → ms 1 → ms 2 → m eq s . To visualize the temperature relaxation, the temperature reduction stages are shown for the mixture model in Fig. 2, similar to the comparison between the single-, two-, and three-temperature single-species diatomic models before. The red temperatures and dashed lines are stages in the temperature equalization, which are omitted in the modeling process. Despite being omitted, species-specific translational temperatures can still be obtained from taking moments of the non-equilibrium distribution function f s. However, these temperatures play no role in the definition of the current model. These red steps are justified through the representation of the model in this section. The equilibration process begins with five temperatures for the two-species The diatomic mixture model can be easily shown to reduce to the single-species single-temperature model 30 for Zr = Zv = 1 and then equal species properties and more specifically ρs = ρ, ms = m, and θv s = θv. This demonstrates the consistency of the diatomic mixture model.
In this section, the governing equations of a new kinetic model for a binary diatomic gas mixture were presented. Furthermore, modeling decisions were discussed and justified. The relaxation process was detailed. Sections III A-III D focus on closing the model by deriving the expressions for the target velocity and species temperatures, demonstrating the compatibility of the model. In this section, the compatibility property of the model is demonstrated. It is shown that taking moments with respect to the collision invariants of the collision integral on the RHS of the governing equation equals zero. This property was also used in monoatomic mixture models 4,5 to provide the expressions for the target species velocity u (g) s and target temperatureT. Similarly, these expressions will be derived for the diatomic mixture model together with the expressions linking the temperatures received by taking moments of the distribution function and the translational-rotational temperature Ttr and the equilibrium temperature T.

A. Mass compatibility
The collision term on the RHS of the governing equation is , and even though the model is for a mixture, the mass is conserved per species s, since the gas is non-reacting. This will differ for both momentum and energy compatibility equations. The proof of the mass compatibility is trivial, since the first moment of all functions is the species density ρs and does not impose new conditions between the macroscopic variables.

B. Momentum compatibility
The momentum of the diatomic mixture model is conserved for the whole mixture and not per species s. Therefore, the sum of the moments of the collision integral Qs with respect to the collision invariant u is assessed as and leads to a requirement for the target velocity u (g) s , C. Definition of the u (g) s velocity So far, we have described the velocity u (g) s as an equilibrium target velocity for each species representing the velocity in the initial equilibrium distribution function ms 1 . However, taking moments of the non-equilibrium distribution function f s will lead to the species velocity us. The gas mixture velocity u 0 can be found using the species velocity us and Eq. (21). Up until this point, we do not have a mathematical expression for the species target velocity u (g) s . Further to Eq. (21), a second velocity relationship is required. This is provided by the same constraint as for the model of Groppi et al. 4 for the species-velocity equalization and includes the introduction of the relaxation parameter η, where Q 1 and Q 2 are the collision terms for species 1 and 2, respectively. This leads to A necessary condition is that Zr ≠ 1 and Zr cannot be speciesspecific for this derivation, while Zv can be. In comparison to the monoatomic mixture models, 4,49 where the expression for the velocity has a dependency on the relaxation parameters η and 1/ν = τ, for the diatomic mixture model, a further dependency on the collision number Zr is present. For typical values of Zr in the interval 1 < Zr ≤ 5, the ratio Zr/(Zr − 1) will always be bigger than 1 and less than 2. It can therefore be argued that the collision number dependency leads to a small weighted advantage toward the species velocity u s rather than u 0 in the diatomic mixture model. In addition, relatively to the monoatomic models, the same coefficient ητ will lead to the value of u (g) s closer to u s . Now, all velocities associated with the presented diatomic mixture model are specified and can be computed.

D. Energy compatibility and temperature definitions
Equivalent to the momentum, the energy of the diatomic mixture is preserved for the sum of all species and not individually. Furthermore, it is preserved for the sum of all types of energies: translational and internal. The sum of the moments of the collision integral Qs with respect to the collision invariants for the translational, rotational, and vibrational energy is where |u| 2 = u 2 + v 2 + w 2 . The only way to guarantee this expression is equal to 0 for different choices of Zr and Zv is if each individual relaxation expression with a different collision number in front of it is 0, leading to three separate conditions to ensure the energy compatibility of the mixture model. After a long but standard integration, we reach the following requirements for theT, Ttr, and T, defined from the species temperatures Ttran s , Tr s , and Tv s , which we find as moments of the non-equilibrium distribution functions: To guarantee positive temperatures, the limit of possible values for the coefficient ητ = η/ν is found, similar to the monoatomic mixture models. 5 For positive macroscopic values for the temperatures, the expressions in Eqs. (25) and (27)

IV. TRANSPORT PROPERTIES FOR THE GAS MIXTURE
The Chapman-Enskog 15 expansion is widely used in the kinetic gas theory, and it was established and derived specifically for monoatomic gases, while it is also often used for diatomic models. Monoatomic kinetic models utilize the expansion to demonstrate the continuum limit and show the corresponding closure of the Navier-Stokes equations. The derived continuum limit allows for the transport coefficients to be evaluated and with the help of relaxation coefficients to be correctly fitted. The transport coefficients for the diatomic mixture model inherit the properties of the single-species diatomic model: viscosity, bulk viscosity, and thermal conductivity due to the internal degrees of freedom. 52 Additional properties arise that are mixture related such as the ordinary diffusion, Soret and Dufour effects. Comparing the continuum equations of a monoatomic mixture model, single-species single-, two-, and three-temperature diatomic models, and diatomic mixture will highlight these properties. The focus of this section is on the novel diatomic mixture model, while the properties and Navier-Stokes equations of models we used to build-up to it are summarized in Appendixes A-C.

A. Diffusion
In comparison to the three single-species diatomic models, where no mass diffusion behavior is observed, a mass flux term is present in the continuum limit of the mixture model, caused by the gradients of the concentration, pressure, and temperature. The mass flux term is in the form which is identical to the mass flux of the monoatomic mixture models. The diffusion effect of the monoatomic mixture and the diatomic mixture is identical. 53 In the monoatomic mixture model, 49 as well as here, we ignore the temperature effect (associated with the Soret coefficient), leading our focus to the ordinary diffusion. Therefore, the binary diffusion coefficient Dsr = Drs, where r ∈ {1, 2} and r ≠ s, is such that

B. Viscosity
For the diatomic mixture model, two viscosity coefficients are present: the ordinary viscosity and the bulk viscosity. The momentum equation can be rewritten as where the ordinary viscosity μ for the diatomic gas mixture considered here is identical to that for the monoatomic model and is expressed as Here, we should take a moment to appreciate that the value of ητ is defined from the relationship of the ordinary diffusion and ordinary viscosity coefficients, formulated as 53 where A * 12 is a non-dimensional coefficient, defined by the ratio of collision integrals, and is, in general, a function of the gas temperature and the force law between molecules. This relationship is identical to the monoatomic gas mixture, 49 as are the ordinary diffusion and viscosity coefficients. Therefore, the expression for ητ is the same and equal to η/ν. Furthermore, the value of A * 12 = 1.10 is found appropriate, 54 This coefficient will vary the most in regions of strong nonequilibrium and for bigger mass ratios.
As expected, we observe changes in the momentum equation. If we compare the monoatomic and diatomic momentum equations, an additional term appears. This additional term is known as the bulk viscosity. The bulk viscosity for a diatomic gas is due to the presence of internal degrees of freedom. 10,52,53,55 The difference of relaxation times between translational, rotational, and vibrational degrees of freedom contributes to the definition of this transport term. Due to this difference, the work done by the pressure on the gas affects the mechanical (translational) energy faster than the thermal (by inelastic collisions) energy, which then gives rise to the bulk viscosity. 15,56 The momentum equation in the x-direction for the monoatomic gas is shown in Appendixes A-C. Splitting the viscous term in the diatomic momentum equation into a normal stress that matches the monoatomic gas momentum and examining the remaining terms lead to the bulk viscosity. This is possible since we know that the ordinary viscosity does not differ for a monoatomic and diatomic gas. 53 After identifying the ordinary viscosity coefficient, the contribution to the normal stress due to the bulk viscosity is We observe an expression similar to the single-species threetemperature diatomic model (see Appendixes A-C). Again, three temperatures (T, Ttr, T) contribute directly to the bulk viscosity expression, and the relaxation phenomenon is presented by the collision numbers Zr and Zv. The difference with the singlespecies expression is only in the last term, which is multiplying the derivative of the mixture velocity u 0 . There, a mass averaged contribution of the different degrees of freedom is present. In the limiting case where s = 1, the normal stress due to the bulk viscosity for the mixture will reduce to the single-species expression (in Appendixes A-C).

C. Heat flux
In the energy equation for the gas mixture (see Appendixes A-C), the order O(ε) expression in the spatial derivative consists of three terms: the work done by the two types of viscosities (ordinary and bulk) and the heat flux, which includes the thermal conductivity and the Dufour effect. The heat flux is expressed as where the term associated with the temperature derivative is the thermal conductivity contribution and the diffusion effect (the term multiplied by the velocity difference) is the Dufour effect.

ARTICLE scitation.org/journal/adv
The Dufour effect qD is where contributions from each species temperatures for the rotation Tr s , the vibration Tv s , and the common target translational temper-atureT are present, while the velocity contribution sums to 0 for the mixture. This expression of the Dufour effects resembles the monoatomic expression, but the additional internal degrees of freedom of the diatomic model give rise to additional terms based on these extra DoFs. Note that for the same species mass, the difference between us and u 0 vanishes, and therefore, the Dufour effect disappears. The other effect in the heat flux expression is the thermal conductivity. The individual contributions of the separate translational, rotational, and vibrational DoFs to the thermal conductivity are These expression are similar to the thermal conductivity for a singlespecies diatomic gas and will, in fact, reduce to it for ms = m, ρs = ρ, and identical characteristic temperature θv s = θv. No Prandtl correction is used in the current model. This is a drawback of the current kinetic model and should be addressed in the future.

V. NUMERICAL METHOD
The diatomic mixture model is tested and evaluated for the profile of a normal shock wave. The computational setup and the problems are detailed in this section.

A. Discrete velocity method
The results of the diatomic kinetic model are acquired using the discrete-velocity method within the Multi-Physics Code (MΦC) in-house framework. 48,49 The method is chosen for its simplicity and accuracy. In the current computational setup, all test cases are considered as two-dimensional (2D) problems. This choice for 2D was based on the planned applications to more complex problems (flow over cylinders, flat plates, and ramps). For the current one-dimensional (1D) shock problems, the extra cost was deemed acceptable, since it avoids the development of a separate 1D version of the solver. Code parallelization is required for all cases. The number of uniformly spaced discrete velocities in the velocity grid is 48 × 48 or 2304 for all test cases. Velocity bounds are between umin = −8 and umax = 12 (i.e., 8 and 12 times the most probable molecular velocity at reference temperature), leading to a velocitystep Δu = 0.417. The local time-step Δt, which is dependent on the local cell size, uses a CFL-based criterion (CFL = 0.5) and is limited by the maximum convection speed. The smallest step is at the shock, where the cells are clustered and Δt = 0.004, which is always smaller than the mean collision time (non-dimensional relaxation time is in the range of 0.1-1). The evaluation of the distribution functions and the fluxes is performed for each velocity in every cell. The macroscopic variables are obtained from taking the corresponding moments of the distribution functions, evaluated using the trapezoidal rule. The finite-volume method used here was derived for multi-block structured, curvilinear meshes with cell-centered data. The second-order accuracy in space is used with a min-mod limiter to suppress oscillations. For the time integration, a second-order accurate Total Variation Diminishing (TVD) time marching method 57,58 is used. This explicit timemarching approach is well suited for the evaluation of the novel kinetic models for the considered supersonic flows. However, for lower-Mach number flows, a change to implicit time-marching would create a more efficient time integration method for the discretized kinetic equations. In terms of relative computational costs of the DVM simulations and the DSMC simulations (detailed in Sec. VI B), a number of observations can be made. First, the considered test case involving the steady flow through a normal shock wave does not constitute a test case where a direct, deterministic method (DVM) has an advantage over DSMC. In fact, the computational resources for the DVM simulation were generally at least one order of magnitude larger than for DSMC, in terms of CPU time and memory. A second important observation is that the DVM implementation used in the present was not optimized in terms of computational efficiency. Instead, the discretization methods were selected on the basis of robustness and relative ease of implementation.

B. Problems setup
The numerical evaluation starts with the study of the profile of a normal shock wave, since it is one of the simplest problems with high levels of non-equilibrium. The normal shock simulations are all performed on a grid with 1152 cells, with x ∈ [−75L ref , 75L ref ]. Note that the domain is selected to avoid the disturbances, downstream of the shock during the simulation, to reflect on the downstream boundary, and to affect the shock structure. The grid is two-dimensional with four blocks and a non-uniform physical space, which is optimized for the smallest in size cells to gather near the shock, where Δx = 0.125.
The considered diatomic gases are nitrogen (N 2 ) and oxygen (O 2 ) with a concentration ratio similar to air for most cases: 80% and 20%, respectively. The normal shock setup is initialized with two Maxwellian equilibrium distributions on each side of the shock. The equilibrium functions are defined based on the macroscopic variables provided by the modified Rankine-Hugoniot conditions for a diatomic gas with rotational and vibrational excitation. In the test cases studied, the pre-shock values are evaluated at a temperature T∞ = 288.15 K for which the level of vibrational excitation is considered negligible and the number of excited vibrational DoFs for each species δs(Tv s ) = 0, while the rotation is taken as fully excited for both gases. Therefore, for the pre-shock values, the ratio of specific heats is taken as γ = 1.4, and the speed of sound is based on the mass mixture mmix = ρ/n. An iterative method is applied to evaluate the jump conditions, described in detail in Ref. 10, for each test case.
The validation of the numerical results for the normal shock is performed by comparing with the results acquired by a DSMC solver with the same macroscopic conditions. There are not many validation data for gas mixtures available in the literature with which ARTICLE scitation.org/journal/adv the accuracy of the model itself can be evaluated and the profile of the shock is a good starting point.

C. Dimensional reduction
The dimensional reduction is a standard and very useful procedure for solvers based on kinetic modeling. The distribution functions are reduced with respect to the rotation and the vibration as follows: leading to a set of governing equations per species s that need to be solved together, where the corresponding equilibrium functions are M tran Further to this reduction and similar to the monoatomic solver for gas mixtures, these functions can be further reduced with respect to the w-velocity, since for the two-dimensional setup, the mean velocity w 0 = 0. This is done by reducing the distribution functions to a translational Gs component and a thermal energy Hs component, where all distribution functions are reduced translationally, while the thermal component is required only for the kinetic part of the energy equation associated with the first of the set of governing equations in Eq. (39). Therefore, the solver uses in total eight governing equations: four for each species. The macroscopic variables are then found by taking moments of the reduced non-equilibrium distribution functions: G tran s , G rot s , G vib s , and H tran s : The macroscopic variables are found as moments of the reduced distribution functions, similarly to the monoatomic solver, finding ρs, us, Ttran s , Tr s , and Tv s from which u 0 , u where the expression for Tv s contains δs(Tv s ) and an iterative procedure is needed to obtain Tv s .

D. Dimensionless form
The non-dimensionalization and reference values are summarized here. For the diatomic mixture model including vibrational excitation, the gases in the mixture need to defined in terms of their respective characteristic temperature, required for the calculation of the number of vibrational degrees of freedom. The focus is on air, and therefore, the mixture consists of nitrogen with characteristic temperature θv 1 = 3371K and oxygen with θv 2 = 2256K. The mass ratio β = m 2 /m 1 = 1.14.
Note that for a hard sphere model, the power ω in the power-law for viscosity used here becomes ω = 0.5 and the expression for the viscosity is where the mass of the mixture mmix = ρ/n and the ratio between the kinetic reference length L ref and the mean free-path λ ref is now easily found. It is an important ratio when comparing with the results from DSMC.

E. Normalized values
Similar to the results from the monoatomic simulations, the macroscopic quantities presented for the shock wave are normalized following Kosuge's approach 59 from both the kinetic model and DSMC, where y is the macroscopic variable, Y− is the pre-shock value of y, and Y+ is post-shock, Note that the tilde is omitted from here on for simplicity. In addition, the origin of the plots (X = 0) is defined at the location where the total number density is exactly half of the sum of the pre-shock and post-shock values.

A. Single species normal shock wave
The numerical evaluation of the gas mixture kinetic model begins with the evaluation of the profile of a normal shock wave. Before conducting a thorough test with a gas mixture, however, we will first focus briefly on inspecting numerically the collapse of the mixture model to a single-species diatomic model. We focus solely on evaluating the model for a single-species but preserving the three-temperatures and the three-step relaxation. The results will demonstrate the separate rotational and vibrational relaxation processes. We show that under the same initial conditions for two identical species and with 50% concentration for each species, the macroscopic variables for both gases are identical throughout the shock wave. Note that further to the identical initial conditions, concentration and species mass, the characteristic vibrational temperature for both species has been set for one type of gas only, i.e., the same value is used since two identical gases are considered. In this case, we have selected the characteristic vibrational temperature for nitrogen. To initialize the problem, further to the concentration and the Mach number, the free-stream temperature is defined dimensionally to evaluate the number of vibrationally excited DoFs. The translational and rotational DoFs are considered fully excited, which is a reasonable assumption for the considered temperature T∞ = 288.15 K. Furthermore, the values of the collision coefficients Zr and Zv need to be specified before each simulation. As described during the development of the kinetic model, the two coefficients are considered constants and are not temperature-dependent, which is a simplification to be addressed in future work. The free-stream conditions are as follows: Mach number M∞ = 3, Zr = 5, Zv = 25, and free-stream temperature T∞ = 288.15 K; the hard-sphere model is employed for both species so that ω is constant and equal for both species at ω = 0.5 for two N 2 gases with a characteristic temperature θv = 3371 K. What we expect to find is all species macroscopic variables equal to each other, i.e., n 1 = n 2 , u 1 = u 2 , Ttran 1 = Ttran 2 , Tr 1 = Tr 2 , and Tv 1 = Tv 2 , since there is no process to separate the two species. In Figs. 3(a) and 3(b), these properties are inspected, where the solid colored lines are used to denote species 1 properties and the black dashed lines with the corresponding colored symbol indicate the results for species 2. In (a), the species number density is in red, while the species velocity is in green. In (b), the translational temperature is in red, the rotational temperature is in blue, and the vibrational temperature in green. This color convention for the species temperatures is kept throughout this section. It is clear that the two species have identical shock profiles and have collapsed into a single species, which was the desired result. This test shows the consistency of the model numerically. Note, in this section, we focus on showing the property of the model to reduce to the same species, while the physical phenomenon observed in the results will be discussed in detail in the section concerning the gas mixture.

B. Gas mixture normal shock wave-comparison with DSMC
In this section, the diatomic kinetic model is tested for a binary gas mixture. In comparison to Sec. VI A, the effect of diffusion between species is present and the species temperatures can deviate from each other. Note that for the practical application of air, the mass ratio between nitrogen and oxygen is quite small (m 2 /m 1 = 1.14) and is not significant enough to create a big difference between the species in terms of species-mean velocities. The non-equilibrium caused by the species mass difference, which we observed in the monoatomic mixture, 49 is not a dominant effect here. Of course, for different gases, e.g., a mixture of hydrogen and nitrogen with a mass ratio of 14, the diffusion effect will become very important. Another case where the species diffusion becomes important is when the mixture of a diatomic and monoatomic gas is considered, which will create in most cases significant mass ratios. For a chemically reacting air, a five-species gas mixture is often considered with elements N 2 , O 2 , NO, N, and O, and therefore, a mass ratio of 2 is present. In this case also, only some of the gases have rotational and vibrational degrees of freedom and temperatures, while the atoms will have all of their energy stored in the kinetic and translational components. Example application is the atmosphere of different planets beyond the Earth, e.g., Saturn and Jupiter with mostly hydrogen and helium constituents. Therefore, there are a number of cases in which the diffusion of the gas mixture is significant. At present, the new model can model the binary mixture of diatomic gases only. However, when modeling air as a gas mixture, the different levels of vibrational excitation for each species has the greatest impact on the flow-field. This is due to the different characteristic vibrational temperatures for nitrogen and oxygen. It leads to different expressions for each species equilibrium function and separation between the species in terms of excitation of internal degrees of freedom. We will focus our attention on demonstrating this difference through the results for the species translational, rotational, and vibrational temperatures. The pre-shock conditions are evaluated at a temperature T∞ for which the vibrational excitation is negligible, i.e., T∞ = 288.15 K, while in comparison, the characteristic temperatures for the two gases are θv 1 = 3371K and θv 2 = 2256K. Throughout the shock, however, the number of excited vibrational DoFs, which is a non-integer value, changes with the rise in temperature. The normal shock, reducing the speed to subsonic flow behind it, leads to a jump in the internal energy of the flow, represented by the jump in temperature. Since vibration is present, it absorbs part of the internal energy and the temperature is not as high as for the perfect gas jump. This leads to an increase in the density (and density jump) behind the shock. 11 The information detailed in Table I describes the initial conditions for the computation of a normal shock for a diatomic gas mixture, including the Mach number, the concentration ratio of the two gases, and the collision coefficients for rotation and vibration. All other variables are not changed through the simulations, e.g., the free-stream temperature T∞ = 288.15 K, the two gas mixture components N 2 and O 2 , and their characteristic temperatures. The results from the kinetic model are compared and validated against DSMC simulations under the same conditions. The DSMC solver used is dsmcFoam+. 51 The hard-sphere potential is used for simplicity with ω = 0.5 in the power law used for viscosity as in the kinetic model. The number of DSMC particles in the steady state is 2.7 × 10 6 for M∞ = 3 and increases to 3.6 × 10 6 for the higher Mach number M∞ = 5, with a minimum of 25 particles in each cell. The number of samples in the steady state is 35 000, and the spatial grid contains 40 000 cells. Each point in the results is an average of 80 cells, bringing the resolution to 500 evenly spaced bins for all simulations. The smallest mean free path and the mean collision time are found behind the Mach M∞ = 5 shock. The cells are at least 36 times smaller than the mean free path in this region and 50 times smaller than the mean collision time, i.e., all of the DSMC simulations have cells much smaller than the mean free path and time steps much smaller than the mean collision time in the entire domain, which meets good DSMC practices.
The validation with DSMC allows for a very accurate assessment of the validity of the kinetic model. This is due to the fact that further to initializing the flow with the same initial conditions, including the free-stream temperature, we can also set up constant and identical collision coefficients Zr and Zv for both solvers. Currently, the available results in the literature do not allow for the same The jump in the species number density and the mixture number density is shown in Fig. 4. The results from DSMC are presented with elements, squares for species 1, triangles for species 2, and circles for the mixture value, while the kinetic model has red lines for species 1, blue for species 2, and green for the mixture values. The number densities are not normalized. It is easier to inspect the results without the normalization due to the small mass ratio between species and consequent small variation of the normalized results. This also allows us to appreciate the different jump in density for initial velocity or concentration variations. Three cases are shown: M∞ = 3, nitrogen to oxygen concentration 80:20 denoted with n 1 /n = 0.8 in Fig. 4(a), M∞ = 5, n 1 /n = 0.8 in Fig. 4(b), and M∞ = 5, n 1 /n = 0.5 in Fig. 4(c). Furthermore, only results for the vibrational collision coefficient Zv = 25 are shown. Note that the different coefficients will not lead to a different jump condition, and even though it influences the vibrational degrees of freedom and temperatures through the internal part of the shock wave, it has a small effect on the number density. What we observe when we plot the same results for Zv = 100 is almost identical number densities, with the profile of the shock slightly less steep. Note that even though the mass ratio between the two species is very small, the nitrogen reacts first to the shock wave, since it is the lighter gas. We observe that the number density n 1 (in red) is the first to start rising through the shock, which is consistent with the monoatomic mixture model results.
The jump for the same Mach number (M∞ = 5) but with different concentrations is of similar value but is not identical as we see when comparing Figs. 4(b) and 4(c). Having less nitrogen present in (c) leads to a higher jump in number density, since more of the  internal energy is absorbed by the easier to excite the vibrational mode for oxygen, which comprises 50% of the gas here. This leads to a lower jump in mixture temperature and higher number density jump, consistent with the discussion earlier in this section and with the literature. 11 Overall, we observe excellent agreement between the results from the diatomic kinetic model and DSMC for all cases. The number density profiles are mostly identical from the two solvers with only a small difference observed at the second part of the shock. We believe that this is caused by the small variation between vibrational temperatures between the results acquired by DSMC and the kinetic model, discussed further when the species temperatures are plotted. Next, the species translational, rotational, and vibrational temperatures are shown and compared with DSMC for all five cases defined in Table I are with solid lines for species 1 (nitrogen) and dashed lines for species 2 (oxygen), while the results acquired with DSMC are presented with symbols and more specifically squares for species 1 and triangles for species 2. All temperatures are normalized, i.e., leading to a value of 0 and 1 under pre-shock and post-shock conditions, respectively. We observe a good agreement between the results from the kinetic model and DSMC for the structure of the shock. Overall, the profile of each temperature jump is steeper when the Mach number is increased, and therefore, the number of mean free paths for the temperatures to equilibrate post-shock is bigger for the smaller Mach number M∞ = 3, similarly to the monoatomic gas mixture. The translational temperatures for both species have longer upstream tails for the kinetic model in comparison to DSMC. This is a classical characteristic of all BGK-based models for both monoatomic and diatomic gases. Another known phenomenon that can also be seen here is the overshoot of the translational temperature in all figures, similar to the results from the monoatomic gas mixture. The rotational temperatures lag behind the translational and are excited more slowly with respect to the translational temperatures. For all cases, the collisional coefficient Zr = 5, and therefore, the rotation's lag is identical for each case. The same is not true for the vibrational temperatures. The change from Zv = 25 in (a) to Zv = 100 in (b) in Figs. 5 and 6 is very significant. The slope of the vibrational temperatures for the bigger Zv is much more gradual in comparison to the smaller value. The change between the two Zv values leads to ∼15 mean free paths (where λ ref is evaluated under free-stream conditions) difference for the vibrational temperature to reach the equilibrium post-shock temperature for the M∞ = 5 case and ∼20 mean free paths for the M∞ = 3 case. Moreover, the slope of the rotational temperatures is steeper than the vibrational temperatures, while the translational temperatures have the steepest slope. All of these phenomena are consistent with non-equilibrium thermodynamics and physical observations. However, we note that both the kinetic model and the results for DSMC show that upstream of the shock, the vibrational temperatures are not lagging behind the other temperatures. In particular, species 1, which is leading even the translational temperatures. Such a result is surprising but can be observed in the results by the more detailed General Boltzmann Equation (GBE). 17 Even though further work is required to justify this phenomenon, it is not considered a drawback introduced by the new kinetic model, since it observed by the DSMC and GBE results.
The assumption for a constant Zv throughout the shock and its value is too big pre-shock, where the free-stream temperature suggests no vibrational excitation. It is also possible that the values of post-shock are too low, especially for the case of Zv = 25. It is clear that the assumption of a constant Zv is a simplification, which is reasonable for initial method/solver development that should be improved upon and addressed in future work.
Until now, we were focused on the different types of temperatures. If we turn our attention to the species differences we observe that species 1 (nitrogen) and 2 (oxygen) for both translational and rotational temperatures differ very little for both the kinetic model and DSMC. This is consistent through the different Mach numbers and concentrations. Certainly, increasing the mass ratio between species will separate the same type of temperature between the two species (specifically the translational and rotational temperatures), as observed for the translational temperature in the monoatomic mixture. 5,49 The difference between species vibrational temperature is bigger, in comparison to the difference between species for the other two types of temperatures, even for the small mass ratio considered here. The number of vibrationally excited DoFs between the oxygen and nitrogen gases is also different, which is a key emphasis of the model. This is due to the different species characteristic vibrational temperatures. It is highly likely that for a species and temperature dependent Zv s , this difference will be more pronounced. Similarly, having a species Zr s for the rotational collision coefficient will separate the two rotational temperatures. Here, we observe the biggest difference between the kinetic model and DSMC: in the vibrational temperatures and even here, the two numerical models still show good agreement.
This discussion of the results shows that the kinetic model and DSMC results follow the same physical trends, as well as some of the limitations of the current assumptions. It demonstrates that the model itself provides the accurate results for the macroscopic variables under strong non-equilibrium conditions, i.e., as created by M∞ = 5 normal shock conditions. The validation and evaluation of the model beyond M∞ = 5 need further study.

C. Mixture properties of the normal shock profile
From here on, the physics of the flow will be investigated further. The vibrational degrees of freedom through the shock are considered. At the considered temperatures, the translational and rotational degrees of freedom are considered fully excited. In contrast, the vibrational degrees of freedom are species-and temperaturedependent. In Fig. 8 At the specified free-stream temperature, the pre-shock values are expected to be negligible. With the rise in temperature through the shock wave, the number of vibrational DoFs increases. The noninteger values of δ(T) in the pre-shock flow-field are close to 0, increase through the shock, and reach the maximum values of postshock. This maximum value depends on the temperature jump but asymptotically approaches the value of 2 for very high vibrational temperatures. Figure 8 shows that the values of δ are non-integer values, as expected, and are in the bounds discussed for δ(T): [0; 2). Furthermore, the values of δ for oxygen are much higher than the values of δ for nitrogen for both vibrational collision numbers Zv. Therefore, we observe that oxygen is more vibrationally excited than nitrogen at the same temperature, which can be explained from the lower characteristic temperature for O 2 . This is consistent with nonequilibrium thermodynamics theory. 10,11,15,53 In Figs. 8(a) and 8(b), the effect of the larger vibrational collision number Zv = 100 is significant. For the two collision numbers, the jump in the value of δ is identical, but for Zv = 100, the slope of the jump is less steep and the number of mean free paths that are required to post-shock equilibrium is much larger than for Zv = 25. Comparing the plots in (a) and (b), it is clear that for larger Mach numbers, the level of excitation of both species is much higher. This is due to the larger temperature ARTICLE scitation.org/journal/adv jump caused by the M∞ = 5. In contrast, the concentration change in (c) has not really affected the number of vibrational degrees of freedom for either species. These figures demonstrate that there is a significant difference between species in the level of excitation and modeling them separately is important.
Inspecting the velocity drift between species through the shock demonstrates the importance of having separate species velocity in the formulation of the kinetic model. In Fig. 9, the difference between velocities u 2 − u 1 is shown for cases when the concentration ratio is 80% N 2 and 20% O 2 . The cases considered are M∞ = 3 and M∞ = 5 with vibrational collision numbers Zv = 25 and Zv = 100. The results for M∞ = 3 are in red and for M∞ = 5 are in green, with Zv = 25 displayed with a solid line and Zv = 100 with a dashed line and a symbol. We observe from this figure that the velocity difference does not change when the vibrational coefficient is modified. However, the higher the Mach number, the bigger the velocity difference. For the diatomic mixture, the maximum values at β = 1.14 are: at M∞ = 3, u 2 − u 1 ∼ 0.08u ref and for M∞ = 5 is u 2 − u 1 ∼ 0.17u ref , where u ref is the most probable speed for the mixture under free-stream (reference) conditions. This drift velocity will increase with increasing mass ratios between species, as it was shown in the parametric study in the monoatomic mixture gas, 49 as well as with the increase in the free-stream velocity. velocity distribution functions for the case at the highest Mach number M∞ = 5 and the larger vibrational collision number Zv = 100 at 80%:20% nitrogen to oxygen concentration ratio in Figs. 10 and 11. In Fig. 10, the focus is on comparing the distributions far upstream and far downstream, while Fig. 11 shows the behavior of the distribution functions of the two species through the shock.

Velocity distribution functions
The pre-shock and post-shock values of the translational, rotational, and vibrational distribution functions are shown for the oxygen in Fig. 10. Both distributions are Maxwellian and are set by the macroscopic variables defining the shock jump. The pre-shock distribution functions have a sharp peak, which is however low in value in comparison with the post-shock Maxwellian functions. This is due to the higher values of the post-shock macroscopic variables. Here, it is important to mention that a higher value of the post-shock number density increases the value of distribution function, while high temperature of post-shock actually has the opposite effect. We will observe that the first effect seems to dominate here. For this reason, the pre-shock values are plotted on a third of the postshock scale. At the same time, the post-shock distribution functions are more spread out over a larger velocity mesh. This is consistent with theoretical observations: a larger temperature leads to a larger spreading of the distribution function. The classical question in applications of kinetic models in DVM solvers for the compromise between velocity bounds (to capture the tails of the post-shock values) and the fine velocity mesh (to capture the sharp spike in the preshock values) for the velocity mesh is demonstrated. Through visual observation of these figures, the mesh demonstrates to be adequately selected.
Furthermore, the pre-shock values are most prominent for the translational distribution. It has a higher maximum than the rotational and vibrational functions. Note that the relationship between the reduced translational distribution function and the rotational distribution function, following the expressions in Eq. (40) and in non-dimensionalized form G rot = (Tr/2β)G tran , where Tr is the nondimensionalized rotational temperature and β = m 2 /m 1 is the mass ratio between species. Since the non-dimensional rotational temperature is taken in the pre-shock as Tr = Tr(dimensional)/T∞ = 1, it follows that Tr < 2 * β in the pre-shock then G rot < G tran , while in post-shock, the value of Tr jumps and G rot > G tran . Meanwhile, there is no vibrational energy in the pre-shock part of the flow domain. Observe that the color scale for the vibration is three orders of magnitude lower for the purposes of presentation of this figure. The values are approaching zero, consistent with the recognition that no vibrational excitation is present at the pre-shock position. Meanwhile, the biggest pre-and post-shock difference between distribution functions is for the vibrational function G vib (O 2 ), with the highest jump in values of three orders of magnitude. The vibrational distribution function is significantly affected by the choice of freestream temperature. Figure 10(e) shows a more pronounced peak and shape in the pre-shock for T∞ = 1000 K, for example. This selection would also affect the post-shock value dramatically. In Fig. 11, the distribution functions for both species: nitrogen (first column) and oxygen (second column) half-way through the shock are shown. The first observation is that the distribution is non-Maxwellian. It has some features of both the pre-shock and post-shock distributions and the transition between them. The values for nitrogen are expectedly higher than the values for the oxygen due to the four-times higher concentration of nitrogen (i.e., number density multiplied of the distribution), while the shape of the distributions is identical between species. Note that due to the lower value of excited vibrational degrees of freedom δ(N 2 ) as shown in Fig. 8(c), the distribution function of the nitrogen is expected to have a lower value of the peak. However, this is compensated by the higher levels of concentration. A clear bimodal behavior can be observed in Figs. 11(e) and 11(f) for the two species, suggesting that for M∞ = 5, two most probable velocities can be identified. These two velocities are usually closely related to the pre-and post-shock mean values. We observe that the maximum values of the distribution functions for the nitrogen gas are higher. This is due to the higher concentration of species 1 in the mixture for the case with the 80%:20% ratio of nitrogen to oxygen concentration. Furthermore, the displayed results for the distribution functions at Zv = 100 are very close to the values for the lower vibrational collision number Zv = 25. The maximum of the vibrational distribution function G vib s is achieved further downstream relative for the other distribution ARTICLE scitation.org/journal/adv functions at Zv = 100, which is consistent with the findings for the macroscopic variables. This is apparent by inspecting the mid-shock value of the distribution functions and compare it with the equivalent post-shock distribution function. The maximum value of the translational distribution function is of the same order of magnitude mid-shock and post-shock. The rotational distribution function is 1 order of magnitude smaller mid-shock than post-shock, while the vibrational distribution is at 2 orders of magnitude smaller mid-shock in comparison to the maximum value of post-shock. Therefore, the effect of the different relaxation rate of translation, rotation, and vibration can be observed in the velocity distribution functions.

VII. CONCLUSION
Three single-species diatomic models with different levels of complexity were presented. Furthermore, a new diatomic model for the binary gas mixture was derived. The model has separate speciestranslational, -rotational, and -vibrational temperatures. Species diffusion is accounted for with a relaxation parameter, allowing the two species velocities to separate. The continuum limit and transport properties are discussed in detail and compared with the properties of the single-species diatomic models and the mixture monoatomic model. The models are found to have good mathematical properties and to be consistent with the previous models. A disadvantage of the presented model is the lack of Prandtl number correction in the continuum limit. However, diatomic mixture models are still a very challenging and active area of research, and the introduced model is a strong step toward a well-understood, derived, and described model for gas mixtures with drift velocity and speciesspecific translational, rotational, and vibrational temperatures. The new kinetic model for a binary mixture of diatomic gases is evaluated numerically for a profile of a normal shock wave and shows good agreement with the DSMC benchmark results. The good agreement of the results is a stepping stone in the development of modeling air flow in rarefied conditions for practical applications using deterministic and efficient kinetic models. Further testing will expand the application possibilities and capabilities of the model.