Flow Simulations Using Two Dimensional Thermal Lattice Boltzmann Method

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.

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 𝐫(𝑥,𝑦) at time 𝑡, and the particle will be displaced by (𝑑𝑥,𝑑𝑦) in time 𝑑𝑡 with the application force 𝐅 on the liquid molecules [9]. The equation governing the distribution function 𝑓(𝐫,𝐜,𝑡) has two terms, the streaming step and the collision term. Here, 𝑥 and 𝑦 are spatial coordinates, 𝑡 is the time, and 𝐜 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 interval 𝑑𝐫𝑑𝐜. The rate of change of the distribution function is expressed as𝑓(𝐫+𝑑𝐫,𝐜+𝐅𝑑𝑡/𝑚,𝑡+𝑑𝑡)−𝑓(𝐫,𝐜,𝑡)=𝑑𝑡𝜕𝑓𝜕𝑡+𝑐_𝑥𝜕𝑓𝜕𝑥+𝑐_𝑦𝜕𝑓+𝐅𝜕𝑦𝑚⋅𝜕𝑓𝑓.𝜕𝐜=𝜑 (2.1) Here, 𝐅 denotes external forces applied, and 𝜑(𝑓) is the source or the collision term. With the absence of the external forces, (2.1) becomes𝜕𝑓_𝑖𝜕𝑡+𝐜⋅∇𝑓_𝑖=𝜑𝑓_𝑖. (2.2) 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𝜑𝑓_𝑖=1𝜏𝑓^eq_𝑖−𝑓_𝑖. (2.3) Here 𝜔=1/𝜏 is the relaxation frequency and the 𝜏 is the relaxation time 𝑓^eq is the equilibrium value of distribution function and is written as𝑓^eq_𝑖=𝑤_𝑖𝜌1+3𝐜_𝑖⋅𝐕𝑐_𝑠²𝐜+4.5^𝑖⋅𝐕²𝑐_𝑠⁴−1.5𝐕⋅𝐕𝑐_𝑠², (2.4) where 𝐜_𝑖 is the discrete velocities vector, 𝐕 is the bulk fluid velocity, 𝑐_𝑠 is the lattice sound speed and 𝑤_𝑖 is the weight factor, one has𝐜_𝑖=()𝐜()𝜋0,0𝑖=1sin𝑖−12,cos(𝑖−1)𝜋2𝐜√𝑖=2,3,4,52cos(2𝑖−11)𝜋4,√2sin(2𝑖−11)𝜋4𝑤𝑖=6,7,8,9,_𝑖=491𝑖=191𝑖=2,3,4,536𝑖=6,7,8,9. (2.5) Equation (2.3) becomes𝑓_𝑖(𝐫+𝑑𝐫,𝑡+𝑑𝑡)=(1−𝜔)𝑓_𝑖(𝐫,𝑡)+𝜔𝑓^eq_𝑖. (2.6)

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 (𝐫), velocity (𝐜), and weight factor (𝑤).

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 𝑥 and 𝑦 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𝜌=⁹_𝑖=1𝑓_𝑖. (4.1) The bulk fluid velocity (𝐕=(𝑢,𝑣)) is the average of microscopic lattice-directional velocity (𝐜=(𝑐_𝑥,𝑐_𝑦)) and the directional density and is governed by conservation of momentum1𝐕=𝜌^9𝑖=1𝑓_𝑖𝐜_𝑖. (4.2) Here 𝑐_𝑥=𝑑𝑥/𝑑𝑡,𝑐_𝑦=𝑑𝑦/𝑑𝑡 are 𝑥 and 𝑦 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 is 𝑈_in, the length of channel is 𝐿, and the gap between plates is 𝐻. 𝑢 is measured in units of 𝑈_in(𝑈=𝑢/𝑈_in),𝑥 and 𝑦 are measured in units of 𝐿 and 𝐻(𝑋=𝑥/𝐿 and 𝑌=𝑦/𝐻), respectively. Scaled inlet velocity becomes 𝑈=1, and the flow domain (𝑋,𝑌) becomes 0≤𝑌≤1, and 0≤𝑋≤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𝑔_𝑖(𝐫+𝑑𝐫,𝑡+𝑑𝑡)=(1−𝜔)𝑔_𝑖(𝐫,𝑡)+𝜔𝑔^eq_𝑖. (5.1) With normalized temperature 𝜃=(𝑇−𝑇_𝑤)/[𝑇_in(𝐻/2)−𝑇_𝑤], the equilibrium value of the thermal distribution function 𝑔^eq is given by𝑔^eq_𝑖=𝜔_𝑖𝜃1+3𝐜_𝑖⋅𝐕𝑐_𝑠²𝐜+4.5^𝑖⋅𝐕²𝑐_𝑠⁴−1.5𝐕⋅𝐕𝑐_𝑠². (5.2) 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()𝑣=𝛼=𝑑𝑥²13𝑑𝑡𝜔−12. (5.3) Here, 𝑇_𝑤 is the wall temperature, and 𝑇_in(𝑦) is the temperature distribution at the inlet. The Prandtl number Pr=𝜐/𝛼=1. The temperature of the fluid is governed by the conservation of energy𝑇=⁹_𝑖=1𝑔_𝑖. (5.4)

