Background theory of fluid dynamics

In classical Newtonian mechanics, the motion of a solid un-deformable body is dictated by the two Newton's laws of motion.

Solving these equations subjected to a set of boundary conditions gives the analytical solution of the motion of the body. We all remember our physics A level where we were asked to solve the motion of a cannon ball subjected to gravity.

This approach could be used to model each molecule of a fluid with the physics of collision implemented on the top of this physics and a simple programme could be very quickly designed. However, a simple calculus show the computational impossibility of this approach.

This glass contains 20 cl of water, that is about 1 mole, ie 10^23 molecules! Even the biggest and fastest computer on earth could not handle a minute fraction of this number of molecules.

Fortunately, 'we are much bigger than what we are made of' (S. Succi) and we can jump to a bigger time and length scale and consider the overall behaviour of the molecules, i.e. the continuum motion of the fluid.

This is what we call fluid dynamics, and the approach is very similar to the classical many-body collisions previously described. We solve the Navier-Stokes equations for an incompressible continuous fluid.

This set of PDE's are the so-called Navier-Stokes equations and define the behaviour of the whole fluid. They can be solved for simple configuration but fall quickly into very complex mathematics and very often no solutions can be found.

An example would be to model the stream lines of a moving formula one racing car.

Too complicated boundaries make this problem impossible to solve analytically.

However, CFD's engineers can solve such a flow (source: www.mclaren.com) by the help of CFD solvers.

There are basically three classic numerical solution strategies: finite difference, finite element and spectral methods. They differ mainly in the way in which the flow variables are approximated and in the discretisation process used. The finite volume method which has its origin in the finite difference method is nowadays the most well-established method. Many commercial solvers have their codes based on these methods.

The following part describes another method, namely lattice Boltzmann, which is the method used in this project.

Background theory of lattice Boltzmann

Lattice Boltzmann has been used as the host of the hydrodynamics simulation. It is a fairly recent method based on the well-known Boltzmann probability distribution function. The following aims to give the reader the essential concept of lattice Boltzmann.

Consider 1000 molecules of a fluid in a volume V =(Dx)^3

Question: how many of them have their velocity between v and v+Dv?
Answer: the number denoted f (x,v,t,Dx,Dv).

Making Dx and Dv very small transforms the number f into a probability: f(x,v,t). In the proper jargon, f is the density of probability associated to the velocity v at the position x and at a time t. It represents the probability of finding a particle at x having a velocity v at a time t. This approach leaves the classical physics and its law of motion and dives into statistical physics dynamics and its probabilistic motion.

