# A direct method for the Boltzmann equation based on a pseudo-spectral velocity space discretization

###### Abstract

A deterministic method is proposed for solving the Boltzmann equation. The method employs a Galerkin discretization of the velocity space and adopts, as trial and test functions, the collocation basis functions based on weights and roots of a Gauss-Hermite quadrature. This is defined by means of half- and/or full-range Hermite polynomials depending whether or not the distribution function presents a discontinuity in the velocity space. The resulting semi-discrete Boltzmann equation is in the form of a system of hyperbolic partial differential equations whose solution can be obtained by standard numerical approaches. The spectral rate of convergence of the results in the velocity space is shown by solving the spatially uniform homogeneous relaxation to equilibrium of Maxwell molecules. As an application, the two-dimensional cavity flow of a gas composed by hard-sphere molecules is studied for different Knudsen and Mach numbers. Although computationally demanding, the proposed method turns out to be an effective tool for studying low-speed slightly rarefied gas flows.

###### keywords:

Boltzmann equation, deterministic solution, Galerkin method, Gaussian quadrature, half- and full-range Hermite polynomials###### Pacs:

02.70.Bf, 47.45.Ab, 51.10.+y, Corresponding author.

## 1 Introduction

The conventional continuum approach to gas dynamics, namely the Navier-Stokes equations with no-slip boundary conditions, is justified when the average distance traveled by molecules between two successive collisions, , is much smaller than a characteristic length, , associated to the flow geometry. This condition breaks down in several physical situations ranging from the re-entry of spacecraft in upper planetary atmospheres, characterized by large , to fluid-structure interaction in small-scale micro-electro-mechanical systems, characterized by small . In such situations, a microscopic description of the gas based on the Boltzmann equation is required [1]. In the absence of external forces, the Boltzmann equation for a gas composed by a single monatomic species whose atoms have mass , takes the form

(1) |

where

(2) |

In Eqs. (1,2), denotes the distribution function of atomic velocities at spatial location and time , whereas gives the collisional rate of change of at the phase space point at time and we have used the shorthand , and . As is clear from Eq. (2), is a non-linear functional of , whose precise structure depends on the assumed atomic interaction forces through the differential cross section . The dynamics of binary encounters determines the differential cross section as a function of the modulus of the relative velocity of two colliding atoms and of the orientation of the unit impact vector with respect to . The collisional dynamics also determines the pre-collisional velocities and which are changed into and by a binary collision

(3) | |||||

(4) |

The prevalent approach for numerically solving Eq. (1) is a stochastic-based method called Direct Simulation Monte Carlo (DMSC) [2]. The distribution function is represented by a number of mathematical particles which move in the computational domain and collide according to stochastic rules derived from Eqs. (1)-(2). Macroscopic flow properties are obtained by time averaging particle properties. If the averaging time is long enough, then accurate flow simulations can be obtained by a relatively small number of particles. Variants of DSMC have been proposed over the years to improve solution accuracy in the presence of high density gradients [3] and for low Mach number flows [4]. Although particle-based methods are by far the most effective tools in describing non-equilibrium gas flows, they are not well suited to simulate unsteady gas flows. Indeed, in this case the possibility of time averaging is lost or reduced. Acceptable accuracy can only be achieved by increasing the number of simulation particles or superposing several flow snapshots obtained from statistically independent simulations of the same flow but, in both cases, the computing effort is considerably increased. The simulation of steady gas flows in the near continuum limit represents an additional challenge since the time scale on which the particle-based methods are forced to operate is much shorter than the characteristic macroscopic time and therefore explicit integration to steady state is computationally demanding. Approaches based on a direct discretization of the Boltzmann equation in the phase space are believed to be a feasible alternative in these cases since they provide solutions with high accuracy even in unsteady conditions and offer the possibility of a direct steady-state formulation [5, 6]. As such, they have been applied to study several problems of both theoretical interest and practical importance, including the viscous gas damping in microfluidic devices [7, 8], the onset of instability in a rarefied gas environment [9] and the investigation of ghost effects [10]. Deterministic methods of solution present some further assets compared to particle-based methods. Firstly, they are more suited to be adopted within a domain decomposition approach since the need to exchange information between kinetic and macroscopic equations requires smooth numerical solutions [12, 13]. Secondly, unlike particle-based methods their implementation on massively parallel computers with SIMD architecture, such as multi-core and GPU processors, can easily realize the full potential of these supercomputers [14, 15, 16]. These aspects also prompted the development of deterministic methods of solution.

