The Lattice Boltzmann Method Biology Essay

Published: Last Edited:

This essay has been submitted by a student. This is not an example of the work written by our professional essay writers.

Blood is one of the most important parts of the human body as it performs several fundamental physiological functions, such as transferring oxygen and other nutrients to tissues, removing any waste from organs, transporting red blood cells, white blood cells, platelets and antibodies for self-vaccination, as well as normalizing the human body's pH and temperature.

However, the flow of the blood can be disturbed and considered abnormal from a broad range of disorders and human diseases, including heart diseases, anemia (several kinds), malaria, diabetes, and atherosclerosis. Therefore, it comes in great need to be able to fully understand the rheological properties of blood, in order to be able to use that knowledge in several sectors such as biomedical or bioengineering. This could lead to innovations such as development of blood substitutes as well as cures for common but lethal diseases.

Visualizing blood at a microscopic depth, blood is a highly concentrated suspension of numerous components such as red blood cells (or erythrocytes), white blood cells (or leukocytes) and platelets. All of the above cells and platelets are suspended in a solution with similar properties as water, which contains a variety of molecules.

Figure 1: The components of the blood: White blood cells, red blood cells and platelets. Penn Medicine,(2007).[1]

Among the components of the blood, it can be said that the RBC's (red blood cells) play the most important role in the mechanics of the blood, mainly due to their significant mechanical properties and their large number of density which is . Blood flow varies in parts of the human body. For instance, in large scale flows such as those occurring in the heart or other large in diameter blood vessels, blood behaves like a fluid whose as the shear rate becomes higher, the viscosity decreases which means that the blood is thinning. On the other hand, in small-scale flows such as the outer branches of the circulatory system (microcirculation), where the diameter is smaller than 100 μm, the components of the blood are more apparent since there are less, and the mechanics of the red blood cell (RBC) plays a significant role in the self-regulation of the blood flow, which leads to a healthy organism.

Several in vitro and in vivo studies have shown that the extremely flexible RBCs tend to migrate to the central region of a capillary or a tube in a normal blood flow. Moreover, a cell-free layer of only plasma appears at the bottom of the capillary walls, which acts as a lubricant for the ease of the flow. Therefore, looking closely at the capillary or the tube, a core in high viscosity of cells appears during a regular flow. This behavior agrees with the Fahraeu-Lindqvist effect, R. Fahraeus and T. Lindqvist ,(1931) [2], where there is a decrease in viscosity as the diameter of the capillary decreases. However, this phenomenon occurs only in vessels which are smaller than 300μm. Additionally, the mechanics of the blood flow in microcapillaries agrees with the Fahraeus effect ,(1931),[2], which states that the hematocrit of the blood decreases as the diameter of the capillary decreases. In other words, the concentration of the RBCs in human blood depends on the diameter of the vessel. It should be noted that the Fahraeus effect influences the Fahraeus-Lindqvist effect but it has not be proven that the former is the only cause of the later.

According to Popel and Johnson ,(2005),[3], the blood flow formation plays an important role in the circulation of blood in the microvascular system of the human body by determining the resistance of the flow and the transfer of the biological properties.

During the 1960s, several vital studies of the human blood flow and microcirculation advanced majorly due the rapid technological innovations in video microscopy, computer software, the development of methods for modeling the motion of the cell as well as the medical gain of knowledge for the exchange of water and macromolecules. The advances in microscopy have allowed researchers to find significant information for the structure and properties of the cells as well as of the vessels of the microcirculatory system. They were able to gain more knowledge about the mechanism that causes biological exchange and find out how all the components of blood, vessels and tissues are linked together to form a perfect working system. However, all of the above developments and applications, not only supply information but also raise even more questions as the procedure becomes more complicated. Therefore, a very deep analysis of the physical processes is required to be able to understand the fundamentals of the human physics. In blood flow, there are many important physical properties such as velocity, pressure, shear stress, molecular concentration and stress of the viscous-elastic membrane of the RBC. All of these properties are hard to measure in such a small scale (capillary) in the laboratory due to the fact that the in vivo process cannot be translated exactly to the in vitro. For example, the contractions and the structure of the vessel is very hard to copy in an in vitro experiment which makes the results doubtful of realism.

In recent years, there has been an immense progression in computational methods that are able to measure properties of the blood by numerical simulation that could not be measured directly the years before. The flow- induced deformation of capsules with several different shapes such as spherical or spheroidal in an unstressed shape and fluids (plasma) with different or equal viscosities with the RBC have become possible to simulate.

