Open Journal of Applied Sciences
Vol.07 No.02(2017), Article ID:74518,14 pages

Estimates of the Fast and Termal Flux in Blanket of Critical Reactors by Using Multi-Group Methods

Aybaba Hançerlioğullari1*, Aslı Kurnaz1, Yosef G. Ali Madee1, Ltfei A. Abdalsmd2, Salem A. A. Shufat2, Khaled M. Elhadad2, Hand Hadia Almezogi2, Mansur Mohamed Ali Mansur2

1Faculty of Arts and Sciences, Kastamonu University, Kastamonu, Turkey

2Engineering Faculty, Kastamonu University, Kastamonu, Turkey

Copyright © 2017 by authors and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

Received: June 28, 2016; Accepted: February 25, 2017; Published: February 28, 2017


In this study, based differential equations methods are used to solve equations because these methods are dependent on boundary value data more than other mathematical equations. We have calculated neutron flux, criticality and geometrical eigenvalue by using multi-group method and solving the neutron diffusion equation for finite and infinite cylindrical and spherical reactors in this study. For the calculation of the total neutron flux cross sections, we need the neutron diffusion equation. Thus, we have established the relationship be- tween neuron flow and cross-section of neuron depending on neutron energy. Critical calculations have been made by comparing the results with MNCP (montecarlo n-partical) simulation methods. For necessary computer calculations, the programme, Wolfram-Matematica-7 has been used.


Critical Reactor, Neutron Diffusion Equation, Mcnp, Multi-Group Method, Simulation

1. Introduction

Bessel differential equations are used for the calculations of neutron flux () and criticality coefficient (K) and cylindrical geometric structure is taken into account as the reactor geometry. Nowadays, the procedure of obtaining electric energy from nuclear energy is supplied by light water or heavy water reactors. Precious fossil fuel energy sources which can be divided with thermal neutrons are used in these reactors. Thermal neutrons are important for the fission reactions. Thermal and delayed neutrons are important for the continuity of the nuclear reaction. During the occurring of these important reactions, the thermal power and the performance of the materials’ structural reactor change significantly.

The Bessel Equation is formulated as follows and this equation is a special case of Bessel’s Equation [1] :


In which n is an integer (n = 0, 1, 2, 3…), if we let and, after multiplication by r. Using these approaches, we can reach the balance equation. From our recent discussions, we recognize this as Bessel’s differential Equation. The geometric eigenvalue, shown by in this equation, is

given as. Leakages, taken into account also in the diffusion equa-

tion don’t allow neutron flux to be zero on the border [2] :


The solutions of this equation are called Bessel Functions of order n. Since Bessel’s differential equation is a second order ordinary differential equation, two sets of functions, the Bessel function of the first kind Y1=A Jn(x)0 and Y2 = CYn(x) are the solutions to the above formulated equation: [1] - [6] .


Y1 and Y2 are respectively called as the functions of the Bessel function of the first kind and the Bessel function of the second kind [7] - [12] .

2. Method

2.1. Reactor Geometry for Neutron Flux

We have calculated neutron flux, criticality and geometrical eigenvalue by using multi-group method and solving the neutron diffusion equation for finite and infinite cylindrical and spherical reactors in this study. The neutron distribution in the reactor can be explained in these ways. At any time, we can specify the neutron angular distribution as Ω by connected to a specified E energy, the number of neutrons in a unit volume [4] [5] [6] . We can write the expanded neutron diffusion equation in homogeneous reactors as. In the diffusion equation, D is used for diffusion coefficient and L stands for diffusion length and is used for the calculation of diffusion length with the help of where, a is absorbtion macroscopic cross section [13] [14] [15] [16] .


If p = 0, the reactor is plane; if p = 1, the reactor is cylinder and if p = 2, it can be thought that the core of the reactor is spherical. In case of critical reactor, considering that the number of neutrons will not be changed by the time, the Boltzmann diffusion equation turns to:


In the equation, the first term stands for leakage neutrons, the second term is for the absorption neutron and the last term S is for the neutron resources in the reactor core. The expression gives the number of absorption in reactor core and the expression is in a unit time and volume.

2.2. The Application of Modified Bessel for NDM

The Bessel function of the second kind of order can be expressed in terms of the Bessel function of the second kind also known as the Weber Function. As it is shown in Figure 1 by arrows, neutron flux must be finite everywhere the diffusion equation is applied. At Figure 1 is given for a reactor with an infinite R radius, the flux change is maximum when R = 0. The size of a reactor with the reflector installed can be much smaller than of a reactor with the same material but without the reflector. For critical reactor P is possibility of not leaking and expressed as following formula [5] [6] (Figure 1).


2.3. The Flux Distribution in the Bare Infinite Cylinder

We can write Boltzmann equation for an infinite cylindrical core with the help of modified Bessel equation. In this equation, we can write diffusion Equation separately from an infinite cylindrical reactor core with a bare R radius .When we apply the expression Laplace in cylindrical coordination to diffusion equation,


it becomes:


For the bare infinite cylindrical reactor, geometrical factor (buckling factor)

can be written as B2 = ()2. The total power of the reactor is written as


Figure 1. The model in fast and termal fluxes in reactor reflector andcore.

2.4. The Flux Distribution in the Bare Finite Cylinder

For a finite cylindrical core with H height and R radius, the diffusion Equation becomes in geometrical eigenvalue:

[2] [4] [6] (9)

Considering that the flux is zero in the expanded radius of the core, the function, Y2(Br), must be C = 0 because it is infinite when r = 0.

In this case the equation, is the essential solution.

is written as


2.5. The Flux Distribution in a Bare Spherical Reactor

If we express the flux distribution for the spherical reactor that we study on, the operator, , can be expressed as it is written below. Hence, a diffusion equation for a spherical reactor is obtained as it’s written below with the help of the expression.

Changing a parameter in the solution of diffusion equation can be written

with the method of instead of [2] .


and written as ,by using the parameter changing method of diffusion equation and changing and by converting into a new Bessel function for a finite

cylindrical reactor core the solution becomes like that.. Solution

of this differential equation,


Considering that the flux is never infinite in anywhere in the reactor, it must be = 0. Therefore, when it is = A, the last version of the flux distribution expression becomes and geometrical eigenvalue

factor becomes B2 =.

3. The Flux Distribution in HAM and MCNP

Homotopy analysis method, based on the homotopy term which is one of the basements of topology and differential geometry was stated by Shi Jun Liao in 1992 and since that, it has been applied in various areas of economics and engineering [17] [18] [19] . HAM is an analytical method providing mathematical, linear, nonlinear, partial differential and differential-integral solutions for the equations. This method can be applied in the boundary conditions of the problem handled in neutron diffusion Equations, this method includes an equation system like j = 1, 2, 3. If an equation when the space boundary conditions of neutron diffusion equation are valid is written, the HAM equation is written as and becomes as N is operator

Considering that the N linear is the operator for the bare finite cylinder with HAM method, the Laplace expression changes as it is written below [20] [21] [22] .


Considering, as a result of the basic HPM operator’s expansion in a series, it becomes


It forms an equation set in the upper series of p. when x=0, the expression is equal to A and finite.

Therefore, the solution of homotopy equation in terms of power series;


And the expression in limiting case;

or (16)

for x, diffusion flux [18] [19] [20] [21] [22] .

This flux must be given in boundary conditions. The Laplace expression for a finite cylindrical flux with the homotopy method, Considering


the solution of diffusion flux distribution equation will be:


with the method of Homotopy. Considering the partial differential solution, the expressions, respectively,


As a result, the finite cylindrical flux distribution function can be:

[18] [19] [20] [21] [22] . (20)

