NUMERICAL MULTIGROUP TRANSIENT ANALYSIS OF SLAB NUCLEAR REACTOR WITH THERMAL FEEDBACK

The paper describes a new numerical code for multigroup transient analyses with a thermal feedback. The code is developed at Institute of Nuclear and Physical Engineering. It is necessary to carefully investigate transient states of fast neutron reactors, due to recriticality issues after accident scenarios. The code solves numerical diffusion equation for 1D problem with possible neutron source incorporation. Crank-Nicholson numerical method is used for the transient states. The investigated cases are describing behavior of PWR fuel assembly inside the spent fuel pool and with the incorporated neutron source for better illustration of the thermal feedback.


Introduction
Nuclear reactors are most reliable source of energy at present time.The advantage of new fast neutron reactor concepts is the possibility of the transmutation of the 238 U into the fissionable 239 Pu, so the worldwide resources of fissionable materials could be enlarged when the fast neutron reactors will be commercialized.However, it is necessary to carefully investigate the transient states of these new concepts, because the recriticality after accident scenario may be achieved [1].
For this purpose, a new numerical code has been developed within the frame of neutronic calculation group at the Institute of Nuclear and Physical Engineering (INPE).The code has been used as the basic step in understanding of transient state and fast reactor R&D and we refer to it as TRANSOS.The TRANSOS code is able to solve a time dependent neutron flux distribution of slab systems with respect to the thermal feedback in multigroup approximation.The power method with infinity norm (Manhattan norm) is used for a calculation of stationary state [2].The transient state is calculated by the Crank-Nicholson numerical method [3].It is possible to define the geometry in the form of coarse mesh intervals with material parameters assigned to each of the mesh interval.Point and linear neutron sources can be incorporated to the design with different fluency for each investigated energy group of neutrons.The boundary conditions have to be specified in conjunction with design parameters such as temperature of the system and power generation rate.Initial state of neutron flux should be given for successful run.The developed code can be used for academic purposes at the faculty and for investigation of simple systems.
This paper follows the previous work [4] and next sections describe calculation methodology of time de-pendent and in-dependent solution, accuracy of the methodology and analysis of two nuclear systems with different 235 U enrichment.Higher impact of the thermal feedback can be observed on microscopic cross section properties for low energy neutron interaction with matter, so the thermal systems are investigated to show this effect.

Stationary state (time in-dependent solution)
The code solves common neutron diffusion equation with internal sources for the stationary state given by [5] Lower index g in all equations represents the neutron energy group and D stands for the associated diffusion coefficient.The first term represents the leakage of neutrons from the system, the second term represents absorptions (radiation capture), the third represents scattering from the group to another group h, the fourth represents contributions of scattering from another group and the last term describes the contribution of fission source to specified energy group (χ g is the proportion of neutrons emitted by fission in group g, ν is an average yield of neutrons from the fission).The Box Scheme was used for the discretization of examined area and therefore (1) can be written in a numerical form (2) for two energy groups of neutrons [2,4], where Σ r represents the removal macroscopic cross section (3), k represents spatial node, h k is the length of the spatial interval, λ stands for the fundamental eigenvalue, Σ f is the fission macroscopic cross section, Σ s is the scattering macroscopic cross section between energy groups of neutrons, and D arbitrary coupling coefficient (3) [2].

Crank-Nicholson numerical method
The Crank-Nicholson method is often applied to the diffusion problems and achieves good numerical stability for the transient states.The finite difference method is used as for the stationary state.This method is able to solve heat equation and similar partial differential equations [3] where M represents material parameters and U is solution of the differential equation.Applying the finite difference method, the central node of the common diffusion equation can be rewritten numerically for the transient state as where the index k represents spatial node, index j represents step time ( U •,j stands for initial condition of the function at present time step, U •,j+1 is initial condition of the function at next time step), the h k length of the spatial interval, h t length of the time interval, f (h k ) is increment of the function U (particular solution).The variables and unknowns from ( 5) can be rewritten in the form (6) for the neutron diffusion equation ( 1) and (2).
where the last term S represents the external neutron source.The heat diffusion equation can be similarly rewritten to [6] where Kk = 2 , and where K represents material conductivity, ρ is material density, c stands for material heat capacity, q is power generation, T stands for temperature, E is an average released energy per fission in eV and parameter e is an electron charge.Accuracy of the above mentioned method is discussed in the next subsection.

