A velocity-space adaptive unified gas kinetic scheme for continuum and rarefied flows

With the notation of cell averaged quantities in the control volume,

where $\Omega_{i,j}$ and $\Omega_{l,m}$ are the cell area in the discretized physical and velocity space, the macroscopic conservative variables and the particle distribution function are updated in the UGKS,

where $\textbf{F}_{r}$ is the flux of conservative variables, $f_{r}$ is the time-dependent gas distribution function at cell interface and $\Delta L_{r}$ is the cell interface length.

In the UGKS, the flux function $\mathbf{F}_{r}$ is evaluated through the particle distribution function at the cell interface, which is constructed from the evolving solution of the Shakhov equation. If the interface location is simplified with $x_{i+1/2}=0$, $y_{j}=0$, and $t^{n}=0$, with a local constant collision time $\tau$, the integral solution in Eq.(4) can be written as,

where $x^{\prime}=-u_{l}(t-t^{\prime})$ and $y^{\prime}=v_{m}(t-t^{\prime})$ are the particle trajectories, and $x^{0},y^{0}$ are the initial locations for the particle which passes through the cell interface at time $t$. Here $f_{0}$ is the gas distribution function at the beginning of $n$-th time step. The internal degree of freedom $\xi$ denotes the random motion in $z$ direction. The above scale-dependent evolving solution provides a multiple scale evolution model, where the contributions from both equilibrium hydrodynamic and non-equilibrium kinetic flow physics are included uniformly.

In the detailed numerical scheme, to the second order accuracy, the initial gas distribution function $f_{0}$ is reconstructed as

where $f_{i+1/2,j,l,m}^{L}$ and $f_{i+1/2,j,l,m}^{R}$ are the reconstructed initial distribution functions at the left and right hand sides of a cell interface, and $\sigma$ and $\theta$ are the slopes of distribution function along $x$ and $y$ directions. In addition, the equilibrium distribution function around a cell interface is constructed as

where $f^{+}_{0}$ is the equilibrium distribution at $(x=0,y=0,t=0)$. The coefficients $(\bar{a}^{L,R},\bar{b},\bar{A})$ can be evaluated from the spatial distribution of conservative variables on both sides of the cell interface and the compatibility condition [19]. After all coefficients are determined, the time dependent interface distribution function becomes

where $\widetilde{f}^{+}_{i+1/2,j,k,l}$ is related to equilibrium state integration and $\widetilde{f}_{i+1/2,j,k,l}$ is related to the initial distribution. With the variation of the ratio between evolving time $t$ (i.e., the time step in the computation) and collision time $\tau$, the interface distribution function above provides self-conditioned multiple scale flow physics across different flow regimes. After the interface distribution function is determined, the corresponding fluxes of conservative variables are evaluated through

Inside each control volume, the collision term $Q(f)$ is to be determined for the update of particle distribution function in Eq.(6). In the unified scheme, the numerical treatments of $Q(f)$ for the exact Boltzmann collision term or the Shakhov model are given in [20, 21, 22, 23].

In continuum flow with intensive intermolecular collisions, the particle distribution function is close to local equilibrium state, and thus the Navier-Stokes equations are valid to describe macroscopic fluid motion. In this case, it is straightforward to apply the Chapman-Enskog expansion to construct the corresponding distribution function at every time step, and thus only macroscopic conservative variables need to be stored and updated. Following this procedure, in the gas kinetic scheme with continuous particle velocity space [24], the gas distribution function $f_{0}$ in Eq.(7) at the beginning of each time step is expanded as following,

where $f_{0}^{+(l)}$ and $f_{0}^{+(r)}$ are the corresponding equilibrium state at left and right sides of the interface, and the coefficients $(a^{l,r},b^{l,r},A^{l,r})$ are related with spatial distribution of macroscopic variables. The equilibrium distribution function $f^{+}$ around a cell interface is constructed in the same way in Eq.(9), and the corresponding interface distribution function becomes,

The interface distribution function here is in the continuous form of particle velocity $(u,v)$, and the fluxes for macroscopic variables can be obtained through the analytical integration of Gaussian distribution.

