Forecast of the thermal regime of an underground storage facility for heat-generating materials under mixed convection conditions

Pavel V. Amosov

doi:10.3897/nucet.7.64361

Forecast of the thermal regime of an underground storage facility for heat-generating materials under mixed convection conditions*

Pavel V. Amosov^‡

‡ Institute of North Industrial Ecology Problems, Apatity, Russia

Corresponding author: Pavel V. Amosov ( p.amosov@ksc.ru )

Academic editor: Yury Kazansky

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 heat-generating materials under mixed convection conditions. Nuclear Energy and Technology 7(1): 1-8. https://doi.org/10.3897/nucet.7.64361

Abstract

The paper presents the results of a study based on numerical simulation methods of the thermal regime of an underground facility for long-term spent nuclear fuel storage in the version of a built-in reinforced concrete structure. A multiphysical computer model was constructed in a two-dimensional setting by means of the COMSOL software. The mathematical model was based on the continuity equations, Navier-Stokes 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 built-in 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 non-stationary 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 built-in structure material. The dynamics of heat flows directed into the rock mass through the base and from the surface of the built-in 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 built-in structure material and the air flow rate on the values of heat flows directed into the air and rock mass was demonstrated.

Keywords

Heat-producing materials, spent nuclear fuel, storage facility, numerical simulation, mixed convection, thermal regime

Introduction

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 design-layout diagram of an underground SNF storage facility (Melnikov et al. 2003, 2006). The facility was proposed to be created in a rock mass at a depth of 100 m. The SNF storage time was assumed to be about 50 years. Along with nuclear and radiation safety, one of the main factors determining the operating conditions of the SNF storage facility is the temperature regime of the irradiated fuel storage module (NP Safety Requirements-035-02 2002, Nuclear Energy Facilities NP-061-05 2005, Kalinkin et al. 2009). For a dry SNF storage facility, along with general requirements for all types of SNF storage facilities, there is a requirement to ensure a safe fuel storage temperature regime (NP Safety Requirements-035-02 2002, Kalinkin et al. 2009). According to the temperature control requirements, it is necessary to remove residual heat release from spent nuclear fuel, which must be organized so as to exclude the possibility of overheating (358 K (85 °C) for the surface of the built-in structure and 373 K (100 °C) for the rock mass (Kalinkin et al. 2009; Naumov and Gusak 2019; Rzhevsky and Novik 1978)).

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), ANSYS 2016, Shchelyaev (accessed Aug 21, 2018) (a full description in English is available for users from the Help menu), which can be used in such evaluations. The results of the studies on the thermal regime of the underground SNF storage facility in versions of either a built-in reinforced concrete structure or reinforced concrete containers and natural (rock) protective barriers are presented in (Amosov and Podshivalova 2010a, b). In these studies, the thermal regime of the facility was simulated using the COMSOL software (COMSOL Documentation (accessed Jan 20, 2020), Yegorov 2006). The mathematical model was based on the continuity equations, Navier-Stokes equations for a viscous incompressible fluid and the general heat transfer equation. The implemented formulation lacked an important term in the Navier-Stokes equations, which took into account the effect of changes in air density with temperature, which undoubtedly reduced the value of the studies performed. As noted by the authors of (Amosov and Podshivalova 2010b), their attempts to include in the Navier-Stokes equations the effect of gravity in the Boussinesq approximation for the initial residual heat release value indicated below were unsuccessful.

Correcting this point, in view of the experience accumulated over the past years in simulating aerothermodynamic processes (Amosov and Novozhilova 2014, Amosov et al. 2018), this article presents the results of a study on the thermal regime and aerothermodynamics of the atmosphere of the underground facility for long-term SNF storage (the option of a built-in reinforced concrete structure) under mixed convection conditions.

Model geometry, governing equations and input data

The focus is on a module, which is a working of chamber-type for storing spent nuclear fuel in a built-in reinforced concrete structure (Amosov and Podshivalova 2010a). The geometric parameters required for constructing a two-dimensional model are shown in Fig. 1. Areas (built-in reinforced concrete structure, working air and enclosing rock mass) are identified, which differ both in the heat transfer mechanisms and in the thermophysical parameters.

Figure 1.

Model geometry of the underground spent nuclear fuel storage facility (dimensions are in meters) (Amosov and Podshivalova 2010a)

In view of the convective and conductive heat transfer mechanisms, the equation is written as follows (COMSOL Documentation (accessed Jan 20, 2020), Yegorov 2006):

