Modeling Neutron Dynamics in Nuclear Reactor using Fractional-order
Point Reactor Kinetics Model with Adiabatic Temperature Feedback
Vishwesh A. Vyawahare
1
and P. S. V. Nataraj
2
1
Department of Electronics Engineering, Ramrao Adik Institute of Technology, Nerul, Navi Mumbai, India
2
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 fission of nuclei of fissile 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 fission 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 efficient 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
1
v
φ (x,t)
t
+ (Σ
a
νΣ
f
)φ(x,t) = D
2
φ(x,t)
x
2
, (1)
where v is the neutron velocity, φ(x,t) is the neutron
flux at location x at time instant t, D is the diffu-
sion coefficient, ν is the average number of neutrons
emitted per fission reaction, and Σ
a
,Σ
f
are the respec-
tive macroscopic cross-sections of absorption and fis-
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 infinite 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.
τ
α
v
α
2α
φ(x,t)
t
2α
+M
1
α
φ(x,t)
t
α
+M
2
φ(x,t) = D
2
φ(x,t)
x
2
(2)
where 0 < α < 1, and M
1
= τ
α
(Σ
a
νΣ
f
)+ 1/v
α
and
M
2
= Σ
a
νΣ
f
. 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
Copyright
c
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):
d
α
dt
α
P(t) =
ρ(t) β
Λ
P(t) + λC(t),
d
dt
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 fission, 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
definition is considered in models (2) and (3) which
is defined as follows (Samko et al., 1997). The Ca-
puto fractional derivative (FD) of order α R
+
of a
causal function f(t) is given by
0
D
α
t
f(t) =
1
Γ(n α)
Z
t
0
f
n
(τ)
(t τ)
αn+1
dτ, (4)
with n N,n 1 < α < n, and f
n
(τ) is the n
th
-order
derivative of the function f(t). It is seen that this
definition requires f(t) to be n-times differentiable
and furthermore this derivative has to be integrable.
This condition makes definition (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-defined 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 flux 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 flux and reactivity: re-
activity affecting power, which in turn affects the re-
activity. The justifications 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 (ρ
0
< β). 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 first 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 first 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 configuration 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
0
. Now there
will be two types of reactivities present in the reactor,
ρ
f
[P
0
] feedback reactivity due to P
0
,
ρ
0
external reactivity. (5)
The feedback reactivity ρ
f
[P
0
] mostly corresponds to
a negativereactivity, trying to reduce the neutron flux,
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 fission 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 ρ
0
must be ap-
plied (like withdrawal of control rods) to balance the
negative reactivity such that
ρ(t) = ρ
0
+ ρ
f
[P
0
] = 0. (6)
Now let the power change to a new value P(t). The
incremental power is defined as the deviation of the
power from the equilibrium value,
p(t) P(t) P
0
. (7)
The corresponding incremental changes in the reac-
tivities can be expressed as
δρ
ext
(t) = ρ
ext
(t) ρ
0
,
δρ
f
[p] = ρ
f
[P] ρ
f
[P
0
], (8)
where, ρ
f
[P] feedback reactivity due to P(t), and
ρ
ext
(t) external reactivity to counterbalance ρ
f
[P].
As a result, the net reactivity input to the reactor is
comprised of two components,
ρ(t) = δρ
ext
(t) + δρ
f
[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):
d
dt
P(t) =
ρ(t) β
Λ
P(t) + λC(t),
d
dt
C(t) =
β
Λ
P(t) λC(t), (10)
d
dt
ρ(t) =
d
dt
γ(t) K
C
α
T
P(t),
where, γ(t) is the impressed reactivity, α
T
is the tem-
perature coefficient of reactivityand K
C
is the recipro-
cal of the reactor heat capacity. The ICs are P
0
= P(0),
C
0
= C(0), and ρ
0
= ρ(0).
Now we derive the FPRK model with reactivity
feedback. The 1G FPRK model (3) is rewritten:
d
α
dt
α
P(t) =
ρ(t) β
Λ
P(t) + λC(t),
d
dt
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) T
0
), (12)
where γ(t) is the impressed time-varying reactivity
and T(t) is the reactor core temperature with T
0
be-
ing the initial temperature at t = 0. Next the adiabatic
model is:
d
dt
T(t) = K
C
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
d
dt
ρ(t) =
d
dt
γ(t) α
T
K
C
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
d
α
dt
α
P(t) =
ρ(t) β
Λ
P(t) + λC(t),
d
dt
C(t) =
β
Λ
P(t) λC(t), (15)
d
dt
ρ(t) =
d
dt
γ(t) K
C
α
T
P(t).
This is the FPRK model with reactivity feedback. It
is also appended with three ICs P
0
, C
0
, and ρ
0
.
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 (ρ
0
< β) 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 difficult to solve analytically, we go for the
numerical solution. The Adams-Bashforth-Moulton
method, used widely in the field 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 ρ
0
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 ρ
0
< β. Thus, γ(t) = ρ
0
d
dt
γ(t) = 0. So the
ODE for reactivity in IPRK and FPRK models (10),
(15) becomes
d
dt
ρ(t) = α
T
K
C
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
0
= 1 watts, C
0
= 93.75, and
ρ
0
= 0.0025. We take K
C
= 0.05 K/MW sec and
α
T
= 5 × 10
5
K
1
. Variation in the reactor temper-
ature T(t) is also studied. It is obtained using (12)
(with γ(t) = ρ
0
) as
T(t) =
ρ
0
ρ(t)
α
T
+ T
0
. (17)
As a convention it is assumed that T
0
= 0
C.
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
0
10
time (sec)
P(t)
IPRK with reactivity feedback (step reactivity)
0 50 100 150 200 250 300 350 400 450 500
0
1000
2000
time (sec)
C(t)
0 50 100 150 200 250 300 350 400 450 500
−5
0
5
x 10
−3
time (sec)
ρ(t)
0 50 100 150 200 250 300 350 400 450 500
0
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
final core temperature is 101.83
C.
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
0
10
20
time (sec)
P(t)
FPRK with reactivity feedback (step reactivity), α = 0.7
0 50 100 150 200 250 300 350 400 450 500
0
1000
2000
time (sec)
C(t)
0 50 100 150 200 250 300 350 400 450 500
−5
0
5
x 10
−3
time (sec)
ρ(t)
0 50 100 150 200 250 300 350 400 450 500
0
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
C.
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
= 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
0
10
time (sec)
P(t)
FPRK with reactivity feedback (step reactivity), α = 0.9
0 50 100 150 200 250 300 350 400 450 500
0
1000
2000
time (sec)
C(t)
0 50 100 150 200 250 300 350 400 450 500
−5
0
5
x 10
−3
time (sec)
ρ(t)
0 50 100 150 200 250 300 350 400 450 500
0
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
C.
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
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)
C
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 final 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 confirm 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 influence 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 confirms 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 flux 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 difficult 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
Identification and Controls. Springer, Germany.
Das, S. and Biswas, B. B. (2007). Fractional divergence for
neutron flux profile 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
λ
D
µ
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 first-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 field 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 modified to handle a system of FDEs.
Consider the single term FDE with Caputo FD
0
D
α
t
y(t) = f(t,y(t)), (18)
where α R
+
and with the appropriate initial condi-
tions:
D
k
t
y(0) = y
(k)
0
, k = 0,1,...,m 1, (19)
where, m = α is the ceil function. The equivalent
Volterra integral equation is
y(t) =
m1
k=0
t
k
k!
D
k
t
y(0)+
1
Γ(α)
Z
t
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
Z
t
k+1
0
(t
k+1
τ)
α1
g(τ)dτ, (21)
as
Z
t
k+1
0
(t
k+1
τ)
α1
g
k+1
(τ)dτ, (22)
where ˜g
k+1
piecewise linear interpolation for g(t)
with grid points at t
j
, j = 0,1,2, . . . , k + 1. Thus we
can write the integral (22) as
Z
t
k+1
0
(t
k+1
τ)
α1
g
k+1
(τ)dτ =
k+1
j=0
a
j,k+1
g(t
j
), (23)
ModelingNeutronDynamicsinNuclearReactorusingFractional-orderPointReactorKineticsModelwithAdiabatic
TemperatureFeedback
359
for the equispaced nodes (t
j
= jh with some fixed
step-size h). The values of a
j,k+1
are given for j = 0
as
h
α
α(α+ 1)
k
α+1
(k α)(k+ 1)
α
,
for 1 j k as
h
α
α(α+ 1)
(d),
where
d = (k j+ 2)
α+1
+ (k j)
α+1
2(k j+ 1)
α+1
,
and for j = k + 1 as
h
α
α(α+ 1)
.
So the corrector formula is
y
k+1
=
m1
j=0
t
j
k+1
j!
y
( j)
0
+
1
Γ(α)
k
j=0
a
j,k+1
f(t
j
,y
j
) (24)
+
1
Γ(α)
a
k+1,k+1
f(t
k+1
,y
P
k+1
)
,
where now the predictor y
P
k+1
is evaluated as
y
P
k+1
=
m1
j=0
t
j
k+1
j!
y
( j)
0
+
1
Γ(α)
k
j=0
b
j,k+1
f(t
j
,y
j
), (25)
with
b
j,k+1
=
h
α
α
((k+ 1 j)
α
(k j)
α
). (26)
For 0 < α < 1, the predictor and corrector expressions
get modified as
y
P
k+1
= y
0
+
1
Γ(α)
k
j=0
b
j,k+1
f(t
j
,y
j
), (27)
and
y
k+1
= y
0
+
1
Γ(α)
k
j=0
a
j,k+1
f(t
j
,y
j
)
+
1
Γ(α)
a
k+1,k+1
f(t
k+1
,y
P
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