The effect of a changing fuel solution composition on a transient in a fissile solution

Abstract This paper presents an extension to a point kinetics model of fissile solution undergoing a transient through the development and addition of correlations which describe neutronics and thermal parameters and physical models. These correlations allow relevant parameters to be modelled as a function of time as the composition of the solution changes over time due to the addition of material and the evaporation of water from the surface of the solution. This allows the simulation of two scenarios. In the first scenario a critical system eventually becomes subcritical through under-moderation as its water content evaporates. In the second scenario an under-moderated system becomes critical as water is added before becoming subcritical as it becomes over-moderated. The models and correlations used in this paper are relatively idealised and are limited to a particular geometry and fissile solution composition. However, the results produced appear physically plausible and demonstrate that simulation of these processes are important to the long term development of transients in fissile solutions and provide a qualitative indication of the types of behaviour that may result in such situations.


Introduction
A fissile solution is an aqueous solution formed of a fissile solute (such as uranyl nitrate) dissolved in water and, potentially, an acid component (such as nitric acid) to increase the solubility of the main solute. Fissile solutions may be used in Aqueous Homogeneous Reactors (AHRs) or as part of fuel fabrication or waste management processes. In the case of AHRs, criticality and a non-zero power is a desirable quality of the system as it allows the functioning of the reactor. In the case of fuel fabrication and waste storage, criticality is to be avoided. However, there have been several accidents involving such solutions such as the Y12 accident (Patton et al., 1958) and the Tokaimura accident (Komura et al., 2000). For either the safe operation of an AHR or the prediction of an accident scenario in a fissile solution it is important to be able to simulate the behaviour of a transient within a fissile solution. Point kinetics codes are commonly used for this purpose (Mather et al., 2002;Mitake et al., 2003;Cooling et al., 2014b) but higher dimensional models which couple neutronics transport and Computational Fluid Dynamics (CFD) have also been produced (Buchan et al., 2013). The purpose of this work is to develop an improved point kinetics model that will track the effects of changing composition of a fissile solution during a criticality accident. This is particularly relevant for accidents such as the Y12 accident (Patton et al., 1958;Zamacinski et al., 2014) where the addition of water caused the solution to become first critical and then subcritical again. The model is very simple and is based upon the models found in , Cooling et al. (2014a) and Zamacinski et al. (2014). The additions to the models presented in those works will concern themselves with the simulation of changing composition due to the addition of material and the evaporation of water and the production of empirical correlations describing key neutronics parameters as a function of the state of the system including the composition of the solution. Although Basoglu et al. (1998) has examined evaporation from the solution surface before, it is the authors' belief that this work represents the first attempt to use a point kinetics model to dynamically simulate the effects of a changing composition caused by dilution or evaporation on a transient as it progresses. It is assumed that few enough fissions will occur during the simulated transients that burnup will not cause the composition of the system to vary significantly or for a significant number of fission products to be created. As a result, simulation of the effects of burnup is neglected. The resulting model is applied to two cases in Section 3. In the first, the system begins with excess reactivity and is initially over-moderated. It is eventually shut down by the evaporation of water from the solution which leads to a reduction in moderation to the point where the system becomes subcritical, causing the fission rate to drop to near zero. In the second, water is added to an initially under-moderated and subcritical system in order to cause the system to become critical and an excursion to occur before the added water eventually leads to the system becoming over-moderated and subcritical one more, halting the reaction.

