Tian Qin, Yingying Luo, Jiemin Zhan, Wenqing Hu^*, Yutian Li, Ziwei Huang
Department of Applied Mechanics and Engineering, School of Aeronautics and Astronautics, Shenzhen Campus of Sun Yat-sen University, Shenzhen, China.
DOI: 10.4236/ojfd.2023.133013   PDF    HTML   XML   114 Downloads   558 Views   Citations

Abstract

This paper conducts a Large Eddy Simulation (LES) of Rayleigh Bénard convection in a cubic cavity based on the WMLES S-Omega subgrid-scale model. For a cubic cavity with a vertical temperature difference of 6.7°C and 20°C, the velocity pulsation profiles and the mean velocity profiles of the vertical section in the middle of the cubic cavity were simulated, respectively. And they are consistent with the experiment results. Furthermore, the mean velocity field of the vertical cross-section in the middle of the cavity was calculated. Structures of the mean velocity field in the two cases are similar. A counterclockwise large vortex is found to occupy the cavity, and there are two small clockwise vortices in the lower left and upper right corners, and the mean velocity fields at two different temperature differences are consistent with the experimental results. The two-dimensional instantaneous temperature field and mean temperature field with different cross-sections in the z-direction, as well as the three-dimensional instantaneous isothermal surface structure, indicate that the large-scale circulation motion within the cubic cavity is moving diagonally. In addition, the structure of the mean streamline also illustrates this viewpoint. For the reverse vortex formed at two corners in the mean streamline structure, we used the Q criterion to identify and obtain two vortex structures similar to boomerangs. The basic turbulent structure in RB thermal convection includes the rising and falling plumes generated by buoyancy effects.

Keywords

Large Eddy Simulation, Rayleigh-Bénard Convection, Cubic Cavity, Plume Flow

Share and Cite:

1. Introduction

Rayleigh-Bénard convection (RBC) is a thermal convection movement occurring in the horizontal fluid layer. In the atmosphere and oceans, thermal convection has a significant impact on global climate; The thermal convection in the Earth’s outer core is considered to be the origin of the Earth’s magnetic field; Besides, heat convection is widely related to industries, such as cooling and heating of air conditions, heat distribution of electronic components, and heat distribution and cooling inside nuclear reactors [1] . According to relevant research, when RBC is above Rayleigh number Ra ≈ 1708 [2] in the cubic cavity, the fluid flows due to buoyancy generated by the part with low density. The flow inside RBC first generates cold and hot plumes and mixes, and its shape is very close to the Prandtl number, Rayleigh number, and initial conditions of RBC. When the Rayleigh number reaches a certain value, RBC will change from stable laminar flow to highly nonlinear turbulent motion.

RBC in a cube cavity is widely used for Computational Fluid Dynamics (CFD) code verification due to its geometric simplicity. The numerical simulation of natural convection has been extensively studied and analyzed in laminar and turbulent states. Many numerical studies have focused on convection in rectangular cavities using horizontal temperature differences (sidewall heating and cooling). The initial numerical research on RBC was conducted through two-dimensional numerical simulations, which were used as a benchmark for evaluating numerical solutions [3] [4] [5] . The experimental implementation of two-dimensional turbulent convection is a complex task because two-dimensional turbulence is more susceptible to perturbation than three-dimensional turbulence [6] . When the Rayleigh number is 1.58 × 10⁹, the convection experiment in the cubic cavity with the aspect ratio of 0.5 is taken as the experimental benchmark for two-dimensional convection in the cubic cavity with an isothermal vertical boundary [7] . In a high rectangular cavity filled with air with an aspect ratio of 0.03 [8] , a nearly two-dimensional convective simulation was achieved. In the air convection in a cavity with an isothermal vertical wall and a length-width ratio of 0.32, numerical studies (DNS and LES) and experimental studies were carried out at the Rayleigh number Ra = 1.58 × 10⁹ [9] . By comparing the simulation results of DNS and LES, it is found that there is a good consistency and consistent with the experimental results. The convection in the cube cavity with vertical temperature difference is characterized by a variety of states. Extensive research on the medium value of Rayleigh number Ra < 10⁸ [10] [11] [12] [13] [14] has found that, even when the Rayleigh number is reduced to close to the critical number, the three-dimensional structure of convection under vertical temperature difference is very similar to the three-dimensional characteristics of high Rayleigh number. For the convection study of cubic cavity with vertical temperature difference and Rayleigh number Ra > 10⁸, it is found [15] [16] [17] that the internal movement of the cavity under a high Rayleigh number is very complex and is characterized by various large-scale flows.

