Lotus Griffin-Pronghorn Model

Mesh

Starting first with the mesh, this block defines the computational domain that the neutronics solve will operate on. Here, a cubit generated 3D mesh is imported as an exodus file. It is correspondingly manipulated to add in a reactor boundary condition, and is scaled to the correct size of the reactor.

[Mesh]
  coord_type = 'XYZ'
  [fmg]
    type = FileMeshGenerator
    file = '../mesh/mcre_mesh.e'
  []
  [reactor_out_bdry]
    type = ParsedGenerateSideset
    input = 'fmg'
    new_sideset_name = 'reactor_out'
    combinatorial_geometry = '(y > 775) | ((y <= 775) & (x > 1150 | x < -250))'
    include_only_external_sides = true
  []
  [scaling]
    type = TransformGenerator
    input = reactor_out_bdry
    transform = 'SCALE'
    vector_value = '0.001 0.001 0.001'
  []
  [delete_bad_sideset]
    type = BoundaryDeletionGenerator
    input = 'scaling'
    boundary_names = 'reactor_out'
  []

  [repair_mesh]
    type = MeshRepairGenerator
    input = delete_bad_sideset
    fix_node_overlap = true
  []
  [recreate_outer]
    type = SideSetsFromNormalsGenerator
    input = repair_mesh
    new_boundary = 'reactor_out'
    # Dummy parameter
    normals = '1 0 0'
    normal_tol = 1
  []
  [recreate_outer_part2]
    type = ParsedGenerateSideset
    input = 'recreate_outer'
    combinatorial_geometry = 'x<=-0.24999999'
    included_subdomains = 'reactor'
    new_sideset_name = 'reactor_out'
    include_only_external_sides = true
  []
[]

Transport Systems

Next, the process of converting the basic conservation equations into MOOSE variables and kernels is automated with the TransportSystems block:

[TransportSystems]
  particle = neutron
  equation_type = eigenvalue

  G = 9

  VacuumBoundary = 'wall-reflector reactor_out'

  [diff]
    scheme = CFEM-Diffusion
    n_delay_groups = 6
    external_dnp_variable = 'dnp'
    family = LAGRANGE
    order = FIRST
    fission_source_aux = true

    # For PJFNKMO
    assemble_scattering_jacobian = true
    assemble_fission_jacobian = true

    initialize_scalar_flux = true
  []
[]

Details about neutron transport equations can be found in Griffin theory manual.

Here we are specifying an eigenvalue neutronics problem using 9 energy groups (G = 9) solved via the diffusion approximation with a continuous finite element discretization scheme (scheme = CFEM-Diffusion).

For brevity, the multi-group diffusion equation for the eigenvalue problem solved at steady-state is given as:

(1)

where the symbols represent the following: diffusion coefficient for energy group , neutron scalar flux of energy group at position , removal cross-section from energy group , eigenvalue representing the effective multiplication factor of the reactor, prompt fission spectrum for neutrons born in energy group , number of neutrons per fission, average fission cross-section of neutrons in energy group , differential scattering cross-section for neutrons scattering from energy group to energy group , delayed neutrons fraction, delayed fission spectrum for neutrons born in energy group due to decay of neutron precursors, decay constant for precursor group , concentration of neutron precursor group , at position , and total number of energy groups.

Note the external_dnp_variable = 'dnp' parameter. This is a special option needed for liquid-fueled MSRs which signals that the conservation equations for delayed neutron precursors (DNPs) will be handled "externally" from the default Griffin implementation which assumes that the precursors do not have the turbulent treatment. DNP treatment will be discussed later in the thermal-hydraulics and DNP distribution input files.

The type of shape function Lagrange and order FIRST are defined here as well as the fission_source_aux which will be used by the sub-apps.

Power Density

Next, the Power Density is set for this specific reactor. Here the neutronics solution is normalized to the rated power of kW.

[PowerDensity]
  power = 25e3
  power_density_variable = power_density
  power_scaling_postprocessor = power_scaling
  family = MONOMIAL
  order = CONSTANT
[]

The power density is evaluated with the normalized neutronics solution. It provides the energy source in the fluid energy conservation equation. Because the fluid energy equation is discretized with FV, we evaluate the power density variable with a constant monomial function.

Auxiliary Variables

Moving on, the AuxVariables are now set which are specific variables that can be passed and read in from the other sub-apps for multiphysics calculations.

[AuxVariables]
  [tfuel]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 900.0 # in degree K
    block = 'reactor pump pipe'
  []
  [trefl]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 760.0 # in degree K
    block = 'reflector'
  []
  [c1]
    order = CONSTANT
    family = MONOMIAL
    block = 'reactor pump pipe'
  []
  [c2]
    order = CONSTANT
    family = MONOMIAL
    block = 'reactor pump pipe'
  []
  [c3]
    order = CONSTANT
    family = MONOMIAL
    block = 'reactor pump pipe'
  []
  [c4]
    order = CONSTANT
    family = MONOMIAL
    block = 'reactor pump pipe'
  []
  [c5]
    order = CONSTANT
    family = MONOMIAL
    block = 'reactor pump pipe'
  []
  [c6]
    order = CONSTANT
    family = MONOMIAL
    block = 'reactor pump pipe'
  []
  [dnp]
    order = CONSTANT
    family = MONOMIAL
    components = 6
    block = 'reactor pump pipe'
  []
[]

Here the tfuel and trefl are the temperature of the fuel salt and reflector, whereas c1 - c6 are the DNP distributions, and dnp is an array auxiliary with 6 component corresponding to the 6 DNP groups used here. This array is then read into Griffin in the external_dnp_variable = 'dnp' parameter within the Transport Systems block.

Auxiliary Kernels

Correspondingly, the AuxKernels are functors which operate on the AuxVariables. Here the six separate DNP groups are compiled into the Aux Variable dnp to be used within Griffin.

[AuxKernels]
  [build_dnp]
    type = BuildArrayVariableAux
    variable = dnp
    component_variables = 'c1 c2 c3 c4 c5 c6'
    execute_on = 'timestep_begin final'
  []
[]

Materials

The next critical need for any neutronics calculation is a set of multigroup cross sections. Generating cross sections is a topic that is left outside the scope of this example. A set has been generated for the LMCR problem and stored in the repository using Griffin's ISOXML format. These cross sections are included by the blocks,

[Materials]
  [reactor_material]
    type = CoupledFeedbackNeutronicsMaterial
    library_name = 'serpent_MCRE_xs'
    library_file = '../mgxs/serpent_MCRE_xs_new.xml'
    grid_names = 'tfuel'
    grid_variables = 'tfuel'
    plus = true
    isotopes = 'pseudo'
    densities = '1.0'
    is_meter = true
    material_id = 1
    block = 'reactor pump pipe'
  []
  [reflector_material]
    type = CoupledFeedbackNeutronicsMaterial
    library_name = 'serpent_MCRE_xs'
    library_file = '../mgxs/serpent_MCRE_xs_new.xml'
    grid_names = 'trefl'
    grid_variables = 'trefl'
    plus = true
    isotopes = 'pseudo'
    densities = '1.0'
    is_meter = true
    material_id = 2
    block = 'reflector'
  []
[]

CoupledFeedbackNeutronicsMaterial is able to use the temperature transferred from the fluid system for evaluating multigroup cross sections based on a table lookup on element quadrature points to bring in the feedback effect.

Executioner

Next, the Executioner block sets up the type of problem, and the numerical methods to solve the neutronics and multiphysics problems.

