Micro Reactor Drum Rotation

Contact: Zachary Prince, zachmprince.at.gmail.com, Javier Ortensi javier.ortensi.at.inl.gov

Model link: Micro Reactor Drum Rotation Model

Reactor Description

This reactor model is based on the prototypical micro reactor presented in Terlizzi and Labouré (2023). The design is based on the Empire reactor (Matthews et al., 2021) with modifications made to ensure a negative overall temperature reactivity coefficient, necessary to simulate realistic transients. This 2 MWth core consists of 18 hexagonal assemblies arranged in two rings and surrounded by 12 control drums contained within a radial reflector. Each assembly is composed of 96 fuel pins, 60 moderator pins, and 61 heat pipes. Detailed dimensions and material composition can be found in Terlizzi and Labouré (2023).

This model has been simplified to a 2D geometry that leverages symmetry to reduce the domain to a 1/12th of the full core. Furthermore, the heat pipes are simulated using a convective boundary conditions, as opposed to using Sockeye to calculate the actual heat removal. Prince et al. (2023) presents the 3D version of this model.

The purpose of this exposition is to show how to perform a multiphysics transient simulation involving control drum rotation using Griffin. Initially, the drum is positioned halfway at 90 degrees. The drum is then rotated outward (inserting reactivity) at 20 degrees per second for two seconds, then rotated inward (removing reactivity) at 20 degrees per second for three seconds. Thus, the total transient duration is five seconds. Due to the symmetry imposed on the model, rotating this drum is equivalent to rotating all the drums in the core.

Model Description

Mesh

Two different meshes were created for this model: a fine-mesh representing the exact geometry and a coarse-mesh for CMFD.

Fine Mesh

The following creates a mesh of the entire core:

Listing 1: Mesh blocks for full-core fine mesh

[Mesh]
  [FUEL_pin]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2'
    background_intervals = 1
    background_block_ids = '8'
    background_block_names = 'MONOLITH'
    polygon_size = 1.075
    polygon_size_style = 'apothem'
    ring_radii = '1.0'
    ring_intervals = '3'
    ring_block_ids = '1 2'
    ring_block_names = 'FUEL FUEL_QUAD'
    preserve_volumes = on
    quad_center_elements = false
  []

  [MOD_pin]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2'
    background_intervals = 1
    background_block_ids = '8'
    background_block_names = 'MONOLITH'
    polygon_size = 1.075
    polygon_size_style = 'apothem'
    ring_radii = '0.975 1'
    ring_intervals = '3 1'
    ring_block_ids = '3 4 5'
    ring_block_names = 'MOD MOD_QUAD MGAP'
    preserve_volumes = on
    quad_center_elements = false
  []

  [HPIPE_pin]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2'
    background_intervals = 1
    background_block_ids = '8'
    background_block_names = 'MONOLITH'
    polygon_size = 1.075
    polygon_size_style = 'apothem'
    ring_radii = '1.0'
    ring_intervals = '3'
    ring_block_ids = '6 7'
    ring_block_names = 'HPIPE HPIPE_QUAD'
    preserve_volumes = on
    quad_center_elements = false
  []

  [AIRHOLE_CELL]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2'
    background_intervals = 3
    background_block_ids = '20 21'
    background_block_names = 'AIRHOLE AIRHOLE_QUAD'
    polygon_size = 1.075
    polygon_size_style = 'apothem'
    preserve_volumes = on
    quad_center_elements = false
  []

  [REFL_CELL]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2'
    background_intervals = 2
    background_block_ids = '14 15'
    background_block_names = 'REFL REFL_QUAD'
    polygon_size = 1.075
    polygon_size_style = 'apothem'
    preserve_volumes = on
    quad_center_elements = false
  []

  [Assembly_1]
    type = PatternedHexMeshGenerator
    inputs = 'MOD_pin HPIPE_pin FUEL_pin'
    # Pattern ID  0        1        2
    background_intervals = 1
    background_block_id = '8'
    background_block_name = 'MONOLITH'
    duct_sizes = '16.0765'
    duct_sizes_style = 'apothem'
    duct_block_ids = '22'
    duct_block_names = 'AIR'
    duct_intervals = 1
    hexagon_size = '16.1765'
    pattern = '1 0 1 0 1 0 1 0 1;
              0 2 2 2 2 2 2 2 2 0;
             1 2 1 0 1 0 1 0 1 2 1;
            0 2 0 2 2 2 2 2 2 0 2 0;
           1 2 1 2 1 0 1 0 1 2 1 2 1;
          0 2 0 2 0 2 2 2 2 0 2 0 2 0;
         1 2 1 2 1 2 1 0 1 2 1 2 1 2 1;
        0 2 0 2 0 2 0 2 2 0 2 0 2 0 2 0;
       1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1;
        0 2 0 2 0 2 0 2 2 0 2 0 2 0 2 0;
         1 2 1 2 1 2 1 0 1 2 1 2 1 2 1;
          0 2 0 2 0 2 2 2 2 0 2 0 2 0;
           1 2 1 2 1 0 1 0 1 2 1 2 1;
            0 2 0 2 2 2 2 2 2 0 2 0;
             1 2 1 0 1 0 1 0 1 2 1;
              0 2 2 2 2 2 2 2 2 0;
               1 0 1 0 1 0 1 0 1'
  []

  [AIRHOLE]
    type = PatternedHexMeshGenerator
    inputs = 'AIRHOLE_CELL'
    # Pattern ID       0
    background_intervals = 1
    background_block_id = '21'
    background_block_name = 'AIRHOLE_QUAD'
    duct_sizes = '16.0765'
    duct_sizes_style = 'apothem'
    duct_block_ids = '22'
    duct_block_names = 'AIR'
    duct_intervals = 1
    hexagon_size = '16.1765'
    pattern = '0 0 0 0 0 0 0 0 0;
              0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0 0;
            0 0 0 0 0 0 0 0 0 0 0 0;
           0 0 0 0 0 0 0 0 0 0 0 0 0;
          0 0 0 0 0 0 0 0 0 0 0 0 0 0;
         0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
       0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
         0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
          0 0 0 0 0 0 0 0 0 0 0 0 0 0;
           0 0 0 0 0 0 0 0 0 0 0 0 0;
            0 0 0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0 0;
              0 0 0 0 0 0 0 0 0 0;
               0 0 0 0 0 0 0 0 0'
  []

  [REFL]
    type = PatternedHexMeshGenerator
    inputs = 'REFL_CELL'
    # Pattern ID    0
    background_intervals = 1
    background_block_id = '15'
    background_block_name = 'REFL_QUAD'
    duct_sizes = '16.0765'
    duct_sizes_style = 'apothem'
    duct_block_ids = '22'
    duct_block_names = 'AIR'
    duct_intervals = 1
    hexagon_size = '16.1765'
    pattern = '0 0 0 0 0 0 0 0 0;
              0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0 0;
            0 0 0 0 0 0 0 0 0 0 0 0;
           0 0 0 0 0 0 0 0 0 0 0 0 0;
          0 0 0 0 0 0 0 0 0 0 0 0 0 0;
         0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
       0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
         0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
          0 0 0 0 0 0 0 0 0 0 0 0 0 0;
           0 0 0 0 0 0 0 0 0 0 0 0 0;
            0 0 0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0 0;
              0 0 0 0 0 0 0 0 0 0;
               0 0 0 0 0 0 0 0 0'
  []

  [cd_0]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Assembly_1 Assembly_1 Assembly_1 Assembly_1 Assembly_1 Assembly_1'
    sides_to_adapt = '0 1 2 3 4 5'
    num_sectors_per_side = '8 8 8 8 8 8'
    hexagon_size = 16.1765
    background_intervals = 3
    background_block_ids = 11
    background_block_names = 'CRREFL_QUAD'
    ring_radii = '13.8 14.8'
    ring_intervals = '5 3'
    ring_block_ids = '10 11 13'
    ring_block_names = 'CRREFL CRREFL_QUAD CRREFL_DYNAMIC'
    preserve_volumes = true
    is_control_drum = true
    # quad_center_elements = true
    # center_quad_factor = 0.1
  []

  [Core]
    type = PatternedHexMeshGenerator
    inputs = 'Assembly_1 Assembly_1 Assembly_1 AIRHOLE REFL cd_0'
    # Pattern ID   0            1          2         3     4    5
    generate_core_metadata = true
    pattern_boundary = none
    pattern = '4 4 4 4 4;
              4 4 5 5 4 4;
             4 5 1 2 1 5 4;
            4 5 2 0 0 2 5 4;
           4 4 1 0 3 0 1 4 4;
            4 5 2 0 0 2 5 4;
             4 5 1 2 1 5 4;
              4 4 5 5 4 4;
               4 4 4 4 4'
  []
