FeCrAlThermal

Computes the specific heat and thermal conductivity for FeCrAl alloys.

Description

The FeCrAlThermal model computes the specific heat and thermal conductivity for a variety of FeCrAl alloys. The alloys to available include the commercial alloys Special Metals Incoloy MA956, Plansee PM2000, Kanthal APMT, Resistalloy International Fecralloy, and Oak Ridge National Laboratory alloys C06M, C35M, and C36M.

The recent FeCrAl Handbook (Field et al., 2018) contains temperature dependent correlations for thermal conductivity and specific heat of four of the alloys listed above: Kanthal APMT, C06M, C35M, and C36M. The thermal conductivity is given as:

$(1)var element = document.getElementById("moose-equation-d058928f-276e-4919-bafc-46e2002f6d6c");katex.render(" k = A_1T^2 + A_2T + A_3", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-df44c8a4-7b83-43fc-b4de-56e11ce541e7");katex.render("A_1", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-62bbcd67-688a-48f4-a979-f164117eaf86");katex.render("A_2", element, {displayMode:false,throwOnError:false});$ , and $var element = document.getElementById("moose-equation-b34936d9-b112-499a-acfd-309b3c464bb9");katex.render("A_3", element, {displayMode:false,throwOnError:false});$ are fitting constants tabulated in Table 1, and $var element = document.getElementById("moose-equation-0ce70e04-2afb-42df-9eff-3d341fe4d910");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the temperature in K. Note, the mean values provided in the handbook are tabulated.

Table 1: Coefficients used in the polynomial fit for thermal conductivity of select FeCrAl alloys

Alloy	A $var element = document.getElementById("moose-equation-10c1710a-f9fe-4e08-965d-9a529f8be2cb");katex.render("_1", element, {displayMode:false,throwOnError:false});$ ( $var element = document.getElementById("moose-equation-ed75be46-da19-4a6e-982a-0d21552348cd");katex.render("\\times 10^{-7}", element, {displayMode:false,throwOnError:false});$ )	A $var element = document.getElementById("moose-equation-a0d6481d-4ffc-4495-9bf0-47615e5374e0");katex.render("_2", element, {displayMode:false,throwOnError:false});$ ( $var element = document.getElementById("moose-equation-718676d6-5e2e-43c8-98b8-11e9b15c3bf1");katex.render("\\times 10^{-2}", element, {displayMode:false,throwOnError:false});$ )	A $var element = document.getElementById("moose-equation-fcf7f023-0b56-4204-aedc-97be7cde48d2");katex.render("_3", element, {displayMode:false,throwOnError:false});$
Kanthal APMT	-7.223	1.5628	6.569
C06M	6.762	1.032	9.956
C35M	-19.86	1.537	8.502
C36M	-9.184	1.368	8.187

The specific heat capacity correlations provided in the handbook are divided into two temperature regimes. Below the Curie Temperature ( $var element = document.getElementById("moose-equation-cc15834f-80e0-43a7-b7fb-882f744bf0ba");katex.render("T_c", element, {displayMode:false,throwOnError:false});$ ) the specific heat correlation is given by:

$var element = document.getElementById("moose-equation-76baa93b-2ad3-4346-b72a-54a3efcec6c4");katex.render("C_p = AT + BT^2 + CT^3", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-7c4cb56d-6c1f-460a-bd96-667a517010cf");katex.render("A", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-84c17256-c66b-43c3-9a97-2b643c0f10d3");katex.render("B", element, {displayMode:false,throwOnError:false});$ , and $var element = document.getElementById("moose-equation-71afe2e7-14ae-4d35-bd80-7fcddb441eb3");katex.render("C", element, {displayMode:false,throwOnError:false});$ are fitting constants, $var element = document.getElementById("moose-equation-5f33f435-b581-4262-9d6f-40b37bd3d120");katex.render("C_p", element, {displayMode:false,throwOnError:false});$ is the specific heat capacity in J/kg-K and $var element = document.getElementById("moose-equation-86280e54-42ae-4f22-82ec-3aa3a1052eac");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the temperature in K. The fitting constants as well as the Curie temperature for the alloys are summarized in Table 2. It should be mentioned that the units provided for these coefficients here are different than the handbook. It appears that the handbook was incorrectly assuming $var element = document.getElementById("moose-equation-967f0d70-168e-4125-b425-147ac9afee3d");katex.render("T", element, {displayMode:false,throwOnError:false});$ in units of K $var element = document.getElementById("moose-equation-a16224b8-1a19-4370-8239-647ed00a17d4");katex.render("^{-1}", element, {displayMode:false,throwOnError:false});$ when deriving the units of the coefficients.

Table 2: Coefficients used in the polynomial fit for specific heat of select alloys below the Curie temperature

Alloy	A (J/kg-K $var element = document.getElementById("moose-equation-713ab2c1-28ba-427c-a95d-dc142fc273e2");katex.render("^{2}", element, {displayMode:false,throwOnError:false});$ )	B (J/kg-K $var element = document.getElementById("moose-equation-66b50f68-b3d1-42d4-a760-f10d0fb5d547");katex.render("^{3}", element, {displayMode:false,throwOnError:false});$ $var element = document.getElementById("moose-equation-6c493146-c999-4301-810b-d433eb8fbe67");katex.render("\\times 10^{-3}", element, {displayMode:false,throwOnError:false});$ )	C (J/kg-K $var element = document.getElementById("moose-equation-e13ccdf5-438f-4f87-bfce-cf76361e55b1");katex.render("^{4}", element, {displayMode:false,throwOnError:false});$ $var element = document.getElementById("moose-equation-cf741da7-8ee2-465a-875e-477252d3cf61");katex.render("\\times 10^{-6}", element, {displayMode:false,throwOnError:false});$ )	T $var element = document.getElementById("moose-equation-83f8a2a5-3ec5-448f-b2ab-b75615467006");katex.render("_c", element, {displayMode:false,throwOnError:false});$ (K)
Kanthal APMT	2.540	-4.311	2.982	852
C06M	2.430	-3.957	2.656	888
C35M	2.450	-4.002	2.720	870
C36M	2.995	-5.953	4.516	771

Above the Curie temperature the coefficients have additional terms to capture the larger increase in specific heat in the 750-900 K range. This increase is attributed to a phase change in the alloys from ferromagnetic to paramagnetic. The correlation is given by:

$var element = document.getElementById("moose-equation-9808fa02-f90a-415e-8e3d-50c418625e3d");katex.render("C_p = AT + BT^2 + CT^3 + \\frac{D}{T} + E \\ln \\left( \\frac{\\left | \\left(T - T_c\\right) \\right |}{T_c} \\right)", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-58955d9e-2c74-4fe1-9770-007db1bc4925");katex.render("D", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-cc13f99a-a069-4cec-a5a5-5c5fb74822b9");katex.render("E", element, {displayMode:false,throwOnError:false});$ are additional fitting constants. The fitting constants in this temperature regime are summarized in Table 3.

Table 3: Coefficients used in the polynomial fit for specific heat of select alloys above the Curie temperature

Alloy	A (J/kg-K $var element = document.getElementById("moose-equation-49ec0258-b9df-445f-bbda-b50e7f430c4c");katex.render("^{2}", element, {displayMode:false,throwOnError:false});$ )	B (J/kg-K $var element = document.getElementById("moose-equation-ff5e7333-f63c-465a-9e0a-8664f30468df");katex.render("^{3}", element, {displayMode:false,throwOnError:false});$ $var element = document.getElementById("moose-equation-6d99b93c-8009-42e5-801f-06f99c2251a6");katex.render("\\times 10^{-3}", element, {displayMode:false,throwOnError:false});$ )	C (J/kg-K $var element = document.getElementById("moose-equation-03b92cc5-ba75-492b-b402-fca00ebf8170");katex.render("^{4}", element, {displayMode:false,throwOnError:false});$ $var element = document.getElementById("moose-equation-ccd036b8-6d4f-4ea7-a045-1c79fd9c9cfc");katex.render("\\times 10^{-6}", element, {displayMode:false,throwOnError:false});$ )	D (J/kg $var element = document.getElementById("moose-equation-e4fa49e6-19da-4b5f-b148-12aafe97fa69");katex.render("\\times 10^{3}", element, {displayMode:false,throwOnError:false});$ )	E (J/kg-K)
Kanthal APMT	1.840	-1.843	0.643	-5.712	-50.38
C06M	1.827	-1.807	0.6134	-9.419	-54.54
C35M	1.946	-2.002	0.698	-1.652	-53.93
C36M	1.456	-1.296	0.438	26.45	-46.89

For the commercial alloys Special Metals MA956 and Plansee PM2000 the temperature dependent thermal conductivity and specific heat capacity are given in tabular form in their data sheets. Since no additional information is given about the behavior as a function of temperature, the material properties are linearly interpolated between the values provided in the table. For temperatures outside of the range provided in the table the thermal conductivity or specific heat is taken as the closest known value to avoid extrapolation into areas where no data is known.

The tabulated data for MA956 (Corporation, 2004) and PM2000 (MatWeb, 2014) are reproduced in Table 4 and Table 5 respectively. The thermal conductivity and specific heat capacity of Fecralloy are given as constant values of 16.0 W/m-K and 460 J/kg-K respectively (MatWeb, 2014).

Table 4: Temperature dependent thermal conductivity and specific heat capacity of MA956 alloy

Temperature (K)	Thermal Conductivity (W/m-K)	Specific Heat Capacity (J/kg-K)
293.15	10.9	469
373.15	12.2	491
473.15	13.9	519
573.15	15.4	547
673.15	16.9	575
773.15	18.4	608
873.15	19.8	630
973.15	21.2	658
1073.15	22.6	686
1173.15	24.1	714
1273.15	25.5	741
1373.15	27	769

$var element = document.getElementById("moose-equation-be9ce7e0-4a59-4875-add3-08826b26c58a");katex.render("~", element, {displayMode:false,throwOnError:false});$

Table 5: Temperature dependent thermal conductivity and specific heat capacity of PM2000 alloy

Temperature (K)	Thermal Conductivity (W/m-K)	Specific Heat Capacity (J/kg-K)
293.15	10.9
373.15		500
473.15	16	480
773.15	21	610
1073.15	22	680
1273.15	25.5	740
1473.15	28