Model
The model assumes a simple cylinder of solution of radius 0.32m and a surface height that is free to move dependent on the total mass and density of the solution. The solution contains water, nitric acid and uranyl nitrate with an enrichment of 20%. As a result the elements present are limited to hydrogen, oxygen, nitrogen, uranium-235 and uranium-238. The neutronics variables of the reactor are described as point values, the temperature of the solution is assumed homogeneous and only the total void volumes are tracked. As a result no parameter discussed has any spatial variation. The power of the system and the concentration of the six groups of delayed neutron precursors are governed by the standard point kinetics equations. The radiolytic gas in the system is modelled to be formed immediately in stoichiometric proportions. This simplification is consistent with the physical case that the system is already fully saturated with radiolytic gas, meaning a more complex model of dissolved gas, such as that found in Zamacinski et al. (2014) is unnecessary. Steam bubbles within the solution are produced at a rate proportionate to the super-heat of the system. This occurs after the creation of radiolytic gas as radiolytic gas is produced in a transient before the solution has warmed sufficiently for boiling to occur which means the radiolytic gas bubbles may act as nucleation sites for the boiling. Both radiolytic gas and steam leave the system as the gas exits the top of the solution as in Zamacinski et al. (2014).  found the characteristic upward velocity for radiolytic gas is approximately 4.35cm/s and this will be used as the upward velocity of the gases in this model. The temperature of the solution is increased by the energy released by fission and reduced by conduction through the sides of the vessel, the addition of new material, the creation of steam and evaporation from the surface of the solution. The resulting expression for the rate of change of temperature is given in Equation (1): where T S (t) is the temperature of the solution, P (t) is the fission power,Ė B (t) is the rate at which energy is removed from the solution for the production of steam,Ė side (t) is the rate of heat loss through the sides of the container to the environment (which is considered to have a constant temperature of 300K),ṁ a (t) is the mass addition rate for material added to the system, c a and T a are the specific heat capacity and temperature of the added material,ṁ e (t) is the rate at which mass is removed from the solution through the evaporation of water at the top surface of the solution, L s is the latent heat of evaporation of water to steam and m S (t) and c S are the mass and specific heat capacity of the solution with the latter assumed constant. Many of the terms in Equation (1) are direct analogues of those used in Zamacinski et al. (2014) but the term relating to the evaporation from the surface is a new addition and is discussed in more detail in Section 2.1. This is the only modification amde to this equation compared to the equivalent presented in Zamacinski et al. (2014). In the interests of creating a simple, abstract model, no assumption is made regarding the environment external to the fuel solution. Instead, it is assumed that the exterior is held at a constant temperature of 300K and the heat transfer coefficient through both the sides and base to this temperature is 100W/K/m 2 .

Evaporation
The model includes several equations meant to model the effects of evaporation of water from the surface of the solution which, in contrast to boiling within the solution, will occur even when the solution is below its saturation temperature. The presence of salts in a solution will reduce the rate of evaporation compared to pure water. However, little data is readily available on the way that uranyl nitrate solute affects the evaporation rate so the model makes the approximation that the evaporation at the surface occurs as if the solution was pure water. This is clearly an assumption which reduces the accuracy of the model and an ambition for the future would be to update the evaporation rate to reflect the effect of the dissolved uranyl nitrate. To evaluate the rate at which mass is removed from the solution surface through the evaporation of waterṁ e a correlation found in Bansal and Xie (1998) is employed (with the assumption that air flow over the surface is negligible): whereṁ e is the rate at which water evaporates from the surface in units of kg/s, r S is the radius of the circular surface in m, p v is the vapour pressure of the liquid in kPa, and p wa is the partial pressure of the water in the air above the surface in kPa. Equation (3) notes the Antoine Equation and is used to find the vapour pressure of the solution p v : where p v is the vapour pressure in kPa, T S (t) is the temperature in Celsius and A, B, and C are constants specific to the evaporating . In this model, A, B, and C depend on the ambient temperature. If T S (t) < 100 o C, A = 8.07131, B = 1730.63, and C = 233.426. Otherwise, A = 8.14019, B = 1810.94, and C = 244.485. For the purposes of this study we will assume an ambient temperature of 300K and an ambient humidity of 50% for the purposes of calculating p wa which is done using Equation (3) and multiplying the resulting value of p v by the humidity resulting in a value for p wa of 1.785kPa.

Solution Density
The density of the solution is used to determine the height of the solution surface. Zamacinski et al. (2014) derived a correlation for the density of uranyl nitrate of a specific concentration of nitric acid. Through the use of experimental data relating to the density of uranyl nitrate found in UKAEA (1975) this correlation has been augmented to include the effect of varying nitric acid concentrations in Equation (4): where ρ S (t) is the density of the solution, T S (t) is the temperature of the solution in K and U S (t) is the uranium mass fraction of the solution and N S,acid is the mass fraction of nitrogen contained in nitric acid (as opposed to the uranyl nitrate). Comparison of the results of this correlation with data found in UKAEA (1975) found agreement to within 5% in all cases across a wide range of conditions and better agreement (∼ 1%) in the majority of cases.

