Assessment of the fractional neutron point kinetic equation to simulate core transients with Newtonian temperature feedback

Annals of Nuclear Energy 138 (2020) 107197

Contents lists available at ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

Assessment of the fractional neutron point kinetic equation to simulate core transients with Newtonian temperature feedback M.A. Polo-Labarrios a,⇑, S. Quezada-García b, G. Espinosa-Paredes c, J. Ortiz-Villafuerte d a

Departamento de Ciencias Naturales, Universidad Autónoma Metropolitana-Cuajimalpa, Avenida Vasco de Quiroga 4871, Santa Fe Cuajimalpa, 05370 Cuajimalpa, CDMX, Mexico Departamento de Sistemas Energéticos, Universidad Nacional Autónoma de México, Avenida Universidad 3000, Ciudad Universitaria, 04510 Coyoacán, CDMX, Mexico c Área de Ingeniería en Recursos Energéticos, Universidad Autónoma Metropolitana-Iztapalapa, Avenida San Rafael Atlixco 186, Vicentina, 09340 Iztapalapa, CDMX, Mexico d Departamento de Sistemas Nucleares, Instituto Nacional de Investigaciones Nucleares, Carretera México Toluca-La Marquesa s/n, Ocoyoacac, Estado de México C.P. 52750, Mexico b

a r t i c l e

i n f o

Article history: Received 4 June 2019 Received in revised form 25 October 2019 Accepted 4 November 2019

Keywords: Reactor dynamics Fractional neutron point kinetic equations Anomalous diffusion exponent Non-local time memory Newtonian temperature feedback reactivity Multi term higher-order linear approximation

a b s t r a c t The model of the Fractional Neutron Point Kinetic equations (FNPK) considers a relaxation time for neutrons, as they travel in surrounding media. Relaxation time is associated with a rapid variation in neutron flux, as reflection of fast change of reactivity, and it can be implemented in the differential operator of fractional order related with non-Fick effects, also known as anomalous diffusion exponent and represents non-local time memory. In this work the model of FNPK extended, by adding the temporal term to including the neutron density derivative, is using with the objective of studying the effect of this added term on the behavior of neutron density during a transient, when there exists Newtonian temperature feedback reactivity, using different values of the anomalous diffusion exponent. Although the one-group delayed neutron precursor is use for this analysis, the results of the neutron density behavior were compared with the classical neutron point kinetics equations. The results show that the additional term implemented in this study led to show how the sub-diffusion phenomena have noticeable relevance in the neutron density behavior in those scenarios where the dynamic of a nuclear system is associated to slow transients. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction The nuclear power distribution in a nuclear reactor implies calculation of neutron transport in heterogeneous media, where some areas are expected to show strong neutron absorption, besides scattering by the target nuclides or escaping from the active part of the reactor. To deal with neutron dynamics in a manner analogous to that of heat diffusion, it requires performing several approximations in the transport equation, which includes a weak angular dependence for the neutron angular distribution, isotropic neutron sources, and disregarding the derivative for the neutron current density respect to time, in comparison with other terms that appear in the neutron transport equation (Duderstadt and Hamilton, 1976). Transient situations occurring in a reactor core can be predicted through approximated solutions to the transport equation, to determine neutron flux and, as a result, it is possible to reach a sufficiently precise prediction of the consequences of the evolving transients. For that prediction, it is enough to relate the magnitude of the neutron flux, which varies in time, with the ⇑ Corresponding author. E-mail address: [email protected] (M.A. Polo-Labarrios). https://doi.org/10.1016/j.anucene.2019.107197 0306-4549/Ó 2019 Elsevier Ltd. All rights reserved.

neutron population in the core of a nuclear reactor (Henry, 1975). Point kinetics equations relate those parameters to describe the temporal behaviour of a nuclear reactor and thus allowing a study of transient situations that may occur in reactor core, by assuming that the duration of the transient the neutron flux spatial shape varies very little. Neutron distribution can be calculated directly from a sequence of approximations of the neutron transport equation, as the neutron diffusion equation, or according to Stacey (2001) and Henry (1975) through a heuristic procedure. All of them lead to the same functional form of the system of nonlinear coupled differential equations, but with different coefficients. The numerical solution of those equations is limited by the stiffness of the mathematical model, that results from the different time scales of the instantaneous and delayed neutrons (Chen et al., 2007, 2006a; Li et al., 2009). The classical neutron point kinetic equation model has been the subject of several studies and applications with different approaches. For example, Espinosa-Paredes et al. (2011) derived a Fractional Neutron Point Kinetics (FNPK) equation model, in which the derivative terms of the point kinetics equations have a noninteger order, that is adopting fractional derivatives (Aboanber and Nahla, 2016; Altahhan et al., 2016; Espinosa-Paredes, 2016).

2

M.A. Polo-Labarrios et al. / Annals of Nuclear Energy 138 (2020) 107197

