Progress of discrete unified gas-kinetic scheme for multiscale flows

Multiscale gas flows appear in many fields and have received particular attention in recent years. It is challenging to model and simulate such processes due to the large span of temporal and spatial scales. The discrete unified gas kinetic scheme (DUGKS) is a recently developed numerical approach for simulating multiscale flows based on kinetic models. The finite-volume DUGKS differs from the classical kinetic methods in the modeling of gas evolution and the reconstruction of interface flux. Particularly, the distribution function at a cell interface is reconstructed from the characteristic solution of the kinetic equation in space and time, such that the particle transport and collision effects are coupled, accumulated, and evaluated in a numerical time step scale. Consequently, the cell size and time step of DUGKS are not passively limited by the particle mean-free-path and relaxation time. As a result, the DUGKS can capture the flow behaviors in all regimes without resolving the kinetic scale. Particularly, with the variation of the ratio between numerical mesh size scale and kinetic mean free path scale, the DUGKS can serve as a self-adaptive multiscale method. The DUGKS has been successfully applied to a number of flow problems with multiple flow regimes. This paper presents a brief review of the progress of this method.


Introduction
Multiscale gas flows appear in many natural and industrial systems, such as nano/micro devices, aerospace vehicles, vacuum techniques, and unconventional natural gas exploitation. Such systems usually involve a large span of length and time scales, which brings challenges in modelling and simulations to capture the flow physics in different scales under a unified framework.
Physically, the transport behaviors are determined by the microscopic dynamics of the underlying gas molecules, which can be further modelled efficiently at different scales. It is well understood that each molecule in a gas system undergoes free streaming and collision (or scattering) dynamics, and a straightforward tracking of these two processes leads to the molecular dynamic (MD) model at the microscopic (molecular) scale. On the other hand, the collective behaviors of the gas molecules at large length (l 0 ) and time (t 0 ) scales can be described macroscopically through phenomenological models, such as the Euler or Navier-Stokes equations. In between the molecular scale and the macroscopic (hydrodynamic) scale, the gas evolution can be described by the Boltzmann equation that models the evolution of velocity distribution function of gas molecules on kinetic scales, i.e., mean-free-path λ and mean collision time τ = λ/c 0 , with c 0 being a typical molecular velocity which is of the order of sound speed [1,2]. The MD simulation tracks the motion of individual molecule in the system and thus is quite computationally intensive, which limits its application to short time dynamics of systems with the resolution of molecular diameter. On the other hand, the macroscopic models, although relatively easier to handle, are limited to large scale systems under the assumptions of continuum mechanics and local thermodynamic equilibrium. Therefore, both the microscopic and macroscopic models are inadequate for modelling the multiscale transport due to the existence of multiple length and time scales.
It is well understood that the Boltzmann equation can lead to hydrodynamic equations in the asymptotic limits via perturbation analysis in terms of the Knudsen number Kn = λ/l 0 (Note that if we choose t 0 as the acoustic time scale t 0 = l 0 /c 0 , Kn can also be expressed as Kn = τ/t 0 [2]). Therefore, the Boltzmann equation provides a solid basis for developing uniformly accurate and stable numerical schemes for gas flows in all regimes from kinetic (Kn 1) to hydrodynamic ones (Kn 1). Actually, a large number of numerical methods for kinetic equations have been developed from different points of view, such as the direct simulation Monte-Carlo (DSMC) [1], discrete velocity method (DVM) [3,4], lattice Boltzmann equation (LBE) [5], gas-kinetic scheme [6], semi-Lagrange method [7], and implicit-explicit (IMEX) method [8]. The progress of the numerical schemes based on kinetic equations can be found in several review papers [9][10][11].
For a kinetic scheme, its capability of simulating multiscale flows closely depends on its asymptotic behavior at small Knudsen numbers, i.e., whether it can capture the hydrodynamic behaviors without resolving the kinetic scale. A scheme with this property is usually called as an "asymptotic preserving" (AP) one [11], which was originally defined for the Euler limit (Kn → 0) and later was used to analyze schemes in the Navier-Stokes limit. A more elaborate concept, unified preserving (UP), was proposed recently [12], with which the detailed asymptotic order of a kinetic scheme can be assessed. In the past years, a number of kinetic schemes with Euler or Navier-Stokes asymptotics have been developed, e.g., [13][14][15][16][17][18][19].
Generally, for a multiscale kinetic scheme it is necessary to preserve the Navier-Stokes asymptotics with a numerical resolution (mesh cell size x and time step t) much larger than kinetic scale, or, it should be at least second-order UP. In the continuum flow regime, the UP schemes should keep the same properties as the shock capturing schemes designed for the Navier-Stokes equations directly in the calculation of hydrodynamic wave structure, such as the boundary layer with the resolution of a few mesh points. From this perspective, the DSMC and classical explicit DVM are not good choices for multiscale flows in that the cell size and/or time step are required to be smaller than the mean-freepath and relaxation time, respectively, which is a severe limitation for near-continuum flow computation. Indeed, both methods are mainly used to simulate rarefied gas flows. On the other hand, the LBE and GKS, with implicit and explicit Chapman-Enskog approximations, respectively, are designed mainly for continuum flows, and therefore are also inadequate for multiscale flow simulations. Some UP schemes, which aim to capture flow behaviors in all regimes, have also been developed in the past decade [18][19][20][21][22][23][24]. Particularly, the finite-volume unified gas-kinetic scheme (UGKS) [18] has gained much attention due to its special reconstruction of cell-interface flux, in which the analytical time evolved integral solution of the kinetic equation is adopted to approximate the distribution function at cell-interface. With similar principle, the discrete UGKS (DUGKS) was developed later [19], in which a simpler numerical characteristic solution of the kinetic equation is employed to reconstruct the cell-interface flux. The discrete evolution of the cell-interface distribution function along the characteristic line resembles that of the LBE, and thus the DUGKS can be viewed as a combination of the LBE and GKS methods.
With the consideration of both free streaming and collision physics in the flux reconstruction, it can be shown that the DUGKS actually solves the collision-less Boltzmann equation as Kn → ∞, and the Navier-Stokes equations as Kn 1 even as x λ 0 and t τ 0 [12]. The transition between flow regimes is realized automatically based on the local flow physics accumulated from the particle transport and collision within a time step, and therefore the DUGKS can be regarded as a self-adaptive multiscale method, which is critical for simulating different regime flow behavior in a single computation. Furthermore, the full temporal-spatial coupling also makes the DUGKS a numerical method respecting the conservation-dissipation mechanics, which is a fundamental requirement for irreversible processes [25]. Even in continuum regime, the finite-volume formulation and the release of tight coupling between time step and mesh size make the DUGKS a competitive tool in comparison with LBE. With these nice properties, the DUGKS has been successfully applied to a variety of flow problems in different flow regimes, such as turbulent flows [26][27][28], micro flows [29][30][31][32], compressible flows [33][34][35], multiphase flows [36,37], gas-solid flows [38,39], and gas mixture systems [40,41]. Besides flow problems, the DUGKS was also extended to multiscale transport problems such as phonon heat transfer [42][43][44] and radiation of photons [45,46]. In this paper, we will give a brief review on the recent progress of the DUGKS. The rest of the paper is organized as follows. In Section 2, the basic structure of DUGKS is presented, together with a brief discussion on its fundamental properties. Section 3 shows a comparison between the DUGKS and LBE for continuum flows, where the DUGKS can be viewed as a special FV-LBE in this regime. In Section 4, a comparison between the DUGKS and UGKS is presented, focusing on the structure of the time averaged interface distribution functions reconstructed in both methods. A number of sample applications of DUGKS to single-phase gas flows in different flow regimes are demonstrated in Section 5, and some extensions of DUGKS to more complex fluid flows are reported in Section 6. In Section 7, some extensions of DUGKS to multiscale transport problems beyond gas flows will be shown. Finally, a brief summary and outlook of the DUGKS are given in Section 8.

Formulation
The DUGKS solves the Boltzmann equation with relaxation models. Without loss of generality, here we take the BGK model for monatomic gases as an example to illustrate the basic idea. The BGK model reads [47], where f = f (x, ξ , t) is the velocity distribution function for molecules moving with velocity ξ at position x and time t, τ is the relaxation time, and f (eq) is the equilibrium distribution function given by, where ρ, u, and T are the density, velocity, and temperature of the gas, respectively, D is the spatial dimension, and R = k B /m is the gas constant with k B the Boltzmann constant and m the molecular mass. The conservative flow variables W = (ρ, ρu, ρE) T are defined as the moments of the distribution function, where ψ = 1, ξ , ξ 2 /2 T is the vector of the elementary collision invariants, and ρE = ρ u 2 + DRT /2 is the total energy. The viscous stress τ and heat flux q can be defined as The DUGKS is a finite-volume discretization of Eq. (1). First, the flow domain is divided into a set of control volumes or cells. Integrating Eq. (1) on a control volume V j centered at x j from time t n to t n+1 = t n + t with a time step t, one can obtain that where the midpoint rule for the time integration of the convection term and trapezoidal rule for the collision term are used, respectively; f n j is the cell-averaged distribution function, with |V j | being the volume of cell V j . The flux F n+1/2 j is evaluated via the midpoint rule for the surface integration, i.e., where ∂V j is the surface of cell V j and n is the outward unit vector normal to the surface, S jk is the surface area of the interface between cell V j and its k-th neighboring cell V k , x jk and n jk are the corresponding face center and unit normal vector of S jk , respectively, as shown in Fig. 1. It is noted that in the volume-average collision term, Q j , the averaged equilibrium distribution function is approximated as f where ρ j , u j , and T j are the cell-averaged density, velocity, and temperature, respectively. This approximation guarantees the conservative properties of the average collision term.
The key feature of DUGKS lies in the reconstruction of the interface flux, which is determined by the distribution function f x jk , ξ , t n+1/2 . As shown in Fig. 1, the intermediate interface distribution function f x jk , ξ , t n+1/2 is obtained by integrating the kinetic Eq. (1) along its characteristic line dx/dt = ξ from x = x jk − ξ h to x jk with a half time step h = t/2, where the trapezoidal rule is again applied to the collision term. Equation (8) can be rewritten in an explicit form, wherē Note that Eq. (9) is just the same as the collision-streaming process in the lattice Boltzmann equation, except that the time step is t/2 now.f + represents the postcollision state at the starting point x of the characteristic line, and then it moves to the face center x jk with a half time step. In Eq. (9),f + x , ξ , t n is reconstructed from the cell-averaged values of the neighboring cells. For smooth flows, it can be approximated as a linear function around x jk , i.e., wheref + x jk and the slope σ jk = ∇f + x jk are approximated by linear interpolations. On the other hand, for flows with discontinuities,f + is assumed to be piecewise linear in each cell, and thusf + x , ξ , t n is determined from the upstream cell, where σ i is the slope ofφ + in cell V i (i = j or k), which can be reconstructed from the cell-averaged values using certain numerical limiters [48].
After determiningf + x , ξ , t n , the half-time distribution functionsf x jk , ξ , t n + h can be obtained from Eq. (9), from which we can then get the conservative variables based on Eq. (10), where we have used the conservative properties of the collision term. Then, the original distribution function can be obtained from Eq. (10), (14) and then the interface flux F n+1/2 j can be determined based on Eq. (7). After obtaining the flux, we can update the cell-averaged distribution function f j . Note that Eq. (5) is implicit due to the involvement of Q n+1 j . Two approaches can be employed to implement Eq. (5) explicitly. The first approach is to introduce two new distribution functions [19,33], Then the evolution Eq. (5) can be rewritten as In practical simulations, we can just trackf instead of the original one since the hydrodynamic variables can be obtained fromf , i.e., where β = τ/(τ + t/2). Another approach for explicit implementation of Eq. (5) is similar with that used in the IMEX [8] and UGKS [18] methods, which is realized by first taking the conservative moments of Eq. (5) to obtain and then the updated distribution function can be calculated as This approach was also adopted in some later conserved DUGKS methods [49,50]. An advantage of this approach is that macroscopic conservation (in terms of W ) from the discrete collision term is satisfied exactly, regardless of the quadrature rule for the integral in velocity space. However, it should be pointed out that the microscopic conservation in terms of f is not necessarily satisfied, as discussed in Section 2.2. In summary, the numerical procedure of each time step in DUGKS can be executed as follows: (1) Reconstructingf + x , ξ , t n from cell-averaged distribution functions for each cell interface S jk and velocity ξ ; (2) Movingf + x jk − ξ h, ξ , t n to the face center along the characteristic line to obtain f x jk , ξ , t n + h ; (3) Calculating W n+1/2 jk fromf x jk , ξ , t n + h , and then obtaining f x jk , ξ , t n + h to give the flux F j ; (4) Updatingf j or f j from t n to t n+1 for each cell.
The above procedure for one-dimensional case can be illustrated in Fig. 2, where the first approach for updating the distribution function is used. It is noted that the BGK model gives a fixed Prandtl number (Pr = 1). Some other relaxation models, such as the ellipsoidal statistical model [51] and the Shakhov model [52], can overcome this deficiency. Later the DUGKS was generalized to variable Prandtl number problems based on the Shakhov model [33], which has a similar structure as the above algorithm.

