# Simulation of Turbulent Convection at High Rayleigh Numbers

Authors: Sergey Dmitriev; Andrey Kozelkov; Andrey Kurkin; Nataliya Tarasova; Valentin Efremov; Vadim Kurulin; Roman Shamin; Maksim Legchanov
Journal: Modelling and Simulation in Engineering
 (2018)

---

## Abstract

The paper considers the possibility of using different approaches to modeling turbulence under conditions of highly developed convection at high Rayleigh numbers. A number of industrially oriented problems with experimental data have been chosen for the study. It is shown that, at Rayleigh numbers from 109 to 1017, the application of the eddy-resolving LES model makes it possible to substantially increase the accuracy of modeling natural convection in comparison with the linear vortex viscosity model SST. This advantage is most pronounced for cases of a vertical temperature difference with the formation of a large zone of convection of strong intensity. The use of the Reynolds stress model EARSM is shown for cases of natural convective flow in domains with dihedral angles in the simulated region and the predominance of secondary currents. When simulating a less intense convective flow, when the temperature difference is reached at one boundary, the differences in the approaches used to model turbulence are less significant. It is shown that, with increasing values of Rayleigh numbers, errors in the determination of thermohydraulic characteristics increase and, for more accurate determination of them, it is expedient to use eddy-resolving approaches to the modeling of turbulence.


---

## Body

## 1. Introduction

Flows of liquids with developed turbulent convection, characterized by a strong influence of gravitational forces, are of particular interest because this phenomenon still is poorly understood. The need to study such currents arises in many branches of science and technology, such as astrophysics, geophysics, geodynamics, atomic energy, and others. The most important tasks, where turbulent convection is necessary, include large-scale currents in the atmosphere and liquid cores of planets, the flow of liquid-metal coolants in reactor installations, convection in microgravity in spacecraft and fuel tanks, and other diverse flows in technical products. Particular interest in the study of turbulent convection is manifested in the nuclear energy industry. Along with the increase of technical and economic parameters of nuclear power plants (NPPs), the security of reactor facilities during severe accidents and minimizing their consequences have been of particular importance [1]. In severe accidents accompanied by loss of the coolant, the reactor core is destroyed. As a result, the molten elements of the reactor and reactor core structure are moved to the bottom of the reactor vessel, which results in the formation of a high-temperature fuel melt in its lower part, where it must be localized. The measure to ensure the safety of the reactor facility is to ensure the necessary heat removal from the boundaries of the melt localization device, where heat exchange takes place in the natural convection mode.Experiments in nuclear power engineering, where the mechanism of convective heat transfer of a liquid in cavities with volume heating was studied at Rayleigh numbers of high orders of 1014 −1017, include the COPO [2–5] and BALI [6, 7]. The complexity of experimental studies of this class of currents actualizes the development of numerical simulation. Since in the course of experimental studies the pronounced turbulent nature of the flow is confirmed, one of the main problems in the numerical modeling of these flows is turbulence. The question of the influence of turbulence on the nature of the flow and the magnitude of the boundary heat transfer is of interest. The existing level of computing power allows modeling convection, at an acceptable time, using the Reynolds-averaged Navier-Stokes equations (RANS approach) [8, 9]. This imposes limitations on the description of flow turbulence, for both natural and forced convection, due to the presence of gravitational forces modeled in the RANS approach by empirical relationships, which leads to a significant error in describing the process as a whole. This problem can be solved by using direct numerical simulation (DNS). However, the currently conducted numerical experiments using DNS for natural convection problems are limited to modest Rayleigh numbers [10] (maximum 107). Interest in the study of turbulent convection is the high Rayleigh numbers. For technical applications, these are of orders 1014–1017. For compulsory convection problems, such as mixing high-speed multitemperature flows and jets characterized by strong buoyancy forces, numerical calculations using DNS are limited to small values of Reynolds numbers (of the order of 105), which leads to constraints on simulated flow regimes. Of interest in scientific and technical applications are problems with Reynolds numbers of the order of 107 and higher. In most cases, to simulate the problems with such values of the Reynolds number an approach based on the solution of the Navier-Stokes equations averaged over the Reynolds numbers (RANS) is used, supplemented by one of the turbulence models to close the system [11]; however, in [12], it is noted that the standard model of RANS models is not always applicable for solving problems of turbulent convection. Numerical modeling of the BALI experiment for the purpose of verification and calibration of various models of turbulence was carried out in [13]. It was shown in [14] that the use AKN turbulence model (the AKN model is described in [15–17]) improves the prediction of the heat flux distribution. Modeling using the standard k-ε model of turbulence shows that its application leads to an underestimated heat flux distribution. At the same time, an increase in the accuracy and informativeness of the results for the problems of natural convection can be achieved by using eddy-resolving turbulence models, for example, LES, the successful application of which can be found in [8, 9].The purpose of this paper is to numerically simulate the processes of highly developed convection of a fuel liquid at high Rayleigh numbers and to study the application of approaches to the description of turbulence for freely convective currents. The results of modeling based on RANS, RSM, and LES approaches are presented.The first part of the paper presents the results of modeling a freely convective flow between two walls of different temperatures in a closed cubic cavity, which represents the most common class of natural convection problems. In the main part of the paper, we present the results of a numerical study of the problems of natural convection at high Rayleigh numbers using the example of simulation of the COPO and BALI experiments. Results of application of various approaches to modeling of turbulence in comparison with results of experimental researches are shown.


## 2. Approaches to Modeling Turbulence