The physical interpretation of this model is: the fractional order is related with non-Fick effects from the neutron diffusion equation point of view, and it is a consequence of anomalous diffusion phenomena processes that occur because of the highly heterogeneous configuration in nuclear reactors (Espinosa-Paredes and PoloLabarrios, 2012). The model considers a relaxation time for neutrons, as they start their motion associated with a rapid variation in the neutron flux, due to the fast variation of reactivity. The FNPK equation constitute a useful tool to provide important information on the dynamics of the reactor and it have been analyzed and applied in different works (Das et al., 2013; Espinosa-Paredes, 2017; Espinosa-Paredes et al., 2017; Hamada and Brikaa, 2017; Nahla, 2017; Nahla and Hemeda, 2017; Nowak et al., 2015, 2014b, 2014a; Ray, 2015; Ray and Patra, 2015, 2013, 2012; Schramm et al., 2013; Vyawahare and Nataraj, 2013b, 2013a; Vyawaharea and Espinosa-Paredes, 2018). A relevant application of the FNPK equation is that by Cázares-Ramírez et al. (2017), who investigate the stability of linear fractional-order neutron point kinetics models in the closed-loop configuration. These authors selected three values of the anomalous diffusion exponent to represent three levels of sub-diffusion in the reactor core. The authors succeed to obtain a best representation against of the plant data of a transient compared to the CNPK equations. There are several methods available in the literature to consider the Newtonian temperature feedback reactivity using classical neutron point kinetic equation (cf. e.g. Chen et al., 2007, 2006b; Hamieh and Saidinezhad, 2012; Li et al., 2007; Nahla, 2011, 2009). In this work, the fractional neutron point kinetics equation formulated by Espinosa-Paredes et al. (2011), extended by adding the temporal term proposed by Hamada (2017) to including the neutron density derivative, is used with the objective of studying the effect of this additional term on the behaviour of neutron density during a transient, when there exists Newtonian temperature feedback reactivity. The one-group delayed neutron approximation is used for this analysis. The results of neutron density obtained from the modified FNPK model are calculate using the numerical solution proposed by Polo-Labarrios et al., (2020). This paper presents a comparison between the solution method used in this work against solution method developed by Hamada (2017), applied to the case of sinusoidal reactivity. Applying the both numerical methods developed with a time step size h ¼ 0:005s, for a ¼ 0:2; 0:3; 0:4, the results obtained shows that a significant difference exists between the neutron density results obtained by the solution method developed by Hamada (2017) and the new solution method applied in this work. This difference is attributed to the fact that Hamada (2017) considered the added term as a fractional differential operator, while in this work is treated as a fractional differential-integral operator. Details of Fractional Calculus fundamentals can be found in several references, as for example (Diethelm, 1997; Gorenflo and Mainardi, 1997; Hilfer, 2000; Podlubny, 1998). In this work, the model of FNPK extended, by adding the temporal term to including the neutron density derivative, is using with the objective of studying the effect of this added term on the behavior of neutron density during a transient, when there exists Newtonian temperature feedback reactivity, using different values of the anomalous diffusion exponent. Although the one-group delayed neutron precursor is use for this analysis, the results of the neutron density behavior were compared with the classical neutron point kinetics equations. The solution of the FNPK extended is by the application of the fractional derivative and the fractional integral (Polo-Labarrios et al., 2020). The results show that the additional term implemented in this study led to show how the sub-diffusion phenomena have noticeable relevance in the neutron density behavior in those scenarios where the dynamic of a nuclear system is associated to slow transients. Physically, the magnitude of the anomalous diffusion

exponent reflects clearly on the dynamics of a system undergoing some perturbation. 2. Modified fractional neutron point kinetic mathematical formulation To considering the reactivity as time dependent function Hamada (2017) extends the Fractional Neutron Point Kinetics (FNPK) model, based in the work of Espinosa-Paredes et al. (2011), to obtain the following Equation:

sa

aþ1

d

nðt Þ

þ

 a
dnðt Þ d PNL ð1  qðt ÞÞ  1 þ b þ sa a nðtÞ dt dt K

dtaþ1 a qðtÞ  b d C ðt Þ nðt Þ þ sa k þ kC ðtÞ ¼ dta K

ð1Þ

where nðt Þ is the neutron density, PNL is the neutron non-leakage probability, K is one generation average lifetime of instantaneous neutron, b is the total fraction of delayed neutron, k is the decay constant of delayed neutron precursor, C ðt Þ is the delayed neutron precursor, and s is the relaxation time defined as

1

s¼

tRtr

¼

3D

ð2Þ

t

where t is the neutron velocity, Rtr ðr; t Þ is the transport crosssection, and D is the diffusion coefficient. Hamada (2017) discussed the applicability of Eq. (1) for finite and infinite media. In this study, the Eq. (1) is intended for a finite medium, where neutrons are generated by an external source or by fission, and lost because of scattering collisions, absorption or leaking from the system. The third term on the left side of Eq. (1) can be rewritten as:

 a
a d PNL ð1  qðt ÞÞ  1 þ b d n ð t Þ ¼ a ½Pðt ÞnðtÞ a dt K dt

ð3Þ

where

Pðt Þ ¼

PNL ½1  qðt Þ  1 þ b K

ð4Þ

Then, the Eq. (1) takes the following form:

sa

aþ1

d

dt

nðt Þ

aþ1

þ

a

dnðt Þ d qðtÞ  b nðt Þ þ sa a ½PðtÞnðtÞ ¼ dt dt K a d C ðt Þ þ kC ðtÞ þ sa k dta

ð5Þ

Therefore, in this way, the time derivative of the reactivity can be taking into a count during the process to obtain a more general FNPK model. However, this consideration is only valid for dynamic systems with feedback effects. For cases where the reactivity inserted is time dependent from the very beginning, it is needed the specific time dependent function, for example, the simple case when reactivity is inserted as a step change and it takes a constant value during the whole transient. Now, considering the constraint of 0 < a < 1, Hamada (2017) applying the generalized Leibniz rule for the fractional order derivative to the product of two functions (Podlubny, 1998): a

1 X d a ½P ðt Þnðt Þ ¼ dt i¼0

 

a i

PðiÞ ðtÞ

d

ai

nðt Þ dt ai

ð6Þ

  a ¼ aða  1Þða  2Þ    binomial coefficients i   a ða  i þ 1Þ=i! (with ¼ 1, aeRþ, and ieN) do not become zero 0 wherein

the

whenever i > a because a R N. P ðiÞ ðt Þ is the i  th derivative of Pðt Þ.

3

M.A. Polo-Labarrios et al. / Annals of Nuclear Energy 138 (2020) 107197

Note that the binomial sum is now of infinite extent, then this equation takes the form a

a

a1

d d nðt Þ dPðtÞ d nðt Þ ½Pðt Þnðt Þ ffi Pðt Þ þa dta dt dt a dta1   2   3 a2 a d PðtÞ d nðtÞ a d P ðtÞ da3 nðtÞ þ þ þ ... 2 a 2 2 3 dt dt dt3 dt a3

d

aþ1

nðtÞ

ð7Þ

dt aþ1   ¼


 a dnðtÞ PNL ð1  qðtÞÞ  1 þ b d nðt Þ þ sa dt dta K

þ

