SAM Generic PBR Model

Contact: Zhiee Jhia Ooi, [email protected]

Model link: Generic PBR SAM Model

The input file (pbr.i) is a model for the generic 200 MW pebble bed reactor developed using SAM (Hu et al., 2021). The model focuses only on the core of the reactor where the auxiliary components and power conversion system are not included. The design and operation information of the core is obtained primarily from the work by Stewart et al. (Stewart et al., 2021), who obtained geometric data of the core from publicly available figures and models, supplemented with information from the PBMR-400 report published by the Nuclear Energy Agency (OECDNEA, 2013).

The schematic of the core model is shown in Figure 1. The core consists of the pebble bed, reflector, core barrel, reactor pressure vessels (RPV), reactor cavity cooling system (RCCS) panel, and helium layers. Note that the angled lower region of the pebble bed is known as the fuel chute. In the model, helium enters from the bottom of the core and flows upward to the upper inlet plenum via the upcomer. It then flows downward through the pebble bed, carrying heat from the pebbles, to the lower outlet plenum before leaving the core.

Figure 1: Schematic of the core model of the generic pebble bed reactor.

The SAM model is developed with the so-called core channel approach where the pebble bed is modeled with a number of PbCoreChannel components while the reflectors and other solid structures are modeled with the PbCoupledHeatStructure component. Along with them are the PbOneDFluidCompoenent used for modeling the flow channels. As shown in the above figure, the pebble bed is partitioned into five individual blocks radially where each block is a PbCoreChannel component. In this model, the number of PbCoreChannel components in the radial direction is chosen to match the partitioning of the reactivity feedback coefficients in the work by Stewart et al. (Stewart et al., 2021). Furthermore, the fuel chute is modeled with multiple smaller blocks of PbCoreChannel components to account for its complex geometry. Note that the PbCoreChannel in SAM is essentially a one-dimensional hydraulic component with built-in multi-dimensional heat structures. This allows the specification of hydraulic parameters such as flow area, hydraulic diameter, and surface roughness, as well as solid parameters such as the geometry and material properties of the heat structure. Furthermore, heat transfer parameters such as the heat transfer area density can be provided to the component. Solid, fluid, and heat transfer behaviors are then calculated internally by the PbCoreChannel component according to the provided information. In this model, the PbCoreChannel component is set to have spherical heat structures where the heat-transfer geometry is set as pebble bed.

The solid structures outside of the pebble bed are modeled as concentric rings with the PbCoupledHeatStructure components and are thermally coupled to the adjacent PbOneDFluidComponents when necessary. To ensure proper heat transfer between adjacent solid components, they are thermally coupled in both the radial and axial directions using the SurfaceCoupling component with an arbitrarily large heat transfer coefficient. Furthermore, similar approach is applied to thermally couple adjacent PbCoreChannels to model the core wide radial heat conduction in the pebble bed where the outermost PbCoreChannel is further coupled to the inner surface of the reflector. It is important that the core wide radial heat conduction is modeled correctly as it is one of the primary heat removal mechanisms during loss of forced flow transient scenarios. Heat transfer between solid structures across the helium layers is modeled as radiation heat transfer. Note that the helium layers are assumed as stagnant with no natural convection.

Point kinetics are included in the model for the transient load-following simulation. The point kinetics parameters such as the delayed neutron fractions, the neutron lifetimes, the delayed neutron precursor decay constant, the prompt neutron lifetime, and the reactivity feedback coefficients are obtained from the work by Stewart et al. (Stewart et al., 2021). In this model only temperature reactivity feedback from the fuel, moderator, and reflector regions are considered. Other mechanisms such as coolant reactivity feedback, xenon reactivity feedback, and core expansion reactivity feedback are not considered both in this work and the work by Stewart et al. (Stewart et al., 2021). The distribution of the local reactivity coefficients are illustrated in Figure 2 where the reactivity of the pebble bed is shown to be influenced by the feedback from the fuel, moderator, and reflector. To properly account for the reactivity feedback from different regions, it is necessary to divide the pebble heat structures into three layers of fuel, moderator, and reflector with the correct reactivity coefficients prescribed to each layer. Furthermore, this step is necessary because SAM treats fuel (Doppler) reactivity coefficients differently than moderator/reflector reactivity coefficients. Note that no such distinction is made by SAM between the moderator and reflector reactivity coefficients. This means that in the reflector region where the fuel (Doppler) reactivity coefficients are absent, the local moderator and reflector reactivity coefficients from Stewart et al. (Stewart et al., 2021) can be summed and prescribed directly to the heat structures.

