# Approximation of the Boltzmann Collision Operator Based on Hermite Spectral Method

###### Abstract

Based on the Hermite expansion of the distribution function, we introduce a Galerkin spectral method for the spatially homogeneous Boltzmann equation with the realistic inverse-power-law models. A practical algorithm is proposed to evaluate the coefficients in the spectral method with high accuracy, and these coefficients are also used to construct new computationally affordable collision models. Numerical experiments show that our method captures the low-order moments very efficiently.

Keywords: Boltzmann equation, Hermite spectral method, inverse power law

## 1 Introduction

Over a century ago, Boltzmann devised a profound equation describing the statistical behavior of gas molecules. A number of interesting theoretical and practical problems emerged due to the birth of this equation, among which the numerical simulation for this six-dimensional Boltzmann equation is one significant topic after the invention of computers. The difficulty comes partly from its high dimensionality, and partly from its complicated integral operator modeling the binary collision of gas molecules. People have been using the Monte Carlo method [3] to overcome the difficulty caused by high dimensionality, but nowadays, a six-dimensional simulation using a deterministic solver is no longer unaffordable due to the fast growth of computer flops. Fully six-dimensional computations are carried out in [19, 10] for the simplified BGK-type collision terms. However, numerical simulation of the original Boltzmann equation with the binary collision operator still requires a large amount of computational resources [11].

Currently, the deterministic discretization of the binary collision operator can be categorized into three types: the discrete velocity method [15, 25], the Fourier spectral method [26, 12], and the Hermite spectral method [17, 13]. The discrete velocity method is hardly used in the numerical simulation due to its low order of convergence [25], whereas the Fourier spectral method is more popular because of its fast convergence rate and high numerical efficiency. Compared with the Fourier spectral method which requires periodization of the distribution function, the Hermite spectral method looks more natural since the basis functions with orthogonality in are employed. In fact, the Hermite spectral method has a longer history and has been known as the moment method since Grad’s work [16]. Grad proposed in [16] a general method to find the expansion of the binary collision term with Hermite basis functions. Later, a similar way to expand the binary collision term using Sonine polynomials (also known as spherical Hermite polynomials) was proposed in [24]. The techniques used in formulating the expansion are also introduced in the book [30].

Despite these works, the Hermite spectral method is used in the numerical simulation only until recently [17, 13]. There are two major difficulties in applying this method: one is the evaluation of the coefficients in the expansion of the collision operator; the other is the huge computational cost due to its quadratic form. Although the general procedure to obtain the coefficients is given in [16, 30], following such a procedure involves expansion of a large number of huge polynomials, which is quite expensive even for a modern computer algebraic system; Kumar [24] provided a formula in his expansion using Sonine polynomials, while the formula involves evaluation of a large number of Talmi coefficients, which is not tractable either. As for the computational cost, the computational time of one evaluation of the collision operator is proportional to the cube of the number of degrees of freedom, while in the Fourier spectral method, the time complexity for a direct Galerkin discretization is only the square of the number of modes.

This work is devoted to both of the aforementioned issues. On one hand, by using a number of properties for relavant polynomials, we provide explicit formulas for all the coefficients appearing in the expansion of the collision operator with the Hermite spectral method. These formulas are immediately applicable in the sense of coding, and the computational cost is affordable for a moderate number of degrees of freedom. On the other hand, we combine the modeling strategy and the numerical technique to form a new way to discretize the collision term, where only a portion in the truncated series expansion is treated “quadratically”, and the remaining part just decays exponentially as in the BGK model. Thus the computational cost is greatly reduced and we can still capture the evolution of low-order moments accurately.

The rest of this paper is organized as follows. In Section 2, we briefly review the Boltzmann equation and the Hermite expansion of the distribution function. In Section 3, we first give an explicit expression for the series expansion of the quadratic collision operator, and then construct approximate collision models based on such an expansion. Some numerical experiments verifying our method are carried out in Section 4 and some concluding remarks are made in Section 5. Detailed derivation of the expansion is given in the Appendix.

## 2 Boltzmann equation and Hermite expansion of the distribution function

This section is devoted to the introduction of existing works needed by our further derivation. We will first give a brief review of the Boltzmann equation and the IPL (Inverse-Power-Law) model, and then introduce the expansion of the distribution function used in the Hermite spectral method.

### 2.1 Boltzmann equation

The Boltzmann equation describes the fluid state using a distribution function , where is the time, is the spatial coordinates, and stands for the velocity of gas molecules. The governing equation of is

(2.1) |

where is the collision operator which has a quadratic form

(2.2) |

where and is a unit vector. Hence is a one-dimensional integration over the unit circle perpendicular to . The post-collisional velocities and are

(2.3) | ||||

and from the conservation of momentum and energy, it holds that

(2.4) |

The collision kernel is a non-negative function determined by the force between gas molecules.

