# 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.


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

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)