LEU Fuel Pulse

Contact: Adam Zabriskie, [email protected]

Model link: LEU Fuel Pulse Model

Model Description

This model represents a greatly simplified version of the Transient Reactor Test Facility (TREAT). The model explores the effect on pulse feedback mechanisms from microscale heterogeneity represented by hypothetical 20 m LEU fuel grains. This model couples the heat equation with the neutron governing equation for temperature feedback during a short pulse. The pulse is produced by inserting a reactivity of 4.56 % k/k. The model's simplified geometry consists of a homogeneous cube of TREAT-like fuel surrounded by a cube shell reflector. The simplified geometry has no cladding, coolant channels, or other structures (Zabriskie, 2019). To reduce computation costs, symmetric boundaries are applied to a 1/8th quadrant of the simplified reactor, which is shown in Figure 1. Outer boundaries on the reflector have vacuum and adiabatic boundary conditions for neutron and heat transport.

Figure 1: TREAT model octant

The dimensions of the macro-scale simulation model are given in Table 1, and the fuel grain characteristics are given in Table 2.

Table 1: Reflected cube reactor dimensions.

CharacteristicValue (cm)
Reactor side length
Reflector thickness
Total side length
Core mesh element cube length
Reflector mesh element cube length

Table 2: 20 m radius LEU fuel grain characteristics.

CharacteristicValue (m)
Fuel grain radius
Outer boundary radius
Moderator transition thickness
Moderator shell total thickness
Moderator shell boundary thickness
Fuel grain mesh element
Moderator transition mesh element
Moderator regular mesh element

The neutronics analysis is performed on a mesh characterized by centimeter-sized elements, while grain treatment is on the micrometer element scale where heat conduction problems are solved on microscopic domains encompassing a single fuel grain and a proportional amount of graphite. This model makes use of the Griffin code for neutron transport and the MOOSE heat transfer module for heat transport and radiation. Serpent 2 was used to prepare cross section libraries (Zabriskie, 2019), (Zabriskie et al., 2019).

The particles are spaced cm apart in a regular grid lattice, and are spaced half that distance when next to a boundary. There are particles in each octant of the full reactor.

The power density, calculated from the sampled flux at each location, is transferred to a single microscale simulation.

The MOOSE MultiApp system handles the multiscale and information transfer needs of the simulation. The macroscale simulation provides the average temperature at each microscale simulation location as an outer boundary condition. With the sampled flux at each location and this outer boundary condition, the thermal solution of each microscale grain provides an average moderator shell temperature and fuel grain temperature to the locations in the macroscale simulation. These average temperatures then adjust the cross section, providing temperature feedback.

Input Files

There are four input files for this model.

Transient Pulse Initial Conditions

This input file sets the initial conditions from which the transient pulse starts.

commentnote

Some code blocks appearing in this file also appear in other files. For brevity, these common blocks are only shown in this first input file, but they will be mentioned in each input file that they appear.

Mesh

This block creates the mesh used for the simplified reactor.

[Mesh]
  # Simple Reflected Cube Reactor
  # Start at zero, half core length reflector thickness; 1/8th symmetric
  # Element size: reflector 5.08 cm and core 4.827616834 cm
  # 14 in half-core, 12 in reflector
  [init_mesh]
    type = CartesianMeshGenerator
    dim = 3
    dx = '67.58663568 60.96'
    dy = '67.58663568 60.96'
    dz = '67.58663568 60.96'
    ix = '14 12'
    iy = '14 12'
    iz = '14 12'
  []
  [set_core_id]
    type = SubdomainBoundingBoxGenerator
    block_id = 10
    input = init_mesh
    bottom_left = '0.0 0.0 0.0'
    top_right = '67.58663568 67.58663568 67.58663568'
  []
[]
(htgr/leu_pulse/init_refcube.i)

TransportSystems

This block sets up the neutron transport governing equation in the diffusion approximation form to be used in the simulation.

[TransportSystems]
  # In 3D, back = 0, bottom = 1, right = 2, top = 3, left = 4, front = 5
  # back is -z, bottom is -y, right is +x
  #Boundary Conditions#
  G = 6
  ReflectingBoundary = '0 1 4'
  VacuumBoundary = '2 3 5'
  equation_type = eigenvalue
  for_adjoint = false
  particle = neutron
  [diffing]
    family = LAGRANGE
    n_delay_groups = 6
    order = FIRST
    scheme = CFEM-Diffusion
  []
[]
(htgr/leu_pulse/init_refcube.i)

AuxVariables

This block sets up the auxilliary variables that will be calculated and tracked throughout the runtime of the simulation.

