A direct-forcing immersed boundary method for the thermal lattice Boltzmann method

A direct-forcing immersed boundary method for the thermal lattice Boltzmann method

Computers & Fluids 49 (2011) 36–45 Contents lists available at ScienceDirect Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i ...

1MB Sizes 0 Downloads 89 Views

Computers & Fluids 49 (2011) 36–45

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

A direct-forcing immersed boundary method for the thermal lattice Boltzmann method Shin K. Kang ⇑, Yassin A. Hassan Department of Nuclear Engineering, Texas A&M University, College Station, TX 77843-3133, USA

a r t i c l e

i n f o

Article history: Received 3 July 2010 Received in revised form 26 March 2011 Accepted 20 April 2011 Available online 20 May 2011 Keywords: Direct-forcing immersed boundary method Thermal lattice Boltzmann method Double-population model Hybrid thermal model Natural convection Moving particle

a b s t r a c t In this study, a direct-forcing immersed boundary method (IBM) for thermal lattice Boltzmann method (TLBM) is proposed to simulate the non-isothermal flows. The direct-forcing IBM formulas for thermal equations are derived based on two TLBM models: a double-population model with a simplified thermal lattice Boltzmann equation (Model 1) and a hybrid model with an advection–diffusion equation of temperature (Model 2). As an interface scheme, which is required due to a mismatch between boundary and computational grids in the IBM, the sharp interface scheme based on second-order bilinear and linear interpolations (instead of the diffuse interface scheme, which uses discrete delta functions) is adopted to obtain the more accurate results. The proposed methods are validated through convective heat transfer problems with not only stationary but also moving boundaries – the natural convection in a square cavity with an eccentrically located cylinder and a cold particle sedimentation in an infinite channel. In terms of accuracy, the results from the IBM based on both models are comparable and show a good agreement with those from other numerical methods. In contrast, the IBM based on Model 2 is more numerically efficient than the IBM based on Model 1. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction The immersed boundary method (IBM) has been widely applied to complex and moving-boundary problems. The IBM is based on the structured (usually Cartesian), non-body conformal grid, thus reducing the computational cost and relieving the burden of meshing in complex geometry. The IBM is more beneficial to modeling moving-boundary problems because time-consuming re-meshing in the body-conformal grid methods is no longer required in the fixed, non-body-conformal grids in the IBM. However, since the grids are not fitted to the boundary in the IBM, the boundary effect, i.e. no-slip conditions, has to be modeled by forcing terms in the momentum equations. For the evaluation of the forcing term, the direct-forcing IBM [1,2] has been prevalently used because the forcing term is naturally evaluated in the conventional calculation process without a stability problem and without arbitrary parameters to be adjusted as in the feedback-forcing IBM [3,4]. A velocity included in the direct-forcing IBM formula, which is not directly given due to mismatch between the boundary and the computational nodes, should be evaluated by adopting a suitable interface (interpolation) scheme [5,6]. In recent years, the concept of forcing in momentum equations under the IBM has been also extended to the energy equation in or-

⇑ Corresponding author. Tel.: +1 979 845 7090; fax: +1 979 845 6443. E-mail address: [email protected] (S.K. Kang). 0045-7930/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2011.04.016

der to satisfy the thermal boundary conditions [7–10]. In this paper, we call the forcing terms in the momentum equations for no-slip boundary conditions ‘‘momentum-forcing’’ and the energy source term in the energy equation for thermal boundary conditions ‘‘energy-forcing’’ [8]. The principles of forcing term evaluation in both momentum and energy equations are basically the same; the only difference is that the no-slip boundary condition is only considered in the momentum-forcing IBM, while more various boundary conditions, such as Dirichlet and Neumann types, can be involved in the energy-forcing IBM. On the other hand, the IBM has also been coupled with the lattice Boltzmann method (LBM) – which is called the immersed boundary-lattice Boltzmann method (IB-LBM) – instead of Navier–Stokes equations (NSE), because the IBM and the LBM have the common feature of using Cartesian grids. In addition, the LBM is simple to implement and easy to parallelize compared to the Navier–Stokes equation method because it does not require the iterative procedure for the pressure–velocity coupling [11,12]. Furthermore, in the direct-forcing IB-LBM based on the split-forcing lattice Boltzmann equation (LBE), the boundary forcing term can be easily and accurately evaluated [13]. It should be noted that for complex boundary problems, the conventional LBM imposes the no-slip boundary condition in mesoscopic level by using the density distribution function interpolation, while the IB-LBM based on the sharp interface scheme imposes the no-slip condition in macroscopic level by using the velocity interpolation. Thus, the IB-LBM with sharp interface scheme has an


S.K. Kang, Y.A. Hassan / Computers & Fluids 49 (2011) 36–45

