An Implicit Unified Gaskinetic Scheme for Unsteady Flow in All Knudsen Regimes
Abstract
The unified gaskinetic scheme (UGKS) is a direct modeling method for multiple scale transports. The modeling scale of the scheme is the mesh size and time step, and the ratios of mesh size over particle mean free path or the time step over the particle collision time determine the local evolution regime of flow physics. For a multiscale flow problem, such as the hypersonic flow around a flying vehicle in near space, the UGKS is able to capture the highly compressed NavierStokes solution in one region and fully expanded free molecular flow in another region, with the significant variations of the ratio of the time step over the local particle collision time around the vehicle. For an explicit UGKS, the time step in the whole computational domain is determined by the CFL condition. For steady state computation, in order to increase the efficiency of the scheme, the implicit and multigrid techniques have been implemented in UGKS zhu2016implicit (); zhu2017multigrid (), where the efficiency can be improved by two orders of magnitude in comparison with the explicit scheme. However, for unsteady flow computation, due to the CFL condition the time step in the explicit UGKS is limited by the smallest cell size in the computational domain. As a result, for a largely stretched nonuniform mesh the global time step becomes very small and the ratio of the time step over the local particle collision time may get a very small value. Under such a circumstance, even though the UGKS is a multiscale method, the physics in explicit UGKS may be constrained by the kinetic scale physics only, such as the overall time step being less than the local particle collision time. Therefore, for unsteady flow the multiscale UGKS may become a single scale discrete velocity method (DVM) where the free transport is used for the flux evaluation. In order to keep the advantages of UGKS and improve the efficiency for unsteady flow computation, the restriction from global CFL condition on the time step has to be removed. In this paper, we will develop an implicit UGKS (IUGKS) for unsteady flow by alternatively solving the macroscopic and microscopic governing equations within a step iteratively. With a predefined uniform large evolution time step, the local CFL number varies greatly in different region, such as on the order in the large numerical cell size region, and in the small cell size region. In order to preserve the coherent flow evolution and multiscale nature, the numerical flux across a cell interface is still evaluated by the explicit UGKS with the local CFL condition. Therefore, the multiscale property of the UGKS modeling is preserved over nonuniform meshes. With improved temporal discretization, the current IUGKS can automatically go back to the explicit UGKS and obtain identical solutions when the time step of the implicit scheme gets to that of an explicit scheme. Many numerical examples are included to validate the scheme for both continuum and rarefied flows with a large variation of artificially generated mesh size. The IUGKS has a second order accuracy and presents reasonably good results for unsteady flow computation and its efficiency has been improved by dozen of times in comparison with the explicit UGKS.
keywords:
unified gaskinetic scheme, implicit method, unsteady flow, multiscale transport1 Introduction
The unified gaskinetic scheme (UGKS) is a direct modeling method on the scale of mesh size and time step xu2010unified (); huang2012unified (); xu2015book (), which is able to capture multiscale flow physics in all Knudsen number regimes. In UGKS, a time evolution solution based on the kinetic relaxation model is used at a cell interface for the flux evaluation, where the ratio of the time step over the particle collision time determines the local flow regime. Different from the gaskinetic scheme (GKS) xu2001gas () for the update of macroscopic flow varibales only for the continuum flow, the UGKS updates both the macroscopic flow variables and the gas distribution function. Therefore, highly nonequilibrium effects can be properly captured by UGKS through the update of the gas distribution function, which has much more degrees of freedom to describe the nonequilibrium property. Due to its unified treatment for the continuum and rarefied, as well as high and low speed flows, the UGKS shows excellent performance in many flow problems, such as the studies of hypersonic flows chen2012adaptive (); liu2014diatomic (), rarefied flows jiang2016thesis (), and microflows huang2013unified (); liu2015micro (). In addition, the methodology of direct modeling xu2015book (); xu2017paradigm () can be extended to other multiscale transport processes, such as plasma liu2017plasma (), radiative transfer sun2015radiative (), and phonon heat transfer guo2016discrete (). The UGKS provides a useful tool to solve complex multiscale problems, and has great potentials in the engineering applications, such as microelectromechanical system and spacecraft design.
One of distinguishable features of the UGKS is to use the scaledependent flow evolution in the numerical flux evaluation, where the ratio of time step over particle mean collision time () identifies the continuum ()and rarefied () flow regimes and this ratio can be changed by several orders of magnitude continuously over the computational domain for a multiscale flow problem. For a single scale method, such as the direct simulation of Monte Carlo (DSMC) method and the discrete velocity method (DVM), it always requires to resolve the flow physics on the scale of (), where the mesh size and time step are less than the mean free path and mean collision time. These constraints on the cell size and time step are not enforced by UGKS because the time evolution solution with coupled collision and transport is used in the interface flux evaluation instead of free transport (upwind) only for the single scale method chenap (). The truly multiscale nature makes the UGKS very efficient in the near continuum flow regime. However, for unsteady flow simulations the explicit UGKS uses a uniform time step which is determined by the CourantFriedrichsLewy (CFL) condition and the smallest cell size may determine the time step. For hypersonic flow with a wide range of cell size, such as the cell size around the corner or the tip of a wing of a flying vehicle, the overall time step can be very small. So, for the explicit scheme the global time step may go to the order of collision time . As a result, the advantage of UGKS may be demolished when using a time step everywhere with the order , even though in most of the computational domain a much large cell size is used and the flow is in the near equilibrium one, such as those far away from the vehicle surface. The above permissible uniform CFL time step for stability is much smaller than the preferred time step for a reasonable accuracy for unsteady flow computation jameson1991time (). Therefore, it is necessary to develop an implicit UGKS (IUGKS) that can break the barrier from the CFL stability condition from these very small cell sizes, and preserve the multiscale nature of the UGKS for simulating different flow regime physics in different locations at the same time.
For unsteady flow simulation in continuum regime, implicit treatments to solve the NavierStokes (NS) equations have already been fully developed using dualtime stepping methods jameson1991time (); pulliam1993time (); tomaro1997implicit (). Recently, the timeaccurate implicit GKS has also been constructed tan2017time (); li2017implementation () to accelerate computational efficiency for highspeed unsteady and turbulent flows. In rarefied flow study, there are several implicit schemes that have been constructed based on the DVM yang1995rarefied (); mieussens2000discrete (); peng2016implicit () and the UGKS mao2015implicit (); zhu2016implicit () for steady state solutions only. In the study of radiative transfer sun2017implicit (), implicit techniques have been used in UGKS to remove the restriction of speed of light on the numerical time step for unsteady radiative transport. In order to fully utilize the multiscale nature of UGKS for unsteady flow simulations, the construction of an implicit UGKS is necessary, especially for mutiscale flow problem with both continuum and rarefied flow regimes and a highly stretched mesh distributions.
The current paper presents an implicit UGKS (IUGKS) for unsteady flows in all Knudsen number regimes. Starting from discretized conservation laws for conservative flow variables and the gas distribution function, we employ the semidiscrete formulation to construct the implicit scheme. To maintain the multiscale transport property, the time derivatives in the semidiscrete UGKS for the flow variables are evaluated from the averaged transport according to the local time step. The temporal derivatives involve the cell size effect on the flow evolution with respect to the local particle mean free path. Based on the implicit discretization for the macroscopic flow variables and the gas distribution function, a matrix connecting the variables in the whole computational domain can be easily solved through the lowerupper symmetric GaussSeidel (LUSGS) method jameson1987lower (); yoon1988lower (), pointrelaxation (PR) method yuan2002comparison (), and the generalized minimal residual (GMRES) algorithm saad1986gmres (). During the inner iterative process within each timemarching step, the implicit governing equations for macroscopic flow variables and microscopic distribution are solved alternatively with high efficiency for both continuum and rarefied flows.
2 Implicit unified gaskinetic scheme
2.1 Unified gaskinetic scheme
The UGKS is constructed based on a direct modeling of the flow evolution on the scale of time step and cell size. For a discrete finite volume cell and discretized time step , the discretized governing equation for the gas distribution function is
(1) 
where denotes the set of neighboring cells of , and cell is one of the neighbors. The interface connecting the cells and is referred to as . is volume of cell , and is area of the cell interface . Here is the normal component of the particles’ microscopic velocity , which is perpendicular to the cell interface. is a timedependent distribution function at the cell interface , which describes the time evolution of flow physics around the cell interface. The collision term takes into account the particles’ interaction in one time step. Eq. (1) is an exact physical law for the discrete cell , which simply describes the evolution of the gas distribution function due to particles’ transport and collision on the discretized level. The corresponding macroscopic governing equations are
(2) 
where for two dimensional cases. and are the particles’ velocities along  and direction, respectively. denotes the particles’ random motion for the internal degrees of freedom. is a vector of the conservative flow variables, i.e., representing the density of mass, momentum and energy. Eq. (2) describes the conservation of macroscopic variables, which is exact as well.
The above two governing equations are the physical laws on a discretized level. The accuracy of the solutions from Eq.(1) and (2) depends on the construction of the collision term and the flux transport, which is similar to the constitutive relationship provided in the NS equations. The special ingredient of the UGKS is to construct the collision term and the flux function according to the scale of cell size and time step. According to the scale of time step relative to the particle collision time , the full Boltzmann collision term can be replaced by the BGKtype kinetic model in the Eq.(1), especially for the cases of liu2016boltzmann (). In the following, the kinetic equation with the BGKtype model will be used only. The multiscale nature of the UGKS comes from the construction of the timedependent distribution function at the cell interface, , which is the evolution solution of the BGK model. The analytic solution gives
(3)  
from the reconstructed initial gas distribution function
(4) 
and an equilibrium state distributed in time and space
(5) 
where and are the spatial derivatives of the initial distribution functions for the left and right hand sides of the interface. and are the spatial and temporal derivatives of the equilibrium state, respectively. is the Maxwellian distribution, corresponding to the local conservative flow variables.
From the first two terms in Eq. (3), we can observe a transition from the initial nonequilibrium distribution function to the equilibrium one with the increment of time scale, which denotes a flow transition from kinetic scale to hydrodynamic one. It describes an accumulating effect from particle collision in the evolution process. The physical picture is illustrated in Fig. 1, where a time evolution of flow physics around the interface within a cellsizedetermined time step is presented. The automatic scale variation in the physical solution depends on the ratio of , which is the key ingredient to construct a multiscale scheme. Details about the implementation of the explicit UGKS can be found in xu2010unified () and xu2015book ().
With a local constant and the assumed local distribution functions in Eq. (4) and Eq. (5), the analytic solution (3) is valid in a local small region on a cell size scale and a small time step from the CFL condition, as illustrated in Fig. 1. Eq. (1), Eq. (2) and Eq. (3) are fundamental physical laws, which presents a local flow evolution in an explicit way. Due to its equationsbased formulation instead of particlebased, Eq. (1) and Eq. (2) can be transformed into semidiscrete forms for the construction of an implicit scheme. The accelerating techniques in CFD can be easily employed here in the UGKS.
2.2 Timeaccurate implicit UGKS
From Eq. (1) and Eq. (2), the semidiscrete macroscopic governing equations for conservative variables can be written as
(6) 
where is the macroscopic flux across the cell interface. The governing equation for the gas distribution function is
(7) 
where is the interface flux of the gas distribution function. The above semidiscrete formulations describe the instantaneous variation of the flow field. It seems to present the flow evolution in a time scale . Although Eq. (6) and Eq. (7) are still consistent with the UGKS formulation, it reduces the timedependent UGKS fluxes into a particle free transport, which makes UGKS as a single scale method within the mean free path and mean collision time. Therefore, in order to incorporate the multiscale property of UGKS into the above formulation, it is better to reconstruct the interface fluxes in Eq. (6) and Eq. (7) as a timeaveraged quantity instead of instantaneous one. Specifically, the above fluxes are averaged over a cellsizedetermined time step, which is given in Fig. 1. Therefore, a cellsizedetermined time scale related flow physics evolution is implicitly included in the semidiscrete governing equations. In order to distinguish this cellsizedetermined time step from the numerical marching time step for an implicit scheme, this time step is given by which determines the local flow physics. is the same as the explicit local time step, which is computed by the local CFL number less than one.
For a large numerical marching time step , we have the discrete macroscopic governing equations
(8) 
and the discrete microscopic equations
(9) 
where and denote the macroscopic fluxes computed by taking moments of the microscopic flux of a gas distribution function. Here and are the timeaveraged fluxes over from the explicit UGKS formulation. For in Eq. (8) and Eq. (9), it corresponds to the CrankNicolson method with secondorder accuracy in time, and for it is a backward Euler scheme with firstorder temporal accuracy.
Eq. (8) and Eq. (9) are fully coupled and in Eq. (9) has onetoone corresponding with in Eq. (8). It is difficult to simultaneously solve these two implicit systems. Instead we can solve Eq. (8) and Eq. (9) in an alternative way. Specifically, we first get an approximate equilibrium state from Eq. (8) from an assumed approximate solution , then obtain a more precise solution of from Eq. (9) using the previously obtained . After several alternations, both and can be solved. However it isn’t straightforward to get solutions from Eq. (8) and Eq. (9) due to the implicit fluxes on the left hand side of these equations. These equations will result in a large nonlinear system, especially for a secondorder accurate numerical fluxes. Therefore, implicit schemes usually adopt a deltaform with implicit parts on left hand side and explicit parts on right hand side of an equation pulliam1993time (). Based on the deltaform equations, the implicit parts on the left hand side can be much simplified using a firstorder flux function, and it becomes much easier to take iterative processes to solve the implicit system tomaro1997implicit (). Details are given in the following.
Given an intermediate approximate solution and , the deltaform governing equations of Eq. (8) and Eq. (9) can be rewritten as
(10) 
and
(11) 
where
(12) 
and
(13)  
where and . The fluxes in the residuals and are fully evaluated by the UGKS with an averaging over a time step , i.e.,
(14) 
and
(15) 
where is given in Eq. (3).
For the terms on the left hand side of Eq. (10), the Euler equationsbased flux
(16) 
is adopted to simplify the implicit macroscopic fluxes, where is the Euler fluxes and is the spectral radius of the Euler flux Jacobian with an additional stable factor related to kinematic viscosity chen2000fast (), i.e.,
(17) 
where is the normal vector of cell interface , is macroscopic velocity and is the speed of sound. For simplifying the numerical fluxes on the left hand side of Eq. (11), we use an upwind approach
(18) 
Substituting Eq. (14) and Eq. (15) into Eq. (10) and Eq. (11), the implicit governing equations for conservative flow variables become
(19) 
where
(20) 
The implicit governing equation for gas distribution function is
(21) 
where
(22)  
For the above deltaform implicit equations, the LUSGS method jameson1987lower (); yoon1988lower (), the pointrelaxation (PR) method yuan2002comparison (), or the GMRES algorithm saad1986gmres () can be applied to solve the variation of conservative variables and the gas distribution function. In the current paper, we employ the PR(2) scheme described in Ref. yuan2002comparison (), for which each inner iterative process is dealt with twice forward and backward sweeps. The details for solving the matrix system will not be further discussed, see zhu2016implicit (); yuan2002comparison (). Therefore, starting with and , Eq. (19) and Eq. (21) can be alternatively solved within several inner iterations.
In summary, the time evolution of the implicit UGKS for unsteady flow simulations from to can be described as
 Step 1