[Executioner]
  type = Eigenvalue
  solve_type = PJFNKMO

  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart '
  petsc_options_value = 'hypre boomeramg 50'
  nl_abs_tol = 1e-9
  #l_max_its = 100

  free_power_iterations = 5 # important to obtain fundamental mode eigenvalue
  normalization = fission_source_integral

  # Parameters for fixed point iteration with MultiApps
  fixed_point_algorithm = picard
  fixed_point_min_its = 2
  fixed_point_max_its = 30
  fixed_point_abs_tol = 1e-3
  relaxation_factor = 0.7
[]

The calculation is driven by Eigenvalue, an executioner available in Griffin. The PJFNKMO option for the solve_type parameter is able to drive the eigenvalue calculation with the contribution of DNP to the neutron transport equation as an external source scaled with -effective.

Additionally, PETSc options and tolerances for the neutronics and multiphysics fixed point iteration coupling method are provided.

Post Processors

The PostProcessors block sets up various calculations of reactor parameters that may be of interest to the user. This can be helpful to ensure the model is implemented correctly. Here the average, maximum, and minimum of various variables (e.g., power, fluxes, and DNPs) can be computed.

[Postprocessors]
  [power]
    type = ElementIntegralVariablePostprocessor
    variable = power_density
    execute_on = 'initial timestep_begin timestep_end transfer'
    block = 'reactor pump pipe'
  []
  [dt_limit]
    type = Receiver
    default = 0.1
  []

  [flux_0_average_reac]
    type = ElementAverageValue
    variable = sflux_g0
    block = 'reactor'
  []
  [flux_1_average_reac]
    type = ElementAverageValue
    variable = sflux_g1
    block = 'reactor'
  []
  [flux_2_average_reac]
    type = ElementAverageValue
    variable = sflux_g2
    block = 'reactor'
  []
  [flux_3_average_reac]
    type = ElementAverageValue
    variable = sflux_g3
    block = 'reactor'
  []
  [flux_4_average_reac]
    type = ElementAverageValue
    variable = sflux_g4
    block = 'reactor'
  []
  [flux_5_average_reac]
    type = ElementAverageValue
    variable = sflux_g5
    block = 'reactor'
  []
  [flux_6_average_reac]
    type = ElementAverageValue
    variable = sflux_g6
    block = 'reactor'
  []
  [flux_7_average_reac]
    type = ElementAverageValue
    variable = sflux_g7
    block = 'reactor'
  []
  [flux_0_average_ref]
    type = ElementAverageValue
    variable = sflux_g0
    block = 'reactor'
  []
  [flux_1_average_ref]
    type = ElementAverageValue
    variable = sflux_g1
    block = 'reflector'
  []
  [flux_2_average_ref]
    type = ElementAverageValue
    variable = sflux_g2
    block = 'reflector'
  []
  [flux_3_average_ref]
    type = ElementAverageValue
    variable = sflux_g3
    block = 'reflector'
  []
  [flux_4_average_ref]
    type = ElementAverageValue
    variable = sflux_g4
    block = 'reflector'
  []
  [flux_5_average_ref]
    type = ElementAverageValue
    variable = sflux_g5
    block = 'reflector'
  []
  [flux_6_average_ref]
    type = ElementAverageValue
    variable = sflux_g6
    block = 'reflector'
  []
  [flux_7_average_ref]
    type = ElementAverageValue
    variable = sflux_g7
    block = 'reflector'
  []
  [flux_0_average_reac_max]
    type = ElementExtremeValue
    variable = sflux_g0
    block = 'reactor'
  []
  [flux_1_average_reac_max]
    type = ElementExtremeValue
    variable = sflux_g1
    block = 'reactor'
  []
  [flux_2_average_reac_max]
    type = ElementExtremeValue
    variable = sflux_g2
    block = 'reactor'
  []
  [flux_3_average_reac_max]
    type = ElementExtremeValue
    variable = sflux_g3
    block = 'reactor'
  []
  [flux_4_average_reac_max]
    type = ElementExtremeValue
    variable = sflux_g4
    block = 'reactor'
  []
  [flux_5_average_reac_max]
    type = ElementExtremeValue
    variable = sflux_g5
    block = 'reactor'
  []
  [flux_6_average_reac_max]
    type = ElementExtremeValue
    variable = sflux_g6
    block = 'reactor'
  []
  [flux_7_average_reac_max]
    type = ElementExtremeValue
    variable = sflux_g7
    block = 'reactor'
  []
  [flux_0_average_ref_max]
    type = ElementExtremeValue
    variable = sflux_g0
    block = 'reactor'
  []
  [flux_1_average_ref_max]
    type = ElementExtremeValue
    variable = sflux_g1
    block = 'reflector'
  []
  [flux_2_average_ref_max]
    type = ElementExtremeValue
    variable = sflux_g2
    block = 'reflector'
  []
  [flux_3_average_ref_max]
    type = ElementExtremeValue
    variable = sflux_g3
    block = 'reflector'
  []
  [flux_4_average_ref_max]
    type = ElementExtremeValue
    variable = sflux_g4
    block = 'reflector'
  []
  [flux_5_average_ref_max]
    type = ElementExtremeValue
    variable = sflux_g5
    block = 'reflector'
  []
  [flux_6_average_ref_max]
    type = ElementExtremeValue
    variable = sflux_g6
    block = 'reflector'
  []
  [flux_7_average_ref_max]
    type = ElementExtremeValue
    variable = sflux_g7
    block = 'reflector'
  []
  [c1_average_reac]
    type = ElementAverageValue
    variable = c1
    block = 'reactor'
  []
  [c2_average_reac]
    type = ElementAverageValue
    variable = c2
    block = 'reactor'
  []
  [c3_average_reac]
    type = ElementAverageValue
    variable = c3
    block = 'reactor'
  []
  [c4_average_reac]
    type = ElementAverageValue
    variable = c4
    block = 'reactor'
  []
  [c5_average_reac]
    type = ElementAverageValue
    variable = c5
    block = 'reactor'
  []
  [c6_average_reac]
    type = ElementAverageValue
    variable = c6
    block = 'reactor'
  []
  [c1_average_pipe]
    type = ElementAverageValue
    variable = c1
    block = 'pipe pump'
  []
  [c2_average_pipe]
    type = ElementAverageValue
    variable = c2
    block = 'pipe pump'
  []
  [c3_average_pipe]
    type = ElementAverageValue
    variable = c3
    block = 'pipe pump'
  []
  [c4_average_pipe]
    type = ElementAverageValue
    variable = c4
    block = 'pipe pump'
  []
  [c5_average_pipe]
    type = ElementAverageValue
    variable = c5
    block = 'pipe pump'
  []
  [c6_average_pipe]
    type = ElementAverageValue
    variable = c6
    block = 'pipe pump'
  []
[]

Outputs

The Outputs block sets up the types of outputs the user would like to visualize or interpret the results. Here both an exodus and csv file are selected. Additionally, a restart file is generated so that transient solutions starting from steady-state can be computed.

[Outputs]
  exodus = true
  csv = true
  checkpoint = true
  [console]
    type = Console
    show = 'power dt_limit eigenvalue power_scaling'
  []
  [restart]
    type = Exodus
    execute_on = 'final'
    file_base = 'run_neutronics_restart'
  []
  print_linear_converged_reason = false
  print_linear_residuals = false
  print_nonlinear_converged_reason = false
[]

MultiApps

Finally, the MultiApps block sets up the sub-app that will be driven by the main Griffin app. Here, the Griffin input is the main-application and includes a single sub-application with the open source MOOSE Navier-Stokes input run_ns_initial.i.

[MultiApps]
  [ns]
    type = FullSolveMultiApp
    input_files = 'run_ns_initial.i'
    execute_on = 'timestep_end'
    # Not restoring means the solution keeps marching in time
    # with an updated power (heat source) distribution
    no_restore = true
    keep_solution_during_restore = true
  []
