Molten Salt Reactor Experiment (MSRE) Multiphysics Model

Problem Parameters

The beginning of the input file lists the problem parameters the user can edit including the model's physical properties and initial conditions.

Global Parameters

Next, Global Parameters are defined in the Global Parameters block which give the macroscopic cross section file path for the MSRE.

[GlobalParams]
  library_file = '../mgxs/xs.xml'
  library_name = 'MSRE-Simplified'
  is_meter = true
[]

Transport Systems

Next, the Transport System block determines what sort of neutronic solve Griffin will perform. Here the methodology to compute the eigenvalue of the MSRE system is computed and corresponding boundary conditions are set — reflecting boundary for the RZ symmetric core and vacuum boundaries elsewhere.

Additionally, the type of transport equation solve is selected. In this case the a multigroup Diffusion equation is sufficient for fast multiphysics coupling. Additionally, options for the Jacobian and fission source auxiliary variable are selected here.

[TransportSystems]
  particle = neutron
  equation_type = eigenvalue
  G = 16
  ReflectingBoundary = 'left'
  VacuumBoundary = 'bottom loop_boundary right top top_core_barrel'
  [transport]
    scheme = CFEM-Diffusion
    family = LAGRANGE
    order = FIRST
    n_delay_groups = 6

    assemble_scattering_jacobian = true
    assemble_fission_jacobian = true

    external_dnp_variable = 'dnp'
    fission_source_aux = true
  []
[]

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

