Higher-order approximate solutions of fractional stochastic point kinetics equations in nuclear reactor dynamics

S.Singh·S.Saha Ray

Abstract Stochastic point kinetics equations(SPKEs)are a system of Itoˆstochastic differential equations whose solution has been obtained by higher-order approximation.In this study,a fractional model of SPKEs has been analyzed.The efficiency of the proposed higher-order approximation scheme has been discussed in the results section.The solutions of SPKEs in the presence of Newtonian temperature feedback have also been provided to further discuss the physical behavior of the fractional model.

Keywords Fractional stochastic point reactor kinetics equations·Fractional calculus·Higher-order approximation·Caputo derivative·Neutron population

1 Introduction

Stochastic point kinetics equations(SPKEs)are a system of coupled nonlinear stochastic differential equations(SDEs)[1]and are important in nuclear engineering.The SPKEs model a system of ItoˆSDEs,specifically,neutron population and delayed neutron precursors.The physical dynamical system has been established to be a population process,and techniques have been employed by Hayes and Allen[2]to transform deterministic point kinetics equations into a system of SDEs.The fractional diffusion model is normally applied to large variations of neutron cross sections,which preclude the use of classical neutron diffusion equations[3–6].Various approximation methods have been developed to improve the control of processes in the nuclear reactor.

In recent developments,dynamic,multiphysics phenomena face many challenges regarding accurate numerical schemes, resulting in severe computational requirements.An approach to reduce the severe computational requirements of standard low-order simulations is to employ higher-order formulations.In the hierarchy of highorder methods,compact schemes represent an attractive choice for reducing dispersion and anisotropy errors.

Nowak et al.[7]presented the numerical solutions of a fractional neutron point kinetics model for a nuclear reactor.Numerical solutions of SPKEs by implementing stochastic piecewise continuous approximation(PCA)have been obtained by Hayes and Allen[2],which provided a very succinct idea about the randomness of neutron density and precursor concentrations.Saha Ray[8]showed that Euler–Maruyama(EM)and Taylor 1.5 strong order numerical schemes are well-founded estimators compared with stochastic PCA.Saha Ray and Patra[9]applied the Grunwald–Letnikov definition of fractional derivative for solving SPKEs.Nahla and Edress[10]showed the efficiency of analytical exponential model(AEM)for obtaining solutions of SPKEs.

In our present work,we demonstrate the fractional numerical method for obtaining the mean neutron population for various reactivities.In Sect.2,we introduce the fractional stochastic nonlinear point reactor kinetics equations(SNPKE).The fractional ItoˆSDE for NPKEs is obtained using the central limit theorem.In Sect.3,we discuss the higher-order approximation method for Caputo derivative,and in Sect.4,we discuss its application.In Sect.5,we discuss the obtained numerical solutions for various reactivities.

2 Framework of fractional SNPKE

We now introduce an elementary definition before stating the central limit theorem.

De finition 2.1 Let¯XMbe the sample mean of M independent samples X1，...，XMof a random variable X such that

where μ is the mean and σ2is the variance for each independent and identically distributed real-valued random variables Xj.

Theorem 2.1 (Central Limit Theorem)If each random variable Xjhas a finite second moment with Var Xj=σ2,the distribution ofconverges to that of a Gaussian random variable with mean μ and variance.Thus,converges in distribution to Z ～N（0，σ2）.

The fractional ItoˆSDE for NPKEs with temperature feedback effects[11–15]obtained from the central limit theorem can be written as follows[16]:

where α is the fractional order of the derivative and 0＜α≤1,

The coefficient matrix A（t）is represented in the fol-

where ρ is the total reactivity,is the total fraction of delayed neutrons,βiis the fraction and λiis the decay constant of i-group of delayed neutrons,and l is the prompt neutron generation time.

The covariance matrix B（t）,which is evaluated in Ref.[16],is represented as follows:

3 Higher-order approximation scheme

3.1 Fractional calculus

The Caputo derivative operator for α ∈ （0，1）is defined as follows:

in which Γ（·）is the Euler gamma function.

