**1.3 System's variables of nuclear reactor**

There are several controlling factors in nuclear reactors such as:


The Theoretical Simulation of a Model by SIMULINK for

Where:

and:

Where:

and:

**1.4 Six factors coefficients** 

The effective multiplying coefficient is [5]:

η υ

As the thermal fission coefficient ( )

Surveying the Work and Dynamical Stability of Nuclear Reactors Cores 181

<sup>1</sup> ( ) () ( ) ( ) *Y s G s C sI A B*

> <sup>1</sup> ( ) () *sI A s* ϕ

*Gs C sB* () () = ϕ

The effective multiplying coefficient is: the ratio of generated neutrons in every generation to generated neutrons in last generation. So to operate the nuclear reactor in steady state, this parameter should be: 1 means the generated neutrons in every generation are equal with neutrons which have absorbed or leaked in last generation that means: critical state.

.... . .... *FNL THNL TNL* = =

235 235 235 238 238 238 . .. .. .. .

*F O O fm aa aa a*

235

238

*U U*

235 235 235 238 238 238

*M rM r M*

.. .. . .. .. . . . *<sup>M</sup>*

<sup>Σ</sup> + + = = Σ +Σ + ++ + (21)

238

*m A*. *<sup>N</sup>*

*m A*. *<sup>N</sup>*

. . 22 2

*F O O a aa aa a F M OO HH O O a a aa aa a a a*

 σ

*N gN gN <sup>f</sup> N gN gN N N* σ

*U mA mA N N*

== =

235 235 235 238 238 238

*O U*

σ

Also the thermal absorption coefficient (*f*) is [5]:

235

*<sup>f</sup> p P P fp P eff*

*f f f*

235

238

υ

 ηε

235 235 235

*N g N gN gN*

 σ

*K* (16)

σσ

*<sup>M</sup>* <sup>=</sup> , (18)

*<sup>M</sup>* <sup>=</sup> (19)

(20)

 σ

<sup>Σ</sup> = = Σ ++ (17)

235 238

σσ

σσ

(1 )

+ −

The minimum value of Keff is: 0 and maximum of it is: υ namely: 2.43.

ηis [5, 6]:

*F*

ηε

σ

<sup>−</sup> = =− (13)

<sup>−</sup> − = (14)

(15)

*u s*


If each system variable as a mathematics variable is considered then can write:

$$
\dot{x}(t) = Ax(t) + Bu(t) \tag{5}
$$

and:

$$y(t) = \mathbb{C}x(t) \tag{6}$$

Where:

*x(t)* is: variable of system, *x t* ( ) is: derivative of system's variable, *y(t)* is: output, *A* is: system's matrix, B is: control's matrix, C is: matrix of output and *u(t)* is: control's variable.

By taking the laplace conversions from two sides of above equations can write:

$$sX(\mathbf{s}) = AX(\mathbf{s}) + Bu(\mathbf{s})\tag{7}$$

and:

$$Y(\mathbf{s}) = \mathbb{C}X(\mathbf{s})\tag{8}$$

So two last equations that are based on Laplace conversions can be converted to following form:

$$(sI - A)X(s) = Bu(s) \tag{9}$$

and:

$$Y(\mathbf{s}) = \mathbb{C}X(\mathbf{s})\tag{10}$$

Therefore:

$$X(\mathbf{s}) = \left(\mathbf{s}I - A\right)^{-1} B u(\mathbf{s})\tag{11}$$

and:

$$Y(s) = \mathcal{C}(sI - A)^{-1}Bu(s) \tag{12}$$

In order to define the transfer function, it will be deduced as shown below:

$$\mathbf{C}(\mathbf{s}) = \frac{\mathbf{Y}(\mathbf{s})}{\mu(\mathbf{s})} = \mathbf{C}(\mathbf{s}I - A)^{-1}B \tag{13}$$

Where:

180 Nuclear Reactors

*x(t)* is: variable of system, *x t* ( ) is: derivative of system's variable, *y(t)* is: output, *A* is: system's matrix, B is: control's matrix, C is: matrix of output and *u(t)* is: control's variable.

So two last equations that are based on Laplace conversions can be converted to following

*x t Ax t Bu t* () () () = + (5)

*y*() () *t Cx t* = (6)

*sX s AX s Bu s* () () () = + (7)

*Y s CX s* () () = (8)

( ) () () *sI A X s Bu s* − = (9)

*Y s CX s* () () = (10)

<sup>1</sup> *X s sI A Bu s* () ( ) () <sup>−</sup> = − (11)

<sup>1</sup> *Y s C sI A Bu s* () ( ) () <sup>−</sup> = − (12)

If each system variable as a mathematics variable is considered then can write:

By taking the laplace conversions from two sides of above equations can write:

In order to define the transfer function, it will be deduced as shown below:

9. Fission poisons like Xe and Sm.

11. Burn up. 12. Power demand. 13. The kind of fuel. 14. Energy of neutrons. 15. Doppler's effect. 16. Value of β.

and:

Where:

and:

form:

and:

and:

Therefore:

10. Fission fragments and fission products.

$$\left(\mathrm{sI} - \mathrm{A}\right)^{-1} = \mathfrak{g}(\mathrm{s})\tag{14}$$

and:

$$G(s) = \mathbb{C}\varphi(s)B\tag{15}$$

#### **1.4 Six factors coefficients**

The effective multiplying coefficient is: the ratio of generated neutrons in every generation to generated neutrons in last generation. So to operate the nuclear reactor in steady state, this parameter should be: 1 means the generated neutrons in every generation are equal with neutrons which have absorbed or leaked in last generation that means: critical state. The minimum value of Keff is: 0 and maximum of it is: υ namely: 2.43.

The effective multiplying coefficient is [5]:

$$K\_{\rm eff} = \mathfrak{η}.f.p.\varepsilon.P\_{\rm FNL}.P\_{\rm THNL} = \mathfrak{η}.f.p.\varepsilon.P\_{\rm TNL} \tag{16}$$

As the thermal fission coefficient ( ) ηis [5, 6]:

$$\eta = \upsilon \frac{\boldsymbol{\Sigma}\_f^F}{\boldsymbol{\Sigma}\_{fn}^F} = \frac{\upsilon \boldsymbol{\nu} N^{235} \boldsymbol{\sigma}\_f^{235} \cdot \mathbf{g}\_f^{235}}{N^{235} \boldsymbol{\sigma}\_a^{235} \cdot \mathbf{g}\_a^{235} + N^{238} \boldsymbol{\sigma}\_a^{238} \cdot \mathbf{g}\_a^{238} + N^O \boldsymbol{\sigma}\_a^O} \tag{17}$$

Where:

$$N^{235} = \frac{m^{235} \, A}{M^{235}} \, \text{\,\, \,} \tag{18}$$

$$\mathcal{N}^{238} = \frac{m^{238} \mathcal{A}}{\mathcal{M}^{238}} \tag{19}$$

and:

$$N^{\circlearrowright} = 2N^{\circlearrowright} = 2\frac{m^{\circlearrowright} \cdot A}{M^{\circlearrowright}} = 2\frac{m^{\circlearrowright} \cdot A}{rM^{\circlearrowright} + (1 - r)M^{\circlearrowright}} \tag{20}$$

Also the thermal absorption coefficient (*f*) is [5]:

$$f = \frac{\boldsymbol{\Sigma}\_a^F}{\boldsymbol{\Sigma}\_a^F + \boldsymbol{\Sigma}\_a^M} = \frac{N^{235}.\sigma\_a^{235}.\boldsymbol{g}\_a^{225} + N^{238}.\sigma\_a^{238}.\boldsymbol{g}\_a^{238} + N^O.\sigma\_a^O}{N^{225}.\sigma\_a^{235}.\sigma\_a^{238}.\sigma\_a^{238} + N^O.\sigma\_a^O + N^H.\sigma\_a^H + N^{O\_M}.\sigma\_a^O} \tag{21}$$

The Theoretical Simulation of a Model by SIMULINK for

( ) ( ). *<sup>T</sup> G G*

β

therefore once a signal with delay time ( *<sup>d</sup>*

**1.5 Neutron point kinetics (NPK)**

and:

around 1 cm.

τ

Where:

and:

neutrons, is given by:

1 1

+

reactivity. When 0 *<sup>k</sup>*

 τ

τ

Surveying the Work and Dynamical Stability of Nuclear Reactors Cores 183

Firstly the transfer function is supposed. Secondly due to mostly reactor's cores are twin

The treatment of the neutron transport as a diffusion process has only been validated. For example, in a Light Water Reactor (LWR) the mean free path of thermal neutrons is typically

The fractional model has been derived for the NPK equations with n groups of delayed

*τ* is the relaxation time, *k* is the anomalous diffusion order (for sub-diffusion process: 0 < *k* < 1; while that for super-diffusion process: 1 < *k* < 2), *n* is the neutron density, Ci is the concentration of delayed neutron precursor, l is the prompt-neutron lifetime for finite media, K is the neutron generation time, β is the fraction of delayed neutrons, and ρ is the

The fractional model includes three additional terms relating to the classic equations which

1 1 *k k d n dt* +

> *k k d n*

*k i k d C*

The physical meaning of above terms suggests that for sub-diffusion processes, the first term has an important contribution for rapid changes in the neutron density (for example in the turbine trip in a BWR nuclear power plant (NPP)), while the second term represents an important contribution when the changes in the neutron density is almost slow, for example during startup in a NPP that involves operational maneuvers due to movement of control rod mechanism. The importance of third term is when the reactor sets in shutdown state, it

 ρρ

= ∇ (32)

) from first part to second part is transmitted then

 λ

<sup>+</sup> , (34)

*dt* (35)

*dt* (36)

1 1

*i i*

− − + + += + + Λ Λ , 0 2 < ≤ *<sup>k</sup>* (33)

λτ

ρ

τ

1 (1 ) *k k m m <sup>k</sup> K K k i k k ii i k*

*d n d n dn d C n C dt l dt dt dt*

<sup>+</sup> = =

→ , the classic NPK equation is recovered.

are contained of fractional derivatives (Gilberto and Espinosa, 2011):

 ρβ

that signal with a same delay time from second part to first part will be transmitted.

The resonance escape probability for fast neutrons (*p*) also is calculated as following [5]:

$$p = e^{-\left(\frac{N^{288}}{\zeta'\Delta\_\*}\right)\times I\_{\mathbb{CP}}}\tag{22}$$

Where:

$$I\_{\rm eff} = 3.9 \times \left(\frac{\Sigma\_s}{N^{238}}\right)^{0.415} \text{ \AA} \tag{23}$$

$$
\xi = \frac{A}{A + \frac{2}{3}} \tag{24}
$$

and:

$$
\Sigma\_s = N^{235}.\sigma\_s^{235} + N^{238}.\sigma\_s^{238} + N^0.\sigma^O + N^H.\sigma^H + N^{O\_M}.\sigma^O + N^{Zr}.\sigma^{Zr} \tag{25}
$$

If the enrichment of fuel is 100% then the resonance escape probability for fast neutrons (*p*) will be maximum value.

According to enrichment of applied fuel and its mass can write [7]:

$$r = \frac{m\_{fl}}{m\_f} \,\text{\,\,\,}\tag{26}$$

$$f\_{f\_{fm}} = \frac{m\_f}{m\_{fm}} = \frac{rM\_{ff} + (1 - r)M\_{nf}}{rM\_{ff} + (1 - r)M\_{nf} + M\_{O\_2}} \tag{27}$$

Where [8]:

$$m\_{l1O\_2} = m\_{fm} = \frac{N\_{fm}M\_{fm}}{A} = \rho\_{fm}.V\_{fm} \,\, \, \, \, \tag{28}$$

$$m\_{\rm LI} = m\_f = m\_{fm}.f\_{fm} = \rho\_{fm}.V\_{fm}.\frac{rM\_{ff} + (1 - r)M\_{nf}}{rM\_{ff} + (1 - r)M\_{nf} + M\_{O\_2}},\tag{29}$$

and:

$$m\_{\rm ul} \coloneqq m\_{\rm ff} = m\_{\rm fm}, f\_{\rm fm}, r = \rho\_{\rm fm} V\_{\rm fm}, \frac{rM\_{\rm ff} + (1 - r)M\_{\rm nf}}{rM\_{\rm ff} + (1 - r)M\_{\rm nf} + M\_{\rm O\_2}}. \tag{30}$$

If the transfer function of G(s) as *V* function is assumes then:

$$\dot{G}(\rho) = \frac{\partial G(\rho)}{\partial \rho\_1} \frac{d\rho\_1}{dt} + \frac{\partial G(\rho)}{\partial \rho\_2} \frac{d\rho\_2}{dt} + \dots + \frac{\partial G(\rho)}{\partial \rho\_n} \frac{d\rho\_n}{dt} \tag{31}$$

and:

182 Nuclear Reactors

238 . *eff s <sup>N</sup> <sup>I</sup>*

238

2 3

σσ

*N* Σ = ×

*A*

+

*A*

235 235 238 238 . . .. . . *OO HH OM O Zr Zr* Σ*<sup>s</sup>* = + ++ + + *N N NN N N*

If the enrichment of fuel is 100% then the resonance escape probability for fast neutrons (*p*)

*ff f*

*f ff nf fm ff nf O*

*m rM r M <sup>f</sup> m rM r M M*

*UO fm fm fm N M m m V <sup>A</sup>* == =

*m m mf V rM r M M* ρ

*m mm f r V r*

1 2

*Gd Gd G <sup>d</sup> <sup>G</sup>*

∂∂ ∂ = + ++ ∂∂ ∂

() () ( ) ( ) *<sup>n</sup>*

 ρρ

ρ

1 2

ρρ

ρρ

(1 ) . .. (1 )

(1 )

. . *fm fm*

ρ

(1 ) .. .. . (1 )

*dt dt dt*

*m r*

(1 ) *fm*

*s s*

0.415

<sup>Σ</sup> <sup>=</sup> (22)

, (23)

 σ

*<sup>m</sup>* <sup>=</sup> , (26)

, (28)

2

2

2

*ff nf*

*rM r M*

+ − == = +− + , (29)

*ff nf O*

*ff nf*

*rM r M*

*rM r M M*

+ − == = +− + (30)

(31)

*ff nf O*

*n*

 ρ

 ρ ρ

+ − = = +− + (27)

(24)

σ

(25)

 − ×

ζ

The resonance escape probability for fast neutrons (*p*) also is calculated as following [5]:

*p e*

3.9 *<sup>s</sup> eff I*

ξ=

 σ

According to enrichment of applied fuel and its mass can write [7]:

2

*U f fm fm fm fm*

*U ff fm fm fm fm*

If the transfer function of G(s) as *V* function is assumes then:

ρ

235

σ

Where:

and:

Where [8]:

and:

will be maximum value.

$$\dot{\mathbf{G}}(\rho) = \nabla \mathbf{G}^{\top}(\rho).\dot{\rho} \tag{32}$$

Firstly the transfer function is supposed. Secondly due to mostly reactor's cores are twin therefore once a signal with delay time ( *<sup>d</sup>* τ ) from first part to second part is transmitted then that signal with a same delay time from second part to first part will be transmitted.

#### **1.5 Neutron point kinetics (NPK)**

The treatment of the neutron transport as a diffusion process has only been validated. For example, in a Light Water Reactor (LWR) the mean free path of thermal neutrons is typically around 1 cm.

The fractional model has been derived for the NPK equations with n groups of delayed neutrons, is given by:

$$\pi^{K}\frac{d^{k+1}n}{dt^{k+1}} + \pi^{K}\left[\frac{1}{l} + \frac{(1-\beta)}{\Lambda}\right]\frac{d^{k}n}{dt^{k}} + \frac{dn}{dt} = \frac{\rho-\beta}{\Lambda}n + \sum\_{i=1}^{m} \lambda\_{i}\mathbb{C}\_{i} + \pi^{k}\sum\_{i=1}^{m} \left(\lambda\_{i}\frac{d^{k}\mathbb{C}\_{i}}{dt^{k}}\right), \ 0 < k \le 2 \tag{33}$$

Where:

*τ* is the relaxation time, *k* is the anomalous diffusion order (for sub-diffusion process: 0 < *k* < 1; while that for super-diffusion process: 1 < *k* < 2), *n* is the neutron density, Ci is the concentration of delayed neutron precursor, l is the prompt-neutron lifetime for finite media, K is the neutron generation time, β is the fraction of delayed neutrons, and ρ is the reactivity. When 0 *<sup>k</sup>* τ→ , the classic NPK equation is recovered.

The fractional model includes three additional terms relating to the classic equations which are contained of fractional derivatives (Gilberto and Espinosa, 2011):

$$\frac{d^{k+1}n}{dt^{k+1}}\text{ .}\tag{34}$$

$$\frac{d^k n}{dt^k} \tag{35}$$

and:

$$\frac{d^k \mathbb{C}\_i}{dt^k} \tag{36}$$

The physical meaning of above terms suggests that for sub-diffusion processes, the first term has an important contribution for rapid changes in the neutron density (for example in the turbine trip in a BWR nuclear power plant (NPP)), while the second term represents an important contribution when the changes in the neutron density is almost slow, for example during startup in a NPP that involves operational maneuvers due to movement of control rod mechanism. The importance of third term is when the reactor sets in shutdown state, it

The Theoretical Simulation of a Model by SIMULINK for

*n s G s*

*s*

Therefore the transfer function will be as shown:

δ

δρ

β

+ −

*n s G s*

Also the equation of zero power transfer function is as following:

[ ]

1 0

δ

δρ

δρ

0

λ

*n s*

2

β

β

*n s n s <sup>s</sup> s s s s s s*

and for variation of fission fragment density as per error reactivity variation, it can write:

( )

 + <sup>Λ</sup>

λ

β

Λ

*s*

+

λ

λ

δ

δρ *s*

[ ]

*s*

1 0

Transfer Function:

Where:

Where:

Surveying the Work and Dynamical Stability of Nuclear Reactors Cores 185

*adj*

+ − <sup>Λ</sup> <sup>Λ</sup> <sup>−</sup> <sup>+</sup> <sup>Λ</sup> = = <sup>=</sup>

*s*

*s n*

<sup>+</sup> <sup>Λ</sup> <sup>+</sup> Λ Λ

λ λ

β

+ − Λ <sup>−</sup> <sup>+</sup> Λ

β

β

*s*

*s*

<sup>0</sup> ( ) ( ) ( ) *<sup>e</sup>*

0

(43)

(44)

(45)

(47)

(46)

*s n*

λ

λ

β

β

β

+ − Λ <sup>−</sup> <sup>+</sup> Λ

β

> β

λ

λ

<sup>2</sup> ( )( )

*s n*

 β

β λ

λ λ

β

<sup>+</sup> <sup>Λ</sup> <sup>+</sup> Λ Λ = =

*s*

+ + Λ β

0

0

λ

 β

 λ

*s s ss*

λ

<sup>Λ</sup> = + + = ++ <sup>−</sup> Λ Λ

[ ]

*<sup>s</sup> s s*

( ) <sup>0</sup> ( ) ( ) ( ) *<sup>e</sup>*

<sup>0</sup> () () *e f*

0 0

*n n*

( ) ( ) ( ) ( ) ( ) ( ) () *<sup>e</sup>*

<sup>+</sup> <sup>Λ</sup> <sup>Λ</sup> <sup>+</sup> = == Λ + ++ +

Λ Λ

*s*

 β

λ

 ρ *s s* = −δρ

1 0

β *s*

λ

λ

*s*

0

0

[ ]

1 0

could also be important to understand the processes in the accelerator driven system (ADS), which is a subcritical system characterized by a low fraction of delayed neutrons and by a small Doppler reactivity coefficient and totally there are many interesting problems to consider under the view point of fractional differential equations (FDEs) (Gilberto and Espinosa, 2011).

The neutron point kinetics (NPK) equations are one of the most important reduced models of nuclear engineering, and they have been the subject of countless studies and applications to understand the neutron dynamics and its effects. These equations are shown below:

$$\frac{dn(t)}{dt} = \dot{n}(t) = \frac{\rho\_e - \beta}{\Lambda} n(t) + \sum\_{i=1}^{6} \lambda\_i \mathbb{C}\_i(t) \tag{37}$$

and:

$$\frac{d\mathbf{C}\_i(t)}{dt} = \dot{\mathbf{C}}\_i(t) = \frac{\beta}{\Lambda} n(t) - \lambda \mathbf{C}\_i(t) \tag{38}$$

Where:

*n t* ( ) is: the variations of neutrons density, *ρe* is: injected reactivity to system's transfer function, *β* is: delayed neutrons fraction, Λ is: neutron generation time, *λi* is: decay coefficient per density of each group of delayed neutrons, *Ci* is: density of each group of fission fragments which are as delayed neutrons generators.

If the partial variations per each system and control variables including: *n t* ( ) *, ρe*, *n t*( ) *, C t*( ) and*C t*( ) are considered the these can write as following:

$$
\delta\dot{n}(t) = \frac{\delta\phi - \beta}{\Lambda}(n\_0 + \delta n(t)) + \mathcal{A}(\mathbb{C}\_0 + \delta\mathbb{C}(t))\tag{39}
$$

and:

$$
\delta \dot{\mathcal{C}}(t) = \frac{\beta}{\Lambda} (n\_0 + \delta n(t)) - \lambda (\mathcal{C}\_0 + \delta \mathcal{C}(t)) \tag{40}
$$

If these both δ*n t* ( ) and ( ) δ*C t* equations to matrix format are written, then the matrix format of them is shown as below:

$$
\begin{bmatrix}
\delta\dot{n}(t) \\
\delta\dot{\mathbf{C}}(t)
\end{bmatrix} = \begin{bmatrix}
\Lambda & \\
\frac{\partial}{\Lambda} & -\lambda
\end{bmatrix} \times \begin{bmatrix}
\delta n(t) \\
\delta \mathbf{C}(t)
\end{bmatrix} + \begin{bmatrix}
n\_0 \\
\Lambda \\
0
\end{bmatrix} \delta\rho\tag{41}
$$

In the linear state can write last equation as the following:

$$
\delta\dot{m}(t) = -\frac{\beta}{\Lambda}\delta n(t) + \lambda\delta\mathcal{C}(t) + \frac{n\_0\delta\rho}{\Lambda} \tag{42}
$$

In case the point kinetic equations of reactor from both neutron density and fission fragments aspects are considered, it can be written:

$$\begin{aligned} \text{Transfer Function:} \qquad \quad G(s) = \frac{\delta n(s)}{\delta \rho\_{\epsilon}(s)} = \frac{\begin{bmatrix} 1 & 0 \end{bmatrix} adj\begin{bmatrix} s + \frac{\beta}{\Lambda} & -\lambda \\ -\frac{\beta}{\Lambda} & s + \lambda \end{bmatrix} \begin{bmatrix} \frac{\eta\_{0}}{\Lambda} \\ 0 \end{bmatrix}}{\begin{bmatrix} s + \frac{\beta}{\Lambda} & -\lambda \\ s + \frac{\beta}{\Lambda} & s + \lambda \end{bmatrix}} = \frac{\begin{bmatrix} \frac{\beta}{\Lambda} & -\frac{\beta}{\Lambda} \end{bmatrix} \begin{bmatrix} \frac{\eta\_{0}}{\Lambda} & -\frac{\beta}{\Lambda} \\ -\frac{\beta}{\Lambda} & s + \lambda \end{bmatrix}}{\begin{bmatrix} -\frac{\beta}{\Lambda} & s + \lambda \end{bmatrix}} \\\\ \begin{bmatrix} \frac{\tau}{s + \lambda} & \frac{\tau}{\Lambda} & \frac{\tau}{\Lambda} \end{bmatrix} \begin{bmatrix} \frac{\tau}{\Lambda} \end{bmatrix} \end{aligned}$$

$$
\begin{array}{c|c}
\begin{bmatrix} 1 & 0 \end{bmatrix} \begin{bmatrix} s+\lambda & \lambda \\ \hline \hline \end{bmatrix} & \begin{bmatrix} \lambda \\ s+\lambda \end{bmatrix} \begin{bmatrix} \frac{n\_0}{\Lambda} \\ 0 \end{bmatrix} \\\hline \begin{array}{c|c} \frac{\beta}{s+\beta} & -\lambda \\ \hline \hline \end{array} \\\hline \begin{array}{c|c} -\beta & s+\lambda \\ \hline \end{array} & \begin{array}{c|c} -\beta \end{array} \end{array} \tag{43}
$$

Where:

184 Nuclear Reactors

could also be important to understand the processes in the accelerator driven system (ADS), which is a subcritical system characterized by a low fraction of delayed neutrons and by a small Doppler reactivity coefficient and totally there are many interesting problems to consider under the view point of fractional differential equations (FDEs) (Gilberto and

The neutron point kinetics (NPK) equations are one of the most important reduced models of nuclear engineering, and they have been the subject of countless studies and applications to understand the neutron dynamics and its effects. These equations are shown below:

( ) () () () *<sup>e</sup> i i*

( ) () () () *<sup>i</sup> i i i dC t C t nt C t*

Λ

*n t* ( ) is: the variations of neutrons density, *ρe* is: injected reactivity to system's transfer function, *β* is: delayed neutrons fraction, Λ is: neutron generation time, *λi* is: decay coefficient per density of each group of delayed neutrons, *Ci* is: density of each group of

If the partial variations per each system and control variables including: *n t* ( ) *, ρe*, *n t*( ) *,* 

0 0 *nt n nt C Ct* ( ) ( ( )) ( ( ))

0 0 *Ct n nt C Ct* ( ) ( ( )) ( ( ))

= + −+

 λδ

( ) *n t C t* δ

+

δ

 δ

λ

<sup>×</sup> ( )

λ

<sup>0</sup> () () () *<sup>n</sup> nt nt Ct*

In case the point kinetic equations of reactor from both neutron density and fission

 λδ=− + +

 λδ

<sup>Λ</sup> (39)

(40)

0

(41)

0 *n* Λ δρ

δρ

Λ Λ (42)

equations to matrix format are written, then the matrix format

δ

<sup>−</sup> = + ++

β== −

*dn t nt nt C t*

ρ β

<sup>−</sup> == +

*dt*

fission fragments which are as delayed neutrons generators.

*C t*( ) and*C t*( ) are considered the these can write as following:

δ

δ

( ) ( ) *n t C t*

δ

In the linear state can write last equation as the following:

δ

δ

fragments aspects are considered, it can be written:

δρ β

> β

Λ

β

−

β

 δ

β

 <sup>Λ</sup> <sup>=</sup> <sup>−</sup> <sup>Λ</sup>

*dt*

6

1

λ

<sup>Λ</sup> (37)

(38)

*i*

=

λ

Espinosa, 2011).

and:

Where:

and:

If these both

δ

of them is shown as below:

*n t* ( ) and ( ) δ*C t*

$$\begin{vmatrix} s + \frac{\beta}{\Lambda} & -\lambda \\ \frac{\Lambda}{\Lambda} & s + \lambda \\ \frac{-\beta}{\Lambda} & s + \lambda \end{vmatrix} = s^2 + s(\lambda + \frac{\beta}{\Lambda}) = s(s + \lambda + \frac{\beta}{\Lambda}) \tag{44}$$

Therefore the transfer function will be as shown:

$$G(s) = \frac{\delta n(s)}{\delta \rho\_c(s)} = \frac{\begin{bmatrix} 1 & 0 \end{bmatrix} \begin{bmatrix} s + \lambda & \lambda \\ \hline \Delta & s + \frac{\beta}{\Lambda} \end{bmatrix} \begin{bmatrix} n\_0 \\ \Lambda \\ 0 \end{bmatrix}}{s(s + \lambda + \frac{\beta}{\Lambda})}\tag{45}$$

Where:

$$
\delta \mathfrak{d} \mathfrak{p}\_{\mathfrak{e}}(\mathbf{s}) = \mathfrak{p}\_{\mathfrak{o}} - \delta \mathfrak{p}\_{f}(\mathbf{s}) \tag{46}
$$

Also the equation of zero power transfer function is as following:

$$
\begin{bmatrix} 1 & 0 \end{bmatrix} \begin{bmatrix} \frac{n\_0(s+\lambda)}{\Lambda} \\ \frac{\beta n\_0}{\Lambda^2} \\ \frac{\beta n\_0}{\Lambda} \end{bmatrix} = \frac{n\_0}{\frac{\Lambda}{\Lambda}(s+\lambda)} = \frac{n\_0}{s(s+\lambda)} = \frac{n\_0(s+\lambda)}{s(s\Lambda+\beta)}\tag{47}
$$

and for variation of fission fragment density as per error reactivity variation, it can write:

The Theoretical Simulation of a Model by SIMULINK for

In large twin reactor's cores, it can be written:

*q1* , *q2* and *ρ*i are respectively as following:

ρ

*<sup>i</sup>* equals:

According to Newton's low of cooling can write:

variables are considered respectively as below:

ρ

So can write:

It is also supposed:

and:

and:

So ( ) *<sup>i</sup> q t* divided by

Surveying the Work and Dynamical Stability of Nuclear Reactors Cores 187

() () *<sup>T</sup> Gs c sb* = φ

ρ β 6

1 () () () () *<sup>e</sup> ii i i nt nt C t q t*

= <sup>−</sup> = ++

1 () ( ) *ij i j ij j i i qt nt T* α

= ≠

α α

12 1 2 12 1 *qt nt T* () ( ) α

21 2 1 21 2 *qt nt T* () ( ) α

> ( ) ( ) *j ij*

> > 1

= ≠

*nt T n t*

( ) ( ) ( ) ( ) ( )

*j ij j i i i i j ij i i j i ij*

<sup>−</sup> <sup>Λ</sup> = = <sup>−</sup> <sup>Λ</sup>

<sup>1</sup> () () () *p T t n t aT t*

If the variations of neutron density than initial neutron numbers, fission fragment density than initial fission fragment density and temperature than initial temperature as system

*nt T q t n t nt T n t*

1

= ≠

*ij*

α

*j i i*

α

*i ij j i i*

 α

ρ

1

= ≠

1

= ≠

λ

<sup>Λ</sup> (54)

= − <sup>Λ</sup> (55)

12 21 = (56)

= − <sup>Λ</sup> , (57)

= − <sup>Λ</sup> (58)

<sup>−</sup> <sup>=</sup> (59)

*mc* = − (61)

(60)

(53)

$$
\frac{\delta \mathcal{C}(s)}{\delta \rho\_c(s)} = \frac{\begin{bmatrix} n\_0(s + \lambda) \\ \Lambda \\ \frac{\beta n\_0}{\Lambda^2} \end{bmatrix}}{s(s + \lambda + \frac{\beta}{\Lambda})} = \frac{\beta n\_0/\Lambda^2}{s(s + \frac{\beta}{\Lambda})} \tag{48}
$$

In case the point kinetic equations of reactor from neutron density and fission fragments aspects and also temperature aspect (according to Newton's low of cooling) are considered, it can be expressed:

$$G(s) = \frac{\delta n(s)}{\delta \rho\_{\varepsilon}(s)} = \frac{\begin{bmatrix} \lambda & 0 & 0 \end{bmatrix} adj}{\begin{bmatrix} -\beta & \lambda & 0 \\ \Lambda & s + \lambda & 0 \end{bmatrix} \begin{bmatrix} n\_0 \\ \Lambda \\ 0 \\ 0 \\ 0 \end{bmatrix}} \begin{bmatrix} n\_0 \\ \Lambda \\ 0 \\ 0 \\ 0 \end{bmatrix} \tag{49}$$
 
$$G(s) = \frac{\delta n(s)}{\delta \rho\_{\varepsilon}(s)} = \frac{\begin{bmatrix} \beta & -\beta & \lambda & -\beta & -\lambda & -\beta \end{bmatrix}}{\begin{bmatrix} \lambda & -\beta & \lambda & 0 \\ -\beta & s + \lambda & 0 \\ -\beta & -s & -\beta \end{bmatrix}} \tag{40}$$

Where:

$$\begin{vmatrix} s + \frac{\beta}{\Lambda} & -\lambda & \frac{\alpha n\_0}{\Lambda} \\ -\frac{\beta}{\Lambda} & s + \lambda & 0 \\ \frac{\Lambda}{\Lambda} & 0 & s + a \end{vmatrix} = s^3 + s^2(a + \lambda + \frac{\beta}{\Lambda}) + s(a\lambda + \frac{\beta a}{\Lambda} + \frac{\alpha n\_0 k}{\Lambda}) + \frac{\alpha n\_0 k \lambda}{\Lambda} \tag{50}$$

In this state the equation of zero power transfer function is as following:

$$G(s) = \frac{\delta n(s)}{\delta \rho\_\epsilon(s)} = \frac{(n\_0/\Lambda)(s+\lambda)(s+a)}{s^3 + s^2(a+\lambda+\frac{\beta}{\Lambda}) + s(a\lambda + \frac{\beta a}{\Lambda} + \frac{\alpha n\_0 k}{\Lambda}) + \frac{\alpha n\_0 k \lambda}{\Lambda}}\tag{51}$$

In this stage according to the main transfer function that has mentioned in the last stage, *G*( ) ρis defined:

$$G(\rho) = \frac{\delta n(\rho)}{\delta \rho\_\epsilon(\rho)} = \frac{(n\_0/\Lambda)(\rho + \lambda)(\rho + a)}{\rho^3 + \rho^2(a + \lambda + \frac{\beta}{\Lambda}) + \rho(a\lambda + \frac{\beta a}{\Lambda} + \frac{\alpha n\_0 k}{\Lambda}) + \frac{\alpha n\_0 k \bar{\lambda}}{\Lambda}} \tag{52}$$

Therefore through *G*( ) ρ , the parts of ( ) *G* ρ will be defined and the stability condition through *V* function is applied.

So can write:

186 Nuclear Reactors

0

*C s n <sup>s</sup> s s s s*

λ

*n s*

0

β

*n*

β

( ) ( ) () *<sup>e</sup>*

In case the point kinetic equations of reactor from neutron density and fission fragments aspects and also temperature aspect (according to Newton's low of cooling) are considered,

( )

 + <sup>Λ</sup> <sup>Λ</sup> <sup>Λ</sup> = = ++ + Λ Λ

λ

2 2 0

100 0 0

β

*adj s*

β

*s n s*

β

+ −

<sup>−</sup> <sup>+</sup> Λ

− + = =

β

0

− +

3 2 0 0 0 ( )( )

3 2 0 0

 λ

3 2 0 0

ρ λ ρ

+ ++ + + + +

 ρλ

λ

 β

Λ ΛΛ Λ

 β

Λ ΛΛ Λ

α

*s a nk nk*

+ ++ + + + +

λλ

β

0

β

In this stage according to the main transfer function that has mentioned in the last stage,

0

*a a*

β

( ) ( )( )( ) ( ) ( ) ( )( ) *<sup>e</sup> n s n s sa G s*

*s s a sa*

Λ+ + = =

λ

( ) ( )( )( ) ( ) ( ) ( )( ) *<sup>e</sup> <sup>n</sup> n a <sup>G</sup>*

Λ+ + = =

 λ

ρ

<sup>−</sup> + = + ++ + + + + Λ Λ ΛΛ Λ

*k sa*

λ

Λ Λ

*s*

β

 β

0

*n*

α λ

*s n*

+ − Λ Λ <sup>Λ</sup> <sup>−</sup>

λ

<sup>+</sup> <sup>Λ</sup>

*k sa*

α λ

0

0

 β

<sup>0</sup> <sup>0</sup>

0

*a nk nk*

 αλ

*a nk nk*

will be defined and the stability condition

 αλ

α

 αλ

α

(48)

(49)

(50)

(51)

(52)

[ ]

[ ]

( )

δ

δρ

( ) ( ) ( )

δ

δρ

*n s G s*

*e*

0

*s s s a sa*

In this state the equation of zero power transfer function is as following:

*n*

α λ

0

− +

*k sa*

δ

δρ

δ ρ

δρ ρ

ρ

ρ ρ

, the parts of ( ) *G*

ρ λ

Λ Λ

*s*

β

+ −

β

it can be expressed:

Where:

*G*( ) ρ

is defined:

Therefore through *G*( )

through *V* function is applied.

0 1

$$\mathbf{G}(\mathbf{s}) = \mathbf{c}^T \boldsymbol{\phi}(\mathbf{s}) \mathbf{b} \tag{53}$$

In large twin reactor's cores, it can be written:

$$
\dot{m}(t) = \frac{\rho\_c - \beta}{\Lambda} n(t) + \sum\_{i=1}^{6} \lambda\_i \mathbb{C}\_i(t) + \eta\_i(t) \tag{54}
$$

It is also supposed:

$$q\_i(t) = \sum\_{j=1 \neq i} \frac{\alpha\_{ij}}{\Lambda\_i} n\_j(t - T\_{ij}) \tag{55}$$

and:

$$
\alpha\_{12} = \alpha\_{21} \tag{56}
$$

*q1* , *q2* and *ρ*i are respectively as following:

$$q\_1(t) = \frac{\alpha\_{12}}{\Lambda\_1} n\_2(t - T\_{12}) \tag{57}$$

$$
\eta\_2(t) = \frac{\alpha\_{21}}{\Lambda\_2} n\_1(t - T\_{21}) \tag{58}
$$

and:

$$\rho\_i = \sum\_{j=1 \neq i} \alpha\_{ij} \frac{n\_j(t - T\_{ij})}{n\_i(t)} \tag{59}$$

So ( ) *<sup>i</sup> q t* divided by ρ*<sup>i</sup>* equals:

$$\frac{q\_i(t)}{\rho\_i} = \frac{\sum\_{j=1 \neq i} \frac{\alpha\_{ij}}{\Lambda\_i} n\_j(t - T\_{ij})}{\sum\_{j=1 \neq i} \alpha\_{ij} \frac{n\_j(t - T\_{ij})}{n\_j(t)}} = \sum\_{j=1 \neq i} \frac{n\_i(t)}{\Lambda\_i} \tag{60}$$

According to Newton's low of cooling can write:

$$
\dot{T}(t) = \frac{1}{mc\_p} n(t) - nT(t) \tag{61}
$$

If the variations of neutron density than initial neutron numbers, fission fragment density than initial fission fragment density and temperature than initial temperature as system variables are considered respectively as below:

The Theoretical Simulation of a Model by SIMULINK for

the perturbation theory (Gupta and Trasi, 1986).

equations has been applied (Merk and Cacuci, 2005).

(Pandey, 1996; Tsuji et al., 1993; Konno et al., 1994).

This function has features as following:

2. The partial differential of *V x*( ) is continued.

4. The *V x*( ) function can be written as following:

1. *V x*( ) has definite positive value.

3. Where the *ρ* is momentum, *V*( )

(Aboanberand, 2002).

used (Nahla, 2009).

*V(x)* Function.

It can also write:

reactors (Ward and Lee, 1987).

for analyzing a point reactor model (Devooght and Smets, 1967).

the effect of time-delay in the automated control system (Konno et al, 1992).

group delayed neutron have also been applied (Munoz and Verdu, 1991).

system the stability condition is when the total energy of a system decreases.

ρ

Surveying the Work and Dynamical Stability of Nuclear Reactors Cores 189

The topological and Lyapunov methods were compared with Aizermann-Rosen methods

The Padé approximations has been used to obtain solutions for point kinetic equations

Perturbation theory has also been widely used in studying reactor dynamics. There has been obtained specific types of steady solutions to study power oscillations in a reactor results from a Hopf bifurcation (Pandey, 1996; Munoz and Verdu, 1991; Tsuji et al., 1993; Konno et al., 1994). The KBM theory has been used for the nonlinear analysis of a reactor model with

A singular perturbation has been used to study relaxation oscillations in typical nuclear

The point kinetic equations in the presence of delayed neutrons with one temperature reactivity coefficient for a step input of reactivity have analytically been solved by applying

The regular perturbations to obtain an analytical solution for general reactivity have been

The multiple time-scales expansion to obtain analytical solutions of the neutron kinetic

The variation methods in conjunction with the Hopf bifurcation theory for a BWR with one

A combination of the center-manifold reduction (CMR) and the method of normal forms have already been applied widely for nonlinear analysis of nuclear reactor dynamics

In the Liapunov model, the dynamically equations may be either differential or non differential equations. But the method of problem evaluation is not based on solving the equations. This method is based on energy in classic mechanical. In one of mechanical

Liapunov applied this property and based the stability function. This function is entitled:

1 2

*V x dx V x dx V x dx V x*

∂∂ ∂ = + ++ ∂∂ ∂

() () ( ) () . . ... . *<sup>n</sup>*

1 2

is negative quasi relatively.

*x dt x dt x dt*

*n*

() () *<sup>T</sup> Vx V xx* = ∇ (70)

(69)

$$\mathbf{x}(t) = \frac{n(t) - n\_0}{n\_0},\tag{62}$$

$$y(t) = \frac{C(t) - C\_0}{C\_0} \tag{63}$$

and:

$$z(t) = \frac{T(t) - T\_0}{T\_0} \tag{64}$$

Then the differentials of these system variables are respectively as following:

$$\dot{\mathbf{x}}\_{i}(t) = \frac{\mathcal{S}\_{i}(\mathbf{x},t)}{\Lambda} [1 + \mathbf{x}\_{i}(t)] - \frac{\mathcal{B}}{\Lambda} \mathbf{x}\_{i}(t) + \frac{\mathcal{B}}{\Lambda} \mathbf{y}\_{i}(t) - \sum\_{j=1 \neq i}^{N} \frac{n\_{j\_{0}}}{n\_{i\_{0}}} \cdot \frac{\alpha\_{\overline{\mathbf{y}}}}{\Lambda} \mathbf{x}\_{i}(t) + \sum\_{j=1 \neq i}^{N} \frac{n\_{j\_{0}}}{n\_{i\_{0}}} \cdot \frac{\alpha\_{\overline{\mathbf{y}}}}{\Lambda} \mathbf{x}\_{j}(t - T\_{\overline{\mathbf{y}}}) \,, \tag{65}$$

$$
\dot{y}\_i(t) = \lambda x\_i(t) - \lambda y\_i(t) \tag{66}
$$

and:

$$
\dot{z}\_i(t) = a x\_i(t) - a z\_i(t) \tag{67}
$$

If these system variables differentials are to matrix format are written, then can write:

$$
\begin{bmatrix}
\dot{x}\_i(t) \\
\dot{y}\_i(t) \\
\dot{z}\_i(t)
\end{bmatrix} = \begin{bmatrix}
\frac{\delta\_i - \beta}{\Lambda} - \sum\_{j=1 \neq i}^{N} \frac{n\_{j\bullet}}{n\_{i\_0}} \cdot \frac{\beta\_0}{\Lambda} & \frac{\beta}{\Lambda} \\
& \lambda & 0 \\
& & a & 0 \\
& & & 0 & -a
\end{bmatrix} \times \begin{bmatrix}
x\_i(t) \\
y\_i(t) \\
\dot{y}\_i(t) \\
z\_i(t)
\end{bmatrix} + \begin{bmatrix}
\frac{\delta\_i}{\Lambda} \\
0 \\
0 \\
\end{bmatrix} \tag{68}
$$

#### **1.6 Stability of reactor**

There are some methods for determination of nuclear reactor stability. Among methods which are applied for this aim are as following:

Liapunov method- Lagrange method- Popov method- Pontryagin method.

There is Lyapunov's method to determine the stability of nonlinear reactor dynamics by constructing certain positive definite functions of the reactor variables and parameters (Pankaj and Vivek, 2011).

There has been calculated the Lyapunov exponents from a time series of the excess neutron population of a boiling water reactor (BWR) and used it to conclude about the stability of the steady state operation of that particular BWR (Munoz et al., 1992).

There has also discussed the application of topological methods in reactor kinetics study (Smets and Giftopoulos, 1959).

( ) ( ) *nt n x t n*

( ) ( ) *Ct C y t <sup>C</sup>*

( ) ( ) *Tt T z t T*

( ,) ( ) [1 ( )] ( ) ( ) . ( ) . ( ) *N N j ij j ij <sup>i</sup> i i ii i j ij*

= +− + − + − Λ ΛΛ Λ Λ , (65)

> () () () *iii y t xt yt* = − λ

 λ

*<sup>N</sup> j ij <sup>i</sup> <sup>i</sup>*

 λ

β

*x t n n x t xt xt y t x t xt T*

If these system variables differentials are to matrix format are written, then can write:

0 1 <sup>0</sup> . ( ) ( )

*j i i i i i i i i*

λ

Liapunov method- Lagrange method- Popov method- Pontryagin method.

the steady state operation of that particular BWR (Munoz et al., 1992).

= ≠

*n x t n x t y t y t z t a a z t*

( ) 0 () 0 ( ) 0 () 0

There are some methods for determination of nuclear reactor stability. Among methods

There is Lyapunov's method to determine the stability of nonlinear reactor dynamics by constructing certain positive definite functions of the reactor variables and parameters

There has been calculated the Lyapunov exponents from a time series of the excess neutron population of a boiling water reactor (BWR) and used it to conclude about the stability of

There has also discussed the application of topological methods in reactor kinetics study

<sup>−</sup> <sup>−</sup> Λ ΛΛ <sup>Λ</sup> <sup>=</sup> × + <sup>−</sup>

α

Then the differentials of these system variables are respectively as following:

 β

β

δ β

 

which are applied for this aim are as following:

**1.6 Stability of reactor** 

(Pankaj and Vivek, 2011).

(Smets and Giftopoulos, 1959).

and:

and:

δ

0 0

> 0 0

0 0

<sup>−</sup> <sup>=</sup> , (62)

<sup>−</sup> <sup>=</sup> (63)

<sup>−</sup> <sup>=</sup> (64)

 α

(66)

δ

(68)

0 0 1 1 0 0

() () () *i ii z t ax t az t* = − (67)

*n n*

*j i i i j i*

= ≠ = ≠

α

The topological and Lyapunov methods were compared with Aizermann-Rosen methods for analyzing a point reactor model (Devooght and Smets, 1967).

The Padé approximations has been used to obtain solutions for point kinetic equations (Aboanberand, 2002).

Perturbation theory has also been widely used in studying reactor dynamics. There has been obtained specific types of steady solutions to study power oscillations in a reactor results from a Hopf bifurcation (Pandey, 1996; Munoz and Verdu, 1991; Tsuji et al., 1993; Konno et al., 1994). The KBM theory has been used for the nonlinear analysis of a reactor model with the effect of time-delay in the automated control system (Konno et al, 1992).

A singular perturbation has been used to study relaxation oscillations in typical nuclear reactors (Ward and Lee, 1987).

The point kinetic equations in the presence of delayed neutrons with one temperature reactivity coefficient for a step input of reactivity have analytically been solved by applying the perturbation theory (Gupta and Trasi, 1986).

The regular perturbations to obtain an analytical solution for general reactivity have been used (Nahla, 2009).

The multiple time-scales expansion to obtain analytical solutions of the neutron kinetic equations has been applied (Merk and Cacuci, 2005).

The variation methods in conjunction with the Hopf bifurcation theory for a BWR with one group delayed neutron have also been applied (Munoz and Verdu, 1991).

A combination of the center-manifold reduction (CMR) and the method of normal forms have already been applied widely for nonlinear analysis of nuclear reactor dynamics (Pandey, 1996; Tsuji et al., 1993; Konno et al., 1994).

In the Liapunov model, the dynamically equations may be either differential or non differential equations. But the method of problem evaluation is not based on solving the equations. This method is based on energy in classic mechanical. In one of mechanical system the stability condition is when the total energy of a system decreases.

Liapunov applied this property and based the stability function. This function is entitled: *V(x)* Function.

This function has features as following:


$$\dot{V}(\mathbf{x}) = \frac{\partial V(\mathbf{x})}{\partial \mathbf{x}\_1} \frac{d\mathbf{x}\_1}{dt} + \frac{\partial V(\mathbf{x})}{\partial \mathbf{x}\_2} \frac{d\mathbf{x}\_2}{dt} + \dots + \frac{\partial V(\mathbf{x})}{\partial \mathbf{x}\_n} \frac{d\mathbf{x}\_n}{dt} \tag{69}$$

It can also write:

$$
\dot{V}(\mathbf{x}) = \nabla V^T(\mathbf{x}) \dot{\mathbf{x}} \tag{70}
$$

The Theoretical Simulation of a Model by SIMULINK for

the secondary value of Keff in the recent position of control rod.

The simulated model is considered according to Fig.3:

Fig. 3. The block diagram of simulation by SIMULINK

**2.2 The simulation by SIMULINK and related block diagram** 

amount of the feedback which has been sampled as follows: *x(t)*.

Where *dx*

So:

Thus:

Supposition: *x*(0) 0 = ,

simulation (Tewari, 2002).

Surveying the Work and Dynamical Stability of Nuclear Reactors Cores 191

<sup>0</sup> <sup>0</sup> () [ () ( ) ] *<sup>t</sup>*

Where *x* is: absolutely descending, *0x* is: the initial value of *x* and *tD* is: the innate delay time. The SIMULINK of MATLAB is an appropriate software to analyze the performance of this

The work processes have been simulated by the SIMULINK of MATLAB software and all responses such as oscillation and transient responses have been analyzed by it as well. The main function *(F)* is *Fcn*. It includes two other functions that are: *u*[1] and *u*[2] which are defined for SIMULINK. *u*[1] is one of two input functions of Mux that has been shown by the input block that is: *H*. This block is presenting amount of the Keff. [2] *u* equals with

sgn( ) *sp dx v FK*

*dt* is: the velocity of control rod from the movement aspect to up and down, *Ksp* :

*dt* = − (76)

*<sup>D</sup> sp x t v kx t H t t K K dt* <sup>=</sup> − +− *sgn* (77)

<sup>0</sup> *x t x vt t t* () ( ) =± − *sgn <sup>D</sup>* (78)

<sup>0</sup> *Fcn F ku u K* == × + [1] [2] (79)

In fact the Liapunov function is a function that considers either state variables or variables which cause to imbalance state. Due to the potential function is able to do a process, so the Liapunov function is as *V* function. One can write:

$$V = f(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_n) \tag{71}$$

If: *V* = 0 then the system will be steady state.

This method revolved about the determination of a *V* function, which satisfies certain requirements of the stability theorem. In an initial reactor model a point reactor with constant power removal and without delayed neutrons is considered.

Due to *ρ* is a variable which causes unbalancing the system, therefore it can be considered as an *x* variable. Then it can be written:

$$\dot{V}(\rho) = \frac{\partial V(\rho)}{\partial \rho\_1} \frac{d\rho\_1}{dt} + \frac{\partial V(\rho)}{\partial \rho\_2} \frac{d\rho\_2}{dt} + \dots + \frac{\partial V(\rho)}{\partial \rho\_n} \frac{d\rho\_n}{dt} \tag{72}$$

and:

$$\dot{V}(\rho) = \nabla V^{T}(\rho).\dot{\rho} \tag{73}$$

#### **2. Methodology**

#### **2.1 Analyzing the theory by mathematical model**

In this work the value of Keff as a comparable value is supposed and attributed to input parameter block (*H*) and then this value with the received feedback value is compared.

The unit of the control rod velocity (*v*) can be mm/s, the rate is steady, and the control rod movement is only to up and down directions, so: x(0) =0 (Shirazi et al., 2010).

Since the *sgn(x)* function is nonlinear; so conversion function can not be calculated; thus in this stage arguing the frequency response is not meaningful. Therefore the steady state must be considered for this nonlinear function; though it is rather complicated (Marie and Mokhtari, 2000).

To analyze the controlling system theory these are assumed:

If: Input=*H*; Output=*x(t)*; in the top of control rod: *x=*0; in the bottom of the control rod: *x=xmax*;

$$F = kH\mathbf{x}(t) + K\_0 \tag{74}$$

Where:

*F* : Function, *k* : constant coefficient, *H* : input parameter, *x t*( ) : the control rod position, *K*<sup>0</sup> : initial value of Keff.

$$
\Delta \mathfrak{x} = \upsilon \operatorname{sgn}(F - K\_{sp}) \Delta t \tag{75}
$$

Where Δ*x* is: the amount of control rod movement, Δ*t* is: time**.**

$$\frac{d\chi}{dt} = \upsilon \operatorname{sgn}(F - K\_{sp}) \tag{76}$$

Where *dx dt* is: the velocity of control rod from the movement aspect to up and down, *Ksp* : the secondary value of Keff in the recent position of control rod.

$$\mathbf{x}(t) = \int\_{0}^{t} \mathbf{v} \text{sgn}[k\mathbf{x}(t)H(t - t\_D) + K\_0 - K\_{sp}]dt\tag{77}$$

Supposition: *x*(0) 0 = ,

So:

190 Nuclear Reactors

In fact the Liapunov function is a function that considers either state variables or variables which cause to imbalance state. Due to the potential function is able to do a process, so the

This method revolved about the determination of a *V* function, which satisfies certain requirements of the stability theorem. In an initial reactor model a point reactor with

Due to *ρ* is a variable which causes unbalancing the system, therefore it can be considered as

1 2

*Vd Vd V <sup>d</sup> <sup>V</sup>*

∂∂ ∂ = + ++ ∂∂ ∂

() () ( ) ( ) *<sup>n</sup>*

 ρρ

In this work the value of Keff as a comparable value is supposed and attributed to input parameter block (*H*) and then this value with the received feedback value is compared.

The unit of the control rod velocity (*v*) can be mm/s, the rate is steady, and the control rod

Since the *sgn(x)* function is nonlinear; so conversion function can not be calculated; thus in this stage arguing the frequency response is not meaningful. Therefore the steady state must be considered for this nonlinear function; though it is rather complicated (Marie and

If: Input=*H*; Output=*x(t)*; in the top of control rod: *x=*0; in the bottom of the control rod:

*F* : Function, *k* : constant coefficient, *H* : input parameter, *x t*( ) : the control rod position,

*dt dt dt*

 ρρ

1 2

ρ

movement is only to up and down directions, so: x(0) =0 (Shirazi et al., 2010).

ρρ

constant power removal and without delayed neutrons is considered.

ρρ

1 2 ( , ,..., ) *V fx x x* = *<sup>n</sup>* (71)

*n*

= ∇ (73)

<sup>0</sup> *F kHx t K* = + ( ) (74)

sgn( ). *sp* Δ= − Δ *xv FK t* (75)

 ρ

(72)

 ρ ρ

Liapunov function is as *V* function. One can write:

ρ

( ) ( ). *<sup>T</sup> V V*

**2.1 Analyzing the theory by mathematical model** 

To analyze the controlling system theory these are assumed:

Where Δ*x* is: the amount of control rod movement, Δ*t* is: time**.**

If: *V* = 0 then the system will be steady state.

an *x* variable. Then it can be written:

and:

**2. Methodology** 

Mokhtari, 2000).

*K*<sup>0</sup> : initial value of Keff.

*x=xmax*;

Where:

$$\mathbf{x}(t) = \mathbf{x}\_0 \pm \upsilon t \mathbf{s} \mathbf{g}(t - t\_D) \tag{78}$$

Where *x* is: absolutely descending, *0x* is: the initial value of *x* and *tD* is: the innate delay time.

The SIMULINK of MATLAB is an appropriate software to analyze the performance of this simulation (Tewari, 2002).

The simulated model is considered according to Fig.3:

Fig. 3. The block diagram of simulation by SIMULINK

#### **2.2 The simulation by SIMULINK and related block diagram**

The work processes have been simulated by the SIMULINK of MATLAB software and all responses such as oscillation and transient responses have been analyzed by it as well. The main function *(F)* is *Fcn*. It includes two other functions that are: *u*[1] and *u*[2] which are defined for SIMULINK. *u*[1] is one of two input functions of Mux that has been shown by the input block that is: *H*. This block is presenting amount of the Keff. [2] *u* equals with amount of the feedback which has been sampled as follows: *x(t)*.

Thus:

$$Fcn = F = ku[1] \times u[2] + \mathbb{K}\_0\tag{79}$$

The Theoretical Simulation of a Model by SIMULINK for

Surveying the Work and Dynamical Stability of Nuclear Reactors Cores 193

Fig. 5. The low oscillation for: Set Point: 100, *v*: 1, H: 100, delay time: 15 ms

Fig. 6. The medium oscillation for: Set Point: 100, v: 5, H: 100, delay time: 15 ms

Because of the control rod movement is steady, in order to calculate the total amount of the discrete movements of control rod, the Discrete-time Integrator block has been used. The *Fcn* produced function has been transferred to Zero-Order Hold block which plays logic converter role. In addition the Transport Delay block is related to the inherent delay time that is: *tD*. The parameters which must be adjusted are: Set Point that is: the default amount of Keff as reference Keff and the meaning of the Set Point=100 is: Keff=1, the velocity of control rod (*v*), recent Keff (block H) and the stop time that is: the innate delay time or *tD*. The graphs can be observed by the oscilloscope.