Velocity discretization
In the above derivation we focus on space and time discretizations. In practical computations, the molecular velocity space should be discretized into a finite set of discrete velocities ξ i |i = 1, 2, · · · , N , like the classical discrete velocity methods. Generally, two types of discrete velocities are used. The first one chooses the abscissas of certain Gaussian quadratures, such as the full or half-range [53] Gauss-Hermite rules [54], as the discrete velocity set, which spans the continuous velocity space nonuniformly; Another type uses a Cartesian velocity grid with a bounded range dependent on temperature and flow velocity. Unstructured velocity mesh was also developed recently [50], and some locally adaptive techniques were proposed to enhance memory and computational efficiency [55,56].
With the discrete velocity space, the velocity moments of the distribution function will be replaced by the corresponding discrete ones, e.g., where w i (i = 1, 2, · · · , N) are the weights of the quadrature. Provided the number of discrete velocities is large enough, the numerical quadrature can be quite accurate. It should be noted that, however, the discrete collision operator is generally not conservative strictly, i.e., k ψ ξ k Q ξ k = 0. This inconsistency comes from the discrete equilibrium distribution function with the original conservative hydrodynamic variables W , because generally This inconsistency is reasonable, since a discrete velocity model evolving in a discrete phase space may have its own equilibrium state that is different from the continuous counterpart. A way is to find the discrete equilibrium based on discrete version of entropy minimization [4], i.e., the regularized discrete equilibrium distribution is defined as f This nonlinear equation can be solved using the Newton iteration method. It has been shown that with the discrete equilibrium distribution defined above, the required number of discrete velocities for DVM can be much reduced to obtain the same accuracy [57]. However, additional computational cost is needed for solving Eq. (22). For low speed or small Mach number flows, a simpler equilibrium can be obtained by means of the Hermite expansion [5,58]. First, the continuous Maxwellian equilibrium distribution is projected onto the space spanned by the Hermite polynomials, where ω ξ = e −ξ 2 /2 /(2πRT 0 ) D/2 ,ξ = ξ /(RT 0 ) 1/2 with T 0 being a reference temperature, H (k) (k = 0, 1, · · · ) are the Hermite polynomials, and P K is a K-th order polynomial ofξ . The expansion coefficients are given by For instance, up to third-order, the expanded equilibrium distribution can be written as whereû = u/(RT 0 ) 1/2 andθ = T/T 0 . Since P K is a K-th order polynomial ofξ , if we choose a Gauss-Hermite quadrature [with weight function ω ξ ] of order higher than K + 2, the conservative moments of this expanded equilibrium distribution can be evaluated exactly, where W i and ξ i /(RT 0 ) 1/2 (i = 1 ∼ N) are the corresponding weights and abscissas of the Gauss-Hermite quadrature, respectively, and w i = W i /ω ξ i (i = 1 ∼ N) are the weight for evaluating the discrete moments (20). The expanded equilibrium distribution is very similar to that used in the LBE, but for the DUGKS off-lattice discrete velocities can be employed naturally, which brings much freedom for the choice of quadrature rules. It should be noticed that although the expanded formulation of equilibrium distribution f (eq) K can ensure the microscopic conservation of the collision operator, it may take negative values and does not minimize any entropies.

Boundary conditions
For flows with solid walls, appropriate boundary conditions should be specified for the discrete distribution functions at the wall surface in the evolution of DUGKS. A general boundary condition is the diffuse-scattering rule, which assumes the distribution function for the reflected molecules follows the Maxwellian one with the wall temperature and velocity. Specifically, at the half-time step t * = t n +h, the unknown distribution functions are given by where x w is the center of the cell interface located at the wall, n is the unit vector normal to the wall pointing to the cell, u w and T w are the velocity and temperature specified at the wall, c i = ξ − u w , and ρ w is the gas density determined by the condition that no molecules can go through the wall, i.e., where the distribution functions f x w , ξ i , t * with c i · n < 0, i.e., for molecules moving towards the wall, can be constructed following the procedure described in Section 2.1. Equation (28) gives that (29) [Note that there is a typo for the diffuse boundary condition in Ref. [19], although it is correctly implemented in the code: the weight w i appearing in each summation in Eqs. (28) and (29) is missing there]. The above diffuse-scattering boundary condition can also be extended to partially diffusive walls with an accommodation coefficient 0 ≤ σ ≤ 1, where ξ i = ξ i − 2n(c i · n) + u w is the velocity of the specularly reflected molecules with incident velocity ξ − u w . For low speed and isothermal flows, the bounce-back method can be employed to realize no-slip velocity boundary condition, which just reverses the molecular velocity after hitting the wall moving with velocity u w , where the weights W i satisfy i W i = 1 and i W i ξ i ξ i = RT 0 I, and ρ w is given by