In this paper,a （3- α）thorder scheme for Caputo derivativewith α ∈ （0，1）has been discussed for obtaining the solutions of fractional SNPKEs.Let 0=t0＜t1＜···＜tM=T and tn=t0+nτ be the equidistant discretized times withfor M∈Z such that τ∈ （0，1）.Now,using the Taylor expansion to Ψ′（s）,Ψ（ti-1）,and Ψ（ti+1）at point t=ti（0 ≤ i＜n）,we get

Hence,we can obtain the following:

Therefore,the Caputo derivative can be discretized as

where n=1，2，...，M and τnis the truncation error.

Therefore,the Caputo derivative has the following numerical approximation[17]:

The right-hand side of Eq.(3.2)can be expressed as follows:

where

Therefore,

Let

It can be proven that|S（n）|is bounded for n≥1[17–19].This proves that the seriesconverges.

On further simplification,we get:

In Eq.(3.2),if i=0,then Ψi-1= Ψ-1,which lies outside of ［0，T］.Various options have been provided to approach Ψ-1.In numerical calculations,the neighboring function values have been used to approximate Ψ-1,that is,

1.When Ψ′（0）= Ψ′′（0）=0,then Ψ-1= Ψ0+O（τ3）;the convergence order is O（τ3-α）.

2.When Ψ′（0）=0, Ψ′′（0）/=0, then Ψ-1= Ψ0the convergence order of Eq.(3.1)is O（τ2）.

3.When Ψ′（0）/=0,then the convergence order is O（τ）.

4 Solution of SPKE by higher-order approximation method

In this section,higher-order approximation to Caputo derivative has been applied to Eq.(2.1)as follows:

Therefore,the above expression can be simplified into an explicit numerical scheme as follows:

Table 1 Mean peak values of N（t）for different step reactivities

Fig.1(Color online)Mean N（t）and two arbitrary sample paths for a step reactivity ρ =0.003 and α =0.96,b step reactivity ρ =0.003 and α =0.98,c step reactivity ρ =0.003 and α =0.99

Table 2 Mean peak values of N（t）for ramp reactivity ρ =0.1βt and different values of fractional order α

Fig.2(Color online)Mean N（t）and two arbitrary sample paths for a ramp reactivity ρ =0.1βt and α =0.96,b ramp reactivity ρ =0.1βt and α =0.98,c ramp reactivity ρ =0.1βt and α =0.99

Table 3 Mean peak values of N（t）for ρex=0.5β,ρex=0.75β,and ρex= β,b comparison of mean peak values of N（t）for ρex=0.5β,ρex=0.75β,and ρex= β for α=1,and α=0.98

Fig.3(Color online)Mean N（t）and two arbitrary sample paths for a step reactivity ρex=0.5β and α =0.96,b step reactivity ρex=0.75β and α =0.96,c step reactivity ρex= β and α =0.96

Fig.4(Color online)Mean N（t）and two arbitrary sample paths for a step reactivity ρex=0.5β and α =0.98,b step reactivity ρex=0.75β and α=0.98,c step reactivity ρex= β and α =0.98

where S0n，S1n，S2n，...，Smnare random variables from N（0，1）with mean equal to 0 and variance equal to 1.

Theorem 4.1 The local truncation error of the scheme is O（τ3-α）.

Proof The local truncation error of the higher-order scheme for Eq.(4.2)can be derived as follows:

Fig.5(Color online)Mean N（t）and two arbitrary sample paths for a step reactivity ρex=0.5β and α =0.99,b step reactivity ρex=0.75β and α =0.99,c step reactivity ρex= β and α =0.99

5 Numerical results and discussions

Here,the solutions of SPKE(i=6)have been obtained using higher-order approximation scheme.

5.1 Step reactivity

In this section,the numerical solutions of the fractional stochastic point kinetic model[10]using the following parameters have been obtained: λ1=0.0127 s-1,λ2=0.0317 s-1, λ3=0.115 s-1, λ4=0.311 s-1, λ5=1.4 s-1, λ6=3.87 s-1, β1=0.000266, β2=0.001491,β3=0.001316, β4=0.002849, β5=0.000896, β6=0.000182,β=0.007,l=2.0×10-5s,v=2.5,and q=0 with N（0）=N0=100(neutron)and Ci（0）=