$var element = document.getElementById("moose-equation-a6d83eff-70e0-4f93-95ac-1bf8ccfdd210");katex.render("~", element, {displayMode:false,throwOnError:false});$

Range of Applicability and Uncertainty

All of the thermal conductivity and specific heat models have an associated range of temperature over which they are valid. The user is cautioned to examine their results when using the correlations or tabulated values outside their applicability ranges. The correlations for Kanthal APMT, C06M, C35M, and C36M are applicable from 300 K to $var element = document.getElementById("moose-equation-7f01f3a4-8342-4d73-a0fa-9fbf8b230760");katex.render("\\sim", element, {displayMode:false,throwOnError:false});$ 1400 K. Recall that the specific heat correlation has two different equations depending upon whether the temperature is above or below the Curie temperature. For the MA956 and PM2000 alloys, which only have tabulated values, the closest known value is taken for thermal conductivity or specific heat for temperatures outside the available range to avoid unknown consequences of extrapolation. The correlations for MA956 are applicable from 293.15 K to 1373.15 K. For PM2000, the thermal conductivity correlation is applicable from 293.15 K to 1473.15 K whereas the specific heat correlation is applicable from 373.15 K to 1473.15 K.

The correlations for Kanthal APMT, C06M, C35M, and C36M from the FeCrAl handbook (Field et al., 2018) are fit to experimental data that has some associated uncertainty. For thermal conductivity, the correlations all have R $var element = document.getElementById("moose-equation-3ff1bf85-20bc-4f32-bd44-f2cbd646a27b");katex.render("^2", element, {displayMode:false,throwOnError:false});$ values greater than 0.96 with the experimental data. Thus, the uncertainty on the experimental data is applied to the models for thermal conductivity. This uncertainty over the entire temperature range of applicability is 7%. For specific heat currently no uncertainty in the experimental data has been provided.

Example Input Syntax

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [fuel_thermalFeCrAl]
    type = FeCrAlThermal<<<{"description": "Computes the specific heat and thermal conductivity for FeCrAl alloys.", "href": "FeCrAlThermal.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 0
    material<<<{"description": "The FeCrAl alloy being used for the cladding. Choices are: APMT C06M C35M C36M MA956 PM2000 FECRALLOY"}>>> = APMT
    scale_factor_k<<<{"description": "Scale factor to be applied to the thermal conductivity"}>>> = 1.0
    scale_factor_cp<<<{"description": "Scale factor to be applied to the specific heat"}>>> = 1.0
    temperature<<<{"description": "Coupled temperature"}>>> = T
  []
[]

Input Parameters

temperatureCoupled temperature
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Coupled temperature

Required Parameters

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.
materialAPMTThe FeCrAl alloy being used for the cladding. Choices are: APMT C06M C35M C36M MA956 PM2000 FECRALLOY
Default:APMT
C++ Type:MooseEnum
Options:APMT, C06M, C35M, C36M, MA956, PM2000, FECRALLOY
Controllable:No
Description:The FeCrAl alloy being used for the cladding. Choices are: APMT C06M C35M C36M MA956 PM2000 FECRALLOY
scale_factor_cp1Scale factor to be applied to the specific heat
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Scale factor to be applied to the specific heat
scale_factor_k1Scale factor to be applied to the thermal conductivity
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Scale factor to be applied to the thermal conductivity

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/thermalFeCrAl/thermalFeCrAl_APMT.i)
(test/tests/solid_mechanics/fecral_elasticity_tensor/elasticity_MA956.i)
(test/tests/thermalFeCrAl/thermalFeCrAl_Fecralloy.i)
(test/tests/solid_mechanics/fecral_elasticity_tensor/elasticity_Fecralloy.i)
(test/tests/fecral_oxidation/corrosion_test_fecral_pwr.i)
(test/tests/thermalFeCrAl/thermalFeCrAl_MA956.i)
(test/tests/solid_mechanics/fecral_elasticity_tensor/elasticity_PM2000.i)
(test/tests/solid_mechanics/fecral_elasticity_tensor/elasticity_APMT.i)
(test/tests/thermalFeCrAl/thermalFeCrAl_C36M.i)
(examples/accident_tolerant_fuel/uo2_fecral/uo2_fecral.i)
(test/tests/solid_mechanics/fecral_elasticity_tensor/elasticity_C06M.i)
(test/tests/solid_mechanics/fecral_eigenstrains/fecral_vswelling/swelling.i)
(test/tests/fecral_oxidation/corrosion_test_fecral_high_temp.i)
(test/tests/thermalFeCrAl/thermalFeCrAl_C06M.i)
(test/tests/thermalFeCrAl/thermalFeCrAl_PM2000.i)
(test/tests/thermalFeCrAl/thermalFeCrAl_C35M.i)
(test/tests/fecral_oxidation/corrosion_test_fecral_bwr.i)
(test/tests/fecral_oxidation/corrosion_test_fecral.i)

References

Special Metals Corporation. Special Metals Incoloy alloy MA956. www.specialmetals.com/documents/Incoloy+alloy+MA956.pdf, 2004.

@MISC{ma956_2004,
    author = "Corporation, Special Metals",
    title = "{Special Metals Incoloy alloy MA956}",
    howpublished = "www.specialmetals.com/documents/Incoloy+alloy+MA956.pdf",
    year = "2004"
}

K. G. Field, M. A. Snead, Y. Yamamoto, and K. A. Terrani. Handbook on the material properties of fecral alloys for nuclear power production applications. Technical Report ORNL/SPR-2018/905 Rev. 1, Oak Ridge National Laboratory, 2018.

@TECHREPORT{fecral_handbook,
    author = "Field, K. G. and Snead, M. A. and Yamamoto, Y. and Terrani, K. A.",
    title = "Handbook on the Material Properties of FeCrAl Alloys for Nuclear Power Production Applications",
    institution = "Oak Ridge National Laboratory",
    number = "ORNL/SPR-2018/905 Rev. 1",
    year = "2018"
}

MatWeb. Resistalloy International Fecralloy Electrical Resistance Steel. http://www.matweb.com/search/datasheet.aspx?MatGUID=c2427c6297594858bedac2a4e5981d2f, 2014.

@MISC{fecralloy_matweb,
    author = "MatWeb",
    title = "{Resistalloy International Fecralloy Electrical Resistance Steel}",
    howpublished = "http://www.matweb.com/search/datasheet.aspx?MatGUID=c2427c6297594858bedac2a4e5981d2f",
    year = "2014"
}

MatWeb. Schwarzkopf Plansee PM 2000. http://www.matweb.com/search/datasheet.aspx?matguid=21e9ec9a0de24b47bcf69ab11c375567, 2014.

@MISC{pm2000_matweb,
    author = "MatWeb",
    title = "{Schwarzkopf Plansee PM 2000}",
    howpublished = "http://www.matweb.com/search/datasheet.aspx?matguid=21e9ec9a0de24b47bcf69ab11c375567",
    year = "2014"
}

(test/tests/thermalFeCrAl/thermalFeCrAl_APMT.i)

