HPMR-H $var element = document.getElementById("moose-equation-cd83828e-bef4-4fe7-9cf1-9c7e15d27cab");katex.render("_2", element, {displayMode:false,throwOnError:false});$ Reactor Model

Contact: Stefano Terlizzi sbt5572.at.psu.edu, Javier Ortensi javier.ortensi.at.inl.gov

Model link: Direwolf Steady State Model

Mesh

The reactor module in MOOSE (Shemon et al., 2022) was used to generate three different meshes for: (a) the Discontinuous Finite Element (DFEM) discrete ordinates (SN) neutronic solver shown in Figure 1 (a), one for the Bison model, shown in Figure 1 (b), and one for the Coarse Mesh Finite Difference (CMFD) acceleration mesh. displayed in Figure 1 (c). A one-twelfth radially reflected geometry was built by exploiting the problem symmetry to minimize the number of degrees of freedom and increase the speed of calculations (Terlizzi and Labouré, 2023).

Figure 1: (a) DFEM-SN mesh, (b) Bison mesh, and (c) CMFD mesh (Terlizzi and Labouré, 2023).

Cross Sections

Serpent (v. 2.1.32) was used to generate the multigroup cross sections for the HPMR-H $var element = document.getElementById("moose-equation-ad849900-b8f5-4934-ae1f-12e6e1ec6ce4");katex.render("_2", element, {displayMode:false,throwOnError:false});$ problem. The ENDF/B-VIII.0 continuous energy library was utilized to leverage the $var element = document.getElementById("moose-equation-cfee6a48-41d7-4b1d-a7e7-bb053eec64df");katex.render("YH_x", element, {displayMode:false,throwOnError:false});$ scattering libraries, which was then converted into an 11-group structure to perform the calculations. The conversion is performed via ISOXML (which is contained within Griffin) by reading the Serpent tallies and converting them into a multigroup XS library readable by Griffin. The group upper boundaries are reported in the Table 1.

Table 1: Energy group (upper) boundaries for the 11-group structure used in the Griffin model (Terlizzi and Labouré, 2023).

Group	Energy (MeV)	Group	Energy (MeV)
1	4.00E+01	7	9.88E-06
2	8.21E-01	8	4.00E-06
3	1.83E-01	9	1.00E-06
4	4.90E-02	10	3.20E-07
5	4.54E-04	11	6.70E-08
6	4.81E-05

The boundaries were set to capture the main resonances, and the thermal peak in the spectra is shown in Figure 2 (Terlizzi and Labouré, 2023).

Figure 2: Spectrum per unit lethargy for the first ring of assemblies, second ring of assemblies, full reactor, and group boundaries (dotted vertical gray lines) (Terlizzi and Labouré, 2023).

The cross sections were parametrized with respect to four variables: (1) fuel temperature, $var element = document.getElementById("moose-equation-19eeb5f7-d82f-4569-aaca-fbedb2975eae");katex.render("T_f", element, {displayMode:false,throwOnError:false});$ ; (2) moderator, monolith, and HP temperature, $var element = document.getElementById("moose-equation-e519d036-a82a-4be9-94ee-e6253f0b127b");katex.render("T_m", element, {displayMode:false,throwOnError:false});$ ; (3) reflector temperature, $var element = document.getElementById("moose-equation-49ae06e8-c20d-4b70-af81-e47d66c15aac");katex.render("T_r", element, {displayMode:false,throwOnError:false});$ ; and (4) the hydrogen-yttrium stoichiometric ratio, $var element = document.getElementById("moose-equation-1f741bd8-d979-416c-81fa-eb1f20a7a837");katex.render("c_H", element, {displayMode:false,throwOnError:false});$ . The choice of capturing the change in hydrogen content by tabulating the cross sections based on the local value of the stoichiometric ratio was based on the analysis performed in (Ni et al., 2021). The state points at which the multigroup cross sections were computed are as follows in Table 2 (Terlizzi and Labouré, 2023):

Table 2: State points at which the multigroup cross sections were computed (Terlizzi and Labouré, 2023).

Variable	1	2	3
$var element = document.getElementById("moose-equation-5c140231-5f92-42a5-8bc6-c4b490ea5e56");katex.render("T_f", element, {displayMode:false,throwOnError:false});$ , K	800	1000	1200
$var element = document.getElementById("moose-equation-b7d7f923-3c68-4c97-885d-cf6af0c50223");katex.render("T_m", element, {displayMode:false,throwOnError:false});$ , K	800	1000	1200
$var element = document.getElementById("moose-equation-1fa166e3-8f1f-4fa5-9d44-0d59334385d4");katex.render("T_r", element, {displayMode:false,throwOnError:false});$ , K	800	1000	N/A
$var element = document.getElementById("moose-equation-1c6cc499-89ce-44e2-87d5-099147415bc3");katex.render("c_H", element, {displayMode:false,throwOnError:false});$	1.7	1.8	1.9

Griffin Input for Neutronics

Griffin is the parent application for this simulation and is the input associated to the neutronic simulation. The Griffin model uses the discontinuous finite element method (DFEM-SN) accelerated through the coarse mesh finite difference (CMFD) acceleration to solve the multigroup neutron transport equation (Wang et al., 2021). This is visible from the transport system block where the type of problem to be solved and the angular and energy approximation used are specified.

[TransportSystems]
  equation_type = eigenvalue
  particle = neutron
  G = 11

  ReflectingBoundary = 'bottom topleft'
  VacuumBoundary = 'right back front'

  [sn]
    scheme = DFEM-SN
    n_delay_groups = 6
    family = MONOMIAL
    order = FIRST
    AQtype = Gauss-Chebyshev
    NPolar = 1
    NAzmthl = 3
    NA = 1
    sweep_type = asynchronous_parallel_sweeper
    using_array_variable = true
    collapse_scattering = true
    hide_angular_flux = true
  []
[]

The fixed point logics based upon the effective multiplication factor convergence within 1 pcm is enforced through "fixed_point_solve_outer = true" and the imposition of the multiplication factor as a custom postprocessor in the executioner block.

[Executioner]
  type = SweepUpdate
  # verbose = true
  # debug_richardson = true
  # Richardson
  richardson_rel_tol = 1e-3
  richardson_max_its = 50
  inner_solve_type = GMRes
  max_inner_its = 2
  # CMFD solver
  cmfd_acceleration = true
  coarse_element_id = coarse_element_id
  diffusion_eigen_solver_type = newton
  prolongation_type = multiplicative
  max_diffusion_coefficient = 1
  # Fixed point iteration with MultiApps
  fixed_point_solve_outer = true
  fixed_point_algorithm = picard
  custom_pp = eigenvalue
  custom_rel_tol = 1e-50
  custom_abs_tol = 1e-5
  force_fixed_point_solve = true
  fixed_point_max_its = 6
[]

The full input is reported below for the sake of completeness.

# ##################################################################################################################
# AUTHORS:            Stefano Terlizzi and Vincent Labouré
# DATE:               June 2023
# ORGANIZATION:       Idaho National Laboratory (INL), C-110
# CITATION:           S. Terlizzi and V. Labouré. (2023).
#                     "Asymptotic hydrogen redistribution analysis in yttrium-hydride-moderated heat-pipe-cooled
#                     microreactors using DireWolf", Annals of Nuclear Energy, Vol. 186, 109735.
# DESCRIPTION:        Main input file containing neutronic model for the SiMBA reactor.
# FUNDS:              This work was supported by INL Laboratory Directed Research&
#                     Development (LDRD) Program under U.S. Department of Energy (DOE) Idaho Operations Office,
#                     United States of America contract DE-AC07-05ID14517.
# ##################################################################################################################
# Block assignment
bison_fuel_blocks = '1 2'
moderator_blocks = '3 4'
monolith_blocks = '8 22' # For now fill the air gap with monolith
reflector_blocks = '10 11 13 14 15 16 17 18 19'
absorber_blocks = '12'
bison_mod_blocks = '${moderator_blocks} ${monolith_blocks}'
bison_ref_blocks = '${reflector_blocks} ${absorber_blocks}'
[Mesh]
  [main]
    type = FileMeshGenerator
    file = 'one_twelfth_simba_core_in.e'
  []
  [coarse_mesh]
    type = FileMeshGenerator
    file = 'one_twelfth_simba_core_CM_in.e'
  []
  [mesh]
    type = CoarseMeshExtraElementIDGenerator
    input = main
    coarse_mesh = coarse_mesh
    extra_element_id_name = coarse_element_id
  []
[]

[PowerDensity]
  power = 1.66666e5 # power: 2e6 / 12 = 1.66666e5 W (per 1/12th core)
  power_density_variable = power_density
  integrated_power_postprocessor = integrated_power
[]

[TransportSystems]
  equation_type = eigenvalue
  particle = neutron
  G = 11

  ReflectingBoundary = 'bottom topleft'
  VacuumBoundary = 'right back front'

  [sn]
    scheme = DFEM-SN
    n_delay_groups = 6
    family = MONOMIAL
    order = FIRST
    AQtype = Gauss-Chebyshev
    NPolar = 1
    NAzmthl = 3
    NA = 1
    sweep_type = asynchronous_parallel_sweeper
    using_array_variable = true
    collapse_scattering = true
    hide_angular_flux = true
  []
[]

[GlobalParams]
  library_file = '../cross_sections_library/simba_fmh_b_fine.xml'
  library_name = simba_fmh_b_fine
  densities = 1.0
  isotopes = 'pseudo'
  dbgmat = false
  grid_names = 'Tfuel Tmod Trefl ch'
  grid_variables = 'Tfuel Tmod Trefl ch'
  is_meter = true
[]

[AuxVariables]
  [Tfuel] # fuel temperature
    initial_condition = 1000
  []
  [Tmod] # moderator + monolith + HP temperature
    initial_condition = 1000
  []
  [Trefl] # radial reflector temperature
    initial_condition = 900
  []
  [ch]
    initial_condition = 1.8
  []
[]

[Materials]
  [fuel]
    type = CoupledFeedbackNeutronicsMaterial
    block = '1 2' # fuel pin with 1 cm outer radius, no gap
    material_id = 1001
    plus = true
  []
  [moderator]
    type = CoupledFeedbackNeutronicsMaterial
    block = '3 4 5' # moderator pin with 0.975 cm outer radius
    material_id = 1002
  []
  [monolith]
    type = CoupledFeedbackNeutronicsMaterial
    block = '8'
    material_id = 1003
  []
  [hpipe]
    type = CoupledFeedbackNeutronicsMaterial
    block = '6 7' # gap homogenized with HP
    material_id = 1004
  []
  [be]
    type = CoupledFeedbackNeutronicsMaterial
    block = '10 11 13 14 15'
    material_id = 1005
  []
  [b4c]
    type = CoupledFeedbackNeutronicsMaterial
    block = '12'
    material_id = 1006
  []
  [air]
    type = CoupledFeedbackNeutronicsMaterial
    block = '20 21 22'
    material_id = 1007
  []
  [axial_refl]
    type = CoupledFeedbackNeutronicsMaterial
    block = '16 17 18 19'
    material_id = 1008
  []
[]

