moltres.tex

\section{Moltres} \label{sec:moltres}

In this section, we describe Moltres \cite{lindsay_introduction_2018}, the
multiphysics reactor simulation tool, and the specific modeling approach for
simulating the CNRS Benchmark cases in Moltres. Much of Moltres' development
focuses on meeting the needs of \gls{MSR} multiphysics.

\subsection{Description of Moltres} \label{sec:description-of-moltres}

Moltres models coupled neutronics and thermal-hydraulics in reactors. While
generally applicable to most reactor concepts, much of
Moltres' development focuses on meeting the needs of \gls{MSR} multiphysics.
Together with \gls{MOOSE}'s \cite{permann_moose_2020} \texttt{Heat}
\texttt{Conduction} and \texttt{Navier-Stokes} \cite{peterson_overview_2018}
modules, Moltres solves the multigroup neutron diffusion
equations, for an arbitrary number of energy and precursor groups, and
thermal-hydraulics equations simultaneously on the same mesh (or separately
through fixed-point iterations if desired).
Previously, Lindsay et al. \cite{lindsay_introduction_2018}
demonstrated Moltres' \gls{MSR} neutronics modeling capabilities with 1D salt
flow in 2D-axisymmetric and 3D models of the \gls{MSRE}. The neutron flux and
temperature distributions agreed qualitatively with legacy
\gls{MSRE} data, albeit with some minor quantitative discrepancies due to
simplifications and assumptions in the reactor geometry. Moltres has
since undergone further development to support 1) the looping of \gls{DNP}
drift back into the reactor core, 2) coupling the aforementioned \gls{DNP}
drift to incompressible Navier-Stokes velocity flows within the reactor core,
and 3) a decay heat model to simulate decay heat from fission products.

To perform neutronics calculations, Moltres requires homogenized group constant
data from dedicated high-fidelity neutronics software such as the NEWT module
in SCALE \cite{dehart_reactor_2011}, Serpent 2 \cite{leppanen_serpent_2014}, or
OpenMC \cite{romano_openmc:_2015}. Users can run a Python script in Moltres'
Github repository, which automatically reads
user-provided SCALE/Serpent 2/OpenMC output data files and creates
Moltres-compatible JSON or text files containing all required group constant
data.

