Layered 1D LWR Rod Tutorial

Layered 1D models, sometimes referred to as 1.5D models, consist of several 1D cylindrical fuel 'slices' that are coupled together with out-of-plane strain variables. In this tutorial, we'll work with the example problem in the examples/tensor_mechanics/1.5D_rodlet_10pellets/ directory. Our tutorial will focus on walking through each section of the existing input file, and we encourage users to make incremental changes to the base example input file as a way to explore the capabilities of Bison.

This example problem uses 10 equally sized axial slices to model a 10 pellet rodlet: one axial slice per fuel pellet. As such, this simulation can be run in a matter of minutes on a desktop or laptop machine.

Layered 1D Model Definition

The term 'Layered 1D' in Bison is defined as a cylindrical fuel geometry model composed of coupled 1D simulations on each axial fuel slice at given axial positions along the rod. The fuel rod is divided into several axial slices, where 1D versions of the thermal and mechanical governing PDEs are solved.

note:Bison Layered 1D Model Overview

See the Introduction to Layered 1D Models for a description of how the Layered 1D capability is implemented in Bison, including a useful schematic which illustrates the key concept of axial slices.

Global Parameters

Global parameters offer users a quick way to define constant parameters that are used several places in an input file. Users must still exercise care to avoid overzealous use of these global parameters which may produce incorrect simulation results.

The Bison team recommends a standard set of global parameters for all LWR simulations, as shown in the input file snippet below:


[GlobalParams]
  density = 10431.0
  initial_porosity = 0.05
  order = SECOND
  family = LAGRANGE
  energy_per_fission = 3.2e-11 # J/fission
  displacements = disp_x
  temperature = temperature
[]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The density, initial_porosity, and energy_per_fission parameters are used to define the characteristics of the UO fuel in this LWR rodlet. The parameters order and family are used to define the shape functions used in the underlying finite element calculations, see the MOOSE documentation on finite element concepts for more information about shape functions. Finally, the displacements and temperature parameters are used to define the primary field variables in the thermo-mechanical problem solved by Bison.

Coordinate System

The Bison Layered 1D models are designed to work only in the cylindrical system which is ideal for fuel rods. We tell Bison to use the cylindrical coordinate system with the setting


[Problem]
  # Specify coordinate system type
  coord_type = RZ
[]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Layered 1D Mesh Generator

A mesh for a Layered 1D LWR fuel model consists of rows of 1D elements at each axial elevation to represent the fuel and cladding, and one more additional row of elements to represent the cladding in the plenum region, as described in Introduction to Layered 1D Models.

This mesh can be created externally to Bison, but to facilitate the setup of the Layered 1D models, the Bison team has developed an internal mesh generator: Layered1DMesh. To use this generator, we need to provide basic information about the fuel and cladding geometry, as shown in the example.


[Mesh]
  type = Layered1DMesh
  patch_update_strategy = auto
  partitioner = centroid
  centroid_partitioner_direction = y
  slices_per_block = 10
  clad_gap_width = 8.0e-5
  clad_thickness = 0.00056
  fuel_height = 0.1186
  plenum_height = 0.027
[]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Since we are modeling the rodlet with 10 equal slices, we set the parameter slices_per_block = 10 to divide the fuel into 10 sections. Generally we recommend using at least 10 slices to model a fuel rod. A user can specify the heights of each axial slice when equally spaced slices are not desired. See the Layered1DMesh documentation for more information.

Fuel Pin Geometry User Object

Along with the internal Bison mesh generator, we will also use the Layered1DFuelPinGeometry user object code class to pass information about the mesh to the function and the Layered1DMaster action, to calculate the axial pressure from the plenum and coolant, and to a variety of postprocessors.

Adding the fuel pin geometry user object is straight forward:


[UserObjects]
  [./pin_geometry]
    type = Layered1DFuelPinGeometry
  [../]
[]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Remember the name that we gave to the fuel pin geometry object, pin_geometry; we'll use this name as an argument to parameters in several other input file sections.

Thermal PDE Components

In this tutorial we'll focus first on adding in the components of the thermal PDE and the associated constitutive material properties before adding in the mechanical analysis components. For a quick refresher on the PDEs solved in a Bison simulation, refer to the Bison Fundamental Equations page.

Temperature Variable

We start by creating the field variable for temperature, as shown below. While it is possible to give the temperature variable any name, we use temperature here as the name of the variable block.

Note that we also set the value of the initial_condition in the temperature block. Here we use a higher value from the inlet coolant temperature instead of the usual room temperature value.

note:Initial Temperature Generally Set to Room Temperature

Generally in LWR analysis, the initial temperature condition is set to the fabrication temperature. Usually the fabrication temperature is assumed to be room temperature, that is, 293K.


[Variables]
  [./temperature]
    initial_condition = 580.0 # set initial temperature to coolant inlet
  [../]
[]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Since this is the first variable, we created the outer [Variable] block at the same time as creating the temperature variable. In fact, the temperature variable will be the only variable we explicitly define in this Layered1D input file: we'll use the Layered1D Master action to create and define the displacement variables.

Temperature PDE Terms (Kernels)

Recall the governing energy balance equation for Bison has three terms: (1) where , and are the temperature, density and specific heat, respectively, is the temperature, is the thermal conductivity, is the energy released in a single fission event, and is the volumetric fission rate. can be prescribed as a function of time and space, or input from a separate neutronics calculation. In the Layered1D input file, each term in Eq. 1 is represented with a single kernel, as shown below.

All three of the kernels we will discuss below should be added to a new input file section called [Kernels]. If you are not sure how to create this new [Kernels] section in the input file, click on the blue path link to the complete input file that is shown at the bottom of all the input file examples we use in this tutorial.

In all three of these kernels, we set variable = temperature because these kernels act on the temperature variable.

The first term in Eq. 1 is modeled with the HeatConductionTimeDerivative kernel:


[./heat_ie] # time term in heat conduction equation
  type = HeatConductionTimeDerivative
  variable = temperature
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The second term in Eq. 1, the gradient of the heat flux, is represented with the HeatConduction in a similar fashion:

[./heat] # gradient term in heat conduction equation
  type = HeatConduction
  variable = temperature
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The final term in Eq. 1 is only applicable for the fuel pellets, where heat from fission contributes as a heat source term. We model this fission heat source with the NeutronHeatSource kernel. Because this kernel relies on the fission rate, we specify the additional parameter burnup_function = burnup. The function for burnup is created by the Burnup Action system which we will discuss in the next section of this tutorial.


[./heat_source] # source term in heat conduction equation
  type = NeutronHeatSource
  variable = temperature
  block = fuel # fission rate applied to the fuel (block 2) only
  burnup_function = burnup
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Note that in NeutronHeatSource input file section we must restrict this kernel to act only on the fuel pellets mesh block; we accomplish this by setting the parameter block = fuel in the input file. This block setting was not necessary in the previous two kernel input file sections because Bison applies the kernels, and all input file sections, to all of the blocks of a mesh by default.

Thermal Materials

The energy balance equation, Eq. 1, requires two constitutive thermal material properties: the thermal conductivity, , and the density, . Because these thermal properties are material dependent, our Layered1D input file must include two sets of thermal material properties: one for the fuel pellets and and one for the cladding mesh block.

All four of the materials we will discuss below should be added to a new input file section called [Materials]. If you are not sure how to create this new [Materials] section in the input file, click on the blue path link to the complete input file that is shown at the bottom of all the input file examples we use in this tutorial.

Fuel Pellet Thermal Materials

To restrict these thermal material properties to the fuel pellet mesh block, we set the parameter block = fuel in each of the two fuel thermal material properties input file sections.

To model the thermal conductivity of the fuel pellets, we use the ThermalFuel material. As discussed on the documentation page, the ThermalFuel has several models for the thermal conductivity of UO; we will use the NFIR model. To select the NFIR model, we set the parameter model = 4 because NFIR is the fourth model. Additionally, we need to provide the ThermalFuel material with the names of the coupled temperature variable (temp = temperature) and the burnup function name (burnup_function = burnup).

[./fuel_thermal]
  type = ThermalFuel
  block = fuel
  model = 4 # NFIR
  temp = temperature
  burnup_function = burnup
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The density of the fuel is defined with the Density material:


[./fuel_density]
  type = Density
  block = fuel
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Because we already defined the fuel density as a global parameter (in the [GlobalParams] section), at the start of this tutorial, we do not need to respecify the value of density for the fuel pellet mesh block.

Cladding Thermal Materials

To restrict these thermal material properties to the cladding mesh block, we set the parameter block = clad in each of the two cladding thermal material properties input file sections.

Since the cladding thermal conductivity is not a function of the fuel burnup, we use the general HeatConductionMaterial to calculate the thermal conductivity of the zirconium cladding. We assume the thermal conductivity and the specific heat are constant by specifying constant values for both of these parameters as shown below:


[./clad_thermal]
  type = HeatConductionMaterial
  block = clad
  thermal_conductivity = 16.0
  specific_heat = 330.0
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

As in the fuel pellet section, we use the same Density material to provide information about the density of the cladding. In this case, we must specify the value of the density for zirconium to overwrite the setting in the [GlobalParams] section of the input file:


[./clad_density]
  type = Density
  block = clad
  density = 6551.0
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Burnup Calculation

In Bison we use the Burnup Action system to determine the values of the fuel burnup and the fission rate across the radial and axial profile of the fuel rod. Because the burnup model is included in the input file through the use of a type of code called an Action, this input file section does not have a type parameter. Instead, the specific input file section name [Burnup] tells Bison to use the Burnup Action. Actions are types of code elements which greatly simplify and shorten the input file by combining common parameter options.

As in other fuel specific sections of this input file, we set the parameter block = fuel to ensure that the Burnup Action is active only on the fuel mesh blocks.

The Burnup Action system creates a secondary burnup-specific grid to over the fuel rod to interpolate the burnup and fission rate at many points across the fuel pellets. This secondary grid is a set of rectangles which we define with the parameters num_radial and num_axial. Since the burnup and fission rates change quickly throughout the fuel pellet radial thickness, we use the relatively high num_radial = 80 to set 80 divisions of the burnup-specific interpolation grid across the radial dimension of the fuel rod. We use the setting num_axial = 11 because this rodlet is short (the default value is 20 axial divisions). The parameters order and family tell Bison how to perform the interpolation for burnup on this secondary grid, as discussed in the MOOSE documentation on finite element shape functions.


[Burnup]
  [./burnup]
    block = fuel
    rod_ave_lin_pow = power_history # using the power function defined above
    axial_power_profile = axial_peaking_factors # using the axial power profile function defined above
    num_radial = 80
    num_axial = 11
    order = CONSTANT
    family = MONOMIAL
    fuel_pin_geometry = pin_geometry
    fuel_volume_ratio = 1.0 # for use with dished pellets (ratio of actual volume to cylinder volume)
    RPF = RPF
    i_enrich = '0.05 0.95 0 0 0 0'
  [../]
[]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Because of this secondary mesh grid, we need to provide the Burnup Action with specific information about the geometry of the fuel rod. We provide this information with the name of the fuel pin geometry user object we already defined: fuel_pin_geometry = pin_geometry. The parameter setting fuel_volume_ratio = 1 simply tells Bison that our fuel pellets are not dished.

