Journal of Applied Sciences

Year: 2010 | Volume: 10 | Issue: 15 | Page No.: 1511-1524
DOI: 10.3923/jas.2010.1511.1524

Mesoscale Numerical Approach to Predict Macroscale Fluid Flow Problems

C.S. Nor Azwadi and M.S. Idris

Abstract: We present a detailed analysis of the lattice Boltzmann method to simulate an incompressible fluid flow problem. Thorough derivation of macroscopic hydrodynamics equations from the continuous Boltzmann equation is performed. After showing how the formulation of the mesoscale particle dynamics fits in to the framework of lattice Boltzmann simulations, numerical results of isothermal, thermal and multiphase fluid flow are presented to highlight the applicability of the approach. The objective of the paper is to gain better understanding of this relatively new approach for applied engineering problems in fluid transport phenomena.

Keywords: lattice model, Distribution function, BGK collision model, boltzmann equation and computational fluid dynamics

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, 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 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²³ 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 to 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 (Wolf-Gladrow, 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 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.

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 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 Chapmann-Enskog procedure. For the rest particle, Eq. 7 becomes:

and

for other particles. This gives:

The symmetric properties of the tensor Σ_ic_σiαc_σiβ... are needed in the derivation and given as follow:

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 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(ε²) results in:

Substituting Eq. 17 and 18 into Eq. 19, the equation to order of ε and ε² 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²₃= 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²₃ρ.

We next multiply Eq. 22 with c_iσ 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 B₁ and B₂:

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 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:

Substituting Eq. 47 into Eq. 4 results in:

where

and

are the equilibrium distribution function for internal energy, the term due to the external force and the heat dissipation term respectively. 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 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 u² 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²) gives:

To recover the macroscopic energy equation, the zeroth-to second-order 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 Gauss-Hermite quadrature. Hence, the Gauss-Hermite quadrature must consistently give accurate result for quadratures of zeroth-to-third-order of velocity moment of g^eq. This implies that the second-order Gauss-Hermite 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, second-order Gauss-Hermite quadrature N = 2, the value for W_k 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).

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 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 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 Ψ(ρ) = ρ₀[1-exp(-ρ/ρ₀)] (Raiskinmaki et al., 2002), Ψ(ρ) = ρ [20], Ψρ₂₀ρ²/[2(ρ₀+ρ)²] (Pan et al., 2004), Ψ(ρ) = Ψ₀exp(-ρ₀/ρ) (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 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 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.

SIMULATION RESULTS

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:

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.

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:

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³ and 10⁵. The plots of streamlines and isotherms are shown in Fig. 5a, b and 6a, b.

At Ra = 10³, 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⁵), 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.

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.

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 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.

REFERENCES