An Efficient NR Method for BoltzmannBGK Equation
Abstract
In [3], we proposed a numerical regularized moment method of arbitrary order (abbreviated as NR method) for BoltzmannBGK 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 13moment system, but the system is problematic: the hyperbolicity can only be obtained in the neighbourhood of Maxwellian, and the structure of shock wave is nonsmooth 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 JinSlemrod 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 ChapmanEnskog 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 BoltzmannBGK 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 stateofart techniques available to us. Precisely, three numerical techniques are employed:

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.

RKC time integration scheme. The RKC time integration scheme [21] is adopted to enlarge the time step sizes. The RKC method is a RungeKutta type method originally designed for diffusive PDEs to provide large stability regions, while the regularized moment equations are convectivediffusive problems, where imaginary parts appear in the eigenvalues of the right hand side of the semidiscrete 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.

Revised Strang splitting method. With enlarged time step, the numerical accuracy can be harmed by convectioncollision 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 BoltzmannBGK 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 BoltzmannBGK 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 BoltzmannBGK equation reads
(2.1) 
where is the distribution function, is the collision frequency, and is the local Maxwellian defined by
(2.2) 
Here , and are local macroscopic variables which represent the density, velocity and temperature respectively. They can be calculated from the distribution function by
(2.3)  
(2.4)  
(2.5) 
2.2 A firstorder 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:
(2.6) 
where is a dimensional multiindex. ’s and ’s act as the basis functions and the corresponding coefficients respectively, and is defined by
(2.7) 
where ’s are the orthogonal Hermite polynomials defined by
(2.8) 
This idea originates from [5] where the 13moment equations are derived. In order to get a finite system, the NR method [4] chooses an positive integer and approximate with by
(2.9) 
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:
(2.10) 
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 :
(2.11) 
where the coefficients ’s with can be worked out by solving the ODE system
(2.12) 
for all until , and setting . In Eq. (2.12), is defined as
(2.13) 
The system (2.12) will be solved by a RungeKutta 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
(2.14) 
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
(2.15)  
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
(2.16) 
Then, by making use of the recurrence relation of the Hermite polynomials, we have
(2.17) 
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 Riemannsolverfree finite volume method, such as LaxFriedrichs scheme and HLL scheme (see e.g. [10]). Suppose the problem is onedimensional, 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 firstorder scheme with HLL numerical flux can be described as follows:

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

Apply (2.9) to obtain the st order moments:
(2.18) 
Solve the convection part with the HLL scheme:
(2.19) where
(2.20) In Eq. (2.20), the signal speed and are approximated by
(2.21) 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.

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 firstorder 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 reformulated as
(3.1) 
where
(3.2)  
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:
(3.3) 
and
(3.4) 
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 :
(3.5) 
Here is a distribution. Obviously, the most convenient representation of is
(3.6) 
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
(3.7) 
In the implementation, we use the simplest minmod slope limiter for reconstruction
(3.8) 
This reconstruction is a conservative reconstruction, since
(3.9) 
where is a linear function of defined on as
(3.10) 
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):
(3.11) 
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 onesided 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 timestepping scheme with 2ndorder time integration will be proposed.
4.1 The RKC timestepping
The RKC method is a series of explicit RungeKutta 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 secondorder RKC scheme is needed. The stage secondorder RKC formula for the ODE system
(4.1) 
was deduced in [21] as
(4.2)  
where , , and
(4.3) 
Here the parameters , , , are relevant to a manually selected damping factor . They are given by
(4.4)  
where is the first kind Chebyshev polynomials
(4.5) 
which can also be defined by the recurrence relation
(4.6) 
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 advectiondiffusion 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 secondorder 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
(4.7) 
In our implementation, the RKC method is only applied to the finite volume scheme (2.19). The time step is determined by
(4.8) 
where is defined as
(4.9) 
Then, we find the smallest positive integer satisfying
(4.10) 
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:
(4.11) 
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 firstorder accuracy due to the convectioncollision splitting, which introduces an error . In order to restore the time integration to a secondorder one, the Strang splitting technique [14] is employed. In our algorithm, a direct usage of the Strang splitting can be described as follows:

Let .

Determine the time step .

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

Solve the convection part using the secondorder RKC scheme over a time step of length .

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

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:

Let and determine the time step .

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

Solve the convection part using the secondorder RKC scheme over a time step of length .

Determine the next time step .

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 .

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

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
(5.1) 
The CFL number is chosen to be . All the computations are performed on the Dell OptiPlex 755 desktop computer with a dualcore 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
(5.2) 
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
(5.3) 
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 firstorder 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.
5.2 The shocktube 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 shocktube problem is considered. It has been studied by Yang and Huang in [23] using the discrete ordinate method. The initial states are
(5.4) 
with
(5.5) 
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 .
6 Concluding remarks
An efficient numerical scheme with high resolution for the NR method has been presented. Since the NR method gives a convectivediffusive 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.
Acknowledgements
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.
References
 [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 onecomponent 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 BoltzmannBGK 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 highorder moment approach for capturing nonequilibrium phenomena in the transition regime. J. Fluid Mech., 636:177–216, 2009.
 [8] W. Hundsdorfer and J. Verwer. Numerical Solution of TimeDependent AdvectionDiffusionReaction 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 20moment equations for nonequilibrium 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. SpringerVerlag, 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 13moment equations. SIAM Multiscale Model. Simul., 5(3):695–728, 2006.
 [19] M. Torrilhon and R. Jeltsch. Essentially optimal explicit RungeKutta methods with application to hyperbolicparabolic equations. Numer. Math., 106(2):303–334, 2007.
 [20] M. Torrilhon and K. Xu. Stability and consistency of kinetic upwinding for advectiondiffusion 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 RungeKutta methods for large values. Z. Angew. Math. Mech., 60(10):479–485, 1980.
 [22] J. G. Verwer, B. P. Sommeijer, and W. Hundsdorfer. RKC timestepping for advectiondiffusionreaction 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.