In a multiscale flow problem, to overcome the computational deficiency and memory burden from a large amount of discretized velocity points, it is feasible to combine both continuum and rarefied flow solvers into a single framework with dynamic continuous-discrete velocity transformation. As shown in Fig.1, in near-equilibrium flow regions, the particle distribution function (PDF) is formulated with a continuous velocity space based on the Chapman-Enskog expansion, and only macroscopic flow variables are updated. In non-equilibrium regions, the AUGKS switches to a discretized velocity space to follow the evolution of particle distribution function. The continuous and discrete velocity space are connected with an adaptation interface, at which the continuous solution of distribution function is sorted into discretized velocity space.

In the detailed numerical scheme, the macroscopic conservative variables are updated in Eq.(5), while in non-equilibrium regions the particle distribution function is updated in Eq.(6). Near the adaptation interface, at every time step $t^{n}$, if there is no recorded discretized distribution function at $t^{n-1}$ in the current ”non-equilibrium” cell $(i,j)$, a local discretized velocity mesh will be generated, where the particle distribution function at velocity point $(l,m)$ is given by the discrete Chapman-Enskog expansion,

where the coefficients $a,b,A$ are related with the spatial distribution of macroscopic variables, and can be determined in the same way in Sec. 3.2. In the current scheme, the velocity mesh is generated within

where $(U,V)$ is macroscopic flow velocity, $T$ is temperature and $R$ is the gas constant.

To update the discretized distribution function in the adjacent cell next to the adaptation interface, the interface distribution function from the continuous GKS solution in Eq.(12) is rewritten into the following discrete form,

and then it can be used to determine the fluxes of macroscopic variables and particle distribution function. In this way, the fluxes at adaptation interface are fully determined and can be used to update the macroscopic variables in Eq.(5) and particle distribution function in Eq.(6).

In the current scheme, the time step is determined by the CFL condition,

where CFL is the CFL number, and $(u_{max},v_{max})$ is the largest particle velocity in $x$ and $y$ directions.

The accuracy and efficiency of the current adaptive scheme are based on a proper choice of location of velocity space adaptation. The transition of discrete to continuous velocity space must be located in the region where the Navier-Stokes solutions provided by the GKS with a continuous velocity space are still valid. For this issue, many empirical parameters for the breakdown of continuum description have been proposed. Bird [4] proposed a parameter $P=D(ln\rho)/Dt/\nu$ for the DSMC simulation of expansion flows, where $\rho$ is gas density and $\nu$ is collision frequency, and the breakdown value of $P$ for translational equilibrium is $0.05$. Boyd et al. [25] extended the above concept to a gradient-length-local Knudsen number $Kn_{GLL}=\ell|\nabla Q|/{Q}$, where $\ell$ is the local mean free path and $Q$ is the macroscopic flow quantity of interest, with the critical value $0.05$. Considering the terms in the Chapman-Enskog distribution function, Garcia et al. [30] proposed a breakdown parameter based on dimensionless stress and heat flux $B=max(|\tau^{*}|,|q^{*}|)$, with the switching criterion of $0.1$.

Since the particle distribution function takes the Chapman-Enskog expansion form in the evolution process of the continuous GKS solver, here we propose an alternative switching criterion of particle velocity space directly from the Chapman-Enskog expansion. For brevity, one-dimensional case will be used for illustration. When there is no discontinuity inside the flow field, the particle distribution function in the Chapman-Enskog form writes,

where the coefficients $a$ and $A$ are the space and time variations of a distribution function, which can be expanded based on the collision invariants,

In the current one-dimensional case, the spatial slope $a$ is determined from the following relation,

where $\psi$ is a vector of collision invariants, $M_{\alpha,\beta}=\int f^{+}_{0}\psi_{\alpha}\alpha\psi_{\alpha}\beta du$ and $\mathbf{a}=(a_{1},a_{2},a_{3})^{T}$. The solution $\mathbf{a}$ can be written as,