asa PNL dqðtÞ da1 nðtÞ K

dta1

dt



dt 2

dt a2

aða  1Þða  2Þsa PNL d3 qðtÞ da3 nðtÞ 6K

qðtÞ  b K

dt a3

dt 3

nðtÞ þ sa k

a

d C ðt Þ þ kC ðt Þ dt a

ð8Þ

for 0 < a < 1. The fractional model given by above Equation includes the time variation of the reactivity terms d qðtÞ=dti for i ¼ 1; :::; 1, which are relevant for those transients when the absorption cross-section parameter Ra ðt Þ is not constant (Ray and Patra, 2015), as in transient with feedback effects or for dynamic systems. It is the main purpose of this study to perform an analysis of the impact of such terms on the dynamics of the reactor core during certain transient scenarios. It is quite important to highlight in previous equation that, as a result of the application of fractional operators, the first and fourth to six terms on the left side, because these terms require special attention during the process of achieving a numerical solution: i

First term :

aþ1

d

nðtÞ

dt aþ1

Fourth to six terms :

¼

 a  1þa d dnðt Þ d nðt Þ ; – a dt dt dt1þa

ai

d

nðt Þ

dtai

¼

a

d dta

ð9Þ

! i a iþa d nðt Þ d d nðtÞ ; ¼ a J i nðt Þ– i dt dt dt iþa ð10Þ

Eq. (8) describe the neutron density dynamic behavior in a nuclear reactor core when there is a transient with variation of reactivity. In particular, the fourth term implies using both a fractional derivative and a fractional integral. To solve this equation, it is necessary to consider the order of the operators in Eqs. (9) and (10) in the numerical method (Polo-Labarrios et al., 2020), to solve the fractional ordinary differential equation system. The particular cases: when the reactivity insertion is a constant, namely, that at time t ¼ 0 a disturbance is inserted and subsequently the reactivity remains constant, Eq. (8) is reduced to:


 a d nðtÞ dnðtÞ PNL ð1  qÞ  1 þ b d nðt Þ þ sa sa aþ1 þ dt dta K dt aþ1

¼

qðtÞ  b K

nðtÞ þ sa k

a

d C ðt Þ þ kC ðt Þ dt a

dt

nðtÞ

aþ1

þ


 a dnðt Þ PNL ð1  qðt ÞÞ  1 þ b d nðtÞ þ sa dt dt a K a1

a

c asa PNL d nðtÞ qðt Þ  b d C ðt Þ ¼ nðtÞ þ sa k þ kC ðt Þ dt a K K dt a1

ð11Þ

In this case q ¼ constant and their subsequences derivatives are cero and no additional terms appear. This mean that the form and quantity of additional terms that will appear in modified FNPK equation will depend on the form of the reactivity parameter function according to each transient event in a nuclear reactor, and the number of additional terms will correspond to the order of the approximation to FNPK equation. For example, the case when the reactivity inserted is time dependent as q(t) = ct, where c is an arbitrary constant, the Eq. (8) take the following form

ð12Þ

where only one additional term appears. So, to the previous equation correspond with a first order approximation to FNPK equation, which corresponding to the existence of the first reactivity function derivative. For the case when reactivity is a sinusoidal function qðtÞ ¼ q0 sinðpt=T Þ, whose i  th derivative exist, if limiting the value of i ¼ 3 it have:

p pt  d2 q p2 pt  d3 q dq ¼ q0 sin ¼ q0 cos 2 T dt T dt3 dt T T   p3 pt ¼ q0 cos : T T

aða  1Þsa PNL d2 qðtÞ da2 nðtÞ 2K

aþ1

d



Now, substituting this approximation into Eq. (5), leads to

sa

sa

ð13Þ

in this case, the modified FNPK equation obtained corresponds to a third order approximation to FNPK equation, and the Eq. (8) takes the following form:

sa


 a dnðt Þ a P NL ð1  qðt ÞÞ  1 þ b d nðt Þ s þ dt dt a K dt aþ1 p pt  da1 nðt Þ aða  1Þ p2 pt  da2 nðtÞ  ax þ x sin cos T 2 T T T dt a1 dt a2   a3   3 aða  1Þða  2Þ p pt d nðtÞ qðtÞ  b þ x cos ¼ nðtÞ 6 T T K dt a3 aþ1

d

nðtÞ

þ sa k

þ

a

d C ðt Þ þ kC ðt Þ dt a

ð14Þ

a

where x ¼ s PKNL q0 . Finally, the differential equation of delayed neutron precursors is:

dC ðt Þ b ¼ nðt Þ  kC ðtÞ: dt K

ð15Þ

The first and second terms on the right side of this equation are the rate of production of the precursors and radioactive decay, respectively. Eqs. (8) and (15) are a system of stiff nonlinear fractional ordinary differential equations, which can be used to analyse the effects of sub-diffusion processes, due the anomalous diffusion exponent and relaxation time in the neutron kinetics for a finite medium. Such equation system incorporates a time varying reactivity, because it depends mainly on the absorption cross section, which in turn is time dependent. Note that when s ¼ 0, the classical neutron point kinetic equation model is obtained directly from Eq. (8), together with Eq. (15). 3. Newtonian temperature feedback reactivity application The reactivity is one of the most important parameters in the operation of a nuclear reactor, because it is directly related to the control of the reactor (Espinosa-Paredes et al., 2014). Reactivity, positive or negative, can be inserted by the reactor operator via movement of the control rods. However, uncontrolled control rod withdrawal or ejection is a common type of initiator event for a reactivity insertion accident. The consequent change of coolant temperature needs to be quantified. One approach is the application of the Newtonian temperature feedback reactivity model, which has been studied in several works (e.g. Duderstadt and Hamilton, 1976; Glasstone and Sesonske, 2012; Palma et al., 2009; Stacey, 2001; Zhang et al., 2008). In this work, Eq. (8) is applied to that same case of Newtonian temperature feedback reactivity, to study the effects of the time variation of the reactivity, by quantifying the effect of the additional terms on the behavior of neutron density.

4

M.A. Polo-Labarrios et al. / Annals of Nuclear Energy 138 (2020) 107197