[]
(microreactors/drum_rotation/empire_2d_CD_fine.i)

Then the mesh is "sliced" to produce a 1/12th geometry and scaled to meters:

Listing 2: Mesh blocks trimming full-core fine mesh to 1/12th geometry and scaling to meters

[Mesh]
  [half]
    type = PlaneDeletionGenerator
    point = '0 0 0'
    normal = '0 -1 0'
    input = Core
  []

  [twelfth]
    type = PlaneDeletionGenerator
    point = '0 0 0'
    normal = '${fparse -cos(pi/3)} ${fparse sin(pi/3)} 0'
    input = half
  []

  [trim]
    type = PlaneDeletionGenerator
    point = '84.0556 48.5295 0'
    normal = '84.0556 48.5295 0'
    input = twelfth
  []

  [scale]
    type = TransformGenerator
    input = trim
    transform = SCALE
    vector_value = '0.01 0.01 0.01'
  []
[]
(microreactors/drum_rotation/empire_2d_CD_fine.i)

However, performing the transient with the mesh so far causes a significant cusping effect (Yamamoto, 2004) in the calculated power and reactivity, even when using the cusping treatment described in Schunert et al. (2019). Investigations indicate that the cusping effect is due to the discretization when the cross sections of the absorbing material and the non-absorbing material are significantly different. It is possible to mitigate the cusping effect by increasing the spatial expansion order in the drum region locally. However, for this model, the mesh in the drum region was refined such that drum's front would align with element edges at every time step. We plan on addressing this cusping issue in future work since this refinement limits the flexibility of time step sizes and unnecessarily increases the number of elements.

The mesh for the drum is created separately with the following input:

Listing 3: Input file generating fine-mesh for control drum

[Mesh]
  [drum]
    type = ConcentricCircleMeshGenerator
    num_sectors = 450
    radii = '13.8 14.8'
    rings = '1 3'
    preserve_volumes = true
    has_outer_square = false
    portion = top_half
  []
  [rename_drum]
    type = RenameBlockGenerator
    input = drum
    old_block = 2
    new_block = 13
  []
  [remove_inner]
    type = BlockDeletionGenerator
    input = rename_drum
    block = 1
    new_boundary = 'boundary_to_fill'
  []
  [inner]
    type = ConcentricCircleMeshGenerator
    num_sectors = 50
    radii = '11 14.8'
    rings = '1 3'
    preserve_volumes = false
    has_outer_square = false
    portion = top_half
  []
  [rename_inner]
    type = RenameBlockGenerator
    input = inner
    old_block = 1
    new_block = 11
  []
  [remove_outer]
    type = BlockDeletionGenerator
    input = rename_inner
    block = 2
    new_boundary = 'boundary_to_fill'
  []
  [fill]
    type = FillBetweenSidesetsGenerator
    input_mesh_1 = remove_inner
    input_mesh_2 = remove_outer
    boundary_1 = boundary_to_fill
    boundary_2 = boundary_to_fill
    num_layers = 5
  []
  [remove_boundaries]
    type = BoundaryDeletionGenerator
    input = fill
    boundary_names = 'outer bottom'
  []
  [add_bottom_boundary]
    type = SideSetsFromNormalsGenerator
    input = remove_boundaries
    normals = '0 -1 0'
    new_boundary = sym
  []
  [mirror]
    type = SymmetryTransformGenerator
    input = add_bottom_boundary
    mirror_normal_vector = '0 -1 0'
    mirror_point = '0 0 0'
  []
  [stitch]
    type = StitchedMeshGenerator
    inputs = 'add_bottom_boundary mirror'
    stitch_boundaries_pairs = 'sym sym'
    clear_stitched_boundary_ids = true
  []
  [rename_block]
    type = RenameBlockGenerator
    input = stitch
    old_block = 1
    new_block = 10
  []
