Corresponding author: Pavel V. Amosov ( p.amosov@ksc.ru ) Academic editor: Yury Kazansky
© 2021 Pavel V. Amosov.
This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Citation:
Amosov PV (2021) Forecast of the thermal regime of an underground storage facility for heatgenerating materials under mixed convection conditions. Nuclear Energy and Technology 7(1): 18. https://doi.org/10.3897/nucet.7.64361

The paper presents the results of a study based on numerical simulation methods of the thermal regime of an underground facility for longterm spent nuclear fuel storage in the version of a builtin reinforced concrete structure. A multiphysical computer model was constructed in a twodimensional setting by means of the COMSOL software. The mathematical model was based on the continuity equations, NavierStokes equations and the general heat transfer equation. The conditions of mixed convection were taken into account in the ‘incompressible ideal gas’ approximation, in which the thermophysical properties of air were a function of temperature only. For two parameters of the model, the following values were taken: the air flow rates providing the velocity at the inflow boundary = 0.01, 0.03 and 0.05 m/s, and the effective heat conductivity coefficients of the material of the builtin structure = 1.0 and 2.0 W/(m×K). Numerical experiments were performed for a period of up to 5 years of fuel storage. Special emphasis was given to the fundamental difference between the nonstationary structure of the velocity fields forecasted in the model of an ‘incompressible ideal gas’ and the ‘frozen’ picture of aerodynamic parameters in the model of an incompressible fluid. An analysis was made of the dynamics of spatial temperature field distributions in different areas of the model. It was shown that the criterion temperature control requirements were met when the facility was operated under conservative ventilation conditions in terms of the air flow rate and the heat conductivity coefficient of the builtin structure material. The dynamics of heat flows directed into the rock mass through the base and from the surface of the builtin structure into the air was analyzed. The heat flow dominance from the structure surface was also noted. Finally, the influence of the effective heat conductivity coefficient of the builtin structure material and the air flow rate on the values of heat flows directed into the air and rock mass was demonstrated.
Heatproducing materials, spent nuclear fuel, storage facility, numerical simulation, mixed convection, thermal regime
One of the results of research carried out by specialists of the Mining Institute of the KSC RAS on the problem of safe management of spent nuclear fuel (SNF) from ship power plants in the European North of Russia was the development of a conceptual designlayout diagram of an underground SNF storage facility (
Forecast calculations of the thermal state of materials of an underground SNF storage facility for such a long period are usually performed on the basis of numerical simulation methods. There is a number of modern verified computer programs, e.g., COMSOL, FlowVision or ANSYS (COMSOL Documentation (accessed Jan 20, 2020),
Correcting this point, in view of the experience accumulated over the past years in simulating aerothermodynamic processes (
The focus is on a module, which is a working of chambertype for storing spent nuclear fuel in a builtin reinforced concrete structure (
In view of the convective and conductive heat transfer mechanisms, the equation is written as follows (COMSOL Documentation (accessed Jan 20, 2020),
$\rho \xb7{C}_{p}(\partial T/\partial t)\nabla \left(k\xb7\nabla +\rho \xb7{C}_{p}\xb7T\xb7\mathbf{u}\right)={q}_{w}$, (1)
where u is the air velocity vector; ρ is the density; C_{p} is the specific heat capacitance at constant pressure; k is the heat conductivity coefficient; q_{w} is the heat source or sink; T is the temperature; t is the time; ∇ is the Hamiltonian operator.
The nonisothermal flow regime is simulated using the continuity and NavierStokes equations, which describe the relationship between the air velocity and pressure (COMSOL Documentation (accessed Jan 20, 2020),
$\rho (\nabla \mathbf{u})\mathbf{u}=\nabla \left[p\mathbf{I}+\eta \left(\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^{T}\right)(2\eta /3k)(\nabla \mathbf{u})\mathbf{I}\right]+\left(\rho {\rho}_{0}\right)\mathbf{g},\nabla \left(\rho \mathbf{u}\right)=0$ (2)
where η is the dynamic viscosity; g is the gravity vector; p is the pressure; I is the identity matrix. Air heating causes deviations in the local density ρ compared to the inflow air density ρ_{0}. The result is a buoyancy force expressed as (ρ – ρ_{0})g. In addition to the wellknown approach used to solve such problems (the socalled Boussinesq approximation), in a number of software products (ANSYS FLUENT, COMSOL, FlowVision), it is possible to use the ‘incompressible ideal gas’ model. In this model, some properties of air are considered as dependent solely on temperature. For example, COMSOL users are advised to use the following ratios(
ρ = p_{0}×μ/(R×T),
lg k = 0,865×lg T – 3,723, (3)
η = 6,0×10^{–6} + 4,0×10^{–8.}T,
where p_{0} = 101325 Pa; μ = 0.0288 kg/mol; R = 8.314 J/(mol×K).
A detailed description of the variables used in nonstationary heat transfer models is given in the COMSOL documentation and reference files (COMSOL Documentation (accessed Jan 20, 2020),
The initial temperatures of the areas of the model and the values of the thermophysical parameters adopted in the calculations (
Initial temperatures T_{0} and physical parameters of the areas of the model.
T_{0}, K  C_{p}, J/(kg×K)  k, W/(m×K)  , kg/m^{3}  , Pa×с  

Air  288  1000  k(T)  (T)  (T) 
Granite  285  740  3  2490  
Heat release area  300  580  1.0 and 2.0.  5860 
In the performed numerical experiments, the values of two parameters were changed ‘manually’, in relation to which the authors of (
Technically, the values of the thermophysical parameters of complex heterogeneous systems can be calculated using the relations given in (
The volumetric power curve of the residual heat release q_{w} as a function of time was adopted as before (
q_{w} = 9.149 – 0.224.t + 0.00277×t^{2} – 1.651×10^{–5}×t^{3},
where t are years.
To implement the model in numerical terms, a triangular computational grid of the ‘extremely fine’ type was used. The total number of elements was more than 6 100 (the number of degrees of freedom was almost 25 000). With the area of the air environment of 657 m^{2} (see Fig.
A highly efficient direct solver UMFPACK was used to solve the combined system of equations (1)–(3). To ensure the stability of the calculations, it was required to involve the Isotropic Diffusion artificial viscosity technique with an adjustment parameter at the level of 0.25–0.35 for the NavierStokes equations and 0.6 for the heat transfer equation. Numerous attempts to obtain stable solutions using other techniques (Streamline Diffusion, Crosswind Diffusion) to combat numerical instability were not successful (COMSOL Documentation (accessed Jan 20, 2020),
The boundary conditions used in the numerical solution of the system of equations (1)–(3) are standard (
Boundary conditions for the numerical solution of the system of equations (1)–(3).
Boundary  Aerodynamic equations  Heat transfer equation 

Inflow  Horizontal component of velocity  Incoming air temperature 
Outflow  No stress  Convective flow 
Storage module walls  No slip  Continuity condition 
External model boundaries  –  Zero flow condition 
The simulation time in all the calculations was five years.
The first fundamental difference in the results of numerical experiments in 2010 (
It is clearly seen that the structure of the velocity fields, calculated in the ‘incompressible ideal gas’ approximation, is significantly different from that used by the authors in the model of an incompressible fluid (
‘Instantaneous’ picture of the structure of the velocity field for fiveyear modeling with variations in the incoming air flow rate (upper 0.06, average 0.18 and lower 0.30 m^{3}/s) and the effective heat conductivity coefficient of the heat release area: a) – 1.0 W / (m×K); b) – 2.0 W / (m×K).
For example, Fig.
The influence of the effective heat conductivity coefficient on the velocity field structure can be conveniently traced in the following examples. For an air flow rate of 0.18 m^{3}/s (middle pictures), Fig.
The velocity field formation dynamics was analyzed using the example of the horizontal component of the velocity vector in three vertical sections of the model (marks 70, 80 and 90 m along the horizontal axis (see Figs
Thus, it becomes obvious that, in the models of an incompressible liquid, on the basis of which the studies of 2010 were carried out (
To answer the question about the conditions for meeting the criterion temperature values, the spacetime temperature distributions for the abovementioned options of the air flow rate and the values of the effective heat conductivity coefficient of the heat release area were analyzed. Let us remember that, first of all, the criterion temperature values are of interest – in the area of the rock mass (≤ 100 °C) and on the surface of the reinforced concrete structure (≤ 85 °C). Analysis of spacetime temperature distributions shows that the search for maximum temperatures in different areas of the model must be performed in different sections (see above):
In the indicated sections, the dynamics of the spatial temperature distribution is analyzed for the option of an air flow rate of 0.18 m^{3}/s and an effective heat conductivity coefficient of the heat release area of 1.0 W/(m×K). The levels of the forecast temperatures in all the areas of the model are fixed (rock mass – the bottom and roof of the working of chambertype, air, builtin reinforced concrete structure).
The most important point of the performed analysis, i.e., the condition for not exceeding the criterion temperature values in the storage facility for unrecycable types of SNF, is confirmed. However, if, in the model of an incompressible fluid, the heating of the rock mass was forecasted to a level of 40 °C and the surface of the engineering structure up to 50 °C (
Analysis of spacetime temperature distributions in different areas of the model shows that the coordinates of points with maximum temperature values change not only in time, which is natural for a nonstationary problem, but also in space when the values of the above parameters (k_{w} and Q) vary.
Let us note several general patterns of behavior of the analyzed graphical dependencies:
Table
Calculated values of the maximum temperatures in different regions of the model with variations in k_{w} and Q.
k_{w}, W/(m×K)  Q, m^{3}/s  Controlled areas of the model  

Heat releases  Surface  Granite_1  Granite_2  Air  
1,0  0,06  90,2  32,4  90,1  33,4  33,4 
0,18  82,7  22,5  82,6  22,3  22,3  
0,30  82,3  21,8  82,2  19,0  19,0  
2,0  0,06  64,2  33,5  64,2  35,2  35,2 
0,18  56,4  24,9  56,4  22,8  22,8  
0,30  55,0  22,3  55,0  19,3  19,3  
Note. Granite_1 – base, Granite_2 – roof 
The maximum heating takes place at the minimum values of the air flow rate and the effective heat conductivity coefficient of the heat release area (conservative conditions). For the estimated time of five years, the forecast temperature of the rock mass reaches 90 °C: this is close to the criterion value. But it should be noted that by this time the heating rate is noticeably reduced. And taking into account the decline in the residual heat generation power, it can be guaranteed that the threshold value will not be exceeded. In addition, there is one more argument in hand for not exceeding the criterion value, i.e., the model geometry used in the calculations is finite. In real conditions, heat outflow will occur in larger volumes of rock, which will further reduce the forecast temperatures.
In addition to the analysis of the temperature distribution, the issue of heat leakage in time from the heat release area was considered. A program option was used that integrates physical parameters either over the boundary or over the volume. In this situation, a point of interest is the heat flow leaving the specified area. An analysis of the heat flows directed through the transverse and longitudinal boundaries of the heat release area showed that the heat flows directed through the lower base and the upper surface of the builtin structure are dominant. The dynamics of these heat flows with variations in both the incoming air flow rate and the effective heat conductivity coefficient of the heat release area is shown in Fig.
Dynamics of the heat power leaving the heat release region through (a) the lower base and (b) the upper surface of the builtin structure: I – 0.06 m^{3}/s, 1.0 W/(m×K); II – 0.06 m^{3}/s, 2.0 W/(м×K); III – 0.18 m^{3}/s, 1.0 W/(m×K); IV– 0.18 m^{3}/s, 2.0 W/(m×K); V – 0.30 m^{3}/s, 1.0 W/(m×K); VI – 0.30 m^{3}/s, 2,0 W/(m×K).
Figures
The qualitative and quantitative indicators of the graphs in Fig.
The behavior of the graphs in Fig.
The graphs in Fig.
Based on the conducted research, the following conclusions can be formulated.
A computer model of an underground facility for longterm spent nuclear fuel storage was constructed in the version of a builtin reinforced concrete structure in a twodimensional setting. The effect of gravity was taken into account in the ‘incompressible ideal gas’ approximation.
A comparative analysis was carried out regarding the calculated aerothermodynamic parameters in the ‘incompressible ideal gas’ approximation with the results of numerical experiments in the model of an incompressible fluid.
The fundamental differences were shown in the structure of the velocity fields forecasted in the storage facility based on the two indicated models, with variations in the incoming air flow rate and the effective heat conductivity coefficient of the builtin structure material.
The maximum temperature values in various areas of the model were analyzed for the considered range of values of the effective heat conductivity coefficient of the heat release area and the flow rate of air supplied to the facility. Consideration was also given to the physical features of the influence of these parameters on the maximum temperatures of the areas of the model.
The dynamics of heat flows directed into the rock mass through the base and from the surface of the builtin structure into the air is analyzed. The heat flow dominance from the structure surface is noted. As the rock mass is gradually heated up, the dominance increases.