LBM - gnomeCreative/HYBIRD GitHub Wiki
Outline of the LBM model
In more traditional CFD solvers like the Finite Element Method or the Finite Volume Method, the conservation of mass and momentum is imposed directly on macroscopic quantities, i.e. velocity and pressure, through a solution of the Navier-Stokes equations. In LBM, the fluid is represented as a collection of particles moving in space, conceptually similar to the real composition of a fluid as a molecular system. Like in a real fluid, mass and momentum conservation are applied at the microscopic level whenever a collision occurs. However, the number of molecules present in any practically relevant volume of fluid is so large, that there is no hope to simulate it using computers. LBM solves this problem by drastically reducing the number of degrees of freedom, describing time, space, and velocity in a discrete form. While time and space are discrete also in traditional CFD solvers, the discretization of the velocity space is peculiar to LBM, and is inherited by Lattice Automata, from which LBM stems McNamara & Zanetti (1988).
The LBM cubic lattice (a) and the set of discrete velocities (b) used in the D3Q19.
The theoretical foundations of LBM lies within the kinetic theory of gases. Imagine a system of gas particles, following a chaotic motion typical of particles. This system can be described by employing a probability density function, $f(x,t,c)$, which represent the probability of finding a particle in the position $x$ at time $t$, moving with speed $c$. As mentioned before, the overwhelming amount of degrees of freedom in this description is drastically reduced by employing a discretized time-space. In the LBM formulations presented here, a cubic lattice is used, as represented in the figure above (a). At every time step, the particles located in a node $x_0$ are allowed to move to a neighbor node $x_i$. Since the time step $\Delta t$ is fixed, every movement corresponds to a specific velocity, calculated as
$c_i=\frac{x_i-x_0}{\Delta t}$.
LBM limits the eligible neighbors to a fixed set. The more neighbors are allowed, the higher the resolution power of LBM becomes, together with an increase in computational time, iterations, and complexity of the code. For the simulation of fluid dynamics, a very common compromise consist in choosing the 19 neighbors depicted in the figure above (b), which correspond to the set of velocities:
$\left(0,0,0\right)~\textrm{for}~i=0$,
$\left(\pm1,0,0\right)~\textrm{for}~i=1,2$,
$\left(0,\pm1,0\right)~\textrm{for}~i=3,4$,
$\left(0,0,\pm1\right)~\textrm{for}~i=5,6$,
$\left(\pm1,\pm1,0\right)~\textrm{for}~i=7\textrm{..}10$,
$\left(0,\pm1,\pm1\right)~\textrm{for}~i=11\textrm{..}14$,
$\left(\pm1,0,\pm1\right) \textrm{for}~i=15\textrm{..}18$.
A lattice defined in this way is commonly named D3Q19 (3 dimensions, 19 velocity vectors). At every node $x$ and time $t$, 19 probabilities distribution functions $f_i$ are defined, each corresponding to one of the discrete velocities: $f_i(x,t)=f(x,t,c_i)$. Each of these functions can be imagined as a population of particles streaming with the same velocity and in the same direction. In the following, time and space discretizations are considered unitary, i.e. $\Delta S=1$ and $\Delta t=1$, which greatly simplifies the notation, since the lattice speed $c=\Delta S / \Delta t$ is also unit. The reference density $\rho_\textrm{ref}$ is also unitary, yielding a system of completely dimensionless variables. In an actual simulation, the equations are indeed solved with unit discretizations, and the results are scaled during post-processing. In the following, we will express all units in dimensionless forms, and a few notes on scaling will be given later on.
Once the lattice is defined, the macroscopic variables density $\rho_\textrm{f}$ and velocity $u_\textrm{f}$ are reconstructed at every node by simple summation:
$\rho_\textrm{f}(x,t)=\sum\limits_i f_i(x,t)$,
$u_\textrm{f}(x,t)=\left(\sum\limits_i f_i(x,t) c_i\right)/\rho_\textrm{f}(x,t)$.
The density of the fluid is therefore treated as a variable, and the medium is granted a limited compressibility. The pressure can be reconstructed as
$p_\textrm{f}(x,t) =c_\textrm{s}^2 \cdot \rho_\textrm{f}(x,t)$,
where $c_\textrm{s}$ is the norm of the lattice speed of sound, defined as
$c_\textrm{s} =\frac{c}{\sqrt{3}}=\frac{1}{\sqrt{3}}$.
At every time step, the populations converging to the same node are colliding at that position. The dynamics of the collisions is governed by the Boltzmann equation, another result from the kinetic theory. The evolution equation, in absence of body forces, reads
$\frac{\partial f}{\partial t}+c\nabla f=\Omega_\textrm{coll}$,
where the term $\Omega_\textrm{coll}$ implements the effect of collisions. In the lattice, the Boltmann equation appears in a discretized form
$f_i (x+ c_i,t+1) = f_i (x,t)+\Omega_{\textrm{coll},i} (x,t)$.
The definition of $\Omega_\textrm{coll}$ is a crucial point in the development of a LBM formulation. Its precise expression, as given by Boltzmann, is of no practical use in this case. Fortunately, Bhatnagar, Gross, and Krook (1954) proposed an alternative approximate form of $\Omega_\textrm{coll}$, as an operator simply driving the system towards thermodynamics equilibrium. Its expression reads
$\Omega_{\textrm{coll},i}=\frac{f^\textrm{eq}_i-f_i}{\tau}$,
where $\tau$, is a time constant that governs how hard the system is pushed towards equilibrium. The expression of $f^\textrm{eq}$ follows a Maxwell-Boltzmann equilibrium distribution.
$f^\textrm{eq}=\frac{\rho_\textrm{f}}{(2\pi R T)^{3/2}}e^{\frac{-(c-u_\textrm{f})^2}{2 R T}}$,
with $R$ denoting the universal gas constant and $T$ the temperature. The discretized form of the equilibrium distribution is obtained after a Taylor expansion on $u_\textrm{f}$ up to the second order. Considering that $c=\sqrt{3RT}$, one obtains
$f_i^\textrm{eq} (u_\textrm{f},\rho_\textrm{f}) = \rho_\textrm{f} w_i \left(1+3 c_i\cdot u_\textrm{f} +\frac{9}{2} \left(c_i\cdot u_\textrm{f}\right)^2 -\frac{3}{2}u_\textrm{f}\cdot u_\textrm{f}\right)$.
The set of weights $w_i$ is chosen so that the collision operator conserves mass and momentum. The values for the D3Q19 are
- $w_i = 1/3$ for $i=1$,
- $w_i = 1/18$ for $i=2\textrm{..}7$,
- $w_i = 1/36$ for $i=8\textrm{..}19$.
With the LBM presented so far, the dynamics of a Newtonian fluid can be solved. A realization of LBM can be proven to be analogous to the solution of the Navier-Stokes equation if the Knudsen number is small, i.e. provided that the size of an element is much smaller than the length scale of the system~ (Chapman & Cowling, 1991). Stability and accuracy are however limited to the range of small Mach numbers, i.e. in the simulation the condition $\textrm{Ma} = u_\textrm{f}/c_s \ll 1$ should always be respected. The viscosity of the fluid is directly related to the relaxation time $\tau$, according to
$\nu(x,t)=\frac{\tau(x,t)-1/2}{3}$.
If a volumetric force is acting on the system (e.g. gravity), an additional operator $\Omega_\textrm{force}$ needs to be added to the system. There exist different techniques for the implementation of such an operator. Commonly employed is the procedure suggested by Guo et al. (2002), which, given a uniform force field $F$, consists of the addition of the following operator to the right hand of the lattice-Boltzmann equation:
$\Omega_{\textrm{force},i}=w_i\left(1-\frac{1}{2\tau}\right) \left[3\left(c_i-u_\textrm{f}\right)+ c_i\left(c_i\cdot u_\textrm{f}\right)\right]F$.
The reconstruction of the macroscopic velocity also needs to be modified, becoming
$\rho_\textrm{f}(x,t)=\sum\limits_i f_i(x,t)$,
$u_\textrm{f}(x,t)=\left(\sum\limits_i f_i(x,t) c_i+F/2\right)/\rho_\textrm{f}(x,t)$.
This basic formulation needs to be extended with a few more terms in order to model debris flow. These represent the effects from the non-Newtonian rheology, the free surface and the interaction with particles, and are addressed in this order in the following.