[]

Transfers

Correspondingly, the Transfers block sets up the Auxiliary Variables that will be passed to and from the thermal hydraulics sub app.

[Transfers]
  [power_density]
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = ns
    source_variable = power_density
    variable = power_density
  []
  [fission_source]
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = ns
    source_variable = fission_source
    variable = fission_source
  []

  [from_ns]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = ns
    source_variable = 'c1 c2 c3 c4 c5 c6 T T_ref'
    variable = 'c1 c2 c3 c4 c5 c6 tfuel trefl'
    source_type = 'variable_default variable_default variable_default variable_default variable_default
                   variable_default variable_default variable_default'
  []
[]

Here the power_density, and fission_sourceis transferred to the Navier-Stokes input file which operate as sources in the energy conservation and DNP advection equations respectively. Lastly, the DNP distributions for groups c1 - c6 and the temperature of the fuel T and the reflector T_ref are passed back to the neutronics solutution for multiphysics convergence.

Problem Parameters

The physical properties of the fuel salt and reflector (e.g., density, viscosity, thermal conductivity, and heat capacity) are defined first. Additionally, other parameters such as a friction force, pump force, porosity, and area are also set. Lastly, numerical interpolation schemes such as upwind and rhie-chow are also selected and a mixing length turbulence calibration informed from high-fidelity CFD is also set here.

################################################################################
## Molten Chloride Reactor - Lotus design                                     ##
## Pronghorn input file to initialize velocity fields                         ##
##                                                                            ##
## This input can be set to run a slow relaxation to steady state while       ##
## ramping down the fluid viscosity.                                          ##
################################################################################

# Notes:
# - These inputs are not using the Physics syntax to be consistent with inputs
#   from another project.
#   Please consider using the Physics syntax for setting up flow simulations for your own models
# - These inputs are deliberately not using later developments of the code.
#   Please contact the model owner if you require an updated version.

# Material properties fuel
rho = 3279. # density [kg / m^3]  (@1000K)
mu = 0.005926 # viscosity [Pa s]
k = 0.38 # [W / m / K]
cp = 640. # [J / kg / K]

# Material properties reflector
k_ref = 30. # [W / m / K]
cp_ref = 880. # [J / kg / K]
rho_ref = 3580. # [kg / m^3]

power = 25e3 # [W]

# alpha_b = '${fparse 1.0/rho}'

# Mass flow rate tuning
friction = 11.0 # [kg / m^4]
pump_force = '${fparse 0.25*4.0e6}' # [N / m^3]
porosity = 1.0
# T_hx = 592
Reactor_Area = '${fparse 3.14159*0.2575*0.2575}'
Pr = '${fparse mu*cp/k}'

# Numerical scheme parameters
advected_interp_method = 'upwind'
velocity_interp_method = 'rc'

# Calibration of the mixing length parameter
mixing_length_pipe_calibrated = '${fparse 0.07 * 0.1 * 0.06}'
mixing_length_reactor_calibrated = '${fparse 0.07 * 0.1 * 2}'

Global Parameters

Next, Global Parameters lists variables and parameters that will be used throughout the entire input file. Here the finite volume interpolation methods are set, and superficial velocities, pressure, porosity, density, viscosity, and the mixing_length model are defined for the entire input file. Some of these values are originally set in the Problem Parameters / header of the input and are referenced here using the ${} notation to be used implicitly for the rest of the input file.

[GlobalParams]
  rhie_chow_user_object = 'pins_rhie_chow_interpolator'

  two_term_boundary_expansion = true
  advected_interp_method = ${advected_interp_method}
  velocity_interp_method = ${velocity_interp_method}
  u = superficial_vel_x
  v = superficial_vel_y
  w = superficial_vel_z
  pressure = pressure
  porosity = porosity

  rho = rho
  mu = mu

  mixing_length = 'mixing_length'
[]

Mesh

This block defines the geometry of the thermal hydraulics computational domain, as shown in Figure 2. This block reads in the same 3D mesh that was generated using CUBIT which is identical for both Griffin and Pronghorn.

Then several boundaries and the mixing plate subdomain are defined by editing the mesh to correctly define the boundary conditions and porosity modeling necessary in the thermal hydraulics model.

[Mesh]
  coord_type = 'XYZ'
  [fmg]
    type = FileMeshGenerator
    file = '../mesh/mcre_mesh.e'
  []
  [scaling]
    type = TransformGenerator
    input = 'fmg'
    transform = 'SCALE'
    vector_value = '0.001 0.001 0.001'
  []
  [new_sideset_wall_reactor]
    type = SideSetsBetweenSubdomainsGenerator
    input = 'scaling'
    primary_block = 'reactor'
    paired_block = 'reflector'
    new_boundary = 'wall-reactor-full'
  []
  [reactor_top]
    type = ParsedGenerateSideset
    input = 'new_sideset_wall_reactor'
    new_sideset_name = 'reactor_top'
    included_subdomains = 'reactor'
    included_neighbors = 'pipe'
    combinatorial_geometry = 'x > 0'
  []
  [reactor_bottom]
    type = ParsedGenerateSideset
    input = 'reactor_top'
    new_sideset_name = 'reactor_bot'
    included_subdomains = 'reactor'
    included_neighbors = 'pipe'
    combinatorial_geometry = 'x < 0'
  []
  [mixing_plate]
    type = ParsedSubdomainMeshGenerator
    input = 'reactor_bottom'
    excluded_subdomains = 'pipe pump reflector'
    combinatorial_geometry = 'x > -0.3 & x < -0.22'
    block_id = '5'
    block_name = 'mixing-plate'
  []
  [delete_old_sidesets]
    type = BoundaryDeletionGenerator
    input = mixing_plate
    boundary_names = 'wall-reactor-caps  wall-reactor-reflector'
  []
  [mix_plate_downstream]
    type = ParsedGenerateSideset
    input = 'delete_old_sidesets'
    new_sideset_name = 'mixing-plate-downstream'
    included_subdomains = ' mixing-plate'
    combinatorial_geometry = ' (x > -0.23 & x < -0.2) & (y*y + z*z<0.12*0.12) '
  []
  [add-cooled-wall-pipes-outshield]
    type = ParsedGenerateSideset
    input = mix_plate_downstream
    combinatorial_geometry = 'y>1.'
    included_boundaries = wall-pipe
    new_sideset_name = 'heat-loss-section-outshield'
  []
  [add-cooled-wall-pipes-inshield]
    type = ParsedGenerateSideset
    input = add-cooled-wall-pipes-outshield
    combinatorial_geometry = 'y<1.'
    included_boundaries = wall-pipe
    new_sideset_name = 'heat-loss-section-inshield'
  []
  [heated-reflector_walls]
    type = ParsedGenerateSideset
    input = 'add-cooled-wall-pipes-inshield'
    new_sideset_name = 'heated-reflector-walls'
    included_subdomains = 'reflector'
    combinatorial_geometry = ' (x < 1.151 & y < -0.7749) | (x < 1.151 & y > 0.7749) |
                               (x < 1.151 & z < -0.5749) | (x < 1.151 & z > 0.5749)'
  []
  [non-heated-reflector_walls]
    type = ParsedGenerateSideset
    input = 'heated-reflector_walls'
    new_sideset_name = 'non-heated-reflector-walls'
    included_subdomains = 'reflector'
    combinatorial_geometry = ' (x > 1.149) | (x < -0.24)'
  []
  [reactor_out_bdry]
    type = ParsedGenerateSideset
    input = 'non-heated-reflector_walls'
    new_sideset_name = 'reactor_out'
    combinatorial_geometry = '(y > 775) | ((y <= 775) & (x > 1150 | x < -250))'
    include_only_external_sides = true
  []
  [regenerate_reactor_wall_to_reflector]
    type = SideSetsBetweenSubdomainsGenerator
    input = 'reactor_out_bdry'
    primary_block = 'reactor'
    paired_block = 'reflector'
    new_boundary = 'wall-reactor-reflector'
  []

  [repair_mesh]
    type = MeshRepairGenerator
    input = regenerate_reactor_wall_to_reflector
    fix_node_overlap = true
  []