The modeling of turbulence is based on several basic approaches. RANS models of turbulence are currently the most popular; however, they are not universal and suitable for a wide range of applications, since they only describe the averaged characteristics, which imposes certain requirements on their applicability in practice [8, 18]. The highest rating of applicability from RANS models is model SST (Shear Stress Transport) [7]. However, as in all RANS models, SST uses the Boussinesq hypothesis on turbulent viscosity [19], which is valid only in the case of isotropic turbulence. For turbulent convection, the anisotropic properties of the flow have a significant effect; this is due to the complex nature of the flow and the presence of secondary currents. To correctly predict the structure of such flows, it is necessary to use RSM models (Reynolds Stress Modeling) or alternative EARSMs (Explicit Algebraic Reynolds Stress Modeling) that take into account the influence on the main flow of all components of the Reynolds stress tensor with the help of nonlinear relationships and thus can increase the accuracy of turbulent calculations, anisotropic flows [19]. In this paper, as an alternative to the SST model, an explicit algebraic EARSM model is considered, which by its qualities is not inferior to differential RSM and allows significantly reducing computational costs [20].The method of large eddy LES (Large Eddy Simulation) makes it possible to obtain good results for both attached and detached flows and to substantially refine the prediction of a number of basic physical processes [8, 9, 21]. This approach can provide detailed information about nonstationary fields of fluctuations of velocity, temperature, and pressure, which in turn can influence the magnitude and nature of the boundary heat transfer. However, the LES method imposes great demands on the quality of discrete models and significantly increases, in comparison with RANS, the amount of required computational resources. Thus, for convective currents, there arises the question of studying the structure of the turbulent flow and the application of eddy-resolving turbulence models for establishing the connection between small-scale oscillations and heat transfer in the zone under consideration.


## 3. Numerical Experiments

All calculations in the present work were carried out with the help of the LOGOS package. This package is the Russian software product of engineering analysis intended for solving conjugate three-dimensional problems of convective heat and mass transfer, aerodynamics, and hydrodynamics on parallel computers. The LOGOS software package has successfully passed the verification and has shown quite good results on a series of various hydrodynamic problems [8, 9, 17, 22, 23], including calculations of turbulent and nonstationary flows [20, 21] and geophysical phenomena [24–26].


### 3.1. The Problem of Developed Turbulent Convection

Simulation of free-convective fluid flow in a cubic cavity between two different-temperature walls is of fundamental interest, since currents with a superimposed vertical temperature difference are characterized by the presence of different-scale vortex structures considering the main features of naturally convective flow and are widely used in engineering applications. The paper presents the results of numerical simulation of the circulation of distilled water in a domain heated from below and cooled from above. The basic research of this issue is carried out for convection in cylindrical and cubic cavities. Experimental studies of convection in a cubic cavity with a vertical temperature difference were performed in a wide range of Rayleigh numbers 103–109.The experimental setup is a cubic region with a sideD = 250 mm [10]. The horizontal walls are made of copper and they act as heat exchangers, and the vertical walls are made of plexiglass and provide isothermal boundary conditions. With the help of thermostats, a thermostatic fluid is passed through the heat exchangers, which heats (cools) the heat exchangers with respect to the average temperature of the liquid. The experiment is described in detail in [10]. In this paper we consider a flow characterized by a difference of temperature of 20°C, which corresponds to the Rayleigh number 6,1 × 109, which is determined by the following formula:(1)Ra=gβΔTH3λν,where g is the acceleration due to gravity, β is the thermal expansion coefficient, H is the characteristic dimension, λ is the thermal conductivity coefficient, ν is the kinematic viscosity, and ΔT is the temperature difference. Water is used as a coolant, the properties of which are expressed in the form of dependencies [27]:(2)μ=0.2574-0.00331T+1.77e-5T2-5.03e-8T3+7.997e-11T4-6.7373e-14T5+2.349e-17T6,Cp=-81756.058+1075.58T-5.325T2+0.013T3-1.57972e-5T4+7.6e-9T5,λ=-0.4937784389+0.005793217949T-7.074592075e-6T2,ρ=754.5044626+1.872414452T-0.003557109557T2.For calculations, an isotropic block-structured grid consisting of 3.375 million cells of regular hexahedral shape with a size of 0.0016 m was used. This grid model corresponds to the requirement of applicability of standard wall functions with parametery+ no higher than 5 and also allows small-scale turbulence to be resolved within the LES calculation.A detailed and isotropic grid allows for RANS calculations to use a countercurrent scheme of the first order. For the LES calculation, a fully implicit scheme was used, increasing the accuracy of the calculation and allowing in a stable mode using a time step corresponding to the Courant number equal to unity. A description of the scheme is given in [9, 21].The nonstationary calculation was carried out up to the time point of 4200 seconds, which corresponds to the estimated time of the main characteristics of the physical fields in the experiment. The averaging of the basic physical fields was carried out from the instant of time of 100 seconds. Figure1 shows calculated and experimental fields of instantaneous and averaged velocity in the central section of the experimental setup.Figure 1
The velocity field distribution: instantaneous velocity field, mm/s (a); averaged velocity field, mm/s (b) (experimental data are presented on the left panel, and calculations on the right panel).
(a)
(b)The maximum value of the velocity in the picture of the experimental instantaneous field is reached near the walls; vortices appear at opposite angles of the region; the numerical field of instantaneous velocity represents a similar picture, describing small-sized vortices localized in the corners. In the experiment and in the calculation, the structure of the average velocity fields is very close. It is clearly seen that a large-scale vortex appears in the region commensurate with its size. In addition to the main vortex in the middle fields, there are also secondary vortices of smaller scale localized in opposite corners of the cubic region.For qualitative and quantitative evaluation, Figure2 shows the profiles of the horizontal and vertical velocity components averaged over time. As can be seen from the figures, the velocity profiles show good agreement, not only repeating the flow structure but also reproducing the velocity value in the boundary layers.Figure 2
Profiles of the time-average profile of the horizontal (a) and vertical (b) velocity components for a fixed coordinatex = 125 mm and z = 125 mm, respectively.
(a)
(b)The velocity profiles obtained with the SST linear vortex viscosity model deviate from the experiment. The root-mean-square deviation for this model was about 6%. Differences in the calculated data can be related to the geometric feature of the region characterized by a large number of dihedral angles. Thus, in [23], the problem of the flow of a liquid in a tube with a square cross section shows that the SST model gives a large error and for such problems the application of the EARSM model allowing modeling of the influence on the main flow of all components of the Reynolds stress tensor is shown. Using this model significantly reduced the computational error, the standard deviation of the EARSM results was about 2.5%.Using the LES model, the best agreement was obtained with the experimental data. The calculated profiles of the averaged velocity according to LES reproduce the general character of the convective heat transfer in a cubic cavity heated from below and cooled from above with a high accuracy; the error of the results was ~1.5%. Note here that the results for all models are obtained on the same calculation grid.The numerical solution of this problem shows that the simulation of convective motion for the numbers Ra ~109 over all models of turbulence yields an acceptable result. However, practically important currents for the safety of nuclear power plants are characterized by higher Rayleigh numbers, up to 1017. The applicability of the presented approaches to turbulence modeling for such Rayleigh numbers can be estimated in the experiments of COPO and BALI, which are presented below.


