3D Multiphysics Simulation of Control Drum Inadvertent Rotation in HP-MR

Contact: Ahmed Amin Abdelhameed (aabdelhameed.at.anl.gov), Yinbin Miao (ymiao.at.anl.gov), Nicolas Stauff (nstauff.at.anl.gov)

Model link: HPMR Model

Overview

In microreactors, modeling events involving rotating control drums can present challenges due to geometry deformation, localized power peaking, and complex multiphysics feedback. These effects must be accurately captured in high-fidelity, three-dimensional multiphysics simulations.

In this work, we investigate inadvertent drum rotation in HP-MR using a fully coupled 3D multiphysics model involving Griffin**, BISON**, and Sockeye**.

The MOOSE MultiApp system is used to couple simulations. Here Griffin is the parent application. The MultiApp system enables seamless data exchange between Griffin, its child application (BISON), and its grandchild application (Sockeye).

Figure 1: MultiApp hierarchy of the HP-MR model.

In this coupling strategy, Griffin** solves the high-fidelity neutronics problem to obtain the power density distribution at each time step of the dynamic transient. This power data is transferred to BISON** using shape-function-preserving mapping to ensure integrated power conservation, allowing BISON to solve the heat conduction equations and compute solid temperatures. The resulting heat flux at the heat pipe surfaces is passed to Sockeye**, which solves the thermal transport in the heat pipes and returns the updated surface temperatures to BISON for feedback. The simulation is tightly coupled using a Picard fixed-point iteration** with defined convergence criteria to ensure stability and accuracy.

Mesh File

In this analysis, a 3D whole-core model with 1/6 symmetry was employed. The MOOSE Reactor Module** provides several tools for generating finite element meshes, enabling rapid construction of detailed heterogeneous reactor geometry—including pins, assemblies, control drums, and peripheral zones—through operations such as extrusion and rotation.

#####################################################################
# This is the MOOSE input file to generate the mesh that can be
# used for the 1/6 core Heat-Pipe Microreactor Multiphysics
# simulations (BISON thermal simulation).
# Running this input requires MOOSE Reactor Module Objects
# Users should use
# --mesh-only HPMR_OneSixth_Core_meshgenerator_tri_rotate_bdry_fine.e
# command line argument to generate
# exodus file for further Multiphysics simulations.
#####################################################################

fine_azi_sectors = 120

