Modeling Neutron Dynamics in Nuclear Reactor using Fractional-order

Point Reactor Kinetics Model with Adiabatic Temperature Feedback

Vishwesh A. Vyawahare

and P. S. V. Nataraj

Department of Electronics Engineering, Ramrao Adik Institute of Technology, Nerul, Navi Mumbai, India

IDP in Systems and Control Engineering, Indian Institute of Technology Bombay, Powai, Mumbai, India

Keywords:

Nuclear Reactor, Point Reactor Kinetics Model, Dynamic Model, Fractional-order Modeling.

Abstract:

This paper deals with the development and analysis of fractional-order (FO) point reactor kinetics (FPRK)

model with reactivity feedback for a nuclear reactor. Incorporation of adiabatic temperature feedback of re-

activity makes this model nonlinear. It basically forms a system of coupled, nonlinear ordinary differential

equations. The nonlinear subprompt critical FPRK model is developed and analyzed in detail. Fundamental

motivation for this model is the fact that neutron transport inside the core of a reactor is truly a subdiffusion.

The work presented here analyzes the effect of temperature feedback on the neutron concentration dynamics

inside reactor core which is modeled using fractional differential equations. The system of nonlinear differ-

ential equations is solved numerically. The analysis clearly establishes the fact that the proposed model is

‘stable’ in the sense that it predicts self-limitting power excursions. The model presented in this paper consti-

tutes an important step in the development of fractional-order model for a nuclear reactor, which can be used

to achieve better control and operation.

1 INTRODUCTION

The heart of a nuclear power plant is the nuclear re-

actor. In this, the heat energy is generated by carrying

out a controlled ﬁssion of nuclei of ﬁssile radioactive

materials with the help of neutrons. Fission reactions

are a result of neutrons moving inside the reactor core

and colliding with the nuclei of core material. Due

to the use of radioactive materials and high probabil-

ity of this ﬁssion chain reaction becoming uncontrol-

lable, utmost care has to be taken to design, construct,

maintain, and operate/control a nuclear reactor. In

view of this, the mathematical modeling of nuclear

reactor is a key step in designing an efﬁcient and safe

reactor. As given in (Duderstadt and Hamilton, 1976;

Glasstone and Sesonske, 2002), this reactor model is

fundamentally based on the model of neutron trans-

port in reactor core. Thus, the validity and applicabil-

ity of reactor model will depend on how perfectly one

models the neutron transport in its core.

In the classical analysis, the diffusion approxi-

mation of neutron transport is used widely. The

integer-order (IO) neutron diffusion equation (in one-

dimension), based on the Fick’s constitutive law, is

∂φ (x,t)

∂t

+ (Σ

− νΣ

)φ(x,t) = D

∂

φ(x,t)

∂x

, (1)

where v is the neutron velocity, φ(x,t) is the neutron

ﬂux at location x at time instant t, D is the diffu-

sion coefﬁcient, ν is the average number of neutrons

emitted per ﬁssion reaction, and Σ

,Σ

are the respec-

tive macroscopic cross-sections of absorption and ﬁs-

sion reactions. But the concept of modeling neutron

transport as diffusion has some problems, viz., the

diffusion model is applicable mainly in the modera-

tor of the core and should not be used to model the

neutron movements near the regions with strong ab-

sorption, next, it predicts inﬁnite speed of propaga-

tion of neutrons (Beckurts and Wirtz, 1964; Meghre-

blian and Holmes, 1960; Espinosa-Paredes et al.,

2008). In an attempt to rectify these shortcomings

and achieve a better representation of neutron move-

ments, an FO neutron telegraph equation model was

proposed in (Vyawahare and Nataraj, 2010; Vyawa-

hare and Nataraj, 2013b) (see (2)). It was developed

using the stochastic technique of continuous-timeran-

dom walk (CTRW) as given in (Compte and Metzler,

1997) for a slab reactor.

∂

2α

φ(x,t)

∂t

2α

∂

φ(x,t)

∂t

φ(x,t) = D

∂

φ(x,t)

∂x

(2)

where 0 < α < 1, and M

= τ

(Σ

− νΣ

)+ 1/v

and

= Σ

− νΣ

. The terms like

∂

∂t

denote the Caputo

352

Vyawahare V. and Nataraj P..

Modeling Neutron Dynamics in Nuclear Reactor using Fractional-order Point Reactor Kinetics Model with Adiabatic Temperature Feedback.

DOI: 10.5220/0005038103520360

In Proceedings of the 4th International Conference on Simulation and Modeling Methodologies, Technologies and Applications (SIMULTECH-2014),

pages 352-360

ISBN: 978-989-758-038-3

 2014 SCITEPRESS (Science and Technology Publications, Lda.)

fractional time-derivatives of order α (see (4)). The

FO point reactor kinetics (FPRK) model was reported

in (Vyawahare and Nataraj, 2012), which is a system

on coupled nonlinear ordinary differential equations

(ODEs):

P(t) =

ρ(t) − β

