INTRODUCTION
Computational Fluid Dynamics (CFD) has emerged as a powerful tool for the analysis
of system involving fluid flow, heat transfer and associated phenomena such
as chemical reactions, evaporation, condensation, etc (Santos
and Costa, 2009; Suresh and Anthony, 1996; Su
and Dong, 1999; Bolokhonov et al., 2006).
From 1960s onwards, the aerospace industry has integrated CFD techniques into
the design, research and development and manufacturing of aircraft and jet engines
(Bianco et al., 2009; Coiro
and Nicolosi, 2001). More recently, CFD has been applied to the design of
internal combustion engines, combustion chambers of gas turbines and furnaces
(Kumaran and Babu, 2009; Kumar and
Kale, 2002; Abbassi and Khoshmanesh, 2008). Furthermore,
motor vehicle manufacturers now routinely predict drag forces, underbonnet
airflows and the incar environment with CFD (Efisio et
al., 2008; Tsubokuraa et al., 2009; Toumi
et al., 2009).
The fundamental law of any fluid flow problems is the NavierStokes equations, which define any singlephase fluid flow. These equations can be simplified by removing terms describing viscosity to yield the Euler equations. Further simplification, by removing terms describing vorticity yields the full potential equations. Finally, these equations can be linearized to yield the linearized potential equations.
Historically, Finite Difference Method (FDM) (Richardson,
1911) was the first computational method used by researchers to solve fluid
flow and heat transfer problem by solving NavierStokes equation. However, due
to the frustration on FDM, which cannot be effectively used on complex geometry,
Finite Element Method (FEM) has been introduced in 1950s (Strang
and Fix, 1973). In 1980s, Finite Volume Method (FVM) was developed at Imperial
College, mainly to solve fluid dynamic problems (Patankar,
1980). Since then the finite volume method is extensively used to solve
transport phenomena problems. Indeed, the FDM, FEM and FVM belong to the same
family of weighted residual method, however, they limit their simulation in
the range of continuum fluid.
There are few numerical methods that simulate the evolution of fluid flow at
particle level. Among them are Direct Simulation Monte Carlo (DSMC) (Bird,
1963) and Molecular Dynamic (MD) (Alder and Wainwright,
1957) methods. In these methods, the trajectories of every particle together
with their position in the system are predicted using the second Newton’s
law. From the knowledge of the forces on each atom, it is possible to determine
the acceleration of each atom in the system. These methods are deterministic;
once the positions and velocities of each atom are known, the state of the system
can be predicted at any time in the future or the past. But remember, a cup
of water contains 10^{23} number of molecules. Even when a gas is being
considered where there are fewer molecules and a larger timestep can be used,
because of the longer mean free path of the molecules, the number of molecules
that can be considered is still limited. However, the question is do we really
need to know the behavior of each molecule or atom? The answer is no. It is
not important to know the behavior of each particle, it is important to know
the function that can represent the behavior of many particles (mesoscale).
As a result, in 1988, the lattice Boltzmann method (LBM) (McNamara
and Zanetti, 1988), a mesoscale numerical method based on statistical distribution
function has been introduced to replace MD and DSMC methods.
Historically, LBM was derived from Lattice Gas Automata (LGA) (Frish
et al., 1986). Consequently, LBM inherits some features from it precursor,
the LGA method. The first LBM model was a floatingpoint version of its LGA
counterpart. Each particle in LGA model (represented by single bit Boolean integer)
was replaced by a single particle distribution function represented by a floatingpoint
number. The lattice structure and the evolution rule remain the same. One important
improvement to enhance the computational efficiency has been made for the LBM
was that the linearization of collision operator (Bhatnagar
et al., 1954). The uniform lattice structure was remaining unchanged.
The starting point in the lattice Boltzmann scheme is by tracking the evolution
of the singleparticle distribution function. The concept of particle distribution
has already well developed in the field of statistical mechanics while discussing
the kinetic theory of gases and liquids (Harris, 1971).
The definition implies the probable number of molecules in a certain volume
at a certain time made from a huge number of particles in a system that travel
freely, without collision, for distance (mean free path) long compared to their
sizes. Once the distribution functions are obtained, the hydrodynamics equations
can be derived.
Although LBM approach treats gases and liquids as systems consisting of individual particles, the primary goal of this approach is to build a bridge between the mesoscopic and macroscopic dynamics, rather than to deal with macroscopic dynamics directly. In other words, the goal is to derive macroscopic equations from mesoscopic dynamics by means of statistic, rather than to solve macroscopic equations.
The LBM has a number of advantages over other conventional computational fluid
dynamics approaches. The algorithm is simple and can be implemented with a kernel
of just a few hundred lines (Azwadi and Tanahashi, 2006).
The algorithm can also be easily modified to allow for the application of other,
more complex simulation components. For example, the LBM can be extended to
describe the evolution of binary mixtures, or extended to allow for more complex
boundary conditions (Abe, 1997; Benzi
and Succi, 1990; Bernsdorf et al., 2000;
Azwadi and Tanahashi, 2007). Thus the LBM is an ideal
tool in fluid simulation.
The objective of present study is to introduce and discuss the formulation of LBM in simulating fluid flow problem. The derivation of macroscopic continuity and momentum equations from mesoscale Boltzmann equation is discussed in details. After showing how the formulation of LBM fits in to the framework of macroscale flow, numerical results of liddriven cavity flow, natural convection in an enclosure and dynamics of droplet on solid surface are presented to highlight the applicability of the approach in various fields of fluid dynamics.
LATTICE BOLTZMANN FORMULATION
Isothermal lattice Boltzmann method: Ludwig Boltzmann (18441906) introduced
a transport equation based on statistical mechanics describing the evolution
of gas particle in a system as:
where, f, c, a and Ω stand for density distribution function, mesoscopic speed, acceleration due to external force, collision function and external force respectively. If there is no external force, Eq. 1 is no more than a hyperbolic wave equation with source term given as:
Any solution of the Boltzmann equation, Eq. 2, requires an expression for the collision operator Ω. If the collision is to conserve mass, momentum and energy, it is required that:
However, the expression for Ω is too complex to be solved. Even if we
only consider twobody collision, the collision integral term needs t
o consider the scattering angle of the binary collision, the speed and direction
before and after the collision, etc. Any replacement of collision must satisfy
the conservation law as expressed in Eq. 3. The idea behind
this replacement is that large amount of detail of twobody interaction is not
likely to influence significantly the values of many experimental measured quantities
(WolfGladrow, 2000).
There are a few version of collision operator published in the literature.
However, the most well accepted version due to its simplicity and efficiency
is the BhatnagarGrossKrook collision model (Bhatnagar
et al., 1954) with a single relaxation time. The equation that represents
this model is given by:equilibrium distribution function and τ is the time
to reach equilibrium condition during collision process and is often called
the relaxation time.
Equation 4 also describes that 1/τ of nonequilibrium
distribution relaxes to equilibrium state within time τ on every collision
process. Substituting Eq. 4 into Eq. 2 gives:
which is known as the Boltzmann BGK equation.
Equation 4 describes two main processes at mesoscale level. The left hand side refers to the propagation of distribution function to the next node in the direction of its probable velocity and the right hand side represents the collision of the particle distribution functions. In lattice Boltzmann formulation, magnitude of c is set up so that in each time step Δt, every distribution function propagates in a distance of lattice nodes spacing Δx. This will ensure that distribution function arrives exactly at the lattice nodes after Δt and collides simultaneously.
In order to apply Eq. 4 into the digital computer, the mesoscopic
velocity space has to be discretised. This can be done by discretising the physical
space into uniform lattice nodes. Every node in the network is then connected
with its neighbours through a number of lattice velocities to be determined
through the model chosen. The general form of the lattice velocity model is
expressed as DnQm where D represents spatial dimension and Q is the number of
connection (lattice velocity) at every node (He and Luo,
1997; He and Doolen, 2002). There are many lattice
velocity models published in the literature, however, the most well used due
to its simplicity is D2Q9 and shown in Fig. 1.