Basic properties of the DUGKS
We now discuss some fundamental properties of the DUGKS, including the consistency, accuracy, and uniform stability in terms of the relaxation time. Particularly, the asymptotic behaviors of DUGKS for small Knudsen number, which is critical for capturing the correct hydrodynamic physics without resolving the kinetic scale, will also be analyzed. The use of the midpoint and trapezoidal rules in Eqs. (5) and (8) for the time integration, as well as the linear spatial reconstruction of the distribution function at the cell interface, ensures the second-order accuracy in both space and time of the DUGKS. This is more evident from the modified equation. For simplicity without loss of generality, we consider one-dimensional case and assume the flow is smooth. Then by applying the Taylor expansion to Eqs. (5) and (8), we can obtain the modified equation of the DUGKS [12], where This confirms the DUGKS is a consistent second-order scheme for the BGK equation in both time and space for a fixed τ .
Regarding the uniform stability of DUGKS, we can see that the collision term is treated semi-implicitly in Eqs. (5) and (8). Therefore, it is expected that the time step t should not be limited by the relaxation time τ . But the convection term is still treated explicitly, and thus t should satisfy the Courant-Friedrichs-Lewy (CFL) condition, where 0 < η ≤ η 0 is the CFL number with η 0 being some constant, |ξ | max and |u| max are the maximum values of the discrete velocities and flow speed, respectively. The time step not only influences numerical stability, but also contributes to numerical dissipation in kinetic schemes, such as the gas-kinetic scheme [6]. The numerical dissipation from the time integration was also analyzed in [59]. Roughly, the asymptotic property of the DUGKS can be demonstrated by analyzing the reconstructed distribution function at cell interface. From Eq. (8), it can be deduced that f x jk , ξ , t n + h = f x jk − ξ h, ξ , t n as τ/ t → ∞, which is just the solution of the collisionless BGK equation. On the other hand, in the continuum limit where τ/ t → 0, it can be shown that [19] which is just the Chapman-Enskog solution at the Navier-Stokes order. These results indicate that the DUGKS can indeed preserve the asymptotic behaviors in both continuum and free molecular limits. In the transitional regime, it is also expected that the DUGKS can give a good approximate solution to the kinetic equation since it is a consistent discretization. In [12], it is shown more rigorously that as x/λ and t/τ are of the order of √ Kn, the DUGKS preserves the Navier-Stokes limit exactly. This means that the DUGKS can yield the Navier-Stokes solutions as x λ and t τ for small Kn. The above arguments suggest that the DUGKS is a self-adaptive multiscale method for flows covering different flow regimes.
In the early FV-LBE methods (e.g., [60,61]), simple interpolations are employed to reconstruct the distribution function at cell-interface and the collision term is treated explicitly, such that the time-step is severely limited by the relaxation time in addition to the CFL condition [62]. Some later FV-LBE methods improved the numerical stability by employing certain upwind interpolations [63,64]. In most of available FV-LBE methods, the construction of flux at cell interfaces ignores the effect of collision term, and thus could lead to large numerical dissipation. We note that a FV-LBE similar to the DUGKS was developed more recently [64], which employs an explicit discrete characteristic solution of the kinetic equation in the flux reconstruction such that the collision effect is taken into account. This FV-LBE has similar accuracy with the DUGKS, but the numerical stability is degenerated due to the explicit treatment in the flux reconstruction [64]. Some comparative studies on the performance of DUGKS and LBE have been conducted in the literature. For instance, Wang et al. made a comparison of accuracy, stability, and efficiency of the DUGKS and two standard LBE models, i.e., the lattice BGK (LBGK) and multiple-relaxation-time (MRT) models [66]. They simulated the two-dimensional (2D) laminar cavity flow and the flow around a square cylinder at different Reynolds numbers. The results showed that the two LBE models and the DUGKS could yield qualitatively similar results in both test cases with sufficient mesh resolutions. However, some unphysical oscillations in pressure field were observed in the results of both LBGK and MRT models, but the pressure field predicted by the DUGKS was smooth. Furthermore, it was found that the LBE methods could yield inaccurate vortex shedding behaviors for the flow over a square cylinder with a blockage ratio of 1/8. Specifically, with a uniform mesh of size 2000×320, the frequencies of vortex shedding predicted by both LBE methods took a single value, and the flow patterns were alternatively symmetric. But as the mesh resolution increased to 5000 × 800, the vortex shedding demonstrated a multiple frequency style and the flow pattern became asymmetric, which was in agreement with results from a Navier-Stokes solver [67]. On the other hand, the multiple vortex shedding frequency and asymmetric flow pattern could be captured by the DUGKS with the uniform coarse mesh as well as a nonuniform one. The flow patterns predicted by the three methods are shown in Fig. 3. The comparison also showed that the numerical stability of the DUGKS and the two LBE methods could be considerably different. Generally, with a same uniform mesh, the DUGKS could simulate flows at much higher Reynolds numbers than the LBE methods. The computational efficiency was also compared in [66]. It was shown that with the same mesh resolution, the LBE methods were more efficient than DUGKS. However, the efficiency of DUGKS could be much improved by using non-uniform meshes.
The performance of DUGKS was also compared with a characteristics-based LBE (C-LBE) in finite-volume formulation [59]. This C-LBE model was originally developed as a characteristic discretization of DVM [68], and was shown to be able to overcome the time step restriction from the relaxation time. The comparison with several other nonstandard LBE methods indicated that the C-LBE is the most stable and accurate one [69]. Zhu et al. showed that this C-LBE could be re-expressed in a finite-volume formulation, where the collision effect was also incorporated in the flux reconstruction. Specifically, the distribution function at a cell interface of the C-LBE is approximated as [59], Comparing this with Eq. (8), we can see that the only difference between the DUGKS and C-LBE lies in the evaluation of the collision term in calculating the characteristic solution at cell interfaces. The C-LBE treats the time integral of the collision term with the explicit one-point quadrature when integrating the BGK equation along the characteristic line, while the DUGKS evaluates the integral of collision term with the semi-implicit trapezoidal quadrature. Therefore, it can be expected in principle that the DUGKS can be more accurate and stable than the C-LBE. Several numerical test results, including the cavity flow, laminar boundary layer flow, and the unsteady Taylor-Green vortex flow, confirmed the above analysis. For instance, for the Taylor-Green vortex flow, the L 2 errors of DUGKS are about 1/4 times of those of C-LBE with the same uniform meshes. The tests of the steady cavity flow and the laminar boundary layer flow over a flat plate also demonstrated that the DUGKS could give more accurate results than the C-LBE. It was found that the flat laminar boundary layer could be captured accurately by the DUGKS and the results were insensitive to the computational meshes. In particular, with a coarse mesh ( y min = 0.1) which has only 4 cells in the boundary layer, the results were still rather satisfactory. On the other hand, the C-LBE could only give satisfactory results with a fine mesh ( y min = 0.01). The C-LBE was also found to be sensitive to grid resolution. With the coarser mesh, the thickness of the boundary layer was over predicted, indicating that the C-LBE is more dissipative. The standard LBE was also employed to simulate the laminar boundary layer flow, and it was shown that the LBE could give an accurate solution with the coarse mesh ( y min = 0.1), but the computational time was much expensive than the DUGKS due to the use of uniform mesh.
The numerical results also showed that the DUGKS is more stable than the C-LBE. It was shown that the computation of the C-LBE was unstable at moderate values of t/τ even though η < 1, but the stability of DUGKS was almost not affected by the CFL number provided η < 1.1. The numerical efficiencies of the two methods were also measured. It was found that the computational time of DUGKS was about twice as that of the C-LBE with the same mesh. However, the DUGKS can achieve an accurate solution with a much coarser mesh, and thus it would be more efficient than C-LBE to obtain a reliable solution.
In a recent study, the numerical stability of DUGKS and another advanced FV-LBE model were also compared by simulating the laminar boundary layer flow [65]. It was found the CFL number required by the FV-LBE was smaller than that of the DUGKS, and the computational efficiency was also lower due to the smaller time step.

Comparison with the UGKS for multiscale flows
The DUGKS can be viewed as a simplified version of UGKS, by taking the advantages of the LBE in the reconstruction of the interface flux. The collision and free streaming processes are taken into account simultaneously in the flux reconstruction in both methods, making them share the UP property for simulating multiscale flows.
The main difference between the two schemes lies in the way to obtain the distribution function at cell interfaces f x jk , t . In UGKS, the distribution function takes the formal integral solution of the kinetic equation [18], It is apparent that this solution consists of a hydrodynamic part and a kinetic part (the first and second terms on the right hand side, respectively). The hydrodynamic part includes the contribution from collisions, and the kinetic part is the transport of the distribution function at time t n , which reflects the contribution from free transport. These two parts contribute to the distribution function dynamically according to the ratio of relaxation time to the numerical time step. With proper approximation can be reconstructed based on Eq. (38). Then the micro flux can be obtained, wheref = 1 t t n+1 t n f (t)dt is the time-averaged interface distribution function.
Compared with the micro flux defined by Eq. (7) together with Eq. (14) in the DUGKS, it can be found that the calculation of UGKS is a bit more complex. We now make an analysis of the difference between the reconstructed time-averaged interface distribution functions in the UGKS and DUGKS. Specifically, if physical quantities are assumed to be linearly distributed around the cell, we can obtain the structure off x jk , ξ in UGKS based on the results in [70] as follows, where the terms on the right hand side are evaluated at x jk , ξ , t n , D = ∂ t + ξ · ∇, and a 0 and a 1 are defined as with α = t/τ . On the other hand, the same structure of time-averaged interface distribution function for the DUGKS can be obtained based on Eq. (8), with As α 1 (i.e., free molecular regime), it can be shown that It is clear that the four coefficients are identical in the leading order, and the difference between a 1 and a 1 is negligible as t/τ 1. On the other hand, in the limit of α 1 (i.e., continuum regime), it can be shown that It seems that the differences between the corresponding parameters are large in this limit. However, in this case the distribution function can be approximated by the first-order Chapman-Enskog expansion at the Navier-Stokes level, i.e., With this approximation, the last two terms on the right hand side of Eq. (40) are of order O(τ 3 / t) and O τ 2 , respectively; while these two terms of Eq. (42) are of order O τ 2 and O(τ t), respectively. Although different, these terms are high-order terms and do not appear at the Navier-Stokes level. Actually, in [70] and [19], it has been shown that bothf ugks andf dugks share the same formulation for continuum flows, namely, In the intermediate regime where α ≈ 1, it is easy to estimate that a 0 ≈ 0.632 and a 1 ≈ 0.264, while a 0 ≈ 0.6 and a 1 ≈ 0.3. The differences between a 0 and a 0 and between a 1 and a 1 are not large, either. Therefore, it can be expected that both the UGKS and DUGKS will behave similarly in this regime. Based on the above analysis, we can conclude that both the UGKS and DUGKS share the same behaviors in the whole flow regimes. This is also confirmed by a number of available numerical results. For instance, the shock structure of argon gas at different Mach numbers has been simulated by both the DUGKS and UGKS, and the results (density, temperature, shear stress, and heat flux) are nearly indistinguishable [33]. The Sod's shock tube problem was also tested by both methods from continuum to free molecular regimes, and again nearly identical results were obtained [33]. Wang simulated the lid-driven cavity flow at Kn = 0.1, and the flow structures predicted by both methods were in excellent agreement [71]. Some comparisons of the DUGKS and UGKS for flows of binary gas mixtures in different flow regimes were also provided in [40], including the shock structure, channel flow, and cavity flow. The results confirm again the similar behaviors of both methods in all flow regimes. The comparison of computational costs for the one-dimensional shock structure problem shows that the DUGKS is approximately 10% ∼ 20% faster than the UGKS [71], which is consistent with the previous estimation.

Turbulent flows in a periodic box
The DUGKS has been applied to turbulent flows as a direct numerical simulation method. Wang et al first studied the decaying homogeneous isotropic turbulence (DHIT) in a periodic cubic box [26]. In the study, simulations at different Taylor Reynold numbers Re λ were performed. Some quantities characterizing the turbulent statistics were compared with those of the standard LBE and pseudo-spectral (PS) methods. It is shown that the instantaneous vorticity field and the spectra of the kinetic energy and dissipation rate at Re λ = 26.06 predicted by the DUGKS with a uniform mesh of size 256 3 agreed well with those of the PS and LBE methods. The velocity derivative skewness and flatness were also compared at different Reynolds numbers and mesh resolutions. The results showed that the DUGKS could adequately resolve the flow when the minimum spatial resolution parameter k max η > 3, where k max is the maximum resolved wave number and η is the flow Kolmogorov length. This resolution requirement can be compared with the requirements of k max η > 1 for the pseudo-spectral method and k max η > 2 for the LBE. This indicates that the DUGKS has a relatively large numerical dissipation compared with the LBE method, although both are of second-order accuracy in time and space. The difference can be attributed to the finite-volume formulation of the DUGKS, while LBE is a Lagrangian scheme.
The DUGKS was also employed to simulate the Kida vortex flow with a relatively low Reynolds number [26] with the minimum k max η = 3.29. The results showed that the DUGKS could accurately predict the low-order statistics (total energy, dissipation rate, enstrophy, and velocity correlations), and could capture the smallest Kolmogorov length scale. The results of the pressure-velocity correlation also demonstrated that the DUGKS could well reproduce the incompressibility behavior of the flow.
Another turbulent flow simulated by the DUGKS is the Taylor-Green vortex in a periodic box [27]. The statistical properties, including the kinetic energy, dissipation rate, skewness, and flatness, were compared with the results of PS and LBE methods as well as a short-time analytical solution. The results confirmed the dissipation property of the DUGKS as observed in [26].
The simulation results of the decaying turbulent flows in a periodic box suggest that the DUGKS is less accurate than LBE, which seems to be contradictory to the findings in the previous comparisons [59,66]. The cause may be related to the different realizations of the no-slip boundary condition on solid walls in the LBE and DUGKS. In LBE the bounceback rule is implemented at grid nodes which are not located at walls, while in DUGKS the rule is implemented at the cell surfaces located on walls. For flows with periodic boundaries, there are no influences from boundary conditions, and thus reflect more clearly the intrinsic dissipation nature of the methods. If this argument is sound, the DUGKS will be more effective for flows involving solid walls, particularly with nonuniform meshes clustered near boundaries.

