# Simulating Turbulent Buoyant Flow by a Simple LES-Based Thermal Lattice Boltzmann Model

Authors: Sheng Chen
Journal: ISRN Thermodynamics
 (2012)

---

## Abstract

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 from104 to 1010 with Prandtl number at 0.7. The advantages of the present model are validated by numerical experiments.


---

## Body

## 1. Introduction

In the last two decades, the lattice Boltzmann method (LBM) has proved its capability to simulate a large variety of fluid flows [1–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–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–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 double-population 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 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–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 LBM-based 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=108 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=1013 [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 primitive-variables 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-relaxation-time 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 primitive-variables 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 advection-diffusion 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 vorticity-streamfunction equations serve as the governing equations such as in geophysical science.The rest of the paper is organized as follows. Vorticity-streamfunction-based governing equations for turbulent buoyant flow are presented in Section2. 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.


## 2. Governing Equations

The turbulent buoyant flow is governed by unsteady Boussinesq equations. The primitive-variables dimensionless equations read [27, 29](1)∇⋅u⃗=0,∂u⃗∂t+u⃗⋅∇u⃗=-∇p+νeΔu⃗+PrTg⃗|g⃗|,∂T∂t+u⃗⋅∇T=DeΔT.
The normalizing characteristic quantities are length with H, velocity with (α/H)Ra0.5, pressure with ρ(α/H)2Ra, and time with H2/α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 De 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=PrRa-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](2)νt=(CΔ)2|S¯|,
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(3)|S̅|=2SαβSαβ.Sαβ is the strain rate tensor.The effective thermal diffusivityDe=D0+Dt, where D0=Ra-0.5 and Dt=υt/Prt. Prt=0.4 is the turbulent Prandtl number.When the vorticityω and the streamfunction ψ are introduced, (1) can be transformed to (4)∂ω∂t+u∂ω∂x+v∂ω∂y=υe(∂2ω∂x2+∂2ω∂y2)+Pr∂T∂x,(5)∂T∂t+u∂T∂x+v∂T∂y=De(∂2T∂x2+∂2T∂y2),(6)∂2ψ∂x2+∂2ψ∂y2=-ω.
The velocities u and v are obtained from(7)u=∂ψ∂y,v=-∂ψ∂x.


## 3. 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(8)gk(x⃗+ce⃗kΔt,t+Δt)-gk(x⃗,t)=-τe-1[gk(x⃗,t)-gk(eq)(x⃗,t)]+ΔtΥo,k,
where e⃗k (k=0···4) are the discrete velocity directions:(9)e⃗k={(0,0):k=0,(cos⁡(k-1)π2,sin⁡(k-1)π2):k=1,2,3,4.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. Υo,k satisfies(10)∑k≥0Υo,k=Υo,Υo=Pr(∂T/∂x),0, for (4) and (5), respectively.The simplest choice satisfying the constraint equation (10) is(11)Υo,k=Υo5.
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 distributiongk(eq) is defined by(12)gk(eq)=δ5[1+2.5e⃗k⋅u⃗c]δ=ω,T, for (4) and (5), respectively, and is obtained by(13)δ=∑k≥0gk.
The dimensionless relaxation time τe for (4) is determined by(14)τe=τ0+(CΔ)2|S¯|cs2Δt,
where τ0=ν0/cs2Δt+0.5 and cs2=(2/5)c2.The dimensionless relaxation timeτe for (5) is determined by(15)τe=τ0+(CΔ)2|S¯|/Prtcs2Δt,
where τ0=D0/cs2Δt+0.5.The complication of calculation ofτe of the primitive-variables 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)|S¯|=|ω|
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 (17)fk(x⃗+ce⃗kΔt,t+Δt)-fk(x⃗,t)=Ωk+Ωk′,
where Ωk=-τψ-1[fk(x⃗,t)-fk(eq)(x⃗,t)],  Ωk′=-ΔtζkωD, and D=(c2/2)(0.5-τψ). τψ>0.5 is the dimensionless relaxation time. fk(eq) is the equilibrium distribution function and defined by(18)fk(eq)={(ξ0-1.0)ψ:k=0,ξkψ:k=1,2,3,4.ξk and ζk are weight parameters given as ξ0=ζ0=0, ξk=ζk=1/4(k=1···4). ψ is obtained by(19)ψ=∑k≥1fk.
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.


## 4. 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 104 ≤ Ra ≤ 1010is 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≥107, the iteration step of 4×105 is used to guarantee the flow field is fully developed.Figure 1
The configuration of computational domain.


### 4.1. The Laminar FlowRa≤106

Figures2–4 show the isothermal and streamfunction contours at 104≤Ra≤106. For low Ra≤104 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=105. When Ra continues increasing, the two vortices move towards the walls, giving space for a third vortex to develop. Even for higher Ra=106, 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].Table 1
Comparison of average Nusselt number at the hot wall of laminar flow with previous works.
RaReference [19]Reference [40]Present1042.2472.2382.22971054.52264.5094.51181068.8058.8178.7776Isothermal (a) and streamfunction (b) contours atRa=104.
(a)(b)Isothermal (a) and streamfunction (b) contours atRa=105.
(a)(b)Isothermal (a) and streamfunction (b) contours atRa=106.
(a)(b)


### 4.2. The Transitional Flow107≤Ra≤108

At higherRa=107 and 108, 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 Figures 5 and 6 illustrate. 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. 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].Table 2
Comparison of average Nusselt number at the hot wall of transitional flow with previous works.
RaReference [19]Reference [41]Present10716.7916.52316.735210830.50630.22530.3177Isothermal (a) and streamfunction (b) contours atRa=107.
(a)(b)Isothermal (a) and streamfunction (b) contours atRa=108.
(a)(b)


