Analysis of Lattice Boltzmann Boundary Conditions

(1)

Analysis of Lattice Boltzmann Boundary Conditions

Dissertation

zur Erlangung des akademischen Grades

des Doktors der Naturwissenschaften (Dr. rer. nat.) an der

Universit¨ at Konstanz

Mathematisch-Naturwissenschaftliche Sektion Fachbereich Mathematik und Statistik

vorgelegt von

M.Sc. Zhaoxia Yang

Tag der m¨undlichen Pr¨ufung : 25. Juli 2007

Referent: Prof. Dr. Michael Junk Referent: Prof. Dr. Vincent Heuveline

Konstanzer Online-Publikations-System (KOPS) URL: http://www.ub.uni-konstanz.de/kops/volltexte/2007/3507/

URN: http://nbn-resolving.de/urn:nbn:de:bsz:352-opus-35079

(2)

(3)

to my parents

(4)

(5)

Preface

This dissertation is the summation of the work at the University of Kaiser- slautern, the University of Saarland and the University of Constance starting from the end of 2001, which is financially supported by the German Research Foundation (Deutsche Forschungsgemeinschaft).

I would like to sincerely express my deepest thanks and gratitude to my su- pervisor Prof. Dr. Michael Junk for providing me the opportunity to achieve this work, and for the guidance and the numerous discussions and suggestions, and for the careful proof reading of this dissertation. This work would not be possible without his constant support and wisdom.

I would also like to express my grateful thanks to Prof. Dr. Vincent Heuveline for refereeing my PhD thesis.

I sincerely thank Prof. Dr. Li-Shi Luo and Prof. Dr. Wen-an Yong for useful discussions and encouragement and sharing useful references. My thanks also go to those people, from whom valuable comments are supplied in ICMMES and different workshops.

I gratefully thank all my friends in Beijing, Shenzhen, Qingdao, Hongkong, Kaiserslautern and Konstanz for their care and help and support. I keep each happy moment with them in my memory forever.

Special thanks to my parents for giving me life and endless love wherever I am.

Thanks to my sister and brothers for their care, understanding, support, and unconditional acceptance to me. I am grateful for that from the deepest of my heart.

(6)

(7)

Abstract

In this dissertation, we investigate a class of standard linear and nonlinear lattice Boltzmann methods from the point of view of mathematical analysis.

First we study the consistency of the lattice Boltzmann method on a bounded domain by means of asymptotic analysis. From the analysis of the lattice Boltz- mann update rule, we find a representation of the lattice Boltzmann solutions in form of truncated regular expansions, which clearly exhibit the relation to solutions of the Navier-Stokes equation. Through the analysis of the initial conditions and the well-known bounce back boundary rule, we demonstrate the general procedure to integrate the boundary analysis process in the whole analysis, and find that our approach can reliably predict the accuracy of the lattice Boltzmann solutions as approximations to Navier-Stokes solutions.

Next, a rigorous convergence proof is achieved for the class of standard linear and nonlinear lattice Boltzmann methods considered in this thesis.

Concentrating on realizations of Dirichlet velocity boundary conditions, we then investigate the consistency of several existing implementations, predict their accuracy, and their advantages and shortcomings. In order to overcome a general drawback of the methods, we construct a class of purely local boundary treatments. All of these methods lead to a second order accurate velocity and a first order accurate pressure. A careful numerical comparison of their properties such as stability, mass conservation and error behavior is presented, as well as a guide for choosing a boundary implementation among the various possibilities.

Regarding Navier-Stokes outflow conditions which are hardly studied in the lattice Boltzmann literature, we deal with three kind of Neumann-type conditions. We have proposed their implementations in the lattice Bolzmann framework, and briefly carry out their consistency analysis. Several numerical results demonstrate the capability of these outflow treatments. For the unsteady benchmark problem like flows around fixed cylinders in an infinitely long channel, the proposed do-nothing and zero normal stress conditions perform very well. For the steady flow, all of the methods produce convincing results.

(8)

(9)

Zussamenfassung

In dieser Arbeit werden Konsistenz, Konvergenz und Randbedingungen für eine Klasse von Gitter-Boltzmann Verfahren behandelt, welche im wesentlichen zur numerischen Lösung strömungsdynamischer Gleichungen wie der Navier- Stokes Gleichung eingesetzt werden. Insbesondere stehen Randbedingungen vom Dirichlet bzw. Neumann Typ im Mittelpunkt.

Die Konsistenzanalyse basiert auf einer asymptotischen Entwicklung der numerischen Gitter-Boltzmann Lösung. Dabei wird die numerische Lösung durch eine abgeschnittene reguläre Entwicklung approximiert, anhand derer sich die Verbindung zur Navier-Stokes Gleichung herstellen läßt. Zunächst wird die Konsistenzanalyse zur Untersuchung von Anfangsbedingungen und der bekan- nten Bounce-back Randbedingung herangezogen und beispielhaft erläutert. Es stellt sich heraus, daß der analytische Zugang eine zuverlässige Vorhersage der Genauigkeit (Konvergenzordnung) des Gitter-Boltzmann Verfahrens ermöglicht.

Die theoretischen Untersuchungen werden durch einen Konvergenzbeweis f¨ur das Gitter-Boltzmannverfahren abgerundet. Dabei werden neben periodischen Randbedingungen auch die klassischen Bounce-back Randbedingungen betrach- tet.