[Executioner]
  type = SweepUpdate
  # verbose = true
  # debug_richardson = true
  # Richardson
  richardson_rel_tol = 1e-3
  richardson_max_its = 50
  inner_solve_type = GMRes
  max_inner_its = 2
  # CMFD solver
  cmfd_acceleration = true
  coarse_element_id = coarse_element_id
  diffusion_eigen_solver_type = newton
  prolongation_type = multiplicative
  max_diffusion_coefficient = 1
  # Fixed point iteration with MultiApps
  fixed_point_solve_outer = true
  fixed_point_algorithm = picard
  custom_pp = eigenvalue
  custom_rel_tol = 1e-50
  custom_abs_tol = 1e-5
  force_fixed_point_solve = true
  fixed_point_max_its = 6
[]

[MultiApps]
  [bison]
    type = FullSolveMultiApp
    input_files = thermal_ss.i
    execute_on = 'timestep_end'
    keep_solution_during_restore = true # to restart from the latest solve of the multiapp (for pseudo-transient)
  []
[]

[Transfers]
  [pdens_to_modules]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    to_multi_app = bison
    source_variable = power_density
    variable = power_density
    from_blocks = '1 2'
    to_blocks = '${bison_fuel_blocks}'
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [tfuel_from_modules]
    type = MultiAppGeneralFieldNearestNodeTransfer
    from_multi_app = bison
    source_variable = Tsolid
    variable = Tfuel
    from_blocks = '${bison_fuel_blocks}'
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [tmod_from_modules]
    type = MultiAppGeneralFieldNearestNodeTransfer
    from_multi_app = bison
    source_variable = Tsolid
    variable = Tmod
    from_blocks = '${bison_mod_blocks}'
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [trefl_from_modules]
    type = MultiAppGeneralFieldNearestNodeTransfer
    from_multi_app = bison
    source_variable = Tsolid
    variable = Trefl
    from_blocks = '${bison_ref_blocks}'
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [ch_from_modules]
    type = MultiAppGeneralFieldNearestNodeTransfer
    from_multi_app = bison
    source_variable = ch
    variable = ch
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
[]

[Postprocessors]
  [fuel_temp_avg]
    type = ElementAverageValue
    variable = Tfuel
    block = '1 2'
  []
  [fuel_temp_max]
    type = ElementExtremeValue
    variable = Tfuel
    block = '1 2'
  []
  [fuel_temp_min]
    type = ElementExtremeValue
    variable = Tfuel
    value_type = min
    block = '1 2'
  []
  [mod_temp_avg]
    type = ElementAverageValue
    variable = Tmod
    block = '3 4'
  []
  [mod_temp_max]
    type = ElementExtremeValue
    variable = Tmod
    block = '3 4'
  []
  [mod_temp_min]
    type = ElementExtremeValue
    variable = Tmod
    value_type = min
    block = '3 4'
  []
  [fuel_volume]
    type = VolumePostprocessor
    block = '1 2'
  []
  [ch_avg]
    type = ElementAverageValue
    variable = ch
    block = '3 4'
  []
  [ch_max]
    type = ElementExtremeValue
    variable = ch
    block = '3 4'
  []
  [ch_min]
    type = ElementExtremeValue
    variable = ch
    value_type = min
    block = '3 4'
  []
[]

[Outputs]
  csv = true
  file_base = 'neutronics_eigenvalue'
  [console]
    type = Console
    outlier_variable_norms = false
  []
  [pgraph]
    type = PerfGraphOutput
    level = 2
  []
[]

Bison Input for Heat Conduction and Hydrogen Redistribution

The Bison input is used to solve the heat conduction problem in the solid and the hydrogen redistribution in the moderator. The hydrogen redistribution problem is modeled through the following equation for the hydrogen stoichiometric ratio, $var element = document.getElementById("moose-equation-f4ae46ad-2b1f-4b0a-ad2c-45ac318d7da2");katex.render("c_H", element, {displayMode:false,throwOnError:false});$ (Matthews et al., 2021):

$var element = document.getElementById("moose-equation-f64f678f-4837-4a14-a2f6-e544932d2f51");katex.render(" \\frac{\\partial {c_H}}{\\partial t} = \\nabla \\cdot [-D(\\nabla {c_H} + \\frac{Q {c_H}}{R T^2} \\nabla T)],", element, {displayMode:true,throwOnError:false});$

where D is the diffusion coefficient of hydrogen in $var element = document.getElementById("moose-equation-4e717ded-db3c-4880-b711-1d21341015c9");katex.render("YH_x", element, {displayMode:false,throwOnError:false});$ , R is the gas constant, T is the temperature, and Q is the heat of transport . In the Bison input, the diffusion coefficient is set to one in order to accelerate the convergence for the steady-state calculations given than the asymptotic hydrogen distribution does not depend on the diffusion coefficient (Terlizzi and Labouré, 2023). The Bison model is coupled to 101 heat pipes (HP) sub-applications that model the flow of sodium within the reactor. The 101 instances of HPs are generated in the correct position through the MultiApps system from the single heat pipe input described in the following section. The coupling between Bison and Sockeye is performed through convective boundary conditions.

The full input is reported below for the sake of completeness.

# ##################################################################################################################
# AUTHORS:            Stefano Terlizzi and Vincent Labouré
# DATE:               June 2023
# ORGANIZATION:       Idaho National Laboratory (INL), C-110
# CITATION:           S. Terlizzi and V. Labouré. (2023).
#                     "Asymptotic hydrogen redistribution analysis in yttrium-hydride-moderated heat-pipe-cooled
#                     microreactors using DireWolf", Annals of Nuclear Energy, Vol. 186, 109735.
# DESCRIPTION:        Model for heat conduction and hydrogen migration.
# FUNDS:              This work was supported by INL Laboratory Directed Research&
#                     Development (LDRD) Program under U.S. Department of Energy (DOE) Idaho Operations Office,
#                     United States of America contract DE-AC07-05ID14517.
# ##################################################################################################################
# Block assignment
fuel_blocks = '1 2'
moderator_blocks = '3 4'
monolith_blocks = '8 22' # For now fill the air gap with monolith
reflector_blocks = '10 11 13 14 15 16 17 18 19'
upper_reflector_blocks = '16 17'
absorber_blocks = '12'
# for radial reflector BC
rrefl_htc = 25 # W/m^2/K arbitrary (but small) HTC
rrefl_Tsink = 800 # K, arbitrary Tsink
T_init = 950
# Parameter to perform sensitivity analysis on Q, heat of transport
qvalue_multiplier = 1.0

[Mesh]
  [main]
    type = FileMeshGenerator
    file = 'one_twelfth_simba_core_in.e'
  []
  # remove blocks and add sidesets to apply boundary conditions
  [add_sideset_hp]
    type = SideSetsBetweenSubdomainsGenerator
    input = main
    primary_block = '8 17' # add 16 so the HP boundary extends into the upper axial reflector
    paired_block = '7'
    new_boundary = 'hp'
  []
  [add_sideset_inner_mod_gap]
    type = SideSetsBetweenSubdomainsGenerator
    input = add_sideset_hp
    primary_block = '4'
    paired_block = '5'
    new_boundary = 'gap_mod_inner'
  []
  [add_sideset_outer_mod_gap]
    type = SideSetsBetweenSubdomainsGenerator
    input = add_sideset_inner_mod_gap
    primary_block = '8'
    paired_block = '5'
    new_boundary = 'gap_mod_outer'
  []
  [add_sideset_central_hole]
    type = SideSetsBetweenSubdomainsGenerator
    input = add_sideset_outer_mod_gap
    primary_block = '22'
    paired_block = '21'
    new_boundary = 'central_hole'
  []
  [remove_blocks] # remove HPs, moderator gap and central shutdown cooling hole
    type = BlockDeletionGenerator
    input = add_sideset_central_hole
    block = '6 7 5 20 21'
  []
[]

[Variables]
  [Tsolid]
    initial_condition = ${T_init} # set all the temperatures to T_sink in Sockeye to allow a smooth ramp-up
  []
  [ch]
    initial_condition = 1.8
    block = ${moderator_blocks}
  []
[]

[Kernels]
  # Heat conduction
  [heat_time_derivative]
    type = HeatConductionTimeDerivative
    variable = Tsolid
  []
  [heat_conduction]
    type = HeatConduction
    variable = Tsolid
  []
  [heat_source_fuel]
    type = CoupledForce
    variable = Tsolid
    block = '${fuel_blocks}'
    v = power_density
  []
  # Hydrogen redistribution
  [ch_time_derivative]
    type = TimeDerivative
    variable = ch
    block = ${moderator_blocks}
  []
  [ch_dxdx]
    type = MatDiffusion
    variable = ch
    diffusivity = D
    block = ${moderator_blocks}
  []
  [soretDiff]
    type = ThermoDiffusion
    variable = ch
    temp = Tsolid
    heat_of_transport = Q
    mass_diffusivity = D
    block = ${moderator_blocks}
  []
[]

[AuxVariables]
  [power_density]
    block = '${fuel_blocks}'
    family = MONOMIAL
    order = FIRST
    initial_condition = 2302579.81614
  []
  # add upper_reflector_blocks blocks because HP cooling occurs there too
  [htc_hp_var]
    block = '${monolith_blocks} ${upper_reflector_blocks}'
    initial_condition = 372
  []
  [Tsink_hp_var]
    block = '${monolith_blocks} ${upper_reflector_blocks}'
    initial_condition = 900
  []
  [Tfuel_layered_average]
    order = CONSTANT
    family = MONOMIAL
    block = '${fuel_blocks}'
  []
[]

[BCs]
  [hp_bc]
    type = CoupledConvectiveHeatFluxBC
    variable = Tsolid
    boundary = 'hp' # inside of heat pipe
    htc = htc_hp_var # given by Sockeye
    T_infinity = Tsink_hp_var # given by Sockeye
  []
  [rrefl_bc]
    type = CoupledConvectiveHeatFluxBC
    variable = Tsolid
    boundary = 'right' # reflector cooling
    htc = ${rrefl_htc}
    T_infinity = ${rrefl_Tsink}
  []
[]

[ThermalContact]
  [moderator_gap]
    type = GapHeatTransfer
    emissivity_primary = 0.8
    emissivity_secondary = 0.8
    variable = Tsolid
    primary = 'gap_mod_inner'
    secondary = 'gap_mod_outer'
    gap_conductivity = 0.08 # W/m/K, typical thermal connectivity for air at 1000C
    quadrature = true
  []
[]

[Materials]
  #### DENSITY #####
  # units of kg/m^3
  [fuel_density]
    type = Density
    block = '${fuel_blocks}'
    density = 14.3e3 # same as in Serpent input
  []
  [moderator_density]
    type = Density
    block = '${moderator_blocks}'
    density = 4.3e3 # same as in Serpent input
  []
  [monolith_density]
    type = Density
    block = '${monolith_blocks}'
    density = 1.8e3 # same as in Serpent input
  []
  [reflector_density]
    type = Density
    block = '${reflector_blocks}'
    density = 1.85376e3 # same as in Serpent input
  []
  [absorber_density]
    type = Density
    block = '${absorber_blocks}'
    density = 2.52e3 # same as in Serpent input
  []
  ### THERMAL CONDUCTIVITY ###
  # units of W/m-K
  [fuel_kappa]
    type = HeatConductionMaterial
    block = '${fuel_blocks}'
    thermal_conductivity = 19 # W/m/K
    specific_heat = 300 # reasonable value
  []
  [moderator_kappa]
    type = HeatConductionMaterial
    block = '${moderator_blocks}'
    thermal_conductivity = 20 # W/m/K
    specific_heat = 500 # random value
  []
  [monolith_thermal_props]
    type = HeatConductionMaterial
    block = '${monolith_blocks}'
    thermal_conductivity = 70 # typical value for G348 graphite
    specific_heat = 1830 # typical value for G348 graphite
  []
  [reflector_kappa]
    type = HeatConductionMaterial
    block = '${reflector_blocks}'
    thermal_conductivity = 200 # W/m/K, typical value for Be
    specific_heat = 1825 # typical value for Be
  []
  [absorber_kappa]
    type = HeatConductionMaterial
    block = '${absorber_blocks}'
    thermal_conductivity = 20 # W/m/K, typical value for B4C
    specific_heat = 1000 # typical value for Be
  []
  # Material properties for migration. D can e adjusted for convergence in steady-state
  [moderator_diffusivity]
    type = GenericConstantMaterial
    block = ${moderator_blocks}
    prop_names = 'D Q'
    prop_values = '1 ${fparse 5.3e3 * qvalue_multiplier}'
  []
