Unlocking the Key to Accelerating Convergence in the Discrete Velocity Method for Flows in the Near Continuous/Continuous Flow Regimes

How to improve the computational efficiency of flow field simulations around irregular objects in near-continuum and continuum flow regimes has always been a challenge in the aerospace re-entry process. The discrete velocity method (DVM) is a commonly used algorithm for the discretized solutions of the Boltzmann-BGK model equation. However, the discretization of both physical and molecular velocity spaces in DVM can result in significant computational costs. This paper focuses on unlocking the key to accelerate the convergence in DVM calculations, thereby reducing the computational burden. Three versions of DVM are investigated: the semi-implicit DVM (DVM-I), fully implicit DVM (DVM-II), and fully implicit DVM with an inner iteration of the macroscopic governing equation (DVM-III). In order to achieve full implicit discretization of the collision term in the Boltzmann-BGK equation, it is necessary to solve the corresponding macroscopic governing equation in DVM-II and DVM-III. In DVM-III, an inner iterative process of the macroscopic governing equation is employed between two adjacent DVM steps, enabling a more accurate prediction of the equilibrium state for the full implicit discretization of the collision term. Fortunately, the computational cost of solving the macroscopic governing equation is significantly lower than that of the Boltzmann-BGK equation. This is primarily due to the smaller number of conservative variables in the macroscopic governing equation compared to the discrete velocity distribution functions in the Boltzmann-BGK equation. Our findings demonstrate that the fully implicit discretization of the collision term in the Boltzmann-BGK equation can accelerate DVM calculations by one order of magnitude in continuum and near-continuum flow regimes. Furthermore, the introduction of the inner iteration of the macroscopic governing equation provides an additional 1–2 orders of magnitude acceleration. Such advancements hold promise in providing a computational approach for simulating flows around irregular objects in near-space environments.


Introduction
The Boltzmann equation, being free from the limitations of continuum assumption, provides a versatile tool for delineating the molecular transport phenomena that form the bedrock of intricate gas dynamics studies [1,2].However, the original collision term in the Boltzmann equation involves complex fivefold integration in both molecular velocity It is evident in the multiscale approaches that the Boltzmann-BGK equation governs the solution in the rarefied flow regime, while the macroscopic governing equation takes precedence in the near-continuum and continuum flow regimes.Building upon this observation, an improved discrete velocity method (IDVM) was developed by integrating the local solution of the Boltzmann-BGK equation into the calculation of the macroscopic flux [34].In the IDVM, the calculation of the numerical flux for the Boltzmann-BGK equation remains the same as in the conventional DVM to preserve its inherent simplicity.Additionally, by applying the LU-SGS method to address both the Boltzmann-BGK equation and the macroscopic governing equation, a fully implicit IDVM was developed and verified in all flow regimes.
To achieve faster convergence in the near-continuum and continuum flow regimes, the concept of inner iteration was recently introduced into the macroscopic governing equation.This approach allows for a more accurate prediction of the equilibrium state at the new time level, thereby accelerating the convergence rate.Several examples of this technique include the general synthetic iterative scheme (GSIS) [35], the multi-prediction implicit scheme [36], the IDVM with inner iteration [37], and the general implicit iterative method for UGKS [38].Due to the significantly fewer number of conservative variables in the macroscopic governing equation compared to the discrete velocity distribution functions in the Boltzmann-BGK equation, the increase in computational cost for each outer loop iteration is minimal.Consequently, the overall computational cost can be substantially reduced by reducing the total number of iterations required for the outer loop.
Although several multiscale approaches have been developed to accelerate the convergence in DVM, a systematic comparative study in this field is still lacking.In this work, we aim to address this gap by focusing on the discrete velocity Boltzmann-BGK equation and undertaking a comparative analysis of three distinct multiscale schemes for its solution.The first scheme is a multiscale semi-implicit DVM (DVM-I).This approach utilizes the local discrete characteristic solution to calculate the numerical flux, addressing the limitation of mesh size being smaller than the mean free path of molecules.Additionally, it adopts a semi-implicit approach to discretize the collision term, thereby loosening the constraint of the time step size.The second scheme is a multiscale fully implicit DVM (DVM-II).Differing from DVM-I, DVM-II employs a fully implicit scheme to discretize the collision term and introduces the corresponding macroscopic governing equation to forecast the equilibrium state at the new time level.The third scheme is a multiscale fully implicit DVM with inner iteration (DVM-III).In this scheme, an inner iteration is incorporated to solve the macroscopic governing equation, aiming to enhance the accuracy of the predicted equilibrium state.To assess the performance of these multiscale schemes, we conduct an asymptotic analysis and numerical experiments for flows in the near-continuous/continuous flow regimes.By doing so, we aim to uncover the key factors that facilitate convergence acceleration in DVM calculations.

Discrete Velocity Boltzmann-BGK Equation
The original Boltzmann equation, an integral-differential equation with a complex collision term, poses significant challenges in practical engineering applications.In this study, our focus is on the Boltzmann equation integrated with the BGK model [3].This amalgamation gives rise to the Boltzmann-BGK equation, which can be expressed as where f represents the distribution function defined in the physical space x, the molecular velocity space ξ, and the time t.Ω denotes the collision operator and τ is the collision time.
The equilibrium state, denoted by g, is defined as: Here, ρ represents the density, T signifies the temperature, c = ξ − u stands for the molecular thermal velocity vector, u denotes the mean flow velocity vector, c = |c| symbolizes the magnitude of c, and R g represents the gas constant.To solve Equation (1) using the DVM, it is necessary to truncate and discretize the molecular velocity space into a series of discrete velocity points.This discretization process results in where, N V represents the total number of discrete velocity points and the subscript α denotes the index in the discrete velocity space.
With Equation (1), we can derive the corresponding macroscopic governing equation of conservation laws through moment integration.By multiplying Equation (1) by the microscopic conservative moment Ψ = 1, ξ, ξ 2 /2 T and subsequently integrating the resulting equation in the molecular velocity space, we can obtain: where, W and F represent the conservative flow variables and fluxes, respectively.The notation f α = ∑ N V f α defines the numerical quadrature of distribution functions across the entire discrete molecular velocity space.As reported in previous studies [8,39], the conservation with respect to the collision integral plays a crucial role in ensuring the stability and accuracy of the DVM.However, the primary focus of this paper is to investigate the impact of different discretization strategies for the collision term on the computational efficiency of the DVM.To address the potential adverse effects caused by numerical quadrature errors, we will employ a relatively large number of points in the discretization of the molecular velocity space during the numerical tests.