Fig. 1: 
D2Q9 lattice Model 
After discretisation in velocity space, Eq. 4 can be rewritten
in the following form:
where, i refers to the number of discrete velocity. σ refers to the direction of mesocopic velocity where σ = 0 when i = 0, σ =1 when i = 1,3 , 5, 7 and σ =2 when i = 2, 4, 6, 8.
The Eulerian expression of the left hand side of Eq. 5 can be transformed into the Lagrangian form. To do this, we take Euler time step in conjunction with an upwind spatial discretization and then setting the grid spacing divided by the time step equal to the velocity. This leads to the wellknown lattice Boltzmann BGK equation:
where, ε is a small lattice time unit in physical unit.
Derivation of macroscopic equations: The equilibrium distribution function
f^{eq }is choosen so that we can reconstruct the hydrodynamics of fluid
flow. The general form of f^{eq }can be written as:
Here A_{σ}., B_{σ}., C_{σ}.. and D_{σ} are the coefficients to be determined based of ChapmannEnskog procedure. For the rest particle, Eq. 7 becomes:
and
for other particles. This gives:
The symmetric properties of the tensor Σ_{i}c_{σiα}c_{σiβ}... are needed in the derivation and given as follow:
• 
The odd orders of tensor are equal to zero 
• 
The second order tensor satisfies where
δ_{αβ} is the Kronecker delta and c_{1} =
1 and 
• 
The fourth order tensor has an expression as: 
for σ = 1 and
for σ = 2 where:
By considering conservation laws of:
and
results in:
These give constraints for coefficients A_{σ}., B_{σ}., C_{σ}..and D_{σ} as:
and
To satisfy Eq. 14 we chose ,
and
We next decompose the timescale into slow and fast timescale. This is to represent two different phenomena occur at different timescale such as advection and diffusion:
where, ε plays the role of Knudsen number (Hou et
al., 1995). We also expand f about f^{eq}:
Here to
imply that the
nonequilibrium distributions do not contribute to the
local values of density and momentum.
Taylor expanding of Eq. 6 and retaining terms up to O(ε^{2}) results in:
Substituting Eq. 17 and 18 into Eq.
19, the equation to order of ε and ε^{2} are:
and
respectively. Equation 21 can be further simplified by using
Eq. 20 gives:
The first order continuity equation can be obtained by taking summation of Eq. 20 respect to σ and i
Taking the same summation of Eq. 22 gives:
Combining Eq. 23 and 43 gives the correct
form of the continuity equation:
We next multiply Eq. 20 with c_{iσ} and taking the summation as above gives:
where, is
the momentum flux
tensor. Substituting the expression of the equilibrium distribution, Π^{eq} can be written as:
which are the pressure term and two nonlinear terms. This gives:
and
where, c^{2}_{3}= 1/3 is speed of sound. In order to obtain a velocity independent pressure and Galilean invariance, we choose:
and
This gives the final expression for Π^{eq} as:
Substituting Eq. 32 into Eq. 26 results in Euler equation as:
and the pressure is given by p = c^{2}_{3}ρ.
We next multiply Eq. 22 with c_{iσ} and taking the summation respect to σ and i:
Substituting the expression of the nonequilibrium distribution and using Eq.
25 and 32 leads to:
To maintain isotropy, we set:
Recalling Eq. 16, gives the expression for B_{1} and B_{2}:
Using Eq. 33 in the form:
Equation 35 can be written as:
Combining Eq. 33, 34 and 39
gives the momentum equation in two dimensions:
where:
For an incompressible fluid, the momentum equation becomes:
The remaining coefficients are determined as follow:
Finally, the equilibrium distribution functions can be written as follow:
and
From the above derivation, we can see that the mesoscale Boltzmann equation can lead to the macroscopic NavierStokes equation by the Chapman Enskog expansion.
Thermal lattice Boltzmann method:In general, the current thermal lattice
Boltzmann models fall into three categories: the multispeed approach (McNamara
and Alder, 1993), the passive scalar approach (Shan,
1997) and the thermal energy distribution model proposed by He
et al. (1998). The multispeed approach uses the same distribution
function in defining the macroscopic velocity, pressure and temperature. In
addition to mass and momentum, in order to preserve the kinetic energy in the
collision on each lattice point, this model requires more variations of speed
than those of the isothermal model and equilibrium distribution function usually
include higher order velocity terms. However, this model is reported to suffer
severe numerical instability and is not computationally efficient (McNamara
and Alder, 1995).
In the passive scalar model, the flow fields (velocity and density) and the
temperature are represented by two different distribution functions. The macroscopic
temperature is assumed to satisfy the same evolution equation as a passive scale,
which is advected by the flow velocity, but does not affect the flow field.
It has been shown that the passive scalar model has better numerical stability
than the multispeed model (Eggels and Somers, 1995).
He et al. (1998) in their model introduce the
internal energy density distribution function, which can be derived from the
Boltzmann equation. This model is shown to be a suitable model for simulating
real thermal problems. However, the complicated gradient operator term appears
in the evolution equation and thus the simplicity property of the lattice Boltzmann
scheme has been lost. To overcome this problem, current author has developed
the simplest lattice structure for internal energy density distribution function
by neglecting the viscous and compressive heating effects. To see this, the
new variable, the internal energy density distribution function is introduced:are
the equilibrium distribution function for internal energy, the term due to the
external force and the heat dissipation term respectively.
Substituting Eq. 47 into Eq. 4 results
in:
where
and
Note that the external force is reintroduced and two different relaxation times
are used to characterize the momentum and energy transport. Equation
47 represents the internal energy carried by the particles and therefore
Eq. 48 can be called as the evolution equation of internal
energy density distribution function. The macroscopic temperature field can
be defined in term of distribution function g as:
As mentioned earlier, for simulating incompressible flow, the viscous heat
dissipation can be neglected. So the evolution equation of internal energy density
distribution function is reduced as follow:
By omitting the dissipation term, the complicated gradient operator in the
evolution equation of internal energy distribution function can be dropped.
In the previous section, the equilibrium function for density distribution
f^{eq} was derived following the same procedure as in Lattice Gas Automata
(LGA) (Frish et al., 1986). In other word, our
understanding of the basis of LBM has been restricted by our knowledge of the
statistical mechanics of LGA. Here, the derivation of equilibrium function for
internal energy density distribution function will be carried out, starting
from the MaxwellBoltzmann equilibrium distribution function written as:
where, R is the ideal gas constant, D is the dimension of the space and ρ,
u and T are the macroscopic density, velocity and temperature respectively.
Recall Eq. 49 and expanding up to u^{2} for the internal
energy distribution function gives:
Regroup to avoid higher order quadrature gives:
It has been proven by Shi et al. (2004) that
the zeroth through second order moments in the last square bracket and the zeroth
and first order moments in the second square bracket in the right hand side
of Eq. 56 vanish. The exclusion of the second order moments
in the second square bracket in Eq. 56 only related to the
constant parameter in the thermal conductivity which can be absorbed by manipulating
the parameter τ_{g} in the computation. Therefore by dropping the
terms in the last two square brackets on the right hand side of Eq.
56 gives:
For low Mach number flows, the internal energy density equilibrium distribution
function can be further simplified by neglecting the terms O(u^{2})
gives:
To recover the macroscopic energy equation, the zerothto secondorder moments
of g^{eq} must be exact. In general:
where, I_{g} should be exact for m equal to zero two respectively.
Now a new variables φ_{g} is introduced defined by:
Rewriting Eq. 59 after substituting Eq. 60
gives:
For simplicity the microscopic velocity is normalized as
Equation 62 can be calculated by using the GaussHermite
quadrature. Hence, the GaussHermite quadrature must consistently give accurate
result for quadratures of zerothtothirdorder of velocity moment of g^{eq}.
This implies that the secondorder GaussHermite quadrature can be applied in
evaluating I_{g}. Therefore Eq. 62 can be evaluated
as follow:
where, N is the number of abscissas, W_{k} is the corresponding weight
coefficient and ζ_{k} is the Gaussian abscissa. For two dimensional,
D = 2, secondorder GaussHermite quadrature N = 2, the value for W_{k}
and ζ_{k} are:
After some modification in order to satisfy macroscopic energy equation via
ChapmannEnskog expansion procedure, the discretised internal energy density
distribution function is obtained as:
This type of lattice structure for internal energy density distribution is
shown in Fig. 2.
From above derivations, we can see that the evolution equation for internal energy distribution function can be directly derived from the Maxwell Boltzmann distribution function.
Multiphase lattice Boltzmann method: Recently, a number of researchers
have used LBM to study multiphase fluid flow in some specific engineering problems
such as granular flow and droplet dynamics (Briant et
al., 2002; Swift et al., 1995).

