Gas‐kinetic unified algorithm for plane external force‐driven flows covering all flow regimes by modeling of Boltzmann equation

The nonequilibrium steady gas flows under the external forces are essentially associated with some extremely complicated nonlinear dynamics, due to the acceleration or deceleration effects of the external forces on the gas molecules by the velocity distribution function. In this article, the gas‐kinetic unified algorithm (GKUA) for rarefied transition to continuum flows under external forces is developed by solving the unified Boltzmann model equation. The computable modeling of the Boltzmann equation with the external force terms is presented at the first time by introducing the gas molecular collision relaxing parameter and the local equilibrium distribution function integrated in the unified expression with the flow state controlling parameter, including the macroscopic flow variables, the gas viscosity transport coefficient, the thermodynamic effect, the molecular power law, and molecular models, covering a full spectrum of flow regimes. The conservative discrete velocity ordinate (DVO) method is utilized to transform the governing equation into the hyperbolic conservation forms at each of the DVO points. The corresponding numerical schemes are constructed, especially the forward‐backward MacCormack predictor‐corrector method for the convection term in the molecular velocity space, which is unlike the original type. Some typical numerical examples are conducted to test the present new algorithm. The results obtained by the relevant direct simulation Monte Carlo method, Euler/Navier‐Stokes solver, unified gas‐kinetic scheme, and moment methods are compared with the numerical analysis solutions of the present GKUA, which are in good agreement, demonstrating the high accuracy of the present algorithm. Besides, some anomalous features in these flows are observed and analyzed in detail. The numerical experience indicates that the present GKUA can provide potential applications for the simulations of the nonequilibrium external‐force driven flows, such as the gravity, the electric force, and the Lorentz force fields covering all flow regimes.


INTRODUCTION
In general, there exist one class of nonequilibrium steady gas flows under the external forces, for example, the atmospheric convection on the planets under the constant gravitational force, 1 and the plasma flows under the electric field force and the magnetic field Lorentz force, 2 which is different from the internal shock flow, the microscale flow, and the reentry vehicle flow. This class of gas flows is driven by the external force ⃗ F, which causes an acceleration or deceleration of the gas molecule and changes the varying rate of molecular velocity ⃗ V in the velocity space. The intensity of the external force ⃗ F affects and embodies the deviation degree from the equilibrium state directly. Besides, in these flows, because of the coupling with particle transport, collision behaviors, and rarefaction effects, extremely complicated nonlinear dynamics will appear along with some peculiar flow features. The nonequilibrium steady gas flows under the external forces have been constantly confronted in various of practical engineering problems, such as the microelectronics, 3 the plasma physics, 4 the hypersonic reentry ionized flow, 5 the miscible displacement processes in porous media flows under gravity, 6 the thermal conduction, 7 the color conductivity problems, 8 the astrophysics, and so on. Hence, the simulations and investigations for this kind of flow problems have been being of great importance and common concerns and interesting issues nowadays.
In terms of the degree of gas rarefaction, the gas flows can be approximately categorized into the four flow regimes, 9 including the continuum flow, the slip flow, the transition flow, and the highly rarefied free-molecule flow. The classification criterion is based on the dimensionless controlling parameter of the flow state, that is, the Knudsen number Kn, which is defined as the ratio of the molecular mean free path and the characteristic physical length scale L. For the continuum flows, the computational fluid dynamics (CFD) techniques based on the Euler or Navier-Stokes (NS) equations have been devoted to obtain convincing and precise simulation results, 10,11 such that the high-resolution well-balanced schemes have been constructed to solve the gravitational fields. 10 For the rarefied nonequilibrium flows with large Kn, due to the advantage of the collision-relaxation evolution feature from the Boltzmann equation, the direct simulation Monte Carlo (DSMC) method based on the particle stochastic statistics theory can present accurate numerical results under sustainable computations. 12 However, for flows in the continuum-transition regimes with small Knudsen numbers, which have been encountered in many practical engineering problems frequently, the NS equations are known to be inadequate due to the invalidation of the continuous fluid hypothesis. 13 In addition, it is quite difficult for the DSMC method to simulate this flow regime due to the limitation to the cell size and the time step required by the method for itself. 14,15 In the past decades, there have been plentiful hybrid methods, which combined NS and kinetic approaches, to be used to model the flows, which are in the continuum-transition regimes or have both continuum and rarefied regimes. 16 These approaches are applied geometrically through the coupling zone, which may sensitively depend on the location of the interface between different methods and hardly be used and developed widely. 17 Hence, it has been becoming more and more necessary and important to develop an integrated simulation method 9,14,18 for high rarefied to continuum flows.
The Boltzmann equation, 19,20 describing the time evolution of the gas molecular velocity distribution function from nonequilibrium to equilibrium state, has been seen to model the flow transport phenomena for the full spectrum of flow regimes and act as the main foundation for the study of complex gas dynamics. 20 The Euler, NS, and Burnett equations in the macroscopic fluid dynamics can be derived from the Boltzmann equation by applying the Chapman-Enskog expansion to the zeroth, first, and second-order approximation of the Maxwellian distribution, respectively. 19 To computably model the nonlinear integral differential equation into partial differential equation, several types of kinetic model equations, resembling to the original Boltzmann equation and concerning the various order of moments, have been presented according to the mass, momentum, and energy conservation laws. 21,22 One of them is the famous Bhatnagar-Gross-Krook (BGK) model, 21 which replaces the collision integro-differential term of the original Boltzmann equation with a simple colliding relaxation model and has been seen to provide a more tractable way to solve comparatively complex gas flow problems. During the last decades, several methods have been developed to simulate the flows based on the kinetic theory of gases theories, for example, the lattice Boltzmann method (LBM), 23 the kinetic flux vector splitting schemes, 24 the asymptotic relaxation schemes, 25 the moment methods, 26 and a series of gas-kinetic BGK-type schemes, such as the gas-kinetic scheme (GKS), 27 the unified gas-kinetic scheme (UGKS), 17 and so on. Specially, the stability of the Boltzmann equation with external forces is studied in Reference 28, and the kinetic models are presented to simulate the homogeneous color diffusion problem in a dilute binary mixture. 29 The LBM is used to compute the fluid motions driven by external forces, 30 and the UGKS is extended to develop well-balanced methods to study the multiple nonequilibrium flow phenomena for the gravitational field. 31 To study aerodynamic problems covering various flow regimes, the gas-kinetic unified algorithm (GKUA) has been applied to the gas flow simulations from the rarefied free-molecule flow to the continuum flow regimes with a wide range of Mach numbers. [32][33][34][35][36][37][38] The GKUA has unified expressions on the molecular collision relaxing parameter and the local equilibrium distribution function, which are associated with the macroscopic flow variables, the gas viscosity transport coefficient, the thermodynamic effect, the molecular power law, molecular models, and the flow state controlling parameter Kn from various flow regimes. In the past, the GKUA has been successfully used to study hypersonic reentry aerothermodynamics around kinds of space vehicles and microscale flows involved in MEMS devices [35][36][37] without the external force. In this article, the GKUA will be extended to develop the computable modeling of the Boltzmann equation with of the external force terms, and the proposed algorithm will be applied to the multiscale nonequilibrium flows under the external force field. This is the first attempt to study and use the direct Boltzmann-type gas-kinetic solver of the GKUA framework in diatomic gas to compute the external force-driven flows covering various flow regimes with a wide range [5 × 10 −5 , 10] of Knudsen numbers. The new DVO technique and numerical quadrature methods, gas-kinetic boundary conditions and numerical schemes are constructed directly to solve the unified Boltzmann model equation in external force term.
Hereby, the remaining part of this article will be structured as follows. The unified Boltzmann model equation for external force-driven flows will be presented in Section 2. Numerical procedures will be constructed in Section 3, including the numerical schemes and the boundary condition procedures for external force-driven flows. Especially, the conservative DVO method (C-DVOM) is first developed and utilized to discretize the velocity distribution function equation into a set of hyperbolic conservation laws with source terms to preserve the conservation of collision term at the discrete levels. Utilizing the new algorithm, some numerical examples are investigated to validate the present algorithm in Section 4, including the Sod shock-tube problem with backward external force, the plane external force-driven Poiseuille flows, and the lid-driven cavity flow under gravitational field, and the effects of the external force and the rarefaction on these flows are discussed as well. Conclusions will be summarized in Section 5.

