Acceleration of a 2-dimensional, 2-energy group neutron noise solver based on a discrete ordinates method in the frequency domain

The acceleration of the convergence rate is studied for a neutron transport solver to simulate 2-D, 2-energy-group neutron noise problems in the frequency domain. The Coarse Mesh Finite Difference (CMFD) method is compared to the Diffusion Synthetic Acceleration (DSA) method. Numerical tests are performed using heterogeneous system configurations with different boundary conditions. The CMFD scheme leads to a better convergence rate. The results also show that CMFD can accelerate neutron noise problems in a wide range of perturbation frequencies with almost equal efficiency. An unstable convergence behavior is nevertheless observed in problems with purely reflective boundary conditions. Stabilization techniques such as performing multiple transport sweeps, underrelaxing the flux update, and using the lpCMFD method are investigated and improvements can be obtained.


INTRODUCTION
Reactor neutron noise is referred to as the fluctuations of the neutron flux around an expected mean value because of perturbations of the nuclear reactor properties.The analysis of the neutron noise is helpful to core monitoring and diagnostics, and it subsequently contributes to enhanced safety.For this purpose, simulations of the system response induced by possible perturbations are often required.Most of the past work in the area of modelling and simulations of neutron noise problems relies on neutron diffusion theory [1].Recent research focus on high order transport methods [2,3,4].These methods can provide more accurate insights and can be used to assess the limitations of the diffusion model.A neutron noise simulator is under development at Chalmers University of Technology and is based on a high order transport method for the solution of the neutron noise balance equations in the frequency domain.The simulator is based on discrete ordinates (SN) method and its current version can deal with two-dimensional, two-energy group problems [5].The overall solution algorithm consists of 2 steps: the first step solves the criticality problem and the second step determines the neutron noise in the frequency domain.The application of the conventional inner-outer iterative algorithm to the frequency domain neutron noise calculations results in very slow convergence rate.Thus, the numerical acceleration of the algorithm is required.The acceleration for static transport calculations have been widely studied; on the other hand, an effective acceleration of the neutron noise calculations in the frequency domain is an open question.A first effort was to test the DSA method.Such a method leads to a significant improvement in the overall convergence rate of the neutron noise simulator.However, the frequency-domain calculation still requires a relatively large number of iterations to converge and its acceleration deteriorates as the frequency of the perturbation becomes lower.In order to overcome these issues, the CMFD method is considered for the acceleration of both the static problem and the neutron noise solution in the frequency domain.
In this paper, the formulation of the CMFD method for the case of frequency-domain neutron noise calculations is presented.Then the numerical performance of this approach is compared to the DSA method.

THE FREQUENCY DOMAIN TRANSPORT NEUTRON NOISE EQUATION
A neutron multiplying system is considered with a perturbation that can be modelled as small, stationary fluctuations of the neutron cross sections.The induced neutron noise can be estimated in the frequency domain from the solution of the following equation: where the noise source term S g �r ⃗, Ω � , ω� is expressed as: and the energy spectrum    ( ⃗, ) is: Eq. ( 1) is derived from the time-dependent, multi-energy group neutron transport equation with a generic number of families of precursors of delayed neutrons (more details are discussed in [5]).The solution of Eq. ( 1) is obtained from a fixed source problem and provides the scalar neutron noise δ  and the angular neutron noise   as complex values.The perturbations are assumed to occur in a critical system and to not affect the value of k eff .