Im anwendungsorientierten Teil werden verschiedene, bereits existierende Um- setzungen der Dirichlet Randbedingung in Bezug auf Konsistenz und sonstige Vor- und Nachteile verglichen. Dies führt zur Konstruktion einer neuen, lokalen (Ein-Knoten) Randbedingung, welche die Schwächen anderer Randbedingun- gen überkommt, ohne wesentliche Vorteile einzubüßen. Konsistenzanalysen und numerische Test zeigen, daß alle Randbedingungen von zweiter Approx- imationsgüte im Geschwindigkeitsfeld sind, während der Druck nur mit erster Genauigkeitsordnung berechnet werden kann. Ergänzend werden Tests zur Sta- bilität und Massenerhaltung der Randbedingungen durchgeführt.

Trotz ihrer praktischen Bedeutung sind Ausflußbedingungen für die Navier- Stokes Gleichung bisher kaum in der Gitter-Boltzmann Literatur diskutiert worden. Zur Umsetzung dieser Neumann-artigen Randbedingungen werden hier drei verschiedene Ansätze verfolgt: die gewöhnliche Neumannbedingung, die Bedingung verschwindender Normalspannung sowie die ’do-nothing’ Bedin- gung. Anhand numerischer Simulationen stationärer und transienter Kanal- strömungen (mit feststehendem Hindernis) werden die Bedingungen erprobt.

Während im stationären Falle alle Randbedingungen zufriedenstellende Resul- tate liefern, überzeugen allein die beiden letzteren bei zeitabhängigen Ström- ungen.

(10)

(11)

Introduction

1.1 Brief history

The lattice Boltzmann method (LBM) originated from the lattice gas automata (LGA) [8, 28, 23, 77, 22, 57, 39, 36, 6, 66, 67, 11, 12] by taking ensemble averages on the related evolution equations. The motivation for the transition from LGA to LBM is to get rid of the noise due to the fluctuation of particle numbers. The lattice Boltzmann method as an independent numerical method is first introduced by McNamara and Zanneti [57] in 1988 to simulate hydrodynamic problems. Later in order to improve the computational efficiency which is mainly pulled down by the cumbersome collision process, the linearized collision operator [36] is preferred. For simplicity, some variations are developed including the well-known single relaxation time BGK model [6] and the MRT (multi-relaxation time) model [17]. Since then the application of the lattice Boltzmann method is simplified and greatly increased.

Without tracing back to its predecessor, the lattice Boltzmann method can also be viewed as a special finite difference scheme for a continuous Boltzmann equation with finite number of velocities by using a small Mach number expansion and by discretizing the phase space and time in a coherent manner.

The detailed derivation can be found in [33, 3]. From this point of view, many numerical techniques for solving PDEs can be used to improve the lattice Boltz- mann method.

Although the lattice Boltzmann method is formulated on the level of particles, its principal application focuses on the macroscopic behavior. Historically the Chapmann-Enskog analysis with a formal multiscale expansion [22] shows that the Euler equations are recovered on the fast convective time scale, and the Navier-Stokes equation from the slower diffusive time scale in the near incompressible limit (due to low Mach number). Because of the kinetic nature, the lattice Boltzmann method has certain features that distinguish it from con- ventional CFD methods to solve the Navier-Stokes equation directly. A great deal of practice has demonstrated the success of the lattice Boltzmann method as well as its variations to meet different applications, particularly to simulate complex fluids such as porous media flow, multiphase flow, multicomponent flow, granular flow, etc..

1

(16)

As for the boundary condition, the bounce back rule [77, 51] is the classical approach in lattice Boltzmann simulations, which is used to simulate the non- slip solid wall and later proved not to be accurate enough [16, 80, 24, 34].

Thereafter much effort has been made to construct more accurate boundary schemes than the bounce back rule. The early boundary treatments appeared as variations of the bounce back rule and were restricted to very simple regular boundary geometries ([71, 62, 63, 40, 54, 14, 82, 81]), for example flat walls.

The velocity and/or density must be known at the node where the boundary scheme is formed. Besides, these methods are more or less dependent on some flow properties too. The so-called hydrodynamic boundary condition in [63]

based on the fixed internal energy being equal to the square of sound speed for the steady-state hydrodynamic field, is such an instance.

Later, more endeavor is dedicated to seeking for methods which deal with curved boundaries. In 1998, Fillipova suggested a purely local method FH ([21, 20]).

However the algorithm involves the factor 1/(τ −1) and is unstable whenτ is close to 1, whereτ is the time relaxation parameter in the collision operator. In [58, 59], Mei et al. suggested an improvement (called MLS) to this method by involving one more neighbor node and thus enlarging the stability region, but without overcoming the inherent drawback of FH. In 2001, Bouzidi et al. [7]

proposed a different link-based method (BFL) by means of linear interpolation with one or two neighbor nodes. Later Ginzburg and d’Humi´eres submitted a more general solution in [25]. Beside the exceptions like MR in [25], the other methods mentioned here can be seen as direct extensions of the bounce back rule which will be demonstrated clearly later.

From the point of view of mathematical analysis, the lattice Boltzmann method is a kinetic relaxation method [10] for the macroscopic equations. The issues of convergence, consistency and stability are as important as for other finite difference methods. A von Neumann stability analysis is carried out in [72, 49]

for linearized lattice Boltzmann methods. In various publications, accuracy and stability are checked numerically [68, 64, 78, 53, 76, 38, 4, 56, 10]. However, rigorous analysis results are very limited until now. For example, Elton [19, 18]

