Cascaded lattice Boltzmann method for incompressible thermal flows with heat sources and general thermal boundary conditions

Cascaded or central-moment-based lattice Boltzmann method (CLBM) is a relatively recent development in the LBM community, which has better numerical stability and naturally achieves better Galilean invariance for a specified lattice compared with the classical single-relation-time (SRT) LBM. Recently, CLBM has been extended to simulate thermal flows based on the double-distribution-function (DDF) approach [L. Fei \textit{et al.}, Int. J. Heat Mass Transfer 120, 624 (2018)]. In this work, CLBM is further extended to simulate thermal flows involving complex thermal boundary conditions and/or a heat source. Particularly, a discrete source term in the central-moment space is proposed to include a heat source, and a general bounce-back scheme is employed to implement thermal boundary conditions. The numerical results for several canonical problems are in good agreement with the analytical solutions and/or numerical results in the literature, which verifies the present CLBM implementation for thermal flows.


Introduction
In the last three decades or so, the lattice Boltzmann method (LBM), which is a mesoscopic numerical method based on the kinetic theory, has been developed to be a widely used numerical method for solving various fluid flows and heat transfer problems [1,2,3,4,5,6,7].In the LBM, a discretized Boltzmann equation, based on a specific discrete velocity set and designed to reproduce the Navier-Stokes (N-S) equations in the macroscopic limit, is solved for the distribution functions (DFs).Generally, the mesoscopic nature of LBM allows its natural incorporation of microscopic and/or mesoscopic physical phenomena, while the highly efficient algorithm makes it affordable computationally [8,9].
In the extensively used algorithm for LBM, the numerical process can be split into two steps [8,10,9]: the "collision" step and the "streaming" step.In the collision step, the single-relaxation-time (SRT) or BGK 10 scheme [3] is the most widely used collision operator.In the BGK-LBM, all the distribution functions are relaxed to their local equilibrium states at an identical rate, where the relaxation rate is related to the kinematic viscosity.The BGK-LBM is quite simple to implement and can simulate low Reynolds number flows, but it may have numerical instability at high Reynolds number or low-viscosity flows, as well as inaccuracy of implementing the boundary conditions [11,12,13,14,15].To overcome these difficulties, the multiple-relaxation-time (MRT) collision operator was proposed in the literature [11,12].In the MRT-LBM, the DF is transformed into a raw moment space, where different raw moments of the DF can be relaxed at different relaxation rates to the local equilibrium raw moments.Compared with the BGK-LBM, the MRT-LBM can enhance numerical stability by carefully separating the time scales among the kinetic modes [12], as well as improve the numerical accuracy for non-slip boundary conditions by choosing a specified relaxation rate for the energy flux [13,14,15].However, Geier et al. argued that the MRT-LBM may also encounter instability for high Reynolds number flows due to the insufficient degree of Galilean invariance and the "cross-talk" effect induced by relaxing the raw moments [16].By relaxing central moments of the DF in the co-moving frame, a cascaded or central-moment-based operator was proposed in 2006 [16].In the cascaded LBM, also known as CLBM, the "cross-talk" effect in the MRT-LBM is eliminated naturally, and a higher degree of Galilean invariance for a specified lattice can be preserved readily by matching the higher order central moments of the continuous Maxwell-Boltzmann distribution.By setting the relaxation rates for high-order central moments to be 1, CLBM has been applied to simulate high Reynolds number (Re = 1400000) turbulent flow using coarse grids without resorting to any turbulence models [16].Recently, CLBM has been extended to simulate multiphase flows coupled with the pseudo-potential model [17] by Lycett-Brown and Luo [18].Compared with the BGK-LBM for multiphase flows, the proposed multiphase CLBM reduces the spurious currents near the phase interface significantly [18], and achieves higher stability range for the Reynolds number [19].As is known, the basic pseudo-potential model has some drawbacks, such as thermodynamic inconsistency, large spurious currents, and suffers from the problem of tuning the surface tension independently of the density ratio [9].More recently, Li et al. proposed an approach of achieving thermodynamic consistency via tuning the mechanical stability condition [20,21], and analyzed the effects of the equation of state on the thermodynamic consistency [22].Inspired by the methods in [20,21,22], an improved forcing scheme in pseudo-potential model was proposed in [23].By coupling the improved forcing scheme with the cascaded operator, Lycett-Brown and Luo achieved very high parameters in the simulation of binary droplet collision [24].
More recently, CLBM was first extended to simulate thermal flows by the present authors [25], where a thermal cascaded lattice Boltzmann method (TCLBM) was proposed based on the double-distributionfunction (DDF) approach.In our TCLBM, the CLBM is used to simulate the flow field and another total energy BGK-LBM is used for the temperature field, where the two fields are coupled by equation of state for the ideal gas.The proposed TCLBM has been proved to be able to simulate low-Mach compressible thermal flows with commendable stability and accuracy.For incompressible thermal flows without viscous dissipation and pressure work, another CLBM has been constructed on a simpler lattice (D2Q5) to solve the passive-scalar temperature field [26].Compared with the D2Q5 MRT-LBM for the temperature equation, the proposed D2Q5 CLBM is shown to be better Galilean invariant.Thus a higher characteristic velocity can be adopted for convection heat transfer problems, which decreases the computational load significantly.
Although CLBM has been applied to several thermal problems [25,26], less attention has been paid to two important factors: temperature field with a heat source and non-isothermal boundary conditions.In this work, we will present the implementation of a heat source and a general bounce-back scheme for the thermal boundary conditions.
The rest of the paper is structured as follows: In Section 2, a brief introduction for the DDF-based CLBM for incompressible thermal flows is given, followed by the implementation of a heat source and general bounce-back scheme for thermal boundary conditions.Numerical experiments are carried out for several benchmark problems to validate the employed method in Section 3. Finally, a brief summary is given in Section 4.