In the presence of Newtonian temperature feedback reactivity, the reactivity is written as (Glasstone and Sesonske, 2012; Hetrick, 1993; Stacey, 2001):

method developed at section 4 of Polo Labarrios et al. (Acepted). This numerical method is based on the numerical algorithm given by Edwards et al. (2002). To simplify notation, the differential

qðtÞ ¼ q0  aq ½T ðtÞ  T 0 ;

operator Da is used instead of d =dta , and J b instead d =dtb . The equations system as a linear multi-term fractional form can be re-written as: Fractional point kinetics

ð16Þ

where q0 is the initial reactivity, aq is the reactivity coefficient, T 0 is initial temperature, and T is the prompt temperature given by:

dT ðt Þ ¼ K c nðtÞ; dt

ð17Þ

where K c is the reciprocal of thermal capacity of the reactor. Differentiating Eq. (16) with respect to time t to take into account the additional terms at Eq. (8), we have: dq dt

¼

d2 q dt 2

dq dT ðtÞ dT dt

¼ ddTq2

ð18Þ

d2 T ðt Þ dt2

2

where, dq dT

Daþ1 nðt Þ þ a3 Dnðt Þ þ a2 Da nðtÞ þ a1 nðt Þ þ a5 Dqðt ÞDa Jnðt Þ ¼ b2 Da C ðt Þ þ b1 C ðtÞ;

ð19Þ

¼ 0:

DC þ b2 C ¼ a0 nðt Þ:

DT ðt Þ ¼ K c nðtÞ:

Dqðt Þ ¼ aq K c nðt Þ;

b ; K

a1 ¼ 

ð20Þ

Introducing Eq. (20) in Eq. (8), it is obtained a first order approximation to FNPK equations with Newtonian temperature feedback reactivity effects in a nuclear reactor:

nðtÞ

dt aþ1 þ

þ


 a dnðtÞ PNL ð1  qðt ÞÞ  1 þ b d nðtÞ þ sa dt dt a K

asa PNL aq K c da1 nðtÞ K

dt

a1

nðt Þ ¼

qðtÞ  b K

nðt Þ þ sa k

¼ 0; t¼0

 d nðt Þ dt a  a

t¼0

 d nðtÞ ¼ 0;  dt a1  a1

¼ 0;

a

d C ðt Þ þ kC ðt Þ: dt a ð21Þ

nðt ¼ 0Þ ¼ n0 ;

t¼0

ð22aÞ Precursor concentration initial conditions:

 d C ðtÞ ¼ 0; a  dt t¼0 a

C ð0Þ ¼

a2 ¼ a3 ¼

The necessary initial conditions by the following equations are given: Fractional neutron point kinetics initial conditions:

 dnðt Þ dt 

ð24Þ

ð25Þ

ð26Þ

where the coefficients in the equations above are:

dqðt Þ ¼ aq K c nðt Þ: dt

aþ1

ð23Þ

Prom temperature

a0 ¼

d

0 < a  1;

Precursor concentration

Then, only one additional term appears at Eq. (8), it is obtained using the Eqs. (17) and (19) into Eq. (18):

sa

b

Newtonian temperature feedback reactivity

¼ aq;

d2 q dT 2

a

b n0 : Kk

ð22bÞ

Newtonian temperature initial conditions:

T ðt ¼ 0Þ ¼ 0:

ð22cÞ

Newtonian temperature feedback reactivity initial conditions:

qð0Þ ¼ 0;

ð22dÞ

and, the initial reactivity q0 is inserted into the reactor at t ¼ 0. The initial conditions given by Eqs. (22) satisfy the static initial condition. The Eqs. (15), (16), (20) and (21) form a system of stiff nonlinear fractional ordinary differential equations to describe the Newtonian temperature feedback reactivity effect in nuclear reactor. 4. Numerical solution of nonlinear multi-term fractional differential equations In order to solve the system of stiff nonlinear multi-term fractional ordinary differential equations, we apply the numerical

a5 ¼ b1 ¼

ð27Þ

qðtÞ  b ; sa K

PNL ð1  qðtÞÞ  1 þ b ; K 1

sa

;

aP NL aq K

k

sa

ð28Þ ð29Þ ð30Þ

;

;

b2 ¼ k:

ð31Þ ð32Þ ð33Þ

The equation system of this problem is transformed in an ordinary and fractional differential equations (OFDE) system. Represented as a multi-term high-order linear fractional differential equation, which is calculated by writing the problem as a system of OFDE, as follows variable change: Fractional point kinetics:

x1 ðtÞ ¼ nðtÞ;

ð34Þ

x2 ðtÞ ¼ Da nðtÞ ¼ Da x1 ðtÞ;

ð35Þ

x3 ðtÞ ¼ DnðtÞ ¼ Dx1 ðtÞ;

ð36Þ

x4 ðtÞ ¼ D1 nðtÞ ¼ Jnðt Þ ¼ Jx1 ðt Þ;

ð37Þ

x5 ðtÞ ¼ Da D1 nðt Þ ¼ Da x4 ðt Þ;

ð38Þ

x6 ðtÞ ¼ Da DnðtÞ ¼ Daþ1 nðtÞ ¼ Da x3 ðtÞ:

ð39Þ

Precursor concentration

y1 ðtÞ ¼ CðtÞ;

ð40Þ

y2 ðtÞ ¼ Da CðtÞ ¼ Da y1 ðtÞ;

ð41Þ

y3 ðtÞ ¼ DCðtÞ ¼ Dy1 ðtÞ:

ð42Þ

5

M.A. Polo-Labarrios et al. / Annals of Nuclear Energy 138 (2020) 107197

And the Prompt temperature equation, Eq. (25), and the Newtonian temperature feedback reactivity equation (Eq. (26)), are now represented as: Prompt temperature

t 1 ðtÞ ¼ T ðt Þ;

ð43Þ

S5;i ¼

i X

a

xp;i x3;ip þ

x3;0

a

xp;i y1;ip þ

y1;0

a

p¼1

S6;i ¼

i X p¼1