P(t) + λC(t),

C(t) =

P(t) − λC(t), (3)

with one delayed neutron group (1G). Here P(t) is

the reactor power, ρ(t) is the reactivity, Λ is the mean

generation time between the birth of neutron subse-

quence absorption inducing ﬁssion, C(t) is the aver-

age concentration of the delayed neutrons, and λ is

the average decay constant of the delayed neutrons.

Various versions of the FPRK model were reported in

(Vyawahare and Nataraj, 2013a). The FPRK model

(or the PRK model in general) forms a set of coupled,

nonlinear differential equations in reactor power P(t)

and delayed neutron concentration C(t).

Note that the Caputo fractional time-derivative

deﬁnition is considered in models (2) and (3) which

is deﬁned as follows (Samko et al., 1997). The Ca-

puto fractional derivative (FD) of order α ∈ R

of a

causal function f(t) is given by

f(t) =

Γ(n− α)

(τ)

(t − τ)

α−n+1

dτ, (4)

with n ∈ N,n− 1 < α < n, and f

(τ) is the n

-order

derivative of the function f(t). It is seen that this

deﬁnition requires f(t) to be n-times differentiable

and furthermore this derivative has to be integrable.

This condition makes deﬁnition (4) quite restrictive.

Nevertheless, it is preferred by engineers and physi-

cists because FO differential equations (FDEs) with

Caputo derivatives have same initial conditions (ICs)

(which have well-deﬁned physical meanings) as that

for the integer-order differential equations. Recently,

fractional derivatives have been extensively used for

modeling a variety of systems and processes (Das,

2011; Magin, 2006), and also in control (Monje et al.,

2010). One of the major applications of FDEs is in

modeling of anomalous diffusion occurring in com-

plex systems (Compte and Metzler, 1997). These

FDE models are found to be more realistic and com-

pact than their counterparts, the classical integer-

order models. For a detailed history and an exhaus-

tive bibliography of fractional calculus and its appli-

cations, see (Machado et al., 2011).

The point reactor kinetics model establishes de-

pendence of the neutron ﬂux or power in reactor core

on reactivity. The remarkable feature about reactor

mechanism is that the reactivity also depends on the

power. So there is an inherent feedback (negative,

in fact) present in the reactor (Duderstadt and Hamil-

ton, 1976; Hetrick, 1993). There is a kind of ‘cyclic’

mechanism related to neutron ﬂux and reactivity: re-

activity affecting power, which in turn affects the re-

activity. The justiﬁcations explaining the dependence

of reactivity on the power can be summarized simplis-

tically as

1. Reactor power depends on the reactivity.

2. Core temperature depends on the reactor power.

3. Reactivity depends on the core temperature.

In this paper, we consider this feedback mechanism

to develop FO point reactor kinetics (FPRK) model

with temperature feedback of reactivity, which mim-

ics the situation of a subprompt critical reactor sub-

jected to a small positive reactivity (ρ

< β). Litera-

ture survey reveals that there have been only two at-

tempts in which the analysis of FPRK model with re-

activity feedback is carried out. The ﬁrst of these ref-

erences (Espinosa-Paredes et al., 2014) analyzes the

FPRK model with Newtonian reactivity feedback and

uses the FPRK model developedin (Espinosa-Paredes

et al., 2011). The second contribution (Vyawa-

hare and Nataraj, 2014) reports the analysis of FO

Nordheim-Fuchs model with adiabatic temperature

feedback of reactivity. To the best of our knowledge,

this is for the ﬁrst time that the development of a sub-

prompt critical nonlinear fractional-order point reac-

tor kinetics model using adiabatic temperature feed-

back of reactivity for a nuclear reactor is reported.

The proposed FPRK model is based on the funda-

mental assumption of considering the neutron trans-

port as anomalous diffusion, particularly subdiffu-

sion (Klages et al., 2008; Metzler and Klafter, 2000).

The literature survey reveals that there have been at-

tempts to developother types of fractional-order mod-

els of the neutron transport and the nuclear reac-

tor, (Espinosa-Paredes et al., 2008; Das and Biswas,

2007; Sardar et al., 2010), (Das et al., 2011; Espinosa-

Paredes et al., 2011), (Kadem, 2009; Kadem and

Baleanu, 2010). For a detailed and rigorous review on

PRK models in general, see (Espinosa-Paredes et al.,

2011).

The paper is organized as follows. Next sec-

tion discusses in brief the inherent reactivity feed-

back mechanism present in nuclear reactor. In section

3, development of the proposed FPRK model is pre-

sented. Analysis of the proposed FO model is given

in section 4. A comparison of the results with the

IO point reactor kinetics (IPRK) model in terms of

time evolution of power, reactivity and reactor core

temperature is presented. Issues related to the use of

various solution methods (both analytical and numer-

ModelingNeutronDynamicsinNuclearReactorusingFractional-orderPointReactorKineticsModelwithAdiabatic

TemperatureFeedback

353

ical) for solving the stiff FPRK model are discussed