[AuxVariables]
  [temperature]
    family = LAGRANGE
    initial_condition = 300.0
    order = FIRST
  []
  [Boron_Conc]
    family = MONOMIAL
    initial_condition = 1.89489259748
    order = CONSTANT
  []
  [PowerDensity]
    block = '10'
    family = MONOMIAL
    order = CONSTANT
  []
  [avg_coretemp]
    block = 0
    family = LAGRANGE
    initial_condition = 300.0
    order = FIRST
  []
[]
(htgr/leu_pulse/init_refcube.i)

AuxKernels

This block defines the equations used to calculate auxilliary values for use in the simulation.

[AuxKernels]
  [PowerDensityCalc]
    type = VectorReactionRate
    block = '10'
    cross_section = kappa_sigma_fission
    dummies = UnscaledTotalPower
    execute_on = 'initial linear'
    scalar_flux = 'sflux_g0 sflux_g1 sflux_g2 sflux_g3 sflux_g4 sflux_g5'
    scale_factor = PowerScaling
    variable = PowerDensity
  []
  [Set_coreT]
    type = SetAuxByPostprocessor
    block = 0
    execute_on = 'linear timestep_end'
    postproc_value = avg_coretemp
    variable = avg_coretemp
  []
[]
(htgr/leu_pulse/init_refcube.i)

Postprocessors

This block processes simulation solution data for output.

[Postprocessors]
  [UnscaledTotalPower]
    type = FluxRxnIntegral
    block = 10
    coupled_flux_groups = 'sflux_g0 sflux_g1 sflux_g2 sflux_g3 sflux_g4 sflux_g5'
    cross_section = kappa_sigma_fission
    execute_on = linear
  []
  [PowerScaling]
    # 100 kW for entire TREAT
    type = PowerModulateFactor
    execute_on = linear
    power_pp = UnscaledTotalPower
    rated_power = 12500.0
  []
  [avg_coretemp]
    type = ElementAverageValue
    block = 10
    execute_on = linear
    outputs = all
    variable = temperature
  []
  [avg_refltemp]
    type = ElementAverageValue
    block = 0
    execute_on = linear
    outputs = all
    variable = temperature
  []
  [avg_powerden]
    type = ElementAverageValue
    block = 10
    execute_on = timestep_end
    outputs = all
    variable = PowerDensity
  []
  [ScaledTotalPower]
    type = ElementIntegralVariablePostprocessor
    block = 10
    execute_on = linear
    variable = PowerDensity
  []
  [delta_time]
    type = TimestepSize
  []
  [nl_steps]
    type = NumNonlinearIterations
  []
  [lin_steps]
    type = NumLinearIterations
  []
  [Eq_TREAT_Power]
    type = ScalePostprocessor
    scaling_factor = 2469860.77609
    value = avg_powerden
  []
[]
(htgr/leu_pulse/init_refcube.i)

Materials

This block sets up the material properties for use in the governing equations. The files leu_20r_is_6g_d.xml and leu_macro_6g.xml contain the cross sections.

