Numerical solution analysis of fractional point kinetics and heat exchange in nuclear reactor

Nuclear Engineering and Design 281 (2015) 121–130

Contents lists available at ScienceDirect

Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes

Numerical solution analysis of fractional point kinetics and heat exchange in nuclear reactor Tomasz Karol Nowak ∗ , Kazimierz Duzinkiewicz, Robert Piotrowski Gdansk University of Technology, Faculty of Electrical and Control Engineering, Narutowicza 11/12, 80-233 Gdansk, Poland

a r t i c l e

i n f o

Article history: Received 18 June 2014 Received in revised form 5 November 2014 Accepted 7 November 2014

a b s t r a c t The paper presents the neutron point kinetics and heat exchange models for the nuclear reactor. The models consist of a nonlinear system of fractional ordinary differential and algebraic equations. Two numerical algorithms are used to solve them. The first algorithm is application of discrete Grünwald–Letnikov definition of the fractional derivative in the model. The second involves building an analog scheme in the FOMCON Toolbox in MATLAB environment where Oustaloup recursive filter is defined. The simulation results for input step reactivity are discussed and compared. The impact of the model parameters is examined. © 2014 Elsevier B.V. All rights reserved.

1. Introduction The neutron kinetics takes an important role in the nuclear reactor. With respect to time of appearance, neutrons are divided into immediate and delayed neutrons. The immediate neutrons (more than 99% of all fission neutrons) are emitted directly in the fission act, while the delayed neutrons (less than 1%) are produced during the radioactive decay of certain fission products – delayed neutron precursors. Research works related to modeling of neutron point kinetics have a long history (Duderstadt and Hamilton, 1976; Hetrick, 1993; Li et al., 2009; Barati and Setayeshi, 2010; Hamieh and Saidinezhad, 2012). In most cases the classic bilinear model is used. The analytical (Nahla, 2010) and numerical (Kim et al., 2014) solution for this model were obtained. In this paper, in opposite of previous works, nonlinear system of fractional and ordinary differential equations for neutron point kinetics model is the subjects of research. In (Espinosa-Paredes et al., 2011) the fractional point kinetics (FPK) model was introduced. It consists of fractional derivative and one group of delayed neutron precursors. The authors extended the model to six groups (Nowak et al., 2014a) and normalized it (Nowak et al., 2014b).

Fractional calculus generalizes traditional calculus to noninteger differential and integral orders. Defining fractional derivatives or integrators in the models may be proposed as a better approximation of many processes. Some of the fractional models were presented in (Petráˇs, 2011b; Podlubny, 1999). The authors’ research concentrate on the nuclear reactor analysis, therefore the FPK model is the subject of their interests. In the paper, the FPK model for six groups of delayed neutron precursors was extended to heat exchange (HE) model. The HE takes place between the core and the coolant. The role of the coolant in the reactor is played by water flowing in the space between bars with the fissile material. The paper is organized as follows. Firstly, the definition of fractional derivative is introduced. Then the FPK model with six groups of delayed neutron precursors and HE model were described and normalized, along with the assumed values of particular constant parameters and the adopted initial conditions. After that, methods to solve the FPK and HE models were presented and described. Finally, test results, including comparative analysis, were illustrated. 2. The mathematical models 2.1. The fractional point kinetics model

∗ Corresponding author. Tel.: +48 583471742; fax: +48 583472487. E-mail addresses: [email protected] (T.K. Nowak), [email protected] (K. Duzinkiewicz), [email protected] (R. Piotrowski). http://dx.doi.org/10.1016/j.nucengdes.2014.11.028 0029-5493/© 2014 Elsevier B.V. All rights reserved.

There are two main approaches for defining a fractional derivative. The first one considers differentiation and integration as limits of finite differences (Grünwald–Letnikov definition). The second approach generalizes a convolution type representation of repeated integration (Riemann–Liouville and Caputo definitions). To solve

122

T.K. Nowak et al. / Nuclear Engineering and Design 281 (2015) 121–130

the fractional derivative models, different analytical and numerical algorithms are used. Their detailed review can be found in Yang (2010). In the paper the Grünwald–Letnikov approach is used. It has the following form (Petráˇs, 2011b; Podlubny, 1999):