t 2 ðtÞ ¼ DT ðtÞ ¼ Dt1 ðt Þ:

ð44Þ S7;i ¼ y1;i1 þ

Newtonian temperature feedback reactivity

z1 ðt Þ ¼ qðt Þ;

ð45Þ

z2 ðt Þ ¼ Dqðt Þ ¼ Dz1 ðtÞ;

ð46Þ

where the term x6 is obtained by substituting Eqs. (34)–(38) and Eq. (41) into Eq. (23); and y3 was obtained by substituting Eqs. (36) and (41) into Eq. (24), that is

x6 ¼ 

5 X

aj xj þ

2 X

j¼1

bj yj ;

ð47Þ

j¼1

y3 ¼ a0 x1  b2 y1 :

ð48Þ

Note that a4 ¼ 0. In similar way, t 2 and z2 are obtained substituting the Eq. (36) into Eqs. (25) and (26), respectively:

t 2 ¼ K c x1 ;

ð49Þ

z2 ¼ aq K c x1 :

ð50Þ

The influence of the order of the operators mentioned in Eq. (8) is taken into account during the process of order reduction of the fractional differential equation through Eqs. (37)–(39), corresponding to the variables x4 , x5 and x6 , respectively. Now, following the procedure developed by Polo Labarrios et al. (Acepted), the numerical solution in matrix form is:

0

a x0;i

ac

i

B 1 0 B B B h=2 0 B B 0 0 B B B a c a1 a ci a2 i B B 0 0 B B B ða0 hÞ=2 0 B B 0 @ hr 0 =2 hT 0 =2

0

i X

a

xp;i x1;ip þ

x1;0

a

;

h S2;i ¼ x1;i1 þ x3;i1 ; 2

S4;i ¼

0

0

0

0

0

0

1

0

0

0

0

0

a x0;i

ac

0

0

0

0 0

a ci a5 0

0

0

0

1 þ ðb2 hÞ=2

0

0

0

0

0

0

0

1

0

0

0

0

0

0

p¼1

a

xp;i x4;ip þ

ð52Þ

ð53Þ

h x1;i ; 2 i X

The first problem to be attacked is described by the one-speed diffusion model of a finite, bare homogenous cylindrical reactor with material composition similar to that of a pressurized water reactor (PWR). This exercise allows studying the stability of the proposed numerical method when no feedback is expected, as in

ac

i a5

ac

i

0 0

0

10

x1

1

0

S1

1

B C B C 0C C B x2 C B S2 C CB C B C B C B C 0C C B x3 C B S3 C CB C B C 0 C B x4 C B S4 C CB C B C B C B C 0C C B x5 C ¼ B S5 C CB C B C 0 C B y 1 C B S6 C CB C B C B C B C 0C C B y 2 C B S7 C CB C B C 0 A @ z 1 A @ S8 A 1 t1 S9

ð51Þ

the case of a sudden reactivity insertion. The data used in this first computation are given next. The fuel is taken as UO2 enriched to

p¼1

S3;i ¼

5.1. Nuclear parameters

h=2

where

S1;i ¼

The numerical solution process proposed for the system of Eqs. (51)–(60) is applied to the case of modified FNPK equations with one-group delayed neutron precursor, when Newtonian temperature feedback reactivity effects are taken into account. The numerical procedure was implemented in the software MATLABÒ. In this work, the Gaussian elimination method to invert the matrix given by Eq. (51) was applied.

0

0;i

ð58Þ

5. Simulations

0

a x

h a0 x1;i1 þ b2 y1;i1 ; 2

ð60Þ

0

i b1

ð57Þ

h S9;i ¼ t 1;i1 þ T 0 x1;i1 : 2

0

ac

;

ð59Þ

0

i

ð56Þ

h S8;i ¼ z1;i1 þ r 0 x1;i1 ; 2

0

ða ci a3 þ a x0;i Þ 0

a

;

ð54Þ x4;0

a

;

ð55Þ

2:78% of U 235 . The other parameter values needed for the solution are: Dc ¼ 9:21 cm, PNL ¼ 0:975, and the neutron speed t ¼ 220000 cms , which corresponds to the energy of the slow (thermal) neutrons (0:025eV) (Hetrick, 1993). The value for the temperature reactivity coefficient is aq ¼ 5  105 K 1 , the average lifetime of instantaneous neutron is K ¼ 0:00005s, the reciprocal of thermal capacity of the reactor is K c ¼ 0:05 K=MW, and for the decay constant of delayed neutron precursor and the total fraction of delayed neutron used in this paper are, respectively:

b k ¼ P6

bi i¼1 ki

¼ 0:0769478

ð61Þ

6

M.A. Polo-Labarrios et al. / Annals of Nuclear Energy 138 (2020) 107197

Table 1 Neutron delayed fractions and decay constants, six precursor group parameters considered (Chao and Attard, 1985; Nobrega, 1971). Decay constants and Yield values  i ki s1

bi

1 2 3 4 5 6

0.00021 0.00141 0.00127 0.00255 0.00074 0.00027

b¼

0.0124 0.0305 0.1110 0.3010 1.1300 3.00

6 X

bi ¼ 0:00645;

ð62Þ

i¼1

where bi and ki are present in Table 1 (Chao and Attard, 1985; Nobrega, 1971). All simulations started from the equilibrium conditions with a normalized neutron density (n0 = 1) and one-group of delayed neutron density (C0 = bn0/Kk). The stability limit for time-step size (h) is determined for a value of the anomalous diffusion exponent

a ¼ 0:5, and the relaxation time s ¼ 1:2559  104 s. Fig. 1 shows the effect of different step sizes h on the neutron density (Fig. 1a) and total reactivity (Fig. 1b) as the transient evolves, obtained from the modified FNPK equation with a step reactivity change of q0 ¼ 0:9b. Fig. 1a shows that only the two largest values of h (h ¼0.1 s, 0.05 s), lead to some noticeable higher values of the peak value of the neutron density, while in Fig. 1b, there is no appreciable change in the total resulting reactivity during the whole transient. It is clear then that the numerical scheme introduced has a very stable performance for time steps in a range of at least two orders of magnitude (0.001 s to 0.1 s). This is a relevant result because of the potential computational time savings, when considering the need of performing many calculations. Then, the time step value h and the total simulation time considered for these numerical experiments is depended of the initial reactivity value, taken as constant during the whole transient scenario.