[Mesh]
  # Moderator pins
  [Mod_hex]
    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 = '10'
    polygon_size = 1.15
    polygon_size_style = 'apothem'
    ring_radii = '0.825 0.92'
    ring_intervals = '2 1'
    ring_block_ids = '103 100 101' # 103 is tri mesh
    preserve_volumes = on
    quad_center_elements = false
  []
  # Heat Pipe pins
  [HP_hex]
    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 = '10'
    # inner background boundary layer, only for BISON mesh
    background_inner_boundary_layer_bias = 1.5
    background_inner_boundary_layer_intervals = 3
    background_inner_boundary_layer_width = 0.03
    polygon_size = 1.15
    polygon_size_style = 'apothem'
    ring_radii = '0.97 1.07'
    ring_intervals = '2 1'
    ring_block_ids = '203 200 201' # 203 is tri mesh
    preserve_volumes = on
    quad_center_elements = false
  []
  # Fuel pins
  [Fuel_hex]
    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 = '10'
    polygon_size = 1.15
    polygon_size_style = 'apothem'
    ring_radii = '1'
    ring_intervals = '2'
    ring_block_ids = '303 301' # 303 is tri mesh
    preserve_volumes = on
    quad_center_elements = false
  []
  # Assembly mesh with all the three types of pins
  [Patterned_o]
    type = PatternedHexMeshGenerator
    inputs = 'Fuel_hex HP_hex Mod_hex'
    hexagon_size = 13.376
    background_block_id = 10
    background_intervals = 1
    pattern = '1 0 1 0 1 0 1;
              0 2 0 2 0 2 0 0;
             1 0 1 0 1 0 1 2 1;
            0 2 0 2 0 2 0 0 0 0;
           1 0 1 0 1 0 1 2 1 2 1;
          0 2 0 2 0 2 0 0 0 0 0 0;
         1 0 1 0 1 0 1 2 1 2 1 2 1;
          0 2 0 2 0 2 0 0 0 0 0 0;
           1 0 1 0 1 0 1 2 1 2 1;
            0 2 0 2 0 2 0 0 0 0;
             1 0 1 0 1 0 1 2 1;
              0 2 0 2 0 2 0 0;
               1 0 1 0 1 0 1'
  []
  [Patterned]
    type = PatternedHexPeripheralModifier
    input = Patterned_o
    input_mesh_external_boundary = 10000
    new_num_sector = ${fine_azi_sectors}
    num_layers = 1
  []
  # Control drum at 12 o'clock
  [cd_12]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Patterned Patterned'
    sides_to_adapt = '3 4'
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = 504
    ring_radii = '12.25 13.25'
    ring_intervals = '2 1'
    ring_block_ids = '500 501 5012'
    preserve_volumes = true
    is_control_drum = true
  []
  # Define the absorber section
  # Control drum at 6 o'clock
  [cd_6]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Patterned Patterned'
    sides_to_adapt = '0 1'
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = 504
    ring_radii = '12.25 13.25'
    ring_intervals = '2 1'
    ring_block_ids = '500 501 5006'
    preserve_volumes = true
    is_control_drum = true
  []
  # Define the absorber section
  # Control drum at 9 o'clock
  [cd_9]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Patterned Patterned Patterned'
    sides_to_adapt = '0 4 5'
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = 504
    ring_radii = '12.25 13.25'
    ring_intervals = '2 1'
    ring_block_ids = '500 501 5009'
    preserve_volumes = true
    is_control_drum = true
  []
  # Define the absorber section
  # Control drum at 3 o'clock
  [cd_3]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Patterned Patterned Patterned'
    sides_to_adapt = '1 2 3'
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = 504
    ring_radii = '12.25 13.25'
    ring_intervals = '2 1'
    ring_block_ids = '500 501 5003'
    preserve_volumes = true
    is_control_drum = true
  []
  # Define the absorber section
  # Control drum at 1 o'clock
  [cd_1]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Patterned Patterned Patterned'
    sides_to_adapt = '2 3 4'
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = 504
    ring_radii = '12.25 13.25'
    ring_intervals = '2 1'
    ring_block_ids = '500 501 5001'
    preserve_volumes = true
    is_control_drum = true
  []
  # Define the absorber section
  # Control drum at 2 o'clock
  [cd_2]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Patterned Patterned'
    sides_to_adapt = '2 3'
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = 504
    ring_radii = '12.25 13.25'
    ring_intervals = '2 1'
    ring_block_ids = '500 501 5002'
    preserve_volumes = true
    is_control_drum = true
  []
  # Define the absorber section
  # Control drum at 4 o'clock
  [cd_4]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Patterned Patterned'
    sides_to_adapt = '1 2'
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = 504
    ring_radii = '12.25 13.25'
    ring_intervals = '2 1'
    ring_block_ids = '500 501 5004'
    preserve_volumes = true
    is_control_drum = true
  []
  # Define the absorber section
  # Control drum at 5 o'clock
  [cd_5]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Patterned Patterned Patterned'
    sides_to_adapt = '0 1 2'
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = 504
    ring_radii = '12.25 13.25'
    ring_intervals = '2 1'
    ring_block_ids = '500 501 5005'
    preserve_volumes = true
    is_control_drum = true
  []
  # Define the absorber section
  # Control drum at 7 o'clock
  [cd_7]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Patterned Patterned Patterned'
    sides_to_adapt = '0 1 5'
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = 504
    ring_radii = '12.25 13.25'
    ring_intervals = '2 1'
    ring_block_ids = '500 501 5007'
    preserve_volumes = true
    is_control_drum = true
  []
  # Define the absorber section
  # Control drum at 8 o'clock
  [cd_8]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Patterned Patterned'
    sides_to_adapt = '0 5'
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = 504
    ring_radii = '12.25 13.25'
    ring_intervals = '2 1'
    ring_block_ids = '500 501 5008'
    preserve_volumes = true
    is_control_drum = true
  []
  # Define the absorber section
  # Control drum at 10 o'clock
  [cd_10]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Patterned Patterned'
    sides_to_adapt = '4 5'
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = 504
    ring_radii = '12.25 13.25'
    ring_intervals = '2 1'
    ring_block_ids = '500 501 5010'
    preserve_volumes = true
    is_control_drum = true
  []
  # Define the absorber section
  # Control drum at 11 o'clock
  [cd_11]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Patterned Patterned Patterned'
    sides_to_adapt = '3 4 5'
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = 504
    ring_radii = '12.25 13.25'
    ring_intervals = '2 1'
    ring_block_ids = '500 501 5011'
    preserve_volumes = true
    is_control_drum = true
  []
  [ref_0]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Patterned'
    sides_to_adapt = '0'
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = '400 401'
  []
  [ref_1]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Patterned'
    sides_to_adapt = '1'
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = '400 401'
  []
  [ref_2]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Patterned'
    sides_to_adapt = '2'
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = '400 401'
  []
  [ref_3]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Patterned'
    sides_to_adapt = '3'
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = '400 401'
  []
  [ref_4]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Patterned'
    sides_to_adapt = '4'
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = '400 401'
  []
  [ref_5]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    meshes_to_adapt_to = 'Patterned'
    sides_to_adapt = '5'
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = '400 401'
  []
  # Central void region
  [air_center]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    meshes_to_adapt_to = 'Patterned Patterned Patterned Patterned Patterned Patterned'
    sides_to_adapt = '0 1 2 3 4 5'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = '600 601'
  []
  # Dummy assembly to help form a regular hexagonal core pattern
  # This will be deleted later
  [dummy]
    type = HexagonConcentricCircleAdaptiveBoundaryMeshGenerator
    num_sectors_per_side = '${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors} ${fine_azi_sectors}'
    hexagon_size = 13.376
    background_intervals = 2
    background_block_ids = '700 701'
    # external_boundary_id = 9998
  []
  # Stitching assemblies together to form the core mesh
  [core]
    type = PatternedHexMeshGenerator
    inputs = 'Patterned cd_1 cd_2 cd_3 cd_4 cd_5 cd_6 cd_7 cd_8 cd_9 cd_10 cd_11 cd_12 ref_0 ref_1 ref_2 ref_3 ref_4 ref_5 dummy air_center'
    # Pattern ID  0       1    2    3    4    5    6    7    8    9    10    11    12    13    14    15    16    17    18   19       20
    pattern_boundary = none
    generate_core_metadata = true
    pattern = '19 17 12 16 19;
             17 11  0  0  1 16;
           10  0  0  0  0  0  2;
          18 0  0   0  0  0  0 15;
        19  9  0  0  20  0  0  3 19;
          18 0  0   0  0  0  0 15;
            8  0  0  0  0  0  4;
             13  7  0  0  5 14;
               19 13  6 14 19'
    rotate_angle = 60
  []
  # Delete the dummy assemblies
  [del_dummy]
    type = BlockDeletionGenerator
    block = '700 701'
    input = core
    new_boundary = 10000
  []
  # Add peripheral circular region
  [add_outer_shield]
    type = PeripheralRingMeshGenerator
    input = del_dummy
    # Use `peripheral_layer_num = 2` for BISON mesh
    peripheral_layer_num = 2
    peripheral_ring_radius = 115.0
    input_mesh_external_boundary = 10000
    peripheral_ring_block_id = 250
    peripheral_ring_block_name = outer_shield
    # Peripheral boundary layer, only for BISON mesh
    peripheral_outer_boundary_layer_bias = 0.625
    peripheral_outer_boundary_layer_intervals = 3
    peripheral_outer_boundary_layer_width = 2
  []
  # trim the full core to 1/6 core
  [del_1]
    type = PlaneDeletionGenerator
    point = '0 0 0'
    normal = '10 17.32 0'
    input = add_outer_shield
    new_boundary = 147
  []
  [del_2]
    type = PlaneDeletionGenerator
    point = '0 0 0'
    normal = '10 -17.32 0'
    input = del_1
    new_boundary = 147
  []
  # extrusion
   [extrude]
     type = AdvancedExtruderGenerator
     input = del_2
     heights = '20 160 20'
     # Use `num_layers = '6 16 6'` for BISON mesh
     num_layers = '1 8 1'
     subdomain_swaps = '101  1000  100  1000  103  1003    201  201     200  200     203  203    250  250  301  1000   303  1003  400    400  401     401    10  1000   600  600    601   601    504    504     500   500   501   501  10000 1009 ;
                        101  101   100  100   103  103     201  201     200  200     203  203    250  250  301  301    303  303   400    400  401     401    10  10     600  600    601   601    504    504     500   500   501   501  10000 10000;
                        101  1000  100  1000  103  1003    201  1000    200  1000    203  1003   250  250  301  1000   303  1003  400    400  401     401    10  1000   600  600    601   601    504    504     500   500   501   501  10000 1008 '

    #subdomain_swaps = '101  1000  100  1000  103  1003    201  201     200  200     203  203    250  250  301  1000   303  1003  400    400  401     401    10  1000  503    503  600  600    601   601    504    504     500   500   501   501  10000 1009  5001   5001  5002   5002  5003   5003   5004   5004   5005  5005   5006  5006   5007  5007   5008  5008   5009  5009   5010  5010    5011  5011    5012  5012;
    #                   101  101   100  100   103  103     201  201     200  200     203  203    250  250  301  301    303  303   400    400  401     401    10  10    503    503  600  600    601   601    504    504     500   500   501   501  10000 10000 5001   5001  5002   5002  5003   5003   5004   5004   5005  5005   5006  5006   5007  5007   5008  5008   5009  5009   5010  5010    5011  5011    5012  5012;
    #                   101  1000  100  1000  103  1003    201  1000    200  1000    203  1003   250  250  301  1000   303  1003  400    400  401     401    10  1000  503    503  600  600    601   601    504    504     500   500   501   501  10000 1008  5001   5001  5002   5002  5003   5003   5004   5004   5005  5005   5006  5006   5007  5007   5008  5008   5009  5009   5010  5010    5011  5011    5012  5012'

     # biased upper and lower reflector mesh, only for BISON mesh
     #biases = '1.6 1.0 0.625'
     direction = '0 0 1'
     top_boundary = 2000
     bottom_boundary = 3000
   []
  ## Define some special reflector blocks
  #[reflector_bottom_quad]
  #  type = ParsedSubdomainMeshGenerator
  #  input = extrude
  #  combinatorial_geometry = 'z<=20'
  #  block_id = 1000
  #  excluded_subdomain_ids = '103 203 303 250 400 401 500 501 502 503 504 600 601'
  #[]
  #[reflector_bottom_tri]
  #  type = ParsedSubdomainMeshGenerator
  #  input = reflector_bottom_quad
  #  combinatorial_geometry = 'z<=20'
  #  block_id = 1003
  #  excluded_subdomain_ids = '1000 250 400 401 500 501 502 503 504 600 601'
  #[]
  #[reflector_top_quad]
  #  type = ParsedSubdomainMeshGenerator
  #  input = reflector_bottom_tri
  #  combinatorial_geometry = 'z>=180'
  #  block_id = 1000
  #  excluded_subdomain_ids = '103 203 303 200 201 250 400 401 500 501 502 503 504 600 601'
  #[]
  #[reflector_top_tri]
  #  type = ParsedSubdomainMeshGenerator
  #  input = reflector_top_quad
  #  combinatorial_geometry = 'z>=180'
  #  block_id = 1003
  #  excluded_subdomain_ids = '1000 200 201 203 250 400 401 500 501 502 503 504 600 601'
  #[]
  ## Assgin block names
  #[rename_blocks]
  #  type = RenameBlockGenerator
  #   old_block_id = '     101       100         103          201        200           203            301       303        400            401         10     600         601        504          500                501          502                  1003'
  #   new_block_name = ' mod_ss moderator_quad moderator_tri hp_ss  heat_pipes_quad heat_pipes_tri fuel_quad fuel_tri reflector_tri reflector_quad monolith  air_gap_tri air_gap_quad air_gap_quad reflector_tri  reflector_quad Drum_channel     reflector'
  #  input = extrude
  #[]
 # Assign boundary names
  [rename_boundaries]
    type = RenameBoundaryGenerator
    input = extrude
    old_boundary = '10000 2000 3000 147'
    new_boundary = 'side top bottom side_mirror'
  []
  # Scale to m unit for general MOOSE application
  # Note that is_meter needs to be true in Griffin
  [scale]
    type = TransformGenerator
    input = rename_boundaries
    transform = SCALE
    vector_value = '1e-2 1e-2 1e-2'
  []
  # Rotate to the desired orientation for simulation
  [rotate_end]
    type = TransformGenerator
    input = scale
    transform = ROTATE
    vector_value = '0 0 -120'
  []
  ## A series of process to prepare the mesh for BISON simulation
  #[split_hp_ss]
  #  type = SubdomainBoundingBoxGenerator
  #  input = rotate_end
  #  block_id = 1201
  #  block_name = 'hp_ss_up'
  #  restricted_subdomains = 'hp_ss'
  #  bottom_left = '-100 -100 1.8'
  #  top_right = '100 100 3.0'
  #[]
  #[add_exterior_ht_low]
  #  type = SideSetsBetweenSubdomainsGenerator
  #  input = split_hp_ss
  #  paired_block = 'hp_ss'
  #  primary_block = 'monolith'
  #  new_boundary = 'heat_pipe_ht_surf_low'
  #[]
  #[add_exterior_ht_up]
  #  type = SideSetsBetweenSubdomainsGenerator
  #  input = add_exterior_ht_low
  #  paired_block = 'hp_ss_up'
  #  primary_block = 'reflector_quad'
  #  new_boundary = 'heat_pipe_ht_surf_up'
  #[]
  #[add_exterior_bot]
  #  type = SideSetsBetweenSubdomainsGenerator
  #  input = add_exterior_ht_up
  #  paired_block = 'heat_pipes_quad heat_pipes_tri hp_ss'
  #  primary_block = 'reflector_quad reflector_tri'
  #  new_boundary = 'heat_pipe_ht_surf_bot'
  #[]
  #[merge_hp_surf]
  #  type = RenameBoundaryGenerator
  #  input = add_exterior_bot
  #  old_boundary = 'heat_pipe_ht_surf_low heat_pipe_ht_surf_up'
  #  new_boundary = 'heat_pipe_ht_surf heat_pipe_ht_surf'
  #[]
  #[remove_hp]
  #  type = BlockDeletionGenerator
  #  input = merge_hp_surf
  #  block = 'hp_ss hp_ss_up heat_pipes_quad heat_pipes_tri'
  #[]
  #
  ## remove extra nodesets to limit the size of the mesh
  #[clean_up]
  #  type = BoundaryDeletionGenerator
  #  input = remove_hp
  #  boundary_names = '1 3'
  #[]
