The description of turbulent flow is one of the major problems in engineering and theoretical fluid mechanics. In turbulent flows, motions with large and small scales of eddies and dissipative structures coexist.

Numerical analysis of turbulent flows has been conducted in many studies [1–3]. Various models have been developed to describe turbulent flows, and one of the most popular models is that of Large-Eddy Simulation (LES) [4], [5] The LES model is widely used in many fields, including the analysis of geophysical phenomena in the atmosphere and ocean, and is considered one of the fundamental models of turbulent flows in fluid mechanics and physics [6–9]. However, the analysis of turbulent flows requires specific information that can be obtained by numerical analysis, which leads to the development of new algorithms and the use of parallel computing tools. The selection of an optimal numerical method that balances efficiency and accuracy is one of the obvious problems.

The Lattice Boltzmann Method (LBM) has gained great popularity over the last two decades. Apart from its simplicity and efficiency, it is used for simulating complex flow problems in different boundary conditions, such as multiphase and multicomponent flows, turbulent flows, and micro-scale flows [10–16], [17–19]. The LBM has shown its potential in simulating turbulent flows efficiently, and its application in computational fluid dynamics has some benefits compared with other methods for simulating turbulence flows [20–24]. Moreover, it has greater generality because of the advantageous balance between efficiency and accuracy. Yu et al. [25] used multiple relaxation times in the lattice Boltzmann equation for simulating turbulence flows by the method of Large Eddy Simulation (LES) and showed that it is a suitable method for simulating turbulence flows. Fernandino et al. [26] performed LES for simulating turbulent open-channel flows by the LBM and found good agreement with experimental results qualitatively. Chen [27] presented a novel and simple LES model based on the LBM for simulating two-dimensional turbulence flows and showed its efficiency, stability, and simplicity in simulating two-dimensional turbulence flows.

Kareem et al. [28] in their research paper proposed an LBM-based framework with a scheme for the simulation of convergent homogeneous turbulent flows. The authors observed that the turbulence characteristics of the flow, regardless of the model used, are achieved by methods similar to those used in the simulation of the flow in previous investigations. In the study by Nee [29], the author presented the application of hybrid lattice Boltzmann models for simulating natural convection heat transfer up to Rayleigh numbers of 1012. Yu et al. [25] presented the application of multiple relaxation times in the lattice Boltzmann equation for simulating LESs of turbulent flows and found it applicable for the simulation of such flows. Chen et al. [30] developed a new thermal lattice Boltzmann model for the simulation of natural convection in fluids with a high Prandtl number.

Chen [30] developed a simple and innovative LES model based on the LBM for the simulation of two-dimensional turbulent flows. The authors proved that the model is efficient, stable, and simple for the simulation of two-dimensional turbulence. Polasanapalli and Anupindi [31] developed a lattice Boltzmann-based solver for the simulation of turbulent natural convection in cylindrical cavities by the LES method. Moreover, Han et al. [32] performed a numerical analysis on the LES in the LBM for the simulation of airflow in a non-isothermal enclosed cavity. Agoujil et al. [33] simulated the process of natural convection with the LBM, considering a localized heating process within a cavity. Moreover, Lin et al. [34] derived a set of scaling laws for the boundary layer of a natural convection process, which can be employed for the validation of turbulence models.

The primary objective of this study is to introduce a turbulence model based on the LBM, considering a simple and explicit form. For this purpose, turbulent natural convection is simulated within a cavity, considering a wide range of Rayleigh numbers. Initially, a brief introduction is made about the LBM, followed by the application of LES within this method, and then the results are compared with those of earlier studies.

The geometric configuration of interest is given in Fig. 1 and consists of a two-dimensional square cavity of height H and width W, where H = W. The fluid inside is air with a Prandtl number of Pr = 0.71 and is considered a Newtonian fluid with constant properties, except for its density in the buoyancy term as given by the Boussinesq approximation .

The Mach number is set at Ma = 0.1 and is used to maintain an incompressible fluid flow regime in the implementation of the LBM. The fluid inside is air with Prandtl number Pr = 0.71, as is common in natural convection problems.

In the lattice Boltzmann model of incompressible thermal flows, two distribution functions, namely f and g, are used to describe the flow and temperature fields, respectively. These functions are used to compute the macroscopic flow variables, which include velocity, temperature, and pressure. An important characteristic of the LBM is that all parameters are made dimensionless. In this study, a square grid and the D2Q9 model, shown in Fig. 2, are used for both flow and temperature distribution functions.

By discretizing the Navier-Stokes equations, the governing equations for the flow and temperature distribution functions are written as follows [25], [35-37]:

For the flow field:

1 eq

fi(x + ciΔT, t + Δt) − fi(x, t) = [fi(x, t) − fi (x, t)] + ΔtciF (1)

τν For the temperature field:

In these relations, cₖ is the discrete velocity vector, Δt is the time step, and τ_f and τ_g are the relaxation times for the flow and temperature fields, respectively. fₖeq and gₖeq are the local equilibrium distribution functions, which are calculated according to the following relations for the flow and temperature fields. Furthermore, F represents the external force.

The discrete velocity vector cₖ for the D2Q9 model is defined as follows:

0i = 0 c (cos[(i − 1) ] , sin[(i − 1) ]) i = 1 − 4

ci = 2 2 (5)

