Lattice Boltzmann simulation of thermal flows with complex geometry using a single-node curved boundary condition

International Journal of Thermal Sciences 146 (2019) 106112

Contents lists available at ScienceDirect

International Journal of Thermal Sciences journal homepage: http://www.elsevier.com/locate/ijts

Lattice Boltzmann simulation of thermal flows with complex geometry using a single-node curved boundary condition Tao Shi*, Baiman Chen, Hanmin Xiao, Simin Huang Key Laboratory of Distributed Energy Systems of Guangdong Province, Dongguan University of Technology, Dongguan, 523808, China

A R T I C L E I N F O

A B S T R A C T

Keywords: Lattice Boltzmann method Thermal flows Curved boundary condition Single-node scheme Second-order accuracy

Thermal curved boundary condition is proposed for lattice Boltzmann simulation of thermal flows with complex geometry. In this scheme, the unknown temperature distribution function at the boundary point is interpolated by the distribution functions at the interface and a further fluid point. The former one is decomposed into the equilibrium and non-equilibrium parts. The equilibrium part can be determined by the boundary conditions, while the anti-bounce back rule is applied to obtain the non-equilibrium part, i.e., the key point for the present scheme. On the other hand, the latter one streams directly from the boundary point. Since information of a single fluid point is needed, the computation involved in the present scheme is completely localized. Another advantage of the method is for the second-order convergence rate in simulations, preserving the overall accuracy of the lattice Boltzmann method. Four thermal problems are used to validate the present scheme. The first two cases, i. e., thermal Couette flow with wall injection, and natural convection in a square cavity are performed to verify the second-order accuracy. It is further applied to the other three thermal flows with complex geometries, including the natural convections in a cavity with a cylinder and in a concentric annulus, and a hot cylinder moving in a channel. Good agreement can be found between the present results and those available in the literature. The five test cases together demonstrate that the present method retains the property of local computation and secondorder accuracy simultaneously in the simulations.

1. Introduction Numerical investigation of the flow with heat transfer is important for a variety of applications in the fields of energy development, chemical engineering, environmental governance and so on. That nonisothermal flow problem is usually described by the Navier-Stokes (NS) coupled with the convection-diffusion (CD) equations. The NS-CD equations used to be solved by the traditional computational fluid dy­ namics (CFD) methods. In the past two decades, the lattice Boltzmann method (LBM) has been advanced into an effective tool for various isothermal flow problems [1,2]. The good parallel performance and the particle-based background together make it promising for the complex non-isothermal flows [3]. There are generally three approaches to incorporate the thermal ef­ fect into the LBM: the multi-speed model [4], the hybrid scheme [5] and the double-population method [6]. The first method uses additional velocity vectors for recovering the macroscopic NS-CD equations, and hence higher order velocity terms are usually included in the equilib­ rium distribution function. The most attractive feature of the

multi-speed scheme is only one set of distribution function needed in the simulations. However, it was reported to possess poor numerical sta­ bility [7]. Such deficiency can be improved in the latter two methods, which belong to the passive scalar model. That is to say, both for the double-population and hybrid schemes, the fluid motion and heat transfer are separately solved by the NS and CD equations. The differ­ ence lies in that, those two methods apply the traditional CFD method and another set of distribution function to handle the CD equation, respectively. For the integrity in taking advantage of the LBM, the double-population scheme is adopted for the simulations of flow with heat transfer in this study. Boundary condition should be specified, i.e., identifying the un­ known distribution function (UDF), in the investigation of the fluid flows using LBM [8]. Before the discussion of the boundary condition for the thermal flows, it is beneficial to first introduce the isothermal counterpart. It is known that one of the most popular boundary schemes is the simple bounce back (BB) rule proposed by Ladd [9]. In such method, the UDF can be directly determined from the distribution function in the opposite direction. However, with a simplification to

* Corresponding author. E-mail address: [email protected] (S. Tao). https://doi.org/10.1016/j.ijthermalsci.2019.106112 Received 19 December 2018; Received in revised form 18 September 2019; Accepted 20 September 2019 Available online 26 September 2019 1290-0729/© 2019 Elsevier Masson SAS. All rights reserved.

S. Tao et al.

International Journal of Thermal Sciences 146 (2019) 106112

Fig. 3. Schematic of the diffusion in the plane shear flow with wall injection. Both walls have a retraction of qδx.