[]
(microreactors/mrad/mesh/HPMR_OneSixth_finercdrum.i)

Mesh Generation Workflow

The meshing process begins by defining various pin types (moderator pins, heat pipe pins, and fuel pins) using the PolygonConcentricCircleMeshGenerator. These pins are then arranged within fuel assemblies using the PatternedHexMeshGenerator. To handle external boundaries of the hexagonal assemblies, the PatternedHexPeripheralModifier is applied. This class replaces the outermost layer of quadrilateral elements with a transitional layer of triangular elements. This ensures node placement at designated boundary positions, facilitating the stitching of adjacent hexagonal assemblies with differing boundary node counts due to interior pin variations or azimuthal discretization.

Control Drum Meshing

Control drums were meshed using the HexagonConcentricCircleAdaptiveBoundaryMeshGenerator. For drum-containing assemblies, the num_sectors_per_side parameter was increased to ensure mesh refinement in the drum region, guaranteeing the drum's actual front face aligns with the mesh's element edges at each time step.

Figure 2: HP-MR mesh used in the analysis.

Reflector and 3D Core Construction

The radial reflector was modeled with the PeripheralRingMeshGenerator to generate the outer shield region. A dummy fuel assembly was temporarily added to complete the hexagonal core pattern and later removed using the BlockDeletionGenerator. Assemblies (including those with control drums) were patterned with the PatternedHexMeshGenerator. A 1/6 symmetry model of the HP-MR core was created by trimming the full hexagonal core using the PlaneDeletionGenerator. The resulting 2D mesh was extruded into 3D using the AdvancedExtruderGenerator, with axial layers defined via the subdomain_swaps parameter.

Griffin Neutronics Simulation

Multi-Group Cross Sections

Multi-group cross sections with 11 energy groups were generated using the Serpent 2 Monte Carlo code**, based on a parametric grid defined by Control drum rotation angle** (4 values), Fuel temperature** (5 values), and Temperature of moderator, reflector, monolith, and heat pipe** (4 values), using a total of 80 Serpent-2 simulations.

################################################################################
## NEAMS Micro-Reactor Application Driver                                     ##
## Heat Pipe Microreactor Control Drum Rotation Transient                     ##
## Griffin Main Application input file                                        ##
## DFEM-SN (1, 3) with CMFD acceleration                                      ##
################################################################################

fuel_blocks = '301 303'
air_blocks = '504 600 601 250'
ring_blocks = '5010 5011 5012'
mono_blocks = '10 10000'
mod_clad_blocks = '101'
yh_blocks = '100 103'
mod_blocks = '${yh_blocks} ${mod_clad_blocks}'
ref_blocks = '1000 1008 1009 400 401 500 501'
hp_blocks = '200 201 203'

non_hp_blocks = '${fuel_blocks} ${air_blocks} ${ring_blocks} ${mono_blocks} ${mod_blocks} ${ref_blocks}'
non_hp_fuel_blocks = '${air_blocks} ${ring_blocks} ${mono_blocks} ${mod_blocks} ${ref_blocks}'
non_hp_yh_blocks = '${fuel_blocks} ${air_blocks} ${ring_blocks} ${mono_blocks} ${mod_clad_blocks} ${ref_blocks}'

angle_step = ${fparse 360 / 720} # Angular distance between elements in drum
dstep = 1                        # Number of elements to pass each timestep
pos_start = 35                   # Starting position
t_out = 1                        # Time moving outward
speed = 5                        # Degrees per second

[Mesh]
  [fmg]
  # If you do not have a presplit mesh already, you should generate it first:
  # 1. Uncomment all the mesh blocks in the !include file
  # 2. Use the exodus file in the fmg block
  # 3. Comment the "parallel_type = distributed" line
  # 4. Use moose executable to presplit the mesh
  # Once you have the presplit mesh
  # 1. Comment all the mesh blocks except the fmg block in the !include file
  # 2. Use the cpr file in the fmg block
  # 3. Uncomment the "parallel_type = distributed" line

    type = FileMeshGenerator
    file = '../mesh/HPMR_OneSixth_finercdrum_in.e'
    # file = 'griffin_mesh.cpr'
  []
  [fmg_id]
    type = SubdomainExtraElementIDGenerator
    input = fmg
    subdomains ='         200 203 100 103 301 303 10  504 600 601 201 101 400 401 250 500 1003 501 1000 10000 1008 1009'
    extra_element_id_names = 'material_id equivalence_id'
    extra_element_ids = ' 815 815 802 802 801 801 803 820 820 820 817 816 805 805 820 805  805 805  805 803    805  805;
                          815 815 802 802 801 801 803 820 820 820 817 816 805 805 820 805  805 805  805 803    805  805'
  []
  [coarse_mesh]
    type = GeneratedMeshGenerator
    dim = 3
    nx = 10
    ny = 10
    nz = 10
    xmin = -0.1
    xmax = 1.1
    ymin = -0.1
    ymax = 1.2
    zmin = -0.0
    zmax = 2.1
  []
  [assign_coarse_id]
    type = CoarseMeshExtraElementIDGenerator
    input = fmg_id
    coarse_mesh = coarse_mesh
    extra_element_id_name = coarse_element_id
  []
  uniform_refine = 0
  parallel_type = distributed
[]

[Executioner]
  type = TransientSweepUpdate

  richardson_abs_tol = 1e-4
  richardson_rel_tol = 5e-5
  richardson_max_its = 1000

  inner_solve_type = GMRes
  max_inner_its = 100

  cmfd_acceleration = true
  coarse_element_id = coarse_element_id
  prolongation_type = multiplicative
  max_diffusion_coefficient = 1

  end_time = 5000
  dt = ${fparse angle_step / speed * dstep}
  dtmin = 0.001
[]