Three Versions of DVM
To investigate the key factors that contribute to convergence acceleration in DVM calculations, this section introduces and compares three versions of DVM.In these schemes, similar to the DUGKS [27][28][29], the local discrete characteristic solution of the Boltzmann-BGK equation is employed to calculate the numerical flux, addressing the limitation associated with mesh size.This local solution is obtained by integrating Equation (1) from t n = 0 to t n + ∆t p /2 along the characteristic line and approximating the collision term using the trapezoidal rule.By doing so, we can obtain where x ij denotes the mid-point of the cell interface.h = ∆t p /2 signifies the half-time step size and ∆t p stands for the physical time step used only for flux reconstruction to avoid extrapolation, which is defined by Here, ξ x,max and ξ y,max represent the maximum molecular velocities in the xand ydirections, respectively.V signifies the volume of the control cell.∆S x and ∆S y are the projected areas of the control volume in the yand x-directions, respectively.σ p is the Courant-Friedrichs-Lewy (CFL) number, which is set as σ p = 0.95.
To fully determine the distribution function at the cell interface f α x ij , h , it is essential to pre-calculate the discrete distribution function at the surrounding point of the cell interface f a x ij − ξ α h, 0 , the equilibrium state at the surrounding point of the cell interface g α x ij − ξ α h, 0 , and the equilibrium state at the mid-point of the cell interface g α x ij , h .For the calculation of f a x ij − ξ α h, 0 and g α x ij − ξ α h, 0 , a second-order interpolation scheme employing van Leer's slope limiter is utilized.Taking φ to represent either f a or g α , we can obtain: where φ L and φ R represent the interfacial states of φ at the left and right sides, respectively.∇φ(x i , 0) and ∇φ x j , 0 denote the derivatives of φ at the left and right cells.L(φ, x i ) and L φ, x j are the corresponding slope limiter functions.
For the calculation of g α x ij , h , the macroscopic flow variables at the cell interface are required.This can be achieved by computing the conservative moments on both sides of Equation (7), which yields: The above equation can be rearranged as follows Since f a x ij − ξ α h, 0 and g a x ij − ξ α h, 0 have been obtained previously, W x ij , h as well as g α x ij , h can be computed explicitly.
Once the local discrete characteristic solution of the Boltzmann-BGK equation at the cell interface is determined, calculating the numerical fluxes for both the Boltzmann-BGK equation and the corresponding macroscopic governing equation becomes straightforward.This facilitates the evolution of both the discrete distribution functions and the conservative variables.In the following subsection, three different strategies are employed to advance the evolution of these quantities.

Semi-Implicit DVM (DVM-I)
The first multiscale approach is the semi-implicit DVM, similar to the work of Yang and Huang [18], where the equilibrium state in the collision term is calculated using the flow variables at the current time level, while other parts of the Boltzmann-BGK equation are discretized implicitly.This leads to the following discretized equation: where ∆t represents the time step size for the evolution of the distribution function, determined by the CFL condition.But different from Equation (8), the CFL number for the calculation of ∆t can be chosen larger than one to attain faster convergence.Integrating Equation ( 12) over a control volume V i , we can obtain a finite volume discretization of the Boltzmann-BGK equation as follows: where N(i) is the set of neighboring cells of the cell i, S ij stands for the area of the interface shared by the cells i and j, and n ij signifies the unit normal vector of the shared interface pointing from cell i to cell j.
By introducing the incremental ∆ f n i,α = f n+1 i,α − f n i,α into Equation ( 13), we can obtain: with the right-hand side of: Since the distribution function at the cell interface has been determined by Equation ( 7), the right-hand side of Equation ( 14) can be calculated explicitly.For the reconstruction of ∆ f n ij,α , the first-order upwind scheme is employed, i.e., where 16) into Equation ( 14) and applying the LU-SGS method to solve the resultant equation [40], the incremental ∆ f n i,α can be obtained by the following forward and backward sweepings: with: where L(i) is the subset of N(i) with elements less than i, while U(i) is the subset of N(i) with elements larger than i. ∆ f n i,α denotes the intermediate result of the forward sweep.Upon obtaining the incremental ∆ f n i,α , the discretized distribution function can be updated by f n+1 i,α = f n i,α + ∆ f n i,α and the macroscopic flow variables can be obtained by Equations ( 5) and (6).

Fully Implicit DVM (DVM-II)
Unlike the semi-implicit DVM, the fully implicit DVM calculates the equilibrium state within the collision term using the predicted solution from the corresponding macroscopic governing equation at the new time level.This strategy has found widespread use in the development of the fully implicit scheme, culminating in the following discretized Boltzmann-BGK equation [32][33][34] Following the solution process of Equation (12), Equation ( 19) can be expressed as: Entropy 2023, 25, 1609 7 of 28 with the right-hand side given by: Equation ( 20) can be easily solved using the LU-SGS method.However, before we proceed to solve this equation, it is necessary to determine the predicted equilibrium state ĝn+1 i,α .As we know, the macroscopic governing equation can be derived from the Boltzmann-BGK equation through moment integration.By multiplying Equation ( 19) by the microscopic conservative moment Ψ = 1, ξ, ξ 2 /2 T and subsequently integrating the resulting equation across the molecular velocity space, we can derive: By introducing the incremental ∆ Ŵn = Ŵn+1 − W n and integrating Equation (22) over a control volume V i , we can obtain: To solve the aforementioned equation, it is necessary to linearize the residual term Rn+1 i into the following form: Since the distribution function at the cell interface has been determined by Equation ( 7), the first term on the right-hand side of Equation ( 24) can be calculated explicitly.
For the calculation of the last term of Equation ( 24), the Euler equations-based flux splitting method is adopted: where, Here, F c = (ρu, ρuu + pI, (ρE + p)u) T denotes the convective flux of the macroscopic governing equation, c s represents the sound speed, γ signifies the specific heat ratio, Pr stands for the Prandtl number, and x i and x j represent the centroids of cell i and cell j, respectively.By substituting Equations ( 24) and (25) into Equation (23) and applying the LU-SGS method to solve the resultant equation [40], we can obtain: where