in detail. Conclusion is given in section 5. The ap-

pendix at the end gives the numerical algorithm used

for solving the fractional-order model.

2 REACTIVITY FEEDBACK

MECHANISM IN NUCLEAR

REACTOR

The feedback mechanism in the reactor can be repre-

sented using the block diagram in Fig. 1 (see (Duder-

stadt and Hamilton, 1976), (Hetrick, 1993)). Let the

Figure 1: Closed-loop conﬁguration with reactivity feed-

back.

reactor be represented by the PRK model. We start

with the assumption that the reactor is operating at a

steady-state equilibrium power level P

. Now there

will be two types of reactivities present in the reactor,

] ≡ feedback reactivity due to P

≡ external reactivity. (5)

The feedback reactivity ρ

] mostly corresponds to

a negativereactivity, trying to reduce the neutron ﬂux,

and ultimately the power. Hence the feedback shown

in Fig. 1 is to be considered as a negative feedback,

and is also known as the power defect in reactivity

(Duderstadt and Hamilton, 1976). If we allow this

process to continue, it will result into the gradual re-

duction in the number of ﬁssion reactions and so in

the number of neutrons. The reactor will gradually

become more and more subcritical and a time will

come when the reactor will eventually shut down. In

order to keep the reactor running and maintain its crit-

icality, an external (positive) reactivity ρ

must be ap-

plied (like withdrawal of control rods) to balance the

negative reactivity such that

ρ(t) = ρ

+ ρ

] = 0. (6)

Now let the power change to a new value P(t). The

incremental power is deﬁned as the deviation of the

power from the equilibrium value,

p(t) ≡ P(t) − P

. (7)

The corresponding incremental changes in the reac-

tivities can be expressed as

δρ

ext

(t) = ρ

ext

(t) − ρ

δρ

[p] = ρ

[P] − ρ

], (8)

where, ρ

[P] ≡ feedback reactivity due to P(t), and

ext

(t) ≡ external reactivity to counterbalance ρ

[P].

As a result, the net reactivity input to the reactor is

comprised of two components,

ρ(t) = δρ

ext

(t) + δρ

[p], (9)

which is depicted in the block diagram of Fig. 1.

3 FPRK MODEL WITH

REACTIVITY FEEDBACK

In this section, we develop the FPRK model with

feedback of reactivity. Note that we consider the sit-

uation of one delayed group (1G) only, although ex-

tending this model to the six delayed group case is

quite trivial and straightforward.

First, we consider the IPRK model with reactivity

feedback (Hetrick, 1993), which is a system of three

nonlinear IO ordinary differential equations (IDEs):

P(t) =

ρ(t) − β

P(t) + λC(t),

C(t) =

P(t) − λC(t), (10)

ρ(t) =

γ(t) − K

P(t),

where, γ(t) is the impressed reactivity, α

is the tem-

perature coefﬁcient of reactivityand K

is the recipro-

cal of the reactor heat capacity. The ICs are P

= P(0),

= C(0), and ρ

= ρ(0).

Now we derive the FPRK model with reactivity

feedback. The 1G FPRK model (3) is rewritten:

P(t) =

ρ(t) − β

P(t) + λC(t),

C(t) =

P(t) − λC(t). (11)

We now append to this system the reactivity feedback

equation. Temperature feedback for time-varying re-

activity is given by (Hetrick, 1993)

ρ(t) = γ(t)− α

(T(t) − T

), (12)

where γ(t) is the impressed time-varying reactivity

and T(t) is the reactor core temperature with T

be-

ing the initial temperature at t = 0. Next the adiabatic

model is:

T(t) = K

P(t). (13)

SIMULTECH2014-4thInternationalConferenceonSimulationandModelingMethodologies,Technologiesand

Applications

354

Differentiating (12) with respect to t, and using (13),

we get the ODE for reactivity, as

ρ(t) =

γ(t) − α

P(t). (14)

Thus we have a system of three nonlinear differential

equations: an FDE for P(t) and two ODEs for C(t)

and ρ(t) as

P(t) =

ρ(t) − β

P(t) + λC(t),

C(t) =

P(t) − λC(t), (15)

ρ(t) =

γ(t) − K

P(t).

This is the FPRK model with reactivity feedback. It

is also appended with three ICs P

, C

, and ρ

In the next section, we analyze and compare the

IO and FO models (10) and (15) for a step-type re-

activity insertion. This will give us a clear picture of

how the reactor power varies for a sudden change in

reactivity input with the negative temperature feed-

back when subdiffusive neutron transport framework

is used.

4 ANALYSIS OF THE PROPOSED

MODEL

In this section, we carry out a thorough analysis of the

proposed FPRK model (15). We make the reactor be-

low prompt critical, that is a step reactivity variation

of much smaller magnitude (ρ

< β) is impressed. It

will be seen that the negative temperature feedback of

reactivity limits the power rise, eventually bringing it

back to zero (Duderstadt and Hamilton, 1976; Het-

rick, 1993; Nahla, 2009). Since this nonlinear model