Fig. 2: 
Lattice structure for internal energy density distribution
function 
Microscopically, the phase segregation and surface tension in multiphase flow
are because of the interparticle forces/interactions. Due to its kinetic nature,
the LBM is capable of incorporating these interparticle interactions, which
are difficult to implement in traditional methods.
In general there are three types of lattice Boltzmann models have been advanced
to simulate multiphase flow systems. The first type is the socalled colored
model for immiscible twophase flow proposed by Gunstensen
et al. (1991) and based on the original lattice gas model by Rothmann
and Keller (1988) and Gunstensen et al. (1991)
used colored particles to distinguish between phases. The color model was further
developed by later studies (Grunau et al., 1993),
but it has serious limitations. One of the most significant problems is that
the model is not rigorously based upon thermodynamics, so it is difficult to
incorporate microscopic physics into the model (Boghosian
and Coveney, 2000).
The second type of LB approach used to model multicomponent fluids was derived
by Shan and Chen (SC model) (1993) and later extended
by others (Shan and Chen, 1994). In the SC model, a
nonlocal interaction force between particles at neighboring lattice sites is
introduced. The net momentum, modified by interparticle forces, is not conserved
by the collision operator at each local lattice node, yet the system’s
global momentum conservation is exactly satisfied when boundary effects are
excluded (Martys and Chen, 1996). Hou
et al. (1997) compared the above two types of models for simulating
a static bubble in a two fluid system and concluded that the SC model is a major
improvement over the colored model.
The main drawback of the SC model, however, is that it is not wellestablished
thermodynamically. One cannot introduce temperature since the existence of any
energylike quantity is not known (Hazi et al., 2002).
The third type of LB model for multiphase flow is based on the FreeEnergy
(FE) approach, developed by Swift et al. (1995),
who imposed an additional constraint on the equilibrium distribution functions.
The FE model conserves mass and momentum locally and globally and it is formulated
to account for equilibrium thermodynamics of nonideal fluids, allowing for the
introduction of well defined temperature and thermodynamics (Swift
et al., 1996). The major drawback of the FE approach is the unphysical
nonGalilean invariance for the viscous terms in the macroscopic NavierStokes
equation. Efforts have been made to restore the Galilean invariance to secondorder
accuracy by incorporating the density gradient terms into the pressure tensor
(Kalarakis et al., 2002).
Present study focuses on the multiphase LBM which is based on the work of ShanChen which involves evolution equation of single particle distribution function f and can be written as:
To simulate multiphase fluids, longrange interactions between particles are
needed. To do this, an interaction force between the nearest neighbors of fluid
particle is incorporated for single component multiphase fluid. Therefore, the
external force F in Eq. 61 is given as:
where, ω is the weight as in equilibrium distribution function. G is the
interaction strength and Ψ is the interaction potential. There are few
types of interaction potential exist in the literature. Some of them are Ψ(ρ)
= ρ_{0}[1exp(ρ/ρ_{0})] (Raiskinmaki
et al., 2002), Ψ(ρ) = ρ [20], Ψρ_{20}ρ^{2}/[2(ρ_{0}+ρ)^{2}]
(Pan et al., 2004), Ψ(ρ) = Ψ_{0}exp(ρ_{0}/ρ)
(Shan and Chen, 1993), etc.
Following Shan and Chen (1994), the nonideal equation
of state from the abovementioned external force and interaction potential can
be written as:
One of the advantages of ShanChen multiphase model over other LBM multiphase
models is its capability to simulate the dynamics of droplet on solid surface.
The study of droplet dynamics is very close to the phenomena in daily lives
such as droplet sliding on car’s front glass or windows, spreading of water
on table, droplet falling and many more to mention. This phenomenon also plays
an important role in real engineering applications. Coating of substances, boiling
water reactor, injection of ink are some of them.
Other than cohesive force described earlier, the interaction between fluid particles and surface is needed to include the adhesive force between these two phases. The interaction force between the fluid particles at site x and the solid wall at site x' is formulated as:
The coefficient to describe the strength of force between solid and fluid is
different to that of between fluid and fluid. Therefore we need to use a different
two coefficients to characterize these two types of forces. At the fluidsolid
interface, the solid is regarded as a phase with constant density s,
which is 1 for a solid and 0 for a pore. Using these definitions, the fluid
momentum is changed at each time step according to:
where, u' is the new fluid velocity.
SIMULATION RESULTS
Earlier the formulation of mesoscale lattice Boltzmann scheme and the derivation of macroscopic NavierStokes equation were demonstrated. Here, the numerical prediction of liddriven cavity flow, natural convection in an enclosure and droplet dynamics on solid surface are performed to demonstrate the applicability of the scheme.
Liddriven cavity flow: The liddriven cavity flow has been used as
a benchmark problem for many numerical methods due to its simple geometry and
complicated flow behaviors. It is usually very difficult to capture the flow
phenomena near the singular points at the corners of the cavity.
In this subsection, the mesoscale LBM scheme is applied to this liddriven cavity flow of height L. The top plate moves from left to right along the x direction with a constant velocity U and the other three walls are fixed. In the simulation, the Reynolds number is defined as:

Fig. 3: 
Streamline plots for Re = 1000 and Re = 5000 
Two different values of Reynolds number, Re = 1000 and Re = 5000 are chosen
and the solutions at steady state are compared with the benchmark results (Ghia
et al., 1982).
Figure 3 show plots of streamline for the Reynolds numbers
considered. They are apparent that the flow structures are in good agreement
with the results published in the literature by previous studies (Ghia
et al., 1982; Sidik et al., 2008,
2009). For Re = 1000, the primary vortex appears at
the cavity center and circular shaped. In addition to the primary, a pair of
counter rotating eddies develop at the lower corners of the cavity. At Re =
5000, a third secondary vortex is evolved in the upper left corner of the cavity
and the size of the secondary vortices become larger.
The velocity components along the vertical and horizontal lines through the cavity center together with the benchmark solutions are shown in Fig. 4a and b. Good agreement between the mesoscale LBM and the benchmark solutions are observed. It is noted that, the mesoscale LBM is able to capture the critical points in the tested problem.
Natural convection in an enclosure: Flow in an enclosure driven by buoyancy
force is a fundamental problem in fluid mechanics. This type of flow can be
found in certain engineering applications within electronic cooling technologies,
in everyday situation such as roof ventilation or in academic research where
it may be used as a benchmark problem for testing newly developed numerical
methods. A classic example is the case where the flow is induced by differentially
heated walls of the cavity boundaries.

