Linear and nonlinear stability in nuclear reactors with delayed effects

Abstract. The research of a nuclear reactor model has been observed, where the system consists of two differential equations with one delay. A linear analysis has been performed, the asymptotic stability model of the area D0 and D2 has been defined, in which a stable periodic solution of one frequency appears. In the nonlinear analysis the analytical expression of the solution is presented with the help of bifurcation theories. In the numerical experiment using the scientific simulation program “Model Maker” numerical Runge–Kutta IV series method asymptotically stable solution and a stable periodic solution has been received and compared to the stable periodic solution received in nonlinear analysis with the help of bifurcation theories.


Introduction
Nuclear (atomic) reactor is a device to initiate and control a sustained nuclear chain reaction.In the nuclear reactor uranium 235, plutonium 239 and uranium-238 are controlled by the time and quantity aspects fission into lighter atoms, releasing large amount of neutrons which have much kinetic energy.These neutrons, gradually impacting the reactor walls or nuclei borders transform kinetic energy into thermal energy.Therefore, the nuclear reactor walls are retarders (moderators) of the atoms in the nuclei fission reactions and transformers of the neutron from kinetic energy into thermal energy [1].
In reactivity transient analysis, the heat transfer from the fuel to the coolant is crucial.This heat transfer is significantly affected by the behaviour of the gas gap between the fuel and the cladding, but also by the heat transfer conditions at the cladding surface.It affects, on one side, the fuel temperature, being a significant safety and neutron kinetics feedback parameter.The heat transfer crisis can occur at the cladding surface leading to the wall superheating, which results in the cladding destruction [2].
When the heat transfer situations include extremely high temperature gradients, extremely large heat fluxes, or an extremely short transient duration, the heat propagation speeds are finite, and the heat conduction mode is propagative and non diffusive.
Although of radial fuel element the thickness of the gap is quite small, the low thermal conductivity of gases causes a large temperature drop across the gap.Because the gap spacing of radial fuel element is not uniform and the heat conduction transfer process is very complex.
For the nuclear fuel case, the physical meaning of this term suggests that there is an effect of heat transfer due to the rapid change of the source term, i.e., the neutron power, and this effect occurs in nuclear reactors, due to changes in the disturbances reactivity of different physical effects, mainly by Doppler [2].
According to the purpose, nuclear reactors are subdivided into energy, isotopic and research [1].
Since the linear theory dass not sufficiently evaluate nuclear reactors stationary mode stability, therefore, adjudicating stability problems are essential to investigate the nonlinear mathematical model of the nuclear reactor [3].
Whereas nuclear reactors are objects with distributed parameters and their sufficiently strict mathematical models are expressed in the nonlinear differential equations with partial derivatives or to certain additional simplifications -equations with a delayed argument [3].Each nuclear reactor unit can be described in various types and complications of differential or integration-differential equations [4].
The difficulty of a nuclear power reactor's stationary modes analysis is associated with a modern reactor design complexity, a large variety of different processes in the reactors and interrelationships.Determining dynamics of the reactors stationary conditions, a research of various theoretical characteristics is relevant, e.g.temperature, power, a nuclear fuel cycle [3].
Instability analysis in boiling water reactor (BWR) plays a central role in the understanding of the physical mechanism that induces the observed power oscillations, which are unstable and occur at low-core flow and relatively high-power conditions [5].
Complex heterogeneous processes tangle in the reactor and form internal reversible links.According to certain equations the feedback connections are described and divided into: mathematical, concentrated, distributed and hybrid.Neutron density changes in the reactor lay by heat amount exude changes, as a result various environments and reactor's design elements of temperature and density change.When these feedbacks are sufficiently strong, reactor parameters change has the required phase offset, hence reactors stationary conditions are unstable.Additional feedbacks, determined by management and regulatory organ also may be the reason for the reactor's unstable operation.Instability of the stationary conditions can contain different forms, where the most common are nonperiodic and periodic [4].
We studied reactor used for scientific research.We are analyzing the point model of a nuclear reactor dynamics, which ignores the effect of delayed neutrons, however considers the power and temperature reactivity coefficients and the delays while changing the reactivity of thermal convection.The system of the model consists of differential equations with one delay, which has been partly analyzed [13].
There N (t) is the power of the reactor.N 0 is its stationary value.T (t) is the temperature deviation from the stationary value.A is the coefficient of the temperature reactivity.B is coefficient of power, k is a positive constant, h is a steady time delay, t is time.We analyze system (1), when k = 1, h = 1.