6. Temperature feedback reactivity

Fig. 1. Effect of the step sizes h ¼ 0:1; 0:05; 0:02; 0:01; 0:005; 0:001s at time on the neutron densities and reactivity obtained by the modified FNPK equation for q0 ¼ 0:9b and a ¼ 0:5.

In order to investigate the joint effect of the anomalous diffusion exponent (a) and the relaxation time (s) on the behavior of the resulting neutron density, temperature, and reactivity, a series of numerical simulations were carried out, with different values of the anomalous diffusion exponent for different values of the initial reactivity q0 . The transient behavior of neutron density, reactivity, and ‘‘Newtonian” temperature is shown in Figs. 2–4 for an initial reactivity values q0 ¼ 0:1b, q0 ¼ 0:5b, and q0 ¼ 0:9b, respectively. Each of these three cases consider six different values for the anomalous diffusion exponent a ¼ 0:9; 0:7; 0:5; 0:3; 0:1; and a step size h ¼ 0:005s; 0:005s; 0:01s; respectively. These time steps were chosen because they were already shown to produce good approximation and the computational cost is low. In Fig. 2a, it is noted that for those values of the anomalous diffusion exponent 0:3  a  0:7, neutron density was under predicted when compared to the results obtained with the classical model. However, when a  0:9, tending to a ! 1, the neutron density profiles now show a similar behavior to that from the classical model. For a ¼ 0:1, the behavior of neutron density predicted by the modified model FNPK equation over predicts the corresponding behavior described by the classical model. In this last comparison, a maximum value of 4:68 is obtained with the modified model FNPK, while with the classic model the maximum value

equals 2:0, that is the difference between those maximum values reaches about 134%. Therefore, for values of the anomalous diffusion exponent greater than 0.7, there is a close agreement between the neutron density results obtained by the modified model FNPK equation and those result obtained the classical point kinetics model. However, values of the anomalous diffusion exponent less than 0.5 lead to clear deviations from expected results, and even to unrealistic results, as shown for the case a ¼ 0:1. Regarding the results obtained for the reactivity profile during the transient, Fig. 2b shows that the modified FNPK equation clearly over predicts the reactivity profiles obtained with the classical point kinetics model when the anomalous diffusion exponent is between 0.3 and 0.5. As expected, the case a ¼ 0:1 leads to unrealistic values, and values greater than 0.7 lead to a close agreement with results from the classical point kinetics model. Regarding the change in temperature, Fig. 2c shows that this variable resulting data can be described in the same terms as discussed for the neutron density, when compared to results from the classical point kinetics model, for the different values of the anomalous diffusion exponent used in the simulation. For the cases q0 ¼ 0:5b and q0 ¼ 0:9b, Figs. 3 and 4 show that for all values of the anomalous diffusion exponent used in the modified FNPK equation, the profiles of neutron density, total

M.A. Polo-Labarrios et al. / Annals of Nuclear Energy 138 (2020) 107197

Fig. 2. Results from the Modified FNPK equations including Newtonian temperature feedback reactivity effects, for an initial reactivity value q0 ¼ 0:1b and step size h ¼ 0:005s. (a) Time dependent neutron density; (b) Time dependent reactivity; and (c) Time dependent temperature difference.

7

Fig. 3. Results from the Modified FNPK equations including Newtonian temperature feedback reactivity effects, for an initial reactivity value q0 ¼ 0:5b and step size h ¼ 0:005s. (a) Time dependent neutron density; (b) Time dependent reactivity; and (c) Time dependent temperature difference.

8

M.A. Polo-Labarrios et al. / Annals of Nuclear Energy 138 (2020) 107197

numerical values. In Fig. 3a, it is noted that for a ¼0.5, 0.3, the behavior of the neutron density obtained with the modified FNPK equation under predicts the result obtained with the classical model, until about 85s of simulation time. After that time, the classical model predicts slightly higher values. For a ¼ 0:1, the difference in neutron density values are noticeable. For this case, the modified FNPK equation over predicts the values obtained by the classical model even at the steady state. The maximum values of neutron density occurred at 59:27s and 67:22s, for the modified FNPK equation and the classical point kinetics model, respectively. After reaching the maximum value, the results from the modified FNPK equation fall under those obtained with the classical model. Regarding the results of the total reactivity, Fig. 3b shows that with exception of a ¼ 0:1, the rest of values of the anomalous diffusion exponent led to reactivity profiles slightly higher that the results by the classical model. The case a ¼ 0:1, on the contrary, led to smaller values of reactivity, but the profile is still consistent with that obtained by the classical model. Regarding the temperature results, Fig. 3c shows that with exception of a ¼ 0:1, the rest of values of the anomalous diffusion exponent led to reactivity profiles slightly lower that the results by the classical model. The case a ¼ 0:1, on the contrary, led to higher values of reactivity, but, as occurred with the reactivity results, the profile is still consistent with that obtained by the classical model. For the case q0 ¼ 0:9b, Fig. 4a, shows that using a ¼ 0:3; 0:1, the behavior of the neutron density obtained with the modified FNPK equation under predicts the profile obtained with the classical model up to 27s. After this time, the behavior obtained by the modified model of the FNPK equation surpassed the profile obtained with the classical model. The classic model yielded a maximum value at 16:15 s, while with the anomalous diffusion exponent value of a ¼ 0:3 the maximum value of neutron density obtained at 16:9s, and the maximum value obtained with for a ¼ 0:1 at 20:45s. The results shown in Figs. 2–4 clearly indicate that slow transients with some feedback effects included, induced by mall amounts of reactivity inserted to start the perturbation, can only be simulated with enough precision if the anomalous diffusion exponent values are above 0.7. For those cases when the reactivity inserted are in a range of medium to high values (q0 ¼ 0:5b and q0 ¼ 0:9b), that is fast transients, the modified FPK equation proposed in this work yields results with comparable precision to those obtained with the classical point kinetics model, with the exception when using small values of the anomalous diffusion exponent. Physically, the magnitude of the anomalous diffusion exponent reflects clearly on the dynamics of a system undergoing some perturbation. For slow transients (q0 ¼ 0:1b), the phenomena of subdiffusion proves being quite relevant in the dynamics of neutron density. However, in fast transients, the change in the absorption cross-section is so rapid that its variation is not noticeably reflected on neutron density, and consequently on reactivity and Newtonian temperature. In any case, this work shows that choosing an appropriate value for the anomalous diffusion exponent is paramount to obtain reliable results.

