1 Introduction

The point reactor kinetics equations are the most essential model in the field of nuclear engineering. This system is a coupled linear differential equation and describes the neutron population density and the precursor concentration of delayed neutrons at the center point of homogenous reactors. A point reactor is a homogenous reactor in which the spatial effects have been eliminated. This is obviously possible if the homogenous reactor’s length is infinite in all spatial dimensions.

The dynamical process explained by the point kinetics equation is stochastic in nature. The neutron population density and delayed neutron precursor concentrations differ randomly with respect to time. At the levels of high power, the random behavior is imperceptible. But at low-power levels, such as at the beginning, random fluctuation in the neutron population density and neutron precursor concentrations can be crucial. The aim of this work presents an accurate method for stochastic point kinetics equations with step, ramp, and sinusoidal reactivities.

There are some techniques which are used for stochastic point reactor kinetics equations. The first of these techniques is the stochastic Piecewise Constant Approximation (stochastic PCA) method and Monte Carlo computations, which were used to calculate the neutron population density and sum of the precursors concentration population density for different values of step reactivity [1, 2]. A simplified stochastic model based on the forward stochastic model in the stochastic kinetics theory and the Itô stochastic differential equations was developed for treating monoenergetic space-time nuclear reactor kinetics in one dimension [3, 4]. Simulation and experimental study of a random neutron analyzing system with a ²⁵²Cf neutron source was presented in Ref. [5]. The Euler-Maruyama and Taylor 1.5 strong order methods were presented for solving stochastic point kinetics equations with step and sinusoidal reactivities [6, 7]. Finally, simplified stochastic point kinetics equations (SSPK) were modeled with a system of Itô stochastic differential equations. This approach does not require computing the square root of a matrix, which is a great computational advantage [8].

In this work, the system of stochastic point kinetics model in the dynamic nuclear reactor is derived in Sec. 2. The analytical exponential model is presented and applied to solve the stochastic point kinetics equations in Sec. 3. The numerical results of the proposed method are discussed and compared with the traditional methods in Sec. 4. Finally, the conclusions and future work are discussed in Sec. 5.

2 Stochastic Model

The one speed neutron diffusion equations with a multi-group of delayed precursor concentrations are written as Ref. [9-11]:

where Φ(r,t) is the neutron population, r is the position (cm), t is the time (s), Ci(r,t) is the i-group of delayed precursor concentration, S₀(r,t) is the external neutron source, D is the diffusion coefficient, ∑a is the absorption cross sections, ∑f is the fission cross sections, ν is the neutron fission, v is the neutron speed, λi is the decay constant of i-group of delayed neutrons, βi is the fraction of i-group delayed neutrons, and $\beta ={\displaystyle \sum _{i=1}^I}{\beta }_i$ is the total fraction of delayed neutrons.

The neutron population and the i-group of the delayed precursor concentration, using separation of variables, can be written as:

where n(t) is the neutron population density as a function of time only and ci(t) is the precursor concentration density of delayed neutrons.

The function $\widehat{\Psi }(r)$ is the fundamental function, which can be determined from the following

where $B_g^2$ is the geometric buckling.

Substituting Eqs. (3) and (4) into Eqs. (1) and (2) we get

According to Ref. [1, 2, 11], Eqs. (5) and (6) can be separated into four terms: deaths, births, decay, and external sources. Therefore,

Equations (7) and (8) can be rewritten as deterministic point kinetics equations using the reactivity $\rho =1-\frac{{\sum }_a+D{B}_g^2}{\nu {\sum }_f}$, which is a function of time, the generation time of neutrons $\Lambda =\frac{1}{\nu {\sum }_{f}v}$, and the constant $\alpha =\frac{1}{\nu }$.

Note that, n(t) is the population size of neutrons. ci(t) is the population size of the i^th group of delayed precursor concentration. Equations (9) and (10) are the deterministic point kinetics equations where terms are separated into births, deaths, decay, and external source q. Indeed, the neutron birth rate due to fission is $b=\frac{1-\alpha -\beta }{\Lambda [(1-\beta )\nu -1]}=\frac{1}{\nu \Lambda }$, where [(1-β)ν-1] is the number of new born neutrons in each fission. The neutron death rate due to captures and leakage is $d=\frac{1-\rho -\alpha }{\Lambda }$. Also, λici is the rate of the i^th group of delayed precursor decay.

