Approaches to radiative heat transfer simulation in a cavity above melt

. The article presents a description of the THERA module, which allows the calculation of radiative heat transfer in the cavity above the melt in the complex geometry with the presence of participating medium in relation to the simulation of severe accidents at nuclear power plants with VVER type reactor. The net-radiation method is implemented in the module as a basic method for calculating radiative heat transfer; the radiation properties of gases are considered using the SLWSGG model. The results of THERA cross-verification on test tasks are presented. The paper presents the results of THERA simulation of radiative heat transfer from the corium surface to the reactor shaft structures at a nuclear power plant with a VVER – 1000 reactor as part of the analysis of the ex-vessel stage of a severe accident in axisymmetric and realistic non-axisymmetric forms.


Introduction
The safety of nuclear power plants with a VVER reactor is based on the concept of defence in depth (DiD), assuming the presence of barriers on the path of radioactive fission products to the environment.Retention of radioactive substances within safety barriers and their integrity can be distinguished as the main functions of the DiD concept.The reactor core or fuel in the pool are heated and melted, followed by the interaction of the melt with the NPP equipment, during a severe accident.Localization and cooling of the formed corium are the main tasks to mitigate the consequences of a severe accident to prevent subsequent failures of safety barriers.
During a hypothetical severe accident, in case of failure of the reactor vessel and the absence of the core catcher, the formed corium goes outside the vessel into the concrete reactor shaft.The melt is a viscous liquid, heated above 2000 K, with volumetric residual heat; it consists of molten core components and structural materials of the reactor vessel.Radiative heat transfer is the main mechanism for removing heat from the corium at such high temperatures and without any additional heat removing methods.Thermal radiation from the corium surface heats up the walls of the reactor shaft located above the surface and the elements of the reactor vessel remaining after the failure.Intense thermal radiation on the reactor shaft structures leads to their melting and moving into the melt.Specialized computational models are required to analyze the effectiveness of measures to mitigate consequences of thermal radiation affecting the shaft construction, as well as to assess the survivability of equipment in such conditions, to determine the * Corresponding author: shm9688@gmail.comdynamics of their heating and melting.These models consider the features of this task: -complex geometry of reactor plant structures and containment; -changes in the geometry of structures during heating and collapsing; -substantial inhomogeneous temperature distribution over the surfaces of the computational domain; -the presence of gases in the cavity above the melt, which are characterized by absorption and emission of radiation.
During the simulation of a severe accident, thermal radiation from the melt surface is often considered simplistically using the Stefan-Boltzmann lawas the radiative heat flux from the melt surface to a fictive surface with desired temperature.This approach does not allow us to estimate the thermal loads on certain structures above the melt level and correctly take into account the change in the geometry over time.Also, it is necessary to consider the presence of gases in the cavity above the corium during the radiative heat transfer calculation in the analysis of severe accidents at nuclear power plants (NPP).At the ex-vessel stage of a severe accident occurs the thermal decomposition of concrete under the influence of a high-temperature melt, what is accompanied, according to [1], by an intense release of steam and carbon dioxide.The oxidation reactions between released gases and the metal components of the melt led to formation of hydrogen and carbon monoxide.Formed gases cause absorption and emission of thermal radiation in the cavity above the melt, which requires correct consideration during calculation of radiative heat transfer.The estimates given in [2] show that even a small mole fraction of steam (0.05 -0.2 mole fraction) E3S Web of Conferences 459, 07011 (2023) https://doi.org/10.1051/e3sconf/202345907011XXXIX Siberian Thermophysical Seminar in the steam-gas mixture led to absorption and emission of radiation in the gas volume.
This article describes approaches for radiative heat transfer simulation in the cavity above the melt with the complex geometry, implemented in the THERA (Thermal Radiation) module of the TSAR (Toolkits for Severe Accident Research) application package.The results of the cross-verification of the THERA module on test tasks in a transparent cavity and in a cavity filled with steam-gas mixture are presented in this article.During operation of the THERA module is calculated the radiative heat transfer from the melt surface to the reactor shaft structures for the analysis of the ex-vessel stage of a hypothetical severe accident at a nuclear power plant with a VVER -1000 reactor.Several cases of melt configuration were simulated in relation to two domains: axisymmetric and realistic, the latter taking into account the lack of axial symmetry of the shaft construction.