is very difﬁcult to solve analytically, we go for the

numerical solution. The Adams-Bashforth-Moulton

method, used widely in the ﬁeld of fractional calculus

and FO control was used to solve these models. The

variations in power, precursor concentration, temper-

ature, and reactivity are obtained for different values

of fractional power α. Some issues related to the con-

vergence of the numerical method are reported and

discussed.

A certain amount ρ

of positive reactivity is sud-

denly inserted into the reactor. We need to keep the re-

actor below prompt critical (Hetrick, 1993) by choos-

ing ρ

< β. Thus, γ(t) = ρ

⇒

γ(t) = 0. So the

ODE for reactivity in IPRK and FPRK models (10),

(15) becomes

ρ(t) = −α

P(t). (16)

We use the data from (Duderstadt and Hamilton,

1976): β = 0.0075, λ = 0.08 sec

−1

, Λ = 10

−3

sec.

The ICs chosen are P

= 1 watts, C

= 93.75, and

= 0.0025. We take K

= 0.05 K/MW sec and

= 5 × 10

−5

−1

. Variation in the reactor temper-

ature T(t) is also studied. It is obtained using (12)

(with γ(t) = ρ

) as

T(t) =

− ρ(t)

+ T

. (17)

As a convention it is assumed that T

= 0

◦

The IPRK system (10) was solved using the MAT-

LAB ODE solver

ode15s

suitable for the stiff ODEs

(Mathworks, 2005). This particular solver was cho-

sen because the PRK models incorporating temper-

ature feedback of reactivity, in general, form a stiff

system of nonlinear ODEs (Espinosa-Paredes et al.,

2011; Aboanber and Nahla, 2004). The step-size was

h = 1 × 10

−3

sec. Time-variation of power, delayed

neutron precursor concentration, reactivity, and tem-

perature are shown in Fig. 2.

0 50 100 150 200 250 300 350 400 450 500

time (sec)

P(t)

IPRK with reactivity feedback (step reactivity)

0 50 100 150 200 250 300 350 400 450 500

1000

2000

time (sec)

C(t)

0 50 100 150 200 250 300 350 400 450 500

−5

x 10

−3

time (sec)

ρ(t)

0 50 100 150 200 250 300 350 400 450 500

100

200

time (sec)

T(t)

Figure 2: IPRK model with reactivity feedback: plots for

P(t), C(t), ρ(t), and T(t) for a step reactivity input.

As we see, the power starts rising due to the in-

sertion of positive reactivity. However, its rate of in-

crease is much slower. This peculiar behaviour is due

to the presence of delayed neutrons as they help in

slowing down the dynamics of the reactor. This in-

crease in power causes the reactor temperature to rise.

The adiabatic negative temperature feedback shows

its effect and reactivity starts decreasing. The power

attains a peak value of P

max

= 14.63 watts at the in-

stant t = 101.75 sec. The reactivity and the precur-

sor concentration at this instant are 0.0004523 and

1289.13 respectively. Finally, the power reduces to

zero with the reactivity settling at −0.002592. The

ﬁnal core temperature is 101.83

◦

Next we solve and analyze the FPRK model (15).

The same data and ICs are considered. As it is im-

ModelingNeutronDynamicsinNuclearReactorusingFractional-orderPointReactorKineticsModelwithAdiabatic

TemperatureFeedback

355

possible to obtain a closed-form solution for the non-

linear FDE system, we opt for other techniques for its

solution. To mention explicitly, the Adomian Decom-

position Method (ADM) (Daftardar-Gejji and Jafari,

2005), and Variational Iteration Method (VIM) (Odi-

bat and Momani, 2006) were tried. These methods

as such don’t come under the category of numerical

methods, because they provide solutions in the form

of a power series with easily computable terms. These

methods are claimed to have many advantages over

the classical numerical methods, viz., no discretiza-

tion, high accuracy, minimal calculations, to name a

few. We tried to implement these methods for our

problem using Mathematica (Ruskeepaa, 2009). But

due to inherent stiff nature of the FPRK model, these

methods did not work and a convergentsolution could

not be achieved.

It was then decided to use an improved version

of the Adams-Bashforth-Moulton (ABM) algorithm

(Diethelm, 2010) which is based on the Predictor-

Corrector scheme for the FDE system (Diethelm

et al., 2002). This method worked perfectly for the

given FDE system. It should be noted that the order of

convergence for the ABM method is a non-decreasing

function of the fractional order α. Only two values

of fractional order α, 0.7, and 0.9 are considered as

we could not make the algorithm convergefor smaller

values of α. Salient steps in the ABM algorithm are

given in the Appendix.

The FPRK model (15) is solved for two values of

α and analysis is carried out in detail.

1. α = 0.7

The step-size used for the ABM method was h =

0.05 sec. The plots for P(t), C(t), ρ(t), and T(t)

are shown in Fig. 3.

0 50 100 150 200 250 300 350 400 450 500

