Generalization of the analytical solution of neutron point kinetics equations with time-dependent external source

Point reactor kinetics equations with one group of delayed neutrons in the presence of the time-dependent external neutron source are solved analytically during the start-up of a nuclear reactor. Our model incorporates the random nature of the source and linear reactivity variation. We establish a general relationship between the expectation values of source intensity and the expectation values of neutron density of the sub-critical reactor by ignoring the term of the second derivative for neutron density in neutron point kinetics equations. The results of the analytical solution are in good agreement with the results obtained with numerical solution.


Introduction
The neutron point kinetics (NPK) equations are a system of coupled nonlinear and stiff equations [1][2][3]. This system describes the nuclear parameters such as neutron density, reactivity and the precursor concentrations of delayed neutrons [4]. Calculations of these parameters are concerned with the reactor dynamics. The neutron density is one of the most important parameters in reactor dynamics [5][6][7]. According to the importance of neutron density in the cold start-up stage, the external neutron source makes a significant contribution to the reactor power [8]. In most reactors, the reactivity is introduced mainly by discontinuous transferring of control rods which is, in practice, a linear function. Applying too much reactivity will cause a high increase in reactor power and can result in an overpressure accident in the reactor core [9]. Therefore, one should carefully study the influence of external neutron source on the reactor power [10,11]. So far, appropriate mathematical models have already been developed to study the sub-critical kinetics, such as including the constant external neutron source [6], approximation models [12][13][14] and prompt jump approximation (PJA) [9,15,16].
In this work we extend the analytical method to determine the effect of external source on the behavior of the neutron density during the start-up of a nuclear power reactor. The source of neutron is considered as a function of time.
This paper is organized as follows. In Sect. ''Start-up process'', a brief description of the start-up process is presented. In Sect. ''External neutron source contribution'', the mechanism of interaction between the control rod and external source is reviewed, and finally the authors discuss and interpret the results. The section is also followed by two appendices.

Start-up process
The start-up process of a reactor deals with the variation of reactivity in the system which occurs by raising the control rods in a discontinuous way [12]. So, the rule of neutron density variation is significant because it not only describes the relationship between lifting speed of the control rod, duration and speed of neutron density response, but also is helpful to the operator to control the reactor start-up from any accident [9,17,18]. The authors modified the analytical expressions for the calculation of neutron density from a set of physical and mathematical approximations [9,12,15]. Since the reactor is in sub-critical state in the cold start-up stage, the external neutron source cannot be ignored. Due to the low average temperature in the reactor core, here, the effect of temperature feedback can be ignored [9,19]; consequently, the NPK equations with the external neutron source term for one group of delayed neutrons can be written as follows [8,11,12,16]: where nðneutron cm À3 Þ is the neutron density, cðneutron cm À3 Þ is the precursor concentrations, b is the total fraction of delayed neutrons, kðs À1 Þ is the decay constant of delayed neutrons, lðsÞ is the prompt neutron generation time, q is the reactivity as a function of time, and qðneutron cm À3 s À1 Þ is the external neutron source.
To have a general solution, the variation of external neutron source with respect to time is considered. Let us first consider the external neutron source defined as the ratios of polynomials of degree N: We denote their reactivity by: where q 0 is the initial subcritical reactivity, rðs À1 Þ represents the ramp reactivity addition rates and t 0 ðsÞ is the time for taking out the rods. Eliminating the dependency of precursor concentration, the differential equation that governs the neutron density during the withdrawal of the control rods is given by [9,12]: Since l d 2 nðtÞ dt 2 is much smaller than the other terms in Eq. (3), we ''Ignore the second term derivative'' (ISD). So, the only approximation which we used in the calculations is ISD [9,12]. According to the ISD, one can write: According to the data of some thermal reactors, i.e. k ¼ 0:001 s À1 ; b ¼ 0:0075; q 0 ¼ À0:006; l ¼ 0:0015 s, the assumption of b À q ) kl for very small ramp rate reactivities is valid [12]. To find the analytical solution, provided n ¼ n 1 , q 1 ¼ q 0 þ rt, k 1 ¼ kq 0 þr r and k 2 ¼ bÀq 0 r , the equivalent form of Eq. (4) is rearranged as follows: Considering (u ¼ t À k 2 , g ¼ k 1 þ kk 2 À 1) [12], Eq. (5) reduces to: n m k m 2 ðkq n þðn þ 1Þq nþ1 Þu nÀm þ AÞ: One can write (for more details see ''Appendix A''): It is possible to determine the integration constant ðAÞ by using initial condition, ðn 1 ð0Þ ¼ n 0 Þ, and so we have: For t ¼ t 0 , the reactivity will reach a constant value ðq 2 ¼ q 0 þ rt 0 Þ. At this time, neutron density is equal to n 2 ðtÞ ¼ n 2 . Therefore, in this stage of introducing reactivity, the differential equation is given by: where s ¼ kðq 0 þrt 0 Þ q 0 þrt 0 Àb and b n ¼ À lðkq n þðnþ1Þq nþ1 Þ q 0 þrt 0 Àb . Finally, the solution of Eq. (8) can be obtained from the integrating factor method (for more details see ''Appendix B''): The integral constant can be calculated by applying the continuity condition (n 2 ðt 0 Þ ¼ n 1 ðt 0 Þ): In summary, utilizing Eqs. (6) and (9), the neutron density was calculated when the core of the reactor was interpolated by an external source at the start-up stage. In the following, the neutron density was calculated for the case of a sinusoidal source, which is very close to a real system. In addition, the constant source was also considered to examine the validity of the obtained results.