[Materials]
  [neut_mix]
    type = CoupledFeedbackNeutronicsMaterial
    block = 10
    densities = '0.998448391539 0.00155160846058'
    grid_names = 'Tfuel Tmod Rod'
    grid_variables = 'temperature temperature Boron_Conc'
    isotopes = 'pseudo1 pseudo2'
    library_file = 'cross_sections/leu_20r_is_6g_d.xml'
    library_name = leu_20r_is_6g_d
    material_id = 1
    plus = true
  []
  [kth]
    # Volume weighted harmonic mean
    # Divided fg_kth by 100 to get it into cm
    type = ParsedMaterial
    coupled_variables = 'temperature'
    block = 10
    constant_expressions = '3.35103216383e-08 1.31125888571e-07 2.14325144175e-05 0.3014 0.01046 1.0 0.05 1.5 1.0'
    constant_names = 'vol_fg vol_fl vol_gr gr_kth fl_kth beta p_vol sigma kap3x'
    property_name = 'thermal_conductivity'
    expression = 'lt := temperature / 1000.0; fresh := (100.0 / (6.548 + 23.533 * lt) + 6400.0 * exp(-16.35 / lt) / pow(lt, 5.0/2.0)) / 100.0; kap1d := (1.09 / pow(beta, 3.265) + 0.0643 * sqrt(temperature) / sqrt(beta)) * atan(1.0 / (1.09 / pow(beta, 3.265) + sqrt(temperature) * 0.0643 / sqrt(beta))); kap1p := 1.0 + 0.019 * beta / ((3.0 - 0.019 * beta) * (1.0 + exp(-(temperature - 1200.0) / 100.0))); kap2p := (1.0 - p_vol) / (1.0 + (sigma - 1.0) * p_vol); kap4r := 1.0 - 0.2 / (1.0 + exp((temperature - 900.0) / 80.0)); fg_kth := fresh * kap1d * kap1p * kap2p * kap3x * kap4r; (vol_fg + vol_fl + vol_gr) / (vol_fg / fg_kth + vol_fl / fl_kth + vol_gr / gr_kth)'
  []
  [rho_cp]
    # Volume weighted arithmetic mean (Irradiation has no effect)
    type = ParsedMaterial
    coupled_variables = 'temperature'
    block = 10
    constant_expressions = '3.35103216383e-08 2.1563640306e-05 0.0018 0.010963'
    constant_names = 'vol_fg vol_gr rho_gr rho_fg'
    property_name = 'heat_capacity'
    expression = 'lt := temperature / 1000.0; gr_rhocp := rho_gr / (11.07 * pow(temperature, -1.644) + 0.0003688 * pow(temperature, 0.02191)); fink_cp := 52.1743 + 87.951 * lt - 84.2411 * pow(lt, 2) + 31.542 * pow(lt, 3) - 2.6334 * pow(lt, 4) - 0.71391 * pow(lt, -2); fg_rhocp := rho_fg * fink_cp / 267.2 * 1000.0; (vol_fg * fg_rhocp + vol_gr * gr_rhocp) / (vol_fg + vol_gr)'
  []
  [neut_refl]
    type = CoupledFeedbackNeutronicsMaterial
    block = 0
    densities = '1'
    grid_names = 'Trefl Tcore Rod'
    grid_variables = 'temperature avg_coretemp Boron_Conc'
    isotopes = 'pseudo'
    library_file = 'cross_sections/leu_macro_6g.xml'
    library_name = leu_macro_6g
    material_id = 2
    plus = true
  []
  [ref_kth]
    type = GenericConstantMaterial
    block = 0
    prop_names = 'thermal_conductivity'
    prop_values = '0.3014'
  []
  [ref_rho_cp]
    type = ParsedMaterial
    coupled_variables = 'temperature'
    block = '0'
    constant_expressions = '0.0018'
    constant_names = 'rho_gr'
    property_name = 'heat_capacity'
    expression = 'rho_gr / (11.07 * pow(temperature, -1.644) + 0.0003688 * pow(temperature, 0.02191))'
  []
[]
(htgr/leu_pulse/init_refcube.i)

Preconditioner

This block describes the preconditioner used by the solver.

[Preconditioning]
  [SMP_full]
    type = SMP
    full = true
    petsc_options = '-snes_ksp_ew -snes_converged_reason'
    petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart -pc_hypre_boomeramg_max_iter -pc_hypre_boomeramg_tol'
    petsc_options_value = 'hypre boomeramg 101 20 1.0e-6'
    solve_type = 'PJFNK'
  []
[]
(htgr/leu_pulse/init_refcube.i)

Executioner

This block describes the simulation solution approach.

[Executioner]
  type = Eigenvalue
  free_power_iterations = 4
  l_max_its = 100
  l_tol = 1e-4
  nl_abs_tol = 1e-8
  nl_max_its = 200
  nl_rel_tol = 1e-7
[]
(htgr/leu_pulse/init_refcube.i)

Adjoint initial conditions

This input file provides the values needed to calculate the point kinetics equation parameters.

This input file is also almost identical to the previous input file, with a few alterations:

  • In the Mesh block, a different name is used for the mesh.

  • In the TransportSystems block, the input for_adjoint is set to true.

  • In the AuxVariables block, [Boron_Conc] sub-block, initial_condition is set to 2.0.

Microscale Particles

This input file describes a microscale particle heat solution. This input is replicated at many locations in the macroscale simulation.

Mesh

This block sets the meshes for the fission damage layer and fuel grain.

[Mesh]
  coord_type = RSPHERICAL
  [leu_mesh]
    type = CartesianMeshGenerator # Remember in cm
    dim = 1
    dx = '0.0020 0.0014 0.0020 0.011875711758' # cm
    # radius particle, damage layer, transition graphite, b_r - dam_lay - 2a_r
    ix = '10 7 10 40' # number of elements
  []
  [set_damlay_id]
    type = SubdomainBoundingBoxGenerator
    input = leu_mesh
    block_id = 11 # Damage Layer id
    bottom_left = '-0.0001 -0.0001 -0.0001'
    top_right = '0.0034 0.0001 0.0001' # x should be a_r
  []
  [set_grain_id]
    type = SubdomainBoundingBoxGenerator
    input = set_damlay_id
    block_id = 10 # Fuel grain id is now 10
    bottom_left = '-0.0001 -0.0001 -0.0001'
    top_right = '0.0020 0.0001 0.0001' # x should be a_r
  []