$ρ \cdot C_{p} (\partial T / \partial t) - \nabla (- k \cdot \nabla + ρ \cdot C_{p} \cdot T \cdot u) = 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 non-isothermal flow regime is simulated using the continuity and Navier-Stokes equations, which describe the relationship between the air velocity and pressure (COMSOL Documentation (accessed Jan 20, 2020), Yegorov 2006):

$ρ (\nabla u) u = \nabla [- p I + η (\nabla u + (\nabla u)^{T}) - (2 η / 3 - k) (\nabla u) I] + (ρ - ρ_{0}) g, \nabla (ρ u) = 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 ρ₀. The result is a buoyancy force expressed as (ρ – ρ₀)g. In addition to the well-known approach used to solve such problems (the so-called 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(Yegorov 2006):

ρ = p₀×μ/(R×T),

lg k = 0,865×lg T – 3,723, (3)

η = 6,0×10^–6 + 4,0×10^–8.T,

where p₀ = 101325 Pa; μ = 0.0288 kg/mol; R = 8.314 J/(mol×K).

A detailed description of the variables used in non-stationary heat transfer models is given in the COMSOL documentation and reference files (COMSOL Documentation (accessed Jan 20, 2020), Yegorov 2006).

The initial temperatures of the areas of the model and the values of the thermophysical parameters adopted in the calculations (Amosov et al. 2018; Earth Science 1969; Vargaftik 1972; Pekhovich and Zhidkikh 1968; Melnikov et al. 1975; Dortman 1984; Horai 1971; Chirkin 1968; Melnikov et al. 2010) are shown in Tab. 1.

Table 1.

Download as

CSV

XLSX

Initial temperatures T₀ and physical parameters of the areas of the model.

	T₀, K	C_p, J/(kg×K)	k, W/(m×K)	, kg/m³	, 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 (Amosov and Podshivalova 2010a) had certain justified doubts: (1) the effective heat conductivity coefficient of the heat release area (k_w), consisting of several dissimilar materials with a large spread in the values of this parameter (SNF, reinforced concrete, air); and (2) the supplied air flow rate (Q).

Technically, the values of the thermophysical parameters of complex heterogeneous systems can be calculated using the relations given in (Pekhovich and Zhidkikh 1968, Melnikov et al. 1975, Melnikov et al. 2010). This is exactly what the authors of publications (Amosov and Podshivalova 2010a, b) did with respect to the density and heat capacity of the heat release area. With account taken of the earlier conclusion about the unconditional criterion temperature observance at high values of the effective heat conductivity coefficient in the heat release area (Amosov and Podshivalova 2010a), the values of this parameter in the area of the built-in structure in numerical experiments were taken as 1.0 and 2.0 W/(m×K). For the flow rate of the air supplied to the facility, the values of 0.06, 0.18 and 0.30 m³/s were used.

The volumetric power curve of the residual heat release q_w as a function of time was adopted as before (Amosov and Podshivalova 2010a):

q_w = 9.149 – 0.224.t + 0.00277×t² – 1.651×10^–5×t³,

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² (see Fig. 1), the number of triangular grid elements was 2609, which corresponds to the average value of the triangle side of about 0.75 m. To solve the non-stationary problem, the time option FREE (COMSOL Documentation (accessed Jan 20, 2020), Yegorov 2006) was used with the lower boundary of the interval time step equal to 3.15 s, which ensures that the Courant criterion is met with a margin at the specified air flow rates. Note that stable calculations turned out to be those with the value of the lower boundary of the time step interval by an order of magnitude larger. At the same time, the differences in the values of the maximum temperature at the same moments of the calculated time were insignificant (the third decimal place).

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 Navier-Stokes 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), Yegorov 2006).

The boundary conditions used in the numerical solution of the system of equations (1)–(3) are standard (Yegorov 2006) and presented in Tab. 2.

Table 2.

Download as

CSV

XLSX

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.

Results of numerical experiments and their analysis

Aerodynamic flow characteristics

