An iterative source correction based immersed boundary-lattice Boltzmann method for thermal flow simulations

An iterative source correction based immersed boundary-lattice Boltzmann method for thermal flow simulations

International Journal of Heat and Mass Transfer 115 (2017) 450–460 Contents lists available at ScienceDirect International Journal of Heat and Mass ...

2MB Sizes 6 Downloads 90 Views

International Journal of Heat and Mass Transfer 115 (2017) 450–460

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage:

An iterative source correction based immersed boundary-lattice Boltzmann method for thermal flow simulations Jiayang Wu a,b,c, Yongguang Cheng a,⇑, Laura A. Miller b a

State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, PR China Department of Mathematics, University of North Carolina, Chapel Hill, NC 27599, United States c Key Laboratory of Hydraulic Machinery Transients, Ministry of Education, Wuhan University, Wuhan 430072, PR China b

a r t i c l e

i n f o

Article history: Received 29 April 2017 Received in revised form 18 July 2017 Accepted 2 August 2017

Keywords: Heat transfer Immersed boundary method Lattice Boltzmann method Non-slip boundary condition Iterative correction Fluid-structure interaction

a b s t r a c t Temperature jump at the boundary occurs when the conventional immersed boundary-lattice Boltzmann (IB-LB) method is applied to simulating the near boundary flows with heat transfer. To remedy this problem, an iterative correction is proposed to modify the heat source term in the IB-LB method. The source term in the LB equation is treated by Cheng’s scheme, in which the heat source at the next timestep is taken as unknowns and iteratively corrected until the resulting boundary temperature matches its desired value. Typical verification cases, including the two-dimensional (2D) heat transfer between two horizontal plates, the natural convection between two concentric circular cylinders, and 2D sedimentation of a single particle with heat convection are simulated to analyze the accuracy of the method. It is shown that the boundary temperature jump can be effectively removed for a certain range of LB relaxation time s, while the first-order spatial convergence of the IB method is still maintained. Also, a theoretical analysis is conducted based on the case of heat transfer between two plates. It is shown that the proposed method outperforms the widely-used direct source method in treating the Dirichlet boundary conditions when s is smaller than 1.624. To further demonstrate its capability for resolving complicated fluid-structure interaction problems, a three-dimensional sedimentation of a single particle in a vertical channel is analyzed. We find that the thermal convection may fundamentally affect the way the particle interacts with the surrounding fluid. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction The simple fact that both the immersed boundary (IB) method and the lattice Boltzmann (LB) method work on the same regular Eulerian grid makes the IB-LB coupling method efficient, which attracts considerable attention in the past decades. Applications such as bio-prosthetic heart valves [1], parachute opening dynamics [2], and flexible plate dynamics [3] have exhibited its effectiveness and great potential in simulating realistic and complex fluidstructure-interaction (FSI) problems. However, the non-slip boundary condition cannot be fully satisfied by the conventional IB method [4,5], in which the direct-forcing [6] or momentum exchange method [7] is applied to computing the boundary force. Great effort have been put into alleviating the boundary velocity slip problem. Wu [8] proposed an implicit velocity correction procedure, in which the force densities at all boundary points are considered as unknowns and calculated by solving a large banded ⇑ Corresponding author. E-mail addresses: [email protected] (J. Wu), [email protected] (Y. Cheng). 0017-9310/Ó 2017 Elsevier Ltd. All rights reserved.

matrix. Based on direct-forcing scheme, Wang [9] introduced a multi-direct forcing technique, in which the velocity field around the IB is corrected multiply times to satisfy the desired boundary condition. Recently, a more accurate IB-LB method was proposed by Zhang [10], in which the Lagrangian body forces are formulated by Cheng’s scheme [11] and corrected in an iterative manner. Up to now, the non-slip problem has been well resolved and many new applications have been reported [12,13], demonstrating IB-LB method’s capability in simulating curved and moving boundaries. However, the previous studies of the IB-LB method were mainly focused on resolving the athermal flows. When it comes to simulating the heat transfer problems by thermal IB-LB method, one may find the temperature errors at the IB also exist, similar to the boundary velocity slip in the athermal cases. Using a standard direct source [14] or penalty IB method [15] would not only cause spurious mass exchange, but also provoke severe non-physical temperature jump across the boundary. The strategies aiming at reducing this non-physical temperature jump can share the same ground as those designed for eliminating the boundary velocity slip. Wang [16] extended their multi-direct forcing technique to


J. Wu et al. / International Journal of Heat and Mass Transfer 115 (2017) 450–460

the thermal flow simulations. With the Navier-Stokes (NS) equations discretized by the finite-difference method, their scheme can achieve second-order spatial convergence. Based on Wu’s implicit velocity correction method, Seta [17] simulated the natural convection between two concentric cylinders and found that the error at each boundary point can be reduced by increasing the number of Lagrangian points. Despite all this progress, less understood is the relation between the temperature jump and the relaxation time in the thermal LB method. For the simple athermal shearing flows, it is well understood that the velocity slip at the IB would be more profound if the LB relaxation time becomes larger, and the non-slip boundary condition can be exactly satisfied only when the relaxation time takes certain value [5]. However, a concrete studies of whether the same conclusion can be applied to the thermal case is imperative. In this paper, following the recently proposed iterative force correction for the IB-LB method, an iterative source correction for the thermal IB-LB is proposed for accurate simulations of heat transfer flows. Several simple natural convection and a complex three-dimensional (3D) particle sedimentation problems are simulated to verify the method. A theoretical analysis is also performed for the heat transfer between two horizontal plates. It is shown that (1) the proposed method can effectively eliminate the temperature jump for a certain range of relaxation time; (2) it outperforms the traditional direct source method when the relaxation time is smaller than 1.624; (3) the temperature gradient will always deviate from its bulk value no matter how careful the relaxation time is chosen, provided that the single-relaxation-time (SRT) LB model is used. The rest of the paper is organized as follows. In Section 2, an overall description of the proposed method, including the governing equations, the thermal IB-LB coupling, and the heat source correcting, is given. In Section 3, the accuracy and ability of eliminating non-physical temperature jump of the method are verified by three typical examples. Then, a 3D case of particle sedimentation is presented to demonstrate the capability in simulating complex thermal FSI problem. Finally a brief conclusion and some discussions are made to end the article. 2. The iterative source correction IB-LB method for thermal flows 2.1. Governing equations for thermal flows with moving rigid body In the IB method, an Eulerian description of the NS and energy equations is used for thermal hydrodynamics, and a Lagrangian description of curvilinear boundaries is used for dynamics of the rigid bodies immersed in the fluid. To make the governing equations dimensionless, the characteristic velocity U 0 and length L are introduced. The definition of these scales are problemdependent and will be determined later. The time and pressure are scaled by L=U 0 . and qf U 20 , respectively, with qf being the fluid density. The dimensionless temperature is defined as ðT  T 0 Þ=jT s  T 0 j ¼ ðT  T 0 Þ=MT, in which both T s and T 0 are characteristic temperatures in the system. Then the governing equations can be recast in dimensionless forms as [18,19]:

r  u ¼ 0;


@u 1 Gr þ u  ru ¼ rp þ Du  2 Teg þ f ; @t Re Re


@T 1 þ u  rT ¼ DT þ q; @t RePr Z Fðs; tÞdðx  Xðs; tÞÞds; f ðx; tÞ ¼ Cb

Z qðx; tÞ ¼


Q ðs; tÞdðx  Xðs; tÞÞds:


Eqs. (1) and (2) are the incompressible NS equations in which the Boussinesq approximation is applied to obtain the buoyancy force, while Eq. (3) is the energy equation in which the pressure work and viscous heat dissipation are neglected. Note that the static gravity term is absorbed into the pressure gradient for convenience. Eqs. (4) and (5) are the boundary reaction equations to the fluid for force and heat, respectively. Here, X; F and Q are the Lagrangian boundary position, force density, and heat source, respectively; x; f ; u; T; q and p are the Eulerian spatial coordinate, external force field, fluid velocity, temperature, heat source and pressure, respectively; eg is the unit vector in the direction of gravity; s is the curvilinear coordinates of the boundary, which is usually chosen as the arc length at the boundary’s zero-stress state. The interaction of the fluid and the immersed boundary is realized by the Dirac delta function dðrÞ, which is used to spread the boundary force F and the heat source Q to the nearby fluid in Eqs. (4) and (5), respectively. The dimensionless parameters are

Reynolds number Re ¼ Grashof number Gr ¼

U0 L




m Prandtl number Pr ¼ ; v



in which m; b and g are the fluid viscosity, thermal expansion coefficient and gravitational acceleration, respectively. v is the thermal diffusivity which can be computed by dividing the thermal conductivity k by fluid density and specific heat capacity cf as k=qf cf . For the ease of elucidating, we will occasionally refer to the Rayleigh number Ra ¼ GrPr as well. For a 3D rigid body immersed in the fluid such as those in Section 3.4, its motion is governed by the Newton-Euler equations in the body-fixed non-inertial reference frame as [20,21]:

@v þ X  ms v ¼ F f þ E; @t @M þ X  M ¼ K f þ T; @t




where ms ; v ; X and M are the mass, translational velocity, angular velocity, and angular momentum of the rigid body, respectively. F f and K f are the force and moment exerted by the surrounding fluid. All other forces and moments such as gravity and collision force are represented by E and T, respectively. Moreover, if the axes of the body-fixed frame coincide with the principal axes of inertia, M can be expressed as M ¼ Is X, where Is is the diagonal inertial matrix. Once the translational and angular velocity are obtained by solving Eqs. (6) and (7), the position of the rigid body can be determined by solving following system of equations for center of mass do and Euler angle Ho

d_ o ¼ Rv ;


_ o ¼ Q X; H


where R and Q are the transformation matrices defined by the Euler angle as follows 2

3 cos w cos h  sin w cos / þ cos w sin h sin / cos w sin h cos / þ sin w cos / 6 7 R ¼ 4 sin w cos h cos w cos / þ sin w sin h sin / sin w sin h cos /  cos w sin / 5;  sin h cos h sin / cos h cos /


ð3Þ ð4Þ


1 sin / sin h= cos h cos / sin h= cos h

6 Q ¼ 40 0

cos /

 sin /

sin /= cos h

cos /= cos h

3 7 5:



J. Wu et al. / International Journal of Heat and Mass Transfer 115 (2017) 450–460

Note that any quantity can be easily transformed from the inertial reference frame to the non-inertial reference frame and vice versa. Assuming X IRF is some vector in inertial reference frame, its counterpart in non-inertial reference frame X NIRF can be calculated by




with R ¼ R . More details about the motion of rigid body can be found in [20]. 2.2. A double distribution LB model for resolving the thermal flow In the incompressible limit, if the transport coefficients such as viscosity or thermal diffusivity are independent of the temperature, the NS equations and energy equation can be solved in a separate manner. In the present study, a double distribution LB model is used to solve Eqs. (1)–(3). Two sets of distribution functions, namely the density distribution functions and the energy distribution functions, are introduced and evolve according to their own LB equations. For two-dimensional (2D) thermal flow, the most common LB model is the one using a square lattice with nine discrete velocity directions (donated as D2Q9 model). The governing equation with SRT collision approximation reads [22]

g a ðx þ ea dt ; t þ dt Þ  g a ðx; tÞ ¼ f a ðx þ ea dt ; t þ dt Þ  f a ðx; tÞ ¼


s 1


ðg eq a ðx; tÞ  g a ðx; tÞÞ þ dt F a ;


ðf eq a ðx; tÞ  f a ðx; tÞÞ þ dt Q a ;


where g a ðx; tÞ and f a ðx; tÞ are the discrete density and energy distribution functions along ath particle velocity direction ea ; g eq a ðx; tÞ and ðx; tÞ are its corresponding equilibrium; s and s are the relaxf eq c a ation times and can be computed from m and v, respectively; dt is the discrete timestep; F a is the LB forcing term which may be calculated by [23]

   1 e u e u  Fa ¼ 1  xa 3 a 2 þ 9 a 4 ea  f ; 2s c c


# 3ea  u 9ðea  uÞ2 3u2 ; þ  c2 2c4 2c2


f eq 0 ¼ 


f eq 14

2T u  u ; 3 c2 " # T 3 3ea  u 9ðea  uÞ2 3u2 ; ¼ þ  þ 9 2 2c2 2c4 2c2


f eq 58

" # T 6ea  u 9ðea  uÞ2 3u2 3þ ¼ þ  2 : 36 2c2 2c4 2c

X f a: a