Based on the experimental model of Vasiliev et al. [18] , the mean velocity profiles, velocity pulsation profiles, and mean velocity field of the middle vertical section obtained by numerical simulation were compared with the experiment in this paper. On this basis, we studied the relevant two-dimensional cross-section and three-dimensional structures in the cubic cavity and determined the motion characteristics of the large-scale circulation.

2. Methods

The flow in a cubic cavity is described by the Boussinesq equations

$\nabla \cdot u=0$ (1)

$\frac{\partial u}{\partial t}+u\cdot \nabla u=-\nabla p+\frac{Pr}{\sqrt{Ra}}{\nabla }^{2}u-PrT\frac{g}{\left|g\right|}$ (2)

$\frac{\partial T}{\partial t}+u\cdot \nabla T=\frac{1}{\sqrt{Ra}}{\nabla }^{2}T$ (3)

where, u is the velocity vector, T is temperature, g is gravity, and p is pressure.

The expressions for the Prandtl number Pr and the Rayleigh number Ra are:

$Pr=\frac{\nu }{\chi }$ (4)

$Ra=\frac{g\beta \Delta T{L}^3}{\nu \chi }$ (5)

where, β is the coefficient of thermal expansion, ν is the kinematic viscosity, and χ is the thermal diffusivity.

In turbulence, the root mean square (RMS) of velocity is generally used to represent velocity pulsation, and the calculation formula is as follows:

$u^{\prime }=\sqrt{\frac{1}{N}{\displaystyle \underset{1}{\overset{N}{\sum }}{\left(u_i-u_m\right)}^2}}$ (6)

where $u_i$ is the instantaneous velocity, $u_m$ is the mean velocity.

LES can directly analyze large vortices and model small vortices, while vortices smaller than the grid size will be directly filtered. The filtered control equations are

$\nabla \cdot \bar{u}=0$ (7)

$\frac{\partial \bar{u}}{\partial t}+\bar{u}\cdot \nabla \bar{u}=-\nabla \bar{p}+\frac{Pr}{\sqrt{Ra}}{\nabla }^2\bar{u}-Pr\bar{T}\frac{g}{\left|g\right|}-\nabla \cdot \tau$ (8)

$\frac{\partial \bar{T}}{\partial t}+\bar{u}\cdot \nabla \bar{T}=\frac{1}{\sqrt{Ra}}{\nabla }^2\bar{T}-\nabla \cdot h$ (9)

where $\tau$ and h are sub-lattice stress and sub-lattice heat flux, respectively, expressed as

${\tau }_{ij}=\overline{u_i{u}_j}-{\bar{u}}_i{\bar{u}}_j$ (10)

$h_i=\overline{u_{i}T}-{\bar{u}}_i\bar{T}$ (11)

The subgrid stress and subgrid heat flux construct the filtered control equations into a closed mode, which is the key to realizing LES. The subgrid stress is proportional to the filtered strain rate tensor on the minimum solution size, and the calculation formula is

${\tau }_{ij}=\frac{{\delta }_{ij}{\tau }_{kk}}{3}-2{\nu }_t{\bar{S}}_{ij}$ (12)

${\bar{S}}_{ij}=\frac{1}{2}\left(\frac{\partial {\bar{u}}_i}{\partial {\bar{x}}_j}+\frac{\partial {\bar{u}}_j}{\partial {\bar{x}}_i}\right)$ (13)

where ${\nu }_t$ is the eddy viscosity coefficient.

In WMLES S-Omega subgrid-scale model, ${\nu }_t$ is defined as

${\nu }_t=\left(C_S\Delta \right)\left(S-\Omega \right)\left\{1-\exp \left[-{\left(y^{+}/25\right)}^3\right]\right\}$ (14)