[TransportSystems]
  particle = neutron
  equation_type = transient

  G = 11
  VacuumBoundary = '10000 2000 3000'
  ReflectingBoundary = '147'

  [sn]
    scheme = DFEM-SN
    family = MONOMIAL
    order = first
    AQtype = Gauss-Chebyshev
    NPolar = 1
    NAzmthl = 3
    NA = 2
    sweep_type = asynchronous_parallel_sweeper
    using_array_variable = true
    collapse_scattering = true
    n_delay_groups = 6
  []
[]

[AuxVariables]
  [Rotation] # drum angle (0 = fully in, 180 = fully out)
  []
  [Tf]
    initial_condition = 873.15
    order = CONSTANT
    family = MONOMIAL
  []
  [Ts]
    initial_condition = 873.15
    order = CONSTANT
    family = MONOMIAL
  []
[]

[GlobalParams]
  library_file = '../isoxml/HP-MR_XS.xml'
  library_name = HP-MR_XS
  isotopes = 'pseudo'
  densities = 1.0
  is_meter = true
  plus = true
  dbgmat = false
  grid_names = 'Rotation Ts Tfuel'
  grid_variables = 'Rotation Ts Tf'
[]

[Functions]
  [drum10_offset]
    type = ConstantFunction
    value = 15
  []
  [drum10_position]
    type = ConstantFunction
    value = 90
  []
  [drum10_fun]
    type = ParsedFunction
    expression =    'drum10_position + drum10_offset'
    symbol_names =  'drum10_position   drum10_offset'
    symbol_values = 'drum10_position   drum10_offset'
  []

  [drum11_offset]
    type = ConstantFunction
    value = 0
  []
  [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}'
  []
  [drum11_position]
    type = ParsedFunction
    expression = 'if(drum_linear < 0, 0, if(drum_linear > 180, 180, drum_linear))'
    symbol_names = 'drum_linear'
    symbol_values = 'drum_linear'
  []
  [drum11_fun]
    type = ParsedFunction
    expression = 'drum11_position + drum11_offset'
    symbol_names = 'drum11_position  drum11_offset'
    symbol_values = 'drum11_position drum11_offset'
  []

  [drum12_offset]
    type = ConstantFunction
    value = 45
  []
  [drum12_position]
    type = ConstantFunction
    value = 0
  []
  [drum12_fun]
    type = ParsedFunction
    expression =    'drum12_position + drum12_offset'
    symbol_names =  'drum12_position   drum12_offset'
    symbol_values = 'drum12_position   drum12_offset'
  []
[]

[AuxKernels]
  [drum10_Rotation_aux]
    type = FunctionAux
    variable = Rotation
    function = drum10_position
    execute_on = 'initial timestep_end'
  []

  [drum11_Rotation_aux]
    type = FunctionAux
    variable = Rotation
    function = drum11_position
    execute_on = 'initial timestep_end'
  []

  [drum12_Rotation_aux]
    type = FunctionAux
    variable = Rotation
    function = drum12_position
    execute_on = 'initial timestep_end'
  []

  [hp_Tf]
    type = SpatialUserObjectAux
    variable = Tf
    block = ${hp_blocks}
    user_object = Tf_avg
  []
  [hp_Ts]
    type = SpatialUserObjectAux
    variable = Ts
    block = ${hp_blocks}
    user_object = Ts_avg
  []
[]

[UserObjects]
  [Tf_avg]
    type = LayeredAverage
    variable = Tf
    direction = z
    num_layers = 100
    block = ${non_hp_fuel_blocks}
    execute_on = 'LINEAR TIMESTEP_END'
  []
  [Ts_avg]
    type = LayeredAverage
    variable = Ts
    direction = z
    num_layers = 100
    block = ${non_hp_yh_blocks}
    execute_on = 'LINEAR TIMESTEP_END'
  []
[]

[Materials]
  [mod]
    type = CoupledFeedbackMatIDNeutronicsMaterial
    block = '200 203 100 103 301 303 10  504 600 601 201 101 400 401 250 500 1003 501 1000 10000 1008 1009'
  []
  [drum_10]
    type = CoupledFeedbackRoddedNeutronicsMaterial
    block = '5010'
    front_position_function = drum10_fun
    rotation_center = '0 0.926716 0'
    rod_segment_length = '270 90'
    segment_material_ids = '805 811'
    isotopes = 'pseudo; pseudo'
    densities = '1.0 1.0'
    mesh_alignment_tolerance=1E-4
  []
   [drum_11]
    type = CoupledFeedbackRoddedNeutronicsMaterial
    block = '5011'
    front_position_function = drum11_fun
    rotation_center = '0.40128 0.695037 0'
    rod_segment_length = '270 90'
    segment_material_ids = '805 811'
    isotopes = 'pseudo; pseudo'
    densities = '1.0 1.0'
    mesh_alignment_tolerance=1E-4
  []
    [drum_12]
    type = CoupledFeedbackRoddedNeutronicsMaterial
    block = '5012'
    front_position_function = drum12_fun
    rotation_center = '0.80256 0.463358 0'
    rod_segment_length = '270 90'
    segment_material_ids = '805 811'
    isotopes = 'pseudo; pseudo'
    densities = '1.0 1.0'
    mesh_alignment_tolerance=1E-4
  []
[]

[PowerDensity]
  power = 345.6e3
  power_density_variable = power_density
  integrated_power_postprocessor = integrated_power
[]

[MultiApps]
  [bison]
    type = TransientMultiApp
    app_type = BisonApp
    input_files = HPMR_thermo_tr.i
    execute_on = 'initial timestep_end'
  []
[]

[Transfers]
  [to_sub_power_density]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    to_multi_app = bison
    source_variable = power_density
    variable = power_density
    from_postprocessors_to_be_preserved = integrated_power
    to_postprocessors_to_be_preserved = power
    execute_on = 'timestep_end'
  []
  [from_sub_temp_fuel]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = bison
    variable = Tf
    source_variable = Tfuel
    execute_on = 'initial timestep_end'
    to_blocks = ${non_hp_blocks}
  []
  [from_sub_temp_mod]
    type = MultiAppGeneralFieldNearestLocationTransfer
    from_multi_app = bison
    variable = Ts
    source_variable = Tmod
    execute_on = 'initial timestep_end'
    to_blocks = ${non_hp_blocks}
  []
[]

[UserObjects]
  [ss]
    type = TransportSolutionVectorFile
    transport_system = sn
    writing = false
    execute_on = initial
    folder = '../3D_core_drum_rotation_ss'
  []
[]

[Postprocessors]
  [scaled_power_avg]
    type = ElementAverageValue
    block = ${fuel_blocks}
    variable = power_density
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_avg]
    type = ElementAverageValue
    variable = Tf
    block = ${fuel_blocks}
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_max]
    type = ElementExtremeValue
    value_type = max
    variable = Tf
    block = ${fuel_blocks}
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_min]
    type = ElementExtremeValue
    value_type = min
    variable = Tf
    block = ${fuel_blocks}
    execute_on = 'initial timestep_end'
  []
  [mod_temp_avg]
    type = ElementAverageValue
    variable = Ts
    block = ${yh_blocks}
    execute_on = 'initial timestep_end'
  []
  [mod_temp_max]
    type = ElementExtremeValue
    value_type = max
    variable = Ts
    block = ${yh_blocks}
    execute_on = 'initial timestep_end'
  []
  [mod_temp_min]
    type = ElementExtremeValue
    value_type = min
    variable = Ts
    block = ${yh_blocks}
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  csv = true
  exodus = false
  perf_graph = true
  checkpoint = false
[]
(microreactors/mrad/3D_core_drum_rotation_tr/HPMR_dfem_griffin_tr.i)

Mesh Split

A split-mesh workflow in Griffin entails first generating the appropriate mesh split configuration before running simulations with MOOSE. This process is guided by the [Mesh] block in Griffin’s input files. If a presplit mesh is not already available, the user must first generate it by uncommenting all mesh blocks, using the Exodus mesh file in the [fmg] block, and commenting out the parallel_type = distributed. Then, the MOOSE executable is used to perform the mesh splitting. Once the pre-split mesh is generated, the user switches to simulation mode by commenting all mesh blocks except the [fmg] block, replacing the Exodus file with the pre-split .cpr file, and uncommenting the parallel_type = distributed line to enable distributed processing during simulation.

Mesh Configuration and Element ID Assignment

In the [Mesh] block, the SubdomainExtraElementIDGenerator was used to assign extra element IDs that were derived from the mesh’s subdomain IDs. The subdomains parameter in this case used numeric identifiers. The extra_element_ids input was specified as a two-dimensional vector, with each row corresponding to an entry in the extra_element_id_names parameter.

