Kinetic Fokker–Planck (FP) methods for modeling rarefied gas flows have received increasing attention over the last few years. However, formulating such models for realistic multi-species gases is still an open subject of research. Therefore, in this letter, we develop a kinetic FP model for describing gas mixtures with particles interacting according to the variable hard-sphere interaction potential. In accordance with the kinetic FP framework, a stochastic solution algorithm is employed in order to solve the model on a particle level. Different test cases are carried out, and the performance of the proposed method is compared with the direct simulation Monte Carlo algorithm.

Aerodynamic analysis of a spacecraft requires modeling of hypersonic, rarefied gas flows. Because of high Mach and Knudsen numbers, such flows are dominated by strong non-equilibrium phenomena, for example, shock waves and expansion fans. Since classical Navier–Stokes solvers assume only slight deviations from local equilibrium, they are not applicable for such cases.

Instead, such flows can be described in general by the highly accurate Boltzmann equation.1 The Boltzmann equation models the evolution of a particle distribution function in phase space and allows describing even strong non-equilibrium effects. However, due to the high dimensionality of the phase space and the complexity of the Boltzmann collision operator, its direct numerical solution becomes a complex task. An alternative approach to model non-equilibrium flows is the direct simulation Monte Carlo (DSMC) method.2 DSMC directly simulates the stochastic motion of representative particles, while macroscopic quantities, like density or temperature, are calculated by averaging the particle ensemble. DSMC is widely used to model rarefied gas flows, but its computational effort becomes significant for high pressure gases, so it is not suitable for handling flows that feature locally small Knudsen numbers. For these reasons, it is the subject of research to find approximation methods for the Boltzmann equation and the DSMC method.

An example features the recently proposed kinetic Fokker–Planck (FP) method.3–6 This method approximates the collision integral of the Boltzmann equation by an FP operator in velocity space. Instead of solving the resulting FP equation directly, an equivalent stochastic simulation algorithm is used to model the gas at a microscopic level. Since no collisions have to be calculated in the kinetic FP method, the computational effort of this method is independent of the Knudsen number. Therefore, kinetic FP has the potential to be more efficient than DSMC for small Knudsen flows. Pioneering work concerning kinetic FP has been carried out by Jenny, Torrilhon, and Heinz3 who introduced a first model for describing single species gases and by Gorji, Torrilhon, and Jenny4 who extended this approach in order to fix a wrong prediction of heat fluxes. Since then, various kinetic FP models have been developed for describing complex gas flows,7–10 and the method has been applied to various test cases.11–16 

Recently, Hepp, Grabe, and Hannemann10 constructed a kinetic FP model for describing gas mixtures, assuming particle interaction according to the hard-sphere (HS) collision potential. In this letter, we use a similar approach to construct a kinetic FP model, assuming the variable hard-sphere17 (VHS) collision potential.

Note that the VHS collision model was originally developed for the DSMC method.17 Because the VHS collision model can represent transport properties of many gases better than the HS collision model,18 it has become the standard model for describing non-equilibrium gas flows.2,19,20 The development of a kinetic FP model based on the VHS collision potential is therefore crucial in order to make the kinetic FP framework applicable for real engineering applications.

The first part of this letter presents the construction of the FP operator. The second part describes a stochastic solution algorithm that solves the previously described FP equation on a microscopic level. In the last part, different test cases are examined in order to check the performance of our model against the DSMC algorithm.

The construction of our model is similar to that of the HS FP model in Ref. 10. The statistical state of a multispecies gas is described by a set of species specific distribution functions f(α)v(α),x,t, representing the number of particles of species α, that can be found in small volume dx around position x, with velocities in a small range dc(α) around velocity c(α) at time t. These distribution functions are described by a set of Fokker–Planck equations,

f(α)t+vi(α)f(α)xi+Fi(α)m(α)f(α)vi(α)  =f(α)Ai(α)vi(α)+2f(α)D(α)vi(α)vi(α)SFP(α),
(1)