[]
(htgr/leu_pulse/ht_20r_leu_fl.i)

Variables

This block sets the temperature variable of the governing equation.

[Variables]
  [temperature]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300.0 # K
  []
[]
(htgr/leu_pulse/ht_20r_leu_fl.i)

Kernels

This block defines the heat conduction governing equation.

[Kernels]
  [HeatConduction]
    type = HeatConduction
    variable = temperature
  []
  [HeatStorage]
    type = HeatCapacityConductionTimeDerivative
    variable = temperature
  []
  [HeatSource_fg]
    type = CoupledForce
    v = PowerDensity
    variable = temperature
    block = 10
    coef = 508.59714978 # Scales to energy in grain
  []
  [HeatSource_dl]
    type = CoupledForce
    v = PowerDensity
    variable = temperature
    block = 11
    coef = 34.7291950039 # Scales to energy in damage layer
  []
[]
(htgr/leu_pulse/ht_20r_leu_fl.i)

AuxVariables

This block sets up the power density and integrated power auxiliary variables.

[AuxVariables]
  [PowerDensity]
    order = CONSTANT
    family = MONOMIAL
  []
  [IntegralPower]
    order = CONSTANT
    family = MONOMIAL
    [InitialCondition]
      type = ConstantIC
      value = 0.0
    []
  []
[]
(htgr/leu_pulse/ht_20r_leu_fl.i)

AuxKernels

This block calculates needed auxiliary values during the simulation.

[AuxKernels]
  [PowerIntegrator]
    type = VariableTimeIntegrationAux
    variable = IntegralPower #J/cm^3
    variable_to_integrate = PowerDensity
    execute_on = timestep_end
    order = 2
  []
  [SetPowerDensity]
    type = SetAuxByPostprocessor
    execute_on = 'timestep_begin timestep_end'
    postproc_value = local_power_density # Not scaled
    variable = PowerDensity
  []
[]
(htgr/leu_pulse/ht_20r_leu_fl.i)

Postprocessors

This block calculates outputs from solution data.

[Postprocessors]
  [local_power_density]
    type = Receiver
  []
  [TotalPower]
    type = ElementIntegralVariablePostprocessor
    variable = PowerDensity
  []
  [IntegratedPower]
    type = ElementIntegralVariablePostprocessor
    variable = IntegralPower
  []
  [avg_graphtemp]
    type = ElementAverageValue
    block = '0 11'
    variable = temperature
  []
  [avg_graintemp]
    type = ElementAverageValue
    block = 10
    variable = temperature
  []
  [avg_gr_rhocp]
    type = ElementAverageMaterialProperty
    block = '0 11'
    mat_prop = 'heat_capacity'
  []
  [avg_gr_kth]
    type = ElementAverageMaterialProperty
    block = '0 11'
    mat_prop = 'thermal_conductivity'
  []
  [avg_fg_rhocp]
    type = ElementAverageMaterialProperty
    block = 10
    mat_prop = 'heat_capacity'
  []
  [avg_fg_kth]
    type = ElementAverageMaterialProperty
    block = 10
    mat_prop = 'thermal_conductivity'
  []
  [delta_time]
    type = TimestepSize
  []
  [nl_steps]
    type = NumNonlinearIterations
  []
  [lin_steps]
    type = NumLinearIterations
  []
[]
(htgr/leu_pulse/ht_20r_leu_fl.i)

Materials

This block contains material properties used in the governing equation.