time (sec)

P(t)

FPRK with reactivity feedback (step reactivity), α = 0.7

0 50 100 150 200 250 300 350 400 450 500

1000

2000

time (sec)

C(t)

0 50 100 150 200 250 300 350 400 450 500

−5

x 10

−3

time (sec)

ρ(t)

0 50 100 150 200 250 300 350 400 450 500

100

200

time (sec)

T(t)

Figure 3: FPRK model with reactivity feedback (α = 0.7):

plots for P(t), C(t), ρ(t), and T(t) for a step reactivity input.

We notice that the behaviour of this FPRK model

in line with the reactor dynamics. The impressed

step reactivity causes the power to shoot up, albeit

at a slower rate. The power reaches to its peak

value, P

max

= 14.2031 watts (which is less than

the P

max

of IPRK model) at t = 103.75 sec. The

reactivity at this peak value of power is 0.000454

and the corresponding precursor concentration is

1255.09. The rest of the dynamics is similar to

the IPRK model. Negative reactivity required to

bring power to zero is −0.002585. The tempera-

ture ultimately settles at 101.7043

◦

2. α = 0.9

Same step-size h = 0.05 sec. is used to solve the

nonlinear FPRK model (15) with α = 0.9. As said

earlier, the behaviour of this FPRK model is very

similar to that of the IPRK model. The effect of

imposing a positive reactivity of ρ

= 0.0025 on

power, precursor concentration, core temperature,

and the reactivity itself is depicted in Fig. 4.

0 50 100 150 200 250 300 350 400 450 500

time (sec)

P(t)

FPRK with reactivity feedback (step reactivity), α = 0.9

0 50 100 150 200 250 300 350 400 450 500

1000

2000

time (sec)

C(t)

0 50 100 150 200 250 300 350 400 450 500

−5

x 10

−3

time (sec)

ρ(t)

0 50 100 150 200 250 300 350 400 450 500

100

200

time (sec)

T(t)

Figure 4: FPRK model with reactivity feedback (α = 0.9):

plots for P(t), C(t), ρ(t), and T(t) for a step reactivity input.

The power rise as a result of achieving the sub-

prompt criticality is slow. The peak power is

14.5434 watts. It occurs at t = 101.75 sec. The

reactivity and precursor concentration at P

max

are

0.0004513 and 1282.094 respectively. The power

excursion is controlled by reduction in reactiv-

ity due to negative temperature feedback. Fi-

nally, a negative reactivity of −0.002591 brings

the power to zero with reactor temperature settling

at T = 101.8262

◦

To carry out a comparative study for the IO and FO

models, these observations are compiled in tables 1

and 2. Survey of these two tables and the plots in

Figs. 2, 3, and 4 bring following observations to our

notice:

1. The power overshoot in all models is small. In

case of FPRK models, for both α = 0.7 and 0.9,

the peak power attained was nearly equal to that

SIMULTECH2014-4thInternationalConferenceonSimulationandModelingMethodologies,Technologiesand

Applications

356

Table 1: IPRK and FPRK models with reactivity feedback:

values of various quantities at P

max

for step reactivity input.

α P

max

ρ(t) C(t)

(watts) at at P

max

at P

max

(sec)

0.7 14.2 103.7 4.5e-4 1255

0.9 14.5 101.7 4.5e-4 1282

IPRK 14.6 101.1 4.5e-4 1289

Table 2: IPRK and FPRK models with reactivity feedback:

values of various quantities at steady-state for step reactivity

input.

α P(t) C(t) ρ(t) T(t)

(watts)

◦

0.7 0.0258 2.9 −0.002585 101.7

0.9 0.0131 1.6 −0.002591 101.8

IPRK 0.0111 1.4 −0.002592 101.8

in the IPRK model. These maximum values in

power occur almost at the same instant for all

these models. Also, values of positive reactiv-

ity and precursor concentration at P

max

are almost

equal.

2. The inherent negative feedback of temperature

helps in limiting the power overshoot and even-

tually brings it down to a negligible value. This

dynamics is observed in both IO and FO models.

The delayed neutron dynamics also works for this

cause.

3. After the power excursion dies out, it settles to a

very low value (almost equal for IPRK and FPRK

models). A very little negative reactivity is re-

quired to bring this power to zero. The reactor

temperature settles at around 101

◦

C for all the

models.

4. However a different scenario is observed in the

dynamics of precursor concentration. The val-

ues of C(t) when power reaches to its peak value

are different for IPRK and FPRK models. Fur-

thermore, the values also depend on the fractional

derivative power α. It is seen that as the derivative

order increases, precursor concentration at P

max

increases (lowest for FPRK model with α = 0.7

and highest for IPRK model). However, exactly

opposite situation is observed during the steady-

state. Now, the ﬁnal value at which C(t) settles

after the power excursion has died out is largest

for FPRK model with derivativeorder 0.7. It grad-

ually decreases and is lowest for the IPRK model.

One may deduce that for the reactivity feedback case,