A coarse mesh was generated using the GeneratedMeshGenerator, with the mesh dimension set to three and the number of elements in the x, y, and z directions set to 10. The mesh boundaries were defined using the xmax, xmin, ymax, ymin, zmax, and zmin parameters. To associate coarse element IDs with the fine mesh, the CoarseMeshExtraElementIDGenerator was applied. This generator used the coarse_mesh parameter to specify the coarse mesh source, and the extra_element_id_name parameter was used to define the name of the element ID field that was added to the fine mesh. This setup ensured that each element in the fine mesh received a coarse ID based on the spatial mapping defined by the coarse mesh.

Griffin Solver

The Griffin simulations used the DFEM-SN transport solver** with CMFD acceleration**. Angular discretization was performed using a Gauss-Chebyshev quadrature. In the multiphysics analysis, the angular quadrature was configured with one polar angle and three azimuthal angles, which produced a total of 24 discrete directions in the three-dimensional model. The finite element shape function family for the primal variables, i.e., the angular fluxes, was set to MONOMIAL, using first-order shape functions. The maximum order of scattering anisotropy was set to 2.

To improve computational efficiency, the using_array_variable setting was enabled. This allowed angular fluxes for each group to be stored using MOOSE’s ArrayVariable system, which reduces the number of computational kernels and lowers the overall cost. In addition, the collapse_scattering option was enabled to allow in-group scattering sources to be formed directly within the scattering kernels, further decreasing simulation runtime.

CMFD acceleration was enabled to assist with convergence. The prolongation_type was set to multiplicative, meaning that the coarse mesh solution was used to update the fine mesh flux in a multiplicative manner. This choice affects convergence behavior and may be tuned based on solver performance. A maximum diffusion coefficient value of 1 was defined using the max_diffusion_coefficient parameter. When this cap is triggered, a warning is issued, and the value is clipped. Users must ensure that the Richardson iteration converges properly when this cap is applied, as limiting the diffusion coefficient can have a negative impact on convergence.

Solver Configuration

The Richardson outer iteration was configured with an absolute tolerance of 1e-4, a relative tolerance of 5e-5, and a maximum of 1000 iterations. The inner linear solver used GMRES, with a maximum of 100 iterations per solve.

Neutronics Material Models

For all material regions except the control drums, the CoupledFeedbackMatIDNeutronicsMaterial was used. This material model performs on-the-fly interpolation of cross sections that are tabulated in the ISOXML XML format. The tabulation is evaluated dynamically using the local values of coupled variables governing the feedback mechanisms, such as temperature.

For control drum regions, the CoupledFeedbackRoddedNeutronicsMaterial was employed. This material type is specifically designed to handle control drum and control rod movements. Drum rotation is modeled using a set of Griffin input parameters including rod_segment_length, front_position_function, rotation_center, and segment_material_ids. These inputs define the poison and non-poison regions and their associated materials. The drum’s rotational position is governed by the front_position_function, which references a specific MOOSE Function that varies with the simulation time.

Heat Transfer Modeling

Heat conduction in solid components of the HP-MR system was modeled using BISON, which employed the HeatConduction object from MOOSE. This simulation included all solid materials except for the heat pipes. Sockeye was used to model the heat pipes, employing a 2D R-Z axisymmetric conduction model to simulate transient temperature profiles.

################################################################################
## NEAMS Micro-Reactor Application Driver                                     ##
## Heat Pipe Microreactor Control Drum Rotation Transient                     ##
## BISON Child Application input file                                         ##
## Thermal Only Physics                                                       ##
################################################################################

R_clad_o = 0.0105 # heat pipe outer radius
R_hp_hole = 0.0107 # heat pipe + gap
num_sides = 12 # number of sides of heat pipe as a result of mesh polygonization
alpha = '${fparse 2 * pi / num_sides}'
perimeter_correction = '${fparse 0.5 * alpha / sin(0.5 * alpha)}' # polygonization correction factor for perimeter
area_correction = '${fparse sqrt(alpha / sin(alpha))}' # polygonization correction factor for area
corr_factor = '${fparse R_hp_hole / R_clad_o * area_correction / perimeter_correction}'

fuel_blocks = 'fuel_tri fuel_quad'
air_blocks = 'air_gap_quad air_gap_tri outer_shield'
b4c_blocks = 'B4C'
mono_blocks = 'monolith'
mod_clad_blocks = 'mod_ss'
yh_blocks = 'moderator_quad moderator_tri'
mod_blocks = '${yh_blocks} ${mod_clad_blocks}'
ref_blocks = 'reflector_quad reflector_tri'

non_fuel_blocks = '${air_blocks} ${b4c_blocks} ${mono_blocks} ${mod_blocks} ${ref_blocks}'
non_yh_blocks = '${fuel_blocks} ${air_blocks} ${b4c_blocks} ${mono_blocks} ${mod_clad_blocks} ${ref_blocks}'

simulation_time = 10000
initial_transient_dt = 1

[GlobalParams]
  flux_conversion_factor = 1
[]

[Problem]
  restart_file_base = '../3D_core_drum_rotation_ss/HPMR_dfem_griffin_ss_out_bison0_cp_cp/LATEST'
[]

[Mesh]
  file = '../3D_core_drum_rotation_ss/HPMR_dfem_griffin_ss_out_bison0_cp_cp/LATEST'
  parallel_type = distributed
[]

[Variables]
  [temp]
    # initial_condition = 800
  []
[]

[Kernels]
  [heat_conduction_td]
    type = HeatConductionTimeDerivative
    variable = temp
  []
  [heat_conduction]
    type = HeatConduction
    variable = temp
  []
  [heat_source_fuel]
    type = CoupledForce
    variable = temp
    block = ${fuel_blocks}
    v = power_density
  []
[]

[AuxVariables]
  [power_density]
    block = ${fuel_blocks}
    family = L2_LAGRANGE
    order = FIRST
    # initial_condition = 3.4e6
  []
  [Tfuel]
    order = CONSTANT
    family = MONOMIAL
    # initial_condition = 800
  []
  [Tmod]
    order = CONSTANT
    family = MONOMIAL
    # initial_condition = 800
  []
  [fuel_thermal_conductivity]
    block = ${fuel_blocks}
    order = CONSTANT
    family = MONOMIAL
  []
  [fuel_specific_heat]
    block = ${fuel_blocks}
    order = CONSTANT
    family = MONOMIAL
  []
  [monolith_thermal_conductivity]
    block = ${mono_blocks}
    order = CONSTANT
    family = MONOMIAL
  []
  [monolith_specific_heat]
    block = ${mono_blocks}
    order = CONSTANT
    family = MONOMIAL
  []
  [flux_uo] #auxvariable to hold heat pipe surface flux from UserObject
    block = 'reflector_quad monolith'
  []
  [flux_uo_corr] #auxvariable to hold corrected flux_uo
    block = 'reflector_quad monolith'
  []
  [hp_temp_aux]
    block = 'reflector_quad monolith'
    # initial_condition = 800
  []
[]

[AuxKernels]
  [assign_tfuel_f]
    type = NormalizationAux
    variable = Tfuel
    source_variable = temp
    execute_on = 'timestep_end'
    block = ${fuel_blocks}
  []
  [assign_tfuel_nf]
    type = SpatialUserObjectAux
    variable = Tfuel
    user_object = Tf_avg
    execute_on = 'timestep_end'
    block = ${non_fuel_blocks}
  []
  [assign_tmod_m]
    type = NormalizationAux
    variable = Tmod
    source_variable = temp
    execute_on = 'timestep_end'
    block = ${yh_blocks}
  []
  [assign_tmod_nm]
    type = SpatialUserObjectAux
    variable = Tmod
    user_object = Tm_avg
    execute_on = 'timestep_end'
    block = ${non_yh_blocks}
  []

  [fuel_thermal_conductivity]
    type = MaterialRealAux
    variable = fuel_thermal_conductivity
    property = thermal_conductivity
    execute_on = timestep_end
  []
  [fuel_specific_heat]
    type = MaterialRealAux
    variable = fuel_specific_heat
    property = specific_heat
    execute_on = timestep_end
  []
  [monolith_thermal_conductivity]
    type = MaterialRealAux
    variable = monolith_thermal_conductivity
    property = thermal_conductivity
    execute_on = timestep_end
  []
  [monolith_specific_heat]
    type = MaterialRealAux
    variable = monolith_specific_heat
    property = specific_heat
    execute_on = timestep_end
  []
  [flux_uo]
    type = SpatialUserObjectAux
    variable = flux_uo
    user_object = flux_uo
  []
  [flux_uo_corr]
    type = NormalizationAux
    variable = flux_uo_corr
    source_variable = flux_uo
    normal_factor = '${fparse corr_factor}'
  []
