A comparative study of boundary conditions for lattice Boltzmann simulations of high Reynolds number ﬂows

Four commonly-used boundary conditions in lattice Boltzmann simulation, i.e. the bounce-back, non- equilibrium bounce-back, non-equilibrium extrapolation, and the kinetic boundary condition, have been systematically investigated to assess their accuracy, stability and eﬃciency in simulating high Reynolds number ﬂows. For the classical lid-driven cavity ﬂow problem, it is found that the bounce-back scheme does not inﬂuence the simulation accuracy in the bulk region if the boundary condition is properly imple- mented to avoid generating non-physical slip velocity. Although the kinetic boundary condition naturally produces physical slip velocity at the wall, it gives overall satisfactory predictions of the center-line ve- locity proﬁle and the vortex center locations for the Reynolds numbers considered. For the cavity ﬂow problem, all four boundary conditions show minimal difference in the computing time needed to reach a steady state. This is surprising because the kinetic boundary condition is signiﬁcantly different from the other three schemes which are designed speciﬁcally for no-slip boundary conditions. The bounce-back scheme is the most computationally eﬃcient in updating boundary points, which is particularly attrac-tive if there are a large number of solid bodies in the ﬂow ﬁeld. For the numerical stability, we further test the pressure-driven channel ﬂow with or without a enclosed square cylinder. Overall, the kinetic boundary condition is the most stable of the four schemes. The non-equilibrium extrapolation scheme presents excellent stability second to the kinetic boundary condition for the lid-driven cavity ﬂow. In comparison with other threes schemes, the stability of non-equilibrium bounce-back scheme appears to be less satisfactory for both ﬂows.


Introduction
The lattice Boltzmann (LB) method has been developed into an efficient mesoscopic simulation tool for fluid dynamics [1] , which has shown its strength in simulating multi-phase and multicomponent flows [2] , and flows in porous media [3] . Its potential for simulating high Reynolds number ( Re ) flows has also attracted significant interest. For example, with a second-order numerical scheme in both time and space, the LB method was shown to perform well for decaying turbulence, although it may require higher spatial resolution, with a correspondingly smaller time step, in comparison to a spectral method [4] . However, a high-order LB model has been successfully applied for simulating the Kida vortex flow [5] . There have been various attempts to simulate turbulence including large eddy simulation and direct numerical simulation [6] .
It is of importance to choose an appropriate boundary condition for simulating high-Re flows as it is often the key to stability, efficiency, and accuracy of simulations, particularly for complex geometry [7] . Conventionally, we need to implement the macroscopic no-slip boundary condition, which has many different implementations in the LB method because of its mesoscopic nature. The most commonly-used implementation is the bounce-back (BB) scheme, which originally describes a stationary wall. However, with a simple modification, it can be used for moving walls as well [8] . By assuming bounce-back of the non-equilibrium part of the distribution function, Zou and He derived the so-called non-equilibrium bounce-back (NEBB) boundary condition [9] for both moving and  [10,11] proposed the non-equilibrium extrapolation (NEEP) scheme, which is also suitable for both moving and stationary walls. There are also other implementations, e.g., the counter-slip boundary [12] .
The kinetic boundary condition (KBC) is of great interest for the LB method [13] . This scheme can induce physically realistic slip velocity at a wall boundary, and has been often used in the discrete velocity method (DVM) [14] for rarefied gas flows, e.g., [15] . As a special form of DVM, the LB method can also use the KBC to simulate rarefied gas flows over a broad range of Knudsen numbers ( Kn ) [16] . Moreover, the KBC has a nice property of retaining the positivity of outgoing particle distribution functions at the boundary nodes if the incoming mass flux is positive, which is helpful for numerical stability. Therefore, it is interesting to investigate its performance for high-Re flows.
The accuracy of the BB scheme has been studied in detail for simple flows, e.g. [17,18] . It was shown that the BB scheme may induce a numerical slip velocity and cause different order of errors depending on its implementation (e.g.,"on-grid" or "halfway"). While this artificial slip velocity can be eliminated for simple geometries [19] , we will investigate the performance of the BB scheme without the deficiency.
Other boundary schemes are also assessed in the literature including those specifically designed for curved boundaries (see e.g., [20][21][22] ). Of particular interest is that the NEEP scheme shows second-order accuracy for "moderate" Reynolds number flows in complex domains [20] , which is consistent with the analysis in [10,11,18] for relatively simple flows. Moreover, the scheme also presents sufficient numerical stability for the investigated Reynolds numbers.
In this work, we will focus on accuracy, efficiency, and stability of the BB, NEEP NEBB and the KBC schemes for high Reynolds number flows. Our benchmark simulations will be conducted using a D2Q9 lattice model for the classical lid-driven cavity flow and the flow around a square cylinder confined in a channel. The results will be particularly focused on the cavity flow where we are able to assess the four boundary schemes for Reynolds numbers up to 7500. For the flow around a square cylinder, we will primarily investigate their numerical stability.

