Cascaded lattice Boltzmann method for thermal flows on standard lattices

In this paper, a thermal cascaded lattice Boltzmann method (TCLBM) is developed in combination with the double-distribution-function (DDF) approach on the standard D2Q9 lattice. A density distribution function relaxed by the cascaded scheme based on central moments is employed to solve the flow field, and a total energy distribution function relaxed by the BGK scheme is used to solve the temperature field. The two distribution functions are coupled naturally to provide a new TCLBM. In this method, the viscous heat dissipation and compression work are taken into account, the Prandtl number and specific-heat ratio are adjustable, and the external force is considered directly without the Boussinesq assumption. The TCLBM is validated by numerical experiments of the thermal Couette flow, low-Mach number shock tube problem, Rayleigh-Bénard convection, and natural convection in a square cavity with a large temperature difference. The simulation results agree well with the analytical solutions and/or results given by previous researchers.


Introduction
The lattice Boltzmann method (LBM), based on the kinetic theory, has achieved remarkable success as an alternative method to conventional computational fluid dynamics (CFD) for thermal flow and heat transfer applications during the past three decades [1][2][3][4][5][6][7][8][9]. Different from solving the discretized Navier-Stokes (N-S) equations in traditional CFD methods, the LBM solves a discrete kinetic equation at the mesoscopic scale, designed to reproduce the N-S equations in the macroscopic limit. The main advantages for LBM over traditional CFD include [10,11]: convenience to deal with complex boundary, easiness of programming, high parallel efficiency, and natural incorporation of micro and meso-scale physics.
The basic algorithm realization of LBM is collision-streaming or streaming-collision, although other time and space evolution schemes can also be used. To be specific, at each time step the collision is first locally executed and followed by streaming the post-collision distributions to their neighbors, or just exchanging the above procedure [12]. Based on this algorithm, various collision operators can be adopted, such as the single-relaxation-time (SRT) or BGK operator [13], tworelaxation-time (TRT) operator [14,15], multiple-relaxation-time (MRT) operator [16,17], and entropic operator [18][19][20]. Compared with these extensively used operators, cascaded or central moment operator, first proposed by Geier et al. [21] in 2006, is more recent. The collision in the cascaded Lattice Boltzmann method (CLBM) is performed by relaxing central moments to their local equilibrium values separately, which is different from MRT LBM where the raw moments are relaxed. As mentioned in Ref. [21], central moments can be expressed as polynomials of raw moments of the same order and below. When a raw moment is relaxed (in MRT), all central moments at the same or higher orders will be changed. This "cross-talk" is a source of instability and can be removed in CLBM. By choosing the relaxation parameters properly, CLBM can be adopted to simulate very high Reynolds number flows using coarse grids without adopting any turbulence models or entropic stability [21]. Recently, Lycett-Brown and Luo [22] extended the CLBM to multiphase flow using the interaction potential method [23] with the EDM force scheme [24]. Compared with the LBGK method, the proposed model provided significant improvement in reducing spurious velocities, and increasing the stability range for the Reynolds number and liquid to gas density ratio. They further extended the model to three dimensions and achieved high Weber number, high Reynolds number and high density ratio simultaneously in binary droplet collision simulations [25,26]. More recently, based on a generalized multiple-relaxation-time (GMRT) framework, we proposed a consistent method to incorporate a force field into CLBM and clarified the relation between CLBM and MRT LBM [27]. Although CLBM has obtained success in high Reynolds number single-phase flows and multiphase flows, its applications are so far limited to incompressible flows. Recently, we showed that for incompressible thermal flows, CLBM can improve the numerical stability significantly compared with the BGK model [9]. The purpose of the present study is to extend CLBM to Low-Mach compressible thermal flows. Generally, there are three feasible ways to construct thermal LBMs. The first one, multispeed approach [28,29], is a straightforward extension of athermal to thermal LBMs, in which more discrete velocities are adopted to match higher-order moment constraints of the density distribution function for recovering the energy equation. In the second one, a density distribution function is still used to simulate the velocity field, while other methods, such as finite difference or finite volume [30,31], are adopted for the temperature field. The doubledistribution-function (DDF) [1,2] approach is the third one, where two different distribution functions are adopted to solve flow and temperature fields, respectively. In DDF-based thermal LBMs, the compression work and heat dissipation can be simply included, and the specific-heat ratio and Prandtl number are adjustable. On the whole, the DDF approach keeps the intrinsic features and simple structures of the standard LBM, and more comparisons and discussions among the three methods can be found in Refs. [2,4,32]. In the history, the first DDF thermal model was proposed by He et al. [1] by using an internalenergy-distribution-function-based DDF approach. Guo et al. [2] then presented another DDF thermal model using a total energy distribution function to solve the energy equation, which is simpler than He and coworkers' model. In Guo and co-authors' model, the local temperature in equilibrium density and energy distribution functions is replaced by the reference temperature, thus it is a decoupling model and is limited to Boussinesq flows. In 2012, Li et al. [4] developed a coupling DDF thermal model which can simulate more general thermal flows, and the model was extended to three-dimensions by Feng et al. [33] recently. Inspired by these works, we construct a thermal cascaded lattice Boltzmann method (TCLBM) in the present work based on the DDF approach. In the TCLBM, a density distribution function is relaxed using the cascaded scheme, a total energy distribution function is relaxed using the SRT scheme, and the external force is considered directly without the Boussinesq assumption.
The rest of the paper is structured as follows: Section 2 briefly introduces the cascaded LBM. Section 3 presents a method to incorporate the force field into cascaded LBM. In Section 4, we extend the athermal CLBM to TCLBM. Numerical experiments are carried out for several benchmark problems to validate the proposed model in Section 5. Finally, conclusions of this work are made in Section 6.