Calculate the numerical fluxes and averaged over using the gas distribution function and equilibrium state by the explicit UGKS fluxes.
 Step 2
 Step 3

Evaluate the macroscopic residual in Eq. (12) using the numerical macroscopic fluxes obtained in Step 1 and Step 2.
 Step 4

Solve Eq. (19) to obtain the equilibrium state .
 Step 5

Evaluate the microscopic residual in Eq. (13) using the equilibrium state recently obtained in Step 4 and the microfluxes and evaluated in Step 1 and Step 2.
 Step 6

Solve Eq. (21) to obtain the gas distribution function , and update the equilibrium state by compatibility condition.
 Step 7

Check the convergence through norm of the macroscopic residual

if convergent state is reached for the current time step, update and into and by and

if not, go to Step 2 and continue the inner iterations.

2.3 Modification for temporal discretization
The timeaveraged fluxes are evaluated over a finite time step instead of at the instant of and . When the evolution time step is comparable with , the size of should be used in the temporal discretization. As illustrated in Fig.2, the fluxes at time can be approximated by those at time and with weighted coefficients
(23)  
Then Eq. (8) and Eq. (9) become
(24) 
and
(25) 
where the modification is only applied to the flux evaluation terms. Then the deltaform of governing equations can be written as
(26) 
and
(27) 
where
and
Particularly, when , equals zero. Then we have
and
which automatically reduce to the explicit UGKS where both macroscopic variables and distribution function are updated within one inner iteration only. For the CN scheme with , the collision term is approximated by trapezoidal rules. It should be noted that without the modified coefficient the early implicit UGKS can not directly go to the explicit scheme. Even though it can still get identical solutions, it may take several inner iterations to obtain the convergent solutions.
3 Numerical validation
3.1 Sod’s shock tube
Sod’s shock tube is computed to validate the IUGKS for unsteady flows in onedimensional case. The initial condition is given by
The time steps are denoted by CFL numbers,
(28) 
For the explicit UGKS, is adopted, which is also used to give the scaledetermining time step for evaluating the timeaveraged fluxes in the implicit scheme.
First of all, we compute the Sod’s shock tube problem at with a small time step of . The computational domain is discretized into 1000 uniform cells. Comparison between the implicit schemes with the original and the modified is carried out when the time step becomes as small as that used in the explicit scheme. For simplification, in the following the IUGKS with will be denoted by IUGKS, and the other one with will be IUGKS. In this case, for the CN method is chosen to give a secondorder accuracy in time, and will become due to . Distributions of density, velocity and temperature at the time are given in Fig. 3. It can be observed that both implicit schemes can give identical solutions to that of explicit UGKS. As to the computational cost, from Fig.3(d) we can find that IUGKS is comparable with the explicit scheme and is much cheaper than IUGKS. This is due to the reason given in Section 2.3, where IUGKS with automatically reduces to the explicit scheme and one inner iteration is needed for each time step. However, for IUGKS several inner iterations are still needed to get a convergent solution in each time step.
Different values of correspond to different temporal discretization methods. Two typical discretization methods for the IUGKS will be presented in the Sod’s test with . Distributions of density, velocity, pressure, and temperature obtained from IUGKS with of the CN scheme, and of the backward Euler method are given in Fig. 4. As observed in the papers luo2001accurate (); tan2017time (), due to the dispersion characteristics of the CN scheme, numerical oscillations will appear around the normal shock wave in the large time step evolution. Regardless of numerical oscillations, the CN scheme can give a sharp discontinuity, while the backward Euler scheme can effectively suppress the numerical oscillations around the smeared normal shock wave. Therefore, for a balance of capturing discontinuity and accuracy preserving, a selection of proper in the interval of is necessary.
In order to figure out how much the implicit scheme can accelerate the numerical simulations, we compute the Sod case using IUGKS with and different time steps. The computational domain is discretized into 400 nonuniform cells with a minimum cell size of and a maximum cell size of . The results in both continuum and rarefied flow regimes are given in Fig. 5, which are compared with the solutions from the exlicit UGKS. For the continuum case at , the velocity space is discretized into 200 velocity points, where the trapezoidal rule is used for the moments integration. The computational cost are listed in Table. 1. It can be seen that generally the IUGKS is more than ten times faster than the explicit UGKS in the continuum regime. For the rarefied flow at , we employ a nonuniform mesh with cells which has a smallest cell size of around the initial discontinuity. In order to get rid of the ray effect and to obtain a smooth solution, we employ 20000 points for velocity space discretization. Actually, using uniform mesh and much fewer velocity points can also give acceptable results for this case, but here we use the much large number of grid points for testing the efficiency of the schemes. The computational cost and acceleration rate from IUGKS are listed in Table. 2 for the rarefied case. It can be seen that for a reasonably good result, e.g. the one with , the computational efficiency can increase by dozens of times.
[c]
Explicit UGKS  IUGKS with  

