半古典晶格波滋曼方法

(1)

1

國立臺灣大學工學院應用力學研究所博士論文

Department of Applied Mechanics College of Engineering National Taiwan University

Doctoral Dissertation

半古典晶格波滋曼方法

Semiclassical Lattice Boltzmann Method

洪立昕 Li-Hsin Hung

指導教授：楊照彥博士

Advisor: Jaw-Yen Yang, Ph.D.

中華民國 99 年 10 月 October, 2010

(2)

終於走到這一步了，四年多前的熱情慢慢累積成一行行的程式、

一篇篇的文章，也差不多消耗殆盡了。很高興也很欣慰能在不算太久的時間內走出台大大門。最感謝的是父母的支持，無論是精神上還是經濟上，都是給我最多幫助的。妹妹來自美國的鼓勵也幫助我不少，

尤其妳比我更累的研究生活，更是鞭策我繼續向上的動力。謝謝珮君和冠孜，妳們的陪伴帶給我豐富而精采的精神生活，方能戰勝這漫長而艱辛的一役。

感謝指導教授楊照彥博士四年多來的提攜，教授不論是嚴謹的研究態度還是充沛的精力跑遍各國都另人景仰和欽羨。

論文定稿之際，承蒙本校許文翰教授、黃美嬌教授，清華大學林昭安教授，中正大學何正榮教授，成功大學陳朝光教授、黃吉川教授以及楊玉姿教授的指導，讓我看到許多不同思考層面，更進一步精進論文內容，令我受益匪淺。

實驗室待了許久，相處過的學長學弟不可勝數，每個人都是我永

遠的好友。特別感謝育炘帶我進入這全新的研究領域、澤揚在博後階

段和我的討論，還有聖鑫及耀天在部分研究內容上的幫忙，讓我能夠

發更多文章，順利畢業。最後，感謝許多我沒提及但仍愛我的人，少

了你(妳)們，這一切也難以成真，謝謝。

(3)

摘要

本論文提出了半古典晶格波滋曼法，此方法是在波滋曼法的基礎下，使用量子統計(Bose-Einstein Statistics 和 Fermi-Dirac Statistics)取代古典統計

(Maxwell-Boltzmann Statistics)後展開所得的方法。此方法可驗證在古典極限(Classical Limit)下可回復原本古典晶格波滋曼法，亦可由數值驗證得知古典極限下可得傳統晶格波滋曼法之結果。本文以此一新提出之模擬方法研究量子氣體運動問題，模擬之問題忽略粒子間交互作用，但因為使用量子統計之故，粒子在量子統計下的特性如庖立不相容原理(Pauli Exclusion Principle)、巨觀上傳輸係數之修正等等均有考慮在內。本文探討一維量子氣體之震波管問題、

二維圓柱流問題、二維微管問題以及三維頂蓋流問題作為驗證半古典晶格波滋曼法之方式，模擬結果指出了量子統計和古典統計結果的主要差異。此外，本文亦提出了另一種雙分布函數熱晶格波滋曼法，可得出一般雙分布函數熱晶格波滋曼法所得結果，也可作為未來更進一步推廣半古典晶格波滋曼法之方向。

關鍵詞:晶格波滋曼法、量子統計、波滋曼方程式

(4)

Unlike describing the physical phenomenon in coordinate or momentum spaces in quantum mechanics, semiclassical Boltzmann equation treats the system in phase space, and it is much easier to describe the dynamics of quantum gases. In this thesis, a class of semiclassical lattice Boltzmann methods is developed for solving quantum hydrodynamics and beyond.

The present method is directly derived by projecting the Uehling-Uhlenbeck Boltzmann-BGK equations onto the tensor Hermite polynomials following Grad's moment expansion method. The intrinsic discrete nodes of the Gauss-Hermite quadrature provide the natural lattice velocities for the semiclassical lattice Boltzmann method. Formulations for the second-order and third order expansion of the semiclassical equilibrium distribution functions are derived and their corresponding hydrodynamics are studied. Gases of particles of arbitrary statistics can be considered. Simulations of one-dimensional compressible gas flow by using D1Q5 lattices, two dimensional microchannel flow, two dimensional flow over cylinder by using D2Q9 lattices and three dimensional lid driven cavity flow by using D3Q19 lattices are provided for validating this method. It is shown that the classical flow patterns such like vortex and vortices shedding in flow over cylinder simulations, temperature and pressure contours together with streamline patterns could be produced from the present method in classical limit. The results also indicate the distinct characteristics of the effects of quantum statistics when they are compared with fluid phenomena in classical statistics.

Keywords: Lattice Boltzmann Method, Semiclassical, Quantum.

(5)

List of Tables

3.1 D1Q5 quadratures and weights . . . 41 3.2 D2Q9 quadratures and weights . . . 41 3.3 D3Q19 quadratures and weights . . . 41 5.1 Comparisons of drag coecient of dierent boundary condi-

tions in LB method . . . 76 5.2 Comparisons of the average Nusselt number and the maximum

velocity components across the cavity center. The data in parentheses are the locations of the maxima. . . 76

(10)

(11)

List of Figures

2.1 Bose and Fermi functions in dierent fugacity z. (a) Bose Function g1.5^Bose(z), g2.5^Bose(z), and ^g_g^Bose^2.5^Bose^(z)

1.5 (z), (b) Fermi Function g_1.5^{F ermi}(z),g^{F ermi}2.5 (z), and ^g_g^{F ermi}^2.5^{F ermi}^(z)

1.5 (z), (c) Comparing Bose ^g_g^Bose^2.5^Bose^(z)

1.5 (z), and Fermi Function ^g_g^{F ermi}^2.5F ermi^(z)

1.5 (z) . . . 19

3.1 Comparisons of streamlines of uniform ow over a circular cylinder in a BE quantum gas in dierent expansion order with z = 0.2 and Re∞ = 40. (a) N = 3, (b) N = 2. . . . 40