Numerical method
The macroscopic governing equations for incompressible thermal flows can be written as: where u, p, ρ 0 , T , ν and α are the velocity, pressure, reference density, temperature, kinematic viscosity, and thermal diffusivity, respectively.The Boussinesq approximation is employed in this work, thus the force field is defined as, where the gravitational acceleration vector g points to the negative direction of y-axis, β is the thermal expansion coefficient, T 0 is the reference temperature, and F v is an external body force.

CLBM for the flow field
In the present work, the D2Q9 discrete velocity model [3] is used to simulate two-dimensional problems.As usual, the lattice spacing ∆x, time step ∆t and lattice speed c = ∆x/∆t are set to be 1.The discrete velocities e i = [|e ix , |e iy ] are defined by where i = 0...8, |• denotes the column vector, and the superscript ⊤ indicates transposition.
For the cascaded collision operator, the collision step is carried out in the central-moment space.The raw moments and central moments of the discrete distribution functions (DFs) f i are defined as: and the equilibrium values k eq mn and keq mn are defined analogously by replacing f i with the discrete equilibrium distribution functions (EDFs) f eq i .In this work, a simplified raw-moment set is adopted [26], and so do the central moments Γi .Specifically, the raw moments can be given from f i through a transformation matrix M by |Γ i = M |f i , and the central moments shifted from raw moments can be performed through a shift matrix N by Γi = N |Γ i .The formulations for M and N can be easily obtained according to the raw-moments set [27].In the present study, M and N are expressed as [26], The post-collision central moments can be obtained by relaxing each of them to their local equilibrium states independently, Γ * where the block-diagonal relation matrix is given by, where ρ is the fluid density, and c s = 1/3 is the lattice sound speed.Consistently, the forcing source terms in central moments space are given by [27], It may be noted that the method of incorporating a force field into the CLBM is the most recently proposed consistent forcing scheme [27] and it shows great advantages over the previous forcing schemes in CLBM.
In the streaming step, the post-collision discrete DFs f * i in space x and time t stream to their neighbors in the next time step as usual, where the post-collision discrete DFs are determined by |f Using the Chapman-Enskog analysis, the incompressible N-S equaltions in Eqs.(1) can be reproduced in the low-Mach number limit [27,29].The hydrodynamics variables are obtained by, It should be noted that the incompressible approximation [30] is employed in the present work.Thus the dynamic variable density ρ can be divided into the reference density ρ 0 and a small density fluctuation δρ.3), {e i = [|e ix , |e iy ] |i = 0, 1, ...4 }) can be used to construct the CLBM for the temperature field [26].Similarly, the raw moments and central moments of the temperature distribution functions g i are defined by [26], In the D2Q5 lattice, the following five raw moments are adopted [26], and so do the central moments ΓT i .Analogously, the raw moments and central moments can be calculated through a transformation matrix M T and a shift matrix N T , respectively [26], Explicitly, the transformation matrix M T is expressed as [26], and the shift matrix N T is given by, The collision in the central-moment space can be written as, where is the diagonal relaxation matrix.The thermal diffusivity is related to the relaxation parameter by α = (1/λ 1 − 0.5)c 2 T ∆t.The equilibrium values of the central moments are given by, ΓT,eq where c 2 T is the sound speed in the D2Q5 lattice.The post-collision temperature distribution functions g * i can be obtained by The streaming step for g * i is also as usual, The temperature T is computed as, Through the Chapman-Enskog analysis, the convection-diffusion equation for the temperature field can be recovered in the macroscopic limit.