[10,11]. The development of thermal CBC is basically the same with the isothermal counterpart outlined above. Zhang et al. [12] first proposed a boundary scheme based on the simple BB, which was used to simulate complex flows with heat and mass transfer in the porous media. How­ ever, since the relative distance from the boundary point to the interface was not involved in the construction, it might lose accuracy for general curved interface. Subsequently, that scheme by Zhang et al. [12] was advanced by Chen et al. [13]. The authors argued that the real tem­ perature at the middle of the lattice link should be used, which was determined by interpolation or extrapolation for different cases. With the parameter of relative distance incorporated in their scheme, improved results had been reported in the simulations [13]. On the other hand, Lin et al. [14] presented a thermal boundary method for complex flows, where the UDF is a combination of the known distribution func­ tion and a corrector. A second fluid node should be used in the linear interpolation for obtaining the correction term. Furthermore, Fattahi et al. [15] and Khazaeli et al. [16] respectively extended the non-equilibrium extrapolation method by Guo et al. [10] and ghost fluid scheme by Tiwari et al. [11] for the isothermal flows to develop the corresponding thermal CBC versions. A common disadvantage of those two methods [15,16] is that they are inconsistent around the threshold value of 0.5 about the relative distance. Recently, there were some other thermal CBCs [17–20] with more complicated interpolation or extrap­ olation computations involved. Those described thermal CBCs [13–20] all dramatically improve the accuracy to be second-order in the treatment of boundaries with com­ plex geometry. However, most of them need to obtain the UDF at the boundary node with the aid of the information of one or two and more neighboring fluid nodes. It is inevitable to undermine the local computation property and further the parallel performance of LBM [21]. However, it should be noted clearly that, the local computation property of the original LBM has been completely preserved in the simple BB scheme by Zhang et al. [12], because there is no extra fluid node involved in that boundary treatment. Therefore, it can be observed that, most methods of boundary treatment cannot integrate the two proper­ ties of the accuracy and local computation simultaneously in the simu­ lation of thermal flows using LBM. Therefore, the development of a single-node and second-order accurate thermal CBC seems is desirable in the thermal flow simulations, which is the motivation of the present study. The remaining part of the paper is organized as follows. A methodology introduction is provided in Section 2. The numerical scheme is validated in Section 3. Finally, conclusions are summarized in Section 4.

Fig. 1. Schematic of the present single-node thermal boundary condition. xA is the boundary node which has UDFs. xW and xC are the intersection point and the nearest fluid point along the intersection direction, respectively.

Fig. 2. Illustration of the evaluation of the Nusselt number for the straight (a) and curved walls (b). xB is a boundary node, and BE is a line perpendicular to the solid surface with xE the foot point.

locate the interface at the center of a lattice link, the scheme inevitably models the original smooth particle surface to a series of zig-zag lines. That renders a collapse of geometry integrity and results in a deficiency of low precision in the boundary treatment, which is insufficient in the simulation of complex flows [8]. To improve the accuracy of the simple BB scheme, many refined methods, called curved boundary condition (CBC) are then developed for considering the real shape of the interface

2. Numerical methodology The incompressible viscous non-isothermal flow is generally gov­ erned by the following NS-CD equations,

2

S. Tao et al.

International Journal of Thermal Sciences 146 (2019) 106112

Fig. 4. The temperature profiles at different Re and Pr with q ¼ 0.2: (a) Pr ¼ 0.7; (b) Re ¼ 10.

Fig. 5. The relative errors of the temperature at different Re and Pr: (a) q ¼ 0.2; (b) q ¼ 0.8.

r:u ¼ 0;

(1)

∂u 1 rp þ νr2 u þ ab ; þ u:ru ¼ ∂t ρ

(2)

∂T þ r:ðuTÞ ¼ κr2 T; ∂t

(3)

where ρ, u, p, T, ν and κ are the density, velocity, pressure, temperature, viscosity and thermal diffusivity of fluid, respectively. ab is the body force in the momentum equation, which denotes the buoyancy force and can be expressed with the Boussinesq approximation [22] as ab ¼ g0 βðT

(4)

T0 Þj;

where g0 is the gravity acceleration, β is the thermal expansion coeffi­ cient, T0 is the average temperature of fluid, and j is the unit vector in the vertical direction. 2.1. Lattice Boltzmann method with body force

Fig. 6. Schematic of the natural convection in a square hollow enclosure, which shifts to the right side with a distance of qδx.

The double-population scheme is applied to solve the above NS-CD equations, which is capable for a wide range of temperature variations and has good numerical stability [23,24]. Two set of distribution func­ tions are adopted in that approach, b 1 X

fi ðx þ ei δt ; t þ δt Þ

fi ðx; tÞ ¼

M 1 SM

� ij

fj

� fjeq þ δt Fi ;

i

j¼1

¼ 0; 1; …; b

3

1:

(5)

S. Tao et al.

International Journal of Thermal Sciences 146 (2019) 106112

Fig. 7. Isotherms (left column) and streamlines (right column) in the square cavity: (a) and (b), Ra ¼ 1.0 � 103; (c) and (d), Ra ¼ 1.0 � 105; (e) and (f), Ra ¼ 1.0 � 106. Table 1 Comparison of the uy,max and Numax for the thermal flow in a hollow cavity. Ra

1.0 � 103

1.0 � 104

uy,maxL/κ

Numax

uy,maxL/κ

Numax

uy,maxL/κ

Numax

uy,maxL/κ

Numax

Davis et al. [35] Present

3.697 3.701

1.505 1.509

19.617 19.632

3.328 3.461

68.59 68.80

7.717 7.683

219.36 217.85

17.925 17.308

b 1 X

gi ðx þ ei δt ; t þ δt Þ

gi ðx; tÞ ¼

M 1 SM

� ij

gj

� geq j ;

1.0 � 105