Considerable progress has been accomplished in developing deterministic method for kinetic model equations [17, 18]. By contrast, an accurate and efficient direct solution of the Boltzmann equation itself remains a challenging problem. A common strategy adopted for solving Eq. (1) consists in decoupling the transport and the collision terms by time-splitting the evolution operator into a transport step and a collisional step. The transport step requires to solve a system of hyperbolic conservation laws coupled at the boundaries. Their discretization can be performed in a variety of ways, including finite-difference, finite-volume, finite-element or spectral methods [19]. The collision step consists of solving a spatially homogeneous relaxation equation. This is the more computationally demanding part since it involves the computation of the bilinear five-fold integral defining the collision operator, Eq. (2). The numerical approaches to evaluate the collision step may be grouped into two broad categories. To the first category belong methods referred to as discrete velocity models (DVM). They make use of a Cartesian grid in velocity space and construct a discrete collision mechanics on the points of the grid that preserve the main physical properties of the collision integral, namely equilibrium states, collision invariants and entropy inequality [20]. DVM methods have high computational cost and low order of accuracy although fast algorithms have been recently developed for a restricted set of collision kernels [21]. To the second category belong methods which adopt a Galerkin discretization of the velocity space. They are based on expanding the velocity dependence of the distribution function in a set of trial functions with expansion coefficients that depend on position and time. The Galerkin ansatz is substituted in the space homogeneous relaxation equation which is subsequently multiplied by test functions and integrated in the velocity space. According to the Galerkin approach, test and trial functions are assumed to be the same. The above procedure yields a system of hyperbolic partial differential equations for the expansion coefficients. Galerkin methods can be further distinguished depending on the basis functions which they employ. In Fourier-Galerkin approach, the distribution function is expanded in trigonometric polynomials and the fast Fourier Transform is used to accelerate the computation of the collision integral in the velocity space. Several different methods have been developed starting from different representation of the collision integral [22, 23, 24, 25]. These methods are generally very efficient and spectrally accurate for smooth solutions. Their major shortcoming is the loss of some of the properties of the solution such as positivity and conservation of momentum and energy. Preservation of collision invariants can be enforced but the use of corrections procedure may limit the accuracy of the solutions. Discontinuous Galerkin methods adopt discontinuous piecewise polynomials as test and trial functions [26, 27, 28, 29]. Although computationally demanding, these methods have the remarkable feature to provide spectral accuracy in the velocity space even for discontinuous solutions which typically occur in the presence of solid surfaces. An hybrid approach is adopted in Refs. [30, 31] where the distribution function is expanded in Laguerre polynomials with respect to the velocity components parallel to solid surfaces whereas quadratic finite element functions have been used for the normal velocity component. This approach permits to explicitly account for the discontinuity of the distribution function and its application to one-dimensional problems has been shown to provide accurate numerical solutions in both linearized and non-linear regime.

The main objective of the present work is to propose a new approach to the deterministic solution of the Boltzmann equation.
The method uses a Galerkin discretization of the velocity space and adopts, as trial and test functions,
the collocation basis functions based on weights and roots of a Gauss-Hermite quadrature.
The Boltzmann equation is thus simplified to a system of hyperbolic conservation laws with non-linear source terms in the physical space.
A fully discrete numerical scheme can then be derived by using standard approaches [19].
The method proposed in Ref. [32] shows strong similarities with the present work but it differs in two significant respects.
In Ref. [32], the method is developed for a system of linear Boltzmann equations and
the collocation basis functions have been defined only on the basis of full-range Hermite polynomials.
By contrast, we are here concerned with the non-linear Boltzmann equation and the roots and weights of the half-range Gauss-Hermite
quadrature formula also enter in the definition of the collocation basis functions.
Half-range Hermite polynomials have been widely used to study flow and heat transfer problems in rarefied gas dynamics [33, 34, 35, 36] and
the Milne problem in radiative transfer [37]. In comparison with full-range Hermite polynomials,
they allow to deal exactly with boundary conditions and to achieve a faster convergence rate in the presence of discontinuous distribution functions.

The rest of the paper is organized as follows.
Section 2 illustrates the method of solution in detail.
Section 3 covers two aspects. Firstly, the spectral accuracy of the proposed velocity space discretization is illustrated by a comparison with the
exact solution of the spatially homogeneous relaxation for Maxwellian molecules.
Secondly, the possibilities of the proposed method are demonstrated by solving the two-dimensional driven cavity flow of a gas composed by hard-sphere molecules for
different Knudsen and Mach numbers.
Section 4 summarizes the main results and the future research directions.

## 2 Mathematical formulation

The numerical scheme to deterministically solve the Boltzmann equation is derived through a two-step procedure. As a first step, Eq. (1) is rewritten in terms of the deviational part of the distribution function, , defined as

