Comparative Study on Flux Solution Methods of Discrete Unified Gas Kinetic Scheme

TL;DR

This paper compares different flux calculation methods in a gas kinetic scheme, focusing on their accuracy and efficiency in various flow simulations.

Contribution

The paper introduces the Simpson method for flux calculation in DUGKS and compares its performance with other variants.

Findings

Simpson–DUGKS reduces computational time by about 59% compared to original DUGKS under certain conditions.

Optimized DUGKS shows the smallest L2 error for velocity in Taylor–Green vortex flow but the largest for density.

Optimized DUGKS results deviate more from reference solutions in lid-driven cavity flow with small grid numbers.

Abstract

In this work, the Simpson method is proposed to calculate the interface flux of a discrete unified gas kinetic scheme (DUGKS) according to the distribution function at the node and the midpoint of the interface, which is noted by Simpson–DUGKS. Moreover, the optimized DUGKS and Simpson–DUGKS considering the force term are derived. Then, the original DUGKS, optimized DUGKS, and Simpson–DUGKS are compared and analyzed in theory. Finally, the numerical tests are performed under different grid numbers (N). In the steady unidirectional flow (Couette flow and Poiseuille flow), the three methods are stable under different Courant–Friedrichs–Lewy (CFL) numbers, and the calculated L2 errors are the same. In the Taylor–Green vortex flow, the L2 error of the optimized DUGKS is the smallest with respect to the analytical solution of velocity, but the L2 error of the optimized DUGKS is the largest with respect to the analytical solution of density. In the lid-driven cavity flow, the results of the optimized DUGKS deviate more from the reference results in terms of accuracy, especially in the case of a small grid number. In terms of computational efficiency, it should be noted that the computational time of optimized DUGKS increases by about 40% compared with the original DUGKS when CFL = 0.1 and N = 16, and the calculation time of Simpson–DUGKS is reduced by about 59% compared with the original DUGKS when CFL = 0.95 and N = 16.

Figures

Full text

1. Introduction

Multiscale flows, characterized by significant spatial and temporal variations, are prevalent in natural phenomena and industrial processes. These flows pose substantial challenges to unified modeling and computational frameworks due to involving a large span of space and time [1]. In recent years, numerical methods for multiscale flow simulation have gained considerable attention across diverse fields, including micro–nano device design, aerospace engineering, and subsurface resource extraction [2,3,4,5,6,7,8,9]. For instance, during spacecraft atmospheric re-entry, the transition from continuum to rarefied flow regimes necessitates multiscale analysis [10,11,12]. While the Navier–Stokes equations, based on continuum assumptions, fail in rarefied regimes, the Boltzmann equation—derived from molecular dynamics principles—provides a more robust foundation. Numerical approaches to solving the Boltzmann equation include the following three categories: Particle-based methods decouple the transport and collision terms of the Boltzmann equation, simulating molecular streaming and collisions through computational particles. The direct simulation Monte Carlo (DSMC) method is a prominent example [13,14,15]. Discrete velocity methods (DVM) include the gas-kinetic unified algorithm (GKUA) [16], unified gas-kinetic scheme (UGKS) [17], and discrete unified gas kinetic scheme (DUGKS) [18]. By concurrently addressing transport and collision effects in flux calculations, these methods eliminate restrictions on the grid size and time steps imposed by molecular mean-free paths and collision times. This enables unified simulations across all flow regimes, overcoming DSMC’s limitations in continuum and near-continuum flows. In addition to stochastic particle methods and deterministic approaches, it would be valuable to include discussions on coupled methods. For example, under the UGKS framework, the unified gas-kinetic wave-particle (UGKWP) method demonstrates notable advantages in both computational efficiency and multiscale simulations across different flow regimes [19,20,21].