distribution functions, with velocity ei at position x and time t. δt denotes the temporal step, and b is the total number of discrete velocities. M is a b � b transform matrix and S a relaxation matrix. f eq and geq are the Maxwellian-type equilibrium distribution functions. The force term of Fi is used to incorporate the buoyancy force ab. Note that equations (5) and (6) are solved generally by the following two steps of collision and

i

j¼1

¼ 0; 1; …; b

1:

1.0 � 106

(6)

where fi(x, t) and gi(x, t) are respectively the fluid and temperature 4

S. Tao et al.

International Journal of Thermal Sciences 146 (2019) 106112

Fig. 8. Schematic of the temperature profile at the horizontal centerline (a) and the errors of temperature (b).

For simplicity and without loss of generality, the D2Q9 velocity model (two dimensions with nine lattice velocities, b ¼ 9) is employed in the present study, where the velocity set and the corresponding weight coefficients are given as 8 8 < ð0; 0Þ; < 4=9; i ¼ 0; 1Þ π =4; sinði 1Þ π =4�; ω ¼ 1=9; i ¼ 1; 2; 3; 4; ei ¼ ½cosði i : pffiffi : 1=36; i ¼ 5; 6; 7; 8: 2½cosði 1Þπ=4; sinði 1Þπ=4�; (10) The transform matrix M is given by Ref. [26] and the relaxation matrix is � � � � � � � � � � (11) S ¼ diag 1 τρ ; 1 τe ; 1 τε ; 1 τd ; 1 τq ; 1 τd ; 1 τq ; 1 τf ; 1 τs : where τ is the relaxation rate and relaxation time; (τρ, τd) and the remaining relaxation times are for conserved and non-conserved mo­ ments. Note that τρ and τd can take any non-zero values, while τe, τε and τq should be set greater than 0.5, and are determined following Ref. [26]. τf should be τT for the temperature field. Through Chapman-Enskog expansion, the macroscopic quantities, i. e., ρ, u, p, ν, T and κ can be respectively obtained as , 8 8 X X ρ¼ fi ; ρu ¼ ci fi þ δt ag 2; p ¼ ρc2s ; ν ¼ c2s ðτs 0:5Þδt ; T

Fig. 9. Schematic of the geometry arrangement and boundary conditions for the thermal flow in a cavity with a concentric cylinder.

i¼0

¼

0:5Þδt : (12)

2.2. Single-node curved boundary condition for thermal flow

identified by the density ρ, velocity u and temperature T of the fluid, " # # " �2 �2 ej :u ej :u ej :u u2 ej :u u2 eq eq fj ¼ ωj ρ 1 þ 2 þ ; gj ¼ ωj T 1 þ 2 þ ; cs cs 2c4s 2c2s 2c4s 2c2s

Thermal curved boundary condition (CBC) should be applied in the simulation of a flow with complex geometry. For LBM, the main concern is to determine the unknown distribution function (UDF) at the boundary point, which has at least a lattice link with the nodes at the opposite side of the interface. Distribution functions in some directions for the boundary point are unknown after the streaming step. Taking point A for example that illustrated in Fig. 1, the f 1, f 2 and f 5 are UDFs. As outlined previously, the simple BB method and the CBC scheme usually have the deficiencies of low precision and non-local computa­ tion, respectively. Therefore, this paper introduces a one-point secondorder CBC for both the fluid and temperature fields, where the UDF is obtained as

(7) pffiffiffiffiffiffi where ωj is the weight coefficient; cs ¼ RT is the lattice sound speed pffiffiffi and usually set to be c/ 3 for isothermal flows (R is the gas constant, and c ¼ δx/δt ¼ 1 in this paper with δx the lattice spacing). For consis­ tency, the forcing term Fi in Eq. (5) should be given by Refs. [25,29]. S=2ÞMF;

gi ; κ ¼ c2s ðτT

i¼0

streaming [7]. Furthermore, the multi-relaxation-time (MRT) collision model is used here instead of the lattice Bhatnagar-Gross-Krook (LBGK) to avoid the unphysical numerical artifact and improve the stability [24]. The equilibrium distribution function fjeq and gjeq are typically

F ¼ M 1 ðI

i¼0 8 X

(8)

where I is the identity matrix, F ¼ (F0, F1, …, Fb-1), and F ¼ ðF0 ; F1 ; …; Fb 1 Þ with �� � ei :ag uag : ei ei c2s I ; i ¼ 0; 1; …; b 1: (9) Fi ¼ w i þ c4s c2s

gi ðxA ; t þ δt Þ ¼

1 geq ðxW ; t þ δt Þ 1þq i

� gneq ðxA ; tÞ þ qg*i ðxA ; tÞ ; i

(13)

where the distance ratio q ¼ |xA – xW|/|xA – xC| denotes the relative distance from the boundary point to the solid surface. 5

International Journal of Thermal Sciences 146 (2019) 106112

S. Tao et al.

Fig. 10. Isotherms (left column) and streamlines (right column) in the enclosure with D/L ¼ 0.4: (a) and (b), Ra ¼ 1.0 � 104; (c) and (d), Ra ¼ 1.0 � 105; (e) and (f), Ra ¼ 1.0 � 106. Table 2 Comparison of the surface-averaged Nu of the inner concentric cylinder. Ra