2.3. The iterative source correction IB method The iterative force correction based IB-LB method [10] is extended to the simulation of thermal flows by using Eq. (20). Unlike other existing correction methods that target at removing boundary velocity slip [16,17], both the heat source and the energy distribution functions will be corrected in the present study. The extended method is called the iterative source correction IB method. First, the evolution equation of f a with the heat source term Eq. (20) is rewritten as

f da ðx þ ea dt ; t þ dt Þ  f a ðx; tÞ ¼ Xa þ

dt xa d ½q ðx; tÞ þ qd ðx þ ea dt ; t þ dt Þ; 2 ð21Þ

where Xa ¼ ½ðf eq a ðx; tÞ  f a ðx; tÞÞ=sc , and the superscript d indicates the desired value when the boundary condition is exactly satisfied. In practice, because the discrete Dirac delta functions appearing in Eqs. (4) and (5) introduce interpolation error, the Eulerian heat source q would be underestimated, causing the computed distribution functions to deviate from its desired value. Here, qd ðx þ ea dt ; t þ dt Þ is taken as unknown and will be determined in such a way that the resulting distribution function f ðx þ ea dt ; t þ dt Þ could match its desired value. Replacing qd ðx þ ea dt ; t þ dt Þ by an intermediate value q ðx þ ea dt ; t þ dt Þ in Eq. (21), the following equation is obtained

f a ðx þ ea dt ; t þ dt Þ  f a ðx; tÞ ¼ Xa þ

dt xa d ½q ðx; tÞ þ q ðx þ ea dt ; t þ dt Þ: 2 ð22Þ

Subtracting it from Eq. (21) leads to the following expression

f a ðx þ ea dt ; t þ dt Þ  f a ðx þ ea dt ; t þ dt Þ d




X f ea g a þ ; u¼ 2 a

1 ½xa qðx; tÞ þ xa qðx þ ea dt ; t þ dt Þ; 2

where qðx; tÞ is the Eulerian heat source in Eq. (3). Through Chapman-Enskog expansion detailed in Appendix, we prove that with Q a given by Eqs. (20), (13) can recover the exact energy equation within the incompressible limit. However, since Eq. (20) involves q evaluated at both t and t þ dt , an iteration procedure is needed within each timestep.


Then the macroscopic density, velocity and temperature are defined as

X q ¼ ga;

Qa ¼


in which xa is the weighting factor [22] and c is the magnitude of discrete velocity ea in either x or y direction. Q a is the heat source term which need special treatment in the present thermal IB-LB method and will be discussed later; For D2Q9 model used here, the equilibrium distribution functions are

g eq 08 ¼ xa q þ

The above choices for the definition of equilibrium and macroscopic variables originate from the incompressible LB model, which decouples the density fluctuation from the momentum equation and in this way improves the numerical stability [22]. As previously pointed out, the heat source term Q a in Eq. (13) should be carefully chosen. Using a regular treatment may result in inaccuracy and introduce extra term to the macroscopic energy equation [11]. Furthermore, in order to apply the iterative source correction to enhance the desired boundary condition, Q a should be formulated in a way that both the present and the next timestep effect are taken into account [10]. Here, the scheme originally proposed by Cheng [11] is extended to the thermal case. The heat source term is formulated as

dt xa d ½q ðx þ ea dt ; t þ dt Þ  q ðx þ ea dt ; t þ dt Þ: 2


The summation of Eq. (23) over a yields

q ðx þ ea dt ; t þ dt Þ  q ðx þ ea dt ; t þ dt Þ d

ð18Þ ð19Þ


2½T d ðx þ ea dt ; t þ dt Þ  T  ðx þ ea dt ; t þ dt Þ ; dt


where T d is the desired temperature which satisfies the non-slip condition, while T  is the intermediate temperature obtained by

J. Wu et al. / International Journal of Heat and Mass Transfer 115 (2017) 450–460


applying the heat source q . Eq. (24) can be viewed as a correction procedure in which the intermediate heat source q is iteratively corrected and forced to match its desired value qd . Assuming that Eq. (24) is also valid at Lagrangian point, one could obtain

Q d ðs; tÞ ¼ Q  ðs; tÞ þ

2½T dl ðs; tÞ  T l ðs; tÞ ; dt


where T dl ðs; tÞ is the desired boundary temperature, and T l ðs; tÞ is the intermediate temperature defined on the IB, which can be interpolated from the nearby flow field. To maintain the symmetry between the spreading and the interpolation process, the Dirac delta function is also used to implement the interpolation. Thus, T l ðs; tÞ can be obtained by

T l ðs; tÞ ¼



T  ðx þ ea dt ; t þ dt Þdðx  Xðs; tÞÞdx:


Fig. 1. Numerical setup for the heat transfer between two horizontal plates.

It should be noted that at the beginning of each timestep, despite that q ðx þ ea dt ; t þ dt Þ is determined by the iterative correction detailed above, qd ðx; tÞ in Eq. (22) needs to be obtained in advance. In the present work, a direct sourcing method which is a natural extension of its original direct forcing method is applied to compute qd ðx; tÞ, and a fourth order Runge-Kutta method is used to solve Eqs. (6) and (7). 3. Verification by typical cases Several benchmark cases are simulated to verify the accuracy and capability of the proposed method in resolving natural convection flows and thermal FSI problems. The accuracy compared with both analytical and previous numerical solutions, the spatial convergence order, the dependence of non-physical temperature jump on the LB relaxation time, and the comparison between present correction method with other existing IB method are studied in detail. Particularly, a theoretical analysis is carried out to investigate the boundary slip observed in the IB-LB framework. If otherwise specified, the Dirichlet boundary condition of velocity or temperature is realized by the non-equilibrium extrapolation method [24], while the Neumann boundary condition is handled by a second-order accurate method described in [25]. 3.1. Heat transfer between two horizontal plates The numerical setup is shown in Fig. 1. A uniform lattice of 201  201 is used to cover the computational domain with periodic boundary conditions applied at four sides. Two horizontal plates are placed at y ¼ 50 and y ¼ 150 with their velocity being zero. The temperatures of the upper and bottom plates are T d ¼ 1 and T d ¼ 1, respectively. The number of Lagrangian points for discretizing each plate is determined by making the distance between two adjacent IB points smaller than half the Eulerian grid spacing. In this case, the gravity term is dropped from the governing equations. The analytical temperature profile in the steady state will assume a linear distribution, which satisfies the exact boundary conditions at two plates.

Fig. 2. Spatial convergence order of present method for the heat transfer between two horizontal plates.