Neutronics Correlations
The wide range of possible states of the system in terms of composition, temperature and geometry led to the construction of correlations for the k eff , generation time Λ, the delayed neutron fractions for the six groups β i and the delayed neutron precursor decay rates for each of the six groups λ i . These correlations were formulated via the construction of MCNP models of the system in a number of different configurations that varied the mass, nitric acid concentration, uranium concentration, voidage and temperature (and hence the solution density and height of the solution surface). The results of these MCNP calculations are found in Appendix A. These correlations may be evaluated in a quasi-static fashion in order to evaluate the neutronics parameters as evaporation, addition of material, heating and so on move the system around the parameter space considered as a simulation progresses. The correlations presented in this section present the types of behaviour one might expect from the system although it would be desirable for future work to include additional scenarios to further improve the correlations. The correlations are only valid for the particular system presented in this paper with the facts that the system is a cylinder with a particular radius, that the enrichment of the uranium is 20% and that there is no reflector (or any other surrounding material) being the primary factors that restricts the applicability of these correlations to the scenario studied here. A more general approach would require dynamically solving the neutron transport equation or some approximation to it for the given arrangement of the system, although this would require a substantially more complex model. The first empirical correlation which is fitted to the data presented in Appendix A is Equation (5) which describes the k eff of the system: where m S is the mass of the solution in kg, M HNO 3 is the concentration of HNO 3 in moles per litre, V F S is the void fraction of the solution/void mixture T S is the solution temperature in K and H U is the ratio of moles of hydrogen to moles of uranium. This expression is an empirical correlation developed here to represent the data in Appendix A and so all terms do not have an obvious physical analogue. However, it can be seen that the k eff increases with mass and tends to an asymptotic value as mass increases. Increasing the concentration of nitric acid slightly decreases the reactivity but the effect is less than that of other parameters for practical values. Increasing the voidage or solution temperature decreases k eff whilst the relationship between k eff and the hydrogen to uranium ratio is more complex. For the range of values studied in this paper, k eff forms a peak at a ratio of around 72 (corresponding to the optimally moderated state) and decreases at a modest pace on either side of this peak as the ratio changes. The generation time is described by the correlation given in Equation (6): where Λ(t) is the generation time in µs and all other variables have the same meaning and units as in Equation (5). This correlation produces generation times which, at worst, differ by around 10% from the MCNP results but are generally accurate to within 5%. This expression is independent of the total mass of the solution as simply extending the extent of the solution will not significantly change the time a neutron takes to be moderated and undergo fission. This is because, all other things being equal, the neutron will have to interact with the same number of nuclei in the slowing down process and the average distance between these nuclei will not have changed. The generation time sees a weak dependence on the nitric acid content and the temperature because both of these influence the average distance between the hydrogen and uranium nuclei which are involved in the slowing down and fission of the neutrons. The relationship with the H U ratio is stronger and more complex as this affects the degree to which a neutron will thermalise before causing fission. However, over the range observed, increasing the ratio always increases the generation time. This is because increasing this ratio means the average neutron undergoing fission will have a higher energy and so have been moderated fewer times by hydrogen nuclei meaning fewer collisions are required. A related reason is that the uranium nuclei have a much higher concentration and so neutrons of a given energy will have less distance to travel before they are captured by a uranium nucleus. The void fraction has a strong influence on the overall result as increasing the voidage increases the average distance between the nuclei the neutrons interact with while the atomic fractions of different isotopes are unchanged. We note that this approximation assumes the mean path length a neutron takes over its lifetime is not very much shorter than the separation between bubbles which make up the void's contribution to the volume. Both the delayed neutron fractions β i and the delayed neutron precursor decay rates λ i are weak functions of the hydrogen to uranium ratio only. This is because the change in moderation affects the energy spectrum of neutrons causing fission which affects the distribution of fission products including isotopes which are represented in the delayed neutron precursor groups. As a result, the dependency of these variables on the state of the system is only on the hydrogen to uranium ratio and then is only significant at high uranium concentrations. Several values of β i do not show any significant variation at all and will be treated as constant. The correlations for these variables are given in Equations (7) to (18): β 4 (t) = 0.00268 + 0.01 where β i is dimensionless and λ i has units of s −1 .