[]

[Functions]
[]

[Executioner]
  type = Transient
  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart'
  petsc_options_value = 'hypre boomeramg 300'
  line_search = 'none'
  automatic_scaling = true

  l_tol = 1e-02
  nl_abs_tol = 1e-7
  nl_rel_tol = 1e-8

  l_max_its = 50
  nl_max_its = 25

  start_time = 0
  end_time = 3e4

  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.5
    dt = 0.01 # the initial timestep should be very small to avoid discrepancy due to lagging at the beginning of each new Richardson iteration
  []
  dtmax = 5e3 # convergence of the global energy balance has less to do with the dt and end_time and more with the actual number of time steps!! (at least after the internal balance in both models is really close to 0)
  dtmin = 1e-3
[]

[UserObjects]
  [average_Twall_hp]
    type = NearestPointLayeredSideAverage
    boundary = 'hp'
    variable = Tsolid
    direction = z
    bounds = '0.2 2.0' # evaporator section - could add more layers but delta T less than 50K (do it after seeing if it messes with Sockeye)
    points_file = 'twelfth_core_hp_centers.txt' # assembly centers
    execute_on = 'initial timestep_end'
  []
[]

[MultiApps]
  [sockeye]
    type = TransientMultiApp
    positions_file = 'twelfth_core_hp_centers.txt' # 101 HPs
    app_type = SockeyeApp
    input_files = 'heatpipe_vapor_only.i'
    execute_on = 'timestep_end'
    max_procs_per_app = 1
    sub_cycling = true
    max_failures = 1000 # number of times the sub-app is allowed to cut its timesteps before the master-app cuts it
    output_in_position = true
    cli_args = 'weight=0.5 weight=1 weight=1 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=1 weight=1
                weight=0.5 weight=1 weight=1 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=1 weight=0.5
                weight=1 weight=1 weight=0.5 weight=1 weight=0.5 weight=0.5 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=1 weight=1 weight=0.5 weight=1
                weight=1 weight=1 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=1
                weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=1
                weight=1 weight=0.5 weight=1 weight=1 weight=1 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=0.5
                weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=1'
  []
[]

[Transfers]
  [to_sockeye_twall]
    type = MultiAppUserObjectTransfer
    to_multi_app = sockeye
    variable = T_ext_var
    user_object = average_Twall_hp
  []
  [from_sub_htc]
    type = MultiAppPostprocessorInterpolationTransfer
    from_multi_app = sockeye
    variable = htc_hp_var
    postprocessor = htc_equiv_cond
    num_points = 1 # interpolate with one point (~closest point)
    power = 0 # interpolate with constant function
  []
  [from_sub_tsink]
    type = MultiAppPostprocessorInterpolationTransfer
    from_multi_app = sockeye
    variable = Tsink_hp_var
    postprocessor = T_secondary #T_evap
    num_points = 1 # interpolate with one point (~closest point)
    power = 0 # interpolate with constant function
  []
  [total_power_from_condensers]
    type = MultiAppPostprocessorTransfer
    from_multi_app = sockeye
    from_postprocessor = 'total_condenser_power'
    to_postprocessor = 'total_condenser_power'
    reduction_type = sum
  []
[]

[Postprocessors]
  [fuel_temp_avg]
    type = ElementAverageValue
    variable = Tsolid
    block = ${fuel_blocks}
  []
  [fuel_temp_max]
    type = ElementExtremeValue
    variable = Tsolid
    block = ${fuel_blocks}
  []
  [fuel_temp_min]
    type = ElementExtremeValue
    variable = Tsolid
    block = ${fuel_blocks}
    value_type = min
  []
  [mod_temp_avg]
    type = ElementAverageValue
    variable = Tsolid
    block = ${moderator_blocks}
  []
  [mod_temp_max]
    type = ElementExtremeValue
    variable = Tsolid
    block = ${moderator_blocks}
  []
  [mod_temp_min]
    type = ElementExtremeValue
    variable = Tsolid
    block = ${moderator_blocks}
    value_type = min
  []
  [monolith_temp_avg]
    type = ElementAverageValue
    variable = Tsolid
    block = ${monolith_blocks}
  []
  [monolith_temp_max]
    type = ElementExtremeValue
    variable = Tsolid
    block = ${monolith_blocks}
  []
  [monolith_temp_min]
    type = ElementExtremeValue
    variable = Tsolid
    block = ${monolith_blocks}
    value_type = min
  []
  [refl_temp_avg]
    type = ElementAverageValue
    variable = Tsolid
    block = ${reflector_blocks}
  []
  [refl_temp_max]
    type = ElementExtremeValue
    variable = Tsolid
    block = ${reflector_blocks}
  []
  [refl_temp_min]
    type = ElementExtremeValue
    variable = Tsolid
    block = ${reflector_blocks}
    value_type = min
  []
  [heatpipe_surface_temp_avg]
    type = SideAverageValue
    variable = Tsolid
    boundary = 'hp'
  []
  [hp_bc_heat_flux]
    type = ConvectiveHeatTransferSideIntegral
    boundary = 'hp'
    T_solid = Tsolid
    T_fluid_var = Tsink_hp_var
    htc_var = htc_hp_var
    execute_on = 'initial timestep_end'
  []
  [rrefl_heat_flux]
    type = ConvectiveHeatTransferSideIntegral
    boundary = 'right'
    T_solid = Tsolid
    T_fluid = ${rrefl_Tsink}
    htc = ${rrefl_htc}
    execute_on = 'initial timestep_end'
  []
  [total_power]
    type = ElementIntegralVariablePostprocessor
    block = '${fuel_blocks}'
    variable = power_density
    execute_on = 'initial timestep_begin timestep_end'
  []
  [fuel_volume]
    type = VolumePostprocessor
    block = '${fuel_blocks}'
  []
  [hp_surface]
    type = AreaPostprocessor
    boundary = 'hp'
  []
  [Tsink_hp_var_avg]
    type = SideAverageValue
    variable = Tsink_hp_var
    # block = ${monolith_blocks}
    boundary = 'hp'
  []
  [htc_hp_var_avg]
    type = SideAverageValue
    variable = htc_hp_var
    boundary = 'hp'
  []
  [total_condenser_power]
    type = Receiver
  []
  [global_energy_balance] # balance performed globally (i.e. using the heat flux out of the condenser, much more difficult to bring to 0 than if using hp_bc_heat_flux)
    type = ParsedPostprocessor
    function = '(total_power - total_condenser_power - rrefl_heat_flux) / total_power'
    pp_names = 'total_power total_condenser_power rrefl_heat_flux'
  []
  [ch_avg]
    type = ElementAverageValue
    variable = ch
    block = ${moderator_blocks}
  []
  [ch_max]
    type = ElementExtremeValue
    variable = ch
    block = ${moderator_blocks}
  []
  [ch_min]
    type = ElementExtremeValue
    variable = ch
    block = ${moderator_blocks}
    value_type = min
  []
[]

[Outputs]
  perf_graph = true
  color = true
  csv = true
  [exodus]
    type = Exodus
    overwrite = true
    show = 'Tsolid power_density Tsink_hp_var htc_hp_var ch'
  []
[]

The Sockeye input contains the model for a single heat pipe. The latter uses the VOFM to model the 1D sodium flow within the core region of the HP and a 2D heat conduction model in the annulus and wick region. The full input is reported below for the sake of completeness.

# ##################################################################################################################
# AUTHORS:            Stefano Terlizzi and Vincent Labouré
# DATE:               June 2023
# ORGANIZATION:       Idaho National Laboratory (INL), C-110
# CITATION:           S. Terlizzi and V. Labouré. (2023).
#                     "Asymptotic hydrogen redistribution analysis in yttrium-hydride-moderated heat-pipe-cooled
#                     microreactors using DireWolf", Annals of Nuclear Energy, Vol. 186, 109735.
# DESCRIPTION:        Input for single heat pipe.
# FUNDS:              This work was supported by INL Laboratory Directed Research&
#                     Development (LDRD) Program under U.S. Department of Energy (DOE) Idaho Operations Office,
#                     United States of America contract DE-AC07-05ID14517.
# ##################################################################################################################
length_evap = 1.8
length_adia = 0.3
length_cond = 0.9

n_elems_evap = 60
n_elems_adia = 10
n_elems_cond = 30
n_elems_wick = 3
n_elems_ann = 3
n_elems_clad = 3

R_clad_o = 0.01
t_clad = 1.15e-03
t_ann = 8.35e-04
t_wick = 1.00e-03
R_pore = 1.5e-05
porosity = 7.06e-01
permeability = 1.80e-10

D_clad_o = '${fparse 2 * R_clad_o}'
D_clad_i = '${fparse D_clad_o - 2 * t_clad}'
D_wick_o = '${fparse D_clad_i - 2 * t_ann}'
D_wick_i = '${fparse D_wick_o - 2 * t_wick}'

magnitude_of_gravity = 9.805

S_evap = '${fparse pi * D_clad_o * length_evap}'
num_sides = 12 # number of sides of heat pipe as a result of mesh polygonization
alpha = '${fparse 2 * pi / num_sides}'
perimeter_correction = '${fparse 0.5 * alpha / sin(0.5 * alpha)}' # polygonization correction factor for perimeter
area_correction = '${fparse sqrt(alpha / sin(alpha))}' # polygonization correction factor for area
surface_correction = '${fparse area_correction / perimeter_correction}' # because mesh is volume preserving

T_sink = 900
htc_sink = 1e3

T_init = 950
T_ext = ${T_init}
htc_evap = 1e6 # chosen very large in lieu of a DirichletBC to obtain conservation

weight = 1 # set to 0.5 for half-HP (i.e. cut by the mesh)

fill_ratio = 1.01
T_ref_fill_ratio = '${fparse T_sink - 50}'

[GlobalParams]
  initial_T = ${T_init}
[]

[Closures]
  [sockeye]
    type = HeatPipeVOClosures
  []
[]

[FluidProperties]
  [fp_2phase]
    type = SodiumIdealGasTwoPhaseFluidProperties
  []
[]

[SolidProperties]
  [clad_mat]
    type = ThermalSS316Properties
  []
[]

