Consistent neutron-physical and thermal-physical calculations of fuel rods of VVER type reactors

For modeling the isotopic composition of fuel, and maximum temperatures at different moments of time, one can use different algorithms and codes. In connection with the development of new types of fuel assemblies and progress in computer technology, the task makes important to increase accuracy in modeling of the above characteristics of fuel assemblies during the operation. Calculations of neutronphysical characteristics of fuel rods are mainly based on models using averaged temperature, thermal conductivity factors, and heat power density. In this paper, complex approach is presented, based on modern algorithms, methods and codes to solve separate tasks of thermal conductivity, neutron transport, and nuclide transformation kinetics. It allows to perform neutron-physical and thermal-physical calculation of the reactor with detailed temperature distribution, with account of temperature-depending thermal conductivity and other characteristics. It was applied to studies of fuel cell of the VVER-1000 reactor. When developing new algorithms and programs, which should improve the accuracy of modeling the isotopic composition and maximum temperature in the fuel rod, it is necessary to have a set of test tasks for verification. The proposed approach can be used for development of such verification base for testing calculation of fuel rods of VVER type reactors


Introduction
Simplified models of calculation are often used in traditional approach to modeling neutron-physical characteristics of fuel rods and isotopic composition in process of nuclear fuel burning [1][2][3][4][5].For instance, in many cases neutron-physical characteristics of fuel depending on temperature are determined by use of average temperature of fuel rod, while the existing dependencies of the neutron cross sections on temperature are neglected.In thermal-physical calculations of the radial distribution of fuel temperature, constant average value of thermal conductivity factor is often used.However this factor depends really both on the temperature and on fuel burnup, which change with radius of fuel rod.
In this paper, mathematical models of consistent neutron-physical and thermal-physical calculation are presented.They are realized in codes designed in MEPhI and "Kurchatov Institute".This approach was applied to calculation investigation of fuel element of the power reactor of VVER-1000 type.Short description of benchmarks for verification of obtained results is given.It is demonstrated how the account of dependence of the thermal conductivity factor on temperature affects on the value of temperature in the center of fuel rod for the VVER type reactor.

Algorithms and methods
In order to obtain the dependence of temperature on radius of cylindrical fuel rod it is necessary to solve the equation of thermal conductivity with specified boundary conditions.If density of energy release and thermal conductivity factor are constant independent on radius of fuel rod, this problem has an analytical solution.However, thermal conductivity factor of oxide nuclear fuel depends on temperature and on fuel burnup.Energy release density and fuel burnup vary nonuniformly with radius of fuel rod in process of fuel burnup.Swelling of fuel pellet under deep burnup results in reduction of the gap between fuel pellet and cladding.
To value the contribution of each of these effects in change of maximum temperature of fuel rod, it is necessary to be able to calculate fuel burnup and change of isotopic composition, as well as distribution of energy release density on radius of fuel rod.For that it is necessary to solve multigroup equation of transport or diffusion of neutrons and system of equations of isotope kinetics.Remember that usually these equations are solved one after another.Firstly transport equation is solved with macroscopic constants calculated on the base of the known isotopic composition.Then one-group cross sections of all isotopes belonging to equations of isotope kinetics are defined.Hereon the system of the equations of isotope kinetics is solved and new isotopic composition for the given moment of time is calculated.After that, we should solve once again transport equation since change of isotopic composition can result in change of macroscopic constants and of neutron spectrum.
Let's consider mathematical settings of the above described tasks.

