UO Pebbles

Heat Conduction Model

The MOOSE heat transfer module is used to solve for steady-state heat conduction,

(1)

where is the solid thermal conductivity, is the solid temperature, and is a volumetric heat source in the solid.

The solid mesh is shown in Figure 1. The pebble surface is sideset 0. The MOOSE solid problem is set up with a length unit of meters, which is convenient for heat transfer applications because material properties for thermal-fluids physics are almost always given in SI units. Further, many MOOSE physics modules inherently assume SI units, such as the fluid property module and solid property module.

Figure 1: Mesh for solid domain

Because heat transfer and fluid flow in the FLiBe is not modeled in this example, heat removal by the fluid is approximated by setting the outer surface of the pebble to a convection boundary condition,

(2)

where W/m/K and is set to a function linearly ranging from 650°C at cm to 750°C at cm. In this example, the MOOSE heat transfer module will be run first. The initial solid temperature is set to 650°C and the heat source to zero.

OpenMC Model

The OpenMC model is built using CSG; we will leverage the lattice feature in OpenMC to repeat a universe in a structured manner throughout the domain. Here, the universe consists of a single pebble and surrounding FLiBe, which is repeated three times throughout the geometry. The overall domain is enclosed in a box with - width of 44 cm and height of 12 cm. All boundaries of this box are reflecting.

OpenMC's Python API is used to create the pebbles model with the script shown below. First, we define materials for the various regions and create the geometry. We create a universe consisting of a single pebble, and then repeat that universe three times in a lattice. Because we have one cell per pebble, each pebble will receive a single temperature from MOOSE. The cell level is 1 because the pebble lattice is nested once with respect to the highest level, and we want to apply feedback to each pebble individually. If we were to specify level zero in the geometry, we would apply feedback to the main_cell, which would apply temperature feedback across the problem globally. The OpenMC geometry as produced via plots is shown in Figure 2.

import openmc
import numpy as np

# Add materials for the pebbles and the FLiBe

uo2 = openmc.Material()
uo2.set_density('g/cc', 10.0)
uo2.add_element('U', 1.0, enrichment=5.0)
uo2.add_element('O', 2.0)

enrichment_li7 = 0.99995
flibe = openmc.Material(name='2LiF-BeF2')
flibe.set_density('kg/m3', 1960)
flibe.add_nuclide('Li7', 2.0*enrichment_li7)
flibe.add_nuclide('Li6', 2.0*(1 - enrichment_li7))
flibe.add_element('Be', 1.0)
flibe.add_element('F', 4.0)

mats = openmc.Materials([uo2, flibe])
mats.export_to_xml()

model = openmc.model.Model()

# create a universe containing the repeatable universe of
# a single pebble plus surrounding flibe
R = 1.5
L = 4.0
sphere_surface = openmc.Sphere(r=R, name='Sphere surface')
flibe_surface = openmc.model.RectangularPrism(L, L, boundary_type='reflective')

sphere_cell = openmc.Cell(fill=uo2, region=-sphere_surface, name='Pebble')
flibe_cell = openmc.Cell(fill=flibe, region=+sphere_surface & -flibe_surface, name='Flibe')
repeatable_univ = openmc.Universe(cells=[sphere_cell, flibe_cell])

outer_cell = openmc.Cell(fill=flibe, region=+sphere_surface & -flibe_surface, name='Outside')
outer_univ = openmc.Universe(cells=[outer_cell])

# create a lattice to repeat the pebble + flibe universe
N = 3
lattice = openmc.RectLattice()
lattice.lower_left = (-L/2.0, -L/2.0, 0)
lattice.pitch = (L, L, L)
lattice.universes = np.full((N, 1, 1), repeatable_univ)
lattice.outer = outer_univ

# create the surfaces that will bound the lattice
top = openmc.ZPlane(z0=N*L, boundary_type='reflective')
bottom = openmc.ZPlane(z0=0.0, boundary_type='reflective')
main_cell = openmc.Cell(fill=lattice, region=-flibe_surface & +bottom & -top, name='Main cell')

model.geometry = openmc.Geometry([main_cell])
model.geometry.export_to_xml()

# Finally, define some run settings
model.settings = openmc.Settings()
model.settings.batches = 1000
model.settings.inactive = 100
model.settings.particles = 10000