[Components]
  [heat_pipe]
    type = HeatPipeVO

    position_evap_end = '0 0 0'
    direction_evap_to_cond = '0 0 1'

    n_elems_evap = ${n_elems_evap}
    n_elems_adia = ${n_elems_adia}
    n_elems_cond = ${n_elems_cond}

    n_elems_wick = ${n_elems_wick}
    n_elems_ann = ${n_elems_ann}
    n_elems_clad = ${n_elems_clad}

    fill_ratio = ${fill_ratio}
    T_ref_fill_ratio = ${T_ref_fill_ratio}

    T_ref_rho_wick = ${T_sink}
    T_ref_rho_clad = ${T_sink}

    gravity_vector = '0 -${magnitude_of_gravity} 0' # horizontal HP

    fp_2phase = fp_2phase
    sp_wick = clad_mat
    sp_clad = clad_mat
    closures = sockeye

    D_clad_o = ${D_clad_o}
    D_clad_i = ${D_clad_i}
    D_wick_i = ${D_wick_i}
    D_wick_o = ${D_wick_o}
    R_pore = ${R_pore}
    porosity = ${porosity}
    permeability = ${permeability}

    L_evap = ${length_evap}
    L_adia = ${length_adia}
    L_cond = ${length_cond}
    slope_reconstruction = NONE
    stop_vapor_at_condenser_pool = false
  []

  [evaporator_boundary]
    type = HSBoundaryExternalAppConvection
    boundary = 'heat_pipe:hs:evap:outer'
    hs = heat_pipe:hs
    T_ext = T_ext_var
    htc_ext = htc_evap_var
  []
  [adiabatic_boundary]
    type = HSBoundaryHeatFlux
    boundary = 'heat_pipe:hs:adia:outer'
    hs = heat_pipe:hs
    q = 0.0
  []
  [condenser_boundary]
    type = HSBoundaryExternalAppConvection
    boundary = 'heat_pipe:hs:cond:outer'
    hs = heat_pipe:hs
    T_ext = T_cond_var
    htc_ext = htc_cond_var
  []
[]

[AuxVariables]
  [T_ext_var]
    initial_condition = ${T_ext}
  []
  [htc_evap_var]
    initial_condition = ${htc_evap}
  []
  [htc_cond_var]
    initial_condition = ${htc_sink}
  []
  [T_cond_var]
    initial_condition = ${T_sink}
  []
[]

[Postprocessors]
  [Text_avg]
    type = SideAverageValue
    variable = T_ext_var
    boundary = 'heat_pipe:hs:evap:outer'
    execute_on = 'INITIAL TIMESTEP_END'
  []

  [htc_equiv_cond]
    type = ParsedPostprocessor
    expression = 'evaporator_boundary_integral / (Text_avg - T_secondary) / ${S_evap} / ${surface_correction}'
    pp_names = 'evaporator_boundary_integral Text_avg T_secondary' # evaporator_boundary_integral is the true numerical flux through the heat_pipe:hs:evap:outer
  []
  [T_secondary]
    type = Receiver
    default = ${T_sink}
  []
  [total_condenser_power]
    type = ParsedPostprocessor
    expression = '-${weight} * condenser_boundary_integral'
    pp_names = 'condenser_boundary_integral'
  []
[]

[Preconditioning]
  [pc]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  scheme = bdf2

  end_time = 3e4
  dtmin = 1e-8
  timestep_tolerance = 1e-6

  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.5
    dt = 1e-3
  []

  steady_state_detection = true
  steady_state_tolerance = 1e-8

  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-8
  nl_max_its = 15

  l_tol = 1e-3
  l_max_its = 10
[]

[Outputs]
  csv = true
  print_linear_converged_reason = false
  print_nonlinear_converged_reason = false
  print_linear_residuals = false

  [console]
    type = Console
    max_rows = 2
  []
[]

Multiphysics Coupling

Figure 3 shows the computational workflow for the coupled simulation. The left box includes the preliminary operations—namely, the mesh preparation and cross-section tabulations performed with Serpent (v. 2.1.32).

Figure 3: Integrated computational workflow (Terlizzi and Labouré, 2023).

Griffin is run to simulate the neutron transport in the core. The power density, $var element = document.getElementById("moose-equation-b0efa59c-29df-4898-abce-005799e9db07");katex.render("P_d", element, {displayMode:false,throwOnError:false});$ , computed from the neutron flux, is then transferred to Bison that is utilized to solve for the temperature in the reflector, $var element = document.getElementById("moose-equation-a3e65af2-b4e7-4f21-8ac5-05722d7665d5");katex.render("T_r", element, {displayMode:false,throwOnError:false});$ , moderator, $var element = document.getElementById("moose-equation-1d7f007c-696f-42e3-ba53-9154c0113fb8");katex.render("T_m", element, {displayMode:false,throwOnError:false});$ , and fuel, $var element = document.getElementById("moose-equation-41c4a1aa-4c87-49ea-8f50-cf452a59c0c1");katex.render("T_f", element, {displayMode:false,throwOnError:false});$ , and the hydrogen-yttrium stoichiometric ratio in the yttrium hydride, here denoted as $var element = document.getElementById("moose-equation-a0d8a398-68c3-4dfa-a2cf-db126463a4a2");katex.render("c_H", element, {displayMode:false,throwOnError:false});$ . Bison is coupled to 101 single heat pipes sub-applications that are used to compute the equivalent heat transfer coefficient, $var element = document.getElementById("moose-equation-ec65a172-23b2-4b1a-98c5-9d7130ac028b");katex.render("h_{eq}", element, {displayMode:false,throwOnError:false});$ , and the secondary-side temperature, $var element = document.getElementById("moose-equation-22e85533-0fbf-4610-a532-44e807f0909a");katex.render("T_{secondary}^i", element, {displayMode:false,throwOnError:false});$ , to obtain the thermal fluids' response. The updated temperature field and the hydrogen-yttrium stoichiometric ratio are then fed back into Griffin, where the macroscopic cross sections are updated based on the new operating conditions,

Running the Model

To run the input on the INL HPC, an interactive node can be requested, e.g., using:


qsub -I -l walltime=1:00:00  -l select=4:ncpus=48:mpiprocs=48 -P project

where 'project' should be replaced with the relevant project name (see here for a list of project names).

Then, load the following modules:


module load use.moose moose-apps direwolf

Finally, to run the input, make sure you are in the correct directory, then use the run command below. The input being ran will be the neutronics file.


mpirun -np 192 dire_wolf-opt -i neutronics_eigenvalue.i

With these commands, the run should take 20-25 minutes. The input can be also run with the BlueCRAB executable as long as it contains Sockeye using:


mpirun -np 192 blue_crab-opt -i neutronics_eigenvalue.i

References

Christopher Matthews, Vincent Laboure, Mark DeHart, Joshua Hansel, David Andrs, Yaqi Wang, Javier Ortensi, and Richard Martineau. Coupled multiphysics simulations of heat pipe microreactors using direwolf. Nuclear Technology, 207(7):1142–1162, 2021.

@article{Empire,
    author = "Matthews, Christopher and Laboure, Vincent and DeHart, Mark and Hansel, Joshua and Andrs, David and Wang, Yaqi and Ortensi, Javier and Martineau, Richard",
    title = "Coupled Multiphysics Simulations of Heat Pipe Microreactors Using DireWolf",
    journal = "Nuclear Technology",
    volume = "207",
    number = "7",
    pages = "1142-1162",
    year = "2021",
    publisher = "Taylor \& Francis"
}

Kan Ni, Yan Cao, Nicolas Stauff, and Jason Hou. Assessment of Griffin Cross-Section Interpolation Capability on TRISO-Fueled Heat-Pipe Micro-Reactor. In International Conference on Physics of Reactors 2022 (PHYSOR 2022). Pittsburgh PA, 2021. American Nuclear Society.

@inproceedings{PhysorANL,
    author = "Ni, Kan and Cao, Yan and Stauff, Nicolas and Hou, Jason",
    address = "Pittsburgh PA",
    booktitle = "International Conference on Physics of Reactors 2022 (PHYSOR 2022)",
    publisher = "American Nuclear Society",
    title = "{Assessment of Griffin Cross-Section Interpolation Capability on TRISO-Fueled Heat-Pipe Micro-Reactor}",
    year = "2021"
}

Emily Shemon, Yinbin Miao, Shikhar Kumar, Kun Mo, Yeon Sang Jung, Aaron Oaks, Guillaume Giudicelli, Logan Harbour, and Roy Stogner. MOOSE Reactor Module Meshing Enhancements to Support Reactor Analysis. Technical Report ANL/NSE-22/65, Argonne National Laboratory and Idaho National Laboratory, 2022.

@techreport{MOOSEReactorModule,
    author = "Shemon, Emily and Miao, Yinbin and Kumar, Shikhar and Mo, Kun and Jung, Yeon Sang and Oaks, Aaron and Giudicelli, Guillaume and Harbour, Logan and Stogner, Roy",
    number = "ANL/NSE-22/65",
    institution = "Argonne National Laboratory and Idaho National Laboratory",
    title = "{MOOSE Reactor Module Meshing Enhancements to Support Reactor Analysis}",
    year = "2022"
}

Stefano Terlizzi and Vincent Labouré. Asymptotic hydrogen redistribution analysis in yttrium-hydride-moderated heat-pipe-cooled microreactors using direwolf. Ann. Nucl. Energy, 186:109735, 2023. doi:https://doi.org/10.1016/j.anucene.2023.109735.

@article{Terlizzi2023,
    author = "Terlizzi, Stefano and Labouré, Vincent",
    title = "Asymptotic hydrogen redistribution analysis in yttrium-hydride-moderated heat-pipe-cooled microreactors using DireWolf",
    journal = "Ann. Nucl. Energy",
    volume = "186",
    pages = "109735",
    year = "2023",
    issn = "0306-4549",
    doi = "https://doi.org/10.1016/j.anucene.2023.109735"
}

Yaqi Wang, Changho Lee, Yeon Sang Jung, Zachary Prince, Joshua Hanophy, Logan Harbour, Hansol Park, and Javier Ortensi. Performance improvements to the griffin transport solvers. Technical Report INL/EXT-21-64272, Idaho National Laboratory, 2021.

@techreport{Wang2021perimp,
    author = "Wang, Yaqi and Lee, Changho and Jung, Yeon Sang and Prince, Zachary and Hanophy, Joshua and Harbour, Logan and Park, Hansol and Ortensi, Javier",
    title = "Performance Improvements to the Griffin Transport Solvers",
    Institution = "Idaho National Laboratory",
    Number = "INL/EXT-21-64272",
    place = "United States",
    Year = "2021"
}

(microreactors/hpmr_h2/steady/neutronics_eigenvalue.i)

