Dynamics of Fluid Fuel Reactors in the Presence of Periodic Perturbations

The appearance of perturbations characterized by a periodic time behavior in fluid fuel reactors is connected to the possible precipitation of fissile compounds which are moved within the primary circuit by the fuel motion. In this paper, the time-dependent response of a critical fluid fuel system to periodic perturbations is analyzed, solving the full neutronic model and comparing the results with approximate methods, such as point kinetics. A fundamental eigenvalue of the problem is defined, characterizing the trend of divergence of the power. Parametric studies on the reactivity insertion, the fuel velocity, and the recirculation time are performed, evidencing the sensitivity of the eigenvalue on typical design parameters. Nonlinear calculations in the presence of a negative feedback term are then performed, in order to assess the possibility to control a fluid fuel system when periodic reactivity perturbations are involved.


INTRODUCTION
Circulating fluid fuel reactors have been studied with alternate fortune during the last forty years. The interest in the development of this reactor design lies in the promising perspectives concerning the management of the fuel cycle and in the possibility to perform the transmutation of long-life radioactive wastes. For these reasons, the molten salt reactor concept meets Generation IV requirements in terms of sustainability [1].
In the years 1965-1969, an extensive experimental activity on fluid fuel reactors was performed at Oak Ridge Laboratories [2,3]. The Molten Salt Reactor Experiment (MSRE) proved the feasibility of the design and allowed to study the peculiar physical phenomena characterizing this kind of systems. Experience in the field of fuel performances, material irradiation and corrosion, operation and maintenance was also acquired.
The experimental results and conclusions drawn from the MSRE program have served as a starting point for the research activities developed in recent years in different countries in Europe. The MOST Project [4], financed by the European Community within the V Framework Program, has taken profit of the MSRE results performing an assessment of the state of the art of the technology and carrying out a benchmark activity on the available computational tools for fluid fuel systems, tested against experimental MSRE data [5]. A Coordinated Research Project supported by IAEA, "Studies of Advanced Reactor Technology Options for Effective Incineration of Radioactive Waste" [6], has developed a research task devoted to the safety assessment of a novel design of molten salt reactor, the MOSART concept [7], performing dynamic calculations [8].
The physics of fluid fuel systems is characterized by some peculiar phenomena, connected to the movement of the multiplying material through the core and in the primary circuit. The motion of the fuel affects the delayed neutron precursors, which are dragged within the system by the fluid flow; as a consequence, the emission of delayed neutrons takes place in different spatial positions with respect to the corresponding prompt emissions. This effect results in a spatial redistribution of delayed neutron precursors, depending on their lifetime, and a general reduction of their role in the fission process, since a part of them is swept outside the reactor core and decays in the primary circuit.
The reduction of importance of delayed emissions associated to the motion of the fluid fuel affects both the static and dynamic behavior of the reactor. The critical condition depends on the velocity regime and the effective fraction of delayed neutrons is reduced, with obvious effects on the time-dependent response of the system.
The physics of such systems also allows the appearance of peculiar perturbations of the reactor material composition which constitutes an interesting aspect of neutron dynamics to be studied. In this paper, the problem of the possible appearance of perturbations with a periodic behavior in time is approached: this aspect is strongly connected to some physical phenomena occurring in fluid fuel reactors.
The operation of a molten salt reactor and the consequent modification of the fuel composition due to burnup can cause the precipitation of fissile compounds. Since the isotopic vector of the fissile components is modified by fissions, the addition of fissile material to keep reactivity can lead to the solubility limit for some fissile nuclides, which as a result precipitate. The consequent lump of fissile material is driven through the system by the fuel motion and represents a localized perturbation of the multiplicativity of the medium with a periodic behavior in time. Analogously, the formation of helium gas bubbles within the reactor during operation constitutes an additional, localized perturbation of the material characteristics of the core which is swept throughout the whole system by the fuel velocity. The presence of gas bubbles was experimentally detected during the operation of MSRE [2] and represents an issue to be solved for the assessment of the reliability of the technology. The amount of reactivity associated to these phenomena can be rather relevant, amounting to hundreds of pcm [9].
The analysis of the time-dependent response of a molten salt reactor to this kind of perturbations can be treated by taking profit of previous studies on the peculiar phenomena associated to periodic reactivity insertions. In previous works [10][11][12] the problem of a periodic reactivity insertion has been approached by means of an analytical solution of the point kinetic model for both critical and subcritical systems. More recently, a preliminary analysis on the precipitation of fissile lumps and the formation of gas bubble in fluid fuel systems has been carried out [13].
In this paper, a comprehensive study of the physical phenomena associated to periodic perturbations in fluid fuel systems is performed. The time-dependent behavior of the system in response to a fissile lump moving within the primary circuit is simulated with point kinetics and solving the full space-time neutronic model, in order to evidence the role of spatial phenomena connected to the movement of a localized perturbation in the reactor. Results for the pure neutronic model are presented, together with calculations in the presence of thermal feedback, in order to determine the stability condition of the reactor.