[]
(microreactors/drum_rotation/drum.i)

This drum mesh is then loaded in by the main mesh input, stitched into the previous control drum assembly, and inserted into current 1/12th model (replacing the previous assembly):

Listing 4: Mesh blocks inserting pregenerated control drum in 1/12th fine mesh

[Mesh]
  [drum]
    type = FileMeshGenerator
    file = drum_in.e
  []

  [drum_insert]
    type = XYDelaunayGenerator
    boundary = cd_0
    holes = 'drum'
    stitch_holes = true
    refine_holes = false
  []

  [drum_scale]
    type = TransformGenerator
    input = drum_insert
    transform = SCALE
    vector_value = '0.01 0.01 0.01'
  []

  [drum_rotate]
    type = TransformGenerator
    input = drum_scale
    transform = ROTATE
    vector_value = '0 0 30'
  []

  [drum_move]
    type = TransformGenerator
    input = drum_rotate
    transform = TRANSLATE
    vector_value = '${x_center} ${y_center} 0'
  []

  [drum_blocks]
    type = RenameBlockGenerator
    input = drum_move
    old_block = '0'
    new_block = '10'
  []

  [drum_boundary]
    type = SideSetsAroundSubdomainGenerator
    input = drum_blocks
    block = 10
    new_boundary = 'drum_outer'
    include_only_external_sides = true
  []

  [drum_remove]
    type = BlockDeletionGenerator
    input = scale
    block = 'CRREFL CRREFL_QUAD CRREFL_DYNAMIC'
    new_boundary = 'drum_inner'
  []

  [stitch]
    type = StitchedMeshGenerator
    inputs = 'drum_boundary drum_remove'
    stitch_boundaries_pairs = 'drum_outer drum_inner'
    clear_stitched_boundary_ids = false
  []
[]
(microreactors/drum_rotation/empire_2d_CD_fine.i)

Finally, side sets are added to the mesh so that boundary conditions can be properly applied:

Listing 5: Adding boundaries to the fine mesh

[Mesh]
  [bottom_boundary]
    type = SideSetsFromNormalsGenerator
    input = stitch
    normals = '0 -1 0'
    new_boundary = 'bottom'
  []

  [topleft_boundary]
    type = SideSetsFromNormalsGenerator
    input = bottom_boundary
    normals = '${fparse -cos(pi/3)} ${fparse sin(pi/3)} 0'
    new_boundary = 'topleft'
  []

  [right_boundary]
    type = SideSetsFromNormalsGenerator
    input = topleft_boundary
    normals = '84.0556 48.5295 0'
    new_boundary = 'right'
  []
[]
(microreactors/drum_rotation/empire_2d_CD_fine.i)

The resulting mesh is shown in Figure 1.

Figure 1: Fine mesh used for control drum rotation transient

A list of all the blocks is shown in Table 1.

Table 1: Description of blocks in geometry

RegionBlock IDs
Fuel1 2
Moderator3 4 5
Heat Pipes6 7
Monolith8
Reflector10 11 14 15
Control Drum13
Air20 21 22

Coarse Mesh

The coarse mesh used for CMFD does not resolve pins, instead it creates quad4 meshes for each hexagonal assembly:

Listing 6: Mesh blocks for full-core coarse mesh

[Mesh]
  [F1_1]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2'
    background_intervals = 8
    background_block_ids = '100 101'
    background_block_names = 'F1 F1_QUAD'
    duct_sizes = '14.8956'
    duct_sizes_style = apothem
    duct_block_ids = '101'
    duct_block_names = 'F1_QUAD'
    duct_intervals = 1
    polygon_size = 16.1765
    polygon_size_style = 'apothem'
    preserve_volumes = on
    quad_center_elements = false
  []

  [F1_2]
    type = PolygonConcentricCircleMeshGenerator
    num_sides = 6 # must be six to use hex pattern
    num_sectors_per_side = '2 2 2 2 2 2'
    polygon_size = 16.1765
    polygon_size_style = 'apothem'
    background_intervals = 1
    background_block_ids = 101
    background_block_names = 'F1_QUAD'
    ring_radii = '13.8'
    ring_intervals = '5'
    ring_block_ids = '100 101'
    ring_block_names = 'F1 F1_QUAD'
    preserve_volumes = off
    quad_center_elements = false
  []

  [Core_CM]
    type = PatternedHexMeshGenerator
    inputs = 'F1_1 F1_2'
    # Pattern ID 0  1
    pattern_boundary = none
    pattern = '0 0 0 0 0;
              0 0 1 1 0 0;
             0 1 0 0 0 1 0;
            0 1 0 0 0 0 1 0;
           0 0 0 0 0 0 0 0 0;
            0 1 0 0 0 0 1 0;
             0 1 0 0 0 1 0;
              0 0 1 1 0 0;
               0 0 0 0 0'
  []
[]
(microreactors/drum_rotation/empire_2d_CD_coarse.i)

Same as the fine mesh, this full core mesh is trimmed to a 1/12th geometry and scaled to meters:

Listing 7: Mesh blocks trimming full-core coarse mesh to 1/12th geometry and scaling to meters

[Mesh]
  [half_CM]
    type = PlaneDeletionGenerator
    point = '0 0 0'
    normal = '0 -1 0'
    input = Core_CM
    new_boundary = bottom
  []

  [twelfth_CM]
    type = PlaneDeletionGenerator
    point = '0 0 0'
    normal = '${fparse -cos(pi/3)} ${fparse sin(pi/3)} 0'
    input = half_CM
    new_boundary = topleft
  []

  [trim_CM]
    type = PlaneDeletionGenerator
    point = '84.0556 48.5295 0'
    normal = '84.0556 48.5295 0'
    input = twelfth_CM
    new_boundary = right
  []

  [scale_CM]
    type = TransformGenerator
    input = trim_CM
    transform = SCALE
    vector_value = '0.01 0.01 0.01'
  []
