Griffin steady state input

Problem Parameters

We first define in the header parameters that can be called throughout the input file:

• geometric parameters for the core height and the active core radius

• core porosity value

• burnup group information

• geometric parameters for pebbles

• parameters for the streamline calculations for pebble cycling

• the core power

• and initial values for the temperature, density, and reference density.

Global Parameters

The next block is a GlobalParams block. These parameters will be inserted in every block in the input, and are used to reduce the size of the input file. Here we specify the group cross section library. We will give additional details about the group cross sections in the Materials block.

[GlobalParams]
  library_file = '../data/gFHR_4g_pebble.xml'
  library_name = 'gFHR'
  is_meter = true
[]

Transport Systems

The [TransportSystems] block is used to specify the solver parameters.

We chose a diffusion solver as accuracy is generally satisfactory with graphite-moderated reactors, as we confirmed by benchmarking with Serpent for a similar reactor (Giudicelli et al., 2021). We select the eigenvalue equation type as the steady state coupling is an eigenvalue calculation for neutron transport. The steady state is the eigenpair of the and the steady flux distribution. We also specify the boundary conditions in this block. This is a 2D RZ model, so the center of the geometry is an axis of symmetry with a reflecting boundary condition on the left of the model. A vacuum boundary condition is placed on the other boundaries. This approximation is appropriate when the boundaries are sufficiently far from the active region.

[TransportSystems]
  particle = neutron
  equation_type = eigenvalue
  G = 4
  ReflectingBoundary = 'left'
  VacuumBoundary = 'right top bottom'

  [diff]
    scheme = CFEM-Diffusion
    family = LAGRANGE
    order = SECOND
    n_delay_groups = 6
    assemble_scattering_jacobian = true
    assemble_fission_jacobian = true
  []
[]

Mesh

The next part of the input specifies the geometry. Here we use the CartesianMeshGenerator to create our own mesh, specifying 5 blocks, including the pebble_bed, reflector, ss stainless steel vessel, downcomer, and cr the control rod. Here we have a 2D, RZ model, and use the subdomain_id to draw our mesh with dx and dy setting the width of the corresponding ix and iy which determine the number of elements in the mesh matching what is drawn in subdomain_id.

The neutronics mesh is shown in Figure 1.

Figure 1: Griffin equilibrium core mesh from (Ortensi et al., 2023).

We then edit the mesh, through the SubdomainExtraElementIDGenerator to define material_ids and subsequently use RenameBlockGenerator to merge block 2 into block 1.

[Mesh]
  block_id = '     1        3          4      5      6'
  block_name = 'pebble_bed  reflector   ss downcomer  cr'
  [cartesian_mesh]
    type = CartesianMeshGenerator
    dim = 2

    dx = '0.9 0.3 0.156 0.444 0.02 0.05 0.04'
    ix = '  6   2    1     3    1     1    1'

    dy = '0.6 0.30947 2.47576 0.30947 0.6'
    iy = '  4    2       16       2     4'

    subdomain_id = '
   3 3 3 3 4 5 4
   2 2 6 3 4 5 4
   1 2 6 3 4 5 4
   2 2 6 3 4 5 4
   3 3 6 3 4 5 4'
  []

  [assign_material_id]
    type = SubdomainExtraElementIDGenerator
    subdomains = '1 2 3 4 5 6 '
    extra_element_id_names = 'material_id'
    extra_element_ids = '1 2 3 4 5 2'
    input = cartesian_mesh
  []
  [merge]
    type = RenameBlockGenerator
    input = assign_material_id
    old_block = '1 2 3 4 5 6'
    new_block = '1 1 3 4 5 6'
  []
  coord_type = RZ
  second_order = true
[]

Auxiliary Variables

The Variables and Kernels need not be defined when using Griffin, as the TransportSystems block is an Action in MOOSE vocabulary that takes care of defining those. However it is helpful to define AuxVariables which can be used to define fields/variables that are not solved for. They may be used for coupling, for computing material properties, for outputting quantities of interests, and many other uses. In this input file, Tsolid and Rho are used to incorporate fields from the thermal hydraulics simulation for a multiphysics coupled solution. These fields are populated by Transfers from the Pronghorn app, and will be discussed further in the input file. We also define other auxiliary variables to track items like burnup, porosity, and power densities or temperatures like pden_avg and Tmod_max.

