The Theoretical Simulation of a Model by SIMULINK for Surveying the Work and Dynamical Stability of Nuclear Reactors Cores

This book presents a comprehensive review of studies in nuclear reactors technology from authors across the globe. Topics discussed in this compilation include: thermal hydraulic investigation of TRIGA type research reactor, materials testing reactor and high temperature gas-cooled reactor; the use of radiogenic lead recovered from ores as a coolant for fast reactors; decay heat in reactors and spent-fuel pools; present status of two-phase flow studies in reactor components; thermal aspects of conventional and alternative fuels in supercritical water?cooled reactor; two-phase flow coolant behavior in boiling water reactors under earthquake condition; simulation of nuclear reactors core; fuel life control in light-water reactors; methods for monitoring and controlling power in nuclear reactors; structural materials modeling for the next generation of nuclear reactors; application of the results of finite group theory in reactor physics; and the usability of vermiculite as a shield for nuclear reactor.


Introduction
According to complexity of nuclear reactor technology, applying a highly developed simulation is necessary for controlling the nuclear reactor control rods, so in this proposal the processes of a controlling model for nuclear reactors have been developed and simulated by the SIMULINK tool kit of MATLAB software and all responses, including oscillation and transient responses, have been analyzed.
In this work an arbitrary value of K eff as a comparable value is purposed and attributed to input block (H) of diagram and then this value with the received feedback value from block diagram is compared. Since the stability of the cited simulation depends on either velocity or delay time values, therefore according to this simulation the best response and operation which a reactor can have from stability aspect, have been derived. Meantime by viewing the results, the best ranges of velocity and delay time of control rod movement (in unit per second and millisecond respectively) for stability a nuclear reactor has been deduced.
Though the highlights of this proposal are respectively the following: • Defining a mathematical model for control rod movement • Simulation of a mathematical model by SIMULINK of MATLAB • Determination of the best ranges for both velocity and delay time of control rod movement (in unit per second and millisecond respectively) based on the obtained results for stability an LWR nuclear reactor In view of the great advancing the nuclear reactors technology, the phenomenal and significant changes in evolution of made nuclear reactors is observed. Since the make of the first nuclear reactor on 1948 until modern reactors, too changes are obvious. The major of these changes to: the kind of reactor design, the percent of fuel enrichment, the kind of coolant and neutron moderator, more safety and the dimensions of core are referred. The power control system is a key control system for a nuclear reactor, which directly affects the safe operation of a nuclear reactor. Much attention has been spent to the power control system performance of nuclear reactor in engineering (Zhao et al., 2003).
High reliability is one of the main objectives of the design and operation of control systems in nuclear power plants (Basu and Zemdegs, 1978;Stark, 1976).
Prototyping a control-rod driving mechanism (CRDM), which is a crucial safety system in the Taiwan Research Reactor (TRR-II) has been implemented, by iterative parallel procedures. Hence to ensure the mechanical integrity and substance of the prototype, a series of performance testing and design improvement have been interactively executed. Functional testing results show that the overall performance of the CRDM meets the specification requirements (Chyou and Cheng, 2004).
Also the SCK.CEN/ININ joint project, which deals with the design and application of modern/expert control and real-time simulation techniques for the secure operation of a TRIGA Mark III research nuclear reactor, has been undertaken (Dong et al., 2009).
This project has been proposed as the first of its kind under a general collaboration agreement between the Belgian Nuclear Research Centre (SCK·CEN) and the National Nuclear Research Institute (ININ) of Mexico (Benítez et al., 2005). In addition to the fuzzy proportional-integral-derivative (fuzzy-PID) control strategy has been applied recently as a nuclear reactor power control system. In the fuzzy-PID control strategy, the fuzzy logic controller (FLC) is exploited to extend the finite sets of PID gains to the possible combinations of PID gains in stable region and the genetic algorithm (Cheng et al., 2009;Park and Cho, 1992).
Until now, manual controlling systems have been used for controlling and tuning the control rods in the core of Gen II and some Gen III reactors (Tachibana et al., 2004).
But by application of this simulation that is the subject of this proposal the best response for operating and the best velocity and delay time of control rod movement in which can be caused to stability and critical state of a nuclear reactor, have been derived.
The safe situation is state in which the reactor stabilizes in the critical situation, meaning that the period is infinite and the Keff is 1 (Lamarsh, 1975).

Accidents
In which two states the positive reactivity overcomes the temperature coefficient of reactivity ( T ): Increasing the power and temperature of reactor core might decrease concentration of boric acid. Accordingly this event might cause to inject positive reactivity.
In addition for the reactors which apply fuels including Pu, because of having a resonance for Pu in thermal neutrons range so through increasing core's temperature the related resonance is broadened and absorbs more neutrons and because of Pu is fissionable therefore fissionable absorption occurs and is caused excess reactivity.
Also either accident or unfavorable issues as a feedback can be considered. Accidents of a nuclear reactor are totally classified based on following: 1. Over power accident.

www.intechopen.com
Each mentioned issues are divided to other sub issues. Over power accident is due factors such as: 1. Control rod withdrawal including uncontrolled rod withdrawal at sub critical power and uncontrolled rod withdrawal at power that will cause power excursion. 2. Control rod ejection. 3. Spent fuel handling. 4. Stem line break. 5. External events such as earthquake, enemy attack and etc.
In each mentioned issues the positive reactivity to nuclear reactor core can be injected. But there is another important accident that is: under cooling accident. Under cooling accident is classified to three sub accident among: LOCA (loss of coolant accident), LOHA (loss of heat sink accident) and LOFA (loss of flow accident).
LOCA accident from loosing coolant is derived. This event in PWR reactor can occur through breaking in primary loop of reactor either hot leg or cold leg.
In this state the existent water in the primary loop along with steam are strongly leaked that blow down event occurs. But when the lost water of primary loop through RHRS (Residual and Heat Removal System) including HPIS (High Pressure Injection System), LPIS (Low Pressure Injection System) and Accumulator (passive system) are filled this process is entitled Refill. When the primary is filled and all the lost water in it is compensated then the reactor sets in the normal status.
In the blow down status the raised steam is due loss of pressure in primary loop and moving the situation of reactor's primary loop from single phase to two phase flow.
LOHA accident from loosing heat sink in reactor is derived. Heat sink is as steam generator in nuclear reactor. This event when occurs heat exchanging between primary and secondary loops are not done.
This event might occur through lacking water circulation in the secondary loop of reactor. Not circulating the water might occur through closing either block valve (which sets after demineralizer tank) or other existent valves in the secondary loop.
LOFA accident from disabling and loosing the pumps in either primary loop or secondary is derived. In case either primary loop's pump or secondary encounter with problems LOFA accident might occur.
In the secondary loop the main feed water pump has duty of circulation of feed water to steam generator and sets after condenser pump.
In view of dynamically stability of a nuclear reactor, there is a stable system so that an excess reactivity is injected to it and it able to be stabled again in shortest time. The stability of linear systems in the field of complex numbers by defining the polarity of closed loop transfer's function is determined.
In case all the polarities are in the left side of imagine page then system will be stabled. In the time field the stability definition means system's response to each input will be definitive. In the matrix form all the Eigen values of system have real negative part. www.intechopen.com

Dynamics of nuclear reactor
There are several methods for investigation of nuclear reactors dynamics.
One of the most important methods to study reactor dynamics and the stability of nuclear reactor is define of transfer's functions and application of it to analyze the closed loop function.
According to following figure a closed loop system including transfer function, feedback and related applied reactivity are shown: is: error reactivity in frequency field that is as input reactivity to transfer function, () Gs is: transfer function, () Hs is: feedback function and () ns is: output of closed loop conversion function that means the density of neutrons.
There is also: According to Fig.1 for both transfer function and feedback function existing in closed loop can write: and it can also be written: In order to survey the stability of a closed loop system the term of [1 ( ). ( )] Gs Hs + must be set zero and by solving this equation, all the roots that are as zero and pole for closes loop system, will be defined. The stability condition of a closed loop system is lack of positive real part of poles. It means all the poles must be the left side of real-imagine graph.
Reactivity feedback causes the steady operation of nuclear reactor and equilibrium of its dynamical system.
A transfer function can be either linear or not. Each system variable can be affected as an input reactivity to transfer function as shown in Fig