The time derivative $A$ is related to the temporal variation of conservative flow variables respectively, and can be evaluated from a successive procedure in the Chapman-Enskog expansion with the help of the Euler equations,

Generally, the Navier-Stokes solutions can be applied when the Chapman-Enskog expansion is a proper approximation of the distribution function in near-equilibrium regime. Therefore, based on the dimensionless collision time and space-time variations, a switching criterion for the velocity space transformation can be defined as

where dimensionless variables are defined as,

where $R$ is gas constant, $L_{0}$ and $T_{0}$ are reference length and temperature. The current switching criterion $B$ for particle velocity space will be tested in numerical experiments.

The numerical algorithm of the adaptive unified gas kinetic scheme is as following. In the AUGKS, we follow the evolution of both both conservative variables and particle distribution function. In near-equilibrium flow regions, the particle distribution function is formulated by the Chapman-Enskog expansion in a continuous velocity space, and macroscopic flow variables are updated in Eq.(5). For non-equilibrium flows, the particle distribution function is updated explicitly in Eq.(6). The scale-dependent flux function is determined by the particle distribution function at the interface, which comes from the integral solution of kinetic model equation in Eq.(4). In each time step, the domain of continuous and discrete velocity space is specified by Eq.(18), and the corresponding interface fluxes are provided by Eq.(10), Eq.(15) and Eq.(12). The detailed numerical procedures for AUGKS are given in Fig. 2.

The first case is the Sod shock tube problem. The flow domain $x\in[0,1]$ is divided into 100 uniform cells. The initial condition is set as

The simulations are performed with reference Knudsen numbers varying from $Kn=0.0001$ to $Kn=1.0$, corresponding to different flow regimes. The current criterion value for velocity space transformation is set as $B=0.0001$. The velocity space is discretized into $80$ uniform points for the update of particle distribution function. In the AUGKS and UGKS, the full Boltzmann collision operator is solved here by the fast spectral method [21]. The numerical solutions at $t=0.2$ are presented in Fig. 3, 4, 5 and 6. The reference solution of continuum flow is calculated by the continuous GKS solver with 1000 cells, and the free molecular flow solution is derived from the collisionless Boltzmann equation.

In the simulation, the region with initial homogeneous spatial distribution of flow variables is calculated with a continuous velocity space, except the central discontinuity simulated with a discretized velocity. As time evolves, the non-equilibrium region enlarges along with the use of discrete velocity space. As presented in Fig. 7a, at $Kn=0.0001$ and $t=0.2$, the flow domain is divided into some subzones, where the non-equilibrium particle distribution function inside rarefaction wave, contact discontinuity, and shock wave is fully resolved with the discretized velocity space, while in the rest near-equilibrium regions the Chapman-Enskog expansion is adopted over a continuous velocity space. The solutions of AUGKS at $Kn=0.0001$ and $t=0.2$ are presented in Fig. 3a, 3b, 3c, which match the benchmark continuum and UGKS solutions accurately. As the Knudsen number gets to $Kn=0.001$ at $t=0.2$, near-equilibrium region confines to a small part near the left tube boundary, where the distribution function has a continuous velocity space, which is shown in Fig. 7b. With increasing rarefaction effect, the distributions of flow variables deviate from the NS solutions gradually and tend to collisionless Boltzmann solutions. As the reference Knudsen number gets to $Kn=0.01$, the Navier-Stokes solutions lose its validity quickly from the initial condition, and the non-equilibrium region occupies the whole tube at $t=0.2$. The numerical solution approaches to the collisionless Boltzmann solution at $Kn=1.0$, as shown in Fig. 5 and 6.

This test case illustrates the capacity of AUGKS to simulate flow in different regimes. The asymptotic preserving (AP) property is confirmed in the two limiting solutions. With increasing reference Knudsen number, there is a smooth transition from the Euler solution of the Riemann problem to collisionless Boltzmann solution. Table 1 presents the CPU time cost from the current adaptive scheme and original UGKS method. As is shown, the AUGKS is about $3.63$ times faster than UGKS at $Kn=0.0001$. Since the computational cost is proportional to the mesh points in the velocity space, it is expected that the computational efficiency is closely related to the size of non-equilibrium region. When the rarefaction increases, the CPU time of AUGKS increases correspondingly, while it is still more efficient than the original UGKS.

