A Truly Second-order and Unconditionally Stable Thermal Lattice Boltzmann Method

An unconditionally stable thermal lattice Boltzmann method (USTLBM) is proposed in this paper for simulating incompressible thermal flows. In USTLBM, solutions to the macroscopic governing equations that are recovered from lattice Boltzmann equation (LBE) through Chapman–Enskog (C-E) expansion analysis are resolved in a predictor–corrector scheme and reconstructed within lattice Boltzmann framework. The development of USTLBM is inspired by the recently proposed simplified thermal lattice Boltzmann method (STLBM). Comparing with STLBM which can only achieve the first-order of accuracy in time, the present USTLBM ensures the second-order of accuracy both in space and in time. Meanwhile, all merits of STLBM are maintained by USTLBM. Specifically, USTLBM directly updates macroscopic variables rather than distribution functions, which greatly saves virtual memories and facilitates implementation of physical boundary conditions. Through von Neumann stability analysis, it can be theoretically proven that USTLBM is unconditionally stable. It is also shown in numerical tests that, comparing to STLBM, lower numerical error can be expected in USTLBM at the same mesh resolution. Four typical numerical examples are presented to demonstrate the robustness of USTLBM and its flexibility on non-uniform and body-fitted meshes.


Introduction
Incompressible thermal flows are frequently encountered in engineering applications [1].Therefore, much research has been conducted on this topic.Numerical simulation is an emerging and popular approach in studying this topic [2][3][4].Among various numerical methods, lattice Boltzmann method (LBM) is attracting more and more interest because of its simplicity and explicitness [5][6][7][8][9].However, with its broader applications, some drawbacks of standard LBM are revealed and reported, which include high cost in virtual memories, inconvenient implementation of physical boundary conditions, lattice uniformity and poor numerical stability in high Reynolds/Rayleigh numbers [10].
To overcome or alleviate these drawbacks, a simplified lattice Boltzmann method (SLBM) was recently developed and extended to simulations of thermal flows [11,12].Simplified thermal lattice Boltzmann method (STLBM) is derived from reconstructing solutions to the macroscopic governing equations recovered from lattice Boltzmann equation and resolved in a predictor-corrector scheme.Final formulations of STLBM only involve the equilibrium and the non-equilibrium distribution functions.The equilibrium distribution functions are evaluated from the macroscopic variables, while the non-equilibrium distribution function is approximated from the difference of two equilibrium distribution functions at different locations and time levels.Consequently, the macroscopic variables are directly updated in STLBM, which brings various merits to this method.On the one hand, virtual memories are reduced dramatically because the distribution functions no longer need to be saved.On the other hand, physical boundary conditions can be directly implemented.However, one flaw in STLBM constrains its robustness in numerical perspective, which is its first-order of accuracy in time in the corrector step.Under the motivation of eliminating this flaw, we present a truly second-order unconditionally stable thermal lattice Boltzmann method (USTLBM) in this paper.Assisted by Chapman-Enskog expansion analysis and Taylor series expansion formulas, higher order of accuracy can be achieved in approximating the non-equilibrium terms.Meanwhile, slight shifting of half time step can be obtained in the time level, which enables us to construct a typical central difference scheme that ensures the second-order of accuracy in time.
While eliminating the flaws, USTLBM maintains all merits of STLBM.The final formulations of USTLBM are very similar to those of STLBM, with slight changes made in the corrector step.Therefore, the simplicity and explicitness are well maintained.At the same time, USTLBM also directly updates the macroscopic variables.Thus, virtual memories are greatly lowered and physical boundary conditions can be directly implemented.Most importantly, through von Neumann stability analysis, USTLBM can be theoretically proven unconditionally stable, which indicates that it inherits the good numerical stability of STLBM and has great advantages over the standard LBM.
The remaining parts of this paper are laid out as follows.Section 0 gives brief introduction to the conventional thermal lattice Boltzmann method and the recently developed simplified thermal lattice Boltzmann method (STLBM).Derivations of the unconditionally stable thermal lattice Boltzmann method (USTLBM), as well as its boundary treatments and stability analysis, are discussed in detail in Section 3. Four typical numerical examples are presented in Section 4 for comprehensive evaluation of the robustness and flexibility of USTLBM.Conclusions are finally summarized in Section 5.