[]
(microreactors/drum_rotation/empire_2d_CD_coarse.i)

The resulting coarse mesh is shown in Figure 2.

Figure 2: Coarse mesh used for control drum rotation transient

Materials

Cross Sections

Serpent (v. 2.1.32) was used to generate the multigroup cross sections for the SiMBA problem. The ENDF/B-VIII.0 continuous energy library was utilized to leverage the scattering libraries, which was then converted into an 11-group structure to perform the calculations. The cross sections were parametrized with respect to three variables: (1) fuel temperature; (2) moderator, monolith, heat pipe, and reflector temperature; and (3) control drum angle. The grid points selected can be seen in the multigroup cross-section library file:

Listing 8: Cross section tabulation grid

      <Tabulation>Tfuel Tmod CD</Tabulation>
      <Tfuel>800 1000 1200 1400 1600</Tfuel>
      <Tmod>800 1000 1200</Tmod>
      <CD>0 60 120 180</CD>
(microreactors/drum_rotation/empire_core_modified_11G_CD.xml)

The cross sections are then loaded and applied to neutronics materials, namely CoupledFeedbackNeutronicsMaterial. Auxiliary variables are also added that represent the grid variables.

Listing 9: Neutronics materials

[AuxVariables]
  [Tfuel]
    # fuel temperature
    initial_condition = 1000
  []
  [Tmod]
    # moderator + monolith + HP + reflector temperature
    initial_condition = 1000
  []
  [CD]
    # drum angle (0 = fully in, 180 = fully out)
  []
[]

[GlobalParams]
  library_file = empire_core_modified_11G_CD.xml
  library_name = empire_core_modified_11G_CD
  densities = 1.0
  isotopes = 'pseudo'
  dbgmat = false
  grid_names = 'Tfuel Tmod CD'
  grid_variables = 'Tfuel Tmod CD'
  is_meter = true

  # Reduces transfers efficiency for now, can be removed once transferred fields are checked
  bbox_factor = 10
[]

[Materials]
  [fuel]
    type = CoupledFeedbackNeutronicsMaterial
    block = '1 2' # fuel pin with 1 cm outer radius, no gap
    material_id = 1001
    plus = true
  []
  [moderator]
    type = CoupledFeedbackNeutronicsMaterial
    block = '3 4 5' # moderator pin with 0.975 cm outer radius
    material_id = 1002
  []
  [monolith]
    type = CoupledFeedbackNeutronicsMaterial
    block = '8'
    material_id = 1003
  []
  [hpipe]
    type = CoupledFeedbackNeutronicsMaterial
    block = '6 7' # gap homogenized with HP
    material_id = 1004
  []
  [be]
    type = CoupledFeedbackNeutronicsMaterial
    block = '10 11 14 15'
    material_id = 1005
  []
  [drum]
    type = CoupledFeedbackRoddedNeutronicsMaterial
    block = '13'
    front_position_function = drum_fun
    rotation_center = '${x_center} ${y_center} 0'
    rod_segment_length = '270 90'
    segment_material_ids = '1005 1006'
    isotopes = 'pseudo; pseudo'
    densities = '1.0 1.0'
  []
  [air]
    type = CoupledFeedbackNeutronicsMaterial
    block = '20 21 22'
    material_id = 1007
  []
[]
(microreactors/drum_rotation/neutronics_eigenvalue.i)

Note the exception for the material type in the drum region. Here, we use CoupledFeedbackRoddedNeutronicsMaterial, which is material specifically for dealing with control drum and rod movement. The position of the drum is controlled by the front_position_function parameter, which accepts the name of a MOOSE Function dependent on time. Typically, three functions are defined: (1) the offset between what the actual position and what the user defines as the position, (2) a user-defined position as a function of time, and (3) combining the offset and position. For the steady-state input, is function is constant:

Listing 10: Defining control drum position

[Functions]
  [offset]
    type = ConstantFunction
    value = 225
  []
  [drum_position]
    type = ConstantFunction
    value = 90
  []
  [drum_fun]
    type = ParsedFunction
    expression = 'drum_position + offset'
    symbol_names = 'drum_position offset'
    symbol_values = 'drum_position offset'
  []
[]

[AuxKernels]
  [CD_aux]
    type = FunctionAux
    variable = CD
    function = drum_position
    execute_on = 'initial timestep_end'
  []
[]
(microreactors/drum_rotation/neutronics_eigenvalue.i)

It is often useful debug the offset by creating a simplified input that has drum position outputted to exodus. This can be done using the following debug option:

Listing 11: Debug option for outputting drum position

[Debug]
  show_rodded_materials_average_segment_in = segment_id
[]
(microreactors/drum_rotation/neutronics_eigenvalue.i)

Thermal Properties

The following lists the thermal properties used for this model:

Table 2: Thermal properties

RegionDensity (kg/m)Conductivity (W/m-K)Specific Heat Capacity (J/kg-K)
Fuel14,30019300
Moderator4,30020300
Monolith1,800701,830
Reflector1,853.762001,825
Absorber2,520201,000

These properties are applied using HeatConductionMaterial and Density:

Listing 12: Thermal materials

