ZircaloyDamage

Computes response of Zircaloy cladding subject to damage and the effect of hydrides.

Description

ZircaloyDamage is a model for the mechanical effect of hydrides in Zircaloy cladding. It is based on Rashid et al. (2009).

The model takes as inputs elastic-plastic properties for hydrides and Zircaloy, the volume fraction of radial and circumferential hydrides, and a damage strain parameter for both Zircaloy and hydride. A combined response is developed based on the behavior of each constituent subject to certain constraints. Hydrides may appear in the circumferential direction, the radial direction, or both.

The basic building block of the model is the arrangement described in Figure 1 with intact and damaged hydride platelet and matrix materials. The fraction of these materials is $var element = document.getElementById("moose-equation-f876f0ab-8087-4c91-80f8-eadc8611b6e5");katex.render(" f_p+f_m=1 \\\\ f_1+f_3=f_p \\\\ f_2+f_4=f_m", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-84bee944-9e4b-4354-9316-c14d15d7669c");katex.render("f", element, {displayMode:false,throwOnError:false});$ is the fractional amount with subscripts $var element = document.getElementById("moose-equation-c133a627-74e4-4bdb-923e-6f17e42b6a5c");katex.render("p", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-86788cb0-efee-427d-aaf1-798af0188ac5");katex.render("m", element, {displayMode:false,throwOnError:false});$ , 1, 2, 3, and 4 representing platelet, matrix, and domains 1 through 4 based on Figure 1. As damage develops, one or both of $var element = document.getElementById("moose-equation-12712495-c18d-4392-9536-baa2f73e220f");katex.render("f_3", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-c5677cb2-af9c-4e25-8c02-ae62478dfc74");katex.render("f_4", element, {displayMode:false,throwOnError:false});$ are greater than zero.

Figure 1: Four domains in the damage model: 1, intact platelet; 2, intact matrix; 3, damaged platelet; and 4, damaged matrix. Based on Fig. 3 in Rashid et al. (2009).

Considering first the undamaged state, the two constituents are constrained to experience the same strains in the plane of the hydride and the same stresses out of the plane of the hydride. For undamaged platelet and matrix, $var element = document.getElementById("moose-equation-009d1c27-d27e-4656-84ad-2700696f7693");katex.render(" \\sigma'_{zz} = f_1 \\sigma'^p_{zz}+(1-f_1)\\sigma'^m_{zz}", element, {displayMode:true,throwOnError:false});$ where the prime indicates that the stress is oriented in the platelet's local coordinate frame.

More generally, $var element = document.getElementById("moose-equation-8ff1c835-a2e6-4bb6-be43-59aed57a5585");katex.render(" \\sigma = \\sum_{i=1}^4 f_i \\sigma^i.", element, {displayMode:true,throwOnError:false});$

The model enforces zero stress across damaged hydride and damaged Zircaloy phases. Any load perpendicular to damaged material is carried by the undamaged material.

Damage is measured by tracking the strain normal to the platelet. Increased strain increases the damage. With increasing damage, more material (hydride and Zircaloy) is considered damaged and less considered undamaged.

The damage function used in the model depends on the maximum normal strain, $var element = document.getElementById("moose-equation-5c8f7022-0f9b-49a1-b4ea-88fd242daa38");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ , the damage parameter $var element = document.getElementById("moose-equation-211e11be-807f-4e2e-80ea-53345f35c501");katex.render("\\epsilon_0", element, {displayMode:false,throwOnError:false});$ , a volume fraction, $var element = document.getElementById("moose-equation-5557f4b7-589a-42f6-b3a5-1a29a6229caa");katex.render("f", element, {displayMode:false,throwOnError:false});$ , and a damage fraction parameter, $var element = document.getElementById("moose-equation-eda1a81d-a378-40f7-bcc1-8ace29dc1b8d");katex.render("f_0", element, {displayMode:false,throwOnError:false});$ . Damage, $var element = document.getElementById("moose-equation-b40dc769-12eb-4dbf-84b6-8acd75952839");katex.render("d", element, {displayMode:false,throwOnError:false});$ , is given as $var element = document.getElementById("moose-equation-e3aebb9c-7894-4b56-9e83-fcffe882dff2");katex.render(" d= \\begin{dcases} 0 & \\text{for } \\epsilon < \\epsilon_0/2\\\\ \\left( \\frac{\\epsilon - \\epsilon_0 / 2}{\\bar\\epsilon - \\epsilon_0 / 2} \\right)^2 & \\text{for } \\epsilon_0/2 < \\epsilon < \\bar\\epsilon\\\\ 1 & \\text{for } \\epsilon > \\bar\\epsilon \\end{dcases}", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-ab9541a1-1666-4dc0-be90-bc62b2ef8217");katex.render(" \\bar\\epsilon = \\epsilon_0(1-f_0+f_0/f).", element, {displayMode:true,throwOnError:false});$

When the volume fraction, $var element = document.getElementById("moose-equation-1e13af8f-0d18-4bdf-9933-a55cadc9bd59");katex.render("f", element, {displayMode:false,throwOnError:false});$ , is one, $var element = document.getElementById("moose-equation-17700dd9-5a23-48d7-8440-42b093074d91");katex.render("\\bar\\epsilon = \\epsilon_0", element, {displayMode:false,throwOnError:false});$ with damage beginning at $var element = document.getElementById("moose-equation-a9d89bd0-3244-4656-92e4-2c216003aefa");katex.render("\\epsilon_0/2", element, {displayMode:false,throwOnError:false});$ and ending at $var element = document.getElementById("moose-equation-ab21b794-eaf2-41db-a423-de50e994dd8a");katex.render("\\epsilon_0", element, {displayMode:false,throwOnError:false});$ through a quadratic transition. With a smaller, more typical volume fraction, damage begins at the same point but reaches one at a larger strain.

Overall damage in the matrix is computed as the sum of the damage in the hydride squared and a separate damage calculation for the matrix. In this way, matrix damage is often lower than the hydride damage but reaches one when the hydride is fully damaged.

By considering radial hydrides alongside a matrix composed of Zircaloy with circumferential hydrides and also circumferential hydrides alongside a matrix composed of Zircaloy with radial hydrides, the interaction between all three is determined.

In this model, damage is driven purely by tensile strains in the radial and circumferential directions. With no tensile strain in these directions, no damage will be computed.

Note that this material is based on an elastic-plastic constitutive relationship and not on creep. The parameters reference_strain_hydride and reference_strain_zircaloy are expected to be greater than zero and default to 0.001. The yield stress is $var element = document.getElementById("moose-equation-444d464d-77ff-4a9a-abf8-710081161599");katex.render(" Y=Y_0*(1+(\\epsilon/\\epsilon_{ref})^n)", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-6ef2d81e-ef19-4962-816d-f5e2723c4326");katex.render("Y", element, {displayMode:false,throwOnError:false});$ is the yield stress, $var element = document.getElementById("moose-equation-4a83095c-07d3-4bb4-ab16-ed5200112414");katex.render("Y_0", element, {displayMode:false,throwOnError:false});$ is the initial yield stress (initial_yield_stress_hydride or initial_yield_stress_zircaloy), $var element = document.getElementById("moose-equation-1e1e335e-dc05-48f6-b656-ad42fefc4eda");katex.render("\\epsilon", element, {displayMode:false,throwOnError:false});$ is the equivalent plastic strain, $var element = document.getElementById("moose-equation-e8eb1172-31ac-4488-8c50-5dd308883073");katex.render("\\epsilon_{ref}", element, {displayMode:false,throwOnError:false});$ is the reference strain, and $var element = document.getElementById("moose-equation-c4149897-9ee9-49a3-aef3-b0262c8d6697");katex.render("n", element, {displayMode:false,throwOnError:false});$ is the hardening exponent.

For a complete description of the model, see Rashid et al. (2009).

Example Input Syntax

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [ZircaloyDamage]
    type = ZircaloyDamage<<<{"description": "Computes response of Zircaloy cladding subject to damage and the effect of hydrides.", "href": "ZircaloyDamage.html"}>>>

    elastic_modulus_hydride<<<{"description": "elastic modulus for hydride inclusions"}>>> = 9.67e10
    poisson_ratio_hydride<<<{"description": "poisson ratio for hydride inclusions"}>>> = 0.389
    K_hydride<<<{"description": "K parameter for hydride inclusions.  Equals initial yield stress if 'reference_strain_hydride' is greater than zero."}>>> = 411e6
    reference_strain_hydride<<<{"description": "reference strain in plasticity model for hydride inclusions"}>>> = 1e-3
    hardening_exponent_hydride<<<{"description": "hardening exponent in plasticity model for hydride inclusions"}>>> = 0
    damage_strain_hydride<<<{"description": "reference normal strain for damage in hydride"}>>> = 0.004

    K_zircaloy<<<{"description": "K parameter for zirconium matrix.  Equals initial yield stress if 'reference_strain_hydride' is greater than zero."}>>> = 4e8
    reference_strain_zircaloy<<<{"description": "reference strain in plasticity model for zirconium matrix"}>>> = 1e-3
    hardening_exponent_zircaloy<<<{"description": "hardening exponent in plasticity model for zirconium matrix"}>>> = 0.05
    damage_strain_zircaloy<<<{"description": "reference normal strain for damage in zirconium"}>>> = 0.035

    volume_fraction_radial<<<{"description": "volume fraction of hydride platelets in the radial direction"}>>> = 0.005
    volume_fraction_circumferential<<<{"description": "volume fraction of hydride platelets in the circumferential direction"}>>> = 0.067
    scale_circumferential<<<{"description": "Whether to scale the distribution of circumferential hydrides, with fewer toward the inner radius and more toward the outer radius."}>>> = true
    clad_inner_radius<<<{"description": "Inner cladding radius (m)"}>>> = 0.004125
    clad_outer_radius<<<{"description": "Outer cladding radius (m)"}>>> = 0.004535
  []
[]

Input Parameters

K_hydrideK parameter for hydride inclusions. Equals initial yield stress if 'reference_strain_hydride' is greater than zero.
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:K parameter for hydride inclusions. Equals initial yield stress if 'reference_strain_hydride' is greater than zero.
K_zircaloyK parameter for zirconium matrix. Equals initial yield stress if 'reference_strain_hydride' is greater than zero.
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:K parameter for zirconium matrix. Equals initial yield stress if 'reference_strain_hydride' is greater than zero.
damage_strain_hydridereference normal strain for damage in hydride
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:reference normal strain for damage in hydride
damage_strain_zircaloyreference normal strain for damage in zirconium
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:reference normal strain for damage in zirconium
elastic_modulus_hydrideelastic modulus for hydride inclusions
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:elastic modulus for hydride inclusions
hardening_exponent_hydridehardening exponent in plasticity model for hydride inclusions
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:hardening exponent in plasticity model for hydride inclusions
hardening_exponent_zircaloyhardening exponent in plasticity model for zirconium matrix
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:hardening exponent in plasticity model for zirconium matrix
poisson_ratio_hydridepoisson ratio for hydride inclusions
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:poisson ratio for hydride inclusions

Required Parameters

base_nameOptional parameter that allows the user to define multiple mechanics material systems on the same block, i.e. for multiple phases
C++ Type:std::string
Controllable:No
Description:Optional parameter that allows the user to define multiple mechanics material systems on the same block, i.e. for multiple phases
blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
clad_inner_radius0Inner cladding radius (m)
Default:0
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Inner cladding radius (m)
clad_outer_radius0Outer cladding radius (m)
Default:0
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Outer cladding radius (m)
computeTrueWhen false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
Default:True
C++ Type:bool
Controllable:No
Description:When false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
constant_onNONEWhen ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
Default:NONE
C++ Type:MooseEnum
Options:NONE, ELEMENT, SUBDOMAIN
Controllable:No
Description:When ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
debug0Debugging output flag: 0 = no output (default), 1 = simple debugging output to console
Default:0
C++ Type:int
Controllable:No
Description:Debugging output flag: 0 = no output (default), 1 = simple debugging output to console
declare_suffixAn optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any declared properties. The suffix will be prepended with a '_' character.
fuel_pin_geometryName of the UserObject that reads the pin geometry from the mesh.
C++ Type:UserObjectName
Controllable:No
Description:Name of the UserObject that reads the pin geometry from the mesh.
normal_circumferential1 0 0unit normal for circumferential hydrides
Default:1 0 0
C++ Type:libMesh::VectorValue<double>
Unit:(no unit assumed)
Controllable:No
Description:unit normal for circumferential hydrides
normal_radial0 0 1unit normal for radial hydrides
Default:0 0 1
C++ Type:libMesh::VectorValue<double>
Unit:(no unit assumed)
Controllable:No
Description:unit normal for radial hydrides
reference_strain_hydride0.001reference strain in plasticity model for hydride inclusions
Default:0.001
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:reference strain in plasticity model for hydride inclusions
reference_strain_zircaloy0.001reference strain in plasticity model for zirconium matrix
Default:0.001
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:reference strain in plasticity model for zirconium matrix
scale_circumferentialFalseWhether to scale the distribution of circumferential hydrides, with fewer toward the inner radius and more toward the outer radius.
Default:False
C++ Type:bool
Controllable:No
Description:Whether to scale the distribution of circumferential hydrides, with fewer toward the inner radius and more toward the outer radius.
volume_fraction_circumferentialvolume fraction of hydride platelets in the circumferential direction
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:volume fraction of hydride platelets in the circumferential direction
volume_fraction_circumferential_propertymaterial property with volume fraction of hydride platelets in the circumferential direction
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:material property with volume fraction of hydride platelets in the circumferential direction
volume_fraction_cross0.05reference volume fraction for cross-damage in zirconium
Default:0.05
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:reference volume fraction for cross-damage in zirconium
volume_fraction_radialvolume fraction of hydride platelets in the radial direction
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:volume fraction of hydride platelets in the radial direction
volume_fraction_radial_propertymaterial property with volume fraction of hydride platelets in the radial direction
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:material property with volume fraction of hydride platelets in the radial direction
volume_fraction_self0.05reference volume fraction for self-damage in hydride
Default:0.05
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:reference volume fraction for self-damage in hydride
zircaloy_modelOptional material model to use to compute the response of the Zircaloy matrix.
C++ Type:MaterialName
Controllable:No
Description:Optional material model to use to compute the response of the Zircaloy matrix.

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.
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator

Advanced Parameters

output_propertiesList of material properties, from this material, to output (outputs must also be defined to an output type)
C++ Type:std::vector<std::string>
Controllable:No
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
outputsnone Vector of output names where you would like to restrict the output of variables(s) associated with this object
Default:none
C++ Type:std::vector<OutputName>
Controllable:No
Description:Vector of output names where you would like to restrict the output of variables(s) associated with this object

Outputs Parameters

prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.

Material Property Retrieval Parameters

Input Files

(test/tests/zircaloy_damage/plastic_independent.i)
(test/tests/zircaloy_damage/radial_crack.i)
(test/tests/zircaloy_damage/burst_one_elem.i)
(test/tests/zircaloy_damage/elastic_rotation.i)
(test/tests/zircaloy_damage/homogeneous_elastic.i)
(test/tests/solid_mechanics/zry_creep/creep_sra_thermal_only_damage.i)
(test/tests/zircaloy_damage/radial_crack_plus_y.i)
(test/tests/zircaloy_damage/plastic.i)
(test/tests/zircaloy_damage/elastic_patch.i)
(test/tests/zircaloy_damage/homogeneous_elastic_2.i)

References

Joe Rashid, Mark Rashid, Albert Machiels, and Robert Dunham. A new material constitutive model for predicting cladding failure. In Proceedings of Top Fuel 2009. Paris, France, September 6-10 2009.

@inproceedings{rashid_hydride_topfuel_2009,
    Author = "Rashid, Joe and Rashid, Mark and Machiels, Albert and Dunham, Robert",
    Address = "Paris, France",
    Booktitle = "Proceedings of Top Fuel 2009",
    Month = "September 6-10",
    Title = "A New Material Constitutive Model for Predicting Cladding Failure",
    Year = "2009"
}

(test/tests/zircaloy_damage/burst_one_elem.i)

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

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

[Mesh]
  coord_type = RZ
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 1
    xmin = 0.004125
    xmax = 0.004535
    ymax = 2.05e-5
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [cladding]
    add_variables = true
    strain = FINITE
    incremental = true
    extra_vector_tags = 'ref'
    generate_output = 'strain_zz stress_zz'
  []
[]

[BCs]
  [y_fixed]
    # plane strain
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom top'
    value = 0.0
  []
  [push]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = left
    function = 1e-7*t
  []
[]

[Materials]
  [zircaloy_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 8.5e10
    poissons_ratio = 0.3
  []
  [ZircaloyDamage]
    type = ZircaloyDamage

    elastic_modulus_hydride = 9.67e10
    poisson_ratio_hydride = 0.389
    K_hydride = 411e6
    reference_strain_hydride = 1e-3
    hardening_exponent_hydride = 0
    damage_strain_hydride = 0.004

    K_zircaloy = 4e8
    reference_strain_zircaloy = 1e-3
    hardening_exponent_zircaloy = 0.05
    damage_strain_zircaloy = 0.035

    volume_fraction_radial = 0.005
    volume_fraction_circumferential = 0.067
    scale_circumferential = true
    clad_inner_radius = 0.004125
    clad_outer_radius = 0.004535
  []
[]

[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  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
  nl_max_its = 50
  nl_rel_tol = 1e-3
  nl_abs_tol = 1e-5
  l_tol = 1e-2
  start_time = 0.0
  end_time = 3.3e3
  dt = 15.0
  dtmin = 1e-6
[]

[Outputs]
  time_step_interval =  1
  exodus = true
  csv = true
[]

[Postprocessors]
  [strain_zz]
    type = ElementalVariableValue
    variable = strain_zz
    elementid = 0
  []
  [stress_zz]
    type = ElementalVariableValue
    variable = stress_zz
    elementid = 0
  []
[]

(test/tests/zircaloy_damage/plastic_independent.i)

# This is a test of the isotropic power law hardening portion of the model.
#
# In this problem, a single Hex 8 element is fixed at the bottom and pulled at the top
# at a constant rate of 0.1.
#
# Before yield, stress = strain (=0.1*t) as youngs modulus is 1.0.
#
# For the hydride, the K value is 0.25, and the reference strain and exponent
# are set to 0.5.  In this case, the model uses
#   Y = K * (1 + pow(eps / eps_0, n))
# to compute the yield stress where K is the K value, eps is equivalent plastic
# strain, eps_0 is the reference strain, and n is the exponent.
#
# A solid_mechanics model is used for the Zircaloy with parameters to match
# the behavior of the hydride.
#
# Therefore, the material should start yielding at t = 2.5 seconds and then follow
# stress = 0.25 * (1 + pow(eps / 0.5, 0.5)).

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

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

[Mesh]
  use_displaced_mesh = false
  [mesh]
    type = GeneratedMeshGenerator
    dim = 3
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [pellets]
    add_variables = true
    strain = SMALL
    incremental = true
    extra_vector_tags = 'ref'
    generate_output = 'stress_xx stress_xy stress_yy stress_zz strain_yy'
  []
[]

[Functions]
  [disp_pull]
    type = ParsedFunction
    expression = 0.1*t
  []
[]

[BCs]
  [x_face]
    type = DirichletBC
    variable = disp_x
    boundary = left
    value = 0.0
  []

  [y_face]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  []

  [y_pull]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top
    function = disp_pull
  []

  [z_face]
    type = DirichletBC
    variable = disp_z
    boundary = back
    value = 0.0
  []
[]

[Materials]
  [zircaloy_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1
    poissons_ratio = 0.3
  []
  [ZircDamage]
    type = ZircaloyDamage

    elastic_modulus_hydride = 1
    poisson_ratio_hydride = 0.3
    K_hydride = 0.25
    reference_strain_hydride = 5e-1
    hardening_exponent_hydride = 0.5

    K_zircaloy = 0.25
    reference_strain_zircaloy = 5e-1
    hardening_exponent_zircaloy = 0.5

    damage_strain_hydride = 2.e2
    volume_fraction_self = 5.e-2
    volume_fraction_cross = 5.e-2
    damage_strain_zircaloy = 2.e2
    volume_fraction_radial = 8.e-2
    volume_fraction_circumferential = 8.e-2
    normal_radial = '0.0 0.0 1.0'
    normal_circumferential = '0.0 1.0 0.0'
    clad_inner_radius = 0
    clad_outer_radius = 1
    debug = 0
    zircaloy_model = flubber
  []
  [flubber]
    type = IsotropicPlasticityStressUpdate
    temperature = 0
    yield_stress = 0.25
    hardening_function = 'if(t<=0,0.25,0.25*(1+pow(t/5e-1,0.5)))'
  []
[]

[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  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
  nl_max_its = 50
  nl_rel_tol = 1e-3
  nl_abs_tol = 1e-6
  l_tol = 1e-2
  start_time = 0.0
  end_time = 5
  dt = 0.25
  dtmin = 1e-6
[]

[Outputs]
  time_step_interval =  1
  exodus = true
  print_linear_residuals = true
[]

[Postprocessors]
  [stress_yy]
    type = ElementAverageValue
    variable = stress_yy
  []
  [strain_yy]
    type = ElementAverageValue
    variable = strain_yy
  []
[]

(test/tests/zircaloy_damage/radial_crack.i)

#
# This test pulls a single element in the z direction until failure.
#
[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

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

[Mesh]
  use_displaced_mesh = false
  [mesh]
    type = GeneratedMeshGenerator
    dim = 3
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [cube]
    add_variables = true
    strain = SMALL
    incremental = true
    extra_vector_tags = 'ref'
    generate_output = 'stress_xx stress_xy stress_yy stress_zz strain_xx strain_yy strain_zz'
  []
[]

[BCs]
  [x_fixed]
    type = DirichletBC
    variable = disp_x
    boundary = 'left right'
    value = 0.0
  []
  [y_fixed]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom top'
    value = 0.0
  []
  [back_fixed]
    type = DirichletBC
    variable = disp_z
    boundary = back
    value = 0.0
  []

  [front_pull]
    type = FunctionDirichletBC
    variable = disp_z
    boundary = front
    function = '0.1*t'
  []
[]

[Materials]
  [zircaloy_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 5e10
    poissons_ratio = 0.25
  []
  [ZircDamage]
    type = ZircaloyDamage
    elastic_modulus_hydride = 10e10
    poisson_ratio_hydride = 0.25
    K_hydride = 1e30
    reference_strain_hydride = 1e-3
    hardening_exponent_hydride = 8e-2
    damage_strain_hydride = 3e-2

    K_zircaloy = 5e8
    reference_strain_zircaloy = 1e-3
    hardening_exponent_zircaloy = 0.05
    damage_strain_zircaloy = 100

    volume_fraction_radial = 0.5
    volume_fraction_circumferential = 0.0
    scale_circumferential = false
    clad_inner_radius = 0
    clad_outer_radius = 1
    debug = 0
  []
[]

[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  solve_type = 'NEWTON'
  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
  nl_max_its = 50
  nl_rel_tol = 1e-10
  nl_abs_tol = 1e-8
  l_tol = 1e-2
  start_time = 0.0
  end_time = 0.7
  dt = 0.01
  dtmin = 1e-6
[]

[Outputs]
  time_step_interval =  1
  exodus = true
  print_linear_residuals = true
[]

[Postprocessors]
  [stress_xx]
    type = ElementAverageValue
    variable = stress_xx
  []
  [strain_xx]
    type = ElementAverageValue
    variable = strain_xx
  []
  [stress_yy]
    type = ElementAverageValue
    variable = stress_yy
  []
  [strain_yy]
    type = ElementAverageValue
    variable = strain_yy
  []
  [stress_zz]
    type = ElementAverageValue
    variable = stress_zz
  []
  [strain_zz]
    type = ElementAverageValue
    variable = strain_zz
  []
[]

(test/tests/zircaloy_damage/burst_one_elem.i)

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

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

[Mesh]
  coord_type = RZ
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 1
    xmin = 0.004125
    xmax = 0.004535
    ymax = 2.05e-5
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [cladding]
    add_variables = true
    strain = FINITE
    incremental = true
    extra_vector_tags = 'ref'
    generate_output = 'strain_zz stress_zz'
  []
[]

[BCs]
  [y_fixed]
    # plane strain
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom top'
    value = 0.0
  []
  [push]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = left
    function = 1e-7*t
  []
[]

[Materials]
  [zircaloy_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 8.5e10
    poissons_ratio = 0.3
  []
  [ZircaloyDamage]
    type = ZircaloyDamage

    elastic_modulus_hydride = 9.67e10
    poisson_ratio_hydride = 0.389
    K_hydride = 411e6
    reference_strain_hydride = 1e-3
    hardening_exponent_hydride = 0
    damage_strain_hydride = 0.004

    K_zircaloy = 4e8
    reference_strain_zircaloy = 1e-3
    hardening_exponent_zircaloy = 0.05
    damage_strain_zircaloy = 0.035

    volume_fraction_radial = 0.005
    volume_fraction_circumferential = 0.067
    scale_circumferential = true
    clad_inner_radius = 0.004125
    clad_outer_radius = 0.004535
  []
[]

[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  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
  nl_max_its = 50
  nl_rel_tol = 1e-3
  nl_abs_tol = 1e-5
  l_tol = 1e-2
  start_time = 0.0
  end_time = 3.3e3
  dt = 15.0
  dtmin = 1e-6
[]

[Outputs]
  time_step_interval =  1
  exodus = true
  csv = true
[]

[Postprocessors]
  [strain_zz]
    type = ElementalVariableValue
    variable = strain_zz
    elementid = 0
  []
  [stress_zz]
    type = ElementalVariableValue
    variable = stress_zz
    elementid = 0
  []
[]

(test/tests/zircaloy_damage/elastic_rotation.i)

#
# Rotation Test
#
# This test is designed to compute a uniaxial stress and then follow that
# stress as the mesh is rotated 90 degrees.
#
# The mesh is composed of one block with a single element.  The nodal
# displacements in the x and y directions are prescribed.  Poisson's
# ratio is zero.
#

[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 3
    xmin = 0
    ymin = 0
    zmin = 0
    xmax = 1
    ymax = 1
    zmax = 1
  []
  [bottom_left]
    type = BoundingBoxNodeSetGenerator
    input = mesh
    bottom_left = '-1 -1 -1'
    top_right = '.5 .5 2'
    new_boundary = ns_100
  []
  [bottom_right]
    type = BoundingBoxNodeSetGenerator
    input = bottom_left
    bottom_left = '.5 -1 -1'
    top_right = '2 .5 2'
    new_boundary = ns_200
  []
  [top_right]
    type = BoundingBoxNodeSetGenerator
    input = bottom_right
    bottom_left = '.5 .5 -1'
    top_right = '2 2 2'
    new_boundary = ns_300
  []
  [top_left]
    type = BoundingBoxNodeSetGenerator
    input = top_right
    bottom_left = '-1 .5 -1'
    top_right = '.5 2 2'
    new_boundary = ns_400
  []
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[Functions]
  [x_200]
    type = ParsedFunction
    symbol_names = 'delta t0'
    symbol_values = '-1e-6 1.0'
    expression = 'if(t<=1.0, delta*t, (1.0+delta)*cos(pi/2*(t-t0)) - 1.0)'
  []
  [y_200]
    type = ParsedFunction
    symbol_names = 'delta t0'
    symbol_values = '-1e-6 1.0'
    expression = 'if(t<=1.0, 0.0, (1.0+delta)*sin(pi/2*(t-t0)))'
  []
  [x_300]
    type = ParsedFunction
    symbol_names = 'delta t0'
    symbol_values = '-1e-6 1.0'
    expression = 'if(t<=1.0, delta*t, (1.0+delta)*cos(pi/2.0*(t-t0)) - sin(pi/2.0*(t-t0)) - 1.0)'
  []
  [y_300]
    type = ParsedFunction
    symbol_names = 'delta t0'
    symbol_values = '-1e-6 1.0'
    expression = 'if(t<=1.0, 0.0, cos(pi/2.0*(t-t0)) + (1+delta)*sin(pi/2.0*(t-t0)) - 1.0)'
  []
  [x_400]
    type = ParsedFunction
    symbol_names = 'delta t0'
    symbol_values = '-1e-6 1.0'
    expression = 'if(t<=1.0, 0.0, -sin(pi/2.0*(t-t0)))'
  []
  [y_400]
    type = ParsedFunction
    symbol_names = 'delta t0'
    symbol_values = '-1e-6 1.0'
    expression = 'if(t<=1.0, 0.0, cos(pi/2.0*(t-t0)) - 1.0)'
  []
[]

[Variables]
  [disp_x]
    order = FIRST
    family = LAGRANGE
  []
  [disp_y]
    order = FIRST
    family = LAGRANGE
  []
  [disp_z]
    order = FIRST
    family = LAGRANGE
  []
[]

[AuxVariables]
  [normal1_x]
    order = CONSTANT
    family = MONOMIAL
  []
  [normal1_y]
    order = CONSTANT
    family = MONOMIAL
  []
  [normal1_z]
    order = CONSTANT
    family = MONOMIAL
  []
  [normal2_x]
    order = CONSTANT
    family = MONOMIAL
  []
  [normal2_y]
    order = CONSTANT
    family = MONOMIAL
  []
  [normal2_z]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [fred]
    generate_output = 'stress_xx stress_yy stress_xy'
    strain = FINITE
  []
[]

[AuxKernels]
  [normal1_x]
    type = MaterialRealVectorValueAux
    property = unit_norm_to_p1p2
    component = 0
    variable = normal1_x
  []
  [normal1_y]
    type = MaterialRealVectorValueAux
    property = unit_norm_to_p1p2
    component = 1
    variable = normal1_y
  []
  [normal1_z]
    type = MaterialRealVectorValueAux
    property = unit_norm_to_p1p2
    component = 2
    variable = normal1_z
  []
  [normal2_x]
    type = MaterialRealVectorValueAux
    property = unit_norm_to_p3p4p5p6
    component = 0
    variable = normal2_x
  []
  [normal2_y]
    type = MaterialRealVectorValueAux
    property = unit_norm_to_p3p4p5p6
    component = 1
    variable = normal2_y
  []
  [normal2_z]
    type = MaterialRealVectorValueAux
    property = unit_norm_to_p3p4p5p6
    component = 2
    variable = normal2_z
  []
[]

[BCs]
  [no_x]
    type = DirichletBC
    variable = disp_x
    boundary = ns_100
    value = 0.0
  []
  [no_y]
    type = DirichletBC
    variable = disp_y
    boundary = ns_100
    value = 0.0
  []
  [x_200]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = ns_200
    function = x_200
  []
  [y_200]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = ns_200
    function = y_200
  []
  [x_300]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = ns_300
    function = x_300
  []
  [y_300]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = ns_300
    function = y_300
  []
  [x_400]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = ns_400
    function = x_400
  []
  [y_400]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = ns_400
    function = y_400
  []
  [no_z]
    type = DirichletBC
    variable = disp_z
    boundary = 'ns_100 ns_200 ns_300 ns_400'
    value = 0.0
  []
[]

[Materials]
  [zircaloy_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e6
    poissons_ratio = 0
  []
  [ZircDamage]
    type = ZircaloyDamage

    elastic_modulus_hydride = 1e6
    poisson_ratio_hydride = 0
    K_hydride = 1e20
    reference_strain_hydride = 5.e-2
    hardening_exponent_hydride = 8.e-2

    K_zircaloy = 1e20
    reference_strain_zircaloy = 5.e-2
    hardening_exponent_zircaloy = 8.e-2

    damage_strain_hydride = 2.e2
    volume_fraction_self = 5.e-2
    volume_fraction_cross = 5.e-2
    damage_strain_zircaloy = 2.e2
    volume_fraction_radial = 8.e-2
    volume_fraction_circumferential = 8.e-2
    normal_radial = '1.0 0.0 0.0'
    normal_circumferential = '0.0 1.0 0.0' # in y direction for testing
    clad_inner_radius = 0
    clad_outer_radius = 1
    debug = 0
  []
[]

[Postprocessors]
  [stress_xx]
    type = ElementIntegralMaterialProperty
    mat_prop = stress_xx
  []
  [stress_yy]
    type = ElementIntegralMaterialProperty
    mat_prop = stress_yy
  []
  [stress_xy]
    type = ElementIntegralMaterialProperty
    mat_prop = stress_xy
  []
  [normal1_x]
    type = ElementAverageValue
    variable = normal1_x
  []
  [normal1_y]
    type = ElementAverageValue
    variable = normal1_y
  []
  [normal1_z]
    type = ElementAverageValue
    variable = normal1_z
  []
  [normal2_x]
    type = ElementAverageValue
    variable = normal2_x
  []
  [normal2_y]
    type = ElementAverageValue
    variable = normal2_y
  []
  [normal2_z]
    type = ElementAverageValue
    variable = normal2_z
  []
[]

[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient

  solve_type = 'NEWTON'

  petsc_options_iname = '-pc_type '
  petsc_options_value = 'lu'
  nl_rel_tol = 1e-30
  nl_abs_tol = 1e-20
  l_max_its = 20
  start_time = 0.0
  dt = 0.01
  end_time = 2.0
[]

[Outputs]
  exodus = true
[]

(test/tests/zircaloy_damage/homogeneous_elastic.i)

#
# This problem is taken from the Abaqus verification manual:
#   "1.5.4 Patch test for axisymmetric elements"
# The stress solution is given as:
#   xx = yy = zz = 2000
#   xy = 400
#
# Since the strain is 1e-3 in all three directions, the new density should be
#   new_density = original_density * V_0 / V
#   new_density = 0.283 / (1 + 1e-3 + 1e-3 + 1e-3) = 0.282153

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  coord_type = RZ
  use_displaced_mesh = false
  [mesh]
    type = ExamplePatchMeshGenerator
    dim = 2
    x_offset = 1
  []
[]

[Physics/SolidMechanics/QuasiStatic/All]
  strain = SMALL
  incremental = true
  add_variables = true
  generate_output = 'stress_xx stress_yy stress_zz stress_xy stress_yz stress_zx'
[]

[Kernels]
  [body]
    type = BodyForce
    variable = disp_y
    value = 1
    function = '-400/x'
  []
[]

[BCs]
  [ur]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = 'left right'
    function = '1e-3*x'
  []
  [uz]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = 'left right'
    function = '1e-3*(x+y)'
  []
[]

[Materials]
  [zircaloy_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e6
    poissons_ratio = 0.25
  []
  [ZircDamage]
    type = ZircaloyDamage

    elastic_modulus_hydride = 1e6
    poisson_ratio_hydride = 0.25
    K_hydride = 1e20
    reference_strain_hydride = 5.e-2
    hardening_exponent_hydride = 8.e-2

    K_zircaloy = 1e20
    reference_strain_zircaloy = 5.e-2
    hardening_exponent_zircaloy = 8.e-2

    damage_strain_hydride = 2.e2
    volume_fraction_self = 5.e-2
    volume_fraction_cross = 5.e-2
    damage_strain_zircaloy = 2.e2
    volume_fraction_radial = 8.e-2
    volume_fraction_circumferential = 8.e-2
    clad_inner_radius = 0
    clad_outer_radius = 1
    debug = 1
  []

  [density]
    type = StrainAdjustedDensity
    strain_free_density = 0.283
    outputs = all
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  end_time = 1.0
[]

[Outputs]
  exodus = true
[]

(test/tests/solid_mechanics/zry_creep/creep_sra_thermal_only_damage.i)

#--------------------------------------------------------------------------------
#
# This test is based on the one described below.  With the addition of hydrides
# and damage, the results no longer match those below, though they are close.
#
#--------------------------------------------------------------------------------
#
#   This test case is prepared to test stress recrystalization annealed creep
#   model in ZryCreepUpdate, with only thermal creep activated.
#
#   - Geometry:
#      Ri = 0.005 m
#      Ro = 0.0055 m
#      H  = 0.01 m
#
#   - Single element
#
#   - Temperature = 650 K
#
#   - Boundary conditions:
#       pressure at inner surface = 40  MPa
#       pressure at outer surface = 0.0 MPa
#       displacement in Z constrained at top and bottom (infinitely long tube)
#
#   - Stresses from original solid mechanics version:
#     sigma_rr    =   -18.88 MPa
#     sigma_theta =   402.2 MPa
#     sigma_zz    =   120.4 MPa
#
#  -  Hand calculation results:
#      effective_creep_strain = thermal_creep_rate * dt
#     The Limback equation for primary creep is
#       thermal_creep_rate = A * E / temperature * stress_term^n * exp(-Q/(R*t)) * 1/3600 [1/sec]
#      where the creep rate is calculated using the material properties A, Q, and n for
#      stress relieved annealed zircaloy. The equations for the young's modulus, E,
#      and the stress_term are given in the documentation for the ZryThermalCreepLimbackUpdate
#      class, and the stress_delta term in that equation is determined by RadialReturnStressUpdate
#
#      Given the geometry of the one element mesh, qps 0 and 2 have the same stress delta, and
#      qps 1 and 3 have the same delta.  The calculated effective thermal strain rate listed
#      below is the average of the effective thermal strain rate calculated at these four
#      individual qps.
#
#     Time Step    Analytical Effective Creep Strain    BISON Effective Creep Strain
#        1                3.694e-5                           3.680e-5
#        2                7.373e-5                           7.346e-5
#        3                1.104e-4                           1.100e-4
#        4                1.469e-4                           1.465e-4
#
#--------------------------------------------------------------------------------
[GlobalParams]
  displacements = 'disp_x disp_y'
  temperature = temp
  volumetric_locking_correction = true
[]

[Mesh]
  coord_type = RZ
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 1
    ny = 1
    xmin = 0.005
    xmax = 0.0055
    ymin = 0.0
    ymax = 0.01
    elem_type = QUAD4
  []
  [sbb1]
    type = SubdomainBoundingBoxGenerator
    input = gmg
    block_id = 1
    bottom_left = '0 0 0'
    top_right = '.1 .1 0'
  []
[]

[Variables]
  [temp]
    initial_condition = 650.0
  []
[]

[AuxVariables]
  [fast_neutron_flux]
  []
  [fast_neutron_fluence]
  []
[]

[Functions]
  [pressure_function]
    type = PiecewiseLinear
    x = '0 100000'
    y = '1      1'
  []
[]

[Physics/SolidMechanics/QuasiStatic/clad]
  strain = FINITE
  add_variables = true
  generate_output = 'stress_xx stress_yy stress_zz elastic_strain_xx elastic_strain_yy elastic_strain_zz creep_strain_xx creep_strain_yy creep_strain_zz'
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temp
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temp
  []
[]

[AuxKernels]
  [fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    factor = 1e18 # n/m^2-sec
    execute_on = 'initial linear'
  []
  [fast_neutron_fluence]
    type = FastNeutronFluenceAux
    variable = fast_neutron_fluence
    fast_neutron_flux = fast_neutron_flux
  []
[]

[BCs]
  [Pressure]
    [outer_surface]
      boundary = 'right'
      factor = 0.0
      function = pressure_function
    []
    [inner_surface]
      boundary = 'left'
      factor = 40.0e6
      function = pressure_function
    []
  []
  [u_bottom_fix]
    type = DirichletBC
    variable = disp_y
    boundary = 'top bottom'
    value = 0.0
  []
  [temp_bc_1]
    type = DirichletBC
    variable = temp
    boundary = 'left top right bottom'
    value = 650.0
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1.0e11
    poissons_ratio = 0.3
  []
  [stess]
    type = ZircaloyDamage

    elastic_modulus_hydride = 9.67e10
    poisson_ratio_hydride = 0.389
    K_hydride = 411e6
    reference_strain_hydride = 1e-3
    hardening_exponent_hydride = 0
    damage_strain_hydride = 0.004

    K_zircaloy = 4e8
    reference_strain_zircaloy = 1e-3
    hardening_exponent_zircaloy = 0.05
    damage_strain_zircaloy = 0.035

    volume_fraction_radial = 0.005
    volume_fraction_circumferential = 0.067
    scale_circumferential = true
    clad_inner_radius = 0.005
    clad_outer_radius = 0.0055
    zircaloy_model = 'zry_thermal_creep'
  []
  [zry_thermal_creep]
    type = ZryCreepLimbackHoppeUpdate
    temperature = temp
    fast_neutron_fluence = fast_neutron_fluence
    model_primary_creep = false
    model_irradiation_creep = false
    zircaloy_material_type = stress_relief_annealed
  []
  [clad_density]
    type = StrainAdjustedDensity
    strain_free_density = 6500
  []
  [thermal]
    type = HeatConductionMaterial
    specific_heat = 1.0
    thermal_conductivity = 100.0
  []
[]

[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 = 20
  nl_max_its = 20
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-6
  l_tol = 1e-5
  start_time = 0.0
  end_time = 400
  dt = 100
[]

[Postprocessors]
  [elastic_strain_xx]
    type = ElementAverageValue
    variable = elastic_strain_xx
  []
  [elastic_strain_yy]
    type = ElementAverageValue
    variable = elastic_strain_yy
  []
  [elastic_strain_zz]
    type = ElementAverageValue
    variable = elastic_strain_zz
  []
  [creep_strain_xx]
    type = ElementAverageValue
    variable = creep_strain_xx
  []
  [creep_strain_yy]
    type = ElementAverageValue
    variable = creep_strain_yy
  []
  [creep_strain_zz]
    type = ElementAverageValue
    variable = creep_strain_zz
  []
  [stress_xx]
    type = ElementAverageValue
    variable = stress_xx
  []
  [stress_yy]
    type = ElementAverageValue
    variable = stress_yy
  []
  [stress_zz]
    type = ElementAverageValue
    variable = stress_zz
  []
[]

[Outputs]
  csv = true
[]

(test/tests/zircaloy_damage/radial_crack_plus_y.i)

#
# This test pulls a single element in the z direction until failure.
#   Shear is also introduced.
#
[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

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

[Mesh]
  use_displaced_mesh = false
  [mesh]
    type = GeneratedMeshGenerator
    dim = 3
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [cube]
    add_variables = true
    strain = SMALL
    incremental = true
    extra_vector_tags = 'ref'
    generate_output = 'stress_xx stress_yz stress_yy stress_zz strain_xx strain_yz strain_yy strain_zz'
  []
[]

[BCs]
  [x_fixed]
    type = DirichletBC
    variable = disp_x
    boundary = 'left right'
    value = 0.0
  []
  [bottom_fixed]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  []
  [back_fixed]
    type = DirichletBC
    variable = disp_z
    boundary = back
    value = 0.0
  []

  [front_pull]
    type = FunctionDirichletBC
    variable = disp_z
    boundary = front
    function = '0.1*t'
  []
  [front_shear]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = front
    function = '0.01*t'
  []
[]

[Materials]
  [zircaloy_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 5e10
    poissons_ratio = 0.25
  []
  [ZircDamage]
    type = ZircaloyDamage

    elastic_modulus_hydride = 10e10
    poisson_ratio_hydride = 0.25
    K_hydride = 1e30
    reference_strain_hydride = 1e-3
    hardening_exponent_hydride = 8e-2
    damage_strain_hydride = 3e-2

    K_zircaloy = 5e8
    reference_strain_zircaloy = 1e-3
    hardening_exponent_zircaloy = 0.05
    damage_strain_zircaloy = 100

    volume_fraction_radial = 0.5
    volume_fraction_circumferential = 0.0
    scale_circumferential = false
    clad_inner_radius = 0
    clad_outer_radius = 1
    debug = 0
  []
[]

[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  solve_type = 'NEWTON'
  petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -snes_force_iteration'
  petsc_options_value = 'lu       superlu_dist                  1'

  line_search = 'none'

  l_max_its = 50
  nl_max_its = 50
  nl_rel_tol = 1e-3
  nl_abs_tol = 1e-6
  l_tol = 1e-2
  start_time = 0.0
  end_time = 0.7
  dt = 0.01
  dtmin = 1e-6
[]

[Outputs]
  time_step_interval =  1
  exodus = true
  print_linear_residuals = true
[]

[Postprocessors]
  [stress_xx]
    type = ElementAverageValue
    variable = stress_xx
  []
  [strain_xx]
    type = ElementAverageValue
    variable = strain_xx
  []
  [stress_yy]
    type = ElementAverageValue
    variable = stress_yy
  []
  [strain_yy]
    type = ElementAverageValue
    variable = strain_yy
  []
  [stress_zz]
    type = ElementAverageValue
    variable = stress_zz
  []
  [strain_zz]
    type = ElementAverageValue
    variable = strain_zz
  []
  [stress_yz]
    type = ElementAverageValue
    variable = stress_yz
  []
  [strain_yz]
    type = ElementAverageValue
    variable = strain_yz
  []
[]

(test/tests/zircaloy_damage/plastic.i)

# This is a test of the isotropic power law hardening portion of the model.
#
# In this problem, a single Hex 8 element is fixed at the bottom and pulled at the top
# at a constant rate of 0.1.
#
# Before yield, stress = strain (=0.1*t) as youngs modulus is 1.0.
#
# The reference strain is set to zero for both hydride and Zircaloy.  In this case, the model uses
#   Y = K * e^n
# to compute the yield stress.
# The initial yield stress is
#   Y0 = K * (K/E)^(n/(1-n))
#
# The initial yield stress for this problem is 0.25 (0.5 * (0.5 / 1.0)^(0.5/(1-0.5)))
# as the strength coefficient is 0.5 and the strain rate exponent is 0.5).
#
# Therefore, the material should start yielding at t = 2.5 seconds and then follow
# stress = K * e^n or stress ~ 0.5 * (0.1*t)^0.5.

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

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

[Mesh]
  use_displaced_mesh = false
  [mesh]
    type = GeneratedMeshGenerator
    dim = 3
  []
[]

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

[Physics/SolidMechanics/QuasiStatic]
  [pellets]
    add_variables = true
    strain = SMALL
    incremental = true
    extra_vector_tags = 'ref'
    generate_output = 'stress_xx stress_xy stress_yy stress_zz'
  []
[]

[AuxKernels]
  [total_strain_yy]
    type = RankTwoAux
    rank_two_tensor = total_strain
    variable = total_strain_yy
    index_i = 1
    index_j = 1
  []
[]

[Functions]
  [disp_pull]
    type = ParsedFunction
    expression = 0.1*t
  []
[]

[BCs]
  [x_face]
    type = DirichletBC
    variable = disp_x
    boundary = left
    value = 0.0
  []

  [y_face]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  []

  [y_pull]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top
    function = disp_pull
  []

  [z_face]
    type = DirichletBC
    variable = disp_z
    boundary = back
    value = 0.0
  []
[]

[Materials]
  [zircaloy_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1
    poissons_ratio = 0.3
  []
  [ZircDamage]
    type = ZircaloyDamage

    elastic_modulus_hydride = 1
    poisson_ratio_hydride = 0.3
    K_hydride = 0.5
    reference_strain_hydride = 0
    hardening_exponent_hydride = 0.5

    K_zircaloy = 0.5
    reference_strain_zircaloy = 0
    hardening_exponent_zircaloy = 0.5

    damage_strain_hydride = 2.e2
    volume_fraction_self = 5.e-2
    volume_fraction_cross = 5.e-2
    damage_strain_zircaloy = 2.e2
    volume_fraction_radial = 8.e-2
    volume_fraction_circumferential = 8.e-2
    normal_radial = '0.0 0.0 1.0'
    normal_circumferential = '0.0 1.0 0.0'
    clad_inner_radius = 0
    clad_outer_radius = 1
    debug = 0
  []
[]

[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]

[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
  nl_max_its = 50
  nl_rel_tol = 1e-3
  nl_abs_tol = 1e-6
  l_tol = 1e-2
  start_time = 0.0
  end_time = 5
  dt = 0.25
  dtmin = 1e-6
[]

[Outputs]
  time_step_interval =  1
  exodus = true
  print_linear_residuals = true
[]

[Postprocessors]
  [stress_yy]
    type = ElementAverageValue
    variable = stress_yy
  []
  [strain_yy]
    type = ElementAverageValue
    variable = total_strain_yy
  []
[]

(test/tests/zircaloy_damage/elastic_patch.i)

# Patch Test

# This test is designed to compute constant xx, yy, zz, xy, yz, and zx
#  stress on a set of irregular hexes.  The mesh is composed of one
#  block with seven elements.  The elements form a unit cube with one
#  internal element.  There is a nodeset for each exterior node.

# The cube is displaced by 1e-6 units in x, 2e-6 in y, and 3e-6 in z.
#  The faces are sheared as well (1e-6, 2e-6, and 3e-6 for xy, yz, and
#  zx).  This gives a uniform strain/stress state for all six unique
#  tensor components.

# With Young's modulus at 1e6 and Poisson's ratio at 0, the shear
#  modulus is 5e5 (G=E/2/(1+nu)).  Therefore,
#
#  stress xx = 1e6 * 1e-6 = 1
#  stress yy = 1e6 * 2e-6 = 2
#  stress zz = 1e6 * 3e-6 = 3
#  stress xy = 2 * 5e5 * 1e-6 / 2 = 0.5
#             (2 * G   * gamma_xy / 2 = 2 * G * epsilon_xy)
#  stress yz = 2 * 5e5 * 2e-6 / 2 = 1
#  stress zx = 2 * 5e5 * 3e-6 / 2 = 1.5

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

[Mesh]
  use_displaced_mesh = false
  [mesh]
    type = ExamplePatchMeshGenerator
    dim = 3
  []
[]

[Functions]
  [rampConstant1]
    type = PiecewiseLinear
    x = '0. 1. 2.'
    y = '0. 1. 1.'
    scale_factor = 1e-6
  []
  [rampConstant2]
    type = PiecewiseLinear
    x = '0. 1. 2.'
    y = '0. 1. 1.'
    scale_factor = 2e-6
  []
  [rampConstant3]
    type = PiecewiseLinear
    x = '0. 1. 2.'
    y = '0. 1. 1.'
    scale_factor = 3e-6
  []
  [rampConstant4]
    type = PiecewiseLinear
    x = '0. 1. 2.'
    y = '0. 1. 1.'
    scale_factor = 4e-6
  []
  [rampConstant6]
    type = PiecewiseLinear
    x = '0. 1. 2.'
    y = '0. 1. 1.'
    scale_factor = 6e-6
  []
[]

[Variables]
  [disp_x]
    order = FIRST
    family = LAGRANGE
  []
  [disp_y]
    order = FIRST
    family = LAGRANGE
  []
  [disp_z]
    order = FIRST
    family = LAGRANGE
  []
[]

[AuxVariables]
  [elastic_energy]
    order = CONSTANT
    family = MONOMIAL
  []
  [direction]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [fred]
    strain = SMALL
    incremental = true
    generate_output = 'stress_xx stress_yy stress_zz stress_xy stress_yz stress_zx vonmises_stress hydrostatic_stress firstinv_stress secondinv_stress thirdinv_stress maxshear_stress intensity_stress max_principal_stress mid_principal_stress min_principal_stress'
  []
[]

[AuxKernels]
  [elastic_energy]
    type = ElasticEnergyAux
    variable = elastic_energy
  []
  [direction]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = direction
    scalar_type = direction
    direction = '1 1 1'
  []
[]

[BCs]
  [node1_x]
    type = DirichletBC
    variable = disp_x
    boundary = bottom_front_left
    value = 0.0
  []
  [node1_y]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = bottom_front_left
    function = rampConstant2
  []
  [node1_z]
    type = FunctionDirichletBC
    variable = disp_z
    boundary = bottom_front_left
    function = rampConstant3
  []

  [node2_x]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = bottom_front_right
    function = rampConstant1
  []
  [node2_y]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = bottom_front_right
    function = rampConstant2
  []
  [node2_z]
    type = FunctionDirichletBC
    variable = disp_z
    boundary = bottom_front_right
    function = rampConstant6
  []

  [node3_x]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = bottom_back_right
    function = rampConstant1
  []
  [node3_y]
    type = DirichletBC
    variable = disp_y
    boundary = bottom_back_right
    value = 0.0
  []
  [node3_z]
    type = FunctionDirichletBC
    variable = disp_z
    boundary = bottom_back_right
    function = rampConstant3
  []

  [node4_x]
    type = DirichletBC
    variable = disp_x
    boundary = bottom_back_left
    value = 0.0
  []
  [node4_y]
    type = DirichletBC
    variable = disp_y
    boundary = bottom_back_left
    value = 0.0
  []
  [node4_z]
    type = DirichletBC
    variable = disp_z
    boundary = bottom_back_left
    value = 0.0
  []

  [node5_x]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = top_front_left
    function = rampConstant1
  []
  [node5_y]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top_front_left
    function = rampConstant4
  []
  [node5_z]
    type = FunctionDirichletBC
    variable = disp_z
    boundary = top_front_left
    function = rampConstant3
  []

  [node6_x]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = top_front_right
    function = rampConstant2
  []
  [node6_y]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top_front_right
    function = rampConstant4
  []
  [node6_z]
    type = FunctionDirichletBC
    variable = disp_z
    boundary = top_front_right
    function = rampConstant6
  []

  [node7_x]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = top_back_right
    function = rampConstant2
  []
  [node7_y]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top_back_right
    function = rampConstant2
  []
  [node7_z]
    type = FunctionDirichletBC
    variable = disp_z
    boundary = top_back_right
    function = rampConstant3
  []

  [node8_x]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = top_back_left
    function = rampConstant1
  []
  [node8_y]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top_back_left
    function = rampConstant2
  []
  [node8_z]
    type = DirichletBC
    variable = disp_z
    boundary = top_back_left
    value = 0.0
  []
[]

[Materials]
  [zircaloy_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e6
    poissons_ratio = 0
  []
  [ZircDamage]
    type = ZircaloyDamage

    elastic_modulus_hydride = 1e6
    poisson_ratio_hydride = 0
    K_hydride = 1e20
    reference_strain_hydride = 5.e-2
    hardening_exponent_hydride = 8.e-2
    K_zircaloy = 1e20
    reference_strain_zircaloy = 5.e-2
    hardening_exponent_zircaloy = 8.e-2

    damage_strain_hydride = 2.e2
    volume_fraction_self = 5.e-2
    volume_fraction_cross = 5.e-2
    damage_strain_zircaloy = 2.e2
    volume_fraction_radial = 8.e-2
    volume_fraction_circumferential = 8.e-2
    clad_inner_radius = 0
    clad_outer_radius = 1
    debug = 0
  []
[]

[Executioner]

  type = Transient

  solve_type = 'PJFNK'

  nl_abs_tol = 1e-10

  l_max_its = 20

  start_time = 0.0
  dt = 1.0
  num_steps = 2
  end_time = 2.0
[]

[Outputs]
  exodus = true
[]

(test/tests/zircaloy_damage/homogeneous_elastic_2.i)

#
# This problem is taken from the Abaqus verification manual:
#   "1.5.4 Patch test for axisymmetric elements"
# The stress solution is given as:
#   xx = yy = zz = 2000
#   xy = 400
#
# Since the strain is 1e-3 in all three directions, the new density should be
#   new_density = original_density * V_0 / V
#   new_density = 0.283 / (1 + 1e-3 + 1e-3 + 1e-3) = 0.282153
#
# This test uses a separate Material for the Zircaloy
#

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  coord_type = RZ
  use_displaced_mesh = false
  [mesh]
    type = ExamplePatchMeshGenerator
    dim = 2
    x_offset = 1
  []
[]

[Physics/SolidMechanics/QuasiStatic/All]
  strain = SMALL
  incremental = true
  add_variables = true
  generate_output = 'stress_xx stress_yy stress_zz stress_xy stress_yz stress_zx'
[]

[Kernels]
  [body]
    type = BodyForce
    variable = disp_y
    value = 1
    function = '-400/x'
  []
[]

[BCs]
  [ur]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = 'left right'
    function = '1e-3*x'
  []
  [uz]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = 'left right'
    function = '1e-3*(x+y)'
  []
[]

[Materials]
  [zircaloy_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1e6
    poissons_ratio = 0.25
  []
  [ZircDamage]
    type = ZircaloyDamage

    elastic_modulus_hydride = 1e6
    poisson_ratio_hydride = 0.25
    K_hydride = 1e20
    reference_strain_hydride = 5.e-2
    hardening_exponent_hydride = 8.e-2

    K_zircaloy = 1e20
    reference_strain_zircaloy = 5.e-2
    hardening_exponent_zircaloy = 8.e-2

    damage_strain_hydride = 2.e2
    volume_fraction_self = 5.e-2
    volume_fraction_cross = 5.e-2
    damage_strain_zircaloy = 2.e2
    volume_fraction_radial = 8.e-2
    volume_fraction_circumferential = 8.e-2
    clad_inner_radius = 0
    clad_outer_radius = 1
    debug = 1
    zircaloy_model = flubber
  []
  [flubber]
    type = IsotropicPlasticityStressUpdate
    temperature = 0
    yield_stress = 1e20
    hardening_constant = 0
  []

  [density]
    type = StrainAdjustedDensity
    strain_free_density = 0.283
    outputs = all
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  end_time = 1.0
[]

[Outputs]
  exodus = true
[]

Description
Example Input Syntax
Input Parameters
Input Files
References