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.