The Burnup Action can also retrieve information about burnup, fission rate, and heavy metal isotope concentrations from the burnup-specific interpolation grid; this information is then saved to a set of auxkernels created by the action. In this case we've elected to have Bison save information about the radial power factor to an auxkernel named RPF with the setting RPF = RPF.

Finally we specify the initial enrichment of the fuel, i_enrich as 5% U-235, 95% U-238 with no Pu-239, Pu-240, Pu-241, nor Pu-242 content. This initial enrichment is the default setting for the Burnup Action system.

Functions Associated with Burnup

As discussed in the Burnup section of this tutorial, the Burnup Action requires information about the power time history and the axial power distribution. We will start a new input file section called [Functions] with these three functions discussed below.

note:Inputs in SI Units Required

The csv files used to provide information about the power history and power axial profile should use SI units: seconds for time, W/m for power, and m for axial position.

We use the PiecewiseLinear function type to read in a csv file containing information about the power over time. This file includes the time (s) and the average linear power (W/m). The csv data file, here named powerhistory.csv must be located in the same directory as the input file. The default value of unity is used for the scale_factor parameter in this example.

[./power_history]
  type = PiecewiseLinear # reads and interpolates an input file containing rod average linear power vs time
  data_file = powerhistory.csv
  scale_factor = 1
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

To provide Bison with information about the distribution of the power across the axial height of the fuel rodlet, we use the PiecewiseBilinear function. This file includes information about the power (W/m) at given axial heights(m), relative to the bottom of the fuel pellet stack, as a function of time(s). Generally these values are provided by the test reactors. The csv data file used in this input file section is named peakingfactors.csv, and this csv file must also be located in the same directory as the input file.

As in the PiecewiseLinear function, a scale_factor of unity is applied. The axis parameter tells Bison that the locations provided in the csv data file are along the y axis.


[./axial_peaking_factors] # reads and interpolates an input file containing the axial power profile vs time
  type = PiecewiseBilinear
  data_file = peakingfactors.csv
  scale_factor = 1
  axis = 1 # (0,1,2) => (x,y,z)
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

These two files can be combined into a single piecewise function using the CompositeFunction to create an axial profile of the power over time. This composite function is used to calculate the fuel relocation, which is discussed under the Mechanical PDE Components section of this tutorial.


[./q]
  type = CompositeFunction
  functions = 'power_history axial_peaking_factors'
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Additional functions are required to correctly generate and apply the boundary conditions associated with the pressure on the fuel rod. These additional functions are discussed under the Functions Associated with Pressure part of this tutorial.

Mechanics PDE Components

Now that we have defined the thermal PDE components, we turn our attention to the mechanical governing PDE, that is, the stress divergence equation: (2) where is the stress, , , and are the elastic, inelastic, and eigen strains, respectively, and is the elasticity tensor.

The TensorMechanics MasterAction system is used to create the displacement nonlinear variables, the stress divergence kernel, and the strain calculator material. The Layered1D Master action sytem, designed specifically for the Layered 1D capability in Bison, adds the scalar out-of-plane variables, out-of-plane strain, and axial pressure required by the Layered 1D model.

As with other actions, the Layered1D Master action simplifies the input file by gathering all of these different input file sections into a single master input file section. The specific input file section name [Modules/TensorMechanics/Layered1DMaster] indicates the use of the Layered 1D master action. Since the mechanical behavior of the fuel and clad are different, we will define two different Layered1DMaster action input file sections. These two input file sections, discussed below, discussed below are both included within the [Modules/TensorMechanics/Layered1DMaster] section.

Fuel Pellet MasterAction

We restrict this version of the Layered1D Master to the fuel mesh block with the parameter setting block = fuel, just as we've done in previous input file sections.

With the next parameter, add_variables = true, we instruct the action to create the nonlinear displacement variables required for this mechanical problem. The Layered1DMaster action refers to the global parameter, in the [GlobalParams] section of the input file, displacements = 'disp_x' to determine that only one nonlinear displacement variable with the name disp_x should be created.

Next we tell the Layered1DMaster action to create the ComputeAxisymmetric1DFiniteStrain strain calculator material with the parameter setting strain = FINITE. The Layered1DMaster action detects that this particular strain calculator material is appropriate based on the strain parameter setting, the mesh geometry, and the number of displacement variables specified with the displacements parameter. The Layered1DMaster action also sets the StressDivergenceRZTensors to solve the governing stress divergence equation, Eq. 2.

We use the Layered1DMaster action to create the scalar variables for the out-of-plane strains for all 10 axial fuel 1D layers with the parameters add_scalar_variables and out_of_plane_strain_name.

Additionally, to properly account for the axial, or out-of-plane, strain due to axial pressure on the fuel from the plenum, we set the out_of_plane_pressure parameter to the name of the fuel axial pressure function we'll discuss in the Boundary Condition part of this tutorial.

As in the Burnup section of the input file, we use the name of the fuel pin geometry user object we already defined to provide the Layered1DMaster action with specific geometry of the fuel by setting fuel_pin_geometry = pin_geometry.

[./fuel]
  block = fuel
  add_variables = true
  strain = FINITE
  add_scalar_variables = true
  out_of_plane_strain_name = strain_yy
  fuel_pin_geometry = pin_geometry
  out_of_plane_pressure = fuel_axial_pressure
  eigenstrain_names = 'fuelthermal_strain swelling_strain fuel_relocation_strain'
  generate_output = 'stress_xx stress_yy stress_zz vonmises_stress strain_xx'
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Next we provide the Layered1DMaster action with a list of all the eigenstrains acting on the fuel block as the argument to the eigenstrain_names parameter. Eigenstrains are strains which do not result directly from an applied force but do change the volume of the material. We'll define the three different eigenstrains acting on the fuel mesh block, the thermal expansion, volumetric swelling, and relocation eigenstrains in the Mechanical Materials part of this tutorial.

Finally, we use the Layered1DMaster action to generate the set of auxkernels and auxvariables necessary to output various scalar quantities derived from the Rank-2 tensors created: the stress and strain tensors. Here we have selected to output the diagonal components of the stress tensor, the von Mises stress measure, and the radial strain by listing these quantities as arguments for the generate_output parameter.

Cladding MasterAction

The settings for the Layered1D Master for the cladding mesh block follow the same pattern as discussed for the fuel mesh block with only a few exceptions. As we'd expect, we will restrict this version of the Layered1DMaster action to operate solely on the clad mesh block with the parameter setting block = clad.

The cladding is subjected to out-of-plane, or axial, pressures from both the plenum, on the inner clad surface, and the coolant, which acts on the outer clad surface. To reflect this different out-of-plane pressure, we set the out_of_plane_pressure parameter to a different function name. As mentioned before, we'll discuss this clad axial pressure function in the Boundary Conditions part of this tutorial.


[./clad]
  block = clad
  add_variables = true
  strain = FINITE
  add_scalar_variables = true
  out_of_plane_strain_name = strain_yy
  fuel_pin_geometry = pin_geometry
  out_of_plane_pressure = clad_axial_pressure
  eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_strain'
  generate_output = 'stress_xx stress_yy stress_zz vonmises_stress strain_xx'
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The final difference in the cladding Layered1DMaster action settings is for the eigenstrain_names parameter: since eigenstrains are material dependent, the cladding has a different set of eigenstrains. We'll create the input file sections for the clad thermal expansion and irradiation swelling eigenstrains in the next part of this tutorial under Mechanical Materials.

Mechanical Materials

Since the Layered1D Master action sections have already defined the strain calculator materials, we only need to focus on providing Bison with information about the elasticity tensor, stress, and eigenstrain materials in the [Materials] section of the input file.

Fuel Pellet Mechanical Materials

We begin by defining the elasticity tensor for the fuel block with the general ComputeIsotropicElasticityTensor. We set the constant elastic constants youngs_modulus and poissons_ratio to the appropriate values shown in our example. In addition we set the block = fuel parameter to restrict this elasticity tensor to the fuel mesh block.


[./fuel_elasticity_tensor]
  type = ComputeIsotropicElasticityTensor
  block = fuel
  youngs_modulus = 2.0e11
  poissons_ratio = 0.345
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

In this example problem, we will model the fuel behavior as simply elastic with the ComputeFiniteStrainElasticStress. In this input file section we only need to add the mesh block restriction with block = fuel.


[./fuel_stress]
  type = ComputeFiniteStrainElasticStress
  block = fuel
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

As described in the the Layered1DMaster action section above, the eigenstrain materials represent changes in the strain that are not directly caused by applied force. In the case of the fuel pellets, we include three separate eigenstrains for thermal expansion, volumetric swelling due to fission products, and a relocation strain which accounts for fuel cracking.

In all of these eigenstrain materials, it is imperative that the name given as the argument to the eigenstrain_name parameter matches the name given in the string for the Layered1DMaster action eigenstrain_names parameter. These eigenstrain names are used by Bison to connect the eigenstrain materials to the Layered1DMaster action.

The first fuel eigenstrain we will add is due to the thermal expansion of the fuel, and we model with eigenstrain with the ComputeThermalExpansionEigenstrain. We restrict this mateiral to the fuel pellets mesh block with the block = fuel parameter.

The coefficient of thermal expansion is material dependent and is set with the thermal_expansion_coeff parameter. The stress_free_temperature is the temperature (in K) at which the material does not change volume due to temperature. Generally, if this value is nor provided by experimental data, we set the stress_free_temperature to the fabrication or room temperature, as we've done in this example. As discussed above, we use the setting eigenstrain_name = fuel_thermal_strain to match the name given in the Layered1DMaster action.


[./fuel_thermal_strain]
  type = ComputeThermalExpansionEigenstrain
  block = fuel
  thermal_expansion_coeff = 10.0e-6
  stress_free_temperature = 295.0
  eigenstrain_name = fuelthermal_strain
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The second fuel eigenstrain, the volumetric swelling of the fuel pellets due to fission products, is modeled with UO2VolumetricSwellingEigenstrain. As before, we restrict this materials to act on only the fuel mesh block with block = fuel.

We use the Bison recommended model for gaseous fission product release, Sifgrs, with the setting gas_swelling_model_type = SIFGRS. We will discuss this fission gas model in the following Fission Gas Production and Release part of this tutorial. Since the fission products depend on burnup, we provide the UO2VolumetricSwellingEigenstrain material with the name of the burnup function as the argument to the burnup_function parameter. Similarly we set the parameter initial_fuel_density to the ceramic fuel density value.

In order to post process more detailed information about the different components of the volumetric swelling, we set both the save_densification and the save_solid_swelling parameters to true. Later in this tutorial we will add these quantities to the simulation data output with sets of AuxKernels, AuxVariables, and PostProcessors. Finally, we set the eigenstrain_name = swelling_strain to match the name given in the Layered1DMaster action.


[./fuel_swelling]
  type = UO2VolumetricSwellingEigenstrain
  block = fuel
  gas_swelling_model_type = SIFGRS
  burnup_function = burnup
  initial_fuel_density = 10431.0
  save_densification = true
  save_solid_swelling = true
  eigenstrain_name = swelling_strain
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The third and final fuel eigenstrain in our example simulation is the relocation eigenstrain, which we model with UO2RelocationEigenstrain. Relocation is an engineering scale approximation of the volume change in fuel pellets due to cracking.