7. Conclusion Fig. 4. Results from the modified FNPK equations including Newtonian temperature feedback reactivity effects, for an initial reactivity value q0 ¼ 0:9b and step size h ¼ 0:005s. (a) Time dependent neutron density; (b) Time dependent reactivity; and (c) Time dependent temperature difference.

reactivity and change in temperature are in close agreement with those profiles obtained with when using the classical model. Only the case when a ¼ 0:1 leads to noticeable difference in the

Recently, Hamada (2017) had proposed a modified to the model of the fractional neutron point kinetic equation, introducing a new term related with the derivative of reactivity. In this sense, this paper a novel numerical solution for the modified fractional neutron point kinetics equations was used. In order to analyze the effect of the additional term over the neutron density behavior descript by modified fractional neutron point kinetic equation we

M.A. Polo-Labarrios et al. / Annals of Nuclear Energy 138 (2020) 107197

the numerical solution obtained is applied for case of Newtonian temperature feedback reactivity effects in a nuclear reactor, for one-group delayed neutron precursor, and for different values of the anomalous diffusion exponent. The results obtained with a novel solution for modified FNPK model was compared with the classical neutron point kinetic equations for simulated cases. We can see that the additional term, is relevant to analyze the information on the dynamic of the nuclear reactor, and such determine the effect on the neutron density behavior associated with slow transients (q0 ¼ 0:1b). Where the dynamic system is affected by the change of the absorption cross-section with time, which in turn modifies the reactivity value (q0 ) due to feedback from the nuclear reactor. With the value of the reactivity close to zero, the results obtained with the modified FNPK equation show greater differences with respect to the results obtained with the classical point kinetic neutron model. With the values of anomalous diffusion exponent that present greater difference are a ¼ 0:1; 0:3; 0:5, this means that in slow transients the phenomena of sub diffusion have greater relevance in the behavior of the neutron density, in comparison with the fast transients, where the differences are smaller. This indicates that the time derivative for the neutron current density should not be disregarded in slow transient analysis, like a several accident situations. In any case, this work shows that choosing an appropriate value for the anomalous diffusion exponent is paramount to obtain reliable results. Acknowledgements The authors greatly acknowledge the support of PRODEP to the project 47410644 through their NPTC Fund, and the support of SENER-CONACYT to the project FSE-2013-04-213519 through their Energy Sustainability Fund. References Aboanber, A., Nahla, A.A., 2016. Comment on the paper: Espinosa-Parrdes, et al., 2011. Fractional neutron point kinetics equations for nuclear reactor dynamics. Ann. Nucl. Energy 38, 307–330. Ann. Nucl. Energy 88, 301–302. Altahhan, M.R., Nagy, M.S., Abou-Gabal, H.H., Aboanber, A.A., 2016. Formulation of a point kinetics model based on the neutron telegraph equation. Ann. Nucl. Energy 91, 176–188. Cázares-Ramírez, R.I., Vyawahare, V.A., Espinosa-Paredes, G., Nataraj, P.S.V., 2017. On the feedback stability of linear FNPK equations. Prog. Nucl. Energy 98, 45– 58. https://doi.org/10.1016/j.pnucene.2017.02.007. Chao, Y.A., Attard, A., 1985. A resolution of the stiffness problem of reactor kinetics. Nucl. Sci. Eng. 90, 40–46. Chen, W.Z., Guo, L.F., Zhu, B., Li, H., 2007. Accuracy of analytical methods for obtaining supercritical transients with temperature feedback. Prog. Nucl. Energy 49, 290–302. Chen, W.Z., Kuang, B., Guo, L.F., Chen, Z.Y., Zhu, B., 2006. New analysis of prompt supercritical process with temperature feedback. Nucl. Eng. Des. 236, 1326– 1329. Das, S., Mukherjee, S., Das, S., Pan, I., Gupta, A., 2013. Continuous order identification of PHWR models under step-back for the design of hyper-damped power tracking controller with enhanced reactor safety. Nucl. Eng. Des. 257, 109–127. Diethelm, K., 1997. An algorithms for numerical solution of differential equations of fractional order. Electron. Trans. Numer. Anal. 5, 97–126. Duderstadt, J.J., Hamilton, L.J., 1976. Nuclear Reactor Analysis. Jhon Wiley and Sons, United States of America. Edwards, J.T., Ford, N.J., Simpson, A.C., 2002. The numerical solution of linear multiterm fractional differential equations. J. Comput. Appl. Math. 148, 401– 418. Espinosa-Paredes, G., 2017. Fractional-space neutron point kinetics (F-SNPK) equations for nuclear reactor dynamics. Ann. Nucl. Energy 107, 136–143. Espinosa-Paredes, G., 2016. A note on ‘‘Comment on the paper: Espinosa-Parrdes, et al., 2011. Fractional neutron point kinetics equations for nuclear reactor dynamics. Ann. Nucl. Energy 38, 307–330”. by Aboanber, A.E., Nahla, A.A. Ann. Nucl. Energy 88, 303–304. Espinosa-Paredes, G., Cázares-Ramírez, R.I., François, J.L., Martin-del-Campo, C., 2017. On the stability of fractional neutron point kinetics (FNPK). Appl. Math. Model. 45, 505–515. Espinosa-Paredes, G., Gallegos, E.D.V., Núñez-Carrera, A., Polo-Labarrios, M.-A., Espinosa-Martínez, E.-G., Vázquez-Rodríguez, A., 2014. Fractional neutron point

