Na-HPMR Multiphysics Model

Cross-Section Generation Using Serpent Code

Similar to the K-HPMR model, the first step involves generating homogenized multi-group cross sections using Serpent-2. The Serpent-2 model for the Na-HPMR is similar to that of the K-HPMR, except that sodium is employed as the working fluid in the heat pipes instead of potassium (K). An 11-group energy structure is employed, with parameters defined over grids of:

1. Fuel temperature (5 values)

2. Temperature of the moderator, reflector, monolith, and heat pipe (5 values)

3. Hydrogen content (7 values)

4. Control drum rotation (1 value: control drums out)

The results obtained from Serpent-2 are consistent with expected reactor physics:

• Increasing the hydrogen content within the tabulated range (from YH to YH) yields a higher value, due to enhanced moderation.

• Increasing the fuel temperature consistently reduces , reflecting the negative reactivity feedback from Doppler broadening.

This multi-grid cross-section library enables analysis of thermal reactivity feedback effects as well as the impacts of hydrogen redistribution within the moderator. At present, simulations are being carried out for the control rod–out configuration.

Finally, the generated cross sections are converted into an XML-file format for compatibility with Griffin.

Sockeye Model

For each heat pipe in the 1/6 HP-MR core, a Sockeye grandchild application is used to calculate its thermal performance.

Instead of the effective thermal conductance model used in the K-HPMR model, the LCVF heat pipe model is adopted in the Na-HPMR model. As the name implies, the advanced heat pipe model uses a mechanistic flow model to predict the heat transfer performance of the heat pipe kernel.

################################################################################
## NEAMS Micro-Reactor Application Driver                                     ##
## Heat Pipe Microreactor with Na Working Fluid Steady State (Na-HPMR) SS     ##
## Sockeye Grandchild Application input file                                  ##
## Liquid-Conductance Vapor-Flow (LCVF) Heat Pipe Model (Vapor-Only)          ##
################################################################################

t_start = -5e4
t_end = 0

length_evap = 1.8
length_adia = 0.3
length_cond = 0.9

n_elems_evap = 360
n_elems_adia = 60
n_elems_cond = 180

n_elems_wick = 8
n_elems_ann  = 4
n_elems_clad = 4

D_clad_o = 2.10e-2
D_clad_i = 1.94e-2
D_wick_o = 1.80e-2
D_wick_i = 1.60e-2

porosity = 0.70
permeability = 2e-9
R_pore = 15e-6

fill_ratio = 1.1
T_ref_fill_ratio = 600

T_ref_rho_wick = 1000.0
T_ref_rho_clad = 1000.0

T_ext_cond = 900
htc_ext_cond = 1.0e6

[FluidProperties]
  [fp_2phase]
    # Better convergence compare to SodiumIdealGasTwoPhaseFluidProperties
    type = SodiumTwoPhaseFluidProperties
    emit_on_nan = none
  []
[]

[Closures]
  [closures]
    type = HeatPipeVOClosures
  []
[]

[SolidProperties]
  [sp_ss316]
    type = ThermalSS316Properties
  []
[]

[Components]
  [heat_pipe]
    type = HeatPipeVO

    position_evap_end = '0 0 0'
    direction_evap_to_cond = '0 0 1' # to be compatible with the core mesh
    gravity_vector = '9.8 0 0' # assume the core is horizontal

    L_evap = ${length_evap}
    L_adia = ${length_adia}
    L_cond = ${length_cond}

    D_clad_o = ${D_clad_o}
    D_clad_i = ${D_clad_i}
    D_wick_o = ${D_wick_o}
    D_wick_i = ${D_wick_i}

    porosity = ${porosity}
    permeability = ${permeability}
    R_pore = ${R_pore}

    n_elems_evap = ${n_elems_evap}
    n_elems_adia = ${n_elems_adia}
    n_elems_cond = ${n_elems_cond}

    n_elems_wick = ${n_elems_wick}
    n_elems_ann  = ${n_elems_ann}
    n_elems_clad = ${n_elems_clad}

    initial_T = ${T_ext_cond}

    fill_ratio = ${fill_ratio}
    T_ref_fill_ratio = ${T_ref_fill_ratio}
    T_ref_rho_wick = ${T_ref_rho_wick}
    T_ref_rho_clad = ${T_ref_rho_clad}

    fp_2phase = fp_2phase
    sp_wick = sp_ss316
    sp_clad = sp_ss316
    closures = closures

    slope_reconstruction = NONE
    stop_vapor_at_condenser_pool = true
    startup_front_option = NONE
    D_collision = 0.362e-9
    Kn_transition = 0.01
  []

  [evaporator_boundary]
    type = HSBoundaryExternalAppConvection
    boundary = 'heat_pipe:hs:evap:outer'
    hs = heat_pipe:hs
    T_ext = virtual_Text
    htc_ext = virtual_htc
  []

  # This works but is simplified
  # Improvement needed in the future
  [condenser_boundary]
    type = HSBoundaryAmbientConvection
    boundary = 'heat_pipe:hs:cond:outer'
    hs = heat_pipe:hs
    T_ambient = ${T_ext_cond}
    htc_ambient = ${htc_ext_cond}
  []