Two very important and many times used methods to describe the deformation of the RBC are the boundary integral method and the immersed boundary method, Pozrikidis (1992),[4] Peskin (1977) [5]. Using the boundary integral method has the advantage of using discretization only in the surface of the membrane of the RBC where it is required. A concern using that method is the restriction of the creeping flow, however, it is not of much importance since the blood flow occurs at very low Reynolds number, which makes it convenient for this otherwise ineffective method. Many successful applications of this method have been reported such as of the RBC motion in shear and channel flow from Pozrikidis and that of the white blood cell migration in a suspension in channel flow by Freund (2007)[6].

The most common way to represent the RBC is the immersed boundary method where the interaction between the fluid(plasma) and the membrane(RBC) is implemented, by membrane forces at the fluid and afterward by updating the membrane configuration grid according to the flow of the plasma. This method comes with several advantages such as that the problem can be solved by standard numerical schemes over a mesh. The greatest advantage of the IBM is that it can be easily combined with a wide variety of fluid formulations such as finite element, finite volume and lattice-Boltzmann method. Other numerical approaches that have been used in the past comprise of multi-particle dynamics, lubrication approximations, Secomb et al. (2001)[7], dissipative particle dynamics, Pivkin & Karniadakis (2008) [8] and liquid drop model.

Especially for the modeling and simulation of red blood cells, vessels, tissues and other biological parts that needed to be studied, the immersed boundary method has been widely used, Sui et al. (2008).[9]. A further important contribution to the study of the RBC was of Eggleton and Popel (1998),[10]where they conducted a research on the red blood cell dynamics in shear flow. Moreover, Bagchi(2007),[11], examined and modeled the flow of suspensions in extremely small vessels of 20 to 300 μm of diameter. He performed a successful numerical simulation of an RBC with a spherical shape in a parabolic vessel flow (Poisseuille flow) and observed the migration of the RBC towards the centerline of the vessel. In the case of the liquid drop model, it was observed that the higher the deformation of the RBC, the faster the migration velocity of it towards the centre. Moreover, there have been many researches for the remarkable deformability of the membrane of the RBC, including that of Zhang et al. (2007)[12] which gave an insight for the cell aggregation. Except from the RBCs, platelets and white blood cells have been given enough attention from authors such as AlMomani(2008)[13].

Despite the evolution in computing, there are still many restrictions that present researchers with computational challenges. Many realistic attempts have been made in order to simulate the RBC properly, as with a biconcave resting shape and high internal viscosity, but unfortunately a few have been successful. Takagi et al.(2009)[14], performed successful simulations of the RBC in shear flow using the immersed boundary method, taking into account the elastic nature of the membrane and its flexural stiffness. Their results showed that the liposomes which are enclosed by membranes deform in a different way than RBCs.

Another efficient ethod to simulate RBC deformation is the lattice-Boltzmann method. Dupn et al.(2007),[15] , managed to simulate the deformation of the cell in three dimensions. The results of this approach are quite similar to the results using traditional methods. However, some parameters which are adjustable have to be chosen. Moreover, it is difficult using this method to be able to extract information from the results such as membrane tension and perturbation flow. Perhaps most important is the fact that Dupn at al. failed to take into account the effect of the membrane viscosity.


The goal of this project is to be able to offer any further documentation and physical approaches in the motion and deformation of an RBC in infinite shear flow using the lattice-Boltzmann method and the immersed boundary method. Attention will be given to the internal viscosity of the membrane and the effects of it in the motion of the cell. Moreover, great emphasis will be given in the spontaneous migration of the cell to the centerline of the flow. However, the simulation of the RBC in three dimensions in a biconcave shape is not possible due to computational restrictions which give persistent challenges related to the adaptive mesh of the membrane. Nonetheless, a two dimensional model can give enough information concerning the dynamics of the cell in a two dimensional flow.


Throughout recent years, the method that has been very promising in the field of fluid flow simulations and fluid modeling, is the lattice Boltzmann method. The scheme has been proven very successful in applications involving fluid flows with interfacial dynamics and complicated boundaries. Traditional methods involved the discretization of macroscopic equations whereas the lattice Boltzmann method is not based on a macroscopic field but on mesoscopic kinetic equations and microscopic models. Not having to turn the equations from a macroscopic to a microscopic field has made this method very attractive to scientists whose research is based on a micro environment. The principal idea of the LBM is to form straightforward kinetic models that consist of the vital physics of microscopic or mesoscopic equations. The ground for using these simple kinetic equations for macroscopic flows, according to Kadanoff (1986),[16] is that the dynamics of a fluid in the macroscopic field comprises of many microscopic particles which work in a collective way and act as a unite system. Moreover it is assumed that the macroscopic dynamics is not affected by the details of the microscopic field. It is an easy method to use, since by developing a simpler version of the kinetic equation, solving complex equations such as the Boltzmann equation is avoided. Moreover, in molecular dynamics simulations, one has to follow each particle separately, which is not necessary in the LBM. However it provides many advantages of molecular dynamics such as clear physical pictures, simple implementation of even complex boundary conditions and parallelization with other algorithms and methods such as the immersed boundary method.