MCNP simulation is randomly number selection technique from one or more probabilistic distribution in a special trial or simulation study. The complexity in the nature of the industrial problems unfortunately makes analytical solution impossible diffusion problems. Monte-Carlo simulation method (MCNP) is randomly number selection technique from one or more probabilistic distribution in a special trial or simulation study. The method was then adopted easily for solution of much more complicated and non-statistical problems such as Integra-differential evaluation problems [19] [20] . It is possible to calculate the multiple integrals on phase transitions by Monte Carlo method [8] . When the integral of an f(x) function between [a, b].It becomes:


In this case, its integral is calculated by multiplication of average value with (b-a). If the arithmetical average of the function on N points, chosen arbitrarily between [a, b] is calculated, it becomes


Therefore, Monte Carlo reaches the integration:


We can benefit from the analyses described in the programme for flux distribution with MCNP method. On the geometrical surface and cells of the reactor , we study on, the code requirement is stated with F1, F2, F4, F5, F6, F7, F8, *F1, *F2, *F3, *F4, *F5 analysis [13] [14] [15] .

The above seven tally categories represent the basic MCNP tally types. To have many tallies of a given type, add multiples of 10 to the tally number. For example, F1, F11, F21, ・・・ , F981, F991 are all type F1 tallies

The quantities actually scored in MCNP before the final normalization per starting particle are presented in Table 1. Note that adding an (*Fn) changes the units and multiples the tally as indicated in the last column of Table 1. The F1 surface current tally estimates the following quantity:


Table 1. Flux distribution with MCNP method.

