U3Si2CreepUpdate

Calculates the thermal creep behavior of U3Si2 fuel. This material must be run in conjunction with ComputeMultipleInelasticStress.

Description

Thermal creep of U $var element = document.getElementById("moose-equation-3078d253-9cb7-41fc-9e2f-3cd515e9be14");katex.render("_3", element, {displayMode:false,throwOnError:false});$ Si $var element = document.getElementById("moose-equation-f136b71c-7385-44c2-b4a6-2a0dfbea488b");katex.render("_2", element, {displayMode:false,throwOnError:false});$ fuel is calculated by the U3Si2CreepUpdate model. This model must be run in conjunction with ComputeMultipleInelasticStress. Three models for creep exist for U $var element = document.getElementById("moose-equation-90d74224-abfc-4d09-8185-11bf1037a175");katex.render("_3", element, {displayMode:false,throwOnError:false});$ Si $var element = document.getElementById("moose-equation-e5953eab-640b-4729-970c-b1d12d60af71");katex.render("_2", element, {displayMode:false,throwOnError:false});$ : Freeman, Metzger, and Yingling.

Freeman Model

The thermal creep correlation is in the form of an Arrhenius law based upon experiments completed at the University of South Carolina. Details of the derivation of the pre-exponential constant, stress exponent, and activation energy can be found in Freeman et al. (2018).