External neutron source contribution
The analytical solution is calculated with two terms: one as a function of sinusoidal fluctuation in the external source, and the other as a function of constant source. Using the sinusoidal solution, the importance of the source term in the reactor kinetics analysis is better manifest.

Constant external neutron source
By calculating the neutron density for constant source, ðq ¼ q 0 Þ, we have tried to show that the obtained result reduces to Palma et al.'s result [12]. From density equation definition of Eq. (6), we can write: where constants a 00 and A are defined by: The complete solution of Eq. (8) with external constant source is given by: In deriving n 2 ðtÞ, we used the continuity condition. In Eq. (12), B is the integration constant given in terms of the incomplete gamma function as: Therefore, Eqs. (11) and (12) are in agreement with previous results [Eqs. (16), (24) [12]]. These solutions are expressed in terms of incomplete gamma function and related to the probability integral of the chi-squared distribution, as discussed and tabulated by Abramowitz and Stegun [20].

Sinusoidal fluctuation
Due to fluctuation of neutron source around a constant value, the source is time dependent. Accounting for the sinusoidal external neutron source, qðtÞ ¼ q 0 sinxt (x ¼ 2p T , where xðs À1 Þ and TðsÞ are angular frequency and period of source, respectively) can be useful for indicating the environment effects in neutron density [16]. Let us recall Eq. (6) for neutron density coefficients: Neutron density can be presented as: After imposing the initial condition, n 1 ð0Þ ¼ n 0 , one can determine the constant A from Eq. (15). So, we have: The solution of Eq. (8), in the second stage from insertion of reactivity, can be calculated from the integration factor method: In Eq. (16), B is the integration constant given in terms of the incomplete gamma functions, which can be obtained from continuity condition ðn 1 ðt 0 Þ ¼ n 2 ðt 0 ÞÞ: Equations (15) and (16) represent the complete solutions of Eqs. (15) and (8) in the presence of the sinusoidal external neutron source.