Cascaded LBM
In this paper, the D2Q9 lattice [13] is adopted, and the discrete velocities here δ x and δ t are the lattice spacing and time step, and = = = c δ δ 1 x t is used in this work. For the derivation of CLBM, we follow Lycett-Brown and Luo [22] and begin with the velocity moments of the discrete distribution function (DF) f a , and then f a and f a eq can be formulated as functions of the corresponding moments and equilibrium moments.
The raw moments are defined as Recombining the second-order moments, the trace of the pressure tensor, the normal stress difference and the off diagonal element of the pressure tensor are given by , . 20 02 20 02 11 (2) According to the definition above, we get the raw moment representation of populations: (3b) (3c) (3d) (3e) [ ] , 5 21 12 22 (3f) [ ] , 6 21 12 22 (3g) [ ] , 7 21 12 22 (3h) 8 21 12 22 (3i) It should be noted that other variables can also be expressed using their moments in this form similarly. Central moments are defined in a reference frame shifted by the local velocity, The transformation between the raw moments and central moments can be expressed using the binomial theorem as given by Lycett-Brown and Luo [22]. To construct a CLBM, we follow the assumption adopted in Ref. [34], by setting the discrete equilibrium central moments equal to the corresponding continuous values, where f eq is the local Maxwell-Boltzmann distribution for athermal fluid at temperature T 0 in continuous particle velocity space ξ ξ ( , ) x y , and the lattice sound speed = c R T s 0 is set to be 1/ 3 in this work. Substituting Eq. (6) into Eq. (5), we can calculate the second order and above central moments, and write them using the combination as done in raw moments: eq eq eq eq eq eq 21 12 The implementation of CLBM is also composed of collision step and streaming step. For the collision step, central moments are relaxed to their equilibrium values, separately: where w 1 and w 2 are dependent on the shear and bulk viscosities, re- ), and the parameters for the third-and fourth-central moments (w 3 and w 4 ) are freely tunable to control the stability. The post-collision raw moments can then be recovered according to the binomial theorem, 3 .
x y x y x y (9f) Then we get the post-collision distribution using Eq. (3), and the streaming step is as usual,