The first fundamental difference in the results of numerical experiments in 2010 (Amosov and Podshivalova 2010a, b) and presented in the article lies in the structure of the air flows of the storage facility. If in the model of an incompressible fluid (Fig. 2), the structure of the air flows at the considered air flow rates is identical throughout the solution of the non-stationary heat transfer problem (the difference is only in the values of the velocity components for different flow rates), then the effect of heat on aerodynamics makes the problem truly unsteady. As an example, Fig. 3 shows ‘instantaneous’ pictures of velocity fields at the simulation time of five years with variations in both the air flow rate Q (upper 0.06, average 0.18 and lower 0.30 m³/s) and the effective heat conductivity coefficient k_w (1.0 (Fig. 3a) and 2.0 (Fig. 3b) W/(m×K)) in the heat release area. To display the velocity field, the option of ‘normalized’ arrows was used, which makes it possible to see the structure of the velocity field in the areas where the absolute values of the velocity components are close to zero and, when the option ‘proportional’ arrows is used, the vectors turn into points.

Figure 2.

The structure of the velocity field in the numerical experiments (2010) (ANSYS 2016).

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 (Amosov and Podshivalova 2010a, b). In addition, a visual comparison of ‘instantaneous’ pictures (see Fig. 3) makes it possible to qualitatively assess the effect of the incoming air flow rate.

Figure 3.

‘Instantaneous’ picture of the structure of the velocity field for five-year modeling with variations in the incoming air flow rate (upper 0.06, average 0.18 and lower 0.30 m³/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. 3a shows that at a flow rate of 0.06 m³/s (upper picture), the reverse air flow in the middle of the facility occupies almost half the height of the space between the surface of the reinforced concrete structure and the roof. As the air flow rate increases to 0.18 m³/s middle picture), the reverse flow area along the roof decreases by about a third. And with a further increase in the flow rate up to 0.30 m³/s (lower picture), the reverse flow at the considered moment of time practically disappears. It seems that in a situation of maximum air flow rate over time, and with a monotonic decrease in the volumetric power of residual heat release, the flow structure will more and more resemble the structure of the velocity field for the model with an incompressible fluid.

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³/s (middle pictures), Fig. 3a, b clearly show the difference in the flow structure in the upper left area of the facility (up to the built-in structure with SNF). At maximum air flow (lower picture), noticeable differences in the velocity field structure are clearly visible under the roof in the center of the facility and on the right, closer to the outflow.

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 2, 3)). The most detailed analysis was carried out for the option of an air flow rate of 0.18 m³/s and an effective heat conductivity coefficient of the heat release area of 1.0 W/(m.K). For almost a year and a half, the air movement was unidirectional towards the outflow boundary. However, over time, a reverse flow begins to form in different vertical sections – in a section of 90 m in about a year and a half, and in sections of 70 and 80 m negative horizontal velocities begin to appear after two years from the beginning of the simulation process.

Thus, it becomes obvious that, in the models of an incompressible liquid, on the basis of which the studies of 2010 were carried out (Amosov and Podshivalova 2010a, b) and, in the model of an ‘incompressible ideal gas’, the convective mechanism of heat transfer, which is dominant in the air, differs significantly. If, in the model of an incompressible fluid (Amosov and Podshivalova 2010a, b), there is practically unidirectional heat transfer to the outflow from the storage facility, in the model of an ‘incompressible ideal gas’, convective heat transfer varies both in time and in space. It is obvious that the result of these differences can be certain changes in the final indicators of the thermal regime of the facility.

Thermal regime of the storage facility

To answer the question about the conditions for meeting the criterion temperature values, the space-time temperature distributions for the above-mentioned 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 space-time temperature distributions shows that the search for maximum temperatures in different areas of the model must be performed in different sections (see above):

for the region of heat release and granite near the base, the cross section is 70 m;
for the surface of the reinforced concrete structure, the cross section is 80 m; and
for granite near the roof and air, the cross section is 90 m.

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³/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 chamber-type, air, built-in 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 (Amosov and Podshivalova 2010a, b), then, in a more complete model of an ‘incompressible ideal gas’, the forecast temperature values are different.

Analysis of space-time 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 non-stationary 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:

a decrease in the air flow rate leads to heating up (to various degrees) of all the regions of the model without exception; and
an increase in the effective heat conductivity coefficient of the heat release region decreases the predicted temperature of the rock mass near the base and the heat release region, but increases the temperature of the surface of the engineering structure, the air of the outgoing jet and the rock mass of the roof.

Table 3 shows the values of the maximum temperatures in various regions of the model with variations in the heat conductivity coefficient of the heat release area and the supplied air flow rate.

Table 3.

Download as

CSV

XLSX

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³/s	Controlled areas of the model
k_w, W/(m×K)	Q, m³/s	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 built-in 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. 4. The ‘minus’ sign for the heat flow directed through the base of the built-in structure is due to the selected direction of the vertical axis.

