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, under-bonnet
airflows and the in-car 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 Navier-Stokes equations, which define any single-phase 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 Navier-Stokes 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 Newtons
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 1023 number of molecules. Even when a gas is being
considered where there are fewer molecules and a larger time-step 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 floating-point 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 floating-point
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 single-particle 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 lid-driven 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 (1844-1906) 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 two-body 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 two-body interaction is not
likely to influence significantly the values of many experimental measured quantities
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 Bhatnagar-Gross-Krook 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 non-equilibrium
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.
|| 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 well-known lattice Boltzmann BGK equation:
where, ε is a small lattice time unit in physical unit.
Derivation of macroscopic equations: The equilibrium distribution function
feq is choosen so that we can reconstruct the hydrodynamics of fluid
flow. The general form of feq can be written as:
Here Aσ., Bσ., Cσ.. and Dσ are the coefficients to be determined based of Chapmann-Enskog procedure. For the rest particle, Eq. 7 becomes:
for other particles. This gives:
The symmetric properties of the tensor Σicσ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 c1 =
||The fourth order tensor has an expression as:
for σ = 1 and
for σ = 2 where:
By considering conservation laws of:
These give constraints for coefficients Aσ., Bσ., Cσ..and Dσ as:
To satisfy Eq. 14 we chose ,
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 feq:
imply that the
non-equilibrium 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:
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 ciσ and taking the summation as above gives:
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:
where, c23= 1/3 is speed of sound. In order to obtain a velocity independent pressure and Galilean invariance, we choose:
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 = c23ρ.
We next multiply Eq. 22 with ciσ and taking the summation respect to σ and i:
Substituting the expression of the non-equilibrium distribution and using Eq.
25 and 32 leads to:
To maintain isotropy, we set:
Recalling Eq. 16, gives the expression for B1 and B2:
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:
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:
From the above derivation, we can see that the mesoscale Boltzmann equation can lead to the macroscopic Navier-Stokes equation by the Chapman Enskog expansion.
Thermal lattice Boltzmann method:In general, the current thermal lattice
Boltzmann models fall into three categories: the multi-speed approach (Mc-Namara
and Alder, 1993), the passive scalar approach (Shan,
1997) and the thermal energy distribution model proposed by He
et al. (1998). The multi-speed 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 (Mc-Namara
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 multi-speed 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
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
feq 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 Maxwell-Boltzmann 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 u2 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.
For low Mach number flows, the internal energy density equilibrium distribution
function can be further simplified by neglecting the terms O(u2)
To recover the macroscopic energy equation, the zeroth-to second-order moments
of geq must be exact. In general:
where, Ig 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
For simplicity the microscopic velocity is normalized as
Equation 62 can be calculated by using the Gauss-Hermite
quadrature. Hence, the Gauss-Hermite quadrature must consistently give accurate
result for quadratures of zeroth-to-third-order of velocity moment of geq.
This implies that the second-order Gauss-Hermite quadrature can be applied in
evaluating Ig. Therefore Eq. 62 can be evaluated
where, N is the number of abscissas, Wk is the corresponding weight
coefficient and ζk is the Gaussian abscissa. For two dimensional,
D = 2, second-order Gauss-Hermite quadrature N = 2, the value for Wk
and ζk are:
After some modification in order to satisfy macroscopic energy equation via
Chapmann-Enskog 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).
||Lattice structure for internal energy density distribution
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 so-called colored
model for immiscible two-phase 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 multi-component fluids was derived
by Shan and Chen (SC model) (1993) and later extended
by others (Shan and Chen, 1994). In the SC model, a
non-local 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 systems
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 well-established
thermodynamically. One cannot introduce temperature since the existence of any
energy-like quantity is not known (Hazi et al., 2002).
The third type of LB model for multiphase flow is based on the Free-Energy
(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
non-Galilean invariance for the viscous terms in the macroscopic Navier-Stokes
equation. Efforts have been made to restore the Galilean invariance to second-order
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 Shan-Chen which involves evolution equation of single particle distribution function f and can be written as:
To simulate multiphase fluids, long-range 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[1-exp(-ρ/ρ0)] (Raiskinmaki
et al., 2002), Ψ(ρ) = ρ , Ψρ20ρ2/[2(ρ0+ρ)2]
(Pan et al., 2004), Ψ(ρ) = Ψ0exp(-ρ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 Shan-Chen 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 cars 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 fluid-solid
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.
Earlier the formulation of mesoscale lattice Boltzmann scheme and the derivation of macroscopic Navier-Stokes equation were demonstrated. Here, the numerical prediction of lid-driven cavity flow, natural convection in an enclosure and droplet dynamics on solid surface are performed to demonstrate the applicability of the scheme.
Lid-driven cavity flow: The lid-driven 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 lid-driven 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:
|| 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.
||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:
|| Streamlines plots for (a) Ra = 103 and (b) Ra
|| Isotherms plots for (a) Ra = 103 and (b) Ra =
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 103 and 105.
The plots of streamlines and isotherms are shown in Fig. 5a,
b and 6a, b.
At Ra = 103, 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 = 105), 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.
||Computed results at 35°, 55° and 150° contact
angles. (a) Contact angle 35°, (b) Contact angle 55° and (c) Contact
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 -46-354 for contact angle
of 0°, 90° and 180°, respectively.
The simulated results at three other values of G' are shown in Fig.
7a-c. The contact angle of the computed figures are then
determined using Auto CAD software.
||Comparison of computed droplets diameter to height ratio
with analytical solution
In present investigation, we also calculate the magnitude of the ratio of droplets
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.
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 lid-driven 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, solid-fluid 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.