[Materials]
  # Pure Graphite
  [graph_kth]
    type = GenericConstantMaterial
    block = 0
    prop_names = 'thermal_conductivity'
    prop_values = '0.3014' # W/cm K
  []
  [graph_rho_cp]
    type = ParsedMaterial
    block = '0 11' # Irradiation does not affect rho * cp
    constant_names = 'rho_gr'
    constant_expressions = '0.0018' # kg/cm3
    coupled_variables = 'temperature' # Variable
    property_name = 'heat_capacity' # Property name
    expression = 'rho_gr / (11.07 * pow(temperature, -1.644) + 0.0003688 * pow(temperature, 0.02191))' # J/cm3 K
  []
  # Fission Damaged Graphite
  [damlay_kth]
    type = GenericConstantMaterial
    block = 11 # Damage layer of graphite
    prop_names = 'thermal_conductivity'
    prop_values = '0.01046' # W/cm K
  []

  # Fuel Grain UO2
  [grain_kth]
    type = ParsedMaterial
    block = 10
    constant_names = 'beta p_vol sigma kap3x'
    constant_expressions = '1.0 0.05 1.5 1.0'
    coupled_variables = 'temperature' # Variable
    property_name = 'thermal_conductivity' # Property name
    expression = 'lt := temperature / 1000.0; fresh := (100.0 / (6.548 + 23.533 * lt) + 6400.0 * exp(-16.35 / lt) / pow(lt, 5.0/2.0)) / 100.0; kap1d := (1.09 / pow(beta, 3.265) + 0.0643 * sqrt(temperature) / sqrt(beta)) * atan(1.0 / (1.09 / pow(beta, 3.265) + sqrt(temperature) * 0.0643 / sqrt(beta))); kap1p := 1.0 + 0.019 * beta / ((3.0 - 0.019 * beta) * (1.0 + exp(-(temperature - 1200.0) / 100.0))); kap2p := (1.0 - p_vol) / (1.0 + (sigma - 1.0) * p_vol); kap4r := 1.0 - 0.2 / (1.0 + exp((temperature - 900.0) / 80.0)); fresh * kap1d * kap1p * kap2p * kap3x * kap4r'
    # Divided fresh by 100 to get it into cm
  []
  [grain_rho_cp]
    type = ParsedMaterial
    block = 10
    constant_names = 'rho_fg'
    constant_expressions = '0.010963' # kg/cm3
    coupled_variables = 'temperature' # Variable
    property_name = 'heat_capacity' # Property name
    expression = 'lt := temperature / 1000.0; fink_cp := 52.1743 + 87.951 * lt - 84.2411 * pow(lt, 2) + 31.542 * pow(lt, 3) - 2.6334 * pow(lt, 4) - 0.71391 * pow(lt, -2); rho_fg * fink_cp / 267.2 * 1000.0' # J/cm3 K
  []
[]
(htgr/leu_pulse/ht_20r_leu_fl.i)

Preconditioning

This block is identical to the [Preconditioning] block in the Transient Pulse Initial Conditions file.

Executioner

This block describes the simulation solution approach.

[Executioner]
  type = Transient

  # Solver parameters
  l_tol = 1e-3
  l_max_its = 100
  nl_max_its = 200
  nl_abs_tol = 1e-8
  nl_rel_tol = 1e-8

  # Transient parameters
  start_time = 0.0
  end_time = 10.0
  [TimeStepper]
    type = ConstantDT
    dt = 0.005
    growth_factor = 1.5
  []
[]
(htgr/leu_pulse/ht_20r_leu_fl.i)

Main Input File

This is the controlling input file that must be run. All other input files are controlled by the MultiApp system in the main input file.

Mesh

This block is almost identical to the Mesh block in the Transient Pulse Initial Conditions file, with a different name used for the mesh.

MultiApps

This block defines the MOOSE MultiApp system. The file refcube_sub_micro.txt is a list of coordinates for the positions of each of the 125 fuel grains.

[MultiApps]
  [initial_solve]
    type = FullSolveMultiApp
    execute_on = initial
    input_files = 'init_refcube.i'
    positions = '0 0 0'
  []
  [init_adj]
    type = FullSolveMultiApp
    execute_on = initial
    input_files = 'adj_refcube.i'
    positions = '0 0 0'
  []
  [micro]
    type = TransientMultiApp
    execute_on = timestep_end
    input_files = 'ht_20r_leu_fl.i'
    max_procs_per_app = 1
    positions_file = 'refcube_sub_micro.txt'
  []
[]
(htgr/leu_pulse/refcube.i)

Transfers

This block defines the passing of solution data between MultiApp simulations.