[]

User Objects

The model now turns to implementing the equations, variables, kernels, and boundary conditions that need to be specified to solve the thermal hydraulic problem. Rather than using a Module where the Navier-Stokes Finite Volume action is used to define the problem, this model opts to explicitly set up the problem for added flexibility.

[UserObjects]
  [pins_rhie_chow_interpolator]
    type = PINSFVRhieChowInterpolator
    block = 'reactor pipe pump mixing-plate'
  []
[]

Here the porous formulation of the incompressible Rhie Chow Interpolator user object is specified and acts on the fluid blocks defined.

Variables

Next, the variables that must be explicitly solved for in the thermal hydraulics model are defined. Here the pressure, superficial velocities, and temperatures are required.

[Variables]
  [pressure]
    type = INSFVPressureVariable
    initial_condition = 0.1
    block = 'reactor pipe pump mixing-plate'
  []
  [superficial_vel_x]
    type = PINSFVSuperficialVelocityVariable
    initial_condition = 1e-6
    block = 'reactor pipe pump mixing-plate'
  []
  [superficial_vel_y]
    type = PINSFVSuperficialVelocityVariable
    initial_condition = 1e-6
    block = 'reactor pipe pump mixing-plate'
  []
  [superficial_vel_z]
    type = PINSFVSuperficialVelocityVariable
    initial_condition = 1e-6
    block = 'reactor pipe pump mixing-plate'
  []
  [lambda]
    family = SCALAR
    order = FIRST
    initial_condition = 1e-6
    block = 'reactor pipe pump mixing-plate'
  []
  [T]
    type = INSFVEnergyVariable
    initial_condition = 900.0
    block = 'reactor pipe pump mixing-plate'
    # Segregating T_fluid for performance
    # This is only OK because we are marching to steady state
    solver_sys = 'heat_transfer'
  []
  [T_ref]
    type = INSFVEnergyVariable
    initial_condition = 900.0
    block = 'reflector'
    scaling = 0.001
    # Segregating T_ref for performance
    solver_sys = 'heat_transfer'
  []
[]

Kernels

Correspondingly, the kernels are the functor or terms that manipulate the variables and form the conservation equations. Here the conservation of mass, momentum, and energy in the fuel salt and reflector are explicitly set.

[FVKernels]
  # The inactive parameter can be used to solve for a steady state
  # inactive = 'u_time v_time w_time heat_time heat_time_ref'

  [mass]
    type = PINSFVMassAdvection
    variable = pressure
    advected_interp_method = 'skewness-corrected'
    velocity_interp_method = 'rc'
    rho = ${rho}
  []
  [mean_zero_pressure]
    type = FVIntegralValueConstraint
    variable = pressure
    lambda = lambda
  []

  [u_time]
    type = PINSFVMomentumTimeDerivative
    variable = superficial_vel_x
    rho = ${rho}
    momentum_component = 'x'
  []
  [u_advection]
    type = PINSFVMomentumAdvection
    variable = superficial_vel_x
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = ${rho}
    momentum_component = 'x'
  []
  [u_viscosity]
    type = PINSFVMomentumDiffusion
    variable = superficial_vel_x
    momentum_component = 'x'
  []
  [u_viscosity_rans_pipe]
    type = INSFVMixingLengthReynoldsStress
    variable = superficial_vel_x
    rho = ${rho}
    momentum_component = 'x'
  []
  [u_pressure]
    type = PINSFVMomentumPressure
    variable = superficial_vel_x
    momentum_component = 'x'
    pressure = pressure
  []
  [u_friction_pump]
    type = PINSFVMomentumFriction
    variable = superficial_vel_x
    momentum_component = 'x'
    Darcy_name = 'DFC'
    block = 'pump'
  []
  [u_friction_pump_correction]
    type = PINSFVMomentumFrictionCorrection
    variable = superficial_vel_x
    momentum_component = 'x'
    rho = ${rho}
    Darcy_name = 'DFC'
    block = 'mixing-plate'
    consistent_scaling = 100.
  []

  [u_friction_mixing_plate]
    type = PINSFVMomentumFriction
    variable = superficial_vel_x
    momentum_component = 'x'
    Darcy_name = 'DFC_plate'
    block = 'mixing-plate'
  []
  [u_friction_mixing_plate_correction]
    type = PINSFVMomentumFrictionCorrection
    variable = superficial_vel_x
    momentum_component = 'x'
    rho = ${rho}
    Darcy_name = 'DFC_plate'
    block = 'mixing-plate'
    consistent_scaling = 100.
  []

  [v_time]
    type = PINSFVMomentumTimeDerivative
    variable = superficial_vel_y
    rho = ${rho}
    momentum_component = 'y'
  []
  [v_advection]
    type = PINSFVMomentumAdvection
    variable = superficial_vel_y
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = ${rho}
    momentum_component = 'y'
  []
  [v_viscosity]
    type = PINSFVMomentumDiffusion
    variable = superficial_vel_y
    momentum_component = 'y'
  []
  [v_viscosity_rans_pipe]
    type = INSFVMixingLengthReynoldsStress
    variable = superficial_vel_y
    rho = ${rho}
    mixing_length = ${mixing_length_pipe_calibrated}
    momentum_component = 'y'
    block = 'pipe pump'
  []
  [v_viscosity_rans]
    type = INSFVMixingLengthReynoldsStress
    variable = superficial_vel_y
    rho = ${rho}
    momentum_component = 'y'
  []
  [v_pressure]
    type = PINSFVMomentumPressure
    variable = superficial_vel_y
    momentum_component = 'y'
    pressure = pressure
  []
  [v_friction_pump]
    type = PINSFVMomentumFriction
    variable = superficial_vel_y
    momentum_component = 'y'
    Darcy_name = 'DFC'
    block = 'pump'
  []
  [v_friction_pump_correction]
    type = PINSFVMomentumFrictionCorrection
    variable = superficial_vel_y
    momentum_component = 'y'
    rho = ${rho}
    Darcy_name = 'DFC'
    block = 'mixing-plate'
    consistent_scaling = 100.
  []
  [v_friction_mixing_plate]
    type = PINSFVMomentumFriction
    variable = superficial_vel_y
    momentum_component = 'y'
    Darcy_name = 'DFC_plate'
    block = 'mixing-plate'
  []
  [pump]
    type = INSFVBodyForce
    variable = superficial_vel_y
    functor = ${pump_force}
    block = 'pump'
    momentum_component = 'y'
  []

  [w_time]
    type = PINSFVMomentumTimeDerivative
    variable = superficial_vel_z
    rho = ${rho}
    momentum_component = 'z'
  []
  [w_advection]
    type = PINSFVMomentumAdvection
    variable = superficial_vel_z
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = ${rho}
    momentum_component = 'z'
  []
  [w_viscosity]
    type = PINSFVMomentumDiffusion
    variable = superficial_vel_z
    momentum_component = 'z'
  []
  [w_viscosity_rans]
    type = INSFVMixingLengthReynoldsStress
    variable = superficial_vel_z
    rho = ${rho}
    momentum_component = 'z'
  []
  [w_pressure]
    type = PINSFVMomentumPressure
    variable = superficial_vel_z
    momentum_component = 'z'
    pressure = pressure
  []
  [w_friction_pump]
    type = PINSFVMomentumFriction
    variable = superficial_vel_z
    momentum_component = 'z'
    Darcy_name = 'DFC'
    block = 'pump'
  []
  [w_friction_mixing_plate]
    type = PINSFVMomentumFriction
    variable = superficial_vel_z
    momentum_component = 'z'
    Darcy_name = 'DFC_plate'
    block = 'mixing-plate'
  []

  ####### FUEL ENERGY EQUATION #######

  [heat_time]
    type = PINSFVEnergyTimeDerivative
    variable = T
    is_solid = false
    cp = ${cp}
  []
  [heat_advection]
    type = PINSFVEnergyAdvection
    variable = T
  []
  [heat_diffusion]
    type = PINSFVEnergyDiffusion
    variable = T
    k = ${k}
  []
  [heat_turbulent_diffusion]
    type = WCNSFVMixingLengthEnergyDiffusion
    variable = T
    schmidt_number = 0.9
    cp = ${cp}
  []
  [heat_src]
    type = FVBodyForce
    variable = T
    function = cosine_guess
    value = '${fparse power/0.21757}' #Volume integral of cosine shape is 0.21757
    block = 'reactor'
  []

  ####### REFLECTOR ENERGY EQUATION #######

  [heat_time_ref]
    type = PINSFVEnergyTimeDerivative
    variable = T_ref
    cp = ${cp_ref}
    rho = ${rho_ref}
    is_solid = true
    porosity = 0
  []
  [heat_diffusion_ref]
    type = FVDiffusion
    variable = T_ref
    coeff = ${k_ref}
  []