those yielded by other correction based IB-LB method. Despite that the LB method is second-order accurate in both time and space, the introduction of the first-order Dirac delta function in Eq. (4) and (5) will jeopardize the global accuracy even though only a finite layer is affected. Next, the relation between the temperature jump and the relaxation time is analyzed. As has been well documented in [5,17], the boundary velocity slip is strongly dependent on the LB relaxation time in the IB-LB method, and a larger relaxation time would cause the boundary velocity slip to be more prominent. The temperature profiles calculated with different sc are shown in Fig. 3. When sc is

3.1.1. Accuracy verification First, the spatial convergence order is accessed. The global relative error is

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uP u x2fluid ðTðx; tÞ  T d ðx; tÞÞ2 ; Error T ¼ t P 2 x2fluid Tðx; tÞ


where T d ðx; tÞ is the analytical solution. Here, both the lattice resolution and the IB grid are refined and the results are shown in Fig. 2. The first-order convergence is evident, which is consistent with

Fig. 3. Temperature distribution for different horizontal plates.

sc in the heat transfer between two


J. Wu et al. / International Journal of Heat and Mass Transfer 115 (2017) 450–460

relatively small, such as sc ¼ 1:0, the computed result matches the exact solution remarkably well. However, when sc is increased to 5.0, the simulated temperatures within the IB supporting layer starts to deviate from the analytical results. The temperature jump becomes even more severe as sc is further increased to 15.0. A close examination of the temperature profile within the IB layer reveals that when compared with the analytical solution, the simulated temperature is larger at the fluid node that coincides with the IB, while at the fluid nodes that are one or two lattices away from the IB, the obtained temperatures are lower. To get a better understanding of the observed phenomenon, a theoretical analysis is conducted. 3.1.2. A theoretical derivation for the temperature jump Following Le [5], when reaching its steady state, the velocity and temperature fields satisfy the following conditions.

u ¼ v ¼ 0;

@T ¼ 0; @x


@T ¼ 0; @t


where u and u are the x and y components of fluid velocity. In what follows, both the timestep dt and lattice spacing ea dt are set to unity for convenience. According to Eq. (13), (16) and (20), the energy distribution functions can be written as

f 0j ¼

4sc q j ; 9 j

f 1j ¼


T sc q þ ; 6 9 T sc q þ ; 6 9


f 2j ¼

T j1 sc  1 j1 1 j þ f þ ðq þ qj1 Þ; 18 6sc sc 2


T jþ1 sc  1 jþ1 1 j ðq þ qjþ1 Þ; ¼ þ f þ 18 6sc sc 4


T j1 sc  1 j1 1 j ðq þ qj1 Þ; þ f þ f 5j ¼ 72 12sc sc 5

T j  T jþ1 ¼


q0 2v


holds for jj  j0 j P 2. Eq. (33) implies that the temperature gradient is balanced by the applied heat source. However, it is only valid in the bulk region. Within the IB supporting layer, the temperature gradient would deviate from its desired value drastically, as indicated by Eqs. (31) and (32). Given that T ¼ 0 on the middle plane of the domain, starting from the middle plane and moving vertically towards the plate, T increases linearly with its gradient given by Eqs. (31)–(33). Thus, combining Eqs. (31)–(33) gives the temperatures at j0 and j0  1 as

  q0 h 1 1 þ 8sc  8s2c   ; v 4 4 24

  q0 h 1 1 þ 8sc  8s2c   ; v 4 2 48



where h is the distance between two plates. However, the exact boundary temperature T j0 can be derived by imposing a constant gradient i.e. Eq. (33) as

T j0 ¼


  q0 1 1 þ 8sc  8s2c  : v 2 48


Beyond the IB supporting layer (jj  j0 j P 2), the heat source is zero everywhere, and the temperature difference will remain constant.

T j0 1 ¼


  q0 1 1 þ 8sc  8s2c  ; v 4 48

T j0 þ1  T j0 þ2 ¼


f 3j ¼

f 4j

T j0  T j0 þ1 ¼

T j0 ¼ j


the total heat source spread by the IB and will be determined later. Using above relations, the temperature difference between adjacent fluid nodes is

hq0 : 4v


Thus, the boundary slip T js0 is given by

T js0 ¼ T j0  T j0 ¼



  1 1 þ 8sc  8s2c   : 4 24


On the other hand, according to Eq. (26), the temperature of IB

f 6j

f 7j

T j1 sc  1 j1 1 j ðq þ qj1 Þ; ¼ þ f þ 72 12sc sc 6 T jþ1 sc  1 jþ1 1 j ðq þ qjþ1 Þ; ¼ þ f þ 72 12sc sc 7

T jþ1 sc  1 jþ1 1 j þ f þ ðq þ qjþ1 Þ: f 8j ¼ 72 12sc sc 8


T l ¼ ð29hÞ


Here, the superscript j indicates the coordinate in y direction. Similarly, the same set of equations as Eq. (29) can also be derived for fluid nodes j þ 1 and j  1. Substituting these equations into Eq. (19) at fluid node j yields

0 ¼ vðT jþ1  2T j þ T j1 Þ þ þq : j

T l can be expressed as a linear combination of T j0 ; T j0 þ1 and T j0 1 as

  T j0 T j0 þ1 T j0 1 q0 h 3 1 þ 8sc  8s2c :   þ þ ¼ 2 4 4 v 4 8 32


In the present method, because the boundary temperature T l is iteratively corrected to match the desired value T d , this gives us a way to obtain q0 as

q0 ¼

32vT d : 8h  13  8sc þ 8s2c


Substituting Eq. (39) into Eq. (34) leads to the final expressions for the dimensionless boundary velocity T j0 =T d and temperature slip T js0 =T d as

1 þ 8sc  8s2c jþ1 ðq  2q j þ qj1 Þ 12 ð30Þ

Eq. (30) can be viewed as a finite difference discretization of the governing equation for our special case [17]. Here, we will focus on the temperature field near the plate where j takes the value j ¼ j0 (j0 ¼ 50 or j0 ¼ 150 in our case). Due to the symmetry, we have T j0 þ1 ¼ T j0 1 . The Dirac delta function in [26] is applied, suggesting that qj0 ¼ q0 =2; qj0 1 ¼ q0 =4 and qjjj0 j>1 ¼ 0, where q0 is

T j0 =T d ¼

4 6h  7  8sc þ 8s2c ; 3 8h  13  8sc þ 8s2c