THERA module
The THERA module is designed for: the calculation of the radiative heat fluxes at the late stage of a severe accident at a nuclear power plant with a VVER reactor in order to determine the dynamics of steel structures melting above the melt level at the in-vessel stage of the accident; the analysis of thermal loads on the equipment of the reactor structures or on the core catcher at the ex-vessel stage of the accident; the evaluation of the effectiveness of measures of cooling water supply to the corium surface to shield the equipment from thermal radiation from the surface of the melt.
The THERA module implements a threedimensional form of the task in order to consider the complex geometry of computational domain.The computational domain is a set of surfaces with specified properties and boundary conditions for which is performed the calculation of heat fluxes and heat sources, associated with the absorption and emission of radiation in the cavity.The discretization of the computational domain is performed in a mutually consistent manner with TSAR code using a non-structural grid of triangular elements with the option of mesh refinement.

THERA description
In the THERA module, the calculation of radiative heat transfer is based on the net-radiation method, first introduced by Hottel [3].The three-dimensional computational domain simulating the cavity above the melt is divided into triangular finite surfaces.The following assumptions are implemented within each surface: -Temperature of the surface is constant; -Surface is considered grey, absorptivity is equal to emissivity and determined only by the surface temperature T: ε(T) = α(T) = 1p.This assumption greatly simplifies the analysis due to the fact that radiation properties of real surfaces and their spectral characteristics are difficult to determine because of many affecting factorsfrom the degree of contamination and surface roughness to the composition of the material, which varies widely in a multicomponent corium; -radiation emitted or reflected by a surface is considered diffuse; -incident and reflected radiation fluxes are considered constant.
The steam-gas mixture in the cavity is considered isothermal, homogeneous, radiation absorbing and emitting are considered, the effects of scattering are not taken into account.
Figure 1 shows a schematic representation of a two-dimensional cavity, which illustrates the difference in the calculation of closed and unclosed computational domains.Fictive surface is added for an open cavity, thus closing the cavity.The heat, leaving the domain in this case, is associated with fictive surface.

Radiation properties of gases
The major difficulty in calculation of radiative heat transfer in a cavity above the melt with steam-gas mixture is to determine the radiation properties of gases due to the discrete nature of the spectral characteristics and the irregular dependence of the absorption coefficient on the wavelength of the spectrum.
In THERA module, the spectral-line-weighted sum of gray gases (SLWSGG) model was chosen for calculation of the radiation properties of gases [4].This model is an improved and more precise version of the classic weighted sum of gray gases (WSGG) model [3] and is a representative of a class of models based on the global absorption-line blackbody distribution function (ALBDF).The main feature of the global models is the transition from integration of the spectral radiative properties of gases over the wavelengths of the gas absorption spectrum to summation over gray gases simulating the absorption spectrum.The transition from integration over absorption bands to summation over a set of gray gases, the absorption coefficients of which are constant values, greatly simplifies the analysis and provides high calculation speed.
The spectrum of the gas absorbtion is divided into M gray gases, for each of which the radiation transfer equation (RTE) in a non-scattering medium is written in the form [5]: Here   ,   ,  radiation intensity, absorption coefficient and weight coefficient for j grey gas,  the blackbody radiation intensity,   =

𝜋𝜋
, the beam length,  the gas temperature.
The absorption coefficient for j grey gas can be written as: , molar concentration and mole fraction of the gas,  ̃−1 ,  ̃absorption cross-section for j-1 and j grey gases.
The values of the absorption cross sections for any of j grey gases are logarithmically distributed between the minimum   and maximum   values of the absorption cross-sections, which represent the entire gas absorption spectrum: The calculation of the weights coefficients is based on the concept of the absorption distribution function which is defined as a fraction of the energy of blackbody, emitting at a temperature   , in the absorption spectrum of a gas with mole fraction   , with pressure P and temperature   for which the absorption cross-section   is less than a chosen C value: The weight coefficient for the j-th gas is defined as the difference between the values of the ALBDF calculated on neighboring regions of the absorption cross-section: In the case of a steam-gas mixture consisting of N components in the cavity above the melt, the multiplication approach from [6] is used to consider the radiation properties of the mixture.This method is based on multiplying the ALBDF of individual components of the mixture: The calculation of the ALBDF is performed according to the correlations given in [4].Parameters required for the calculation of correlations are taken from [7].They were obtained based on the HITEMP-2010 high-resolution spectroscopy database for H 2 O, CO 2 и CO.