[AuxVariables]
  [Tsolid]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = '${initial_temperature}'
  []
  [Rho]
    order = CONSTANT
    family = MONOMIAL
    initial_condition = '${Rho}'
  []
  [Tfuel_avg]
    order = CONSTANT
    family = MONOMIAL
  []
  [Tfuel_max]
    order = CONSTANT
    family = MONOMIAL
  []
  [Tmod_avg]
    order = CONSTANT
    family = MONOMIAL
  []
  [Tmod_max]
    order = CONSTANT
    family = MONOMIAL
  []
  [Burnup]
    order = CONSTANT
    family = MONOMIAL
    components = 9
    # use burnup midpoints and add one additional for discard group
    initial_condition = ${burnup_group_avg}
    # outputs = none
  []
  [Burnup_avg]
    order = CONSTANT
    family = MONOMIAL
  []
  [porosity]
    order = CONSTANT
    family = MONOMIAL
  []
  [pden_avg]
    order = CONSTANT
    family = MONOMIAL
  []
  [pden_max]
    order = CONSTANT
    family = MONOMIAL
  []
  [power_peaking]
    order = CONSTANT
    family = MONOMIAL
  []
  [rod_position]
    order = CONSTANT
    family = MONOMIAL
    block = 'cr'
  []
[]

Auxiliary Kernels

Next, we define AuxKernels which operate on AuxVariables. They may scale, multiply, add and perform many other operations. They may be block restricted, as some AuxVariables are not defined over the entire domain, or they may inherit the same block restriction as the variable. Some examples here are Auxiliary Variables which find the max or average value pden_max and Tfuel_avg using the ArrayVarExtremeValueAux and PebbleAveragedAux respectively. We also set the porosity auxiliary variable through a FunctionAux which reads in a function called porosity_f which will be discussed next.

[AuxKernels]
  [pden_max]
    type = ArrayVarExtremeValueAux
    block = 'pebble_bed'
    variable = pden_max
    array_variable = partial_power_density
    value_type = MAX
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [pden_avg]
    type = PebbleAveragedAux
    block = 'pebble_bed'
    variable = pden_avg
    array_variable = partial_power_density
    pebble_volume_fraction = pebble_volume_fraction
    n_fresh_pebble_types = 1
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [power_peaking_aux]
    type = ParsedAux
    block = 'pebble_bed'
    variable = power_peaking
    coupled_variables = 'total_power_density pden_avg'
    expression = 'total_power_density / pden_avg'
  []
  [Tfuel_avg]
    type = PebbleAveragedAux
    block = 'pebble_bed'
    variable = Tfuel_avg
    array_variable = triso_temperature
    pebble_volume_fraction = pebble_volume_fraction
    n_fresh_pebble_types = 1
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tfuel_max]
    type = ArrayVarExtremeValueAux
    block = 'pebble_bed'
    variable = Tfuel_max
    array_variable = triso_temperature
    value_type = MAX
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tmod_avg]
    type = PebbleAveragedAux
    block = 'pebble_bed'
    variable = Tmod_avg
    array_variable = graphite_temperature
    pebble_volume_fraction = pebble_volume_fraction
    n_fresh_pebble_types = 1
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tmod_max]
    type = ArrayVarExtremeValueAux
    block = 'pebble_bed'
    variable = Tmod_max
    array_variable = graphite_temperature
    value_type = MAX
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Burnup_avg]
    type = PebbleAveragedAux
    block = 'pebble_bed'
    variable = Burnup_avg
    array_variable = Burnup
    pebble_volume_fraction = pebble_volume_fraction
    n_fresh_pebble_types = 1
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [porosity]
    type = FunctionAux
    block = 'pebble_bed'
    variable = porosity
    function = porosity_f
    execute_on = 'INITIAL'
  []
[]

Functions

Functions may be defined in MOOSE in a Functions block. Here a function to define the porosity is read in using the PiecewiseMulticonstant function to read a data file defined in the text file gFHR_porosity.txt. Additionally the control rod position is set in this input using a ConstantFunction.