# ##################################################################################################################
# AUTHORS:            Stefano Terlizzi and Vincent Labouré
# DATE:               June 2023
# ORGANIZATION:       Idaho National Laboratory (INL), C-110
# CITATION:           S. Terlizzi and V. Labouré. (2023).
#                     "Asymptotic hydrogen redistribution analysis in yttrium-hydride-moderated heat-pipe-cooled
#                     microreactors using DireWolf", Annals of Nuclear Energy, Vol. 186, 109735.
# DESCRIPTION:        Main input file containing neutronic model for the SiMBA reactor.
# FUNDS:              This work was supported by INL Laboratory Directed Research&
#                     Development (LDRD) Program under U.S. Department of Energy (DOE) Idaho Operations Office,
#                     United States of America contract DE-AC07-05ID14517.
# ##################################################################################################################
# Block assignment
bison_fuel_blocks = '1 2'
moderator_blocks = '3 4'
monolith_blocks = '8 22' # For now fill the air gap with monolith
reflector_blocks = '10 11 13 14 15 16 17 18 19'
absorber_blocks = '12'
bison_mod_blocks = '${moderator_blocks} ${monolith_blocks}'
bison_ref_blocks = '${reflector_blocks} ${absorber_blocks}'
[Mesh]
  [main]
    type = FileMeshGenerator
    file = 'one_twelfth_simba_core_in.e'
  []
  [coarse_mesh]
    type = FileMeshGenerator
    file = 'one_twelfth_simba_core_CM_in.e'
  []
  [mesh]
    type = CoarseMeshExtraElementIDGenerator
    input = main
    coarse_mesh = coarse_mesh
    extra_element_id_name = coarse_element_id
  []
[]

[PowerDensity]
  power = 1.66666e5 # power: 2e6 / 12 = 1.66666e5 W (per 1/12th core)
  power_density_variable = power_density
  integrated_power_postprocessor = integrated_power
[]

[TransportSystems]
  equation_type = eigenvalue
  particle = neutron
  G = 11

  ReflectingBoundary = 'bottom topleft'
  VacuumBoundary = 'right back front'

  [sn]
    scheme = DFEM-SN
    n_delay_groups = 6
    family = MONOMIAL
    order = FIRST
    AQtype = Gauss-Chebyshev
    NPolar = 1
    NAzmthl = 3
    NA = 1
    sweep_type = asynchronous_parallel_sweeper
    using_array_variable = true
    collapse_scattering = true
    hide_angular_flux = true
  []
[]

[GlobalParams]
  library_file = '../cross_sections_library/simba_fmh_b_fine.xml'
  library_name = simba_fmh_b_fine
  densities = 1.0
  isotopes = 'pseudo'
  dbgmat = false
  grid_names = 'Tfuel Tmod Trefl ch'
  grid_variables = 'Tfuel Tmod Trefl ch'
  is_meter = true
[]

[AuxVariables]
  [Tfuel] # fuel temperature
    initial_condition = 1000
  []
  [Tmod] # moderator + monolith + HP temperature
    initial_condition = 1000
  []
  [Trefl] # radial reflector temperature
    initial_condition = 900
  []
  [ch]
    initial_condition = 1.8
  []
[]

[Materials]
  [fuel]
    type = CoupledFeedbackNeutronicsMaterial
    block = '1 2' # fuel pin with 1 cm outer radius, no gap
    material_id = 1001
    plus = true
  []
  [moderator]
    type = CoupledFeedbackNeutronicsMaterial
    block = '3 4 5' # moderator pin with 0.975 cm outer radius
    material_id = 1002
  []
  [monolith]
    type = CoupledFeedbackNeutronicsMaterial
    block = '8'
    material_id = 1003
  []
  [hpipe]
    type = CoupledFeedbackNeutronicsMaterial
    block = '6 7' # gap homogenized with HP
    material_id = 1004
  []
  [be]
    type = CoupledFeedbackNeutronicsMaterial
    block = '10 11 13 14 15'
    material_id = 1005
  []
  [b4c]
    type = CoupledFeedbackNeutronicsMaterial
    block = '12'
    material_id = 1006
  []
  [air]
    type = CoupledFeedbackNeutronicsMaterial
    block = '20 21 22'
    material_id = 1007
  []
  [axial_refl]
    type = CoupledFeedbackNeutronicsMaterial
    block = '16 17 18 19'
    material_id = 1008
  []
[]

[Executioner]
  type = SweepUpdate
  # verbose = true
  # debug_richardson = true
  # Richardson
  richardson_rel_tol = 1e-3
  richardson_max_its = 50
  inner_solve_type = GMRes
  max_inner_its = 2
  # CMFD solver
  cmfd_acceleration = true
  coarse_element_id = coarse_element_id
  diffusion_eigen_solver_type = newton
  prolongation_type = multiplicative
  max_diffusion_coefficient = 1
  # Fixed point iteration with MultiApps
  fixed_point_solve_outer = true
  fixed_point_algorithm = picard
  custom_pp = eigenvalue
  custom_rel_tol = 1e-50
  custom_abs_tol = 1e-5
  force_fixed_point_solve = true
  fixed_point_max_its = 6
[]

[MultiApps]
  [bison]
    type = FullSolveMultiApp
    input_files = thermal_ss.i
    execute_on = 'timestep_end'
    keep_solution_during_restore = true # to restart from the latest solve of the multiapp (for pseudo-transient)
  []
[]

[Transfers]
  [pdens_to_modules]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    to_multi_app = bison
    source_variable = power_density
    variable = power_density
    from_blocks = '1 2'
    to_blocks = '${bison_fuel_blocks}'
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [tfuel_from_modules]
    type = MultiAppGeneralFieldNearestNodeTransfer
    from_multi_app = bison
    source_variable = Tsolid
    variable = Tfuel
    from_blocks = '${bison_fuel_blocks}'
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [tmod_from_modules]
    type = MultiAppGeneralFieldNearestNodeTransfer
    from_multi_app = bison
    source_variable = Tsolid
    variable = Tmod
    from_blocks = '${bison_mod_blocks}'
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [trefl_from_modules]
    type = MultiAppGeneralFieldNearestNodeTransfer
    from_multi_app = bison
    source_variable = Tsolid
    variable = Trefl
    from_blocks = '${bison_ref_blocks}'
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [ch_from_modules]
    type = MultiAppGeneralFieldNearestNodeTransfer
    from_multi_app = bison
    source_variable = ch
    variable = ch
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
[]

[Postprocessors]
  [fuel_temp_avg]
    type = ElementAverageValue
    variable = Tfuel
    block = '1 2'
  []
  [fuel_temp_max]
    type = ElementExtremeValue
    variable = Tfuel
    block = '1 2'
  []
  [fuel_temp_min]
    type = ElementExtremeValue
    variable = Tfuel
    value_type = min
    block = '1 2'
  []
  [mod_temp_avg]
    type = ElementAverageValue
    variable = Tmod
    block = '3 4'
  []
  [mod_temp_max]
    type = ElementExtremeValue
    variable = Tmod
    block = '3 4'
  []
  [mod_temp_min]
    type = ElementExtremeValue
    variable = Tmod
    value_type = min
    block = '3 4'
  []
  [fuel_volume]
    type = VolumePostprocessor
    block = '1 2'
  []
  [ch_avg]
    type = ElementAverageValue
    variable = ch
    block = '3 4'
  []
  [ch_max]
    type = ElementExtremeValue
    variable = ch
    block = '3 4'
  []
  [ch_min]
    type = ElementExtremeValue
    variable = ch
    value_type = min
    block = '3 4'
  []
[]

[Outputs]
  csv = true
  file_base = 'neutronics_eigenvalue'
  [console]
    type = Console
    outlier_variable_norms = false
  []
  [pgraph]
    type = PerfGraphOutput
    level = 2
  []
[]

(microreactors/hpmr_h2/steady/neutronics_eigenvalue.i)

# ##################################################################################################################
# AUTHORS:            Stefano Terlizzi and Vincent Labouré
# DATE:               June 2023
# ORGANIZATION:       Idaho National Laboratory (INL), C-110
# CITATION:           S. Terlizzi and V. Labouré. (2023).
#                     "Asymptotic hydrogen redistribution analysis in yttrium-hydride-moderated heat-pipe-cooled
#                     microreactors using DireWolf", Annals of Nuclear Energy, Vol. 186, 109735.
# DESCRIPTION:        Main input file containing neutronic model for the SiMBA reactor.
# FUNDS:              This work was supported by INL Laboratory Directed Research&
#                     Development (LDRD) Program under U.S. Department of Energy (DOE) Idaho Operations Office,
#                     United States of America contract DE-AC07-05ID14517.
# ##################################################################################################################
# Block assignment
bison_fuel_blocks = '1 2'
moderator_blocks = '3 4'
monolith_blocks = '8 22' # For now fill the air gap with monolith
reflector_blocks = '10 11 13 14 15 16 17 18 19'
absorber_blocks = '12'
bison_mod_blocks = '${moderator_blocks} ${monolith_blocks}'
bison_ref_blocks = '${reflector_blocks} ${absorber_blocks}'
[Mesh]
  [main]
    type = FileMeshGenerator
    file = 'one_twelfth_simba_core_in.e'
  []
  [coarse_mesh]
    type = FileMeshGenerator
    file = 'one_twelfth_simba_core_CM_in.e'
  []
  [mesh]
    type = CoarseMeshExtraElementIDGenerator
    input = main
    coarse_mesh = coarse_mesh
    extra_element_id_name = coarse_element_id
  []
[]

[PowerDensity]
  power = 1.66666e5 # power: 2e6 / 12 = 1.66666e5 W (per 1/12th core)
  power_density_variable = power_density
  integrated_power_postprocessor = integrated_power
[]

[TransportSystems]
  equation_type = eigenvalue
  particle = neutron
  G = 11

  ReflectingBoundary = 'bottom topleft'
  VacuumBoundary = 'right back front'

  [sn]
    scheme = DFEM-SN
    n_delay_groups = 6
    family = MONOMIAL
    order = FIRST
    AQtype = Gauss-Chebyshev
    NPolar = 1
    NAzmthl = 3
    NA = 1
    sweep_type = asynchronous_parallel_sweeper
    using_array_variable = true
    collapse_scattering = true
    hide_angular_flux = true
  []
[]

[GlobalParams]
  library_file = '../cross_sections_library/simba_fmh_b_fine.xml'
  library_name = simba_fmh_b_fine
  densities = 1.0
  isotopes = 'pseudo'
  dbgmat = false
  grid_names = 'Tfuel Tmod Trefl ch'
  grid_variables = 'Tfuel Tmod Trefl ch'
  is_meter = true
[]

[AuxVariables]
  [Tfuel] # fuel temperature
    initial_condition = 1000
  []
  [Tmod] # moderator + monolith + HP temperature
    initial_condition = 1000
  []
  [Trefl] # radial reflector temperature
    initial_condition = 900
  []
  [ch]
    initial_condition = 1.8
  []
[]

[Materials]
  [fuel]
    type = CoupledFeedbackNeutronicsMaterial
    block = '1 2' # fuel pin with 1 cm outer radius, no gap
    material_id = 1001
    plus = true
  []
  [moderator]
    type = CoupledFeedbackNeutronicsMaterial
    block = '3 4 5' # moderator pin with 0.975 cm outer radius
    material_id = 1002
  []
  [monolith]
    type = CoupledFeedbackNeutronicsMaterial
    block = '8'
    material_id = 1003
  []
  [hpipe]
    type = CoupledFeedbackNeutronicsMaterial
    block = '6 7' # gap homogenized with HP
    material_id = 1004
  []
  [be]
    type = CoupledFeedbackNeutronicsMaterial
    block = '10 11 13 14 15'
    material_id = 1005
  []
  [b4c]
    type = CoupledFeedbackNeutronicsMaterial
    block = '12'
    material_id = 1006
  []
  [air]
    type = CoupledFeedbackNeutronicsMaterial
    block = '20 21 22'
    material_id = 1007
  []
  [axial_refl]
    type = CoupledFeedbackNeutronicsMaterial
    block = '16 17 18 19'
    material_id = 1008
  []
