Modeling incompressible thermal flows using a central-moment-based lattice Boltzmann method

In this paper, a central-moment-based lattice Boltzmann (CLB) method for incompressible thermal flows is proposed. In the method, the incompressible Navier-Stokes equations and the convection-diffusion equation for the temperature field are sloved separately by two different CLB equations. Through the Chapman-Enskog analysis, the macroscopic governing equations for incompressible thermal flows can be reproduced. For the flow field, the tedious implementation for CLB method is simplified by using the shift matrix with a simplified central-moment set, and the consistent forcing scheme is adopted to incorporate forcing effects. Compared with several D2Q5 multiple-relaxation-time (MRT) lattice Boltzmann methods for the temperature equation, the proposed method is shown to be better Galilean invariant through measuring the thermal diffusivities on a moving reference frame. Thus a higher Mach number can be used for convection flows, which decreases the computational load significantly. Numerical simulations for several typical problems confirm the accuracy, efficiency, and stability of the present method. The grid convergence tests indicate that the proposed CLB method for incompressible thermal flows is of second-order accuracy in space.


Introduction
As a mesoscopic numerical method based on the kinetic theory, the lattice Boltzmann method (LBM) [1,2,3] has obtained remarkable success in the applications to fluid flows and heat transfer problems during the past three decades [4,5,6,7,8,9]. The LBM solves a discreted Boltzmann equation, designed to recover the Navier-Stokes (N-S) equations in the macroscopic limit. The highly efficient and easy algorithm of LBM makes it at afforadable computatinal cost, while the meso-scale nature allows its natural incorporation of micro and/or meso-scale physics [7].
In the standard collision-streaming algorithm for LBM, the simplest collision operator is the singlerelaxation-time (SRT) or BGK operator, in which all the distirbution functions are relaxed to their local equilibrium values at a common rate [1]. However, the BGK-LBM may meet troubles of inaccuracy in im- 10 plementing the boundary conditions [10,11,12], as well as numerical instability at high Reynolds numebr or low-viscosity flows [13]. To overcome these difficulties, other collision operators, such as multiple-relaxationtime (MRT) operator [14,13], entropic operator [15,16,17] and two-relaxation-time (TRT) operator [18,19] have been developed. Recently, a cascaded or central-moment-based operator was proposed by Geier et al. [20]. The collision in the central-moment-based Lattice Boltzmann method (CLBM) is carried out by relaxing central moments of the discrete distribution functions separately, rather than raw moments as in the MRT-LBM. By matching the higher order central moments of the continuous Maxwell-Boltzmann distribution naturally, CLBM achieves a higher order of Galilean invariance [20]. Meanwhile, as mentioned by Geier et al. [20], central moments can be expressed as polynomials of raw moments at the same order and below, which means relaxing raw moments (in MRT) affects the independent relaxation for central moments 20 at higher orders. The "cross-talk" is a source of numerical instability, and can be removed through relaxing central moments in CLBM. By setting high-oder central moments to their equilibrim values, CLBM has been used to simulate turbulence flow at Re = 1400000 using coarse grids without resorting to any turbulence models or entropic stabilization [20]. More recently, CLBM has been extended to multiphase flows based on the interaction potential method [21] by Lycett-Brown and Luo [22] . Compared with the BGK-LBM for multiphase flows, the proposed multiphase CLBM enables significant improvment in reducing spurious currents near the phase interface [22], and achieving higher stability range for the Reynolds number [23]. They further extended the model with an improved forcing scheme [24], and made a breakthrough for large density ratio multiphase flow with high Reynolds and Weber numbers simultaneously [25]. Although CLBM has gained success in high Reynolds number single-phase flows and multiphase flows, 30 its applications are still mainly limited to isothermal flows [26]. The motivation of the present work is to extend CLBM to incompressible thermal flows. Up to now, the double-distribution function (DDF) approach has been extensively used for constructing thermal LBMs [8,27,28,29,30,31,32,33]. In DDFbased incompressible thermal LBMs, a density or pressure distribution function is used to solve the velocity field, with another distribution for temperature field, where the two fields are usually coupled through the Boussinesq approximation [29,30]. Due to the simplicity of the convection diffusion equation, a simpler lattice is often used for the temperature field to decrease the computational load and save the memory storage [29,34]. The MRT operator is also widely used for the temperature field to achieve better numerical stability and more accuarte boundary conditions [35,36,37]. Inspired by these studies, we try to extend the CLBM to incompressible thermal flows in the present study based on the DDF approach. Specifically, 40 a density distribution function and a temperature distribution function are used to simulate the velocity field and the temperature field, respectively, and both of them are relaxed using the central-moment-based operater. The rest of the paper is structured as follows: In Section 2, the 2D CLBM for incompressible thermal flows is presented in detail. Numerical experiments are carried out for several benchmark problems to validate the proposed method in Section 3. Finally, concluding remarks are given in Section 4.