[]

[UserObjects]
  [surf_T]
    type = LayeredSideAverage
    direction = z
    num_layers = 180
    variable = T_solid
    boundary = 'heat_pipe:hs: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'
    block = '3 4 5'
  []
[]

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

[Postprocessors]
  [refill]
    type = FunctionValuePostprocessor
    function = '1'
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

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

# Added dampers as some NAN behavior was observed
[Dampers]
  [max_T]
    type = MaxIncrement
    increment_type = fractional
    max_increment = 0.5
    variable = T_solid
  []
  [bd_rhoA]
    type = BoundingValueElementDamper
    variable = rhoA
    min_value = 0.0
    min_damping = 0.001
  []
  [bd_rhoEA]
    type = BoundingValueElementDamper
    variable = rhoEA
    min_value = 0.0
    min_damping = 0.001
  []
[]

[Executioner]
  type = Transient

  scheme = bdf2
  start_time = ${t_start}
  end_time = ${t_end}

  dtmin = 1e-4
  dtmax = 5000

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.1
    optimal_iterations = 15
    iteration_window = 2
    growth_factor = 2
    cutback_factor = 0.8
  []

  [Quadrature]
    type = GAUSS
    order = SECOND
  []

  ## PJFNK + superlu seems to make the problem converge faster
  solve_type = PJFNK
  petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -ksp_gmres_restart'
  petsc_options_value = 'lu       superlu_dist                  101'

  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-8
  nl_max_its = 25

  l_tol = 1e-3
  l_max_its = 10

  automatic_scaling = true
  compute_scaling_once = false
[]

[Outputs]
  exodus = false
  [console]
    type = Console
    execute_postprocessors_on = 'NONE'
  []
  [csv]
    type = CSV
    show = 'evaporator_boundary_integral condenser_boundary_integral'
    execute_on = 'INITIAL FINAL'
    enable = false
  []
[]

Load Following Transient Simulation

The load following transient simulation for the Na-HPMR is performed using the same approach as the counterpart simulation for the K-HPMR. The transient is initiated by a significant reduction in the heat removal capability of the secondary coolant loop through alterning the condensor envelope surface boundary condition in the Sockeye model. The heat transfer coefficient (HTC) of that convective boundary condition is reduced by 99.99% (i.e., from 10 W/m-K to 100 W/m-K) to mimic the loss of heat removal capability in the secondary coolant loop. The transient simulation is performed for 2,000 seconds.

Startup Transient Simulation

The startup transient simulation is conducted using a different and simplified approach compared to the load following transient simulation. The reactivity startup of the microreactor is achieved by precise control of the control drum rotation to compensate the reactivity feedback from the increasing temperature to maintain a designated power ramping profile. Such a startup procedure would be expensive to perform if a high-fidelity neutronics modeling approach is used. Alternatively, the startup neutronics can be modeled using point kinetics with a much lower computational cost. However, as the focus of the startup transient simulation is to demonstrate the capability of the Sockeye model in a multiphysics scheme instead of the neutronics behavior during startup, it is not essential to include a low-fidelity neutronics model such as point kinetics in the model. Instead, the startup transient simulation is performed by directly controlling the power ramping profile in the BISON model, while the Griffin model is bypassed. A linear power ramping profile which reaches the nominal power within 3,600 seconds is used for the startup transient simulation. The temperature evolution in different reactor components and the activation of the heat pipes are the focus of the analysis. It is also worth noting that the startup transient involved here is not a frozen startup. Instead, the initial temperature is set at 700 K, which is lower than the activation temperature of the sodium work fluid but higher than the melting point of sodium.

To establish this initial state, the external boundary temperature at the condenser was reduced from 900 K to 700 K. To maintain a steady state condenser temperature of ~900 K at full power despite the 700 K external boundary, the heat transfer coefficient was adjusted to 152 W/m²K. This configuration holds the condenser near 900 K at approximately 1800 W per heat pipe (the design power level). While this simplified boundary condition effectively captures the desired startup temperature behavior, a coupled, high-fidelity secondary system model will be integrated in the near future for more realistic simulations.