As in previous input file sections, we set restrict this material with the block = fuel parameter, and provide the burnup and geometry information with the burnup_function and fuel_pin_geometry parameters. We will also use the name of the CompositeFunction for the axial profile of power applied to the rod over time, which we defined earlier. Recall from the Functions Associated with Burnup part of this tutorial that we gave this composite function the name q. We tell Bison to use this composite function in the UO2RelocationEigenstrain material by setting the parameter linear_heat_rate_function = q.

We use the recommended default parameters setting relocation_activation1 = 5000.0 for the first linear power bound in the calculation of the relocation strain coefficient. The burnup_relocation_stop is the burnup value (in FIMA) at which relocation is considered to stop, when the clad comes into contact with the fuel; this parameter setting is specific to each individual fuel analysis. For this example we have selected the value 0.24. Finally we select the relocation_model = ESCORE_modified. As with the other eigenstrains, we set the parameter eigenstrain_name = fuel_relocation_strain, taking care to ensure the name used here matches the Layered1DMaster action.


[./fuel_relocation]
  type = UO2RelocationEigenstrain
  block = fuel
  burnup_function = burnup
  fuel_pin_geometry = pin_geometry
  linear_heat_rate_function = q
  relocation_activation1 = 5000.0
  burnup_relocation_stop = 0.024
  relocation_model = ESCORE_modified
  eigenstrain_name = fuel_relocation_strain
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Cladding Mechanical Materials

As we did for the fuel materials, we begin the section of cladding mechanical materials by creating the elasticity tensor section with the zirconium specific ZryElasticityTensor material. This material model contains the defined elastic constants for zirconium. We also restrict this material to act only on the clad mesh block by setting block = clad.


[./clad_elasticity_tensor]
  type = ZryElasticityTensor
  block = clad
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Unlike in the fuel, we elect to use an inelastic material model, creep, to simulate the stress response of the zirconium cladding. This choice requires the use of two stress materials in the input file. In both stress input file sections we will restrict these materials to the clad mesh block with the setting block = clad.

The first stress material, ComputeMultipleInelasticStress is used to manage the nonlinear stress response iterations. We instruct Bison to use the more computationally efficient (for well behaved problems) tangent_operator = elastic. For more information on the tangent operator, see the ComputeMultipleInelasticStress documentation page. As we did with the eigenstrains and the Layered1DMaster action, we tell the ComputeMultipleInelasticStress material the names of the inelastic stress materials to call with the inelastic_models = zrycreep.


[./clad_stress]
  type = ComputeMultipleInelasticStress
  block = clad
  tangent_operator = elastic
  inelastic_models = 'zrycreep'
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

In the next input file section, we create the zrycreep material. It is important that the name of this block matches the argument given for the inelastic_models parameter in the ComputeMultipleInelasticStress input file section.

To model the zirconium creep we use the ZryCreepLimbackHoppeUpdate material model. Because the creep behavior depends on the irradiation environment, we pass the names of the neutron flux and neutron fluence auxvariables, which we will discuss in the AuxVariables & AuxKernels section of this tutorial, as the arguments to the fast_neutron_flux and fast_neutron_fluence parameters.

We elect to include contributions from irradiation, and primary and secondary thermal creep by setting all of the corresponding parameters to true. Note that this is the default behavior to include all three of these creep contributions. Finally we set the zircaloy_material_type = stress_relief_annealed to use the material properties for Zr-2 which has been annealed.


[./zrycreep]
  type = ZryCreepLimbackHoppeUpdate
  block = clad
  fast_neutron_flux = fast_neutron_flux
  fast_neutron_fluence = fast_neutron_fluence
  model_irradiation_creep = true
  model_primary_creep = true
  model_thermal_creep = true
  zircaloy_material_type = stress_relief_annealed
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

In the clad we will consider the effects of two eigenstrains. As in the fuel section above, it is critical to ensure that the names given to the individual eigenstrains with the eigenstrain_name parameter matches the names used in the Layered1DMaster action eigenstrain_names parameter.

The first cladding eigenstrain we add is also a thermal expansion eigenstrain, but in the case of the cladding we will use the material specific ZryThermalExpansionMatproEigenstrain. We restrict this thermal expansion eigenstrain to just the clad with the block = clad parameter setting. We use the same fabrication temperature value for the stress_free_temperature parameter as we did in the fuel thermal expansion section. We set eigenstrain_name = clad_thermal_eigenstrain.


[./clad_thermal_expansion]
  type = ZryThermalExpansionMatproEigenstrain
  block = clad
  stress_free_temperature = 295.0
  eigenstrain_name = clad_thermal_eigenstrain
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The second cladding eigenstrain in our example simulation is due to irradiation induced swelling, and we use the ZryIrradiationGrowthEigenstrain material. After restricting this material to the clad mesh block with block = clad, we pass the name of the fast neutron fluence auxvariable that we used in the creep stress input file section above to the fast_neutron_fluence parameter. We also use the same zircaloy_material_model_type = stress_relief_annealed as in the creep stress input file section to ensure consistency in the material coefficients. Finally we set the parameter eigenstrain_name = clad_irradiation_strain, ensuring that the name used in the eigenstrain matches the name used in the Layered1DMaster action.


[./clad_irradiation_swelling]
  type = ZryIrradiationGrowthEigenstrain
  block = clad
  fast_neutron_fluence = fast_neutron_fluence
  zircaloy_material_model_type = stress_relief_annealed
  eigenstrain_name = clad_irradiation_strain
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Fission Gas Production and Release

The last material model we will use in our Layered 1D example problem is the material model for the generation and release of fission gas from the fuel into the plenum space. As a material model, it will be added to the [Materials] section of the input file. We use the model Sifgrs to calculate the amount of fission gas generated as a result of fission in the fuel and the amount of fission gas released to the plenum.

Since only the fuel pellets undergo fission, we restrict the Sifgrs material to the fuel mesh block with the familiar block = fuel setting. We provide the name of the nonlinear temperature variable as the argument to the temperature parameter and the name of the burnup function to the corresponding burnup_function parameter, as shown below.


[./fission_gas_release]
  type = Sifgrs
  block = fuel
  temp = temperature
  burnup_function = burnup
  gbs_model = true
  grain_radius = grain_radius
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

We elect to couple the fission gas behavior to the growth of the ceramic fuel crystalline grains by setting the parameter gbs_model = true. In order to properly use the grain growth model, we must provide Sifgrs with the name of the auxvariable storing the grain radius size: we set the parameter grain_radius to this auxvariable name. We will discuss the creation of the grain radius auxvariable in the AuxVariables & AuxKernels section of this tutorial.

Displacement AuxVariables for Visualization

In order to view the simulation results with Paraview, the visualization tool recommended for use with Bison, we need to add two additional auxvariables to create a complete set of three displacement quantities for Paraview. Since the radial displacement, disp_x, is defined as a nonlinear variable by the Layered1D Master action, we add two auxvariables for the disp_y and disp_z displacements in the axial and azmuthial directions, respectively. We will put these two auxvariables in the [AuxVariables] section of the input file.

[./disp_y] ## Required for easier visualization in Paraview
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

[./disp_z] ## Required for easier visualization in Paraview
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Boundary Conditions

The boundary conditions, often abbreviated as BCs, help satisfy the essential and natural conditions required to solve a finite element problem. In our Layered1D model we apply boundary conditions on both the displacement variables (with a displacement or with a stress traction (e.g. pressure)) and on the temperature (with the CoolantChannel action). We'll discuss first the mechanical displacements boundary conditions before the thermal boundary conditions; the next three input file sections will go into the [BCs] section of the input file.

Displacements

Since we have only one nonlinear displacement variable, the radial displacement disp_x, in our Layered1D simulation, we only need one displacement boundary condition. We select the PresetBC to set the displacements along the given boundary to a constant value for the duration of the simulation.

In our example problem we use the setting boundary = 12 to apply this boundary condition to the centerline of the fuel pin. Since we are assuming an axisymmetric model, this centerline is shown as the left most boundary of our Layered1D model. For more information on the boundary settings used in Bison, refer to the Mesh Script documentation page.

To prevent the radial displacement, disp_x, from moving from the centerline, we set value = 0.0 to prescribe a value of zero displacement.

[./no_x_all] # pin pellets and clad along axis of symmetry (y)
  type = PresetBC
  variable = disp_x
  boundary = 12
  value = 0.0
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Pressures

We use two sets of pressure boundary conditions: one from the coolant which applies pressure to the outside of the cladding, and a second which applies pressure from the plenum pressure on both the fuel exterior and the cladding interior surfaces. Since the Layered 1D model is a series of individual 1D axial slices models, these pressures are only applied in the radial directions.

The first pressure we will add to our Layered1D example problem is from the coolant exerting pressure on the cladding exterior. We use the Pressure action system to create this pressure traction force. To apply this boundary condition to the cladding exterior, we use the setting boundary = 2. For more information on the boundary settings used in Bison, refer to the Mesh Script documentation page.

The evolution of the pressure on the cladding from the coolant is defined with the function named pressure_ramp; we will discuss the creation of this function in the next part of this tutorial, Functions Associated with Pressure. We use the pressure ramp function name as the argument to the function parameter here. The last parameter factor we use as a scaling factor for the pressure ramp function.


[./Pressure] # apply coolant pressure on clad outer walls
  [./coolantPressure]
    boundary = 2
    function = pressure_ramp # use the pressure_ramp function defined above
    factor = 15.5e6
  [../]
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The second pressure input file section we will add to our example Layered1D input file models the pressure from the plenum that acts on both the cladding interior and the exterior of the fuel pellet surfaces. We use the PlenumPressure action system. As shown in the example below, we use the parameter setting boundary = 9 to apply this pressure from the plenum on the combination of the pellet exterior and the cladding interior surfaces. For more information on the boundary settings used in Bison, refer to the Mesh Script documentation page.

The initial_pressure parameter is set to the value of the pressure in the plenum at the start time of the reactor operation. Generally the plenum is filled with an inert gas such that the inital_pressure is nonzero, but the value differs for each fuel rod. Following the general procedure for Bison fuel rod analyses, we set starup_time = 0 to indicate that the reactor begins operation at a time of 0 seconds. The parameter R is the universal gas constant, and in Bison we recommend the use of SI units such that R = 8.314 is the most common setting.


