Hybrid LES/URANS Simulation of Rayleigh-Bénard Convection Using BEM

In this paper, we develop and test a unified hybrid LES/URANS turbulence model with two different Large Eddy Simulation (LES) turbulence models. The numerical algorithm is based on the Boundary Element Method. In the existing hybrid LES/URANS turbulence model we implemented a new Smagorinsky LES turbulence model. The hybrid LES/URANS turbulence model is unified, which means that the LES/URANS interface is changed dynamically during simulation using a physical quantity. In order to define the interface between LES and unsteady Reynolds Averaged Navier Stokes (URANS) zones during the simulation, we use the Reynolds number based on turbulent kinetic energy as a switching criterion. This means that the flow characteristics define where the sub-grid scale or URANS effective viscosity and thermal conductivity are used in the governing equations in the next time step. In unified hybrid turbulence models, only one set of governing equations is used for LES and URANS regions. The developed hybrid LES/URANS model was tested on non-isothermal, unsteady and turbulent Rayleigh-Bénard Convection and compared with an existing model, where LES is based on turbulent kinetic energy. The hybrid turbulence model was implemented within a numerical algorithm based on the Boundary-Domain Integral Method, where a single domain and sub-domain approaches were used. The numerical algorithm uses governing equations written in a velocity-vorticity form. The false transient time scheme is used for the kinematics equation.


