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 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 (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):

Q (r, z) = Q₀J₀(2,405r/R_p)cos(πz/H_p), (1)

where Q₀ is the maximum value of the power density at the reactor center, MW; J₀ 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).

Figure 1. 

Power density distribution radially Q_r and axially Q_z in a homogeneous cylindrical core.

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.

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.

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³; and q is the specific power density of one rod, W/m³.

The problem has been solved for a cylindrically-shaped fuel element (Fig. 2).

Figure 2. 

Geometrical parameters of the rod and the boundary conditions influencing the solution.

Geometrical parameters:

rod half-length l = 3.68/2 = 1.84 m;

rod radius R₀ = 0.0076/2 = 0.0038 m. (3)

Initial and boundary conditions influencing the solution:

• initial rod temperature

Т (r, z, 0) = T₀ = 592 K, rϵ[0, R₀], 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;

• ambient air temperature Taa = 592 K;
• third-order boundary conditions at the rod ends (Fig. 2)

q ₁(r, –l, τ) = 0, rϵ[0, R₀], (5)

q ₂(r, l, τ) = 0, rϵ[0, R₀], (6)

where q₁(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²; and q₂(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²;

• third-order boundary conditions on the rod’s side surface

q ₃(R₀, z, τ) = α·(Т (R₀, z, τ) – T_аa), zϵ[–l, l], (7)

where q₃(R₀, z, t) is the heat flux on the rod’s side surface, W/m² (see Fig. 2); a is the coefficient of heat exchange with the fluid, W/(m²·K); T (R₀, 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₀, (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² is the thermal conductivity, m²/s; R₀ 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²·K); l is the conductivity coefficient, W/(m·K); f₀ = q₀/(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₀(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₀·h·J₀(z) – x·J₁(z) = 0; µ_0i = x_i²/R₀² is the ratio of the equation root x_i to the fuel radius, 1/m²; and J₀, J₁ 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).

Figure 3. 

Calculated value of the maximum fuel temperature versus the number of the series terms for the system of equations solved.

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

Figure 4. 

Fuel temperature field determination grid for the axisymmetrical model solution.

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.

Table 1.Download as CSV 

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

Table 2.Download as CSV 

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

Figure 5. 

A triangular 2D grid for the fuel element.

Figure 6. 

A 3D grid for the fuel element.

Figure 7. 

Uncertainty of the result versus the number of the grid layers.

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