studied the convergence, stability and consistency of lattice Boltzmann methods for the viscous Burger´s equation and advection-diffusion systems. For a class of lattice Boltzmann methods solving for Navier-Stokes equation on a periodic domain, the consistency [42, 45] and stability have also been achieved.

However, the analysis of the lattice Boltzmann method on general bounded domains with general boundary conditions is still under consideration. Besides, most of already matured boundary conditions reflect only the Dirichlet case for velocity at the macroscopic level. Other boundary treatments involving derivatives of the fluid velocity are much less studied.

The aim of this thesis is to fill some of these gaps.

(17)

1.2. SCOPE OF THIS WORK 3

1.2 Scope of this work

1.2.1 Consistency

The classical Chapman Enskog analysis with two time scales, which is usually taken as basis for the analysis of lattice Boltzmann schemes (see for example [22, 5, 12, 13, 31, 33]), leads to a compressible Navier-Stokes equation. The incompressible Navier-Stokes equation is obtained consecutively by another limit process.

In [43], we have successfully applied the regular expansion and multiscale expansion in asymptotic analysis to finite difference methods and meanwhile compared with other classical consistency analysis methods such as the truncation error analysis and the modified equation approach. The merits of the approach are: (i) The asymptotic analysis predicts and analyzes the behavior of finite difference solutions order by order very precisely. Usually the lower order coefficients produce the exact solution of the partial difference equations, the higher order coefficients show the numerical error. (ii) Asymptotic analysis also makes it possible to investigate the behavior of initial layers, boundary layers and oscil- lation, by changing the scale appropriate to the layer thickness. (iii) Boundary conditions can be directly introduced into the analysis process and the overall consistency accuracy can be predicted exactly.

Compared to other finite difference methods for the incompressible Navier- Stokes equation, the lattice Boltzmann method is very particular, since the hydrodynamic macroquantities are connected to the microquantities by means of averages according to kinetic theory. More precisely, from a mathematical point of view, the solution of the Navier-Stokes equation results from an asymptotic singular limit of the lattice Boltzmann solution. We show that the chosen approach can simply and easily provide a consistency analysis of boundary algorithms and simultaneously of the whole scheme. The success and advantages of the asymptotic analysis approach indicates that it is a very appropriate tool for the lattice Boltzmann method on a bounded domain.

In [42], the authors also advocate the regular expansion by illustrating the consistency of a wide class of lattice Boltzmann equation on the whole space or on a periodic domain. A comparison with the Chapman-Enskog expansion is also given at last. Comparatively, the asymptotic analysis relates the numerical solution to the exact solution and the error in a more straightforward manner.

In this work, we continue the asymptotic analysis as in [42] and apply it to the consistency analysis of lattice Boltzmann methods on general bounded domains.

We will see that it plays an important role in the rigorous convergence proof too.

1.2.2 Convergence

Very few results about convergence have been presented so far. For instance, El- ton [18, 19] established the convergence theory for nonlinear convective-diffusive

(18)

lattice Boltzmann methods, but his results apply to schemes with non-linear collision operators satisfying an H-theorem. Moreover, his work is focused on the initial value problem, i.e. the spatial domain is either the whole space or periodic.

The use of truncated regular asymptotic expansions in the proof of convergence can be dated back to [73], and has been successfully applied to the lattice Boltzmann case in [18, 19]. We continue the idea in this work to prove the convergence for a class of standard linear and nonlinear lattice Boltzmann methods based on a rigorous stability result [44].

1.2.3 Boundary conditions

Implementing boundary conditions in the lattice Boltzmann setting is diffi- cult because there is no one-to-one mapping between the variables of the algorithm, the so-called particle distributions, and the predetermined hydrodynamic macroquantities given at the boundary. For example, if Dirichlet boundary conditions are to be implemented, one cannot directly prescribe the velocity at boundary nodes but one has to set the particle distributions in such a way that the average velocity satisfies the required conditions. Typically, the required number of conditions on the kinetic level exceeds the available conditions from the Navier-Stokes problems. This indicates that the kinetic conditions have to be chosen carefully in order to avoid the appearance of extra conditions on the Navier-Stokes level which would render the problem ill posed (leading to an unwanted behavior on the grid scale like boundary layers, oscillations etc.). For other cases like outflow conditions or conditions at free surfaces, the situation is additionally complicated by the fact that the relation between the particle distributions and the quantities specified in the boundary conditions is less evident than in the Dirichlet case.

In this work the Dirichlet and outflow boundary conditions are studied in detail. Our aims are: (i) to analyze the consistency of existing lattice Boltzmann boundary conditions, (ii) to construct new accurate boundary conditions and carry out the related consistency analysis, and (iii) to offer a practical guide for choosing a boundary implementation among all the possibilities.

1.2.4 Publications

This work is supported by the Deutsche Forschungsgemeinschaft (DFG). Several results achieved in this thesis have already been published in collaboration with Prof. Junk. However, the presentation here is generally more detailed with more numerical examples and a uniform use of MRT models in the derivation.

1. Analysis of lattice Boltzmann boundary conditions, Proc. Appl.

Math. Mech. 3(2003)7679.

2. Asymptotic analysis of finite difference methods, Appl. Math. Com- put. 158(2004)267301.

(19)