PHYSICAL-MATHEMATICAL MODEL FOR THE NEUTRONICS OF FLUID FUEL REACTORS
The study of the neutronics of fluid fuel reactors requires the proper modification of the neutron and precursor balance equations, in order to take into account the presence of a velocity field within the reactor core. Prompt fission emissions are not affected, due to their different time scale, while a con-vective term is introduced in the equation for precursors. The general model for neutrons and precursors can be written in terms of neutron density n and delayed emissivities E i as [14] ∂n(r, E, Ω, t) ∂t R families of delayed neutron precursors are considered and the delayed emissivity is related to the precursor density C i (r, t) by the expression In model (1), the leakage L, prompt multiplication M p and delayed multiplication M i operators are introduced, while the velocity field, represented by the vector u(r, t), is taken as an input function. This system of equations requires initial and boundary conditions for both the neutron density and the delayed emissivities. Boundary conditions for delayed emissivities need to properly represent the correlation between the exiting flow of precursors at the outlet of the reactor core and the undecayed fraction of precursors reentering the system after a certain time interval. The general mathematical formulation of this condition can be written as where the recirculation time in the external circuit τ is introduced. Condition (3) equates the incoming flow of precursor into the reactor core at time t, left-hand side of the equation, to the undecayed fraction of the outgoing flow of precursor evaluated at a preceding instant in time t − τ, related to the transit time in the recirculation loop. A geometrical redistribution function F is also defined, expressing the probability for a precursor exiting at point r to re-enter the core at point r [15].
The solution of the full model (1) in the case of the precipitation of a fissile compound can be modeled as a localized perturbation of the multiplicativity of the medium moving at velocity u. This solution approach allows to correctly describe all the spatial and spectral effects associated to the presence of the moving perturbation, but it obviously requires a large computational effort.
An alternative procedure to simulate this kind of transients implies the adoption of a point-like kinetic model.
Since the physical-mathematical model has been modified for molten salt reactors, also the corresponding approximate models, such as point kinetics, need to be extended and adapted. A consistent formulation of point kinetics can be obtained by making use of the standard factorizationprojection technique, and a generalized formulation of the time-dependent amplitude equations is obtained [14] dP dt where P and G i are the amplitude functions for the neutron distribution and for the ith precursor family, respectively. The kinetic parameters in system (4) have a different definition from the standard formulation and depend on the shape of neutrons and precursors, since both unknowns are factorized and projected in the procedure. Additional terms, such as ζ i and η i , appear as a consequence of the fuel motion, allowing this reformulated point model to consistently describe the dynamics of a fluid fuel reactor, although all spatial and spectral effects are neglected. In the point kinetic framework the simulation of the transit of a localized perturbation affects the value of the reactivity ρ and the value of the effective delayed neutron fraction β when perturbations of multiplicativity are concerned.
The study of periodic perturbations in molten salt reactors is here performed solving the full model (1) in cylindrical (r −z) geometry and comparing it to point kinetic results. The perturbation is moved within the system along the axial direction from the core inlet to the core outlet in the time interval τ c , while its residence time in the external circuit is τ l . The time-dependence of the reactivity ρ during τ c can be qualitatively assimilated to the function sin 2 (πt/τ c ), according to elementary perturbation theory in diffusion approximation [16], while in the interval τ l the kinetic parameters of the reactor are left unperturbed. Also the effective delayed neutron fraction β experiences a perturbation with a periodic behavior in time.
The analysis is focused on the stability characteristics of the power oscillations appearing as a consequence of the periodic reactivity introduction and on the role of spatial effects, which are neglected when point-like calculations are performed. The model (1) is solved by numerical inversion of the full operator, while the point kinetic solution is produced adopting a semianalytical approach [12], trying to evidence the asymptotic trend in the power evolution. In a previous work [11] it has been proved that a periodic perturbation with zero mean value results in a divergence of the power in a critical system, due to the accumulation of precursors. In the case of the precipitation of a fissile lump a reactivity oscillation with a positive mean value is involved, and consequently the power is expected to be diverging.
At last, results in the presence of thermal feedback are presented. The point kinetic model is coupled to a zerodimensional thermal model for the nuclear fuel, performing a thermal balance over the full core. The fuel temperature is updated on the basis of the power production, the feedback on neutronics is performed by means of an integral parameter α, such as and the total reactivity ρ tot is inserted into (4). The implementation of a thermal module allows to perform a preliminary analysis of the effect of thermal feedback on the stability of the reactor in the presence of these perturbations, evidencing the influence of various physical parameters.

