Mass-Conserved Wall Treatment of the Non-Equilibrium Extrapolation Boundary Condition in Lattice Boltzmann Method

Abstract

In lattice Boltzmann simulations, the widely used non-equilibrium extrapolation method for velocity and pressure boundary conditions can cause a constant mass leakage under certain circumstances, particularly when an external force field is imposed on the fluid domain. The non-equilibrium distribution function at the boundary uses a first-order extrapolation method on the corresponding data of adjacent fluid nodes. In addition, based on this extrapolation method, the macroscopic velocity and density at the boundary nodes are obtained. Therefore, the corresponding equilibrium component of the distribution function can be calculated explicitly. Regarding the no-slip wall boundary condition, we found that the mass leakage primarily results from the extrapolation scheme for the density term in the equilibrium component of the distribution function at the boundary node. In this study, a mass-conserved wall treatment method is developed to correct the existing density term for guaranteeing the conservation of mass. Several benchmark test cases were simulated and compared to prove the justification of the newly developed mass-conserved boundary condition, and the results show a good agreement with those in the existing literature.

1. Introduction

In recent decades, the lattice Boltzmann method (LBM) [1,2] has been developed into a promising and successful alternative to the conventional computational fluid dynamics (CFD) method, particularly for incompressible flows [3], porous-media flows [4], magnetohydrodynamics [5], and single- and multi-phase fluid flows [6]; more recently see [7,8,9,10,11]. Different from the traditional CFD, the LBM models the fluid as a group of fictive particles by using discrete distribution functions [1]. These particles perform consecutive collision and streaming processes over a discrete lattice mesh. Owing to its particulate nature and local dynamics, the LBM has several advantages over the traditional CFD methods, particularly the intrinsic parallelism of the algorithm and data structure [12,13], and the treatment of complex boundaries [14,15,16]. Moreover, the lattice Boltzmann equation would have a second-order accuracy in the time and space discretisation [2], which is comparable to the accuracy of the finite-volume numerical solutions of the Navier-Stokes equations.