1.3. SYNOPSIS OF THIS WORK 5 3. Asymptotic Analysis of lattice Boltzmann boundary conditions, J.

Stat. Phys. 121(2005), 335.

4. One-point boundary condition for the lattice Boltzmann method, Phys. Rev. E, 72,1(2005).

5. Convergence of lattice Boltzmann methods for Stokes flows in periodic and bounded domains, Int. J. of Comp. Fluid Dynamics, Vol.

20 No. 6(2006).

6. Convergence of lattice Boltzmann methods for Navier-Stokes flows in periodic and bounded domains, submitted to Numer. Math.

7. Outflow conditions for the lattice Boltzmann method, Progress in Computational Fluid Dynamics, to appear

1.3 Synopsis of this work

Since the analysis of lattice Boltzmann methods together with the various boundary treatments is the main task, the layout of this work is ordered according to the macroscopic boundary conditions.

Inchapter 2, we describe the lattice Boltzmann methods studied in this work.

In particular, the required properties of the collision operators, discrete velocity sets and the related weight functions are discussed. It is stressed that the considered class of lattice Boltzmann methods contains the generally used BGK [6] and MRT [17] models as well asD2Q9,D3Q15,D3Q19 andD3Q27 discrete velocity sets and weight functions [67, 58].

In section 2.3, a feasible initial condition is proposed followed by a list of boundary conditions for well-posed Navier-Stokes problems. To some of them, the corresponding treatments in the lattice Boltzmann framework are also introduced.

The Dirichlet condition and Neumann-type outflow conditions are carefully studied in chapters 5 and 6.

In chapter 3, the main task is to demonstrate the procedure of consistency analysis of lattice Boltzmann methods on general bounded domains by means of asymptotic analysis.

The first step is to obtain the explicit expressions for the coefficientsf_i^(k)in the regular asymptotic expansion of the lattice Boltzmann solution

f_i(n,j) =f_i⁽⁰⁾(t_n,x_j) +hf_i⁽¹⁾(t_n,x_j) +h²f_i⁽²⁾(t_n,x_j) +. . .

at least to fifth order in h. This is done in section 3.1 for both linear and nonlinear lattice Boltzmann methods by reviewing the results of [42]. For linear flows and the polynomial Poiseuille flow, the coefficients f_i^(k) are computed explicitly. Based on these known coefficients we obtain a truncated expansion, which is also a representation of the lattice Boltzmann solutions to some extent.

(20)

This truncated expansion clearly shows a relation between the Navier-Stokes solutions and the lattice Boltzmann solutions, as well as the numerical error terms governed by several Oseen-type equations. The first order error terms can disappear if the corresponding initial and boundary values are zero. The second order error terms, however, usually can not be zero.

Next, the analysis of initial conditions [42, 44] is recalled in section 3.3. Three initial treatments are explained and numerically tested. The most important information is that the first order and second order accurate treatments can be generally implemented. On the contrary, the third order accurate approach is generally not feasible because the required values depend on the initial time derivative of fluid pressure. It is also remarked that the second order accurate initial condition generally leads to a first order accurate pressure and second order accurate velocity.

Last we come to the analysis of boundary conditions which is the core part of this chapter. In section 3.5 we take the bounce back rule as an example to demonstrate the boundary scheme analysis in detail. We insert the regular expansion in the bounce back rule, and do Taylor expansion at a suitable point on the boundary, then collect terms according to the order of h. By imposing the residue of the boundary scheme to be of high order in h, we achieve the accuracy of the boundary scheme as well as of the entire algorithm.

In the case of the bounce back rule we encounter conditions at second order which have a structure that is very typical and common in the residue for other lattice Boltzmann boundary schemes and which cannot be satisfied by smooth functions. This indicates irregular behavior which is actually observed in the error of the velocity (see the following figure for the cases of lattice Boltzmann solutions of the circular flow (Left) in a unit disk and the Taylor vortex flow in a unit square).

−1 −0.5 0 0.5 1

−0.04

−0.02 0 0.02 0.04

0 0.5 1

−3

−2

−1 0 1 2 3x 10⁻³

Figure 1.1: Left: Plot of velocity error (horizontal component (∗) and vertical component (◦)) at t = 0.6 for the circular flow. Right: Plot of horizontal component of velocity error along cutsx= 0.1 (◦),x= 0.3 (∗) for the decaying Taylor vortex flow.

(21)

1.3. SYNOPSIS OF THIS WORK 7 In chapter 4, a rigorous convergence proof for linear and nonlinear lattice Boltzmann methods is given in the case of periodic boundaries as well as the Dirichlet case with bounce back rule at half links.

Summarizing the consistency and convergence results, we have found that, in periodic cases, the regular expansion representation is generally consistent and convergent to the lattice Boltzmann solution in the same order, which is the order of the residue concerned with the initial treatments. This conclusion co- incides with the analytical results in section 3.3 and is verified by the numerical tests in section 3.4. It showed that the convergence order of moments to the Stokes or Navier-Stokes solutions on a periodic domain is generally determined only by the initial error, and is up to two for both the velocity and the pressure of the fluid, provided that the solutions are sufficiently regular. In Dirichlet cases, due to the coarse estimate, we do not achieve as good consistency and convergence results as in the periodic cases. Nevertheless, from the convergence proof of the nonlinear lattice Boltzmann method, we see that, the grid size is required to be suitably small if the Navier-Stokes solution has strong gradient and the Reynolds number is high.

