Print
Development of the model to determine the fuel temperature field in a two-dimensional problem statement*
expand article infoVladimir A. Gorbunov, Natalya B. Ivanova, Nikita A. Lonshakov, Yaroslav V. Belov
‡ Ivanovo State Power University, Ivanovo, Russia
Open Access

Abstract

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.

Keywords

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

Introduction

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) = Q0J0(2,405r/Rp)cos(πz/Hp), (1)

where Q0 is the maximum value of the power density at the reactor center, MW; J0 is the zero-order Bessel function; Re is the effective radius, m; He 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 Qr and axially Qz 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.

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/() + q/(), (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/m3; and q is the specific power density of one rod, W/m3.

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 R0 = 0.0076/2 = 0.0038 m. (3)

Initial and boundary conditions influencing the solution:

  • initial rod temperature

Т (r, z, 0) = T0 = 592 K, rϵ[0, R0], 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 1(r, –l, τ) = 0, rϵ[0, R0], (5)

q 2(r, l, τ) = 0, rϵ[0, R0], (6)

where q1(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/m2; and q2(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/m2;

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

q 3(R0, z, τ) = α·(Т (R0, z, τ) – Tаa), zϵ[–l, l], (7)

where q3(R0, z, t) is the heat flux on the rod’s side surface, W/m2 (see Fig. 2); a is the coefficient of heat exchange with the fluid, W/(m2·K); T (R0, z, t) is the temperature of the rod side surface points at the time t, K; and Taa 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) = T0, (9)

|T (0, z, t)| < ∞, (10)

(11)

Tz (r, –l, t) = Tz (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; a2 is the thermal conductivity, m2/s; R0 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/(m2·K); l is the conductivity coefficient, W/(m·K); f0 = q0/(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 Re 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/Rp changes slightly within one rod, the value of the function J0(2,405·R/Rp) 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 xi is the roots of the equation R0·h·J0(z) – x·J1(z) = 0; µ0i = xi2/R02 is the ratio of the equation root xi to the fuel radius, 1/m2; and J0, J1 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.

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

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.

Conclusion

  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.

References

  • Batenin VM, Belyaev IA, Listratov YaI, Pyatnitskaya NYu, Razuvanov NG, Sviridov VG, Sviridov EV, Frik PG (2017) Peculiarities of the heavy liquid metal heat exchange in nuclear power plants of a new generation. Proc. of the 13th International Scientific and Practical Conference on Nuclear Power. Safety, Efficiency, Endurance. Sevastopol. SevNTU Publ., 79–81. [in Russian]
  • Dementyev BA (1990) Nuclear Power Reactors. Moscow. Energoatomizdat Publ., 352 pp. [in Russian]
  • Dolgov A (2016) Effective fuel solutions using SNF reprocessing. International Forum “AtomExpo”, Moscow, JSC Tvel, May 31.
  • Goltsev AO, Davydova GB, Davidenko VD (2009) Influence of the neutron flux depression in the RBMK cell on the magnitude of the maximum and average fuel temperature. Izvestiya Tomskogo Politekhnicheskogo Universiteta 314(4): 5–7. [in Russian]
  • Gorbunov VA (2011a) Experience in using the software complex at Ivanovo State Power Engineering University named after V.I. Lenin. ANSYS Advantage (Russkaya Redaktsiya) 15: 38–39. [in Russian]
  • Gorbunov VA (2011b) Predicting the accuracy of results in solving heat transfer problems based on neural network technologies. Promyshlennaya Energetika, 8: 48–52. [in Russian]
  • Kartashov KV, Bogoslovskaya GP (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. Yadernaya Energetika 2: 3–11.
  • Kolpakov GN, Selivanikova OV (2009) Designs of Fuel Elements, Channels and Cores of Power Reactors. Tomskiy Politekhnicheskiy Universitet Publ., Tomsk, 118 pp. [in Russian]
  • Leskin ST, Shelegov AS, Slobodchuk VI (2011) Physical Features and Design of the VVER-1000 Reactor. NIYaU MIFI Publ., Moscow, 116 pp. [in Russian]
  • Loginov VS (2009) Approximate Methods of Thermal Calculation for Active Elements of Electrophysical Plants. Fizmatlit Publ., Moscow, 273 pp. [in Russian]
  • Palmer LD, Swanson LL (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.
  • Perimov RR, Sorokin GA, Sorokina TV (2004) Modeling of the thermal-engineering reliability of a fuel rod with different options of the power density and temperature variation. Promyshlennaya Teplotekhnika, 26(5): 150–153. [in Russian]
  • Shcherbakova TS, Gorbunov VA (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]
  • Starkov VA, Marikhin NYu (2013) Procedures and program of the stationary temperature field calculation in the system of multi-zone cylindrical fuel rods. Izvestiya Vysshykh Uchebnykh Zavedeniy. Yadernaya Energetika 1: 54–62. [in Russian]
  • Ushakov PA, Subbotin VI (1972) Approximate calculations of the hydrodynamic characteristics of a turbulent fluid flow in ring channels. Teplofizika Vysokikh Temperatur 10(5): 1025–1030. [in Russian]
  • Velesyuk A, Morgunov I (1990) CFD codes: problems and prospects in nuclear power engineering. Atomny Ekspert. Supplement to Atomnaya Energetika 8: 37–39. [in Russian]
  • Vishnyakova AD, Gulina OM, Salnikov NL (2015) The possibility of using the neural network apparatus for predicting erosive and corrosive wear of the NPP equipment. Izvestiya Vysshykh Uchebnykh Zavedeniy. Yadernaya Energetika 4: 61–71. [in Russian]
  • Zhukov AV, Kuzina YuA, Sorokin AP (2009) Systematization of the studies on heat exchange in fuel assemblies and selected problems of liquid metal cooling. Izvestiya Vysshykh Uchebnykh Zavedeniy. Yadernaya Energetika 4: 95–108. [in Russian]

* Russian text published: Izvestiya vuzov. Yadernaya Energetika (ISSN 0204-3327), 2019, n. 2, pp. 174–184.