A nonstiff solution for the stochastic neutron point kinetics equations
Abstract
We propose an approach to solve the stochastic neutron point kinetics equations using an adaptation of the diagonalizationdecomposition method (DDM). This new approach (DoubleDDM) yields a nonstiff solution for the stochastic formulation, allowing the calculation of the neutron and precursor densities at any time of interest without the need of using progressive time steps. We use DoubleDDM to compute results for stochastic problems with constant, linear, and sinusoidal reactivities. We show that these results strongly agree with those obtained by other approaches established in the literature. We also compute and analyze the first four statistical moments of the solutions.
keywords:
point kinetics, stochastic moments, stiffness, decomposition methodauthoryear
1 Introduction
The neutron point kinetics equations (Hetrick, 1971; Kinard and Allen, 2004; Hayes and Allen, 2005) are the coupled differential equations for the neutron density and for the delayed neutron precursor concentrations. These equations model the timedependent behavior of a nuclear reactor and provide insight into the dynamics of its operation. The timedependent parameters in this system are the reactivity function and the neutron source term.
The neutron density and delayed neutron precursor concentrations vary randomly with time; however, the point kinetics equations are deterministic and can only be used to estimate average values. Random fluctuations in the neutron density and precursor concentrations can be significant at low power levels (Hurwitz Jr et al., 1963a, b), which points to the importance of estimating these variations.
Hayes and Allen (2005) have generalized the standard deterministic point kinetics equations, deriving a system of stochastic differential equations that model the random behavior of the neutron density and the precursor concentrations in a point reactor. Due to the issue of stiffness, this system was implemented numerically using a stochastic piecewise constant approximation method (Stochastic PCA). Work performed by Saha Ray (2012) has shown that order 1.5 strong Taylor and EulerMaruyama numerical methods are valid computational alternatives to Stochastic PCA in solving the stochastic point kinetics equations. However, with the exception of cases modeled with either none or only one group of neutron precursors, the stiffness of the problem remains.
In this paper we propose to solve this stochastic formulation using a double decomposition approach based on the diagonalizationdecomposition method (DDM) decribed by Wollmann da Silva et al. (2014). This proposed method is the major novelty and principal contribution of this work, yielding a nonstiff solution for the stochastic point kinetics equations. Specifically, this approach allows the calculation of the neutron and precursor densities at any time of interest without the need of using progressive time steps. This solution is obtained with a minimal amount of numerical approximations of the model; the largest numerical effort lies in the truncation of the decomposition and the integrations required by DDM.
The major caveat in this approach is that convergence of DDM is yet to be proven. For this reason, a Lyapunov criterion (Boichenko et al., 2005) is used to guarantee convergence (cf. Petersen et al., 2011; Wollmann da Silva et al., 2014). We present computational results for problems with constant, linear, and sinusoidal reactivities. The results of the proposed method are compared against those of other approaches established in the literature, showing strong agreement. We also compute the first four statistical moments of the solutions.
This work is an expanded version of a recent conference paper (Wollmann da Silva et al., 2015). The remainder of this paper is organized as follows. In Section 2 we present a brief review on the key aspects of DDM. In Section 3 we formulate the stochastic point kinetics equations. We introduce the proposed double decomposition approach in Section 4. Numerical results are given in Section 5 for problems with constant (Section 5.1) and timedependent (Section 5.2) reactivities. The paper concludes in Section 6 with a discussion of the work presented.
2 The DiagonalizationDecomposition Method (DDM)
Following Wollmann da Silva et al. (2014), one can obtain an analytical representation for the solution of the neutron point kinetics equations that is free of stiffness. The neutron point kinetics equations with six groups of precursors and timedependent reactivity are written as: \cref@addtoresetequationparentequation
(2.1a)  
(2.1b) 
for . Here, is the neutron density; is the density of the i delayed neutron precursor group; is the decay constant for a specific group ; represents the neutron mean generation time; and represents the delayedneutron fraction in a specific group . The total fraction of delayed neutrons is given by .
A recursive scheme with finite recursive depth is used to obtain a solution. The truncation index is determined with exponential convergence by the Lyapunov criterion (Boichenko et al., 2005; Petersen et al., 2011), evaluated after each recursion step. The neutron population and the precursors concentration are written in terms of the solution from a recursion initialization () and the respective correction terms () for an appropriate : \cref@addtoresetequationparentequation
(2.2a)  
(2.2b) 
The combination of Eq. 2.1 and (2.2) yields a system with unknowns. We define \cref@addtoresetequationparentequation
(2.3a)  
(2.3b)  
(2.3c)  
and  
(2.3d) 
where (constant) and are such that . Given the recursive system \cref@addtoresetequationparentequation
(2.4a)  
(2.4b) 
the solution of Eq. 2.4a is \cref@addtoresetequationparentequation
(2.5a)  
with . Equation 2.4b may be formally solved by the Laplace transform:  
(2.5b) 
since the initial condition from Eq. 2.1 is fully absorbed in Eq. 2.5a. The integral in Eq. 2.5b is evaluated using the GaussLegendre method.
A flowchart describing the implementation of this method is given in Fig. 1. The solution is obtained in an analytical representation that may be evaluated for any time value (free of stiffness).
3 The Stochastic Formulation
Hayes and Allen (2005) derived a system of Itô stochastic differential equations that model the dynamics of the neutron density and the delayed neutron precursors in a point nuclear reactor. This formulation describes the variation of the population and can be interpreted as a balance between deaths, births, and transformations of neutrons in the system. The probabilities of these events are determined by the physical parameters of the model, such as the total and partial delayed neutron fraction; the fraction of delayed neutrons of each precursor group; the decay constant of each group; and the average number of neutrons produced in each fission.
Assuming a time interval small enough such that only one event occurs, one can write \cref@addtoresetequationparentequation
(3.1a)  
where  
(3.1b)  
(3.1c)  
(3.1d)  
(3.1e)  
(3.1f)  
Here, and are given by Eqs. 2.3d and 2.3c, are Wiener processes, and  
(3.1g)  
(3.1h)  
(3.1i) 
with number of neutrons per fission. Note that if = 0, then Eq. 3.1a (with ) reduces to the deterministic problem discussed in Section 2.
4 The Proposed Method (DoubleDDM)
We propose to solve the stochastic formulation in Section 3 by adapting the recursive method described in Section 2. This approach yields a nonstiff solution to the stochastic system in Eq. 3.1.
Since the matrix depends on both the neutron populations and the delayed neutron concentrations, we resort to a double decomposition to obtain a solution for this problem:

