## Abstract

By means of the continuity equation of the incompressible Navier–Stokes equations, additional physical arguments for the derivation of a formulation of the no-slip boundary condition for the lattice Boltzmann method for straight walls at rest are obtained. This leads to a boundary condition that is second-order accurate with respect to the grid spacing and conserves mass. In addition, the boundary condition is stable for relaxation frequencies close to two.

## 1. Introduction

The lattice Boltzmann method has become increasingly popular for the simulation of incompressible fluid flow (e.g. [1–3]). It is based on a ‘particle’ approach, meaning that the fluid is considered as an ensemble of particles rather than continuous matter. These particles are, however, bound to travel only from one grid node to the next grid node where they undergo collisions between each other. In the method itself, this is manifested in the collision and streaming step. Opposed to the earlier lattice gas automata, cf. [2], the lattice Boltzmann method does not handle each particle individually. Instead, the central object of the lattice Boltzmann method is the particle distribution function *f*(** x**,

**,**

*c**t*), which can be seen as a kind of probability density of finding a particle at a certain position with a certain velocity at a certain point in time. Since we restrict the possible velocities for each particle to a finite set of nine velocities

*c*_{i},

*i*=0,…,8, in two dimensions, such that the particles can only move from one grid node to the next one during one time step, cf. figure 1 and table 1, the particle distribution function at a node can be reduced to nine discrete values

*f*

_{i}(

**,**

*x**t*)=

*f*(

**,**

*x*

*c*_{i},

*t*),

*i*=0,…,8, which are called the populations. The grid consisting of nodes

**=(**

*x**α*,

*β*)

^{T}, , and the interconnecting velocity vectors

*c*_{i}is called a lattice. Therefore, the resulting lattice in two dimensions with nine velocities is called the D2Q9 lattice. Accounting for the difference in length, lattice weights

*t*

_{i}, cf. table 1, are associated to each of the lattice directions.

The collision step at a node ** x** is then modelled as follows:
1.1
where

*f*

^{*}

_{i}is the post-collision value of the population

*f*

_{i}and

*Ω*

_{i}is the collision operator that is a function of all the populations at the node in question. The streaming step on the other hand shifts the post-collision values to the neighbouring grid nodes: 1.2 The collision operator

*Ω*is modelled as a relaxation towards an equilibrium distribution , the Bhatnagar–Gross–Krook collision operator (e.g. [3]): 1.3 where

*ρ*is the local density and

**is the local (mean) velocity. The equilibrium distribution used in the present discussion is given by 1.4 where is the speed of sound of the lattice and**

*u***Q**

_{i}=

*c*_{i}

*c*_{i}−(1/3)

*c*

^{2}

_{s}

**I**is a second-order tensor. The density

*ρ*and the velocity

**at a node are obtained as statistical moments of the populations. The density**

*u**ρ*is the zeroth moment of the populations: 1.5 Whereas the velocity

**is given by the first statistical moment: 1.6 The density**

*u**ρ*and the velocity

**are the actual quantities of interest, computed indirectly by the populations. It can be shown by a procedure called the Chapman–Enskog multiscale expansion (e.g. [3]) that the so-computed macroscopic quantities obey the incompressible Navier–Stokes equations up to second order with respect to the grid spacing if compressibility effects can be made negligible: 1.7 and 1.8 where the pressure**

*u**p*is given by

*p*=

*c*

^{2}

_{s}

*ρ*. An additional result of this expansion is that we can relate the leading order of the non-equilibrium part of the populations to the rate of strain tensor

**S**=(1/2)(∇

**+(∇**

*u***)**

*u*^{T}): 1.9 On the other hand, we can find an inverse relation by 1.10 Equation (1.10) or similar equations are used by different authors [4,5] in order to derive boundary conditions by first determining the density

*ρ*, the velocity

**and the rate of strain tensor**

*u***S**at the boundary node and then computing the populations before collision by 1.11 Opposed to ‘classical’ solvers based on the incompressible Navier–Stokes equations, the implementation of the no-slip boundary condition is not straightforward. For classical solvers, the value of the velocity is known for a no-slip boundary condition

**=0. However, for the lattice Boltzmann method, parts of the populations are unknown and a closure needs to be found for them in order to impose a no-slip boundary condition. In the example of a wall situated in the south of the computational domain as depicted in figure 1, after streaming, the populations**

*u**f*

_{1},

*f*

_{7},

*f*

_{8}are unknown and need be given a value in order to apply the collision step. Different closures have been proposed in the literature (e.g. [5]). The earliest use a particle–wall-collision-type argument where at the wall particles are reflected back into the directions where they were coming from, hence the name bounce back (e.g. [4]). The collision step is left out for this type of boundary condition. This boundary condition is mass conserving and has programming advantages over all other boundary conditions, since it does not depend on the orientation of the wall. However, it does not represent the position of the wall very accurately. The resulting wall position happens to be not at the boundary nodes but somewhere between the boundary nodes and the first fluid nodes. The idea of reflection plays also a central role in the boundary conditions of Inamuro

*et al*. [6]. However, they borrowed an argument from the kinetic theory of gases, the diffuse reflection. This assumption states that after contact with the wall, the particles are reflected with a probability density function that has the form of an equilibrium distribution centred around the wall velocity plus a slip velocity. Inamuro

*et al*. [6] used this argument to derive a boundary condition for the lattice Boltzmann method, which is second order with respect to the grid spacing.