In the LBM, the definition of proper boundary conditions has always been an important issue to simulate a fluid flow accurately and stably in the domain. Different from the conventional CFD methods, the boundary condition of the LBM is implemented in the form of distribution functions. Therefore, an appropriate formulation needs to be proposed to obtain the unknown distribution functions at the boundary node based on other known distribution functions at the boundary or fluid nodes. The most commonly used boundary condition in the LBM is the bounce-back method [17], which was first proposed for lattice gas automata [18]. In this boundary condition, the particle distribution scatters back to the previous fluid node when it streams to a solid wall node, which can be directly formulated and implemented on the boundary of the calculation domain. However, regrettably, the bounce-back boundary condition only has a first-order accuracy [19], which is not consistent with the order of the LBM in interior points. To improve its inherent accuracy, many other boundary conditions have been developed in recent decades (Noble et al. [20], Inamuro et al. [21], Latt and Chopard [22], such as the halfway bounce-back boundary method [23], non-equilibrium bounce-back scheme [19], and extrapolation scheme [24].

In the extrapolation boundary method proposed by Chen et al. [24], the unknown distributions on the boundary are obtained from the distributions in the fluid through a second-order scheme. Although this method is rather straightforward and can be extended to other boundary conditions, its numerical stability was found to be unsatisfactory. Regarding this issue, Guo et al. [25] proposed a new extrapolation scheme to increase the stability and accuracy. In their approach, the distribution function is divided into two different parts: equilibrium and non-equilibrium. In particular, at the boundary nodes, the equilibrium distribution function was calculated from the macroscopic velocity and density field. The non-equilibrium distribution function was obtained by using a first-order extrapolation of the non-equilibrium component of the distributions in the flow domain. This non-equilibrium extrapolation (NEE) boundary method is also proven to be of second order. Moreover, it can be easily extended to design schemes for other boundary conditions, such as mass, momentum, and heat flux. In this study, the NEE boundary condition was analysed and discussed in detail. We found that it cannot meet the requirement of mass conservation, especially when an external force field such as gravitational force is imposed on the flow domain. This mass leakage (ML) can cause serious problems in a variety of circumstances. Therefore, this paper proposes a new implementation scheme for the NEE boundary method to guarantee mass conservation in the whole flow domain.

Although abundant literature has discussed the mass conservation of many other boundary methods, the consequences of ML in the flow domain have rarely been analysed in detail. In this study, we aimed to observe the generation and aftermaths of ML as well as the effect of its inherent problem on the flow domain. Undertaking this issue, this paper is organised as follows: Section 2 outlines the basic governing equation and computational methodology. Section 3 compares and analyses the numerical results. Finally, Section 4 presents the major conclusions.

2. Design of LBM Simulation

2.1. Single-Relaxation-Time LBM

In the single-relaxation-time lattice Boltzmann method (SRT-LBM), the collision model proposed by Bhatnagar, Gross, and Krook [26] is applied to the lattice Boltzmann equation (LBE), which is given by

$$f_i(\mathbf{x}+{\mathbf{c}}_i\Delta t,t+\Delta t)-f_i(\mathbf{x},t)=-\frac{1}{\tau }(f_i(\mathbf{x},t)-f_i^{eq}(\mathbf{x},t))$$

(1)

where $f_i(\mathbf{x},t)$ and $f_i^{eq}(\mathbf{x},t)$ are the particle and equilibrium distribution functions at $(\mathbf{x},t)$, ${\mathbf{c}}_{\mathbf{i}}$ is the particle velocity along the ith direction, and $\tau$ is the relaxation time that manages the rate of the distribution function’s approach to equilibrium state. The implementation of the controlling equation can be divided into two processes, the collision and streaming processes, which are given by

$$\mathrm{Collision} \mathrm{process}:{\tilde{f}}_i(\mathbf{x},t)=f_i(\mathbf{x},t)-\frac{1}{\tau }(f_i(\mathbf{x},t)-f_i^{eq}(\mathbf{x},t))$$

(2)

$$\mathrm{Streaming} \mathrm{process}:f_i(\mathbf{x}+{\mathbf{c}}_i\Delta t,t+\Delta t)={\tilde{f}}_i(\mathbf{x},t)$$

(3)

where ${\tilde{f}}_i(\mathbf{x},t)$ refers to the post-collision distribution function. In the LBM, the relaxation time should be linked to the viscosity of the fluid through the Chapman-Enskog expansion in such a way that the macroscopic variables such as velocity, pressure, and density computed from LBM can satisfy the Navier-Stokes equations with second-order accuracy [27]. It is denoted as

$$\tau =\frac{6\nu +1}{2}$$

(4)

The equilibrium distribution function is a type of polynomial equation simplified from the exponential form of the Maxwellian distribution function (He et al. [28]) as shown below:

$$f_i^{eq}=\rho w_i\left[1+\frac{{\mathbf{c}}_{\mathbf{i}}·\mathbf{u}}{c_s^2}+\frac{{({\mathbf{c}}_{\mathbf{i}}·\mathbf{u})}^2}{2c_s^4}-\frac{\mathbf{u}·\mathbf{u}}{2c_s^2}\right]$$

(5)

where $c_s=c/\sqrt{3}$ is the lattice sound speed and $w_i$ is the lattice weight. The macroscopic quantities such as density, velocity, and pressure can be obtained from the distribution functions as follows:

$$\rho =\sum _{i=0}^{N-1}{f}_i$$

(6)

$$\rho \mathbf{u}=\sum _{i=0}^{N-1}{f}_i{\mathbf{c}}_{\mathbf{i}}$$

(7)

$$p=\rho c_s^2$$

(8)

where N is the total number of discrete velocity directions. In the LBM, there are several lattice models. For 2D flows, the 2D nine-velocity model (D2Q9) is widely used and has nine discrete velocity vectors of ${\mathbf{c}}_{\mathbf{i}}$, as shown in Figure 1.

For the D2Q9 lattice model, the discrete velocity and lattice weights are written as

$${\mathbf{c}}_{\mathbf{i}}=\left\{\begin{array}{cc}0\hfill & i=0,\hfill \\ c(cos(\frac{(i-1)\pi }{4}),sin(\frac{(i-1)\pi }{4})\hfill & i=1,2,3,4,\hfill \\ \sqrt{2}c(cos(\frac{(i-1)\pi }{4}),sin(\frac{(i-1)\pi }{4})\hfill & i=5,6,7,8\hfill \end{array}\right.$$

(9)

$$w_i=\left\{\begin{array}{cc}4/9\hfill & i=0,\hfill \\ 1/9\hfill & i=1,2,3,4,\hfill \\ 1/36\hfill & i=5,6,7,8\hfill \end{array}\right.$$

(10)

As the above Equation (1) does not consider the force term, this would cause a variety of fluid problems where an external or internal force is involved, such as in a pressure-driven or multi-phase flows. In this study, the lumped-forcing LBE [28,29] is used, which is expressed as

$$f_i(\mathbf{x}+{\mathbf{c}}_i\Delta t,t+\Delta t)-f_i(\mathbf{x},t)=-\frac{1}{\tau }(f_i(\mathbf{x},t)-f_i^{eq}(\mathbf{x},t))+\Delta t{F}_i(\mathbf{x},t)$$

(11)

where $F_i$ denotes the discrete external force, which can be obtained from

$$F_i(\mathbf{x},t)=\frac{w_i}{c_s^2}{\mathbf{c}}_{\mathbf{i}}·\mathbf{F}(\mathbf{x},t)$$

(12)

where $\mathbf{F}$ denotes the external force vector. The external force effect can be achieved by updating the post-collision distribution function, which is given as

$${\tilde{f}}_i(\mathbf{x},t)=f_i(\mathbf{x},t)-\frac{1}{\tau }(f_i(\mathbf{x},t)-f_i^{eq}(\mathbf{x},t))+\frac{w_i}{c_s^2}{\mathbf{c}}_{\mathbf{i}}·\mathbf{F}(\mathbf{x},t)\Delta t$$

(13)

2.2. Existing NEE Boundary Condition

In the NEE boundary method, the distribution function can be decomposed into its equilibrium and non-equilibrium components, as shown below:

$$f_i(\mathbf{x},t)=f_i^{eq}(\mathbf{x},t)+f_i^{neq}(\mathbf{x},t))$$

(14)

where $f_i^{eq}(\mathbf{x},t)$ and $f_i^{neq}(\mathbf{x},t))$ are the equilibrium and non-equilibrium components of $f_i(\mathbf{x},t)$, respectively. Many boundary conditions in the LBM are implemented based on post-collision distribution functions. Therefore, based on this decomposition scheme, the post-collision distribution function can be rewritten as

$${\tilde{f}}_i(\mathbf{x},t)=f_i^{eq}(\mathbf{x},t)+(1-\omega )f_i^{neq}(\mathbf{x},t))$$

(15)

where $\omega =1/\tau$ is the relaxation frequency. To obtain a complete post-collision distribution function ${\tilde{f}}_i({\mathbf{x}}_{\mathbf{b}},t)$ at the boundary node, both $f_i^{eq}(\mathbf{x},t)$ and $f_i^{neq}(\mathbf{x},t))$ at the boundary node need to be known functions in advance. As shown in Figure 2, the fluid node ${\mathbf{x}}_{\mathbf{f}\mathbf{1}}$ is the closest adjacent to the boundary node ${\mathbf{x}}_{\mathbf{b}}$, whereas $x_{f2}$ and $x_{f3}$ are the neighbouring nodes in the fluid.

Based on the Chapman-Enskog analysis, the non-equilibrium distribution function at the boundary node is approximated by a first-order extrapolation, as follows:

$$f_i^{neq}({\mathbf{x}}_{\mathbf{b}},t)=f_i^{neq}({\mathbf{x}}_{\mathbf{f}\mathbf{1}},t)=f_i({\mathbf{x}}_{\mathbf{f}\mathbf{1}},t)-f_i^{eq}({\mathbf{x}}_{\mathbf{f}\mathbf{1}},t)$$

(16)

After completing the streaming process, the density and velocity, along with the distribution function $f_i({\mathbf{x}}_{\mathbf{f}\mathbf{1}},t)$, are calculated at the fluid node ${\mathbf{x}}_{\mathbf{f}}$. Based on Equation (5), the equilibrium distribution function can be retrieved from the known density and velocity. Finally, the non-equilibrium distribution function can be obtained at the boundary node.

The equilibrium distribution function can be determined according to the velocity and pressure boundary conditions. For the velocity boundary condition, the velocity can be obtained, but the pressure is still unknown. From Equation (8), it is found that the pressure can be directly calculated from the density. The density at the boundary node can be extrapolated from the adjacent fluid node. Therefore, the equilibrium distribution function at the boundary node can be given as

$$f_i^{eq}({\mathbf{x}}_{\mathbf{b}},t)=\rho ({\mathbf{x}}_{\mathbf{f}\mathbf{1}},t)w_i\left[1+\frac{{\mathbf{c}}_{\mathbf{i}}·\mathbf{u}({\mathbf{x}}_{\mathbf{b}},t)}{c_s^2}+\frac{{({\mathbf{c}}_{\mathbf{i}}·\mathbf{u}({\mathbf{x}}_{\mathbf{b}},t))}^2}{2c_s^4}-\frac{\mathbf{u}({\mathbf{x}}_{\mathbf{b}},t)·\mathbf{u}({\mathbf{x}}_{\mathbf{b}},t)}{2c_s^2}\right]$$

(17)

Regarding the pressure boundary condition, the density is known, but the velocity is unknown at the boundary node. Similarly to the velocity boundary condition, the velocity can be extrapolated from the adjacent fluid node. Therefore, the equilibrium distribution function at the boundary node can be expressed as

$$f_i^{eq}({\mathbf{x}}_{\mathbf{b}},t)=\rho ({\mathbf{x}}_{\mathbf{b}},t)w_i\left[1+\frac{{\mathbf{c}}_{\mathbf{i}}·\mathbf{u}({\mathbf{x}}_{\mathbf{f}\mathbf{1}},t)}{c_s^2}+\frac{{({\mathbf{c}}_{\mathbf{i}}·\mathbf{u}({\mathbf{x}}_{\mathbf{f}\mathbf{1}},t))}^2}{2c_s^4}-\frac{\mathbf{u}({\mathbf{x}}_{\mathbf{f}\mathbf{1}},t)·\mathbf{u}({\mathbf{x}}_{\mathbf{f}\mathbf{1}},t)}{2c_s^2}\right]$$

(18)

After obtaining both the equilibrium and non-equilibrium distribution functions at the boundary node, the post-collision distribution function can be calculated based on Equation (15).

2.3. Source of ML in the NEE Boundary Condition

The post-collision distribution functions can be obtained at both the boundary and fluid nodes after implementation of the collision process and boundary condition at the boundary nodes. Then, the post-collision distribution functions at the boundary nodes will propagate into the fluid domain, but they move towards the boundary nodes at the neighbouring fluid nodes. In this study, as we only focus on the discussion of mass conservation on fluid domain, ML is defined as the exchange surplus of post-collision distribution functions between the neighbouring fluid and boundary nodes. The internal fluid will move out when ML has a positive value; however, the external fluid will move into the computational domain with a negative value. As shown in Figure 2, on a bottom boundary, it can be easily found that the number of post-collision distribution functions used for calculating mass leakage would vary on different positions. On the corner and boundary nodes adjacent to the corner, 2 and 4 distribution functions would be needed respectively, while 6 distribution functions for those boundary nodes which are far away from the corner. Because the same scheme has been implemented for calculating the mass leakage and implementing the newly developed mass-conserved boundary method, the detailed derivation process is only described here for those boundary nodes far away from the corner. Based on above definition, the mass leakage on a single boundary node can be given as

$$\begin{array}{ccc}\hfill M_l({\mathbf{x}}_{\mathbf{b}},t)& =& \left[{\tilde{f}}_4({\mathbf{x}}_{\mathbf{f}\mathbf{1}},t)+{\tilde{f}}_7({\mathbf{x}}_{\mathbf{f}\mathbf{2}},t)+{\tilde{f}}_8({\mathbf{x}}_{\mathbf{f}\mathbf{3}},t)\right]\hfill \\ & & -\left[{\tilde{f}}_2({\mathbf{x}}_{\mathbf{b}},t)+{\tilde{f}}_5({\mathbf{x}}_{\mathbf{b}},t)+{\tilde{f}}_6({\mathbf{x}}_{\mathbf{b}},t)\right]\hfill \end{array}$$

(19)

where $M_l({\mathbf{x}}_{\mathbf{b}},t)$ refers to ML on the boundary node ${\mathbf{x}}_{\mathbf{b}}$, and the three fluid nodes adjacent to the boundary node ${\mathbf{x}}_{\mathbf{b}}$ are denoted by ${\mathbf{x}}_{\mathbf{f}\mathbf{1}}$, ${\mathbf{x}}_{\mathbf{f}\mathbf{2}}$, and ${\mathbf{x}}_{\mathbf{f}\mathbf{3}}$, respectively. Based on the definition of the post-collision distribution function in Equation (2), it can also be calculated from the equilibrium and non-equilibrium distribution functions, as follows:

$${\tilde{f}}_i=f_i-\frac{1}{\tau }\left(f_i-f_i^{eq}\right)=f_i^{eq}+{\omega }^{\prime }{f}_i^{neq}$$

(20)

where ${\omega }^{\prime }=(1-1/\tau )$. Based on the definition of equilibrium distribution function above, it can be denoted as the multiplication of the density and velocity-relevant terms, as shown below:

$$f_i^{eq}=\rho Z_i(\mathbf{u}),Z_i(\mathbf{u})=w_i\left[1+\frac{{\mathbf{c}}_{\mathbf{i}}·\mathbf{u}}{c_s^2}+\frac{{({\mathbf{c}}_{\mathbf{i}}·\mathbf{u})}^2}{2c_s^4}-\frac{\mathbf{u}·\mathbf{u}}{2c_s^2}\right]$$

(21)

As shown in Figure 2, when the gravitational force field exists in the flow domain, the post-collision distribution function of fluid nodes should be calculated by using Equation (13). In this analysis of ML, the relaxation time was set to $\tau =1$, and the post-collision distribution function can be denoted by the corresponding equilibrium distribution function based on Equation (20). In the no-slip boundary condition, the velocity on the boundary is known and fixed to zero. Substituting the post-collision distribution functions in Equation (19), the ML can be given as

$$\begin{array}{ccc}\hfill M_l({\mathbf{x}}_{\mathbf{b}},t)& =& [\rho ({\mathbf{x}}_{\mathbf{f}\mathbf{1}})Z_4(\mathbf{u}({\mathbf{x}}_{\mathbf{f}\mathbf{1}}))+\rho ({\mathbf{x}}_{\mathbf{f}\mathbf{2}})Z_7(\mathbf{u}({\mathbf{x}}_{\mathbf{f}\mathbf{2}}))\hfill \\ & & +\rho ({\mathbf{x}}_{\mathbf{f}\mathbf{3}})Z_8(\mathbf{u}({\mathbf{x}}_{\mathbf{f}\mathbf{3}}))+3g(w_4+w_7+w_8)]\hfill \\ & & -\rho ({\mathbf{x}}_{\mathbf{b}})(w_2+w_5+w_6)\hfill \end{array}$$

(22)

As the fluid nodes ${\mathbf{x}}_{\mathbf{f}\mathbf{1}}$, ${\mathbf{x}}_{\mathbf{f}\mathbf{2}}$, and ${\mathbf{x}}_{\mathbf{f}\mathbf{3}}$ are adjacent to the no-slip boundary, in certain cases, their velocity can be very close to each other and approach to zero, which is given as

$$\mathbf{u}({\mathbf{x}}_{\mathbf{f}\mathbf{1}})\approx \mathbf{u}({\mathbf{x}}_{\mathbf{f}\mathbf{2}})\approx \mathbf{u}({\mathbf{x}}_{\mathbf{f}\mathbf{3}})\approx (0,0)$$

(23)

With the approximation in Equation (23), the velocity term $Z_i(\mathbf{u})$ in the equilibrium distribution function can be simplified as $Z_4(\mathbf{u}({\mathbf{x}}_{\mathbf{f}\mathbf{1}}))=w_4$, $Z_7(\mathbf{u}({\mathbf{x}}_{\mathbf{f}\mathbf{2}}))=w_7$, $Z_8(\mathbf{u}({\mathbf{x}}_{\mathbf{f}\mathbf{3}}))=w_8$. In addition, the density on the boundary node $\rho ({\mathbf{x}}_{\mathbf{b}})$ is extrapolated from the adjacent fluid node $\rho ({\mathbf{x}}_{\mathbf{f}\mathbf{1}})$ for the NEE boundary method, and finally, the ML can be expressed as

$$\begin{array}{ccc}\hfill M_l({\mathbf{x}}_{\mathbf{b}},t)& =& \frac{1}{36}\left[\rho ({\mathbf{x}}_{\mathbf{f}\mathbf{2}})+\rho ({\mathbf{x}}_{\mathbf{f}\mathbf{3}})-2\rho ({\mathbf{x}}_{\mathbf{f}\mathbf{1}}))\right]+\frac{g}{2}\hfill \\ & =& \frac{1}{36}\left[\rho ({\mathbf{x}}_{\mathbf{f}\mathbf{2}})+\rho ({\mathbf{x}}_{\mathbf{f}\mathbf{3}})-2\rho ({\mathbf{x}}_{\mathbf{b}}))\right]+\frac{g}{2}\hfill \end{array}$$

