79urn:lsid:arphahub.com:pub:D015427E-53F6-553C-9279-2E51CD684756Nuclear Energy and TechnologyNUCET2452-3038National Research Nuclear University MEPhI (Moscow Engineering Physics Institute)10.3897/nucet.5.3931839318Research ArticleHeat transfer and Fluid dynamicsDevelopment of the model to determine the fuel temperature field in a two-dimensional problem statementGorbunovVladimir A.1gorbynov.w@mail.ruIvanovaNatalya B.1LonshakovNikita A.1BelovYaroslav V.1Ivanovo State Power University n.a. V.I. Lenin, 34 Rabfakovskaya Str., Ivanovo, 153003 Russian FederationIvanovo State Power UniversityIvanovoRussia
Corresponding author: Vladimir A. Gorbunov (gorbynov.w@mail.ru)
Academic editor: Boris Balakin
2019250920195321922456B51DFB-8C7B-5BB6-BD8F-8C92E0331BB334704931006201920082019Vladimir A. Gorbunov, Natalya B. Ivanova, Nikita A. Lonshakov, Yaroslav V. BelovThis 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.
Water-cooled water-moderated reactors (VVER) are widely used at Russian nuclear power plants. The VVER reactor core is formed by fuel assemblies consisting of fuel rods. The fuel in fuel rods is uranium dioxide. The safety of the reactor operation is ensured through stringent requirements for the maximum nuclear fuel temperature. Calculation of temperature fields within the reactor core requires associated problems to be solved to determine the internal energy release in fuel based on neutronic characteristics. Dedicated software for such calculations is not available to a broad range of users. At the present time, there are numerical thermophysical modeling packages available for training or noncommercial applications which are used extensively, including Elcut, Flow Vision, Ansys Fluent, and Comsol Multiphysics. Verification of the obtained results is becoming an important issue in building models using these calculation packages.
An analytical solution was obtained as part of the study for the fuel temperature field determination. A program was developed in MathCAD based on this solution. A model was developed in Comsol Multiphysics to determine the fuel temperature field with constant thermophysical properties in a two-dimensional problem statement. The numerical model was verified using the analytical solution. The influence of the number of the grid nodes on the solution accuracy was established. The analytical solution can be used to determine the fuel temperature field at any radial coordinate of the reactor. The temperature field determination model developed in MathCAD can be used to verify numerical models of the fuel temperature field determination developed in dedicated packages.
Fuel elementfuel temperature fieldanalytical solutionmathematical model to determine fuel temperature fieldsverification of numerical calculationsinfluence of tuning coefficientssafety of fuel element heatingIntroduction
Operating safety of a nuclear reactor requires stringent limits to be placed on the maximum temperature of fuel (Perimov et al. 2004). The yield of uranium fission gaseous products rises sharply at a temperature of over 1600 °C. This leads to the requirement for an operating safety procedure that sets the maximum permissible temperature in the middle of the fuel column at 1690 °C (Kolpakov and Selivanikova 2009).
The purpose of the nuclear reactor core is to generate thermal energy and transfer it from the fuel surface to the coolant. VVER-1000 reactors use cylindrically shaped bare-rod fuel elements. It has been found that the reactor volume power density, in the event of a homogeneous reactor, is proportional to the neutron flux and can be determined by the following expression (Dementyev 1990):
where Q_{0} is the maximum value of the power density at the reactor center, MW; J_{0} is the zero-order Bessel function; R_{e} is the effective radius, m; H_{e} is the effective height, m; and r, z are the coordinates of the calculation point radially and axially respectively, m (Fig. 1).
Power density distribution radially Q_{r} and axially Q_{z} in a homogeneous cylindrical core.
https://binary.pensoft.net/fig/341993
The fuel element is a pressurized tube of a zirconium-niobium mixture which is filled with uranium dioxide pellets (Dolgov 2016). Thermal energy is released within uranium dioxide. Since uranium dioxide has rather low conductivity, the temperature field in the column has a highly non-uniform distribution (Leskin et al. 2011).
To determine the temperature field of a fuel element, one can use dedicated commonly available numerical software products which make it possible to build a computer model of heat fluxes inside of the reactor core. After the model is built and calculated, it requires verification based on experimental data which is rather hard to find. An important task to be solved in the process of computer simulation is to determine the influence of adjustment coefficients (number of grid nodes, number of time steps, etc.) on the accuracy of the model operation results.
Therefore, one of the paper’s objectives is to obtain an accurate analytical solution for the fuel temperature field determination (Goltsev et al. 2009). The accurate analytical solution shall be used to verify the numerical model of the fuel temperature field determination.
Investigation techniques
Numerical simulation methods are used for the fuel temperature field investigation. An axially symmetrical model has been built and solved based on a conductivity equation in a 2D problem statement with an internal heat source with the second-order boundary conditions on the fuel element’s end faces and the third-order boundary conditions on the rod’s side surface (Zhukov et al. 2009). Software not available for a broad range of users (BERKUT (Batenin et al. 2017), SOKRAT, HYDRA, YEVKLID, KORSAR, RELAP5, MELCOR, ASTEC (New Generation Codes 2018, Vishnyakova et al. 2015)) are used as tools for such calculations. Commonly available software products, such as STAR-CCM+, STAR-CD, ANSYS CFX, and ANSYS Fluent, are used for CFD calculations at OKB Gidropress (Velesyuk and Morgunov 1990, Agranat et al. 2015, Kartashov and Bogoslovskaya 2012); apart from these codes, Russian-made software products (Elcud, Flow Vision) and the Comsol multi-physics package are used at universities.
Problem statement and solution
The purpose of the investigation is to determine the temperature field of an axially symmetrical fuel element model in a 2D problem statement with constant thermal-physical coefficients with regard for the internal heat source (Ushakov and Subbotin 1972, Loginov 2009). A Comsol Multiphisics software package is used for the problem solution by finite-difference element method. The obtained results of the fuel temperature fuel calculation need to be further verified using the accurate analytical solution of the heat conductivity equation for the fuel element (Starkov and Marikhin 2013). Third-order boundary conditions on the fuel element’s side surface are used to solve the problem, and the conditions on the end faces are assumed to be adiabatic.
The following assumptions are used for the problem solution:
the thermal-physical properties of uranium dioxide are constant;
the uranium-dioxide rod does not have a zirconium alloy cladding;
when the fuel element is heated, no effects of it being flown over by the fluid current are taken into account.
We shall consider the fuel conductivity equation
δT/δt = λΔT/(cρ) + q/(cρ), (2)
where T is the temperature inside of the rod, K; t is time, s; l is the rod conductivity factor, W/(m×K); c is the specific heat capacity, J/(kg·K); r is the rod material density, kg/m^{3}; and q is the specific power density of one rod, W/m^{3}.
The problem has been solved for a cylindrically-shaped fuel element (Fig. 2).
where T (r, z, 0) is the temperature of the rod points with the coordinates (r, z) at the time τ = 0;
ambient air temperature Taa = 592 K;
third-order boundary conditions at the rod ends (Fig. 2)
q_{1}(r, –l, τ) = 0, rϵ[0, R_{0}], (5)
q_{2}(r, l, τ) = 0, rϵ[0, R_{0}], (6)
where q_{1}(r, –l, t) is the heat flux at the rod’s lower end at the point with the coordinate r at the time t, W/m^{2}; and q_{2}(r, l, t) is the heat flux at the rod’s upper end at the point with the coordinate r at the time t, W/m^{2};
third-order boundary conditions on the rod’s side surface
where q_{3}(R_{0}, z, t) is the heat flux on the rod’s side surface, W/m^{2} (see Fig. 2); a is the coefficient of heat exchange with the fluid, W/(m^{2}·K); T (R_{0}, z, t) is the temperature of the rod side surface points at the time t, K; and T_{aa} is the ambient temperature, K.
We shall finally formulate in a cylindrical system of coordinates the initial boundary value problem required to obtain the analytical solution:
(8)
T (r, z, 0) = T_{0}, (9)
|T (0, z, t)| < ∞, (10)
(11)
T_{z} (r, –l, t) = T_{z} (r, l, t) = 0, (12)
(13)
where T = T (r, z, t) is the temperature of the rod points with the coordinates (r, z) at the time t, K;
f (r, z) = q (r, z)/(c·r) is the ratio of the specific power density depending on the coordinates r, z to the specific heat capacity and density of uranium dioxide, K/s; a^{2} is the thermal conductivity, m^{2}/s; R_{0} is the fuel element radius, m; l is the fuel element half-length, m; h = a/l is the ratio of the coefficient of heat exchange with the fluid to the conductivity coefficient, 1/m; a is the coefficient of heat exchange with the fluid, W/(m^{2}·K); l is the conductivity coefficient, W/(m·K); f_{0} = q_{0}/(c·r) is the ratio of the maximum specific power density in the fuel element to the specific heat capacity and density of uranium dioxide, K/s; d is the effective addition along the reactor height, m; R is the current radius of the fuel element position relative to the core center, m; and R_{e} is the reactor radius with regard for the effective addition, m.
Investigation results
The following assumptions were used as part of the analytical calculation.
Since the ratio R/R_{p} changes slightly within one rod, the value of the function J_{0}(2,405·R/R_{p}) can be conveniently presented as the parameter K, constant for each rod, the value of which is defined by the rod position inside of the core.
Since the radius of the fuel element is much smaller than its height and the heat transfer along the axis z (see Fig. 2) is negligibly small as compared with the heat release over the fuel element’s side surface, it can be assumed that the thermal energy distribution along the rod height is time-constant and agrees with the distribution defined at the initial calculation stage.
With regard for these assumptions, the analytical solution was obtained in the following form:
(14)
where x_{i} is the roots of the equation R_{0}·h·J_{0}(z) – x·J_{1}(z) = 0; µ_{0}_{i} = x_{i}^{2}/R_{0}^{2} is the ratio of the equation root x_{i} to the fuel radius, 1/m^{2}; and J_{0}, J_{1} are the zero- and first-order Bessel functions.
Based on this solution, a program was compiled in the MathCAD package. The test program has shown that the maximum temperature value for the given fuel element changes insignificantly with seven and more terms in the series. To obtain the temperature field with an accuracy of up to one degree, it is enough to take nine series terms, and 15 terms are enough when the accuracy is up to 0.1 degree.
The program compiled based on the analytical solution of the conductivity equation can be used to determine the temperature field in steady-state and non-steady-state modes of the fuel element operation. The convergence of the solution is presented for the steady-state fuel operating mode which occurs at t > 100 s (Fig. 3).
Calculated value of the maximum fuel temperature versus the number of the series terms for the system of equations solved.
https://binary.pensoft.net/fig/341995
The proposed analytical solution, with regard for the assumptions made, is accurate and can be therefore used to verify numerical programs built in other numerical simulation packages.
To investigate the influence of adjustment coefficients in numerical packages on the fuel temperature field determination accuracy, we shall solve this problem in the Comsol Multiphysics package. The number of the grid nodes, the grid node irregularity, and the grid type (a triangular or rectangular one for 2D problems, and a tetrahedron, hexahedron, pyramid or prism for 3D problems) can be used as adjustment coefficients when solving a steady-state problem for the temperature field determination (Shcherbakova and Gorbunov 2008, Gorbunov 2011a, 2011b, Palmer and Swanson 1961). The number of time steps is used as one of the adjustment coefficients in the solution of time-dependent problems. The predetermined values of adjustment coefficients influence the accuracy of the obtained result and the computer resources in the calculation.
We shall consider an axisymmetrical model of a fuel element with a triangular grid (Fig. 4).
Fuel temperature field determination grid for the axisymmetrical model solution.
https://binary.pensoft.net/fig/341996
The temperature field calculated in the Comsol Multiphysics package agrees with the results of the analytical solution in MathCAD (see Tables 1, 2). The results of the fuel temperature field determination were analyzed using a triangular grid which had 6961 nodes. A 2D axisymmetrical model, in which a one fourth part of the fuel element was built, was considered as the model.
Results of comparing the temperature field calculation with the radial coordinate r = 0 m.
Current axial coordinate z, m
Maximum fuel temperature calculated in Comsol, K
Maximum fuel temperature obtained in MathCAD, K
0
1038.5
1038.4
0.3
1025.74
1025.3
0.6
986.67
986.79
0.9
925.50
925.06
1.20
843.20
843.78
1.50
748.40
747.71
1.84
628.80
628.11
Results of comparing the temperature field calculation with the axial coordinate z = 0 m.
Current radial coordinate r, m
Maximum fuel temperature calculated in Comsol, K
Maximum fuel temperature obtained in MathCAD, K
0
1038.5
1038.4
0.0006
1029.74
1029.4
0.0012
1002.32
1002.4
0.0018
957.97
957.42
0.0024
893.77
894.42
0.003
813.80
813.42
0.0038
677.50
677.41
The influence of the number of the grid nodes on the maximum fuel temperature value determination has been established. The numerical experiment, in which the number of the nodes varied in a range of 21 to 94657, has shown that the number of the nodes in this problem does not influence the accuracy of the maximum temperature calculation result.
To take into account the influence of thermal hydraulics on the fuel temperature fields requires the use of 3D fuel models which have a prism-shaped grid. Comsol Multiphysics makes it possible to transform a 2D- model into a 3D model through the grid extrusion. An example of a 2D grid is shown in Fig. 5.
When a 3D model is developed, the grid extends in the radial direction (Fig. 6). In this case, the number of the grid elements influences the accuracy of the maximum fuel temperature determination result: with the coarsest grid division into 40 elements, it is enough to divide the grid into 30 axial layers to get a result with an uncertainty of one degree. If the number of the layers is greater, the uncertainty is larger due to the computational error accumulation.
With more grid nodes in the fuel element cross-section, it is enough to have 10 axial layers to get a result with an uncertainty of one degree (Fig. 7).
Uncertainty of the result versus the number of the grid layers.
https://binary.pensoft.net/fig/341999Conclusion
An accurate analytical solution of an axisymmetrical model has been obtained in a 2D statement of the problem to determine the fuel element temperature field with the assumptions that the physical properties of fuel do not depend on temperature and that the uranium-dioxide rod does not have a cladding.
The mathematical model based on the proposed analytical solution has been compiled in the MathCAD package and can be used to verify numerical models of the fuel temperature field determination.
ReferencesAgranatVMalinMPioroIAbdullahRPerminovVA (2015) BateninVMBelyaevIAListratovYaIPyatnitskayaNYuRazuvanovNGSviridovVGSviridovEVFrikPG (2017) Peculiarities of the heavy liquid metal heat exchange in nuclear power plants of a new generation. Proc. of the 13^{th} International Scientific and Practical Conference on Nuclear Power. Safety, Efficiency, Endurance. Sevastopol. SevNTU Publ., 79–81. [in Russian]DementyevBA (1990) Nuclear Power Reactors. Moscow. Energoatomizdat Publ., 352 pp. [in Russian]DolgovA (2016) Effective fuel solutions using SNF reprocessing. International Forum “AtomExpo”, Moscow, JSC Tvel, May 31.GoltsevAODavydovaGBDavidenkoVD (2009) Influence of the neutron flux depression in the RBMK cell on the magnitude of the maximum and average fuel temperature.GorbunovVA (2011a) Experience in using the software complex at Ivanovo State Power Engineering University named after V.I. Lenin.GorbunovVA (2011b) Predicting the accuracy of results in solving heat transfer problems based on neural network technologies.KartashovKVBogoslovskayaGP (2012) Calculations for the optimization of geometrical and operating parameters of VVER-SKD fuel assemblies for different modes of the reactor operation with supercritical water parameters. Izvestiya Vysshykh Uchebnykh Zavedeniy.KolpakovGNSelivanikovaOV (2009) LeskinSTShelegovASSlobodchukVI (2011) LoginovVS (2009) New Generation Codes (2018) New Generation Codes. http://www.ibrae.ac.ru/contents/68/ [accessed Jun 10, 2018] [in Russian]PalmerLDSwansonLL (1961) Measurements of Heat Transfer coefficients, Friction Factors and Velocity Profiles for Air Flowing Parallel to Closely Spaced Rods. Proc. of the Conf. of ASME, AIChE and Inst. Ch. E. on Internat. Developments in Heat Transfer. ASME, part 3, August: 535–542.PerimovRRSorokinGASorokinaTV (2004) Modeling of the thermal-engineering reliability of a fuel rod with different options of the power density and temperature variation.ShcherbakovaTSGorbunovVA (2008) A study into the solution of accuracy problems with second-order boundary conditions of heating in the Phoenics multipurpose computing complex. Proc. of the All-Russian Student Olympiad, Scientific and Practical Conference and Exhibition of Student, Postgraduate and Young Scientist Works “Energy and Resource Saving, Power Supply, Nonconventional and Renewable Energy Sources”. UGTU-UPI Publ., Yekaterinburg, 174–175. [in Russian]StarkovVAMarikhinNYu (2013) Procedures and program of the stationary temperature field calculation in the system of multi-zone cylindrical fuel rods. Izvestiya Vysshykh Uchebnykh Zavedeniy.UshakovPASubbotinVI (1972) Approximate calculations of the hydrodynamic characteristics of a turbulent fluid flow in ring channels.VelesyukAMorgunovI (1990) CFD codes: problems and prospects in nuclear power engineering. Atomny Ekspert.VishnyakovaADGulinaOMSalnikovNL (2015) The possibility of using the neural network apparatus for predicting erosive and corrosive wear of the NPP equipment. Izvestiya Vysshykh Uchebnykh Zavedeniy.ZhukovAVKuzinaYuASorokinAP (2009) Systematization of the studies on heat exchange in fuel assemblies and selected problems of liquid metal cooling. Izvestiya Vysshykh Uchebnykh Zavedeniy.
* Russian text published: Izvestiya vuzov. Yadernaya Energetika (ISSN 0204-3327), 2019, n. 2, pp. 174–184.