1.0 � 104

D/L Shu et al. [30] Lin et al. [14] Present

0.2 2.082 2.081 2.090

1.0 � 105 0.4 3.245 3.229 3257

0.6 5.395 5.183 5.077

0.2 3.786 3.801 3.796

1.0 � 106 0.4 4.861 4.924 4.946

6

0.6 6.214 6.147 6.281

0.2 6.106 6.108 6.112

0.4 8.898 9.680 8.938

0.6 12.000 11.954 12.002

S. Tao et al.

International Journal of Thermal Sciences 146 (2019) 106112

It can be obviously found that, there is only one fluid node involved in the present boundary scheme. Furthermore, the two approximations involved in the present approach, i.e., Eq. (18) and Eq. (19) had been demonstrated through Chapman-Enskog analysis to be first- and secondorder accurate in the literature [12,27,28]. As a result, the non-equilibrium distribution functions are obtained being at least first-order accuracy, which is enough to derive the second-order con­ struction of UDFs. Therefore, the present single-node thermal CBC is of second-order accuracy in space theoretically. Note that the no-slip boundary condition should also be imple­ mented at the solid surface in the flow with heat transfer. The methods proposed by Guo et al. [10] and Tiwari et al. [11], among others [36,37] can be applied for such purpose. In this study, the single-node CBC [38] for the isothermal flows is adopted for consistency. 3. Results and discussions The present thermal CBC introduced above obviously indicates itself as a single-node boundary scheme. And theoretically, it is second-order accurate spatially, which is to be validated numerically in this section with the thermal Couette flow with wall injection, and natural convec­ tion in a square cavity. The boundary scheme is further verified in the applications to the complex thermal flows, including the natural con­ vection in a cavity with a cylinder, concentrically or eccentrically, the convective heat transfer in a concentric annulus, and a hot cylinder moving in a channel. The thermal flow is generally controlled by several parameters, such as the Reynolds, Prandtl and Rayleigh numbers, which are defined respectively as

Fig. 11. Schematic of the geometry arrangement and boundary conditions for the thermal flow in a cavity with an eccentric cylinder.

For the derivation of the formula given by Eq. (13), it is clear that the UDFs are the combination of distribution functions at the node xC and intersection point xW in the present boundary scheme, gi ðxA ; t þ δt Þ ¼ as

1 ðgi ðxW ; t þ δt Þ þ qgi ðxC ; t þ δt ÞÞ: 1þq

Re ¼

(14)

(15)

Therefore, the remaining difficulty is to determine the thermal dis­ tribution function at xW. In this study, inspired by Refs. [10,27], the distribution function is divided into two parts of equilibrium and non-equilibrium, neq gi ðxW Þ ¼ geq i ðxW Þ þ gi ðxW Þ:

∂

(16)

Nu ¼

t þ δt Þ ¼

gneq i ðxA ;

(17)

t þ δt Þ:

gneq ðxA ; tÞ; i

(18)

gneq ðxA ; tÞ: i

(21)

T Tc Th Tc

∂n

jwall ;

(23)

3.1. Thermal Couette flow with wall injection

(19)

The first test case is the thermal Couette flow with wall injection. As illustrated in Fig. 3, the flow is in a rectangular domain with width L. Since periodic boundary condition is specified in the horizontal direc­ tion, the length of the domain should have negligible influence on the simulation results. Furthermore, the upper wall of channel with a higher temperature T ¼ Th moves constantly at the velocity u ¼ Ux. The cold lower wall with temperature T ¼ Tc is always stationary. Fluid with vertical velocity u ¼ Uy is injected into and extracted from the channel at the lower and upper walls, respectively. Note that Uy is taken as the

2.3. Substituting Eq. (19) into Eq. (18), there has gneq i ðxW ; t þ δt Þ ¼

g0 βΔTL3 ; νκ

Here, xB is a boundary node, and BE is a line perpendicular to the solid surface with xE the foot point. The temperature at xE is determined by the known boundary conditions, and TxB is obtained by the zeroth order moment of the temperature distribution function. Note that the following simulations are all implemented under the limitation of low Mach number (Ma ¼ U/cs < 0.3).

However, the non-equilibrium distribution functions at xA are usu­ ally unknown after the streaming step. To solve that problem, the antiBB rules [12,28] for the non-equilibrium distribution function is applied here, which has gneq i ðxA ; t þ δt Þ ¼

Ra ¼

where n is the normal direction with respect to the wall. Specifically, as shown in Fig. 2, Nu is obtained in the simulations according to Tx TxB �ΔT : (24) Nu ¼ E jxE xB j L

On the other hand, the non-equilibrium part can be approximated as those of xA by extrapolation, gneq i ðxW ;

ν

Pr ¼ ; κ

;

Furthermore, to evaluate the degree of convective heat transfer, the Nusselt number will be calculated, which is defined as � �

According to Eq. (7), the equilibrium part at the intersection node xW can be determined by the known interface velocity and interface tem­ perature as eq geq i ðxW Þ ¼ gi ðuW ; TÞ:

ν