### 3.2. Modeling the COPO Experiment

In the COPO experiment, the problem of developed turbulent convection of an aqueous solution of salt in a closed region is considered. The COPO experiment simulates an emergency in a water-water power reactor (WWPR) at which the core is destroyed, which leads to the formation of a high-temperature fuel melt in the lower part of the reactor vessel [1–5]. To prevent further destruction of the reactor facility, heat is removed from the reactor vessel by cooling its lower part with boiling water, the upper part of the reactor is cooled by radiation heat exchange, and the flow within the region occurs in the natural convection mode.The experiment was carried out at a scale of 1 : 2. In slice geometry, the installation length was 1.77 m, the thickness of the layer was 0.094, the large horizontal axis of the half-ellipse wasa = 0.89 m, and the vertical axis was b=1/2a. The general view of the experimental setup is shown in Figure 3; it is a thin vertical slice corresponding to the center of the lower part of the WWPR-440 reactor shell, with an elliptical bottom.Figure 3
COPO setup scheme.The liquid simulating the melt was a solution of salt in water. The volumetric heating of which provided distributed electrodes throughout the all domain. Vertical boundaries were thermally insulated. The lateral and upper boundaries were cooled by liquid nitrogen through the interlayer of the heat insulator, which led to the formation of a relatively thin crust of ice and provided isothermal boundary conditions. In one of the modes, the upper bound of the experimental setup was thermally insulated, which is analogous to the case in which the upper oxidized basin is covered in the melt localization device and the heat release is not large. For the remaining series of experiments, the upper bound was cooled similarly to a cylindrical one, which provided a constant temperature. This mode corresponds to the strongest heat sink from the upper boundary of the melt.The main characteristics of the flow with internal heat sources are the parameters of the heat flux density at the boundary and the maximum excess of the liquid temperature in the volume above the boundary temperature. The dimensionless numbers characterizing these quantities are the Nusselt and Rayleigh criteria, respectively. For a heat-generating melt, a modified Rayleigh number is used:(3)Rai=gβqνR5λaTν,where g is the acceleration of gravity, β is the coefficient of thermal expansion, qV is the volume volumetric energy release, R is the height of the half-ellipse, λ is the coefficient of thermal conductivity of the medium, aT is the thermal diffusivity of the medium, and ν is the kinematic viscosity. The value of the modified Rayleigh number in this experiment was ~1015.The indicator of effective heat transfer through the volume boundary is the Nusselt number:(4)Nu=qHλΔT.The all thermophilically quantities are given by temperature-dependent (2), and the coefficient of temperature expansion, depending on the temperature, is given by the following:(5)β=0.007114T-1.87241445754.5045+1.8724T-0.00355711T2.In numerical experiments, we consider a mode with a heat-insulated top wall (mode 1) and also three regimes with different energy release and different height of the region (aspect ratio). Table1 presents the main characteristics of the modes.Table 1
Parameters of experimental regimes.
Mode
Height of the area, m
Volumetric energy release, kW/m2
1
0.847
46.5
2
0.639
153.2
3
0.73
161.5
4
0.847
104.3The calculated region corresponded to the complete model of the experimental setup. A block-structured grid consisting of 2.5 million cells of a hexahedral shape was used and the cell size in the main region was 0.004 m thick, with concentration to all walls (the size of the wall cell was 0.0002 m). The grid model corresponds to the requirement of applicability of standard wall functions with parametery+ no higher than 5. This allows considering the small-scale turbulence in the LES models.A detailed and isotropic grid allows for RANS calculations to use a countercurrent scheme of the first order. The LES as scheme discretization in space and the equation of conservation of momentum was used with the BCD scheme calibrated for this scheme constant Smagorinsky. The BCD scheme is a scheme of limited central differences, where a diagram of normalized variable is used to formulate the restrictive condition, similarly to the GAMMA scheme [28]. The description of the scheme BCD is given in [9, 21].The time for setting the heat-hydraulic process in the calculations is about 5000 seconds. The averaging of the basic physical fields was carried out from the time point of 2000 seconds. A series of nonstationary calculations was performed; the time step for RANS models was 1 s, for LES – 0.1 s. It has been experimentally confirmed that, for regime number 1 with a heat-insulated top wall, a stratified temperature distribution is typical [5]. Figure 4 shows the calculated temperature field, it can be seen from the figure that the liquid with a higher temperature (minimum density) is at the top of the region, there are practically no disturbances in the flow, and a clear temperature stratification is observed.Figure 4
Temperature distribution.The results of calculating the experimental regime 4 with a cooled upper boundary are below. For this regime, the stratified temperature distribution is observed only in the lower part of the calculated region, and most of the upper region is occupied by the Rayleigh-Bener convection zone, where large vortex structures are observed. Figures5 and 6 show the velocity and temperature distributions using various approaches to modeling turbulence. As can be seen from the figures, the temperature fields for SST and EARSM have the same structure, and a similar distribution is observed in the fields of the averaged temperature for LES. On the fields of instantaneous velocity for all models of turbulence in the upper part of the region, two large vortices are observed and the maximum velocity is observed near the cooled cylindrical boundary.Figure 5
The temperature distribution in the central section of the region: instantaneous temperature (a) and average temperature (b).
(a)
(b)Figure 6
The velocity distribution in the central section of the region: instantaneous velocity, m/s (a) and average velocity, m/s (b).
(a)
(b)For the experimental mode number 1 in Figure7(a), the distribution of the heat flux on the cooled elliptic boundary is compared with the experimental data. The experiment contains data on the distribution of heat flux density, both on the left side of the body and on the right side. The difference in these graphs seems to be due to the purity of the experiment. In the case of modeling, the difference in the flow between the right and left parts can be due only to an insufficient averaging of the results over time. The SST and EARSM models represent a similar result for the distribution of heat flux density, the root-mean-square error is about 6%, and the use of EARSM gives a slight advantage, reducing the error by 0.5%. The LES model refines the prediction; the mean square deviation error of the results is about 3.5%.Figure 7
Heat flow distribution: mode with a heat-insulated upper boundary (a) and mode with a cooled one (b).
(a)
(b)For the experimental regime with a cooled upper boundary (Figure7(b)), the results have a higher error. All models overestimate the heat flow in the upper part of the calculation area. The root-mean-square deviation of the EARSM and SST results is 11% and 13%, respectively, but the LES as a whole is closer to the experiment and the error was about 6%. The results using the SST and EARSM models relate giving approximately the same result, but the EARSM model still slightly clarifies the distribution of the heat flux on the boundary.The main integral indicator for tasks with internal energy release and heat dissipation from boundaries is the distribution of the boundary heat transfer, characterized by the Nusselt number. Based on the results of recent experimental studies, the following dependence of the Nu number on the modified Ra number was obtained (the dependence was established for Ra ~1013–1017 [3]):(6)Nu=0.385Rai0.233.In Figure8, in accordance with Table 1, the results of numerical simulation of three experimental regimes are presented. For all variants, the dependence of the Nusselt number on the modified Rayleigh number for various approaches to the modeling of turbulent flows is presented.Figure 8
Dependence of the Nusselt number from the modified Rayleigh number.Analyzing the results, we can assume that a sufficiently large deviation is associated with a understating of the heat flux at the upper boundary, while at the cylindrical boundary the heat flux is overestimated. The heat dissipation from the upper boundary is less intensive: in comparison with the experimental data, on the elliptical boundary it is overestimated by an average of 30%. The maximum deviation from the experimental dependence reaches 33%; for the LES, the maximum deviation was 25%. Thus, the use of the LES approach refines the forecasting, but the error remains substantial, which may be due to the overestimation of the heat flow through the upper boundary.In the framework of the COPO experiment, a convective flow (with Rayleigh number ~1015 [2–5]) is considered. Since the currents accompanying accidents in reactor installations are characterized by Rayleigh numbers of even higher orders, considerable interest is in the study of convection at Rayleigh numbers ~1017. The experiment that covers the region most important for the safety problems of nuclear reactors is the BALI experiment.