To derive the stochastic dynamical system, let us take a small time interval size h where the probability of more than one event occurring during this time interval is small. During time interval h, there are (I+3) different possibilities for an event. The change in populations n and ci during this time interval is

In the present investigation let us assume that the changes are approximately normally distributed. The (I+3) possibilities and its probabilities are Ref. [1, 2]:

In the present analysis, it is assumed that the extraneous source produces neutrons randomly following a Poisson process with intensity q. According to these assumptions, the changes in neutron population and precursor concentration are approximately normally distributed with a mean of

According to Ref. [8], the variance takes the form

where B is the diagonal matrix as

Using central limit theorem, the random variate $\frac{|\Delta P\rangle -E(|\Delta P\rangle )}{\sqrt{Var(|\Delta P\rangle )}}=|\eta \rangle$ follows standard normal distribution [1, 2, 12]. This implies

where $|\eta \rangle =\left(\begin{array}{c}\eta 1\\ \eta 2\\ ⋮\\ {\eta }_{I+1}\end{array}\right)$, and η₁,η₂, η_I+1 N(0,1).

Substituting Eqs. (18) and (19) into Eq. (20), we have

(21)

This Eq. (21) gives, as h→0, the following Itô stochastic differential equation system

where,

(23)

and $|\Delta W(t)\rangle =\sqrt{h}|\eta \rangle$ where, W₁(t), W₂(t),, W_I+1(t) are Wiener processes [6, 7].

Equation (22) represents the stochastic point reactor kinetics model. Since B=0, it reduces to the deterministic point kinetics model. This model was solved analytically and numerically in many references, for example, Ref. [13-18].

3 Analytical Exponential Model

Using the integration factor, differential Eq. (22) in matrix form is rewritten as

Let us divide the time into M, very small time intervals with step size h. The matrices A and B are constants during the time interval [tm,t_m+1] where t_m+1=tm+h, and m=0,1,2,⋅s,M-1. Equation (24) can then take the following form

and consequently, we get

The mathematical treatment of this system can be found by calculating all the eigenvalues and corresponding eigenvectors of matrix A and performing straightforward computations. However, this is an expensive scheme, especially when the reactivity varies with time, since the eigenvalues of matrix A are calculated by solving the Inhour equation, a (I+1)th-order algebraic equation, at each time step.

The eigenvectors of A, denoted by ket vectors |U_j⟩ and the corresponding eigenvalues denoted by ωj, obey the relation

The eigenvalues ωj of matrix A are the roots of the inhour formula

For j, an arbitrary function f(hA) satisfies the following expression [13-16]

and, consequently, we get

with the properties

The ket eigenvectors |U_j⟩ and the bra eigenvectors U_j| are calculated analytically [13-16]

For the normalized condition ⟨U_j|U_j⟩=1, we get ${\Omega }_j=\frac{1}{\sqrt{1+{\displaystyle \sum _{i=1}^I}\frac{{\beta }_i{\lambda }_i}{\Lambda {({\omega }_j+{\lambda }_i)}^2}}}\mathrm{,}=\mathrm{0,1,2,}\cdots \mathrm{,}I$.

Hence, the general solution of the stochastic point kinetics equation takes the following form

where the initial condition is $|P(0)\rangle =|P(t_0)\rangle ={\left(n_0\frac{{\beta }_1{n}_0}{\Lambda {\lambda }_1}\frac{{\beta }_2{n}_0}{\Lambda {\lambda }_2}\cdots \frac{{\beta }_I{n}_0}{\Lambda {\lambda }_I}\right)}^T$.

4 Numerical Results and Discussion

In this section, many of the examples are presented to measure the accuracy of the analytical exponential model (AEM) for stochastic point kinetics equations. The mean and standard deviation of the neutron and precursor population are calculated by solving the stochastic point kinetics equations with three different cases: step, ramp, and sinusoidal reactivities, and are compared with the published stochastic methods.