Inchapter 5, the consistency and accuracy of boundary schemes simulating velocity Dirichlet conditions is studied by using the asymptotic analysis developed in chapter 3.

First, we investigate several existing boundary schemes and display their advantages and drawbacks. It turns out that the finite difference techniques, FH, MLS, and BFL can be considered as improvements of the bounce back rule, in the sense that the leading order error terms of the bounce back method are removed, and generally lead to 2nd order accurate velocity and 1st order accurate pressure. Multi-reflection methods (MR) try to use the information at three nodes to get higher accuracy up to 3rd order. However,

• the finite difference technique, MLS, BFL and MR are not local and em- ploy two or three nodes. For nodes where the required neighbors are not available, these methods are no more applicable. Usually the bounce back rule is suggested as a remedy at these points. However we find that the low order error of bounce back rule is transported gradually everywhere into the domain.

• the FH method is local, but leads to unwanted velocity condition at certain nodes.

• the finite difference technique and the MR method show very bad stability in our numerical tests. FH and MLS methods depend on 1/(τ −1) and 1/(τ−2) respectively which eventually leads to instability when τ ≈1 or 2. Only BFL shows satisfying stability.

In order to overcome the shortcomings of the above-mentioned methods while remaining the advantages, a class of boundary treatments POP_θ (θ∈[0,1]) is developed in section 5.3 which has the following characteristics:

(22)

• POPθ are local. POP₀ is explicit and POP₁ is purely implicit.

• the numerical stability of POP_θ is similar to the one of BFL.

• POP_θleads to 2nd order accurate velocity and 1st order accurate pressure.

In addition, in all of our numerical tests POPθ leads tosmaller error and less total mass variations compared to BFL. The following figures are examples to illustrate their behaviors.

0.00183 0.00157 0.0013 0.00104 0.000508

0.2 0.4 0.6 0.8 0.1

0.2 0.3 0.4 0.5 0.6 0.7 0.8

0.9 −8.87e−05

0.000184 0.000162 0.000139 7.06e−05

0.2 0.4 0.6 0.8 0.1

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Figure 1.2: Error contour lines of x-component velocity for the stationary linear flow (3.82) withν = 0.01 and boundary conditions: BFL (left) and POP1(right).

0.5 0.55 0.6

−1 0 1 2 3

x 10⁻⁷

log10(h)

log△˜ρ(t)

0.5 0.55 0.6

0 1 2 3

x 10⁻⁷

log10(h)

log△˜ρ(t)

Figure 1.3: Deviation of the average density from its initial value (△ρ(t)) versus time˜ on a grid of size h = 1/10(⋄), 1/20(◦), 1/30(∇), 1/40(+) and 1/50(∗) for a circular flow. Right: POP¹; Left: BFL.

Inchapter 6, lattice Boltzmann algorithms are constructed to deal with Neumann- type outflow boundary conditions including the Neumann boundary condition (NBC), the Do-nothing condition (DNT) and the zero normal stress condition (ZNS). The related consistency is also presented which shows that all these

(23)

1.3. SYNOPSIS OF THIS WORK 9 treatments lead to 1st order accuracy with respect to the corresponding outflow conditions.

As a benchmark problem, we use the flow around a fixed cylinder in an infinite long channel to evaluate these outflow boundary conditions. One aspect is to check the drag coefficientC_d, the lift coefficient C_l and the pressure difference

△P between the front and end point of the cylinder, which can reflect how the outflow boundary condition influences the inner flow relatively far away from the outflow boundary. Another aspect is to check the behavior of the flow at the outflow boundary itself.

For the steady flow around the cylinder, all these three outflow conditions work very well and produce results forC_d,C_l and△P in the reference interval if the grid is sufficiently fine. Besides, the length of channel does not have a strong effect on these quantities.

However, for the unsteady flow around the cylinder, we find that the length of the channel has obvious effects on the values of Cd,Cl and △P. The NBC outflow condition has a strong impact on the inner flow. On the contrary, ZNS and DNT lead to less variations of the velocity and pressure. Only a slight phase difference occurs when the channel length is varied. Figure 1.4 shows the comparison of the pressure along the center line of the channel. Figures 1.5–1.6 show plots of the fluid streamlines. When the channel is too short, the streamlines are deformed apparently if NBC is employed. On the contrary, ZNS and DNT show a very similar flow structure even in short channels.

0 50 100 150

−1

−0.5 0 0.5 1 1.5 2 2.5

x

0 50 100 150

−1

−0.5 0 0.5 1 1.5 2 2.5

x

0 50 100 150

0 1 2 3 4 5

x

Figure 1.4: The pressure along the center line of the channel for the unsteady flow with ZNS (left), DNT (middle) and NBC (right) outflow conditions respectively. The symbols +,◦and▽stand for the length/width ratio of the channel of 2, 3 and 5. The comparison is restricted to the commonx-region of the three channels.

(24)

50 100 150 50 100 150 50 100 150

Figure 1.5: The truncated streamlines with NBC in channels of different length/width ratio 2 (left), 3 (middle) and 5 (right). A clear channel length dependence is visible.

100 150 100 150

Figure 1.6: The truncated streamlines with ZNS (left) and DNT (right) in channels with length/width ratio 2 (upper row), 3 (middle row) and 5 (lower row). The channel length dependence is much weaker than that for NBC condition.

(25)

Chapter 2

Lattice Boltzmann method