[Transfers]
  #Below are communication for adjoint IQS for PKE
  #Below are communication for Multiscale with micro-subs
  [copy_flux]
    type = TransportSystemVariableTransfer
    execute_on = initial
    from_transport_system = diffing
    from_multi_app = initial_solve
    to_transport_system = diffing
  []
  [copy_pp]
    # Scales everything to the initial power.
    type = MultiAppPostprocessorTransfer
    execute_on = initial
    from_postprocessor = UnscaledTotalPower
    from_multi_app = initial_solve
    reduction_type = maximum
    to_postprocessor = UnscaledTotalPower
  []
  [copy_sf]
    # Scales factor
    type = MultiAppPostprocessorTransfer
    execute_on = initial
    from_postprocessor = PowerScaling
    from_multi_app = initial_solve
    reduction_type = maximum
    to_postprocessor = PowerScaling
  []
  [copy_adjoint_vars]
    type = MultiAppCopyTransfer
    execute_on = initial
    from_multi_app = init_adj
    source_variable = 'sflux_g0 sflux_g1 sflux_g2 sflux_g3 sflux_g4 sflux_g5'
    variable = 'adjoint_flux_g0 adjoint_flux_g1 adjoint_flux_g2 adjoint_flux_g3 adjoint_flux_g4 adjoint_flux_g5'
  []
  [powden_down]
    type = MultiAppVariableValueSamplePostprocessorTransfer
    to_multi_app = micro
    postprocessor = local_power_density
    source_variable = PowerDensity
  []
  [modtemp_up]
    type = MultiAppPostprocessorInterpolationTransfer
    from_multi_app = micro
    num_points = 6
    postprocessor = avg_graphtemp
    power = 2
    variable = temp_ms
  []
  [graintemp_up]
    type = MultiAppPostprocessorInterpolationTransfer
    from_multi_app = micro
    num_points = 6
    postprocessor = avg_graintemp
    power = 2
    variable = temp_fg
  []
[]
(htgr/leu_pulse/refcube.i)

TransportSystems

This block is almost identical to the TransportSystems block in the Transient Pulse Initial conditions file, with the exception being that equation_type is set to transient.

Variables

This block sets the temperature variable of the governing equation.

[Variables]
  [temperature]
    family = LAGRANGE
    initial_condition = 300.0
    order = FIRST
    scaling = 1e-8
  []
[]
(htgr/leu_pulse/refcube.i)

Kernels

This block defines the heat conduction governing equation.

[Kernels]
  [HeatConduction]
    type = HeatConduction
    variable = temperature
  []
  [HeatStorage]
    type = HeatCapacityConductionTimeDerivative
    variable = temperature
  []
  [HeatSource]
    block = 10
    type = CoupledForce
    v = PowerDensity
    variable = temperature
  []
[]
(htgr/leu_pulse/refcube.i)

Functions

This block is used to insert reactivity to start the pulse from the initial conditions.

[Functions]
  [boron_state]
    type = PiecewiseLinear
    x = '0.0 0.005 30'
    y = '1.89489259748 2.0 2.0'
  []
[]
(htgr/leu_pulse/refcube.i)

AuxVariables

This block sets up the power density, integrated power, and other auxilliary variables

[AuxVariables]
  [Boron_Conc]
    family = MONOMIAL
    initial_condition = 1.89489259748
    order = CONSTANT
  []
  [PowerDensity]
    block = '10'
    family = MONOMIAL
    order = CONSTANT
  []
  [IntegralPower]
    block = 10
    family = MONOMIAL
    initial_condition = 0.0
    order = CONSTANT
  []
  [avg_coretemp]
    block = 0
    family = LAGRANGE
    initial_condition = 300.0
    order = FIRST
  []
  [temp_fg]
    #Fuel Grain
    block = 10
    family = LAGRANGE
    initial_condition = 300.0
    order = FIRST
  []
  [temp_ms]
    #Moderator Shell
    block = 10
    family = LAGRANGE
    initial_condition = 300.0
    order = FIRST
  []
[]
(htgr/leu_pulse/refcube.i)

AuxKernels

This block defines the equations used to calculate auxilliary values for use in the simulation.

[AuxKernels]
  [pulse_boron]
    type = FunctionAux
    execute_on = timestep_end
    function = boron_state
    variable = Boron_Conc
  []
  [PowerDensityCalc]
    type = VectorReactionRate
    block = '10'
    cross_section = kappa_sigma_fission
    dummies = UnscaledTotalPower
    execute_on = 'initial linear'
    scalar_flux = 'sflux_g0 sflux_g1 sflux_g2 sflux_g3 sflux_g4 sflux_g5'
    scale_factor = PowerScaling
    variable = PowerDensity
  []
  [Powerintegrator]
    type = VariableTimeIntegrationAux
    block = 10
    execute_on = timestep_end
    variable = IntegralPower
    variable_to_integrate = PowerDensity
  []
  [Set_coreT]
    type = SetAuxByPostprocessor
    block = 0
    execute_on = 'linear timestep_end'
    postproc_value = avg_coretemp
    variable = avg_coretemp
  []
[]
(htgr/leu_pulse/refcube.i)

Postprocessors

This block calculates outputs from solution data.