Thermal Lattice Boltzmann Method (TLBM) and Chapman-Enskog (C-E) Expansion Analysis
In incompressible thermal fluids, effects of temperature and velocity are coupled.Within LBM framework, one popular approach to interpret such coupling effect is to introduce internal energy distribution function [13].This approach is named as the thermal lattice Boltzmann method (TLBM) and the evolutions of the distribution functions can be expressed as g α (r + e α δ t , t + δ t ) = g α (r, t) + g eq α (r, t) − g α (r, t) τ c , α = 0, 1, ... , N where f α and g α represent the density distribution function and the internal energy distribution function along the α direction, respectively, with the superscript "eq" denoting their equilibrium state; τ υ and τ c are single relaxation parameters related to the kinematic viscosity and the thermal diffusivity, respectively; e α and δ t represent the lattice velocity and the time step applied in the model; and M and N are the number of the lattice velocities used in the density and the internal energy distribution functions, respectively.The equilibrium density distribution function is given by Appl.Sci.2017, 7, 277 3 of 21 while the equilibrium internal energy distribution functions in two and three dimensions are, respectively, expressed as 2 evaluated from the conservation laws of mass, momentum and energy: Symmetric lattice velocity directions are adopted in TLBM.The lattice velocity model is named as DaQb, where a and b stand for the dimension and the number of lattice velocity directions, respectively.In the present work, typical D2Q9 model and D3Q15 model (see Figure 1) are applied in two-dimensional and three-dimensional simulations, respectively.
e u e u r u e u e u (5) where c = δx/δt is set to be 1 in this work; and δx is the lattice spacing.
Although TLBM directly updates distribution functions, macroscopic properties can still be evaluated from the conservation laws of mass, momentum and energy: , , Symmetric lattice velocity directions are adopted in TLBM.The lattice velocity model is named as DaQb, where a and b stand for the dimension and the number of lattice velocity directions, respectively.In the present work, typical D2Q9 model and D3Q15 model (see Figure 1) are applied in two-dimensional and three-dimensional simulations, respectively.Chapman-Enskog (C-E) expansion analysis is a multi-scale expansion method which is capable to link TLBM to macroscopic governing equations.It can be proven that [14,15], to the second-order of accuracy both in time and in space, Equations (1) and ( 2) are equivalent to the following macroscopic expressions: where Π and Q are the momentum flux tensor and the thermal flux vector in the following form where the non-equilibrium distribution functions have the expressions of The single relaxation parameters are related to kinematic viscosity and thermal diffusivity by

Simplified Thermal Lattice Boltzmann Method (STLBM)
The simplified thermal lattice Boltzmann method (STLBM) is a novel method that is recently developed within LBM framework [12].Its derivation starts from the macroscopic governing equations recovered from TLBM, i.e., Equations ( 8)- (10).Assisted by fractional step technique, these equations can be resolved in a predictor-corrector scheme Predictor step: ∂ρe ∂t Corrector step: According to the Boussinesq approximation [16], the external force term can be regarded as a buoyance force related to the temperature profile, i.e., where g is the gravitational acceleration; β is the thermal expansion coefficient; and T m is the average temperature.By using lattice properties and the relationships given by C-E analysis, the above Equations ( 17)-( 21) can be reconstructed within the LBM framework to obtain the formulations of simplified thermal lattice Boltzmann method (STLBM) Predictor step: Corrector step: ρ(r, t)e(r, t) = ρ * e * + 1 − 1 Inspired by the treatment in lattice Boltzmann flux solver (LBFS), the non-equilibrium terms involved in Equations ( 28) and (29) are approximated by the difference between two equilibrium distribution functions at different locations and time levels, i.e., According to Shu et al. [17] and Wang et al. [18], the equilibrium distribution functions at time level of t is evaluated from the intermediate macroscopic variables obtained in the predictor step.
Compared to the standard TLBM, simplified thermal lattice Boltzmann method (STLBM) possesses many advantages.Although it is developed within LBM framework, STLBM directly updates macroscopic variables instead of the distribution functions.In this way, much fewer virtual memories are required in the computation, and physical boundary conditions can be directly and flexibly implemented without transforming into conditions for the distribution functions.Moreover, by performing von Neumann stability analysis, STLBM can be theoretically proven to be unconditionally stable, which makes it suitable for computations at high Reynolds/Rayleigh numbers.
Despite of the above merits, a serious flaw in STLBM theory is its accuracy.By performing Taylor-series expansion analysis and conservation laws, it can be proven that Equations ( 24)- (26) in the predictor step can recover corresponding macroscopic Equations ( 17)- (19) with the second-order of accuracy both in time and in space.However, in the corrector step, it can be shown that when recovering macroscopic Equations ( 21) and (22), Equations ( 28) and ( 29) maintain the second-order of accuracy in space, but present a Euler scheme in time marching, which is only in the first-order of accuracy (see [12] for detail).Since, in LBM framework, the time interval and the mesh spacing are tied up, low order of accuracy in time may have intrinsic effect on overall accuracy of numerical solution.