in velocity space. Here, F(α) refers to an external force, m(α) denotes the particle mass, t is the time, and SFP(α) is the FP operator. The species specific drift coefficients A(α) and the diffusion coefficients D(α) are model parameters, and they are constructed so that the Fokker–Planck operator reproduces lower order Boltzmann production terms,

PBoltz(α)X=!PFP(α)X,Xci(α),ci(α)ci(α),c<i(α)cj>(α)ci(α)cj(α)cj(α),
(2)

in the limit of small Knudsen numbers. In the above-mentioned equation, the quantity c(α) = v(α)u(α) denotes the thermal particle velocity of species α with respect to the species flow velocity u(α), while the production term PΨ(α) with respect to a collision operator SΨ(α) is defined as

PΨ(α)χ=SΨ(α)χdc.
(3)

Note that Eq. (2) leads the FP model to conserve mass, momentum, and energy and to correctly predict transport coefficients for the mixture in the continuum limit.21,22

In order to construct an FP operator that satisfies Eq. (2), production terms must be evaluated. Production terms for the FP operator can be calculated independently of a distribution function.10 Relevant results are

PFP(α)ci(α)=m(α)Ai(α)f(α),
(4)
PFP(α)c(α)c(α)=2m(α)Ai(α)ci(α)f(α)+6ρ(α)D(α),
(5)
PFP(α)ci(α)cj(α)=2δijρ(α)D(α)+m(α)Ai(α)cj(α)f(α)+m(α)Aj(α)ci(α)f(α),
(6)
PFP(α)ci(α)c(α)c(α)=m(α)Ai(α)cj(α)cj(α)f(α)+2m(α)Aj(α)cj(α)ci(α)f(α).
(7)

Here, ρ(α) denotes the mass density of species α, and g|f(α)=f(α)x,v(α),tgv(α)dv(α) denotes a velocity moment with respect to the distribution function f(α).

Production terms for the Boltzmann collision operator must be calculated with respect to the VHS collision potential since our model is intended to describe particle interaction according to the VHS potential. Such production terms can be extracted, for example, from the work of Gupta and Torrilhon23 who derived Boltzmann production terms for binary gas mixtures assuming a general interaction potential and a distribution function according to Grad’s 13 moment method. These results can be easily adapted to the VHS interaction potential and mixtures with any number of species. The calculations are straightforward, and the details can be found in the supplementary material.

The final results for lower order Boltzmann production terms read as

PBoltz(α)ci(α)=β=1Nν|VHS(αβ)μ(β)53ρ(α)VHS[1]ui(α)ui(β)+16θ^(αβ)VHS[2]ĥi(α)ρ(α)ρ(β)ĥi(β),
(8)
P(α)c(α)c(α)=10β=1Nν(αβ)μ(β)ρ(α)VHS[1]×kBm(α)+m(β)T(α)T(β)μ(β)3u(α)u(β)2.
(9)

Here, N denotes the number of species in the mixture, and the quantities μ(α), θ^(α), θ^(αβ), Δθ^(αβ), h^(α), and ν(αβ) are defined as

μ(α)m(α)m(α)+m(β),
(10)
θ^(α)kBT^(α)m(α),
(11)
θ^(αβ)θ^(α)+θ^(β)2,
(12)
Δθ^(αβ)μ(α)θ^(α)μ(β)θ^(β)θ^(αβ),
(13)
h^(α)q^(α)52ρ(α)θ^(α)ud(α),
(14)
ν(αβ)=165VHS[6]πn(β)dref(αβ)2θ^(αβ).
(15)

In the above-mentioned equations, the species flow velocity is given by

u(α)x,t=1ρ(α)x,tm(α)v(α)|f(α),
(16)

the mixture velocity is given by

umixx,t=1ρα=1Nu(α)x,tρ(α),
(17)

with the mass density ρ(α) = m(α) · n(α), and the diffusion velocity is given by

udx,t=umixx,tu(α)x,t.
(18)

Note, that this letter distinguishes between thermal particle velocities c^(α)=v(α)umix and c(α) = v(α)u(α). For instance, the temperature T^(α) is given by

T^(α)x,t=m(α)3kBρ(α)m(α)c^(α)c^(α)|f(α),
(19)

and the heat flux is given by

