Corresponding author: Vladimir A. Gorbunov ( gorbynov.w@mail.ru ) Academic editor: Boris Balakin
© 2019 Vladimir A. Gorbunov, Natalya B. Ivanova, Nikita A. Lonshakov, Yaroslav V. Belov.
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:
Gorbunov VA, Ivanova NB, Lonshakov NA, Belov YV (2019) Development of the model to determine the fuel temperature field in a twodimensional problem statement. Nuclear Energy and Technology 5(3): 219224. https://doi.org/10.3897/nucet.5.39318

Watercooled watermoderated 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 twodimensional 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 element, fuel temperature field, analytical solution, mathematical model to determine fuel temperature fields, verification of numerical calculations, influence of tuning coefficients, safety of fuel element heating
Operating safety of a nuclear reactor requires stringent limits to be placed on the maximum temperature of fuel (
The purpose of the nuclear reactor core is to generate thermal energy and transfer it from the fuel surface to the coolant. VVER1000 reactors use cylindrically shaped barerod 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 (
Q (r, z) = Q_{0}J_{0}(2,405r/R_{p})cos(πz/H_{p}), (1)
where Q_{0} is the maximum value of the power density at the reactor center, MW; J_{0} is the zeroorder 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.
The fuel element is a pressurized tube of a zirconiumniobium mixture which is filled with uranium dioxide pellets (
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 (
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 secondorder boundary conditions on the fuel element’s end faces and the thirdorder boundary conditions on the rod’s side surface (
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 thermalphysical coefficients with regard for the internal heat source (
The following assumptions are used for the problem solution:
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 cylindricallyshaped fuel element (Fig.
Geometrical parameters:
rod halflength l = 3.68/2 = 1.84 m;
rod radius R_{0} = 0.0076/2 = 0.0038 m. (3)
Initial and boundary conditions influencing the solution:
Т (r, z, 0) = T_{0} = 592 K, rϵ[0, R_{0}], zϵ[–l, l], (4)
where T (r, z, 0) is the temperature of the rod points with the coordinates (r, z) at the time τ = 0;
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};
q _{3}(R_{0}, z, τ) = α·(Т (R_{0}, z, τ) – T_{аa}), zϵ[–l, l], (7)
where q_{3}(R_{0}, z, t) is the heat flux on the rod’s side surface, W/m^{2} (see Fig.
We shall finally formulate in a cylindrical system of coordinates the initial boundary value problem required to obtain the analytical solution:
T (r, z, 0) = T_{0}, (9)
T (0, z, t) < ∞, (10)
T_{z} (r, –l, t) = T_{z} (r, l, t) = 0, (12)
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 halflength, 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.
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.
With regard for these assumptions, the analytical solution was obtained in the following form:
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 firstorder 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 steadystate and nonsteadystate modes of the fuel element operation. The convergence of the solution is presented for the steadystate fuel operating mode which occurs at t > 100 s (Fig.
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 steadystate problem for the temperature field determination (
We shall consider an axisymmetrical model of a fuel element with a triangular grid (Fig.
The temperature field calculated in the Comsol Multiphysics package agrees with the results of the analytical solution in MathCAD (see Tables
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 prismshaped 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.
When a 3D model is developed, the grid extends in the radial direction (Fig.
With more grid nodes in the fuel element crosssection, it is enough to have 10 axial layers to get a result with an uncertainty of one degree (Fig.