[./PlenumPressure] # apply plenum pressure on clad inner walls and pellet surfaces
  [./plenumPressure]
    boundary = 9
    initial_pressure = 2.0e6
    startup_time = 0
    R = 8.314
    output_initial_moles = initial_moles # coupling to post processor to get initial fill gas mass
    temperature = ave_temp_interior # coupling to post processor to get gas temperature approximation
    volume = gas_volume # coupling to post processor to get gas volume
    material_input = fis_gas_released # coupling to post processor to get fission gas added
    output = plenum_pressure # coupling to post processor to output plenum/gap pressure
  [../]
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Additionally the PlenumPressure action couples to many additional PostProcessors, which we will discuss in the PostProcessors part of this tutorial, to include information which impact the plenum. In these cases, the name of the postprocessor is the argument to the PlenumPressure action parameters. The output_initial_moles parameter is used to store the value of the initial moles of gas. This option is highly recommended as a debugging tool when comparing simulations of different geometries. The temperature parameter does not take the name of the nonlinear temperature variable but instead is intended for the name of the post processor that calculates the average temperature of the entire plenum space. The volume parameter argument is the name of the postprocessor that determines the current volume of the plenum cavity. The material_input parameter takes the name of the postprocessor that stores the amount of gaseous fission products released into the plenum. The last parameter, output, is given the name of the postprocessor which stores the current pressure in the plenum. This parameter setting is required for the proper functioning of Layered 1D models. All of these postprocessors will be discussed under the Postprocessors part of this tutorial.

Functions Associated with Pressure

We will add these additional pressure functions to the [Functions] section of the input file which we started when adding our burnup functions.

The first function we will add to our input file is the function which defines the change in pressure from the simulation start time to the time when the reactor is considered to start operation. We name this function pressure_ramp to match the name the PlenumPressure action is expecting. We use the PiecewiseLinear function here. We use the x parameter to specify the times and the y parameter to specify the different values of the pressure at those specific times. Following the convention, we set the simulation start time to -200 with zero pressure, and then at the reactor operation start time, 0 we consider the plenum pressure to be at the full capacity. Even though the simulation runs well past the 0 seconds specified by the parameters, the PiecewiseLinear function will continue to hold the last given value until the simulation ends.


[./pressure_ramp] # reads and interpolates input data defining amplitude curve for fill gas pressure
  type = PiecewiseLinear
  x = '-200 0'
  y = '0 1'
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The fuel_axial_pressure function is a Layered 1D specific requirement, and is required by the Layered1D Master action to determine the amount of axial pressure acting on the fuel. We use the particular name fuel_axial_pressure for this function to match the inputs we gave earlier for the Layered1DMaster action. We use the ParsedFunction and set all three parameters equal to the plenum_pressure postprocessor that is created by the PlenumPressure action.


[./fuel_axial_pressure]
  type = ParsedFunction
  value = plenum_pressure
  vars = plenum_pressure
  vals = plenum_pressure
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The clad_axial_pressure function is a Layered 1D specific requirement, and is used by the Layered1D Master action to determine the amount of axial pressure acting on the cladding from both the plenum and the coolant. We use the specific name clad_axial_pressure for this function to match the inputs we gave earlier for the Layered1DMaster action.

We use the Layered 1D specific CladdingAxialPressureFunction to create this function of combined pressures. The parameter plenum_pressure takes the name of the plenum_pressure postprocessor that is created by the PlenumPressure action. We use the name of the pressure_ramp function we just defined as the argument for the coolant_pressure parameter; the value of the coolant_pressure_scaling_factor should match the value of the factor parameter from the coolantPressure boundary condition we created with the Pressure action.


[./clad_axial_pressure]
  type = CladdingAxialPressureFunction
  plenum_pressure = plenum_pressure
  coolant_pressure = pressure_ramp
  coolant_pressure_scaling_factor = 15.5e6
  fuel_pin_geometry = pin_geometry
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Finally we provide the CladdingAxialPressureFunction with specific geometry information with the nearly ubiquitous Layered1DFuelPinGeometry with the parameter setting fuel_pin_geometry = pin_geometry.

CoolantChannel Model

We apply boundary conditions on the temperature nonlinear variable with the CoolantChannel action system. We set the parameter variable = temperature to tell Bison to apply the temperature profiles calculated with the CoolantChannel to the temperature variable. We set boundary = 2 to set temperature boundary condition on the outer cladding surface. For more information on the boundary settings used in Bison, refer to the Mesh Script documentation page.

The values for the inlet_temperature, inlet_pressure are specific to the individual rod rod; however, we generally set the inlet_temperature = 580. The value of the inlet_pressure parameter argument should correspond to the value used in creating pressure boundary conditions.

note:Match Inlet Pressure and Coolant Pressure Scaling Factors

Note that the value of the inlet_pressure matches the value of the scaling factor we used in defining the coolant pressure boundary conditions with the Pressure action and the CladdingAxialPressureFunction.


[CoolantChannel]
  [./convective_clad_surface] # apply convective boundary to clad outer surface
    variable = temperature
    boundary = 2
    inlet_temperature = 580 # K
    inlet_pressure = 15.5e6 # Pa
    inlet_massflux = 3800 # kg/m^2-sec
    rod_diameter = 0.948e-2 # m
    rod_pitch = 1.26e-2 # m
    linear_heat_rate = power_history
    axial_power_profile = axial_peaking_factors
  [../]
[]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The values of the rod_diameter and rod_pitch are specific to each fuel rod being analyzed and will change from Bison simulation to simulation. We use the settings shown above as an example in our Layered 1D tutorial input file.

The parameters which provide information about the power applied to the rod, linear_heat_rate and axial_power_profile take the same function names for arguments as we gave in the Burnup action input file section.

Thermal Contact

Within a fuel rod, heat transfers from the fuel pellets to the cladding first across the plenum, when there is a gap between the fuel and the cladding, and later directly from the fuel to the cladding once the two materials are in contact with each other. We model this transfer of heat with the GapHeatTransferLWR boundary condition. Bison knows to apply this contact boundary condition to the nonlinear temperature variable because we set the parameter variable = temperature.

The boundary condition differs from the conditions we have previously discussed because it borrows the concept of contact surfaces from contact theory. Traditionally the surface which acts on the second surface is termed the master surface and while the second surface in the contact pair is termed the slave surface. Using the boundaries defined by Bison's Mesh Script we set master = 5 for the interior surface of the cladding and slave = 10 for the exterior surface of all the fuel pellets.


[ThermalContact]
  [./thermal_contact]
    type = GapHeatTransferLWR
    variable = temperature
    master = 5
    slave = 10
    initial_moles = initial_moles # coupling to a postprocessor which supplies the initial plenum/gap gas mass
    gas_released = fis_gas_released # coupling to a postprocessor which supplies the fission gas addition
    contact_pressure = contact_pressure
  [../]
[]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The GapHeatTransferLWR boundary condition is coupled to thePlenumPressure action input file section through the number of moles of gas in the initial plenum; the parameter initial_moles takes the name of the postprocessor specified in the PlenumPressure. Similarly, the GapHeatTransferLWR boundary condition is connected to the Sifgrs fission gas release model with the parameter gas_released that takes the name of the postprocessor storing the current amount of fission gas released into the plenum. Finally the optional contact_pressure parameter is used to create the name of a postprocessor to store the value of the contact pressure for use by postprocessors.

Mechanical Contact

To enforce mechanical contact between the fuel pellet and the cladding, we use the Contact Action system. As with other actions previously discussed in this tutorial, the action simplifies the input file section and does not have a type parameter. The specific input file section name [Contact] indicates that Bison will use the Contact Action.

As in the thermal contact input file section, the recommended settings for the master and slave parameters are given below. The clad block is generally set as the master surface which applies pressure onto the slave fuel pellet surface.


[Contact]
  [./pellet_clad_mechanical]
    master = 5
    slave = 10
    formulation = kinematic
    model = frictionless
    system = Constraint
    penalty = 1e7
  [../]
[]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Since the axial slices of fuel and clad in the mesh block are not allowed to move vertically, the system = Constraint with the formulation = kinematic and model = frictionless are recommended in Layered 1D simulations. The penalty parameter is required by the kinematic formulation, and the setting recommended for Layered1D simulations is shown above.

AuxVariables & AuxKernels

Because AuxKernels and AuxVariables rely on each other strongly (the auxvariable stores the data calculated by the auxkernel), we will discuss and add these two different input file sections in pairs. Remember that the AuxKernels should be added into the [AuxKernels] section of the input file while the AuxVariables should be included in the [AuxVariables] input file section.

Fast Neutron Flux and Fluence

Within the cladding, the flow of fast neutrons from the pellets into the clad affects the mechanical material behavior. The fast neutron flux and fluence auxvariables are used by the cladding creep ZryCreepLimbackHoppeUpdate and eigenstrain ZryIrradiationGrowthEigenstrain.

To model the fast neutron flux, we use the FastNeutronFluxAux auxkernel. We provide the auxkernel with the name of the corresponding auxvariable and restrict this auxkernel to act on only the clad mesh block with block = clad. As in the Burnup action, the parameters rod_ave_lin_pow and axial_power_profile are set to the names of the power history and peaking factors functions with we discussed in the Functions Associated with Burnup part of this tutorial. The parameter factor is used to calculate the fast neutron flux from the linear heat rate and generally has a value of 3e13 (n/(ms)/(W/m)). Finally, we instruct Bison to calculate the current fast neutron flux at the start of each time step with the execute_on = timestep_begin setting.

[./fast_neutron_flux]
  type = FastNeutronFluxAux
  variable = fast_neutron_flux
  block = clad
  rod_ave_lin_pow = power_history
  axial_power_profile = axial_peaking_factors
  factor = 3e13
  execute_on = timestep_begin
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The corresponding auxvariable, fast_neutron_flux is restricted to the cladding mesh block with block = clad. The shape functions for this auxvariable are defined with the order and family settings in the [GlobalParams] section at the top of our input file.

In a similar manner, we determine the fast neutron fluence with the FastNeutronFluenceAux auxkernel. Since we only use the fast neutron fluence in the cladding, we keep the auxkernel on the clad mesh block with block = clad. FastNeutronFluenceAux takes as an argument the name of the fast_neutron_flux auxvariable we just defined. As in the FastNeutronFluxAux auxkernel, we have Bison calculate this fast neutron fluence value at the start of each time step with the parameter setting execute_on = timestep_begin.


[./fast_neutron_fluence]
  type = FastNeutronFluenceAux
  variable = fast_neutron_fluence
  block = clad
  fast_neutron_flux = fast_neutron_flux
  execute_on = timestep_begin
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

We create the corresponding fast_neutron_fluence auxvariable in the same manner as the fast_neutron_flux auxvariable we previously discussed.


[./fast_neutron_fluence]
  block = clad
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Grain Radius

The grain radius is used by the fission gas release Sifgrs material to better capture the interaction among the fission gas release and grain boundaries in the fuel. We use the GrainRadiusAux auxkernel to compute the change of the grain radius within the ceramic fuel. Since this auxkernel is only relevant for the fuel pellets, we use the setting block = fuel. GrainRadiusAux takes the name of the grain radius auxvariable. We also couple to the nonlinear temperature variable with the setting temp = temperature. We also instruct Bison to compute this auxkernel on each linear iteration of the time step with execute_on = linear.


[./grain_radius]
  type = GrainRadiusAux
  block = fuel
  variable = grain_radius
  temp = temperature
  execute_on = linear
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

As with the neutron auxvariables, we create the corresponding grain_radius auxvariable and restrict it to the fuel mesh block, block = fuel and use the shape function order and family parameters from the [GlobalParams] section of input file. For the grain_radius auxvariable we also set the initial_condition to 1 m, or 10e-6 m.


