Multi-species modeling in the particle-based ellipsoidal statistical Bhatnagar-Gross-Krook method including internal degrees of freedom

The implementation of the ellipsoidal statistical Bhatnagar-Gross-Krook (ESBGK) method in the open-source particle code PICLas is extended for multi-species modeling of polyatomic molecules, including internal energies with multiple vibrational degrees of freedom. For this, the models of Mathiaud, Mieussens, Pfeiffer, and Brull are combined. In order to determine the transport coefficients of the gas mixture, Wilke's mixing rules and collision integrals are compared. The implementation is verified with simulation test cases of a supersonic Couette flow as well as a hypersonic flow around a 70{\deg} blunted cone. The solutions of the ESBGK method are compared to the Direct Simulation Monte Carlo (DSMC) method to assess the accuracy, where overall good agreement is achieved. In general, collision integrals should be preferred for the determination of the transport coefficients, since the results using Wilke's mixing rules show larger deviations in particular for large mass ratios. Considering the computational efficiency of the methods, a considerable reduction in computational time is possible with the ESBGK model.


Introduction
Numerical simulations of fluid dynamics in space applications, micro and nano flows, as well as vacuum technology present significant challenges for applied numerical methods.These applications involve large density gradients, ranging from the continuum to the free molecular flow regime, and require versatile numerical approaches due to their multiscale and non-equilibrium nature.Although a highly accurate solution of these flows can be achieved by using the well-established Direct Simulation Monte Carlo (DSMC) [5] method, the latter requires excessive computational effort in the transition and continuum regimes.Thus, the coupling of DSMC with a computationally efficient method is desirable.
For continuum flows, computational fluid dynamics (CFD) are typically used.However, many problems arise in the coupling with the DSMC method due to the very different underlying approaches of both methods.Especially the statistical noise of the DSMC method is problematic at the boundaries between DSMC and CFD [6,7,8,9].
Particle-based continuum methods have emerged as promising alternatives, bridging the transition regime.Recent developments include but are not limited to the Fokker-Planck approach [10,11,12,13] and the Bhatnagar-Gross-Krook (BGK) [14,15,6,16] model, which is the focus of this paper.The latter approximates the collision integral of the Boltzmann equation by a relaxation process.With this, the mean free path and the collision frequency do not need to be resolved, which means that the choice of time step and particle weighting factor for particle simulations can be less restrictive.The necessary local resolution is rather defined by the gradients of the macroscopic moments of the distribution function present in the flow, such as density, velocity or temperature, similar to that described in Pfeiffer and Gorji [11].The time step is then determined by the stiffness of the BGK collision term, i.e. the relaxation frequency to be resolved.The advantage over DSMC, however, is that the BGK equation no longer represents a jump process but can use different time integration and space interpolation methods to achieve significantly coarser resolutions.While this has been a focus of research in connection with BGK in the context of discrete velocity methods [17,18] for a long time, progress has also been made in recent years on the particle methods side [19,20,21,22].Thus, less computational costs are expected compared to DSMC simulations, in particular for low-Kn regimes.
The particle-based ellipsoidal statistical BGK (ESBGK) [23] and the Shakhov BGK [24] methods were investigated in particular [25,26], both producing the correct Prandtl number of the gas.It was showed that in some cases the ESBGK model is more robust and more efficient in the particle context [25,8].In addition, the ESBGK method fulfills the H-theorem [27], which is why this approach is followed here.
An overview of different BGK models for gas mixtures is available from Pirner [28].A rough distinction is made between two types of mixture models.The first type are models with a sum of several relaxation terms, so that each species combination has its own relaxation term [29,30,31,32,33,34].The advantage here is that the correct collision rates and thus energy/moment exchange rates between the individual species can theoretically be reproduced correctly according to the Boltzmann equation.The disadvantage, however, is that the models are significantly more complex and therefore require more computing time.In addition, the choice of free parameters such as the different relaxation frequencies is challenging from a physical point of view, for example if certain viscosities or Prandtl numbers are to be achieved for the entire mixture.This is especially true for multi-species mixtures with more than two components and when additional exchange terms for internal degrees of freedom are added in the polyatomic case.On the other hand, there are models with only one relaxation term on the right hand side [35,3,4,36,37].The obvious disadvantage is that the correct energy and momentum exchange times cannot be modeled with this approach.However, there are some advantages for the practical application of the models.The models are less complex and therefore require less computing time.In addition, the collective behavior of the gas mixture, such as viscosity, heat flow, Prandtl number, etc., can be modeled correctly with significantly fewer parameters, whereby established mixing rules can be used.Since the left-hand side in the form of advection is still modeled correctly, the error due to the use of only one relaxation term can be small, depending on the application.
The focus of this paper is the extension of the particle-based ESBGK model to gas mixtures including non-equilibrium states of internal degrees of freedom of di-and polyatomic molecules, including quantized vibrational states by using a model with only one relaxation term.The work is based on the models of Mathiaud et al. [1], Pfeiffer [2] and Brull [3,4].In the applications presented here, this model using only one relaxation term is already capable of accurately representing non-equilibrium effects in complex 2D flows.This was already demonstrated in Pfeiffer et al. [38] using pure atomic mixtures and has been successfully extended here to mixtures with any number of atomic and polyatomic species, including non-equilibrium effects of internal excitation energies.In the models with one relaxation term mentioned above, either binary molecular mixtures [37], atomic mixtures [35,36,3], or equilibrium in the internal excitation states [4] have been assumed so far.In the specific case of stochastic particle methods, models with only one relaxation term have the additional advantage that solely the moments of the entire mixture, not the individual species, are required.Thus, a significantly smaller particle representation of the distribution function, which means fewer simulation particles, can be used compared to the case where all individual moments of the distribution function are needed, as it is in the case of multi-relaxation term models.

