CompositeSiCElasticityTensor

Computes the orthotropic elasticity tensor for composite (CVI) SiC-SiC

Description

The CompositeSiCElasticityTensor material model determines the elasticity tensor for composite SiC-SiC. Composite SiC-SiC for cladding application is generally orthotropic: 3D braided architecture with fibers along the axial direction as well as at angles $var element = document.getElementById("moose-equation-998dd69a-246d-4cc3-a190-db1aeb544a8f");katex.render("+\\theta", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-3329ac1e-257f-4947-8411-f34f8ef86e87");katex.render("-\\theta", element, {displayMode:false,throwOnError:false});$ from the axial direction of the cladding. The angle $var element = document.getElementById("moose-equation-ece6f0f1-7846-4b8f-bd03-9ccc9da8759d");katex.render("\\theta", element, {displayMode:false,throwOnError:false});$ is generally in the range 40 $var element = document.getElementById("moose-equation-bff67e7a-7f79-4f1b-a4ab-065bdd0bbb1a");katex.render("^\\circ", element, {displayMode:false,throwOnError:false});$ - 60 $var element = document.getElementById("moose-equation-a33bdb37-777c-42d9-8ec2-bf98e4e94cc4");katex.render("^\\circ", element, {displayMode:false,throwOnError:false});$ .

There are three models available for the user: 1) A general model based on the properties presented in Koyanagi et al. (2017) and assumed properties 2) Model based on the properties from GA (2020) and assumed properties and 3) Model based on the properties from Singh et al. (2019). Either of these models can be invoked by providing the parameter elasticity_model with a value from GENERAL, GA and SINGH. It should be noted that the elastic properties of SiC-SiC composite depends on several factors including the fiber orientation and porosity. These models are based on tests conducted on SiC-SiC samples that had a specific orientation of fibers and porosity. The details of these physical properties can be found in respective references.

Note that for a 3D simulation, this model assumes that the cladding's axis is aligned along the Y-axis and continues to be so during the operation. If the axis bends during operation then the elasticity tensor obtained will be erroneous, with the error increasing with the extent of bending.

The variation of the stiffness ( $var element = document.getElementById("moose-equation-31f75322-7dba-4684-bc52-36a1fd0379b6");katex.render("D", element, {displayMode:false,throwOnError:false});$ ) of the material with irradiation dose is considered in all of these models as:

$var element = document.getElementById("moose-equation-b34f9328-6c0d-4b36-b31b-14eebacc1c9e");katex.render("D_{irradiated} = D_{unirradiated} (1 - 18\\epsilon^s)", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-cafdd2b3-45b6-469d-a9b4-9e67ce0a3e96");katex.render("\\epsilon^s", element, {displayMode:false,throwOnError:false});$ (unitless) is the swelling strain Singh and Wirth (2022). Note that $var element = document.getElementById("moose-equation-e1f76207-f351-4bf0-ac71-ce65ebc269d8");katex.render("\\epsilon^s", element, {displayMode:false,throwOnError:false});$ is linear swelling strain and not the volumetric swelling strain.

Figure 1: Schematic of cladding with Young's moduli along different directions.

The GENERAL model is the default model with following values of elastic constants:

Table 1: Elastic constants for the GENERAL model

Elastic constant	Value	Remarks/Reference
$var element = document.getElementById("moose-equation-1b42cf9f-d1d1-4e6b-b974-2e184fb77a16");katex.render("E_r", element, {displayMode:false,throwOnError:false});$ (GPa)	100	assumed
$var element = document.getElementById("moose-equation-60f22ce5-327a-4b5d-a640-49408f13d6c1");katex.render("E_y", element, {displayMode:false,throwOnError:false});$ (GPa)	302.24 - 1.9148 $var element = document.getElementById("moose-equation-fa2e1baf-b909-4fcb-8642-2aafac883ed8");katex.render("\\beta", element, {displayMode:false,throwOnError:false});$	estimated from Koyanagi et al. (2017)
$var element = document.getElementById("moose-equation-71ea38a7-16fa-4ff3-9d89-c11f55c9f3db");katex.render("E_\\theta", element, {displayMode:false,throwOnError:false});$ (GPa)	160	Koyanagi et al. (2017)
$var element = document.getElementById("moose-equation-dbca7d54-2457-4287-8073-3470ac338b32");katex.render("\\nu_{y\\theta}", element, {displayMode:false,throwOnError:false});$	0.25	assumed based on Nozawa et al. (2012)
$var element = document.getElementById("moose-equation-6377da0c-a9dd-48f8-b2d2-60806dd2f36f");katex.render("\\nu_{\\theta r}", element, {displayMode:false,throwOnError:false});$	0.13	assumed based on Nozawa et al. (2012)
$var element = document.getElementById("moose-equation-1695bbdd-6aa2-43e2-86d6-2eafa3acfbee");katex.render("\\nu_{ry}", element, {displayMode:false,throwOnError:false});$	0.13	assumed based on Nozawa et al. (2012)
$var element = document.getElementById("moose-equation-186e23f5-6a92-4431-b7b9-432660714391");katex.render("G_{y\\theta}", element, {displayMode:false,throwOnError:false});$ (GPa)	110	assumed
$var element = document.getElementById("moose-equation-21f2a14b-d731-4d33-90f9-3667290a8349");katex.render("G_{\\theta r}", element, {displayMode:false,throwOnError:false});$ (GPa)	90	assumed
$var element = document.getElementById("moose-equation-de9f32a5-534b-4a62-b492-5edf9bacb7a1");katex.render("G_{ry}", element, {displayMode:false,throwOnError:false});$ (GPa)	90	assumed

where $var element = document.getElementById("moose-equation-58efe526-b02a-4d5a-9c53-4af4e4408fdf");katex.render("\\beta", element, {displayMode:false,throwOnError:false});$ is the angle between fiber orientation and cladding axis. Note that the value of $var element = document.getElementById("moose-equation-67f64e84-2b4f-4fe4-a426-7f609b0727eb");katex.render("E_y", element, {displayMode:false,throwOnError:false});$ is valid only for $var element = document.getElementById("moose-equation-116786c9-e0d2-4cf4-af69-2e3a13498de2");katex.render("\\beta", element, {displayMode:false,throwOnError:false});$ in the range 30 $var element = document.getElementById("moose-equation-2ce3a1af-9c18-4680-a99c-396a662ec003");katex.render("^\\circ", element, {displayMode:false,throwOnError:false});$ - 55 $var element = document.getElementById("moose-equation-1dbbaa16-7b7e-49b6-8593-9b2f107e970a");katex.render("^\\circ", element, {displayMode:false,throwOnError:false});$ . For unknown elastic constants (Young's and shear moduli) in GENERAL and GA models, values close to those reported in Singh et al. (2019) are assumed.

The GA model has the following elastic properties:

Table 2: Elastic constants for the GA model

Elastic constant	Value	Remarks/Reference
$var element = document.getElementById("moose-equation-d208b69a-ef53-45cb-9f13-48ccf0c67109");katex.render("E_r", element, {displayMode:false,throwOnError:false});$ (GPa)	100	assumed
$var element = document.getElementById("moose-equation-9bb5c65f-450c-40f6-b8ad-38c634b8d631");katex.render("E_y", element, {displayMode:false,throwOnError:false});$ (GPa)	171	GA (2020)
$var element = document.getElementById("moose-equation-dd45be31-eb6e-4832-8fe1-b19a605813a5");katex.render("E_\\theta", element, {displayMode:false,throwOnError:false});$ (GPa)	207	GA (2020)
$var element = document.getElementById("moose-equation-28838d73-c0be-4a10-9f41-43817b866326");katex.render("\\nu_{y\\theta}", element, {displayMode:false,throwOnError:false});$	0.12	GA (2020)
$var element = document.getElementById("moose-equation-397d8c07-d556-4736-ab08-17ee9553ccec");katex.render("\\nu_{\\theta r}", element, {displayMode:false,throwOnError:false});$	0.13	assumed based on Nozawa et al. (2012)
$var element = document.getElementById("moose-equation-1e8ebec7-9422-4de0-80cc-238b665b3e02");katex.render("\\nu_{ry}", element, {displayMode:false,throwOnError:false});$	0.13	assumed based on Nozawa et al. (2012)
$var element = document.getElementById("moose-equation-0c94beb2-0795-4e93-9f06-486b282417e1");katex.render("G_{y\\theta}", element, {displayMode:false,throwOnError:false});$ (GPa)	101	GA (2020)
$var element = document.getElementById("moose-equation-125f5177-144f-43d0-9885-4d37c2dc1897");katex.render("G_{\\theta r}", element, {displayMode:false,throwOnError:false});$ (GPa)	90	assumed
$var element = document.getElementById("moose-equation-25ec8bb3-0027-46a6-92a0-e86fd911a6ed");katex.render("G_{ry}", element, {displayMode:false,throwOnError:false});$ (GPa)	90	assumed

The SINGH model has the following elastic properties:

Table 3: Elastic constancts for the SINGH model

Elastic constant	Value	Remarks/Reference
$var element = document.getElementById("moose-equation-4cbd3e20-1eee-415d-ac87-0a3ba2eb35f6");katex.render("E_r", element, {displayMode:false,throwOnError:false});$ (GPa)	79.9	Singh et al. (2019)
$var element = document.getElementById("moose-equation-aa45f0ac-50e6-4f7d-be85-dc0b643f4414");katex.render("E_y", element, {displayMode:false,throwOnError:false});$ (GPa)	263.9	Singh et al. (2019)
$var element = document.getElementById("moose-equation-5d6b22a6-9ed2-4957-9f79-3992b3373cf5");katex.render("E_\\theta", element, {displayMode:false,throwOnError:false});$ (GPa)	248.8	Singh et al. (2019)
$var element = document.getElementById("moose-equation-14c2820a-b902-469e-8e49-ccbeb0c31626");katex.render("\\nu_{y\\theta}", element, {displayMode:false,throwOnError:false});$	0.170	Singh et al. (2019)
$var element = document.getElementById("moose-equation-beba7ebd-a58f-4960-8a33-7f73fd995ddf");katex.render("\\nu_{\\theta r}", element, {displayMode:false,throwOnError:false});$	0.191	Singh et al. (2019)
$var element = document.getElementById("moose-equation-8dadc58d-4890-42a5-819e-3a6e4e2da40c");katex.render("\\nu_{ry}", element, {displayMode:false,throwOnError:false});$	0.204	Singh et al. (2019)
$var element = document.getElementById("moose-equation-defa8036-49cf-47cb-a012-9947a52fcd1a");katex.render("G_{\\theta y}", element, {displayMode:false,throwOnError:false});$ (GPa)	109.5	Singh et al. (2019)
$var element = document.getElementById("moose-equation-38b0557f-0a1b-4b4e-b9b4-5056266e8b79");katex.render("G_{r\\theta}", element, {displayMode:false,throwOnError:false});$ (GPa)	70.5	Singh et al. (2019)
$var element = document.getElementById("moose-equation-11fbc98e-2574-49c2-b9b6-7701d2c741a0");katex.render("G_{ry}", element, {displayMode:false,throwOnError:false});$ (GPa)	85.7	Singh et al. (2019)

Example Input Syntax

[Materials<<<{"href": "../../../syntax/Materials/index.html"}>>>]
  [composite_elasticity_tensor]
    type = CompositeSiCElasticityTensor<<<{"description": "Computes the orthotropic elasticity tensor for composite (CVI) SiC-SiC", "href": "CompositeSiCElasticityTensor.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = clad
    elasticity_model<<<{"description": "Composite SiC elasticity model"}>>> = GA
  []
[]

Input 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
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
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.
elasticity_modelGENERALComposite SiC elasticity model
Default:GENERAL
C++ Type:MooseEnum
Options:GENERAL, GA, SINGH
Controllable:No
Description:Composite SiC elasticity model
elasticity_tensor_prefactorOptional function to use as a scalar prefactor on the elasticity tensor.
C++ Type:FunctionName
Unit:(no unit assumed)
Controllable:No
Description:Optional function to use as a scalar prefactor on the elasticity tensor.
fiber_winding_angle45The angle at which the fibers in the composite are wound relative to the axial direction.
Default:45
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The angle at which the fibers in the composite are wound relative to the axial direction.
swellingswellingIrradiation induced swelling strain in SiC composite
Default:swelling
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Irradiation induced swelling strain in SiC composite

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
use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

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/solid_mechanics/compositeSiC_mechanics/axial_3D.i)
(test/tests/solid_mechanics/compositeSiC_mechanics/axial_2D_ad.i)
(test/tests/solid_mechanics/compositeSiC_mechanics/hoop_2D.i)
(test/tests/solid_mechanics/compositeSiC_mechanics/axial_2D.i)
(test/tests/solid_mechanics/compositeSiC_mechanics/swelling_dependence_2D.i)
(assessment/LWR/validation/ATF_SiC_HFIR/analysis/base.i)
(test/tests/solid_mechanics/compositeSiC_mechanics/hoop_3D.i)
(test/tests/solid_mechanics/compositeSiC_mechanics/hoop_2D_ad.i)
(examples/accident_tolerant_fuel/u3si2_sic/u3si2_outer_monolith_1.5D.i)

References

GA. Material properties manual: general atomics silicon carbide cladding. Technical Report GA-A28712X, General Atomics Inc., 2020.

@TECHREPORT{ga-sic-manual,
    author = "GA",
    title = "Material Properties Manual: General Atomics Silicon Carbide Cladding",
    institution = "General Atomics Inc.",
    year = "2020",
    number = "GA-A28712X"
}

T. Koyanagi, Y. Katoh, G. Singh, and M. Snead. SiC/SiC cladding materials properties handbook. Technical Report ORNL/TM-2017/385, Oak Ridge National Laboratory, 2017.

@TECHREPORT{koyanagi2017,
    author = "Koyanagi, T. and Katoh, Y. and Singh, G. and Snead, M.",
    title = "{SiC/SiC} Cladding Materials Properties Handbook",
    year = "2017",
    number = "ORNL/TM-2017/385",
    institution = "Oak Ridge National Laboratory"
}

Takashi Nozawa, Kazumi Ozawa, Yong-Bum Choi, Akira Kohyama, and Hiroyasu Tanigawa. Determination and prediction of axial/off-axial mechanical properties of SiC/SiC composites. Fusion Engineering and Design, 87(5-6):803–807, 2012. doi:10.1016/j.fusengdes.2012.02.026.

@article{NOZAWA2012803,
    author = "Nozawa, Takashi and Ozawa, Kazumi and Choi, Yong-Bum and Kohyama, Akira and Tanigawa, Hiroyasu",
    title = "Determination and prediction of axial/off-axial mechanical properties of {SiC/SiC} composites",
    journal = "Fusion Engineering and Design",
    volume = "87",
    number = "5-6",
    pages = "803--807",
    year = "2012",
    doi = "10.1016/j.fusengdes.2012.02.026",
    publisher = "Elsevier"
}

G. Singh, T. Koyanagi, C. Petrie, C. Deck, K. Terrani, J.D. Arregui-Mena, and Y. Katoh. Elastic moduli reduction in SiC-SiC tubular specimen after high heat flux neutron irradiation measured by resonant ultrasound spectroscopy. Journal of Nuclear Materials, 523:391–401, 2019. doi:10.1016/j.jnucmat.2019.06.026.

@article{SINGH2019391,
    author = "Singh, G. and Koyanagi, T. and Petrie, C. and Deck, C. and Terrani, K. and Arregui-Mena, J.D. and Katoh, Y.",
    title = "Elastic moduli reduction in {SiC-SiC} tubular specimen after high heat flux neutron irradiation measured by resonant ultrasound spectroscopy",
    journal = "Journal of Nuclear Materials",
    volume = "523",
    pages = "391--401",
    doi = "10.1016/j.jnucmat.2019.06.026",
    year = "2019"
}

Gyanender Singh and Brian D Wirth. Impact of circumferential variation in power, neutron flux and spacer grids on structural behavior of sic-sic cladding. Journal of Nuclear Materials, 565:153726, 2022. doi:10.1016/j.jnucmat.2022.153726.

@article{SINGH2022565,
    author = "Singh, Gyanender and Wirth, Brian D",
    title = "Impact of circumferential variation in power, neutron flux and spacer grids on structural behavior of SiC-SiC cladding",
    journal = "Journal of Nuclear Materials",
    volume = "565",
    pages = "153726",
    year = "2022",
    doi = "10.1016/j.jnucmat.2022.153726",
    publisher = "Elsevier"
}

(test/tests/solid_mechanics/compositeSiC_mechanics/hoop_3D.i)

# In a tube (ri = 1 m, ro = 1.03 m, h = 2 m) an internal pressure (10 MPa) is
# applied. The hoop strain calculations are shown below:
#
# Subscript r: radial    y: axial    t: hoop
#
# epsilon_tt = term1 + term2 + term3
#              term1 = -nu_rt * sigma_rr / E_r
#              term2 = -nu_yt * sigma_yy / E_y
#              term3 = sigma_tt / E_t
#
# nu_rt = nu_tr * E_r / E_t = 0.0628
# nu_yr = nu_ry * E_y / E_r = 0.2223
#
# epsilon_tt = term1 + term2 + term3
#            = -6.2801e-6 + 0.0 + 1.6344e-3
#            = 1.6282e-3
#
# The results of BISON compared to the analytical solution:
#
#          BISON               Analytical
#       strain_hoop (m/m)    strain_hoop (m/m)
#        1.6220e-03           1.6282e-3
#
# The minor difference between analytical solution and BISON prediction
# is due to a) thin cylinder approximation in the analytical calculation
# (calculation of hoop stress sigma_tt) and b) coarse mesh used in the
# simulation.
#
[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

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

[Mesh]
  coord_type = XYZ
  [mesh]
    type = FileMeshGenerator
    file = quarter_cylinder_coarse.e
  []
[]

[Functions]
  [swelling_func]
    type = ParsedFunction
    expression = 0.0
  []
[]

[AuxVariables]
  [hoop_stress]
    block = clad
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [hoop_stress]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = hoop_stress
    scalar_type = HoopStress
    execute_on = timestep_end
  []
[]

[BCs]
  [x_symm]
    type = DirichletBC
    variable = disp_x
    boundary = x_plane
    value = 0.0
  []
  [z_symm]
    type = DirichletBC
    variable = disp_z
    boundary = z_plane
    value = 0.0
  []
  [no_y]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  []
  [Pressure]
    [innerPressure]
      boundary = inner
      function = 10e6
    []
  []
[]

[Materials]
  [composite_elasticity_tensor]
    type = CompositeSiCElasticityTensor
    block = clad
    elasticity_model = GA
  []
  [composite_stress]
    type = ComputeStrainIncrementBasedStress
    block = clad
  []
  [irradiation_swelling]
    type = GenericFunctionMaterial
    block = clad
    prop_names = swelling
    prop_values = swelling_func
    outputs = all
  []
[]

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

[Executioner]
  solve_type = Newton
  type = Steady
  l_tol = 1E-5
  nl_abs_tol = 1E-3
  nl_rel_tol = 1E-8
  l_max_its = 200
  nl_max_its = 400

  petsc_options_iname = '-pc_type -pc_asm_overlap -sub_pc_type -ksp_type -ksp_gmres_restart'
  petsc_options_value = ' asm      2              lu            gmres     200'
[]

[Postprocessors]
  [hoop_stress]
    type = ElementalVariableValue
    variable = hoop_stress
    elementid = 231
  []
  [hoop_strain]
    type = ElementalVariableValue
    variable = strain_zz
    elementid = 231
  []
[]

[Outputs]
  csv = true
  color = false
[]

(test/tests/solid_mechanics/compositeSiC_mechanics/axial_3D.i)

# In a tube (ri = 1 m, ro = 1.03 m, h = 2 m) pressure (-10 MPa) is applied at
# its cross-section. The axial strain calculations are shown below:
#
# Subscript r: radial    y: axial    t: hoop
#
# epsilon_yy = term1 + term2 + term3
#              term1 = -nu_ry * sigma_rr / E_r
#              term2 = sigma_yy / E_y
#              term3 = -nu_ty * sigma_tt / E_t
##
# epsilon_tt = term1 + term2 + term3
#            = 0.0 + 3.7893e-5 + 0
#            = 3.7893e-5
#
# The results of BISON compared to the analytical solution:
#
#          BISON               Analytical
#       strain_hoop (m/m)    strain_hoop (m/m)
#        3.7883e-5           3.7893-5
#

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

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

[Mesh]
  coord_type = XYZ
  [mesh]
    type = FileMeshGenerator
    file = quarter_cylinder_coarse.e
  []
[]

[Functions]
  [swelling_func]
    type = ParsedFunction
    expression = 0.0
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = clad
    add_variables = true
    strain = FINITE
    generate_output = 'stress_yy strain_yy'
    extra_vector_tags = 'ref'
  []
[]

[BCs]
  [x_symm]
    type = DirichletBC
    variable = disp_x
    boundary = x_plane
    value = 0.0
  []
  [z_symm]
    type = DirichletBC
    variable = disp_z
    boundary = z_plane
    value = 0.0
  []
  [no_y]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  []
  [Pressure]
    [topPressure]
      boundary = top
      function = -10e6
    []
  []
[]

[Materials]
  [composite_elasticity_tensor]
    type = CompositeSiCElasticityTensor
    block = clad
    elasticity_model = SINGH
    swelling = swelling
  []
  [composite_stress]
    type = ComputeStrainIncrementBasedStress
    block = clad
  []
  [irradiation_swelling]
    type = GenericFunctionMaterial
    block = clad
    prop_names = swelling
    prop_values = swelling_func
    outputs = all
  []
[]

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

[Executioner]
  solve_type = Newton
  type = Steady
  l_tol = 1E-5
  nl_abs_tol = 1E-3
  nl_rel_tol = 1E-8
  l_max_its = 200
  nl_max_its = 400

  petsc_options_iname = '-pc_type -pc_asm_overlap -sub_pc_type -ksp_type -ksp_gmres_restart'
  petsc_options_value = ' asm      2              lu            gmres     200'
[]

[Postprocessors]
  [axial_strain]
    type = ElementalVariableValue
    variable = strain_yy
    elementid = 231
    execute_on = 'timestep_end'
  []
[]

[Outputs]
  csv = true
  color = false
[]

(test/tests/solid_mechanics/compositeSiC_mechanics/axial_2D_ad.i)

# This test case is prepared to test the elastic moduli of composite (CVI)
# SiC/SiC using the correlation provided in:
#
# L. L. Snead, T. Nozawa, Y. Katoh, T-S. Byun, S. Kondo, D. A. Petti,
# "Handbook of SiC propeties for fuel performance modeling," Journal of Nuclear
# Materials, 371, 2007, p. 329-377.
#
# The geometry represents a ring of monolothic SiC with inner radius of 0.005 m,
# outer radius of 0.0055 m, and a height of 0.01 m.
#
# A pressure (-100 MPa) is applied to the top of the ring (longitudinal cross-section).
# The strain calculation in the radial direction is shown below.
#
# Subscript r: radial    y: axial    t: circumferential
#
# epsilon_rr = term1 + term2 + term3
#              term1 = sigma_rr / E_r
#              term2 = -nu_yr * sigma_yy / E_y
#              term3 = -nu_tr * sigma_tt / E_t
#
# term1 and term3 are very small for axial loading as sigma_rr and sigma_tt are
# small and can be neglected.
#
# epsilon_rr = -nu_yr * sigma_yy / E_y
#            = -(nu_ry * E_y / E_r) * sigma_yy / E_y
#            = -nu_ry * sigma_yy / E_r
#
# epsilon_rr = 0.13 * 100e6 / 100e9 = 0.00013
#
#
# epsilon_yy = term1 + term2 + term3
#              term1 = -nu_ry * sigma_rr / E_r
#              term2 = sigma_yy / E_y
#              term3 = -nu_ty * sigma_tt / E_t
#
# term1 and term3 are very small for axial loading as sigma_rr and sigma_tt are
# small and can be neglected.
#
# epsilon_yy = sigma_yy / E_y
#            = -100e6 / 216.074e9
#            = -4.6280e-4
#
#
# The results of BISON compared to the analytical solution:
#
#          BISON             Analytical          BISON             Analytical
#       strain_rr (m/m)    strain_rr (m/m)    strain_yy (m/m)     strain_yy (m/m)
#         0.00013             0.00013          -4.6280e-4          -4.6280e-4

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  coord_type = RZ
  use_displaced_mesh = false
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0.005
    xmax = 0.0055
    ymin = 0
    ymax = 0.01
  []
[]

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = SMALL
        block = 0
        add_variables = true
        generate_output = 'stress_xx stress_yy stress_zz elastic_strain_xx elastic_strain_yy elastic_strain_zz'
        use_automatic_differentiation = true
      []
    []
  []
[]

[Functions]
  [pressure_function]
    type = PiecewiseLinear
    x = '0  1'
    y = '1  1'
  []
  [swelling_func]
    type = ParsedFunction
    expression = 0.0
  []
[]

[BCs]
  [top_pressure] # Simplified BC to apply pressure in only disp_y
    type = ADPressure
    boundary = top
    variable = disp_y
    function = pressure_function
    factor = 100e6
  []

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

[Materials]
  [elasticity_tensor]
    type = ADCompositeSiCElasticityTensor
  []
  [stress]
    type = ADComputeLinearElasticStress
  []
  [clad_density]
    type = ADStrainAdjustedDensity
    strain_free_density = 2700.0
  []
  [irradiation_swelling]
    type = ADGenericFunctionMaterial
    prop_names = swelling
    prop_values = swelling_func
    outputs = all
  []
[]

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

[Executioner]
  type = Transient

  solve_type = Newton

  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-8
  nl_abs_tol = 1e-8
  l_tol = 1e-8
  start_time = 0.0
  end_time = 1
  dt = 1
[]

[Postprocessors]
  [elastic_strain_rr]
    type = ElementAverageValue
    variable = elastic_strain_xx
  []
  [elastic_strain_yy]
    type = ElementAverageValue
    variable = elastic_strain_yy
  []
[]

[Outputs]
  csv = true
[]

(test/tests/solid_mechanics/compositeSiC_mechanics/hoop_2D.i)

# This test case is prepared to test the elastic moduli of composite (CVI)
# SiC/SiC using the correlation provided in:
#
# L. L. Snead, T. Nozawa, Y. Katoh, T-S. Byun, S. Kondo, D. A. Petti,
# "Handbook of SiC propeties for fuel performance modeling," Journal of Nuclear
# Materials, 371, 2007, p. 329-377.
#
# The geometry represents a ring of monolothic SiC with inner radius of 0.005 m,
# outer radius of 0.0055 m, and a height of 0.01 m.
#
# A radial pressure (10 MPa) is applied to the inner surface of the ring.
# The strain calculations in the hoop direction are shown below.
#
# Subscript r: radial    y: axial    t: hoop
#
# epsilon_tt = term1 + term2 + term3
#              term1 = -nu_rt * sigma_rr / E_r
#              term2 = -nu_yt * sigma_yy / E_y
#              term3 = sigma_tt / E_t
#
# nu_rt = nu_tr * E_r / E_t = 0.08125
# nu_yr = nu_ry * E_y / E_r = 0.25
#
# epsilon_tt = term1 + term2 + term3
#            = -8.125e-6 + 0.0 + 0.00065625
#            = 0.000648125
#
# The results of BISON compared to the analytical solution:
#
#          BISON               Analytical
#       strain_hoop (m/m)    strain_hoop (m/m)
#        6.28499e-04           6.48125e-04
#
# The small difference is due to using approximation of thin cylinder for
# sigma_tt calculation

[GlobalParams]
  displacements = 'disp_x disp_y'
  volumetric_locking_correction = true
[]

[Mesh]
  coord_type = RZ
  use_displaced_mesh = false
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0.005
    xmax = 0.0055
    ymin = 0
    ymax = 0.01
  []
[]

[Functions]
  [pressure_function]
    type = PiecewiseLinear
    x = '0  1'
    y = '1  1'
  []
  [swelling_func]
    type = ParsedFunction
    expression = 0.0
  []
[]

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = SMALL
        block = 0
        add_variables = true
        generate_output = 'elastic_strain_xx elastic_strain_yy elastic_strain_zz stress_xx stress_yy stress_zz'
      []
    []
  []
[]

[BCs]
  [Pressure]
    [inner_pressure]
      boundary = left
      function = pressure_function
      factor = 10e6
    []
  []

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

[Materials]
  [elasticity_tensor]
    type = CompositeSiCElasticityTensor
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [clad_density]
    type = StrainAdjustedDensity
    strain_free_density = 2700.0
  []
  [irradiation_swelling]
    type = GenericFunctionMaterial
    prop_names = swelling
    prop_values = swelling_func
    outputs = all
  []
[]

[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-8
  nl_abs_tol = 1e-8
  l_tol = 1e-8
  start_time = 0.0
  end_time = 1
  dt = 1
[]

[Postprocessors]
  [elastic_strain_hoop]
    type = ElementAverageValue
    variable = elastic_strain_zz
  []
[]

[Outputs]
  csv = true
[]

(test/tests/solid_mechanics/compositeSiC_mechanics/axial_2D.i)

# This test case is prepared to test the elastic moduli of composite (CVI)
# SiC/SiC using the correlation provided in:
#
# L. L. Snead, T. Nozawa, Y. Katoh, T-S. Byun, S. Kondo, D. A. Petti,
# "Handbook of SiC propeties for fuel performance modeling," Journal of Nuclear
# Materials, 371, 2007, p. 329-377.
#
# The geometry represents a ring of monolothic SiC with inner radius of 0.005 m,
# outer radius of 0.0055 m, and a height of 0.01 m.
#
# A pressure (-100 MPa) is applied to the top of the ring (longitudinal cross-section).
# The strain calculation in the radial direction is shown below.
#
# Subscript r: radial    y: axial    t: circumferential
#
# epsilon_rr = term1 + term2 + term3
#              term1 = sigma_rr / E_r
#              term2 = -nu_yr * sigma_yy / E_y
#              term3 = -nu_tr * sigma_tt / E_t
#
# term1 and term3 are very small for axial loading as sigma_rr and sigma_tt are
# small and can be neglected.
#
# epsilon_rr = -nu_yr * sigma_yy / E_y
#            = -(nu_ry * E_y / E_r) * sigma_yy / E_y
#            = -nu_ry * sigma_yy / E_r
#
# epsilon_rr = 0.13 * 100e6 / 100e9 = 0.00013
#
#
# epsilon_yy = term1 + term2 + term3
#              term1 = -nu_ry * sigma_rr / E_r
#              term2 = sigma_yy / E_y
#              term3 = -nu_ty * sigma_tt / E_t
#
# term1 and term3 are very small for axial loading as sigma_rr and sigma_tt are
# small and can be neglected.
#
# epsilon_yy = sigma_yy / E_y
#            = -100e6 / 216.074e9
#            = -4.6280e-4
#
#
# The results of BISON compared to the analytical solution:
#
#          BISON             Analytical          BISON             Analytical
#       strain_rr (m/m)    strain_rr (m/m)    strain_yy (m/m)     strain_yy (m/m)
#         0.00013             0.00013          -4.6280e-4          -4.6280e-4

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  coord_type = RZ
  use_displaced_mesh = false
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0.005
    xmax = 0.0055
    ymin = 0
    ymax = 0.01
  []
[]

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = SMALL
        block = 0
        add_variables = true
        generate_output = 'stress_xx stress_yy stress_zz elastic_strain_xx elastic_strain_yy elastic_strain_zz'
      []
    []
  []
[]

[Functions]
  [pressure_function]
    type = PiecewiseLinear
    x = '0  1'
    y = '1  1'
  []
  [swelling_func]
    type = ParsedFunction
    expression = 0.0
  []
[]

[BCs]
  [top_pressure] # Simplified BC to apply pressure in only disp_y
    type = Pressure
    boundary = top
    variable = disp_y
    function = pressure_function
    factor = 100e6
  []

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

[Materials]
  [elasticity_tensor]
    type = CompositeSiCElasticityTensor
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [clad_density]
    type = StrainAdjustedDensity
    strain_free_density = 2700.0
  []
  [irradiation_swelling]
    type = GenericFunctionMaterial
    prop_names = swelling
    prop_values = swelling_func
    outputs = all
  []
[]

[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-8
  nl_abs_tol = 1e-8
  l_tol = 1e-8
  start_time = 0.0
  end_time = 1
  dt = 1
[]

[Postprocessors]
  [elastic_strain_rr]
    type = ElementAverageValue
    variable = elastic_strain_xx
  []
  [elastic_strain_yy]
    type = ElementAverageValue
    variable = elastic_strain_yy
  []
[]

[Outputs]
  csv = true
[]

(test/tests/solid_mechanics/compositeSiC_mechanics/swelling_dependence_2D.i)

# This test case is prepared to test the elastic moduli of composite (CVI)
# SiC/SiC using the correlation provided in:
#
# L. L. Snead, T. Nozawa, Y. Katoh, T-S. Byun, S. Kondo, D. A. Petti,
# "Handbook of SiC propeties for fuel performance modeling," Journal of Nuclear
# Materials, 371, 2007, p. 329-377.
#
# The geometry represents a ring of monolothic SiC with inner radius of 0.005 m,
# outer radius of 0.0055 m, and a height of 0.01 m.
#
# A pressure (-100 MPa) is applied to the top of the ring (longitudinal cross-section).
# The swelling strain is considered to be 0.067 (~2% volumetric strain). The strain
# calculation in the radial direction is shown below.
#
# Subscript r: radial    y: axial    t: circumferential
#
# epsilon_rr = term1 + term2 + term3
#              term1 = sigma_rr / E_r
#              term2 = -nu_yr * sigma_yy / E_y
#              term3 = -nu_tr * sigma_tt / E_t
#
# term1 and term3 are very small for axial loading as sigma_rr and sigma_tt are
# small and can be neglected.
#
# epsilon_rr = -nu_yr * sigma_yy / E_y
#            = -(nu_ry * E_y / E_r) * sigma_yy / E_y
#            = -nu_ry * sigma_yy / E_r
#
# epsilon_rr = 0.13 * 100e6 / (100e9 * (1 - 18.0 *0.0067))= 1.4782e-4
#
#
# epsilon_yy = term1 + term2 + term3
#              term1 = -nu_ry * sigma_rr / E_r
#              term2 = sigma_yy / E_y
#              term3 = -nu_ty * sigma_tt / E_t
#
# term1 and term3 are very small for axial loading as sigma_rr and sigma_tt are
# small and can be neglected.
#
# epsilon_yy = sigma_yy / E_y
#            = -100e6 / (216.074e9 * (1 - 18.0 * 0.0067))
#            = -5.2627e-4
#
# The results of BISON compared to the analytical solution:
#
#          BISON             Analytical          BISON             Analytical
#      strain_rr (m/m)    strain_rr (m/m)    strain_yy (m/m)     strain_yy (m/m)
#         1.4780e-4          1.4782e-4         -5.2627e-04         -5.2627e-4

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  coord_type = RZ
  use_displaced_mesh = false
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0.005
    xmax = 0.0055
    ymin = 0
    ymax = 0.01
  []
[]

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = SMALL
        block = 0
        add_variables = true
        generate_output = 'stress_xx stress_yy stress_zz elastic_strain_xx elastic_strain_yy elastic_strain_zz'
      []
    []
  []
[]

[Functions]
  [pressure_function]
    type = PiecewiseLinear
    x = '0  1'
    y = '1  1'
  []
  [swelling_func]
    type = ParsedFunction
    expression = 0.0067
  []
[]

[BCs]
  [top_pressure] # Simplified BC to apply pressure in only disp_y
    type = Pressure
    boundary = top
    variable = disp_y
    function = pressure_function
    factor = 100e6
  []

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

[Materials]
  [elasticity_tensor]
    type = CompositeSiCElasticityTensor
  []
  [stress]
    type = ComputeLinearElasticStress
  []
  [clad_density]
    type = StrainAdjustedDensity
    strain_free_density = 2700.0
  []
  [irradiation_swelling]
    type = GenericFunctionMaterial
    prop_names = swelling
    prop_values = swelling_func
    outputs = all
  []
[]

[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-8
  nl_abs_tol = 1e-8
  l_tol = 1e-8
  start_time = 0.0
  end_time = 1
  dt = 1
[]

[Postprocessors]
  [elastic_strain_rr]
    type = ElementAverageValue
    variable = elastic_strain_xx
  []
  [elastic_strain_yy]
    type = ElementAverageValue
    variable = elastic_strain_yy
  []
[]

[Outputs]
  csv = true
[]

(assessment/LWR/validation/ATF_SiC_HFIR/analysis/base.i)

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

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

[Mesh]
  coord_type = RZ
  [clad_mesh]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 100
    xmin = 3.51e-3
    xmax = 4.27e-3
    ymin = 0.0
    ymax = 16.04e-3
  []
[]

[Variables]
  [temperature]
    initial_condition = 295.0 # set initial temperature to coolant inlet
  []
[]

[AuxVariables]
  [fast_neutron_flux]
    order = CONSTANT
    family = MONOMIAL
  []
  [fast_neutron_fluence]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 1e4 2.16e6 2.3328e6'
    y = '0 333333.0 333333.0 0'
    scale_factor = 1
  []
  [axial_peaking_factors]
    type = ParsedFunction
    expression = 1
  []
  [pressure_ramp]
    type = PiecewiseLinear
    x = '-200 0 1e8'
    y = '6.537e-3 1 1'
  []
  [max_time_step_size_func]
    type = ParsedFunction
    expression = 'if(t < 2.16e6, 5e4, 1e4)'
  []
  [heat_flux_func]
    type = ParsedFunction
    expression = 'if(t < 2.16e6, 6e5, 0)'
  []
[]

[Constraints]
  [top]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    secondary = 'top'
    penalty = 1e+11
    formulation = penalty
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    add_variables = true
    strain = FINITE
    eigenstrain_names = 'composite_thermal_strain composite_swelling_strain'
    automatic_eigenstrain_names = true
    generate_output = 'hoop_stress hoop_strain stress_xx stress_yy stress_zz vonmises_stress strain_zz strain_yy'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    extra_vector_tags = 'ref'
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temperature
    extra_vector_tags = 'ref'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temperature
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [fast_neutron_flux]
    type = MaterialRealAux
    variable = fast_neutron_flux
    property = fast_neutron_flux
    execute_on = timestep_begin
  []
  [fast_neutron_fluence]
    type = MaterialRealAux
    variable = fast_neutron_fluence
    property = fast_neutron_fluence
    execute_on = timestep_begin
  []
[]

[BCs]
  [no_y]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  []
  [heat_flux_clad_inner_surface]
    type = FunctionNeumannBC
    boundary = left
    variable = temperature
    function = heat_flux_func
  []
  [right]
    type = ConvectiveHeatFluxBC
    variable = temperature
    boundary = 'right'
    T_infinity = 295.0
    heat_transfer_coefficient = 12000
    heat_transfer_coefficient_dT = 0
  []
[]

[Materials]
  [flux]
    type = FastNeutronFlux
    calculate_fluence = true
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    factor = 3e13
  []
  [composite_thermal]
    type = CompositeSiCThermal
    thermal_conductivity_model = STONE
    temperature = temperature
    outputs = 'all'
  []
  [composite_density]
    type = StrainAdjustedDensity
    strain_free_density = 2700.0
  []
  [composite_elasticity_tensor]
    type = CompositeSiCElasticityTensor
    elasticity_model = SINGH
  []
  [composite_stress]
    type = ComputeStrainIncrementBasedStress
  []
  [composite_thermal_expansion]
    type = CompositeSiCThermalExpansionEigenstrain
    stress_free_temperature = 295.0
    temperature = temperature
    eigenstrain_name = composite_thermal_strain
    outputs = 'all'
  []
  [composite_irradiation_swelling]
    type = CompositeSiCVolumetricSwellingEigenstrain
    temperature = temperature
    fast_neutron_fluence = fast_neutron_fluence
    swelling_model = KATOH
    number_of_substeps = 1000
    eigenstrain_name = composite_swelling_strain
    outputs = 'all'
  []
[]

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

[PerformanceMetricOutputs]
[]

[Executioner]
  type = Transient
  solve_type = 'NEWTON' #PJFNK'

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

  line_search = 'none'

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

  start_time = -200
  n_startup_steps = 1
  end_time = 2.3328e6
  dtmin = 1

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 2e2
    optimal_iterations = 25
    iteration_window = 5
    growth_factor = 2
    cutback_factor = .5
    timestep_limiting_postprocessor = max_time_step_size_pp
  []
[]

[Postprocessors]
  [hoop_stress_outer]
    type = ElementalVariableValue
    variable = hoop_stress
    elementid = 499
  []
  [axial_stress_outer]
    type = ElementalVariableValue
    variable = stress_yy
    elementid = 499
  []
  [hoop_stress_inner]
    type = ElementalVariableValue
    variable = hoop_stress
    elementid = 490
  []
  [axial_stress_inner]
    type = ElementalVariableValue
    variable = stress_yy
    elementid = 490
  []
  [hoop_strain_inner]
    type = ElementalVariableValue
    variable = strain_zz
    elementid = 490
  []
  [axial_strain_inner]
    type = ElementalVariableValue
    variable = strain_yy
    elementid = 490
  []
  [hoop_strain_outer]
    type = ElementalVariableValue
    variable = strain_zz
    elementid = 499
  []
  [axial_strain_outer]
    type = ElementalVariableValue
    variable = strain_yy
    elementid = 499
  []
  [max_time_step_size_pp]
    type = FunctionValuePostprocessor
    function = max_time_step_size_func
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [outlet_heat_flux]
    type = SideDiffusiveFluxIntegral
    variable = temperature
    boundary = left
    diffusivity = thermal_conductivity
  []
  [swelling_strain_outer]
    type = ElementalVariableValue
    variable = composite_swelling_strain_22
    elementid = 499
  []
  [swelling_strain_inner]
    type = ElementalVariableValue
    variable = composite_swelling_strain_22
    elementid = 490
  []
  [thermal_strain_outer]
    type = ElementalVariableValue
    variable = composite_thermal_strain_22
    elementid = 499
  []
  [thermal_strain_inner]
    type = ElementalVariableValue
    variable = composite_thermal_strain_22
    elementid = 490
  []
[]

[Outputs]
  perf_graph = true
  exodus = true
  [csv]
    type = CSV
    execute_on = final
  []
  color = false
[]

(test/tests/solid_mechanics/compositeSiC_mechanics/hoop_3D.i)

# In a tube (ri = 1 m, ro = 1.03 m, h = 2 m) an internal pressure (10 MPa) is
# applied. The hoop strain calculations are shown below:
#
# Subscript r: radial    y: axial    t: hoop
#
# epsilon_tt = term1 + term2 + term3
#              term1 = -nu_rt * sigma_rr / E_r
#              term2 = -nu_yt * sigma_yy / E_y
#              term3 = sigma_tt / E_t
#
# nu_rt = nu_tr * E_r / E_t = 0.0628
# nu_yr = nu_ry * E_y / E_r = 0.2223
#
# epsilon_tt = term1 + term2 + term3
#            = -6.2801e-6 + 0.0 + 1.6344e-3
#            = 1.6282e-3
#
# The results of BISON compared to the analytical solution:
#
#          BISON               Analytical
#       strain_hoop (m/m)    strain_hoop (m/m)
#        1.6220e-03           1.6282e-3
#
# The minor difference between analytical solution and BISON prediction
# is due to a) thin cylinder approximation in the analytical calculation
# (calculation of hoop stress sigma_tt) and b) coarse mesh used in the
# simulation.
#
[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

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

[Mesh]
  coord_type = XYZ
  [mesh]
    type = FileMeshGenerator
    file = quarter_cylinder_coarse.e
  []
[]

[Functions]
  [swelling_func]
    type = ParsedFunction
    expression = 0.0
  []
[]

[AuxVariables]
  [hoop_stress]
    block = clad
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [hoop_stress]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = hoop_stress
    scalar_type = HoopStress
    execute_on = timestep_end
  []
[]

[BCs]
  [x_symm]
    type = DirichletBC
    variable = disp_x
    boundary = x_plane
    value = 0.0
  []
  [z_symm]
    type = DirichletBC
    variable = disp_z
    boundary = z_plane
    value = 0.0
  []
  [no_y]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  []
  [Pressure]
    [innerPressure]
      boundary = inner
      function = 10e6
    []
  []
[]

[Materials]
  [composite_elasticity_tensor]
    type = CompositeSiCElasticityTensor
    block = clad
    elasticity_model = GA
  []
  [composite_stress]
    type = ComputeStrainIncrementBasedStress
    block = clad
  []
  [irradiation_swelling]
    type = GenericFunctionMaterial
    block = clad
    prop_names = swelling
    prop_values = swelling_func
    outputs = all
  []
[]

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

[Executioner]
  solve_type = Newton
  type = Steady
  l_tol = 1E-5
  nl_abs_tol = 1E-3
  nl_rel_tol = 1E-8
  l_max_its = 200
  nl_max_its = 400

  petsc_options_iname = '-pc_type -pc_asm_overlap -sub_pc_type -ksp_type -ksp_gmres_restart'
  petsc_options_value = ' asm      2              lu            gmres     200'
[]

[Postprocessors]
  [hoop_stress]
    type = ElementalVariableValue
    variable = hoop_stress
    elementid = 231
  []
  [hoop_strain]
    type = ElementalVariableValue
    variable = strain_zz
    elementid = 231
  []
[]

[Outputs]
  csv = true
  color = false
[]

(test/tests/solid_mechanics/compositeSiC_mechanics/hoop_2D_ad.i)

# This test case is prepared to test the elastic moduli of composite (CVI)
# SiC/SiC using the correlation provided in:
#
# L. L. Snead, T. Nozawa, Y. Katoh, T-S. Byun, S. Kondo, D. A. Petti,
# "Handbook of SiC propeties for fuel performance modeling," Journal of Nuclear
# Materials, 371, 2007, p. 329-377.
#
# The geometry represents a ring of monolothic SiC with inner radius of 0.005 m,
# outer radius of 0.0055 m, and a height of 0.01 m.
#
# A radial pressure (10 MPa) is applied to the inner surface of the ring.
# The strain calculations in the hoop direction are shown below.
#
# Subscript r: radial    y: axial    t: hoop
#
# epsilon_tt = term1 + term2 + term3
#              term1 = -nu_rt * sigma_rr / E_r
#              term2 = -nu_yt * sigma_yy / E_y
#              term3 = sigma_tt / E_t
#
# nu_rt = nu_tr * E_r / E_t = 0.08125
# nu_yr = nu_ry * E_y / E_r = 0.25
#
# epsilon_tt = term1 + term2 + term3
#            = -8.125e-6 + 0.0 + 0.00065625
#            = 0.000648125
#
# The results of BISON compared to the analytical solution:
#
#          BISON               Analytical
#       strain_hoop (m/m)    strain_hoop (m/m)
#        6.28499e-04            6.48125e-04
#
# The small difference is due to using approximation of thin cylinder for
# sigma_tt calculation

[GlobalParams]
  displacements = 'disp_x disp_y'
  volumetric_locking_correction = true
[]

[Mesh]
  coord_type = RZ
  use_displaced_mesh = false
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0.005
    xmax = 0.0055
    ymin = 0
    ymax = 0.01
  []
[]

[Functions]
  [pressure_function]
    type = PiecewiseLinear
    x = '0  1'
    y = '1  1'
  []
  [swelling_func]
    type = ParsedFunction
    expression = 0.0
  []
[]

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = SMALL
        block = 0
        add_variables = true
        generate_output = 'elastic_strain_xx elastic_strain_yy elastic_strain_zz stress_xx stress_yy stress_zz'
        use_automatic_differentiation = true
      []
    []
  []
[]

[BCs]
  [Pressure]
    [inner_pressure]
      boundary = left
      function = pressure_function
      factor = 10e6
      use_automatic_differentiation = true
    []
  []

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

[Materials]
  [elasticity_tensor]
    type = ADCompositeSiCElasticityTensor
  []
  [stress]
    type = ADComputeLinearElasticStress
  []
  [clad_density]
    type = ADStrainAdjustedDensity
    strain_free_density = 2700.0
  []
  [irradiation_swelling]
    type = ADGenericFunctionMaterial
    prop_names = swelling
    prop_values = swelling_func
    outputs = all
  []
[]

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

[Executioner]
  type = Transient

  solve_type = Newton

  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-8
  nl_abs_tol = 1e-8
  l_tol = 1e-8
  start_time = 0.0
  end_time = 1
  dt = 1
[]

[Postprocessors]
  [elastic_strain_hoop]
    type = ElementAverageValue
    variable = elastic_strain_zz
  []
[]

[Outputs]
  csv = true
[]

(examples/accident_tolerant_fuel/u3si2_sic/u3si2_outer_monolith_1.5D.i)

# Model is of a 10 pellet fuel rodlet modeled in 1.5D.  The rodlet contains
# U3Si2 fuel and a multilayer silicon carbide cladding (an inner composite
# winding layer) and an outer monolithic layer. The inner composite layer is
# 0.75 mm thick and the outer monolithic layer is 0.25 mm thick. The internal
# layered1D mesh generator can model a clad an arbitrary number of additional blocks.
# Therefore, to create the multilayer SiC clad the composite layer is assigned to the
# clad block and the monolithic_layer is assigned to the monolithic_layer block
# as specified in the additional_block_names parameter in the Mesh block.


initial_fuel_density = 11590.0

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

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

[Mesh]
  coord_type = RZ
  [layered1D_mesh]
    type = Layered1DMeshGenerator
    slices_per_block = 10
    clad_gap_width = 8.0e-5
    clad_mesh_density = customize
    clad_thickness = 0.00075
    nx_c = 5
    additional_block_names = 'monolithic_layer'
    additional_elements_per_ring = '3'
    additional_ring_thicknesses = '0.00025'
    fuel_height = 0.1186
    plenum_height = 0.027
  []
  patch_update_strategy = auto
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[UserObjects]
  [pin_geometry]
    type = Layered1DFuelPinGeometry
    mesh_generator = layered1D_mesh
  []
[]

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

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

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 1e4 1e8'
    y = '0 25000 25000'
    scale_factor = 1
  []
  [axial_peaking_factors]
    type = ParsedFunction
    expression = 1
  []
  [pressure_ramp]
    type = PiecewiseLinear
    x = '-200 0 1e8'
    y = '6.537e-3 1 1'
  []
  [q]
    type = CompositeFunction
    functions = 'power_history axial_peaking_factors'
  []
  [clad_axial_pressure]
    type = CladdingAxialPressureFunction
    plenum_pressure = plenum_pressure
    coolant_pressure = pressure_ramp
    coolant_pressure_scaling_factor = 15.5e6
    fuel_pin_geometry = pin_geometry
  []
  [fuel_axial_pressure]
    type = ParsedFunction
    expression = plenum_pressure
    symbol_names = plenum_pressure
    symbol_values = plenum_pressure
  []
[]

[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
    extra_vector_tags = 'ref'
  []
[]

[Physics]
  [SolidMechanics]
    [Layered1D]
      [fuel]
        block = fuel
        add_variables = true
        strain = SMALL
        incremental = true
        add_scalar_variables = true
        out_of_plane_strain_name = strain_yy
        fuel_pin_geometry = pin_geometry
        out_of_plane_pressure_function = fuel_axial_pressure
        eigenstrain_names = 'fuel_thermal_strain fuel_swelling_strain'
        generate_output = 'stress_xx stress_yy stress_zz vonmises_stress strain_xx'
        extra_vector_tags = 'ref'
        group_scalar_vars_in_reference_residual = true
        mesh_generator = layered1D_mesh
      []
      [composite]
        block = clad
        add_variables = true
        strain = SMALL
        incremental = true
        add_scalar_variables = true
        out_of_plane_strain_name = strain_yy
        fuel_pin_geometry = pin_geometry
        out_of_plane_pressure_function = clad_axial_pressure
        eigenstrain_names = 'composite_thermal_strain composite_swelling_strain'
        generate_output = 'stress_xx stress_yy stress_zz vonmises_stress strain_xx'
        extra_vector_tags = 'ref'
        group_scalar_vars_in_reference_residual = true
        mesh_generator = layered1D_mesh
      []
      [monolith]
        block = monolithic_layer
        add_variables = true
        strain = SMALL
        incremental = true
        add_scalar_variables = true
        out_of_plane_strain_name = strain_yy
        fuel_pin_geometry = pin_geometry
        out_of_plane_pressure_function = clad_axial_pressure
        eigenstrain_names = 'monolith_thermal_strain monolith_swelling_strain'
        generate_output = 'stress_xx stress_yy stress_zz vonmises_stress strain_xx'
        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
    order = CONSTANT
    family = MONOMIAL
    fuel_pin_geometry = pin_geometry
    fuel_volume_ratio = 1.0
    RPF = RPF
    fuel_type = U3Si2
  []
[]

[AuxKernels]
  [fast_neutron_flux]
    type = MaterialRealAux
    variable = fast_neutron_flux
    property = fast_neutron_flux
    block = 'clad monolithic_layer'
    execute_on = timestep_begin
  []
  [fast_neutron_fluence]
    type = MaterialRealAux
    variable = fast_neutron_fluence
    property = fast_neutron_fluence
    block = 'clad monolithic_layer'
    execute_on = timestep_begin
  []
  [grain_radius]
    type = GrainRadiusAux
    block = fuel
    variable = grain_radius
    temperature = temperature
    execute_on = linear
  []
  [solid_swell]
    type = MaterialRealAux
    variable = solid_swell
    property = solid_swelling
    execute_on = timestep_end
    block = fuel
  []
  [gas_swell]
    type = MaterialRealAux
    variable = gaseous_swelling
    property = gaseous_swelling
    execute_on = timestep_end
    block = fuel
  []
  [densification]
    type = MaterialRealAux
    variable = densification
    property = densification
    execute_on = timestep_end
    block = fuel
  []
  [volumetric_swelling_strain]
    type = MaterialRealAux
    variable = volumetric_swelling_strain
    property = volumetric_swelling_strain
    execute_on = timestep_end
    block = fuel
  []
[]

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

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

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

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

[Materials]
  [flux]
    type = FastNeutronFlux
    calculate_fluence = true
    block = 'clad monolithic_layer'
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    factor = 3e13
  []
  ### U3Si2 Fuel
  [fuel_thermal]
    type = SilicideFuelThermal
    block = fuel
    thermal_conductivity_model = WHITE
    temperature = temperature
  []
  [fuel_elasticity_tensor]
    type = U3Si2ElasticityTensor
    block = fuel
  []
  [fuel_stress]
    type = ComputeMultipleInelasticStress
    block = fuel
    tangent_operator = elastic
    inelastic_models = 'fuel_creep'
  []
  [fuel_creep]
    type = U3Si2CreepUpdate
    block = fuel
    temperature = temperature
  []
  [fuel_thermal_expansion]
    type = U3Si2ThermalExpansionEigenstrain
    block = fuel
    temperature = temperature
    stress_free_temperature = 295.0
    eigenstrain_name = fuel_thermal_strain
  []
  [fuel_volumetric_swelling]
    type = U3Si2VolumetricSwellingEigenstrain
    block = fuel
    gaseous_swelling_type = FINLAY
    temperature = temperature
    burnup_function = burnup
    eigenstrain_name = fuel_swelling_strain
  []
  [fission_gas_release]
    type = U3Si2Sifgrs
    block = fuel
    temperature = temperature
    burnup_function = burnup
    grain_radius_const = 2.5e-05
  []

  [fuel_density]
    type = StrainAdjustedDensity
    block = fuel
    strain_free_density = ${initial_fuel_density}
  []

  ### Composite SiC
  [composite_thermal]
    type = CompositeSiCThermal
    thermal_conductivity_model = STONE
    temperature = temperature
    block = clad
  []
  [composite_density]
    type = StrainAdjustedDensity
    block = clad
    strain_free_density = 2700.0
  []
  [composite_elasticity_tensor]
    type = CompositeSiCElasticityTensor
    block = clad
  []
  [composite_stress]
    type = ComputeStrainIncrementBasedStress
    block = clad
  []
  [composite_thermal_expansion]
    type = CompositeSiCThermalExpansionEigenstrain
    block = clad
    stress_free_temperature = 295.0
    temperature = temperature
    eigenstrain_name = composite_thermal_strain
  []
  [composite_irradiation_swelling]
    type = CompositeSiCVolumetricSwellingEigenstrain
    block = clad
    temperature = temperature
    fast_neutron_fluence = fast_neutron_fluence
    swelling_model = KATOH
    number_of_substeps = 1000
    eigenstrain_name = composite_swelling_strain
  []

  ### Monolithic SiC
  [monolith_thermal]
    type = MonolithicSiCThermal
    temperature = temperature
    thermal_conductivity_model = STONE
    block = monolithic_layer
  []
  [monolith_density]
    type = StrainAdjustedDensity
    block = monolithic_layer
    strain_free_density = 3120.0
  []
  [monolith_elasticity_tensor]
    type = MonolithicSiCElasticityTensor
    block = monolithic_layer
  []
  [monolith_stress]
    type = ComputeMultipleInelasticStress
    block = monolithic_layer
    tangent_operator = elastic
    inelastic_models = 'monolith_creep'
  []
  [monolith_creep]
    type = MonolithicSiCCreepUpdate
    block = monolithic_layer
    fast_neutron_flux = fast_neutron_flux
    temperature = temperature
    k_function = 2e-37
  []
  [monolith_thermal_expansion]
    type = MonolithicSiCThermalExpansionEigenstrain
    block = monolithic_layer
    stress_free_temperature = 295.0
    temperature = temperature
    eigenstrain_name = monolith_thermal_strain
  []
  [monolith_irradiation_swelling]
    type = CompositeSiCVolumetricSwellingEigenstrain
    block = monolithic_layer
    temperature = temperature
    fast_neutron_fluence = fast_neutron_fluence
    swelling_model = KATOH
    number_of_substeps = 1000
    eigenstrain_name = monolith_swelling_strain
  []
[]

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

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

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

  line_search = 'none'

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

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

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

[Postprocessors]
  [ave_temp_interior] # average temperature of the cladding interior and all pellet exteriors
    type = LayeredSideAverageValuePostprocessor
    boundary = 9
    variable = temperature
    execute_on = 'initial linear'
    fuel_pin_geometry = pin_geometry
  []
  [clad_inner_vol] # volume inside of cladding
    type = LayeredInternalVolumePostprocessor
    boundary = 7
    component = 0
    fuel_pin_geometry = pin_geometry
    out_of_plane_strain = strain_yy
    execute_on = 'initial linear'
  []
  [pellet_volume] # fuel pellet total volume
    type = LayeredInternalVolumePostprocessor
    boundary = 8
    component = 0
    fuel_pin_geometry = pin_geometry
    out_of_plane_strain = strain_yy
    execute_on = 'initial linear'
  []
  [avg_clad_temp] # average temperature of cladding interior
    type = LayeredSideAverageValuePostprocessor
    boundary = 7
    variable = temperature
    fuel_pin_geometry = pin_geometry
    execute_on = 'initial linear'
  []
  [fis_gas_produced] # fission gas produced (moles)
    type = LayeredElementIntegralFisGasGeneratedSifgrsPostprocessor
    block = fuel
    fuel_pin_geometry = pin_geometry
  []
  [fis_gas_released] # fission gas released to plenum (moles)
    type = LayeredElementIntegralFisGasReleasedSifgrsPostprocessor
    block = fuel
    fuel_pin_geometry = pin_geometry
  []
  [fis_gas_grain]
    type = LayeredElementIntegralFisGasGrainSifgrsPostprocessor
    block = fuel
    outputs = exodus
    fuel_pin_geometry = pin_geometry
  []
  [fis_gas_boundary]
    type = LayeredElementIntegralFisGasBoundarySifgrsPostprocessor
    block = fuel
    outputs = exodus
    fuel_pin_geometry = pin_geometry
  []
  [fission_gas_release]
    type = FGRPercent
    fission_gas_released = fis_gas_released
    fission_gas_generated = fis_gas_produced
  []

  [gas_volume]
    type = LayeredInternalVolumePostprocessor
    boundary = 9
    execute_on = 'initial linear'
    component = 0
    out_of_plane_strain = strain_yy
    fuel_pin_geometry = pin_geometry
  []
  [flux_from_clad] # area integrated heat flux from the cladding
    type = LayeredSideFluxIntegralPostprocessor
    variable = temperature
    boundary = 5
    diffusivity = thermal_conductivity
    fuel_pin_geometry = pin_geometry
  []
  [flux_from_fuel] # area integrated heat flux from the fuel
    type = LayeredSideFluxIntegralPostprocessor
    variable = temperature
    boundary = 10
    diffusivity = thermal_conductivity
    fuel_pin_geometry = pin_geometry
  []
  [_dt] # time step
    type = TimestepSize
  []
  [num_lin_it]
    type = NumLinearIterations
  []
  [num_nonlin_it]
    type = NumNonlinearIterations
  []
  [tot_lin_it]
    type = CumulativeValuePostprocessor
    postprocessor = num_lin_it
  []
  [tot_nonlin_it]
    type = CumulativeValuePostprocessor
    postprocessor = num_nonlin_it
  []
  [alive_time]
    type = PerfGraphData
    section_name = Root
    data_type = TOTAL
  []

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

  [solid_swelling]
    type = ElementAverageValue
    variable = solid_swell
    block = fuel
  []
  [gaseous_swelling]
    type = ElementAverageValue
    variable = gaseous_swelling
    block = fuel
  []
  [densification]
    type = ElementAverageValue
    variable = densification
    block = fuel
  []
  [volumetric_swelling]
    type = ElementAverageValue
    variable = volumetric_swelling_strain
    block = fuel
  []
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  color = false
[]

Description
Example Input Syntax
Input Parameters
Input Files
References