DESCRIPTION OF THE UNIFIED BOLTZMANN MODEL EQUATION FOR EXTERNAL FORCE-DRIVEN FLOWS
Considering the effects of the external forces, the governing equation of the flow transport phenomena from the full spectrum of flow regimes in the view of statistical mechanics has been transformed from the original Boltzmann equation of the external force terms into the following form: 20 where f is the molecular velocity distribution function, which depends on space ⃗ r, molecular velocity ⃗ V and the time t. ⃗ F is the vector of the external force. Q(f, f 1 ) is the collision term describing the binary collision between particles, which is a complex integro-differential term. The gas transport properties and the governing equations describing macroscopic gas flows can be obtained from Equation (1) or its model equations. 21,22 Applying the BGK model, 21 which has a relaxation collision term and is robust for a wide range of numerical schemes, Equation (1) can be transformed as the following form: where is a constant parameter called the collision frequency, which is determined by the viscosity, and f M is a standard Maxwellian equilibrium distribution function of the following form: 21,27 Herein, n and T denote the density and temperature of gas flow, respectively, R is the gas constant, ⃗ C represents the molecule thermal velocity defined as ⃗ where ⃗ U is the mean flow velocity. D is the number of dimensions, K is the number of internal degree of freedom and ⃗ = ( 1 , 2 , … , K ) T . The relation between K and the ratio of specific heat is: For example, a monatomic gas at one-dimensional problem has K = 2, which is accounted for the motion in y, z direction, and ⃗ = ( V y , V z ) T , while K = 0 for a monatomic gas at three-dimensional problem and ⃗ = 0.
Based on the gas-kinetic Shakhov model, 22 the local Maxwellian f M in the BGK Equation (2) can be replaced by the local equilibrium distribution function f N of the following form to obtain the correct Prandtl number and the properties of viscosity and heat conduction: temperature exponent of the coefficient of viscosity transport, ref is the reference mean free path, and other parameters have the same meanings as mentioned before. Generally, we set the free-stream values as the reference values. Obviously, the local equilibrium distribution function f N and the collision relaxing parameter can be determined by the macroscopic flow variables, the viscosity transport coefficient, the thermodynamic effect, the molecular power law models, and the flow state controlling parameter from various flow regimes.
Following the characteristic variables in Table 1, the nondimensional unified simplified velocity distribution function equation with the external force term is presented in the Cartesian coordinates as follows: with The constant is introduced as the dimensionless adjustment parameter in order to unify different relations of the dimensionless density, temperature, and pressure in some previous literatures. Generally, = 1, and Equation (9) can be transformed into the same forms in References 32-38. All of the variables in Equation (9) and in all following equations are nondimensionalized by using Table 1, if it has not been mentioned specifically.
Equation (9) describes the time dependence of the molecular velocity distribution function on the position space ⃗ r and velocity space ⃗ V along with providing the statistical description of the gas flow from the level of the kinetic theory of gases. All of the macroscopic flow variables of gas dynamics under consideration, such as the density of the gas n, the flow velocity ⃗ U, the temperature T, the pressure P, the viscous stress tensor , and the heat flux vector ⃗ q, can be integrated from the moments of the velocity distribution function over the molecular velocity space:

GAS-KINETIC NUMERICAL PROCEDURES IN SOLVING BOLTZMANN MODEL VELOCITY DISTRIBUTION FUNCTION EQUATION
In this section, the numerical procedures for the unified velocity distribution function Equation (9) will be presented in detail for external force-driven flows at the first time, including the C-DVOM, the improved Euler method for the time evolution, the NND-4(a) splitting scheme for the convection in the physical space, and the forward-backward MacCormack predictor-corrector method for the convection in the molecular velocity space.

Conservative DVO method
Due to the multidimensional complexity of the molecular velocity distribution function f , the discrete velocity ordinate method (DVOM) [32][33][34][35][36][37][38] in the gas-kinetic theory has been introduced to get rid of the continuous dependency of g and h on the velocity space. The DVOM is based on the representation of functions by a set of DVO points and coincides with the evaluated points of the quadrature rule for the computation of the macroscopic flow moments of the velocity distribution function. Herein, we use the two-dimensional case as an example, that is, D = 2. Utilizing the DVOM in the V x and V y velocity space, Equation (9) for the two-dimensional flows can be transformed into hyperbolic conservation forms at each of the DVO points V x , V y , ( = 1, 2, … , N ; = 1, 2, … , N ) as the following forms: with where g , (x, y, t), h , (x, y, t), G N , (x, y, t), and H N , (x, y, t) correspond to the values of g, h, G N , and H N at the DVO points (V x , V y ), and the subscript , represents the discrete velocity grid indexes in the V x -direction and V y -direction, respectively.
Once the discrete velocity distribution functions g , (x, y, t) and h , (x, y, t) are solved, the macroscopic flow moments at any time in each point of the physical space can be obtained and updated by the appropriate discrete velocity quadrature method. For example, using the discrete velocity quadrature method, the gas density by Equation (10a) can be calculated as: where the coefficients A , A and V x , V y are the weights and the DVOs, respectively, which are determined by the sets of the quadrature rules. In our previous studies, 32-36 the multiple quadrature rules have been adopted for vast flow simulations successfully, including the Gauss-Hermite, the composite Newton-Cotes, the multisubinterval Gauss-Legendre and the Gauss-Chebyshev quadrature rules, whose weights and DVOs have been given and omitted here for the conciseness. However, due to the errors caused by the accuracy problems of the computational schemes and the quadrature rules, the macroscopic flow variables obtained by the numerical procedures are not highly accurate. Increasing the precision of the schemes, or enlarging the integral domain of the discrete velocity space, may cure this problems, which increases computing demand and has high computational cost. In order to analyze and solve this problem, the C-DVOM 41,42 is utilized in this article, which is on the basis of the conservative conditions of the Boltzmann model equation and increases the accuracy of the computational macroscopic flow variables in an acceptable computing cost scope.
It is well known that the conservation laws for mass, momentum, and energy of the gas can be obtained by multiplying (1) by collision invariants = )T and integrating over the entire molecular velocity space: where Q(f, f 1 ) is the exact collision integral term. It is natural to require that this property be maintained at the discrete level by the numerical method. Since the differential parts in the left side of the exact Equation (1) and the Boltzmann model Equation (9) on the molecular velocity distribution function are the same as the external force term, the approximation condition means that the first few moments of the exact collision integral Q(f, f 1 ) and the collision relaxation term in the right side of Equation (9) have the following relations: Recalling the expressions of f N , one can obtain the following integration results precisely: However, in the actual numerical simulations, the molecular velocity space is interceptive and discretized, which leads to the moment of f N calculated by the quadrature rules. Besides, f N is constructed by the macroscopic flow variables , which were updated and calculated using the numerical quadrature rules so that the DVOM might not satisfy the constraints (14). If the quadrature rule used in the DVOM is of the rth-order accuracy, it can obtain the following results, Combining Equation (10) and the above results (16), we have Theoretically, the value of | | should be zero. However, | | ≠ 0 due to the computational errors caused by the numerical integration or the schemes. From Equation (17), it can be found that the value of | | will decrease following the increase of the quadrature order when the state of the incoming flow remains unchanged, while the higher order quadrature needs more DVOs, which brings huge and unbearable calculation amount. On the other hand, for simulating the continuum flow, that is, Kn → 0, it can be seen that | | ≫ 0 even when we apply higher order quadrature rules. Hence, finding a method to minimize | |, which is not dependent on the Knudsen number Kn and the integral accuracy, will increase the computational efficiency and enhance the robustness of the algorithm.
Considering the heat flux ⃗ q has appeared in the expression of the local equilibrium distribution function f N , we introduce new variables , and can calculate the following relation: In combination with Equations (10f), (14), and (18), we have the following system of nonlinear equations: Applying some type of quadrature rule, Equation (19) can be rewritten as follows: where (20) can be solved by using the Newton-Raphson method for nonlinear system equations with multivariables, whose iterative formula can be expressed as the following form: where the matrix M, defined as M = H∕ W, is the Jacobi matrix that contains the differential and numerical integration calculations, and is difficult to obtain its exact expressions. Recalling Equation (6) and using the nondimensional variables, the local equilibrium distribution function G N and H N have the following identity integral formulas: Correspondingly, Equation (20) for the two-dimensional flows can be rewritten specifically as: Especially, the approximate Jacobi matrix  (2) and its inverse matrix  −1 (2) in the two-dimensional flows have the following forms:

Numerical scheme
In order to treat arbitrary geometry configuration, the body fitted coordinate is introduced. In the transformed coordinates ( , ), Equation (11) can be transformed into hyperbolic conservation equations as follows with nonlinear source terms at each of discrete space grid points: with where J is the Jacobian of the general transformation, whose metric forms are denoted as and U, V are the so-called contravariant molecular velocity defined as The Jacobian coefficient matrices A =  ∕  and B = ∕  of the transformed Equation (25) are diagonal and have the real eigenvalues a = U and b = V.
The time-splitting numerical method 38 is adopted to solve Equation (25) simultaneously in physical space and velocity space as the first attempt. The discrete velocity distribution function Equation (25) can be decomposed into colliding relaxation equations with nonlinear source terms and convective equations in the physical and velocity space, respectively: The above equations can be solved by the following subsequent four steps: (1) Using the improved Euler method to solve the nonlinear source terms Equation (28a), which represents the variation of the colliding relaxation: One can prove the above difference equation  n+1 = L s (Δt) n is compatible with Equation (28a) in the second-order accuracy.
(2) The convective Equations (28b) and (28c) in the -and -direction can be numerically solved by using the NND-4(a) scheme 43 based on primitive variables with second-order accuracy in time and space. For example, the convective Equation (28b) in the -direction can be solved as: where One can prove the above difference equation  n+1 = L (Δt) n is compatible with Equation (28b) at the (Δt 2 , Δ 2 ) accuracy. (3) For Equations (28d) and (28e), the improved Euler method is applied to discretize the time derivative term ∕ t.
As DVOM have been applied to discretize the molecular velocity spaces, leading to the finite difference method as the natural dealing approach for these two equations. Hence, various finite difference methods can be considered to wield to discretize the convective terms ∕ V x and ∕ V y in Equation (25).
However, the meshes from the nodes of the most quadrature rules are nonuniform including the Gauss-Hermite, the multisubinterval Gauss-Legendre and the Gauss-Chebyshev quadrature rules, and only uniform in the composite Newton-Cotes rules. The nodes of these quadrature rules in one case are shown in Figure 1. Hence, the finite difference methods on nonuniform meshes of the discretized velocity space should be presented initially in the procedures for the external force-driven flows with serious nonequilibrium effects, considering the universality in the future.
The central finite difference method is not suitable to be applied here due to the nonuniform meshes. The forward or backward finite difference can be used to approximate the first-order derivative discretely. Based on the explicit Mac-Cormack operator splitting concept, the forward-backward MacCormack predictor-corrector method, where the forward difference method is used in the predictor step and the backward difference method is applied in the corrector step, is utilized for solving Equations (28d) and (28e).
For the convective Equation (28d) in the V x -direction, the specific expressions are given as follows: For the convective Equation (28e) in the V y -direction, one can calculate as: For the boundary grid points of the discretized velocity space, the linear interpolation is applied to obtain the variation. For example, in the predictor step of Equation (31), the forward difference is used, which makes t  * at the point = 1 incomputable. Assuming the variations t  * at the points = 1 to = 3 are linear to the nonuniform mesh spacing, as the following relation: the variations t  * at the point = 1 can be expressed by Similarly, in the corrector step, the variations t  * * at the point = N can be expressed by For Equation (32), the same approaches are applied for the variations at the points = 1 and = N, whose specific expressions are omitted here. One can prove the difference equation In conclusion, the finite difference scheme for the hyperbolic conservation equations with nonlinear source terms (25) can be written as follows: Considering simultaneously proceeding on the convective movement and colliding relaxation in real gas, the computing order of the previous and subsequent time steps is interchanged to eliminate and control the numerical discretized errors caused by the sequential direction of splitting computation. The second-order finite difference scheme directly solving the discrete velocity distribution functions (25) is constructed as: In addition, the time step size Δt in the computation should be controlled only by the stability condition of the scheme: where CFL is the adjusting coefficient of the time step in the scheme.

Boundary condition procedures for the evolution of molecular velocity distribution function
There exist four type boundary conditions, that is, the wall surface condition, the entrance and exit condition, the periodic condition, and the symmetrical condition, which should be considered in this article. Since our present algorithm evaluates the time evolution of the molecular velocity distribution function explicitly at each of discrete grid points from the physical and velocity space, all kinds of boundary conditions should be numerically implemented by acting on the velocity distribution function directly, instead of only using the macroscopic flow variables.

Solid wall boundary
In order to specify the interaction between the gas molecules and the solid wall surface, it is assumed that the gas molecules striking the surface are subsequently emitted with the Maxwellian velocity distribution partly accommodating to the wall temperature T w and the velocity ⃗ U w = (U w , V w ), coupling with the specular reflection, which is characterized by the accommodation coefficient . Define C n = ⃗ C ⋅ ⃗ n, and ⃗ n is the outward unit vector normal to the wall surface.
When C n > 0 in the wall cells of two-dimensional flows, for the Maxwell's fully accommodation boundary, the velocity distribution functions reflected from the wall surface can be written as follows in discretized form: where n w represents the number density of molecules diffusing from the solid surface. n w can be determined by the mass balance condition on the surface, which assumes that the mass flux normal to the wall surface for the incident molecules equals that of the reflective molecules at any time, of the following form: From Equation (40), we have where C − n = (C n − |C n |)∕2. Once the number density of molecules diffusing from the solid surface n w has been calculated, the reflected velocity distribution functions can be obtained. For the specular reflection wall boundary, it has the following expression, From the above, the reflected velocity distribution functions can be obtained by We set = 1 at the computations in the remaining part of this article, if it has not been mentioned specifically. When C n ⩽ 0, the discrete velocity distribution functions at the wall cells are solved by using second-order upwind-difference approximations from the adjacent grids.