Modifications on STLBM
The initial formulations in the corrector step are shown in Equations ( 27)- (29).Apparently, Equation (27) recovers Equation (20) exactly.In the present work, Equations ( 28) and ( 29) are, respectively, replaced by the following formulas To present theoretical analysis on the accuracy, we start from non-equilibrium distribution functions involved in the above formulations.Performing Taylor-series expansion analysis on Equation (30) with third-order of truncated error gives Similarly, we have Equations ( 34) and ( 35) indicate that, with desired order of accuracy, the approximated non-equilibrium distribution functions in Equations ( 30) and (31) are at space/time level of (r − 0.5e α δ t , t − 0.5δ t ).Such shifting in time and space favors our construction of a central-difference scheme both in space and in time, with which the second-order of accuracy can be guaranteed.Maintaining its location in time level, we further expand Equation (34) with respect to the space/time levels of (r, t − 0.5δ t ), which gives Substituting Equations ( 36) and (37) into Equation ( 32) leads to Apparently, the first two terms on right hand side of Equation ( 38) represent the approximation of the space derivative term in Equation ( 21) with the second-order of accuracy in space.Meanwhile, the approximated space derivatives are defined at the time level of (t − 0.5δ t ), which represents a typical central difference scheme that ensures the second-order of accuracy in time marching.Through similar processes, it can also be proven that Equation (33) recovers Equation ( 22) with the second-order of accuracy both in time and in space.

Unconditionally Stable Thermal Lattice Boltzmann Method (USTLBM)
With the modifications made on previous simplified thermal lattice Boltzmann method (STLBM), formulations of the present unconditionally stable thermal lattice Boltzmann method (USTLBM) are shown as follows: Predictor step: Corrector step: where the non-equilibrium terms are evaluated by To ensure the second-order of accuracy in time, the external force term F E should be defined on (r, t − 0.5δ t ).According to Boussinesq approximation shown in Equation (23), this term is approximated by Since T(r,t) can be solved firstly and then substituted into Equation (47) in the same time step, the equations are decoupled and the explicitness of the scheme is ensured.
As discussed above, USTLBM is an improved version of STLBM and is a truly second-order method.Nevertheless, the basic idea of USTLBM is similar to that of STLBM, which is to reconstruct solutions to macroscopic governing equations within LBM framework.Such characteristics grant USTLBM advantages from both STLBM and standard LBM.On one hand, USTLBM maintains the simplicity and explicitness, and ensures the second-order of accuracy in time and space as standard LBM does.On the other hand, USTLBM directly updates macroscopic variables rather than distribution functions, which greatly lowers virtual memories and facilitates implementation of physical boundary conditions.Moreover, as its name indicates, USTLBM can be theoretically proven to be unconditionally stable, which will be discussed in detail in following sections.

Boundary Conditions
In its computational process, USTLBM calls for boundary values of both macroscopic variables and non-equilibrium distribution functions.The physical boundary condition, which is related to macroscopic variables, can be directly implemented.Specifically, Dirichlet boundary condition can be imposed by artificially assigning boundary values of macroscopic variables on boundary meshes; and Neumann boundary condition can be fulfilled by adopting information on several inner layers of meshes to achieve desired order of accuracy.
Boundary values of non-equilibrium distribution functions are required in the corrector step of USTLBM, as shown in Equations ( 43) and (44).It has been proven in previous studies that the non-equilibrium extrapolation scheme can ensure the second order of accuracy in space [17].In this treatment, the boundary values of non-equilibrium distribution functions are set equal to those on the first inner layer of mesh points.Such boundary treatment can be flexibly applied in uniform and non-uniform meshes, and will be used in the present method in determining boundary values of non-equilibrium terms.