Fully Implicit DVM with Inner Iteration (DVM-III)
In DVM-II, it is observed that both the Boltzmann-BGK equation and the corresponding macroscopic governing equation iterate only once at each time step.In the nearcontinuum and continuum flow regimes, the macroscopic governing equation tends to approach the Navier-Stokes equation.However, the flux Jacobian of the macroscopic governing equation is approximated using the Euler equations-based flux splitting method, which might lead to inaccuracies in predicting the equilibrium state.To enhance the accuracy of the predicted equilibrium state in the near-continuum and continuum flow regimes, similar to prior studies [35][36][37][38], an inner iteration is introduced to solve the macroscopic governing equation.To achieve this goal, a time derivative term with respect to the pseudo-time variable η is added on the left-hand side of Equation ( 23) as follows: In DVM-III, the flux Jacobian of the macroscopic governing equation ∆ Rn i is calculated by taking the difference between the fluxes of the Navier-Stokes equation at the current time level and the new time level.This can be expressed as: where F n ij,NS and Fn+1 ij,NS represent the numerical fluxes of the Navier-Stokes equation at the current time level and the new time level, respectively.This equation signifies that the accurate difference in the numerical fluxes of the macroscopic governing equation, computed using the distribution function, is estimated by the disparity in numerical fluxes of the Navier-Stokes equation.This approximation is considered more reasonable than the Euler equations-based flux splitting method, particularly in the context of near-continuum and continuum flow regimes [36][37][38].
By discretizing the time derivative of Equation ( 28) with the backward Euler scheme and substituting Equation ( 29) into the resultant equation, we can obtain: where ∆ Fn+1,m c,j . Substituting Equation (32) into Equation (30), we can derive: with: Equation ( 33) can be easily resolved through the LU-SGS method [40].At the beginning of the inner iteration, the initial values of the conservative variables and fluxes are set as Ŵn+1,m=1 , where M signifies the iteration number of the inner iteration and it is chosen as 50 in this work.In addition, for simplicity, the time step size in the pseudo-time level is chosen as ∆η n+1,m i = ∆t n i .

Comparison of Three Schemes
The three schemes mentioned above all entail calculating the numerical flux using the local discrete characteristic solution of the Boltzmann-BGK equation.Additionally, they employ either a semi-implicit or fully implicit scheme to discretize the collision term.This allows for a mesh size unrestricted by the mean free path of molecules and a time step size not limited by the collision time.The primary distinction among these schemes lies in how they incorporate the solution of the macroscopic governing equation.In DVM-II, the solution of the macroscopic governing equation is introduced, while in DVM-III, an inner iteration is used to solve the macroscopic governing equation.However, considering that the computational cost of solving the macroscopic governing equation is significantly lower than that of the Boltzmann-BGK equation, the increased computational cost of DVM-II and DVM-III in each time step for evolving the discrete distribution function is negligible.In this subsection, we will qualitatively compare the three schemes in terms of the collisionless limit and continuum limit.For convenient comparison, let us rewrite the evolution equations of the three schemes as follows: In DVM-I, only the evolution equation of the distribution function is involved.This equation can be written as: In DVM-II, both the evolution equations of the distribution function and the macroscopic conservative variables are involved.These equations can be expressed as: Entropy 2023, 25, 1609 10 of 28 In DVM-III, the evolution equations of the distribution function and the macroscopic conservative variables maintain the same forms as Equations ( 35) and (36), with the distinction that the flux Jacobian is calculated by the disparity of fluxes of Navier-Stokes equation: Note that, in DVM-II and DVM-III, the predicted equilibrium state ĝn+1 α is calculated from the solution of the macroscopic governing equation Ŵn+1 .
(1) Collisionless limit In the collisionless limit, where the collision time is significantly larger than the time step size, τ ∆t, the evolution equation of the discrete distribution function for all three schemes can be simplified as follows: This simplification disregards any equilibrium state and can be considered as a solution to the following collisionless Boltzmann equation.
For DVM-II and DVM-III, even though the evolution of macroscopic conservative variables is involved, the predicted equilibrium state does not directly impact the evolution of the discrete distribution function.Hence, all three schemes can provide reasonable results in the collisionless limit.
(2) Continuum limit In the continuum limit, where the collision time is significantly less than the time step size, τ ∆t, the evolution equations of the discrete distribution function for DVM-I, DVM-II, and DVM-III can be respectively simplified as follows: For DVM-I: For DVM-II and -III: It can be observed that the effective time step size for the evolution of the discrete distribution function is the collision time.However, the initial values (the predicted equilibrium state) differ among the three schemes.A more accurate initial value leads to a faster convergence rate.
In order to achieve a more precise predicted equilibrium state, the evolution equation of the predicted macroscopic conservative variables (as shown in Equation ( 36)) is introduced in DVM-II and DVM-III.It is evident that the effective time step size of Equation ( 36) is much larger than the collision time in this scenario.This implies that the macroscopic governing equation takes precedence in the continuum limit.Consequently, the flow field can be accurately described by the Navier-Stokes equation in this regime.Naturally, DVM-III yields a more accurate predicted equilibrium state compared to DVM-II since the flux Jacobian is calculated using the Navier-Stokes equation in this method.