Theory
The BGK operator is an approximation of the Boltzmann collision integral as relaxation process of the particle distribution function f s (x, v, t) of a species s at position x and with velocity v towards a target distribution f t s [14] using the relaxation frequency ν.In a mixed model of atoms and molecules, the target distribution function differs between the two because the internal degrees of freedom (DOF) for the molecules must be taken into account.While the atoms only carry translational DOFs, the molecules carry vibrational and rotational DOFs additionally.Therefore, the target distribution function of atoms only consists of a translational part depending on the particle velocity v, whereas the distribution function for molecules can be separated into the translational, the rotational and the vibrational part as demonstrated in Mathiaud et al. [1], Dauvois et al. [39] depending also on the rotational energy E rot and the vibrational quantum number i vib .

Macroscopic Flow Values
The macroscopic flow values particle density n, flow velocity u and temperature T tr of each species s are defined as the moments of the particle distribution function: with the thermal particle velocity c = v − u from the particle velocity v, the total energy E and the thermal energy E. Furthermore, the macroscopic mean values of the flow are given by: For molecular species, it is also necessary to define the mean rotational energy ⟨E⟩ rot and the mean vibrational energy ⟨E⟩ vib , which are directly linked to the rotational temperature T rot and the vibrational temperature T vib : ξ are the degrees of freedom of rotation and vibration.In the diatomic case, ξ rot = 2 applies, while in the polyatomic case, ξ rot = 2 applies for linear molecules and ξ rot = 3 for non-linear molecules [40,41].For the vibration, the simple harmonic oscillator (SHO) model is assumed with the characteristic vibrational temperature Θ vib .Due to the quantization, the degrees of freedom are not fixed but temperature-dependent.In the polyatomic case, γ is the number of vibration modes as a function of the total number of atoms a in the molecule [40]: In the case of a diatomic molecule, γ = 1 follows.With the assumption of the SHO, an expression for the mean energy and vibration degrees of freedom as a function of temperature T vib,s can be found:

ESBGK mixture model
The model proposed here combines the model of Mathiaud et al. [1], Pfeiffer [2] and the model of Brull [3,4] into a mixture model that allows for non-equilibrium effects in the internal degrees of freedom.Different than the standard BGK model using a Maxwellian target distribution function, the ellipsoidal statistical BGK (ESBGK) model produces the correct Prandtl number by using an anisotropic Gaussian distribution [23] f ES,tr with the mass density of the considered species ρ s and the average flow velocity u.The anisotropic matrix A s reads: T tr,rel T tr P − T tr,rel I and consists of the identity matrix I and the pressure tensor which are both symmetric [3].The parameter α is a variable of the model depending on mass fraction, density fraction and internal degrees of freedom [4]: The internal degrees of freedom are calculated as the sum of the vibrational and rotational degrees of freedom, respectively: The model used and presented in this paper uses one relaxation term per species s.It fulfills the indifferentiability principle so that the model reduces to a single species model for identical species.Additionally, the Maxwellian distribution is produced in the equilibrium state.The relaxation frequency ν of the model is defined by with the viscosity of the mixture µ and the translational temperature T tr .Pr is the targeted Prandtl number of the gas mixture, which is calculated using Wilke's mixing rules [42] or collision integrals [43] (see Section 2.3).
For molecular species, the rotational and vibrational components of the distribution function must also be defined, which are taken directly from Mathiaud et al. [1]: The temperatures T tr,rel , T rot,rel,s and T vib,rel,s introduced here are chosen so that the relaxation corresponds to the Landau-Teller relaxation: Here Z r,s is the collision number of the internal degrees of freedom r per species s and τ s,C = 1/ν s,C is the collision time of the gas per species.The latter is different from τ = 1/ν, the relaxation time of the BGK equation (1), whereas a detailed discussion can be found in Mathiaud et al. [1].The collision frequency for the variable hard sphere model (VHS) of each species s in a mixture of M species is calculated by [44]: with If direct VHS data are available for each collision combination, these can also be used accordingly here [45,46,47].With these definitions, T tr,rel , T rot,rel,s and T vib,rel,s can now be defined with r = rot, vib: ⟨E⟩ r,s (T r,rel,s ) = ⟨E⟩ r,s (T r,s ) ⟨E⟩ tr (T tr,rel ) = ⟨E⟩ tr (T tr )