T js0 =T d ¼

4 7  8sc þ 8s2c : 3 8h  13  8sc þ 8s2c


Now, it is evident that the non-physical temperature jump introduced by present thermal IB-LB method is dependent on sc . The only

sc that satisfies the equations T js0 =T d ¼ 0 and sc > 0:5 is

J. Wu et al. / International Journal of Heat and Mass Transfer 115 (2017) 450–460

Fig. 4. Analytical and simulated results with respect to non-physical temperature jump.


sc for the heat transfer between two horizontal plates. (a) Dimensionless boundary temperature. (b) Dimensionless

1.561. Also, the boundary slip would be more prominent as sc starts from 1.561 and approaches either 0.5 or infinite. Since in most cases sc works in the region ½0:5; 1:561, the above result shows that the dimensionless boundary slip could not be larger than 1:53% for present case. It should be pointed out that even if we can make sc close enough to 1:561 so that the desired boundary condition is approximated (T js0 < 10%), the temperature gradient within the IB layer would still deviate from its bulk value, since T j0  T j0 þ1 – q0 =2v and T j0 þ1  T j0 þ2 – q0 =2v. Fig. 4 compares the analytical solutions with the numerical results with respect to sc . We also include the results of direct sourcing method [17]. We can see that for both methods, the simulated curves match the theoretical predictions remarkably well, which validates our theoretical analysis and numerical implementation. Compared with the direct sourcing method, the nonphysical temperature jump can be reduced by present method when sc < 1:624. The two method yield almost identical results when sc approaches 0.5, and when it goes to infinite, the direct

sourcing method however is obviously the better choice. Moreover, Fig. 4 suggests that to make the present method accurate enough in approximating the boundary condition (T js0 < 10%), the relaxation time should not be larger than 3.528. 3.2. Natural convection between two concentric circular cylinders Whether the desired boundary condition can be enhanced by the IB-LB method may be affected by the orientation of the IB grid with respect to the LB lattice [27]. The natural convection between two concentric cylinders is simulated to verify its accuracy when the lattice orientation relative to the IB grid varies along the boundary. Here, a uniform 501  501 lattice is used. The radii of the inner and outer circular cylinders are set to Ri ¼ 62:5 and Ro ¼ 162:5, while their temperatures are T i ¼ 1 and T o ¼ 0, respectively. For comparison, all parameters are chosen so as to match those used in [18,28], namely the Rayleigh number Re ¼ 103  5  104 and the Prandtl number Pr ¼ 0:7. To make sure

Fig. 5. Temperature distributions for the natural convection between two concentric circular cylinders. (a) Re ¼ 103 , (b) Re ¼ 104 , (c) Re ¼ 105 . The top and bottom panels are from the present simulation and previous experiment [28], respectively.


J. Wu et al. / International Journal of Heat and Mass Transfer 115 (2017) 450–460

the fluid is incompressible, the Mach number is set to Ma ¼ U 0 =c ¼ 0:1. The characteristic scales are chosen as pffiffiffiffiffiffiffiffiffiffiffiffiffiffi L ¼ Ro  Ri ; DT ¼ Ri  Ro , and U 0 ¼ bg DTL. The temperature distributions under different Rayleigh numbers are shown in Fig. 5. To help elucidating, the polar angle h starting from the positive y direction and increasing clockwise is introduced. The simulated temperature distribution agrees well with previous study. At a given Rayleigh number, the temperature contours resemble those of the natural convection near a single horizontal cylinder. At the bottom of the inner cylinder (h ¼ 180 ), the thickness of thermal boundary layer is the smallest, indicating that the largest temperature gradient and heat flux occur at the bottom. The thickness of boundary layer increases gradually from the bottom to the top, where the buoyant plume above the inner cylinder impinges upon the outer cylinder, creating a high temperature gradient around the top surface. When Rayleigh number is small (Ra ¼ 103 ), the temperature variation with respect to h is small, indicating that the flow field is dominated by pure conduction. As Rayleigh number increases (Ra ¼ 104 ), the temperature distribution will get distorted, and the heat convection takes over. With further increase in Rayleigh number (Ra ¼ 5  104 ), the convection becomes more enhanced, resulting in the emerging of inversion region, in which the temperature near the outer cylinder is higher than that closer to the inner cylinder. The above results agree fairly well with those reported in [18,28]. Fig. 6 compares the simulated temperature and angular velocity profiles with those obtained by experiment [28] at Ra ¼ 5  104 . Good agreement can be found although they are not identical. The thermal boundary layer near both inner and outer cylinders are clearly seen, and for h ¼ 30 and h ¼ 90 , the temperature inversion mentioned above forces the curves to increase in the middle area of annulus. In Table 1, the distribution of local equiv-

alent conductivity keq on the inner cylinder is presented with respect to h. The local equivalent conductivity is defined as the actual heat flux divided by the heat flux that would occur by pure conduction in the absence of fluid motion [28]. For the inner cylinder, the highest value of keq occurs at the bottom, where the temperature gradient reaches its maximum as indicated in Fig. 5. The value decreases from the bottom to the top. Owing to the impinging plume, the temperature gradient on the inner cylinder reaches its minimum at the top, which results in a much lower keq . Generally speaking, the present results agree with the numerical solution of [28] and is consistent with Fig. 5. 3.3. Sedimentation of a 2D particle in an enclosure The sedimentation of a hot or cold particle is a very interesting case, because it provides abundant hydrodynamic and thermal flow structures depending on the Reynolds number and Grashof number. In this section, a test case originally elaborated by Feng et al. [29] is considered for the study of drag and terminal velocity of a 2D particle sedimentation with heat transfer. Initially, a circular particle with diameter d and Gr ¼ 100 is placed at 4d below the top of a rectangular enclosure with dimension ð16d; 40dÞ. The particle and boundary temperatures are both fixed throughout the simulation. According to Feng, the Reynolds number and the drag coefficient are defined as

Re ¼




CD ¼

FD ; 0:5qf V 2T d


where V T and F D are the simulated particle terminal velocity and drag force, which is measured when the particle reaches its steady state. Here, the Prandtl number and the fluid viscosity are set to 1 and 0.02, respectively. The particle to fluid density ratio .r is

Fig. 6. Dimensionless temperature and angular velocity in the annulus gap for the natural convection between two concentric circular cylinders. (a) Temperature. (b) Angular velocity.

