Rod Ejection and Drop Accident Analysis of Aqueous Homogeneous Solution Reactor

Aqueous Homogeneous Solution Reactor concept has been proposed for producing medical isotopes (Medical Isotope Production Reactor-MIPR). However, there are several difficulties in transient calculation of aqueous homogeneous solution reactors. First, there are no assemblies in the core which is different from the traditional reactor core. Second, the operation of aqueous solution reactor at a power of 200 kW will generate radiolytic-gas bubbles. The void volume created by these bubbles in the solution core will introduce a strong negative reactivity feedback. Third, the complex structure of the coolant pipes immersed in fuel solution requires unstructured neutron diffusion calculation methods. Therefore, analytic basis functions expansion nodal method for arbitrary triangular-z node is established to solve the complex structure geometry neutron diffusion equation. Based on this, a software named TABFEN-K has been developed to solve the three-dimensional space-time neutron kinetic equations. Then, TABFEN-K code is used for typical accident analysis of a solution reactor. A simplified geometry model, bubbles generation model, thermal conduction model, and cross section feedback model are given in this paper. A software called TABFEN-MIPR is developed and used for the simulations of the control rod ejection and drop. The same characteristics in the transient process with the results from literatures are obtained.

Introduction

Other than the irradiation of uranium targets in heterogeneous reactors, aqueous homogeneous solution reactors is an alternative way to produce medical isotopes (Ball, 1997). Based on many potential advantages of aqueous homogeneous solution reactors, the Medical Isotope Production Reactor (MIPR) concept has been proposed for producing medical isotopes.

However, many difficulties exist in transient calculation of aqueous homogeneous solution reactors (Briesmeister, 2002). First, there are no assemblies in the core which is very different from the traditional reactor core. Second, the operation of aqueous solution reactor at a power of 200 kW will generate radiolytic-gas bubbles. The void volume created by these bubbles in the solution core will introduce a strong negative reactivity feedback (Dunenfeld and Stitt, 1963). Third, the complex structure of the coolant pipes immersed in fuel solution requires unstructured neutron diffusion calculation methods (Lecorche and Seale, 1973).

Therefore, analytic basis functions expansion nodal method for arbitrary triangular-z node is established to solve the complex structure geometry neutron diffusion equation. The distribution of detailed neutron flux within each node is expanded into the sum of a set of analytic basis functions by accurately formulating the multi-group matrix form neutron diffusion equation and appropriately choosing the expansion order and characteristically directions. To balance between accuracy and efficiency, adjacent nodes are coupled with each other through the zero- and first-order partial current moments defined at the interface. Coordinate conversion technique is used to transform arbitrary triangle into a regular triangle to simplify the derivation of the response matrix between the incoming and outgoing partial currents. A new sweeping scheme is designed for the triangular-z node to solve the nodal diffusion equation iteratively. Then, this analytic basis function expansion nodal method is extended for solving the space-time neutron kinetic equations. The efficient fully implicit difference method is applied to discretize the time variable of the prompt neutron diffusion equation together with the time integration method employed to handle the delayed neutron precursor equations. At each step after the time variable discretization, a fixed source neutron diffusion equation is obtained. Different from the steady-state case, an intranodal neutron flux consists of two components, they are an analytic basis function expanded homogeneous terms and a heterogeneous term which can be approximated into a second-order orthogonal polynomial. Similar to the steady-state, nodal averaged volume flux and surface partial current moments are calculated by introducing the coordinate conversion scheme. Based on this, a code named TABFEN-K has been developed to solve the three-dimensional space-time neutron kinetic equations.

At last, the TABFEN-K code is used for typical accident analysis of a solution reactor. The simplified geometry model, bubbles generation model and the thermal conduction model of solution reactor are given together with the cross section feedback model. A code called TABFEN-MIPR is developed and used for the simulations of the control rod ejection and drop accident. The same characteristics in the transient process with the results from literatures are obtained.

Calculation Model

The simplified geometry model, bubbles generation model, and the thermal conduction model of solution reactor are given together with the cross section feedback model.

Solution of the Diffusion Equation

In a particular homogeneous mesh, the matrix form of the standard multigroup neutron diffusion equation in triangular geometry can be written as (Noh and Cho, 1993, 1994; Nam Zin and Noh, 1995):

$\begin{array}{ll}-{\nabla }^2\Phi (\text{r})+\Sigma (k_{\text{eff}})\Phi (\text{r})=0& (1)\end{array}$

where Φ(r) = col{Φ₁(r), Φ₂(r), …, Φ_G(r)} is the neutron flux and Φ_G(r) denotes the G-th group neutron flux; the matrix Σ (k_eff) is a square matrix of order G with its elements ${\sum }_{g{g}^{\prime }}$ as follows:

$\begin{array}{ll}{\sum }_{g{g}^{\prime }}=\left({\delta }_{g{g}^{\prime }}{\displaystyle \sum _t^g}-{\displaystyle \sum _s^{g^{\prime }-g}}-\frac{{\chi }^g}{k_{\text{eff}}}v{\displaystyle \sum _f^{g^{\prime }}}\right)/D^g& (2)\end{array}$

where ${\delta }_{g{g}^{\prime }}$ is the Kronecker symbol and the other notations are as in conventional manner.

The solution of the diffusion Equation (1) depends on whether the eigenvalues λ_m and corresponding eigenvectors u_m of the matrix Σ(k_eff) are real or complex. For simplicity, we present our derivation here only for the case in which the eigenvalues are real. In fact, the eigenvalues are exactly real for two-group problem. Then, Ψ(r) is introduced as:

$\begin{array}{ll}\text{U}=\left[{\text{u}}_1,{\text{u}}_2,\dots ,{\text{u}}_G\right]& (3)\end{array}$

$\begin{array}{ll}\text{Ψ}(\text{r})={\text{U}}^{-1}\Phi (\text{r})& (4)\end{array}$

Substitute Equation (4) to Equation (1) to reduce Equation (1) to a decoupled form

$\begin{array}{ll}-{\nabla }^2{\text{Ψ}}_m(\text{r})+{\lambda }_m{\text{Ψ}}_m(\text{r})=0,m=1\dots G& (5)\end{array}$

Equation (5) is a typical Helmholz equation with analytical solution

$\begin{array}{ll}{\text{Ψ}}_m(\text{r})={\displaystyle \sum _{l=1}^{\infty }}\left[A_{ml}{e}^{\sqrt{{\lambda }_m}{\text{e}}_l\text{r}}+B_{ml}{e}^{-\sqrt{{\lambda }_m}{\text{e}}_l\text{r}}\right]& (6)\end{array}$

where e_l is an arbitrary unit vector. In practice, the number of terms l in Equation (6) mainly depends on the number of the available nodal boundary conditions. Generally speaking, the higher of the terms result in higher accuracy but lower calculation speed, vice versa. For triangular geometry, 3 interfaces' averaged partial current and first-order partial current moment are chosen as the coupling condition to determine the expansion coefficients. Therefore, Equation (6) can be written as:

$\begin{array}{ll}{\Psi }_m\left(\text{r}\right)=\{\begin{array}{c}{\displaystyle \sum _{l=1}^3{A}_{ml}Sinh\left(k_m{\text{e}}_l\text{r}\right)}+{\displaystyle \sum _{l=1}^3{B}_{ml}Cosh\left(k_m{\text{e}}_l\text{r}\right)},{\text{λ}}_m>0\\ {\displaystyle \sum _{l=1}^3{A}_{ml}Sin\left(k_m{\text{e}}_l\text{r}\right)}+{\displaystyle \sum _{l=1}^3{B}_{ml}Cos\left(k_m{\text{e}}_l\text{r}\right)},{\text{λ}}_m<0\end{array}& (7)\end{array}$

where $k_m=\sqrt{|{\lambda }_m|}$, by using Equation (5) and Equation (7), the intranodal flux distribution can be eventually obtained as:

$\begin{array}{ll}{\Phi }_g\left(\text{r}\right)=\{\begin{array}{l}{\displaystyle \sum _{m=1}^G{u}_{gm}\left({\displaystyle \sum _{l=1}^3{A}_{ml}Sinh\left(k_m{\text{e}}_l\text{r}\right)}+{\displaystyle \sum _{l=1}^3{B}_{ml}Cosh\left(k_m{\text{e}}_l\text{r}\right)}\right)},\hfill \\ {\text{λ}}_m>0\hfill \\ {\displaystyle \sum _{m=1}^G{u}_{gm}\left({\displaystyle \sum _{l=1}^3{A}_{ml}Sin\left(k_m{\text{e}}_l\text{r}\right)}+{\displaystyle \sum _{l=1}^3{B}_{ml}Cos\left(k_m{\text{e}}_l\text{r}\right)}\right),}\hfill \\ {\text{λ}}_m<0\hfill \end{array}& (8)\end{array}$

where u_gm is the element of the matrix U defined in Equation (3).

Based on this, a code named TABFEN-K has been developed for solving the three-dimensional space-time neutron kinetic equations. The code is verified by 3D-LMW, 3D-LRA benchmarks. Both the relative power and power distribution are in good agreement with the benchmarks (Ban et al., 2012).

Void Bubbles Model

In aqueous homogeneous solution reactors, void bubbles affect the core state by changing the density and volume of fuel solution. The change of fuel solution density was considered by interpolating macroscopic cross-sections in terms of void fraction, while the volume change of fuel solution was considered by modifying the height of fuel solution surface.

In the volume calculation of void bubbles, the total void volume is given by (Souto et al., 2005).

$\begin{array}{ll}{\bar{V}}_B=\left(1+\frac{1}{\xi }\right)\frac{G(H_2)}{N_A}{R}_g{T}_g\left(\frac{r_B}{2\sigma }\right){\tau }_B·\bar{n}& (9)\end{array}$

where

${\bar{V}}_B$ : total void volume (m³);

$\bar{n}$ : the reactor power (MW);

G(H₂) : hydrogen yield in fuel solution (J⁻¹);

ξ : fraction of H₂ molecules per O₂ molecule produced by water radiolysis;

N_A : Avogadro's constant (mol⁻¹);

R_g : gas constant (J·mol⁻¹·K⁻¹);

T_g: void bubble temperature (K);

r_B : void bubble characteristic radius (um);

σ : fuel solution surface tension (Kg·s⁻¹);

τ_B : the characteristic time for void bubble release from the fuel solution (s).

The distributions of void bubbles in radial and axial direction are different. In radial direction, void distribution is decided by radial power distribution; while in axial direction by the law of axial void movement. In this paper, it assumes that void distribution in axial direction obeys integral power distribution.

Fuel solution volume expansion is caused by generation of void bubbles and thermal expansion and contraction of the fuel solution. The increment of the fuel solution volume can be determined by

$\begin{array}{ll}\Delta V=V_{Void}+\Delta V_{Hot}=V_{Void}+\left(\frac{{\rho }_{Cold}}{{\rho }_{Power}}-1\right)·V_{Cold}& (10)\end{array}$

where

V_Void : total void volume (cm³);

ΔV_Hot : solution volume increment caused by thermal expansion and contraction (cm³);

ρ_Cold : density of cold solution (g/l);

ρ_Power : density of hot solution (g/l);

V_Cold : volume of cold solution (cm³).

Temperature Model

The solution of the fluidity of liquid fuel, combined with the solution of smaller volume, so that the temperature of the solution of liquid fuel is equal, then according to the quantity of heat production and export, the calculation formula of fuel solution can be written as:

$\begin{array}{ll}\frac{dT(t)}{dt}=\frac{1}{c_f{m}_f}p(t)-\frac{1}{{\tau }_T}(T(t)-T_0)& (11)\end{array}$

where

c_f : specific heat capacity (J/Kg/K);

m_f : Solution fuel weight (kg);

p(t) : Power (W);

τ_t : Characteristics of the time (s);

T(t) : volume of cold solution (K).

Few-Group Constants Fitting

The few-group constants are fitted with the following form

$\begin{array}{ll}{\Sigma }_x=f_1(Bu)f_2(v)f_3(p)+f_4(Bu)·\delta (Cr)& (12)\end{array}$

where

f_i(x) (i = 1,2,3,4): a polynomial of variable x;

δ(Cr) : a dirac function which equals to 1 or 0 when control rods exist or not;

Bu : burnup (MW·d/tU);

v : void fraction;

p : relative power.

At last, the TABFEN-K code is used for typical accident analysis of a solution reactor. The simplified geometry model, bubbles generation model and the thermal conduction model of solution reactor are given together with the cross section feedback model. A code called TABFEN-MIPR is developed.

Numerical Results

The structure and material of the aqueous homogeneous solution reactor are given in Figures 1, 2 gives the radial mesh of the core used in the calculation, respectively.

FIGURE 1

Figure 1. Structure of MIPR.

FIGURE 2

Figure 2. The radial meshing map of MIPR.

The code TABFEN-MIPR is used for the simulations of the control rod ejection and drop accident. The same characteristics in the transient process with the results from literatures are obtained. At the beginning, the power of MIPR is 0.033MW, the temperature of the fuel solution is 344.4164K, the volume of bubble is 376.2742 cm³. The control rods are at 25.68 cm from the bottom of the core.

Rods Ejection Accident

Rods ejection accident belongs to the limit accident, accident control rod in the pop-up core in the process of introducing positive reactivity, which can cause power has arisen sharply, the safety of the pile is affected by the great, usually is the focus of the nuclear safety analysis of the research object. For rods ejection accident, the control rods eject to the top of the core at 32.05 cm, in this case 1,143 pcm reactivity are introduced. Figure 3 gives the solution of ejection rod in the change of the power, Figure 4 gives the corresponding solution temperature and bubble content changes over time. We can find that, after the control rod pop-up, due to the introduction of larger positive reactivity, power rises rapidly in a short time, because of the heat transfer lag, the temperature of the fuel and air bubble content is slowly rising, then began to rapid increase in power; The increment of fuel temperature and bubble content, in turn, has affected the nuclear reaction cross section, the comprehensive result of negative feedback effect to power up to a certain value began to decline, with declining power, and the solution temperature and content of bubble growth began to slow, finally the change of power tend to flat, and finally stability on above the level of the power of the initial time, temperature, and bubble content at this time also started to become stable, and keep on the above the level of the initial time. Because of solution temperature and bubble content, are under the influence of power, has obvious lag in time, Figure 4 also shows that the change of the solution temperature and air bubbles curve shape is consistent, and lagging behind the growth of the power, hence TABFEN-MIPR calculation results conform to play well the transient progress regularity of the accident.

FIGURE 3

Figure 3. Relative power variation over time.

FIGURE 4

Figure 4. Temperature variation over time.

Rod Drop Accident

For rods drop accident, the control rods drop to the core at 20 cm, in this case −1,733 pcm reactivity are introduced. The negative reactivity caused power drops rapidly, but due to the solution of temperature and air bubbles in the feedback, is likely to lead to power fluctuations, may also cause reactor shutdown, therefore must also be analyzed. The results show that, at the beginning of the rod drop, core into larger negative reactivity, core power fells sharply; With declining power, temperature, and decreases in the bubble content and temperature of the solution and the bubble content is introduced to reduce the reactivity of the reactivity makes for control rod falling into negative reactivity decreases, so power after plunging began to slow down. As the accumulation of feedback effects, stable power eventually below the initial power of the power level. And solution temperature and bubble content change lag, which makes them relatively flat. Therefore, TABFEN—MIPR calculation result is reasonable and accords with the transient progress law of rod drop accident.

Conclusions

In this study, TABFEN-K software is used for typical accident analysis of a solution reactor. The simplified geometry model, bubbles generation model and the thermal conduction model of solution reactor are given together with the cross section feedback model. The code TABFEN-MIPR is developed and used for the simulations of the control rod ejection and drop accident. The same characteristics in the transient process with the results from literatures are obtained as well.

Data Availability Statement

All datasets generated for this study are included in the article/supplementary material.

Author Contributions

KW developed the code and wrote the manuscript. CX performed the calculation. CL edited the manuscript. YL and JJ analyzed the data. HW gave advice and comments. CX, CL and YL gave the technical support. HW and JJ gave advice for this project.

Funding

The research was supported by National Science and Technology Major Project of China (2019ZX06005001).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Ball, R. M. (1997). Medical Isotope Production Reactor [P]. US: 5596611

Google Scholar

Ban, Y., Endo, T., and Yamamoto, A. (2012). A unified approach for numerical calculation of space-dependent kinetic equation[J]. J. Nucl. Sci. Technol. 49, 496–515. doi: 10.1080/00223131.2012.677126

CrossRef Full Text | Google Scholar

Briesmeister, J. F. (2002). MCNP-A General Monte Carlo N-Particle Transport Code. Version 4C, LA-13709–M.

Google Scholar

Dunenfeld, M. S., and Stitt, R. K. (1963). Summary Review of the Kinetics Experiments on Water Boilers. NAA-SR-7087. Atomics International.

Google Scholar

Lecorche, P., and Seale, R. L. (1973). Review of the Experiments Performed to Determine the Radiological Consequences of a Criticality Accident. Y-CDC-12. Oak Ridge National Laboratory.

Google Scholar

Nam Zin, C., and Noh, J. M. (1995). Analytic function expansion nodal method for hexagonal geometry. Nuclear Sci. Eng. 121, 245–253. doi: 10.13182/NSE95-A28561

CrossRef Full Text | Google Scholar

Noh, J. M., and Cho, N. Z. (1993). A new diffusion nodal method based on analytic basis function expansion. Transac. Am. Nucl. Soc. 69, 462–463.

Google Scholar

Noh, J. M., and Cho, N. Z. (1994). A new approach of analytic basis function expansion to neutron diffusion nodal calculation. Nucl. Sci. Eng. 116, 165–180. doi: 10.13182/NSE94-A19811

CrossRef Full Text | Google Scholar

Souto, F. J., Kimpland, R. H., and Heger, A. S. (2005). Analysis of the effects of radiolytic-gas bubbles on the operation of solution reactors for the production of medical isotopes. Nucl. Sci. Eng. 150, 322–335. doi: 10.13182/NSE05-A2519

CrossRef Full Text | Google Scholar