The first benchmark problem does not model an actual physical nuclear reactor problem, but this problem provides a simple computational solution for comparing the stochastic model. This model simulates a step reactivity insertion and assumes one neutron precursor [1]. The parameters of this benchmark are: the neutron generation time $\Lambda =\frac{2}{3}$ s, reactivity $\rho =-\frac{1}{3}$, decay constant λ₁=0.1 s^-1, fraction delayed neutron β₁=β=0.05, number of neutrons per fission ν=2.5, external source q=200 s^-1, and the initial condition assumes n(0)=400 and c(0)=300. Table 1 presents the mean and standard deviations of neutron and precursor populations using 40 time intervals for a time interval of length of 2. The results of the AEM are compared with the results of the Monte Carlo, stochastic PCA method [1], Euler-Maruyama [6], Taylor 1.5 strong order [6] and the deterministic point kinetics model (DPKM) i.e B=0 at time t=2 s. All the stochastic methods use 5000 trails. Table 1 shows the accuracy of the analytical exponential model such that the results of the AEM are in good agreement with the deterministic point kinetics model, more than other methods under the same conditions. This example confirms that the AEM is an efficient method for solving stochastic point kinetics equations in nuclear reactor dynamics.

The second benchmark problem simulates a step reactivity insertion for an actual nuclear reactor with six groups of delayed neutrons [1]. The parameters of this reactor are as follows: βi=[0.000266, 0.001491, 0.001316, 0.002849, 0.000896, 0.000182], β=0.007, λi=[0.0127, 0.0317, 0.115, 0.311, 1.4, 3.87] s^-1, Λ=0.00002 s, ν=2.5, n₀=100, and no external source q=0 s^-1. Two cases of step reactivity insertion, ρ=0.003 and ρ=0.007, are presented. Calculation results of the mean and standard deviations can be seen in Table 2 at time 0.1 s for the first case ρ=0.003 and at time 0.001 s for the second case ρ=0.007. The results of the AEM are compared with the results of Monte Carlo, Stochastic PCA method [1], Euler-Maruyama [6], Taylor 1.5 Strong order [6], simplified stochastic point kinetics (SSPK) method [8], and DPKM (B=0). All the stochastic methods use 5000 trails. The agreement is seen between the results of the AEM, deterministic point kinetics model, Monte Carlo, Stochastic PCA, and SSPK methods in Table 2 for two cases of reactivity.

The mean neutron population and two individual neutron sample paths are given in Fig. 1 for a step reactivity insertion (ρ=0.003) and in Fig. 2 for a step reactivity insertion (ρ=0.007).

Fig. 1.Full-screen View

(Color online) Mean neutron population and two individual neutron sample paths for reactivity ρ=0.003.

Fig. 2.Full-screen View

(Color online) Mean neutron population and two individual neutron sample paths for reactivity ρ=0.007.

The third benchmark problem simulates a ramp reactivity insertion for the above actual nuclear reactor with six groups of delayed neutrons. The parameters of this reactor are the same as in the second benchmark problem, except the reactivity is ρ=0.1βt. The mean neutron population and two individual neutron sample paths are shown in Fig. 3 using a time interval size of h=0.01 s after 5000 trails.

Fig. 3.Full-screen View

(Color online) Mean neutron population and two individual neutron sample paths for a ramp reactivity.

The final benchmark problem simulates a sinusoidal reactivity insertion, $\rho ={\rho }_0\sin \left(\frac{\pi t}{T}\right)$ [7]. The parameters of the reactor with one group of delayed neutrons are as follows: ρ₀=0.005333 (0.68 $), λ₁=0.077 s^-1, β₁=β=0.0079, Λ=10^-3 s, q=0 s^-1, n₀=1, and a half period time of T=50 s. The mean neutron population and two individual neutron sample paths are seen in Fig. 4 using the time interval size h=0.1 s after 5000 trails.

Fig. 4.Full-screen View

(Color online) Mean neutron population and two individual neutron sample paths for a sinusoidal reactivity.

5 Conclusion

In this work, stochastic point kinetics equations were introduced, which represent the generalization of the deterministic point kinetics equations. The analytical exponential model was described, which is based on the eigenvalues and corresponding eigenvectors. This method is an efficient approximate for the stiff system of stochastic point kinetics differential equations. The mean and standard deviations of the neutron and precursor populations using the analytical exponential model are in good agreement with these by the deterministic point kinetics equations, Monte Carlo, and stochastic PCA methods more than other references methods. This agreement confirms that the analytical exponential model is an efficient method for solving stochastic point kinetics equations in nuclear reactor dynamics using a step, ramp, and sinusoidal reactivities.

Possible future work may include the derivation and study of stochastic multi-energy group reactor kinetics equations where spatial effects are included.