is used to determine the matrix for a sequence of discrete time steps (its components are constant in each time step);

is obtained through diagonalization ( is symmetric);

Since is known for a specific time interval, a decomposition scheme analogue to DDM is applied: \cref@addtoresetequationparentequation
(4.1a) (4.1b) where are constants known in each time step.
The total number of stochastic components “drawn” in this approach is defined by the Central Limit Theorem (CLM) to guarantee a small statistical error for the first four moments. For all results shown in this paper, is large enough to guarantee a statistical error smaller than 0.05% with 95% confidence. The flowchart in Fig. 2 describes the implementation of the proposed method.
5 Numerical Results
In this section we present numerical results for the DoubleDDM approach proposed in Section 4 for examples with (i) constant and (ii) timedependent reactivities. We compare these results against those obtained with other approaches established in the literature. We calculate the expected value and variance Var), given by \cref@addtoresetequationparentequation
(5.1a)  
(5.1b)  
where the index represents the different choices of stochastic components (histories). We also calculate skewness and excess kurtosis for the neutron density , defined as  
(5.1c)  
(5.1d) 
which gives us further insight on the behavior of the stochastic solutions.
5.1 Constant Reactivity
In the following examples we present the results of four methods established in the literature: Monte Carlo and Stochastic PCA (Hayes and Allen, 2005); order 1.5 strong Taylor and EulerMaruyama (Saha Ray, 2012). These results are reproduced here as they were reported in the aforementioned references.
We compare these results with those obtained with the deterministic diagonalizationdecomposition method (DDM) and with the DoubleDDM approach. The solutions of the deterministic DDM are obtained by solving Eq. 3.1a with ; we point out that this is also the first step of DoubleDDM.
5.1.1 Stepreactivity insertion: one precursor
This example does not model an actual physical nuclear reactor problem. Nevertheless, considering only one group of precursors implies that the stiffness of the problem disappears; this provides a simple computational solution that is useful for a first comparison with other methods.
The physical parameters are , , , , , and for . The initial condition is . The expected values and standard deviations of and at time are presented in Table 1 for each of the methods.
Monte  Stochastic  Euler  Taylor  DDM  Double  