CFL number  0.5  50  250  450  750 
CPU time (sec)  159.6  23.9  9.2  6.6  5.6 
Speedup  1.0  6.7  17.4  24.3  28.7 
[c]
Explicit UGKS  IUGKS with  

CFL number  0.5  50  200  400  800 
CPU time (sec)  4405.4  270.4  77.7  46.0  29.8 
Speedup  1.0  16.3  56.7  95.8  147.7 
3.2 Couette flow with a temperature gradient
In the incompressible limit, the Couette flow with a temperature gradient has a steady state analytic solution under the assumption of constant viscosity and heat conduction coefficients. Such a steady state problem provides a good test case to validate the spatial accuracy of the IUGKS, where the influence on the time accuracy can be ignored. This case describes the flow between two parallel solid boundaries. The setup is as follows. The top solid boundary has a constant temperature and moves horizontally at a velocity . The bottom wall is fixed and has a lower temperature . If the distance between these two walls is , it gives an analytic solution for the normalized temperature at steady state
(29) 
where is the Eckert number, and is the specific heat ratio at constant pressure.
This case has , , and for the incompressible limit. The argon gas with molecular mass is employed, and . Uniform meshes with , , and discrete cells are used. The results obtained by the IUGKS with in continuum regime are plotted in Fig. 6.
The errors of norms are calculated by
(30) 
where denotes numerical distribution of a specific flow variable, and is the corresponding exact solution. For this case, a normalized temperature distribution is used to evaluate the errors. The norms with respect to mesh size are plotted in Fig. 8(a). It can be seen that the IUGKS has a secondorder accuracy in space for the continuum flow.
In order to measure the spatial convergence of the IUGKS in rarefied regime, the same case at is computed on a group of meshes of , , and cells. The results are given in Fig. 7. Since there is no analytic solution for the Couette flow in rarefied regime, a relative method to measure the error is adopted, i.e.,
(31) 
where and are the numerical solutions on the coarse and the fine meshes, respectively. The results given in Fig. 8(b) show that a secondorder accuracy in space can be obtained as well for the case in the collisionless limit.
3.3 Advection of density perturbation
In order to test the temporal accuracy of the IUGKS for unsteady flow simulations, the onedimension case of advection of density perturbation pan2016efficient () is employed. The initial condition is set as follows
(32) 
The periodic boundary condition is implemented and it gives an analytic solution
(33) 
Since the IUGKS is designed to be secondorder accurate both in space and time, it is supposed to have an error on order of . When is dominant, we can capture the convergence order of the IUGKS with respect to time step. Therefore, in this case we adopt a very fine uniform mesh with cells and use very large time steps. A small mean collision time is used here to drive the IUGKS in a continuum limit to the inviscid case.
Fig. 9 shows the density distribution at time that obtained by the IUGKS with different time steps. It can be predicted that the IUGKS with for the CN temporal discretization has a secondorder accuracy in time because it has been proved that the IUGKS can automatically reduce into the secondorder accurate explicit UGKS for small time step cases. In order to further validate this, the errors for different time step cases are calculated by both Eq. (30) and Eq. (31) to measure the temporal accuracy of the IUGKS. From Fig. 10, it shows that the CN scheme with achieves the expected temporal accuracy.
The accuracy of IUGKS with a temporal discretization of has also been tested. The results obtained by both Eq. (30) and Eq. (31) are given in Fig. 11. As expected, the IUGKS for the case with has a firstorder accuracy in time.
Generally, the stability of IUGKS can be achieved by adjusting the parameter so that the IUGKS is able to evolve in time using any large time step. However, large time steps may lead to unsatisfactory results due to the decrease of time accuracy. Even though large time steps are permissible, the time resolution of flow physics evolutions should be considered. In order to explore the applicability of the IUGKS, we utilize this periodic case on a uniform mesh with cells to further validate the scheme. In Fig. 12(a), we give the results obtained with different time steps, i.e., , , and , where is the period of this sine wave. With respect to the wave length, the relative error of the location of extremum value is about , , and . From Fig. 12(a), it can be found that the shape of the sine wave is well preserved even for the case with , and the solutions can coincide to the exact one after a shift, see Fig. 12(b). Based on this case, we can roughly draw a conclusion that for the requirement of time accuracy, in order to capture the flow variable fluctuations, the numerical time step chosen to evolve the unsteady flow field should not be larger than onetenth of the period of fluctuations.
3.4 Rayleigh flow
The Rayleigh flow is an unsteady gas flow around a vertical plate with infinite length. Initially, the argon gas with molecular mass is stationary and has a temperature of , and suddenly the plate obtains a constant vertical velocity of and gets a higher temperature of . The computational domain is long, which is the characteristic length to define the Knudsen number by the variable hard sphere (VHS) model. The dynamic viscosity is computed by with . Results at time are discussed here.
The target of the IUGKS is to release the time step restriction from smallsize cells on nonuniform mesh and thus to accelerate overall computational efficiency. For this test, we employ nonuniform mesh with the minimum cell size of near the plate. The results at different Knudsen numbers are plotted in Fig. 13, where the density and temperature are normalized by and . In comparison with the DSMC results obtained from the reference paper huang2013unified () and those from the explicit UGKS simulation, IUGKS with can give satisfactory results with a relative large time step. With the continuous increment of time step, the solutions gradually deviate from the reference ones due to the increase of temporal discretization error. For the collisionless limit which requires the interval of discrete velocities satisfy , in order to get smooth solutions velocity points uniformly in direction and 100 points in direction to cover a range from to are used in both directions. Extra more discrete points are employed for discretization of the physical and velocity space in this case, so that it gives more reliable results in the efficiency testing due to long time program running. We present the computational cost in Table. 3 at different Knudsen numbers with various time steps. Generally more inner iterations are required in the small Knudsen number case. The increase of computational efficiency for near continuum flows would not be as much as that in rarefied cases, but it still can be about ten times faster than that of the explicit scheme with . Since the IUGKS is stable enough, we also give the computational cost for case with as a reference even though the solutions under such a condition may not make any sense for a satisfactory solution.
[c]
Kn  Explicit UGKS  IUGKS with  