Incorporating forcing terms into cascaded LBM
To include the force effect on the flow field, we define f a changes due to this force field by a source term S a . To match the overall accuracy in LBM, one way to add the source term in CLBM is to employ the second-order trapezoidal rule along the characteristic line, To remove the implicitness in Eq. (11), the transformation method in Ref. [34] is adopted, He et al. [1] proposed that the presence of the force field = F F F ( , ) x y changes the continuous distribution function as follows: We then follow the assumption in Ref. [34] that the discrete central moments of S a are equal to the continuous central moments of f Δ F : Substituting Eq. (13) into the integral, we get the nonzero central moments of S a , where a x and a y are horizontal and vertical components of the acceleration. As suggested by Premnath and Banerjee [34], we remove thirdorder moments in Eq. (15) henceforward for convenience. Using the binomial theorem once again, we yield analytical raw moments of S a , Thus, the analytical expressions of S a can be written in the same form as Eq. (3). From the definition in Eq. (12), the conserved raw moments of the , respectively. The corresponding non-conserved raw and central moments can then be obtained straightforwardly, With Eqs. (7) and (15), the non-conserved equilibrium central moments will remain the same as the ones before transformed, thus the collision step for the central moments will not be affected in Eq. (8). According to the relationship between raw moments mentioned above, the postcollision raw moments are slightly different from Eq. (9), x y x y y x ) , x y x x y y * * 2 2 (18c) 3 .
x y x y x y x y x x y y y x ) to get f a * , the streaming step is then given as: The hydrodynamic variables are then obtained as:

Coupling DDF cascaded LBM for thermal flows
For thermal flows, the temperature T is now a function of space and time, not a constant value T 0 . The equilibrium distribution function f eq is given by and then the reference temperature T 0 in Eq. (7) should be replaced by local temperature T. Inspired by the total-energy-based DDF models [2,4], we adopt a density distribution function relaxed by the cascaded scheme to solve the flow field, together with a total energy distribution function using the BGK relaxation scheme to simulate the temperature field, and the two fields are coupled through the ideal gas equation of state (EOS, = p ρRT ). The discrete total energy distribution function has the kinetic equation [2], To recover the compressible N-S equations, discrete raw moments of f a eq should be consistent with the continuous raw moments of f eq from the zeroth-to third-order. As mentioned in Sec. 2, two of the third-order raw moments are not independent due to the lack of symmetry in D2Q9 lattice, The last term at the RHS is a deviation from the continuous moments for f eq , where = θ T T / 0 , and = δ 1 ijkl . This means that the diagonal elements ( = δ 1 ijkl ) for the third-order velocity moments deviate from the needed relationship. As pointed out by Prasianakis and Karlin [35], the deviation can be removed only by adding a correction term C a into the evolution equation for standard lattices. According to Li and co-workers' work [4], the raw moments for C a should satisfy, The other raw moments can be chosen as: In the above, all the third-order velocity terms have been neglected because of the low Mach number limit. Then all the central moments for C a are zero except: In the simulation, the derivative terms can be evaluated using a secondorder central difference. Then the analytical expressions of C a can be written in the same form as Eq. (3).
By using a second-order trapezoidal rule, the evolution equation can be written analogously: Due to the non-zero second-order central moments ( ∼ E c and ∼ N c ) for C a , the equilibrium central moments for the transformed distribution f a should be: The cascaded relaxation for central moments is in the same form as Eq. (8), while the dynamic viscosity μ and bulk viscosity μ B are: Because the conserved raw moments for C a are zero, the calculation for post-collision raw moments is the same as Eq. (18) In the same manner, we use a transformed total energy distribution function: Then the time-discrete form of Eq. (22) is, The macroscopic temperature is obtained by, If the compression work and viscous heat dissipation are negligible, we can adopt an analogous internal energy (temperature) distribution function [4] to solve the temperature field. This concludes the development of a new thermal cascaded lattice Boltzmann method (TCLBM).

Numerical experiments
In this section, a series of numerical experiments are conducted to verify the developed model. Unless otherwise specified, is adopted in simulations.