Although it can be said that the LBM is based on a particle frame, its main focus is the mean macroscopic behavior. In current years, very fast computing machines are available which makes sense to employ codes that can use the intrinsic features or parallelism. The LBM fulfills these requirements in an uncomplicated manner.

There are three very important features of the LBM that make it stand apart from the rest of the numerical methods. First and most important is that the streaming operator, or convection operator, in velocity space is linear. It is essential to try to turn the system its nonlinear nature to a linear environment. This characteristic is borrowed from kinetic theory and is a huge advantage compared to other traditional methods that use a macroscopic approach with nonlinearity. The simple streaming operator is combined with a collision operator and gives the recovery of the nonlinear advection through multi-scale expansions.

Second feature is that the incompressible Navier-Stokes equations can be retrieved in the limit of the LBM which is nearly incompressible. This means that the pressure of the lattice is computed with the use of an equation of state. On the other hand, in traditional direct numerical methods that are used to simulate the incompressible Navier-Stokes equations, the pressure has to satisfy the Poisson formulae with the velocity strains acting as the basis. This equation describing the pressure of the flow can be very challenging to solve, as it comes with several numerical difficulties and requires iterations or relaxations. Using the LBM, this time consuming and complex iteration is avoided.

The third feature is that the LBM uses a certain amount of velocities in space. In another traditional kinetic theories such as the Maxwell-Boltzmann equilibrium distribution, the space is utterly functional. This means that when averaging the macroscopic field, the whole velocity space is involved which means that complex arithmetical calculations are necessary. On the other hand, in the LBM this averaging is greatly simplified because of the minimized velocity directions. The transformation from macroscopic in mesoscopic field consists of easier numerical calculations than other methods.


The lattice Boltzmann method originated from lattice gas automata (LG), which is a discrete particle kinetics that utilizes distinct lattice and distinct time. Another way to view LBM is like a kind of finite difference scheme for the kinetic equation of the velocity distribution formulae. Broadwell (1964),[17] was one of the first to simulate fluid flows in order to examine shock structures, using the simplified kinetic equation with a single particle velocity . As it was a beginner's steps, Broadwell's model can be viewed as a one-dimensional lattice Boltzmann based method. Another method that has been used to study shock structures is the multispeed discrete particle velocity models which was used innovatively by Inammuro and Sturtevant (1990) [18]. Common variable in these models is that although the particle velocity in the distribution function was discretized, time and space were continuous, unlike today's lattice Boltzmann model. It was not until 1976 when Hardy et al. (1976) used the fully discrete particle velocity model where all three components, velocity, time and space were discrete, in order to examine the transport properties of fluids. Later, Frisch et al. (1986) [19], influenced the way the lattice gas automaton method was used for two-dimensional hydrodynamics. They acknowledged the heavy importance of the symmetry of the lattice for the recovery of the Navier-Stokes equation. Remarkably, for the first time they managed to obtain the correct form of the Navier-Stokes equation with the use of the lattice gas automata on a hexagonal lattice, as a starting point. After that breakthrough, more ideas came by such as the cellular automaton model by Wolfram (1986) [20] and the three dimensional model using four dimensional face centered hyper cubic lattice by d'Hummieres et al (1986)[21].

The lattice gas automaton is structured as a simple imaginary molecular dynamic in which the velocity of the particles, space and time are all discrete which gives the method a great advantage in numerical computations. For this reason, the lattice gas method can be also referred to as the lattice gas cellular automata. The general structure of the lattice gas automaton is a regular lattice with particles residing on the nodes of it. The straightforward operation behind the lattice is that a set of Boolean variables that indicate that the particle occupation is defined. The development function of the lattice gas automata is as follows:

Where M is the number of directions of the velocities of the particles at each node and are the local velocities of the particles. The process in which the method works is quite simple to understand. In the initial stage, the configurations of all the particles at each time step progress in two sub-steps, streaming and collision. In the streaming step, each particle moves to the nearest node in the direction of its given velocity. The collision step occurs when particles that arrive at the next node (from the streaming step) interact and change their velocity directions according to scattering rules. An important feature of this method is that only one particle is allowed to exist at a discrete time, with a discrete velocity at a discrete node. This is called the exclusion principle and it is imposed to this method for memory efficiency and simplicity. According to Frisch et al. (1987),[22] it can be said that this leads to a Fermi-Dirac local equilibrium distribution.