Entrance and exit boundaries
For the entrance and exit boundary conditions, it should be treated using characteristics-based boundary conditions, due to the finite distance among the entrance boundary, the exit boundary, and the body. The distribution functions for outgoing molecules through the entrance or exit boundary are numerically discretized by the second-order difference approximation in terms of the condition, which is in accord with the upwind nature of the interior point scheme.
For molecules incoming from outside, the following conditions are used for the calculation to the boundary cells: when it is along upstream boundary ahead of the body, the distribution functions of the ingoing molecules from upstream are set as an equilibrium distribution with prescribed free stream properties; and it is assumed that there is no gradient along the outward direction for the distribution function in the downstream boundary.

Periodic boundary
For the periodic condition, it assumes that the entrance and exit boundaries are connected, which means the molecular velocity distribution functions of the outside fictitious nodes are copied from the functions of the interior nodes, which is shown in Figure 2.

NUMERICAL SIMULATIONS AND DISCUSSION
In this section, we present results of some selected computational examples, including the Sod shock-tube problem with backward external force, the external force-driven Poiseuille flow, and the lid-driven cavity flow under gravitational field, to validate the GKUA for flows under external force terms covering the whole range of flow regimes.

Sod shock-tube problem with backward external force
The first test case is the standard Sod shock-tube problem coupled with backward external force, which has a diaphragm located at x = 0.5 and contains two different constant equilibrium states at t = 0. Following the problem setup in References 10,31, the computational domain is set as [0, 1], and the initial conditions are given by: The backward external force is set as F x = −1.0 in the negative x-direction. In Reference 31, the UGKS is used to simulate this problem in monatomic gas with = 5∕3. Herein, the simulation uses diatomic gas with = 1.4 and performs on a uniform mesh with N = 1001 until final time t = 0.2. The Prandtl number is set as Pr = 1, and other parameters are set as = 1, = 0.5, and CFL = 0.9. The modified Gauss-Hermite quadrature rule with 32 DVO points is applied for the simulation to evaluate the macroscopic flow moments over the velocity space. The specular reflecting boundary condition is applied at the left and right ends. In this case, the Knudsen number is defined as the ratio of the reference mean free path and the length of the tube, which can be varied with the reference number density. The time step is decided by Equation (38) only from the stability condition of the scheme, which can be found that a larger time step can be set in solving the flow with higher Knudsen number, leading to the less computational time. Figure 3 presents the density, flow velocity, energy, and pressure profiles obtained by the present GKUA at different Knudsen numbers of Kn = 5 × 10 −5 to 10, as well as the Euler solver results with the high-order well-balanced WENO scheme in Reference 10. It can be observed that the GKUA results at Kn = 5 × 10 −5 match well with the Euler solver results, which demonstrates the feasibility and accuracy of the present algorithm in solving the continuum gas flows under the impacts of the external force. Besides, it is shown that the distinct shock and rarefaction wave is gradually forming following the decrease of the Knudsen number. There is no shock and expanding wave, but strong disturbance for the fully rarefied flow related to the state Kn = 10.0 only exists. The above-mentioned phenomena are only the representation of the notable differences between continuum and rarefied flows, which are aroused by the difference of molecular transport capability from various flow regimes. For the case F x = −1.0, the results of computed density, velocity, and temperature profiles (denoted as symbols) by using the DVOM and C-DVOM are shown in Figures 4 to 6, respectively. The solid line denotes the solution of the finest 201 velocity nodes in velocity space, that is, considered as the reference solutions. It can be observed that the C-DVOM results agree very well among the cases of different velocity nodes, while the cases of 31 velocity nodes are not accurate for the DVOM results. This phenomenon is caused by that the velocity grids are not fine enough to accurately evaluate the moment integrals of the distribution function for flow variables and the compatibility condition of the discrete collision integral is not well preserved by the nonconservative calculations of collision integral. Although the C-DVOM contains steps of solving nonlinear system of equations and calculating the matrixes, less number of velocity nodes are used in the C-DVOM, leading to less computation time and more accurate results than the DVOM. Hereby, the C-DVOM is suitable to be utilized in the present algorithm.
Comparing with the simulation results of the original Sod shock-tube problem without the external force, as shown in References 32-38, it can be observed that the density distribution is pulling toward the left direction, and negative velocity appears in some regions due to the negative external force. It also can be indicated that external force causes more serious nonequilibrium flow effects, and the C-DVOM can improve the robustness of the present GKUA by satisfying the conservation condition in the real-time computational process with less DVO points.