Thermal Couette flow
To check the capability of describing viscous heat dissipation by the present TCLBM, two-dimensional thermal Couette flow is simulated. We consider the viscous fluid between two infinite parallel plates, in which the upper one is moving at speed U with temperature T 0 and the lower adiabatic plate is fixed. With the assumption that = Pr μc λ / p is constant and = μ μ T T ( / ) ( / ) 0 0 , there is an analytical solution [36], 2)/ is the specific-heat ratio, = Ma U γRT / 0 is the Mach number, and H is the distance between the two plates.
In our simulations, we set x y is employed. For the top and bottom walls, the non-equilibrium bounce-back method [37] and nonequilibrium extrapolation method [2,38] are adopted for velocity and temperature boundary conditions, respectively, while the periodic boundary condition is imposed along the x direction. The upper wall temperature and the reference temperature are set to unity, with a reference dynamic viscosity = μ 0.35 0 . The relaxation parameter w 1 is a field variable related to the local dynamic viscosity as given in Eq. (29).  Table 1. The relative error is defined as 0 . As presented in Table 1, the relative errors are less than 1% in all the cases.

Low-Mach shock tube problem
To check the present model's ability of simulating Low-Mach number compressible flow, a shock tube problem is studied in this section. The construction of this problem is that a long tube containing the same gas is separated by a barrier in the middle into two parts with different pressures, densities and temperatures. At the moment of removing the barrier, a complex flow is set up. The initial conditions for our simulations are, 0 are the reference density, temperature, pressure and velocity, respectively, and L 0 is the length of the tube.
In the simulation, a × 1000 5 lattice is used, the periodic boundary condition is imposed along the y direction, while EDFs are used in = x 0 and = x 1000. The specific heat ratio γ and Prandtl number Pr are set to 1.4 and 0.71, with = w 1.89 1 . Simulation results are compared with the analytical ones in Fig. 2. The four plots present dimensionless density, pressure, horizontal velocity, and temperature profiles, respectively, at time = t δ 520 t . It can be observed that numerical results are in good agreement with the theoretical ones.

