Multiphysics for a TRISO Gas-Cooled Compact

OpenMC Input Files

The neutronics physics is solved over the entire domain with OpenMC. The OpenMC wrapping used for the OpenMC-NekRS-MOOSE coupling is described in the openmc_nek.i input file. We begin by defining a number of constants and by setting up the mesh mirror on which OpenMC will receive temperature and density from THM-MOOSE, and on which OpenMC will write the fission heat source. Because the coupled applications use length units of meters, the mesh mirror must also be in units of meters. For simplicity, we just use the same mesh as generated with solid_mesh.i earlier, though this it not necessary.

!include common_input.i

# This input file runs coupled OpenMC Monte Carlo transport, MOOSE heat
# conduction, and NekRS fluid flow and heat transfer.
# This input should be run with:
#
# cardinal-opt -i openmc_nek.i

num_layers_for_THM = 150
density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3
N = 1000

[Mesh]
  [solid]
    type = FileMeshGenerator
    file = solid_mesh_in.e
  []
[]

Next, we define a number of auxiliary variables to query the OpenMC solution and set up multiphysics coupling. We add cell_temperature and cell_density in order to read the cell temperatures and densities directly from OpenMC in order to visualize the temperatures and densities ultimately applied to OpenMC's cells. In order to compute density using the ideal gas EOS given a temperature and a fixed pressure, we use a FluidDensityAux to set density.

[AuxVariables<<<{"href": "../syntax/AuxVariables/index.html"}>>>]
  [cell_temperature]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [cell_density]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
[]

