Hermite spectral method for multi-species Boltzmann equation
Abstract
We introduce a numerical scheme for the full multi-species Boltzmann equation based on Hermite spectral method. With the proper choice of expansion centers for different species, a practical algorithm is derived to evaluate the complicated multi-species binary collision operator. New collision models are built by combining the quadratic collision model and the simple BGK collision model under the framework of the Hermite spectral method, which enables us to balance the computational cost and accuracy.
Several numerical experiments are implemented to validate the dramatic efficiency of this new Hermite spectral method. Moreover, we can handle the problems with as many as 100 species, which is far beyond the capability of the state-of-art algorithms.
Keywords: multi-species Boltzmann equation;
multi-species binary collision operator; Hermite spectral method
1 Introduction
For rarefied gas flows where the molecular mean free path is comparable to the characteristic length, the continuum-fluid models such as Euler equations and Navier-Stokes equations are no longer accurate. The Boltzmann equation is an important model to study rarefied gas dynamics, which uses the distribution function to govern the dilute gas behavior at the mesoscopic scale. Gas mixtures with species diffusion, which play an essential role in the evaporation/condensation process, could also be described by the multi-species Boltzmann equation. For each species of the gases, it is represented by a distribution function , where is time, is the physical space, is the particle velocity space, and represents the kind of species. Thus, the multi-species Boltzmann equation, which describes a many-particle system comprised of a mixture of species, has the form below [7, 17]
| (1.1) |
where is the collision operator that models binary collisions between -th and -th species:
| (1.2) |
During the collisions, the momentum and energy are conserved. and denote the pre-collision and post-collision velocity pairs satisfying
| (1.3) |
where is a unit vector indicating the direction of the post-collision relative velocity, and are the mass of particles of -th and -th species respectively.
The collision kernel in (1.1) characterizes the interaction between the particles, which could be determined by the classical scattering theory based on the interaction potential as
| (1.4) |
with the deviation angle between and , and the impact factor [21]. For the multi-species Boltzmann equation, the most difficult part lies in the binary collision term due to its complex quadratic form and high-dimensional integral. Several collision models are introduced, such as the variable soft sphere model (VSS) [24] by assuming
| (1.5) |
where is the diameter and is the scattering parameter. For the diameter , we choose the same as in [21], which is borrowed from [2, Eq. (4.79)] as
| (1.6) |
with , the Gamma function, the viscosity index, and are the reference diameter and temperature. Particularly, we could obtain the variable hard sphere model (VHS) when and (especially, the Maxwell molecules and hard sphere model (HS) are corresponding to and respectively); and the VSS model when and . It is full of difficulty in numerically dealing with the binary collision term, especially for the VSS model, due to its extremely complicated form of integral.
For now, several numerical methods have been proposed to solve the multi-species Boltzmann equation. By stochastically modelling the binary interactions, the Direct Simulation Monte Carlo (DSMC) method [2] is widely used in the simulation of the Boltzmann equation. DSMC method is efficient in the highly-rarefied gas flows but does not work well in the unsteady-state problems due to its stochastical noise. In recent years, research on the deterministic methods has made great progress. Discrete velocity method (DVM) [13, 28] is a classical deterministic method. However, low efficiency of DVM has greatly limited its usage [32]. Spectral methods such as Fourier spectral method [33, 31, 10] have made great success in the simulation of the Boltzmann equation, and then, it has been extended to the multi-species Boltzmann equations [21, 36]. Another important spectral method is the Hermite spectral method [12, 35], where the Hermite polynomials are utilized as the basis functions to approximate the distribution function. An asymptotic-preserving scheme for the two-species binary collisional kinetic system was proposed in [11]. A comprehensive review of these methods can be referred to [8]. Another important method is the moment method, which was firstly proposed by Grad [14] in 1949. However, the non-hyperbolicity (even near the Maxwellian) of the moment equation led by Grad’s method, has blocked the its development. Recently, several regularization methods [3, 23] have been proposed to derive the globally hyperbolic moment systems.
In this paper, we are focusing on the Hermite spectral method. Since the Maxwellian is the steady-state solution of the Boltzmann equation, it is natural to consider expanding the distribution function with Hermite polynomials, whose weight function is just the Maxwellian. Another advantage of the Hermite spectral method is that several macroscopic variables such as the density, velocity, temperature, shear stress, and heat flux could be expressed by the expansion coefficients of the first few orders under the framework of the Hermite expansions. Recently, an algorithm to reduce the computational cost of approximating the binary collision term is proposed in [35], which is achieved by exactly calculating the expansion coefficients of the quadratic collision model and building new collision models through combining the quadratic collision term and the simplfied collision models such as the BGK model. This method is verified to make great success in the simulation of rarefied gas [20], and is soon extended to the field of the collisional plasma [27, 26].
Following similar lines, we develop the Hermite spectral method for the multi-species Boltzmann equation. Distribution functions of each species are expanded with the Hermite basis functions, with the expansion center in the basis function chosen differently for different species. We first simplify the complicated quadratic collision models, where the expansion centers are decided by the macroscopic velocity, the average temperature, and the mass of each species. Following the techniques in [35, 27], we will derive the algorithm to calculate the expansion coefficients of the binary collision models under the framework of the Hermite expansions. The different expansion centers for different species will make this calculation much easier. Moreover, the choice of the expansion centers makes the expansion coefficients of the collision models only depend on the mass ratio of the two particles, which also could be pre-computed and only need to be computed once. To reduce the computational cost of this complicated binary collision model, new collision models are built by combining the quadratic collision models and the simple multi-species BGK model [22], which enables us to obtain a high-order approximation with low computational cost and this advantage is also verified in the numerical experiments. The time-splitting scheme is utilized to deal with the convection term and the collision term separately, where the convection term is solved similarly as in [20, 4] with the standard finite volume method. With the new collision model in the collision step, the total computational cost to approximate one collision term is greatly reduced, which is about , where and are the length of the binary collision part and the order of Hermite expansion, respectively.
In the numerical experiments, the test for the exact Krook-Wu solution [25] validates the correctness of this new numerical method. The simulations are implemented for one-dimensional benchmark problems, including Couette flow and Fourier heat transfer, and the two-dimensional lid-driven cavity flow. To further verify the accuracy and efficiency of this Hermite spectral method, we finally implement a 100-species Krook-Wu solution with tolerable costs of time and memory space, where the number of species can be hardly handled by most of the methods nowadays.
The rest of this paper is organized as follows. In Sec. 2, the multi-species Boltzmann equation and the Hermite spectral method are introduced. The numerical method to calculate the expansion coefficients of the quadratic collision term in the framework of the Hermite spectral method is shown in Sec. 3, and the algorithm to build the new collision model is also presented in this section. The discussion of the moment equations and the description of the whole numerical scheme is given in Sec. 4. The numerical experiments are presented in Sec. 5, with some concluding remarks in Sec. 6. Some supplementary contents are presented in App. 7.
2 Preliminaries
For easier computation, we follow the similar nodimensionalization as in [21] to scale the variables as
| (2.1) |
where and are the characteristic length, number density, temperature and mass, with the character velocity defined as .
2.1 Multi-species Boltzmann equation
Substituting the nondimensionalization into the Boltzmann equation (1.1), we derive the dimensionless Boltzmann equation as (dropping for simplicity)
| (2.2) |
where the collision term is simplified as
| (2.3) |
and is the Knudsen number, which is defined as the ratio of the mean free path and characteristic length. For the VSS model, the dimensionless collision kernel has the form
| (2.4) |
From the nondimensionalization, the Knudsen number can be computed as
| (2.5) |
The macroscopic variables such as the number density , density , macroscopic velocity and temperature for each species could be derived from the distribution function as
| (2.6) | ||||
The stress tensor , heat flux and total energy could also be derived from the distribution function as
| (2.7) | ||||
The macroscopic variables for the whole system are defined as
| (2.8) |
Define , then the steady-state solution to the multi-species Boltzmann equation (2.2) has the form as [21]
| (2.9) |
with
| (2.10) |
which is named as the local Maxwellian. Here, we define the total average temperature as
| (2.11) |
which could express the total average temperature of the system. We also want to emphasize that due to the nonlinearity of the total energy respected to the macroscopic velocity , is different from in (2.8). Based on the conservation of mass, momentum, and energy during the collision, the properties hold for the collision term as
| (2.12) |
Due to the complex form of the quadratic collision term, several simplified collision models are introduced, such as the BGK collision model [22]
| (2.13) |
where is the collision frequency, and are decided by the macroscopic variables of species and , which should satisfy the conversation relationship (2.12). Here, we want to mention if it holds that
| (2.14) |
then , and . Once is decided, could be calculated with the conservation law (2.12), so as and . There are several types of BGK models to decide and , such as the work in [22, 15], which we will not discuss in detail.
For now, we have introduced the properties of the multi-species Boltzmann equation. Due to its high dimensionality and the complex collision terms, several numerical methods have been proposed, such as the stochastic DSMC method [2], the discrete velocity method [13, 32], and the Fourier spectral method [21, 36]. So far, the Hermite spectral method has been successfully utilized to solve Boltzmann equation with single species and the plasma with collisions [26, 20]. In this work, we will extend the methodology therein to solve the multi-species Boltzmann equation, which we will discuss in detail in the next sections.
2.2 Hermite spectral method
In the framework of the Hermite spectral method, the distribution function is approximated by a series of expansions. Following the method in [20], we choose the weighted Hermite polynomials as the basis function as
Definition (Hermite Polynomials).
For , the three-dimensional Hermite polynomial is defined as
| (2.15) |
with and . Here and are two parameters which are chosen problem-dependently. is the Maxwellian defined in (2.10). The Hermite polynomials have several useful properties when approximating the complex collision term, which are listed in Appendix 7.3.
With the Hermite polynomials, the distribution function is expanded as
| (2.16) |
where are the basis functions, is the problem-dependent expansion center. One of the advantages of this expansion is that the problem-dependent expansion center could also be interpreted as the ”moments” of the distribution function. With a properly chosen expansion center, the computational cost may be greatly reduced. More details will be revealed in the successive sections. By the orthogonality of the basis function (7.13), it holds that
| (2.17) |
The macroscopic variables (2.6) could also be expressed using as
| (2.18) |
Sometimes, the expansion center is chosen according to the numerical scheme, such as dealing with the force term when simulating the plasma with collision [26], or according to the average temperature in the whole space for the Fourier flow problem [20]. In this work, this expansion center is also chosen differently for the convection part and the collision term.
3 Modelling the quadratic collision term with Hermite spectral method
In this section, we will model the quadratic collision term using the Hermite spectral method. The influence of the convection term will be ignored temporarily, and we only discuss the algorithm to approximate the collision term in (2.2).
3.1 Hermite expansion of the collision term
For the homogeneous Boltzmann equation, the total momentum and the kinetic energy are conserved. With a proper frame of reference, we can suppose that the total macroscopic velocity in (2.8) is . With (2.11), the total average temperature is constant and solely decided by the total energy . Therefore, we choose the macroscopic velocity and total average temperature to decide the expansion center. For the th species, the expansion center is set as with and to approximate the collision term.
Thus, the basis function of is chosen as , and the variable is also omitted in the distribution function for simplicity. The Hermite expansion (2.16) is then written as
| (3.1) |
In what follows, we denote and as the Hermite polynomial and Maxwellian corresponding to . Substituting (3.1) in the the collision term (2.3), the collision could be expanded as
| (3.2) |
Now we consider how to obtain for fixed and . From the orthogonality (7.13), we can obtain that
| (3.3) |
with . Here, the second equality can be obtained by inserting (3.1) into (2.3), and
| (3.4) |
where is the post-collision velocity pair determined by (1.3). Here, and are functions related to the expansion center . For simplicity, from now on, the coefficients and will be shortened as and if the expansion center is set as ; and this parameter will be explicitly written out if expansion center other than is utilized.
In the computation of the expansion coefficients , it has a complicated high-dimensional integral. Moreover, for the multi-species particles, this should be calculated for each pair of , which has greatly increased the computational cost compared to the single species case. Here, the properties of the Hermite polynomials and the conservation of the collision model are all utilized to compute .
Denoting as the mass ratio, with the transitivity of Hermite polynomials (7.14) and the definition of Maxwellian (2.10), it holds that
| (3.5) |
where we omit the superscript of and short it as for simplicity. Besides, from the conversation of energy, it also holds that
| (3.6) |
Based on (2.10), (3.5) and (3.6), we can derive that
| (3.7) |
Then, (3.4) can be simplified as
| (3.8) |
With the properties of the basis functions, the calculation of could be greatly simplified, and the final result is listed in the theorem below.
Theorem 1.
The coefficients could be simplified as the form below:
| (3.9) |
where and , with . Besides, the symbol ‘’ means
The coefficients are given by
| (3.10) |
with the generalized combination number which is regarded as zero when or . The coefficients are defined as
| (3.11) |
where is the relative velocity and is the post-collision relative velocity.
The proof of Theorem 1 could be finished with the following steps. First, we introduce the lemma below
Lemma 1.
The proof of Lemma 1 could be found in App. 7.1. Based on Lemma 1, we could derive the corollary below
Corollary 1.
For any , define , then it holds that
| (3.13) |
The proof is straightforward by letting in (3.12). Based on Lemma 1 and Corollary 1, the proof of Theorem 1 is as below.
Proof of Theorem 1.
First, we denote , in (3.8), then it holds that
| (3.14) |
Similarly, with the conservation of the collision term, if , , then it holds that
| (3.15) |
Thus, following the transformation method in [35], we can derive that
| (3.16) |
Combining Lemma 1 and Corollary 1, we rewrite the above equality as an integral of and
| (3.17) |
where the coefficients are integrals with respect to defined in (3.11), and are integrals with respect to defined by
| (3.18) |
The second equality can be obtained by orthogonality (7.13). With (3.17) and (3.18), we complete the proof of Theorem 1. ∎
For some special collision kernels, the calculation of could be further simplified. For example, for the VSS collision kernel (1.5), where could be rewritten in the form
| (3.19) |
where is some constant and is a function of . Then for the VSS model, with the transitivity property (7.14), the coefficient in (3.11) could be further simplified as
| (3.20) |
Here, a further simplification of is introduced in [35], and we list the result in App. 7.2. For the VSS model, the ultimate form of only needs a two-dimensional integration, and the computational cost has been greatly reduced.
Remark 1.
From the calculation of , we want to emphasize that the expansion center is not unique. As long as the same expansion center and scaling factor is chosen for each species, (i.e. the expansion center is for species ), (3.4) could be computed in a similar way. It is easy to verify that for the VSS model, it holds that
| (3.21) |
and only depends on the mass ratio . Therefore, in the real simulation, we only need to compute and store with each different mass ratio in the system.
For now, we have introduced the algorithm to calculate the expansion coefficients for the Boltzmann collision term. Especially for the VSS model, the final computation is reduced into some elementary operations and a two-dimensional integral. Moreover, the coefficients could be precomputed offline and saved for the simulation. The computational cost of the coefficients is listed in Tab. 1, where the computational cost for the related coefficients is also listed.
| Coefficients | Formula | Computational cost |
|---|---|---|
| (7.9) | ||
| (7.12) | ||
| (7.11) | ||
| (3.20) and (7.8) | ||
| (3.10) | ||
| (3.8) |
Even though the expansion coefficients could be precomputed, the memory required to store is extremely large. Thus, in the framework of the Hermite spectral method to solve the Boltzmann equation, the direct simulation of the collision term is still quite expensive. Following the thought in [35], a new collision model will be built to approximate the quadratic collision term in the framework of Hermite expansion.
3.2 Approximation to the Boltzmann collision term
Supposing that the truncation order of the expansion approximation in (3.1) is , then the computational cost to compute the collision terms in (3.3) is , where is the number of species. Such computational cost is unacceptable for a large . For the spatially non-homogeneous problems, the computational cost will be much larger. Therefore, we want to build new collision models following the method proposed in [35, 20].
3.2.1 The BGK collision term
Before introducing the new collision model, we will discuss the BGK model (2.13) firstly. For the BGK model in the multi-species Boltzmann equation, other than that in the single-species case where the BGK model has only one form, there are several variants here. The essential part is how to choose , and collision frequency in (2.13), with different models discussed in related literature. It is assumed that is a linear combination of and in [22], and an estimation depending on the macroscopic variables is provided in [15]. In this work, for simplicity, we assume that and . Then, as is discussed in [16], it holds for and that
| (3.22) |
To approximate the BGK collision term under the framework of the Hermite expansion, we choose the same expansion center as in (3.1). Thus, the BGK collision term (2.13) could be expanded as
| (3.23) |
where
| (3.24) |
Here is the operator defined as below.
Definition.
First, define the function space, and as
| (3.25) |
Assume that for , then the projection operator is defined as
| (3.26) |
3.2.2 Building new collision model
In the problem that is far from equilibrium, the BGK model could not describe the movements of the particles well, so the quadratic collision model should be utilized. However, the quadratic collision is quite expensive to simulate. Thus, we want to build a new collision model combining the quadratic and BGK collision model.
For the single-species Boltzmann equation, the BGK collision model could not give the correct Prandtl number [34], and the Shakhov model has fixed this problem by modifying the Maxwellian term in the BGK model. To build the new collision model, we borrow this idea from the Shakhov model. Instead of modifying the Maxwellian, we will revise the expansion coefficients of the collision model under the framework of the Hermite spectral method.
In the new collision model, we first choose to decide the length of the expansion coefficients which are derived from the quadratic collision term, then the rest are obtained from the BGK collision model. Precisely, the new collision model could be written as
| (3.30) |
with
| (3.31) |
Here, are some positive constants deciding the damping rate, where we follow the same idea in [20] to decide this parameter. Since could be regarded as a matrix respect to and for each fixed , and then, we set as the minimum of all eigenvalues of among . Building the new collision model (3.30) here could also be understood as combining the quadratic and BGK collision models as in [20, 35]. In this way, the computational cost for the quadratic collision part in the new collision model will be , which will be reduced greatly compared to (3.2).
For the new collision model, the length of the quadratic collision model is problem-dependent as the further the system from the equilibrium, the larger is needed. However, there is no standard principle to decide . There is an adaptive method in [6] to choose , which is still an initial work, and we have not adopted the algorithm therein.
Remark 2.
In this section, we are focusing on approximating the quadratic collision term in the homogeneous problem. Therefore, the expansion center to calculate the expansion coefficients and build the new collision model (3.30) is set as . As is stated in Remark 1, the expansion center could be chosen arbitrarily, so as the expansion of the BGK collision term (3.23). The transition between different expansion centers could be derived through Theorem 2, and we refer the readers to [20] for more details.
4 Moment equations and numerical algorithm
In the last section, we have explained the approximation of the quadratic collision term for the homogeneous problem under the framework of the Hermite expansion. In this section, the numerical scheme to solve the whole Boltzmann equation (2.2) will be introduced. We will first introduce the moment systems in the framework of the Hermite expansions. Since the Strang splitting is adopted to solve the convection and collision part separately, the moment equations will also be derived for the convection and collision step. Precisely, the Boltzmann equation (2.2) is split into the convection step and the collision step as
-
•
convection step:
(4.1) -
•
collision step:
(4.2)
4.1 Moment equations
To obtain the finite moment system, we first conduct a truncation of the expansion (2.16) as
| (4.3) |
where the expansion center is chosen differently for the convection and collision step, and is the functional space
| (4.4) |
For the convection step, we follow the framework in [20], and the identical expansion center all over the spatial space is utilized for each species. Precisely, the choice of is based on a priori estimation of the problem such that the expansion center is not too far away from the average macroscopic variables. Then, substituting (4.3) into (4.1), we could derive the moment equations for the convection part as
| (4.5) |
where is a column vector of all for and is a constant matrix made up by the convection term, the detailed form of which can be referred to [20] and will be omitted here. Moreover, without loss of generality, we will record as if the corresponding distribution function . For the collision term, the expansion center could be the same as the convection term, the same expansion center as in Sec. 3, or some other choices could also be utilized. For any expansion center, substituting (4.3) into (4.2) and matching the basis functions on both sides, we can derive the moment equations for the collision step as
| (4.6) |
To solve the moment equations of the convection step (4.5) and the collision step (4.6), the standard finite volume method will be adopted, which we will introduce in detail in the following sections.
Here, we also want to emphasize that if different expansion centers are utilized for the convection and the collision steps, the transition of the distribution functions between different expansion centers can be achieved with Theorem 2.
4.2 Numerical scheme
Supposing the spatial domain is discretized by a uniform grid with the cell size and the cell centers . Using to approximate the average of over the th grid cell at time , then we will introduce the numerical scheme for the Boltzmann equation. It will split into two steps, which we will begin from the convection step. The system (4.5) can be solved by the forward-Euler method with time step length as follows:
| (4.7) |
where the finite volume method is employed for spatial discretization and is the result of convection step (4.5). is the HLL flux [18] given by
| (4.8) |
where and can be computed with the linear reconstruction (also used in [20])
| (4.9) |
or the 3-order WENO reconstruction [29]
| (4.10) |
where the square of in (4.10) indicates squaring by each expansion coefficient, and the superscript of is omitted in (4.9) and (4.10) for simplicity.
In (4.8), and are the minimum and maximum eigenvalue of respectively. Due to the same form of here as in [20], and are set as and . The time step is chosen to satisfy the CFL condition as
| (4.11) |
After obtaining the results at each cell in the convection step, we will update the collision step. Since the expansion center of the convection step and the collision step may be different, we will apply the algorithm in Theorem 2 to obtain the new expansion coefficients in the expansion center of the collision step.
For the collision step, we choose a similar expansion center as in Sec. 3. Precisely, the expansion center for species is , where is the the average macroscopic velocity in (2.8) and is the average temperature in (2.11) in the -th cell at time step . Since the collision term does not change the average macroscopic velocity and the average temperature in each cell, the expansion center here is computed using the distribution function after the convection step. Thus, with the projection operator defined in (3.26), we first compute as
| (4.12) |
Then, at the -th cell, we compute the collision term in (3.31) with the distribution function . The forward-Euler scheme is adopted here to update the distribution function in the collision step (4.6) as
| (4.13) |
The final step is to project the distribution function to the expansion center of the convection step as
| (4.14) |
For the forward Euler scheme in (4.14), the Runge-Kutta scheme could also be utilized, and the time step may also be restricted by the collision term. A small enough and problem-dependent should be adopted, which we will not discuss in detail here. The Maxwell boundary condition [30] is adopted here and the detailed implementation can be referred to [20].
For now, we have finished the whole numerical scheme, and it will be summarized as Algorithm 1.
5 Numerical experiments
In this section, several numerical examples are studied to validate the numerical algorithm proposed in the last section. Precisely, the spatial homogeneous test with Krook-Wu solution, two spatially one-dimensional and a spatially two-dimensional problem are tested. For the spatially non-homogeneous problems, the SPARTA [9] is utilized to obtain the reference solution with DSMC method, and the mixture of Argon and Krypton is taken as the working gas, the parameters of which are listed in Tab. 9 in App. 7.4. In the final example, we study the Krook-Wu solution with species to show the superiority of this numerical method.
5.1 Spatially homogeneous case: Krook-Wu solution
For the spatially homogeneous problem, we will consider a constant collision kernel. In this case, the spatially homogeneous two-species Boltzmann equation can be simplified as
| (5.1) |
where is the number of species. and are some positive constants. This problem is also studied in [21, 35]. There exists an exact solution to this problem [25] as
| (5.2) |
where
| (5.3) |
with
| (5.4) |
and is a constant. The number density should be chosen such that , while the mass can be arbitrary with the following condition satisfied
| (5.5) |
Here, we choose the same parameters as in [35] that , , and the constant . Moreover, the constants are chosen as , and .
It can be directly computed from (5.2) that the average temperature is , and therefore we choose the expansion center as Based on this expansion center, we can obtain the expansion coefficient of the exact solution (5.2) as
| (5.6) |
for all and .
In the computation, the time step length is fixed as , and the fourth-order Runge-Kutta scheme is adopted. The total expansion order is , and two cases with different length of the quadratic collision term are tested, where is set as and respectively. The numerical results of from to are given in Fig. 1, since the expansion coefficients for the lower orders are all constant with respect to . We can find that for the two cases, the numerical results of both species all match well with the exact solution at any time.
To see the numerical results of the distribution function more clearly, we randomly select several vectors for each species and show the results of along at time and , which are shown in Fig. 2. In Fig. 2, three cases are plotted with respectively. For the case , there is an obvious distance between the numerical solution and the exact solution, and for the case , the numerical solution is still far away from the numerical solution, while the results for the case are almost on top of the exact solution. From this, we can see that with the increase of , the numerical solution could describe the distribution function better. The results in Fig. 2 also show that the -order approximation is much better than those with -order, which exhibits the importance of the expansion coefficients with the order higher than and validates the advantage of this numerical algorithm as well.
5.2 1D case: Couette flow
In this section, we will consider the Couette flow, which is also tested in [21, 20]. The Ar-Kr mixture between two infinite parallel plates with a distance of m will be studied when they arrive at the steady state. Both the left and right plates have the temperature K and move in the opposite direction along with the plate with the speed m/s. Besides, both walls are purely diffusive. The VSS collision model (1.5) is utilized here, the detailed parameters of which are shown in Tab. 9. In this simulation, a uniform grid with cells is utilized for the spatial discretization with the WENO reconstruction adopted. The CFL condition is set as . Three tests with different number densities are carried out, which are corresponding to different Knudsen numbers. The detailed parameters for the initial conditions are listed in Tab. 2, where the parameters for the nondimensionalization are also listed. With (2.5), we can obtain the corresponding Knudsen number for the three cases as in Tab. 2.
| Case 1 | Case 2 | Case 3 | |
| Initial conditions: | |||
| Temperature (K) | 273 | 273 | |
| Velocity (m/s) | (0, 0, 0) | (0, 0, 0) | (0, 0, 0) |
| Number density (m-3) | 1.68e21 | 5.6e20 | 1.68e20 |
| Number density (m-3) | 8.009e20 | 2.667e20 | 8.009e19 |
| Characteristic variables: | |||
| Length (m) | 1e-3 | 1e-3 | 1e-3 |
| Mass (kg) | 6.63e-27 | 6.63e-27 | 6.63e-27 |
| Number density (m-3) | 1.68e21 | 5.6e20 | 1.68e20 |
| Temperature (K) | 273 | 273 | 273 |
| Velocity (m/s) | 238.377 | 238.377 | 238.377 |
| Knudsen number , | (0.793, 0.804) | (2.379, 2.411) | (7.931, 8.036) |
| Knudsen number , | (0.606, 0.555) | (1.819, 1.664) | (6.065, 5.548) |
For the Couette flow, the average macroscopic velocity of the whole domain is zero. Therefore, we choose the expansion center for the convection step as . The average temperature for the initial condition after non-dimensionalization is . Therefore, in the convection step, the expansion center of temperature is set as and , where and are the mass after non-dimensionalization. In the collision step, the expansion center will be decided by local macroscopic variables. The details can be referred to Algorithm 1.
In the simulation, the macroscopic velocity in the -direction , temperature , stress tensor , and the heat flux are studied. In the first case, since the Knudsen number is relatively small, we set the length of the quadratic collision term as , and the total moment number as . The numerical results are shown in Fig. 3, from which we can see that for these both species Ar and Kr, the numerical solutions of these four variables match well with the reference solutions obtained by the DSMC method. Though there is a little difference in the temperature between the numerical solution and the reference, the relative error is less than . For the second case, we still set and . The behavior of the numerical solution is similar to that of the first case, as is shown in Fig. 4, where the numerical solution and the reference solution are almost on top of each other. For the third case, since the Knudsen number is relatively large, we increase to and keep due to the limitation of the computational cost. From Fig. 5, we can see that for the two species, though there is a little disparity between the numerical solution and the reference solution for the velocity and the heat flux on the boundary, the relative error is quite small. We also find that in this case, there are small oscillations for the reference solution in temperature , while the numerical results are still quite smooth.
| Case 1 | Case 2 | Case 3 | |
| Run-time data: | |||
| Total CPUtime (s): | 15264 | 17217 | 27273 |
| Elapsed time(Wall time) (s): | 4923.69 | 6600.08 | 10164.6 |
The simulation of the Couette flow is done on the CPU model Intel Xeon E5-2697A V4 @ 2.6GHz, and 4 threads are used in this test. We implement the simulation for a long enough time to obtain the steady-state solution. In the real simulation, the final time for the simulation is for all three cases. The total CPU time and wall time are recorded in Tab. 3, from which we can see that the computational time is increasing with the increase of . The computational time also shows the high efficiency of this Hermite spectral method, even for the quite rarefied cases.
5.3 1D case: Fourier heat transfer
The second example of the 1D case is the Fourier heat transfer problem, which is also widely studied [21, 20, 19]. In this example, we also consider the Ar-Kr mixture between two infinitely large plates. Similar to the Couette flow, the distance between the walls is m, and both boundaries are assumed to be purely diffusive. Different from the Couette flow, the Fourier heat transfer problem considers a distinction of temperature instead of velocity on two walls. In this example, a different collision model, the VHS model, is adopted for the collision term, the detailed parameters of which are listed in Tab. 9. For the Fourier heat transfer problem, the steady-state solution is studied and the DSMC method is utilized to obtain the reference solution.
In the simulation, the same WENO reconstruction with uniform grid cells and is utilized as the Couette flow problem. Two different number densities and two different boundary conditions are tested. The detailed parameters for the initial and boundary conditions are listed in Tab. 4, where the parameters for the nondimensionalization are also listed. The same expansion centers and are adopted here for the convection step while the expansion centers of the collision step are decided with local macroscopic variables as in Algorithm 1.
| Case 1 | Case 2 | Case 3 | Case 4 | |
| Initial conditions | ||||
| Temperature (K) | 273 | 273 | 273 | 273 |
| Velocity (m/s) | (0, 0, 0) | (0, 0, 0) | (0 ,0, 0) | (0, 0, 0) |
| Number density (m-3) | 1.68e21 | 5.6e20 | 1.68e21 | 5.6e20 |
| Number density (m-3) | 8.009e20 | 2.667e20 | 8.009e20 | 2.667e20 |
| Characteristic variables | ||||
| Length (m) | 1e-3 | 1e-3 | 1e-3 | 1e-3 |
| Mass (kg) | 6.63e-27 | 6.63e-27 | 6.63e-27 | 6.63e-27 |
| Number density (m-3) | 1.68e21 | 5.6e20 | 1.68e21 | 5.6e20 |
| Temperature (K) | 273 | 273 | 273 | 273 |
| Velocity (m/s) | 238.377 | 238.377 | 238.377 | 238.377 |
| Knudsen number , | (0.77, 0.782) | (2.311, 2.346) | (0.77, 0.782) | (2.311, 2.346) |
| Knudsen number , | (0.54, 0.591) | (1.62, 1.774) | (0.54, 0.591) | (1.62, 1.774) |
| Boundary conditions | ||||
| Left wall temperature (K) | 223 | 223 | 109.2 | 109.2 |
| Right wall temperature (K) | 323 | 323 | 436.8 | 436.8 |
For the Fourier heat transfer problem, the density , temperature , stress tensor , and heat flux are studied. In the first and second cases in Tab. 4, since that difference between the two boundary conditions is small, which are , we set the order of the total moment number and the length of the quadratic collision term as . The numerical results are shown in Fig. 6 and 7. In Fig. 6, the numerical solutions match well with the reference solutions, and there are some oscillations in the reference solutions in the stress tensor, while those of the numerical solutions are still smooth. In Fig. 7, for the temperature , stress tensor , and heat flux all match well with the reference solution, although there is a little discrepancy for the density on the left boundary. However, the relative error is still quite small, which is about .
In the third and fourth cases in Tab. 4, the boundary condition is chosen as , where the ratio of the temperature is . In case 3, since the Knudsen number is still small, we set , the results of which are shown in Fig. 8. Most of the numerical solutions and the reference solutions are on top of each other, excepting the little difference in the stress tensor. In case 4, we increase to due to the large Knudsen number and keep due to the limitation of the computational cost. Fig. 9 shows that even for this large Knudsen number and the ratio of temperature on the boundary condition, we can still catch the behavior of the two species.
The Fourier heat transfer problem is simulated using the same machine as the Couette flow problem. Here, the final computation time is still , and the total CPU time and wall time are shown in Tab. 5. We can find that with the same , the total run times for the first three cases are almost the same, while that of the fourth case is much longer, which may be due to the large Knudsen number and the high ratio of the temperature on the boundary conditions.
| Case 1 | Case 2 | Case 3 | Case 4 | |
| Run-time data: | ||||
| Total CPUtime (s): | 16384 | 17599 | 16978 | 43887 |
| Elapsed time(Wall time) (s): | 5238.95 | 5839.7 | 5366.48 | 18839.9 |
5.4 2D case: Lid-driven cavity flow
In this section, the two-dimensional lid-driven cavity flow is studied, which is also tested in [20, 5, 28]. In this test, the Ar-Kr mixture is also taken as the working gas, and the HS collision kernel is utilized here in the collision model, the detailed parameters of which are listed in Tab. 9. For the lid-driven cavity flow, the gas is confined in a square cavity with side length . All walls are purely diffusive and have temperature . The top lid moves to the right with a constant speed , while the other three boundaries are stationary.
In this simulation, a uniform grid with cells is utilized for the spatial discretization with the linear reconstruction adopted. The CFL condition is set as . Two tests with different number densities are carried out, which are corresponding to different Knudsen numbers. The detailed parameters for the initial conditions are listed in Tab. 6, where the parameters for the non-dimensionalization are also listed. With (2.5), we can obtain the corresponding Knudsen number for the two cases as in Tab. 6. The same expansion centers and as the Couette flow and Fourier heat transfer problem are adopted here for the convection step. For the collision step, the expansion centers are also decided with local macroscopic variables as in Algorithm 1.
| Case 1 | Case 2 | |
| Initial conditions | ||
| Temperature (K) | 273 | 273 |
| Velocity (m/s) | (0, 0, 0) | (0, 0, 0) |
| Number density (m-3) | 1.708e22 | 1.708e21 |
| Number density (m-3) | 8.141e21 | 8.141e20 |
| Characteristic variables | ||
| Mass (kg) | 6.63e-27 | 6.63e-27 |
| Number density (m-3) | 1.68e21 | 5.6e20 |
| Temperature (K) | 273 | 273 |
| Velocity (m/s) | 238.377 | 238.377 |
| Knudsen number (, ) | (0.1, 0.101) | (1, 1.014) |
| Knudsen number (, ) | (0.07, 0.076) | (0.7, 0.762) |
In the simulation, the temperature and the stress tensor of the two species are studied. For case 1 in Tab. 6, since the Knudsen number is small, we set the total expansion order and the length of the quadratic collision term as . The numerical results are plotted in Fig. 10, where we find that the numerical solution matches well with the reference solution obtained by the DSMC method. Moreover, small oscillations exist in the reference solution of the temperature, while the numerical solution keeps smooth. For case 2 in Tab. 6, we set as , and the numerical results are shown in Fig. 11. We find that for the case with increased Knudsen number, the numerical results and the reference solution are also almost the same, except that the reference solution has some oscillations and the numerical solutions are all smooth.
5.5 100-species KW solution
To verify the superiority of the numerical methods proposed in Algorithm 1 in the cost of time and space, we again implement the KW solution test in 5.1. The difference is that we will consider the 100-species case, which is much more challenging in the required memory and running time. The numerical simulation with such a quantity of species is unreachable by most methods with acceptable consumption of time and space. The general -species Boltzmann equation with constant collision kernel has the same form as in (5.1), and it also has the exact solution [25], the form of which is similar to (5.2). Precisely, the exact solution has the form below
| (5.7) |
where
| (5.8) |
Here is a constant which is large enough to ensure is positive for each , and we set . The same expansion center is utilized here for each species as in Sec. 5.1. The corresponding coefficients can be computed as
| (5.9) |
In the simulation, the fourth-order Runge-Kutta scheme with the time step length is utilized. The total expansion order and the length of the quadratic collision is set as . To show the efficiency of this new method, instead of showing the numerical solutions, we evaluate the relative error and weighted error
| (5.10) |
for each species. The integral will be evaluated with the 10-point Gauss-Hermite quadrature in each dimension. Two times and are recorded, the detailed of which is shown in Tab. 7. We can see that the two errors are all quite small, which are almost at the order . The results clearly show the accuracy of this method. Moreover, with time evolution, the distribution function is approximating the Maxwellian, and that is why the error reduces with time increasing.
| Average of 100 species: | ||
|---|---|---|
| Maximum: | ||
| Average of 100 species: | ||
| Maximum: |
The computational cost including the total computational time and the total memory used is shown in Tab. 8. The simulation is run on the CPU model Intel Xeon E5-2697A V4 @ 2.6GHz with 8 threads. Tab. 8 shows that the computing time for each collision term is only s, which is quite short for such a problem with species. Moreover, the total memory to store the expansion coefficients , which are represented in the double-precision floating-point format, is GB, which is also quite easy to accomplish for computers nowadays. This is due to the special algorithm to calculate , which is also the unique point in this new method. For Maxwell molecules (i.e., constant collision kernel ), it has been proved in [35] that can be nonzero only when . Combining (3.9) and (3.20), we know that is nonzero only when . Therefore, we only need to store the coefficients of . The memory for each can be reduced to . Besides, since the set of expansion coefficients is only decided by the mass ratio , we only need to store one group of coefficients for duplicate , which will further reduce the storage cost. Furthermore, the new collision model combined by the quadratic collision term and the BGK collision term will reduce the computational cost and greatly improve the efficiency of this new method. All of these special properties make the simulation of this 100-species KW solution problem possible, while similar simulations are never seen before.
| Total CPU time(s): | 2102 | 10642 |
|---|---|---|
| Collision terms computed: | ||
| CPU time per collision term(s): | ||
| Elapsed time(Wall time)(s): | 302.02 | 1759.32 |
| Total memory of (Gigabytes): | 18.85 | 18.85 |
6 Conclusion
In this paper, we have developed a numerical scheme of multi-species Boltzmann equation based on the Hermite spectral method. A new collision model is built by combing the quadratic collision term and the BGK collision model to balance the computational cost and the accuracy. The expansion coefficients of the quadratic collision model are almost computed with the properties of the Hermite basis functions. Several numerical examples are implemented to validate this new numerical method. Even for problems with 100-species, this new method could capture the behavior of the particles well with a low computational cost.
Such numerical results make this new numerical method promising when applied to high-dimensional and more complicated problems. Moreover, the GPU may be adopted to further speed up and the research on the quantum Boltzmann equations are ongoing.
Acknowledgements
We thank Prof. Zhenning Cai from NUS for his valuable suggestions. The work of Yanli Wang is partially supported by the National Natural Science Foundation of China (Grant No. 12171026, U1930402 and 12031013).
7 Appendix
7.1 Proof of Lemma 1
Now we provide the proof of Lemma 1. For convenience, we will omit the superscript of Hermite polynomials in the proof.
proof of Lemma 1.
First, it is easy to verify that , = and
| (7.1) |
By the Leibniz rule, we have the following relation:
| (7.2) |
where . Denote
| (7.3) |
Based on the orthogonality of Hermite polynomials, we only need to prove that
| (7.4) |
With , we can simplify (7.3) as
| (7.5) |
where the second equality uses the definition (2.15) and the third equality is obtained using integration by parts. From the differentiation relation (7.16)
| (7.6) |
and the orthogonality (7.13), it holds that (7.5) is nonzero only when , , , which implies . When holds, we can apply (7.6) to (7.5) and get
| (7.7) |
Thus (7.4) holds, which completes the proof of the lemma. ∎
7.2 The algorithm of
In this section, we will briefly introduce the algorithm to obtain , and readers can be referred to [35, Theorem 2] for the detailed algorithm. As is shown in [35], it holds that
| (7.8) |
where
| (7.9) |
and is the coefficient of in the polynomial
| (7.10) |
and
| (7.11) |
where are the Laguerre polynomials and are the Legendre polynomials. Moreover, is in fact decided by based on (1.4). As for , we can obtain with the recursion formula:
| (7.12) |
then we can derive the expression of based on (7.12).
7.3 Properties of Hermite polynomials
For the Hermite polynomials (2.15), several important properties are listed below
Property 1.
(Orthogonality)
| (7.13) |
Property 2.
(Transitivity)
| (7.14) |
Property 3.
(Recursion)
| (7.15) |
Property 4.
(Differential of Hermite polynomial)
| (7.16) |
7.4 Solver configurations for spatially inhomogeneous cases
In the numerical simulation, SPARTA [9] is employed to carry out the DSMC simulation as the reference. In spatially inhomogeneous cases, the mixture of Argon and Krypton will be taken as the working gas. Three collision models, the VSS, VHS, and HS collision kernels, are utilized in the simulation. The parameters of these three collision models between Argon and Krypton are shown in Tab. 9, where the reference diameter is derived with [2, Eq. (4.62)].
| Collision model for Ar-Kr mixture | VSS | VHS | HS |
| Molecular mass: (kg) | 6.63 | 6.63 | 6.63 |
| Molecular mass: (kg) | 13.91 | 13.91 | 13.91 |
| Ref. viscosity: (Pa s) | 2.117 | 2.117 | 2.117 |
| Ref. viscosity: (Pa s) | 2.328 | 2.328 | 2.328 |
| Viscosity index: () | (0.81, 0.805) | (0.81, 0.805) | (0.5, 0.5) |
| Viscosity index: () | (0.805, 0.8) | (0.805, 0.8) | (0.5, 0.5) |
| Scattering parameter: () | (1.4, 1.36) | (1, 1) | (1, 1) |
| Scattering parameter: () | (1.36, 1.32) | (1, 1) | (1, 1) |
| Ref. diameter: ()() | (4.11, 4.405) | (4.17, 4.465) | (3.63, 3.895) |
| Ref. diameter: ()() | (4.405, 4.7) | (4.465, 4.76) | (3.895, 4.16) |
| Ref. temperature: ()(K) | 273 | 273 | 273 |
| Ref. temperature: ()(K) | 273 | 273 | 273 |
References
- [1] M. Abramowitz and I. Stegun. Handbook of mathematical functions with formulas, graphs, and mathematical tables. New York: Dover, 1964.
- [2] G. Bird. Molecular Gas Dynamics and the Direct Simulation of Gas Flows. Oxford: Clarendon Press, 1994.
- [3] Z. Cai, Y. Fan, and R. Li. Globally hyperbolic regularization of Grad’s moment system in one dimensional space. Comm. Math. Sci., 11(2):547–571, 2013.
- [4] Z. Cai, Y. Fan, and R. Li. A framework on moment model reduction for kinetic equation. SIAM J. Appl. Math., 75(5):2001–2023, 2015.
- [5] Z. Cai and M. Torrilhon. Numerical simulation of microflows using moment methods with linearized collision operator. J. Sci. Comput., 74:336–374, 2018.
- [6] Z. Cai and Y. Wang. Numerical solver for the Boltzmann equation with self-adaptive collision operators. arXiv:2102.08559, 2021.
- [7] C. Cercignani. The Boltzmann Equation and Its Applications. Springer-Verlag, New York, 1988.
- [8] G. Dimarco and L. Pareschi. Numerical methods for kinetic equations. Acta Numerica, 23:369–520, 2014.
- [9] M. Gallis, J. Torczynski, S. Plimpton, D. Rader, T. Koehler, and J. Fan. Direct simulation monte carlo: the quest for speed. AIP Conference Proceedings, 1628(1):27–36, 2014.
- [10] I. Gamba, J. Haack, C. Hauck, and J. Hu. A fast spectral method for the Boltzmann collision operator with general collision kernels. SIAM J. Sci. Comput., 39(14):B658–B674, 2017.
- [11] I. Gamba, S. Jin, and L. Liu. Asymptotic-preserving schemes for two-species binary collisional kinetic system with disparate masses i: time discretization and asymptotic analysis. Comm. Math. Sci., 17:1257–1289, 2019.
- [12] I. Gamba and S. Rjasanow. Galerkin-Petrov approach for the Boltzmann equation. J. Comput. Phys., 366:341–365, 2018.
- [13] D. Goldstein, B. Sturtevant, and J. E. Broadwell. Investigations of the motion of discrete-velocity gases. Progress in Astronautics and Aeronautics, 117:100–117, 1989.
- [14] H. Grad. On the kinetic theory of rarefied gases. Comm. Pure Appl. Math., 2(4):331–407, 1949.
- [15] J. Haack, C. Hauck, C. Klingenberg, M. Pirner, and S. Warnecke. A consistent BGK model with velocity-dependent collision frequency for gas mixtures. J. Stat. Phys., 184(31), 2021.
- [16] J. Haack, C. Hauck, and M. Murillo. A conservative entropic multispecies BGK model. J. Stat. Phys., 168:826–856, 2017.
- [17] S. Harris. An introduction to the theory of the Boltzmann equation. Holt, Rinehart and Winston, Inc., 1971.
- [18] A. Harten, P. Lax, and B. Leer. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Rev., 25(1):35–61, 1983.
- [19] Z. Hu and Z. Cai. Burnett spectral method for high-speed rarefied gas flows. SIAM J. Sci. Comput., 42(5):B1193–B1226, 2020.
- [20] Z. Hu, Z. Cai, and Y. Wang. Numerical simulation of microflows using Hermite spectral methods. SIAM J. Sci. Comput., 42(1):B105–B134, 2020.
- [21] S. Jaiswal, A. Alexeenko, and J. Hu. A discontinuous Galerkin fast spectral method for the multi-species Boltzmann equation. Comput. Methods Appl. Mech. Engrg., 352:56–84, 2019.
- [22] C. Klingenberg, M. Pirner, and G. Puppo. A consistent kinetic model for a two-component mixture with an application to plasma. Kinet. Relat. Models, 10(2):445–465, 2017.
- [23] J. Koellermeier and M. Torrilhon. Numerical study of partially conservative moment equations in kinetic theory. Commun. Comput. Phys., 21(4):981–1011, 2017.
- [24] K. Koura and H. Matsumoto. Variable soft sphere molecular model for inverse-power-law or Lennard-Jones potential. Phys. Fluids A, 3(10), 1991.
- [25] M. Krook and T. T. Wu. Exact solutions of the Boltzmann equation. Phys. Fluids, 20(10):1589–1595, 1977.
- [26] R. Li, Y. Ren, and Y. Wang. Hermite spectral method for Fokker-Planck-Landau equation modeling collisional plasma. J. Comput. Phys., 434:110235, 2021.
- [27] R. Li, Y. Wang, and Y. Wang. Approximation to singular quadratic collision model in Fokker-Planck-Landau equation. SIAM J. Sci. Comput., 42(3):B792–B815, 2020.
- [28] C. Liu and K. Xu. A unified gas-kinetic scheme for micro flow simulation based on linearized kinetic equation. Adv. Aerodyn., 2(21), 2020.
- [29] X. Liu, S. Osher, and T. Chan. Weighted essentially non-oscillatory schemes. J. Comput. Phys., 115:200–212, 1994.
- [30] J. Maxwell. On stresses in rarefied gases arising from inequalities of temperature. Proc. R. Soc. Lond., 27:304–308, 1878.
- [31] C. Mouhot and L. Pareschi. Fast algorithms for computing the Boltzmann collision operator. Math. Comp., 75(256):1833–1852, 2006.
- [32] A. Panferov and A. Heintz. A new consistent discrete-velocity model for the Boltzmann equation. Math. Method Appl. Sci., 25(7):571–593, 2002.
- [33] L. Pareschi and B. Perthame. A Fourier spectral method for homogeneous Boltzmann equa- tions. Transport Theor. Stat., 25:369–382, 1996.
- [34] H. Struchtrup. Macroscopic Transport Equations for Rarefied Gas Flows: Approximation Methods in Kinetic Theory. Springer, 2005.
- [35] Y. Wang and Z. Cai. Approximation of the Boltzmann collision operator based on hermite spectral method. J. Comput. Phys., 397:108815, 2019.
- [36] L. Wu, J. Zhang, J. Reese, and Y. Zhang. A fast spectral method for the Boltzmann equation for monatomic gas mixtures. J. Comput. Phys., 298:602–621, 2015.