Stability Analysis
Theoretical analysis for the stability of lattice Boltzmann method (LBM) was firstly conducted by Sterling and Chen [18].It was found that standard LBM suffers from numerical instability when the single relaxation parameter τ is close to 0.5, which indicates that standard LBM is not suitable to problems at high Reynolds numbers.The recently developed simplified thermal lattice Boltzmann method (STLBM), on the other hand, performs very well in numerical stability.Through von Neumann stability analysis, it can be proven that STLBM is unconditionally stable [12].The USTLBM proposed in this paper can be regarded as an improved version of STLBM that ensures higher order of accuracy in time.In this section, we use von Neumann stability analysis to study its stability performance.For the ease of analysis, we take the two-dimensional formulations as the example here.Analysis on the three-dimensional formulations can be carried out through similar processes, and same conclusions can be obtained.
To begin with, the following macroscopic variables are defined: The vector of these variables is then written as The equilibrium distribution functions Equations ( 3) and ( 4) are functions of the macroscopic variables, which can be written as In von-Neumann stability analysis, the variables are approximated by the Fourier serial expansions [2].Consider a random component of the variables encountered here, which can be written as In addition, note that the following isotropic properties of the lattice tensors are generally applicable within LBM framework: In the corrector step, it is shown in Equation (42) that no correction is performed on the density.Therefore, we can simply derive that dy Next, we still take the variable y 2 to illustrate the analyzing process in the corrector step.According to Equations (40), (43), and (45), we have where the subscript "−α" represents the space level of (r + e α δ t ).Similar to the procedures conducted in the corrector step, here we linearize the formula by differentiating it, which gives dS = ∂S ∂y (64) Thus, the overall stability analysis of USTLBM can be written in the matrix form of Apparently, all eigenvalues in the characteristic matrix are 1.Therefore, it is proven that USTLBM proposed in this paper is theoretically unconditionally stable.
Through the above von Neumann stability analysis, it turns out that USTLBM maintains the merit of STLBM in numerical stability; therefore, both of them are advantageous over the standard LBM.Meanwhile, considering its advantage in temporal accuracy, USTLBM is more competitive than STLBM in the perspective of numerical accuracy.

Computational Procedure
The computational sequence of the present unconditionally stable thermal lattice Boltzmann method (USTLBM) can be summarized as follows: (1) Specify the streaming distance δ x (δ x = δ t ).Determine the single relaxation parameters τ υ and τ c according to Equations ( 15) and ( 16).

Numerical Tests
In this part, we present four numerical examples to demonstrate the robustness and the flexibility of USTLBM in modeling incompressible thermal flows.Numerical accuracy and stability are evaluated through typical numerical experiments to validate the theoretical conclusions made in previous sections (i.e., USTLBM ensures the second-order of accuracy in space and is unconditionally stable).Combining with Lagrangian interpolation algorithms, the present method can be flexibly extended to non-uniform meshes and body-fitted meshes to achieve higher efficiency and to study problems with curved boundaries.