(24)

From Equation (24), it is found that the gravitational force effect on ML cannot be negligible when the external force is perpendicular to the boundary. In addition, it will increase with the increasing magnitude of gravitational force. When there is no external force field on the flow domain, if the density on the adjacent fluid nodes (i.e., ${\mathbf{x}}_{\mathbf{f}\mathbf{1}}$, ${\mathbf{x}}_{\mathbf{f}\mathbf{2}}$, and ${\mathbf{x}}_{\mathbf{f}\mathbf{3}}$) varies linearly along the boundary, then there is little ML on the boundary node ${\mathbf{x}}_{\mathbf{b}}$. However, in the whole flow domain, owing to the complicated flow structures, it is not realistic to keep this relationship at every position of the boundary. Moreover, without the approximation in Equation (23), this ML would be stronger.

2.4. Newly Developed Mass-Conserved Boundary Condition

To guarantee mass conservation strictly at the local boundary nodes, the ML must be zero. As discussed above, the ML is relevant to the exchange of post-collision distribution functions between the neighbouring fluid and boundary nodes. After making the collision process, all the post-collision distribution functions at the fluid nodes are known.

On the specific velocity boundary, the values of $Z_i(\mathbf{u})$ at different boundary nodes are all the same. To enforce local mass conservation, the ML on a single boundary node should be zero, as shown below:

$${\tilde{f}}_2({\mathbf{x}}_{\mathbf{b}},t)+{\tilde{f}}_5({\mathbf{x}}_{\mathbf{b}},t)+{\tilde{f}}_6({\mathbf{x}}_{\mathbf{b}},t)={\tilde{f}}_4({\mathbf{x}}_{\mathbf{f}\mathbf{1}},t)+{\tilde{f}}_7({\mathbf{x}}_{\mathbf{f}\mathbf{2}},t)+{\tilde{f}}_8({\mathbf{x}}_{\mathbf{f}\mathbf{3}},t)$$

(25)

where the post-collision distribution functions have been updated by including external force effects. Substituting the post-collision distribution functions at the boundary node with Equation (20), it becomes as follows:

$$\begin{array}{c}\hfill {\tilde{f}}_2^{eq}({\mathbf{x}}_{\mathbf{b}},t)+{\tilde{f}}_5^{eq}({\mathbf{x}}_{\mathbf{b}},t)+{\tilde{f}}_6^{eq}({\mathbf{x}}_{\mathbf{b}},t)=[{\tilde{f}}_4({\mathbf{x}}_{\mathbf{f}\mathbf{1}},t)+{\tilde{f}}_7({\mathbf{x}}_{\mathbf{f}\mathbf{2}},t)+{\tilde{f}}_8({\mathbf{x}}_{\mathbf{f}\mathbf{3}},t)]\\ \hfill -{\omega }^{\prime }\left[{\tilde{f}}_2^{neq}({\mathbf{x}}_{\mathbf{b}},t)+{\tilde{f}}_5^{neq}({\mathbf{x}}_{\mathbf{b}},t)+{\tilde{f}}_6^{neq}({\mathbf{x}}_{\mathbf{b}},t)\right]\end{array}$$