[]

[BCs]
  [outside_bc]
    type = ConvectiveFluxFunction # (Robin BC)
    variable = temp
    boundary = 'bottom top side'
    coefficient = 1e3 # W/K/m^2
    T_infinity = 800 # K air temperature at the top of the core
  []
  [hp_temp]
    type = CoupledConvectiveHeatFluxBC
    boundary = 'heat_pipe_ht_surf'
    variable = temp
    T_infinity = hp_temp_aux
    htc = 750
  []
[]

[Materials]
  [fuel_matrix_thermal]
    type = GraphiteMatrixThermal
    block = ${fuel_blocks}
    # unirradiated_type = 'A3_27_1800'
    packing_fraction = 0.4
    specific_heat_scale_factor = 1.0
    thermal_conductivity_scale_factor = 1.0
    fast_neutron_fluence = 0 #6.75E+24 # this value is nuetron fluence over 0.1MeV
    temperature = temp
  []
  [monolith_matrix_thermal]
    type = GraphiteMatrixThermal
    block = ${mono_blocks}
    # unirradiated_type = 'A3_27_1800'
    packing_fraction = 0
    specific_heat_scale_factor = 1.0
    thermal_conductivity_scale_factor = 1.0
    fast_neutron_fluence = 0 #6.75E+24 # this value is nuetron fluence over 0.1MeV
    temperature = temp
  []
  [moderator_thermal]
    type = HeatConductionMaterial
    block = ${yh_blocks}
    temp = temp
    thermal_conductivity = 20 # W/m/K
    specific_heat = 500 # random value
  []
  [airgap_thermal]
    type = HeatConductionMaterial
    block = ${air_blocks} # Helium gap
    temp = temp
    thermal_conductivity = 0.15 # W/m/K
    specific_heat = 5197 # random value
  []
  [axial_reflector_thermal]
    type = HeatConductionMaterial
    block = ${ref_blocks}
    temp = temp
    thermal_conductivity = 199 # W/m/K
    specific_heat = 1867 # random value
  []
  [B4C_thermal]
    type = HeatConductionMaterial
    block = ${b4c_blocks}
    temp = temp
    thermal_conductivity = 92 # W/m/K
    specific_heat = 960 # random value
  []
  [SS_thermal]
    type = SS316Thermal
    temperature = temp
    block = mod_ss
  []
  [fuel_density]
    type = Density
    block = ${fuel_blocks}
    density = 2276.5
  []
  [moderator_density]
    type = Density
    block = ${yh_blocks}
    density = 4.3e3
  []
  [monolith_density]
    type = Density
    block = ${mono_blocks}
    density = 1806
  []
  [airgap_density]
    type = Density
    block = ${air_blocks} #helium
    density = 180
  []
  [axial_reflector_density]
    type = Density
    block = ${ref_blocks}
    density = 1848
  []
  [B4C_density]
    type = Density
    block = B4C
    density = 2510
  []
  [SS_density]
    type = Density
    density = 7990
    block = mod_ss
  []
[]

[MultiApps]
  [sockeye]
    type = TransientMultiApp
    app_type = SockeyeApp
    positions_file = 'hp_centers.txt'
    input_files = 'HPMR_sockeye_tr.i'
    execute_on = 'timestep_begin' # execute on timestep begin because hard to have a good initial guess on heat flux
    max_procs_per_app = 1
    output_in_position = true
    sub_cycling = true
  []
[]

[Transfers]
  # HP Transfers
  [from_sockeye_temp]
    type = MultiAppNearestNodeTransfer
    from_multi_app = sockeye
    source_variable = hp_temp_aux
    variable = hp_temp_aux
    execute_on = 'timestep_begin'
    fixed_meshes = true
  []
  [to_sockeye_flux]
    type = MultiAppNearestNodeTransfer
    variable = master_flux
    source_variable = flux_uo_corr
    to_multi_app = sockeye
    execute_on = 'timestep_begin'
    fixed_meshes = true
  []
[]

[UserObjects]
  [flux_uo]
    type = NearestPointLayeredSideDiffusiveFluxAverage
    direction = z
    num_layers = 100
    points_file = 'hp_centers.txt'
    variable = temp
    diffusivity = thermal_conductivity
    execute_on = linear
    boundary = 'heat_pipe_ht_surf'
  []
  [Tf_avg]
    type = LayeredAverage
    variable = temp
    direction = z
    num_layers = 100
    block = ${fuel_blocks}
  []
  [Tm_avg]
    type = LayeredAverage
    variable = temp
    direction = z
    num_layers = 100
    block = ${yh_blocks}
  []
[]

[Executioner]
  type = Transient

  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart '
  petsc_options_value = 'hypre boomeramg 100'

  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-8

  start_time = 0
  end_time = ${simulation_time}
  dtmin = 1
  dtmax = 5000
  # dt = ${initial_transient_dt}

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = ${initial_transient_dt}
  []
[]

[Postprocessors]
  [hp_heat_integral]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 'heat_pipe_ht_surf'
    diffusivity = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [hp_end_heat_integral]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 'heat_pipe_ht_surf_bot'
    diffusivity = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [ext_side_integral]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 'side'
    diffusivity = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [mirror_side_integral]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 'side_mirror'
    diffusivity = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [tb_integral]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 'top bottom'
    diffusivity = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [total_sink]
    type = SumPostprocessor
    values = 'hp_heat_integral hp_end_heat_integral ext_side_integral mirror_side_integral tb_integral'
    execute_on = 'initial timestep_end'
  []
  [fuel_temp_avg]
    type = ElementAverageValue
    variable = temp
    block = ${fuel_blocks}
  []
  [fuel_temp_max]
    type = ElementExtremeValue
    variable = temp
    block = ${fuel_blocks}
  []
  [fuel_temp_min]
    type = ElementExtremeValue
    variable = temp
    block = ${fuel_blocks}
    value_type = min
  []
  [mod_temp_avg]
    type = ElementAverageValue
    variable = temp
    block = ${yh_blocks}
  []
  [mod_temp_max]
    type = ElementExtremeValue
    variable = temp
    block = ${yh_blocks}
  []
  [mod_temp_min]
    type = ElementExtremeValue
    variable = temp
    block = ${yh_blocks}
    value_type = min
  []
  [monolith_temp_avg]
    type = ElementAverageValue
    variable = temp
    block = monolith
  []
  [monolith_temp_max]
    type = ElementExtremeValue
    variable = temp
    block = monolith
  []
  [monolith_temp_min]
    type = ElementExtremeValue
    variable = temp
    block = monolith
    value_type = min
  []
  [heatpipe_surface_temp_avg]
    type = SideAverageValue
    variable = temp
    boundary = heat_pipe_ht_surf
  []
  [power]
    type = ElementIntegralVariablePostprocessor
    block = ${fuel_blocks}
    variable = power_density
    execute_on = 'initial timestep_end transfer'
  []
  [fuel_cP]
    type = ElementExtremeValue
    value_type = 'max'
    variable = fuel_specific_heat
    block = ${fuel_blocks}
    execute_on = 'initial timestep_end'
  []
  [fuel_k]
    type = ElementExtremeValue
    value_type = 'max'
    variable = fuel_thermal_conductivity
    block = ${fuel_blocks}
    execute_on = 'initial timestep_end'
  []
  [monolith_cP]
    type = ElementExtremeValue
    value_type = 'max'
    variable = monolith_specific_heat
    block = ${mono_blocks}
    execute_on = 'initial timestep_end'
  []
  [monolith_k]
    type = ElementExtremeValue
    value_type = 'max'
    variable = monolith_thermal_conductivity
    block = ${mono_blocks}
    execute_on = 'initial timestep_end'
  []
  [flux_uo_avg]
    type = ElementAverageValue
    variable = flux_uo_corr
    block = 'reflector_quad monolith'
  []
[]