where ΔT is the difference of temperature between the hot and cold boundaries. L is the characteristic length. U is the characteristic velocity, and if not given explicitly, it can be obtained as pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi U ¼ g0 βΔTL: (22)

Fortunately, the distribution functions at xC stream directly from xA

gi ðxC ; t þ δt Þ ¼ g*i ðxA ; tÞ:

UL

(20)

The thermal distribution function at xW can be eventually deter­ mined by Eqs. (17) and (20), according to Eq. (16). With the distribution functions of xC and xW at hands, the UDFs at xA can be obtained by Eq. (14) (also Eq. (13)). 7

S. Tao et al.

International Journal of Thermal Sciences 146 (2019) 106112

Fig. 12. Isotherms (left column) and streamlines (right column) in the cavity with ε/L ¼ 0.25 and D/L ¼ 0.3846: (a) and (b), θ ¼ 0; (c) and (d), θ ¼ 3π/4; (e) and (f), θ ¼ π.

characteristic velocity in the simulation, and the buoyancy force is not taken into account in this problem. In order to take advantage of the CBC, the two walls of channel both have a retraction of qδx. For this flow problem, the temperature can be obtained analytically as � y=LPrRe � e 1 TðyÞ ¼ Tc þ (25) ðTh Tc Þ: ePrRe 1

Table 3 Comparison of the surface-averaged Nu of the inner eccentric cylinder. θ

0

π/4

π/2

3π/4

π

Shu et al. [30] Lin et al. [14] Present

7.95 7.98 8.06

7.06 6.85 6.87

7.40 7.43 7.51

6.23 6.04 6.04

6.90 6.97 7.02

The L2-norm is calculated to evaluate the numerical error as

8

S. Tao et al.

International Journal of Thermal Sciences 146 (2019) 106112

which preserves well the precision of the LBM. 3.2. Natural convection in a square hollow cavity To test the convergence rate of the present scheme for the flow with thermal buoyancy, the simulation of natural convection in a square cavity is further performed. The gravity acceleration g0 points down vertically. As illustrated in Fig. 6, the cavity moves to the right side with an offset of qδx for taking advantage of the CBC. The length of the square enclosure is L. The outer boundaries are all static walls, while the left wall is hot and the right is cold, and the upper and lower walls are both set in adiabatic. The major control parameter, Ra, ranges from 1 � 103 to 1 � 106. The Pr is fixed at 0.71. The length of the cavity is resolved by 100 lattice if not stated. Fig. 7 shows the isotherms and streamlines in the square hollow cavity at different Ra. It is clearly that with the increase of Ra, the iso­ therms are distorted clockwise and the orientation changes from vertical to horizontal, indicating the variation of heat transfer from the diffusion to convection. For the feature of streamlines, only a bubble emerges in the cavity center at small Ra (Fig. 7(b)). As the disturbance of the flow gradually strengthens with larger Ra, it is torn into two (Fig. 7(d)) and three ones (Fig. 7(f)). Furthermore, the maximum y-component velocity uy,max (normalized by κ/L) at the horizontal centerline of the cavity and the maximum Nusselt number Numax along the cold wall are respectively listed in Table 1. Good agreement can be found when compared with those by Davis et al. [35]. The temperature profiles at the horizontal centerline of the enclosure are presented in Fig. 8(a). It can be observed that for both the cases of q ¼ 0.2 and 0.4, the temperature profile is able to converge, respectively. The result at the finest grid resolution, i.e., δx ¼ 1/200 is then taken as a standard solution, since there is no analytical solution available for this flow problem. The results of the error of the temperature are shown in Fig. 8(b). It is found that the present CBC is generally second-order ac­ curate in space, consistent with the simulation results of the pure diffusion case outlined in the above section.

Fig. 13. Schematic of the geometry configuration and boundary conditions for the thermal flow between two concentric cylinders.

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi, ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uX i u N h N X ðTLBM Ta Þ2 ErrorðTÞ ¼ t ðTa Þ2 i¼1

;

(26)

i¼1

where T and Ta are the temperatures obtained respectively by the pre­ sent CBC and the analytical solution derived from Eq. (25), and N is the total number of the fluid points along the channel length. Fig. 4 presents the velocity profiles at different Re and Pr, where the retraction ratio q is fixed at 0.2. The results obtained by the present thermal CBC are in good agreement with those of the analytical data by Eq. (25), indicating both the accuracy and the feasibility to the moving boundary. Furthermore, the L2-norm error of the temperature under different grid resolutions is exhibited in Fig. 5. It can be clearly observed that the convergence rate of the present thermal boundary scheme is generally second-order for different sets of Pr, Re, and q. In general, the present single-node thermal CBC is second-order accurate in space,

3.3. Natural convection in a cavity with a cylinder 3.3.1. The concentric case The third problem that considered is the natural convection in an enclosure with a circular cylinder, which has been frequently studied in the literature [14,30–32]. The simulation of the concentric case is firstly performed, as shown in Fig. 9. The temperature of the cylinder is high with T ¼ Th, and the cavity is cold with T ¼ Tc. Hence, the average

Fig. 14. Isotherms for the natural convection in the concentric annulus: (a) Ra ¼ 1.0 � 103; (b) Ra ¼ 1.0 � 104. 9