Numerical Examples
In this section, we perform a numerical investigation to evaluate the performance of the three schemes and uncover the key to accelerating convergence in DVM calculations.The evaluation encompasses six diverse test examples with varying Knudsen/Reynolds numbers, effectively covering a broad spectrum of flow regimes.These examples include Couette flow, heat transfer between parallel plates, force-driven Poiseuille flow, lid-driven cavity flow, flow around a NACA0012 airfoil, and flow in a planar microchannel.All computations are carried out on a PC equipped with an Intel(R) Xeon(R) Gold 6226R CPU operating at 2.9 GHz.No parallelization techniques are employed during these computations.

Case 1: Couette Flow
The first test case involves the Couette flow, consisting of two vertically positioned plates situated at x L = 0 and x R = 1.Both plates have the same temperature of T 0 = 1 and different constant vertical velocities: v L = −0.25 for the left plate and v R = 0.25 for the right plate.The Maxwellian diffuse boundary condition is employed for the left and right boundaries, while the periodic boundary condition is adopted for the upper and lower boundaries.The Knudsen number for this case is defined as: where ρ 0 = 1 is the reference density, L 0 = x R − x L denotes the reference length, µ 0 represents the reference viscosity, and R g = 0.5 signifies the specific gas constant.The dynamic viscosity is determined using the following intermolecular interaction model: with the viscosity index w = 0.81.In the simulation, we consider five different Knudsen numbers: Kn = 0.001, 0.01, 0.1, 1, and 10.The physical space is discretized non-uniformly into 100 cells in the x-direction and 5 cells in the y-direction.For the discretization of the molecular velocity space, we use the Gauss-Hermite quadrature with 28 × 28 mesh points for the cases of Kn = 0.001, 0.01, and 0.1.In the cases of Kn = 1 and 10, we utilize the Newton-Cotes quadrature with 101 × 101 mesh points uniformly distributed in the range of The converged velocity and temperature distributions along the x-direction are shown in Figure 1, demonstrating good agreement among DVM-I, DVM-II, and DVM-III.Figure 2 compares the convergence history of the three schemes.It is evident that, in the rarefied flow regime, all three schemes converge similarly.However, in the near-continuum and continuum flow regimes, the DVM-III achieves the fastest convergence, followed by the DVM-II, while the DVM-I exhibits slower convergence.This indicates that the fully implicit discretization of the collision term is pivotal for accelerating convergence in DVM calculations, and a more accurate prediction of the equilibrium state leads to a faster convergence rate.The quantitative comparison of the computational cost of the three schemes is tabulated in Table 1.Since the computational cost of solving the macroscopic governing equation is significantly lower than that of the Boltzmann-BGK equation, DVM-II achieves approximately one order of magnitude acceleration compared to DVM-I.Additionally, DVM-III achieves an additional order of magnitude acceleration compared to DVM-II.
calculations, and a more accurate prediction of the equilibrium state leads to a faster convergence rate.The quantitative comparison of the computational cost of the three schemes is tabulated in Table 1.Since the computational cost of solving the macroscopic governing equation is significantly lower than that of the Boltzmann-BGK equation, DVM-II achieves approximately one order of magnitude acceleration compared to DVM-I.Additionally, DVM-III achieves an additional order of magnitude acceleration compared to DVM-II.

Case 2: Heat Transfer between Two Parallel Plates
The second test case involves heat transfer between two stationary parallel plates.The geometry configuration and mesh in the physical space for this test case are the same as those in the first test case.The left plate is maintained at a fixed temperature of 0.75 L T = , while the right plate is kept at a fixed temperature of 1.25 R T = .The Knudsen number is determined by Equation ( 43) with the reference density of 0 1 ρ = and refer- ence temperature of 0 1 T = , and the dynamic viscosity is calculated using Equation (44) with the viscosity index 0.81 w = .In the simulation, we consider five different Knudsen numbers: Kn = 0.001, 0.01, 0.1, 1, and 10, which cover all the flow regimes.For each Knud-

Case 2: Heat Transfer between Two Parallel Plates
The second test case involves heat transfer between two stationary parallel plates.The geometry configuration and mesh in the physical space for this test case are the same as those in the first test case.The left plate is maintained at a fixed temperature of T L = 0.75, while the right plate is kept at a fixed temperature of T R = 1.25.The Knudsen number is determined by Equation ( 43) with the reference density of ρ 0 = 1 and reference temperature of T 0 = 1, and the dynamic viscosity is calculated using Equation ( 44) with the viscosity index w = 0.81.In the simulation, we consider five different Knudsen numbers: Kn = 0.001, 0.01, 0.1, 1, and 10, which cover all the flow regimes.For each Knudsen number, the molecular velocity space discretization follows the same strategy as that in the first test case.
Figure 3 illustrates a comparison of the converged density and temperature distributions along the x-direction obtained using DVM-I, DVM-II, and DVM-III.The results from all three schemes show good agreement with each other, indicating that these multiscale approaches can provide reasonable results across different flow regimes.To further analyze the convergence behavior, Figure 4 presents a comparison of the convergence history.Similar to the first test case, DVM-III demonstrates the fastest convergence among the three schemes.This can be attributed to its ability to provide more accurate predictions of the equilibrium state in the discretization of the collision term.As a result, DVM-III achieves approximately 1-2 orders of magnitude acceleration compared to DVM-I and about one order of magnitude acceleration compared to DVM-II, as shown in Table 2.