In 1872, Boltzmann derived the evolution of such probabilities (f's) (refs: [Chapman], [Cercignani], [Harris]):

Consider a gas in which an external force mF acts (where F is the vectorial sum of the external forces such as gravity, electric, etc.) and assume initially that no collisions take place between the gas molecules (i.e. perfect gas). In time dt the velocity c of any molecule will change to c+Fdt and its position x will change to x+cdt. Thus the number of molecules f(x,c,t)dxdc is equal to the number of molecules f(r+cdt,c+Fdt,t+dt)drdc.

If, however, collisions do occur between the molecules there will be a net difference between these numbers. This can be written W(f)drdcdt where W(f) is the collision operator. This gives the following equation describing the evolution of the distribution function : the Boltzmann equation.

The fluid density, velocity and internal energy can be found from the distribution function f. Any solution of Boltzmann's equation (2.3) requires that an expression is found for the collision operator. Without knowing the form of W(f) there are however several properties which can be deduced. The collision is set to conserve mass, momentum and energy. In the simplest approach, the collision term is substituted with a phenomenological term (Bhatnagar-Gross-Krook approach, source : Phys, Rev. 94, 511, 1954): 

where feq indicates the (local) equilibrium distribution function, and is t a microscopic relaxation time.

The Boltzmann equilibrium distribution function

The previous part introduced the notion of local equilibrium with the simple BGK collision operator. This idea is central in lB and differs from global equilibrium where the revolutionary potential of the system has been exhausted (entropy of the system is maximal). Mathematically, it is defined by an exact balance between gain and loses for the velocities so that the collision term is annihilated (it is thermodynamically equivalent to the first case of previous part where the molecules where supposed not to undergo any collision): W(feq)=0.

It can be shown that  where c is the peculiar speed c=v-u, namely the relative speed of the molecules with respect to the fluid, the quantity VT is the thermal speed associated with the fluid temperature T.

It should be noted at this point that this equilibrium function represents a gaussian distribution function centred on u in the velocity space and with T as band width.

To help visualise this equilibrium distribution function, consider the following: each of these fluids has the same uniform velocity u, the graph show the probability of finding a particle with its velocity around u, the animation's show the associated movement of the particles.

Note that this distribution accords with the so-called thermal agitation.

The model we use is supposed to be isothermal, which simplifies considerably the algorithm and parametrisation.

Discretisation: from BE to LBE (2nd order accurate for NS)

This continuum distribution is not very handy to use as such in a computer model as the distribution function assumes an infinite number of different velocities. It is therefore a necessity to discretise the possible velocities of the particles. By doing so, time or space has to be also discretised ... or both. The latter approach is taken as it simplifies greatly the parameter space and stay coherent with the continuum approach. Below is the so called D2Q9 (2 dimension, 9 velocities) scheme.

Each node of the lattice has Q-1 links where densities (particles) propagate each time-step. The D2Q9 model is said to have three levels of energy (three set of velocities): the rest velocity (energy=0, velocity=0), the short links (energy=1/2, velocity=1) and the long links (energy=1, velocity=2^(1/2)).

Each particle is only allowed to reside on the nodes of a lattice (yellow dots) and can only move onto a neighbouring node (via the links in blue). A lattice based automaton also discretises time: the particles only interact (collide) at very precise and regular instant of time and do not see each other for the rest of the time. The 'life' of a particle on a lattice based automaton (such as this one) could be summarised as

The evolution equation of the system does not depend on the symmetry of the lattice and can be written as: 

where I denotes the link (velocity) in the discretised space (lattice). This approach has been directly inspired by the BGK's previously described and has appropriately been named Lattice BGK, or LBGK for short. It is difficult to imagine a simpler scheme that recovers the Navier-Stokes equation (to second order).

LBGK equilibria

Qian et al. (1992) have shown that a whole range of geometry's was possible for the LBGK model. They dubbed them DnQm for m speed model in n dimensions. Most popular examples are D1Q3, D1Q5, D2Q9. An important asset of this scheme is that there is no problem in going into three dimensions. Straight forward in other words. In fact, there are in fact a number of robust three-dimensional LBGK's such as the D3Q15 and D3Q19.

The generic family of LBGK equilibria can be expressed as where r is the fluid's nodal density, wi is the weight associated to the i-th direction (velocity) for the given geometry and Qiab is the projector along the i-th discrete direction.

Restrictions on mass and momentum conservation as well as isotropy impose a set of constrains on the weights wi. Under-constrain due to the extra freedom provided by the multi-energy (multi-speed) lattice (more discrete speed than constraints) allows several solutions (geometry's). The weights can either be interpreted as degeneracy numbers from higher dimensions or simply as different masses of the particles moving along the different directions (from [Succi]). Below is a table depicting some of the geometry's and the associated weights:

Succi notes that there is a real analogy between the LBGK weights and the numerical stencils (sometimes called 'computational molecules') of finite difference and finite element calculations.

Note that to populate the fluid at rest in the lattice, we set the densities to equilibrium term with zero velocity and the targeted density.

To find more: a very good reading is the book from S. Succi: The Lattice Boltzmann Equation, Oxford University Press, Clarendon Press,ISBN: 0198503989 (available in any good retailer shop) which gives an overview of the lattice Boltzmann method and place it within the CFD's world.