DUGKS distinguishes itself from the lattice Boltzmann method (LBM) and UGKS in several ways. While LBM primarily targets Navier–Stokes solutions for continuum flows and struggles in rarefied regimes, DUGKS achieves full flow-regime coverage across all Knudsen (Kn) numbers, capturing phenomena inaccessible to LBM [22,23]. As a finite-volume scheme, DUGKS naturally accommodates unstructured grids, enhancing geometric flexibility. Compared with UGKS, which relies on analytical time-integrated flux solutions, DUGKS employs a simplified numerical characteristic approach, reducing computational costs to approximately 70% of UGKS [24], while maintaining accuracy. DUGKS has been successfully applied to a wide range of multiscale flow problems, including turbulent flows [25,26,27,28], microscale flows (e.g., Knudsen layer effects, velocity slip) [29,30,31,32,33], compressible flows (shock waves, high-speed aerodynamics) [34,35,36,37], multiphase and particle-laden flows [38,39,40,41,42], complex geometries (3D unstructured meshes, porous media) [43,44,45,46,47,48,49], plasma dynamics and gas mixtures [45,46,50,51], and multiscale transport phenomena (phonon heat transfer [52,53,54], radiative transport [55,56]).

However, Wang et al. [57] think that the computational efficiency of DUGKS is not satisfactory for the following reasons: (1) The time step is determined by a small Courant–Friederichs–Lewy (CFL) number, which limits the overall computational efficiency. (2) The DUGKS framework is somewhat complicated due to the distribution function at the cell interface, which reduces the computational efficiency. Wang et al. [22] found that DUGKS is four times slower than LBM in the simulation of two-dimensional incompressible flow with the same grid resolution. Therefore, Wang et al. [57] proposed an optimized DUGKS to improve the efficiency of the original DUGKS. The evolution equation of the optimized DUGKS is approximately the same as that of the original DUGKS. Different from the original DUGKS, the optimized DUGKS uses a trapezoidal method to calculate the interface flux according to the distribution function at the nodes. However, there are higher-order methods to calculate interface flux than the trapezoidal method. Therefore, the Simpson method is proposed to calculate the interface flux according to the distribution function at the node and the midpoint of the interface, which is called Simpson–DUGKS in this paper. Moreover, the optimized DUGKS and Simpson–DUGKS considering force terms are derived. Then, the original DUGKS, optimized DUGKS, and Simpson–DUGKS are simply compared and analyzed in theory. Finally, the performance of original DUGKS, optimized DUGKS, and Simpson–DUGKS in computational efficiency, accuracy, and stability is tested by numerical examples, and a conclusion is drawn.

2. Numerical Method

2.1. Simplified Governing Equations of DUGKS with Force Term

Based on kinetic theory of gases, Guo et al. proposed the original DUGKS [18]. The basic variable is particle distribution function $[eqn]$ for molecules moving with velocity $[eqn]$ at position $[eqn]$ and time $[eqn]$ . Derived from the Boltzmann equation, the governing equation with force term DUGKS is as follows [23]:

$$[eqn]$$

Equation (1) describes the spatial and temporal evolution process of particle distribution function $[eqn]$ . The collision term Ω approximated by the Bhatnagar–Gross–Krook (BGK) model is as follows:

$$[eqn]$$

The relaxation time $[eqn]$ is computed from $[eqn]$ , where $[eqn]$ , R, and T represent the kinematic viscosity, the gas constant, and the temperature, respectively.

The force term S is approximately [23]:

$$[eqn]$$

where $[eqn]$ is the acceleration caused by the force term. In isothermal flow, the temperature is constant with $[eqn]$ .

At a low Mach number, Taylor expansion is carried out around zero particle velocity, and the Maxwell equilibrium distribution function f^eq^ is approximately [18]:

$$[eqn]$$

The symbol $[eqn]$ represents the lattice speed of sound, which is a constant depending on the lattice arrangement ( $[eqn]$ for the D2Q9 model). The symbol ρ represents the macroscopic density.

In the D2Q9 (two-dimensional and nine velocity) model [55], the velocity space set { $[eqn]$ } and the corresponding weight set { $[eqn]$ } are as follows:

$$[eqn]$$ $$[eqn]$$

The symbol c is the basic speed on the lattice (c = 1 usually).

In DUGKS with the D2Q9 model, the computational domain is divided into many cells by a Cartesian grid, and each cell has the size V_C_ = ΔxΔy, whose center is x_C_, as shown in Figure 1. As the finite volume scheme, the cell-average value of the distribution function and the force term are introduced.