2D Natural Convection in a Square Cavity
Natural convection in a square cavity is a classic benchmark test to validate numerical methods for thermal flows.In this section, we use the example in two-dimensional scenario to assess accuracy, convergence, and stability performance of USTLBM.Physical configuration of this problem is shown in Figure 2. Key non-dimensional parameters that affect the flow pattern are the Rayleigh number and the Prandtl number defined as where ∆T is the temperature difference between the left hot wall and the right cold wall; and V c is the characteristic velocity defined as V c = gβL • ∆T.In our simulations, Pr is set to be 0.71, and V c = 1 to fulfill the low Mach number assumption.Various Nusselt numbers are defined to quantify the heat transfer rate.Specifically, the volume averaged Nusselt number and the average Nusselt number on the hot wall are Initially, we use the case with Ra = 10 4 to validate spatial accuracy of USTLBM.Since analytical solution to this problem is unavailable, the benchmark result of the volume averaged Nusselt number obtained by the mesh size of 301 × 301 is used as the reference data Nu a .Various mesh sizes of 11 × 11, 21 × 21, 41 × 41, 61 × 61 and 81 × 81 are adopted and corresponding relative errors are calculated through the following formula The relative errors obtained at different mesh sizes are depicted in Figure 3.As can be seen, the linearly fitted line of the data reveals a slope of 1.956 in log scale, which indicates that USTLBM can ensure the second-order of accuracy in space.The result obtained by STLBM is also included in this figure.The fitted line of the data given by STLBM also has a slope that is very near to 2.0, which validates that STLBM can ensure the second-order of accuracy in space.However, it should be noted that the order of accuracy and the numerical error (e.g., the relative error to the reference data in this example) are two different concepts.The numerical error is usually written as Error = O(h n ) = C•h n , where h is the mesh spacing, n is the order of accuracy, and C is a constant.Clearly, the order of accuracy is the slope of straight line of log (error) vs. log(h).Even if the order of accuracy n is the same, the relative numerical error of one method can still be smaller than that of another method at the same mesh size (i.e., C value of two methods may not be the same).Detailed comparisons show that, at the same mesh spacing, relative error of the results given by USTLBM is lower than that obtained by STLBM.Therefore, it is reasonable to expect more accurate results at the same mesh resolution when applying USTLBM.L g u = 0, v = 0, ∂T/∂y = 0 The relative errors obtained at different mesh sizes are depicted in Figure 3.As can be seen, the linearly fitted line of the data reveals a slope of 1.956 in log scale, which indicates that USTLBM can ensure the second-order of accuracy in space.The result obtained by STLBM is also included in this figure .The fitted line of the data given by STLBM also has a slope that is very near to 2.0, which validates that STLBM can ensure the second-order of accuracy in space.However, it should be noted that the order of accuracy and the numerical error (e.g., the relative error to the reference data in this example) are two different concepts.The numerical error is usually written as Error = O(h n ) = C•h n , where h is the mesh spacing, n is the order of accuracy, and C is a constant.Clearly, the order of accuracy is the slope of straight line of log (error) vs. log(h).Even if the order of accuracy n is the same, the relative numerical error of one method can still be smaller than that of another method at the same mesh size (i.e., C value of two methods may not be the same).Detailed comparisons show that, at the same mesh spacing, relative error of the results given by USTLBM is lower than that obtained by STLBM.Therefore, it is reasonable to expect more accurate results at the same mesh resolution when applying USTLBM.More cases in high Rayleigh numbers of Ra = 10 5 , 10 6 , 10 7 and 10 8 are simulated by both USTLBM and STLBM with different mesh resolutions.Note that, in cases of Ra = 10 7 and 10 8 , non-uniform meshes are adopted to capture severe changes of flow fields near boundaries.Lagrange interpolation algorithm is applied to determine functional values on locations where the streaming process initiates.Representative flow parameters are tabulated in Tables 1-4.It is noteworthy that, in all cases, both USTLBM and STLBM do not diverge at very coarse mesh size of 11 × 11, and reasonable results of flow parameters can be achieved.Such performance demonstrates that both USTLBM and STLBM have very good numerical stability, which makes them superior to standard LBM that can easily diverge in cases with high Reynolds/Rayleigh numbers.Convergence study is also carried out with various mesh sizes adopted in these cases.It turns out that the calculated flow parameters are approaching to the reference data by refining the meshes, which verifies that both USTLBM and STLBM can converge to the correct answer.More cases in high Rayleigh numbers of Ra = 10 5 , 10 6 , 10 7 and 10 8 are simulated by both USTLBM and STLBM with different mesh resolutions.Note that, in cases of Ra = 10 7 and 10 8 , non-uniform meshes are adopted to capture severe changes of flow fields near boundaries.Lagrange interpolation algorithm is applied to determine functional values on locations where the streaming process initiates.Representative flow parameters are tabulated in Tables 1-4.It is noteworthy that, in all cases, both USTLBM and STLBM do not diverge at very coarse mesh size of 11 × 11, and reasonable results of flow parameters can be achieved.Such performance demonstrates that both USTLBM and STLBM have very good numerical stability, which makes them superior to standard LBM that can easily diverge in cases with high Reynolds/Rayleigh numbers.Convergence study is also carried out with various mesh sizes adopted in these cases.It turns out that the calculated flow parameters are approaching to the reference data by refining the meshes, which verifies that both USTLBM and STLBM can converge to the correct answer.Another notable issue reflected in these tables is that at the same mesh size, the results obtained by USTLBM show smaller relative errors to the reference data.At high Rayleigh numbers, STLBM usually requires finer meshes to achieve comparable results.Here, we take the maximum vertical velocity for illustration.At Ra = 10 6 , results obtained by USTLBM are within the range given by reference data at the mesh size of 201 × 201, while STLBM can only get comparable results at the mesh size of 251 × 251.In the case of Ra = 10 8 , the relative error of results obtained by USTLBM at the mesh size of Appl.Sci.2017, 7, 277 14 of 21 301 × 301 is within 0.02%, while the relative error of STLBM at the same mesh size is 12.66%.These results evidently prove that, although in the same order of accuracy, USTLBM is more accurate than STLBM at the same mesh size.