(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.

Mesh

The Mesh block can either generate a mesh using the MOOSE meshing capabiility or reads in a previously generated mesh. Here, the Mesh block reads in a MOOSE generated mesh and specifies that the model is RZ.

[Mesh]
  coord_type = 'RZ'
  [fmg]
    type = FileMeshGenerator
    file = '../mesh/mesh_in.e'
  []
[]

Aux Variables

The AuxVariables block is specificed next. Through the Auxiliary Variables, developers can define external variables that are solved or used in the primary simulation. The variable types, computation, and scaling factors are explained in detail here. Auxiliary variables can be set explicitly in the AuxKernel section, or set from another app using Transfers to perform multiphysics coupling.

In this case, the velocity, fuel salt and solid structures temperatures, and delayed neutron precursor group distributions will be informed by the Pronghorn sub app, whereas the volumetric power and fission rate fields will be calculated by Griffin and passed to the Pronghorn sub app.

[AuxVariables]
  [vel_x]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 0.000
    block = 'core lower_plenum upper_plenum down_comer riser pump elbow'
  []
  [vel_y]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 0.000
    block = 'core lower_plenum upper_plenum down_comer riser pump elbow'
  []
  [T_salt]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = ${T_Salt_initial}
    block = 'core lower_plenum upper_plenum down_comer riser pump elbow'
  []
  [T_solid]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = ${T_Salt_initial}
    block = 'core core_barrel'
  []
  [fission_rate]
    family = MONOMIAL
    order = CONSTANT
  []
  [Absorption_rate]
    family = MONOMIAL
    order = CONSTANT
  []
  [Leakage_rate]
    family = MONOMIAL
    order = CONSTANT
  []
  [c1]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 0.000
    block = 'core lower_plenum upper_plenum down_comer riser pump elbow'
  []
  [c2]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 0.000
    block = 'core lower_plenum upper_plenum down_comer riser pump elbow'
  []
  [c3]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 0.000
    block = 'core lower_plenum upper_plenum down_comer riser pump elbow'
  []
  [c4]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 0.000
    block = 'core lower_plenum upper_plenum down_comer riser pump elbow'
  []
  [c5]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 0.000
    block = 'core lower_plenum upper_plenum down_comer riser pump elbow'
  []
  [c6]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 0.000
    block = 'core lower_plenum upper_plenum down_comer riser pump elbow'
  []
  [dnp]
    order = CONSTANT
    family = MONOMIAL
    components = 6
    block = 'core lower_plenum upper_plenum down_comer riser pump elbow'
  []

  [ad_U235]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 9.692763E-05
  []
  [ad_U238]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 1.967924E-04
  []
  [ad_Be9]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 9.496943E-03
  []
  [ad_Li7]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 2.121311E-02
  []
  [ad_F9]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 4.790899E-02
  []
  [ad_Zr90]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 8.395505E-04
  []
  [ad_Zr91]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 1.830855E-04
  []
  [ad_Zr92]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 2.798495E-04
  []
  [ad_Zr94]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 2.836027E-04
  []
  [ad_Zr96]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 4.568967E-05
  []
  [ad_C12]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = 7.247622E-02
    block = 'core core_barrel'
  []
[]

Aux Kernels

Correspondingly, the AuxKernels block specifies Kernels which act on the auxiliary variables. Here an array of delayed neutron precursors concentration is built to be read in by Griffin. Additionally, the fuel salt nuclide constituent concentrations are updated due to fuel salt density changes.

[AuxKernels]
  [build_dnp]
    type = BuildArrayVariableAux
    variable = dnp
    component_variables = 'c1 c2 c3 c4 c5 c6'
    execute_on = 'initial timestep_begin final'
  []
  [update_ad_f_U235]
    type = ParsedAux
    block = 'lower_plenum upper_plenum down_comer riser pump elbow'
    variable = ad_U235
    coupled_variables = 'T_salt'
    expression = '9.692763E-05*(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_f_U238]
    type = ParsedAux
    block = 'lower_plenum upper_plenum down_comer riser pump elbow'
    variable = ad_U238
    coupled_variables = 'T_salt'
    expression = '1.967924E-04 *(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_f_Be9]
    type = ParsedAux
    block = 'lower_plenum upper_plenum down_comer riser pump elbow'
    variable = ad_Be9
    coupled_variables = 'T_salt'
    expression = '9.496943E-03*(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_f_Li7]
    type = ParsedAux
    block = 'lower_plenum upper_plenum down_comer riser pump elbow'
    variable = ad_Li7
    coupled_variables = 'T_salt'
    expression = '2.121311E-02 *(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_f_F9]
    type = ParsedAux
    block = 'lower_plenum upper_plenum down_comer riser pump elbow'
    variable = ad_F9
    coupled_variables = 'T_salt'
    expression = '4.790899E-02*(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_f_Zr90]
    type = ParsedAux
    block = 'lower_plenum upper_plenum down_comer riser pump elbow'
    variable = ad_Zr90
    coupled_variables = 'T_salt'
    expression = '8.395505E-04*(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_f_Zr91]
    type = ParsedAux
    block = 'lower_plenum upper_plenum down_comer riser pump elbow'
    variable = ad_Zr91
    coupled_variables = 'T_salt'
    expression = '1.830855E-04*(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_f_Zr92]
    type = ParsedAux
    block = 'lower_plenum upper_plenum down_comer riser pump elbow'
    variable = ad_Zr92
    coupled_variables = 'T_salt'
    expression = '2.798495E-04*(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_f_Zr94]
    type = ParsedAux
    block = 'lower_plenum upper_plenum down_comer riser pump elbow'
    variable = ad_Zr94
    coupled_variables = 'T_salt'
    expression = '2.836027E-04*(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_f_Zr96]
    type = ParsedAux
    block = 'lower_plenum upper_plenum down_comer riser pump elbow'
    variable = ad_Zr96
    coupled_variables = 'T_salt'
    expression = '4.568967E-05*(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_c_U235]
    type = ParsedAux
    block = 'core'
    variable = ad_U235
    coupled_variables = 'T_salt'
    expression = '2.159856E-05*(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_c_U238]
    type = ParsedAux
    block = 'core'
    variable = ad_U238
    coupled_variables = 'T_salt'
    expression = '4.385162E-05*(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_c_Be9]
    type = ParsedAux
    block = 'core'
    variable = ad_Be9
    coupled_variables = 'T_salt'
    expression = '2.116221E-03*(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_c_Li7]
    type = ParsedAux
    block = 'core'
    variable = ad_Li7
    coupled_variables = 'T_salt'
    expression = '4.726958E-03*(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_c_F9]
    type = ParsedAux
    block = 'core'
    variable = ad_F9
    coupled_variables = 'T_salt'
    expression = '1.067565E-02*(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_c_Zr90]
    type = ParsedAux
    block = 'core'
    variable = ad_Zr90
    coupled_variables = 'T_salt'
    expression = '1.870786E-04*(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_c_Zr91]
    type = ParsedAux
    block = 'core'
    variable = ad_Zr91
    coupled_variables = 'T_salt'
    expression = '4.079728E-05*(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_c_Zr92]
    type = ParsedAux
    block = 'core'
    variable = ad_Zr92
    coupled_variables = 'T_salt'
    expression = '6.235937E-05*(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_c_Zr94]
    type = ParsedAux
    block = 'core'
    variable = ad_Zr94
    coupled_variables = 'T_salt'
    expression = '6.319571E-05*(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
  [update_ad_Zr96]
    type = ParsedAux
    block = 'core'
    variable = ad_Zr96
    coupled_variables = 'T_salt'
    expression = '1.018111E-05*(1.0-0.4798*(T_salt-${T_Salt_initial})/${Salt_Density_initial})'
    execute_on = 'INITIAL timestep_end'
  []
[]

User Objects

Next, the UserOjbects block specifies User Objects that can be used by other MOOSE applications. Here both the transport solutions and aux variable solutions are stored in their corresponding user objects. This user object is necessary for restarting a transient simulation from a steady-state solution.

[UserObjects]
  [transport_solution_s1]
    type = TransportSolutionVectorFile
    transport_system = transport
    writing = true
    execute_on = 'FINAL'
  []
  [auxvar_solution_s1]
    type = SolutionVectorFile
    var = 'vel_x  vel_y  T_salt  T_solid  c1  c2  c3  c4  c5  c6
           ad_C12  ad_U235 ad_U238 ad_Be9  ad_Li7  ad_F9
           ad_Zr90 ad_Zr91 ad_Zr92 ad_Zr94 ad_Zr96'
    writing = true
    execute_on = 'FINAL'
  []
[]

Power Density

The PowerDensity block specifies how the Power Density is calculated in Griffin. Here the total power is used to scale the power density calculated during the eigenvalue solution.

[PowerDensity]
  power = '${fparse total_power}'
  power_density_variable = power_density
  integrated_power_postprocessor = total_power
  family = L2_LAGRANGE
  order = FIRST
[]

Materials

The Materials block specifies what materials will be used in the neutronic calculation. Each material links to the corresponding material IDs from the macroscopic ISOXML cross section library.

Here the nuclide densities are updated through the Aux Kernel operations, and transfered to the Materials block via the auxvar_solution_s1 user object.

[Materials]
  [activeCore]
    type = CoupledFeedbackNeutronicsMaterial
    grid_names = ' Tfuel'
    grid_variables = ' T_salt'
    plus = true
    isotopes = ' C12   U235   U238    BE9      LI7    F19
                ZR90   ZR91   ZR92   ZR94   ZR96'
    densities = 'ad_C12   ad_U235  ad_U238  ad_Be9   ad_Li7   ad_F9
                 ad_Zr90  ad_Zr91  ad_Zr92  ad_Zr94  ad_Zr96'
    material_id = 1
    block = 'core'
  []
  [flow_loop]
    type = CoupledFeedbackNeutronicsMaterial
    grid_names = ' Tfuel '
    grid_variables = ' T_salt'
    plus = true
    isotopes = 'U235   U238    BE9      LI7    F19
                ZR90   ZR91   ZR92   ZR94   ZR96'
    densities = 'ad_U235  ad_U238  ad_Be9   ad_Li7   ad_F9
                 ad_Zr90  ad_Zr91  ad_Zr92  ad_Zr94  ad_Zr96'
    material_id = 2
    block = 'lower_plenum upper_plenum down_comer riser pump elbow'
  []
  [core_barrel]
    type = CoupledFeedbackNeutronicsMaterial
    grid_names = ' Tfuel '
    grid_variables = ' T_solid'
    plus = true
    isotopes = 'C12'
    densities = 'ad_C12'
    material_id = 1
    block = 'core_barrel'
  []
[]

Postprocessors

The Postprocessors block allows for the user to calculate additional items of interest and have them included in the output.

Here the max and average temperatures for both the fuel and moderator, as well as reaction rates and neutron leakage are calculated in this block.

[Postprocessors]
  [Fuel_max_Temp]
    type = ElementExtremeValue
    variable = T_salt
    block = 'core lower_plenum upper_plenum down_comer riser pump elbow'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Fuel_avg_Temp]
    type = ElementAverageValue
    variable = T_salt
    block = 'core lower_plenum upper_plenum down_comer riser pump elbow'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Fuel_core_max_Temp]
    type = ElementExtremeValue
    variable = T_salt
    block = 'core'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Fuel_core_avg_Temp]
    type = ElementAverageValue
    variable = T_salt
    block = 'core'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Mod_max_Temp]
    type = ElementExtremeValue
    variable = T_solid
    block = 'core core_barrel'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Mod_avg_Temp]
    type = ElementAverageValue
    variable = T_solid
    block = 'core core_barrel'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [nufission_RR]
    type = FluxRxnIntegral
    block = 'core lower_plenum upper_plenum down_comer riser pump elbow'
    cross_section = nu_sigma_fission
    coupled_flux_groups = ' sflux_g0    sflux_g1    sflux_g2    sflux_g3 sflux_g4    sflux_g5    sflux_g6    sflux_g7 sflux_g8    sflux_g9    sflux_g10   sflux_g11 sflux_g12   sflux_g13   sflux_g14   sflux_g15'
    execute_on = 'transfer timestep_end'
  []
  [absorption_RR]
    type = FluxRxnIntegral
    block = 'core lower_plenum upper_plenum down_comer riser pump elbow'
    cross_section = sigma_absorption
    coupled_flux_groups = ' sflux_g0    sflux_g1    sflux_g2    sflux_g3 sflux_g4    sflux_g5    sflux_g6    sflux_g7 sflux_g8    sflux_g9    sflux_g10   sflux_g11 sflux_g12   sflux_g13   sflux_g14   sflux_g15'
    execute_on = 'transfer timestep_end'
  []
  [fission_RR]
    type = FluxRxnIntegral
    block = 'core lower_plenum upper_plenum down_comer riser pump elbow'
    cross_section = sigma_fission
    coupled_flux_groups = ' sflux_g0    sflux_g1    sflux_g2    sflux_g3 sflux_g4    sflux_g5    sflux_g6    sflux_g7 sflux_g8    sflux_g9    sflux_g10   sflux_g11 sflux_g12   sflux_g13   sflux_g14   sflux_g15'
    execute_on = 'transfer timestep_end'
  []
  [ngamma_RR]
    type = FluxRxnIntegral
    block = 'core lower_plenum upper_plenum down_comer riser pump elbow'
    cross_section = sigma_ngamma
    coupled_flux_groups = ' sflux_g0    sflux_g1    sflux_g2    sflux_g3 sflux_g4    sflux_g5    sflux_g6    sflux_g7 sflux_g8    sflux_g9    sflux_g10   sflux_g11 sflux_g12   sflux_g13   sflux_g14   sflux_g15'
    execute_on = 'transfer timestep_end'
  []
  [Leakage]
    type = PartialSurfaceCurrent
    boundary = 'right bottom loop_boundary pump_outlet downcomer_inlet top_core_barrel'
    transport_system = transport
  []