The creep rate is given by: $var element = document.getElementById("moose-equation-2f96156f-64f2-4c1e-b310-4f3680a154d7");katex.render(" \\dot{\\epsilon} = 2.0386 \\times 10^{-4} \\sigma^{1.2063} \\exp\\left( \\frac{295550}{RT}\\right)", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-e4307866-dbd1-4129-b0de-68e1e2bd3a6d");katex.render("R", element, {displayMode:false,throwOnError:false});$ is the ideal gas constant with a value of 8.314 J/mol-K and $var element = document.getElementById("moose-equation-03488309-0426-4777-a4ec-b42cf04b96a4");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the temperature in Kelvin.

Metzger Model

The mechanism for creep can fall into one of three regimes based upon the temperature and stress state within the fuel. This model is developed by Metzger (2016). Below a homologous temperature of 0.45 T $var element = document.getElementById("moose-equation-293cd8d3-39ea-489b-9906-200451c03992");katex.render("_m", element, {displayMode:false,throwOnError:false});$ (872.0 K) the creep is said to be athermal and irradiation induced. This is known as the Nabarro-Herring regime with a creep rate given by:

$var element = document.getElementById("moose-equation-dcfa525f-9807-4b5e-a54f-3c23df930d15");katex.render("\\dot{\\epsilon}_{NH} = \\frac{A_{NH}D_{irr}b^3\\sigma}{kTd^2}", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-1286031c-951a-41a0-8cc1-c5fa66da97cc");katex.render("A_{NH}", element, {displayMode:false,throwOnError:false});$ is Nabarro-Herring creep coefficient, $var element = document.getElementById("moose-equation-eef6b7b1-56d7-4275-94ab-dd950ab15c5e");katex.render("D_{irr}", element, {displayMode:false,throwOnError:false});$ is the diffusion coefficient under athermal conditions, $var element = document.getElementById("moose-equation-852d9b8e-f365-4d34-b518-2dbf8a8fb7a8");katex.render("b", element, {displayMode:false,throwOnError:false});$ is the burgers vector, $var element = document.getElementById("moose-equation-0a0f7088-c30d-4994-a751-896145ed343c");katex.render("\\sigma", element, {displayMode:false,throwOnError:false});$ is the Von Mises stress in Pa, $var element = document.getElementById("moose-equation-21e51704-ed55-4d55-bf48-4b418b47ea6f");katex.render("k", element, {displayMode:false,throwOnError:false});$ is the Boltzmann constant, $var element = document.getElementById("moose-equation-e012932f-c994-4bef-a572-eef408b7cbd8");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the temperature in K, and $var element = document.getElementById("moose-equation-fd0231f4-cef1-415d-b2f8-42761fe80875");katex.render("d", element, {displayMode:false,throwOnError:false});$ is the grain size. Above a homologous temperature of 0.45 T $var element = document.getElementById("moose-equation-760e1a84-3726-444f-bb46-dec08379a471");katex.render("_m", element, {displayMode:false,throwOnError:false});$ the creep can either driven by a Coble grain boundary creep mechanism or a dislocation creep mechanism depending upon the stress state. If $var element = document.getElementById("moose-equation-92586fd9-c58e-4270-ba56-19e458474aff");katex.render("\\sigma/G > 10^{-4}", element, {displayMode:false,throwOnError:false});$ where $var element = document.getElementById("moose-equation-fba7a971-0582-451b-8d24-24cea0816dfe");katex.render("G", element, {displayMode:false,throwOnError:false});$ is the shear modulus of U $var element = document.getElementById("moose-equation-a875caa7-7146-4e0b-ac11-242c02de6323");katex.render("_3", element, {displayMode:false,throwOnError:false});$ Si $var element = document.getElementById("moose-equation-dabe8917-4fe4-4f2d-9b1b-7c2e6a675612");katex.render("_2", element, {displayMode:false,throwOnError:false});$ , which is assumed to be constant at 50 GPa in this model, then the creep mechanism is governed by dislocations through:

$var element = document.getElementById("moose-equation-706a855d-3bb5-4862-9573-16dd321802eb");katex.render("\\dot{\\epsilon}_{dis} = \\frac{A_{dis}D_{L}b\\sigma^5}{kTG^4}", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-0bd0531f-4710-4139-b459-4668e20e441a");katex.render("A_{dis}", element, {displayMode:false,throwOnError:false});$ is the dislocation creep coefficient, and $var element = document.getElementById("moose-equation-396ebfc9-94f9-46c2-8b7c-8890c87e0fa0");katex.render("D_{L}", element, {displayMode:false,throwOnError:false});$ is the lattice diffusion coefficient. if $var element = document.getElementById("moose-equation-c655e5b0-a6ce-4209-8702-953f20afa0a4");katex.render("\\sigma/G \\leq 10^{-4}", element, {displayMode:false,throwOnError:false});$ then creep occurs due to a Coble grain boundary mechanism given by:

$var element = document.getElementById("moose-equation-eb2cc71d-7f44-4eca-95d9-ac5687c8c4da");katex.render("\\dot{\\epsilon}_{co} = \\frac{A_{co}D_{gb}b^4\\sigma}{kTd^3}", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-851a5b53-b3e7-4d11-9b51-de5fb82d9ef0");katex.render("A_{co}", element, {displayMode:false,throwOnError:false});$ is the Coble creep coefficient, and $var element = document.getElementById("moose-equation-71736e7b-a67f-4dc7-8e51-58241bef5f6c");katex.render("D_{gb}", element, {displayMode:false,throwOnError:false});$ is the grain boundary coefficient. The diffusion coefficients in the three different regimes are given by:

$var element = document.getElementById("moose-equation-17171cb6-24cd-48b8-8c3b-bfa9270e5013");katex.render("D_{irr} = D_{o_{gb}} \\exp\\left(-\\frac{Q_{gb}}{872.0k}\\right)", element, {displayMode:true,throwOnError:false});$ $var element = document.getElementById("moose-equation-67e51db7-8a18-4707-95ac-b60d0016815a");katex.render("D_{gb} = D_{o_{gb}} \\exp\\left(-\\frac{Q_{gb}}{kT}\\right)", element, {displayMode:true,throwOnError:false});$ $var element = document.getElementById("moose-equation-eb47ecc0-6c04-4a39-be85-4bc56f44de8f");katex.render("D_{L} = D_{o_{L}} \\exp\\left(-\\frac{Q_{L}}{kT}\\right)", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-501399d3-e483-4fcb-b47a-c4e7e303abb4");katex.render("D_{o_{gb}}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-b5199303-7bb5-4ee9-976f-d964421337d2");katex.render("D_{o_{L}}", element, {displayMode:false,throwOnError:false});$ are coefficients, and $var element = document.getElementById("moose-equation-1a784953-c28b-4a35-a626-6753a7d61909");katex.render("Q_{gb}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-cdfc3863-689b-409e-991f-eedde6018d48");katex.render("Q_{L}", element, {displayMode:false,throwOnError:false});$ are activation energies. The values assumed by the model for all the parameters are summarized in Table 1.

Table 1: Values for different parameters used in the Metzger creep model.

Parameter	Value	Units
$var element = document.getElementById("moose-equation-aac316ce-4ec1-4455-8d35-609b7b07fb01");katex.render("A_{NH}", element, {displayMode:false,throwOnError:false});$	12.5	dimensionless
$var element = document.getElementById("moose-equation-66521925-1c0f-4dfc-b731-838c0c8ee7b9");katex.render("A_{dis}", element, {displayMode:false,throwOnError:false});$	6 $var element = document.getElementById("moose-equation-6f9d14ee-d721-4bdc-9f09-07eb963cfe44");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-a5df9b0a-9c55-4f25-80c9-8bbc9003417b");katex.render("^7", element, {displayMode:false,throwOnError:false});$	dimensionless
$var element = document.getElementById("moose-equation-534c7d4e-68e9-49ea-8d0c-ca967d3f3574");katex.render("A_{co}", element, {displayMode:false,throwOnError:false});$	40.0	dimensionless
$var element = document.getElementById("moose-equation-d809ea9a-79de-4338-8c79-533b04fa5c71");katex.render("D_{o_{gb}}", element, {displayMode:false,throwOnError:false});$	2365.0	m $var element = document.getElementById("moose-equation-2a1edd0e-bb39-461c-8b69-f769485762f7");katex.render("^2", element, {displayMode:false,throwOnError:false});$ /s
$var element = document.getElementById("moose-equation-4603bbd8-563c-431e-ae29-f8870dc856c8");katex.render("D_{o_{L}}", element, {displayMode:false,throwOnError:false});$	6.86 $var element = document.getElementById("moose-equation-726b2bcb-bc7a-45d9-b390-4dc97c8c4fca");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-e7f5d07f-c74f-40a0-bd1d-31059cfbdd4a");katex.render("^{24}", element, {displayMode:false,throwOnError:false});$	m $var element = document.getElementById("moose-equation-bd6a2fef-d91d-4b78-85f2-d0144d316725");katex.render("^2", element, {displayMode:false,throwOnError:false});$ /s
$var element = document.getElementById("moose-equation-69220995-3022-4ccc-b764-06ff806737d7");katex.render("Q_{gb}", element, {displayMode:false,throwOnError:false});$	9.97 $var element = document.getElementById("moose-equation-e79ab38c-1cbf-4b1c-82ac-949176a6ab0b");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-8c95c58b-d76d-4fb5-9c2b-c6142b9372f7");katex.render("^{-19}", element, {displayMode:false,throwOnError:false});$	J
$var element = document.getElementById("moose-equation-e62ff996-df4f-4cac-af6b-1d701d115a04");katex.render("Q_{L}", element, {displayMode:false,throwOnError:false});$	2.0 $var element = document.getElementById("moose-equation-5f208793-2b47-4788-aecb-5ad12082a086");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-52fb76ef-4c51-4b97-9448-58bc09c6feec");katex.render("^{-18}", element, {displayMode:false,throwOnError:false});$	J
$var element = document.getElementById("moose-equation-4ad5adea-f15a-468a-b563-0084e362d3e2");katex.render("b", element, {displayMode:false,throwOnError:false});$	0.56 $var element = document.getElementById("moose-equation-a95e3162-f049-4cd4-b1cf-d53ad39e7bcc");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-a925c7ec-23ab-474a-8272-92572b89f886");katex.render("^{-9}", element, {displayMode:false,throwOnError:false});$	m
$var element = document.getElementById("moose-equation-8506885a-8511-4814-a606-6a204c2aff33");katex.render("k", element, {displayMode:false,throwOnError:false});$	1.38064852 $var element = document.getElementById("moose-equation-8af42fc8-9747-493b-a1f4-3a8bec39c52d");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-4da4343e-d795-4fdc-96fd-8405e5ee125b");katex.render("^{-23}", element, {displayMode:false,throwOnError:false});$	J/K
$var element = document.getElementById("moose-equation-a5112dbe-4033-4992-8d79-bb4b9437701a");katex.render("d", element, {displayMode:false,throwOnError:false});$	20 $var element = document.getElementById("moose-equation-4b51e776-d1e2-47da-b19e-3951ab6e9bfd");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-edb9226e-fa34-4cc3-a268-4c4d41cc0534");katex.render("^{-6}", element, {displayMode:false,throwOnError:false});$	m

Yingling Model

The thermal creep correlation is the Mukherjee-Bird-Dorn empirical equation based upon experiments completed at the University of South Carolina. Details of the derivation of the pre-exponential constant, stress and grain size exponents, and activation energy can be found in Yingling (2019).

The creep rate is given by: $var element = document.getElementById("moose-equation-68d5b1b8-13b7-4c41-acbb-7e7e3d17272b");katex.render(" \\dot{\\epsilon} = 4.841 \\times 10^{-19} \\sigma^{1.936} d^{-1.86} \\exp\\left( \\frac{223100}{RT}\\right)", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-c203c719-1a61-47fd-a04e-2073fd50af7b");katex.render("R", element, {displayMode:false,throwOnError:false});$ is the ideal gas constant with a value of 8.314 J/mol-K, $var element = document.getElementById("moose-equation-0f18e644-c48f-43a5-a893-3efb4fdbde25");katex.render("d", element, {displayMode:false,throwOnError:false});$ is the average grain size in meters, and $var element = document.getElementById("moose-equation-c443b723-77e3-4c85-a120-fcc888f76a0c");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the temperature in Kelvin.

Example Input Syntax

An example of using the Freeman model:

[Materials<<<{"href": "../../../syntax/Materials/index.html"}>>>]
  [u3si2creep]
    type = U3Si2CreepUpdate<<<{"description": "Calculates the thermal creep behavior of U3Si2 fuel.  This material must be run in conjunction with ComputeMultipleInelasticStress.", "href": "U3Si2CreepUpdate.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 1
    temperature<<<{"description": "The coupled temperature (K)"}>>> = temp
  []
[]

An example of using the Metzger model

[Materials<<<{"href": "../../../syntax/Materials/index.html"}>>>]
  [u3si2creep]
    type = U3Si2CreepUpdate<<<{"description": "Calculates the thermal creep behavior of U3Si2 fuel.  This material must be run in conjunction with ComputeMultipleInelasticStress.", "href": "U3Si2CreepUpdate.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 1
    creep_model<<<{"description": "The model to be used for thermal creep."}>>> = METZGER
    temperature<<<{"description": "The coupled temperature (K)"}>>> = temp
  []
[]

U3Si2CreepUpdate must be run in conjunction with the inelastic strain return mapping stress calculator as shown below:

[Materials<<<{"href": "../../../syntax/Materials/index.html"}>>>]
  [stress]
    type = ComputeMultipleInelasticStress<<<{"description": "Compute state (stress and internal parameters such as plastic strains and internal parameters) using an iterative process.  Combinations of creep models and plastic models may be used.", "href": "../ComputeMultipleInelasticStress.html"}>>>
    tangent_operator<<<{"description": "Type of tangent operator to return.  'elastic': return the elasticity tensor.  'nonlinear': return the full, general consistent tangent operator."}>>> = elastic
    inelastic_models<<<{"description": "The material objects to use to calculate stress and inelastic strains. Note: specify creep models first and plasticity models second."}>>> = 'u3si2creep'
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 1
  []
[]

Input Parameters

temperatureThe coupled temperature (K)
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The coupled temperature (K)

Required Parameters

absolute_tolerance1e-11Absolute convergence tolerance for Newton iteration
Default:1e-11
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Absolute convergence tolerance for Newton iteration
acceptable_multiplier10Factor applied to relative and absolute tolerance for acceptable convergence if iterations are no longer making progress
Default:10
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Factor applied to relative and absolute tolerance for acceptable convergence if iterations are no longer making progress
adaptive_substeppingFalseUse adaptive substepping, where the number of substeps is successively doubled until the return mapping model successfully converges or the maximum number of substeps is reached.
Default:False
C++ Type:bool
Controllable:No
Description:Use adaptive substepping, where the number of substeps is successively doubled until the return mapping model successfully converges or the maximum number of substeps is reached.
automatic_differentiation_return_mappingFalseWhether to use automatic differentiation to compute the derivative.
Default:False
C++ Type:bool
Controllable:No
Description:Whether to use automatic differentiation to compute the derivative.
base_nameOptional parameter that defines a prefix for all material properties related to this stress update model. This allows for multiple models of the same type to be used without naming conflicts.
C++ Type:std::string
Controllable:No
Description:Optional parameter that defines a prefix for all material properties related to this stress update model. This allows for multiple models of the same type to be used without naming conflicts.
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
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
creep_modelYINGLINGThe model to be used for thermal creep.
Default:YINGLING
C++ Type:MooseEnum
Options:FREEMAN, METZGER, YINGLING
Controllable:No
Description:The model to be used for thermal creep.
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.
grain_size2e-05Constant U3Si2 grain size.
Default:2e-05
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Constant U3Si2 grain size.
grain_size_exponent_scale_factor1Scale factor to be applied to the thermal creep grain size exponent. Used for calibration and sensitivity studies.
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Scale factor to be applied to the thermal creep grain size exponent. Used for calibration and sensitivity studies.
max_inelastic_increment0.0001The maximum inelastic strain increment allowed in a time step
Default:0.0001
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The maximum inelastic strain increment allowed in a time step
maximum_number_substeps25The maximum number of substeps allowed before cutting the time step.
Default:25
C++ Type:unsigned int
Controllable:No
Description:The maximum number of substeps allowed before cutting the time step.
relative_tolerance1e-08Relative convergence tolerance for Newton iteration
Default:1e-08
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Relative convergence tolerance for Newton iteration
use_substep_integration_errorFalseIf true, it establishes a substep size that will yield, at most,the creep numerical integration error given by substep_strain_tolerance.
Default:False
C++ Type:bool
Controllable:No
Description:If true, it establishes a substep size that will yield, at most,the creep numerical integration error given by substep_strain_tolerance.
use_substeppingNONEWhether and how to use substepping
Default:NONE
C++ Type:MooseEnum
Options:NONE, ERROR_BASED, INCREMENT_BASED
Controllable:No
Description:Whether and how to use substepping

Optional Parameters

activation_energy_scale_factor1Scale factor to be applied to the thermal creep activation energy. Used for calibration and sensitivity studies.
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Scale factor to be applied to the thermal creep activation energy. Used for calibration and sensitivity studies.
pre_exponential_scale_factor1Scale factor to be applied to the thermal creep exponential prefactor. Used for calibration and sensitivity studies.
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Scale factor to be applied to the thermal creep exponential prefactor. Used for calibration and sensitivity studies.
stress_exponent_scale_factor1Scale factor to be applied to the thermal creep stress exponent. Used for calibration and sensitivity studies.
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Scale factor to be applied to the thermal creep stress exponent. Used for calibration and sensitivity studies.

Advanced: Scaling Factors Parameters

apply_strainTrueFlag to apply strain. Used for testing.
Default:True
C++ Type:bool
Controllable:No
Description:Flag to apply strain. Used for testing.
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.
effective_inelastic_strain_nameeffective_creep_strainName of the material property that stores the effective inelastic strain
Default:effective_creep_strain
C++ Type:std::string
Controllable:No
Description:Name of the material property that stores the effective inelastic strain
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
substep_strain_tolerance0.1Maximum ratio of the initial elastic strain increment at start of the return mapping solve to the maximum inelastic strain allowable in a single substep. Reduce this value to increase the number of substeps
Default:0.1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Maximum ratio of the initial elastic strain increment at start of the return mapping solve to the maximum inelastic strain allowable in a single substep. Reduce this value to increase the number of substeps
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

internal_solve_full_iteration_historyFalseSet true to output full internal Newton iteration history at times determined by `internal_solve_output_on`. If false, only a summary is output.
Default:False
C++ Type:bool
Controllable:No
Description:Set true to output full internal Newton iteration history at times determined by `internal_solve_output_on`. If false, only a summary is output.
internal_solve_output_onon_errorWhen to output internal Newton solve information
Default:on_error
C++ Type:MooseEnum
Options:never, on_error, always
Controllable:No
Description:When to output internal Newton solve information

Debug 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

(examples/accident_tolerant_fuel/u3si2_sic/u3si2_outer_monolith_1.5D.i)
(test/tests/solid_mechanics/u3si2_creep/thermal_creep_u3si2_metzger.i)
(test/tests/solid_mechanics/u3si2_creep/thermal_creep_u3si2.i)
(examples/accident_tolerant_fuel/u3si2_zircaloy/u3si2_zircaloy.i)

References

R. A. Freeman, T. Martin, E. Roberts, and T. W. Knight. Analysis of thermal creep for uranium silicide fuel using Bison. In Proceedings of the 2018 International Congress on Advances in Nuclear Power Plants (ICAPP 18). Charlotte, NC, 2018.

@inproceedings{freeman2018,
    author = "Freeman, R. A. and Martin, T. and Roberts, E. and Knight, T. W.",
    title = "Analysis of thermal creep for uranium silicide fuel using {B}ison",
    year = "2018",
    booktitle = "Proceedings of the 2018 International Congress on Advances in Nuclear Power Plants (ICAPP 18)",
    address = "Charlotte, NC"
}

K. E. Metzger. Analysis of Pellet Cladding Interaction and Creep of U$_3$Si$_2$ Fuel for use in Light Water Reactors. PhD Dissertation, University of South Carolina, Department of Nuclear Engineering, 2016.

@PHDTHESIS{metzger2016,
    author = "Metzger, K. E.",
    title = "Analysis of Pellet Cladding Interaction and Creep of {U$\_3$Si$\_2$} Fuel for use in Light Water Reactors",
    school = "University of South Carolina",
    year = "2016",
    type = "{PhD} Dissertation",
    address = "Department of Nuclear Engineering"
}

J. A. Yingling. BISON simulation-based identification of important design criteria for U$_3$Si$_2$ fuels with composite-monolithic duplex SiC cladding. Master's thesis, University of South Carolina, December 2019.

@mastersthesis{yingling2019,
    author = "Yingling, J. A.",
    title = "{BISON} Simulation-Based Identification of Important Design Criteria for {U$\_3$Si$\_2$} Fuels With Composite-Monolithic Duplex {SiC} Cladding",
    school = "University of South Carolina",
    year = "2019",
    month = "December"
}

(test/tests/solid_mechanics/u3si2_creep/thermal_creep_u3si2.i)

#--------------------------------------------------------------------------------
#  This test verifies the implementation of the thermal creep model that has been
#  added to BISON for U3Si2. The model was taken from Freeman et al. 2018.
#  The correlation is given by:
#
#  e_dot = 2.0386e-4 * sigma^(1.2063) * exp(-295550 / RT)                (s^-1)
#
#  where R is the ideal gas constant taken as 8.314 J/mol-K. In this test a
#  uniaxial specimen is pulled in the axial direction by
#  a varying stress at 1000 K.  Results are taken from different times in the
#  simulation: 2e5, 4e5, 6e5, 8e5 and 1e6 seconds. The creep rate predicted by
#  BISON is compared to the analytical solution for these stresses for the central
#  element in the 2D patch test.
#
#    Time (s)    Von Mises Stress     Analytical Creep       BISON Creep
#                     (MPa)             Rate (s^-1)          Rate (s^-1)
#
#      2e5           17.1306            3.9530e-11           3.9534e-11
#      4e5           33.8741            8.9972e-11           8.9984e-11
#      6e5           50.2010            1.4761e-10           1.4763e-11
#      8e5           66.1119            5.0157e-10           5.0161e-10
#      1e6           81.6204            2.5991e-10           2.5997e-10
#--------------------------------------------------------------------------------

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  coord_type = RZ
  [mesh]
    type = FileMeshGenerator
    file = 1x1_rz.e
  []
[]

[Variables]
  [disp_x]
    order = FIRST
    family = LAGRANGE
  []
  [disp_y]
    order = FIRST
    family = LAGRANGE
  []
  [temp]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300.0
  []
[]

[AuxVariables]
  [vonmises]
    order = CONSTANT
    family = MONOMIAL
  []
  [elastic_strain_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [pressure_pull]
    type = PiecewiseLinear
    x = '0 1e6'
    y = '0 1'
  []
  [temperature_function]
    type = PiecewiseLinear
    x = '0 10000'
    y = '300 1000.0'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = finite
  []
[]

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


[AuxKernels]
  [vonmises]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = vonmises
    scalar_type = VonMisesStress
    execute_on = timestep_end
  []
  [elastic_strain_yy]
    type = RankTwoAux
    rank_two_tensor = elastic_strain
    variable = elastic_strain_yy
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  []

  [creep_strain_yy]
    type = RankTwoAux
    rank_two_tensor = creep_strain
    variable = creep_strain_yy
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  []
  [creep_rate_aux]
    type = MaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
[]

[BCs]
  [Pressure]
    [u_top_pull]
      boundary = 4
      factor = 1e8
      function = pressure_pull
    []
  []
  [bottom_fix_y]
    type = DirichletBC
    variable = disp_y
    boundary = 3
    value = 0.0
  []
  [bottom_fix_x]
    type = DirichletBC
    variable = disp_x
    boundary = 3
    value = 0.0
  []
  [temp_bc]
    type = FunctionDirichletBC
    variable = temp
    boundary = '1 2 3 4'
    function = temperature_function
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1.4e11
    poissons_ratio = 0.17
    block = 1
  []
  [stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'u3si2creep'
    block = 1
  []
  [u3si2creep]
    type = U3Si2CreepUpdate
    block = 1
    temperature = temp
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = 1
    strain_free_density = 11590.0
  []
  [thermal]
    type = HeatConductionMaterial
    block = 1
    thermal_conductivity = 1
    specific_heat = 1
  []
[]

[Executioner]
  type = Transient

  solve_type = 'PJFNK'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -ksp_gmres_restart'
  petsc_options_value = 'lu       superlu_dist                  51'

  line_search = 'none'

  l_max_its  = 100
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-10
  l_tol      = 1e-5
  start_time = 0.0
  end_time   = 1e6

  dt = 1e4
[]

[Postprocessors]
  [creep_rate]
    type = ElementalVariableValue
    variable = creep_rate_aux
    elementid = 2
  []
  [vonmises]
    type = ElementalVariableValue
    variable = vonmises
    elementid = 2
  []
  [fuel_temp]
    type = ElementalVariableValue
    elementid = 2
    variable = temp
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/solid_mechanics/u3si2_creep/thermal_creep_u3si2_metzger.i)

#--------------------------------------------------------------------------------
#  This test verifies the implementation of the thermal creep model from Metzger 2016.
#  The model consists of three different creep regimes. This input file is setup
#  to investigate the Nabarro-Herring low temperature regime (T < 872.0 K).
#  Additional tests are created to test the high temperature Coble and Dislocation
#  regimes using command line arguments in the tests file.
#
#  In this test a uniaxial specimen is pulled in the axial direction by
#  a varying stress at 800 K.  Results are taken from different times in the
#  simulation: 2e5, 4e5, 6e5, 8e5 and 1e6 seconds. The creep rate predicted by
#  BISON is compared to the analytical solution for these stresses for the central
#  element in the 2D patch test.
#
#    Time (s)    Von Mises Stress     Analytical Creep       BISON Creep
#                     (MPa)             Rate (s^-1)          Rate (s^-1)
#
#      2e5           17.2926            2.2027e-23           2.2027e-23
#      4e5           34.5857            4.4054e-23           4.4054e-23
#      6e5           51.8795            6.6082e-23           6.6082e-23
#      8e5           69.1738            8.8111e-23           8.8111e-23
#      1e6           86.4687            1.1014e-22           1.1014e-22
#--------------------------------------------------------------------------------

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  coord_type = RZ
  [mesh]
    type = FileMeshGenerator
    file = 1x1_rz.e
  []
[]

[Variables]
  [disp_x]
    order = FIRST
    family = LAGRANGE
  []
  [disp_y]
    order = FIRST
    family = LAGRANGE
  []
  [temp]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300.0
  []
[]

[AuxVariables]
  [vonmises]
    order = CONSTANT
    family = MONOMIAL
  []
  [elastic_strain_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [pressure_pull]
    type = PiecewiseLinear
    x = '0 1e6'
    y = '0 1'
  []
  [temperature_function]
    type = PiecewiseLinear
    x = '0 10000'
    y = '300 800.0'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = finite
  []
[]

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

[AuxKernels]
  [vonmises]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = vonmises
    scalar_type = VonMisesStress
    execute_on = timestep_end
  []
  [elastic_strain_yy]
    type = RankTwoAux
    rank_two_tensor = elastic_strain
    variable = elastic_strain_yy
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  []
  [creep_strain_yy]
    type = RankTwoAux
    rank_two_tensor = creep_strain
    variable = creep_strain_yy
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  []
  [creep_rate_aux]
    type = MaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
[]

[BCs]
  [Pressure]
    [u_top_pull]
      boundary = 4
      factor = 1e8
      function = pressure_pull
    []
  []
  [bottom_fix_y]
    type = DirichletBC
    variable = disp_y
    boundary = 3
    value = 0.0
  []
  [bottom_fix_x]
    type = DirichletBC
    variable = disp_x
    boundary = 3
    value = 0.0
  []
  [temp_bc]
    type = FunctionDirichletBC
    variable = temp
    boundary = '1 2 3 4'
    function = temperature_function
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1.4e11
    poissons_ratio = 0.17
    block = 1
  []
  [stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'u3si2creep'
    block = 1
  []
  [u3si2creep]
    type = U3Si2CreepUpdate
    block = 1
    creep_model = METZGER
    temperature = temp
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = 1
    strain_free_density = 11590.0
  []
  [thermal]
    type = HeatConductionMaterial
    block = 1
    thermal_conductivity = 1
    specific_heat = 1
  []
[]

[Executioner]
  type = Transient

  solve_type = 'PJFNK'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -ksp_gmres_restart'
  petsc_options_value = 'lu       superlu_dist                  51'

  line_search = 'none'

  l_max_its  = 100
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-10
  l_tol      = 1e-5
  start_time = 0.0
  end_time   = 1e6

  dt = 1e4
[]

[Postprocessors]
  [creep_rate]
    type = ElementalVariableValue
    variable = creep_rate_aux
    elementid = 2
  []
  [vonmises]
    type = ElementalVariableValue
    variable = vonmises
    elementid = 2
  []
  [fuel_temp]
    type = ElementalVariableValue
    elementid = 2
    variable = temp
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/solid_mechanics/u3si2_creep/thermal_creep_u3si2.i)

#--------------------------------------------------------------------------------
#  This test verifies the implementation of the thermal creep model that has been
#  added to BISON for U3Si2. The model was taken from Freeman et al. 2018.
#  The correlation is given by:
#
#  e_dot = 2.0386e-4 * sigma^(1.2063) * exp(-295550 / RT)                (s^-1)
#
#  where R is the ideal gas constant taken as 8.314 J/mol-K. In this test a
#  uniaxial specimen is pulled in the axial direction by
#  a varying stress at 1000 K.  Results are taken from different times in the
#  simulation: 2e5, 4e5, 6e5, 8e5 and 1e6 seconds. The creep rate predicted by
#  BISON is compared to the analytical solution for these stresses for the central
#  element in the 2D patch test.
#
#    Time (s)    Von Mises Stress     Analytical Creep       BISON Creep
#                     (MPa)             Rate (s^-1)          Rate (s^-1)
#
#      2e5           17.1306            3.9530e-11           3.9534e-11
#      4e5           33.8741            8.9972e-11           8.9984e-11
#      6e5           50.2010            1.4761e-10           1.4763e-11
#      8e5           66.1119            5.0157e-10           5.0161e-10
#      1e6           81.6204            2.5991e-10           2.5997e-10
#--------------------------------------------------------------------------------

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  coord_type = RZ
  [mesh]
    type = FileMeshGenerator
    file = 1x1_rz.e
  []
[]

[Variables]
  [disp_x]
    order = FIRST
    family = LAGRANGE
  []
  [disp_y]
    order = FIRST
    family = LAGRANGE
  []
  [temp]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300.0
  []
[]

[AuxVariables]
  [vonmises]
    order = CONSTANT
    family = MONOMIAL
  []
  [elastic_strain_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [pressure_pull]
    type = PiecewiseLinear
    x = '0 1e6'
    y = '0 1'
  []
  [temperature_function]
    type = PiecewiseLinear
    x = '0 10000'
    y = '300 1000.0'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = finite
  []
[]

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


[AuxKernels]
  [vonmises]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = vonmises
    scalar_type = VonMisesStress
    execute_on = timestep_end
  []
  [elastic_strain_yy]
    type = RankTwoAux
    rank_two_tensor = elastic_strain
    variable = elastic_strain_yy
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  []

  [creep_strain_yy]
    type = RankTwoAux
    rank_two_tensor = creep_strain
    variable = creep_strain_yy
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  []
  [creep_rate_aux]
    type = MaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
[]

[BCs]
  [Pressure]
    [u_top_pull]
      boundary = 4
      factor = 1e8
      function = pressure_pull
    []
  []
  [bottom_fix_y]
    type = DirichletBC
    variable = disp_y
    boundary = 3
    value = 0.0
  []
  [bottom_fix_x]
    type = DirichletBC
    variable = disp_x
    boundary = 3
    value = 0.0
  []
  [temp_bc]
    type = FunctionDirichletBC
    variable = temp
    boundary = '1 2 3 4'
    function = temperature_function
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1.4e11
    poissons_ratio = 0.17
    block = 1
  []
  [stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'u3si2creep'
    block = 1
  []
  [u3si2creep]
    type = U3Si2CreepUpdate
    block = 1
    temperature = temp
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = 1
    strain_free_density = 11590.0
  []
  [thermal]
    type = HeatConductionMaterial
    block = 1
    thermal_conductivity = 1
    specific_heat = 1
  []
[]

[Executioner]
  type = Transient

  solve_type = 'PJFNK'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -ksp_gmres_restart'
  petsc_options_value = 'lu       superlu_dist                  51'

  line_search = 'none'

  l_max_its  = 100
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-10
  l_tol      = 1e-5
  start_time = 0.0
  end_time   = 1e6

  dt = 1e4
[]

[Postprocessors]
  [creep_rate]
    type = ElementalVariableValue
    variable = creep_rate_aux
    elementid = 2
  []
  [vonmises]
    type = ElementalVariableValue
    variable = vonmises
    elementid = 2
  []
  [fuel_temp]
    type = ElementalVariableValue
    elementid = 2
    variable = temp
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  exodus = 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
[]

(test/tests/solid_mechanics/u3si2_creep/thermal_creep_u3si2_metzger.i)

#--------------------------------------------------------------------------------
#  This test verifies the implementation of the thermal creep model from Metzger 2016.
#  The model consists of three different creep regimes. This input file is setup
#  to investigate the Nabarro-Herring low temperature regime (T < 872.0 K).
#  Additional tests are created to test the high temperature Coble and Dislocation
#  regimes using command line arguments in the tests file.
#
#  In this test a uniaxial specimen is pulled in the axial direction by
#  a varying stress at 800 K.  Results are taken from different times in the
#  simulation: 2e5, 4e5, 6e5, 8e5 and 1e6 seconds. The creep rate predicted by
#  BISON is compared to the analytical solution for these stresses for the central
#  element in the 2D patch test.
#
#    Time (s)    Von Mises Stress     Analytical Creep       BISON Creep
#                     (MPa)             Rate (s^-1)          Rate (s^-1)
#
#      2e5           17.2926            2.2027e-23           2.2027e-23
#      4e5           34.5857            4.4054e-23           4.4054e-23
#      6e5           51.8795            6.6082e-23           6.6082e-23
#      8e5           69.1738            8.8111e-23           8.8111e-23
#      1e6           86.4687            1.1014e-22           1.1014e-22
#--------------------------------------------------------------------------------

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  coord_type = RZ
  [mesh]
    type = FileMeshGenerator
    file = 1x1_rz.e
  []
[]

[Variables]
  [disp_x]
    order = FIRST
    family = LAGRANGE
  []
  [disp_y]
    order = FIRST
    family = LAGRANGE
  []
  [temp]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300.0
  []
[]

[AuxVariables]
  [vonmises]
    order = CONSTANT
    family = MONOMIAL
  []
  [elastic_strain_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [pressure_pull]
    type = PiecewiseLinear
    x = '0 1e6'
    y = '0 1'
  []
  [temperature_function]
    type = PiecewiseLinear
    x = '0 10000'
    y = '300 800.0'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = finite
  []
[]

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

[AuxKernels]
  [vonmises]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = vonmises
    scalar_type = VonMisesStress
    execute_on = timestep_end
  []
  [elastic_strain_yy]
    type = RankTwoAux
    rank_two_tensor = elastic_strain
    variable = elastic_strain_yy
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  []
  [creep_strain_yy]
    type = RankTwoAux
    rank_two_tensor = creep_strain
    variable = creep_strain_yy
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  []
  [creep_rate_aux]
    type = MaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
[]

[BCs]
  [Pressure]
    [u_top_pull]
      boundary = 4
      factor = 1e8
      function = pressure_pull
    []
  []
  [bottom_fix_y]
    type = DirichletBC
    variable = disp_y
    boundary = 3
    value = 0.0
  []
  [bottom_fix_x]
    type = DirichletBC
    variable = disp_x
    boundary = 3
    value = 0.0
  []
  [temp_bc]
    type = FunctionDirichletBC
    variable = temp
    boundary = '1 2 3 4'
    function = temperature_function
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1.4e11
    poissons_ratio = 0.17
    block = 1
  []
  [stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'u3si2creep'
    block = 1
  []
  [u3si2creep]
    type = U3Si2CreepUpdate
    block = 1
    creep_model = METZGER
    temperature = temp
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = 1
    strain_free_density = 11590.0
  []
  [thermal]
    type = HeatConductionMaterial
    block = 1
    thermal_conductivity = 1
    specific_heat = 1
  []
[]

[Executioner]
  type = Transient

  solve_type = 'PJFNK'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -ksp_gmres_restart'
  petsc_options_value = 'lu       superlu_dist                  51'

  line_search = 'none'

  l_max_its  = 100
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-10
  l_tol      = 1e-5
  start_time = 0.0
  end_time   = 1e6

  dt = 1e4
[]

[Postprocessors]
  [creep_rate]
    type = ElementalVariableValue
    variable = creep_rate_aux
    elementid = 2
  []
  [vonmises]
    type = ElementalVariableValue
    variable = vonmises
    elementid = 2
  []
  [fuel_temp]
    type = ElementalVariableValue
    elementid = 2
    variable = temp
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/solid_mechanics/u3si2_creep/thermal_creep_u3si2.i)

#--------------------------------------------------------------------------------
#  This test verifies the implementation of the thermal creep model that has been
#  added to BISON for U3Si2. The model was taken from Freeman et al. 2018.
#  The correlation is given by:
#
#  e_dot = 2.0386e-4 * sigma^(1.2063) * exp(-295550 / RT)                (s^-1)
#
#  where R is the ideal gas constant taken as 8.314 J/mol-K. In this test a
#  uniaxial specimen is pulled in the axial direction by
#  a varying stress at 1000 K.  Results are taken from different times in the
#  simulation: 2e5, 4e5, 6e5, 8e5 and 1e6 seconds. The creep rate predicted by
#  BISON is compared to the analytical solution for these stresses for the central
#  element in the 2D patch test.
#
#    Time (s)    Von Mises Stress     Analytical Creep       BISON Creep
#                     (MPa)             Rate (s^-1)          Rate (s^-1)
#
#      2e5           17.1306            3.9530e-11           3.9534e-11
#      4e5           33.8741            8.9972e-11           8.9984e-11
#      6e5           50.2010            1.4761e-10           1.4763e-11
#      8e5           66.1119            5.0157e-10           5.0161e-10
#      1e6           81.6204            2.5991e-10           2.5997e-10
#--------------------------------------------------------------------------------

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  coord_type = RZ
  [mesh]
    type = FileMeshGenerator
    file = 1x1_rz.e
  []
[]

[Variables]
  [disp_x]
    order = FIRST
    family = LAGRANGE
  []
  [disp_y]
    order = FIRST
    family = LAGRANGE
  []
  [temp]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300.0
  []
[]

[AuxVariables]
  [vonmises]
    order = CONSTANT
    family = MONOMIAL
  []
  [elastic_strain_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_yy]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [pressure_pull]
    type = PiecewiseLinear
    x = '0 1e6'
    y = '0 1'
  []
  [temperature_function]
    type = PiecewiseLinear
    x = '0 10000'
    y = '300 1000.0'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = finite
  []
[]

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


[AuxKernels]
  [vonmises]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = vonmises
    scalar_type = VonMisesStress
    execute_on = timestep_end
  []
  [elastic_strain_yy]
    type = RankTwoAux
    rank_two_tensor = elastic_strain
    variable = elastic_strain_yy
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  []

  [creep_strain_yy]
    type = RankTwoAux
    rank_two_tensor = creep_strain
    variable = creep_strain_yy
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  []
  [creep_rate_aux]
    type = MaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
[]

[BCs]
  [Pressure]
    [u_top_pull]
      boundary = 4
      factor = 1e8
      function = pressure_pull
    []
  []
  [bottom_fix_y]
    type = DirichletBC
    variable = disp_y
    boundary = 3
    value = 0.0
  []
  [bottom_fix_x]
    type = DirichletBC
    variable = disp_x
    boundary = 3
    value = 0.0
  []
  [temp_bc]
    type = FunctionDirichletBC
    variable = temp
    boundary = '1 2 3 4'
    function = temperature_function
  []
[]

[Materials]
  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1.4e11
    poissons_ratio = 0.17
    block = 1
  []
  [stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'u3si2creep'
    block = 1
  []
  [u3si2creep]
    type = U3Si2CreepUpdate
    block = 1
    temperature = temp
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = 1
    strain_free_density = 11590.0
  []
  [thermal]
    type = HeatConductionMaterial
    block = 1
    thermal_conductivity = 1
    specific_heat = 1
  []
[]

[Executioner]
  type = Transient

  solve_type = 'PJFNK'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -ksp_gmres_restart'
  petsc_options_value = 'lu       superlu_dist                  51'

  line_search = 'none'

  l_max_its  = 100
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-10
  l_tol      = 1e-5
  start_time = 0.0
  end_time   = 1e6

  dt = 1e4
[]

[Postprocessors]
  [creep_rate]
    type = ElementalVariableValue
    variable = creep_rate_aux
    elementid = 2
  []
  [vonmises]
    type = ElementalVariableValue
    variable = vonmises
    elementid = 2
  []
  [fuel_temp]
    type = ElementalVariableValue
    elementid = 2
    variable = temp
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  exodus = true
[]

(examples/accident_tolerant_fuel/u3si2_zircaloy/u3si2_zircaloy.i)


initial_fuel_density = 11590.0

[GlobalParams]
  # Set initial fuel density, other global parameters
  density = ${initial_fuel_density}
  order = SECOND
  family = LAGRANGE
  energy_per_fission = 3.2e-11 # J/fission
  volumetric_locking_correction = false
  displacements = 'disp_x disp_y'
[]

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

[Mesh]
  coord_type = RZ
  # Import mesh file
  patch_size = 10 # For contact algorithm
  patch_update_strategy = auto
  partitioner = centroid
  centroid_partitioner_direction = y
  [mesh]
    type = FileMeshGenerator
    file = u3si2_zircaloy_smeared.e
  []
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [temp]
    initial_condition = 293.0
  []
[]

[UserObjects]
  [pin_geometry]
    type = FuelPinGeometry
    clad_inner_wall = 5
    clad_outer_wall = 2
    clad_top = 3
    clad_bottom = 1
    pellet_exteriors = 8
  []
[]

[AuxVariables]
  # Define auxilary variables
  [fast_neutron_flux]
    block = clad
  []
  [fast_neutron_fluence]
    block = clad
  []
  [gap_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [hoop_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxide_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_hoop_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_rate]
    order = CONSTANT
    family = MONOMIAL
  []
  [densification]
    order = CONSTANT
    family = MONOMIAL
  []
  [solid_swell]
    order = CONSTANT
    family = MONOMIAL
  []
  [gaseous_swell]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 1e4 1e8'
    y = '0 2.5e4 2.5e4'
    scale_factor = 1
  []
  [axial_peaking_factors]
    type = ParsedFunction
    expression = 1
  []
  [pressure_ramp]
    type = PiecewiseLinear
    x = '-200 0 1e8'
    y = '6.537e-3 1 1'
    scale_factor = 15.5e6
  []
  [mass_flux_func]
    type = PiecewiseLinear
    x = '-200 0 1e8'
    y = '3800 3800 3800'
  []
  [q]
    type = CompositeFunction
    functions = 'power_history axial_peaking_factors'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [pellets]
    block = pellet_type_1
    strain = FINITE
    eigenstrain_names = 'fuel_thermal_strain fuel_volumetric_strain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
    extra_vector_tags = 'ref'
  []
  [clad]
    block = clad
    strain = FINITE
    eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_strain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
    extra_vector_tags = 'ref'
  []
[]

[Kernels]
  [gravity]
    type = Gravity
    variable = disp_y
    value = -9.81
    extra_vector_tags = 'ref'
  []
  [heat]
    type = HeatConduction
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temp
    extra_vector_tags = 'ref'
  []
  [heat_source]
    type = NeutronHeatSource
    variable = temp
    block = pellet_type_1
    burnup_function = burnup
    extra_vector_tags = 'ref'
  []
[]

[Burnup]
  [burnup]
    block = pellet_type_1
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    num_radial = 81
    num_axial = 11
    fuel_pin_geometry = pin_geometry
    fuel_volume_ratio = 1.0
    RPF = RPF
    fuel_type = U3Si2
  []
[]

[AuxKernels]
  [fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    block = clad
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    factor = 3e13
    execute_on = timestep_begin
  []
  [fast_neutron_fluence]
    type = FastNeutronFluenceAux
    variable = fast_neutron_fluence
    block = clad
    fast_neutron_flux = fast_neutron_flux
    execute_on = timestep_begin
  []
  [hoop_stress]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = hoop_stress
    scalar_type = HoopStress
    execute_on = timestep_end
  []
  [total_hoop_strain]
    type = RankTwoScalarAux
    rank_two_tensor = total_strain
    variable = total_hoop_strain
    scalar_type = HoopStress
    execute_on = timestep_end
  []
  [conductance]
    type = MaterialRealAux
    property = gap_conductance
    variable = gap_cond
    boundary = 10
  []
  [coolant_htc]
    type = MaterialRealAux
    property = coolant_channel_htc
    variable = coolant_htc
    boundary = 2
  []
  [oxide]
    type = MaterialRealAux
    variable = oxide_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [creep_rate]
    type = MaterialRealAux
    variable = creep_rate
    property = creep_rate
    execute_on = timestep_end
    block = clad
  []
  [densfication]
    type = MaterialRealAux
    property = densification
    variable = densification
    block = pellet_type_1
  []
  [solid_swell]
    type = MaterialRealAux
    property = solid_swelling
    variable = solid_swell
    block = pellet_type_1
  []
  [gaseous_swell]
    type = MaterialRealAux
    property = gaseous_swelling
    variable = gaseous_swell
    block = pellet_type_1
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    normal_smoothing_distance = 0.1
    penalty = 1e7
  []
[]

[ThermalContact]
  [thermal_contact]
    type = GasGapHeatTransfer
    variable = temp
    primary = 5
    secondary = 10
    initial_moles = initial_moles
    gas_released = fis_gas_released
    contact_pressure = contact_pressure
    quadrature = true
  []
[]

[BCs]
  [no_x_all] # pin pellets and clad along axis of symmetry (y)
    type = DirichletBC
    variable = disp_x
    boundary = 12
    value = 0.0
  []
  [no_y_clad_bottom] # pin clad bottom in the axial direction (y)
    type = DirichletBC
    variable = disp_y
    boundary = '1'
    value = 0.0
  []
  [no_y_fuel_bottom] # pin fuel bottom in the axial direction (y)
    type = DirichletBC
    variable = disp_y
    boundary = '1020'
    value = 0.0
  []
  [Pressure]
    [coolantPressure]
      boundary = '1 2 3'
      function = pressure_ramp
    []
  []
  [PlenumPressure]
    [plenumPressure]
      boundary = 9
      initial_pressure = 2.0e6
      startup_time = 0
      R = 8.3143
      output_initial_moles = initial_moles
      temperature = ave_temp_interior
      volume = gas_volume
      material_input = fis_gas_released
      output = plenum_pressure
    []
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = '1 2 3'
    variable = temp
    inlet_temperature = 580 # K
    inlet_pressure = pressure_ramp # Pa
    inlet_massflux = mass_flux_func # kg/m^2-sec
    rod_diameter = 9.4996e-3 # m
    rod_pitch = 1.26e-2 # m
    linear_heat_rate = power_history
    axial_power_profile = axial_peaking_factors
  []
[]

[Materials]
  [fuel_thermal]
    type = SilicideFuelThermal
    block = pellet_type_1
    thermal_conductivity_model = WHITE
    silicon_mole_fraction = 0.4
    temperature = temp
  []
  [fuel_elasticity_tensor]
    type = U3Si2ElasticityTensor
    block = pellet_type_1
  []
  [fuel_stress]
    type = ComputeMultipleInelasticStress
    block = pellet_type_1
    tangent_operator = elastic
    inelastic_models = 'fuel_creep'
  []
  [fuel_creep]
    type = U3Si2CreepUpdate
    block = pellet_type_1
    temperature = temp
  []
  [fuel_thermal_expansion]
    type = U3Si2ThermalExpansionEigenstrain
    block = pellet_type_1
    temperature = temp
    stress_free_temperature = 295.0
    eigenstrain_name = fuel_thermal_strain
  []
  [fuel_volumetric_swelling]
    type = U3Si2VolumetricSwellingEigenstrain
    block = pellet_type_1
    gaseous_swelling_type = U3SI2FG
    temperature = temp
    burnup_function = burnup
    eigenstrain_name = fuel_volumetric_strain
  []

  [clad_thermal]
    type = ZryThermal
    temperature = temp
    block = clad
  []
  [clad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 7.5e10
    poissons_ratio = 0.3
    block = clad
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep clad_plasticity'
    relative_tolerance = 1e-5
    block = clad
  []
  [clad_zrycreep]
    type = ZryCreepLimbackHoppeUpdate
    block = clad
    temperature = temp
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    model_irradiation_creep = true
    model_primary_creep = true
    model_thermal_creep = true
    relative_tolerance = 1e-5
    max_inelastic_increment = 1e-4
    zircaloy_material_type = stress_relief_annealed
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = clad
    temperature = temp
    stress_free_temperature = 295.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [irradiation_swelling]
    type = ZryIrradiationGrowthEigenstrain
    block = clad
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = stress_relief_annealed
    eigenstrain_name = clad_irradiation_strain
  []
  [clad_plasticity]
    type = ZryPlasticityUpdate
    block = clad
    temperature = temp
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    relative_tolerance = 1e-5
    cold_work_factor = 0.5
    plasticity_model_type = MATPRO
    zircaloy_alloy_type = 4
  []
  [fission_gas_behavior]
    type = U3Si2Sifgrs
    block = pellet_type_1
    temperature = temp
    burnup_function = burnup
    saturation_coverage = 0.5
  []
  [clad_density]
    type = StrainAdjustedDensity
    block = clad
    strain_free_density = 6511.0
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = pellet_type_1
    strain_free_density = ${initial_fuel_density}
  []
  [ZryOxidation]
    type = ZryOxidation
    boundary = 2
    clad_inner_radius = 4.1783e-3
    clad_outer_radius = 4.7498e-3
    normal_operating_temperature_model = epri_kwu_ce
    temperature = temp
    fast_neutron_flux = fast_neutron_flux
    use_coolant_channel = true
    oxygen_weight_fraction_initial = 0.0012
  []
[]

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

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -ksp_gmres_restart'
  petsc_options_value = 'lu       superlu_dist                  51'

  line_search = 'none'

  l_max_its = 100
  l_tol = 8e-3

  nl_max_its = 25
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-10

  start_time = -200
  n_startup_steps = 1
  end_time = 1e8

  dtmax = 1e6
  dtmin = 1e-3

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 2.0e2
    force_step_every_function_point = true
    timestep_limiting_function = power_history
    max_function_change = 3e20
    optimal_iterations = 10
    iteration_window = 2
    linear_iteration_ratio = 100
    timestep_limiting_postprocessor = material_timestep
  []
  [Quadrature]
    order = FIFTH
    side_order = SEVENTH
  []
[]

[Postprocessors]
  [avg_fuel_surface]
    type = SideAverageValue
    boundary = 10
    variable = temp
  []
  [ave_temp_interior]
    type = SideAverageValue
    boundary = 9
    variable = temp
    execute_on = 'initial linear'
  []
  [clad_inner_vol]
    type = InternalVolume
    boundary = 7
  []
  [pellet_volume]
    type = InternalVolume
    boundary = 8
  []
  [avg_clad_temp]
    type = SideAverageValue
    boundary = 7
    variable = temp
  []
  [fis_gas_produced]
    type = ElementIntegralFisGasGeneratedSifgrs
    block = pellet_type_1
  []
  [fis_gas_released]
    type = ElementIntegralFisGasReleasedSifgrs
    block = pellet_type_1
  []
  [fis_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    block = pellet_type_1
  []
  [fis_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    block = pellet_type_1
  []
  [gas_volume]
    type = InternalVolume
    boundary = 9
    execute_on = 'initial linear'
  []
  [flux_from_clad]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 5
    diffusivity = thermal_conductivity
  []
  [flux_from_fuel]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 10
    diffusivity = thermal_conductivity
  []
  [_dt]
    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 = ElementIntegralPower
    variable = temp
    burnup_function = burnup
    block = pellet_type_1
  []
  [rod_input_power]
    type = FunctionValuePostprocessor
    function = power_history
    scale_factor = 0.1186
  []
  [average_burnup]
    type = ElementAverageValue
    block = pellet_type_1
    variable = burnup
  []
  [oxide_thickness]
    type = ElementExtremeValue
    block = clad
    variable = oxide_thickness
  []
  [fis_gas_percent]
    type = FGRPercent
    fission_gas_released = fis_gas_released
    fission_gas_generated = fis_gas_produced
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = clad
  []
[]

[Outputs]
  perf_graph = true
  time_step_interval =  1
  exodus = true
  color = false
  csv = true
  print_linear_residuals = true
  [console]
    type = Console
    max_rows = 25
  []
[]

Description
Freeman Model
Metzger Model
Yingling Model
Example Input Syntax
Input Parameters
Input Files
References