Novel Schemes of No-Slip Boundary Conditions for the Discrete Unified Gas Kinetic Scheme Based on the Moment Constraints

The boundary conditions are crucial for numerical methods. This study aims to contribute to this growing area of research by exploring boundary conditions for the discrete unified gas kinetic scheme (DUGKS). The importance and originality of this study are that it assesses and validates the novel schemes of the bounce back (BB), non-equilibrium bounce back (NEBB), and Moment-based boundary conditions for the DUGKS, which translate boundary conditions into constraints on the transformed distribution functions at a half time step based on the moment constraints. A theoretical assessment shows that both present NEBB and Moment-based schemes for the DUGKS can implement a no-slip condition at the wall boundary without slip error. The present schemes are validated by numerical simulations of Couette flow, Poiseuille flow, Lid-driven cavity flow, dipole–wall collision, and Rayleigh–Taylor instability. The present schemes of second-order accuracy are more accurate than the original schemes. Both present NEBB and Moment-based schemes are more accurate than the present BB scheme in most cases and have higher computational efficiency than the present BB scheme in the simulation of Couette flow at high Re. The present Moment-based scheme is more accurate than the present BB, NEBB schemes, and reference schemes in the simulation of Poiseuille flow and dipole–wall collision, compared to the analytical solution and reference data. Good agreement with reference data in the numerical simulation of Rayleigh–Taylor instability shows that they are also of use to the multiphase flow. The present Moment-based scheme is more competitive in boundary conditions for the DUGKS.