Table 1 Comparison of local equivalent conductivity on inner cylinder at different Ra numbers. 0







T.H. Kuehn

















T.H. Kuehn present
















5  104

T.H. Kuehn
















Ra 3



J. Wu et al. / International Journal of Heat and Mass Transfer 115 (2017) 450–460

varying in the range ½1:001; 1:1 to yield a qr -dependent particle terminal velocity. Based on our simulations, the Reynolds number so defined is in the range ½5; 120. In the present simulation, both the hot and cold particle cases are considered. Initially, the temperature of the fluid and boundaries are set to 0, while the particle temperature is 1. The computation domain is discretized using a uniform lattice 321  801. In Fig. 7, the Drag coefficient are plotted against Reynolds number while the dimensionless terminal velocity V T =V ref is plotted against pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi the reference velocity V ref ¼ 0:5pdð.r  1Þg . For comparison, the results from [29] are included. At relatively low Reynolds number, due to the natural convection, the drag coefficient of the cold particle is larger than that of the hot particle, as expected. However, it is interesting to note that the C D difference between the hot and cold particles dies out gradually as Re increases. A similar trend can also be found in the variation of particle terminal velocity. The fluid around the cold particle has lower temperature and is thus heavier than that around the hot one, assisting in the particle downward motion and causing a larger terminal velocity. In general, our results are in good agreement with those reported by Feng. 3.4. Sedimentation of a 3D particle in an infinite channel In an original paper of Gan et al. [18], the sedimentation of a 2D particle with thermal convection is firstly simulated using an Arbitrary Lagrangian-Eulerian (ALE) based finite-element method. Up to five regimes have been detected for Gr < 104 . Different particle dynamics and wake structures are reported as Grashof number varies from one regime to another. In this section, we extend Gan’s work to 3D using the proposed thermal IB-LB method. The same set of characteristic scales as in [18] is adopted, namely the particle diameter d as characteristic length, the temperatures of the particle and the side walls as characteristic temperatures. The characpffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi teristic velocity is defined as U 0 ¼ 0:5pdð.r  1Þg with .r being the density ratio of particle to fluid. In the present study, the parameters are chosen as .r ¼ 1:00232; Re ¼ 40:5; Pr ¼ 0:7, and the channel width W ¼ 3d. The temperatures of the particle and the side walls are T s ¼ 1 and T 0 ¼ 0, respectively. The channel is long enough to eliminate possible effects of top and bottom boundaries. Also, to keep the particle within the computational domain, as the particle’s centroid reaches below one third of the channel height, both the flow field (distribution functions, velocity, temperature, etc.) and the particle position are shifted one lattice upward.


First, the sedimentation of either a cold or hot particle T s ¼ 1 at fixed Grashof number is analyzed. Fig. 8 shows the dimensionless particle settling velocity versus time at Gr ¼ 1100. As the simulation begins, the particle starts to move downward in the accelerated manner until it reaches the point where the gravity is balanced by the buoyancy force. Because the buoyancy force is dependent on the fluid density which in turn is a function of its temperature, the cold, hot and athermal particles behave in quite different ways. Compared with the cold particle case, the fluid temperature near the hot particle is much higher, which hinders the particle’s downward motion. As a result, when tU 0 =d < 5, the curve of the cold particle is steeper than that of the hot one, and the athermal curve hovers somewhere between them. The terminal settling velocity for the cold particle is the highest among all three cases, while the hot particle is the lowest. It is also interesting to note that for a hot particle, the settling velocity will undergo a slight increase right after its magnitude peaks. Such behavior could not be seen in the athermal or cold particle case, in which the settling velocity decreases monotonically until stabilizing at its steady equilibrium. When reaching its steady state, the distributions of dimensionless velocity and temperature are shown in Fig. 9. The thermal convection has a significant influence on the velocity field. For the cold particle, the velocity field resembles that of an athermal case, in which a pair of standing vortices trails in the wake. However, the

Fig. 8. Histories of dimensionless settling velocity for the sedimentation of a single particle with heat transfer.

Fig. 7. (a) Drag coefficient versus Reynolds number, and (b) dimensionless terminal velocity versus reference velocity for the 2D sedimentation of a single particle at Pr ¼ 1 and Gr ¼ 100.


J. Wu et al. / International Journal of Heat and Mass Transfer 115 (2017) 450–460

Fig. 9. Dimensionless velocity and temperature fields for the single particle sedimentation at Re ¼ 40:5; Gr ¼ 100 and Pr ¼ 0:7. (a) Velocity field for T s ¼ 1, (b) velocity field for T s ¼ 1, (c) temperature field for T s ¼ 1, (d) temperature field for T s ¼ 1.

Fig. 10. (a) Dimensionless velocity, and (b) temperature distributions for the sedimentation of a particle with natural convection between two side walls.

trailing vortices disappear for the hot particle case. Instead, another two pairs of vortices emerge, with one attached to the particle side and the other being at some distance behind the particle.

The disappearing of the trailing vortices may result from the fact that for the hot particle case, the natural and forced convection are in the same direction. The detachment of the boundary layer from the particle surface is suppressed by the hot fluid that carries upward momentum. With help of the upgoing hot current, the boundary layer eventually overcomes the reverse pressure gradient in the wake and remains attached to the particle surface. Next, the sedimentation of a particle with natural convection between side walls is simulated. To the best of our knowledge, similar study has not been reported yet. Such case can find its applications in chemical and food industries and also in solar collector [30]. Besides, Chiba et al. [31] have shown that natural convection of emulsion, which is a mixture of fluid and particles, has significant effects on the melting process of clathrate and is different from natural convection of pure substance. The simulation setup is the same as the previous case, except that the temperatures at the left and right walls are 0:5 and 0:5, respectively. The initial temperatures of the flow field and the particle are set to zero. Unlike the previous case, the thermal convection between two side walls will cause the flow field to be asymmetrical, and as a result, the particle would drift away from the centerline. Fig. 10 plots the velocity vectors and temperature distribution when the particle reaches its quasi-static state. Now, it is clear why the particle would migrate towards the wall with higher temperature. Due to the gravity, the fluid near the warmer

Fig. 11. Histories of x (a) and z (b) components of particle dimensionless velocities under different Reynolds numbers.


J. Wu et al. / International Journal of Heat and Mass Transfer 115 (2017) 450–460

Fig. 12. Temperature profiles for the heat transfer between two horizontal plates calculated by the IB-MRT-LB model. (a) Direct sourcing method, (b) iterative source correction method.