External force-driven Poiseuille flow
The Poiseuille flow, which is confined between two rigid parallel plates that are stationary and act as thermal reservoirs, displays complex fluid mechanical phenomena with multiple scales and covers the most basic and important flow properties, even under simple geometric conditions. 44 Generally, the Poiseuille flow is driven by the pressure gradient, which has been studied extensively by both CFD and the kinetic approaches. 27,45,46 There exist another type of Poiseuille flows which is driven by the external force in the flow direction, 45,46 for example, the gravity. The external force-driven Poiseuille flow can be confined as follows: considering an ideal gas between two parallel infinite plates at rest with a common uniform temperature, a steady unidirectional flow of the gas is caused between the plates, when the gas is subject to a uniform external force in the direction parallel to the plates. The planar external force-driven Poiseuille flow exhibits exceptional features when the gas is rarefied, for example, a nonconstant pressure profile and a bimodal temperature distribution with a local minimum lying at the center of the channel. 26 In order to understand these phenomena, many analyses have been done. These features were originally predicted in Reference 47 using a perturbative solution to the BGK model and further confirmed in Reference 46 using DSMC method. Both the Navier-Stokes-Fourier (NSF) and Burnett equations are unable to capture the observed bimodal temperature profile, even in the slip regime, although the Burnett equations do recover the nonconstant pressure profile. 48 In addition, the temperature minimum at the center is explained only through the super-Burnett solution or the kinetic theory. 27 This flows can be a suitable and referential example to verify the accuracy and efficiency of the present GKUA with the external force terms.
The parameters are given in Reference 45 as follows. The simulated fluid is a hard sphere gas with particle mass m = 1 and diameter d = 1. The Boltzmann's constant is taken as k = 0.5, and the gas constant = 5∕3 for a monatomic gas. At From Figure 7, the present GKUA results are able to reproduce the anomalous features predicted in the DSMC simulations fairly accurately, such as the bimodal temperature, nonconstant pressure, and heat flux profiles, while the slip NS results are quite different from other two results. It can be seen that the density is lower in the middle part of the channel, and increases sharply in the parts near the up and down walls. In comparison to the DSMC data, the temperature at the center of the channel is underpredicted slightly by GKUA, and the lowest temperature is arisen at the up and down walls. The temperature shown in Figure 7B includes energy components in the three spatial dimensions, whose definition can be seen in Reference 49. The nondimensionalized pressure distribution along the perpendicular wall direction is shown in Figure 7C, from which we can find that the GKUA results' tendency agrees with the DSMC results basically, while the slip NS solver obtains different results and appears as constant. The pressure is at the bottom in the middle part of the channel, and increases away from the center symmetrically. Figure 7D gives the distribution of the shear stress along the perpendicular wall direction. It can be found that all three methods obtain almost the same results, which can demonstrate the correctness of the shear stress profile and the feasibility of the slip NS solver in predicting some flow variables of the continuum-transition flow regime such as a fairly large Knudsen number, that is, Kn = 0.1. The theoretical reason for the slip NS solver to correctly solve the linear distribution of the shear stress along with the perpendicular wall direction is that the shear stress roots in the first-order crossderivatives derived from flow velocity vs position space to be solved accurately in the NS equation, although it is more difficult for the NS solver to solve the macroscopic flow variables of higher order derivatives. And it is difficult and invalid for the slip NS method to solve more nonequilibrium flow distributions. The profiles of the heat flux in the axial and lateral components are shown in Figure 7E,F, respectively. It can be found that the GKUA results match the DSMC simulations well, only existing slight deviation at part near the wall. The heat flux in the axial component is in symmetrical distribution, which gets the minimum value at the center of The velocity profiles along the perpendicular wall direction for the first case of the force-driven Poiseuille flow are shown in Figure 8, from which we can find that the axial velocity profiles obtained by these three methods are almost the same, which are parabolic, while the crosswise velocity profiles have obvious distinction. The crosswise velocity profile obtained by the GKUA is presented as a sine-like wave, which is about zero at the center, the up and down wall, and reaches the maximal absolute value at the positions H∕4 and 3H∕4. The sinusoidal-like distribution of lateral velocity V is due to zero momentum in this Y -direction and nonlinear distribution, and needs to maintain a similar distribution to the lateral heat flux distribution shown in Figure 7F. The maximal absolute value of the crosswise velocity is not exceeding 6 × 10 −5 , whereas the DSMC results are almost zero, with wider margin of oscillation. We think Furthermore, the reason on the deviation of the bimodal part of the temperature is that the rarefied effects of temperature jump and velocity slip nearby wall surface are revealed by the present GKUA different from the DSMC. The higher temperature jump and the zero of lateral velocity from Figure 8B are resulted from the tremendous statistical fluctuations of the DSMC method in simulating the weak flow of the external force-driven Poiseuille flow of Kn = 0.1 and F x = 6.868 × 10 −5 , however, the lower temperature jump and the weak waveform structure of lateral velocity are obtained by the GKUA, which brings about the slight difference of the bimodal distribution of the temperature at the center of the channel existed in the results of the GKUA and DSMC, although the slip N-S solver cannot solve this nonlinear flow phenomena.
The temperature, the pressure, the axial velocity, the shear stress, the axial heat flux, and the lateral heat flux profiles for the second case of the force-driven Poiseuille flow are shown in Figure 9, where the black solid line represents the numerical simulation results obtained by the present GKUA, the blue delta symbols stand for the results obtained by solving the 13 moment equations, 26 and the red circle symbols represent the results obtained by solving the 26 moment equations, 26 which are presented for comparison. The varying rules and tendencies of all flow variable distributions in this case are the same as the ones in the first case, including the bimodal temperature, nonconstant pressure, and heat flux profiles. In addition, the simulated results obtained by the present algorithm are in close agreement with the 26 moment results, which verify the feasibility of the present GKUA in solving the Poiseuille channel flows from the near-continuum slip flow regime of Kn = 0.05 under the external force field. By comparing the varying features of nonconstant pressure profiles shown in Figure 9B, it is clear that the 26-moment method better reveals the nonequilibrium flow phenomena than the 13-moment method, while the GKUA most clearly represents the nonconstant pressure distribution among the three approaches. By the simulation shown in Figures 7 to 9, the validity and accuracy of the present GKUA is verified in solving external force-driven flows covering various flow regimes, and the nonequilibrium flow phenomena with different Knudsen numbers from near-continuum and transitional flow regimes are revealed and caused under constant external forces.