[Materials]
  #### DENSITY #####
  # units of kg/m^3
  [fuel_density]
    type = Density
    block = '${fuel_blocks}'
    density = 14.3e3 # same as in Serpent input
  []
  [moderator_density]
    type = Density
    block = '${moderator_blocks}'
    density = 4.3e3 # same as in Serpent input
  []
  [monolith_density]
    type = Density
    block = '${monolith_blocks}'
    density = 1.8e3 # same as in Serpent input
  []
  [reflector_density]
    type = Density
    block = '${reflector_blocks}'
    density = 1.85376e3 # same as in Serpent input
  []

  ### THERMAL CONDUCTIVITY ###
  # units of W/m-K
  [fuel_kappa]
    type = HeatConductionMaterial
    block = '${fuel_blocks}'
    # temp = Tsolid
    thermal_conductivity = 19 # W/m/K
    specific_heat = 300 # reasonable value
  []
  [moderator_kappa]
    type = HeatConductionMaterial
    block = '${moderator_blocks}'
    # temp = Tsolid
    thermal_conductivity = 20 # W/m/K
    specific_heat = 500 # random value
  []
  [monolith_thermal_props]
    type = HeatConductionMaterial
    block = '${monolith_blocks}'
    # temp = Tsolid
    thermal_conductivity = 70 # typical value for G348 graphite
    specific_heat = 1830 # typical value for G348 graphite
  []
  [reflector_kappa]
    type = HeatConductionMaterial
    block = '${reflector_blocks}'
    # temp = Tsolid
    thermal_conductivity = 200 # W/m/K, typical value for Be
    specific_heat = 1825 # typical value for Be
  []

  [drum_material]
    type = GenericFunctionMaterial
    block = '${drum_blocks}'
    prop_names = 'thermal_conductivity    specific_heat density   thermal_conductivity_dT'
    prop_values = 'drum_k                  drum_cp        drum_rho 0'
  []
[]
(microreactors/drum_rotation/thermal_ss.i)

Since the materials in the drum region change based on the position of the control drum, GenericFunctionMaterial is used to manually set these properties. Based on the position (provided by a post-processor values), these functions define whether a supplied (x,y) location is outside or inside the absorber region and output the appropriate property value:

Listing 13: Functions for thermal properties in drum region

[Functions]
  [cos_theta_hat]
    type = ParsedFunction
    expression = '(cos(theta*pi/180) * (xc - x) + sin(theta*pi/180) * (yc - y)) / (sqrt((xc-x)^2+(yc-y)^2))'
    symbol_names = 'theta xc yc'
    symbol_values = 'drum_position ${x_center} ${y_center}'
  []
  [drum_k]
    type = ParsedFunction
    expression = 'if(cost < cos45, 200, 20)'
    symbol_names = 'cost cos45'
    symbol_values = 'cos_theta_hat 0.7071067811865476'
  []
  [drum_cp]
    type = ParsedFunction
    expression = 'if(cost < cos45, 1825, 1000)'
    symbol_names = 'cost cos45'
    symbol_values = 'cos_theta_hat 0.7071067811865476'
  []
  [drum_rho]
    type = ParsedFunction
    expression = 'if(cost < cos45, 1.85376e3, 2.52e3)'
    symbol_names = 'cost cos45'
    symbol_values = 'cos_theta_hat 0.7071067811865476'
  []
[]
(microreactors/drum_rotation/thermal_ss.i)

Physics

Neutronics Physics

For this model, DFEM-SN discretization is used for the neutronics. The angular quadrature uses Gauss-Chebyshev collocation with one polar angle and three azimuthal. There are six delayed neutron precursors (evident from multigroup cross-section library) and first-order scattering. The physics for the neutronics is setup using the TransportSystems syntax in Griffin:

Listing 14: Defining the physics in the neutronics input

[TransportSystems]
  equation_type = eigenvalue
  particle = neutron
  G = 11

  ReflectingBoundary = 'bottom topleft'
  VacuumBoundary = 'right'

  [sn]
    scheme = DFEM-SN
    n_delay_groups = 6
    family = MONOMIAL
    order = FIRST
    AQtype = Gauss-Chebyshev
    NPolar = 1 # use >=2 for final runs (4 sawtooth nodes sufficient)
    NAzmthl = 3 # use >=6 for final runs (4 sawtooth nodes sufficient)
    NA = 1
    sweep_type = asynchronous_parallel_sweeper
    using_array_variable = true
    collapse_scattering = true
    hide_angular_flux = true
    hide_higher_flux_moment = 0
  []
[]
(microreactors/drum_rotation/neutronics_eigenvalue.i)

The initial power for this model is 2 MWth. Due to the 1/12th 2D geometry, this power equations to 83 kW/m in the simulation. The PowerDensity syntax in griffin adds an auxiliary variable for the computed power density field:

Listing 15: Power density evaluation

[PowerDensity]
  power = '${fparse 2e6 / 12 / 2}' # power: 2e6 W / 12 / 2 m (axial)
  power_density_variable = power_density
  integrated_power_postprocessor = integrated_power
[]
(microreactors/drum_rotation/neutronics_eigenvalue.i)

Thermal Physics

The heat conduction problem is only solved on the solid portions of the domain, i.e. excluding heat pipes, air, and moderator gaps. As such the system variable is restricted to these blocks, which requires removing some default checks in the problem:

Listing 16: Block restricted variable in thermal input

[Variables]
  [Tsolid]
    initial_condition = 950
    block = '${solid_blocks}'
  []
[]

[Problem]
  kernel_coverage_check = false
  material_coverage_check = false
[]
(microreactors/drum_rotation/thermal_ss.i)

The volumetric kernels are quite simple, including heat conduction in the solid regions and a source applied in the fuel region from the power density:

Listing 17: Thermal input volumetric kernels

[Kernels]
  [heat_conduction]
    type = HeatConduction
    variable = Tsolid
    block = '${solid_blocks}'
  []
  [heat_source_fuel]
    type = CoupledForce
    variable = Tsolid
    block = '${fuel_blocks}'
    v = power_density
  []
[]
(microreactors/drum_rotation/thermal_ss.i)

In order to capture the heat transfer into the heat pipes and between the moderator gaps, a few side sets are added to the mesh:

Listing 18: Side sets added to thermal mesh for boundary and gap heat transfer

[Mesh]
  [main]
    type = FileMeshGenerator
    file = empire_2d_CD_fine_in.e
  []
  [add_sideset_hp]
    type = SideSetsBetweenSubdomainsGenerator
    input = main
    primary_block = '8' # add 16 so the HP boundary extends into the upper axial reflector
    paired_block = '7'
    new_boundary = 'hp'
  []
  [add_sideset_inner_mod_gap]
    type = SideSetsBetweenSubdomainsGenerator
    input = add_sideset_hp
    primary_block = '4'
    paired_block = '5'
    new_boundary = 'gap_mod_inner'
  []
  [add_sideset_outer_mod_gap]
    type = SideSetsBetweenSubdomainsGenerator
    input = add_sideset_inner_mod_gap
    primary_block = '8'
    paired_block = '5'
    new_boundary = 'gap_mod_outer'
  []