Carlo  PCA  Maruyama  1.5  DDM  
400.03  395.32  412.23  412.10  412.13  402.35  
27.311  29.411  34.391  34.519  –  28.610  
300.00  300.67  315.96  315.93  315.93  305.84  
7.8073  8.3564  8.2656  8.3158  –  7.9240 
A total of histories were accumulated for the DoubleDDM approach. Skewness and excess kurtosis for the neutron density were found to be and . The fact that and are nearly zero implies that the distribution of stochastic solutions is symmetric and has Gaussianlike peak and tail.
It can be seen that there exists a close agreement between DoubleDDM and the results obtained with Monte Carlo and Stochastic PCA. EulerMaruyama and order 1.5 strong Taylor (Taylor 1.5) yield slightly higher results, very close to those obtained with the deterministic DDM.
5.1.2 Stepreactivity insertion: six precursors
The following two examples model stepreactivity insertions for a thermal nuclear reactor. In this case, we consider a stiff system of equations with six precursor groups. The set of physical parameters is taken from Kinard and Allen (2004): , , , and , with and given in Table 2. The initial condition is given by
(5.2) 
Group  1  2  3  4  5  6 

0.266  1.491  1.316  2.849  0.896  0.182  
0.0127  0.0317  0.115  0.311  1.4  3.87 
We compute results for a prompt subcritical stepreactivity insertion at time , and for a prompt critical stepreactivity insertion at time . We define
(5.3) 
and present the expected values and standard deviations for each of the methods in Table 3.
Monte  Stochastic  Euler  Taylor  DDM  Double  

Carlo  PCA  Maruyama  1.5  DDM  
183.04  186.31  208.60  199.41  200.01  187.05  
168.79  164.16  255.95  168.55  –  167.83  
4.478  4.491  4.498  4.497  4.497  4.488  
1495.7  1917.2  1233.4  1218.8  –  1475.6  
135.67  134.55  139.57  139.57  139.61  135.86  
93.376  91.242  92.042  92.047  –  93.210  
4.464  4.464  4.463  4.463  4.463  4.463  
16.226  19.444  6.071  18.337  –  17.845 
We collected histories for the subcritical stepreactivity . Skewness was found to be , and excess kurtosis . For the critical stepreactivity , histories were collected. Skewness and excess kurtosis were computed respectively as and . These results for the third and fourth moments indicate that the distribution of the stochastic solutions is symmetric and has neither a sharp peak nor a heavy tail.
As in the previous example, the results obtained with DoubleDDM are in close agreement to the results from Monte Carlo and Stochastic PCA. The results from EulerMaruyama and order 1.5 strong Taylor are closer to those of deterministic DDM.
5.2 TimeDependent Reactivity
The current literature lacks numerical results for the stochastic system in Eq. 3.1a with timedependent reactivities. For this reason, the results collected from the literature and reproduced next represent only the deterministic solution (with ). Although not ideal, this approach allows us to verify that the expected value obtained with DoubleDDM closely agrees with well established models for problems with timedependent reactivities. All the following examples take into account six precursor groups.
5.2.1 Linear reactivity
The following two examples model a ramp reactivity for a thermal nuclear reactor. The physical parameters considered are: , , and , with and taken from Lewins (1978) and given in Table 4.
Group  1  2  3  4  5  6 

0.246  1.363  1.203  2.605  0.819  0.167  
0.0127  0.0317  0.115  0.311  1.4  3.87 
We compute the neutron density for two different choices of constant : and . The results obtained with DoubleDDM are given in Table 5 for times and . We compare these results with those obtained with the Padé approximation (Aboanber and Nahla, 2002) and the generalization of the analytical exponential model (GAEM), as reported by Nahla (2008).
Time  Padé  GAEM  DDM  Double  