[]

Executioner

The Executioner block sets up how the neutronics solve is performed. Here the eigenvalue solve is selected, and the numerical method of solving it — the Preconditioned Jacobian Free Newton Krylov - Matrix Only (PJFNKMO) — is specified.

Additionally, PETSc options and tolerances for internal convergence and Picard iterations for multiphysics convergence are selected.

[Executioner]
  type = Eigenvalue
  solve_type = PJFNKMO

  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart'
  petsc_options_value = 'hypre boomeramg 50'
  l_max_its = 5000
  nl_abs_tol = 1e-10

  free_power_iterations = 4 # important to obtain fundamental mode eigenvalue

  # MultiApp
  fixed_point_min_its = 2
  fixed_point_max_its = 50
  fixed_point_rel_tol = 1e-6
  fixed_point_abs_tol = 1e-6
[]

Multi Apps

Finally, the MultiApps block sets up any sub apps that will be used in a multiphysics coupling. Here we just have one multi app which is the Pronghorn solve of the thermal hydraulics and delayed neutron precursor group distributions.

[MultiApps]
  [flow_dnp]
    type = FullSolveMultiApp
    input_files = 'th.i'
    execute_on = 'timestep_end'
    max_procs_per_app = 48
    keep_solution_during_restore = true
  []
[]

