LWR Pincell

MOOSE Heat Conduction Model

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

(2)

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

MeshGenerators are used to construct the solid mesh. Figure 1 shows the solid mesh with block IDs and sidesets. Different block IDs are used for the hexahedral and prism elements in the pellet region because libMesh does not allow different element types to exist on the same block ID. Because this mesh is generated using the MeshGenerator system in MOOSE, the mesh is created at runtime. If you want to generate a mesh file, you can do so by running the solid input file in mesh generation mode:

cardinal-opt -i mesh.i --mesh-only

Figure 1: Mesh for the solid portions of an LWR pincell; boundary IDs are shown to the right

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

(3)

where W/m/K and is a function linearly increasing from °C at the bottom of the pincell to °C at the top of the pincell. The top and bottom of the solid pincell are assumed insulated.

The gap region is unmeshed, and a quadrature-based thermal contact model is applied based on the sum of thermal conduction and thermal radiation (across a transparent medium). For a paired set of boundaries, each quadrature point on boundary is paired with the nearest quadrature point on boundary . Then, the sum of the radiation and conduction heat fluxes imposed between quadrature point pairs is

(4)

where is the Stefan-Boltzmann constant, is the temperature at a quadrature point, is the temperature of the nearest quadrature point across the gap, and are emissivities of boundaries and , respectively, and is the conduction resistance. For cylindrical geometries, the conduction resistance is given as

(5)

where are the radial coordinates associated with the outer and inner radii of the cylindrical annulus, is the height of the annulus, and is the thermal conductivity of the annulus material.

In this example, the MOOSE heat transfer module will run first. The initial solid temperature is 280°C and the initial power is zero.

OpenMC Model

The OpenMC model is built using CSG, which are cells created from regions of space formed by half-spaces of various common surfaces. When creating the OpenMC geometry, there are three aspects that you must pay attention to when using Cardinal -

1. The resolution of the temperature (and density, for fluid feedback) feedback to impose in OpenMC

2. The "level" of the cells with which you want to perform feedback

3. The manner in which temperature feedback is applied to the cross section data

The temperature in OpenMC is stored on the cells. One constant temperature can be set for each cell. Therefore, the resolution of the temperature feedback received from MOOSE is determined during the OpenMC model setup. Each OpenMC cell will receive a unique temperature, so the number of OpenMC cells dictates the temperature feedback resolution.

The second consideration is slightly more subtle, but allows great flexibility for imposing multiphysics feedback for very heterogeneous geometries, such as TRISO fuels. The "level" of a cell refers to the number of nested universes (relative to the root universe) at which you would like to impose feedback. If you construct your geometry without filling any OpenMC cells with other universes, then all your cells are at level zero - i.e. the highest level in the geometry. But if your model contains lattices or filled universes, the level that you want to perform multiphysics coupling on is most likely not the highest level in the geometry. For instance, LWR cores are comprised of hundreds of assemblies. The simplest approach to creating this geometry would be to make a single fuel assembly universe, then repeat that universe several times throughout the geometry. If you wanted to apply a single temperature to the entire assembly (pins plus coolant), then you would set a level of 1. If your assembly itself was formed as a lattice of single-pin unit cells, then to apply unique temperatures to each pin, you would set a level of 2. In other words, when we establish a mapping to MOOSE, we do not simply use the lowest cell level in the geometry because it is sometimes desirable to set the temperature/density for all cells contained in a particular cell.

OpenMC's Python API is used to create the pincell model with the script shown below. First, we define materials and create the geometry. We add 40 cells to receive temperature feedback by dividing the entire axial height by 41 axial planes. Note that this particular choice of axial cells has no relationship to the solid mesh. That is, MOOSE elements can span more than one OpenMC cell; when mapping by centroid, however, each MOOSE element will then be associated with one OpenMC cell (which may be an imperfect alignment, but represents a form of discretization error in the MOOSE element refinement). For simplicity here we select the same axial discretization in the OpenMC model as used in the [Mesh]. The OpenMC geometry as produced via plots is shown in Figure 2.

import openmc
import numpy as np
from argparse import ArgumentParser

ap = ArgumentParser()
ap.add_argument('-n', dest='n_axial', type=int, default=40,
                help='Number of axial cell divisions')
ap.add_argument('-s', '--entropy', action='store_true',
                help='Whether to add a Shannon entropy mesh')
args = ap.parse_args()

N = args.n_axial # Number of axial cells to build in the solid to receive feedback
height = 300.0   # Total height of pincell