[Functions]
  [porosity_f]
    type = PiecewiseMulticonstant
    direction = 'right right'
    data_file = '../data/gFHR_porosity.txt'
  []
  [CR_pos_f]
    type = ConstantFunction
    value = '3.38523' #Note that when position is 3.6947 CR is withdrawn
  []
[]

Pebble Depletion

Now we introduce a new capability developed in Griffin for this specific model to do PBFHR depletion with streamlines for pebble shuffling and achieve an equilibrium core. This PebbleDepletion action is quite extensive and the theory can be found here.

Therefore, we will very briefly cover this block at a high level. First the Pebble Depletion action reads in variables and Auxiliary variables such as the power, burnup, Tfuel, Tmod, and Rho. The coolant and pebble compositions are read in from the Compositions block which will be covered shortly. Transmutation data in the ISOXML format is read in for this problem from the file DRAGON5_DT.xml. Specific Pebble and Tabulation options are also specified.

[PebbleDepletion]
  block = 'pebble_bed'

  power = ${total_power}
  integrated_power_postprocessor = total_power
  power_density_variable = total_power_density
  family = MONOMIAL
  order = CONSTANT

  burnup_grid_name = 'Burnup'
  fuel_temperature_grid_name = 'Tfuel'
  moderator_temperature_grid_name = 'Tmod'
  additional_grid_name_variable_mapping = 'Rho Rho'

  # transmutation data
  dataset = ISOXML
  isoxml_data_file = '../data/DRAGON5_DT.xml'
  isoxml_lib_name = 'DRAGON'
  dtl_physicality = DISABLE

  # Isotopic options
  n_fresh_pebble_types = 1
  fresh_pebble_compositions = pebble
  track_isotopes = 'U235 U236 U238 PU238 PU239 PU240 PU241 PU242 AM241 AM242M CS135 CS137 XE135 XE136 I131 I135 SR90'
  decay_heat = true

  # Pebble options
  porosity_name = porosity
  burnup_group_boundaries = ${burnup_group_boundaries}

  # Tabulation data.
  initial_moderator_temperature = '${initial_temperature}'
  initial_fuel_temperature = '${initial_temperature}'

  [DepletionScheme]
    type = ConstantStreamlineEquilibrium

    # Streamline definition
    major_streamline_axis = y
    streamline_points = '0.152  0.60 0  0.152  3.6947 0;
                         0.450  0.60 0  0.450  3.6947 0;
                         0.750  0.60 0  0.750  3.6947 0;
                         1.050  0.60 0  1.050  3.6947 0'
    streamline_segment_subdivisions = '20; 20; 20; 20' # align with the neutronics mesh

    # Pebble options
    pebble_unloading_rate = ${pebble_unloading_rate}
    pebble_flow_rate_distribution = ${streamline_flow_fraction}
    pebble_diameter = '${fparse pebble_radius * 2.0}'
    burnup_limit = '${burnup_limit}'

    # Solver parameters
    sweep_tol = 1e-8
    sweep_max_iterations = 100

    # Output
    exodus_streamline_output = false
  []

  # pebble conduction
  pebble_conduction_input_file = ${pebble_conduction_input_file}
  pebble_positions_file = '../data/pebble_heat_pos_8r_20z.txt'
  surface_temperature_sub_app_postprocessor = T_surface
  surface_temperature_main_app_variable = Tsolid
  power_sub_app_postprocessor = pebble_power_density_pp
  fuel_temperature_sub_app_postprocessor = T_fuel
  moderator_temperature_sub_app_postprocessor = T_mod

  # coolant settings
  coolant_composition_name = coolant
  coolant_density_variable = 'Rho'
  coolant_density_ref = 1973.8 # kg/m^3
  # TODO: add material id 2 with coolant as well
  coolant_material_id = '1'
[]

More specifically, the specific DepletionScheme that is used here is the ConstantStreamLineEquilibrium option. Here the streamlines of the shuffling pebbles are used to correctly determine the equilibrium depletion.