Lattice Boltzmann equation and D2Q9 lattice model
The LB method can be considered as an approximation to the Boltzmann-BGK equation [23][24][25] . The governing equation can be written as This equation describes the evolution of the distribution f α (x , t ) for the discrete velocity c α at position x and time t . The complex molecular interaction is simplified as a relaxation scheme towards the discrete equilibrium distribution f eq α (x , t ) . In order to simulate incompressible and isothermal flows, it is common to use an equilibrium function with second-order velocity terms which is determined by the density, ρ, the fluid velocity, u , , and the reference temperature, T 0 . For gas flows, the constant, R , can be conveniently understood as the gas constant. If a liquid is involved, R , together with T 0 , can be considered as a reference quantity. For convenience, the sound speed c s is often considered to be RT 0 , although there is a constant factor √ γ of difference where γ is the ratio of specific heats. The weighting factor is denoted by w α for the discrete velocity c α . Although the LB method appears to be a primary tool for modeling gas flow, based on its origin from kinetic theory, it can also be used for liquid flows. This utilizes the fact that, under appropriate conditions, the Navier-Stokes (NS) equations can be recovered, which is valid for both gas and liquid flows.
The relaxation time τ is related to the fluid viscosity via the Chapman-Enskog expansion, i.e., µ = pτ, where µ is the viscosity and p the pressure. Hence, for isothermal and incompressible flows, the Reynolds number becomes Re = ρ 0 u 0 L/µ = u 0 L/ (τ RT 0 ) , where the subscript 0 denotes the reference value and L is the characteristic length of the system. It is worth noting that the important non-dimensional number, Kn , can be defined as µ 0 RT 0 / (p 0 L ) . So, τ is also related to the Knudsen number through viscosity, i.e., Kn = τ RT 0 /L where p 0 = ρ 0 RT 0 . We can therefore obtain Kn × Re = u 0 / RT 0 = Ma . For gas flows, the Knudsen number is the ratio of the gas mean free path and the characteristic length, which measures the rarefied level of gas flows. If a liquid flow is simulated, we may not be able to define a physically meaningful Knudsen number. However, the above defined Knudsen number can be regarded as an important "numerical" number for LB simulation due to its kinetic nature. In this case, this "numerical" Knudsen number will influence the dynamical behavior of modelling liquid flow. When Kn becomes large, the LB method may deviate from the NS dynamics significantly, even in the bulk region away from any surface [26] . Therefore, for liquid flows, it is important to ensure a small "numerical" Kn in any LB simulation in order to avoid deviating from the NS hydrodynamics. However, this is not a problem for gas flows. Instead, it provides an opportunity to model rarefied gas flows with suitable ranges of Knudsen numbers (see e.g. [16] ).
For the isothermal case, only density and velocity are of interest and can be obtained as A smart scheme can be utilized to numerically solve Eq. (1) and achieve the particle-jump like simulation [27] . By applying the trapezoidal integration rule for both sides of Eq. (1) and introducing a transformation we obtain the discretized form for ˜ f as which provides the stream-collision algorithm. The last key step is to choose an appropriate set of discrete velocities. The D2Q9 lattice is the most popular choice for twodimensional flows, which employs nine discrete velocities The corresponding weights are w 0 = 4 / 9 , w 1 −4 = 1 / 9 , w 5 −8 = 1 / 36 .
As discussed, the stream-collision algorithm can readily be implemented now and the only requirement is to tie the space and time step together. If one is determined, the other is determined also. For instance, assuming the system length is L , we may set the space step dx and dy as dx = dy = L/ 10 0 and then dt = L/ (100 3 RT 0 ) with 10 0 cells of spatial resolution. This insures the "particles" are jumping on a uniform grid system. Alternatively, Illustration of south (bottom) wall boundary.