[Outputs]
  perf_graph = true
  color = true
  csv = false
  [exodus]
    type = Exodus
    execute_on = 'FINAL'
    enable = false
  []
  [cp]
    type = Checkpoint
    additional_execute_on = 'FINAL'
    enable = false
  []
[]
(microreactors/mrad/3D_core_drum_rotation_tr/HPMR_thermo_tr.i)
################################################################################
## NEAMS Micro-Reactor Application Driver                                     ##
## Heat Pipe Microreactor Control Drum Rotation Transient                     ##
## Sockeye Grandchild Application input file                                  ##
## Effective Heat Conduction Model with Operation Limits                      ##
################################################################################

# Total heat removed/added to heat pipe
# Q_hp = 1800.

# Wick characteristics
R_pore = 15.0e-6
D_h_pore = '${fparse 2.0 * R_pore}'
permeability = 2e-9
porosity = 0.70

# Envelope ("env")
# SS316. Incropera & DeWitt, 3rd ed, Table A.1 @ 900K (627C)
# Density (kg/m3)
rho_env = 8238.
# Thermal conductivity (W/m-K)
k_env = 23.05
# Specific heat capacity (J/kg-K)
cp_env = 589.

# Potassium vapor
# From Appendix C of text book "Heat Pipe Design and Technology"
# Sat. Vapor Potassium at T = 650C (923 K)
# Density (kg/m3)
rho_vapor = 0.193
# Effective "super" conductivity (W/m-K)
#k_vapor = 1.0e5 # Reaches 600 seconds
k_vapor = 1.0e6
# Specific heat capacity (J/kg-K)
cp_vapor = 5320.

# Potassium liquid
# From Appendix C of text book "Heat Pipe Design and Technology"
# Sat. Liquid Potassium at T = 650C (923 K)
# Density (kg/m3)
rho_liquid = 695.4
# Thermal conductivity (W/m-K)
k_liquid = 40.08
# Specific heat capacity (J/kg-K)
# From Table 1.1, no temperature data given
cp_liquid = 810.
# Melting point, Table 3.1, lists 62C, rounding up
# T_melting = 340.

# Wick, homogenize envelope and fluid
# Density (kg/m3)
rho_wick = '${fparse porosity * rho_liquid + (1.0 - porosity) * rho_env}'
# Thermal conductivity (W/m-K)
k_wick = '${fparse porosity * k_liquid + (1.0 - porosity) * k_env}'
# Specific heat capacity (J/kg-K)
# From Table 1.1, no temperature data given
cp_wick = '${fparse porosity * cp_liquid + (1.0 - porosity) * cp_env}'

# Elevations and lengths
# Note: For blackbox model -- manually update "length" input
length_evap = 180.0e-2
length_adia = 30.0e-2
length_cond = 90.0e-2

# Mesh density
# The dimensions are nicely divisible by 3 cm mesh.
nelem_base_evap = 50
nelem_base_adia = 10
nelem_base_cond = 30
mesh_density = 3
nelem_evap = '${fparse mesh_density*nelem_base_evap}'
nelem_adia = '${fparse mesh_density*nelem_base_adia}'
nelem_cond = '${fparse mesh_density*nelem_base_cond}'

# Envelope thickness
t_env = 0.08e-2
# Liquid annulus thickness
t_ann = 0.07e-2
# Wick thickness
t_wick = 0.1e-2

# Radial geometry
# Envelope outer
R_hp_o = 1.05e-2
D_hp_o = '${fparse 2.0 * R_hp_o}'
# Inner Envelope/outer annulus
R_hp_i = '${fparse R_hp_o - t_env}'
D_hp_i = '${fparse 2.0 * R_hp_i}'
# Inner annulus/wick outer
R_wick_o = '${fparse R_hp_i - t_ann}'
D_wick_o = '${fparse 2.0 * R_wick_o}'
# Inner wick/vapor core outer
R_wick_i = '${fparse R_wick_o - t_wick}'
D_wick_i = '${fparse 2.0 * R_wick_i}'

# BCs for condenser
T_ext_cond = 800.
htc_ext_cond = 1.0e6

# Evaporator parameters
# S_evap = '${fparse pi * D_hp_o * length_evap}'
# q_evap = '${fparse Q_hp / S_evap}'

simulation_time = 10000
initial_transient_dt = 1

[GlobalParams]
  scaling_factor_temperature = 1e-2
  fp_2phase = fp_2phase
[]

[FluidProperties]
  [fp_2phase]
    type = PotassiumTwoPhaseFluidProperties
    emit_on_nan = none
  []
[]

[Components]
  [hp]
    type = HeatPipeConduction

    # Common to both HeatPipe2Phase and HeatPipeBlackbox
    position = '0 0 0'
    orientation = '0 0 1'
    length = '${length_evap} ${length_adia} ${length_cond}'
    n_elems = '${nelem_evap} ${nelem_adia} ${nelem_cond}'
    gravity_vector = '0 0 -9.8'
    D_wick_i = ${D_wick_i}
    D_wick_o = ${D_wick_o}
    R_pore = ${R_pore}
    porosity = ${porosity}
    permeability = ${permeability}

    # HeatPipeConduction only
    # Axial dimensions (for heat transfer & analytic limits)
    axial_region_names = 'evap adia cond'
    L_evap = ${length_evap}
    L_adia = ${length_adia}
    L_cond = ${length_cond}
    # Radial dimensions, mesh, and materials for heat transfer problem
    D_clad_o = ${D_hp_o}
    D_clad_i = ${D_hp_i}
    D_h_pore = ${D_h_pore}
    # Mesh
    n_elems_clad = 4
    n_elems_wick = 8
    n_elems_core = 10

    # Thermal properties
    sp_vapor = sp_vapor
    sp_liquid = sp_wick_ann
    sp_wick = sp_wick_ann
    sp_clad = sp_clad
    k_core = ${k_vapor}
    k_eff = ${k_wick}

    fp_2phase = fp_2phase
    evaporator_at_start_end = true
    # Initial temperature of block
    # initial_T = ${T_ext_cond}
    T_ref = T_inner_avg

    # To evaluate the constant properties
    T_ref_density = 1000

    make_pressure_corrections = true
  []

  [condenser_boundary]
    type = HSBoundaryAmbientConvection
    boundary = 'hp:cond:outer'
    hs = hp
    T_ambient = ${T_ext_cond}
    htc_ambient = ${htc_ext_cond} #large value to approach an effective DirichletBC
    scale = bc_scale_pp
  []
  [evaporator_boundary]
    type = HSBoundaryExternalAppConvection
    boundary = 'hp:evap:outer'
    hs = hp
    T_ext = virtual_Text
    htc_ext = virtual_htc
  []
[]

[SolidProperties]
  [sp_vapor]
    type = ThermalFunctionSolidProperties
    rho = ${rho_vapor}
    cp = ${cp_vapor}
    k = ${k_vapor}
  []
  [sp_wick_ann]
    type = ThermalFunctionSolidProperties
    rho = ${rho_wick}
    cp = ${cp_wick}
    k = ${k_wick}
  []
  [sp_clad]
    type = ThermalFunctionSolidProperties
    rho = ${rho_env}
    cp = ${cp_env}
    k = ${k_env}
  []
[]

[UserObjects]
  [surf_T]
    type = LayeredSideAverage
    direction = z
    num_layers = 100
    variable = T_solid
    boundary = 'hp:evap:outer'
  []
[]

[AuxKernels]
  [hp_var]
    type = SpatialUserObjectAux
    variable = hp_temp_aux
    user_object = surf_T
  []
  [virtual_Text]
    type = ParsedAux
    variable = virtual_Text
    coupled_variables = 'T_solid master_flux virtual_htc'
    expression = 'master_flux/virtual_htc + T_solid'
  []
[]

[Functions]
  [scale_fcn]
    type = ParsedFunction
    symbol_names = 'catastrophic_pp recoverable_pp operational_pp'
    symbol_values = 'catastrophic_pp recoverable_pp operational_pp'
    expression = 'catastrophic_pp*recoverable_pp*operational_pp'
  []
[]

[AuxVariables]
  [T_wall_var]
    # initial_condition = ${T_ext_cond}
  []
  [operational_aux]
    # initial_condition = 1
  []
  [master_flux]
    # initial_condition = ${q_evap}
  []
  [hp_temp_aux]
    # initial_condition = ${T_ext_cond}
  []
  [virtual_Text]
    # initial_condition = ${T_ext_cond}
  []
  [virtual_htc]
    # initial_condition = 1.0
  []