$$[eqn]$$

The subscript C represents the location in the cell center, and the superscript n represents the time t^n^. At each cell V_C_, Equation (1) is integrated from time t^n^ = nΔt to time t^n^^+1^ = (n + 1)Δt. The midpoint method is used to calculate the time integrals of the convection term, and the trapezoidal method is used to calculate the time integrals of the collision term and force term. The simplified governing equation of DUGKS is as follows

$$[eqn]$$

The following new transformation functions are introduced:

$$[eqn]$$

Expanding the collision term in the BGK model, Equation (9) becomes as follows:

$$[eqn]$$ $$[eqn]$$

The flux $[eqn]$ across the cell interface is calculated as follows:

$$[eqn]$$

It should be pointed out that in this paper, the calculation of interface flux is based on the distribution function at the time $[eqn]$ .

2.2. Three Methods for Calculating Interface Flux

2.2.1. Original DUGKS

For the original DUGKS [18], the flux is calculated according to the midpoint method, as follows:

$$[eqn]$$

n is the outward unit vector normal to the cell interface ∂V_C_. $[eqn]$ represents the distribution function of midpoint x**b in the cell interface at the time $[eqn]$ .

Within half a time step (t^n^ to t^n^^+1/2^), Equation (1) is integrated along the characteristic line (the midpoint x**b of the cell interface is regarded as the end point of the characteristic line), and the collision term and force term are treated by the trapezoidal method, thus obtaining the following:

$$[eqn]$$

Similarly, the following new transformation functions are introduced:

$$[eqn]$$ $$[eqn]$$

Then, Equation (14) can be simplified as follows:

$$[eqn]$$

Using Taylor expansion, the right term of Equation (17) can be approximated as follows:

$$[eqn]$$

The auxiliary function $[eqn]$ at the midpoint of the interface and its gradient $[eqn]$ can be calculated by linear interpolation approximation, for example for the following directional uniform grid:

$$[eqn]$$

The function $[eqn]$ at the cell center x**C can be obtained by the following transformation:

$$[eqn]$$

$[eqn]$ can be calculated according to Equations (17)–(20). Based on the conservation of mass and momentum, the density and velocity at the center of the cell interface are as follows:

$$[eqn]$$

$[eqn]$ is obtained according to Equation (4). $[eqn]$ is determined according to Equation (15), and then, the flux $[eqn]$ can be obtained according to Equation (13). Then, according to Equation (8), the distribution function $[eqn]$ tracked in the DUGKS can be updated to the next time step. Since $[eqn]$ has been obtained according to Equation (20), $[eqn]$ is obtained by the following equation:

$$[eqn]$$

The evolution steps of the original DUGKS to update to the next time step are as follows:

$$[eqn]$$

In this paper, for the original DUGKS, optimized DUGKS, and Simpson–DUGKS, the macroscopic quantity (density and velocity) at the cell center can be obtained by tracking the following distribution function $[eqn]$ :

$$[eqn]$$

In all flow regimes, the time step $[eqn]$ is independent of the relaxation time $[eqn]$ and is determined by the Courant–Friedrichs–Lewy (CFL) condition [18]:

$$[eqn]$$

$[eqn]$ represents the minimum grid spacing. The maximum discrete speed $[eqn]$ is set to $[eqn]$ = $[eqn]$ . The CFL number is defined as CFL = CΔt/Δxmin. The CFL condition ensures that information propagating through the computational domain (e.g., waves or disturbances) does not travel further than one grid cell per time step. If the CFL condition is violated (i.e., the CFL number is too large), numerical instabilities can arise, leading to non-physical oscillations or the divergence of the solution. As demonstrated by the authors of [58], the stability of DUGKS is not affected by the CFL number apparently as long as CFL < 1. So, the CFL number often lies in 0 and 1.

Taking the two-dimensional uniform grid as an example, as shown in Figure 1, the evolution process of the distribution function of the original DUGKS can be written as follows:

$$[eqn]$$

where $[eqn]$ , $[eqn]$ . The distribution function $[eqn]$ is abbreviated as $[eqn]$ . The distribution function $[eqn]$ at the midpoint x_bx_ of the cell interface can be obtained by interpolating the distribution functions of the surrounding six cell centers:

$$[eqn]$$

where

$$[eqn]$$

Then, the evolution steps of $[eqn]$ updating to the next time step in the original DUGKS in the two-dimensional uniform grid are as follows:

$$[eqn]$$

2.2.2. Optimize DUGKS

For the optimized DUGKS [57], the trapezoidal method is used to calculate the interface flux:

$$[eqn]$$

$[eqn]$ represents the distribution function at the intersection point (node) x**N of the cell interface at the time $[eqn]$ , which is abbreviated as $[eqn]$ .

Taking the two-dimensional uniform grid as an example, as shown in Figure 2, the calculation formula of interface flux $[eqn]$ in the optimized DUGKS is as follows:

$$[eqn]$$

where $[eqn]$ , $[eqn]$ .

Next, the calculation of distribution function $[eqn]$ is introduced.

Using Taylor expansion, $[eqn]$ can be obtained from Equation (1) (see the Appendix A) as follows:

$$[eqn]$$

$[eqn]$ and $[eqn]$ can be obtained by the following distribution function of the surrounding four cell centers:

$$[eqn]$$

where $[eqn]$ .

$[eqn]$ can be calculated after $[eqn]$ is known, as follows:

$$[eqn]$$

$[eqn]$ is obtained from Equation (4), where the density and velocity are as follows:

$$[eqn]$$

The evolution process of the distribution function of optimization DUGKS in the two-dimensional uniform grid can be written as follows:

$$[eqn]$$

Then, the evolution steps of updating to the next time step in the optimization DUGKS are as follows:

$$[eqn]$$

2.2.3. Simpson–DUGKS

For Simpson–DUGKS, the Simpson method (three-point method) is proposed to calculate the interface flux, as follows:

$$[eqn]$$

Taking the two-dimensional uniform grid as an example, as shown in Figure 2, the calculation formula of interface flux $[eqn]$ in the Simpson–DUGKS is as follows:

$$[eqn]$$

where $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ .

The evolution process of the distribution function of Simpson–DUGKS in the two-dimensional uniform grid can be written as follows:

$$[eqn]$$

Then, the evolution steps of updating to the next time step in the Simpson–DUGKS are as follows:

$$[eqn]$$

2.3. Comparative Analysis

As shown in Figure 1 and Figure 2, the number of the nodes x_N_ is (Nx+1) (Ny+1), the number of x_bx_ in the x direction is Ny (Nx+1), the number of x_by_ in the y direction is Nx (Ny+1), and the number of x_C_ is NxNy. In the original DUGKS, the flux is calculated according to the distribution function at the cell interface center x**b. In the optimized DUGKS, the flux is calculated according to the distribution function at the intersection point (node) x**N. In the Simpson–DUGKS, the flux is calculated according to the distribution function of x**b and x**N. Therefore, the total number of points is NxNy + Ny (Nx+1) + Nx (Ny+1) in the original DUGKS, the total number of points is NxNy + (Nx+1) (Ny+1) in the optimized DUGKS, and the total number of points is 4NxNy+2 (Nx + Ny) + 1 in the Simpson–DUGKS. It can be seen that the optimized DUGKS has the minimum computational cost, so the optimized DUGKS has the highest computational efficiency in theoretical analysis.

There are three methods to evaluate the interface flux, as follows: the original DUGKS adopts the midpoint method; the optimize DUGKS adopts the trapezoidal method; and Simpson–DUGKS adopts the Simpson method (three-point method). The truncation errors of these three integration methods are as follows:

$$[eqn]$$

From Equation (40), it can be seen that the local truncation error of midpoint method is twice that of trapezoidal method. The algebraic accuracy of midpoint method and trapezoidal method is 1, and the algebraic accuracy of Simpson method is 3. Both the midpoint method and the trapezoidal method introduce spurious dissipation due to integration errors, smoothing out small-scale vortex structures. Moreover, they may also distort the phase characteristics of waves, affecting the simulation of transient flows. The Simpson method reduces these errors by means of high-order polynomial interpolation and lower truncation errors, and it has a significant effect, especially when simulating problems, such as unsteady and transient flows.

