An Efficient NRxx Method for Boltzmann-BGK Equation

An Efficient NR Method for Boltzmann-BGK Equation

Zhenning Cai,   Ruo Li   and Yanli Wang School of Mathematical Sciences, Peking University, Beijing, China, email:, LMAM & School of Mathematical Sciences, Peking University, Beijing, China, email: of Mathematical Sciences, Peking University, Beijing, China, email:

In [3], we proposed a numerical regularized moment method of arbitrary order (abbreviated as NR method) for Boltzmann-BGK equation, which makes numerical simulation using very large number of moments possible. In this paper, we are further exploring the efficiency of NR method with techniques including the 2nd order HLL flux with linear reconstruction to improve spatial accuracy, the RKC schemes to relieve the time step length constraint by the regularization terms, and the revised Strang splitting to calculate convective and diffusive terms only once without loss of accuracy. It is validated by the numerical results that the overall efficiency is significantly improved and the convergence order is kept well.

1 Introduction

In 1949, the moment method was introduced by Grad [5] as a technique to approximate the Boltzmann equation in a macroscopic view. With proficient mathematical skills, he derived the famous 13-moment system, but the system is problematic: the hyperbolicity can only be obtained in the neighbourhood of Maxwellian, and the structure of shock wave is non-smooth when the Mach number is large (see e.g. [13, 6]). During a long time, the moment method suffers lots of criticism, and very few progresses are made before 1990s. In the recent twenty years, a number of new ideas based on the moment method have come out, such as the Jin-Slemrod regularization of the Burnett equations [9], the COET (Consistently Ordered Extended Thermodynamics) method [12], the order of magnitude approach [15], and et al. In [17], regularization based on 1st order Chapman-Enskog expansion was considered to remedy the defects of Grad’s 13 moment equations. The regularized system obtained therein was referred as R13 equations. Soon, the R20 and R26 equations are studied [11, 7] and the aspiration to extend this method for system with more moments [16] was called as R by the authors of [17]. In [3], we proposed a numerical regularized moment method of arbitrary order for Boltzmann-BGK equation. The method is abbreviated as “NR method” later on for convenience following the tradition. The NR method makes it possible to investigate large moment systems, by direct numerical approximation of the regularized moment systems without deriving the explicit forms of the moment equations. The regularization in [3] is similar to the original derivation of the R13 equations [17]. And later in [4], the idea of COET [12] is adopted to revise the regularization terms, which results in a parabolic system.

Though one can escape from deriving the complex moment system, and the developing of the simulation program is greatly simplified by the NR method, the numerical efficiency of the method in [3] should be further explored. It has been verified that the computational time is linear in the number of moments. However, a 1st order HLL numerical flux, which is over diffusive, was used in [3] for the transportation part, such that we have to use quite fine spatial grids to reduce the numerical error to a moderate level. Noting that the regularized moment system is parabolic, the time step is quadratic in the grid size. It turns out that in the 1D case, the total computational cost is cubic in the number of spatial grids. With a quite fine spatial grid to achieve enough accuracy, the numerical simulation therein is rather time consuming.

In this paper, we focus on improving the efficiency of the NR method in [3] with some state-of-art techniques available to us. Precisely, three numerical techniques are employed:

  1. Linear spatial reconstruction. The piecewise linear spatial reconstruction is able to provide a 2nd order HLL flux and greatly reduces the numerical diffusion. Owing to the absence of analytical expressions of the moment equations, the reconstruction needs to be done carefully. A conservative reconstruction is proposed, and it is numerically verified to have achieved high resolution. Although the final scheme is still of first order, the numerical error is greatly reduced. Thus much less grids are used in the computation.

  2. RKC time integration scheme. The RKC time integration scheme [21] is adopted to enlarge the time step sizes. The RKC method is a Runge-Kutta type method originally designed for diffusive PDEs to provide large stability regions, while the regularized moment equations are convective-diffusive problems, where imaginary parts appear in the eigenvalues of the right hand side of the semi-discrete system. As is well known, the RKC schemes contain a parameter called as the damping factor, which allows a small imaginary perturbation of eigenvalues and is usually selected as a small positive value. In order to make the RKC scheme compatible with the convection terms, we use a large damping factor in our scheme to ensure stability. The final time step size is equivalent to the grid size, and the number of internal time steps is equivalent to the square root of the grid number. Thus the total number of time steps is essentially reduced, and no instability phenomenon apprears in our numerical test.

  3. Revised Strang splitting method. With enlarged time step, the numerical accuracy can be harmed by convection-collision splitting. Therefore, the Strang splitting is utilized to win higher order of accuracy in the time direction. Usually, the implementation of the Strang splitting requires twice calculations of the collision term in one time step. For this problem, we combine the two steps of collision in the successive time steps, so only once calculation of the collision term is performed in a time step.