3.2 D1Q5 lattice structure . . . 41

4.1 Layout of the curved wall boundary in the regularly spaced lattices. . . 57

4.2 Periodic boundary conditions. . . 57

4.3 Layout of the curved wall boundary in the regularly spaced lattices. . . 58

4.4 Bounce back boundary conditions (a) q < 1/2, (b) q ≥ 1/2. . . 59

4.5 The inuences of dierent accommodation coecients on velocity prole based on two dierent boundary conditions. (a) Bounce-back specular reection boundary condition, (b) Dis- crete maxwellian boundary condition. . . 60

4.6 The inuences of dierent boundary conditions on velocity pro- le based on the same accommodation r = 0.8. (a) Kn=0.5, (b) Kn=2. . . 61

(12)

4.7 The velocity prole of dierent Knudsen number based on the same accommodation coecient r = 0.8 for two dierent bound- ary conditions. (a) Bounce-back specular reection boundary condition, (b) Discrete maxwellian boundary condition. . . 62

5.1 Convergence of solution with rened grids, the cells number is 100, 200, and 400. (a) Velocity, (b) Density, (a) Pressure, (d) Temperature. . . 78 5.2 Shock tube ow in a quantum gas. The eect due to dier-

ent particle statistics; Bose-Einstein: BE, Fermi-Dirac: FD, Maxwell-Boltzmann: MB. (a) Density, (b) Velocity, (c) Pres- sure, (d) Temperature. . . 80 5.3 Shock tube ow in a quantum gas. The eect in high tem-

perature limit; Bose-Einstein: BE, Fermi-Dirac: FD, Maxwell- Boltzmann: MB. (a) Density, (b) Velocity, (c) Pressure, (d) Temperature. . . 82 5.4 Velocity proles in a channel ow of gases of arbitrary statistics

(gas with z = 0.2). (a) Kn=0.05, (b) Kn=0.5, (c) Kn=1. (d) Normalized mass ux in a channel ow as a function of Kn number. . . 83

5.5 Drag and lift coecients in classical limit. (a) Re = 100, (b) Re = 500. . . 84 5.6 Drag and lift coecients in dierent boundary conditions (im-

mersed boundary condition and bounce back boundary condi- tion). (a) Re = 100, (b) Re = 200. . . . 85

(13)

List of Figures ix

5.7 Streamlines of uniform ow over a circular cylinder in a quan- tum gas with z = 0.2 and Re_∞ = 20. (a) BE gas, (b) MB gas, (c) FD gas. . . 86 5.8 Streamlines of uniform ow over a circular cylinder in a quan-

tum gas with z = 0.2 and Re_∞ = 40. (a) BE gas, (b) MB gas, (c) FD gas. . . 87 5.9 Streamlines of uniform ow over a circular cylinder in a quan-

tum gas with z = 0.2 and Re∞= 200. (a) BE gas, (b) MB gas, (c) FD gas. . . 88 5.10 Dierent contours and streamlines of uniform ow over a circu-

lar cylinder in a quantum gas in a classical limit with Re∞ = 200, steps=94200. (a) X-direction velocity(ux), (b) y-direction velocity(uy), (c) fugacity(z), (d) density(ρ), (e) temperature(θ), (f) streamlines. . . 89 5.11 Streamlines for dierent Rayleigh numbers of natural convec-

tion. (a) Ra = 10³, (b) Ra = 10⁴, (c) Ra = 10⁵. . . 90 5.12 Isotherm lines for dierent Rayleigh numbers of natural con-

vection. (a) Ra = 10³, (b) Ra = 10⁴, (c) Ra = 10⁵. . . 90 5.13 Streamlines of Rayleigh-Benard convection for dierent Rayleigh

numbers. (a) Ra = 4000, (b) Ra = 10000, (c) Ra = 50000. . . 91 5.14 Isotherms of Rayleigh-Benard convection for dierent Rayleigh

numbers. (a) Ra = 4000, (b) Ra = 10000, (c) Ra = 50000. . . 92 5.15 Grid renement of three dimensional lid driven simulation, Re =

400 for two dierent centerlines. (a) (x,0,0), (b) (0,0,z). . . 93 5.16 Grid renement of three dimensional lid driven simulation, Re =

1000 for two dierent centerlines. (a) (x,0,0), (b) (0,0,z). . . . 93

(14)

5.17 Streamlines of MB statistics in three mid-planes, Re = 100. (a) xy planes, (b) zy planes, (c) xz planes. . . 94 5.18 Streamlines of MB statistics in three mid-planes, Re = 400. (a)

xy planes, (b) yz planes, (c) xz planes. . . 94 5.19 Streamlines of MB statistics in three mid-planes, Re = 1000.

(a) xy planes, (b) yz planes, (c) xz planes. . . 94 5.20 Centerlines velocity of classical and semiclassical model. (a)

Classical, (b) semiclassical . . . 95 5.21 Velocity vectors of BE statistics in three mid-planes, Re = 100,

z = 0.2. (a) xy planes, (b) yz planes, (c) xz planes. . . 96 5.22 Pressure contours of BE statistics in three mid-planes, Re =

100, z = 0.2. (a) xy planes, (b) yz planes, (c) xz planes. . . . . 96 5.23 Streamlines patters of BE statistics in three mid-planes, Re =

100, z = 0.2. (a) xy planes, (b) yz planes, (c) xz planes. . . . . 96 5.24 Velocity vectors of FD statistics in three mid-planes, Re = 100,

z = 0.2. (a) xy planes, (b) yz planes, (c) xz planes. . . 97 5.25 Pressure contours of FD statistics in three mid-planes, Re =

100, z = 0.2. (a) xy planes, (b) yz planes, (c) xz planes. . . . . 97 5.26 Streamlines patters of FD statistics in three mid-planes, Re =