###############################################################################
# Create materials for the problem

uo2 = openmc.Material(name='UO2 fuel at 2.4% wt enrichment')
uo2.set_density('g/cm3', 10.29769)
uo2.add_element('U', 1., enrichment=2.4)
uo2.add_element('O', 2.)

helium = openmc.Material(name='Helium for gap')
helium.set_density('g/cm3', 0.001598)
helium.add_element('He', 2.4044e-4)

zircaloy = openmc.Material(name='Zircaloy 4')
zircaloy.set_density('g/cm3', 6.55)
zircaloy.add_element('Sn', 0.014  , 'wo')
zircaloy.add_element('Fe', 0.00165, 'wo')
zircaloy.add_element('Cr', 0.001  , 'wo')
zircaloy.add_element('Zr', 0.98335, 'wo')

borated_water = openmc.Material(name='Borated water')
borated_water.set_density('g/cm3', 0.740582)
borated_water.add_element('B', 4.0e-5)
borated_water.add_element('H', 5.0e-2)
borated_water.add_element('O', 2.4e-2)
borated_water.add_s_alpha_beta('c_H_in_H2O')

# Collect the materials together and export to XML
materials = openmc.Materials([uo2, helium, zircaloy, borated_water])
materials.export_to_xml()

###############################################################################
# Define problem geometry

# Create cylindrical surfaces
fuel_or = openmc.ZCylinder(r=0.39218, name='Fuel OR')
clad_ir = openmc.ZCylinder(r=0.40005, name='Clad IR')
clad_or = openmc.ZCylinder(r=0.45720, name='Clad OR')

# Create a region represented as the inside of a rectangular prism
pitch = 1.25984
box = openmc.model.RectangularPrism(pitch, pitch, boundary_type='reflective')

# Create cells, mapping materials to regions - split up the axial height
planes = np.linspace(0.0, height, N + 1)
plane_surfaces = []
for i in range(N + 1):
  plane_surfaces.append(openmc.ZPlane(z0=planes[i]))

# set the boundary condition on the topmost and bottommost planes to vacuum
plane_surfaces[0].boundary_type = 'vacuum'
plane_surfaces[-1].boundary_type = 'vacuum'

fuel_cells = []
clad_cells = []
gap_cells = []
water_cells = []
all_cells = []
for i in range(N):
  layer = +plane_surfaces[i] & -plane_surfaces[i + 1]
  fuel_cells.append(openmc.Cell(fill=uo2, region=-fuel_or & layer, name='Fuel{:n}'.format(i)))
  gap_cells.append(openmc.Cell(fill=helium, region=+fuel_or & -clad_ir & layer, name='Gap{:n}'.format(i)))
  clad_cells.append(openmc.Cell(fill=zircaloy, region=+clad_ir & -clad_or & layer, name='Clad{:n}'.format(i)))
  water_cells.append(openmc.Cell(fill=borated_water, region=+clad_or & layer & -box, name='Water{:n}'.format(i)))
  all_cells.append(fuel_cells[i])
  all_cells.append(gap_cells[i])
  all_cells.append(clad_cells[i])
  all_cells.append(water_cells[i])

# Create a geometry and export to XML
geometry = openmc.Geometry(all_cells)
geometry.export_to_xml()

###############################################################################
# Define problem settings

# Indicate how many particles to run
settings = openmc.Settings()
settings.batches = 1500
settings.inactive = 500
settings.particles = 20000

# Create an initial uniform spatial source distribution over fissionable zones
lower_left = (-pitch/2, -pitch/2, 0.0)
upper_right = (pitch/2, pitch/2, height)
uniform_dist = openmc.stats.Box(lower_left, upper_right, only_fissionable=True)
settings.source = openmc.IndependentSource(space=uniform_dist)

if (args.entropy):
  entropy_mesh = openmc.RegularMesh()
  entropy_mesh.lower_left = lower_left
  entropy_mesh.upper_right = upper_right
  entropy_mesh.dimension = (10, 10, 20)
  settings.entropy_mesh = entropy_mesh

settings.temperature = {'default': 280.0 + 273.15,
                        'method': 'interpolation',
                        'range': (294.0, 3000.0),
                        'tolerance': 1000.0}

settings.export_to_xml()

# create some plots to look at the geometry for the sake of the tutorial
plot1          = openmc.Plot()
plot1.filename = 'plot1'
plot1.width    = (pitch, pitch)
plot1.basis    = 'xy'
plot1.origin   = (0.0, 0.0, height/2.0)
plot1.pixels   = (1000, 1000)
plot1.color_by = 'cell'