Several numerical examples are carried out to show the efficiency of our algorithm. Our prediction of the computational cost together with the order of convergence is validated by numerical experiments. We also compare the current method with the scheme without linear reconstruction to demonstrate much higher resolution for the one with linear reconstruction.

The rest of this paper is arranged as follows: in Section 2, we give a brief review of the Boltzmann-BGK equation and the NR method. In Section 3, the linear reconstruction is added to the HLL numerical flux. In Section 4, the time step size is enlarged by the RKC method, and the accuracy is improved by the Strang splitting method. In Section 5, numerical examples are carried out to make illustrations. Some concluding remarks are given in Section 6.

2 The NR method

2.1 The Boltzmann-BGK model

In the kinetic theory, it is generally accepted that the Boltzmann equation is able to describe the fluids accurately. Due to the complexity of the collision term, in the computational field, several simplified collision operators are adopted, among which the BGK model [1] is the simplest but useful. The Boltzmann-BGK equation reads


where is the distribution function, is the collision frequency, and is the local Maxwellian defined by


Here , and are local macroscopic variables which represent the density, velocity and temperature respectively. They can be calculated from the distribution function by


2.2 A first-order scheme for the NR method

The NR method is raised in [3, 4] as a new tool for the computation of large moment systems. It is based on an Hermite expansion of the distribution function:


where is a -dimensional multi-index. ’s and ’s act as the basis functions and the corresponding coefficients respectively, and is defined by


where ’s are the orthogonal Hermite polynomials defined by


This idea originates from [5] where the 13-moment equations are derived. In order to get a finite system, the NR method [4] chooses an positive integer and approximate with by


Now that all ’s with have been related with ’s with , we can get a closed moment system of all ’s with by putting the expansion (2.6) into the BGK equation (2.1).

The explicit expressions of such moment systems can be written in a uniform style for any choice of , but those expressions are not convenient for computation, since they are not in a conservative form. In order to avoid an intricate process to obtain balance laws, the NR method treats the distribution function (2.6) as a whole, instead of considering each moment as an individual variable. For the construction of an applicable scheme, it is necessary to provide a method which is able to apply addition or subtraction on the following two distributions:


Additionally, due to the cutoff, only ’s and ’s with are known. In [3], the authors proposed a homotopic method to calculate a new representation of :


where the coefficients ’s with can be worked out by solving the ODE system


for all until , and setting . In Eq. (2.12), is defined as


The system (2.12) will be solved by a Runge-Kutta type scheme. Once the form (2.11) is obtained, or can be calculated naturally.

By solving (2.12), we are able to represent for any and as


once there is one such representation known for some particular and . Here the ellipsis means the remaining coefficients are unknown. For any and and the associated representation of (2.14), we have


When and , (2.14) is called as the standard representation of . Note that the coefficients in (2.6) and (2.9) are in the sense of standard representation.

Another important technique is to calculate the flux based on the Hermite expansion of . Since only the moments with are known, can only be accurately given upto the -th order moments. Suppose is presented as (2.14) for some and , and we let


Then, by making use of the recurrence relation of the Hermite polynomials, we have


Here is taken to be zero if .

Once the linear operations between discrete distributions and the calculation of the flux are applicable, the convection progress can be simulated by a Riemann-solver-free finite volume method, such as Lax-Friedrichs scheme and HLL scheme (see e.g. [10]). Suppose the problem is one-dimensional, and the spatial mesh is uniform with mesh size . We denote by the distribution function on the -th grid at the time step . Other symbols such as and are defined similarly. Then, a first-order scheme with HLL numerical flux can be described as follows:

  1. Let and set to be the initial value, which is in its standard representation.

  2. Apply (2.9) to obtain the -st order moments:

  3. Solve the convection part with the HLL scheme:




    In Eq. (2.20), the signal speed and are approximated by


    where is the maximal root of Hermite polynomial . Eq. (2.21) is also used to determine the time step . We refer to [3] for details. Note that (2.19) is only accurate up to the th order moments.

  4. Solve the collision part analytically:

    1. For each , get its standard representation using (2.12) and (2.15). We suppose the result is

    2. Multiply all coefficients , by . The result is



  5. Increase by and return to step 2.