The main different feature of the lattice Boltzmann method to the lattice gas automata is that the particle occupation variables (Boolean variables) , which are present in equation 1, are replaced by single particle distribution functions , where < > denoted the ensemble average. Another important feature is that the individual particle motion and particle with particle correlations in the kinetic equations are being neglected. This simple procedure helps to reduce the statistical noise. In this method, the primal variables are the averaged particle distributions which fall in the mesocsopic category of variables. However, as it was mentioned before, the locality of the equation and its advantages still remain since the kinetic form of the LG automata is the same. This makes it a unique feature since locality is important in order for the method to be able to be parallelized with other methods.

A significant problem with the LBM is that the collision operator was non-linear. However, Higuera and Jimenez (1989)[23], assumed that the distribution is close to the local equilibrium state and proposed an approach for the collision operator which is linearly stable. This was an important step for the further simplification of the LBM, since this better version of the collision operator makes use of the relaxation time towards the local equilibrium using a single time relaxation.

Bhatnagar-Gross-Krook (BGK),(1954)[24], are famous for the relaxation term which was given their name (BGK collision operator) and has been used by several different authors when using lattices for researches. In this model, the local equilibrium distribution is chosen to form the Navier- Stokes macroscopic equations which makes the numerical computations far more resourceful and allows flexibility of the transferring coefficients.


LBE: An Extension of LG Automata

There are many ways to retrieve the lattice Boltzmann equation (LBE) and also several ways to obtain the macroscopic Navier-Stokes equations.

Navier-Stokes equations


Discrete velocity model

Boltzmann kinetic equation

As it can be seen from the graph, two ways to obtain LBE is from either a discrete velocity model or the Boltzmann kinetic equation. Once having the LBE, there are several straightforward ways to obtain the Navier-Stokes equations. Additionally, since the LBM is an advanced lattice gas method, an LBE will be introduced which is derived from a discrete kinetic equation for the particle distribution formula which is similar to equation 1 shown before.

Where is the velocity distribution function of the particle in the ith direction, Δt and Δx are time and space increments and is the collision operator which is the rate of change of the particle velocity distribution function which results from the collision process. An interesting fact is that when the space over the time increment equals the local velocity of a particle, , both of equations 1 and 2 have the exact same discretizations. Moreover, it is essential to point out that depends only on the local distribution function. In the LB method, space in conveniently discretized in a certain way that it is consistent with the kinetic equation, a feature that makes calculations simpler. This means that the coordinates of the nearest points around a position of a particle, are .

Other properties such as density ρ and momentum density ρu are parts of the distribution function and can be characterized as particle velocity moments. The following equations represent the density and the momentum density respectively.

Moreover, is required to satisfy the conservation of mass and conservation of momentum at each lattice, meaning:

Performing a Taylor series expansion in the space and time increment, the kinetic equation in second order in continuum form is obtained:

Furthermore in order to derive the macroscopic hydrodynamic equation, the Chapman-Enskog expansion is being employed which is mainly a multi-scaling expansion, Frisch et al. (1987).

From the above equation it can be assumed that the convection time scale is much faster than the diffusion time scale . Moreover, in the same expansion way, thevelocity distribution function can be formed to include the local equilibrium distribution function . The expansion gives the following result:

Where is the non equilibrium distribution function and can be described by:

In this case the local equilibrium function depends only on the local macroscopic variables of the density and the momentum density. However there are two constrains that should be held:

Moreover, there are two constrains for the non equilibrium function :

For k=1 and k=2.

Now by substituting into the collision operator , the Taylor expansion gives:


Higuera and Jimenez (1989), noticed that in equation 5, when ε tends to zero, then This provides us with a linearized collision operator:

Where is the collision matrix that is rensponsible for the scattering rate of the directions i and j.

The collision matrix depends solely on the angle between the directions i and j. However this matrix has to satisfy the mass and momentum conservation as shown by Benzi et al. (1992)[25]:

Furthermore, if it can be assumed that there is a single rate Ï„ in which the local particle distribution relaxes to equilibrium, then the collision matrix becomes:

Bhatnagar at al.(1954)[27], managed to conclude to the LBGK collision term by substitution of equations 6, 10 and 12:

And furthermore the complete LBGK equation:

The BGK collision term has been previously used from others in Boltzmann simulations in order to examine shock dynamics and sometimes with finite-volume methods.