Linear analysis
We consider the stability of stationary solution For this purpose we use D-decomposition method [25].
We performed an amendment have in the system (1) We obtain a system of differential equations where the linear part is system of differential equations Let us note a = AN 0 , b = BN 0 .Then the characteristic quasipolynomial of system ( 4) is defined as Equation ( 5) has a zero root λ = 0, when Line ( 6) becomes one of the D-decomposition curves on the plane ab.Let us conceive that quasipolynomial (5) has a purely imaginary root λ = iσ.After the set of λ = iσ in Eq. ( 5) and the separation of real and imaginary parts, we get the remaining parametric equations of D-decomposition curves The analysis of formulas ( 6) and ( 7) allows us estimate, that when σ → 0, gets coordinates of reversible point (0.5; −0.5).Furthermore, the second equation of system (7) follows, that b = 0, when σ sin σ − cos σ = 0, σ = 0.As a result we obtain the equation Equation ( 8) is valid only when σ ∈ (πk; π/2 + πk), k = 0, 1, 2, . . . .With a direct calculation we estimate that ( 8) is legitimate at σ ∈ (0; π/2], when σ ≈ 0.86 radians, then a ≈ 1.13.Taking into account the periodicity of functions sin σ and cos σ, we obtain values of parameter a, where b = 0. Further, considering equality (7) we estimate, that increasing the parameter a, the parameter b also increases, when σ ∈ (0; π/2).
Furthermore it is necessary to emphasize that: Subsequent calculations pointed out that π/2 + α < σ π, a < 0, b < 0, parameters a and b are growing together.However, when π < σ 3π/2, then negative value a matches positive value b, and a decreases, when b increases.
Taking into account the periodicity of functions sin σ and cos σ, also considering conditions (8), as well as the previous reasoning we estimate curves of D-decomposition.
Figure 1 shows D-decomposition curves on the parameters a and b in the plane.We are to find the area of the asymptotic stability D 0 , which corresponds to quasipolynomial, where real parts of all roots are negative [25].For the set of area D 0 , it is sufficient to find at least one point that corresponds to the quasipolynomial, where real parts of all roots are negative [26].
From solution (2) of the stability point a = 0, b = b * > 0 follows that the area D 0 is the asymptotic stability domain, because quasipolynomial (5) at this point has only two roots in λ 1 = −b * , λ 2 = −1, which real parts are negative.
We are to analyze the characteristic quasipolynomial (5) in the case when b is fixed and the parameter a > 0. We are to ascertain the case, when the parameter a increases two purely imaginary roots ±iσ 0 that appear in the characteristic quasipolynomial (5) for the first time, while the remaining roots have negative real parts.It is, when where σ 0 is only root of equation ctg σ = (σ 2 − b)/((1 + b)σ), that belongs to the interval (0; π/2 + α), where α < 0.43 is radians.This follows from the preceding reasoning.
Lemma has been substantiated.
Thus, transitioning parameters a and b from the area D 0 to area D 2 , two complex conjoint roots λ(ε) and λ(ε) with positive real parts appear.