t − a k a Dt f (t)

= lim

1

h→0 hk

h 

(−1)p

p=0

p

f (t − ph)

dk+1

d  n(t) +  k n(t) + dt dt k+1 =

6 



j=1

k

dk dt k

Cj (t) +

6 



dk+1 dt k+1

1 1−ˇ + l 



dk dt k

n(t) +

6 

ˇ − (t) n(t) 

j Cj (t)

(2)

(3)

6

ˇ, neutrons of the jth group in the population of neutrons, ˇ = j=1 j j is the decay constant of delayed-neutron precursors, l is the immediate neutron lifetime, (t) is the reactivity,  k is the relaxation time, and  is the effective lifetime of the neutrons. Neutron density was examined in the paper. It is proportional to nuclear reactor’s power (Duderstadt and Hamilton, 1976; Duzinkiewicz et al., 1989). Therefore, the neutron density will be normalized further in the paper for purposes of testing nuclear reactor’s power. Reactivity is the relative deviation of the relation between the sizes of successive generations of neutrons. If nuclear reactor power is greater than 0 and if reactivity is negative – the power decreases (and also neutron density), if positive – the power increases and if equal to 0 – the power is at a constant level (Duderstadt and Hamilton, 1976; Hetrick, 1993; Li et al., 2009; Hamieh and Saidinezhad, 2012; Duzinkiewicz et al., 1989). Relaxation time is the parameter specifying the rate of system return to the state of equilibrium and is time-constant. Greater values of this parameter mean that the system (in this case: the density of neutrons in the nuclear reactor core) returns to the state of equilibrium slower. The impact of relaxation time on the neutron density of the FPK model (without HE) was examined in Nowak et al. (2014a). The fractional derivative and relaxation time will be examined. The discussion of results will be further in the paper. The analysis was performed using the normalized values of n(t) and Cj (t). The normalized values of neutron density are denoted by n, whilst the normalized concentrations of delayed neutron precursors nuclei of the jth group are denoted by Cj . In order to perform normalization, the normalized quantities and the relation between them are defined: Cj =

Cj (t) Cj,B

,

Cj,B =

ˇj j

nB , j = 1, 2, . . ., 6

ˇj

Cj +

1 1−ˇ + l  6  ˇj



dk dt k

n+

ˇ − (t) 

d C = j n − j Cj for j = 1, . . ., 6 dt j

(6)

dt k j

j=1



The critical state of equilibrium was adopted as the initial state. In such case, the ordinary and fractional derivatives are equal zero (Hilfer, 1995). Therefore the initial conditions for this state were: n0 = Cj0 = 1.

The point heat exchange model can be written in a system of differential equations (Duzinkiewicz et al., 1989): Mfe cep

where n0 = n(0) and Cj0 = Cj (0) are initial conditions, n(t) is the density of neutrons, Cj (t) is the concentration of delayed neutron precursor nuclei of the jth group, ˇj is the fraction of the delayed

n(t) , nB

dk



(5)

Mc cc

n=

k

d n + k dt

2.2. The heat exchange model

j=1

ˇj d C (t) = n(t) − j Cj (t) for j = 1, . . ., 6 dt j 

n+

Cj

j=1

(1)

where k is the order of fractional derivative, h is step-size, a and t are the bounds relating to the above operation. The FPK model with six groups of delayed neutron precursors is given by fractional differential equations (Nowak et al., 2014a,b): k

k

=

  k

where nB and Cj,B are base values of above quantities. The FPK model after normalization is given by formula:

(4)

1 dTfe (t) = Q˙ r (t) − (Tfe (t) − Tc (t)) R dt

(7)

dTc (t) 1 = (Tfe (t) − Tc (t)) − 2w(t)cp (Tc (t) − Tin (t)) R dt

(8)