Gas mixture properties
The ESBGK model requires the calculation of the correct Prandtl number of the gas mixture to assess the relaxation frequency from Equation (21).Here, µ is the viscosity, K is the thermal conductivity and c p is the specific heat of the considered mixture.The latter is calculated as: For the calculation of the mixture viscosity and thermal conductivity, any source can be used in the presented model.Basically, the same problems also arise with multi-species CFD codes and therefore the same solutions can also be applied.Typically, the transport properties are calculated on the basis of Wilke's mixture rules [42] or the first approximation of the transport properties using collision integrals [43], whereby various libraries are also available that can be used directly [48,49,50].However, in order to achieve the same transport coefficients as far as possible when coupling with DSMC, the procedure for the Variable Hard Sphere (VHS) potential model is presented in this paper.

Wilke's mixing rules
Wilke's mixing rules [42] are used for the calculation of the transport coefficients of the mixture.The well-known exponential ansatz is used for the viscosity of each species µ s , which can be calculated with: Here, µ ref,s is the reference dynamic viscosity at a reference temperature T ref,s for each species s, and ω VHS,s is a speciesspecific parameter of the Variable Hard Sphere (VHS) model [5,46,45,47].Using the VHS reference diameter d ref,s , the reference dynamic viscosity for a VHS gas is defined as [51]: The mixture viscosity of M species is then calculated by [42]: The thermal conductivity of each species K s can be calculated using the Eucken's relation [52] and the correction by Hirschfelder [53] with the viscosity: Here, c v,s = ξ s k B /2m s is the specific heat capacity at constant volume for species s, with the corresponding degrees of freedom for translation ξ tr = 3 for each species and the internal degrees of freedom ξ int,s from Equation (20).In general, c v = c p − k B /m applies.The prefactors used in Equation ( 36) are defined as f tr = 5/2 and f int = 1.328.With that, the thermal conductivity of the mixture K is calculated:

First approximation of transport properties
Another calculation possibility is the first approximation to the viscosity of species s depending on the collision integral Ω (2)  s (2), which is given by [43] .
The mixture viscosity is determined by where b s is the contribution of each species to the total mixture viscosity and is determined by solving the system with the mole fraction χ, the density ρ ′ k of species k when pure at pressure and temperature of the actual gas mixture, the parameter A sk , defined by sk (2) 5Ω (1)  sk ( 1) and the binary diffusion coefficient with the reduced mass m * sk : .
The mixture thermal conductivity K is calculated by with a s being the translational species contribution to the total mixture thermal conductivity and K rot,s as well as K vib,s being the rotational and vibrational contributions of the species s [44].The factors a s are determined by solving the system Here, K s is the first approximation of the thermal conductivity of species s: and the parameter B sk is defined by sk (3) 5Ω (1)  sk ( 1) The collision integrals for the VHS model are given by [54] Ω VHS,( 1) with the VHS parameters d ref , T ref and ω.Generally, for s = k, the translational temperature of the species T tr,s is used, while for s k, the translational temperature in the cell T tr is used as an approximation in order to reduce noise.The rotational and vibrational contributions K r,s with r = rot, vib to the thermal conductivity of a species s are calculated as [44]:

Implementation
The proposed ESBGK method gas mixtures with internal degrees of freedom is implemented in the open-source particle code PICLas [55].A stochastic particle approach is chosen for the implementation [25,2,38,56], although the model can of course also be used in discrete velocity (DVM) codes.Especially to save simulation time in DVM multi-species simulations, the model should be extended to a reduced model in the future [57, 39, 1].