Nonlinear analysis
We are to consider the system of nonlinear equations (1), when a = a 0 + ε.After the use of substitution (3), we will replace the system of differential equation with one of the second series with a delay: where Let us suppose, that Since τ 0 > 0, differential equations with a delay solution construction method can be applied to the differential equation ( 9) [22].
On the ground she, in Eq. ( 9) the substitution t = (1 + c)τ , (|c| < 1), we be fixed by the time and record lines ( 10)- (12) in Eq. ( 9), we obtain, that the differential equation ( 9) will become a formal identity The obtained identities of the left and right sides located ξ degrees and equated the coefficients to equal ξ the degrees, we get differential equations, which after appropriate modifications are: Equation ( 14) is a second series linear of a not homogeneous differential equation with steady coefficients.Therefore, we are looking for the solution in the form After inserting solution (16) to Eq. ( 14) we obtain where Alternatives of Fredholm [26,27] follows, that the differential equation ( 16) has a periodic solution when these equalities are satisfied There f 3 (τ ) is equal to Eqs. (15) right side.
From conditions (17) we get a system with two linear equations in regard to variables b 2 and c 2 : where

Numerical experiment
When the qualitative analysis of model ( 1) had been performed, we obtained the results using the numerical experiment.We are to search the differential equation with one delay numerical solution using Runge-Kutta IV series method.Asymptotically stable solution of the differential equation ( 9) obtained by the numerical method is presented in Fig. 2, under the following parameter values: a = 0.9; b = 0.34; h = 1; x 0 = 0.9; y 0 = 0.8, where a and b values are taken from the domain D 0 .
In Fig. 2 we observe, that the temperature coefficient of reactivity a and coefficient of power b selected from domain D 0 , is the temperature deviation of the nuclear reactor, where the stationary value is asymptotically stable.Moving neutrons at different reactor's height are increasing its temperature.The reactivity of temperature together with the power of reactor asymptotically approach to the equilibrium state.
A stable periodic solution of the differential equation (9) obtained by the numerical method is represented in Fig. 3, under the following parameter values: a = 1.2; b = 0.04; h = 1; x 0 = 0.9; y 0 = 0.8, where a and b values are taken from the area D 2 .
Figure 3.We observe that the temperature coefficient of reactivity a and coefficient of power b selected from domain D 2 , is the temperature deviation of the nuclear reactor where the stationary value is a stable periodic solution.The negative temperature of reactivity affects the power o the reactor's stabilization [16].As the temperature solution is periodically stable, the reactor's power has a stable periodic mode.
We represent the stable periodic solution (continuous line) of the differential equation (9) obtained by a numerical method and a stable periodic solution obtained (21) by the formula (dotted line).
Figure 5.We compared the stable periodic solution of the temperature obtained by a numerical method and the stable periodic solution of the temperature obtained by analytical methods.

Conclusions
After performing the linear analysis of the reactor mathematical model described by two differential equations with one delay (1) by the D-decomposition method were obtain the asymptotic stability area D 0 and area D 2 , which appears a stable periodic solution of one frequency.Nonlinear analysis of model has been performed using the bifurcation theory and the analytical expression of a stable periodic solution has been constructed in the area D 2 .Sustaining obtained results with a simulation package "Model Maker" in the linear analysis of the numerical experiment, a stable asymptomatic solution and a stable periodic solution have been received with the help of a numerical method, which acceptably very well with the obtained stable periodic solution in the nonlinear analysis constructed with the help of bifurcation theories.
www.mii.lt/NA 2 Dynamics model of the nuclear reactor

Fig. 2 .
Fig.2.The asymptotic stable solution received with the help of numerical method.

Figure 4
Figure 4. We obtained a stable periodic solution of the temperature derived from formula(21).We represent the stable periodic solution (continuous line) of the differential equation (9) obtained by a numerical method and a stable periodic solution obtained (21) by the formula (dotted line).Figure5.We compared the stable periodic solution of the temperature obtained by a numerical method and the stable periodic solution of the temperature obtained by analytical methods.

Fig. 3 .
Fig.3.The stable periodic solution received with the help of numerical method.

Fig. 4 .
Fig. 4. The stable periodic solution received with the help of equation (21).