For different step reactivities,i.e.,ρ=0.003 and ρ =0.007,the mean peak values of N（t）with respect to time for fractional orders α=0.96,0.98,and 0.99 at step size h=0.001 s and for 500 trials are presented in Table 1.In addition,the results have been compared with those obtained by other methods,namely EM[8],Taylor 1.5 strong order[8],AEM[20],and efficient stochastic model(ESM)[10]with graphical representation in Fig.1.The solutions of SPKE obtained using the above-discussed fractional scheme have been tabulated to establish the efficiency of the higher-order approximation method.

5.2 Ramp reactivity

The numerical solutions of the fractional SNPKE have been obtained using the same parameters as those in Sect.5.1.Here,the reactivity can be represented as ρ =0.1βt.

For ramp external reactivity,i.e.,ρ =0.1βt,the mean peak values of N（t）with respect to time for fractional orders α=0.96,0.98,and 0.99 at step size h=0.001 s and for 500 trials are listed in Table 2.In addition,the results have been compared with those obtained using other methods,namely AEM[20]and ESM[10]with graphical representation in Fig.2.The solutions of SPKE obtained using the above-discussed fractional scheme have been tabulated to establish the efficiency of the higher-order approximation method.

5.3 Temperature feedback reactivity

In this section,solutions of fractional stochastic point kinetic model with i-group of delayed neutrons(i=6)in the presence of Newtonian temperature feedback have been obtained using the higher-order approximation method.

Table 4 Mean peak values of N（t）for=0.1t and 0.01t,b comparison of mean peak values of N（t）for=0.1t and 0.01t,for α =1 and α =0.98

Table 4 Mean peak values of N（t）for=0.1t and 0.01t,b comparison of mean peak values of N（t）for=0.1t and 0.01t,for α =1 and α =0.98

a σ α=0.96 α=0.98 α=0.99 Peak Time(s) Peak Time(s) Peak Time(s)(a)0.01 10-11 1.69869E+10 1.084 1.76031E+10 1.103 1.73211E+10 1.112 10-13 2.0663E+12 1.132 2.18262E+12 1.151 2.13422E+12 1.161 0.1 10-11 1.78236E+11 0.223 1.90873E+11 0.232 2.01965E+11 0.235 10-13 2.34211E+13 0.238 2.37627E+13 0.248 2.41795E+13 0.252 a σ SSFEMM[15](α =1) DFMM[15](α=1) AEM[10](α =1) α =0.98 Peak Time(s) Peak Time(s) Peak Time(s) Peak Time(s)(b)0.01 10-11 1.68604E+10 1.118 1.69492E+10 1.118 1.673436E+10 0.854 1.76031E+10 1.103 10-13 2.12034E+12 1.169 2.12802E+12 1.168 2.082531E+12 0.877 2.18262E+12 1.151 0.1 10-11 1.88849E+11 0.235 1.89642E+11 0.235 1.790577E+11 0.142 1.90873E+11 0.232 10-13 2.24448E+13 0.251 2.26026E+13 0.251 2.143778E+13 0.150 2.37627E+13 0.248

The total reactivity of the reactor in the presence of temperature feedback[21]is represented in the following form:

where σ = αKc,ρex（t）is the external reactivity,T（t）is the temperature,T0is the initial temperature,α is the temperature coefficient,and Kcis the reciprocal of the thermal capacity.

In the interval ［tk，tk+1］,the total reactivity can be expressed as follows[15]:

5.3.1 Step external reactivity

In this section,the numerical solutions of the SPKE of U235nuclear reactor[3]using the following parameters have been obtained: λ1=0.0124 s-1, λ2=0.0305 s-1,λ3=0.111 s-1, λ4=0.301 s-1, λ5=1.13 s-1, λ6=3.0 s-1, β1=0.00021, β2=0.00141, β3=0.00127,β4=0.00255,β5=0.00074,β6=0.00027,β =0.00645,l=5.0×10-5s, α=5.0×10-5K-1, and Kc=0.05 K/MWswith N（0）=N0=1(neutron)and Ci（0）=

