GKS and UGKS for High-Speed Flows

The gas-kinetic scheme (GKS) and the unified gas-kinetic scheme (UGKS) are numerical methods based on the gas-kinetic theory, which have been widely used in the numerical simulations of high-speed and non-equilibrium flows. Both methods employ a multiscale flux function constructed from the integral solutions of kinetic equations to describe the local evolution process of particles’ free transport and collision. The accumulating effect of particles’ collision during transport process within a time step is used in the construction of the schemes, and the intrinsic simulating flow physics in the schemes depends on the ratio of the particle collision time and the time step, i.e., the so-called cell’s Knudsen number. With the initial distribution function reconstructed from the Chapman–Enskog expansion, the GKS can recover the Navier–Stokes solutions in the continuum regime at a small Knudsen number, and gain multi-dimensional properties by taking into account both normal and tangential flow variations in the flux function. By employing a discrete velocity distribution function, the UGKS can capture highly non-equilibrium physics, and is capable of simulating continuum and rarefied flow in all Knudsen number regimes. For high-speed non-equilibrium flow simulation, the real gas effects should be considered, and the computational efficiency and robustness of the schemes are the great challenges. Therefore, many efforts have been made to improve the validity and reliability of the GKS and UGKS in both the physical modeling and numerical techniques. In this paper, we give a review of the development of the GKS and UGKS in the past decades, such as physical modeling of a diatomic gas with molecular rotation and vibration at high temperature, plasma physics, computational techniques including implicit and multigrid acceleration, memory reduction methods, and wave–particle adaptation.