In the next section, the computational efficiency and error of these three methods are compared by numerical tests, as well as the stability.

3. Numerical Test

The Couette flow, Poiseuille flow, Taylor–Green vortex flow, and lid-driven cavity flow are numerically tested. In the simulations, $[eqn]$ , $[eqn]$ , $[eqn]$ , and ν = HU0/Re. The characteristic height H is fixed as H = 1, and the reference velocity U0 is fixed as U0 = 0.1. The numerical experiments demonstrate that the maximum CFL number available for original and optimized DUGKS are around 1.0 and 1.7 [57], and the CFL number often lies in 0 and 1.0. Therefore, the CFL number is selected as CFL = 0.1, 0.5, 0.95, 1.0, 1.3, 1.5, and 1.7 in this work. Considering the computational cost, the grid size including the low-resolution grid, medium-resolution grid, and high-resolution grid is selected as N = 16, 32, 64, and 128. The grid number N, CFL number, and Re number are adjusted to compare the computational efficiency, accuracy, and numerical stability under different conditions. To compare the numerical stability, the ratio of time step to relaxation time (∆t/τ = 5√2CFLRe/(3N)) is considered. So, a high Re number (Re = 10,000) is selected for the numerical tests, and some commonly used Re number (Re = 100 and 1000) is also selected.

The convergence criteria for attaining the steady-state solution are defined as follows:

$$[eqn]$$

where $[eqn]$ represents the velocity in the fluid domain at time t = nΔt. error represents the convergence error. Unless stated otherwise, the convergence criterion is set to be error < 10^−6^.

To assess the accuracy, the L2 error of macroscopic physical quantity is estimated as follows:

$$[eqn]$$

$[eqn]$ is the analytical solution.

To keep the boundary conditions consistent when comparing the three methods, the non-slip solid wall boundary condition adopts a bounce-back scheme to deal with the distribution function $[eqn]$ at the center of the ghost cells (i = 0 or Nx+1, j = 0 or Ny+1), as shown in Figure 1 and Figure 2.

It should be noted that all the simulations are performed on an Intel(R) Core(TM) i9-10900K CPU and 64-bit operating system.

3.1. Couette Flow

In the Couette flow, the lower wall is fixed, and the upper wall moves with horizontal velocity U0 = 0.1. Periodic boundary conditions are adopted for inlet and outlet (left and right boundaries). The analytical solution of horizontal velocity u(y) in the Couette flow is as follows:

$$[eqn]$$

The CPU time, L2 error, and the number of steps to reach convergence calculated by the three methods are shown in Table 1. Note that the symbol “-” indicates that the calculation does not converge, Δ1 represents the percentage increase in the calculation time of Simpson–DUGKS compared with the original DUGKS, and Δ2 represents the percentage decrease in the calculation time of optimized DUGKS compared with the original DUGKS. In terms of accuracy, the L2 errors calculated by the three methods according to Equation (42) are the same. In terms of computational efficiency, although the three methods have the same number of steps to achieve the convergence criterion of Equation (41), the computational time of optimized DUGKS is reduced by about 18–47% compared with the original DUGKS, and the computational time of Simpson–DUGKS is increased by about 27–120% compared with the original DUGKS. In terms of computational stability, the three methods have the same stability under different CFL numbers and grid numbers in the test. All three methods can keep the computational convergence, except for the divergence at Nx = Ny = 128 with CFL = 1.7.

Furthermore, the numerical performance of the three methods with different CFL numbers (CFL = 0.95, 1.5, 1.7) and different grid numbers (Nx = Ny = 16, 32, 64, 128) is tested at a high Re number. The results are shown in Table 2. With CFL = 0.95 and 1.5, the three methods do not converge with the grid numbers Nx = Ny = 16, 32, and 64 but converge stably only with the grid number Nx = Ny = 128. With CFL = 1.7, the three methods cannot work with the grid number Nx = Ny = 16, 32, and 64. Only with the grid number Nx = Ny = 128, the calculation can converge stably. In terms of computational efficiency, although the steps to reach the convergence criterion of Equation (41) are the same among the three methods, the computational time of the optimized DUGKS is reduced by about 14–39% compared with that of the original DUGKS.