Case 3: Force-Driven Poiseuille Flow
This test case involves the flow of a fluid between two infinite parallel plates separated by a distance of L 0 = 1 in the y-direction.The plates maintain a temperature of T 0 = 1 at the reference temperature.In the force-driven Poiseuille flow, an external force is imposed in the x-direction.Consequently, the Boltzmann-BGK equation is modified as follows: where F α,x denotes the force term.In the context of the Maxwellian distribution function, F α,x can be expressed as: Here, ξ α,x stands for the x-component of the discrete molecular velocity, u x denotes the x-component of the mean flow velocity, and G represents the magnitude of the external acceleration.With the external force, the corresponding modified macroscopic governing equation is as follows: where S = (0, G, 0, uG) T is the source term.In this test example, the Knudsen number is calculated by Equation ( 43) with the reference density of ρ 0 = 1 and the dynamic viscosity is determined using Equation (44) with the viscosity index w = 0.5.The apparent gas permeability is introduced to quantify the simulation results, which is defined as: In the simulation, we consider Knudsen numbers ranging from Kn = 10 −4 to 10.The magnitude of the external acceleration is set as follows: G = 10 −5 for the cases of 10 −4 ≤ Kn < 10 −3 , G = 10 −4 for the cases of 10 −3 ≤ Kn < 10 −2 , G = 10 −3 for the cases 10 −2 ≤ Kn < 10 −1 , and G = 10 −2 for the cases 10 −1 ≤ Kn.The physical space is discretized uniformly into 40 cells in the y-direction and 5 cells in the x-direction.For the discretization of the molecular velocity space, we utilize the Gauss-Hermite quadrature with 8 × 8 mesh points and 28 × 28 mesh points for the cases of Kn < 10 −2 and 10 −2 ≤ Kn < 1, respectively.For the cases of 1 ≤ Kn ≤ 10, we employ the Newton-Cotes quadrature with 101 × 101 mesh points uniformly distributed in the range of [−4, 4] Figure 5 presents a comparison of the apparent gas permeability for force-driven Poiseuille flow with different Knudsen numbers using DVM-I, DVM-II, DVM-III, and DUGKS [41] as a reference.The results demonstrate that all three schemes yield consistent results with the reference data, indicating their accuracy as multiscale approaches.To evaluate the computational efficiency of the schemes, we analyze the convergence history of the apparent gas permeability for Kn = 0.0001, 0.001, 0.01, and 0.1, as shown in Figure 6.
It is evident that DVM-III exhibits faster convergence compared to DVM-I and DVM-II.Specifically, for the case of Kn = 0.0001, DVM-III achieves approximately five and two orders of magnitude faster convergence compared to DVM-I and DVM-II, respectively.Even in the slip flow regime (Kn = 0.01), DVM-III still demonstrates an approximate two orders of magnitude faster convergence compared to DVM-II.These findings confirm that a more accurate prediction of the equilibrium state in the collision term discretization leads to a faster convergence rate in DVM calculations.

Case 4: Lid-Driven Cavity Flow
In this subsection, we will simulate the two-dimensional lid-driven cavity flow, which introduces a non-trivial dimension compared to the previous test examples.The square cavity with an edge length of L 0 = 1 remains stationary, except for the top wall that moves with a velocity of u W = 0.15.All walls are maintained at a fixed temperature equal to the reference temperature of T 0 = 1.The simulation considers three different Knudsen numbers: Kn = 0.075, 1 and 10, as well as two different Reynolds numbers: Re = 100 and 1000.For the cases with Kn = 0.075, 1, and 10, the dynamic viscosity is computed using Equation ( 44) with the viscosity index w = 0.81 and the reference dynamic viscosity µ 0 is determined by the Knudsen number: For the cases with Re = 100 and 1000, the dynamic viscosity µ is directly calculated by the Reynolds number: Additionally, the computational domain is uniformly discretized into 50 × 50 cells for the cases with Kn = 0.075, 1 and 10, while it is divided into 150 × 150 cells for the cases with Re = 100 and 1000.Regarding the discretization of the molecular velocity space, the Newton-Cotes quadrature with 101 × 101 points uniformly distributed in the range of [−4, 4] × [−4, 4] is employed for the test cases with Kn = 1 and 10.The Gauss-Hermite quadrature with 28 × 28 points is utilized for the test case with Kn = 0.075, and the Gauss-Hermite quadrature with 8 × 8 points is used for the test cases with Re = 100 and 1000.
The comparisons of density, temperature, the x-component of heat flux, and the ycomponent of heat flux contours for lid-driven cavity flow with Knudsen numbers of 0.075, 1, and 10 are displayed in Figures 7-9, respectively.The results obtained from all three schemes (DVM-I, DVM-II, and DVM-III) exhibit good agreement, indicating their capability to accurately capture the flow behavior across different regimes.Figure 10 provides a quantitative comparison of the velocity profiles along the vertical and horizontal central lines of the cavity for the three schemes, along with the results of UGKS [42] and the numerical results of Ghia et al. [43] as reference data.Once again, the results obtained from all three schemes closely align with the reference data, further supporting the applicability of these multiscale approaches in the different flow regimes.To assess the convergence performance, Figure 11 compares the convergence history for lid-driven cavity flow with different Knudsen/Reynolds numbers.Additionally, Table 3 quantitatively compares the computational cost among the three schemes.Clearly, DVM-III exhibits the fastest convergence rate among the three schemes, achieving approximately 1-2 orders of magnitude faster convergence compared to DVM-I.This observation confirms that the accuracy in predicting the equilibrium state for the collision term's discretization plays a critical role in accelerating convergence in DVM calculations.

Case 5: Flow around a NACA0012 Airfoil
This test case involves the flow around a NACA0012 airfoil with a Mach number of Ma = 0.2 and an angle of attack of AoA = 10 degrees.The temperature of the airfoil is maintained at the reference temperature of 0 1 T = .The Knudsen number is defined by Equa-