Table 1
Parameters for two grids of a physical system defined by u 0 , ν 0 and T 0 , where we may transform Eq. (1) to its non-dimensional form first and then apply the scheme given by Eq. (5) , as done in [26] . Obviously, the non-dimensional transformation will not alter the final simulation results.
It is common to choose the space step dx or dy and time step dt as reference values, i.e., the "lattice" unit. For the D2Q9 model, the transformation is where the hat symbol is used to denote normalized variables in the lattice unit. Also, we note dx = dt 3 RT 0 . With the transformation, Eq.
For convenience, a new relaxation time ˆ In addition, it is also common practice to calculate the important non-dimensional variables in the lattice unit, i.e., where ν denotes the kinematic viscosity and we have ν = ˆ νdx 3 RT 0 . The length ˆ L in the lattice unit becomes the cell number of one coordinate and is commonly denoted as N . Finally, the relation between the viscosity and relaxation time in the lattice unit is While any proper transformations will not influence the final results, we find that an issue may occur with the lattice unit when conducting grid independent tests. For convenience, we consider isothermal and incompressible flow. Following the aforementioned procedure, we list parameters in both the physical unit and lattice unit of two different grid systems in Table 1 . The physical system is defined by a set of fixed parameters u 0 , ν 0 and T 0 . It shows that, to keep all the parameters exactly the same in physical units, the relaxation time and viscosity in the lattice unit should depend linearly on the mesh size. It is different from the practice where the viscosity ˆ ν in the lattice unit is fixed. From the above transformation and Table 1 , if ˆ ν and Re are fixed, the corresponding phys-ical system needs to have either a different viscosity or a different characteristic length while the characteristic velocity must be changed. This turns out to be, however, inconsistent with common CFD practice. In fact, since both Ma and Kn are changed in this way, there is a risk that Kn for the coarsest mesh becomes so large (cf. Eq. 9 ) that the LB method may deviate from the NS dynamics, even in the bulk [26] . For this reason, we shall keep all the parameters unchanged in physical units for grid independent tests or scheme accuracy analysis.

Boundary conditions
We choose to discuss the KBC and the BB, NEBB, NEEP schemes, and their implementation is illustrated in Fig. 1 . For any boundary scheme, the aim is to determine the outgoing distribution func- where the superscript b is used to denote boundary point, according to the incoming distribution functions f b 7 , f b 4 and f b 8 known from the stream step. It is worth noting that we will use the "on-grid" implementation for all four boundary conditions.
The main idea of the BB scheme is that an incoming particle distribution directing to the wall bounces back into the fluid domain in the opposite direction. Therefore, for the D2Q9 model, we As discussed in [19] , the nonphysical slip velocity can be easily eliminated for this cavity flow by taking care of f b 1 and f b 3 . The KBC assumes that an outgoing particle completely forgets its history and its velocity is renormalized by the Maxwellian distribution. Its implementation for the LB method has been discussed not only for standard models [13] but also for multi-speed models [28] . For the D2Q9 model as shown in Fig. 1 , the rule can be simplified as, where the wall density ρ w is determined according to mass conversation and is actually canceled in the above formula.
In the NEBB scheme [9] , the unknown distribution functions are solved by setting the constraints of The equations are indefinite since there are more unknowns than equations. To resolve this, the distribution is first divided into the equilibrium and non-equilibrium parts and the bounce-back scheme is assumed to be valid for the non-equibrilium part of distribution function normal to the wall, i.e., With this assumption, all unknown distribution functions can be determined.
The last boundary condition is the non-equilibrium extrapolation scheme [10,11,29] , which divides all the unknown boundary distribution functions into two parts: the equilibrium part f eq,b α and the non-equilibrium part f For the D2Q9 model, the equilibrium part is The non-equilibrium part is approximated as where ρ 1 and u 1 denote the quantities at the node ( i , 1), as shown in Fig. 1 . Obviously, a first-order extrapolation scheme is used for both the density and the non-equilibrium part. However, it has been shown that the NEEP scheme can give second-order accuracy for incompressible flows [11] .