In terms of accuracy, the L2 errors calculated by the three methods according to Equation (42) are the same. The corresponding velocity distribution is shown in Figure 3, which is in good agreement with the analytical solution.

3.2. Poiseuille Flow

The Poiseuille flow is driven by external force ρa_x_ (a_y_ = 0), and the inlet and outlet (left and right boundaries) adopt periodic boundary conditions. The analytical solution of the horizontal velocity u(y) in the Poiseuille flow can be expressed as follows:

$$[eqn]$$

where U0 = a_x_H^2^/8ν.

The CPU time, L2 error, and the number of steps to reach convergence calculated by the three methods are shown in Table 3. Note that Δ1 represents the percentage increase in the calculation time of Simpson–DUGKS compared with the original DUGKS, and Δ2 represents the percentage decrease in the calculation time of optimized DUGKS compared with the original DUGKS. In terms of accuracy, the L2 errors calculated by the three methods according to Equation (42) are the same. In terms of computational efficiency, although the three methods have the same number of steps to achieve the convergence criterion of Equation (41), the computational time of optimized DUGKS is reduced by about 20–39% compared with the original DUGKS, and the computational time of Simpson–DUGKS is increased by about 35–62% compared with the original DUGKS. In terms of computational stability, the three methods have the same stability under different CFL numbers and grid numbers in the test. The three methods can keep the computational convergence, except for the divergence with the grid number Nx = Ny = 128 and CFL = 1.7.

Furthermore, the numerical performance of the three methods with different grid numbers is tested at high Re number and large CFL number, as shown in Table 4. With the grid number Nx = Ny = 16, 32, and 64, the three methods cannot work. Only with the grid number Nx = Ny = 128, the calculation can converge stably. In terms of computational efficiency, although the steps to reach the convergence criterion of Equation (41) are the same among the three methods, the computational time of the optimized DUGKS is reduced by about 44% compared with that of the original DUGKS. In terms of accuracy, the L2 errors calculated by the three methods according to Equation (42) are the same. The corresponding velocity distribution is shown in Figure 4, which is in agreement with the analytical solution. Zooming in on this area can show that the predictions of the three methods noticeably deviate from the theoretical values near the center.

3.3. Taylor–Green Vortex Flow

The Taylor–Green vortex flow is characterized by turbulence and unsteady in macroscopic view. It is a simple model for studying the decomposition process of larger-scale vortexes gradually declining into smaller-scale vortexes and the resulting uniform isotropic turbulence. It is widely used as a standard example of a numerical algorithm for fluid solutions. In the two-dimensional Taylor–Green vortex flow with periodic boundary conditions, the analytical solutions of its velocity and density are as follows:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where u and v denote horizontal and vertical velocities, respectively. ρ denotes the density. $[eqn]$ and $[eqn]$ represent the wave numbers in the horizontal and vertical direction after expansion through Fourier space, respectively. $[eqn]$ represents the wave number in Fourier space. υ denotes the kinematic viscosity of fluid. U0 denotes the initial velocity peak. ρ0 denotes the average density, and it is set as ρ0 = 1.

The convergence error, L2 error, and CPU time (s) calculated by the three methods to 10,000 steps at Re = 1000 are tested, as shown in Table 5, Table 6, Table 7 and Table 8. Note that Δ1 represents the percentage increase in the calculation time of Simpson–DUGKS compared with the original DUGKS, and Δ2 represents the percentage decrease in the calculation time of optimized DUGKS compared with the original DUGKS. In terms of computational stability, all three methods can calculate stably under different grid numbers in the test. In terms of accuracy, all three methods are in good agreement with the analytical solution, as shown in Figure 5 and Figure 6. However, as shown in Table 5, Table 6 and Table 7, the convergence error calculated by Equation (41) and the L2 error calculated by Equation (42) are different. The convergence error of optimized DUGKS is the smallest under the same number of grids. Compared with the analytical solution of velocity, the L2 error of the optimized DUGKS is the smallest, and the L2 error of the original DUGKS is the largest under the same grid number. Compared with the analytical solution of density, the L2 error of optimized DUGKS is the largest under the same number of grids. As shown in Table 8, in terms of computational efficiency, when the number of calculation steps is 10,000, the calculation time of optimized DUGKS is reduced by about 12–24% compared with the original DUGKS, and the calculation time of Simpson–DUGKS is increased by about 22–30% compared with the original DUGKS. As shown in Figure 7, all three methods are of second-order accuracy in space and time.