Turbulent channel flows
The first attempt to apply the DUGKS to wall-bounded turbulent flows was attributed to Bo et al. [27], in which a turbulent channel flow at Re τ = 180 was simulated with a nonuniform mesh with minimum mesh spacing y + = 0.425 in wall unit, which is sufficient to resolve the wall layer. For comparison, the standard LBE was also applied to this flow, but with a smaller computational domain due to the expensive computational costs with a uniform mesh.
The time-averaged velocity and turbulent Reynolds stress predicted by the DUGKS were compared with the PS and LBE results. It was observed that the DUGKS results were in good agreement with the predictions of the PS and LBE methods. Particularly, the velocity profile fits well with the standard linear viscous sublayer scaling for y + < 5 and the inertial sublayer scaling starting at y + > 30. The averaged Reynolds stresses predicted by the three methods were also in good agreement. However, it was reported that for the statistic stationary mean flow speed (averaged over the whole flow domain), the DUGKS gave a better prediction than LBE in comparison with the PS result, with relative errors 0.5% and 1.1%, respectively. The root-mean-square (rms) fluctuation velocities were also measured. Overall, the DUGKS gave better predictions than the LBE, particularly in the near-wall region. It should be noted that in the simulations the grid resolution used in the DUGKS was rather coarse in comparison with the PS method, and the grid cell aspect ratios near the wall were quite large. Actually, the maximum aspect ratio of streamwise to transverse grid spacing near the wall reached 41.62. Even with this grid ratio, the DUGKS was still numerically stable. It was also pointed out that the LBE with a uniform cubic grid could not handle the same domain size with limited computational resources.
Finally, it was shown that although the time step in DUGKS was small due to the use of non-uniform grid, the reduced number of grid points made the overall computational cost of DUGKS comparable to LBE [27]. The results also suggested while the DUGKS was used as a direct numerical simulation (DNS) tool since the local grid spacing was sufficiently small to adequately resolve all scales of the flow, it had the potential to be used as an implicit large-eddy simulation (LES) tool for high-Reynolds number flows at a given grid resolution due to its preferred numerical stability.

Turbulent natural convection
Besides isothermal turbulent problems, the DUGKS was also applied to turbulent flows with heat transfer. Wang et al. studied the natural convection of air in a three-dimensional (3D) cavity using a simplified version of the DUGKS under the Boussinesq assumption [72], which is constructed based on two weakly coupled reduced kinetic equations for the velocity and temperature evolutions, respectively [73].
In the study, the flow characteristics were analyzed as the Rayleigh number (Ra) ranged from 10 3 to 10 10 . From the instantaneous velocity and temperature patterns on the symmetry planes for 10 7 Ra 10 10 , it was found that both the temperature and velocity boundary layers near the two isothermal side walls became thinner with increasing Ra, while no apparent boundary layers were developed near the adiabatic walls. Two flow tunnels appeared near the center of the isothermal walls, and moved to the corners in the joint of isothermal and adiabatic walls with increasing Ra. In addition, as Ra increased above a critical value, the flow developed from steady to unsteady states and eventually turned to be turbulent as Ra ≈ 10 9 . It was also found that the lateral adiabatic walls had an inhibition effect on the temperature field, and the distribution of the local Nusselt number on the hot wall suggested that the lateral adiabatic walls could suppress the heat transfer. Generally, the convective intensity in the cavity center was observed to be stronger than that close to the adiabatic walls, while the effect of adiabatic walls on the overall heat transfer decreases with increasing Ra.
The time-averaged characteristic quantities of interest on the symmetry plane for Ra up to 10 10 were measured, and a scaling law between the Nusselt numbers (local and overall) and the Rayleigh number were identified, i.e., Nu m = 0.1522Ra 0.2942 , 10 3 Ra < 10 7 , 0.3533Ra 0.2395 , 10 7 Ra 10 10 , and Nu o = 0.1270Ra 0.3052 , 10 3 Ra < 10 7 , 0.3408Ra 0.2410 , 10 7 Ra 10 10 , where Nu m and Nu o are the mean Nusselt numbers on the symmetry plane and the overall one, respectively. These correlations were shown to agree well with the available numerical and experimental data at several specific values of Ra.

Compressible flows
Besides nearly incompressible flows with low Mach numbers, the DUGKS has also been used to simulate a number of high Mach number (Ma) compressible flows. For instance, the one-dimensional shock structure of argon gas was tested by the DUGKS based on the Shakhov model under different Mach numbers (1.2 Ma 8.0) [33]. It was found that the density, temperature, shear stress, and heat flux were nearly indistinguishable from the data of UGKS for all cases, and both were consistent with the Boltzmann solutions or DSMC results for Ma = 1.2. For cases of Ma = 3 and 8, the predicted density and stress still agreed well with the Boltzmann/DSMC solutions, but some discrepancies in temperature and heat flux appeared in the upstream, although they were still in good agreement with the benchmark data in the downstream. The inconsistency could be attributed to the use of single relaxation time in the Shakhov model. The simulation results also demonstrated that the DUGKS can serve as an effective shock-capturing scheme as the numerical cell size is much larger than mean-free-path.
The DUGKS was also employed to simulate the standard shock tube problem under different mean-free-paths. It was shown that the DUGKS could give satisfactory predictions from continuum to free molecular flow regimes. This property was also confirmed by simulating the 2D Riemann problem by comparing the numerical results with the Euler solution in the continuum limit and the analytical solution of the collision-less Boltzmann equation in the free molecular limit [33].
The supersonic flow over a circular cylinder at Ma = 5 was simulated via the DUGKS with an unstructured mesh to investigate the influence of Knudsen number [48]. At Kn = 0.1, it was found the DUGKS results, including heat flux, pressure, and shear stress along the cylinder surface, all agreed well with those of the DSMC results overall. But some small discrepancies in the temperature field in front of the bow shock were observed, which was attributed again to the difference between the Shakhov model and DSMC. For the case of Kn = 1.0, the discrepancies in the front shock were more obvious due to the increase of non-equilibrium effects. However, the heat flux, normal pressure, and shear stress along the cylinder surface were still quite satisfactory in comparison with the DSMC data, which are critical quantities of interest in engineering applications.
The flow over two side-by-side identical circular cylinders with a gap at Ma = 2 and Kn = 0.1 was further simulated. Again the temperature and Mach number distributions were found to be nearly identical to the DSMC results, and it was observed that the gas could be chocked in the gap between the cylinders. The high temperature region in front of the cylinders covered the upstream of the gap, and the temperature dropped gradually at the downstream of the gap, where the gas was accelerated to supersonic speed. The heat flux, pressure, and shear stress on the cylinder were measured to quantify the results, and the results were in good agreement with the DSMC data.
A more challenging multiscale compressible flow was investigated using the DUGKS in [48]. In this problem, two cavities (A and B) filled with gas at different pressures were connected by a thin tube. Initially, a diaphragm was put at the middle of the channel, and the gas temperatures in the two cavities were set to be 273K. The initial pressures in cavities A and B were P A = 48.78Pa and P B = 0.004878Pa, respectively, such that the corresponding Knudsen numbers were Kn A = 0.001 and Kn B = 10, respectively. At time t = 0, the diaphragm was removed suddenly, and the gas then started to expand from cavity A to cavity B. The dynamic behavior of the gas during the expansion process was then measured. It was found that the shock wave developed at an early stage and reached the center of cavity B. At this stage, the gas in cavity B was still very rarefied and underwent ballistically. The pressure in cavity B increased continuously with time as the gas moved in. The pressure ratio between the two cavities was high enough to form a supersonic jet at the outlet of the channel. At a later time, the initial shock wave disappeared and two symmetric vortexes appeared in cavity B. The detailed distributions of temperature, velocity and pressure during the expansion process showed that the shock wave was strong in the early stage, and became weaker gradually with time. It was also observed that the flow in cavity A changed slightly in the time interval, and the temperature and pressure were nearly uniform, but the temperature decreased with time as the international energy converted into kinetic energy during the expansion process.
Finally, we note that some compressible turbulent flows were investigated using the DUGKS by Wang and his colleagues [34,35]. A 5th-order weighted essentially nonoscillatory (WENO) scheme was incorporated into the DUGKS to better reconstruct the distribution functions at cell interfaces, although the overall accuracy is still second order. The new DUGKS was tested by simulating the compressible DHIT problem with low and high initial turbulent Mach numbers and Taylor Reynolds numbers. The turbulence properties, such as turbulent kinetic energy, rms values of the fluctuations, and the probability density function (PDF) of the local Mach number, were compared with those of other high-order DNS results, and good agreement was observed. Particularly, the distribution of the local Knudsen number was measured as a typical feature of compressible turbulence, and some non-continuum regions were identified in the flow field. Via the DUGKS, the influence of bulk viscosity on turbulence statistics and flow structures was also numerically studied. Overall, the results show that the DUGKS scheme can serve as a reliable tool for simulating compressible turbulence at low and moderate turbulent Mach numbers.

Micro flows
For micro flows of gases, a typical feature is the flow speed is usually small and the Knudsen number may vary in a wide range. Consequently, the DSMC method needs a long time averaging to lower the statistical noise, and it becomes computational expensive as Kn is small due to constraints on time step and mesh size. The classical DOM is free of statistical noise, but still suffers from expensive computational costs due to the time-splitting treatment of the collision and streaming processes as in DSMC. It is noted that some improved DSMC have been developed to reduce the noise, e.g., [74]. Particularly, a particle version of the DUGKS was proposed recently which can improve the performance of particle simulation greatly [75].
Unlike the DSMC and classical DOM, the DUGKS provides an efficient tool for the simulation of micro flows with the nice UP properties. Actually, a number of such systems have been investigated numerically by the DUGKS. The applications can be classified into two categories according to the devices. The first type is for devices with moving parts, and the second type is for those without any moving parts.