6. Results and Discussion

The discretized equations (2.6) and (5.1) for momentum and thermal distribution 𝑓 and 𝑔 for 𝑀×𝑁 nodes in 𝑥 and 𝑦 directions, respectively, are solved by employing Gauss-Siedel iterations. For boundary nodes the detailed discretized equations for 𝑓 and 𝑔 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 (𝑈_in=0.02m/sec), and the temperature of the fluid satisfies 𝜃_in=4×(𝑌−𝑌²). Boundary conditions imposed on the velocity field at 𝑌=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/m³ and 𝜇=855×10⁻³N⋅s/m²). For the example illustrated in this paper, the flow rate considered is 0.4?kg/s and the corresponding Reynolds number Re=𝜌2𝐻𝑈_in/𝜇=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 (𝑀×𝑁) are placed uniformly in the direction of 𝑥(𝑀 nodes) and 𝑦(𝑁 nodes). The convergence test is displayed in Figure 4 as the velocity and temperature profile at 𝑋=0.2 plotted for various 𝑀 and 𝑁. 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)

(b)

Figure 4 

(a) Velocity and (b) temperature profiles at 𝑋=0.2 plotted for various 𝑁 and 𝑀.

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 Figure 5. 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 𝑥/𝐻=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 𝑋=0.725 also agrees very well with each other. They also agree well with the analytical solution, 𝑈=6(𝑌−𝑌²), obtained by the self-similar solution for fully developed laminar flow in a channel, as shown in Figure 5(b).

(a)

(b)

Figure 5 

(a) Velocity profiles at various cross-sections predicted by LBM and FLUENT. (b) Velocity profile at 𝑋=0.725 predicted by LBM and FLUENT and the self-similar solution for the fully developed laminar channel flow.

The temperature profiles predicted by LBM and FLUENT at various cross-sections are shown in Figure 6. 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, 𝐶_𝑓, and the Nusselt number, Nu, are determined from the numerical solution as𝐶_𝑓(𝑋)=𝜇1/2𝜌𝑈²_in𝜕𝑢()=4𝜕𝑦𝑥,0Re𝜕𝑈()(𝑋)𝑇𝜕𝑌𝑋,0,Nu=2ⁱⁿ−𝑇_𝑤𝑇_𝑚−𝑇_𝑤𝜕𝜃(),𝜕𝑌𝑋,0 (6.1) where 𝑇_𝑚 is the bulk temperature of the fluid and is calculated at each cross-section as𝑇_𝑚=1𝑇∫𝑈𝑑𝐴^𝑤+𝑇_in(0.5𝐻)−𝑇_𝑤𝜃𝑈𝑑𝐴. (6.2) The skin friction and the Nusselt number predicted by LBM are plotted in Figure 7 as a function 𝑥/2𝐻. Re𝐶_𝑓 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 𝐶_𝑓 and Nu as well documented in the literature.

(a)

(b)

Figure 7 

(a) Skin friction and (b) Nusselt number plotted as a function of 𝑥/2𝐻.

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.

Appendices

A. Discretized LB Equations for the Velocity Field at the Inlet, the Outlet, and at the Surface of Plates