We present the governing equations for the various physics models implemented
in Moltres. Moltres solves for the neutron fluxes governed by
the multigroup neutron diffusion equations given by:
%
\begin{align}
    \frac{1}{v_g} \frac{\partial \phi_g}{\partial t} =& \nabla \cdot D_g
    \nabla \phi_g - \Sigma^r_g \phi_g +
    \sum^G_{g' \neq g} \Sigma^s_{g' \rightarrow g} \phi_{g'} \nonumber \\
    &+ \chi^p_g \sum^G_{g'=1} \left( 1-\beta \right) \nu \Sigma^f_{g'}
    \phi_{g'} + \chi^d_g \sum^I_i \lambda_i C_i \label{eq:neut} \\
    %
    \shortintertext{where}
    v_g =& \text{ average speed of neutrons in group $g$,} 
    \nonumber \\
    \phi_g =& \text{ neutron flux in group $g$,}
    \nonumber \\
    t =& \text{ time,} \nonumber \\
    D_g =& \text{ diffusion coefficient of neutrons in} \nonumber \\
    &\text{ group $g$,} \nonumber \\
    \Sigma^r_g =& \text{ macroscopic cross section for removal of} \nonumber \\
    &\text{ neutrons from group $g$,} \nonumber \\
    \Sigma^s_{g' \rightarrow g} =& \text{ macroscopic cross section of
    scattering from} \nonumber \\
    &\text{ groups $g'$ to $g$,} \nonumber \\
    \chi^p_g =& \text{ prompt fission spectrum for neutrons in} \nonumber \\
    &\text{ group $g$,} \nonumber \\
    G =& \text{ total number of discrete neutron groups,} \nonumber \\
    \nu =& \text{ average number of neutrons produced per} \nonumber \\
    &\text{ fission,} \nonumber \\
    \Sigma^f_{g} =& \text{ macroscopic fission cross section for neutron}
    \nonumber \\
    &\text{ in group $g$,} \nonumber \\
    \chi^d_g =& \text{ delayed fission spectrum for neutrons in} \nonumber \\
    &\text{ group $g$,} \nonumber \\
    I =& \text{ total number of delayed neutron precursor} \nonumber \\
    &\text{ groups,} \nonumber \\
    \beta =& \text{ total delayed neutron fraction.} \nonumber
\end{align}

The delayed neutron precursor concentrations are
governed in Moltres by the following equation:
%
\begin{align}
    \frac{\partial C_i}{\partial t} =& \beta_i \sum^G_{g'=1} \nu \Sigma^f_{g'}
    \phi_{g'} - \lambda_i C_i - \vec{u} \cdot \nabla C_i + \nabla \cdot
    D_{\text{P}} \nabla C_i \label{eq:dnp} \\
    %
    \shortintertext{where}
    \beta_i =& \text{ delayed neutron fraction of precursor group $i$,}
    \nonumber \\
    \lambda_i =& \text{ average decay constant of delayed neutron} \nonumber \\
    &\text{ precursors in precursor group $i$,} \nonumber \\
    C_i =& \text{ concentration of delayed neutron precursors in}
    \nonumber \\
    &\text{ precursor group $i$,} \nonumber \\
    \vec{u} =& \text{ molten salt flow velocity vector,}
    \nonumber \\
    D_{\text{P}} =& \text{ diffusion coefficient of the delayed}
    \nonumber \\
    &\text{ neutron precursors.} \nonumber
\end{align}

The last two terms in Equation \ref{eq:dnp} represent the advection and
diffusion terms, respectively, to model the movement of \gls{DNP} in
liquid-fuel \glspl{MSR}.

The governing equation for temperature in Moltres is an advection-diffusion
equation with a fission heat source term $Q_f$ given by:
%
\begin{align}
    \rho c_{p} \frac{\partial T}{\partial t} =& - \rho c_p \vec{u}
    \cdot \nabla T + \nabla \cdot \left(k \nabla T \right) + Q_f
    \label{eq:temp}
    \shortintertext{and}
    Q_f =& \sum^G_{g=1} \epsilon_g \Sigma_g^f \phi_g
    \shortintertext{where}
    \rho =& \text{ density of the molten salt,}
    \nonumber \\
    c_p =& \text{ specific heat capacity of molten salt,} \nonumber \\
    T =& \text{ temperature of molten salt,} \nonumber \\
    k =& \text{ effective thermal conductivity of molten salt,} \nonumber \\
    Q_f =& \text{ fission heat source,} \nonumber \\
    \epsilon_g =& \text{ average fission energy released by neutrons in group
    $g$.} \nonumber
\end{align}

Lastly, the governing equations for the incompressible Navier-Stokes flow in
Moltres are given by:
%
\begin{align}
    \rho \frac{\partial \vec{u}}{\partial t} =&
    -\rho (\vec{u}
    \cdot \nabla) \vec{u} - \nabla p + \mu \nabla^2 \vec{u}
    + \rho \alpha \vec{g} \left(T - T_{\text{ref}} \right)
    \label{eq:momemtum}
    \shortintertext{and}
    \nabla \cdot \vec{u} =& 0
    \label{eq:divergence}
    \shortintertext{where}
    p =& \text{ pressure,} \nonumber \\
    \mu =& \text{ dynamic viscosity,} \nonumber \\
    \alpha =& \text{ coefficient of thermal expansion,} \nonumber \\
    \vec{g} =& \text{ gravitational force vector,} \nonumber
    \\
    T_{\text{ref}} =& \text{ reference temperature at which the nominal}
    \nonumber \\
    &\text{ density is provided.} \nonumber
    \nonumber
\end{align}

The velocity, temperature, and delayed neutron precursor variables are all
susceptible to numerical node-to-node oscillations near discontinuous boundary
conditions commonly observed when resolving advection-dominated transport using
continuous \gls{FEM} \cite{kuhlmann_lid-driven_2018}. MOOSE's
\texttt{Navier-Stokes} module provides the \gls{SUPG} stabilization scheme
\cite{brooks_streamline_1982} for the velocity and temperature variables to
combat these oscillations. We refer readers to \cite{peterson_overview_2018}
for details on the implementation of these methods in the
\texttt{Navier-Stokes} module. For the delayed neutron precursor variables, we
discretized them using discontinuous shape functions supported by \gls{MOOSE}'s
\gls{DGFEM} solver. The cell-centered \gls{DGFEM} scheme does not produce the
aforementioned numerical oscillations observed with node-centered continuous
\gls{FEM}.

\subsection[Modeling approach]{Modeling approach\footnote{The input files for
all benchmark
cases are available on the Moltres GitHub repository at 
\url{https://github.com/arfc/moltres/tree/devel/problems/2021-cnrs-benchmark}.
}} \label{sec:model}

For this work, we ran the benchmark cases on a uniformly-spaced mesh consisting
of 200$\times$200 elements (0.01m$\times$0.01m each). This mesh resolution was
sufficient for
mesh convergence as further refinement to 0.005m$\times$0.005m produced a
small 0.1 pcm increase in the reactivity for Step 0.2. We adopted the group
constant data provided by Tiberga et al. \cite{tiberga_results_2020}. Next, we
discretized most of the relevant variables, i.e., neutron fluxes, velocity
components, pressure, and temperature, using continuous, first-order Lagrange
shape functions. The only exceptions are the precursor concentration variables,
which we discretized using piecewise constant shape functions for the
\gls{DGFEM} solver mentioned in Section 3.1. In terms of solver accuracy,
\gls{DGFEM} with piecewise constant discretization is similar to first-order
\gls{FVM} because they share the same
number of degrees of freedom. We interpolated the resulting discontinuous,
cell-centered precursor values to obtain the nodal values for results
analysis. Given the high Schmidt number ($2\times10^8$)
\cite{tiberga_results_2020} of the salt, we neglected the precursor diffusion
term in Equation \ref{eq:dnp} as it has no observable effect on the
distribution.

As mentioned in Section \ref{sec:description-of-moltres}, the
\texttt{Navier-Stokes} and \texttt{Heat} \texttt{Conduction} modules from
\gls{MOOSE} provide some of the capabilities for
modeling incompressible flow and heat transfer. In particular, we stabilized
the incompressible flow and temperature governing equations using the
\gls{SUPG} stabilization method implemented in \gls{MOOSE}
\cite{peterson_overview_2018}. Without \gls{SUPG} stabilization, we
observed spurious numerical oscillations in the velocity and temperature near
the top boundary due to the singularity on the top left corner where different
velocity boundary conditions meet. We also applied the \gls{PSPG} stabilization
scheme \cite{hughes_new_1986} from the Navier-Stokes module
\cite{peterson_overview_2018}
which enables equal-order discretizations in the velocity and pressure
variables. Equal-order discretizations with \gls{PSPG} are computationally
cheaper and more convenient to work with than implementing higher-order
velocity discretizations for stability without \gls{PSPG}
\cite{chapelle_inf-sup_1993}.

We performed all neutronics criticality calculations in Steps 0.2, 1.1, 1.2,
1.3, and 1.4 using the inverse power method solver in \gls{MOOSE} and other
calculations in Steps 0.1, 0.3, and 2.1 using the Preconditioned Newton-Krylov
solver \cite{gaston_physics-based_2015}. The coupled steady-state problems in
Steps 1.2, 1.3, and 1.4 required segregated solvers due to the unique problem
setups involving time-independent neutronics criticality calculations
and pseudo-transient thermal-hydraulics calculations.

\begin{table}[tb]
    \caption{Timestep sizes used for the time-dependent cases in
    Step 2.1, corresponding to 1/200th of the perturbation period.}
	\centering
	\setlength\tabcolsep{2.5pt}
	\begin{tabular}{l l l l l l l l}
	    \toprule
	    Frequency [Hz] & 0.0125 & 0.025 & 0.05 & 0.1 & 0.2 & 0.4 & 0.8 \\
	    \midrule
	    Timestep size [s] & 0.2 & 0.2 & 0.1 & 0.05 & 0.025 & 0.0125 & 0.00625
	    \\
	    \bottomrule
	\end{tabular}
	\label{table:timestep}
\end{table}

For the time-dependent cases in Step 2.1, we employed fully coupled solves with
a second-order implicit Backward Differential Formula (BDF2) time-stepping
scheme. We set the timestep sizes to 1/200th of the perturbation period for
each driving frequency in the heat transfer coefficient. Table
\ref{table:timestep} shows the timestep sizes. We assumed the
systems reached asymptotic behavior when the magnitudes of neighboring power
peaks differed by less than 0.001\% for at least ten wavelengths. Under this
assumption, the phase shift measurements between adjacent waves always
converged before the magnitude measurements of the power peaks.

Table \ref{table:software} compares the numerical methods, meshing schemes, and
neutronics models of Moltres and the four participating software packages in
the CNRS benchmark paper \cite{tiberga_results_2020}. The $SP_N$ and
$S_N$ neutronics models refer to the simplified $P_N$ spherical harmonics and
$S_N$ discrete ordinates neutron transport models, respectively. Based on the
solvers and methods of solution, Moltres is most similar to the
PHANTOM-$S_N$ + DGFlows \cite{tiberga_discontinuous_2019} multiphysics package
from \gls{TUD} with the $S_2$ neutron transport model. Participants from
\gls{CNRS} and \gls{PSI} employed non-uniform meshes with smaller mesh elements
near the boundaries, while we and
the participants from \gls{PoliMi} and \gls{TUD} employed uniform meshes.

\FloatBarrier

\begin{landscape}
\begin{table*}[p]
    \caption{List of software packages and their corresponding model
    specifications for the CNRS Benchmark simulations
    \cite{tiberga_results_2020}.}
    \centering
    \begin{tabular}{p{4.2cm} p{7cm} p{3.3cm} p{2cm} p{2.7cm}}
        \toprule
        Software & Institute & Numerical method & Mesh & Neutronics model \\
        \midrule
        OpenFOAM & Centre national de la recherche scientifique (CNRS) & Finite volume & 200$\times$200 \newline Non-uniform & $SP_1$ \& $SP_3$ \\
        OpenFOAM & Politecnico di Milano (PoliMi) & Finite volume & 400$\times$400 \newline Uniform & Neutron diffusion \\
        GeN-Foam & Paul Scherrer Institute (PSI) & Finite volume & 200$\times$200 \newline Non-uniform & Neutron diffusion \\
        PHANTOM-$S_N$+DGFlows & Delft University of Technology (TUD) & Discontinuous finite \newline element & 50$\times$50 \newline Uniform & $S_2$ \& $S_6$ \\
        Moltres (This work) & University of Illinois at Urbana-Champaign (UIUC) & Continuous \& discontinuous finite element & 200$\times$200 \newline Uniform & Neutron diffusion \\
        \bottomrule
    \end{tabular}
    \label{table:software}
\end{table*}
\end{landscape}

\FloatBarrier