Case 5: Flow around a NACA0012 Airfoil
This test case involves the flow around a NACA0012 airfoil with a Mach number of Ma = 0.2 and an angle of attack of AoA = 10 degrees.The temperature of the airfoil is maintained at the reference temperature of T 0 = 1.The Knudsen number is defined by Equation ( 43) using the reference density of ρ 0 = 1 and the normalized airfoil chord length of L 0 = 1, and the dynamic viscosity is determined using Equation (44) with the viscosity index w = 0.81.Consequently, the relationship between the Knudsen number, Mach number, and Reynolds number can be expressed as: In the simulation, we consider three different Reynolds numbers: Re = 5, 50, and 500.The computational domain extends up to a distance of 25L 0 from the center of the airfoil and is discretized using an unstructured mesh comprising 16,042 triangular cells.On the airfoil surface, we place 200 discrete points.For the discretization of the molecular velocity space, we employ the Gauss-Hermite quadrature with a mesh of 28 × 28 mesh points.
Figures 12-14 present the density, u-velocity, v-velocity, and pressure contours for flows around a NACA0012 airfoil at Reynolds numbers of 5, 50, and 500, respectively.The results obtained by the DVM-I, DVM-II, and DVM-III are found to be consistent with each other.Figures 15-17 compare the pressure coefficient and skin friction coefficient distributions along the airfoil surface for the cases of Re = 5, 50, and 500.These figures also include numerical results obtained using the Navier-Stokes equation with a conventional CFD scheme [44].It is observed that the results obtained by DVM-I, DVM-II, and DVM-III align well with each other and show good agreement with the results from the Navier-Stokes equation for the case of Re = 500.However, for the cases of Re = 5 and 50, there are deviations between the results obtained by the three schemes and those from the Navier-Stokes equation due to the rarefaction effect.In Figure 18, we compare the convergence history of the three schemes at different Reynolds numbers.It is evident that DVM-III exhibits faster convergence compared to the other two schemes.Specifically, for the case of Re = 500, there is an acceleration of approximately one order of magnitude compared to DVM-I.

Case 6: Flow in a Planar Microchannel
The last test example focuses on the flow in a planar microchannel, which has been previously investigated by Titarev [45,46].The configuration of this test case is illustrated in Figure 19 and consists of a left reservoir with pressure p 1 and temperature T 1 , a right reservoir with pressure p 2 and temperature T 2 , and a channel with the aspect ratio of l/a = 10, where l represents half of the length and a represents half of the width.In our simulation, p 1 /p 2 = 1.1 and T 1 /T 2 = 1 are considered.To compare our results with those obtained by Titarev [45,46], we calculate the collision time using τ = 1/δ, where δ is a rarefaction parameter that is inversely proportional to the Knudsen number.In the simulation, we examine six different values of δ: 0.01, 0.1, 1, 10, 100, and 1000.Note that δ = 1000 corresponds to the continuum flow regime, which poses a challenge for DVM-I to efficiently achieve convergence and was not tested in the work of Titarev [45,46]. .Consequently, the relationship between the Knudsen number, Mach number, and Reynolds number can be expressed as: In the simulation, we consider three different Reynolds numbers: Re = 5, 50, and 500.
The computational domain extends up to a distance of 0 25L from the center of the airfoil and is discretized using an unstructured mesh comprising 16,042 triangular cells.On the airfoil surface, we place 200 discrete points.For the discretization of the molecular velocity space, we employ the Gauss-Hermite quadrature with a mesh of 28 × 28 mesh points.To quantitatively compare the simulation results with the reference data provided by Titarev [45,46], the non-dimensional mass flow rate is introduced with: where p 0 = (p 1 + p 2 )/2 and T 0 = (T 1 + T 2 )/2 represent the reference pressure and ref- erence temperature, respectively.In the simulation, the physical space is discretized by 5501 quadrilateral cells, and the molecular velocity space is discretized by the Gauss-Hermite quadrature with a mesh of 28 × 28 mesh points for δ ≤ 100 and 8 × 8 mesh points for δ ≤ 1000.

Case 6: Flow in a Planar Microchannel
The last test example focuses on the flow in a planar microchannel, which has been previously investigated by Titarev [45,46].The configuration of this test case is illustrated in Figure 19 and consists of a left reservoir with pressure 1 p and temperature 1 T , a right reservoir with pressure 2 p and temperature 2 T , and a channel with the aspect ratio of corresponds to the continuum flow regime, which poses a challenge for DVM-I to efficiently achieve convergence and was not tested in the work of Titarev [45,46].
To quantitatively compare the simulation results with the reference data provided by Titarev [45,46], the non-dimensional mass flow rate is introduced to DVM-I.Furthermore, Table 4 presents the non-dimensional mass flow rates obtained by the three schemes, alongside the reference data provided by Titarev [46].Basically, the simulation results demonstrate good agreement with the reference data, validating the accuracy of all three schemes.4 presents the non-dimensional mass flow rates obtained by the three schemes, alongside the reference data provided by Titarev [46].Basically, the simulation results demonstrate good agreement with the reference data, validating the accuracy of all three schemes.