Rayleigh-Bénard convection
To check the ability of simulating thermal flow with external force by the present TCLBM, the numerical experiment of the Rayleigh-Bénard convective flow is conducted in this section. The Rayleigh-Bénard instability is one of the classical thermal instability phenomena, in which the fluid is enclosed between two parallel stationary walls, cold at the top and hot at the bottom, and experiences the gravity force. Linear stability theory has proven that convection develops most readily when the wave number is at the critical value 3.117 [40], which is approximately corresponding to length-width ratio 2:1 in the flow domain.
Since the present model is a coupling model, we can implement the force by means of central moments with = a 0 x and = − a g y (g is the gravity acceleration), without using the Boussinesq assumption. We conduct the experiment in the weakly compressible regime, with a 6% temperature difference of the reference temperature = T 1.0 0 . To delete the heat dissipation term, we adopt the internal energy distribution function to simulate temperature field [4,41]. The non-equilibrium extrapolation scheme [2,38] is used to treat the upper and lower wall boundaries for both velocities and temperatures, whereas the periodic boundary scheme is applied along the horizontal direction. The Prandtl number corresponds to air, = Pr 0.71. Then the flow is characterized by   should be set to be an appropriate value, for example 0.08 in our simulation, to keep the flow in the low-Mach number regime. And β is the thermal expansion coefficient, which is the reciprocal of reference temperature for the ideal gas considered here.
For this kind of instability phenomenon, the driven force by the density variations induced by the temperature variations will balance with the viscosity force at the critical Rayleigh number Ra c , while if the Rayleigh number is increased above the threshold, the driving force will dominate and convection will start. First, we use a × = × N N 200 100 x y grid to calculate the critical Rayleigh number. We initialize the temperature field with a linear distribution in the y direction and give a small perturbation for density around the reference density = ρ 1.0 0 . It is noted that the total kinetic energy will keep increasing/decreasing lineally after the initial unsteady period around the critical Rayleigh number. For that, the total kinetic energy increment e Δ every 10000 time steps in the domain is measured, where e Δ is measured by the slope of the total kinetic energy change with time, not at a certain time step. The critical Rayleigh number extrapolated is = Ra 1707.07 c (see Fig. 3   Bénard convection at = Ra 5000, 10000 and 50000 are shown in Fig. 4. 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. To quantify this, the Nusselt number in the system is calculated, where the square bracket represents the average over the whole system and k is the thermal conductivity of the fluid. The obtained values of Nusselt number Nu n at various Rayleigh numbers are compared with the reference data in Table 2, and plotted in Fig. 5. The simulation results are in good agreement with those of Ref. [40] in a wide range of Ra as given in Table 2. During the small Rayleigh number range ( < Ra 5000), convection is suppressed so that the Nusselt number decreases rapidly to 1.0 at = Ra Ra c , where the empirical formula loses efficacy. At very high Rayleigh numbers, the numerical results slightly underestimate the heat transfer, while this trend was also observed in other LBM studies [1,33,35].

Natural convection in a square cavity
In this subsection, the present TCLBM is employed to simulate the natural convection in a square cavity. In this case, the large temperature difference condition is considered, in which the left and right walls are maintained at = T 960 h K and = T 240 l K respectively, and the horizontal sidewalls are adiabatic. The dynamic viscosity is defined by the Sutherland's Law: where = T 273 * K, = S 110.5 K. The μ * is the dynamic viscosity at T * , which is determined from the reference viscosity μ r at the reference temperature = T 600 r K. In the simulation. The Prandtl number is taken as 0.71, μ r is set to be 0.14, 0.1, 0.06, and 0.055 for = Ra 10 3 , 10 4 , 10 5 , and 10 6 , and the grid sizes are 100 × 100, 150 × 150, 250 × 250 and 650 × 650, respectively. The isotherms and streamlines are presented in Figs. 6 and 7. The streamlines are drawn through the numerical integration of the steady compressible stream function y y x ), where the partial derivative is discretized by the second-order central difference. The asymmetry feature of the convection with a large temperature difference is clearly shown in Figs. 6 and 7: For = Ra 10 3 and 10 4 , the center of the primary vortex shifts to the lower right side of the cavity center, and the two small vortexes at higher Rayleigh numbers are not similar in size and shape, which agree well with the benchmark solutions given in Ref. [42].
To  [4,42]. The results predicted by the present TCLBM are compared with the previous benchmark solutions [42] and LB results [4] in Table 3. From the table we can see that the present results are in good agreement with the data in previous works.

Conclusions
In this paper, we developed a thermal cascaded lattice Boltzmann method (TCLBM) for low-Mach compressible thermal flows on standard lattices. Considering the proven numerical performances of DDF-based thermal LBMs, we constructed the TCLBM in this framework. Firstly, the reference temperature in equilibrium central moments for the density DF was replaced by the local temperature. Secondly, a correction term was introduced similarly by means of central moments to remove the derivation of two diagonal elements for the third-order raw moments. Then a total energy EDF was introduced according to the required raw moments. Finally, by relaxing the density DF and total energy DF using the cascaded and BGK schemes respectively, a DDF thermal CLBM was constructed, where the density DF solves the flow field and the total energy DF solves the temperature field and they are coupled naturally by EOS for the ideal gas.
To verify the proposed model, a thermal Couette flow was simulated first, and the simulation results agreed well with the analytical solutions in different cases, which demonstrated the present TCLBM can include viscous heat dissipation with different Prandtl number and specific-heat ratio. The TCLBM's ability of simulating low-Mach compressible flows was then verified by simulating a low-Mach shock tube problem. The numerical results for the Rayleigh-Bénard and square cavity convections confirmed that the model can simulate thermal problems with force field without invoking the Boussinesq assumption.
In summary, the present TCLBM retains the simplicity and numerical efficiency of the DDF-LB method while the numerical stability of CLBM is preserved. In this method, the viscous heat dissipation and compression work are taken into account, the Prandtl number and specific-heat ratio are adjustable, and the external force is considered directly without the Boussinesq assumption. Finally, the present method can be extended to three dimensions (3D) readily based on the 3D cascaded LBM [21,25,26].

Acknowledgments
Support from the MOST National Key Research and Development Programme (Project No. 2016YFB0600805) and the Center for Combustion Energy at Tsinghua University is gratefully acknowledged.