Description of the lattice Boltzmann method

To begin with, we discretize the space R^d by a regular cubic lattice using a uniform grid size h. The grid points are xj = hj with j ∈Z^d. Collecting all the grid points in ¯Ω, we get a discretization of the domain Ω. Likewise, we discretize the time domain similarly by placing a grid on the temporal interval [0,∞) with grid spacing △t and grid pointstn=n△t,n∈N₀.

In the corresponding lattice Boltzmann setup to simulate (2.1), the time step is △t=h². This ratio between time and space step is related to the diffusive scaling of the kinetic equations. For details we refer to [42].

On this lattice, the update rule of the lattice Boltzmann method is described by

f_i(n+ 1,j+c_i) =f_i(n,j) +J(f)_i(n,j) +g_i(n,j), (2.20) where V ={c₁, . . . ,cN} ⊂R^d is the finite discrete velocity set. The numbers fi(n,j) represent the particle distributions related to velocity c_i at time level tnand nodexj. The functiongi models the body force (in appendix D, the set V,f^∗ andcs for several well known models are given)

g_i(n,j) =c⁻²_s h³f_i^∗c_i·G(t_n,x_j).

J is the collision operator, and chosen in this work to be of general relaxation type

J(f) =A(f^eq(f)−f), (2.21) where A is a linear mapping and f^eq is a so-called equilibrium function. The common used models for this kind of relaxation type collision include the fa-miliar BGK model and MRT. Obviously, the collision is determined by several basic ingredients: the discrete velocity setV, the equilibriumf_i^eqand the linear

2.2. DESCRIPTION OF THE LATTICE BOLTZMANN METHOD 15 mappingA. Without restriction to some specific model, we state our assump-tions containing the widely used models such as BGK and MRT.

The velocity setV admits a symmetry property, i.e.,

V=−V. (2.22)

Those velocity sets for the well known models D2Q9, D3Q15, D3Q19 and D3Q27, which are displayed in appendix D, all possess this property. We use the equilibrium function recommended in [32] for the incompressible lattice Boltzmann model, which is based on the assumption that the fluid density slightly fluctuates around a constant ¯ρ. In this work we set, without loss of generality,

¯

ρ= 1. (2.23)

Hence the equilibrium function is of a polynomial form:

f_i^eq(f) =F_i^eq(ˆρ,u) =ˆ f_i^∗ with respect to the total mass density ˆρ and the average velocity ˆu (or more precisely, the average momentum ¯ρu),ˆ

ˆ called constant equilibrium distribution, obeys the symmetry property

f_i^∗=f_i^∗∗, if c_i =−c_i∗, (2.26) and is identified by the following constraints

N For the convenience in the later use, we split the equilibrium function into two parts

f_i^eq(f) =f_i^L(f) +f_i^Q(f,f), (2.29) namely a linear part f_i^L and a quadratic part f_i^Q, here f is the vector with components offi. Defining them in a concise way, the linear part is written as

f^L(f) =F^L(ˆρ,u),ˆ F^L(ˆρ,u) = (ˆˆ ρ+c⁻²_s uˆ·V)f^∗,

wheres∈ F andh1,Vsi= ˆw. Due to the properties off_i^∗ in (2.27) (2.28), we can find out

h1,F^L(ˆρ,u)ˆ i= ˆρ, (2.30) h1,VF^L(ˆρ,u)ˆ i= ˆu, (2.31) h1,V⊗VF^L(ˆρ,u)ˆ i=c²_sρI,ˆ (2.32) and

h1,F^Q(ˆu,w)ˆ i= 0, (2.33)

h1,VF^Q(ˆu,w)ˆ i=0, (2.34)

h1,V⊗VF^Q(ˆu,w)ˆ i= ˆu⊗w.ˆ (2.35) After giving the equilibrium function f_i^eq, conditions to determine the linear mappingAare required in order to completely fix the collision operator. These are:

(i) A is symmetric;

(ii) Ais positive semi-definite;

(iii) K ={1,v₁, . . . ,v_d} generates the kernel of A.

(iv) AΛf^∗ =c²_s/µΛf^∗ with Λ =V⊗V−1/d|V|²I and µ=ν+c²_s/2,

where condition (iv) means that the components of Λf^∗ are eigenvectors of the collision matrixA with eigenvaluec²_s/µ.

Remark 1 Let Q be the orthogonal projection onto the kernel of A and P :=

I−Qthe projection on the complement. Then we define A^†= (A|PRN)⁻¹P to be the peseudoinverse of A, and A^† has the following properties:

QA^†=A^†Q= 0, P A^†=A^†P =A^†, AA^†=A^†A=P (2.36) Remark 2 The so-called BGK collision operator J(f) = _τ¹(f^(eq) −f) is a special case considered here. A = ¹_τP with τ = µ/c²_s is a particular choice which satisfies all the above conditions (i) to (iv), and A^†=τ P. Moreover

P(f^(eq)−f) = (P +Q)(f^(eq)−f) = (f^(eq)−f), (2.37) since f^(eq)−f is orthogonal to the kernel of A which is easily checked by ob-serving (2.30) to (2.35).

Since the equilibrium functions have a linear and quadratic part and A is a constant matrix, the collision operator also consists of two parts,

J(f) =J^L(f) +J^Q(f,f), with a linear collision operator

J^L(f) =A(f^L(f)−f), (2.38)

2.2. DESCRIPTION OF THE LATTICE BOLTZMANN METHOD 17 and a quadratic operator

J^Q(f,f) =Af^Q(f,f). (2.39) Besides, observing that the term concerned with the quadratic collision operator is the only nonlinear part in the lattice Boltzmann scheme (2.20), hereafter the lattice Boltzmann method with only a linear collision operator is called the linear lattice Boltzmann method, whereas the lattice Boltzmann method with both linear and quadratic collision operators is called nonlinear lattice Boltzmann method. In this work, if not particularly pointed out, the results hold for the nonlinear case usually and for the linear case by dropping the quadratic equilibrium function.

While equipped with proper initial valuesfi(0,j), xj ∈Ω and boundary con-¯ ditions at boundary nodes, the lattice Boltzmann method becomes a complete system. The evolution consists of two processes, one is a collision process which is described by the right hand side of (2.20), i.e.,

f_i^c(n,j) =fi(n,j) +J(f)i(n,j) +gi(n,j). (2.40) This process models the local interaction among particles at a node. Second, the transport process realizes the advection of particles in one time interval,

fi(n+ 1,j+c_i) =f_i^c(n,j). (2.41) At an ordinary node, particles simply move to one of their neighbors with a certain velocity inV. When a nodex_j is next to the boundary and thus some of its neighbors are out of the domain (see figure 2.2),

∂Ω xj

h hqji

cⁱ c^i∗

xji=xj−hqjici

Figure 2.2: Intersection of links and boundary give rise toxji∈∂Ω.

for example x_j −hc_i, then a particular treatment for f_i(n+ 1,x_j) must be introduced. This treatment depends on the geometry of Ω and the predeter-mined boundary settings, for instance, the prescribed average velocity along the boundary. Therefore the transport process differs respectively to the various boundary treatments. Since analysis of the boundary conditions is one of our main goals in this work we describe the frequently used initial conditions, and give a list of existing boundary conditions in computational fluid dynamics with their possible treatments in the lattice Boltzmann framework.