DDM  
0.25  1.069840  1.069541  1.069542  1.069763  
0.50  1.157065  1.156694  1.156695  1.157867  
0.75  1.265795  1.265331  1.265332  1.269374  
1.0  1.402562  1.401981  1.401982  1.403561  
0.25  1.149544  1.149200  1.149210  1.137216  
0.50  1.369438  1.368927  1.368928  1.356934  
0.75  1.708411  1.707600  1.707601  1.695607  
1.0  2.276692  2.275271  2.275272  2.263278 
The number of histories collected for the case was . We computed the higher moments for time , finding the standard deviation , skewness , and excess kurtosis . For the case , we collected histories, and found , , . This shows that, in both cases, the distribution of stochastic solutions is nearly normal.
DoubleDDM shows good agreement with the other methods shown in Table 5. The results obtained with DoubleDDM are slightly larger for the first case (), and slightly smaller for the second case. The relative differences are around or smaller, being well within 1 standard deviation.
5.2.2 Sinusoidal reactivity
The last example simulates a case with sinusoidal reactivity , with , , and . The total fraction of delayed neutrons is given by , with and shown in Table 6.
Group  1  2  3  4  5  6 

0.214  1.423  1.247  2.568  0.748  0.273  
0.0124  0.0305  0.111  0.301  1.14  3.01 
The results obtained with DoubleDDM are presented in Table 7 for every whole second up until . We compare these results with those reported in (Wollmann da Silva et al., 2014). These were obtained with the method introduced by Kang and Hansen (1973), referred to as K & H in Table 7, and with the method of Taylor series (cf. Nahla, 2011).
Time  K & H  Taylor  DDM  Double 