### 3.3. Modeling the BALI Experiment

Experiments of BALI were carried out in order to investigate the convection of the fuel at a serious accident of the PWR reactor in the second half of the 1990s. In contrast to the COPO experiment, BALI experiments simulate the processes occurring during a severe accident in a reactor with a spherical bottom. Four series of experiments were performed under different conditions and configurations on a 1 : 1 scale. The general view of the experimental setup is shown in Figure9 and represents a cutout along the axial plane [6]. The liquid simulating the melt was a solution of salt in water. The volumetric heating provided by electrodes was distributed throughout the region. The side walls are thermally insulated; the cylindrical boundary was cooled by liquid nitrogen through the interlayer of the heat insulator, which led to the formation of a relatively thin crust of ice and provided isothermal boundary conditions. In the first series of experiments, the upper boundary was thermally insulated; in the second series of experiments, it was cooled similarly to the cylindrical boundary. The Nusselt number (4) is an indicator of the effective heat transfer across the volume boundary.Figure 9
BALI installation scheme.The calculated region corresponded to the complete model of the experimental setup. To simulate the BALI experiment, a block-structured grid was constructed, consisting of 4.2 million cells of a hexahedral shape, the cell size in the main region 0.006 m thick with concentration to all walls, and the size of the wall cell 0.0001 m. As a space-sampling scheme and in the pulse-conservation equation for the LES model, we used the BCD scheme with the Smagorinsky constant. The time step was chosen in such a way as to ensure the value of the Courant number of the order of unity. The thermophilically properties of the coolant used in modeling the experiments are given by temperature-dependent (2), and the coefficient of variable temperature expansion, depending on the temperature, is determined by relation (5).It should be noted that in [6] the results of the experiment are given in dimensionless form in relative units and there are no exact data on volumetric energy release, cooling conditions, and about ice crust on cooled surfaces. It was shown in [6, 7] that the heat flux profile obtained in various experiments is universal and does not depend on the Rayleigh number; however, the shape of the basin and the aspect ratio greatly influence the heat transfer distribution [1, 6].The output parameters for the heat flux and temperature are given in a dimensionless form. The physical picture for the first calculation case, with the thermally insulated top wall, is similar to the COPO problem, the temperature field is clearly stratified, the liquid with the higher temperature (minimum density) is distributed at the top of the region, there is practically no disturbance in the flow, and the greatest value of the velocity arises near the cooled cylindrical wall.In calculations of the current with the adiabatic upper boundary, volumetric energy release was varied, it was selected in such a way as to provide the following values of the Rayleigh criterion:Rai = 2 × 1016 (Ra1), Rai = 3 × 1016 (Ra2). Figure 10(a) shows the distribution of the dimensionless temperature along the vertical boundary and the distribution of the dimensionless heat flux along the cylindrical boundary. It can be seen from the graphs that the profiles of temperature and density of the heat flux are universal and almost independent of the Rayleigh number. The results of the first series of numerical experiments show that the EARSM and SST turbulence models accurately reproduced the general character of convective heat transfer in the experimental BALI setup. The distribution of the dimensionless temperature on the vertical wall obtained during the calculation agrees with the experimental data for all turbulence models and for various Rayleigh numbers. The distribution of the dimensionless heat flux density on the cooled cylindrical boundary is in good agreement with the experimental data in the lower part of the experimental setup, and its magnitude is underestimated in its upper part.Figure 10
The distribution of the time-averaged dimensionless temperature along the vertical boundary (a) and the dimensionless density of the heat flux at the cylindrical boundary (b).
(a)
(b)The SST and EARSM models represented a similar result of the distribution of heat flow and temperature; the root-mean-square error for the heat flux is about 6%. The use of EARSM gives a slight advantage, reducing the error by 1%. The LES model refines the prediction; the rms error of the results is about 4%.The results of the numerical simulation of the problem of BALI for the regime with a cooled upper boundary are below. Several flow zones are observed: the upper unstable layer of a homogeneous temperature, where the Rayleigh-Boehner convection zone occupies the most part; a lower zone with a temperature stratification at which the liquid rises to the center at a low velocity; boundary layer where the liquid in contact with ice cools and returned at a high speed to the bottom of the basin. A similar picture is observed in numerical simulation; a large part of the domain is occupied by a zone of large-scale convection where large vortex structures are localized; heat exchange occurs in the regime of highly developed turbulent convection. For this mode, Figure11 shows the instantaneous velocity field and the instantaneous distribution of the dimensionless temperature, with a Rayleigh criterion of 3.5 × 1016 for the LES turbulence simulation.Figure 11
Instantaneous temperature field (a) and instantaneous velocity field (b).
(a)
(b)The stratified temperature distribution is observed only in the lower part of the calculation area. The greater part of the upper domain is occupied by the convection zone of Rayleigh-Bener and it is much larger than the zone observed in the COPO simulation. This may be due to both the geometry of the installation and the increased value of the number Rayleigh.The experiment contains data on the distribution of temperature and heat flux for the flow regime with cooling walls. The results are presented for the averaged density of the heat flux (Figure12(b)) and the distribution of the dimensionless temperature (Figure 12(a)).Figure 12
The distribution of the time-averaged temperature along the vertical boundary (a) and the distribution of the heat flux density at the cylindrical boundary (b).
(a)
(b)Analysis of the calculation results shows that the RANS models underestimate the heat flow in the upper part of the calculation area. The LES results are closer to the experiment. The root-mean-square deviation of the EARSM and SST results is 11% and 13%, respectively. The LES results are 6%. The results using the SST and EARSM models correlated, giving approximately the same result, but the EARSM model lies closer to the experimental data.Figure13 shows the results of numerical simulation of three experimental regimes with a cooled upper wall for the numbers Ra1 = 6.82 × 1015, Ra2 = 3.16 × 1016, and Ra3 = 1.07 × 1017 in comparison with the obtained Nusselt number dependence on the modified Rayleigh number (6) for EARSM and LES turbulence models.Figure 13
Dependencies of the number Nu on the modified number Ra.The maximum error for the Nusselt number at the upper boundary was about 36%, the vortex-resolving LES Smagorinsky model refines the prediction, and the maximum error using this model is about 33% in the case of the maximum Rayleigh number. Numerical modeling qualitatively corresponds to the physical picture of natural convection. Due to the high-temperature difference, large vortex flows occur at the upper wall. In the lower part of the basin, there is a stable zone of temperature stratification. This picture corresponds to the flow in the experimental setup. The distribution of the dimensionless temperature at the vertical boundary is described with sufficient accuracy, with the greatest deviations being observed in the central area. The calculated density distribution of the dimensionless heat flux at the cooled upper boundary is described with some deviation from the experimental data.