Accuracy, stability and efficiency for lid-driven cavity flow
In the following we will compare accuracy, numerical stability and efficiency of these four schemes for simulating the classical lid-driven cavity flow problem. The fluid is confined in a square box but the top (north) lid of the box is moving with a constant speed. The Reynolds number and the top wall speed (i.e., Ma ) will be specified for each of the following simulation. We will mainly consider steady flows with Re smaller than 7500, so that we can focus on the effect of boundary conditions without interference of physical instability. Hence, the criterion is used to judge whether the steady state is reached, where u x and u y are the x -component and y -component of the velocity, u , respectively. In our simulations, the largest value of E is 10 −6 .
For a no-slip boundary condition, we need first to ensure noslip can be correctly achieved. For example, we need to eliminate any non-physical slip velocity at a boundary for the BB scheme [19] . Therefore, accuracy analysis for a scheme can focus on numerical error related to mesh size in a common convergence test.
It is found that predictions for all boundary schemes are close to the benchmark results, given the fact that the KBC will produce slip velocity at a wall boundary, see Fig. 2 for the velocity profiles through centerlines and Table 2 for predictions of the vortex center location. This confirms, as expected, that the KBC can be used for flows with very low Knudsen and high Reynolds numbers. In fact, Fig. 3 shows decreasing velocity slip with increasing Re (and decreasing Kn ) while Ma is fixed, which indicates that the impact of the slip velocity caused by KBC will become negligible for high Reynolds numbers and supports the findings in Fig. 2 and Table 2 . It is also worthwhile to note that the slip velocity produced by the KBC represents physical effects for gas flows. The BB boundary scheme also works well since the non-physical slip is removed here.
We now perform common convergence tests for accuracy analysis. Following the discussion in Section 2 , we keep all parameters unchanged in physical units. We check the results at a series of centerline points adjacent to the wall boundary, namely x = 0 . 5 and y = 1 / 32 , 1 / 16 , 1 / 8 , 1 / 4 . A moderate Reynolds number of 400 is considered and the top wall speed is set to be 0.1 in the lattice unit. Five sets of grids are investigated for each boundary condition, namely 96 × 96, 128 × 128, 192 × 192, 256 × 256 and 512 × 512. The results with the finest grid are considered to be "correct" to calculate the actual simulation accuracy, which is shown in Fig. 4 .
All the simulations show similar accuracy, although those using the NEEP scheme show slightly higher accuracy. It is interesting to see that the simulations with the BB scheme can achieve accuracy better than first-order. This reflects the fact that the convergence test is actually measuring the accuracy of second order of Eq. (5) (cf. Ref. [17] , particularly Eq. (10) in Page 121.) once the nonphysical slip is carefully removed for the BB scheme, even with an on-grid implementation. Moreover, in comparison with the other three schemes, the first-order extrapolation of the NEEP scheme appears to enhance the simulation accuracy in the bulk region.
To study the numerical stability, we carry out two different tests. For the first test, we use a fixed number of grid points to simulate different Reynolds number, and record the maximum Re that can be achieved stably; For the second, the Reynolds number is fixed with a difference number of grid points for the simulations. We then record the minimum grid point number needed for a stable simulation. The results are listed in Tables 3 and 4 . The KBC shows excellent stability, which can be attributed to the favorable property of positivity. The NEEP scheme also has good stability while the NEBB scheme is the less satisfactory.
The numerical efficiency of the boundary conditions is also measured by two criteria: the computational time needed for approaching the steady state determined by setting E = 10 −6 in Eq. (17) and the computational time for updating the wall boundary points. For this purpose, we run simulations on an Intel Core i5 3.20GHZ CPU. The results for the first criterion are listed in Table 5 . There is no significant difference between the four schemes. This is particularly interesting for the KBC because previous simulations, using a discrete velocity model, became increasingly slower as the Knudsen number decreased [15] for rarefied gas flows. If we consider the lattice Boltzmann model to be a special form of the discrete velocity model, the evidence here shows that it presents little difference to those specially designed boundary conditions for continuum flows.
In Table 6 , we list the computational time for updating all boundary points at the four walls during one time step. It is found that the NEEP boundary condition needs the longest time, which is around 40 times more expensive than the most efficient BB scheme. This shows that the BB scheme should be an efficient choice for simulating flows with a large percentage of solid boundary points, such as particle-fluid multiphase flows.