[./grain_radius]
  block = fuel
  initial_condition = 10e-6
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Mechanical Material Data Output

Any material property declared as a Property in Bison can be added to the postprocessors for data output with the use of auxkernels and auxvariables; however, we must know the name of the material property that is used within the source code of Bison to use this approach. To output all of the following six mechanical material properties, we use the MaterialRealAux auxkernel to compute the average of the material property value over the entire element.

Since we are averaging over the element, we can use a lower order finite element shape function; we set order = CONSTANT and family = MONOMIAL for each auxvariable created below.

To store the zirconium creep strain rate in the cladding, we tell Bison to retrieve the material property with the setting property = creep_rate from the ZryCreepLimbackHoppeUpdate material. We also restrict this auxkernel to the cladding with block = clad.


[./creep_strain_rate]
  type = MaterialRealAux
  property = creep_rate
  variable = creep_strain_rate
  block = clad
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The corresponding auxvariable is named creep_strain_rate and is also restricted to the cladding mesh block.


[./creep_strain_rate]
  order = CONSTANT
  family = MONOMIAL
  block = clad
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

In setting up the parameters for the UO2VolumetricSwellingEigenstrain, we elected to save the solid swelling fission products. This material property is connected to this auxkernel with the setting property = solid_swelling. This auxkernel is restricted to the fuel mesh block with block = fuel. Since we use this auxkernel to collect a computed material property, we set execute_on = timestep_end to run this auxkernel at the end of the time step.


[./solid_swell]
  type = MaterialRealAux
  variable = solid_swell
  property = solid_swelling
  execute_on = timestep_end
  block = fuel
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The corresponding auxvariable, solid_swell, is defined only on the fuel mesh block with block = fuel.


[./solid_swell]
  order = CONSTANT
  family = MONOMIAL
  block = fuel
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The UO2VolumetricSwellingEigenstrain also creates the material property of gas fission products. We tell Bison to retrieve this material property with the setting property = gas_swelling, restrict the auxkernel to operate on the fuel with block = fuel, and run the auxkernel at the end of the time step with the execute_on = timestep_end setting.


[./gas_swell]
  type = MaterialRealAux
  variable = gas_swell
  property = gas_swelling
  execute_on = timestep_end
  block = fuel
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

And the corresponding auxvariable is defined as gas_swell with the same settings as the solid_swell auxvariable.


[./gas_swell]
  order = CONSTANT
  family = MONOMIAL
  block = fuel
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

In the same fashion we create an auxkernel to retrieve the mechanical material property with the setting property = densification


[./densification]
  type = MaterialRealAux
  variable = densification
  property = densification
  execute_on = timestep_end
  block = fuel
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

and create the densification auxvariable to complete the pair:


[./densification]
  order = CONSTANT
  family = MONOMIAL
  block = fuel
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Finally, we instruct Bison to retrieve the total volumetric eigenstrain material property with the auxkernel setting property = volumetric_swelling_strain. We use the same execute_on and block settings as for the other auxkernels storing UO2VolumetricSwellingEigenstrain material properties


[./volumetric_swelling_strain]
  type = MaterialRealAux
  variable = volumetric_swelling_strain
  property = volumetric_swelling_strain
  execute_on = timestep_end
  block = fuel
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

and the corresponding volumetric_swelling_strain auxvariable.


[./volumetric_swelling_strain]
  order = CONSTANT
  family = MONOMIAL
  block = fuel
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The UO2RelocationEigenstrain creates the material property relocation eigenstrain. We use the MaterialRealAux setting property = relocation_strain to retrieve this eigenstrain material property. As with the volumetric eigenstrain auxkernels, we restrict this auxkernel to only the fuel pellets with block = fuel, and we instruct Bison to run this auxkernel at the end of the timesteps with execute_on = timestep_end.


[./relocation_strain]
  type = MaterialRealAux
  variable = relocation
  property = relocation_strain
  execute_on = timestep_end
  block = fuel
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The corresponding relocation auxvariable is created in a similar manner and is also restricted to the fuel mesh block with block = fuel.


[./relocation]
  order = CONSTANT
  family = MONOMIAL
  block = fuel
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Postprocessors

Postprocessors are used to output simulation data to an Exodus file or a CSV file, as specified in the [Outputs] input file section. The example input file we are using in this tutorial includes too many postprocessors to work through each in this tutorial. We cover three different sets of postprocessors in this tutorial below.

Standard Bison PostProcessors

The average temperature within the plenum space, between the fuel pellet exterior and the cladding interior, is used by PlenumPressure. We use the specialized SideAverageValueLayered1D postprocessor to compute this value. As shown in the example below, we use the parameter setting boundary = 9 to calculate the average from the plenum on the combination of the pellet exterior and the cladding interior surfaces. For more information on the boundary settings used in Bison, refer to the Mesh Script documentation page.

To find the average temperature value, we set variable = temperature, and we provide the required geometry information to the postprocessor with the Layered1DFuelPinGeometry by setting fuel_pin_geometry = pin_geometry.

Because the ave_temp_interior is used in the PlenumPressure input file section, we must tell Bison to compute the average plenum temperature at both the start of the time step and on every linear iteration with the setting execute_on = initial linear.

[./ave_temp_interior] # average temperature of the cladding interior and all pellet exteriors
  type = SideAverageValueLayered1D
  boundary = 9
  variable = temperature
  execute_on = 'initial linear'
  fuel_pin_geometry = pin_geometry
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The volume of the plenum space, or the volume of the plenum gas, is also used by the PlenumPressure. We use another specialized postprocessor, InternalVolumeLayered1D to calculate this volume value. As with the ave_temp_interior postprocessor, we use the boundary = 9 and execute_on = 'initial linear' parameter settings.

Because we have a single nonlinear displacement variable, we set the parameter to component = 0. We also need to provide the InternalVolumeLayered1D postprocessor with information about the out-of-plane scalar strain variables with the parameter out_of_plane_strain = strain_yy setting.


[./gas_volume]
  type = InternalVolumeLayered1D
  boundary = 9
  execute_on = 'initial linear'
  component = 0
  out_of_plane_strain = strain_yy
  fuel_pin_geometry = pin_geometry
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The amount of fission gas products released into the PlenumPressure and the GapHeatTransferLWR boundary condition sections of the input file. We use the ElementIntegralFisGasReleasedSifgrs postprocessor to collect the amount of fission gas released from the Sifgrs.md material model. We couple the nonlinear temperature variable, restrict the post processor to the fuel mesh block, and provide information about the fuel rod geometry with the fuel_pin_geometry user object.


[./fis_gas_released] # fission gas released to plenum (moles)
  type = ElementIntegralFisGasReleasedSifgrs
  variable = temperature
  block = fuel
  fuel_pin_geometry = pin_geometry
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The next four postprocessors are also associated with the Sifgrs material. Search for the postprocessor name given as the argument to the type parameters in our Bison documentation list.


[./fis_gas_produced] # fission gas produced (moles)
  type = ElementIntegralFisGasGeneratedSifgrs
  variable = temperature
  block = fuel
  fuel_pin_geometry = pin_geometry
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

[./fis_gas_grain]
  type = ElementIntegralFisGasGrainSifgrs
  variable = temperature
  block = fuel
  outputs = exodus
  fuel_pin_geometry = pin_geometry
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

[./fis_gas_boundary]
  type = ElementIntegralFisGasBoundarySifgrs
  variable = temperature
  block = fuel
  outputs = exodus
  fuel_pin_geometry = pin_geometry
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

[./fis_gas_released] # fission gas released to plenum (moles)
  type = ElementIntegralFisGasReleasedSifgrs
  variable = temperature
  block = fuel
  fuel_pin_geometry = pin_geometry
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Finally, we use the ElementIntegralPower postprocessor to verify the calculation of power by Bison. We pass information about the nonlinear temperature variable, burnup, and fuel rod geometry with the variable, burnup_function, and fuel_pin_geometry parameters. We also restrict this postprocessor to operate only on the fuel mesh block.


[./rod_total_power]
  type = ElementIntegralPower
  variable = temperature
  burnup_function = burnup
  block = fuel
  fuel_pin_geometry = pin_geometry
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Mechanical Materials PostProcessors

To add all of the mechanical material property auxvariables we created in the AuxKernels & AuxVariables part of the tutorial. For all of these material properties, we use the ElementAverageValue postprocessor. In each of the six examples shown below, we use the name of the auxvariable created above as the argument to the variable parameter. We use the same mesh block restriction as set in the creation of the corresponding auxkernal and auxvariable pairs.


[./effective_creep_strain_rate]
  type = ElementAverageValue
  variable = creep_strain_rate
  block = clad
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

[./solid_swelling]
  type = ElementAverageValue
  variable = solid_swell
  block = fuel
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

[./gas_swelling]
  type = ElementAverageValue
  variable = gas_swell
  block = fuel
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

[./densification]
  type = ElementAverageValue
  variable = densification
  block = fuel
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

[./volumetric_swelling]
  type = ElementAverageValue
  variable = volumetric_swelling_strain
  block = fuel
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

[./relocation]
  type = ElementAverageValue
  variable = relocation
  block = fuel
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Additional PostProcessors Examples

This input file has many more postprocessors than a standard Bison fuel rod analysis simulation because this input file is used to compare to a 2D-RZ geometry where more detailed information is required. In this part of the tutorial, we will go through only three additional postprocessors as examples. Refer to the complete input file for the complete set of postprocessors used in our Layered 1D example.

We use ElementAverageValue postprocessor to compute and store the average value of the temperature in the fuel. Restricting this postprocessor to the fuel mesh block with block = fuel, we tell Bison to compute the average value of the nonlinear variable temperature across the entire element.


[./ave_fuel_temp]
  type = ElementAverageValue
  block = fuel
  variable = temperature
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

The NodalVariableValue is used in this example postprocessor section to output the value of the nonlinear variable temperature at a specific node. We tell Bison what the number of this specific node is with the node_id parameter.

note:Bison Node ID is One Less than Paraview GlobalID value

Recall, if you are using Paraview to determine the ID numbers of specific nodes, the Bison mesh numbering scheme is a value one less than the value of the Paraview Global ID value.


[./central_fuel_temp]
  type = NodalVariableValue
  nodeid = 262 #Mesh dependent (0.0041, 0.05661)
  variable = temperature
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Our final example uses the NodalExtremeValue postprocessor to output the maximum value of the nonlinear variable temperature on the fuel mesh block. We tell Bison to output the maximum value with the setting value_type = max.


[./max_fuel_temp]
  type = NodalExtremeValue
  block = fuel
  value_type = max
  variable = temperature
[../]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

VectorPostprocessors

The NodalValueSampler is used to output the values of the radial displacement as a vector along the length of the fuel rod. The boundary settings are used to output the displacement values along the outer cladding surface (2) and the exterior fuel pellet surface (10). For more information on the boundary settings used in Bison, refer to the Mesh Script documentation page.

To sort the radial displacements by axial position of the fuel rod, we set the parameter sort_by = y. We give each vectorpostprocessor a unique name with the outputs parameters in the two input file sections.