CLBM for incompressible thermal flows
In this section, details for the construction of the CLBM for incompressible thermal flows are given. The macroscopic govering equations for the flow fields are where u, p, ρ 0 , F and ν are the velocity, pressure, reference density, force field and kinematic viscosity, respectively. The convection-diffusion equation for a scalar variable φ with diffusion coefficient D can be written as ∂φ ∂t For the incompressible thermal flows considered in this study, the scalar variable and diffusion coefficient are specified as temperature T and thermal diffusivity α, respectively. To include the effect of temperature field on the flow field, the Boussinesq assumption is used and the force field is defined by 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 50
The two-dimensional (2D) problems are considered in this study, and the D2Q9 lattice [1] is used for the flow field. The lattice speed c = ∆x/∆t = 1 is adopted, in which ∆x and ∆t are the lattice spacing and time step. The discrete velocities e i = [|e ix , |e iy ] are defined by where i = 0...8, |· indicates the colunm vector, and the superscript ⊤ indicates the transposition.
To construct the central-moment-based collision operator, raw moments and central moments for the discreted distribution functions (DFs) f i are introduced, and the equilibrium values k eq mn andk eq mn are defined analogously by replacing f i with the discrete equilibrium distribution functions (EDFs) f eq i . In the literature, many researchers [20,22,23,25,26,38,39] using the recombined raw moments k 20 + k 02 , k 20 − k 02 , to treat the trace of the pressure tensor and the normal stress difference independently. To simplify the tedious implementation of using the recombined raw-moment set, a simplified method was proposed in [40], where the raw moments k 20 , k 02 were used and some modifications were made in the relaxation matrix. In this work, the simplified raw-moment set is adopted, and so do the recombined central momentsΓ i . To be more specific, 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 realized through a shfit matrix N by Γ i = N |Γ i . Similar to the expressions in [38], M and N are written as, By relaxing each central moment to its equilibrium counterpart independently, the post-collision central moments are given by where C i are the forcing source terms in central moments space, and the block-diagonal relation matrix is given by, with s + = (s b + s ν )/2 and s − = (s b − s ν )/2. The equalibrium central moments of the f eq i are set equal to the continuous central moments of the Maxwellian-Boltzmann distribution in continuous velocity space, where ρ is the fluid density, and c s = 1/3 is the lattice sound speed. The corresponding EDF is in fact a generalized local equilibrium [40,39]. Consistently, the forcing source terms in central moments space are given by [38], In the streaming step, the post-collision discrete DFs 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 Eq. (1) can be reproduced in the low-Mach number limit [39]. The hydrodynamics variables are obtained by, The kinematic and bulk viscosities are related to the relaxation parameters by ν = (1/s ν − 0.5)c 2 s ∆t and ξ = (1/s b − 0.5)c 2 s ∆t, respectively.