DDM  
0  1.00000  1.00000  1.00000  1.0000000 
1  1.12397  1.12394  1.12396  1.1119659 
2  1.16881  1.16884  1.16889  1.1568959 
3  1.07443  1.07442  1.07448  1.0624859 
4  0.95381  0.95380  0.95382  0.9418259 
5  0.90737  0.90737  0.90735  0.8953559 
6  0.96151  0.96158  0.96153  0.9495359 
7  1.08748  1.08749  1.08745  1.0754559 
8  1.17168  1.17164  1.17167  1.1596759 
9  1.11128  1.11124  1.11130  1.0993059 
10  0.98464  0.98464  0.98468  0.9726859 
We collected histories to achieve the requirement imposed for the statistical error (less than 0.05% with 95% confidence). The higher moments for time yield , , and . These results indicate (i) a very small asymmetry in the distribution of stochastic solutions, with a slightly larger left tail; and (ii) a very small yet noticeable “flatter” peak when compared to a normal.
In general, DoubleDDM closely agrees with the other methods presented here for comparison. Results displayed in Table 7 show that DoubleDDM yields slightly smaller results than those attained with the other methods. This can be confirmed in Fig. 3, which depicts the average behavior of the stochastic solution compared to the deterministic solution (DDM). These differences are, once again, very small (), and well within 1 standard deviation.
6 Conclusions
In this paper, we propose an approach to solve the stochastic neutron point kinetics equations with a solution procedure that is free of stiffness. This is achieved through an adaptation of the diagonalizationdecomposition method (DDM) introduced in Wollmann da Silva et al. (2014), wich provides a nonstiff solution for the deterministic point kinetics equations. DDM uses the Laplace transform to obtain a formal solution, then applies a decomposition into a recursive scheme, using GaussLegendre integration.
The essential steps of the proposed approach (DoubleDDM) are: (i) the deterministic problem is solved with DDM; (ii) the deterministic solution is used to build the stochastic component; (iii) another decomposition scheme analogue to DDM is used to solve the stochastic system. This allows the calculation of the neutron density and precursor concentrations at any time of interest, without the need to resort to progressive time steps. The elimination of stiffness comes from the fact that the evolution of the solution by recursion adds correction terms to the whole time interval of interest in each step, and simultaneously for each term that depends on a specific time scale.
Since convergence of DDM is yet to be proven, a Lyapunov criterion is used to guarantee convergence. The results of the proposed method are compared against results obtained through other approaches established in the literature. This comparison shows close agreement for problems with constant stepreactivity insertions, as well as timedependent ramp reactivity and sinusoidal reactivity insertions.
In the current literature, the stochastic problem is mainly solved for constant reactivities; numerical solutions are limited to feasible time intervals due to the stiffness inherent to the problem. The mitigation of the stiffness character in solving the stochastic formulation is the major novelty and principal contribution of this work. Moreover, the analysis of the third and fourth moments of the stochastic solutions is, to the best of our knowledge, new.
The analyzed moments still depend on the size of the sample set and on the frequency with which the stochastic fluctuations are applied. In principle, an adjustment such as variance reduction and its consequences for higher moments could yield more realistic results. It would be necessary to find a reference scale in order to obtain such results independently of the sample size or frequency of application. This still needs to be identified, as well as a necessary ingredient to mimic reactor fluctuations. These tasks are left for future work.
Acknowledgments
Milena Wollmann da Silva would like to thank the Coordination for the Improvement of Higher Education Personnel (CAPES, Brazil) for financial support. Bardo E.J. Bodmann and Marco Túllio Vilhena thank the Federal University of Rio Grande do Sul (UFRGS) and also acknowledge financial support from the National Council of Scientific and Technological Development (CNPq, Brazil). Richard Vasques prepared this paper under award number NRCHQ8414G0052 from the Nuclear Regulatory Commission. The statements, findings, conclusions, and recommendations are those of the authors and do not necessarily reflect the view of the U.S. Nuclear Regulatory Commission.
References
References
 Aboanber and Nahla (2002) Aboanber, A. E., Nahla, A. A., 2002. Solution of the point kinetics equations in the presence of newtonian temperature feedback by padÃ© approximations via the analytical inversion method. Journal of Physics A: Mathematical and General 35 (45), 9609–9627.
 Boichenko et al. (2005) Boichenko, V. A., Leonov, G. A., Reitmann, V., 2005. Dimension Theory for Ordinary Differential Equations. Teubner Verlag.
 Hayes and Allen (2005) Hayes, J. G., Allen, E. J., 2005. Stochastic pointkinetics equations in nuclear reactor dynamics. Annals of Nuclear Energy 32 (6), 572–587.
 Hetrick (1971) Hetrick, D. L., 1971. Dynamics of Nuclear Reactors, 1st Edition. University of Chicago Press, Chicago, U.S.A.
 Hurwitz Jr et al. (1963a) Hurwitz Jr, H., MacMillian, D. B., Smith, J. H., Storm, M. L., 1963a. Kinetics of low source reactor startups part i. Nuclear Science and Engineering 15 (2), 166–186.
 Hurwitz Jr et al. (1963b) Hurwitz Jr, H., MacMillian, D. B., Smith, J. H., Storm, M. L., 1963b. Kinetics of low source reactor startups part ii. Nuclear Science and Engineering 15 (2), 187–196.
 Kang and Hansen (1973) Kang, C. M., Hansen, K. F., 1973. Finite element methods for reactor analysis. Nuclear Science and Engineering 51 (4), 456–495.
 Kinard and Allen (2004) Kinard, M., Allen, E., 2004. Efficient numerical solution of the point kinetics equations in nuclear reactor dynamics. Annals of Nuclear Energy 31 (9), 1039–1051.
 Lewins (1978) Lewins, J., 1978. Nuclear Reactor Kinetics and Control. Pergamon Press, Oxford, UK.
 Nahla (2008) Nahla, A. A., 2008. Generalization of the analytical exponential model to solve the point kinetics equations of be and d2omoderated reactors. Nuclear Engineering and Design 238 (10), 2648–2653.
 Nahla (2011) Nahla, A. A., 2011. Taylor series method for solving the nonlinear point kinetics equations. Nuclear Engineering and Design 241 (5), 1592–1595.
 Petersen et al. (2011) Petersen, C. Z., Dulla, S., Vilhena, M. T., Ravetto, P., 2011. An analytical solution of the point kinetics equations with timevariable reactivity by the decomposition method. Progress in Nuclear Energy 53 (8), 1091–1094.
 Saha Ray (2012) Saha Ray, S., 2012. Numerical simulation of stochastic point kinetic equation in the dynamical system of nuclear reactor. Annals of Nuclear Energy 49, 154–159.
 Wollmann da Silva et al. (2014) Wollmann da Silva, M., Leite, S. B., Vilhena, M. T., Bodmann, B. E. J., 2014. On an analytical representation for the solution of the neutron point kinetics equation free of stiffness. Annals of Nuclear Energy 71, 97–102.
 Wollmann da Silva et al. (2015) Wollmann da Silva, M., Vilhena, M. T., Bodmann, B. E. J., Vasques, R., 2015. The solution of the neutron point kinetics equation with stochastic extension: an analysis of two moments. In: Proc. International Nuclear Atlantic Conference.