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 Heinz^{3} who introduced a first model for describing single species gases and by Gorji, Torrilhon, and Jenny^{4} 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 Hannemann^{10} 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-sphere^{17} (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(\alpha )v(\alpha ),x,t$, representing the number of particles of species *α*, that can be found in small volume d**x** around position **x**, with velocities in a small range d**c**^{(α)} around velocity **c**^{(α)} at time *t*. These distribution functions are described by a set of Fokker–Planck equations,

in velocity space. Here, **F**^{(α)} refers to an external force, *m*^{(α)} denotes the particle mass, *t* is the time, and $SFP(\alpha )$ 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,

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\Psi (\alpha )$ with respect to a collision operator $S\Psi (\alpha )$ is defined as

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

Here, *ρ*^{(α)} denotes the mass density of species *α*, and $g|f(\alpha )=\u222bf(\alpha )x,v(\alpha ),tgv(\alpha )dv(\alpha )$ 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 Torrilhon^{23} 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

Here, *N* denotes the number of species in the mixture, and the quantities *μ*^{(α)}, $\theta ^(\alpha )$, $\theta ^(\alpha \beta )$, $\Delta \theta ^(\alpha \beta )$, $h^(\alpha )$, and *ν*^{(αβ)} are defined as

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

the mixture velocity is given by

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

Note, that this letter distinguishes between thermal particle velocities $c^(\alpha )=v(\alpha )\u2212umix$ and **c**^{(α)} = **v**^{(α)} − **u**^{(α)}. For instance, the temperature $T^(\alpha )$ is given by

and the heat flux is given by

The temperature *T*^{(α)} is defined similar to that in Eq. (19), where velocities $c^(\alpha )$ are replaced by velocities **c**^{(α)}.

The VHS-scaling parameters are defined as

where Γ denotes the gamma function. Note, that the velocity exponent *ν*^{(αβ)} and the reference diameter $dref(\alpha \beta )$ 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(\alpha )$ and $X=ci(\alpha )ci(\alpha )$ 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,

can be found. For the drift coefficient, the ansatz

is made. The frequencies *s*^{(α)} are given by

while the constants $Ki(\alpha )$ are defined as

The nine model parameters $\psi ij(\alpha )$ and $\gamma i(\alpha )$ in Eq. (24) are chosen such that the equality,

holds. Production terms $PBoltz(\alpha )c<i(\alpha )cj>(\alpha ),PBoltz(\alpha )cl(\alpha )cl(\alpha )$ and $PBoltz(\alpha )ci(\alpha )c(\alpha )c(\alpha )$ 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 $\psi ij(\alpha )$ and $\gamma i(\alpha )$. 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}

Here, **X** and **V**^{(α)} represent particle velocity and position, respectively, d**W**^{(α)} refers to a Wiener process with zero expectation, and $dWi(\alpha )dWj(\alpha )=\delta 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

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

while the scaling factor *ϵ* is given by

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

Similar to the HS FP model,^{10} the temperatures *T*^{(α),n+1} and flow velocities $ui(\alpha ),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.

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 $\psi ij(\alpha )$ and $\gamma i(\alpha )$. |

3. Use the model parameters $\psi ij(\alpha )$ and $\gamma i(\alpha )$ to calculate the new thermal particle velocities $Ci(\alpha ),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(\alpha ),n+1$ at the end of the integration step. |

5. Use the temperatures T^{(α),n+1}, the flow velocities $ui(\alpha ),n+1$, and the thermal particle velocities $Ci(\alpha ),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 $\psi ij(\alpha )$ and $\gamma i(\alpha )$. |

3. Use the model parameters $\psi ij(\alpha )$ and $\gamma i(\alpha )$ to calculate the new thermal particle velocities $Ci(\alpha ),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(\alpha ),n+1$ at the end of the integration step. |

5. Use the temperatures T^{(α),n+1}, the flow velocities $ui(\alpha ),n+1$, and the thermal particle velocities $Ci(\alpha ),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.

Species α
. | Species β
. | $dref(\alpha \beta )(10\u221210m)$ . | ω^{(αβ)}
. | $Tref(\alpha \beta )$ . |
---|---|---|---|---|

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(\alpha \beta )(10\u221210m)$ . | ω^{(αβ)}
. | $Tref(\alpha \beta )$ . |
---|---|---|---|---|

He | He | 2.33 | 0.66 | 273.15 |

He | Ar | 3.25 | 0.735 | 273.15 |

Ar | Ar | 4.17 | 0.81 | 273.15 |

Simulation case . | N_{p}
. | N_{c}
. | Δt (s)
. | n_{0} (1/m^{3})
. |
---|---|---|---|---|

Supersonic Couette flow | 76 000 | 100 | 10^{–6} | 3 · 10^{19} |

Diffusion flow | 20 000 | 100 | 2 · 10^{−6} | 6 · 10^{19} |

Simulation case . | N_{p}
. | N_{c}
. | Δt (s)
. | n_{0} (1/m^{3})
. |
---|---|---|---|---|

Supersonic Couette flow | 76 000 | 100 | 10^{–6} | 3 · 10^{19} |

Diffusion flow | 20 000 | 100 | 2 · 10^{−6} | 6 · 10^{19} |

All simulations are performed using the SPARTA code^{26} 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 *T*_{w} = 300 K. The higher wall moves in the *y*-direction with a velocity of *v*_{w} = 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 $\chi 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.

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.

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 *T*_{0} = 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 $\chi f(He)=\chi 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.

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.

Simulation case . | Time DSMC . | Time FP . |
---|---|---|

Supersonic Couette flow | 1 | 5.0 |

Diffusion flow | 1 | 8.3 |

Simulation case . | Time DSMC . | Time FP . |
---|---|---|

Supersonic Couette flow | 1 | 5.0 |

Diffusion flow | 1 | 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 mixtures^{10} 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 VSS^{27} or M1^{28} 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 AVAILABILITY

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