the behaviours of IPRK and FPRK models (with α =

0.7,0.9) nearly coincide. However, it should be noted

that a nuclear reactor is very critical and complex sys-

tem and its safe operation is a very crucial issue. A

slight variation in any parameter may have severe ef-

fect on performance of the system. In view of this, the

results predicted by the proposed FPRK model may

yield a better understanding of the nuclear reactor dy-

namics.

Furthermore, we think that the effect of consid-

ering subdiffusive neutron transport would be more

saliently visible for lower values of α. But as men-

tioned earlier, we found it almost impossible to get the

ABM method converged for smaller values of frac-

tional differentiation order. Nevertheless, the results

presented here conﬁrm the validity of the developed

FPRK model. It faithfully captures the behaviour of a

reactor subjected to a subprompt step in reactivity un-

der the inﬂuence of adiabatic temperature feedback.

Thus, in this section we analyzed the FPRK model

with reactivity feedback. As the aim was to study the

reactor with below prompt critical situation, the dy-

namics of precursor concentration was also consid-

ered in the model. The derived model was subjected

to a detailed analysis. As a part of the study, we ex-

amined the behaviour of FPRK model when subjected

to the step reactivity. The reactivity insertion was

chosen in such a way so as to keep the reactor sub-

prompt critical. This exercise conﬁrms that the devel-

oped fractional-order nonlinear model with reactivity

feedback faithfully represents the reactor dynamics.

A comparison of this FO model with the classical IO

model is also carried out and various results are pre-

sented.

5 CONCLUSIONS

To get a more reliable and realistic model of the nu-

clear reactor, it is necessary to consider processes by

which the reactivity is affected by the neutron ﬂux or

power. Using this fact, in this paper, we proposed

a new version of the fractional-order point reactor

kinetics model considering the effect of temperature

feedback on the reactivity. The adiabatic model of

temperature dependence on power is used.

The proposed FPRK model is used to analyze

small perturbations (below prompt critical) in the

power. This model is a system of three coupled, non-

linear ODEs (one FDE + two IDEs). Various stan-

dard methods like VIM, ADM, were tried. But it was

found that the stiff nature of the PRK model in gen-

eral proved as a big hurdle in achieving the conver-

gence for these methods. Finally, the ABM method

was used to numerically solve the FO model. ‘It be-

comes more and more difﬁcult to make a numerical

ModelingNeutronDynamicsinNuclearReactorusingFractional-orderPointReactorKineticsModelwithAdiabatic

TemperatureFeedback

357

method applied to an FDE converge for smaller val-

ues of α’: this frequently observed phenomenon was

experienced for the model under consideration.

The behaviour predicted by the FO model was

found to be in line with reactor physics. Each time,

the power excursion was found to be self-limiting and

therefore stable. Thus this paper presents a major step

in the development and analysis of fractional-order

model for a nuclear reactor under the consideration

of reactivity feedback. The developed FO model very

faithfully mimics the actual behaviour of the reactor

in these situations. Also the developed FO model has

broader applicability, and is easy to derive and solve.

The classical integer-order model forms a special case

of the proposed FO model.

The analysis carried out in this paper can be made

more exhaustive by studying the FPRK models with

other types of reactivity feedback mechanisms and

carrying out a comparative study of these models.

Further, it has been proved in literature that the heat

transfer mechanism is better represented using frac-

tional dynamics. Hence a more detailed study of the

reactivity feedback in nuclear reactor can be achieved

by additionally considering a fractional-order model

for the temperature dependence of reactivity.

REFERENCES

Aboanber, A. E. and Nahla, A. A. (2004). On pade’ approx-

imations to the exponential function and application to

the point kinetics equations. Progress in Nuclear En-

ergy, 44(4):347–368.

Beckurts, K. H. and Wirtz, K. (1964). Neutron Physics.

Springer-Verlag, Germany.

Compte, A. and Metzler, R. (1997). The generalized Cat-

taneo equation for the description of anomalous trans-

port processes. Journal of Physics A: Mathematical

and General, 30:7277–7289.

Connolly, J. A. (2004). The numerical solution of fractional

and distributed order differential equations. PhD the-

sis, University of Liverpool, UK.

Daftardar-Gejji, V. and Jafari, H. (2005). Adomian decom-

position: a tool for solving a system of fractional dif-

ferential equations. Journal of Mathematical Analysis

and Applications, 301:508–518.

Das, S. (2011). Functional Fractional Calculus for System

Identiﬁcation and Controls. Springer, Germany.

Das, S. and Biswas, B. B. (2007). Fractional divergence for

neutron ﬂux proﬁle in nuclear reactor. International

Journal of Nuclear Energy Science and Technology,

3(2):139–159.

Das, S., Das, S., and Gupta, A. (2011). Fractional or-

der modeling of a PHWR under step-back condition

and control of its global power with a robust PI

controller. IEEE Transactions on Nuclear Science,

58(5):2431–2441.

Diethelm, K. (2010). The Analysis of Fractional Differ-

