main.tex

\documentclass[letterpaper]{mc2025}

%%% Packages Required by Class (already included)
% fancyhdr
% lastpage
% titling
% titlesec
% ragged2e
% enumitem
% amsmath
% graphicx
% geometry
% newtxtext
% newtxmath
% hyperref
% cleveref
% caption
% authblk
% apptools
% appendix
% ifpdf
% epstopdf

%%% Some other useful packages
% \usepackage{tikz}
% \usepackage{color}
% \usepackage{subcaption}
% \usepackage{algcompatible}
% \usepackage{bm}
% \usepackage{array}

\usepackage[acronym]{glossaries}
%\makeglossaries
\include{acros}

\usepackage{xspace}
%\usepackage{graphicx}
\graphicspath{{images/}}

\usepackage{placeins}
\usepackage{booktabs} % nice rules (thick lines) for tables
\usepackage{array}
\usepackage{microtype} % improves typography for PDF

%\usepackage[hidelinks]{hyperref}
%\usepackage{caption}
\usepackage{subcaption}
\usepackage{hhline}
%\usepackage{amsmath}
%\usepackage{amssymb}
\usepackage{mathtools}
\allowdisplaybreaks
\usepackage{color}
\usepackage{multirow}
\usepackage{siunitx}
\usepackage{xfrac}
\usepackage{bm}
\usepackage{stmaryrd}

\usepackage{threeparttable, tablefootnote}

\usepackage{environ}
\makeatletter

\usepackage{tabularx}
\usepackage{float}
\usepackage{enumitem}
\setlist[itemize]{noitemsep}
\usepackage{diagbox}
\usepackage{courier}
\usepackage{pdflscape}

%\usepackage{cleveref}
\usepackage{datatool}
% \usepackage[numbers]{natbib}
\usepackage{notoccite}

\usepackage{nomencl} % If needed
\makenomenclature

\usepackage{tikz}
\usepackage{tkz-euclide}
\usetikzlibrary{positioning, arrows, decorations, shapes}
\usetikzlibrary{shapes.geometric,arrows}
\definecolor{illiniblue}{HTML}{B1C6E2}
\definecolor{illiniorange}{HTML}{f8c2a2}
\definecolor{green}{HTML}{c2e2b1}
\tikzstyle{process} = [rectangle, rounded corners, thick, minimum width=7cm, minimum height=1cm, text centered, draw=black, fill=illiniblue, text width=18em]
\tikzstyle{object} = [ellipse, minimum width=2cm, thick, minimum height=2.2cm, text centered, draw=black, fill=illiniorange, text width=12em]
\tikzstyle{decision} = [diamond, thick, aspect=2, minimum width=2cm, minimum height=2cm, text centered, draw=black, fill=green, text width=6em]
\tikzstyle{arrow} = [thick,->,>=stealth]

\title{A Hybrid $S_N$-Diffusion Method for Molten Salt Reactor Control Rod Modeling}

%%% Authors (use arabic numbers: 1, 2, 3, etc. for affiliationNumber)
%%% \addAuthor{GivenName MiddleInitial. FamilyName}{affiliationNumber}
\addAuthor[smpark3@illinois.edu]{Sun Myung Park}{1}
\addAuthor{Kathryn D. Huff}{1}
\addAuthor{Madicken Munk}{1,2}
% To move the * for the corresponding author move the Email address for primary/corresponding author only. Move the * next to the appropriate name. Do not use all capital letters for any part of the author's name] and insert it before the Author addition ex. if corresponding author is A then adjust \addAuthor[email]{Author A}{1}

%%% Affiliations (from authblk)
%%% \addAffiliation{affiliationNumber}{Name of Institute, City, State/Country}
\addAffiliation{1}{University of Illinois Urbana-Champaign, Urbana, IL}
\addAffiliation{2}{Oregon State University, Corvallis, OR}