Stability for flow around a square cylinder
A two-dimensional channel flow around a square cylinder is utilized here to further assess the numerical stability of the four boundary schemes. The channel length and height is set to be L and H , respectively. The square is of width D and its center is located at the point (4.5 D , H /2). For convenience, we choose the channel height H as the reference length to define the Reynolds number as well as the Knudsen number. The flow is driven by the pressure difference between the inlet and outlet of the channel. The pressure profiles at the inlet and outlet are specified using a first-order extrapolation scheme which is implemented in the way described in, e.g., [32] . While both the NEEP scheme and the NEBB scheme can be used as the pressure flow boundary, they will not be employed in our comparisons to avoid the unnecessary coupling effect. Also, to eliminate the impact of transient effects on the comparison, steady cases will be considered here so that the Reynolds number will not be as high as those in the lid-driven cavity flows. However, to make the stability tests sensible, we will This observation is generally consistent with previous findings. In [11] , it is found that ˆ τ o ≈ 0 . 51 is the threshold if the pressure inlet/outlet is also implemented using the NEEP scheme, while the critical value becomes approximately 0.7 if extrapolation schemes are used for both wall and flow boundary. Therefore, the less satisfactory stability behavior of the NEEP scheme may be attributed to the mixing with the extrapolation pressure flow boundary. We also note that the NEEP scheme appears to work well with a fixed velocity inlet boundary condition for similar flow problems [20] . The NEBB scheme also fails to achieve a stable simulation with 20 0 0 0 × 200 cells for Re = 2 . 16 .   For the flow around a square cylinder, we will mainly investigate the stability of the KBC and the BB scheme based on the observation of the pressure-driven flow. For this purpose, the channel length L is set to be 15 D and the height H equal to 2.5 D where the flow will feel the blockage effects. To test the stability, we setup a case with 300 × 50 cells and ˆ τ o = 0 . 500866 ( Kn = 0 . 00001 ). The  Table 5 Computational time (h) needed for approaching the steady state determined by setting the criterion, E = 10 −6 , given by Eq. (17) . The NEBB scheme is not tested for Re = 50 0 0 since the simulation is not stable. pressure difference is set to be 10 −8 ( × ρ 0 RT 0 ). Both schemes can achieve a stable simulation, where "stable" means the simulation lasts at least 60 0 0 0 iterations without the occurrence of "NAN" and afterwards we stop monitoring. This observation further verifies the stability of both the KBC and the BB scheme.

Re
To examine the impact of the slip velocity, we set the pressure to be 1.0 0 02 ( × ρ 0 RT 0 ) at the inlet and 1 ( × ρ 0 RT 0 ) at the outlet.    5 . The Reynolds number is estimated to be about 28 using the average speed at the inlet. The simulation with the KBC produces qualitatively similar results to the one with the BB scheme, which again indicates that the slip velocity will play a minor role if the Knudsen number is small and the Reynolds number is sufficient large.

Concluding remarks
We have investigated the accuracy, stability, and efficiency of four popular boundary conditions for LB simulation of flows with high Reynolds numbers. For the present study, the classical liddriven cavity problem is the primary benchmark application while the channel flow around a square cylinder is used to supplement the observation of numerical stability.
It is found that all four schemes give reasonable predictions of the centerline velocity profile and the vortex position for the liddriven cavity flow. In particular, the slip velocity given by the KBC has no significant effect on the Reynolds numbers considered (low Knudsen number). Once the non-physical slip is removed, the BB boundary condition also performs well. Moreover, the discretization accuracy is not influenced by the BB scheme in the bulk region. In fact, simulations using any of the four boundary conditions give similar orders of accuracy for the tested bulk points. The KBC shows the best stability for not only the cavity flow but also the flow around a confined square cylinder, possibly because it is able to maintain its positivity of distribution function. The NEEP scheme shows excellent stability for the cavity flow. The BB scheme shows consistent stability behavior for both flows. In the present study, we find that the computational time required to reach the steady state for the lid-driven cavity flow is similar for all four boundary conditions considered. Interestingly, the KBC scheme, which is often used for rarefied gas flows, shows no significant difference in comparison to the other three boundary conditions specially designed for continuum flows. However, both the KBC and NEEP schemes need a significantly longer time to update the boundary points while the BB boundary condition, as expected, is the most computationally efficient.