Doppler Slab Analytic Benchmark

Problem Description

This test case models the 1D Doppler slab benchmark from Remley and Griesheimer (2019), which couples neutron transport and heat conduction in a semi-infinite slab with microscopic absorption cross section , number density , and thermal conductivity . A beam of neutrons (intensity ) is incident on the left face () of the slab, which is held at a fixed temperature of . The geometry of the benchmark problem can be found in Figure 1.

Doppler slab geometry, taken from [!cite](doppler_slab_benchmark).

Figure 1: Doppler slab geometry, taken from Remley and Griesheimer (2019).

The macroscopic absorption cross section in the slab takes one of two functional forms. The first is a linear dependence on temperature:

(1)

The second is an inverse square-root dependence on temperature:

(2)

Absorption reactions heat the slab with an energy released per reaction of . The analytic solution for the temperature (Eq. (3)) and flux (Eq. (4)) in the slab when the cross section is represented by Eq. (1) can be obtained:

(3)(4)

is defined as:

(5)

and is:

(6)

For the case where Eq. (2) is used for the cross section, the following solutions for the temperature (Eq. (7)) and flux (Eq. (8)) can be obtained:

(7)

and

(8)

is defined as:

(9)

and is:

(10)

is the Lambert W function (also known as the product logarithm). When using linear temperature feedback, the benchmark parameters must be selected such that to prevent non- physical solutions (negative cross sections). There are no such constraints on the inverse-root feedback. The cannonical benchmark parameters can be found in Table 1.

Table 1: Doppler slab benchmark parameters proposed in Remley and Griesheimer (2019).

ParameterValue
(n cm-2 s-1)
(atom barn-1 cm-1)
(barn)
(eV)
(K)
0.6 W (m-1 K-1)
(cm-1 K-1)

Cardinal Multi-Group Monte Carlo Solutions

Cardinal's coupling between MOOSE's heat conduction solver, OpenMC, and NekRS running as a solid heat conduction solver was verified in Hegazy and Novak (2023) with the Doppler slab problem. 200 mesh elements with linear Lagrange basis functions were used in the heat conduction solution. The OpenMC simulation used 200 cells and a mesh tally with 200 bins to ensure temperature feedback and flux tallies matched the heat conduction mesh. As the benchmark problem is semi-infinite, a length must be chosen to truncate the problem at in order to apply a symmetry boundary condition (heat conduction) and vacuum condition (neutronics). This was selected to be cm. Results for Eq. (1) can be found in Figure 2, and results for Eq. (2) can be found in Figure 3.

Cardinal simulations of the Doppler slab problem with linear temperature feedback, taken from [!cite](doppler_slab_mc_cardinal).

Figure 2: Cardinal simulations of the Doppler slab problem with linear temperature feedback, taken from Hegazy and Novak (2023).

Cardinal simulations of the Doppler slab problem with inverse root temperature feedback, taken from [!cite](doppler_slab_mc_cardinal).

Figure 3: Cardinal simulations of the Doppler slab problem with inverse root temperature feedback, taken from Hegazy and Novak (2023).

The maximum relative error for linear feedback is (temperature) and (flux). The maximim relative error for inverse root feedback is (temperature) and (flux). Error in Cardinal results is dominated by (i) statistical noise near the right boundary due to neutron attenuation, and (ii) boundary condition effects at the right boundary caused by the finite (as opposed to semi-infinite) slab. More information regarding the problem setup and results can be found in Hegazy and Novak (2023).

Input files for the linear temperature feedback case are included in Cardinal to verify the heat conduction module, Cardinal's OpenMC coupling, and different temperature interpolation schemes present in OpenMC. These input files can be found in /test/tests/neutronics/mg/doppler_slab_lin. The OpenMC and heat conduction inpute files for this problem can be found below.

joule_per_ev = 1.60218e-19

q      = ${fparse 1.0e6}       # Energy release per absorption (eV)
Sigma0 = ${fparse 4.0 * 0.025} # Initial macroscopic XS, 1/cm
T0     = 293.6                 # Surface temperature (K)
Y0     = 6.65e11               # Source intensity (n / cm2-s)
alpha  = -0.0001               # Linear doppler Coefficient (1 / K)
k      = 0.006                 # Thermal conductivity (W / cm-K)