Figure 2: Distribution of local reactivity coefficients from Griffin/Pronghorn by Stewart et al. (Stewart et al., 2021).

During steady-state, the one-dimensional fluid components of the model has an inlet velocity 7.17 m/s whose value is determined from the work by Stewart et al. (Stewart et al., 2021) and a pressure outlet boundary condition of 6 MPa. For the solid structures, the top and bottom surfaces are adiabatic while the outer surface of the RCCS panel is given a constant temperature boundary condition of 293.15 K. The total power of the reactor is set at 200 MW. During the transient, the coolant flow rate is reduced to 25% of its nominal value in a span of 900 s, then kept constant at that level for 1800 s, and lastly raised back to the nominal value over 900 s and is held constant to the end of the simulation. Note that other than the coolant flow rate, the other boundary conditions are kept unchanged.

The output files consist of: (1) a csv file that writes all user-specified variables and postprocessors at each time step; (2) a checkpoint folder that saves the snapshots of the simulation data including all meshes, solutions, and stateful object data. They are saved for restarting the run if needed; and (3) an ExodusII file that also has all mesh and solution data. Users can use Paraview to visualize the solution and analyze the data. This tutorial describes the content of the input file, the output files and how the model can be run using the SAM code.

Input Description

SAM uses a block-structured input syntax. Each block begins with square brackets which contain the type of input and ends with empty square brackets. Each block may contain sub-blocks. The blocks in the input file are described in the order as they appear. Before the first block entries, users can define variables and specify their values which are subsequently used in the input model. For example,

emissivity   = 0.8.             # Material emissivity
T_RCCS       = 293.15           # Temperature of RCCS panel
V_in         = 7.186055         # Coolant inlet velocity

In SAM, comments are entered after the # sign

Global parameters

This block contains the parameters such as global initial pressure, velocity, and temperature conditions, the scaling factors for primary variable residuals, etc. For example, to specify global pressure of 6.0e6 Pa, the user can input



global_init_P	 = 	6.0e6

This block also contains a sub-block PBModelParams which specifies the modeling parameters associated with the primitive-variable based fluid model. New users should leave this sub-block unchanged.

EOS

This block specifies the Equation of State. The user can choose from built-in fluid library for common fluids like air, nitrogen, helium, sodium, molten salts, etc. The user can also input the properties of the fluid as constants or function of temperature. For example, the built-in eos for helium can be input as

[EOS]
  [eos]
    type = HeEquationOfState
  []
[]
(htgr/generic-pbr/pbr.i)

Functions

Users can define functions for parameters used in the model. These include temporal, spatial, and temperature dependent functions. For example, users can input enthalpy as a function of temperature, power history as a function of time, or power profile as a function of position. Users should be cautious when time-dependent and temperature-dependent functions are in use to avoid confusion between time (t) and temperature (T). Currently, only certain parameters can be made temperature-dependent. These include material properties and the h_gap for the SurfaceCoupling component. The input below specifies the thermal conductivity of fuel pebbles as a function of temperature (in K).

[Functions]
  [k_TRISO]
    # Effective fuel thermal conductivity calculated with Chiew & Glandt model: kp = 4.13 W/(m K)
    type = PiecewiseLinear
    x = '600                700             800             900             1000            1100            1200            1300            1400            1500'
    y = '55.6153061 51.02219975 47.11901811 43.95203134 41.16924224 38.85202882 36.89323509 35.04777834 33.20027175 31.3520767'
  []
[]
(htgr/generic-pbr/pbr.i)

MaterialProperties

Material properties are input in this block. The values can be constants or temperature dependent as defined in the Functions block. For example, the properties of fuel pebbles are input as

[MaterialProperties]
  [fuel-mat]
    type = SolidMaterialProps
    k = k_TRISO
    Cp = 1697 # Specific heat
    rho = 1780 # Density assumed the same as normal graphite as suggested by NEA report
  []
[]
(htgr/generic-pbr/pbr.i)

The thermal conductivity is defined by the function k_TRISO which appears under the Functions block.

ComponentInputParameters

This block is used to input common features for Components (section below) so that these common features do not need to be repeated in the inputs for Components later on. For example, if pipes are used in various parts of the model and the pipes all have the same diameter, then the diameter can be specified in ComponentsInputParameters and it applies to all pipes used in the model.