The lattice Boltzmann method has many variational models due to different collision operators, various discrete velocity sets and the involved numerical techniques after realizing that the lattice Boltzmann method is a finite difference approximation of the continuous discrete Boltzmann equation [31, 33, 41], for example lattice Boltzmann model on non-uniform lattice. There are also a lot of generation models of the lattice Boltzmann method to meet various applications such as multiphase flows, subgrid scale modeling, flows involving energy, etc.

Naturally not all kinds of lattice Boltzmann models could be studied in one thesis.

This chapter is dedicated to highlight a standard class of lattice Boltzmann models used in this text, which is supposed to be a numerical solver for the nondimensional incompressible Navier-Stokes equation:

∇ ·u= 0, ∂tu+∇p+ (u· ∇)u=ν∇²u+G, (2.1) with an initially divergence free velocity field

u(0,x) =ψ(x), x∈Ω, (2.2) and a compatible boundary condition

u(t,x) =φ(t,x), x∈∂Ω (2.3)

satisfying

Z

∂Ω

φ·n= 0.

Here Ω∈R^d is the flow domain,u(t,x) is the velocity field of the fluid,p(t, x) represents the flow pressure,ν is the fluid viscosity. G(t,x) stands for the body force.

Before the detailed description some basic notations and definition are given.

2.1 Preliminaries and Notations

In this section we introduce the notation which is used in the later parts of this work.

11

(26)

Throughout this text, the domain Ω denotes a bounded, non-empty, open set in space R^d (d = 2,3). ∂Ω is its boundary and ¯Ω its closure. The lattice Boltzmann method is formulated at the discrete nodes in the domain ¯Ω. At each node there is a finite number of possible velocity directions. For example, the following figure depicts the node distribution in a two dimensional domain with a velocity set of D2Q9 type (see appendix D.1 for details).

0000 1111 0000 1111 0000 1111 0000 1111

0000 1111 0000 1111 0000 1111

0000 1111 0000 1111

0000 1111 0000 1111 0000 1111

Ω ∂Ω

Figure 2.1: Node distribution plot in a 2D case with velocity directions of the D2Q9 model. Fluid nodes in Ω are marked by filled circles, the non-fluid nodes are marked by hollow circles.

The connection between two nodes along some velocity direction is called link.

The nodes are classified into two types, fluid nodes and non-fluid nodes (see the above figure). A fluid node is called an ordinary node if all of its closest neighbors along the links are fluid nodes. Correspondingly, a fluid node called boundary node has at least one non-fluid node as neighbor. In figure 2.1, the black circles are ordinary nodes. The grey nodes are boundary nodes and the white circles are non-fluid nodes. In addition, a link is called boundary link (see dashed line in figure 2.1) if it connects a boundary node and a non-fluid node.

For the variables in our work we adopt the following conventions,

• The bold face symbols such asf,u,v always represent vectors.

• Greek subscripts stand for the spatial coordinates.

• Roman subscripts represent the indices of discrete velocities.

Let V ={c₁, . . . ,c_N} ⊂R^d denote a finite discrete velocity set, and 1 be the vector inR^N with all entries one, we introduce dvectors inR^N based onV,

vα= [c_1α, . . . , cN α]^T (2.4) wherecjα is theαth component of vector cj, and α is taken from 1 tod. Next we defined operatorsVα:R^N →R^N by

Vαf = diag(vα)f = [c_1αf₁, . . . , cN αfN]^T. (2.5)

(27)

2.1. PRELIMINARIES AND NOTATIONS 13 AmongVα, any two of them commute,

VαVβ =VβVα, ∀α, β ∈ {1,2, . . . , d}. (2.6) ObviouslyV_α is a linear mapping on the vector space R^N and V_α1=v_α. A vector operator is hence generated by setting V = [V₁, . . . , V_d]^T, and V : R^N →R^d×N,

Vf = [V₁f, . . . , Vdf]^T, ∀f ∈R^N. (2.7) Then|V|²=V₁²+V₂²+. . .+V_d²is an operator fromR^N toR^N. Let·denote the inner product in the spaceR^d. Observing that∇= [∂₁, . . . , ∂_d]^T is also a vector operator with the same dimension asV, we can thus define a composition by

V· ∇=

d

X

α=1

Vα∂α, (2.8)

so that V· ∇ is an operator from F = C¹(R^d,R^N) to R^N. Since V_α and ∂_α commute,V and ∇also commute with respect to·, i.e., V· ∇=∇ ·V. Iff is an arbitrary function inF, then

(∇ ·V)f = [(c₁· ∇)f₁, . . . ,(c_N · ∇)f_N]^T. (2.9) Moreover, for an arbitrary vector π in R^d, (π·V) is again an operator from R^N to R^N. It turns out to be

(π·V)f =

d

X

α=1

παVαf = [(π·c₁)f₁, . . . ,(π·cN)fN]^T. (2.10) Further leth·,·idenote the inner product on the vector spaceR^N. Iff,g∈R^N are two arbitrary vectors, then the following products will occur frequently in the lattice Boltzmann context,

h1,fi=

N

X

i=1

fi, (2.11)

h1, Vαfi=hvα,fi=

N

X

i=1

ciαfi, (2.12)

hVβg, Vαfi=hg, VβVαfi=

N

X

i=1

ciαciβfigi. (2.13) Since V is a vector operator, we can apply the inner product to each of its component, for example,