Flows in devices with moving parts
For problems of the first type, the micro Couette flow between two moving planar plates and the micro lid-driven cavity flow have been investigated under isothermal condition in the first paper of the DUGKS series [19]. It was shown the nonlinear velocity near the walls could be successfully captured in the Couette flow at finite Knudsen numbers, and the shear stress was in excellent agreement with the results of DSMC and Boltzmann equation in whole flow regimes.
For the cavity flow, a typical feature is that the distribution function can be highly irregular with discontinuities induced by the walls, particularly around the corners, and the deviation from the local equilibrium increases with the Knudsen number. In the simulations, it was found that the DUGKS could yield satisfactory predictions that agreed excellently with the DSMC solutions as Kn changes from 0.1 to 8, and the slip velocity on the walls approached a finite value at each wall. However, vortex center was found to be close to the vertical center line of the cavity in all cases, and moved downward with increasing Kn. The micro cavity flow with the consideration of thermal effects was investigated later by the DUGKS based on the Shakhov model [48]. The flow characteristics ranging from continuum to free molecular regimes were analyzed and compared to the DSMC or Navier-Stokes results. It was observed the heat flux was not consistent with the Fourier law for noncontinuum flows, even as Kn = 0.075. This observation is interesting, since generally it is believed that the Navier-Stokes-Fourier equations can still be used for simulating micro flows provided suitable slip and temperature jump boundary conditions are imposed. This test case clearly demonstrates this is not necessarily the case in some problems and we must be careful to use such an approach in the study of micro flows. A 3D micro cavity flow was also simulated with the DUGKS and compared with the DSMC results [76], which confirmed again the advantages of the DUGKS in simulating low speed micro flows.
Some unsteady micro flows in devices with moving parts were also investigated by the DUGKS. For instance, Wang and his colleagues conducted a number of studies on oscillating flows. In [77], the unsteady flow in a 2D rectangular cavity with an oscillating top lid was numerically investigated. To demonstrate the general flow characteristics of this problem, a number of simulations were performed under different conditions. Specifically, a wide range of parameters, including the aspect ratio of cavity width to length (A), the oscillation Mach number, the Strouhal number (St), and the Knudsen number, were considered. It was found that the flow properties, including the flow velocity, temperature, shear stress, and heat flux, were symmetrical about the vertical centerline of the cavity for Ma = 0.01 as Kn changed from 0.001 to 10. Actually, in this case the magnitude of the oscillation was so small that the system was nearly linear. On the contrary, these flow properties became asymmetric at Ma = 1.2 when the system was fully nonlinear, and the strength of shear stress near the top-right corner of the cavity could be much larger than that at the top-left corner, while the temperature at the top-right corner could be significantly higher than the wall temperature. An interesting phenomenon was reported that the heat transfer behavior could be altered by the oscillation. Previously, it was shown that heat could be transferred from the cold to hot regions for low speed cavity flow with constant driven velocity in non-continuum regimes [48]. However, as the oscillation frequency was above a certain value, it was found that heat was still transferred from hot to cold regions for both Ma = 0.01 and 1.2, although the heat flux could be non-parallel to the temperature gradient. The damping force exerted on the oscillating lid was found to depend on the Strouhal number non-monotonously generally, which was attributed to the anti-resonance and resonance of rarefied gas flows, respectively. Then scaling laws for the anti-resonance frequency and the inverse aspect ratio of the cavity were established for cases from near continuum to highly rarefied regimes.
The study of 2D oscillating cavity flow was later extended to 3D case with a small Mach number [78], focusing on the confinement effects of the additional lateral walls on the damping force on the lid. It was found the damping force in a 3D cavity could even be smaller than that in a 2D one for highly rarefied flows with a low oscillation frequency. Furthermore, the damping force was found to increase with decreasing depth of the cavity due to the effect of gas anti-resonance as the frequency was high for highly rarefied flow. It was also shown that the gas resonance and anti-resonance phenomena appeared in 2D cavity also occurred in 3D case, but the presence of the lateral walls tended to suppress their formation. Similar to the 2D cavity case, linear scaling laws for the anti-resonance and resonance frequencies were also obtained for flows from near continuum to free molecular regimes, which suggested that reducing the Knudsen number and increasing the aspect ratio of the cavity could enhance the 3D effects on formation of gas resonance and anti-resonance.
Further study on the effects of oscillation on heat transfer in a 2D square cavity was conducted [79]. It was found that the thermal convection could be dramatically enhanced under oscillation conditions even at moderate Knudsen numbers, which played a dominant role in the heat transfer. The average Nusselt number on the lid was found to decrease with increasing Kn, and could vary non-monotonically with the oscillation frequency, with the maximum occurring at nearly the anti-resonance frequency.

Flows in devices without moving parts
For micro flows in devices without moving parts, Zhu and his colleagues studied a number of micro flows driven by temperature difference [29,30,32]. Under rarefied conditions, gas flows can be induced by nonuniform temperature fields or temperature changes at solid boundaries. Generally, different types of such flows can be identified, e.g., thermal creep (or thermal transpiration) flow in a channel as the temperature gradient is applied along the channel walls, radiometric flow around a solid plate with different surface temperatures immersed in a vacuum, flow between a heated cantilever next to a cold substrate inside a vacuum enclosure, and flow induced by temperature discontinuity. Although the flow velocity of such flows is usually very low, they have potential applications such as Knudsen pump and mesosphere flight vehicle. Zhu and Guo [32] explored these four types of thermally induced flows in 2D geometries for a wide range of Knudsen numbers. Specifically, the thermal creep flow in a closed channel at a length-towidth ratio of 5 was simulated, where the temperatures at two ends were maintained at two different constants while the temperatures on the two lateral side walls were both linearly distributed between the two end ones. It was observed that the flow patterns could be quite different at different Knudsen numbers (0.01 Kn 10).
The flow induced by a hot micro beam with square cross section immersed in a square cavity was further numerically studied [32]. It is known that a flow would develop around a heated object as the separation between it and a cold substrate is comparable to the mean free path of the gas, and hence a net force (Knudsen force) will be exerted on the object. Therefore, as the micro beam and the outer cavity wall serve as the heated object and the cold substrate, respectively, a flow will also be induced. Such flow is fundamentally different from the thermal creep flow. It is actually caused by the combined effects of thermal stress flow and thermal edge flow due to the inhomogeneity of the temperature gradient and the sharp corners of the micro beam. Three cases, i.e., Kn = 0.1, 1, and 10, were considered in the study. It was found that some complicated vortexes were developed at the corners of the micro beam in each case, and the non-uniformities of both the normal stress and the shear stress along the micro beam surface tended to be more obvious as the Knudsen number decreased, and the stresses were accumulated near the corners of the micro beam.
Another thermally induced flow considered in [32] was the radiometric flow generated by a small fixed plate with differentially heated sides placed in a chamber with cold walls. Two Knudsen numbers were considered, i.e., Kn = 0.1 and 1. It was found four vortexes appeared at the corners of the plate in each case, and the strength and sizes of the vortex were nearly the same at Kn = 1. But at Kn = 0.1, the sizes of the two vortexes near the cold surface of the plate were smaller than the two ones near the hot surface. The pressure (normal stress) difference between the hot and cold surfaces of the plate along the vertical direction, which is the main contribution to the radiometric force, was also measured. It was shown that the pressure difference was nearly uniform along the plate surface in the vertical direction at Kn = 1, but was nonuniform and became larger near the top and bottom surfaces of the plate at Kn = 0.1, and the smallest pressure difference was still larger than that in the case of Kn = 1.
The micro flow in a fixed square cavity induced by temperature discontinuities was also investigated in [32]. In this problem, the top wall maintained at a higher temperature (T h ), and the other walls maintained at a lower one (T c ). The temperature discontinuities at the two upper corners could induce a flow in the cavity. The flow and heat transfer behaviors from continuum to free molecular regimes were analyzed. It was observed that in the continuum limit (Kn = 0.001), the velocity was rather weak and no vortex appeared, and the temperature could be well described by the heat conduction equation as (T h − T c )/T c 1, which has an analytical solution. For the slip (Kn = 0.1) and transitional (Kn = 1) cases, two pairs of vortexes appeared symmetrically along the vertical center line of the cavity. As Kn = 10, a number of secondary vortexes appeared and the flow pattern became more complicated. Some of the DUGKS results for the four cases were confirmed by DSMC [32]. It was shown that generally the predictions of both methods were in good agreement. But some discrepancies could be observed for large Kn, which could be attributed to the differences between the relaxation model employed in the DUGKS and the collision model in DSMC. However, some statistical noises were observed in the DSMC results as expected.
Recently, the DUGKS, together with the DSMC, was employed to study the flow and heat transfer behaviors in three categories of radiometric pumps consisting of channels with bottom or top surfaces periodically patterned with different structures [80]. The dominant mechanism in the radiometric force production was analyzed and operational performance of the devices was evaluated based on the numerical results. The analysis showed that the zigzag channel could generate maximum velocity with a parabolic profile, but the net radiometric force was rather weak in this case. For other configurations, the flow exhibited a linear profile in the open section of the channel. It was revealed that the radiometric force was produced due to the difference in particle momentum on both sides of the fins, and the mechanism was different from the standard Crookes radiometer due to the presence of the bottom surfaces. The mass flux was found to be a nonlinear function of the temperature difference, and its dependence on Kn was mainly caused by the structures of temperature field under different configurations.

Flows in devices with moving parts and temperature differences
Micro flows in devices with both moving parts and temperature differences were also studied, e.g., [30]. In such problems the overall flow was generated by both the forced motion of the moving part and the inhomogeneity of temperature. In [30], the nonequilibrium gas flow in a two-dimensional microchannel with a fixed hot ratchet surface and a moving cold wall was investigated via the DUGKS. With the periodic asymmetrical ratchet structures on the bottom wall and the temperature difference between the channel walls, a flow could be induced and a tangential propelling force could be generated on the wall.
Such thermally induced propelling mechanism could be utilized as a model heat engine. In the study, the flow field and propelling force were measured under different wall velocities and Knudsen numbers. Particularly, the flow fields at the critical wall velocity, at which the thermally induced force just balanced the drag force due to the active motion of the top wall, were analyzed. It was found the force changed linearly with the wall velocity, and the forces on the static wall and the top wall velocity at shear-free state achieved their maximum values as 0.1 ≤ Kn ≤ 1. The magnitude of the force was linearly dependent on the wall temperature difference, while the relative height of the ratchet structures affected the thermal driving effect significantly. A counter-intuitive relation between the flow direction and the shear force was also observed in highly rarefied conditions.
The output power and thermal efficiency of the system working as a model heat engine were also analyzed based on the momentum and energy transfer between the walls. The effects of Knudsen number, temperature difference, and geometric configurations were investigated. It was shown that the power output and the thermal efficiency achieved the maximum values in the early transition regime and were significantly larger than those in the free molecular regime. These results were helpful for improving the mechanical performance of the device.
In summary, the available results show that the DUGKS can be faithfully used for lowspeed micro flows, and the deterministic nature of the DUGKS makes it much more efficient than DSMC for such systems.