where Tfe (t) – averaged temperature of the fuel element core volume, Tc (t) – averaged temperature of the refrigerant, Tin (t) – coolant temperature at the inlet to the core, Mfe – total mass of the fuel elements in the core volume, Mc – total mass of the coolant and components in the core volume, cfe – weighted average specific heat of the fuel assemblies in the core, cc – weighted average specific heat of the coolant (at constant pressure) and materials of the core, cp – coolant’s specific heat (at constant pressure), R – core’s thermal resistance, w – coolant mass flow rate (at constant pressure), Q˙ r – thermal nuclear’s reactor power. In order to obtain Tfe (0) = 1 and Tc (0) = 1, the HE model has to be normalized in the same way as the FPK model. Obtained relations for normalized core temperature Tfe and coolant temperature Tc are as follows: Tfe =

Tfe (t) , Tfe,B

Tc =

Tc (t) Tc,B

(9)

where Tfe,B – base value of core temperature, Tc,B – base value of coolant temperature. Their adoption allows to normalize model (7) and (8): Tfe − (Tc,B /Tfe,B )Tc dTfe Q˙ r = − Tfe,B Mfe cfe RMfe cfe dt 1 dTc = RMc cc dt



Tfe,B T − Tc Tc,B fe



−

2wcp Mc cc

(10)

 Tc −

Tin Tc,B

 (11)

The heat exchange model takes into account reactivity effects – the variations of reactivity caused by change of core and coolant temperature. The reactivity effect can be written as follows (Duzinkiewicz et al., 1989):  = (t) − 0 = ˛c Tc,B (Tc − 1) + ˛fe Tfe,B (Tfe − 1)

(12)

where: ˛c – reactivity coefficient of coolant temperature change, ˛fe – reactivity coefficient of core temperature change, 0 – reactivity input. Neutron density n is proportional to nuclear reactor’s power Q˙ r (Duderstadt and Hamilton, 1976; Duzinkiewicz et al., 1989): Q˙ r = n where:  – proportional coefficient.

(13)

T.K. Nowak et al. / Nuclear Engineering and Design 281 (2015) 121–130

Hence, in FPK and HE models can be used both n or Q˙ r .

By substituting the right hand side of (19) into Dk+1 n, Dn and Dk n in (14):

3. Two approaches to solve the nonlinear fractional differential system of equations

i 

1 h1+k

The FPK and HE models consist of nonlinear systems of fractional and ordinary differential and algebraic equations. Nonlinearity is caused by multiplying neutron density and reactivity in (3). Reactivity is an algebraic function of core and coolant temperature. Fractional derivative is present in the neutron point kinetic equation. The solution of the FPK model was obtained applying the numerical solution of linear multi-term fractional differential system of equations proposed in (Edwards et al., 2002). It was used in the previous papers (Nowak et al., 2014a,b). That method could not be applied, because the considered model is nonlinear (see (5) and (12)). The Adomian decomposition method (Jafari and DaftardarGejji, 2006) was attempted to solve this model. The method turned out to be ineffective. Computed absolute values of coefficients to solve the model are increasing every iteration for specified problem. Eventually, they were too large to define in PC. That is why other methods to solve the FPK and HE model were proposed. The first entails substituting the discrete Grünwald–Letnikov definition of fractional derivative into the FPK and HE derivatives (Petráˇs, 2011a). The second is to use the FOMCON Toolbox in MATLAB to build an analog scheme in Simulink, where the system of equations could be solved (Tepljakov et al., 2011). The following symbols were introduced to simplify the notation: Dk instead of dk /dtk , D instead of d/dt, Dk+1 instead of dk+1 /dtk+1 and  instead of (t). After applying the new nomenclature, (5), (6) and (10)–(13) become as follows: D

k+1

=

n + a3 Dn + a2 D n + a11 n + a12 n b2j Dk Cj +

6 

j=1

b1j Cj

(14)

j=1

DCj = a0j n − a0j Cj for j = 1, . . ., 6

(15)

DTfe = m1 n + m2 Tfe + m3 Tc

(16)

DTc = m4 Tfe + m5 Tc + F

(17)

 = 0 + ˛c Tc,B (Tc − 1) + ˛fe Tfe,B (Tfe−1 )

(18)

According to (5), (6) and (10)–(12) the coefficients take the folˇ 1−ˇ 1 lowing form: a0j = j , a11 = 1k  , a12 = 1k  , a2 = 1l +  , a3 = 1k , b1j = m4 =

ˇj k 

ˇj , j