(26)

As with the NEE method above, the non-equilibrium distribution functions at the boundary node ${\mathbf{x}}_{\mathbf{b}}$ are approximated directly with first-order accuracy from the non-equilibrium distribution functions of the adjacent fluid node ${\mathbf{x}}_{\mathbf{f}\mathbf{1}}$. Substituting this approximation and Equation (21) into Equation (26), the subsequent equation can be obtained as follows:

$$\begin{array}{c}\hfill \rho ({\mathbf{x}}_{\mathbf{b}},t)[Z_2(\mathbf{u}({\mathbf{x}}_{\mathbf{b}},t))+Z_5(\mathbf{u}({\mathbf{x}}_{\mathbf{b}},t))+Z_6(\mathbf{u}({\mathbf{x}}_{\mathbf{b}},t))]={\tilde{f}}_4({\mathbf{x}}_{\mathbf{f}\mathbf{1}},t)+{\tilde{f}}_7({\mathbf{x}}_{\mathbf{f}\mathbf{2}},t)\\ \hfill +{\tilde{f}}_8({\mathbf{x}}_{\mathbf{f}\mathbf{3}},t)-{\omega }^{\prime }\left[{\tilde{f}}_2^{neq}({\mathbf{x}}_{\mathbf{b}},t)+{\tilde{f}}_5^{neq}({\mathbf{x}}_{\mathbf{b}},t)+{\tilde{f}}_6^{neq}({\mathbf{x}}_{\mathbf{b}},t)\right]\end{array}$$