Lid-driven cavity flow under gravitational field
The lid-driven cavity flow problem is a complex system including boundary effect, shearing structure, heat transfer, and nonequilibrium thermodynamics, which is often considered to be a benchmark problem for validation of algorithm. Gaseous flow and heat transfer in square cavity enclosures that represent vacuum-packaged MEMS devices have been studied in Reference 50, and the driven cavity problem was also examined with the aid of the LBM in The stationary profile contours of density, temperature, tangential x-direction, and vertical y-direction velocity in the lid-driven cavity flow without the external force are shown in Figure 10. Moreover, the distribution contours of these four flow variables in the cavity flow with F y = −1.0 and −2.0 are shown in Figures 11 and 12, respectively. In addition, the streamlines and velocity vector fields of each case are shown at their corresponding velocity contours, and the heat flux stream traces in each case are plotted overlaid on each temperature contour as well.
In the case without the external force, as shown in Figure 10, it is observed that a significant drop in density is seen toward the top left corner of the cavity, whereas a significant rise in density is observed toward the top right corner of the cavity, as well as the contour of the temperature profile. This can be attributed to the expansion of gas at the top left corner as it is driven by the moving lid, whereas it is compressed at the opposite end. In addition, the modes of the heat flux propagations are quite unordered in this case, including the cold to cold regions, the hot to hot regions, the cold to hot regions, and the hot to cold regions, which can be seen in Figure 10B. The tangential x-direction velocity contours are symmetric along the vertical centerline, while the vertical y-direction velocity contours are antisymmetric along this line, as shown in Figure 10C,D.
However, in the case with the external force, the movement of upper wall and the external force are two driving sources for the flow motion in the cavity system, which lead to quite different distributions of the flow variables. From Figures 11 to 12, it can be observed that the density distribution changes significantly along the vertical direction, and a great scale drop in temperature is seen toward the top left corner of the cavity, whereas the top right corner and the bottom part of the cavity get rise in temperature. The heat propagation is also different from the case without the external force, which is from the top (cold) to the bottom (hot) regions, as shown in Figures 11B and 12B. This phenomenon is against Fourier's law and also different from the nonequilibrium heat flux due to the rarefaction effect only, as presented in Reference 12. The results show that the heat flux is transferred from the cold to the hot region due to the external force effect, or the gravitational effect, which may be applied to explain the gravity-thermal instability phenomena arising in the field of astrophysics. The tangential x-direction velocity contours are also symmetric along the vertical centerline, while the vertical y-direction velocity contours are not antisymmetric along this line, which has a large upward or downward vertical velocity. Comparing Figures 11 and 12, it can be found that the distinction of the density and temperature profiles between the top and bottom regions of the cavity becomes larger with the increase of the external force. Besides, the region of the temperature drop gets extended following the increase of the external force. The velocity components, U and V, in the X-and Y -directions at the vertical and horizontal central lines of the cavity, respectively, are shown in Figure 13. The mark lines represent the results obtained by the present GKUA with different external forces, and the mark symbols in Figure 13A are from the results given in Reference 31 for comparison with uniform cells in the physical space. It is observed that the present GKUA results match well with Reference 31 results for different external forces, which verify the feasibility of the present GKUA in computable modeling of the Boltzmann-type velocity distribution function equation to simulate the lid-driven cavity flow with the downward gravitational external forces. Nevertheless, some difference of the tangential velocity U profile exists nearby the highly nonequilibrium region of the moving upper wall for higher external force of F y = −2. The reasons causing this disparity may be further investigated in the future work.
From Figure 13A, it can be seen that the U-velocity profile is characterized by the velocity of the upper wall. It is interesting to note that the nondimensional slip velocity becomes to be smaller from 0.4 to 0.53 with the increase of the external force from 0 to 2 for the three cases of varying external force, this is due to the external forces in the negative Y -direction although there is a tangential velocity of U w = 0.15. In the case without the external force, the V-velocity profile is symmetrical about a vertical line crossing the cavity center. However, with the increment of the external force, the V-velocity profile is no longer symmetrical and is shifted to the right, which is consistent with the vortex center shifting toward the right of the cavity at larger external force.
The impact of rarefaction on the flow and heat transfer in the lid-driven cavity flow with downward external force can be assessed by carrying out numerical simulations for a wide range of Knudsen numbers from the continuum to the free-molecule flow regime. Only selected results will be shown here, that is, Kn = 0.0005, 0.001, 0.01, 0.5, 5.0, and 10.0. In this part of the study, the lid velocity is fixed at U w = 0.15, and the external force is fixed at F y = −1.0. Figure 14 shows the stationary profile contours of the temperature in the lid-driven cavity flow with external force at different Knudsen numbers, as well as the heat flux traces in each case. From Figure 14, the core regions of the observed fall and rise in temperature at either end of the cavity enlarge as the Knudsen number increases. At Kn = 0.0005, which represents the continuum flow regime, the heat flux propagations are from hot to cold regions, as shown in Figure 14A. While at Kn = 0.001, the direction of heat transfer becomes astatic, which include the hot to cold and cold to hot regions, and there exist one source point and one sink line in the heat transfer, as shown in Figure 14B. With the increase of the Knudsen number, from Figure 14C with Kn = 0.01, the sink line disperses into a series of downward traces, and the source point shifts up and disperses into lines. In the rarefied and transitional flow regimes, the direction of heat transfer is from the cold region to the hot region, as shown in Figure 14D-F. Besides, for higher values of Knudsen number, temperature variations are observed throughout the cavity. With increasing Kn, the region of expansion cooling on the heat transfer mechanism increases, as T < T w in a larger region near the left wall. For the case Kn = 5.0, the gas temperature is less than the wall temperature in the entire region near the left wall. This governs the heat transfer for the flows of higher Knudsen numbers. The heat transfer plots indicate that with an increase in rarefaction, both expansion cooling and downward gravitational external force play a greater role in determining the heat transfer characteristics, as the heat flux vectors become more predominantly from the stationary top left wall to the bottom right wall. Figure 15A shows the tangential U-velocity component along a vertical line crossing the cavity center as a function of Kn. It can be observed that the velocity slip increases with any increase in the Knudsen number. The vertical V-velocity component along a horizontal line crossing the center of the cavity at different values of Kn is shown in Figure 15B. The V-velocity is expected to be close to zero at both the left and right stationary walls of the cavity for flow in the continuum regime. And as observed in the previous description, the V-velocity profile is no symmetrical and is shifted to the right with the nonequilibrium wave-like phenomena, which is consistent with the vortex center shifting toward the right of the cavity due to the effects of the external force. Besides, with the decrease of Knudsen numbers, the V-velocity slips at the walls become smaller, while the maximum absolute values of tangential U-velocity profile become larger, in which the rarefaction effect weakens, and the viscous transport and convective diffusion effects strengthen and approach to the dominant position of the whole flow field. An interesting discovery is that when the flow change from the rarefied transition to continuum flow regimes, the maximum absolute values of V profile dive dramatically. This may be caused by the fact that the rarefied effect dominates gradually the whole flow field and needs larger computational region for the Knudsen numbers greater than 0.01. On the other hand, the viscous dissipation of the gas and viscous heat generation increase to consume part of the kinetic energy from instable velocity variations in the y-direction with further decrease of the Knudsen numbers than 0.01. The more reasons causing these phenomena may be investigated by enlarging computational domain for higher Knudsen number flow and developing the present GKUA to solve three-dimensional lid-driven cavity flow under gravitational field in the future work as well.