Transfers

Correspondingly, the Transfers block determines how the main and sub apps will communicate with each other. These transfers send auxiliary variables to and from different apps. Here the power_density and fission_source calculated by Griffin is sent to the sub-app Pronghorn to inform the thermal hydraulic solution.

In return, Pronghorn sends the delayed neutron precursor group distributions c1-c6, the solid and liquid temperatures, and the velocities to Griffin to determine the steady state multiphysics solution of the reactor.

[Transfers]
  [power_density]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    to_multi_app = flow_dnp
    source_variable = power_density
    variable = power_density
    execute_on = 'timestep_end'
  []
  [fission_source]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    to_multi_app = flow_dnp
    source_variable = fission_source
    variable = fission_source
    execute_on = 'timestep_end'
  []
  [c1]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    from_multi_app = flow_dnp
    source_variable = 'c1'
    variable = 'c1'
    execute_on = 'timestep_end'
  []
  [c2]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    from_multi_app = flow_dnp
    source_variable = 'c2'
    variable = 'c2'
    execute_on = 'timestep_end'
  []
  [c3]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    from_multi_app = flow_dnp
    source_variable = 'c3'
    variable = 'c3'
    execute_on = 'timestep_end'
  []
  [c4]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    from_multi_app = flow_dnp
    source_variable = 'c4'
    variable = 'c4'
    execute_on = 'timestep_end'
  []
  [c5]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    from_multi_app = flow_dnp
    source_variable = 'c5'
    variable = 'c5'
    execute_on = 'timestep_end'
  []
  [c6]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    from_multi_app = flow_dnp
    source_variable = 'c6'
    variable = 'c6'
    execute_on = 'timestep_end'
  []
  [T_salt]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    from_multi_app = flow_dnp
    source_variable = 'T_fluid'
    variable = 'T_salt'
    execute_on = 'timestep_end'
  []
  [T_graph]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    from_multi_app = flow_dnp
    source_variable = 'T_solid'
    variable = 'T_solid'
    execute_on = 'timestep_end'
  []
  [vel_x]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    from_multi_app = flow_dnp
    source_variable = 'superficial_vel_x'
    variable = 'vel_x'
    execute_on = 'timestep_end'
  []
  [vel_y]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    from_multi_app = flow_dnp
    source_variable = 'superficial_vel_y'
    variable = 'vel_y'
    execute_on = 'timestep_end'
  []
[]

Outputs

Lastly, the Outputs block specifies what type of output (e.g. exodus and CSV) will be created during the simulation.

[Outputs]
  file_base = msre_neutronics_ss_s2_out
  exodus = true
  csv = true
  perf_graph = true
  execute_on = 'INITIAL FINAL TIMESTEP_END'
[]

Problem Parameters

The physical properties are defined first. The core porosity is defined at a ratio of 0.222831, calculated as the quotient of flow area by total core area. The porosity for the rest of the components is set to 1, full fluid region, but is editable by the user. The thermophysical properties of the solid structures are also defined here.

# Properties -------------------------------------------------------------------
core_porosity = 0.222831853 # core porosity salt VF=0.222831853, Graphite VF=0.777168147
down_comer_porosity = 1.0 # downcomer porosity
lower_plenum_porosity = 0.5 # lower pelnum porosity
upper_plenum_porosity = 1.0 # upper pelnum porosity
riser_porosity = 1.0 # riser porosity
pump_porosity = 1.0 # pump porosity
elbow_porosity = 1.0 # elbow porosity

cp_steel = 500.0 # (J/(kg.K)) specific heat of steel
rho_steel = 8000.0 # (kg/(m3)) density of steel
k_steel = 15.0 # # (W/(m.k)) density of steel

The following section focuses on the operational parameters. The mass flow rate was defined at 191.19 kilograms per second, with a core outlet pressure approximately atmospheric at 101.325 kiloPascals. The salt core inlet temperature is defined to be 908.15 Kelvin (K) with an ambient room temperature defined to be 300K. Finally, this section defines the pump force scaling for the centrifugal pump functor to obtain a specific pump force of 1800.0 kiloNewton per meter cube (kN/m). The pump force was adapted to produce a loop circulation time of ~25 seconds.