Fig. 4: 
Velocity components along the vertical and horizontal lines
through the cavity center for (a) Re = 1000 and (b) Re = 5000 (lines: LBM,
symbol: Ghia et al., 1982) 
Two vertical walls with constant hot and cold temperature is the most well
defined geometry and was studied extensively in the literature.
In this study, we reproduce the phenomenon of natural convection in differentially heated walls by using LBM scheme due to vast numerical and experimental data can be obtained for the sake of comparison. There are two dimensionless parameters which govern the characteristic of thermal and fluid flow in the enclosure; the Prandtl and Rayleigh numbers defined as follow:

Fig. 5: 
Streamlines plots for (a) Ra = 10^{3} and (b) Ra
= 10^{5} 

Fig. 6: 
Isotherms plots for (a) Ra = 10^{3} and (b) Ra =
10^{5} 
where, χ, L and ΔT are the thermal diffusivity, width of the cavity
and temperature different between left and right walls respectively. In present
study, the Prandlt number of 0.71 was used to represent the circulation of air
in the system. The Rayleigh number is set at 10^{3 }and 10^{5}.
The plots of streamlines and isotherms are shown in Fig. 5a,
b and 6a, b.
At Ra = 10^{3}, streamlines are those of a single vortex, with its
center in the center of the system. The corresponding isotherms are parallel
to the heated walls, indicating that most of the heat transfer is by heat conduction.
As the Rayleigh number increases (Ra = 10^{5}), the central streamline
is elongated and two secondary vortices appear inside it. The isotherms initially
parallel to the differentially heated walls at low Re become horizontal at the
center of the cavity at high Re indicating that the dominant of heat transfer
mechanism is convection. All of these findings are in good agreement with previous
studies (Azwadi and Tanahashi, 2006, 2007,
2008; Azwadi and Syahrullail, 2009;
Azwadi et al., 2010; Azwadi
and Irwan, 2010; Davis, 1983; Fusegi
et al., 1991).
Droplet dynamics on solid surface: The droplet initially positioned
just touching the solid surface. In this case, the magnitude of G' has to be
determined before the calculation begins.