q^(α)x,t=m(α)c^(α)c^(α)c^(α)|f(α).
(20)

The temperature T(α) is defined similar to that in Eq. (19), where velocities c^(α) are replaced by velocities c(α).

The VHS-scaling parameters are defined as

VHS[1]=3Γ3ν(αβ)Γ4ν(αβ),
(21)
VHS[2]=615Γ3ν(αβ)Γ4ν(αβ),
(22)

where Γ denotes the gamma function. Note, that the velocity exponent ν(αβ) and the reference diameter dref(αβ) indicate species specific model parameters related to the VHS collision potential.17 

It is worth noting that Boltzmann production terms for the VHS interaction potential (8) and (9) feature the same structure as the corresponding production terms for the HS collision potential.24 This is because the HS collision potential can be expressed as a special case of the VHS collision potential. In this case, the velocity exponent becomes ν(αβ) = 0, the parameters (21) and (22) become one, and VHS production terms (8) and (9) reduce to HS production terms, as employed in Ref. 10.

Similar to the HS FP model,10 first, an FP operator with a drift coefficient that is linear in the thermal particle velocities and that obeys requirement (2) for X=ci(α) and X=ci(α)ci(α) is constructed. Therefore, drift and diffusion coefficients are chosen so that FP production terms (4) and (5) reproduce corresponding VHS Boltzmann production terms (8) and (9). Afterward, the drift coefficient is extended by a second order expression in order to obey the full requirement (2). As described above, VHS production terms (8) and (9) feature the same structure as the corresponding HS production terms. Therefore, drift and diffusion coefficients of our model can be derived similar to the HS FP model described in Ref. 10. For brevity, this derivation is not repeated. Instead, only the final results are presented.

For the diffusion coefficient,

D(α)2=103VHS[1]β=1Nν(αβ)μ(β)μ(β)×kBT(α)m(α)kBT(β)m(β)+13u(α)u(β)2
(23)

can be found. For the drift coefficient, the ansatz

Ai(α)s(α)ci(α)+Ki(α)+ψij(α)cj(α)+γi(α)cj(α)cj(α)cj(α)cj(α)|f(α)
(24)

is made. The frequencies s(α) are given by

s(α)=53VHS[1]β=1Nν(αβ)μ(β),
(25)

while the constants Ki(α) are defined as

Ki(α)=β=1Nν(αβ)μ(β)53VHS[1]ui(α)ui(β)+16θ^(α)VHS[2]1ρ(α)ĥi(α)1ρ(β)ĥi(β).
(26)

The nine model parameters ψij(α) and γi(α) in Eq. (24) are chosen such that the equality,

PFP(α)ci(α)cj(α)=!PBoltz(α)c<i(α)cj>(α)+13δijPBoltz(α)cl(α)cl(α),
(27)
  PFP(α)ci(α)c(α)c(α)=!PBoltz(α)ci(α)c(α)c(α),
(28)

holds. Production terms PBoltz(α)c<i(α)cj>(α),PBoltz(α)cl(α)cl(α) and PBoltz(α)ci(α)c(α)c(α) are given by Eqs. (27)–(29) in the supplementary material. Equations (27) and (28) can be transformed together with Eqs. (24), (6), and (7) to a system of nine linear equations that have to be solved in order to determine the nine model parameters ψij(α) and γi(α). This system is presented in the supplementary material. It is noteworthy here that the drift and diffusion coefficients of our model feature a similar structure to the drift and diffusion coefficients of the HS FP model.10 In particular, for the special case of HS molecules, our model reduces to the HS FP model. This is due to the similarity between the HS and VHS production terms, as described above.

A main motivation for introducing our kinetic FP model is to efficiently simulate non-equilibrium gas flows. Therefore, in the following, we will discuss the applicability of our model to high Knudsen number flows.