[]

[Postprocessors]
  [Integral_BC_Total]
    type = SumPostprocessor
    values = 'condenser_boundary_integral evaporator_boundary_integral'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [ZeroPP]
    type = EmptyPostprocessor
  []
  [Integral_BC_Cond]
    type = DifferencePostprocessor
    value1 = ZeroPP
    value2 = condenser_boundary_integral
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Integral_BC_RelErr]
    type = RelativeDifferencePostprocessor
    value1 = evaporator_boundary_integral
    value2 = Integral_BC_Cond
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [bc_scale_pp]
    type = FunctionValuePostprocessor
    function = 1.0
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [operational_pp]
    type = ElementAverageValue
    variable = operational_aux
    execute_on = 'initial timestep_begin TIMESTEP_END'
  []
  # set `catastrophic_pp` as if it is recoverable for a solving to steady-state simulation
  # this MUST be changed back to catestrophic in transient simulations
  [catastrophic_pp]
    type = HeatRemovalRateLimitScale
    heat_addition_pps = 'evaporator_boundary_integral'
    limit_condenser_side = false
    recoverable_heat_removal_limit_pps = 'hp_boiling_limit hp_capillary_limit hp_entrainment_limit'
    catastrophic_heat_removal_limit_pps = ''
    T = T_inner_avg
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [recoverable_pp]
    type = HeatRemovalRateLimitScale
    heat_addition_pps = 'evaporator_boundary_integral'
    limit_condenser_side = false
    catastrophic_heat_removal_limit_pps = ''
    recoverable_heat_removal_limit_pps = 'hp_sonic_limit hp_viscous_limit'
    T = T_inner_avg
    execute_on = 'INITIAL linear nonlinear TIMESTEP_END'
  []
  [T_evap_inner]
    type = SideAverageValue
    boundary = hp:evap:inner
    variable = T_solid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_cond_inner]
    type = SideAverageValue
    boundary = hp:cond:inner
    variable = T_solid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_evap_outer]
    type = SideAverageValue
    boundary = hp:evap:outer
    variable = T_solid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_cond_outer]
    type = SideAverageValue
    boundary = hp:cond:outer
    variable = T_solid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_wall_var_avg]
    type = ElementAverageValue
    variable = T_wall_var
    execute_on = 'Initial timestep_end'
  []
  [T_inner_avg]
    type = SideAverageValue
    variable = T_solid
    boundary = hp:inner
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_inner_max]
    type = NodalExtremeValue
    variable = T_solid
    boundary = hp:inner
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_inner_min]
    type = NodalExtremeValue
    variable = T_solid
    boundary = hp:inner
    execute_on = 'INITIAL TIMESTEP_END'
    value_type = min
  []
  [DT_outer]
    type = DifferencePostprocessor
    value1 = T_evap_outer
    value2 = T_cond_outer
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [DT_inner]
    type = DifferencePostprocessor
    value1 = T_evap_inner
    value2 = T_cond_inner
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [scale_pp]
    type = FunctionValuePostprocessor
    function = scale_fcn
  []
  [A_int_master_flux]
    type = SideIntegralVariablePostprocessor
    variable = master_flux
    boundary = 'hp:evap:inner'
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [A_int_T_ext]
    type = SideIntegralVariablePostprocessor
    variable = virtual_Text
    boundary = 'hp:evap:inner'
    execute_on = 'INITIAL LINEAR'
  []
  [A_avg_T_aux]
    type = AverageNodalVariableValue
    variable = hp_temp_aux
    boundary = 'hp:evap:inner'
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[VectorPostprocessors]
  [env_vpp]
    type = NodalValueSampler
    variable = T_solid
    block = 'hp:clad'
    sort_by = x
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[Preconditioning]
  [pc]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient

  solve_type = NEWTON
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  scheme = bdf2
  line_search = none

  nl_abs_tol = 1e-6
  nl_rel_tol = 1e-8
  nl_max_its = 100

  l_tol = 1e-3
  l_max_its = 100

  start_time = 0
  end_time = ${simulation_time}
  dtmin = 1

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = ${initial_transient_dt}
  []
[]

[Outputs]
  [console]
    type = Console
    max_rows = 5
    execute_postprocessors_on = 'INITIAL TIMESTEP_END FINAL FAILED'
  []
  [csv]
    type = CSV
    execute_on = 'INITIAL TIMESTEP_END FINAL FAILED'
    execute_vector_postprocessors_on = 'FINAL FAILED'
    enable = false
  []
[]
(microreactors/mrad/3D_core_drum_rotation_tr/HPMR_sockeye_tr.i)

Neutronic Transient Results

Before proceeding to the full multiphysics analysis, it is beneficial to first examine the neutronics-only results to establish a baseline for comparison. In this study, a dynamic simulation was performed on a 1/6 core model, where the middle control drum (the one not intersected by any symmetry lines) was rotated, while the other two drums remained stationary. Due to the 1/6 core symmetry, rotating this single drum represents the synchronized rotation of six control drums in the full core. Reflective boundary conditions were applied along the symmetry planes, and vacuum boundary conditions were imposed on the top, bottom, and outer radial surfaces. At the start of the simulation, the rotating drum was inserted by 35 degrees, while the two stationary drums remained in the fully withdrawn position throughout the simulation.

The neutronics-only transient simulations were carried out by disabling the [MultiApps] and [Transfers] blocks, as well as any postprocessors related to heat transfer. The initial condition for the transient run was obtained from a previously computed steady-state solution. This steady-state solution is saved to a file and then imported into the transient simulation using the TransportSolutionVectorFile user object. To generate the steady-state solution file, appropriate UserObjects are added to the steady-state input file, with the execute_on parameter set to final to ensure the solution is captured at the end of the simulation.

The transient simulation employed the DFEM-SN method with CMFD acceleration, using an SN(2,3) with NA=2. Two scenarios were analyzed. In Scenario I, the control drums accidentally rotate outward at a rate of 5 degrees per second for 1 second, followed by an inward rotation at the same rate for 2 seconds due to reinserting the drums. In Scenario II, the drums rotate outward at a faster rate of 20 degrees per second for 1 second, then inward at the same speed for 2 seconds. As the drums move outward from their initial position, positive reactivity is introduced, causing a rise in reactor power. This increase is subsequently reduced as the drums are reinserted during the corrective action.

Figure 3: Scenario I: 3D Visualization of Power Density (Griffin-Only).

Figure 4: Scenario I: Normalized power (Griffin-Only).

Figure 5: Scenario II: Normalized power (Griffin-Only).

Multiphysics Transient Results

The control drum rotation in Griffin was modeled using a multiphysics simulation that coupled Griffin, BISON, and Sockeye, as previously described, to capture the reactor’s reactivity response to drum movement while accounting for temperature effects. Each time step in the simulation required approximately 15 Richardson iterations involving both the BISON child application and the Sockeye grandchild application, making the multiphysics control drum rotation transient computationally expensive. To reduce computational demands, the neutronics model used for the multiphysics simulation was simplified from SN(2,3) NA=2 to SN(1,3) NA=2. Compared to the neutronics-only model, the multiphysics simulation predicted a slightly lower power increase during control drum rotation. This difference is primarily attributed to temperature feedback effects, which are not captured in the single-physics (neutronics-only) model. As the control drums rotate outward, they insert positive reactivity, resulting in a rapid increase in reactor power. This power spike leads to a prompt rise in both fuel and moderator temperatures. Due to the inherently negative temperature reactivity feedback of the core, this temperature rise introduces negative reactivity, which counteracts the initial positive reactivity and ultimately suppresses the peak power. After one second of outward drum rotation, the control drums are reinserted, leading to a quick drop in both reactor power and temperatures. By the third second of the transient, both the multiphysics and neutronics-only models converge to nearly identical power levels. At this stage, the moderator and fuel temperatures are below their nominal values, and the control drums are inserted deeper than their original positions prior to the transient.

Figure 6: Multiphysics 3D visualization of HP-MR power density during the power increase period caused by control drum rotation.

Figure 7: Time evolution of HP-MR power due to control drum rotation: multiphysics simulation compared to single-physics neutronics

Figure 8: Time evolution of the average fuel and moderator temperature from multiphysics simulation due to control drum rotation.