Heat source and boundary conditions
The DDF-based CLBM introduced above has been proved to be able to simulate several incompressible thermal flows with isothermal boundary condition.However, it can hardly simulate convective heat transfer problems with a heat source.Inspired by the previous method to include the heat source in the BGK and MRT LBM [31,32], here we present a CLBM for the temperature equation with a generalized heat source term.Similar to the consistent forcing scheme in CLBM, a heat source Q can be incorporated into Eq.( 18) by means of central moments, where R i correspond to the central moments of the heat source, Analogously, the calculation of temperature is modified, To implement thermal boundary conditions, a general bounce-back scheme is adopted in this work.After the collision step, the post-collision temperature distribution functions are obtained by Eq. ( 20).In the streaming step, the distribution functions entering from "outside" of the boundary g i (x f , t + ∆t) are determined by, where e i = −e i , and T w is the temperature at the boundary.For the general thermal boundary conditions, b 1 ∂T w /∂n+b 2 T w = b 3 , the boundary temperature T w can be solved using a finite-difference scheme.Different from the method in [33], a second-order finite-difference scheme is adopted for the temperature gradient, where T 1 and T 2 are temperatures at the first and second layer nodes neighboring the boundary, and n is the boundary normal vector.The boundary temperature can be calculated as, After obtaining T w , the unknown distribution functions g i (x f , t + ∆t) can be calculated using Eq. ( 26).

Numerical experiments
In this section, several benchmark problems are conducted to verify our implementation of the heat source and boundary conditions.In the present CLBM for the temperature field, the value of c T can be independent of c s , and is set to be c T = 2/5 in this work.Unless otherwise specified, the half-way bounceback boundary scheme is used for both velocity and temperature boundary conditions, while s 3 in Eq. ( 8) is chosen according to the non-slip rule s 3 = (16 − 8s ν )/(8 − s ν ) [27].

Time-independent diffusion problem
The first tested problem is a time-independent diffusion problem, which can be described by the following simplified equation and boundary conditions, where T 0 and T L are the temperatures at the bottom and the top of a straight channel.The heat source is Q = 2α∆T /L 2 , with ∆T = (T L − T 0 ), and the exact solution is, Due to the simple flow configuration, only 6 nodes are used to cover the channel width (L = 6∆x).The simulation results are compared with the analytical solution in Fig. 1.Two cases with α = [1/10, 1/3] are considered, where the boundary conditions are T 0 = 0 and T L = 1, respectively.The corresponding relaxation rates are chosen as: (1) λ 1 = 4/3 and λ 2 = 3/4 for the first case; (2) λ 1 = 3/4 and λ 2 = 4/3 for the second case.It is seen that the simulation results are in very good agreement with the analytical solution.
As analyzed by Cui et al. [32], when the relaxation rate λ 2 is specified as λ 2 = 12(λ 1 − 2)/(λ 1 − 12), the numerical slip in the D2Q5 MRT can be eliminated.To check its applicability in the present D2Q5 CLBM, a series of simulations are carried out with λ 2 changing from 0.2 to 1.8.As shown in Fig. 2, the global relative error E 2 , defined as E 2 = (T − T a ) 2 / T a 2 , reaches the minimum values at λ 2 = 3/4 and λ 2 = 4/3 for α = 1/10 and α = 1/3, respectively.Thus the non-slip rule in the D2Q5 MRT is also suitable for the present D2Q5 CLBM, which further verifies our previous analysis that the MRT-LBM and CLBM can be put into a unified general framework [27].

