Finite Element Method with Spectral Green's Function in Slab Geometry for Neutron Diffusion in Multiplying Media and One Energy Group

ABSTRACT

The physical phenomenon of neutrons transport associated with eigenvalue problems appears in the criticality calculations of nuclear reactors and can be treated as a diffusion process. This paper presents a method to solve eigenvalue problems of neutron diffusion in slab geometry and one energy group. This formulation combines the Finite Element Method, considered an intermediate mesh method, with the Spectral Green's Function Method, which is free of truncation errors, and it is considered a coarse mesh method. The novelty of this formulation is to approach the spatial moments of the neutron flux distribution by the first-order polynomials obtained from the spectral analysis of diffusion equation. The approximations provided by the proposed formulation allow obtaining accurate results in coarse mesh calculations. To validate the proposed method, we compared its results with methods described in the literature. The accuracy and computational performance of our formulation were characterized by solving a benchmark problem with a high degree of heterogeneity.

Keywords:
Eigenvalue problems; Neutron diffusion equation; Spectral Green's function

RESUMO

A utilização de centrais nucleares para geração de energia elétrica é uma das principais alternativas para atender a demanda energética nos próximos anos sem agravar os problemas ambientais associados ao aquecimento global. O uso de métodos e técnicas de simulação que estimem a população de nêutrons no núcleo é fundamental para garantir a operação segura e confiável do reator. Neste trabalho apresentamos uma formulação híbrida para resolver problemas de autovalor de difusão de nêutrons em domínios unidimensionais e aproximação de uma velocidade. Esta formulação combina a simplicidade do Método deElementos Finitos com uma aproximação Espectro-Nodal e permite obter resultados precisos em cálculos de malha grossa. A precisão e o desempenho computacional do métodoproposto são caracterizados através da solução de um problema modelo com alto grau de heterogeneidade.

Palavras-chave:
Problemas de autovalor; Equação de difusão de nêutrons; Função espectral de Green

1 INTRODUCTION

Population growth and economic development require the continuous increase of electricity generation capacity. To attend this demand without aggravation of the environmental problems associated with the greenhouse effect is a major challenge for scientists, economists and politicians nowadays. In this context, nuclear energy appears as one of the generation alternatives without carbon dioxide emissions. In Brazil, the energy development plan anticipates, in the next years, an increase of nuclear power generation in the energy matrix¹³.

Among the main concerns related to the use of nuclear energy there is the safety and the radioactive materials contention in the reactor core. These issues are associated with the monitoring and control of the neutron population in the reactor. For this reason, during the design and operation of nuclear plants, we need accurate and efficient numerical methods to assess the neutron flux distribution and the effective multiplication factor. This problem is associated with the nuclear reactor criticality calculations, and is also known as a neutron transport eigenvalue problem¹⁶.

The physical modeling of the phenomenon of neutron transport involves the migration of neutrons in a host medium and the probability of interaction between neutrons and atoms⁵. Many scientists treat the neutron transport phenomenon as a diffusion process, thus, one of the mathematical models that allow us to describe this physical phenomenon is the neutron diffusion equation (DE). The DE is a simple model that provides accurate results for the distribution of neutron flux and the effective multiplication coefficient ¹². Many neutron transport problems are well modelled in one-dimensional Cartesian geometry. That is a good approximation when the neutron flux has rough changes in one spatial direction x (axial) and smoother changes in the symmetrical plane y - z. For that reason, many authors use one-dimensional Cartesian geometry models in recent neutron transport works^{9), (19), (20), (21}.