A Rayleigh flow forms over a plate which suddenly acquires a constant parallel velocity and temperature. In this test case, we follow the setup by Sun [26]. As shown in Fig. 8, the argon gas is at rest and has a unit temperature initially. When $t>0$, the plate suddenly moves with a constant velocity $U_{w}=0.0296$ and temperature $T_{w}=1.36$. The momentum and energy are transported into the flow field through a shearing effect in the unsteady process. A physical domain $y\in[0,1]$ with 100 uniform cells are set up for the simulation, and the 32 uniform points are used in the velocity space where the particle distribution function is updated directy. In this case, the full Boltzmann collision operator in the AUGKS and UGKS is calculated by the fast spectrum method [21]. The current switching criterion of velocity space is set as $B=0.0001$.

Numerical simulations are performed with a series of reference Knudsen number, and solutions at same output times are plotted in Fig. 9, 10, 11 and 12. Besides AUGKS solutions, the UGKS and DSMC results are also provided as benchmarks. With Maxwell’s fully accommodation boundary condition, Bird [4] proposed an analytical solution from the collisionless Boltzmann equation when the time is much less than the reference mean collision time $\tau_{0}=\ell_{0}/v_{0}$, where $\ell_{0}$ is particle mean free path and $v_{0}$ is the mean molecular speed. The analytical collisionless solution is also plotted in figures.

As presented in Fig. 9, 10, 11 and 12, for the case at $Kn=2.66$ and $t=0.1\tau_{0}$, the AUGKS recovers exact collisionless Botlzmann solution. In the transition regime $Kn=0.266$ and $Kn=0.0266$ at $t=\tau_{0}$ and $t=10\tau_{0}$, the numerical solutions deviate from collisionless solutions gradually due to increasing intermolecular collisions. At $t=100\tau_{0}$ and $Kn=0.00266$ corresponding to a near-continuum regime, the current adaptive scheme recovers the Navier-Stokes solutions with intensive intermolecular collisions. As plotted, in all cases the AUGKS solutions match up well with the benchmark solutions from DSMC and UGKS. It is worth mentioning that in comparison with DSMC method, the current Boltzmann-equation-based adaptive unified scheme has no statistical scattering, which is beneficial in low speed simulations.

Fig. 13 presents the velocity space adaptation inside the flow domain at the output time. In the case with $Kn=0.00266$, in the near-wall region with large slope of macroscopic variables, the AUGKS uses a discrete velocity space, while a continuous velocity space is used in the outer domain. As the Knudsen number increases, the enhanced dimensionless viscosity and heat conductivity lead to a large non-equilibrium region. As a result, the non-equilibrium region enlarges faster. For the case $Kn=0.266$, in all flow region the distribution function deviates from the Chapman-Enskog solution and its evolution must be followed with a discretized velocity space. Table 2 presents the computational cost of the AUGKS and UGKS. When $Kn=0.00266$, the AUGKS is $7.45$ times faster than the original UGKS. When the Knudsen number increases, the enhanced non-equilibrium effects increase the computational cost in AUGKS. From the current numerical experiments, it is clear that the AUGKS provides a self-adjusted algorithm from continuum to rarefied flow simulation with the consideration of both accuracy and efficiency.