FREQUENCY DOMAIN NEUTRON NOISE SOLVER WITH CMFD ACCELERATION
The transport neutron noise solver under development at Chalmers University of Technology is based on the 2-energy-group, 2-dimensional version of Eq. ( 1) and is arranged in a static and a dynamic module.In both modules, the transport equations are discretized angularly with the discrete ordinates method and spatially with the diamond finite difference scheme.The conventional inner with-in group scattering sweeps and outer fission source iterations are used.
For the acceleration of the static module, both DSA [6] and CMFD [7] methods are implemented.In the dynamic module, the DSA method was implemented as a first attempt to accelerate the dynamic calculations [5].The analysis of the DSA-based scheme showed that the convergence rate can be improved, although the number of iterations is still high.In addition, the acceleration was found to degrade with the decrease of frequency.Thus, the CMFD method, which is widely used in accelerating static [7] and timedependent [8] calculations, is investigated.
The formulation of the CMFD scheme for neutron noise calculations in the frequency domain follows the formulation used for the static case.The CMFD equations for the dynamic calculations are derived through a homogenization process that results in lower-order diffusion-like equations which preserve the different terms in the original transport problem.This process is briefly described in the following.
At the ( + 1)-th outer iteration, the transport calculation provides the scalar and angular neutron noise fluxes ϕ ,, ,(+1/2) and  ,,, ,(+1/2) for all the energy groups  = 1, … , .The indices  and  are used for the fine transport mesh along the  and  directions, respectively.The index  denotes the discrete ordinates which is defined by the cosine values   and   for the  and  directions, respectively, with weights   . Then are operated on the fully discretized 2-dimensional version of Eq. ( 1) successively.The volume  , of the fine cell (, ) equals to ∆  ∆  .The spatial indices ,  are associated with the coarse CMFD mesh.Such operations lead to: According to the CMFD approach, the current terms are approximated by Fick's Law together with a transport correction.The current in the  direction is approximated as: where the normal coupling coefficient  � and the correction coefficient  � ,+1/2, are defined by: and ,(+1/2) − Φ ,, ,(+1/2) � �Φ ,+1, ,(+1/2) + Φ ,, ,(+1/2) � In Eq. ( 6), the diffusion coefficients of the coarse cells are calculated using the total cross-sections of the coarse cells that are homogenized with the static fluxes: In Eq. ( 7), δψ ,+1/2, and Φ ,+1, ,(+1/2) are, respectively, the surface averaged and cell averaged (homogenized) angular and scalar neutron noise calculated from the transport sweep.In the calculation of the coupling coefficients, the treatment of different boundary conditions is similar to the static case [7].Similar expressions are used for the current in the  direction.
In the derivation of Eq. ( 4) the term related to the time derivative is added to the total cross-section term, i.e.
The same homogenization procedure is applied to the cross-sections  , ′ →,, and νΣ f,g ′ ,i,j .The term χ � g,p,q dyn in the fission source term of Eq. ( 3) is homogenized according to: The noise source term is homogenized with only the volume of the mesh.
Eq. ( 4) also represents a fixed source problem, and this linear system of equations is solved in matrix form for the coarse mesh complex valued scalar neutron noise flux Φ ,, ,(+1) .Since the system considered in this work is relatively small, the corresponding matrix is simply inverted by applying an LU factorization.
The CMFD method is known to be unstable for some static problems [9], and especially for those having high scattering ratios and coarse cells with optical large thicknesses.In the CMFD algorithm for dynamic calculations, unstable behavior is also observed under some circumstances and therefore three stabilizing techniques are tested.The first technique is simply performing more than one transport sweeps in each overall iteration before passing the flux and current information to the CMFD calculation.The increase of transport sweep can be applied by performing either more inner scattering sweeps per outer iteration or more outer iterations where each outer iteration includes one inner sweep.The second technique tested is the Linear Prolongation CMFD (lpCMFD) method [10], in which the flux update equation is modified as follows: ϕ ,, ,(+1) = ϕ ,, ,(+1/2) + ϕ ,, ,(+1/2) The correction quantity ϕ ,, ,(+1/2) is obtained by first calculating the difference between the coarse mesh CMFD and transport fluxes at the vertices and then interpolating them to the centers of the fine mesh cells.The third stabilizing technique relies on flux update relaxation.In this case, the flux update equation (Eq.( 11)) is replaced by: The value of the underrelaxation parameter θ varies between zero and unity.The unaccelerated scheme is obtained when θ = 0.The standard CMFD update and thus Eq. ( 11) are recovered when θ = 1.For the underrelaxation and lpCMFD method, only one transport sweep is performed per overall iteration.