(5) |

where is the Maxwellian at equilibrium with uniform and constant density and temperature , i.e.,

(6) |

As discussed in Refs. [15, 26, 28], the use of permits to greatly improve the accuracy of the method of solution when the gas approaches equilibrium since it reduces the truncation errors which arise in evaluating Eq. (2) for distribution functions close to Maxwellian. In this respect, an equilibrium distribution function which depends on space and time, , may be also adopted in Eq. (5) to further increase accuracy of the computations. The applicability of the method, however, is not affected by this choice. Substituting Eq. (5) into Eq. (1) yields

(7) |

It is worth noticing that if the perturbation is sufficiently small, i.e., ,
the quadratic terms in become negligible and Eq. (7)
simplifies to the linearized Boltzmann equation.
The present formulation in terms of the deviational part of the distribution function, however,
is not restricted to a vanishing perturbation and it is valid in the non-linear case as well.

As a second step, we seek an approximate solution of Eq. (7) by using a Galerkin approach.
A suitable finite dimensional linear space of orthonormal trial functions
is chosen and the deviational part of the distribution function is expanded in a series of the form

(8) |

Inserting Eq. (8) into Eq. (7) implies a residual depending on the expansion coefficients, . These are determined by requiring that the residual is orthonormal with respect to some test functions. According to the Galerkin method, trial and test functions are assumed to be equal. Therefore, we substitute the spectral Galerkin ansatz, Eq. (8), into Eq. (1), multiply both side by and integrate with respect to the velocity variable . We thus obtain a closed first-order hyperbolic system of conservation laws

(9) |

In Eq. (9), the coefficient matrices , and are given, respectively, by

(10) | |||||

(11) | |||||

(12) |

The time-independent kernel of the integral operators in Eqs. (11) and (12), , read,

(13) |

Boundary conditions can be treated similarly. For the sake of simplicity, let us consider a fixed plane wall at with temperature . Let us further suppose that the gas fills the half space and molecules which strike the wall are re-emitted according to the Maxwell’s scattering kernel with complete accommodation [1]. In terms of the deviational distribution function, the boundary condition reads

(14) |

Substituting Eq. (8) into Eq. (14), multiplying both side by and integrating with respect to the velocity variable , yields

(15) |

where

(16) |

Several methods of solution for the Boltzmann equation employ a Galerkin discretization of the velocity space but very different numerical schemes are derived due to the different trial and test functions [22, 23, 24, 25, 26, 27, 28, 29, 30, 31]. In the present work, we adopt the collocation basis functions based on weights, , and roots, , of the Gaussian quadrature formula defined by means of half- and/or full-range Hermite polynomials, i.e., with as given by Eq. (28) (see Appendix A for further details). A first consequence of this choice is that, by virtue of Eq. (34), expansion coefficients are proportional to the values of the deviational part of the distribution function at the quadrature roots , that is . The proposed method of solution is thus closely related to discrete velocity methods. Furthermore, computing the integrals over the velocity space in Eqs. (10)-(12) as finite summations according to the Gauss-Hermite formula for numerical quadrature, Eq. (26), yields

(17) | |||

(18) | |||

(19) |

Likewise, integrals in Eqs. (16) can be rewritten as

(20) |

Matrices given by Eqs. (18),(19) and (20), can be computed with good accuracy once and for all and
used in many individual simulations as long as their velocity space discretization is the same.
It is worth noticing that Gauss-Hermite quadrature formula can integrate exactly polynomials of degree at most in each velocity component .
Since integrals in Eqs. (11) and (12) may involve polynomials of higher order depending on the differential cross section,
their numerical evaluation is approximated.
The proposed method for the numerical solution of the Boltzmann equation has therefore a pseudo-spectral nature.

The collocation basis functions introduced above hold several advantages compared with alternative choices of functions .
Firstly, as shown by Eq. (17), they lead to a diagonal system of conservation laws. This simplifies the problem of solving the
advection step since boundary conditions are explicitly identified.
Secondly, using half-range Hermite polynomials permits to exactly compute the half-space integrals appearing in the boundary conditions,
Eq. (16). This permits to greatly improve accuracy of the solution for boundary value problems [33, 35].
It is worth pointing out that the use of half-range Hermite polynomials is not restricted to a cartesian geometry since, by making use of a
body fitted coordinate system, arbitrary geometries can be dealt with as well [38].
Finally, as it will be shown in Subsec. 3.1, collision invariants are preserved during the relaxation step
within machine precision and no correction method to enforce conservation of momentum and energy is thus required.

## 3 Results and discussion