In step 3, since diffusion terms exist in the equations, the time step satisfies . The whole scheme is of first order.

3 A high resolution scheme

As is well known, the first-order HLL flux (2.20) adds excessive diffusion to the numerical solution in a general case. In order to reduce numerical diffusion, the technique of linear reconstruction is introduced below to the finite volume method. Suppose the boundary between the -th and -st cells is located at . Our aim is to construct two distributions and based on all cell averages to approximate the left and right limit of at point . In this section, the superscript will be omitted since the reconstruction is only applied on the same time step. Thus the numerical flux in (2.19) can be re-formulated as




During the reconstruction, different methods are applied to the convection part () and the diffusion part ().

3.1 Reconstruction for the convection part

For the convection part, it is important that the quantities used in reconstruction are conservative variables. Before reconstruction, according to the algorithm described in section 2.2, we already have the standard representations for all distributions , and the coefficients are assumed to be , . The simplest idea is to use together with and to make linear reconstruction:




where and are constants, and is a constant vector. However, this method leads to incorrect numerical results since none of the variables in Eq. (3.3) is conservative.

We take the original idea of the NR method, and consider the distribution function as a whole. Thus, linear reconstruction means taking the following approximation of and :


Here is a distribution. Obviously, the most convenient representation of is


Thus no ODEs are to be solved during the calculation of (3.5). Now, the coefficients ’s can be naturally given if we represent and as


In the implementation, we use the simplest minmod slope limiter for reconstruction


This reconstruction is a conservative reconstruction, since


where is a linear function of defined on as


However, this condition is not satisfied by the reconstruction (3.3).

3.2 Reconstruction for the diffusion part

The diffusion terms (2.9) provide approximation to the -st order moments. Since (2.9) is in the sense of standard representation, we first need to get the standard representations of and , and we use (3.4) to denote the results. The reconstruction of and with , which are involved in the ellipsis of Eq. (3.4), is a direct discretization of (2.9):


Now we give a general discussion on the order of accuracy. The time splitting introduces an error of magnitude . With (3.5), (3.6) and (3.8), the numerical flux (3.1) turns out to be a second order numerical flux. However, as discussed in [18, 20], due to the one-sided approximation of the diffusive gradients (3.11), the final accuracy only appears to be the first order. However, comparing with the original scheme described in Section 2.2, numerical error is significantly reduced by the linear reconstruction.

4 Enlarging the time step

The regularization of the moment method introduces diffusion terms into the system, which yields a relatively small time step . In order to enlarge the time step length, we use the RKC time stepping in the temporal discretization. In this section, a large time-stepping scheme with 2nd-order time integration will be proposed.

4.1 The RKC time-stepping

The RKC method is a series of explicit Runge-Kutta schemes for parabolic problems with large stability region and good internal stability. Our aim is to use in our algorithm without producing errors larger than the original method. Thus a second-order RKC scheme is needed. The -stage second-order RKC formula for the ODE system


was deduced in [21] as


where , , and


Here the parameters , , , are relevant to a manually selected damping factor . They are given by


where is the first kind Chebyshev polynomials


which can also be defined by the recurrence relation


The damping factor is often chosen as a small positive value such that a small imaginary perturbation of the eigenvalues of is allowable. In [8], the authors suggest that be chosen as , which results in a reduction in the stability boundary of about 2%.

For advection-diffusion problems, the imaginary part of the eigenvalues of may be large. Some analysis of the eigenvalue structure of upwind finite volume methods can be found in [19], where the authors show that the imaginary parts of the eigenvalues are less than for the second-order scheme if the time step satisfies the CFL condition of a pure advection problem. In order to ensure the stability, we follow [22] and use a large damping factor . Thus the stability boundary is approximately given by


In our implementation, the RKC method is only applied to the finite volume scheme (2.19). The time step is determined by


where is defined as


Then, we find the smallest positive integer satisfying


and use the -stage RKC method instead of the forward Euler scheme in (2.19). The inequalities (4.8) and (4.10) give and , which leads to the following estimation of the total computational time:


where is the length of the computational domain, and is the finishing time of the computation.

4.2 The Strang splitting

