CompositeSiCThermal

Computes thermal conductivity and specific heat of composite (CVI) SiC/SiC cladding.

Description

The CompositeSiCThermal model computes the thermal conductivity and specific heat capacity of composite silicon carbide.

Thermal Conductivity Models

Three models exist for the thermal conductivity of composite silicon carbide; Koyanagi, Stone, and the Combined.

Koyanagi Model

Koyanagi et al. (2017) developed a model for unirradiated composite SiC. The through thickness (of tube specimens) thermal conductivity correlation is calculated from the measured thermal diffusivity data through: $var element = document.getElementById("moose-equation-84d4477c-9af8-4843-a43e-a8825f60585e");katex.render("k = \\alpha c_P \\rho", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-46a4a91e-3d55-4c52-a3bc-be019caac70c");katex.render("k", element, {displayMode:false,throwOnError:false});$ is the thermal conductivity (W/m-K), $var element = document.getElementById("moose-equation-03e57c60-83b3-4c65-ac77-916f986bfd6b");katex.render("\\alpha", element, {displayMode:false,throwOnError:false});$ is the thermal diffusivity, $var element = document.getElementById("moose-equation-2ac86f19-71b6-4aed-88a3-c65a324b002e");katex.render("c_P", element, {displayMode:false,throwOnError:false});$ is the specific heat at constant pressure, and $var element = document.getElementById("moose-equation-ad5ec17e-f47f-4267-abcf-16877140708d");katex.render("\\rho", element, {displayMode:false,throwOnError:false});$ is the density of composite SiC. The density of composite SiC is taken as the average of the reported value by Koyanagi et al. (2017) ( $var element = document.getElementById("moose-equation-72173554-531f-4ab2-a9c0-50ef7c4d5a6b");katex.render("\\rho", element, {displayMode:false,throwOnError:false});$ = 2700 kg/m $var element = document.getElementById("moose-equation-e33182cf-2522-4ed6-b5b5-3aaf52e4a0e6");katex.render("^3", element, {displayMode:false,throwOnError:false});$ ). The specific heat is calculated by Eq. (1) or Eq. (2). The range of thermal diffusivity measurements for full SiC/SiC composite as a function of temperature are shown in Table 1.

Table 1: Ranges of temperature dependent thermal diffusivity of Full SiC/SiC composite (Koyanagi et al., 2017).

Temperature (K)	Thermal Diffusivity Range (mm $var element = document.getElementById("moose-equation-8f354ce4-adea-423f-ad2a-ca0a7f05b2e8");katex.render("^2", element, {displayMode:false,throwOnError:false});$ /s)
298.0	2-7
573.15	1.25-4.5
1073.15	1-3.25

For the computation of the thermal conductivity, the thermal diffusivity is calculated from a piecewise linear function generated from the mean of the ranges given in Table 1. After converting to SI units for the thermal diffusivity (m $var element = document.getElementById("moose-equation-b530e63c-2ec9-41d5-904e-b4e1ab398dda");katex.render("^2", element, {displayMode:false,throwOnError:false});$ /s) the tabulated values used in the piecewise linear function are shown in Table 2.

Table 2: Mean temperature dependent thermal diffusivity used in the BISON piecewise linear function.

Temperature (K)	Thermal Diffusivity (m $var element = document.getElementById("moose-equation-52a2327a-b141-4dcf-938d-34014dd4efaa");katex.render("^2", element, {displayMode:false,throwOnError:false});$ /s)
298.0	4.5 $var element = document.getElementById("moose-equation-fd1e95fd-8427-4c4c-9843-ae0a912571cf");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-09850370-a2f8-4603-83d1-243504b3c734");katex.render("^{-6}", element, {displayMode:false,throwOnError:false});$
573.15	2.875 $var element = document.getElementById("moose-equation-b4276568-1a6b-4285-a467-c65d3451766a");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-db20795d-e375-4ac4-8717-e5995009c368");katex.render("^{-6}", element, {displayMode:false,throwOnError:false});$
1073.15	2.125 $var element = document.getElementById("moose-equation-7ef6ef19-1b6a-4d3f-9026-c0a6fd4eadd5");katex.render("\\times", element, {displayMode:false,throwOnError:false});$ 10 $var element = document.getElementById("moose-equation-6a1240f1-a442-4a7c-b1dc-ac73f5f0d619");katex.render("^{-6}", element, {displayMode:false,throwOnError:false});$

Stone Model

Stone et al. (2015) developed a model that applies to both irradiated and unirradiated composite SiC. The unirradiated thermal conductivity (W/m-K) is given by: $var element = document.getElementById("moose-equation-2b060315-2c1a-4102-835b-b69903fabc41");katex.render(" k_{unirradiated} = -1.71 \\times 10^{-11} T^4 + 7.35 \\times 10^{-8} T^3 - 1.10 \\times 10^{-4} T^2 + 0.061 T + 7.97", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-2f11b569-e6e9-4bb3-a855-c53b50f7d91d");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the temperature in K. Irradiation effects are taken into account by assuming a resistance network where the total thermal conductivity is given by: $var element = document.getElementById("moose-equation-935fcd92-50c2-4d65-840b-d017452444d4");katex.render(" \\frac{1}{k_{irradiated}} = \\frac{1}{k_{unirradiated}} + R_{irradiation}", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-ac32b653-3a3b-414b-8aa0-bb6ae0f0f3f7");katex.render("R_{irradiation}", element, {displayMode:false,throwOnError:false});$ is the resistance due to irradiation effects given by: $var element = document.getElementById("moose-equation-c319666d-516a-4b1c-a220-1a2be58e2895");katex.render(" R_{irradiation} = 15.11 S", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-bee76787-de84-4249-b905-f160e312a64a");katex.render("S", element, {displayMode:false,throwOnError:false});$ is the total volumetric swelling strain calculated by CompositeSiCVolumetricSwellingEigenstrain.

Combined Model

Since Koyanagi's model is preferred for including tube specimen data, but lacks an irradiation damage description, the combined model incorporates Stone's irradiation damage resistivity term to Koyanagi's model: $var element = document.getElementById("moose-equation-f8911b32-99e2-4a4c-a0fc-7143fa641f61");katex.render(" \\frac{1}{k_{irradiated}} = \\frac{1}{k_{koyanagi.unirradiated}} + R_{stone.irradiation}", element, {displayMode:true,throwOnError:false});$

Specific Heat Capacity Models

Two models exist for the specific heat of composite silicon carbide; Snead and the GA.

Snead Model

Snead et al. (2007) used the same correlation for specific heat (J/kg-K) as for monolithic SiC (Snead et al., 2007):

$(1)var element = document.getElementById("moose-equation-9dc65793-46e1-4768-bd6c-b759e7f4b2a8");katex.render(" c_P = 925.65 + 0.3772T - 7.9259\\times10^{-5}T^2 - \\frac{3.1946\\times10^7}{T^2}", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-fe0b9e53-72c2-45f0-81c4-a901cbf1bf4c");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the temperature in K.

GA Model

GA (2020) developed a correlation to mechanically represent the measured specific heat (J/kg-K) of General Atomics (GA) prototype cladding:

$(2)var element = document.getElementById("moose-equation-862a0468-cbe5-4eec-861c-840de65d0255");katex.render(" c_P = 4.66\\times10^{-13}T^5 - 2.96\\times10^{-9}T^4 + 7.47\\times10^{-6}T^3 - 9.41\\times10^{-3}T^2 + 6.17T - 5.38\\times10^{2}", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-86cd728c-23dc-4441-91aa-6f9f74c475b9");katex.render("T", element, {displayMode:false,throwOnError:false});$ is the temperature in K.

Range of Applicability and Uncertainty

Thermal Conductivity

The Koyanagi model was derived from experiments in the temperature range 298.0 - 1073.15 K. The experimental uncertainty was not published.

The Stone model was derived from a 4th order polymonial fit from temperatures 273 - 1573 K. The uncertainty in the calculation was not published.

Specific Heat

The Snead model applies to temperatures 200 - 2400 K with uncertainty of $var element = document.getElementById("moose-equation-062ea77d-3de9-4801-8113-1da0a5f619a4");katex.render("\\pm", element, {displayMode:false,throwOnError:false});$ 7% for 200 $var element = document.getElementById("moose-equation-2f4dbf9b-8a9c-4179-9a86-4dbb3392c70d");katex.render("\\le", element, {displayMode:false,throwOnError:false});$ T $var element = document.getElementById("moose-equation-3828d307-bc90-4e8c-8ffe-8693e620c312");katex.render("\\le", element, {displayMode:false,throwOnError:false});$ 1000 K, and $var element = document.getElementById("moose-equation-04fd91d4-f277-4dd4-bdb3-5b940309516a");katex.render("\\pm", element, {displayMode:false,throwOnError:false});$ 4% for 1000 $var element = document.getElementById("moose-equation-9eb79013-3a83-41e3-993a-fb92adaf956e");katex.render("\\le", element, {displayMode:false,throwOnError:false});$ T $var element = document.getElementById("moose-equation-36612179-e563-4481-a373-e957fc84026d");katex.render("\\le", element, {displayMode:false,throwOnError:false});$ 2400 K.

The GA model was derived from experiments in the temperature range 150 - 2000 K. The experimental uncertainty was not published.

Example Input Syntax

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [thermalCompositeSiC]
    type = CompositeSiCThermal<<<{"description": "Computes thermal conductivity and specific heat of composite (CVI) SiC/SiC cladding.", "href": "CompositeSiCThermal.html"}>>>
    temperature<<<{"description": "Coupled Temperature"}>>> = temperature
    thermal_conductivity_model<<<{"description": "Options for the correlation used to calculate thermal conductivity"}>>> = STONE # This is the default
  []
[]

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.
specific_heat_modelSNEADOptions for the correlation used to calculate specific heat
Default:SNEAD
C++ Type:MooseEnum
Options:SNEAD, GA
Controllable:No
Description:Options for the correlation used to calculate specific heat
thermal_conductivity_modelSTONEOptions for the correlation used to calculate thermal conductivity
Default:STONE
C++ Type:MooseEnum
Options:KOYANAGI, COMBINED, STONE
Controllable:No
Description:Options for the correlation used to calculate thermal conductivity
thermal_conductivity_scale_factor1Thermal conductivity scale factor
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Thermal conductivity scale factor

Optional Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters

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

Outputs Parameters

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

Material Property Retrieval Parameters

Input Files

(test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_mieloszyk_ad.i)
(test/tests/thermalCompositeSiC/thermal_koyanagi.i)
(test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_katoh_ad.i)
(test/tests/thermalCompositeSiC/specific_heat_ga.i)
(assessment/LWR/validation/ATF_SiC_HFIR/analysis/base.i)
(test/tests/thermalCompositeSiC/thermal_stone_irradiated.i)
(test/tests/thermalCompositeSiC/thermal_koyanagi_ad.i)
(test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_mieloszyk.i)
(examples/accident_tolerant_fuel/u3si2_sic/u3si2_outer_monolith_1.5D.i)
(test/tests/thermalCompositeSiC/thermal_koyanagi_stone_combined.i)
(test/tests/thermalCompositeSiC/thermal_stone_unirradiated.i)
(test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_katoh.i)

References

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

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

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

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

L. L. Snead, T. Nozawa, Y. Katoh, T.-S. Byun, S. Kondo, and D. A. Petti. Handbook of sic properties for fuel performance modeling. Journal of Nuclear Materials, 371:329–377, 2007.

@article{snead2007,
    author = "Snead, L. L. and Nozawa, T. and Katoh, Y. and Byun, T.-S. and Kondo, S. and Petti, D. A.",
    title = "Handbook of SiC properties for fuel performance modeling",
    journal = "Journal of Nuclear Materials",
    volume = "371",
    year = "2007",
    pages = "329-377"
}

J. G. Stone, R. Schleicher, C. P. Deck, G. M. Jacobsen, H. E. Khalifa, and C. A. Back. Stress analysis and probabilistic assessment of multi-layer sic-based accident tolerant nuclear fuel cladding. Journal of Nuclear Materials, 466:682–697, 2015.

@article{stone2015,
    author = "Stone, J. G. and Schleicher, R. and Deck, C. P. and Jacobsen, G. M. and H. E. Khalifa and Back, C. A.",
    title = "Stress analysis and probabilistic assessment of multi-layer SiC-based accident tolerant nuclear fuel cladding",
    journal = "Journal of Nuclear Materials",
    volume = "466",
    pages = "682-697",
    year = "2015",
    publisher = "Elsevier"
}

(test/tests/thermalCompositeSiC/thermal_stone_irradiated.i)

# The mesh is a 1x1 square (single element).
# The temperature is held constant at 1073.15 K with a varying fast fluence
# from 0 to 1e25 n/m^2 over 1 second. The swelling is calculated by the Katoh
# volumetric swelling model. At 1073.15 K the calculated unirradiated thermal
# conductivity using the model of:
#
# J. G. Stone, R. Schleicher, C. P. Deck, G. M. Jacobsen, H. E. Khalifa, C. A.
# Back, "Stress analysis and probablistic assessment of multi-layer SiC-based
# accident tolerant nuclear fuel cladding," Journal of Nuclear Materials, 466,
# 2015, p. 682-697.
#
# is 14.9090 W/m-K. Then, the degradation due to swelling for composite silicon
# carbide is given by:
#
# 1 / k_irradiated = R_0 + R_irr
#
# where R_0 = 1 / k_unirradiated and R_irr = 15.11 * S where S is the swelling
# strain. The BISON results are compared to the analytical using the BISON
# values for swelling.
#
#  Fluence     Swelling      SiC/SiC k        SiC/SiC k
#   (dpa)     Strain (-)    BISON(W/m-K)   analytical(W/m-K)
#   1.0       5.91955e-3      6.38904          6.38904
#   2.0       7.20081e-3      5.68577          5.68577
#   3.0       7.70534e-3      5.44955          5.44955
#   4.0       7.90061e-3      5.36332          5.36332
#   5.0       7.97852e-3      5.32967          5.32967

[GlobalParams]
  displacements = 'disp_x disp_y'
  order = FIRST
  family = LAGRANGE
[]

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

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [temperature]
    initial_condition = 1073.15
  []
[]

[AuxVariables]
  [fast_neutron_fluence]
  []
  [swell]
    order = CONSTANT
    family = MONOMIAL
  []
  [thermal_conductivity]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [square]
    block = 0
    add_variables = false
    strain = FINITE
    eigenstrain_names = 'swelling_strain'
  []
[]

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

[AuxKernels]
  [thermal_conductivity]
    type = MaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 0
    execute_on = 'initial linear'
  []
  [fast_neutron_fluence]
    type = FunctionAux
    variable = fast_neutron_fluence
    function = fast_neutron_fluence
    block = 0
    execute_on = 'timestep_begin'
  []
  [total_swelling]
    type = MaterialRealAux
    variable = swell
    property = swelling
    block = 0
  []
[]

[Functions]
  [fast_neutron_fluence]
    type = PiecewiseLinear
    x = '0 1.0'
    y = '0 5e25'
  []
[]

[BCs]
  [temperature_all]
    type = DirichletBC
    boundary = 'left right bottom top'
    variable = temperature
    value = 1073.15
  []
  [leftX]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0.0
  []
  [bottomY]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom'
    value = 0.0
  []
[]

[Materials]
  [thermalCompositeSiC]
    type = CompositeSiCThermal
    temperature = temperature
    thermal_conductivity_model = STONE # This is the default
  []
  [density]
    type = StrainAdjustedDensity
    block = 0
    strain_free_density = 2700.0
  []
  [swelling]
    type = CompositeSiCVolumetricSwellingEigenstrain
    block = 0
    temperature = temperature
    swelling_model = KATOH
    number_of_substeps = 100
    fast_neutron_fluence = fast_neutron_fluence
    eigenstrain_name = swelling_strain
  []
  [elasticity_tensor]        # isotropic elasticity tensor
    type = ComputeIsotropicElasticityTensor
    block = 0
    youngs_modulus = 3e11
    poissons_ratio = 0.17
  []
  [elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 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'

   nl_rel_tol = 1e-5
   nl_abs_tol = 1e-10
   l_tol = 1e-5

   start_time = 0.0
   dt = 0.01
   end_time = 1.0
[]

[Postprocessors]
  [avg_temperature]
    type = SideAverageValue
    boundary = 'left right top bottom'
    execute_on = 'initial timestep_end'
    variable = temperature
  []
  [avg_th_cond]
    type = ElementAverageValue
    block = 0
    variable = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [swelling_out]
    type = ElementAverageValue
    variable = swell
    execute_on = 'initial timestep_end'
    block = 0
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_mieloszyk_ad.i)

# The purpose of this test is to verify the implementation of the Mieloszyk
# swelling model for composite and monolithic silicon carbide. The algorithm
# used is from:
#
# A. J. Mieloszyk, "Assessing Thermo-Mechanical Performance of ThO2 and SiC Clad
# Light Water Reactor Fuel Rods with a Modular Simulation Tool," PhD Dissertation,
# Massachusetts Institute of Technology, 2015.
#
# In the above reference an example is provided that applies a temperature and
# fluence profile over 300 days (here we simplify that to 300 seconds with the
# same profiles).  The reference also provides an expected profile of the
# volumetric swelling as a function of time given those profiles.  The BISON
# results from this test are compared to the results provided in Mieloszyk's
# dissertation in the swelling_mieloszyk_testing.xlsx Excel file contained
# within this directory.
#
# The small discrepancies between the BISON and Mieloszyk curves can be
# attributed to the number of points used in the digitization of
# Mieloszyk's data.
#
# The results of predicted swelling at selected times used in the temperature
# function of the test are tabulated below:
#
# Time (s)    Temperature (K)     Fluence (dpa)    Swelling Strain (-)
#   0             650.0                0                   0
#  100            650.0               2.0              0.014267
#  120            350.0               2.0              0.014267
#  200            750.0               4.2              0.011403
#  250            650.0               4.9              0.014604
#  300            550.0               5.6              0.017932

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

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

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

[AuxVariables]
  [fast_neutron_fluence]
    order = FIRST
    family = LAGRANGE
  []
  [swell]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [square]
    block = 0
    add_variables = false
    strain = FINITE
    eigenstrain_names = 'swelling_strain'
    use_automatic_differentiation = true
  []
[]

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

[Functions]
  [fast_neutron_fluence_func]
    type = PiecewiseLinear
    x = '0 100.0 120.0 200.0  300.0'
    y = '0 2e25  2e25  4.2e25 5.6e25'
  []
  [temperature]
    type = PiecewiseLinear
    x = '0   100 100.01 120 120.01 200.0 200.01 250 250.01 300'
    y = '650 650 350    350 675    750   600    650 625    550'
  []
[]

[AuxKernels]
  [fast_neutron_fluence]
    type = FunctionAux
    variable = fast_neutron_fluence
    function = fast_neutron_fluence_func
    execute_on = timestep_begin
  []
  [total_swelling]
    type = ADMaterialRealAux
    variable = swell
    property = swelling
    block = 0
  []
[]

[BCs]
  [leftX]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0.0
  []
  [bottomY]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom'
    value = 0.0
  []
  [rightX]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = right
    function = '1e-6*t'
  []
  [topY]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top
    function = '1e-6*t'
  []
  [all_temp]
    type = FunctionDirichletBC
    variable = temp
    boundary = 'bottom right top left'
    function = temperature
  []
[]

[Materials]
  [swelling]
    type = ADCompositeSiCVolumetricSwellingEigenstrain
    block = 0
    temperature = temp
    swelling_model = MIELOSZYK
    fast_neutron_fluence = fast_neutron_fluence
    eigenstrain_name = swelling_strain
  []
  [elasticity_tensor] # isotropic elasticity tensor
    type = ADComputeIsotropicElasticityTensor
    block = 0
    youngs_modulus = 3e11
    poissons_ratio = 0.17
  []
  [elastic_stress]
    type = ADComputeFiniteStrainElasticStress
    block = 0
  []
  [thermalCompositeSiC]
    type = ADCompositeSiCThermal
    temperature = temp
    thermal_conductivity_model = STONE # This is the default
  []
  [density]
    type = ADStrainAdjustedDensity
    block = 0
    strain_free_density = 3100.0
  []
[]

[Executioner]
  type = Transient

  solve_type = Newton

  petsc_options = '-snes_ksp_ew '
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu       superlu_dist'
  line_search = none
  l_max_its = 60
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-3
  l_tol = 1e-5

  start_time = 0.0
  dt = 1
  end_time = 300.0

[]

[Postprocessors]
  [avgFuelTemp]
    type = ElementAverageValue
    variable = temp
    execute_on = 'initial timestep_end'
    block = 0
  []
  [swelling_out]
    type = ElementAverageValue
    variable = swell
    execute_on = 'initial timestep_end'
    block = 0
  []
  [average_fluence]
    type = ElementAverageValue
    variable = fast_neutron_fluence
    execute_on = 'initial timestep_end'
    block = 0
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/thermalCompositeSiC/thermal_koyanagi.i)

# The mesh is a 1x1x1 cube (single element).
# The temperature is ramped on all faces of the cube from 300 K to 1800 K
# over 15.0 s.
#
# Specific heat capacity is computed using the correlation from:
#
# L. L. Snead, T. Nozawa, Y. Katoh, T-S. Byun, S. Kondo, D. A. Petti,
# "Handbook of SiC propeties for fuel performance modeling," Journal of Nuclear
# Materials, 371, 2007, p. 329-377.
#
# and thermal conductivity is computed using the thermal diffusivity data from:
#
# T. Koyanagi, Y. Katoh, G. Singh, M. Snead, "SiC/SiC Cladding Materials
# Properties Handbook," Technical Report ORNL/TM-2017/385, 2017.
#
# The thermal conductivity and specific heat computed by BISON for the
# temperatures shown below and compared with analytical solutions.
# The results are as follows:
#
#              SiC/SiC k       SiC/SiC k        SiC/SiC Cp      SiC/SiC Cp
# Temp(K)     BISON(W/m-K)  analytical(W/m-K)  BISON(J/kg-K)  analytical(J/kg-K)
#  300          8.20058         8.20058          676.721          676.721
#  600          7.91933         7.91933         1034.679         1034.679
#  900          7.47856         7.47856         1161.491         1161.491
# 1200          7.12582         7.12582         1241.972         1241.972
# 1500          7.45255         7.45255         1298.919         1298.919
# 1800          7.67649         7.67649         1337.951         1337.951

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

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

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

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

[AuxKernels]
  [thermal_conductivity]
    type = MaterialRealAux
    variable = thermal_conductivity
    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 15.0'
    y = '300 1800'
  []
[]

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

[Materials]
  [thermalCompositeSiC]
    type = CompositeSiCThermal
    temperature = temperature
    thermal_conductivity_model = KOYANAGI
  []
  [density]
    type = ParsedMaterial
    block = 0
    property_name = density
    expression = 2700.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-6
   nl_abs_tol = 1e-10
   l_tol = 1e-5
   start_time = 0.0
   num_steps = 15
   dt = 1.0
[]

[Postprocessors]
   [avg_temperature]
     type = SideAverageValue
     boundary = 'left right top bottom front back'
     execute_on = 'initial timestep_end'
     variable = temperature
   []
   [avg_th_cond]
     type = ElementAverageValue
     block = 0
     variable = thermal_conductivity
     execute_on = 'initial timestep_end'
   []
   [avg_specific_heat]
     type = ElementAverageValue
     block = 0
     variable = specific_heat
     execute_on = 'initial timestep_end'
   []
[]

[Outputs]
  exodus = true
[]

(test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_katoh_ad.i)

# This test tests the Katoh swelling model used for composite and monolithic
# silicon carbide.  The model is obtained from:
#
# Y. Katoh, K. Ozawa, C. Shih, T. Nozawa, R. J. Shinavski, A. Hasegawa, and
# L. L. Snead. Continuous SiC fiber, CVI SiC matrix composites for nuclear
# applications: properties and irradiation effects. Journal of Nuclear Materials,
# 448, pp. 448-476, 2014.
#
# The swelling rate is a function of the fast neutron fluence:
#
# dS/dgamma = k(T) * gamma^(-1/3) * exp(-gamma/gamma_sc(T))
#
# where gamma is the fluence in dpa, S is the swelling strain, and k(T) and
# gamma_sc(T) are temperature dependent constants.  The above equation cannot
# be analytically integrated and an incremental approach is used to calculate the
# swelling strain. The "analytical" solution is taken as the results of the
# numerical integral calculated using Wolfram Alpha or Mathematica. Depending
# upon the number of substeps used the BISON simulations approach this analytical
# solution.
#
# Here, a single 1x1 m 2D element is subjected to a constant temperature of
# 1073.15 K with a varying fast neutron fluence from 0 to 1e25 n/m^2 (0 to 5 dpa)
# over 1 second. Below BISON results are included for 100, 1000 and 10000 substeps to
# illustrate that the results approach that of the analytical solution. This test
# is for the 100 substep case. Tests for the 1000 and 10000 substeps at high
# temperature as well as low and mid range temperatures of 473.15 K
# and 773.15 K (with 100 substeps) are included using cli_args in the tests file.
#
# Fluence               BISON Swelling Strain (-)                 Analytical
#  (dpa)        100 substeps   1000 substeps  10000 substeps       Swelling
#   1.0          5.91955e-3      5.83353e-3     5.83899e-3        5.84045e-3
#   2.0          7.20081e-3      7.22617e-3     7.23163e-3        7.23270e-3
#   3.0          7.70534e-3      7.73070e-3     7.73617e-3        7.73713e-3
#   4.0          7.90061e-3      7.92597e-3     7.93143e-3        7.93235e-3
#   5.0          7.97852e-3      8.00388e-3     8.00934e-3        8.01025e-3
#

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

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

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

[AuxVariables]
  [fast_neutron_fluence]
    order = FIRST
    family = LAGRANGE
  []
  [swell]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [square]
    block = 0
    add_variables = false
    strain = FINITE
    eigenstrain_names = 'swelling_strain'
    use_automatic_differentiation = true
  []
[]

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

  [heat_ie]
    type = ADHeatConductionTimeDerivative
    variable = temp
  []
[]

[Functions]
  [fast_neutron_fluence_func]
    type = PiecewiseLinear
    x = '0 1.0'
    y = '0 5e25'
  []
[]

[AuxKernels]
  [fast_neutron_fluence]
    type = FunctionAux
    variable = fast_neutron_fluence
    function = fast_neutron_fluence_func
    execute_on = timestep_begin
  []
  [total_swelling]
    type = MaterialRealAux
    variable = swell
    property = swelling
    block = 0
  []
[]

[BCs]
  [leftX]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0.0
  []
  [bottomY]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom'
    value = 0.0
  []
  [all_temp]
    type = DirichletBC
    variable = temp
    boundary = 'bottom left right top'
    value = 1073.15
  []
[]

[Materials]
  [swelling]
    type = ADCompositeSiCVolumetricSwellingEigenstrain
    block = 0
    temperature = temp
    swelling_model = KATOH
    number_of_substeps = 100
    fast_neutron_fluence = fast_neutron_fluence
    eigenstrain_name = swelling_strain
  []
  [elasticity_tensor] # isotropic elasticity tensor
    type = ADComputeIsotropicElasticityTensor
    block = 0
    youngs_modulus = 3e11
    poissons_ratio = 0.17
  []
  [elastic_stress]
    type = ADComputeFiniteStrainElasticStress
    block = 0
  []
  [thermalCompositeSiC]
    type = ADCompositeSiCThermal
    temperature = temp
    thermal_conductivity_model = STONE # This is the default
  []
  [density]
    type = ADStrainAdjustedDensity
    block = 0
    density = 3100.0
  []
[]

[Executioner]
  type = Transient

  solve_type = Newton

  petsc_options = '-snes_ksp_ew '
  petsc_options_iname = '-ksp_gmres_restart -pc_type -pc_hypre_type'
  petsc_options_value = '70 hypre boomeramg'
  l_max_its = 60
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-10
  l_tol = 1e-5

  start_time = 0.0
  dt = 0.01
  end_time = 1.0

[]

[Postprocessors]
  [avgFuelTemp]
    type = ElementAverageValue
    variable = temp
    execute_on = 'initial timestep_end'
    block = 0
  []
  [swelling_out]
    type = ElementAverageValue
    variable = swell
    execute_on = 'initial timestep_end'
    block = 0
  []
  [average_fluence]
    type = ElementAverageValue
    variable = fast_neutron_fluence
    execute_on = 'initial timestep_end'
    block = 0
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/thermalCompositeSiC/specific_heat_ga.i)

# The mesh is a 1x1x1 cube (single element).
# The temperature is ramped on all faces of the cube from 300 K to 1800 K
# over 15.0 s.
#
# Specific heat capacity is computed using the correlation from:
#
# G. Jacobsen, H. Chiger, K. Shapovalov, "Materials property manual: General
# Atomics silicon carbide cladding," GA Report GA-A28712X, Revision B1, 2020.
#
# and thermal conductivity is computed using the thermal diffusivity data from:
#
# T. Koyanagi, Y. Katoh, G. Singh, M. Snead, "SiC/SiC Cladding Materials
# Properties Handbook," Technical Report ORNL/TM-2017/385, 2017.
#
# The thermal conductivity and specific heat computed by BISON for the
# temperatures shown below and compared with analytical solutions.
# The results are as follows:
#
#              SiC/SiC k       SiC/SiC k        SiC/SiC Cp      SiC/SiC Cp
# Temp(K)     BISON(W/m-K)  analytical(W/m-K)  BISON(J/kg-K)  analytical(J/kg-K)
#  300          7.81553         7.81553          644.946          644.946
#  600          7.97935         7.97935         1042.540         1042.540
#  900          7.54392         7.54392         1171.642         1171.642
# 1200          7.14583         7.14583         1245.461         1245.461
# 1500          7.51290         7.51290         1309.438         1309.438
# 1800          7.90129         7.90129         1377.131         1377.131

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

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

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

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

[AuxKernels]
  [thermal_conductivity]
    type = MaterialRealAux
    variable = thermal_conductivity
    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 15.0'
    y = '300 1800'
  []
[]

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

[Materials]
  [thermalCompositeSiC]
    type = CompositeSiCThermal
    temperature = temperature
    thermal_conductivity_model = KOYANAGI
    specific_heat_model = GA
  []
  [density]
    type = ParsedMaterial
    block = 0
    property_name = density
    expression = 2700.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-6
   nl_abs_tol = 1e-10
   l_tol = 1e-5
   start_time = 0.0
   num_steps = 15
   dt = 1.0
[]

[Postprocessors]
   [avg_temperature]
     type = SideAverageValue
     boundary = 'left right top bottom front back'
     execute_on = 'initial timestep_end'
     variable = temperature
   []
   [avg_th_cond]
     type = ElementAverageValue
     block = 0
     variable = thermal_conductivity
     execute_on = 'initial timestep_end'
   []
   [avg_specific_heat]
     type = ElementAverageValue
     block = 0
     variable = specific_heat
     execute_on = 'initial timestep_end'
   []
[]

[Outputs]
  exodus = true
[]

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[PerformanceMetricOutputs]
[]

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

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

  line_search = 'none'

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

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

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

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

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

(test/tests/thermalCompositeSiC/thermal_stone_irradiated.i)

# The mesh is a 1x1 square (single element).
# The temperature is held constant at 1073.15 K with a varying fast fluence
# from 0 to 1e25 n/m^2 over 1 second. The swelling is calculated by the Katoh
# volumetric swelling model. At 1073.15 K the calculated unirradiated thermal
# conductivity using the model of:
#
# J. G. Stone, R. Schleicher, C. P. Deck, G. M. Jacobsen, H. E. Khalifa, C. A.
# Back, "Stress analysis and probablistic assessment of multi-layer SiC-based
# accident tolerant nuclear fuel cladding," Journal of Nuclear Materials, 466,
# 2015, p. 682-697.
#
# is 14.9090 W/m-K. Then, the degradation due to swelling for composite silicon
# carbide is given by:
#
# 1 / k_irradiated = R_0 + R_irr
#
# where R_0 = 1 / k_unirradiated and R_irr = 15.11 * S where S is the swelling
# strain. The BISON results are compared to the analytical using the BISON
# values for swelling.
#
#  Fluence     Swelling      SiC/SiC k        SiC/SiC k
#   (dpa)     Strain (-)    BISON(W/m-K)   analytical(W/m-K)
#   1.0       5.91955e-3      6.38904          6.38904
#   2.0       7.20081e-3      5.68577          5.68577
#   3.0       7.70534e-3      5.44955          5.44955
#   4.0       7.90061e-3      5.36332          5.36332
#   5.0       7.97852e-3      5.32967          5.32967

[GlobalParams]
  displacements = 'disp_x disp_y'
  order = FIRST
  family = LAGRANGE
[]

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

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [temperature]
    initial_condition = 1073.15
  []
[]

[AuxVariables]
  [fast_neutron_fluence]
  []
  [swell]
    order = CONSTANT
    family = MONOMIAL
  []
  [thermal_conductivity]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [square]
    block = 0
    add_variables = false
    strain = FINITE
    eigenstrain_names = 'swelling_strain'
  []
[]

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

[AuxKernels]
  [thermal_conductivity]
    type = MaterialRealAux
    variable = thermal_conductivity
    property = thermal_conductivity
    block = 0
    execute_on = 'initial linear'
  []
  [fast_neutron_fluence]
    type = FunctionAux
    variable = fast_neutron_fluence
    function = fast_neutron_fluence
    block = 0
    execute_on = 'timestep_begin'
  []
  [total_swelling]
    type = MaterialRealAux
    variable = swell
    property = swelling
    block = 0
  []
[]

[Functions]
  [fast_neutron_fluence]
    type = PiecewiseLinear
    x = '0 1.0'
    y = '0 5e25'
  []
[]

[BCs]
  [temperature_all]
    type = DirichletBC
    boundary = 'left right bottom top'
    variable = temperature
    value = 1073.15
  []
  [leftX]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0.0
  []
  [bottomY]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom'
    value = 0.0
  []
[]

[Materials]
  [thermalCompositeSiC]
    type = CompositeSiCThermal
    temperature = temperature
    thermal_conductivity_model = STONE # This is the default
  []
  [density]
    type = StrainAdjustedDensity
    block = 0
    strain_free_density = 2700.0
  []
  [swelling]
    type = CompositeSiCVolumetricSwellingEigenstrain
    block = 0
    temperature = temperature
    swelling_model = KATOH
    number_of_substeps = 100
    fast_neutron_fluence = fast_neutron_fluence
    eigenstrain_name = swelling_strain
  []
  [elasticity_tensor]        # isotropic elasticity tensor
    type = ComputeIsotropicElasticityTensor
    block = 0
    youngs_modulus = 3e11
    poissons_ratio = 0.17
  []
  [elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 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'

   nl_rel_tol = 1e-5
   nl_abs_tol = 1e-10
   l_tol = 1e-5

   start_time = 0.0
   dt = 0.01
   end_time = 1.0
[]

[Postprocessors]
  [avg_temperature]
    type = SideAverageValue
    boundary = 'left right top bottom'
    execute_on = 'initial timestep_end'
    variable = temperature
  []
  [avg_th_cond]
    type = ElementAverageValue
    block = 0
    variable = thermal_conductivity
    execute_on = 'initial timestep_end'
  []
  [swelling_out]
    type = ElementAverageValue
    variable = swell
    execute_on = 'initial timestep_end'
    block = 0
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/thermalCompositeSiC/thermal_koyanagi_ad.i)

# The mesh is a 1x1x1 cube (single element).
# The temperature is ramped on all faces of the cube from 300 K to 1800 K
# over 15.0 s.
#
# Specific heat capacity is computed using the correlation from:
#
# L. L. Snead, T. Nozawa, Y. Katoh, T-S. Byun, S. Kondo, D. A. Petti,
# "Handbook of SiC propeties for fuel performance modeling," Journal of Nuclear
# Materials, 371, 2007, p. 329-377.
#
# and thermal conductivity is computed using the thermal diffusivity data from:
#
# T. Koyanagi, Y. Katoh, G. Singh, M. Snead, "SiC/SiC Cladding Materials
# Properties Handbook," Technical Report ORNL/TM-2017/385, 2017.
#
# The thermal conductivity and specific heat computed by BISON for the
# temperatures shown below and compared with analytical solutions.
# The results are as follows:
#
#              SiC/SiC k       SiC/SiC k        SiC/SiC Cp      SiC/SiC Cp
# Temp(K)     BISON(W/m-K)  analytical(W/m-K)  BISON(J/kg-K)  analytical(J/kg-K)
#  300          8.20058         8.20058          676.721          676.721
#  600          7.91933         7.91933         1034.679         1034.679
#  900          7.47856         7.47856         1161.491         1161.491
# 1200          7.12582         7.12582         1241.972         1241.972
# 1500          7.45255         7.45255         1298.919         1298.919
# 1800          7.67649         7.67649         1337.951         1337.951

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

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

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

[Kernels]
  [heat]
    type = ADHeatConduction
    variable = temperature
  []
  [heat_ie]
    type = ADHeatConductionTimeDerivative
    variable = temperature
  []
[]

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

[Functions]
  [temp_ramp]
    type = PiecewiseLinear
    x = '0.0 15.0'
    y = '300 1800'
  []
[]

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

[Materials]
  [thermalCompositeSiC]
    type = ADCompositeSiCThermal
    temperature = temperature
    thermal_conductivity_model = KOYANAGI
  []
  [density]
    type = ADParsedMaterial
    block = 0
    property_name = density
    expression = 2700.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-6
   nl_abs_tol = 1e-10
   l_tol = 1e-5
   start_time = 0.0
   num_steps = 15
   dt = 1.0
[]

[Postprocessors]
   [avg_temperature]
     type = SideAverageValue
     boundary = 'left right top bottom front back'
     execute_on = 'initial timestep_end'
     variable = temperature
   []
   [avg_th_cond]
     type = ElementAverageValue
     block = 0
     variable = thermal_conductivity
     execute_on = 'initial timestep_end'
   []
   [avg_specific_heat]
     type = ElementAverageValue
     block = 0
     variable = specific_heat
     execute_on = 'initial timestep_end'
   []
[]

[Outputs]
  exodus = true
[]

(test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_mieloszyk.i)

# The purpose of this test is to verify the implementation of the Mieloszyk
# swelling model for composite and monolithic silicon carbide. The algorithm
# used is from:
#
# A. J. Mieloszyk, "Assessing Thermo-Mechanical Performance of ThO2 and SiC Clad
# Light Water Reactor Fuel Rods with a Modular Simulation Tool," PhD Dissertation,
# Massachusetts Institute of Technology, 2015.
#
# In the above reference an example is provided that applies a temperature and
# fluence profile over 300 days (here we simplify that to 300 seconds with the
# same profiles).  The reference also provides an expected profile of the
# volumetric swelling as a function of time given those profiles.  The BISON
# results from this test are compared to the results provided in Mieloszyk's
# dissertation in the swelling_mieloszyk_testing.xlsx Excel file contained
# within this directory.
#
# The small discrepancies between the BISON and Mieloszyk curves can be
# attributed to the number of points used in the digitization of
# Mieloszyk's data.
#
# The results of predicted swelling at selected times used in the temperature
# function of the test are tabulated below:
#
# Time (s)    Temperature (K)     Fluence (dpa)    Swelling Strain (-)
#   0             650.0                0                   0
#  100            650.0               2.0              0.014267
#  120            350.0               2.0              0.014267
#  200            750.0               4.2              0.011403
#  250            650.0               4.9              0.014604
#  300            550.0               5.6              0.017932

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

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

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

[AuxVariables]
  [fast_neutron_fluence]
    order = FIRST
    family = LAGRANGE
  []
  [swell]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [square]
    block = 0
    add_variables = false
    strain = FINITE
    eigenstrain_names = 'swelling_strain'
  []
[]

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

[Functions]
  [fast_neutron_fluence_func]
    type = PiecewiseLinear
    x = '0 100.0 120.0 200.0  300.0'
    y = '0 2e25  2e25  4.2e25 5.6e25'
  []
  [temperature]
    type = PiecewiseLinear
    x = '0   100 100.01 120 120.01 200.0 200.01 250 250.01 300'
    y = '650 650 350    350 675    750   600    650 625    550'
  []
[]

[AuxKernels]
  [fast_neutron_fluence]
    type = FunctionAux
    variable = fast_neutron_fluence
    function = fast_neutron_fluence_func
    execute_on = timestep_begin
  []
  [total_swelling]
    type = MaterialRealAux
    variable = swell
    property = swelling
    block = 0
  []
[]

[BCs]
  [leftX]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0.0
  []
  [bottomY]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom'
    value = 0.0
  []
  [rightX]
    type = FunctionDirichletBC
    variable = disp_x
    boundary = right
    function = '1e-6*t'
  []
  [topY]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top
    function = '1e-6*t'
  []
  [all_temp]
    type = FunctionDirichletBC
    variable = temp
    boundary = 'bottom right top left'
    function = temperature
  []
[]

[Materials]
  [swelling]
    type = CompositeSiCVolumetricSwellingEigenstrain
    block = 0
    temperature = temp
    swelling_model = MIELOSZYK
    fast_neutron_fluence = fast_neutron_fluence
    eigenstrain_name = swelling_strain
  []
  [elasticity_tensor] # isotropic elasticity tensor
    type = ComputeIsotropicElasticityTensor
    block = 0
    youngs_modulus = 3e11
    poissons_ratio = 0.17
  []
  [elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 0
  []
  [thermalCompositeSiC]
    type = CompositeSiCThermal
    temperature = temp
    thermal_conductivity_model = STONE # This is the default
  []
  [density]
    type = StrainAdjustedDensity
    block = 0
    strain_free_density = 3100.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'
  line_search = none
  l_max_its = 60
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-3
  l_tol = 1e-5

  start_time = 0.0
  dt = 1
  end_time = 300.0

[]

[Postprocessors]
  [avgFuelTemp]
    type = ElementAverageValue
    variable = temp
    execute_on = 'initial timestep_end'
    block = 0
  []
  [swelling_out]
    type = ElementAverageValue
    variable = swell
    execute_on = 'initial timestep_end'
    block = 0
  []
  [average_fluence]
    type = ElementAverageValue
    variable = fast_neutron_fluence
    execute_on = 'initial timestep_end'
    block = 0
  []
[]

[Outputs]
  exodus = true
[]

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

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


initial_fuel_density = 11590.0

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

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

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

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

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

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

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

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

[Physics]
  [SolidMechanics]
    [Layered1D]
      [fuel]
        block = fuel
        add_variables = true
        strain = SMALL
        incremental = true
        add_scalar_variables = true
        out_of_plane_strain_name = strain_yy
        fuel_pin_geometry = pin_geometry
        out_of_plane_pressure_function = fuel_axial_pressure
        eigenstrain_names = 'fuel_thermal_strain fuel_swelling_strain'
        generate_output = 'stress_xx stress_yy stress_zz vonmises_stress strain_xx'
        extra_vector_tags = 'ref'
        group_scalar_vars_in_reference_residual = true
        mesh_generator = layered1D_mesh
      []
      [composite]
        block = clad
        add_variables = true
        strain = SMALL
        incremental = true
        add_scalar_variables = true
        out_of_plane_strain_name = strain_yy
        fuel_pin_geometry = pin_geometry
        out_of_plane_pressure_function = clad_axial_pressure
        eigenstrain_names = 'composite_thermal_strain composite_swelling_strain'
        generate_output = 'stress_xx stress_yy stress_zz vonmises_stress strain_xx'
        extra_vector_tags = 'ref'
        group_scalar_vars_in_reference_residual = true
        mesh_generator = layered1D_mesh
      []
      [monolith]
        block = monolithic_layer
        add_variables = true
        strain = SMALL
        incremental = true
        add_scalar_variables = true
        out_of_plane_strain_name = strain_yy
        fuel_pin_geometry = pin_geometry
        out_of_plane_pressure_function = clad_axial_pressure
        eigenstrain_names = 'monolith_thermal_strain monolith_swelling_strain'
        generate_output = 'stress_xx stress_yy stress_zz vonmises_stress strain_xx'
        extra_vector_tags = 'ref'
        group_scalar_vars_in_reference_residual = true
        mesh_generator = layered1D_mesh
      []
    []
  []
[]

[Burnup]
  [burnup]
    block = fuel
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    num_radial = 80
    num_axial = 11
    order = CONSTANT
    family = MONOMIAL
    fuel_pin_geometry = pin_geometry
    fuel_volume_ratio = 1.0
    RPF = RPF
    fuel_type = U3Si2
  []
[]

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

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

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

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

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

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

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

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

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

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

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

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

  line_search = 'none'

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

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

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

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

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

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

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

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

(test/tests/thermalCompositeSiC/thermal_koyanagi_stone_combined.i)

# The mesh is a 1x1x1 cube (single element).
# The temperature is ramped on all faces of the cube from 300 K to 1800 K
# over 15.0 s.
#
#
# Unirradiated thermal conductivity is computed using the thermal diffusivity data from:
#
# T. Koyanagi, Y. Katoh, G. Singh, M. Snead, "SiC/SiC Cladding Materials
# Properties Handbook," Technical Report ORNL/TM-2017/385, 2017.
#
# Irradiation resistance is computed according to:
#
# J. G. Stone, R. Schleicher, C. P. Deck, G. M. Jacobsen, H. E. Khalifa, C. A.
# Back, "Stress analysis and probablistic assessment of multi-layer SiC-based
# accident tolerant nuclear fuel cladding," Journal of Nuclear Materials, 466,
# 2015, p. 682-697.
#
# 1 / k_irradiated = R_0 + R_irr
#
# where R_0 = 1 / k_unirradiated and R_irr = 15.11 * S where S is the swelling
# strain.
#
# The thermal conductivity and specific heat computed by BISON for the
# temperatures shown below and compared with analytical solutions.
# The results are as follows:
#
#              SiC/SiC k       SiC/SiC k        SiC/SiC Cp      SiC/SiC Cp
# Temp(K)     BISON(W/m-K)  analytical(W/m-K)  BISON(J/kg-K)  analytical(J/kg-K)
# 300         8.20058       8.20058            676.721         676.721
# 600         1.67342       1.66730            1034.698        1034.697
# 900         1.62824       1.63774            1161.491        1161.491
# 1200        1.59471       1.60440            1241.972        1241.972
# 1800        1.44797       1.45651            1337.951        1337.9509

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

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

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

[AuxVariables]
  [fast_neutron_fluence]
  []
  [swell]
    order = CONSTANT
    family = MONOMIAL
  []
  [specific_heat]
    order = CONSTANT
    family = MONOMIAL
  []
  [thermal_conductivity]
    order = CONSTANT
    family = MONOMIAL
  []
[]
[Physics/SolidMechanics/QuasiStatic]
  [block]
    block = 0
    add_variables = false
    strain = FINITE
    eigenstrain_names = 'swelling_strain'
  []
[]

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

[AuxKernels]
  [fast_neutron_fluence]
    type = FunctionAux
    variable = fast_neutron_fluence
    function = fast_neutron_fluence
    block = 0
    execute_on = 'timestep_begin'
  []
  [total_swelling]
    type = MaterialRealAux
    variable = swell
    property = swelling
    block = 0
  []
  [thermal_conductivity]
    type = MaterialRealAux
    variable = thermal_conductivity
    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]
  [fast_neutron_fluence]
    type = PiecewiseLinear
    x = '0 5.0'
    y = '0 5e25'
  []
  [temp_ramp]
    type = PiecewiseLinear
    x = '0.0 5.0'
    y = '300 1800'
  []
[]

[BCs]
  [temperature_all]
    type = FunctionDirichletBC
    boundary = 'left right front back bottom top'     # All cube faces
    variable = temperature
    function = temp_ramp
  []
  [leftX]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0.0
  []
  [frontY]
    type = DirichletBC
    variable = disp_y
    boundary = 'front'
    value = 0.0
  []
  [bottomZ]
    type = DirichletBC
    variable = disp_z
    boundary = 'bottom'
    value = 0.0
  []
[]

[Materials]
  [thermalCompositeSiC]
    type = CompositeSiCThermal
    temperature = temperature
    thermal_conductivity_model = COMBINED
  []
  [density]
    type = StrainAdjustedDensity
    block = 0
    strain_free_density = 2700.0
  []
  [swelling]
    type = CompositeSiCVolumetricSwellingEigenstrain
    block = 0
    temperature = temperature
    swelling_model = KATOH
    number_of_substeps = 100
    fast_neutron_fluence = fast_neutron_fluence
    eigenstrain_name = swelling_strain
  []
  [elasticity_tensor]        # isotropic elasticity tensor
    type = ComputeIsotropicElasticityTensor
    block = 0
    youngs_modulus = 3e11
    poissons_ratio = 0.17
  []
  [elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 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-5
   nl_abs_tol = 1e-10
   l_tol = 1e-5
   start_time = 0.0
   num_steps = 10
   dt = 0.5
[]

[Postprocessors]
   [avg_temperature]
     type = SideAverageValue
     boundary = 'left right top bottom front back'
     execute_on = 'initial timestep_end'
     variable = temperature
   []
   [avg_th_cond]
     type = ElementAverageValue
     block = 0
     variable = thermal_conductivity
     execute_on = 'initial timestep_end'
   []
   [avg_specific_heat]
     type = ElementAverageValue
     block = 0
     variable = specific_heat
     execute_on = 'initial timestep_end'
   []
   [swelling_out]
     type = ElementAverageValue
     variable = swell
     execute_on = 'initial timestep_end'
     block = 0
   []
[]

[Outputs]
  exodus = true
[]

(test/tests/thermalCompositeSiC/thermal_stone_unirradiated.i)

# The mesh is a 1x1x1 cube (single element).
# The temperature is ramped on all faces of the cube from 300 K to 1500 K
# over 12.0 s.
#
# Specific heat capacity is computed using the correlation from:
#
# L. L. Snead, T. Nozawa, Y. Katoh, T-S. Byun, S. Kondo, D. A. Petti,
# "Handbook of SiC propeties for fuel performance modeling," Journal of Nuclear
# Materials, 371, 2007, p. 329-377.
#
# and thermal conductivity is computed using the thermal conductivity correlation
# from:
#
# J. G. Stone, R. Schleicher, C. P. Deck, G. M. Jacobsen, H. E. Khalifa, C. A.
# Back, "Stress analysis and probablistic assessment of multi-layer SiC-based
# accident tolerant nuclear fuel cladding," Journal of Nuclear Materials, 466,
# 2015, p. 682-697.
#
# The specific heat correlation is tested elsewhere.  The results of the BISON
# simulation for unirradiatedthermal conductivity is compared to the analytical
# solution below.
#
#              SiC/SiC k       SiC/SiC k
# Temp(K)     BISON(W/m-K)  analytical(W/m-K)
#  300          18.21599        18.21599
#  600          18.62984        18.62984
#  900          16.13219        16.13219
# 1200          14.31944        14.31944
# 1500          13.46375        13.46375

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

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

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

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

[AuxKernels]
  [thermal_conductivity]
    type = MaterialRealAux
    variable = thermal_conductivity
    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 12.0'
    y = '300 1500'
  []
[]

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

[Materials]
  [thermalCompositeSiC]
    type = CompositeSiCThermal
    temperature = temperature
    thermal_conductivity_model = STONE # This is the default
  []
  [density]
    type = ParsedMaterial
    block = 0
    property_name = density
    expression = 2700.0
  []
  [swelling]
    type = CompositeSiCVolumetricSwellingEigenstrain
    block = 0
    temperature = temperature
    fast_neutron_fluence = 0.0
    eigenstrain_name = swelling_strain
  []
[]

[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-6
   nl_abs_tol = 1e-10
   l_tol = 1e-5
   start_time = 0.0
   num_steps = 12
   dt = 1.0
[]

[Postprocessors]
   [avg_temperature]
     type = SideAverageValue
     boundary = 'left right top bottom front back'
     execute_on = 'initial timestep_end'
     variable = temperature
   []
   [avg_th_cond]
     type = ElementAverageValue
     block = 0
     variable = thermal_conductivity
     execute_on = 'initial timestep_end'
   []
   [avg_specific_heat]
     type = ElementAverageValue
     block = 0
     variable = specific_heat
     execute_on = 'initial timestep_end'
   []
[]

[Outputs]
  exodus = true
[]

(test/tests/solid_mechanics/compositeSiC_eigenstrains/volumetric_swelling/swelling_katoh.i)

# This test tests the Katoh swelling model used for composite and monolithic
# silicon carbide.  The model is obtained from:
#
# Y. Katoh, K. Ozawa, C. Shih, T. Nozawa, R. J. Shinavski, A. Hasegawa, and
# L. L. Snead. Continuous SiC fiber, CVI SiC matrix composites for nuclear
# applications: properties and irradiation effects. Journal of Nuclear Materials,
# 448, pp. 448-476, 2014.
#
# The swelling rate is a function of the fast neutron fluence:
#
# dS/dgamma = k(T) * gamma^(-1/3) * exp(-gamma/gamma_sc(T))
#
# where gamma is the fluence in dpa, S is the swelling strain, and k(T) and
# gamma_sc(T) are temperature dependent constants.  The above equation cannot
# be analytically integrated and an incremental approach is used to calculate the
# swelling strain. The "analytical" solution is taken as the results of the
# numerical integral calculated using Wolfram Alpha or Mathematica. Depending
# upon the number of substeps used the BISON simulations approach this analytical
# solution.
#
# Here, a single 1x1 m 2D element is subjected to a constant temperature of
# 1073.15 K with a varying fast neutron fluence from 0 to 1e25 n/m^2 (0 to 5 dpa)
# over 1 second. Below BISON results are included for 100, 1000 and 10000 substeps to
# illustrate that the results approach that of the analytical solution. This test
# is for the 100 substep case. Tests for the 1000 and 10000 substeps at high
# temperature as well as low and mid range temperatures of 473.15 K
# and 773.15 K (with 100 substeps) are included using cli_args in the tests file.
#
# Fluence               BISON Swelling Strain (-)                 Analytical
#  (dpa)        100 substeps   1000 substeps  10000 substeps       Swelling
#   1.0          5.91955e-3      5.83353e-3     5.83899e-3        5.84045e-3
#   2.0          7.20081e-3      7.22617e-3     7.23163e-3        7.23270e-3
#   3.0          7.70534e-3      7.73070e-3     7.73617e-3        7.73713e-3
#   4.0          7.90061e-3      7.92597e-3     7.93143e-3        7.93235e-3
#   5.0          7.97852e-3      8.00388e-3     8.00934e-3        8.01025e-3
#

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

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

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

[AuxVariables]
  [fast_neutron_fluence]
    order = FIRST
    family = LAGRANGE
  []
  [swell]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [square]
    block = 0
    add_variables = false
    strain = FINITE
    eigenstrain_names = 'swelling_strain'
  []
[]

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

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

[Functions]
  [fast_neutron_fluence_func]
    type = PiecewiseLinear
    x = '0 1.0'
    y = '0 5e25'
  []
[]

[AuxKernels]
  [fast_neutron_fluence]
    type = FunctionAux
    variable = fast_neutron_fluence
    function = fast_neutron_fluence_func
    execute_on = timestep_begin
  []
  [total_swelling]
    type = MaterialRealAux
    variable = swell
    property = swelling
    block = 0
  []
[]


[BCs]
  [leftX]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0.0
  []
  [bottomY]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom'
    value = 0.0
  []
  [all_temp]
    type = DirichletBC
    variable = temp
    boundary = 'bottom left right top'
    value = 1073.15
  []
[]

[Materials]
  [swelling]
    type = CompositeSiCVolumetricSwellingEigenstrain
    block = 0
    temperature = temp
    swelling_model = KATOH
    number_of_substeps = 100
    fast_neutron_fluence = fast_neutron_fluence
    eigenstrain_name = swelling_strain
  []
  [elasticity_tensor]        # isotropic elasticity tensor
    type = ComputeIsotropicElasticityTensor
    block = 0
    youngs_modulus = 3e11
    poissons_ratio = 0.17
  []
  [elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 0
  []
  [thermalCompositeSiC]
    type = CompositeSiCThermal
    temperature = temp
    thermal_conductivity_model = STONE # This is the default
  []
  [density]
    type = StrainAdjustedDensity
    block = 0
    strain_free_density = 3100.0
  []
[]

[Executioner]
  type = Transient

  solve_type = 'PJFNK'

  petsc_options = '-snes_ksp_ew '
  petsc_options_iname = '-ksp_gmres_restart -pc_type -pc_hypre_type'
  petsc_options_value = '70 hypre boomeramg'
  l_max_its = 60
  nl_rel_tol = 1e-5
  nl_abs_tol = 1e-10
  l_tol = 1e-5

  start_time = 0.0
  dt = 0.01
  end_time = 1.0

[]

[Postprocessors]
  [avgFuelTemp]
    type = ElementAverageValue
    variable = temp
    execute_on = 'initial timestep_end'
    block = 0
  []
  [swelling_out]
    type = ElementAverageValue
    variable = swell
    execute_on = 'initial timestep_end'
    block = 0
  []
  [average_fluence]
    type = ElementAverageValue
    variable = fast_neutron_fluence
    execute_on = 'initial timestep_end'
    block = 0
  []
[]

[Outputs]
  exodus = true
[]

Description
Thermal Conductivity Models
Specific Heat Capacity Models
Range of Applicability and Uncertainty
Example Input Syntax
Input Parameters
Input Files
References