Selection of a turbulence model to calculate the temperature profile near the surface of VVER-1000 fuel assemblies in the NPP spent fuel pool*

The paper considers the problem of selecting a turbulence model to simulate natural convection near the surface of a VVER-1000 fuel assembly unloaded from the reactor by computational fluid dynamics (CFD simulation) methods. The turbulence model is selected by comparing the calculated data obtained using the Ansys Fluent software package with the results of experimental studies on the natural convection near the surface of a heated vertical plate immersed in water, which simulates the side face of the VVER-1000 fuel assembly in a first approximation. Two-parameter semi-empirical models of turbulence, k-ε and k-ω, are considered as those most commonly used in engineering design. The calculated and experimental data were compared based on the excessive temperature of the plate surface and the water temperature profiles in the turbulent boundary layer for convection modes with a Rayleigh number of 8∙1013 to 3.28∙1014. It has been shown that the best agreement with experimental data, with an average deviation not exceeding ~ 8%, is provided by the RNG k-ε model which is recommended to be used to simulate natural convection near the surface of VVER-1000 FAs in the NPP spent fuel storage pool.


Introduction
A long-term experience in testing fuel assemblies (FA) for pressurized-water power reactors (PWR and BWR) on inspection benches in the NPP spent fuel storage pools has proved demonstratively its efficiency in such missions as providing the R&D basis for the FA operation, introducing new FA and fuel element designs, and justifying new fuel cycles (Pavlov 2013(Pavlov , 2014. The efficiency of pool inspections is defined primarily by the timely acquisition of information right after the FA withdrawal from the reactor core, and by the possibility of testing more FAs involving much smaller costs than through tests in safety chambers at research centers (Pavlov 2014).
A repair and inspection bench for the VVER-1000 reactor FAs has been developed and is deployed at NPPs which makes it possible to examine closely the FAs, identify their form change, detect and remove leaky fuel elements, and perform eddy-current radioscopy of the fuel cladding (Ivanov et al. 2017). Essential in improving further the procedural guidance for inspection benches and in developing new nondestructive inspection techniques are the characteristics of tested FAs capable to affect the performance of sensors or the accuracy of measurements. Such characteristics include the FA radiological parameters (gamma and neutron dose rates) and thermophysical parameters (temperature of the fuel element and FA components).
Right after the VVER-1000 FA is withdrawn from the reactor core, its decay heat is rather high (reaching hundreds of kilowatts) due to which natural convection of water develops near the FA surface in the spent fuel pool. The convective layer is characterized by the thickness and the water temperature and velocity distribution in it, the convection parameters varying lengthwise the FA in a range from a laminar mode in the assembly's bottom part to full turbulence at the top (Voronina and Pavlov 2019).
The convective layer formed near the FA surface with a noticeable water temperature gradient should be taken into account, specifically when using ultrasonic pulse-echo techniques to measure the FA dimensions (Pavlov et al. 1991). To estimate the accuracy of measuring the distance between the sensors and the FA surface, bearing in mind that the ultrasonic wave propagation velocity in water depends on temperature, one needs to know how the sound velocity changes in water along the ultrasonic wave propagation path and, accordingly, how the temperature changes in the boundary layer.
Numerical simulation based on CFD (Computational Fluid Dynamics) techniques can be used to estimate the parameters of the convection near the FA surface (Raman et al. 2018). In a general form, the convection near the FA surface is described by Navier-Stokes equations (McLean 2013), which are solved in engineering design by the RANS (Reynolds-Averaged Navier-Stokes) method (McLean 2013). The RANS method suggests that Navier-Stokes equations are represented as a system of Reynolds-averaged equations with different semi-empirical turbulence models used to close it (McLean 2013).
Hundreds of turbulence models have currently been developed of which none is however versatile. When there is a new item to be numerically simulated for convection, a problem therefore arises as to what particular turbulence model to select and how to estimate the confidence of the results obtained by simulation.
As a rule, the rationale for the turbulence model selection is provided through comparing the calculated data obtained using different models with experimental results. And the experiment can be conducted using the tested item as such or its model or the item's convection processes which are physically close (identical) to same.
Literature does not contain experimental data on the parameters of natural convection near the surface of a VVER-1000 FA with decay heat immersed in water. This leads to the need to search for examples of experimental studies on natural convection for items which are close, as far as possible, to the heated VVER-1000 FA in terms of thermal physics. The convective layer along the FA's outer surface results primarily from deposition of heat from the fuel elements in the peripheral row which can be simulated in a first approximation by a heated vertical plate with the width equal to the VVER-1000 FA flat (~ 135 mm) and the height equal to the fuel element height (~ 4000 mm), and with the heat flux equal to the flux from the fuel elements and directed into the environment perpendicularly to the plate. In other words, the VVER-1000 FA is simulated by a heated hexagon.
The paper presents the results of a CFD simulation using the Ansys Fluent software package (ANSYS Fluent) for the natural convection along a vertical heated plate immersed in water. Two-parameter semi-empirical models of the k-ε and k-ω families were used as the turbulence models. The calculation results were compared with experimental results (Vliet andLiu 1969, Qureshi andGebhart 1978).

CFD simulation of convection
There is a great number of theoretical and experimental studies dealing with investigating the natural convection near a vertical heated surface. Since the phenomenon of natural convection is of most interest to scientists in calculation of structural elements and utility systems regarding the removal or supply of heat, many papers are concerned with studying the natural convection near a vertical surface uniformly heated in an air environment. Studies to investigate convection near a non-isothermal surface are much fewer. One of the earliest papers with published results from experiments to study laminar convection near a vertical surface with a constant heat flux in water is that by Lock and Trotter (Lock and Trotter 1968). Vliet and Liu (Vliet and Liu 1969) continued investigating natural convection near an immersed plate with a constant heat flux. Experiments were conducted for different water temperature and heat flux values. Similar studies were undertaken by Gebhart (Qureshi and Gebhart 1978).
The experiments published in (Vliet andLiu 1969, Qureshi andGebhart 1978) were selected to compare the results obtained by CFD simulation with the experimental data. These studies provide detailed data on the excessive wall temperature and the temperature profile in the boundary layer.
In a general form, the flow of a viscous liquid is described by a system of Navier-Stokes equations consisting of a mass conservation equation and the law of conservation of momentum (Lecheler 2014): where v is the liquid velocity vector; μ is the dynamic viscosity; F is the vector of volume forces; р is static pres-sure; ρ is the density;▽ is the Hamilton operator; and Δ is the Laplace operator.
There are three approaches to simulation of turbulent flows: • direct numerical simulation (DNS); • solution of systems of Reynolds-averaged equations (RANS); • large eddy simulation (LES) method.
DNS suggests solving complete unsteady-state Navier-Stokes equations. This method is hard to use to simulate actual physical flows, since it requires computational resources expected to become available by only 2080 (Spalart 2000). LES requires a quality description of the turbulent flow with a smaller power as compared with the previous method and makes it possible to simulate not only minor eddies but also to calculate directly large eddies (Belov 2001). Due to LES requiring the use of parallel computation systems, the RANS method was used for turbulence simulation. Many RANS turbulence models have been developed for closing Reynolds-averaged equations. This method allows one to obtain a physically correct result with a relatively small resource capacity (Spalart et al. 2014).
Many studies were undertaken to simulate a free convective flow near a vertical heating surface for isothermal conditions (Bagayev and Syraleva 2018, Yingchun 2014, Ziganshin et al. 2012. Various turbulence models were used to close system (1). These studies have shown that the most common turbulence models, k-ε and k-ω, out of the family of two-parameter RANS models can be used to calculate natural convective flows and heat exchange. It is however required to assess their applicability for investigating a free convective flow near a non-isothermal surface.
The problem was simulated numerically in an unsteady two-dimensional statement. A flow was studied along a heat emitting plate of the height 1.5 m immersed in a sta-tic liquid (water) with the constant temperature T ∞ . The diagram of the computational region is presented in Fig. 1.
The initial conditions have the following form: u(x, y, 0) = v (x, y, 0) = k(x, y, 0) = ε(x, y, 0) = ω(x, y, 0) = 0; where u, v are the projections of the velocity vector on the x and y axes respectively; x, y are the coordinates, m; k is the kinetic energy of turbulence, m 2 /s 2 ; ε is the dissipation rate of the turbulence kinetic energy, m 2 /s 3 ; ω is the specific dissipation rate of the turbulence energy, s -1 ; and T ∞ is the water temperature, °С. Temperature conditions of the second kind were defined as the boundary conditions on the plate surface FE T q w O wn as well as no-slip and impermeability conditions for the velocity v = 0, where λ is the wall material's thermal conductivity; n is the vector of the normal to the surface; q is the heat flux density, W/m 2 ; and v is the velocity vector.
A symmetry condition meaning the equality to zero of the boundary-normal liquid motion velocity component and of the temperature and pressure gradients has been defined for areas AE and FB: v(x, y, t) = 0, The following boundary conditions are defined on boundaries AD, CD and BC for velocity, pressure and temperature: where p ∞ and T ∞ are the water pressure and temperature respectively.
The following k-ε and k-ω turbulence models implemented in Ansys Fluent were used for the turbulence simulation: Standard k-ε, RNG k-ε, Realizable k-ε, Standard k-ω, and SST k-ω. Due to a limited size, the paper does not describe in detail each of the models. The mathematical formulation and the peculiarities of the models are presented in (Lecheler 2014).
Models of the k-ε family are not capable to simulate the effects taking place in the near-wall region. Semi-empirical formulas, referred to as wall functions, are used for the quality simulation of the flow in the boundary layer using k-ε models. Ansys Fluent has the following models for near-wall turbulence: Standard Wall Functions, Non-equilibrium Wall Functions, and Enhanced Wall Treatment (ANSYS Fluent). Enhanced Wall Treatment was used in the study. The use of the Standard Wall Function and the Non-equilibrium Wall Function has shown a poor agreement in the simulation of a similar phenomenon (convection near an isothermal plate) (Ziganshin et al. 2012).
Boussinesq approximation is used not infrequently to describe convective motion (Bykalyuk et al. 2015). This approximation can be however used only with minor temperature drops (ΔT < 2 °С for water) (Ferziger and Peric 2002). In this connection, all thermophysical properties of water, including density, are assumed to be polynomial functions of temperature.
An equation solver, called Pressure-Based Solver, was used for solving. The Simple (Semi-Implicit Method for Pressure-Linked Equations) algorithm was used to calculate the velocity field coupling with pressure involving a counter-flow pattern of the second order accuracy (Second Order Up-wind) for the convective members in the momentum conservation equation, for the kinetic turbulent energy equation, and for the turbulent energy dissipation equation. The convergence of the solution was regarded to be achieved in events when the difference between iterations in solving continuity, momentum, energy, turbulent kinetic energy, and dissipation rate equations reached 1×10 -6 .
Regular meshes condensing towards the heated wall were used for the calculations. Ansys Meshing was used as the mesh generator. The mesh model was built with regard for the key criteria that characterize its quality (ANSYS Fluent).

Simulation results and comparison with experimental data
Convection was simulated for the conditions of two experiments the results of which are published in (Vliet and Liu 1969), and for one experiment in (Qureshi and Gebhart 1978). In experiment 1 in (Vliet and Liu 1969), the heat flux from the plate was q = 19497 W/m 2 ; the water temperature was T ∞ = 33 °С; and the temperature profile in the boundary layer was recorded at a height of h = 1.067 m from the plate's bottom edge. The conditions of experiment 2 in (Vliet and Liu 1969) were as follows: q = 28661 W/m 2 ; T ∞ = 23 °С; h = 0.762 m, and those of the experiment in (Qureshi and Gebhart 1978) were q = 4488 W/m 2 , T ∞ = 23 °С; h = 1.194 m.
Figs 2, 3 present variations in the excessive wall temperature ΔT = T w -T ∞ (where T w is the plate surface temperature) along the plate height for experiments 1 and 2 in (Vliet and Liu 1969). The horizontal lines in the figures show the elevations for the end of the convective laminar flow with Ra = 3•10 12 and for the full turbulence onset, Ra = 4•10 13 .
It can be seen from the figures that the results for all of the turbulence models used are close to the experimental results in the laminar region and predict satisfactorily the laminar-turbulent transition region (3•10 12 < Ra < 4•10 13 ).  In the developed turbulence region, Ra > 4•10 13 , the results of the calculations based on the Standard k-ω model stand out sharply against the overall picture giving the plate wall temperature overestimated by nearly 10 °C (Fig. 3). This model was not therefore further used to calculate the water temperature profile in the boundary layer. Fig. 4 through 6 present the numerical and experimental investigation results for the temperature profile in a boundary turbulent layer. Basically, all four remaining turbulence models describe satisfactorily the temperature profile for the three experiments under consideration.
The greatest discrepancies between the experimental and calculated data fall on the region near the plate wall being in the limits of ~ 0.1 mm. Here, the calculated data exceeds that obtained in the experiment. The possible reason for such discrepancies is that the water temperature is hard to measure correctly in the immediate vicinity of the plate, since the temperature sensors perturb the structure of the convective flow viscous sublayer. Table 1 presents the results of calculating the deviation of the numerical simulation results obtained using different turbulence models from the experimental data in Fig. 4 through 6. No experimental data obtained nearer to the plate surface than 0.1 mm was used to calculate the deviations.
It can be seen from the table that the best result is achieved in the event of the RNG k-ε model with the average deviation from the experimental data being ~ 8%, as the Standard k-ε, Realizable k-ε, and SST k-ω models give 11.6, 14.7, and 14.3% respectively.

Conclusions
The Ansys Fluent code was used for the CFD simulation of the natural convection near a submerged vertical heated plate with a constant heat flux density. Five two-parameter semi-empirical turbulence models, most often applied in engineering design, were used for the simulation, including Standard k-ε, RNG k-ε, Realizable k-ε, Standard k-ω, and SST k-ω.
The simulation results were compared with the experimental data obtained earlier for the excessive temperature of the plate surface and the water temperature profiles in the turbulent boundary layer for the three heat flux density values: ~ 4.5, 19.5, and 28.7 kW/m 2 . The comparison results have shown the following.
All turbulence models used give the calculation results which are close to the experimental values in the laminar region and predict satisfactorily the laminar-turbulent transition region.
To the exception of Standard k-ω, all models under consideration calculate satisfactorily the plate surface temperature in the developed turbulence area. As far as the temperature profile is concerned, the best result is provided by the RNG k-ε model with the average deviation from the experimental data being ~ 8%.
Therefore, it is recommended to use a two-parameter semi-empirical turbulence model, RNG k-ε, as the turbu-    lence model to calculate, using the Ansys Fluent code, the temperature profile in the laminar and turbulent boundary layers of natural convection near the outer surface of the irradiated VVER-1000 FA immersed in water.