Relaxation and sampling
Regarding movement and boundary conditions, particles are treated in the same manner as in the DSMC method.But instead of performing binary collisions between the particles, all particles in a cell relax with a probability P in each time step ∆t: which corresponds to Equation (1) in a stochastic interpretation as described in Pfeiffer [25].The first order method from Pfeiffer [25] is used here.In the future, however, this model will be extended to the second-order particle method from Pfeiffer et al. [19].The relaxation frequency ν, which is the same as in Equation ( 1), is evaluated in each time step for each cell according to Equation ( 21) using the transport coefficients of the mixture calculated in the code.
The particles that are chosen to relax are then sampled from the target distribution (see Equation ( 15)).Here, an approximation of the transformation matrix S [58] with A = SS is used to transform a random vector chosen from a Maxwellian to a velocity vector chosen from the ESBGK distribution: For a detailed description of the sampling process, the reader is referred to Ref. [25], in which also alternative methods for sampling are described.At this point, a special solution for the stochastic particle approach must be chosen.Even though the presented model is both momentum and energy conserving, additional effort must be made for both in a stochastic approach.The basic problem is that the new states of the particles (velocities, internal energies) are chosen randomly from the distribution functions.Thus, in this random step, conservation of energy and momentum is not automatically given and one would need a large number of particles for stable simulations.To construct energy conservation for the stochastic particle method, the energy difference between the old and new sampled energy must be distributed over the particles in such a way that energy conservation is given.The idea is to distribute this energy difference evenly according to the respective degrees of freedom to the translational energy of all particles in the cell and the internal degrees of freedom of the particles relaxing in the degrees of freedom within the cell, as this has led to the best results in previous work [2].An artificial equilibrium temperature T equi is introduced for this purpose.Starting from this point, molecules are chosen to relax in the rotation and vibration towards the equilibrium temperature T equi with the probability so that Equations ( 1), ( 22) and ( 23) are fulfilled.To save calculation time, instead of using the probability P to relax the particles to T vib,rel,s and T rot,rel,s corresponding to Equation ( 22) and ( 23), only the fraction τ Z r,s τ s,C relaxes directly to T tr according to Equation (29), which leads to the same result but with significantly fewer relaxing particles.Additionally, a correction parameter β r,s for the rotation and vibration, respectively, is defined per species in order to follow the idea of equal distribution of energies for the conservation of energy as described: Since T equi directly depends on β rot,s and β vib,s of each species, the solution of this equation system is obtained numerically.For this purpose, the following system must be solved: with the starting values ξ vib,s (T 0 equi ) = ξ vib,s (T vib,s ), (56) and iteration step n, until the accuracy ϵ is reached with the condition T n+1 equi − T n equi < ϵ.N is hereby defined as the total number of particles in the gas mixture, while N s is the particle number per species.
The new rotational energy of a molecule p, which is chosen for a rotational relaxation according to the probability in Equation (51), is calculated by: using a random number R p .In the same manner, the new vibrational energy of a vibrational relaxing particle is calculated: The post-relaxation energies are hereby denoted with * .To fulfill the energy conservation, these energies need to be scaled additionally.Note: At this point the vibration is still continuous but already with the correct quantized degrees of freedom.The transition to quantized states happens with the conservation of energy.

Energy and momentum conservation
The energy and momentum conservation is based on Ref.
[2] and extended to polyatomic molecules and mixtures of different gas species.In the implemented scheme, only relaxations in translation-vibration (T − V) and translation-rotation (T − R) are allowed directly.In the following, the process is described in greater detail.Here, the parameters after energy conservation are denoted with ′ .
For the energy conservation process, the vibrational energy is scaled first.The energy E T −V is the sum of the total translational energy of the gas mixture and the vibrational energy of all relaxing particles in vibration, N vib , before the relaxation process: It should be distributed equally over all translational and vibrational DOFs after the relaxation process to fulfill the energy conservation.This is reached with a scaling factor α vib,s , that is calculated per species as: With this, the new vibrational energy of each particle p after the relaxation and energy conservation processes is for the continuous handling of the vibrational energies and including the zero-point energy corresponding to the second sum.
For quantized vibrational energy states, additional steps are to be performed.The energy per vibrational mode α vib,s E * vib,p, j is reformulated to a quantum number i p, j with a random number R p, j ∈ [0, 1): If the condition is fulfilled, is the new energy of this vibrational mode, including the zeropoint energy.Otherwise, a new quantum number of this mode is calculated with a new random number R p, j as: until Equation ( 63) is fulfilled.E ′ vib,p, j is then calculated with this new quantum number using Equation (64) subsequently.
The energy E T −V is subsequently updated with before the algorithm is repeated for each vibrational mode j of each vibrational relaxing particle p.In the end, the new energy of the each particle p is: Afterwards, the remaining energy for translation and rotation is the sum of the remaining E T −V and the rotational energy of all rotational relaxing particles before the relaxation process: Here, N rot is the number of the rotational relaxing particles.E T −R should be distributed equally over all translational and rotational DOFs similar to E T −V .For this, a scaling factor α rot,s is calculated per species as: Using this, the new rotational energy of each particle p after the relaxation and energy conservation processes is calculated with: Different than for the internal DOFs, and as described earlier, all particles in the gas mixture are scaled in the translation, whether they relax or not.The scaling factor α tr is derived to be Finally, the new particle velocities are calculated using with the average flow velocities If no relaxation occurs for a particle p, v * p = v p applies.Using this approach, energy conservation is guaranteed and momentum conservation is ensured due to as described in Pfeiffer [2].

