An Efficient NR Method for Boltzmann-BGK Equation
In , 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.
In 1949, the moment method was introduced by Grad  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 , the COET (Consistently Ordered Extended Thermodynamics) method , the order of magnitude approach , and et al. In , 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  was called as R by the authors of . In , 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  is similar to the original derivation of the R13 equations . And later in , the idea of COET  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  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  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  with some state-of-art 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  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.
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  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
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
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 , 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
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
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. ). 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:
Let and set to be the initial value, which is in its standard representation.
Apply (2.9) to obtain the -st order moments:
Solve the convection part with the HLL scheme:
In Eq. (2.20), the signal speed and are approximated by
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  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 , 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 , 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  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
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  is employed. In our algorithm, a direct usage of the Strang splitting can be described as follows:
Determine the time step .
Solve the collision part over a half time step of length .
Solve the convection part using the second-order 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 second-order 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.
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  (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
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 , 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 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  using the discrete ordinate method. The initial states are
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 . 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 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.
-  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.
-  G. A. Bird. Molecular Gas Dynamics and the Direct Simulation of Gas Flows. Oxford: Clarendon Press, 1994.
-  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.
-  Z. Cai, R. Li, and Y. Wang. Numerical regularized moment method for high Mach number flow. arXiv:1101.5787, 2010.
-  H. Grad. On the kinetic theory of rarefied gases. Comm. Pure Appl. Math., 2(4):331–407, 1949.
-  H. Grad. The profile of a steady plane shock wave. Comm. Pure Appl. Math., 5(3):257–300, 1952.
-  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.
-  W. Hundsdorfer and J. Verwer. Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations. Springer, 2003.
-  S. Jin and M. Slemrod. Regularization of the Burnett equations via relaxation. J. Stat. Phys, 103(5–6):1009–1033, 2001.
-  R. J. Leveque. Finite Volume Methods for Hyperbolic Problems. Cambridge, 2002.
-  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.
-  I. Müller, D. Reitebuch, and W. Weiss. Extended thermodynamics – consistent in order of magnitude. Continuum Mech. Thermodyn., 15(2):113–146, 2002.
-  I. Müller and T. Ruggeri. Rational Extended Thermodynamics, Second Edition, volume 37 of Springer tracts in natural philosophy. Springer-Verlag, New York, 1998.
-  G. Strang. On the construction and comparison of difference schemes. SIAM J. Numer. Anal., 5(3):506–517, 1968.
-  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.
-  H. Struchtrup. How many moments do we need, really? Presentation on Workshop on Moment Methods in Kinetic Theory, ETH Zürich, Switzerland, November 2008.
-  H. Struchtrup and M. Torrilhon. Regularization of Grad’s 13 moment equations: Derivation and linear analysis. Phys. Fluids, 15(9):2668–2680, 2003.
-  M. Torrilhon. Two dimensional bulk microflow simulations based on regularized Grad’s 13-moment equations. SIAM Multiscale Model. Simul., 5(3):695–728, 2006.
-  M. Torrilhon and R. Jeltsch. Essentially optimal explicit Runge-Kutta methods with application to hyperbolic-parabolic equations. Numer. Math., 106(2):303–334, 2007.
-  M. Torrilhon and K. Xu. Stability and consistency of kinetic upwinding for advection-diffusion equations. IMA J. Numer. Analy., 26(4):686–722, 2006.
-  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.
-  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.
-  J. Y. Yang and J. C. Huang. Rarefied flow computations using nonlinear model Boltzmann equations. J. Comput. Phys., 120(2):323–339, 1995.