From equation 5 it can be obtained that:


The first order equation above, using equation 15, can be rewritten in the form of:

Afterwards the mass and momentum equations can be obtained with the use of equations 15 and 17:


Where Π is the momentum flux tensor and can be described as:

Where is the component of the velocity vector in α coordinate direction and is the component of the velocity vector in the β coordinate direction.

The detailed form of the momentum flux tensor can be described by the specification of the lattice structure and the corresponding equilibrium distribution (number of velocities and directions). In order to simplify the problem, a two dimensional square lattice with 9 velocities is considered (Figure 2). The velocities can be described by:

Figure 2: A two dimensional square lattice with nine directions. For i=0 the velocity is considered to be zero. (2008)[26]

According to Frisch et al (1986), the choice of that specific lattice is not random. It is necessary to sustain the lattice symmetry in order to derive the Navier-Stokes equations correctly. Thus choosing a simpler five velocity square lattice is not going to provide the symmetry needed in this situation. However, the simplest lattice in two dimensions where the symmetry is guaranteed is the nine velocity square lattice for the LBE to recover the Navier-Stokes equations.

The problem although rises from the fact that the Navier-Stokes equations have a second-order nonlinearity. Chen et al.(1992)[27], noticed that the general form of the equilibrium distribution function can be written as :

Where α, b, c and d are constants. This expansion by Chen at al. can only be used for small velocities such as in microfluidics, or even small Mach numbers.

In order to obtain the above lattice constants (a,b,c,d) ,according to Qian et al.(1992)[28] the constrains in equation 7 have to be utilized. Therefore:

with and and

Afterwards, by substituting the above equation in equation 20 the result is:

Where p is the pressure and it is equal to ρ/3 and ν is the kinematic viscosity which can be described by ν = (2τ-1)/6. As a result of the above, the momentum equation becomes:

Qian and Orszag(1993)[29], showed that the above momentum equation is exactly the same with the Navier-Stokes equations if the density variation δρ is really small.

LBE: Approximation to the Continuum Boltzmann Equation

Recently it has been proven by He and Luo(1997)[30], that despite the fact that the LBE evolved from the LG automaton method, it can also be obtained from the continuum Boltzmann Equation for discrete velocities, by using a very small Mach number expansion which agrees with the fact mentioned above, that the density cariation δρ is really small in the Navier-Stokes equations of momentum. In these deriviations, the source is the Boltzmann Bhatnagar Gross Krook equation:

Where is the single-particle distribution function in continuum phase space and can be described by = (x,,t) and is the Maxwell-Boltzmann equilibrium distribution function and can be described as:

Where D is the spatial dimension. Moreover, ξ which is the particle velocity and u which is the fluid velocity have been normalized by , where T is the temperature. The velocity moments of the function g, are the macroscopic fluid variables:

Where ε is the internal energy and it can be described as: ε=DT/2. Koelman(1991)[31] assumed that the fluid velocity in equation 26, is a small parameter and showed that the equilibrium distribution can be described by:

The kinetic evolution in the BGK equation 25, will only need the solution of . In addition, in discrete velocity models, only a few particle velocities and their distribution functions at these velocities are used. Moreover, the two definitions of equation 27, ρ and ρu, can be approximated with the help of the discrete velocities in a Gaussian-type quadrature, which are the following:

An effective distribution function must be constructed that will satisfy the conservation relations in equation 3. This distribution function is:

It is assumed that and are constant and that the function satisfies the following equation which g satisfies as well:


He and Luo(1997) [32], in order to obtain the weight ,use the Hermite formula which is of third order, in order to estimate the integrals in equation 27. Moreover, Abe(1997) [33] assumes that is of straightforward form which is based on . In both papers, the square lattice with nine velocities was used (D2Q9) and it was shown in both that and and which are the same as in the equilibrium distribution function in equation 22 which confirms the results of Quian et al.(1992).

It was observed from Cao et al, (1997) [34], that in equation 31, if the time derivative is replaced by a first order time difference and the fist order discretization is used as well as the collision term is used, then the finite difference equation for is obtained:

Where and . It can be observed that in the above equation, if and then it becomes the standard LBE which was introduced in equation 2.

From this discretization process it can be seen that equation 32 has solely first order convergence in space and time ( and ) to equation 31. On the other hand, Sterling and Chen(1996) [35] explained that because equation 32 has a Langrangian nature, a second order accuracy in space and time can be achieved with the addition of a special form of the discretization error in the viscous term.