[PebbleDepletion]
  [DepletionScheme]
    type = ConstantStreamlineEquilibrium

    # Streamline definition
    major_streamline_axis = y
    streamline_points = '0.152  0.60 0  0.152  3.6947 0;
                         0.450  0.60 0  0.450  3.6947 0;
                         0.750  0.60 0  0.750  3.6947 0;
                         1.050  0.60 0  1.050  3.6947 0'
    streamline_segment_subdivisions = '20; 20; 20; 20' # align with the neutronics mesh

    # Pebble options
    pebble_unloading_rate = ${pebble_unloading_rate}
    pebble_flow_rate_distribution = ${streamline_flow_fraction}
    pebble_diameter = '${fparse pebble_radius * 2.0}'
    burnup_limit = '${burnup_limit}'

    # Solver parameters
    sweep_tol = 1e-8
    sweep_max_iterations = 100

    # Output
    exodus_streamline_output = false
  []
[]

Lastly, the pebble conduction problem is solved through a sub-app gFHR_pebble_triso_ss.i where the positions are loaded in from the text data file pebble_heat_pos_16r_40z.txt. Specific postprocessors from that sub-app are then read in and used in the PebbleDepletion action. See this page for more information on the multiscale pebble - triso fuel performance model.

Compositions

Next, specific Compositions are defined that will be read in by the PebbleDepletion actions. Here IsotopeComposition is used to define the coolant and pebble isotopic compositions. These compositions are then loaded into the PebbleDepletion action previously described.

Materials

The group cross sections are distributed through the geometry using the Materials block. Each block in this section is a Material object defining the group cross section in a block.

Here there are two types of neutronics materials including the CoupledFeedbackRoddedNeutronicsMaterial and the CoupledFeedbackMatIDNeutronicsMaterial. Here the CR_dynamic material is used to define the control rod, and the downcomer and non-pebble-bed-regions correspondingly identify the material and the auxiliary variables that are read in from the Pronghorn sub-app.

Specifically, the CoupledFeedbackMatIDNeutronicsMaterial is used to take into account the effect of the fuel and salt temperature on the group cross section. The cross section library name and the file containing it are specified using the GlobalParams block. The temperature dependence is captured using a tabulation. The names used in the tabulation are matched to the variables in the simulation using grid_names and grid_variables. The correct entry in the library is selected using the material_id, except this was read in previously from material IDs in the mesh in assign_material_id.

[Materials]
  [CR_dynamic]
    type = CoupledFeedbackRoddedNeutronicsMaterial
    block = 'cr'
    grid_names = 'Tsolid Tsolid Tsolid'
    grid_variables = 'Tsolid Tsolid Tsolid'
    isotopes = 'pseudo; pseudo; pseudo'
    densities = '1.0 1.0 1.0'
    rod_segment_length = 3.0947
    front_position_function = 'CR_pos_f'
    segment_material_ids = '7 6 7'
    rod_withdrawn_direction = 'y'
    plus = true
  []
  [downcomer]
    type = CoupledFeedbackMatIDNeutronicsMaterial
    block = 'downcomer'
    grid_names = 'Rho'
    grid_variables = 'Rho'
    isotopes = 'pseudo'
    densities = '1.0'
  []
  [non-pebble-bed-regions]
    type = CoupledFeedbackMatIDNeutronicsMaterial
    block = 'reflector ss'
    grid_names = 'Tsolid'
    grid_variables = 'Tsolid'
    isotopes = 'pseudo'
    densities = '1.0'
  []
[]

Executioner

The Executioner block specifies how to solve the equation system. We choose an eigenvalue executioner, Eigenvalue, as a criticality calculation is an eigenvalue problem. The solve type uses the PJFNKMO or the Preconditioned Jacobian-Free Newton_krylov - Matrix Only methods of solving the eigenvalue problem.

The number of non-linear iterations and the non-linear relative and absolute convergence criteria may be respectively reduced and loosened to reduce the computational cost of the solution at the expense of its convergence. Additionally, for multiphysics coupled solutions, the fixed point iteration convergence criteria can also be specified.

[Executioner]
  type = Eigenvalue
  solve_type = PJFNKMO

  #Linear/nonlinear iterations.
  nl_max_its = 50
  nl_abs_tol = 1e-6

  fixed_point_force_norms = true
  fixed_point_max_its = 50
  fixed_point_abs_tol = 1e-6
[]

MultiApps