Neutron diffusion equation
To obtain the one-grup cross sections and neutron fluxes, it is necessary to simulate neutron spectrum and hence to solve the neutron transport equation or diffusion equation.Usually multigroup approach is used for this purpose, in which one must solve a system of multigroup equations like (1).All symbols are traditional, g is the number of the group while G is a full number of groups.
The solution of these equations in particular geometry with specified boundary conditions allows to determine the neutron spectrum that can be used to obtain the one-group constants.To determine one-group neutron fluxes it is necassaru to calculate additionaly the spatial distribution of energy release density and know the full energy release in a simulated system.
We should keep in mind that diffusion coefficients ) (r D & and macroscopic cross-sections ) (r g & ¦ depend not only on the spatial variable, but also on the temperature in the spatial domain.In order to make correct system of multigroup equations with a purpose to prepare prepare one-group cross sections and fluxes, one need to set or be able to calculate the temperature in spatial areas.

Equation of isotope kinetics
System of equations of isotope kinetics is formulated for vector of isotope concentrations of nuclides ) (t i ρ .Isotope concentrations in each spatial zone with number i are assumed to be the same in all points of spatial sone.
Each equation of the system of equations of isotope kinetics include the following values in addition to concentrations.
i l σ -one-group cross section of incineration of isotope with number l in spatial zone with number i due to nuclear reactions on l-th isotipe.Vector of concentrations changes in process of burning of nuclear fuel.Concentrations of isotopes may vary in different spatial locations of the fuel core.For example, area with modified microstructure, so-called rim-zone appeares in the outer layer of the fuel rod of thermal-neutron reactors.Its formation is associated with high concentration of fission products [6][7][8].Therefore, application of constant concentrations over the whole fuel rod can lead to distortion of the simulation results.Along with isotope concentrations, other properties of fuel composition can vary with radius within fuel rod in the process of burnup and can influence on space distribution of the temperature.These variation of temperature don't take into account in the "classic" approach, and the process of formation of rim-zone can't be calculated.
In order to solve a system of equations ( 2 (2), the new value of the concentration vector ) (t i ρ is obtained, which corresponds to a certain moment of time.

Equation of thermal conductivity
In "classic" approach, the factor of thermal conductivity is considered to be independant on the radius within fuel rod, that allows to obtain an analytical decision, if radial distribution of the energy release density qv(r) is known.
However, really thermal conductivity factor of oxide nuclear fuel depends on the temperature T of this fuel and fuel burnup B [9] (Fig. 1).
where λ(T(r),ȼ) is thermal conductivity factor of fuel matherial depending on temperature and burnup, qv is energy release density.
Several methods can be used to solve equation ( 3).One of them consists in dividing the fuel rod into annular zones.Factors of thermal conductivity remain constant within each zone : λ(T(r)) = {λ (Ti)} = {λi}.
Then the thermal conductivity equation becomes linear within each zone.However, in this case we need to use iterative procedure.Each iteration specifies the value of thermal conductivity factor in each layer.This allows to consider the dependence of thermal conductivity factor on temperature.Similarly, we can proceed with the dependence of thermal conductivity factor on fuel burnup : in each step of burnup to take from the table the new value of the factor of thermal conductivity and insert it into the equation to find the temperature.

Complex model of neutron-physics calculation of fuel cell
All of the above described relations and methods of their descriptions can be represented in a model of consistent thermal physical and neutron-physical calculations.In Fig. 2 the following designations are used: Eq.1 -transport equation, Eq.2 -equation of isotope kinetics, Eq.3 -equation of thermal conductivity.DB1 (data base) -dependence of microscopic cross sections from neutron energy, account of shielding, account of the temperature dependence of microscopic cross sections, DB2 -library of decay constant.Ȉ, ı are macroscopic and microscopic cross sections of interaction of neutrons with nuclei of the medium.Φabsolute neutron flux (matched with the power and onegroup cross sections), ĳ -multi-group neutron flux (the result of solving the transport equation).ȡ is the vector of the nuclide concentration (the result of solving the system of equations of isotopic kinetics).qv -energy release density.B -fuel burnup.T -temperature (the solution to the thermal conductivity equation).A1, A2, A3, A4, A5 -the algorithms that implement the ratio between the required values based on the model and experimental data : A1.

Temperature calc.
Neutron physics calc. 1.We assign initial value of energy release density (qv) and thermal conductivity factor Ȝ. Then iterative procedure is made.Firstly, thermal conductivity equation (Eq.3) is solved.The result is temperature value in the cell.In case the cell is divided into zones, we obtain values of the temperature in each zone of the cell.Secondly we calculate new value of the thermal conductivity factor Ȝ, corresponding to new value of the temperature.These two steps can be repeated many times to obtain required accuracy of the temperature.Simultaneously we get more accurate Ȝ value.
2. This temperature value is used as input parameter for database (DB1) of the temperature dependencies of the microscopic cross sections ı.Using these data of ı together with initial values of nuclide concentrations (ρ0), we obtain value of macroscopic cross section (Σg) by means of algorithm (A1).
3. With this macroscopic cross section we solve the transport equation (Eq.1).As a result, due to algorithm A2, we obtain group neutron fluxes in zones (Φg), that is to say spatial dependency of group neutron fluxes.At the same time, we save new values of the microscopic cross sections (Φ σ).We use the group neutron fluxes to calculate new value of energy release.We repeat the steps 1-3 as iterations over the energy release density (algorithm A3).Herewith we make more accurate the temperature distribution and take into account the temperature dependency of the thermal conductivity factor.
4. Using group neutron fluxes along with decay constants (DB2) and new values of the microscopic cross sections, we solve the equation of isotope kinetics (Eq.2).As a result we obtain isotope concentrations on the next step of burnup (ρk).With these new values of isotope concentrations and values of microscopic cross sections we obtain new value of macroscopic cross sections (algorithm A1).The steps 3-4 can be repeated several times as iterations.
We can perform iterations over the steps 1-4 to achieve required accuracy.Herewith we shall take into account the dependency of Ȝ on fuel burnup.
Thereby, new complex approach allows to take into account: -spatial distribution of energy release to determine neutron flux; -dependence of the diffusion factor and macroscopic cross section not only on spatial value, but also on the temperature in given spatial area; -presence of spatial-nonuniform structure in fuel rod (rim-zone); -dependence of thermal conductivity factor of oxide fuel on the temperature of fuel; -reduction of thermal conductivity factor of fuel in process of burnup.
In order to appreciate efficiency of the presented approach to calculation of temperature distribution over the radius of fuel rod, it is necessary to answer the following questions: 1. How do estimations of the maximum temperature of fuel change when turning from previous approaches to new approaches?
2. Is it necessary to change traditional schemes of calculation of isotopic composition in process of fuel burnup by including in them block of calculation of the temperature distribution?
3. How should we calculate fuel temperature if thermal conductivity factor depends on the temperature and fuel burnup?

Algorithms and codes used for calculations
Thermal-physical calculations in this paper were made using codes HEATING7, UNK_teplo, TempR_5.
HEATING7 is multipurpose code for solution of thermal conductivity tasks [10,11].It is included into code complex SCALE [12].However it is not integrated code of the SCALE because it doesn't have standard format of input data and doesn't use in functional sequences of the SCALE.
Heatinh7 can solve stationary and nonstationary thermal conductivity tasks in one-, two-or threedimensional square, cylindrical or spherical coordinates.Calculating model can include different materials.Thermal conductivity, density and energy release of each material can depend on time and temperature.Energy release can depend on time, temperature, space.Thermal conductivity can be anisotropic.Phase change can be taken into account.Thermal properties of materials can be included by user or can be taken from the own library of material thermal properties.
The energy release may depend on time, temperature, position, and boundary temperatures, which in turn may depend on time and coordinates.The boundary conditions may be assigned on the outer surface of the model or on the surface between parts of the model.They can be specified by temperatures or by any combination of heat flux, forced convection, natural convection and radiation.The parameters of the boundary conditions can depend on time and temperature.The spatial lattice for calculations can be arbitrary for any variable.HEATING uses the scheme of memory distribution during runtime, which allows to avoid the need of re-calculations and to meet the memory requirements for any tasks.The code UNK_teplo can solve the stationary thermal conductivity tasks in one-dimensional cylindrical coordinates with different materials.Density and thermal conductivity materials should be assigned for each layer.Numerical methods are used for calculation of temperature in each layer.
The code TempR_5 allows to calculate temperatures in cylindrical fuel element with piecewise constant thermal conductivity factor.Temperature distribution in each layer is obtained from analytical decision of thermal conductivity equation with assigned constant values of radius, thermal conductivity factor, and temperature.The results of calculation are average values of temperature over the layer, temperatures in points between layers and in the center of fuel element.
UNK_teplo is code for solution of thermal conductivity tasks.Code is included into program complex UNK [10,13] UNK_teplo can solve the stationary thermal conductivity tasks in one-dimensional cylindrical coordinates.Calculation model may include various materials.Thermal conductivity and density of each material is specified separately for each layer.Thermal properties of materials can be entered by the user.
Numerical methods are used to obtain the temperature values in each layer.The system of equations compiled according to the finite-difference scheme, then the equations are solved by sweep method.
The calculation of temperatures in the UNK complex is made in coordination with calculation of fuel burnup.Firstly, macro constants are prepared and the temperatures are calculated in different zones of the given geometry.After that, standard module of UNK starts to operate.It calculates the value of Keff , the new macro constants, energy release in each zone.Then temperatures in zones are recalculated.
The code TempR_5 allows to calculate temperatures in cylindrical fuel element with piecewise constant thermal conductivity factor.Temperature distribution in each layer is obtained from analytical decision of thermal conductivity equation with assigned constant values of radius, thermal conductivity factor, and temperature.The results of calculation are average values of temperature over the layer, temperatures in points between layers and in the center of fuel element.

Results of test tasks calculations
Application of the proposed approach is illustrated by calculations of the temperature distribution along the radius of the fuel rod of the VVER-1000 reactor.These calculations can be considered as test tasks (benchmarks) to verify the developed algorithms.Test tasks will help to answer the above questions and to indicate the expediency of the complexity of the model for calculation of the fuel rod.We present below several results illustrating the capabilities of the algorithms.The tasks were solved in one-dimensional cylindrical geometry.
Model of the elementary cell of the VVER-1000 is taken as base model.The diameter of the central hole of the fuel pellet was taken as 1.5 mm, the external diameter of the pellet 7.55 mm; internal and external diameters of the cladding 7.72 and 9.17 mm; the pitch of the triangular lattice of rods 12.75 mm, temperature of coolant 575.7 K.For calculations, fuel rod was divided on 11 radial layers.
The temperature of the outer surface of the cladding was the same for all calculations.It was equal to the temperature of the coolant 575,7 K.The thermal conductivity factor was taken as 0.35 W/mK for the center hole and the gas gap between the fuel rod and cladding, and as 20.42 W/mK for cladding.These values were the same for all calculations.
First of all it is necessary to illustrate the advantage of the developed technique, which allows to take into account the dependence of thermal conductivity factor on temperature and on fuel burnup.Fig. 3 shows the results of calculations of the temperature inside the fuel element for four cases.The case 1 (curve 1) corresponds to the traditional approach, in which λ was taken to be constant throughout the fuel rod.The λ value was calculated for the average temperature over the rod as 1000 K.The case 2 (curve 2) corresponds to the application of the developed technique, however only dependence λ (T) for burnup B = 0 was considered.In cases 3 and 4, the dependence λ(T, B) was takes into account.Calculations were performed for B = 60 (curve 3) and 120 MWÂday/kgU (curve 4).For all four cases, it was assumed that no gas gap was there between fuel rod and cladding, and energy release density qv = 388 W/cm 3 was constant over the fuel rod.
The results presented in Fig. 3 can be easily interpreted using data on λ(T, B) from Fig. 1.In case 1 for λ = const, we obtain curve 1 in Fig. 3.The temperature decreases from the center to the periphery of the fuel rod.Peripheral temperature is determined by the coolant temperature of 575 K.
In case 2 when we consider λ(T) at B = 0 (the upper curve 1 in Fig. 1), the λ increases sharply due to temperature decrease in going from the center of the fuel rod to its periphery.Since the peripheral temperature is specified, the temperature at the center is lower than in case 1.The curve 2 lies below the curve 1.
In case 3 (the dependence λ(T, B) corresponds to the curve 3 in Fig. 1) the λ values are much less than in cases 1 and 2. The increase of the λ in going from the center to the periphery of fuel rod is slower than in cases 1 and 2. Since the energy release density is the same in all cases and peripheral temperature is specified, the temperature in the center of the rod is higher than in case 1.The curve 3 in Fig. 3 lies above the curve 1.
In case 4 (curve λ (T, B) is more flatter than curve 4 in Fig. 1) similar arguments show that the temperature in the center of the rod is higher than in case 3. The curve 4 in Fig. 3  These data demonstrate the efficiensy of the developed methods of calculation of the temperature in the fuel element of the VVER reactor.Corrections to the value of the maximum temperature can reach 10-15%.
Let's now to turn to the results showing the capabilities of the algorithms.
In tasks 1 and 2 (Fig. 4 and 5), the uniform energy release density of 388 W/cm 3 was taken over the whole fuel rod.Thermal conductivity factor of the fuel rod λ was also taken constant.We considered three values of λ = 1.7 (curves 1 in figures 4 and 5), 2.1 (curves 2) and 3.0 W/mK (curves 3).The difference between the results shown in Figs. 4 and 5 was that Fig. 4 corresponded to fresh fuel with the gap between fuel rod and cladding, while Fig. 5 corresponded to the "swollen" fuel.There was no gap in this case, and fuel was in thermal contact with cladding.
Results presented in Figs. 4 and 5 demonstrate apparent dependence of the temperature curves for different values of λ.The lower λ, the higher is the curve T(R).All three curves converge on the outer surface of the fuel rod.Fig. 4 clearly demonstrates sudden change in temperature in the gap between the outer surface of the fuel rod and cladding.In task 3 (Fig. 7), thermal conductivity factor of fuel λ depended on the radius within the fuel rod.This dependence was the same for curves 1 and 2. The λ value increased smoothly from 2.1 W/mK in the central layer to 6.5 W/mK in the outer layer of fuel rod.The energy release density was taken to be constant over the volume of the rod and was equal to 388 W/m 3 for curve 1 and 370 W/m 3 for curve 2. Curve 1 corresponded to the presence of the gap between fuel rod and cladding, while curve 2 corresponded to the absence of the gap, and fuel was in thermal contact with cladding.The result is similar to the curves of Fig. 6, however in Fig. 7, the curves are more flat due to the increase of λ from the center to the periphery.

Conclusions
Thus, we present a consistent approach to complex calculation of nuclear fuel for VVER reactors.It is based on modern software and allows to obtain detailed distribution of power density and temperatures [14] in the fuel rod.The proposed approach can serve further for development of verification base for testing of calculations of the fuel rods of VVER reactors.

)
group cross section of appearance of isotope with number l in spatial zone with number i due to nuclear reactions on m-th isotipe.i φ -one-group neutron flux in i-th spatial zone.l λ -decay constant of l-th isotope l k → λ -decay constant of k-th isotope resulting in appearance of l-th isotope.
. It is also necessary to know the decay constants l λ and l k→ λ for each nuclide and initial concentrations of isotopes i l ρ .After the solution of the equation system

Fig. 2 .
Fig. 2. Diagram of consistent thermal physical and neutronphysical calculations Below we propose 4-steps sequence of calculations of neutron-physical characteristics of fuel rods based on diagram Fig. 2. It allows to simulate physical and neutron-physical characteristic of fuel more correctly and increases accuracy of the results.1.We assign initial value of energy release density (qv) and thermal conductivity factor Ȝ. Then iterative procedure is made.Firstly, thermal conductivity equation (Eq.3) is solved.The result is temperature value in the cell.In case the cell is divided into zones, we obtain values of the temperature in each zone of the cell.Secondly we calculate new value of the thermal conductivity factor Ȝ, corresponding to new value of the temperature.These two steps can be repeated many times to obtain required accuracy of the temperature.Simultaneously we get more accurate Ȝ value.