It is known that the truncuation error of the LBM is second order in space. Aside from that, the accuracy of the solution depends greatly on the boundary conditions that are applied on the problem and it is found to be only first order in many cases. It is very important to be able to understand the effect of the boundary conditions applied since their use can affect the results of modeling and simulations. In this section, a boundary condition is discussed, the bounce back

boundary condition.

The origin of the wall boundary conditions in the LBM was the LG method. An efficient example of boundary condition is that from Wolfram(1986)[36], where he used a particle distribution function bounce back scheme at walls in order to gain no-slip conditions. It is important to notice that one of the great advantages of the LBM is the easy implementation of boundary conditions, especially no-slip velocity conditions, which makes it ideal for simulating fluid flows in complicated geometries.


One of the simplest way to create solid walls in simulations involving LBM is the bounce-back boundary condition. Essentially, in this condition particles that fall at a point in the wall bounce back with the same velocity to the point they came from. This regulation create stick boundary conditions at one and a half distance along the vector that joints the solid and the fluid nodes. This ensures that the velocity of the fluid that contacts the solid is equal to the velocity of the latter. However, a common case that appears in problems is that the zero-velocity plane must be located exactly inside the imposed boundary layer in order to achieve better accuracy. In this case an interpolation scheme must be applied.

Figure 3: Bounce back boundary condition for no-slip(rigid wall) and slip(movable wall).(2008)[37]

Moreover, it should be noticed that an alternative condition that is used by many authors is to place the boundary nodes midway between the solid wall nodes and the fluid nodes, thus three of the velocity nodes of the lattice will be touching the solid wall inside the boundary layer.

It is understood that the bounce-back boundary condition is similar to a no-slip boundary condition where the velocity of the fluid relative to the wall is zero. However, to achieve better accuracy in simulations, many authors have imposed much more sophisticated and complicated boundary conditions such as second-order boundaries, where the no-slip boundary is exactly at the wall nodes. The disadvantage unfortunately of such boundaries is that they can only be applied to regular geometries such as flat walls, which makes it unappealing to many authors that look to research problems with more complex geometries. Therefore, in order to achieve practical simulations with small errors, the bounce-back boundary is ideal since it can be computationally efficient in several problems.

Previous studies including the bounce-back boundary condition, were of simple geometry of the walls such as a flat wall. Recently, evaluations of the bounce-back rule have been reported from Gallivan et. al (1997)[38]and compared to results of evaluations made with the finite difference method for flow around octagonal and circular objects. In this particular research, the no-slip boundary was taken to be at the wall node and the error in the results was first-order convergent in space, as usual. The results showed that is can be very difficult to predict the location of the no-slip wall for specifically arbitrary geometries. Moreover, an additional research that compared the lattice Boltzmann method, the finite element method and experimental results of a fluid flow in chemical reactor showed that the lattice Boltzmann results were satisfactory. To conclude, the results of the researches showed that the bounce-back boundary is very useful in complicated geometries. However, in certain applications, satisfactory accuracy can only be gained on large lattices. In cases like so, better accuracy can be achieved for the smaller lattices by refining the grid (or mesh) in the vicinity of walls.

It should be noted that there are certain surface forces that is a result from the bounce back rule of the particles which are calculated from the momentum transfer from each node and summed up to obtain the force and torque on each obstacle object (like the RBC). A great advantage over the finite element and finite difference methods where the local surface normals are required in order to integrate the stresses over the obstacle's surface, in the bounce back rule, these difficulties are removed since the surface forces are directly summed up.

Numerical Efficiency, Stability and Nonuniform Grid