When the RKC time stepping is used, according to (4.8), the time step has the same magnitude as the grid size. Thus, when the gas is dense, we can still have only the first-order accuracy due to the convection-collision splitting, which introduces an error . In order to restore the time integration to a second-order one, the Strang splitting technique [14] is employed. In our algorithm, a direct usage of the Strang splitting can be described as follows:

  1. Let .

  2. Determine the time step .

  3. Solve the collision part over a half time step of length .

  4. Solve the convection part using the second-order RKC scheme over a time step of length .

  5. Solve the collision part again over a half time step of length .

  6. Increase by and return to Step 2.

This scheme requires twice calculation of the collision operator in one time step, but decrease the time integration error to . Here we introduce a method equivalent to the above Strang splitting scheme, but only once calculation of collision term is needed per time step.

According to (4.8) and (2.21), the time step determined in Step 2 is only relevant to and . Since the collision steps (Step 3 and Step 5) do not change these two quantities, we can actually calculate before Step 5. Thus the above algorithm can be rearranged as follows:

  1. Let and determine the time step .

  2. Solve the collision part over a half time step of length .

  3. Solve the convection part using the second-order RKC scheme over a time step of length .

  4. Determine the next time step .

  5. Solve the collision part over a half time step of length , and then solve the collision part again over a half time step of length .

  6. Increase by and return to Step 3.

Recall that for the BGK model, the collision part can be solved analytically as (2.23) and (2.24). Therefore, we can replace Step 5 by

  1. Solve the collision part over a time step of length .

Thus, the numerical result is identical to the original Strang splitting scheme, but only one collision step is performed in a complete time step.

5 Numerical examples

In this section, two numerical examples are presented to validate the efficiency and accuracy of our algorithm. In all the tests, the collision frequency is given by a simple form , which is corresponding to the Maxwell molecules. Here is the global Knudsen number, and it is slightly different from the Knudsen number defined in [2] (denoted by ) by


The CFL number is chosen to be . All the computations are performed on the Dell OptiPlex 755 desktop computer with a dual-core processor and CPU speed 2.33GHz.

5.1 An example with smooth solution

The first example is a repetition of the 1D periodic problem in [18, 3]. The computational domain is and the boundary condition is assumed to be periodic. The initial condition is


and the distribution is the local Maxwellian everywhere. The computation is stopped at . We set and . The numerical results for the density and temperature are plotted in Figure 1. Since the Knudsen number is large, the profiles of density and temperature differ a lot for different moment systems.

In order to examine the efficiency of our algorithm, we discretize the problem on a series of spatial grids with grid numbers ranging from to . Using , the analysis at the end of Section 4.1 gives


These results are validated by Figure 2. We can see that the large time step sizes are achieved, and the results are still stable. This leads to a remarkable reduction of the computational time.

Also, the first-order convergence rate is illustrated in Figure 3, where the “exact” solution is obtained on a mesh with grids and the errors are shown. We would like to remark that although the scheme is still of first order, the magnitude of error is much smaller than the original method. In [3], where the HLL scheme without reconstruction was utilized, the authors used grids to obtain numerical results with similar resolution as those in Figure 1. Note that small grid size implies smaller time steps, which leads to a huge computational cost.

(a) , density
(b) , temperature
(c) , density
(d) , temperature
(e) , density
(f) , temperature
Figure 1: Numerical solutions of Problem 5.1 using spatial grids
(a) Time cost
(b) Average step size
(c) Average effective step size
Figure 2: Computational costs and time step sizes for Problem 5.1 on different spatial grids. The -axis is the logarithm of the grid number .
Figure 3: Convergence rates for different moment systems

5.2 The shock-tube test

In this example, we show that our method is able to achieve high resolution when sharp layers exist in the numerical solution. Here a Riemann shock-tube problem is considered. It has been studied by Yang and Huang in [23] using the discrete ordinate method. The initial states are




The Knudsen number are selected by setting in (5.1). The computational domain is set to be , and we solve the problem until in order to match the results in [23].

Since is small, only the case is considered here. If large is used, the results are nearly identical to the current case. Some results are listed in the left column of Figure 4, whose validity can be confirmed by comparing them with those in [23]. In order to see the effects of reconstruction, we set in (3.8) and rerun the program. The results are in the right column of Figure 4. It is obvious that the left column provides much higher resolution near the shock wave, while the right column is even unable to achieve the correct peak value for and .

(b) without reconstruction
(d) , without reconstruction
(f) , without reconstruction
Figure 4: Density plots of the numerical solutions of Problem 5.2. is the grid number.

6 Concluding remarks