# Operational Parameters --------------------------------------------------------
#p_outlet = 1.01325e+05 # Reactor outlet pressure (Pa)
p_outlet = 1.50653e+05 # Reactor outlet pressure (Pa)
T_inlet = 908.15 # Salt inlet temperature (K).
T_Salt_initial = 923.0 # inital salt temperature (will change in steady-state)

pump_force = -1.3e6 # pump force functor (set to get a loop circulation time of ~25 seconds)
vol_hx = 1e10 # (W/(m3.K)) volumetric heat exchange coefficient for heat exchanger
# Note: vol_hx need to be tuned to match intermediate HX performance for transients
bulk_hx = 100.0 # (W/(m3.K)) core bulk volumetric heat exchange coefficient (already callibrated)

Next, the hydraulic diameters of the flow channels are set respectively. The fluid blocks are defined to indicate the fluid regions.

# Thermal-Hydraulic diameters ----------------------------------------------------
D_H_fuel_channel = 0.0191334114 # Hydraulic diameter of bypass
D_H_downcomer = 0.045589414 # Hydraulic diameter of riser
D_H_pipe = '${fparse 5*0.0254}' # Riser Hydraulic Diameter
D_H_plena = '${fparse 2*core_radius}' # Hydraulic diameter of riser

Then, first section defines the delayed neutron group properties, both their production rate due to fission and decay rate. It also define the turbulent Schmidt number and the initial temperature of the fuel salt.

# Delayed neutron precursors constants ------------------------------------------
lambda1 = 0.0133104
lambda2 = 0.0305427
lambda3 = 0.115179
lambda4 = 0.301152
lambda5 = 0.879376
lambda6 = 2.91303
beta1 = 8.42817e-05
beta2 = 0.000684616
beta3 = 0.000479796
beta4 = 0.00103883
beta5 = 0.000549185
beta6 = 0.000184087

Sc_t = 1 # turbulent Schmidt number

The end of the first section defines the set of blocks where the fluid flow and solid heat conduction equations will be solved.

# Utils -------------------------------------------------------------------------
# fluid blocks define fluid vars and solve for them
fluid_blocks = 'core lower_plenum upper_plenum down_comer riser pump elbow'
solid_blocks = 'core core_barrel'

Global Parameters

Next, Global Parameters are defined in the Global Parameters block. This section is used to define variables and methods that are required by many methods in the input file.

[GlobalParams]
  fp = fluid_properties_obj
  porosity = 'porosity'
  rhie_chow_user_object = 'pins_rhie_chow_interpolator'

  u = superficial_vel_x
  v = superficial_vel_y

  advected_interp_method = 'upwind'
  velocity_interp_method = 'rc'
  mixing_length = 'mixing_length'
[]

Mesh

This block defines the geometry of the core, as shown in Figure 1. This block reads in the same mesh that was generated using MOOSE which is identical for both Griffin and Pronghorn.

[Mesh]
  coord_type = 'RZ'
  [fmg]
    type = FileMeshGenerator
    file = '../mesh/mesh_in.e'
  []
[]

Modules

This block defines the Module that will be used in this problem - specifically the Navier Stokes Finite Volume action.

The Navier Stokes Finite Volume action is used to define fluid properties, compute the "weakly-compressible" flow, and set up boundary conditions.

Here, the action operates on the fluid blocks which are defined at the top of the header and are substituted using the $ sign. Recall that the $ sign refers to variable substitutions. "Weakly-compressible" with porous_medium_treatment is selected for this model and the energy equation is also included.

Numerical schemes can be selected for solving the Navier Stokes equations with different methods of pressure and velocity interpolation.

Porosity, friction, and turbulence treatment can be modified here using various correlations and models. Fluid properties are also set here and can be modified.

The action also incorporates external physics like a volumetric heat source due to the power_density calculated by Griffin. Also, boundary conditions are set to constrain the different Navier Stokes equations.

The scalar equations for delayed neutron precursor groups are not included here in the action since they are defined externally elsewhere in the input file

Lastly, scaling factors can be tuned to increase convergence speed if needed.

The action allows users to simply modify the equations that are solved, choose numerical schemes, define the porosity and friction treatment, define fluid properties, couple with other physics for energy deposition, and set boundary conditions.

[Modules]
  [NavierStokesFV]
    # Basic settings - weakly-compressible, turbulent flow with buoyancy
    block = ${fluid_blocks}
    compressibility = 'weakly-compressible'
    porous_medium_treatment = true
    add_energy_equation = true
    gravity = '0.0 -9.81 0.0'

    # Variable naming
    velocity_variable = 'superficial_vel_x superficial_vel_y'
    pressure_variable = 'pressure'
    fluid_temperature_variable = 'T_fluid'

    # Numerical schemes
    momentum_advection_interpolation = upwind
    mass_advection_interpolation = upwind
    energy_advection_interpolation = upwind
    velocity_interpolation = rc

    # Porous & Friction treatement
    use_friction_correction = true
    friction_types = 'darcy forchheimer'
    friction_coeffs = 'Darcy_coefficient Forchheimer_coefficient'
    consistent_scaling = 100.0
    porosity_smoothing_layers = 2
    turbulence_handling = 'mixing-length'

    # fluid properties
    density = 'rho'
    dynamic_viscosity = 'mu'
    thermal_conductivity = 'kappa'
    specific_heat = 'cp'

    # Energy source-sink
    external_heat_source = 'power_density'
    #energy_scaling = 2.0

    # Boundary Conditions
    wall_boundaries = 'left      top      bottom   right    loop_boundary '
    momentum_wall_types = 'symmetry  slip     noslip   noslip   noslip'
    energy_wall_types = 'heatflux  heatflux heatflux heatflux heatflux'
    energy_wall_function = '0        0        0        0        0'

    # Constrain Pressure
    pin_pressure = true
    pinned_pressure_value = ${p_outlet}
    pinned_pressure_point = '0.0 2.13859 0.0'
    pinned_pressure_type = point-value-uo

    # Passive Scalar -- solved separetely to integrate porosity jumps
    add_scalar_equation = false
  []