3.4. Lid-Driven Cavity Flow

As a standard example, the lid-driven cavity flow contains complex flows, and the results of Ghia et al. [59] are usually used as reference data. The upper wall moves with horizontal velocity U0 = 0.1, and the lower wall and left- and right-side walls are fixed.

The CPU time calculated by the three methods is shown in Table 9. Note that the symbol “-” indicates that the calculation does not converge, Δ1 represents the percentage increase in the calculation time of Simpson–DUGKS compared with the original DUGKS, and Δ2 represents the percentage decrease in the calculation time of optimized DUGKS compared with the original DUGKS. In terms of computational efficiency, compared with the original DUGKS, the computational time of the optimized DUGKS can be reduced by about 77% at most. It should be noted that, when CFL = 0.1 and N = 16, the calculation time of optimized DUGKS increases by about 40% compared with the original DUGKS. When CFL = 0.95 and N = 16, the calculation time of Simpson–DUGKS is reduced by about 59% compared with the original DUGKS. In terms of stability, the optimized DUGKS can still be computationally stable at large CFL numbers (CFL = 1.5 and 1.7). In terms of accuracy, the horizontal velocity of the vertical central axis obtained by the three methods is compared with the reference data, as shown in Figure 8, Figure 9 and Figure 10. It can be seen that the result of optimized DUGKS deviates more from the result of Ghia et al. [59], especially in the case of a small grid number. With CFL = 0.95, the horizontal velocity of the vertical central axis is shown in Figure 11, Figure 12, Figure 13 and Figure 14. With N = 64 and 128, the distribution of horizontal velocity obtained by the three methods is roughly similar, but with N = 16 and 32, the distribution of horizontal velocity obtained by optimized DUGKS is different from the other two methods, especially when N = 16, the error of horizontal velocity obtained by optimized DUGKS is obviously large. As shown in Figure 15, Figure 16, Figure 17 and Figure 18, the results of the optimized DUGKS contain some oscillations at the corner, while the solution of the original DUGKS and Simpson–DUGKS is smoother. What is more, the average density obtained by the Simpson–DUGKS is closest to the initial density (ρ0 = 1). The density error exists in the result obtained by the optimized DUGKS. The total number of points in the optimized DUGKS is the least. For the lid-driven cavity flow, it seems that each point plays an important role in the performance of numerical method. So, the minimum number of points may result in the large error.

Furthermore, the performance of the three methods is tested under CFL = 0.1, 0.95, 1, 1.5, and 1.7 at Re = 10,000. The three methods cannot work when CFL = 1, 1.5, and 1.7. The three methods can converge when CFL = 0.1 and 0.95, and the calculation time to reach the convergence criterion is shown in Table 10. In terms of computational efficiency, compared with the original DUGKS, the computational time of the optimized DUGKS can be reduced by about 55% at most. In terms of stability, when CFL = 0.1 and N = 32, the optimized DUGKS cannot converge. As shown in Figure 19 and Figure 20, no non-physical oscillation phenomenon is found at the high Re number, and the velocity contours obtained by original DUGKS and Simpson–DUGKS are almost consistent, but the velocity contours obtained by optimized DUGKS are different, which seems to deviate from the reference results. As shown in Figure 21 and Figure 22, the results of the optimized DUGKS contain some oscillations at the corner, while the solution of the original DUGKS and Simpson–DUGKS is smoother.