plot2          = openmc.Plot()
plot2.filename = 'plot2'
plot2.width    = (pitch, height)
plot2.basis    = 'xz'
plot2.origin   = (0.0, 0.0, height/2.0)
plot2.pixels   = (100, int(100 * (height/2.0/pitch)))
plot2.color_by = 'cell'

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

Figure 2: OpenMC geometry colored by cell ID shown on the - and - planes

The top and bottom of the pincell are vacuum boundaries. The four lateral faces of the pincell are reflective. Because we are not filling any universes with other universes or lattices, all of the cells in this problem are at the highest level of the geometry - i.e. the root universe.

We use a linear-linear stochastic interpolation between the two cross section data sets nearest to the imposed temperature by setting the method parameter on settings.temperature to interpolation. Finally, when OpenMC is initialized, cross section data is only loaded for a range of temperatures to save on memory usage; this range is specified with the temperature range element. You should set the range to be larger than the temperature range you expect in your coupled multiphysics problem - otherwise, you will encounter a runtime error that data is not available at the requested temperatures.

Because OpenMC runs at the end of each timestep, the initial fluid temperature is set to 280°C. And as 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 generate the XML files needed to run OpenMC, you can run the following:

python make_openmc_model.py

or simply use the XML files checked in to the tutorials/lwr_solid directory.

Solid Input Files

The solid phase is solved with the MOOSE heat transfer module, and is described in the solid.i input. The solid mesh is created using mesh generators in the mesh.i input:

r_pin = 0.39218
clad_ir = 0.40005
clad_or = 0.45720
L = 300.0

[Mesh]
  [clad] # This makes a circular annulus that will represent the clad
    type = AnnularMeshGenerator
    nr = 3
    nt = 20
    rmin = ${clad_ir}
    rmax = ${clad_or}
    quad_subdomain_id = 1
    tri_subdomain_id = 0
  []
  [extrude_clad] # this extrudes the circular annulus in the axial direction
    type = AdvancedExtruderGenerator
    input = clad
    heights = '${L}'
    num_layers = '40'
    direction = '0 0 1'
  []

  # A sideset in MOOSE is both an ID and a name.
  # This renames the sideset numbers on the clad inner surface (0) and outer surface (1)
  # to (5) and (4), respectively. This is done to avoid name collisions with another
  # AnnularMeshGenerator we use for the fuel pellet. After we rename the sideset IDs, we
  # need another RenameBoundaryGenerator to rename the sideset names.
  [rename_clad]
    type = RenameBoundaryGenerator
    input = extrude_clad
    old_boundary = '1 0' # outer surface, inner surface
    new_boundary = '5 4'
  []
  [rename_clad_names]
    type = RenameBoundaryGenerator
    input = rename_clad
    old_boundary = 'rmax rmin' # outer surface, inner surface
    new_boundary = 'rmax_c rmin_c'
  []

  [fuel] # this makes a circle that will represent the fuel
    type = AnnularMeshGenerator
    nr = 10
    nt = 20
    rmin = 0
    rmax = ${r_pin}
    quad_subdomain_id = 2
    tri_subdomain_id = 3
    growth_r = -1.2
  []
  [extrude] # this extrudes the circle in the axial direction
    type = AdvancedExtruderGenerator
    input = fuel
    heights = '${L}'
    num_layers = '40'
    direction = '0 0 1'
  []
  [combine]
    type = CombinerGenerator
    inputs = 'rename_clad_names extrude'
  []

  # one of the mesh generators does not work with distributed mesh
  parallel_type = replicated
[]

We generate the mesh by running cardinal-opt -i mesh.i --mesh-only to create the mesh_in.e file, which we then use in the solid input file.

L = 300.0
T_fluid = ${fparse 280.0 + 273.15}

[Mesh]
  [file]
    type = FileMeshGenerator
    file = mesh_in.e
  []
[]

The heat transfer module will solve for temperature, which we define as a nonlinear variable and apply a simple uniform initial condition.

[Variables]
  [temp]
    initial_condition = ${T_fluid}
  []
[]

The Transfer system in MOOSE is used to communicate variables across applications; a heat source will be computed by OpenMC and applied as a source term in MOOSE. In the opposite direction, MOOSE will compute a temperature that will be applied to the OpenMC geometry. Because we already solve for temperature, we simply need to add an auxiliary variable to receive the heat source from OpenMC.

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

