Numerical modeling of neutron transport in SP3 approximation by finite element method
Alexander V. Avvakumov Valery F. Strizhov Petr N. Vabishchevich Alexander O. Vasilev






ki]Alexander V. Avvakumov
Avvakumov2009@rambler.ru

nsi]Valery F. Strizhov
vfs@ibrae.ac.ru

nsi,univ]Petr N. Vabishchevichcor
vabishchevich@gmail.com

univ]Alexander O. Vasilev
haska87@gmail.com

[ki]National Research Center Kurchatov Institute,  1, Sq. Academician Kurchatov, Moscow, Russia
[nsi]Nuclear Safety Institute, Russian Academy of Sciences, 52, B. Tulskaya, Moscow, Russia
[univ]North-Eastern Federal University, 58, Belinskogo, Yakutsk, Russia

[cor]Corresponding author


The SP_3 approximation of the neutron transport equation allows improving the accuracy for both static and transient simulations for reactor core analysis compared with the neutron diffusion theory. 
Besides, the SP_3 calculation costs are much less than higher order transport methods (S_N or P_N). 
Another advantage of the SP_3 approximation is a similar structure of equations that is used in the diffusion method. 
Therefore, there is no difficulty to implement the SP_3 solution option to the multi-group neutron diffusion codes. 
In this work, the application of the SP_3 methodology based on solution of the λ- and α-spectral problems has been tested for the IAEA-2D and HWR reactor benchmark tests. 
The FEM is chosen to achieve the 3D geometrical generality, using GMSH as a generic mesh generator. 
The results calculated with the diffusion and SP_3 methods are compared with the reference transport calculation results. 
It was found for the HWR reactor test that some eigenvalues are complex when calculating using both diffusion and SP_3 options.



neutron transport equation, diffusion theory, SP_3 approximation, reactor core, spectral problems, eigenvalues.







§ INTRODUCTION

The diffusion approximation of the neutron transport equation is widely used in nuclear reactor analysis allowing whole-core calculations with reasonable accuracy. 
The main feature of neutron diffusion equation is following: it is assumed that the neutron current is proportional to the neutron flux gradient (the Fick’s law). 
There are also three assumptions: neutron absorption much less likely than scattering, linear spatial variation of the neutron distribution and isotropic scattering <cit.>. 
To provide the validity of diffusion theory, the modern diffusion codes use, as a rule, assembly-by-assembly coarse-mesh calculation scheme with effective homogenized cross sections, prepared by more accurate transport approximations. 
To improve the diffusion code restrictions related with limitations on mesh spacing, different approaches are used including nodal and finite element methods <cit.>.

For many situations of interest (for instance, the pin-by-pin calculation taking account strongly absorbing control rods), the applicability of neutron diffusion theory are limited. 
Therefore, a more rigorous approximation for the neutron transport is required.

The solution of the neutron transport equation is very complicated problem because of seven independent variables: five for space-angular description, one for energy and one for time. 
To simplify the transport problem, different approaches are used such as the spherical harmonics (P_N) approximation <cit.>. 
The P_N approximation of the neutron transport equation is derived by expansion of the angular dependence of the neutron flux in the N spherical harmonics. 
During the last time, the simplest version of the P_N method, namely the simplified P_N approximation became widespread <cit.>. 
The major feature of the SP_N method is following: the three-dimensional neutron transport equation is transformed to a set of one-dimensional equations. 
The number of the SP_N trial functions is equal to 2(N+1) compared with the P_N method which uses (N+1)^2 trial functions. 
This leads to significant reduce in the computation time for typical whole-core calculations. 

The SP_N approximation was first derived by Gelbard <cit.> in the early 1960s. 
He replaced the spatial derivatives with Laplacian and divergence operators in a one-dimensional planar geometry. 
The resulting SP_N equations are elliptic, for example, the SP_3 equations consist of two equations of diffusion type with two unknown fluxes: the scalar flux and the second angular flux moment. 
More rigorous theoretical foundation of the SP_3 methodology has been derived by Brantley and Larsen <cit.> on the basis of variational methods.

The SP_3 method, as expected, can provide accuracy improvement compared with the common used diffusion method. 
Besides, implementation of the SP_3 equations into the diffusion code is not difficult because of the similar structure of the SP_3 and diffusion equations. 
For this reason the SP_3 method was adopted in different whole-core calculation codes, such as DYN3D <cit.>, PARCS <cit.> and others. According to <cit.>, application of the SP_3 theory to the pin-by-pin calculation for BWR geometry resulted in remarkable improvement in the calculation accuracy compared with the diffusion method. 
Another report <cit.> shows the comparison of the diffusion and  SP_3 methods to calculate the control rod reactivity  in a light-water reactor. As compared with the Monte Carlo reference calculation, the  SP_3 method gives twice as accurate result compared with the diffusion method. Besides, as it turned out, the computation time using the SP_3 method is only 1.5 times longer than that using the diffusion method <cit.>.