For computational performance reasons it is desirable to obtain accurate results for the distribution of neutron flux and the effective multiplication coefficient in a coarse mesh. Among the coarse mesh methods, we can highlight the Spectral Green's Function Nodal Methods (SGF)^{2), (10}. The SGF method had its origins in the Spectral-Nodal Method with constant approximation for the discrete ordinates model in source-fixed problems³. Since then, several other formulations have been developed, for example, the SGF method for discrete ordinates⁴ and the SGF formulation for diffusion approximation¹¹, both of them, in Cartesian geometry for neutron multiplying systems. The SGFs methods have high algebraic cost and complex computational algorithm and, for that reason, the linear Finite Element Method (FEM) has been used to solve some of the nuclear engineering problems. The FEM method offers appropriated results for fine and intermediate meshes^{22), (23}. The use of analytical solutions, enriching the finite element methods is well established in the solution of several engineering problems^{14), (15), (18}. However, this approach is less known in neutron transport applications.

The formulation proposed in this work combines the linear approximation of FEM method with a quasi-analytic approach of SGF to solve neutron diffusion eigenvalue problems in slab geometry and one energy group. The combination of these methods aims to generate accurate results with high computational performance, preserving the analytical spectral solutions inside the elements. We validate this formulation by solving benchmark problems and compare the obtained results with other methods.

In the next section we discuss the difussion mathematical model and the discretization scheme.In the section 3, we present the hybrid formulation fundamentals. The numerical results are offered in section 4 for a benchmark problem with high degree of heterogeneity. Finally, in section 5, we give the final comments and suggestions for future works.

2 THE MATHEMATICAL MODEL

In this paper, the deterministic mathematical modeling of the neutron transport in multiplying medium is based on the diffusion equation. We consider a static model in slab geometry, disregarding the energy dependence. This model^{16), (22} is a system composed by two equations, the balance equation

(2.1)

and, the Fick's law,

(2.2)

The equation (2.1) represents the mono-energetic diffusion equation, which describes the gain and loss of neutrons, while (2.2) refers to the Fick's law, which describes the neutron's difussion, that is, the neutrons migration from high density regions to low density regions¹². In (2.1) and (2.2), J(x) is the neutron current, Φ(x) is the neutron flux, Σ_α (x) is the macroscopic absorption cross section, k_eff is the effective multiplication coefficient, ν is the average number of neutrons released by fission event, Σ_f (x) is the macroscopic fission cross section, D(x) is the diffusion coefficient and L is the length of the one-dimensional domain.

The boundary conditions considered in this work are the Albedo boundary conditions ¹⁷, which are represented by

(2.3)

the terms α_LR are Albedo proportionality constants.

2.1 Spatial Discretization Scheme

In general, to obtain the numerical solution for differential equations systems, the domain discretization is required. At this point, we define a spatial mesh with total length L and I elements called Ω_i , i=1:I, each one with length h_i . The representation of this spatial mesh can be seen in Figure 1.

Figure 1
Spatial mesh in the domain of length L.

Now, we write the equations (2.1) and (2.2) in discretized form for each element Ω_i in the domain

(2.4)

(2.5)

then, the Albedo boundary conditions, equations (2.3), appear in the following discretized form

J_1/2 = αLΦ_1/2, J_{I + 1/2} = αRΦ_{I + 1/2.}

The next step is to obtain the spatial balance equations for the zero and first order moments in each element Ω_i . Thus, we need to apply the operator

in the equations (2.4) and (2.5). This operator is based in the Legendre polynomials¹, where P_l represents the polynomial of l degree. For the zero order balance equations, l=0, the Legendre polynomial is a constant, P ₀(x) =1. For first order, l=1, we have a linear polynomial, P ₁(x) =x.

The balance equations, for the zero order moment appear in the following form

where, the element average flux is

(2.6)

and the element average current is

For the first order balance equations we have the following expressions

where, the flux first order moment is

(2.7)

and the the neutron current first order moment is defined as

3 THE FEM-SGF FORMULATION