Results
Two scenarios are modelled in this paper. The first sees a supercritical over-moderated system undergoing a transient which evaporates a substantial amount of water from the solution, eventually causing the solution to become subcritical and halting the reaction.
In the second a subcritical under-moderated system has water added until the system becomes supercritical and a transient ensues. Further addition of water causes the system to eventually become over-moderated and the system eventually becomes subcritical. In both cases the longer term changes in reactivity occur due to the changing H U ratio. This effect is discussed in Thomas (1978) which shows how there is optimal value for this ratio in fissile solutions in terms of maximising reactivity, with k ef f decreasing as the H U ratio deviates further from this optimal ratio in either direction, as shown in Figure 1. This occurs because water acts as both a moderator and as an absorber. When the H U ratio is low the addition of more water causes increased moderation which is more important than the increased absorption but when the H U ratio is high there is ample hydrogen to moderate the system efficiently and adding more water does not cause significantly more efficient moderation but does cause an increase in absorption. In the first case the system begins with a H U ratio above optimal before it decreases to optimal and then to below optimal. In the second case the ratio H U begins below optimal before increasing to optimal and ends above optimal. Figure 1: A qualitative representation of the relationship between k ef f and H U for a fissile solution and the way the two simulated cases presented in this paper move through this space as the simulated time progresses.

Case 1: Step Reactivity Insertion
The first scenario to be studied with the model described in Section 2 is the case where the system begins at t = 0 with a significant positive reactivity due to the composition, mass and temperature of the system at this time, zero power and zero gas content (in terms of radiolytic gas and steam) and is in thermal equilibrium with its environment. This approximates the case where a large positive reactivity step is inserted into a previously subcritical cold system. A small source is present in this simulation and there is no addition of material once the simulation begins such thatṁ a (t) = 0. The simulated response to such a scenario is found in Figure 2. Initially the neutrons injected by the source begin to increase sharply in number due to the high reactivity. The power rises to a maximum of 1.25×10 9 W at 0.036s before the production of radiolytic gas reduces the reactivity of the system and causes the power to drop to around 1×10 6 W. At this power level the decay of delayed neutrons produced in the initial power peak produces enough neutrons to balance the neutron losses through the sub-criticality of the system and so the power holds relatively steady, decreasing only as the number of delayed neutron precursors decrease. On the time-scale of seconds the radiolytic gas produced in the initial power peak begins to leave the solution, increasing reactivity, and by 25s the system is critical again and the power has increased. The solution increased in temperature by approximately 10K in the first power peak and the elevated power after 25s causes significant heating to resume, which slowly reduces the reactivity and power. At 317s the solution temperature is above the saturation temperature of the solution and rapid steam production occurs. This causes a reduction of reactivity and power, which causes the steam production rate and therefore steam volume to drop after a few seconds. At this stage the power and temperature are fairly stable and the solution begins to evaporate, causing a reduction in mass and pH and an increase in uranium concentration. This causes a slow increase in reactivity as the system was initially over-moderated and the power peaks at 15.1kW at around 70,000s (compared to 12.6kW just after the onset of boiling). At around this time the evaporation of more water reduces the reactivity of the system as the system is now under-moderated. In the time up to 200,000s the steam and radiolytic gas content and the temperature all fall as the power slowly drops. This keeps the reactivity near zero and limits the rate at which the power may fall but, after the radiolytic gas and steam content have reached zero and temperature has reached 300K there is no more negative reactivity which can be removed from the system and the reactivity declines quickly as more water evaporates from the solution (this continues to occur because the air is modelled as having 50% humidity and, as a result, evaporation still occurs even when the solution is the same temperature as the air above it).

Case 2: Under-moderated Solution
The second scenario studied is that of an initially under-moderated subcritical solution to which water is steadily added. The aim of this simulation is to form a case analogous to the Y12 accident (Patton et al., 1958) where such an influx of water causes a uranyl nitrate solution to become supercritical and a criticality excursion to occur until the continued water addition caused over-moderation and the system became sub-critical again. It is stressed that this scenario is not intended to provide a simulation of the Y-12 accident itself but it is noted that there are strong qualitative similarities between this scenario and the accident. Again, the system initially begins at zero power and in thermal equilibrium with its environment and a small source is present. The initial mass of the solution is 137.5kg and water at room temperature (300K) is added at a rate of 0.05kg/s until the mass of the solution reaches 780kg such thatṁ a (t) is described by the equation: The initial reactivity of the system is -3.7$ but this soon rises as water is added until the system becomes critical at 5.3s. At this point the power begins to increase with the rate of increase rising substantially at 6.9s when the system becomes prompt supercritical. As the reactivity increase is a ramp instead of a step there is no power peak formed and the power rises fairly smoothly. The temperature and radiolytic gas content also rise slowly until 50s when the solution temperature exceeds the saturation temperature and steam begins to form. This causes a sudden reduction in power. Over the next 350s the steam content rises and then falls. This is because enough steam must be present in the system for the reactivity to be near zero and, following Equation (5), an increasing H U ratio causes the reactivity first to rise and then to fall as first the − At approximately 390s the power drops low enough that it cannot maintain the temperature of the solution at the saturation temperature against the dominant cooling effect of the influx of cold water. At this time the power begins to slowly decline as the increasing H U ratio reduces the reactivity faster than the cooling of the solution through the added material can raise it. The power is still substantial, however, and a significant amount of radiolytic gas is produced. There is more radiolytic gas present than earlier in the simulation because the value of k eff in Equation (5) is dependent on the void fraction not the actual volume of void and, as shown by Figure 3d the surface height has increased substantially, reflecting that the overall volume of the fuel solution/void mixture has increased. The power continues to fall at a rate governed by the decay of delayed neutron precursors until the end of the simulation. The inflow of water stops at around 715s and the temperature of the system begins to fall more slowly as the main medium of cooling has been removed and the temperature begins to tend towards the environment temperature as energy is lost through the sides of the system.