CONCLUSION
In this article, based on the kinetic theory and the previous works, the GKUA has been developed into new version for some nonequilibrium steady gas flows under the external forces, covering the rarefied free-molecule flow to the continuum flow regimes. The elements of the unified Boltzmann model equation have been given briefly by introducing the gas molecular collision relaxing parameter and the local equilibrium distribution function integrated in unified expression with the flow state controlling parameter, the macroscopic flow variables, the gas viscosity transport coefficient, the thermodynamic effect, the molecular power law, and molecular models for external force-driven flows covering a full spectrum of flow regimes. Utilizing the conservative DVOM, the governing equation has been transformed into the hyperbolic conservation forms at each of the DVO points. Their corresponding numerical schemes have been constructed in two-dimensional detail, including the improved Euler method for the time evolution, the NND-4(a) splitting scheme for the convection in the physical space, and the forward-backward MacCormack predictor-corrector method for the convection in the molecule velocity space.
Applying the new version of the GKUA, some numerical examples have been investigated to validate the present algorithm, including the one-dimensional Sod shock-tube problem with backward external force, the plane external force-driven Poiseuille flows, and the two-dimensional lid-driven cavity flow under gravitational field. Comparing the results presented in the previous articles, it has been observed that the simulation results obtained by the present algorithm match well with them, which demonstrates the validity and accuracy in solving external force-driven flows covering various flow regimes. These results illustrate that the new version of the GKUA can provide valuable insight and potential applications into all aspects of the nonequilibrium steady gas flows under the external forces, such as the gravity, the electric field force, and the Lorentz force.
Furthermore, through comparison, it has been found that the conservative DVOM is more suitable than the original DVOM to be utilized in the present GKUA, due to the conservative calculations of the collision integral of the Boltzmann equation.
In addition, some anomalous features in these numerical examples have been observed and investigated in detail by the current algorithm, which can be summarized as follows: (1) For the one-dimensional Sod shock-tube problem with backward external force, it has been revealed that the density distribution is pulling toward the left direction, and negative velocity appears in some regions due to the negative external force. It also has been found that the distinct shock and rarefaction wave is gradually forming following the decrease of the Knudsen number, and there exist no shock and expanding wave, but strong disturbance for the fully rarefied flow. (2) For the plane external force-driven Poiseuille flows, the present GKUA is validated with some comparative confirmation of the DSMC, slip NS solver, moment methods, and so on in solving external force-driven nonequilibrium flow phenomena with different Knudsen numbers of [0.05, 0.1] from near-continuum and transitional regimes. Some anomalous features have been captured in high resolution and analyzed by the present GKUA, including the bimodal temperature, nonconstant pressure, and heat flux profiles. (3) For the lid-driven cavity flow under gravitational field, after the present GKUA is verified and confirmed with the UGKS in solving the near-continuum lid-driven cavity flows of Kn = 0.01 with/without the external force, the GKUA has been applied in solving and analyzing the effects of rarefaction, external force, viscous transport, and convective diffusion on the flow variable distribution and heat transfer in the lid-driven cavity flow for a wide range of Knudsen numbers of Kn = 5 × 10 −5 to 10 covering all flow regimes. It has been indicated that the heat propagation is from the cold to the hot regions in rarefied transitional flow regimes, due to the dual effects of the movement of upper wall in the tangential direction and the external force in the negative Y -direction. Besides, it has been revealed that with the increase of the external force, the distinction of the density and temperature profiles between the top and bottom regions of the cavity becomes larger. In addition, with an increase in rarefaction, both expansion cooling and gravitational external force play a greater role in determining the heat transfer characteristics, as the heat flux vectors become more predominant from the stationary top left wall to the bottom right wall.
As this work is only the beginning of a gas-kinetic numerical study of external force-driven flows by directly solving the Bolzmann-type velocity distribution function equation in the frame of the present GKUA, further investigations on the kinetic models for three-dimensional complex real gas flows, involving internal energy and external forces, and so on, need to be carried out in more detail.