[]

The porous flow equations for incompressible flow read as follows:

(2)

(3)

(4)

(5)

where is the superficial velocity defined as and is the interstitial or physical velocity, is the density, is the pressure, is the internal energy, is the enthalpy, is the fluid temperature, is the solid temperature, is the solid density, is the specific heat of the solid phase, is the gravity vector, is the pressure drop coefficient, is the momentum source that is used to model the pump, is the effective thermal conductivity of the fluid, is the turbulent heat diffusivity, is the effective solid thermal conductivity, is the heat source being deposited directly in liquid fuel (e.g., fission heat source), and is the heat source in the solid (e.g., residual power production in structures). The effective thermal conductivities and are in general diagonal tensors.

The fluid domain comprises porous regions with and free flow regions with ; in the free-flow region, is not solved and we have , , and (where is the thermal conductivity of the fluid). Similarly, is not solved in solid-only regions where and with as the solid conductivity.

Finite Volume Interface Kernels

In addition to the Kernels that operate on the variables in the bulk domain, Interface Kernels operate at an interface. In this case the convection Interface Kernel is selected to account for heat transfer from the fuel salt to the reflector thorugh the vessel boundary.

[FVInterfaceKernels]
  [convection]
    type = FVConvectionCorrelationInterface
    variable1 = T
    variable2 = T_ref
    boundary = 'wall-reactor-reflector'
    h = htc
    T_solid = T_ref
    T_fluid = T
    subdomain1 = reactor
    subdomain2 = reflector
    wall_cell_is_bulk = true
  []
[]

Finite Volume Boundary Conditions

Next, the Finite Volume Boundary Conditions must be set for the conservation equations. Here no slip boundary conditions for the velocities for the momentum equation, and various heat transfer boundary conditions for the energy equation are specified.

[FVBCs]
  [no-slip-u]
    type = INSFVNoSlipWallBC
    boundary = 'wall-reactor wall-pipe wall-pump wall-reactor-full'
    variable = superficial_vel_x
    function = 0
  []
  [no-slip-v]
    type = INSFVNoSlipWallBC
    boundary = 'wall-reactor wall-pipe wall-pump wall-reactor-full'
    variable = superficial_vel_y
    function = 0
  []
  [no-slip-w]
    type = INSFVNoSlipWallBC
    boundary = 'wall-reactor wall-pipe wall-pump wall-reactor-full'
    variable = superficial_vel_z
    function = 0
  []
  [heat-losses-outshield]
    type = FVFunctorConvectiveHeatFluxBC
    variable = T
    T_bulk = T
    T_solid = 300.
    is_solid = false
    heat_transfer_coefficient = 3. #50.
    boundary = 'heat-loss-section-outshield'
  []
  [heated-outshield-pipe]
    type = FVFunctorNeumannBC
    variable = T
    boundary = 'heat-loss-section-outshield'
    functor = heat-input-pipe
  []
  [heat-losses-reflector]
    type = FVFunctorConvectiveHeatFluxBC
    variable = T_ref
    T_bulk = 350.
    T_solid = T_ref
    is_solid = true
    heat_transfer_coefficient = htc_rad_ref
    boundary = 'wall-reflector'
  []
  [heated-reflector-walls]
    type = FVFunctorNeumannBC
    variable = T_ref
    boundary = 'heated-reflector-walls'
    functor = heat-input-ref
  []
  [heated-inshield-pipe]
    type = FVNeumannBC
    variable = T
    boundary = 'heat-loss-section-inshield'
    value = 0.0 #500.
  []
[]

Auxiliary Variables

The Auxiliary Variables block defines variables of interest that can be operated on and transferred to sub-apps below or main-apps above. Here the power_density and fission_source is read in from the main-app neutronics solve, and c1-c6 is read in from the precursor advection sub-app. Additionally, auxiliary variables that will be transferred to the precursor advection sub-app are also defined here as a_u_var, a_v_var, and a_w_var respectively.

[AuxVariables]
  [h_DeltaT]
    type = MooseVariableFVReal
    block = 'reactor pipe pump mixing-plate'
  []
  [DeltaT_rad_aux]
    type = MooseVariableFVReal
    block = 'reflector'
  []
  [h_DeltaT_rad]
    type = MooseVariableFVReal
    block = 'reflector'
  []
  [a_u_var]
    type = MooseVariableFVReal
    block = 'reactor pipe pump mixing-plate'
  []
  [a_v_var]
    type = MooseVariableFVReal
    block = 'reactor pipe pump mixing-plate'
  []
  [a_w_var]
    type = MooseVariableFVReal
    block = 'reactor pipe pump mixing-plate'
  []
  [power_density]
    type = MooseVariableFVReal
  []
  [fission_source]
    type = MooseVariableFVReal
  []
  [c1]
    type = MooseVariableFVReal
    initial_condition = 1e-6
    block = 'reactor pipe pump mixing-plate'
  []
  [c2]
    type = MooseVariableFVReal
    initial_condition = 1e-6
    block = 'reactor pipe pump mixing-plate'
  []
  [c3]
    type = MooseVariableFVReal
    initial_condition = 1e-6
    block = 'reactor pipe pump mixing-plate'
  []
  [c4]
    type = MooseVariableFVReal
    initial_condition = 1e-6
    block = 'reactor pipe pump mixing-plate'
  []
  [c5]
    type = MooseVariableFVReal
    initial_condition = 1e-6
    block = 'reactor pipe pump mixing-plate'
  []
  [c6]
    type = MooseVariableFVReal
    initial_condition = 1e-6
    block = 'reactor pipe pump mixing-plate'
  []
[]

Auxiliary Kernels

Correspondingly, the Auxiliary Kernels operate on the Auxiliary Variables. In this case the auxiliary variables needed by the precursor advection sub-app are generated by copying the superficial velocities given by the thermal hydraulics solution.