Computational method
A closed cavity consisting of N triangular surfaces is filled with absorbing and emitting gas with temperature   .Expression for resulting radiation heat flux density of the k surface for j grey gas can be written: ,,density of outgoing radiation heat flux from the surface k or radiosity,  ,,density of incident radiation heat flux into k surface or irradiation,  ,, the density of the resulting radiation flux.The density of radiative flux emitted by the surface k can be written [5] using the grey gas coefficient   of SLWSGG model.Than the expression of radiosity for surface k and j grey gas is: ε emissivity of the k surface, the Stefan-Boltzmann constant,  temperature of the k surface,  ,,density of incident radiation heat flux on the surface k for j-th grey gas.Solving equation (1) analytically and obtaining the integral-differential form of the RTE for the j-th gas, we obtain an expression for determining the density of the incident radiation flux for the j-th gray gas on the surface k from all N surfaces, considering the absorption and emission of radiation in the gas volume with beam length S:  ,, = ∑  ,,  ,  , +  ,  ,     4  =1 (9)  ,the view factor between surfaces k and i.In equation ( 9)  , ,  ,transitivity and absorptivity of jth grey gas: the distance between surfaces k and j,  absorption coefficient of jth grey gas, calculating by SLWSGG model.
The resulting system of linear algebraic equations with respect to the density of effective radiation heat fluxes for j-th grey gas is obtained using ( 8) and ( 9).The values of heat fluxes densities   supplied from the outside are set as boundary conditions for some surfaces, temperatures T are set as boundary conditions for other surfaces.The final system of equations has the form: δ ki -the Kronecker delta, the first equation of the system is written for surfaces with temperature boundary conditions   , the second equationfor heat flux boundary condition with value   .
As a result of solving (12) is obtained a set of densities of effective radiation heat fluxes for N surfaces, that provides the relevant densities of the resulting radiation heat fluxes for the j-th gray gas.According to [5], the values of the integral densities of the resulting radiation fluxes are obtained by summing over M gray gases: In the absence of an absorbing and radiating medium in the cavity: the absorption coefficient of the medium is zero, the weight coefficients are assumed to be equal to one and the system of equations (12) describes the radiative heat transfer in an optically transparent cavity.

View factors calculation
The net-radiation method is based on the calculation of the view factors of surfaces in the computational domain.The view factor is a geometric characteristic of the computational domain, which determines the proportion of radiation emitted by surface 1 and incident on surface 2.
In THERA, the following expression is used to calculate view factor between two elementary surfaces i and j: where   ,  surface areas of i and j,  vector connecting the barycenters of surfaces i and j,     are angles between the normals to surfaces i and j and the vector  ⃗  .
Figure 2 shows a schematic representation of the values involved in calculating of the view factor between two surfaces involved in heat transfer.The calculation of view factors in a complex geometry implies correct consideration of possible shadow effect between surfaces.The Moller-Trumbore algorithm [8] is used in THERA to consider shadow effect of the surfaces.This algorithm searches for a possible intersection point of the beam and the shading surface on the path of the beam to the target surface.As can be seen from Figure 3, if Surface 2 occurs on the path of the beam emitted from Surface 1 to Surface 3, then it is assumed that surfaces 1 and 3 clearly do not participate in radiative heat transfer, and the corresponding view factor is equal to zero.The Moller-Trumbore algorithm is characterized by its high computational efficiency, as well as the possibility of additional acceleration due to parallelization of GPU (CUDA) calculation.The Moller-Trumbore algorithm is based on searching of the intersection point of the beam emitted from the barycentre of surface 1 and the surface 2, the vertices of the surface 2 are set by coordinates V 1 , V 2 and V 3 .The equation of the beam line is given in parametric form, the coordinates of the surface 2 and the intersection points are given in barycentric coordinates.
The back-face culling method [9] is used in the THERA module to improve the computational time of the shadow effect algorithm.This method determines the position of the surface, receiving the emitted beam from the observer surface, relative to the observer surface itself.The position of a surface is determined by comparing the triangle vertices winding order of the two surfaces.The receiving surface is considered culled relative to the observer surface in the case when they have different triangle vertices winding order.