(27)

After completing the collision process, both the equilibrium and non-equilibrium distribution functions at the fluid nodes are known. Therefore, for the velocity boundary, the revised density at the boundary nodes can be finally calculated as shown in Equation (27). This is substantially different from the original NEE boundary method in which the density at the boundary node ${\mathbf{x}}_{\mathbf{b}}$ is directly extrapolated from the neighbouring fluid node ${\mathbf{x}}_{\mathbf{f}\mathbf{1}}$. As an application of this revised density term, the complete post-collision distribution function could be obtained at the boundary nodes, and the mass conservation of the whole flow domain was also validated.

3. Results and Discussion

3.1. ML with the NEE Boundary Condition

The NEE boundary condition was applied for the simulation of several fluid flows. Although reasonable results could be obtained, it was accidentally found that mass conservation was not attained and maintained. In other words, once an external force field exists in the flow domain, ML will increase and the total mass decreases drastically. In order to prove the ML and the effects of an external force field on the variation of mass flux, a steady incompressible flow in a multiple right bending pipe was simulated and the force imposing this flow was provided by the gravitational force. As shown in Figure 3, periodic boundary conditions were applied to both the inlet and outlet of the pipe, and the bold solid lines represent the boundary walls on which the NEE boundary condition is applied. In order to compare the ML on each wall, Figure 3b indicates the labels that enclose the boundary walls. In addition, to observe the effect of body force inside the flow domain, the vertical direction ’down’ is assumed to be the direction of gravitational force, which has a magnitude of ${10}^{-4}$. The initial density in the fluid domain is set to 1.0. The Reynolds number used in the simulation is 37, which is based on the relaxation time 0.92 and reference vertical velocity 0.12.

Figure 4 shows the contour plots of horizontal and vertical velocity using the NEE boundary condition. In the upper domain of the multiple right bending pipe, the flow first moves down and turns right towards the positive horizontal direction (x-direction), whereas the flow direction is the opposite in the lower domain. As earlier indicated, the mass conservation should be checked and verified under the condition having the body force, i.e., the gravitational force. Therefore, the variation of the non-dimensional mass inside the domain and effects of the external force magnitude are shown in Figure 5. In this figure, $M_0$ indicates the initial mass in the domain as reference value, whereas M the total mass. The magnitude of the external force varies around a factor of 2.5, which includes 0.00004 and 0.0001. As shown in the figure, it is found that the decrease in mass is almost linearly proportional to the time step. In addition, with a larger magnitude of gravitational force (i.e., $g={10}^{-4}$), ML on the boundary becomes even stronger, which agrees well with the former estimation of ML in Equation (24). More importantly, the above facts have clearly demonstrated the existence of ML in the NEE boundary condition.

In the flow simulation inside the multiple bending pipe, it would be advantageous to know the places where ML occurs, such as the vertical and horizontal boundaries. Among them, some salient planes were chosen: B1, B2, B3, B4 and B5 (see Figure 3b). Figure 6 shows the variation of the non-dimensional ML with time step on these boundaries. Interestingly, it has been found that the ML is almost linearly proportional to the time step, which indicates a constant ML per time step. In addition, the ML on the horizontal boundaries B2 and B4 is larger than that on the vertical boundaries B1, B3 and B5. On the horizontal boundary B2, the ML is negative, which indicates that the external fluid would move into the flow domain of the pipe. Regarding the vertical boundaries, the ML on B1 is larger than that on the boundaries B3 and B5.

Therefore, when an external force is imposed on the flow domain based on the NEE boundary condition, it can be concluded that the ML on the boundaries having a perpendicular force is larger than that on other directional boundaries. In the bending pipe, those vertical boundaries are parallel to the gravitational force, which indicates that this force does not have any effect on the ML. However, it is found that ML still exists on the vertical boundaries B1, B3 and B5.