In this paper, we are mainly concerned with the IPL model, for which the force between two molecules is always repulsive and proportional to a negative power of their distance. In this case, the kernel in (2.2) has the form

(2.5) |

where is the index in the power of distance. The case corresponds to the “hard potential”, and the case corresponds to the “soft potential”. When , the collision kernel is independent of , and in this model the gas molecules are called “Maxwell molecules”. The dimensionless impact parameter is related to the angle by

(2.6) |

and is a positive real number satisfying

(2.7) |

It can be easily shown that the above equation of admits a unique positive solution when and .

Apparently, the quadratic collision term is the most complicated part in the Boltzmann equation. In this paper, we will focus on the numerical approximation of . For simplicity, we assume that the gas is homogeneous in space, and thus we can remove the variable in the distribution function to get the spatially homogeneous Boltzmann equation

(2.8) |

It is well known that the steady state of this equation takes the form of the Maxwellian:

(2.9) |

where the density , velocity and temperature can be obtained by

(2.10) |

These quantities are invariant during the evolution, and therefore (2.10) holds for any . By selecting proper frame of reference and applying appropriate non-dimensionalization, we can obtain

(2.11) |

and thus the Maxwellian (2.9) is reduced to

(2.12) |

Hereafter, the normalization (2.11) will always be assumed.

In the literature, people have been trying to avoid the complicated form of the collision operator by introducing simpler approximations to it. For example, the BGK collision model

(2.13) |

was proposed in [2]. Here is the mean relaxation time, which is usually obtained from the first approximation of the Chapman-Enskog theory [8]. When (2.13) is used to approximate the IPL model,

(2.14) |

where . With replaced by in (2.8), the collision process becomes an exponential convergence to the Maxwellian. Such a simple approximation provides incorrect Prandtl number . Hence some other models such as the Shakhov model [28] and ES-BGK model [18] are later proposed to fix the Prandtl number by changing the Maxwellian in (2.13) to a non-equilibrium distribution function. We will call these models “BGK-type models” hereafter.

Numerical evaluations on these BGK-type models can be found in [14, 9], where one can find that these approximations are not accurate enough when the non-equilibrium is strong. Hence the study on efficient numerical methods for the original Boltzmann equation with the quadratic collision operator is still necessary.

### 2.2 Series expansion of the distribution function

Our numerical discretization will be based on the following series expansion of the distribution function in the weighted space :

(2.15) |

where is the Maxwellian, and we have used the abbreviation

(2.16) |

In (2.15), are the Hermite polynomials defined as follows:

###### Definition 1 (Hermite polynomials).

The expansion (2.2) was proposed by Grad in [16], where such an expansion was used to derive moment methods. The relation between the coefficients and the moments can be seen from the orthogonality of Hermite polynomials

(2.18) |

For example, by the above orthogonality, we can insert the expansion (2.15) into the definition of in (2.10) to get . In our case, the normalization (2.11) gives us . Similarly, it can be deduced from the other two equations in (2.10) and (2.11) that

(2.19) |

Other interesting moments include the heat flux and the stress tensor , which are defined as

They are related to the coefficients by

and

## 3 Approximation of the collision term

To get the evolution of the coefficients in the expansion (2.15), we need to expand the collision term using the same basis functions. The expansions of the BGK-type collision operators are usually quite straightforward. For instance, the series expansion of the BGK collision term (2.13) is given in [6] as

(3.1) |

where

The expansions for the ES-BGK and Shakhov operators can be found in [5, 4]. In this section, we will first discuss the series expansion of the quadratic collision term defined in (2.2), and then mimic the BGK-type collision operators to construct collision models with better accuracy.

### 3.1 Series expansions of general collision terms

Suppose the binary collision term can be expanded as

(3.2) |

By the orthogonality of Hermite polynomials, we get

(3.3) |

where the second equality can be obtained by inserting (2.15) into (2.2), and

(3.4) | ||||

It can be seen from (3.4) that the evaluation of every coefficient requires integration of an eight-dimensional function. In principle, this can be done by numerical quadrature; however, the computational cost for obtaining all these coefficients would be huge. Actually, in [16, 30], a strategy to simplify the above integral has been introduced, and for small indices, the values are given in the literature. However, when the indices are large, no explicit formulae are provided in [16, 30], and the procedure therein is not easy to follow. Inspired by these works, we give in this paper explicit equations of the coefficients for any collision kernel, except for an integral with respect the two parameters in the kernel function . The main results are summarized in the following two theorems:

###### Theorem 1.

###### Theorem 2.

###### Definition 2 (Legendre functions).

For , the Legendre polynomial is defined as

###### Definition 3 (Laguerre polynomials).

For , let . For , define the Laguerre polynomial as

Through these two theorems, the eight-dimensional integration in (3.4) has been reduced into a series of summations and a two-dimensional integration. Among all the coefficients introduced in these theorems, and can be computed directly. As for , we need to expand polynomial , which can be done recursively using the following recursion formula:

(3.13) |

This recursion formula can be derived from the recursion relation of Legendre polynomials, and it also shows that for every monomial in the expansion of , the degree of equals the degree of . Therefore is nonzero only when . This means in (3.9), the summand is nonzero only when

(3.14) |

Consequently, when evaluating defined in (3.12), we only need to take into account the case . Generally, can be computed by numerical quadrature; for the IPL model, the integral with respect to can be written explicitly, which will be elaborated in the following section.

### 3.2 Series expansion of collision operators for IPL models

The formulae given in the previous section are almost ready to be coded, except that specific collision models are needed to calculate the integral defined in (3.12). This section is devoted to further simplifying this integral for IPL models, which completes the algorithm for computing the coefficients .

For the IPL model (2.5), we first consider the integral with respect to in (3.12). To this aim, we extract all the terms related to from (3.12), and define as

To evaluate the above integral, we follow the method introduced in [7] and apply the change of variable

to get

(3.15) |

Below we write the above equation as

(3.16) |

where denotes the integral in (3.15). In general, we need to evaluate by numerical quadrature. In our implementation, the adaptive integrator introduced in [27, Section 3.3.7] is used to compute this integral.

Now we consider the integral with respect to . Using the result (3.16), we can rewrite (3.12) as

(3.17) |

where , and we have applied the change of variable , and taken into account the relation . In general, we can adopt the formula

(3.18) |

introduced in [29, eq. (10)] to calculate (3.17). Specially, when , which corresponds to the model of Maxwell molecules, we can use the orthogonality of Laguerre polynomials to get

(3.19) |

and thus the computational cost can be further reduced. In fact, Grad has already pointed out in [16] that for Maxwell molecules, is nonzero only when

(3.20) |

This can also be seen from our calculation: from (3.19), we can find that only when , the coefficient given in (3.9) is nonzero; therefore in (3.5), if the summand is nonzero, the sum of , and must equal the sum of , and , which is equivalent to (3.20) due to (3.6).

The above analysis shows that for the IPL model, we only need to apply the numerical quadrature to the one-dimensional integrals , which makes it easier to obtain the coefficients with high accuracy.

### 3.3 Approximation of the collision term

Until now, we already have a complete algorithm to calculate the coefficients . These coefficients can be used either to discretize the collision term or to construct new collision models. We will discuss both topics in this section.

#### 3.3.1 Discretization of the homogeneous Boltzmann equation

Based on the expansion of the distribution function (2.15), the most natural discretization of the homogeneous Boltzman equation is to use the Galerkin spectral method. From this point of view, for any positive integer , we define the space of the numerical solution

(3.21) |

where is the index set

Then the semi-discrete distribution function satisfies

(3.22) |

Suppose

(3.23) |

The equations (3.2) and (3.3) show that the variational form (3.22) is equivalent to the following ODE system:

(3.24) |

It is easy to see that the time complexity for the computation of all the right hand sides is , where is the number of elements in :

(3.25) |

To fully formulate the ODE system (3.24), we need the coefficients for all . When the collision kernel is chosen and is fixed, we only need to compute these coefficients once, and then they can be used repeatedly. For a given , the algorithm for computing these coefficients is summarized in Table 1. The general procedure is to sequentially compute the coefficients in the first column, with indices described in the third column, and the equations to follow are given in the second column. For IPL models, we can use (3.17) and (3.18) instead to obtain the values of . In the third column of Table 1, it is worth mentioning that some indices are in the index set instead of , as is due to the equation (3.6), which shows that

Therefore the corresponding indices for and must lie in . Similar arguments hold for the coefficients .

The last column in Table 1 shows an estimation of the computational cost for each coefficient, from which one can see that the total cost for getting is . Now we compare this with the numerical cost by applying numerical integration directly to (3.4). We assume the number of quadrature points on is , and the number of quadrature points on the unit sphere (domain for and ) is . Thus using numerical integration to evaluate all the coefficients has time complexity . In most cases, we will choose to get accurate results. Hence our method listed in Table 1 is significantly faster.

#### 3.3.2 Approximation of the collision operator

In the previous section, a complete numerical method has been given to solve the spatially homogeneous Boltzmann equation. However, due to the rapid growth of the number of coefficients as increases, the storage requirement of this algorithm is quite strong. Table 2 shows the memory required to store the coefficients , where we assume that the coefficients are represented in the double-precision floating-point format, whose typical size is bytes per number. It can be seen that the case has already exceeded the memory caps of most current desktops. Although the data given in Table 2 can be reduced by taking the symmetry of the coefficients into consideration, it can still easily hit our memory limit by increasing slightly. Even if the memory cost is acceptable for large , the computational cost becomes an issue especially when solving the spatially non-homogeneous problems.

Memory (Gigabytes) | Memory (Gigabytes) | ||
---|---|---|---|

5 | 25 | ||

10 | 30 | ||

15 | 35 | ||

20 | 40 |