System's variables of nuclear reactor
There are several controlling factors in nuclear reactors such as: If each system variable as a mathematics variable is considered then can write: and: Where: By taking the laplace conversions from two sides of above equations can write: and: So two last equations that are based on Laplace conversions can be converted to following form: and: Therefore: and: In order to define the transfer function, it will be deduced as shown below: Where: and:

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 K eff is: 0 and maximum of it is: υ namely: 2.43.
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]: and: If the transfer function of G(s) as V function is assumes then: www.intechopen.com and: () () .
Firstly the transfer function is supposed. Secondly due to mostly reactor's cores are twin therefore once a signal with delay time ( d τ ) 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.

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: 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 k τ → , 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): and: 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 www.intechopen.com 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: and: Where: () nt  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, C i 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: () nt  , ρ e , () nt , () Ct  and () Ct are considered the these can write as following: and: If these both () nt δ  and ( ) Ct δ  equations to matrix format are written, then the matrix format of them is shown as below: 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: Therefore through () G ρ , the parts of ( ) G ρ  will be defined and the stability condition through V function is applied.
www.intechopen.com In large twin reactor's cores, it can be written: It is also supposed: and: 12 21 αα = (56) q 1 , q 2 and ρ i are respectively as following: According to Newton's low of cooling can write: 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: www.intechopen.com and: Then the differentials of these system variables are respectively as following: and: If these system variables differentials are to matrix format are written, then can write:

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).
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: This function has features as following: 1.
() Vx  has definite positive value. 2. The partial differential of () Vx is continued. 3. Where the ρ is momentum, () V ρ is negative quasi relatively.

The ()
Vx function can be written as following: It can also write: 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: If: 0 V =  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: and:

Analyzing the theory by mathematical model
In this work the value of K eff 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: Supposition Where x is: absolutely descending, 0 x is: the initial value of x and t D is: the innate delay time.
The SIMULINK of MATLAB is an appropriate software to analyze the performance of this simulation (Tewari, 2002 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: t D . The parameters which must be adjusted are: Set Point that is: the default amount of K eff as reference K eff and the meaning of the Set Point=100 is: K eff =1, the velocity of control rod (v), recent K eff (block H) and the stop time that is: the innate delay time or t D . The graphs can be observed by the oscilloscope.

Results
In this simulation the input parameter value (H) is attributed to arbitrary K eff . So if the favorite K eff is: 1 then the value of H will be defined: 100 and it is also for Set point (reference K eff ). So Set Point: 50 means the reference K eff is: 0.50 and in this situation this K eff is enumerated as the arbitrary and favorite K eff that stability of reactor in this situation is based on it. This arbitrary K eff with output of Zero-Order Hold block is compared for performance of simulation by SIMULINK software. Also for example velocity: 3 means the velocity of control rod in this simulation is: 3 units per second (for example: 3mm/s). In this state the velocity of control rod is increased comparing to the last stage which was: 1 unit per second. The velocity of control rod belongs to Speed block which is an input to Tsp Sum block in the block diagram.
By changing the values of the cited parameters (which were: Set Point, the velocity of control rod (v), H and the delay time), the different states of the graphs can be shown according to Figs.4, 5, 6 and 7:   The Figs.4, 5, 6 and 7 show if the velocity of control rod for upward and downward movements is increased then the K eff will be pendulous surrounds of defined Set Point (reference K eff ) and also the oscillation amplitude will be more than lower velocity situations. For low velocities the oscillation amplitude is slight and acceptable. Another effective factor is inherent delay (such as derived delay of control rod mechanism) and its inordinate increasing can cause unstable states.

Conclusion
By this simulation the best response and operation which a reactor can have from stability aspect, according to its control rod velocity is derived.
According to Figs.4-8 the best status in the Fig.4 is observed in which there is no vibration in response. In this simulation the stability of reactor depends on either velocity or delay time values directly because delay time plays a key role. Therefore in this simulation the admissible ranges of velocity and delay time which can be caused to stable the reactor are respectively: low velocity of control rod around 1-3 units per second and short delay time (10ms). However in this case reach critical state (Keff=1) for nuclear reactor will be taken more than modes Figs.5-8. As in all the figures is observed, so can deduce the velocity of control rod plays the more important role than delay time of mechanism in stability of nuclear reactor. Whereas for minor and major changes of reactivity and shut down of reactor in emergency states, there are some kinds of control rods at nuclear reactors cores such as regulating rods, safety rods and shim rods, therefore this simulation can be applied for each control rod in either LWR nuclear reactor or research reactor cores which have vertical control rods. Also this simulation can be applied for each batch control rods which