Accuracy of the Crank-Nicholson method
The accuracy of the method for the initial and boundary conditions (more in Section 3) is presented in Fig. 1.The analytical solution of the differential equation for a homogenous slab system is sinus, so the problem can be approximated by function f (x) f (x) = a sin(bx) + c.In this case, the whole system is homogeneous subcritical nuclear system.It is necessary to insert the external neutron source because the neutron flux will never stabilize on non-zero value without it.The external neutron source is linear, located between 217-227 mm, with fluency 1.0 × 10 6 s −1 in thermal energy spectrum.Deviations between the modeled problem and analytical solution are presented in Table 1.Amplitude (a) and shape parameter (b) are in accordance with the analytical solution (Table 1).Deviance in location parameter (c) is caused by the external neutron source and by zero flux boundary condition.

Geometry of the Nuclear Slab System
The investigated geometry represents middle cut of one heterogenic VVER440 fuel assembly that is placed in the middle.One edge of the investigated geometry is represented by homogenous representation of the heterogeneous fuel assembly.Control rod assembly is placed on the opposite edge.This case was introduced due to simulation of the transient state during operational conditions.Geometry is also simplified into one dimension and some lengths are changed for faster convergence of numerical calculation (Fig. 2).Red material represents UO 2 (dimension d1 is 7 mm and is the same for every red and orange material).Orange material is the representation of UO 2 with 3.35 % of 153 Gd.
The two cases of the enrichment are investigated (Table 2).Green is B 4 C (dimension d5 is 142 mm),  blue stands for H 2 O (dimension d2 is 5 mm, d3 = 19 mm, d4 = 10 mm).Purple is a mixture of the whole area between x-y points (d5 is 142 mm).All cross sections are obtained from ENDF/B-VII.1 library for the temperatures of 296.3 K, 600 K and 1800 K, for the discrete energies of 0.0253 eV and 2 MeV [7].Isotropic scattering is assumed and the diffusion coefficient can be calculated according to (9) [8], where Σ t stands for the total macroscopic cross section, Σ tr is the transport macroscopic cross section, and μ0 represents unsymmetrical parameter of neutron scattering (if isotropic scattering is applied μ0 = 0).In this case, a molecular cross section of H 2 O is not considered, even if the molecular cross section is higher than calculated cross sections of particular nuclides [8].
More details about geometry input can be found in [4].It is necessary to note that this is the only theoretical approximation and that it is mainly because the values of the macroscopic scattering cross section were simplified with certain degree of uncertainty.

Initial and Boundary Conditions
The zero incoming current boundary condition is used in all cases of the neutron diffusion where J is the neutron current and φ is the neutron flux.This can be interpreted as that the environment behind the boundary of the nuclear system is vacuum.
According to (2), the stationary state is represented by proportional neutron flux distribution that is normed to maximal value 1.It is necessary to set the initial condition in the form of total power generation The absolute value of the neutron flux distribution is calculated from the proportional neutron flux distribution (2) and from the initial power generation condition (10).The power generation rate is calculated for the case 1 from the typical power generation of VVER reactor (1375 MW th ) and it is approximately 1.576 kW/m per the assembly where ASSMEBLY NR is the total number of fuel assemblies in the reactor and HEIGHT represents height of the core.The power generation rate is the same for the case 2 with higher enrichment to lower impact of the thermal feedback to ensure the nuclear system's subcritical state.Temperature is kept at constant 296.3K at the boundaries that correspond to Dirichlet boundary condition.Temperature of the whole system is also set to 296.3 K before the transient for all cases.No thermo-mechanical physics is applied (the densities of the particular materials are constant).Also no fluid mechanics is applied (the heat diffusion equation solves just thermal conduction equation) and this simplification results in total temperature increase of the whole nuclear system.