Figure 4.

Dynamics of the heat power leaving the heat release region through (a) the lower base and (b) the upper surface of the built-in structure: I – 0.06 m³/s, 1.0 W/(m×K); II – 0.06 m³/s, 2.0 W/(м×K); III – 0.18 m³/s, 1.0 W/(m×K); IV– 0.18 m³/s, 2.0 W/(m×K); V – 0.30 m³/s, 1.0 W/(m×K); VI – 0.30 m³/s, 2,0 W/(m×K).

Figures 4 clearly shows that the effect of the effective heat conductivity coefficient of the heat release area is more significant than that of the air flow rate. The graphs are well grouped in three, following the variable air flow rate, according to the value of the effective heat conductivity coefficient of the heat release area. At the same time, for air flow rates of 0.18 and 0.30 m³/s, the graphs of heat flows practically merge and, for a flow rate of 0.06 m³/s, the graph is slightly lower than the other two. It seems that, since the heat flows directed through the considered boundaries are interrelated, the probable reason for this behavior of the graphs can be precisely the structure of the velocity field above the surface of the built-in structure (see Fig. 3). Let us remember that it is for the minimum air flow rate of 0.06 m³/s that the structure of the velocity field differs significantly from the structure of velocity fields for other situations.

The qualitative and quantitative indicators of the graphs in Fig. 4 show that we are dealing with monotone functions. In terms of absolute values, at the initial stage of the heating process, the derivatives of the analyzed functions are significant but with time the rate of growth of heat flows gradually decreases. It is clearly seen that the dominating heat flow is directed into the air, where the convective transfer mechanism ensures heat transfer and mixing from the surface of the structure and its further removal from the SNF storage facility.

The behavior of the graphs in Fig. 4b shows that with an effective heat conductivity coefficient of the structural material of 1.0 W/(m×K), the heat flow reaches its maximum value (about 780 W) within five years and, with a doubled value of the effective thermal conductivity, the process of reaching a maximum value of about 850 W occurs earlier (3.0–3.5 years). These moments are in full agreement with the previously noted features of the spatial temperature distribution dynamics. In particular, this confirms the forecast decrease in the maximum temperatures of the rock mass with an increase in the effective heat conductivity coefficient of the built-in structure material.

The graphs in Fig. 4a illustrate a gradual decrease in absolute value of the heat flow, which is associated with gradual heating of the enclosing rock mass. For the maximum effective heat conductivity coefficient of the heat release area (2.0 W/(m×K)), the mass is heated up much faster, which is expressed by the outflow of the curves in Fig. 4a practically on the ‘plateau’ for five years. With the minimum thermal conductivity coefficient of the heat release region (1.0 W/(m×K)), the rock mass receives a greater amount of heat, which is especially typical for a situation with an air flow rate of 0.06 m³/s. As the graph in Fig. 4a shows, for the described situation the rock mass heating will continue after five years but with a slowing down rate.

Conclusion

Based on the conducted research, the following conclusions can be formulated.

A computer model of an underground facility for long-term spent nuclear fuel storage was constructed in the version of a built-in reinforced concrete structure in a two-dimensional 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 built-in 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 built-in 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.

References

Amosov PV, Kozyrev SA, Nazarchuk OV (2018) Development of a Computer Model of Aerothermodynamics of the Pit Open Atmosphere in ANSYS FLUENT. Bulletin of the St. Petersburg State Technological Institute (Technical University) 44(70): 121–125. https://doi.org/10.15217/issn1998984-9.2018.44.121 [in Russian]

Amosov PV, Novozhilova NV (2014) Two-Dimensional Numerical Simulation of Atmosphere Aerothermodynamic of Open Pit. Proc. of the All-Russian Sci.-Tech. Conf. Environmental Strategy for the Development of the Mining Industry – the Formation of a New Worldview in the Development of Natural Resources. Mining Institute of the Kola Scientific Center RAS: 153–159. [in Russian]

Amosov PV, Podshivalova AV (2010a) Modeling of the Thermal Regime of an Underground Storage Facility for Heat-Emitting Radiation-Hazardous Materials. Vestnik MGTU [MSTU Bulletin] 13(3): 562–566. [in Russian]

Amosov PV, Podshivalova AV (2010b) Features of the Thermal Regime of an Underground Object of Insulation of Heat-Generating Materials (Container Version). Izvestia vysshikh uchebnykh zawedeniy. Yadernaya energetika [News of Higher Educational Institutions. Nuclear Power Enguneering] 3: 133–140. [in Russian]