Conclusions
In this study, three versions of multiscale DVM are compared to investigate the key factor in accelerating convergence in DVM calculations.To ensure accuracy across different flow regimes, particularly in the near-continuum and continuum flows, these approaches use the local discrete characteristic solution of the Boltzmann-BGK equation to calculate the numerical flux at the cell interface, similar to the DUGKS.The first version, DVM-I, employs a semi-implicit scheme to discretize the collision term.It approximates the equilibrium state using its current time step value.On the other hand, DVM-II and DVM-III utilize a fully implicit scheme to discretize the Boltzmann-BGK equation, including the collision term.In these versions, the equilibrium state is predicted based on the solution of the corresponding macroscopic governing equation.Notably, DVM-III introduces an inner iteration of the macroscopic governing equation between adjacent DVM steps, leading to a more accurate prediction of the equilibrium state for the fully implicit discretization of the collision term.
To investigate the key factor in accelerating convergence in DVM calculations, simulations have been conducted on six benchmark cases, including Couette flow, heat transfer between two parallel plates, force-driven Poiseuille flow, lid-driven cavity flow, flow around a NACA0012 airfoil, and flow in a planar microchannel.By considering the collisional effect in the calculation of the numerical flux at the cell interface, all three multiscale approaches (DVM-I, DVM-II, and DVM-III) yield reasonable results in good agreement.Concerning computational efficiency, DVM-III exhibits the highest efficiency in near-continuum and continuum flow regimes, followed by DVM-II, while DVM-I shows a lower efficiency.This indicates that the fully implicit discretization of the collision term plays a crucial role in accelerating convergence in DVM computations.Additionally, the more precise prediction of the equilibrium state resulted in a higher convergence rate.These findings provide valuable insights into the development of efficient and accurate multiscale approaches for simulating flows around irregular objects in near-space environments [47,48].However, it is crucial to acknowledge that for practical applications, parallel implementation becomes necessary due to the substantial computational cost associated with the DVM.Employing multicore calculations alongside the LU-SGS scheme can significantly impact convergence.Furthermore, implementing parallelization for DVM-II and DVM-III, where resolving the macroscopic governing equation is required, poses greater challenges compared to DVM-I.Developing parallel versions for these schemes will be addressed in future work.
NS represent the increments in the predicted conservative variable vector and the flux vector of the Navier-Stokes equation at the pseudo-time level m of the inner iteration.During the inner iteration process, as the solution converges, the values of ∆ Ŵn+1,m i and ∆ Fn+1,m ij,NS tend to approach zero.Thus, the Euler equations-based flux splitting method can be employed to calculate the flux Jacobian in the inner iteration, ij,NS .At the conclusion of the inner iteration, we can set Ŵn+1 i = Ŵn,m=M i

Figure 1 .
Figure 1.Velocity (a) and temperature (b) distributions along the x-direction for Couette flow with different Knudsen numbers.

Figure 1 . 29 Figure 2 .
Figure 1.Velocity (a) and temperature (b) distributions along the x-direction for Couette flow with different Knudsen numbers.Entropy 2023, 25, x FOR PEER REVIEW 13 of 29

Figure 2 .
Figure 2. Convergence history of DVM-I, DVM-II, and DVM-III for Couette flow with different Knudsen numbers.

Figure 3 .
Figure 3. Density (a) and temperature (b) distributions along the x-direction for heat transfer between two parallel plates with different Knudsen numbers.

Figure 4 .
Figure 4. Convergence history of DVM-I, DVM-II, and DVM-III for heat transfer between two parallel plates with different Knudsen numbers.

Figure 3 .
Figure 3. Density (a) and temperature (b) distributions along the x-direction for heat transfer between two parallel plates with different Knudsen numbers.

Figure 3 .
Figure 3. Density (a) and temperature (b) distributions along the x-direction for heat transfer between two parallel plates with different Knudsen numbers.

Figure 4 .
Figure 4. Convergence history of DVM-I, DVM-II, and DVM-III for heat transfer between two parallel plates with different Knudsen numbers.

Figure 4 .
Figure 4. Convergence history of DVM-I, DVM-II, and DVM-III for heat transfer between two parallel plates with different Knudsen numbers.

Figure 6 .
Figure 6.Convergence history of the apparent gas permeability for force-driven Poiseuille flow with different Knudsen numbers (a-d).

Figure 6 .
Figure 6.Convergence history of the apparent gas permeability for force-driven Poiseuille flow with different Knudsen numbers (a-d).

Figure 6 .
Figure 6.Convergence history of the apparent gas permeability for force-driven Poiseuille flow with different Knudsen numbers (a-d).

Figure 7 .
Figure 7.Comparison of (a) density, (b) temperature, (c) x-component of heat flux, and (d) y-component of heat flux contours for lid-driven cavity flow with Kn = 0.075 (DVM-I: white dash line; DVM-II: red dash-dot line; and DVM-III: colored background).

Figure 8 .
Figure 8.Comparison of (a) density, (b) temperature, (c) x-component of heat flux, and (d) y-component of heat flux contours for lid-driven cavity flow with Kn = 1 (DVM-I: white dash line; DVM-II: red dash-dot line; and DVM-III: colored background).

Figure 9 .
Figure 9.Comparison of (a) density, (b) temperature, (c) x-component of heat flux, and (d) y-component of heat flux contours for lid-driven cavity flow with Kn = 10 (DVM-I: white dash line; DVM-II: red dash-dot line; and DVM-III: colored background).

Figure 9 .
Figure 9.Comparison of (a) density, (b) temperature, (c) x-component of heat flux, and (d) y-component of heat flux contours for lid-driven cavity flow with Kn = 10 (DVM-I: white dash line; DVM-II: red dash-dot line; and DVM-III: colored background).

Figure 10 .
Figure 10.Velocity profiles along the vertical and horizontal central lines of the cavity with different Knudsen (a)/Reynolds (b) numbers [43].

4. 5 .
Case 5: Flow around a NACA0012 Airfoil This test case involves the flow around a NACA0012 airfoil with a Mach number of Ma = 0.2 and an angle of attack of AoA = 10 degrees.The temperature of the airfoil is main-

Figure 10 . 29 Figure 10 .
Figure 10.Velocity profiles along the vertical and horizontal central lines of the cavity with different Knudsen (a)/Reynolds (b) numbers [43].

Figure 11 .
Figure 11.Convergence history for lid-driven cavity flow with different Knudsen (a)/Reynolds (b) numbers.

Figure 11 .
Figure 11.Convergence history for lid-driven cavity flow with different Knudsen (a)/Reynolds (b) numbers.

Entropy 2023 , 1 L
25, x FOR PEER REVIEW 21 of 29 0 = , and the dynamic viscosity is determined using Equation (44) with the viscosity index 0.81 w =