9

kinetics equation with Newtonian temperature feedback effects. Prog. Nucl. Energy 73, 96–101. Espinosa-Paredes, G., Polo-Labarrios, M.-A., 2012. Time-Fractional telegrapher´s equation (P1) approximation for the transport equation. Nucl. Sci. Eng. 171, 258–264. Espinosa-Paredes, G., Polo-Labarrios, M.-A., Espinosa-Martínez, E.-G., Del ValleGallegos, E., 2011. Fractional neutron point kinetics equations for nuclear reactor dynamics. Ann. Nucl. Energy 38, 307–330. Glasstone, S., Sesonske, A., 2012. Nuclear reactor engineering: reactor systems engineering. Springer Science & Business Media. Gorenflo, R., Mainardi, F., 1997. Fractional calculus: integral and differential equations of fractional order. In: Carpinteri, A., Mainardi, F. (Eds.), Fractals and Fractional Calculus in Continuum Mechanics. Springer, Berlin. Hamada, Y.M., 2017. Modified fractional neutron point kinetics equations for finite and infinite medium of bar reactor core. Ann. Nucl. Energy 106, 118–126. Hamada, Y.M., Brikaa, M.G., 2017. Applied Mathematical Modelling. Ann. Nucl. Energy 102, 359–367. Hamieh, S.D., Saidinezhad, M., 2012. Analytical solution of the point reactor kinetics equations with temperature feedback. Ann. Nucl. Energy 42, 148–152. Henry, A.F., 1975. Nuclear-Reactor Analysis. Cambrige, Massachusetts and London. Hetrick, D.L., 1993. Dynamics of Nuclear Reactor. Illinois, USA. Hilfer, R., 2000. Applications of fractional calculus in physics. World Scientific, Singapore. Li, H., Chen, W., Luo, L., Zhu, Q., 2009. A new integral method for solving the point reactor neutron kinetics equations. Ann. Nucl. Energy 36, 427–432. Li, H., Chen, W., Zhang, F., Luo, Z., 2007. Approximate solutions of point kinetics equations with one delayed neutron group and temperature feedback during delayed supercritical process. Ann. Nucl. Energy 34, 521–526. Nahla, A.A., 2017. Analytical solution of the fractional point kinetics equations with multi-group of delayed neutrons during start-up of a nuclear reactor. Ann. Nucl. Energy 99, 247–252. Nahla, A.A., 2011. An efficient technique for the point reactor kinetics equations with Newtonian temperature feedback effects. Ann. Nucl. Energy 38, 2810– 2817. Nahla, A.A., 2009. An analytical solution for the point reactor kinetics equations with one group of delayed neutrons and the adiabatic feedback model. Prog. Nucl. Energy 51, 124–128. Nahla, A.A., Hemeda, A., 2017. Picard iteration and Padé approximations for stiff fractional point kinetics equations. Appl. Math. Comput. 293, 72–80. Nobrega, J., 1971. A new solution of the point kinetics equations. Nucl. Sci. Eng. 46, 366–375. Nowak, T.K., Duzinkiewicz, K., Piotrowski, R., 2015. Numerical solution analysis of fractional point kinetics and heat exchange in nuclear reactor. Nucl. Eng, Des, p. 281. Nowak, T.K., Duzinkiewicz, K., Piotrowski, R., 2014a. Fractional neutron point kinetics equations for nuclear reactor dynamics – numerical solution investigations. Ann. Nucl. Energy 73, 317–329. Nowak, T.K., Duzinkiewicz, K., Piotrowski, R., 2014b. Numerical solution of fractional neutron point kinetics model in nuclear reactor. Arch. Control Sci. 24, 129–154. Palma, D.A.P., Martinez, A.S., Gonçalves, A.C., 2009. Analytical solution of point kinetics equations for linear reactivity variation during the start-up of a nuclear reactor. Ann. Nucl. Energy 36, 1469–1471. Podlubny, I., 1998. Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications. Academic Press, New York. Polo-Labarrios, M.A., Quezada-García, S., Espinosa-Paredes, G., Franco-Pérez, L., Ortiz-Villafuerte, J., 2020. Novel numerical solution to the fractional neutron point kinetic equation in nuclear reactor dynamics. Ann. Nucl. Energy 137, 107173. Ray, S.S., 2015. Fractional calculus with applications for nuclear reactor dynamics. CRC Press. Ray, S.S., Patra, A., 2015. On the solution of the nonlinear fractional neutron pointkinetics equation with newtonian temperature feedback reactivity. Nucl. Technol. 189, 103–109. Ray, S.S., Patra, A., 2013. Numerical solution of fractional stochastic neutron point kinetic equation for nuclear reactor dynamics. Ann. Nucl. Energy 54, 154–161. Ray, S.S., Patra, A., 2012. An explicit finite difference scheme for numerical solution of fractional neutron point kinetic equation. Ann. Nucl. Energy 41, 61–66. Schramm, M., Petersen, C.Z., Vilhena, M.T., Bodmann, B.E.J., Alvim, A.C.M., 2013. On the fractional neutron point kinetics equations. In: Integral Methods in Science and Engineering. Springer, New York. Stacey, W.M., 2001. Nuclear Reactor Physics. John Wiley and Sons, Inc., USA. Vyawahare, V.A., Nataraj, P.S.V., 2013a. Development and analysis of some versions of the fractional-order point reactor kinetics model for a nuclear reactor with slab geometry. Commun. Nonlinear Sci. Numer. Simul. 18, 1840–1856. Vyawahare, V.A., Nataraj, P.S.V., 2013b. Fractional-order modeling of neutron transport in a nuclear reactor. Appl. Math. Model. 37, 9747–9767. Vyawaharea, V.A., Espinosa-Paredes, G., 2018. On the stability of linear fractionalspace neutron point kinetics (F-SNPK) models for nuclear reactor dynamics. Ann. Nucl. Energy 111, 12–21. Zhang, F., Chen, W.Z., Gui, X.W., 2008. Analytic method study of point-reactor kinetic equation when cold start-up. Ann. Nucl. Energy 35, 746–749.