Results and discussions
Neutron external source play an important role during the start-up of a nuclear reactor. Therefore, the analytical solutions of neutron point kinetics equations in the presence of an external source are important in predicting the variation of neutron population during the start-up of a nuclear reactor. Here, the analyses are presented for the thermal reactor with the following parameters [9,12]: k ¼ 0:001s À1 , l ¼ 0:0015s, b ¼ 0:0075, r ¼ 0:0001s À1 , q 0 ¼ À0:006, and t 0 ¼ 10s. Equations (6) and (9) are the solutions of the neutron point kinetics equations with onegroup delayed neutron during the start-up of the thermal reactor in a two-stage insertion of reactivity. They are valid for any function that can be expressed by a power series. Due to fluctuation of the neutron source around a mean value, the source is actually time dependent. The general solutions of NPK equations can well describe neutron density response to any external neutron source, both quantitatively and qualitatively. To verify the validity of the analytical solutions, the variation of neutron density with constant and sinusoidal external neutron source was calculated. Specially, for n ¼ 0 and m ¼ 0; the neutron density was obtained in the presence of a constant source [see Eqs. (11), (12)]. This result is in line with Palma et al.'s work [12]. In this work, numerical calculations were carried out with generalized Runge-Kutta (GRK) method. The relative percentage errors (RPEs) of the neutron density using analytical and numerical solutions are defined as follows: ðRPEsÞ 12 ¼ j n 1 ðtÞ À n 2 ðtÞ n 1 ðtÞ j Â 100 ðRPEsÞ 13 ¼ j n 1 ðtÞ À n 3 ðtÞ n 1 ðtÞ j Â 100;  wherethe neutron densities of n 1 ; n 2 and n 3 are the solutions of numerical calculations without ISD, with ISD and analytical calculation with ISD, respectively. According to Table 1, the RPEs are smaller than 1 %, and the resultant RPEs are acceptable for engineering applications. So, the results of the analytical calculations are in a good agreement with the numerical results (see Fig. 1). Studying neutron density with conditions close to the critical state and small ramp rate reactivities ðr ( 1Þ is important for reactor control and safety. Effective multiplication factor ðk eff Þ values, close to the critical state is shown in Fig. 2.
As an example, Fig. 3 shows the numerical and analytical solution results for r ¼ 10 À5 s À1 ; q 0 ¼ À0:003 and k eff ¼ 0:9965, which are in good agreement with each other. Data in Table 2 show that the RPEs of the numerical and analytical solution are less than 0.5 % and also the RPEs decrease for k eff À! 1 and r ( 1. So, the ISD approximation is valid for neutronic calculations very close to criticality and r ( 1. For n ¼ 2k þ 1 and k ¼ 0; 1; 2; ::: the power series external source is: Such a phenomenon is also explained by Hetrick [16]. As a result, in the case of sinusoidal source with a small amplitude of oscillations, the constant source approximation is applicable at the start-up process [21,22]. As an example, the neutron density using analytical solution for n ¼ 20, m ¼ 0; ::9, T ¼ 50 s and q 0 ¼ 10 8 (blue dashed-dotted line) was compared with numerical results, with ISD and without ISD in Fig. 4. After raising the control rods and during sturt-up, 10 s t 200 s, the neutron density increases with time, as shown in Fig. 5. Regarding Table 3, the RPEs do not reach 1 %. Therefore, the results of numerical and analytical methods are the same.

Conclusion
Since the external neutron source term is very important at the start-up stage, it cannot be ignored. Due to the random nature of the external neutron source, we can take a timedependent source. Therefore, in this work, an analytical solution is presented for any neutron external source that can be found in terms of a power series. Analytical solutions in a general form are as in Eqs. (6) and (9). In particular, analytical solutions are obtained for the constant [12] and sinusoidal sources [16]. Numerical and analytical results in both cases are in a good agreement with each other.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.
Appendix A Using transformation of variable u ¼ t À k 2 and replacing in Eq. (5), we have: where g ¼ k 1 þ kk 2 À 1. Also using binomial expansion as follows: one can write: where constant dnm are defined as follows: Equation (19) can be solved through the use of the integrating factor method. Therefore, we get: Using the transformation of variable y ¼ Àku and a ¼ g þ n À m þ 1, we have: The solution of Eq. (21) can be written as: Appendix B Incomplete gamma function is defined as: Equation (9) can be solved through the use of the integrating factor method. Therefore we get: Using y ¼ Àst, one can write: