Rarefied gas flow simulation is always a research hotspot of computational fluid dynamics (CFD). Recent years, multiscale gas-kinetic methods based on the discrete velocity method (DVM, [1, 2, 3, 4, 5]) framework for nonequilibrium rarefied flow simulation have been developed, like the unified gas-kinetic scheme (UGKS)  by Xu and Huang, the discrete unified gas-kinetic scheme (DUGKS) [7, 8]
by Guo et al. These multiscale methods overcome the time step and cell size restrictions of the original DVM method which requires time step and cell size of the order of mean collision time and mean free path, and thus have attracted more and more researchers’ attention. It is worth pointing out that although UGKS and DUGKS can adopt time step and cell size comparable to the traditional macroscopic Navier-Stokes (NS) method, they still involve large amount of computation due to the curse of dimensionality. Hence, many researches on the acceleration of these multiscale methods have been carried out, including Mao et al.’s implicit UGKS, Zhu et al.’s prediction based implicit UGKS [10, 11], Zhu et al.’s implicit multigrid UGKS algorithm , Yang et al.’s memory saving implicit multiscale scheme , Pan et al.’s implicit DUGKS , etc. Following these previous works, it is quite valuable to further develop the fast algorithm for the multiscale method.
In this paper, a multiple prediction implicit multiscale method for steady state calculation of gas flow in all flow regimes is proposed. The idea of macroscopic prediction presented by Zhu et al.  is further developed. A prediction solver is used to predict the macroscopic variable based on the macroscopic residual, and a multiple prediction procedure is constructed. The prediction solver is designed to ensure the accuracy of the predicted macroscopic variable in the continuum flow regime and the stability of the numerical system in all flow regimes, which makes the method very efficient in the continuum flow regime and stable in all flow regimes. Our test cases show that the present method is one order of magnitude faster than the previous implicit multiscale method in the continuum flow regime.
2 Numerical method
In this paper, the monatomic gas is considered and the governing equation is BGK-type equation ,
where is the gas particle velocity distribution function, is the particle velocity, is the relaxation time calculated as ( and are the viscosity and pressure). is the equilibrium state which has a form of Maxwellian distribution,
or if the Shakhov model  is used
where is the peculiar velocity and is the macroscopic gas velocity, is the heat flux, is a variable related to the temperature by . Pr is the Prandtl number and has a value of for monatomic gas. is related to the macroscopic variables by
is the vector of the macroscopic conservative variables,
is the vector of moments,
is the velocity space element. The stress tensorand the heat flux can also be calculated by as
Moreover, and obey the conservation law,
Adopting the integral form, the steady state of the governing equation Eq. 1 is
where is the control volume, is the volume element, is the surface area element and is the outward normal unit vector. Take the moment of Eq. 8 for , the corresponding macroscopic governing equation can be written as
where the flux has the relation with the distribution function by
This paper is about the numerical method of determining the steady state defined by Eq. 8. It is time-consuming to directly solve Eq. 8 through a microscopic scheme involving discretization in both physical space and velocity space. The main idea of the present method is summarized as that, using the accurate but expensive scheme to calculate the residual of the system deviating from the steady state, then utilizing this residual and using the less accurate but efficient scheme to do the evolution. More distinctly, an accurate multiscale microscopic scheme based on the DVM framework is used to handle the microscopic numerical system with Eq. 8, and a fast prediction scheme is used to do the evolution of the macroscopic variables. The prediction scheme can be some kind of macroscopic scheme based on macroscopic variables or even a scheme based on the DVM framework but with less velocity points. The schematic of the general algorithmic framework for the present method is shown in Fig. 1. The method consists of several loops in different layers. The outermost loop is denoted by . One iteration of the loop includes a loop denoted by and a loop denoted by , where the macroscopic variable and the residual are given as the input, the new and are the output. In the loop, the predicted macroscopic variable is determined by the prediction scheme and by the numerical smoothing process. In the loop, the microscopic variable is calculated and the new and are obtained. The present method is a development of the prediction method of Zhu et al. , and has a structure similar to the multigrid method of Zhu et al. , therefore we call it as “multiple prediction method”. The method is detailed in following paragraphs.
2.1 Construction of the loop
In the loop, residuals of the numerical system are evaluated through the microscopic scheme and the microscopic variables (the discrete distribution function) are updated through an implicit method (the numerical smoothing process). The microscopic scheme is very important because it determines the final steady state of the whole numerical system and thus determines the nature of the present numerical method.
The microscopic scheme is based on Eq. 8. Discretizing the physical space by finite volume method and discretizing the velocity space into discrete velocity points, the microscopic governing equation Eq. 8 can be expressed as
where the signs correspond to the discretizations in physical space and velocity space respectively. denotes the neighboring cell of cell and is the set of all of the neighbors of . Subscript denotes the variable at the interface between cell and . is the interface area, is the outward normal unit vector of interface relative to cell , and is the volume of cell . The loop aims to find the solution of Eq. 11 with the input predicted variable , therefore Eq. 11 can be written more exactly as
where the symbol denotes the predicted variables at the th step. and can be directly calculated from the input variable . The distribution function at the interface is very important to ensure the multiscale property of the scheme. In this paper, following the idea of DUGKS [7, 8], the construction of in reference  is adopted, i.e.
In above equations, and can be obtained through the reconstruction of the distribution function data. and are calculated by the same way as the method of GKS  and they can be both calculated from the predicted macroscopic variable . For , it is determined by the interface macroscopic variables , which can be calculated as
where the superscripts and denote variables at the left and right sides of the interface, and can be determined after the spacial reconstruction of . For , it is calculated as
where the pressure , at two sides of the interface can be obtained from the reconstruction and the second term on the right is for artificial viscosity. in above equations is calculated from the physical local time step
The physical local time step for the cell is determined by the local CFL condition as
where is the Heaviside function defined as
For more details about the construction of the interface distribution function please refer to reference .
Eq. 12 is solved by iterations. The microscopic residual at the th iteration can be defined as
According to the previous descriptions, can be calculated from and through the spatial data reconstruction. The increment equation to get the microscopic variable at the iteration is constructed by backward Euler method,
For the increment of the interface distribution function , it is simply handled by a modified upwind scheme
where the coefficient is the corresponding coefficient multiplied by in Eq. 13. This coefficient is multiplied because during the whole loop the term in Eq. 13 is calculated by the predicted macroscopic variable and therefore is an invariant, so the variation of the microscopic variable only influences the term , which is multiplied by the coefficient . It is worth noting that in an actual implementation of the method, the interface distribution function at the th step may be taken as the initial value at for the step to reduce computation cost, in this situation the variation should also account for the variation of , and the coefficient shouldn’t be multiplied at the first iteration of the loop, i.e.
In this situation, after the first iteration of the loop, the interface distribution function will be calculated with the newly predicted and Eq. 23 is used to handle again. Without loss of generality, substituting Eq. 23 into Eq. 22 will yield
where is the set of ’s neighboring cells satisfying while for it satisfies . For simplicity, Eq. 25 is solved by the Symmetric Gauss-Seidel (SGS) method, or also known as the Point Relaxation Symmetric Gauss-Seidel (PRSGS) method [19, 20]. In each time of the SGS iteration, a forward sweep from the first to the last cell and a backward sweep from the last to the first cell are implemented, during which the data of a cell is always updated by the latest data of its adjacent cells through Eq. 25. Such a SGS iteration procedure is totally matrix-free and easy to implement.
After several times of SGS iterations for solving Eq. 25, an evaluation of with a certain precision can be obtained. Then the residual at the th iteration of the loop can be computed from and , and a new turn of the loop will be performed. After several iterations of the loop, an evaluation of with a certain precision can be obtained, and the interface distribution function can be calculated by Eq. 13. Then the macroscopic numerical flux at the interface can be got by numerical integral in the discrete velocity space
and the macroscopic residual defined by the macroscopic governing equation Eq. 9 at the th step can be calculated from the flux by
Note that in the loop we solve the microscopic system Eq. 12 which is under the condition of the predicted variable , so is not zero even if the microscopic system is solved sufficiently accurately. Finally, the macroscopic variable is calculated by numerical integral as
where the term is the integral error compensation term to make the scheme conservative, more details about this term please refer to reference .
The iteration of the loop is similar to the numerical smoothing process in multigrid method . The computation procedure of the loop is listed as follows:
- Step 1.
Set the initial value .
- Step 2.
- Step 3.
Make judgement: if the residual meets the convergence criterion, or if the iteration number of the loop meets the maximum limit, break out of the loop and go to Step 5.
- Step 4.
Solve Eq. 25 by several times of SGS iterations, obtain , and go to Step 2.
- Step 5.
2.2 Construction of the loop
In the loop, based on the macroscopic variable and the macroscopic residual at the
th step, a reasonable estimation for the macroscopic variableis obtained through a fast prediction scheme to accelerate convergence. Theoretically speaking, the prediction scheme can be either a macroscopic scheme based on macroscopic variables or a microscopic scheme based on the DVM framework but with less velocity points. In this paper, a macroscopic scheme is designed to do the prediction. The process of the loop has certain similarity to the coarse grid correction in the multigrid method .
The macroscopic residual has the form of Eq. 27. To reduce the residual, a prediction equation is constructed by backward Euler formula
is the local prediction time step, the purpose of this time step is to constrain the marching time depth of the prediction process to make the scheme stable in the extreme case. The predicted residual increment is calculated by
where and are fluxes calculated by the prediction solver from and the predicted with data reconstruction. This prediction solver is well-designed to balance between accuracy and stability, and will be presented later in the next section.
The aim of the loop is to solve Eq. 29 and give an estimation for with a certain precision. Like what we do in the loop, Eq. 29 is also solved by iterations. The residual at the th iteration can be defined by Eq. 29 and expressed as
and the corresponding increment equation for is
where is the pseudo time step. Considering Eq. 31, the increment of the residual can be expressed as
the variation of the flux is further approximated by
where has the form  of the well-known Roe’s flux function
Here is the Euler flux
Eq. 38 is solved by several times’ SGS iterations. An estimation of with a certain precision can be obtained from Eq. 38, then and the residual at the th iteration of the loop can be calculated. After several turns of the loop, the predicted macroscopic variable with a certain precision can be determined.
In fact, utilizing and , a prediction for the microscopic variable can also be obtained to accelerate the convergence of the microscopic numerical system (i.e. the loop). The increment of the distribution function can be calculated from the Chapman-Enskog expansions  based on macroscopic variables and . This strategy will increase the complexity of the algorithm and thus is not adopted in the present method.
Likewise, as one can see, the process of the loop is similar to the numerical smoothing process in multigrid method . The computation procedure of the loop is listed as follows:
- Step 1.
Set the initial value .
- Step 2.
Calculate the residual by Eq. 31 from , and (data reconstruction is implemented).
- Step 3.
Make judgement: if the residual meets the convergence criterion, or if the iteration number of the loop meets the maximum limit, break out of the loop and the predicted macroscopic variable is determined.
- Step 4.
Solve Eq. 38 by several times of SGS iterations, obtain , and go to Step 2.
2.2.2 Prediction solver
The prediction solver used to calculate the fluxes and in Eq. 30 requires careful design. For the continuum flow, the prediction solver should be as accurate as a traditional NS solver. For the rarefied flow, it’s unrealistic for a fast solver based on macroscopic variables to provide a very precise flux, but the solver should be stable so that the present method can be applied to all flow regimes. Thus, there are two principles for the prediction solver: accurate in the continuum flow regime, stable in all flow regimes.
We start constructing the solver from the view of gas kinetic theory. Based on the famous Chapman-Enskog expansion , the distribution function obtained from the BGK equation Eq. 1 to the first order of is
Suppose there is an interface in direction. If the interface distribution function has the form of Eq. 39, take moments of to Eq. 39 and ignore second (and higher) order terms of , we can get the NS flux [18, 23], where the term corresponds to the Euler flux and terms with (i.e. terms except ) correspond to viscous terms in the NS flux. Flux directly calculated from Eq. 39 will lead to divergence in many cases, and we introduce some modifications below.
The Euler flux often causes stability issue. Inspiring by gas-kinetic scheme (GKS) or also known as BGK-NS scheme , we replace it by a weighting of Euler flux and the flux of kinetic flux vector splitting (KFVS) . That is, we replace the term in Eq. 39 and the interface distribution function is expressed as
which is determined by the reconstructed macroscopic variables on the two side of the interface. The interface macroscopic variable is calculated as
where is for artificial viscosity. is the local CFL time step and is equal to in Eq. 13. Eq. 40 has a form similar to the interface distribution function of GKS , except that the viscous term is not upwind split and the weight factor is constructed following the thought of DUGKS [7, 8]. Because the KFVS scheme is very robust, the flux obtained from Eq. 40 makes the numerical system more stable than directly using Eq. 39. In the continuum flow regime, , if the flow is continuous the term for artificial viscosity will be negligible and Eq. 40 will recover the NS flux, while if the flow is discontinuous the term will be activated and Eq. 40 will work as a stable KFVS solver. In the rarefied flow simulation, and the inviscid part of Eq. 40 generally provides a KFVS flux, which can increase the stability of the scheme.
The flux obtained from Eq. 40 works well in the continuum flow regime. However, in the case of large Kn number, the numerical system based on Eq. 40 is very stiff due to the large NS-type linear viscous term and the scheme is easy to blow up. Therefore, we multiply the viscous term by a limiting factor and Eq. 40 is transformed into
Here we emphasize that the limiting factor aims not to accurately calculate the flux, but to increase the stability in the case of large Kn number. One can view it as an empirical parameter. The limiting factor is constructed considering the form of nonlinear coupled constitutive relations (NCCR) [25, 26], and is expressed as
which has and . is related to the viscous term and calculated as
where and correspond to the heat flux and stress in NS equation, is the specific heat at constant pressure. is a molecular model coefficient  involved in the variable soft sphere (VSS) model [27, 28] and is calculated as
where the molecular scattering factor and the heat index depend on the type of gas molecule. The limiting factor is constructed to weaken the viscous term in large Kn number case to make the scheme stable. It can be seen from Eq. 45 and Eq. 46 that, when the stress and heat flux are small, is approaching to and we can get the NS viscous term in Eq. 44, when the stress and heat flux are large, is approaching to and the viscous term is weakened. Here we further reveal the mechanism of through a simple one-dimensional case where there is no stress but only heat flux, i.e. and . In this case is
and the heat flux from the viscous term of Eq. 44 is