The hybrid formulation proposed (FEM-SGF) in this paper combines two methods found in the literature. One of these methods is the FEM that obtains the solution of the diffusion model for eigenvalue problems using approximation functions. In this work, we choose two linear approximation functions, one for the neutron flux and the other for neutron current, such functions satisfy the conditions described by integral equations in the domain. The other method that composes the hybrid formulation is the one-dimensional SGF², this method, reaches the solution of the diffusion equation with two spectral para meters α_i and β_i , obtained by spectral analysis of the mathematical model. Thus, the hybrid formulation proposed is based on the combination of linear approximation of the FEM with the approach quasi-analytic of the SGF.

The approximation functions of the FEM-SGF method appear in the following form

(3.1)

(3.2)

where, we introduce the spectral parameters of the SGF in the linear functions of the FEM.

To obtain the discretization of the approximation functions for the FEM-SGF, we need to evaluate those functions into the element edges. For the elements belonging to left half of the domain, we evaluate the functions (3.1) and (3.2) in the element right edge, x_{i + 1/2}, these functions discretized appear as

(3.3)

(3.4)

For the elements in the right half of domain, the functions (3.1) and (3.2) are evaluated in the element left edge and the result is

(3.5)

(3.6)

When the mesh has an odd number of elements, we must equate the average flux in equations (3.3) and (3.5), and the average current in equations (3.4) and (3.6), inside the central element to ensure the symmetry conditions in the domain. The resulting equations are

Now, we calculate the spectral parameters, α_i and β_i , of the approximation functions in the FEM-SGF. For that, we need to obtain general analytical solutions for the diffusion equation. These solutions can be calculated by the spectral analysis, a technique proposed by Case and Zweifel⁸. The general analytical solutions can be write as

where k ₁ and k ₂ are arbitrary constants, and μ_i are the eigenvalues of the problem. The arbitrary constants must be used to preserve the DE general analytical solution inside each spatialelement Ω_i .

Next, we can choose, arbitrarily, any of the functions in the general solution, because both functions carry the arbitrary constants k ₁ and k ₂, then, we select the function Φ(x), equation (3). The next step is to obtain the fluxes on the edges Φ_{i ± 1/2}, by evaluating (3) in the element edges (x_{i ± 1/2}), which appears as

(3.7)