S. Tao et al.

International Journal of Thermal Sciences 146 (2019) 106112

Fig. 15. The profiles of the temperature (a) and angular velocity (b) in the annulus gap at Ra ¼ 5.0 � 104.

different Ra are presented in Fig. 10. Duo to the symmetry of geometry arrangement and the boundary condition, the flow and temperature files are all symmetrical about the vertical centerline of cavity. Furthermore, the fluid heated by the inner cylinder moves upwards, and it has a lower density near the cold wall and falls downwards. That eventually leads to the formation of two large recirculation zones in the enclosure. Partic­ ularly, the features of the figures vary greatly with the increase of the Ra. When it is small (Ra ¼ 1 � 104), the heat transfer between the cavity and the cylinder is dominated by diffusion, and the thermal convection is weak. Therefore, even though the fluid field emerges two bubbles, which has little influence on the temperature field. However, at larger Ra, the lower bubble disappears and the upper one becomes bigger. The isotherm is gradually distorted by the flow, and so the convection dominates the heat transfer in the annulus. Note that the qualitative results mentioned above are in good agreement with those in the liter­ ature [31,32]. Table 2 further exhibits the quantitative results of the surfaceaveraged Nusselt number (Nu) of the inner cylinder, which is define as Z S ∂T Nu ¼ (27) j dθ; 2π ∂n wall

Fig. 16. Schematic of a hot circular cylinder moving in a channel.

where S is half of the perimeter of the inner cylinder. The data from Shu et al. [30] and Lin et al. [14] are also included for comparison, and those results agree well with each other. 3.3.2. The eccentric case For the thermal flow in a cavity with a cylinder, the location of the cylinder has significant influence on the fluid and temperature fields. The flow can be complicated by the eccentric arrangement of the cyl­ inder, which is considered in this subsection. As plotted in Fig. 11, the boundary conditions and the grid resolution are the same with those in the concentric case. Furthermore, two parameters are introduced in the description of the eccentric, i.e., ε and θ, which denote the length be­ tween the two centers of cavity and cylinder, and the location angle between the center link and the vertical line, respectively. In the simu­ lations, Ra, D/L and ε/L are fixed at 1.0 � 106, 0.3846 and 0.25, respectively, and θ ranges from 0 to π . Fig. 12 presents the isotherms and streamlines in the enclosure under different θ. It can be observed that the thermal plum induced by hot cylinder develops well at θ ¼ 0. With the increase of θ, the cylinder approaches to the upper wall of cavity and the thermal plum is sup­ pressed. Therefore, the gap between the isotherms near the lower wall is large, indicating a small temperature gradient and then a poor perfor­ mance of heat transfer. The present results agree well with those re­ ported by Shu et al. [30] And Lin et al. [14]. Furthermore, the quantitative results of the surface-averaged Nusselt number of the inner cylinder are listed in Table 3. Good agreement can be found again when

Fig. 17. Isotherms around the cylinder (a) and the temperature profile along the channel centerline (b). xc(t) is the center of the cylinder.

temperature of the flow is obtained as T0 ¼ (Th þ Tc)/2. The outer boundaries of the enclosure are all stationary walls. The cavity length is L, which is used in the definition of Ra. The diameter of the cylinder is D, and three sets of length ratios, i.e., D/L ¼ 0.2, 0.3 and 0.4 are studied in the present work. The Ra ranges from 1 � 104 to 1 � 106, the Pr is fixed at 0.71, and the length of the enclosure is resolved by 300 lattices. The isotherms and streamlines in the cavity at D/L ¼ 0.4 and 10

S. Tao et al.

International Journal of Thermal Sciences 146 (2019) 106112

compared with those in the literature [14,30].

Acknowledgments This work was supported by the Nationa Natural Science Foundation of China (51906044), and the Doctoral Start-up Foundation of Dong­ guan University of Technology, China (Grant No. GC300502-39).