[]

The porous flow equations for weakly-compressible 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.

Fluid Properties

This block contains parameters applied to all the core components. The fluid_properties_obj refers to the primary salt F-Li-Be (FLiBe), already defined within MOOSE.

[FluidProperties]
  [fluid_properties_obj]
    type = FlibeFluidProperties
  []
[]

Fluid properties are further defined within the Navier Stokes Finite Volume action and the Materials block.

    # fluid properties
    density = 'rho'
    dynamic_viscosity = 'mu'
    thermal_conductivity = 'kappa'
    specific_heat = 'cp'

  # Setting up fluid properties at blocks material blocks
  [fluid_props_to_mat_props]
    type = GeneralFunctorFluidProps
    pressure = 'pressure'
    T_fluid = 'T_fluid'
    speed = 'speed'
    characteristic_length = characteristic_length
    block = ${fluid_blocks}
  []

Variables

This block contains the variables that are explicitly solved for in this model. These include the velocity, pressure, temperature, and six delayed neutron precursor groups.

[Variables]
  [superficial_vel_x]
    type = PINSFVSuperficialVelocityVariable
    initial_condition = 1e-8
    block = ${fluid_blocks}
    scaling = 1e-3
  []
  [superficial_vel_y]
    type = PINSFVSuperficialVelocityVariable
    initial_condition = 1e-8
    block = ${fluid_blocks}
    scaling = 1e-3
  []
  [pressure]
    type = INSFVPressureVariable
    initial_condition = ${p_outlet}
    block = ${fluid_blocks}
    face_interp_method = average
    scaling = 10
  []
  [T_fluid]
    type = INSFVEnergyVariable
    initial_condition = ${T_Salt_initial}
    block = ${fluid_blocks}
  []
  [T_solid]
    type = INSFVEnergyVariable
    initial_condition = ${T_Salt_initial}
    block = ${solid_blocks}
  []
  [c1]
    type = MooseVariableFVReal
    block = ${fluid_blocks}
  []
  [c2]
    type = MooseVariableFVReal
    block = ${fluid_blocks}
  []
  [c3]
    type = MooseVariableFVReal
    block = ${fluid_blocks}
  []
  [c4]
    type = MooseVariableFVReal
    block = ${fluid_blocks}
  []
  [c5]
    type = MooseVariableFVReal
    block = ${fluid_blocks}
  []
  [c6]
    type = MooseVariableFVReal
    block = ${fluid_blocks}
  []
[]

Finite Volume Kernels

Furthermore, corresponding to the variables are the finite volume kernels which operate on these variables. Most of the kernels are set implicitly within the Navier Stokes action. However Kernel tweaking is possible and here the pump, delayed neutron precursor group advection equations, and solid heat conduction are explicitly set.