, b2j =

Tfe,B 1 , RMc cc Tc,B



m5 = −



1 , Tfe,B Mfe cfe

m1 =

(1/R)+2wcp , Tc,B Mc cc

m2 =

F=

1 , RMfe cfe

2wcp Tin . Tc,B Mc cc

m3 =

 Tc,B /Tfe,B , RMfe cfe

The above symbols

will be used both in the both methods.

  (−1)p

p=0

k

p

f (ti−p ) = h−k

hk

p=0

 6

(k)

cp Cj,i−p +

p=0

(21)

b1j Cj,i

j=1

To get the neutron density in the ith iteration one can write it on the left hand side of the equation, while the remaining part of the equation on the right hand side. It results in:

6

i

b2j j=1 h

ni =

1 1

h 1+k

+

+

1 h1+k

i



i (1+k) (1) a c ni−p + h3 c ni−p p=1 p p=1 p a a 1 + h3 + 2k + a12 i + a11 h1+k h

+ 1 h1+k



6 (k) c cj,i−p + b C n p=0 p j=1 1j j,i i−p a3 a2 + h + k + a12 i + a11 h

+

a2 hk a3 h

i +

(k) c ni−p p=1 p a2 + a12 i + a11 hk

(22)

Using the same transformation, values for other quantities in ith iteration can be calculated from following equations as: Cj,i =

Tc,i =

a0j nj −

1 h

i

(1) c Cj,i−p p=1 p

1 h

− b2j

m1 ni + m3 Tc,i − 1 h

m4 Tfe,i + F − 1 h

1 h

1 h

for j = 1, . . ., 6

(23)

i

(1) c Tfe,i−p p=1 p

(24)

− m2

i

(1) c Tc,i−p p=1 p

(25)

− m5

i = 0 + ˛c Tc,B (Tc,i − 1) + ˛fe Tfe,B (Tfe,i − 1)

(26)

It is worth noting that both on the left and right hand side of the above equations are quantities in ith iteration. i can be found on the right hand side of (22). In (23) it is ni , in (24) – ni and Tc,i , in (25) – Tfe,i and in (26) – Tep,i and Tc,i . Because it is impossible to calculate the value of one quantity in ith iteration while not having the value of another quantity in this iteration, the mentioned quantities on the right hand side of the equations in ith iteration were replaced by i − 1th iteration. The application of discrete Grünwald–Letnikov definition of the fractional derivative is hereinafter referred to as the first method. In figures it appears as GL.

In order to make a clear presentation of an analog scheme, the auxiliary variables were defined:

To numerically solve nonlinear fractional derivative equation in ith iteration one can use the discrete expression of (1) (Petráˇs, 2011a): i 

 b2j  i

j=1

i

3.2. Application of FOMCON toolbox

3.1. Application of discrete Grünwald–Letnikov definition of the fractional derivative

Dk f (t) ≈ h−k

ni−p +

p=0

6

=

a3  (1) a2  (k) cp ni−p + cp ni−p +(a11 +a12 i )ni h hk i

(1+k)

cp

p=0

Tfe,i =

k

6 

123

i 

(k)

cp f (ti−p )

(19)

p=0

x1 = n, x2 = Dk n = Dk x1 , x3 = Dn = Dx1 , x4 = Dk+1 n = Dk+1 x1 = Dk x3

y1j = Cj ,

y2j = Dk Cj = Dk y1j ,

(27)

y3j = DCj = Dy1j

for

j = 1, . . ., 6

(k)

(28)

where: ti = ih and cp (p = 1, 2, . . ., i) – binomial coefficients that can be written as (Petráˇs, 2011a): (k)

c0 = 1,

(k)

cp =



1−

1+k p



(k)

cp−1

(20)

z11 = Tfe , z12 = DTfe = Dz11 , z21 = Tc , z22 = DTc = Dz21 , z3 = 

(29)

124

T.K. Nowak et al. / Nuclear Engineering and Design 281 (2015) 121–130

integrator. For sk where 0 < k < 1 it is defined by (Tepljakov et al., 2011):

It results from (14)–(18) and (27)–(29) that: 3 

x4 = −

ai xi − a11 x1 − a12 z3 x1 +