where constant C_S = 0.2, Δ is the subgrid scale, S is the strain rate, Ω is the vortex strength, and y⁺ is the dimensionless distance from the wall surface to the center of the first layer of grid.

3. Numerical Model

3.1. Geometric Model and Boundary

According to the experimental model proposed by A. Vasiliev et al. [18] , experimental studies on turbulent convection were conducted in two essentially identical cubic cavities , namely cubic cavity A and cubic cavity B. The two experimental closed cavities are both cubic cavities with a side length L of 250 mm, and the main difference between the two experimental units (A and B) is related to the material selection of the heat exchanger. The experimental plans are Experiment I and Experiment II, and their basic parameters are shown in Table 1.

The closed cubic cavity structure is shown in Figure 1, with a calculation domain of 250 × 250 × 250 mm³, the lower side plate is hot, and the upper side plate is cold. The upper and lower side plates have a vertical temperature difference, and the remaining four sides are adiabatic wall surfaces. Taking the middle vertical section (y = 125 mm) as the research object. The water in the cubic cavity is regarded as an incompressible, continuous medium with constant physical properties. The proportion of radiation heat transfer is very low, so it can be ignored. The gravitational acceleration direction is the negative z-axis direction. The initial water temperature is 25˚C, the specific heat capacity (C_P) is 4179 J/(kg∙K), the dynamic viscosity (μ) is 9.076 × 10⁻⁴ N∙s/m², and the thermal conductivity is 0.6185 W/(m∙K). According to Formula (4), Pr of the water in the cubic cavity is 6.1. The numerical simulation uses two operating conditions, where the vertical temperature differences are 6.7˚C and 20˚C, respectively, and the corresponding Ra = 2 × 10⁹, 6 × 10⁹. The parameters related to numerical simulation are consistent with the parameter settings in experimental I and II in the reference [18] .

The RBC is solved using the WMLES S-Omega model of LES based on Fluent 15.0 transient solution. Pressure-velocity coupling algorithm of SIMPLE is adopted, the finite central difference scheme is used for momentum, and the second-order upwind difference scheme is used for other equations. During the numerical simulation, the flow field reached dynamic stability after 2000 seconds of computational time. It continued with 2000 seconds of monitoring statistics.

3.2. Grid Settings

In general, accurate numerical results were obtained by the fine enough mesh grids. However, an excessive number of grids will consume a lot of computing resources. Therefore, it is necessary to select an appropriate number of grids for the model by combining the two factors. For LES, based on verifying the grid independence, this paper adopted a uniform grid with a quantity of 10⁶. The grid of the middle vertical section is shown in Figure 2.

4. Results and Analysis

4.1. Verification of Mean Velocity Field

In the case of the vertical temperature difference of 6.7˚C, the numerical result of the mean velocity field in the middle vertical section is shown in Figure 3(a). The overall mean velocity field in the middle vertical section of the cubic cavity is an anticlockwise vortex that occupies the entire cavity. There are two small reverse-rotating eddies in the lower left corner and the upper right corner that are opposite to the mainstream vortex. As the center approaches, the mean velocity value decreases, and the mean velocity increases as the wall approaches. The maximum mean velocity occurred on all four walls. The numerical simulation agrees with the results of Exp. I-A (Figure 3(b)) and Exp. I-B (Figure 3(c)) [18] . In the case of the vertical temperature difference of 20˚C, the numerical result of the mean velocity field in the middle vertical section is shown in Figure 4(a) and consistent with the Exp. II-A (Figure 4(b)) [18] . It is found that when the vertical temperature difference increases by nearly three times, the structural distribution of the mean velocity field in the middle section is similar, and generally presents a structure of a large anticlockwise vortex plus two small inverted vortices at the corners.

4.2. Verification of Mean Velocity Profiles and Velocity Pulsation Profiles

The experimental data from the reference can be used as a benchmark for CFD code, and the accuracy of the selected model can be verified through numerical simulation by setting relevant conditions. Figure 5 shows the comparison between the simulation of mean velocity profiles and Exp. I at the vertical temperature difference of 6.7˚C. Figure 6 shows the comparison between the simulation of velocity pulsation profiles (RMS profiles) and Exp. I at the vertical temperature difference of 6.7˚C. Figure 7 shows the comparison between the simulation