[Postprocessors]
  [UnscaledTotalPower]
    type = Receiver
    outputs = none
  []
  [PowerScaling]
    type = Receiver
    outputs = none
  []
  [avg_coretemp]
    type = ElementAverageValue
    block = 10
    execute_on = linear
    outputs = all
    variable = temp_ms
  []
  [avg_refltemp]
    type = ElementAverageValue
    block = 0
    execute_on = linear
    outputs = all
    variable = temperature
  []
  [max_tempms]
    type = ElementExtremeValue
    block = '10'
    outputs = all
    value_type = max
    variable = temp_ms
  []
  [max_tempfg]
    type = ElementExtremeValue
    block = '10'
    outputs = all
    value_type = max
    variable = temp_fg
  []
  [avg_ms_temp]
    type = ElementAverageValue
    block = 10
    outputs = all
    variable = temp_ms
  []
  [avg_fg_temp]
    type = ElementAverageValue
    block = 10
    outputs = all
    variable = temp_fg
  []
  [avg_powerden]
    type = ElementAverageValue
    block = 10
    execute_on = timestep_end
    outputs = all
    variable = PowerDensity
  []
  [peak_avgPD]
    type = TimeExtremeValue
    execute_on = timestep_end
    postprocessor = avg_powerden
    value_type = max
  []
  [powden_ratio]
    type = PostprocessorRatio
    denominator = peak_avgPD
    execute_on = timestep_end
    numerator = avg_powerden
  []
  [der1_avgPD]
    type = ElementAverageTimeDerivative
    block = 10
    execute_on = timestep_end
    variable = PowerDensity
  []
  [peak_der1avgPD]
    type = TimeExtremeValue
    execute_on = timestep_end
    postprocessor = der1_avgPD
    value_type = max
  []
  [der1PD_ratio]
    type = PostprocessorRatio
    denominator = peak_der1avgPD
    execute_on = timestep_end
    numerator = der1_avgPD
  []
  [ScaledTotalPower]
    type = ElementIntegralVariablePostprocessor
    block = 10
    execute_on = linear
    variable = PowerDensity
  []
  [IntegratedPower]
    type = ElementIntegralVariablePostprocessor
    block = 10
    execute_on = timestep_end
    variable = IntegralPower
  []
  [delta_time]
    type = TimestepSize
  []
  [nl_steps]
    type = NumNonlinearIterations
  []
  [lin_steps]
    type = NumLinearIterations
  []
  [Eq_TREAT_Power]
    type = ScalePostprocessor
    scaling_factor = 2469860.77609
    value = avg_powerden
  []
[]
(htgr/leu_pulse/refcube.i)

UserObjects

The Terminator user object calculates when the simulation concludes.

[UserObjects]
  [der_pulse_end]
    type = Terminator
    execute_on = timestep_end
    expression = '(powden_ratio < 0.25) & (abs(der1PD_ratio) < 0.04)'
  []
[]
(htgr/leu_pulse/refcube.i)

Materials

This block contains material properties used in the governing equation.