Mixed Convection in a Lid-Driven Cavity
In this section, we present the example of mixed convection in a lid-driven cavity.The physical configuration of the problem is a square cavity with moving top lid and three stationary walls.The top and the bottom walls are maintained at constant cold and hot temperatures (T c and T h ), respectively (see Figure 4).Neumann boundary condition for the temperature with its normal gradient equal to zero is implemented on the right/left walls.Typical non-dimensional parameters defined for this problem include Reynolds number, Grashof number, Prandtl number and Richardson number, which are expressed as  Referring to set-ups in previous studies [26,27], here we set Gr = 10 6 and Pr = 0.71, and three cases with different Richardson numbers of 10, 1.0 and 0.1 are simulated.Uniform mesh with the size of 251 × 251 is adopted in our simulations.
Figure 5 presents isotherms and streamlines at steady state obtained by USTLBM.As can be seen, smoothness of the isotherms and streamlines is well maintained by the present USTLBM.It is also observed that, at larger Ri number (e.g., Ri = 10), two major vortices show up in the flow pattern.With the decreasing of the Ri number, the lower vortex is gradually compressed.When Ri number is sufficiently small (e.g., Ri = 0.1), only one major vortex exists and the flow pattern is quite similar to the classic lid-driven cavity flow at the same Reynolds number.The physical mechanism of such phenomenon is the relative importance of natural convection to the forced convection, which is the essence of Richardson number.At low Ri number, the effect of the natural convection is quite small or even negligible, due to which the final flow pattern is mainly shaped by the shear force induced by the moving lid.Therefore, the streamlines at low Ri number are very much like those of isothermal lid-driven cavity flows at the same Reynolds number.In high Ri number scenario, natural convection has comparable or even larger effect on the flow field as the forced convection does.Subsequently, the buoyance force induced by the temperature gradient enlarges the vortex on the right corner and eventually makes its size comparable to that of the main vortex.Referring to set-ups in previous studies [26,27], here we set Gr = 10 6 and Pr = 0.71, and three cases with different Richardson numbers of 10, 1.0 and 0.1 are simulated.Uniform mesh with the size of 251 × 251 is adopted in our simulations.
Figure 5 presents isotherms and streamlines at steady state obtained by USTLBM.As can be seen, smoothness of the isotherms and streamlines is well maintained by the present USTLBM.It is also observed that, at larger Ri number (e.g., Ri = 10), two major vortices show up in the flow pattern.With the decreasing of the Ri number, the lower vortex is gradually compressed.When Ri number is sufficiently small (e.g., Ri = 0.1), only one major vortex exists and the flow pattern is quite similar to the classic lid-driven cavity flow at the same Reynolds number.The physical mechanism of such phenomenon is the relative importance of natural convection to the forced convection, which is the essence of Richardson number.At low Ri number, the effect of the natural convection is quite small or even negligible, due to which the final flow pattern is mainly shaped by the shear force induced by the moving lid.Therefore, the streamlines at low Ri number are very much like those of isothermal lid-driven cavity flows at the same Reynolds number.In high Ri number scenario, natural convection has comparable or even larger effect on the flow field as the forced convection does.Subsequently, the buoyance force induced by the temperature gradient enlarges the vortex on the right corner and eventually makes its size comparable to that of the main vortex.
smoothness of the isotherms and streamlines is well maintained by the present USTLBM.It is also observed that, at larger Ri number (e.g., Ri = 10), two major vortices show up in the flow pattern.With the decreasing of the Ri number, the lower vortex is gradually compressed.When Ri number is sufficiently small (e.g., Ri = 0.1), only one major vortex exists and the flow pattern is quite similar to the classic lid-driven cavity flow at the same Reynolds number.The physical mechanism of such phenomenon is the relative importance of natural convection to the forced convection, which is the essence of Richardson number.At low Ri number, the effect of the natural convection is quite small or even negligible, due to which the final flow pattern is mainly shaped by the shear force induced by the moving lid.Therefore, the streamlines at low Ri number are very much like those of isothermal lid-driven cavity flows at the same Reynolds number.In high Ri number scenario, natural convection has comparable or even larger effect on the flow field as the forced convection does.Subsequently, the buoyance force induced by the temperature gradient enlarges the vortex on the right corner and eventually makes its size comparable to that of the main vortex.For quantitative comparison, the average Nusselt numbers on the hot wall at different Richardson numbers are presented in Table 5. Apparently, the computational results in our simulations agree well with the reference data in the literature.Therefore, the robustness of USTLBM in this mixed convection problem is validated.For quantitative comparison, the average Nusselt numbers on the hot wall at different Richardson numbers are presented in Table 5. Apparently, the computational results in our simulations agree well with the reference data in the literature.Therefore, the robustness of USTLBM in this mixed convection problem is validated.