Because of condition (2), our model captures the correct continuum limit. However, for the moderate Knudsen number regime, higher order production terms, which are not included in condition (2), become relevant. In combination with the polynomial approach for the drift coefficient (24), this might cause our model to become invalid for simulating non-equilibrium flows. On the other hand, the FP Eq. (1) also reproduces the collision-less Boltzmann equation in the limit of high Knudsen numbers. This is because the drift coefficient (24) and the diffusion coefficient (23) include the collision frequencies ν(αβ) in every term, which vanish with decreasing particle density. This behavior indicates that the model can capture at least some rarefaction effects, for instance, free molecular transport at low particle densities. A similar behavior has also been reported in the open literature. For instance, the HS FP model,10 which uses a similar approach to the model presented in this study, has been successfully applied to a Couette flow up to a Knudsen number of 0.5.

The authors, therefore, believe that the model presented above can also be applied to the moderate Knudsen number regime. However, further studies, that go beyond the scope of this work, are necessary to check this assumption.

The high dimension of the FP Eq. (1) makes its direct solution complex. Therefore, the FP Eq. (1) is not solved directly. Instead, a kinetic particle scheme is employed to simulate the corresponding microscopic particle motion. Therefore, the distribution functions f(α) are approximated by an ensemble of computational particles. The motion of particles is described by a set of stochastic processes, consistent with the FP Eq. (1), while macroscopic moments, e.g., flow velocity or temperature, are calculated by averaging the particle ensemble. Note, that this approach is equivalent to directly solving Eq. (1) in the limit of an infinite number of simulated particles.3,10 For more details on the kinetic FP framework, please refer to Refs. 4, 9, 10, and 25.

If external forces are neglected, the particle path can be described in accordance with the FP model described above by the following stochastic equations of motion:3,4,10

dVi(α)dt=Ai(α)+2D(α)dWi(α)dt,
(29)
dXidt=Vi(α).
(30)

Here, X and V(α) represent particle velocity and position, respectively, dW(α) refers to a Wiener process with zero expectation, and dWi(α)dWj(α)=δij. Note that X and V(α) denote random variables that must be distinguished from phase space coordinates x and v(α).

In order to solve systems (29) and (30), the scheme derived in Ref. 10 is applied. In particular, particle velocities and positions are updated as

Vi(α),n+1=ϵCi(α),n+1+ui(α),n+1,
(31)
Xin+1=Xin+VinΔt.
(32)

The index n refers to the time step, and Δt denotes the time step size. Thermal particle velocities C(α) = V(α)u(α) are defined with respect to the species flow velocity and are updated as

Ci(α),n+1=Ci(α),nexps|VHS(α),nΔt+ΔtNi(α)C(α),n+D(α),ns(α),n1exp2s(α),nΔtξi(α),
(33)

while the scaling factor ϵ is given by

ϵ=3kBT(α),n+1/m(α)Ci(α),n+1Ci(α),n+1/n(α).
(34)

In the above-mentioned equations, ξi(α) denotes independent standard normal variates, Ci(α),n+1Ci(α),n+1 means an ensemble average, and the parameter Ni(α) is defined as

Ni(α)C(α),n=ψij(α)Cj(α),n+γi(α)×Cj(α),nCj(α),nCj(α),nCj(α),n|f(α).
(35)

Similar to the HS FP model,10 the temperatures T(α),n+1 and flow velocities ui(α),n+1 can be extracted as a solution of a system of differential equations. This system is described in more detail in the supplementary material. In summary, particle positions and velocities are updated as described in Algorithm 1.

ALGORITHM 1

: Update of particle positions and velocities.

1. Evaluate the required statistical moments for every grid cell. 
2. Solve the linear system that is described by Eqs. (31) and (32) in the supplementary material. This yields the model parameters ψij(α) and γi(α)
3. Use the model parameters ψij(α) and γi(α) to calculate the new thermal particle velocities Ci(α),n+1 using Eq. (33)
4. Solve the system of differential equations that is described by Eq. (33) in the supplementary material. This yields the temperatures T(α),n+1 and flow velocities ui(α),n+1 at the end of the integration step. 
5. Use the temperatures T(α),n+1, the flow velocities ui(α),n+1, and the thermal particle velocities Ci(α),n+1 to calculate the final particle positions and velocities by Eqs. (31) and (32)
1. Evaluate the required statistical moments for every grid cell. 
2. Solve the linear system that is described by Eqs. (31) and (32) in the supplementary material. This yields the model parameters ψij(α) and γi(α)
3. Use the model parameters ψij(α) and γi(α) to calculate the new thermal particle velocities Ci(α),n+1 using Eq. (33)
4. Solve the system of differential equations that is described by Eq. (33) in the supplementary material. This yields the temperatures T(α),n+1 and flow velocities ui(α),n+1 at the end of the integration step. 
5. Use the temperatures T(α),n+1, the flow velocities ui(α),n+1, and the thermal particle velocities Ci(α),n+1 to calculate the final particle positions and velocities by Eqs. (31) and (32)