Two-phase flows
The DUGKS has been extended to isothermal two-phase flows based on some phase-filed models. In such cases, the DUGKS can be viewed as a special finite-volume counterpart of the LBE method. In phase-field theory for a binary fluid system, the thermodynamic behavior is described by a free-energy function related to an order parameter φ and its spatial derivatives. The order parameter is used to distinguish the fluid phase and varies continuously from one phase to the other, and the thickness of the fluid interface is finite. A widely used model in phase-field theory is the Landau free-energy function defined as where ψ(φ) is the bulk free-energy density, which usually takes a double-well formulation; κ is a parameter related to surface tension, and V is the system volume. When coupled with a flow field, the evolution of order parameter can be described by certain convectiondiffusion equations, such as the Cahn-Hilliard equation and the Allen-Cahn equation. The first DUGKS for immiscible two-phase flows was developed by Zhang et al. based on a quasi-incompressible Navier-Stokes equations and the Cahn-Hilliard equation [36], in which the governing equations read ∇ · u = −γ ∇ · λ∇μ φ , ρ ∂u ∂t + u · ∇u = −∇p + ∇ · ρν ∇u + (∇u) T + F, where the parameter γ is related to the density ratio of the two fluids, λ is the mobility for the order parameter, and F s = −φ∇μ φ is the interfacial force with μ φ = δF /δφ being the chemical potential. The DUGKS was constructed based on two discrete velocity kinetic equations,  (52) following the DUGKS framework, and simulated a number of two-phase fluid systems. For the problem of a stationary droplet immersed in another fluid, it was demonstrated that the DUGKS could predict the Laplace law correctly, and showed better conservation property than the corresponding LBE. A layered two-fluid flow in a channel with large viscosity ratio (up to 10 3 ) was also investigated, and the velocity profiles were well predicted. Particularly, by taking the advantage of nonuniform mesh, the numerical accuracy near the fluid interface was much improved. The test of a rising bubble under gravity showed that the shapes and positions of the bubble in the time history were well captured at different density ratios, and the numerical stability of DUGKS at large density ratio seemed to be better than the corresponding LBE. The Rayleigh-Taylor instability phenomenon of two fluids under gravity was further simulated. The results at different Atwood and Reynolds numbers showed that the positions of both bubble front and spike tip, as well as the interface growth amplitude, agreed well quantitatively with the previous numerical results and existing experimental data.
Chen et al. later employed the above DUGKS to study [37] the 3D Rayleigh-Taylor instability and two-phase homogeneous isotropic decaying turbulence. The results from the DUGKS were compared to those of the corresponding LBE and/or ARCHER code based on the Coupled Level Set-Volume of Fluid (CLSVOF) method [81]. The results demonstrated the reliability of DUGKS. For the Rayleigh-Taylor instability problem, a quantitative comparison was performed by tracking the position of the bubble, spike, and saddle points in the time history. It was found that the bubble front and saddle locations predicted by both DUGKS and ARCHER were in excellent agreement. For the spike, the results from both DUGKS and ARCHER were also in excellent agreement at early times. While at later times, the DUGKS predicted a lower spike location, which could be attributed to the numerical dissipation in both methods. For the two-phase decaying homogeneous isotropic turbulence, a droplet was initially positioned in a turbulent flow field. A key issue of this problem is to set up a consistent initial velocity field across the droplet interface. In [37], a forced homogeneous isotropic turbulence was first run using the ARCHER code to create a developed single phase turbulent velocity field. Then a solid particle, which was treated by an immersed boundary method, was put into the flow field. After several large eddy turnover times, the velocity field was used to initialize the two-phase flow by replacing the solid particle by the droplet with the same size. The time evolutions of the velocity and vorticity fields predicted by both DUGKS and ARCHER were compared, and it was observed that the complexity of the interface topology increased with time, and the results of both methods were in excellent agreement.
To further improve the capability in simulating two fluid flows with large density ratio, a DUGKS based on the incompressible Navier-Stokes equations coupled with a conservative Allen-Cahn phase-field model was developed in [82]. Two discrete-velocity kinetic equations in the formulation of (52) were again adopted as the starting point, but with different definitions of equilibrium distribution functions and source terms. The method was then employed to simulate a number of stationary and dynamic problems with density ratio up to 1000, and reliable solutions were obtained. But in comparison with the corresponding LBE, it was found that the numerical dissipation of DUGKS was relatively large and could destroy the interface structures, which suggested that a high-order DUGKS is preferred for two-phase systems involving complex interface changes.
The DUGKS method was also applied to solid-liquid phase change problems [83]. The flow field is governed by the Navier-Stokes equations, and the energy field is described by an equation for the enthalpy, where H is the total enthalpy. Two discrete velocity kinetic equations with the formulation of (52) were designed in which the equilibrium distribution functions and source terms were chosen such that the governing equations could be recovered in the Champan-Enskog analysis. Then the two kinetic equations were discretized following the DUGKS procedure. A number of phase change problems, including the two-region phase change, phase change under constant heat flux and natural convection with phase change, were simulated. The results predicted by the DUGKS were compared with benchmark solutions and satisfactory agreement was observed. The dynamic movement of the solidliquid interface during phase change processes was accurately captured. Particularly, for the natural convection problem, it was found that the temperature field near the interface could be predicted by the DUGKS and was more accurate than LBE.

Gas-solid flows
Flows with solid particles are another type of two-phase system. A number of direct numerical simulation methods for continuum particulate systems with resolved particle shapes have been developed, such as LBE, immersed-boundary method (IBM) [84], and fictitious domain method (FDM) [85]. Recently, the DUGKS was applied to particulate flows by treating the solid body with the immersed-boundary (IB) technique [86][87][88].
In the IB-DUGKS method for isothermal particulate flows [86], the entire domain, including the interior of solid particles, is assumed to be filled with fluid, and the surface of a particle is discretized into a set of Lagrangian points X l (l = 1, 2, · · · , N L ). The interaction between fluid and solid body is realized by adding an immersed boundary force to the fluid, which is distributed from the Lagrangian points to the Eulerian mesh with cell centers x j (j = 1, 2, · · · , N E ) employed by the DUGKS. The acceleration due to the IB force, a l at the Lagrangian point X l is determined by the velocity difference between the fluid and solid body at the point, where U l and u l are the desired velocity and the intermediate fluid velocity (without IB force) at X l , respectively. For a solid particle, U l = u p + ω p × X l − x p , with u p and ω p the translational and rotational velocities of the particle, and x p the position of particle center. The intermediate fluid velocity u l is interpolated from the fluid field using certain local weight functions such as the smoothed Dirac δ function. The calculated force a l is then distributed back to the Eulerian points using the same interpolation function to obtain the body force a for the fluid. This procedure can be iterated to ensure the no-slip boundary condition accurately [86].
Once the IB force is determined, the fluid can feel the existence of immersed boundary of the solid particle, and the BGK equation for the fluid can be written as where the forcing term is F = −a · ∇ ξ f , with a being the acceleration due to the IB force, can be simplified for continuum flows, In order to solve Eq. (55) with the DUGKS, the Strang splitting scheme is employed in the IB-DUGKS, where the evolution of the distribution function at time step t n follows three steps, (1) Pre-forcing step: Advance f n in each cell to obtain the first intermediate value, (2) DUGKS step without force: Advancef * following the standard DUGKS without forcing term to obtain the second intermediate valuef * * ; (3) Post-forcing step: Advance f * * to obtain the distribution function at next step, The dynamic of the solid particle in the IB-DUGKS follows the Newton's law, where M p and I p are the mass and moment of inertia of the particle, respectively, F p is the total force on the particle including the counter-acting IB force, and T p is the corresponding total torque. The accuracy of the IB-DUGKS was tested by simulating several 2D and 3D particulate flows, including the sedimentation of a particle and the drafting-kissing-tumbling (DKT) dynamics of two particles in a channel, and a group of particles settling in an enclosure. The numerical results predicted by the IB-DUGKS were found to be in good agreement with benchmark data.
Recently, the IB-DUGKS was further extended to systems with heat transfer between fluid and solid bodies with fixed temperature following similar idea [87,88], in which the energy equation was solved by another BGK equation for temperature distribution function. Specifically, a non-iterative technique was proposed to realize the exact velocity and temperature boundary conditions on the solid surface in the method reported in [87]. Both IB-DUGKS methods were tested by several thermal flows involving stationary solid bodies, but no results were reported for flow problems with moving bodies.

Gas-mixture flows
Multiscale flows of gas mixtures are widely encountered in many industrial and natural processes, and the transport of each gas species can be described by a Boltzmann equation with inter-molecular collisions between the same and different species. A number of simplified kinetic models have been proposed, among which the McCormack model [89] and the Andries-Aoki-Perthame (AAP) [90] model are widely used. The former is a model with a linearized collision operator under the assumption of slight deviation from equilibrium, while the AAP model uses a single BGK operator in which both self and cross-collision effects are incorporated.
Zhang et al. developed a DUGKS for multiscale binary mixture flows based on the following AAP model [40], where f α is the distribution function for species α, and the equilibrium distribution function is defined as where m α and ρ α are the molecular mass and density of species α, u * α and T * α are two parameters related to the hydrodynamic variables W α = (ρ α , ρ α u α , ρ α E α ) for each species, which are defined as the moments of f α similar to that for a single gas given by Eq. (3). Specifically, u * α and T * α are defined such that the total mass, momentum, and energy are conserved [90]. Generally, u * α (W ) and T * α (W ) are nonlinear functions and depend on interaction potential between molecules, which can be simplified for Maxwell molecules.
The DUGKS can then be constructed based on Eq. (60), similar to that for the single gas BGK Eq. (1). But a nonlinear system must be solved to determine u * α (W ) and T * α (W ) in the calculation of equilibrium distribution functions, which appear in both the update of cell-averaged variables and flux evaluation at cell interfaces.
The above AAP-DUGKS was applied to several 1D and 2D flows of binary mixtures with different mass ratios at different regimes, including the shock structure problem, the channel flows driven by a small pressure, temperature, or concentration gradient, the 2D plane Couette flow, and the cavity flow. Overall the DUGKS results agreed well with benchmark data obtained by the UGKS, DSMC, and/or the linearized Boltzmann equation. However, it was found that there existed some deviations for the light species as the difference in molecular masses was large, particularly at large Knudsen numbers. The discrepancies could be attributed to the limitations of the AAP model such as the incompatible transport coefficients.
In order to overcome the deficiency of the AAP-based DUGKS, a DUGKS based on the linear McCormack model was later developed [41]. In the McCormack model, it is assumed that the perturbations in concentration, pressure, and/or temperature are small (and thus the velocity is also small), then the system deviates slightly from the global equilibrium such that the collision term in the Boltzmann equation can be linearized [89]. The velocity moments of the collision term of the McCormack model match those of the Boltzmann equation up to third order, thus leading to the same transport coefficients as the latter.
The McCormack model can be expressed as where h α is the perturbation defined by f α = f Here ρ α0 and T 0 are the reference density of species α and reference temperature, respectively. The linearized collision operator can be expressed explicitly [41,89]. Based on Eq. (62), the DUGKS was constructed for solving the perturbation distribution functions h α following the same idea as the original DUGKS [40]. A number of low-speed binary mixture flows covering a wide range of flow regimes were simulated to test the DUGKS. Since the McCormack model can reproduce all transport coefficients, it is not surprising that the results of the McCormack-DUGKS were in better agreement with those of the linearized and full Boltzmann equations than the AAP-DUGKS, particularly in the transitional and near-continuum regimes for systems with large molecular mass difference. However, it should be noted that the McCormack-DUGKS can only be applied to the low speed flows that slightly deviate from equilibrium and is unsuitable for nonlinear problems. Therefore, there is still much room to improve DUGKS for more general multiscale mixture flows.