advantage of involving less number of interpolations than the conventional LBM with interpolation treatments. For a complex 2D geometry, the IB-LBM involves the maximum two-directional (x and y) interpolations, whereas the conventional LBM requires the maximum five-directional interpolations. The difference in the number of interpolations involved can increase in the 3D case. Although numerous IB-LBMs have been applied to isothermal flows, the coupling between the IBM and thermal lattice Boltzmann equation (TLBE) for non-isothermal flows was not documented until the recent work of Jeong et al. [14]. They called the method the immersed boundary-thermal lattice Boltzmann method (IB-TLBM). However, they adopted the feedback-forcing method to evaluate the momentum and energy forcing, thus causing the stability problem and arbitrariness in selecting the parameter. In addition, their method is based on the double-population model with a complex TLBE [15]. TLBMs can be classified into three categories: the single-population model, the double-population model, and the hybrid model. In the single-population model (or multi-speed model), only the density distribution function was used, but additional discrete velocities were introduced to obtain the energy equation, and the equilibrium distribution functions usually include higher order velocity terms [16,17]. Earlier single-population models suffered from a severe numerical instability and the range of temperature variation was limited [18]. However, recently, several improved single-population models to overcome such problems have been developed [19–21]. On the other hand, in the double-population models, distribution functions for temperature (or internal energy) were introduced in addition to the original density-distribution function, so that the athermal LBE of density distribution functions for momentum and the thermal LBE of temperature (or internal energy) distribution functions for energy were separately solved. As a result, this kind of model could effectively overcome two limitations of the multi-speed models; namely, severe numerical instability and narrow range of temperature variation [15]. The correct double-population model, including terms of the viscous heat dissipation and the compression work done by pressure, was derived by He et al. [15]. However, this method is too complicated to use because it contains a complex gradient operator term in the thermal LBE. Hence, several double-population models with a simplified thermal LBE have been proposed [22–24] and applied to various heat transfer problems successfully [25–29]. The simplified thermal LBEs corresponded to the energy equation without terms of the viscous dissipation and compression work done by pressure, i.e. advection–diffusion equation of temperature. However, even double-population models with such simplified TLBEs are still numerically inefficient because they utilize a full set of distribution functions to calculate the temperature, although a reduced set of distribution functions [23,30] can slightly improve the numerical inefficiency [31]. Thus, the hybrid TLBM [31–34] was proposed in which the mass and momentum conservation are solved by the usual athermal LBE, while the advection–diffusion equation satisfied by the temperature is solved separately by a finite difference technique. The hybrid methods could effectively overcome both the instability of single-population models and the numerical inefficiency of double-population models. The objective of this study is to extend the direct-forcing IBM for the isothermal LBM to the thermal LBM. Our strategy is as follows. We consider two thermal LBM models: the double-population model with a simplified thermal LBE, and the hybrid model with an advection–diffusion equation of temperature. We introduce the energy source term for each thermal equation and then derive the direct-forcing IBM formula for both equations, as done in the isothermal IBM. To obtain more accurate boundary effects,

we adopt the sharp interface scheme based on bilinear and linear interpolations instead of the diffuse interface scheme based on discrete delta functions. The proposed methods are tested through convective heat transfer problems with not only stationary but also moving boundaries. The remaining part of this paper is organized as follows. The direct-forcing IBM for the isothermal LBE is briefly reviewed in Section 2.1. Then, the direct-forcing method is extended to a simplified internal energy LBE with an energy source term under the double-population TLBM model in Section 2.2 and then to an advection–diffusion equation of temperature under the frame of the hybrid TLBM model in Section 2.3. In Section 2.4, Newton’s equations of motion are discussed for the moving particle application. To validate the LBM code for two TLBM models, the natural convection in a square cavity is first considered in Section 3.1. Then, the present IBM is applied to a complex, fixed-boundary heat transfer problem (the natural convection in a square cavity with an eccentrically located cylinder) in Section 3.2, and then to a complex, moving-boundary heat transfer problem (a cold particle sedimentation in a hot infinite channel) in Section 3.3. In Section 4, summary and conclusions are addressed. 2. Numerical methods 2.1. Immersed boundary method for the isothermal lattice Boltzmann equation The single-relaxation-time (SRT) LBE with an unsteady, nonuniform, body-force density term can be expressed as [35]:

fi ðx þ ei Dt; t þ DtÞ ¼ fi ðx; tÞ 




fi ðx; tÞ  fi

 i  1 F i Dt ðx; tÞ þ 1  2s ð1Þ

where fi(x, t) is the density distribution function for the discrete velocity ei and s is dimensionless relaxation time. The equilibrium ðeqÞ distribution function fi is given as ðeqÞ


  3 4:5 1:5 ¼ wi q 1 þ 2 ðei  uÞ þ 2 ðei  uÞ2  2 u2 c c c


where the lattice speed c = Dx/Dt, and Dx and Dt are the lattice constant and the time step size, respectively, and weighting coefficient wi depends on the discrete velocity set. The nine-velocity LBE model on a square lattice, denoted as D2Q9 model, is commonly used in two dimensional flows. In D2Q9 model, discrete velocity vectors are defined as

8 0; i¼0 > > h i h i > <  ði1Þp ði1Þp ; i ¼ 1; 2; 3; 4 ei ¼ c cos 2 ; sin 2 > h i h i pffiffiffi  > > ði1Þ p ði1Þ p : 2c cos ; sin 2 ; i ¼ 5; 6; 7; 8 2


and the corresponding weighting coefficients wi are

8 i¼0 > < 4=9; wi ¼ 1=9; i ¼ 1; 2; 3; 4 > : 1=36; i ¼ 5; 6; 7; 8


The discrete forcing term Fi(x,t) is defined as:

  ei  uðx; tÞ ei  uðx; tÞ F i ðx; tÞ ¼ wi 3 þ9 ei  Fðx; tÞ 2 4 c c


where F is the force density. The fluid velocity u is determined by

qu ¼

X i

ei fi þ

Dt F 2



S.K. Kang, Y.A. Hassan / Computers & Fluids 49 (2011) 36–45

Using the Chapman-Enskog multi-scale analysis, we can show that this LBE correctly recovers the following Navier–Stokes equation [35]:

@q þ r  qu ¼ 0 @t @ ðquÞ þ r  ðquuÞ ¼ rp þ mr  ½qru þ ðruÞT  þ F @t

Ud ¼ ð7Þ ð8Þ

where the viscosity is defined by

1 3


 1 2 c Dt: 2

to the boundary. To evaluate Ud on the forcing node, the following bilinear interpolation is basically used:


In the actual calculation, the LBE is solved through the following four steps.