[]
(microreactors/drum_rotation/thermal_ss.i)

The convective heat transfer into the heat pipes is applied as a boundary condition, while the moderator gaps use GapHeatTransfer:

Listing 19: Side-set heat transfer

[BCs]
  [hp_bc]
    type = CoupledConvectiveHeatFluxBC
    variable = Tsolid
    boundary = 'hp' # inside of heat pipe
    htc = htc_hp_var
    T_infinity = Tsink_hp_var # eventually, will be given by Sockeye
  []
  [rrefl_bc]
    type = CoupledConvectiveHeatFluxBC
    variable = Tsolid
    boundary = 'right' # reflector cooling
    htc = 50 # arbitrary (but small) HTC
    T_infinity = 500 # arbitrary Tsink
  []
[]

[AuxVariables]
  [power_density]
    block = '${fuel_blocks}'
    family = MONOMIAL
    order = FIRST
    initial_condition = 2.3e6
  []
  [htc_hp_var]
    block = '${monolith_blocks}'
    initial_condition = 3e2
  []
  [Tsink_hp_var]
    block = '${monolith_blocks}'
    initial_condition = 900
  []
[]

[ThermalContact]
  [RPV_gap]
    type = GapHeatTransfer
    gap_geometry_type = 'PLATE'
    emissivity_primary = 0.8
    emissivity_secondary = 0.8 # varies from 0.6 to 1.0 in RPV paper, 0.79 in NRC paper
    variable = Tsolid # graphite -> rpv gap
    primary = 'gap_mod_inner'
    secondary = 'gap_mod_outer'
    gap_conductivity = 0.08 # W/m/K, typical thermal connectivity for air at 1000C
    quadrature = true
  []
[]
(microreactors/drum_rotation/thermal_ss.i)

Coupling

The coupling between the neutronics and thermal inputs utilize the MultiApps and Transfers systems. For this model, the neutronics serves as the main application. For steady-state, the thermal sub-application is defined as a FullSolveMultiApp:

Listing 20: Steady-state multiapp block

[MultiApps]
  [bison]
    type = FullSolveMultiApp
    input_files = thermal_ss.i
    execute_on = 'timestep_end'
    no_restore = true
  []
[]
(microreactors/drum_rotation/neutronics_eigenvalue.i)

Power density and drum position are transferred to the sub-application and temperatures are retrieved. Power density simply uses a MultiAppCopyTransfer, which directly copies the degrees of freedom of the variable from the main application to the sub-application. The drum position is computed as a post-processor, where the data is transferred via MultiAppReporterTransfer. Since the cross sections are tabulated with separate temperatures for fuel and non-fuel. However, these variables must be defined everywhere in order to evaluate the cross sections. To do this, we use MultiAppGeneralFieldNearestLocationTransfer that is block restricted on the sub-application side. This transfer will then do a direct degree a freedom transfer for Tfuel in the fuel region and extrapolate from the nearest node in the non-fuel region. The same goes for Tmod in the non-fuel vs. fuel regions.

Listing 21: Variable transfers

[Transfers]
  [pdens_to_modules]
    type = MultiAppCopyTransfer
    to_multi_app = bison
    source_variable = power_density
    variable = power_density
  []
  [dp_to_modules]
    type = MultiAppReporterTransfer
    to_multi_app = bison
    from_reporters = dp/value
    to_reporters = drum_position/value
  []
  [tfuel_from_modules]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = bison
    variable = Tfuel
    source_variable = Tsolid
    from_blocks = '1 2'
    # Values are transferred outside the block restriction of Tfuel, leading to some indetermination
    search_value_conflicts = false
  []
  [tmod_from_modules]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = bison
    variable = Tmod
    source_variable = Tsolid
    from_blocks = '3 4 8 10 11 13 14 15 22'
    search_value_conflicts = false
  []
[]
(microreactors/drum_rotation/neutronics_eigenvalue.i)

Solvers

Griffin utilizes a customized executioner for DFEM-SN systems known as SweepUpdate. See this Griffin tutorial for more details on the methodology.

Listing 22: Griffin executioner for steady-state DFEM-SN with CMFD

[Mesh]
  [main]
    type = FileMeshGenerator
    file = empire_2d_CD_fine_in.e
  []
  [coarse_mesh]
    type = FileMeshGenerator
    file = empire_2d_CD_coarse_in.e
  []
  [assign_coarse]
    type = CoarseMeshExtraElementIDGenerator
    input = main
    coarse_mesh = coarse_mesh
    extra_element_id_name = coarse_element_id
  []
[]

[Executioner]
  type = SweepUpdate

  richardson_rel_tol = 1e-10
  richardson_abs_tol = 1e-8
  richardson_max_its = 100
  richardson_value = eigenvalue
  inner_solve_type = GMRes
  max_inner_its = 2

  cmfd_acceleration = true
  coarse_element_id = coarse_element_id
  prolongation_type = multiplicative
  max_diffusion_coefficient = 1
[]
(microreactors/drum_rotation/neutronics_eigenvalue.i)

The thermal application simply uses the traditional Steady executioner, which specifies using AMG to speed up convergence.

Listing 23: Executioner for steady-state thermal problem

[Executioner]
  type = Steady
  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart'
  petsc_options_value = 'hypre boomeramg 300'
  line_search = 'none'

  l_tol = 1e-02
  nl_abs_tol = 1e-7
  nl_rel_tol = 1e-8

  l_max_its = 50
  nl_max_its = 25
[]
(microreactors/drum_rotation/thermal_ss.i)

Transient

Initial Conditions

For this model, the initial condition comes from the steady-state solution. This solution is written to a file and loaded in the transient simulation using SolutionVectorFile and TransportSolutionVectorFile. The steady-state solution is written by adding the user-objects to the steady-state inputs, with the important parameter specification: execute_on = FINAL.

Listing 24: Writing steady-state neutronics solution to a file