[VectorPostprocessors]
  [./clad]
    type = NodalValueSampler
    variable = disp_x
    boundary = 2
    sort_by = y
    outputs = 'clad_radial_displacement'
  [../]
  [./pellet]
    type = NodalValueSampler
    variable = disp_x
    boundary = 10
    sort_by = y
    outputs = 'fuel_radial_displacement'
  [../]
[]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Outputs

The standard Bison output settings are shown in the example below. Note the match in the names of the sub-outputs with the arguments of the vectorpostprocessor outputs parameters.


[Outputs]
  exodus = true
  csv = true
  color = false
  print_perf_log = true
  [./clad_radial_displacement]
    type = CSV
    execute_on = 'FINAL'
  [../]
  [./fuel_radial_displacement]
    type = CSV
    execute_on = 'FINAL'
  [../]
[]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Solver Options

Nearly all Bison simulations used the Transient option since the simulations involve the evolution of temperature and displacements over time. We also use the default solve_type = PJFNK.

There are numerous PETSc options available, and discussing all of the options is beyond the scope of this tutorial. The Bison team has observed good convergence behavior when using the superlu_dist package for simulations involving contact but run on a small number of processors. Since the Layered 1D problems run quickly on a single processor, they are an ideal use case for the superlu_dist option.

The Transient checks both the linear and nonlinear tolerances for a time step to determine if the solve is convergence. The parameters l_tol and l_max_its control the linear iterations, and the nl_rel_tol, nl_abs_tol, and nl_max_its control the nonlinear iterations.

As we discussed in the Functions Associated with Pressure part of this tutorial, when building the pressure_ramp PiecewiseLinear function, we set this Bison simulation to begin at start_time = -200 and to take a single start up step with n_startup_steps = 1. The end_time parameter tells Bison for how long the simulation should run.

The parameters dtmax and dtmin provide Bison with maximum and minimum limits on the size of the time step size. While Bison will prevent the simulation from taking time step sizes larger than dtmax, the simulation will fail if Bison attempts to take a time step smaller than dtmin.


[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'

  line_search = 'none'

  l_max_its = 50
  l_tol = 8e-3
  nl_max_its = 25
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-7

  start_time = -200
  n_startup_steps = 1
  end_time = 8.0e7
  dtmax = 2e6
  dtmin = 1

  [./TimeStepper]
    type = IterationAdaptiveDT
    dt = 2e2
    optimal_iterations = 8
    iteration_window = 2
    growth_factor = 2
    cutback_factor = .5
  [../]
[]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

In addition to the limits on the size of the time step, we also impose limitations on the size of the nonlinear variable changes during a single time step. We apply these limits with the MaxIncrement damper as shown below. The value of the max_increment parameter depends on which nonlinear variable the MaxIncrement damper is used on.


[Dampers]
  [./limitT]
    type = MaxIncrement
    max_increment = 100.0
    variable = temperature
  [../]
  [./limitX]
    type = MaxIncrement
    max_increment = 1e-5
    variable = disp_x
  [../]
[]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)

Conclusion

Within this tutorial we have constructed the our input file in an order that best matches the natural discussion of the physics of nuclear fuel rod analysis; however, this tutorial order deviates from the standard input file order recommended by the Bison team. Refer to the complete example input file for the standard Bison input file order.

# Model is of a 10 pellet stack of fuel modeled in 1.5d

[GlobalParams]
  density = 10431.0
  initial_porosity = 0.05
  order = SECOND
  family = LAGRANGE
  energy_per_fission = 3.2e-11  # J/fission
  displacements = disp_x
  temperature = temperature
[]

[Problem]
  # Specify coordinate system type
  coord_type = RZ
[]

[Mesh]
  type = Layered1DMesh
  patch_update_strategy = auto
  partitioner = centroid
  centroid_partitioner_direction = y
  slices_per_block = 10
  clad_gap_width = 8.0e-5
  clad_thickness = 0.00056
  fuel_height = 0.1186
  plenum_height = 0.027
[]

[UserObjects]
  [./pin_geometry]
    type = Layered1DFuelPinGeometry
  [../]
[]

[Variables]
  [./temperature]
    initial_condition = 580.0     # set initial temperature to coolant inlet
  [../]
[]

[AuxVariables]
  [./disp_y] ## Required for easier visualization in Paraview
  [../]
  [./disp_z] ## Required for easier visualization in Paraview
  [../]
  [./fast_neutron_flux]
    block = clad
  [../]
  [./fast_neutron_fluence]
    block = clad
  [../]
  [./grain_radius]
    block = fuel
    initial_condition = 10e-6
  [../]
  [./creep_strain_rate]
    order = CONSTANT
    family = MONOMIAL
    block = clad
  [../]
  [./creep_strain]
    order = CONSTANT
    family = MONOMIAL
    block = clad
  [../]
  [./solid_swell]
    order = CONSTANT
    family = MONOMIAL
    block = fuel
  [../]
  [./gas_swell]
    order = CONSTANT
    family = MONOMIAL
    block = fuel
  [../]
  [./densification]
    order = CONSTANT
    family = MONOMIAL
    block = fuel
  [../]
  [./volumetric_swelling_strain]
    order = CONSTANT
    family = MONOMIAL
    block = fuel
  [../]
  [./relocation]
    order = CONSTANT
    family = MONOMIAL
    block = fuel
  [../]
[]

[Functions]
  [./power_history]
    type = PiecewiseLinear   # reads and interpolates an input file containing rod average linear power vs time
    data_file = powerhistory.csv
    scale_factor = 1
  [../]
  [./axial_peaking_factors]      # reads and interpolates an input file containing the axial power profile vs time
    type = PiecewiseBilinear
    data_file = peakingfactors.csv
    scale_factor = 1
    axis = 1 # (0,1,2) => (x,y,z)
  [../]
  [./q]
    type = CompositeFunction
    functions = 'power_history axial_peaking_factors'
  [../]
  [./pressure_ramp]              # reads and interpolates input data defining amplitude curve for fill gas pressure
    type = PiecewiseLinear
    x = '-200 0'
    y = '0 1'
  [../]
  [./clad_axial_pressure]
    type = CladdingAxialPressureFunction
    plenum_pressure = plenum_pressure
    coolant_pressure = pressure_ramp
    coolant_pressure_scaling_factor = 15.5e6
    fuel_pin_geometry = pin_geometry
  [../]
  [./fuel_axial_pressure]
    type = ParsedFunction
    value = plenum_pressure
    vars = plenum_pressure
    vals = plenum_pressure
  [../]
[]

[Kernels]
  [./heat]         # gradient term in heat conduction equation
    type = HeatConduction
    variable = temperature
  [../]
  [./heat_ie]       # time term in heat conduction equation
    type = HeatConductionTimeDerivative
    variable = temperature
  [../]
  [./heat_source]  # source term in heat conduction equation
    type = NeutronHeatSource
    variable = temperature
    block = fuel # fission rate applied to the fuel (block 2) only
    burnup_function = burnup
  [../]
[]

[Modules]
  [./TensorMechanics]
    [./Layered1DMaster]
      [./fuel]
        block = fuel
        add_variables = true
        strain = FINITE
        add_scalar_variables = true
        out_of_plane_strain_name = strain_yy
        fuel_pin_geometry = pin_geometry
        out_of_plane_pressure = fuel_axial_pressure
        eigenstrain_names = 'fuelthermal_strain swelling_strain fuel_relocation_strain'
        generate_output = 'stress_xx stress_yy stress_zz vonmises_stress strain_xx'
      [../]
      [./clad]
        block = clad
        add_variables = true
        strain = FINITE
        add_scalar_variables = true
        out_of_plane_strain_name = strain_yy
        fuel_pin_geometry = pin_geometry
        out_of_plane_pressure = clad_axial_pressure
        eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_strain'
        generate_output = 'stress_xx stress_yy stress_zz vonmises_stress strain_xx'
      [../]
    [../]
  [../]
[]

[Burnup]
  [./burnup]
    block = fuel
    rod_ave_lin_pow = power_history          # using the power function defined above
    axial_power_profile = axial_peaking_factors     # using the axial power profile function defined above
    num_radial = 80
    num_axial = 11
    order = CONSTANT
    family = MONOMIAL
    fuel_pin_geometry = pin_geometry
    fuel_volume_ratio = 1.0 # for use with dished pellets (ratio of actual volume to cylinder volume)
    RPF = RPF
    i_enrich = '0.05 0.95 0 0 0 0'
  [../]
[]