of mean velocity profiles and Exp. II at the vertical temperature difference of 20˚C. In general, the numerical simulation results of velocity profiles and mean velocity profiles are agree with the experimental results.

4.3. Large-Scale Circulation (LSC) Dynamics

Slices of the cubic cavity along the z-direction are used to study the dynamic characteristics of large-scale circulation, with horizontal cross-sections (parallel to the bottom) of z = 10 mm, z = 125 mm, and z = 240 mm, respectively. From Figure 8, it can be observed from the instantaneous temperature contour of the z = 10 mm cross-section that the higher temperature values in both temperature difference conditions are concentrated in the upper right part of the diagonal line passing through points (250, 0) and (0, 250) on the cross-section, indicating

that the fluid near the bottom heating surface first tends to move towards the corner after being heated. The temperature distribution is relatively uniform at z = 125 mm, but the minimum and maximum temperature distributions are concentrated in the upper left and lower right corners of the section respectively. At z = 240 mm, it can be observed that the lower temperature values are mainly concentrated in the lower left part of the diagonal between points (250, 0) and (0, 250), indicating that when the hot plume rises to the top, it is cooled by the upper wall and then gradually decreases as it moves along the upper wall. After moving diagonally to the upper wall, it becomes a cold plume flowing downwards. Figure 9 is the mean temperature contour of the corresponding horizontal section, and it can be observed that the temperature distribution is more uniform. The high temperature of the z = 10 mm section is concentrated in the

upper right part of the diagonal and near the upper and right walls. The low temperature and high temperature of the z = 125 mm section are concentrated in two diagonal parts, respectively. The high temperature of the z = 240 mm section is concentrated in the lower left part of the diagonal, mainly concentrated in the left and lower walls. In addition, it can be seen that the average temperature of the horizontal section presents an approximately symmetrical structure with the diagonal axis of the crossing points (250, 0) and (0, 250) as the symmetry axis.

Heat flux refers to the amount of heat transferred by conduction, convection, and radiation through a certain area of an object in a unit of time when there is a temperature difference on both sides. Figure 10 shows the average heat flux distribution on the top and bottom surfaces of the cubic cavity. The bottom plates in Figure 10(a) and Figure 10(c) are used as heat plates higher than the initial water temperature, and the temperature gradient between the bottom plate and the water is negative, so the heat flux on the bottom surface is always positive; The dark red area of heat flux is caused by the decrease of the cold plume, increasing the absolute value of temperature gradient. Similarly, Figure 10(b) and Figure 10(d) show that the top surface of the cold plate absorbs heat from the water, while the temperature gradient between the bottom plate and water is positive, resulting in a constant negative value of heat flux. The dark blue area of heat flux is caused by the rise of the hot plume, increasing the absolute value of the temperature gradient.

Through the analysis of the two-dimensional structure, it can be seen that the fluid in the cubic cavity tends to move along the cross-section where the diagonal lines passing through points (0, 0, 0) and (250, 250) are located.

An instantaneous feature of RBC in turbulent flow regimes is the presence of LSC, which is a large-scale identifiable convective structure composed of cold and hot plumes maintained by its dynamics and depends not only on the state of turbulence but also on geometric constraints. Figure 11(a) and Figure 11(b) are instantaneous iso-surface of temperature when the vertical temperature difference is 6.7˚C and 20˚C, respectively. It is observed that a hot plume flows upward in one corner of the cubic cavity, and a cold plume flows downward in the other corner, indicating that the fluid in the cubic cavity is driven by temperature to form an LSC along the diagonal. By analyzing the motion trend of three-dimensional cold and hot plumes, it can be seen that this corresponds to the two-dimensional structure we previously studied.