The proposed kinetic FP model is applied to different test cases in order to check its performance. Simulations are performed for a He–Ar mixture. The associated collision parameters are listed in Table I. For all test cases, reference DSMC simulations are performed, assuming the same VHS collision model. For additional comparison, results of simulations assuming an equivalent HS collision model are also shown. Therefore, collision parameters similar to the ones listed in Tab. I are used, but the viscosity exponent is set to ω(ij) = 0.5. To avoid different results due to numerical discretization, the same spatial and temporal resolution is used for the DSMC and the Fokker–Planck simulations. Detailed information on the simulation settings can be found in Table II.

TABLE I.

Interspecies VHS model parameter.2 Velocity and viscosity exponents are related by ω(αβ) = ν(αβ) + 1/2.

Species αSpecies βdref(αβ)(1010m)ω(αβ)Tref(αβ)
He He 2.33 0.66 273.15 
He Ar 3.25 0.735 273.15 
Ar Ar 4.17 0.81 273.15 
Species αSpecies βdref(αβ)(1010m)ω(αβ)Tref(αβ)
He He 2.33 0.66 273.15 
He Ar 3.25 0.735 273.15 
Ar Ar 4.17 0.81 273.15 
TABLE II.

Number of simulation particles Np, number of cells Nc, time step size Δt, and reference number density n0 for both simulation cases. The reference number density n0 is used to normalize simulation results given in Figs. 2 and 3. Note that the diffusion test case allows using a smaller number of simulation particles because it features a larger signal-to-noise ratio.

Simulation caseNpNcΔt (s)n0 (1/m3)
Supersonic Couette flow 76 000 100 10–6 3 · 1019 
Diffusion flow 20 000 100 2 · 10−6 6 · 1019 
Simulation caseNpNcΔt (s)n0 (1/m3)
Supersonic Couette flow 76 000 100 10–6 3 · 1019 
Diffusion flow 20 000 100 2 · 10−6 6 · 1019 

All simulations are performed using the SPARTA code26 that has been extended by the proposed kinetic Fokker–Planck model.

To check the performance of the model in predicting shear stresses and heat fluxes, a supersonic Couette flow is investigated. The left side of Fig. 1 illustrates the numerical setup. Domain boundaries in the y- and z-direction are assumed to be periodic while boundaries in the x-direction are modeled as fully diffusive walls with a temperature of Tw = 300 K. The higher wall moves in the y-direction with a velocity of vw = 1000 m/s while the lower wall is stationary. Both walls are separated by a distance of 1 m, which is also used as a reference length for defining the Knudsen number. Simulations are performed for a Knudsen number of Kn = 0.05, assuming a number fraction of χf(He)=n(He)/n(He)+n(Ar)=0.5. The spatial and temporal discretization is set to resolve the particle mean free path and mean collision time.

FIG. 1.

Numerical setup for the test cases that are investigated. Left side: Supersonic Couette flow. Right side: Diffusion flow.

FIG. 1.

Numerical setup for the test cases that are investigated. Left side: Supersonic Couette flow. Right side: Diffusion flow.

Close modal