The governing equation solved by MOOSE is specified in the [Kernels] block with the HeatConduction and CoupledForce kernels. The heat source provided by OpenMC is given by the receiver heat_source auxiliary variable.

[Kernels]
  [hc]
    type = HeatConduction
    variable = temp
  []
  [heat]
    type = CoupledForce
    variable = temp
    v = heat_source
  []
[]

We also set thermal conductivity values in the pellet and clad. Constant values are used for simplicity.

[Materials]
  [k_clad]
    type = GenericConstantMaterial
    prop_values = '0.5'
    prop_names = 'thermal_conductivity'
    block = '1'
  []
  [k_fuel]
    type = GenericConstantMaterial
    prop_values = '0.05'
    prop_names = 'thermal_conductivity'
    block = '2 3'
  []
[]

Next, the boundary conditions on the solid are applied, including the thermal contact model between the pellet and the clad.

[Functions]
  [T_fluid]
    type = ParsedFunction
    expression = '573.0 + 50.0 * (z / 300.0)'
  []
[]

[BCs]
  [surface]
    type = ConvectiveFluxFunction
    T_infinity = T_fluid

    # convert from W/m2/K to W/cm2/K
    coefficient = ${fparse 1000.0/100.0/100.0}
    variable = temp
    boundary = 'rmax_c'
  []
[]

[ThermalContact]
  # This adds boundary conditions bewteen the fuel and the cladding, which represents
  # the heat flux in both directions as
  # q''= h * (T_1 - T_2)
  # where h is a conductance that accounts for conduction through a material and
  # radiation between two infinite parallel plate gray bodies.
  [one_to_two]
    type = GapHeatTransfer
    variable = temp
    primary = 'rmax'
    secondary = 'rmin_c'

    # we will use a quadrature-based approach to find the gap width and cross-side temperature
    quadrature = true

    # emissivity of the fuel
    emissivity_primary = 0.8

    # emissivity of the clad
    emissivity_secondary = 0.8

    # thermal conductivity of the gap material
    gap_conductivity = 1.0

    # geometric terms related to the gap
    gap_geometry_type = CYLINDER
    cylinder_axis_point_1 = '0 0 0'
    cylinder_axis_point_2 = '0 0 ${L}'
  []
[]

In this example, the overall calculation workflow is as follows:

1. Run MOOSE heat conduction with a given power distribution from OpenMC.

2. Send temperature to OpenMC by modifying the temperature of OpenMC cells. A volume average is performed over all the MOOSE elements that mapped to a particular OpenMC cell.

3. Run OpenMC with an updated temperature distribution.

4. Extract the kappa-fission distribution (the recoverable fission energy) computed by OpenMC and map in the opposite direction from OpenMC cells to all the MOOSE elements that mapped to each cell.

The above sequence is repeated until desired convergence of the coupled domain is achieved. The MultiApps and Transfers blocks describe the interaction between Cardinal and MOOSE. The MOOSE heat conduction application is run as the main application, with OpenMC run as the sub-application. We specify that MOOSE will run first on each time step.

Two transfers are required to couple OpenMC and MOOSE for heat source and temperature feedback. The first is a transfer of heat source from Cardinal to MOOSE. The second is transfer of temperature from MOOSE to Cardinal.

[MultiApps]
  [openmc]
    type = TransientMultiApp
    input_files = 'openmc.i'
    execute_on = timestep_end
  []
[]

[Transfers]
  [heat_source_from_openmc]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    from_multi_app = openmc
    variable = heat_source
    source_variable = heat_source
    from_postprocessors_to_be_preserved = heat_source
    to_postprocessors_to_be_preserved = source_integral
  []
  [temp_to_openmc]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    to_multi_app = openmc
    variable = temp
    source_variable = temp
  []
[]

For the heat source transfer from OpenMC, we ensure conservation by requiring that the integral of heat source computed by OpenMC (in the heat_source postprocessor) matches the integral of the heat source received by MOOSE (in the source_integral postprocessor). We also add a postprocessor to evaluate the maximum solid temperature.

[Postprocessors]
  [source_integral]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
    execute_on = transfer
  []
  [max_T]
    type = NodalExtremeValue
    variable = temp
  []
[]

Because we did not specify sub-cycling in the [MultiApps] block, this means that OpenMC will run for exactly the same number of time steps (but the actual time step size used by the OpenMC wrapping is of no consequence because OpenMC is run in -eigenvalue mode). By setting a fixed number of time steps, this example will simply run a fixed number of Picard iterations.