The hyperbolic equation which is an approximation of the Navier-Stokes equation,in this chapter shown as equation 31, in the nearly incompressible limit. According to Reider and Sterling(1995),[39] the numerical accuracy of the equation depends on the Mach number. The velocity is a constant vector which is not the case with the dependent velocity in compressible hydrodynamic equations. This dependence prevents the LBM from solving nonlinear Riemann problems. The LBM s considered to be a relaxation method for the macroscopic equations. However, it is very similar to the "pseudocompressibility method" according to Sterling and Chen (1996). This view, by Ancona (1994) [40], was used to construct the LBM to include Lagrangian methods for the solution of partial differential equations. Jin and Xin (1995) [41] proposed a way to solve a hyperbolic conservation system, the kinetic relaxation method. This method, in order to model the nonlinear flux terms in the equations in the macroscopic field, it uses the relaxation approach and Riemann solvers are not needed. Furthermore, Jin and Katsoulakis (1997) [42], using the above relaxation method managed to calculate curvature-dependent front propagation. In principle the kinetic relaxation was developed for shock capture. On the other hand, LBM is used mainly on viscous nearly incompressible flows. However, the relaxation method and the LBM are very alike. A simulation method was proposed by Nadiga and Pullin (1994)[43] for discrete velocity gases. This method was proved to be more flexible and simple to use as well as more accurate. In 1995 this technique was more developed in order to solve the compressible Euler equation. Furthermore, matters of convergence, stability and numerical efficiency were studied by Elton et al. (1995) for lattice Boltzmann models. Succi et al. (1991) [44] showed that the number of calculations would be for the LBM, while for the pseudospectral method it would be , where N is the number of lattice points along one direction and r is 150 for the 2D square lattice. Comparing further the two methods, Chen at al. (1992)[45] reported that for three-dimensionalsimulations, the LBM can be almost 2.5 times faster than the pseudospectral method, which is a remarkable difference considering the simplicity and efficiency of the LBM. Noble et al (1996) [46] studied the scalability of the LBM. Their results showed that there is a much greater liner scalability when using the LBM, when while increasing the size of the computational domain, at the same time the number of processors is increased, which makes perfect sense. Noble et al. found that for a two dimensional problem which was modeled with lattices, they accomplished 13.9 gigaflops with the use of 512 processors of CM-5 machines. However, the computational capabilities have changed dramatically since then and computational efficiency can be achieved easier.

Equation 26, which is the equilibrium distribution function in the Boltzmann equation, can be translated as the maximum entropy state. This is called the Boltzann's H-theorem and means that any initial state will continue to a state of higher entropy. This theorem guarantees stability as it ensures an increase in entropy. Moreover, Frisch et al.(1987) developed an H-theorem for the LG method. On the other hand, in the LBM method, only a few discrete velocities of the particles are used and the equilibrium distribution is reduced to . As a result, an H-theorem cannot ensure that when put together with the LBM it will produce the correct form of the macroscopic equations. Thus, the LBM often encounters numerical instability.

It is known that the original LBM equation (equation 2) is a finite-difference solution of the Boltzmann equation (equation 31) and the discretization process of equation 32 with a=1, which satisfies the Courant-Friedrichs-Lewy condition. The CFLC is a required condition when solving partial differential equations, especially hyperbolic ones, using the method of finite differences. It is very important in computational simulations as it means that the time step must be less than a certain time set by the scientist, in many explicit time simulations or else the results will be incorrect.

Furthermore, Sterling and Chen (1996) proposed a von Neumann linear stability analysis for the LBM. Basically, they developed to become , where are the constant global populations which depend on the mean density and constant velocity. Afterwards, they expanded equation 2 using Taylor series and developed the following linear equation:

Where is the independent of time and location Jacobian matrix that corresponds to the coefficient of the linear term in the original LBM equation (equation 2). Taking out the Fourier transformation of equation 33, in order to analyze the stability, and resolving the eigenvalue problem for . Substituting k with zero, the results showed the linear stability condition which is correct considering the positivity of the viscosity. Chen et al. (1992) carried out a stability analysis for the hexagonal two-dimensional LBM model and the square two-dimensional LBM model. Moreover, Qian et al. (1992) also carried out a stability analysis for the D3Q15 (three-dimensional with 15 velocities) model. Both of them concluded that in order for the LBM to be stable, it is required that the mean velocity of the flow should be belowa maximum velocity which is dependent of several parameters such as the time of the relaxation, the sound speed and the wave number. Ancona(1994) and McNamara et al.(1995)[47] showed that an improvement of the numerical stability of the LBM can be made by starting from the kinetic equation (equation 31) and discretizing using several finite difference methods.

For years, a major obstacle that stood in the way of achieving better numerical efficiency in the LBM, is that the discretization of the Boltzmann equation (equation 31) was constrained in a form of uniform, regular lattices. This was a restriction especially to flows with complex geometries. In order to overcome this problem, Koelman (1991) proposed the inclusion of a spatial non-uniformity in the structure of the lattice. This idea was developed by Nannelli and Succi (1992)[48] and Amati et al.(1997) [49] where they created the finite volume lattice Boltzmann method (FVLBM) which is based on the Boltzmann equation (equation 31). They created a non-unifrom lattice whose cell enclosed many original lattices. Thus they created the equation where C is the cell and is the volume of the cell. This equation requires the estimation of the flux across the boundaries of the cell C which could be approximated by a constant or a linear interpolation. This method was extremely useful in studies involving flows past "fake" bodies. Amati et al. (1997) simulated a three-dimensional turbulent channel flow using the FVLBM model. The results showed that this method is a lot faster than any of the traditional non-uniform CFD methods. Most simulation results agree with the results from the traditional methods, however, this comparison is not satisfactory, partially due to the fact that the interpolation schemes are of low order, which makes the method less successful. A non-uniform finite difference LBM for the Boltzmann equation (equation 31) was developed by Cao et al. (1997) where he modeled annular flows in polar coordinates. Moreover, He et al. (1996)[50] developed a time-dependent interpolation scheme which increases the density of the grid in regions where the shear tension is higher than other regions. This interpolation scheme is used in many simulations including finite element and finite difference as well as LBM, as they managed to simulate successfully flows in a two-dimensional channel with sudden expansion and a flow around a circular cylinder.

