Spatial Reactor Dynamics and Thermo Hydraulic Behavior Simulation of a Large AGR Nuclear Power Reactor in Response to a Reactivity Step Change Disturbance

In this article, two-dimensional partial differential equations with time representation of nuclear power reactor kinetics are considered for spatial reactor dynamics and thermo hydraulic behavior analysis of a large thermal advanced gas cooled reactor (AGR) type used for nuclear power generation. The equations include the neutron flux equation and delayed neutron precursor concentration, together with taking into account the equations to represent the thermo hydraulic behavior of the fuel, coolant and moderator temperatures. These equations are solved numerically using the finite difference method. For time propagation, an implicit method is applied. The desired initial condition for the reactor to stay at stable critical condition is established by finding the correct value of reactivity. The reactivity disturbance effect in the reactor is studied for different cases and presented for high reactivity values. The model was developed for the analysis of a large AGR with 2000 MWe for future power generation. The results show that the model not only behaves stably but also predicts the results physically for all the various parameters.


Introduction
Nuclear power stations play an important role in electricity generation in industrial countries.Nuclear energy is a clean and cheap energy source after renewable energy sources, such as solar and hydraulic energy.Thus, for industrial development, progress in nuclear power reactors is essential.Advanced gas cooled reactors (AGR) and the next generation high temperature gas cooled reactors (HTGR) are safe with respect to their moderator and cooling materials and can generate a high amount of energy; these types of nuclear reactors are desired for high capacity nuclear power plants.
Nuclear reactor simulations are conducted for analysis of reactor behavior.There are some computer codes that have been developed by universities or research establishments in order to predict flow regimes, hydro dynamical instability, reactor kinetics, etc.These computer codes (analyzers) include MINCS, PHOENICS, TRAC, RELAP5, and others, and they not only predict the optimized criteria for the thermo hydraulic condition for nuclear reactors but also fill in the gaps between the scarce experimental results that have been found by different research sources and obtained under substantially different conditions.These codes have some capabilities and limitations that have been discussed by Wulff [1] and Physical Benchmarking Exercise [2].RELAP5 has been used massively for nuclear reactor modeling.Comparison of results obtained by RELAP5 modeling and experimental data collected by Groudev [3] for pressurized water reactors (PWR) of VVER440/V230 type showed that the analysis using this computer code is valid.Because such codes are not easily accessible, especially for developing countries, modeling of the nuclear reactors with more accurate mathematics and higher physics is required to predict more precise simulation results.
Arab-Alibeik et al. [4,5] modeled a PWR reactor.Their model is one-dimensional and includes five differential equations.They include the neutron flux equation without a diffusion term, the six group delayed neutron equation, energy conservation equations for fuel temperature and an energy conservation equation for the cooling flow.They also used fuel and moderator feedback coefficients in the model.Marseguerra et al. [6] used the point kinetics neutron flux equation with thermal equilibrium equations for the reactor core.In their model, the neutron diffusion term was omitted from the neutron flux equation.They considered six groups for the delayed neutron equations.Weyfeng et al. [7] presented a three-dimensional model that based on a 200 MW mina tour reactor with a thermal pool reactor.Their model included the neutron flux equation, one group equation including a neutron diffusion term, delayed neutrons, the cooling flow thermal equilibrium equation and one state equation for the effect of control rod replacement on the neutron flux.Sadek et al. [8] used a one-dimensional second order parabolic equation to model the nuclear reactor, including neutron flux but without delayed neutrons and thermo hydraulic terms.
In this article, the two-dimensional neutron flux equation is considered for the neutron spatial distribution together with equations to represent the dynamic behavior of the fuel, coolant and moderator temperatures.In this model, which has been applied to an AGR type reactor, the cooling temperature is different than the moderator temperature and is included in the equations definition.This model has one more differential equation with respect to the PWR model that was presented by other researchers.The present model will consider symmetric and anti-symmetric reactivity inlet disturbances and their effect on nuclear reactor behaviors.