Simulation Results
The previously described multi-species ESBGK implementation including molecules with internal degrees of freedom is verified with simulation test cases of a supersonic Couette flow as well as a hypersonic flow around a 70°blunted cone.The different approaches for the calculation of the transport coefficients, Wilke's mixing rules (denoted by Wilke) and the direct use of collision integrals (denoted by CollInt), are compared using the VHS model.The VHS species-specific parameters are listed in Table 1.The parameters for a collision pair of unlike species are determined as an average (referred to as collisionaveraged).

Supersonic Couette flow
The first test case is a three-dimensional simulation of a supersonic Couette flow.A mesh with 100 cells in y direction and a single cell in x and z directions is used.For the second test case of N 2 -N with Kn VHS ≈ 0.121 and Kn VHS ≈ 1.21 as well as for test cases four and five, the mesh is refined and thus contains 400 cells in y direction.The time steps and particle weighting factors are chosen so that the mean free path and the collision frequency are resolved.The boundaries in y direction have a velocity of v wall,1 = 350 m s −1 and v wall,2 = −350 m s −1 , respectively, assuming diffuse reflection and complete thermal accommodation at a constant wall temperature of T wall = 273 K.In x and z directions, all boundaries are periodic, leading to particles reappearing on a side after leaving the opposite side.Different gas mixtures are simulated, all initialized at v 0 = 0 m s −1 and T 0 = 273 K.An overview of all Couette flow test cases is given in Table 2.
First, a N 2 -He mixture is simulated, which represents a challenging case due to the relatively large mass differences.Mixing ratios of 50 %-50 %, 20 %-80 % and 80 %-20 % with a initial particle density of n 0 = 1.3 • 10 20 m −3 are tested to estimate the influence of composition on the simulation results.The correct Prandtl number approximation is crucial for an accurate temperature representation in the Couette case.However, since it is not a linear function of the mixing ratio [59], these three ratios were chosen to capture the resulting Prandtl numbers as a function of the mixing ratio.In Figure 1   temperature is displayed for the different mixing ratios.In general, the ESBGK results using collision integrals are in excellent agreement with the DSMC results with a maximum error of 3.6 %, while larger deviations up to 8.5 % are visible using Wilke's mixing rules.This is probably due to the large mass ratio which leads to a poorer approximation of the Prandtl number.The results slightly differ with the mixing ratio of N 2 and He.In Figure 2, the rotational temperature is compared for the different mixing ratios.The results are similar to the results of the translational temperature due to the thermal equilibrium in the Couette flow.Therefore, only the translational temperature plots are shown for the following test cases.
The second test mixture is a 50 %-50 % N 2 -N mixture, initialized with the same particle density as for the first case, n 0 = 1.3 • 10 20 m −3 , as well as lower densities by one and two orders of magnitude, respectively.This mixture was chosen because it represents a typical mass difference of species in air mixtures (N 2 ,N,O 2 ,O, NO).The initial densities correspond to Knudsen numbers of Kn VHS ≈ 0.0121 − 1.21, whereby the implemented ESBGK method can be verified across different flow regimes.The results of these simulations for the translational temperatures are shown in Figure 3. Again, the results of the ESBGK simulations using collision integrals are in very good agreement with the DSMC simulations.Due to the lower mass ratios, also Wilke's mixing rules achieve good agreement with the DSMC results.For higher Knudsen numbers, the errors for both ESBGK methods are increasing which is also consistent with observations for the one species ESBGK model [25].
Here, deviations of up to 32 % occur.The performance of the ESBGK method thus decreases for increasing Knudsen numbers.
As a third test case, a mixture with two molecular species (polyatomic and diatomic) is chosen to examine the behavior of the model with two molecular species: 50 %-50 % CO 2 -N 2 .Again, n 0 = 1.3 • 10 20 m −3 is the initial particle density.The results for the translational temperature are depicted in Figure 4.All simulation methods show very good agreement with errors below 3.8 %.
With test cases four and five, two more realistic atmospheric gas mixtures with a larger number of constituents are simulated to test the behavior of the model with multi-species mixtures.From now on, only the DSMC and ESBGK solutions using collision integrals for the calculation of the transport coefficients are compared since the results using Wilke's mixing rules show larger errors, in particular when considering mixtures with large mass ratios.
Case four shows the results for a 78 %-21 %-1 % N 2 -O 2 -Ar mixture which corresponds to a typical air composition close to the ground, again initialized with n 0 = 1.0 • 10 20 m −3 .In   mimic a gas composition after dissociation processes in the air.The ESBGK results show very good agreement with the DSMC results, both illustrated in Figure 6, with maximum deviations of 2.6%.Both case 4 and case 5 emphasize that good results can be achieved with the model even with multi-species mixtures.