To further study the large-scale flow structure, the three-dimensional mean velocity field was evaluated. Figure 12 shows a 3D perspective view of a mean streamline color coded with mean velocity magnitude. The LSC is oriented along a black arrow indicating the direction of flow. As can be seen from Figure 12(a) and Figure 12(c), the flow line shows that the hot fluid rises through the cubic cavity and strikes the left corner of the top plate. Near the top plate, the flow occurs horizontally, along the diagonal plane direction where the arrow is located. Similarly, the cold fluid descends along the right side and strikes the bottom plate. In the other diagonal plane, the fluids that rise and fall from the bottom and top plates converge with each other at intermediate heights, and flow occurs from the corners inward into the interior of the cubic cavity. In addition, when the hot fluid reached the roof, a small reverse vortex can be seen near the upper left corner of the cubic cavity. As shown in Figure 12, this interaction between the rising hot fluid and the cold plume generated near the roof results in a recirculation vortex. Similar reverse corner vortices will occur near the lower right side of the box near the bottom plate. To quantify these vortices, the mean vortex structure is shown by the Q criterion, $Q=1/2\left({\Omega }_{ij}{\Omega }_{ij}-S_{ij}{S}_{ij}\right)$ [19] , which is based on a second invariant based on the velocity gradient tensor. As shown in Figure 13, under the vertical temperature difference of 6.7˚C and 20˚C, vortices

shaped as a boomerang appeared at the two corners of the cubic cavity. In this cubic cavity with an aspect ratio of 1, the mean flow structure exhibits an approximately symmetric structure with a diagonal plane as the symmetry plane and is significantly affected by temperature deflection.

5. Conclusion

This paper conducted a series of numerical simulations on the Rayleigh-Bénard convection in the cubic cavity and compared the results with those of references. Based on the verification of the grid independence, we obtained the mean velocity profiles and velocity pulsation profiles of the vertical section in the middle of the cubic cavity, and the numerical results are consistent with the experimental results. It is found that the mean velocity fields of the middle vertical cross-section are very similar in both cases of the vertical temperature difference of 6.7˚C and 20˚C. A counterclockwise large eddy current occupied the entire cavity, and two small clockwise eddies appeared in the lower left and upper right corners, which look like boomerangs. Moreover, we studied the instantaneous temperature fields and mean temperature fields of cross-sections, and three-dimensional isothermal surface structure as well, and found the overall flow is along the diagonal direction, which is also confirmed by time-averaged streamline. And the basic turbulent structure in RB thermal convection includes the rising and falling plumes generated by buoyancy effects.

Acknowledgements