[AuxKernels]
  [./fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    block = clad
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    factor = 3e13
    execute_on = timestep_begin
  [../]
  [./fast_neutron_fluence]
    type = FastNeutronFluenceAux
    variable = fast_neutron_fluence
    block = clad
    fast_neutron_flux = fast_neutron_flux
    execute_on = timestep_begin
  [../]
  [./grain_radius]
    type = GrainRadiusAux
    block = fuel
    variable = grain_radius
    temp = temperature
    execute_on = linear
  [../]
  [./creep_strain]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = creep_strain
    block = clad
    execute_on = timestep_end
  [../]
  [./creep_strain_rate]
    type = MaterialRealAux
    property = creep_rate
    variable = creep_strain_rate
    block = clad
  [../]
  [./solid_swell]
    type =  MaterialRealAux
    variable = solid_swell
    property = solid_swelling
    execute_on = timestep_end
    block = fuel
  [../]
  [./gas_swell]
    type = MaterialRealAux
    variable = gas_swell
    property = gas_swelling
    execute_on = timestep_end
    block = fuel
  [../]
  [./densification]
    type = MaterialRealAux
    variable = densification
    property = densification
    execute_on = timestep_end
    block = fuel
  [../]
  [./volumetric_swelling_strain]
    type = MaterialRealAux
    variable = volumetric_swelling_strain
    property = volumetric_swelling_strain
    execute_on = timestep_end
    block = fuel
  [../]
  [./relocation_strain]
    type = MaterialRealAux
    variable = relocation
    property = relocation_strain
    execute_on = timestep_end
    block = fuel
  [../]
[]

[Contact]
  [./pellet_clad_mechanical]
    master = 5
    slave = 10
    formulation = kinematic
    model = frictionless
    system = Constraint
    penalty = 1e7
  [../]
[]

[ThermalContact]
  [./thermal_contact]
    type = GapHeatTransferLWR
    variable = temperature
    master = 5
    slave = 10
    initial_moles = initial_moles       # coupling to a postprocessor which supplies the initial plenum/gap gas mass
    gas_released = fis_gas_released     # coupling to a postprocessor which supplies the fission gas addition
    contact_pressure = contact_pressure
  [../]
[]

[BCs]
  [./no_x_all] # pin pellets and clad along axis of symmetry (y)
    type = PresetBC
    variable = disp_x
    boundary = 12
    value = 0.0
  [../]
  [./Pressure] # apply coolant pressure on clad outer walls
    [./coolantPressure]
      boundary = 2
      function = pressure_ramp   # use the pressure_ramp function defined above
      factor = 15.5e6
    [../]
  [../]
  [./PlenumPressure] # apply plenum pressure on clad inner walls and pellet surfaces
    [./plenumPressure]
      boundary = 9
      initial_pressure = 2.0e6
      startup_time = 0
      R = 8.314
      output_initial_moles = initial_moles       # coupling to post processor to get initial fill gas mass
      temperature = ave_temp_interior            # coupling to post processor to get gas temperature approximation
      volume = gas_volume                        # coupling to post processor to get gas volume
      material_input = fis_gas_released          # coupling to post processor to get fission gas added
      output = plenum_pressure                   # coupling to post processor to output plenum/gap pressure
    [../]
  [../]
[]

[CoolantChannel]
  [./convective_clad_surface] # apply convective boundary to clad outer surface
    variable = temperature
    boundary = 2
    inlet_temperature = 580      # K
    inlet_pressure = 15.5e6   # Pa
    inlet_massflux = 3800     # kg/m^2-sec
    rod_diameter = 0.948e-2 # m
    rod_pitch = 1.26e-2  # m
    linear_heat_rate  = power_history
    axial_power_profile = axial_peaking_factors
  [../]
[]

[Materials]
  [./fuel_thermal]
    type = ThermalFuel
    block = fuel
    model = 4 # NFIR
    temp = temperature
    burnup_function = burnup
  [../]
  [./fuel_density]
    type = Density
    block = fuel
  [../]

  [./fuel_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = fuel
    youngs_modulus = 2.0e11
    poissons_ratio = 0.345
  [../]
  [./fuel_stress]
    type = ComputeFiniteStrainElasticStress
    block = fuel
  [../]
  [./fuel_thermal_strain]
    type = ComputeThermalExpansionEigenstrain
    block = fuel
    thermal_expansion_coeff = 10.0e-6
    stress_free_temperature = 295.0
    eigenstrain_name = fuelthermal_strain
  [../]
  [./fuel_swelling]
    type = UO2VolumetricSwellingEigenstrain
    block = fuel
    gas_swelling_model_type = SIFGRS
    burnup_function = burnup
    initial_fuel_density = 10431.0
    save_densification = true
    save_solid_swelling = true
    eigenstrain_name = swelling_strain
  [../]
  [./fuel_relocation]
    type = UO2RelocationEigenstrain
    block = fuel
    burnup_function = burnup
    fuel_pin_geometry = pin_geometry
    linear_heat_rate_function = q
    relocation_activation1 = 5000.0
    burnup_relocation_stop = 0.024
    relocation_model = ESCORE_modified
    eigenstrain_name = fuel_relocation_strain
  [../]
  [./fission_gas_release]
    type = Sifgrs
    block = fuel
    temp = temperature
    burnup_function = burnup
    gbs_model = true
    grain_radius = grain_radius
  [../]

  [./clad_thermal]
    type = HeatConductionMaterial
    block = clad
    thermal_conductivity = 16.0
    specific_heat = 330.0
  [../]
  [./clad_density]
    type = Density
    block = clad
    density = 6551.0
  [../]

  [./clad_elasticity_tensor]
    type = ZryElasticityTensor
    block = clad
  [../]
  [./clad_stress]
    type = ComputeMultipleInelasticStress
    block = clad
    tangent_operator = elastic
    inelastic_models = 'zrycreep'
  [../]
  [./zrycreep]
    type = ZryCreepLimbackHoppeUpdate
    block = clad
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    model_irradiation_creep = true
    model_primary_creep = true
    model_thermal_creep = true
    zircaloy_material_type = stress_relief_annealed
  [../]
  [./clad_thermal_expansion]
    type = ZryThermalExpansionMatproEigenstrain
    block = clad
    stress_free_temperature = 295.0
    eigenstrain_name = clad_thermal_eigenstrain
  [../]
  [./clad_irradiation_swelling]
    type = ZryIrradiationGrowthEigenstrain
    block = clad
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_model_type = stress_relief_annealed
    eigenstrain_name = clad_irradiation_strain
  [../]
[]

[Dampers]
  [./limitT]
    type = MaxIncrement
    max_increment = 100.0
    variable = temperature
  [../]
  [./limitX]
    type = MaxIncrement
    max_increment = 1e-5
    variable = disp_x
  [../]
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu superlu_dist'

  line_search = 'none'

  l_max_its = 50
  l_tol = 8e-3
  nl_max_its = 25
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-7

  start_time = -200
  n_startup_steps = 1
  end_time = 8.0e7
  dtmax = 2e6
  dtmin = 1

  [./TimeStepper]
    type = IterationAdaptiveDT
    dt = 2e2
    optimal_iterations = 8
    iteration_window = 2
    growth_factor = 2
    cutback_factor = .5
  [../]
[]

[Postprocessors]
  [./ave_temp_interior]            # average temperature of the cladding interior and all pellet exteriors
    type = SideAverageValueLayered1D
    boundary = 9
    variable = temperature
    execute_on = 'initial linear'
    fuel_pin_geometry = pin_geometry
  [../]
  [./clad_inner_vol]              # volume inside of cladding
    type = InternalVolumeLayered1D
    boundary = 7
    component = 0
    fuel_pin_geometry = pin_geometry
    out_of_plane_strain = strain_yy
    execute_on = 'initial linear'
    #outputs = exodus
  [../]
  [./pellet_volume]               # fuel pellet total volume
    type = InternalVolumeLayered1D
    boundary = 8
#    scale_factor = -1
    component = 0
    fuel_pin_geometry = pin_geometry
    out_of_plane_strain = strain_yy
    execute_on = 'initial linear'
    #outputs = exodus
  [../]
  [./avg_clad_temp]               # average temperature of cladding interior
    type = SideAverageValueLayered1D
    boundary = 7
    variable = temperature
    fuel_pin_geometry = pin_geometry
    execute_on = 'initial linear'
  [../]
  [./fis_gas_produced]           # fission gas produced (moles)
    type = ElementIntegralFisGasGeneratedSifgrs
    variable = temperature
    block = fuel
    fuel_pin_geometry = pin_geometry
  [../]
  [./fis_gas_released]           # fission gas released to plenum (moles)
    type = ElementIntegralFisGasReleasedSifgrs
    variable = temperature
    block = fuel
    fuel_pin_geometry = pin_geometry
  [../]
  [./fis_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    variable = temperature
    block = fuel
    outputs = exodus
    fuel_pin_geometry = pin_geometry
  [../]
  [./fis_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    variable = temperature
    block = fuel
    outputs = exodus
    fuel_pin_geometry = pin_geometry
  [../]
  [./fisson_gas_release]
    type = FGRPercent
    value1 = fis_gas_released
    value2 = fis_gas_produced
  [../]

  [./gas_volume]
    type = InternalVolumeLayered1D
    boundary = 9
    execute_on = 'initial linear'
    component = 0
    out_of_plane_strain = strain_yy
    fuel_pin_geometry = pin_geometry
  [../]
  [./flux_from_clad]           # area integrated heat flux from the cladding
    type = SideFluxIntegralLayered1D
    variable = temperature
    boundary = 5
    diffusivity = thermal_conductivity
    fuel_pin_geometry = pin_geometry
  [../]
  [./flux_from_fuel]          # area integrated heat flux from the fuel
    type = SideFluxIntegralLayered1D
    variable = temperature
    boundary = 10
    diffusivity = thermal_conductivity
    fuel_pin_geometry = pin_geometry
  [../]
  [./_dt]                     # time step
    type = TimestepSize
  [../]
  [./num_lin_it]
    type = NumLinearIterations
  [../]
  [./num_nonlin_it]
    type = NumNonlinearIterations
  [../]
  [./tot_lin_it]
    type = CumulativeValuePostprocessor
    postprocessor = num_lin_it
  [../]
  [./tot_nonlin_it]
    type = CumulativeValuePostprocessor
    postprocessor = num_nonlin_it
  [../]
  [./alive_time]
    type = PerformanceData
    event = ALIVE
  [../]

  [./rod_total_power]
    type = ElementIntegralPower
    variable = temperature
    burnup_function = burnup
    block = fuel
    fuel_pin_geometry = pin_geometry
  [../]
  [./rod_input_power]
    type = FunctionValuePostprocessor
    function = power_history
    scale_factor = 0.1186 # rod height
  [../]

  [./ave_fuel_temp]
    type = ElementAverageValue
    block = fuel
    variable = temperature
  [../]
  [./central_fuel_temp]
    type = NodalVariableValue
    nodeid = 262 #Mesh dependent (0.0041, 0.05661)
    variable = temperature
  [../]
  [./max_fuel_temp]
    type = NodalExtremeValue
    block = fuel
    value_type = max
    variable = temperature
  [../]
  [./max_clad_temp]
    type = NodalExtremeValue
    block = clad
    value_type = max
    variable = temperature
  [../]

  ### Comparisons for 1.5D work, mesh specific ####################    # von Mises Stress
  [./top_vonMises_fuel]
    type = ElementalVariableValue
    elementid = 171 # mesh dependent (contains pt. 0.0041, 0.09219)
    variable = vonmises_stress
  [../]
  [./center_vonMises_fuel]
    type = ElementalVariableValue
    elementid = 123 # mesh dependent (contains pt. 0.0041, 0.05661)
    variable = vonmises_stress
  [../]
  [./bottom_vonMises_fuel]
    type = ElementalVariableValue
    elementid = 75 # mesh dependent (contains pt. 0.0041, 0.02103)
    variable = vonmises_stress
  [../]
  [./average_vonMises_fuel]
    type = ElementAverageValue
    variable = vonmises_stress
    block = fuel
  [../]
  [./top_vonMises_clad_inner]
    type = ElementalVariableValue
    elementid = 28 # mesh dependent (contains pt. 0.00418, 0.09219)
    variable = vonmises_stress
  [../]
  [./top_vonMises_clad_outer]
    type = ElementalVariableValue
    elementid = 31 # mesh dependent (contains pt. 0.00474, 0.09219)
    variable = vonmises_stress
  [../]
  [./center_vonMises_clad_inner]
    type = ElementalVariableValue
    elementid = 16 # mesh dependent (contains pt. 0.00418, 0.05661)
    variable = vonmises_stress
  [../]
  [./center_vonMises_clad_outer]
    type = ElementalVariableValue
    elementid = 19 # mesh dependent (contains pt. 0.00474, 0.05661)
    variable = vonmises_stress
  [../]
  [./bottom_vonMises_clad_inner]
    type = ElementalVariableValue
    elementid = 4 # mesh dependent (contains pt. 0.00418, 0.02103)
    variable = vonmises_stress
  [../]
  [./bottom_vonMises_clad_outer]
    type = ElementalVariableValue
    elementid = 7 # mesh dependent (contains pt. 0.00474, 0.02103)
    variable = vonmises_stress
  [../]
  [./average_vonMises_clad]
    type = ElementAverageValue
    variable = vonmises_stress
    block = clad
  [../]

  # radial stress
  [./top_stress_rr_fuel]
    type = ElementalVariableValue
    elementid = 171 # mesh dependent (contains pt. 0.0041, 0.09219)
    variable = stress_xx
  [../]
  [./center_stress_rr_fuel]
    type = ElementalVariableValue
    elementid = 123 # mesh dependent (contains pt. 0.0041, 0.05661)
    variable = stress_xx
  [../]
  [./bottom_stress_rr_fuel]
    type = ElementalVariableValue
    elementid = 75 # mesh dependent (contains pt. 0.0041, 0.02103)
    variable = stress_xx
  [../]
  [./average_stress_rr_fuel]
    type = ElementAverageValue
    variable = stress_xx
    block = fuel
  [../]
  [./top_stress_rr_clad_inner]
    type = ElementalVariableValue
    elementid = 28 # mesh dependent (contains pt. 0.00418, 0.09219)
    variable = stress_xx
  [../]
  [./top_stress_rr_clad_outer]
    type = ElementalVariableValue
    elementid = 31 # mesh dependent (contains pt. 0.00474, 0.09219)
    variable = stress_xx
  [../]
  [./center_stress_rr_clad_inner]
    type = ElementalVariableValue
    elementid = 16 # mesh dependent (contains pt. 0.00418, 0.05661)
    variable = stress_xx
  [../]
  [./center_stress_rr_clad_outer]
    type = ElementalVariableValue
    elementid = 19 # mesh dependent (contains pt. 0.00474, 0.05661)
    variable = stress_xx
  [../]
  [./bottom_stress_rr_clad_inner]
    type = ElementalVariableValue
    elementid = 4 # mesh dependent (contains pt. 0.00418, 0.02103)
    variable = stress_xx
  [../]
  [./bottom_stress_rr_clad_outer]
    type = ElementalVariableValue
    elementid = 7 # mesh dependent (contains pt. 0.00474, 0.02103)
    variable = stress_xx
  [../]
  [./average_stress_rr_clad]
    type = ElementAverageValue
    variable = stress_xx
    block = clad
  [../]

  # radial strain
  [./top_strain_rr_fuel]
    type = ElementalVariableValue
    elementid = 171 # mesh dependent (contains pt. 0.0041, 0.09219)
    variable = strain_xx
  [../]
  [./center_strain_rr_fuel]
    type = ElementalVariableValue
    elementid = 123 # mesh dependent (contains pt. 0.0041, 0.05661)
    variable = strain_xx
  [../]
  [./bottom_strain_rr_fuel]
    type = ElementalVariableValue
    elementid = 75 # mesh dependent (contains pt. 0.0041, 0.02103)
    variable = strain_xx
  [../]
  [./average_strain_rr_fuel]
    type = ElementAverageValue
    variable = strain_xx
    block = fuel
  [../]
  [./top_strain_rr_clad_inner]
    type = ElementalVariableValue
    elementid = 28 # mesh dependent (contains pt. 0.00418, 0.09219)
    variable = strain_xx
  [../]
  [./top_strain_rr_clad_outer]
    type = ElementalVariableValue
    elementid = 31 # mesh dependent (contains pt. 0.00474, 0.09219)
    variable = strain_xx
  [../]
  [./center_strain_rr_clad_inner]
    type = ElementalVariableValue
    elementid = 16 # mesh dependent (contains pt. 0.00418, 0.05661)
    variable = strain_xx
  [../]
  [./center_strain_rr_clad_outer]
    type = ElementalVariableValue
    elementid = 19 # mesh dependent (contains pt. 0.00474, 0.05661)
    variable = strain_xx
  [../]
  [./bottom_strain_rr_clad_inner]
    type = ElementalVariableValue
    elementid = 4 # mesh dependent (contains pt. 0.00418, 0.02103)
    variable = strain_xx
  [../]
  [./bottom_strain_rr_clad_outer]
    type = ElementalVariableValue
    elementid = 7 # mesh dependent (contains pt. 0.00474, 0.02103)
    variable = strain_xx
  [../]
  [./average_strain_rr_clad]
    type = ElementAverageValue
    variable = strain_xx
    block = clad
  [../]

  # effective creep strain
  [./top_creep_strain_clad_inner]
    type = ElementalVariableValue
    elementid = 28 # mesh dependent (contains pt. 0.00418, 0.09219)
    variable = creep_strain
  [../]
  [./top_creep_strain_clad_outer]
    type = ElementalVariableValue
    elementid = 31 # mesh dependent (contains pt. 0.00474, 0.09219)
    variable = creep_strain
  [../]
  [./center_creep_strain_clad_inner]
    type = ElementalVariableValue
    elementid = 16 # mesh dependent (contains pt. 0.00418, 0.05661)
    variable = creep_strain
  [../]
  [./center_creep_strain_clad_outer]
    type = ElementalVariableValue
    elementid = 19 # mesh dependent (contains pt. 0.00474, 0.05661)
    variable = creep_strain
  [../]
  [./bottom_creep_strain_clad_inner]
    type = ElementalVariableValue
    elementid = 4 # mesh dependent (contains pt. 0.00418, 0.02103)
    variable = creep_strain
  [../]
  [./bottom_creep_strain_clad_outer]
    type = ElementalVariableValue
    elementid = 7 # mesh dependent (contains pt. 0.00474, 0.02103)
    variable = creep_strain
  [../]
  [./average_creep_strain_clad]
    type = ElementAverageValue
    variable = creep_strain
    block = clad
  [../]

  ### Nodal displacements
  [./top_disp_r_fuel]
    type = NodalVariableValue
    variable = disp_x
    nodeid = 361 # mesh dependent, at (0.0041, 0.09219)
  [../]
  [./center_disp_r_fuel]
    type = NodalVariableValue
    variable = disp_x
    nodeid = 262 # mesh dependent, at (0.0041, 0.05661)
  [../]
  [./bottom_disp_r_fuel]
    type = NodalVariableValue
    variable = disp_x
    nodeid = 163 # mesh dependent, at (0.0041, 0.02103)
  [../]
  [./top_disp_r_clad_inner]
    type = NodalVariableValue
    variable = disp_x
    nodeid = 63 #mesh dependent, at (0.00418, 0.09219)
  [../]
  [./top_disp_r_clad_outer]
    type = NodalVariableValue
    variable = disp_x
    nodeid = 68 #mesh dependent, at (0.00474, 0.09219)
  [../]
  [./center_disp_r_clad_inner]
    type = NodalVariableValue
    variable = disp_x
    nodeid = 36 #mesh dependent, at (0.00418, 0.05661)
  [../]
  [./center_disp_r_clad_outer]
    type = NodalVariableValue
    variable = disp_x
    nodeid = 43 #mesh dependent, at (0.00474, 0.05661)
  [../]
  [./bottom_disp_r_clad_inner]
    type = NodalVariableValue
    variable = disp_x
    nodeid = 9 #mesh dependent, at (0.00418, 0.02103)
  [../]
  [./bottom_disp_r_clad_outer]
    type = NodalVariableValue
    variable = disp_x
    nodeid = 16 #mesh dependent, at (0.00418, 0.02103)
  [../]

  ### Nodal temperatures
  [./top_temp_fuel]
    type = NodalVariableValue
    variable = temperature
    nodeid = 361 # mesh dependent, at (0.0041, 0.09219)
  [../]
  [./center_temp_fuel]
    type = NodalVariableValue
    variable = temperature
    nodeid = 262 # mesh dependent, at (0.0041, 0.05661)
  [../]
  [./bottom_temp_fuel]
    type = NodalVariableValue
    variable = temperature
    nodeid = 163 # mesh dependent, at (0.0041, 0.02103)
  [../]
  [./top_temp_clad_inner]
    type = NodalVariableValue
    variable = temperature
    nodeid = 63 #mesh dependent, at (0.00418, 0.09219)
  [../]
  [./top_temp_clad_outer]
    type = NodalVariableValue
    variable = temperature
    nodeid = 68 #mesh dependent, at (0.00474, 0.09219)
  [../]
  [./center_temp_clad_inner]
    type = NodalVariableValue
    variable = temperature
    nodeid = 36 #mesh dependent, at (0.00418, 0.05661)
  [../]
  [./center_temp_clad_outer]
    type = NodalVariableValue
    variable = temperature
    nodeid = 43 #mesh dependent, at (0.00474, 0.05661)
  [../]
  [./bottom_temp_clad_inner]
    type = NodalVariableValue
    variable = temperature
    nodeid = 9 #mesh dependent, at (0.00418, 0.02103)
  [../]
  [./bottom_temp_clad_outer]
    type = NodalVariableValue
    variable = temperature
    nodeid = 16 #mesh dependent, at (0.00418, 0.02103)
  [../]

  ### Nodal penetration
  [./top_penetration_fuel]
    type = NodalVariableValue
    variable = penetration
    nodeid = 361 # mesh dependent, at (0.0041, 0.09219)
  [../]
  [./center_penetration_fuel]
    type = NodalVariableValue
    variable = penetration
    nodeid = 262 # mesh dependent, at (0.0041, 0.05661)
  [../]
  [./bottom_penetration_fuel]
    type = NodalVariableValue
    variable = penetration
    nodeid = 163 # mesh dependent, at (0.0041, 0.02103)
  [../]

  ### Nodal contact pressure
  [./top_contact_pressure_fuel]
    type = NodalVariableValue
    variable = contact_pressure
    nodeid = 361 # mesh dependent, at (0.0041, 0.09219)
  [../]
  [./center_contact_pressure_fuel]
    type = NodalVariableValue
    variable = contact_pressure
    nodeid = 262 # mesh dependent, at (0.0041, 0.05661)
  [../]
  [./bottom_contact_pressure_fuel]
    type = NodalVariableValue
    variable = contact_pressure
    nodeid = 163 # mesh dependent, at (0.0041, 0.02103)
  [../]
### End of 1.5D comparisons

  [./center_effective_creep_rate_inner]
    type = ElementalVariableValue
    elementid = 16 # mesh dependent
    variable = creep_strain_rate
  [../]
  [./center_effective_creep_rate_outer]
    type = ElementalVariableValue
    elementid = 19 # mesh dependent
    variable = creep_strain_rate
  [../]
  [./effective_creep_strain_rate]
    type = ElementAverageValue
    variable = creep_strain_rate
    block = clad
  [../]

  [./solid_swelling]
    type = ElementAverageValue
    variable = solid_swell
    block = fuel
  [../]
  [./gas_swelling]
    type = ElementAverageValue
    variable = gas_swell
    block = fuel
  [../]
  [./densification]
    type = ElementAverageValue
    variable = densification
    block = fuel
  [../]
  [./volumetric_swelling]
    type = ElementAverageValue
    variable = volumetric_swelling_strain
    block = fuel
  [../]
  [./relocation]
    type = ElementAverageValue
    variable = relocation
    block = fuel
  [../]
[]

[VectorPostprocessors]
  [./clad]
    type = NodalValueSampler
    variable = disp_x
    boundary = 2
    sort_by = y
    outputs = 'clad_radial_displacement'
  [../]
  [./pellet]
    type = NodalValueSampler
    variable = disp_x
    boundary = 10
    sort_by = y
    outputs = 'fuel_radial_displacement'
  [../]
[]

[Outputs]
  exodus = true
  csv = true
  color = false
  print_perf_log = true
  [./clad_radial_displacement]
    type = CSV
    execute_on = 'FINAL'
  [../]
  [./fuel_radial_displacement]
    type = CSV
    execute_on = 'FINAL'
  [../]
[]
(examples/tensor_mechanics/1.5D_rodlet_10pellets/Smeared_1.5D.i)