[AuxKernels]
  [h_DeltaT_out]
    type = ParsedAux
    variable = h_DeltaT
    coupled_variables = 'T'
    expression = '3.*(T-300.)'
  []
  [DeltaT_rad_out]
    type = ParsedAux
    variable = DeltaT_rad_aux
    coupled_variables = 'T_ref'
    expression = 'T_ref-350.'
  []
  [h_DeltaT_rad_out]
    type = FunctorAux
    functor = 'DeltaT_rad_aux'
    variable = h_DeltaT_rad
    factor = 'htc_rad_ref'
  []
  [comp_a_u]
    type = FunctorAux
    functor = 'ax'
    variable = 'a_u_var'
    block = 'reactor pipe pump mixing-plate'
    execute_on = 'timestep_end'
  []
  [comp_a_v]
    type = FunctorAux
    functor = 'ay'
    variable = 'a_v_var'
    block = 'reactor pipe pump mixing-plate'
    execute_on = 'timestep_end'
  []
  [comp_a_w]
    type = FunctorAux
    functor = 'az'
    variable = 'a_w_var'
    block = 'reactor pipe pump mixing-plate'
    execute_on = 'timestep_end'
  []
[]

Functions

Next, various functions are defined to be used throughout the input file. Here heat transfer coefficients, reynolds numbers, power-density functions, and cosine guesses can all be generated and used for initial conditions, kernels, etc.

[Functions]
  [heat-input-ref]
    type = ParsedFunction
    expression = '0.'
  []
  [heat-input-pipe]
    type = ParsedFunction
    expression = '0.'
  []
  [Re_reactor]
    type = ParsedFunction
    expression = 'flow_hx_bot/Reactor_Area * (2*0.2575) / mu '
    symbol_names = 'flow_hx_bot mu Reactor_Area'
    #symbol_values = 'flow_hx_bot ${mu} ${Reactor_Area}'
    symbol_values = '25.2 ${mu} ${Reactor_Area}'
  []
  [htc]
    type = ParsedFunction
    expression = 'k/0.2575* 0.023 * Re_reactor^0.8 * Pr^0.3'
    symbol_names = 'k Re_reactor Pr'
    symbol_values = '${k} Re_reactor ${Pr}'
  []
  [htc_rad_ref]
    type = ParsedFunction
    # See https://web.mit.edu/16.unified/www/FALL/thermodynamics/notes/node136.html
    expression = '(T_ref_wall^2+350.^2)*(T_ref_wall+350.)*5.67e-8 / (1/0.18+1/0.35-1.)' #e_mgo=0.18 #e_steel=0.35
    symbol_names = 'T_ref_wall'
    symbol_values = 'T_ref_wall'
  []
  [power-density-func]
    type = ParsedFunction
    expression = '${power}/0.21757 * max(0, cos(x*pi/2/1.5))*max(0, cos(y*pi/2/1.5))*max(0, cos(z*pi/2/1.5))'
  []
  [ad_rampdown_mu_func]
    type = ParsedFunction
    # This expression can be set to impose a slowly decreasing profile in mu
    expression = 'mu' # *(0.1*exp(-3*t)+1)
    symbol_names = 'mu'
    symbol_values = ${mu}
  []
  [cosine_guess]
    type = ParsedFunction
    expression = 'max(0, cos(x*pi/2/1.5))*max(0, cos(y*pi/2/1.5))*max(0, cos(z*pi/2/1.5))'
  []
  [dts]
    type = PiecewiseLinear
    x = '0   1e7 ${fparse 1e7+0.05} ${fparse 1e7+0.2}'
    y = '0.1 0.1 1e7                1e7'
  []
[]

Materials

Next, the Materials block specifies material properties of interest in the fuel salt, reactor vessel, and reflector. Here, viscosity, wall frictions, and turbulent mixing length models can be implemented.

[FunctorMaterials]
  [generic]
    type = ADGenericFunctorMaterial
    prop_names = 'rho porosity'
    prop_values = '${rho} ${porosity}'
  []
  [interstitial_velocity_norm]
    type = PINSFVSpeedFunctorMaterial
    superficial_vel_x = superficial_vel_x
    superficial_vel_y = superficial_vel_y
    superficial_vel_z = superficial_vel_z
    porosity = porosity
  []
  [mu_spatial]
    type = ADPiecewiseByBlockFunctorMaterial
    prop_name = 'mu'
    subdomain_to_prop_value = 'pipe ad_rampdown_mu_func
                               pump ad_rampdown_mu_func
                               mixing-plate ad_rampdown_mu_func
                               reactor ad_rampdown_mu_func'
  []
  # TODO: remove this, fix ParsedFunctorMaterial
  [constant]
    type = GenericFunctorMaterial
    prop_names = 'constant_1'
    prop_values = '${fparse 2 * friction}'
  []
  [D_friction]
    type = ADParsedFunctorMaterial
    property_name = D_friction
    # This convoluted expression is because the model was developed back when
    # Pronghorn used a linear friction force as W * rho * v.
    # Now Pronghorn has explicit Darcy and Forchheimer coefficients
    expression = '2 * friction * rho / porosity / mu'
    functor_symbols = 'friction rho porosity mu'
    functor_names = 'constant_1 rho porosity mu'
  []
  [friction_material_pump]
    type = ADGenericVectorFunctorMaterial
    prop_names = 'DFC FFC'
    prop_values = 'D_friction D_friction D_friction
                   0 0 0'
  []
  [friction_material_mixing_plate]
    type = ADGenericVectorFunctorMaterial
    prop_names = 'DFC_plate FFC_plate'
    prop_values = 'D_friction D_friction D_friction
                   0 0 0'
  []
  [ins_fv]
    type = INSFVEnthalpyFunctorMaterial
    rho = ${rho}
    cp = ${cp}
    temperature = 'T'
  []
  [mixing_length_mat]
    type = ADPiecewiseByBlockFunctorMaterial
    prop_name = 'mixing_length'
    subdomain_to_prop_value = 'reactor      ${mixing_length_reactor_calibrated}
                               mixing-plate ${mixing_length_reactor_calibrated}
                               pipe         ${mixing_length_pipe_calibrated}
                               pump         ${mixing_length_pipe_calibrated}'
  []
[]

Executioner

The Executioner block defines how the problem will be run. Here the type of solution is a Transient and a TimeStepper function is defined. Additionally, the numerical tolerances for the solution are selected to control the level of convergence required.

[Executioner]
  type = Transient

  # Time-stepping parameters
  start_time = 0
  end_time = 1e7
  # dt = 1e6
  steady_state_detection = true
  steady_state_start_time = 5
  steady_state_tolerance = 1e-5

  [TimeStepper]
    type = IterationAdaptiveDT
    optimal_iterations = 10
    dt = 0.005
    timestep_limiting_postprocessor = 'dt_limit'
  []

  # Solver parameters
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -pc_factor_shift_type -ksp_gmres_restart -pc_factor_mat_solver_package'
  petsc_options_value = 'lu NONZERO 20 superlu_dist'
  line_search = 'none'
  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-4
  nl_max_its = 10 # fail early and try again with a shorter time step
  l_max_its = 80
  automatic_scaling = true

  # MultiApp iteration parameters
  relaxation_factor = 1.0

  # Multi-system parameter
  system_names = 'nl0 heat_transfer'
[]

Outputs

After computing the solution, the Outputs block defines how the user would like to process the final solution data. Here a CSV file is generated containing all of the calculations performed in the Postprocessors block.

[Outputs]
  [out]
    type = CSV
  []
  [restart]
    type = Exodus
    execute_on = 'timestep_end final'
  []
  # Reduce base output
  print_linear_converged_reason = false
  print_linear_residuals = false
  print_nonlinear_converged_reason = false
  hide = 'dt_limit'