The key feature of the method described in Section 2 is the introduction of a grid in the velocity space
whose points are the roots of the Gaussian quadrature formula based on half- and/or full-range Hermite polynomials.
This discretization of the velocity space has major impact on the evaluation of the collision integral.
In Subsec. 3.1 the proposed method is thus assessed by computing spatially homogeneous relaxation for Maxwell molecules.
Most of the spectral approaches for solving the Boltzmann equation have been developed for
a restricted set of collision kernels which do not include Maxwell molecules in the three-dimensional velocity space. By contrast, the
method proposed here can be carried over to arbitrary particle interactions and therefore we can consider
the fully three-dimensional relaxation process.
The comparison with an exact solution will show that mass, momentum and energy are conserved and numerical results are spectrally accurate
in the velocity space.

In Subsec. 3.2 the possibilities of the proposed method are demonstrated by solving the driven cavity flow of a gas composed by
hard-sphere molecules. This is a classical benchmark problem which it is often used for the validation of numerical codes.
In spite of its simple geometry, in fact, it contains most of the features of more complicated problems described by kinetic equations.
A comprehensive study of this problem has been carried out in Refs. [44] for the Bhatnagar-Gross-Krook-Welander kinetic model
equation in the linearized regime. By contrast, only few results are available for the hard-sphere Boltzmann equation [15, 45].
Subsec. 3.2 is thus of interest in that it provides accurate reference solutions albeit for a restricted range of Knudsen and Mach numbers.

### 3.1 Spatially homogeneous relaxation

The method described in Sec. 2 is here assessed by comparing the numerical solution of the three-dimensional spatially homogeneous relaxation of Maxwell molecules with the exact Bobylev-Krook-Wu (BKW) distribution [40, 41, 42]. The BKW distribution is radially symmetric and reads

(21) |

where and are the equilibrium density and temperature, respectively, and

(22) |

being the constant which enters in the definition of the Maxwell’s molecules differential cross section, i.e., . The distribution function given by Eq. (21) has physical validity only for since it is negative for smaller times. The three-dimensional velocity space is discretized using discrete velocities. Since the BKW distribution function is continuous in the velocity space, the discrete velocities are the quadrature nodes only based on the full-range Hermite polynomial of degree . A simple first order explicit Euler method is employed to solve Eq. (9) in which transport terms are neglected, , and the coefficient matrices, Eqs. (18) and (19), are computed using the differential cross section of Maxwell’s molecules. The initial condition is determined by expanding Eq. (21) at according to Eq. (8). The ability of the polynomial expansion to reproduce the initial condition is shown in Fig. 1 where a comparison is made between the BKW distribution evaluated at (solid line) and the distribution function as given by Eq. (8) (dashed lines) for different number of expansion polynomials used in each direction of the velocity space . Because of the symmetry of the BKW distribution, odd full-range Hermite polynomials do not effectively contribute to the expansion defined in Eq. (8) and thus only results for the even values of have been reported. A relatively high number of polynomials is required to obtain an accurate representation of the initial condition. This is reasonable since the BKW solution at is an highly non-equilibrium distribution function. However, we notice that even if a low number of expansion polynomials is used, i.e., , the distribution functions given by Eqs. (8) and (21) provide the same mass, momentum and energy. As a matter of fact, by construction, the first moments of and are the same. In order to quantitatively determine the convergence rate of the polynomial expansion, we report in Fig. 2 the error between BGW and the numerical distribution function for two different initial instants of time, and , versus the total number of discrete velocities, . The error has been computed according to the formula

(23) |

where is the BKW distribution function, Eq. (21), and is the numerical solution. As Figure 2 clearly shows, for a given , the error at is lower than the error at , which is not unexpected. Indeed, the higher is the closer is the distribution function to a Maxwellian and therefore the smaller is the number of expansion polynomials required to represent the distribution function with a given accuracy. Furthermore, independently of the initial instant of time, as increases, the error decreases with a rate which is higher than any algebraic order. The proposed velocity space discretization is therefore spectrally accurate. In Fig. (3) the error is reported as a function of time for and with . The most important feature to highlight is that the error asymptotically tends to zero up to machine precision whatever is the value of . This behavior indicates that distribution function approaches the Maxwellian distribution while conserving mass, momentum and energy.

### 3.2 Two dimensional driven cavity flow

Ma | |||||

The method described in Sec. 2 is here applied to solve the two-dimensional driven cavity flow of a gas composed
by hard-sphere molecules. The problem geometry consists of a square enclosure with length .
All the walls are kept at uniform and constant temperature .
Initially, the gas is in uniform equilibrium with density and temperature .
The flow is driven by the translation of the lid of the cavity with constant velocity [44].

