Axial Relocation Action System

Description

The AxialRelocation action simplifies the input file of a LWR simulation that includes axial relocation by creating the necessary AuxVariables, AuxKernels, Materials, UserObjects, and Postprocessors required for the Axial Relocation Model.

warning:For LWR Use Only

This action is designed for use with only Light Water Reactor simulations run with a mesh generated by Layered1DMeshGenerator or Layered2DMeshGenerator.

Constructed Objects

The objects constructed by the AxialRelocation action are summarized in Table 1. The Generated Name column refers to the name of the object constructed by the action. Depending on the functionality of the object, the name could represent an aux variable, material, postprocessor or userobject.

Table 1: Correspondence Among Action Functionality and Moose/BISON Objects for the AxialRelocation Action

Functionality	Generated Name	Created Classes	Associated Parameters
Inner clad radius aux variable	`inner_clad_radius`	Radius
Outer fuel radius aux variable	`outer_fuel_radius`	Radius
Pulverized aux variable	`pulverized`	MaterialRealAux
Axial relocation strain aux variable	`axial_relocation_strain`	MaterialRealAux
Gap conductivity aux variable	`gap_conductivity`	MaterialRealAux
Axial relocation output aux variables	`layered_` + settings of `axial_relocation_output_options` parameter	AxialRelocationOutputAux	`axial_relocation_output_options`: the properties to output from the `axial_relocation` userobject
UserObject output aux variables	`layered_average_contact_pressure`	SpatialUserObjectAux	`output_userobjects`: whether or not to create aux variables for the properties calculated in the user objects used in the axial relocation calculation
	`layered_pulverized_fuel_volume`	SpatialUserObjectAux
	`layered_average_burnup`	SpatialUserObjectAux
	`layered_maximum_clad_radius`	SpatialUserObjectAux
	`layered_average_gap_conductivity`	SpatialUserObjectAux
UO $var element = document.getElementById("moose-equation-d047ca05-5e08-4d37-9ce2-dc96023b68e4");katex.render("_2", element, {displayMode:false,throwOnError:false});$ pulverization material model		UO2Pulverization or UO2PulverizationMesoscale	`burnup_variable`: the burnup variable
			`temperature_variable`: the temperature variable
UO $var element = document.getElementById("moose-equation-1504258d-7a63-43f8-981b-969ab058db50");katex.render("_2", element, {displayMode:false,throwOnError:false});$ axial relocation eigenstrain material model		UO2AxialRelocationEigenstrain
Linear heat rate postprocessor	`rod_input_power`	FunctionValuePostprocessor	`rod_ave_lin_pow`: function containing the linear heat generation rate
Maximum linear heat rate as a function of time postprocessor	`maximum_power`	TimeExtremeValue
Clad internal volume postprocessor	`clad_volume`	LayeredInternalVolumePostprocessor	`out_of_plane_strain_variable`: the out-of-plane strain variable
			`clad_inner_volume_addition`: addition to the internal clad volume. Typically used for refabricated rods.
Fuel volume postprocessor	`pellet_volume`	LayeredCrumbledInternalVolumePostprocessor	`out_of_plane_strain_variable`: the out-of-plane strain variable
Gas volume postprocessor	`gas_volume`	LinearCombinationPostprocessor
Fuel pin geometry user object	`fuel_pin_geometry`	Layered1DFuelPinGeometry or Layered2DFuelPinGeometry	`formulation`: the type of layered mesh being used in the analysis
Layered average contact pressure user object	`layered_average_contact_pressure`	LayeredSideAverage	`contact_pressure_variable`: the contact pressure variable (default name of `contact_pressure` comes from the contact modules)
			`axial_direction`: the axial direction of the layers
			`num_layers`: the number of layers
Layered axial relocation eigenstrain user object	`layered_eigenstrain`	LayeredAxialRelocationEigenstrainUserObject	`penetration_variable`: the penetration variable (default name of `penetration` comes from the contact module)
			`axial_direction`: the axial direction of the layers
			`num_layers`: the number of layers
Layered pulverized volume user object	`layered_pulverized_fuel_volume`	LayeredVariableIntegral	`axial_direction`: the axial direction of the layers
			`num_layers`: the number of layers
Layered average burnup user object	`layered_average_burnup`	LayeredAverage	`axial_direction`: the axial direction of the layers
			`num_layers`: the number of layers
Layered maximum inner cladding radius user object	`layered_maximum_clad_radius`	LayeredNodalExtremeValue	`axial_direction`: the axial direction of the layers
			`num_layers`: the number of layers
Layered clad internal volume user object	`layered_clad_internal_volume`	LayeredInternalVolume	`out_of_plane_strain_variable`: the out-of-plane strain variable
			`axial_direction`: the axial direction of the layers
			`num_layers`: the number of layers
Layered average gap conductivity user object	`layered_average_gap_conductivity`	LayeredNodalExtremeValue	`axial_direction`: the axial direction of the layers
			`num_layers`: the number of layers
Axial relocation user object	`axial_relocation`	AxialRelocationUserObject	`axial_direction`: the axial direction of the layers
			`fragmentation_model`: the fragmentation model to use
			`num_layers`: the number of layers
			`density`: as-fabricated fuel density
			`gap_thickness_threshold`: the threshold gap thickness required for axial relocation
			`nonrelocatable_fuel_fraction`: the fraction of the initial fuel in a layer that cannot relocate
			`fragment_packing_fraction`: the packing fraction of large fragments of fuel
			`pulver_packing_fraction`: the packing fraction of the smaller pulvers of fuel
			`pulver_characteristic_length`: the characteristic length of the pulvers
			`tolerance`: the tolerance used in the sub-newton iteration loop for determining the packing fraction
			\| `maximum_iterations`: the maximum number of sub-newton iterations allowed when determining the packing fraction
			`output_iteration_info`: whether or not to output sub-newton iteration information

If the full thermo-mechanical coupling during axial relocation is desired, the axial_relocation_object parameter must be set equal to axial_relocation in NeutronHeatSource, UO2Thermal, and UO2ElasticityTensor. The gap_thermal_conductivity parameter in UO2Thermal must also be set to layered_average_gap_conductivity.

Figure 1 provides an illustration on how all of the different components (AuxKernels, Materials, and UserObjects) of the axial relocation system created by the AxialRelocationAction interact to predict the redistribution of fuel during Loss of Coolant Accidents.

Figure 1: A flowchart that illustrates how the different components created by the AxialRelocationAction interact with one another.

Example Input File Syntax

An example of using the AxialRelocationAction is given by:

[AxialRelocation<<<{"href": "index.html"}>>>]
  [relocation]
    rod_ave_lin_pow<<<{"description": "Rod average linear power function."}>>> = power
    axial_direction<<<{"description": "The axial direction of the model."}>>> = y
    fuel_blocks<<<{"description": "The list of names of the blocks (subdomains) for the fuel pellets."}>>> = fuel
    clad_blocks<<<{"description": "The list of names of the block (subdomain) for the cladding."}>>> = clad
    contact_pressure_variable<<<{"description": "The contact pressure variable."}>>> = contact_pressure
    out_of_plane_strain_variable<<<{"description": "The out-of-plane strain nodal variable."}>>> = strain_yy_0
    penetration_variable<<<{"description": "The penetration variable."}>>> = penetration
    clad_inner_volume_addition<<<{"description": "An additional volume to be included in the internal volume calculation. A time-dependent function is expected."}>>> = 0
    burnup_variable<<<{"description": "The burnup variable."}>>> = burnup
    temperature<<<{"description": "Coupled temperature"}>>> = temperature
    axial_relocation_output_options<<<{"description": "The layered quantities to be output from the axial relocation userobject."}>>> = 'MASS_FRACTION PACKING_FRACTION'
    mesh_generator<<<{"description": "The name of the generator to use as the prefix for mesh meta data properties."}>>> = layered1D_mesh
  []
[]

An example of NeutronHeatSource when using full thermo-mechanical coupling is given by:

[Kernels<<<{"href": "../Kernels/index.html"}>>>]
  [heat_source]
    type = NeutronHeatSource<<<{"description": "Compute heat generation due to fission.", "href": "../../source/kernels/NeutronHeatSource.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = fuel
    burnup_function<<<{"description": "Burnup function"}>>> = burnup
    axial_relocation_object<<<{"description": "Name of the AxialRelocationUserObject that determines whether the fuel has crumbled."}>>> = axial_relocation
    extra_vector_tags<<<{"description": "The extra tags for the vectors this Kernel should fill"}>>> = 'ref'
  []
[]

An example of UO2Thermal when using full thermo-mechanical coupling is given by:

[Materials<<<{"href": "../Materials/index.html"}>>>]
  [fuel_thermal]
    type = UO2Thermal<<<{"description": "Computes thermal conductivity and specific heat capacity of uranium dioxide fuel.", "href": "../../source/materials/UO2Thermal.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = fuel
    thermal_conductivity_model<<<{"description": "The thermal conductivity model."}>>> = STAICU
    hbs_porosity_correction<<<{"description": "The porosity correction method for the thermal conductivity model in the high burnup structure."}>>> = KAMPF
    model_hbs_formation<<<{"description": "Flag on whether to take into account the volume fraction of high burnup structure that has formed on the thermal conductivity calculation."}>>> = true
    temperature<<<{"description": "Coupled temperature"}>>> = temperature
    burnup_function<<<{"description": "Burnup function"}>>> = burnup
    axial_relocation_object<<<{"description": "Name of the AxialRelocationUserObject that determines whether the fuel has crumbled."}>>> = axial_relocation
    gap_thermal_conductivity<<<{"description": "The layered average thermal conductivity across the gas gap."}>>> = layered_average_gap_conductivity
  []
[]

An example of UO2ElasticityTensor when using full thermo-mechanical coupling is given by:

[Materials<<<{"href": "../Materials/index.html"}>>>]
  [fuel_elasticity_tensor]
    type = UO2IsotropicDamageElasticityTensor<<<{"description": "Computes the isotropic elastic constants for UO2 fuel as a scaled function of the number of cracks in the fuel.", "href": "../../source/materials/solid_mechanics/UO2IsotropicDamageElasticityTensor.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = fuel
    fragmentation_model<<<{"description": "The model used to calculate the number of cracks in the fuel."}>>> = BARANI
    temperature<<<{"description": "Coupled temperature"}>>> = temperature
    rod_ave_lin_pow<<<{"description": "Average linear power function"}>>> = power_history
    axial_relocation_object<<<{"description": "Name of the AxialRelocationUserObject that determines whether the fuel has crumbled."}>>> = axial_relocation
    crumbling_scale_factor<<<{"description": "Scaling factor to apply to the Young's modulus in layers that have crumbled during axial relocation."}>>> = 0.0001
  []
[]

Input Parameters

axial_directionThe axial direction of the model.
C++ Type:MooseEnum
Options:x, y, z
Controllable:No
Description:The axial direction of the model.
axial_relocation_output_optionsThe layered quantities to be output from the axial relocation userobject.
C++ Type:MultiMooseEnum
Options:MASS_FRACTION, MASS, PULVER_MASS_FRACTION, PACKING_FRACTION, RELOCATABLE_MASS, CRUMBLED_MASS, LOWER_MASS_CONSTRAINT, UPPER_MASS_CONSTRAINT, CRUMBLED_STATUS
Controllable:No
Description:The layered quantities to be output from the axial relocation userobject.
contact_pressure_variableThe contact pressure variable.
C++ Type:VariableName
Unit:(no unit assumed)
Controllable:No
Description:The contact pressure variable.
fragmentation_modelCOINDREAUThe model used to calculate the number of large fragments.
Default:COINDREAU
C++ Type:MooseEnum
Options:COINDREAU, WALTON, BARANI
Controllable:No
Description:The model used to calculate the number of large fragments.
mesh_generatorThe name of the generator to use as the prefix for mesh meta data properties.
C++ Type:MeshGeneratorName
Controllable:No
Description:The name of the generator to use as the prefix for mesh meta data properties.
out_of_plane_strain_variableThe out-of-plane strain nodal variable.
C++ Type:VariableName
Unit:(no unit assumed)
Controllable:No
Description:The out-of-plane strain nodal variable.
penetration_variableThe penetration variable.
C++ Type:VariableName
Unit:(no unit assumed)
Controllable:No
Description:The penetration variable.
rod_ave_lin_powRod average linear power function.
C++ Type:FunctionName
Unit:(no unit assumed)
Controllable:No
Description:Rod average linear power function.
temperatureCoupled temperature
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Coupled temperature

Required Parameters

active__all__ If specified only the blocks named will be visited and made active
Default:__all__
C++ Type:std::vector<std::string>
Controllable:No
Description:If specified only the blocks named will be visited and made active
burnup_functionBurnup function
C++ Type:FunctionName
Unit:(no unit assumed)
Controllable:No
Description:Burnup function
burnup_variableThe burnup variable.
C++ Type:VariableName
Unit:(no unit assumed)
Controllable:No
Description:The burnup variable.
clad_blocksThe list of names of the block (subdomain) for the cladding.
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of names of the block (subdomain) for the cladding.
clad_inner_volume_addition0An additional volume to be included in the internal volume calculation. A time-dependent function is expected.
Default:0
C++ Type:FunctionName
Unit:(no unit assumed)
Controllable:No
Description:An additional volume to be included in the internal volume calculation. A time-dependent function is expected.
critical_gb_stress1.3e+08Critical stress at which a grain boundary fractures for the pulverization model
Default:1.3e+08
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Critical stress at which a grain boundary fractures for the pulverization model
density10431As-fabricated fuel pellet density.
Default:10431
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:As-fabricated fuel pellet density.
evolve_bubble_pressure_hbsTrueAllow bubble pressure to evolve from initial conditions based on evolution equations.
Default:True
C++ Type:bool
Controllable:No
Description:Allow bubble pressure to evolve from initial conditions based on evolution equations.
fission_rateCoupled fission rate
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Coupled fission rate
formulationLayered1DThe layered formulation being used in the simulation.
Default:Layered1D
C++ Type:MooseEnum
Options:Layered1D, Layered2D
Controllable:No
Description:The layered formulation being used in the simulation.
fragment_packing_fraction0.69The packing fraction of large fragments of fuel.
Default:0.69
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The packing fraction of large fragments of fuel.
fragment_particle_shapeTriangular_PrismThe larger particle shape used to calculate the pulverized binary packing fraction: Sphere Cube Octahedron Ideal_Cylinder Triangular_Prism
Default:Triangular_Prism
C++ Type:MooseEnum
Options:Sphere, Cube, Octahedron, Ideal_Cylinder, Triangular_Prism
Controllable:No
Description:The larger particle shape used to calculate the pulverized binary packing fraction: Sphere Cube Octahedron Ideal_Cylinder Triangular_Prism
fuel_blocksThe list of names of the blocks (subdomains) for the fuel pellets.
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of names of the blocks (subdomains) for the fuel pellets.
gap_thickness_threshold0.0002The threshold gap thickness required for axial relocation.
Default:0.0002
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The threshold gap thickness required for axial relocation.
hbs_materialhbs_formationName of the HighBurnupStructureFormation block.
Default:hbs_formation
C++ Type:MaterialName
Controllable:No
Description:Name of the HighBurnupStructureFormation block.
id_offset0This offset is added to the generated boundary IDs
Default:0
C++ Type:short
Controllable:No
Description:This offset is added to the generated boundary IDs
inactiveIf specified blocks matching these identifiers will be skipped.
C++ Type:std::vector<std::string>
Controllable:No
Description:If specified blocks matching these identifiers will be skipped.
maximum_iterations10Maximum number of sub-newton iterations allowed when calculating the packing fraction.
Default:10
C++ Type:unsigned int
Controllable:No
Description:Maximum number of sub-newton iterations allowed when calculating the packing fraction.
name_prefixIf provided, prefix the built in boundary names with this string
C++ Type:std::string
Controllable:No
Description:If provided, prefix the built in boundary names with this string
nonrelocatable_fuel_fraction0.01The fraction of the initial fuel in a layer that cannot relocate.
Default:0.01
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The fraction of the initial fuel in a layer that cannot relocate.
output_iteration_infoFalseSet true to output sub-newton iteration information when calculating the packing fraction.
Default:False
C++ Type:bool
Controllable:No
Description:Set true to output sub-newton iteration information when calculating the packing fraction.
output_userobjectsTrueWhether or not to create AuxVariables and AuxKernels to output the results of the userobjects used in the axial relocation calculation.
Default:True
C++ Type:bool
Controllable:No
Description:Whether or not to create AuxVariables and AuxKernels to output the results of the userobjects used in the axial relocation calculation.
pulver_characteristic_length0.0001The characteristic length of the pulvers.
Default:0.0001
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The characteristic length of the pulvers.
pulver_packing_fraction0.72The packing fraction of the smaller pulvers of fuel.
Default:0.72
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The packing fraction of the smaller pulvers of fuel.
pulver_particle_shapeOctahedronThe smaller particle shape used to calculate the pulverized binary packing fraction: Sphere Cube Octahedron Ideal_Cylinder Triangular_Prism
Default:Octahedron
C++ Type:MooseEnum
Options:Sphere, Cube, Octahedron, Ideal_Cylinder, Triangular_Prism
Controllable:No
Description:The smaller particle shape used to calculate the pulverized binary packing fraction: Sphere Cube Octahedron Ideal_Cylinder Triangular_Prism
pulverization_thresholdTURNBULLThe model to use for the pulverization threshold.
Default:TURNBULL
C++ Type:MooseEnum
Options:TURNBULL, ANALYTICAL, PHASE_FIELD_2D, PHASE_FIELD_3D, PHASE_FIELD_2D_MULTISCALE
Controllable:No
Description:The model to use for the pulverization threshold.
residual_gap_size2e-06The residual gap remaining in layers that have crumbled.
Default:2e-06
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The residual gap remaining in layers that have crumbled.
scalar_critical_pressure1Multiplicative scalar to adjust the value of the critical fragmentation pressure in the pulverization model
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Multiplicative scalar to adjust the value of the critical fragmentation pressure in the pulverization model
tolerance1e-10Absolute convergence tolerance for the sub-newton iteration used to calculate the packing fraction.
Default:1e-10
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Absolute convergence tolerance for the sub-newton iteration used to calculate the packing fraction.
translation_vector0 0 0The value to use for the translation. The xyz coordinates are applied in each direction respectively.
Default:0 0 0
C++ Type:libMesh::VectorValue<double>
Unit:(no unit assumed)
Controllable:No
Description:The value to use for the translation. The xyz coordinates are applied in each direction respectively.
use_automatic_differentiationFalseFlag to use automatic differentiation (AD) objects when possible
Default:False
C++ Type:bool
Controllable:No
Description:Flag to use automatic differentiation (AD) objects when possible
use_axial_gas_communicationFalseUsed when the axial gas communication model is run with relocation.
Default:False
C++ Type:bool
Controllable:No
Description:Used when the axial gas communication model is run with relocation.

Optional Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.

Advanced Parameters

Associated Actions

Available Actions

Bison App
AxialRelocationActionSets up AuxVariables, AuxKernels, Materials, Postprocessors, and Userobjects for simulations containing axial fuel relocation.

(test/tests/axial_relocation/axial_relocation_eigenstrain_action.i)

# The purpose of this test is to verify the calculation of the eigenstrain
# applied to the fuel in the event that an axial layer has crumbled (accommodated
# additional mass) during axial relocation in Layered 1D. This eigenstrain is used
# to move the outer radius of the fuel to within a small residual gap size of the
# inside surface of the cladding.
#
# The eigenstrain is calculated as:
#
# epsilon = ln(1.0 + (R_fcurr + g_width - R_fo - g^r) / R_fo)
#
# where R_fcurr is the current outer radius of the pellet in the layer,
# g_width is the current fuel-to-clad gap width, g^r is the residual gap size
# (default of 2.0e-6 m), and R_fo is the fuel radius at the time the fuel begins
# to move to that layer.  In this test case since there is no thermal expansion
# or other radial displacements of the fuel, R_fo is the as-fabricated fuel radius.
#
# In this test crumbling of the middle layer of 5 (layer 2 since indexing of
# the layers begins at zero) occurs at the end of the simulation. Since only this
# layer is deemed crumbled it is only this layer to which the eigenstrain is applied.
# At this time the inner cladding radius is calculated to be 4.5e-3 + 0.00052 =
# 5.02e-3 m. Therefore, the axial relocation eigenstrain is calculated to be:
#
# epsilon = ln(1.0 + (5.02e-3 - 4.5e-3 - 2.0e-6) / 4.5e-3)
#
# epsilon = 0.108954
#
# This is verified by checking the value of the axial_relocation_strain column in the
# outputs vector postprocessor. The axial variation of packing fraction is also output
# as this input is used with cli_args to test different particle shape options and the
# impact on the packing fraction.

[GlobalParams]
  order = FIRST
  family = LAGRANGE
  displacements = disp_x
  temperature = temperature
[]

[Mesh]
  coord_type = RZ
  [layered1D_mesh]
    type = Layered1DMeshGenerator
    elem_type = EDGE2
    slices_per_block = 5
    pellet_outer_radius = 4.5e-3
    fuel_height = 0.5
    include_plenum = false
    nx_p = 10
    clad_gap_width = 0.0
    pellet_mesh_density = customize
    pellet_bottom_coor = 0.0
  []
  patch_update_strategy = auto
[]

[Variables]
  [temperature]
    order = FIRST
    family = LAGRANGE
    initial_condition = 1200
  []
[]

[AuxVariables]
  [disp_y]
  []
  [disp_z]
  []
  [burnup]
    order = FIRST
    family = LAGRANGE
  []
  [strain_yy_0]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [burnup_function]
    type = ParsedFunction
    expression = 'x*100.0/4.275'
  []
  [power]
    type = PiecewiseLinear
    x = '0 100'
    y = '15000 15000'
  []
  [clad_displacement_function]
    type = ParsedFunction
    expression = '2.0e-5 * t * sin(pi * y / 0.5)'
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temperature
  []
[]

[Physics]
  [SolidMechanics]
    [Layered1D]
      [gps_fuel]
        add_scalar_variables = true
        add_variables = true
        out_of_plane_strain_name = strain_yy
        fuel_pin_geometry = fuel_pin_geometry
        strain = finite
        block = fuel
        eigenstrain_names = 'axial_relocation_eigenstrain'
        decomposition_method = EigenSolution
        mesh_generator = layered1D_mesh
      []
      [gps_clad]
        add_scalar_variables = true
        add_variables = true
        out_of_plane_strain_name = strain_yy
        fuel_pin_geometry = fuel_pin_geometry
        strain = finite
        block = clad
        decomposition_method = EigenSolution
        mesh_generator = layered1D_mesh
      []
    []
  []
[]

[AuxKernels]
  [burnup]
    type = FunctionAux
    variable = burnup
    function = burnup_function
    execute_on = 'initial linear'
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    penalty = 1e7
    formulation = kinematic
    model = frictionless
  []
[]

[ThermalContact]
  [thermal_contact]
    type = GasGapHeatTransfer
    variable = temperature
    primary = 5
    secondary = 10
    gap_conductivity = 1
  []
[]

[BCs]
  [temperature]
    type = DirichletBC
    boundary = '10 12 5 2'
    variable = temperature
    value = 1200
  []
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 12
    value = 0.0
  []
  [inner_clad_displacement]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = '5'
    function = clad_displacement_function
  []
[]

[Materials]
  [fuel_thermal]
    type = HeatConductionMaterial
    block = fuel
    thermal_conductivity = 1.0
    specific_heat = 1.0
  []
  [fuel_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = fuel
    youngs_modulus = 2.e11
    poissons_ratio = .345
  []
  [fuel_stress]
    type = ComputeFiniteStrainElasticStress
    block = fuel
  []
  [clad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = clad
    youngs_modulus = 7.5e10
    poissons_ratio = 0.3
  []
  [clad_stress]
    type = ComputeFiniteStrainElasticStress
    block = clad
  []
  [clad_thermal]
    type = HeatConductionMaterial
    block = clad
    thermal_conductivity = 16.0
    specific_heat = 330.0
  []
[]

[AxialRelocation]
  [relocation]
    rod_ave_lin_pow = power
    axial_direction = y
    fuel_blocks = fuel
    clad_blocks = clad
    contact_pressure_variable = contact_pressure
    out_of_plane_strain_variable = strain_yy_0
    penetration_variable = penetration
    clad_inner_volume_addition = 0
    burnup_variable = burnup
    temperature = temperature
    axial_relocation_output_options = 'MASS_FRACTION PACKING_FRACTION'
    mesh_generator = layered1D_mesh
  []
[]

[Executioner]
  type = Transient

  solve_type = PJFNK

  petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-ksp_gmres_restart'
  petsc_options_value = '101'

  line_search = 'none'

  l_max_its = 100
  nl_max_its = 100
  nl_rel_tol = 1e-10
  nl_abs_tol = 1e-10
  l_tol = 1e-5
  start_time = 0.0
  num_steps = 14
  dt = 2
[]

[VectorPostprocessors]
  [vpp]
    type = LineValueSampler
    start_point = '2.25e-3 0.05 0'
    end_point = '2.25e-3 0.45 0'
    num_points = 5
    sort_by = y
    variable = 'axial_relocation_strain layered_packing_fraction'
    outputs = results
  []
[]

[Outputs]
  hide = penetration
  [results]
    type = CSV
    execute_on = final
    create_final_symlink = true
  []
[]

(assessment/LWR/validation/LOCA_IFA_650/analysis/IFA_650_9/IFA_650_9_part1.i)


initial_fuel_density = 10430.0

[GlobalParams]
  density = ${initial_fuel_density}
  initial_porosity = 0.048
  order = SECOND
  family = LAGRANGE
  displacements = disp_x
  temperature = temperature
  energy_per_fission = 3.2e-11 #J/fission
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
  acceptable_multiplier = 10
[]

[Mesh]
  coord_type = RZ
  [layered1D_mesh]
    type = Layered1DMeshGenerator
    slices_per_block = 30
    pellet_outer_radius = 4.565e-3
    clad_gap_width = 0.085e-3
    clad_thickness = 0.725e-3
    fuel_height = 0.480
    plenum_height = 0.262416
    pellet_mesh_density = customize
    clad_mesh_density = customize
    nx_p = 11
    nx_c = 5
  []
  patch_update_strategy = auto
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[Variables]
  [disp_x]
  []
  [temperature]
    initial_condition = 295.0
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    data_file = power_history.csv
    scale_factor = 1.0
    format = columns
  []
  [axial_peaking_factors]
    type = PiecewiseBilinear
    data_file = axial_peaking_factors.csv
    axis = 1
    scale_factor = 1
  []
  [pressure_ramp]
    type = PiecewiseLinear
    data_file = coolant_pressure.csv
    scale_factor = 1
    format = columns
  []
  [average_htc]
    type = PiecewiseBilinear
    data_file = average_coolant_htc.csv
    axis = 1
    scale_factor = 1
  []
  [forced_times]
    type = PiecewiseLinear
    data_file = timestep_limiting.csv
    scale_factor = 1
    format = columns
  []
  [heat_sink_temperature]
    type = PiecewiseBilinear
    data_file = heater_temp.csv
    scale_factor = 1
    axis = 1
  []
  [clad_outer_temperature]
    type = PiecewiseBilinear
    data_file = clad_surface_temp.csv
    scale_factor = 1
    axis = 1
  []
  [heat_transfer_mode]
    type = PiecewiseConstant
    x = '-200  200412461 200413048'
    y = '9  9          8       '
    direction = 'right'
  []
  [clad_axial_pressure]
    type = CladdingAxialPressureFunction
    plenum_pressure = plenum_pressure
    coolant_pressure = pressure_ramp
    coolant_pressure_scaling_factor = 1.0
    fuel_pin_geometry = fuel_pin_geometry
  []
  [fuel_axial_pressure]
    type = ParsedFunction
    expression = plenum_pressure
    symbol_names = plenum_pressure
    symbol_values = plenum_pressure
  []
[]

[AuxVariables]
  [disp_y]
  []
  [disp_z]
  []
  [fast_neutron_flux]
    block = clad
  []
  [fast_neutron_fluence]
    block = clad
  []
  [grain_radius]
    block = fuel
    initial_condition = 5.0e-6
  []
  [effective_creep_strain]
    block = clad
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxide_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
  [gap_conductance]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_rate]
    order = CONSTANT
    family = MONOMIAL
  []
  [layered_maximum_clad_radius]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_source]
    type = NeutronHeatSource
    variable = temperature
    block = fuel
    burnup_function = burnup
    axial_relocation_object = axial_relocation
    extra_vector_tags = 'ref'
  []
[]

[Physics]
  [SolidMechanics]
    [Layered1D]
      [fuel]
        add_scalar_variables = true
        out_of_plane_strain_name = strain_yy
        fuel_pin_geometry = fuel_pin_geometry
        out_of_plane_pressure_function = fuel_axial_pressure
        strain = finite
        block = fuel
        eigenstrain_names = 'fuel_thermal_strain fuel_swelling_strain
          fuel_relocation_strain axial_relocation_eigenstrain'
        decomposition_method = EigenSolution
        cylindrical_axis_point1 = '0 0 0'
        cylindrical_axis_point2 = '0 1 0'
        generate_output = 'stress_xx stress_yy stress_zz vonmises_stress
          hoop_stress'
        extra_vector_tags = 'ref'
        group_scalar_vars_in_reference_residual = true
        mesh_generator = layered1D_mesh
      []
      [clad]
        add_scalar_variables = true
        out_of_plane_strain_name = strain_yy
        fuel_pin_geometry = fuel_pin_geometry
        strain = finite
        out_of_plane_pressure_function = clad_axial_pressure
        block = clad
        cylindrical_axis_point1 = '0 0 0'
        cylindrical_axis_point2 = '0 1 0'
        eigenstrain_names = 'clad_thermal_strain clad_irradiation_strain'
        generate_output = 'stress_xx stress_yy stress_zz vonmises_stress
          strain_zz creep_strain_zz hoop_stress'
        decomposition_method = EigenSolution
        extra_vector_tags = 'ref'
        group_scalar_vars_in_reference_residual = true
        mesh_generator = layered1D_mesh
      []
    []
  []
[]

[Burnup]
  [burnup]
    block = fuel
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    num_radial = 80
    num_axial = 11
    fuel_pin_geometry = fuel_pin_geometry
    fuel_volume_ratio = 1.0
    order = CONSTANT
    family = MONOMIAL
    RPF = RPF
    isotopes = 'U235 U238 Pu239 Pu240 Pu241 Pu242'
    isotope_fractions = '0.035 0.965 0 0 0 0'
  []
[]

[AuxKernels]
  [fast_neutron_flux]
    type = FastNeutronFluxAux
    block = clad
    variable = fast_neutron_flux
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    factor = 3e13
    execute_on = timestep_begin
  []
  [fast_neutron_fluence]
    type = FastNeutronFluenceAux
    block = clad
    variable = fast_neutron_fluence
    fast_neutron_flux = fast_neutron_flux
    execute_on = timestep_begin
  []
  [grain_radius]
    type = GrainRadiusAux
    block = fuel
    variable = grain_radius
    temperature = temperature
    execute_on = linear
  []
  [effective_creep_strain]
    type = MaterialRealAux
    block = clad
    variable = effective_creep_strain
    property = effective_creep_strain
    execute_on = 'timestep_end'
  []
  [fract_bphase]
    type = MaterialRealAux
    block = clad
    variable = fract_beta_phase
    property = fract_beta_phase
    execute_on = 'initial linear'
  []
  [oxide_thickness]
    type = MaterialRealAux
    boundary = 2
    variable = oxide_thickness
    property = oxide_scale_thickness
    execute_on = 'initial linear'
  []
  [hasburst]
    type = MaterialRealAux
    boundary = 2
    variable = burst
    property = failed
    execute_on = timestep_end
  []
  [conductance]
    type = MaterialRealAux
    boundary = 10
    property = gap_conductance
    variable = gap_conductance
    execute_on = 'initial linear'
  []
  [coolant_htc]
    type = MaterialRealAux
    property = coolant_channel_htc
    variable = coolant_htc
    boundary = 2
    execute_on = 'initial linear'
  []
  [creep_rate]
    type = MaterialRealAux
    block = clad
    variable = creep_rate
    property = creep_rate
    execute_on = timestep_end
  []
[]

[AxialRelocation]
  [rel]
    rod_ave_lin_pow = power_history
    axial_direction = y
    fuel_blocks = fuel
    clad_blocks = clad
    contact_pressure_variable = contact_pressure
    out_of_plane_strain_variable = strain_yy
    penetration_variable = penetration
    clad_inner_volume_addition = 0
    burnup_variable = burnup
    temperature = temperature
    gap_thickness_threshold = 0.00039
    axial_relocation_output_options = 'MASS_FRACTION'
    mesh_generator = layered1D_mesh
  []
[]

[CoolantChannel]
  [convective_clad_surface] # apply convective boundary to clad outer surface
    boundary = 2
    variable = temperature
    heat_transfer_mode = heat_transfer_mode
    heat_transfer_coefficient = average_htc # Calculated from an initial simulation of the base irradiation using the inlet_pressure, inlet_massflux, and inlet_temperature commented out below.
    inlet_temperature = heat_sink_temperature # K
    effective_emissivity = 0.75
    # inlet_temperature = 580
    # inlet_pressure    = 15.3e6   # Pa
    # inlet_massflux    = 3800     # kg/m^2-sec
    rod_diameter = 0.01075 # m
    rod_pitch = 1.26e-2 # m
    compute_enthalpy = false
    linear_heat_rate = power_history
    axial_power_profile = axial_peaking_factors
    output_properties = 'coolant_channel_htype coolant_channel_hmode'
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    penalty = 1e7
    formulation = kinematic
    model = frictionless
  []
[]

[ThermalContact]
  [thermal_contact]
    type = GasGapHeatTransfer
    variable = temperature
    primary = 5
    secondary = 10
    initial_moles = initial_moles
    gas_released = fis_gas_released
    plenum_pressure = plenum_pressure
    contact_pressure = contact_pressure
    jump_distance_model = LANNING
    roughness_coef = 3.2
    refab_gas_types = 'He Ar'
    refab_fractions = '0.05 0.95'
    refab_time = 199159200
    refab_type = 0
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 12
    value = 0.0
  []
  [Pressure]
    [coolantPressure]
      boundary = 2
      function = pressure_ramp
      factor = 1.0
    []
  []
  [clad_outer_temp]
    type = FunctionDirichletBC
    boundary = 2
    variable = temperature
    function = clad_outer_temperature
  []
  [PlenumPressure]
    [plenumPressure]
      boundary = 9
      initial_pressure = 2.0e6
      startup_time = 0
      R = 8.3143
      output_initial_moles = initial_moles
      temperature = plenum_temp
      volume = 'clad_volume pellet_volume'
      material_input = fis_gas_released
      output = plenum_pressure
      refab_time = 199159200
      refab_pressure = 4.0e6
      refab_temperature = 295.0
      refab_volume = 1.9e-05
    []
  []
[]

[LayeredPlenumTemperature]
  [plenum_temp]
    boundary = 5
    fuel_pin_geometry = fuel_pin_geometry
    out_of_plane_strain = strain_yy
    inner_surfaces = '5'
    outer_surfaces = '10'
    temperature = temperature
  []
[]

[Controls]
  [period0]
    type = TimePeriod
    disable_objects = 'BCs/clad_outer_temp'
    start_time = -200.0
    end_time = 199159200.0
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
    execute_on = timestep_end
  []
[]

[Materials]
  [fuel_thermal]
    type = UO2Thermal
    block = fuel
    thermal_conductivity_model = STAICU
    hbs_porosity_correction = KAMPF
    model_hbs_formation = true
    temperature = temperature
    burnup_function = burnup
    axial_relocation_object = axial_relocation
    gap_thermal_conductivity = layered_average_gap_conductivity
  []
  [relocation]
    type = UO2RelocationEigenstrain
    block = fuel
    burnup_function = burnup
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    fuel_pin_geometry = fuel_pin_geometry
    burnup_relocation_stop = 0.024
    relocation_activation1 = 5000.0
    relocation_model = ESCORE_modified
    eigenstrain_name = fuel_relocation_strain
  []
  [fuel_thermal_strain]
    type = ComputeThermalExpansionEigenstrain
    block = fuel
    thermal_expansion_coeff = 10.0e-6
    stress_free_temperature = 295.0
    eigenstrain_name = fuel_thermal_strain
  []
  [fuel_swelling]
    type = UO2VolumetricSwellingEigenstrain
    gas_swelling_model_type = SIFGRS
    block = fuel
    burnup_function = burnup
    initial_fuel_density = 10430.0
    eigenstrain_name = fuel_swelling_strain
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = fuel
    temperature = temperature
    burnup_function = burnup
    grain_radius = grain_radius
    transient_option = MICROCRACKING_BURNUP
    diff_coeff_option = TURNBULL_D1_D2
    gbs_model = true
  []
  [fuel_elasticity_tensor]
    type = UO2IsotropicDamageElasticityTensor
    block = fuel
    fragmentation_model = BARANI
    temperature = temperature
    rod_ave_lin_pow = power_history
    axial_relocation_object = axial_relocation
    crumbling_scale_factor = 0.0001
  []
  [fuel_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'fuel_creep'
    block = fuel
  []

  [fuel_creep]
    type = UO2CreepUpdate
    block = fuel
    temperature = temperature
    burnup_function = burnup
    initial_grain_radius = 5.0e-6
  []

  [HBS]
    type = HighBurnupStructureFormation
    block = fuel
    burnup_function = burnup
    temperature = temperature
    output_properties = 'hbs_volume_fraction'
    outputs = 'exodus'
  []

  [clad_elasticity_tensor]
    type = ZryElasticityTensor
    block = clad
  []
  [stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'zrycreep'
    block = clad
  []
  [zrycreep]
    type = ZryCreepLOCAUpdate
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    model_irradiation_creep = true
    model_primary_creep = true
    model_thermal_creep = true
    max_inelastic_increment = 5e-4
    zircaloy_material_type = stress_relief_annealed
    block = clad
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = clad
    stress_free_temperature = 295.0
    eigenstrain_name = clad_thermal_strain
  []
  [irradiation_swelling]
    type = ZryIrradiationGrowthEigenstrain
    block = clad
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = stress_relief_annealed
    eigenstrain_name = clad_irradiation_strain
  []
  [clad_phase]
    type = ZrPhase
    block = clad
    temperature = temperature
    numerical_method = 2
  []
  [clad_oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 4.65e-03
    clad_outer_radius = 5.375e-03
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = cathcart
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
  []
  [clad_thermal]
    type = ZryThermal
    block = clad
    temperature = temperature
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = fuel
    strain_free_density = ${initial_fuel_density}
  []
  [clad_density]
    type = StrainAdjustedDensity
    block = clad
    strain_free_density = 6551.0
  []
[]

[Postprocessors]
  [ave_temp_interior]
    type = LayeredSideAverageValuePostprocessor
    boundary = 9
    variable = temperature
    execute_on = 'initial linear'
    fuel_pin_geometry = fuel_pin_geometry
  []
  [pellet_volume_2]
    type = LayeredInternalVolumePostprocessor
    boundary = 8
    component = 0
    fuel_pin_geometry = fuel_pin_geometry
    out_of_plane_strain = strain_yy
    execute_on = 'initial linear'
  []
  [avg_clad_temp]
    type = LayeredSideAverageValuePostprocessor
    boundary = 7
    variable = temperature
    fuel_pin_geometry = fuel_pin_geometry
    execute_on = 'initial linear'
  []
  [fis_gas_produced]
    type = LayeredElementIntegralFisGasGeneratedSifgrsPostprocessor
    block = fuel
    fuel_pin_geometry = fuel_pin_geometry
  []
  [fis_gas_released]
    type = LayeredElementIntegralFisGasReleasedSifgrsPostprocessor
    block = fuel
    fuel_pin_geometry = fuel_pin_geometry
  []
  [fis_gas_grain]
    type = LayeredElementIntegralFisGasGrainSifgrsPostprocessor
    block = fuel
    outputs = exodus
    fuel_pin_geometry = fuel_pin_geometry
  []
  [fis_gas_boundary]
    type = LayeredElementIntegralFisGasBoundarySifgrsPostprocessor
    block = fuel
    outputs = exodus
    fuel_pin_geometry = fuel_pin_geometry
  []
  [fission_gas_release]
    type = FGRPercent
    fission_gas_released = fis_gas_released
    fission_gas_generated = fis_gas_produced
  []
  [average_coolant_htc]
    type = LayeredSideAverageValuePostprocessor
    boundary = 2
    variable = coolant_htc
    execute_on = 'initial linear'
    fuel_pin_geometry = fuel_pin_geometry
  []
  [average_burnup]
    type = RodAverageBurnup
    burnup_function = burnup
  []
  [temp_clad_max]
    type = NodalExtremeValue
    block = clad
    value_type = max
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [betaph_fract_max]
    type = ElementExtremeValue
    value_type = max
    variable = fract_beta_phase
    block = clad
    execute_on = 'initial timestep_end'
  []
  [burst]
    type = ElementExtremeValue
    value_type = max
    variable = burst
    block = clad
    execute_on = 'initial timestep_end'
  []
  [timestep_material]
    type = MaterialTimeStepPostprocessor
    block = clad
    execute_on = 'initial timestep_end'
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = clad
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geometry
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
[]

[Dampers]
  [limitT]
    type = BoundingValueNodalDamper
    variable = temperature
    max_value = 3200.0
    min_value = 0.0
  []
  [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 -ksp_gmres_restart'
  petsc_options_value = 'lu       superlu_dist                  51'

  line_search = 'none'

  l_max_its = 50
  l_tol = 1e-3
  nl_max_its = 30
  nl_rel_tol = 1e-4
  nl_abs_tol = 1e-8
  dtmax = 5e5
  dtmin = 1e-5
  start_time = -200.0
  end_time = 199159200 # End base irradiation
  #  end_time = 200412431 # Begin Blowdown
  # end_time = 200413048 # End

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 200
    timestep_limiting_postprocessor = timestep_material
    optimal_iterations = 20
    iteration_window = 4
    linear_iteration_ratio = 100
    timestep_limiting_function = forced_times
    force_step_every_function_point = true
    max_function_change = 2000
    time_t = '199159200 200312431 200411431 200412431 200412461 200413048'
    time_dt = '1.0e04    1.0e04    10.0        5.0         0.5        5.0'
  []
[]

[VectorPostprocessors]
  [clad_radial_disp]
    type = NodalValueSampler
    variable = disp_x
    boundary = 2
    sort_by = y
    outputs = 'outfile_1'
  []
  [clad_out_temp]
    type = NodalValueSampler
    variable = temperature
    boundary = 2
    sort_by = y
    outputs = 'outfile_temp_1'
  []
  [mass_fraction]
    type = LineValueSampler
    start_point = '0 0.01124 0'
    end_point = '0 0.47524 0'
    num_points = 30
    sort_by = y
    variable = layered_mass_fraction
    outputs = 'outfile_mass_1'
  []
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  csv = true
  color = false
  perf_graph = true
  [exodus]
    type = Exodus
    file_base = IFA_650_9_part1_out
    execute_on = 'initial timestep_end'
  []
  [checkpoint]
    type = Checkpoint
    time_step_interval =  1
    num_files = 1
  []
  [outfile_1]
    type = CSV
    execute_on = 'FINAL'
    create_final_symlink = true
  []
  [outfile_temp_1]
    type = CSV
    execute_on = 'FINAL'
    create_final_symlink = true
  []
  [outfile_mass_1]
    type = CSV
    execute_on = 'FINAL'
    create_final_symlink = true
  []
[]

(assessment/LWR/validation/LOCA_IFA_650/analysis/IFA_650_9/IFA_650_9_part1.i)


initial_fuel_density = 10430.0

[GlobalParams]
  density = ${initial_fuel_density}
  initial_porosity = 0.048
  order = SECOND
  family = LAGRANGE
  displacements = disp_x
  temperature = temperature
  energy_per_fission = 3.2e-11 #J/fission
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
  acceptable_multiplier = 10
[]

[Mesh]
  coord_type = RZ
  [layered1D_mesh]
    type = Layered1DMeshGenerator
    slices_per_block = 30
    pellet_outer_radius = 4.565e-3
    clad_gap_width = 0.085e-3
    clad_thickness = 0.725e-3
    fuel_height = 0.480
    plenum_height = 0.262416
    pellet_mesh_density = customize
    clad_mesh_density = customize
    nx_p = 11
    nx_c = 5
  []
  patch_update_strategy = auto
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[Variables]
  [disp_x]
  []
  [temperature]
    initial_condition = 295.0
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    data_file = power_history.csv
    scale_factor = 1.0
    format = columns
  []
  [axial_peaking_factors]
    type = PiecewiseBilinear
    data_file = axial_peaking_factors.csv
    axis = 1
    scale_factor = 1
  []
  [pressure_ramp]
    type = PiecewiseLinear
    data_file = coolant_pressure.csv
    scale_factor = 1
    format = columns
  []
  [average_htc]
    type = PiecewiseBilinear
    data_file = average_coolant_htc.csv
    axis = 1
    scale_factor = 1
  []
  [forced_times]
    type = PiecewiseLinear
    data_file = timestep_limiting.csv
    scale_factor = 1
    format = columns
  []
  [heat_sink_temperature]
    type = PiecewiseBilinear
    data_file = heater_temp.csv
    scale_factor = 1
    axis = 1
  []
  [clad_outer_temperature]
    type = PiecewiseBilinear
    data_file = clad_surface_temp.csv
    scale_factor = 1
    axis = 1
  []
  [heat_transfer_mode]
    type = PiecewiseConstant
    x = '-200  200412461 200413048'
    y = '9  9          8       '
    direction = 'right'
  []
  [clad_axial_pressure]
    type = CladdingAxialPressureFunction
    plenum_pressure = plenum_pressure
    coolant_pressure = pressure_ramp
    coolant_pressure_scaling_factor = 1.0
    fuel_pin_geometry = fuel_pin_geometry
  []
  [fuel_axial_pressure]
    type = ParsedFunction
    expression = plenum_pressure
    symbol_names = plenum_pressure
    symbol_values = plenum_pressure
  []
[]

[AuxVariables]
  [disp_y]
  []
  [disp_z]
  []
  [fast_neutron_flux]
    block = clad
  []
  [fast_neutron_fluence]
    block = clad
  []
  [grain_radius]
    block = fuel
    initial_condition = 5.0e-6
  []
  [effective_creep_strain]
    block = clad
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxide_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
  [gap_conductance]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_rate]
    order = CONSTANT
    family = MONOMIAL
  []
  [layered_maximum_clad_radius]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_source]
    type = NeutronHeatSource
    variable = temperature
    block = fuel
    burnup_function = burnup
    axial_relocation_object = axial_relocation
    extra_vector_tags = 'ref'
  []
[]

[Physics]
  [SolidMechanics]
    [Layered1D]
      [fuel]
        add_scalar_variables = true
        out_of_plane_strain_name = strain_yy
        fuel_pin_geometry = fuel_pin_geometry
        out_of_plane_pressure_function = fuel_axial_pressure
        strain = finite
        block = fuel
        eigenstrain_names = 'fuel_thermal_strain fuel_swelling_strain
          fuel_relocation_strain axial_relocation_eigenstrain'
        decomposition_method = EigenSolution
        cylindrical_axis_point1 = '0 0 0'
        cylindrical_axis_point2 = '0 1 0'
        generate_output = 'stress_xx stress_yy stress_zz vonmises_stress
          hoop_stress'
        extra_vector_tags = 'ref'
        group_scalar_vars_in_reference_residual = true
        mesh_generator = layered1D_mesh
      []
      [clad]
        add_scalar_variables = true
        out_of_plane_strain_name = strain_yy
        fuel_pin_geometry = fuel_pin_geometry
        strain = finite
        out_of_plane_pressure_function = clad_axial_pressure
        block = clad
        cylindrical_axis_point1 = '0 0 0'
        cylindrical_axis_point2 = '0 1 0'
        eigenstrain_names = 'clad_thermal_strain clad_irradiation_strain'
        generate_output = 'stress_xx stress_yy stress_zz vonmises_stress
          strain_zz creep_strain_zz hoop_stress'
        decomposition_method = EigenSolution
        extra_vector_tags = 'ref'
        group_scalar_vars_in_reference_residual = true
        mesh_generator = layered1D_mesh
      []
    []
  []
[]

[Burnup]
  [burnup]
    block = fuel
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    num_radial = 80
    num_axial = 11
    fuel_pin_geometry = fuel_pin_geometry
    fuel_volume_ratio = 1.0
    order = CONSTANT
    family = MONOMIAL
    RPF = RPF
    isotopes = 'U235 U238 Pu239 Pu240 Pu241 Pu242'
    isotope_fractions = '0.035 0.965 0 0 0 0'
  []
[]

[AuxKernels]
  [fast_neutron_flux]
    type = FastNeutronFluxAux
    block = clad
    variable = fast_neutron_flux
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    factor = 3e13
    execute_on = timestep_begin
  []
  [fast_neutron_fluence]
    type = FastNeutronFluenceAux
    block = clad
    variable = fast_neutron_fluence
    fast_neutron_flux = fast_neutron_flux
    execute_on = timestep_begin
  []
  [grain_radius]
    type = GrainRadiusAux
    block = fuel
    variable = grain_radius
    temperature = temperature
    execute_on = linear
  []
  [effective_creep_strain]
    type = MaterialRealAux
    block = clad
    variable = effective_creep_strain
    property = effective_creep_strain
    execute_on = 'timestep_end'
  []
  [fract_bphase]
    type = MaterialRealAux
    block = clad
    variable = fract_beta_phase
    property = fract_beta_phase
    execute_on = 'initial linear'
  []
  [oxide_thickness]
    type = MaterialRealAux
    boundary = 2
    variable = oxide_thickness
    property = oxide_scale_thickness
    execute_on = 'initial linear'
  []
  [hasburst]
    type = MaterialRealAux
    boundary = 2
    variable = burst
    property = failed
    execute_on = timestep_end
  []
  [conductance]
    type = MaterialRealAux
    boundary = 10
    property = gap_conductance
    variable = gap_conductance
    execute_on = 'initial linear'
  []
  [coolant_htc]
    type = MaterialRealAux
    property = coolant_channel_htc
    variable = coolant_htc
    boundary = 2
    execute_on = 'initial linear'
  []
  [creep_rate]
    type = MaterialRealAux
    block = clad
    variable = creep_rate
    property = creep_rate
    execute_on = timestep_end
  []
[]

[AxialRelocation]
  [rel]
    rod_ave_lin_pow = power_history
    axial_direction = y
    fuel_blocks = fuel
    clad_blocks = clad
    contact_pressure_variable = contact_pressure
    out_of_plane_strain_variable = strain_yy
    penetration_variable = penetration
    clad_inner_volume_addition = 0
    burnup_variable = burnup
    temperature = temperature
    gap_thickness_threshold = 0.00039
    axial_relocation_output_options = 'MASS_FRACTION'
    mesh_generator = layered1D_mesh
  []
[]

[CoolantChannel]
  [convective_clad_surface] # apply convective boundary to clad outer surface
    boundary = 2
    variable = temperature
    heat_transfer_mode = heat_transfer_mode
    heat_transfer_coefficient = average_htc # Calculated from an initial simulation of the base irradiation using the inlet_pressure, inlet_massflux, and inlet_temperature commented out below.
    inlet_temperature = heat_sink_temperature # K
    effective_emissivity = 0.75
    # inlet_temperature = 580
    # inlet_pressure    = 15.3e6   # Pa
    # inlet_massflux    = 3800     # kg/m^2-sec
    rod_diameter = 0.01075 # m
    rod_pitch = 1.26e-2 # m
    compute_enthalpy = false
    linear_heat_rate = power_history
    axial_power_profile = axial_peaking_factors
    output_properties = 'coolant_channel_htype coolant_channel_hmode'
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    penalty = 1e7
    formulation = kinematic
    model = frictionless
  []
[]

[ThermalContact]
  [thermal_contact]
    type = GasGapHeatTransfer
    variable = temperature
    primary = 5
    secondary = 10
    initial_moles = initial_moles
    gas_released = fis_gas_released
    plenum_pressure = plenum_pressure
    contact_pressure = contact_pressure
    jump_distance_model = LANNING
    roughness_coef = 3.2
    refab_gas_types = 'He Ar'
    refab_fractions = '0.05 0.95'
    refab_time = 199159200
    refab_type = 0
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 12
    value = 0.0
  []
  [Pressure]
    [coolantPressure]
      boundary = 2
      function = pressure_ramp
      factor = 1.0
    []
  []
  [clad_outer_temp]
    type = FunctionDirichletBC
    boundary = 2
    variable = temperature
    function = clad_outer_temperature
  []
  [PlenumPressure]
    [plenumPressure]
      boundary = 9
      initial_pressure = 2.0e6
      startup_time = 0
      R = 8.3143
      output_initial_moles = initial_moles
      temperature = plenum_temp
      volume = 'clad_volume pellet_volume'
      material_input = fis_gas_released
      output = plenum_pressure
      refab_time = 199159200
      refab_pressure = 4.0e6
      refab_temperature = 295.0
      refab_volume = 1.9e-05
    []
  []
[]

[LayeredPlenumTemperature]
  [plenum_temp]
    boundary = 5
    fuel_pin_geometry = fuel_pin_geometry
    out_of_plane_strain = strain_yy
    inner_surfaces = '5'
    outer_surfaces = '10'
    temperature = temperature
  []
[]

[Controls]
  [period0]
    type = TimePeriod
    disable_objects = 'BCs/clad_outer_temp'
    start_time = -200.0
    end_time = 199159200.0
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
    execute_on = timestep_end
  []
[]

[Materials]
  [fuel_thermal]
    type = UO2Thermal
    block = fuel
    thermal_conductivity_model = STAICU
    hbs_porosity_correction = KAMPF
    model_hbs_formation = true
    temperature = temperature
    burnup_function = burnup
    axial_relocation_object = axial_relocation
    gap_thermal_conductivity = layered_average_gap_conductivity
  []
  [relocation]
    type = UO2RelocationEigenstrain
    block = fuel
    burnup_function = burnup
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    fuel_pin_geometry = fuel_pin_geometry
    burnup_relocation_stop = 0.024
    relocation_activation1 = 5000.0
    relocation_model = ESCORE_modified
    eigenstrain_name = fuel_relocation_strain
  []
  [fuel_thermal_strain]
    type = ComputeThermalExpansionEigenstrain
    block = fuel
    thermal_expansion_coeff = 10.0e-6
    stress_free_temperature = 295.0
    eigenstrain_name = fuel_thermal_strain
  []
  [fuel_swelling]
    type = UO2VolumetricSwellingEigenstrain
    gas_swelling_model_type = SIFGRS
    block = fuel
    burnup_function = burnup
    initial_fuel_density = 10430.0
    eigenstrain_name = fuel_swelling_strain
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = fuel
    temperature = temperature
    burnup_function = burnup
    grain_radius = grain_radius
    transient_option = MICROCRACKING_BURNUP
    diff_coeff_option = TURNBULL_D1_D2
    gbs_model = true
  []
  [fuel_elasticity_tensor]
    type = UO2IsotropicDamageElasticityTensor
    block = fuel
    fragmentation_model = BARANI
    temperature = temperature
    rod_ave_lin_pow = power_history
    axial_relocation_object = axial_relocation
    crumbling_scale_factor = 0.0001
  []
  [fuel_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'fuel_creep'
    block = fuel
  []

  [fuel_creep]
    type = UO2CreepUpdate
    block = fuel
    temperature = temperature
    burnup_function = burnup
    initial_grain_radius = 5.0e-6
  []

  [HBS]
    type = HighBurnupStructureFormation
    block = fuel
    burnup_function = burnup
    temperature = temperature
    output_properties = 'hbs_volume_fraction'
    outputs = 'exodus'
  []

  [clad_elasticity_tensor]
    type = ZryElasticityTensor
    block = clad
  []
  [stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'zrycreep'
    block = clad
  []
  [zrycreep]
    type = ZryCreepLOCAUpdate
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    model_irradiation_creep = true
    model_primary_creep = true
    model_thermal_creep = true
    max_inelastic_increment = 5e-4
    zircaloy_material_type = stress_relief_annealed
    block = clad
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = clad
    stress_free_temperature = 295.0
    eigenstrain_name = clad_thermal_strain
  []
  [irradiation_swelling]
    type = ZryIrradiationGrowthEigenstrain
    block = clad
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = stress_relief_annealed
    eigenstrain_name = clad_irradiation_strain
  []
  [clad_phase]
    type = ZrPhase
    block = clad
    temperature = temperature
    numerical_method = 2
  []
  [clad_oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 4.65e-03
    clad_outer_radius = 5.375e-03
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = cathcart
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
  []
  [clad_thermal]
    type = ZryThermal
    block = clad
    temperature = temperature
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = fuel
    strain_free_density = ${initial_fuel_density}
  []
  [clad_density]
    type = StrainAdjustedDensity
    block = clad
    strain_free_density = 6551.0
  []
[]

[Postprocessors]
  [ave_temp_interior]
    type = LayeredSideAverageValuePostprocessor
    boundary = 9
    variable = temperature
    execute_on = 'initial linear'
    fuel_pin_geometry = fuel_pin_geometry
  []
  [pellet_volume_2]
    type = LayeredInternalVolumePostprocessor
    boundary = 8
    component = 0
    fuel_pin_geometry = fuel_pin_geometry
    out_of_plane_strain = strain_yy
    execute_on = 'initial linear'
  []
  [avg_clad_temp]
    type = LayeredSideAverageValuePostprocessor
    boundary = 7
    variable = temperature
    fuel_pin_geometry = fuel_pin_geometry
    execute_on = 'initial linear'
  []
  [fis_gas_produced]
    type = LayeredElementIntegralFisGasGeneratedSifgrsPostprocessor
    block = fuel
    fuel_pin_geometry = fuel_pin_geometry
  []
  [fis_gas_released]
    type = LayeredElementIntegralFisGasReleasedSifgrsPostprocessor
    block = fuel
    fuel_pin_geometry = fuel_pin_geometry
  []
  [fis_gas_grain]
    type = LayeredElementIntegralFisGasGrainSifgrsPostprocessor
    block = fuel
    outputs = exodus
    fuel_pin_geometry = fuel_pin_geometry
  []
  [fis_gas_boundary]
    type = LayeredElementIntegralFisGasBoundarySifgrsPostprocessor
    block = fuel
    outputs = exodus
    fuel_pin_geometry = fuel_pin_geometry
  []
  [fission_gas_release]
    type = FGRPercent
    fission_gas_released = fis_gas_released
    fission_gas_generated = fis_gas_produced
  []
  [average_coolant_htc]
    type = LayeredSideAverageValuePostprocessor
    boundary = 2
    variable = coolant_htc
    execute_on = 'initial linear'
    fuel_pin_geometry = fuel_pin_geometry
  []
  [average_burnup]
    type = RodAverageBurnup
    burnup_function = burnup
  []
  [temp_clad_max]
    type = NodalExtremeValue
    block = clad
    value_type = max
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [betaph_fract_max]
    type = ElementExtremeValue
    value_type = max
    variable = fract_beta_phase
    block = clad
    execute_on = 'initial timestep_end'
  []
  [burst]
    type = ElementExtremeValue
    value_type = max
    variable = burst
    block = clad
    execute_on = 'initial timestep_end'
  []
  [timestep_material]
    type = MaterialTimeStepPostprocessor
    block = clad
    execute_on = 'initial timestep_end'
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = clad
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geometry
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
[]

[Dampers]
  [limitT]
    type = BoundingValueNodalDamper
    variable = temperature
    max_value = 3200.0
    min_value = 0.0
  []
  [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 -ksp_gmres_restart'
  petsc_options_value = 'lu       superlu_dist                  51'

  line_search = 'none'

  l_max_its = 50
  l_tol = 1e-3
  nl_max_its = 30
  nl_rel_tol = 1e-4
  nl_abs_tol = 1e-8
  dtmax = 5e5
  dtmin = 1e-5
  start_time = -200.0
  end_time = 199159200 # End base irradiation
  #  end_time = 200412431 # Begin Blowdown
  # end_time = 200413048 # End

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 200
    timestep_limiting_postprocessor = timestep_material
    optimal_iterations = 20
    iteration_window = 4
    linear_iteration_ratio = 100
    timestep_limiting_function = forced_times
    force_step_every_function_point = true
    max_function_change = 2000
    time_t = '199159200 200312431 200411431 200412431 200412461 200413048'
    time_dt = '1.0e04    1.0e04    10.0        5.0         0.5        5.0'
  []
[]

[VectorPostprocessors]
  [clad_radial_disp]
    type = NodalValueSampler
    variable = disp_x
    boundary = 2
    sort_by = y
    outputs = 'outfile_1'
  []
  [clad_out_temp]
    type = NodalValueSampler
    variable = temperature
    boundary = 2
    sort_by = y
    outputs = 'outfile_temp_1'
  []
  [mass_fraction]
    type = LineValueSampler
    start_point = '0 0.01124 0'
    end_point = '0 0.47524 0'
    num_points = 30
    sort_by = y
    variable = layered_mass_fraction
    outputs = 'outfile_mass_1'
  []
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  csv = true
  color = false
  perf_graph = true
  [exodus]
    type = Exodus
    file_base = IFA_650_9_part1_out
    execute_on = 'initial timestep_end'
  []
  [checkpoint]
    type = Checkpoint
    time_step_interval =  1
    num_files = 1
  []
  [outfile_1]
    type = CSV
    execute_on = 'FINAL'
    create_final_symlink = true
  []
  [outfile_temp_1]
    type = CSV
    execute_on = 'FINAL'
    create_final_symlink = true
  []
  [outfile_mass_1]
    type = CSV
    execute_on = 'FINAL'
    create_final_symlink = true
  []
[]

(assessment/LWR/validation/LOCA_IFA_650/analysis/IFA_650_9/IFA_650_9_part1.i)


initial_fuel_density = 10430.0

[GlobalParams]
  density = ${initial_fuel_density}
  initial_porosity = 0.048
  order = SECOND
  family = LAGRANGE
  displacements = disp_x
  temperature = temperature
  energy_per_fission = 3.2e-11 #J/fission
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
  acceptable_multiplier = 10
[]

[Mesh]
  coord_type = RZ
  [layered1D_mesh]
    type = Layered1DMeshGenerator
    slices_per_block = 30
    pellet_outer_radius = 4.565e-3
    clad_gap_width = 0.085e-3
    clad_thickness = 0.725e-3
    fuel_height = 0.480
    plenum_height = 0.262416
    pellet_mesh_density = customize
    clad_mesh_density = customize
    nx_p = 11
    nx_c = 5
  []
  patch_update_strategy = auto
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[Variables]
  [disp_x]
  []
  [temperature]
    initial_condition = 295.0
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    data_file = power_history.csv
    scale_factor = 1.0
    format = columns
  []
  [axial_peaking_factors]
    type = PiecewiseBilinear
    data_file = axial_peaking_factors.csv
    axis = 1
    scale_factor = 1
  []
  [pressure_ramp]
    type = PiecewiseLinear
    data_file = coolant_pressure.csv
    scale_factor = 1
    format = columns
  []
  [average_htc]
    type = PiecewiseBilinear
    data_file = average_coolant_htc.csv
    axis = 1
    scale_factor = 1
  []
  [forced_times]
    type = PiecewiseLinear
    data_file = timestep_limiting.csv
    scale_factor = 1
    format = columns
  []
  [heat_sink_temperature]
    type = PiecewiseBilinear
    data_file = heater_temp.csv
    scale_factor = 1
    axis = 1
  []
  [clad_outer_temperature]
    type = PiecewiseBilinear
    data_file = clad_surface_temp.csv
    scale_factor = 1
    axis = 1
  []
  [heat_transfer_mode]
    type = PiecewiseConstant
    x = '-200  200412461 200413048'
    y = '9  9          8       '
    direction = 'right'
  []
  [clad_axial_pressure]
    type = CladdingAxialPressureFunction
    plenum_pressure = plenum_pressure
    coolant_pressure = pressure_ramp
    coolant_pressure_scaling_factor = 1.0
    fuel_pin_geometry = fuel_pin_geometry
  []
  [fuel_axial_pressure]
    type = ParsedFunction
    expression = plenum_pressure
    symbol_names = plenum_pressure
    symbol_values = plenum_pressure
  []
[]

[AuxVariables]
  [disp_y]
  []
  [disp_z]
  []
  [fast_neutron_flux]
    block = clad
  []
  [fast_neutron_fluence]
    block = clad
  []
  [grain_radius]
    block = fuel
    initial_condition = 5.0e-6
  []
  [effective_creep_strain]
    block = clad
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxide_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
  [gap_conductance]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_rate]
    order = CONSTANT
    family = MONOMIAL
  []
  [layered_maximum_clad_radius]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_source]
    type = NeutronHeatSource
    variable = temperature
    block = fuel
    burnup_function = burnup
    axial_relocation_object = axial_relocation
    extra_vector_tags = 'ref'
  []
[]

[Physics]
  [SolidMechanics]
    [Layered1D]
      [fuel]
        add_scalar_variables = true
        out_of_plane_strain_name = strain_yy
        fuel_pin_geometry = fuel_pin_geometry
        out_of_plane_pressure_function = fuel_axial_pressure
        strain = finite
        block = fuel
        eigenstrain_names = 'fuel_thermal_strain fuel_swelling_strain
          fuel_relocation_strain axial_relocation_eigenstrain'
        decomposition_method = EigenSolution
        cylindrical_axis_point1 = '0 0 0'
        cylindrical_axis_point2 = '0 1 0'
        generate_output = 'stress_xx stress_yy stress_zz vonmises_stress
          hoop_stress'
        extra_vector_tags = 'ref'
        group_scalar_vars_in_reference_residual = true
        mesh_generator = layered1D_mesh
      []
      [clad]
        add_scalar_variables = true
        out_of_plane_strain_name = strain_yy
        fuel_pin_geometry = fuel_pin_geometry
        strain = finite
        out_of_plane_pressure_function = clad_axial_pressure
        block = clad
        cylindrical_axis_point1 = '0 0 0'
        cylindrical_axis_point2 = '0 1 0'
        eigenstrain_names = 'clad_thermal_strain clad_irradiation_strain'
        generate_output = 'stress_xx stress_yy stress_zz vonmises_stress
          strain_zz creep_strain_zz hoop_stress'
        decomposition_method = EigenSolution
        extra_vector_tags = 'ref'
        group_scalar_vars_in_reference_residual = true
        mesh_generator = layered1D_mesh
      []
    []
  []
[]

[Burnup]
  [burnup]
    block = fuel
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    num_radial = 80
    num_axial = 11
    fuel_pin_geometry = fuel_pin_geometry
    fuel_volume_ratio = 1.0
    order = CONSTANT
    family = MONOMIAL
    RPF = RPF
    isotopes = 'U235 U238 Pu239 Pu240 Pu241 Pu242'
    isotope_fractions = '0.035 0.965 0 0 0 0'
  []
[]

[AuxKernels]
  [fast_neutron_flux]
    type = FastNeutronFluxAux
    block = clad
    variable = fast_neutron_flux
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    factor = 3e13
    execute_on = timestep_begin
  []
  [fast_neutron_fluence]
    type = FastNeutronFluenceAux
    block = clad
    variable = fast_neutron_fluence
    fast_neutron_flux = fast_neutron_flux
    execute_on = timestep_begin
  []
  [grain_radius]
    type = GrainRadiusAux
    block = fuel
    variable = grain_radius
    temperature = temperature
    execute_on = linear
  []
  [effective_creep_strain]
    type = MaterialRealAux
    block = clad
    variable = effective_creep_strain
    property = effective_creep_strain
    execute_on = 'timestep_end'
  []
  [fract_bphase]
    type = MaterialRealAux
    block = clad
    variable = fract_beta_phase
    property = fract_beta_phase
    execute_on = 'initial linear'
  []
  [oxide_thickness]
    type = MaterialRealAux
    boundary = 2
    variable = oxide_thickness
    property = oxide_scale_thickness
    execute_on = 'initial linear'
  []
  [hasburst]
    type = MaterialRealAux
    boundary = 2
    variable = burst
    property = failed
    execute_on = timestep_end
  []
  [conductance]
    type = MaterialRealAux
    boundary = 10
    property = gap_conductance
    variable = gap_conductance
    execute_on = 'initial linear'
  []
  [coolant_htc]
    type = MaterialRealAux
    property = coolant_channel_htc
    variable = coolant_htc
    boundary = 2
    execute_on = 'initial linear'
  []
  [creep_rate]
    type = MaterialRealAux
    block = clad
    variable = creep_rate
    property = creep_rate
    execute_on = timestep_end
  []
[]

[AxialRelocation]
  [rel]
    rod_ave_lin_pow = power_history
    axial_direction = y
    fuel_blocks = fuel
    clad_blocks = clad
    contact_pressure_variable = contact_pressure
    out_of_plane_strain_variable = strain_yy
    penetration_variable = penetration
    clad_inner_volume_addition = 0
    burnup_variable = burnup
    temperature = temperature
    gap_thickness_threshold = 0.00039
    axial_relocation_output_options = 'MASS_FRACTION'
    mesh_generator = layered1D_mesh
  []
[]

[CoolantChannel]
  [convective_clad_surface] # apply convective boundary to clad outer surface
    boundary = 2
    variable = temperature
    heat_transfer_mode = heat_transfer_mode
    heat_transfer_coefficient = average_htc # Calculated from an initial simulation of the base irradiation using the inlet_pressure, inlet_massflux, and inlet_temperature commented out below.
    inlet_temperature = heat_sink_temperature # K
    effective_emissivity = 0.75
    # inlet_temperature = 580
    # inlet_pressure    = 15.3e6   # Pa
    # inlet_massflux    = 3800     # kg/m^2-sec
    rod_diameter = 0.01075 # m
    rod_pitch = 1.26e-2 # m
    compute_enthalpy = false
    linear_heat_rate = power_history
    axial_power_profile = axial_peaking_factors
    output_properties = 'coolant_channel_htype coolant_channel_hmode'
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    penalty = 1e7
    formulation = kinematic
    model = frictionless
  []
[]

[ThermalContact]
  [thermal_contact]
    type = GasGapHeatTransfer
    variable = temperature
    primary = 5
    secondary = 10
    initial_moles = initial_moles
    gas_released = fis_gas_released
    plenum_pressure = plenum_pressure
    contact_pressure = contact_pressure
    jump_distance_model = LANNING
    roughness_coef = 3.2
    refab_gas_types = 'He Ar'
    refab_fractions = '0.05 0.95'
    refab_time = 199159200
    refab_type = 0
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 12
    value = 0.0
  []
  [Pressure]
    [coolantPressure]
      boundary = 2
      function = pressure_ramp
      factor = 1.0
    []
  []
  [clad_outer_temp]
    type = FunctionDirichletBC
    boundary = 2
    variable = temperature
    function = clad_outer_temperature
  []
  [PlenumPressure]
    [plenumPressure]
      boundary = 9
      initial_pressure = 2.0e6
      startup_time = 0
      R = 8.3143
      output_initial_moles = initial_moles
      temperature = plenum_temp
      volume = 'clad_volume pellet_volume'
      material_input = fis_gas_released
      output = plenum_pressure
      refab_time = 199159200
      refab_pressure = 4.0e6
      refab_temperature = 295.0
      refab_volume = 1.9e-05
    []
  []
[]

[LayeredPlenumTemperature]
  [plenum_temp]
    boundary = 5
    fuel_pin_geometry = fuel_pin_geometry
    out_of_plane_strain = strain_yy
    inner_surfaces = '5'
    outer_surfaces = '10'
    temperature = temperature
  []
[]

[Controls]
  [period0]
    type = TimePeriod
    disable_objects = 'BCs/clad_outer_temp'
    start_time = -200.0
    end_time = 199159200.0
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
    execute_on = timestep_end
  []
[]

[Materials]
  [fuel_thermal]
    type = UO2Thermal
    block = fuel
    thermal_conductivity_model = STAICU
    hbs_porosity_correction = KAMPF
    model_hbs_formation = true
    temperature = temperature
    burnup_function = burnup
    axial_relocation_object = axial_relocation
    gap_thermal_conductivity = layered_average_gap_conductivity
  []
  [relocation]
    type = UO2RelocationEigenstrain
    block = fuel
    burnup_function = burnup
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    fuel_pin_geometry = fuel_pin_geometry
    burnup_relocation_stop = 0.024
    relocation_activation1 = 5000.0
    relocation_model = ESCORE_modified
    eigenstrain_name = fuel_relocation_strain
  []
  [fuel_thermal_strain]
    type = ComputeThermalExpansionEigenstrain
    block = fuel
    thermal_expansion_coeff = 10.0e-6
    stress_free_temperature = 295.0
    eigenstrain_name = fuel_thermal_strain
  []
  [fuel_swelling]
    type = UO2VolumetricSwellingEigenstrain
    gas_swelling_model_type = SIFGRS
    block = fuel
    burnup_function = burnup
    initial_fuel_density = 10430.0
    eigenstrain_name = fuel_swelling_strain
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = fuel
    temperature = temperature
    burnup_function = burnup
    grain_radius = grain_radius
    transient_option = MICROCRACKING_BURNUP
    diff_coeff_option = TURNBULL_D1_D2
    gbs_model = true
  []
  [fuel_elasticity_tensor]
    type = UO2IsotropicDamageElasticityTensor
    block = fuel
    fragmentation_model = BARANI
    temperature = temperature
    rod_ave_lin_pow = power_history
    axial_relocation_object = axial_relocation
    crumbling_scale_factor = 0.0001
  []
  [fuel_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'fuel_creep'
    block = fuel
  []

  [fuel_creep]
    type = UO2CreepUpdate
    block = fuel
    temperature = temperature
    burnup_function = burnup
    initial_grain_radius = 5.0e-6
  []

  [HBS]
    type = HighBurnupStructureFormation
    block = fuel
    burnup_function = burnup
    temperature = temperature
    output_properties = 'hbs_volume_fraction'
    outputs = 'exodus'
  []

  [clad_elasticity_tensor]
    type = ZryElasticityTensor
    block = clad
  []
  [stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'zrycreep'
    block = clad
  []
  [zrycreep]
    type = ZryCreepLOCAUpdate
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    model_irradiation_creep = true
    model_primary_creep = true
    model_thermal_creep = true
    max_inelastic_increment = 5e-4
    zircaloy_material_type = stress_relief_annealed
    block = clad
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = clad
    stress_free_temperature = 295.0
    eigenstrain_name = clad_thermal_strain
  []
  [irradiation_swelling]
    type = ZryIrradiationGrowthEigenstrain
    block = clad
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = stress_relief_annealed
    eigenstrain_name = clad_irradiation_strain
  []
  [clad_phase]
    type = ZrPhase
    block = clad
    temperature = temperature
    numerical_method = 2
  []
  [clad_oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 4.65e-03
    clad_outer_radius = 5.375e-03
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = cathcart
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
  []
  [clad_thermal]
    type = ZryThermal
    block = clad
    temperature = temperature
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = fuel
    strain_free_density = ${initial_fuel_density}
  []
  [clad_density]
    type = StrainAdjustedDensity
    block = clad
    strain_free_density = 6551.0
  []
[]

[Postprocessors]
  [ave_temp_interior]
    type = LayeredSideAverageValuePostprocessor
    boundary = 9
    variable = temperature
    execute_on = 'initial linear'
    fuel_pin_geometry = fuel_pin_geometry
  []
  [pellet_volume_2]
    type = LayeredInternalVolumePostprocessor
    boundary = 8
    component = 0
    fuel_pin_geometry = fuel_pin_geometry
    out_of_plane_strain = strain_yy
    execute_on = 'initial linear'
  []
  [avg_clad_temp]
    type = LayeredSideAverageValuePostprocessor
    boundary = 7
    variable = temperature
    fuel_pin_geometry = fuel_pin_geometry
    execute_on = 'initial linear'
  []
  [fis_gas_produced]
    type = LayeredElementIntegralFisGasGeneratedSifgrsPostprocessor
    block = fuel
    fuel_pin_geometry = fuel_pin_geometry
  []
  [fis_gas_released]
    type = LayeredElementIntegralFisGasReleasedSifgrsPostprocessor
    block = fuel
    fuel_pin_geometry = fuel_pin_geometry
  []
  [fis_gas_grain]
    type = LayeredElementIntegralFisGasGrainSifgrsPostprocessor
    block = fuel
    outputs = exodus
    fuel_pin_geometry = fuel_pin_geometry
  []
  [fis_gas_boundary]
    type = LayeredElementIntegralFisGasBoundarySifgrsPostprocessor
    block = fuel
    outputs = exodus
    fuel_pin_geometry = fuel_pin_geometry
  []
  [fission_gas_release]
    type = FGRPercent
    fission_gas_released = fis_gas_released
    fission_gas_generated = fis_gas_produced
  []
  [average_coolant_htc]
    type = LayeredSideAverageValuePostprocessor
    boundary = 2
    variable = coolant_htc
    execute_on = 'initial linear'
    fuel_pin_geometry = fuel_pin_geometry
  []
  [average_burnup]
    type = RodAverageBurnup
    burnup_function = burnup
  []
  [temp_clad_max]
    type = NodalExtremeValue
    block = clad
    value_type = max
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [betaph_fract_max]
    type = ElementExtremeValue
    value_type = max
    variable = fract_beta_phase
    block = clad
    execute_on = 'initial timestep_end'
  []
  [burst]
    type = ElementExtremeValue
    value_type = max
    variable = burst
    block = clad
    execute_on = 'initial timestep_end'
  []
  [timestep_material]
    type = MaterialTimeStepPostprocessor
    block = clad
    execute_on = 'initial timestep_end'
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = clad
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geometry
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
[]

[Dampers]
  [limitT]
    type = BoundingValueNodalDamper
    variable = temperature
    max_value = 3200.0
    min_value = 0.0
  []
  [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 -ksp_gmres_restart'
  petsc_options_value = 'lu       superlu_dist                  51'

  line_search = 'none'

  l_max_its = 50
  l_tol = 1e-3
  nl_max_its = 30
  nl_rel_tol = 1e-4
  nl_abs_tol = 1e-8
  dtmax = 5e5
  dtmin = 1e-5
  start_time = -200.0
  end_time = 199159200 # End base irradiation
  #  end_time = 200412431 # Begin Blowdown
  # end_time = 200413048 # End

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 200
    timestep_limiting_postprocessor = timestep_material
    optimal_iterations = 20
    iteration_window = 4
    linear_iteration_ratio = 100
    timestep_limiting_function = forced_times
    force_step_every_function_point = true
    max_function_change = 2000
    time_t = '199159200 200312431 200411431 200412431 200412461 200413048'
    time_dt = '1.0e04    1.0e04    10.0        5.0         0.5        5.0'
  []
[]

[VectorPostprocessors]
  [clad_radial_disp]
    type = NodalValueSampler
    variable = disp_x
    boundary = 2
    sort_by = y
    outputs = 'outfile_1'
  []
  [clad_out_temp]
    type = NodalValueSampler
    variable = temperature
    boundary = 2
    sort_by = y
    outputs = 'outfile_temp_1'
  []
  [mass_fraction]
    type = LineValueSampler
    start_point = '0 0.01124 0'
    end_point = '0 0.47524 0'
    num_points = 30
    sort_by = y
    variable = layered_mass_fraction
    outputs = 'outfile_mass_1'
  []
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  csv = true
  color = false
  perf_graph = true
  [exodus]
    type = Exodus
    file_base = IFA_650_9_part1_out
    execute_on = 'initial timestep_end'
  []
  [checkpoint]
    type = Checkpoint
    time_step_interval =  1
    num_files = 1
  []
  [outfile_1]
    type = CSV
    execute_on = 'FINAL'
    create_final_symlink = true
  []
  [outfile_temp_1]
    type = CSV
    execute_on = 'FINAL'
    create_final_symlink = true
  []
  [outfile_mass_1]
    type = CSV
    execute_on = 'FINAL'
    create_final_symlink = true
  []
[]

Description
Constructed Objects
Example Input File Syntax
Input Parameters
Associated Actions
Available Actions