[Executioner]
  type = Transient
  nl_abs_tol = 1e-8
  num_steps = 5
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

[Outputs]
  exodus = true
  print_linear_residuals = false
[]

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 on which OpenMC will receive temperature from the coupled MOOSE application, and on which OpenMC will write the fission heat source. In this example, we use exactly the same solid mesh as the coupled MOOSE application; this is not a requirement, and is relaxed in other tutorials.

[Mesh]
  [file]
    type = FileMeshGenerator
    file = mesh_in.e
  []
[]

Next, the Problem and Tallies blocks describe all objects necessary for the actual neutronics solve. To replace MOOSE finite element calculations with OpenMC particle transport calculations, the OpenMCCellAverageProblem class is used.

[Problem]
  type = OpenMCCellAverageProblem
  verbose = true
  power = '${fparse 3000e6 / 273 / (17 * 17)}'
  temperature_blocks = '1 2 3'
  cell_level = 0

  volume_calculation = vol

  [Tallies]
    [heat_source]
      type = CellTally
      blocks = '2 3'
      name = heat_source
    []
  []
[]

For this example, we first start by specifying that we wish to add a CellTally in [Tallies]. The blocks are then used to indicate which OpenMC cells to add tallies to (as inferred from the mapping of MOOSE elements to OpenMC cells). If not specified, we add tallies to all OpenMC cells. But for this problem, we already know that the cladding doesn't have any fissile material, so we can save some effort with the tallies by skipping tallies in those regions by setting blocks to blocks 2 and 3. OpenMCCellAverageProblem will then take the information provided in the [Tallies] block and add the necessary OpenMC tally.

For this example, we specify the total fission power by which to normalize OpenMC's tally results (because OpenMC's tally results are in units of eV/source particle). Next, we indicate which blocks in the [Mesh] should be considered for temperature feedback using the temperature_blocks parameter. Here, we specify temperature feedback for the pellet (blocks 2 and 3) and the cladding (block 1). During the initialization, OpenMCCellAverageProblem will automatically map from MOOSE elements to OpenMC cells, and store which MOOSE elements are providing feedback. Then when temperature is sent into OpenMC, that mapping is used to compute a volume-average temperature to apply to each OpenMC cell. We specify the level in the geometry on which the cells exist. Because we don't have any lattices or filled universes in our OpenMC model, the cell level is zero.

Finally, we add a stochastic volume calculation in order to compare the actual volumes of OpenMC cells against the [Mesh] elements to which they map. A good mapping should show close agreement between these two values.

[UserObjects]
  [vol]
    type = OpenMCVolumeCalculation
    n_samples = 100000
  []
[]

Next, we add a series of auxiliary variables for solution visualization (these are not requried for coupling). To help with understanding how Cardinal volume-averages temperature over the mesh, we add a CellTemperatureAux.

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

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

[Problem]
  type = OpenMCCellAverageProblem
  verbose = true
  power = ${fparse 3000e6 / 273 / (17 * 17)}
  temperature_blocks = '1 2 3'
  cell_level = 0

  volume_calculation = vol

Next, we specify an executioner and output settings. Even though OpenMC technically performs a criticality calculation (with no time dependence), we use the transient executioner so that if we wanted to run OpenMC more times than the coupled main application via subcycling, we would have a way to control that.

[Executioner]
  type = Transient
[]

[Outputs]
  exodus = true
  csv = true
[]

Finally, we add a postprocessor to evaluate the total heat source computed by OpenMC. We also include a TallyRelativeError postprocessor to evaluate the maximum relative error of the cell tally and a third postprocessor to evaluate the maximum heat source.

[Postprocessors]
  [heat_source]
    type = ElementIntegralVariablePostprocessor
    variable = heat_source
  []
  [max_tally_rel_err]
    type = TallyRelativeError
  []
  [max_heat_source]
    type = ElementExtremeValue
    variable = heat_source
  []
[]

You will likely notice that many of the always-included MOOSE blocks are not present in the openmc.i input. OpenMCCellAverageProblem automatically adds the heat_source, temp, and density (if density is coupled to OpenMC) variables in the openmc.i input, so these will never appear in the OpenMC wrapper file explicitly. It is as if the following is included in the input file:

[AuxVariables]
  [heat_source]
    family = MONOMIAL
    order = CONSTANT
  []
  [temp]
    family = MONOMIAL
    order = CONSTANT
  []
  [density] # only present if density_blocks was provided (not this example)
    family = MONOMIAL
    order = CONSTANT
  []
[]