70 degree blunted cone
To test the behavior of the ESBGK mixture model in particular the presence of strong gradients such as across shock waves, simulations of a hypersonic flow around a 70°blunted cone are done.The geometry of this cone is illustrated in Figure 7.The simulations are performed two-dimensional axisymmetrical.At the surface of the cone diffuse reflection and complete thermal accommodation at a constant wall temperature of T wall = 300 K is assumed.The inflow conditions are listed in Table 3.A 50 %-25 %-25 % N 2 -O 2 -NO gas mixture is chosen for Cases 1-3, while a 50 %-50 % CO 2 -N 2 mixture is used for   Case 4. The first three cases were chosen to test the behavior of the model for different Knudsen numbers on the more complex case and the fourth case was selected to test the behavior with a polyatomic molecule.The mesh contains 5064 cells for Cases 1 and 4 with 1.06 • 10 6 total simulation particles and 16320 cells for Cases 2 and 3 with larger particle densities in the free-stream.For Case 2, 3.13 • 10 6 particles are in the simulation, while for Case 3, there are 7.38 • 10 6 ESBGK simulation particles and 1.45 • 10 7 DSMC simulation particles, respectively.For the DSMC method, a time step of 1 • 10 −8 s is used for Cases 1, 2 and 4, and a smaller time step is chosen for Case 3 with 2 • 10 −9 s, while for the ESBGK method, the time step is 2.5 • 10 −8 s.The results using collision integrals directly and Wilke's mixing rules are almost identical for all simulations due to the small mass differences between the gas constituents and thus good approximation with Wilke's mixing rules, which is why the latter are not shown separately.

Case 1
For the first test case, a density in the free-stream of n ∞ = 3.7 • 10 20 m −3 is assumed, corresponding to Kn VHS,∞ = 0.147, to show the abilities of the ESBGK model for nonequilibrium flows.
The simulation results of the mean flow variables of the mixture are shown in Figure 8. Overall good agreement is visible between the results of both methods.
In Figure 9, the mean translational, rotational and vibrational temperatures of the gas mixture are compared for the DSMC and ESBGK simulations.The ESBGK model predicts a earlier onset of the translational temperature increase compared to DSMC, which results in a wider shock region, also visible in Figure 10, where the translational temperature in the flow field is shown in comparison.However, this behavior is typical for the ESBGK method, even for just one species [25,26,38].The agreement in the post-shock region is excellent.Figure 9 also shows the vibrational temperatures of each species at y = 20 mm.Even in the non-equilibrium region in the wake of the 70°cone, the agreement between the results of the ESBGK and DSMC methods is very good.Furthermore, a relatively large deviation in the vibrational temperature is visible in the front of the shock in Figure 9.The reason for this is stochastic noise.Due to the very low inflow temperature in order to achieve the high Mach number and due to the quantized treatment of the vibration, only very few particles are vibrational excited in this area, which leads to this high statistical noise.However, using significantly more particles only makes limited sense here, as this region in terms of vibrational excitation has no influence on the behavior of the flow.Comparing the vibrational temperatures is only meaningful closer to the shock, above temperatures of around 400 K, as a sufficient number of particles is excited here.
The heat flux on the cone surface is compared in Figure 11.For both the flow-facing cone surface as well as the rear cone part, excellent agreement between both methods is shown.

Case 2
As a second test case, a particle density in the free-stream of n ∞ = 1.48 • 10 21 m −3 is assumed, corresponding to Kn VHS,∞ = 0.0366, and representing a transition regime case.
The simulation results of the mean values of the velocity in x direction and the particle density of the mixture are shown in Figure 12 with overall good agreement between the results of both methods.
In Figure 13, the mean translational, rotational and vibrational temperatures of the gas mixture are compared for the DSMC and ESBGK simulations.Similar to Case 1, the ES-BGK model predicts an earlier onset of the translational temperature increase compared to DSMC, resulting in a wider shock region.Excellent agreement is visible in the post-shock region.Figure 13 shows the vibrational temperatures of each species at y = 20 mm additionally.The noise behavior is also visible here, similar to Case 1. Again, there is very good agreement between the results of the ESBGK and DSMC methods even in the non-equilibrium region in the wake of the 70°cone.
A comparison of the heat flux on the cone surface is shown in Figure 14.For both the flow-facing cone surface as well as the rear cone part, both the ESBGK and the DSMC simulation results are in excellent agreement.

Case 3
The third test case represents a continuum case with a particle density in the free-stream of n ∞ = 5.92 • 10 21 m −3 , corresponding to Kn VHS,∞ = 0.0092.
The the mean values of the velocity in x direction and the particle density of the mixture are compared in Figure 15 for the ESBGK and the DSMC models, both showing excellent agreement.
In Figure 16, the mean translational, rotational and vibrational temperatures of the gas mixture are shown.The earlier onset of the translational temperature increase of the ESBGK model is again visible compared to DSMC, but the difference in the translational temperature profile is smaller compared to Cases 1 and 2. This is due to the higher density, leading to a narrower shock region in general.Again, excellent agreement is visible in the post-shock region.Furthermore, Figure 16 illustrates the vibrational temperatures of each species    The simulation results for the heat flux on the cone surface is compared in Figure 17 for both of the models.There is good agreement of the results for the flow-facing cone surface as well as only little deviations for the rear cone part.