Figures 12 -
14 present the density, u-velocity, v-velocity, and pressure contours for flows around a NACA0012 airfoil at Reynolds numbers of 5, 50, and 500, respectively.The results obtained by the DVM-I, DVM-II, and DVM-III are found to be consistent with each other.Figures 15-17 compare the pressure coefficient and skin friction coefficient distributions along the airfoil surface for the cases of Re = 5, 50, and 500.These figures also include numerical results obtained using the Navier-Stokes equation with a conventional CFD scheme[44].It is observed that the results obtained by DVM-I, DVM-II, and DVM-III align well with each other and show good agreement with the results from the Navier-Stokes equation for the case of Re = 500.However, for the cases of Re = 5 and 50, there are deviations between the results obtained by the three schemes and those from the Navier-Stokes equation due to the rarefaction effect.In Figure18, we compare the convergence history of the three schemes at different Reynolds numbers.It is evident that DVM-III exhibits faster convergence compared to the other two schemes.Specifically, for the case of Re = 500, there is an acceleration of approximately one order of magnitude compared to DVM-I.

Figure 15 .
Figure 15.Comparison of pressure coefficient (a) and skin friction coefficient (b) distributions along the airfoil surface for flow around a NACA0012 airfoil with Re = 5.

Figure 15 .
Figure 15.Comparison of pressure coefficient (a) and skin friction coefficient (b) distributions along the airfoil surface for flow around a NACA0012 airfoil with Re = 5.

Figure 15 .
Figure 15.Comparison of pressure coefficient (a) and skin friction coefficient (b) distributions along the airfoil surface for flow around a NACA0012 airfoil with Re = 5.

Figure 16 .
Figure 16.Comparison of pressure coefficient (a) and skin friction coefficient (b) distributions along the airfoil surface for flow around a NACA0012 airfoil with Re = 50.

Figure 17 .
Figure 17.Comparison of pressure coefficient (a) and skin friction coefficient (b) distributions along the airfoil surface for flow around a NACA0012 airfoil with Re = 500.

Figure 16 . 29 Figure 16 .
Figure 16.Comparison of pressure coefficient (a) and skin friction coefficient (b) distributions along the airfoil surface for flow around a NACA0012 airfoil with Re = 50.

Figure 17 .
Figure 17.Comparison of pressure coefficient (a) and skin friction coefficient (b) distributions along the airfoil surface for flow around a NACA0012 airfoil with Re = 500.

Figure 17 .
Figure 17.Comparison of pressure coefficient (a) and skin friction coefficient (b) distributions along the airfoil surface for flow around a NACA0012 airfoil with Re = 500.

Figure 18 .
Figure 18.Convergence history for flow around a NACA0012 airfoil with different Reynolds numbers (a-c).
10 l a = , where l represents half of the length and a represents half of the width.In our simulation, 1 2 1.1 p p = and 1 2 1 T T = are considered.To compare our results with those obtained by Titarev [45,46], we calculate the collision time using 1 τ δ = , where δ is a rarefaction parameter that is inversely proportional to the Knudsen number.In the simulation, we examine six different values of δ : 0.01, 0.1, 1, 10, 100, and 1000.Note that 1000 δ =

Figure 18 .
Figure 18.Convergence history for flow around a NACA0012 airfoil with different Reynolds numbers (a-c).

Figure 19 .
Figure 19.Geometry and computational mesh in the physical space for flow in a planar microchannel.

Figure 19 .
Figure 19.Geometry and computational mesh in the physical space for flow in a planar microchannel.

Figure 20 illustrates
Figure 20 illustrates the convergence history for the flow in a planar microchannel with different rarefaction parameters.In this figure, we only display the convergence history obtained by DVM-I, DVM-II, and DVM-III for the cases of δ ≥ 1, since the convergence rates of the three schemes are identical for the cases of δ = 0.01 and 0.1.It is evident that DVM-III exhibits the fastest convergence, followed by DVM-II, while DVM-I converges at a slower rate.Particularly for the case of δ = 1000, both DVM-II and DVM-III achieve a significant acceleration of convergence by two orders of magnitude compared to DVM-I.Furthermore, Table4presents the non-dimensional mass flow rates obtained by the three schemes, alongside the reference data provided by Titarev[46].Basically, the simulation

Figure 19 .
Figure 19.Geometry and computational mesh in the physical space for flow in a planar microchannel.

Figure 20 .
Figure 20.Convergence history for flow in a planar microchannel with different rarefaction parameters (a-d).

Figure 20 .
Figure 20.Convergence history for flow in a planar microchannel with different rarefaction parameters (a-d).

Table 1 .
Computational cost (hours) of DVM-I, DVM-II, and DVM-III for Couette flow with different Knudsen numbers.

Table 1 .
Computational cost (hours) of DVM-I, DVM-II, and DVM-III for Couette flow with different Knudsen numbers.

Table 2 .
Computational cost (hours) of DVM-I, DVM-II, and DVM-III for heat transfer between two parallel plates with different Knudsen numbers.

Table 2 .
Computational cost (hours) of DVM-I, DVM-II, and DVM-III for heat transfer between two parallel plates with different Knudsen numbers.

Table 2 .
Computational cost (hours) of DVM-I, DVM-II, and DVM-III for heat transfer between two parallel plates with different Knudsen numbers.
Note: "Ratio 1" and "Ratio 2" have the same definitions as provided in Table1.

Table 3 .
Computational cost (hours) of the DVM-I, DVM-II, and DVM-III for lid-driven cavity flow with different Knudsen/Reynolds numbers.
Note: "Ratio 1" and "Ratio 2" have the same definitions as provided in Table1.

Table 3 .
Computational cost (hours) of the DVM-I, DVM-II, and DVM-III for lid-driven cavity flow with different Knudsen/Reynolds numbers.
Note: "Ratio 1" and "Ratio 2" have the same definitions as provided in Table1

Table 3 .
Computational cost (hours) of the DVM-I, DVM-II, and DVM-III for lid-driven cavity flow with different Knudsen/Reynolds numbers.
Note: "Ratio 1" and "Ratio 2" have the same definitions as provided in Table1.

Table 4 .
Mass flow rate for flow in a planar microchannel with different rarefaction parameters.