The cavity flow problem has been solved for Mach numbers in the range and for two values of the Knudsen number, and .
The Mach and Knudsen numbers are defined, respectively, as with the speed of sound, and with
and the viscosity of the hard sphere gas [1].
Numerical results have been obtained by solving Eq. (9) where the coefficient matrices, Eqs. (17)-(19),
are computed with the hard-sphere differential cross section, i.e., , being the molecules’ diameter.
The solution in one time step has been obtained by first integrating the transport equation, Eq. (9) with ,
and afterwards the space homogeneous equation, Eq. (9) with , using the output of the previous step as initial condition.
Iterating the process provides the numerical solution at later times.
The three-dimensional velocity space is discretized using discrete velocities.
Because of the presence of solid walls, the discrete velocities along the and directions are the quadrature nodes of the positive and negative
-order half-range Hermite polynomials whereas discrete velocities along the direction are the the quadrature nodes of the -order
full-range Hermite polynomial.
The two-dimensional physical space has been divided into rectangular cells. Their lengths have been uniformly stretched with
a progression ratio for and for .
The smaller cells are located close to the upper corners of the cavity where severe gradients are anticipated.
The collision and transport steps have been solved by using simple first-order explicit Euler and upwind schemes.

We firstly establish the convergence rate of the method by computing two global flow field properties, namely the drag coefficient, , and the flow rate of the main vortex, . We refer to Refs. [15, 46] for precise definitions of all the macroscopic quantities used in the present Section. The absolute relative error in the long-term dimensionless values of and are shown in Figs 4a and 4b, respectively, versus the dimensionless spatial grid size, , for two values of the Knudsen number, (solid symbols) and (empty symbols), and two values of the Mach number, (squares) and (circles). It is worth noticing that the scheme which has been used for numerically solving the Boltzmann equation in the linearized regime, that is for , is provided by Eq. (9) in which the second term on the right hand side is disregarded. The exact values of and , which are referred to as and , have been extrapolated from the linear fit of the results when . The linear behavior of the absolute relative errors demonstrates that the results are in the asymptotic range of convergence and the method is overall first order accurate in the physical space [47]. Table 1 shows the numerical predictions of the volume flow rate of the main vortex , the drag coefficient, , the mean temperature, , and the -component of the mean heat flux, , along the moving wall. The error bound for and is also reported for and . The largest percentage errors in and are about and and they are attained at in the linearized regime. Results for intermediate Mach numbers are expected to retain a similar level of accuracy.

Figures 5 and 6 show the profiles of the dimensionless horizontal component of the velocity, , along the vertical line crossing the center of the cavity and the dimensionless vertical component of the velocity, , along the horizontal line crossing the center of the cavity, respectively, for and Mach numbers in the range . The possibilities of the proposed method of solution are further illustrated by the computations reported in Figs. 7 and 8. Fig. 7 shows snapshots of the heat flux lines superimposed to the temperature field for two values of the Mach number, (left panel) and (right panel) whereas in Fig. 8 the time evolution of the drag coefficient, , is reported during the simulation for and Mach numbers in the range . These results would be difficult to obtain with a particle-based method since computationally expensive time and/or ensemble averaging are needed to provide such smooth macroscopic fields.

We conclude with a few remarks on the parallel implementation of the proposed method of solution. Profiling of the serial code has highlighted that the collision step has a computational cost which is four order of magnitude greater than the computational cost of the streaming step and it takes about of the total simulation time. Therefore, only the collision step has been parallelized leaving the transport step performed by all MPI processes. Matrix-vector operations involved in the collision step are independent along , and indices. As parallelization strategy we decided to implement an hybrid solution decomposing the simulation along the index with OpenMP and the index with MPI. Simulations have been performed on a cluster which is comprised by nodes, each node having cores Intel Xeon X5660 and of RAM. The nodes are interconnected by an Infiniband QDR/DDR Voltaire network. Benchmarks showed that the parallelization strategy which has been adopted yields an efficiency greater than up to 128 cores. The results reported in this Section were obtained using 64 cores, with a memory requirement for each simulation of about of RAM and a computational cost of the order of floating point operations per time iteration, being and the total numbers of cells in the physical and velocity space, respectively. The time required to execute one time iteration was about .

## 4 Conclusion

