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

Academic editor: Yury Kazansky

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.

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 (

Forecast calculations of the thermal state of materials of an underground

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 chamber-type for storing spent nuclear fuel in a built-in reinforced concrete structure (

Model geometry of the underground spent nuclear fuel storage facility (dimensions are in meters) (

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

where _{p}_{w}

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),

where η is the dynamic viscosity; _{0}. The result is a buoyancy force expressed as (ρ – ρ_{0})

ρ = _{0}×μ/(

lg

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

where _{0} = 101325 Pa; μ = 0.0288 kg/mol;

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),

The initial temperatures of the areas of the model and the values of the thermophysical parameters adopted in the calculations (

Initial temperatures _{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 (_{w}

Technically, the values of the thermophysical parameters of complex heterogeneous systems can be calculated using the relations given in (^{3}/s were used.

The volumetric power curve of the residual heat release _{w}

_{w}^{2} – 1.651×10^{–5}×^{3},

where

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 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),

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 (^{3}/s) and the effective heat conductivity coefficient _{w}

The structure of the velocity field in the numerical experiments (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 five-year 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. ^{3}/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^{3}/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^{3}/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^{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 ^{3}/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 (

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^{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 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

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 (_{w}

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

Calculated values of the maximum temperatures in different regions of the model with variations in _{w}

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 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.

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^{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 ^{3}/s, the graphs of heat flows practically merge and, for a flow rate of 0.06 m^{3}/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}/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.

The behavior of the graphs in Fig.

The graphs in Fig. ^{3}/s. As the graph in Fig.

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.

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