lower_left = (-L, -L, 0)
upper_right = (L, L, N*L)
uniform_dist = openmc.stats.Box(lower_left, upper_right, only_fissionable=True)
model.settings.source = openmc.IndependentSource(space=uniform_dist)
model.settings.temperature = {'default': 650.0 + 273.15,
                              'method': 'interpolation',
                              'multipole': False,
                              'tolerance': 1000.0,
                              'range': (294.0, 1600.0)}
model.settings.export_to_xml()

# create some plots for making images in the tutorial

height = N*L
plot1          = openmc.Plot()
plot1.filename = 'plot1'
plot1.width    = (L, height)
plot1.basis    = 'xz'
plot1.origin   = (0.0, 0.0, height/2.0)
plot1.pixels   = (100, 300)
plot1.color_by = 'cell'

plot2          = openmc.Plot()
plot2.filename = 'plot2'
plot2.width    = (L, L)
plot2.basis    = 'xy'
plot2.origin   = (0.0, 0.0, 2.0)
plot2.pixels   = (100, 100)
plot2.color_by = 'cell'

plots = openmc.Plots([plot1, plot2])
plots.export_to_xml()

Figure 2: OpenMC geometry (colored by cell ID)

As you can see, there are only two unique cell IDs. This problem is built using distributed cells, meaning a cell with the same ID appears multiple times in the geometry; each time the cell is repeated, we assign a new instance to that cell. So, the three pebbles are represented as:

• Cell ID 1, instance 0

• Cell ID 1, instance 1

• Cell ID 1, instance 2

and the FLiBe region is represented as:

• Cell ID 2, instance 0

• Cell ID 2, instance 1

• Cell ID 2, instance 2

Because OpenMC runs after the MOOSE heat transfer module, initial conditions are only required for the FLiBe temperature, which is set to 650°C. Because there is no density feedback in this example, the densities initially imposed in the OpenMC model remain fixed at the values set in the OpenMC input files.

To create the XML files required to run OpenMC, run the script:

python make_openmc_model.py

Solid Input Files

The solid phase is solved with the MOOSE heat transfer module, and is described in the solid.i input. We set up the mesh using the CombinerGenerator to translate a single sphere mesh to multiple locations.

[Mesh<<<{"href": "../syntax/Mesh/index.html"}>>>]
  [pebble]
    type = SphereMeshGenerator<<<{"description": "Generate a 3-D sphere mesh centered on the origin", "href": "../source/meshgenerators/SphereMeshGenerator.html"}>>>
    nr<<<{"description": "Number of radial elements"}>>> = 2
    radius<<<{"description": "Sphere radius"}>>> = 0.015
  []
  [repeat]
    type = CombinerGenerator<<<{"description": "Combine multiple meshes (or copies of one mesh) together into one (disjoint) mesh.  Can optionally translate those meshes before combining them.", "href": "../source/meshgenerators/CombinerGenerator.html"}>>>
    inputs<<<{"description": "The input MeshGenerators.  This can either be N generators or 1 generator.  If only 1 is given then 'positions' must also be given."}>>> = pebble
    positions<<<{"description": "The (optional) position of each given mesh.  If N 'inputs' were given then this must either be left blank or N positions must be given.  If 1 input was given then this MUST be provided."}>>> = '0 0 0.02
                 0 0 0.06
                 0 0 0.10'
  []
[]

Next, we define the temperature variable, temp, and specify the governing equations and boundary conditions that we will apply. We set the thermal conductivity of the pebbles to 50 W/m/K.

[Variables<<<{"href": "../syntax/Variables/index.html"}>>>]
  [temp]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = ${T_fluid}
  []
[]