%%% Write text for abstract
%%% Most text modifying commands will work in abstract
\Abstract{%
This paper presents a hybrid $S_N$-diffusion method for accurate modeling of control rods in
time-dependent molten salt reactor simulations. Neutron diffusion methods struggle with strong neutron
absorptivity of control rods, while neutron transport methods are computationally prohibitive
for time-dependent reactor analyses. The proposed hybrid method combines the strengths of both
approaches by applying transport corrections generated using the $S_N$ method near control rods.
Six 1-D test cases with increasing complexity are used to verify the implementation and assess
performance. Results show that the hybrid method provides improved rod worth estimates compared to
pure diffusion models. The hybrid method iterative scheme demonstrates robustness in convergence
and converges rapidly within two outer iterations. Additionally, the method accurately reproduces
much of the reference transport correction distribution from a full $S_N$ calculation despite the
$P_1$ approximation for $S_N$-diffusion boundary coupling. Overall, the hybrid $S_N$-diffusion
method offers a promising avenue for accurate control rod modeling in time-dependent reactor
analyses. Further work is ongoing to extend the hybrid method to higher dimensional models and
time-dependent analyses.
}

%%% List up to 5 keywords separated by a comma
\keywords{control rod modeling, molten salt reactor, neutron transport, reactor physics}

%%% Provide a short title for the header on odd pages
\shortTitle{A Hybrid $S_N$-Diffusion Method for Molten Salt Reactor Control Rod Modeling}

%%% Provide a short author listing for the header on even pages
\authorHead{S. Park, et al.}

%%% If LaTeX reports the line number of an error at \begin{document} it
%%%   is most likely due to an error in one of the commands above
\begin{document}

\section{Introduction}\label{sec:1}

\glspl*{MSR} are advanced reactors noted for strong passive safety due to their liquid fuel form,
negative temperature reactivity feedback, and low operating pressure. While some \gls*{MSR} concepts
completely eliminated control rods from their designs, many still retain them.

Routine time-dependent and spatially-resolved reactor analyses remain in the domain of neutron
diffusion methods due to the large computational expense of neutron transport methods. Several
transport correction techniques exist in the literature for improving neutron flux and
$k_\text{eff}$ estimates of neutron diffusion methods augmented by correction parameters in
various forms \cite{kavenoky_sph_1978, koebke_new_1980, gross_comprehensive_2023,
tamang_multilevel_2014, stehle_hybrid_2014}.

In this paper, we propose a hybrid $S_N$-diffusion neutronics method to improve control rod
modeling in neutron diffusion solvers. This iterative method adaptively applies the $S_N$ discrete
ordinates neutron transport method in subregions containing the control rod to obtain pointwise
transport corrections for the diffusion method.
This new method is unique in its approach towards adaptive domain decomposition into transport and
diffusion problem subdomains through on-the-fly evaluations of the transport correction
parameter distributions.

\section{Theory} \label{sec:theory}

\subsection{Multigroup \Gls{SAAF} $S_N$ Neutron Transport Method} \label{sec:saaf}

We implemented the hybrid method using finite element method numerical solver capabilities available through
Moltres \cite{lindsay_introduction_2018} and the \gls*{MOOSE} framework \cite{giudicelli_30_2024}.
The derivation for the \gls*{SAAF} formulation of the $S_N$ method and its
implementation in Moltres follows closely the derivation developed by Wang et al.\
\cite{wang_diffusion_2014} in Griffin (formerly Rattlesnake). This formulation reformulates
the first-order $S_N$ equations as second-order equations, enabling the use of fast and scalable
HYPRE multigrid preconditioners \cite{hypre_hypre_2022}.