# The mesh is a 1x1x1 cube (single element).
# The temperature is ramped on all faces of the cube from 300 K to 1500K
# over 100 Ms. The cladding alloy is Kanthal APMT.
#
# The thermal conductivity and specific heat computed by BISON
# for the temperatures shown are compared with the analytical solution below.
# The results are as follows:
#
# Temp(K)    k BISON     k analytical      cp BISON      cp analytical
#            (W/m-K)       (W/m-K)         (J/kg-K)         (J/kg-K)
#
#   300     11.19239      11.19239         454.5240         454.5240
#   408     12.82499      12.82499         521.2231         521.2231
#   504     14.26204      14.26204         566.8648         566.8648
#   600     15.68577      15.68577         616.1520         616.1520
#   708     17.27156      17.27156         695.6675         695.6675
#   804     18.66701      18.66701         805.2611         805.2611
#   900     20.04914      20.04914         770.4826         770.4826
#  1008     21.58812      21.58812         720.5348         720.5348
#  1104     22.94196      22.94196         706.4823         706.4823
#  1200     24.28249      24.28249         705.5334         705.5334
#  1308     25.77467      25.77467         719.6348         719.6348
#  1404     27.08690      27.08690         747.7660         747.7660
#  1500     28.38583      28.38583         793.3558         793.3558
#
#------------------------------------------------------------------------
[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 3
  []
[]

[Variables]
  [T]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300     # set initial T to 500 K
  []
[]

[AuxVariables]
  [th_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [specific_heat]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [th_cond]
    type = MaterialRealAux
    variable = th_cond
    property = thermal_conductivity
    execute_on = 'initial linear'
    block = 0
  []
  [specific_heat]
    type = MaterialRealAux
    variable = specific_heat
    property = specific_heat
    execute_on = 'initial linear'
    block = 0
  []
[]

[Functions]
  [temp_ramp]
    type = PiecewiseLinear
    x = '0.0 100.0e6'
    y = '300 1500'
  []
[]

[BCs]
   [VariableT]
     type = FunctionDirichletBC
     boundary = 'top bottom front back left right'      # All cube faces
     variable = T
     function = temp_ramp
   []
[]

[Materials]
  [fuel_thermalFeCrAl]
    type = FeCrAlThermal
    block = 0
    material = APMT
    scale_factor_k = 1.0
    scale_factor_cp = 1.0
    temperature = T
  []
  [density]
    type = ParsedMaterial
    block = 0
    property_name = density
    expression = 7250.0
  []
[]

[Executioner]
   type = Transient

   solve_type = PJFNK

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

   line_search = 'none'

   l_max_its = 100
   nl_max_its = 100
   nl_rel_tol = 1e-4
   nl_abs_tol = 1e-6
   l_tol = 1e-5
   start_time = 0.0
   num_steps = 100
   dt = 1.0e6
[]

[Postprocessors]
  [ave_temp_interior]
    type = SideAverageValue
    boundary = 'top bottom front back left right'
    variable = T
    execute_on = 'initial linear'
  []
  [avg_th_cond]
    type = ElementAverageValue
    block = 0
    variable = th_cond
    execute_on = 'initial linear'
  []
  [average_clad_T]
    type = ElementAverageValue
    block = 0
    variable = T
    execute_on = 'initial linear'
  []
  [avg_specific_heat]
    type = ElementAverageValue
    block = 0
    variable = specific_heat
    execute_on = 'initial linear'
  []
[]

[Outputs]
  [out]
    type = Exodus
  []
[]

(test/tests/thermalFeCrAl/thermalFeCrAl_APMT.i)

# The mesh is a 1x1x1 cube (single element).
# The temperature is ramped on all faces of the cube from 300 K to 1500K
# over 100 Ms. The cladding alloy is Kanthal APMT.
#
# The thermal conductivity and specific heat computed by BISON
# for the temperatures shown are compared with the analytical solution below.
# The results are as follows:
#
# Temp(K)    k BISON     k analytical      cp BISON      cp analytical
#            (W/m-K)       (W/m-K)         (J/kg-K)         (J/kg-K)
#
#   300     11.19239      11.19239         454.5240         454.5240
#   408     12.82499      12.82499         521.2231         521.2231
#   504     14.26204      14.26204         566.8648         566.8648
#   600     15.68577      15.68577         616.1520         616.1520
#   708     17.27156      17.27156         695.6675         695.6675
#   804     18.66701      18.66701         805.2611         805.2611
#   900     20.04914      20.04914         770.4826         770.4826
#  1008     21.58812      21.58812         720.5348         720.5348
#  1104     22.94196      22.94196         706.4823         706.4823
#  1200     24.28249      24.28249         705.5334         705.5334
#  1308     25.77467      25.77467         719.6348         719.6348
#  1404     27.08690      27.08690         747.7660         747.7660
#  1500     28.38583      28.38583         793.3558         793.3558
#
#------------------------------------------------------------------------
[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 3
  []
[]

[Variables]
  [T]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300     # set initial T to 500 K
  []
[]

[AuxVariables]
  [th_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [specific_heat]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [th_cond]
    type = MaterialRealAux
    variable = th_cond
    property = thermal_conductivity
    execute_on = 'initial linear'
    block = 0
  []
  [specific_heat]
    type = MaterialRealAux
    variable = specific_heat
    property = specific_heat
    execute_on = 'initial linear'
    block = 0
  []
[]

[Functions]
  [temp_ramp]
    type = PiecewiseLinear
    x = '0.0 100.0e6'
    y = '300 1500'
  []
[]

[BCs]
   [VariableT]
     type = FunctionDirichletBC
     boundary = 'top bottom front back left right'      # All cube faces
     variable = T
     function = temp_ramp
   []
[]

[Materials]
  [fuel_thermalFeCrAl]
    type = FeCrAlThermal
    block = 0
    material = APMT
    scale_factor_k = 1.0
    scale_factor_cp = 1.0
    temperature = T
  []
  [density]
    type = ParsedMaterial
    block = 0
    property_name = density
    expression = 7250.0
  []
[]

[Executioner]
   type = Transient

   solve_type = PJFNK

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

   line_search = 'none'

   l_max_its = 100
   nl_max_its = 100
   nl_rel_tol = 1e-4
   nl_abs_tol = 1e-6
   l_tol = 1e-5
   start_time = 0.0
   num_steps = 100
   dt = 1.0e6
[]

[Postprocessors]
  [ave_temp_interior]
    type = SideAverageValue
    boundary = 'top bottom front back left right'
    variable = T
    execute_on = 'initial linear'
  []
  [avg_th_cond]
    type = ElementAverageValue
    block = 0
    variable = th_cond
    execute_on = 'initial linear'
  []
  [average_clad_T]
    type = ElementAverageValue
    block = 0
    variable = T
    execute_on = 'initial linear'
  []
  [avg_specific_heat]
    type = ElementAverageValue
    block = 0
    variable = specific_heat
    execute_on = 'initial linear'
  []
[]

[Outputs]
  [out]
    type = Exodus
  []
[]

(test/tests/solid_mechanics/fecral_elasticity_tensor/elasticity_MA956.i)

# This test case is prepared to test the elastic moduli of MA956 FeCrAl alloy.
#
# The geometry represents a ring of MA956 with inner radius of 0.005 m,
# outer radius of 0.0055 m, and a height of 0.01 m.
#
# An axial load is applied to the top of the ring with a value of 100 MPa.
# The strain in the radial direction can be calculated using the Young's modulus
# and the Poisson's ratio from the applied 100.0 MPa pressure.
#
# epsilon_{rr} = nu * sigma_{axial} / E
#
# The Young's modulus (E) for MA956 is calculated via a piecewse linear function.
#
# The Poisson's ratio of MA956 is calculated is taken as 0.3.
#
# The results of BISON compared to the analytical solution for an increase in
# temperature under constant pressure are summarized:
#
#                          BISON             Analytical
#    Temperature (K)    strain_rr (m/m)    strain_rr (m/m)
#        420              1.794104e-4        1.794104e-4
#        540              1.900288e-4        1.900288e-4
#        660              2.057089e-4        2.057089e-4
#        780              2.202522e-4        2.202522e-4
#        900              2.477721e-4        2.477721e-4
#       1020              2.703530e-4        2.703530e-4
#       1140              2.999621e-4        2.999621e-4
#       1260              3.807386e-4        3.807386e-4
#       1380              4.126547e-4        4.126547e-4
#       1500              4.126547e-4        4.126547e-4
#

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

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

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

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = FINITE
        block = 0
        add_variables = true
        generate_output = 'elastic_strain_xx'
      []
    []
  []
[]

[Functions]
  [temperature_function]
    type = PiecewiseLinear
    x = '0       10'
    y = '300   1500'
  []
  [pressure_function]
    type = PiecewiseLinear
    x = '0  10'
    y = '1   1'
  []
[]

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

[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
  []

  [temp_bc]
    type = FunctionDirichletBC
    variable = temperature
    boundary = 'bottom top left right'
    function = temperature_function
  []
[]

[Materials]
  [elasticity_tensor]
    type = FeCrAlElasticityTensor
    temperature = temperature
    fecral_material_type = MA956
  []
  [stress]
    type = ComputeFiniteStrainElasticStress
  []
  [clad_density]
    type = StrainAdjustedDensity
    strain_free_density = 1.0
  []
  [thermal]
    type = FeCrAlThermal
    temperature = temperature
    material = MA956
  []
[]

[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   = 10
  dt         = 1
[]

[Postprocessors]
  [elastic_strain_rr]
    type = ElementAverageValue
    variable = elastic_strain_xx
  []
  [temp]
    type = AverageNodalVariableValue
    variable = temperature
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/thermalFeCrAl/thermalFeCrAl_Fecralloy.i)

# The mesh is a 1x1x1 cube (single element).
# The temperature is ramped on all faces of the cube from 300 K to 1500K
# over 100 Ms. The cladding alloy is Resistalloy International Fecralloy.
#
# The thermal conductivity and specific heat computed by BISON
# for the temperatures shown are compared with the analytical solution below.
# The results are as follows:
#
# Temp(K)    k BISON     k analytical      cp BISON      cp analytical
#            (W/m-K)       (W/m-K)         (J/kg-K)         (J/kg-K)
#     16
#   300        16            16              460              460
#   408        16            16              460              460
#   504        16            16              460              460
#   600        16            16              460              460
#   708        16            16              460              460
#   804        16            16              460              460
#   900        16            16              460              460
#  1008        16            16              460              460
#  1104        16            16              460              460
#  1200        16            16              460              460
#  1308        16            16              460              460
#  1404        16            16              460              460
#  1500        16            16              460              460
#
#------------------------------------------------------------------------
[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 3
  []
[]

[Variables]
  [T]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300     # set initial T to 500 K
  []
[]

[AuxVariables]
  [th_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [specific_heat]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = T
  []
[]

[AuxKernels]
  [th_cond]
    type = MaterialRealAux
    variable = th_cond
    property = thermal_conductivity
    block = 0
    execute_on = 'initial linear'
  []
  [specific_heat]
    type = MaterialRealAux
    variable = specific_heat
    property = specific_heat
    block = 0
    execute_on = 'initial linear'
  []
[]

[Functions]
  [temp_ramp]
    type = PiecewiseLinear
    x = '0.0 100.0e6'
    y = '300 1500'
  []
[]

[BCs]
   [VariableT]
     type = FunctionDirichletBC
     boundary = 'top bottom front back left right'       # All cube faces
     variable = T
     function = temp_ramp
   []
[]

[Materials]
  [fuel_thermalFeCrAl]
    type = FeCrAlThermal
    block = 0
    material = FECRALLOY
    scale_factor_k = 1.0
    scale_factor_cp = 1.0
    temperature = T
  []
  [density]
    type = ParsedMaterial
    block = 0
    property_name = density
    expression = 7250.0
  []
[]

[Executioner]
  type = Transient

  solve_type = PJFNK

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

  line_search = 'none'

  l_max_its = 100
  nl_max_its = 100
  nl_rel_tol = 1e-4
  nl_abs_tol = 1e-6
  l_tol = 1e-5
  start_time = 0.0
  num_steps = 100
  dt = 1.0e6
[]

[Postprocessors]
  [ave_temp_interior]
    type = SideAverageValue
    boundary = 'top bottom front back left right'
    variable = T
    execute_on = 'initial linear'
  []
  [avg_th_cond]
    type = ElementAverageValue
    block = 0
    variable = th_cond
    execute_on = 'initial linear'
  []
  [average_clad_T]
    type = ElementAverageValue
    block = 0
    variable = T
    execute_on = 'initial linear'
  []
  [avg_specific_heat]
    type = ElementAverageValue
    block = 0
    variable = specific_heat
    execute_on = 'initial linear'
  []
[]

[Outputs]
  [out]
    type = Exodus
  []
[]

(test/tests/solid_mechanics/fecral_elasticity_tensor/elasticity_Fecralloy.i)

# This test case is prepared to test the elastic moduli of Fecralloy FeCrAl alloy.
#
# The geometry represents a ring of Fecralloy with inner radius of 0.005 m,
# outer radius of 0.0055 m, and a height of 0.01 m.
#
# An axial load is applied to the top of the ring with a value of 100 MPa.
# The strain in the radial direction can be calculated using the Young's modulus
# and the Poisson's ratio from the applied 100.0 MPa pressure.
#
# epsilon_{rr} = nu * sigma_{axial} / E
#
# The Young's modulus (E) for Fecralloy is taken as 180 GPa.
#
# The Poisson's ratio of Fecralloy is calculated is taken as 0.3.
#
# The results of BISON compared to the analytical solution for an increase in
# temperature under constant pressure are summarized:
#
#                          BISON             Analytical
#    Temperature (K)    strain_rr (m/m)    strain_rr (m/m)
#        420              1.666667e-4        1.666667e-4
#        540              1.666667e-4        1.666667e-4
#        660              1.666667e-4        1.666667e-4
#        780              1.666667e-4        1.666667e-4
#        900              1.666667e-4        1.666667e-4
#       1020              1.666667e-4        1.666667e-4
#       1140              1.666667e-4        1.666667e-4
#       1260              1.666667e-4        1.666667e-4
#       1380              1.666667e-4        1.666667e-4
#       1500              1.666667e-4        1.666667e-4
#

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

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

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

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = FINITE
        block = 0
        add_variables = true
        generate_output = 'elastic_strain_xx'
      []
    []
  []
[]

[Functions]
  [temperature_function]
    type = PiecewiseLinear
    x = '0       10'
    y = '300   1500'
  []
  [pressure_function]
    type = PiecewiseLinear
    x = '0  10'
    y = '1   1'
  []
[]

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

[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
  []

  [temp_bc]
    type = FunctionDirichletBC
    variable = temperature
    boundary = 'bottom top left right'
    function = temperature_function
  []
[]

[Materials]
  [elasticity_tensor]
    type = FeCrAlElasticityTensor
    temperature = temperature
    fecral_material_type = FECRALLOY
  []
  [stress]
    type = ComputeFiniteStrainElasticStress
  []
  [clad_density]
    type = StrainAdjustedDensity
    strain_free_density = 1.0
  []
  [thermal]
    type = FeCrAlThermal
    temperature = temperature
    material = FECRALLOY
  []
[]

[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   = 10
  dt         = 1
[]

[Postprocessors]
  [elastic_strain_rr]
    type = ElementAverageValue
    variable = elastic_strain_xx
  []
  [temp]
    type = AverageNodalVariableValue
    variable = temperature
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/fecral_oxidation/corrosion_test_fecral_pwr.i)

################################################################################
#
# The test case is a plate with an axial length of 0.01 meter and
# a thickness of 0.05cm; cladding OD temperature is 500 K;
# EOL time = 277.778 hours (1000000 seconds).
#
# The model calculates the oxide (oxygen) mass gain and oxide layer thickness
# based upon the latest data of Terrani et al.:
#
# [1] Terrani et al. "Uniform corrosion of FeCrAl alloys in LWR Coolant
#     environments, Journal of Nuclear Materials, p. 36-47, 2016"
#
# For PWR environments the parabolic rate constant for the oxidation correlation
# is given by: k_1 = 3.49e-3 mg/cm^2-h^(1/2) which is the average of all alloys
# studied.  The analytical oxygen mass gain is then determined by:
#
# w = k_1 * sqrt(t)                                              (1)
#
# Where t is the time within the oxidizing environment.  Then the oxide
# thickness is determined from w and the oxygen density in the oxide, d:
#
# x = w / d                                                      (2)
#
# The oxygen mass gain calculated using equation (1) and converting to SI units
# gives w = 5.81667e-4 kg/m^2.  BISON calculates an oxygen mass gain of
# 5.81667e-4 kg/m^2 as well.
#
# Then in SI units of meters the oxide thickness layer is determined to be
# 4.03935e-7 m (0.403935 um). BISON calculates an oxide thickness of 4.03935e-7 m.
#
# Mesh:
# Number of elements in radial direction = 5
# Number of elements in axial  direction = 1
#
################################################################################
[Mesh]
  [mesh]
    type = FileMeshGenerator
    file = clad_rz_short_rev4.e
  []
[]

[AuxVariables]
   [oxide_thickness]
     order  = CONSTANT
     family = MONOMIAL
   []
   [mass_gain]
     order = CONSTANT
     family = MONOMIAL
   []
[]
[Variables]
   [temp]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300
  []
[]

[AuxKernels]
  [oxide]
    type = MaterialRealAux
    variable = oxide_thickness
    property = scale_thickness
    boundary = 3
  []
  [mass_gain]
    type = MaterialRealAux
    variable = mass_gain
    property = oxide_mass_gain
    boundary = 3
  []
[]

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

 [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temp
  []
[]

[BCs]
  [clad_outer_surface]
    type =  DirichletBC
    boundary =  3
    value =  500
    variable =  temp
  []
[]

[Materials]
  [clad_thermal]
    type = FeCrAlThermal
    material = APMT
    block = 1
    temperature = temp
  []
  [density]
    type = ParsedMaterial
    block = 1
    property_name = density
    expression = 7250.0
  []
  [oxidation]
    type = FeCrAlOxidation
    reactor_type = PWR
    boundary = 3
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  l_max_its = 200
  l_tol = 8e-3

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

  start_time = 0.0
  dt = 10000
  end_time = 1000000
[]

[Outputs]
  exodus = true
  print_linear_residuals = true
[]

(test/tests/thermalFeCrAl/thermalFeCrAl_MA956.i)

# The mesh is a 1x1x1 cube (single element).
# The temperature is ramped on all faces of the cube from 300 K to 1500K
# over 100 Ms. The cladding alloy is INCOLOY MA956.
#
# The thermal conductivity and specific heat computed by BISON
# for the temperatures shown are compared with the analytical solution below.
# The results are as follows:
#
# Temp(K)    k BISON     k analytical      cp BISON      cp analytical
#            (W/m-K)       (W/m-K)         (J/kg-K)         (J/kg-K)
#
#   300     11.01131      11.01131         470.8838         470.8838
#   408     12.79245      12.79245         500.7580         500.7580
#   504     14.90275      14.90275         527.6380         527.6380
#   600     15.80275      15.80275         554.5180         554.5180
#   708     17.42275      17.42275         586.5005         586.5005
#   804     18.83190      18.83190         614.7870         614.7870
#   900     20.17590      20.17590         637.5180         637.5180
#  1008     21.68790      21.68790         667.7580         667.7580
#  1104     23.06275      23.06275         694.6380         694.6380
#  1200     24.47590      24.47590         721.2495         721.2495
#  1308     26.02275      26.02275         750.7580         750.7580
#  1404     27.00000      27.00000         769.0000         769.0000
#  1500     27.00000      27.00000         769.0000         769.0000
#
#------------------------------------------------------------------------

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

[Variables]
  [T]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300     # set initial T to 500 K
  []
[]

[AuxVariables]
  [th_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [specific_heat]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = T
  []
[]

[AuxKernels]
  [th_cond]
    type = MaterialRealAux
    variable = th_cond
    property = thermal_conductivity
    block = 0
    execute_on = 'initial linear'
  []
  [specific_heat]
    type = MaterialRealAux
    variable = specific_heat
    property = specific_heat
    block = 0
    execute_on = 'initial linear'
  []
[]

[Functions]
  [temp_ramp]
    type = PiecewiseLinear
    x = '0.0 100.0e6'
    y = '300 1500'
  []
[]

[BCs]
   [VariableT]
     type = FunctionDirichletBC
     boundary = 'top bottom front back left right'       # All cube faces
     variable = T
     function = temp_ramp
   []
[]

[Materials]
  [fuel_thermalFeCrAl]
    type = FeCrAlThermal
    block = 0
    material = MA956
    scale_factor_k = 1.0
    scale_factor_cp = 1.0
    temperature = T
  []
  [density]
    type = ParsedMaterial
    block = 0
    property_name = density
    expression = 7250.0
  []
[]

[Executioner]
  type = Transient

  solve_type = PJFNK

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

  line_search = 'none'

  l_max_its = 100
  nl_max_its = 100
  nl_rel_tol = 1e-4
  nl_abs_tol = 1e-6
  l_tol = 1e-5
  start_time = 0.0
  num_steps = 100
  dt = 1.0e6
[]

[Postprocessors]
  [ave_temp_interior]
    type = SideAverageValue
    boundary = 'top bottom front back left right'
    variable = T
    execute_on = 'initial linear'
  []
  [avg_th_cond]
    type = ElementAverageValue
    block = 0
    variable = th_cond
    execute_on = 'initial linear'
  []
  [average_clad_T]
    type = ElementAverageValue
    block = 0
    variable = T
    execute_on = 'initial linear'
  []
  [avg_specific_heat]
    type = ElementAverageValue
    block = 0
    variable = specific_heat
    execute_on = 'initial linear'
  []
[]

[Outputs]
  [out]
    type = Exodus
  []
[]

(test/tests/solid_mechanics/fecral_elasticity_tensor/elasticity_PM2000.i)

# This test case is prepared to test the elastic moduli of PM2000 FeCrAl alloy.
#
# The geometry represents a ring of PM2000 with inner radius of 0.005 m,
# outer radius of 0.0055 m, and a height of 0.01 m.
#
# An axial load is applied to the top of the ring with a value of 100 MPa.
# The strain in the radial direction can be calculated using the Young's modulus
# and the Poisson's ratio from the applied 100.0 MPa pressure.
#
# epsilon_{rr} = nu * sigma_{axial} / E
#
# The Young's modulus (E) for PM2000 is calculated via a piecewse linear function.
#
# The Poisson's ratio of PM2000 is calculated is taken as 0.33.
#
# The results of BISON compared to the analytical solution for an increase in
# temperature under constant pressure are summarized:
#
#                          BISON             Analytical
#    Temperature (K)    strain_rr (m/m)    strain_rr (m/m)
#        420              2.062500e-4        2.062500e-4
#        540              2.129221e-4        2.129221e-4
#        660              2.260487e-4        2.260487e-4
#        780              2.456911e-4        2.456911e-4
#        900              2.697952e-4        2.697952e-4
#       1020              2.929232e-4        2.929232e-4
#       1140              3.194115e-4        3.194115e-4
#       1260              3.473684e-4        3.473684e-4
#       1380              3.473684e-4        3.473684e-4
#       1500              3.473684e-4        3.473684e-4
#

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

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

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

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = FINITE
        block = 0
        add_variables = true
        generate_output = 'elastic_strain_xx'
      []
    []
  []
[]

[Functions]
  [temperature_function]
    type = PiecewiseLinear
    x = '0       10'
    y = '300   1500'
  []
  [pressure_function]
    type = PiecewiseLinear
    x = '0  10'
    y = '1   1'
  []
[]

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

[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
  []

  [temp_bc]
    type = FunctionDirichletBC
    variable = temperature
    boundary = 'bottom top left right'
    function = temperature_function
  []
[]

[Materials]
  [elasticity_tensor]
    type = FeCrAlElasticityTensor
    temperature = temperature
    fecral_material_type = PM2000
  []
  [stress]
    type = ComputeFiniteStrainElasticStress
  []
  [clad_density]
    type = StrainAdjustedDensity
    strain_free_density = 1.0
  []
  [thermal]
    type = FeCrAlThermal
    temperature = temperature
    material = PM2000
  []
[]

[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   = 10
  dt         = 1
[]

[Postprocessors]
  [elastic_strain_rr]
    type = ElementAverageValue
    variable = elastic_strain_xx
  []
  [temp]
    type = AverageNodalVariableValue
    variable = temperature
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/solid_mechanics/fecral_elasticity_tensor/elasticity_APMT.i)

# This test case is prepared to test the elastic moduli of APMT FeCrAl alloy.
#
# The geometry represents a ring of APMT with inner radius of 0.005 m,
# outer radius of 0.0055 m, and a height of 0.01 m.
#
# An axial load is applied to the top of the ring with a value of 100 MPa.
# The strain in the radial direction can be calculated using the Young's modulus
# and the Poisson's ratio from the applied 100.0 MPa pressure.
#
# epsilon_{rr} = nu * sigma_{axial} / E
#
# The Young's modulus (E) for APMT is calculated via a piecewse linear function.
#
# The Poisson's ratio of APMT is calculated is taken as 0.3.
#
# The results of BISON compared to the analytical solution for an increase in
# temperature under constant pressure are summarized:
#
#                          BISON             Analytical
#    Temperature (K)    strain_rr (m/m)    strain_rr (m/m)
#        420              1.444687e-4        1.444687e-4
#        540              1.500103e-4        1.500103e-4
#        660              1.570794e-4        1.570794e-4
#        780              1.673033e-4        1.673033e-4
#        900              1.793025e-4        1.793025e-4
#       1020              1.931558e-4        1.931558e-4
#       1140              2.093291e-4        2.093291e-4
#       1260              2.284583e-4        2.284583e-4
#       1380              2.307692e-4        2.307692e-4
#       1500              2.307692e-4        2.307692e-4
#

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

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

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

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = FINITE
        block = 0
        add_variables = true
        generate_output = 'elastic_strain_xx'
      []
    []
  []
[]

[Functions]
  [temperature_function]
    type = PiecewiseLinear
    x = '0       10'
    y = '300   1500'
  []
  [pressure_function]
    type = PiecewiseLinear
    x = '0  10'
    y = '1   1'
  []
[]

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

[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
  []

  [temp_bc]
    type = FunctionDirichletBC
    variable = temperature
    boundary = 'bottom top left right'
    function = temperature_function
  []
[]

[Materials]
  [elasticity_tensor]
    type = FeCrAlElasticityTensor
    temperature = temperature
    fecral_material_type = APMT
  []
  [stress]
    type = ComputeFiniteStrainElasticStress
  []
  [clad_density]
    type = StrainAdjustedDensity
    strain_free_density = 1.0
  []
  [thermal]
    type = FeCrAlThermal
    temperature = temperature
    material = APMT
  []
[]

[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   = 10
  dt         = 1
[]

[Postprocessors]
  [elastic_strain_rr]
    type = ElementAverageValue
    variable = elastic_strain_xx
  []
  [temp]
    type = AverageNodalVariableValue
    variable = temperature
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/thermalFeCrAl/thermalFeCrAl_C36M.i)

# The mesh is a 1x1x1 cube (single element).
# The temperature is ramped on all faces of the cube from 300 K to 1500K
# over 100 Ms. The cladding alloy is C36M.
#
# The thermal conductivity and specific heat computed by BISON
# for the temperatures shown are compared with the analytical solution below.
# The results are as follows:
#
# Temp(K)    k BISON     k analytical      cp BISON      cp analytical
#            (W/m-K)       (W/m-K)         (J/kg-K)         (J/kg-K)
#
#   300     12.20834      12.20834         484.6620         484.6620
#   408     13.61556      13.61556         537.7144         537.7144
#   504     14.84843      14.84843         575.4794         575.4794
#   600     16.06438      16.06438         629.3760         629.3760
#   708     17.41208      17.41208         739.1408         739.1408
#   804     18.59205      18.59205         741.1624         741.1624
#   900     19.75510      19.75510         693.1644         693.1644
#  1008     21.04329      21.04329         680.9782         680.9782
#  1104     22.17036      22.17036         680.5240         680.5240
#  1200     23.28050      23.28050         687.3541         687.3541
#  1308     24.50918      24.50918         704.5105         704.5105
#  1404     25.58336      25.58336         729.8179         729.8179
#  1500     26.64060      26.64060         766.5099         766.5099
#
#------------------------------------------------------------------------
[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 3
  []
[]

[Variables]
  [T]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300     # set initial T to 500 K
  []
[]

[AuxVariables]
  [th_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [specific_heat]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [th_cond]
    type = MaterialRealAux
    variable = th_cond
    property = thermal_conductivity
    execute_on = 'initial linear'
    block = 0
  []
  [specific_heat]
    type = MaterialRealAux
    variable = specific_heat
    property = specific_heat
    execute_on = 'initial linear'
    block = 0
  []
[]

[Functions]
  [temp_ramp]
    type = PiecewiseLinear
    x = '0.0 100.0e6'
    y = '300 1500'
  []
[]

[BCs]
   [VariableT]
     type = FunctionDirichletBC
     boundary = 'top bottom front back left right'      # All cube faces
     variable = T
     function = temp_ramp
   []
[]

[Materials]
  [fuel_thermalFeCrAl]
    type = FeCrAlThermal
    block = 0
    material = C36M
    scale_factor_k = 1.0
    scale_factor_cp = 1.0
    temperature = T
  []
  [density]
    type = ParsedMaterial
    block = 0
    property_name = density
    expression = 7250.0
  []
[]

[Executioner]
   type = Transient

   solve_type = PJFNK

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

   line_search = 'none'

   l_max_its = 100
   nl_max_its = 100
   nl_rel_tol = 1e-4
   nl_abs_tol = 1e-6
   l_tol = 1e-5
   start_time = 0.0
   num_steps = 100
   dt = 1.0e6
[]

[Postprocessors]
  [ave_temp_interior]
    type = SideAverageValue
    boundary = 'top bottom front back left right'
    variable = T
    execute_on = 'initial linear'
  []
  [avg_th_cond]
    type = ElementAverageValue
    block = 0
    variable = th_cond
    execute_on = 'initial linear'
  []
  [average_clad_T]
    type = ElementAverageValue
    block = 0
    variable = T
    execute_on = 'initial linear'
  []
  [avg_specific_heat]
    type = ElementAverageValue
    block = 0
    variable = specific_heat
    execute_on = 'initial linear'
  []
[]

[Outputs]
  [out]
    type = Exodus
  []
[]

(examples/accident_tolerant_fuel/uo2_fecral/uo2_fecral.i)


initial_fuel_density = 10431.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
  displacements = 'disp_x disp_y'
[]

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

[Mesh]
  coord_type = RZ
  displacements = 'disp_x disp_y'
  patch_size = 10 # For contact algorithm
  patch_update_strategy = auto
  partitioner = centroid
  centroid_partitioner_direction = y
  [mesh]
    type = FileMeshGenerator
    file = uo2_fecral_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]
  [fast_neutron_flux]
    block = clad
  []
  [fast_neutron_fluence]
    block = clad
  []
  [grain_radius]
    block = pellet_type_1
    initial_condition = 10e-6
  []
  [gap_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxide_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [total_hoop_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_hoop]
    order = CONSTANT
    family = MONOMIAL
  []
  [hoop_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_rate]
    order = CONSTANT
    family = MONOMIAL
  []
  [mass_gain]
    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.'
  []
[]

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

[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
  []
[]

[AuxKernels]
  [fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    block = clad
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    factor = 3e13
    execute_on = timestep_begin
  []
  [fast_neutron_fluence]
    type = FastNeutronFluenceAux
    variable = fast_neutron_fluence
    block = clad
    fast_neutron_flux = fast_neutron_flux
    execute_on = timestep_begin
  []
  [grain_radius]
    type = GrainRadiusAux
    block = pellet_type_1
    variable = grain_radius
    temperature = temp
    execute_on = linear
  []
  [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
  []
  [creep_strain_hoop]
    type = RankTwoScalarAux
    rank_two_tensor = creep_strain
    variable = creep_strain_hoop
    scalar_type = HoopStress
    block = clad
    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
  []
  [creep_rate]
    type = MaterialRealAux
    variable = creep_rate
    property = creep_rate
    execute_on = timestep_end
    block = clad
  []
  [oxide]
    type = MaterialRealAux
    variable = oxide_thickness
    property = scale_thickness
    boundary = 2
  []
  [mass_gain]
    type = MaterialRealAux
    variable = mass_gain
    property = oxide_mass_gain
    boundary = 2
  []
[]

[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]
    type = DirichletBC
    variable = disp_x
    boundary = 12
    value = 0.0
  []

  [no_y_clad_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 1
    value = 0.0
  []

  [no_y_fuel_bottom]
    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.5e-3 # m
    rod_pitch = 1.26e-2 # m
    linear_heat_rate = power_history
    axial_power_profile = axial_peaking_factors
    oxide_thickness = oxide_thickness
  []
[]

[Materials]
  [fuel_thermal]
    type = UO2Thermal
    block = pellet_type_1
    thermal_conductivity_model = NFIR
    temperature = temp
    burnup_function = burnup
  []
  [fuel_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = pellet_type_1
    youngs_modulus = 2.0e11
    poissons_ratio = 0.345
  []
  [elastic_stress]
    type = ComputeSmearedCrackingStress
    block = pellet_type_1
    cracking_stress = 1.68e8
    inelastic_models = 'fuel_creep'
    softening_models = exponential_softening
    shear_retention_factor = 0.1
    max_stress_correction = 0
    cracked_elasticity_type = DIAGONAL
    output_properties = crack_damage
    outputs = exodus
  []
  [exponential_softening]
    type = ExponentialSoftening
  []
  [fuel_creep]
    type = UO2CreepUpdate
    block = pellet_type_1
    burnup_function = burnup
    temperature = temp
  []
  [fuel_relocation]
    type = UO2RelocationEigenstrain
    block = pellet_type_1
    burnup_function = burnup
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    fuel_pin_geometry = 'pin_geometry'
    relocation_activation1 = 5000
    relocation_model = ESCORE_modified
    eigenstrain_name = fuel_relocation_strain
  []
  [fuel_thermal_expansion]
    type = UO2ThermalExpansionMATPROEigenstrain
    block = pellet_type_1
    temperature = temp
    stress_free_temperature = 295.0
    eigenstrain_name = fuel_thermal_strain
  []
  [fuel_volumetric_swelling]
    type = UO2VolumetricSwellingEigenstrain
    gas_swelling_model_type = SIFGRS
    block = pellet_type_1
    temperature = temp
    burnup_function = burnup
    initial_fuel_density = 10431.0
    eigenstrain_name = fuel_volumetric_strain
  []
  [clad_thermal]
    type = FeCrAlThermal
    material = C35M
    block = clad
    temperature = temp
  []
  [clad_elasticity_tensor] # isotropic elasticity tensor for Zry cladding
    type = FeCrAlElasticityTensor
    temperature = temp
    fecral_material_type = C35M
    block = clad
  []
  [clad_stress] # stress update class to govern the return mapping algorithm for creep
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_creep clad_plasticity'
    block = clad
  []
  [clad_creep]
    type = FeCrAlCreepUpdate
    block = clad
    temperature = temp
    fecral_material_type = C35M
    fast_neutron_flux = fast_neutron_flux
    model_irradiation_creep = true
    model_thermal_creep = true
    max_inelastic_increment = 1e-4
  []
  [thermal_expansion]
    type = FeCrAlThermalExpansionEigenstrain
    block = clad
    temperature = temp
    fecral_material_type = C35M
    stress_free_temperature = 295.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [irradiation_swelling]
    type = FeCrAlVolumetricSwellingEigenstrain
    block = clad
    temperature = temp
    fast_neutron_fluence = fast_neutron_fluence
    eigenstrain_name = clad_volumetric_strain
  []
  [clad_plasticity]
    type = FeCrAlPlasticityUpdate
    block = clad
    hardening_constant = 2.5e9
    temperature = temp
    yield_stress = 500.0
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = pellet_type_1
    temperature = temp
    burnup_function = burnup
    grain_radius = grain_radius
    gbs_model = true
  []
  [clad_density]
    type = StrainAdjustedDensity
    block = clad
    strain_free_density = 7250.0
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = pellet_type_1
    strain_free_density = ${initial_fuel_density}
  []

  [failure_criterion]
    type = FeCrAlCladdingFailure
    boundary = '2 5'
    hoop_stress = hoop_stress
    failure_criterion = UTS
    temperature = temp
  []
  [oxidation]
    type = FeCrAlOxidation
    reactor_type = PWR
    boundary = 2
  []
[]

[Dampers]
  [limitT]
    type = BoundingValueNodalDamper
    max_value = 3200.0
    min_value = 293.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 = 1e5
  dtmin = 1

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 2.0e2
    force_step_every_function_point = true
    timestep_limiting_function = power_history
    max_function_change = 5e5
    optimal_iterations = 10
    iteration_window = 2
    linear_iteration_ratio = 100
    growth_factor = 2.0
    timestep_limiting_postprocessor = material_timestep
  []

  [Quadrature]
    order = FIFTH
    side_order = SEVENTH
  []
[]

[Postprocessors]
  [ave_temp_interior]
    type = SideAverageValue
    boundary = 9
    variable = temp
    execute_on = 'initial linear'
  []

  [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'
  []

  [_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
  []

  [alhr_input]
    type = FunctionValuePostprocessor
    function = power_history
  []
  [average_burnup]
    type = ElementAverageValue
    block = pellet_type_1
    variable = burnup
  []
  [oxide_thickness]
    type = ElementExtremeValue
    block = clad
    variable = oxide_thickness
  []
  [mass_gain]
    type = ElementExtremeValue
    block = clad
    variable = mass_gain
  []
  [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
  csv = true
  print_linear_residuals = true
  color = false

  [console]
    type = Console
    max_rows = 25
  []
[]

(test/tests/solid_mechanics/fecral_elasticity_tensor/elasticity_C06M.i)

# This test case is prepared to test the elastic moduli of C06M FeCrAl alloy.
#
# The geometry represents a ring of C06M with inner radius of 0.005 m,
# outer radius of 0.0055 m, and a height of 0.01 m.
#
# An axial load is applied to the top of the ring with a value of 100 MPa.
# The strain in the radial direction can be calculated using the Young's modulus
# and the Poisson's ratio from the applied 100.0 MPa pressure.
#
# epsilon_{rr} = nu * sigma_{axial} / E
#
# The Young's modulus (E) for C06M is calculated via:
#
# E = -5.46e-5 * (T - 273.15)^2 - 3.85e-2 * (T - 273.15) + 1.99e2) * 1e9
#
# where E is in Pa and T is the temperature in K.
#
# The Poisson's ratio of C06M is calculated via:
#
# nu = 4.46e-5 * (t - 273.15) + 0.27
#
# The results of BISON compared to the analytical solution for an increase in
# temperature under constant pressure are summarized:
#
#                          BISON             Analytical
#    Temperature (K)    strain_rr (m/m)    strain_rr (m/m)
#        420              1.439097e-4        1.439097e-4
#        540              1.525125e-4        1.525125e-4
#        660              1.632723e-4        1.632723e-4
#        780              1.768440e-4        1.768440e-4
#        900              1.942209e-4        1.942209e-4
#       1020              2.169732e-4        2.169732e-4
#       1140              2.477254e-4        2.477254e-4
#       1260              2.912040e-4        2.912040e-4
#       1380              3.568533e-4        3.568533e-4
#       1500              4.666521e-4        4.666521e-4
#
# The answers are the same for the FeCrAlloy alloys C35M and C36M.

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

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

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

[Physics]
  [SolidMechanics]
    [QuasiStatic]
      [all]
        strain = FINITE
        block = 0
        add_variables = true
        generate_output = 'elastic_strain_xx'
      []
    []
  []
[]

[Functions]
  [temperature_function]
    type = PiecewiseLinear
    x = '0       10'
    y = '300   1500'
  []
  [pressure_function]
    type = PiecewiseLinear
    x = '0  10'
    y = '1   1'
  []
[]

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

[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
  []

  [temp_bc]
    type = FunctionDirichletBC
    variable = temperature
    boundary = 'bottom top left right'
    function = temperature_function
  []
[]

[Materials]
  [elasticity_tensor]
    type = FeCrAlElasticityTensor
    temperature = temperature
    fecral_material_type = C06M
  []
  [stress]
    type = ComputeFiniteStrainElasticStress
  []
  [clad_density]
    type = StrainAdjustedDensity
    strain_free_density = 1.0
  []
  [thermal]
    type = FeCrAlThermal
    temperature = temperature
    material = C06M
  []
[]

[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   = 10
  dt         = 1
[]

[Postprocessors]
  [elastic_strain_rr]
    type = ElementAverageValue
    variable = elastic_strain_xx
  []
  [temp]
    type = AverageNodalVariableValue
    variable = temperature
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/solid_mechanics/fecral_eigenstrains/fecral_vswelling/swelling.i)

# Test of volumetric swelling of FeCrAl alloys (1000K).
#
# The model used is based upon the upper limit of swelling suggested by:
#
# K. A. Terrani, T. M. Karlsen, and Y. Yamamoto, "Input Correlations for
# Irradiation Creep of FeCrAl and SiC based on in-pile Halden test results,"
# Technical Report ORNL/TM-2016/191, Oak Ridge National Laboratory, May 2016
#
# The suggested upper limit is 0.05% per dpa.  Using a conversion factor of
# 1x10^25 neutrons/m^2-s = 0.9 dpa as proposed by:
#
# K. G. Field, X. Hu, K. C. Littrell, Y. Yamamoto, and L. L. Snead, "Radiation
# tolerance of neutroni-irradiated model Fe-Cr-Al alloys," Journal of Nuclear
# Materials, Vol. 465, pp. 746-755, 2015.
#
# A creep strain rate due to swelling is obtained and given by:
#
# e_dot_sw = 0.0005 per dpa
#
# e_dot_sw = 0.0005 * 0.9 / 1e25 [m^2-s/n]
#
# To get the correct units multiply the factor above by flux
#
# e_dot_sw = 4.5e-29 * fast_flux
#
# To get the volumteric strain the above equation must be integrated over time:
#
# e_sw = 4.5e-29 * fast_fluence
#
# With a fast fluence of 1e24, the volumetric strain should be 4.5e-5.
# This is easily checked with the 'volume' postprocessor value.

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

[Mesh]
  [mesh]
    type = FileMeshGenerator
    file = patch.e
  []
[]

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

[AuxVariables]
  [flux]
    order = FIRST
    family = LAGRANGE
  []
  [fluence]
    order = FIRST
    family = LAGRANGE
  []
  [swelling]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [all]
    strain = finite
    eigenstrain_names = 'fuelthermal_strain swell'
  []
[]

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

[AuxKernels]
  [flux_neutron_flux]
    type = FastNeutronFluxAux
    variable = flux
    factor = 1e18
    execute_on = 'initial timestep_begin'
  []
  [fast_neutron_fluence]
    type = FastNeutronFluenceAux
    variable = fluence
    block = '1 2 3 4 5 6 7'
    fast_neutron_flux = flux
    execute_on = timestep_begin
  []
  [swelling]
    type = MaterialRealAux
    variable = swelling
    property = volumetric_swelling_strain
    execute_on = timestep_end
  []
[]

[BCs]
  [bottom_x]
    type = DirichletBC
    variable = disp_x
    boundary = 10
    value = 0.0
  []
  [bottom_y]
    type = DirichletBC
    variable = disp_y
    boundary = 9
    value = 0.0
  []
  [bottom_z]
    type = DirichletBC
    variable = disp_z
    boundary = 14
    value = 0.0
  []
  [bottom_temp]
    type = DirichletBC
    variable = temp
    boundary = 3
    value = 1000.0
  []
  [top_temp]
    type = DirichletBC
    variable = temp
    boundary = 5
    value = 1000.0
  []
[]

[Materials]
  [clad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = '1 2 3 4 5 6 7'
    youngs_modulus = 1.0
    poissons_ratio = 0.3
  []
  [clad_swelling]
    type = FeCrAlVolumetricSwellingEigenstrain
    block = '1 2 3 4 5 6 7'
    temperature = temp
    fast_neutron_fluence = fluence
    eigenstrain_name = swell
  []
  [clad_thermal_strain]
    type = ComputeThermalExpansionEigenstrain
    block = '1 2 3 4 5 6 7'
    thermal_expansion_coeff = 1.0e-5
    temperature = temp
    stress_free_temperature = 1000.0
    eigenstrain_name = fuelthermal_strain
  []
  [clad_stress]
    type = ComputeFiniteStrainElasticStress
    block = '1 2 3 4 5 6 7'
  []
  [thermal]
    type = FeCrAlThermal
    block = '1 2 3 4 5 6 7'
    material = APMT
    temperature = temp
  []
  [density]
    type = StrainAdjustedDensity
    block = '1 2 3 4 5 6 7'
    strain_free_density = 7225.0
  []
[]

[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'
  l_max_its = 60
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-10
  l_tol = 1e-6

  start_time = 0.0
  dt = 10000
  num_steps = 100
[]

[Postprocessors]
  [volume]
    type = InternalVolume
    boundary = 15
    execute_on = 'initial linear'
  []
  [ave_fluence]
    type = ElementAverageValue
    variable = fluence
    block = '1 2 3 4 5 6 7'
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/fecral_oxidation/corrosion_test_fecral_high_temp.i)

################################################################################
#
# The test case is a plate with an axial length of 0.01 meter and
# a thickness of 0.05cm; cladding OD temperature is 1473 K;
# EOL time = 2.778 hours (10000 seconds).
#
# The model calculates the oxide (oxygen) mass gain and oxide layer thickness
# based upon the high-temperature data from:
#
# [2] K.G. Field, M.A. Snead, Y. Yamamoto, K.A. Terrani, Handbook on the Material
#     Properties of FeCrAl Alloys for Nuclear Power Production Applications,
#     ORNL/SPR-2018/905 (2018).
#
# The model considers the growth of a scale of alumina (Al2O3) above 1173 K.
# The parabolic rate constant for the oxidation correlation is given by:
# k_1 = sqrt(1.803e3*exp(-3.440e5/(RT))) kg^2/m^4-s, where the pre-factor
# 1.803e3 corresponds to 2.3 times that of the APMT alloy [2].
# The analytical oxygen mass gain is then determined by:
#
# w = k_1 * sqrt(t)                                              (1)
#
# where t is the time within the oxidizing environment.  Then the oxide
# thickness is determined from w and the oxygen density in the oxide, d:
#
# x = w / d                                                      (2)
#
# The oxygen mass gain calculated using equation (1) gives w = 3.38e-3 kg/m^2.
# BISON calculates an oxygen mass gain of 3.38e-3 kg/m^2 as well.
#
# Then the oxide thickness layer is determined by equation (2) to be 1.79e-6 m
# (1.79 um). BISON calculates an oxide thickness of 1.79e-6 m as well.
#
# Mesh:
# Number of elements in radial direction = 5
# Number of elements in axial  direction = 1
#
################################################################################
[Mesh]
  [mesh]
    type = FileMeshGenerator
    file = clad_rz_short_rev4.e
  []
[]

[AuxVariables]
   [oxide_thickness]
     order  = CONSTANT
     family = MONOMIAL
   []
   [mass_gain]
     order = CONSTANT
     family = MONOMIAL
   []
[]
[Variables]
   [temp]
    order = FIRST
    family = LAGRANGE
    initial_condition = 1473
  []
[]

[AuxKernels]
  [oxide]
    type = MaterialRealAux
    variable = oxide_thickness
    property = scale_thickness
    boundary = 3
  []
  [mass_gain]
    type = MaterialRealAux
    variable = mass_gain
    property = oxide_mass_gain
    boundary = 3
  []
[]

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

 [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temp
  []
[]

[BCs]
  [clad_outer_surface]
    type =  DirichletBC
    boundary =  3
    value =  1473
    variable =  temp
  []
[]

[Materials]
  [clad_thermal]
    type = FeCrAlThermal
    material = APMT
    block = 1
    temperature = temp
  []
  [density]
    type = ParsedMaterial
    block = 1
    property_name = density
    expression = 7250.0
  []
  [oxidation]
    type = FeCrAlOxidation
    reactor_type = PWR
    temperature = temp
    boundary = 3
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  l_max_its = 200
  l_tol = 8e-3

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

  start_time = 0.0
  dt = 20
  end_time = 10000
[]

[Outputs]
  exodus = true
  print_linear_residuals = true
[]

(test/tests/thermalFeCrAl/thermalFeCrAl_C06M.i)

# The mesh is a 1x1x1 cube (single element).
# The temperature is ramped on all faces of the cube from 300 K to 1500K
# over 100 Ms. The cladding alloy is C06M.
#
# The thermal conductivity and specific heat computed by BISON
# for the temperatures shown are compared with the analytical solution below.
# The results are as follows:
#
# Temp(K)    k BISON     k analytical      cp BISON      cp analytical
#            (W/m-K)       (W/m-K)         (J/kg-K)         (J/kg-K)
#
#   300     13.11286      13.11286         444.5820         444.5820
#   408     14.27912      14.27912         513.1303         513.1303
#   504     15.32905      15.32905         559.6106         559.6106
#   600     16.39143      16.39143         607.1760         607.1760
#   708     17.60151      17.60151         679.5392         679.5392
#   804     18.69039      18.69039         776.2241         776.2241
#   900     19.79172      19.79172         852.0768         852.0768
#  1008     21.04562      21.04562         733.6445         733.6445
#  1104     22.17344      22.17344         708.5530         708.5530
#  1200     23.31373      23.31373         699.4732         699.4732
#  1308     24.61145      24.61145         704.4914         704.4914
#  1404     25.77822      25.77822         723.6579         723.6579
#  1500     26.95745      26.95745         758.9976         758.9976
#
#------------------------------------------------------------------------
[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 3
  []
[]

[Variables]
  [T]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300     # set initial T to 500 K
  []
[]

[AuxVariables]
  [th_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [specific_heat]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [th_cond]
    type = MaterialRealAux
    variable = th_cond
    property = thermal_conductivity
    execute_on = 'initial linear'
    block = 0
  []
  [specific_heat]
    type = MaterialRealAux
    variable = specific_heat
    property = specific_heat
    execute_on = 'initial linear'
    block = 0
  []
[]

[Functions]
  [temp_ramp]
    type = PiecewiseLinear
    x = '0.0 100.0e6'
    y = '300 1500'
  []
[]

[BCs]
   [VariableT]
     type = FunctionDirichletBC
     boundary = 'top bottom front back left right'      # All cube faces
     variable = T
     function = temp_ramp
   []
[]

[Materials]
  [fuel_thermalFeCrAl]
    type = FeCrAlThermal
    block = 0
    material = C06M
    scale_factor_k = 1.0
    scale_factor_cp = 1.0
    temperature = T
  []
  [density]
    type = ParsedMaterial
    block = 0
    property_name = density
    expression = 7250.0
  []
[]

[Executioner]
   type = Transient

   solve_type = PJFNK

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

   line_search = 'none'

   l_max_its = 100
   nl_max_its = 100
   nl_rel_tol = 1e-4
   nl_abs_tol = 1e-6
   l_tol = 1e-5
   start_time = 0.0
   num_steps = 100
   dt = 1.0e6
[]

[Postprocessors]
  [ave_temp_interior]
    type = SideAverageValue
    boundary = 'top bottom front back left right'
    variable = T
    execute_on = 'initial linear'
  []
  [avg_th_cond]
    type = ElementAverageValue
    block = 0
    variable = th_cond
    execute_on = 'initial linear'
  []
  [average_clad_T]
    type = ElementAverageValue
    block = 0
    variable = T
    execute_on = 'initial linear'
  []
  [avg_specific_heat]
    type = ElementAverageValue
    block = 0
    variable = specific_heat
    execute_on = 'initial linear'
  []
[]

[Outputs]
  [out]
    type = Exodus
  []
[]

(test/tests/thermalFeCrAl/thermalFeCrAl_PM2000.i)

# The mesh is a 1x1x1 cube (single element).
# The temperature is ramped on all faces of the cube from 300 K to 1500K
# over 100 Ms. The cladding alloy is Plansee PM2000.
#
# The thermal conductivity and specific heat computed by BISON
# for the temperatures shown are compared with the analytical solution below.
# The results are as follows:
#
# Temp(K)    k BISON     k analytical      cp BISON      cp analytical
#            (W/m-K)       (W/m-K)         (J/kg-K)         (J/kg-K)
#
#   300     11.09408      11.09408         500.0000         500.0000
#   408     14.15408      14.15408         493.0230         493.0230
#   504     16.51416      16.51416         493.3683         493.3683
#   600     18.11412      18.11412         534.9683         534.9683
#   708     19.91417      19.91417         581.7683         581.7683
#   804     21.10284      21.10284         618.6378         618.6378
#   900     21.46283      21.46283         645.5180         645.5180
#  1008     21.78283      21.78283         675.7580         675.7580
#  1104     22.53988      22.53988         699.4040         699.4040
#  1200     24.21988      24.21988         722.4440         722.4440
#  1308     25.93563      25.93563         740.0000         740.0000
#  1404     27.13563      27.13563         740.0000         740.0000
#  1500     28.00000      28.00000         740.0000         740.0000
#
#------------------------------------------------------------------------

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

[Variables]
  [T]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300     # set initial T to 500 K
  []
[]

[AuxVariables]
  [th_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [specific_heat]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [th_cond]
    type = MaterialRealAux
    variable = th_cond
    property = thermal_conductivity
    block = 0
    execute_on = 'initial linear'
  []
  [specific_heat]
    type = MaterialRealAux
    variable = specific_heat
    property = specific_heat
    block = 0
    execute_on = 'initial linear'
  []
[]

[Functions]
  [temp_ramp]
    type = PiecewiseLinear
    x = '0.0 100.0e6'
    y = '300 1500'
  []
[]

[BCs]
  [VariableT]
    type = FunctionDirichletBC
    boundary = 'top bottom front back left right'       # All cube faces
    variable = T
    function = temp_ramp
  []
[]

[Materials]
  [fuel_thermalFeCrAl]
    type = FeCrAlThermal
    block = 0
    material = PM2000
    scale_factor_k = 1.0
    scale_factor_cp = 1.0
    temperature = T
  []
  [density]
    type = ParsedMaterial
    block = 0
    property_name = density
    expression = 7250.0
  []
[]

[Executioner]
   type = Transient

   solve_type = PJFNK

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

   line_search = 'none'

   l_max_its = 100
   nl_max_its = 100
   nl_rel_tol = 1e-4
   nl_abs_tol = 1e-6
   l_tol = 1e-5
   start_time = 0.0
   num_steps = 100
   dt = 1.0e6
[]

[Postprocessors]
  [ave_temp_interior]
    type = SideAverageValue
    boundary = 'top bottom front back left right'
    variable = T
    execute_on = 'initial linear'
  []
  [avg_th_cond]
    type = ElementAverageValue
    block = 0
    variable = th_cond
    execute_on = 'initial linear'
  []
  [average_clad_T]
    type = ElementAverageValue
    block = 0
    variable = T
    execute_on = 'initial linear'
  []
  [avg_specific_heat]
    type = ElementAverageValue
    block = 0
    variable = specific_heat
    execute_on = 'initial linear'
  []
[]

[Outputs]
  [out]
    type = Exodus
  []
[]

(test/tests/thermalFeCrAl/thermalFeCrAl_C35M.i)

# The mesh is a 1x1x1 cube (single element).
# The temperature is ramped on all faces of the cube from 300 K to 1500K
# over 100 Ms. The cladding alloy is C35M.
#
# The thermal conductivity and specific heat computed by BISON
# for the temperatures shown are compared with the analytical solution below.
# The results are as follows:
#
# Temp(K)    k BISON     k analytical      cp BISON      cp analytical
#            (W/m-K)       (W/m-K)         (J/kg-K)         (J/kg-K)
#
#   300     12.93426      12.93426         448.2600         448.2600
#   408     14.44236      14.44236         518.1462         518.1462
#   504     15.74400      15.74400         566.4534         566.4534
#   600     17.00904      17.00904         616.8000         616.8000
#   708     18.38845      18.38845         693.8556         693.8556
#   804     19.57570      19.57570         796.4774         796.4774
#   900     20.72634      20.72634         818.3847         818.3847
#  1008     21.97706      21.97706         739.9534         739.9534
#  1104     23.04991      23.04991         716.8472         716.8472
#  1200     24.08616      24.08616         709.3671         709.3671
#  1308     25.20818      25.20818         717.9575         717.9575
#  1404     26.16665      26.16665         742.7319         742.7319
#  1500     27.08850      27.08850         786.5558         786.5558
#
#------------------------------------------------------------------------
[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 3
  []
[]

[Variables]
  [T]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300     # set initial T to 500 K
  []
[]

[AuxVariables]
  [th_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [specific_heat]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [th_cond]
    type = MaterialRealAux
    variable = th_cond
    property = thermal_conductivity
    execute_on = 'initial linear'
    block = 0
  []
  [specific_heat]
    type = MaterialRealAux
    variable = specific_heat
    property = specific_heat
    execute_on = 'initial linear'
    block = 0
  []
[]

[Functions]
  [temp_ramp]
    type = PiecewiseLinear
    x = '0.0 100.0e6'
    y = '300 1500'
  []
[]

[BCs]
   [VariableT]
     type = FunctionDirichletBC
     boundary = 'top bottom front back left right'      # All cube faces
     variable = T
     function = temp_ramp
   []
[]

[Materials]
  [fuel_thermalFeCrAl]
    type = FeCrAlThermal
    block = 0
    material = C35M
    scale_factor_k = 1.0
    scale_factor_cp = 1.0
    temperature = T
  []
  [density]
    type = ParsedMaterial
    block = 0
    property_name = density
    expression = 7250.0
  []
[]

[Executioner]
   type = Transient

   solve_type = PJFNK

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

   line_search = 'none'

   l_max_its = 100
   nl_max_its = 100
   nl_rel_tol = 1e-4
   nl_abs_tol = 1e-6
   l_tol = 1e-5
   start_time = 0.0
   num_steps = 100
   dt = 1.0e6
[]

[Postprocessors]
  [ave_temp_interior]
    type = SideAverageValue
    boundary = 'top bottom front back left right'
    variable = T
    execute_on = 'initial linear'
  []
  [avg_th_cond]
    type = ElementAverageValue
    block = 0
    variable = th_cond
    execute_on = 'initial linear'
  []
  [average_clad_T]
    type = ElementAverageValue
    block = 0
    variable = T
    execute_on = 'initial linear'
  []
  [avg_specific_heat]
    type = ElementAverageValue
    block = 0
    variable = specific_heat
    execute_on = 'initial linear'
  []
[]

[Outputs]
  [out]
    type = Exodus
  []
[]

(test/tests/fecral_oxidation/corrosion_test_fecral_bwr.i)

################################################################################
#
# The test case is a plate with an axial length of 0.01 meter and
# a thickness of 0.05cm; cladding OD temperature is 500 K;
# EOL time = 277.778 hours (1000000 seconds).
#
# The model calculates the oxide (oxygen) mass gain and oxide layer thickness
# based upon the latest data of Terrani et al.:
#
# [1] Terrani et al. "Uniform corrosion of FeCrAl alloys in LWR Coolant
#     environments, Journal of Nuclear Materials, p. 36-47, 2016"
#
# For normal water chemistry BWR environments the parabolic rate constant for
# the oxidation correlation is given by: k_1 = 3.45e-4 mg/cm^2-h^(1/2) for the
# alloy with a composition very close to C35M. (Fe-13Cr-4Al).  The analytical
# oxygen mass gain is then determined by:
#
# w = k_1 * sqrt(t)                                              (1)
#
# Where t is the time within the oxidizing environment.  Then the oxide
# thickness is determined from w and the oxygen density in the oxide, d:
#
# x = w / d                                                      (2)
#
# The oxygen mass gain calculated using equation (1) and converting to SI units
# gives w = 5.75e-5 kg/m^2.  BISON calculates an oxygen mass gain of
# 5.75e-5 kg/m^2 as well.
#
# Then in SI units of meters the oxide thickness layer is determined to be
# 3.99306e-8 m (0.0399306 um). BISON calculates an oxide thickness of
# 3.99306e-8 m.
#
# Mesh:
# Number of elements in radial direction = 5
# Number of elements in axial  direction = 1
#
################################################################################
[Mesh]
  [mesh]
    type = FileMeshGenerator
    file = clad_rz_short_rev4.e
  []
[]

[AuxVariables]
   [oxide_thickness]
     order  = CONSTANT
     family = MONOMIAL
   []
   [mass_gain]
     order = CONSTANT
     family = MONOMIAL
   []
[]
[Variables]
   [temp]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300
  []
[]

[AuxKernels]
  [oxide]
    type = MaterialRealAux
    variable = oxide_thickness
    property = scale_thickness
    boundary = 3
  []
  [mass_gain]
    type = MaterialRealAux
    variable = mass_gain
    property = oxide_mass_gain
    boundary = 3
  []
[]

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

 [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temp
  []
[]

[BCs]
  [clad_outer_surface]
    type =  DirichletBC
    boundary =  3
    value =  500
    variable =  temp
  []
[]

[Materials]
  [clad_thermal]
    type = FeCrAlThermal
    material = APMT
    block = 1
    temperature = temp
  []
  [density]
    type = ParsedMaterial
    block = 1
    property_name = density
    expression = 7250.0
  []
  [oxidation]
    type = FeCrAlOxidation
    reactor_type = BWR
    boundary = 3
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  l_max_its = 200
  l_tol = 8e-3

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

  start_time = 0.0
  dt = 10000
  end_time = 1000000
[]

[Outputs]
  exodus = true
  print_linear_residuals = true
[]

(test/tests/fecral_oxidation/corrosion_test_fecral.i)

################################################################################
#
# The test case is a plate with an axial length of 0.01 meter and
# a thickness of 0.05cm; cladding OD temperature is 1473 K;
# EOL time = 1000 hours.
#
# The oxide thickness layer is analytically calculated to be 22.6162e-5
# (22.6162 um).BISON calculates an oxide thickness of 22.6162e-5 m as well.
#
# Mesh:
# Number of elements in radial direction = 5
# Number of elements in axial  direction = 1
#
#--------------------------------------------------------------------------------
[Mesh]
  [mesh]
    type = FileMeshGenerator
    file = clad_rz_short_rev4.e
  []
[]

# Define dependent variables, element order and shape function family,
# and initial conditions
[AuxVariables]
  [oxide_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
[]
[Variables]
  [temp]
    order = FIRST
    family = LAGRANGE
    initial_condition = 300 # K,  initial temperature
  []
[]

[AuxKernels]
  [oxide]
    type = FeCrAlOxideAux
    variable = oxide_thickness
    temperature = temp
    boundary = 3
    execute_on = timestep_end
  []
[]
[Kernels]
  [heat]
    type = HeatConduction
    variable = temp
  []

  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temp
  []
[]

# Define boundary conditions
[BCs]
  [clad_outer_surface]
    type = DirichletBC
    boundary = 3
    value = 1473 # K
    variable = temp
  []
[]

[Materials]
  [clad_thermal]
    type = FeCrAlThermal
    material = APMT
    block = 1
    temperature = temp
  []
  [density]
    type = ParsedMaterial
    block = 1
    property_name = density
    expression = 7250.0 # kg/m^3
  []
[]

[Executioner]
  type = Transient
  # PETSC options

  solve_type = 'PJFNK'

  # controls for linear iterations
  l_max_its = 200
  l_tol = 8e-3

  # controls for nonlinear iterations
  nl_max_its = 15
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-10

  # time control
  start_time = 0.0
  dt = 36000
  end_time = 3600000
[]

# Define output file(s)
[Outputs]
  exodus = true
  print_linear_residuals = true
[]

Description
Range of Applicability and Uncertainty
Example Input Syntax
Input Parameters
Input Files
References