Introduction
Heat transfer via natural convection is used widely in the area of HVAC, electronics' cooling, convection in solar panels, process and power engineering. One of the typical natural convection problems in engineering is the Rayleigh-Bénard Convection (RBC), where the fluid develops a regular pattern of convection cells. These occur in a plane horizontal layer of fluid heated from below. In industrial applications, RBC occurs in parallel piped shaped mobile tanks, which are filled with hydraulic oil, air, or other liquid [Ayed, Živković and Tomić (2017)]. Furthermore, RBC occurs in the rotating cavities in the compressor rotor of a gas turbine [Owen (2010)]. RBC represents a heat transfer mechanism, which is achieved by vortical turbulent flow developed by thermal buoyancy.
Numerical analysis are widely used in research and development in several industrial branches.
Employing CFD numerical algorithm Wang et al. [Wang, Yang and Huang (2019)] successfully establish diffusion model of chaff cloud under airflow. Aerodynamic coefficients and factors calculated using different turbulence models (Shear Stress Transport (SST), Reynolds Stress Model (RSM) and k − ) were compared with wind tunnel results.
Tian et al. [Tian, Liu, Wang et al. (2019)] performed a study of underwater explosion bubble using Boundary integral equation (BIE) based 2D Numerical model. He et al. [He and Ma (2019)] developed eigenstrain BIE model for solids with fluid-filled pores analysis.
Fluid and heat transfer engineering problems are, nowadays, solely, or in combination with experiments, solved using Computational Fluid Dynamics (CFD) programs employing different turbulence models. Despite disadvantages in comparison with Direct Numerical Simulation (DNS) and Large Eddy Simulation (LES), the most widely spread turbulence models in the engineering industry are steady RANS and unsteady URANS respectively. From the practical point of view, they are still the most robust and versatile, and provide assessable computational times [Fröhlich and von Terzi (2008)]. In steady RANS or unsteady URANS, the Reynolds Averaged Navier Stokes equations are used for the simulation of the fluid flow. In comparison with other approaches of turbulence modelling LES and DNS, the main advantage of (U)RANS modelling is the attached and fully developed boundary layer. In the near wall region, where, due to small Kolmogorov length, LES requires a very fine mesh, (U)RANS requires an affordable mesh density.
When solving near-wall flows by LES the smallest flow structures must be resolved. The main restriction of the LES is the prediction of the near-wall flow due to the requirement of a very fine grid resolution in the near-wall region. On the other hand, when dealing with flows where large turbulent scales are dominant and the mesh resolution requirements are not too restrictive, e.q. the bluff body flows, LES produces excellent results. URANS faces difficulties when the bulk flow structure is strongly affected by the dynamic of large-scale turbulent eddies. URANS simulations, based on statistical averaging, average the effect of the unsteady large scales and thus lose the information of the flow.
Comparing both turbulence modelling approaches in the bulk flow, the main advantage of the LES is that it solves large-scale eddies directly, while (U)RANS provides only statistical information of the turbulent flow, and, thus, does not provide information of structures, frequencies, etc. Thus, one of the main reasons for hybrid LES/URANS development was to overcome those limitations of both turbulence models.
The basic idea of a hybrid LES/URANS turbulence model is to combine the advantages of both LES and URANS turbulence models. In the near wall region, where LES requires dense meshes, the URANS model is used, and in the bulk flow, where the URANS model lacks information about the flow, the LES model is used [Spalart and Shur (1997)]. That means that one model considers the attached flow in the boundary layer using URANS, and other resolves detached eddies in the bulk flow using LES.
We can divide LES/URANS hybrid turbulence models based on the fact of how they divide both regions in the computational domain. One approach is represented by segregated hybrid models, where the computational domain is divided prior to simulation [Xiao and Jenny (2012)]. In segregated hybrid models, LES and URANS are solved separately. The other approach is unified hybrid models, where the same set of governing equations is used for fluid flow calculation by LES and URANS. Besides unified hybrid models with a hard interface, where the LES/URANS interface is defined geometrically, we know models with dynamical interfaces, where the switching criterion, based on physical quantity, defines the LES/URANS interface dynamically during the simulation [Breuer, Aybay and Jaffrezic (2010)]. Hybrid models with a dynamic interface depend on flow conditions. The advantage of such hybrid models is that they are adapted to the flow structure in every time step and vortex position. In this manner we achieve a smoother transition between the URANS and LES regions. After the basic concept of hybrid models was known, there were still open topics, like coupling techniques between regions, which quantity is the most suitable for the switching criterion, and also which LES and URANS are the most suitable for particular application of fluid flow.
Ravnik et al. [Ravnik, Škerget and Hriberšek (2006)] employed a numerical solver with a combination of a BEM based solution of the kinematics equation and FEM-based solution of the kinetics equation. The numerical algorithm with the LES model was tested on non-isothermal turbulent simulation of natural convection. Thermal turbulent flow employing BEM was also studied further by Ramšak et al. [Ramšak and Škerget (2008)]. A two-equation k − turbulence model was used in the numerical algorithm with the multidomain Boundary Element Method using mixed boundary elements and a subdomain technique. Lupše et al. [Lupše, Škerget and Ravnik (2014)] compared Spalart-Allmaras, Chien and Abe-Kondoh-Nagano turbulence models employing a BEM based algorithm and velocity-vorticity form of governing equation. Simulations were performed on a turbulent channel flow and backward facing step.
The novelty of this work is a new Smagorinsky LES turbulence model implemented in a verified hybrid LES/URANS turbulence model with dynamic LES/URANS interface. We developed the hybrid LES/URANS turbulence model and used a numerical algorithm based on BEM using a combination of single-domain and sub-domain approaches. We employed a hybrid LES/URANS turbulence model using BEM and a combination of single-domain and sub-domain approaches. The LES/URANS hybrid model was used for non-stationary non-isothermal turbulent flow simulation. The main characteristic of a unified hybrid model is that LES and URANS models are unified within the transport equation for turbulent kinetic energy k, and are being solved simultaneously. Using a LES/URANS hybrid model means that the LES and URANS regions are not defined before simulation, but they are defined dynamically during simulation itself. During the numerical simulation, the algorithm switches automatically to URANS within the near wall region, and to the LES in the bulk flow. Sub-grid viscosity is used in LES mode, while, URANS modelled viscosity in URANS mode.
A 2D benchmark simulation of Rayleigh-Bénard Convection in a rectangular cavity with length/height ratio 4/1 was used for testing the present BEM implementation of the hybrid model. We chose the benchmark simulation of Rayleigh-Bénard Convection because it requires a reasonable increase of computational resources in comparison with the previous study of natural convection in a square cavity. The benchmark case has a simple geometry, well defined boundary conditions, and the same time features complex turbulent flow with more main big vortexes, small eddies, thin boundary layer and challenging high temperature gradients.
In the past studies, various experimental and numerical approaches were used to examine the Rayleigh-Bénard Convention for a wide range of Rayleigh numbers and enclosure dimension ratios.
Kenjereš et al. [Kenjereš and Hanjalić (2000)] studied two-dimensional Rayleigh-Bénard Convection employing an algebraic model for turbulent heat flux. They studied the existence, creation, and behaviour of convective rolls. A wide range of Rayleigh numbers was performed, from Ra = 10 5 up to Ra = 10 12 , and enclosure aspect ratios from 4 : 1 up to 32 : 1. A comparison was made for full three-dimensional direct numerical DNS simulations, large eddy simulations, and three-dimensional transient Reynolds-averaged Navier-Stokes approaches. They have shown that the 2D approach is an acceptable method for average analysis of fully 3D flows with at least one homogeneous direction. A numerical study of Rayleigh-Bénard Convection was recently performed by Dabbagh et al. [Dabbagh, Trias, Gorobets et al. (2016a]. Dabbagh et al. [Dabbagh, Trias, Gorobets et al. (2016a)] employed Direct Numerical Simulation to study small scale motions of Rayleigh-Bénard Convection. The fluid used in the simulation was air. Furthermore, Dabbagh et al. [Dabbagh, Trias, Gorobets et al. (2016a)], used a new subgrid-scale model S3PQR for simulation of Rayleigh-Bénard Convection. The developed SGS model was tested on Rayleigh number Ra = 10 8 , and came in good agreement with the prediction of turbulent kinetics, although it showed higher heat flux within the boundary layer.
Rayleigh-Bénard Convection was also studied experimentally. Chu et al. [Chu and Goldstein (1973)] first found the existence of convection roll cells. They also found that the thermal vortexes were released periodically from relatively fixed locations, which becomes even more stable when the Ra number is increasing. Chu and Goldstein concluded their investigation with a Nusselt and Rayleigh number correlation with power law: N u = 0.123 · Ra 0.294 .
Although vast numerical and experimental research has been done, there is still some discrepancy between heat transfer scaling. For full understanding of the details and an overview of the whole physical phenomena of Rayleigh-Bénard Convection, DNS is still the most widely used.

Velocity-vorticity formulation
In the present work, we used a unified hybrid turbulence model where only one set of governing equation is solved. The velocity-vorticity formulation is used for governing equations, and a false transient time scheme for the kinematics equation. The fluid used in the simulation is considered as incompressible, with constant material properties: Molecular viscosity ν m and heat diffusivity a m . The turbulent quantities are calculated by URANS and SGS models and further, according to the switching criterion used in the governing equations for kinetics and kinematics. The governing equations have to be averaged for URANS and filtered for LES simulation. In time averaged/filtered governing equations, we use effective kinematic viscosity ν ef and effective heat diffusivity a ef . Both quantities are the sum of the molecular and turbulent parts. Effective kinematic viscosity is the sum of molecular and turbulent kinematic viscosity defined as ν ef = ν m + ν t , and effective heat diffusivity is the sum of molecular and turbulent heat diffusivity defined as The kinematics equation is the curl of the vorticity vector, written as a pseudo transient problem that allows an under-relaxation approach. The kinetics equation is the vorticity equation obtained using the curl of the momentum equation. Employing false diffusivity α, the kinematic equation for two-dimensional incompressible flow in the velocity-vorticity formulation is written as Škerget et al. [Škerget, Hriberšek and Kuhn (1999)]: The kinetics equation is the vorticity equation obtained using the curl of the momentum equation. The only acting force on the fluid in natural convection driven flows is buoyancy resulting from density differences due to the temperature. For simulating buoyancy in incompressible fluid, we introduce Boussinesq's approximation in the following kinetic equation. For solving non-isothermal turbulent flow, we also have to introduce the correlation between turbulent viscosity ν t and turbulent thermal diffusivity a t ; we used correlation using turbulent Prandtl number P r t defined as: a t = νt P rt . With the curl of the second extended form of the momentum conservation equation, the equation for kinetics can be written as: where the parameters of the equation are the reference temperature T 0 , thermal volumetric expansion coefficient β T , gravity force g and the non-zero component of the vorticity vector ω. The term f m i in the Eq.
(2) is defined as: For the simulation of the non-isothermal turbulent flow we have to write the energy equation for temperature T [Lupše, Škerget and Ravnik (2014) and the transport equation for turbulent kinetic energy k, written as Breuer et al. [Breuer, Jaffrezic and Arora (2008)]: In the Eq. (5), the P is the production part, D the dissipation part, its sum is the source term I, ν m is molecular and ν t is turbulent viscosity.
For the fluid flow simulation in the present research, we used a hybrid turbulence model based on the transport equation for turbulent kinetic energy (5). The main characteristic is that, in the URANS region, turbulent viscosity ν t and dissipation part D, calculated using a URANS model, ν t,U RAN S and D U RAN S are used, and in the LES region, turbulent viscosity and the dissipation part, calculated using SGS model ν t,SGS and D SGS .

Smagorinsky LES model
For LES subgrid-scale viscosity ν t,sgs and D sgs we used Smagorinsky LES model [Pope (2000)]. The filter width is calculated as: Subgrid-scale viscosity ν t,sgs is calculated as where l S is the Smagorinsky length scale, C S the Smagorinsky coefficient andS is defined as: The dissipation part of the Smagorinsky model is defined as:

Turbulence kinetic energy based turbulence model for LES
The LES subgrid-scale viscosity used in the turbulence model, based on the transport-equation for turbulent kinetic energy k, is calculated as Pope [Pope (2000)]. Filter width ∆, dissipation part D sgs and turbulent viscosity ν t,sgs for the sub-grid model can be written as: where C d = 1.0, C µ = 0.05 [Schumann (1975)] .

Turbulence model for URANS
For the URANS part of the hybrid model, we used the one-equation turbulence model based on the transport equation for turbulent kinetic energy [Breuer, Jaffrezic and Arora (2008)]. Turbulent viscosity ν t,U RAN S is defined as: where l µ,v is the characteristic length defined as l µ,v = C l,µ · y, C l,µ = 0.33 and v 2 normal velocity fluctuations defined as: Dissipation part D U RAN S is defined empirically as: is dissipation length defined as: l D,v = 1.3 · y/ 1 + 2.12 ν y(v 2 ) 1/2 . For the switching criterion we calculate Reynolds number defined turbulent kinetic energy k:

Switching criterion between LES/URANS region
For a dynamical definition of the LES and URANS area in the computational domain during the calculation, we used the switching criterion based on the physical quantity Reynolds number defined by turbulent kinetic energy k; Re k . Re k can be written as: where k is the turbulent kinetic energy, ν the viscosity of the fluid, and y the normal distance from the wall.
Thus, the switching criterion between the LES and URANS areas can be defined as Breuer et al. [Breuer, Jaffrezic and Arora (2008)]: Re k ≤ C switch → RANS mode,

Summary
Characteristics and labels of each LES/URANS turbulence model are summarised in the Tab. 1. In models 1 and 2, the same URANS model and switching criterion were used, but

Numerical method
In the present study of fluid flow, we used the subdomain Boundary Element Method based numerical algorithm. The Navier-Stokes equations were written in a velocity-vorticity formulation [Kocutar, Škerget and Ravnik (2015)]. For the BEM formulation, the effective viscosity is divided into constant and variable parts ν m + ν t = ν ef = ν 0 + ν, in the same manner, we divide thermal diffusivity a m + a t = a ef = a 0 + a.

Integral formulation of parabolic-diffusive fundamental solution
Employing a parabolic-diffusion fundamental solution where (ξ, t F ) is a source point, (s, t) a reference point within a domain, r a vector from source point ξ to reference point s, and τ a time step, the kinetic Eq.
(2) for two-dimensional flow in integral form can be written as: In the same way, we can write the kinematics Eq. (1). For further simplification of the numerical scheme, we select the under relaxation term of false transient time scheme α equal to the viscosity term ν 0 and write: Considering the Eq. (18) and the kinematic Eq.
(1) is also used to determine boundary vorticity value.
Before implementing the governing transport equations (velocity, vorticity, temperature, and turbulence kinetic energy) into the numerical algorithm, they were written in the discrete form. A collocation strategy was employed to derive the final system of linear equations.

Discrete form
Before implementing the governing transport equations (velocity, vorticity, temperature and turbulence kinetic energy) into the numerical algorithm they have to be written in the discrete form.
The computation domain consists of C internal cells and E boundary elements, thus we can write: Γ = E e=1 Γ e , Ω = C c=1 Ω c . For calculation of the sum of integrals over boundary domain we interpolate field function over every boundary element E using Φ and every computational cell using φ. In the numerical algorithm, we used quadratic interpolation function for computational nodes and constant function for fluxes.
In the following Eqs. (21) to (23) the number of nodes in each internal cell or boundary element is noted with index p and the degree of time polynomial ψ with index m. Employing parabolic-diffusion fundamental solution u * and time step F the equation for kinematics (20) for two-dimensional flow can be written as: where transposition is labeled with T and vectors of nodal values are labeled with {}.
In the same matter the equation for kinetics (17) can be written as: For solving non-isothermal fluid flow the energy Eq. (4) has to be written as: and for the purpose of modeling the turbulent fluid flow the transport equation for turbulent kinetic energy k (5) as: For solving the system of equations written in matrix notation we employed the Sleijpen et al. [Sleijpen and Fokkema (1993)].

Numerical algorithm
The dual reciprocity method is usually applied when dealing with non-linear problems such as fluid flow, to remove the need for performing domain integration. In our work, we use a domain decomposition approach instead. We mesh the entire domain using a standard meshing tools steaming from computational fluid dynamics packages and preformed domain decomposition in such a way that each domain element represents a subbomain. Within each subdomain we apply BEM for produce linear equations for each of the function and flux nodes within the subdomain. Next we use continuity of function and compatibility of flux for join subdomain equations into a large system of liner equations. The resulting system is sparse, and thus enables the use of dense meshes. This domain-decomposition approach is used to discretizise transport equations (velocity, vorticity, temperature, and turbulence kinetic energy). The determination of boundary vorticity values (20) is done using standard BEM with calculation of the domain integrals. Internal mesh used in simulation is refined in y-direction towards the upper and lower wall using a geometrical series with the ratio between longest and shortest elements R = 16.
Since our method is based on domain decomposition, we can use standard CFD meshing tools to produce a domain mesh, which is capable of capturing the small flow structures required by the LES. We are proposing a hybrid RANS-LES approach, since the LES demands a very fine grid close to the walls, which is currently too costly of our BEM based algorithm. Thus, using the switching criterion we are able to use RANS in parts of the domain, where the grid is fine enough for LES.
For the transient simulation, the numerical algorithm starts with a time loop, and continues with a global iteration loop. Then, the algorithm starts a kinematics loop, where the boundary vorticity ω and then internal velocity v are calculated. When the convergence for ω and v is fulfilled, the algorithm continues on to the next step. In the next step, the energy equation for kinetics is solved, and, finally, the vorticity ω for kinetics. Handling turbulence starts with the outer-loop for effective viscosity. In the LES/URANS hybrid turbulence model, the turbulence viscosity and dissipation part are calculated for sub-grid ν sgs , D sgs , Sections 3.1.1 and 3.1.2, and for URANS model ν U RAN S , D U RAN S , Sections 3.2. Calculating the Re k 3.3 the interface is defined between the LES and URANS regions. By defining the interface, we choose which D and η ef are used for the turbulent kinetic energy k equation, and, in the following time step, for solving the governing equations for the fluid flow. After that, the transport equation for turbulent kinetic energy k is solved. After the convergence for ν ef is fulfilled, the algorithm checks the convergence for vorticity ω and temperature T .

Simulation of natural Rayleigh-Bénard Convection within a square cavity
The developed turbulence model was tested on an unsteady turbulent Rayleigh-Bénard Convection with the ratio of height and length of H : L = 1 : 4. The Rayleigh-Bénard Convection is a buoyancy-driven flow, where the temperature difference between the heated lower plate and cooled upper plate acts as a major generator for heat and fluid momentum transfer. When the Rayleigh number is increased the flow becomes more turbulent, and the arrays of roll cells are becoming more unsteady and chaotic. Due to buoyancy, the flow is rising from the lower hot plate and penetrating the bulk flow up to the upper cold plate, where flow mixing leads to intensive vortexes and an almost uni-thermal middle area. Rayleigh-Bénard Convection within the cavity is steady up to a Rayleigh number around Ra = 10 7 .
The boundary conditions for Rayleigh-Bénard Convection are defined in Fig. 1 in a non-dimensional setting: Temperature T h = 1 on the uniformly heated lower wall, temperature T c = 0 on the uniformly cooled upper wall, adiabatic boundary condition q = 0 on the left and right walls, no-slip velocity ∂ v ∂n = 0 boundary conditions on all four walls, and gravity g in the opposite direction of coordinate y. The problem is defined by the Rayleigh number. The Rayleigh number and the Prandtl number, which we choose to be the Prandtl number of air P r = 0.71, are defined as The cavity of Rayleigh-Bénard Convection is assumed to be L = 0.8m wide and H = 0.2m high. The fluid in the cavity is air with kinematic viscosity ν = 17.6 · 10 −6 m 2 /s, thermal volumetric expansion coefficients β = 1/313.15K −1 , and thermal conductivity λ = 0.0265W/mK. That means that the difference between heated and cooled walls is ∆T = 174.08K at Ra = 10 8 .
For numerical simulations, we used a computational mesh with 160 × 80 domain elements and 321 × 161 nodes, respectively. The mesh was compressed in the normal wall direction in the y coordinate, with the geometrical series having the ratio between the longest and shortest elements R = 16. The results of Rayleigh-Bénard Convection simulation are evaluated using the averaged Nusselt number. The Nusselt number measures the heat flux through a solid wall, and can be defined as: 6 Results 6.1 Direct Numerical Simulation of Rayleigh-Bénard Convection at Ra = 10 4 − 10 7 Prior to using the LES/URANS hybrid model for turbulent flow, the domain mesh and numerical algorithm were validated using a steady-state direct numerical simulation of Rayleigh-Bénard Convection for Rayleigh numbers from Ra = 10 4 to Ra = 10 7 . We tested our numerical algorithm on three different meshes, 140 × 70R16, 160 × 80R16, and 180 × 90R16, which means 281 × 141, 321 × 161 and 361 × 181 nodes. The results for Rayleigh numbers from Ra = 10 4 to Ra = 10 7 were compared with the results of other authors, and are shown in the Tab. 2: Kenjereš et al. [Kenjereš and Hanjalić (2000)] and Peng et al. [Peng, Hanjalić and Davidson (2006)]. The results with all meshes 140×70R16, 160×80R16 and 180×90R16 show good agreement with benchmark results. The computational mesh 140 × 70R16 was too coarse for turbulent flow at Ra = 10 8 due to a thin boundary layer and high temperature gradient. Since the values of Nusselt number N u/L compare well with published values and since only small differences were observed between different meshes, we concluded that the numerical algorithm is appropriate for further simulations using the Hybrid LES/URANS turbulence model using the 321 × 161 nodes mesh. Fig. 2 shows the temperature field of laminar Rayleigh-Bénard Convection for Rayleigh numbers from Ra = 10 4 up to Ra = 10 7 . We found that, with the increasing of the Rayleigh number, the temperature gradient near the upper and lower wall was increasing, and the area of unithermal bulk flow (green color) was getting nearer to the horizontal walls. In the Ra = 10 7 case, we can see perturbations in the temperature peaks, which indicate the transition to unsteady turbulent flow. Fig. 3 shows the temperature profiles of Rayleigh-Bénard Convection using computational mesh 160x80R16 for the laminar steady-state simulations for Rayleigh numbers from Ra = 10 4 to Ra = 10 7 at height y = L/2 and length x = L/2 respectively.

Turbulent Rayleigh-Bénard Convection at Ra = 10 8
Simulation of turbulent Rayleigh-Bénard Convection at Ra = 10 8 was performed using the computational mesh 160x80R16 and by using the hybrid LES/URANS turbulence model. For the switching criterion between LES and URANS area we defined C switch = 20. and chaotic. From the streamlines on the temperature field, Figs. 6 and 7 above, we noticed that the vortexes occur in the corners and near connection lines (0.25, 0.5 and 0.75 · L). As mentioned before, the driving force of Rayleigh-Bénard Convection is the difference in temperature between the upper and lower walls, which, due too gravity and different density of the fluid, generates buoyancy. The fluid is accelerating in two areas, partly along the heated and cooler walls in the x-direction, and mainly in the y-direction near the walls and connection lines (0.25, 0.5 and 0.75·L), where the vortexes collide in the vertical walls and each other respectively. Also, the majority of turbulent eddies occur in areas where flows collide in the wall or in each other. From the temperature and vorticity fields, Figs. 6 and 7 above and middle, it can be seen that the boundary layer near the horizontal is very thin and temperature gradients are high. The bulk flow between the plates is almost unithermal. Due to a thin boundary layer and high temperature gradients the mesh should be compressed near the walls. When comparing the Hybrid LES/URANS region with the temperature and vorticity fields, we noticed that the URANS is being used near the fix-wall regions where attached flow occurs, and LES in the bulk where eddies detach.

Figs. 6 and 7 show that Rayleigh-Bénard Convection is already turbulent, thus unsteady
From temperature-vorticity phase portraits for Ra = 10 8 for point A (0.05, 0.0498), Fig. 4, we see that the flow is unsteady and chaotic, thus turbulent. At Rayleigh-Bénard Convection the transition from laminar to turbulent flow occurs between Ra = 10 7 −10 8 . The unsteady and turbulent characteristics of Rayleigh-Bénard Convection can be seen in Figs. 6 and 7, where the temperature and the vorticity fields for Ra = 10 8 are shown. Fig. 5 shows the temperature time series power spectra at a point 0.050, 0.4983 for both hybrid turbulence models. Power spectra show the presence of eddies of the different sizes and frequencies of its forming oscillations. The power spectra were calculated using the Fast Fourier Transform. Fig. 8 shows the Nusselt number N u/L at Ra = 10 8 for both hybrid turbulence models 1 and 2 are shown. The ratio of the area where the LES turbulence model is used is shown on the Fig. 9.
According to the Kolmogorov hypothesis, both turbulence models were compared for a one-dimensional longitudinal energy spectrum of turbulence. Fig. 10 shows that both turbulence models start with increasing of E(k) in the energy containing range of eddies and continue with inertial subrange, where bigger eddies forward the energy to smaller scale eddies, and end up with a dissipation subrange. It can be seen that the spectrum for both models in the inertial subrange is parallel to the −5/3 slope used in the Kolmogorov hypothesis.
Tab. 3 shows the averaged Nusselt number Nu/L of the present LES/URANS hybrid turbulence models 1 and 2 and compared the results of other researchers. Kenjereš et al. [Kenjereš and Hanjalić (2000)] performed numerical simulation based on URANS, while Dabbagh et al. performed simulations using DNS [Dabbagh, Trias, Gorobets et al. (2016b)] and LES [Dabbagh, Trias, Gorobets et al. (2016a)], Peng et al. used LES [Peng, Hanjalić and Davidson (2006)] and Choi et al. elliptic blending second-moment closure [Choi and Kim (2008)]. Grossmann et al. [Grossmann and Lohse (2000)] derived a scaling theory based on experimental data, and proposed correlation between Rayleigh and Nusselt numbers for Rayleigh-Bénard Convection defined as Nu = 0 .27Ra 1/4 + 0.038Ra 1/3 . For the further evaluation of both turbulence hybrid models and the influence of implemented Smagorinsky LES model into Hybrid model 2, the results of Hybrid models Comparing the results for both Hybrid models we noticed that the new Hybrid model 2 comes in better agreement with the average value of Numerical simulations and also with the average value of LES simulations.

Conclusions
In the present work, we used the developed hybrid LES/URANS turbulence models successfully using two different approaches for LES for solving unsteady non-isothermal turbulent Rayleigh-Bénard Convection. We employed a BEM based numerical algorithm in combination with the velocity-vorticity formulation of governing equations.
In the LES/URANS models we used two different approaches for the LES model, a one-equation LES model based on the transport equation for turbulent kinetic energy k and basic Smagorinsky LES model, calculated from the current velocity field. The numerical algorithm was tested using air Rayleigh-Bénard Convection with the ratio length/height 4 : 1. For the steady flow up to Rayleigh numbers Ra = 10 7 , we used direct numerical simulation. The transition from laminar to turbulent flow occurs at Rayleigh numbers Ra = 10 7 − 10 8 , where flow becomes unsteady and turbulent. For developed turbulent flow above Ra = 10 8 , we used a unified one-equation LES/URANS hybrid turbulence model. In the LES/URANS hybrid models we employed two different LES models, a one-equation LES model based on the transport equation for turbulent kinetic energy k and a Smagorinsky LES model. Both turbulence models use the same one-equation k RANS model and same switching criterion, and Reynolds number based on turbulent kinetic energy Re tot . Results from the hybrid models using both LES approaches came in with good agreement with results from Kenjereš at al. [Kenjereš and Hanjalić (2000)], Peng et al. [Peng, Hanjalić and Davidson (2006)], Choi et al. [Choi and Kim (2008)], Dabbagh et al. [Dabbagh, Trias, Gorobets et al. (2016b,a)] and Grossmann et al. [Grossmann and Lohse (2000)].
Funding Statement: The authors acknowledge the financial support from the Slovenian Research Agency (research core funding No. P2-0196).

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.