## 3.1. The Problem of Developed Turbulent Convection

Simulation of free-convective fluid flow in a cubic cavity between two different-temperature walls is of fundamental interest, since currents with a superimposed vertical temperature difference are characterized by the presence of different-scale vortex structures considering the main features of naturally convective flow and are widely used in engineering applications. The paper presents the results of numerical simulation of the circulation of distilled water in a domain heated from below and cooled from above. The basic research of this issue is carried out for convection in cylindrical and cubic cavities. Experimental studies of convection in a cubic cavity with a vertical temperature difference were performed in a wide range of Rayleigh numbers 103–109.The experimental setup is a cubic region with a sideD = 250 mm [10]. The horizontal walls are made of copper and they act as heat exchangers, and the vertical walls are made of plexiglass and provide isothermal boundary conditions. With the help of thermostats, a thermostatic fluid is passed through the heat exchangers, which heats (cools) the heat exchangers with respect to the average temperature of the liquid. The experiment is described in detail in [10]. In this paper we consider a flow characterized by a difference of temperature of 20°C, which corresponds to the Rayleigh number 6,1 × 109, which is determined by the following formula:(1)Ra=gβΔTH3λν,where g is the acceleration due to gravity, β is the thermal expansion coefficient, H is the characteristic dimension, λ is the thermal conductivity coefficient, ν is the kinematic viscosity, and ΔT is the temperature difference. Water is used as a coolant, the properties of which are expressed in the form of dependencies [27]:(2)μ=0.2574-0.00331T+1.77e-5T2-5.03e-8T3+7.997e-11T4-6.7373e-14T5+2.349e-17T6,Cp=-81756.058+1075.58T-5.325T2+0.013T3-1.57972e-5T4+7.6e-9T5,λ=-0.4937784389+0.005793217949T-7.074592075e-6T2,ρ=754.5044626+1.872414452T-0.003557109557T2.For calculations, an isotropic block-structured grid consisting of 3.375 million cells of regular hexahedral shape with a size of 0.0016 m was used. This grid model corresponds to the requirement of applicability of standard wall functions with parametery+ no higher than 5 and also allows small-scale turbulence to be resolved within the LES calculation.A detailed and isotropic grid allows for RANS calculations to use a countercurrent scheme of the first order. For the LES calculation, a fully implicit scheme was used, increasing the accuracy of the calculation and allowing in a stable mode using a time step corresponding to the Courant number equal to unity. A description of the scheme is given in [9, 21].The nonstationary calculation was carried out up to the time point of 4200 seconds, which corresponds to the estimated time of the main characteristics of the physical fields in the experiment. The averaging of the basic physical fields was carried out from the instant of time of 100 seconds. Figure1 shows calculated and experimental fields of instantaneous and averaged velocity in the central section of the experimental setup.Figure 1
The velocity field distribution: instantaneous velocity field, mm/s (a); averaged velocity field, mm/s (b) (experimental data are presented on the left panel, and calculations on the right panel).
(a)
(b)The maximum value of the velocity in the picture of the experimental instantaneous field is reached near the walls; vortices appear at opposite angles of the region; the numerical field of instantaneous velocity represents a similar picture, describing small-sized vortices localized in the corners. In the experiment and in the calculation, the structure of the average velocity fields is very close. It is clearly seen that a large-scale vortex appears in the region commensurate with its size. In addition to the main vortex in the middle fields, there are also secondary vortices of smaller scale localized in opposite corners of the cubic region.For qualitative and quantitative evaluation, Figure2 shows the profiles of the horizontal and vertical velocity components averaged over time. As can be seen from the figures, the velocity profiles show good agreement, not only repeating the flow structure but also reproducing the velocity value in the boundary layers.Figure 2
Profiles of the time-average profile of the horizontal (a) and vertical (b) velocity components for a fixed coordinatex = 125 mm and z = 125 mm, respectively.
(a)
(b)The velocity profiles obtained with the SST linear vortex viscosity model deviate from the experiment. The root-mean-square deviation for this model was about 6%. Differences in the calculated data can be related to the geometric feature of the region characterized by a large number of dihedral angles. Thus, in [23], the problem of the flow of a liquid in a tube with a square cross section shows that the SST model gives a large error and for such problems the application of the EARSM model allowing modeling of the influence on the main flow of all components of the Reynolds stress tensor is shown. Using this model significantly reduced the computational error, the standard deviation of the EARSM results was about 2.5%.Using the LES model, the best agreement was obtained with the experimental data. The calculated profiles of the averaged velocity according to LES reproduce the general character of the convective heat transfer in a cubic cavity heated from below and cooled from above with a high accuracy; the error of the results was ~1.5%. Note here that the results for all models are obtained on the same calculation grid.The numerical solution of this problem shows that the simulation of convective motion for the numbers Ra ~109 over all models of turbulence yields an acceptable result. However, practically important currents for the safety of nuclear power plants are characterized by higher Rayleigh numbers, up to 1017. The applicability of the presented approaches to turbulence modeling for such Rayleigh numbers can be estimated in the experiments of COPO and BALI, which are presented below.