NUMERICAL RESULTS
The CMFD and DSA accelerated versions of the dynamic solver are tested using the heterogeneous C3 and C4V benchmark configurations reported in [10].Both systems consist of a 2 × 2 squared array of fuel assemblies with two UO2 assemblies and two MOX assemblies.Each fuel assembly contains 17 × 17 homogenized fuel cells.The C3 problem has reflective boundary condition on all the 4 sides of the system, while the C4V case has reflective and vacuum boundary conditions depending on the side of the system.
The neutron noise source is assumed to be a stationary fluctuation of the capture cross-section in both energy groups, placed in a fuel cell of one of the MOX fuel assemblies.The amplitude of the perturbation is taken to be 5% of the nominal values of the capture cross-sections for each group.The frequency of the neutron noise source is set to 1 Hz.
In the calculations presented in this section, the transport sweeps are carried out over a fine mesh in which each fuel cell of size 1.26 × 1.26 cm is discretized with 3 × 3 equally sized cells, and a Level-symmetric S8 quadrature set is applied.Numerical tests with other orders of discrete ordinates such as S16 and S20, lead to similar results.For the CMFD calculations, the coarse mesh size corresponds to the size of each fuel cell.The convergence is checked for the real part, imaginary part, amplitude and phase of the scalar neutron noise flux in a point wise manner.The calculation is stopped when the relative differences between two iterations, for all the quantities, are lower than 1E-6.The solution is then used to evaluate the residuals.For each of the cases discussed below, the Euclidean norm of the point-wise residuals, normalized with respect to the noise strength, is about 1E-5 or lower.
For the C4V problem the dynamic solver is accelerated by the CMFD method with 1 SN outer iteration and with 2 SN outer iterations (where only 1 inner iteration is run per each outer iteration), by the lpCMFD method, and by the DSA method.The neutron noise calculated with the accelerated methods are compared to the one calculated with the unaccelerated solver taken as reference, where the relative differences are evaluated in a pointwise manner for the amplitude and phase of the scalar neutron noise flux.The solutions from the DSA and CMFD methods are almost identical to the unaccelerated solution, with a maximum difference of 0.036%.In Table I the performances of the different methods are reported.The total number of space-angle sweeps needed for each energy group to reach the criterion for stopping the calculation, shows that the CMFD-based algorithms are numerically more efficient than the DSA scheme.For the different methods, the relative changes in the real and imaginary part with respect to the iteration number are plotted in Fig. 1.The standard CMFD with 1 or 2 outer iterations have a similar convergence behavior.
The lpCMFD method requires a slightly more computational time.Since the results are satisfactory, no underrelaxation is applied to this problem.In Table II the results for the C3 configuration (where all the boundary conditions are reflective) are given.
In this case the CMFD method with under-relaxion factor  = 0.5 or 0.18 are also considered.The most efficient convergence rate is obtained from the lpCMFD method.The CMFD methods becomes unstable as shown in Fig, 2, where the changes of the relative differences in the real and imaginary part are plotted with respect to the iteration number.Since the C3 and C4V configurations differ in only their boundary conditions, the unstable behavior may thus be related to the treatment of the reflective boundary condition in the algorithms.All the methods provide similar results in terms of calculated neutron noise, since the maximum relative differences with respect to the unaccelerated solution are lower than 0.04%.The efficiency of the CMFD acceleration method is also investigated for different frequencies of the noise source.For the C4V configuration the CMFD method with 1 transport sweep is chosen because it is performing better than the other methods considered.For the C3 configuration, the lpCMFD method which lead to fastest convergence rate, is selected.The total numbers of transport sweeps necessary for convergence are summarized in Table III.For the C4V case, the standard CMFD method with 1 transport sweep provides excellent acceleration effect at any frequency, even at very low frequencies which are numerically challenging for other methods.For example, DSA method would require more than 3000 iterations at 0.001 Hz.In the C3 case, the lpCMFD performs well in the region where the neutron noise amplitude does not depend on the frequency.For very low frequencies the performance deteriorates.For very high frequencies (e.g.1000 Hz) large oscillations in the change of the relative difference with respect to the number of iterations, for the real part, are observed and the solver cannot converge.This unstable behavior is also found when other stabilized CMFD methods were used.

CONCLUSIONS
The acceleration of a frequency domain neutron noise solver based on a discrete ordinates method was investigated using several variants of the CMFD method.The convergence properties of the CMFD accelerated algorithms were compared with the DSA accelerated algorithm.Two heterogeneous benchmark problems that differ in the boundary conditions were used for the study.In the C4V problem (where reflective and vacuum boundary conditions are applied), the CMFD methods outperform the DSA method in terms of convergence rate, while predicting very similar values of the calculated neutron noise.In the calculations with the C3 configuration (where all the boundary conditions are reflective), the CMFD methods provide faster convergence than the DSA method, but issues on the convergence stability are observed.For extremely low and high frequencies, the CMFD methods may become unstable and even fail at high frequencies.Future work is required to improve the CMFD-based scheme for problems with pure reflective boundaries.Moreover, the noise solver is planned to be extended so that the energy variable can be discretized with a generic number of energy groups.

Figure 1
Figure 1 C4V problem -Convergence of the real (on the left) and the imaginary (on the right) part of the calculated neutron noise

Table I .
C4V problem -Acceleration of the dynamic solver with CMFD and DSA methods

Table II .
C3 problem -Acceleration of the dynamic solver with CMFD and DSA methods

Table III
Performance of the CMFD method for different frequencies of the noise source