[FVKernels]
  # Extra kernels for the thermal-hydraulics solve in the fluid
  [pump_x]
    type = INSFVPump
    momentum_component = x
    rhie_chow_user_object = 'pins_rhie_chow_interpolator'
    variable = superficial_vel_x
    block = 'pump'
    pump_volume_force = ${pump_force}
  []
  [pump_y]
    type = INSFVPump
    momentum_component = y
    rhie_chow_user_object = 'pins_rhie_chow_interpolator'
    variable = superficial_vel_y
    block = 'pump'
    pump_volume_force = ${pump_force}
  []
  [convection_fluid_hx]
    type = NSFVEnergyAmbientConvection
    variable = T_fluid
    T_ambient = ${T_inlet}
    alpha = ${vol_hx}
    block = 'pump'
  []

  # Kernels for solve in the solid blocks
  [heat_time_solid]
    type = INSFVEnergyTimeDerivative
    variable = T_solid
    dh_dt = dh_dt
    rho = ${rho_steel}
  []
  [heat_diffusion_solid]
    type = FVDiffusion
    variable = T_solid
    coeff = ${k_steel}
  []
  [convection_core]
    type = PINSFVEnergyAmbientConvection
    variable = T_solid
    T_fluid = T_fluid
    T_solid = T_solid
    is_solid = true
    h_solid_fluid = ${bulk_hx}
    block = 'core'
  []
  [convection_core_completmeent]
    type = PINSFVEnergyAmbientConvection
    variable = T_fluid
    T_fluid = T_fluid
    T_solid = T_solid
    is_solid = false
    h_solid_fluid = ${bulk_hx}
    block = 'core'
  []

  # Kernels for solve of delayed neutron precursor transport
  [c1_time]
    type = FVFunctorTimeKernel
    variable = 'c1'
  []
  [c2_time]
    type = FVFunctorTimeKernel
    variable = 'c2'
  []
  [c3_time]
    type = FVFunctorTimeKernel
    variable = 'c3'
  []
  [c4_time]
    type = FVFunctorTimeKernel
    variable = 'c4'
  []
  [c5_time]
    type = FVFunctorTimeKernel
    variable = 'c5'
  []
  [c6_time]
    type = FVFunctorTimeKernel
    variable = 'c6'
  []
  [c1_advection]
    type = PINSFVMassAdvection
    variable = c1
    rho = 'c1_porous'
    block = ${fluid_blocks}
  []
  [c2_advection]
    type = PINSFVMassAdvection
    variable = c2
    rho = 'c2_porous'
    block = ${fluid_blocks}
  []
  [c3_advection]
    type = PINSFVMassAdvection
    variable = c3
    rho = 'c3_porous'
    block = ${fluid_blocks}
  []
  [c4_advection]
    type = PINSFVMassAdvection
    variable = c4
    rho = 'c4_porous'
    block = ${fluid_blocks}
  []
  [c5_advection]
    type = PINSFVMassAdvection
    variable = c5
    rho = 'c5_porous'
    block = ${fluid_blocks}
  []
  [c6_advection]
    type = PINSFVMassAdvection
    variable = c6
    rho = 'c6_porous'
    block = ${fluid_blocks}
  []
  [c1_turb_diffusion]
    type = INSFVMixingLengthScalarDiffusion
    schmidt_number = ${Sc_t}
    variable = c1
    block = ${fluid_blocks}
  []
  [c2_turb_diffusion]
    type = INSFVMixingLengthScalarDiffusion
    schmidt_number = ${Sc_t}
    variable = c2
    block = ${fluid_blocks}
  []
  [c3_turb_diffusion]
    type = INSFVMixingLengthScalarDiffusion
    schmidt_number = ${Sc_t}
    variable = c3
    block = ${fluid_blocks}
  []
  [c4_turb_diffusion]
    type = INSFVMixingLengthScalarDiffusion
    schmidt_number = ${Sc_t}
    variable = c4
    block = ${fluid_blocks}
  []
  [c5_turb_diffusion]
    type = INSFVMixingLengthScalarDiffusion
    schmidt_number = ${Sc_t}
    variable = c5
    block = ${fluid_blocks}
  []
  [c6_turb_diffusion]
    type = INSFVMixingLengthScalarDiffusion
    schmidt_number = ${Sc_t}
    variable = c6
    block = ${fluid_blocks}
  []
  [c1_src]
    type = FVCoupledForce
    variable = c1
    v = fission_source
    coef = ${beta1}
    block = ${fluid_blocks}
  []
  [c2_src]
    type = FVCoupledForce
    variable = c2
    v = fission_source
    coef = ${beta2}
    block = ${fluid_blocks}
  []
  [c3_src]
    type = FVCoupledForce
    variable = c3
    v = fission_source
    coef = ${beta3}
    block = ${fluid_blocks}
  []
  [c4_src]
    type = FVCoupledForce
    variable = c4
    v = fission_source
    coef = ${beta4}
    block = ${fluid_blocks}
  []
  [c5_src]
    type = FVCoupledForce
    variable = c5
    v = fission_source
    coef = ${beta5}
    block = ${fluid_blocks}
  []
  [c6_src]
    type = FVCoupledForce
    variable = c6
    v = fission_source
    coef = ${beta6}
    block = ${fluid_blocks}
  []
  [c1_decay]
    type = FVReaction
    variable = c1
    rate = ${lambda1}
    block = ${fluid_blocks}
  []
  [c2_decay]
    type = FVReaction
    variable = c2
    rate = ${lambda2}
    block = ${fluid_blocks}
  []
  [c3_decay]
    type = FVReaction
    variable = c3
    rate = ${lambda3}
    block = ${fluid_blocks}
  []
  [c4_decay]
    type = FVReaction
    variable = c4
    rate = ${lambda4}
    block = ${fluid_blocks}
  []
  [c5_decay]
    type = FVReaction
    variable = c5
    rate = ${lambda5}
    block = ${fluid_blocks}
  []
  [c6_decay]
    type = FVReaction
    variable = c6
    rate = ${lambda6}
    block = ${fluid_blocks}
  []
[]

Here the delayed neutron precursor distribution is modeled via an advection diffusion equation as follows:

(6)

where interstitial advection velocity vector at position , average molecular diffusion for neutron precursors of type , turbulent kinematic viscosity, and the turbulent Schmidt number. There are as many equations of type Eq. (6) as the number of neutron precursor groups. Since we are adding these equations in manually, the FVKernels entered correspond to the associated terms in the equation for the six delayed neutron precursor groups.

Finite Volume Interface Kernels

Additionally, finite volume interface Kernels can be deployed to account for specific interfacial modeling at boundaries. Here an explicit heat transfer correlation is implimented at the core_downcomer_boundary.

[FVInterfaceKernels]
  # Conjugated heat transfer with core barrel
  [convection]
    type = FVConvectionCorrelationInterface
    variable1 = T_fluid
    variable2 = T_solid
    boundary = 'core_downcomer_boundary'
    h = ${bulk_hx}
    T_solid = T_solid
    T_fluid = T_fluid
    subdomain1 = 'core down_comer lower_plenum upper_plenum'
    subdomain2 = 'core_barrel'
    wall_cell_is_bulk = true
  []