Conclusion
This paper has presented a model which allows evaporation of the system or the addition of material to change the chemical composition of a fissile solution undergoing a criticality excursion and has used correlations informed by MCNP simulations to simulate the effect of this changing composition on the transient. The examples of a system losing enough moderator through evaporation to cause it to become subcritical and the addition of water causing an initially under-moderated system to become critical and then sub-critical have been simulated. In both cases the results produced appeared physically plausible although no direct comparison to a physical system has been made. The effect of evaporation on the system becomes important for the evolution of the system between 1,000s and 10,000s as the rate of evaporation is fairly low, although modelling the effects considered in this paper are shown to be very important at all timescales when the addition of material is an important part of a scenario being simulated. This work has shown the feasibility and value of modelling the effect of changing solution composition over both short and long timescales in simulations of fissile solutions. Future work in this area could include the comparison of this model to accident scenarios or experiments, such as the CRAC or SILENE experiments, to verify the results of this model. The correlations used for the neutronics parameters and the evaporation rate could also be refined, particularly the correlation for the evaporation rate which currently has no dependence on the salt concentration. The addition of other physical processes important to the long term development of a transient, such as the production and decay of Xenon, would also make a valuable addition to this model.

Acknowledgements
The authors would like to thank EPSRC for their support through the following grants: Adaptive Hierarchical Radiation Transport Methods to Meet Future Challenges in Reactor Physics (EPSRC grant number: EP/J002011/1) and Nuclear Reactor Kinetics Modelling and Simulation Tools for Small Modular Reactor (SMR) Start-up Dynamics and Nuclear Critically Safety Assessment of Nuclear Fuel Processing Facilities (EPSRC grant number: EP/K503733/1).

Appendix A. MCNP Simulations
This appendix details the MCNP simulations performed to construct correlations for various neutronics parameters in Section 2.3. Note that the number of temperatures at which the simulations could be performed was limited by the number of temperatures the S(α, β) libraries were available within MCNP. Simulations at 293.6K were performed using the MCNP S(α, β) library lwtr.10, the 350K simulation using lwtr.11t and the 400K simulation with lwtr.12. The relatively small number of temperatures available is not expected to cause a significant error because, as discussed in , there is good indication that the key parameters such as the value of k eff are well approximated by linear functions of temperature.    The specific heat capacity of the material being added to the system c S The specific heat capacity of the solutioṅ E B (t) The rate at which energy is being used to create steam within the solutioṅ E side (t) The rate at which energy is lost through the sides of the container k ef f (t) The effective neutron multiplication factor of the system H U (t) The atomic ratio of hydrogen and uranium in the solution L s The latent heat of evaporation of water to steaṁ m a (t) The mass addition rate for material added to the systeṁ m e (t) The rate at which mass evaporates at the solution surface M HNO 3 (t) The concentration of the nitric acid m S (t) The mass of the solution N S,acid The mass fraction of the nitrogen contained in the nitric acid only P (t) The power produced by the system p v (t) The vapour pressure of the solution p wa The partial pressure of water in the air above the solution T a The temperature of the material being added to the system T S (t) The temperature of the solution (assumed homogeneous) U S (t) The uranium mass fraction of the solution V F S (t) The void fraction of the solution/void mixture β i The delayed neutron fraction relating to the ithe precursor group λ i The decay rate of a delayed neutron precursor in the ith precursor group Λ(t) The generation time of the system ρ S (t) The density of the solution Table B.5: A description of the variable and parameters