Strongly inhomogeneous flows
As a fluid is confined in a nanoscale structure, the solid-fluid interaction is significant and fluid properties may become strongly inhomogeneous. For such systems, the molecular size σ may be comparable to the device size L, and the confinement effect cannot be ignored. Therefore, besides the Knudsen number, the ratio σ/L is also a key dimensionless parameter, which is usually taken to be zero in the classical kinetic theories for homogeneous or weak inhomogeneous fluids. Consequently, the classical Boltzmann equation for ideal gases and Enskog equation for dense gases, are not adequate to describe such systems. Some modified kinetic models have been developed and employed to study static or dynamic behaviors of nano-confined fluids. For instance, a tractable kinetic model was proposed based on the local average density approximation and density-functional theory [91]. Based on this model, a DUGKS was developed for nano-confined strongly inhomogeneous fluid systems recently [92].
The tractable inhomogeneous kinetic model reads [91], where φ e is the external potential such as gravity or that exerted by the confined walls, φ m (n) is the mean field potential induced by other fluid molecules dependent on the local number density n, and J ex (n,n) is the excess part of the collision beyond the hard-core Boltzmann one, withn being a weighted average number density. The detailed expressions of J ex andn can be found in [91]. The relaxation time τ also depends on both n and n to account for the inhomogeneity, i.e., τ = μ(n)/nk B T, where μ(n) is the viscosity for a dense gas with number densityn derived from the Enskog theory for homogeneous fluids. The kinetic Eq. (64) can be regrouped as where G is the collection of the potential term and J ex . Then the DUGKS can be constructed to solve Eq. (65) by employing the Strang splitting technique as shown in Section 6.2. Note that significant density oscillations may occur in a strongly inhomogeneous fluid system, so the van Leer limiter was employed in the flux reconstruction during the DUGKS step without force.
The static structures and dynamic behaviors of dense fluids in slits were simulated by the DUGKS. The effects of fluid density, slit size, fluid-fluid and fluid-wall interactions, were investigated. The DUGKS results were in good agreement with those of the Monte Carlo and/or molecular dynamic simulations. Particularly, the layered structures due to the competition between fluid-fluid and fluid-wall interactions were successfully captured. Furthermore, in order to evaluate the local nonequilibrium effects, the local effective Knudsen number Kn e was measured, which is defined as where χ(n) is the radial distribution function evaluated with the average number density. It was found that Kn e could vary from 10 −4 to the order of 100 in the system. The simulation results also showed that the velocity profiles of the Couette flow and the Poiseuille flow deviated from the Navier-Stokes solutions clearly, as a result of the inhomogeneous nature of the dense fluids in nano confined geometries.

Multiscale transports beyond fluid flows
Besides gas molecules, other energy carriers such as phonons, photons, electrons, and plasma, can also undergo multiscale transports. The Boltzmann equations with different equilibrium statistics can also be employed to describe the transports of such particles [93]. Based on certain kinetic models, the DUGKS has been extended to such systems recently.
A few examples will be given briefly in this section.

Phonon heat transfer
A phonon is a quantum of lattice vibrational energy in crystalline solids and has particlelike properties [94]. Phonon transport is the main mechanism for heat transfer in semiconductors and dielectrics. Heat transfer process in nanostructures of such materials usually involves multiple temporal and spatial scales, characterized by the Knudsen number defined as the ratio of phonon mean free path and typical length of the material. The Boltzmann equation can be used to describe the phonon transport when the phase coherence effects are negligible, which can be expressed as where f p = f p (x, k p , s, t) = f p (x, ω p , s, t) is the distribution function for phonons of polarization branch p with wave number k p (or frequency ω p ) at position x and time t; v p = ∂ω p /∂k p = v p s is the group velocity with which the phonon travels along the unit direction s, and Q p represents the rate of change due to phonon scatterings. The scattering between phonons is quite complicated, and a widely used model is the BGK-like one, where f (eq) p (T * ) is the equilibrium distribution following the Bose-Einstein statistics, where is the Planck's constant divided by 2π, and T * is a local pseudo-temperature which is introduced to ensure the energy conservation of the scattering term, and may be different from the thermodynamic temperature T defined below. The effective relaxation time τ p (ω p , T) reflects the combination effects of all scattering processes and is usually estimated using the Matthiessen's rule if the individual scattering processes are independent of each other [94]. The mean free path of phonons is related to the relaxation time, i.e., λ = τ p v 0 , with v 0 being a typical value of the phonon group velocity. It should be noted that the mean free path depends on the frequency, polarization, and temperature, and can change over several orders of magnitude. Therefore, the phonon transport exhibits an intrinsic multiscale nature. The energy and heat flux carried by phonons are defined from the distribution function, where D p (ω) is the density of state and is the solid angle. The thermodynamic temperature T and pseudo-temperature T * are defined from the following constraints, It is noted that as τ p is frequency-and polarization-independent, the above two constraints are identical and T = T * . The temperature T is related to the energy as T = E/C v , with C v being the volume specific heat capacity. Even with the relaxation time approximation, the kinetic equation for phonon transport is still difficult to solve due to the dependence of frequency and polarization. Some further simplified models have been proposed to reduce the complicity [94]. For instance, the gray model assumes phonons of all polarizations and frequencies are same and the group speed v = v g s, with v g a constant. In this case we can introduce an energy distribution function independent of frequency and polarization, Then the total phonon energy E and heat flux q can be determined from e , The transport equation for the energy distribution function can be written as where τ is now a singlet frequency-independent relaxation time, but it may still depend on temperature. The equilibrium energy distribution e (eq) is given by which is just the angular average of the total energy. Apparently, the scattering term in the gray model defined above conserves the total energy. The gray model (74) has the same formulation as the BGK Eq. (1) for gas molecules, and a DUGKS could be constructed straightforwardly, which was reported in [42]. It was proved rigorously that the DUGKS is a Lax-Wendroff discretization of the kinetic equation without the scattering term, while becomes a finite-volume scheme for the diffusion equation in diffusive limit. This analysis suggests that the DUGKS has the unified preserving properties [12] and can serve as an ideal method for multiscale phonon transport problems. The property was also numerically verified by simulating several steady and unsteady heat transfer problems at different regimes. Particularly, a multiscale heat transfer across an inhomogeneous film with Kn changes from 10 −4 to 10 was investigated, and the multiscale transient heat transfer behaviors were successfully captured even with a coarse mesh and large time step. On the contrary, the classical DOM failed to capture the behaviors unless a very fine mesh was employed. Consequently, the DUGKS was much more computational efficient than DOM for this multiscale problem.
Although the gray model is simple and can provide some insightful predictions on the phonon transport behaviors, the energy exchange between different modes cannot be clearly identified in this model. Consequently, the DUGKS based on the gray model is limited to systems where the phonon relaxation mechanisms are unimportant. A DUGKS for phonon transport with the consideration of dispersion and polarization effects was later proposed [44]. The model was based on the assumption that the temperature difference in the system is small enough, i.e., T = |T − T 0 | T 0 (T 0 is the reference temperature). Under this assumption, the relaxation time and specific heat can be regarded as temperature-independent, and the thermodynamic temperature T and pseudo-temperature in f (eq) p can be regarded to be indifferentiable. Then the kinetic equation can be linearized and can be written in terms of the deviational distribution function where with C p (ω) = ω p D p (ω)∂f (eq) p /∂T being the mode specific heat capacity. The DUGKS was then constructed based on Eq. (76), and was verified by simulating several phononmediated heat conduction problems from ballistic to diffusive regimes [44].
Although the frequency and polarization effects are considered in the above improved DUGKS, it is limited to systems with small temperature difference and thus linear phonon transports. A DUGKS for multiscale heat transfer with arbitrary temperature difference was developed recently [43], which was based on the nonlinear relaxation model (67) together with (68). Different from the two DUGKS methods based on the gray model and the linear relaxation model in which only one temperature (T) is involved, two temperatures (T and T * ) appeared in the nonlinear DUGKS, which must be computed in the evolution. This was achieved by solving the nonlinear system Eq. (71) by the Newton iteration procedure. The numerical results of some steady and transient heat transfer problems showed that this nonlinear DUGKS could capture the multiscale phonon transport accurately for systems with both large and small temperature differences. Particularly, it was observed that the thermal transport phenomena with large temperature differences could behave quite differently from those with small ones, owing to the temperature-dependent relaxation time. Generally, the average temperature increased in ballistic regime but decreased in diffusive regime with the increase of temperature difference. Therefore, it is necessary to simultaneously consider the multiscale effects in both spatial and frequency spaces for problems with large temperature differences.
Overall, the DUGKS's based on different kinetic models for phonon transport provide efficient tools for predicting multiscale heat transfer problems. It is also noted that some special techniques were developed to accelerate the convergence for steady problems [95,96].

Radiation heat transfer
Radiative heat transfer caused by electromagnetic waves (or photons) is another type of transport phenomena, which can undergo multiscale behaviors with the change of optical thickness of the medium. The transport of photons can be described by the Boltzmann equation for the distribution function of photons f (x, s, ν, t), where s is the unit direction vector of photon propagation and ν is the photon frequency. However, radiative heat transfer is usually described in terms of the radiative intensity I = hνcf , where c is the light speed and h is the Planck constant. For gray media where the frequency dependence is ignored, the radiative heat transfer equation (RTE) can be written as where Q is the collision operator accounting for the interactions between photons and surrounding matter, which can be expressed as where β is the extinction coefficient that is the inverse of the local mean free path of photos, ω is the scattering albedo, I b is the blackbody intensity, and s , s is the scattering phase function which describes the fraction of the radiative energy scattered into the outgoing direction s from the incoming direction s, with being the corresponding solid angle.
For isotropic scattering problems, = 1; while for anisotropic radiation, depends on the scattering angle and can be approximated by a finite series of Legendre polynomials. For equilibrium radiation, the blackbody intensity I b is determined by energy conservation, i.e., while for nonequilibrium systems, it is given by the Stefan-Boltzmann law, where σ is the Stefan-Boltzmann constant and T is the local temperature of the medium. The radiative energy G and heat flux q are defined as the moments of I, The RTE (77) can also be written in other formulations. For example, in order to reduce the ray effect from boundaries, the RTE can be rewritten in a splitting formulation [97], where I = I c + I d , I c (x, s, t) is the direct intensity from boundaries and I d (x, s, t) is the diffuse one resulting from emission within the medium. The solution of I c can be given analytically with prescribed boundary conditions, and I d can be obtained numerically by solving Eq. (83b). Based on Eq. (83), Luo et al. developed the first DUGKS for isotropic radiative heat transfer problems [45]. Equation (83b), together with the analytical solution of I c , is solved with the procedure of DUGKS. The method was validated by simulating a number of transient and steady problems, including the transient radiative transfer in a plane-parallel slab filled with an absorbing-scattering medium, the equilibrium and nonequilibrium steady-state radiative heat transfer in a two-dimensional square medium, and a three-dimensional multiscale radiative transfer problem in a cube filled with an inhomogeneous absorbing-scattering medium. It was confirmed that the DUGKS exhibited the asymptotic preserving property such that it could give accurate solutions for systems with different optical thicknesses with relatively coarse meshes. Particularly, the DUGKS was found to be more computational efficient than conventional numerical methods in the optically thick regime. Furthermore, the results also demonstrated the good capability of the DUGKS in capturing the sharp spatial discontinuities in the radiation field and modeling multiscale radiative transfer.
The DUGKS developed in [45] is for isotropic media with = 1. Recently, a DUGKS for anisotropic radiative heat transfer was developed based on the RTE (77) directly [46]. A difficulty in this DUGKS arises in the relations between the original radiation intensity and the transformed ones due to the nonlinear anisotropic phase function . In isotropic case, the transformations are explicit, but for anisotropic media the relations are implicit. A simple iterative method was suggested in the calculations, which was shown to be able to give satisfied predictions [46]. But more efficient methods are desirable to improve the computational efficiency. The DUGKS was tested by simulating several 1D and 2D anisotropic radiative transfer problems with different optical thicknesses. The results confirmed the DUGKS exhibited the nice unified preserving properties and could serve as an efficient and accurate tool for radiative heat transfer in multiscale anisotropic media.