400, z = 0.2. (a) xy planes, (b) yz planes, (c) xz planes. . . . . 98 5.28 Streamlines patters of BE statistics in three mid-planes, Re =

1000, z = 0.2. (a) xy planes, (b) yz planes, (c) xz planes. . . . 98

(15)

List of Figures 1

5.30 Streamlines patters of BE statistics in three mid-planes, Re = 1000, z = 0.2. (a) xy planes, (b) yz planes, (c) xz planes. . . . 99 5.31 Centerlines velocity of dierent quantum statistics, z = 0.1 and

z = 0.2 . . . 100

(16)

(17)

Chapter 1

Introduction

1.1 Overview

After more than twenty years developments on lattice Boltzmann (LB) method since its introduction [1][2], the LB method is not just the extension of its ancestor lattice-gas automaton [3], and even not just a hydrodynamical equations solver. It has been another approach analyzing Boltzmann equation in compared with Chapman-Enskog expansion [4] and Grad moment method [5].

LB method is based on the kinetic equations for simulating uid ow, see [6][7][8]. Over the past two decades, signicant advances in the development of the LB method [2][9][10][11] based on classical Boltzmann equations with the relaxation time approximation of Bhatnagar, Gross and Krook (BGK) [12] have been achieved. The LB method has demonstrated its ability to simulate hydrodynamic systems, magnetohydrodynamic systems, multi-phase

(18)

and multi-component uids, multi-component ow through porous media, and complex uid systems, see [13]. The LB equations can also be directly derived in a priori manner from the continuous Boltzmann equations [14][15][16]. Most of the classical LB methods are accurate up to the second order, i.e., Navier- Stokes hydrodynamics and have not been extended beyond the level of the Navier-Stokes hydrodynamics. A systematical method [17][16] was proposed for kinetic theory representation of hydrodynamics beyond the Navier-Stokes equations using Grad's moment expansion method [5][18]. The use of Grad's moment expansion method in other kinetic equations such as quantum kinetic equations and Enskog equations can be found in [19][20].

Despite their great success, however, most of the existing LB methods are limited to hydrodynamics of classical particles. Modern development in nanoscale transport requires carriers of particles of arbitrary statistics, e.g., phonon Boltzmann transport in nanocomposite and carrier transport in semiconductors. The extension and generalization of the successful classical LB method to quantum LB method for quantum particles is desirable. Analogous to the classical Boltzmann equations, a semiclassical Boltzmann equations for transport phenomenon in quantum gases has been developed by Uehling and Uhlenbeck [21][22]. Following the work of Uehling and Uhlenbeck based on the Chapman-Enskog procedure [4], the hydrodynamic equations of a trapped dilute Bose gas with damping have been derived [23]. In [19], the quantum Grad expansion using tensor Hermite polynomials has been applied to obtain the non-equilibrium density matrix which reduces to the classical Grad moment expansion if the gas obeys the Boltzmann statistics. The full Boltzmann equations is mathematically dicult to handle due to the collision integral in dierent types of collisions. To avoid the complexity of the the collision term, the relaxation time model originally proposed by BGK model for the clas-

(19)

1.2. Conventional LB method 3

sical non-relativistic neutral and charged gases has been widely used. Also, BGK-type relaxation time models to capture the essential properties of carrier scattering mechanisms can be similarly devised for the UUB equations for var- ious carriers and have been widely used in carrier transports [24]. Recently, kinetic numerical methods for ideal quantum gas dynamics based on Bose- Einstein(BE) and Fermi-Dirac(FD) statistics have been presented [25][26]. A gas-kinetic method for the semiclassical Boltzmann-BGK equations for non- equilibrium transport has been devised [27].

In this thesis, a new semiclassical LB method for the UUB-BGK equations based on Grad's moment expansion method by projecting the UUB-BGK equations onto Hermite polynomial basis has been derived. The relations between the relaxation time, viscosity and thermal conductivity are obtained by applying the Chapman-Enskog method [4] to the UUB-BGK equations for providing the basis for determining relaxation time used in the present semiclassical LB method. Hydrodynamics based on moments up to second and third order expansions are presented. Computational examples to illustrate the methods are given and the eects due to quantum statistics are delineated.

In the rest of this dissertation, the conventional LB method and quantum LB method will be introduced, then comes the basics of semiclassical LB method which was rst proposed in [28].

1.2 Conventional LB method

This section will describe the basic algorithm of the LB method. First, an overview of the historical development of LB method including recently developments and reviews will be given. Next, the method itself will be described.

Finally, a derivation of the necessary equations and the lattice structures will

(20)

be discussed. The conventional LB method originated from its predecessor, the lattice gas cellular automata models [3] has become a competitive numerical tool for simulating uid ows over a wide range of complex physical problems [2] [29] [30] [31] [32] [33]. The theoretical background of LB method is the kinetic theory and Boltzmann equation, which are connected with the macroscopic Navier-Stokes equation by the Chapman-Enskog expansion. The LB method can be regarded as a simplied kinetic scheme by using a nite set of discrete velocities and a simplied collision integral. Both the algorithm and the boundary conditions are easy to implement in LB method. As the LB method computes macroscopic behavior, such as the motion of a uid, with equations describing microscopic scales, it operates on a mesoscopic level in between those two extremes. When compared to the traditional computational uid dynamics techniques, the advantages of LB method are mostly on its clear physics and simple algorithm. While conventional solvers directly discretize the Navier-Stokes equations, the LB method is essentially a rst order explicit discretization of the Boltzmann equation in a discrete phase-space. It can also be shown, that the LB method approximates the Navier-Stokes equations with good accuracy. A good review of LB method can be found, e.g., in [6][7][8]. In LB method, the simulation region is divided into a cartesian grid of cells, each of which only interacts with cells in its direct neighbor- hood. The LB method consists of two steps, the stream step, and the collide step. These are usually combined with no-slip boundary conditions for the domain boundaries or obstacles. The simplicity of the algorithm is especially evident when implementing it. Using a LB method, the particle movement is restricted to a limited number of directions. As shown in Fig. 3.3 and Fig. 3.4, a two-dimensional model with 9 velocities (commonly denoted as D2Q9), and three-dimensional model with 19 velocities (commonly denoted

(21)

1.3. Basics of Semiclassical LB method 5

as D3Q19) are provided.