ential Equations: An Application-Oriented Exposi-

tion Using Differential Operators of Caputo Type.

Springer, Germany.

Diethelm, K., Ford, N. J., and Freed, A. D. (2002). A

Predictor-Corrector approach for the numerical solu-

tion of fractional differential equations. Nonlinear Dy-

namics, 29:3–22.

Diethelm, K., Ford, N. J., Freed, A. D., and Luchko, Y.

(2005). Algorithms for the fractional calculus: A se-

lection of numerical methods. Computer Methods in

Applied Mechanics and Engineering, 194:743–773.

Duderstadt, J. J. and Hamilton, L. J. (1976). Nuclear Reac-

tor Analysis. John Wiley & Sons, USA.

Espinosa-Paredes, G., del Valle-Gallegos, E., N´u¯nez-

Carrera, A., Polo-Labarrios, M.-A., Espinosa-

Mart´ınez, E.-G., and V´azquez-Rodr´ıguez, R. (2014).

Fractional neutron point kinetics equation with new-

tonian temperature feedback effects. Progress in Nu-

clear Energy, 73:96–101.

Espinosa-Paredes, G., Morales-Sandoval, J. B., V´azquez-

Rodr´ıguez, R., and Espinosa-Mart´ınez, E.-G. (2008).

Constitutive laws for the neutron transport current.

Annals of Nuclear Energy, 35:1963–1967.

Espinosa-Paredes, G., Polo-Labarrios, M.-A., Espinosa-

Mart´ınez, E.-G., and del Valle-Gallegos, E. (2011).

Fractional neutron point kinetics equations for nuclear

reactor dynamics. Annals of Nuclear Energy, 38:307–

330.

Glasstone, S. and Sesonske, A. (2002). Nuclear Reactor

Engineering: Vol. 1. CBS Publishers & Distributors,

India.

Hetrick, D. L. (1993). Dynamics of Nuclear Reactors.

American Nuclear Society, USA.

Kadem, A. (2009). The fractional transport equation: an

analytical solution and a spectral approximation by

Chebyshev polynomials. Applied Sciences, 11:78–90.

Kadem, A. and Baleanu, D. (2010). Analytical method

based on Walsh function combined with orthogonal

polynomial for fractional transport equation. Commu-

nications in Nonlinear Science and Numerical Simu-

lation, 15(3):491–501.

Klages, R., Radons, G., and Sokolov, I. M., editors (2008).

Anomalous Transport. WILEY-VCH Verlag GmbH &

Co.

Li, C. and Peng, G. (2004). Chaos in Chen’s system

with a fractional order. Chaos, Solitons & Fractals,

22(2):443–450.

Machado, J. T., Kiryakova, V., and Mainardi, F. (2011).

Recent history of fractional calculus. Communica-

tions in Nonlinear Science and Numerical Simulation,

16(3):1140–1153.

Magin, R. L. (2006). Fractional Calculus in Bioengineer-

ing. Begell House Publishers, USA.

Mathworks (2005). MATLAB Manual. The Mathworks

Inc., MATLAB version 7.1 (R14), USA.

Meghreblian, R. V. and Holmes, D. K. (1960). Reactor

Analysis. McGraw-Hill Book Company, USA.

SIMULTECH2014-4thInternationalConferenceonSimulationandModelingMethodologies,Technologiesand

Applications

358

Metzler, R. and Klafter, J. (2000). The random walk’s guide

to anomalous diffusion: A fractional dynamics ap-

proach. Physics Reports, 339:1–77.

Monje, C. A., Chen, Y. Q., Vinagre, B. M., Xue, D.,

and Feliu, V. (2010). Fractional-order Systems and

Control: Fundamentals and Applications. Springer-

Verlag London Limited, UK.

Nahla, A. A. (2009). An analytical solution for the point

reactor kinetics equations with one group of delayed

neutrons and the adiabatic feedback model. Progress

in Nuclear Energy, 51:124–128.

Odibat, Z. and Momani, S. (2006). Application of Varia-

tional Iteration Method to nonlinear differential equa-

tions of fractional order. International Journal of Non-

linear Sciences and Numerical Simulation, 1(7):15–

27.

Ruskeepaa, H. (2009). Mathematica Navigator: Mathemat-

ics, Statistics and Graphics. Academic Press, USA.

Samko, S. G., Kilbas, A. A., and Marichev, O. I. (1997).

Fractional Integrals and Derivatives. Gordon and

Breach Science Publishers, Netherlands.

Sardar, T., Ray, S. S., Bera, R., Biswas, B., and Das, S.

(2010). The solution of coupled fractional neutron dif-

fusion equations with delayed neutrons. International

Journal of Nuclear Energy Science and Technology,

5(2):105–113.

Tavazoei, M. S. and Haeri, M. (2007). A necessary condi-

tion for double scroll attractor existence in fractional-

order systems. Physics Letters A, 367:102–113.

Vyawahare, V. A. and Nataraj, P. S. V. (2010). Modeling