CFL=0.5  40  100  400  1200  

2.66  437.8  22.5  10.3  4.0  1.9  
0.266  435.0  25.9  11.5  4.8  1.8  
0.0266  436.7  31.0  14.3  7.1  3.3  
0.00266  439.8  48.2  28.2  13.7  11.1  
Speedup  2.66  1.0  19.5  42.3  108.7  236.6  
0.266  1.0  16.8  37.7  91.0  236.8  
0.0266  1.0  14.1  30.6  61.1  132.4  
0.00266  1.0  9.1  15.6  32.1  39.5 
In order to check the sensitivity of the IUGKS to the mesh size, the Rayleigh problems are computed on two kinds of meshes, whose cell size distributions are given in Fig. 14. The mesh denoted by the green diamonds is used in previous calculations, and the one denoted by red circles has a cell size shrinking at . Comparison between results obtained from these two meshes are given in Fig. 15, from which it can be observed that flow can efficiently spread through the small cell size region, and identical solutions to the previous ones can be obtained. A fixed time step used corresponds to the previous case at , which is shown in Fig. 13.
3.5 Wall bounded Rayleigh flow
The Rayleigh problems tested in Section 3.4 are onedimensional, where periodic boundary conditions are imposed in the direction. Here we added two parallel solid walls along xdirection at two locations separated by in the axis, which restrict the movement of the flow in the vertical direction. The schematic for the setting is shown in Fig. 16. The lower and upper solid walls have a length of and a constant temperature of . The argon gas has molecular mass , which corresponds to a specific gas constant of . The initial gas temperature is . The side plates move vertically at with a higher temperature of . The Maxwellian diffusive boundary condition is imposed on the solid isothermal walls. In addition, symmetric boundary condition is imposed at , so that half domain can be used in computation. The case of is considered here.
In order to check the convergence of solutions with respect to the cell size and the time step, numerical simulations are carried out on different meshes using different time steps. The spatial discretizations used in the test case are , , and . In the following, these meshes will be denoted by , , and for simplicity, where means the number of discrete cells along the vertical plate. In this case, the nondimensional time and time step are used, and the normalization is based on a reference time scale , where and . The time steps used in these cases for the IUGKS are , and . As a reference, the cases with are computed as well using the explicit UGKS. Flow variables along the central vertical and horizontal lines ( and ) at are monitored to illustrate the mesh convergence solutions of the unsteady flow, where the results are plotted in Fig. 17 and Fig. 18. It shows that the results on the mesh of have little difference from those on the mesh of , which can be regarded as a mesh convergent solution. Flow variables along central lines obtained using different time steps are given in Fig. 19 and Fig. 20 which show the convergence of the evolving solutions using different numerical time steps.
The flow variables such as magnitude of velocity and temperature at different times are plotted in Fig. 21 and Fig. 22, which clearly present the flow evolution. The mechanism for the flow evolution comes from plate’s shearing and thermal heating effect. Based on the data on the central lines, the rarefied gas effect appears, such as the velocity slip and temperature jump near the solid walls.
The surface quantities, such as pressure, shear stress, and heat flux at four different instants are given in Fig. 23, Fig. 24, and Fig. 25. Detailed data for this case have been listed in the Appendix for future references.
3.6 Hypersonic flow past a square cylinder
This highspeed rarefied flow past a square cylinder is computed to validate the efficiency and stability of the current implicit scheme. The freestream is hypersonic flow of the argon gas at a Mach number with a temperature of . The computational domain is , and the square center locates at the origin. The spatial discretization is the same as that described in chen2017memory (). The density in the freestream is , resulting in a global Knudsen number which is defined with respect to the square side length . The dynamic viscosity is calculated by where .
The density, pressure and velocity distributions at the steady state are given in Fig. 26, where the results obtained from the IUGKS are compared with those from the DSMC method. Generally, the IUGKS results agree well with the DSMC data. For the velocity along direction, the IUGKS gives a straight line on the symmetric plane for the zero contour line while the DSMC gives some noises in the freestream due to the statistic error. The surface quantities, such as pressure, shear stress, and heat flux on the upper half walls of the cylinder are given in Fig. 27 and compared with the DSMC data. Good agreement has been obtained. For unsteady flow simulation, we show the time evolution of the temperature distribution at several instants in Fig. 28 to illustrate the capability of the IUGKS for capturing unsteady flow evolution.
The final steady state solution is obtained at the physical time , which takes about hours using a serial computation. In comparison with the DSMC method, for evolving the solution by the same amount of time, the IUGKS basically takes the same computational cost as DSMC, where the DSMC targets on the steady state solution and no sampling is conducted in the unsteady process. The IUGKS has no noise and doesn’t need any additional sampling for unsteady solution. If only the final steady state solution is needed, a fully implicit with multigrid technique can be used in UGKS to get a much higher efficiency, where the UGKS may become hundreds times faster than DSMC zhu2017multigrid ().
4 Discussions and conclusions
In this paper, an implicit unified gaskinetic scheme (IUGKS) is proposed for unsteady flow simulations in all Knudsen number regimes. By solving the timeaccurate implicit governing equations of the macroscopic variables and the gas distribution function within several inner iterations, the flow field can be updated with a large time step, which is no longer restricted by the CFL stability condition.
For unsteady flow simulation, nonuniform meshes are usually used in many practical engineering applications, such as very fine mesh for viscous flow near solid boundaries and the mesh around the airfoil tip. For the explicit scheme, due to the CFL condition the numerical global time step can become very small, which destroys the advantage of the multiscale nature of UGKS. Since the UGKS employs the analytic solution of the kinetic relaxation model, the numerical flux covers different flow regimes with the changing of the local ratio of the time step over the local mean collision time, i.e., . A uniform time step determined by the smallest cell size will make the UGKS flux to the kinetic scale everywhere, such as . Consequently, a free transport dynamics will be used for the flux evaluation even for these cells with much larger cell size away from the boundaries. Therefore, the multiscale property of the UGKS is not fully utilized. In the current IUGKS, the implicit treatment releases the CFL time step constraint, and the flow field can be updated by a numerical time step which is much larger than the global minimum one. In other words, for the unsteady flow with a large variation of cell size the flow dynamics from IUGKS in different regions can be in different flow regimes. In addition, according to the ratio between the numerical time step and the cellsizedetermined time step , the temporal discretization for the implicit scheme has been improved. As a result, when the global numerical time step in IUGKS becomes as small as the one used in the explicit UGKS, the IUGKS can automatically recover the explicit method, where only one inner iteration for each step is needed. Based on numerical tests, the IUGKS has been validated in terms of accuracy, stability, computational efficiency. The IUGKS has secondorder accuracy both in space and time in both continuum and rarefied flow regimes. Generally, for unsteady flow the IUGKS is at least one order of magnitude more efficient than the explicit UGKS. Even for the hypersonic rarefied flow as shown in this paper, the IUGKS can obtain the time accurate solutions as efficient as DSMC.
Acknowledgment
The work of Zhu and Zhong was supported by the National Natural Science Foundation of China (Grant No. 11472219) and the National PreResearch Foundation of China, as well as the 111 Project of China (B17037). The research work of Xu is supported by the Hong Kong research grant council (16206617,16207715,16211014) and NSFC (91530319,11772281).
Appendix A Data for wall bounded Rayleigh flow
The wall bounded Rayleigh flow tested in Section 3 is a lowspeed rarefied flow, which may become a benchmark for validation of numerical schemes. Here we list the data for at some points along the central lines and on solid walls at time .
(m)  (m/s)  (m/s)  (K)  