[]

Postprocessors

The Postprocessors block is used to defined to calculate specific items of interest to the user. In this case, velocities, volumetric flow rate, temperatures, pressure drop, heat losses, and DNP concentrations can be calculated and printed in the Outputs CSV file.

[Postprocessors]
  [max_v]
    type = ElementExtremeValue
    variable = superficial_vel_x
    value_type = max
    block = 'reactor pump pipe'
  []
  [flow]
    type = VolumetricFlowRate
    boundary = 'reactor_bot'
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    vel_z = superficial_vel_z
    advected_quantity = ${rho}
  []
  [T_reac_bot]
    type = SideAverageValue
    boundary = 'reactor_bot'
    variable = T
  []
  [T_reac_top]
    type = SideAverageValue
    boundary = 'reactor_top'
    variable = T
  []
  [T_ref_wall]
    type = SideAverageValue
    boundary = 'wall-reflector'
    variable = T_ref
  []
  [pdrop]
    type = PressureDrop
    pressure = pressure
    upstream_boundary = 'reactor_bot'
    downstream_boundary = 'mixing-plate-downstream'
    boundary = 'mixing-plate-downstream reactor_bot'
  []
  [heat-loss-pipe]
    type = SideIntegralVariablePostprocessor
    boundary = 'heat-loss-section-outshield'
    variable = h_DeltaT
  []
  [heat-loss-reflector]
    type = SideIntegralVariablePostprocessor
    boundary = 'wall-reflector'
    variable = h_DeltaT_rad
  []
  [heat-input-ref-pp]
    type = FunctionSideIntegral
    boundary = 'heated-reflector-walls'
    function = heat-input-ref
  []
  [heat-input-pipe-pp]
    type = FunctionSideIntegral
    boundary = 'heat-loss-section-outshield'
    function = heat-input-pipe
  []
  [reactor-power]
    type = ElementIntegralFunctorPostprocessor
    block = 'reactor'
    functor = 'power-density-func'
  []
  [T_reac_ave]
    type = ElementAverageValue
    variable = T
    block = 'reactor'
  []
  [T_reac_max]
    type = ElementExtremeValue
    variable = T
    block = 'reactor'
  []
  [T_pipe_ave]
    type = ElementAverageValue
    variable = T
    block = 'pipe'
  []
  [T_ref_ave]
    type = ElementAverageValue
    variable = T_ref
    block = 'reflector'
  []
  [T_ref_max]
    type = ElementExtremeValue
    variable = T_ref
    block = 'reflector'
  []

  [c1_outlet]
    type = VolumetricFlowRate
    boundary = 'reactor_top'
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    vel_z = superficial_vel_z
    advected_quantity = c1
  []
  [c2_outlet]
    type = VolumetricFlowRate
    boundary = 'reactor_top'
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    vel_z = superficial_vel_z
    advected_quantity = c2
  []
  [c3_outlet]
    type = VolumetricFlowRate
    boundary = 'reactor_top'
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    vel_z = superficial_vel_z
    advected_quantity = c3
  []
  [c4_outlet]
    type = VolumetricFlowRate
    boundary = 'reactor_top'
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    vel_z = superficial_vel_z
    advected_quantity = c4
  []
  [c5_outlet]
    type = VolumetricFlowRate
    boundary = 'reactor_top'
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    vel_z = superficial_vel_z
    advected_quantity = c5
  []
  [c6_outlet]
    type = VolumetricFlowRate
    boundary = 'reactor_top'
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    vel_z = superficial_vel_z
    advected_quantity = c6
  []

  [c1_inlet]
    type = VolumetricFlowRate
    boundary = 'reactor_bot'
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    vel_z = superficial_vel_z
    advected_quantity = c1
  []
  [c2_inlet]
    type = VolumetricFlowRate
    boundary = 'reactor_bot'
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    vel_z = superficial_vel_z
    advected_quantity = c2
  []
  [c3_inlet]
    type = VolumetricFlowRate
    boundary = 'reactor_bot'
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    vel_z = superficial_vel_z
    advected_quantity = c3
  []
  [c4_inlet]
    type = VolumetricFlowRate
    boundary = 'reactor_bot'
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    vel_z = superficial_vel_z
    advected_quantity = c4
  []
  [c5_inlet]
    type = VolumetricFlowRate
    boundary = 'reactor_bot'
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    vel_z = superficial_vel_z
    advected_quantity = c5
  []
  [c6_inlet]
    type = VolumetricFlowRate
    boundary = 'reactor_bot'
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    vel_z = superficial_vel_z
    advected_quantity = c6
  []

  [dt_limit]
    type = Receiver
    default = 1e6
  []
[]

MultiApps

Moving next to the MultiApps block, this block connects the present app with any sub-apps. In this case the prec_transport sub-app is defined using the run_prec_transport.i input file.

[MultiApps]
  [prec_transport]
    type = TransientMultiApp
    input_files = 'run_prec_transport.i'
    execute_on = 'timestep_end'
    # no_restore = true
    # Allow smaller time steps by the child applications
    sub_cycling = true
  []
[]

Transfers

Lastly, the Transfers block is specified to showcase how Auxiliary Variables should be shared between the current application and the sub-app. Here power_density, fission_source, and superficial velocities, are sent to the precursor transport sub-app, whereas c1 - c6 will be returned back to the current app after solving the DNP distributions.

Then the neutronics application above calls for c1-c6 from the current application in it's Transfers block and uses these values to solve the multiphysics coupled neutronics eigenvalue problem.

[Transfers]
  [power_density]
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = prec_transport
    source_variable = power_density
    variable = power_density
  []
  [fission_source]
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = prec_transport
    source_variable = fission_source
    variable = fission_source
  []
  [u_x]
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = prec_transport
    source_variable = superficial_vel_x
    variable = superficial_vel_x
  []
  [u_y]
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = prec_transport
    source_variable = superficial_vel_y
    variable = superficial_vel_y
  []
  [u_z]
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = prec_transport
    source_variable = superficial_vel_z
    variable = superficial_vel_z
  []
  [a_u]
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = prec_transport
    source_variable = a_u_var
    variable = a_u_var
  []
  [a_v]
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = prec_transport
    source_variable = a_v_var
    variable = a_v_var
  []
  [a_w]
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = prec_transport
    source_variable = a_w_var
    variable = a_w_var
  []

  [c1]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = prec_transport
    source_variable = 'c1'
    variable = 'c1'
  []
  [c2]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = prec_transport
    source_variable = 'c2'
    variable = 'c2'
  []
  [c3]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = prec_transport
    source_variable = 'c3'
    variable = 'c3'
  []
  [c4]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = prec_transport
    source_variable = 'c4'
    variable = 'c4'
  []
  [c5]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = prec_transport
    source_variable = 'c5'
    variable = 'c5'
  []
  [c6]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = prec_transport
    source_variable = 'c6'
    variable = 'c6'
  []
[]

Problem Parameters

Similar to the thermal hydraulics input file run_ns_initial.i, the physical properties of the problem are defined first. Uniquely, the decay constants and production fractions are defined for the six DNP groups.

################################################################################
## Molten Chloride Reactor - Lotus design                                     ##
## Pronghorn input file to initialize DNP fields                              ##
################################################################################

# Mass flow rate tuning
porosity = 1.0
mu = 0.0166

# Numerical scheme parameters
advected_interp_method = 'upwind'
velocity_interp_method = 'rc'

# Dynamic scaling paramters
mixing_length_pipe_calibrated = '${fparse 0.07 * 0.1 * 0.06}'
mixing_length_reactor_calibrated = '${fparse 0.07 * 0.1 * 2}'
Sc_t = 0.9