[AuxKernels<<<{"href": "../syntax/AuxKernels/index.html"}>>>]
  [cell_temperature]
    type = CellTemperatureAux<<<{"description": "OpenMC cell temperature (K), mapped to each MOOSE element", "href": "../source/auxkernels/CellTemperatureAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = cell_temperature
  []
  [cell_density]
    type = CellDensityAux<<<{"description": "OpenMC fluid density (kg/m$^3$), mapped to each MOOSE element", "href": "../source/auxkernels/CellDensityAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = cell_density
  []
  [density]
    type = FluidDensityAux<<<{"description": "Computes density from pressure and temperature", "href": "../source/auxkernels/FluidDensityAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = density
    p<<<{"description": "Pressure"}>>> = ${outlet_P}
    T<<<{"description": "Temperature"}>>> = nek_temp
    fp<<<{"description": "The name of the user object for fluid properties"}>>> = helium
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_begin'
  []
[]

[FluidProperties<<<{"href": "../syntax/FluidProperties/index.html"}>>>]
  [helium]
    type = IdealGasFluidProperties<<<{"description": "Fluid properties for an ideal gas", "href": "../source/fluidproperties/IdealGasFluidProperties.html"}>>>
    molar_mass<<<{"description": "Constant molar mass of the fluid (kg/mol)"}>>> = 4e-3
    gamma<<<{"description": "gamma value (cp/cv)"}>>> = 1.668282 # should correspond to  Cp = 5189 J/kg/K
    k<<<{"description": "Thermal conductivity, W/(m-K)"}>>> = 0.2556
    mu<<<{"description": "Dynamic viscosity, Pa.s"}>>> = 3.22639e-5
  []
[]

Next, we define the necessary initial conditions using functions; all temperatures (nek_temp and solid_temp are to be discussed shortly) are set to a linear variation from the inlet to outlet fluid temperatures.

[ICs<<<{"href": "../syntax/Cardinal/ICs/index.html"}>>>]
  [temp]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = nek_temp
    function<<<{"description": "The initial condition function."}>>> = temp_ic
  []
  [solid_temp]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = solid_temp
    function<<<{"description": "The initial condition function."}>>> = temp_ic
  []
[]

[Functions<<<{"href": "../syntax/Functions/index.html"}>>>]
  [temp_ic]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / (${mdot} / ${n_bundles} / ${n_coolant_channels_per_block}) / ${fluid_Cp}'
  []
[]

The wrapping of OpenMC is specified in the [Problem] block and the addition of tallies is done in the [Tallies] block. Here, we indicate that we will provide both temperature and density feedback to OpenMC. In order to visualize the tally standard deviation, we output the fission tally standard deviation using the output parameter. The heat source from OpenMC will be relaxed using Robbins-Monro relaxation.

By default, OpenMC will try to read temperature from the temp variable. However, in this case we have multiple applications (NekRS/THM for the fluids, MOOSE for the solids) which want to send temperatures into OpenMC. To be sure one data transfer does not overwrite the second, we need to tell OpenMC the names of the variables to read temperature from. Cardinal contains convenient syntax to automatically set up the necessary receiver variables for temperature, by using the temperature_variables and temperature_blocks parameters.

Finally, a number of "triggers" are used to automatically terminate OpenMC's active batches once reaching desired uncertainties in and the fission tally. The number of batches here is terminated once both of the following are satisfied:

• Standard deviation in is less than 75 pcm

• Maximum fission tally relative error is less than 1%

These criteria are checked every batch_interval, up to a maximum number of batches.

[Problem<<<{"href": "../syntax/Problem/index.html"}>>>]
  type = OpenMCCellAverageProblem

  power = ${unit_cell_power}
  scaling = 100.0
  density_blocks = ${density_blocks}
  cell_level = 1

  relaxation = robbins_monro

  temperature_variables = 'solid_temp;   nek_temp'
  temperature_blocks = '${temperature_blocks}; ${density_blocks}'

  k_trigger = std_dev
  k_trigger_threshold = 7.5e-4
  batches = 40
  max_batches = 100
  batch_interval = 5

  [Tallies<<<{"href": "../syntax/Problem/Tallies/index.html"}>>>]
    [heat_source]
      type = CellTally<<<{"description": "A class which implements distributed cell tallies.", "href": "../source/tallies/CellTally.html"}>>>
      block<<<{"description": "Subdomains for which to add tallies in OpenMC. If not provided, tallies will be applied over the entire domain corresponding to the [Mesh] block."}>>> = ${fuel_blocks}
      name<<<{"description": "Auxiliary variable name(s) to use for OpenMC tallies. If not specified, defaults to the names of the scores"}>>> = heat_source

      check_equal_mapped_tally_volumes<<<{"description": "Whether to check if the tallied cells map to regions in the mesh of equal volume. This can be helpful to ensure that the volume normalization of OpenMC's tallies doesn't introduce any unintentional distortion just because the mapped volumes are different. You should only set this to true if your OpenMC tally cells are all the same volume!"}>>> = true

      trigger<<<{"description": "Trigger criterion to determine when OpenMC simulation is complete based on tallies. If multiple scores are specified in 'score, this same trigger is applied to all scores."}>>> = rel_err
      trigger_threshold<<<{"description": "Threshold for the tally trigger"}>>> = 1e-2

      output<<<{"description": "UNRELAXED field(s) to output from OpenMC for each tally score. unrelaxed_tally_std_dev will write the standard deviation of each tally into auxiliary variables named *_std_dev. unrelaxed_tally_rel_error will write the relative standard deviation (unrelaxed_tally_std_dev / unrelaxed_tally) of each tally into auxiliary variables named *_rel_error. unrelaxed_tally will write the raw unrelaxed tally into auxiliary variables named *_raw (replace * with 'name')."}>>> = 'unrelaxed_tally_std_dev'
    []
  []
[]

Next, we define a transient executioner - while OpenMC is technically solving a steady -eigenvalue problem, using a time-dependent executioner with the notion of a "time step" will allow us to control the frequency with which OpenMC sends data to/from it's sub-app (MOOSE heat conduction). Here, we set the time step to be 1000 times the NekRS fluid time step. We will terminate the solution once reaching the specified steady_state_tolerance, which by setting check_aux = true will terminate the entire solution once reaching less than 1% change in all auxiliary variables in the OpenMC input file (which includes temperatures and power).

[Executioner<<<{"href": "../syntax/Executioner/index.html"}>>>]
  type = Transient
  dt = '${fparse N * nek_dt * t0}'

  # This is somewhat loose, and you should use a tighter tolerance for production runs;
  # we use 1e-2 to obtain a faster tutorial
  steady_state_detection = true
  check_aux = true
  steady_state_tolerance = 1e-2
[]

Next, we define the MOOSE heat conduction sub-application and data transfers to/from that application. Most important to note is that while the MOOSE heat transfer module does not itself compute fluid temperature, we see a transfer getting the fluid temperature that has been transferred to the MOOSE heat transfer module by the doubly-nested NekRS sub-application.

[MultiApps<<<{"href": "../syntax/MultiApps/index.html"}>>>]
  [bison]
    type = TransientMultiApp<<<{"description": "MultiApp for performing coupled simulations with the parent and sub-application both progressing in time.", "href": "../source/multiapps/TransientMultiApp.html"}>>>
    input_files<<<{"description": "The input file for each App.  If this parameter only contains one input file it will be used for all of the Apps.  When using 'positions_from_file' it is also admissable to provide one input_file per file."}>>> = 'solid_nek.i'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
    sub_cycling<<<{"description": "Set to true to allow this MultiApp to take smaller timesteps than the rest of the simulation.  More than one timestep will be performed for each parent application timestep"}>>> = true
  []
[]

[Transfers<<<{"href": "../syntax/Transfers/index.html"}>>>]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer<<<{"description": "Transfers the value to the target domain from a combination/interpolation of the values on the nearest nodes in the source domain, using coefficients based on the distance to each node.", "href": "../source/transfers/MultiAppGeometricInterpolationTransfer.html"}>>>
    source_variable<<<{"description": "The variable to transfer from."}>>> = T
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = solid_temp
    from_multi_app<<<{"description": "The name of the MultiApp to receive data from"}>>> = bison
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer<<<{"description": "Transfers field data at the MultiApp position using the finite element shape functions from the origin application.", "href": "../source/transfers/MultiAppGeneralFieldShapeEvaluationTransfer.html"}>>>
    source_variable<<<{"description": "The variable to transfer from."}>>> = heat_source
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = power
    to_multi_app<<<{"description": "The name of the MultiApp to transfer the data to"}>>> = bison
    from_postprocessors_to_be_preserved<<<{"description": "The name of the Postprocessor in the from-app to evaluate an adjusting factor."}>>> = heat_source
    to_postprocessors_to_be_preserved<<<{"description": "The name of the Postprocessor in the to-app to evaluate an adjusting factor."}>>> = power
  []
  [temp_from_nek]
    type = MultiAppGeneralFieldShapeEvaluationTransfer<<<{"description": "Transfers field data at the MultiApp position using the finite element shape functions from the origin application.", "href": "../source/transfers/MultiAppGeneralFieldShapeEvaluationTransfer.html"}>>>
    source_variable<<<{"description": "The variable to transfer from."}>>> = nek_bulk_temp
    from_multi_app<<<{"description": "The name of the MultiApp to receive data from"}>>> = bison
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = nek_temp
  []
[]

Next, we define several postprocessors for querying the solution. The heat_source postprocessor will be used to ensure conservation of power when sent to the MOOSE heat conduction application. All other postprocessors are used for general solution monitoring purposes.

[Postprocessors<<<{"href": "../syntax/Postprocessors/index.html"}>>>]
  [heat_source]
    type = ElementIntegralVariablePostprocessor<<<{"description": "Computes a volume integral of the specified variable", "href": "../source/postprocessors/ElementIntegralVariablePostprocessor.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = heat_source
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'transfer initial timestep_end'
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = ${fuel_blocks}
  []
  [max_tally_err]
    type = TallyRelativeError<<<{"description": "Maximum/minimum tally relative error", "href": "../source/postprocessors/TallyRelativeError.html"}>>>
    value_type<<<{"description": "Whether to give the maximum or minimum tally relative error"}>>> = max
  []
  [k]
    type = KEigenvalue<<<{"description": "k eigenvalue computed by OpenMC", "href": "../source/postprocessors/KEigenvalue.html"}>>>
  []
  [k_std_dev]
    type = KEigenvalue<<<{"description": "k eigenvalue computed by OpenMC", "href": "../source/postprocessors/KEigenvalue.html"}>>>
    output<<<{"description": "The value to output. Options are $k_{eff}$ (mean), the standard deviation of $k_{eff}$ (std_dev), or the relative error of $k_{eff}$ (rel_err)."}>>> = 'std_dev'
  []
  [min_power]
    type = ElementExtremeValue<<<{"description": "Finds either the min or max elemental value of a variable over the domain.", "href": "../source/postprocessors/ElementExtremeValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on"}>>> = heat_source
    value_type<<<{"description": "Type of extreme value to return. 'max' returns the maximum value. 'min' returns the minimum value. 'max_abs' returns the maximum of the absolute value."}>>> = min
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = ${fuel_blocks}
  []
  [max_power]
    type = ElementExtremeValue<<<{"description": "Finds either the min or max elemental value of a variable over the domain.", "href": "../source/postprocessors/ElementExtremeValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on"}>>> = heat_source
    value_type<<<{"description": "Type of extreme value to return. 'max' returns the maximum value. 'min' returns the minimum value. 'max_abs' returns the maximum of the absolute value."}>>> = max
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = ${fuel_blocks}
  []
[]

For postprocessing, we also compute the average power distribution in a number of layers using a LayeredAverage and output to CSV using a SpatialUserObjectVectorPostprocessor, in combination with a CSV output.

[UserObjects<<<{"href": "../syntax/UserObjects/index.html"}>>>]
  [average_power_axial]
    type = LayeredAverage<<<{"description": "Computes averages of variables over layers", "href": "../source/userobjects/LayeredAverage.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = heat_source
    direction<<<{"description": "The direction of the layers."}>>> = z
    num_layers<<<{"description": "The number of layers."}>>> = ${num_layers_for_plots}
    block<<<{"description": "The list of block ids (SubdomainID) that this object will be applied"}>>> = ${fuel_blocks}
  []
[]

[VectorPostprocessors<<<{"href": "../syntax/VectorPostprocessors/index.html"}>>>]
  [power_avg]
    type = SpatialUserObjectVectorPostprocessor<<<{"description": "Outputs the values of a spatial user object in the order of the specified spatial points", "href": "../source/vectorpostprocessors/SpatialUserObjectVectorPostprocessor.html"}>>>
    userobject<<<{"description": "The userobject whose values are to be reported"}>>> = average_power_axial
  []
[]

[Outputs<<<{"href": "../syntax/Outputs/index.html"}>>>]
  [out]
    type = Exodus<<<{"description": "Object for output data in the Exodus format", "href": "../source/outputs/Exodus.html"}>>>
    hide<<<{"description": "A list of the variables and postprocessors that should NOT be output to the Exodus file (may include Variables, ScalarVariables, and Postprocessor names)."}>>> = 'solid_temp nek_temp'
  []

  [csv]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../source/outputs/CSV.html"}>>>
    file_base<<<{"description": "The desired solution output name without an extension. If not provided, MOOSE sets it with Outputs/file_base when available. Otherwise, MOOSE uses input file name and this object name for a master input or uses master file_base, the subapp name and this object name for a subapp input to set it."}>>> = 'csv_nek/openmc_nek'
  []
[]

Solid Input Files

The solid heat conduction physics is solved over the solid regions of the unit cell using the MOOSE heat transfer module. The input file for this portion of the physics is the solid_nek.i input. We begin by defining a number of constants and by setting up the mesh for solving heat conduction.

channel_diameter = 0.016                 # diameter of the coolant channels (m)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected
kernel_radius = 214.85e-6                # fissile kernel outer radius (m)
buffer_radius = 314.85e-6                # buffer outer radius (m)
iPyC_radius = 354.85e-6                  # inner PyC outer radius (m)
SiC_radius = 389.85e-6                   # SiC outer radius (m)
oPyC_radius = 429.85e-6                  # outer PyC outer radius (m)
buffer_k = 0.5                           # buffer thermal conductivity (W/m/K)
PyC_k = 4.0                              # PyC thermal conductivity (W/m/K)
SiC_k = 13.9                             # SiC thermal conductivity (W/m/K)
kernel_k = 3.5                           # fissil kernel thermal conductivity (W/m/K)
matrix_k = 15.0                          # graphite matrix thermal conductivity (W/m/K)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
height = 6.343                           # height of the full core (m)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
triso_pf = 0.15                          # TRISO packing fraction (%)

density_blocks = 'coolant'
temperature_blocks = 'graphite compacts compacts_trimmer_tri'
fuel_blocks = 'compacts compacts_trimmer_tri'

unit_cell_power = ${fparse power / (n_bundles * n_coolant_channels_per_block) * unit_cell_height / height}
unit_cell_mdot = ${fparse mdot / (n_bundles * n_coolant_channels_per_block)}

# compute the volume fraction of each TRISO layer in a TRISO particle
# for use in computing average thermophysical properties
kernel_fraction = ${fparse kernel_radius^3 / oPyC_radius^3}
buffer_fraction = ${fparse (buffer_radius^3 - kernel_radius^3) / oPyC_radius^3}
ipyc_fraction = ${fparse (iPyC_radius^3 - buffer_radius^3) / oPyC_radius^3}
sic_fraction = ${fparse (SiC_radius^3 - iPyC_radius^3) / oPyC_radius^3}
opyc_fraction = ${fparse (oPyC_radius^3 - SiC_radius^3) / oPyC_radius^3}

N = 50

U_ref = ${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}
t0 = ${fparse channel_diameter / U_ref}
nek_dt = 6e-3

[Mesh<<<{"href": "../syntax/Mesh/index.html"}>>>]
  [solid]
    type = FileMeshGenerator<<<{"description": "Read a mesh from a file.", "href": "../source/meshgenerators/FileMeshGenerator.html"}>>>
    file<<<{"description": "The filename to read."}>>> = solid_mesh_in.e
  []
[]

Because we have a block in the problem that we don't need to define any material properties on, we technically need to turn off a material coverage check, or else we're going to get an error from MOOSE. FEProblem is just the default problem, which we need to list in order to turn off the material coverage check.

[Problem<<<{"href": "../syntax/Problem/index.html"}>>>]
  type = FEProblem
  material_coverage_check = false
[]

Next, we define the nonlinear variable that this application will solve for (T), the solid temperature. On the solid blocks, we solve the heat equation, but on the fluid blocks that exist exclusively for transferring data, we add a NullKernel to skip the solve in those regions. On the channel wall, the temperature from NekRS is applied as a Dirichlet condition.

[Variables<<<{"href": "../syntax/Variables/index.html"}>>>]
  [T]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = ${inlet_T}
  []
[]

[Kernels<<<{"href": "../syntax/Kernels/index.html"}>>>]
  [diffusion]
    type = HeatConduction<<<{"description": "Diffusive heat conduction term $-\\nabla\\cdot(k\\nabla T)$ of the thermal energy conservation equation", "href": "../source/kernels/HeatConduction.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = ${temperature_blocks}
  []
  [source]
    type = CoupledForce<<<{"description": "Implements a source term proportional to the value of a coupled variable. Weak form: $(\\psi_i, -\\sigma v)$.", "href": "../source/kernels/CoupledForce.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
    v<<<{"description": "The coupled variable which provides the force"}>>> = power
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = ${fuel_blocks}
  []
  [null] # The value of T on the fluid blocks is meaningless
    type = NullKernel<<<{"description": "Kernel that sets a zero residual.", "href": "../source/kernels/NullKernel.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = ${density_blocks}
  []
[]

[BCs<<<{"href": "../syntax/BCs/index.html"}>>>]
  [pin_outer]
    type = MatchedValueBC<<<{"description": "Implements a NodalBC which equates two different Variables' values on a specified boundary.", "href": "../source/bcs/MatchedValueBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
    v<<<{"description": "The variable whose value we are to match."}>>> = nek_temp
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'fluid_solid_interface'
  []
[]

[ICs<<<{"href": "../syntax/Cardinal/ICs/index.html"}>>>]
  [nek_temp]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = nek_temp
    function<<<{"description": "The initial condition function."}>>> = axial_fluid_temp
  []
[]

Next, we define a number of functions to set the solid material properties and define an initial condition for the wall temperature. The solid material properties are then applied with a HeatConductionMaterial.

[ICs<<<{"href": "../syntax/Cardinal/ICs/index.html"}>>>]
  [nek_temp]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = nek_temp
    function<<<{"description": "The initial condition function."}>>> = axial_fluid_temp
  []
[]

[Functions<<<{"href": "../syntax/Functions/index.html"}>>>]
  [k_graphite]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'k_TRISO k_graphite'
    symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = 'k_TRISO k_graphite'
  []
  [axial_fluid_temp]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / ${unit_cell_mdot} / ${fluid_Cp}'
  []
[]

[Materials<<<{"href": "../syntax/Materials/index.html"}>>>]
  [graphite]
    type = HeatConductionMaterial<<<{"description": "General-purpose material model for heat conduction", "href": "../source/materials/HeatConductionMaterial.html"}>>>
    thermal_conductivity_temperature_function<<<{"description": "Thermal conductivity as a function of temperature."}>>> = k_graphite
    temp<<<{"description": "Coupled Temperature"}>>> = T
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 'graphite'
  []
  [compacts]
    type = HeatConductionMaterial<<<{"description": "General-purpose material model for heat conduction", "href": "../source/materials/HeatConductionMaterial.html"}>>>
    thermal_conductivity_temperature_function<<<{"description": "Thermal conductivity as a function of temperature."}>>> = k_compacts
    temp<<<{"description": "Coupled Temperature"}>>> = T
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = ${fuel_blocks}
  []
[]

[AuxVariables<<<{"href": "../syntax/AuxVariables/index.html"}>>>]
  [nek_temp]
    block = ${density_blocks}
  []
  [nek_bulk_temp]
    block = ${density_blocks}
  []
  [flux]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [power]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
[]

[AuxKernels<<<{"href": "../syntax/AuxKernels/index.html"}>>>]
  [flux]
    type = DiffusionFluxAux<<<{"description": "Compute components of flux vector for diffusion problems $(\\vec{J} = -D \\nabla C)$.", "href": "../source/auxkernels/DiffusionFluxAux.html"}>>>
    diffusion_variable<<<{"description": "The name of the variable"}>>> = T
    component<<<{"description": "The desired component of flux."}>>> = normal
    diffusivity<<<{"description": "The name of the diffusivity material property that will be used in the flux computation."}>>> = thermal_conductivity
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = flux
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 'fluid_solid_interface'
  []
[]

Next, we define a number of auxiliary variables - power will receive the heat source from OpenMC, nek_temp will receive the wall temperature from NekRS, nek_bulk_temp will receive the volumetric fluid temperature from NekRS, and finally flux will be used to compute the wall heat flux to send to NekRS.

[AuxVariables<<<{"href": "../syntax/AuxVariables/index.html"}>>>]
  [nek_temp]
    block = ${density_blocks}
  []
  [nek_bulk_temp]
    block = ${density_blocks}
  []
  [flux]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [power]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
[]

[AuxKernels<<<{"href": "../syntax/AuxKernels/index.html"}>>>]
  [flux]
    type = DiffusionFluxAux<<<{"description": "Compute components of flux vector for diffusion problems $(\\vec{J} = -D \\nabla C)$.", "href": "../source/auxkernels/DiffusionFluxAux.html"}>>>
    diffusion_variable<<<{"description": "The name of the variable"}>>> = T
    component<<<{"description": "The desired component of flux."}>>> = normal
    diffusivity<<<{"description": "The name of the diffusivity material property that will be used in the flux computation."}>>> = thermal_conductivity
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = flux
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 'fluid_solid_interface'
  []
[]

We then add a NekRS sub-application and define the transfers to/from NekRS, including the transfer of fluid temperature from the NekRS volume to the dummy coolant blocks in the MOOSE solid model. Because the NekRS mesh mirror will be a volume mirror (in order to extract volumetric temperatures for the neutronics feedback), a significant cost savings can be obtained by using the "minimal transfer" feature of NekRSProblem (which requires sending a dummy Receiver postprocessor, here named synchronization_to_nek, to indicate when data is to be exchanged). For more information, please consult the documentation for NekRSProblem.

[MultiApps<<<{"href": "../syntax/MultiApps/index.html"}>>>]
  [nek]
    type = TransientMultiApp<<<{"description": "MultiApp for performing coupled simulations with the parent and sub-application both progressing in time.", "href": "../source/multiapps/TransientMultiApp.html"}>>>
    input_files<<<{"description": "The input file for each App.  If this parameter only contains one input file it will be used for all of the Apps.  When using 'positions_from_file' it is also admissable to provide one input_file per file."}>>> = 'nek.i'
    sub_cycling<<<{"description": "Set to true to allow this MultiApp to take smaller timesteps than the rest of the simulation.  More than one timestep will be performed for each parent application timestep"}>>> = true
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
[]

[Transfers<<<{"href": "../syntax/Transfers/index.html"}>>>]
  [heat_flux_to_nek]
    type = MultiAppGeneralFieldNearestLocationTransfer<<<{"description": "Transfers field data at the MultiApp position by finding the value at the nearest neighbor(s) in the origin application.", "href": "../source/transfers/MultiAppGeneralFieldNearestLocationTransfer.html"}>>>
    source_variable<<<{"description": "The variable to transfer from."}>>> = flux
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = flux
    from_boundaries<<<{"description": "The boundary we are transferring from (if not specified, whole domain is used)."}>>> = 'fluid_solid_interface'
    to_boundaries<<<{"description": "The boundary we are transferring to (if not specified, whole domain is used)."}>>> = '3'
    to_multi_app<<<{"description": "The name of the MultiApp to transfer the data to"}>>> = nek
  []
  [flux_integral_to_nek]
    type = MultiAppPostprocessorTransfer<<<{"description": "Transfers postprocessor data between the master application and sub-application(s).", "href": "../source/transfers/MultiAppPostprocessorTransfer.html"}>>>
    to_postprocessor<<<{"description": "The name of the Postprocessor in the MultiApp to transfer the value to.  This should most likely be a Reporter Postprocessor."}>>> = flux_integral
    from_postprocessor<<<{"description": "The name of the Postprocessor in the Master to transfer the value from."}>>> = flux_integral
    to_multi_app<<<{"description": "The name of the MultiApp to transfer the data to"}>>> = nek
  []
  [temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer<<<{"description": "Transfers field data at the MultiApp position by finding the value at the nearest neighbor(s) in the origin application.", "href": "../source/transfers/MultiAppGeneralFieldNearestLocationTransfer.html"}>>>
    source_variable<<<{"description": "The variable to transfer from."}>>> = temperature
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = nek_temp
    from_boundaries<<<{"description": "The boundary we are transferring from (if not specified, whole domain is used)."}>>> = '3'
    to_boundaries<<<{"description": "The boundary we are transferring to (if not specified, whole domain is used)."}>>> = 'fluid_solid_interface'
    from_multi_app<<<{"description": "The name of the MultiApp to receive data from"}>>> = nek
  []
  [bulk_temperature_to_bison]
    type = MultiAppGeneralFieldNearestLocationTransfer<<<{"description": "Transfers field data at the MultiApp position by finding the value at the nearest neighbor(s) in the origin application.", "href": "../source/transfers/MultiAppGeneralFieldNearestLocationTransfer.html"}>>>
    source_variable<<<{"description": "The variable to transfer from."}>>> = temperature
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = nek_bulk_temp
    from_multi_app<<<{"description": "The name of the MultiApp to receive data from"}>>> = nek
  []
  [synchronization_to_nek]
    type = MultiAppPostprocessorTransfer<<<{"description": "Transfers postprocessor data between the master application and sub-application(s).", "href": "../source/transfers/MultiAppPostprocessorTransfer.html"}>>>
    to_postprocessor<<<{"description": "The name of the Postprocessor in the MultiApp to transfer the value to.  This should most likely be a Reporter Postprocessor."}>>> = transfer_in
    from_postprocessor<<<{"description": "The name of the Postprocessor in the Master to transfer the value from."}>>> = synchronization_to_nek
    to_multi_app<<<{"description": "The name of the MultiApp to transfer the data to"}>>> = nek
  []
[]

We add several postprocessors to facilitate the data transfers as well as to query the solution.

[Postprocessors<<<{"href": "../syntax/Postprocessors/index.html"}>>>]
  [flux_integral]
    # evaluate the total heat flux for normalization
    type = SideDiffusiveFluxIntegral<<<{"description": "Computes the integral of the diffusive flux over the specified boundary", "href": "../source/postprocessors/SideDiffusiveFluxIntegral.html"}>>>
    diffusivity<<<{"description": "The name of the diffusivity material property that will be used in the flux computation. This must be provided if the variable is of finite element type"}>>> = thermal_conductivity
    variable<<<{"description": "The name of the variable which this postprocessor integrates"}>>> = T
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'fluid_solid_interface'
  []
  [max_fuel_T]
    type = ElementExtremeValue<<<{"description": "Finds either the min or max elemental value of a variable over the domain.", "href": "../source/postprocessors/ElementExtremeValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on"}>>> = T
    value_type<<<{"description": "Type of extreme value to return. 'max' returns the maximum value. 'min' returns the minimum value. 'max_abs' returns the maximum of the absolute value."}>>> = max
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = ${fuel_blocks}
  []
  [max_block_T]
    type = ElementExtremeValue<<<{"description": "Finds either the min or max elemental value of a variable over the domain.", "href": "../source/postprocessors/ElementExtremeValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on"}>>> = T
    value_type<<<{"description": "Type of extreme value to return. 'max' returns the maximum value. 'min' returns the minimum value. 'max_abs' returns the maximum of the absolute value."}>>> = max
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 'graphite'
  []
  [synchronization_to_nek]
    type = Receiver<<<{"description": "Reports the value stored in this processor, which is usually filled in by another object. The Receiver does not compute its own value.", "href": "../source/postprocessors/Receiver.html"}>>>
    default<<<{"description": "The default value"}>>> = 1.0
  []
  [power]
    type = ElementIntegralVariablePostprocessor<<<{"description": "Computes a volume integral of the specified variable", "href": "../source/postprocessors/ElementIntegralVariablePostprocessor.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = power
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = ${fuel_blocks}
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = transfer
  []
[]

We add a transient executioner - again, even though the MOOSE heat conduction module in this tutorial solves steady equations, a transient executioner allows us to control the frequency with which MOOSE and NekRS iterate the CHT physics. Here, we use a time step that is 50 times bigger than the NekRS time step.

[Executioner<<<{"href": "../syntax/Executioner/index.html"}>>>]
  type = Transient
  dt = '${fparse N * nek_dt * t0}'
  nl_abs_tol = 1e-5
  nl_rel_tol = 1e-16
  petsc_options_value = 'hypre boomeramg'
  petsc_options_iname = '-pc_type -pc_hypre_type'
[]

Finally, we add a number of LayeredAverage user objects to compute averages of the fuel and graphite temperatures in the axial direction, which are output to CSV using a SpatialUserObjectVectorPostprocessor, in combination with a CSV output.

[UserObjects<<<{"href": "../syntax/UserObjects/index.html"}>>>]
  [average_fuel_axial]
    type = LayeredAverage<<<{"description": "Computes averages of variables over layers", "href": "../source/userobjects/LayeredAverage.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = T
    direction<<<{"description": "The direction of the layers."}>>> = z
    num_layers<<<{"description": "The number of layers."}>>> = ${num_layers_for_plots}
    block<<<{"description": "The list of block ids (SubdomainID) that this object will be applied"}>>> = ${fuel_blocks}
  []
  [average_block_axial]
    type = LayeredAverage<<<{"description": "Computes averages of variables over layers", "href": "../source/userobjects/LayeredAverage.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = T
    direction<<<{"description": "The direction of the layers."}>>> = z
    num_layers<<<{"description": "The number of layers."}>>> = ${num_layers_for_plots}
    block<<<{"description": "The list of block ids (SubdomainID) that this object will be applied"}>>> = 'graphite'
  []
[]

[VectorPostprocessors<<<{"href": "../syntax/VectorPostprocessors/index.html"}>>>]
  [fuel_axial_avg]
    type = SpatialUserObjectVectorPostprocessor<<<{"description": "Outputs the values of a spatial user object in the order of the specified spatial points", "href": "../source/vectorpostprocessors/SpatialUserObjectVectorPostprocessor.html"}>>>
    userobject<<<{"description": "The userobject whose values are to be reported"}>>> = average_fuel_axial
  []
  [block_axial_avg]
    type = SpatialUserObjectVectorPostprocessor<<<{"description": "Outputs the values of a spatial user object in the order of the specified spatial points", "href": "../source/vectorpostprocessors/SpatialUserObjectVectorPostprocessor.html"}>>>
    userobject<<<{"description": "The userobject whose values are to be reported"}>>> = average_block_axial
  []
[]

[Outputs<<<{"href": "../syntax/Outputs/index.html"}>>>]
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
  hide<<<{"description": "A list of the variables and postprocessors that should NOT be output to the Exodus file (may include Variables, ScalarVariables, and Postprocessor names)."}>>> = 'synchronization_to_nek'

  [csv]
    file_base<<<{"description": "The desired solution output name without an extension. If not provided, MOOSE sets it with Outputs/file_base when available. Otherwise, MOOSE uses input file name and this object name for a master input or uses master file_base, the subapp name and this object name for a subapp input to set it."}>>> = 'csv/solid_nek'
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../source/outputs/CSV.html"}>>>
  []
[]

Fluid Input Files

The fluid mass, momentum, and energy transport physics are solved using NekRS. The input file for this portion of the physics is the nek.i input. We begin by defining a number of constants and by setting up the NekRSMesh mesh mirror. Because we are coupling via boundary CHT to MOOSE, we set boundary = '3' so that we will be able to extract the boundary temperature from boundary 3 (the wall). We are also coupling via volumes to OpenMC higher in the MultiApp hierarchy. In order to extract volume representations of the fluid temperature, we also set volume = true.

# copy-pasta from common_input.i

channel_diameter = 0.016                 # diameter of the coolant channels (m)
height = 6.343                           # height of the full core (m)
inlet_T = 598.0                          # inlet fluid temperature (K)
power = 200e6                            # full core power (W)
mdot = 117.3                             # fluid mass flowrate (kg/s)
fluid_density = 5.5508                   # fluid density (kg/m3)
fluid_Cp = 5189.0                        # fluid isobaric specific heat (J/kg/K)
n_bundles = 12                           # number of bundles in the full core
n_coolant_channels_per_block = 108       # number of coolant channels per assembly
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected

[Mesh<<<{"href": "../syntax/Mesh/index.html"}>>>]
  type = NekRSMesh
  boundary = '3'
  volume = true
  scaling = ${channel_diameter}
[]

The bulk of the NekRS wrapping occurs in the [Problem] block with NekRSProblem. The NekRS input files are in non-dimensional form, whereas all other coupled applications use dimensional units. The various *_ref and *_0 parameters define the characteristic scales that were used to non-dimensionalize the NekRS input. We add the NekBoundaryFlux and NekFieldVariable transfers in order to write heat flux and read temperature from NekRS. Finally, we indicate that we will be minimizing the data transfers in/out of NekRS unless new data is actually available from the MOOSE heat transfer module with the synchronization_interval = parent_app parameter.

[Problem<<<{"href": "../syntax/Problem/index.html"}>>>]
  type = NekRSProblem
  casename<<<{"description": "Case name for the NekRS input files; this is <case> in <case>.par, <case>.udf, <case>.oudf, and <case>.re2."}>>> = 'ranstube'
  n_usrwrk_slots<<<{"description": "Number of slots to allocate in nrs->usrwrk to hold fields either related to coupling (which will be populated by Cardinal), or other custom usages, such as a distance-to-wall calculation"}>>> = 2

  [Dimensionalize<<<{"href": "../syntax/Problem/Dimensionalize/index.html"}>>>]
    U<<<{"description": "Reference velocity"}>>> = '${fparse mdot / (n_bundles * n_coolant_channels_per_block) / fluid_density / (pi * channel_diameter * channel_diameter / 4.0)}'
    T<<<{"description": "Reference temperature"}>>> = ${inlet_T}
    dT<<<{"description": "Reference temperature difference"}>>> = '${fparse power / mdot / fluid_Cp * unit_cell_height / height}'
    L<<<{"description": "Reference length"}>>> = ${channel_diameter}
    rho<<<{"description": "Reference density"}>>> = ${fluid_density}
    Cp<<<{"description": "Reference isobaric specific heat"}>>> = ${fluid_Cp}
  []

  [FieldTransfers<<<{"href": "../syntax/Problem/FieldTransfers/index.html"}>>>]
    [flux]
      type = NekBoundaryFlux<<<{"description": "Reads/writes boundary flux data between NekRS and MOOSE."}>>>
      direction<<<{"description": "Direction in which to send data"}>>> = to_nek
      usrwrk_slot<<<{"description": "When 'direction = to_nek', the slot(s) in the usrwrk array to write the incoming data; provide one entry for each quantity being passed"}>>> = 0
    []
    [temperature]
      type = NekFieldVariable<<<{"description": "Reads/writes volumetric field data between NekRS and MOOSE."}>>>
      direction<<<{"description": "Direction in which to send data"}>>> = from_nek
    []
  []

  synchronization_interval = parent_app
[]

Next, we will allow NekRS to select its own time step using the NekTimeStepper, combined with a transient executioner.

[Executioner<<<{"href": "../syntax/Executioner/index.html"}>>>]
  type = Transient

  [TimeStepper<<<{"href": "../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = NekTimeStepper
  []
[]

We also add a number of postprocessors to query the Nek solution.

[Postprocessors<<<{"href": "../syntax/Postprocessors/index.html"}>>>]
  [inlet_T]
    type = NekSideAverage<<<{"description": "Average of a field over a boundary of the NekRS mesh", "href": "../source/postprocessors/NekSideAverage.html"}>>>
    field<<<{"description": "Field to apply this object to"}>>> = temperature
    boundary<<<{"description": "Boundary ID(s) for which to compute the postprocessor"}>>> = '1'
  []
  [outlet_T]
    type = NekSideAverage<<<{"description": "Average of a field over a boundary of the NekRS mesh", "href": "../source/postprocessors/NekSideAverage.html"}>>>
    field<<<{"description": "Field to apply this object to"}>>> = temperature
    boundary<<<{"description": "Boundary ID(s) for which to compute the postprocessor"}>>> = '2'
  []
  [max_T]
    type = NekVolumeExtremeValue<<<{"description": "Extreme value (max/min) of a field over the NekRS mesh", "href": "../source/postprocessors/NekVolumeExtremeValue.html"}>>>
    field<<<{"description": "Field to apply this object to"}>>> = temperature
    value_type<<<{"description": "Whether to give the maximum or minimum extreme value"}>>> = max
  []
[]

Finally, we define the output formats and hide the automatically-created flux_integral and transfer_in postprocessors from the screen (console) to have neater output.

[Outputs<<<{"href": "../syntax/Outputs/index.html"}>>>]
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true

  [screen]
    type = Console<<<{"description": "Object for screen output.", "href": "../source/outputs/Console.html"}>>>
    hide<<<{"description": "A list of the variables and postprocessors that should NOT be output to the Exodus file (may include Variables, ScalarVariables, and Postprocessor names)."}>>> = 'flux_integral transfer_in'
  []

  [csv]
    file_base<<<{"description": "The desired solution output name without an extension. If not provided, MOOSE sets it with Outputs/file_base when available. Otherwise, MOOSE uses input file name and this object name for a master input or uses master file_base, the subapp name and this object name for a subapp input to set it."}>>> = 'csv/nek'
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../source/outputs/CSV.html"}>>>
  []
[]

OpenMC Input Files

The neutronics physics is solved over the entire domain with OpenMC. The OpenMC wrapping used for the OpenMC-THM-MOOSE coupling is described in the openm_thm.i input file. Relative to the OpenMC input files used for the OpenMC-NekRS-MOOSE coupling shown previously, we now need to add variables to hold the wall temperature, thm_temp_wall, that THM passes to the MOOSE heat transfer module, and the wall heat flux, flux, that the MOOSE heat transfer module passes to THM. Both of these are simply receiver variables that are written into by the two sub-applications (MOOSE heat conduction and THM).

[AuxVariables<<<{"href": "../syntax/AuxVariables/index.html"}>>>]
  [cell_temperature]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [cell_density]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [thm_temp_wall]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
    block = ${density_blocks}
  []
  [flux]
  []

  [q_wall]
  []
[]

Next, when using the "tree" MultiApp hierarchy, OpenMC runs after the MOOSE heat transfer module, but before THM. Therefore, initial conditions are required for the fluid temperature (for OpenMC) as well as for the fluid wall temperature (which will be sent to the MOOSE heat transfer module on the first time step). Because OpenMC sends its heat source to the MOOSE heat transfer module at the beginning of a time step, we also need to set an initial condition on the heat source.

[ICs<<<{"href": "../syntax/Cardinal/ICs/index.html"}>>>]
  [fluid_temp_wall]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = thm_temp_wall
    function<<<{"description": "The initial condition function."}>>> = temp_ic
  []
  [fluid_temp]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../source/ics/FunctionIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = thm_temp
    function<<<{"description": "The initial condition function."}>>> = temp_ic
  []
  [heat_source]
    type = ConstantIC<<<{"description": "Sets a constant field value.", "href": "../source/ics/ConstantIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = heat_source
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = ${fuel_blocks}
    value<<<{"description": "The value to be set in IC"}>>> = ${fparse unit_cell_power / (2.0 * pi * compact_diameter * compact_diameter / 4.0 * unit_cell_height)}
  []
[]

[Functions<<<{"href": "../syntax/Functions/index.html"}>>>]
  [temp_ic]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = '${inlet_T} + z / ${unit_cell_height} * ${unit_cell_power} / (${mdot} / ${n_bundles} / ${n_coolant_channels_per_block}) / ${fluid_Cp}'
  []
[]

The [Problem] block is identical to that shown earlier for the OpenMC-NekRS-MOOSE coupling. Next, we define the two sub-applications and the transfers. The transfers of wall heat flux and wall temperature between MOOSE and THM first pass up through OpenMC to one of the sub-applications, giving four additional transfers in the main application's input file than we saw for the "single-stack" hierarchy.

[MultiApps<<<{"href": "../syntax/MultiApps/index.html"}>>>]
  [bison]
    type = TransientMultiApp<<<{"description": "MultiApp for performing coupled simulations with the parent and sub-application both progressing in time.", "href": "../source/multiapps/TransientMultiApp.html"}>>>
    input_files<<<{"description": "The input file for each App.  If this parameter only contains one input file it will be used for all of the Apps.  When using 'positions_from_file' it is also admissable to provide one input_file per file."}>>> = 'solid_thm.i'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_begin
  []
  [thm]
    type = FullSolveMultiApp<<<{"description": "Performs a complete simulation during each execution.", "href": "../source/multiapps/FullSolveMultiApp.html"}>>>
    input_files<<<{"description": "The input file for each App.  If this parameter only contains one input file it will be used for all of the Apps.  When using 'positions_from_file' it is also admissable to provide one input_file per file."}>>> = 'thm.i'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
    max_procs_per_app<<<{"description": "Maximum number of processors to give to each App in this MultiApp.  Useful for restricting small solves to just a few procs so they don't get spread out"}>>> = 1
    bounding_box_padding<<<{"description": "Additional padding added to the dimensions of the bounding box. The values are added to the x, y and z dimension respectively."}>>> = '0.1 0.1 0'
  []
[]

[Transfers<<<{"href": "../syntax/Transfers/index.html"}>>>]
  [solid_temp_to_openmc]
    type = MultiAppGeometricInterpolationTransfer<<<{"description": "Transfers the value to the target domain from a combination/interpolation of the values on the nearest nodes in the source domain, using coefficients based on the distance to each node.", "href": "../source/transfers/MultiAppGeometricInterpolationTransfer.html"}>>>
    source_variable<<<{"description": "The variable to transfer from."}>>> = T
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = solid_temp
    from_multi_app<<<{"description": "The name of the MultiApp to receive data from"}>>> = bison
  []
  [heat_flux_to_openmc]
    type = MultiAppGeneralFieldNearestLocationTransfer<<<{"description": "Transfers field data at the MultiApp position by finding the value at the nearest neighbor(s) in the origin application.", "href": "../source/transfers/MultiAppGeneralFieldNearestLocationTransfer.html"}>>>
    source_variable<<<{"description": "The variable to transfer from."}>>> = flux
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = flux
    from_multi_app<<<{"description": "The name of the MultiApp to receive data from"}>>> = bison
    from_boundaries<<<{"description": "The boundary we are transferring from (if not specified, whole domain is used)."}>>> = 'fluid_solid_interface'
    to_boundaries<<<{"description": "The boundary we are transferring to (if not specified, whole domain is used)."}>>> = 'fluid_solid_interface'
    from_postprocessors_to_be_preserved<<<{"description": "The name of the Postprocessor in the from-app to evaluate an adjusting factor."}>>> = flux_integral
    to_postprocessors_to_be_preserved<<<{"description": "The name of the Postprocessor in the to-app to evaluate an adjusting factor."}>>> = flux_integral
  []
  [source_to_bison]
    type = MultiAppGeneralFieldShapeEvaluationTransfer<<<{"description": "Transfers field data at the MultiApp position using the finite element shape functions from the origin application.", "href": "../source/transfers/MultiAppGeneralFieldShapeEvaluationTransfer.html"}>>>
    source_variable<<<{"description": "The variable to transfer from."}>>> = heat_source
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = power
    to_multi_app<<<{"description": "The name of the MultiApp to transfer the data to"}>>> = bison
    from_postprocessors_to_be_preserved<<<{"description": "The name of the Postprocessor in the from-app to evaluate an adjusting factor."}>>> = heat_source
    to_postprocessors_to_be_preserved<<<{"description": "The name of the Postprocessor in the to-app to evaluate an adjusting factor."}>>> = power
  []
  [thm_temp_to_bison]
    type = MultiAppGeometricInterpolationTransfer<<<{"description": "Transfers the value to the target domain from a combination/interpolation of the values on the nearest nodes in the source domain, using coefficients based on the distance to each node.", "href": "../source/transfers/MultiAppGeometricInterpolationTransfer.html"}>>>
    source_variable<<<{"description": "The variable to transfer from."}>>> = thm_temp_wall
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = fluid_temp
    to_multi_app<<<{"description": "The name of the MultiApp to transfer the data to"}>>> = bison
  []

  [q_wall_to_thm]
    type = MultiAppGeneralFieldUserObjectTransfer<<<{"description": "Transfers user object spatial evaluations from an origin app onto a variable in the target application.", "href": "../source/transfers/MultiAppGeneralFieldUserObjectTransfer.html"}>>>
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = q_wall
    to_multi_app<<<{"description": "The name of the MultiApp to transfer the data to"}>>> = thm
    source_user_object<<<{"description": "The UserObject you want to transfer values from. It must implement the SpatialValue() class routine"}>>> = q_wall_avg
  []
  [T_wall_from_thm]
    type = MultiAppGeneralFieldNearestLocationTransfer<<<{"description": "Transfers field data at the MultiApp position by finding the value at the nearest neighbor(s) in the origin application.", "href": "../source/transfers/MultiAppGeneralFieldNearestLocationTransfer.html"}>>>
    source_variable<<<{"description": "The variable to transfer from."}>>> = T_wall
    from_multi_app<<<{"description": "The name of the MultiApp to receive data from"}>>> = thm
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = thm_temp_wall
  []
  [T_bulk_from_thm]
    type = MultiAppGeneralFieldNearestLocationTransfer<<<{"description": "Transfers field data at the MultiApp position by finding the value at the nearest neighbor(s) in the origin application.", "href": "../source/transfers/MultiAppGeneralFieldNearestLocationTransfer.html"}>>>
    source_variable<<<{"description": "The variable to transfer from."}>>> = T
    from_multi_app<<<{"description": "The name of the MultiApp to receive data from"}>>> = thm
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = thm_temp
  []
[]

Finally, because we are coupling to THM, we need to average the wall heat flux around the coolant channel into a number of layers to transfer from the 3-D MOOSE heat conduction model to the 1-D fluid flow model. Heat flux is averaged using a LayeredSideAverage.

[UserObjects<<<{"href": "../syntax/UserObjects/index.html"}>>>]
  [q_wall_avg]
    type = LayeredSideAverage<<<{"description": "Computes side averages of a variable storing partial sums for the specified number of intervals in a direction (x,y,z).", "href": "../source/userobjects/LayeredSideAverage.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = fluid_solid_interface
    variable<<<{"description": "The name of the variable that this boundary condition applies to"}>>> = flux

    # Note: make this to match the num_elems in the channel
    direction<<<{"description": "The direction of the layers."}>>> = z
    num_layers<<<{"description": "The number of layers."}>>> = ${num_layers_for_THM}

    direction_min<<<{"description": "Minimum coordinate along 'direction' that bounds the layers"}>>> = 0.0
    direction_max<<<{"description": "Maximum coordinate along 'direction' that bounds the layers"}>>> = ${unit_cell_height}
  []
  [average_power_axial]
    type = LayeredAverage<<<{"description": "Computes averages of variables over layers", "href": "../source/userobjects/LayeredAverage.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = heat_source
    direction<<<{"description": "The direction of the layers."}>>> = z
    num_layers<<<{"description": "The number of layers."}>>> = ${num_layers}
    block<<<{"description": "The list of block ids (SubdomainID) that this object will be applied"}>>> = ${fuel_blocks}
  []
[]

All other aspects of the input file are the same as for the OpenMC-NekRS-MOOSE case.

Solid Input Files

The solid heat conduction input file is the solid_thm.i input. This input file is identical to the solid_nek.i input file for the OpenMC-NekRS-MOOSE case except for:

• There are no MultiApps, and therefore no transfers, to a thermal-fluid code

• There are no initial conditions in the input file, since OpenMC sends the MOOSE heat conduction file all necessary initial conditions for wall temperature and heat source

• There is no need to explicitly turn off the material coverage check or add a NullKernel because there is not a dummy fluid block to receive fluid temperature from a sub-application

Because all of these differences are simply omissions, this concludes the discussion of the solid input file. For reference, the full file is below.

# copy-pasta from common_input.i
inlet_T = 598.0                          # inlet fluid temperature (K)
buffer_k = 0.5                           # buffer thermal conductivity (W/m/K)
PyC_k = 4.0                              # PyC thermal conductivity (W/m/K)
SiC_k = 13.9                             # SiC thermal conductivity (W/m/K)
kernel_k = 3.5                           # fissil kernel thermal conductivity (W/m/K)
matrix_k = 15.0                          # graphite matrix thermal conductivity (W/m/K)
num_layers_for_plots = 50                # number of layers to average fields over for plotting
triso_pf = 0.15                          # TRISO packing fraction (%)
kernel_radius = 214.85e-6                # fissile kernel outer radius (m)
buffer_radius = 314.85e-6                # buffer outer radius (m)
iPyC_radius = 354.85e-6                  # inner PyC outer radius (m)
SiC_radius = 389.85e-6                   # SiC outer radius (m)
oPyC_radius = 429.85e-6                  # outer PyC outer radius (m)

# compute the volume fraction of each TRISO layer in a TRISO particle
# for use in computing average thermophysical properties
kernel_fraction = ${fparse kernel_radius^3 / oPyC_radius^3}
buffer_fraction = ${fparse (buffer_radius^3 - kernel_radius^3) / oPyC_radius^3}
ipyc_fraction = ${fparse (iPyC_radius^3 - buffer_radius^3) / oPyC_radius^3}
sic_fraction = ${fparse (SiC_radius^3 - iPyC_radius^3) / oPyC_radius^3}
opyc_fraction = ${fparse (oPyC_radius^3 - SiC_radius^3) / oPyC_radius^3}

fuel_blocks = 'compacts compacts_trimmer_tri'

[Mesh<<<{"href": "../syntax/Mesh/index.html"}>>>]
  [mesh]
   type = FileMeshGenerator<<<{"description": "Read a mesh from a file.", "href": "../source/meshgenerators/FileMeshGenerator.html"}>>>
   file<<<{"description": "The filename to read."}>>> = solid_mesh_in.e
  []
  [delete_coolant]
    type = BlockDeletionGenerator<<<{"description": "Mesh generator which removes elements from the specified subdomains", "href": "../source/meshgenerators/BlockDeletionGenerator.html"}>>>
    input<<<{"description": "The mesh we want to modify"}>>> = mesh
    block<<<{"description": "The list of blocks to be deleted"}>>> = 'coolant'
  []
[]

[Variables<<<{"href": "../syntax/Variables/index.html"}>>>]
  [T]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = ${inlet_T}
  []
[]

[Kernels<<<{"href": "../syntax/Kernels/index.html"}>>>]
  [diffusion]
    type = HeatConduction<<<{"description": "Diffusive heat conduction term $-\\nabla\\cdot(k\\nabla T)$ of the thermal energy conservation equation", "href": "../source/kernels/HeatConduction.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
  []
  [source]
    type = CoupledForce<<<{"description": "Implements a source term proportional to the value of a coupled variable. Weak form: $(\\psi_i, -\\sigma v)$.", "href": "../source/kernels/CoupledForce.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
    v<<<{"description": "The coupled variable which provides the force"}>>> = power
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = ${fuel_blocks}
  []
[]

[BCs<<<{"href": "../syntax/BCs/index.html"}>>>]
  [pin_outer]
    type = MatchedValueBC<<<{"description": "Implements a NodalBC which equates two different Variables' values on a specified boundary.", "href": "../source/bcs/MatchedValueBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
    v<<<{"description": "The variable whose value we are to match."}>>> = fluid_temp
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'fluid_solid_interface'
  []
[]

[Functions<<<{"href": "../syntax/Functions/index.html"}>>>]
  [k_graphite]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = '${matrix_k}'
  []
  [k_TRISO]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = '${kernel_fraction} * ${kernel_k} + ${buffer_fraction} * ${buffer_k} + ${fparse ipyc_fraction + opyc_fraction} * ${PyC_k} + ${sic_fraction} * ${SiC_k}'
  []
  [k_compacts]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = '${triso_pf} * k_TRISO + ${fparse 1.0 - triso_pf} * k_graphite'
    symbol_names<<<{"description": "Symbols (excluding t,x,y,z) that are bound to the values provided by the corresponding items in the vals vector."}>>> = 'k_TRISO k_graphite'
    symbol_values<<<{"description": "Constant numeric values, postprocessor names, function names, and scalar variables corresponding to the symbols in symbol_names."}>>> = 'k_TRISO k_graphite'
  []
[]

[Materials<<<{"href": "../syntax/Materials/index.html"}>>>]
  [graphite]
    type = HeatConductionMaterial<<<{"description": "General-purpose material model for heat conduction", "href": "../source/materials/HeatConductionMaterial.html"}>>>
    thermal_conductivity_temperature_function<<<{"description": "Thermal conductivity as a function of temperature."}>>> = k_graphite
    temp<<<{"description": "Coupled Temperature"}>>> = T
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 'graphite'
  []
  [compacts]
    type = HeatConductionMaterial<<<{"description": "General-purpose material model for heat conduction", "href": "../source/materials/HeatConductionMaterial.html"}>>>
    thermal_conductivity_temperature_function<<<{"description": "Thermal conductivity as a function of temperature."}>>> = k_compacts
    temp<<<{"description": "Coupled Temperature"}>>> = T
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = ${fuel_blocks}
  []
[]

[Postprocessors<<<{"href": "../syntax/Postprocessors/index.html"}>>>]
  [flux_integral] # evaluate the total heat flux for normalization
    type = SideDiffusiveFluxIntegral<<<{"description": "Computes the integral of the diffusive flux over the specified boundary", "href": "../source/postprocessors/SideDiffusiveFluxIntegral.html"}>>>
    diffusivity<<<{"description": "The name of the diffusivity material property that will be used in the flux computation. This must be provided if the variable is of finite element type"}>>> = thermal_conductivity
    variable<<<{"description": "The name of the variable which this postprocessor integrates"}>>> = T
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'fluid_solid_interface'
  []
  [max_fuel_T]
    type = ElementExtremeValue<<<{"description": "Finds either the min or max elemental value of a variable over the domain.", "href": "../source/postprocessors/ElementExtremeValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on"}>>> = T
    value_type<<<{"description": "Type of extreme value to return. 'max' returns the maximum value. 'min' returns the minimum value. 'max_abs' returns the maximum of the absolute value."}>>> = max
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = ${fuel_blocks}
  []
  [max_block_T]
    type = ElementExtremeValue<<<{"description": "Finds either the min or max elemental value of a variable over the domain.", "href": "../source/postprocessors/ElementExtremeValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on"}>>> = T
    value_type<<<{"description": "Type of extreme value to return. 'max' returns the maximum value. 'min' returns the minimum value. 'max_abs' returns the maximum of the absolute value."}>>> = max
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 'graphite'
  []
  [power]
    type = ElementIntegralVariablePostprocessor<<<{"description": "Computes a volume integral of the specified variable", "href": "../source/postprocessors/ElementIntegralVariablePostprocessor.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = power
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = ${fuel_blocks}
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'transfer'
  []
[]

[AuxVariables<<<{"href": "../syntax/AuxVariables/index.html"}>>>]
  [fluid_temp]
  []
  [flux]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
  [power]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
[]

[AuxKernels<<<{"href": "../syntax/AuxKernels/index.html"}>>>]
  [flux]
    type = DiffusionFluxAux<<<{"description": "Compute components of flux vector for diffusion problems $(\\vec{J} = -D \\nabla C)$.", "href": "../source/auxkernels/DiffusionFluxAux.html"}>>>
    diffusion_variable<<<{"description": "The name of the variable"}>>> = T
    component<<<{"description": "The desired component of flux."}>>> = normal
    diffusivity<<<{"description": "The name of the diffusivity material property that will be used in the flux computation."}>>> = thermal_conductivity
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = flux
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 'fluid_solid_interface'
  []
[]

[Executioner<<<{"href": "../syntax/Executioner/index.html"}>>>]
  type = Transient
  nl_abs_tol = 1e-5
  nl_rel_tol = 1e-16
  petsc_options_value = 'hypre boomeramg'
  petsc_options_iname = '-pc_type -pc_hypre_type'
[]

[UserObjects<<<{"href": "../syntax/UserObjects/index.html"}>>>]
  [average_fuel_axial]
    type = LayeredAverage<<<{"description": "Computes averages of variables over layers", "href": "../source/userobjects/LayeredAverage.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = T
    direction<<<{"description": "The direction of the layers."}>>> = z
    num_layers<<<{"description": "The number of layers."}>>> = ${num_layers_for_plots}
    block<<<{"description": "The list of block ids (SubdomainID) that this object will be applied"}>>> = ${fuel_blocks}
  []
  [average_block_axial]
    type = LayeredAverage<<<{"description": "Computes averages of variables over layers", "href": "../source/userobjects/LayeredAverage.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = T
    direction<<<{"description": "The direction of the layers."}>>> = z
    num_layers<<<{"description": "The number of layers."}>>> = ${num_layers_for_plots}
    block<<<{"description": "The list of block ids (SubdomainID) that this object will be applied"}>>> = 'graphite'
  []
[]

[VectorPostprocessors<<<{"href": "../syntax/VectorPostprocessors/index.html"}>>>]
  [fuel_axial_avg]
    type = SpatialUserObjectVectorPostprocessor<<<{"description": "Outputs the values of a spatial user object in the order of the specified spatial points", "href": "../source/vectorpostprocessors/SpatialUserObjectVectorPostprocessor.html"}>>>
    userobject<<<{"description": "The userobject whose values are to be reported"}>>> = average_fuel_axial
  []
  [block_axial_avg]
    type = SpatialUserObjectVectorPostprocessor<<<{"description": "Outputs the values of a spatial user object in the order of the specified spatial points", "href": "../source/vectorpostprocessors/SpatialUserObjectVectorPostprocessor.html"}>>>
    userobject<<<{"description": "The userobject whose values are to be reported"}>>> = average_block_axial
  []
[]

[Outputs<<<{"href": "../syntax/Outputs/index.html"}>>>]
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
  print_linear_residuals<<<{"description": "Enable printing of linear residuals to the screen (Console)"}>>> = false

  [csv]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../source/outputs/CSV.html"}>>>
    file_base<<<{"description": "The desired solution output name without an extension. If not provided, MOOSE sets it with Outputs/file_base when available. Otherwise, MOOSE uses input file name and this object name for a master input or uses master file_base, the subapp name and this object name for a subapp input to set it."}>>> = 'csv_thm/solid_out'
  []
[]

Fluid Input Files

The fluid mass, momentum, and energy transport physics are solved using THM. The input file for this portion of the physics is the thm.i input. The THM input file is built using syntax specific to THM - we will only briefly cover the syntax, and instead refer users to the THM documentation for more information. First, we define a number of constants at the beginning of the file and apply some global settings. We set the initial conditions for pressure, velocity, and temperature and indicate the fluid EOS object using IdealGasFluidProperties.

# copy-pasta from common.i
inlet_T = 598.0                          # inlet fluid temperature (K)
mdot = ${fparse 117.3 / 12 / 108}        # fluid mass flowrate (kg/s)
outlet_P = 7.1e6                         # fluid outlet pressure (Pa)
channel_diameter = 0.016                 # diameter of the coolant channels (m)
unit_cell_height = 1.6                   # unit cell height - arbitrarily selected

num_layers_for_THM = 150

[GlobalParams<<<{"href": "../syntax/GlobalParams/index.html"}>>>]
  initial_p = ${outlet_P}
  initial_T = ${inlet_T}
  initial_vel = ${fparse mdot / outlet_P / 8.3144598 * 4.0e-3 / inlet_T / (pi * channel_diameter * channel_diameter / 4.0)}

  rdg_slope_reconstruction = full
  closures = none
  fp = helium
[]

[Closures<<<{"href": "../syntax/Closures/index.html"}>>>]
  [none]
    type = WallTemperature1PhaseClosures<<<{"description": "Adds wall temperature material for single-phase flow.", "href": "../source/closures/WallTemperature1PhaseClosures.html"}>>>
  []
[]

[FluidProperties<<<{"href": "../syntax/FluidProperties/index.html"}>>>]
  [helium]
    type = IdealGasFluidProperties<<<{"description": "Fluid properties for an ideal gas", "href": "../source/fluidproperties/IdealGasFluidProperties.html"}>>>
    molar_mass<<<{"description": "Constant molar mass of the fluid (kg/mol)"}>>> = 4e-3
    gamma<<<{"description": "gamma value (cp/cv)"}>>> = 1.668282 # should correspond to  Cp = 5189 J/kg/K
    k<<<{"description": "Thermal conductivity, W/(m-K)"}>>> = 0.2556
    mu<<<{"description": "Dynamic viscosity, Pa.s"}>>> = 3.22639e-5
  []
[]

Next, we define the "components" in the domain. These components essentially consist of the physics equations and boundary conditions solved by THM, but expressed in THM-specific syntax. These components define single-phase flow in a pipe, an inlet mass flowrate boundary condition, an outlet pressure condition, and heat transfer to the pipe wall.

[Components<<<{"href": "../syntax/Components/index.html"}>>>]
  [channel]
    type = FlowChannel1Phase<<<{"description": "1-phase 1D flow channel", "href": "../source/components/FlowChannel1Phase.html"}>>>
    position<<<{"description": "Start position of axis in 3-D space [m]"}>>> = '0 0 0'
    orientation<<<{"description": "Direction of flow channel from start position to end position (no need to normalize). For curved flow channels, it is the (tangent) direction at the start position."}>>> = '0 0 1'

    A<<<{"description": "Area of the flow channel, can be a constant or a function"}>>> = '${fparse pi * channel_diameter * channel_diameter / 4}'
    D_h<<<{"description": "Hydraulic diameter [m]"}>>> = ${channel_diameter}
    length<<<{"description": "Length of each axial section [m]"}>>> = ${unit_cell_height}
    n_elems<<<{"description": "Number of elements in each axial section"}>>> = ${num_layers_for_THM}
  []

  [inlet]
    type = InletMassFlowRateTemperature1Phase<<<{"description": "Boundary condition with prescribed mass flow rate and temperature for 1-phase flow channels.", "href": "../source/components/InletMassFlowRateTemperature1Phase.html"}>>>
    input<<<{"description": "Name of the input"}>>> = 'channel:in'
    m_dot<<<{"description": "Prescribed mass flow rate [kg/s]"}>>> = ${mdot}
    T<<<{"description": "Prescribed temperature [K]"}>>> = ${inlet_T}
  []

  [outlet]
    type = Outlet1Phase<<<{"description": "Boundary condition with prescribed pressure for 1-phase flow channels.", "href": "../source/components/Outlet1Phase.html"}>>>
    input<<<{"description": "Name of the input"}>>> = 'channel:out'
    p<<<{"description": "Prescribed pressure [Pa]"}>>> = ${outlet_P}
  []

  [ht_ext]
    type = HeatTransferFromExternalAppHeatFlux1Phase<<<{"description": "Heat transfer specified by heat flux provided by an external application going into 1-phase flow channel.", "href": "../source/components/HeatTransferFromExternalAppHeatFlux1Phase.html"}>>>
    flow_channel<<<{"description": "Name of flow channel component to connect to"}>>> = channel
    P_hf<<<{"description": "Heat flux perimeter [m]"}>>> = '${fparse channel_diameter * pi}'
  []
[]

Associated with these components are a number of closures, defined as materials. We set up the Churchill correlation for the friction factor and the Dittus-Boelter correlation for the convective heat transfer coefficient. For the Dittus-Boelter correlation, we use a corrected version of the closure (with the leading coefficient changed from 0.023 to 0.021) based on the NekRS simulations described in an earlier tutorial. Additional materials are created to represent dimensionless numbers and other auxiliary terms, such as the wall temperature. As can be seen here, the Material system is not always used to represent quantities traditioanlly thought of as "material properties."

[Materials<<<{"href": "../syntax/Materials/index.html"}>>>]
  # wall friction closure
  [f_mat]
    type = ADWallFrictionChurchillMaterial<<<{"description": "Computes the Darcy friction factor using the Churchill correlation.", "href": "../source/materials/ADWallFrictionChurchillMaterial.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = channel
    D_h<<<{"description": "hydraulic diameter"}>>> = D_h
    f_D<<<{"description": "Darcy friction factor material property"}>>> = f_D
    vel<<<{"description": "x-component of the velocity"}>>> = vel
    rho<<<{"description": "Density"}>>> = rho
    mu<<<{"description": "Dynamic viscosity material property"}>>> = mu
  []

  # Wall heat transfer closure (all important is in Nu_mat)
  [Re_mat]
    type = ADReynoldsNumberMaterial<<<{"description": "Computes Reynolds number as a material property", "href": "../source/materials/ADReynoldsNumberMaterial.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = channel
    Re<<<{"description": "Reynolds number property name"}>>> = Re
    D_h<<<{"description": "Hydraulic diameter"}>>> = D_h
    mu<<<{"description": "Dynamic viscosity of the phase"}>>> = mu
    vel<<<{"description": "Velocity of the phase"}>>> = vel
    rho<<<{"description": "Density of the phase"}>>> = rho
  []
  [Pr_mat]
    type = ADPrandtlNumberMaterial<<<{"description": "Computes Prandtl number as material property", "href": "../source/materials/ADPrandtlNumberMaterial.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = channel
    cp<<<{"description": "Constant-pressure specific heat"}>>> = cp
    mu<<<{"description": "Dynamic viscosity"}>>> = mu
    k<<<{"description": "Thermal conductivity"}>>> = k
  []
  [Nu_mat]
    type = ADParsedMaterial<<<{"description": "Parsed expression Material.", "href": "../source/materials/ParsedMaterial.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = channel
    # Dittus-Boelter
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '0.021 * pow(Re, 0.8) * pow(Pr, 0.4)'
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'Nu'
    material_property_names<<<{"description": "Vector of material properties used in the parsed function"}>>> = 'Re Pr'
  []
  [Hw_mat]
    type = ADConvectiveHeatTransferCoefficientMaterial<<<{"description": "Computes convective heat transfer coefficient from Nusselt number", "href": "../source/materials/ADConvectiveHeatTransferCoefficientMaterial.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = channel
    D_h<<<{"description": "Hydraulic diameter"}>>> = D_h
    Nu<<<{"description": "Nusselt number"}>>> = Nu
    k<<<{"description": "Thermal conductivity"}>>> = k
  []
  [T_wall]
    type = ADTemperatureWall3EqnMaterial<<<{"description": "Computes the wall temperature from the fluid temperature, the heat flux and the heat transfer coefficient", "href": "../source/materials/ADTemperatureWall3EqnMaterial.html"}>>>
    Hw<<<{"description": "Heat transfer coefficient"}>>> = Hw
    T<<<{"description": "Fluid temperature"}>>> = T
    q_wall<<<{"description": "Wall heat flux"}>>> = q_wall
  []
[]

THM computes the wall temperature to apply a boundary condition in the MOOSE heat transfer module. To convert the T_wall material into an auxiliary variable, we use the ADMaterialRealAux.

[AuxVariables<<<{"href": "../syntax/AuxVariables/index.html"}>>>]
  [T_wall]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
[]

[AuxKernels<<<{"href": "../syntax/AuxKernels/index.html"}>>>]
  [Tw_aux]
    type = ADMaterialRealAux<<<{"description": "Outputs element volume-averaged material properties", "href": "../source/auxkernels/MaterialRealAux.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = channel
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = T_wall
    property<<<{"description": "The material property name."}>>> = T_wall
  []
[]

Finally, we set the preconditioner, a Transient executioner, and set an Exodus output. We will run THM to convergence based on a tight steady state relative tolerance of .

[Preconditioning<<<{"href": "../syntax/Preconditioning/index.html"}>>>]
  [pc]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
  []
[]

[Executioner<<<{"href": "../syntax/Executioner/index.html"}>>>]
  type = Transient
  dt = 0.1

  steady_state_detection = true
  steady_state_tolerance = 1e-08

  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-6

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  solve_type = NEWTON
  line_search = basic
[]

[Outputs<<<{"href": "../syntax/Outputs/index.html"}>>>]
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
  print_linear_residuals<<<{"description": "Enable printing of linear residuals to the screen (Console)"}>>> = false

  [console]
    type = Console<<<{"description": "Object for screen output.", "href": "../source/outputs/Console.html"}>>>
    outlier_variable_norms<<<{"description": "If true, outlier variable norms will be printed after each solve"}>>> = false
  []

  [csv]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../source/outputs/CSV.html"}>>>
    file_base<<<{"description": "The desired solution output name without an extension. If not provided, MOOSE sets it with Outputs/file_base when available. Otherwise, MOOSE uses input file name and this object name for a master input or uses master file_base, the subapp name and this object name for a subapp input to set it."}>>> = 'csv_thm/thm'
  []
[]