The upper left side of Fig. 2 shows the number density distribution plots. As described in Ref. 10, thermodiffusion induces species separation. Compared to the HS model, the separation effect is slightly smaller for the VHS model, indicating a smaller thermal-diffusion ratio.1 The upper right side of Fig. 2 shows the temperature distributions. For both collision models, a significant temperature slip occurs at walls. In addition, species temperatures are separated, indicating a strong degree of thermal non-equilibrium. The lower side of Fig. 2 shows shear stress and heat flux distributions along the simulation domain. While the distributions are nearly constant and linear in the bulk flow, they adopt a non-constant and non-linear shape at the vicinity of the walls, indicating the boundary layer. As expected, the collision model strongly influences shear stress and heat flux distributions.

FIG. 2.

Species density, temperature, shear stress, and heat flux distributions for a hypersonic Couette flow. The reference number density n0 is given in Table II. Lines: DSMC/VHS results. Circles: FP/VHS results. Dashed lines: FP/HS results. Black: He. Red: Ar.

FIG. 2.

Species density, temperature, shear stress, and heat flux distributions for a hypersonic Couette flow. The reference number density n0 is given in Table II. Lines: DSMC/VHS results. Circles: FP/VHS results. Dashed lines: FP/HS results. Black: He. Red: Ar.

Close modal

In general, very good agreement between the results of our FP model and reference DSMC simulations can be found. Separation effects in the species density and temperature distributions as well as shear stresses and heat fluxes are correctly predicted.

To check the capability of the model in prediction of diffusion effects, a diffusion test case is investigated. The setup is illustrated at the right side of Fig. 1. The lower x-boundary is modeled as the reservoir of Ar particles while the higher x-boundary is modeled as the reservoir of He particles. For both reservoirs, the same temperature T0 = 300 K is assumend. The distance between the x-boundaries is set to 1 m. This length scale is also used as the reference length to define the Knudsen number. The number densities in the reservoir are set corresponding to a Knudsen number of Kn = 0.05, assuming a constant density fraction χf(He)=χf(Ar)=0.5. The spatial and temporal discretization is set to resolve the particle mean free path and the mean collision time.

The left side of Fig. 3 shows the species density distribution plots, while the right side of Fig. 3 shows species velocity distributions. Interestingly, no significant differences can be found between the HS and VHS collision model, indicating a weak dependence of the diffusion coefficient on the VHS velocity exponent. In general, our FP model accurately reproduces DSMC results.

FIG. 3.

Species density and velocity distributions for a diffusion flow. The reference number density n0 is given in Table II. Lines: DSMC/VHS results. Circles: FP/VHS results. Dashed lines: FP/HS results.

FIG. 3.

Species density and velocity distributions for a diffusion flow. The reference number density n0 is given in Table II. Lines: DSMC/VHS results. Circles: FP/VHS results. Dashed lines: FP/HS results.

Close modal

The main motivation for introducing our kinetic FP model is to reduce the computational effort compared to a conventional DSMC solver. Therefore, in the following, we discuss the computational efficiency of our kinetic FP method compared to the DSMC algorithm.

As a main advantage over DSMC, the kinetic FP model does not have to resolve molecular length scales. Therefore, kinetic FP has the potential to be more efficient than DSMC in the small Knudsen number regime since simulations can be carried out with a smaller number of cells and particles and a larger time step than required for the DSMC algorithm. In general, this efficiency gain increases as the spatial and temporal numerical resolution for the FP simulation increases compared to molecular length scales.

On the other hand, the efficiency gain of the FP model decreases when the spatial and temporal resolution becomes small compared to molecular length scales. At some point, when using a very small numerical discretization for kinetic FP, DSMC will become more efficient. This is because the computing time per time step and cell for modeling DSMC collisions decreases with a decreasing time step and cell size, while the computing time per time step and cell remains constant for the FP algorithm.

A similar behavior can also be observed when considering the simulation cases discussed above. Table III shows the required computing times. The DSMC algorithm is clearly more efficient in both cases. However, it should be noted that the FP simulations use the same numerical resolution as that used for the DSMC simulations. For instance, the cell size Δx is chosen to be five times smaller than the average mean free path λ. As described above, our model cannot be expected to be more efficient than DSMC in such a case. On the other hand, the numerical resolution for the kinetic FP simulations was chosen to be, only for the sake of simplicity, as small as that in the DSMC simulations and is clearly not representative for a real application. However, a full efficiency study, which takes into account the effect of numerical resolution, should be discussed for more complex test cases and goes beyond the scope of this work.