## 3.2. Modeling the COPO Experiment

In the COPO experiment, the problem of developed turbulent convection of an aqueous solution of salt in a closed region is considered. The COPO experiment simulates an emergency in a water-water power reactor (WWPR) at which the core is destroyed, which leads to the formation of a high-temperature fuel melt in the lower part of the reactor vessel [1–5]. To prevent further destruction of the reactor facility, heat is removed from the reactor vessel by cooling its lower part with boiling water, the upper part of the reactor is cooled by radiation heat exchange, and the flow within the region occurs in the natural convection mode.The experiment was carried out at a scale of 1 : 2. In slice geometry, the installation length was 1.77 m, the thickness of the layer was 0.094, the large horizontal axis of the half-ellipse wasa = 0.89 m, and the vertical axis was b=1/2a. The general view of the experimental setup is shown in Figure 3; it is a thin vertical slice corresponding to the center of the lower part of the WWPR-440 reactor shell, with an elliptical bottom.Figure 3
COPO setup scheme.The liquid simulating the melt was a solution of salt in water. The volumetric heating of which provided distributed electrodes throughout the all domain. Vertical boundaries were thermally insulated. The lateral and upper boundaries were cooled by liquid nitrogen through the interlayer of the heat insulator, which led to the formation of a relatively thin crust of ice and provided isothermal boundary conditions. In one of the modes, the upper bound of the experimental setup was thermally insulated, which is analogous to the case in which the upper oxidized basin is covered in the melt localization device and the heat release is not large. For the remaining series of experiments, the upper bound was cooled similarly to a cylindrical one, which provided a constant temperature. This mode corresponds to the strongest heat sink from the upper boundary of the melt.The main characteristics of the flow with internal heat sources are the parameters of the heat flux density at the boundary and the maximum excess of the liquid temperature in the volume above the boundary temperature. The dimensionless numbers characterizing these quantities are the Nusselt and Rayleigh criteria, respectively. For a heat-generating melt, a modified Rayleigh number is used:(3)Rai=gβqνR5λaTν,where g is the acceleration of gravity, β is the coefficient of thermal expansion, qV is the volume volumetric energy release, R is the height of the half-ellipse, λ is the coefficient of thermal conductivity of the medium, aT is the thermal diffusivity of the medium, and ν is the kinematic viscosity. The value of the modified Rayleigh number in this experiment was ~1015.The indicator of effective heat transfer through the volume boundary is the Nusselt number:(4)Nu=qHλΔT.The all thermophilically quantities are given by temperature-dependent (2), and the coefficient of temperature expansion, depending on the temperature, is given by the following:(5)β=0.007114T-1.87241445754.5045+1.8724T-0.00355711T2.In numerical experiments, we consider a mode with a heat-insulated top wall (mode 1) and also three regimes with different energy release and different height of the region (aspect ratio). Table1 presents the main characteristics of the modes.Table 1
Parameters of experimental regimes.
Mode
Height of the area, m
Volumetric energy release, kW/m2
1
0.847
46.5
2
0.639
153.2
3
0.73
161.5
4
0.847
104.3The calculated region corresponded to the complete model of the experimental setup. A block-structured grid consisting of 2.5 million cells of a hexahedral shape was used and the cell size in the main region was 0.004 m thick, with concentration to all walls (the size of the wall cell was 0.0002 m). The grid model corresponds to the requirement of applicability of standard wall functions with parametery+ no higher than 5. This allows considering the small-scale turbulence in the LES models.A detailed and isotropic grid allows for RANS calculations to use a countercurrent scheme of the first order. The LES as scheme discretization in space and the equation of conservation of momentum was used with the BCD scheme calibrated for this scheme constant Smagorinsky. The BCD scheme is a scheme of limited central differences, where a diagram of normalized variable is used to formulate the restrictive condition, similarly to the GAMMA scheme [28]. The description of the scheme BCD is given in [9, 21].The time for setting the heat-hydraulic process in the calculations is about 5000 seconds. The averaging of the basic physical fields was carried out from the time point of 2000 seconds. A series of nonstationary calculations was performed; the time step for RANS models was 1 s, for LES – 0.1 s. It has been experimentally confirmed that, for regime number 1 with a heat-insulated top wall, a stratified temperature distribution is typical [5]. Figure 4 shows the calculated temperature field, it can be seen from the figure that the liquid with a higher temperature (minimum density) is at the top of the region, there are practically no disturbances in the flow, and a clear temperature stratification is observed.Figure 4
Temperature distribution.The results of calculating the experimental regime 4 with a cooled upper boundary are below. For this regime, the stratified temperature distribution is observed only in the lower part of the calculated region, and most of the upper region is occupied by the Rayleigh-Bener convection zone, where large vortex structures are observed. Figures5 and 6 show the velocity and temperature distributions using various approaches to modeling turbulence. As can be seen from the figures, the temperature fields for SST and EARSM have the same structure, and a similar distribution is observed in the fields of the averaged temperature for LES. On the fields of instantaneous velocity for all models of turbulence in the upper part of the region, two large vortices are observed and the maximum velocity is observed near the cooled cylindrical boundary.Figure 5
The temperature distribution in the central section of the region: instantaneous temperature (a) and average temperature (b).
(a)
(b)Figure 6
The velocity distribution in the central section of the region: instantaneous velocity, m/s (a) and average velocity, m/s (b).
(a)
(b)For the experimental mode number 1 in Figure7(a), the distribution of the heat flux on the cooled elliptic boundary is compared with the experimental data. The experiment contains data on the distribution of heat flux density, both on the left side of the body and on the right side. The difference in these graphs seems to be due to the purity of the experiment. In the case of modeling, the difference in the flow between the right and left parts can be due only to an insufficient averaging of the results over time. The SST and EARSM models represent a similar result for the distribution of heat flux density, the root-mean-square error is about 6%, and the use of EARSM gives a slight advantage, reducing the error by 0.5%. The LES model refines the prediction; the mean square deviation error of the results is about 3.5%.Figure 7
Heat flow distribution: mode with a heat-insulated upper boundary (a) and mode with a cooled one (b).
(a)
(b)For the experimental regime with a cooled upper boundary (Figure7(b)), the results have a higher error. All models overestimate the heat flow in the upper part of the calculation area. The root-mean-square deviation of the EARSM and SST results is 11% and 13%, respectively, but the LES as a whole is closer to the experiment and the error was about 6%. The results using the SST and EARSM models relate giving approximately the same result, but the EARSM model still slightly clarifies the distribution of the heat flux on the boundary.The main integral indicator for tasks with internal energy release and heat dissipation from boundaries is the distribution of the boundary heat transfer, characterized by the Nusselt number. Based on the results of recent experimental studies, the following dependence of the Nu number on the modified Ra number was obtained (the dependence was established for Ra ~1013–1017 [3]):(6)Nu=0.385Rai0.233.In Figure8, in accordance with Table 1, the results of numerical simulation of three experimental regimes are presented. For all variants, the dependence of the Nusselt number on the modified Rayleigh number for various approaches to the modeling of turbulent flows is presented.Figure 8
Dependence of the Nusselt number from the modified Rayleigh number.Analyzing the results, we can assume that a sufficiently large deviation is associated with a understating of the heat flux at the upper boundary, while at the cylindrical boundary the heat flux is overestimated. The heat dissipation from the upper boundary is less intensive: in comparison with the experimental data, on the elliptical boundary it is overestimated by an average of 30%. The maximum deviation from the experimental dependence reaches 33%; for the LES, the maximum deviation was 25%. Thus, the use of the LES approach refines the forecasting, but the error remains substantial, which may be due to the overestimation of the heat flow through the upper boundary.In the framework of the COPO experiment, a convective flow (with Rayleigh number ~1015 [2–5]) is considered. Since the currents accompanying accidents in reactor installations are characterized by Rayleigh numbers of even higher orders, considerable interest is in the study of convection at Rayleigh numbers ~1017. The experiment that covers the region most important for the safety problems of nuclear reactors is the BALI experiment.


