On the Zero-Neutron Density in Stochastic Nuclear Dynamics

: In this short paper, we compare the deterministic model for the nuclear reactor dynamic (Hetrick, 1993) with the stochastic model (Kinard and Allen, 2004). Our numerical results show coincidences between the deterministic model and the mean of the stochastic paths, although, as already observed by other authors, there is alarge amount of dispersion between the individual paths. Notably, we always observe that the neutron density approaches zero within a short time. In this paper, we investigate this question; more concretely, we study the mean-extinction of the neutron density. The technique used here ﬁrst builds the backward Kolmogorov differential equation and then solves it numerically using the ﬁnite-element method with F REE F EM ++. Our results conﬁrm that in a very short time the neutrons disappear although later they recover probably due to the external source.

In this paper, we research the time-dependency of a nuclear reactor. Previous models in nuclear reactor dynamics are investigated in [1,2]. We aim to determinate the time behavior of the power level of a nuclear reactor to obtain better control over it, for example, by using the control rod position. The deterministic models explained in [3] or [4] build onegroup, time-dependent diffusion equations. The literature that has studied their solutions and their numerical simulation is very extensive, as can be seen in the references of [4]. A few years later, J.G Hayes and E.J. Allen in [5] built a stochastic model, which is explained below. J.G Hayes and E.J. Allen consider the simplest case of a single precursor, while the determinist models assume m type of precursors; moreover, they assume that the nuclear reactor is large and homogeneous, so that spatial effects can be ignored. Then, they have two variables, n(t) and c(t), representing neutron and precursor densities, respectively.
In this way, considering a very small time-interval t, such that the probability of more than one event occurring during t is very small, they assume that the changes and their probabilities regarding the first order are given in Table 1 with x = ( n, c) T . Change Probability The physical parameters are as follows: • b the neutron birth rate due to fission; • d the neutron death rate due to fission; • ν the total number of neutrons per fission: (1 − β)ν prompt and βν delayed; • λ the constant of fission product c(t); • βν the number of atoms of fission product c(t) produced per fission; • q is the extraneous neutron sources.
Fixing x(t) at time t, we calculate the expected change for the change x = (n, c) T where and the covariance matrix such that Finally we obtain the following Stochastic Differential System (SDE) The deterministic nuclear model for two variables n(t) and c(t) is the linear Ordinary Differential Equations (EDOs) with a fix-point in whose stability depends on its Jacobin's eigenvalues. We can find these using the symbolic manipulator of MATLAB, although it has a complicated expression. Let us consider the two following examples.

Case 1:
As in the first numerical case, we consider the example in [6], p. 163, with which are values for a nuclear reactor involving the thermal fission of uranium-235 in [3]. For these values, the deterministic model (6) has a fixed point at 10 5 × (−0.0033, −3.3772) and its Jacobian has eigenvalues of 0.0480 and −48.1301, so the fixed point is unstable. In Figure 1, we plot the numerical solution with Matlab for n(0) = c(0) = 10 and n(0) = 250, c(0) = 10, 0 ≤ t ≤ 0.1. (We refer to [7][8][9][10] for a review of the numerical solutions of ODEs).

Case 2:
In the second example, we change the parameter ν = 2.432/2 to modify the asymptotic behavior of (6). When the fixed point is 10 3 × (0.0020, 1.0191) and the spectrum of the Jacobian is −5009.0 and −1.01, then the fixed point is an attractor of the linear system. In Figure 1, we plot the numerical solution for n(0) = 1, c(0) = 10 and n(0) = 3, c(0) = 10, 0 ≤ t ≤ 0.001. Our numerical simulations for stochastic model (5) were performed using the classic Euler-Maruyama numerical method, although it has a strong order 1/2 and weak order 1 (we refer to [11][12][13][14] and, more recently, Ref. [15] for a review on the numerical solutions of SDEs). Thus, we have where η 1 , η 2 ∼ N (0, 1). For each step, we have to compute the matrix D 1/2 . At this point, we have two options: it is well-known that a positive definite matrix has a unique positive definite matrix square root, and can be caluclated for a 2 × 2 matrix ([6], Remark 5.3). The other option uses the command sqrtm of MATLAB. The results are the same, although the second can be applied to any dimension. In Figure 2 we plot case 1. The deterministic solution is shown in red, the mean of 1000 trials is given in black, and tree sample paths are shown in green, blue and magenta. The step size is t = 10 −4 for n(0) = 10, c(0) = 10 on the left and n(0) = 250, c(0) = 10 on the right. These numerical results confirm the observations in Figure 5.7 in ( [6], p. 164), Figure 1 in [5] and Figure 1 in [13]: the means are close to the deterministic model, but with very important deviations; in other words, a sample path can be very far from that mean. One of the two variables could even reach very small values, close to zero. This same situation is observed in case 2, on which we do not believe it is necessary to comment. In the previous figures, we were surprised to see that the paths became very close to zero in a very short time. To understand this behavior, we researched a new random variable T defining the persistence time, i.e., the time it takes for the size of either variables to reach zero.
T ≡ inf{ t ≥ 0 : n(t) = 0 or c(t) = 0}, obviously, T depends on the initial values n(0)) and c(0), although this is not explicitly indicated. As the mean persistence-time τ ≡ E(T) for (5) satisfies the stationary backward Kolmogorov equation (see, for example, ( [6], p. 150), [16][17][18]), then, if we name n(0) = x and c(0) = y, this partial differential equation is and the boundary conditions are To resolve the boundary problems (9) and (10) using the FEM, let us multiply (9) by a regular function φ(x, y), satisfying the homogeneous Dirichlet boundary conditions. After integrating these over the domain Ω = [0, M x ] × [0, M y ], the following terms will appear: such that from (10) We computed the FEM solution of (9) and (10) with M x = 300, M y = 100 using FREEFEM++, with the result shown in Figure 3. From these results, we highlight the following: For this set of initial values: 0 < n(0) < 300, 0 < c(0) < 100, the mean time for the number of neutrons to reach zero is less than a tenth; more specifically, as we can see

Conclusions
Our main conclusion is that the deterministic solutions and mean of stochastic solution averages are very close, but, in a very short time, the neutrons disappear. However, they recover quickly, probably due to the external source. To test this assumption, we consider q = 100. In Figure 4, on the left, we plot n(0) = 10; c(0) = 10: the deterministic solution is given in red, the mean of 1000 trials in black, and tree sample paths in green, blue and magenta in 0 ≤ t ≤ 0.1. On the right, the numerical solution of (9) and (10) with M x = 100, M y = 100 is shown. These results suggest that when the source is reduced, the stochastic model moves away from the deterministic one and, above all, the neutrons disappear much earlier. However, in [20,21] an alternative stochastic formulation was proposed based on the classical Chemical Calgevin Method (e.g., [22][23][24][25]). However, in my opinion this statement is not correct. In remark 5.4 from ( [6], p. 143) or [26], it is shown that the both models are the same. Data Availability Statement: The numerical methods were implemented in MATLAB© and FREEFEM++ which is freely available [27]. The experiments were carried out in an Intel(R) Core(TM)i7-8665U CPU @ 1.90 GHz, 16.0 GB of RAM. The codes for the numerical tests are available on request.