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 from to with Prandtl number at 0.7. The advantages of the present model are validated by numerical experiments.
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 numbers by the LBM [19]. More recently, Kuznik and his cooperators used the LBM to simulate natural convection up to 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 into two parts and . is the molecular viscosity and the eddy viscosity [30]. The effective lattice relaxation time depends on . 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 [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 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 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 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 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.
2. 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 , velocity with , pressure with , and time with . , , , and are the characteristic length, density, pressure, and coefficient of thermal expansion, respectively. The direction of the gravity is along the negative -axis. is the velocity vector. and 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 . The molecular viscosity . The eddy viscosity is computed from the local shear rate and a length scale when the Smagorinsky model is used [30, 37, 38] where the constant is called the Smagorinsky constant and is adjustable. In our simulations, we take . is the filter width . is the lattice grid spacing. The local intensity of the strain rate tensor is defined as is the strain rate tensor.
The effective thermal diffusivity , where and . is the turbulent Prandtl number.
When the vorticity and the streamfunction are introduced, (1) can be transformed to The velocities and are obtained from
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 where () are the discrete velocity directions: is the fluid particle speed. and are the time step and the dimensionless relaxation time, respectively. is the discrete form of the source term . satisfies, 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 is defined by, for (4) and (5), respectively, and is obtained by The dimensionless relaxation time for (4) is determined by where and .
The dimensionless relaxation time for (5) is determined by where .
The complication of calculation of of the primitive-variables LES-based LB models results from the complication of calculation of (cf. [30, equation (12)], [37, equation (9)] and [38, equation (22)]). Fortunately, thanks to the intrinsic features of vorticity-streamfunction equations, the calculation of is very simple in the present model, just 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 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 , , and . is the dimensionless relaxation time. is the equilibrium distribution function and defined by and are weight parameters given as , . is obtained by 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 for all cases. Because the flow field becomes unsteady when , the iteration step of is used to guarantee the flow field is fully developed.
4.1. The Laminar Flow
Figures 2–4 show the isothermal and streamfunction contours at . For low a central vortex appears as the dominant characteristic of the flow. As increases, the vortex tends to become elliptic and finally breaks up into two vortices at . When continues increasing, the two vortices move towards the walls, giving space for a third vortex to develop. Even for higher , 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 increases. For low 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 () 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].
(a)
(b)
(a)
(b)
(a)
(b)
4.2. The Transitional Flow
At higher and , 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 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].
(a)
(b)
(a)
(b)
4.3. The Fully Turbulent Flow
When , 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).
(a)
(b)
When the is increased to , 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 between different models result from the effect of the eddy viscosity becoming significant when .
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 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 from to with . The results obtained by the present model are in excellent agreement with those in previous publications. When , the average Nusselt numbers at the hot wall obtained by the present model agree well with previous studies. With the increasing of , the effect of the eddy viscosity becomes significant so that the deviations of between different models on the boundaries set in.
Appendices
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 . Substituting (A.1) into (A.2), we get where . And then, we can obtain the following equations in consecutive order of the parameter : Equation (A.6) can be simplified by (A.5): Because we can obtain where where Combining (A.9) and (A.10), we can obtain (4) and (5) if , respectively.
The detailed derivation from (17) to (6) can be found in our previous work [34].
B. Calculation of in the Present Model
For two-dimensional problems, , so With the aid of the continuum equation and the definition of vorticity Equation (B.1) is reduced as For incompressible flow, there exists [42] With the aid of (B.5), (B.4) is further reduced as Consequently, with second-order accuracy, consistent with the numerical accuracy of the LB method.
Acknowledgments
This work was partially supported by the Alexander von Humboldt Foundation, Germany, the Research Fund for the Doctoral Program of Higher Education of China (Grant no. 20100142120048), and the National Natural Science Foundation of China (Grant nos. 51006043 and 51176061).