i=2

6 2  

bij yij

(30) sk = Gf (s) = K

i=1 j=1

(32)

(33)z22 = m4 z11 + m5 z21 + F z3 = 0 + ˛fe Tfe,B (z11 − 1) + ˛c Tc,B (z21 − 1)

k

l=−N

(31)y3j = a0j x1 − b2j y1j for j = 1, . . ., 6 z12 = m1 x1 + m2 z11 + m3 z21

N

s + ω

(34)

Eqs. (27)–(34) were modeled as an analog scheme (see Fig. 1). The presented analog scheme was built using blocks from the FOMCON Toolbox. Those blocks are fractional derivative sk and fractional integrator s−k . They consist of the Oustaloup recursive filter that gives a very good approximation of fractional derivative and

(35)

s + ωk

where N is the approximation order and the ωk , ωk , K are obtained from: ωk = ωb ωk = ωb

 ω (l+N+(1/2)(1−k))/2N+1 h

ωb

 ω (l+N+(1/2)(1−k))/2N+1 h

ωb

, ,

K = whk

(36)

where (ωb , ωh ) is the specified frequency range. The default frequency range used in toolbox blocks is (ωb , ωh ) = (10−3 , 103 ), whilst the default approximation order is N = 5.

Fig. 1. FPK and HE models – analog scheme.

T.K. Nowak et al. / Nuclear Engineering and Design 281 (2015) 121–130

Default values returned results inconsistent with basic laws of physics taking place in a nuclear reactor. For this reason, the frequency range was changed to (ωb , ωh ) = (10−1 , 1). The transfer function for the Oustaloup recursive filter was calculated for two values of k (0.25, 0.75) for which the tests were

125

performed. For all mentioned k: K = 1. The transfer functions of Oustaloup recursive filter are: s0.25 =

(s + 0.1082)(s + 0.1334)· · ·(s + 0.8774) (s + 0.1140)(s + 0.1405)· · ·(s + 0.9245)

Fig. 2. Normalized neutron density for 0 = 0.0005,  k = 10−4 sk and k = 0.25.

Fig. 3. Normalized neutron density for 0 = 0.0005,  k = 10−4 sk and k = 0.25.

Fig. 4. Normalized neutron density for 0 = 0.0005,  k = 10−4 sk and k = 0.75.

(37)

126

T.K. Nowak et al. / Nuclear Engineering and Design 281 (2015) 121–130

s0.75 =

(s + 0.1027)(s + 0.1266)· · ·(s + 0.8326) (s + 0.1201)(s + 0.1481)· · ·(s + 0.9742)

(38)

The Oustaloup recursive filter for s−k , which is the fractional integration of order k, can be obtained from: s−k =

1 sk

=

1 Gf (s)

(39)

It is suggested to use ode23tb solver (Tepljakov et al., 2011), which is the implicit Runge–Kutta method. It is a two-stage method. Its first stage is a trapezoidal rule. The second is a second-order backward differentiation formula. The application of FOMCON Toolbox was hereinafter referred to as the second method. In figures it appears as FOMCON. 4. Results and discussion The numerical results were obtained for FPK and HE models. Two methods to solve nonlinear fractional models were applied. The obtained results were shown for input step reactivity 0 = 0.0005, selected values of relaxation time  k (10−4 sk and

10−5 sk ), fractional order of derivative k (0.25, 0.75, 0.96 and 0.99) and step-size h (0.025 s, 0.05 s and 0.1 s). The results for order of derivative equal 0.96 and 0.99 were shown in order to compare them with the results obtained in (Espinosa-Paredes et al., 2011). In (Espinosa-Paredes et al., 2011) the fractional neutron point kinetics model with one group of delayed-neutron precursors was examined for k equal 0.96, 0.97, 0.98 and 0.99. In the paper, two extreme values were used. It should be noted, that in the paper the FPK model was built for six groups of delayed-neutrons precursors and that the HE model was added. The discussion was made. The following parameters were adopted in (5)–(6): ˇj = [0.000231, 0.001533, 0.001327, 0.002765, 0.000805, 0.000294], ˇ = 0.07, j = [0.0124, 0.0305, 0.111, 0.301, 1.13, 3] (Kinard and Allen, 2004),  = l = 0.00003, whilst in (10)–(13): Mfe = 60,645 kg, R = 2.325 × 10–7 J/◦ C K, Mc = 10,573 kg, cc cfe = 313.07 J/kg K, = 3402.9 J/kg K, cp = 5460 J/kg K, Q˙ r = 1375 MW, Tin = 269 ◦ C, Tfe,B = 603.16 ◦ C, Tc,B = 283.47 ◦ C, ˛c = –1.2 × 10–4 ◦ C−1 , ˛fe = –3.4 × 10–5 ◦ C−1 , w = 8700 kg/s,  = Q˙ r /nB 1375 × 106 /122.42 × 106 = 11.232 W cm−3 . The initial state of the reactor was assumed by the critical state, which is state of equilibrium. For the first test, the results are presented for simulation time equal