hg,Vfi= [hg, V₁fi, . . . ,hg, Vdfi]^T . (2.14) It is easy to prove thathg,Vfi=hVg,fiand in particular,

h1,Vfi= [hv₁,fi, . . . ,hv_d,fi]^T =

N

X

i=1

c_ifi. (2.15)

(28)

In addition if a tensor product between two vectors a and bis defined by (a⊗b)ij = 1

2(aibj+ajbi), (2.16) we can define a matrix operatorV⊗Vand calculate

hg,V⊗Vfi^αβ = (hVαg, Vβfi)_αβ. (2.17) Moreover, a related calculation

h1,V(V· ∇)fi=∇ · h1,V⊗Vfi (2.18) is frequently used in this work. Finally the : − product between matrices A and B,

A:B =X

i,j

AijBij (2.19)

is introduced.

2.2 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êq(f)−f), (2.21) where A is a linear mapping and fêq 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êqand the linear

(29)

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^∗

ˆ

ρ+c⁻²_s uˆ·c_i+c⁻⁴_s

2 (ˆu·c_i)²− c⁻²_s 2 |uˆ|²

(2.24) with respect to the total mass density ˆρ and the average velocity ˆu (or more precisely, the average momentum ¯ρu),ˆ

ˆ ρ=

N

X

i=1

f_i, uˆ =

N

X

i=1

c_if_i. (2.25)

Related to the velocity setV, the weight functionf_i^∗=F_i^eq(1,0), which is also 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

X

i=1

f_i^∗= 1,

N

X

i=1

ciαciβf_i^∗ =c²_sδαβ, (2.27)

N

X

i=1

ciαciβciγciδf_i^∗=c⁴_s(δαβδγδ+δαδδβγ+δαγδβδ). (2.28) 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^∗,

for anyf ∈ F with ˆρand ûdefined by (2.25). The quadratic part is defined by f^Q(f,s) =F^Q(û,w),ˆ F_i^Q(û,w) =ˆ c⁻⁴_s

2

(ˆu⊗w) : (cˆ i⊗ci−c²_sI)f_i^∗ ,

(30)

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(û,w)ˆ i= û⊗w.ˆ (2.35) After giving the equilibrium function f_iêq, 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 observing (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)

(31)

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ⁱ^∗

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 predetermined 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.

(32)

2.3 Initial condition

A feasible initialization proposed in [71] is

fi(0,j) =f_i^eq(1 +c⁻²_s h²p(0,x_j), hψ(xj))−c⁻²_s h²A^†V· ∇(V·ψ)(xj) (2.42) where p(0,x) is the pressure corresponding to the initial velocity fieldψ. It is obtained by solving the Poisson equation

∆p=−∇·(ψ· ∇ψ) +∇·G.

2.4 Boundary condition

Here we list some typical boundary conditions for the velocity field u and pressurep. A similar list can be found in [75] and [54].

1. periodic boundary

Assume the problem is periodic along the coordinate axes with a period π∈R^d, whereπiis the period inith direction. Hence Ω =Qd

i=1[0, πi), the solutions are completely characterized by their values on the periodicity cell. This implies that, for example, the velocity field satisfies

u(t, x₁, . . . ,0, . . . , xd) =u(t, x₁, . . . , πi, . . . , xd),

and similar conditions follow for the pressure and their derivatives. In this work we take πi = 1 without loss of generality.

2. symmetric boundary

The symmetric boundary is usually realized by setting the component of velocity normal to the boundary to be zero,

u(t,x)·n(x) = 0, x∈∂Ω,

Here n(x) denotes the unit outer normal vector at the pointx∈∂Ω.

3. Dirichlet boundary

In many applications, the fluid moves in a bounded domain. We treat the behavior of flow along the boundary by the Dirichlet boundary condition, i.e., the velocity components are prescribed.

For most solid surfaces which are impermeable to fluid, the fluid sticks to their surfaces. Hence, there is no slip and no penetration, and the fluid particles on the wall move with the velocity of the wall. Namely, the velocity component normal to the wall is set to be zero, the tangential velocity is identical to the velocity of the wall. In the special case of a stationary wall, the velocity is zero.

(33)

2.4.1 Periodic boundary 19 4. outflow boundary

Outflow conditions are often used at artificial boundaries to simulate a larger flow domain. One possibility is to use the so-called zero Neumann condition ^∂u_∂n(t,x) = 0, x ∈ ∂Ω. There are more outflow conditions mentioned later in chapter 6.

5. pressure boundary

This kind of boundary condition allows to stimulate the flow experiments where a flow is induced by applying a pressure difference. From the mathematical analysis in [35], it is pointed out that the incompressible Navier- Stokes equation (2.1) with the prescribed average pressure condition

p(t,x)n−ν∂_nu(t,x) = ¯P_Γn, x∈Γ (2.43) is a well-posed problem. Where ¯P_Γ is the known average pressure on the boundary segment Γ.

6. free surface

Free surfaces occur at the interface between two fluids, for example water and oil or air. Such interfaces require a kinematic condition which relates the motion of the free interface to the fluid velocities at the free surface and a dynamic condition which is concerned with the force balance at the free surface.

In the following we address the implementations for two of the above mentioned boundary conditions in the lattice Boltzmann setting, which will be analyzed in the next chapter. Additional implementations will be discussed in chapter 5 and chapter 6.

2.4.1 Periodic boundary