CLBM for the temperature field
A new D2Q5 (the five discrete velocity set is defined in Eq. (4), {e i = [|e ix , |e iy ] |i = 0, 1, ...4 }) CLBM is proposed to solve the convection-diffusion equation for the temperature field. Similarly, the raw moments and central moments of the temperature distribution functions g i can be defined by In the D2Q5 lattice, the following five raw moments are considered, and so do the recombined 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 Explicitly, the transformation matrix M T is expressed as and the shift matrix N T is given by The collision in central moments can also be written as where The equibrium values of the central moments are given by where c 2 T,s 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 are also as usual The temperature T is computed as As shown in Appendix A, the convection-diffusion Eq. (2) can be recovered by the D2Q5 CLBM presented in this subsection, and the thermal diffusivity is related to the relaxation parameter by α = (1/λ 1 − 0.5)c 2 T,s ∆t.

Numerical experiments
In this section, we conduct several benchmark cases to verify the effectiveness and accuracy of the proposed CLBM for incompressible thermal flows. Unless otherwise specified, the lattice sound speed for the D2Q5 CLBM is set to c sT = 1/2, the tunable relaxation parameters for high-order central moments are 60 set to 1.0, and the the non-equilibrium bounce-back method [41] and non-equilibrium extrapolation method [42] are adopted for velocity and temperature boundary conditions in simulations, respectively.

The decay of a temperature wave
Firstly, the decay of a temperature wave on a moving frame is considered. The problem is specified by the following velocity and temperature fields: (24a) where φ = 2π/L, A represents the vertical reference velocity component, and B is the initial amplitude of the temperature wave. Periodic boundary conditions are used along the x and y axes. In the simulations, T 0 = 1.0, L = 100, and the initial amplitude is set to B = 0.01. The velocity field is given, thus only the D2Q5 CLBM is adopted to solve the temperature field. Another two D2Q5 MRT LBMs in [36] and [37] are also used for comparison, and they are denoted by MRT-LBM1 and MRT-LBM2, respectively.  present D2Q5 CLBM, the measured thermal diffusivity is independent of the reference velocity (or M a) and always agrees with the given value. For the other two D2Q5 MRT LBMs, the measured diffusivities decrease with the increase of the reference velocity. To be specific, the relative errors at M a = 0.3 are around 12% and 8% for MRT-LBM1 and MRT-LBM2, respectively.

Normal plate velocity problem with a temperature difference
The normal plate velocity problem is a fully developed channel flow, where the upper plate moves with a uniform velocity u 0 , and a uniform normal flow with velocity v 0 is injected through the bottom plate and withdrawn from the upper plate. The analytical solution of the flow is given by [43], where the Reynolds number is based on the width of the channel, L, and defined by Re = v 0 L/ν. For the present study, a temperature difference ∆T = T H − T L is considered, where T H and T L are the temperatures at upper hot plate and bottom cold plate, respectively. The steady temperature profile satisfies [29], where P r = ν/α is the Prandtl number.
In the simulations, T H = 1.0 and T L = 0 are used, and the viscosity is set to ν = 0.1. Firstly, we set the Reynolds number Re = 10, u 0 = v 0 = Reν/L with L = 50. Periodic boundary conditions are adopted at the inlet and outlet of channel, and the length of the channel is covered by 5 grids to save the computational cost. Three simulation cases with P r = [0.1, 1, 10] are conducted to verify the numerical performace of the present method at a wide range of the Prandtl number. The tunable parameter λ 2 is set according to the zero-numerical-slip condition [37], λ 2 = 12(λ 1 − 2)/(λ 1 − 12). The residual error E R < 1 × 10 −9 is used as the convergent criterion, the defination of E R can be seen in [38]. As shown in Fig. 3, the numerical results for the non-dimensional temperature T * = (T − T 0 )/∆T agree very well with the analytical solutions.
Then we set the Prandtl number corresponding to air, P r = 0.71, with different values of the Reynolds number, Re = [10,20,30]. The width of the channel covered by a series of grid nodes, L = [30,60,90,120,150], are considered to validate the convergence rate in space. The relative errors of temperature and velocity are calculated according to the following definations, The relationships between grid size and relative errors of the present method are plotted in Fig. 4, and 90 the slopes of each fitting lines are very close to 2.0. This demonstrates that the method proposed is of second-order convergence rate in space.