Viscous dissipation in Poiseuille flow
To validate the implementation of a spatially variable heat source, viscous dissipation in Poiseuille flow is simulated.The flow is driven by a constant body force along x direction, F = [F x , 0], while the walls are at constant temperature T w .The viscous dissipation is considered by adding a heat source, Q = ν(∂u/∂y) 2 , in Eq. (1c).By using the non-slip rule for the velocity field, a very accurate velocity profile can be provided by the D2Q9 CLBM in Section.2.1.The analytical temperature field is [31], where h is the half-width of the channel.The simulation result is compared with the analytical solution in Fig. 3, where the dimensionless temperature is defined as, It is clearly shown that the simulation result agrees well with the analytical solution.The global relative errors at different grid sizes are shown in Fig. 4. A very good linear fit is seen in the simulation results, and the slop is 2.05.It indicates that the implementation of the boundary conditions and heat source for the present problem has second-order accuracy in space.

Natural convection in a square cavity
The natural convection driven by the buoyancy force in a square cavity is simulated to validate the implementation of complex thermal boundary conditions.This problem has been widely examined in the literatures [25,34,35,36].The left and right walls of a square cavity are at constant temperature T h = 1 and T l = 0, respectively, while the top and the bottom walls are adiabatic.The problem is characterized by the Prandtl number P r = ν/a and Rayleigh number Ra = gβ(T h − T l )H 3 /(νa), where H is the cavity hight.In the present paper, P r is set to be 0.71, and the characteristic velocity U = gβ(T h − T l )H is set to be 0. and 6, respectively.Qualitatively, all the characteristics in both temperature and flow fields agree well with the results in previous studies [34,35,36].To be more quantitative, data of the present work are listed in Table 1, compared with those reported in previous studies [36,35].The following quantities are compared: the maximum horizontal velocity component u max at x = H/2 and its location y max , the maximum vertical 120 velocity component v max at y = H/2 and its location x max , and the average Nusselt number N u along the cold wall.There is an excellent agreement between the present results and the benchmark solutions in the previous studies [35,36].

Conclusions
In this work, we extend previous DDF-based thermal CLBM to simulate more general incompressible thermal flows with heat sources and thermal boundary conditions.To include a heat source in the temperature equation, a discrete source term R i is added to the collision step in central-moment space.To deal with thermal boundary conditions, the general bounce-back boundary scheme in MRT-LBM is modified and adopted in the present D2Q5 CLBM.Through numerical simulations of several benchmark cases, very good accuracy of the proposed implementation for the heat source and boundary conditions are confirmed.In with s + = (s b + s ν )/2 and s − = (s b − s ν )/2 [28, 26].The kinematic ν and bulk viscosities ν b are related to the relaxation parameters by ν = (1/s ν − 0.5)/3 and ν b = (1/s b − 0.5)/3, respectively.The equilibrium central moments Γeq i are defined equal to the continuous central moments of the Maxwellian-Boltzmann distribution in continuous velocity space,

70 2 . 2 .
CLBM for the temperature field Due to the simplicity of convection-diffusion equation, a D2Q5 discrete velocity model (the five discrete velocity set is defined in Eq. (

Figure 1 :
Figure 1: Comparison of temperature profiles predicted by the D2Q5 CLBM simulation and the analytical solution.

Figure 3 :
Figure 3: Comparison of the dimensionless temperature profiles predicted by the D2Q5 CLBM simulation and the analytical solution.

1 . 5 ∆Figure 4 :
Figure 4: Global relative errors E 2 change with grid sizes for viscous dissipation problem in Poiseuille flow.