ANSYS (2016) Release 17.0 Documentation. Swanson Analysis Systems.

Chirkin VS (1968) Thermophysical Properties of Nuclear Engineering Materials. Moscow. Atomizdat Publ., 485 pp. [in Russian]

COMSOL Documentatio (2020) COMSOL Documentation. https://www.comsol.ru/documentation [accessed Jan 20, 2020] [in Russian]

Dortman NB (1984) Physical Properties of Rocks and Minerals (Petrophysics): Handbook of Geophysics. Moscow. Nedra Publ., 455 pp. [in Russian]

Earth Science (1969) Vol. 21: Handbook of Physical Constants of Rocks. Moscow. Mir Publ., 543 pp. [in Russian]

Horai K (1971) Thermal Conductivity of Rock-Forming Minerals. Journal of Geophysical Research 76(5): 1278–1308. https://doi.org/10.1029/JB076i005p01278

Kalinkin VI, Kritsky VG, Tokarenko AI, Tikhonov NS, Razmashkin NV, Serova AL, Balitskaya AN (2009) Storage of Spent Nuclear Fuel of Power Reactors: Preprint. St. Petersburg: VNIPIET Publ., 124 pp. [in Russian]

Melnikov NN, Konukhin VP, Naumov VA, Amosov PV, Gusak SA, Naumov AV, Katkov YR (2003) Spent Nuclear Fuel of Ship Power Plants in the European North of Russia. Part II. Apatity. KNTS RAN Publ., 209 pp. [in Russian]

Melnikov NN, Konukhin VP, Naumov VA, Amosov PV, Gusak SA, Naumov AV, Katkov YR, Smirnov YuG, Orlov AO, Rybin YuYu (2006) Concept of Underground Storage of Spent Nuclear Fuel for Marine Nuclear Power Plants on the Kola Peninsula. Vestnik MGTU [MSTU Bulletin] 9(3): 408–418. [in Russian]

Melnikov NN, Konukhin VP, Naumov VA, Amosov PV, Gusak SA, Naumov AV, Orlov AO, Smirnov YuG, Karavaeva EV, Novozhilova NV, Klimin SG (2010) Scientific and Engineering Aspects of Safe Storage and Disposal of Radiation-Hazardous Materials in the European North of Russia. Apatity. KNTS RAN Publ., 305 pp. [in Russian]

Melnikov NV, Rzhevsky VV, Protodyakonov MM (1975) Handbook (Cadastre) of Physical Properties of Rocks. Moscow. Nedra Publ., 279 pp. [in Russian]

Naumov VA, Gusak SA (2019) A Study on Formation Regularities of Heat Sources in the Cask Type Storage Facilities for Spent Nuclear Fuel of the Low-Power Rectors. Vestnik KSC RAS [KSC RAS Bulletin] 2(11): 105–115. https://doi.org/10.25702/KSC.2307-5228.2019.11.2.105-115 [in Russian]

NP Safety Requirements-035-02 (2002) Dry Storage Facilities for Spent Nuclear Fuel. Vestnik Gosatomnadzora Rossii [Gosatomnadzor Bulletin] 3: 51–58. [in Russian]

Nuclear Energy Facilities NP-061-05 (2005) Safety Rules for Storage and Transportation of Nuclear Fuel. https://files.stroyinf.ru/Data1/47/47340/ [accessed Jan 20, 2020] [in Russian]

Pekhovich AI, Zhidkikh VM (1968) Calculations of the Thermal Regime of Solids. Leningrad. Energia Publ., 304 pp. [in Russian]

Rzhevsky VV, Novik GYa (1978) Fundamentals of Rock Physics. Moscow. Nedra Publ., 390 pp. [in Russian]

Shchelyaev A FlowVision-modern Russian mathematical modeling tool. https://sapr.ru/article/21879 [accessed Aug 21, 2018] [in Russian]

Vargaftik NB (1972) Handbook of Thermophysical Properties of Gases and Liquids. Moscow. Nauka Publ., 720 pp. [in Russian]

Yegorov VI (2006) Application of Computers for Solving Problems of Thermal Conductivity. Tutotial. St. Petersburg: SP GU ITMO Publ., 77 pp. [in Russian]

* Russian text published: Izvestiya vuzov. Yadernaya Energetika (ISSN 0204-3327), 2020, n. 3, pp. 80–92.