Fig. 5. Normalized neutron density for 0 = 0.0005,  k = 10−4 sk and k = 0.96.

Fig. 6. Normalized neutron density for 0 = 0.0005,  k = 10−4 sk and k = 0.99.

T.K. Nowak et al. / Nuclear Engineering and Design 281 (2015) 121–130

to 150 s (which show the steady state) and 5 s (which show the biggest changes taking place in the first seconds after reactivity step) (see Figs. 2–6). Other results are illustrated only for time 5 s (see Figs. 7 and 8). Normalized neutron density for ␶k = 10−4 sk is shown in Figs. 2–6. The neutron density responses can be divided into two areas. The first one is neutron density increment caused by positive step reactivity. The second one is neutron density decrement caused by reactivity effects. Reactivity effect, caused by changes in core and coolant temperature, affects the reactivity. This results in reactivity decrement to 0. It causes the stabilization of neutron density. The neutron density is supposed to be greater after stabilization than its initial value (before tests) (Duzinkiewicz et al., 1989). The neutron density progress was different for tested values of fractional derivative. The higher value of k, the maximum neutron density value is higher. This relationship can be seen for both methods. Higher k causes faster intersection of neutron density for both methods.

127

The results for k equal 0.96 and 0.99 are quantitatively the same as for other examined values of k. Comparing the results with obtained in (Espinosa-Paredes et al., 2011) it can be written that the reactivity effects causes stabilizing of neutron density and that the value of k does not have the impact on the nature of the results. Neutron density when applying the discrete Grünwald– Letnikov definition increases slower than when using the FOMCON Toolbox, but it also decreases at a lower rate. After about one second neutron density obtained for first method is greater than for the second. In Figs. 7 and 8 neutron density for  k = 10−5 sk is presented (Figs. 9 and 10). . There is a difference between results obtained for  k = 10−5 sk and  k = 10−4 sk . The smaller  k , the smaller the difference between first and second method. It can be explained by the fact, that  k → 0 causes change in (2), which becomes classic point kinetic

Fig. 7. Normalized neutron density for 0 = 0.0005,  k = 10−5 sk and k = 0.25.

Fig. 8. Normalized neutron density for 0 = 0.0005,  k = 10−5 sk and k = 0.75.

128

T.K. Nowak et al. / Nuclear Engineering and Design 281 (2015) 121–130

model described in Duderstadt and Hamilton (1976). The parts of equation containing fractional derivative are becoming no longer relevant. Due to this fact, the difference between results for tested value of fractional derivative is smaller for  k = 10−5 sk than for  k = 10−4 sk . The neutron density stabilization for results obtained with  k = 10−5 sk can be explained in the same way as it was with  k = 10−4 sk . Obtained results for positive input reactivity show that concerned models, where the fractional derivative is defined in FPK, may be used to approximate neutron kinetics and heat exchange in nuclear reactor. The results are consistent with basic laws of physics taking place in nuclear reactor dynamics. The impact of selected parameters was examined. 5. Comparative analysis There is no analytical solution for nonlinear system of fractional and ordinary differential and algebraic equations consisting of (5),

(6) and (10)–(13). In order to check which of two methods is more accurate for FPK and HE model, a comparison to a classic model was drawn. The classic model does not consist of fractional derivatives. After normalization the equation for neutron density is given by Duderstadt and Hamilton (1976):

 ˇj d (t) − ˇ n(t) = n(t) + C (t) dt   j 6