Introduction
High-speed flows are usually involved in aeronautical and aerospace engineering, such as in the re-entry of spacecraft, launch of rockets, and near-space vehicle cruising. As computational fluid dynamics (CFD) calculations start to play important roles in many fields, numerical algorithms for high-speed flow simulations become crucial as well, especially for the cases where it is difficult for experiments to be conducted. However, there are still great challenges for the CFD method in the aspects of reliability, accuracy, and robustness in the simulation of high-speed flows. When the speed of flight gets very high, the flow physics becomes very complicated with the emergence of highly non-equilibrium phenomena. Multiple-scale transport processes in space and time may occur in high-speed flight, where both continuum and rarefied flows would appear around a vehicle. Moreover, due to high temperature, real gas effects such as molecular rotation, vibration, chemical reactions, ionization, and plasma physics have to be taken into account. The gas-kinetic scheme (GKS) [1][2][3] and the unified gas-kinetic scheme (UGKS) [4][5][6] are constructed from the gas-kinetic theory and attempt to recover the flow dynamics with a solid physical foundation. With the incorporation of direct physical modeling in the algorithms, realistic flow physics can be recovered in the schemes, such as the excellent performance in the capturing of multiscale transport in the continuum and rarefied flows.
The GKS proposed by Xu [1,2] aims to obtain the Navier-Stokes (NS) solutions in the continuum regime. Based on the finite volume framework for the macroscopic conservation laws, the GKS adopts the time evolution solution of the kinetic model equation to construct the flux function across the cell interface. For the continuum flow the initial gas distribution function at the beginning of each time step is constructed from the Chapman-Enskog (CE) expansion [7]; the time accurate solution at a cell interface combines the inviscid and viscous terms in the flux evaluation and recovers a physical evolution process from the initial non-equilibrium to the equilibrium one. Moreover, the upwind scheme and a generalized Lax-Wendroff-type central difference method are dynamically coupled in the flux function, and the corresponding limiting schemes could be achieved respectively in the discontinuous and continuous flow regions. This makes the GKS robust and accurate for capturing both shock discontinuities and smooth flow structures. Another distinguishable feature of the GKS is the multi-dimensional property of the flux function [3,8], which incorporates both normal and tangential variations of the flow field in the flux computation, and makes the GKS suitable for flow simulations on unstructured mesh. In order to simulate high-speed flow, the GKS has been further developed in the aspects of physical modeling and numerical algorithms, such as the kinetic boundary condition [9,10], multiple temperature model [11][12][13][14][15][16][17][18][19][20], multi-component and reactive flow [21][22][23][24][25], and implicit scheme and multigrid acceleration [3,[26][27][28][29][30][31][32][33]. By employing a modified relaxation time [10,13,34], the corresponding constitutive relationship in the GKS have been improved so that the GKS can be applied in the near continuum regime. Due to the assumption of a small relaxation time and the use of CE expansion for the distribution function, a continuous particle velocity space is used in GKS and the moments of a gas distribution function can be integrated analytically. However, the applicable regime of the GKS is limited to describing the equilibrium and near-equilibrium flow.
With the adoption of a discrete particle velocity space, the UGKS [4][5][6] is a much enhanced GKS for continuum and rarefied flow simulation. Highly non-equilibrium flow physics can be described and captured in the UGKS by employing a discrete distribution function [35][36][37][38]. In other words, the GKS can be regarded as a limiting case of UGKS in the NS regimes, where under intensive particle collisions the distribution function in UGKS will automatically reach the near-equilibrium CE expansion. The UGKS takes direct modeling on the numerical discretization scales, i.e., the mesh size and time step, aiming to construct the corresponding flow physics on the observational scale for a better description of a multiscale transport process with high efficiency. The key ingredients of the UGKS are to follow the basic conservation laws of the macroscopic flow variables and the microscopic gas distribution function in a discretized space, and to construct a multiscale flux function from the integral solution of the kinetic equation, taking into account the accumulating effects of particles' collision during the transport process in the scale of a numerical time-step.
For high-speed flow simulation, the UGKS can include real gas effects, such as molecular rotation and vibration [39][40][41], multi-component gas mixture, and plasma physics [42,43]. In contrast to the GKS, for which there is a continuous velocity space in which the distribution function can be expressed by the macroscopic flow variables and their gradients, sufficient discrete velocity points should be provided in order to capture the local non-equilibrium distribution in the UGKS. Thus, the computational cost and memory consumption are huge in the UGKS, especially for high-speed rarefied flow simulation with discrete points to cover a six-dimensional physical and velocity space. Therefore, many numerical techniques have been developed and implemented in the UGKS to increase the computational efficiency and reduce memory cost, such as unstructured mesh computation [44,45], moving grids [46,47], velocity space adaptation [47,48], memory reduction [49,50], wave-particle adaptation [51,52], implicit algorithms [53][54][55][56], and fur-ther simplification and modification [57,58]. With these treatments, the UGKS becomes a powerful tool to solve multiscale problems, and shows great advantages in the simulations of high-speed and non-equilibrium flow with a large variation of local Knudsen number in a single computation.
In this paper, we will give a review of the GKS and the UGKS for high-speed flow simulations. The development of the GKS and the UGKS in the aspects of physical modeling and numerical algorithms will be introduced in Sections 2 and 3, respectively. A summary of the GKS and the UGKS for high-speed flow simulation will be presented in the Section 4.

Gas-Kinetic Scheme and Kinetic Boundary Condition
In rarefied gas dynamics, the fundamental governing equation is the Boltzmann equation [59][60][61][62], which employs the seven-dimensional gas distribution function f = f (u, v, w, x, y, z, t) to describe the time evolution of the probability density of gas particles around a microscopic velocity u = (u, v, w) T , at physical location x = (x, y, z) T and time t. Due to the complexity of the full Boltzmann collision term, the kinetic models, such as the Bhatnagar-Gross-Krook (BGK) [63], Shakhov [64], and ellipsoidal statistical BGK (ES-BGK) [65] models, are usually employed. The BGK model equation can be written as where the collision term is approximated by a relaxation process of the distribution function approaching a local equilibrium state. The macroscopic flow variables can be computed from the microscopic distribution function, such as the conservative flow variables, i.e., the densities of mass, momentum, and energy where dΞ = dudvdw and ψ = (1, u, v, w, 1 2 (u 2 + v 2 + w 2 )) T . In the BGK model (1), τ is the relaxation time or particles' mean collision time, which denotes the average time interval between two successive collisions. The equilibrium state g is a Maxwellian distribution where λ is related to the temperature T by where m 0 denotes the molecular mass of gas particle and k B is the Boltzmann constant. Based on the gas-kinetic theory, the NS equations can be derived from the Boltzmann equation or the simplified kinetic model equations using CE expansion [7] when τ is small in the continuum regime. The GKS proposed in [1,2] is a hydrodynamic NS solver based on the kinetic theory in the framework of the finite volume method.
For a control volume i during the time step ∆t = t n+1 − t n , the conservation law of macroscopic flow variables gives where V i is the volume of cell i, N(i) denotes the set of the face-bordered neighbors of cell i, and j is one of the neighboring cells. The interface between cells i and j is denoted by ij, so S ij and F ij are the interface area and time step-averaged fluxes across cell interface ij, respectively. The interface flux in GKS is constructed from the integral solution of the BGK model equation along the characteristic line is the initial distribution function around cell interface x ij at the beginning of each time step t 0 , and it adopts a linear distribution in space and the first-order CE expansion Here, l, r denotes different values on the left and right sides of the cell interface. The initial distribution function is fully determined by the macroscopic flow variables and their gradients. The equilibrium state in space and time is expanded as which can be computed from the spatial and temporal gradients of the macroscopic flow variables at the cell interface. The integral solution describes the evolution process of an initial non-equilibrium distribution approaching a local equilibrium one with the time increment. Since the transport from the initial distribution function f 0 has an upwind mechanism, such as a kinetic flux vector splitting scheme, and the integral evolution of the equilibrium state g has a Lax-Wendroff-type central difference formulation, the GKS automatically provides a transition between the upwind and central difference according to the ratio of numerical time step ∆t over the relaxation time τ. For example, in the discontinuous shock region, the numerical relaxation time τ is large and the GKS becomes an upwind scheme to capture the discontinuity; in the smooth region at a small τ relative to the time step ∆t, the GKS recovers the central difference scheme with high-order accuracy. Another distinguishable feature of the GKS is that both the normal and tangential gradients along the cell interface participate in the flux evaluation, which results in a multi-dimensional hydrodynamic NS solver [3]. The multi-dimensional property makes the GKS suitable for flow simulation on unstructured mesh. In the GKS, the kinetic boundary condition can easily be applied for a better simulation of high-speed flow in the near-continuum regime, where the slip boundary condition can be automatically recovered [9,10]. Without loss of generality, assume that an isothermal wall is on the right side of the inner flow field. A time-accurate gas distribution function for incident particles can be evaluated by one-side interpolation and the incident particles will be reflected from the wall with a Maxwellian distribution where V w and W w are the tangential moving velocities of the wall and λ w is computed from the given wall temperature. The only unknown variable, which is the reflected density ρ w , can be calculated from the non-penetration condition of a solid wall u>0 u n f in dΞdt + u<0 u n g out dΞdt = 0 In the study of paper [66], the hypersonic flow around a blunt cone is computed at the freestream condition of Mach number Ma ∞ = 10.6, Reynolds number defined with respect to the head radius Re = 1.1 × 10 5 , and angle of attack AoA = 20 • . The study provides a detailed comparison of the heat flux on the wall between the GKS solution and experimental data, which show excellent agreement. Many numerical cases have been computed for validating the GKS in high-speed flow simulations, such as supersonic flow past a cylinder [3,[67][68][69], type-IV shock-shock interaction [3], flow past a double cone [3,29], shock boundary layer interaction, double Mach reflection [68,69], hypersonic ramp [68,70], flow around a hollow cylinder flare [9,10], transonic flow around a Falcon Business Jet configuration [71], and a viscous shock tube [72][73][74]. The GKS shows excellent performance in the aspects of robustness and accuracy for simulating high-speed viscous flows.
The kinetic model equation with multiple translational temperatures can be written as where g M is a Maxwellian distribution with one temperature as given in Equation (3) and g is an intermediate state with multiple translational temperature where λ x , λ y , and λ z are related to the x-, y-, and z-directional translational temperatures T x , T y , and T z , respectively. The average temperature in the equilibrium state g M is evaluated from 3 Z t is a parameter denoting the ratio of the relaxation time from intermediate state g to g M over that from non-equilibrium distribution f to the intermediate state g. When Z t equals one, Equation (12) becomes the standard BGK model equation with one translational temperature; for Z t → ∞, the distribution function f approaches a Gaussian distribution with different translational temperatures in the x-, y-, and z-directions [12]. In the original study of [11], a weight function 1/2 is used with the consideration that the average temperature should return to the same temperature as that in the single temperature BGK model, which is equivalent to the case with Z t = 2 in Equation (12). In [18,75], Z t = 1 is employed to re-write the collision process of the BGK model into two physical subprocesses. Several numerical cases, including shock wave, Couette flow, Poiseuille flow, and lid-driven cavity flow, have been carried out to show the capability of the multiple temperature kinetic model in capturing non-equilibrium physics in the near-continuum regime [11,12,75]. In general, the translational temperature should be constructed as a second-order symmetric tensor, which deserves further investigation. Besides the translational temperature, the kinetic model in Equation (12) can be extended to rotational and vibrational degrees of freedom as well to simulate diatomic gas flows. Adding two internal degrees of freedom into Equation (12), the kinetic model equation for diatomic gas can be written as where the equilibrium state becomes Here, ξ 2 = ξ 2 1 + ξ 2 2 denotes two rotational degrees of freedom. The intermediate state has different translational and rotational temperatures Z r is the so-called rotational collision number.
In the studies of [13][14][15], slightly different from the above equation with two relaxation processes, the BGK kinetic model with different translational and rotational temperatures is used, and the rotational collision number Z r is artificially added into the source term of the rotational energy moment equation, which takes into account the longer relaxation time of the rotational energy to reach equilibrium. In [16,76], the two-stage relaxation model is more consistent with the well-known Rykov model [77] for diatomic gas simulation.
In addition, the vibrational energy has been included as well in the GKS to compute high Mach number shocks and type IV shock-shock interaction [17,20]. It should be pointed out that in order to capture non-equilibrium in a near-continuum regime, a generalized collision time τ * has been proposed [34], and it was used in the previously mentioned studies [10,11,13,14,16]. The generalized collision time is computed by where D = ∂ t + u · ∂ x and . . . = (. . . )φdΞ. Here, different moments φ can be used to provide different values of τ * for the viscosity term and heat conduction coefficient [14]. τ * depends not only on the macroscopic flow variables but also their gradients. It can be regarded as a generalization of constitutive relationship with the inclusion of higher-order terms resulting in better capability of capturing rarefied non-equilibrium effects.

Multi-Component and Reactive Flow
For high-speed and high-temperature flows, a multi-component gas mixture is always involved. Numerical modeling and algorithms for multi-component and reactive gas flows have also been developed in the GKS [21][22][23][24][25].
For two-species gas flow, the kinetic equation can be written as where f 1 and f 2 are gas distribution functions for the components 1 and 2, and g 1 and g 2 are the corresponding equilibrium states. In the studies of [21,22], the equilibrium state for each component is a Maxwellian distribution, i.e., where ξ 2 1 = ξ 2 1,1 + · · · + ξ 2 1,K 1 , and ξ 2 2 = ξ 2 2,1 + · · · + ξ 2 2,K 2 . K 1 and K 2 denote the internal degrees of freedom for two components. It is assumed that g 1 and g 2 share the same temperature and velocity. For a two-component gas mixture, the compatibility condition gives where 2 )) T . The common temperature and velocity can be obtained from the conservation law In the study of Xu [21], it is assumed that the relaxation time τ is much smaller than the time step, and the exchange of momentum and energy can be finished instantaneously, which results in both species following the common temperature and velocity at any time.
Similar to the GKS with single-component gas, the integral solution of the kinetic model is used as well in the multi-component GKS to construct the interface flux for each component where the initial gas distribution function f 0 is constructed by interpolation from each side of the cell interface. g 1 and g 2 are the equilibrium states for both species with the same temperature and velocity computed from the compatibility condition. After obtaining the fluxes across the cell interface for both species, the flow variables can be updated for each species.
In the studies of Lian [22] and Pan et al. [25], the chemical reaction is considered. The ZND model proposed by Zel'dovich, von Neumann, and Doering has been incorporated in the multi-component GKS in a splitting way, where an ordinary differential equation (ODE) is solved to take into account the source term from the chemical reaction, i.e., where K(T) depends on the temperature T and is computed by Here, K 0 is a positive constant and E + denotes the activation energy. α = 0 is used in the studies of [22,25]. Equation (24) can be solved by the one-step forward Euler method or the multi-step Runge-Kutta method. In the framework of the two-stage fourth-order method, high-order GKS with multi-component gases has been developed by Pan et al. [25].
In the study of [25], the case of Rayleigh-Taylor instability has been computed to validate the multicomponent GKS with gravitational force. The development of the instability from the initial condition is well captured by the multi-component GKS. In addition, the numerical test cases, including the Sod shock tube, gas-vacuum expansion, ZND detonation, 2D viscous reactive flow, and shock-bubble interaction, have also been computed [21][22][23][24][25] to show the capability of the GKS for simulation of multi-component and reactive gas flows.