1.3 Basics of Semiclassical LB method

Although Boltzmann Equation has been successfully applied to many eld like dilute gas dynamics, multi-scale simulations, however, in recent years, micro and nano technology has been emerged quickly, and the transport phenomena in semiconductors at low temperatures is very important. There have been a successful theory in statistical mechanics which can predict the transport coecients like shear viscosity and thermal conductivity of ideal quantum

uid such like electrons in the metal. The questions arise of whether quantum systems like that can be described similar to the one developed for the classical counterpart. When solving these kinds of problems, classical Boltzmann Equation is not enough and it require quantum mechanical treatments. In quantum mechanics identical particles are absolutely indistinguishable from one another and N-particle system can be described by a wave function with permutation symmetry. In nature, it is found that particles with antisymmet- ric wave functions are called fermions which obey FD statistics and particles with symmetric wave functions are called bosons which obey BE statistics.

The statistical properties of fermion and boson systems are profoundly dierent at low temperature. However, in the classical limit, both quantum distributions reduce to the Maxwell-Boltzmann(MB) distribution. Boltzmann equation describes the dynamic behavior of ideal gas by a single-particle distribution function.

In general, there are three strategies to take for statistically treating a quantum system [34]. One is to use a kinetic equation governing the density matrix, another one is to use a kinetic equation with the Wigner distribution

(22)

function, the last one, which is also the method used in this thesis, is to assume a semiclassical kinetic equation such as the UUB [21][22] kinetic equation as a generalization of the classical one. In UUB equation, the collision term of the Boltzmann equation is rewritten with the quantum particle scattering form.

The equilibrium distribution to semiclassical Boltzmann equation is the BE or FD distributions. The particles obey BE distribution are called bosons and FD distribution for fermions, no third category has yet been found.

Considering of the recently successful developing on LB method and well deriving semiclassical Boltzmann equation, it is natural to extend the conventional LB method to the semiclassical LB method for dealing with dierent quantum statistics. Although UUB equation introduced above has considered the quantum statistics and would reveal quantum eects in specic problems.

However, dealing with UUB equation is still a challenge. In this dissertation, we proposed a new semiclassical LB method which extends well known LB method to solve semiclassical Boltzmann equation with BGK approximation.

In [28], a new semiclassical LB-BGK method had been developed, and this method would describe quantum systems in dierent approaches. The following works about rareed channel ow [35] and axisymmetric ow of quantum statistics [36] present dierent applications of this new method. The idea of extending the conventional LB method to semiclassical LB method is to adopt the Grad's moment method to nd solutions to semiclassical Boltz- mann equation by expanding f(x, ζ, t) in terms of Hermite polynomials. And in simulations, the N-th nite order truncated distribution function f^N was considered. the details will be shown in chapter 3.

It is worth mentioning here the dierences between the present method and quantum LB method which is proposed by Succi and Benzi [37][38]. In the quantum LB method, Dirac equation which is the most general equation de-

(23)

1.3. Basics of Semiclassical LB method 7

scribing single particle motion in compliance with quantum theory and special relativity is solved by LB method. In [37], the procedure builds on a formal analogy between the Dirac equation and a special discrete kinetic equation known as LB method, it was then shown that the non-relativistic Schr¨odinger equation ensues from the Dirac equation under an adiabatic assumption that is formally similar to the one which takes the Boltzmann equation to the Navier-Stokes equations in kinetic theory. In [38], it was further shown that by a proper resort to operator splitting methods, the Dirac equation can be in- tegrated as a sequence of three one dimensional LB equation evolving complex valued distribution function. In these works, LB method with complex distribution function is treated as a numerical tool for solving a complex equation.

By using multiscale technique and the Chapman-Enskog expansion on complex variables, the complex partial dierential equations could be recovered.

Recent years, this procedure is applied for solving many dierent equations, for examples, the one-dimensional nonlinear Dirac equation[39], the nonlinear convection-diusion equations in [40] and the complex Ginzburg-Landau equation in [41]. It should be also noticed that several quantum lattice gas cellular automata methods [37][42][43][44] have been recently presented which applying and extending the concept of classical lattice gas cellular automata models to treat the time evolution of wave functions for spinning particles and the Schrödiger equation or the Dirac equation directly. For a more detailed review, see [45]. However, present semiclassical LB method associated with the works presented in [25][26][27] are based on the semiclassical kinetic de- scription. i.e., the particle motion (velocity or momentum) and position are treated in classical mechanics manner while the particles can be of quantum statistics. The procedure and physical meanings are much dierent.

(24)

1.4 Contents of the Dissertation

This Dissertation is organized as following.

Chapter 1 gives an overview of this dissertation, briey describes the basic ideas of semiclassical LB method which combines the elements of semiclassical kinetic theory and LB method. History and recent developments of LB method are also given.

Chapter 2 describes the semiclassical kinetic theory. First, some important concepts of quantum kinetic theory are introduced, then comes the semiclassical Boltzmann equation which has been developed for several years since UUB equation [21][22]. Following the previous works, the corresponding semiclassical hydrodynamic equations are developed and compared with the classical hydrodynamic equations.

In chapter 3, the semiclassical LB method is derived based on expanding the BE and FD equilibrium distribution function onto the Hermite polynomials. The analysis of semiclassical LB method is also given. In chapter 4, the initial condition and the boundary conditions of the semiclassical LB method based on conventional LB method are introduced. Since LB method is a mature uid simulating method, extending the boundary conditions from conventional LB method to semiclassical LB method is straightforward.