wall goes up, while the fluid near the cooler wall goes down. Thus, a simple shear flow is formed in the channel, which makes the particle rotate in the negative y direction. According to Bernoulli’s principle, a pressure difference across particle’s horizontal diameter is generated, and the particle is pushed away from centerline towards the hotter wall. To the study the influence of Reynolds number on the particle settling behavior, both x and z component of particle velocities are plotted in Fig. 11. Since the Grashof number is fixed in our case, a lower Reynolds number leads to a higher thermal expansion coefficient, which in turn results in a stronger shearing flow. As a consequence, the curve corresponding to a lower Reynolds number is steeper, indicating that a stronger net force is applied. As the particle gets close to the wall, it slows down and finally settles after a couple of oscillations. The histories of particle settling velocity under different Reynolds numbers are quite different. Similarly, a lower Reynolds number gives a stronger buoyancy force when the particle enters into the warmer area. Thus, the lower the Reynolds number, the sooner the settling velocity increases after hitting the bottom, and the lower the final settling velocity. Note that when Re is below 40:5, an overshoot can be found right before the curve becomes horizontal. This may correspond to the oscillations observed in the histories of uxðtÞ under the same Re.

4. Conclusion and discussion In this paper, an iterative source correction based IB-LB method is proposed to simulate the FSI with thermal convection. The external heat source is incorporated into the LB equation through Cheng’s scheme, in which the heat source is iteratively corrected to satisfy the desired Dirichlet boundary condition. In the three benchmark simulations, good agreement between our results and those in previous works verifies the accuracy of the proposed method. In addition, a theoretical analysis is carried out to find the intrinsic relation between the observed boundary temperature jump and the LB relaxation time. Our result indicates that for a certain range of relaxation time, our method outperforms the traditional direct source method in term of enhancing the boundary condition. Finally, a 3D sedimentation of a single particle with heat transfer is simulated to demonstrate its capability in simulating complicated FSI problem. As pointed out in [4], introducing the two-relaxation-time (TRT) or multi-relaxation-time (MRT) collision technique into the LB model could significantly reduce the boundary velocity slip. Their conclusion has been numerically and theoretically verified by sev-

eral simple cases, including symmetric shear flow and Poiseuille flow. Thus, we may intuitively conjecture that applying the TRT or MRT LB model could erase the non-physical temperature jump as well. To check the veracity of our conjecture, the case of heat transfer between two horizontal plates is simulated again with the MRT-LB model given in [32]. Since we just need to solve the energy equation in this case, the MRT collision technique is applied to Eq. (14). In the MRT-LB model, instead of the velocity space spanned by the discrete velocity set, the collision process happens in the moment space according to

Sðmeq f a ðx þ ea dt ; t þ dt Þ  f a ðx; tÞ ¼ M1 b a ðx; tÞ  ma ðx; tÞÞ þ dt Q a ; ð43Þ where ma ðx; tÞ is the moments of the distribution functions, and M is the transformation matrix which can be used to compute the moments from the distribution functions. b S is the diagonal relaxation matrix which is related to the thermal diffusivity. Details about the thermal MRT-LB model can be found in [32], and will not be elaborated here. Fig. 12 shows the temperature profiles for the two-plates problem under different sc . The results yielded either by direct sourcing method or the iterative source correction method are greatly improved compared with those obtained by SRT-LB model (Fig. 3). Although, the direct sourcing method does not eliminate the non-physical temperature jump completely, the proposed iterative source correction method gives fairly good results. Moreover, the proposed method does not produce a sc dependent result any more, and the temperature jump at high sc has been effectively erased. Thus, it may be concluded that only when the correction-based IB method is combined with the MRTLB model, can the non-physical temperature jump be significantly reduced for any relaxation time. Therefore, it is imperative to derived the theoretical formulae for the correction based IB-MRTLB model to understand how exactly the temperature jump is reduced even for high sc . This may be the focus of our future work. Conflict of interest The authors declared that they have no conflicts of interest to this work. Acknowledgements This work was supported by the National Natural Science Foundation of China (NSFC, Grant Nos. 11172219 and 51579187) and


J. Wu et al. / International Journal of Heat and Mass Transfer 115 (2017) 450–460

the Specialized Research Fund for the Doctoral Program of Higher Education of China (Grant No. 20130141110013). Appendix A. Appendix: Chapman-Enskog derivation of the energy equation With e representing the Knudsen number and assuming that e ¼ dt , the following multiscale expansions are introduced [22], ð1Þ 2 ð2Þ 3 f a ¼ f ð0Þ a þ ef a þ e f a þ Oðe Þ;


@ @ @ @ ¼ þe þ e2 þ Oðe3 Þ: @t @t0 @t 1 @t2


By expanding both f a ðx þ ea dt ; t þ dt Þ in Eq. (14) and qðx þ ea dt ; t þ dt Þ in Eq. (20) through the Taylor series and substituting Eqs. (44) and (45) into Eq. (14), one can find the equations corresponding to different orders of e as follows: The equation to order of e0 : eq f ð0Þ a ¼ fa :


The equation to order of



 @ 1 ð1Þ þ ea  r f ð0Þ f þ xa q: a ¼ @t0 sc a The equation to order of


e2 :

   2 @ ð0Þ @ 1 @ fa þ þ ea  r f ð1Þ þ ea  r f ð0Þ a þ a @t 1 @t 0 2 @t 0   1 1 @ ¼  f ð2Þ þ þ e  r xa q: sc a 2 @t0 a


Substituting Eq. (47) into Eq. (48) yields

   @ ð0Þ 1 @ 1 ð2Þ fa þ 1  þ ea  r f ð1Þ f : a ¼ @t 1 2s @t 0 sc a


Keeping in mind that with the equilibrium defined by Eq. (16), P ð0Þ P ð0Þ a xa ¼ 1; a f a ¼ T and a ea f a ¼ uT. The summation of Eq. (47) over a by using Eq. (20) gives the macroscopic equation on e1 time scale


@T þ u  rT ¼ q: @t 0


Similarly, the summation of Eq. (48) over a yields the macroscopic equation on e2 time scale

@T 2 ¼ @t 1 3


 1 MT: 2


Therefore, the combination of Eq. (50) + e  Eq. (51) results in

@T þ u  rT ¼ vMT þ q þ Oðd2t Þ; @t


which also gives the expression for the thermal diffusivity as

