% ****** Start of file aipsamp.tex ****** % % This file is part of the AIP files in the AIP distribution for REVTeX 4. % Version 4.1 of REVTeX, October 2009 % % Copyright (c) 2009 American Institute of Physics. % % See the AIP README file for restrictions and more information. % % TeX'ing this file requires that you have AMS-LaTeX 2.0 installed % as well as the rest of the prerequisites for REVTeX 4.1 % % It also requires running BibTeX. The commands are as follows: % % 1) latex aipsamp % 2) bibtex aipsamp % 3) latex aipsamp % 4) latex aipsamp % % Use this file as a source of example code for your aip document. % Use the file aiptemplate.tex as a template for your document. \documentclass[% aip, %jmp,% %bmf,% sd,% %rsi,% amsmath,amssymb, %preprint,% reprint,% %author-year,% author-numerical,% twoside ]{revtex4-1} \usepackage{graphicx}% Include figure files \usepackage{dcolumn}% Align table columns on decimal point \usepackage{bm}% bold math \usepackage[mathlines]{lineno}% Enable numbering of text and display math %\linenumbers\relax % Commence numbering lines \usepackage{hyperref} \usepackage{fancyhdr} \usepackage{appendix} \usepackage{ifthen} \usepackage{float} \usepackage{multirow} \usepackage{tabularx} \usepackage{array} \pagestyle{fancy} \fancyhf{} \lhead{\textit{Control rod worth prediction in a Materials Test Reactor}} \fancyfoot[RO]{\textit{\small{McMaster Journal of Engineering Physics}, 2017, \bf{[2]}, \thepage}} \fancyfoot[LE]{\textit{\small{McMaster Journal of Engineering Physics}, 2017, \bf{[v2]}, \thepage}} \renewcommand{\footrulewidth}{0.4pt} %\fancyfoot{\ifthenelse{\value{page}=1}{\small{\textit{McMaster Journal of Engineering Physics}, 2016, \bf{[vol]}, \thepage}}} \newcolumntype{P}[1]{>{\centering\arraybackslash}p{#1}} \begin{document} \preprint{AIP/123-QED} \title[\textit{Sample manuscript short title}]{Control rod worth prediction in a Materials Test Reactor by finite difference solution of the two-group point reactor kinetics equations} \author{J. Saraco} %only include a), b) etc if there are footnotes on the authors (see section III.3 Footnotes, and Notes and References at end \affiliation{ McMaster University, 1280 Main St W., Hamilton, Ontario, Canada. %\\This line break forced with \textbackslash\textbackslash } \author{D. FitzGreen} \date{\textbf{\today}}% It is always \today, today, % but any date may be explicitly specified \begin{abstract} A finite difference model was developed in order to predict control rod worth from a set of power measurements following an insertion of reactivity. The model is based on the solution of the point reactor kinetics equations in the two energy group formalism where thermal reactivity is varied to model the gravitational (and thus quadratic) insertion of a control rod in a scram scenario. The model was tested using power measurements from the McMaster Nuclear Reactor during the insertion of two different control rods to shut down the reactor from a low initial power. The model was found to agree with the worths of the two control rods within an order of magnitude, and fail to reproduce the decay rate in power observed for sufficiently high reactivity insertions, leading to the belief that lack of feedback effects constitutes the biggest failure in the model currently. %Valid PACS numbers may be entered using the \verb+\pacs{#1}+ command. \end{abstract} \pacs{28.50.Dr, 28.52.Av, 02.60.Cb, 02.60.Dc, 02.60.Lj, 02.70.Bf}% PACS, the Physics and Astronomy % Classification Scheme. \keywords{nuclear, reactor, control rod, worth, finite difference}%Use showkeys class option if keyword %display desired \maketitle \thispagestyle{fancy} \section{\label{sec:level1}Introduction} In order for a nuclear reactor to be ran safely and efficiently, there must be extensive knowledge of the negative reactivity at an operator's disposal. The most valuable source of negative reactivity comes from control rods usually meant to shutdown the reactor, or routinely offset the initial positive reactivity of the reactor. The worth of these control rods must be known accurately to ensure a cold shutdown at every stage of reactor operation, including during an unexpected insertion of positive reactivity due to an accident scenario. \\ I will define the following reactor concepts that are required in this research. First, reactivity of a reactor core. Reactivity is a measure of the growth or decay in the nuclear chain reaction after subsequent neutron generations. If a fission neutron, on average, produces more than one fission neutron, then the reactor is said to be supercritical, and thus has a reactivity greater than zero. The nuclear reactor, if allowed to continue in this state, will exponentially grow in power. If a fission neutron, on average, produces one fission neutron, then the reactor is critical and the power is stable. If a fission neutron produces less than one neutron in the chain reaction, the reactor is said to be subcritical and the power will decay over time. Most reactors are designed to have excess positive reactivity in the core, leading to a supercritical reactor. To fix this, there are control rods that are inserted into the core that are made out of neutron absorbing materials. The insertion of these will add negative reactivity, and allow a critical state to be reached. As reactor fuel is burned, the positive reactivity they contribute to the reactor decreases. This also contributes to designing an initial state of positive reactivity. An initially critical reactor core, as it's run, would very quickly reach a subcritical state and would require positive reactivity insertions to keep power stable or perform fine adjustments. Since it is safer to be capable of shutting down a reactor from the insertion of a control rod, modern reactor designs have settled on fuel bundles that have excess positive reactivity that are then adjusted with negative reactivity control rods that can be adjusted to keep the reactor stable. \\ For most research or commercial reactors, control rod worth is measured using computer systems that are integrated with the reactor control systems during routine tests, or approximated by computer codes such as RELAP5. These methods work well with the required resources, however there are cases where this may be difficult to achieve, or are disallowed from use. Within countries that have not ratified nuclear non-proliferation treaties (Pakistan and Israel, for example), reactor operators are disallowed from using reactor codes and do not receive support from the International Atomic Energy Agency (IAEA). Political reasons aside, it is in international interest for these reactors to have the resources to run safely and efficiently. \\ For this reason, a small, efficient mathematical model for measuring control rod worth with a dataset of power measurements during a shutdown has been developed. The model only requires a set of power measurements, and as such does not need to be built into the larger reactor system. The aim of the model is to be accessible, easy to use, and with the least mathematical complexity possible. \section{\label{sec:mod}Mathematical Model} \subsection{Summary of Finite Difference Method} To model the reactor under study, the set of power data measured is fit against a solution of the point reactor kinetics equations. These are usually expressed in the one group formalism, where all neutrons are constrained to a single value of energy. The rod worth estimation model, in this case, uses the reactor kinetics equations in the two-group formalism. \\ Neutrons are modelled as being born from fissions in the fast group (having an energy greater than $0.625$ eV in most derivations), they scatter down to the thermal group (having energy less than $0.625$ eV), and then induce another fission, continuing the chain reaction. The set of differential equations is given by:\cite{aboanber_novel_2014} \begin{align} \frac{dn_1(t)}{dt} &= -\left(\frac{\sum_{i=1}^{8}\beta_i}{l_1}+\kappa\right)n_1(t) \notag \\ &+\left(\frac{1-\sum_{i=1}^{8}\beta_i}{l_2}\right)n_2(t)+\sum_{i=1}^{8}\lambda_iC_i(t) \end{align} \begin{align} \frac{dn_2(t)}{dt} = \kappa n_1(t)+\left(\frac{\rho_2(t)-1}{l_2}\right)n_2(t) \end{align} \begin{align} \frac{dC_i(t)}{dt} = \frac{\beta_i}{l_1}n_1(t)+\frac{\beta_i}{l_2}n_2(t)-\lambda_iC_i(t) \end{align} Where $n_1(t)$ is the fast neutron concentration, $n_2(t)$ is the thermal neutron concentration, and $C_i(t)$ is the ith precursor concentration. \\ The solution of these equations are obtained by performing an implicit finite difference discretization, and solving the resulting system of equations at each time step. Upon discretizing the system separating points in time by $\Delta t$, the following matrix is obtained, simplified for only two precursor groups: \begin{widetext} \[ \begin{bmatrix} 1+\frac{\Delta t}{l_1}\left(\sum_{i=1}^{2}\beta_i+\kappa l_1\right) & \frac{\Delta t}{l_2}\left(\sum_{i=1}^{2}\beta_i+1\right) & -\Delta t\lambda_1 & -\Delta t\lambda_2 \\ -\Delta t\kappa & 1+\frac{\Delta t}{l_2}\left[1-\rho_2(t_i)\right] & 0 & 0 \\ -\frac{\Delta t\beta_1}{l_1} & -\frac{\Delta t\beta_1}{l_2} & 1+\Delta t\lambda_1 & 0 \\ -\frac{\Delta t\beta_2}{l_1} & -\frac{\Delta t\beta_2}{l_2} & 0 & 1+\Delta t\lambda_2 \end{bmatrix} \begin{bmatrix} n_{1(i+1)} \\ n_{2(i+1)} \\ C_{1(i+1)} \\ C_{2(i+1)} \end{bmatrix} = \begin{bmatrix} n_{1(i)} \\ n_{2(i)} \\ C_{1(i)} \\ C_{2(i)} \end{bmatrix} \] \end{widetext} By iterating through the time steps until the final time is reached, the functions are built up and can be plotted or analyzed after this point. \subsection{Reactor Parameters} For this model to approximate the behaviour of a Materials Test Reactor, the reactor parameters must be accurate to this system. The parameters have the following values: \begin{table}[H] \centering \caption{Two group reactor parameters} \label{paramTable} \def\arraystretch{1.3} \scalebox{1}{ \begin{tabular}{|l|l|l|} \hline \textbf{Reactor Parameter} & \textbf{Group} & \textbf{Value} \\ \hline \multirow{2}{*}{Diffusion Coefficient (D)} & Fast & 1.35 cm \\ \cline{2-3} & Thermal & 1.08 cm \\ \hline \multirow{2}{*}{Absorption Cross Section ($\Sigma_a$)} & Fast & 0.001382 $\text{cm}^{-1}$ \\ \cline{2-3} & Thermal & 0.0054869 $\text{cm}^{-1}$ \\ \hline \multirow{2}{*}{Fission Cross-Section ($\Sigma_f$)} & Fast & 0.000242 $\text{cm}^{-1}$ \\ \cline{2-3} & Thermal & 0.00408 $\text{cm}^{-1}$ \\ \hline \multirow{2}{*}{Scattering Cross-Section ($\Sigma_{s}$)} & Fast & 0.0023 $\text{cm}^{-1}$ \\ \cline{2-3} & Thermal & 0 $\text{cm}^{-1}$ \\ \hline \multirow{2}{*}{Neutron Speed} & Fast & $3.0\cdot10^7\text{m/s}$ \\ \cline{2-3} & Thermal & $2.2\cdot10^5\text{m/s}$ \\ \hline \end{tabular}} \end{table} Since the Material Test Reactor runs with uranium-235 as fuel, the eight sets of precursor constants are the following: \begin{table}[H] \centering \caption{Eight precursor constants for uranium-235} \label{precurTable} \def\arraystretch{1.3} \scalebox{0.87}{ \begin{tabular}{|c|c|c|} \hline Precursor Group & Delayed Neutron Fraction ($\beta$) & Decay Constant ($\lambda$) \\ \hline 1 & 0.0218 & 0.012467 s$^{-1}$ \\ 2 & 0.1023 & 0.028292 s$^{-1}$ \\ 3 & 0.0605 & 0.042524 s$^{-1}$ \\ 4 & 0.131 & 0.133042 s$^{-1}$ \\ 5 & 0.22 & 0.292467 s$^{-1}$ \\ 6 & 0.06 & 0.666488 s$^{-1}$ \\ 7 & 0.054 & 1.634781 s$^{-1}$ \\ 8 & 0.0152 & 3.5546 s$^{-1}$ \\ \hline \end{tabular}} \end{table} \subsection{Modelling a Reactivity Insertion} Reactor control systems often place their emergency control rods above the reactor core. This allows the system to keep the control rods lifted with power, and in the event of loss of power, the actuators disengage and the control rods fall into the core. For this reason, a reactivity insertion will be modelled as a negative quadratic from 0 reactivity to the full worth of the rod. It is defined that it takes the control rod two seconds to be fully inserted in the reactor. If the control rod is inserted at time $t_0$ when the reactor is critical, the reactivity over time will be given by: \begin{align} \rho(t) = \frac{\rho_{rod}}{4}(t-t_0)^2 \end{align} \begin{figure} \includegraphics[width=260pt]{reactivityPlot.png}% Here is how to import PDF art \caption{\label{fig:reactPlot} Reactivity over time for a simulated rod insertion due to gravity} \end{figure} Figure \ref{fig:reactPlot} is a graph of the reactivity over time in the reactor for a modelled quadratic insert of a control rod. With power measurements of a reactor subject to a scram, the relative power is calculated (power divided by the power at the time of the scram), and $rho_2(t)$ is varied to find the solution of the relative thermal neutron density that fits the relative power curve the best way. Thermal flux is used because a Materials Test Reactor is a thermal reactor, and it is approximated that thermal flux is proportional to power, and thus relative thermal flux is equal to relative power. \subsection{Curve Fit to Power Data} Following an insertion of reactivity, there are two regimes of transience that occur. There are two sources of neutrons in a reactor, prompt neutrons that come directly from fission and delayed neutrons that come from the decay of daughter nuclei. In a very short time following an insertion of positive or negative reactivity, the number of prompt neutrons (which account for a large percentage of total neutrons) experience a large, and rapid change. This is called the prompt jump.\\ To predict rod worth, the variable $\rho_{rod}$ in the expression for $\rho_2(t)$ is varied until the calculated prompt jump most accurately follows the observed prompt jump.\\ While testing the model, the decision to use the prompt jump only to fit and not the entire dataset came from an inability to replicate the decay rate of the observed power data (this will be seen in the Model Testing section). This could be explained by inaccurate reactor parameters, or feedback effects that dampen the decay of the power and are not present in the idealized point reactor kinetics equations. \subsection{Avoiding Stiffness Effects} The point reactor kinetics equations, derived in any number of energy or precursor groups, are known as a "stiff" system of differential equations. A stiff differential equation is one that, when solving with computational methods, produces unstable solutions with sufficiently low time step.\cite{shampine_users_1979} To prove that the model developed does not produce unstable time steps with solution, the square of the error between two solutions set at different time steps was plotted. Figure \ref{fig:reactPlot} shows that the highest error that occurs is on the order of 1\%, and is thus stable for time steps fewer than $10^{-1}$ s. \begin{figure} \includegraphics[width=260pt]{errorPlot.png}% Here is how to import PDF art \caption{\label{fig:errorPlot} Relative error between solutions for a sinusoidal reactivity for $\Delta t=4\cdot10^{-6}$ s and $\Delta t=10^{-1}$ s.} \end{figure} \section{\label{sec:res}Model Testing} Two experiments were performed at the McMaster Nuclear Reactor. From an initial power of about 500 W thermal, two scrams were done with two different control rods, and power measurements were taken. With the power data sets obtained, the relative power was calculated by visually inspecting the plot and determining the point in time of control rod insertion, and dividing all power values by the power at this point. \\ Then, the model was tested against these two relative power data sets by iterating through values of $\rho_{rod}$ inserted quadratically. The value that produced the most accurate agreement between the prompt jump observed and the prompt jump calculated in the solution of the system of differential equations is the estimated control rod worth. \begin{figure} \includegraphics[width=260pt]{rod3plot.png} \caption{\label{fig:rod3Plot} Curve fit of the solution of the two-group PRK model for $\rho_{rod} = -4.28$\text{ mk} with scram data measured from the McMaster Nuclear Reactor from an initial power of about 500 W-thermal} \end{figure} \begin{figure} \includegraphics[width=260pt]{rod5plot.png} \caption{\label{fig:rod5Plot} Curve fit of the solution of the two-group PRK model for $\rho_{rod} = -2.15$\text{ mk} with scram data measured from the McMaster Nuclear Reactor from an initial power of about 500 W-thermal} \end{figure} Using the model to predict rod worth of control rod 3 in the McMaster Nuclear Reactor produces figure \ref{fig:rod3Plot}, and results in an estimated control rod worth of $-4.28$ mk. A plot for control rod 5 is found in figure \ref{fig:rod5Plot}, and results in an estimated control rod worth of $-2.15$ mk. \section{\label{sec:conc}Conclusions} Even though we were not given the true worths of control rods 3 and 5 in the MNR, we were told by operators that our estimates were the right order of magnitude, and that rod 3 has more negative reactivity than rod 5 which was predicted.\\ Theoretically, this method for predicting rod worth should be as accurate as analytic equations that predict rod worth using measured constants like reactor period. This assumes that the error introduced from the solution method being computational and not analytic is negligible, which is not proven in this paper. In figure \ref{fig:rod3Plot} in comparison to figure \ref{fig:rod5Plot}, it can be seen that the decay rate after the prompt jump is much more accurately predicted. Without any study into this matter, it is predicted that Material Test Reactors have feedback effects that are more pronounced for higher insertions of reactivity. Since the model has no feedback effects accounted for, it would make sense that more prevalent feedback effects would cause higher deviation and thus a lower accuracy in the model. If this conclusion is true, then the biggest improvement to this model would come in the form of using higher order reactor kinetics equations, such as incorporating xenon transience or feedback effects (power coefficient, temperature coefficient, effects etc.). \begin{acknowledgments} Thanks to Daniel FitzGreen, lab technician for the nuclear engineering stream at of Engineering Physics and McMaster University, for guiding my research into this problem and performing the power measurements at the McMaster Nuclear Reactor that were used to test the model. \end{acknowledgments} \appendix \subsubsection*{\label{sec:notes}Notes and References} %note that the section numbering is suppressed by the use of "*" before the name of the section \footnotesize \noindent Author to whom correspondence should be addressed by alectronic mail: saracojc@mcmaster.ca %\nocite{*} \bibliographystyle{unsrtnat} \bibliography{MJEP}% Produces the bibliography via BibTeX. \end{document} % % ****** End of file aipsamp.tex ******