[ComponentInputParameters]
  [CoolantChannel]
    type = PBOneDFluidComponentParameters
    eos = eos
    HTC_geometry_type = Pipe
  []
  [Graphite-Fuel]
    type = HeatStructureParameters
    hs_type = cylinder
    dim_hs = 2
    material_hs = 'fuel-mat'
  []
  [Graphite-Reflector]
    type = HeatStructureParameters
    hs_type = cylinder
    dim_hs = 2
    material_hs = 'graphite-mat'
  []
  [SS-Barrel]
    type = HeatStructureParameters
    hs_type = cylinder
    dim_hs = 2
    material_hs = 'ss-mat'
  []
  [SS-RPV]
    type = HeatStructureParameters
    hs_type = cylinder
    dim_hs = 2
    material_hs = 'ss-rpv-mat'
  []
  [Graphite-Reflector-Porous]
    type = HeatStructureParameters
    hs_type = cylinder
    dim_hs = 2
    material_hs = 'graphite-porous-mat'
  []
[]
(htgr/generic-pbr/pbr.i)

Components

This block provides the specifications for all components that make up the core. The components consist of: reactor, core channels, coolant channels, heat structures (fuel, reflectors, core barrel, RPV and RCCS), junctions for connecting components, and point kinetics. In the reactor component, the reactor power is an input and this includes normal operating power and decay heat. Point kinetics can also be enabled for the core channels through the reactor component.

[Components]
  [reactor]
    type = ReactorPower
    initial_power = 200e6 # Initial total reactor power
    operating_power = 200e6
    pke = 'pke'
    enable_decay_heat = true
    isotope_fission_fraction = '1.0 0.0 0.0 0.0'
  []
[]
(htgr/generic-pbr/pbr.i)

Additional point kinetics information is needed by the model and can be provided through the PKE component. Such information includes the delayed neutron fractions, the prompt neutron lifetimes, the delayed neutron precursor decay constant, the names of components with reactivity feedback, and the time for reactivity feedback to start.

[Components]
  [pke]
    type = PointKinetics
    lambda = '1.334e-2 3.273e-2 1.208e-1 3.029e-1 8.501e-1 2.855' # From Griffin/Pronghorn (Table 9 of INL Report [1])
    LAMBDA = 6.519e-4 # From Griffin (Table 9 of INL Report [1])
    betai = '2.344e-4 1.210e-3 1.150e-3 2.588e-3 1.070e-3 4.486e-4' # From Griffin/Pronghorn (Table 9 of INL Report [1])
    feedback_start_time = -200
    feedback_components = 'F-1 F-2 F-3 F-4 F-5 F-6 F-7 F-8 F-9 F-10 F-11 F-12 F-13 F-14 F-15
                               BR-1 BR-2 BR-3 BR-4 BR-5 BR-6 BR-7 BR-8 BR-9 BR-10 BR-11 BR-12 BR-13 BR-14 BR-15 BR-16 BR-17 BR-18 BR-19 BR-20
                               R-4 R-5 R-6 R-7 R-8 R-9 R-10 R-11 R-12 R-13 R-14 R-15 R-16 R-17 R-18 R-19 R-20 R-21 R-22 R-23'
    irk_solver = true # Turn on IRK solver to speed up the simulation for PKE calculation
  []
[]
(htgr/generic-pbr/pbr.i)

The coolant channels are modeled as 1-D fluid flow components, and heat structures are modeled as 2-D components. In the bottom reflector region, each coolant channel communicates with its adjacent heat structures through the name_comp_left and name_comp_right variables such as shown below

[Components]
  [BR-6]
    type = PBCoupledHeatStructure
    input_parameters = Graphite-Reflector-Porous
    radius_i = 0.99
    HS_BC_type = 'Coupled Coupled'
    position = '0.0 0 2.065'
    orientation = '0 0 1'
    length = 0.135
    elem_number_axial = 3
    HT_surface_area_density_left = ${aw_BR6_left}
    name_comp_left = CH-5
    HT_surface_area_density_right = ${aw_BR6_right}
    name_comp_right = CH-6
    width_of_hs = '0.1610973'
    elem_number_radial = 5
    moderator_reactivity_feedback = true
    n_layers_moderator = 1
    moderator_reactivity_coefficients = '${moderator_coef_BR6}'
  []
[]
(htgr/generic-pbr/pbr.i)

Adjacent heat structures are thermally coupled using SurfaceCoupling to ensure temperature continuity in the radial direction

