High Temperature Engineering Test Reactor (HTTR) Null Transient Model Description

Contact: Javier Ortensi javier.ortensi.at.inl.gov

Model link: HTTR Null Transient

A detailed description of the steady-state model is available here. Only modifications necessary to set up a null-transient state are explained below.

For all null-transient input files, initial conditions ([ICs] block and initial_condition parameters) have been removed. All child applications (thermal hydraulics and fuel element heterogeneous heat conduction) are restarted using the parent full core heat conduction checkpoint. The full-core neutronics model uses an Exodus file for restart for the mesh as well as Griffin's binary restart.

Neutronics

Mesh

Parameter use_for_exodus_restart has been set to true to enable the [Mesh] block to restart the fluxes from the exodus file populated from the steady state neutronics simulation.

[Mesh]
  [fmg]
    type = FileMeshGenerator
    file = 'neutronics_eigenvalue${extension}.e'
    exodus_extra_element_integers = 'material_id equivalence_id'
    use_for_exodus_restart = true
  []
  # Restart relies on the ExodusII_IO functionality, which only works with ReplicatedMesh.
  parallel_type = replicated
[]
(htgr/httr/steady_state_and_null_transient/neutronics_null.i)

Transport System

The [TransportSystems] block restarts with the eigenvalue solved from the steady state neutronics calculation, to scale the fission cross sections so the core is critical. The parameter equation_type has been changed to transient.

[TransportSystems]
  particle = neutron
  equation_type = transient
  restart_transport_system = true
  G = 10
  VacuumBoundary = '1 2 3'

  [diff]
    scheme = CFEM-Diffusion
    n_delay_groups = 6
    assemble_scattering_jacobian = true
    assemble_fission_jacobian = true
    assemble_delay_jacobian = true
  []
[]
(htgr/httr/steady_state_and_null_transient/neutronics_null.i)

User Objects

Values not included in the exodus file are restarted from binary in the [UserObjects] block using a XXXXXXXXXX

[UserObjects]
  [ss_temperatures]
    type = SolutionVectorFile
    var = 'Tfuel Tmod'
    writing = false
    execute_on = 'initial'
  []
  [restart_poison_densities]
    type = SolutionVectorFile
    var = 'NI NXe'
    writing = false
    execute_on = 'initial'
  []
[]
(htgr/httr/steady_state_and_null_transient/neutronics_null.i)

MultiApps

The [MultiApps] block has been modified to utilize a TransientMultiApp type. This is because all the physics simulations will progress synchronized in time, there won't be an eigenvalue neutronics calculation as there was previously for the steady state.

[MultiApps]
  [moose_modules]
    type = TransientMultiApp
    positions = '0 0 0'
    input_files = full_core_ht_null.i
    execute_on = 'timestep_end'
  []
[]
(htgr/httr/steady_state_and_null_transient/neutronics_null.i)

Heat Transfer

Heat Transfer - Homogenized Full Core

Mesh

The [Mesh] block now loads the mesh the checkpoint created from the steady simulation homogenized simulation.

[Mesh]
  file = neutronics_eigenvalue${extension}_moose_modules0_checkpoint_cp/LATEST
[]
(htgr/httr/steady_state_and_null_transient/full_core_ht_null.i)

Executioner

The Executioner block has been updated to use a Transient type and include time steps for a 50 second transient run.

[Executioner]
  type = Transient
  automatic_scaling = true
  compute_scaling_once = false
  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart '
  petsc_options_value = 'hypre boomeramg 100'

  start_time = 0
  end_time = 50
  dt = 10

  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-9
  fixed_point_rel_tol = 1e-3
  fixed_point_abs_tol = 1e-7

  line_search = none # default seems bad for convergence
[]
(htgr/httr/steady_state_and_null_transient/full_core_ht_null.i)

Heat Transfer - Single Pin Heterogeneous

Mesh

The Mesh block has been restarted with the checkpoint created from the steady state single pin heterogeneous calculation.

[Mesh]
  [fmg]
    type = FileMeshGenerator
    file = '../mesh/fuel_element/HTTR_fuel_pin_2D_refined_m_5pins_axialref.e'
  []
  uniform_refine = 0 # do not modify if using ThermalContact
  coord_type = RZ
  rz_coord_axis = X
[]
(htgr/httr/steady_state_and_null_transient/fuel_elem_null.i)

Kernels

The Kernels block has time derivatives added back into the nested heat_ie block.

[Kernels]
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temp
  []
[]
(htgr/httr/steady_state_and_null_transient/fuel_elem_null.i)

Executioner

The Executioner block has been updated to use a Transient type and include time steps for a 50 second transient run.

[Executioner]
  type = Transient

  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-6

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu       superlu_dist'

  start_time = 0
  end_time = 50
  dt = 10
[]
(htgr/httr/steady_state_and_null_transient/fuel_elem_null.i)

Thermal-Hydraulics

Thermal-Hydraulics - Control Rods

Executioner

The Executioner block has been updated to use a Transient type and include time steps for a 50 second transient run.

[Executioner]
  type = Transient
  scheme = 'bdf2'
  solve_type = 'NEWTON'
  line_search = 'basic'

  dtmin = 1.e-8

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e-2
  []

  nl_rel_tol = 1e-6
  nl_abs_tol = 5e-8
  nl_max_its = 20

  l_tol = 1e-3
  l_max_its = 100

  start_time = 0
  end_time = 50

  [Quadrature]
    type = GAUSS
    order = SECOND
  []
[]
(htgr/httr/steady_state_and_null_transient/thermal_hydraulics_CR_null.i)

Thermal-Hydraulics - Fuel Pins

Executioner

The Executioner block has been updated to use a Transient type and include time steps for a 50 second transient run.

[Executioner]
  type = Transient
  scheme = 'bdf2'
  solve_type = 'NEWTON'
  line_search = 'basic'

  dtmin = 1.e-8

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e-2
  []

  nl_rel_tol = 1e-6
  nl_abs_tol = 5e-8 # original: 1e-8 (but can result in solve failure in rare cases)
  nl_max_its = 20 # original: 5 (but can result in solve failure in rare cases)

  l_tol = 1e-3
  l_max_its = 100

  start_time = 0
  end_time = 50

  [Quadrature]
    type = GAUSS
    order = SECOND
  []
[]
(htgr/httr/steady_state_and_null_transient/thermal_hydraulics_fuel_pins_null.i)

Null Transient Model input summary:

The inputs can be found below or on the relevant folder on the Github repository.

Griffin model

# ==============================================================================
# Full-core homogenized neutronics model (null)
# Application: Griffin (Sabertooth with child applications)
# ------------------------------------------------------------------------------
# Idaho Falls, INL, April 2023
# Author(s): Vincent Laboure
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
# ==============================================================================

# Defining i as the axial layer index (1 being the top and 18 the bottom), we have the following block IDs:
# 1050 + i: fuel columns of type 1
# 1150 + i: fuel columns of type 2
# 1250 + i: fuel columns of type 3
# 1350 + i: fuel columns of type 4
# 1450 + i: RR columns
# 1550 + i: PR columns
# 2050 + i: CR column, type C
# 2150 + i: CR columns, type R1 without neutron sources for i = 5-6
# 2170 + i: CR columns, type R1 with neutron sources for i = 5-6 (only modeled during transient)
# 2250 + i: CR columns, type R2
# 2350 + i: I columns,
# 2450 + i: CR columns, type R3

# Block IDs of fuel regions (to be SPH corrected)
fuel_blocks = '
1055 1056 1057 1058 1059 1060 1061 1062 1063 1064
1155 1156 1157 1158 1159 1160 1161 1162 1163 1164
1255 1256 1257 1258 1259 1260 1261 1262 1263 1264
1355 1356 1357 1358 1359 1360 1361 1362 1363 1364
'

# Block IDs of non-fuel (CR, RR, I) stack regions to be SPH corrected (active core + one layer of axial reflector on top and bottom)
sph_nonfuel_blocks = '
1054                                                   1065
1154                                                   1165
1254                                                   1265
1354                                                   1365
1454 1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 1465
2054 2055 2056 2057 2058 2059 2060 2061 2062 2063 2064 2065
2154 2155 2156 2157 2158 2159 2160 2161 2162 2163 2164 2165
2174 2175 2176 2177 2178 2179 2180 2181 2182 2183 2184 2185
2254 2255 2256 2257 2258 2259 2260 2261 2262 2263 2264 2265
2354 2355 2356 2357 2358 2359 2360 2361 2362 2363 2364 2365
2454 2455 2456 2457 2458 2459 2460 2461 2462 2463 2464 2465
'

# Block IDs of the part of the axial reflector to not be SPH corrected (3 layers on top and bottom)
axial_reflector_blocks = '
1051 1052 1053 1066 1067 1068 2051 2052 2053 2066 2067 2068
1151 1152 1153 1166 1167 1168 2151 2152 2153 2166 2167 2168 2171 2172 2173 2186 2187 2188
1251 1252 1253 1266 1267 1268 2251 2252 2253 2266 2267 2268
1351 1352 1353 1366 1367 1368 2351 2352 2353 2366 2367 2368
1451 1452 1453 1466 1467 1468 2451 2452 2453 2466 2467 2468
'

# Block IDs of the permanent reflector (PR)
radial_reflector_blocks = '
1551 1552 1553 1554 1555 1556 1557 1558 1559 1560 1561 1562 1563 1564 1565 1566 1567 1568
'

initial_decay_heat = 5.7472e5 # decay heat for 9MW steady-state, in W
total_power = 9e6 # total power (including decay heat), in W
fission_power = '${fparse total_power - initial_decay_heat}'

extension = '_9MW'

[Mesh]
  [fmg]
    type = FileMeshGenerator
    file = 'neutronics_eigenvalue${extension}.e'
    exodus_extra_element_integers = 'material_id equivalence_id'
    use_for_exodus_restart = true
  []
  # Restart relies on the ExodusII_IO functionality, which only works with ReplicatedMesh.
  parallel_type = replicated
[]

[TransportSystems]
  particle = neutron
  equation_type = transient
  restart_transport_system = true
  G = 10
  VacuumBoundary = '1 2 3'

  [diff]
    scheme = CFEM-Diffusion
    n_delay_groups = 6
    assemble_scattering_jacobian = true
    assemble_fission_jacobian = true
    assemble_delay_jacobian = true
  []
[]

[PoisonTracking]
  decay_chains = 'XE135'
  micro_library_file = '../cross_sections/HTTR_5x5${extension}_profiled_VR_kappa_adjusted.xml'
  micro_library_name = 'HTTR_5x5${extension}'
  grid_names = 'Tfuel Tmod'
  grid_variables = 'Tfuel Tmod'
[]

[Equivalence]
  type = SPH
  transport_system = diff
  compute_factors = false
  equivalence_library = '../cross_sections/HTTR_5x5${extension}_equiv_corrected.xml'
  library_name = 'HTTR_5x5_equiv${extension}'

  equivalence_grid_names = 'Tfuel Tmod'
  equivalence_grid_variables = 'Tfuel Tmod'
  sph_block = '${fuel_blocks} ${sph_nonfuel_blocks}'
[]

[PowerDensity]
  power = ${fission_power} # initial power
  power_density_variable = power_density
  integrated_power_postprocessor = integrated_power
[]

[GlobalParams]
  # for materials
  library_file = '../cross_sections/HTTR_5x5${extension}_profiled_VR_kappa_adjusted.xml'
  library_name = 'HTTR_5x5${extension}'
  isotopes = 'pseudo'
  densities = '1.0'
  is_meter = true
  grid_names = 'Tfuel Tmod'
  grid_variables = 'Tfuel Tmod'

  constant_on = ELEMENT