[]

[Executioner]
  type = SweepUpdate
  # verbose = true
  # debug_richardson = true
  # Richardson
  richardson_rel_tol = 1e-3
  richardson_max_its = 50
  inner_solve_type = GMRes
  max_inner_its = 2
  # CMFD solver
  cmfd_acceleration = true
  coarse_element_id = coarse_element_id
  diffusion_eigen_solver_type = newton
  prolongation_type = multiplicative
  max_diffusion_coefficient = 1
  # Fixed point iteration with MultiApps
  fixed_point_solve_outer = true
  fixed_point_algorithm = picard
  custom_pp = eigenvalue
  custom_rel_tol = 1e-50
  custom_abs_tol = 1e-5
  force_fixed_point_solve = true
  fixed_point_max_its = 6
[]

[MultiApps]
  [bison]
    type = FullSolveMultiApp
    input_files = thermal_ss.i
    execute_on = 'timestep_end'
    keep_solution_during_restore = true # to restart from the latest solve of the multiapp (for pseudo-transient)
  []
[]

[Transfers]
  [pdens_to_modules]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    to_multi_app = bison
    source_variable = power_density
    variable = power_density
    from_blocks = '1 2'
    to_blocks = '${bison_fuel_blocks}'
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [tfuel_from_modules]
    type = MultiAppGeneralFieldNearestNodeTransfer
    from_multi_app = bison
    source_variable = Tsolid
    variable = Tfuel
    from_blocks = '${bison_fuel_blocks}'
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [tmod_from_modules]
    type = MultiAppGeneralFieldNearestNodeTransfer
    from_multi_app = bison
    source_variable = Tsolid
    variable = Tmod
    from_blocks = '${bison_mod_blocks}'
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [trefl_from_modules]
    type = MultiAppGeneralFieldNearestNodeTransfer
    from_multi_app = bison
    source_variable = Tsolid
    variable = Trefl
    from_blocks = '${bison_ref_blocks}'
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [ch_from_modules]
    type = MultiAppGeneralFieldNearestNodeTransfer
    from_multi_app = bison
    source_variable = ch
    variable = ch
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
[]

[Postprocessors]
  [fuel_temp_avg]
    type = ElementAverageValue
    variable = Tfuel
    block = '1 2'
  []
  [fuel_temp_max]
    type = ElementExtremeValue
    variable = Tfuel
    block = '1 2'
  []
  [fuel_temp_min]
    type = ElementExtremeValue
    variable = Tfuel
    value_type = min
    block = '1 2'
  []
  [mod_temp_avg]
    type = ElementAverageValue
    variable = Tmod
    block = '3 4'
  []
  [mod_temp_max]
    type = ElementExtremeValue
    variable = Tmod
    block = '3 4'
  []
  [mod_temp_min]
    type = ElementExtremeValue
    variable = Tmod
    value_type = min
    block = '3 4'
  []
  [fuel_volume]
    type = VolumePostprocessor
    block = '1 2'
  []
  [ch_avg]
    type = ElementAverageValue
    variable = ch
    block = '3 4'
  []
  [ch_max]
    type = ElementExtremeValue
    variable = ch
    block = '3 4'
  []
  [ch_min]
    type = ElementExtremeValue
    variable = ch
    value_type = min
    block = '3 4'
  []
[]

[Outputs]
  csv = true
  file_base = 'neutronics_eigenvalue'
  [console]
    type = Console
    outlier_variable_norms = false
  []
  [pgraph]
    type = PerfGraphOutput
    level = 2
  []
[]

(microreactors/hpmr_h2/steady/neutronics_eigenvalue.i)

# ##################################################################################################################
# AUTHORS:            Stefano Terlizzi and Vincent Labouré
# DATE:               June 2023
# ORGANIZATION:       Idaho National Laboratory (INL), C-110
# CITATION:           S. Terlizzi and V. Labouré. (2023).
#                     "Asymptotic hydrogen redistribution analysis in yttrium-hydride-moderated heat-pipe-cooled
#                     microreactors using DireWolf", Annals of Nuclear Energy, Vol. 186, 109735.
# DESCRIPTION:        Main input file containing neutronic model for the SiMBA reactor.
# FUNDS:              This work was supported by INL Laboratory Directed Research&
#                     Development (LDRD) Program under U.S. Department of Energy (DOE) Idaho Operations Office,
#                     United States of America contract DE-AC07-05ID14517.
# ##################################################################################################################
# Block assignment
bison_fuel_blocks = '1 2'
moderator_blocks = '3 4'
monolith_blocks = '8 22' # For now fill the air gap with monolith
reflector_blocks = '10 11 13 14 15 16 17 18 19'
absorber_blocks = '12'
bison_mod_blocks = '${moderator_blocks} ${monolith_blocks}'
bison_ref_blocks = '${reflector_blocks} ${absorber_blocks}'
[Mesh]
  [main]
    type = FileMeshGenerator
    file = 'one_twelfth_simba_core_in.e'
  []
  [coarse_mesh]
    type = FileMeshGenerator
    file = 'one_twelfth_simba_core_CM_in.e'
  []
  [mesh]
    type = CoarseMeshExtraElementIDGenerator
    input = main
    coarse_mesh = coarse_mesh
    extra_element_id_name = coarse_element_id
  []
[]

[PowerDensity]
  power = 1.66666e5 # power: 2e6 / 12 = 1.66666e5 W (per 1/12th core)
  power_density_variable = power_density
  integrated_power_postprocessor = integrated_power
[]

[TransportSystems]
  equation_type = eigenvalue
  particle = neutron
  G = 11

  ReflectingBoundary = 'bottom topleft'
  VacuumBoundary = 'right back front'

  [sn]
    scheme = DFEM-SN
    n_delay_groups = 6
    family = MONOMIAL
    order = FIRST
    AQtype = Gauss-Chebyshev
    NPolar = 1
    NAzmthl = 3
    NA = 1
    sweep_type = asynchronous_parallel_sweeper
    using_array_variable = true
    collapse_scattering = true
    hide_angular_flux = true
  []
[]

[GlobalParams]
  library_file = '../cross_sections_library/simba_fmh_b_fine.xml'
  library_name = simba_fmh_b_fine
  densities = 1.0
  isotopes = 'pseudo'
  dbgmat = false
  grid_names = 'Tfuel Tmod Trefl ch'
  grid_variables = 'Tfuel Tmod Trefl ch'
  is_meter = true
[]

[AuxVariables]
  [Tfuel] # fuel temperature
    initial_condition = 1000
  []
  [Tmod] # moderator + monolith + HP temperature
    initial_condition = 1000
  []
  [Trefl] # radial reflector temperature
    initial_condition = 900
  []
  [ch]
    initial_condition = 1.8
  []
[]

[Materials]
  [fuel]
    type = CoupledFeedbackNeutronicsMaterial
    block = '1 2' # fuel pin with 1 cm outer radius, no gap
    material_id = 1001
    plus = true
  []
  [moderator]
    type = CoupledFeedbackNeutronicsMaterial
    block = '3 4 5' # moderator pin with 0.975 cm outer radius
    material_id = 1002
  []
  [monolith]
    type = CoupledFeedbackNeutronicsMaterial
    block = '8'
    material_id = 1003
  []
  [hpipe]
    type = CoupledFeedbackNeutronicsMaterial
    block = '6 7' # gap homogenized with HP
    material_id = 1004
  []
  [be]
    type = CoupledFeedbackNeutronicsMaterial
    block = '10 11 13 14 15'
    material_id = 1005
  []
  [b4c]
    type = CoupledFeedbackNeutronicsMaterial
    block = '12'
    material_id = 1006
  []
  [air]
    type = CoupledFeedbackNeutronicsMaterial
    block = '20 21 22'
    material_id = 1007
  []
  [axial_refl]
    type = CoupledFeedbackNeutronicsMaterial
    block = '16 17 18 19'
    material_id = 1008
  []
[]

[Executioner]
  type = SweepUpdate
  # verbose = true
  # debug_richardson = true
  # Richardson
  richardson_rel_tol = 1e-3
  richardson_max_its = 50
  inner_solve_type = GMRes
  max_inner_its = 2
  # CMFD solver
  cmfd_acceleration = true
  coarse_element_id = coarse_element_id
  diffusion_eigen_solver_type = newton
  prolongation_type = multiplicative
  max_diffusion_coefficient = 1
  # Fixed point iteration with MultiApps
  fixed_point_solve_outer = true
  fixed_point_algorithm = picard
  custom_pp = eigenvalue
  custom_rel_tol = 1e-50
  custom_abs_tol = 1e-5
  force_fixed_point_solve = true
  fixed_point_max_its = 6
[]

[MultiApps]
  [bison]
    type = FullSolveMultiApp
    input_files = thermal_ss.i
    execute_on = 'timestep_end'
    keep_solution_during_restore = true # to restart from the latest solve of the multiapp (for pseudo-transient)
  []
[]

[Transfers]
  [pdens_to_modules]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    to_multi_app = bison
    source_variable = power_density
    variable = power_density
    from_blocks = '1 2'
    to_blocks = '${bison_fuel_blocks}'
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [tfuel_from_modules]
    type = MultiAppGeneralFieldNearestNodeTransfer
    from_multi_app = bison
    source_variable = Tsolid
    variable = Tfuel
    from_blocks = '${bison_fuel_blocks}'
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [tmod_from_modules]
    type = MultiAppGeneralFieldNearestNodeTransfer
    from_multi_app = bison
    source_variable = Tsolid
    variable = Tmod
    from_blocks = '${bison_mod_blocks}'
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [trefl_from_modules]
    type = MultiAppGeneralFieldNearestNodeTransfer
    from_multi_app = bison
    source_variable = Tsolid
    variable = Trefl
    from_blocks = '${bison_ref_blocks}'
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
  [ch_from_modules]
    type = MultiAppGeneralFieldNearestNodeTransfer
    from_multi_app = bison
    source_variable = ch
    variable = ch
    # Reduces transfers efficiency for now, can be removed once transferred fields are checked
    bbox_factor = 10
  []
[]

[Postprocessors]
  [fuel_temp_avg]
    type = ElementAverageValue
    variable = Tfuel
    block = '1 2'
  []
  [fuel_temp_max]
    type = ElementExtremeValue
    variable = Tfuel
    block = '1 2'
  []
  [fuel_temp_min]
    type = ElementExtremeValue
    variable = Tfuel
    value_type = min
    block = '1 2'
  []
  [mod_temp_avg]
    type = ElementAverageValue
    variable = Tmod
    block = '3 4'
  []
  [mod_temp_max]
    type = ElementExtremeValue
    variable = Tmod
    block = '3 4'
  []
  [mod_temp_min]
    type = ElementExtremeValue
    variable = Tmod
    value_type = min
    block = '3 4'
  []
  [fuel_volume]
    type = VolumePostprocessor
    block = '1 2'
  []
  [ch_avg]
    type = ElementAverageValue
    variable = ch
    block = '3 4'
  []
  [ch_max]
    type = ElementExtremeValue
    variable = ch
    block = '3 4'
  []
  [ch_min]
    type = ElementExtremeValue
    variable = ch
    value_type = min
    block = '3 4'
  []
[]

