# Flow Simulations Using Two Dimensional Thermal Lattice Boltzmann Method

Authors: Saeed J. Almalowi; Alparslan Oztekin
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

Lattice Boltzmann method is implemented to study hydrodynamically and thermally developing steady laminar flows in a channel. Numerical simulation of two-dimensional convective heat transfer problem is conducted using two-dimensional, nine directional D2Q9 thermal lattice Boltzmann arrangements. The velocity and temperature profiles in the developing region predicted by Lattice Boltzmann method are compared against those obtained by ANSYS-FLUENT. Velocity and temperature profiles as well as the skin friction and the Nusselt numbers agree very well with those predicted by the self-similar solutions of fully developed flows. It is clearly shown here that thermal lattice Boltzmann method is an effective computational fluid dynamics (CFD) tool to study nonisothermal flow problems.


---

## Body

## 1. Introduction

Historically, the Lattice Boltzman method (LBM) evolves from Lattice Gas Cellular Automata. In 1988, LBM is proposed to be used to simulate flows for the first time. The LBM is a branch of statistics of mechanics which is an ideal approach to simulate flows in simple or complex geometries. Recently, LBM has been modified to solve nonlinear partial differential equations to model complex fluid flows. Different approaches of the LBM have been discussed by several investigators [1]. However a successfully LBM simulation rests on the correct implementation of the boundary conditions, where unknown distribution function originated from the operation. As it is stated in several literatures, the implementation of the boundary conditions in LBM is the key to successfully model flow problems. Each type of boundary condition requires different technique and has different degree of accuracy as well.LBM has been used to model flows in various geometries by several research groups, but only few investigators compare LBM model against other numerical methods for computational fluid dynamics (CFD). Recently Chen et al. and Begum and Basit [2, 3] had proposed models to be applied to simulate complex flows. It has been shown that LBM can easily be implemented to study single- and multiphase flows. The equations governing conservation of mass and momentum are satisfied at each lattice nodes. The LBM approximation is a linear discretized equation which has two terms: streaming and collision terms. The collision term in the LBM has been approximated to be a linear term using models introduced by Bhatngar, Gross, and Krook (BGK) [4]. Recent study by Shi et al. [5] has shown that the LBM is a promising tool to study microscopic flows. Present work is to illustrate that LBM can be an effective CFD tool and in order to demonstrate that 2D developing nonisothermal flows in a channel is studied by implementing LBM method. The results predicted by LBM have been compared against those obtained using ANSYS-FLUENT and those obtained by self-similar solutions in the developed region for validation.


## 2. Lattice Boltzmann Governing Equation