(40)

j=1

The above equation can be obtained by assuming  k = 0 in (5), but one has to have in mind that (5) is proposed extension of (40) and that the above equation was derived earlier than the one consisting of fractional derivatives. The classic model also consists of (6) and (10)–(13). In this paper, model in (40) is considered to be classic model. The difference between results obtained for the classic model and results obtained for the first and second method will be compared.

Fig. 9. Normalized neutron density for 0 = 0.0005,  k = 10−5 sk and k = 0.96.

Fig. 10. Normalized neutron density for 0 = 0.0005,  k = 10−5 sk and k = 0.99.

T.K. Nowak et al. / Nuclear Engineering and Design 281 (2015) 121–130

129

Fig. 11. Difference between neutron density obtained by classic model and by first and second method for 0 = 0.0005 and  k = 0.

Fig. 12. Difference between neutron density obtained by classic model and by first and second method for 0 = 0.0005 and  k = 0.

For the first method the neutron density in ith iteration for  k = 0 was obtained from (40):

6 ni =



i (1) 1 c ni−p p=1 p h ˇ− 1 + i h

C − j=1 j j,i

(41)

Equations to solve other quantities were already derived (see (23)–(26)). The model in Fig. 1 was simplified for the second method where  k = 0. The fractional integrators, derivatives and blocks multiplied by them were removed. The integrator of order 1 was replaced by two integrators and fractional derivative. The order of the fractional derivative block was set to 1. It is not possible to set initial conditions in FOMCON blocks. Therefore, it was not possible to replace the integrator with a fractional integrator. The initial condition was set in an integrator block. The difference between classic model and the first or second method is shown in Figs. 11 and 12. Fig. 12 shows the difference within the first 0.5 s of simulation. It was caused by the existence of the highest differences at the beginning of the test. In order to verify which one of the methods is more accurate to approximate FPK and HE models, the comparison was made. Three indicators were used. The first one was the maximum value of the difference between the classic model and the first or the second method (MAX). The second one was the integral of this difference

Table 1 Indicators to compare first and second method. Indicator

GLh = 0.1 a

GLh = 0.05 a

GLh = 0.025 a

FOMCONb

MAX ID IAD

0.069893 −0.010043 0.020701

0.062261 −0.0066771 0.011621

0.050105 −0.0050268 0.0073061

0.017888 −0.0033438 0.0034594

a The difference between results obtained from the classic model and the first method. b The difference between results obtained from the classic model and the second method.

(ID). The third one was the absolute value of this difference (IAD). The values of the indicators are presented in Table 1. It can be seen that the values of indicators for the second method are the most satisfying. As it could be expected, the value of stepsize for the first method determines the values of indicators. 6. Conclusions The numerical solution of FPK and HE models were analyzed in the paper. Results obtained by two methods for typical input were shown and compared using three proposed indicators. The first method was application of the discrete Grünwald–Letnikov definition of the fractional derivative. The second method was application

130

T.K. Nowak et al. / Nuclear Engineering and Design 281 (2015) 121–130