[Outputs]
  csv = true
  file_base = 'neutronics_eigenvalue'
  [console]
    type = Console
    outlier_variable_norms = false
  []
  [pgraph]
    type = PerfGraphOutput
    level = 2
  []
[]

(microreactors/hpmr_h2/steady/thermal_ss.i)

# ##################################################################################################################
# AUTHORS:            Stefano Terlizzi and Vincent Labouré
# DATE:               June 2023
# ORGANIZATION:       Idaho National Laboratory (INL), C-110
# CITATION:           S. Terlizzi and V. Labouré. (2023).
#                     "Asymptotic hydrogen redistribution analysis in yttrium-hydride-moderated heat-pipe-cooled
#                     microreactors using DireWolf", Annals of Nuclear Energy, Vol. 186, 109735.
# DESCRIPTION:        Model for heat conduction and hydrogen migration.
# FUNDS:              This work was supported by INL Laboratory Directed Research&
#                     Development (LDRD) Program under U.S. Department of Energy (DOE) Idaho Operations Office,
#                     United States of America contract DE-AC07-05ID14517.
# ##################################################################################################################
# Block assignment
fuel_blocks = '1 2'
moderator_blocks = '3 4'
monolith_blocks = '8 22' # For now fill the air gap with monolith
reflector_blocks = '10 11 13 14 15 16 17 18 19'
upper_reflector_blocks = '16 17'
absorber_blocks = '12'
# for radial reflector BC
rrefl_htc = 25 # W/m^2/K arbitrary (but small) HTC
rrefl_Tsink = 800 # K, arbitrary Tsink
T_init = 950
# Parameter to perform sensitivity analysis on Q, heat of transport
qvalue_multiplier = 1.0

[Mesh]
  [main]
    type = FileMeshGenerator
    file = 'one_twelfth_simba_core_in.e'
  []
  # remove blocks and add sidesets to apply boundary conditions
  [add_sideset_hp]
    type = SideSetsBetweenSubdomainsGenerator
    input = main
    primary_block = '8 17' # add 16 so the HP boundary extends into the upper axial reflector
    paired_block = '7'
    new_boundary = 'hp'
  []
  [add_sideset_inner_mod_gap]
    type = SideSetsBetweenSubdomainsGenerator
    input = add_sideset_hp
    primary_block = '4'
    paired_block = '5'
    new_boundary = 'gap_mod_inner'
  []
  [add_sideset_outer_mod_gap]
    type = SideSetsBetweenSubdomainsGenerator
    input = add_sideset_inner_mod_gap
    primary_block = '8'
    paired_block = '5'
    new_boundary = 'gap_mod_outer'
  []
  [add_sideset_central_hole]
    type = SideSetsBetweenSubdomainsGenerator
    input = add_sideset_outer_mod_gap
    primary_block = '22'
    paired_block = '21'
    new_boundary = 'central_hole'
  []
  [remove_blocks] # remove HPs, moderator gap and central shutdown cooling hole
    type = BlockDeletionGenerator
    input = add_sideset_central_hole
    block = '6 7 5 20 21'
  []
[]

[Variables]
  [Tsolid]
    initial_condition = ${T_init} # set all the temperatures to T_sink in Sockeye to allow a smooth ramp-up
  []
  [ch]
    initial_condition = 1.8
    block = ${moderator_blocks}
  []
[]

[Kernels]
  # Heat conduction
  [heat_time_derivative]
    type = HeatConductionTimeDerivative
    variable = Tsolid
  []
  [heat_conduction]
    type = HeatConduction
    variable = Tsolid
  []
  [heat_source_fuel]
    type = CoupledForce
    variable = Tsolid
    block = '${fuel_blocks}'
    v = power_density
  []
  # Hydrogen redistribution
  [ch_time_derivative]
    type = TimeDerivative
    variable = ch
    block = ${moderator_blocks}
  []
  [ch_dxdx]
    type = MatDiffusion
    variable = ch
    diffusivity = D
    block = ${moderator_blocks}
  []
  [soretDiff]
    type = ThermoDiffusion
    variable = ch
    temp = Tsolid
    heat_of_transport = Q
    mass_diffusivity = D
    block = ${moderator_blocks}
  []
[]

[AuxVariables]
  [power_density]
    block = '${fuel_blocks}'
    family = MONOMIAL
    order = FIRST
    initial_condition = 2302579.81614
  []
  # add upper_reflector_blocks blocks because HP cooling occurs there too
  [htc_hp_var]
    block = '${monolith_blocks} ${upper_reflector_blocks}'
    initial_condition = 372
  []
  [Tsink_hp_var]
    block = '${monolith_blocks} ${upper_reflector_blocks}'
    initial_condition = 900
  []
  [Tfuel_layered_average]
    order = CONSTANT
    family = MONOMIAL
    block = '${fuel_blocks}'
  []
[]

[BCs]
  [hp_bc]
    type = CoupledConvectiveHeatFluxBC
    variable = Tsolid
    boundary = 'hp' # inside of heat pipe
    htc = htc_hp_var # given by Sockeye
    T_infinity = Tsink_hp_var # given by Sockeye
  []
  [rrefl_bc]
    type = CoupledConvectiveHeatFluxBC
    variable = Tsolid
    boundary = 'right' # reflector cooling
    htc = ${rrefl_htc}
    T_infinity = ${rrefl_Tsink}
  []
[]

[ThermalContact]
  [moderator_gap]
    type = GapHeatTransfer
    emissivity_primary = 0.8
    emissivity_secondary = 0.8
    variable = Tsolid
    primary = 'gap_mod_inner'
    secondary = 'gap_mod_outer'
    gap_conductivity = 0.08 # W/m/K, typical thermal connectivity for air at 1000C
    quadrature = true
  []
[]

[Materials]
  #### DENSITY #####
  # units of kg/m^3
  [fuel_density]
    type = Density
    block = '${fuel_blocks}'
    density = 14.3e3 # same as in Serpent input
  []
  [moderator_density]
    type = Density
    block = '${moderator_blocks}'
    density = 4.3e3 # same as in Serpent input
  []
  [monolith_density]
    type = Density
    block = '${monolith_blocks}'
    density = 1.8e3 # same as in Serpent input
  []
  [reflector_density]
    type = Density
    block = '${reflector_blocks}'
    density = 1.85376e3 # same as in Serpent input
  []
  [absorber_density]
    type = Density
    block = '${absorber_blocks}'
    density = 2.52e3 # same as in Serpent input
  []
  ### THERMAL CONDUCTIVITY ###
  # units of W/m-K
  [fuel_kappa]
    type = HeatConductionMaterial
    block = '${fuel_blocks}'
    thermal_conductivity = 19 # W/m/K
    specific_heat = 300 # reasonable value
  []
  [moderator_kappa]
    type = HeatConductionMaterial
    block = '${moderator_blocks}'
    thermal_conductivity = 20 # W/m/K
    specific_heat = 500 # random value
  []
  [monolith_thermal_props]
    type = HeatConductionMaterial
    block = '${monolith_blocks}'
    thermal_conductivity = 70 # typical value for G348 graphite
    specific_heat = 1830 # typical value for G348 graphite
  []
  [reflector_kappa]
    type = HeatConductionMaterial
    block = '${reflector_blocks}'
    thermal_conductivity = 200 # W/m/K, typical value for Be
    specific_heat = 1825 # typical value for Be
  []
  [absorber_kappa]
    type = HeatConductionMaterial
    block = '${absorber_blocks}'
    thermal_conductivity = 20 # W/m/K, typical value for B4C
    specific_heat = 1000 # typical value for Be
  []
  # Material properties for migration. D can e adjusted for convergence in steady-state
  [moderator_diffusivity]
    type = GenericConstantMaterial
    block = ${moderator_blocks}
    prop_names = 'D Q'
    prop_values = '1 ${fparse 5.3e3 * qvalue_multiplier}'
  []
[]

[Functions]
[]

[Executioner]
  type = Transient
  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart'
  petsc_options_value = 'hypre boomeramg 300'
  line_search = 'none'
  automatic_scaling = true

  l_tol = 1e-02
  nl_abs_tol = 1e-7
  nl_rel_tol = 1e-8

  l_max_its = 50
  nl_max_its = 25

  start_time = 0
  end_time = 3e4

  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.5
    dt = 0.01 # the initial timestep should be very small to avoid discrepancy due to lagging at the beginning of each new Richardson iteration
  []
  dtmax = 5e3 # convergence of the global energy balance has less to do with the dt and end_time and more with the actual number of time steps!! (at least after the internal balance in both models is really close to 0)
  dtmin = 1e-3
[]

[UserObjects]
  [average_Twall_hp]
    type = NearestPointLayeredSideAverage
    boundary = 'hp'
    variable = Tsolid
    direction = z
    bounds = '0.2 2.0' # evaporator section - could add more layers but delta T less than 50K (do it after seeing if it messes with Sockeye)
    points_file = 'twelfth_core_hp_centers.txt' # assembly centers
    execute_on = 'initial timestep_end'
  []
[]

[MultiApps]
  [sockeye]
    type = TransientMultiApp
    positions_file = 'twelfth_core_hp_centers.txt' # 101 HPs
    app_type = SockeyeApp
    input_files = 'heatpipe_vapor_only.i'
    execute_on = 'timestep_end'
    max_procs_per_app = 1
    sub_cycling = true
    max_failures = 1000 # number of times the sub-app is allowed to cut its timesteps before the master-app cuts it
    output_in_position = true
    cli_args = 'weight=0.5 weight=1 weight=1 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=1 weight=1
                weight=0.5 weight=1 weight=1 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=1 weight=0.5
                weight=1 weight=1 weight=0.5 weight=1 weight=0.5 weight=0.5 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=1 weight=1 weight=0.5 weight=1
                weight=1 weight=1 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=1
                weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=1
                weight=1 weight=0.5 weight=1 weight=1 weight=1 weight=1 weight=1 weight=0.5 weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=0.5
                weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=1 weight=1'
  []
[]

[Transfers]
  [to_sockeye_twall]
    type = MultiAppUserObjectTransfer
    to_multi_app = sockeye
    variable = T_ext_var
    user_object = average_Twall_hp
  []
  [from_sub_htc]
    type = MultiAppPostprocessorInterpolationTransfer
    from_multi_app = sockeye
    variable = htc_hp_var
    postprocessor = htc_equiv_cond
    num_points = 1 # interpolate with one point (~closest point)
    power = 0 # interpolate with constant function
  []
  [from_sub_tsink]
    type = MultiAppPostprocessorInterpolationTransfer
    from_multi_app = sockeye
    variable = Tsink_hp_var
    postprocessor = T_secondary #T_evap
    num_points = 1 # interpolate with one point (~closest point)
    power = 0 # interpolate with constant function
  []
  [total_power_from_condensers]
    type = MultiAppPostprocessorTransfer
    from_multi_app = sockeye
    from_postprocessor = 'total_condenser_power'
    to_postprocessor = 'total_condenser_power'
    reduction_type = sum
  []
[]

