On the stability of the point reactor kinetics equations

Nuclear Engineering and Design 240 (2010) 1443–1449 Contents lists available at ScienceDirect Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes On the stability of the point reactor kinetics equations Germán G. Theler a,b,∗ , Fabián J. Bonetto b a b Instituto Balseiro, Centro Atómico Bariloche, Río Negro, Argentina TECNA S.A., Buenos Aires, Argentina a r t i c l e i n f o Article history: Received 22 October 2009 Received in revised form 30 January 2010 Accepted 1 March 2010 a b s t r a c t As the basic neutronic problem is unstable by nature, maintaining a reactor critical is a task that requires a lot of effort. This work presents and discusses some aspects related to the stability of the basic physics behind a nuclear reactor core based on the well-known point reactor kinetics equations. First, the linear non-feedback case is studied where differences between Lyapunov and BIBO stability are found. These differences are shown both numerically and analytically, and explained using a reactor physics based reasoning. Finally, a simple model is used to analyse the intrinsic stability of the point reactor equations when reactivity feedbacks are taken into account. A method for constructing conceptual stability design maps is proposed and a basic interpretation of the simple results obtained is given. © 2010 Elsevier B.V. All rights reserved. 1. Introduction The subject of the stability of a nuclear reactor is of particular concern in a conceptual engineering stage, as shown for example by the state-of-the-art generation IV design paths (Zanocco et al., 2003). The mechanisms of instability appearance are diverse, being the most common the coupling between neutron dynamics and the thermalhydraulics of the coolant loop (Theler, 2008; Doutta and Doshi, 2009). A frequency domain analysis can be made in order to study theoretical instabilities (Peng et al., 1984, 1986; Garea et al., 1999), but other linear and non-linear methods also exist (Lahey and Podowski, 1989). On the other hand, best estimate codes may be used to simulate reactor transients. Not only are these neutronic-thermalhydraulic coupled codes used nowadays to evaluate the plant state after postulated accidents (Salah and D’Auria, 2007, 2009) but also as a tool to assist the designer in the decision making process. In this paper we propose a method for analysing the stability of the point reactor kinetics equations, in order to grasp the general behaviour of the system and to generate some conceptual design maps. The development is based on approximations and general facts that are aimed at maintaining the mathematics fairly simple without entering into detailed and precise calculations. Particularly, we address and study the stability of the wellknown point reactor kinetics equations. This model might be a very crude simplification of a real reactor core, but it has two main advantages that make it suitable for our purposes. The first is that point kinetics equations maintain the underlying physics ∗ Corresponding author at: TECNA S.A., Encarnacion Ezcurra 365, Buenos Aires, Argentina. Tel.: +54 (11) 4347 0300. E-mail address: [email protected] (G.G. Theler). 0029-5493/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2010.03.007 of the general transient neutron problem, and the second is that they have a reasonable and simple mathematical formulation. We ought to remark that the facts and results shown here only apply to the mathematical model and equations discussed and they are not directly applicable to real reactors. 2. The point kinetics equations Every stability study presented in this work is based on the point reactor kinetics equations, whose derivation and particular analytical solutions may be found in the standard reactor physics literature. In particular, Duderstadt and Hamilton (1976) provide a rather intuitive approach to the equations while Henry (1975) shows how they can be rigorously obtained from the timedependant diffusion equation. The resulting equations are  dn −ˇ i ci + q = n+ dt  I dci ˇ = i n − i ci dt  i=1 (1) The variables and parameters that appear in these equations are the result of integrating the original partial derivatives formulation of the neutron conservation problem with respect to both the spacial and angular coordinates and to all available neutron energies. Nevertheless, these particular details are not of interest for the analysis that follows. The nomenclature used is the one that may usually be found in the literature and a brief description of each kinetic parameter is shown in Table 1. We take each one of these as a constant and a priori known parameter. To complete the formulation of the problem, initial conditions n(t = 0) and ci (t = 0) should be given. In the absence of an exter- 1444 G.G. Theler, F.J. Bonetto / Nuclear Engineering and Design 240 (2010) 1443–1449 Table 1 Description of the variables and parameters appearing in the point kinetics equations. n ı ın ˇ ı + − ın + ıc    ˇ ıċ = ın − ıc  Neutron population a Number of neutron precursor groups i th precursor group concentration −1 Total reactivity (keff − 1) · keff i th precursor group fission yield Total  fraction of delayed neutrons ˇi ) ( i th precursor group decay constant Average time between the emission of neutrons born in two consecutive generations b Independent neutron source intensity n I ci  ˇi ˇ i  q ıṅ = nal neutron source, the equilibrium condition for a steady-state population n is ˇi  n i  (2) Note that if q = 0, an homogeneous system of equations is obtained, and n and ci may be scaled by an arbitrary factor without affecting the solution form. As the thermal power is proportional to the neutron flux, we may write the kinetics equations as a function of the power P(t) instead of the neutron population in some sections of this paper. 3. Stability studies without reactivity feedback Again, with the objective of looking at the physical background of the problem and not solving for a particular detailed solution, we take only one group of precursors instead of the traditional six groups. This approximation simplifies the associated mathematics and at the same time maintains the physical basis of the original problem. Let us consider the point kinetics equations with one neutron precursor group and with no external source, written as a dynamical system with a state vector x   ẋ = ṅ ċ  = ( − ˇ)/ · n + c ˇ/ · n − c  (3)  − ˇ  n + c   ˇ  n − c  0=  0= (4) As both n and  are positive, the reactivity must be zero for a fixed point to exist. This is the well-known result that a reactor is critical if and only if keff = 1 and, on the other hand, justifies the equilibrium condition shown in Eq. (2). Let us now consider what happens around the fixed point x re-writing the problem using the perturbation variables ın = n − n ıc = c − c  ı =  −  =  There are two different stability concepts that may be of interest for the nuclear reactor core designer. One is the so-called Lyapunov stability which refers to what happens when the vector state is perturbed near a fixed point in phase space. The other one is the bounded-input bounded-output (BIBO) stability concept that copes with perturbations on the external input of a given dynamical system. In order to proceed with the standard linear approach, we may neglect ıın with respect to the other terms—just for this section—and consider ı(t) as an external control action. Thus, the equations are casted in a matrix based form as  ıṅ ıċ   = −ˇ/ ˇ/ −   ın ıc   + n / 0  ı (6) To study the Lyapunov stability of the fixed point x we may compute the eigenvalues of the system matrix. They are s2 = − s1 = 0 ˇ −  (7) As there is one eigenvalue equal to zero, according to the Hartman–Grobman theorem (Solari et al., 1996), nothing can be said about the linear Lyapunov stability of the fixed point x of the dynamical system (6). We will return to this subject on Section 3.3. On the other hand, to study how the system responds to perturbations on the external input, we may compute the transfer function of the neutron population n with respect to the reactivity . Laplace-transforming both Eq. (6) we get   n ˇ ın + ıc + ı   n ˇ ıR(s) s ıN(s) − ın(0) = − N(s) +  ıC(s) +   (8) ˇ ın +  ıc L{ıċ} = L  ˇ s ıC(s) − ıc(0) = ıN(s) −  ıC(s)  (9) −  By summing both sides of these equations we obtain the necessary condition for the system to be in equilibrium n =0  3.1. Intrinsic linear stability L{ıṅ} = L There is a fixed point x = x when the derivatives vanish  · (5) The system (5) is mathematically equivalent to Eq. (1) with q = 0. There is only one explicit non-linearity in the crossed term (ı ın)/. Thus, the point kinetics equations are linear only for a constant reactivity. In the following two sections we take the usual linear analysis approach of neglecting this second order term, while we will come back to study its influence in Sections 3.3 and 4.1. a All variables come from a weighted integral of the continuous diffusion equation. Only when the weight function is identically one, n is equal to the total number of neutrons. b If the reactor is critical,  is the average lifetime of a single neutron. ci (0) = Replacing in the system of Eq. (1) and taking into account that x is a fixed point, we obtain  To derive the open loop transfer function we set the initial conditions ın(0) = ıc(0) = 0 and solve for ıC(s) in Eq. (9) ˇ ıN(s) −  ıC(s)  ˇ C(s) = N(s) (s + ) s ıC(s) = By replacing C(s) in Eq. (8) we arrive at s ıN(s) = − ˇ ˇ n ıN(s) + ıN(s) + ıR(s)  (s + )  G.G. Theler, F.J. Bonetto / Nuclear Engineering and Design 240 (2010) 1443–1449 1445 In Fig. 2 we show the Nyquist plot that results when mapping the imaginary axis s = iω excluding the origin using s = re−i for −/2 <  < /2 with r = 10−3 for the test function KG(s) = Fig. 1. Closed control loop with a proportional compensator action on the linearised point reactor equations with a group of precursors. so we can solve for the linear transfer function ıN(s)/ıR(s) between the reactivity and the neutron population G(s) = n (s + ) ıN(s) = ·  s(s +  + (ˇ/)) ıR(s) (10) whose poles are equal to the eigenvalues (7). Nevertheless, in this case we may conclude that the point kinetic equations are not BIBOstable as the pole in the origin of the complex plane acts as an integrator of the excitation (t), modulated by the factor n /. For a positive reactivity the neutron population would monotonically increase in time, and vice versa. The same conclusion may be reached by noting that the term (n ı)/ in the first equation of the system (5) has the functional form of an external neutron source—or sink if negative—that is proportional to the reactivity perturbation. Therefore, as there is at least one bounded external perturbation—namely a constant finite positive reactivity—that produces an unbounded output, the system is not BIBO-stable. 3.2. Linear stability with a control loop In the previous section we studied the intrinsic stability of the point reactor equations without any feedback on the external reactivity. Without any of these feedback effects, the condition of criticality  = 0 would never be obtained as the probability of the matrix of the multigroup diffusion problem having an eigenvalue keff exactly equal to one is zero (Henry, 1975; Stacey, 2001). So we now turn out attention to study the linear stability of the point reactor equations with a closed control loop. Consider the control scheme depicted in Fig. 1. The open loop transfer function is KG(s) and the closed negative feedback loop transfer function is KG(s) Gcl (s) = 1 + KG(s) (s + ) s(s +  + (ˇ/)) with ˇ = 7.65 × 10−3  = 7.59 × 10−2 s−1  = 1.76 × 10−4 s that approximately represent the kinetic parameters of an open pool research reactor operating at full power with an equilibrium core. We see that the Nyquist plot does not circle the point (−1,0), and as G(s) does not have any zeros in the right half plane, the system is stable for the proposed gain. Moreover, the system is stable for any positive K—provided the feedback loop is negative—as the gain effect is just to scale the diagrams shown in Fig. 2. Thus, even though the point reactor equations are not BIBO-stable, the simplest continuous compensator possible acting on the reactivity is able to render it stable (Lewins, 1978). However, we must keep in mind that the unconditionally stability we found applies only to a linearised version of the point reactor equations that does not take into account any temperature, void or external reactivity feedback, and that the (continuous) controller is applied instantly and directly to . In any real case, there may be feedbacks, non-linear effects (real control rods), discrete sampling times and phase lags that would reduce the stability margin, possibly rendering the system unstable. 3.3. Non-linear effects In Section 3.1 we stated that nothing can be concluded about the Lyapunov stability of the fixed point x = x using the linear theory because the system matrix had a null eigenvalue. We now try to determine the Lyapunov stability by numerically solving the kinetics equations and interpreting the results using reactor physics theory. To see if a small disturbance in phase space makes the system remain close and eventually converge back to the point x , we consider again the equations with  = 0:   ẋ = ṅ ċ  = −ˇ/ · n + c ˇ/ · n − c  (11) Fig. 2. Nyquist plot of the linearised kinetics equations transfer function G(iω). The infinite closure was done excluding the origin making s = re−i for −/2 <  < /2 with r = 10−3 . 1446 G.G. Theler, F.J. Bonetto / Nuclear Engineering and Design 240 (2010) 1443–1449 Fig. 4. Qualitative illustration of the dynamics of an extra deterministic neutron inserted at time t = 0 in a critical core. We may estimate the time derivative of ın and integrate to get Fig. 3. Lyapunov stability. The straight line c/c  = n/n is the kernel of the system matrix, so any point that belongs to it is a fixed point. Initial conditions (spots) were taken on a circumference of a radius equal to 5% of the normalised variables. All cases converge to a new fixed point almost horizontally and, in particular, they never leave the shown square. Thus, the fixed point (1,1) is Lyapunov-stable. that is a linear dynamical system with a singular matrix. For a given n0 = / n , the vector  x0 = n0 ˇ n0  (12) is also a fixed point of the system (11). In particular, we call x = [n c  ]T the design point, whose Lyapunov stability we want to address. We define the non-dimensional state vector as  x̃ = n/n  c/c  By virtue of Eqs. (2) and (12), any point located in the first quadrant diagonal of the non-dimensional phase space is also a fixed point. This result is consistent with the fact that the matrix kernel has dimension one. To study the stability of the design point x̃ = [1 1]T we numerically integrate the differential equations using a few different initial conditions (marked as spots in Fig. 3) located on a circumference centred in x with a non-dimensional radius equal to 5 × 10−2 . In all cases, the state vector x(t) evolves to a point located in the c/c  = n/n straight line (marked as crosses in Fig. 3). It is apparent that the design point x = [n c  ]T is Lyapunovstable, as we see that some perturbations ı( ) (in particular the circumference radius) do not make the resulting phase space trajectories leave a given environment of the fixed point (in particular the square enclosing Fig. 3). However, this stability is not asymptotic as the state vector does not go back to the design point but to another different fixed point. Both this effect and the fact that the trajectories are almost horizontal in Fig. 3 may be explained with some reactor physics reasoning (Gho, 2005). When  = 0, summing Eq. (11) gives ṅ + ċ = 0 so the sum n + c is always constant. Suppose we artificially introduce an extra deterministic neutron into a critical reactor core, that lives exactly a time , gets instantly multiplied by a factor of (1 − ˇ) and that a fraction ˇ of a delayed neutron appears exactly a time 1/ after the fission occurred. For times much smaller than 1/, the extra neutron transforms into successive fractions ın(t + ) = (1 − ˇ)ın(t) ıṅ(t  1/) ≈ (1 − ˇ)ın(t) − ın(t) ˇ = − ın(t) ⇒ ın(t) = e−(ˇ/)t   As at this time the precursor has not yet decayed and ın + ıc equals one, for we introduced only one single extra neutron, then ıċ(t  1/) ≈ ˇ ın(t)  ⇒ ıc(t) = 1 − e−(ˇ/)t This is, for t  1/ the extra neutron exponentially disappears and the precursor concentration rises, both with a time constant ˇ/. If /ˇ  1/, then for times t  1/ the extra neutron may have virtually vanished. When t = 1/, the precursor decays, a new neutron appears and the whole process is repeated. Fig. 4 qualitatively shows this reasoning. The result is that the average increase of the neutron population is the mean value of ın(t) ın = 1 1/ 1/ exp − 0 ˇ ·t  dt (13) that for   ˇ/—as usual for common kinetic parameters—gives ın ≈  ˇ (14) The system is still critical but reaches a new equilibrium point given by the new neutron population and the corresponding precursor concentration. It is interesting to note that even though one neutron is inserted, the population rises only by a small fraction /ˇ. The remainder of the extra neutron appears as “stored” into (1 − /ˇ) precursors. This reasoning also explains Fig. 3, where not only was the initial neutron population suddenly changed with respect to the equilibrium, but also the precursor concentration. Moreover, we can easily prove that the trajectories are always straight lines and for common kinetics parameters, essentially horizontal in the non-dimensional phase space. The neutron population derivative may be written using the chain rule as ṅ = ∂n ċ ∂c If we replace Eq. (11), we get ∂n = −1 ∂c (15) thus trajectories in phase space are straight lines. Furthermore, if we transform Eq. (15) into a non-dimensional expression, we obtain: ∂ñ c ˇ =−  =− ≈ 0− n  ∂c̃ which explains the results obtained numerically in Fig. 3. The minus sign and the fact that the diagonal is an equilibrium subspace prove that the design point x is Lyapunov-stable. G.G. Theler, F.J. Bonetto / Nuclear Engineering and Design 240 (2010) 1443–1449 1447 4. Reactivity feedback analysis The multiplication factor keff is a property of the core as a whole, and it depends on a lot of factors (such as the excess initial reactivity, control rods position, fuel burn-up, poisons distribution, and temperatures) in a rather complex way. In particular, it is a function of the thermal power itself so  = (n) and the point kinetics Eq. (3) are not linear anymore. Thus, the stability analysis cannot be addressed with the standard linear theory—at least not directly. Feedback effects may be classified and analysed in a variety of ways—such as those proposed by Lewins (1978) or Duderstadt and Hamilton (1976) for example. One useful approach to compute the change in reactivity due to variations in the conditions of the core is the perturbation formula: ∞ =− V † (r, E) · 0 ∞ V 0  ı a † (r, E) ·  +ı f f · · (r, E) dE dr (r, E) dE dr that may be even further simplified when using a one-group approximation, as the adjoint flux equals the direct flux so Fig. 5. Lumped parameter thermal model of a research reactor. Both the core inlet mass flow ṁcool and temperature Tin are fixed. groups of delayed neutron precursors may be written as =− V 2 (r) · (ı V a 2 (r) · +ı f f ) dr dr However, to correctly compute the reactivity, one must know how the macroscopic cross-sections depend on the thermalhydraulic parameters or the position of the control rods. As stability is a conceptual design concern, this information may not be available at this engineering stage. Thus, we utilise the standard concept of feedback coefficients in this section. The reactivity feedback coefficient due to the effect of the variation of a certain variable x is defined as ˛x = ∂ ∂x (16) so the change in the reactivity introduced by a variation ıx is approximated by ıx ≈ ˛x · ıx  )+˛  ext +˛Tfuel (Tfuel −Tfuel Tcool (Tcool − Tcool ) − ˇ Ṗ =  ċj = ˇj P − j cj  4.1. Study case In this section we show how a very simple time-domain analysis may be used to construct stability maps that could be used to help the core designer during the conceptual engineering of a nuclear reactor. Consider the zero-dimensional model of a research reactor shown in Fig. 5 where coolant enters at a constant rate ṁcool with a fixed temperature Tin . We take a single fuel temperature and a single coolant temperature—equal to the average of the coolant inlet and outlet temperatures—as lumped parameters for representing the thermal state of the core. Treating the reactivity feedback effects as separable with constant coefficients, the core kinetics equations with six 6  i=1 i ci (18) for j = 1 . . . 6 where ext represents the external reactivity perturbations and the coolant temperature coefficient takes into account the change in reactivity due to both coolant temperature and density variations. Using constant parameters for the lumped parameters thermal model, the coolant and fuel temperatures time derivatives may be written as Ṫfuel = 1 [P − (hA)core · (Tfuel − Tcool )] (mc)fuel Ṫcool = (17) Even though the feedback coefficient (16) is well defined, the actual coefficient varies with time and Eq. (17) applies only for small changes in x. Moreover, when there are many variables whose effect one wants to take into account, the usual practise is to sum up every single contribution to the total reactivity as a linear combination, neglecting correlations between these variables and other non-linear effects. Furthermore, reactivity coefficients play a major role in defining the design basis of a nuclear reactor design (Florido et al., 2007) and are extensively used in nuclear regulations and international safety guides (International Atomic Energy Agency, 1999, 2005). P+ (19) 1 [(hA)core (Tfuel − Tcool ) − ṁcool · ccool · 2(Tcool − Tin )] (mc)cool (20) where (mc)fuel and (mc)cool are the total heat capacity of the fuel and the coolant inside the core respectively, ccool is the coolant specific heat capacity and (hA)core is the product of the heat transfer coefficient and the exchange area between the fuel and the coolant. Eqs. (18)–(20) constitute a ninth order dynamical system ẋ = F(x) with ⎡P ⎤ ⎢ c1 ⎢. ⎢. x = ⎢. ⎢ c6 ⎣ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ Tcool Tfuel For a fixed thermal power P  , the other components of the design point x are ˇj  P j  P + 2 · ṁcool · ccool · Tin  Tcool = 2 · ṁcool · ccool  P  + (hA)core · Tcool  Tfuel = (hA)core cj = (21) 1448 G.G. Theler, F.J. Bonetto / Nuclear Engineering and Design 240 (2010) 1443–1449 Fig. 6. Stability maps in the ˛Tcool vs. ˛Tfuel parameter plane. The perturbation in the Lyapunov study is a perturbation of 2 ◦ C in the inlet temperature of the coolant. The perturbation in the BIBO study is a positive 10 pcm external reactivity insertion. Fig. 7. Sensitivity analysis of the Lyapunov stability map with respect to coolant mass flow. (a) 2.00ṁcool , (b) 0.75ṁcool , (c) 0.50ṁcool , (d) 0.25.ṁcool . 4.2. Stability maps When all the reactivity feedback coefficients are negative, it is clear that the reactor is stable. But when there is at least one positive coefficient, a question arises about the overall stability. Given the fact that the model shown in Fig. 5 has only two feedback coefficients, a two-dimensional stability map in the space ˛Tfuel − ˛Tcool may be constructed as follows. We numerically integrate the point kinetics equations starting from a slightly perturbed equilibrium initial condition—either a perturbation in phase space or in the external reactivity—for a definite amount of time. We then repeat this process for different combinations of the two reactivity coefficients. A two-dimensional contour plot of the final thermal power as a function of both coefficients plays the role of a ˛Tfuel − ˛Tcool stability map. In Fig. 6, the resulting stability maps—in both Lyapunov and BIBO senses—for the study case are shown. In the first case, the initial perturbation corresponds to a step increase of 2 ◦ C of the nominal coolant inlet temperature at t = 0. In the BIBO case, the proposed perturbation is a positive external absolute reactivity insertion of 10−4 or 10 pcm. The shape of both maps is essentially the same, showing that the fuel coefficient is the one that dominates the stability of the kinetics equations a as very positive coolant coefficient is needed in order to render the reactor unstable for negative fuel coefficients. The Lyapunov stability map shows a zone labelled “safe-unstable” where for very negative coolant coefficients, the thermal power goes to zero with the proposed perturbation. This behaviour may be defined unstable with respect to the Lyapunov definition, but nevertheless the power deviates to a safe value. Fig. 8. Sensitivity analysis of the Lyapunov stability map with respect to fuel heat capacity. (a) 2.00ṁfuel , (b) 0.75ṁfuel , (c) 0.50ṁfuel , (d) 0.25ṁfuel . G.G. Theler, F.J. Bonetto / Nuclear Engineering and Design 240 (2010) 1443–1449 The fact that the stability border shown does not pass through the origin has to do with the finite time adopted for the numerical integration. Had we provided a longer end time, the unstable zone would have moved to the left tangentially approaching the origin in both cases. As we already proved, there is a difference between Lyapunov stability and BIBO stability when both coefficients are zero. This also explains why the instability zone in Fig. 6(a) is slightly larger than the one in Fig. 6(b). 4.3. Sensitivity analysis The proposed method of constructing these conceptual stability maps may be used to evaluate the impact of modifying a certain design parameter. For example, Fig. 7 shows how the stability limit changes as the coolant mass flow is reduced to different fractions of the design value ṁcool and Fig. 8 the influence of the fuel total calorific capacity (mc)fuel . The coolant mass flow modifies the stability border and for high mass rates, the fuel temperature reactivity coefficient dominates. On the other hand, the fuel calorific heat capacity does not appreciably affect the stability zone in the proposed model. These results are applicable only to the lumped parameter model and constant feedback coefficients proposed in Section 4.1, but even with the simple analysis shown, the designer may reach some interesting conclusions. First, that the stability topography does not depend heavily on the fuel calorific capacity. On the other hand, the coolant flow rate does change the stability border, but apparently not the safe-unstable zone. For high cooling flows, the coolant temperature is almost constant and the effect on the neutron kinetics is negligible. Therefore, if for some reason a core design has a positive coolant temperature reactivity coefficient, the reactor may be rendered stable by having a good cooling system thus driving the stability essentially with the usually negative fuel temperature coefficient. 5. Conclusions The point reactor kinetics equations are a useful simplification of the time-dependant neutronic Boltzmann problem that allow studying the stability of a reactor core using dynamical systems theory. From a strictly mathematical point of view, we proved that the point reactor equations—in the absence of an external neutron source and without any feedback effects—are not BIBO-stable but are, nevertheless, Lyapunov-stable. Moreover, we were able to show that a simple proportional controller acting on the external reactivity makes the equations BIBO-stable. However, modifying directly the reactivity may not be possible and thus delays and other non-linearities may arise, changing the base problem. The result of unconditional linear stability applies only to the situation depicted in Fig. 1 as a pure theoretical analysis. It should be 1449 remarked again that all the studies and conclusions apply directly only to the mathematical models presented. The maps shown in Figs. 6–8 are just the results of the numerical integration of Eqs. (18), (20) and (19). The applicability of the methods and results to a real nuclear reactor is not straightforward. However, the basic physics of the feedback effects between the thermalhydraulics and the neutronics are retained. This simple method—namely solving a dynamical system for different parameters and analysing the results in the time domain—may be extended to handle more complex cases. On the one hand, existing codes may be used to solve the governing equations thus avoiding a deep mathematical development, at least in early stages of the conceptual engineering. On the other hand, as each system to be solved is independent from the others, the process can be easily parallelised in modern multi-core computers and preliminary results may be obtained quickly. This kind of non-linear stability analysis is a powerful aid for nuclear reactor engineering, and the basic physics underlying the neutronic core design can be better understood by means of these tools. References Doutta, G., Doshi, J., 2009. Nonlinear analysis of nuclear coupled density wave instability in time domain for a boiling water reactor core undergoing corewide and regional modes of oscillations. Progress in Nuclear Energy 51 (8). Duderstadt, J.J., Hamilton, L.J., 1976. Nuclear Reactor Analysis. Wiley, New York. Florido, P.C., Brasnarof, D., Lestani, H., 2007. Base de diseño del coeficiente de vacío del combustible CARA. In: Actas de la XXXIV Reunión Anual de la Asociación Argentina de Tecnología Nuclear. Garea, V., Drew, D., Lahey, R.T., 1999. A moving-boundary nodal model for the analysis of the stability of boiling channels. International Journal of Heat and Mass Transfer 42 (19). Gho, C., 2005. Reactor Physics lecture notes, Undergraduate Nuclear Engineering Course at Instituto Balseiro, Argentina. Henry, A.F., 1975. Nuclear Reactor Analysis. MIT, Cambridge. International Atomic Energy Agency, 1999. Basic Safety Principles for Nuclear Power Plants 1. INSAG 12. IAEA. International Atomic Energy Agency, 2005. Design of the Reactor Core for Nuclear Power Plants. Safety Guide No. NS-G-1.12. IAEA. Lahey, R.T., Podowski, M., 1989. On the analysis of various instabilities in two-phase flows. Multiphase Science and Technology, 4. Lewins, J., 1978. Nuclear Reactor Kinetics and Control. Pergamon Press. Peng, S., Becker, M., Lahey, R.T., Podowski, M., 1984. NUFREQ-NP: a computer code for the stability analysis of boiling water nuclear reactors. Nuclear Science and Engineering 88 (3). Peng, S., Podowski, M., Lahey, R.T., 1986. BWR linear stability analysis. Nuclear Engineering and Design 93 (1). Salah, A.B., D’Auria, F., 2007. Use of coupled code technique for best estimate safety analysis of nuclear power plants. Progress in Nuclear Energy 49 (1). Salah, A.B., D’Auria, F., 2009. Application of coupled code method for the analysis of VVER-1000 coolant transient benchmark. Progress in Nuclear Energy 51 (1). Solari, H.G., Natiello, M.A., Mindlin, G.B., 1996. Nonlinear Dynamics. Institute of Physics Publishing, Bristol/Philadelphia. Stacey, W.M., 2001. Nuclear Reactor Physics. Wiley–VCH. Theler, G., 2008. Non-linear analysis of the coupled thermalhydraulicneutronic problem. Master’s Thesis, Instituto Balseiro, Argentina. Zanocco, P., Giménez, M., Delmastro, D., 2003. Safety design maps: an early evaluation of safety to support reactor design. Nuclear Engineering and Design 225 (2–3).