Case 1
k eff = 0.282 601 is calculated from the stationary state of TRNASOS code for this scenario.This condition may be achieved in spent fuel pool.The worst case is simulated by the neutron flux corresponding to the operational state before the transient (the power generation is 1.576 kW/m).It is possible to observe (Fig. 3) that the thermal energy spectrum of neutrons decreases rapidly due to deep subcritical state of the nuclear system.The same evolution can be seen in case of fast energy spectrum of neutrons in Fig. 4.
The temperature is set to 296.3 K at the beginning of the transient state.The power is generated by (10) according to the neutron flux and the temperature evolution is shown in Fig. 5.The temperature increases up to 460 K and according to the temperature change, the microscopic cross sections are modified.
The whole system is below melting point temperature of the fuel.

Case 2
The k eff is equal to 0.982 693 in this case, due to higher enrichment of 235 U (Table 2).The neutron source is incorporated within the core and its total emission of neutrons is 1.1 × 10 6 s −1 in the area between 217-227 mm on x-axis.It corresponds to the concept of super-thermal liquid-helium ( 4 He) source (UCN) [9] and portion of 1.0 × 10 6 s −1 is in the thermal energy spectrum.Subcritical multiplication can be derived from The result of subcritical multiplication formula is consistent with the calculated results that are shown in Fig. 6 and Fig. 7.The whole starts from zero power level.The transient state of the neutron flux is stabilized approximately after 50 000 s.However, the transient of the temperature evolution is represented by higher inertia and the transient state is stabilized in the region of 140 000 s (Fig. 8).The temperature does not exceed 380 K.

Conclusions
The successful application of the TRANSOS code was demonstrated.It has been found that the maximum value of the neutron flux and its temperature is located in the middle of the system within the heterogeneous region.On the other hand, the diffusion theory is more appropriate for the homogenous problems and calculation error increases with greater heterogeneity.In all cases, the temperature does not exceed melting point of fuel material.The code that was written in C++ can be used for academic purposes at the institute and for solution of simple nuclear systems.
It has been confirmed that the diffusion theory is suitable for fast calculations and it is used during development of theoretical reactor design.Therefore, development of simple codes increases the knowledge level of neutron physics and it can result in enhancing of nuclear safety.
The contribution of 238 U thermal feedback with thermo-hydraulic properties of the system has to be investigated in the future.

Figure 2 .
Figure 2. Graphical representation of investigated nuclear reactor core.

Figure 3 .
Figure 3. Neutron flux distribution for energy group 0.0253 eV in Case 1 (the neutron flux decreases with time evolution due to deep subcritical state).

Figure 4 .
Figure 4. Neutron flux distribution for energy group 2 MeV in Case 1 (the neutron flux decreases with time evolution due to deep subcritical state).

Figure 5 .
Figure 5. Temperature evolution in Case 1 with k eff = 0.282601 (the temperature decreases with time evolution -the nuclear system is brought into safeconditions).

Figure 6 .
Figure 6.Neutron flux distribution for energy group 0.0253 eV in Case 2 (nuclear system is stabilized after subcritical multiplication for k eff = 0.982 693).

Figure 7 .
Figure 7. Neutron flux distribution for energy group 2 MeV in Case 2 (nuclear system is stabilized after subcritical multiplication for k eff = 0.982 693).

Figure 8 .
Figure 8. Temperature evolution in Case 2 (temperature is stabilized below melting point of fuel material for k eff = 0.982 693).

Table 1 .
Deviation between numerical and analytical solution.