# Simulation of Thermal Flow Problems via a Hybrid Immersed Boundary-Lattice Boltzmann Method

Authors: J. Wu; C. Shu; N. Zhao
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

A hybrid immersed boundary-lattice Boltzmann method (IB-LBM) is presented in this work to simulate the thermal flow problems. In current approach, the flow field is resolved by using our recently developed boundary condition-enforced IB-LBM (Wu and Shu, (2009)). The nonslip boundary condition on the solid boundary is enforced in simulation. At the same time, to capture the temperature development, the conventional energy equation is resolved. To model the effect of immersed boundary on temperature field, the heat source term is introduced. Different from previous studies, the heat source term is set as unknown rather than predetermined. Inspired by the idea in (Wu and Shu, (2009)), the unknown is calculated in such a way that the temperature at the boundary interpolated from the corrected temperature field accurately satisfies the thermal boundary condition. In addition, based on the resolved temperature correction, an efficient way to compute the local and average Nusselt numbers is also proposed in this work. As compared with traditional implementation, no approximation for temperature gradients is required. To validate the present method, the numerical simulations of forced convection are carried out. The obtained results show good agreement with data in the literature.


---

## Body

## 1. Introduction

As a long-standing challenge in the computational fluid mechanics, the flow problems with complex geometry have been widely studied. In terms of grid applied, the numerical methods can be generally classified into two categories. In the traditional approach, the body-fitted mesh is employed to discretize the governing equation [1]. Thereafter, the boundary condition could be implemented directly. As a result, the solution of problems strongly depends on the quality of generated mesh. Moreover, the mesh regeneration procedure is usually unavoidable when the object is not fixed. To relieve this cumbersome requirement, the alternative choice is to decouple the governing equation from the computational mesh. In this category, the most famous algorithm is the immersed boundary method (IBM) proposed by Peskin [2].In IBM, the discretization of governing equation is performed on the Cartesian (Eulerian) mesh, and the boundary of object is represented through a set of Lagrangian points. Different from the body-fitted solver, the boundary condition in this method is depicted by introducing a body force (restoring force) into the governing equation, which replaces the effect of boundary on the surrounding flow field. Since the governing equation is independent of the boundary, IBM is extremely suitable for handling various flow problems with complex geometry.To acquire accurate solution via IBM, the appropriate calculation of force density is of great importance. Up till now, there are several ways to fulfill this task. By treating the boundary as deformable with high stiffness, Peskin [2] firstly applied the Hook’s law, associated with the use of free parameter, to compute the restoring force (named penalty method). Fadlun et al. [3] viewed the momentum equations as to be satisfied at boundary points and hence the force density could be subsequently computed (named direct forcing method). This approach has been widely used in current IBM application. Different from two methods introduced above, which are based on the Navier-Stokes (N-S) solver, another force calculation technique was proposed by Niu et al. [4] in the framework of lattice Boltzmann method (LBM) [5]. The concept of momentum exchange at the boundary is utilized in this approach (named momentum exchange method). Nonetheless, it should be stressed that the force density is calculated explicitly in the conventional methods. Consequently, the nonslip boundary condition cannot be guaranteed. Recently, we presented a boundary condition-enforced IB-LBM [6], in which the force density is raised from the velocity correction and is computed implicitly by enforcing the boundary condition. As compared to the conventional IBM, the nonslip boundary condition can be guaranteed in this method.The idea of IBM has been extensively introduced into numerous isothermal flow problems. On the other hand, its application in thermal flow problems is limited and is still in progress. A few attempts have been made in this aspect [7–9]. Similar to the use of force density in isothermal flow problems, a heat source term is employed to meet the influence of boundary in thermal flow problems. However, same as the force calculation in conventional methods, the heat source term is computed explicitly in recent studies [7–9]. In this way, the thermal boundary condition cannot be accurately satisfied, which may affect the accuracy of solution. To address this problem, in this work, we develop a hybrid IB-LBM for simulation of thermal flow problems following the idea of velocity correction in [6]. In this method, the flow field is solved based on the IB-LBM in [6]. In the meanwhile, the temperature field is obtained by solving the conventional energy equation with additional heat source term, which is equivalent to temperature correction. With the temperature correction evaluated implicitly, the thermal boundary condition can be satisfied strictly. In addition, the Nusselt number is a crucial parameter in thermal flow problems. Based on the established temperature correction, in this paper, both the local and average Nusselt numbers can be efficiently computed. In this manner, the troublesome temperature gradient calculation at the boundary points could be eliminated. To validate the proposed algorithm, the forced convection from a stationary isothermal circular cylinder is simulated. The numerical results show good agreement with available data in the literature.


## 2. Methodology

### 2.1. Boundary Condition-Enforced Immersed Boundary-Lattice Boltzmann Method

For viscous incompressible flow problems with an immersed boundary, the governing equations in the framework of lattice Boltzmann equation (LBE) can be written as(2.1)fα(x+eαδt,t+δt)-fα(x,t)=-1τ(fα(x,t)-fαeq(x,t))+Fαδt,(2.2)Fα=(1-12τ)wα(eα-ucs2+eα⋅ucs4eα)⋅f,(2.3)ρu=∑αeαfα+12fδt,
where fα is the distribution function, fαeq is its corresponding equilibrium state, τ is the single relaxation parameter, eα is the lattice velocity, wα are coefficients which depend on the selected lattice velocity model, and f is the force density which is distributed from the boundary force density. Here, the D2Q9 lattice velocity model is used, which is(2.4)eα={(0,0),α=0,(±1,0),(0,±1),α=1~4,(±1,±1),α=5~8,
and the corresponding equilibrium distribution function is(2.5)fαeq(x,t)=ρwα[1+eα⋅ucs2+(eα⋅u)2-(cs|u|)22cs4]
with w0=4/9, wα=1/9, (α=1~4), and wα=1/36, (α=5~8). cs=1/3 is the sound speed of this model. The relationship between the relaxation time and the kinematic viscosity of fluid is υ=(τ-0.5)cs2δt.In order to satisfy the nonslip boundary condition, the force densityf in (2.2) and (2.3) should be set as unknown [6]. It can be calculated from the fluid velocity correction δu. Furthermore, δu can be obtained from the boundary velocity correction δuB. The final expression for δuB is(2.6)AX=B,
where
(2.7)X={δuB1,δuB2,…,δuBm}T,(2.8)A=(δ11δ12⋯δ1nδ21δ22⋯δ2n⋮⋮⋱⋮δm1δm2⋯δmn)(δ11Bδ12B⋯δ1mBδ21Bδ22B⋯δ2mB⋮⋮⋱⋮δn1Bδn2B⋯δnmB),B=(UB1UB2⋮UBm)-(δ11δ12⋯δ1nδ21δ22⋯δ2n⋮⋮⋱⋮δm1δm2⋯δmn)(u1*u2*⋮un*).
Here, m is the number of Lagrangian (boundary) points, and n is the number of surrounding Eulerian points. δuBl(l=1,2,…,m) is the unknown velocity correction vector at the boundary points. δij=Dij(xij-XBl)ΔxΔy and δijB=Dij(xij-XBl)Δsl. Dij(xij-XBl) is the delta function, which is expressed as(2.9)Dij(xij-XBl)=δ(xij-XBl)δ(yij-YBl),δ(r)={14h(1+cos⁡(π|r|2h)),|r|≤2h,0,|r|>2h,
and Δsl is the arc length of boundary element. UBl is the boundary velocity. The intermediate fluid velocity u* is calculated by(2.10)u*=1ρ∑αeαfα.
By solving equation system (2.6), the unknown variables δuBl can be obtained. After that, the fluid velocity correction δu can be calculated by(2.11)δu=∑lδuBlDij(xij-XBl)Δsl.
As shown in (2.3), the relationship between the force density f and the fluid velocity correction δu can be written as(2.12)f=2ρδuδt.
In the LBM, the macroscopic variables such as density and pressure are calculated by(2.13)ρ=∑αfα,P=cs2ρ.


### 2.2. Implicit Temperature Correction for Energy Equation

After the flow field is accurately solved through the boundary condition-enforced IB-LBM [6], the temperature field can be obtained by using the conventional energy equation with additional heat source term, which can be written as(2.14)ρcp(∂T∂t+(u⋅∇)T)=k∇2T+q,
where T is the temperature, cp is the specific heat of the fluid, k is the thermal conductivity of the fluid, and q is the heat source density distributed from the boundary heat flux.To solve (2.14), the fractional step technique is used. First, the energy equation without heat source density is solved, which is called Predictor step:(2.15)ρcp(∂T∂t+(u⋅∇)T)=k∇2T.
The solution of (2.15) is regarded as intermediate temperature T*. In the second step, or named Corrector step, the following equation is solved:(2.16)ρcp∂T∂t=q.
Here, (2.16) means the effect of heat source on temperature field. Using the Euler explicit scheme, ∂T/∂t can be rewritten as Tn+1-T*/δt=δT/δt, where Tn+1 is the required temperature at next time step. Therefore, (2.16) becomes(2.17)q=ρcpδTδt,
where δT is the temperature correction which determines the unknown variable q. It means that to calculate heat source term is equivalent to correcting the temperature field near the boundary. Similar to the implementation of velocity correction in IB-LBM [6], δT can further be calculated from the boundary temperature correction δTB. By enforcing the constant temperature boundary condition, δTB can computed from the following equation system:(2.18)CY=D,
where
(2.19)Y={δTB1,δTB2,…,δTBm}T,(2.20)C=(δ11δ12⋯δ1nδ21δ22⋯δ2n⋮⋮⋱⋮δm1δm2⋯δmn)(δ11Bδ12B⋯δ1mBδ21Bδ22B⋯δ2mB⋮⋮⋱⋮δn1Bδn2B⋯δnmB),D=(TB1TB2⋮TBm)-(δ11δ12⋯δ1nδ21δ22⋯δ2n⋮⋮⋱⋮δm1δm2⋯δmn)(T1*T2*⋮Tn*).
Here, TBl is the specified boundary temperature. The parameters m, n, δij, and δijB have the same meaning as those in equation system (2.6). After δTB is obtained from equation system (2.18), the temperature correction δT can be calculated through(2.21)δT=∑lδTBlDij(xij-XBl)Δsl.
Consequently, the corrected temperature can be computed by(2.22)T=T*+δT.


### 2.3. Local and Average Nusselt Number Evaluation

In thermal flow problems, the Nusselt number is an important measurement parameter. The local Nusselt numberNu on the surface of boundary is defined as(2.23)Nu(XB,t)=h(XB,t)Lck,
where h is the local convective heat transfer coefficient and Lc is the characteristic length. According to Newton’s cooling law and Fourier’s law, the heat conducted away from the boundary should be equal to the heat convection from the boundary, that is,(2.24)-k∂T∂n(XB,t)=h(XB,t)(TB-T∞),
where T∞ is the free stream temperature. Substituting (2.24) into (2.23) gives(2.25)Nu(XB,t)=-LcTB-T∞∂T∂n(XB,t).
Further averaging over the entire boundary, we can have the surface overall mean Nusselt number Nu¯(2.26)Nu¯=1L∫Nuds=-Lc(TB-T∞)L∫∂T∂n(XB,t)ds,
where L is the total length of boundary. Usually, the average Nusselt number can be used to estimate the rate of heat transfer from the heated surface.As shown in (2.25) and (2.26), to evaluate the local and average Nusselt numbers, the calculation of temperature gradient at the boundary point is required. Since the temperature is located at the Eulerian mesh which is not coincided with boundary, we cannot calculate the boundary temperature derivatives directly. On the other hand, according to the Fourier’s law, we can have the following expression for heat flux on the boundary Q:(2.27)Q(XB,t)=-k∂T∂n(XB,t).
At the same time, similar to (2.17), the boundary heat flux can be calculated using the boundary temperature correction(2.28)Q=ρcpδTBδt.
Substituting (2.27) and (2.28) into (2.25) and (2.26), the expressions for the local and average Nusselt numbers can be rewritten as(2.29)Nu=ρcpLckδt(TB-T∞)δTB,(2.30)Nu¯=ρcpLckδt(TB-T∞)L∑lδTBlΔsl.
As can be seen from (2.29) and (2.30), the local and average Nusselt numbers can be easily calculated by using the solved boundary temperature correction, which avoids the temperature gradient computation on the boundary.


## 2.1. Boundary Condition-Enforced Immersed Boundary-Lattice Boltzmann Method

For viscous incompressible flow problems with an immersed boundary, the governing equations in the framework of lattice Boltzmann equation (LBE) can be written as(2.1)fα(x+eαδt,t+δt)-fα(x,t)=-1τ(fα(x,t)-fαeq(x,t))+Fαδt,(2.2)Fα=(1-12τ)wα(eα-ucs2+eα⋅ucs4eα)⋅f,(2.3)ρu=∑αeαfα+12fδt,
where fα is the distribution function, fαeq is its corresponding equilibrium state, τ is the single relaxation parameter, eα is the lattice velocity, wα are coefficients which depend on the selected lattice velocity model, and f is the force density which is distributed from the boundary force density. Here, the D2Q9 lattice velocity model is used, which is(2.4)eα={(0,0),α=0,(±1,0),(0,±1),α=1~4,(±1,±1),α=5~8,
and the corresponding equilibrium distribution function is(2.5)fαeq(x,t)=ρwα[1+eα⋅ucs2+(eα⋅u)2-(cs|u|)22cs4]
with w0=4/9, wα=1/9, (α=1~4), and wα=1/36, (α=5~8). cs=1/3 is the sound speed of this model. The relationship between the relaxation time and the kinematic viscosity of fluid is υ=(τ-0.5)cs2δt.In order to satisfy the nonslip boundary condition, the force densityf in (2.2) and (2.3) should be set as unknown [6]. It can be calculated from the fluid velocity correction δu. Furthermore, δu can be obtained from the boundary velocity correction δuB. The final expression for δuB is(2.6)AX=B,
where
(2.7)X={δuB1,δuB2,…,δuBm}T,(2.8)A=(δ11δ12⋯δ1nδ21δ22⋯δ2n⋮⋮⋱⋮δm1δm2⋯δmn)(δ11Bδ12B⋯δ1mBδ21Bδ22B⋯δ2mB⋮⋮⋱⋮δn1Bδn2B⋯δnmB),B=(UB1UB2⋮UBm)-(δ11δ12⋯δ1nδ21δ22⋯δ2n⋮⋮⋱⋮δm1δm2⋯δmn)(u1*u2*⋮un*).
Here, m is the number of Lagrangian (boundary) points, and n is the number of surrounding Eulerian points. δuBl(l=1,2,…,m) is the unknown velocity correction vector at the boundary points. δij=Dij(xij-XBl)ΔxΔy and δijB=Dij(xij-XBl)Δsl. Dij(xij-XBl) is the delta function, which is expressed as(2.9)Dij(xij-XBl)=δ(xij-XBl)δ(yij-YBl),δ(r)={14h(1+cos⁡(π|r|2h)),|r|≤2h,0,|r|>2h,
and Δsl is the arc length of boundary element. UBl is the boundary velocity. The intermediate fluid velocity u* is calculated by(2.10)u*=1ρ∑αeαfα.
By solving equation system (2.6), the unknown variables δuBl can be obtained. After that, the fluid velocity correction δu can be calculated by(2.11)δu=∑lδuBlDij(xij-XBl)Δsl.
As shown in (2.3), the relationship between the force density f and the fluid velocity correction δu can be written as(2.12)f=2ρδuδt.
In the LBM, the macroscopic variables such as density and pressure are calculated by(2.13)ρ=∑αfα,P=cs2ρ.


## 2.2. Implicit Temperature Correction for Energy Equation

After the flow field is accurately solved through the boundary condition-enforced IB-LBM [6], the temperature field can be obtained by using the conventional energy equation with additional heat source term, which can be written as(2.14)ρcp(∂T∂t+(u⋅∇)T)=k∇2T+q,
where T is the temperature, cp is the specific heat of the fluid, k is the thermal conductivity of the fluid, and q is the heat source density distributed from the boundary heat flux.To solve (2.14), the fractional step technique is used. First, the energy equation without heat source density is solved, which is called Predictor step:(2.15)ρcp(∂T∂t+(u⋅∇)T)=k∇2T.
The solution of (2.15) is regarded as intermediate temperature T*. In the second step, or named Corrector step, the following equation is solved:(2.16)ρcp∂T∂t=q.
Here, (2.16) means the effect of heat source on temperature field. Using the Euler explicit scheme, ∂T/∂t can be rewritten as Tn+1-T*/δt=δT/δt, where Tn+1 is the required temperature at next time step. Therefore, (2.16) becomes(2.17)q=ρcpδTδt,
where δT is the temperature correction which determines the unknown variable q. It means that to calculate heat source term is equivalent to correcting the temperature field near the boundary. Similar to the implementation of velocity correction in IB-LBM [6], δT can further be calculated from the boundary temperature correction δTB. By enforcing the constant temperature boundary condition, δTB can computed from the following equation system:(2.18)CY=D,
where
(2.19)Y={δTB1,δTB2,…,δTBm}T,(2.20)C=(δ11δ12⋯δ1nδ21δ22⋯δ2n⋮⋮⋱⋮δm1δm2⋯δmn)(δ11Bδ12B⋯δ1mBδ21Bδ22B⋯δ2mB⋮⋮⋱⋮δn1Bδn2B⋯δnmB),D=(TB1TB2⋮TBm)-(δ11δ12⋯δ1nδ21δ22⋯δ2n⋮⋮⋱⋮δm1δm2⋯δmn)(T1*T2*⋮Tn*).
Here, TBl is the specified boundary temperature. The parameters m, n, δij, and δijB have the same meaning as those in equation system (2.6). After δTB is obtained from equation system (2.18), the temperature correction δT can be calculated through(2.21)δT=∑lδTBlDij(xij-XBl)Δsl.
Consequently, the corrected temperature can be computed by(2.22)T=T*+δT.


## 2.3. Local and Average Nusselt Number Evaluation

In thermal flow problems, the Nusselt number is an important measurement parameter. The local Nusselt numberNu on the surface of boundary is defined as(2.23)Nu(XB,t)=h(XB,t)Lck,
where h is the local convective heat transfer coefficient and Lc is the characteristic length. According to Newton’s cooling law and Fourier’s law, the heat conducted away from the boundary should be equal to the heat convection from the boundary, that is,(2.24)-k∂T∂n(XB,t)=h(XB,t)(TB-T∞),
where T∞ is the free stream temperature. Substituting (2.24) into (2.23) gives(2.25)Nu(XB,t)=-LcTB-T∞∂T∂n(XB,t).
Further averaging over the entire boundary, we can have the surface overall mean Nusselt number Nu¯(2.26)Nu¯=1L∫Nuds=-Lc(TB-T∞)L∫∂T∂n(XB,t)ds,
where L is the total length of boundary. Usually, the average Nusselt number can be used to estimate the rate of heat transfer from the heated surface.As shown in (2.25) and (2.26), to evaluate the local and average Nusselt numbers, the calculation of temperature gradient at the boundary point is required. Since the temperature is located at the Eulerian mesh which is not coincided with boundary, we cannot calculate the boundary temperature derivatives directly. On the other hand, according to the Fourier’s law, we can have the following expression for heat flux on the boundary Q:(2.27)Q(XB,t)=-k∂T∂n(XB,t).
At the same time, similar to (2.17), the boundary heat flux can be calculated using the boundary temperature correction(2.28)Q=ρcpδTBδt.
Substituting (2.27) and (2.28) into (2.25) and (2.26), the expressions for the local and average Nusselt numbers can be rewritten as(2.29)Nu=ρcpLckδt(TB-T∞)δTB,(2.30)Nu¯=ρcpLckδt(TB-T∞)L∑lδTBlΔsl.
As can be seen from (2.29) and (2.30), the local and average Nusselt numbers can be easily calculated by using the solved boundary temperature correction, which avoids the temperature gradient computation on the boundary.


## 3. Results and Discussion

As a benchmark case, the forced convection from a stationary circular cylinder has been extensively simulated to validate the numerical methods. Both the experimental and numerical results are available in the literature [10–12]. To characterize this thermal flow problem, two nondimensional parameters are used: the Reynolds number Re=ρU∞D/μ and Prandtl number Pr=μcp/k. Here, U∞ is the free stream velocity, D is the cylinder diameter, and μ is the dynamic viscosity. In current simulation, the Prandtl number is fixed at Pr=0.7. Meanwhile, several low Reynolds numbers are selected: Re=20,40 (steady flow), and 80, 150 (unsteady flow). Since this is a one-way interaction problem, the velocity field can influence the temperature field while cannot be affected by the temperature field. As the developed IB-LBM can accurately simulate velocity field [6], we only focus on the temperature field in this case.For the steady flow simulation, the computational domain is60D×40D with the nonuniform mesh size of 321×241. The cylinder is located at (20D,20D). The region around the cylinder is 1.2D×1.2D with a finest and uniform mesh size of 49×49. To apply the IB-LBM on nonuniform mesh, the Taylor series expansion and least-squares-based LBM [13] are utilized. Figure 1 shows the temperature contours at Re=20 and 40. It can be seen from the figure that the isotherms cluster heavily in the front part of cylinder. When the Reynolds number is increased, the behavior of isotherm clustering is strengthened. The clustering of isotherm indicates that the temperature gradient is large there. As a result, the heat transfer rate near the front of cylinder is much larger than other regions. This phenomenon can be further verified via the local Nusselt number distribution on the surface of cylinder, as shown in Figure 2 for the case of Re=20. Without evaluation of temperature gradient on the boundary point, the local Nusselt number can be easily computed from the solved boundary temperature correction, as shown in (2.29). To make comparison, the result of Bharti et al. [12] is also involved. In this figure, θ=0∘ means the front stagnation point of cylinder. From Figure 2, it can be found that the value of local Nusselt number located at θ=0∘ is maximum and it decreases monotonously with respect to θ. Also can be seen from the figure is that the result of current simulation shows good agreement with that of Bharti et al. [12].Figure 1
Temperature contours for steady forced convection from a stationary circular cylinder.Figure 2
Local Nusselt number distribution on the surface of cylinder atRe=20.For the unsteady flow simulation, the computational domain is50D×40D with the mesh size of 401×301. The region around the cylinder keeps 1.2D×1.2D with a uniform mesh size of 97×97. The cylinder is still located at (20D,20D). Figure 3 plots the isotherms as well as vorticity contours at Re=80 and 150. The regular vortex shedding occurs at these two Reynolds numbers. Concurrently, the isotherms lose symmetry and start to display shedding behavior, which synchronously moves downstream with vortex.Instantaneous isotherms (colorful) and vorticity contours (gray) for unsteady forced convection from a stationary cylinder.
(a)(b)To illustrate the performance of average Nusselt number calculation by using (2.30), the obtained numerical results for all Reynolds numbers considered here are shown in Table 1. To compare with available results in the literature, the numerical results of Lange et al. [11] and Bharti et al. [12] are also listed in the table. It is clear that the results from current method basically agree well with the reference data. Besides, the value of average Nusselt number increases with Reynolds number as expected.Table 1
Nusselt number comparison at different Reynolds numbers.
ReferencesCasesNuLange et al. [11]Re=202.409Re=403.281Re=804.571Re=1506.396Bharti et al. [12]Re=202.465Re=403.283Re=80—Re=150—PresentRe=202.472Re=403.281Re=804.611Re=1506.432


## 4. Conclusion

In this paper, a hybrid immersed boundary-lattice Boltzmann method is developed to simulate the heat transfer problems. Employing the newly proposed IB-LBM [6], the nonslip boundary condition is enforced in simulation, and thus the velocity field can be accurately simulated. In the meanwhile, the temperature field is resolved by using the conventional energy equation with additional heat source term which mimics the influence of boundary on temperature field. Similar to the idea of velocity correction in [6], the heat source term is equivalent to the temperature correction and is set as unknown. By enforcing the thermal boundary condition, this unknown variable can be determined. Furthermore, utilizing the relationship between the temperature correction and heat flux, a simple approach for local and average Nusselt number evaluation is proposed. Different from the conventional way, there is no requirement for temperature gradient calculation.The efficiency and capability of developed method as well as way to calculate the local and average Nusselt numbers are illustrated by the simulation of forced convection problems. The obtained numerical results show good agreement with available data in the literature. It is demonstrated that the present method has a promising potential for solving thermal flow problems with complex geometry.



---

*Source: 161484-2012-03-15.xml*