Zou & He [7], on the other hand, developed an argument based on the Chapman–Enskog expansion. The Chapman–Enskog expansion indicates that to leading order the non-equilibrium parts of the populations along the grid axes are equal. This argument leads to a closure formulation for the unknown populations, resulting in a second-order accurate boundary condition.

Apart from the fact that both boundary conditions, the one of Inamuro *et al*. and of Zou & He, are not mass conserving, they show the appearance of numerical instabilities for values of the relaxation frequency *ω* close to two. More recently, an attempt has been made to unify both boundary conditions by use of the continuity equation of the incompressible Navier–Stokes equations (cf. [8]). This allowed in addition to track down the origin of the numerical instabilities. In the present discussion, we focus on some details of this derivation, explaining more carefully relevant steps of the derivation. The present discussion is organized as follows. Section 2 presents the derivation of the present boundary condition, which is then verified in §3. The present discussion is concluded in §4.

## 2. Derivation

A necessary condition for mass conservation during collision is that the sum of the non-equilibrium parts of the populations vanishes:
2.1
Similarly, in order to conserve momentum, we have
2.2
Equations (1.9), (2.1) and (2.2) lead to a system of six equations for the non-equilibrium parts of the populations:
2.3
This system has been found by Halliday *et al*. [9]. However, their conclusions and further steps differ from the further derivation in the present discussion (cf. [8]). For a southern wall at rest, as depicted in figure 1, the velocity is zero, *u*_{x}=0, *u*_{y}=0. In addition, this implies that the derivatives with respect to *x* vanish:
2.4
Owing to the continuity equation (1.8), the normal derivative of the normal component vanishes equally, ∂*u*_{y}/∂*y*=0. Since the mixed derivative ∂*u*_{x}/∂*y* is an unknown, the off-diagonal element of the rate of strain tensor *S*_{xy} is unknown. Therefore, we remove the fifth equation from the system (2.3) to get
2.5
A consequence of this system is the following system:
2.6
2.7
2.8
2.9
2.10
and
2.11
The rule behind the first three equations (2.6)–(2.8) can be illustrated by arranging the non-equilibrium parts of the populations in a 3×3 tableau (figure 2*a*) and requiring that for all populations the sum along the *y*-axis is zero. The last three equations (2.9)–(2.11) are then associated with a vanishing sum along the *x*-axis (figure 2*b*). For the southern wall in question, cf. figure 1, the populations *f*_{3}, *f*_{4} and *f*_{5} are among the incoming populations after streaming. The populations *f*_{1}, *f*_{7} and *f*_{8} are unknown. Although the populations along the wall, *f*_{0}, *f*_{2} and *f*_{6}, are given after streaming, we recompute them in order to satisfy the system of equations (2.6)–(2.11). A consequence of taking only *f*_{3}, *f*_{4} and *f*_{5} as initially given is that the density *ρ* at the wall is unknown. However, equation (2.11) implies that
2.12
Thus, using this condition, we can define the unknown density at the wall *ρ* as follows:
2.13
since the velocity at the wall vanishes. For the non-equilibrium part of the populations along the wall, *f*_{0}, *f*_{2} and *f*_{6}, we use the value of the leading-order term given by equation (1.10), using the elements of the rate of strain tensor computed before (cf. [8]), i.e.
2.14
Finally, the closure formulation of the populations *f*_{0}, *f*_{1}, *f*_{2}, *f*_{6}, *f*_{7} and *f*_{8} as functions of the incoming populations *f*_{3}, *f*_{4} and *f*_{5} for the present no-slip boundary condition is
2.15
2.16
2.17
2.18
2.19
2.20
and
2.21
After having computed *f*_{0}, *f*_{1}, *f*_{2}, *f*_{6}, *f*_{7} and *f*_{8} by equations (2.15)–(2.21), the usual collision and streaming steps are applied. The mass streamed out after collision is given by
2.22
Therefore, the total mass streamed out after collision equals the total mass streamed in before collision, leading to exact mass conservation at the node. As explained in Verschaeve [8], the condition of mass conservation can also be used as a starting point for the present derivation. In this case, equation (2.11) would be satisfied as a consequence of the mass conservation condition. We remark that the present boundary condition is mass conserving in the sense defined by Chopard & Dupuis [10], which includes only the bulk fluid nodes in the mass balance. Hollis *et al*. [11] formulated an alternative mass balance including the boundary nodes. However, their mass balance and the system (2.6)–(2.11) are difficult to satisfy simultaneously.

In the next section we present the results of benchmark tests, which were also performed by Verschaeve [8].

## 3. Numerical verification

We set up a two-dimensional channel flow using the present no-slip boundary condition for the channel walls. For different resolutions *N* of the channel width, we measured the *L*^{2} error using the analytical solution of Poiseuille flow and oscillating Poiseuille flow at a Reynolds number of 10 (e.g. [12]). The results can be seen in figures 3 and 4, which show that the present boundary condition is indeed second-order accurate. The time step was reduced proportional to the square of the grid spacing,
3.1
In addition, we performed simulations of Poiseuille flow with an increasing relaxation frequency *ω* until the simulations became unstable. This happened for values of *ω* larger than 1.995. See [8] for details and other benchmark tests.

## 4. Conclusions

We derived a no-slip boundary condition for straight walls, which is second-order accurate, mass conserving and stable for relaxation frequencies close to two.

## Footnotes

One contribution of 26 to a Theme Issue ‘Discrete simulation of fluid dynamics: methods’.

- This journal is © 2011 The Royal Society