We developed a deterministic method of solution for the Boltzmann equation based on a pseudo-spectral Galerkin discretization of the velocity space.
The novel aspect of the proposed method is the choice of trial and test functions, namely
the collocation basis functions related to weights and roots of a Gauss-Hermite quadrature formula defined on the basis of half- and/or full-range Hermite polynomials,
Eqs. (28).
The use of half-range Hermite polynomials is an important feature of the method since it permits a straightforward and exact formulation of boundary
conditions.
Due to this discretization of the velocity space, the Boltzmann equation simplifies to a system of hyperbolic conservation laws with
non-linear source terms in the physical space, Eqs. (9). A fully discrete numerical scheme is derived by a first-order time splitting
of the evolution operator which couples a first-order finite-volume scheme for the transport step with a first-order explicit Euler scheme for the relaxation step.
However, more sophisticated methods can be used for coupling and solving both transport and relaxation steps.
Numerical results have been presented for the space homogeneous relaxation to equilibrium of Maxwell molecules
and the two-dimensional cavity flow of a gas composed by hard-sphere molecules.
Overall, the numerical results are affected by an error of less than a few percent.

The proposed method of solution has a number of attracting features such as explicit identification and exact treatment of boundary conditions,
spectral accuracy in the velocity space and preservation of the collision invariants.
While the method can in principle handle a wide range of Knudsen and Mach numbers, it turns out to be feasible
only in dealing with slightly rarefied low Mach number gas flows.
The main obstacle is the large number of operations required to evaluate the collision integral.
However, a preliminary study indicates that the numerical code possesses very good scalability
on a memory distributed cluster of GPUs and supercomputing architectures like IBM Blue Gene/Q.
A throughly discussion of these parallel implementations will be the subject of a future publication.

## 5 Acknowledgments