The nozzle flow connecting different flow regimes is an ideal case to test the capacity of AUGKS in capturing multiple scale flow dynamics. The schematic of the nozzle problem is presented in Fig. 14. The argon gas is enclosed in a rectangular box $x\in[0,2.2],y\in[-0.5,0.5]$. The velocity space is discretized into $28\times 28$ uniform points for the update of particle distribution function. The switching criterion of velocity space in this case is set as $B=0.0005$. In this case, the collision term in AUGKS and UGKS is the Shakhov model. The computational domain is divided into two parts which are connected through a nozzle. The gas density in the left is $100$ times higher than taht in the right part, and the initial Knudsen numbers are $Kn_{L}=0.0001$ and $Kn_{R}=0.01$ respectively. Except that, the initial gas is at rest and has the same temperature inside two subdomains, same as the cavity wall. All the walls are set up as Maxwell’s diffusive boundary. The nozzle has two entrances with cross sections $\Delta y_{L}=0.13$ and $\Delta y_{R}=0.33$ along a length $\Delta x=0.14$, from which a jet flow will be formed. The simulation is performed till $t=50\tau_{0}$, where $\tau_{0}=\ell_{0}/v_{0}$ is the mean collision time of initial argon gas in the right domain.

Fig. 15 and 16 presents the solution contours of $U-$velocity and temperature at two times $t=5\tau_{0}$ and $t=20\tau_{0}$. The upper part of color contours are the results calculated by the AUGKS (flood) and UGKS (lines), while the lower part is the Navier-Stokes solutions provided by the GKS with a continuous velocity space. As is shown, the bow shock and expansion cooling region behind shock are captured by the two methods. However, it is clear that at $Kn=0.01$, the Navier-Stokes equations lose validity to quantitatively describe the flow evolution in the right domain, and it is necessary to use kinetic method to get accurate solutions here. Fig. 17 and 18 presents the solutions along the horizontal center line of the box at times $t=5\tau_{0}$ and $t=20\tau_{0}$. It is clear that the AUGKS provide equivalent solutions with the NS ones in the near-equilibrium left region, and with the Boltzmann solutions in the non-equilibrium right region. This test demonstrates the multiscale capability of the adaptive method to get the physical solutions in the corresponding flow regimes.

Fig. 19, 20 and 21 present the components of spatial slope $a$ used in the velocity space switching criterion, the mean collision time and the corresponding velocity space adaptation of the flow domain at three times $t=\tau_{0},\ 5\tau_{0},\ 20\tau_{0}$. As shown, the shock and expansion waves are the major sources for large flow gradients inside the domain. Accompanying with the high-density jet into the right domain, the mean collision time decreases in the jet region. With time increasing, the local flow structure becomes more complicated, leading to a large non-equilibrium region. As a result, the Chapman-Enskog expansion fails in the places where the strong non-equilibrium effects appear, and the discretized velocity space has to be used in AUGKS. Table 3 presents the computational cost of AUGKS, UGKS and GKS in this case. With the current setup of physical mesh and velocity space, the continuous GKS solver is about 25 times faster than the UGKS with discretized velocity space. The AUGKS is $1.7$ times faster than the original UGKS. At the initial stage in the simulation, the memory cost in the AUGKS is on the same order as the GKS, which is about $1/30$ of the UGKS. As flow evolves, the number of cells associated with discretized velocity space increases, and the corresponding memory cost becomes larger. At the final time $t=50\tau_{0}$, there are about 429 out of 914 total cells using continuous velocity space, and the corresponding memory burden is $53\%$ of the original UGKS.

In the previous studies, the AUGKS in the simulation of unsteady flows is well validated. In this case, the flow passing through a circular cylinder is used to test the performance of current adaptive scheme for steady flow. The incoming gas has a uniform velocity with Mach number $Ma=5$ and the same temperature as the cylinder wall. The reference Knudsen number is set up as $Kn=0.001$ and $Kn=0.01$ relative to cylinder radius, and the corresponding dynamic viscosity is $\mu_{ref}=7.313\times 10^{-4}$ and $\mu_{ref}=7.313\times 10^{-4}$. In the calculation, $60$ cells in radial direction and $100$ cells in circumferential direction are used in physical domain. The velocity space is discretized into $41\times 41$ velocity points for the update of particle distribution function. In this case, the collision term for particle distribution function in AUGKS and UGKS is modeled as the Shakhov model. The cylinder’s wall is set up as Maxwell’s diffusive boundary. The switching criterion of particle velocity space is set as $B=0.0005$.