[UserObjects]
  [neutronics_initial]
    type = TransportSolutionVectorFile
    transport_system = sn
    folder = 'binary_90'
    execute_on = final
  []
  [neutronics_thermal_initial]
    type = SolutionVectorFile
    var = 'Tfuel Tmod'
    folder = 'binary_90'
    execute_on = final
  []
[]
(microreactors/drum_rotation/neutronics_eigenvalue.i)

Listing 25: Writing steady-state thermal solution to a file

[UserObjects]
  [thermal_initial]
    type = SolutionVectorFile
    var = 'Tsolid power_density'
    folder = 'binary_90'
    execute_on = final
  []
[]
(microreactors/drum_rotation/thermal_ss.i)

The files are loaded with the same objects (must also be the same name) with parameters: writing = false and execute_on = INITIAL.

Listing 26: Loading steady-state neutronics solution from file

[UserObjects]
  [neutronics_initial]
    type = TransportSolutionVectorFile
    transport_system = sn
    folder = 'binary_90'
    writing = false
    execute_on = initial
  []
  [neutronics_adjoint]
    type = TransportSolutionVectorFile
    folder = 'binary_90'
    transport_system = sn
    writing = false
    load_to_adjoint = true
    disable_eigenvalue_transfer = true
    execute_on = initial
  []
  [neutronics_thermal_initial]
    type = SolutionVectorFile
    var = 'Tfuel Tmod'
    folder = 'binary_90'
    writing = false
    execute_on = initial
  []
[]
(microreactors/drum_rotation/neutronics_transient.i)

Listing 27: Loading steady-state thermal solution from file

[UserObjects]
  [thermal_initial]
    type = SolutionVectorFile
    var = 'Tsolid power_density'
    folder = 'binary_90'
    writing = false
    execute_on = initial
  []
[]
(microreactors/drum_rotation/thermal_transient.i)

Adjoint Problem

In order to compute point-kinetics parameters, like dynamic reactivity and mean generation time, Griffin requires the computation of an adjoint solution. This is done by running the adjoint version of the steady-state neutronics input. There adjoint input is basically the same as the forward input, with a couple key difference. First, the problem is de-coupled, so there are MultiApps or Transfers and the thermal solution is loaded from forward solution file. Second, TransportSystems must have the parameter for_adjoint = true set.

Listing 28: Neutronics adjoint problem

[UserObjects]
  [neutronics_adjoint]
    type = TransportSolutionVectorFile
    transport_system = sn
    folder = 'binary_90'
    execute_on = final
  []
  [neutronics_thermal_initial]
    type = SolutionVectorFile
    folder = 'binary_90'
    var = 'Tfuel Tmod'
    writing = false
    execute_on = initial
  []
[]

[TransportSystems]
  equation_type = eigenvalue
  particle = neutron
  G = 11

  ReflectingBoundary = 'bottom topleft'
  VacuumBoundary = 'right'

  for_adjoint = true
  [sn]
    scheme = DFEM-SN
    n_delay_groups = 6
    family = MONOMIAL
    order = FIRST
    AQtype = Gauss-Chebyshev
    NPolar = 1 # use >=2 for final runs (4 sawtooth nodes sufficient)
    NAzmthl = 3 # use >=6 for final runs (4 sawtooth nodes sufficient)
    NA = 1
    sweep_type = asynchronous_parallel_sweeper
    using_array_variable = true
    collapse_scattering = true
    hide_angular_flux = true
  []
[]
(microreactors/drum_rotation/neutronics_adjoint.i)

Neutronics Transient

The transient input for neutronics is largely identical to the steady-state version. The following presents the key differences. First, equation_type = transient in the TransportSystems block must be set. Second, function describing the control drum position now depends on time. Third, the multi-app for the thermal sub-application is now a TransientMultiApp.

Listing 29: Changes for neutronics transient input

[TransportSystems]
  equation_type = transient
  particle = neutron
  G = 11

  ReflectingBoundary = 'bottom topleft'
  VacuumBoundary = 'right'

  [sn]
    scheme = DFEM-SN
    n_delay_groups = 6
    family = MONOMIAL
    order = FIRST
    AQtype = Gauss-Chebyshev
    NPolar = 1 # use >=2 for final runs (4 sawtooth nodes sufficient)
    NAzmthl = 3 # use >=6 for final runs (4 sawtooth nodes sufficient)
    NA = 1
    dnp_amp_scheme = quadrature
    sweep_type = asynchronous_parallel_sweeper
    using_array_variable = true
    collapse_scattering = true
    hide_angular_flux = true
    hide_higher_flux_moment = 0
  []
[]

[Functions]
  [offset]
    type = ConstantFunction
    value = 225
  []
  [drum_linear]
    type = ParsedFunction
    expression = 'if(t <= t_out, pos_start + speed * t, (pos_start + speed * t_out) - speed * (t - t_out))'
    symbol_names = 't_out pos_start speed'
    symbol_values = '${t_out} ${pos_start} ${speed}'
  []
  [drum_position]
    type = ParsedFunction
    expression = 'if(drum_linear < 0, 0, if(drum_linear > 180, 180, drum_linear))'
    symbol_names = 'drum_linear'
    symbol_values = 'drum_linear'
  []
  [drum_fun]
    type = ParsedFunction
    expression = 'drum_position + offset'
    symbol_names = 'drum_position offset'
    symbol_values = 'drum_position offset'
  []
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    input_files = thermal_transient.i
    execute_on = 'timestep_end'
  []
[]
(microreactors/drum_rotation/neutronics_transient.i)

Finally, the executioner is changed to IQSSweepUpdate, which utilized the improved quasi-static method (IQS). The algorithm this executioner implements is described in detail in Prince et al. (2023). The pke_param_csv = true parameter triggers an output of point kinetics parameters and output_micro_csv = true outputs the detailed power profile of the transient.

Listing 30: Neutronics solver using IQS