Natural Convection in a Concentric Annulus
To test the flexibility of USTLBM on body-fitted mesh and on problems with curved boundaries, numerical test on natural convection in a concentric annulus is carried out in this section.Set-up of the problem is depicted in Figure 6.The gas filled between the concentric annulus is driven by the buoyance force that is induced by the temperature gradient between the hot inner circular cylinder and cold outer circular cylinder.Non-dimensional parameters in this problem include the aspect ratio Ar = R o /R i , the Prandtl number Pr = υ/χ and the Rayleigh number Ra = gβ(T i −T o )L 3 /υχ.In our simulations, we set Ar = 2.6 and Pr = 0.71; and four cases with different Rayleigh numbers of 10 3 , 6 × 10 3 , 10 4 and 5 × 10 4 are considered.To accurately describe the curved boundaries encountered in this problem, body fitted mesh with the mesh size of 251 × 81 is adopted in all simulations.Lagrangian interpolation algorithm is introduced in computational domain to evaluate functional values on locations where particles start the streaming process.problem, body fitted mesh with the mesh size of 251 × 81 is adopted in all simulations.Lagrangian interpolation algorithm is introduced in computational domain to evaluate functional values on locations where particles start the streaming process.Isotherms at steady state of different cases are presented in Figure 7.It can be seen that our simulations keep the smoothness of the temperature field, which validates the stability of USTLBM.It is also observed that, with the increase of Rayleigh number, the isotherms are more squeezed to the lower part of the inner circular cylinder and the upper part of the outer circular cylinder.The reason is that, at higher Rayleigh number, the buoyance force induced by the temperature gradient has larger effect on the flow field, due to which the gas is pushed upward and compressed to the lower part of the inner cylinder and the upper part of the outer cylinder.Isotherms at steady state of different cases are presented in Figure 7.It can be seen that our simulations keep the smoothness of the temperature field, which validates the stability of USTLBM.It is also observed that, with the increase of Rayleigh number, the isotherms are more squeezed to the lower part of the inner circular cylinder and the upper part of the outer circular cylinder.The reason is that, at higher Rayleigh number, the buoyance force induced by the temperature gradient has larger effect on the flow field, due to which the gas is pushed upward and compressed to the lower part of the inner cylinder and the upper part of the outer cylinder.Isotherms at steady state of different cases are presented in Figure 7.It can be seen that our simulations keep the smoothness of the temperature field, which validates the stability of USTLBM.It is also observed that, with the increase of Rayleigh number, the isotherms are more squeezed to the lower part of the inner circular cylinder and the upper part of the outer circular cylinder.The reason is that, at higher Rayleigh number, the buoyance force induced by the temperature gradient has larger effect on the flow field, due to which the gas is pushed upward and compressed to the lower part of the inner cylinder and the upper part of the outer cylinder.6 and compared with reference data in the literature [28,29].It is noted that our computational results are in good agreement with the reference data.The relative errors are within 3%, which demonstrates the robustness of USTLBM in the present example on body-fitted mesh and with curved boundaries.To quantify the correctness of USTLBM in this problem, the following equivalent heat conductivities are defined Calculated results in our simulations are tabulated in Table 6 and compared with reference data in the literature [28,29].It is noted that our computational results are in good agreement with the reference data.The relative errors are within 3%, which demonstrates the robustness of USTLBM in the present example on body-fitted mesh and with curved boundaries.

3D Natural Convection in a Cubic Cavity
To further test the performance of USTLBM in three-dimensional problems, here we present the example of 3D natural convection in a cubic cavity.The physical set-up of the problem is a cubic cavity with its left and right walls, respectively, fixed on constant temperatures of T h and T c .On the other walls, Neumann boundary condition with normal temperature gradient equal to zero is implemented.The non-dimensional parameters that determine the converged flow pattern are still the Rayleigh number and Prandtl number as defined in Equations ( 66) and (67).The average Nusselt number on the wall is defined as Considering the numerical efficiency, non-uniform Cartesian mesh is adopted in this example.Two cases with high Rayleigh numbers of Ra = 10 5 and 10 6 are simulated.
Tables 7 and 8 give the average Nusselt number on the wall, the maximum horizontal velocity along the vertical centerline and the maximum vertical velocity along the horizontal centerline at different mesh sizes and Rayleigh numbers.It is observed that, be refining the mesh resolution, both USTLBM and STLBM can converge to the reference data [30][31][32][33][34].Such performance indicates the robustness of USTLBM and STLBM on non-uniform mesh in three-dimensional scenarios.Careful comparisons reveal that, compared to our two-dimensional test of the natural convection in a square cavity (Section 4.1), the advantage of USTLBM over STLBM is not so obvious in this three-dimensional example.The possible reason is that, since the non-uniform mesh is adopted here, errors are inevitably introduced in the interpolation process, due to which the advantage of USTLBM in lowering numerical errors is counteracted and the results of USTLBM and STLBM are quite close.If uniform mesh is applied, the present USTLBM performs better than STLBM under the same mesh size, especially on coarser meshes (see Table 9).Such evidence indicates that USTLBM can lower the numerical error of STLBM while maintaining the same order of accuracy in space.