Plasma
A plasma is a partially or fully ionized gas containing electrons, ions, and neutral atoms or molecules [98]. The transport of different particles in a plasma can be modelled at different scales. Particularly, certain kinetic models of Boltzmann type have been developed to describe the evolution of a plasma at the kinetic scale, such as the Vlasov equation, the Bhatnagar-Gross-Krook equation, and the Fokker-Planck equation, coupled with the Poisson equation for electric potential, or more generally with the Maxwell equations to include magnetic and electromagnetic effects.
Like the usual gas systems, plasma dynamics is characterized by a wide range of spatial and temporal scales. Furthermore, in addition to the mean free path λ and collision time, velocities span a wide range. This reflects the common limitations of discrete velocity methods for problems involving a wide velocity space.
Very recently, Liu et al. further developed an improved DUGKS which has the asymptotic properties in terms of both Kn and D , based on the BGK-Vlasov equation and a reformulated Poisson equation (RPE) [100]. The RPE reads in dimensionless form as, which was then solved by a finite-element method. It was shown that this DUGKS based on the BGK-Vlasov and RPE equations could preserve the hydrodynamic and quasineutral limits as Kn → 0 and D → 0, which was also confirmed by a number of numerical simulations.

Summary and outlook
It is a challenging task for developing efficient and robust numerical methods for multiscale flows due to the involved large span of spatial and temporal scales. Kinetic schemes based on the Boltzmann or model equations have the potential to serve this purpose, but it is non-trivial to design a kinetic scheme which can capture the hydrodynamics without resolving the kinetic scale, i.e., exhibit the UP properties. The DUGKS is one such kinetic scheme with the desired properties, and its merit lies in the reconstruction of the numerical flux at cell interfaces, which is based on the numerical solution of the kinetic equation itself. This means that the spatial-temporal coupling dynamics is taken into account simultaneously, which is necessary for developing reliable numerical methods [101]. Furthermore, with the coupling of collision and free transport, the underlying physics of the kinetic equation is respected by the DUGKS, and thus it exhibits some nice properties, such as low numerical dissipation and robustness, that distinguish it from other kinetic schemes. It is interesting to note that the standard LBE could also be viewed as spatial-temporal coupling method if one tracks the transformed distribution functionf on a regular lattice [102], which follows the simple collision-streaming procedure, although it is generally considered as a decoupling scheme. The LBE method has shown many distinguish advantages in simulating complex flow problems, mainly lying in continuum regimes. However, some intrinsic limitations still exist in LBE, as pointed out in the recent perspective article [103]. Specifically, four inconveniences were identified: the use of uniform mesh, slow convergence for steady flows, limitation of Mach number, and difficulties in actual coding for the implementation of complex boundary conditions. In addition to these points, it is also a non-trivial task to design a LBE model for strong non-equilibrium flows which usually requires a large number of discrete velocities coupled with a regular lattice. Although some progresses on these subjects have been achieved separately, some special techniques are usually involved and the advantages of LBE will be lost more or less. On the other hand, the DUGKS does not suffer from these limitations and provides a unified numerical tool. Actually, an open source code is available to be used to simulate low and high Mach number flows in different regimes with unstructured meshes [76].
Indeed, after years of development the DUGKS has demonstrated its advantages in simulating multiscale flows, including turbulent flows, particulate flows, two-phase flows, flows of gas mixtures, and micro/nano flows. It has also been extended to study transport phenomena of other energy carriers, such as phonon, photon, and plasmon. The available results present a clear evidence that the DUGKS can serve as an appealing numerical tool for the calculation of multiscale transport problems. Despite the success of DUGKS, there is still much space for improvements. Specifically, the following critical problems should be considered to improve the performance of DUGKS.

Memory reduction
As a deterministic discrete velocity method, a large number of discrete velocities may be required, particularly for 3D highly compressible or highly rarefied flows, to capture the irregular distribution function. For such problems, the memory requirement and computational cost will be quite demanding. Therefore, developing memory reduction techniques is critical for DUGKS in simulating such flows. One natural choice is using the velocity adaption technique, which has been employed in some discrete velocity methods (e.g., [55,56]). In this method, fine and coarse velocity grids are used adaptively in different flow regions, depending on the regularity of the distribution function. However, conservation property should be carefully addressed during the particle transport on different velocity grids.
Another more elaborate memory-reduction technique, i.e., reduced-order-modeling (ROM), was introduced into DUGKS recently [104]. With the ROM, a reduced discrete velocity space can be selected from the original discrete velocity space to represent important dynamical characteristics. Consequently, a large number of grid points in the original discrete velocity space, which contribute little to the dynamics, can be removed in practical computations. The numerical tests showed that the memory in the ROM DUGKS could be significantly reduced and the computational efficiency could be greatly improved [104].
Beside the above two deterministic memory reduction techniques, stochastic particle methods can also be combined into discrete velocity methods. A successful example is the recent unified gas-kinetic wave-particle (UGKWP) method [23]. In this method, both the distribution function and simulated particles are used to describe flow physics, where the particles are sampled only for capturing the local non-equilibrium part caused by free transport. In each control volume, the evolution of kinetic equation is coupled with the macroscopic conservation equations. With the introduction of particles, the velocity space can be sampled adequately such that the memory requirement as well as computational cost could be reduced significantly. This technique could also be incorporated into the framework of DUGKS in principle, but the detailed implementation needs further investigation.
It is also noted that a very efficient memory reduction technique was developed for solving stationary kinetic model equations [105], in which velocity distribution function was reconstructed from macroscopic variables. Therefore, only the macroscopic quantities appearing in the collision term are stored, and the memory requirement for enormous discrete velocities is totally relieved. Furthermore, with the introduction of a prediction step to estimate the equilibrium state by solving the macroscopic governing equations, the efficiency and accuracy of this method were greatly improved in the near continuum and continuum regimes [106]. This technique may also be employed in DUGKS. However, it is still not clear how to apply this technique to transient kinetic models.
Finally, it is remarked that for continuum flow the memory can be greatly reduced by removing the evolution of discrete distribution functions, since in this case the distribution functions can be approximated with the Chapman-Enskog solution at the Navier-Stokes equations, which depends only on the equilibrium distribution function and its gradient. As such, only the conserved variables are required to evolve, and the memory requirement can be much reduced. This idea was originally developed in the lattice Boltzmann flux solver [107,108], and can be used to improve the memory and computational efficiency of DUGKS for continuum flows.

Implicit schemes
With the memory reduction technique, the computational efficiency could be improved. However, as an semi-implicit time-marching method, the time step of the standard DUGKS is limited by the CFL number, which may be rather small and thus controls the overall computational efficiency. In this regard, many techniques widely used in classical computational fluid dynamics (CFD) can be employed to release this restriction. For instance, an implicit DUGKS was designed for steady flows [109], in which the macroscopic equations for the conservative variables were solved iteratively as a prediction step for the iteration procedure of the kinetic equation. The delta formulations of the conservative variables and distribution function were employed in the iterations and the matrix free Lower-Upper Symmetric Gauss Seidel (LU-SGS) was used to solve the implicit equations. Furthermore, a pseudo time step with a large CFL number was adopted to accelerate the convergence. Some numerical tests of (nearly) incompressible and compressible flows under different Knudsen numbers were performed. The results showed that the computational efficiency could be improved by one or two orders of magnitude in comparison with the explicit DUGKS.
The implicit DUGKS developed in [109] was based on the transient kinetic equation, although a large pseudo time step could be used. Recently, a DUGKS was developed to solve directly the steady linear kinetic equation for neutron transport [110], where the distribution function at a cell interface is obtained by integrating the steady kinetic equation along the neutron transport direction. This idea could also be used to design implicit DUGKS for steady flow problems.
The above two schemes were designed for steady problems. More efficient DUGKS algorithms for unsteady problems are still desired. The strategies used in CFD, such as dual time stepping, multi-grid, and temporal adaption techniques, could be adopted to this end. We note particularly that an implicit UGKS for unsteady flows was developed recently [111], which solves the macroscopic conservative equations in delta form iteratively with a large numerical time step. The flux for macroscopic equations is obtained from a local time-averaged one given by the distribution function, which is calculated iteratively from the kinetic equation also in delta form. The cell size effect is incorporated into the local time averaged flux in this method, such that the numerical time step is not restricted by the CFL condition. It is expected that the techniques used in this implicit UGKS could also be employed to design implicit DUGKS for unsteady flows.

High-order schemes
The available DUGKS is of second-order accuracy in both space and time. In some cases, higher order accuracy would be necessary (e.g., [82]). Some efforts have been made to design high-order DUGKS's. For instance, a third-order DUGKS has been developed based on a two-stage time-stepping scheme and a third-order flux reconstruction [112]. Numerical results demonstrated that the scheme is of third order accuracy in both space and time. On the other hand, it was claimed that it was difficult to design a DUGKS of order higher than three in time if the trapezoidal rule is employed in determining the interface distribution function along the characteristic line.
It is noted that even in the second-order DUGKS, some partially high-order techniques could also be employed. For example, one can use high-order interpolations (such as weighted essentially non-oscillatory schemes) to construct the distribution function f + x at the starting point of the characteristic line in (9). As such, although the overall accuracy is still of second-order, the absolute error can be reduced.
Another appealing technique for improving the accuracy is the two-stage Lax-Wendroff time stepping method [101], which has been applied to the gas-kinetic scheme [113]. It is worthwhile to try this technique in developing high-order DUGKS methods in the future.
In summary, the DUGKS has gained much success in simulating multiscale flows and demonstrated great potentials in simulating other transport phenomena, but it is still far from maturity and needs further improvement in many aspects such as memory reduction, implicit discretization, convergence acceleration, and high-order spatial/temporal discretizations. Further applications of DUGKS to multiscale flow physics are also desired.