at this point, we calculate the average flux 2.6), theresult appears as, by substituting (3) in the definition (

(3.8)

finally, we calculate the first order moment of neutron flux 2.7). The result appears in the following form, by substituting (3) in the definition (

(3.9)

Now, we can substitute the expressions for Φi ± 1/2, 3.7), (3.8) and (3.9), in (3.3) or (3.5) and write these in the form and , equations (

(3.10)

where,

But, the general analytical solutions must be preserved for any k1 and k2 values, thus, (3.10) will be satisfied if A=C and B=D. These equalities provide a linear system that lets us calculate the parameters αi and βi. The eigenvalue μi can be real or pure imaginary, thus, we can have two possible solutions,

(i) If

(3.11)

(ii) If

(3.12)

where λi =hi/.2μi The equations (3.11) and (3.12) used to calculate the parameters of auxiliary equations, ensure, that the general analytical solutions are preserved inside each spatialelement Ωi.

The algebraic linear equation system obtained by FEM-SGF is formed by (6I+2) unknowns and the same number of discretized equations. In this system we can eliminate the average element flux and the current (I+2) order. The reduced system, represents an eigenvalue problem, then, we can use direct or iterative techniques⁷ to obtain the neutron flux distribution and the power method⁷ to estimate the effective multiplication coefficient.). Thus, we obtain an equivalent and reduced equation system of (2 and and ). Also, we can remove the first order moment of the flux and the current (

4 RESULTS AND DISCUSSION

In order to quantify the accuracy and the computational performance of the developed method, we compare the proposed formulation with conventional methods, specifically, the Diamond Difference method (DD)¹⁷ and FEM method. DD is considered a fine mesh method and it has a low computational performance.

The benchmark chosen for the computational experiments has a domain with 150 cm of length, divided into 6 regions, whose materials and geometric parameters are shown in Table 1. The boundary conditions are reflexive in the left boundary and vacuum (α_R ⁼ 0.5) in the right boundary. The convergence criterion requires that the relative deviation between two estimates of the effective multiplication factor to be less than or equal to 10^-4%. The nominal power considered in the reactor is 10 MW.

Table 1
Materials and geometric parameters of benchmark problem.

The numerical results for neutron flux distribution appear in Figures 2 and 3. The reference solution was generated using a fine mesh composed of 960 elements with the FEM method. The Figure 2 shows the behavior of the neutron flux in the domain and the Figure 3 presents the relative deviations of the three methods with respect to reference solution. As we see, in this benchmark problem appear strong neutron flux gradients in some regions. Consequently, most of the numerical methods fail to obtain accurate results.

Figure 2
Neutron flux distribution in the domain for DD, FEM and FEM-SGF methods using 30 elements in the spatial grid.

Figure 3
Relative deviations in neutron flux distribution for DD, FEM and FEM-SGF methods in coarse mesh.

Figure 2 shows that the FEM-SGF method, with a coarse mesh (30 elements), keeps the neutron flux close to reference solution, while the FEM, with the same mesh, moves away from the reference solutions in the regions with strong gradients. On the other hand, the DD method generates meaningless results for a coarse mesh calculations.

The differences between numerical results obtained by the three methods can be more noticeable in Figure 3. The strongest variation of the neutron flux occurs in the region between x=50 and x=75. In this region the numerical methods find the major difficulty to estimate theneutron flux. The FEM-SGF presents the smallest relative deviations, 0.01% approximately, while in the FEM the relative errors vary between 1% to 100%. The relative errors associated to DD are even higher.

To analyze the reactor power density results, we choose the region between x=60 and x=100 cm, which presents a strong gradient in the neutron flux and a high neutron density.The Table 2 shows the results for the power density and the relative deviations for each method.

Table 2
Numerical results for power density in the region (60 ≤ x ≤ 100) for DD, FEM and FEM-SGF methods in coarse mesh.

Method (elements)	Power (MW)	Rel. deviation (%)
DD (30)	6.208733E-02	9.769013E+01
FEM (30)	2.523533E+00	6.115466E+00
FEM-SGF (30)	2.692115E+00	1.563635E-01
Reference: FEM (960)	2.687912E+00	...

The power density obtained using the FEM-SGF method presented a good agreement when we compare with the reference solution, its relative deviation is about 0.1%. However, the FEM method shows a relative deviation of 6%, approximately. The DD method was more sensitive to coarse mesh and the strong flux gradient in this region, presenting a relative deviation of 97.6%.

In the Figure 4, we present the relative deviation behavior, in power density calculations, as a function of the element number for all the discussed formulations. In general, we can see that the FEM-SGF reached the best results if compared with the conventional methods, for all spatial meshes considered. Specifically, when we compare the FEM-SGF with the FEM, we notice that it shows similar results with about 1/3 of the elements of the FEM. Therefore, the FEM-SGF is as accurate, using a coarse mesh, as the DD in fine mesh and the FEM in intermediate mesh.To complete this analysis, we compare the behavior of the relative deviation with the error analytical estimates for the analyzed methods (Fig. 4). The convergence order of the DD, FEM and FEM-SGF methods, in asymptotic notation is O(h ²)^{6), (17}. We can also see in Figure 4, that the curves of each method have similar behaviors modified by some multiplicative constant.

Figure 4
Relative deviation as a function of elements number for DD, FEM and FEM-SGF methods.

The numerical results for the effective multiplication coefficient (k_eff ) appear in Table 3. The k_eff value obtained by FEM-SGF have the smallest relative deviation, 10^-5%, approximately. The results for the conventional methods (FEM and DD), suffer a bigger influence of the spatial grid dimensions. They present relative errors of 10^-3% and 10^-2%, respectively.

Table 3
Numerical results for effective multiplication coefficient for DD, FEM and FEM-SGF methods in coarse mesh.

Method (elements)	keff	Rel. deviation (%)
DD (30)	1.140051E+00	3.218296E-02
FEM (30)	1.140309E+00	9.573804E-03
FEM-SGF (30)	1.140419E+00	7.634042E-05
Reference: FEM (960)	1.140418E+00	...

To characterize the FEM-SGF computational performance, we consider the dominant operation in the algorithmic implementation of the numerical formulations. This dominant operation is associated with the solution of the algebraic linear system generated by the discretization procedure in each method. The resulting matrices in the spatial discretization process are sparse, specifically, they are banded symmetric⁷. The bandwidth for DD, FEM and FEM-SGF methods are 3, 4 and 4, respectively. Thus, considering the dominant operation, the computational cost associated to the numerical formulations DD, FEM and FEM-SGF, in asymptotic notation, is O(M)O(I), where M is the iteration number of power method for eigenvalue convergence.It is noteworthy that the DD generate an algebraic linear system with (I+1) variables and equations, and the order of algebraic linear system generated by FEM is equivalent to FEM-SGF. Moreover, due to the specific characteristic of the FEM-SGF, the matrix coefficients depend on problem eigenvalue and must be updated on each iteration.

5 CONCLUSIONS AND FUTURE WORKS

We presented in this work a hybrid formulation to solve the one-dimensional neutron diffusion equation in multiplying media and one energy group. The hybrid character of this formulation is based on the combination of the FEM linear approximations with the SGF spectral parameters. These parameters allow us to preserve the general analytic solutions inside each element. The FEM-SGF formulation, permits to obtain accurate numerical results in coarse mesh, suggesting a high computational performance.

The FEM-SGF generates satisfactory numerical results for the neutron flux distribution and the effective multiplication coefficient in a coarse spatial grid. Considering accuracy and computational performance, the proposed method presents better results than conventional methods, FEM and DD, in most of the cases. The computational experiments with coarse meshes showed that, in regions with strong neutron flux gradients, the FEM-SGF keeps the numerical results close to the reference solution unlike FEM and DD methods that generate results far from it. The hybrid formulation using a coarse mesh presents the same precision than the compared methods with finer meshes and, for that reason, it has the better computational performance.

As future works, we suggest to reduce the FEM-SGF algebraic linear system order to (I+1).We also suggest to extend the hybrid formulation to solve multi-dimensional problems, andconsider multi-group approximations.

ACKNOWLEDGMENTS

Authors acknowledge to Conselho Nacional de Desenvolvimento Científico e Tecnológico(CNPq, BRAZIL) and Fundação de Amparo à Pesquisa do Estado da Bahia (FAPESB,BRAZIL) for the financial support provided. Also, the authors acknowledge the Universidade Estadual de Santa Cruz, for supporting this research and the Núcleo de Biologia Computacional e Gestão de Informações Biotecnológicas for offering the infrastructure for the development of the computational experiments.

REFERENCES

• ¹
  R.E. Attar. "Legendre Polynomials and Functions". CreateSpace, (2009).
• ²
  R.C. Barros & E.W. Larsen. A Numerical Method for Multigroup Slab-Geometry Discrete Ordinates Problems with No Spatial Truncation Error. Transport Theory and Statistical Physics, 20 (1991), 441-462.
• ³
  R.C. Barros & E.W. Larsen. A spectral nodal method for one-group X,Y-geometry discrete ordinates problems. Nuclear Science and Engineering, 111(1) (1992), 34-45.
• ⁴
  R.C. Barros, H.A. Filho, E.T.V. Orellana, F.C. Silva, N. Couto, D.S. Dominguez & C.R.G. Hernandez. The Application of Spectral Nodal Methods to Discrete Ordinates and Diffusion Problems in Cartesian Geometry for Neutron Multiplying Systems. Progress in Nuclear Energy, 42 (2003),385-426.
• ⁵
  G.I. Bell & S. Glasstone. "Nuclear reactor theory". Van Nostrand Reinhold Co., (1970).
• ⁶
  S.C. Brenner & L.R. Scott. "The Mathematical Theory of Finite Element Methods". Springer-Verlag New York, Inc., (1996).
• ⁷
  R.L. Burden & D.J. Faires. "Numerical Analysis". Ninth Edition, BROOKS/COLE - CENGAGE Learning, (2011).
• ⁸
  K.M. Case & P.F. Zweifel. "Linear transport theory". Addison-Wesley Pub. Co., (1967).
• ⁹
  C. Ceolin, M. Schramm, B.E.J. Bodmann, M.T. Vilhena & S.B. Leite. On an analytical evaluation of the flux and dominant eigenvalue problem for the steady state multi-group multi-layer neutron diffusion equation. Kerntechnik, 79 (2014), 430-435.
• ¹⁰
  D.S. Dominguez & R.C. Barros. The spectral Green's function linear-nodal method for one-speed X,Y-geometry discrete ordinates deep penetration problems. Annals of Nuclear Energy, 34 (2007), 958.
• ¹¹
  D.S. Dominguez, C.R.G. Hernandez & R.C. Barros. Spectral nodal method for numerically solving two-energy group X,Y geometry neutron diffusion eigenvalue problems. International Journal of Nuclear Energy, Science and Technology (Print), 5 (2010), 66.
• ¹²
  J.J. Duderstadt & L.J. Hamilton. "Nuclear Reactor Analysis". John Wiley & Sons Inc, (1975).
• ¹³
  Empresa de Pesquisa Energética. "Plano Nacional de Energia 2030", Ministério de Minas e Energia, Rio de Janeiro, Brasil, (2007).
• ¹⁴
  T. Hayashi & D. Inoue. Calculation of leaky Lamb waves with a semi-analytical finite element method. Ultrasonics, 54(6) (2014), 1460-1469.
• ¹⁵
  J.D. Jung & W. Becker. Semi-analytical modeling of composite beams using the scaled boundary finite element method. Composite Structures, 137 (2016), 121-129.
• ¹⁶
  J.R. Lamarsh & A.J. Baratta. "Introduction to Nuclear Engineering". Prentice Hall, (2001).
• ¹⁷
  E.E. Lewis & W.F. Miller Jr. "Computational Methods of Neutron Transport". American Nuclear Society, Illinois, USA, (1993).
• ¹⁸
  P. Liu, D. Wang & M. Oeser. Application of semi-analytical finite element method coupled with infinite element for analysis of asphalt pavement structural response. Journal of Traffic and Transportation Engineering (English Edition), 2(1) (2015), 48-58.
• ¹⁹
  E. Sauter, F.S. Azevedo, M. Thompson & M.T. Vilhena. Eigenvalues of the Anisotropic Transport Equation in a Slab. Transport Theory and Statistical Physics, (2012), 448-472.
• ²⁰
  E. Sauter, F.S. Azevedo, M. Thompson & M.T.M.B. Vilhena. Solution of the one-dimensional transport equation by the vector Green function method: Error bounds and simulation. Applied Mathematics and Computation, 219, (2013), 11291-11301.
• ²¹
  A.C. da Silva, A.S. Martinez & A. da C. Gonçalves. Reconstruction of the Flux in a Slab Reactor. World Journal of Nuclear Science and Technology, (2012), 181-186.
• ²²
  W.M. Stacey. "Nuclear Reactor Physics", Wiley-VCH, (2007).
• ²³
  O.C. Zienkiewicz. "The Finite Element Methods in Engineering Science". McGraw-Hill, (1971).

Publication Dates

• Publication in this collection
  May-Aug 2016

History

• Received
  18 Nov 2015
• Accepted
  01 Aug 2016