Implicit Scheme and Multigrid Acceleration
Implicit scheme and multigrid methods are the commonly used acceleration techniques for solving partial differential equations (PDEs), which have been widely applied in the CFD calculations. The implicit scheme can not only increase the convergence efficiency but also improve the numerical stability, which is also vital for high-speed flow simulations. The implicit and multigrid acceleration techniques have been developed in the GKS as well for fast solving of both steady and unsteady state problems [3,[26][27][28][29][30][31][32][33].
For unsteady flow evolution in a large time step ∆t = t n+1 − t n , the discrete governing equations can be written as where = 0.5 gives a second-order accurate Crank-Nicolson (CN) scheme and = 1 corresponds to the first-order accurate backward Euler method. We can re-write the governing equations in the ∆-form with an intermediate approximate solution for inner iteration, ij . The residual on the right hand side of the equation gives It can be found that when the residual R (s) approximates zero after several iterations, w (s) converges to w n+1 and becomes the solution to Equation (26). Since the residual determines the convergent solution, the fluxes on the right hand side, i.e., F n ij and F (s) , are computed by the GKS as the same as that in the explicit scheme; for terms on the left hand side, an arbitrary approximation can be carried out as long as it can result in a convergent solution. Therefore, ∆t can adopt any positive value, and the flux variation on left hand side can be approximated by the first-order Euler flux [78,79], where the Euler-equation-based flux T is U n = U · n ij is the projected macroscopic velocity along the normal direction of the cell interface ij. Γ ij denotes the spectral radius of the flux Jacobian matrix. For viscous flow, an additional term Γ ν is required for stability considerations and it is computed by Then, the implicit system formed by Equation (27) can be solved by the numerical algorithms commonly used in CFD computation, such as the lower-upper symmetric Gauss-Seidel (LU-SGS) method [80][81][82] and generalized minimal residual (GMRES) method [83]. From the study of Tan and Li [30], it is suggested to use the Jacobian matrix of kinetic flux vector splitting (KFVS) as the pre-conditioner, which can provide better performance for steady high-speed flow simulations in the aspects of robustness and convergence. In the study of Li et al. [31], the temporal discretization involving solutions on three time levels has been utilized in the development of implicit GKS for unsteady flow simulation, by adopting the commonly used dual time-stepping strategy in conventional CFD methods.
For a steady-state solution, the implicit scheme can be further simplified, and the algorithm would be equivalent to the case that = 1, ∆t → ∞ in Equation (28). One evolution step from t n = 0 to t n+1 → ∞ would provide the convergent solution satisfying The multigrid method accelerates the solution for a fine grid by computing corrections on a coarse grid to eliminate low-frequency errors efficiently [84][85][86]. It achieves a fast development [87][88][89][90][91] and has been widely used in CFD calculations to solve the Euler and NS equations [80,92,93]. In the studies of [26,28,94], the multigrid method has been used in the computation of GKS. With the implicit scheme and multigrid method, the computational efficiency of the GKS can usually be improved by one or two orders of magnitude.

Unified Gas-Kinetic Scheme for High-Speed Flows
With the discrete velocity distribution function, the unified gas-kinetic scheme (UGKS) proposed by Xu et al. [4,5] can be regarded as a highly enhanced version of the GKS, which is capable of simulating continuum and rarefied flows in all Knudsen regimes. By adopting the integral solution of the kinetic equation along a characteristic line, the UGKS considers the coupling effect of particles' free streaming and collision in the transport process, and can well describe multiscale flow transport with no need to restrict the numerical resolution on the scales of mean free path and particle collision time. Similarly, the discrete unified gas-kinetic scheme (DUGKS) proposed by Guo et al. [95,96] is capable of multiscale flow simulations in all Knudsen regimes by adopting the discrete solution of the kinetic equation along a characteristic line. After a decade of development and improvement, the UGKS and DUGKS have gained great success in solving multiscale transport problems, such as radiative transfer [97][98][99][100], phonon transport [101,102], plasma physics [42,103,104], neutron transport [105,106], multicomponent and multiphase flow [107][108][109][110][111], granular flow [112][113][114], and turbulent flow [115][116][117].
For high-speed flow simulation, besides the physical modeling of real gas effects such as molecular rotation and vibration, chemical reaction, ionization, and plasma, another big challenge for the UGKS is the huge computational cost resulting from the massive discrete velocity points. In the following, we will give a brief review of the basic algorithm of the UGKS, and introduce the development of the UGKS in physical modeling and numerical treatment for high-speed flow simulations.

Basic Algorithm
For a control volume i during the time step ∆t = t n+1 − t n , the governing equations of the gas distribution function and conservative flow variables are (32) and which describe the distribution function change due to particle transport and collision, and the flow variable variation due to flux transport across the surfaces of the control volume.f ij is the time-averaged distribution function at cell interface ij over time step ∆t and u n,ij = u · n ij denotes the normal component of microscopic velocity u across the interface ij. J ( f , f ) represents particle collision and in the UGKS the trapezoidal rule is used In Equation (33), we re-write the macroscopic conservation law the same as that in the GKS. The main difference here is that the time-averaged macroscopic flux F ij is supported by the underlying microscopic physics, and it is computed by taking moments of microscopic fluxes F ij with the compatibility condition of the collision term With the discrete velocity distribution function, the compatibility condition cannot be naturally satisfied due to the numerical quadrature error, so several conservative kinetic methods [54,[118][119][120][121][122][123] have been constructed, which could maintain the compatibility condition in the discrete form and show that the conservation of collision term can slightly release the strict requirement of velocity space discretization, and gain better convergence for implicit calculations.
Equations (32) and (33) directly describe the conservation of gas distribution function and macroscopic flow variables, which are supposed to be valid at arbitrary discretization scales regardless of variations in mesh size and time step. Following the basic conservation laws is the fundamental step in the construction of a multiscale UGKS on the discretization scales. With this consideration, a new conserved DUGKS has been proposed [124] as well, where the conservative flow variables are updated from the macroscopic conservation law instead of by taking moments of the auxiliary distribution function as in the original DUGKS. Regardless of the quadrature rule in velocity space, the conservation law is always satisfied on the macroscopic level in the conserved DUGKS, which could benefit in the numerical calculations where conservation is of much importance, such as for highly nonequilibrium flow and plasma physics [48,104,125]. Following these conservation laws, we can see from Equation (33) that the evolution mainly depends on the microscopic flux F ij across cell interfaces.
In the UGKS, a multiscale flux is constructed from the integral solution of the BGK model equation, i.e., Equation (6). Different from that in the GKS, f 0 ( x) in the UGKS is the initial distribution function obtained from the previous step evolution, which could be a highly non-equilibrium distribution instead of the near-equilibrium state from the CE expansion. With the distribution function described at discrete velocity points, the UGKS can make full use of the advantage of the multiscale properties of the integral solution. Besides providing a transition between the upwind scheme and central difference discretization, the integral solution describes the evolution of the distribution function from the initial state to the local equilibrium with the accumulation of particle collision. When the discrete time-scale is much larger than particle mean collision time ∆t τ, the local equilibrium state dominates in the flux function and the flow behavior performs as macroscopic wave propagation and interaction; when the time scale is small ∆t < τ, the initial distribution function is dominant and the flow physics is the free transport of microscopic particles. Moreover, the integral solution provides a transition from kinetic scale to hydrodynamic scale, where the observed local physics is determined by the cell's Knudsen number, i.e., τ/∆t. For a well-resolved solution, the flux function is able to describe the flow physics on the corresponding modeling scale of the numerical time step. The scale adaptation through the ratio of physical scale and numerical resolution (τ and ∆t) is the key to construct a truly multiscale numerical method. The traditional PDE approach usually requires the mesh size and the time step to be small enough to match the physical modeling scale in the construction of the governing equations, while the UGKS has no specific requirement for the numerical discretization, except for the stability condition and necessary spatial resolution for well-resolving local flow structure. Owing to the adaptive multiscale flux function and direct modeling concept, the UGKS has better efficiency in the numerical simulation of multiscale problems with large variations in Knudsen numbers.

Diatomic Gas with Molecular Rotation and Vibration
For diatomic molecules, e.g., N 2 and O 2 , the characteristic temperature of rotation is around 3 K [59], so the rotational degrees of freedom have already been excited at room temperature. However, for the vibrational degrees of freedom, the characteristic temperatures of N 2 and O 2 are 3371 K and 2256 K [59], respectively. Therefore, for hightemperature flows at high Mach numbers, both the rotational and vibrational degrees of freedom should be taken into account.
Liu et al. [39] employed the Rykov model to construct the time evolution solution of the gas distribution function in the UGKS for simulation of diatomic gases. The rotational degree of freedom is included and modeled by controlling the energy exchange between translational and rotational energy through the relaxation model. The Rykov kinetic model equation [77] gives The equilibrium states are expressed as whereg and θ is defined by The macroscopic flow variables, such as pressure, temperature, and heat flux for translational and rotational degrees are defined by The parameters used in the above model are ω 0 = 0.2354, ω 1 = 0.3049, and 1/σ = 1.55 for nitrogen gas, and ω 0 = 0.5, ω 1 = 0.286, and 1/σ = 1.55 for oxygen gas.
In order to construct the diatomic UGKS, the Rykov model can be re-written in the same form as the BGK model where the equivalent equilibrium state g eq is Then the UGKS with the Rykov model can be constructed analogous to that for the BGK model. The updating of the macroscopic flow variables can be carried out by In the calculation, the energy relaxation term in the Rykov equation is modeled using a Landau-Teller-Jeans relaxation. The particle collision time multiplied by rotational collision number defines the relaxation rate for the rotational energy equilibrating with the translational energy. The rotational collision number Z r is computed by Z r = Z * r 1 + (π 3/2 /2) T /T tr + (π + π 2 /4)(T/T tr ) (45) whereT is the characteristic temperature of intermolecular potential. For N 2 within a temperature range from 30 K to 3000 K, the values of Z * r = 23.0 andT = 91.5 K are used. In the study of [39], the normal shock has been computed for nitrogen gas at different Mach numbers. The collision number is set at a constant value Z r = 2.4. The comparisons of the density and temperature between the UGKS and DSMC at Ma = 1.53, 4.0, 5.0, and 7.0 are provided. With the differences of kinetic models, the UGKS still obtains acceptable solutions of translational and rotational temperatures in comparison with the data obtained by the DSMC method.
Taking the molecular vibration into account, a kinetic model [40,41] has been proposed to simulate the diatomic gases with activated vibrational degrees of freedom. The kinetic model equation can be written as where Z r and Z v denote the rotational and vibrational collision numbers. The equilibrium states are g trv = ρ λ π 3 2 e −λc 2 λ π e −λ(ξ 2 Kv v g tr,v = ρ λ t,r π 3 2 e −λ t,r c 2 λ t,r π e −λ t,r (ξ 2 Kv v Here, λ, λ t , λ r , λ v , and λ t,r are related to the fully relaxed, translational, rotational, vibrational, and partially relaxed temperatures, respectively. Specifically, we have where K v is the number of vibrational degrees of freedom. This collision model consists of three terms, including the relaxation processes of elastic collision, inelastic collision between molecular translation and rotation, and the energy exchange from vibrational degrees of freedom to the translational and rotational ones. Equation (46) is more like a BGK equation with a fixed Prandtl number, which can be extended to a Shakhov-like model by using orthogonal polynomials in order to obtain the correct relaxation rate of high-order moments [40,41]. With the provided kinetic model, the UGKS for diatomic gases with rotational and vibrational degrees of freedom can be constructed. Different from the translational and rotational degrees of freedom, the vibration degree of freedom K v is a temperature-dependent variable. According to the harmonic oscillator model, the specific vibrational energy associated with a characteristic vibrational temperature Θ v is then according to the equal partition to each degree of freedom, the vibrational degree of freedom K v can be determined by In the study of [40], flow around the Apollo re-entry capsule has been computed at Ma = 10.2, Kn = 0.067 with a freestream temperature 142.2 K and a fixed wall temperature 300 K. Nitrogen gas is considered and the collision numbers adopt Z r = 5 and Z v = 10. The comparison of flow fields between the diatomic UGKS with the Rykov model and the current model with a vibrational degree of freedom is presented. It is shown that due to the vibrational degree of freedom, the translational and rotational temperature distributions obtained by the current model are lower than those computed by the Rykov model. A more realistic solution has been obtained.

Multi-Component Gas Mixture and Plasma
The direct modeling methodology of the UGKS has also been applied in the plasma simulation [42,43]. The dynamics of a fully ionized plasma are modeled by the Fokker-Planck-Landau (FPL) equations on the kinetic level where f α ( x, u, t) is the velocity distribution function of species α (α = i for ion and α = e for electron). F α = e( E + u α × B) is the averaged electromagnetic force. The collision term Q α,β ( f α , f β ) describes the binary collisions between charged particles with long-range Coulomb interactions In order to overcome the complexity and high computational cost of the nonlinear Landau collision term, a single BGK-type collision operator proposed by Andries, Aoki, and Perthame (AAP model) is employed to model the collision process [126].
In the AAP model, one global collision operator is used for each component to take into account both self-collision and cross-collision, and the kinetic equations read where post-collision distribution f + is a Maxwellian distribution The parameters T α and U α are connected to the macroscopic properties of individual components by where ν αr denotes the interaction coefficients that measure the strength of intermolecular collision. The relaxation time is determined by τ α = 1/ ∑ r ν αr . The parameter ν α r is determined by molecular models and the hard sphere model is used [127]. The average electric field E and magnetic field B follow the Maxwell equations, where c is the speed of light and 0 is the vacuum permittivity, which is related to the vacuum permeability ν 0 by c = (ν 0 0 ) −1/2 . The electromagnetic field satisfies the divergence constraints where e is the charge of a proton. Theoretically, the divergence constraints will always hold if they are initially satisfied. However numerical techniques are needed to make sure that the divergence constraints are satisfied by numerical solutions. The perfectly hyperbolic Maxwell equations (PHM) [128] are used to evolve the electromagnetic field, which preserve the divergence constraints, where j is the total electric current density j = e(n i U i − n e U e ), and φ and ψ are artificial correction potentials. Munz proved that the propagation speed of magnetic field divergence error and electric field divergence error are γc and χc [128]. This scheme is built on the BGK-Maxwell system Equations (54), (57), and (58), which are able to cover the flow regimes of plasma from the collisionless Vlasov regime to the continuum MHD regime.
In the study of Liu et al. [42], the asymptotic limits of the BGK-Maxwell system have been analyzed. In the continuum regime, when the interspecies collision is strong, the gas mixture behaves like dielectric material, and the BGK-Maxwell equations become the Euler equations. For a conductive plasma, the BGK-Maxwell equations can span the complete range from the neutral two-fluid system to resistive-MHD, Hall-MHD, and ideal MHD equations. An implicit-explicit method is employed to renew the velocity distribution function so that the constraint of a small time step due to the large electromagnetic acceleration can be removed. Moreover, the flow physics covered by the UGKS are more general than those from either the collisionless Vlasov equation or MHD equations in the corresponding kinetic or hydrodynamic regimes alone. The UGKS can provide a reliable physical solution in the transitional regime as well, which has not been fully explored before from the particle-based and MHD-based numerical methods.
Based on the conserved DUGKS, a finite-volume direct kinetic method for electrostatic plasma has been presented [103], where the BGK-Vlasov-Poisson system is solved with un-splitting treatment of particle transport and collision. In addition to the non-quasineutral plasma simulation, based on the reformulated BGK-Vlasov-Poisson system, a novel DUGKS has been presented [104] for multiscale plasma simulation with a wide range of Knudsen numbers and normalized Debye length covering hydrodynamic, kinetic, quasi-neutral, and non-quasi-neutral regimes. Numerical test cases, including plasma sheath, linear and nonlinear Landau damping, two-stream instability, Brio-Wu shock tube, Orszag-Tang vortex, magnetic reconnection, and bump-on-tail instability, have been calculated in different flow regimes, which show the multiscale properties of the UGKS and DUGKS in plasma physics simulation [42,43,103,104,129].

Computational Techniques
In this section, the computational techniques for improving the efficiency of the UGKS in numerical simulations of high-speed gas flow will be introduced, including the implicit scheme, multigrid method, parallel computation, adaptive physical mesh and velocity space, memory reduction techniques, and the novel wave-particle adaptation.

Implicit UGKS and Multigrid Acceleration
The implicit UGKS (IUGKS) has been developed [53,55,130] by transforming the UGKS on a discretization scale into a semi-discrete form as is usually carried out for the traditional CFD methods. The semi-discrete governing equations of the macroscopic flow variables and the distribution function are and The governing equations in the semi-discrete form usually describe the instant variation of the flow field, which implies a time scale of t → 0. Specifically, F ij and F ij denote the instant macroscopic and microscopic fluxes across the cell interface ij. However, the UGKS is constructed based on the integral solution of the kinetic model, describing the multiscale transport process in a finite time-step, and the fluxes F ij and F ij are not only related to the local physical state τ, but also depend on the mesh-size-determined time step ∆t. Therefore, the fluxes F ij and F ij in Equations (60) and (61) should employ the time-averaged fluxes over a finite time-step ∆t s instead of the instant ones for constructing the implicit UGKS, where the time step ∆t s is used to average the time-dependent fluxes in the explicit scheme. Considering the flow physics in a local region, ∆t s is determined by the resolution of computational mesh and the maximum particle speed, e.g., for one-dimensional case ∆t s = CFL ∆x max |u| (62) In order to distinguish the time step ∆t s to average fluxes from the following numerical time-marching step ∆t, it should be pointed out that ∆t s is only used to determine the coefficients of the initial distribution function and equilibrium state in the flux function.
For unsteady flow evolution in a time step ∆t = t n+1 − t n , in addition to the macroscopic governing equation of conservative flow variables (26), the discrete governing equations can be written as Since the equilibrium state g n+1 i has one-to-one correspondence to w n+1 i and depends on the distribution function f n+1 i , Equations (26) and (63) result in a coupled implicit nonlinear system. It is very difficult to directly solve the large implicit system coupling macroscopic and microscopic governing equations, and the treatment of implicit equilibrium state g n+1 in the collision term is important for the convergent efficiency in the continuum regime. Zhu et al. [53,131] solved the Equations (26) and (63) in an alternative way and used the most recent solved conservative variables to discretize the equilibrium state in the microscopic implicit equation. Specifically, for an intermediate solution during the alternative solving process, the governing equation can be written as where the quantities in the ∆ form are given by ∆Q (s) = Q n+1 − Q (s) , denoting the correction of a specific variable Q. The residual on the right hand side of Equation (64) is It should be noted that the fluxes in the residuals (28) and (65) are computed the same as that in the explicit UGKS with a time step ∆t s .g in Equation (65) is the equilibrium state computed from the most recently updated conservative variables ∆ w (s+1) i in the inner loops of solving Equation (27).
In the study of Zhu et al. [130], the one-dimensional case of advection of density perturbation has been employed to test the temporal accuracy of the IUGKS. The time accuracy test is conducted for the cases at Kn = 0.001, 0.01, 0.1, 1, and 10 on a uniform mesh with 400 cells. The IUGKS achieves second-order accuracy for all cases from continuum to free molecular flows. The unsteady Rayleigh flow induced by a moving plate with higher temperature has been computed to validate the efficiency of the implicit scheme. It is shown that the IUGKS can be ten times faster than the explicit scheme [130].
For steady-state solutions, temporal accuracy can be ignored, and the time discretization of the implicit scheme could use = 1 for a backward Euler scheme. In [53,130], in order to overcome the stiffness problem in the continuum flow, the macroscopic equations are solved with the implicit Euler flux to drive the convergence of the IUGKS; see Equations (28) and (29). This makes the IUGKS efficient in the low Knudsen number cases. In order to further improve the performance of the IUGKS for viscous flows, Yuan et al. [132] took into account the NS terms during the macroscopic iterations. A similar idea of driving the convergence in the continuum regimes can be found in the development of the general synthetic iterative scheme (GSIS) [133,134]. From the study of Yuan et al. [132], it can be found that the efficiency can be further improved by one order of magnitude on the basis of the original IUGKS for NS solutions in the continuum limit.
In addition, Zhu et al. [135] developed the IUGKS with multigrid acceleration, which further improves the convergence of the IUGKS for steady-flow simulations. It is commonly known that an iterative algorithm can reduce the high-frequency errors faster than the low-frequency ones. The multigrid method can reduce the low-frequency errors on a fine mesh more efficiently by transitioning them to a coarser mesh, where the errors become high-frequency ones with respect to the coarse mesh size, and can then be eliminated faster. In the IUGKS, both the implicit systems of macroscopic and microscopic equations are solved by the multigrid method, which consists of pre-smoothing, coarse grid correction, and post-smoothing processes, and involves numerical interpolation between fine and coarse grids [135]. Since the governing equation of macroscopic flow variables for implicit iteration is nonlinear, and that of the microscopic distribution function is linear, the full approximation storage (FAS) scheme and the correction scheme (CS) [88,89] are used to solve the macroscopic and microscopic implicit systems, respectively.
In the study of [135], simulations of lid-driven cavity flows have been carried out at different Knudsen numbers to validate the efficiency of the multigrid IUGKS. Cases at three different Knudsen numbers, i.e., Kn = 10, 1, 0.075, have been tested, from which obvious accelerating effects of the multigrid method on the IUGKS can be observed. In the high Knudsen number cases at Kn = 10 and 1, the multigrid method is about 3 times faster than the original implicit scheme and in the case at Kn = 0.075 the acceleration rate can be increased by up to 8 times.

Parallel Strategy
In the CFD algorithms for the Navier-Stokes equations, the domain decomposition method is commonly used in the parallel computations. Since the kinetic solvers take numerical discretization on both the physical space and the velocity space, the parallel strategy is more flexible for large-scale simulations.
For the UGKS on Cartesian grids, Ragta et al. [136] adopt the parallel strategy based on the domain decomposition, and investigate the parallelization performance up to thousands of cores. Chen et al. [49] developed an implicit kinetic method with memory reduction techniques, which significantly reduces the memory consumption, and the parallel computation based on the discrete velocity distribution function is employed. Li et al. [137] and Tan et al. [105] implemented a hybrid MPI strategy for parallel algorithms in both physical space and velocity space, where a two-dimensional Cartesian topology is used to arrange the physical and velocity blocks. Similarly, Jiang et al. [138] took parallel computing with MPI for the decomposed physical mesh, and used several threads with OpenMP in every MPI process for parallel computation of a discrete distribution function. Parallel algorithms have also been implemented on GPU devices for a two-dimensional UGKS [139].
In the study of Li et al. [137], the parallel speedup ratio has been tested based on a supersonic flow over a sphere. It is shown that the parallel efficiency is around 1 up to 5832 processors, which reveals the good scalability of the parallel UGKS.

Adaptive Mesh
Adaptive mesh refinement (AMR) can provide better spatial discretization with fewer grid points. For the kinetic solvers with discrete phase space, the AMR technique can be applied in the velocity space as well [140][141][142]. In the study of Chen et al. [47], an adaptive UGKS was proposed through the introduction of an adaptive quadtree structure in the velocity space. Together with a moving mesh in the physical space, the UGKS is able to simulate moving fluid-body interaction in the high-speed flows with a better efficiency. Qin et al. [143] presented a simple local discrete velocity adaptation in the UGKS under a uniform background velocity mesh, which avoids the interpolation between different levels of velocity grids. Xiao et al. [144] replaced the discrete velocity points in the UGKS with a continuous velocity space in the region where the NS solutions provided by the GKS are valid. In the studies of [48,125], the unstructured mesh is adopted in the velocity space for the DUGKS, where more discrete velocity points are arranged in the velocity regions that enclose a large number of molecules, and fewer points are needed in the regions with a small amount of molecules, so that the memory consumption can be reduced and the computational efficiency can be increased. Recently, the reduced order modeling (ROM)based approach was introduced in the DUGKS [145] to reduce the discrete velocity space by removing the mesh points that have a negligible effect on the calculation of moments. With these approaches, the UGKS and DUGKS can obtain multiscale solutions with the same accuracy and less computational cost.
With an adaptive mesh in velocity space, Chen et al. [47] computed a moving nozzle flow with the co-existence of both continuum and rarefied regimes in a single computation. This multiscale problem is very challenging for any single-scale-based numerical scheme. However, this is exactly the place where the UGKS can be faithfully applied to obtain accurate solutions in different regimes. For such a multiscale problem, a hybrid approach with domain decomposition will usually be adopted [146][147][148]. However, for the unsteady case, the hybrid approach cannot generally be applied due to the absence of a valid interface.

Memory Reduction
For rarefied flow study, a direct Boltzmann solver based on the discrete velocity distribution function is a common approach to simulate the non-equilibrium flow. However, it requires a huge memory to follow the evolution of a discrete distribution function in a computational space with seven dimensions. This effect would become more severe for numerical computation of high-speed rarefied flows, where a very large region in the velocity space has to be covered. The adaptive methods mentioned in the previous section [47,48,145] could reduce the memory consumption and computational time through optimization of the distribution of discrete velocity points in the velocity space. Furthermore, for steady-state simulation, a memory reduction technique was proposed by Chen et al. [49,149] for solving the Boltzmann kinetic model equations. Yang et al. [50] applied the memory reduction method in the implicit UGKS. Due to the fact that only the spatial connection and variation take effect in a steady-state solution, the non-equilibrium distribution function can be obtained from spatial iterations at each velocity point, therefore only the memory for the distribution function at one discrete velocity point in the physical domain is required, which reduces the memory cost to the same level as that in the hydrodynamic solvers for the Euler and Navier-Stokes equations.
In the study of [49], Chen et al. computed a thermal transpiration flow in a threedimensional square cube with closed walls. In their computation, the spatial domain is discretized by a structured mesh with 60 × 60 × 120 cells, and the velocity space employs 101 × 101 × 101 discrete velocity points. If the distribution function at all points in the discrete velocity space is stored, the memory requirement would be approximately 13 terabytes (TB) with double precision calculation. With the memory reduction technique, the observed data storage for a parallel code is only 181 megabytes (MB) for each process in a total of 128 MPI processes. The memory technique makes it possible to calculate the three-dimensional problems using a discrete-velocities-based kinetic method on a small cluster. However, with the significant memory reduction, the computational cost will not be reduced in the above scheme.

Wave-Particle Adaptation
For rarefied gas simulations, the DSMC method is more preferable for engineering applications than the discrete velocity method, due to the fact that the stochastic particle method takes much less computational cost than the discrete velocity method, especially for high-speed calculations [150][151][152]. Compared to the discrete velocity method, where the discrete velocities are predetermined, the stochastic particle method can be regarded as an optimal velocity space adaptation method, because the particles collide and change their velocity, and appear in the phase space where there should be a non-zero distribution function.
In order to incorporate the adaptive properties of the stochastic methods, the unified gas-kinetic particle (UGKP) method has been developed [153][154][155][156] according to the direct modeling methodology, where stochastic particles are employed to record the nonequilibrium distribution, and the evolution processes are carried out on both the macroscopic and microscopic levels. The UGKP method recovers the UGKS solutions by following the local evolution solution of the kinetic models, which is the key to construct a multiscale method for all Knudsen number flows. Specifically, the particle evolution described by the gas distribution function follows where According to Equation (67), only a portion of particles can move freely in a time step, and the others will encounter collisions with other particles during the transport process. The free transport process can be fully resolved by tracking all the particles, and the contribution of the particles' free motion to the macroscopic fluxes across the cell interface can also be recorded. Instead of moving freely within a whole time step, the collided particles will encounter collision during a time step evolution. Since a portion of particles will suffer many collisions during the collision process, the colliding particles will not be followed. Fortunately, the collective behavior of the colliding particles in the collision process can be well modeled and followed on the macroscopic level based on the evolution of the equilibrium state, which is only related to the macroscopic flow variables and their gradients. Similar to the UGKS, the conservative variables can be updated first with the fluxes contributing from the free transport and collision processes, and then the colliding particles will be re-sampled from a local equilibrium state according to their corresponding conservative variables.
Considering the fact that the colliding particles are re-sampled from a known equilibrium state, there is no need to describe the equilibrium state with discrete particles, but by using the analytic expression itself with high efficiency and low memory-consumption. Therefore, a novel adaptive wave-particle formulation is proposed to describe the multiscale flow physics, resulting in the unified gas-kinetic wave-particle (UGKWP) method [51,52,[157][158][159]. It can be found that if the colliding particles are expressed by an analytic formulation, their contribution to the free transport and collision can be analytically computed as well. The only required treatment for the collided particles is to sample a portion of these particles to take into account the non-equilibrium generating from the free transport mechanism of the equilibrium distribution. There is a competitive balance between the analytic waves and stochastic particles according to local physics. Particles become waves due to collisions and waves transform into particles in order to capture the partially free transport mechanism.
With the adaptive wave-particle formulation, the cost of the UGKWP method can be reduced to a particle method in the highly rarefied flow regime, and it becomes a hydrodynamic NS solver in the continuum regime automatically. It is a novel technique that combines and optimizes the adaptive velocity space method [47] and GKS/UGKS hybrid algorithm [144]. Different from the hybrid methods [146,160,161], which are based on the domain decomposition and solver hybridization, the UGKWP method describes the physical state by an adaptive wave-particle decomposition in each cell with a unified treatment in the whole computational domain.
In the study of [162], hypersonic flow past a circular cylinder has been simulated to show the capability of the UGKWP method for high-speed rarefied flow simulations. The free stream is initialized with the monatomic gas of argon. The UGKWP method obtains consistent solutions with the UGKS for the case at a relatively low Mach number Ma = 5, but is one order of magnitude lower in computational cost and memory consumption. For the case at a very high Mach number Ma = 30, since the memory requirement of the discrete velocity points for the UGKS is unaffordable, only the comparison between the UGKWP method and the particle method is provided. In this computation, the memory cost of the UGKWP method is only 375 MB. The advantage of the particle method with a natural adaptivity in the phase space is well inherited by the UGKWP method for highspeed rarefied flow computations. In the study of Chen et al. [158], the high speed flow at Ma = 10 over a space vehicle in the transition regime at Kn = 10 −3 has been calculated to show the efficiency and capability of the UGKWP method for simulating three-dimensional hypersonic flow over a complex geometry configuration. The UGKWP method shows great power to solve engineering problems involving high-speed non-equilibrium flows, such as flow around near-space or re-entry vehicles.

Conclusions
In this paper, we have reviewed the development of the gas-kinetic scheme (GKS) and the unified gas-kinetic scheme (UGKS) in the aspects of physical modeling and numerical algorithms. Non-equilibrium flow physics, including diatomic gas with molecular rotation and vibration, multi-component and reactive gas, and plasma, have been modeled in order to take the real gas effects into account. The numerical efficiency and robustness of the GKS and the UGKS have been much improved by adopting implicit scheme and multigrid acceleration, using parallel computation and an adaptive mesh.
Specifically, the GKS is a hydrodynamic NS-solver based on the gas-kinetic theory. With the numerical flux constructed from the integral solution of the kinetic model, it achieves adaptive transition between the upwind and central difference discretization in the flux evaluation according to the local physics. Derived from the gas kinetic theory with a gas distribution function, the GKS has a solid foundation and is more flexible for dealing with non-equilibrium flow. With multi-dimensional properties, it can handle complex gas flow and shows excellent performance in the high speed flow simulation. The multitemperature model, multiple distribution function for multi-component and reactive gas flow, or multi-stage relaxation model for diatomic gases have been developed in the GKS for better capturing real gas flow physics. In addition, the constitutive relationship could be improved directly by modifying the relaxation time τ * , so that the GKS can be enhanced to compute low-density non-equilibrium in the near-continuum regime. However, the GKS adopts the distribution function from the CE expansion, which is analytically expressed by the macroscopic flow variables and their gradients, and a small deviation from the Maxwellian equilibrium is assumed, so the capability of the GKS to capture highly nonequilibrium flows is still limited.
By using a discrete velocity space, the UGKS can capture highly non-equilibrium distributions, and it makes full use of the advantages of the integral solution of the kinetic model in describing a multiscale transport process. In the UGKS, physical modeling has been carried out to consider a diatomic gas with molecular rotation and vibration, a multicomponent gas, and plasma physics, which could provide more accurate solutions for high-speed and high-temperature flows. With the implicit discretization and multigrid method, the computational efficiency of the UGKS has been much improved. Parallel strategies based on domain decomposition have been developed in both the physical and velocity spaces, which enables the UGKS to handle large-scale computations. In order to reduce the memory cost of discrete particle velocities, an adaptive mesh is adopted in the velocity space. With the memory reduction technique, the memory requirement of the UGKS can be much reduced to the level of a macroscopic solver, even without reducing computational cost. Moreover, a novel wave-particle adaptation has been proposed, which incorporates the advantages of the stochastic particle method and a deterministic hydrodynamic solver. Multiscale flows at high Mach numbers with large variations of Knudsen number in computational domain can easily be dealt with by the UGKWP.
For high-speed flow simulation, the development of the GKS and the UGKS in the future may focus on these subjects, such as multi-component flow with large numbers of species; non-equilibrium flow physics, including ionization and chemical reactions; and large-scale computations using hybrid CPU-GPU systems.
Author Contributions: Investigation, Y.Z.; writing, original draft preparation, Y.Z. and K.X.; conceptualization, supervision, reviewing, and editing, C.Z. and K.X. All authors have read and agreed to the published version of the manuscript.