TABLE III.

Computational times for DSMC and kinetic FP simulations. For both test cases, times are normalized to the computational time required for the DSMC simulation.

Simulation caseTime DSMCTime FP
Supersonic Couette flow 5.0 
Diffusion flow 8.3 
Simulation caseTime DSMCTime FP
Supersonic Couette flow 5.0 
Diffusion flow 8.3 

In summary, we developed a kinetic FP model for describing monatomic gas mixtures with particles interacting according to the VHS collision potential. The model has been formulated for any number of species and has been tested for a binary He–Ar mixture. It is worth mentioning here that our model includes a recently proposed kinetic FP model for HS gas mixtures10 as a special case. We investigated representative test cases and checked the capability of our model to describe non-equilibrium flows.

In future work, the application of our model to more complex test cases is planned. In addition, it is worth noting that a framework similar to that presented in this letter could be used to construct a kinetic FP model for any other collision potential, provided that appropriate production terms can be calculated. For instance, as an extension of the VHS model, the VSS27 or M128 collision potential could be integrated into future FP models.

See the supplementary material for the derivation of VHS production terms, the linear system of equations that is solved in point 2 of Algorithm 1, and the system of differential equations that is solved in point 4 of Algorithm 1.

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

1.
G. M.
Kremer
,
An Introduction to the Boltzmann Equation and Transport Processes in Gases
(
Springer Berlin Heidelberg
,
2010
).
2.
G. A.
Bird
,
Molecular Gas Dynamics and the Direct Simulation of Gas Flows (Oxford Engineering Science Series)
(
Clarendon Press
,
1994
).
3.
P.
Jenny
,
M.
Torrilhon
, and
S.
Heinz
, “
A solution algorithm for the fluid dynamic equations based on a stochastic model for molecular motion
,”
J. Comput. Phys.
229
,
1077
1098
(
2010
).
4.
M. H.
Gorji
,
M.
Torrilhon
, and
P.
Jenny
, “
Fokker–Planck model for computational studies of monatomic rarefied gas flows
,”
J. Fluid Mech.
680
,
574
601
(
2011
).
5.
J.
Mathiaud
and
L.
Mieussens
, “
A Fokker–Planck model of the Boltzmann equation with correct Prandtl number
,”
J. Stat. Phys.
162
,
397
414
(
2015
).
6.
S. K.
Singh
and
S.
Ansumali
, “
Fokker–Planck model of hydrodynamics
,”
Phys. Rev. E
91
,
033303
(
2015
).
7.
M. H.
Gorji
and
P.
Jenny
, “
A Fokker–Planck based kinetic model for diatomic rarefied gas flows
,”
Phys. Fluids
25
,
062002
(
2013
).
8.
J.
Mathiaud
and
L.
Mieussens
, “
A Fokker–Planck model of the Boltzmann equation with correct Prandtl number for polyatomic gases
,”
J. Stat. Phys.
168
,
1031
1055
(
2017
).
9.
H.
Gorji
and
P.
Jenny
, “
A kinetic model for gas mixtures based on a Fokker–Planck equation
,”
J. Phys.: Conf. Ser.
362
,
012042
(
2012
).
10.
C.
Hepp
,
M.
Grabe
, and
K.
Hannemann
, “
A kinetic Fokker–Planck approach to model hard-sphere gas mixtures
,”
Phys. Fluids
32
,
027103
(
2020
).
11.
M.
Pfeiffer
,
A.
Mirza
, and
P.
Nizenkov
, “
Evaluation of particle-based continuum methods for a coupling with the direct simulation Monte Carlo method based on a nozzle expansion
,”
Phys. Fluids
31
,
073601
(
2019
).
12.
M.
Pfeiffer
and
M. H.
Gorji
, “
Adaptive particle-cell algorithm for Fokker–Planck based rarefied gas flow simulations
,”
Comput. Phys. Commun.
213
,
1
8
(
2017
).
13.
E.
Jun
,
M. H.
Gorji
,
M.
Grabe
, and
K.
Hannemann
, “
Assessment of the cubic Fokker–Planck-DSMC hybrid method for hypersonic rarefied flows past a cylinder
,”
Comput. Fluids
168
,
1
13
(
2018
).
14.
E.
Jun
, “
Cubic Fokker–Planck-DSMC hybrid method for diatomic rarefied gas flow through a slit and an orifice
,”
Vacuum
159
,
125
133
(
2019
).
15.
Y.
Jiang
,
Z.
Gao
, and
C.-H.
Lee
, “
Particle simulation of nonequilibrium gas flows based on ellipsoidal statistical Fokker–Planck model
,”
Comput. Fluids
170
,
106
120
(
2018
).
16.
B.
Xiao
and
R.
Li
, “
Ellipsoidal statistical Fokker–Planck simulation of thermally induced gas flow in a ratchet microchannel
,”
J. Phys.: Conf. Ser.
1168
,
032037
(
2019
).
17.
G. A.
Bird
, “
Monte-Carlo simulation in an engineering context
,”
Prog. Astronaut. Aeronaut.
74
,
239
255
(
1981
).
18.
G. A.
Bird
, “
Recent advances and current challenges for DSMC
,”
Comput. Math. Appl.
35
,
1
14
(
1998
).
19.
R.
Prakash
,
S.
Gai
, and
S.
O’Byrne
, “
A direct simulation Monte Carlo study of hypersonic leading-edge separation with rarefaction effects
,”
Phys. Fluids
30
,
063602
(
2018
).
20.
M. S.
Ivanov
,
A. V.
Kashkovsky
,
S. F.
Gimelshein
,
G. N.
Markelov
,
A. A.
Alexeenko
,
Y. A.
Bondar
,
G. A.
Zhukova
,
S.
Nikiforov
, and
P. V.
Vaschenkov
, “
SMILE system for 2D/3D DSMC computations
,” in
Proceedings of 25th International Symposium on Rarefied Gas Dynamics
(
Siberian Branch of the Russian Academy of Sciences
,
Saint Petersburg, Russia
,
2006
), pp.
21
28
, available at www.worldcat.org/title/rarefied-gas-dynamics-proceedings-of-25th-international-symposium-on-rarefied-gas-dynamics-saint-petersburg-russia-july-21-18-2006-monograph/oclc/862147291
21.
H.
Struchtrup
,
Macroscopic Transport Equations for Rarefied Gas Flows
(
Springer Berlin Heidelberg
,
2005
), pp.
145
160
.
22.
V. K.
Gupta
,
H.
Struchtrup
, and
M.
Torrilhon
, “
Regularized moment equations for binary gas mixtures: Derivation and linear analysis
,”
Phys. Fluids
28
,
042003
(
2016
).
23.
V. K.
Gupta
and
M.
Torrilhon
, “
Automated Boltzmann collision integrals for moment equations
,”
AIP Conf. Proc.
1501
,
67
74
(
2012
).
24.
V. K.
Gupta
, “
Mathematical modeling of rarefied gas-mixtures
,” Ph.D. thesis,
RWTH Aachen University
,
2015
.
25.
P.
Jenny
,
S.
Küchlin
, and
H.
Gorji
, “
Controlling the bias error of Fokker–Planck methods for rarefied gas dynamics simulations
,”
Phys. Fluids
31
,
062005
(
2019
).
26.
M. A.
Gallis
,
J. R.
Torczynski
,
S. J.
Plimpton
,
D. J.
Rader
, and
T.
Koehler
, “
Direct simulation Monte Carlo: The quest for speed
,”
AIP Conf. Proc.
1628
,
27
36
(
2014
).
27.
K.
Koura
and
H.
Matsumoto
, “
Variable soft sphere molecular model for inverse-power-law or Lennard-Jones potential
,”
Phys. Fluids A
3
,
2459
2465
(
1991
).
28.
A.
Kersch
,
W.
Morokoff
, and
C.
Werner
, “
Selfconsistent simulation of sputter deposition with the Monte Carlo method
,”
J. Appl. Phys.
75
,
2278
2285
(
1994
).
All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Supplementary Material