For steady problem, the computational time can be further reduced with the help of the GKS with a continuous velocity space. A convergent coarse flow field can be first obtained by the GKS, and then used as the initial state in the subsquent adaptive method. The method for the computation of steady flow is the following.

1. 1.

   From initial setup, use the GKS solver in the entire domain and obtain a convergent flow field.

2. 2.

   Use the calculated macroscopic flow variables as the initial flow condition to get the particle distribution function with the discretized Chapman-Enskog form by Eq.(13).

3. 3.

   Adapt the velocity space based on the current switching criterion in Eq.(18).

4. 4.

   Use continuous velocity space in near-equilibrium flow region and discretized velocity space in non-equilibrium one, and redo the velocity adaptation after some iterations until a convergent flow field is obtained.

Fig. 22 and Fig. 23 presents the solution contours of $U-$velocity and temperature calculated by the AUGKS, UGKS and GKS methods respectively around the cylinder. The upper part of contours is the results of AUGKS (flood) and UGKS (solid line) solutions, and the lower part is the GKS solutions. As shown, the bow shock and expansion cooling region behind shock are well captured by all methods.

Fig. 24, 25, 26, 27 present the solutions along the horizontal center line in front of and behind the cylinder. At $Kn=0.001$, the cell size and time step in the computation are much larger than particle mean free path and collision time. Due to the limited time-space resolution, all the three methods become shock-capturing schemes, and a sharp shock profile is obtained in front of the cylinder in Fig. 24. Near the cylinder wall, due to the non-equilibrium gas dynamics in gas-surface interaction, there is a slight difference in the solutions provided by UGKS and GKS. At the same time, the gas density reduces much in the wake region behind cylinder with emerging rarefied effects, and there is a significant difference between UGKS and GKS solutions in Fig. 25.

When the reference Knudsen number gets to $Kn=0.01$, the increasing particle mean free path collision leads to a wider shock structure. This non-equilibrium evolution is provided in the scale-dependent interface solution used in AUGKS and UGKS. However, the Chapman-Enskog expansion can only provide incomplete information about this process in continuous GKS solver. As a result, the GKS presents a narrower shock profile than that in AUGKS and UGKS in Fig. 26. In the wake region, due to the increasing collision time at $Kn=0.01$, the results provided by GKS differ significantly from AUGKS and UGKS solutions in Fig. 27. The GKS with a continuous velocity space fails to predict physical solutions in these regions, and particle distribution function is updated explicitly in the AUGKS with a discretized particle velocity space. It is clear that the current velocity adaptive unified scheme captures the equivalent physical solutions as the NS in the near-equilibrium region and the UGKS ones in the non-equilibrium region.

Fig. 28 and 29 presents two components of spatial slope $\mathbf{a}$ used in the velocity-space switching criterion, the mean collision time and the adaptation of velocity space. As can be seen, the shock wave and boundary are two sources for high gradients of flow variables, leading to the failure of Chapman-Enskog expansion and Navier-Stokes solutions. Behind the cylinder, the low density wake leads to an increasing collision time, which is shown in Fig. 28c and Fig. 29c. Therefore, a velocity adaptation is determined as shown in Fig. 28d and Fig. 29d . The incoming flow as well as a small region between the bow shock and cylinder is computed with continuous velocity space, while the rest non-equilibrium region are simulated by the UGKS with discrete velocity space. Due to increasing collision time and non-equilibrium flow dynamics, the region under discrete particle velocity space enlarges at $Kn=0.01$ than that at $Kn=0.01$.

Table 4 and 5 presents the computational cost of AUGKS, UGKS and GKS at $Kn=0.001$ and $Kn=0.01$ respectively. With the current setup of physical mesh and velocity space, the continuous GKS solver is about 30 times faster than the UGKS, and the AUGKS is about $3.3$ times faster than the original UGKS in this steady flow problem. In the convergent steady state, there are about 3196 at $Kn=0.001$ and 1992 at $Kn=0.01$ out of 6000 total cells using continuous velocity space, and the corresponding memory burden is about $47\%$ and $67\%$ of the original UGKS.