Simulating Turbulent Buoyant Flow by a Simple LES-Based Thermal Lattice BoltzmannModel

To simulate turbulent buoyant flow in geophysical science, where usually the vorticity-streamfunction equations instead of the primitive-variables Navier-Stokes equations serve as the governing equations, a novel and simple thermal lattice Boltzmann model is proposed based on large eddy simulation (LES). Thanks to its intrinsic features, the present model is efficient and simple for thermal turbulence simulation. Two-dimensional numerical simulations of natural convection in a square cavity were performed at high Rayleigh number varying from 10 to 10 with Prandtl number at 0.7. The advantages of the present model are validated by numerical experiments.


Introduction
In the last two decades, the lattice Boltzmann method (LBM) has proved its capability to simulate a large variety of fluid flows [1][2][3][4][5][6][7][8].Unlike the conventional numerical methods based on discretizations of macroscopic continuum equations, the LBM is based on microscopic models and mesoscopic kinetic equations for fluids.The kinetic nature of LBM makes it very suitable for fluid systems [3].In particular, the LBM has been compared favourably with the spectral method [9], the artificial compressibility method [10], the finite volume method [11,12], and finite difference method [13]; all quantitative results further validate excellent performances of the LBM not only in computational efficiency but also in computational accuracy [14].
The LBM originally was designed to simulate isothermal flows [1][2][3].Later it was extended for thermal flows.Some of the earliest works in thermal LBM include those of Massaioli et al. [15], Alexander et al. [16], Chen et al. [17], and Eggels and Somers [18], to cite only a few.Thermal lattice Boltzmann (LB) models can be classified into three categories based on their approach in solving the Boltzmann equation, namely, the multispeed approach, the passive scalar approach, and the double-population approach [19][20][21].The multispeed approach is an extension of the common LB equation isothermal models in which only the density distribution function is used.To obtain the temperature evolution equation at the macroscopic level, additional speeds are necessary and the equilibrium distribution must include the higher-order velocity terms.However, the inclusion of higher-order velocity terms leads to numerical instabilities and hence the temperature variation is limited to a narrow range [22].For the passive scalar approach, it consists in solving the velocity by the LBM and the macroscopic temperature equation independently [23].The macroscopic temperature equation is similar to a passive scalar evolution equation if the viscous heat dissipation and compression work done by pressure are negligible.The coupling to LBM is made by adding a potential to the distribution function equation.This approach has attracted much attention compared to the multispeed approach on account of its numerical stability.But the main disadvantage is that viscous heat dissipation and compression work done by pressure cannot be incorporated in this model.The doublepopulation approach introduces an internal energy density distribution function in order to simulate the temperature field, the velocity field still being simulated using the density distribution function [24].Compared with the multispeed thermal LB approach, the numerical scheme is numerically more stable.Moreover, this method is able to include the ISRN Thermodynamics viscous heat dissipation and compression work done by pressure.
Turbulent buoyant flow is of fundamental interest and practical importance [25,26].However, to date the studies using the LBM on turbulent buoyant flow are quite few [19,20,[27][28][29][30].Most of them took the LBM as a direct numerical simulation (DNS) tool to investigate turbulent buoyant flow.For instance, Shi and Guo employed the LBM to simulate turbulent natural convection due to internal heat generation [27].Zhou and his cooperators used an LBMbased commercial software "PowerFLOW" to simulate natural and forced convection [28].Dixit and Babu simulated the fully turbulent natural convection with very high Rayleigh Ra numbers by the LBM [19].More recently, Kuznik and his cooperators used the LBM to simulate natural convection up to Ra = 10 8 with nonuniform mesh [20].But the DNS is too expensive for available computer capability to simulate practical problems.Due to the balance between accuracy and efficiency, several LB models based on large eddy simulation (LES) were developed as an alternative [29,30].The spirit of LES-based LB models is to split the effective fluid viscosity ν e into two parts ν 0 and ν t .ν 0 is the molecular viscosity and ν t the eddy viscosity [30].The effective lattice relaxation time τ e depends on ν e .Treeck and his cooperators combined the LBM and the LES to investigate turbulent convective flows [29] and Liu et al. designed a thermal Bhatnagar-Gross-Krook (BGK) model based on the LES to simulate turbulent natural convection due to internal heat generation up to Ra = 10 13 [30].
All existing LES-based LB models for turbulent buoyant flow are based on primitive-variables Navier-Stokes (NS) equations.Although they have been successfully used in many applications, the disadvantages of the primitivevariables LES-based LB models are obvious.First, in order to recover the correct macroscopic governing equations, in the lattice evolving equation the treatment of buoyant forcing due to temperature field is complicated: one has to redefine the fluid velocity as well as the equilibrium velocity and expands the forcing in a power series in the particle velocity [31].Second, the calculation of τ e of the primitive-variables LES-based LB models is generally complicated, not only for the BGK approximation but also for multiple-relaxationtime models [29], and in some cases to obtain the exact value of τ e is extremely difficult [30,32].In particular, they are inconvenient for geophysical flow simulations where the vorticity-streamfunction equations, instead of the primitivevariables NS equations, serve as the governing equations [33].
To overcome the above disadvantages, in this paper a novel and simple LES-based thermal LB model, which is an extension of the model designed in our previous work [34,35], is proposed to simulate turbulent buoyant flow.Unlike existing primitive-variables LES-based LB models, for the flow field, the target macroscopic equations of the present model are vorticity-streamfunction equations.Because the vorticity-streamfunction equations consist in an advectiondiffusion equation and a Poisson equation, for which the source terms in the LBM can be treated more simply than the force strategy for the primitive-variables-based LB equation with additional forcing [31,34,35], therefore the first disadvantage of the primitive-variables LES-based LB models is overcame.In addition, the calculation of τ e of the present model is much simpler and more efficient than that of primitive-variables LES-based LB models due to the intrinsic features of the present model.We will show them in detail in Section 3.More important, the present model can be used directly to simulate the problems where the vorticitystreamfunction equations serve as the governing equations such as in geophysical science.
The rest of the paper is organized as follows.Vorticitystreamfunction-based governing equations for turbulent buoyant flow are presented in Section 2. In Section 3, a novel and simple LES-based thermal LB model is introduced.In Section 4, numerical experiments are performed to validate the present model.Conclusion is presented in the last section.

Governing Equations
The turbulent buoyant flow is governed by unsteady Boussinesq equations.The primitive-variables dimensionless equations read [27,29] The normalizing characteristic quantities are length with H, velocity with (α/H)Ra 0.5 , pressure with ρ(α/H) 2 Ra, and time with H 2 /αRa −0.5 .H, ρ, p, and α are the characteristic length, density, pressure, and coefficient of thermal expansion, respectively.The direction of the gravity g is along the negative y-axis.u = (u, v) is the velocity vector.ν e and D e are the effective viscosity and the effective thermal diffusivity, respectively.A complete description of the scaling procedure could be found in [36].
The effective viscosity ν e = ν 0 + ν t .The molecular viscosity ν 0 = Pr Ra −0.5 .The eddy viscosity ν t is computed from the local shear rate and a length scale when the Smagorinsky model is used [30,37,38] where the constant C is called the Smagorinsky constant and is adjustable.In our simulations, we take C = 0.1.Δ is the filter width Δ = Δx.Δx is the lattice grid spacing.The local intensity of the strain rate tensor is defined as S αβ is the strain rate tensor.The effective thermal diffusivity D e = D 0 + D t , where D 0 = Ra −0.5 and D t = υ t /Pr t .Pr t = 0.4 is the turbulent Prandtl number.
When the vorticity ω and the streamfunction ψ are introduced, (1) can be transformed to The velocities u and v are obtained from

LES-Based Thermal LB Model
Equation (4) (governing equation for the flow field) and ( 5) (governing equation for the temperature field) are nothing but advection-diffusion equations with/without source terms.There are many matured efficient lattice Boltzmann models for this type of equation [34].In this paper a D2Q5 model is employed to solve these equations.It reads where e k (k = 0 • • • 4) are the discrete velocity directions: c = Δx/Δt is the fluid particle speed.Δt and τ e are the time step and the dimensionless relaxation time, respectively.Υ o,k is the discrete form of the source term Υ o = Pr(∂T/∂x), 0, for (4) and ( 5), respectively.The simplest choice satisfying the constraint equation ( 10) is One can see that in the present model redefining the fluid velocity as well as the equilibrium velocity and expanding the forcing in a power series in the particle velocity, which result from the treatment of buoyant forcing due to temperature field, both are avoided.The equilibrium distribution g (eq) k is defined by δ = ω, T, for (4) and ( 5), respectively, and is obtained by The dimensionless relaxation time τ e for (4) is determined by where τ 0 = ν 0 /c 2 s Δt + 0.5 and c 2 s = (2/5)c 2 .The dimensionless relaxation time τ e for ( 5) is determined by where τ 0 = D 0 /c 2 s Δt + 0.5.The complication of calculation of τ e of the primitivevariables LES-based LB models results from the complication of calculation of |S| (cf.[30, equation (12)], [37, equation (9)] and [38, equation (22)]).Fortunately, thanks to the intrinsic features of vorticity-streamfunction equations, the calculation of |S| is very simple in the present model, just (16) please see Appendix B for the details.Because the value of vorticity ω is already known at every grid point, therefore compared with primitive-variables LES-based LB models [29,30] no extra computational cost is needed for |S| in the present model.Equation ( 6) is just the Poisson equation, which also can be solved by the LB method efficiently.In the present study, the D2Q5 model used in our previous work [34] is employed because this model is more efficient and more accurate than others to solve the Poisson equation.The evolution equation for (6) reads where is the equilibrium distribution function and defined by ξ k and ζ k are weight parameters given as The detailed derivation from ( 8) and ( 17) to ( 4)-( 6) can be found in Appendix A.
In the present model (7) are solved by the central finite difference scheme.

Results and Discussions
Because it is very difficult to find a simple benchmark test in geophysical science, in the present study, the natural convection in a square cavity with Pr = 0.7 and 10 4 ≤ Ra ≤ 10 10 is invoked to validate the present model.The computational domain is illustrated in Figure 1.The nonequilibrium extrapolation method [30,34] is used for boundary treatment.The grid resolution is 256 × 256 for all cases.Because the flow field becomes unsteady when Ra ≥ 10 7 , the iteration step of 4 × 10 5 is used to guarantee the flow field is fully developed.

The Laminar Flow
Ra ≤ 10 6 .Figures 2-4 show the isothermal and streamfunction contours at 10 4 ≤ Ra ≤ 10 6 .For low Ra ≤ 10 4 a central vortex appears as the dominant characteristic of the flow.As Ra increases, the vortex tends to become elliptic and finally breaks up into two vortices at Ra = 10 5 .When Ra continues increasing, the two vortices move towards the walls, giving space for a third vortex to develop.Even for higher Ra = 10 6 , the third vortex is very weak in comparison with the other two.The shape of the isotherms shows how the dominant heat transfer mechanism changes as Ra increases.For low Ra almost vertical isotherms appear, because heat is transferred by conduction between hot and cold walls.As the isotherms depart from the vertical position, the heat transfer mechanism changes from conduction to convection.Table 1 reports the average Nusselt (Nu) number at the hot wall, together with that obtained in previous studies [19,40].The thermal and the flow fields agree very well with those reported in the literature [19,20,39].The transitional flow features reported by previous studies [20,39] are well captured by our model.Table 2 reports the average Nusselt number at the hot wall, together with that obtained in previous studies [19,41].The present results are in excellent agreement with those in [19,41].

4.3.
The Fully Turbulent Flow 10 9 ≤ Ra ≤ 10 10 .When Ra = 10 9 , the flow becomes fully turbulent.The flow structure in entire simulation domain becomes irregular and chaotic.However, the isothermal curves become almost straight at the center and very sharp inside the very thin boundary layers, as shown in Figure 7(a).When the Ra is increased to Ra = 10 10 , the flow becomes even more turbulent.The small-scale structures become increasing and much finer.And the isotherms except those which are very close to the midplane of the cavity are significantly deformed by the turbulent flow and not straight longer, as Figure 7(b) illustrates.
Table 3 reports the average Nusselt number at the hot wall, together with that obtained in previous studies [19,26].The deviations of Nu between different models result from the effect of the eddy viscosity ν t becoming significant when Ra ≥ 10 9 .

Conclusion
In order to simulate turbulent buoyant flow in geophysical science, where usually the vorticity-streamfunction equations instead of the primitive-variables NS equations serve as the governing equations, we propose a novel and simple thermal LB model based on LES.Unlike existing LES-based LB models for thermal turbulent flow simulation, which were based on primitive-variables NS equations, the target macroscopic equations of the present for flow field model are   vorticity-streamfunction equations.Thanks to its intrinsic features, the buoyant forcing term due to the temperature in the present model can be treated more simply than that in the existing LES-based LB models for thermal turbulence, avoiding re-defining the fluid velocity as well as the equilibrium velocity and avoiding expanding the forcing in a power series in the particle velocity.Moreover, the calculation of τ e of the present model is much simpler and more efficient than that of primitive-variables LES-based LB models due to the intrinsic features of the present model.Therefore, the disadvantages of the existing LES-based LB models are overcome by the present model.Furthermore, the present model can be used directly to simulate the problems where the vorticity-streamfunction equations instead of primitivevariables NS equations serve as the governing equations, such as geophysical flow simulations.
The present model is validated by natural convection in a square cavity at Ra from 10 4 to 10 10 with Pr = 0.7.The results obtained by the present model are in excellent agreement with those in previous publications.When Ra < 10 9 , the average Nusselt numbers at the hot wall obtained by the present model agree well with previous studies.With the increasing of Ra, the effect of the eddy viscosity ν t becomes significant so that the deviations of Nu between different models on the boundaries set in.

A. Recovery of the Governing Equations
The recovery of the governing equations is aided by the Chapman-Enskog expansion.Expanding the distribution   functions and the time and space derivatives in terms of a small quantity To perform the Chapman-Enskog expansion we must first Taylor expand (8): where where And then, we can obtain the following equations in consecutive order of the parameter : g (2)  k .(A.6)   Combining (A.9) and (A.10), we can obtain (4) and ( 5) if δ = ω, T, respectively.The detailed derivation from ( 17) to (6) can be found in our previous work [34].with second-order accuracy, consistent with the numerical accuracy of the LB method.

Table 1 :
Comparison of average Nusselt number at the hot wall of laminar flow with previous works.

Table 2 :
Comparison of average Nusselt number at the hot wall of transitional flow with previous works.

Table 3 :
Comparison of average Nusselt number at the hot wall of turbulent flow with previous works.The Transitional Flow 10 7 ≤ Ra ≤ 10 8 .At higher Ra = 10 7 and 10 8 , the velocities at the center of the cavity are very small compared with those at the boundaries where the fluid is moving fast, forming vortices at the lower right and top left corner of the cavity, destabilizing the laminar flow, as Figures5 and 6illustrate.The vortices become narrow when Ra increases, improving the stratification of the flow at the central part of the cavity.The isotherms at the center of the cavity are horizontal and become vertical near the walls.