An efficient numerical scheme with high resolution for the NR method has been presented. Since the NR method gives a convective-diffusive system, we not only perform the linear reconstruction to gain a high spatial resolution, but also use the RKC schemes and the Strang splitting method to enlarge the time step while maintaining the order of accuracy in . In the future work, we are extending it into the 2D case with unstructured grids, together with the specularly reflective boundary conditions.


The research of the second author was supported in part by the National Basic Research Program of China under the grant 2010CBxxxxx and the National Science Foundation of China under the grant 10771008 and grant 10731060.


  • [1] P. L. Bhatnagar, E. P. Gross, and M. Krook. A model for collision processes in gases. I. small amplitude processes in charged and neutral one-component systems. Phys. Rev., 94(3):511–525, 1954.
  • [2] G. A. Bird. Molecular Gas Dynamics and the Direct Simulation of Gas Flows. Oxford: Clarendon Press, 1994.
  • [3] Z. Cai and R. Li. Numerical regularized moment method of arbitrary order for Boltzmann-BGK equation. SIAM J. Sci. Comput., 32(5):2875–2907, 2010.
  • [4] Z. Cai, R. Li, and Y. Wang. Numerical regularized moment method for high Mach number flow. arXiv:1101.5787, 2010.
  • [5] H. Grad. On the kinetic theory of rarefied gases. Comm. Pure Appl. Math., 2(4):331–407, 1949.
  • [6] H. Grad. The profile of a steady plane shock wave. Comm. Pure Appl. Math., 5(3):257–300, 1952.
  • [7] X. J. Gu and D. R. Emerson. A high-order moment approach for capturing non-equilibrium phenomena in the transition regime. J. Fluid Mech., 636:177–216, 2009.
  • [8] W. Hundsdorfer and J. Verwer. Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations. Springer, 2003.
  • [9] S. Jin and M. Slemrod. Regularization of the Burnett equations via relaxation. J. Stat. Phys, 103(5–6):1009–1033, 2001.
  • [10] R. J. Leveque. Finite Volume Methods for Hyperbolic Problems. Cambridge, 2002.
  • [11] S. Mizzi, X. J. Gu, D. R. Emerson, R. W. Barber, and J. M. Reese. Computational framework for the regularized 20-moment equations for non-equilibrium gas flows. Int. J. Num. Meth. Fluids, 56(8):1433–1439, 2008.
  • [12] I. Müller, D. Reitebuch, and W. Weiss. Extended thermodynamics – consistent in order of magnitude. Continuum Mech. Thermodyn., 15(2):113–146, 2002.
  • [13] I. Müller and T. Ruggeri. Rational Extended Thermodynamics, Second Edition, volume 37 of Springer tracts in natural philosophy. Springer-Verlag, New York, 1998.
  • [14] G. Strang. On the construction and comparison of difference schemes. SIAM J. Numer. Anal., 5(3):506–517, 1968.
  • [15] H. Struchtrup. Derivation of 13 moment equations for rarefied gas flow to second order accuracy for arbitrary interaction potentials. Multiscale Model. Simul., 3(1):221–243, 2005.
  • [16] H. Struchtrup. How many moments do we need, really? Presentation on Workshop on Moment Methods in Kinetic Theory, ETH Zürich, Switzerland, November 2008.
  • [17] H. Struchtrup and M. Torrilhon. Regularization of Grad’s 13 moment equations: Derivation and linear analysis. Phys. Fluids, 15(9):2668–2680, 2003.
  • [18] M. Torrilhon. Two dimensional bulk microflow simulations based on regularized Grad’s 13-moment equations. SIAM Multiscale Model. Simul., 5(3):695–728, 2006.
  • [19] M. Torrilhon and R. Jeltsch. Essentially optimal explicit Runge-Kutta methods with application to hyperbolic-parabolic equations. Numer. Math., 106(2):303–334, 2007.
  • [20] M. Torrilhon and K. Xu. Stability and consistency of kinetic upwinding for advection-diffusion equations. IMA J. Numer. Analy., 26(4):686–722, 2006.
  • [21] P. J. van der Houwen and B. P. Sommeijer. On the internal stability of explicit, -stage Runge-Kutta methods for large -values. Z. Angew. Math. Mech., 60(10):479–485, 1980.
  • [22] J. G. Verwer, B. P. Sommeijer, and W. Hundsdorfer. RKC time-stepping for advection-diffusion-reaction problems. J. Comput. Phys, 201(1):61–79, 2004.
  • [23] J. Y. Yang and J. C. Huang. Rarefied flow computations using nonlinear model Boltzmann equations. J. Comput. Phys., 120(2):323–339, 1995.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description