Conclusions
In this paper, we present an unconditionally stable thermal lattice Boltzmann method (USTLBM) for numerical simulation of 2D/3D incompressible thermal flows.This method is developed from the recently proposed simplified thermal lattice Boltzmann method (STLBM) by truly ensuring the second-order of accuracy in time.Specifically, USTLBM is derived by reconstructing solutions to the macroscopic governing equations recovered from lattice Boltzmann equation and resolved in a predictor-corrector scheme.Final formulations of USTLBM only involve the equilibrium and the non-equilibrium distribution functions.In practical computation, the equilibrium distribution function is evaluated by the macroscopic variables, while the non-equilibrium term is approximated by the difference between two equilibrium distribution functions at different locations and time levels.Such treatment means that USTLBM can directly track the evolution of macroscopic variables rather than the distribution functions, which lowers the cost in virtual memories and facilitates implementation of physical boundary conditions.
Von Neumann stability analysis can be performed to prove that USTLBM is an unconditionally stable scheme.Numerical experiments are carried out on cases with coarse meshes and in high Rayleigh numbers to validate this conclusion.Accuracy test of USLTBM reveals that it ensures the second-order of accuracy in space while being able to achieve more accurate results than STLBM at the same mesh size.
Four numerical examples, including 2D natural convection in a square cavity, mixed convection in a lid-driven cavity, natural convection in a concentric annulus and 3D natural convection in a cubic cavity, are presented to validate the robustness and the flexibility of USTLBM.It turns out that the results obtained by USTLBM agree well with reference data in the literature; and the method can be flexibly coupled with Lagrangian interpolation algorithm to handle problems on non-uniform or body-fitted meshes and with curved boundaries.

( 7 )
Repeat Steps (2)-(6) until the computation converges or the prescribed maximum time step is reached.

Figure 2 .
Figure 2. Physical configuration of 2D natural convection in a square cavity.

Figure 2 .
Figure 2. Physical configuration of 2D natural convection in a square cavity.

Figure 3 .
Figure 3. Relative errors of volume averaged Nusselt numbers at different mesh spacing.

Figure 3 .
Figure 3. Relative errors of volume averaged Nusselt numbers at different mesh spacing.

Figure 4 .
Figure 4. Physical configuration of mixed convection in a lid-driven cavity.

Figure 4 .
Figure 4. Physical configuration of mixed convection in a lid-driven cavity.

Figure 6 .
Figure 6.Illustration of the setup for natural convection in a concentric annulus.

Figure 6 .
Figure 6.Illustration of the setup for natural convection in a concentric annulus.

Figure 6 .
Figure 6.Illustration of the setup for natural convection in a concentric annulus.

21 Figure 7 .
Figure 7. Isotherms for natural convection in a concentric annulus.

Figure 7 .
Figure 7. Isotherms for natural convection in a concentric annulus.
The weighting coefficient ωα and the sound speed cs are αx y 2α + e αy y 3α + s e T; and the superscript n represents the time level.For illustration purposes, we take the term [∂y 2 * /∂y 2 ] n to elaborate the analyzing process: * = F(Y n ).To linearize the problem to perform von Neumann stability analysis, the obtained equations are differentiated with respect to y 1 , y 2 , y 3 and y 4 , which leads to dY * = * )

Table 1 .
Comparisons of flow parameters in Ra = 10 5 at different mesh sizes.

Table 2 .
Comparisons of flow parameters in Ra = 10 6 at different mesh sizes.

Table 3 .
Comparisons of flow parameters in Ra = 10 7 at different mesh sizes.

Table 4 .
Comparisons of flow parameters in Ra = 10 8 at different mesh sizes.

Table 5 .
Comparison of average Nusselt number on hot wall for the mixed convection in lid-driven cavity.

Table 5 .
Comparison of average Nusselt number on hot wall for the mixed convection in lid-driven cavity.

Table 6 .
Comparison of average equivalent heat conductivity at different Rayleigh numbers.

Table 6 .
Comparison of average equivalent heat conductivity at different Rayleigh numbers.

Table 7 .
Comparisons of flow parameters for 3D natural convection in Ra = 10 5 .

Table 8 .
Comparisons of flow parameters for 3D natural convection in Ra = 10 6 .

Table 9 .
Average Nusselt number on the wall for 3D natural convection in Ra = 10 5 .