An Implicit Unified Gas-kinetic Scheme for Unsteady Flow in All Knudsen Regimes

An Implicit Unified Gas-kinetic Scheme for Unsteady Flow in All Knudsen Regimes

Yajun Zhu Chengwen Zhong Kun Xu National Key Laboratory of Science and Technology on Aerodynamic Design and Research, Northwestern Polytechnical University, Xi’an, Shaanxi 710072, China Department of Mathematics, Hong Kong University of Science and Technology, Hong Kong, China

The unified gas-kinetic 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 Navier-Stokes 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 non-uniform 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 pre-defined 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 non-uniform 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.

unified gas-kinetic scheme, implicit method, unsteady flow, multiscale transport

1 Introduction

The unified gas-kinetic 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 gas-kinetic 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 non-equilibrium 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 non-equilibrium 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 micro-flows huang2013unified (); liu2015micro (). In addition, the methodology of direct modeling xu2015book (); xu2017paradigm () can be extended to other multi-scale transport processes, such as plasma liu2017plasma (), radiative transfer sun2015radiative (), and phonon heat transfer guo2016discrete (). The UGKS provides a useful tool to solve complex multi-scale problems, and has great potentials in the engineering applications, such as micro-electro-mechanical system and spacecraft design.

One of distinguishable features of the UGKS is to use the scale-dependent 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 chen-ap (). The truly multi-scale 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 Courant-Friedrichs-Lewy (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 multi-scale 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 Navier-Stokes (NS) equations have already been fully developed using dual-time stepping methods jameson1991time (); pulliam1993time (); tomaro1997implicit (). Recently, the time-accurate implicit GKS has also been constructed tan2017time (); li2017implementation () to accelerate computational efficiency for high-speed 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 multi-scale 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 semi-discrete formulation to construct the implicit scheme. To maintain the multiscale transport property, the time derivatives in the semi-discrete 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 lower-upper symmetric Gauss-Seidel (LU-SGS) method jameson1987lower (); yoon1988lower (), point-relaxation (PR) method yuan2002comparison (), and the generalized minimal residual (GMRES) algorithm saad1986gmres (). During the inner iterative process within each time-marching step, the implicit governing equations for macroscopic flow variables and microscopic distribution are solved alternatively with high efficiency for both continuum and rarefied flows.

This paper is organized as follows. Section 2 gives a brief introduction to the explicit UGKS and presents details in the construction of IUGKS. In Section 3, numerical validations and discussions are carried out. A conclusion is given in the last section.

2 Implicit unified gas-kinetic scheme

2.1 Unified gas-kinetic 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


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 time-dependent 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


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 BGK-type kinetic model in the Eq.(1), especially for the cases of liu2016boltzmann (). In the following, the kinetic equation with the BGK-type model will be used only. The multi-scale nature of the UGKS comes from the construction of the time-dependent distribution function at the cell interface, , which is the evolution solution of the BGK model. The analytic solution gives


from the reconstructed initial gas distribution function


and an equilibrium state distributed in time and space


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 non-equilibrium 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 cell-size-determined 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 multi-scale scheme. Details about the implementation of the explicit UGKS can be found in xu2010unified () and xu2015book ().

Figure 1: Physical picture illustrating the time evolution of flow physics around a cell interface in GKS and UGKS.

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 equations-based formulation instead of particle-based, Eq. (1) and Eq. (2) can be transformed into semi-discrete forms for the construction of an implicit scheme. The accelerating techniques in CFD can be easily employed here in the UGKS.

2.2 Time-accurate implicit UGKS

From Eq. (1) and Eq. (2), the semi-discrete macroscopic governing equations for conservative variables can be written as


where is the macroscopic flux across the cell interface. The governing equation for the gas distribution function is


where is the interface flux of the gas distribution function. The above semi-discrete 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 time-dependent 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 multi-scale property of UGKS into the above formulation, it is better to reconstruct the interface fluxes in Eq. (6) and Eq. (7) as a time-averaged quantity instead of instantaneous one. Specifically, the above fluxes are averaged over a cell-size-determined time step, which is given in Fig. 1. Therefore, a cell-size-determined time scale related flow physics evolution is implicitly included in the semi-discrete governing equations. In order to distinguish this cell-size-determined 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


and the discrete microscopic equations


where and denote the macroscopic fluxes computed by taking moments of the microscopic flux of a gas distribution function. Here and are the time-averaged fluxes over from the explicit UGKS formulation. For in Eq. (8) and Eq. (9), it corresponds to the Crank-Nicolson method with second-order accuracy in time, and for it is a backward Euler scheme with first-order temporal accuracy.

Eq. (8) and Eq. (9) are fully coupled and in Eq. (9) has one-to-one 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 second-order accurate numerical fluxes. Therefore, implicit schemes usually adopt a delta-form with implicit parts on left hand side and explicit parts on right hand side of an equation pulliam1993time (). Based on the delta-form equations, the implicit parts on the left hand side can be much simplified using a first-order 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 delta-form governing equations of Eq. (8) and Eq. (9) can be rewritten as








where and . The fluxes in the residuals and are fully evaluated by the UGKS with an averaging over a time step , i.e.,




where is given in Eq. (3).

For the terms on the left hand side of Eq. (10), the Euler equations-based flux


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.,


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


Substituting Eq. (14) and Eq. (15) into Eq. (10) and Eq. (11), the implicit governing equations for conservative flow variables become




The implicit governing equation for gas distribution function is




For the above delta-form implicit equations, the LU-SGS method jameson1987lower (); yoon1988lower (), the point-relaxation (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

Compute the intermediate fluxes and averaged over using the intermediate solution and by Eq. (14) and Eq. (15). For the first inner step, and are adopted.

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 micro-fluxes 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

Figure 2: Structure of temporal discretization.

The time-averaged 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


Then Eq. (8) and Eq. (9) become




where the modification is only applied to the flux evaluation terms. Then the delta-form of governing equations can be written as






Particularly, when , equals zero. Then we have


which automatically reduce to the explicit UGKS where both macroscopic variables and distribution function are updated within one inner iteration only. For the C-N 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 one-dimensional case. The initial condition is given by

The time steps are denoted by CFL numbers,


For the explicit UGKS, is adopted, which is also used to give the scale-determining time step for evaluating the time-averaged 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 C-N method is chosen to give a second-order 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.

(a) Density
(b) Velocity
(c) Temperature
(d) CPU time
Figure 3: Numerical results and computational cost of the Sod’s shock tube problems at obtained from the explicit UGKS, IUGKS- and IUGKS- with . is adopted for the implicit schemes and it gives for IUGKS-.

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 C-N 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 C-N scheme, numerical oscillations will appear around the normal shock wave in the large time step evolution. Regardless of numerical oscillations, the C-N 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.

(a) Density
(b) Velocity
(c) Pressure
(d) Temperature
Figure 4: Comparison of two different temporal discretization methods, i.e., the Crank-Nicolson method and the backward Euler scheme for IUGKS- in the Sod test case with .

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 non-uniform 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.

(a) Density
(b) Density
(c) Velocity
(d) Velocity
(e) Temperature
(f) Temperature
Figure 5: Numerical results of the Sod test case from IUGKS- with different choices of time step. The left ones are the results at , and the right ones are those at .


Table 1: Computational cost of the implicit UGKS and the explicit UGKS for the Sod test case at with different time steps.
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


Table 2: Computational cost of the implicit UGKS and the explicit UGKS for Sod test case at with different time steps.
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


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.

(a) Velocity
(b) Temperature
Figure 6: The normalized velocity and temperature distributions obtained on different meshes by IUGKS with .

The errors of norms are calculated by


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 second-order accuracy in space for the continuum flow.

(a) Velocity
(b) Temperature
Figure 7: The normalized velocity and temperature distribution obtained by IUGKS on different meshes at .
Figure 8: Spatial accuracy of the IUGKS for (a) in continuum flow, and (b) in rarefied 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.,


where and are the numerical solutions on the coarse and the fine meshes, respectively. The results given in Fig. 8(b) show that a second-order 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 one-dimension case of advection of density perturbation pan2016efficient () is employed. The initial condition is set as follows


The periodic boundary condition is implemented and it gives an analytic solution


Since the IUGKS is designed to be second-order 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 C-N temporal discretization has a second-order accuracy in time because it has been proved that the IUGKS can automatically reduce into the second-order 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 C-N scheme with achieves the expected temporal accuracy.

Figure 9: The density distribution at time obtained by the IUGKS with different time steps.
Figure 10: Temporal accuracy of the IUGKS with measured by (a) Eq. (30) and (b) Eq. (31).

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 first-order accuracy in time.

Figure 11: Temporal accuracy of the IUGKS with measured by (a) Eq. (30) and (b) Eq. (31).

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 one-tenth of the period of fluctuations.

Figure 12: Density distribution obtained with , , and . (a) Phase error and (b) shifted results.

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.

(a) Kn=2.66
(b) Kn=0.266
(c) Kn=0.0266
(d) Kn=0.00266
Figure 13: Rayleigh problems at different Knudsen numbers which cover the free molecular, transition and continuum regimes. IUGKS here denotes the IUGKS- with .

The target of the IUGKS is to release the time step restriction from small-size cells on non-uniform mesh and thus to accelerate overall computational efficiency. For this test, we employ non-uniform 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.


Table 3: Computational cost for Rayleigh Problem with different time steps.
Kn Explicit UGKS IUGKS- with
CFL=0.5 40 100 400 1200
CPU time
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
Figure 14: The distribution of the cell size used in simulation of Rayleigh problems.

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.

(a) Kn=2.66
(b) Kn=0.266
(c) Kn=0.0266
(d) Kn=0.00266
Figure 15: Comparison of the results for Rayleigh problems computed on different meshes with a fixed time step .

3.5 Wall bounded Rayleigh flow

The Rayleigh problems tested in Section 3.4 are one-dimensional, where periodic boundary conditions are imposed in the -direction. Here we added two parallel solid walls along x-direction 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.

Figure 16: The schematic for wall bounded Rayleigh flow.

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 non-dimensional 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.

Figure 17: Flow variables along the central horizontal line obtained on different meshes. denotes the discrete cell number along the vertical plate.
Figure 18: Flow variables along the central vertical line obtained on different meshes. denotes the discrete cell number along the vertical plate.
Figure 19: Flow variables along the central horizontal line obtained with different time steps.
Figure 20: Flow variables along the central vertical line obtained with different 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.

Figure 21: Magnitude of velocity distribution for wall bounded Rayleigh flows at different times. Background: explicit UGKS with ; dotted line: IUGKS with time step ; red solid line: streamlines.
Figure 22: Temperature distribution for wall bounded Rayleigh flows at different times. Background: explicit UGKS with ; solid line: IUGKS with time step ; red solid line: heat flux.

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.

(a) Pressure
(b) Shear stress
(c) Heat flux
Figure 23: Surface quantities along the lower wall changing with time. The pressure and shear stress are normalized by where , and the heat flux is normalized by .
(a) Pressure
(b) Shear stress
(c) Heat flux
Figure 24: Surface quantities along the upper wall changing with time. The pressure and shear stress are normalized by where , and the heat flux is normalized by .
(a) Pressure
(b) Shear stress
(c) Heat flux
Figure 25: Surface quantities along the side plate changing with time. The pressure and shear stress are normalized by where , and the heat flux is normalized by .

3.6 Hypersonic flow past a square cylinder

This high-speed 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 .

(a) Normalized pressure
(b) Temperature (K)
(c) X velocity (m/s)
(d) Y velocity (m/s)
Figure 26: Flow variable distributions at steady state for hypersonic flow past a square cylinder at and . The background with dashes lines present the results from the IUGKS, and solid lines in the upper half domain denote the DSMC results.
(a) Pressure
(b) Shear stress
(c) Heat flux
Figure 27: Surface quantities on the upper half walls of the square cylinder. The pressure and shear stress are normalized by where , and the heat flux is normalized by .

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.

Figure 28: Temperature distributions at different times for hypersonic flow past a square cylinder at and ; (a) , (b) , (c) and (d) .

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 gas-kinetic scheme (IUGKS) is proposed for unsteady flow simulations in all Knudsen number regimes. By solving the time-accurate 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, non-uniform 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 multi-scale 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 cell-size-determined 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 second-order 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.


The work of Zhu and Zhong was supported by the National Natural Science Foundation of China (Grant No. 11472219) and the National Pre-Research 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 low-speed 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)