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

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.

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𝑓_𝛼𝐱+𝐞_𝛼𝛿𝑡,𝑡+𝛿𝑡−𝑓_𝛼()1𝐱,𝑡=−𝜏𝑓_𝛼(𝐱,𝑡)−𝑓^eq_𝛼(𝐱,𝑡)+𝐹_𝛼𝐹𝛿𝑡,_𝛼=11−𝑤2𝜏^𝛼𝐞_𝛼−𝐮𝑐²_𝑠+𝐞_𝛼⋅𝐮𝑐⁴_𝑠𝐞_𝛼⋅𝐟,𝜌𝐮=_𝛼𝐞_𝛼𝑓_𝛼+12𝐟𝛿𝑡, (2.1)(2.2)(2.3) where 𝑓_𝛼 is the distribution function, 𝑓^eq_𝛼 is its corresponding equilibrium state, 𝜏 is the single relaxation parameter, 𝐞_𝛼 is the lattice velocity, 𝑤_𝛼 are coefficients which depend on the selected lattice velocity model, and 𝐟 is the force density which is distributed from the boundary force density. Here, the D2Q9 lattice velocity model is used, which is𝐞_𝛼=()(),()()0,0,𝛼=0,±1,00,±1,𝛼=1∼4,±1,±1,𝛼=5∼8, (2.4) and the corresponding equilibrium distribution function is𝑓^eq_𝛼(𝐱,𝑡)=𝜌𝑤_𝛼𝐞1+^𝛼⋅𝐮𝑐²_𝑠+𝐞_𝛼⋅𝐮²−𝑐_𝑠|𝐮|²2𝑐⁴_𝑠 (2.5) with 𝑤₀=4/9, 𝑤_𝛼=1/9, (𝛼=1∼4), and 𝑤_𝛼=1/36, (𝛼=5∼8). 𝑐_𝑠=1/√3 is the sound speed of this model. The relationship between the relaxation time and the kinematic viscosity of fluid is 𝜐=(𝜏−0.5)𝑐²_𝑠𝛿𝑡.

In order to satisfy the nonslip boundary condition, the force density 𝐟 in (2.2) and (2.3) should be set as unknown [6]. It can be calculated from the fluid velocity correction 𝛿𝐮. Furthermore, 𝛿𝐮 can be obtained from the boundary velocity correction 𝛿𝐮_𝐵. The final expression for 𝛿𝐮_𝐵 is𝐀𝐗=𝐁, (2.6) where 𝐗=𝛿𝐮¹_𝐵,𝛿𝐮²_𝐵,…,𝛿𝐮^𝑚_𝐵^𝑇,𝛿𝐀=¹¹𝛿₁₂⋯𝛿_1𝑛𝛿₂₁𝛿₂₂⋯𝛿_2𝑛𝛿⋮⋮⋱⋮_𝑚1𝛿_𝑚2⋯𝛿_𝑚𝑛𝛿^𝐵₁₁𝛿^𝐵₁₂⋯𝛿^𝐵_1𝑚𝛿^𝐵₂₁𝛿^𝐵₂₂⋯𝛿^𝐵_2𝑚𝛿⋮⋮⋱⋮_𝐵𝑛1𝛿^𝐵_𝑛2⋯𝛿^𝐵_𝑛𝑚,𝐔𝐁=^1𝐵𝐔²_𝐵⋮𝐔^𝑚_𝐵−𝛿₁₁𝛿₁₂⋯𝛿_1𝑛𝛿₂₁𝛿₂₂⋯𝛿_2𝑛𝛿⋮⋮⋱⋮_𝑚1𝛿_𝑚2⋯𝛿_𝑚𝑛𝐮^∗₁𝐮^∗₂⋮𝐮^∗_𝑛. (2.7)(2.8) Here, m is the number of Lagrangian (boundary) points, and 𝑛 is the number of surrounding Eulerian points. 𝛿𝐮^𝑙_𝐵(𝑙=1,2,…,𝑚) is the unknown velocity correction vector at the boundary points. 𝛿_𝑖𝑗=𝐷_𝑖𝑗(𝐱_𝑖𝑗−𝐗^𝑙_𝐵)Δ𝑥Δ𝑦 and 𝛿^𝐵_𝑖𝑗=𝐷_𝑖𝑗(𝐱_𝑖𝑗−𝐗^𝑙_𝐵)Δ𝑠_𝑙. 𝐷_𝑖𝑗(𝐱_𝑖𝑗−𝐗^𝑙_𝐵) is the delta function, which is expressed as𝐷_𝑖𝑗𝐱_𝑖𝑗−𝐗^𝑙_𝐵=𝛿𝑥_𝑖𝑗−𝑋^𝑙_𝐵𝛿𝑦_𝑖𝑗−𝑌^𝑙_𝐵,𝛿(𝑟)=1𝜋|𝑟|4ℎ1+cos,|𝑟||𝑟|2ℎ≤2ℎ,0,>2ℎ, (2.9) and Δ𝑠_𝑙 is the arc length of boundary element. 𝐔^𝑙_𝐵 is the boundary velocity. The intermediate fluid velocity 𝐮^∗ is calculated by𝐮^∗=1𝜌_𝛼𝐞_𝛼𝑓_𝛼. (2.10) By solving equation system (2.6), the unknown variables 𝛿𝐮^𝑙_𝐵 can be obtained. After that, the fluid velocity correction 𝛿𝐮 can be calculated by𝛿𝐮=_𝑙𝛿𝐮^𝑙_𝐵𝐷_𝑖𝑗𝐱_𝑖𝑗−𝐗^𝑙_𝐵Δ𝑠_𝑙. (2.11) As shown in (2.3), the relationship between the force density 𝐟 and the fluid velocity correction 𝛿𝐮 can be written as𝐟=2𝜌𝛿𝐮𝛿𝑡. (2.12) In the LBM, the macroscopic variables such as density and pressure are calculated by𝜌=_𝛼𝑓_𝛼,𝑃=𝑐²_𝑠𝜌. (2.13)

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) where 𝑇 is the temperature, 𝑐_𝑝 is the specific heat of the fluid, 𝑘 is the thermal conductivity of the fluid, and 𝑞 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) The solution of (2.15) is regarded as intermediate temperature 𝑇^∗. In the second step, or named Corrector step, the following equation is solved:𝜌𝑐_𝑝𝜕𝑇𝜕𝑡=𝑞. (2.16) Here, (2.16) means the effect of heat source on temperature field. Using the Euler explicit scheme, 𝜕𝑇/𝜕𝑡 can be rewritten as 𝑇^𝑛+1−𝑇^∗/𝛿𝑡=𝛿𝑇/𝛿𝑡, where 𝑇^𝑛+1 is the required temperature at next time step. Therefore, (2.16) becomes𝑞=𝜌𝑐_𝑝𝛿𝑇𝛿𝑡, (2.17) where 𝛿𝑇 is the temperature correction which determines the unknown variable 𝑞. 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], 𝛿𝑇 can further be calculated from the boundary temperature correction 𝛿𝑇_𝐵. By enforcing the constant temperature boundary condition, 𝛿𝑇_𝐵 can computed from the following equation system:𝐂𝐘=𝐃, (2.18) where 𝐘=𝛿𝑇¹_𝐵,𝛿𝑇²_𝐵,…,𝛿𝑇^𝑚_𝐵^𝑇,𝛿𝐂=¹¹𝛿₁₂⋯𝛿_1𝑛𝛿₂₁𝛿₂₂⋯𝛿_2𝑛𝛿⋮⋮⋱⋮_𝑚1𝛿_𝑚2⋯𝛿_𝑚𝑛𝛿^𝐵₁₁𝛿^𝐵₁₂⋯𝛿^𝐵_1𝑚𝛿^𝐵₂₁𝛿^𝐵₂₂⋯𝛿^𝐵_2𝑚𝛿⋮⋮⋱⋮_𝐵𝑛1𝛿^𝐵_𝑛2⋯𝛿^𝐵_𝑛𝑚,𝑇𝐃=^1𝐵𝑇²_𝐵⋮𝑇^𝑚_𝐵−𝛿₁₁𝛿₁₂⋯𝛿_1𝑛𝛿₂₁𝛿₂₂⋯𝛿_2𝑛𝛿⋮⋮⋱⋮_𝑚1𝛿_𝑚2⋯𝛿_𝑚𝑛𝑇^∗₁𝑇^∗₂⋮𝑇^∗_𝑛. (2.19)(2.20) Here, 𝑇^𝑙_𝐵 is the specified boundary temperature. The parameters 𝑚, 𝑛, 𝛿_𝑖𝑗, and 𝛿^𝐵_𝑖𝑗 have the same meaning as those in equation system (2.6). After 𝛿𝑇_𝐵 is obtained from equation system (2.18), the temperature correction 𝛿𝑇 can be calculated through𝛿𝑇=_𝑙𝛿𝑇^𝑙_𝐵𝐷_𝑖𝑗𝐱_𝑖𝑗−𝐗^𝑙_𝐵Δ𝑠_𝑙. (2.21) Consequently, the corrected temperature can be computed by𝑇=𝑇^∗+𝛿𝑇. (2.22)

2.3. Local and Average Nusselt Number Evaluation

In thermal flow problems, the Nusselt number is an important measurement parameter. The local Nusselt number Nu on the surface of boundary is defined asNu𝐗_𝐵=ℎ𝐗,𝑡^𝐵,𝑡𝐿_𝑐𝑘, (2.23) where ℎ is the local convective heat transfer coefficient and 𝐿_𝑐 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) where 𝑇_∞ is the free stream temperature. Substituting (2.24) into (2.23) givesNu𝐗_𝐵𝐿,𝑡=−^𝑐𝑇_𝐵−𝑇_∞𝜕𝑇𝐗𝜕𝑛^𝐵,𝑡. (2.25) Further averaging over the entire boundary, we can have the surface overall mean Nusselt number Nu1Nu=𝐿𝐿Nu𝑑𝑠=−^𝑐𝑇_𝐵−𝑇_∞𝐿𝜕𝑇𝐗𝜕𝑛^𝐵,𝑡𝑑𝑠, (2.26) where 𝐿 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 𝑄:𝑄𝐗_𝐵,𝑡=−𝑘𝜕𝑇𝐗𝜕𝑛^𝐵,𝑡. (2.27) At the same time, similar to (2.17), the boundary heat flux can be calculated using the boundary temperature correction𝑄=𝜌𝑐_𝑝𝛿𝑇_𝐵𝛿𝑡. (2.28) Substituting (2.27) and (2.28) into (2.25) and (2.26), the expressions for the local and average Nusselt numbers can be rewritten asNu=𝜌𝑐_𝑝𝐿_𝑐𝑘𝛿𝑡𝑇_𝐵−𝑇_∞𝛿𝑇_𝐵,Nu=𝜌𝑐_𝑝𝐿_𝑐𝑘𝛿𝑡𝑇_𝐵−𝑇_∞𝐿_𝑙𝛿𝑇^𝑙_𝐵Δ𝑠_𝑙. (2.29)(2.30) 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=𝜌𝑈_∞𝐷/𝜇 and Prandtl number Pr=𝜇𝑐_𝑝/𝑘. Here, 𝑈_∞ is the free stream velocity, 𝐷 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 is 60𝐷×40𝐷 with the nonuniform mesh size of 321×241. The cylinder is located at (20𝐷,20𝐷). The region around the cylinder is 1.2𝐷×1.2𝐷 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 at Re=20.

For the unsteady flow simulation, the computational domain is 50𝐷×40𝐷 with the mesh size of 401×301. The region around the cylinder keeps 1.2𝐷×1.2𝐷 with a uniform mesh size of 97×97. The cylinder is still located at (20𝐷,20𝐷). 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.

(a)

(b)

Figure 3 

Instantaneous isotherms (colorful) and vorticity contours (gray) for unsteady forced convection from a stationary cylinder.

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.

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.