THERA validation
Experiments with radiative heat transfer as a predominant phenomenon are mainly limited to the problems of furnaces and spacecraft.Due to the lack of experimental data for THERA validation, the module was cross-verified on test tasks.Verification of THERA on the problems of calculating radiative heat transfer in an optically transparent medium was performed for a task with a standard geometry characteristic of severe accident simulation objects, the comparison was made with the calculation of CFD.A comparison of the THERA calculation with the work of Liu [10] was performed in order to validate THERA module on radiative heat transfer simulation in gas cavity.

Optical transparent medium
In the case of radiative heat transfer calculation in the optical transparent medium is considered cylindrical geometry, a similar geometry is characteristic of the reactor vessel, the core catcher and the concrete shaft of the reactor.The radius of the cylinder is equal to 2 m, height of the cylinder is 3 m, emissivity of all surfaces in the domain set to ε = 0.4.Temperature boundary conditions  3 = 3000 К,  2 = 1500 К are set on the bases of cylinder 3 and side surface 2, respectively, adiabatic boundary conditions are set on the base of cylinder 1.
The temperature distribution along the diameter of the cylinder 1 base was chosen as a parameter for comparison with CFD, which is shown in Figure 4. Radiative heat transfer was simulated by THERA for meshes with different densities: 2794, 3836, 5390 elements, meshes of border areas were refined for all cases.The differences in the results between THERA and CFD in the areas close to the side surface of the cylinder are explained by interpolating the values of temperature inside the boundary cells in CFD simulation.The error dependence on the mesh density of the computational domain of the calculation results by the THERA module was compared to the CFD code results.The mesh partition in the central part of the cylinder base makes the main contribution to the error.Figure 5 shows that with an increase in the number of elements of the computational domain, there is a better agreement of the THERA module with CFD.The error does not exceed 1% for all cases.

Gaseous medium
A three-dimensional rectangular computational domain was chosen for the test task, the square in the base of domain has its side equal to 2 meters, and the height of the rectangular area is 4 meters.The rectangular cavity is filled with steam at a temperature   = 1000 К, temperature boundary conditions of the walls are set as   = 300 К, the steam pressure is 0.1 MPa.The results of simulation by THERA were compared with the Liu's results [10], in which the calculation of radiative heat transfer with similar conditions was performed by discrete ordinates method, the statistical narrow-band model was used to obtain the radiation properties of steam.Statistical narrow-band model is more accurate than global models [11], thus Liu's work can be used as a benchmark.
The mesh consisting of 2000 surfaces was chosen to match the calculation conditions of Liu's work.Figures 5 and 6 show the calculated values of the densities of radiative heat fluxes on the side and upper surfaces, respectively, in comparison with the Liu's work.As can be seen, the calculation by THERA with calculating radiation properties of gas by SLWSGG demonstrates good agreement with Liu's work, slight differences in the results may be caused by differences in approaches to the mesh partition of the computational domainin Liu's work was used a structured rectangular mesh, while in THERA an unstructured triangular one.
An analysis of the influence of the grey gases amount on the error of the simulation result relative to the Liu's values was performed.It was found that the increase in number of grey gases provides the better agreement of THERA simulation with the benchmark, and it can be noted that dividing the spectrum on 10 -15 gray gases gives the best result in terms of error with respect to Liu's work and computational efficiency.