1 fub  ½Dx ð1  Dy Þu2 þ ð1  Dx Þð1  Dy Þu3 þ ð1  Dx ÞDy u4 g Dx Dy ð15Þ

where subscript b indicates the point on the boundary closest to the forcing node and subscripts 2, 3, and 4 are fluid nodes translated by one grid in y-direction, x-direction, and both directions, and Dx and Dy are distances between points b and 2 in x-direction, and between points b and 4 in y-direction. However, in the case that one of neighboring three points (excluding the boundary point) is another forcing node, the following linear interpolation is employed:

( d

U ¼


ei fi þ


Dt F 2


 1DD uf

if D P 12


if D < 12

2ub  2Duf  ð1  2DÞuff

First-forcing step:

qu ¼

1 u D b

where f and ff indicate the fluid nodes translated from the forcing nodes by one and two grid sizes, either x- or y-direction, and D is the distance between the boundary point b and the fluid node f.

Collision step:

1 ðeqÞ fi0 ðx; tÞ ¼ fi ðx; tÞ  ½fi ðx; tÞ  fi ðq; uÞ



2.2. Immersed boundary method for the thermal lattice Boltzmann equation The simplified TLBE can be expressed as [22]:

Second-forcing step:

  1 fi00 ðx; tÞ ¼ fi0 ðx; tÞ þ 1  DtF i 2sg


g i ðx þ ei Dt; t þ DtÞ ¼ g i ðx; tÞ 

ð10dÞ ðeqÞ

For the simulation of the non-isothermal flows with significant buoyancy force effect, we adopt the Boussinesq approximation and then the body force term in Eq. (8) becomes

F ¼ qbðT  T 0 Þg


where T0 is the average temperature and b is the thermal expansion coefficient at T0. The direct forcing formula in the LBE can be directly derived [13]. Suppose that x is the forcing node. After the streaming step, the velocity under no external force, unoF, at (x, t + Dt) can be expressed as:

qðx; t þ DtÞunoF ðx; t þ DtÞ ¼


ei fi ðx; t þ DtÞ:



If the desired velocity, U , which satisfy the no-slip condition on the boundary, is given, then from Eq. (6)




8 2 wi qeð1:5 uc2 Þ; i¼0 > > < 2 ei u ðei uÞ2 u ¼ wi qeð1:5 þ 1:5 c2 þ 4:5 c4  1:5 c2 Þ; i ¼ 1; 2; 3; 4 > > 2 : 2 wi qeð3 þ 6 eci 2u þ 4:5 ðeicuÞ  1:5 uc2 Þ; i ¼ 5; 6; 7; 8 4 ð18Þ 2

where e = RT with the gas constant R and c = 3RT0 with the mean temperature T0. The energy distribution functions satisfy the following condition:


g i ¼ qe:


In the IBM, we use the energy source term to model the thermal boundary effect. Thus, we need to additionally consider the energy source term to the simplified TLBE. The TLBE with an energy source term can be expressed as [26]:


qðx; t þ DtÞU ¼

ðx; tÞ ðeqÞ

fi ðx þ ei Dt; t þ DtÞ ¼ fi00 ðx; tÞ



½g i ðx; tÞ  g i

where gi is the energy distribution function and g i is the equilibrium energy distribution function, which is determined by:

Streaming step:




Dt ei fi ðx; t þ DtÞ þ Fb ðx; t þ DtÞ: 2


Subtracting Eq. (11) from Eq. (11), we can obtain the following direct-forcing formula: d

Fb ðx; t þ DtÞ ¼ 2qðx; t þ DtÞ

U u


ðx; t þ DtÞ : Dt

g i ðx þ ei Dt; t þ DtÞ ¼ g i ðx; tÞ 


g i ðx; tÞ  g i

  1 Qi þ Dt 1  2sg

i ðx; tÞ ð19Þ

where the discrete energy source function is

Q i ¼ wi Q


and Q is the energy source density term. The macroscopic internal energy can be calculated by


However, since the computational nodes, in general, do not match the boundary in complex boundary problems, we need to adopt the interface scheme. For more accurate boundary treatment, instead of using the diffuse interface scheme based on the discrete delta function, we employ the exterior sharp interface scheme [2,13] where the forcing nodes are located inside the solid closest



qe ¼

X i

gi þ

Dt Q: 2


Applying the Chapman–Enskog multi-scale analysis, we can show that Eq. (5) recovers the following energy equation:

@ ðqeÞ þ r  ðqeuÞ ¼ ar2 ðqeÞ þ Q @t



S.K. Kang, Y.A. Hassan / Computers & Fluids 49 (2011) 36–45

where the thermal diffusivity a is expressed as

2 a¼ 3


 1 2 c Dt:  2



It should be noted that in Eqs. (19) and (22), the compressible work and the viscous heat dissipation terms are neglected. In the same manner as the IB-LBM, we can apply the IBM to the non-isothermal lattice Boltzmann equation. After the streaming step, the temperature under no external energy source, TnoE, at (x, t + Dt) can be expressed as: 2

qðx; t þ DtÞ

c T noE ðx; t þ DtÞ ¼ 3T 0


g i ðx; t þ DtÞ:



If the desired temperature, Td, which satisfies the thermal boundary condition on the boundary, is given, then from Eq. (21) we can obtain

c2 d X Dt qðx; t þ DtÞ T ¼ g i ðx; t þ DtÞ þ Q b ðx; t þ DtÞ: 2 3T 0 i


Subtracting Eq. (24) from Eq. (25), we can obtain the following direct-forcing formula for the boundary energy-forcing term:

Q b ðx; t þ DtÞ ¼ 2qðx; t þ DtÞ

c2 T d  T noE ðx; t þ DtÞ : 3T 0 Dt


To evaluate Td on the energy-forcing node for complex boundary problems, where the boundary does not match computational nodes, we can directly use the same interface scheme as in the isothermal IB-LBM in case the Dirichlet-type boundary condition is imposed. Although we only consider the Dirichlet-type boundary condition in this paper, if the Neumann boundary is imposed, we can use the procedure of transferring it into the Dirichlet-type boundary condition [7,8] and then apply Eq. (26).

    1 1 RHSn ¼ 4aT ni;j þ a þ uniþ1;j T niþ1;j þ a þ uni1;j T ni1;j 2 2     1 n 1 n n þ a þ v i;jþ1 T i;jþ1 þ a þ v i;j1 T ni;j1 : 2 2


Here, Dx = Dy = 1 as in the standard LBM and T is the normalized temperature by the temperature difference. A von Neumann stability analysis of this discretized equation provides the following stability constraint of Model 2 [37]:

u2 þ v 2 6 a 6 0:25: 2


The direct-forcing IBM can be easily applied to this equation. If there were no energy-forcing, the temperature at the next time step would become

T nþ1 ¼ T ni;j þ RHSn Dt: i;j


Using T noF to discern the temperature from that under external i;j force, we can rewrite Eq. (29) as: n T nþ1 ¼ T noF i;j i;j þ qi;j Dt


From Eqs. (32) and (33), the energy-forcing term can be evaluated as:

qni;j ¼

T di;j  T noE i;j Dt


It should be pointed out that in the actual calculation there is no need to explicitly evaluate the energy-forcing term because the explicit time-advancement scheme is adopted [38]. For the interface scheme, we adopt the same method as in the IB-LBM. 2.4. Newton’s equation of motion for a moving particle

2.3. Immersed boundary method for the finite difference energy equation In the hybrid lattice-Boltzmann finite-difference method (LB FDM), the non-isothermal LBE is solved for velocities and pressures and the finite difference advection–diffusion equation for temperature. For the boundary conditions of the momentum equation – i.e. for no-slip boundary conditions – we adopt the IB-LBM as discussed in Section 2.1, while for the thermal boundary conditions we utilize the IBM for the advection–diffusion equation instead of the simplified thermal LBE in Section 2.2. The energy equation neglecting terms of the viscous dissipation and the compression work done by the pressure can be expressed as:

@T þ r  ðuTÞ ¼ ar2 T þ q @t



Q q¼ qcp

For the simulation of a moving particle, like sedimentation of a cold particle in Section 3.3, we have to consider the motion equations of the particle. The Newton’s equation of translational motion is


dUc ¼ dt


r  dS þ ðqs  qf ÞV s g


where Uc is the center-of-mass velocity of the particle; M, S, V, and q are mass, surface, volume, and density, respectively; and subscripts f and s indicate fluid and solid, respectively. The first term on the right-hand side of Eq. (35) indicates the force from fluid to solid, which consists of stationary surface force and added mass force due to acceleration. The stationary surface force can be expressed in terms of the boundary momentum-forcing based on Eq. (14). Hence,


r  dS ¼ 







Fb dV þ

@ @t

Fb dV þ Mf


qf udV

dUs : dt


where cp is the specific heat. We explicitly discretize Eq. (27) using the first-order forward difference scheme in time and the secondorder central difference scheme in space as in [36]. The resulting discretized equation (Model 2) can be expressed as:

On the other hand, the Newton’s equation of angular motion is

T nþ1 ¼ T ni;j þ ðRHSn þ qni;j ÞDt i;j

where Xc is the angular velocity of the particle, Is is the moment of inertia, and Xw and Xc are position vectors of a wall surface and the



dXc ¼ dt


ðXw  Xc Þ  r  dS




S.K. Kang, Y.A. Hassan / Computers & Fluids 49 (2011) 36–45

where g is the gravitational constant, L is the height or the width of the square cavity, and DT is the temperature difference between hot and cold walls, i.e. DT = Thot  Tcold. Therefore, viscosity and the corresponding non-dimensional relaxation time, and thermal diffusivity and the corresponding non-dimensional thermal relaxation time, can be written as

rffiffiffiffiffiffi Pr U c L; m¼ Ra Uc L ffi; a ¼ pffiffiffiffiffiffiffiffiffiffi RaPr

rffiffiffiffiffiffi Pr 1 Uc L þ ; s¼3 Ra 2 3U c L 1 sg ¼ pffiffiffiffiffiffiffiffiffiffiffi þ ; 2 RaPr 2

ð44aÞ ð44bÞ

where Uc is the characteristic velocity, which is defined by

Uc ¼ Fig. 1. Geometry and boundary conditions of the natural convection in a square cavity.

center as shown in Fig. 1. Eq. (37) can be rewritten in terms of boundary momentum-forcing as c Is dX dt

¼ ¼


ðxb  Xc Þ  Fb dV þ @[email protected] V V

ðxb  Xc Þ  Fb dV þ If


ðXw  Xs Þ  qf udV V

@Xc @t


whereIf ¼ Mf R2s . As a result, the following discretized Newton’s equations of motion corresponding to Eqs. (36) and (38) can be written as

Unþ1 ¼ Unc þ c

" # X n Mf n 1  Fb DV b þ ðMs  Mf Þg Dt þ ðUc  Un1 Þ c Ms M s b

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi bg DTL:


The characteristic velocity should be selected to be small so that the compressibility error remains small. For the athermal and thermal boundary conditions, the following non-equilibrium bounce-back scheme [15] was adopted

fineq ¼ fineq


g neq i



gneq i

where ei ¼ ei . Here, to apply this scheme to the adiabatic boundary, the temperature on the wall, T(y0), should be predetermined. For this, the following second-order accurate discretization was adopted:

@T 1 ¼ ½4Tðy0 þ DyÞ  Tðy0 þ 2DyÞ  3Tðy0 Þ þ OðDy2 Þ ¼ 0: @y y0 2Dy ð47Þ

ð39Þ 3



Xnþ1 ¼ Xnc þ c


X If 1  ðxb  Xc Þ  Fnb DV b Dt þ ðXnc  Xn1 Þ; c Is Is b


respectively. Here, acceleration and angular acceleration terms are discretized based on current and previous time steps, as in Feng and Michaelides [39]. The center position at n + 1 time step can be expressed as:

Xnþ1 c



1 þ ðUnþ1 þ Unc ÞDt: 2 c


Thus, the wall velocity on the boundary point Xw closest to a forcing node can be evaluated as

Unþ1 ¼ Unþ1 þ Xnþ1  ðXw  Xc Þ w s s


and using this velocity and direct-forcing formula, we can obtain the next time step boundary force Fnþ1 b .




Simulations were implemented for Ra = 10 , 10 , 10 , and 10 with the present method, and in all simulations the Prandtl number (Pr) was set to be 0.71. The following convergence criteria [22] were used:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi maxðj ðunþ1 Þ2 þ ðv nþ1 Þ2  ðun Þ2 þ ðv n Þ2 jÞ < 109 ; maxðjT nþ1  T n jÞ < 109


where n and n + 1 indicate the old and the new time steps, respectively. Quantitative comparisons of the maximum horizontal velocity (umax) and its vertical position (ymax) on the vertical centerline (x = L/2), the maximum vertical velocity (vmax) and its horizontal position (xmax) on the horizontal centerline (y = L/2), and the average Nusselt number were performed. The average Nusselt number is defined by

Nuav g ¼

L 1 aDT L2

Z 0




qx ðx; yÞdxdy



3. Results and discussions

with the local heat flux in horizontal direction

3.1. Natural convection in a square cavity

qx ðx; yÞ ¼ uTðx; yÞ  að@[email protected]ÞTðx; yÞ

At first, to test the present LBM code, we considered the natural convection in a square cavity, which has been used as a typical benchmark problem for the validation of the code capability to simulate the natural convection. The square cavity has hot and cold isothermal boundary conditions at left and right vertical sides, respectively, and adiabatic boundary conditions at top and bottom horizontal sides, as shown in Fig. 1. This natural convection is characterized by two non-dimensional numbers: the Rayleigh number (Ra) and the Prandtl number (Pr), which are defined by

where u is x-direction velocity. Table 1 shows the effect of the characteristic velocity on the results from both thermal LBM models under Ra = 104 and the grid size of L/150. It is observed that as the characteristic velocity decreases, the calculation results approach the reference data [40]. Under Uc 6 0.1, the Nusselt number does not show much improvement. Thus, in the remaining calculations, we selected the characteristic velocity as Uc = 0.05. Fig. 2 presents errors of the average Nusselt numbers when changing grid sizes from L/50 to L/200 under Ra = 104 and Uc = 0.05. The error is defined by






m a


Error ¼ Nuav g  NuL=400 av g



S.K. Kang, Y.A. Hassan / Computers & Fluids 49 (2011) 36–45 Table 1 The effect of the characteristic velocity. Model





Model 1

0.01 0.05 0.1 0.2 0.5

2.243 2.243 2.242 2.242 2.215

16.164 16.161 16.157 16.114 15.842

19.610 19.604 19.620 19.565 19.355

Model 2

0.01 0.05 0.1 0.12 0.14

2.243 2.243 2.243 2.242 2.241

16.168 16.165 16.152 16.152 16.147

19.613 19.614 19.596 19.624 19.628






However, in terms of efficiency, the two models showed some differences. Although the simplified thermal LBE model (Model 1) greatly reduced the CPU time compared to He et al.’s original double-population thermal LBE model [22], it still spends about 50% more CPU time than the hybrid thermal LBE model (Model 2) as shown in Table 2. This can be attributed to the use of a larger number of energy distribution functions at collision and streaming steps, and the complicated boundary condition calculation in Model 1. It should be also pointed out that according to Eq. (31), Model 2 has the stability constraint of

u2 þ v 2 Uc L 6 pffiffiffiffiffiffiffiffiffiffiffi 6 0:25: 2 RaPr


For low Rayleigh number conditions, for example, under Ra = 103 and grid size of L/100 and Ra = 104 and grid size of L/150, Eq. (52) provides Uc 6 0.0666 and Uc 6 0.140435, respectively. Therefore, when we selected Uc > 0.0666 under the former condition and Uc > 0.140435 under the latter condition, we could not obtain the converged solution, and thus this range was not included in Table 1. 3.2. Natural convection in a square cavity with an eccentric cylinder

Fig. 2. Accuracy of the two thermal LBM models.

NuL=400 av g

where indicates the average Nusselt number under the grid size of L/400, which is adopted as an exact solution because the analytical solution does not exist. It was observed that the two models have the second-order accuracy in space. Based on such grid sensitivity studies, the grid sizes were selected as L/100, L/150, L/200, and L/250 for Ra = 103, 104, 105, and 106, respectively. Table 2 presents the final results for Ra = 103, 104, 105, and 106. In terms of accuracy, all results showed good agreement between the two models and the reference calculation [40] within 1%.

To check the applicability of the present IBM to complex, fixed boundary, we considered the natural convection of fluid in a square cavity with an eccentric circular cylinder. In a square cavity with a width of L, a circular cylinder with a diameter of D = 0.4L is eccentrically located in the cavity by 0.1L upward from the center, as shown in Fig. 3. Hot and cold temperatures were imposed on the circular cylinder boundary and vertical side walls of the cavity boundary, respectively, and adiabatic conditions are imposed on the horizontal top and bottom walls. This problem has been simulated by various numerical methods under Ra = 106 and Pr = 10 [8,10,41,42]. In the present simulation, the direct-forcing IBMs based on Model 1 and Model 2 were adopted for the boundary conditions on the circular cylinder wall. For the square cavity boundary conditions, the non-equilibrium bounce-back scheme was used as in Section 3.1. The characteristic velocity of 0.2 was adopted. Figs. 4 and 5 present the isotherms and the streamlines form the IBM based on Model 1 and Model 2 under the grid size of L/ 200. Due to the buoyancy force, the heated flows around the hot cylinder move upward and the flows cooled by the cold walls move downward along the cold side walls, thus forming the two

Table 2 Comparison of velocities and Nusselt numbers. Ra = 103

Ra = 104

Ra = 105

Ra = 106


Model 1 Model 2 Reference

3.645 3.646 3.649

16.161 16.165 16.178

34.679 34.680 34.73

64.553 64.596 64.63


Model 1 Model 2 Reference

0.810 0.810 0.813

0.820 0.820 0.823

0.855 0.855 0.855

0.848 0.848 0.850


Model 1 Model 2 Reference

3.694 3.695 3.697

19.604 19.614 19.617

68.527 68.545 68.59

219.670 219.593 219.36


Model 1 Model 2 Reference

0.180 0.180 0.178

0.120 0.120 0.119

0.065 0.065 0.066

0.036 0.036 0.0379


Model 1 Model 2 Reference

1.118 1.118 1.118

2.243 2.243 2.243

4.514 4.514 4.519

8.798 8.794 8.800

CPU time (s)

Model 1 Model 2

396 263

2539 1733

16,182 11,096

47,741 27,151

Fig. 3. Geometry and boundary conditions of the natural convection in a square cavity with an eccentric cylinder.


S.K. Kang, Y.A. Hassan / Computers & Fluids 49 (2011) 36–45

Fig. 4. Isotherms obtained from the IBM based on (a) Model 1 and (b) Model 2.

Fig. 5. Streamlines obtained from the IBM based on (a) Model 1 and (b) Model 2.

symmetric free circulations. The isotherms and streamlines show very good agreement with those from Ref. [8]. To quantitatively validate the present method, we compared the local Nusselt numbers along one of the cold side walls with different grid sizes of L/100 and L/200 with those from other numerical methods [41,42] as shown in Fig. 6. It was observed that the results from the IBM based on both models – even under the grid size of L/100 – showed good agreement with those under high resolution although there were some discrepancies near the top wall, which are removed under the grid size of L/200. In terms of numerical efficiency, the IBM based on Model 1 spent 50% more time than the IBM based on Model 2. 3.3. Cold particle sedimentation As mentioned in the Introduction, one of great advantages of the IBM is that it can be easily applied to the moving-boundary problem without re-meshing. To validate the applicability of the present IBM to the heat transfer problem involving moving, complex boundaries, we considered a cold-particle sedimentation in a hot, infinite channel. This problem is challenging because it involves complicated mechanisms between the forced and natural thermal convections and strong wall confinement effect [42]. Firstly, Gan et al. [43] used the Arbitrary Lagrangian–Eulerian (ALE) method to simulate this problem where the particle had been

initially released at the centerline. They suggested the five flow regimes according to different Grashof numbers based on the particle equilibrium positions and the wake structures. Subsequently, Yu et al. [42] used the fictitious domain method to simulate the problem. However, they initially located the particle off the centerline by one particle radius, because the trajectories in some regimes where the particle migrates away from the centerline are not deterministic since the migrations depend on the random numerical disturbances. They also found the different regimes from Gan et al.’s at high Grashof numbers over Gr > 4000. Especially at Gr = 4500, it was observed that flows become turbulent-like and the particle oscillates violently but still regularly. They attributed the difference to the use of fine meshes only for the region in the vicinity of the particle boundary, resulting in a lack of description of far field in the ALE method. This modified Grashof regime at high Grashof numbers was confirmed by the IBM based on the NSE [10]. In the present simulation, we followed the conditions in [42] so that the direct comparison is possible. As shown in Fig. 7, a circular cylinder with a diameter of D was initially located off the centerline by D/2. The Prandtl number was 0.7 and the density ratio of the solid particle to the fluid (qr) was 1.00232. The cylinder had the constant cold temperature and the side walls had the constant hot temperature. Fluid is initially at rest with the hot temperature. The characteristic velocity as adopted in previous sections is no longer valid in this problem because the flow involves forced

S.K. Kang, Y.A. Hassan / Computers & Fluids 49 (2011) 36–45

Fig. 6. Local Nusselt number variation along the cold wall for L = 200Dx from the IBM based on (a) Model 1 and (b) Model 2.


Fig. 8. Horizontal position evolutions obtained from the IBM based on (a) Model 1 and (b) Model 2.

convection due to particle movement. The relaxation time was selected as 0.65, as in [28]. To determine the gravitation constant in the LBM frame, we used the reference Reynolds number of 40.5, as in [10,42]. Here, the reference Reynolds number is defined by [42]

Reref ¼

U ref D



where the reference velocity Uref is defined by

U ref ¼

Fig. 7. Geometry and boundary conditions in a cold particle sedimentation problem.

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pðD=2Þðqr  1Þg :


To simulate the infinite channel, the computational domain of 4D  160D was adopted. The present IBMs based on Model 1 and Model 2 were applied to the cold wall particle boundary. The motion of the particle is calculated using Eqs. (39) and (46). For isothermal case (i.e. Gr = 0), the terminal Reynolds numbers (ReT = uTD/m) with the terminal velocity (uT) obtained from the IBLBM based on both models was 21.2, which is the same as that from Ref. [42]. Fig. 8 presents the time evolution of the particle horizontal positions under different Grashof numbers and Fig. 9 isotherms and vorticity contours at time t = 129.6, where the time t is defined as t = D/Uref. The results follow the modified flow regimes


S.K. Kang, Y.A. Hassan / Computers & Fluids 49 (2011) 36–45

Fig. 9. Isotherms (left) and vorticity contours (right) at time t = 129.6 under (a) Gr = 100, (b) Gr = 564, (c) Gr = 1000, (d) Gr = 2000, (e) Gr = 2500, and (f) Gr = 4500.

Table 3 Equilibrium positions at Gr = 1000 and 2000 and amplitudes at Gr = 4500. Gr

Present (Model 1)

Present (Model 2)

Yu et al. [42]

Feng and Michaelides [39]

1000 2000 4500

2.89 2.74 1.32

2.91 2.74 1.32

2.89 2.74 1.33

2.90 2.73 1.35

well. At Gr = 100 (Regime A in [43]), the particle settles steadily along the centerline and the wake vortices are steady and symmetric (Figs. 8 and 9a). At Gr = 564 (Regime B in [43]), vortex shedding occurs from the particle and the particle oscillates regularly about the centerline (Figs. 8 and 9b). At Gr = 1000 and 2000 (Regime C in [43]), two types of migrations are observed: one with oscillation as a nature extension of Regime B (Figs. 8 and 9c), and the other without oscillation (Figs. 8 and 9d). At Gr = 2500 (Regime D in [43]), the centerline becomes again a stable equilibrium position and vortex shedding is absent (Figs. 8 and 9e). At Gr = 4500 (Regime F in [42]), the flows become turbulent-like and the particle oscillates violently but still regularly (Figs. 8 and 9f). For the quantitative comparison, we compared the particle horizontal positions between the present results and the previous results. Table 3 presents a comparison of equilibrium positions at Gr = 1000 and 2000, and the amplitude at Gr = 4500 between the present scheme, the fictitious domain method of [42] and the IBM based on Naiver–Stokes equation of [10]. Results from both models show good agreement with others. From this simulation, it is confirmed that the present directforcing IBM based on the two thermal LBM models can produce accurate results even for the heat transfer problem with moving boundaries. Again, the IBM based on Model 2 is numerically 50% more efficient than that based on Model 1.

4. Conclusions In this study, we proposed the direct-forcing IBMs coupled with two thermal LBM models: the double-population model with simplified internal energy distribution functions (Model 1), and the hybrid model with a finite difference advection–diffusion equation (Model 2). We showed that the IBMs based on both models had good accuracies for heat transfer problems with not only stationary but also moving complex boundaries where the Boussinesq approximation is valid, and that the viscous heat dissipation and the compression work done by the pressure was negligible. However, the IBM based on Model 2 was faster than that based on Model 1 by 50%. Therefore, the IBM based on the hybrid thermal LBM model is recommended for the actual calculation in terms of efficiency as well as accuracy. In this study, we only considered the SRT lattice Boltzmann equation for velocity and pressure. However, the multi-relaxation-time (MRT) lattice Boltzmann equation could be adopted for the better stability. In addition, the MRT-model could also consider temperature-dependent transport coefficients without explicitly using the Boussinesq approximation [31]. We predict that the direct-forcing formula can be coupled with other simplified double-population TLBE models [23,24]. Also, applying the present method to an improved single-population model, such as [21], is worth considering as future work. In addition, it is expected that we can easily extend the present method to 3D problems by adopting the D3Q19 thermal LBE and the interface scheme based on tri-linear, bi-linear, and linear interpolations. References [1] Fadlun EA, Verzicco R, Orlandi P, Mohd-Yusof J. Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations. J Comput Phys 2000;161(1):35–60.

S.K. Kang, Y.A. Hassan / Computers & Fluids 49 (2011) 36–45 [2] Kim J, Kim D, Choi H. An immersed-boundary finite-volume method for simulations of flow in complex geometries. J Comput Phys 2001;171(1):132–50. [3] Goldstein D, Handler R, Sirovich L. Modeling a no-slip flow boundary with an external force-field. J Comput Phys 1993;105(2):354–66. [4] Saiki EM, Biringen S. Numerical simulation of a cylinder in uniform flow: application of a virtual boundary method. J Comput Phys 1996;123(2):450–65. [5] Iaccarino G, Verzicco R. Immersed boundary technique for turbulent flow simulations. Appl Mech Rev 2003;56:331–47. [6] Balaras E. Modeling complex boundaries using an external force field on fixed Cartesian grids in large-eddy simulations. Comput Fluids 2004;33(3):375–404. [7] Kim J, Choi HC. An immersed-boundary finite-volume method for simulation of heat transfer in complex geometries. KSME Int J 2004;18(6):1026–35. [8] Pacheco JR, Pacheco-Vega A, Rodic T, Peck RE. Numerical simulations of heat transfer and fluid flow problems using an immersed-boundary finite-volume method on nonstaggered grids. Numer Heat Transfer B – Fund 2005;48(1):1–24. [9] Zhang N, Zheng ZC, Eckels S. Study of heat-transfer on the surface of a circular cylinder in flow using an immersed-boundary method. Int J Heat Fluid Flow 2008;29(6):1558–66. [10] Feng ZG, Michaelides EE. Heat transfer in particulate flows with direct numerical simulation (DNS). Int J Heat Mass Transfer 2009;52:777–86. [11] Chen S, Doolen GD. Lattice Boltzmann method for fluid flows. Annu Rev Fluid Mech 1998;30:329–64. [12] Yu DZ, Mei RW, Luo LS, Shyy W. Viscous flow computations with the method of lattice Boltzmann equation. Prog Aerosp Sci 2003;39(5):329–67. [13] Kang SK, Hassan YA. A comparative study of direct-forcing immersed boundary-lattice Boltzmann methods for stationary complex boundaries. Int J Numer Methods Fluids 2010. doi:10.1002/fld.2304. [14] Jeong HK, Yoon HS, Ha MY, Tsutahara M. An immersed boundary-thermal lattice Boltzmann method using an equilibrium internal energy density approach for the simulation of flows with heat transfer. J Comput Phys 2010;229(7):2526–43. [15] He X, Chen S, Doolen GD. A novel thermal model for the lattice Boltzmann method in incompressible limit. J Comput Phys 1998;146(1):282–300. [16] Alexander FJ, Chen S, Sterling JD. Lattice Boltzmann thermohydrodynamics. Phys Rev E 1993;47(4):R2249–52. [17] Mcnamara G, Alder B. Analysis of the lattice Boltzmann treatment of hydrodynamics. Physica A 1993;194(1–4):218–28. [18] Mcnamara GR, Garcia AL, Alder BJ. Stabilization of thermal lattice Boltzmann models. J Stat Phys 1995;81(1–2):395–408. [19] Shan XW, Yuan XF, Chen HD. Kinetic theory representation of hydrodynamics: a way beyond the Navier–Stokes equation. J Fluid Mech 2006;550:413–41. [20] Philippi PC, Hegele LA, dos Santos LOE, Surmas R. From the continuous to the lattice Boltzmann equation: the discretization problem and thermal models. Phys Rev E 2006;73(5):056702. [21] Prasianakis NI, Karlin IV. Lattice Boltzmann method for thermal flow simulation on standard lattices. Phys Rev E 2007;76(1):016702. [22] Peng Y, Shu C, Chew YT. Simplified thermal lattice Boltzmann model for incompressible thermal flows. Phys Rev E 2003;68(2):026701. [23] Guo ZL, Shi BC, Zheng CG. A coupled lattice BGK model for the Boussinesq equations. Int J Numer Methods Fluids 2002;39(4):325–42.


[24] Li Q, He YL, Wang Y, Tang GH. An improved thermal lattice Boltzmann model for flows without viscous heat dissipation and compression work. Int J Mod Phys C 2008;19(1):125–50. [25] Yuan P, Schaefer L. A thermal lattice Boltzmann two-phase flow model and its application to heat transfer problems – Part 2. Integration and validation. J Fluid Eng – Trans ASME 2006;128(1):151–6. [26] Wang JK, Wang M, Li ZX. A lattice Boltzmann algorithm for fluid-solid conjugate heat transfer. Int J Therm Sci 2007;46(3):228–34. [27] Zu YQ, Yan YY, Shi WP, Ren LQ. Numerical method of lattice Boltzmann simulation for flow past a rotating circular cylinder with heat transfer. Int J Numer Method H 2008;18(5–6):766–82. [28] Mandujano F, Rechtman R. Thermal levitation. J Fluid Mech 2008;606:105–14. [29] Barrios G, Rechtman R, Rojas J, Tovar R. The lattice Boltzmann equation for natural convection in a two-dimensional cavity with a partially heated wall. J Fluid Mech 2005;522:91–100. [30] Xuan YM, Yu K, Li Q. Investigation on flow and heat transfer of nanofluids by the thermal lattice Boltzmann model. Prog Comput Fluid Dynam 2005;5(1– 2):13–9. [31] Lallemand P, Luo LS. Theory of the lattice Boltzmann method: acoustic and thermal properties in two and three dimensions. Phys Rev E 2003;68(3):036706. [32] Filippova O, Hanel D. Lattice-BGK model for low Mach number combustion. Int J Mod Phys C 1998;9(8):1439–45. [33] Filippova O, Hanel D. A novel lattice BGK approach for low Mach number combustion. J Comput Phys 2000;158(2):139–60. [34] Mezrhab A, Bouzidi W, Lallemand P. Hybrid lattice-Boltzmann finitedifference simulation of convective flows. Comput Fluids 2004;33(4):623–41. [35] Guo ZL, Zheng CG, Shi BC. Discrete lattice effects on the forcing term in the lattice Boltzmann method. Phys Rev E 2002;65(4):046308. [36] Moussaoui MA, Jami M, Mezrhab A, Naji H. MRT-lattice Boltzmann simulation of forced convection in a plane channel with an inclined square cylinder. Int J Therm Sci 2010;49(1):131–42. [37] Hindmarsh AC, Gresho PM, Griffiths DF. The stability of explicit Euler timeintegration for certain finite-difference approximations of the multidimensional advection diffusion equation. Int J Numer Methods Fluids 1984;4(9):853–97. [38] Yang JM, Balaras E. An embedded-boundary formulation for large-eddy simulation of turbulent flows interacting with moving boundaries. J Comput Phys 2006;215(1):12–40. [39] Feng ZG, Michaelides EE. Robust treatment of no-slip boundary condition and velocity updating for the lattice-Boltzmann simulation of particulate flows. Comput Fluids 2008;38(2):370–81. [40] Davis GD. Natural-convection of air in a square cavity – a bench-mark numerical-solution. Int J Numer Methods Fluids 1983;3(3):249–64. [41] Demirdzic I, Lilek Z, Peric M. Fluid-flow and heat-transfer test problems for nonorthogonal grids – bench-mark solutions. Int J Numer Methods Fluids 1992;15(3):329–54. [42] Yu ZS, Shao XM, Wachs A. A fictitious domain method for particulate flows with heat transfer. J Comput Phys 2006;217(2):424–52. [43] Gan H, Chang JZ, Feng JJ, Hu HH. Direct numerical simulation of the sedimentation of solid particles with thermal convection. J Fluid Mech 2003;481:385–411.