!include mesh.i

[AuxVariables]
  [cell_temperature]
    family = MONOMIAL
    order = CONSTANT
  []

  # Variables added to convert into proper heat source term for the power in the solid
  [heat_source] # W / cm^3
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [cell_temperature]
    type = CellTemperatureAux
    variable = cell_temperature
  []

  # Change the unit of the 'flux' (neutrons / m^2 / s) into a volumetric power. This is
  # based on Eq. (5b) in the paper, which shows the volumetric power is q / k * Sigma * flux
  [compute_power]
    type = ParsedAux
    variable = heat_source
    coupled_variables = 'flux temp'
    expression = 'flux * ${q} * ${joule_per_ev} * ${Sigma0} * (1+(${alpha}/${Sigma0})*(temp-${T0}))'
    execute_on = 'timestep_end'
  []
[]

[ICs]
  [temp]
    type = ConstantIC
    variable = temp
    value = ${T0}
  []
[]

[Problem]
  type = OpenMCCellAverageProblem
  verbose = true
  source_strength = ${fparse Y0 * 2.0 * 2.0} # multiply by the left area to get units of neutrons / s

  cell_level = 0
  temperature_blocks = '0'
  initial_properties = xml

  particles = 1000
  batches = 100
  inactive_batches = 0

  relaxation = 'constant'

  [Tallies]
    [heating]
      type = CellTally
      score = 'flux'
      output = 'unrelaxed_tally_std_dev'
      check_tally_sum = false
    []
  []
[]

[Executioner]
  type = Transient
  dt = 1.0
  steady_state_detection = true
  check_aux = true
  num_steps = 3
[]

[MultiApps]
  [solid]
    type = TransientMultiApp
    app_type = CardinalApp
    input_files = 'solid.i'
    execute_on = timestep_end
    sub_cycling = true
  []
[]

[Transfers]
  [solid_temp]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = T
    variable = temp
    from_multi_app = solid
  []
  [source_to_solid]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    source_variable = heat_source
    variable = power
    to_multi_app = solid
  []
[]

[Postprocessors]
  [max_heat_source]
    type = ElementExtremeValue
    variable = heat_source
    execute_on = 'timestep_begin'
  []
  [max_flux]
    type = ElementExtremeValue
    variable = flux
  []
  [max_temp]
    type = ElementExtremeValue
    variable = temp
  []
[]

[Outputs]
  csv = true
  execute_on = 'TIMESTEP_END'
[]
(test/tests/neutronics/mg/doppler_slab_lin/openmc_base.i)
!include mesh.i

[Variables]
  [T]
    initial_condition = 293.6
  []
[]

[Kernels]
  [diffusion]
    type = HeatConduction
    variable = T
    block = '0'
  []
  [source]
    type = CoupledForce
    variable = T
    v = power
    block = '0'
  []
[]

[BCs]
  [left_boundary]
    type = DirichletBC
    variable = T
    boundary = 'left'
    value = 293.6
  []
  [insulated]
    type = NeumannBC
    variable = T
    boundary = 'right top bottom front back'
    value = 0
  []
[]

[Materials]
  [absorbing]
    type = HeatConductionMaterial
    thermal_conductivity = ${fparse 0.6 / 100.0} # W/cm-K
    block = '0'
  []
[]

[Postprocessors]
  [T_solid]
    type = ElementExtremeValue
    variable = T
  []
[]

[AuxVariables]
  [power]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[Executioner]
  type = Transient
  nl_abs_tol = 1e-5
  nl_rel_tol = 1e-16
  petsc_options_value = 'hypre boomeramg'
  petsc_options_iname = '-pc_type -pc_hypre_type'
[]
(test/tests/neutronics/mg/doppler_slab_lin/solid.i)

References

  1. Aya H. Hegazy and April J. Novak. Verification of the Cardinal Multiphysics Solver with the Doppler Slab Benchmark. In Proceedings of M&C. 2023.[BibTeX]
  2. Kyle E. Remley and David P. Griesheimer. A Fully Analytic Coupled Thermal-Neutronics Benchmark and its Application to Monte Carlo Simulation. In Proceedings of M&C. 2019.[BibTeX]