In case that a flow has a periodic structure in the whole space, the computational domain Ω can be restricted to a region of the size of one period. To be specific, let us assume that the periodicity cell is given by Ω = [0,1)^d and that h= 1/m for some m >1. Defining the modulo addition (j+mk)i = (ji+ki) modmfor vectorsj,k∈Z^d, we can formulate the advection step in the periodic domain uniformly by

fi(n+ 1,j+mci) =f_i^c(n,j), i= 1, . . . , N, x_j ∈Ω, (2.44) no matter ifxj is an ordinary node or a boundary node.

2.4.2 Bounce back boundary

In order to approximate the Dirichlet boundary condition (2.3) on a general bounded domain, the bounce back rule is applied at every boundary nodexj to

(34)

the components off which belong to incoming velocities ci. The bounce back rule has the form

fi(n+ 1,j) =f_i^c∗(n,j) + 2hc⁻²_s f_i^∗φ(tn,xji)·ci, (2.45) where x_ji is the intersection point of the link along c_i and the boundary ∂Ω.

The remaining components are treated with the usual update rule (2.41).

(35)

Chapter 3

Consistency of LBM on a general bounded domain

The consistency analysis of lattice Boltzmann methods (including those models in chapter 2) on a periodic domain has been intensively studied in the last few years [22, 31, 33, 42]. The advantage of applying asymptotic analysis to the lattice Boltzmann method is given in [42] where the periodic case is discussed.

In this chapter, we intend to extend the consistency analysis to the lattice Boltzmann methods introduced in chapter 2 on a general bounded domain by means of asymptotic analysis.

As explained in chapter 2, a typical lattice Boltzmann method consists of two parts. One part is the usual update rule which holds at the majority of the nodes in Ω. The other part is the special treatment at the boundary nodes.

Correspondingly, the consistency analysis splits into two parts as well, the analysis of the interior algorithm and the analysis of the boundary condition.

The asymptotic expansion is taken as the regular expansion suggested in [42]:

fi(n,j) =f_i⁽⁰⁾(tn,xj) +hf_i⁽¹⁾(tn,xj) +h²f_i⁽²⁾(tn,xj) +. . . (3.1) with smooth andh independent functions f_i^(k).

According to [43], the analysis of the update rule yields explicit expressions of the coefficients in the asymptotic expansion, which are determined by a series of partial differential equations (among them, the one for the coefficients in the lowest relevant order is usually the target problem). We insert the expansion (3.1) into the update rule (2.20), do Taylor expansion at node (t_n,x_j), and define f_i^(s) (s ≤k) in such a way that f_i⁽⁰⁾+. . .+h^kf_i^(k) satisfies (2.20) with high accuracy. Then we substitute the asymptotic expansion into the boundary scheme, do Taylor expansion at suitable boundary points, and check which boundary condition leads to a residue of high order inh. In this way, consistency information about the full algorithm is achieved.

This chapter gives an outline of a typical analysis:

1. Section 3.1 derives expressions for the coefficients in (3.1) from the update rule of the algorithm, the relation between the coefficients f_i^(k) and the

21

(36)

solution of the incompressible Navier-Stokes equation, and the partial differential equations for the higher order error terms. The approach is based on the assumption that the numerical solution can be described with smooth functions up to high orders in h. Since the update rule is the same as in the periodic case, the investigation parallels the one in [42].

2. In section 3.3, while summarizing the consistency analysis of the existing initial conditions in [42, 44], a practically feasible initial condition with possible higher accuracy is proposed so that we can concentrate our attention particularly on the boundary conditions in the later parts.

3. In section 3.5, taking the well known bounce back rule as example, a careful description of the boundary consistency analysis is illustrated.

For the lattice Boltzmann D2Q9 model with BGK [6] type collision operator, a similar result has already been presented in our article [44].

3.1 Formal asymptotic expansion

The aim of this section is to compute the coefficients in (3.1). First we see that the asymptotic expansion of the moments ˆρ and ˆu are correspondingly defined by

ˆ

ρ=ρ₀+hρ₁+h²ρ₂+. . . , (3.2) ˆ

u=u₀+hu₁+h²u₂+. . . , (3.3) where the kth order coefficient of the moments

ρk=X

i

f_i^(k), uk=X

i

f_i^(k)ci

are derived from f^(k). From the knowledge of f^(k) we expect to derive the relation between the numerical values ˆρ,uˆand the solution of theincompressible Navier-Stokes equation.

To begin with we insert (3.1) into the left hand side of the lattice Boltzmann algorithm (2.20) and find expressions of the form

f_i^(k)(tn+1,x_j+c_i) =f_i^(k)(tn+h²,x_j+hci).

Since the functions f_i^(k) are assumed to be smooth, we can perform a Taylor expansion around the point (tn,x_j). This process yields

f_i^(k)(tn+h²,x_j+hci) =f_i^(k)+h(ci· ∇)f_i^(k)

+h²(∂t+ (ci· ∇)²/2)f_i^(k)+h³(ci· ∇)(∂t+ (ci· ∇)²/6)f_i^(k)+. . . where the right hand side is evaluated at (t_n,x_j). Generalizing this expansion to arbitrary orders, we formally obtain an infinite series

f_i^(k)(tn+h²,xj+hci) =f_i^(k)(tn,xj) +

∞

X

r=1

h^rDr(∂t,ci· ∇)f_i^(k)(tn,xj) (3.4)

Analysis of Lattice Boltzmann Boundary Conditions