Developed mathematical technique for fractional stochastic point kinetics model in nuclear reactor dynamics

NUCLEAR ENERGY SCIENCE AND ENGINEERING

Developed mathematical technique for fractional stochastic point kinetics model in nuclear reactor dynamics

Ahmed E. Aboanber，

Abdallah A. Nahla，

Adel M. Edress

Nuclear Science and Techniques

Vol.29, No.9

Article number 132

Published in print 01 Sep 2018

Available online 27 Jul 2018

DOI：10.1007/s41365-018-0467-0

35400

Fractional stochastic kinetics equations have proven to be valuable tools for the point reactor kinetics model, where the nuclear reactions are not fully described by deterministic relations. A fractional stochastic model for the point kinetics system with multi-group of precursors, including the effect of temperature feedback, has been developed and analyzed. A major mathematical and inflexible scheme to the point kinetics model is obtained by merging the fractional and stochastic technique. A novel split-step method including mathematical tools of the Laplace transforms, Mittage-Leffler function, eigenvalues of the coefficient matrix and its corresponding eigenvectors have been used for the fractional stochastic matrix differential equation. The validity of the proposed technique has been demonstrated via calculations of the mean and standard deviation of neutrons and precursor populations for various reactivity: step, ramp, sinusoidal, and temperature reactivity feedback. The results of the proposed method agree well with the conventional one of the deterministic point kinetics equations.

Itĉ stochastic point kinetics equationsTemperature feedback effectsWiener processFractional calculusMittage-Leffler function

1 Introduction

The dynamical processes described by the linear or nonlinear point kinetics equations are random processes in nature, that is due to the neutron population and precursor concentrations of delayed neutrons varying randomly with time. At the levels of high power, the stochastic manner is imperceptible, while at low-power levels, for example, at the start up of reactor operation, random fluctuation in the neutron population density and neutron precursor concentrations can be useful. The behavior variations of neutron population and precursor concentrations for nuclear reactors have been described by several innovators through employing stochastic models. Hayes and Allen [1] are the first authors who derived the stochastic model of the linear point reactor kinetics equations. They introduced a simplified stochastic model based on the Itô stochastic differential equations. The numerical results of this model using stochastic piecewise constant approximation (SPCA) have been compared with the Monte Carlo (MC) calculations and the experimental measurements [2]. Ha and Kim [3, 4] have presented the stochastic space-dependent kinetics model (SSKM) to solve the one dimensional monoenergetic space-time reactor kinetics. In 2012, Ray [5] developed Taylor 1.5 strong order methods and the Euler–Maruyama to solve the stochastic point kinetics equations with step reactivity, while in 2013, Ray and Patra [6] presented the same techniques with sinusoidal reactivity. The stochastic partial differential equation and stochastic difference equations have been presented by Ref. [7] for neutron transport equation. Furthermore, the power doubling time for a subcritical reactor is identified through the point kinetics system by Allen as a stochastic first-passage time problem [8]. In 2013, Ray and Patra [9] were the first authors who presented a numerical solution of fractional stochastic neutron point kinetic equations. Ayyoubzadeh and Vosoughi [10] simplified the system of Itô stochastic differential equations via alternative derivation of the stochastic differential equations. In 2016, Nahla and Edress [11] utilized the analytical exponential model (AEM) for the simplest formula of the stochastic point reactor kinetics system with various reactivity. They also proposed an efficient stochastic model (ESM) for the point kinetics model in Ref. [12]. da Silva et al. [13] presented a solution for the stochastic neutron point kinetics model. In 2017, Nahla [14] developed the analytical exponential technique (AET) to solve a stochastic nonlinear system of the point reactor kinetics equations with Newtonian temperature feedback reactivity. Finally, Singh and Ray [15] presented a comparison of two split-step methods for the numerical simulation of stochastic point kinetics equations in the presence of Newtonian temperature feedback effects.

The fundamental objective of this work is speculation of the stochastic point kinetics system to a fractional stochastic point kinetics system including multi-group of the precursor concentration. To overcome the difficulty arising from the merging of the fractional and stochastic techniques, a developed mathematical technique is presented for solving the equation in the matrix form of the proposed fractional stochastic model. This technique is based on split-step technique, Laplace transforms, the Mittage-Leffler function, eigenvalues of the coefficient matrix and its corresponding eigenvectors. The proposed method is applied to the fractional stochastic point kinetics system with various reactivity and different fractional order.

The paper is sorted out as follows: the preliminaries of the stochastic model, as well as the definitions of the Wiener process and fractional calculus, are introduced in Sect. 2. The solution of the fractional stochastic point kinetics equations with multi-group of a delayed precursor is derived in Sect. 3. The computational numerical results of the proposed system are discussed and compared with various stochastic techniques in Sect. 4. General conclusions including future work are presented in Sect. 5.

2 Preliminaries

In the following subsections, preliminaries of the Itô stochastic model, Wiener process, and fractional calculus are introduced briefly.

2.1 Itô stochastic model

Let us consider the following Itô stochastic differential model: [16]