SEMI-ANALYTICAL SOLUTION OF THE POINT KINETIC MODEL
The interest in the study of periodic perturbations relies in the possibility to evidence a peculiar behavior connected to the different time scales characterizing the evolution of the neutron and precursor populations. Moreover, in some specific and physically meaningful cases, like the problem of a square wave reactivity insertion, the point kinetic problem can be treated fully analytically evidencing the asymptotic trend of evolution of the point-like reactor, in both critical and subcritical configurations [11,12].
In the case of a circulating lump or a gas bubble the associated reactivity has a time behavior which cannot allow a totally analytical approach. The localized modifications of the cross sections in the core during the time interval τ c can be approximated by a sin 2 (πt/τ c ) function, while during τ l the inserted reactivity vanishes. In the calculation performed the movement of the perturbation has been simulated on the discrete spatial mesh in the cylindrical (r − z) reactor considered, and the supposition of a"sin 2 " behavior based on perturbation theory has been proved. Analogously, the effect on the value of β has been evaluated.
In the adopted semianalytical approach the full period of the oscillation T = τ c + τ l is discretized in N time steps Δt, retaining the kinetic parameters constant over each interval. The differential problem over the ith time interval in a critical source-free system can be written as where K i is the point kinetic matrix with kinetic parameters associated to the ith time step. In the following of this section only one family of delayed neutron precursors is considered to simplify notation, but the procedure is applicable to the general case of R families. System (6) is discretized on the same time step Δt using an implicit Euler algorithm for the time derivative. During the first time step, the system of algebraic equations which is obtained can be written as where |X(Δt) is the state of the system after the first time step, |X(0) is the initial condition of the system, and the matrix ϑ 1 is defined as The solution of problem (7) can be obtained by analytical tools preliminarily evaluating the real eigenvalues ( (1) 1 and (2) 1 ) and eigenvectors (|Γ (1)   1 and |Γ (2) 1 ) of the matrix ϑ 1 and of its adjoint ϑ † 1 , Γ (1) 1 | and Γ (2) 1 |, solution of the following homogeneous algebraic problems: Due to equality (8), these eigenvectors turn out to be also the eigenvectors of the point kinetic matrix K 1 ; the established relation for the eigenvalues is where ω (1) 1 and ω (2) 1 are the eigenvalues of the point kinetic matrix K 1 . The solution after the first time step can thus be expressed as The solution of the problem over the whole period T = NΔt is obtained by applying formula (11) repetitively over all time steps; the heavy formulation of the final solution can be simplified by introducing some assumptions related to the well-known characteristics of the point kinetic eigenvalues The second eigenvalue of the point kinetic matrix ω (2) j can be approximated as (ρ − β)/ p and is largely negative; the corresponding eigenvalue of matrix ϑ i can be reduced below unity by satisfying the assumption In the cases considered for this analysis ω (2) j is of the order of −10 1 s −1 and from condition (12) this circumstance implies Δt < 10 −1 s. Satisfying condition (12) allows to introduce the following simplification: Other two approximations are introduced: requiring a "pseudo-orthogonality" of the eigenvectors associated to subsequent time steps which is exactly verified for Δt approaching zero. Consequently, the choice of a sufficiently small Δt justifies all the assumptions (13) and (14) and allows to write the solution over the period T as In (15), ξ is the fundamental eigenvalue, while |Ψ is the fundamental state of the system, defined as Applying recursively (15), it is possible to describe the asymptotic evolution of the system for each transit of the fissile lump: where the exponent m indicates the number of cycles.
Observing (17), it appears that the asymptotic evolution of a multiplying system excited through a periodic perturbation is governed by a fundamental eigenvalue that operates as an amplification factor for each cycle. The real eigenstate is obviously not the one calculated by means of point kinetics, but it can be evaluated in general by performing the inversion of the full space-angle-energy-time model.
When performing the inversion of the full model (1) an analytical approach cannot be followed, however the fundamental eigenvalue for the asymptotic evolution of the system can be defined as In real calculations, the evaluation of formula (18) is performed after few cycles owing to the large computational effort associated to the inversion of the full model, trying to evidence the appearance of a stable value of ξ.
The assumption of the existence of a "fundamental state" for the full model (1), appearing when the asymptotic regime is established, can be understood on the basis of well-known physical properties of linear systems in the presence of a periodic input. After an initial transient a dynamic oscillating equilibrium condition is reached for both the populations of neutrons and precursors over a full period; as a consequence, the asymptotic evolution can be described by an integral parameter as in formula (18).