[Materials]
  [neut_mix]
    type = CoupledFeedbackNeutronicsMaterial
    block = 10
    densities = '0.998448391539 0.00155160846058'
    grid_names = 'Tfuel Tmod Rod'
    grid_variables = 'temp_fg temp_ms Boron_Conc'
    isotopes = 'pseudo1 pseudo2'
    library_file = 'cross_sections/leu_20r_is_6g_d.xml'
    library_name = leu_20r_is_6g_d
    material_id = 1
    plus = true
  []
  [kth]
    # Volume weighted harmonic mean
    # Divided fg_kth by 100 to get it into cm
    type = ParsedMaterial
    coupled_variables = 'temp_fg'
    block = 10
    constant_expressions = '3.35103216383e-08 1.31125888571e-07 2.14325144175e-05 0.3014 0.01046 1.0 0.05 1.5 1.0'
    constant_names = 'vol_fg vol_fl vol_gr gr_kth fl_kth beta p_vol sigma kap3x'
    property_name = 'thermal_conductivity'
    expression = 'lt := temp_fg / 1000.0; fresh := (100.0 / (6.548 + 23.533 * lt) + 6400.0 * exp(-16.35 / lt) / pow(lt, 5.0/2.0)) / 100.0; kap1d := (1.09 / pow(beta, 3.265) + 0.0643 * sqrt(temp_fg) / sqrt(beta)) * atan(1.0 / (1.09 / pow(beta, 3.265) + sqrt(temp_fg) * 0.0643 / sqrt(beta))); kap1p := 1.0 + 0.019 * beta / ((3.0 - 0.019 * beta) * (1.0 + exp(-(temp_fg - 1200.0) / 100.0))); kap2p := (1.0 - p_vol) / (1.0 + (sigma - 1.0) * p_vol); kap4r := 1.0 - 0.2 / (1.0 + exp((temp_fg - 900.0) / 80.0)); fg_kth := fresh * kap1d * kap1p * kap2p * kap3x * kap4r; (vol_fg + vol_fl + vol_gr) / (vol_fg / fg_kth + vol_fl / fl_kth + vol_gr / gr_kth)'
  []
  [rho_cp]
    # Volume weighted arithmetic mean (Irradiation has no effect)
    type = ParsedMaterial
    coupled_variables = 'temp_fg temp_ms'
    block = 10
    constant_expressions = '3.35103216383e-08 2.1563640306e-05 0.0018 0.010963'
    constant_names = 'vol_fg vol_gr rho_gr rho_fg'
    property_name = 'heat_capacity'
    expression = 'lt := temp_fg / 1000.0; gr_rhocp := rho_gr / (11.07 * pow(temp_ms, -1.644) + 0.0003688 * pow(temp_ms, 0.02191)); fink_cp := 52.1743 + 87.951 * lt - 84.2411 * pow(lt, 2) + 31.542 * pow(lt, 3) - 2.6334 * pow(lt, 4) - 0.71391 * pow(lt, -2); fg_rhocp := rho_fg * fink_cp / 267.2 * 1000.0; (vol_fg * fg_rhocp + vol_gr * gr_rhocp) / (vol_fg + vol_gr)'
  []
  [neut_refl]
    type = CoupledFeedbackNeutronicsMaterial
    block = 0
    densities = '1'
    grid_names = 'Trefl Tcore Rod'
    grid_variables = 'temperature avg_coretemp Boron_Conc'
    isotopes = 'pseudo'
    library_file = 'cross_sections/leu_macro_6g.xml'
    library_name = leu_macro_6g
    material_id = 2
    plus = true
  []
  [ref_kth]
    type = GenericConstantMaterial
    block = 0
    prop_names = 'thermal_conductivity'
    prop_values = '0.3014'
  []
  [ref_rho_cp]
    type = ParsedMaterial
    coupled_variables = 'temperature'
    block = '0'
    constant_expressions = '0.0018'
    constant_names = 'rho_gr'
    property_name = 'heat_capacity'
    expression = 'rho_gr / (11.07 * pow(temperature, -1.644) + 0.0003688 * pow(temperature, 0.02191))'
  []
[]
(htgr/leu_pulse/refcube.i)

Preconditioning

This block is identical to the [Preconditioning] block in the Transient Pulse Initial Conditions file.

Executioner

This block describes the simulation solution approach.

[Executioner]
  type = IQS
  do_iqs_transient = false
  pke_param_csv = true

  # Time evolution parameters
  dtmin = 1e-7
  start_time = 0.0
  end_time = 10.0
  [TimeStepper]
    type = ConstantDT
    dt = 0.005
    growth_factor = 1.5
  []

  # Solver parameters
  l_max_its = 100
  l_tol = 1e-3
  nl_abs_tol = 1e-7
  nl_max_its = 200
  nl_rel_tol = 1e-7

  # Multiphysics iteration parameters
  fixed_point_abs_tol = 1e-7
  fixed_point_max_its = 10
  fixed_point_rel_tol = 1e-7
[]
(htgr/leu_pulse/refcube.i)

Run Command

Use this command to run the main input file.


mpirun -np 48 /path/to/griffin-opt -i refcube.i

Output Files

When the main input file is run, many output files will be generated.

Most of these have the format out~refcube_micro_<number> as either CSV (.csv) or Exodus (.e) files. These files are output by the main input file. These files contain the postprocessor values specified in the [Postprocessor] block for each of the fuel grains. The main input file will also output a CSV file containing the PKE parameters, called out~refcube_pke_params.csv.

Each of the four input files will output a CSV file and an Exodus file in the format <input_file_name>_out.csv or <input_file_name>_exodus.e. These files will contain the postprocessor values specified in the corresponding [Postprocessor] block.

References

  1. Adam Zabriskie, Sebastian Schunert, Daniel Schwen, Javier Ortensi, Benjamin Baker, Yaqi Wang, Vincent Laboure, Mark DeHart, and Wade Marcum. A coupled multiscale approach to treat leu feedback modeling using a binary-collision monte-carlo–informed heat source. Nuclear Science and Engineering, 193(4):368–387, 2019. URL: https://doi.org/10.1080/00295639.2018.1528802, arXiv:https://doi.org/10.1080/00295639.2018.1528802, doi:10.1080/00295639.2018.1528802.[BibTeX]
  2. Adam X. Zabriskie. Multi-Scale, Multi-Physics Reactor Pulse Simulation Method with Macroscopic and Microscopic Feedback Effects. PhD thesis, Oregon State University, Corvallis OR, 2019.[BibTeX]