## 3.3. Modeling the BALI Experiment

Experiments of BALI were carried out in order to investigate the convection of the fuel at a serious accident of the PWR reactor in the second half of the 1990s. In contrast to the COPO experiment, BALI experiments simulate the processes occurring during a severe accident in a reactor with a spherical bottom. Four series of experiments were performed under different conditions and configurations on a 1 : 1 scale. The general view of the experimental setup is shown in Figure9 and represents a cutout along the axial plane [6]. The liquid simulating the melt was a solution of salt in water. The volumetric heating provided by electrodes was distributed throughout the region. The side walls are thermally insulated; the cylindrical boundary was cooled by liquid nitrogen through the interlayer of the heat insulator, which led to the formation of a relatively thin crust of ice and provided isothermal boundary conditions. In the first series of experiments, the upper boundary was thermally insulated; in the second series of experiments, it was cooled similarly to the cylindrical boundary. The Nusselt number (4) is an indicator of the effective heat transfer across the volume boundary.Figure 9
BALI installation scheme.The calculated region corresponded to the complete model of the experimental setup. To simulate the BALI experiment, a block-structured grid was constructed, consisting of 4.2 million cells of a hexahedral shape, the cell size in the main region 0.006 m thick with concentration to all walls, and the size of the wall cell 0.0001 m. As a space-sampling scheme and in the pulse-conservation equation for the LES model, we used the BCD scheme with the Smagorinsky constant. The time step was chosen in such a way as to ensure the value of the Courant number of the order of unity. The thermophilically properties of the coolant used in modeling the experiments are given by temperature-dependent (2), and the coefficient of variable temperature expansion, depending on the temperature, is determined by relation (5).It should be noted that in [6] the results of the experiment are given in dimensionless form in relative units and there are no exact data on volumetric energy release, cooling conditions, and about ice crust on cooled surfaces. It was shown in [6, 7] that the heat flux profile obtained in various experiments is universal and does not depend on the Rayleigh number; however, the shape of the basin and the aspect ratio greatly influence the heat transfer distribution [1, 6].The output parameters for the heat flux and temperature are given in a dimensionless form. The physical picture for the first calculation case, with the thermally insulated top wall, is similar to the COPO problem, the temperature field is clearly stratified, the liquid with the higher temperature (minimum density) is distributed at the top of the region, there is practically no disturbance in the flow, and the greatest value of the velocity arises near the cooled cylindrical wall.In calculations of the current with the adiabatic upper boundary, volumetric energy release was varied, it was selected in such a way as to provide the following values of the Rayleigh criterion:Rai = 2 × 1016 (Ra1), Rai = 3 × 1016 (Ra2). Figure 10(a) shows the distribution of the dimensionless temperature along the vertical boundary and the distribution of the dimensionless heat flux along the cylindrical boundary. It can be seen from the graphs that the profiles of temperature and density of the heat flux are universal and almost independent of the Rayleigh number. The results of the first series of numerical experiments show that the EARSM and SST turbulence models accurately reproduced the general character of convective heat transfer in the experimental BALI setup. The distribution of the dimensionless temperature on the vertical wall obtained during the calculation agrees with the experimental data for all turbulence models and for various Rayleigh numbers. The distribution of the dimensionless heat flux density on the cooled cylindrical boundary is in good agreement with the experimental data in the lower part of the experimental setup, and its magnitude is underestimated in its upper part.Figure 10
The distribution of the time-averaged dimensionless temperature along the vertical boundary (a) and the dimensionless density of the heat flux at the cylindrical boundary (b).
(a)
(b)The SST and EARSM models represented a similar result of the distribution of heat flow and temperature; the root-mean-square error for the heat flux is about 6%. The use of EARSM gives a slight advantage, reducing the error by 1%. The LES model refines the prediction; the rms error of the results is about 4%.The results of the numerical simulation of the problem of BALI for the regime with a cooled upper boundary are below. Several flow zones are observed: the upper unstable layer of a homogeneous temperature, where the Rayleigh-Boehner convection zone occupies the most part; a lower zone with a temperature stratification at which the liquid rises to the center at a low velocity; boundary layer where the liquid in contact with ice cools and returned at a high speed to the bottom of the basin. A similar picture is observed in numerical simulation; a large part of the domain is occupied by a zone of large-scale convection where large vortex structures are localized; heat exchange occurs in the regime of highly developed turbulent convection. For this mode, Figure11 shows the instantaneous velocity field and the instantaneous distribution of the dimensionless temperature, with a Rayleigh criterion of 3.5 × 1016 for the LES turbulence simulation.Figure 11
Instantaneous temperature field (a) and instantaneous velocity field (b).
(a)
(b)The stratified temperature distribution is observed only in the lower part of the calculation area. The greater part of the upper domain is occupied by the convection zone of Rayleigh-Bener and it is much larger than the zone observed in the COPO simulation. This may be due to both the geometry of the installation and the increased value of the number Rayleigh.The experiment contains data on the distribution of temperature and heat flux for the flow regime with cooling walls. The results are presented for the averaged density of the heat flux (Figure12(b)) and the distribution of the dimensionless temperature (Figure 12(a)).Figure 12
The distribution of the time-averaged temperature along the vertical boundary (a) and the distribution of the heat flux density at the cylindrical boundary (b).
(a)
(b)Analysis of the calculation results shows that the RANS models underestimate the heat flow in the upper part of the calculation area. The LES results are closer to the experiment. The root-mean-square deviation of the EARSM and SST results is 11% and 13%, respectively. The LES results are 6%. The results using the SST and EARSM models correlated, giving approximately the same result, but the EARSM model lies closer to the experimental data.Figure13 shows the results of numerical simulation of three experimental regimes with a cooled upper wall for the numbers Ra1 = 6.82 × 1015, Ra2 = 3.16 × 1016, and Ra3 = 1.07 × 1017 in comparison with the obtained Nusselt number dependence on the modified Rayleigh number (6) for EARSM and LES turbulence models.Figure 13
Dependencies of the number Nu on the modified number Ra.The maximum error for the Nusselt number at the upper boundary was about 36%, the vortex-resolving LES Smagorinsky model refines the prediction, and the maximum error using this model is about 33% in the case of the maximum Rayleigh number. Numerical modeling qualitatively corresponds to the physical picture of natural convection. Due to the high-temperature difference, large vortex flows occur at the upper wall. In the lower part of the basin, there is a stable zone of temperature stratification. This picture corresponds to the flow in the experimental setup. The distribution of the dimensionless temperature at the vertical boundary is described with sufficient accuracy, with the greatest deviations being observed in the central area. The calculated density distribution of the dimensionless heat flux at the cooled upper boundary is described with some deviation from the experimental data.