of the FOMCON Toolbox. The impact of relaxation time, fractional order and step-size were examined. Obtained results were consistent with the basic laws of physics taking place in nuclear reactor. It leads to conclusion that both methods can be used to approximate FPK and HE models. The comparison of methods showed that the second method is quantitatively better. There is a big question for researchers what values should be assumed by relaxation time, fractional order and step-size in order to obtain a satisfactory approximation for the nuclear reactor processes. Acknowledgements This work was supported by the National Centre for Research and Development under Strategic Research Project No. SP/J/10/176450/12. The authors wish to express their thanks for the support. References Barati, R., Setayeshi, S., 2010. A model for nuclear research reactor dynamics. Nucl. Eng. Des. 262 (9), 251–263. Duderstadt, J.J., Hamilton, L.J., 1976. Nuclear Reactor Analysis. John Wiley & Sons, New York. Duzinkiewicz, K., Baum, G., Michalak, A., 1989. Simulation Model of Basic Dynamic Processes in Nuclear Power Plant WWER based on Lumped Parameters. Technical Report. Electric Power and Control Engineering Institute, Gdansk University of Technology (in Polish). Edwards, J.T., Ford, N.J., Simpson, A.Ch., 2002. The numerical solution of linear multiterm fractional differential equations: systems of equations. J. Comput. Appl. Math. 148 (2), 401–418. Espinosa-Paredes, G., Polo-Labarrios, M.-A., Espinosa-Martinez, E.-G., del ValleGallegos, E., 2011. Fractional neutron point kinetics equations for nuclear reactor dynamics. Ann. Nucl. Energy 38 (2–3), 307–330. Hamieh, S.D., Saidinezhad, M., 2012. Analytical solution of the point reactor kinetics equations with temperature feedback. Ann. Nucl. Energy 42 (4), 148–152. Hetrick, D.L., 1993. Dynamics of Nuclear Reactors. American Nuclear Society, La Grange Park, III, USA. Hilfer, I., 1995. Foundation of fractional dynamics. Fractals 3 (3), 549–556. Jafari, H., Daftardar-Gejji, V., 2006. Solving a system of nonlinear fractional differential equations using Adomian decomposition. J. Comput. Appl. Math. 196 (2), 644–651.

Kim, H.T., Park, Y., Kazantzis, N., Parlos, A.G., Vista, I.V., Chong, F.P.K.T., 2014. A numerical solution to the point kinetic equations using Taylor–Lie series combined with a scaling and squaring technique. Nucl. Eng. Des. 272 (6), 1–10. Kinard, M., Allen, E.J., 2004. Efficient numerical solution of the point kinetics equations in nuclear reactor dynamics. Ann. Nucl. Energy 31 (9), 1039–1051. 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 (4), 427–432. Nahla, A.A., 2010. Analytical solution to solve the point reactor kinetics equations. Nucl. Eng. Des. 240 (6), 1622–1629. 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 in nuclear reactor. Arch. Control Sci. 24 (2), 129–154. Petráˇs, I., 2011a. Fractional derivatives, fractional integrals, and fractional differential equations in MATLAB. In: Assi, A.H. (Ed.), Engineering Education and Research using MATLAB. InTech, Rijeka, pp. 239–264. Petráˇs, I., 2011b. Fractional-order nonlinear systems. In: Modelling, Analysis and Simulation. Springer, New York. Podlubny, I., 1999. Fractional Differential Equations. Academic Press, San Diego. Tepljakov, A., Petlenkov, E., Belikov, J., 2011. FOMCON: a MATLAB toolbox for fractional-order system identification and control. Int. J. Microelectron. Comput. Sci. 2 (2), 51–62. Yang, Q., (PhD thesis) 2010. Novel Analytical and Numerical Methods for Solving Fractional Dynamical Systems. Queensland University of Technology, Australia. Tomasz Karol Nowak is a Ph.D. student in the Faculty of Electrical and Control Engineering (Gdansk University of Technology, Poland). He received M.Sc. from Automation and Robotics at the Gdansk University of Technology in 2011. His research interests include: mathematical modeling, fractional calculus, linear and nonlinear dynamical systems, decision support systems. Kazimierz Duzinkiewicz received M.Sc. degree in Electrical Engineering and Ph.D. in Control Engineering from Faculty of Electrical and Control Engineering at Gdansk University of Technology, in 1973 and 1982, respectively. He has been employed as a university teacher starting his work in 1973 from the post of Assistant to the current position of Assistant Professor in Department of Control Engineering. His research interests include: production scheduling and operational control of technological systems with switchable processes (refinery type systems), computer control of electric power station in emergency conditions, safety and reliability analysis of hazardous systems, mathematical modeling of complex systems, multihorizon and multilevel optimization and control structures and algorithms with applications to petroleum industry. Robert Piotrowski is a researcher in the Faculty of Electrical and Control Engineering (Gdansk University of Technology, Poland). He received M.Sc. and Ph.D. from Automation and Robotics at the Gdansk University of Technology in 2001 and 2005, respectively. His research interests include: mathematical modeling, identification, analysis and control design of nonlinear dynamical systems and design of computer control systems.