Flow Around a Circular Cylinder

Since this project focuses on blood flow inside a capillary and the effects of the flow on an erythrocyte, it can be assumed that the flow is almost same with that of a flow around a circular cylinder with the only difference being that the cylinder is not moving. However, the flow should have similar effects on the moving erythrocyte since the velocity is not high and the Reynolds number is low. Many studies have been made about the flow around a circular cylinder using the LBM. A very important study of a flow around an octagonal cylinder was made by Higuera and Succi (1989), where they studied the dynamics of the flow with Reynolds number up to 80. Their results showed that at Reynolds number of 52.8, the flow became periodic after initially being transient for a large amount of time. At a Reynolds number of 77.8, a shedding of the periodic flow appeared. Moreover, they compared their results from the study including Strouhal number, lift coefficients, drag coefficients and flow-separation angle with previous studies both in vitro experimental and simulation results. The comparison showed that their results agreed with the previous researches except from a small number of errors. Wagner (1994), [51] studied the pressure distribution on the surface of the cylinder due to steady flow. His results showed that there is a less of 3.5% of a difference between his Strouhal numbers and those from past simulations. It is very important to note that the LBM scheme has a compressible nature which makes Wgner's result extremely significant since it shows that in a nearly incompressible limit, the LBM gives very accurate results for the pressure distribution in incompressible fluids , which was demonstrated as well from Martinez et al. (1994)[52]. These two simulations of flow around a circular cylinder were performed using uniform lattices which unfortunately do not approximate well the circular geometry. Another drawback was that the no-slip boundary was placed on the rough mesh point and not on the cylinder. The problem of a flow around a circular cylinder was revised by He and Doolen(1997) [53], where they based their study on an interpolation LBM. In the above simulation, the lattices are square and a spatial interpolation was used which was previously used by He et al. on a previous research .Unlike previous researches the domain of the simulation was circular and they a polar coordinate was used. Thus they managed to apply the no-slip wall boundary condition at exactly on the circular domain of the cylinder. It was clear that the accuracy of the simulation was far better than any previous results that used LBM. The grid lines showed that vortex shedding was appearing which was very similar to experimental results. Along with Strouhal numbers, the LBM calculated lift and drag coefficients at different Reynolds numbers and were compared with different simulations. The computational error was small and similar to errors among other ways of simulating flows around a circular cylinder, in a circular domain.


Concluding Remarks

The LBM is a different numerical method for the computation of incompressible or nearly incompressible fluid dynamics. It has the potential to be an efficient solver for flows with low Reynolds number, flows with complex geometries and multiphase flows.

The LBM can be assumed that it belongs to a category of pseudo compressible solvers of the Navier-Stokes equations for flows in the nearly incompressible limit. Therefore, the method has the advantages and the limitations such as restrictive time step of the pseudo compressible methods.

LBE is capable of mimicking complex fluid dynamics which are similar to those of the continuum based models in fluid dynamics. It should be noted that the LBM models have a limited application to multiphase hydrodynamics due to the restriction in implementing fluid properties such as surface tension, density ratio, kinematic viscosity and others.

It should be noted that the LBM is different than methods which are directly based on the Navier-Stokes equations as far as computational aspects concern. These different features which make the method unique are:

The non-linearity of the Navier-Stokes equation is limited since the non-linear terms are included in the quadratic velocity terms of the (equilibrium distribution function). Therefore the second order differential equations of the Navier-Stokes are replaced by the first order partial differential equations of the lattice-Boltzmann equation.

The LBM although it can inefficient for time-dependent problems, due to the complexity and computational time especially when the resolution is high, it still is a very competitive CFD technique as its simplicity is very attractive to researchers.

The convergence rate can be extremely improved both for steady and unsteady flows.

Other CFD techniques wich are based on the Navier-Stokes equations, need to solve the Poisson equation in order to compute the pressure. However, the LBM calculates the pressure through the equation of state. This advantage offers the researcher an easier parallelization and faster performance.