3.4. Natural convection in a concentric annulus In order to verify the present thermal CBC for more complex flows, the natural convection between two concentric cylinders are simulated. Compared with problems considered above, this case is complicated as the outer boundary of the computation domain is also curved in space. The geometry configuration and boundary conditions are illustrated in Fig. 13. The two sets of diameter and temperature for the inner and outer cylinders are (Di, Th) and (Do, Tc), respectively. In the simulations, Do/Di is set to be 2.6 and Th is bigger than Tc. The Pr is fixed at 0.7. Note that here L ¼ (Do-Di)/2 is considered as the reference length. The diameter of the outer cylinder is resolved by 200 lattices, and the angle θ is measured clockwise starting from the vertical line. The isotherms between the two concentric cylinders with Ra ¼ 1.0 � 103 and 1.0 � 104 are presented in Fig. 14(a) and (b), respectively. Those figures also illustrate the transition of the heat transfer from diffusion to convection when increasing the Ra, which are in good agreement with those in the literature [32–34]. For further verify the present scheme, simulation with Ra ¼ 5.0 � 104 is performed. Fig. 15 compares the temperature and angular velocity profiles obtained by the present thermal CBC with those in the literature [33,34]. It is observed that the results are in good agreement.

Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.ijthermalsci.2019.106112. Conflicts of interest The authors declare that there is no conflict of interest. References [1] C.K. Aidun, J.R. Clausen, Lattice-Boltzmann method for complex flows, Annu. Rev. Fluid Mech. 42 (2010) 439–472. [2] L. Fei, K.H. Luo, Cascaded lattice Boltzmann method for thermal flows on standard lattices, Int. J. Therm. Sci. 132 (2018) 368–377. [3] Y.L. He, Q. Liu, Q. Li, W.Q. Tao, Lattice Boltzmann methods for single-phase and solid-liquid phase-change heat transfer in porous media: a review, Int. J. Heat Mass Transf. 129 (2019) 160–197. [4] V.E. Ambrus¸, V. Sofonea, High-order thermal lattice Boltzmann models derived by means of Gauss quadrature in the spherical coordinate system, Phys. Rev. E 86 (1) (2012), 016708. [5] Q. Li, K.H. Luo, Q.J. Kang, Y.L. He, Q. Chen, Q. Liu, Lattice Boltzmann methods for multiphase flow and phase-change heat transfer, Prog. Energy Combust. Sci. 52 (2016) 62–105. [6] M. Jami, F. Moufekkir, A. Mezrhab, J.P. Fontaine, M.H. Bouzidi, New thermal MRT lattice Boltzmann method for simulations of convective flows, Int. J. Therm. Sci. 100 (2016) 98–107. [7] Z. Guo, C. Shu, Lattice Boltzmann Method and its Applications in Engineering, vol. 3, World Scientific, 2013. [8] L. Jahanshaloo, N.A.C. Sidik, A. Fazeli, M.P. HA, An overview of boundary implementation in lattice Boltzmann method for computational heat and mass transfer, Int. Commun. Heat Mass Transf. 78 (2016) 1–12. [9] A.J. Ladd, Lattice-Boltzmann methods for suspensions of solid particles, Mol. Phys. 113 (17–18) (2015) 2531–2537. [10] Z. Guo, C. Zheng, B. Shi, An extrapolation method for boundary conditions in lattice Boltzmann method, Phys. Fluids 14 (6) (2002) 2007–2010. [11] A. Tiwari, S.P. Vanka, A ghost fluid Lattice Boltzmann method for complex geometries, Int. J. Numer. Methods Fluids 69 (2) (2012) 481–498. [12] T. Zhang, B. Shi, Z. Guo, Z. Chai, J. Lu, General bounce-back scheme for concentration boundary condition in the lattice Boltzmann method, Phys. Rev. E 85 (1) (2012), 016701. [13] Q. Chen, X. Zhang, J. Zhang, Improved treatments for general boundary conditions in the lattice Boltzmann method for convection-diffusion and heat transfer processes, Phys. Rev. E 88 (3) (2013), 033304. [14] K.H. Lin, C.C. Liao, S.Y. Lien, C.A. Lin, Thermal lattice Boltzmann simulations of natural convection with complex geometry, Comput. Fluids 69 (2012) 35–44. [15] E. Fattahi, M. Farhadi, K. Sedighi, H. Nemati, Lattice Boltzmann simulation of natural convection heat transfer in nanofluids, Int. J. Therm. Sci. 52 (2012) 137–144. [16] R. Khazaeli, S. Mortazavi, M. Ashrafizaadeh, Application of a ghost fluid approach for a thermal lattice Boltzmann method, J. Comput. Phys. 250 (2013) 126–140. [17] I. Ginzburg, Generic boundary conditions for lattice Boltzmann models and their application to advection and anisotropic dispersion equations, Adv. Water Resour. 28 (11) (2005) 1196–1216. [18] L. Li, R. Mei, J.F. Klausner, Boundary conditions for thermal lattice Boltzmann equation method, J. Comput. Phys. 237 (2013) 366–395. [19] L. Li, R. Mei, J.F. Klausner, Lattice Boltzmann models for the convection-diffusion equation: D2Q5 vs D2Q9, Int. J. Heat Mass Transf. 108 (2017) 41–62. [20] J. Huang, Z. Hu, W.A. Yong, Second-order curved boundary treatments of the lattice Boltzmann method for convection–diffusion equations, J. Comput. Phys. 310 (2016) 26–44. [21] J.W.S. McCullough, C.R. Leonardi, B.D. Jones, S.M. Aminossadati, J.R. Williams, Investigation of Local and Non-local Lattice Boltzmann Models for Transient Heat Transfer between Non-stationary, Disparate Media. Computers & Mathematics with Applications, 2018. [22] C.S. Dai, M. Li, H.Y. Lei, S.X. Wang, Numerical simulation of natural convection between hot and cold microtubes in a cylinder enclosure, Int. J. Therm. Sci. 95 (2015) 115–122. [23] M. Torabi, A. Keyhani, G.P. Peterson, A comprehensive investigation of natural convection inside a partially differentially heated cavity with a thin fin using twoset lattice Boltzmann distribution functions, Int. J. Heat Mass Transf. 115 (2017) 264–277. [24] Z. Li, M. Yang, Y. Zhang, Lattice Boltzmann method simulation of 3-D natural convection with double MRT model, Int. J. Heat Mass Transf. 94 (2016) 222–238.