2 3


 1 dt : 2


From Eq. (52), it can be seen that the present thermal LB method can approximate the exact energy equation and is second-order accurate in time. References [1] I. Borazjani, Fluid-structure interaction, immersed boundary-finite element method simulations of bio-prosthetic heart valves, Comput. Meth. Appl. Mech. Eng. 257 (2013) 103–116.

[2] Y. Kim, C.S. Peskin, 3D parachute simulation by the immersed boundary method, Comput. Fluids 38 (6) (2009) 1080–1090. [3] R. Hua, L. Zhu, X. Lu, Dynamics of fluid flow over a circular flexible plate, J. Fluid Mech. 759 (2014) 56–72. [4] T. Seta, R. Rojas, K. Hayashi, A. Tomiyama, Implicit-correction-based immersed boundary-lattice Boltzmann method with two relaxation times, Phys. Rev. E 89 (2) (2014) 023307. [5] G. Le, J. Zhang, Boundary slip from the immersed boundary lattice Boltzmann models, Phys. Rev. E 79 (2) (2009) 026701. [6] A. Dupuis, P. Chatelain, P. Koumoutsakos, An immersed boundary-latticeBoltzmann method for the simulation of the flow past an impulsively started cylinder, J. Comput. Phys. 227 (9) (2008) 4486–4498. [7] X. Niu, C. Shu, Y. Chew, Y. Peng, A momentum exchange-based immersed boundary-lattice Boltzmann method for simulating incompressible viscous flows, Phys. Lett. A 354 (3) (2006) 173–182. [8] J. Wu, C. Shu, Implicit velocity correction-based immersed boundary-lattice Boltzmann method and its applications, J. Comput. Phys. 228 (6) (2009) 1963– 1979. [9] Z. Wang, J. Fan, K. Luo, Combined multi-direct forcing and immersed boundary method for simulating flows with moving particles, Int. J. Multiphase Flow 34 (3) (2008) 283–302. [10] C. Zhang, Y. Cheng, L. Zhu, J. Wu, Accuracy improvement of the immersed boundary-lattice Boltzmann coupling scheme by iterative force correction, Comput. Fluids 124 (2016) 246–260. [11] Y. Cheng, J. Li, Introducing unsteady non-uniform source terms into the lattice Boltzmann model, Int. J. Numer. Meth. Fluids 56 (6) (2008) 629–641. [12] J. Wu, C. Shu, N. Zhao, F. Tian, Numerical study on the power extraction performance of a flapping foil with a flexible tail, Phys. Fluids 27 (1) (2015) 013602. [13] J. Wu, C. Shu, Particulate flow simulation via a boundary condition-enforced immersed boundary-lattice Boltzmann scheme, Commun. Comput. Phys. 7 (4) (2010) 793. [14] E. Fadlun, R. Verzicco, P. Orlandi, J. Mohd-Yusof, Combined immersedboundary finite-difference methods for three-dimensional complex flow simulations, J. Comput. Phys. 161 (1) (2000) 35–60. [15] Z. Feng, E.E. Michaelides, The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems, J. Comput. Phys. 195 (2) (2004) 602–628. [16] Z. Wang, J. Fan, K. Luo, K. Cen, Immersed boundary method for the simulation of flows with heat transfer, Int. J. Heat Mass Transfer 52 (19) (2009) 4510– 4518. [17] T. Seta, Implicit temperature-correction-based immersed-boundary thermal lattice Boltzmann method for the simulation of natural convection, Phys. Rev. E 87 (6) (2013) 063304. [18] H. Gan, J. Chang, J.J. Feng, H.H. Hu, Direct numerical simulation of the sedimentation of solid particles with thermal convection, J. Fluid Mech. 481 (2003) 385–411. [19] Z. Feng, E.E. Michaelides, Heat transfer in particulate flows with direct numerical simulation (DNS), Int. J. Heat Mass Transfer 52 (3) (2009) 777–786. [20] L. Landau, E. Lifshitz, Mechanics, third ed., Buttersworth-Heinemann, Burlington, MA, 2003. [21] T.I. Fossen, Guidance and Control of Ocean Vehicles, John Wiley & Sons Inc., 1994. [22] P. Lallemand, L. Luo, Theory of the lattice Boltzmann method: dispersion, dissipation, isotropy, Galilean invariance, and stability, Phys. Rev. E 61 (6) (2000) 6546. [23] Z. Guo, C. Zheng, B. Shi, Discrete lattice effects on the forcing term in the lattice Boltzmann method, Phys. Rev. E 65 (4) (2002) 046308. [24] Z. Guo, C. Zheng, B. Shi, Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method, Chin. Phys. 11 (4) (2002) 366. [25] M. Junk, Z. Yang, Outflow boundary conditions for the lattice Boltzmann method, Prog. Comput. Fluid Dyn. Int. J. 8 (1–4) (2008) 38–48. [26] J. Wu, Y. Cheng, C. Zhang, et al., Three-Dimensional simulation of balloon dynamics by the immersed boundary method coupled to the multiplerelaxation-time lattice Boltzmann method, Commun. Comput. Phys. 17 (05) (2015) 1271–1300. [27] Y. Cheng, L. Zhu, C. Zhang, Numerical study of stability and accuracy of the immersed boundary method coupled to the lattice Boltzmann BGK model, Commun. Comput. Phys 16 (1) (2014) 136–168. [28] T. Kuehn, R. Goldstein, An experimental and theoretical study of natural convection in the annulus between horizontal concentric cylinders, J. Fluid Mech. 74 (04) (1976) 695–719. [29] Z. Feng, E.E. Michaelides, Inclusion of heat transfer computations for particle laden flows, Phys. Fluids 20 (4) (2008) 040604. [30] Y. Itaya, N. Arai, M. Hasatani, Characteristics of radiative energy collection by a fine-particle semitransparent suspension (FPSS) used as a heat vehicle and storage dedium of volume heat-trap type solar collector, Trans. Chem. Eng. 11 (6) (1985) 662–667. [31] T. Chiba, M. Okada, K. Matsumoto, Melting process of clathrate in a rectangular cell, Trans. Jpn. Soc. Refrig. Air Condit. Eng. 9 (2011) 169–180. [32] T. Zhang, D. Che, Double MRT thermal lattice Boltzmann simulation for MHD natural convection of nanofluids in an inclined cavity with four square heat sources, Int. J. Heat Mass Transfer 94 (2016) 87–100.