This tally is the number of particles (quantity of energy for *F1 crossing a surface. The scalar current is related to the flux as

The range of integration over area, energy, time, aand angle can be controlled by FS, E, T, and C cards, respectively. The FT card can be used to change the vector relative to which is calculated (FRV option) or to segregate electron current tallies by charge

The F2, F4 and F5 flux tallies are estimates:


4. Track Lenth Estimate of Cell Flux (F4)

The definition of particle flux is, where = particle velocity, N = particle density = particle weight/unit volume. Roughly speaking, the time integrated flux:


More precisely, let ds = vdt. Then the time integrated flux:


Beacuse is a track length density, MCNP estimates this integral by summing for all particle tracks in the cell, time range and energy range. Because of the track length term in the numerator, this tally is known as a track length estimate of the flux. It is generally quite reliable because there are frequently many tracks in a cell (compared to the number of collisions), leading to many contributions to this tally.

5. Surface Flux (F2)

The surface flux is a surface estimator but can be thought of as the limiting case of the cell flux or track length estimator when the cell becomes [10] [11] .


As the cell thickness approaches zero, the volume approaches and the track length approaches, where, the angle between the surface normal and the particle trajectory. This definition of flux also follows directly from the relation between flux and


MCNP sets..

The F2 tally is essential for stochastic calculation of surface areas when the normal analytic procedure fails. Surface areas when the normal analytic procedure fails [20] [21] .

6. Calculations

We have calculated a specific for the flux distribution expression written below that we obtained for spherical reactors, which stands for the r val-

ues for each different expansion values of r =. We have obtained the re-

sults in Table 2 by calculating neutron flux manually for different r values in this formula. As we stated before, (expansion radius) is the distance where neutron flux is zero out of the reactor. We can see that in case of, is zero. is the flux in the reactor center, namely is the maximum flux when r = 0.

In Table 3 microscopic influence lines of some of the materials which are used as fuel in fission reactor has shown.

Table 2. The flux distribution in a spherical reactor (Re = 5 m).

Table 3. Thermal neutron microscopic influence lines of the fuel materials [6] .

In Table 2 shows the flux distribution depending on the distance and theneutronic data of flux distribution for a homogeneous spherical core in case of R.

Connections which are obtained from the result of neutron diffusion equations and their boundary-value conditions obtained by using modify Bessel equations are shown in Table 4 in detail. In Table 4 maximum flux ratio on average flux, geometric buckling (B2) and flux distribution equations of some reactor geometries are shown. Similarly the datum given for the expansion radius for the same spherical reactor in Table 5 and the flux distribution is given in Figure 2 and Figure 3 after obtained by using Matematica-7 programme. Similarly, in Table 6 the flux distribution in a finite cylindrical reactor for R = 2.405 and H = 4.81 is given in Figure 3.

Table 4. NEM and MCNP calculated flux values for the infinite bare cylinder (R = 2.405 m).

Table 5. The flux distribution in a spherical reactor for Re = 2 m.

Figure 2. The flux distribution in a spherical reactor for Re = 2 m.

Figure 3. The flux in infinite cylindrical reactor for R = 2 m and R = 2.405 m (H = 4.81 m).

Table 6. Flux distributions in a finite bare and in an infinite bare reactor.

In Figure 3 is show that flux distribution in a spherical the reactor for radius R = 2 m. It appears to be greater in a small radius of the neutron flux. In Figure 3 and Figure 4, core flux distribution are shown for the radius R = 2 and R = 2.405 in finite and the infinite cylindrical. By using a modified Bessel function for infinite cylindrical reactor R = 2 m of the radius 0, 1, 2 state of flux distributions is shown in Figure 5. In both of the graphics, a view where neutron flux decreases from the center of reactor till the expansion radius is drawn and takes the value 0. After this point, it is seen that it takes negative values. One of the boundary conditions for reactors shows that the neutron flux can never be zero where the diffusion equation is applied. Actually we cannot say that these negative results are in contradiction with this boundary condition.

Because we applied the diffusion equation till the point from r = 0 to r = Re. We should analyze the results we found in spherical reactors for neutron flux in order to be able to compare. In the solution of diffusion equation in cylindrical reactors, a different way is used than we used for the solution of diffusion equation in spherical reactors (see in Table 7).

Figure 4. The flux distribution in a finite cylindrical reactor for R = 2 m.

Figure 5. The flux distribution in case of n = 0, 1, 2 for an infinite cylindrical reactor.

Table 7. Flux distributions finite and infinite height 2H flux distributions for cylinder radius R finite.

In Table 2 we can see the geometric eigenvalue flux distribution expressions found for finite cylindrical and infinite cylindrical reactors. We have chosen the flux distribution obtained for a finite cylindrical reactor to compare. Because the spherical reactor that we study on is handled as a finite system. As it is seen from the graphic, there is no sinusoidal change in the flux distribution obtained for cylindrical reactors. We can say that it draws a picture where the flux decreases till a specific distance from the maximum flux and after reaching the zero it increases. We can also say that the flux in spherical flux takes negative values as much as the twice of the expansion radius after decreasing till the expansion radius from maximum flux.

7. Conclusion and Discussion

Nuclear reactors are the complex machine-equipment systems constructed through the use of advanced engineering technologies. Fission-type reactors are devices developed to generate energy at a stable power by taking the chain reaction under control [15] . Bessel differential equations are second order ordinary differential equations and they offer solutions in the cylindrical, spherical and polar coordinates easily and also required physical parameters in the reactor can easily be obtained through the use of Bessel differential equations [8] . Neutron flux () in the reactor changes according to the geometry of the reactor. Three of the geometries in the field have a smallest critical size and the cube is the largest [8] . Spherical geometry is very difficult to build. For this reason, most of the reactors are constructed in the shape of a cylinder. The size of the reactor must be greater than the minimum critical size to allow for combustion of fuel and accumulation of fission products that are absorbing neutrons for an infinite reactor which can be defined k∞ = ε η fp, as well as, define keff for a finite reactor with leakage terms. At this here, keff; reactor critical coefficient, B; buckling coefficient, Ls; fast neutron moderation length, L; diffusion length and Pf ; reactor shape factors and Pt; reactor size factors. In the cylindrical reactor case, for height and radius, the following is given [15] .

Cite this paper

Hançerlioğullari, A., Kurnaz, A., Madee, Y.G.A., Abdalsmd, L.A., Shufat, S.A.A., Elhadad, K.M., Almezogi, H.H. and Mansur, M.M.A. (2017) Estimates of the Fast and Termal Flux in Blanket of Critical Reactors by Using Multi-Group Methods. Open Journal of Applied Sciences, 7, 68-81.


  1. 1. Gray, A. (1922) A Treatise on Bessel Functions and Their Applications to Physics. Macmillan, London.

  2. 2. White, J.R. (1998) Mathematical Methods the Strum Lioville Problem Neutron Diffusion in Nuclear Lecture Notes Mass.

  3. 3. Watson, G.N. (1922) A Treatise on the Theory of Bessel Functions. Cambridge University Press, Cambridge.

  4. 4. Bataineh, A.S., Noorani, M.S.M. and Hashim, I. (2009) On a New Reliable Modification of Hootopy Analysis Method. Communications in Nonlinear Science and Numerical Simulation, 14, 409-423.

  5. 5. Lamarsh, J.R. (1983) Introduction to Nuclear Engineering. 2nd Edition, Addison-Wesley, Reading, 102-274.

  6. 6. Aybers, N. (1990) Nuclear Reactor Engineering—I. 65-117.

  7. 7. Byerly, W.E. (1895) Fourier’s Series and Spherical. Ginn & Company, Boston.

  8. 8. Kuzmin, R.O. (1930) Bessel Functions.

  9. 9. Lebedev, N.N. (1972) Special Functions and Their Applications. Physico-Technical Institute Academy of Sciences, USSR.

  10. 10. Horie, J. (2012) Numerical Solution to Critical Problem of Finite Cylindrical Reactor by Variation Method. Journal of Nuclear Science and Technology, 11, 359-368.

  11. 11. Rsic Computer Code Collection Mcnp 4B (1995) London.

  12. 12. Horie, J. and Nishihara, H. (1975) Reduced Two-Dimensional Critical Problem of Finite Cylindrical Reactor. Journal of Nuclear Science and Technology, 12, 3531-3542.

  13. 13. Thomas, W.J. (2007) Reactor Physics Simulation with Coupled Monte Carlo Calculation and Computational Fluid Dynamics. International Conference on Emerging Nuclear Energy Systems, Istanbul, 3-8 June 2007, 36-46.

  14. 14. Shayesteh, M. and Hahriari, M.S. (2009) Calculation of Time-Dependent Neutronic Parameters Using Montecarlo Method. Annals of Nuclear Energy, 36, 901-909.

  15. 15. Hançerliogullari, A. (2012) Neutronic Calculations at Uranium Powered Cylindrical Reactor by Using Bessel Differential Equation. Nuclear Science and Technology, 15-24.

  16. 16. Sekimoto, H. (2007) Nuclear Reactor Theory. COE-INES Tokyo Institute of Technology, Tokyo.

  17. 17. Çavdar, S. (2010) Application of the Adomian Decomposition Method to One Group Neutron Diffusion Equation. 5th International Ege Energy Symposium and Exhibitition, Denizli, 27-30 June 2010,1021-1022.

  18. 18. Cassell, J.S. and Williams, M.M.R. (2004) A Solution of Neutron Diffusion Equation for a Hemispherical with Mixed Boundary Conditions. Annals of Nuclear Energy, 31, 1987-2004.

  19. 19. Khasawneh, K., Dababneh, S. and Odibat, Z. (2009) A Solution of the Neutron Diffusion Equation in Hemispherical Symmetry Using the Homotopy Perturbation Method. Annals of Nuclear Energy, 36, 1711-1717.

  20. 20. Bremister, J. (1993) Mcnp-4a General Monte Carlo Code N-Particle Transport Co-de. Version 4a, La-12625, Los Alamos National Laboratory, New Mexico.

  21. 21. Johston, R. (1963) A General Monte Carlo Neutronics Code. LAMS-2856, Los Alamos.

  22. 22. Dababneh, S., Khasawneh, K. and Odibat, Z. (2011) An Alternative Solution of the Neutron Diffusion Equation in Cylindrical Symmetry. Annals of Nuclear Energy, 38, 1140-1143.