[]

Functions

The next block looks at setting functions which are then used in functors or initial conditions in the model. Functions give users added flexibility to edit the model and use it for their own needs.

[Functions]
  [cosine_guess]
    type = ParsedFunction
    expression = 'max(0, cos(x*pi/2/1.0))*max(0, cos((y-1.0)*pi/2/1.1))'
  []
[]

Aux Variables

Through the Auxiliary Variables, developers can define the variable type for porosity, power density, and the fission sources, along with their domain. The variable types, functions, and scaling factors are explained in detail here. Aux variables can be set explicitly in this block, or passed from other apps to perform multiphysics coupling.

The porosity variable is a constant field defined in the fluid blocks across the domain. The fluid blocks are defined at the top of the header and are substituted using the $ sign. Recall that the $ sign refers to variable substitutions.

The power density is a constant field defined in the core and plena blocks. The initial conditions for both the power density and the fission source are a cosine guess, defined below. These initial condition guesses are then udpated with the Griffin solution that is passed to Pronghorn.

[AuxVariables]
  [porosity_var]
    type = MooseVariableFVReal
    block = ${fluid_blocks}
  []
  [power_density]
    type = MooseVariableFVReal
    [FVInitialCondition]
      type = FVFunctionIC
      function = 'cosine_guess'
      scaling_factor = '${fparse 2.9183E+6}'
    []
  []
  [fission_source]
    type = MooseVariableFVReal
    [FVInitialCondition]
      type = FVFunctionIC
      function = 'cosine_guess'
      scaling_factor = '${fparse 1.0}'
    []
  []
  [rho_var]
    type = MooseVariableFVReal
    initial_condition = 1.0
    block = ${fluid_blocks}
  []
[]

Aux Kernels

Correspondingly, the Auxiliary Kernels operate on the Auxiliary Variables. Here the porosity and density auxiliary variable are incorporated into the automatic differentiation process through the ADFunctorElementalAux

[AuxKernels]
  [porosity_var_aux]
    type = FunctorAux
    variable = porosity_var
    functor = 'porosity'
    block = ${fluid_blocks}
  []
  [rho_var_aux]
    type = FunctorAux
    variable = 'rho_var'
    functor = 'rho'
    block = ${fluid_blocks}
  []
[]

Postprocessors

The postprocessor block allows for the user to calculate additional items of interest and have them included in the output. These can be useful to check assumptions that the user makes when setting up the model to ensure the model is performing as the user intends.

Here the inlet and outlet pressures and temperatures are determined as well as the volumetric flow rate and area in the downcomer.

[Postprocessors]
  [outlet_p]
    type = SideAverageValue
    variable = pressure
    boundary = 'pump_outlet'
  []
  [outlet_T]
    type = SideAverageValue
    variable = 'T_fluid'
    boundary = 'pump_outlet'
  []
  [inlet_p]
    type = SideAverageValue
    variable = 'pressure'
    boundary = 'downcomer_inlet'
  []
  [inlet_T]
    type = SideAverageValue
    variable = 'T_fluid'
    boundary = 'downcomer_inlet'
  []
  [vfr_downcomer]
    type = VolumetricFlowRate
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    advected_quantity = 1.0
    boundary = 'downcomer_inlet'
  []
  [area_pp_downcomer_inlet]
    type = AreaPostprocessor
    boundary = 'downcomer_inlet'
    execute_on = 'INITIAL'
  []
  [vfr_pump]
    type = VolumetricFlowRate
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    advected_quantity = 1.0
    boundary = 'pump_outlet'
  []
[]

# ==============================================================================
# EXECUTION PARAMETERS
# ==============================================================================
[Debug]
  show_var_residual_norms = true
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -sub_pc_factor_shift_type'
  petsc_options_value = ' lu       NONZERO'
  automatic_scaling = true
  nl_abs_tol = 1e-6
  #line_search = l2
  nl_max_its = 100

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.1
    optimal_iterations = 20
    iteration_window = 2
    growth_factor = 2
    cutback_factor = 0.5
  []

  end_time = 1e10
  steady_state_detection = true
[]

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

Materials

The Materials block allows for users to set various material parameters and correlations for the model. Here the porosity, hydraulic diameter, fluid properties, drag correlations and porosity delayed neutron precursor distributions are set.

Executioner

Finally, the Executioner block determines how the model is run. The type of solve, method, PETSc options and convergence tolerances are set in this block. Additionally the time stepping method, end time, and steady state detection options are defined in this block.

[Debug]
  show_var_residual_norms = true
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -sub_pc_factor_shift_type'
  petsc_options_value = ' lu       NONZERO'
  automatic_scaling = true
  nl_abs_tol = 1e-6
  #line_search = l2
  nl_max_its = 100

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.1
    optimal_iterations = 20
    iteration_window = 2
    growth_factor = 2
    cutback_factor = 0.5
  []

  end_time = 1e10
  steady_state_detection = true
[]

Output

Lastly, the Output block specifies what type of output (e.g. exodus and CSV) will be created and what items should be printed during the simulation.

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