[Executioner]
  type = IQSSweepUpdate
  pke_param_csv = true
  output_micro_csv = true

  end_time = 5
  dt = '${fparse angle_step / speed * dstep}'
  dtmin = 0.001

  richardson_rel_tol = 1e-4
  richardson_abs_tol = 5e-5
  richardson_max_its = 200
  richardson_value = integrated_power
  inner_solve_type = GMRes
  max_inner_its = 2

  cmfd_acceleration = true
  coarse_element_id = coarse_element_id
  prolongation_type = multiplicative
  max_diffusion_coefficient = 1
[]
(microreactors/drum_rotation/neutronics_transient.i)

Thermal Transient

Again, the transient input of the thermal model is largely the same as the steady-state input. First, the time derivative is added to the volumetric kernels. Second, the executioner type is changed to Transient. Note that the time stepping scheme is not included in this input, since it will be defined by the main neutronics application.

Listing 31: Changes for thermal transient input

[Kernels]
  [time_derivative]
    type = HeatConductionTimeDerivative
    variable = Tsolid
    block = '${solid_blocks}'
  []
  [heat_conduction]
    type = HeatConduction
    variable = Tsolid
    block = '${solid_blocks}'
  []
  [heat_source_fuel]
    type = CoupledForce
    variable = Tsolid
    block = '${fuel_blocks}'
    v = power_density
  []
[]

[Executioner]
  type = Transient

  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart'
  petsc_options_value = 'hypre boomeramg 300'
  line_search = 'none'

  l_tol = 1e-02
  nl_abs_tol = 1e-7
  nl_rel_tol = 1e-8

  l_max_its = 50
  nl_max_its = 25
[]
(microreactors/drum_rotation/thermal_transient.i)

Running Model

This model can be run using a Griffin executable, or any application built with Griffin (BlueCRAB, Direwolf, etc.). The simulations scale decently well, so either a workstation or HPC can be used. The following presents the commands used to run the model and produce results. They must be executed in this order and with the same number of processors. Run times are presented at various processor counts, which were obtained by running the simulations on the INL Sawtooth cluster.

First is building the meshes, which must be run in this order, which produces drum_in.e, empire_2d_CD_fine_in.e, and empire_2d_CD_coarse_in.e:

griffin-opt -i drum.i --mesh-only
griffin-opt -i empire_2d_CD_fine.i --mesh-only
griffin-opt -i empire_2d_CD_coarse.i --mesh-only

Second is the coupled steady-state calculation, which produces neutronics_eigenvalue_out.e and neutronics_eigenvalue_out_bison0.e:

mpiexec -n <n> griffin-opt -i neutronics_eigenvalue.i

Table 3: Run times for steady-state simulation on various processors

Processors (<n>)Run-time (s)
4572
8324
16184
32122
4894

Third is the adjoint neutronics calculation:

mpiexec -n <n> griffin-opt -i neutronics_adjoint.i

Table 4: Run times for adjoint simulation on various processors

Processors (<n>)Run-time (s)
496
857
1634
3222
4817

Finally, the coupled transient can be run, which produces several output files which can be postprocessed:

  • neutronics_transient_out.e: Contains neutronics-evaluated quantities like power density and scalar flux shape. Note that the scalar flux shape must be scaled by the power scaling factor and IQS amplitude to retrieve the actual flux.

  • neutronics_transient_out_bison0.e: Contains thermal-evaluated quantities such as the temperature field.

  • neutronics_transient_out.csv: Contains global neutronics-evaluated quantities, like power.

  • neutronics_transient_out_bison0.csv: Contains global thermal-evaluated quantities, like average and max temperatures.

  • neutronics_transient_out_pke_params.csv: Contains point-kinetics parameters, like reactivity.

  • neutronics_transient_out_micro_power.csv: Contains the flux amplitude, which can be used obtain a detailed temporal profile of the power.

mpiexec -n <n> griffin-opt -i neutronics_transient.i

Table 5: Run times for adjoint simulation on various processors

Processors (<n>)Run-time (min)
4405
8216
16121
3275
4859

Results

The initial temperature and power profile are shown in Figure 3. The full transient profile of power, average and maximum temperatures, and reactivity are shown in Figure 4. These plots show that during the first second, the transient is purely kinetics driven as reactivity from rotating the control drum is inserted into the core. Maximum power is reached around 1.25 seconds. The core temperatures rise significantly at this point, causing the power to drop due to the negative temperature feedback in the fuel, even while the drum is still rotating outward. The power drops exponentially as the drums are rotated back inward, causing the temperatures to flatten and eventually drop due to heat removal. The spatial profile of the power density and temperature during the transient are shown in Figure 5.

Figure 3: Initial power and temperature profile

Figure 4: Values of selected quantities over simulated transient for control drum rotation.

Figure 5: Power density and temperature fields over transient

References

  1. Christopher Matthews, Vincent Laboure, Mark DeHart, Joshua Hansel, David Andrs, Yaqi Wang, Javier Ortensi, and Richard C Martineau. Coupled multiphysics simulations of heat pipe microreactors using direwolf. Nuclear Technology, 207(7):1142–1162, 2021.[BibTeX]
  2. Zachary Prince, Joshua Hanophy, Vincent Labouré, Yaqi Wang, Logan Harbour, and Namjae Choi. Neutron transport methods for multiphysics heterogeneous reactor core simulation in griffin. Submitted to Annals of Nuclear Energy, 2023.[BibTeX]
  3. Sebastian Schunert, Yaqi Wang, Javier Ortensi, Vincent Laboure, Frederick Gleicher, Mark DeHart, and Richard Martineau. Control rod treatment for FEM based radiation transport methods. Annals of Nuclear Energy, 127:293–302, May 2019. doi:10.1016/j.anucene.2018.11.054.[BibTeX]
  4. Stefano Terlizzi and Vincent Labouré. Asymptotic hydrogen redistribution analysis in yttrium-hydride-moderated heat-pipe-cooled microreactors using direwolf. Annals of Nuclear Energy, 186:109735, 2023. URL: https://www.sciencedirect.com/science/article/pii/S0306454923000543, doi:https://doi.org/10.1016/j.anucene.2023.109735.[BibTeX]
  5. A. Yamamoto. A simple and efficient control rod cusping model for three-dimensional pin-by-pin core calculations. Nuclear Technology, 2004.[BibTeX]