1 Introduction
Rarefied gas flow simulation is always a research hotspot of computational fluid dynamics (CFD). Recent years, multiscale gaskinetic 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 gaskinetic scheme (UGKS) [6] by Xu and Huang, the discrete unified gaskinetic 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 NavierStokes (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
[9], Zhu et al.’s prediction based implicit UGKS [10, 11], Zhu et al.’s implicit multigrid UGKS algorithm [12], Yang et al.’s memory saving implicit multiscale scheme [13], Pan et al.’s implicit DUGKS [14], 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. [10] 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 BGKtype equation [15],
(1) 
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,
(2) 
or if the Shakhov model [16] is used
(3) 
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
(4) 
where
is the vector of the macroscopic conservative variables,
is the vector of moments
,is the velocity space element. The stress tensor
and the heat flux can also be calculated by as(5) 
(6) 
Moreover, and obey the conservation law,
(7) 
Adopting the integral form, the steady state of the governing equation Eq. 1 is
(8) 
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
(9) 
where the flux has the relation with the distribution function by
(10) 
This paper is about the numerical method of determining the steady state defined by Eq. 8. It is timeconsuming 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. [10], and has a structure similar to the multigrid method of Zhu et al. [12], 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
(11) 
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
(12) 
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 [17] is adopted, i.e.
(13) 
where
(14) 
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 [18] 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
(15) 
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
(16) 
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
(17) 
The physical local time step for the cell is determined by the local CFL condition as
(18) 
where is the Heaviside function defined as
(19) 
For more details about the construction of the interface distribution function please refer to reference [17].
Eq. 12 is solved by iterations. The microscopic residual at the th iteration can be defined as
(20) 
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,
(21) 
where is the pseudo time step and is always set to be in the present study. Combined with the residual expression Eq. 20, Eq. 21 can be written as
(22) 
For the increment of the interface distribution function , it is simply handled by a modified upwind scheme
(23) 
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.
(24) 
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
(25)  
where is the set of ’s neighboring cells satisfying while for it satisfies . For simplicity, Eq. 25 is solved by the Symmetric GaussSeidel (SGS) method, or also known as the Point Relaxation Symmetric GaussSeidel (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 matrixfree 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
(26) 
and the macroscopic residual defined by the macroscopic governing equation Eq. 9 at the th step can be calculated from the flux by
(27) 
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
(28) 
where the term is the integral error compensation term to make the scheme conservative, more details about this term please refer to reference [17].
The iteration of the loop is similar to the numerical smoothing process in multigrid method [12]. 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 variable
is 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 [12].2.2.1 Framework
The macroscopic residual has the form of Eq. 27. To reduce the residual, a prediction equation is constructed by backward Euler formula
(29) 
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
(30) 
where and are fluxes calculated by the prediction solver from and the predicted with data reconstruction. This prediction solver is welldesigned 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
(31)  
and the corresponding increment equation for is
(32) 
where is the pseudo time step. Considering Eq. 31, the increment of the residual can be expressed as
(33) 
the variation of the flux is further approximated by
(34) 
where has the form [21] of the wellknown Roe’s flux function
(35) 
Here is the Euler flux
(36) 
and is
(37) 
where is the acoustic speed at the interface and is the distance between cell center and . Substitute Eq. 33, Eq. 34 and Eq. 35 into Eq. 32, approximate by , and note that holds, we can get
(38)  
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 ChapmanEnskog expansions [22] 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 [12]. 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 ChapmanEnskog expansion [22], the distribution function obtained from the BGK equation Eq. 1 to the first order of is
(39) 
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 gaskinetic scheme (GKS) or also known as BGKNS scheme [18], we replace it by a weighting of Euler flux and the flux of kinetic flux vector splitting (KFVS) [24]. That is, we replace the term in Eq. 39 and the interface distribution function is expressed as
(40) 
where is
(41) 
which is determined by the reconstructed macroscopic variables on the two side of the interface. The interface macroscopic variable is calculated as
(42) 
and is obtained from . The weight factors and share the same forms as those in Eq. 13 (for how these weight factors are constructed please refer to [17]), and is calculated by
(43) 
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 [18], 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 NStype 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
(44) 
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
(45) 
which has and . is related to the viscous term and calculated as
(46) 
where and correspond to the heat flux and stress in NS equation, is the specific heat at constant pressure. is a molecular model coefficient [26] involved in the variable soft sphere (VSS) model [27, 28] and is calculated as
(47) 
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 onedimensional case where there is no stress but only heat flux, i.e. and . In this case is
(48) 
and the heat flux from the viscous term of Eq. 44 is