[Components]
  [coupling_radial_R9_R10]
    type = SurfaceCoupling
    coupling_type = GapHeatTransfer
    h_gap = ${h_gap}
    surface1_name = R-9:outer_wall
    surface2_name = R-10:inner_wall
    radius_1 = 1.72
  []
[]
(htgr/generic-pbr/pbr.i)

and axial direction

[Components]
  [coupling_axial_R5_R8]
    type = SurfaceCoupling
    coupling_type = GapHeatTransfer
    h_gap = ${h_gap}
    surface1_name = R-8:top_wall
    surface2_name = R-5:bottom_wall
    radius_1 = 1.9
  []
[]
(htgr/generic-pbr/pbr.i)

Note that the coupling in the radial direction is performed on the inner_wall and outer_wall of a heat structure while the coupling in the axial direction is performed on the top_wall and bottom_wall. In both cases the heat transfer coefficient, h_gap, is set to an arbitrarily large value to ensure temperature continuity.

Similarly, SurfaceCoupling is used to thermally couple adjacent core channels to allow for radial heat transfer between the core channels

[Components]
  [coupling_radial_F1_F2]
    type = SurfaceCoupling
    coupling_type = PebbleBedHeatTransfer
    surface1_name = F-1:solid:outer_wall
    surface2_name = F-2:solid:outer_wall
    h_gap = ${h_ZBS_1}
    radius_1 = ${F1_radius}
    radius_2 = ${F2_radius}
    area_ratio = ${F1_F2_ratio}
  []
[]
(htgr/generic-pbr/pbr.i)

Note the difference in coupling_type used for the heat structures and the core channels. The former uses GapHeatTransfer and the latter uses PebbleBedHeatTransfer.

In SAM, 1-D fluid components are connected using PBSingleJunction. The following input is used to connect the outlet of component CH-1 to the inlet of component CH-3

[Components]
  [junc_CH1_CH3]
    type = PBSingleJunction
    eos = eos
    inputs = 'CH-1(out)'
    outputs = 'CH-3(in)'
  []
[]
(htgr/generic-pbr/pbr.i)

Note that core channels can also be connected with the same approach.

For the components with point kinetics, additional information is required. As mentioned previously, the pebbles in the core channels are divided into three layers of fuel, moderator, and reflector. Using the F-6 core channel block as an example,

    n_heatstruct = 3
    elem_number_of_hs = '5 5 5'
    power_fraction = '${power_fraction_6} 0 0'
    pke_material_type = 'FuelDoppler Moderator Moderator'
    fuel_doppler_coef = ${fuel_coef_F6}
    n_layers_doppler = 1
    n_layers_moderator = 1
    moderator_reactivity_coefficients = '${moderator_coef_F6}; ${reflector_coef_F6}'
(htgr/generic-pbr/pbr.i)

the number of layers can be set using the n_heatstruct option. The names, materials, widths, numbers of elements in the radial direction, and power fractions can be provided as a string using the name_of_hs, material_hs, width_of_hs, elem_number_of_hs, and power_fraction options, respectively. Note that because heat is generated only in the fuel region, the power fraction of the this region is set to the desired value with those of the moderator and reflector regions are set to zero. The types of reactivity feedback of each region is specified with the pke_material_type option. The fuel (Doppler) reactivity coefficients are provided with the fuel_doppler_coef option while the moderator/reflector reactivity coefficients are provided with the moderator_reactivity_coefficients option. The numbers of axial elements for point kinetics calculation in the fuel (Doppler) and moderator/reflector regions are provided with the n_layers_doppler and n_layers_moderator options, respectively. Note that the number of PKE axial elements must match the length of the provided reactivity coefficients.

By default, the reactivity feedback in the heat structures of a core channel is calculated using the layered-average solid temperature. However, users can choose to calculate the reactivity feedback in the fuel region with the fuel kernel temperature instead. This can be done by setting

    fuel_kernel_temperature = true
    use_kernel_temperature_doppler = true
(htgr/generic-pbr/pbr.i)

Currently, the fuel kernel temperature is calculated using the analytical model implemented in the TINTE code (Gerwin et al., 2010). Additional information is needed as,

    molar_mass_uranium = 0.238
    molar_mass_kernel_material = 0.27
    heavy_metal_loading = 0.009
    local_power_fraction = 1.0
    heat_capacity_kernel_material = 300.0
    modified_heat_flux_resistance = 2.6e-6
  []