d X (t) = g (t, X (t)) d t + f (t, X (t)) d W (t)

(1)

where, $g : ℝ^{+} \times ℝ^{m} \to ℝ^{m}$ and $f : ℝ^{+} \times ℝ^{m} \to ℝ^{m \times m}$ are locally bounded and measurable functions, W(t), is an m-dimensional Wiener process which is defined as the diffusion term.

The general solution of Eq. (1) can be written as:

\begin{array}{l} X (t) = X (0) & + \int_{0}^{t} g (u, X (u)) d u \\ + \int_{0}^{t} f (u, X (u)) d W (u) . \end{array}

(2)

Here, the integral $\int_{0}^{t} g (u, X (u)) d u$ is typically the Riemann-Lebesgue integral, while the integral $\int_{0}^{t} f (u, X (u)) d W (u)$ is considered as an Itô integral. Generally, the analytical solution of the multi-dimensional Eq.(2) to a great or significant extent is not possible, and numerical techniques are robustly utilized.

2.2 Wiener process

Recall that the standard Wiener process is a continuous-time stochastic process which is also called the standard Brownian motion. The Wiener process W(t) over [0,T] is a random variable, W(t), that depends on a continuous time, t, and is characterized by three conditions as follows [17]

1. W(t)=0 for t=0

2. W(t) has independent increments with $W (t) - W (s) \sim \sqrt{t - s} N (0,1)$ for $s \in [0, t]$ ,where 𝒩(0,1) is the normal distribution with zero mean and unit variance.

3. For 0 ≤ s₁ <t₁<s₂<t₂≤ T the increments W(t₁)-W(s₁) and W(t₂)-W(s₂) are independent random variables.

For numerical computational object, it is helpful to assume the discretized Brownian motion, where W(t) is specialized at discrete values, t. Accordingly, for some positive integer, N, consider h = T/N and let Wi denote W(ti) with ti=i h. The first condition states that W₀=0 with the probability equal to one and the other conditions state that Wi=Wi-1+dWi, i=1, 2,..., N, where each ∂Wi is an independent random variable of the form $\sqrt{h} N (0,1)$ .

2.3 Fractional calculus

The Riemann-Liouville fractional integral and Caputo derivative are defined respectively as follows:

Definition 1:[18] Let g(t) be a continues function, α≥ 0 and tє R. The Riemann-Liouville integral can be defined as:

I_{x}^{α} g (x) = \frac{1}{Γ (α)} \int_{0}^{x} {(x - t)}^{α - 1} g (t) d t,

(3)

where, Γ(α) is the gamma function of the fractional order.

Definition 2[19]: Let n be an integer number, where n-1<α<n and x>0. The Caputo fractional derivative of order α for a function g(x) is defined as

D_{x}^{α} g (x) = \frac{1}{Γ (n - α)} \int_{0}^{x} {(x - t)}^{n - α - 1} g^{(n)} (t) d t

(4)

where $g^{(n)} (t) = \frac{\partial^{n} g (t)}{\partial t^{n}}$ .

The relation between the Caputo fractional derivative [20] and the Riemann-Liouville fractional integral is introduced in the following formula:

I_{x}^{α} D_{x}^{α} g (x) = g (x) - \sum_{r = 0}^{n - 1} g^{(r)} {(0}_{+}) \frac{x^{r}}{r!}, n - 1 < α \leq n .

(5)

3 Developed mathematical technique

The fractional stochastic model of the point reactor kinetics equation can be written as: [1, 2, 21, 22]

D_{t}^{α} | P (t) 〉 = A | P (t) 〉 + | Q 〉 + B^{\frac{1}{2}} D_{t} | W (t) 〉,

(6)

where,

| P (t) 〉 = (\begin{matrix} n (t) \\ c_{1} (t) \\ c_{2} (t) \\ ⋮ \\ c_{I} (t) \end{matrix}), A = (\begin{matrix} \frac{ρ - β}{Λ} & λ_{1} & λ_{2} & \dots & λ_{I} \\ \frac{β_{1}}{Λ} & - λ_{1} & 0 & \dots & 0 \\ \frac{β_{2}}{Λ} & 0 & - λ_{2} & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋱ & ⋮ \\ \frac{β_{I}}{Λ} & 0 & \dots & 0 & - λ_{I} \end{matrix}), | Q 〉 = (\begin{matrix} q \\ 0 \\ 0 \\ ⋮ \\ 0 \end{matrix}),

(7)

t is the time, n(t) is the neutron population, ρ is the total reactivity, ci(t) is the i-group of delayed precursor concentration, λi is the decay constant of i-group of delayed neutrons, βi is the fraction of i-group delayed neutrons, $β = \sum_{i = 1}^{I} β_{i}$ is the total fraction of delayed neutrons, I is the total number of delayed neutron groups, Λ is the prompt neutron generation time, q is the external source of neutrons,

(8)

and $| Δ W (t) 〉 = \sqrt{h} | η 〉$ such that W₁(t), W₂(t),..., W_I+1(t) are the Wiener processes [5, 17].