Case 4
For Case 4, a 50 %-50 % CO 2 -N 2 is chosen to evaluate the accuracy of the ESBGK model for a polyatomic-diatomic mixture with a higher number of degrees of internal freedom in particular in vibration.A particle density in the free-stream of n ∞ = 3.7 • 10 20 m −3 is chosen similarly to Case 1, which corresponds to Kn VHS,∞ = 0.108 for this mixture.
In Figure 18, the simulation results of the mean values of the velocity in x direction and the particle density of the mixture are compared for the DSMC and ESBGK models.Same as for the 50 %-25 %-25 % N 2 -O 2 -NO gas mixture from Case 1, there is overall good agreement.
The mean translational, rotational and vibrational temperatures of the gas mixture are illustrated in Figure 19.Again, a wider shock region is visible for the ESBGK results due to an earlier onset of the translational temperature increase compared to DSMC, while there is very good agreement in the post-shock region.Also, the vibrational temperatures of each species at y = 20 mm are shown.While there is excellent agreement between the results of the ESBGK and DSMC methods in the non-equilibrium region in the wake of the 70°cone for N 2 , a small deviation can be seen for CO 2 .It is difficult to pinpoint the reason for this.Perhaps the behavior of the polyatomic vibratory excitation with the very large Knudsen numbers in the wake can simply no longer be reproduced with the model.
The comparison of the heat flux on the cone surface in Figure 20 shows very good agreement.

Computational performance
It is expected that the resolution requirements for the ESBGK method are less restrictive than for the DSMC method, depending on the Knudsen number, which allows for computational time savings as briefly described in the introduction.Thus, the basic idea is to apply the ESBGK method in regions of low Knudsen numbers to gain computational efficiency compared to the DSMC method.The aim is then to use this particle-based continuum method for a coupling with the DSMC method in order to solve multi-scale problems.In the following, the resolution requirements of both of the methods are discussed in greater detail.
The spatial resolution of the BGK particle method is essentially determined by the gradients in the flow and thus by the decoupling of the particle movement and the collision operator.However, this is not necessarily coupled to the mean free path as in the case of DSMC, which is why in denser cases, significantly fewer particles can be used in a simulation.To achieve an optimal spatial resolution for the BGK method, appropriately adapted grids would be needed, as described in Pfeiffer and Gorji [11].However, all simulations for both the ESBGK and the DSMC methods presented in this paper were performed on the same grid to be able to test the physical accuracy of the ESBGK model.Therefore, the computational optimum of both methods could not be achieved here.Nevertheless, regarding for example the simulations of Case 3, half as many particles could be used in the ESBGK case compared to DSMC, lead-     ing to a higher computational efficiency.This difference in the particle number would become even more pronounced in a 3D simulation case, since the particle number required to resolve the mean free path in the DSMC method increases with the power of the dimension, while the BGK method instead only needs to resolve the gradients in flow quantities.
Regarding temporal resolution, the mean collision time needs to be resolved for the DSMC method while for BGK, the stiff relaxation frequency needs to be resolved instead, at least when an explicit first-order time discretization is used.A few methods exist for particle BGK methods to circumvent this resolution criterion [19,21,20].Nevertheless, in the simulations here, a simple explicit time integration was used.However, the collision frequency in the DSMC method and the relaxation frequency in the BGK method differ.As described in Mathiaud et al. [1], both vary depending on the collision model used.The relaxation frequency of the ESBGK model scales with the Prandtl number, often allowing for larger time steps even with explicit time integrations, which in the end leads to a higher computational efficiency.Considering all the above-mentioned arguments, one concludes that the ESBGK method should become computationally more efficient than the DSMC method when solving smaller Knudsen numbers flows.This is also reflected in the simulations presented above, where in Case 1, the computational time could be reduced by a factor of 10 with the ESBGK method compared to the DSMC method, while in Case 3, due to the smaller Knudsen number, a factor of 28 could be saved.However, for larger Knudsen numbers, the DSMC method will be faster than the ESBGK method at some point, because only few or even no collisions need to be computed anymore.A more detailed comparison of computational times goes beyond the scope of this paper, as it mainly focuses on the theoretical foundations and physical results of the model.