Thus, the SP_3 method can be considered as an improved approximation of the neutron transport equation compared with the diffusion method. 
In this regard, it will be very useful to compare the spectral parameters, calculated by both the diffusion and SP_3 methods. 
To characterize the reactor steady-state conditions or dynamic behavior, some spectral problems are considered <cit.>. 
The steady-state condition is usually described by solution of a spectral problem (λ-eigenvalue problem); the fundamental eigenvalue (the largest eigenvalue) is called k-effective of the reactor core <cit.>. 
The reactor dynamic behavior can naturally be described on the basis of the approximate solution expansion in time-eigenvalue of α-eigenvalue problem <cit.>. 
At large times, one can talk about the asymptotic behavior of a neutron flux, whose amplitude is exp(αt). 
Because the operator matrix is nonsymmetric, some of its
eigenvalues can be complex <cit.>.
Previously the complex eigenvalues and eigenfunctions were found in the spectral problems for some numerical tests <cit.>.

In this paper we consider the SP_3 approximation for the steady-state multi-group neutron transport problem.
To solve spectral problems with nonsymmetrical matrices we use well-designed algorithms and relevant free software including the library SLEPc (Scalable Library for Eigenvalue Problem Computations, http://slepc.upv.es/). 
We use a Krylov-Schur algorithm, a variation of Arnoldi method, described in <cit.>.

The paper is organized as follows. 
The steady-state and dynamic models of a nuclear reactor based on the multigroup SP_3 equations are given in Section 2. 
In Section 3 we discuss various spectral problems. 
Some numerical examples of calculation of spectral characteristics of two-dimensional test problems for IAEA-2D benchmark problem and HWR reactor using the two-group system of diffusion and SP_3 equations is discussed in Section 4. 
The results of the work are summarized in Section 5.



§ PROBLEM STATEMENT


Let’s consider the symmetric form of the SP_3 equation for the neutron flux <cit.>.
The neutron dynamics is considered in the limited convex two-dimensional or three-dimensional area  Ω (x = {x_1, ..., x_d}∈Ω,  d = 2,3) with boundary ∂Ω. 
The neutron transport is described by the system of equations

    1/v_g∂ϕ_0,g/∂ t - 2/v_g∂ϕ_2,g/∂ t    -  ∇· D_0,g∇ϕ_0,g + Σ_r,gϕ_0,g -  2Σ_r,gϕ_2,g
        =  (1-β)χ_n,g S_n + S_s,g + χ_d,g S_d, 
    
     -2/v_g∂ϕ_0,g/∂ t + 9/v_g∂ϕ_2,g/∂ t    - ∇· D_2,g∇ϕ_2,g + (5Σ_t,g + 4Σ_r,g) ϕ_2,g - 2Σ_r,gϕ_0,g
        =  -2(1-β)χ_n,g S_n - 2S_s,g - 2χ_d,g S_d,

where

    S_n =  ∑_g'=1^GνΣ_f,g'ϕ_g', 
      
    S_s,g = ∑_g≠ g'=1^GΣ_s,g'→ gϕ_g',
      
    S_d = ∑_m=1^Mλ_m c_m,


    ϕ_0,g=ϕ_g + 2ϕ_2,g, 
      
    D_0,g = 13Σ_tr,g, 
      
    D_2,g = 97Σ_t,g, 
       g=1,2,...,G.

Here G — number of energy groups,
ϕ_g( x, t) — scalar flux,
ϕ_0,g( x, t) — pseudo 0th moment of angular flux,
ϕ_2,g( x, t) — second moment of angular flux,
Σ_t,g( x, t) — total cross-section, 
Σ_tr,g( x, t) — transport cross-section, 
Σ_r,g( x, t) — removal cross-section,
Σ_s,g'→ g( x, t) — scattering cross-section,
χ_g  — spectra of neutrons, 
νΣ_f,g( x, t) — generation cross-section,
c_m( x, t) — density of sources of delayed neutrons,
λ_m — decay constant of sources of delayed neutrons,
M — number of types of delayed neutrons.

The density of sources of delayed neutrons is described by the equations

    ∂ c_m/∂ t + λ_m c_m = β_m S_n,
        m = 1,2, ..., M,

where β_m is the fraction of delayed neutrons of m-type, and

    β = ∑_m=1^Mβ_m.
 
The Marshak-type conditions are set at the boundary of the area ∂Ω:

    [ J_0,g( x); J_2,g( x);           ]    =
    [  -12  -38;  -38 -218;      ][ ϕ_0,g( x); ϕ_2,g( x);           ],
    
    
    J_i,g( x)     = -D_i,g∇ϕ_i,g( x), 
      
    i = 0, 2.


System of equations (<ref>) and (<ref>) is supplemented with boundary conditions (<ref>) and corresponding initial conditions:

    ϕ_g( x,0) = ϕ_g^0( x), 
        g = 1,2, ..., G,
        c_m( x,0) = c_m^0( x), 
        m = 1,2, ..., M.


Let's write the boundary problem (<ref>)–(<ref>) in operator form. 
The vectors u_1 = {ϕ_0,1, ϕ_0,2, ⋯, ϕ_0,G}, u_2 = {ϕ_2,1, ϕ_2,2, ⋯, ϕ_2,G}, c = {c_1, c_2, ..., c_M} and matrices are defined as follows

    V     = (v_gg'),
       
    v_gg' = 1/v_gδ_gg', 
    
    B     = (b_gg'),
       
    b_gg = -2Σ_r,g,
       
    b_gg' = 2Σ_s, g' → g, 
    
    A_1     = (a_gg'),
       
    a_gg = -∇· D_0,g∇ + Σ_r,g,
       
    a_gg' = -Σ_s, g' → g, 
    
    A_2     = (a_gg'),
       
    a_gg = -∇· D_2,g∇ + 5Σ_tr,g + 4Σ_r,g,
       
    a_gg' = -4Σ_s, g' → g, 
    
    F     = (f_gg'),
       
    f_gg' = χ_n,gνΣ_f,g', 
    
    E     = (e_gm),
       
    e_gm = χ_d,gλ_m, 
    Λ    = (λ_mm'), 
      λ_mm' = δ_mm'λ_m, 
    
    Q     = (q_mg),
     
    q_mg =β_m νΣ_f,g,

where

    δ_g g' = {
     1,     g = g', 
    
     0,     g ≠  g',
     .
  
is the Kronecker symbol.
We shall use the set of vectors u, whose components satisfy the boundary conditions (<ref>). 
Using the set definitions, the system of equations (<ref>) and (<ref>) can be written as following

    V  (d  u_1/d t - 2 d  u_2/d t ) + A_1  u_1 + B  u_2    =(1-β) F ( u_1 - 2 u_2) + E c,
    
    
    V  (- 2 d  u_1/d t + 9 d  u_2/d t )  + B  u_1 + A_2  u_2    =-2(1-β) F ( u_1 - 2 u_2) - 2E c,
    
    d  c/d t + Λ c    = Q ( u_1 - 2 u_2).

Without taking into account delayed neutrons (all neutrons are considered as prompt), we have

    V  (d  u_1/d t - 2 d  u_2/d t ) + A_1  u_1 + B  u_2    = F ( u_1 - 2 u_2),
    
    
    V  ( - 2 d  u_1/d t + 9 d  u_2/d t ) + B  u_1 + A_2  u_2    =-2 F ( u_1 - 2 u_2).

The Cauchy problem is formulated for equations  (<ref>) when

    u_1(0) =  u_1^0,    u_2(0) =  u_2^0,
 
where u_1^0 = {ϕ_0,1^0,  ϕ_0,2^0, ...,  ϕ_0,G^0 }, 
u_2^0 = {ϕ_2,1^0,  ϕ_2,2^0, ...,  ϕ_2,G^0 }. 
Для уравнений (<ref>)  задается также начальное условие

    c(0) =  c^0,
 
where c^0 = { c_1^0,  c_2^0, ...,  c_M^0 }.



§ SPECTRAL PROBLEMS


To characterize the reactor dynamic processes described by
Cauchy problem (<ref>)-(<ref>), let’s consider some spectral problems <cit.>.

The spectral problem, which is known as the λ-spectral problem, is usually considered.
For the system of equations (<ref>), (<ref>), we have

    L φ = λ^(k) M φ,

where

    φ = {φ_1, φ_2},
      
    L = [ A_1   B;   B A_2;     ],
      
    M = [   F -2F; -2F  4F;     ].

The minimal eigenvalue is used for characterisation of neutron field, thus

    k = 1/λ^(k)_1
 
is the effective multiplication factor (k-effective).
The value k = λ^(k)_1 = 1 is related to the critical state of the reactor, and the corresponding eigenfunction φ^(1)( x) is the stationary solution of the Eq (<ref>), (<ref>).
At k > 1, one can speak about supercriticality, at k < 1 — about subcriticality.

The spectral problem (<ref>) cannot directly be connected with the
dynamic processes in a nuclear reactor. The eigenvalues of the multiplication factor of the reactor and the corresponding eigenfunctions do not depend on the time delay for the emission of delayed neutrons. 
The reason is that the problem (<ref>) on eigenvalues is the problem of finding time-independent solutions of the neutron transport equation, and the term describing the contribution of fission to the neutron balance is equal to the total number of fission neutrons, both instantaneous and delayed divided by k.
At the best, we can get only the limiting case — the stationary critical state.
The more acceptable spectral characteristics for the non-stationary equation (<ref>) are related the spectral problem

    L φ - (1 - β) M φ - I  s    = λ^(α) W φ, 
    Λ s - R φ   = λ^(α) s.

The  α-eigenvalue problem without delayed neutrons (<ref>) can be written as follows

    L φ - M φ = λ^(α) W φ,

where

    I = [   E; -2E;     ],
      
    R = [   Q -2Q;     ],
      
    W = [   V -2V; -2V  9V;     ]

The fundamental eigenvalue

    α = λ^(α)_1

is called <cit.> the α-eigenvalue or the period eigenvalue, because it is inversely related to the reactor period.
The problem of the period eigenvalues essentially takes into account the contribution of delayed neutrons.
In particular, the long lifetime of the predecessors of delayed neutrons makes a large contribution to the slowly decreasing eigenfunctions of the reactor period, and this does not occur when only instantaneous neutrons are taken into account.

The asymptotic behaviour of Cauchy problem solution (<ref>)-(<ref>) at large times can be connected with the eigenvalue α.
In this regular mode, the reactor behaviour is described by the function exp(-α t) φ^(1)( x).
If α = 0, then the reactor is critical; if α > 0, then we get the neutron flux decreasing (subcritical state), and if α <  0, then we get the neutron flux increasing (supercritical state).

There is a simple approximate relationship between the dominant eigenvalues of the k-spectral (λ-spectral) and α-spectral problems without delayed neutrons <cit.>:

    1 - kk≈αΛ_pr,  
     1 - k ≈αł_pr,

where Λ_pr —  the prompt neutron generation time and  l_pr —  the prompt neutron lifetime. The prompt neutron generation time is determined as the neutron flux functional using adjoint fundamental α-eigenfunction rather then the k-eigenfunction <cit.>.

For the spectral problems with delayed neutrons, the relationship is as folows

    1 - kk≈α_totΛ_pr + α_tot∑_m=1^Mβ_mλ_m - α_tot,

which corresponds to inhour's equation. 
It is obvious that if we consider small perturbations of reactivity (| k - 1 | < β), then the dominant α_tot-eigenvalue does not practically depend on Λ_pr and is related  only with delayed neutrons parameters λ_m and β_m. Under such assumptions one can expect that the shapes of the fundamental eigenfunctions of different spectral problems are similar and, therefore, the relationships (<ref>) and (<ref>) are close to exact solution. In this case, we can derive the following relationship between α and α_tot.

    α≈α_totΛ_pr∑_m=1^Mβ_mλ_m - α_tot + α_tot,
 
where the second term in the rihght part of the equation (<ref>) can be neglected.

It is supposed that the relationships (<ref>)-(<ref>) are valid for the rest eigenvalues of the k-spectral and α-spectral problems <cit.>.



§ NUMERICAL EXAMPLES


To study the properties of the eigenvalues and eigenfunctions of dufferent types, several benchmarks are studied.
The first benchmark is the IAEA-2D two dimensional hexagonal problem of VVER-type without reflector <cit.>.
The second benchmark is the similar test but with radial reflector <cit.>. 
The aim of the third test is to investigate azimutally non-symmetric effects on the eigenvalues and eigenfunctions 
(the changed IAEA-2D test with reflector).
Finally, we studied the complex eigenvalues and eigenfunctions for the heavy water hexagonal reactor test HWR <cit.>.

The two-group model (G = 2) is used in all tests. 
The method of finite elements <cit.> on triangular calculation grids is used for the approximate solution of the spectral problem. 
The standard Lagrangian finite elements are used.
The software has been developed using the engineering and scientific calculation library FEniCS <cit.>.
SLEPc has been used for numerical solution of the spectral problems.
We used a Krylov-Schur algorithm with an accuracy of 10^-15.
The following parameters were varied in the calculations:
1pt 0pt 0pt

  * n — the number of triangles per one assembly (Fig. <ref>); 

  * p — the order of finite element.

In this work we compare the SP_3 calculations with the previous diffusion model calculations <cit.>.






 §.§ IAEA-2D without reflector


The geometrical model of the IAEA-2D reactor core <cit.> consists of a set of hexagonal assemblies and is presented in Fig. <ref>, where the assemblies of various types are marked with various digits. 
The total size of assembly equals 20 cm. 



Diffusion neutronics constants in the common units are given in Table <ref>.
The following delayed neutrons parameters are used: one group of delayed neutrons with effective fraction β_1 = 6.5·10^-3 and decay constant λ_1 = 0.08 s^-1. 
Neutron velocity  v_1 = 1.25 · 10^7 cm/s and v_2 = 2.5 · 10^5 cm/s.





  §.§.§ Solution of Lambda Modes spectral problem


As a reference solution for the diffusion model, we used the previous results <cit.>; for the SP_3 model — the solution obtained using very fine mesh (p = 3,   n = 96). The maximum difference in assembly power between two models is about 2 percent for the rodded assemblies (material 3, see Fig. <ref>).

The results of the solution of the effective multiplication factor for test IAEA-2D without a reflector are shown in Table <ref>.
Hereinafter, for λ-spectral problems, the following notation is used: k_dif — effective multiplication factor for the diffusion model; k_sp_3 — effective multiplication factor for the SP_3 model; Δ — absolute deviation from the reference value in pcm (10^-5); δ — the standard deviation of the relative power in percent. 
These data demonstrate the convergence of the computed eigenvalues with refinement of the calculation mesh and increase in polynomial degree.

The results of the first 10 eigenvalues for p = 3,   n = 96 are shown in Table <ref>.
The power and error distributions for the diffusion and SP_3 models are presented in Figs <ref> and <ref> for p = 2,  n = 24.
Hereinafter, for each assembly the following data are given: the reference solution (the diffusion or SP_3 model), the solution for p = 2,  n = 24 and the relative error from the reference solution.










  §.§.§ Solution of α-spectral problem without delayed neutrons


As a reference solution for both the diffusion and transport SP_3 models, we used the solutions obtained using very fine mesh (p = 3,   n = 96).
Hereinafter, for α-spectral problems, the following notation is used: α_dif — α-eigenvalue by diffusion model; α_sp_3 — α-eigenvalue by SP_3 model; Δ — absolute deviation from the reference value.

The calculation results for the α-spectral problem without  delayed neutrons using the different meshes and the finite element orders are shown in Table <ref>.
These data demonstrate the convergence of approximate computed eigenvalues with refinement of the calculation mesh and increase in polynomial degree.



The spectral problem results for the first 10 eigenvalues are shown in Table <ref>.
The eigenvalues λ_1^(α)≤λ_2^(α)≤ ... are well separated. 
In this example, the fundamental eigenvalue is less compared the rest and therefore the main harmonic  will attenuate more slowly.
A regular mode of the reactor is thereby defined.
The value α = λ_1^(α) determines the amplitude of neutron flux and is connected directly with reactor period in the regular mode.



The eigenfunctions for fundamental eigenvalue (i=1) of the α-spectral problem without delayed neutron are shown in Fig. <ref>. 
Due to the fact that the reactor state is close to critical (k=k_1≈ 0.97935), the fundamental eigenfunctions of the λ-spectral problem are close to the fundamental eigenfunctions of the α-spectral problem.
The eigenfunctions ϕ_1^(i), i=2,3,4,5 are shown in Fig. <ref>, Fig. <ref>.








  §.§.§ Solution of α-spectral problem with delayed neutrons


As a reference solution for both the diffusion and transport SP_3 models, we used the solutions obtained using very fine mesh (p = 3,   n = 96). 
The  α-spectral problem results with delayed neutrons using the different meshes and finite element orders are shown in Table <ref>.
Compared with the previous case without delayed neutrons, these data demonstrate the similar convergence of the computed eigenvalues.





The first 10 spectral problem eigenvalues are shown in Table <ref>.
Due to the contribution of delayed neutrons, the fundamental eigenvalue is much smaller compared with the case without delayed neutrons.

The eigenfunctions for fundamental eigenvalue (i=1) of the α-spectral problem with delayed neutron are shown in Fig. <ref>. 
The eigenfunctions ϕ_1^(i), i=2,3,4,5 are shown in Fig. <ref>, Fig. <ref>.
The eigenfunctions of the problems without and with delayed neutrons are close to each other in topology.





According to (<ref>)-(<ref>), the prompt neutron generation time Λ_pr = 4.920 · 10^-5 s  for both the diffusion and SP_3 fundamental eigenvalues.
Assuming the same value of Λ_pr for both α-spectral problems, the uncertainty in the fundamental eigenvalues is less than 0.1 percent.



 §.§ IAEA-2D with reflector


This test differs from the previous one only additional external row of reflector assemblies (material 4, see Table <ref>). 



  §.§.§ Solution of Lambda Modes spectral problem


As a reference solution for the diffusion model, we used the previous results <cit.>; for the SP_3 model — the solution obtained using the MCNP4C code <cit.>. As well as in the previous benchmark calculations, the maximum difference in assembly power between two models is about 2 percent for the rodded assemblies (material 3, see Fig. <ref>).

The comparison of the calculated effective multiplication factors is shown in Table <ref>. 
The results of the first 10 eigenvalues for p = 3,   n = 96 are presented in Table <ref>.
One can see that there are several eigenvalues of multiplicity two.

The power distributions and calculation errors for p = 2, n = 24 using the diffusion model are shown in Fig <ref> and for p = 3,   n = 96 using the SP_3 model are shown Fig <ref>.











  §.§.§ Solution of α-spectral problem without delayed neutrons


As a reference solution, we use the fine mesh solutions obtained using the diffusion or transport SP_3 model (p = 3,   n = 96). 
The α-spectral problem results are shown in Table <ref>.



The results of the first 10 eigenvalues for p = 3,   n = 96 are presented in Table <ref>.
As before, the eigenvalues are well separated.
In this example, the fundamental eigenvalue is negative and therefore the main harmonic will increase, while all others will attenuate. 



The eigenfunctions for fundamental eigenvalue (i=1) of the α-spectral problem without delayed neutrons are shown in Fig. <ref>. 
Due to the fact that a state of the reactor is close to critical (k=k_1≈ 1.00625), the fundamental eigenfunctions of the λ-spectral problem are close to the fundamental eigenfunctions of the α-spectral problem.
The eigenfunctions ϕ_1^(i), i=2,3,4,5 are shown in Fig. <ref>, Fig. <ref>.







  §.§.§ Solution of α-spectral problem with delayed neutrons


As a reference solution, we use the fine mesh solutions obtained using the diffusion or transport SP_3 model (p = 3,   n = 96). 
The α-spectral problem results are shown in
Table <ref>.





The calculation results for the first 10 eigenvalues are shown in Table <ref>.
Due to the contribution of delayed neutrons, the fundamental eigenvalue is much smaller than in the case without delayed neutrons. Again the fundamental eigenvalue is negative and therefore the main harmonic will increase, while all others will attenuate.

The eigenfunctions for fundamental eigenvalue (i=1) of the α-spectral problem with delayed neutrons are presented in Fig. <ref>. 
The eigenfunctions ϕ_1^(i), i=2,3,4,5 are shown in Fig. <ref>, Fig. <ref>.
The eigenfunctions of the problems without and with delayed neutrons are close to each other in topology.





According to (<ref>)-(<ref>), the prompt neutron generation time Λ_pr = 5.436 · 10^-5 s  for  the diffusion fundamental eigenvalue and Λ_pr = 5.416 · 10^-5 s  for the SP_3 fundamental eigenvalue. 
Assuming the same value of Λ_pr for both α-spectral problems, the uncertainty in the fundamental eigenvalues is less than 0.1 percent.



 §.§ Azimutally non-symmetric test IAEA-2D with reflector


To investigate azimutally non-symmetric geometry effects on the eigenfunction behaviour, we replaced two unrodded assemblies in the north-east part of the core by rodded ones (material 3, see Fig. <ref>). 





  §.§.§ Solution of Lambda Modes spectral problem


As a reference solution for both the diffusion and transport SP_3 models, we used the solutions obtained using very fine mesh (p = 3,   n = 96). 
The effective multiplication factors are shown in Table <ref>. 
The results of the first 10 eigenvalues for p = 3,   n = 96 are presented in Table <ref>.
As can be see from Table <ref> compared with Table <ref> for unperturbed case, all eigenvalues became well separated 
(pairs of the eigenvalues are vanished).







  §.§.§ Solution of α-spectral problem without delayed neutrons


As a reference solution for both the diffusion and  SP_3 models, we used the solutions obtained using very fine mesh (p = 3,   n = 96).
The α-spectral problem results are shown in Table <ref>.
The results of the first 10 eigenvalues for p = 3,   n = 96 are presented in Table <ref>.





Let's consider changes in the eigenfunctions due to the rodded assembly insertion.
The eigenfunctions for fundamental eigenvalue (i=1) of the α-spectral problem without delayed neutrons are shown in Fig. <ref>. 
The eigenfunctions ϕ_1^(i), i=2,3,4,5 are shown in  Fig. <ref>, Fig. <ref>.
As can be seen from Fig. <ref> to Fig. <ref>, the overal structure of the eigenfunctions is preserved taking into account the neutron flux perturbations.





According to (<ref>)-(<ref>), the prompt neutron generation time Λ_pr = 5.474 · 10^-5 s  for  the diffusion fundamental eigenvalue and Λ_pr = 5.451 · 10^-5 s  for the SP_3 fundamental eigenvalue. Thus, for the non-symmetric test we obtained similar neutronic properties compared with the symmetric test.
The corresponding fundamental eigenvalues of the α-spectral problem with delayed neutrons are calculated using (<ref>):
α_dif = -0.0679 and
α_sp_3 = -0.1078.




 §.§ HWR test problem


This benchmark is a model of large heavy-water reactor HWR <cit.>. 
The geometry of the HWR test is presented in Fig. <ref>. Fuel assemblis (1, 2, 3 and 6 in Fig. <ref>) located in the central part of the core, are surrouded by the target zone and reflector layer (7 and 9 in Fig. <ref>). There are two types of rodded assemblies (4 and 8).
The assembly size is equal to 17.78 cm. 

Diffusion constants  are given in Table <ref>. 
The following delayed neutrons parameters are used: one group of delayed neutrons with effective fraction β_1 = 6.5·10^-3 and decay constant λ_1 = 0.08 s^-1. 
Neutron velocity  v_1 = 1.25 · 10^7 cm/s and v_2 = 2.5 · 10^5 cm/s.



 



  §.§.§ Solution of Lambda Modes spectral problem



As a reference solution for the diffusion model, we used the  results obtained by <cit.>; for the SP_3 model — the solution obtained using very fine mesh (p = 3,   n = 96).

The effective multiplication factors for the HWR test are shown in Table <ref>. 



The results of the first 10 eigenvalues for p = 3,   n = 96 are presented in Table <ref>.
The eigenvalues k_2, k_3, k_4, k_5, k_9, k_10 of the λ-spectral problem are the complex values with small imaginary parts, and the eigenvalues k_1, k_6, k_7, k_8 are the real values. One can see pairs of the complex
conjugate values.



  §.§.§ Solution of α-spectral problem without delayed neutrons


As a reference solution for the diffusion and SP_3 models we use the solutions obtained using very fine mesh (p = 3,   n = 96).

The α-spectral problem results at different computational parameters are shown in Table <ref>. 
The results of the first 10 eigenvalues for p = 3,   n = 96 are presented in Table <ref>.
The eigenvalues are well separated.
The eigenvalues α_2, α_3, α_4, α_5, α_9, α_10 of the α-spectral problem, like for the λ-spectral problem, are the complex values with small imaginary parts, and the eigenvalues α_1, α_6, α_7, α_8 are the real values.

The eigenfunctions for fundamental eigenvalue (n=1) of the α-spectral problem  are shown in Fig. <ref>.
The real part of the eigenfunctions ϕ^(n)_1,  n = 2,3,4,5 is shown in Fig. <ref>.
Fig. <ref> shows the imaginary part of these eigenfunctions.
The eigenfunctions of the λ-spectral and α-spectral problems are close to each other in topology.











  §.§.§ Solution of α-spectral problem with delayed neutrons


As a reference solution for the diffusion and SP_3 models we use the solutions obtained using very fine mesh (p = 3,   n = 96).

The α-spectral problem results are shown in Table <ref>. 
Due to the contribution of delayed neutrons, the fundamental eigenvalue is much smaller than in the case without delayed neutrons.
The results of the first 10 eigenvalues for p = 3,   n = 96 is shown in Table <ref>.
The eigenvalues α_2, α_3, α_4, α_5, α_9, α_10 of the α-spectral problem, like as before, are the complex values with small imaginary parts, and the eigenvalues α_1, α_6, α_7, α_8 are the real values.





According to (<ref>)-(<ref>), the prompt neutron generation time Λ_pr = 1.917 · 10^-4 s  for both the diffusion and SP_3 fundamental eigenvalues.
Assuming the same value of Λ_pr for both α-spectral problems, the uncertainty in the fundamental eigenvalues is less than 0.1 percent.



§ CONCLUSION


Simulation of reactor dynamic processes is basically considered on the basis of multigroup diffusion approximation of the neutron transport equation. 
To improve the calculation accuracy for some situations of interest, including pin-by-pin calculations, the SP_3 approximation was adopted in different whole-core diffusion codes as an improved option compared with the diffusion approximation. 
In this regard, it will be very useful to compare the spectral parameters, calculated by both the diffusion and SP_3 options using the finite element method. 

Solution of the k-spectral (λ-spectral)  and α-spectral problems has been tested for the IAEA-2D and HWR reactor benchmark tests. 
The classical Lagrange finite elements are used for the spatial approximation. 
Accuracy control is performed using condensed grids. 
Spectral problems are solved numerically using well-developed free software SLEPc and using GMSH as a generic mesh generator.
 
Of particular interest is the problem associated with appearance of complex eigenvalues and eigenfunctions. 
It was found that this tendency occurs for both the diffusion and SP_3 solutions of the HWR reactor test. 
On the contrary, for the IAEA-2D benchmark tests we have not  found complex eigenvalues and eigenfunctions, even for the non-symmetric variant of the benchmark.

There were found some features in the k-spectral and α-eigenvalues, which were confirmed by calculations.
Using the known fundamental k-eigenvalue and the evaluated value of the prompt neutron generation time Λ_pr one can estimate any fundamental α-eigenvalue with and without delayed neutrons. The same behaviour occurs for the minor eigenvalues.



§ ACKNOWLEDGEMENTS


This work was supported by the grant of the Russian Federation Government (# 14.Y26.31.0013) and by the Russian Foundation for Basic Research (# 18-31-00315).



27
natexlab#1#1url<#>1urlprefixURL 

[Aragones et al.(2007)Aragones, Ahnert, and
  Garcia-Herran]Aragones2007FiniteDifference
Aragones, J. M., Ahnert, C., Garcia-Herran, N., 2007. The analytic coarse-mesh
  finite difference method for multigroup and multidimensional diffusion
  calculations. Nuclear Science and Engineering 157, 1–15.

[Avvakumov et al.(2017)Avvakumov, Strizhov, Vabishchevich, and
  Vasilev]avvakumov2017spectral
Avvakumov, A. V., Strizhov, V. F., Vabishchevich, P. N., Vasilev, A. O., 2017.
  Spectral properties of dynamic processes in a nuclear reactor. Annals of
  Nuclear Energy 99, 68–79.

[Avvakumov et al.(2014)Avvakumov, Vabishchevich, and
  Vasilev]avvakumov2014
Avvakumov, A. V., Vabishchevich, P. N., Vasilev, A. O., 2014. Metod konechnyh
  ehlementov dlya uravneniya diffuzii nejtronov v geksagonal'noj geometrii.
  Vestnik of North-Eastern Federal University 11 (5), 7–18, in Russian.

[Avvakumov et al.(2015)Avvakumov, Vasilev, and
  Zakharov]avvakumov2015
Avvakumov, A. V., Vasilev, A. O., Zakharov, P. E., 2015. Programmnaya
  realizaciya metoda konechnyh ehlementov dlya uravneniya diffuzii nejtronov.
  Vestnik of North-Eastern Federal University (4 (48)).

[Azmy and Sartori(2010)]azmy2010nuclear
Azmy, Y., Sartori, E., 2010. Nuclear computational science: a century in
  review. Springer.

[Bahabadi et al.(2016)Bahabadi, Pazirandeh, and Mitra]Bahabadi2016
Bahabadi, M. H., Pazirandeh, A., Mitra, A., 12 2016. Analytic function
  expansion nodal (AFEN) method for solving multigroup neutron simplified
  P_3 (SP_3) equations in hexagonal-z geometry. Annals of
  Nuclear Energy 98, 74–80.

[Beckert and Grundmann(2008)]beckert2008development
Beckert, C., Grundmann, U., 2008. Development and verification of a nodal
  approach for solving the multigroup SP_3 equations. Annals of
  Nuclear Energy 35 (1), 75–86.

[Bell and Glasstone(1970)]bell1970
Bell, G. I., Glasstone, S., 1970. Nuclear Reactor Theory. Van Nostrand Reinhold
  Company.

[Brantley and Larsen(2000)]brantley2000simplified
Brantley, P. S., Larsen, E. W., 2000. The simplified P_3
  approximation. Nuclear Science and Engineering 134 (1), 1–21.

[Brenner and Scott(2008)]brenner2008
Brenner, S. C., Scott, L. R., 2008. The mathematical theory of finite element
  methods. Springer.

[Brewster(2018)]Brewster2018
Brewster, W. J., 2018. Development and Monte Carlo validation of a finite
  element reactor analysis framework. Masters Theses 7756, Missouri University
  of Science and Technology.

[Chao and Shatilla(1995)]chao1995
Chao, Y. A., Shatilla, Y. A., 1995. Conformal mapping and hexagonal nodal
  methods-II: Implementation in the ANC-H Code. Nuclear Science and
  Engineering 121, 210–225.

[Downar et al.(2010)Downar, Xu, Seker, and Hudson]downar2010theory
Downar, T., Xu, Y., Seker, V., Hudson, N., 2010. Theory manual for the PARCS
  kinetics core simulator module. Department of Nuclear Engineering and
  Radiological Sciences, University of Michigan.

[Gelbard(1960)]gelbard1960application
Gelbard, E. M., 1960. Application of spherical harmonics method to reactor
  problems. Bettis Atomic Power Laboratory, West Mifflin, PA, Technical Report
  No. WAPD-BT-20.

[Gelbard(1961)]gelbard1961simplified
Gelbard, E. M., 1961. Simplified spherical harmonics equations and their use in
  shielding problems. Tech. rep., Westinghouse Electric Corp. Bettis Atomic
  Power Lab., Pittsburgh.

[Gelbard(1962)]gelbard1962applications
Gelbard, E. M., 1962. Applications of the simplified spherical harmonics
  equations in spherical geometry. Tech. rep., Westinghouse Electric Corp.
  Bettis Atomic Power Lab., Pittsburgh.

[Ginestar et al.(2002)Ginestar, Miro, Verdu, and
  Hennig]ginestar2002transient
Ginestar, D., Miro, R., Verdu, G., Hennig, D., 2002. A transient modal analysis
  of a BWR instability event. Journal of Nuclear Science and Technology
  39 (5), 554–563.

[Lawrence(1986)]lawrence1986progress
Lawrence, R. D., 1986. Progress in nodal methods for the solution of the
  neutron diffusion and transport equations. Progress in Nuclear Energy 17 (3),
  271–301.

[Logg et al.(2012)Logg, Mardal, and Wells]logg2012
Logg, A., Mardal, K.-A., Wells, G., 2012. Automated solution of differential
  equations by the finite element method: The FEniCS book. Vol. 84. Springer
  Science & Business Media.

[McClarren(2010)]mcclarren2010theoretical
McClarren, R. G., 2010. Theoretical aspects of the simplified P_N
  equations. Transport Theory and Statistical Physics 39 (2-4), 73–109.

[Quarteroni and Valli(2008)]quarteroni2008
Quarteroni, A., Valli, A., 2008. Numerical approximation of partial
  differential equations. Springer.

[Ryu and Joo(2013)]ryu2013fem
Ryu, E. H., Joo, H. G., 2013. Finite element method solution of the simplified
  P_3 equations for general geometry applications. Annals of Nuclear
  Energy 56, 194–207.

[Stacey(2007)]stacey2007
Stacey, W. M., 2007. Nuclear Reactor Physics. Wiley.

[Stewart(2001)]stewart2002krylov
Stewart, G. W., 2001. A Krylov–Schur algorithm for large eigenproblems. SIAM
  Journal on Matrix Analysis and Applications 23 (3), 601–614.

[Tada et al.(2008)Tada, Yamamoto, Yamane, and
  Kitamuray]tada2008applicability
Tada, K., Yamamoto, A., Yamane, Y., Kitamuray, Y., 2008. Applicability of the
  diffusion and simplified P_3 theories for pin-by-pin geometry of
  BWR. Journal of nuclear science and technology 45 (10), 997–1008.

[Verdu and Ginestar(2014)]verdu2014modal
Verdu, G., Ginestar, D., 2014. Modal decomposition method for BWR stability
  analysis using alpha-modes. Annals of Nuclear Energy 67, 31–40.

[Verdu et al.(2010)Verdu, Ginestar, Roman, and Vidal]verdu20103d
Verdu, G., Ginestar, D., Roman, J., Vidal, V., 2010. 3D alpha modes of a
  nuclear power reactor. Journal of Nuclear Science and Technology 47 (5),
  501–514.