Introduction
Discrete Unified Gas Kinetic Scheme (DUGKS) for all Knudsen Number Flows was proposed [1], which combined the advantages of the Lattice Boltzmann method (LBM) and the Unified Gas Kinetic Scheme (UGKS). The DUGKS and LBM can share the same equilibrium distribution function, and both start from the Boltzmann equation. Although the computational time cost of DUGKS is slightly higher than that of LBM in a continuous flow, the numerical stability of DUGKS is better than that of LBM, and it can capture flow characteristics that LBM is not competent to [2,3]. Both DUGKS and UGKS have asymptotic preserving properties, where the time step is not limited by the collision time. For the reconstruction of cell-interface flux, the UGKS adopts the analytical time-evolved integral solution and the DUGKS adopts a simpler numerical characteristic solution, in which the calculation time cost is approximately 70% of that in the UGKS [4]. In short, the DUGKS has a simple framework and the powerful ability to capture the flow characteristics in a wide flow regime, which can make the DUGKS a competitive tool in comparison with LBM and UGKS. Recently, the DUGKS has been successfully applied to a variety of flow problems in different flow regimes, such as turbulent flows [5][6][7][8], micro flows [9][10][11][12][13], compressible flows [4,[14][15][16][17], multiphase flows [18][19][20], fluid-particle flows [21,22], flows with complex geometries [23,24], gas mixture systems [25,26], pore-scale porous media flows [27], and plasma physics [28,29]. In addition to flow problems, the DUGKS was also extended to multiscale transport problems such as phonon heat transfer [30][31][32] and the radiation of photons [33,34].
The boundary conditions are crucial for numerical methods [35][36][37]. If boundary conditions are not properly introduced in the numerical methods, severe problems may arise, such as the divergence of numerical solutions. In short, boundary conditions should be treated with great care. Because the LBM and DUGKS share a common kinetic origin, researchers tried to introduce the boundary conditions from the LBM into the DUGKS. For example, two kinds of no-slip boundary conditions are widely used: The bounce-back (BB) scheme [1] and the non-equilibrium bounce-back (NEBB) scheme [38]. They have received considerable critical attention. For the LBM and DUGKS, there are distinctive modeling differences in the particle evolution process. The LBM separates the particle streaming and collision process. However, particle transport and collision are fully coupled in DUGKS. This dynamic difference can result in the deviation of the boundary condition in the LBM and DUGKS. For example, the boundary condition is processed at time t + 0.5∆t in the DUGKS but at time t in the LBM. Yang et al. [39] have analyzed and assessed the BB scheme and the NEBB scheme for the DUGKS. Although the BB scheme is a simple and common method to apply the no-slip boundary condition, it introduces an additional error (a purely artificial numerical slip error) into the DUGKS. The NEBB scheme eliminates the error of the numerical slip, but the closure to find the unknowns at the wall is somewhat arbitrary [40]. There are only a few studies on boundary conditions for the DUGKS, which have not been studied thoroughly. So, it can bring renewed interest and value to introduce and test other boundary conditions for the DUGKS.
Recently, the moment-based boundary condition has received increased interest and attention [40][41][42][43][44][45], which is based on the moments of the LBM [40] and has not yet been introduced to the DUGKS. Numerical simulations show that the moment-based boundary condition converges with second-order accuracy using dipole-wall collisions [40], natural convection in the square cavity [41], and lid-driven cavity flow [43][44][45]. A numerical simulation of the Poiseuille flow shows that the moment-based boundary condition can eliminate the spurious oscillations seen in solutions using other boundary conditions, considering the nonzero deviatoric stress [42]. The moment-based boundary condition can be parallelized easily for efficient computing, and it is conceptually simpler and more straightforward than other methods that involve a mixture of the bounce-back rule, hydrodynamic moments, momentum corrections, and other modifications to the distribution functions [44].
To the best of the authors' knowledge, the moment-based boundary condition for the DUGKS has not been analyzed. So, the present work introduces the moment-based boundary condition to the DUGKS, and its numerical performance is studied in comparison with the BB scheme and the NEBB scheme. Generally, for boundary conditions in the DUGKS, the boundary nodes are located at the cell interface, and the unknown original distribution functions are treated at time t n+1/2 = t n + 0.5∆t [39]. In the DUGKS, the original distribution functions are obtained after the transformed distribution functions are calculated. So, in this study, novel schemes of the BB, NEBB, and Moment-based boundary conditions are proposed for the DUGKS, which treat the unknown transformed distribution functions at time t n+1/2 = t n + 0.5∆t. That is, the implementation of the boundary condition shifts from the original f (x w , ξ α , t + h) to the transformed f (x w , ξ α , t + h), which may improve the numerical performance.
Assessment and validation of the novel schemes are necessary and important to their application. To test boundary conditions for the DUGKS, numerical simulations have only focused on the Couette flow and the Poiseuille flow in Ref. [39]. Considering the generality, more complex flows should be adopted to test boundary conditions, such as the multiphase flow. The multiphase flows are of both academic and industrial interest [46][47][48][49][50][51][52], such as the unsteady Rayleigh-Taylor instability. The benchmark problem of Rayleigh-Taylor instability has been used to validate the newly developed numerical approaches [18,19,50,51]. The bounce-back boundary conditions are applied to the bottom and top wall boundaries in the DUGKS simulations of Rayleigh-Taylor instability [18]. However, to the best of the authors' knowledge, the effect of varied wall boundary conditions on the numerical simulation of Rayleigh-Taylor instability has been not investigated. Therefore, in this work, the unsteady Rayleigh-Taylor instability is used to test and validate the proposed schemes of boundary conditions, in addition to the Couette flow, the Poiseuille flow, the dipole-wall collision, and the lid-driven cavity flow.
The rest of the paper is organized as follows. In Section 2, the DUGKS with a force term, the original and novel schemes of BB, NEBB, and Moment-based boundary conditions for the DUGKS are present. In Section 3, the boundary conditions are tested and assessed using Couette flow, Poiseuille flow, Lid-driven cavity flow, vortex dipole-wall collision, and Rayleigh-Taylor instability. In Section 4, the conclusions are drawn.

Numerical Method
It is noted that the present work is for isothermal continuum flow.

DUGKS with a Force Term
In the LBM, a modeled gas, which is composed of identical particles whose velocities are restricted to a finite set of vectors, is considered. Similar to the LBM, the DUGKS follows the lattice Boltzmann equation. It is inevitably necessary to introduce a force term to the DUGKS in some cases, such as the external force driven by Poiseuille flow. The governing equation of the DUGKS with a force term [3]: Equation (1) describes the spatial and temporal evolution of the distribution f (ξ α , x, t) of particles with velocity ξ α at position x and time t. The fundamental variable in gas kinetic theory is the particle distribution function, which simultaneously represents the density of mass in both physical space and velocity space.
The force term is approximated as [3] where a denotes the acceleration due to the external force, and R and T represent the gas constant and the temperature, respectively. For simplicity, the particle velocity ξ and the fluid velocity u have been normalized by √ 3RT, giving a sound speed of c s = 1/ √ 3 [53]. So, in the following simulations of isothermal continuum flow, the temperature T is constant with c s 2 = RT = 1/3. With the Taylor expansion at approximately zero particle velocity at a low Mach number, the Maxwellian equilibrium distribution function f eq is approximated as [1]  In the widely used D2Q9 model [5], the velocity space {ξ α } and the corresponding weights {W α } are given as The computational domain is divided into the cells V i,j = ∆x i ∆y j centered at (x i , y j ) in the DUGKS with the D2Q9 model. As a finite volume scheme, the cell-averaged values of the distribution function and the force term are introduced, denoted as f n and S n , Integrating Equation (1) into each cell V i,j from time t n = n∆t to time t n+1 = (n + 1)∆t, and using the midpoint rule and trapezoidal rule, the evolution equation of the DUGKS is advanced from t n to t n+1 as, The scalar variable F n+1/2 represents the flux across the cell interface, which is computed as (9) where f (ξ α , x b , t n+1/2 ) represents the distribution at the cell interface x b at the time t n+1/2 = t n + h (h = ∆t/2), and n is the outward unit vector normal to the surface ∂V i,j . For the purpose of eliminating the simplicity, new distribution functions are introduced: The collision term can be expanded in the BGK collision model, and Equation (10) can be converted to the following equations: The evolution equation of the DUGKS is simplified as All other forms of the distribution function can be expressed in terms of f in the DUGKS. So, we focus on the distribution function f .
As can be seen from Equation (13), the critical step is to evaluate the interface flux F n+1/2 . For DUGKS with higher-order accuracy, more intermediate time steps can be selected, for example, the flux at the cell interface at t * = t n + ∆t/6 and t * * = t n + 3∆t/4 need calculating [54]. Considering the easy implementation and fast computation in the present study, we select the intermediate time step t n + ∆t/2 in the present study. Equation (1) is integrated within a half time step (h = ∆t/2) along the characteristic line using the trapezoidal rule to treat the collision and force terms, and we have where x b at the cell interface denotes an end point along the characteristic line. Similar to the above discussion, new distribution functions are introduced: With new distribution functions, Equation (14) is simplified as With the Taylor expansion [1], the right term of Equation (17) can be approximated as where f + (ξ α , x b , t n ) and its gradients ∇ f + (ξ α , x b , t n ) can be approximated by linear interpolations. For the uniform mesh in the x-direction, e.g., The distribution function f + (ξ α , x, t n ) at the cell center can be calculated as Then, f (ξ α , x b , t n+1/2 ) can be obtained by using Equations (17) and (18). Similarly, the original distribution function f (ξ α , The equilibrium function f eq (ξ α , x b , t n+1/2 ) at the cell interface can be determined by the conserved variables. Based on the conservation of mass and momentum, the density and velocity at the cell interface center can be given as Finally, the flux F n+1/2 is obtained by Equation (9) after f (ξ α , x b , t n+1/2 ) is determined by Equation (15). Then, the tracked distribution function f can be updated to the next time step by Equation (13).
Particularly, the following equation will be used in the DUGKS [3], The density ρ and velocity u at the cell center can be calculated from f : In the present DUGKS, the relaxation time τ is computed from τ = ν/RT, where ν is the kinematic viscosity. The time step ∆t is independent of the relaxation time τ for all flow regimes, and it is determined by the Courant-Friedrichs-Lewy (CFL) condition [55]: where χ is the CFL number, ∆x is the minimum grid spacing, and C is the maximum discrete velocity. We set C to be √ 2 in the present work.

Original Schemes of No-Slip Boundary Conditions
The original schemes of BB, NEBB, and Moment-based boundary conditions are introduced, and they are analyzed theoretically from the view of the moment.
The boundary conditions are applied to solve the unknown distributions that have velocities pointing into the fluid domain. So, we need to find the unknown distributions.

Moment-Based Scheme
In brief, the Moment-based scheme states that the unknown distribution functions at the boundary can be determined by the linearly independent moments. For the D2Q9 model, there are nine discrete velocity moments, including the 6 hydrodynamic moments (the density ρ, momentum ρu, and momentum flux Π) and 3 non-hydrodynamic moments (Q xxy , Q xyy , and S xxyy ) [40]: For the no-slip wall boundary, we pick three linearly independent moments (momentum Π x , Π y , and momentum flux Π xx ), and impose constraints upon them, i.e., u x = u wx (const), u y = 0, ∂ x u x = 0 for the south boundary (u wy = 0). Then, by Chapman-Enskog, We obtain the following three conditions for three unknowns: For the south boundary (u wy = 0), solving these equations yields For the south boundary (u wy = 0), the Non-equilibrium Bounce-Back (NEBB) scheme can be expressed as [38]: where u wx denotes the wall velocity in the x-direction.
It is found that the calculation of the density in Equation (28) is the same as that in Equation (29) when u wy = 0. Then, from a fresh perspective based on the moment, it is found that the conditions of the NEBB scheme are That is, for the NEBB scheme, three linearly independent moments are picked for the south boundary: Momentum Π x , Π y , and momentum flux Q xxy . Π x and Π y are useful in defining the no-slip condition, but the condition on a component of the third-order moment Q xxy seems somewhat arbitrary.

Bounce-Back Scheme
The Bounce-Back (BB) scheme assumes that a particle just reverses its velocity after hitting the wall. The unknown distribution functions f (x w , ξ α , t n+1/2 ) for particles leaving the wall are determined as [1], where u w denotes the wall velocity and ρ w denotes the density at the wall boundary, which can be approximated well by the constant average density for nearly incompressible flows [1]. For the south boundary (the wall velocity at y-direction u wy = 0), the unknown original distribution functions f 2 , f 5 , and f 6 are given as, where u wx denotes the wall velocity in the x-direction.
For the BB scheme, three linearly independent moments are picked for the south boundary: Π y , Q xyy , and Q xxy , leading to the conditions on a moment basis, Two components of the third-order moment Q xyy and Q xxy result in the numerical slip error.

Novel Schemes of No-Slip Boundary Conditions
In the DUGKS, new transformed distribution functions are introduced to evaluate the interface flux. So, the original scheme of boundary conditions should be transformed. The original distribution functions f (x w , ξ α , t + h) are obtained after the transformed distribution functions f (x w , ξ α , t + h) are calculated. So, the implementation of the boundary condition shifts from the original f (x w , ξ α , t + h) to the transformed f (x w , ξ α , t + h), which may improve the numerical performance. It is noted that the particle streaming and collision are fully coupled in the DUGKS.
As shown in Figure 1, in the present work, the boundary is located at a cell interface, and we treat the unknown distribution functions f (x w , ξ α , t + h) at the cell interface center x w . Owing to new distribution functions and force terms in the DUGKS, we should convert the original boundary conditions into a new format. The boundary conditions will be translated into constraints on f (x w , ξ α , t + h). It should noted that more details about derivations for the novel schemes are shown in Appendix A.

Analysis of Numerical Slip Errors of Novel Schemes
In the following, the numerical slip errors (u error = u x − u wx ) of the novel schemes of no-slip boundary conditions are theoretically analyzed.
Inspired by Guo et al. [56], He et al. [57], Wang et al. [58], and Yang et al. [59], the unidirectional steady flow is adopted to analyze the numerical error of boundary conditions. The assumptions in the unidirectional steady flow are written as where φ denotes an arbitrary flow variable. a x and a y denote the component of the acceleration a in the x-direction and y-direction, respectively. Considering Equations (14) and (35), we can obtain With the Maxwellian equilibrium distribution function known, as shown in Equation (4), we can obtain Yang et al. [39] treated the unknown distribution functions f (x w , ξ α , t + h) at the cell interface x w , and the numerical slip error of the BB scheme and NEBB scheme generated in the DUGKS can be written as However, in the present work, we treat the unknown distribution functions f (x w , ξ α , t + h) at the cell interface center x w .
For the Moment-based scheme, we substitute Equation (37) into Equation (35), and obtain For the NEBB scheme, we substitute Equation (39) into Equation (35), and obtain For the BB scheme, we substitute Equation (40) into Equation (35), and obtain With Equations (15) and (43), we can obtain It is found that the numerical slip error of the BB scheme in the present work is also zero (u error,BB = u x − u wx = 0) when the acceleration a x is equal to 0.
As shown in Equation (34), the only hydrodynamic condition is that the flow rate through the wall is zero in the BB scheme, and the velocity along the wall is undefined. The BB scheme does not impose a boundary condition on the tangential velocity at the wall, resulting in the numerical slip velocity when the external force term exists. Both the NEBB scheme and Moment-based scheme select momentum Π x and Π y , which can satisfy the no-slip condition without generating numerical slip error. So, momentum Π x and Π y are useful in defining the no-slip condition. It seems sensible to choose the hydrodynamic moments instead of the non-hydrodynamic ones that do not appear in these equations of motion [40].

Explore Extending Present Schemes to Curved Boundaries
To extend the current boundary treatment to curved boundaries, this work also tries to introduce the methods from the LBM into the DUGKS, including the link method (Ladd [60,61]), interpolated methods (Bouzidi et al. [62], Yu et al. [63]), and single-node schemes (Zhao&Yong [64], Tao et al. [65], Zhao et al. [66], Chen et al. [67]). Although the LBM and DUGKS share a common kinetic origin, we need to be aware of distinctive features in the LBM and DUGKS. It should be noted that the present boundary conditions treat the unknown transformed distribution functions f at the cell interface in the DUGKS.

Link Method
As shown in Figure 2, the link method approximates the curved solid boundary to a stair-like zigzag line (dash line), which is set in the center of the solid and fluid nodes (i.e., the cell interface). So, it is quite convenient for the DUGKS to apply the link method to treat the curved boundaries by treating the unknown transformed distribution functions f (x w , ξ α , t + h) at the cell interface center x w and at time t n+1/2 = t n + 0.5∆t (h = 0.5∆t). If the grid mesh is fine enough, it is easy and efficient to treat the curved boundaries by adopting the current boundary treatment (as shown in Section 2.3) directly. It is found that the numerical slip error of the BB scheme in the present work is also zero ( error,BB 0 x wx u u u = − = ) when the acceleration ax is equal to 0.
As shown in Equation (34), the only hydrodynamic condition is that the flow rate through the wall is zero in the BB scheme, and the velocity along the wall is undefined. The BB scheme does not impose a boundary condition on the tangential velocity at the wall, resulting in the numerical slip velocity when the external force term exists. Both the NEBB scheme and Moment-based scheme select momentum Πx and Πy, which can satisfy the no-slip condition without generating numerical slip error. So, momentum Πx and Πy are useful in defining the no-slip condition. It seems sensible to choose the hydrodynamic moments instead of the non-hydrodynamic ones that do not appear in these equations of motion [40].

Explore Extending Present Schemes to Curved Boundaries
To extend the current boundary treatment to curved boundaries, this work also tries to introduce the methods from the LBM into the DUGKS, including the link method (Ladd [60,61]), interpolated methods (Bouzidi et al. [62], Yu et al. [63]), and single-node schemes (Zhao&Yong [64], Tao et al. [65], Zhao et al. [66], Chen et al. [67]). Although the LBM and DUGKS share a common kinetic origin, we need to be aware of distinctive features in the LBM and DUGKS. It should be noted that the present boundary conditions treat the unknown transformed distribution functions f at the cell interface in the DUGKS.

Link Method
As shown in Figure 2, the link method approximates the curved solid boundary to a stair-like zigzag line (dash line), which is set in the center of the solid and fluid nodes (i.e., the cell interface). So, it is quite convenient for the DUGKS to apply the link method to treat the curved boundaries by treating the unknown transformed distribution functions f (xw, ξα, t + h) at the cell interface center xw and at time tn+1/2 = tn + 0.5Δt (h = 0.5Δt). If the grid mesh is fine enough, it is easy and efficient to treat the curved boundaries by adopting the current boundary treatment (as shown in Section 2.3) directly.

Figure 2.
Link method approximates curved boundary to zigzag one. The white square means fluid node xij, black square means solid node, solid line means real curved boundary, and dashed line Figure 2. Link method approximates curved boundary to zigzag one. The white square means fluid node x ij , black square means solid node, solid line means real curved boundary, and dashed line means numerical boundary. Gray circle denotes the cell interface center x b . Gray diamond denotes the numerical boundary node x w , which is located at the cell interface center on the numerical boundary.

Interpolated Method
In actuality, the link method makes the smooth curved boundary into the rough boundary. To improve the curved boundary conditions, some interpolated methods are developed. It should be noted that strategies to design interpolated schemes are not unique. For example, in the LBM, two representative interpolated bounce-back (IBB) schemes are the conditional scheme proposed by Bouzidi et al. [62] and the unified scheme by Yu et al. [63]. It is noted that the interpolated methods are based on the bounce-back scheme in the LBM. Similarly, the present BB scheme for the DUGKS can be applied to the interpolated methods. However, it needs more implementation effort or other strategies to extend the present NEBB scheme and Moment-based scheme to curved boundaries.

Conditional Scheme
Inspired by the conditional scheme [62], the IBB scheme in the DUGKS is proposed by using linear (first-order treatment) or quadratic (second-order treatment) interpolation formulas involving values at two or three nodes. The conditional scheme uses separate treatments for q < 0.5 and q ≥ 0.5 according to the scaled distance.
Taking the boundary configuration in Figure 3 as an example, we use quadratic interpolation to obtain, R PEER REVIEW 13 of 44 Figure 3. Sketch of a curved boundary located between two lattice nodes arbitrarily. Grey circles represent the fluid nodes (i.e., xf, xff, xfff), black circle represents the solid node (xs), and square box represents the intersection (xw) of the boundary and the grid line. ξα defines the lattice velocity of the particle, which travels from xf to xff and α ξ denotes the opposite direction of ξα

Unified Scheme
Inspired by the unified scheme [63], the present work treats the curved boundaries in three main steps: . Sketch of a curved boundary located between two lattice nodes arbitrarily. Grey circles represent the fluid nodes (i.e., x f , x ff , x fff ), black circle represents the solid node (x s ), and square box represents the intersection (x w ) of the boundary and the grid line. ξ α defines the lattice velocity of the particle, which travels from x f to x ff and ξ α denotes the opposite direction of When q = 1 2 , the conditional interpolation formulas are reduced to the "bounce-back" scheme. Since a distribution function travels precisely one grid spacing from t to t + ∆t, particles that start from x f can end precisely at the same location only when x f is half a grid spacing from the wall.

Unified Scheme
Inspired by the unified scheme [63], the present work treats the curved boundaries in three main steps: A virtual distribution function f (x w , ξ α , t + h) is interpolated from the distribution functions at x f , x ff and x fff . We use quadratic interpolation to obtain, All the present BB, NEBB, and Moment-based schemes (as shown in Section 2.3) are

Single-Node Scheme
The interpolated methods [62,63] improve the accuracy to be second-order. However, they involve at least two neighboring fluid nodes, e.g., x f , x ff , and x fff . It is inevitable to undermine the local computation property and further the parallel performance of LBM [68,69]. Moreover, the required number of fluid nodes may not be available for some cases, e.g., dense particle suspensions. Zhao and Yong (2017) [64] developed a singlenode second-order bounce-back scheme for curved boundaries, which used pre-and postcollision density distributions to determine the unknown density distributions at the boundary. An alternative single-node second-order scheme was proposed by   [65]. Zhao et al. (2020) [66] derived a general single-node second-order scheme for curved boundaries. However, the free parameter is limited within the range of [max(0, 2q − 1), 2q]. Later, it was reported that the free parameter γ can be arbitrarily chosen within the range of [0, 2q] in a general single-node second-order scheme proposed by Chen et al. (2021) [67], but γ should be normalized by the gird spacing ∆x. It should be noted that the single-node schemes [64][65][66][67] do not involve the present Moment-based scheme.
Inspired by Zhao and Yong (2017) [64], the present single-node second-order BB scheme in the DUGKS can be similarly expressed as, where ρ 0 represents the mean density. Inspired by   [65], the present single-node second-order NEBB scheme in the DUGKS can be similarly expressed as, f As shown in Equation (55), it is known that the distribution function at x ff streams directly from x f . The distribution function at x w is divided into two parts of equilibrium and non-equilibrium, which are represented by the superscripts eq and neq, respectively.
For the equilibrium part at x w , it can be determined by the known velocity and the approximate fluid density, according to the Maxwellian equilibrium distribution function (Equation (4) It has been demonstrated that for low-speed flow, using ρ f (t) and ρ f (t + h) to approximate ρ f (t + h) and ρ w (t + h) have second-and third-order accuracies, respectively [70].
For the non-equilibrium part, it can be approximated and calculated by the nonequilibrium distribution function at x f , with the idea of non-equilibrium bounce back [38,71].
is obtained having at least first-order accuracy, which is enough for deriving a second-order construction of f (x w , ξ α , t n+1/2 ) [65]. Hence, the present singlenode second-order NEBB scheme is second-order accurate in space theoretically. Inspired by Chen et al. (2021) [67], taking the boundary configuration in Figure 4 as an example, a general single-node second-order scheme for the DUGKS is expressed as The constraint 1 < γ ≤ 2q is applied to ensure that both 1 + γ − 2q and 2q − γ are nonnegative.
In future work, the present schemes for curved boundaries in the DUGKS will be further validated and analyzed by theoretical analysis and numerical tests.

Numerical Tests
We perform the numerical tests of the Couette flow, the Poiseuille flow, the Liddriven cavity flow, and the Rayleigh-Taylor instability to further assess the proposed scheme of boundary conditions for the DUGKS. In our simulations, the CFL number is set to be 0.95, and Re = HU0/ν.
The convergence criterion for attaining the steady-state solution is 1000 6 , , Δ represents the velocity in the fluid domain.
To assess the accuracy, we measure the L2 errors of steady velocity fields, The constraint 1 < γ ≤ 2q is applied to ensure that both 1 + γ − 2q and 2q − γ are non-negative.
In future work, the present schemes for curved boundaries in the DUGKS will be further validated and analyzed by theoretical analysis and numerical tests.

Numerical Tests
We perform the numerical tests of the Couette flow, the Poiseuille flow, the Lid-driven cavity flow, and the Rayleigh-Taylor instability to further assess the proposed scheme of boundary conditions for the DUGKS. In our simulations, the CFL number is set to be 0.95, and Re = HU 0 /ν. The convergence criterion for attaining the steady-state solution is where u n ij = u(x i , y j , n∆t) represents the velocity in the fluid domain. To assess the accuracy, we measure the L 2 errors of steady velocity fields, where u is the analytical solution.

The Couette Flow
In the Couette flow, the top wall moves with the horizontal velocity U 0 = 0.1, and the bottom wall is fixed. We apply the proposed scheme of boundary conditions to the top and bottom walls. The inlet and outlet adopt the periodic boundary condition. The gap between the top and bottom wall is set as H = 1. For the no-slip condition, the analytical solution of horizontal velocity in the Couette flow can be written as, It is noted that the convergence criterion for simulating the Couette flow follows Equation (60). To test the accuracy, the L 2 errors of horizontal velocity along the vertical center line are measured at Re = 100, 1000, and 10,000. The results are shown in Tables 1-3    As shown in Figures 5-7, the present results agree very well with the analytic solution. As shown in Tables 1-3, the numerical errors are almost negligible, even with a coarse mesh (N = 8). It is shown that the BB scheme, NEBB scheme, and Moment-based scheme can accurately simulate the Couette flow with the no-slip condition. As shown in Tables 1-3, the L 2 errors in the Moment-based scheme are equal to those in the NEBB scheme, which are a little less than those in the BB scheme at Re = 100.     To test the stability of the proposed scheme of boundary conditions at high Re, some simulations are performed at Re = 10 5 , 10 6 , and 10 7 with N = 16. The steps of reaching the steady state are shown in Table 4. It is found that both the NEBB scheme and thMomentbased scheme converge to the steady state faster than the BB scheme. As shown in Table 5, the L 2 error of the Moment-based scheme is equal to that of the NEBB scheme, and the L 2 errors of the Moment-based scheme and the NEBB scheme are a little less than those in the BB scheme. The L 2 errors are approximately 0.23%, 2.25%, and 18.6% at Re = 10 5 , 10 6 , and 10 7 , respectively. It is shown that the schemes can predict acceptable results for the simulation of the Couette flow at high Re.

The Lid-Driven Cavity Flow
In the lid-driven cavity flow, the top wall moves with the horizontal velocity U 0 = 0.1. The bottom wall and left-and right-side walls are fixed. The walls adopt the present boundary conditions. The cavity length L is set to be 1. It is noted that the convergence criterion for simulating the lid-driven flow follows Equation (60).
We perform some numerical simulations using the Moment-based scheme at Re = 400, 1000, 3200, 5000, 7500, and 10,000 with the different uniform meshes (N × N, N = 32, 64, 128, 256, and 512). The results of velocity profiles are present in Figures 8 and 9. It is found that the results of the Moment-based scheme are in good agreement with the reference data [72] with a relatively fine mesh (N = 128, 256, 512).
To compare the Moment-based scheme with the BB and NEBB scheme, some simulations are performed at Re = 100, 400, 1000, 3200, 5000, 7500, and 10,000 with the fine mesh (N = 256). The results of velocity profiles are present in Figures 10 and 11. It is shown that the results of the Moment-based, BB, and NEBB schemes are all in good agreement with the reference data [72]. Based on a comparison with the reference data, the proposed schemes of the boundary conditions for the DUGKS are valid and suitable for simulating the lid-driven cavity flow.

The Lid-Driven Cavity Flow
In the lid-driven cavity flow, the top wall moves with the horizontal velocity U0 = 0.1. The bottom wall and left-and right-side walls are fixed. The walls adopt the present boundary conditions. The cavity length L is set to be 1. It is noted that the convergence criterion for simulating the lid-driven flow follows Equation (60).
We perform some numerical simulations using the Moment-based scheme at Re = 400, 1000, 3200, 5000, 7500, and 10,000 with the different uniform meshes (N × N, N = 32, 64, 128, 256, and 512). The results of velocity profiles are present in Figures 8 and 9. It is found that the results of the Moment-based scheme are in good agreement with the reference data [72] with a relatively fine mesh (N = 128, 256, 512).
To compare the Moment-based scheme with the BB and NEBB scheme, some simulations are performed at Re = 100, 400, 1000, 3200, 5000, 7500, and 10,000 with the fine mesh (N = 256). The results of velocity profiles are present in Figures 10 and 11. It is shown that the results of the Moment-based, BB, and NEBB schemes are all in good agreement with the reference data [72]. Based on a comparison with the reference data, the proposed schemes of the boundary conditions for the DUGKS are valid and suitable for simulating the lid-driven cavity flow.

The Poiseuille Flow
The Poiseuille flow is driven by an external force ρa x with periodic boundary conditions at the entrance and exit. Both the top and bottom walls adopt the proposed scheme of boundary conditions. The gap between the top and bottom wall is set as H = 1. For the no-slip condition, the analytical solution of horizontal velocity in the Poiseuille flow can be expressed as, u(y)/Umax = 4y/H(1 -y/H), where Umax = a x H 2 /8ν. It is noted that a y = 0 and the convergence criterion for simulating the Poiseuille flow follows Equation (60). In this subsection, it is noted that the "Original" scheme represents the boundary condition with a x = a y = 0 in Equation (37) or Equation (39), and the "Present" scheme represents the proposed scheme of boundary condition with a x = g, a y = 0 in Equation (37) or Equation (39).
To compare the present NEBB scheme (such as Equation (39) with a x = g, a y = 0) with the original NEBB scheme (such as Equation (39) with a x = a y = 0), the original NEBB scheme is also applied to both the top and bottom walls. The results are shown in Figures 12-14 and Table 6.
As shown in Figures 12-14, it seems that both the present and original NEBB schemes agree very well with the analytical solution with N = 32, 64, and 128 at Re = 100, 1000, and 10,000. Furthermore, the L 2 errors of horizontal velocity are measured. As shown in Table 6, the L 2 error of the present NEBB scheme is less than the original NEBB scheme, which shows the present NEBB scheme is more accurate than the original NEBB scheme. It is found that the present and original NEBB schemes are almost second-order accurate.      resents the boundary condition with ax = ay = 0 in Equation (37) or Equation (39), and the "Present" scheme represents the proposed scheme of boundary condition with ax = g, ay = 0 in Equation (37) or Equation (39). To compare the present NEBB scheme (such as Equation (39) with ax = g, ay = 0) with the original NEBB scheme (such as Equation (39) with ax = ay = 0), the original NEBB scheme is also applied to both the top and bottom walls. The results are shown in Figures 12-14 and Table 6.     To compare the present Moment-based scheme (such as Equation (37) with a x = g, a y = 0) with the original Moment-based scheme (such as Equation (37) with a x = a y = 0), the original Moment-based scheme is also applied to both top and bottom walls. The results are shown in   Table 7.  As shown in Figures 12-14, it seems that both the present and original NEBB schemes agree very well with the analytical solution with N = 32, 64, and 128 at Re = 100, 1000, and 10,000. Furthermore, the L2 errors of horizontal velocity are measured. As shown in Table 6, the L2 error of the present NEBB scheme is less than the original NEBB scheme, which shows the present NEBB scheme is more accurate than the original NEBB scheme. It is found that the present and original NEBB schemes are almost second-order accurate.
To compare the present Moment-based scheme (such as Equation (37) with ax = g, ay = 0) with the original Moment-based scheme (such as Equation (37) with ax = ay = 0), the original Moment-based scheme is also applied to both top and bottom walls. The results are shown in Figures 15-17 and Table 7. As shown in Figures 15-17, it seems that both the present and original Momentbased schemes agree very well with the analytical solution with N = 32, 64, and 128 at Re = 100, 1000, and 10,000. Furthermore, the L 2 errors of horizontal velocity are measured. As shown in Table 7, the L 2 error of the present Moment-based scheme is less than the original Moment-based scheme, which shows the present Moment-based scheme is more accurate than the original Moment-based scheme. It is found that the present and original Moment-based schemes are almost second-order accurate.
Comparing Table 6 with Table 7, it is found that the original Moment-based scheme is more accurate than the original NEBB scheme under different meshes, and the present Moment-based scheme is more accurate than the present NEBB scheme with N = 16 and 32, in contrast to the cases with N = 64 and 128.
To compare the present BB and original and present NEBB and Moment-based schemes with the BB and NEBB schemes proposed by Yang et al. [39], we show the L 2 errors in Table 8.   As shown in Figures 15-17, it seems that both the present and original Momentbased schemes agree very well with the analytical solution with N = 32, 64, and 128 at Re = 100, 1000, and 10,000. Furthermore, the L2 errors of horizontal velocity are measured. As shown in Table 7, the L2 error of the present Moment-based scheme is less than the original Moment-based scheme, which shows the present Moment-based scheme is As shown in Table 8, the original and present schemes can predict more accurate results with a finer mesh. The L 2 errors of the present BB and original and present NEBB schemes are not very sensitive to the meshes, but the L 2 errors of the original and present Momentbased schemes are sensitive to the meshes. The L 2 errors of the present BB and original and present Moment-based schemes are not very sensitive to the kinematic viscosity, but the L 2 errors of the original and present NEBB schemes are very sensitive to the kinematic viscosity. The L 2 errors of the original NEBB scheme are more than those of the present NEBB scheme, which are more than those of the NEBB scheme in Ref. [39]. With a fixed kinematic viscosity and grid number, the L 2 errors of the present Moment-based scheme are minimal.

Normal Dipole-Wall Collision
Two counter-rotating vortices are propelled towards a solid boundary, and they collide with the no-slip boundary. The vortex dipole-wall collision is found in nature, such as the effect of the ground on the formulation of secondary vortices when an airplane takes off or lands [73]; another natural phenomenon is the formulation of large-scale vortices in geophysical turbulence on the coasts of seas and oceans [74]. Therefore, the vortex dipole-wall collision is an important problem.  Comparing Table 6 with Table 7, it is found that the original Moment-based scheme is more accurate than the original NEBB scheme under different meshes, and the present Moment-based scheme is more accurate than the present NEBB scheme with N = 16 and 32, in contrast to the cases with N = 64 and 128.
To compare the present BB and original and present NEBB and Moment-based schemes with the BB and NEBB schemes proposed by Yang et al. [39], we show the L2 errors in Table 8.  In this study, two counter-rotating vortices are confined to a square box with a size of [−1, 1] × [−1, 1]. The present boundary conditions are applied to the no-slip walls. The initial vortex is located at positions (x 1 , y 1 ) and (x 2 , y 2 ). The initial velocities are u x0 = −0.5|We|(y − y 1 ) exp −(r 1 /r 0 ) 2 + 0.5|We|(y − y 2 ) exp −(r 2 /r 0 ) 2 , The symbol r 0 denotes the radius of the monopoles, and We represents the strength of the vortices. To compare the numerical results, the total kinetic energy E(t) and the total enstrophy Ω(t) are calculated where ω is the vorticity.
Clercx and Bruneau [75] simulated a dipole-wall collision using a Finite Difference Method (FDM) and the Chebyshev Pseudospectral Method (CPM). Mohammed et al. [40] simulated the dipole-wall collision using the LBM with two relaxation time models (TRT-LBM). Their authoritative data will be used as benchmark numerical results. For a direct comparison between the present work and the work of Ref. [40] and Ref. [75], we use the As shown in Table 8, the original and present schemes can predict more accurate results with a finer mesh. The L2 errors of the present BB and original and present NEBB schemes are not very sensitive to the meshes, but the L2 errors of the original and present Moment-based schemes are sensitive to the meshes. The L2 errors of the present BB and original and present Moment-based schemes are not very sensitive to the kinematic viscosity, but the L2 errors of the original and present NEBB schemes are very sensitive to the kinematic viscosity. The L2 errors of the original NEBB scheme are more than those of the present NEBB scheme, which are more than those of the NEBB scheme in Ref. [39]. With a fixed kinematic viscosity and grid number, the L2 errors of the present Momentbased scheme are minimal.

Normal Dipole-Wall Collision
Two counter-rotating vortices are propelled towards a solid boundary, and they collide with the no-slip boundary. The vortex dipole-wall collision is found in nature, such as the effect of the ground on the formulation of secondary vortices when an airplane takes off or lands [73]; another natural phenomenon is the formulation of large-scale vortices in geophysical turbulence on the coasts of seas and oceans [74]. Therefore, the vortex dipole-wall collision is an important problem.
In this study, two counter-rotating vortices are confined to a square box with a size of [−1, 1] × [−1, 1]. The present boundary conditions are applied to the no-slip walls. The initial vortex is located at positions (x1, y1) and (x2, y2). The initial velocities are In the dipole-wall collision benchmark test for the normal case, two monopoles are located at positions (x 1 , y 1 ) = (0, 0.1) and (x 2 , y 2 ) = (0, −0.1) initially. Then they are propelled towards the right wall. To test the effect of the present schemes on the vortices after the dipole collides, vorticity contour plots are present in Figures 18-21. As shown in vorticity contour plots, the present schemes are effective to simulate the dipole-wall collision, but the results with BB, NEBB, and Moment schemes are almost indistinguishable. To analyze the present schemes quantitatively, the values of the first and second maxima enstrophy of dipoles will be grouped for comparisons.
The present results for the first and second local maxima of the enstrophy are shown in Tables 9 and 10, and we compare them with the results in Refs. [40,75]. The results computed using the moment-based boundary condition with DUGKS and TRT-LBM are in good agreement. However, the data obtained by the DUGKS appear to be more accurate than the data obtained by the TRT-LBM in the sense that the data are in closer agreement with the data obtained by FDM and CPM. The data obtained using the present momentbased scheme appear to be more accurate than the data obtained using the present BB and NEBB schemes in the sense that they are in closer agreement with the data obtained by FDM and CPM. These show that the proposed moment-based scheme can be a competitive method and gives us the confidence to use it to impose physically more complex conditions. ter the dipole collides, vorticity contour plots are present in Figures 18-21. As shown in vorticity contour plots, the present schemes are effective to simulate the dipole-wall collision, but the results with BB, NEBB, and Moment schemes are almost indistinguishable. To analyze the present schemes quantitatively, the values of the first and second maxima enstrophy of dipoles will be grouped for comparisons.     The present results for the first and second local maxima of the enstrophy are shown in Tables 9 and 10, and we compare them with the results in Refs. [40,75]. The results computed using the moment-based boundary condition with DUGKS and TRT-LBM are in good agreement. However, the data obtained by the DUGKS appear to be more accurate than the data obtained by the TRT-LBM in the sense that the data are in closer agreement with the data obtained by FDM and CPM. The data obtained using the present moment-based scheme appear to be more accurate than the data obtained using   For further comparison with reference data [76], we present the values of the energy E(t) and enstrophy Ω(t) at different times, as shown in Tables 11 and 12. The results obtained by the D3Q19-CM-LBM [76] show a slight mismatch (up to 3%) with respect to the LBM study by Mohammed et al. [40]. It should be noted that the results obtained by the DUGKS with the present Moment scheme are closer to the reference ones by Clercx and Bruneau [75] than those obtained by the TRT-LBM with the Moment scheme [40]. We address this behavior regarding the adoption in the present work of a more accurate boundary condition.

Rayleigh-Taylor Instability
Rayleigh-Taylor instability can occur when a layer of heavy fluid descends as light fluid below it rises. We perform a simulation using the same parameters as the first case of Re = 256, as shown in Figure 7.4 in Ref. [77]. The left and right boundaries adopt the periodic boundary conditions. For the upper and lower boundaries, the mentioned boundary conditions are applied. A computational domain X × Y = 128 × 512 is employed.
Initially, there is a zero-velocity field, and the location of the perturbed interface is set as y = 0.5Y + 0.1Xcos(2πx/X), where x, y, X, and Y are all in lattice units. We set the densities of heavy fluid and light fluid to be ρ h = 0.12 and ρ l = 0.04, respectively, so that At = (ρ h − ρ l )/(ρ h + ρ l ) = 0.5. We set U = 0.04, and the gravitational acceleration is g = U 2 /X = 1.25 × 10 −5 in the -y direction.
In this study, the HCZ model [52] is introduced for the DUGKS to simulate the Rayleigh-Taylor instability. In the model, two distribution functions satisfy the following equations: The distribution functions g and f, corresponding to the density ρ and phase order φ, denote the particle distribution function in terms of position x, discrete particle velocity ξ α , and time t. The relaxation times τ g and τ f are related to the viscosity and the mobility coefficient in the Cahn-Hilliard equation. Usually, we set τ g = τ f .
S g and S f represent the source terms, which can be written as F s represents the force associated with surface tension, where the parameter κ determines the strength of surface tension and κ = 0.01 in our simulations.
ψ represents the function of the density ρ or phase order φ. ψ(ρ) and ψ(φ) are determined as [78] ψ(ρ) = p − c 2 s ρ (72) In our simulations, a = 12RT, b = 4. The equilibrium distributions g eq and f eq can be calculated as Since Equations (66) and (67) share the same format as Equation (1), we will not repeat the detailed depiction of evolution in the DUGKS. In this subsection, we study the effect of varied boundary conditions on the numerical simulation results of Rayleigh-Taylor instability.
The In this subsection, it is noted that the "Original" scheme represents the boundary condition with a x = a y = 0 in Equation (37) or Equation (39), and the "Present" scheme represents the proposed scheme of the boundary condition with a x = 0, a y = -g in Equation (37) or Equation (39).
As shown in Figures 22-27, the fluid interface does not diffuse in our simulations by applying the mentioned boundary conditions. Compared to Figure 7.4 in Ref. [77], the evolution of the interface is indistinguishable from the results computed on the same mesh. It is shown that the present schemes are accurate and stable to simulate the Rayleigh-Taylor instability. Figures 28 and 29 show the density and vertical velocity profiles across the spike, respectively. As shown in Figures 28 and 29, the interface thickness takes approximately four grid spacings (y = 31, 32, 33, and 34). There exist some "jiggles" near the interface, and the "jiggles" are different due to using different boundary conditions. From the magnified view in Figures 28 and 29, the density and absolute value of vertical velocity near the interface using present NEBB and Moment-based schemes are less than that using BB, original NEBB, and original Moment-based schemes. Although the interface thickness is not affected by different boundary conditions, choosing a wall boundary condition can influence the density and velocity near the interface. Usually, numerical instability occurs near the interface. So, the finding reminds us that wall boundary conditions should be treated with care.
The simulation is executed on a Win10 X64 system with Intel(R) Core(TM) i5-8250U CPU (1.60 GHz)s. For comparison, the convergence error and CPU time at the last time step are recorded in Table 13. As shown in Table 13, the error and CPU time of the present Moment-based scheme are minimal, although they are almost indistinguishable. sents the proposed scheme of the boundary condition with ax = 0, ay = -g in Equation (37) or Equation (39).   sents the proposed scheme of the boundary condition with ax = 0, ay = -g in Equation (37) or Equation (39).                near the interface using present NEBB and Moment-based schemes are less than that using BB, original NEBB, and original Moment-based schemes. Although the interface thickness is not affected by different boundary conditions, choosing a wall boundary condition can influence the density and velocity near the interface. Usually, numerical instability occurs near the interface. So, the finding reminds us that wall boundary conditions should be treated with care.   The simulation is executed on a Win10 X64 system with Intel(R) Core(TM) i5-8250U CPU (1.60 GHz)s. For comparison, the convergence error and CPU time at the last time step are recorded in Table 13. As shown in Table 13, the error and CPU time of the present Moment-based scheme are minimal, although they are almost indistinguishable.

Conclusions
Owing to the DUGKS with transformed distribution functions and force terms, we should convert the original boundary conditions into a new format. In this study, the proposed schemes of the BB, NEBB, and Moment-based boundary conditions are proposed for the DUGKS. The boundary conditions will be translated into constraints on the unknown transformed distribution functions at a half time step (t + 0.5Δt). The present work tests and analyzes, for the first time, the Moment-based scheme for the

Conclusions
Owing to the DUGKS with transformed distribution functions and force terms, we should convert the original boundary conditions into a new format. In this study, the proposed schemes of the BB, NEBB, and Moment-based boundary conditions are proposed for the DUGKS. The boundary conditions will be translated into constraints on the unknown transformed distribution functions at a half time step (t + 0.5∆t). The present work tests and analyzes, for the first time, the Moment-based scheme for the DUGKS. The mentioned boundary conditions are evaluated through theoretical analysis and numerical tests.
Using the steady unidirectional flow, we theoretically analyze the numerical slip error of the present BB, NEBB, and Moment-based schemes. It is found that the numerical slip errors of both NEBB and Moment-based schemes are equal to zero, which implements the no-slip condition at the wall boundary.
The accuracy and stability of the present schemes are validated by numerical tests of Couette flow, Poiseuille flow, Lid-driven cavity flow, normal dipole-wall collision, and Rayleigh-Taylor instability. The following conclusions are obtained: (1) Couette flow The present schemes can predict accurate results of simulating the Couette flow, even with a coarse mesh for a large Reynolds number (Re = 10 6 ). The L 2 errors in both NEBB and Moment-based schemes are equal, which are a little less than those in the BB scheme at Re = 100 under different meshes. Both the NEBB and Moment-based schemes are more accurate and converge to the steady state faster than the BB scheme at high Re (Re = 100,000, 1,000,000, and 10,000,000).
(2) Lid-driven cavity flow The present schemes can predict accurate results by simulating the lid-driven cavity flow. It is found that the results of the BB scheme are in better agreement with the reference data than those of Moment and NEBB schemes at Re = 10,000.
(3) Poiseuille flow The results of the present schemes agree very well with the analytical solution under different meshes (N = 32, 64, 128) at Re = 100, 1000, and 10,000. The L 2 errors of the present NEBB and Moment-based schemes are less than the original NEBB and Moment-based schemes, respectively, which shows the present schemes are more accurate than the original schemes. It is found that the present schemes for the DUGKS are second-order accurate. The original Moment scheme is more accurate than the original NEBB scheme under different meshes. The present Moment scheme is more accurate than the present NEBB scheme with N = 16 and 32, in contrast to the cases with N = 64 and 128. Compared with the BB and NEBB schemes in Ref. [39], the L 2 errors of the present Moment-based scheme are minimal.

(4) Normal dipole-wall collision
With the present moment-based scheme, the DUGKS appears to be more accurate than the TRT-LBM. The data obtained by using the present moment-based scheme appears to be more accurate than the data obtained by using the present BB and NEBB schemes in the sense that they are in closer agreement with the benchmark data. These show that the proposed moment-based scheme can be a competitive method.

(5) Rayleigh-Taylor instability
The present results agree well with the reference results, which show the present schemes can effectively capture the evolution of the interface. It is found that the error and CPU time of the present Moment-based scheme are minimal among the mentioned boundary conditions. The Rayleigh-Taylor instability simulation shows that choosing a wall boundary condition can influence the density and velocity near the interface slightly. Institutional Review Board Statement: Not applicable.

Data Availability Statement:
The data that support the findings of this study are available within the article.

Conflicts of Interest:
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.