Conclusion
An ellipsoidal statistical BGK mixture method is proposed and implemented in the open-source particle code PICLas based on the model of Mathiaud et al. [1], Pfeiffer [2] and Brull [3,4].Multi-species gas flows are simulated while allowing for the treatment of internal energies of di-and polyatomic molecules, aiming for a coupling with the DSMC method to solve multi-scale problems.As test cases for verification, different supersonic Couette flows and flows around a 70°blunted cone are used, where in general, overall good agreement between the ESBGK and DSMC methods is achieved for different type of mixtures and Knudsen numbers.Deviations occur concerning the onset of the temperature increase in the shock region in front of the 70°blunted cone, which is a known behavior for the ESBGK method.Visible deviations in the translational temperature profile in the shock region become smaller looking at smaller Knudsen number flows with narrower shock regions.The agreement of the simulation results in both the post-shock region and on the cone surface i.e. the heat flux are excellent.Minor deviations are visible in the wake of the 70°blunted cone for the vibrational temperature of CO 2 which is expected to be related to the polyatomic vibratory excitation in this high-Knudsen regime.Additionally, two methods for the determination of the transport coefficients of a gas mixture, Wilke's mixing rules and collision integrals, are compared, whereby the latter shows larger errors for mixtures with large mass ratios.A considerable reduction of the computational duration in transition and continuum regimes is possible with the ESBGK method compared to DSMC due to less restricted resolution requirements of the BGK method.
In future work, the implementation of chemical reactions into the model is envisioned.Furthermore, an extension of the already implemented atomic Shakhov model to polyatomic gas mixtures will be done.

Figure 1 :
Figure 1: Comparison of the stationary solution for the translational temperature in a supersonic Couette flow for a N 2 -He mixture.The 50 %-50 % mixture results are indicated with solid lines, the 80 %-20 % mixture with dashed lines and the 20 %-80 % mixture with dotted lines.

Figure 2 :
Figure 2: Comparison of the stationary solution for the rotational temperature in a supersonic Couette flow for a N 2 -He mixture.The 50 %-50 % mixture results are indicated with solid lines, the 80 %-20 % mixture with dashed lines and the 20 %-80 % mixture with dotted lines.

Figure 3 :
Figure 3: Comparison of the stationary solution for the temperature in a supersonic Couette flow for a 50 %-50 % N 2 -N mixture with Kn VHS ≈ 0.0121 indicated with solid lines, Kn VHS ≈ 0.121 with dashed lines and Kn VHS ≈ 1.21 with dotted lines.

Figure 4 :Figure 5 :Figure 6 :
Figure 4: Comparison of the stationary solution for the temperature in a supersonic Couette flow for a 50 %-50 % CO 2 -N 2 mixture.

Figure 7 :
Figure 7: Geometry of the 70 • blunted cone.R b = 25.0 mm, R c = 1.25 mm, R j = 2.08 mm, R n = 12.5 mm, R s = 6.25 mm.S denotes the arc length along the surface.

s − 1 Figure 8 :
Figure 8: 70°blunted cone, Case 1: Mixture mean values of velocity in x direction (dashed) and number density (solid) along stagnation stream line using DSMC and ESBGK.

Figure 9 :
Figure 9: 70°blunted cone, Case 1: Mixture mean values of the translation, rotational and vibrational temperatures along the stagnation stream line and speciesspecific values of the vibrational temperature at y = 20 mm using DSMC and ESBGK.

Figure 10 :Figure 11 :
Figure 10: 70°blunted cone, Case 1: Comparison of the mixture mean value of the translational temperature in the flow field using DSMC and ESBGK.Characteristic points A-D on the cone surface are indicated.

Figure 12 :
Figure 12: 70°blunted cone, Case 2: Mixture mean values of velocity in x direction (dashed) and number density (solid) along stagnation stream line using DSMC and ESBGK.
Figure 13: 70°blunted cone, Case 2: Mixture mean values of the translation, rotational and vibrational temperatures along the stagnation stream line and speciesspecific values of the vibrational temperature at y = 20 mm using DSMC and ESBGK.

s − 1 Figure 15 :
Figure 15: 70°blunted cone, Case 3: Mixture mean values of velocity in x direction (dashed) and number density (solid) along stagnation stream line using DSMC and ESBGK.
Figure 16: 70°blunted cone, Case 3: Mixture mean values of the translation, rotational and vibrational temperatures along the stagnation stream line and speciesspecific values of the vibrational temperature at y = 20 mm using DSMC and ESBGK.

s − 1 Figure 18 :
Figure 18: 70°blunted cone, Case 4: Mixture mean values of velocity in x direction (dashed) and number density (solid) along stagnation stream line using DSMC and ESBGK.

2 (
Figure 19: 70°blunted cone, Case 4: Mixture mean values of the translation, rotational and vibrational temperatures along the stagnation stream line and speciesspecific values of the vibrational temperature at y = 20 mm using DSMC and ESBGK.

Table 2 :
Overview of Couette flow test cases.

Table 3 :
Inflow conditions of hypersonic flow around a 70°blunted cone.