With indices 𝑘 denoting 9 directions of LB D2Q9 arrangements, 𝑖 denoting the nodes placed in the 𝑥-direction and 𝑗 denoting the nodes placed in the 𝑦-directions, the distribution functions for velocity, 𝑓, and temperature field, 𝑔, are represented by three-dimensional arrays 𝑓(𝑘,𝑖,𝑗), and 𝑔(𝑘,𝑖,𝑗) for 𝑘=1 to 9, 𝑖=1,𝑀 and 𝑗=1,𝑁.

At the inlet (𝑖=1;𝑗=2 to 𝑁−1), the conservation of mass and momentum yield 𝜌=𝑓1,𝑖,𝑗(𝑓1,𝑖,𝑗+𝑓5,𝑖,𝑗+𝑓3,𝑖,𝑗+2×4,𝑖,𝑗+𝑓7,𝑖,𝑗+𝑓8,𝑖,𝑗,𝑓1−𝑢1,𝑖,𝑗2,𝑖,𝑗−𝑓^eq2,𝑖,𝑗=𝑓4,𝑖,𝑗−𝑓^eq,𝑓+24,𝑖,𝑗2,𝑖,𝑗=𝑓4,𝑖,𝑗3𝜌,𝑓+11,𝑖,𝑗×𝑢1,𝑖,𝑗9,𝑖,𝑗=𝑓6,𝑖,𝑗6𝜌,𝑓+11,𝑖,𝑗×𝑢1,𝑖,𝑗5,𝑖,𝑗=𝑓7,𝑖,𝑗6𝜌1,𝑖,𝑗×𝑢1,𝑖,𝑗. (A.1) The boundary conditions at the surface of the lower plate (𝑖=1,𝑀−1;𝑗=1) imposed on the velocity field give𝜌=𝑓1,𝑖,𝑗(𝑓1,𝑖,𝑗+𝑓5,𝑖,𝑗+𝑓2,𝑖,𝑗+2×3,𝑖,𝑗+𝑓9,𝑖,𝑗+𝑓8,𝑖,𝑗,𝑓+21−𝑣1,𝑖,𝑗5,𝑖,𝑗=𝑓3,𝑖,𝑗3𝜌1,𝑖,𝑗×𝑣1,𝑖,𝑗,𝑓3,𝑖,𝑗−𝑓^eq3,𝑖,𝑗=𝑓5,𝑖,𝑗−𝑓^eq,𝑓+15,𝑖,𝑗6,𝑖,𝑗=𝑓8,𝑖,𝑗2𝑓+14,𝑖,𝑗−𝑓2,𝑖,𝑗6𝜌+11,𝑖,𝑗×𝑣1,𝑖,𝑗2𝜌,𝑓1,𝑖,𝑗×𝑢1,𝑖,𝑗₇=𝑓₉−12𝑓₄−𝑓₂+16𝜌−11,𝑖,𝑗×𝑣1,𝑖,𝑗2𝜌1,𝑖,𝑗×𝑢1,𝑖,𝑗. (A.2) At the top surface (𝑖=1,𝑀−1;𝑗=𝑁), the velocity boundary conditions yield𝜌=𝑓1,𝑖,𝑗(𝑓1,𝑖,𝑗+𝑓2,𝑖,𝑗+𝑓4,𝑖,𝑗+2×5,𝑖,𝑗+𝑓7,𝑖,𝑗+𝑓6,𝑖,𝑗,𝑓−21+𝑣1,𝑖,𝑗3,𝑖,𝑗=𝑓5,𝑖,𝑗3𝜌_𝑢,𝑓−11,𝑖,𝑗×𝑣1,𝑖,𝑗8,𝑖,𝑗=𝑓6,𝑖,𝑗2𝑓−14,𝑖,𝑗−𝑓2,𝑖,𝑗6𝜌−11,𝑖,𝑗×𝑣1,𝑖,𝑗2𝜌,𝑓+11,𝑖,𝑗×𝑢1,𝑖,𝑗9,𝑖,𝑗=𝑓7,𝑖,𝑗2𝑓−14,𝑖,𝑗−𝑓2,𝑖,𝑗6𝜌+11,𝑖,𝑗×𝑣1,𝑖,𝑗2𝜌1,𝑖,𝑗×𝑢1,𝑖,𝑗. (A.3) The outlet conservation of mass and momentum is approximated by using second-order extrapolation which yields𝑓𝑘,𝑀,𝑗=2𝑓𝑘,𝑀−1,𝑗−𝑓𝑘,𝑀−2,𝑗where𝑘=2,6,9,𝑗=2,𝑁−1. (A.4)

B. Discretized LB Equations for the Temperature Field at the Inlet, the Outlet, and at the Surface of Plates

The inlet (𝑖=1;𝑗=2 to 𝑁−1) conservation of energy gives𝑔3,𝑖,𝑗−𝑔^eq3,𝑖,𝑗=𝑔5,𝑖,𝑗−𝑔^eq,𝑔5,𝑖,𝑗2,𝑖,𝑗=𝜃_in(𝑤(4)(2,𝑔1,𝑖,𝑗+𝑤))−𝑔4,𝑖,𝑗6,𝑖,𝑗=𝜃_in(𝑤(6)(8,𝑔1,𝑖,𝑗+𝑤))−𝑔8,𝑖,𝑗9,𝑖,𝑗=𝜃_in1,𝑖,𝑗(𝑤(7)+𝑤(9))−𝑔7,𝑖,𝑗. (B.1) With 𝜃_𝑤=0 the thermal boundary condition at the surface of the lower plate (𝑖=1,𝑀−1;𝑗=1) gives𝑔,𝑔,𝑔,𝑔.7,𝑖,𝑗=−𝑔7,𝑖,𝑗5,𝑖,𝑗=−𝑔3,𝑖,𝑗6,𝑖,𝑗=−𝑔4,𝑖,𝑗2,𝑖,𝑗=−𝑔4,𝑖,𝑗 (B.2) With 𝜃_𝑤=0 the thermal boundary condition at the surface of the upper plate (𝑖=1,𝑀−1;𝑗=𝑁) gives𝑔3,𝑖,𝑗−𝑔^eq3,𝑖,𝑗=𝑔5,𝑖,𝑗−𝑔^eq,𝑔,𝑔,𝑔,𝑔.5,𝑖,𝑗9,𝑖,𝑗=−𝑔7,𝑖,𝑗8,𝑖,𝑗=−𝑔6,𝑖,𝑗3,𝑖,𝑗=−𝑔5,𝑖,𝑗2.𝑖.𝑗=−𝑔4,𝑖,𝑗 (B.3) The outlet conservation of energy is approximated by the second-order extrapolation which yields𝑔𝑘,𝑀,𝑗=2𝑔𝑘,𝑀−1,𝑗−𝑔𝑘,𝑀−2,𝑗where𝑘=2,6,9,𝑗=2,𝑁−1. (B.4)

Nomenclature

𝑓:	Density distribution function
𝑓eq:	Local equilibrium density distribution function
𝑔:	Temperature distribution function
𝑔eq:	Local equilibrium temperature distribution function
𝐜=(𝑐𝑥,𝑐𝑦):	Lattice discrete velocity, 𝑐𝑥 and 𝑐𝑦 are 𝑥 and 𝑦 components
𝐕=(𝑢,𝑣):	Bulk velocity of the fluid, 𝑢 and 𝑣 are 𝑥 and 𝑦 components
𝑈:	Normalized 𝑥 component of the fluid velocity
𝜃:	Normalized temperature
𝜔:	Dimensionless relaxation frequency
𝑤:	Weight factor
𝑇𝑤:	Wall temperature in (°C)
𝑐𝑠:	Lattice sound speed
𝜏:	Dimensionless collision relaxation time
𝐫:	Position vector
(𝑋,𝑌):	Dimensionless 𝑥 and 𝑦 coordinate
𝑈in:	Fluid speed at the inlet
𝑇in:	Temperature profile at the inlet in (°C), 𝑇in(𝑦)
𝜌:	Density of fluid, kg/m3
𝜐:	Kinematic viscosity of fluid, m2/sec
𝛼:	Thermal diffusivity of fluid, m2/sec
Pr:	Prandtl number
Re:	Reynolds number
𝑇:	Bulk temperature in (°C).

Acknowledgment

One of the authors of this paper is grateful to Saudi Arabia government for its support toward his graduate education.