NUMERICAL RESULTS
Calculations are performed using the numerical code DY-NAMOSS [15], solving the multigroup diffusion equations for a fluid fuel system in cylindrical (r − z) geometry and   imposing an axial velocity field without dependence on the radial variable. Two different reactor configurations are considered: the first one assumes geometrical dimensions and material data typical of MSRE with two energy groups and one family of precursors [3]; the second one reproduces the main characteristics of the MOSART design [7] with four energy groups and six precursor families. In both systems, the reactor core has homogenized cross sections and the MOSART-like geometry includes a radial reflector. In Table 1 the main data for the two systems are summarized.
The first set of results concerns the MSRE system in the presence of a fissile lump moving through the system and a parametric analysis is performed on physically meaningful quantities influencing the power evolution. In Figure 1, the flux perturbation associated to the movement of the fissile solid component is shown at different instants during the circulation period in the core τ c . The maximum reactivity associated to this perturbation is about 250 pcm . The appearance of spatial and spectral effects is clearly visible: the increased multiplicativity in this localized region causes a larger absorption of thermal neutrons with a consequent additional production of fast neutrons by fission.
The power evolution following the transit of the fissile lump is then simulated using the point-like approach described in the previous section and it is compared to the results obtained by inversion of the full space-time operator, in order to compare the capability to predict the trend of the power divergence by simplified models, with respect to more accurate ones. In this case the reactivity associated to the localized perturbation amounts to 25 pcm at its maximum. In  The first half period is characterized by an insertion of reactivity into the system, with a consequent increase of the power, which is reproduced by both the dynamic models adopted. During the second part of the period the non-equilibrium condition between neutrons and precursors causes a decrease of the power which is predicted with large differences between the two options. The importance of spatial effects can also be clearly evidenced by observing the modification of the neutron flux along the axial coordinate during the transit of the lump. The spatial modifications introduced by the perturbation on the fission cross section are neglected by the point model, thus giving large discrepancies in the prediction of the power evolution. Figure 2(a) also reports the asymptotic evolution of the power for the full and point model, according to formulae (17) and (18). The agreement of the curves for point kinetics is very good, while in the full case the asymptotic establishment of a "fundamental mode" is not reached during the first oscillation periods, thus giving larger discrepancies. In Table 2, the relation between maximum reactivity insertion and value of the fundamental eigenvalue ξ, evaluated through model (17), is reported. The power starts rising very quickly, since for each cycle it is amplified by a factor ξ. The last example regards a super-prompt-critical configuration, since the maximum reactivity insertion is larger than the value of β.
A square-wave problem producing the same results in terms of fundamental eigenvalue and eigenstate is now considered. Such a problem allows a fully analytical treatment. In Figure 3 the previous graphs are compared to the response for the corresponding square-wave problem. The half-period of the oscillation with a positive insertion of reactivity for both ρ and β corresponds to the time of transit of the lump in the core, while the amplitude is slightly larger (few percents) than the mean value of ρ(t) and β(t) on the interval in which a positive reactivity is inserted. The behavior of the state of the system is quite different during the oscillation period, but the reactor is exactly in the same conditions at the end of each oscillation.
Some additional parametric studies on important quantities are also performed, solving the time-dependent problem with the point kinetic model in order to easily evidence the main integral parameters coming into play. At first, the influence of the recirculation time in the external circuit τ l is analyzed. The geometry of the reactor core is left unchanged in terms of composition and velocity field, while the value of τ l is modified. The critical condition is determined and then the transit of the fissile lump is simulated by perturbing locally the multiplicativity, keeping the relative variation with respect to the critical value constant in all calculations.
In Figure 4, the behaviors of the main parameters of the problem as functions of τ l are reported. The influence on the reactivity insertion is very small, as could be expected since the relative perturbation on cross sections is kept constant, while the increase of the recirculation time causes a corresponding increase of the fundamental eigenvalue. This results in a faster divergence of the power, due to the fact that, with an increased τ l , the role of delayed neutrons is reduced and the reactivity insertion in dollars turns out to be larger. This effect is confirmed by observing Figures 4(c) and 4(d reporting the mean value of β during the transit time in the core τ c (β + ) and the value experienced during the recirculation time τ l (β − ). It must be noted also that the effect on these two values of the effective delayed neutron fraction is comparable, meaning that the change of τ l does not affect the shape in time of the oscillation of β. As a further comment, the appearance of a saturation value for all the integral parameters can be observed, as the increased recirculation time approaches the asymptotic condition in which no delayed neutron precursors re-enter the core [15]. An analogous set of calculations is performed changing parametrically the fuel velocity. The geometry of the system is left unchanged, thus an increase of the fuel velocity implies a reduction of both transit times τ c and τ l . In Figure 5 the results obtained are presented. The value of the reactivity is again basically unperturbed, while some peculiarities can be observed concerning the values of ξ and β. As the fuel velocity increases the effective delayed neutron fraction is reduced, as it is already observed in the parametric study on τ l and as could be expected. On the other side, this reduction of β does not imply an increase of the fundamental eigenvalue, which on its turn decreases. This effect can be explained by the fact that, with the increased velocity, the time spent by the lump on each cycle is reduced as well as its effect in increasing the power level and accumulating precursors. The appearing of a saturation value for all the parameters as the velocity increases and approach to the "infinite velocity" condition can also be observed [15].
At last, results in the presence of a simplified zerodimensional thermal model are presented. The aim of this analysis is to enlighten the possibility to stabilize the reactor behavior and prevent the power divergence. The feedback model adopted introduces an additional reactivity term in the point kinetic equations (4) connected to the modification of the fuel temperature, as in formula (5). The reactor design considered is of MOSART-type: the starting reactor power is 2400 MW, the core inlet temperature T in is 600 • C and the average fuel temperature T fuel,0 is 660 • C. Parametric studies on the value of the feedback coefficient α and on the recirculation time τ l are performed. The possibility to control the reactor with a negative feedback effect is connected to the saturation of the precursor concentration in the system. In Figure 6 the evolution of the precursor amplitude for the first family, G 1 , is plotted adopting different values of the feedback coefficient α. Reference values for the velocity and the recirculation time are adopted, see Table 1. The stabilization of the precursor concentration is clearly visible, implying the stabilization of the system power which keeps oscillating with a constant mean value. The increase of the role of thermal feedback obviously reduces the amplitude of the excursion from the initial value to the asymptotic one. In Figure 7, the response of three different precursor families is plotted during the first oscillation period, showing different behaviors connected to their decay constants. The role of thermal feedback is also evidenced.
The effect of thermal feedback on the system power and on the amplitude for all six precursor families is shown in Figure 8, where the maximum values of these quantities over each oscillation are plotted as functions of the number of periods m. The comparison of the fastest decaying family G 6 to the one with the longest mean life G 1 shows that the number of periods needed in order to reach a stable asymptotic level is higher for the families with smaller decay constants. Consequently, the number of oscillations required to establish the asymptotic equilibrium condition is regulated mainly by the first family.   A parametric study on the recirculation time τ l is then performed, as shown in Figure 9, where the maximum and minimum value of the stable oscillation for power and precursors is plotted as function of the recirculation time. An increased value of τ l results in larger oscillations of both the precursor concentration and the power, since an increased recirculation time reduces the value of the effective delayed neutron fraction. This phenomenon affects mainly the maximum value reached during the oscillation; the importance of feedback effects in order to reduce the amplitude of the power oscillations is clearly visible.

CONCLUSIONS
The problem of the precipitation of solid fissile salts in circulating fuel reactors has been thoroughly analyzed, enlightening its effects on the dynamic behavior of the system. The power oscillations induced by a localized perturbation of multiplicativity moving through the reactor have been simulated by means of the solution of the full space-time model and using a simplified point-like approach, evidencing the appearance of an asymptotically diverging fundamental state. The differences in the evaluation of the timedependent response of the system are discussed and show the importance of full space dynamic tools for a realistic prediction of this type of transients in molten salt reactors. The results of a parametric analysis on important physical parameters such as the fuel velocity and recirculation time have been discussed. Calculations in the presence of a thermal feedback model have been carried out, enlightening the stabilizing effect of negative feedback and the establishment, after an initial transient, of a stable oscillation around a constant power level. These numerical evaluations constitute a preliminary assessment on the safety issues connected to periodic reactivity perturbations in a reactor operating at full power.
Further developments on this research subject include the implementation of more detailed thermal models in order to take into account spatial effects, as well as improvements in the neutronic description. The introduction of more complex fluid-dynamic models and the development of a fully coupled neutronic-thermal-fluid-dynamic module is envisaged, in order to be able to perform a comprehensive analysis of the reactor design.