Governing Equations
Nuclear reactor kinetic equations are declared by neutronic and thermo hydraulic equations to describe the dynamic behavior of the reactor core.Nuclear reactor neutronic equations include prompt and delayed neutrons.The neutron flux differential equation is as follows (Stacey [9] and Glasstone, et al.
The one group delayed neutron precursor concentration equations are In the one group delayed neutrons equations, the six groups of delayed neutrons are averaged, and only one variable is used instead of six variables.The average one group delayed neutron and one-energy group equation is written as To solve the above equations, the variables are normalized as The normalized equations of the neutron flux and delayed neutron become as follows: P is the normalized power, and the time constant of normal power and migration area are obtained as Reactor reactivity will be affected by the fuel and moderator temperature feedback coefficients.The reactivity feedback formula is as follows: Fuel, cooling and moderator temperatures can be found from the following differential equations: The above equations can be derived by considering the thermal equilibrium between the fuel, cooling and moderator temperatures.The fuel and moderator temperatures are varied by the power that is released from the nuclear energy generated from neutron production.The fuel rods and moderator portion from the generated power E 1 and E 2 , respectively, as follows: of the inlet and the outlet temperature of the reactor The mass flow rate is determined in such a way that the cooling temperature is constant, which can be obtained by assuming that the time variation of the cooling temperature is equal to zero.The following equation can be obtained for the cooling flow rate:

Initial and Boundary Conditions
To solve equations, the initial condition for different variables must be determined over the whole reactor cross-section area.The power distribution in the reactor at the initial condition is assumed to be a sinusoidal form, as shown on Figure 1.The reactivity distribution is found using the following equation:  is a constant and is called the convergence coefficient.Its value can be found by trial and error.When the right values of reactivity are found, it should be checked in the model that the power stays at the initial condition without changing when the computer program is run for a long period of time.See Figure 2 for the normalized power distribution after a long period of time (steadystate condition).The other variable distributions were found with respect to the power at the steady-state condition and are presented in Figures 3 to 7. The variables values are high at the reactor center but decline to their minimum values near the reactor boundaries.

Solution Technique
To solve the governing equations numerically, the im- plicit method was applied for differentiation because it has higher numerical stability.The total differentiated equations were solved using the Three Diagonal Matrix Algorithm (TDMA) by considering the Alternative Direction Implicit (ADI) method.Two set of equations, including neutronic and thermo hydraulic equations, were considered.The neutronic equations include the normal power equation, delayed neutrons and temperature feedback on reactivity.The thermo hydraulic equations include the fuel, cooling and moderator temperatures and the cooling flow rate.The two sets of equations are affected mutually by the reactivity feedback equation, and all of the equations must be solved together simultaneously.A square cross-section nuclear reactor is assumed with 10 by 10-meter dimensions.The number of meshes in the x and y directions are 21.Mesh independency was also considered with a different number of meshes and increased up to 101 meshes in each direction.
The results generated from different runs of the computer program showed that the differences are negligible.Thus, the results obtained in the rest of the article depend on 21 meshes in each direction.
The second order partial derivatives of the neutron flux equations were converted into differentiated form.Each step of time was conducted in two half steps.In the first half step of the time, the neutron flux equation was solved implicitly in the x direction but explicitly in the y direction.In the second half of the time step, the same neutron flux equation was solved implicitly in the y direction and explicitly in the x direction.Because of the long mathematical manipulation and shortage of space, only the final form of the obtained differentiated equations is presented below.
The neutron flux equation in the first half of the step of time is The neutron flux equation in the second half of the step of time is The normalized delayed neutron precursor concentra- The equations representing fuel, coolant and moderator temperatures are, respectively, , 1 The required coolant flow varies as Equations ( 13) to (17) were solved for each step of tim .Results

.1. Whole Reactor Cross-Section Disturbed by
nce the correct values of the different variables were

Reactivity
O established, the reactivity value was increased by 0.01 mN (mN is the reactivity unit used in the British nuclear industry) over the entire reactor area after 0.1 seconds.Figure 8 shows the transient response of the reactor with respect to the reactivity change.The reactor power at the middle point increased to about 25 times of the initial state in just 0.7 seconds, then declined to 1.93 times the initial value after 1.7 seconds and then increased slightly.This behavior is the transient response of the reactor with respect to the disturbances.The stable behavior of the reactor starts after about 4 seconds.Figure 8 is plotted for 5 seconds to demonstrate the reactor behavior at an early time at the center of the reactor, and Figure 9 is the continuation of the Figure 8 behavior at a longer time.It also shows that the power increased gradually with time and reached 8.2 times from 4.6 times the initial values in 45 seconds.The sudden increase of the reactor power in Figure 8 was caused by prompt neutrons that were generated quickly by the reactivity disturbance response.The second increase in power was because of the delayed neutron generation.Figure 10 shows same behavior at the center cross-section of the reactor core.presents the reactor power after 50 seconds.It is interesting that, because of the sharp power profile in the initial conditions, the early time power profile is sharp (Figure 10 at a time of 0.5 sec), but, because of the high difference of neutron flux between the center and the boundaries, the power profile becomes almost uniform over the whole reactor cross-section (Figures 10 and 11 at 50 seconds).Figure 12 for an early time and Figure 13 for a longer period of time show the reactivity transient behavior after it was increased intentionally.The reactor power decreases and tends to stabilize the system because of the negative reactivity feedback coefficient effect.Figure 14 demonstrates the transient response of the fuel, coolant and moderator temperatures on the reactivity disturbance at an early time.The fuel temperature increment is greater than the other temperatures, and Copyright © 2011 SciRes.EPE ust .

One Node at Center of Reactor tivity
uring the previous section of the present paper, the from 430˚C in equilibrium it increased to 685˚C in j 2.5 seconds.The temperature behavior is also presented for longer periods of time in Figure 15.Because of the higher power production, the fuel temperature increased more with respect to moderator temperature.Because of the greater time constant of the moderator, the moderator temperature increase was lower than the fuel temperature.The cooling temperature is a function of the fuel temperature increment, but the cooling temperature increment was low because the cooling flow rate also increased.The cooling temperature reached 507˚C from 380˚C for the initial value in 50 seconds.The fuel temperature distribution in the reactor cross-section is presented in Figure 16 for different times, and the cooling mass flow rate response to the uniform reactivity disturbance is presented in Figure 17.
4 Cross-Section Disturbed by Reac D whole cross-section of the reactor core was disturbed by the reactivity.In this section, only one node was disturbed at the center of the core by a reactivity value of 0.1 mN (which is added to the initial value of the reactivity).The reactivity disturbance starts after 0.1 seconds.The reactor transient response for power is presented in Figure 18.Because of the high value of reactivity, the reactor power increased sharply, but after about 0.78 seconds its value decreased to 0.16 times the initial value.
The oscillating behavior of the reactor decreased to a stable value after around 15 seconds.The long-term behavior of the reactor power is presented in Figure 19.
The reactor power oscillation in Figure 18 was caused by the fast generation of the prompt neutrons generated  reactivity disturbance.However, by the high amount of the high value power production was damped and decreased by delayed neutron production.The reactor power behavior was not only affected by the prompt and delayed neutron productions and leakage but also controlled by positive and negative reactivity feedback coefficients due to the fuel and moderator temperatures.It should be mentioned that the large value of the reactivity introduced into reactor system is dangerous to the reactor because it will generate very high power and fuel, moderator and cooling temperature increases, as shown in Figure 19.Such a high value temperatures will melt down the reactor and possibly create a reactor explosion.Our model cannot predict the melt down or reactor explosion (the equations do not contain the complicated condition of the melt down or reactor explosion correlations).However, it is intended to demonstrate and show the responses of different parameters for a high value of reactivity when put into operation intentionally.Then, the reactor responses can be analyzed, and the model capabilities and limitations can be considered.The reactor responses were not only characterized physically (qualitatively) but also checked mathematically without any failure of the model developed (i.e., it remained stable numerically).The results showed that, even for a high reactivity, the model stayed stable numerically without facing hard failure of the solution.

Conclusions
wo-dimensional t T a were nalyzed for nuclear power reactor kinetics and thermo hydraulic behavior for large AGR nuclear power generation.The governing equations used in modeling of the reactor included neutron diffusion, delayed neutron precursor concentration and thermo-hydraulic equations for fuel, moderator and cooling temperatures.The equations were solved numerically using the finite difference method.The reactor transient response was demonstrated for reactivity disturbances in different conditions.For this reason, the initial reactivity was established for the steady-state condition.The initial reactivity values are the reactivity amounts that do not change the reactor power and maintain all the variables at steady-state values.Higher values of the reactivity were imposed as follows.
1) A nuclear reactor AGR type with 2000 MWe was develo ved well physically (qualitatively) and numerically during critical conditions without failure.
2) The reactivity of all nodes was increased by 0.01.The reactor response showed that even th les increased, however, the negative feedback coefficient and delayed neutron generation prevented the reactor system from the critical condition.Nevertheless, if this condition were continued, for a reactor without a control system, the fuel rod temperature would increase, and the reactor would end with crisis.
Only the center node reactivity disturbance increased to a value of 0.1.The reactor response mperature increased to very high values.This behavior could lead to fuel rod melting and explosion, even though this model is not prepared to predict these crises. Fission constant of one-group delayed nu  i Fission constant of delayed nucleons of group i (s −1 )

Figure 1 .
Figure 1.Initial distribution of normalized power at the reactor core cross-section.

Figure 2 .Figure 3 .
Figure 2. Normalized power distribution at the reactor cross-section after using the right values of reactivity.

Figure 4 .
Figure 4. Initial values of fuel temperature distribution across the reactor core.

Figure 14 .Figure 15 .
Figure 14.Transient behavior of fuel, coolant and moderator to reactivity step disturbance at early time.

Figure 16 .
Figure 16.Fuel temperature distribution across the reactor core due to reactivity step change at different times.