neutron transport in a nuclear reactor as subdiffusion:

The neutron fractional-order telegraph equation. In

The 4th IFAC Workshop on Fractional Differentiation

and its Applications, Badajoz, Spain.

Vyawahare, V. A. and Nataraj, P. S. V. (2012). Develop-

ment and analysis of the fractional point reactor kinet-

ics model for a nuclear reactor with slab geometry. In

The 5th IFAC Workshop on Fractional Differentiation

and its Applications, Nanjing, China.

Vyawahare, V. A. and Nataraj, P. S. V. (2013a). Develop-

ment and analysis of some versions of the fractional-

order point reactor kinetics model for a nuclear reac-

tor with slab geometry. Communications in Nonlinear

Science and Numerical Simulation, 18:1840–1856.

Vyawahare, V. A. and Nataraj, P. S. V. (2013b). Fractional-

order Modeling of Neutron Transport in a Nuclear

Reactor. Applied Mathematical Modelling, 37:9747–

9767.

Vyawahare, V. A. and Nataraj, P. S. V. (2014). Devel-

opment and analysis of fractional-order Nordheim-

Fuchs model for nuclear reactor. In Daftardar-Gejji,

V., editor, Fractional Calculus: Theory and Applica-

tions. Narosa Publishing House, India.

APPENDIX: FRACTIONAL

SECOND-ORDER

ADAMS-BASHFORTH-MOULTON

METHOD

Here the fractional second- order Adams-Bashforth-

Moulton (ABM) method which is used in section 4 is

explained in brief. The main computational steps in-

volved in the algorithm are presented here for the eq-

uispaced grid points. For details, refer to (Diethelm,

2010; Diethelm et al., 2005; Connolly, 2004). It is

an extension of the classical ABM method used to

numerically solve the ﬁrst-order ODEs. It comes in

the category of the so-called PECE (Predict, Evaluate,

Correct, Evaluate) type since it involves calculation of

the predictor value which is in turn used to compute

the corrector value. This method and its variants are

very popular in the ﬁeld of fractional calculus and ap-

plied areas (Li and Peng, 2004; Tavazoei and Haeri,

2007). The algorithm explained below is for a sin-

gle fractional differential equation. However, it can

be easily modiﬁed to handle a system of FDEs.

Consider the single term FDE with Caputo FD

y(t) = f(t,y(t)), (18)

where α ∈ R

and with the appropriate initial condi-

tions:

y(0) = y

(k)

, k = 0,1,...,m− 1, (19)

where, m = ⌈α⌉ is the ceil function. The equivalent

Volterra integral equation is

y(t) =

m−1

∑

k=0

y(0)+

Γ(α)

(t−τ)

α−1

f(τ,y(τ))dτ

(20)

The integration limits from 0 to t imply the nonlocal

structure of the fractional derivatives.

The next step is to use the product trapezoidal

quadrature formula to replace the integral in (20). We

approximate the following integral

k+1

− τ)

α−1

g(τ)dτ, (21)

k+1

− τ)

α−1

k+1

(τ)dτ, (22)

where ˜g

k+1

≡ piecewise linear interpolation for g(t)

with grid points at t

, j = 0,1,2, . . . , k + 1. Thus we

can write the integral (22) as

k+1

−τ)

α−1

k+1

(τ)dτ =

k+1

∑

j=0

j,k+1

g(t

), (23)

ModelingNeutronDynamicsinNuclearReactorusingFractional-orderPointReactorKineticsModelwithAdiabatic

TemperatureFeedback

359

for the equispaced nodes (t

= jh with some ﬁxed

step-size h). The values of a

j,k+1

are given for j = 0

α(α+ 1)



α+1

− (k− α)(k+ 1)



for 1 ≤ j ≤ k as



α(α+ 1)



(d),

where

d = (k− j+ 2)

α+1

+ (k− j)

α+1

− 2(k− j+ 1)

α+1

and for j = k + 1 as

α(α+ 1)

So the corrector formula is

k+1

m−1

∑

j=0

k+1

( j)

Γ(α)

∑

j=0

j,k+1

f(t

) (24)

Γ(α)



k+1,k+1

f(t

k+1

)



where now the predictor y

k+1

is evaluated as

k+1

m−1

∑

j=0

k+1

( j)

Γ(α)

∑

j=0

j,k+1

f(t

), (25)

with

j,k+1

((k+ 1− j)

− (k− j)

). (26)

For 0 < α < 1, the predictor and corrector expressions

get modiﬁed as

k+1

= y

Γ(α)

∑

j=0

j,k+1

f(t

), (27)

and

k+1

= y

Γ(α)

∑

j=0

j,k+1

f(t

)

Γ(α)



k+1,k+1

f(t

k+1

)



.(28)

As already mentioned in section 4, the convergence of

this algorithm deteriorates as α → 0. This algorithm

was coded in MATLAB.

SIMULTECH2014-4thInternationalConferenceonSimulationandModelingMethodologies,Technologiesand

Applications

360