As discussed above, when the velocity on adjacent fluid nodes is not zero, the ML on the boundary would be stronger. Therefore, the perpendicular velocity profiles and ML at 10 time steps on the vertical boundaries B1 and B5 are analysed and shown in Figure 7. In the figure, the axis of the abscissa (the boundary sections B1 and B5, i.e., L) is all non-dimensionalised by the length of each section (i.e., $L_0$), whereas the axis of the ordinate (ML) is non-dimensionalised by the loss of total mass in the domain (i.e., $M_0-M_c$), where $M_c$ is the total mass after making the calculation converge. In the figure, $V_0$ indicates the reference mean velocity at the inlet section. As shown in the figure, it is found that the horizontal velocity components, which are the perpendicular velocities on boundaries B1 and B5, do not become zero, particularly at the corner edge. In addition, the magnitude of the non-dimensional ML increases with the increasing magnitude of horizontal velocity. Owing to the larger magnitude of the horizontal velocity on the edges, its ML also becomes higher, particularly at $L/L_0=1$ of section B1 and $L/L_0=0$ of section B5.

In addition, to present the explicit relationship between horizontal velocity and ML on the boundary, Figure 7 clearly shows that the magnitude of horizontal velocity around the boundary edge becomes larger than that at other positions. This fact can be explained with Figure 8, which shows the streamline pattern inside the flow domain having almost the same for the original NEE and our newly developed boundary methods. This figure indicates the general flow and wake structures, including the recirculation regions, while the flow passes through the multiple bending pipe. Not surprisingly, the recirculation regions are generated by the flow separation owing to the negative pressure gradient around edges and corners of the bending pipe. The recirculation flow forms a vortex region behind these separation points, which is believed to induce the increase in perpendicular velocity on the boundary, particularly on the boundary edge. In addition, the ML around the boundary edge should be larger than at other positions.

Figure 9 shows the ML along the left boundaries of the bending pipe. Comparing it with the streamline pattern, it also shows the clear variation of a higher magnitude of ML around the boundary edge.

3.2. Mass-Conserved Boundary Condition

In the newly developed mass-conserved boundary condition, only the density term of the equilibrium distribution function has been revised. Therefore, it is rather straightforward and clear to implement this method in the condition of velocity boundary. In this section, first, the steady incompressible flow in the multiple right bending pipe was simulated with the newly developed mass-conserved boundary condition. Then several benchmark test problems were validated to evaluate the accuracy of this newly developed boundary condition.

3.2.1. Steady Flow in the Multiple Right Bending Pipe

As discussed above, in the NEE boundary condition, ML appears on the boundaries of the multiple right bending pipe for a steady incompressible fluid flow, in particular when an external force field exists in the flow domain. The magnitude of the ML at the boundary nodes differs from each other, which results from the variation in the perpendicular velocity profiles along the boundary. To keep the mass conservation in the flow domain, the newly developed mass-conserved boundary condition was applied to simulate a steady flow in the multiple right bending pipe. The magnitude of the gravitational force is ${10}^{-4}$. As shown in Figure 10, it is found that the mass conservation is kept with the newly developed boundary condition, whereas the total mass decreases almost linearly with time step in the NEE boundary condition.

In addition to the comparison of mass conservation, the perpendicular velocity profiles at the fluid nodes adjacent to the boundary were also compared. In Figure 11, $V_0$ indicates the reference mean velocity at the inlet section. As shown in Figure 11, it is found that the magnitude of the perpendicular velocity with the newly developed boundary condition is smaller than that in the NEE boundary condition. Moreover, in the newly developed mass-conserved boundary condition, the perpendicular velocity profiles are much closer to zero, which indicates that there is no mass flux through the boundary.

Figure 12 shows the contours of non-dimensionalised pressure ($p/p_0$) in the flow domain for both the newly developed and NEE boundary conditions. As shown in Figure 12, it is found that the overall characteristics of pressure distribution are almost the same in both cases. However, in the NEE boundary condition, owing to the definition of pressure in Equation (8) and ML in the flow domain, the magnitude of the pressure field is decreased (by a factor of two, in this case). Regarding the particular features of pressure distribution, the magnitude of pressure around the convex corners is smaller than in other regions. On the contrary, the magnitude of pressure around the concave corners of the perpendicular boundaries B4–B5 and B7–B8 is the largest in the flow domain, which results from the main vertical fluid flow in the multiple right bending pipe.

3.2.2. Poiseuille Flow

The Poiseuille flow is a representative laminar flow based on the pressure drop in an incompressible, pressure-driven, and Newtonian flow in a long channel. In the LBM simulation, a periodic boundary condition is applied to both edges (i.e., inlet and outlet of the channel), whereas the top and bottom solid walls adopt the newly developed mass-conserved boundary condition. The analytical solution of the Poiseuille velocity profile can be expressed as a formula for the parabolic shape, as a function of the distance from the centre of the channel,

$$u=U_0\frac{H-y}{H^2/4}$$

(28)

where H is height of the channel and $U_0$ is the maximum streamwise velocity along the centreline (also called a centreline velocity), which could be calculated from the magnitude of force [30].

$$U_0=\frac{F{H}^2}{8{\rho }_0\nu }$$

(29)

where F is the external force for driving the flow in the channel. The kinematic viscosity can be calculated from the Reynolds number $Re=H{U}_0/\nu$. In the simulation, the relaxation time and Reynolds number were set to be 1.1 and 10, respectively, for all cases; then, the force and kinematic viscosity can be obtained from the equations as described above. To assess the convergence of a steady state in calculating the Poiseuille flow, the following criterion was used:

$$\sum _{i,j}\left|{\mathbf{u}}_{ij}^n-{\mathbf{u}}_{ij}^{n-100}\right|/\sum _{i,j}\left|{\mathbf{u}}_{ij}^n\right|\le 1.0\times {10}^{-8}$$

(30)

where ${\mathbf{u}}_{ij}^n=\mathbf{u}(x_i,y_i,n\Delta t)$, and the summation was implemented for the whole flow domain.

Figure 13 shows the results of the velocity profiles in a channel with height of 100. The LBM results are obtained based on the newly developed mass-conserved boundary condition. The maximum streamwise velocity can be obtained from the above definition of Reynolds number, and it is found that $U_0=0.02$. It is also found that all the calculations converged to steady state, and the streamwise velocity profile simulated by the LBM is in good agreement with the analytical solutions at steady state.

In order to evaluate the accuracy of the new mass-conserved boundary condition for the steady Poiseuille flow, the relative $L_2$-norm is also calculated. As shown in Figure 14, the axis of the abscissa denotes the channel height in log scale, whereas the ordinate shows the accuracy of the discrepancy in log between the streamwise velocity profiles for the LBM and the analytical solution (i.e., the relative global error), which is defined as follows:

$$E=\frac{\sqrt{{\sum }_y{(u(y)-u_e(y))}^2}}{\sqrt{{\sum }_y{u}_e{(y)}^2}}$$

(31)

where $u(y)$ is the streamwise velocity obtained from the LBM simulation and $u_e(y)$ refers to the analytical solution. As shown in Figure 14, it is found that the relative global error decreases with the increase in channel height, which perfectly fits the logarithmic curve of relative global errors. In the curve fitting, the slope of the fitted line is about 2.1. This indicates that the mass-conserved boundary condition owns second-order accuracy in space, which is consistent with the lattice Boltzmann model.

3.2.3. Lid-Driven Flow in Square Cavity

In order to validate the newly developed mass-conserved boundary condition, the steady Poiseuille flow was not only generated and tested, but the 2D lid-driven flow in a square cavity was also evaluated, which is a classical benchmark problem for numerical models and boundary conditions. The simulation was carried out for two different Reynolds numbers, 400 and 1000. The kinematic viscosity was also obtained in the definition of the Reynolds number. The adopted mesh size was $257\times 257$ in the cavity domain. The newly developed mass-conserved boundary condition was applied to the four enclosed boundaries of the square cavity, including the top moving wall and the other three static walls. In the boundary condition of the top moving wall, the streamwise surface velocity was set to $U_0=0.1$. As shown in Figure 15, not surprisingly, the current simulation is in good agreement with the existing results of Ghia et al. [31].

In addition to the application of the newly developed boundary condition, the NEE boundary condition was used to validate the cavity flow for the Reynolds number of 400, and the velocity profiles were obtained. Figure 16 shows the variation in the streamwise velocity profiles between the newly developed and NEE boundary conditions. As shown in the figure, the difference in the streamwise velocity profiles was only around 0.67%. This fact would imply that the difference is not substantial and can be negligible. Hence, the revision of the density term in the equilibrium distribution function can not only improve the accuracy at the boundary node but also retain most advantages of the original NEE BC.

The variation of the total mass in the whole flow domain was also evaluated with the aim to find the effect of the mass-conserved boundary condition, as shown in Figure 17. It is found that the newly developed boundary condition can keep the mass conservation, whereas the NEE boundary condition demonstrates a constant ML, which also indicates that the ML still exists in the NEE boundary condition even without the external force.

4. Conclusions

In this paper, a newly developed second-order mass-conserved boundary condition is proposed for the LBM. First, we analysed the source of ML in the NEE boundary condition, which results from the density term on the boundary node. From the investigation of a steady flow driven by gravitational force in a multiple right bending pipe, it is found that the ML on the boundaries directly against the flow is larger than that on the other boundaries. Although the external force is equivalent at each node in the flow domain, the ML at single nodes along the boundaries is different. This fact results from the different density distributions and velocity profiles along the boundary. The nonlinear variation of the density and velocity profiles along the boundary can induce a substantial ML. To understand the effect of the gravitational force on the ML in the flow domain, two different gravitational force fields (i.e., ${10}^{-4}$ and $4\times {10}^{-5}$) were applied. It is found that a larger external force has greater effect on the decrease in total mass in the flow domain.

Furthermore, a newly developed mass-conserved boundary condition is proposed, and the steady flow in a multiple right bending pipe is simulated with this new boundary condition. In the calculation, the magnitude of velocity profile on the fluid nodes adjacent to the boundary is reduced and approaches to zero, which meets the requirement of mass conservation. In addition, several benchmark test problems are coaxially validated for its accuracy and robustness. In this newly developed boundary method, the density term is revised for the equilibrium distribution function on the boundary nodes, and thus, it can be straightforwardly implemented. Except for a discrepancy in the mass conservation with the NEE boundary condition, the results based on the newly developed boundary condition show a good agreement with other existing results. The requirement of mass conservation in the whole flow domain should be kept for obtaining high-quality flow characteristics.