[Postprocessors]
  [fuel_temp_avg]
    type = ElementAverageValue
    variable = Tsolid
    block = ${fuel_blocks}
  []
  [fuel_temp_max]
    type = ElementExtremeValue
    variable = Tsolid
    block = ${fuel_blocks}
  []
  [fuel_temp_min]
    type = ElementExtremeValue
    variable = Tsolid
    block = ${fuel_blocks}
    value_type = min
  []
  [mod_temp_avg]
    type = ElementAverageValue
    variable = Tsolid
    block = ${moderator_blocks}
  []
  [mod_temp_max]
    type = ElementExtremeValue
    variable = Tsolid
    block = ${moderator_blocks}
  []
  [mod_temp_min]
    type = ElementExtremeValue
    variable = Tsolid
    block = ${moderator_blocks}
    value_type = min
  []
  [monolith_temp_avg]
    type = ElementAverageValue
    variable = Tsolid
    block = ${monolith_blocks}
  []
  [monolith_temp_max]
    type = ElementExtremeValue
    variable = Tsolid
    block = ${monolith_blocks}
  []
  [monolith_temp_min]
    type = ElementExtremeValue
    variable = Tsolid
    block = ${monolith_blocks}
    value_type = min
  []
  [refl_temp_avg]
    type = ElementAverageValue
    variable = Tsolid
    block = ${reflector_blocks}
  []
  [refl_temp_max]
    type = ElementExtremeValue
    variable = Tsolid
    block = ${reflector_blocks}
  []
  [refl_temp_min]
    type = ElementExtremeValue
    variable = Tsolid
    block = ${reflector_blocks}
    value_type = min
  []
  [heatpipe_surface_temp_avg]
    type = SideAverageValue
    variable = Tsolid
    boundary = 'hp'
  []
  [hp_bc_heat_flux]
    type = ConvectiveHeatTransferSideIntegral
    boundary = 'hp'
    T_solid = Tsolid
    T_fluid_var = Tsink_hp_var
    htc_var = htc_hp_var
    execute_on = 'initial timestep_end'
  []
  [rrefl_heat_flux]
    type = ConvectiveHeatTransferSideIntegral
    boundary = 'right'
    T_solid = Tsolid
    T_fluid = ${rrefl_Tsink}
    htc = ${rrefl_htc}
    execute_on = 'initial timestep_end'
  []
  [total_power]
    type = ElementIntegralVariablePostprocessor
    block = '${fuel_blocks}'
    variable = power_density
    execute_on = 'initial timestep_begin timestep_end'
  []
  [fuel_volume]
    type = VolumePostprocessor
    block = '${fuel_blocks}'
  []
  [hp_surface]
    type = AreaPostprocessor
    boundary = 'hp'
  []
  [Tsink_hp_var_avg]
    type = SideAverageValue
    variable = Tsink_hp_var
    # block = ${monolith_blocks}
    boundary = 'hp'
  []
  [htc_hp_var_avg]
    type = SideAverageValue
    variable = htc_hp_var
    boundary = 'hp'
  []
  [total_condenser_power]
    type = Receiver
  []
  [global_energy_balance] # balance performed globally (i.e. using the heat flux out of the condenser, much more difficult to bring to 0 than if using hp_bc_heat_flux)
    type = ParsedPostprocessor
    function = '(total_power - total_condenser_power - rrefl_heat_flux) / total_power'
    pp_names = 'total_power total_condenser_power rrefl_heat_flux'
  []
  [ch_avg]
    type = ElementAverageValue
    variable = ch
    block = ${moderator_blocks}
  []
  [ch_max]
    type = ElementExtremeValue
    variable = ch
    block = ${moderator_blocks}
  []
  [ch_min]
    type = ElementExtremeValue
    variable = ch
    block = ${moderator_blocks}
    value_type = min
  []
[]

[Outputs]
  perf_graph = true
  color = true
  csv = true
  [exodus]
    type = Exodus
    overwrite = true
    show = 'Tsolid power_density Tsink_hp_var htc_hp_var ch'
  []
[]

(microreactors/hpmr_h2/steady/heatpipe_vapor_only.i)

# ##################################################################################################################
# AUTHORS:            Stefano Terlizzi and Vincent Labouré
# DATE:               June 2023
# ORGANIZATION:       Idaho National Laboratory (INL), C-110
# CITATION:           S. Terlizzi and V. Labouré. (2023).
#                     "Asymptotic hydrogen redistribution analysis in yttrium-hydride-moderated heat-pipe-cooled
#                     microreactors using DireWolf", Annals of Nuclear Energy, Vol. 186, 109735.
# DESCRIPTION:        Input for single heat pipe.
# FUNDS:              This work was supported by INL Laboratory Directed Research&
#                     Development (LDRD) Program under U.S. Department of Energy (DOE) Idaho Operations Office,
#                     United States of America contract DE-AC07-05ID14517.
# ##################################################################################################################
length_evap = 1.8
length_adia = 0.3
length_cond = 0.9

n_elems_evap = 60
n_elems_adia = 10
n_elems_cond = 30
n_elems_wick = 3
n_elems_ann = 3
n_elems_clad = 3

R_clad_o = 0.01
t_clad = 1.15e-03
t_ann = 8.35e-04
t_wick = 1.00e-03
R_pore = 1.5e-05
porosity = 7.06e-01
permeability = 1.80e-10

D_clad_o = '${fparse 2 * R_clad_o}'
D_clad_i = '${fparse D_clad_o - 2 * t_clad}'
D_wick_o = '${fparse D_clad_i - 2 * t_ann}'
D_wick_i = '${fparse D_wick_o - 2 * t_wick}'

magnitude_of_gravity = 9.805

S_evap = '${fparse pi * D_clad_o * length_evap}'
num_sides = 12 # number of sides of heat pipe as a result of mesh polygonization
alpha = '${fparse 2 * pi / num_sides}'
perimeter_correction = '${fparse 0.5 * alpha / sin(0.5 * alpha)}' # polygonization correction factor for perimeter
area_correction = '${fparse sqrt(alpha / sin(alpha))}' # polygonization correction factor for area
surface_correction = '${fparse area_correction / perimeter_correction}' # because mesh is volume preserving

T_sink = 900
htc_sink = 1e3

T_init = 950
T_ext = ${T_init}
htc_evap = 1e6 # chosen very large in lieu of a DirichletBC to obtain conservation

weight = 1 # set to 0.5 for half-HP (i.e. cut by the mesh)

fill_ratio = 1.01
T_ref_fill_ratio = '${fparse T_sink - 50}'

[GlobalParams]
  initial_T = ${T_init}
[]

[Closures]
  [sockeye]
    type = HeatPipeVOClosures
  []
[]

[FluidProperties]
  [fp_2phase]
    type = SodiumIdealGasTwoPhaseFluidProperties
  []
[]

[SolidProperties]
  [clad_mat]
    type = ThermalSS316Properties
  []
[]

[Components]
  [heat_pipe]
    type = HeatPipeVO

    position_evap_end = '0 0 0'
    direction_evap_to_cond = '0 0 1'

    n_elems_evap = ${n_elems_evap}
    n_elems_adia = ${n_elems_adia}
    n_elems_cond = ${n_elems_cond}

    n_elems_wick = ${n_elems_wick}
    n_elems_ann = ${n_elems_ann}
    n_elems_clad = ${n_elems_clad}

    fill_ratio = ${fill_ratio}
    T_ref_fill_ratio = ${T_ref_fill_ratio}

    T_ref_rho_wick = ${T_sink}
    T_ref_rho_clad = ${T_sink}

    gravity_vector = '0 -${magnitude_of_gravity} 0' # horizontal HP

    fp_2phase = fp_2phase
    sp_wick = clad_mat
    sp_clad = clad_mat
    closures = sockeye

    D_clad_o = ${D_clad_o}
    D_clad_i = ${D_clad_i}
    D_wick_i = ${D_wick_i}
    D_wick_o = ${D_wick_o}
    R_pore = ${R_pore}
    porosity = ${porosity}
    permeability = ${permeability}

    L_evap = ${length_evap}
    L_adia = ${length_adia}
    L_cond = ${length_cond}
    slope_reconstruction = NONE
    stop_vapor_at_condenser_pool = false
  []

  [evaporator_boundary]
    type = HSBoundaryExternalAppConvection
    boundary = 'heat_pipe:hs:evap:outer'
    hs = heat_pipe:hs
    T_ext = T_ext_var
    htc_ext = htc_evap_var
  []
  [adiabatic_boundary]
    type = HSBoundaryHeatFlux
    boundary = 'heat_pipe:hs:adia:outer'
    hs = heat_pipe:hs
    q = 0.0
  []
  [condenser_boundary]
    type = HSBoundaryExternalAppConvection
    boundary = 'heat_pipe:hs:cond:outer'
    hs = heat_pipe:hs
    T_ext = T_cond_var
    htc_ext = htc_cond_var
  []
[]

[AuxVariables]
  [T_ext_var]
    initial_condition = ${T_ext}
  []
  [htc_evap_var]
    initial_condition = ${htc_evap}
  []
  [htc_cond_var]
    initial_condition = ${htc_sink}
  []
  [T_cond_var]
    initial_condition = ${T_sink}
  []
[]

[Postprocessors]
  [Text_avg]
    type = SideAverageValue
    variable = T_ext_var
    boundary = 'heat_pipe:hs:evap:outer'
    execute_on = 'INITIAL TIMESTEP_END'
  []

  [htc_equiv_cond]
    type = ParsedPostprocessor
    expression = 'evaporator_boundary_integral / (Text_avg - T_secondary) / ${S_evap} / ${surface_correction}'
    pp_names = 'evaporator_boundary_integral Text_avg T_secondary' # evaporator_boundary_integral is the true numerical flux through the heat_pipe:hs:evap:outer
  []
  [T_secondary]
    type = Receiver
    default = ${T_sink}
  []
  [total_condenser_power]
    type = ParsedPostprocessor
    expression = '-${weight} * condenser_boundary_integral'
    pp_names = 'condenser_boundary_integral'
  []
[]

[Preconditioning]
  [pc]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  scheme = bdf2

  end_time = 3e4
  dtmin = 1e-8
  timestep_tolerance = 1e-6

  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.5
    dt = 1e-3
  []

  steady_state_detection = true
  steady_state_tolerance = 1e-8

  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-8
  nl_max_its = 15

  l_tol = 1e-3
  l_max_its = 10
[]

[Outputs]
  csv = true
  print_linear_converged_reason = false
  print_nonlinear_converged_reason = false
  print_linear_residuals = false

  [console]
    type = Console
    max_rows = 2
  []
[]

Mesh
Cross Sections
Griffin Input for Neutronics
Bison Input for Heat Conduction and Hydrogen Redistribution
Multiphysics Coupling
Running the Model
References

HPMR-Hvar element = document.getElementById("moose-equation-cd83828e-bef4-4fe7-9cf1-9c7e15d27cab");katex.render("_2", element, {displayMode:false,throwOnError:false}); Reactor Model

Mesh

Cross Sections

Griffin Input for Neutronics

Bison Input for Heat Conduction and Hydrogen Redistribution

Multiphysics Coupling

Running the Model

References

(microreactors/hpmr_h2/steady/neutronics_eigenvalue.i)

(microreactors/hpmr_h2/steady/neutronics_eigenvalue.i)

(microreactors/hpmr_h2/steady/neutronics_eigenvalue.i)

(microreactors/hpmr_h2/steady/thermal_ss.i)

(microreactors/hpmr_h2/steady/heatpipe_vapor_only.i)

HPMR-H $var element = document.getElementById("moose-equation-cd83828e-bef4-4fe7-9cf1-9c7e15d27cab");katex.render("_2", element, {displayMode:false,throwOnError:false});$ Reactor Model