For different step external reactivities,i.e.,ρex=0.5β,ρex=0.75β,and ρex= β,the mean peak values of N（t）for fractional orders α=0.96,0.98,and 0.99 and for different step external reactivities 0.5β,0.75β,and β using 500 trials are presented in Table 3a.These results have been compared with previously obtained results of split-step forward EM method(SSFEMM)and derivative-free Milstein method(DFMM)[15]in Table 3b with graphical representation in Figs.3,4,and 5.The solutions of SPKE obtained using the above-discussed fractional scheme have been tabulated to establish the efficiency of the higherorder approximation scheme.

For α =0.96,0.98,and 0.99,the mean N（t）with two arbitrary sample paths has been shown for ρex=0.5β,ρex=0.75β,and ρex= β.

5.3.2 Ramp external reactivity

The numerical solutions of the fractional SNPKE of U235nuclear reactor have been obtained using parameters similar to those in Sect.5.3.1.Here,the external reactivity is represented as ρex=0.1t and ρex=0.01t,the nonlinear coefficient σ takes 10-11or 10-13.

For different ramp external reactivities,i.e.,ρex=0.1t and 0.01t,the mean peak values of N（t）with respect to time for fractional orders α=0.96,0.98,and 0.99 at step size h=0.001 s and for 500 trials are listed in Table 4a.These results have been compared with the previously obtained results[10,15]in Table 4b with graphical representation in Figs.6a–c,7a–c,and 8a–c.The solutions of SPKE obtained using the above-discussed fractional scheme have been tabulated to establish the efficiency of the higher-order approximation method.

Fig.6(Color online)Mean N（t）and two arbitrary sample paths for a ramp reactivity ρex=0.01t, σ =10-11,and α =0.96,b ramp reactivity ρex=0.01t, σ =10-13,and α =0.96,c ramp reactivity ρex=0.1t, σ =10-11,and α =0.96,d ramp reactivity ρex=0.1t,σ=10-13,and α=0.96

Fig.7(Color online)Mean N（t）and two arbitrary sample paths for a ramp reactivity ρex=0.01t, σ =10-11,and α =0.98,b ramp reactivity ρex=0.01t, σ =10-13,and α =0.98,c ramp reactivity ρex=0.1t, σ =10-11,and α =0.98,d ramp reactivity ρex=0.1t,σ=10-13,and α=0.98

Fig.8(Color online)Mean N（t）and two arbitrary sample paths for a ramp reactivity ρex=0.01t,σ =10-11,and α=0.99,b ramp reactivity ρex=0.01t,σ =10-13,and α=0.99,c ramp reactivity ρex=0.1t,σ =10-11,and α=0.99,d ramp reactivity ρex=0.1t,σ =10-13,and α=0.99

Table 5 Mean peak values of N（t）for sinusoidal reactivity ρ =0.005333for different values of fractional order α

Table 5 Mean peak values of N（t）for sinusoidal reactivity ρ =0.005333for different values of fractional order α

ρ α=0.96 α=0.98 α=0.99 Peak Time(s) Peak Time(s) Peak Time(s)0.005333 sin 38.0005 38.28 45.8029 38.18 49.1345 38.99 πt T

5.4 Sinusoidal reactivity

In this section,the reactivity is in the form of sinusoidal change,i.e.,.The numerical solution for this reactivity has been obtained using the following parameters: ρ0=0.005333, β1= β=0.0079, λ1=0.077,Λ =0.001,q=0,N（0）=N0=1,and time period T=100 s.The mean peak values of N（t）with respect to time for fractional orders α=0.96,0.98,and 0.99 are presented in Table 5 with graphical representation in Fig.9.

Fig.9(Color online)Mean N（t）and two arbitrary sample paths for sinusoidal reactivity ρ=0.005333 sinwith a fractional order α =0.96,b fractional order α =0.98,and c fractional order α =0.99

6 Conclusion

In this study,fractional SPKEs have been solved using the higher-order approximation scheme with different fractional orders α.The obtained numerical solutions for mean N（t）have been presented in the tables and graphically demonstrated to justify the efficiency of the proposed higher-order approximation method.In addition,the results have been compared to some previous works such as[10,15,17].The results obtained by the implemented fractional model are in good agreement with previous results,which further establishes the efficiency of our proposed scheme.The graphical representation for different reactivities shows the behavior of the mean neutron population.The random fluctuations at low power levels and achievement of equilibrium state after reaching its peak value provide us with a succinct idea about the behavior of N（t）for different reactivities.