3.5. A hot cylinder moving in a channel The present method is further validated in the simulation of a hot cylinder moving in a channel. As shown in Fig. 16, the length and width of the channel respectively are L ¼ 12.5D and H ¼ 2D, where D is the diameter of the circular cylinder. The cylinder has a high temperature Th ¼ 1.1, which migrates at a constant velocity Uc ¼ 5 � 10 3 along the centerline of the channel. The upper and lower boundaries of the channel are both cold walls with temperature Tc ¼ 1.0, and the left and right are periodic boundaries. In the simulations, the cylinder is resolved by 80 lattices. The above mentioned parameters are the same with those in the literature [13,39]. Note that in this test case, the state of some points can be changed from fluid to solid and vice versa. The averaged extrapolation method is used to obtain the distribution function at the newly released fluid nodes [13,40]. Fig. 17 shows the isotherms around the cylinder and the temperature profile along the centerline of the channel. The temperature distribution is asymmetrical about the cylinder, indicating the convection effect on the thermal diffusion. The temperature profile by the present method agrees well with those by Chen et al. [13] and Hu et al. [39]. 4. Conclusions A thermal boundary condition is proposed for the lattice Boltzmann simulation of flows with heat transfer in this paper. The key point in the present scheme is to determine the non-equilibrium distribution func­ tion, which is consistent with the isothermal counterpart. However, the anti-bounce back rule should be used for the thermal problems. With only a single fluid node involved, the present method has the property of local computation. That feature can contribute to the overall high per­ formance computing, giving the present scheme great potential for application in related fields such as flow and heat transfer in dense porous media. The second-order accuracy of the present scheme is first validated in the thermal Couette flow with wall injection, and the nat­ ural convection in a square cavity. Subsequently, simulations of three numerical cases, including the natural convections in a cavity with a cylinder and in a concentric annulus, and a hot cylinder moving in a channel are performed to test the feasibility and accuracy of the present scheme. The results are in good agreement with the data in the literature. 11

S. Tao et al.

International Journal of Thermal Sciences 146 (2019) 106112 [33] J. Wu, Y. Cheng, L.A. Miller, An iterative source correction based immersed boundary-lattice Boltzmann method for thermal flow simulations, Int. J. Heat Mass Transf. 115 (2017) 450–460. [34] T.H. Kuehn, R.J. Goldstein, An experimental and theoretical study of natural convection in the annulus between horizontal concentric cylinders, J. Fluid Mech. 74 (4) (1976) 695–719. [35] G. de Vahl Davis, Natural convection of air in a square cavity: a bench mark numerical solution, Int. J. Numer. Methods Fluids 3 (3) (1983) 249–264. [36] P. Lallemand, L.S. Luo, Lattice Boltzmann method for moving boundaries, J. Comput. Phys. 184 (2) (2003) 406–421. [37] X. Yin, J. Zhang, An improved bounce-back scheme for complex boundary conditions in lattice Boltzmann method, J. Comput. Phys. 231 (11) (2012) 4295–4303. [38] S. Tao, Q. He, B. Chen, X. Yang, S. Huang, One-point second-order curved boundary condition for lattice Boltzmann simulation of suspended particles, Comput. Math. Appl. 76 (7) (2018) 1593–1607. [39] J. Hu, Z. Guo, Numerical study on mass transfer from a composite particle settling in a vertical channel, Int. J. Heat Mass Transf. 117 (2018) 132–142. [40] S. Tao, J. Hu, Z. Guo, An investigation on momentum exchange methods and refilling algorithms for lattice Boltzmann simulation of particulate flows, Comput. Fluids 133 (2016) 1–14.

[25] 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. [26] L.S. Luo, W. Liao, X. Chen, Y. Peng, W. Zhang, Numerics of the lattice Boltzmann method: effects of collision models on the lattice Boltzmann simulations, Phys. Rev. E 83 (5) (2011), 056710. [27] B. Chun, A.J.C. Ladd, Interpolated boundary condition for lattice Boltzmann simulations of flows in narrow gaps, Phys. Rev. 75 (6) (2007), 066705. [28] Q. Zou, X. He, On pressure and velocity boundary conditions for the lattice Boltzmann BGK model, Phys. Fluids 9 (6) (1997) 1591–1598. [29] Z. Wang, Y. Liu, J. Zhang, N. Dang, Study of laminar natural convection in a vertical annulus with inner wall covered by a porous layer by using lattice Boltzmann method, Int. J. Therm. Sci. 135 (2019) 386–397. [30] C. Shu, H. Xue, Y.D. Zhu, Numerical study of natural convection in an eccentric annulus between a square outer cylinder and a circular inner cylinder using DQ method, Int. J. Heat Mass Transf. 44 (17) (2001) 3321–3333. [31] R. Roslan, H. Saleh, I. Hashim, A.S. Bataineh, Natural convection in an enclosure containing a sinusoidally heated cylindrical source, Int. J. Heat Mass Transf. 70 (2014) 119–127. [32] 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.

12