## 4. Conclusion

Numerical simulation of the processes of turbulent natural convection is demonstrated for different values of the Rayleigh number. The application of different approaches to the description of turbulence is provided. Analysis of the experimental data and the results of calculations show that the qualitative picture of the flow with the values of the Rayleigh criterion to 1017 is modeled correctly. In the case that modeling convective heat transfer in a cubic cavity at Ra~109 coincides with the experimental data, the qualitative and quantitative parameters are coinciding.It is shown that the distribution of the dimensionless temperature and the heat flux density along the boundaries does not depend on the turbulence model and is not dependent on Rayleigh criterion. This is in good agreement with the experimental data, but the quantitative characteristics for some regimes have significant differences. For experiments with a heat-generating liquid (COPO and BALI) with a heat-insulated upper boundary, the effect of applying approaches to turbulence modeling is not particularly significant. This may be due to a less intense internal flow, since heat exchange occurs with only one boundary. Within the framework of these regimes, the results of numerical simulation are in good agreement with the results of experimental studies. For problems with two cooled walls, where there is intense circulation of liquid in the closed region, the value of the heat flux at the lower boundary is underestimated and the heat transfer at the upper boundary is increased. In various sources [3, 19, 20], a high error in the heat transfer value is also observed in the simulation of the convection. This error is increased with an increase in the Rayleigh number. This may be due to the increasing influence of the assumptions made during the calculation, since the set of initial data is insufficient to accurately repeat the experiment using CFD modeling. Numerical modeling requires multiple assumptions. These assumptions include differences in the specification of fluid properties, which affects the calculated integral characteristics. A uniform volumetric energy release is accepted, although in the experiment it was inhomogeneous [20]. In addition, there is no data on the concentration of the salt solution in water (the properties of pure water are used). The lack of accurate data on the cooling conditions at the boundaries, as well as data on the geometric dimensions of the formed ice crust, can introduce a high error in the final result. The modeling of the ice crust is a serious problem and does not allow taking into account the surface roughness, which can strongly affect the heat transfer coefficient. The rough structure contributes to an additional flow turbulence and intensification of heat transfer, which is not taken into account when performing numerical experiments.The LOGOS software at a sufficiently high level reproduces the heat exchange in the heat-generating fluid with natural convection under conditions of high Rayleigh numbersRai~1017. An important role is played by the choice of the approach to modeling turbulence. Thus, for regimes with a cooled upper limit, the intensive application of the vortex-resolving LES model both qualitatively and quantitatively clarifies the description of the processes of natural convection and heat and mass transfer.



---

*Source: 5781602-2018-01-09.xml*