Notice that the initial condition can be determined as n(0)=n₀, $\frac{d n (t)}{d t} |_{t = 0} = 0$ , and $\frac{d c_{i} (t)}{d t} |_{t = 0} = 0, i = 1,2, \dots, I$ . This means that

| P (0) 〉 = (\begin{matrix} n_{0} \\ \frac{β_{1} n_{0}}{Λ λ_{1}} \\ \frac{β_{2} n_{0}}{Λ λ_{2}} \\ ⋮ \\ \frac{β_{I} n_{0}}{Λ λ_{I}} \end{matrix}) .

(9)

Of course, when the variance matrix B=0 and the fractional derivative order α=1, Eq. (6) is reduced to the standard point reactor kinetics model.

In what follows, we aim to solve the fractional stochastic differential Eq. (6). The technique of split-step was utilized, that is, Eq. (6) is separated into deterministic and stochastic parts, followed by solving each of them separately. Applying the Laplace transformation on the deterministic part, i.e. consider B=0, of Eq. (6) is as follows: [23-25]

s^{α} | P (s) 〉 - s^{α - 1} | \bar{P} (0) 〉 - A | P (s) 〉 = \frac{1}{s} | Q 〉

(10)

where $| P (s) 〉 = L [| \bar{P} (t) 〉]$ , $| \bar{P} (t) 〉$ is the solution of the deterministic part for Eq. (6). Consequently, we get:

| P (s) 〉 = s^{α - 1} {[s^{α} I - A]}^{- 1} | \bar{P} (0) 〉 + \frac{1}{s} {[s^{α} I - A]}^{- 1} | Q 〉 .

(11)

Using the inverse Laplace transformation, we have:

| \bar{P} (t) 〉 = E_{α,1} (A t^{α}) | \bar{P} (0) 〉 + \int_{0}^{t} {(t - τ)}^{α - 1} E_{α, α} (A {(t - τ)}^{α}) d τ | Q 〉 .

(12)

Let us introduce the parameter z=t-τ into Eq.(12) to have:

| \bar{P} (t) 〉 = E_{α,1} (A t^{α}) | \bar{P} (0) 〉 + \int_{0}^{t} z^{α - 1} E_{α, α} (A z^{α}) d z | Q 〉 .

(13)

Using the integration property of the Mittag-Leffler function which is:

\int_{0}^{t} z^{b - 1} E_{a, b} (A z^{α}) d z = t^{b} E_{a, b + 1} (A t^{α}),

(14)

for b>0, then Eq. (13) reads as follows:

| \bar{P} (t) 〉 = E_{α,1} (A t^{α}) | \bar{P} (0) 〉 + t^{α} E_{α, α + 1} (A t^{α}) | Q 〉 .

(15)

Equation (15) introduces the general solution of the fractional stochastic point reactor kinetics equations, which depends on the stiff coefficient matrix A. To overcome the stiffness of this matrix, the coefficient matrix A was changed by its eigenvalues, ωj, [26-28] and the corresponding eigenvectors, |Vi⟩, of the matrix A [29-32]. Furthermore, over a small time interval with step size h, the matrices B and A are considered constant over the specified time interval [tm,t_m+1] where tm+1=tm+h and m=0,1,2,...,M-1. As a results of this substitution, the Mittag-Leffler function can be written as:

E_{a, b} (A z^{α}) = \sum_{j = 0}^{I} E_{a, b} (ω_{j} z^{α}) | V_{j} 〉 〈 U_{j} |,

(16)

where, the coefficient matrix A satisfies the following property in bra-ket space [33]

A | V_{j} 〉 = ω_{j} | V_{j} 〉, 〈 U_{j} | A^{T} = 〈 U_{j} | ω_{j},

(17)

and