Radiation heat transfer in the reactor shaft 4.1 Features of ex-vessel stage of a severe accident
The simulation of radiative heat transfer is performed as part of the analysis of the ex-vessel stage of a severe accident at a nuclear power plant with a VVER -1000 reactor.As a representative scenario for the calculation of a hypothetical severe accident was taken the generalized scenario of the guillotine break of main circulation pump DN 850 on the inlet of the reactor with simultaneous station blackout.This scenario is characterized by the earliest release of the melt from the reactor vessel (about one and a half hours after the accident onset), which causes the greatest residual heat release in the melt, in comparison with other scenarios.The early movement of the melt into the reactor shaft corresponds to a lower degree of oxidation of zirconium, which increases the intensity of chemical oxidation reactions in the melt during molten corium -concrete interaction and leads to an increase in heat release.The initial data on the corium composition entering the reactor shaft are obtained as a result of simulation with SOCRAT B1/B2 code of the in-vessel stage of a severe accident.Data are shown in Table 1.During the in-vessel stage, the melt in the reactor vessel is stratified into a metal layer containing molten steel and zirconium, and an oxide phase composed of uranium and zirconium dioxides.The melting of the reactor vessel occurs at the level of the metal layer, which determines the specifics of the corium entering the shaft: first of all, the metal phase moves, then, as the vessel melts, the oxide components enter.A molten pool is formed on the bottom of the reactor shaft, the corium begins to interact with the concrete structures of the shaft.Simulation of the molten corium -concrete interaction in the reactor shaft is performed with the HEFEST-ULR code in an axisymmetric geometry.After the significant mass of oxide components enters the non-stratified melt, direct stratification of the melt occurs -the heavier oxide components sink down, the light metal layer floats up.As a result of concrete ablation, the melt is saturated with light oxides of the concrete components, which leads to decrease in the density of the oxide layer and inversion of the meltthe oxide phase floats, the metal phase moves down.These features of the melt behavior affect the radiative heat transfer from the surface of the corium to the shaft structures due to the dependence of emissivity on the composition of the melt.Figure 7 shows the change in the temperature of the melt surface during the ex-vessel stage of a severe accident.The surface temperature of corium during interaction with concrete does not fall below 2000 K, radiative heat transfer is the main mechanism of heat removal from the melt.
The molten corium -concrete interaction is accompanied by the release of decomposition gasessteam and carbon dioxide, as well as the formation of combustible gases because of oxidation reactions of the metal components of the melthydrogen and carbon monoxide.According to [12], gases with symmetric molecules H 2 , O 2 , N 2 hardly absorb or emit radiation in the gas volume, thus they can be omitted during simulation.

Simulation conditions
The reactor shaft of the VVER -1000 is a cylindrical concrete structure with incomplete axial symmetry; it belongs to the accident localization zone.The walls of the shaft are made of silicate concrete, the shaft structure is characterized by the presence of a rectangular compartment with two protective doors: the first is made of concrete with lining of carbon steel, the second is entirely made of carbon steel.
The calculation with THERA module of radiative heat transfer is performed for three corium's configurations to estimate the influence of corium surface emissivity on simulation results: -Case №1: non-stratified melt, consists of metal components, the composition of the vapour-gas mixture in mole fractions: 0.12 air, 0.48 H 2 O, 0.02 CO 2 , 0.36 H 2 , 0.02 CO, the pressure of gas mixture is 0.31 MPa, the mixture temperature is 413 K; -Case №2: stratified melt with upper metal layer, the composition of the vapour-gas mixture in mole fractions: 0.10 air, 0.52 H 2 O, 0.02 CO 2 , 0.32 H 2 , 0.026 CO, the pressure of gas mixture is 0.30 MPa, the mixture temperature is 530 K; -Case №3: inversion corium, with upper oxide layer, the composition of the vapour-gas mixture in mole fractions: 0.04 air, 0.40 H 2 O, 0.025 CO 2 , 0.50 H 2 , 0.027 CO.The pressure of gas mixture is 0.27 MPa, the mixture temperature is 515 K.
The results of simulation of the ex-vessel stage of a severe accident with the HEFEST-ULR code were obtained within the axisymmetric approximation of the computational domain, without considering a rectangular compartment.As a result, at the first stage of the work, the calculation of radiative heat transfer was performed in an axisymmetric geometry, the scheme of which is shown in Figure 8 with an indication of the surfaces involved in heat transfer.The boundary conditions of variant calculations are presented in Table 2.The highest density of the radiation heat flux on the wall is observed in the case 2, which is associated with the highest surface temperature of the melt -3192 K.It should be noted, that at similar temperatures of the melt surfaces for the case 1 and case 3, the highest heat flux density on the wall is observed in the third scenario, which is due to the high value of the emissivity of the main component of the oxide layer (uranium dioxide) in comparison with the emissivity of the metal layer.With an increase in the height of the reactor shaft wall, the values of the radiative heat flux densities become similar regardless of the melt configuration and surface temperature.This can be explained by the effects of radiation absorption in a sufficiently dense steam-gas mixture for all three cases of calculation.According to the calculation results, the greatest heat release in the steam-gas mixture equal to 155 kW/m 3 was obtained in the case 2.