This work is supported by Shenzhen Science and Technology Program (No. JCYJ20220818102012024).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Xie, Y.C., Zhang, L., Ding, G.Y., Chen, X., Xi, H.D. and Xia, K.Q. (2023) Progress in Turbulent Thermal Convection in the Past Decade and Outlook. Advances in Mechanics, 53, 1-47.
[2]	Di Prima, R.C. (1961) Some Variational Principles for Problems in Hydrodynamic and Hydromagnetic Stability. Quarterly of Applied Mathematics, 18, 375-385. https://doi.org/10.1090/qam/116767
[3]	Travis, B.J. anderson, C., Baumgardner, J., Gable, C.W., Hager, B.H., O’Connell, R.J., et al. (1990) A Benchmark Comparison of Numerical Methods for Infinite Prandtl Number Thermal Convection in Two-Dimensional Cartesian Geometry. Geophysical Fluid Dynamics, 55, 137-160. https://doi.org/10.1080/03091929008204111
[4]	Wakashima, S. and Saitoh, T.S. (2004) Benchmark Solutions for Natural Convection in a Cubic Cavity Using the High-Order Time-Space Method. International Journal of Heat and Mass Transfer, 47, 853-864. https://doi.org/10.1016/j.ijheatmasstransfer.2003.08.008
[5]	Michalek, T., Kowalewski, T.A. and Sarler, B. (2005) Natural Convection for Anomalous Density Variation of Water: Numerical Benchmark. Progress in Computational Fluid Dynamics, 5, 158-170. https://doi.org/10.1504/PCFD.2005.006751
[6]	Penot, F., N’Dame, A. and Le Quere, P. (1990) Investigation of the Route to Turbulence in a Differentially Heated Cavity. International Heat Transfer Conference, Jerusalem, 19-24 August1990, 417-422. https://doi.org/10.1615/IHTC9.2410
[7]	Barannikov, V.A., Frik, P.G. and Shaidurov, V.G. (1988) Spectral Characteristics of Two-Dimensional Turbulent Convection in a Vertical Slot. Journal of Applied Mechanics and Technical Physics, 29, 196-200. https://doi.org/10.1007/BF00908581
[8]	Betts, P.L. and Bokhari, I.H. (2000) Experiments on Turbulent Natural Convection in an Enclosed Tall Cavity. International Journal of Heat and Fluid Flow, 21, 675-683. https://doi.org/10.1016/S0142-727X(00)00033-3
[9]	Salat, J., Xin, S., Joubert, P., Sergent, A., Penot, F. and Le Quere, P. (2004) Experimental and Numerical Investigation of Turbulent Natural Convection in a Large Air-Filled Cavity. International Journal of Heat and Fluid Flow, 25, 824-832. https://doi.org/10.1016/j.ijheatfluidflow.2004.04.003
[10]	Valencia, L., Pallares, J., Cuesta, I. and Grau, F.X. (2004) Turbulent Rayleigh-Bénard Convection of Water in Cubical Cavities: A Numerical and Experimental Study. International Journal of Heat and Mass Transfer, 50, 3203-3215. https://doi.org/10.1016/j.ijheatmasstransfer.2007.01.013
[11]	Worner, M. and Grotzbach, G. (1998) Pressure Transport in Direct Numerical Simulations of Turbulent Natural Convection in Horizontal Fluid Layers. International Journal of Heat and Fluid Flow, 19, 150-158. https://doi.org/10.1016/S0142-727X(97)10019-4
[12]	Bairi, A. and de María, J.M.G. (2013) Numerical and Experimental Study of Steady State Free Convection Generated by Constant Heat Flux in Tilted Hemispherical Cavities. International Journal of Heat and Mass Transfer, 66, 355-365. https://doi.org/10.1016/j.ijheatmasstransfer.2013.07.038
[13]	Bondareva, N.S. and Sheremet, M.A. (2017) 3D Natural Convection Melting in a Cubical Cavity with a Heat Source. International Journal of Thermal Sciences, 115, 43-53. https://doi.org/10.1016/j.ijthermalsci.2017.01.021
[14]	Wei, Y., Wang, Z., Yang, J., Dou, H.S. and Qian, Y. (2015) A Simple Lattice Boltzmann Model for Turbulence Rayleigh-Bénard Thermal Convection. Computers & Fluids, 118, 167-171. https://doi.org/10.1016/j.compfluid.2015.06.003
[15]	Peng, S.H., Hanjalic, K. and Davidson, L. (2006) Large-Eddy Simulation and Deduced Scaling Analysis of Rayleigh-Bénard Convection up to Ra = 109. Journal of Turbulence, 7, N66. https://doi.org/10.1080/14685240600953462
[16]	Bosshard, C., Dehbi, A., Deville, M., Leriche, E., Puragliesi, R. and Soldati, A. (2013) Large Eddy Simulation of the Differentially Heated Cubic Cavity Flow by the Spectral Element Method. Computers & Fluids, 86, 210-227. https://doi.org/10.1016/j.compfluid.2013.07.007
[17]	Huang, X.J., Hu, Y.P. and Li, Y.R. (2019) Aspect Ratio Dependence of Rayleigh-Bénard Convection of Cold Water near Its Maximum Density in Box-Shaped Containers. Physics of Fluids, 31, Article ID: 075107. https://doi.org/10.1063/1.5097964
[18]	Vasiliev, A., Sukhanovskii, A., Frick, P., Fomichev, V., Bolshukhin, M. and Romanov, R. (2016) High Rayleigh Number Convection in a Cubic Cell with Adiabatic Sidewalls. International Journal of Heat and Mass Transfer, 102, 201-212. https://doi.org/10.1016/j.ijheatmasstransfer.2016.06.015
[19]	Dubief, Y. and Delcayre, F. (2000) On Coherent-Vortex Identification in Turbulence. Journal of Turbulence, 1, N11. https://doi.org/10.1088/1468-5248/1/1/011