[]

[Materials]
  [fuel_sph]
    type = CoupledFeedbackMatIDNeutronicsMaterial
    block = '${fuel_blocks}'
    plus = true
  []
  [nonfuel_sph]
    type = CoupledFeedbackMatIDNeutronicsMaterial
    block = '${sph_nonfuel_blocks}'
  []
  [nonsph]
    type = CoupledFeedbackMatIDNeutronicsMaterial
    block = '${axial_reflector_blocks} ${radial_reflector_blocks}'
  []
[]

[AuxVariables]
  [Tfuel]
    order = CONSTANT
    family = MONOMIAL
  []
  [Tmod]
    order = FIRST
    family = L2_LAGRANGE
  []
[]

[Postprocessors]
  [Tm_avg]
    type = ElementAverageValue
    variable = Tmod
    execute_on = 'initial timestep_end'
  []
  [Tf_avg]
    type = ElementAverageValue
    variable = Tfuel
    execute_on = 'initial timestep_end'
  []
  [Tm_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tmod
    execute_on = 'initial timestep_end'
  []
  [Tf_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tfuel
    execute_on = 'initial timestep_end'
  []
  [NI_avg]
    type = ElementAverageValue
    variable = NI
    block = '${fuel_blocks}'
    execute_on = 'initial timestep_end'
  []
  [NXe_avg]
    type = ElementAverageValue
    variable = NXe
    block = '${fuel_blocks}'
    execute_on = 'initial timestep_end'
  []
[]

[UserObjects]
  [ss_temperatures]
    type = SolutionVectorFile
    var = 'Tfuel Tmod'
    writing = false
    execute_on = 'initial'
  []
  [restart_poison_densities]
    type = SolutionVectorFile
    var = 'NI NXe'
    writing = false
    execute_on = 'initial'
  []
[]

[Preconditioning]
  [SMP]
    type = SMP
    off_diag_row = '
                       sflux_g1
                       sflux_g2
                       sflux_g3
                       sflux_g4
                       sflux_g5 sflux_g5
                       sflux_g6 sflux_g6 sflux_g6 sflux_g6
                       sflux_g7 sflux_g7 sflux_g7
                       sflux_g8 sflux_g8 sflux_g8
                       sflux_g9
                      '
    off_diag_column = '
                       sflux_g0
                       sflux_g1
                       sflux_g2
                       sflux_g3
                       sflux_g4 sflux_g6
                       sflux_g5 sflux_g7 sflux_g8 sflux_g9
                       sflux_g6 sflux_g8 sflux_g9
                       sflux_g6 sflux_g7 sflux_g9
                       sflux_g8
                      '
  []
[]

[Executioner]
  type = Transient

  start_time = 0
  end_time = 50
  dt = 10

  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart '
  petsc_options_value = 'hypre boomeramg 100'
  line_search = none # default seems bad for convergence

  automatic_scaling = true
  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-8
  l_tol = 1e-2

  fixed_point_rel_tol = 1e-7
  fixed_point_abs_tol = 1e-4
  fixed_point_max_its = 10
  accept_on_max_fixed_point_iteration = true
  disable_fixed_point_residual_norm_check = false
[]

[MultiApps]
  [moose_modules]
    type = TransientMultiApp
    positions = '0 0 0'
    input_files = full_core_ht_null.i
    execute_on = 'timestep_end'
  []
[]

[Transfers]
  [pdens_to_modules]
    type = MultiAppProjectionTransfer
    to_multi_app = moose_modules
    source_variable = power_density
    variable = power_density
    execute_on = 'timestep_end'
  []
  [tmod_from_modules]
    type = MultiAppProjectionTransfer
    from_multi_app = moose_modules
    source_variable = Tmod
    variable = Tmod
    execute_on = 'timestep_end'
  []
  [tfuel_from_modules]
    type = MultiAppProjectionTransfer
    from_multi_app = moose_modules
    source_variable = Tfuel
    variable = Tfuel
    execute_on = 'timestep_end'
  []
[]

[Outputs]
  [exodus]
    type = Exodus
    overwrite = true
  []

  csv = true
  perf_graph = true
[]
(htgr/httr/steady_state_and_null_transient/neutronics_null.i)

Full-core heat transfer model

# ==============================================================================
# Full-core homogenized heat transfer model (null)
# Application: BISON (Sabertooth with child applications)
# ------------------------------------------------------------------------------
# Idaho Falls, INL, April 2023
# Author(s): Vincent Laboure
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
# ==============================================================================

# Defining i as the axial layer index (1 being the top and 18 the bottom), we have the following block IDs:
# 1050 + i: fuel columns of type 1
# 1150 + i: fuel columns of type 2
# 1250 + i: fuel columns of type 3
# 1350 + i: fuel columns of type 4
# 1450 + i: RR columns
# 1550 + i: PR columns
# 2050 + i: CR column, type C
# 2150 + i: CR columns, type R1 without neutron sources for i = 5-6
# 2170 + i: CR columns, type R1 with neutron sources for i = 5-6 (only modeled during transient)
# 2250 + i: CR columns, type R2
# 2350 + i: I columns,
# 2450 + i: CR columns, type R3

# Block IDs of fuel regions
fuel_blocks = '
1055 1056 1057 1058 1059 1060 1061 1062 1063 1064
1155 1156 1157 1158 1159 1160 1161 1162 1163 1164
1255 1256 1257 1258 1259 1260 1261 1262 1263 1264
1355 1356 1357 1358 1359 1360 1361 1362 1363 1364
'

# fuel columns including the axial reflectors (since cooling will occur there)
full_fuel_blocks = '
1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068
1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168
1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268
1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368
'

fuel_blocks_33pin = '
1055 1056 1057 1058 1059 1060 1061 1062 1063 1064
1155 1156 1157 1158 1159 1160 1161 1162 1163 1164
'

fuel_blocks_31pin = '
1255 1256 1257 1258 1259 1260 1261 1262 1263 1264
1355 1356 1357 1358 1359 1360 1361 1362 1363 1364
'

rr_fuel_blocks_33pin = '
1051 1052 1053 1054 1065 1066 1067 1068
1151 1152 1153 1154 1165 1166 1167 1168
'

rr_fuel_blocks_31pin = '
1251 1252 1253 1254 1265 1266 1267 1268
1351 1352 1353 1354 1365 1366 1367 1368
'

# CR columns including the axial reflectors (since cooling will occur there)
full_cr_blocks = '
2051 2052 2053 2054 2055 2056 2057 2058 2059 2060 2061 2062 2063 2064 2065 2066 2067 2068
2151 2152 2153 2154 2155 2156 2157 2158 2159 2160 2161 2162 2163 2164 2165 2166 2167 2168
2171 2172 2173 2174 2175 2176 2177 2178 2179 2180 2181 2182 2183 2184 2185 2186 2187 2188
2251 2252 2253 2254 2255 2256 2257 2258 2259 2260 2261 2262 2263 2264 2265 2266 2267 2268
2451 2452 2453 2454 2455 2456 2457 2458 2459 2460 2461 2462 2463 2464 2465 2466 2467 2468
'

# Remaining blocks containing the removable, permanent reflectors, and instrumentation
full_pr_blocks = '
1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 1465 1466 1467 1468
1551 1552 1553 1554 1555 1556 1557 1558 1559 1560 1561 1562 1563 1564 1565 1566 1567 1568
2351 2352 2353 2354 2355 2356 2357 2358 2359 2360 2361 2362 2363 2364 2365 2366 2367 2368
'

rpv_blocks = '2'

Tinlet = 453.15 # fluid inlet temperature

rpv_outer_radius = 2.872 # m
vcs_radius = 4.00 # m

stefan_boltzmann_constant = 5.670367e-8
vcs_emissivity = 0.95
rpv_emissivity = 0.95

D_ext_fuel = 4.1e-2 # external diameter in m
npins_fuel33 = 33
npins_fuel31 = 31

D_ext_cr = 1.23e-1 # external diameter in m

p = 0.36 # hexagonal pitch in m
Sblock = '${fparse 0.5 * sqrt(3) * p * p}' # Cross-section area of the hexagonal block in m^2

# volume fractions of graphite for the 31 fuel, 33 fuel and CR blocks
vol_frac_fuel33 = '${fparse 1 - 0.25 * npins_fuel33 * pi * D_ext_fuel * D_ext_fuel / Sblock}'
vol_frac_fuel31 = '${fparse 1 - 0.25 * npins_fuel31 * pi * D_ext_fuel * D_ext_fuel / Sblock}'
vol_frac_cr = '${fparse 1 - 0.25 * 3            * pi * D_ext_cr * D_ext_cr / Sblock}' # use 3 instead of 2 because 3 holes (but only 2 CRs)

# weighting factor of sleeve temperature [see HTTR_XS.xlsx for computation]
# Tmod is defined, in the fuel blocks, as: xi * Tsleeve + (1 - xi) * Tsolid
xi31 = 0.142146 # 31-pin fuel blocks
xi33 = 0.154815 # 33-pin fuel blocks

extension = '_9MW'

[Problem]
  restart_file_base = neutronics_eigenvalue${extension}_moose_modules0_checkpoint_cp/LATEST
  force_restart = true
[]

[Mesh]
  file = neutronics_eigenvalue${extension}_moose_modules0_checkpoint_cp/LATEST
[]

[Variables]
  [Tsolid] # also wall temperature (moderator side) from Griffin to RELAP-7
  []
[]

[Functions]
  [T_init]
    type = ParsedFunction
    expression = '${Tinlet}'
  []
  [T_outside]
    type = ParsedFunction
    # set to 300 K to match temp with VCS operation (concrete bio shield temp)
    expression = 300 # change to piecewise function for transient without VCS
  []
  [decay_heat_power_density_func]
    # x in s
    # y in W/m^3
    type = PiecewiseLinear
    x = '0' # 8.64 86.4 259.2 432 864 1296 1728 2160 2592 3024 3456 3888 4320 8640 12960 17280 21600 25920 30240 34560 38880 43200 47520 51840 56160 60480 64800 69120 73440 77760 82080 86400'
    y = '5.7472e+5' # 4.2455e+5 2.8708e+5 2.3162e+5 2.0972e+5 1.8030e+5 1.6245e+5 1.4962e+5 1.3976e+5 1.3190e+5 1.2547e+5 1.2011e+5 1.1555e+5 1.1163e+5 8.9431e+4 7.9024e+4 7.2585e+4 6.8042e+4 6.4565e+4 6.1753e+4 5.9397e+4 5.7373e+4 5.5606e+4 5.4042e+4 5.2646e+4 5.1388e+4 5.0247e+4 4.9205e+4 4.8248e+4 4.7366e+4 4.6547e+4 4.5785e+4 4.5073e+4'
  []
  [ssteel_304_k]
    type = PiecewiseLinear # [W/mK]
    x = '300.0 400.0 500.0 600.0 800.0 1000.0'
    y = '14.9   16.7  18.3  19.7  22.6  25.4'
  []
  [IG110_k] # thermal conductivity divided by the temperature (dirty trick to make AnisoHeatConductionMaterial define a temperature-independent anistropic thermal conductivity)
    type = ParsedFunction
    symbol_values = '6.632e+01 -4.994e-02 1.712e-05' # [W/m/K]
    symbol_names = 'a0        a1         a2       '
    expression = '((a0 + a1 * t + a2 * t * t) - 1) / t'
  []
  [IG110_cp] # assumed to be the same as for H-451 graphite, taken from MHTGR-350-Appendices.r0.pdf
    type = ParsedFunction
    symbol_values = '0.54212 -2.42667e-6 -90.2725 -43449.3 1.59309e7 -1.43688e9 4184' # [J/kg/K]
    symbol_names = 'a0     a1          a2       a3       a4        a5         b'
    expression = 'if(t < 300, 712.76, (a0 + a1 * t + a2 / t + a3 / t / t + a4 / t / t / t + a5 / t / t / t / t) * b)'
  []
[]

[Kernels]
  [heat_conduction]
    type = HeatConduction
    variable = Tsolid
    block = '${rpv_blocks}'
  []
  [aniso_heat_conduction]
    type = AnisoHeatConduction
    variable = Tsolid
    thermal_conductivity = aniso_thermal_conductivity
    block = '${full_fuel_blocks} ${full_cr_blocks} ${full_pr_blocks}' # all blocks but RPV
  []
  [heat_ie] # time term in heat conduction equation, no need for Steady-state
    type = HeatConductionTimeDerivative
    variable = Tsolid
    block = '${rpv_blocks}'
  []
  [aniso_heat_ie] # time term in heat conduction equation, no need for Steady-state
    type = HeatConductionTimeDerivative
    variable = Tsolid
    block = '${full_fuel_blocks} ${full_cr_blocks} ${full_pr_blocks}' # all blocks but RPV
    specific_heat = aniso_specific_heat
  []

  # conductance kernels (homogenized conduction + radiation through block-sleeve gap)
  [heat_loss_conductance_active_fuel]
    type = Removal
    variable = Tsolid
    block = '${fuel_blocks_31pin} ${fuel_blocks_33pin}'
    sigr = gap_conductance_mat
  []
  [heat_gain_conductance_active_fuel]
    type = MatCoupledForce
    variable = Tsolid
    v = inner_Twall
    block = '${fuel_blocks_31pin} ${fuel_blocks_33pin}'
    material_properties = gap_conductance_mat
  []

  # convection kernels (homogenized convection extracted by fuel and CR cooling channels)
  [heat_loss_convection_outer]
    type = Removal
    variable = Tsolid
    block = '${full_fuel_blocks} ${full_cr_blocks}'
    sigr = Hw_outer_homo_mat
  []
  [heat_gain_convection_outer]
    type = MatCoupledForce
    variable = Tsolid
    v = Tfluid
    block = '${full_fuel_blocks} ${full_cr_blocks}'
    material_properties = Hw_outer_homo_mat
  []
[]

[BCs]
  # for now, do not use BC out of the core to not lose heat (most of it should go back in the core since it pre-cools the helium)

  [RPV_in_BC] # FIXME: this BC adds energy to the RPV that is not removed from the fluid! (but only ~40kW for 9MW so about 0.6K error on the inlet fluid temperature)
    type = ConvectiveFluxFunction # (Robin BC)
    variable = Tsolid
    boundary = '200' # inner RPV
    coefficient = 30.7 # calculated value for 9MW forced convection # old value: 1e2 W/K/m^2 - arbitrary value
    T_infinity = T_init
  []
  [RPV_out_BC] # k \nabla T = h (T- T_inf) at RPV outer boundary
    type = ConvectiveFluxFunction # (Robin BC)
    variable = Tsolid
    boundary = '100' # outer RPV
    coefficient = 2.58 # calculated value for natural convection following NRC report method
    T_infinity = T_outside # matches Figure 13 value concrete shield
  []
  [radiative_BC] # radiation from RPV to bio shield
    type = InfiniteCylinderRadiativeBC
    boundary = '100' # outer RPV
    variable = Tsolid
    boundary_radius = ${rpv_outer_radius}
    boundary_emissivity = ${rpv_emissivity}
    cylinder_radius = ${vcs_radius}
    cylinder_emissivity = ${vcs_emissivity} # matching NRC value, 0.79 for inner RPV
    Tinfinity = T_outside # matches Figure 13 value concrete shield
  []
[]

[ThermalContact]
  [RPV_gap]
    type = GapHeatTransfer
    gap_geometry_type = 'CYLINDER'
    emissivity_primary = 0.8
    emissivity_secondary = 0.8 # varies from 0.6 to 1.0 in RPV paper, 0.79 in NRC paper
    variable = Tsolid # graphite -> rpv gap
    primary = '200'
    secondary = '1'
    gap_conductivity = 0.2091 # for He at 453K, 2.8 MPa
    quadrature = true
    cylinder_axis_point_1 = '0 0 0'
    cylinder_axis_point_2 = '5.22 0 0'
  []
[]

[AuxVariables]
  [power_density]
    block = '${fuel_blocks}'
    order = FIRST
    family = L2_LAGRANGE
  []
  [scaled_initial_power_density] # necessary for the transient to make the decay heat proportional to the initial power
    block = '${fuel_blocks}'
    order = FIRST
    family = L2_LAGRANGE
  []
  [decay_heat_power_density]
    block = '${fuel_blocks}'
    order = FIRST
    family = L2_LAGRANGE
  []
  [heat_source_tot]
    block = '${fuel_blocks}'
    order = FIRST
    family = L2_LAGRANGE
  []
  [Tfuel] # fuel temperature to transfer back to Griffin
    order = CONSTANT
    family = MONOMIAL
  []
  [Tsleeve] # sleeve temperature from BISON (only for fueled blocks)
    block = '${fuel_blocks}'
    order = CONSTANT
    family = MONOMIAL
  []
  [Tmod] # moderator temperature: xi * Tsleeve + (1 - xi) * Tsolid (in the fuel blocks), Tsolid otherwise
    block = '${fuel_blocks} ${rr_fuel_blocks_33pin} ${rr_fuel_blocks_31pin} ${full_cr_blocks} ${full_pr_blocks}'
    order = FIRST
    family = L2_LAGRANGE
  []
  [inner_Twall] # wall temperature (fuel side) from BISON to RELAP-7
    block = '${fuel_blocks}'
    order = CONSTANT
    family = MONOMIAL
  []
  [Tfluid] # fluid temperature from RELAP-7
    order = CONSTANT
    family = MONOMIAL
  []
  [average_axial_tsolid]
    block = '${full_fuel_blocks}' #' ${full_cr_blocks}'
    order = CONSTANT
    family = MONOMIAL
  []
  [Tinit_var]
    block = '${rpv_blocks}'
  []
  [T_inf_outside]
    block = '${rpv_blocks}'
  []
  [RPV_convection_removal]
    block = '${rpv_blocks}'
  []
  [layered_averaged_power_density]
    block = '${fuel_blocks}'
    order = CONSTANT
    family = MONOMIAL
  []
  [heat_balance]
    block = '${full_fuel_blocks} ${full_cr_blocks}'
    order = FIRST
    family = L2_LAGRANGE
  []
  [gap_conductance] # homogenized gap conductance for the inner radius of fuel cooling channels [W/K/m^3]
    block = '${fuel_blocks}'
    order = CONSTANT
    family = MONOMIAL
  []
  [Hw_inner] # heat transfer coefficient for the inner radius of fuel cooling channels [W/K/m^2] (not homogenized because only meant to be transferred to BISON)
    block = '${fuel_blocks}'
    order = CONSTANT
    family = MONOMIAL
  []
  [Hw_outer] # heat transfer coefficient for the inner radius of fuel cooling channels [W/K/m^2] (not homogenized because only meant to be transferred to BISON)
    block = '${full_fuel_blocks}' # No need to include the CR blocks because only meant to be transferred to BISON for a homogeneous full-core heat transfer model
    order = CONSTANT
    family = MONOMIAL
  []
  [Hw_outer_homo] # homogenized heat transfer coefficient for the outer radius of cooling channels [W/K/m^3]
    # Homogenization performed in the RELAP inputs  because scaling depends on the number of pins
    block = '${full_fuel_blocks} ${full_cr_blocks}'
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [assign_initial_power_density]
    type = NormalizationAux
    block = '${fuel_blocks}'
    variable = scaled_initial_power_density
    source_variable = power_density
    normalization = ptot
    execute_on = 'initial'
  []
  [assign_decay_heat_power_density]
    type = ScaleAux
    block = '${fuel_blocks}'
    variable = decay_heat_power_density
    source_variable = scaled_initial_power_density
    multiplying_pp = decay_heat_pp
    execute_on = 'initial timestep_begin timestep_end'
  []
  [assign_heat_source_tot]
    type = ParsedAux
    block = '${fuel_blocks}'
    variable = heat_source_tot
    expression = 'power_density + decay_heat_power_density'
    coupled_variables = 'power_density decay_heat_power_density'
    execute_on = 'initial timestep_begin timestep_end'
  []
  [assign_T_inf_outside] # needs to be time dependent for the case without VCS
    type = FunctionAux
    block = '${rpv_blocks}'
    variable = T_inf_outside
    function = T_outside
    execute_on = 'initial timestep_begin timestep_end'
  []
  [assign_average_axial_tsolid]
    type = SpatialUserObjectAux
    variable = average_axial_tsolid
    user_object = average_axial_temperature_UO
    block = '${full_fuel_blocks}' #' ${full_cr_blocks}'
    execute_on = 'timestep_end'
  []
  [assign_layered_averaged_power_density]
    type = SpatialUserObjectAux
    variable = layered_averaged_power_density
    user_object = average_power_UO
    block = '${fuel_blocks}'
    execute_on = 'timestep_end'
  []

  [heat_balance_active_fuel]
    type = ParsedAux
    variable = heat_balance
    block = '${fuel_blocks}'
    expression = 'gap_conductance * (inner_Twall - Tsolid) + Hw_outer_homo * (Tfluid - Tsolid)'
    coupled_variables = 'gap_conductance inner_Twall Tsolid Hw_outer_homo Tfluid'
  []
  [heat_balance_convection_rr_cr]
    type = ParsedAux
    variable = heat_balance
    block = '${rr_fuel_blocks_31pin} ${rr_fuel_blocks_33pin} ${full_cr_blocks}'
    expression = 'Hw_outer_homo * (Tfluid - Tsolid)'
    coupled_variables = 'Hw_outer_homo Tfluid Tsolid'
  []

  [RPV_convection_removal_aux]
    type = ParsedAux
    variable = RPV_convection_removal
    block = '${rpv_blocks}'
    expression = '${stefan_boltzmann_constant} * ${rpv_emissivity} * ${vcs_emissivity} * ${vcs_radius} / (${vcs_emissivity} * ${vcs_radius} + ${rpv_emissivity} * ${rpv_outer_radius} * (1 - ${vcs_emissivity})) * (pow(Tsolid, 4) - pow(T_inf_outside,4))'
    coupled_variables = 'T_inf_outside Tsolid'
  []

  [assign_Tmod_fuelblocks31]
    type = ParsedAux
    variable = Tmod
    block = '${fuel_blocks_31pin}'
    expression = '${xi31} * Tsleeve + (1 - ${xi31}) * Tsolid'
    coupled_variables = 'Tsleeve Tsolid'
    execute_on = 'timestep_end'
  []
  [assign_Tmod_fuelblocks33]
    type = ParsedAux
    variable = Tmod
    block = '${fuel_blocks_33pin}'
    expression = '${xi33} * Tsleeve + (1 - ${xi33}) * Tsolid'
    coupled_variables = 'Tsleeve Tsolid'
    execute_on = 'timestep_end'
  []
  [assign_Tmod_elsewhere]
    type = ParsedAux
    variable = Tmod
    block = ' ${rr_fuel_blocks_33pin} ${rr_fuel_blocks_31pin} ${full_cr_blocks} ${full_pr_blocks}'
    expression = 'Tsolid'
    coupled_variables = 'Tsolid'
    execute_on = 'timestep_end'
  []
[]

[UserObjects]
  [average_power_UO]
    type = NearestPointLayeredAverage
    block = '${fuel_blocks}'
    variable = heat_source_tot
    direction = x
    points_file = '../mesh/centers_fuel_blocks.txt' # fuel column centers
    bounds = '1.16 1.74 2.32 2.90 3.48 4.06' # axial boundaries of each block
    execute_on = 'timestep_begin timestep_end'
  []
  [average_Tsolid_UO]
    type = NearestPointLayeredAverage
    block = '${full_fuel_blocks} ${full_cr_blocks}'
    variable = Tsolid
    direction = x
    points_file = '../mesh/centers_relap.txt' # fuel and CR column centers
    bounds = '0 0.58 1.16 1.74 2.32 2.90 3.48 4.06 4.64 5.22' # axial boundaries of each block
    execute_on = 'timestep_begin timestep_end'
  []
  [average_inner_Twall_UO]
    type = NearestPointLayeredAverage
    block = '${fuel_blocks}'
    variable = inner_Twall
    direction = x
    points_file = '../mesh/centers_fuel_blocks.txt' # fuel and CR column centers
    bounds = '0 0.58 1.16 1.74 2.32 2.90 3.48 4.06 4.64 5.22' # axial boundaries of each block
    execute_on = 'timestep_begin timestep_end'
  []
  [average_htc_inner_UO]
    type = NearestPointLayeredAverage
    block = '${fuel_blocks}'
    variable = Hw_inner
    direction = x
    points_file = '../mesh/centers_fuel_blocks.txt' # fuel and CR column centers
    bounds = '0 0.58 1.16 1.74 2.32 2.90 3.48 4.06 4.64 5.22' # axial boundaries of each block
    execute_on = 'timestep_begin timestep_end'
  []
  [average_htc_outer_UO]
    type = NearestPointLayeredAverage
    block = '${full_fuel_blocks}'
    variable = Hw_outer
    direction = x
    points_file = '../mesh/centers_fuel_blocks.txt' # fuel and CR column centers
    bounds = '0 0.58 1.16 1.74 2.32 2.90 3.48 4.06 4.64 5.22' # axial boundaries of each block
    execute_on = 'timestep_begin timestep_end'
  []
  [average_Tfluid_UO]
    type = NearestPointLayeredAverage
    block = '${full_fuel_blocks} ${full_cr_blocks}'
    variable = Tfluid
    direction = x
    points_file = '../mesh/centers_relap.txt' # fuel and CR column centers
    bounds = '0 0.58 1.16 1.74 2.32 2.90 3.48 4.06 4.64 5.22' # axial boundaries of each block
    execute_on = 'timestep_begin timestep_end'
  []
  [average_axial_temperature_UO] # for postprocessing only (e.g. to get an idea of the axial average temperature)
    type = LayeredAverage
    block = '${full_fuel_blocks}' #' ${full_cr_blocks}'
    variable = Tsolid
    direction = x
    bounds = '0 0.145 0.29 0.435 0.58 0.725 0.87 1.015 1.16 1.305 1.45 1.595 1.74 1.885 2.03 2.175 2.32 2.465 2.61 2.755 2.9 3.045 3.19 3.335 3.48 3.625 3.77 3.915 4.06 4.205 4.35 4.495 4.64 4.785 4.93 5.075 5.22' # for 36 axial elements
    execute_on = 'timestep_begin timestep_end'
  []
[]

[Materials]
  [graphite_moderator]
    type = AnisoHeatConductionMaterial
    block = '${full_fuel_blocks} ${full_cr_blocks} ${full_pr_blocks}' # all blocks but RPV
    temperature = Tsolid
    base_name = 'aniso'
    reference_temperature = 0.0
    thermal_conductivity_temperature_coefficient_function = IG110_k
    thermal_conductivity = '1    0    0
                            0 0.54    0
                            0    0 0.54' # isotropic for testing only, otherwise one in the x-direction
    specific_heat = IG110_cp
  []
  # The base graphite densities come from table 1.27 of the NEA/NSC/DOC(2006)1 report and are then corrected with the volume fraction of graphite in each block
  [graphite_density_fuel33]
    type = Density
    block = '${fuel_blocks_33pin}'
    density = '${fparse 1770 * vol_frac_fuel33}'
  []
  [graphite_density_fuel31]
    type = Density
    block = '${fuel_blocks_31pin}'
    density = '${fparse 1770 * vol_frac_fuel31}'
  []
  [graphite_density_rr_fuel33]
    type = Density
    block = '${rr_fuel_blocks_33pin}'
    density = '${fparse 1760 * vol_frac_fuel33}'
  []
  [graphite_density_rr_fuel31]
    type = Density
    block = '${rr_fuel_blocks_31pin}'
    density = '${fparse 1760 * vol_frac_fuel31}'
  []
  [graphite_density_cr]
    type = Density
    block = '${full_cr_blocks}'
    density = '${fparse 1770 * vol_frac_cr}'
  []
  [graphite_density_rest]
    type = Density
    block = '${full_pr_blocks}'
    density = '${fparse 1760}' # In reality, the actual PR blocks are made of PGX graphite of density 1.732 g/cc but here we assume IG110 graphite for now.
  []
  [ss304_rpv]
    type = HeatConductionMaterial
    block = '${rpv_blocks}'
    thermal_conductivity_temperature_function = ssteel_304_k
    specific_heat = 500 # J/kg/K
    temp = Tsolid
  []
  [ss304_density]
    type = Density
    block = '${rpv_blocks}'
    density = 8000 # typical value for ss304 in kg/m^3
  []
  [gap_conductance_mat]
    type = Variable2Material
    block = '${fuel_blocks}'
    variables = 'gap_conductance'
    prop_names = 'gap_conductance_mat'
  []
  [Hw_outer_mat]
    type = Variable2Material
    block = '${full_fuel_blocks} ${full_cr_blocks}'
    variables = Hw_outer_homo
    prop_names = Hw_outer_homo_mat
  []
[]

[Postprocessors]
  [Tsolid_avg]
    type = ElementAverageValue
    variable = Tsolid
    block = '${full_fuel_blocks} ${full_cr_blocks}'
    execute_on = 'initial linear'
  []
  [Tsolid_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tsolid
    block = '${full_fuel_blocks} ${full_cr_blocks}'
    execute_on = 'initial linear'
  []
  [Tsolid_min]
    type = ElementExtremeValue
    value_type = min
    variable = Tsolid
    block = '${full_fuel_blocks} ${full_cr_blocks}'
    execute_on = 'initial linear'
  []
  [Tsolid_fuel_avg]
    type = ElementAverageValue
    variable = Tsolid
    block = '${fuel_blocks}'
    execute_on = 'initial linear'
  []
  [Tsolid_fuel_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tsolid
    block = '${fuel_blocks}'
    execute_on = 'initial linear'
  []
  [Tsolid_fuel_min]
    type = ElementExtremeValue
    value_type = min
    variable = Tsolid
    block = '${fuel_blocks}'
    execute_on = 'initial linear'
  []
  [Tsolid_pr_avg]
    type = ElementAverageValue
    variable = Tsolid
    block = '${full_pr_blocks}'
    execute_on = 'initial linear'
  []
  [Tsolid_pr_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tsolid
    block = '${full_pr_blocks}'
    execute_on = 'initial linear'
  []
  [Tsolid_pr_min]
    type = ElementExtremeValue
    value_type = min
    variable = Tsolid
    block = '${full_pr_blocks}'
    execute_on = 'initial linear'
  []
  [Tfuel_avg]
    type = ElementAverageValue
    variable = Tfuel
    block = '${fuel_blocks}'
    execute_on = 'initial linear'
  []
  [Tfuel_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tfuel
    block = '${fuel_blocks}'
    execute_on = 'initial linear'
  []
  [Tfuel_min]
    type = ElementExtremeValue
    value_type = min
    variable = Tfuel
    block = '${fuel_blocks}'
    execute_on = 'initial linear'
  []
  [Tsleeve_avg]
    type = ElementAverageValue
    variable = Tsleeve
    block = '${fuel_blocks}'
    execute_on = 'initial linear'
  []
  [Tsleeve_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tsleeve
    block = '${fuel_blocks}'
    execute_on = 'initial linear'
  []
  [Tsleeve_min]
    type = ElementExtremeValue
    value_type = min
    variable = Tsleeve
    block = '${fuel_blocks}'
    execute_on = 'initial linear'
  []
  [Tmod_avg]
    type = ElementAverageValue
    variable = Tmod
    block = '${fuel_blocks}'
    execute_on = 'initial linear'
  []
  [Tmod_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tmod
    block = '${fuel_blocks}'
    execute_on = 'initial linear'
  []
  [Tmod_min]
    type = ElementExtremeValue
    value_type = min
    variable = Tmod
    block = '${fuel_blocks}'
    execute_on = 'initial linear'
  []
  [Trpv_avg]
    type = ElementAverageValue
    variable = Tsolid
    block = '${rpv_blocks}'
    execute_on = 'initial linear'
  []
  [Trpv_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tsolid
    block = '${rpv_blocks}'
    execute_on = 'initial linear'
  []
  [Trpv_min]
    type = ElementExtremeValue
    value_type = min
    variable = Tsolid
    block = '${rpv_blocks}'
    execute_on = 'initial linear'
  []
  [ptot]
    type = ElementIntegralVariablePostprocessor
    variable = power_density
    block = '${fuel_blocks}'
    execute_on = 'initial linear'
  []
  [ptot_decay]
    type = ElementIntegralVariablePostprocessor
    variable = decay_heat_power_density
    block = '${fuel_blocks}'
    execute_on = 'initial linear'
  []
  [decay_heat_pp]
    type = FunctionValuePostprocessor
    function = decay_heat_power_density_func
    execute_on = 'initial timestep_begin timestep_end'
  []
  [pdens_avg]
    type = ElementAverageValue
    variable = power_density
    block = '${fuel_blocks}'
    execute_on = 'initial linear'
  []
  [net_heat_src_avg]
    type = ElementAverageValue
    variable = heat_balance
    block = '${full_fuel_blocks} ${full_cr_blocks}'
    execute_on = 'initial linear'
  []
  [net_heat_src_max]
    type = ElementExtremeValue
    value_type = max
    variable = heat_balance
    block = '${full_fuel_blocks} ${full_cr_blocks}'
    execute_on = 'initial linear'
  []
  [net_heat_src_min]
    type = ElementExtremeValue
    value_type = min
    variable = heat_balance
    block = '${full_fuel_blocks} ${full_cr_blocks}'
    execute_on = 'initial linear'
  []
  [net_heat_src_tot]
    type = ElementIntegralVariablePostprocessor
    variable = heat_balance
    block = '${full_fuel_blocks} ${full_cr_blocks}'
    execute_on = 'initial linear'
  []
  [Tfluid_avg]
    type = ElementAverageValue
    variable = Tfluid
    block = '${full_fuel_blocks} ${full_cr_blocks}'
    execute_on = 'initial linear'
  []
  [Tfluid_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tfluid
    block = '${full_fuel_blocks} ${full_cr_blocks}'
    execute_on = 'initial linear'
  []
  [rpv_convective_in]
    type = ConvectiveHeatTransferSideIntegral
    T_solid = Tsolid
    boundary = '200'
    T_fluid_var = Tinit_var
    htc = 30.7 # calculated value for 9MW forced convection
  []
  [rpv_convective_out]
    type = ConvectiveHeatTransferSideIntegral
    T_solid = Tsolid
    boundary = '100'
    T_fluid_var = T_inf_outside
    htc = 2.58 # calculated value for natural convection following NRC report method
  []
  [rpv_radiative_out]
    type = SideIntegralVariablePostprocessor
    variable = RPV_convection_removal
    boundary = '100'
  []
  [bison_total_power]
    type = Receiver
  []
  [bison_total_bdy_heat_flux]
    type = Receiver
  []
  [Tfluid_out]
    type = Receiver
  []
  [delta_H]
    type = Receiver
  []
  [gap_cond_max]
    type = ElementExtremeValue
    value_type = max
    variable = gap_conductance
    block = '${fuel_blocks}'
    execute_on = 'initial linear'
  []
  [gap_cond_min]
    type = ElementExtremeValue
    value_type = min
    variable = gap_conductance
    block = '${fuel_blocks}'
    execute_on = 'initial linear'
  []
[]

[Executioner]
  type = Transient
  automatic_scaling = true
  compute_scaling_once = false
  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart '
  petsc_options_value = 'hypre boomeramg 100'

  start_time = 0
  end_time = 50
  dt = 10

  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-9
  fixed_point_rel_tol = 1e-3
  fixed_point_abs_tol = 1e-7

  line_search = none # default seems bad for convergence
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    positions_file = '../mesh/centers_relap_33pins.txt
                      ../mesh/centers_relap_31pins.txt'
    input_files = 'fuel_elem_null.i'
    cli_args = 'npins=33 npins=33 npins=33 npins=33 npins=33 npins=33 npins=33 npins=33 npins=33 npins=33 npins=33 npins=33
                npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31'
    execute_on = 'timestep_begin'
    output_in_position = true
  []
  [relap]
    type = TransientMultiApp
    positions_file = '../mesh/centers_relap_33pins.txt
                      ../mesh/centers_relap_31pins.txt
                      ../mesh/centers_relap_CR.txt'

    # the first 12 positions are fuel with 33 pins, the next 18 are fuel with 31 pins, the last 16 are CR
    input_files = 'thermal_hydraulics_fuel_pins_null.i
                   thermal_hydraulics_fuel_pins_null.i
                   thermal_hydraulics_CR_null.i'

    cli_args = 'npins=33 npins=33 npins=33 npins=33 npins=33 npins=33 npins=33 npins=33 npins=33 npins=33 npins=33 npins=33
                npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31 npins=31
                npins=2  npins=2  npins=2  npins=2  npins=2  npins=2  npins=2  npins=2  npins=2  npins=2  npins=2  npins=2  npins=2  npins=2  npins=2  npins=2'

    execute_on = 'timestep_end'
    max_procs_per_app = 1
    sub_cycling = true
  []
[]

[Transfers]
  [pdens_to_bison]
    type = MultiAppUserObjectTransfer
    to_multi_app = bison
    user_object = average_power_UO
    variable = power_density
  []
  [tmod_to_bison]
    type = MultiAppUserObjectTransfer
    to_multi_app = bison
    user_object = average_Tsolid_UO
    variable = Tmod
  []
  [tfluid_to_bison]
    type = MultiAppInterpolationTransfer
    to_multi_app = bison
    source_variable = Tfluid
    variable = Tfluid
  []
  [htc_inner_to_bison]
    type = MultiAppInterpolationTransfer
    to_multi_app = bison
    source_variable = Hw_inner
    variable = Hw_inner
  []
  [htc_outer_to_bison]
    type = MultiAppInterpolationTransfer
    to_multi_app = bison
    source_variable = Hw_outer
    variable = Hw_outer
  []
  [inner_twall_to_relap]
    type = MultiAppInterpolationTransfer
    to_multi_app = relap
    source_variable = inner_Twall
    variable = T_wall:1
  []
  [tmod_to_relap]
    type = MultiAppUserObjectTransfer
    to_multi_app = relap
    user_object = average_Tsolid_UO
    variable = T_wall:2
  []
  [tmod_to_relap_topbottom]
    type = MultiAppUserObjectTransfer
    to_multi_app = relap
    user_object = average_Tsolid_UO
    variable = T_wall
  []
  [tfuel_from_bison]
    type = MultiAppUserObjectTransfer
    from_multi_app = bison
    user_object = average_Tfuel_UO
    variable = Tfuel
    nearest_sub_app = true
  []
  [tsleeve_from_bison]
    type = MultiAppUserObjectTransfer
    from_multi_app = bison
    user_object = average_Tsleeve_UO
    variable = Tsleeve
    nearest_sub_app = true
  []
  [twall_from_bison]
    type = MultiAppUserObjectTransfer
    from_multi_app = bison
    user_object = inner_wall_temp_UO
    variable = inner_Twall
    nearest_sub_app = true
  []
  [heat_source_from_bison]
    type = MultiAppUserObjectTransfer
    from_multi_app = bison
    user_object = gap_conductance_UO
    variable = gap_conductance
    nearest_sub_app = true
  []
  [total_power_from_bison]
    type = MultiAppPostprocessorTransfer
    from_multi_app = bison
    from_postprocessor = 'fuel_block_total_power'
    to_postprocessor = 'bison_total_power'
    reduction_type = sum
  []
  [bison_heat_flux]
    type = MultiAppPostprocessorTransfer
    from_multi_app = bison
    from_postprocessor = 'bdy_heat_flux_tot'
    to_postprocessor = 'bison_total_bdy_heat_flux'
    reduction_type = sum
  []
  [tfluid_from_relap]
    type = MultiAppUserObjectTransfer
    from_multi_app = relap
    variable = Tfluid
    user_object = avg_Tfluid_UO
    nearest_sub_app = true
  []
  [hw_inner_from_relap]
    type = MultiAppUserObjectTransfer
    from_multi_app = relap
    variable = Hw_inner
    user_object = avg_Hw_inner_UO
    nearest_sub_app = true
  []
  [hw_outer_from_relap]
    type = MultiAppUserObjectTransfer
    from_multi_app = relap
    variable = Hw_outer
    user_object = avg_Hw_outer_UO
    nearest_sub_app = true
  []
  [hw_outer_homo_from_relap]
    type = MultiAppUserObjectTransfer
    from_multi_app = relap
    variable = Hw_outer_homo
    user_object = avg_Hw_outer_homo_UO
    nearest_sub_app = true
  []
  [tout_from_relap]
    type = MultiAppPostprocessorTransfer
    from_multi_app = relap
    from_postprocessor = Tout_weighted
    to_postprocessor = Tfluid_out
    reduction_type = sum
  []
  [deltaH_from_relap]
    type = MultiAppPostprocessorTransfer
    from_multi_app = relap
    from_postprocessor = delta_H
    to_postprocessor = delta_H
    reduction_type = sum
  []
[]

[Outputs]
  [exodus]
    type = Exodus
    overwrite = true
  []
  csv = true
[]
(htgr/httr/steady_state_and_null_transient/full_core_ht_null.i)

BISON 5-pin model

# ==============================================================================
# Fuel element heat transfer model (null-transient)
# Application: BISON (Sabertooth with child applications)
# ------------------------------------------------------------------------------
# Idaho Falls, INL, April 2023
# Author(s): Vincent Laboure
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
# ==============================================================================

fuel_blocks = '1'
sleeve_blocks = '2'
mod_blocks = '4'

h = 0.577 # height of a block in m
p = 0.36 # hexagonal pitch in m
Sblock = '${fparse 0.5 * sqrt(3) * p * p}' # Cross-section area of the hexagonal block in m^2
Vblock = '${fparse h * Sblock}' # volume of a block in m^3

npins = 31 # number of fuel channels in the current block
h_fuel = '${fparse 0.56 - 0.014}' # height of a fuel pin in m
inner_radius = 0.005 # inner radius of a fuel pin in m
outer_radius = 0.013 # outer radius of a fuel pin in m
Vfuel = '${fparse npins * h_fuel * pi * (outer_radius * outer_radius - inner_radius * inner_radius)}' # volume of fuel compacts in a block

multiplier = '${fparse Vblock / Vfuel}' # multiplier to obtain the scaled power density. Roughly 8.0 for 33 pins, 8.5 for 31 pins [units = none]

sleeve_outer_radius = 0.017 # outer radius of the sleeve in m
inner_outer_radius = 0.0205 # outer radius of the mesh in m

htc_homo_scaling = '${fparse 2 * pi * inner_outer_radius * npins / Sblock}'

stefan_boltzmann_constant = 5.670367e-8
eps1 = 0.8
eps2 = 0.8
cond = 0.2309 # He at 523K, 2.8MPa

[Mesh]
  [fmg]
    type = FileMeshGenerator
    file = '../mesh/fuel_element/HTTR_fuel_pin_2D_refined_m_5pins_axialref.e'
  []
  uniform_refine = 0 # do not modify if using ThermalContact
  coord_type = RZ
  rz_coord_axis = X
[]

[Variables]
  [temp]
  []
[]

[AuxVariables]
  [power_density]
    block = '${fuel_blocks}'
    order = CONSTANT
    family = MONOMIAL
  []
  [scaled_power_density]
    block = '${fuel_blocks}'
    order = CONSTANT
    family = MONOMIAL
  []
  [Tmod]
    order = CONSTANT
    family = MONOMIAL
  []
  [Tfluid] # fluid temperature from RELAP-7
    order = CONSTANT
    family = MONOMIAL
  []
  [Tfuel] # Tfuel variable for XS interpolation, no longer includes the sleeve!
    block = '${fuel_blocks}'
  []
  [Tsleeve] # Tsleeve variable for XS interpolation (weighted summation with graphite block temperature)
    block = '${sleeve_blocks}'
  []
  [inner_Twall]
    block = '${fuel_blocks} ${sleeve_blocks}'
    order = CONSTANT
    family = MONOMIAL
  []
  [bdy_heat_flux_aux]
    order = FIRST # not CONSTANT MONOMIAL! (or weird effects when evaluating the boundary integral) but has to be discontinuous
    family = L2_LAGRANGE
  []
  [bdy_heat_flux_layered_integral]
    order = CONSTANT
    family = MONOMIAL
  []
  [gap_conductance]
    block = '${fuel_blocks} ${sleeve_blocks}' # because we only want it defined over the active fuel region for proper averaging (won't change radially)
    order = CONSTANT
    family = MONOMIAL
  []
  [Hw_inner] # heat transfer coefficient for the inner radius of fuel cooling channels [W/K/m^2]
    block = '${fuel_blocks} ${sleeve_blocks}' # because we only want it defined over the active fuel region for proper averaging (won't change radially)
    order = CONSTANT
    family = MONOMIAL
  []
  [Hw_outer] # heat transfer coefficient for the outer radius of fuel cooling channels [W/K/m^2]
    block = '${mod_blocks} ${fuel_blocks} ${sleeve_blocks}' # also define on fuel pins because needed for equivalent gap conductance
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temp
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temp
  []
  [heat_source]
    type = CoupledForce
    variable = temp
    block = '${fuel_blocks}'
    v = scaled_power_density
  []
[]

[AuxKernels]
  [GetPowerDensity]
    type = ScaleAux
    variable = scaled_power_density
    source_variable = power_density
    multiplier = ${multiplier}
  []
  [setTfuel]
    type = CoupledAux
    variable = Tfuel
    coupled = temp
    block = '${fuel_blocks}'
  []
  [setTsleeve]
    type = CoupledAux
    variable = Tsleeve
    coupled = temp
    block = '${sleeve_blocks}'
  []
  [set_inner_Twall]
    type = SpatialUserObjectAux
    variable = inner_Twall
    user_object = inner_wall_temp_UO
    block = '${fuel_blocks} ${sleeve_blocks}'
    execute_on = 'nonlinear timestep_end'
  []
  [set_bdy_heat_flux_aux]
    type = ParsedAux
    variable = bdy_heat_flux_aux
    block = '${mod_blocks}'
    expression = '1e5 * (temp - Tmod)'
    coupled_variables = 'temp Tmod'
    execute_on = 'nonlinear timestep_end'
  []
  [set_bdy_heat_flux_layered_integral]
    type = SpatialUserObjectAux
    variable = bdy_heat_flux_layered_integral
    user_object = bdy_heat_flux_UO
    execute_on = 'nonlinear timestep_end'
  []
  [set_gap_conductance]
    type = ParsedAux
    variable = gap_conductance
    block = '${fuel_blocks} ${sleeve_blocks}' # because we only want it defined over the active fuel region for proper averaging (won't change radially)
    expression = '(${cond} / ${inner_outer_radius} / log(${inner_outer_radius} / ${sleeve_outer_radius}) + ${stefan_boltzmann_constant} * (inner_Twall * inner_Twall + Tmod * Tmod) * (inner_Twall + Tmod) / (1 / ${eps1} + 1 / ${eps2} - 1) * 0.5 * (${inner_outer_radius} + ${sleeve_outer_radius}) / ${inner_outer_radius}) * ${htc_homo_scaling}'
    coupled_variables = 'Tmod inner_Twall'
    execute_on = 'nonlinear timestep_end'
  []
[]

[ThermalContact]
  [gap_top]
    type = GapHeatTransfer
    gap_geometry_type = 'PLATE'
    emissivity_primary = 0.8 # MHTGR: 0.85, Shi PhD: 0.8
    emissivity_secondary = 0.8 # MHTGR: 0.85, Shi PhD: 0.8
    # coarser side should be primary
    # (in theory using quadrature = true should make it independent but other errors appear and not exactly conservative)
    primary = Sleeve_top
    secondary = Fuel_top
    gap_conductivity = 0.2735 # He at 668K, 2.8 MPa - Original: 1e-5 W/m/K, Shi PhD: 0.35 (but no gap on the inside)
    variable = temp
    quadrature = true
  []

  [gap_side]
    type = GapHeatTransfer
    gap_geometry_type = 'CYLINDER'
    emissivity_primary = 0.8 # MHTGR: 0.85, Shi PhD: 0.8
    emissivity_secondary = 0.8 # MHTGR: 0.85, Shi PhD: 0.8
    primary = Sleeve_side
    secondary = Fuel_side
    gap_conductivity = 0.2735 # He at 668K, 2.8 MPa - Original: 1e-5 W/m/K, Shi PhD: 0.35 (but no gap on the inside)
    variable = temp
    quadrature = true
  []

  # only use during transient (numerical issues during steady-state and amount of energy through gap should be very small in steady-state)
  [cooling_channel]
    type = GapHeatTransfer
    gap_geometry_type = 'CYLINDER'
    emissivity_primary = 0.8 # MHTGR: 0.85, Shi PhD: 0.8
    emissivity_secondary = 0.8 # MHTGR: 0.85, Shi PhD: 0.8
    primary = Inner_outer_ring
    secondary = Relap7_Tw
    gap_conductivity = 0.2309 # He at 523K, 2.8MPa
    variable = temp
    quadrature = true
  []
[]

[BCs]
  [Neumann]
    type = NeumannBC
    variable = temp
    value = 0.0
    boundary = 'Neumann Sleeve_bot'
  []

  [outside_bc]
    type = CoupledConvectiveHeatFluxBC
    variable = temp
    boundary = 'Outer_outer_ring'
    T_infinity = Tmod
    htc = 1e5
  []

  # For RELAP-7
  [convective_inner]
    type = CoupledConvectiveHeatFluxBC
    boundary = 'Relap7_Tw'
    variable = temp
    T_infinity = Tfluid
    htc = Hw_inner
  []
  [convective_outer]
    type = CoupledConvectiveHeatFluxBC
    boundary = 'Inner_outer_ring'
    variable = temp
    T_infinity = Tfluid
    htc = Hw_outer
  []
[]

[Functions]
  # Functions for the compact k and Cp (and rho) are obtained using the tpmain Fortran script with a packing fraction of 0.3, burnup = 0, fluence = 0
  # FIXME: recompute with real burnup and fluence
  [compact_k]
    type = ParsedFunction
    symbol_values = '-4.438e-09 2.569e-05 -5.186e-02 5.111e+01' # [W/m/K]
    symbol_names = 'a3         a2         a1        a0'
    expression = 'a0 + a1 * t + a2 * t * t + a3 * t * t * t'
  []
  [compact_cp] # assumed to be the same as for H-451 graphite, taken from MHTGR-350-Appendices.r0.pdf
    type = ParsedFunction
    symbol_values = '-2.408e-10 1.468e-06 -3.379e-03 3.654e+00 -2.875e+02' # [J/kg/K]
    symbol_names = 'a4         a3        a2         a1        a0'
    expression = 'a0 + a1 * t + a2 * t * t + a3 * t * t * t + a4 * t * t * t * t'
  []
  [IG110_k]
    type = ParsedFunction
    symbol_values = '6.632e+01 -4.994e-02 1.712e-05' # [W/m/K]
    symbol_names = 'a0        a1         a2'
    expression = 'a0 + a1 * t + a2 * t * t'
  []
  [IG110_cp] # assumed to be the same as for H-451 graphite, taken from MHTGR-350-Appendices.r0.pdf
    type = ParsedFunction
    symbol_values = '0.54212 -2.42667e-6 -90.2725 -43449.3 1.59309e7 -1.43688e9 4184' # [J/kg/K]
    symbol_names = 'a0     a1          a2       a3       a4        a5         b'
    expression = 'if(t < 300, 712.76, (a0 + a1 * t + a2 / t + a3 / t / t + a4 / t / t / t + a5 / t / t / t / t) * b)'
  []
[]

[Materials]
  [Fuel_compact]
    type = HeatConductionMaterial
    block = '${fuel_blocks}'
    temp = temp
    thermal_conductivity_temperature_function = compact_k
    specific_heat_temperature_function = compact_cp
  []
  [Fuel_sleeve]
    type = HeatConductionMaterial
    block = '${sleeve_blocks} ${mod_blocks}'
    temp = temp
    thermal_conductivity_temperature_function = IG110_k
    specific_heat_temperature_function = IG110_cp
  []
  [Fuel_compact_density]
    type = Density
    block = '${fuel_blocks}'
    density = 2573 # kg/m^3 - obtained using the tpmain Fortran script with a packing fraction of 0.3
  []
  [Fuel_sleeve_density]
    type = Density
    block = '${sleeve_blocks} ${mod_blocks}'
    density = 1770 # kg/m^3 - comes from table 1.27 of the NEA/NSC/DOC(2006)1 report
  []
[]

[Executioner]
  type = Transient

  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-6

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu       superlu_dist'

  start_time = 0
  end_time = 50
  dt = 10
[]

[UserObjects]
  [inner_wall_temp_UO]
    type = LayeredSideAverage
    boundary = 'Relap7_Tw'
    variable = temp
    direction = x
    bounds = '1.16 1.74 2.32 2.9 3.48 4.06' # for 9 axial elements
    execute_on = 'nonlinear timestep_end'
  []
  [bdy_heat_flux_UO]
    type = LayeredSideIntegral
    boundary = 'Outer_outer_ring'
    variable = bdy_heat_flux_aux
    direction = x
    bounds = '0 0.58 1.16 1.74 2.32 2.9 3.48 4.06 4.64 5.22' # for 9 axial elements
    execute_on = 'nonlinear timestep_end'
  []
  [gap_conductance_UO]
    type = LayeredAverage
    block = '${fuel_blocks} ${sleeve_blocks}'
    variable = gap_conductance
    direction = x
    bounds = '0 0.58 1.16 1.74 2.32 2.9 3.48 4.06 4.64 5.22' # for 9 axial elements
    execute_on = 'nonlinear timestep_end'
  []
  [average_Tfuel_UO]
    type = LayeredAverage
    block = '${fuel_blocks}'
    variable = Tfuel
    direction = x
    bounds = '0 0.58 1.16 1.74 2.32 2.9 3.48 4.06 4.64 5.22' # for 9 axial elements
    execute_on = 'nonlinear timestep_end'
  []
  [average_Tsleeve_UO]
    type = LayeredAverage
    block = '${sleeve_blocks}'
    variable = Tsleeve
    direction = x
    bounds = '0 0.58 1.16 1.74 2.32 2.9 3.48 4.06 4.64 5.22' # for 9 axial elements
    execute_on = 'nonlinear timestep_end'
  []
[]

[Postprocessors]
  [fuel_k_average]
    type = ElementAverageMaterialProperty
    block = '${fuel_blocks}'
    mat_prop = 'thermal_conductivity'
  []
  [graphite_k_average]
    type = ElementAverageMaterialProperty
    block = '${sleeve_blocks} ${mod_blocks}'
    mat_prop = 'thermal_conductivity'
  []
  [TotalPower]
    type = ElementIntegralVariablePostprocessor
    block = '${fuel_blocks}'
    variable = scaled_power_density
    execute_on = 'initial timestep_end'
  []
  [fuel_block_total_power]
    type = ParsedPostprocessor
    function = '${npins} * TotalPower'
    pp_names = 'TotalPower'
    execute_on = 'initial timestep_end'
  []
  [avg_homo_gap_conductance]
    type = ElementAverageValue
    variable = gap_conductance
    block = '${fuel_blocks} ${sleeve_blocks}'
    execute_on = 'initial timestep_end'
  []
  [avg_gap_conductance]
    type = ParsedPostprocessor
    function = 'avg_homo_gap_conductance / ${htc_homo_scaling}'
    pp_names = 'avg_homo_gap_conductance'
    execute_on = 'initial timestep_end'
  []
  [avg_htc_inner]
    type = ElementAverageValue
    variable = Hw_inner
    block = '${fuel_blocks} ${sleeve_blocks}'
    execute_on = 'initial timestep_end'
  []
  [avg_htc_outer]
    type = ElementAverageValue
    variable = Hw_outer
    block = '${mod_blocks}' #  True average should be computed on mod_blocks to avoid bias towards active fuel regions
    execute_on = 'initial timestep_end'
  []
  [avg_Tfuel]
    type = ElementAverageValue
    variable = temp
    block = '${fuel_blocks}'
    execute_on = 'initial timestep_end'
  []
  [max_Tfuel]
    type = ElementExtremeValue
    variable = temp
    value_type = max
    block = '${fuel_blocks}'
    execute_on = 'initial timestep_end'
  []
  [avg_Tsleeve]
    type = ElementAverageValue
    variable = temp
    block = '${sleeve_blocks}'
    execute_on = 'initial timestep_end'
  []
  [max_Tsleeve]
    type = ElementExtremeValue
    variable = temp
    value_type = max
    block = '${sleeve_blocks}'
    execute_on = 'initial timestep_end'
  []
  [avg_scaled_pdensity]
    type = ElementAverageValue
    variable = scaled_power_density
    block = '${fuel_blocks}'
    execute_on = 'initial timestep_end'
  []
  [avg_pdensity]
    type = ElementAverageValue
    variable = power_density
    block = '${fuel_blocks}'
    execute_on = 'initial timestep_end'
  []
  [avg_Twall]
    type = SideAverageValue
    variable = temp
    boundary = Relap7_Tw
    execute_on = 'initial timestep_end'
  []
  [avg_Tfluid]
    type = ElementAverageValue
    variable = Tfluid
    execute_on = 'initial timestep_end'
  []
  [avg_Tmod]
    type = ElementAverageValue
    variable = Tmod
    execute_on = 'initial timestep_end'
  []
  [convection_sleeve_fluid]
    type = ConvectiveHeatTransferSideIntegral
    T_fluid_var = Tfluid
    htc_var = Hw_outer
    T_solid = temp
    boundary = 'Inner_outer_ring'
  []
  [convection_mod_fluid]
    type = ConvectiveHeatTransferSideIntegral
    T_fluid_var = Tfluid
    htc_var = Hw_inner
    T_solid = temp
    boundary = 'Relap7_Tw'
  []
  [bdy_heat_flux]
    type = ConvectiveHeatTransferSideIntegral
    boundary = 'Outer_outer_ring'
    T_solid = temp
    T_fluid_var = Tmod
    htc = 1e5
  []
  [bdy_heat_flux_tot]
    type = ParsedPostprocessor
    function = '${npins} * bdy_heat_flux'
    pp_names = 'bdy_heat_flux'
  []
[]

[Outputs]
  csv = true
[]
(htgr/httr/steady_state_and_null_transient/fuel_elem_null.i)

RELAP-7 Control Rod Assembly thermal-hydraulics model

# ==============================================================================
# Control Rod Assembly thermal-hydraulics model (null)
# Application: RELAP-7 (Sabertooth with child applications)
# ------------------------------------------------------------------------------
# Idaho Falls, INL, April 2023
# Author(s): Vincent Laboure
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
# ==============================================================================

# - the number of pins should be 33 for fuel blocks of type 1 and 2, 31 for type 3 and 4.
#   For the control rods, the component below should be adjusted and the number of pins is set to 2
#
# - mdot should be computed from the total mass flow rate Mdot with:
#        Mdot = 12.4 kg/s for the 30 MW with Tout = 850C and for the 9 MW
#        Mdot = 10.2 kg/s for the 30 MW with Tout = 950C
#   Then, an assumption must be made to know how much of the helium cools the CR, which is chosen to be
#   5.5% < (1 - alpha) < 8.8% of Mdot [see Maruyama 1994, Shi 2015]. For now choose alpha = 92%
#   For the fuel regions, the mass flow rate is then:
#       mdot = alpha * Mdot / 954, 954 being the total number of fuel channels (33*12+31*18)
#   For the CR regions, the mass flow rate is then:
#       mdot = (1 - alpha) * Mdot / 32, 32 being the total number of CR cooling channels (2*16)

# For now, this assumes that there is no CR impeding the flow (although most of them have partially inserted CRs)
# Thus, we only have one component

pressure = 2.8e6 # pressure in Pa

Mdot = 12.4 # total mass flow rate of coolant in kg/s
Npins = 32 # total number of CR channels
npins = 2 # number of CR channels in the current assembly
alpha = 0.92 # (1-alpha) is the by-pass flow proportion

mdot = '${fparse Mdot * (1-alpha) / Npins }' # kg/s
mdot_npins_over_Mdot = '${fparse (1-alpha) * npins / Npins }'

D_ext = 1.23e-1 # external diameter in m

p = 0.36 # hexagonal pitch in m
Sblock = '${fparse 0.5 * sqrt(3) * p * p}' # Cross-section area of the hexagonal block in m^2
htc_homo_scaling = '${fparse pi * D_ext * npins / Sblock}'

Tinlet = 453.15 # fluid inlet temperature

[GlobalParams]
  # HTTR system at 2.8 MPa
  # Tinlet 180C for 9MW

  scaling_factor_1phase = '1 1e-2 1e-5' # scaling mass, mommentum, and energy
  closures = closure_high_temp_gas
[]

[Closures]
  [closure_high_temp_gas]
    type = Closures1PhaseHighTempGas
  []
[]

[AuxVariables]
  [T_wall:1]
    # dummy variable to make the transfers not crash
  []
  [T_wall:2]
    # dummy variable to make the transfers not crash
  []

  [Hw_inner] # heat transfer coefficient for the inner radius of fuel cooling channels [W/K/m^2] (not homogenized because only meant to be transferred to BISON)
    # dummy variable to make the transfers not crash
    family = MONOMIAL
    order = CONSTANT
  []
  [Hw_outer] # heat transfer coefficient for the outer radius of cooling channels [W/K/m^2]
    family = MONOMIAL
    order = CONSTANT
  []
  [Hw_outer_homo] # homogenized heat transfer coefficient for the outer radius of cooling channels [W/K/m^3]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [Hw_outer_aux]
    type = ADMaterialRealAux
    variable = Hw_outer
    property = Hw
  []
  [Hw_outer_homo_aux]
    type = ParsedAux
    variable = Hw_outer_homo
    expression = 'Hw_outer * ${htc_homo_scaling}'
    coupled_variables = 'Hw_outer'
  []
[]

[UserObjects]
  [avg_Tfluid_UO]
    type = LayeredAverage
    variable = T
    direction = x
    bounds = '0 0.58 1.16 1.74 2.32 2.90 3.48 4.06 4.64 5.22' # axial boundaries of each block
    execute_on = 'timestep_end'
  []
  [avg_Hw_inner_UO]
    type = LayeredAverage
    variable = Hw_inner
    direction = x
    bounds = '0 0.58 1.16 1.74 2.32 2.90 3.48 4.06 4.64 5.22' # axial boundaries of each block
    execute_on = 'timestep_end'
  []
  [avg_Hw_outer_UO]
    type = LayeredAverage
    variable = Hw_outer
    direction = x
    bounds = '0 0.58 1.16 1.74 2.32 2.90 3.48 4.06 4.64 5.22' # axial boundaries of each block
    execute_on = 'timestep_end'
  []
  [avg_Hw_outer_homo_UO]
    type = LayeredAverage
    variable = Hw_outer_homo
    direction = x
    bounds = '0 0.58 1.16 1.74 2.32 2.90 3.48 4.06 4.64 5.22' # axial boundaries of each block
    execute_on = 'timestep_end'
  []
[]

[FluidProperties]
  [helium]
    type = HeliumSBTLFluidProperties
  []
[]

[Components]
  [CR_cooling_channel]
    type = FlowChannel1Phase
    fp = helium
    # Geometry
    position = '5.22 0 0'
    orientation = '-1.0 0 0'
    A = '${fparse 0.25 * pi * D_ext * D_ext}' # Area of the pipe m^2
    D_h = ${D_ext} # Hydraulic diameter A = Dout
    f = 0 # Wall friction

    length = 5.22 # m
    n_elems = 300
  []

  [inlet]
    type = InletMassFlowRateTemperature1Phase
    input = 'CR_cooling_channel:in'
    m_dot = ${mdot}
    T = ${Tinlet} # K
  []

  [outlet]
    type = Outlet1Phase
    input = 'CR_cooling_channel:out'
    p = ${pressure}
  []

  [ht_ext]
    type = HeatTransferFromExternalAppTemperature1Phase
    flow_channel = CR_cooling_channel
    P_hf = '${fparse pi * D_ext}' # Heat flux perimeter m
  []
[]

[Preconditioning]
  [SMP_PJFNK]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = 'lu       mumps'
  []
[]

[Executioner]
  type = Transient
  scheme = 'bdf2'
  solve_type = 'NEWTON'
  line_search = 'basic'

  dtmin = 1.e-8

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e-2
  []

  nl_rel_tol = 1e-6
  nl_abs_tol = 5e-8
  nl_max_its = 20

  l_tol = 1e-3
  l_max_its = 100

  start_time = 0
  end_time = 50

  [Quadrature]
    type = GAUSS
    order = SECOND
  []
[]

[Postprocessors]
  [avg_T]
    type = ElementAverageValue
    variable = T
    execute_on = 'initial timestep_end'
  []
  [max_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    execute_on = 'initial timestep_end'
  []
  [avg_Twall]
    type = ElementAverageValue
    variable = T_wall
    execute_on = 'initial timestep_end'
  []
  [min_Twall]
    type = ElementExtremeValue
    variable = T_wall
    value_type = min
    execute_on = 'initial timestep_end'
  []
  [avg_Hw_outer]
    type = ElementAverageValue
    variable = Hw_outer
    execute_on = 'initial timestep_end'
  []
  [mdot_in]
    type = PointValue
    variable = rhouA
    point = '5.22 0 0'
    execute_on = 'initial timestep_end'
  []
  [mdot_out]
    type = PointValue
    variable = rhouA
    point = '0 0 0'
    execute_on = 'initial timestep_end'
  []
  [H_in]
    type = PointValue
    variable = H
    point = '5.22 0 0'
    execute_on = 'initial timestep_end'
  []
  [H_out]
    type = PointValue
    variable = H
    point = '0 0 0'
    execute_on = 'initial timestep_end'
  []
  [delta_H]
    type = ParsedPostprocessor
    function = '${npins} * ${mdot} * (H_out - H_in)'
    pp_names = 'H_in H_out'
    execute_on = 'initial timestep_end'
  []
  [Tout]
    type = SideAverageValue
    boundary = 'CR_cooling_channel:out'
    variable = T
    execute_on = 'initial timestep_end'
  []
  [Tout_weighted]
    type = ScalePostprocessor
    value = Tout
    scaling_factor = ${mdot_npins_over_Mdot}
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  [console]
    type = Console
    max_rows = 1
  []
[]
(htgr/httr/steady_state_and_null_transient/thermal_hydraulics_CR_null.i)

RELAP-7 Fuel Assembly thermal-hydraulics model

# ==============================================================================
# Fuel Assembly thermal-hydraulics model (null)
# Application: RELAP-7 (Sabertooth with child applications)
# ------------------------------------------------------------------------------
# Idaho Falls, INL, April 2023
# Author(s): Vincent Laboure
# If using or referring to this model, please cite as explained in
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
# ==============================================================================

# - the number of pins should be 33 for fuel blocks of type 1 and 2, 31 for type 3 and 4.
#   For the control rods, the component below should be adjusted and the number of pins is set to 2
#
# - mdot should be computed from the total mass flow rate Mdot with:
#        Mdot = 12.4 kg/s for the 30 MW with Tout = 850C and for the 9 MW
#        Mdot = 10.2 kg/s for the 30 MW with Tout = 950C
#   Then, an assumption must be made to know how much of the helium cools the CR, which is chosen to be
#   5.5% < (1 - alpha) < 8.8% of Mdot [see Maruyama 1994, Shi 2015]. For now choose alpha = 92%
#   For the fuel regions, the mass flow rate is then:
#       mdot = alpha * Mdot / 954, 954 being the total number of fuel channels (33*12+31*18)
#   For the CR regions, the mass flow rate is then:
#       mdot = (1 - alpha) * Mdot / 32, 32 being the total number of CR cooling channels (2*16)

pressure = 2.8e6 # pressure in Pa

Mdot = 12.4 # total mass flow rate of coolant in kg/s
Npins = 954 # total number of fuel channels
npins = 33 # number of fuel channels in the current assembly
alpha = 0.92 # (1-alpha) is the by-pass flow proportion

mdot = '${fparse Mdot * alpha / Npins }' # kg/s
mdot_npins_over_Mdot = '${fparse alpha * npins / Npins }'

D_ext = 4.1e-2 # external diameter in m
D_int = 3.4e-2 # internal diameter in m

p = 0.36 # hexagonal pitch in m
Sblock = '${fparse 0.5 * sqrt(3) * p * p}' # Cross-section area of the hexagonal block in m^2
htc_homo_scaling = '${fparse pi * D_ext * npins / Sblock}'

Tinlet = 453.15 # fluid inlet temperature

[GlobalParams]
  # HTTR system at 2.8 MPa
  # Tinlet 180C for 9MW
  scaling_factor_1phase = '1 1e-2 1e-5' # scaling mass, mommentum, and energy
  closures = closure_high_temp_gas
[]

[Closures]
  [closure_high_temp_gas]
    type = Closures1PhaseHighTempGas
  []
[]

[AuxVariables]
  [Hw_inner] # heat transfer coefficient for the inner radius of fuel cooling channels [W/K/m^2]
    block = 'pipe1'
    family = MONOMIAL
    order = CONSTANT
  []
  [Hw_outer] # heat transfer coefficient for the outer radius of cooling channels [W/K/m^2]
    family = MONOMIAL
    order = CONSTANT
  []
  [Hw_outer_homo] # homogenized heat transfer coefficient for the outer radius of cooling channels [W/K/m^3]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [Hw_inner_aux]
    type = ADMaterialRealAux
    block = 'pipe1'
    variable = Hw_inner
    property = Hw:1
  []
  [Hw_outer_aux]
    type = ADMaterialRealAux
    block = 'pipe1'
    variable = Hw_outer
    property = Hw:2
  []
  [Hw_outer_aux_topbottom]
    type = ADMaterialRealAux
    block = 'pipe_top pipe_bottom'
    variable = Hw_outer
    property = Hw
  []
  [Hw_outer_homo_aux]
    type = ParsedAux
    variable = Hw_outer_homo
    expression = 'Hw_outer * ${htc_homo_scaling}'
    coupled_variables = 'Hw_outer'
  []
[]

[UserObjects]
  [avg_Tfluid_UO]
    type = LayeredAverage
    variable = T
    direction = x
    bounds = '0 0.58 1.16 1.74 2.32 2.90 3.48 4.06 4.64 5.22' # axial boundaries of each block
    execute_on = 'timestep_end'
  []
  [avg_Hw_inner_UO]
    type = LayeredAverage
    variable = Hw_inner
    direction = x
    block = 'pipe1'
    bounds = '0 0.58 1.16 1.74 2.32 2.90 3.48 4.06 4.64 5.22' # axial boundaries of each block
    execute_on = 'timestep_end'
  []
  [avg_Hw_outer_UO]
    type = LayeredAverage
    variable = Hw_outer
    direction = x
    bounds = '0 0.58 1.16 1.74 2.32 2.90 3.48 4.06 4.64 5.22' # axial boundaries of each block
    execute_on = 'timestep_end'
  []
  [avg_Hw_outer_homo_UO]
    type = LayeredAverage
    variable = Hw_outer_homo
    direction = x
    bounds = '0 0.58 1.16 1.74 2.32 2.90 3.48 4.06 4.64 5.22' # axial boundaries of each block
    execute_on = 'timestep_end'
  []
[]

[FluidProperties]
  [helium]
    type = HeliumSBTLFluidProperties
  []
[]

[Components]
  [pipe_top]
    type = FlowChannel1Phase
    fp = helium
    # Geometry
    position = '5.22 0 0'
    orientation = '-1.0 0 0'
    A = '${fparse 0.25 * pi * D_ext * D_ext}' # Area of the pipe m^2
    D_h = ${D_ext} # Hydraulic diameter A = Dout
    f = 0 # Wall friction

    length = 1.16 # m
    n_elems = 60
  []
  [pipe1]
    type = FlowChannel1Phase
    fp = helium
    # Geometry
    position = '4.06 0 0'
    orientation = '-1.0 0 0'
    A = '${fparse 0.25 * pi * (D_ext * D_ext - D_int * D_int)}' # Area of the pipe m^2
    D_h = '${fparse D_ext - D_int}' # Hydraulic diameter A = Dout-Din
    f = 0 # Wall friction

    length = 2.9 # m
    n_elems = 180
  []
  [pipe_bottom]
    type = FlowChannel1Phase
    fp = helium
    # Geometry
    position = '1.16 0 0'
    orientation = '-1.0 0 0'
    A = '${fparse 0.25 * pi * D_ext * D_ext}' # Area of the pipe m^2
    D_h = ${D_ext} # Hydraulic diameter A = Dout
    f = 0 # Wall friction

    length = 1.16
    n_elems = 60
  []

  [inlet]
    type = InletMassFlowRateTemperature1Phase
    input = 'pipe_top:in'
    m_dot = ${mdot}
    T = ${Tinlet} # K
  []

  [junction1]
    type = JunctionOneToOne1Phase
    connections = 'pipe_top:out pipe1:in'
  []
  [junction2]
    type = JunctionOneToOne1Phase
    connections = 'pipe1:out pipe_bottom:in'
  []

  [outlet]
    type = Outlet1Phase
    input = 'pipe_bottom:out'
    p = ${pressure}
  []

  [ht_int]
    type = HeatTransferFromExternalAppTemperature1Phase
    flow_channel = pipe1
    P_hf = '${fparse pi * D_int}' # Heat flux perimeter m
  []

  [ht_ext_top]
    type = HeatTransferFromExternalAppTemperature1Phase
    flow_channel = pipe_top
    P_hf = '${fparse pi * D_ext}' # Heat flux perimeter m
  []
  [ht_ext]
    type = HeatTransferFromExternalAppTemperature1Phase
    flow_channel = pipe1
    P_hf = '${fparse pi * D_ext}' # Heat flux perimeter m
  []
  [ht_ext_bottom]
    type = HeatTransferFromExternalAppTemperature1Phase
    flow_channel = pipe_bottom
    P_hf = '${fparse pi * D_ext}' # Heat flux perimeter m
  []
[]

[Preconditioning]
  [SMP_PJFNK]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = 'lu       mumps'
  []
[]

[Executioner]
  type = Transient
  scheme = 'bdf2'
  solve_type = 'NEWTON'
  line_search = 'basic'

  dtmin = 1.e-8

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e-2
  []

  nl_rel_tol = 1e-6
  nl_abs_tol = 5e-8 # original: 1e-8 (but can result in solve failure in rare cases)
  nl_max_its = 20 # original: 5 (but can result in solve failure in rare cases)

  l_tol = 1e-3
  l_max_its = 100

  start_time = 0
  end_time = 50

  [Quadrature]
    type = GAUSS
    order = SECOND
  []
[]

[Postprocessors]
  [avg_T]
    type = ElementAverageValue
    variable = T
    execute_on = 'initial timestep_end'
  []
  [max_T]
    type = ElementExtremeValue
    variable = T
    value_type = max
    execute_on = 'initial timestep_end'
  []
  [avg_Twall]
    type = ElementAverageValue
    variable = T_wall
    block = 'pipe_top pipe_bottom'
    execute_on = 'initial timestep_end'
  []
  [avg_Twall1]
    type = ElementAverageValue
    variable = T_wall:1
    block = 'pipe1'
    execute_on = 'initial timestep_end'
  []
  [max_Twall1]
    type = ElementExtremeValue
    variable = T_wall:1
    value_type = max
    block = 'pipe1'
    execute_on = 'initial timestep_end'
  []
  [min_Twall1]
    type = ElementExtremeValue
    variable = T_wall:1
    value_type = min
    block = 'pipe1'
    execute_on = 'initial timestep_end'
  []
  [avg_Twall2]
    type = ElementAverageValue
    variable = T_wall:2
    block = 'pipe1'
    execute_on = 'initial timestep_end'
  []
  [max_Twall2]
    type = ElementExtremeValue
    variable = T_wall:2
    value_type = max
    block = 'pipe1'
    execute_on = 'initial timestep_end'
  []
  [min_Twall2]
    type = ElementExtremeValue
    variable = T_wall:2
    value_type = min
    block = 'pipe1'
    execute_on = 'initial timestep_end'
  []
  [avg_Hw_inner]
    type = ElementAverageValue
    variable = Hw_inner
    block = 'pipe1'
    execute_on = 'initial timestep_end'
  []
  [avg_Hw_outer]
    type = ElementAverageValue
    variable = Hw_outer
    execute_on = 'initial timestep_end'
  []
  [mdot_in]
    type = PointValue
    variable = rhouA
    point = '5.22 0 0'
    execute_on = 'initial timestep_end'
  []
  [mdot_out]
    type = PointValue
    variable = rhouA
    point = '0 0 0'
    execute_on = 'initial timestep_end'
  []
  [H_in]
    type = PointValue
    variable = H
    point = '5.22 0 0'
    execute_on = 'initial timestep_end'
  []
  [H_out]
    type = PointValue
    variable = H
    point = '0 0 0'
    execute_on = 'initial timestep_end'
  []
  [delta_H]
    type = ParsedPostprocessor
    function = '${npins} * ${mdot} * (H_out - H_in)'
    pp_names = 'H_in H_out'
    execute_on = 'initial timestep_end'
  []
  [Tout]
    type = SideAverageValue
    boundary = 'pipe_bottom:out'
    variable = T
    execute_on = 'initial timestep_end'
  []
  [Tout_weighted]
    type = ScalePostprocessor
    value = Tout
    scaling_factor = ${mdot_npins_over_Mdot}
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  [console]
    type = Console
    max_rows = 1
  []
[]
(htgr/httr/steady_state_and_null_transient/thermal_hydraulics_fuel_pins_null.i)

Running the simulation

To apply for access to Sabertooth and access to INL High Performance Computing, please visit NCRC. These inputs can be run on your local machine or on a computing cluster, depending on your level of access to the simulation software and the computing power available to you.

Idaho National Laboratory HPC

In all cases (steady-state and null-transient - which should be executed in this order), run the Griffin input with Sabertooth. With at least level 1 access to Sabertooth, the sabertooth module can be loaded and from the sabertooth-opt executable can be used. If you have already run the steady state calculation, you can skip the second step below.



module load use.moose moose-apps sabertooth


mpirun -n 6 sabertooth-opt -i neutronics_eigenvalue.i

mpirun -n 6 sabertooth-opt -i neutronics_null.i

Local Device

In all cases (steady-state and null-transient - which should be executed in this order), run the Griffin input with Sabertooth.

For instance, if you have level 2 access to the binaries, you can download sabertooth through miniconda and use them as follows:



conda activate sabertooth

mpirun -n 6 sabertooth-opt -i neutronics_eigenvalue.i

mpirun -n 6 sabertooth-opt -i neutronics_null.i

Citing and credits

If using or modifying this model, please cite:


Vincent Labouré, Javier Ortensi, Nicolas Martin, Paolo Balestra, Derek Gaston, Yinbin Miao, Gerhard Strydom, "Improved multiphysics model of the High Temperature Engineering Test Reactor for the simulation of loss-of-forced-cooling experiments", Annals of Nuclear Energy, Volume 189, 2023, 109838, ISSN 0306-4549, https://doi.org/10.1016/j.anucene.2023.109838. (https://www.sciencedirect.com/science/article/pii/S0306454923001573)

Vincent M Laboure, Matilda A Lindell, Javier Ortensi, Gerhard Strydom and Paolo Balestra, "FY22 Status Report on the ART-GCR CMVB and CNWG International Collaborations." INL/RPT-22-68891-Rev000 Web. doi:10.2172/1893099.

Documentation written by Kylee Swanson.