Non-Axisymmetric computational domain
Figure 12 shows a three-dimensional non-axisymmetric computational domain on which the simulation of radiation heat transfer by the THERA module was carried out.The mesh consists of 9000 elements, and refinement of the mesh was carried out at the base of the domain.Figures 13 and 14 show the distribution of the radiation heat flux density along the wall of the reactor shaft (surface 2) and along the door (surface 4), obtained during the simulation of radiative heat transfer by the THERA module for the reactor shaft in the nonaxisymmetric formulation.As in the simulation in the axisymmetric formulation, the greatest thermal load due to radiation heat transfer is observed in the case 2 at the maximum temperature of the melt surface.The values of the heat flux densities on the shaft's wall are less than in the case of an axisymmetric formulation due to the redistribution of radiation heat fluxes inside the computational domain.As in the simulation in the axisymmetric formulation, there is a strong influence of the melt configuration on the distribution of heat flux densities across the shaft's structures due to the significant difference in the emissivity of the melt surface for the metal and oxide layers.According to the calculation results, the greatest heat release in the steamgas mixture equal to 142 kW/m 3 was obtained in the case 2.  Figure 15 shows the distribution of heat flux density over one of the surfaces in a "cold" space between the safety doors.The occurrence of heat flux density on the walls of the room isolated from the melt is explained by the presence of a vapor-gas mixture due to the fact that the protective door is not sealed.The higher the temperature of the steam-gas mixture, the higher the density of thermal radiation fluxes on the wallthe temperature in the case 1 is less than in the cases 2 and 3 more than 100 K.

Conclusions
High accuracy of radiative heat transfer simulation at the late phase of a severe accident after molten pool formation is necessary for: realistic estimation of the development of the emergency process, the development and evaluation of the effectiveness of measures aimed at preserving the last safety barrier.To improve the accuracy of simulation and the realism of calculated safety analyses was developed the THERA module of the TSAR application software package in the C++ programming code, designed to calculate radiative heat transfer in the cavity above the melt, which allows the calculation of radiative heat fluxes in case of tasks with complex geometry in transparent media and media filled with steam-gas mixture.The THERA module is based on the net-radiation method adapted for a task with a three-dimensional form and the presence of shading surfaces.The calculation is performed on a nonstructural mesh of triangular elements.The cross-verification of the THERA module was performed on test tasks in an optically transparent medium and in the presence of steam in the cavity to check the operability of the THERA module.For the case of the transparent medium was performed the comparison with a CFD code.Liu's work was used as a benchmark for calculation of radiative heat transfer in the presence of a steam-gas mixture in the cavity.The results of cross-verification demonstrated good convergence of the THERA module with CFD and benchmark.
Simulation of radiative heat transfer from the corium surface on the reactor shaft structures was calculated with THERA module in the relation to the analysis of a hypothetical ex-vessel stage of the severe accident at nuclear power plant with a VVER-1000 reactor.This article presents the results of radiative heat transfer simulation in a simplified axisymmetric form and in a realistic non-asymmetrical geometry of the reactor shaft constructions.The distribution of the radiative heat flux density on the protective door of a rectangular compartment of the reactor shaft room is presented, which is important from the point of view of melting and further spreading of the melt.In both cases of domains, simulation was performed for three configurations of the melt surface in the conditions of the ex-vessel stage of the accident in order to analyse the effect of the surface emissivity on the thermal load of the walls of the reactor shaft.The inverse corium configuration leads to an increase in the heat flux density on the shaft's walls, which is explained by the high value of emissivity of the main component of the oxide layer of the melt (uranium dioxide).

Fig. 4 .
Fig. 4. Distribution of the temperature along the base diameter.

Fig. 5 .
Fig. 5. Distribution of heat flux densities along the side of domain.

Fig. 6 .
Fig. 6.Distribution of heat flux densities along the upper base of domain.

Fig. 15 .
Fig. 15.Distribution of heat flux densities along the wall in the cold space between doors.

Table 1 .
Parameters of corium entering the shaft.