The time-independent, multigroup neutron transport equations and boundary conditions defined on the
3-D spatial domain $\mathcal{D}$ and 2-D unit sphere angular domain $\mathcal{S}$ is:
%
\begin{multline}
  \hat{\Omega}\cdot\nabla\Psi_g(\vec{r},\hat{\Omega},t) + \Sigma_{t,g}
  \Psi_g(\vec{r},\hat{\Omega},t) = \\
  \sum^G_{g'=1}\int_\mathcal{S} \Sigma_s^{g'\rightarrow g}(\hat{\Omega}'\rightarrow\hat{\Omega})
  \Psi_{g'}(\vec{r},\hat{\Omega}',t)d\hat{\Omega}'
  + \frac{1}{4\pi}\frac{\chi_{g}}{k}\sum^G_{g'=1} \nu\Sigma_{f,g'} \phi_{g'}(\vec{r},t)
  \label{eq:mg-nte}
\end{multline}
%
\begin{gather}
  \Psi_g(\vec{r},\hat{\Omega}) = \Psi^\text{inc}_g(\vec{r},\hat{\Omega}) +
  \alpha^s_g\Psi_g(\vec{r},\hat{\Omega}_r)
  \mbox{ on } \vec{r} \in \partial\mathcal{D} \mbox{ and } \hat{\Omega}\cdot\hat{n}_b < 0.
  \label{eq:mg-nte-bc}
\end{gather}
%
We define inner products involving volume and surface integrals over the spatial domain for
expressing the weak formulation of multigroup $S_N$ equations for finite element analysis:
%
\begin{gather}
  \left(\Psi^*, \Psi\right)_\mathcal{D} \equiv \int_\mathcal{D}d\vec{r}
  \sum_{i\in E_e}\Psi^*_i b_i(\vec{r})\sum_{j\in E_e}\Psi_jb_j(\vec{r}), \\
  \left(\Psi^*, \Psi\right)_{\partial\mathcal{D}} \equiv
  \sum_{s\in\partial\mathcal{D}}\int_s d\vec{r}\sum_{i\in E_s}\Psi^*_i b_i(\vec{r})\sum_{j\in E_s}
  \Psi_j b_j(\vec{r}),
\end{gather}
%
where $b_i$ is a set of shape functions such that any function can be approximated as:
%
\begin{gather}
  f(\vec{r}) = \sum^{N_b}_{k=1}f_k b_k(\vec{r}),
\end{gather}
%
with $N_b$ being the dimension of the function space. $E_e$ and $E_s$ are index sets of non-zero
shape functions within each element $e$ or surface $s$.

In addition to the \gls*{SAAF} formulation, a void treatment scheme \cite{wang_diffusion_2014}
provides numerical stabilization in near-void regions where the $1/\Sigma_{t,g}$ factors in the
streaming and source terms evaluate to very large values. The final weak formulations of the
individual terms in the \gls*{SAAF} multigroup $S_N$ equations are given as:
%
\begin{gather}
  \shortintertext{Streaming term:}
%  \left(\mathbb{L}_1\bm{\Psi}^*,
%  \left(\bm{\tau}\mathbb{L}_1-\mathbb{I}+\bm{\tau}\mathbb{L}_2\right)\bm{\Psi}\right) =
  \sum^G_{g=1}\sum^{N_d}_{d=1}w_d\left(\hat{\Omega}_d\cdot\nabla\Psi^*_{g,d},\tau_g\hat{\Omega}
  \cdot\nabla\Psi_{g,d}-(1-\tau_g\Sigma_{t,g})\Psi_{g,d}\right)_\mathcal{D} \\
  \shortintertext{Collision term: }
%  \left(\mathbb{L}_2\bm{\Psi}^*,\bm{\Psi}\right) =
  \sum^G_{g=1}\sum^{N_d}_{d=1}w_d\left(\Psi^*_{g,d},\Sigma_{t,g}\Psi_{g,d}\right)_\mathcal{D} \\
  \shortintertext{Scattering term:}
%  \left(\left(\mathbb{I}+\bm{\tau}\mathbb{L}_1\right)\bm{\Psi}^*,\mathbb{S}\bm{\Psi}\right) =
  \sum^G_{g=1}\sum^{N_d}_{d=1}w_d\left(\Psi^*_{g,d}+\tau_g\hat{\Omega}_d\cdot\nabla\Psi^*_{g,d},
  \sum^G_{g'=1}\sum^L_{l=0}\Sigma^{g'\rightarrow g}_{s,l}\sum^l_{m=-l}
  \frac{2l+1}{w}Y_{l,m}(\hat{\Omega}_d)\phi_{g',l,m}\right)_\mathcal{D} \\
  \shortintertext{Fission source term: }
%  \left(\left(\mathbb{I}+\bm{\tau}\mathbb{L}_1\right)\bm{\Psi}^*,\mathbb{F}\bm{\Psi}\right) =
  \sum^G_{g=1}\sum^{N_d}_{d=1}w_d\left(\Psi^*_{g,d}+\tau_g\hat{\Omega}_d\cdot\nabla\Psi^*_{g,d},
  \frac{1}{w}\frac{\chi_g}{k}\sum^G_{g'=1}\nu\Sigma_{f,g'}\phi_{g'}\right)_\mathcal{D}
\end{gather}
%
$\tau_g$ is a numerical stabilization parameter borne from the void treatment scheme. $w_d$ are
weight values from the angular discretization scheme such that $w=\sum^{N_d}_{d=1} w_d$.
The boundary source and reflecting boundary terms are given as:
%
\begin{gather}
%  \langle\bm{\Psi}^*,\bm{\Psi}\rangle^+ - \langle\bm{\Psi}^*,\bm{\Psi}^\text{inc}\rangle^- =
  \mbox{Boundary source term: }
  \begin{cases}
    \sum^G_{g=1}\sum^{N_d}_{d=1}w_d\left(\Psi^*_{g,d},
    \hat{\Omega}_d\cdot\hat{n}_b\Psi_{g,d}\right)_{\partial\mathcal{D}},
    & \hat{\Omega}\cdot\hat{n}_b>0,\vec{r}\in\partial\mathcal{D} \\
    \sum^G_{g=1}\sum^{N_d}_{d=1}w_d\left(\Psi^*_{g,d},
    \hat{\Omega}_d\cdot\hat{n}_b\Psi^\text{inc}_{g,d}\right)_{\partial\mathcal{D}},
    & \hat{\Omega}\cdot\hat{n}_b<0,\vec{r}\in\partial\mathcal{D}
  \end{cases}, \label{eq:boundary-source} \\
%  \langle\bm{\Psi}^*,\bm{\Psi}\rangle^+ - \langle\bm{\Psi}^*,\mathbb{B}\bm{\Psi}\rangle^- =
  \mbox{Reflecting boundary term: }
  \begin{cases}
    \sum^G_{g=1}\sum^{N_d}_{d=1}w_d\left(\Psi^*_{g,d},
    \hat{\Omega}_d\cdot\hat{n}_b\Psi_{g,d}\right),
    & \hat{\Omega}\cdot\hat{n}_b>0,\vec{r}\in\partial\mathcal{D}_s \\
    \sum^G_{g=1}\sum^{N_d}_{d=1}w_d\left(\Psi^*_{g,d},
    \hat{\Omega}_d\cdot\hat{n}_b\Psi_{g,d_r}\right),
    & \hat{\Omega}\cdot\hat{n}_b<0,\vec{r}\in\partial\mathcal{D}_s
  \end{cases}, \label{eq:reflecting-bc}
\end{gather}
%
where $\Psi^\text{inc}_{g,d}=0$ for vacuum boundaries and
$\hat{\Omega}_{d_r} = \hat{\Omega}_d - 2(\hat{\Omega}_d\cdot\hat{n}_b)\hat{n}_b$ is the neutron
reflection direction.

\subsection{Transport-Corrected Multigroup Neutron Diffusion Method} \label{sec:transport-correction}

For transport corrections to the multigroup neutron diffusion equations, we adopt drift correction
terms which are expressed as $\nabla\cdot \vec{D}_g\phi_g$. Drift
coefficient values correct for the inaccuracy of the diffusion approximation on the streaming term
and the omission of anisotropic effects in neutron diffusion theory.
Derivations for the drift terms depend on the discretization schemes of the high- and low-order
equations. For this work, we adopted drift terms derived by Wang et al.\
\cite{wang_diffusion_2014}. The neutron balance equation derived from
integrating the \gls*{SAAF} $S_N$ equation over the angular variable is:
%
\begin{multline}
  \left(\bm{\Phi}^*,\frac{\partial}{\partial t}\left(\frac{\bm{\Phi}}{\bm{v}}\right)\right)_\mathcal{D}
  + \left(\nabla\bm{\Phi}^*, \mathbb{D}\nabla\bm{\Phi}\right)_\mathcal{D}
  + \left(\mathbb{L}_2\bm{\Phi}^*,\bm{\Phi}\right)_\mathcal{D}
  + \left(\bm{\Phi}^*,\frac{\bm{\Phi}}{2}\right)_{\mathcal{D}_v}
  + \left(\nabla\bm{\Phi}^*,\bm{\tau}\vec{\bm{r}}_1-\vec{\bm{J}}-\mathcal{D}\bm{\Phi}\right)_\mathcal{D}
  + \\
  \left(\bm{\Phi}^*,\bm{J}^\text{out}-\frac{\bm{\Phi}}{2}\right)_{\partial\mathcal{D}_v}
  = \left(\bm{\Phi}^*,\mathbb{S}_0\bm{\Phi}\right)_\mathcal{D}
  + \left(\bm{\Phi}^*,\bm{Q}_0\right)_\mathcal{D}. \label{eq:modified-diff}
\end{multline}
%
We can interpret Eq.\ \ref{eq:modified-diff} as a modified neutron diffusion equation with
transport corrections provided by a drift term with the drift vector defined as
%
\begin{gather}
  \vec{\bm{D}} \equiv \frac{\bm{\tau}\vec{\bm{r}}_1-\vec{\bm{J}}-\mathbb{D}\nabla\bm{\Phi}}{\bm{\Phi}}
\end{gather}
%
and a vacuum boundary correction term with the boundary coefficient vector defined as
%
\begin{gather}
  \bm{\gamma} \equiv \frac{\bm{J}^\text{out}}{\bm{\Phi}}-\frac{1}{2}\bm{I}.
\end{gather}
%
With $S_N$ angular discretization, the drift and boundary correction vector components are
evaluated as
%
\begin{gather}
  \vec{D}_g = \frac{\sum^{N_d}_{d=1}w_d\left(\tau_g\hat{\Omega}_d\hat{\Omega}_d\cdot\nabla\Psi_{g,d}
  + \left(\tau_g\Sigma_{t,g}-1\right)\hat{\Omega}_d\Psi_{g,d}
  - \tau_g\sum^G_{g'=1}\Sigma^{g'\rightarrow g}_{s,1}\hat{\Omega}_d\Psi_{g',d}
  - D_g\nabla\Psi_{g,d}\right)}{\sum^{N_d}_{d=1}w_d\Psi_{g,d}}, \label{eq:drift} \\
  \gamma_g =
  \frac{\sum_{\hat{\Omega}_d\cdot\hat{n}_b > 0}w_d |\hat{\Omega}_d\cdot\hat{n}_b |
  \Psi_{g,d}}{\sum^{N_d}_{d=1}w_d\Psi_{g,d}}. \label{eq:bound-coef}
\end{gather}

\subsection{$S_N$-Diffusion Iteration Algorithm \& Boundary Coupling} \label{sec:hybrid-algorithm}

To reduce the computational cost of the high-level $S_N$ calculation in a reactor
simulation, we propose reducing the problem domain of the $S_N$ method to a
\textit{correction subregion} containing the control rod
and its vicinity. Consequently, the hybrid $S_N$-diffusion method may retain accurate neutron flux
and current estimates around the control rod region from the $S_N$ method while making significant
computational cost savings by treating most of the reactor geometry with the neutron diffusion
method alone. Henceforth, we will refer to the $S_N$ calculation in the correction
region as the $S_N$ \textit{subproblem} or \textit{subsolver}. We define the full problem
domain and the correction subregion as $V_0$ and $V_1$, respectively, where
$V_1\subseteq V_0$. The iterative algorithm for the hybrid $S_N$-diffusion method is as follows:
%
\begin{enumerate}
  \item Start with an initial neutron diffusion calculation in $V_0$ using the standard neutron
    diffusion method.
  \item Pass the neutron diffusion current estimates along
    $\partial V_1$ to the $S_N$ subsolver to evaluate the boundary conditions for the $S_N$
    subproblem.
  \item With the $S_N$ subsolver, evaluate transport correction terms in $V_1$ using Eq.\
    \ref{eq:drift}.
  \item Pass the transport correction terms to the neutron diffusion solver to apply corrections
    within $V_1$ while continuing to apply the standard neutron diffusion solver
    in the rest of $V_0$.
  \item Start the next iteration by running a neutron diffusion calculation with transport
    corrections in $V_1$.
  \item Repeat Steps 2-6 until convergence is reached by meeting desired convergence tolerance
    values.
\end{enumerate}

For the hybrid $S_N$-diffusion method to converge, it requires appropriate boundary conditions for
the $S_N$ subproblem. Given that we want to limit the coverage of $V_1$ to the control rod region
and its immediate vicinity, $V_1$ should be sufficiently smaller than $V_0$, but large enough for
the boundary $\partial V_1$ to be far enough away from the rod that the $S_N$-diffusion interface
has physics that are favorable to diffusion. Therefore, the boundary $\partial V_1$ 
lies within $V_0$. We may obtain the best boundary source estimate for the $S_N$ subproblem by
applying the $P_1$ approximation for evaluating the neutron angular flux along
the discrete ordinates $\hat{\Omega}_d$ of the $S_N$ angular quadrature set as follows
%
\begin{align}
  \Psi_{g,d} \approx \frac{1}{4\pi}\left(\phi^\text{diff}_g+3\hat{\Omega}_d\cdot
  \vec{J}^\text{diff}_g\right)
  =\frac{1}{4\pi}\left(\phi^\text{diff}_g-3\hat{\Omega}_d\cdot D_g\nabla\phi^\text{diff}_g\right)
\end{align}
%
where the diff superscript denotes flux estimates from the latest neutron diffusion calculation.
With this relation, we can formulate the boundary source term for the $S_N$ subsolver by setting
$\Psi^\text{inc}_{g,d}$ in Eq.\ \ref{eq:boundary-source} to
%
\begin{gather}
  \Psi^\text{inc}_{g,d} = \frac{1}{w}
  \left(\phi^\text{diff}_g-3\hat{\Omega}_d\cdot D_g\nabla\phi^\text{diff}_g\right)
\end{gather}

The angular flux solution to the $S_N$ subproblem defined with these boundary sources typically
produces accurate drift correction parameters within $V_1$ but tend to be inaccurate
near $\partial V_1$. Consequently, we define a \textit{buffer region}
$V_1'$, within $V_1$ adjacent $\partial V_1$, where inaccurate drift correction parameters are
discarded. Right before each
neutron diffusion iteration, the solver adaptively determines a cutoff boundary for $V_1'$ which
determines the extent in $V_1$ where the drift correction terms are applied.
A natural/intuitive criterion for the location of the cutoff boundary
would be wherever the components of the drift correction variable $\vec{D}_g$ is zero (i.e.,
wherever the components change signs). The reasons for this criterion are twofold. Firstly, at
points where the $\vec{D}_g$ components are zero, the flux is approximately isotropic along the
axes corresponding to the components because value of $\vec{D}_g$ represents how much correction a
standard neutron diffusion model requires. Secondly, this choice preserves the smoothness of the
neutron scalar flux gradient.

\section{Model Setup}

We modeled the 1-D models in this work after the \gls*{MSRE} design \cite{robertson_msre_1965}. We
employed OpenMC \cite{romano_openmc:_2015} for generating group constants
for the hybrid method and reference solutions for verification. All multigroup
neutronics simulations ran with an eight neutron energy group structure developed by Jaradat
\cite{jaradat_development_2021-1} for \gls*{MSRE} analyses. The OpenMC models ran on continuous
energy (OpenMC-CE) and multigroup (OpenMC-MG) modes and utilized ENDF/B-VII.1.

This work evaluates six 1-D test cases with increasing complexity to verify the $S_N$ and hybrid
$S_N$-diffusion implementations and to test the performance of the hybrid
method in response to various geometrical features. The last two cases resemble the
reference \gls*{MSRE} design, which has centrally located control rods
and air-filled rod guide tubes. Figure \ref{fig:case-geom} shows the geometries of Cases 1a to
3b. All geometries have reflective boundary conditions at $x=0$ cm, thereby
forming half-core or repeating unit cell models. Cases 1a and 1b are repeating unit
cell models with reflecting boundaries on the right-side boundaries. Cases 2a, 2b, 3a, and 3b
are 1-D half-core models with vacuum boundaries on the right-side boundaries.

All material compositions and densities, except the control rod, follow
reference data at 900 K from \gls*{MSRE} benchmark specifications.
We reduced the Gd$_2$O$_3$ content in the control rod region from 70 wt\% to 0.35 wt\% to bring the
rod worth down from approximately 55000 pcm to 20000 pcm in the 1-D models.
%
\begin{figure}[htb!]
  \centering
  \includegraphics[width=0.75\columnwidth]{case-geometry}
  \caption{Geometries of the 1-D test cases. The material labeled ``mixture'' represents a
    homogeneous mixture of fuel and graphite at a ratio of 22.5\%-77.5\% by volume. All geometries
    have reflective boundary conditions on the boundary at $x=0$ cm. The right-side boundaries are
    reflective for Cases 1a and 1b, and vacuum for Cases 2a, 2b, 3a, and 3b.}
  \label{fig:case-geom}
\end{figure}

\section{Numerical Results \& Discussion}

The hybrid method shares some similarities with diffusion-based acceleration schemes for the
neutron transport methods which do not require the neutron transport side of the calculations to
fully converge \cite{wang_diffusion_2014}. The typical convergence threshold value required for
converged 1-D neutron flux calculations is $10^{-8}$. We varied the convergence tolerance value for
the $S_N$ subsolver for maintaining accuracy while minimizing computational costs.
Table \ref{table:sn-tol} lists the number of outer iterations required by the hybrid method for
a given set of $S_8$ subsolver convergence tolerance values in Case 3b.

\begin{table}[t]
  \begin{minipage}{0.49\columnwidth}
    \centering
    \small
    \caption{Number of outer iterations in hybrid method calculations of Case 3b for a given set of
    convergence tolerance values imposed on the $S_8$ subsolver.}
    \begin{tabular}{S S}
      \toprule
      {$S_8$ subsolver tolerance} & {No.\ of outer iterations} \\
      \midrule
      {$10^{-8}$} & 3 \\
      {$10^{-7}$} & 3 \\
      {$10^{-6}$} & 3 \\
      {$10^{-5}$} & 2 \\
      {$10^{-4}$} & 2 \\
      {$10^{-3}$} & 1 \\
      \bottomrule
    \end{tabular}
    \label{table:sn-tol}
  \end{minipage}
  \begin{minipage}{0.49\columnwidth}
    \centering
    \includegraphics[width=0.8\columnwidth]{sn-tol}
    \captionof{figure}{Error estimates of $k_\text{eff}$ in Case 3b as a function of
    convergence tolerance. $q=1.333$ represents the approximate rate of
    convergence.}
    \label{fig:sn-tol}
  \end{minipage}
\end{table}

Figure \ref{fig:sn-tol} shows the $k_\text{eff}$ error estimates relative to the reference value
when $\epsilon_\text{tol}=10^{-8}$. The hybrid method exhibits superlinear convergence ($q=1.333$)
with respect to the $S_N$ subsolver convergence tolerance value. Based on the results, setting
$\epsilon_\text{tol}$ to $10^{-5}$ is the best balance between accuracy and computational
cost because the $k_\text{eff}$ error is less than 0.01 pcm and the number of outer iterations
would increase from 2 to 3 if the tolerance value is further tightened.

Figure \ref{fig:1d-rho} shows the difference in reactivities of OpenMC-MG, $S_8$, neutron
diffusion, and hybrid methods relative to OpenMC-CE for all test cases. The standard deviations of
reactivity values from OpenMC-CE are approximately 40 pcm for Cases 1a and 1b and 60 pcm for the
rest as depicted by the blue highlighted area in Figure \ref{fig:1d-rho}. For Cases 1a and 1b, the
OpenMC-MG, $S_8$, and neutron diffusion methods show excellent agreement with OpenMC-CE as the
reactivity values fall within the one or two standard deviation range.
%
\begin{figure}[htb!]
  \centering
  \includegraphics[width=0.75\columnwidth]{rho}
  \caption{Difference in reactivity $\rho$ of all neutronics methods investigated relative
  to OpenMC-CE.}
  \label{fig:1d-rho}
  \includegraphics[width=0.75\columnwidth]{worth}
  \caption{Percentage difference in rod worth for Cases 2 and 3 of all neutronics methods
  investigated relative to OpenMC-CE.}
  \label{fig:1d-worth}
\end{figure}

For all cases, the OpenMC-MG and the $S_8$ methods show consistent agreement with one another. 
While they deviate from OpenMC-CE by approximately 350 pcm for Cases 2a and 3a, they fall within
two standard deviations of OpenMC-CE for Cases 2b and 3b. When increasing the maximum Legendre
polynomial order to approximate the angular dependence in $\Sigma_s^{g'\rightarrow g}$ from
2nd-order to 3rd-order, the reactivity changed by only 1 pcm. Therefore, we attribute the reactivity
discrepancies to the eight neutron energy group structure which remains the
only significant difference between OpenMC-CE and the multigroup methods.

The neutron diffusion and hybrid methods agree closely with one another for Cases 2a and 3a which
exclude the control rod region. This shows that the hybrid method provides similar $k_\text{eff}$
estimates to the neutron diffusion method in 1-D simulations without highly neutron-absorbing
regions. Compared to the OpenMC-MG and $S_8$ reactivity values, the neutron diffusion and hybrid
method reactivity values agree closer with the reference OpenMC-CE value. However, this is likely
due to favorable error cancellation since diffusion theory is an approximation to neutron
transport.

In Cases 2b and 3b, the neutron diffusion method largely fails to accurately capture the
effect of introducing the control rod region as evidenced by the -1500 pcm and -1150 pcm
discrepancies relative to OpenMC-CE. The hybrid method fares better with -400 pcm and -300 pcm
discrepancies.

Moving the discussion to control rod worths, Figure \ref{fig:1d-worth} shows the percentage
difference in rod worths for Cases 2 and 3 of OpenMC-MG, $S_8$, neutron diffusion, and hybrid
methods relative to OpenMC-CE. We calculated the rod worth estimates from each method by taking the
difference in reactivity between the ``rod out'' (a) and ``rod in'' (b) configurations of Cases 2
and 3. All neutronics methods overestimate the rod worth relative to OpenMC-CE, resulting in
positive percentage difference values observed in the figure.
The neutron diffusion methods clearly show up as outliers with rod worth estimates in excess
of 8\%. The remaining methods cluster around 2.5\% and 3\% percentage difference for Cases
2 and 3.

\begin{figure}[p]
  \centering
  \includegraphics[width=0.98\columnwidth]{case-3b-flux-diff}
  \caption{Absolute difference in neutron group flux distributions for Case 3b from $S_8$,
  neutron diffusion, and hybrid methods relative to OpenMC-MG.}
  \label{fig:3b-flux-diff}
\end{figure}

Given that the hybrid methods rely on transport corrections derived from the $S_N$ method, the
$S_8$ rod worth estimates serve as the
reference point for hybrid method verification. Additionally, errors arising from the multigroup
approximation affect both hybrid and $S_N$ simulations to similar degrees. The hybrid method
provides significant improvements in rod worth estimates over the neutron diffusion method.
%
\begin{figure}[t]
  \centering
  \includegraphics[width=\columnwidth]{case-3b-drift}
  \caption{Group 1 and 8 drift correction ($\vec{D}_g$) $x$-component distributions from the
  Moltres-hybrid and Moltres-$S_8$ solvers. They show nearly perfect overlap in the interior of the
  correction subregion ($x=0$ cm to $8$ cm).}
  \label{fig:3b-drift}
\end{figure}

Figure \ref{fig:3b-flux-diff} shows the absolute
difference in Case 3b neutron group flux distributions of the $S_8$, neutron diffusion, and
hybrid methods relative to OpenMC-MG, instead of OpenMC-CE, to eliminate discrepancies arising from
the multigroup approximation. Case 3b most closely resembles the actual \gls*{MSRE} design with an
inserted rod.
The neutron diffusion and hybrid methods fare worse than the $S_8$ method at capturing the
oscillatory pattern in groups 1, 2, 5, 7, and 8 arising from the fuel-graphite lattice geometry
along $x=5$ cm to 10 cm. The neutron diffusion method
exhibits significant flux deviations in groups 1, 2, 7, and 8 near $x=0$ cm. The hybrid method
provides significant improvements in neutron flux distributions compared to the neutron diffusion
method as shown by the generally smaller flux difference values particularly near $x=0$ cm.

Figure \ref{fig:3b-drift} shows the drift correction parameter distribution in neutron energy
group 1 and 8 for Case 3b. We computed the reference drift distributions from reference flux
solutions of the standard $S_8$ method.
The drift correction parameters generated by the correction scheme match their respective reference
values within the correction subregion from $x=0$ cm to 10 cm except near the subregion boundary at
$x=10$ cm as previously discussed in Section \ref{sec:hybrid-algorithm}. The adaptively defined
buffer region $V_1'$ completely encompasses the region near $x=10$ cm where the drift correction
distributions deviate from the reference distribution and are thus discarded.

\section{CONCLUSIONS}

This paper presents a hybrid $S_N$-diffusion method aimed at control rod modeling in time-dependent
simulations. The method involves generating drift correction parameters within a control rod and
its vicinity, using the $S_N$ method for the neutron diffusion method. The transport-corrected
subregion size adaptively changes in response to the drift correction distributions. We implemented
this method alongside the preexisting neutron diffusion solver in Moltres, a \gls*{MOOSE}-based
application. The reactivity and rod worth results in \gls*{MSR}-based 1-D test cases show the hybrid
method accurately reproduces rod worth estimates. The drift correction parameter
distributions indicate that the influence of the rod on the drift, and
consequently the flux, extend to approximately $x=5$ cm before the drift distribution settles on a
regular repeating pattern influenced by the fuel-graphite lattice structure. Further work is
ongoing for the extension of the hybrid method to 2-D and 3-D models and its application to
time-dependent simulations.

\section*{ACKNOWLEDGEMENTS}
Sun Myung Park was supported by a research assistantship provided by the Department of Nuclear,
Plasma, \& Radiological Engineering (NPRE) at the University of Illinois Urbana-Champaign (UIUC).

\bibliographystyle{mc2025}
\bibliography{bibliography}

\end{document}