Ludwig Eduard Boltzmann (1844–1906), the Austrian physicist, had the greatest achievement in the development of statistical mechanics. This approach has been used to predict macroscopic properties of matter such as the viscosity, thermal conductivity, and diffusion coefficient from the microscopic properties of atoms and molecules [6–8]. The probability of finding particles within certain range of velocities at a certain range of locations replaces tagging each particle as in molecular dynamics simulation. The lattice Boltzmann transportation can be governed by distribution function which represents particles at location r(x,y) at time t, and the particle will be displaced by (dx,dy) in time dt with the application force F on the liquid molecules [9]. The equation governing the distribution function f(r,c,t) has two terms, the streaming step and the collision term. Here, x and y are spatial coordinates, t is the time, and c is the lattice discrete velocity.The collision takes place between the molecules; there will be a net difference between the numbers of molecules in the intervaldrdc. The rate of change of the distribution function is expressed as(2.1)f(r+dr,c+Fdt/m,t+dt)-f(r,c,t)dt=∂f∂t+cx∂f∂x+cy∂f∂y+Fm⋅∂f∂c=φ(f).
Here, F denotes external forces applied, and φ(f) is the source or the collision term. With the absence of the external forces, (2.1) becomes(2.2)∂fi∂t+c⋅∇fi=φ(fi).
Equation (2.2) is known as the lattice Boltzmann governing equation. The right hand side of equation is called a source and is approximated by BGK as(2.3)φ(fi)=1τ(fieq-fi).
Here ω=1/τ is the relaxation frequency and the τ is the relaxation time feq is the equilibrium value of distribution function and is written as(2.4)fieq=wiρ[1+3ci⋅Vcs2+4.5(ci⋅V)2cs4-1.5V⋅Vcs2],
where ci is the discrete velocities vector, V is the bulk fluid velocity, cs is the lattice sound speed and wi is the weight factor, one has(2.5)ci={(0,0)i=1c([sin⁡(i-1)π2],[cos⁡(i-1)π2])i=2,3,4,5c(2[cos⁡(2i-11)π4],2[sin⁡(2i-11)π4])i=6,7,8,9,wi={49i=119i=2,3,4,5136i=6,7,8,9.
Equation (2.3) becomes(2.6)fi(r+dr,t+dt)=(1-ω)fi(r,t)+ωfieq.


## 3. Lattice Boltzmann Arrangements (D2Q9)

Lattice Boltzmann is relatively recent technique that has been shown to be as accurate as traditional CFD methods having ability to integrate arbitrarily complex geometries. LBM can be used for different arrangements such as D1Q2, D2Q4, D2Q9, D3Q15, D3Q19, or D3Q27 [10]. However, in this paper, we only use D2Q9 which implies the two-dimensional and nine velocity components as shown in Figure 1.Figure 1
D2Q9 LBM  arrangement.Each distribution function has position(r), velocity (c), and weight factor (w).


## 4. Momentum Lattice Boltzmann Model

The momentum LBM represents the particles velocity [4]. For instance, for D2Q9 lattice arrangements, the particle at the origin is at rest and the remaining particles move in different directions with different speed. Each velocity vector denotes a lattice per unit step. These velocities are exceptionally convenient in that all x and y components are either 0 or ±1. Mass of particle is taken as unity uniformly throughout the flow domain. The macroscopic fluid density, ρ, is governed by conservation of mass(4.1)ρ=∑i=19fi.
The bulk fluid velocity (V=(u,v)) is the average of microscopic lattice-directional velocity (c=(cx,cy)) and the directional density and is governed by conservation of momentum(4.2)V=1ρ∑i=19fici.
Here cx=dx/dt,cy=dy/dt are x and y components of the lattice directional velocity. Conservation of mass and momentum is also applied at each boundary, at the inlet and the outlet as shown in Figure 2 for the D2Q9 lattice arrangement for nodes placed on boundaries.Figure 2
D2Q9 momentum lattice arrangements at boundaries and inside the flow domain.The uniform inlet velocity isUin, the length of channel is L, and the gap between plates is H. u is measured in units of Uin(U=u/Uin),x and y are measured in units of L and H(X=x/L and Y=y/H), respectively. Scaled inlet velocity becomes U=1, and the flow domain (X,Y) becomes 0≤Y≤1, and 0≤X≤1.


## 5. Thermal Lattice Boltzmann Model

There has been rapid progress in developing the construction of stable thermal lattice Boltzmann equation models to study heat transfer problems. McNamara and Zanetti successfully applied multispeed thermal fluid lattice Boltzmann method to solve heat transfer problems [6]. At the outlet, bounce back or extrapolation boundary conditions are considered as the thermal and flow boundary conditions. Bounce back type boundary conditions are proven to provide more accurate numerical approximations [11] and are used by the present work. The temperature at each wall is specified; however, the temperatures which are pointing to the flow domain are unknowns, and they can be evaluated from streaming and collision steps. The thermal lattice arrangement is illustrated in Figure 3.Figure 3
D2Q9 thermal lattice arrangements at boundaries and inside the flow domain.The rate of change of the thermal distribution function is written as(5.1)gi(r+dr,t+dt)=(1-ω)gi(r,t)+ωgieq.
With normalized temperature θ=(T-Tw)/[Tin(H/2)-Tw], the equilibrium value of the thermal distribution function geq is given by(5.2)gieq=ωiθ[1+3ci⋅Vcs2+4.5(ci⋅V)2cs4-1.5V⋅Vcs2].
For simplicity the relaxation frequency of the thermal and momentum distribution function is selected as the same. Hence the kinematic viscosity, υ, and the thermal diffusivity, α, are the same and are expressed by(5.3)v=α=(dx)23dt(1ω-12).
Here, Tw is the wall temperature, and Tin(y) is the temperature distribution at the inlet. The Prandtl number Pr=υ/α=1. The temperature of the fluid is governed by the conservation of energy(5.4)T=∑i=19gi.


## 6. Results and Discussion

The discretized equations (2.6) and (5.1) for momentum and thermal distribution f and g for M×N nodes in x and y directions, respectively, are solved by employing Gauss-Siedel iterations. For boundary nodes the detailed discretized equations for f and g are described in Appendices A and B. The results are presented for steady incompressible two-dimensional laminar flows in an entrance region of a channel. Flow develops hydrodynamically and thermally in a 1 m long and 0.02 m height channel with the aspect ratio AR of 50. At the inlet the flow is uniform (Uin=0.02m/sec), and the temperature of the fluid satisfies θin=4×(Y-Y2). Boundary conditions imposed on the velocity field at Y=0 and 1 are no-slip and no-penetration, and the thermal boundary conditions applied on each surface are θ=0. The physical properties are considered to be constant and are determined for water at 300 K—(ρ=999.1kg/m3 and μ=855×10-3N·s/m2). For the example illustrated in this paper, the flow rate considered is 0.4 kg/s and the corresponding Reynolds number Re=ρ2HUin/μ=800.Spectral convergence is checked for LBM to ensure that the results predicted by the LBM are not dependent on the number of nodes selected for the numerical simulations. Nodes(M×N) are placed uniformly in the direction of x(M nodes) and y(N nodes). The convergence test is displayed in Figure 4 as the velocity and temperature profile at X=0.2 plotted for various M and N. It is shown that the (50×1000) mesh provides satisfactory spectral convergence and numerical accuracy and is thereby chosen for the numerical simulation results predicted by LBM.(a) Velocity and (b) temperature profiles atX=0.2 plotted for various N and M.
(a)(b)The velocity and temperature profiles are displayed at various cross-sections in the developing region. The profiles that are obtained by LBM are compared against those obtained by ANSYS-FLUENT at the same conditions. The boundary conditions at the inlet and the outlet and on the surface of the plates are selected as the same for both methods. The velocity and temperature field are considered to be converged when the error tolerance is less than 10−3.The velocity profiles predicted by LBM and FLUENT at various cross-sections are shown in Figure5. The solid lines denote the prediction obtained by LBM while the symbols denote the predictions obtained by FLUENT. The velocity profiles at all cross-section predicted by LBM agree very well with those predicted by FLUENT, as shown in Figure 5(a). The development length for velocity field at Re=800 in this flow is expected to be x/H=47.4. The thermal field has the same development length as the hydrodynamic field since Pr is selected to be unity. The nearly fully developed velocity profile obtained by both methods at X=0.725 also agrees very well with each other. They also agree well with the analytical solution, U=6(Y-Y2), obtained by the self-similar solution for fully developed laminar flow in a channel, as shown in Figure 5(b).(a) Velocity profiles at various cross-sections predicted by LBM and FLUENT. (b) Velocity profile atX=0.725 predicted by LBM and FLUENT and the self-similar solution for the fully developed laminar channel flow.
(a)(b)The temperature profiles predicted by LBM and FLUENT at various cross-sections are shown in Figure6. The solid lines denote the prediction obtained by LBM, and the symbols denote the predictions obtained by FLUENT. The temperature profiles predicted by LBM agree very well with those predicted by FLUENT.Figure 6
Temperature profile at various cross-sections predicted by LBM and FLUENT.Wall shear stress and heat transfer coefficient are predicted at various cross-sections in the developing region. The local value of skin friction,Cf, and the Nusselt number, Nu, are determined from the numerical solution as(6.1)Cf(X)=μ1/2ρUin2∂u∂y(x,0)=4Re∂U∂Y(X,0),Nu(X)=2Tin-TwTm-Tw∂θ∂Y(X,0),
where Tm is the bulk temperature of the fluid and is calculated at each cross-section as(6.2)Tm=1∫UdA∫[Tw+(Tin(0.5H)-Tw)θ]UdA.
The skin friction and the Nusselt number predicted by LBM are plotted in Figure 7 as a function x/2H. ReCf tends to 24 as the fully developed region is approached. Similarly, Nusselt number tends to 7.54 as the thermally fully developed region is approached, as shown in Figure 7(b). These values are in perfect agreement with the fully developed values of Cf and Nu as well documented in the literature.(a) Skin friction and (b) Nusselt number plotted as a function ofx/2H.
(a)(b)


## 7. Conclusion

Hydrodynamically and thermally developing laminar steady flow in a channel is considered as an example to illustrate that Lattice Boltzmann method is a promising computational fluid dynamics tool. D2Q9 lattice arrangement is used to predict both velocity and temperature field. Profiles obtained by LBM-D2Q9 and ANSYS-FLUENT agree very well. Away from the inlet as the fully developed region is approached, and the profiles tend to self-similar solutions of developed laminar channel flows. That is confirmed by prediction of the skin friction and Nusselts number. As full-developed region is approached away from the inlet both the skin friction coefficient and the Nusselt number tend to well-documented values in the literature. Extension of LBM method to three-dimensional unsteady complex multiphase flows is natural, and the implementation of LBM to tackle such problems is underway.



---

*Source: 135173-2012-06-13.xml*