π π π π c√2 (cos[(i − 5) + ] , sin[(i − 5) + ]) i = 5 − 8 { 2 4 2 4 In these relations, c = Δx/Δt = Δy/Δt, where Δx and Δy are the grid spacings in the x and y directions. To improve numerical stability, Tm, the average temperature value, is used to calculate cₛ. Using the Chapman-Enskog expansion, the Navier-Stokes equations can be recovered with the proposed model. The kinematic viscosity and thermal diffusion coefficient are related to the relaxation times as follows:

1 2 1 2

ϑ = [τν − ] cs Δtα = [τc − ] cs Δt (6)

In these relations, cₛ is the speed of sound and equals c/√3. In the present simulation, the Boussinesq approximation is used to model the buoyancy force. The external force F is defined as follows:

Fi = 3ωigyβΔT (7)

In this relation, β is the thermal expansion coefficient, gy is the gravitational acceleration, and ΔT is the temperature difference. Finally, the macroscopic quantities of density, velocity, and temperature are calculated using the following relations:

In the LES model, the main objective is to calculate the turbulent viscosity νt such that the total viscosity

νt

equals αt = pr in this relation, Prt is the turbulent Prandtl number, which is considered to be 0.4. The relevant calculations are performed as follows:

In this relation, C is the Smagorinsky constant, assumed to be 0.1. Δ = √(Δx)2 + (Δy)2is the filter width, and Δx and Δy are the grid spacings in the x and y directions. The term |S̄| is defined as follows:

S ̄αβ = (∂αū β + ∂βū α) (13)

Table 1 shows a comparison between the average Nusselt numbers for different Rayleigh numbers. It is evident that the results are in good agreement with previous studies.

The results obtained are presented in this section. Figs. 3 and 4 show the temperature contours and streamlines for Rayleigh numbers ranging from 106 to 109. The results obtained are in good agreement with the numerical analysis of turbulent flow within a cavity. The results clearly show the development of flow and thermal structures with increasing Rayleigh number.

For Ra = 106, the flow field is dominated by a single clockwise rotating cell that covers most of the enclosure. The contours of the stream function are smooth and well-organized, indicating a stable laminar flow field. As such, the isotherms are thicker and gently curved, indicating that heat transfer is dominated by conduction. For Ra = 107, the buoyancy-induced flow field is much stronger than at Ra = 106. As such, the flow field has developed secondary vortical structures at the corners of the enclosure. The streamlines are compressed at the vertical walls of the enclosure, indicating higher velocity gradients at these walls. The isotherms are significantly distorted, especially at the hot and cold walls, indicating the development of thermal boundary layers and a transition towards convection-dominated heat transfer.

At Ra = 108, there is significant unsteadiness and complexity in the flow. There are multiple vortices and distorted streamline patterns, especially in the upper and lower areas of the cavity. The isotherms are nearly horizontal in the core region, indicating strong thermal mixing and a good stratification of the temperature field. The thermal boundary layers near the vertical walls are also reduced, suggesting that the heat transfer by convection is increased.

In the case of the highest Rayleigh number, i.e., Ra = 109, the flow field is complex and turbulent. There are multiple interacting vortices of different scales, especially in the upper part of the enclosure, according to the contours of the stream function. This suggests that there is significant buoyancy-driven instability in this flow field. The isotherms show significant distortion and clustering near the walls, while they remain nearly uniform in the core region. This suggests that convection is dominant over conduction, and there is a significant enhancement of heat transfer.

Thus, an increase in the value of the Rayleigh number causes a significant transition from laminar to turbulent natural convection with smaller thermal boundary layers, more complex flow structures, and an appreciable improvement in the performance of natural convection. The above results are in accordance with the numerical and experimental studies of turbulent natural convection in closed enclosures.

Furthermore, Variation of the dimensionless temperature along the horizontal direction for three different values of the Rayleigh number is shown in Fig. 5. It is clear that near the vertical walls on the left and right sides, there are significant temperature gradients, which confirm the presence of boundary layers. Also, with an increase in the value of the Rayleigh number, the boundary layer thickness decreases, showing an improvement in the rate of convective heat transfer.

In the central region of the enclosure, the temperature profile seems almost uniform for all values of the Rayleigh number, which confirms the well-mixed thermal field due to the presence of a strong convective heat transfer mechanism. However, with an increase in the value of the Rayleigh number, the dimensionless temperature in the central region of the enclosure decreases, which might be due to the improvement in the mixing of fluids due to an increase in the value of the Rayleigh number. From the above results, it might be concluded that with an increase in the value of the Rayleigh number, the rate of convective heat transfer improves, showing sharper temperature gradients near the vertical walls and a uniform temperature profile in the central region of the enclosure.

The numerical analysis of turbulent natural convection within a square enclosure subjected to differentially heated vertical walls has been carried out using the LES model and the LBM. The analysis was performed by considering different Rayleigh numbers from 106 to 109. The numerical results show that the buoyancy-induced flow increases with the increase of the Rayleigh number, resulting in the reduction of the thermal boundary layers and the development of more complex flow structures. The heat transfer within the core of the enclosure is found to be mainly conductive at lower Rayleigh numbers, whereas the role of the convective mechanism increases as the Rayleigh numbers increase. At Ra = 109, the flow structure shows characteristics of turbulent flow, including distorted isotherms and the presence of multiple, interacting vortices. The uniformity of the temperature field increases within the core of the enclosure as the Rayleigh numbers increase, which indicates the enhancement of the heat transfer due to the turbulent flow. The numerical results show good agreement with the benchmark solutions, which indicates the accuracy and reliability of the LES-LBM method.

The authors declare no conflict of interest.

All data are included in the text.

This research received no specific grant from funding agencies in the public, commercial, or not-for-profit sectors.