### 4.3. The Fully Turbulent Flow109≤Ra≤1010

WhenRa=109, 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).Isothermal contours at (a) Ra = 109 and (b) Ra = 1010.
(a)(b)When theRa is increased to Ra=1010, 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.Table3 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≥109.Table 3
Comparison of average Nusselt number at the hot wall of turbulent flow with previous works.
RaReference [19]Reference [26]Present10957.350—49.52351010103.663156.8598.4100


## 4.1. The Laminar FlowRa≤106

Figures2–4 show the isothermal and streamfunction contours at 104≤Ra≤106. For low Ra≤104 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=105. When Ra continues increasing, the two vortices move towards the walls, giving space for a third vortex to develop. Even for higher Ra=106, 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].Table 1
Comparison of average Nusselt number at the hot wall of laminar flow with previous works.
RaReference [19]Reference [40]Present1042.2472.2382.22971054.52264.5094.51181068.8058.8178.7776Isothermal (a) and streamfunction (b) contours atRa=104.
(a)(b)Isothermal (a) and streamfunction (b) contours atRa=105.
(a)(b)Isothermal (a) and streamfunction (b) contours atRa=106.
(a)(b)


## 4.2. The Transitional Flow107≤Ra≤108

At higherRa=107 and 108, 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 Figures 5 and 6 illustrate. 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. 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].Table 2
Comparison of average Nusselt number at the hot wall of transitional flow with previous works.
RaReference [19]Reference [41]Present10716.7916.52316.735210830.50630.22530.3177Isothermal (a) and streamfunction (b) contours atRa=107.
(a)(b)Isothermal (a) and streamfunction (b) contours atRa=108.
(a)(b)


## 4.3. The Fully Turbulent Flow109≤Ra≤1010

WhenRa=109, 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).Isothermal contours at (a) Ra = 109 and (b) Ra = 1010.
(a)(b)When theRa is increased to Ra=1010, 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.Table3 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≥109.Table 3
Comparison of average Nusselt number at the hot wall of turbulent flow with previous works.
RaReference [19]Reference [26]Present10957.350—49.52351010103.663156.8598.4100


## 5. 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 primitive-variables 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 from104 to 1010 with Pr=0.7. The results obtained by the present model are in excellent agreement with those in previous publications. When Ra<109, 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.



---

*Source: 253972-2012-02-29.xml*