The numerical validation and examples are given in chapter 5, including one dimensional shock tube simulation, two dimensional ow over cylinder, natural convection ow, and three dimensional lid-driven cavity simulation.

In chapter 6 the conclusions and future works are given. In Appendix 1, the Chapman-Enskog analysis of semiclassical LB method is given. Finally, a new thermal LB model is proposed for the future development of semiclassical LB method.

(25)

Chapter 2

Semiclassical Kinetic Theory

2.1 Overview

In quantum mechanics, the term semiclassical has dierent meanings and all refer to some approximations or situations that combine quantum and classical properties, for example, Wentzel-Kramers-Brillouin(WKB) approximation or UUB equation. In this chapter, from the basic concepts to quantum gases, the collision term of semiclassical Boltzmann equation is derived. Then, the approximation method, BGK method is introduced, and the equilibrium distribution function is proposed based on Boltzmann H theory. Finally, semiclassical hydrodynamic equations are derived from moment integration of the semiclassical Boltzmann equation with BGK approximation and the results are compared with their classical counterparts.

(26)

2.2 Introduction to Quantum Gases

Comparing with classical mechanics statistics and quantum mechanics statistics, distinguishability is a very important concept. In classical mechanics statistics, particles in a system are distinguishable even if they are identical. It is possible to label particles and track their phase space trajectories with cer- tainty. Conversely, in quantum mechanics statistics the precise descriptions of particles motions are limited by Heisenberg uncertainty relation where the si- multaneous measurements on momentum and position are not allowed. Under these fuzzy descriptions, the particles become untrackable and the property of indistinguishability is inherent. The distinguishable property makes the conguration dierent and the corresponding microstate distinct as particle in- terchanges. The number of microstates are usually very large and microstates will continually change as time goes on. The time-average behavior of the system will be identical to the average behavior of all sort microstates corresponding to the macrostates at time t. The average behaviors on all the collection of microstates in any systems are called ensembles of the systems.

In statistical theory, there are three dierent statistics which are microcanoni- cal ensemble, canonical ensemble and grand canonical ensemble. In canonical ensemble the macrostates are characterized by (N, V, T ), and the partition function is written as QN(V, T ) =∑

Ee^−βE, wherein β = 1/kbT. In the ideal quantum gas system, the energy eigenvalues E are composed of single particle states in occupation set {nε}, ε is the single-particle energy, and E can be expressed in terms of ε as E = ∑

εn_εε. The values of the numbers nε must satisfy the condition

∑

ε

n_ε = N (2.1)

According to the discussions of distinguishable property given above, we

(27)

2.2. Introduction to Quantum Gases 11

can write down the statistical weight factor for FD, BE and MB statistics as:

WF D{nε} = { 1 if all ε = 0 or 1

0 otherwise }, (2.2)

W_BE{nε} = 1, (2.3)

W_{M B} =∏

ε

1

n_ε!. (2.4)

We can further write the partition function of a N particle system Q_N(V, T ) =

∑′ {nε}

W{nε}e^−β^∑^εⁿ^ε^ε (2.5)

∑_′

means the summation over all distribution sets that satises (2.1). To work out the partition function within ∑_′

is relatively complicated since the summation is under some constraints. Now substituting (2.3) and (2.2) into (2.5), one has

Q_N(V, T ) =

∑′ {nε}

e^−β^∑^εⁿ^ε^ε.

The dierences between BE and FD statistics arises from the values that the number nε can take. However, the restriction (2.1) makes the canonical ensemble cumbersome here. The grand partition function seems to provide an ecient way

ϑ(z, V, T )≡

∑∞ N =0

z^NQ_N(V, T ) =

∑∞ N =0

[

∑′ {n^ε}

∏

ε

(ze^−βε)ⁿ^ε]

= ∑

{n0,n1,··· }

(ze^−βε⁰)ⁿ⁰(ze^−βε¹)ⁿ¹(· · · )

=







∏

ε

(1− ze^−βε)⁻¹ BE case, ze^−βε < 1

∏

ε

(1 + ze^−βε) FD case

(28)

The average occupation number in state ε

< nε> = 1 ϑ[−1

β(∂ϑ

∂ε)z,T ,all other ε]

= 1

z⁻¹e^βε+ η

where η = −1 for bosons, η = −1 for fermions, η = 0 for MB gas. It is convenient to dene the thermodynamic potential q to nd the average number and energy

q(z, V, T )≡ lnϑ = P V

kT = η∑

ε

ln(1 + ηze^βε) (2.6)

N = z(¯ ∂q

∂z)_V,T =∑

ε

1

z⁻¹e^βε+ η (2.7)

E =¯ −(∂q

∂β)_z,V =∑

ε

z⁻¹e^βε+ η (2.8)

2.3 Semiclassical Boltzmann Equation(UUB Equa- tion)

Boltzmann equation is successful in describing dilute classical monatomic gases. People showed the derivation of the Boltzmann equation from the Newtonian of motion for the many particle system under suitable assump- tions. Hereinafter, the semiclassical Botlzmann equation (also called UUB equation) is used for describing the ow of fermions and bosons. The generalized Boltzmann equation is

∂_tf + p

m · ∇xf + F · ∇ξf = R(f ), (2.9) where F = ma is the external force. The quantum particles are assumed to obey the streaming mythology like the left hand side of the Boltzmann

(29)

2.3. Semiclassical Boltzmann Equation(UUB Equation) 13

equation. It should be emphasized here that the major and only dierence between semiclassical Boltzmann equation and classical Boltzmann equation is the collision term R(f). In the classical one, collision term is derived based on classical mechanics [46] while it is derived by quantum scattering theory in the semiclassical one [23]. The collision term for semiclassical Boltzmann equation is described below [21].

R(f ) =

∫ dp₁

∫

dΩK(p, p₁, Ω){[1 + ηf(p, t)][1 + ηf(p1)]f (p^∗, t)f (p^∗₁, t) (2.10)

−[1 + ηf(p^∗, t)][1 + ηf (p^∗₁)]f (p, t)f (p1, t)},