[]
(htgr/generic-pbr/pbr.i)

Detailed descriptions of these parameters are available in the SAM Theory Manual (Hu et al., 2021) and the TINTE report (Gerwin et al., 2010).

Postprocessors

This block is used to specify the output quantities written to a csv file that can be further processed in Excel or other processing tools such as Python and Matlab. For example, to output the maximum fuel temperature in F-1:

[Postprocessors]
  [max_Tsolid_F1]
    type = NodalExtremeValue
    block = 'F-1:solid:fuel'
    variable = T_solid
  []
[]
(htgr/generic-pbr/pbr.i)

and the mass flow rate in F-1

[Postprocessors]
  [TotalMassFlowRateInlet_1]
    type = ComponentBoundaryFlow
    input = F-1(in)
  []
[]
(htgr/generic-pbr/pbr.i)

Preconditioning

This block describes the preconditioner to be used by the preconditioned JFNK solver (available through PETSc). Two options are currently available, the single matrix preconditioner (SMP) and the finite difference preconditioner (FDP). The theory behind the preconditioner can be found in the SAM Theory Manual (Hu et al., 2021). New users can leave this block unchanged.

[Preconditioning]
  [SMP_PJFNK]
    type = SMP
    full = true
    solve_type = 'PJFNK'
    petsc_options_iname = '-pc_type -ksp_gmres_restart'
    petsc_options_value = 'lu 101'
  []
[]
(htgr/generic-pbr/pbr.i)

Executioner

This block describes the calculation process flow. The user can specify the start time, end time, time step size for the simulation. The user can also choose to use an adaptive time step size with the IterationAdaptiveDT time stepper. Other inputs in this block include PETSc solver options, convergence tolerance, quadrature for elements, etc. which can be left unchanged.

[Executioner]
  type = Transient
  dt = 1
  dtmin = 0.001
  start_time = -172800
  end_time = -210
  dtmax = 10000
  nl_rel_tol = 1e-8
  l_tol = 1e-8
  nl_abs_tol = 1e-8
  nl_max_its = 30
  l_max_its = 200

  # The TimeStepper block can be activated if users
  # wish to enable adaptive time-stepper.
  # [./TimeStepper]
  #     type = IterationAdaptiveDT
  #     growth_factor = 1.25
  #     optimal_iterations = 8
  #     linear_iteration_ratio = 150
  #     dt = 1
  #     cutback_factor = 0.8
  #     cutback_factor_at_failure = 0.8
  # [../]
  [Quadrature]
    type = SIMPSON
    order = SECOND
  []
[]
(htgr/generic-pbr/pbr.i)

Running the Simulation

SAM can be run on Linux, Unix, and MacOS. Due to its dependence on MOOSE, SAM is not compatible with Windows. SAM can be run from the shell prompt as shown below



sam-opt -i pbr.i

Acknowledgements

This work is supported by U.S. DOE Office of Nuclear Energy’s Nuclear Energy Advanced Modeling and Simulation (NEAMS) program. The submitted manuscript has been created by UChicago Argonne LLC, Operator of Argonne National Laboratory (“Argonne”). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The authors would like to acknowledge the support and assistance from Dr. Ryan Stewart and Dr. Paolo Balestra of Idaho National Laboratory in the completion of this work.

References

  1. H. Gerwin, W. Scherer, A. Lauer, and I. Clifford. TINTE – Nuclear Calculation Theory Description Report. Technical Report Jul-4317, Safety Research and Reactor Technology, Institute for Energy Research, 2010.[BibTeX]
  2. R. Hu, L. Zou, G. Hu, D. Nunez, T. Mui, and T. Fei. SAM Theory Manual. Technical Report ANL/NE-17/4 Rev. 1, Argonne National Laboratory, Lemont, IL, 2021.[BibTeX]
  3. OECDNEA. PBMR Coupled Neutronics/ Thermal-hydraulics Transient Benchmark the PBMR-400 Core Design. Technical Report NEA/NSC/DOC(2013)10, Organisation for Economic Co-operation and Development, Nuclear Energy Agency, Paris, France, 2013.[BibTeX]
  4. R. Stewart, D. Reger, and P. Balestra. Demonstrate Capability of NEAMS Tools to Generate Reactor Kinetics Parameters for Pebble-Bed HTGRs Transient Modeling. Technical Report INL/EXT-21-64176, Idaho National Laboratory, Idaho Falls, ID, 2021.[BibTeX]