The lattice BGK (LBGK) model and the multiple relaxation time (MRT) model, as popular standard LBM, have been successfully applied and well-accepted for incompressible NS solutions [60]. The pure stability of the LBGK, MRT, and Simpson–DUGKS method is evaluated without considering the accuracy of the results. To reach a steady-state solution at Re = 1000, the minimum required mesh resolution of the LBGK, MRT, and Simpson–DUGKS method is the 70 × 70, 20 × 20, 10 × 10 uniform mesh, respectively. On a uniform 128 × 128 mesh, the computations from the LBGK and MRT blow up at Re = 1900 and Re = 8000, respectively. However, the Simpson–DUGKS method works even at Re = 100,000. Clearly, in comparison with standard LBM, the Simpson–DUGKS method has super performance in stability.

Moreover, the average density obtained by the Simpson–DUGKS is closest to the initial density (ρ0 = 1). The density error exists in the result obtained by the optimized DUGKS. The reasons for the large density error can be explained as follows:

• (1)Spatial and temporal truncation errors: optimized DUGKS can introduce truncation errors, which may accumulate and lead to significant density errors in specific flow regimes, especially when the flow has complex structures or high gradients.
• (2)Dispersion and dissipation errors: Optimized DUGKS can suffer from dispersion and dissipation errors. In some flows, these errors might be amplified, causing the density waves to propagate at incorrect speeds or with distorted shapes. This can lead to a discrepancy between the numerical solution and the exact solution, manifesting as a large density error.
• (3)Complex flow structures: flow with complex structures, such as those involving multiple vortices, can pose challenges for optimized DUGKS. The method might have difficulties in accurately resolving these intricate structures, leading to errors in the density field. For example, in a flow with contact discontinuities, the optimized DUGKS may not capture the sharp density gradients at the discontinuities, resulting in a smeared density profile and larger errors.

4. Conclusions

In this paper, the Simpson method is proposed to calculate the interface flux according to the distribution function at the node and the midpoint of the cell interface, and the optimized DUGKS and Simpson–DUGKS, considering the force term, are derived. Then, the original DUGKS, optimized DUGKS, and Simpson–DUGKS are simply compared and analyzed in theory. Finally, through numerical tests, the following conclusions are obtained:

In the numerical test of steady unidirectional flow (Couette flow and Poiseuille flow), the L2 errors calculated by the three methods are consistent, and the stability of the three methods is consistent under different CFL numbers. In terms of computational efficiency, among the three methods, optimized DUGKS has the least computational time, and Simpson–DUGKS has the most computational time, which also verifies the theoretical analysis and comparison. In more complex flows, some findings should be noted:

• (1)In the numerical test of Taylor–Green vortex flow, the stability of the three methods is consistent under different grid numbers. Compared with the analytical solution of velocity, the L2 error of optimized DUGKS is the smallest under the same number of grids, but the L2 error of optimized DUGKS is the largest under the same number of grids compared with the analytical solution of density. In terms of computational efficiency, among the three methods, optimized DUGKS has the least computational time, while Simpson–DUGKS has the most computational time.
• (2)In the numerical test of the lid-driven cavity flow, the calculation time of the optimized DUGKS can be reduced by about 77% at most, compared with the original DUGKS. It should be noted that, when CFL = 0.1 and N = 16, the calculation time of optimized DUGKS increases by about 40% compared with the original DUGKS. When CFL = 0.95 and N = 16, the calculation time of Simpson–DUGKS is reduced by about 59% compared with the original DUGKS. In terms of stability, the optimized DUGKS can still be computationally stable at large CFL numbers (CFL = 1.5 and 1.7). In terms of accuracy, the results of optimized DUGKS deviate more from the reference results, especially in the case of a small grid number. For example, when CFL = 0.95 and N = 16, the error of horizontal velocity contours obtained by optimized DUGKS is obviously large. This may be explained by the fact that the optimized DUGKS requires the fewest number of auxiliary points.

The guidelines on how to select appropriate DUGKS flux calculation methods are shown as follows: The optimized DUGKS is recommended for steady flow (Couette/Poiseuille flow, etc.) where efficiency is critical, which achieves 40–77% faster computation at large CFL. The Simpson–DUGKS is preferred for unsteady/vortical flow (Taylor–Green/cavity flow, etc.) requiring a high resolution of gradients, which reduces L2 errors by 30–50% with low grid resolution at moderate CFL, ensuring excellent stability.