where the function K is the collision kernel, Ω is the solid angle. η = +1 de- notes the case of BE statistics and η = −1 denotes the case of FD statistics.

We observe that the Pauli exclusion principle is included in the case of fermion and a population factor for boson is included. It is noted that the Boltzmann equation of the classical statistics is included in (2.10) as a special case when η = 0. The population in nal state is neglected. This shows the classical collision term can be recovered when considering dilute, non-degenerate classical gases.

Although semiclassical Boltzmann equation has been derived from considering quantum particles collisions, the complexity of the collision term is still a obstacle for computation. One straightforward way to solve this problem is to approximate the complex collision term R(f) with BGK approximation

f−f^eq,Q

τ which was rst proposed in [12]. BGK approximation makes the direct solution of the BGK-Boltzmann equation tractable. After introducing the BGK approximation, (2.9) could be rewritten as

∂

∂tf + ξ· ∇xf + a· ∇ξf =−f − f^eq

τ . (2.11)

(30)

In this BGK-Boltzmann equation, fêq could be quantum equilibrium distri- bution function fêq,Q or classical equilibrium distribution function fêq,C. fêq means the probability of particles (electrons, phonons, photons, etc...) occu- pying a specic quantum state in equilibrium, such that, the physical meaning of fêq is the same as < n > which has been derived in last section. The gen- eralized equilibrium distribution function is

f^eq,Q(ε) =< n_ε >= 1

z⁻¹e^βε+ η (2.12)

where z ≡ e^−βµ is the fugacity and ε = m(ξ − u)²/2is the energy of particle.

The parameter η determines which statistics f^eq,Q belong to.

2.4 Semiclassical Hydrodynamic Equations

The semiclassical hydrodynamic equations are obtained by taking moments ψ = [1, mξ, mξ²/2]on the semiclassical Boltzmann equation of (2.9) with the collision term (2.10), then integrating the resulting equations over all ξ.

∂_t

∫

ψf dΞ +∇x

∫

ψξf dΞ− a

∫

ψ∇ξf dΞ =

∫

R(f )ψdΞ (2.13) The integrals of the collision terms in all three cases should preserve the conservation property. That means the conservation of mass, momentum and energy need to be satised all the time, which is called the compatibility condition,

∫

R(f )ψdΞ = 0, ψ =







1 mξ

1

2m(ξ_x²+ ξ_y²+ ξ_z²)





, (2.14)

where dΞ = m³dξ/h³ is the innitesimal volume in momentum space. We note that the collision in BGK model should also preserve the compatibility

(31)

2.4. Semiclassical Hydrodynamic Equations 15

condition,

1 τ

∫

(f − f^eq)ψdΞdt = 0, ψ =







1 ξi 1

2m(ξ_i²+ ξ_j²+ ξ_k²)





. (2.15)

The denitions of the number density, number density ux, and energy density are given, respectively, by

n(x, t) =

∫ f dΞ, j(x, t) = m

∫

ξf dΞ, E(x, t) = m

∫ ξ² 2f dΞ.

Other denitions of higher order moments such as internal energy density ε(x, t), stress tensor Pij(x, t), heat ux vector qi(x, t) are also given,

ε(x, t) = m 2

∫

|ξ − u|²f (ξ, x, t)dΞ

P_ij(x, t) = m

∫

(ξ_i− ui)(ξ_j − uj)f (ξ, x, t)dΞ q_i(x, t) = m

2

∫

|ξ − u|²(ξ_i− ui)f (ξ, x, t)dΞ.

We can have the familiar form of hydrodynamic equations,

∂n

∂t + ∂nu_i

∂x_i = 0 (2.16)

ρdu_i

dt + ∂P_ij

∂x_j = ρai (2.17)

ρdε

dt + P_ij∂u_i

∂x_j = ∂q_i

∂x_j, (2.18)

where, _dt^d = _∂t^∂ + u_i_∂x^∂

i. We have briey described the equilibrium distribution including MB, FD, and BE distributions. The next step is to derive the corresponding hydrodynamic equations under these statistics. Although the

(32)

hydrodynamic equations are the same in both classical and semiclassical ver- sion, the transport coecients are dierent due to the dierences between f^eq,Q and f^eq,C. After some straightforward algebraic manipulations, the equation of state of pressure p, the number density n and total energies E in classical and quantum statistics as:

p_c = nk_BT n_c = z

Λ³ Ec = 3

2nkBT + D 2mnu², and

p_q = nk_BTg_5/2(z) g_3/2(z) nq = g_3/2(z)

Λ³ E_q = 3

2nk_BTg_5/2(z) g_3/2(z)+ D

2mnu².

g_ν is the generalized FD/BE function which is dened as

g_ν(z)≡ 1 Γ(ν)

∫ _∞

0

x^ν⁻¹

z⁻¹e^x+ ηdx, (2.19) wherein Γ(ν) is the Gamma function. Notice that in classical limit, gν(z) will approach z no matter what the number ν is. We can check that the pressure pq, number density nqand total energies Eqwill approach the classical counterpart in classical limit. We can also check the specic heat γ = Cp/C_v, in ideal classical gas,γc = 5/3 for monatomic gas, but γq = ⁵₃^g^5/2_g^(z)g2 ^1/2^(z)

3/2(z)

in ideal quantum gas, which depends on the fugacity z and not constant.

Obviously, γqapproaches γcin classical limit. Moreover, more hydrodynamical

(33)

coecients like viscosity and thermal conductivity are given by

µq = τ nkBTg_5/2(z) g_3/2(z) = τ P κ_q = τ5k_B

2m[7g_7/2(z)

2g_3/2(z)− 5g_5/2(z)

2g_3/2² (z)]nk_BT, (2.20) and the classical coecients are also listed for comparing:

µ_c = τ nk_BT = τ P κ_c = τ5k_B

2mnk_BT = τ5k_B

2mP. (2.21)

All the semiclassical coecients will approach their classical counterparts and Prandtl number equals 1 in classical coecients. However, in semiclassical model, the Prandtl number depends on the fugacity z and is no more a con- stant number. Similar results could be found in linearized semiclassical Boltz- mann equations [23]. Until now, the semiclassical Boltzmann Equation and the semiclassical hydrodynamical equations have been introduced and compared with their classical counterparts. In those descriptions there appears an important parameter z which is not shown in classical Boltzmann equation or MB statistics. The physical meaning of the fugacity z is described below:

Recall number density in semiclassical representation nq = ^g^3/2_Λ3^(z), if we con- sider 0 < z < 1, then g3/2(z) ≃ z, we have z = nq(h²/2πmk_BT )^3/2 = n_qΛ³. From this equation, z can be interpreted as the ratio of Λ³ to average volume occupied by particles. In other words, it is the ratio of occupied length to particles thermal wavelength. When z is small, it means the order of spa- tial dimension is larger than the thermal wavelength and we can neglect the degeneracy eect of particles (highly non-degenerate gas). For large z, the degeneracy eect is important (highly degenerate gas) because the order of thermal wavelength is comparable to the spatial dimension. Moreover, the thermal wavelength is proportional to T^−1/2, that means when the number

(34)

density is in normal condition this eect will be obvious especially in the low temperature. In summary, when z → 0 the quantum distribution will coin- cide with the classical one, and the physical explanation is that the length dimension of particle is larger than the particle de Broglie wavelength. The wave property will not be important. When z is considerably large, two length scales become comparable and one cannot omit the quantum eect anymore.

So, we can think the fugacity z is the index of the degree of degeneracy. The fugacity z has some restrictions in two dierent quantum distributions. In the case of Boson, z should not exceed 1 because of the non-negative density, and in the Fermion case there are no such restrictions on z.

Finally, the actual correction values of the generalized Fermi function is shown in Fig. 2.1. One nds the BE and FD curves overlap MB curve in z → 0 limit, that means the classical statistics could only work in the classical limit (z → 0) of quantum statistics. Moreover, z is restricted in 0 ≤ z ≤ 1 for Bose gas, the function g3/2(z) increases monotonically with z and is bounded, its largest value being g3/2(1) = 2.612.

(35)

Fugacity

0.2 0.4 0.6 0.8

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Bose(2.5) Bose(1.5) Bose(2.5)/Bose(1.5) (a)

Fugacity

0.2 0.4 0.6 0.8

0.2 0.4 0.6 0.8 1

Fermi(2.5) Fermi(1.5) Fermi(2.5)/Fermi(1.5) (b)

Fugacity

0.2 0.4 0.6 0.8

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Bose(2.5)/Bose(1.5) Fermi(2.5)/Fermi(1.5) (c)

Figure 2.1: Bose and Fermi functions in dierent fugacity z. (a) Bose Func- tion g1.5^Bose(z), g2.5^Bose(z), and ^g_g^2.5^Bose^Bose^(z)

1.5 (z), (b) Fermi Function g^{F ermi}1.5 (z),g2.5^{F ermi}(z), and ^g_g^{F ermi}^2.5^{F ermi}^(z)

1.5 (z), (c) Comparing Bose ^g_g^Bose^2.5^Bose^(z)

1.5 (z), and Fermi Function ^g_g^{F ermi}^2.5^{F ermi}^(z)

1.5 (z)

(36)

(37)

Chapter 3

Semiclassical Lattice Boltzmann Method

3.1 Overview

In this chapter, the semiclassical LB method will be derived from semiclassical Boltzmann equation with BGK approximation (2.11). In general, there are two steps for deriving LB method from the continuous Boltzmann equation. First, the time discretization is achieved by integrating the Boltzmann

(38)

equation along characteristic line. Second, the discretization of space is done by low Mach expansion or Hermite polynomials expansion. The derivations of semiclassical LB method is following the Hermite expansion procedures and the results are extended to multiple relaxation time version. The derived equations could be validated by reducing the SLB equation to classical one in classical limit.

3.2 The Derivations of Conventional LB method

Just like other numerical methods, discretization of time and space before simulating continuous models or equations on a discrete digital computer is necessary, the detailed procedures are listed below.

3.2.1 Time Discretization

Consider the semiclassical Boltzmann equation with BGK approximation (2.11), and neglect the forcing term here:

∂_tf + ξ· ∇f = RBGK(f ) =−1

τ[f− f^eq] (3.1) Rewrite the Boltzmann BGK Equation (3.1) in the form:

D_tf + 1 τf = 1

τf^eq (3.2)

where Dt ≡ ∂t + ξ · ∇. Then, integrating (3.2) over a time step δt along characteristics:

f (x + ξδ_t, ξ, t + δ_t) = e^−δ^t^/τf (x, ξ, t) + 1 τe^−δ^t^/τ

∫ _δ_t

0

e^t^′^/τf^eq(x + ξδ_t^′, ξ, t + δ_t^′)dt^′ (3.3)

(39)

3.2. The Derivations of Conventional LB method 23

By using Taylor expansion, and neglecting the high order terms O(δt2), also with τ^∗ = τ /δ_t, the standard LB equation is derived as:

f (x + ξδ_t, ξ, t + δ_t)− f(x, ξ, t) = RBGK(f ) =−1

τ^∗[f (x, ξ, t)− f^eq(x, ξ, t)], (3.4) This equation indicates the basic procedures of LB method: streaming and collision. Since we have got the time discretization LB equation, the following processes will focus on discretizing the equilibrium distribution function f^eq(x, ξ, t) on velocity space .

3.2.2 Space and Velocity Discretization

In LB method, computational domain is discretized on a regular lattice with