The Pronghorn sub-application is created by the MultiApps block with the pronghorn input file being gFHR_pronghorn_ss.i. Since we are seeking a steady state solution, we want the thermal hydraulics problem to be fully solved at every multiphysics iteration. This is accomplished by a FullSolveMultiApp. The execute_on field is set to timestep_end and FINAL to have the thermal hydraulics solve be performed after the neutronics solve. This means that for the first multiphysics iteration, the temperature fields in the neutronics solve will be set to an initial guess, while the power distribution in the thermal hydraulics solve will be updated once already. 0 -0.09 0 is specified for the position since the two meshes are slightly different and not aligned.

[MultiApps]
  [flow]
    type = FullSolveMultiApp
    input_files = ${flow_subapp_input_file}
    keep_solution_during_restore = true
    positions = '0 -0.09 0'
    execute_on = 'TIMESTEP_END FINAL'
    max_procs_per_app = 12
  []
[]

Transfers

The coupling to the thermal hydraulics simulation is done by using a MultiAppGeneralFieldNearestLocationTransfer of the power density to the pronghorn sub-app. This transfer is conservative, in that the total power in both applications is preserved. This is important for the multiphysics coupling, as the power will not fluctuate in the thermal hydraulics simulation based on the power density profile and the difference between the neutronics and thermal hydraulics meshes.

The conservation of the transfer is achieved through the specification of two postprocessors, one in the source and one in the receiving simulations, which are used to compute a scaling factor.

Additionally, the pronghorn sub-application sends steady state thermal hydraulics solutions of the temperature Tsolid and density Rho back to the neutronics solve to update the new flux, and eigenvalue until the multiphysics solution converges within tolerance set in the Executioner.

[Transfers]
  #Pronghorn flow simulation transfer
  [power_density_to_flow]
    type = MultiAppGeneralFieldNearestLocationTransfer
    to_multi_app = flow
    source_variable = total_power_density
    variable = power_density
    from_postprocessors_to_be_preserved = total_power
    to_postprocessors_to_be_preserved = total_power
  []
  [Tsolid_from_flow]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = flow
    source_variable = T_solid
    variable = Tsolid
  []
  [Rho_from_flow]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = flow
    source_variable = rho_var
    variable = Rho
  []
[]

Postprocessors

The second to last block is the Postprocessors block which specifies calculations that should be made regarding the solution so that data analysis and interpretation of results can be performed. Key items of interest here are the ElementAverageValue for variables of interest like the temperature in the core, or density of the coolant in the core. Similarly ElementExtremeValue finds the maximum value of certain variables on the mesh and reports them. These postprocessors will be displayed in the output to terminal, but will also be stored in the CSV if that option is selected in the Outputs block.

[Postprocessors]
  [power_peak]
    type = ElementExtremeValue
    variable = power_peaking
    value_type = max
    block = 'pebble_bed'
  []
  [Tsolid_core]
    type = ElementAverageValue
    block = 'pebble_bed'
    variable = Tsolid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Rho_core]
    type = ElementAverageValue
    block = 'pebble_bed'
    variable = Rho
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tfuel_avg]
    type = ElementAverageValue
    block = 'pebble_bed'
    variable = Tfuel_avg
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tfuel_max]
    type = ElementExtremeValue
    block = 'pebble_bed'
    variable = Tfuel_max
    value_type = max
    execute_on = 'INITIAL TIMESTEP_END TRANSFER'
  []
  [Tmod_avg]
    type = ElementAverageValue
    block = 'pebble_bed'
    variable = Tmod_avg
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tmod_max]
    type = ElementExtremeValue
    block = 'pebble_bed'
    variable = Tmod_max
    value_type = max
    execute_on = 'INITIAL TIMESTEP_END TRANSFER'
  []
  [Burnup_avg]
    type = ElementAverageValue
    block = 'pebble_bed'
    variable = Burnup_avg
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tsolid_reflector]
    type = ElementAverageValue
    block = 'reflector'
    variable = Tsolid
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

Outputs

Finally, the Outputs block indicates which types of outputs the simulation should return. We specify here two major types:

• exodus. Exodus files contain the mesh, the (aux)variables and global quantities such as postprocessors. We can use Paraview to view the spatial dependence of each field. They may be used to restart simulations.

• comma-separated values, or csv. This file reports the values of the postprocessors at each time step. They are useful for easily importing then plotting the simulation results using Python or Matlab.

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