Rayleigh-Bénard convection
In this section, the Rayleigh-Bénard convective flow is conducted to check the ability of simulating incompressible thermal flows with an external force field. In the 2D Rayleigh-Bénard convective flow, the fluid is enclosed between two parallel stationary walls, with high temperature T H at the bottom and low temperature T L at the top, and experiences the gravity. The gravity field is incorporated by Eq.
where ρ 0 = 1.0, while other points are initialized as ρ 0 . The characteristic time of the system can be expressed by t c ∼ H/u c = Hc s /M a. It is known that the iterations needed for the convergence are proportional to t c . Firstly, we change 1/M a from 0.1 to 0.31 with 100 a 0.03 interval, and calculte the needed time steps until convergence. From Fig. 5, we certainly confirm the linear relation between the needed time steps and 1/M a. Throughout the variation range of M a, the relative changes of the Nusselt number (defined in Eq. 29) are only 0.11%, 0.10% and 0.08% for the Rayleigh numbers at 2500, 3000, and 5000, respectively. Based on this, a higher Mach number can be used in the present method to reduce the computational cost, but not affects the numerical accuracy. In the following simulations, the Mach number is set to M a = 0.3.
According to the linear stability theory, the driven force by the density variations induced by the temperature variations will be balanced by the viscous force when Rayleigh number is lower than a critical value Ra c , while if the Rayleigh number is increased above the threshold, the driving force will dominate and convection will be induced. To determine the critical Rayleigh number, we measure the evolution of the maximun 110 vertical velocity in the system V max for a series of Rayleigh numbers, Ra = [1702, 1704, 1706, 1708, 1710]. It can be seen in Fig 6 that V max will keep increasing/decreasing approximately linearly in the early period, depending on Ra. The critical Rayleigh number is determined by solving zero value for the grow rate of V max with the least-squre method [35]. Compared with the exact value in the linear stability theory, Ra c = 1707.76, the value based on our method, 1706.82, is satisfying.
Flows at different Rayleigh numbers are then simulated. Figs. 7 display the normalized temperature (T − T 0 ) /∆T at Ra = 2000, 10000 and 50000. When the Rayleigh number increases, we can see two clear trends in the figures: the mixing of the hot and cold fluids is enhanced, and the temperature gradients near the bottom and top walls are increased, both of which mean the convective heat transfer is enhanced in the domain. In the meantime, as shown in Figs .8, the vortex is gradually distorted with the increase of the Rayleigh number, which also means the enhancement of convection.To quantify this, the Nusselt number in the system is calculated, where the square bracket represents the average over the whole system. Nusselt numbers obtained at various Rayleigh numbers are compared with the reference data in Table 1. The simulation results are in good agreement with the analytical values in Ref . [44]. On the whole, the relative errors for the present method are smaller than the method in [45], while only 30 nodes are used in the present method rather than 50 nodes as in [45].

120
In the end, it is interesting to find that the proposed method is stable for Ra reaching up to 10 9 with only 220 nodes in the vertical direction as shown in Fig .9, which confirms the commendable stability of the method. However, the study of high-Rayleigh-number thermal convection is beyond the scope of the current paper.

Conclusions
In this paper, we have developed a central-moment-based lattice Boltzmann (CLB) model for simulation of incompressible thermal flows. Combined with the D2Q9 CLB equation for the flow field, another D2Q5 CLB equation is designed to reproduce the temperature equation. Through the Chapman-Enskog analysis, the macroscopic govering equations for incompressible thermal flows can be recovered. By matching the higher order central moments of continuous Maxwell-Boltzmann distribution naturally, the proposed thermal 130 model achieves better Galilean invariance compared with some existing thermal LBMs. Thus a higher Mach number can be adopted for convection flows in the present model, which reduces the computational cost significantly. Through the simulation of several canonical problems, the accuray, efficiency, stability and the second-order convergence rate in space for the present model are verified. The method developed retains the simplicty and numerical efficiency of the standard LBM. In principle, the method can be extended to three-dimensional readily. Besides, the model developed can also be applied to other convection-diffusion problems directly.  [44,45].