[Kernels<<<{"href": "../syntax/Kernels/index.html"}>>>]
  [hc]
    type = HeatConduction<<<{"description": "Diffusive heat conduction term $-\\nabla\\cdot(k\\nabla T)$ of the thermal energy conservation equation", "href": "../source/kernels/HeatConduction.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
  []
  [heat]
    type = CoupledForce<<<{"description": "Implements a source term proportional to the value of a coupled variable. Weak form: $(\\psi_i, -\\sigma v)$.", "href": "../source/kernels/CoupledForce.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
    v<<<{"description": "The coupled variable which provides the force"}>>> = heat_source
  []
[]

[Functions<<<{"href": "../syntax/Functions/index.html"}>>>]
  [T_fluid]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = '${T_fluid}+z*1000'
  []
[]

[BCs<<<{"href": "../syntax/BCs/index.html"}>>>]
  [surface]
    type = ConvectiveFluxFunction<<<{"description": "Determines boundary value by fluid heat transfer coefficient and far-field temperature", "href": "../source/bcs/ConvectiveFluxFunction.html"}>>>
    T_infinity<<<{"description": "Function describing far-field temperature"}>>> = T_fluid
    coefficient<<<{"description": "Function describing heat transfer coefficient"}>>> = 1000.0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = '0'
  []
[]

[Materials<<<{"href": "../syntax/Materials/index.html"}>>>]
  [k]
    type = GenericConstantMaterial<<<{"description": "Declares material properties based on names and values prescribed by input parameters.", "href": "../source/materials/GenericConstantMaterial.html"}>>>
    prop_values<<<{"description": "The values associated with the named properties"}>>> = '50.0'
    prop_names<<<{"description": "The names of the properties this material will have"}>>> = 'thermal_conductivity'
  []
[]

The heat source received from OpenMC is stored in the heat_source auxiliary variable.

[AuxVariables<<<{"href": "../syntax/AuxVariables/index.html"}>>>]
  [heat_source]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
[]

Finally, we specify that we will run a OpenMC as a sub-application, with data transfers of heat source from OpenMC and temperature to OpenMC and use a transient executioner.

[MultiApps<<<{"href": "../syntax/MultiApps/index.html"}>>>]
  [openmc]
    type = TransientMultiApp<<<{"description": "MultiApp for performing coupled simulations with the parent and sub-application both progressing in time.", "href": "../source/multiapps/TransientMultiApp.html"}>>>
    input_files<<<{"description": "The input file for each App.  If this parameter only contains one input file it will be used for all of the Apps.  When using 'positions_from_file' it is also admissable to provide one input_file per file."}>>> = 'openmc.i'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
[]

[Transfers<<<{"href": "../syntax/Transfers/index.html"}>>>]
  [heat_source_from_openmc]
    type = MultiAppGeneralFieldNearestLocationTransfer<<<{"description": "Transfers field data at the MultiApp position by finding the value at the nearest neighbor(s) in the origin application.", "href": "../source/transfers/MultiAppGeneralFieldNearestLocationTransfer.html"}>>>
    from_multi_app<<<{"description": "The name of the MultiApp to receive data from"}>>> = openmc
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = heat_source
    source_variable<<<{"description": "The variable to transfer from."}>>> = heat_source
    from_postprocessors_to_be_preserved<<<{"description": "The name of the Postprocessor in the from-app to evaluate an adjusting factor."}>>> = heat_source
    to_postprocessors_to_be_preserved<<<{"description": "The name of the Postprocessor in the to-app to evaluate an adjusting factor."}>>> = source_integral
  []
  [temp_to_openmc]
    type = MultiAppGeneralFieldShapeEvaluationTransfer<<<{"description": "Transfers field data at the MultiApp position using the finite element shape functions from the origin application.", "href": "../source/transfers/MultiAppGeneralFieldShapeEvaluationTransfer.html"}>>>
    to_multi_app<<<{"description": "The name of the MultiApp to transfer the data to"}>>> = openmc
    variable<<<{"description": "The auxiliary variable to store the transferred values in."}>>> = temp
    source_variable<<<{"description": "The variable to transfer from."}>>> = temp
  []
[]

[Postprocessors<<<{"href": "../syntax/Postprocessors/index.html"}>>>]
  [source_integral]
    type = ElementIntegralVariablePostprocessor<<<{"description": "Computes a volume integral of the specified variable", "href": "../source/postprocessors/ElementIntegralVariablePostprocessor.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = heat_source
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = transfer
  []
  [max_T]
    type = NodalExtremeValue<<<{"description": "Finds either the min or max elemental value of a variable over the domain.", "href": "../source/postprocessors/NodalExtremeValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on"}>>> = temp
  []
[]

[Executioner<<<{"href": "../syntax/Executioner/index.html"}>>>]
  type = Transient
  petsc_options_iname = '-pc_type -pc_hypre_type'
  num_steps = 3
  petsc_options_value = 'hypre boomeramg'
  nl_abs_tol = 1e-10
[]

[Outputs<<<{"href": "../syntax/Outputs/index.html"}>>>]
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
[]

Neutronics Input Files

The neutronics physics is solved over the entire domain using OpenMC. The OpenMC wrapping is described in the openmc.i input file. We begin by defining a mesh mirror on which OpenMC will receive temperature from the coupled MOOSE application, and on which OpenMC will write the fission heat source.

[Mesh<<<{"href": "../syntax/Mesh/index.html"}>>>]
  [pebble]
    type = SphereMeshGenerator<<<{"description": "Generate a 3-D sphere mesh centered on the origin", "href": "../source/meshgenerators/SphereMeshGenerator.html"}>>>
    nr<<<{"description": "Number of radial elements"}>>> = 2
    radius<<<{"description": "Sphere radius"}>>> = 0.015
  []
  [repeat]
    type = CombinerGenerator<<<{"description": "Combine multiple meshes (or copies of one mesh) together into one (disjoint) mesh.  Can optionally translate those meshes before combining them.", "href": "../source/meshgenerators/CombinerGenerator.html"}>>>
    inputs<<<{"description": "The input MeshGenerators.  This can either be N generators or 1 generator.  If only 1 is given then 'positions' must also be given."}>>> = pebble
    positions<<<{"description": "The (optional) position of each given mesh.  If N 'inputs' were given then this must either be left blank or N positions must be given.  If 1 input was given then this MUST be provided."}>>> = '0 0 0.02
                 0 0 0.06
                 0 0 0.10'
  []
[]

Cardinal will also automatically output a variable named cell_id (CellIDAux) and a variable named cell_instance ( CellInstanceAux) to show the spatial mapping. We also define a CellTemperatureAux to show the OpenMC volume-averaged temperatures as they map to the [Mesh].

[AuxVariables<<<{"href": "../syntax/AuxVariables/index.html"}>>>]
  [cell_temperature]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
  []
[]

[AuxKernels<<<{"href": "../syntax/AuxKernels/index.html"}>>>]
  [cell_temperature]
    type = CellTemperatureAux<<<{"description": "OpenMC cell temperature (K), mapped to each MOOSE element", "href": "../source/auxkernels/CellTemperatureAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = cell_temperature
  []
[]

The [Problem] and [Tallies] blocks are then used to specify the OpenMC wrapping. We define a total power of 1500 W, and indicate that we'd like to add a CellTally on block 0, which corresponds to the pebbles. The cell tally setup in Cardinal will then automatically add a tally for each unique cell ID+instance combination. Because the repeated pebble cells we'd like to tally are repeated in the lattice nested one level below the root universe, we set the cell_level = 1.

[Problem<<<{"href": "../syntax/Problem/index.html"}>>>]
  type = OpenMCCellAverageProblem
  verbose = true
  power = 1500.0
  temperature_blocks = '0'
  cell_level = 1
  scaling = 100.0

  [Tallies<<<{"href": "../syntax/Problem/Tallies/index.html"}>>>]
    [heat_source]
      type = CellTally<<<{"description": "A class which implements distributed cell tallies.", "href": "../source/tallies/CellTally.html"}>>>
      block<<<{"description": "Subdomains for which to add tallies in OpenMC. If not provided, cell tallies will be applied over the entire mesh."}>>> = '0'
      name<<<{"description": "Auxiliary variable name(s) to use for OpenMC tallies. If not specified, defaults to the names of the scores"}>>> = heat_source
    []
  []
[]

The scaling parameter is used to indicate a multiplicative factor that should be applied to the [Mesh] in order to get to units of centimeters. Because the [Mesh] is in units of meters, we set scaling = 100.0. This scaling factor is applied within the OpenMCCellAverageProblem::findCell routine that maps MOOSE elements to OpenMC cells - no actual changes are made to the mesh in the [Mesh] block. The scaling is also applied to ensure that the heat source is on the correct per-unit-volume scale that is expected by the [Mesh].

Finally, we set a transient executioner, specify an Exodus output, and define several postprocessors.

[Executioner<<<{"href": "../syntax/Executioner/index.html"}>>>]
  type = Transient
[]

[Outputs<<<{"href": "../syntax/Outputs/index.html"}>>>]
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
[]

[Postprocessors<<<{"href": "../syntax/Postprocessors/index.html"}>>>]
  [heat_source]
    type = ElementIntegralVariablePostprocessor<<<{"description": "Computes a volume integral of the specified variable", "href": "../source/postprocessors/ElementIntegralVariablePostprocessor.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = heat_source
  []
  [max_tally_rel_err]
    type = TallyRelativeError<<<{"description": "Maximum/minimum tally relative error", "href": "../source/postprocessors/TallyRelativeError.html"}>>>
  []
  [k]
    type = KEigenvalue<<<{"description": "k eigenvalue computed by OpenMC", "href": "../source/postprocessors/KEigenvalue.html"}>>>
  []
[]