# Delayed neutron precursor parameters. Lambda values are decay constants in
# [1 / s]. Beta values are production fractions [-].
lambda1 = 0.0124667
lambda2 = 0.0282917
lambda3 = 0.0425244
lambda4 = 0.133042
lambda5 = 0.292467
lambda6 = 0.666488
beta1 = 0.000250489
beta2 = 0.00104893
beta3 = 0.000726812
beta4 = 0.00143736
beta5 = 0.0022503
beta6 = 0.000680667

User Objects

Since this is a sub-app nested underneath the thermal hydraulics application, the rhie chow interpolator user object is updated with the superficial velocity aux variables a_u, a_v, and a_w transferred from the Thermal Hydraulics solve.

[UserObjects]
  [pins_rhie_chow_interpolator]
    type = PINSFVRhieChowInterpolator
    a_u = a_u_var
    a_v = a_v_var
    a_w = a_w_var
  []
[]

Variables

The next difference is that the c1 - c6 groups are no longer AuxVariables but are the finite volume variables that this problem is explicitly solving for.

[Variables]
  [c1]
    type = MooseVariableFVReal
  []
  [c2]
    type = MooseVariableFVReal
  []
  [c3]
    type = MooseVariableFVReal
  []
  [c4]
    type = MooseVariableFVReal
  []
  [c5]
    type = MooseVariableFVReal
  []
  [c6]
    type = MooseVariableFVReal
  []
[]

Auxiliary Variables

In order for this sub-app to take advantage of the transfers from applications above it, it must define the Auxiliary Variables that it will receive and send.

[AuxVariables]
  [power_density]
    type = MooseVariableFVReal
  []
  [fission_source]
    type = MooseVariableFVReal
  []
  [superficial_vel_x]
    type = PINSFVSuperficialVelocityVariable
  []
  [superficial_vel_y]
    type = PINSFVSuperficialVelocityVariable
  []
  [superficial_vel_z]
    type = PINSFVSuperficialVelocityVariable
  []
  [pressure]
    type = INSFVPressureVariable
  []
  [a_u_var]
    type = MooseVariableFVReal
  []
  [a_v_var]
    type = MooseVariableFVReal
  []
  [a_w_var]
    type = MooseVariableFVReal
  []
[]

Here the power_density, fission_source, superficial velocities, pressure, and auxiliary velocities are defined so they can be populated via transfers from the neutronics and thermal hydraulics applications above this sub-app.

Kernels

Correspondingly, the kernels are the functor or terms that manipulate the variables and form the conservation equations. Here the conservation of DNPs as seen in Eq. (6) is implemented by explicitly setting each term for each DNP group.

[FVKernels]
  [c1_advection]
    type = INSFVScalarFieldAdvection
    variable = c1
  []
  [c2_advection]
    type = INSFVScalarFieldAdvection
    variable = c2
  []
  [c3_advection]
    type = INSFVScalarFieldAdvection
    variable = c3
  []
  [c4_advection]
    type = INSFVScalarFieldAdvection
    variable = c4
  []
  [c5_advection]
    type = INSFVScalarFieldAdvection
    variable = c5
  []
  [c6_advection]
    type = INSFVScalarFieldAdvection
    variable = c6
  []
  [c1_diffusion]
    type = FVDiffusion
    coeff = '${fparse mu}'
    variable = c1
  []
  [c2_diffusion]
    type = FVDiffusion
    coeff = '${fparse mu}'
    variable = c2
  []
  [c3_diffusion]
    type = FVDiffusion
    coeff = '${fparse mu}'
    variable = c3
  []
  [c4_diffusion]
    type = FVDiffusion
    coeff = '${fparse mu}'
    variable = c4
  []
  [c5_diffusion]
    type = FVDiffusion
    coeff = '${fparse mu}'
    variable = c5
  []
  [c6_diffusion]
    type = FVDiffusion
    coeff = '${fparse mu}'
    variable = c6
  []
  [c1_diffusion_turb]
    type = INSFVMixingLengthScalarDiffusion
    variable = c1
  []
  [c2_diffusion_turb]
    type = INSFVMixingLengthScalarDiffusion
    variable = c2
  []
  [c3_diffusion_turb]
    type = INSFVMixingLengthScalarDiffusion
    variable = c3
  []
  [c4_diffusion_turb]
    type = INSFVMixingLengthScalarDiffusion
    variable = c4
  []
  [c5_diffusion_turb]
    type = INSFVMixingLengthScalarDiffusion
    variable = c5
  []
  [c6_diffusion_turb]
    type = INSFVMixingLengthScalarDiffusion
    variable = c6
  []
  [c1_src]
    type = FVCoupledForce
    variable = c1
    v = fission_source
    coef = ${beta1}
  []
  [c2_src]
    type = FVCoupledForce
    variable = c2
    v = fission_source
    coef = ${beta2}
  []
  [c3_src]
    type = FVCoupledForce
    variable = c3
    v = fission_source
    coef = ${beta3}
  []
  [c4_src]
    type = FVCoupledForce
    variable = c4
    v = fission_source
    coef = ${beta4}
  []
  [c5_src]
    type = FVCoupledForce
    variable = c5
    v = fission_source
    coef = ${beta5}
  []
  [c6_src]
    type = FVCoupledForce
    variable = c6
    v = fission_source
    coef = ${beta6}
  []
  [c1_decay]
    type = FVReaction
    variable = c1
    rate = ${lambda1}
  []
  [c2_decay]
    type = FVReaction
    variable = c2
    rate = ${lambda2}
  []
  [c3_decay]
    type = FVReaction
    variable = c3
    rate = ${lambda3}
  []
  [c4_decay]
    type = FVReaction
    variable = c4
    rate = ${lambda4}
  []
  [c5_decay]
    type = FVReaction
    variable = c5
    rate = ${lambda5}
  []
  [c6_decay]
    type = FVReaction
    variable = c6
    rate = ${lambda6}
  []
[]

It should be noted that the fission_source variable used as the source term for generating the DNPs was transfered from the Neutronics -> Thermal Hydraulics -> this precursor advection sub-app.

Finite Volume Boundary Conditions

Next, no-slip boundary conditions for the superficial velocities u, v, and w are defined at the boundaries.

[FVBCs]
  [no-slip-u]
    type = INSFVNoSlipWallBC
    boundary = 'wall-reactor wall-pipe wall-pump wall-reactor-full'
    variable = superficial_vel_x
    function = 0
  []
  [no-slip-v]
    type = INSFVNoSlipWallBC
    boundary = 'wall-reactor wall-pipe wall-pump wall-reactor-full'
    variable = superficial_vel_y
    function = 0
  []
  [no-slip-w]
    type = INSFVNoSlipWallBC
    boundary = 'wall-reactor wall-pipe wall-pump wall-reactor-full'
    variable = superficial_vel_z
    function = 0
  []
[]

Materials

Here, porosity and mixing length are incorporated as Materials into the model to correctly modify the DNP conservation equations depending on the fluid blocks porosity and turbulent mixing length.

[FunctorMaterials]
  [porous_mat]
    type = ADGenericFunctorMaterial
    prop_names = 'porosity'
    prop_values = '${porosity}'
    block = 'reactor pipe pump mixing-plate reflector'
  []
  [mixing_length_mat]
    type = ADPiecewiseByBlockFunctorMaterial
    prop_name = 'mixing_length'
    subdomain_to_prop_value = 'reactor      ${mixing_length_reactor_calibrated}
                               mixing-plate ${mixing_length_reactor_calibrated}
                               pipe         ${mixing_length_pipe_calibrated}
                               pump         ${mixing_length_pipe_calibrated}'
  []
[]