〈 U_{l} | V_{j} 〉 = δ_{l, j} = {\begin{cases} 1, l = j \\ 0, l \neq j \end{cases} .

(18)

The ket eigenvectors, |Vj⟩, can be written by the following analytical form [33]

| V_{j} 〉 = σ_{j} (\begin{matrix} 1 \\ \frac{β_{1}}{Λ (ω_{j} + λ_{1})} \\ \frac{β_{2}}{Λ (ω_{j} + λ_{2})} \\ ⋮ \\ \frac{β_{I}}{Λ (ω_{j} + λ_{I})} \end{matrix})

(19)

and the bra eigenvectors, 〈Uj|, is

〈 U_{j} | = σ_{j} (1 \frac{λ_{1}}{(ω_{j} + λ_{1})} \frac{λ_{2}}{(ω_{j} + λ_{2})} \dots \frac{λ_{I}}{(ω_{j} + λ_{I})}) .

(20)

From the normalization ⟨Uj|Vj⟩=1, we can deduce:

σ_{j} = {(1 + \sum_{i = 1}^{I} \frac{β_{i} λ_{i}}{Λ {(ω_{j} + λ_{i})}^{2}})}^{- \frac{1}{2}}, \forall j = 0,1,2, \dots, I .

(21)

After introducing Eq. (16) into Eq. (15) we get:

\begin{array}{l} | \bar{P} (t_{m + 1}) 〉 = & \sum_{j = 0}^{I} [E_{α,1} (ω_{j} h^{α}) | V_{j} 〉 〈 U_{j} | \bar{P} (t_{m}) 〉 \\ + h^{α} E_{α, α + 1} (ω_{j} h^{α}) | V_{j} 〉 〈 U_{j} | | Q 〉 . \end{array}

(22)

The obtained results, $| \bar{P} (t_{m + 1}) 〉$ , from Eq. (22), were used to evaluate the variance matrix $\bar{B} (t_{m + 1})$ . Furthermore, according to the split-step method, the general solution of Eq. (6) is given by:

| P (t_{m + 1}) 〉 = | \bar{P} (t_{m + 1}) 〉 + \sqrt{h} {\bar{B}}^{\frac{1}{2}} (t_{m + 1}) | η 〉,

(23)

where $| Δ W (t) 〉 = \sqrt{h} | η 〉$ .

Equation (23) represents the general solution of the fractional stochastic point kinetics model with a multi-group of delayed precursor concentration.

4 Computational Results

In order to affirm the exactness and validity of the proposed technique over the traditional methods, the developed mathematical technique (DMT) for the fractional stochastic point kinetics model with a multi-group of delayed precursor concentration is tested through Matlab code. The mean and standard deviation of the neutron populations and precursor concentrations in different cases of reactivities, step, ramp, sinusoidal, and also in the presence of temperature feedback, are calculated. The results of the proposed method are compared to five different values of fractional order with the results of the traditional stochastic methods.

4.1 Step reactivity

The first computational example simulates two cases of step reactivity ρ=0.007 and ρ=0.003 for an actual reactor with six groups of delayed precursors [1]. The parameters of the two problems are taken from the following references [1, 5, 9, 11, 12] as follows: ν=2.5, Λ=0.00002 s, n₀=100 (neutrons), λi=[0.0127, 0.0317, 0.115, 0.311, 1.4, 3.87] s^-1, βi=[0.000266, 0.001491, 0.001316, 0.002849, 0.000896, 0.000182], β=0.007, and q=0. The mean and standard deviation at step size h=0.001 s after 500 trails are tabulated in Table 1 for reactivity ρ=0.007 and time T=0.001 s. In the calculation of the developed mathematical technique, the partial requests of the fractional order are taken respectively as: α=1, α=0.98, α=0.99, α=1.01, and α=1.02. The Matlab code is tested with the conventional results under the same conditions. The obtained results of the developed mathematical technique (DMT) are compared with the familiar Monte Carlo (MC) [1], stochastic piecewise constant approximation (SPCA) [1], Euler-Maruyama (EM) [5], Taylor 1.5 strong order (T1.5SO) [5], analytical exponential model (AEM) [12], efficient stochastic model (ESM) [11], and the deterministic point kinetics model (DPKM) [33].

Mean and standard deviation of the neutron population and the sum of the precursor population for step reactivity

	Method	Step size	α	E[n(t)]	σ[n(t)]	$E [\sum_{i = 1}^{6} c_{i} (t)]$	$σ [\sum_{i = 1}^{6} c_{i} (t)]$
ρ=0.007 t=0.001
	MC		1.0	135.67	93.376	4.464×10⁵	7.8073
	SPCA		1.0	134.55	91.242	4.464×10⁵	19.444
	EM		1.0	139.57	92.042	4.463×10⁵	6.071
	T1.5SO		1.0	139.57	92.047	4.463×10⁵	18.337
	AEM	h=0.00001	1.0	134.54	91.234	4.464×10⁵	19.235
	ESM	h=0.00001	1.0	134.96	6.8527	4.464×10⁵	2.529
	DPKM	h=0.000025	1.0	135.0		4.464×10⁵
	DMT	h=0.001	0.98	140.122221	7.872980	446362.78	3.794234
	DMT	h=0.001	0.99	137.268228	7.746689	446361.57	3.657256
	DMT	h=0.001	1.0	134.613600	7.627336	446360.52	3.525117
	DMT	h=0.001	1.01	132.144579	7.514619	446359.63	3.397625
	DMT	h=0.001	1.02	129.848341	7.408245	446358.87	3.274607
	DMT	h=0.0005	0.98	141.287614	7.427251	446362.79	3.262914
	DMT	h=0.0005	0.99	138.120575	7.318368	446361.49	3.137738
	DMT	h=0.0005	1.0	135.195181	7.216352	446360.39	3.017724
	DMT	h=0.0005	1.01	132.493180	7.120848	446359.45	2.902478
	DMT	h=0.0005	1.02	129.997682	7.031509	446358.66	2.791808
	DMT	h=0.0001	0.98	142.454816	7.029539	446363.37	2.946374
	DMT	h=0.0001	0.99	138.562403	6.927727	446361.78	2.796684
	DMT	h=0.0001	1.0	135.024891	6.833936	446360.47	2.656922
	DMT	h=0.0001	1.01	131.810129	6.747624	446359.39	2.525245
	DMT	h=0.0001	1.02	128.888862	6.668273	446358.50	2.401280
ρ=0.003 t=0.1
	MC		1.0	183.04	168.79	4.478×10⁵	1495.7
	SPCA		1.0	186.31	164.16	4.491×10⁵	1917.20
	EM		1.0	208.60	255.95	4.498×10⁵	1233.38
	T1.5SO		1.0	199.41	168.547	4.497×10⁵	1218.82
	AEM	h=0.001	1.0	186.30	164.14	4.490×10⁵	1911.91
	ESM	h=0.001	1.0	179.93	10.555	4.489×10⁵	94.75
	DPKM	h=0.0025	1.0	179.95		4.489×10⁵
	DMT	h=0.001	0.98	180.922922	12.847182	449304.10	94.953085
	DMT	h=0.001	0.99	180.466893	13.190199	449083.17	94.453581
	DMT	h=0.001	1.0	180.036327	13.547699	448878.15	93.938513
	DMT	h=0.001	1.01	179.629458	13.919605	448687.88	93.403821
	DMT	h=0.001	1.02	179.244639	14.305915	448511.28	92.845214
	DMT	h=0.0005	0.98	181.064354	11.889735	449347.96	95.974066
	DMT	h=0.0005	0.99	180.562722	12.307332	449103.51	95.397178
	DMT	h=0.0005	1.0	180.093120	12.742279	448878.30	94.806707
	DMT	h=0.0005	1.01	179.653284	13.194207	448670.78	94.198372
	DMT	h=0.0005	1.02	179.241129	13.662739	448479.56	93.567468
	DMT	h=0.0001	0.98	180.596300	10.464344	449450.09	88.949000
	DMT	h=0.0001	0.99	179.988863	10.964650	449149.09	88.353100
	DMT	h=0.0001	1.0	179.435711	11.490519	448876.38	87.735258
	DMT	h=0.0001	1.01	178.932909	12.042888	448629.26	87.089442
	DMT	h=0.0001	1.02	178.476566	12.622634	448405.28	86.409072

In a similar manner, the results of the second case ρ=0.003 are given in Table 1. The numerical results, in this case, are calculated at time t=0.1 s. Moreover, the pattern of the two individual neutron sample paths and the mean neutron population are shown in Fig. 1 at different values of the fractional order α=0.98, α=1, and α=1.02, symbolized a, b, and c respectively. The red solid curve is the mean neutron population, while the blue dashed curve and black dot-dashed curve represent the two individual neutron sample paths, which show approximately a real behavior for neutron density into the actual reactors.

Fig. 1.

Mean neutron population and two individual neutron sample paths for various values of fractional order and step reactivity ρ =0.003

4.2 Ramp reactivity

The second example simulates a ramp reactivity, where the same parameters are taken from the previous example. The function of reactivity is taken as: ρ=0.1β t and T=1 s. The result of a comparison of the proposed method (DMT), efficient stochastic model (ESM) [11], and the deterministic point kinetics model (DPKM) [33] is given in Table 2. To show the effect of the time interval, the code for the developed method was run at different time step as well as at different values of the parameter α. The relative percentage errors for E[n(t)] with the DMT is recorded as 0.025, 0.026, and 0.159 in compared with the conventional methods SPCA, AEM, and ESM respectively under the same conditions at time step h=0.01 s, Table 2 for ramp reactivity, while the relative percentage error is 0.05 for the DPKM at time step h=0.001 s. Furthermore, in the same table, the relative percentage errors for the DMT method at different step time interval (h=0.01 and h=0.001) are recorded as 0.126, while at h=0.01 and h=0.0001 it is 0.322. In Table 1 for step reactivity the effect of the time step interval appears from the relative percentage error at different time steps, e.g. at h=0.001 and h=0.0005 the RPE is -0.432, while at h=0.001 and h=0.0001 it is -0.306. The previous analysis confirms the stability of the DMT method and the effect of time step is acceptable. The behavior of the results show that an increase with the decreasing value of α and vise versa. Finally, we conclude that the validity of the proposed method shows a high agreement with the deterministic method DPKM as well as with the conventional methods (e.g. MC, SPCA, EM, T1.5SO, AEM, ESM, and DMT).

Mean and standard deviation of the neutron population and the sum of the precursor population for ramp reactivity ρ=0.1β t

Method	Step size	α	E[n(1)]	σ[n(1)]	$E [\sum_{i = 1}^{6} c_{i} (1)]$	$σ [\sum_{i = 1}^{6} c_{i} (1)]$
SPCA	h=0.01	1.0	113.268077	13.330142	448239.846	3009.93141
AEM	h=0.01	1.0	113.267707	13.327291	448239.798	3002.68282
ESM	h=0.01	1.0	113.116433	4.111150	448253.780	47.203115
DPKM	h=0.001	1.0	113.091124		448236.26
DMT	h=0.01	0.98	113.474956	9.476546	448434.483	46.871815
DMT	h=0.01	0.99	113.384203	9.473676	448331.161	46.840134
DMT	h=0.01	1.0	113.296719	9.470503	448232.846	46.811703
DMT	h=0.01	1.01	113.212358	9.468112	448139.299	46.782352
DMT	h=0.01	1.02	113.131034	9.466234	448050.294	46.753327
DMT	h=0.001	0.98	113.415094	3.971593	448535.015	43.751183
DMT	h=0.001	0.99	113.280768	4.058934	448378.858	43.695078
DMT	h=0.001	1.0	113.153337	4.150857	448233.774	43.632678
DMT	h=0.001	1.01	113.032287	4.247468	448098.976	43.555908
DMT	h=0.001	1.02	112.917106	4.348945	447973.736	43.348216
DMT	h=0.0001	0.98	113.305430	3.263194	448638.295	47.552005
DMT	h=0.0001	0.99	113.112671	3.411900	448425.709	47.488998
DMT	h=0.0001	1.0	112.932052	3.567299	448232.798	47.438442
DMT	h=0.0001	1.01	112.762814	3.729720	448057.756	47.366197
DMT	h=0.0001	1.02	112.604169	3.899621	447898.902	47.283241

For various values of fractional orders (α=0.98, α=1, α=1.02), the pattern of the two individual neutron sample paths and the mean neutron are shown in Figs. 2, 3, 4. Furthermore, the intensity of fluctuations for neutron sample paths increases with the mean neutron population. This phenomenon arises from the fact that the variance matrix including are stochastic part is dependent on the mean neutron population and precursors population.

Fig. 2.

Mean neutron population and two individual neutron sample paths for fractional order 0.98 and ramp reactivity ρ=0.1βt

Fig. 3.

Mean neutron population and two individual neutron sample paths for fractional order 1.0 and ramp reactivity ρ=0.1 βt

Fig. 4.

Mean neutron population and two individual neutron sample paths for fractional order 1.02 and ramp reactivity ρ=0.1β t

4.3 Sinusoidal reactivity

In the third example, the variation of reactivity insertion in the form of a sinusoidal change as $ρ = ρ_{0} \sin (\frac{π t}{T})$ [6] is considered. The parameters of the nuclear reactor with one group of delayed neutrons are as follows: ρ₀=0.005333 (0.68 $), β₁=β=0.0079, λ₁=0.077, Λ=10^-3, q=0, n₀=1 (neutrons) with period time of T=100 s. The two individual neutron sample paths and the mean neutron are shown in Figs. 5, 6, and 7 using different values of fractional orders (α=0.98, α=1, α=1.02) respectively, where the time step interval is h=0.01 s after 500 trails. We compared the stochastic and deterministic solutions to deduce the fact that the deterministic solution represents the average of the stochastic approach. Moreover, the individual sample paths oscillates around the deterministic curve. The intensity of fluctuations in these classes of figures arise from the variations of the mean neutron population.

Fig. 5.

Mean neutron population and two individual neutron sample paths for fractional order 0.98 and sinusoidal reactivity

ρ = 0.68 β \sin (\frac{π t}{50})

Fig. 6.

Mean neutron population and two individual neutron sample paths for fractional order 1.0 and sinusoidal reactivity

ρ = 0.68 β \sin (\frac{π t}{50})

Fig. 7.

Mean neutron population and two individual neutron sample paths for fractional order 1.02 and sinusoidal reactivity

ρ = 0.68 β \sin (\frac{π t}{50})

4.4 Temperature feedback reactivity

In the most nuclear literature, there are two cases for the external reactivity, step, and ramp external reactivities. The means of the neutron population are calculated for a U²³⁵ nuclear reactor with step and ramp external reactivities. In what follows, the effect of Newtonian temperature feedback introduced into the reactivity is analyzed. The new reactivity form in the presence of temperature feedback is given by:

ρ (t) = ρ_{ex} (t) - σ \int_{0}^{t} N (τ) d τ

(24)

where, ρ_ex(t) represents the external reactivity and σ is the nonlinear coefficient part which represents the product of the reciprocal of the thermal capacity and the temperature coefficient.

4.4.1 Step external reactivity

To check the developed mathematical technique for the nonlinear fractional stochastic model in the presence of temperature feedback and step external reactivity, let us take the parameters of the U²³⁵ nuclear reactor as follows [14, 15, 34]: λi=[0.0124, 0.0305, 0.111, 0.301, 1.13, 3.0] s^-1, βi=[0.00021, 0.00141, 0.00127, 0.00255, 0.00074, 0.00027, 0.00645], β=0.00645, Λ=5.0× 10^-5 s, σ=2.5× 10^-6 MWs)^-1, and N(0)=1 (neutrons).

In Table 3, the peak of the mean neutron population with the corresponding time at various step external reactivity, ρ_ex=1.0$, ρ_ex=1.5$, and ρ_ex=2$, are tabulated for different values of fractional order 0.98, 0.99, 1.0, 1.01, and 1.02. Using time step h=0.001 s and after 500 trails, the peak of the mean neutron population is compared with the peak of the mean neutron population using the analytical exponential technique (AET) [14]. In addition, the mean neutron population and the two individual neutron sample paths are plotted in Figs. 8, 9, and 10 for the step external reactivity (ρ_ex=0.5$) and Figs. 11, 12, and 13 for the step external reactivity (ρ_ex=0.75$) using different values of fractional derivative order 0.98,1.0, 1.02, respectively. The most important notice is that the mean neutron population increases with time until it reaches the maximum value due to the positive external reactivity. Thereafter, the mean neutron tends to zero with the increasing time due to the effect of temperature reactivity feedback. Again, these figures confirm that the amplitude of fluctuations for neutron sample paths is affected by the direct variations of the mean neutron population.

Peak of the mean neutron population and its time for the nonlinear fractional stochastic model at step external reactivity

Method	α	ρ_ex($)	0.5	0.75	1.0	1.5	2.0
DMT	0.98	Peak	53.0279	189.727	895.547	41083.125	158278.694
		Time	25.378	7.732	0.942	0.151	0.088
	0.99	Peak	49.565	176.089	826.687	36557.357	140582.052
		Time	26.284	8.181	1.025	0.162	0.094
	1.0	Peak	46.289	163.986	768.739	34289.333	132756.323
		Time	28.142	8.762	1.09	0.173	0.101
	1.01	Peak	43.137	152.776	712.311	30306.984	117251.348
		Time	30.303	9.429	1.211	0.186	0.108
	1.02	Peak	40.239	142.525	662.795	28353.797	111465.515
		Time	30.708	9.711	1.245	0.199	0.116
AET	1.0	Peak	46.235334	164.209516	770.5188	33119.58	128083.1
		Time	28.142	8.867	1.040	0.174	0.101

Fig. 8.

Mean neutron population and two individual neutron sample paths for fractional order 0.98 and temperature feedback reactivity ρ_ex=0.5$ and σ=2.5×10^-6

Fig. 9.

Mean neutron population and two individual neutron sample paths for fractional order 1.0 and temperature feedback reactivity ρ_ex=0.5$ and σ=2.5× 10^-6

Fig. 10.

Mean neutron population and two individual neutron sample paths for fractional order 1.02 and temperature feedback reactivity ρ_ex=0.5$ and σ=2.5×10^-6

Fig. 11.

Mean neutron population and two individual neutron sample paths for fractional order 0.98 and temperature feedback reactivity ρ_ex=0.75$ and σ=2.5×10^-6

Fig. 12.

Mean neutron population and two individual neutron sample paths for fractional order 1.0 and temperature feedback reactivity ρ_ex=0.75$ and σ=2.5× 10^-6

Fig. 13.

Mean neutron population and two individual neutron sample paths for fractional order 1.02 and temperature feedback reactivity ρ_ex=0.75$ and σ=2.5×10^-6

4.4.2 Ramp external reactivity

In the test example for ramp external reactivity, the parameters for the U²³⁵ nuclear reactor is taken as the same values from the nonlinear case, where Λ=10×10^-5 s, σ=10^-11, or 10^-13 (MW s)^-1 and the external reactivity is ramp (ρ_ex= 0.01 t, 0.1 t). The peak of the mean neutron population and its time with various values of fractional order (0.98,1.0, 1.02) are compared with the results of the analytical exponential technique (AET) [14] in Table 4. The results are calculated with time step h=0.001 s and the number of trails is 500. In addition, two individual neutron sample paths and the mean neutron population are drawn for external reactivity ρ_ex(t)=0.003t and the nonlinear coefficient σ=10^-11 in Figs. 14, 15, and 16 and σ=10^-13 in Figs. 17, 18, and 19. In Figs. 14, 15, 16, 17, 18 and 19, the mean neutron increases with time until it reaches the peak due to the external reactivity increasing. After that the mean neutron decreases due to the effect of temperature feedback. Therefore, the mean neutron population is almost stable due to the effect of external reactivity, which is equivalent with the effect of temperature feedback. Moreover, the fluctuations of neutron sample paths disappear approximately around the sharp peak due to the logarithmic scale and a large increase in the neutron population at a very small time. Figure 14 shows a pattern of the effect of the neutron population and fluctuations at different sections with time. A sharp increase of fluctuation with a slight increase of the neutron population in the ramp section is followed by a slight increase of the intensity of fluctuations compared with a sharp increase of the neutron population. Moreover, the same remark is observed at the peak and the remainder section of the figures.

Peak of the mean neutron population and its time for the nonlinear fractional stochastic model at reactivity

ρ (t) = a t - σ \int_{0}^{t} N (τ) d τ

α			a=0.003		a=0.01		a=0.1
			σ=10^-11	σ=10^-13	σ=10^-11	σ=10^-13	σ=10^-11	σ=10^-13
DMT
	0.98	Peak	4.409101×10⁹	5.771654×10¹¹	1.701579×10¹⁰	2.144910×10¹²	1.876758×10¹¹	2.242436×10¹³
		Time	2.446	2.488	0.831	0.852	0.130	0.136
	0.99	Peak	4.452107×10⁹	5.813629×10¹¹	1.704601×10¹⁰	2.124988×10¹²	1.918749×10¹¹	2.290324×10¹³
		Time	2.461	2.504	0.839	0.861	0.132	0.138
	1.0	Peak	4.477012×10⁹	5.828427×10¹¹	1.689973×10¹⁰	2.102489×10¹²	1.849657×10¹¹	2.146883×10¹³
		Time	2.475	2.520	0.847	0.869	0.135	0.142
	1.01	Peak	4.505250×10⁹	5.847340×10¹¹	1.673285×10¹⁰	2.081893×10¹²	1.860123×10¹¹	2.233038×10¹³
		Time	2.491	2.537	0.856	0.879	0.138	0.144
	1.02	Peak	4.514868×10⁹	5.844036×10¹¹	1.656222×10¹⁰	2.059339×10¹²	1.786014×10¹¹	2.131587×10¹³
		Time	2.507	2.554	0.864	0.888	0.141	0.148
AET
	1.0	Peak	4.482853×10⁹	5.833649×10¹¹	1.687565×10¹⁰	2.099857×10¹²	1.853996×10¹¹	2.217193×10¹³
		Time	2.476	2.521	0.847	0.870	0.135	0.141

Fig. 14.

Mean neutron population and two individual neutron sample paths for fractional order 0.98, temperature feedback reactivity ρ_ex=0.003 t and σ=10^-11

Fig. 15.

Mean neutron population and two individual neutron sample paths for fractional order 1.0, temperature feedback reactivity ρ_ex=0.003 t and σ=10^-11

Fig. 16.

Mean neutron population and two individual neutron sample paths for fractional order 1.02, temperature feedback reactivity ρ_ex=0.003 t and σ=10^-11

Fig. 17.

Mean neutron population and two individual neutron sample paths for fractional order 0.98, temperature feedback reactivity ρ_ex=0.003 t and σ=10^-13

Fig. 18.

Mean neutron population and two individual neutron sample paths for fractional order 1.0, temperature feedback reactivity ρ_ex=0.003 t and σ=10^-13

Fig. 19.

Mean neutron population and two individual neutron sample paths for fractional order 1.02, temperature feedback reactivity ρ_ex=0.003 t and σ=10^-13

General description and analysis for the results of the developed method can be summarized as follows. Figures 1, 2, 3, 4, 5, 6, and 7 show a large class of changes consists of small variations of the cross sections, positive or negative changes of reactivity attributed to the neutron population, around an expected value E[n(t)], corresponding to a critical case in which the perturbation can be induced flux fluctuations, where: δ n(t)=n(t)-E[n(t)] and δρ(t)=ρ(t)-E[ρ(t)].

Strong-frequency power fluctuations caused by increasing the reactivity (the fission rate is sufficiently large for the average value equations) of the reactor, as shown in Figs. 1, 2, 3, and 4. In the opposite, low-frequency power fluctuations can arise from decreasing the reactivity (i.e. the fission rate) as shown in Figs. 5, 6, 7, 8, 9, 10, 11, 12, and 13. Furthermore, Figs. 14, 15, 16, 17, 18, and 19 deal with the variation of the neutron population, which obtained by studying the fluctuations arise from the temperature feedback. These figures conclude that the randomness of the input is communicated to the output via the response characteristics of the system, where the fluctuation is above or below the mean value at instant time.

5 Conclusion

The developed mathematical technique was presented for a linear/nonlinear fractional stochastic model of the point reactor kinetics system with a multi-group of delayed neutron precursors. This system is characterized by its stochastic behavior and can offer the only average or mean values of the modeled populations. However, the neutron population and the delayed neutron precursor concentrations vary randomly with time, meaning the real dynamical process is stochastic. This system was numerically implemented using a stochastic piecewise constant approximation (SPCA) due to the stiffness of these equations. In this paper, the matrix formula for this fractional stochastic model is solved through a developed mathematical technique, which is based on the split-step method, Laplace transforms, Mittage-Leffler function, eigenvalues, and eigenvectors. The mean and standard deviation of the neutron population and the sum of the precursor population were calculated for step, ramp, sinusoidal, and the temperature feedback reactivity insertion which represents the nonlinear fractional stochastic model. Moreover, this fractional differential system was calculated with different values of the fractional derivative order. In order to validate the proposed method (DMT), we present a comparison with the conventional results in the literature of the stochastic model and the deterministic point kinetics model, showing that the method is in agreement with those already established. The future work will be included the derivation and the study of a fractional stochastic model for the time-space kinetics equations.

References

[1]

J.G. Hayes, E.J. Allen,

Stochastic point kinetics equations in nuclear reactor dynamics

, Ann. Nucl. Energ. 32, 572-587 (2005). doi: 10.1016/j.anucene.2004.11.009