xed grid spacing which makes this method simple. Although original LB method is empirically derived from lattice gas automata, however, modern LB method is derived from continuous Boltzmann equation and has a concrete physical interpretation compared with those come from empirical derivations.

The two dierent approaches to discretizing the computational domain are low mach number expansion and hermite polynomials expansion. In the following derivations, particles are in high temperature, low density and under classical statistics are assumed. That means, the distribution follows the classical MB statistics. And the MB equilibrium distribution function is:

f^eq,C = n

(2πθ)^D/2exp[−c²

2θ] (3.5)

wherein n is the number density, D is the space dimension, θ = RT the temperature and the peculiar velocity c = ξ − u. ξ is the molecular velocity and u the mean velocity. C means classical equilibrium distribution function here.

(40)

3.2.2.1 Low Mach Number Expansion

Expanding (3.5) by Taylor expansion, we can get the discrete equilibrium distribution function as:

f^eq,C ≈ n

(2πθ)^D/2exp(−ξ²

2θ)[1 +ξ· u

θ +(ξ· u)² 2θ² − u²

2θ] (3.6) After choosing some suitable weighting and quadrature points to recover the original moment integration, it will become the frequently used LB model. We use D2Q9 model for example for presenting these procedures. First, we need to discretize the momentum space ξ properly. The integration is replaced by the summations as:

∫

ψ(ξ)f^eq,C(x, ξ, t) =∑

a

W_aψ(ξ_a)f^eq,C(x, ξ_a, t), (3.7) where ψ(ξ) is the polynomial of ξ. Then, the above integral is evaluated by a Gaussian-type quadrature:

I =

∫

exp(−ξ²

2θ)ψ(ξ)dξ =∑

a

W_aexp(−ξ²

2θ)ψ(ξ_a), (3.8) A cartesian coordinate system is used to recover the D2Q9 model. Therefore, we set ψ(ξ) = ξx^mξⁿ_y, where ξx and ξy are the x and y components of ξ. Such that, the integrand I will be

I = (√

θ)^(m+n+2)I_mI_n, (3.9)

where Im =∫₊_∞

−∞ e^−ς and ς = ξx/√

2θ or ς = ξy/√

2θ. For D2Q9 model, the third-order Hermite formula for evaluating Im is chosen, i.e. Im =∑3

j=1ω_jς_j^m. The three abscissas of the quadrature are:

ς₁ =−√ 3/2 ς₂ = 0 ς₃ =√

3/2, (3.10)

(41)

3.2. The Derivations of Conventional LB method 25

and the corresponding weight coecients are:

ω₁ =√ π/6 ω2 = 2√

π/3 ω₃ =√

π/6. (3.11)

After this integrating, the equilibrium distribution function of the D2Q9 model is:

f_a^(0),C = ω_an[1 +3ξ_a· u

c² +9(ξ_a· u)²

2c⁴ − 3u²

2c²], (3.12) where c =√

3θ≡ ^δx_δt. δx is the space between two adjacent lattices. Similarly, other models like D2Q7 and D3Q27 can be derived in a similar manner.

3.2.2.2 Hermite Polynomials Expansion

Hermite polynomials expanding of the distribution function for LB method was rstly presented in [17], then well developed in [47] and [16]. This method is inspired by [5], [18] and [48]. The brief descriptions are described below.

Following the approaches in [17][16][47], we adopt the Grad's moment approach and seek solutions to (3.4) by expanding f(x, ξ, t) in terms of Hermite polynomials,

f (x, ξ, t) = ω(ξ)

∑∞ n=0

1

n!a⁽ⁿ⁾(x, t)H⁽ⁿ⁾(ξ) (3.13) And the expansion coecients a⁽ⁿ⁾ are given by

a⁽ⁿ⁾(x, t) =

∫

f (x, ξ, t)H⁽ⁿ⁾(ξ)dξ (3.14) where ω(ξ) = _(2π)¹3/2e^−ξ²^/2is the weighting function, a⁽ⁿ⁾and H⁽ⁿ⁾(ξ)are rank- n tensors and the product on the right-hand side denotes full contraction. Here and throughout the dissertation, the shorthand notations of Grad [5] for fully symmetric tensors have been adopted.

(42)

Some of the rst few tensor Hermite polynomials are given here, H⁽⁰⁾(ξ) = 1

H⁽¹⁾i (ξ) = ξ_i H⁽²⁾ij (ξ) = ξ_iξ_j− δij

H⁽³⁾_ijk(ξ) = ξ_iξ_jξ_k− ξiδ_jk − ξjδ_ik− ξkδ_ij

The rst few expansion coecients can be easily identied with the familiar hydrodynamic variables:

a⁽⁰⁾ =

∫

f dξ = n a⁽¹⁾ =

∫

f ξdξ = nu a⁽²⁾ =

∫

f (ξ²− δ)dξ = P + n(u²− δ) a⁽³⁾ =

∫

f (ξ³− ξδ)dξ = Q + u(a⁽²⁾− 2nu²) (3.15) Where the momentum ux tensor is dened as P = Pij ≡∫

f c_ic_jdc, the heat ux tensor is dened as Q = gijk ≡ ∫

f c_ic_jc_kdc, and We also have ε =

1

2[a⁽²⁾_ii −n(uu−3)]. The orthogonality of Hermite polynomials implies that the leading moments of a distribution function up to the Nth order are preserved by truncations of the higher-order terms in its Hermite expansion. Thus, a distribution function of the Boltzmann-BGK equation can be approximated by its projection onto a Hilbert space spanned by the rst N Hermite polynomials without aecting the rst N moments. Here, up to Nth order, f^N(x, ξ, t) has exactly the same velocity moments as the original f(x, ξ, t) does. This guaranties that a LB method gas dynamic system can be constructed by a

nite set of macroscopic variables. As a partial sum of Hermite series with

nite terms, the truncated distribution function f^N can be completely and uniquely determined by its values at a set of discrete abscissae in the velocity