Fig. 7: 
Computed results at 35°, 55° and 150° contact
angles. (a) Contact angle 35°, (b) Contact angle 55° and (c) Contact
angle 150° 
The value of G' at special contact angles (0°, 90° and 180°) can
be easily determined by balancing the cohesive and adhesive in different ways.
This gives the value of G' of 327.79, 187.16 and 46354 for contact angle
of 0°, 90° and 180°, respectively.
The simulated results at three other values of G' are shown in Fig.
7ac. The contact angle of the computed figures are then
determined using Auto CAD software.

Fig. 8: 
Comparison of computed droplet’s diameter to height ratio
with analytical solution 
In present investigation, we also calculate the magnitude of the ratio of droplet’s
wetting diameter to the height of droplet for every contact angle. Then the
computed values are compared with the analytical solution and shown in Fig.
8. As can be seen from the Fig. 8, these ratios are in
excellent agreement when compared with the analytical solutions.
CONCLUSSIONS
This study discussed the theory of mesoscale numerical approach namely the lattice Boltzmann method in prediction of various types of fluid flow problems ranging from isothermal to multiphase cases. The derivations of macroscopic governing fluid flow equations were demonstrated is detail based on the kinetic theory of gases. After showing how the formulation of the mesocale particle dynamics fits in to the framework of lattice Boltzmann simulations, numerical results of flow in a liddriven cavity, natural convection in an enclosure and droplet dynamics on solid surface were conducted to demonstrate the applicability of the present method. All of the predicted results were found to be in close agreement with the previous similar study reported in the literature.
The lattice Boltzmann method is still undergoing development. Many models including simulation of flow in porous media, solidfluid flow, were recently proposed. Although innovative and promising, these existing LBM methods, require additional benchmarking and verification. This would open many new areas of application.