The authors would like to thank Prof. A. Frezzotti for several helpful discussions and for critically reading the paper. Moreover, they gratefully acknowledge the support received from Regione Lombardia and SAES Getters within the framework of the projects “Modelli Matematici per la Progettazione di Barriere Attive (MOBA)”. This work has also been supported by Regione Lombardia and CILEA Consortium through a LISA Initiative (Laboratory for Interdisciplinary Advanced Simulation) 2011 grant [link:http://lisa.cilea.it].

## Appendix: Half-range Gauss-Hermite quadrature

We introduce polynomials in three variables defined as the Cartesian product of polynomials in one variable and we organize them according to the lexicographical order, i.e., where , with and the number of polynomials in the velocity component. By applying the Gram-Schmidt procedure to the dimensionless monomial powers , polynomials can be constructed that are orthonormal in the range with respect to the weighted inner product defined as

(24) |

where denotes the Kronecker delta function which is for and otherwise. In Eq. (24), is the Maxwellian weight function with unit mass which reads

(25) |

As a matter of fact, polynomials can be regarded as defined over the whole real axis by setting their values to zero
outside the interval .

The polynomials afore introduced permit to define the following quadrature formula

(26) |

where are the roots of the polynomial and are quadrature weights which are determined by solving the linear systems of equations

(27) |

where .
The quadrature formula given by Eq. (26) can be shown to provide the correct result
for the integral of functions that are products of polynomials in each velocity component up to degree ,
as long as their support is contained in .

In the present work, three polynomial sets are used that satisfy Eq. (24) in a different interval .
The support of full-range Hermite polynomials, , is
whereas the support of negative, , and positive, , half-range Hermite polynomials are and
, respectively.
Zeros and weigths of and are referred to as
and , respectively.
It is not difficult to prove that the zeros of full-range Hermite polynomials are symmetric
about the zero velocity, i.e., ,
whereas the zeros of positive and negative half-range Hermite polynomials have equal magnitude but different sign, i.e., .

The collocation basis functions introduced in Section 2 are
based on roots and weights of the Gauss-Hermite quadrature formula defined by Eq. (26) for the polynomial set
.
More precisely, the collocation basis functions read

(28) |

where and are the nodes and weights of the quadrature formula given by Eq. (26) based on the polynomials . These are defined differently depending whether or not the deviational part of the distribution function is expected to have a discontinuity in the velocity component. More precisely, if is continuous in the velocity space then and , . Instead, if presents a discontinuity in the velocity component then

(29) |

Likewise

(30) |

and

(31) |

Equation (28) can be rewritten with a more compact notation as follows

(32) |

By using the quadrature formula given by Eq. (26), it is possible to prove that the collocation basis functions defined by Eq. (28) are orthonormal with respect to the Maxwellian weight factor,

(33) |

and that the -th basis function vanishes on all nodes but node

(34) |

For illustrative purposes, in Figs 9, we report the collocation basis functions in one-variable based on the quadrature nodes of the full-range Hermite polynomials (left panel) and positive half-range Hermite polynomials (right panel).

## References

- [1] C. Cercignani, The Boltzmann Equation and Its Applications, Springer-Verlag, New Yor1Gk, 1988.
- [2] G. A. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Oxford University Press, 1994.
- [3] S. Rjasanow, W. Wagner, “A stochastic weighted particle method for the Boltzmann equation”, J. Comput. Phys. 124 (1996) 243-253.
- [4] T. M. M. Homolle , N. G. Hadjiconstantinou, “A low-variance deviational simulation Monte Carlo for the Boltzmann equation”, J. Comput. Phys. 226 (2007) 2341-2358.
- [5] V. V. Aristov, Direct Methods for Solving the Boltzmann Equation and Study of Nonequilibrium Flows, Kluwer Academic Publishers, 2001.
- [6] P. Degond, L. Pareschi and G. Russo, Modeling and Computational Methods for Kinetic Equations, Series: Modeling and Simulation in Science, Engineering and Technology, Birkhäuser, 2004.
- [7] S. Lorenzani, L. Gibelli, A. Frezzotti, A. Frangi, C. Cercignani, “Kinetic approach to gas flows in microchannels”, Nanos. Micros. Therm. 11 (2007) 211-226.
- [8] C. Cercignani, A. Frangi, A. Frezzotti, G. P. Ghiroldi, L. Gibelli, and S. Lorenzani, “On the application of the Boltzmann equation to the simulation of fluid structure interaction in micro-electro-mechanical-systems”, Sens. Lett. 6 (1) (2008) 121-129.
- [9] V. V. Aristov, O. I. Rovenskaya, “Application of the Boltzmann kinetic equation to the eddy problems”, Comput. Fluids. 50 (2011) 189-198.
- [10] Y. Sone, T. Doi, “Ghost effect of infinitesimal curvature of in the plane Couette flow of a gas in the continuum limit”, Phys. Fluids 16 (2004) 952-971.
- [11] Y. Sone, Kinetic Theory and Fluid Dynamics (Birkhäuser, Boston, 2002).
- [12] P. Degond, G. Dimarco, L. Mieussens, “A multiscale kinetic-fluid solver with dynamic localization of kinetic effects”, J. Comput. Phys. 229 (2010) 4907-4933.
- [13] A. Alaia, G. Puppo, “A hybrid method for hydrodynamic-kinetic flow. Part I: A particle-grid method for reducing stochastic noise in kinetic regimes”, J. Comput. Phys. 230 (2011) 5660-5683.
- [14] A. Frezzotti, G. P. Ghiroldi, L. Gibelli, “Solving model kinetic equations on GPUs”, Comput. Fluids. 50 (2011) 136-146.
- [15] A. Frezzotti, G. P. Ghiroldi, L. Gibelli, “Solving the Boltzmann equation on GPUs”, Comput. Phys. Commun. 182 (2011) 2445-2453.
- [16] J. R. Haack, I. M. Gamba, “High performance computing with a conservative spectral Boltzmann solver”, 28th International Symposium on Rarefied Gas Dynamics, 9-13 July 2012, Zaragoza, Spain, AIP Conference Proceedings American Institute of Physics, 1501, 334-341 (2012).
- [17] L. Mieussens, “Discrete-velocity models and numerical schemes for the Boltzmann-BGK equation in plane and axisymmetric geometries”, J. Comput. Phys. 162 (2000) 429-466.
- [18] A. M. Alekseenko, “Numerical properties of high order velocity solutions to the BGK kinetic equation”, App. Num. Math. 61 (2011) 410-427.
- [19] R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002.
- [20] V. A. Panferov, A. G. Heintz, “A new consistent discrete-velocity model for the Boltzmann equation”, Math. Method. Appl. Sci. 25 (7) (2002) 571-593.
- [21] C. Mouhot, L. Pareschi, T. Rey, “Convolutive decomposition and fast summation methods for discrete-velocity approximations of the Boltzmann equation”, arxiv:1021.3986v2
- [22] L. Pareschi, G. Russo, “Numerical solution of the Boltzmann equation. I. Spectrally accurate approximation of the collision operator”, SIAM J. Numer. Anal. 37 (4) (2000) 1217-1245.
- [23] C. Mouhot, L. Pareschi, “Fast algorithms for computing the Boltzmann collision operator”, C. R. Acad. Sci. Paris Sér I Math 339 (1) (2004) 71-76.
- [24] F. Filbet, C. Mouhut and L. Pareschi, “Solving the Boltzmann equation in ” SIAM J. Sci. Comput. 28 (3) (2006) 1029-1053.
- [25] I. M. Gamba, S. H. Tharkabhushanam, “Spectral-Lagrangian methods for collisional models of non-equilibrium statistical states”, J. Comput. Phys. 228 (2009) 2012-2036.
- [26] L. L. Baker, N. G. Hadjiconstantinou, “Variance-reduced Monte Carlo solutions of the Boltzmann equation for low-speed gas flows: A discontinuous Galerkin formulation”, Int. J. Numer. Meth. Fl. 58 (2008) 381-402.
- [27] A. Majorana, “A numerical model of the Boltzmann equation related to the discontinuous Galerkin method”, Kin. Rel. Mod. 4 (1) (2011) 139-151.
- [28] A. Alekssenko, E. Josyula, “Deterministic solution of the Boltzmann equation using a discontinuous Galerkin velocity discretization”, 28th International Symposium on Rarefied Gas Dynamics, 9-13 July 2012, Zaragoza, Spain, AIP Conference Proceedings American Institute of Physics, 1501, pp. 279-286 (2012).
- [29] A. Alekssenko, E. Josyula, “Deterministic solution of the Boltzmann equation using discontinuous Galerkin discretization in velocity space”, airxiv:1301.1099v1
- [30] T. Ohwada, “Structure of normal shock wave: Direct numerical analysis of the Boltzmann equation for hard-sphere molecules”, Phys. Fluids 5 (1) (1992) 217-234.
- [31] Y. Sone, T. Ohwada, K. Aoki, “Temperature jump and Knudsen layer in a rarefied gas over a plane wall: Numerical analysis of the linearized Boltzmann equation for hard-sphere molecules”, Phys. Fluids 1 (1989) 363-370.
- [32] M. K. Gobbert, S. G. Webster and T. S. Cale, “A Galerkin method for the simulation of the transient 2-D/2-D and 3-D/3-D linear Boltzmann equation”, J. Sci. Comput. 30 (2) (2007) 237-273.
- [33] E. P. Gross, E. A. Jackson, S. Ziering, “Boundary value problem in kinetic theory of gases”, Ann. Phys. 1 (1957) 141-167.
- [34] C. E. Siewert, “The linearized Boltzmann equation: Concise and accurate solutions to basic flow problems”, Z. Angew. Math. Phys. 54 (2003) 273-303.
- [35] A. Frezzotti, L. Gibelli, B. Franzelli, “A moment method for low speed flows”, Cont. Mech. Thermodyn. 21 (6) (2009) 495-509.
- [36] L. Gibelli, “Velocity slip coefficients based on the hard-sphere Boltzmann equation”, Phys. Fluids 24 (2012) 022001.
- [37] M. J. Lindenfeld, B. D. Shizgal, “The Milne problem: A study of the mass dependence”, Phys. Rev. A 27 (3) (1983) 1657-1670.
- [38] Z. Li, H. Zhang, “Gas kinetic algorithm using Boltzmann model equation”, Computers & Fluids 33 (2004) 967â991.
- [39] J. Fan and C. Shen, “Statistical simulation of low-speed rarefied gas flows”, J. Comput. Phys. 167 (2001) 393-412.
- [40] A. V. Bobylev, “Exact solutions of the Boltzmann equation”, Dokl. Akad. Nauk. S. S. S. R. 225 (1975) 1296-1299 (Russian).
- [41] A. V. Bobylev, “The theory of the nonlinear spatially uniform Boltzmann equation for Maxwell molecules”, Math. Phys. Reviews 7 (1988) 111â233.
- [42] M. Krook, T. T. Wu, “Exact solutions of the Boltzmann equation”, Phys. Fluids 20 (10) (1977) 1589-1595.
- [43] F. Filbet, G. Russo, “High order numerical methods for the space non homogeneous Boltzmann equation”, J. Comput. Phys. 186 (2003) 457-480.
- [44] S. Naris, D. Valougeorgis, “The driven cavity flow over the whole range of the Knudsen number”, Phys. Fluids 17 (2005) 097106.
- [45] S. Mizzi, D. R. Emerson, S. K. Stefanov, R. W. Barber, J. M. Reese, “Effects of Rarefaction on Cavity Flow in the Slip Regimeâ, J. Comput. Theor. Nano. 4 (4) (2007) 817-822.
- [46] G. P. Ghiroldi, L. Gibelli, P. Dagna, A. Invernizzi, “Linearized Boltzmann Equation: A Preliminary Exploration of its Range of Applicability”, 28th International Symposium on Rarefied Gas Dynamics, 9-13 July 2012, Zaragoza, Spain, AIP Conference Proceedings American Institute of Physics, 1501, pp. 735-741 (2012).
- [47] M. D. Salas, “Some observations on grid convergence”, Comp. Fluids 35 (2006) 688-692.