ZirconiumDiffusion

Computes the amount of zirconium that is transported across the mesh.

Description

After all the diffusivities are set, ZirconiumDiffusion is used to transport the zirconium across the mesh.

A pseudo-binary model for U-19Pu-10Zr fuel was implemented into BISON with two essential numerical modifications: addition of an artificial diffusion term in the 2-phase region to stabilize the algorithm, and smoothing of the model coefficients near phase diagram curves. We solve the coupled system of equations, $(1)var element = document.getElementById("moose-equation-c81ec518-1372-40cd-9ba5-31ab10f92df9");katex.render(" \\begin{gathered} \\rho c_p \\frac{\\partial T}{\\partial t} = \\nabla \\cdot k(x,T) \\nabla T + q ( \\mathbf{r}, x) \\\\ \\frac{\\partial x}{\\partial t} = \\nabla \\cdot D(x,T) \\nabla x + \\nabla \\cdot S(x,T) \\nabla T \\end{gathered}", element, {displayMode:true,throwOnError:false});$ for an evolving temperature distribution $var element = document.getElementById("moose-equation-546cddd7-86bd-4d7d-8309-d9153f179200");katex.render("T(\\mathbf{r},t)", element, {displayMode:false,throwOnError:false});$ and Zr atom fraction distribution $var element = document.getElementById("moose-equation-57fa9258-4f32-4ff8-8d13-0232f4d786fd");katex.render("{{x}_{\\text{Zr}}}=x(\\mathbf{r},t)", element, {displayMode:false,throwOnError:false});$ . This is a departure from the model considered in Kim et al. (2004) and Kim et al. (2006) which used a fixed temperature profile. We use a variable heat source $var element = document.getElementById("moose-equation-43a80d64-56a8-498f-8d95-741791439431");katex.render("q", element, {displayMode:false,throwOnError:false});$ that will depend on the axial position in the fuel rod and the local Zr fraction (actually the local fraction of actinides, which is the complement of the Zr fraction). Additionally a thermal conductivity model that accounts for constituent migration as well as porosity effects was used (Galloway model described in UPuZrThermal).

In a single-phase region of the phase diagram (here $var element = document.getElementById("moose-equation-f209f200-fea4-41a4-aaa5-0a5c2fee8578");katex.render("\\pi", element, {displayMode:false,throwOnError:false});$ denotes the arbitrary phase) we have $(2)var element = document.getElementById("moose-equation-78418f35-aa0c-4752-9e55-d8ccc1fdd3e5");katex.render(" \\begin{gathered} D(x,T)={{D}_{\\pi}}(x,T)\\\\ S(x,T)={{D}_{\\pi}}(x,T)\\left( \\frac{x(1-{{x}_{\\text{Pu}}}-x)}{1-{{x}_{\\text{Pu}}}} \\right)\\frac{Q_{\\pi}^{\\star}}{R{{T}^{2}}} \\end{gathered}", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-0ca9a807-e749-4de2-be3e-70bdaa509738");katex.render("{{D}_{\\pi }}", element, {displayMode:false,throwOnError:false});$ is the diffusivity and $var element = document.getElementById("moose-equation-b5cf1ab2-3f74-4818-9135-730573424a86");katex.render("Q_{\\pi }^{\\star }", element, {displayMode:false,throwOnError:false});$ the heat of transport of Zr in the phase $var element = document.getElementById("moose-equation-040aa488-9573-4ee4-8a71-dc117e7954e9");katex.render("\\pi", element, {displayMode:false,throwOnError:false});$ . In the 2-phase $var element = document.getElementById("moose-equation-1c679e6a-1965-46f2-b534-86e4e38ff7c6");katex.render("\\beta + \\gamma", element, {displayMode:false,throwOnError:false});$ region bounded by the solubility limit curves $var element = document.getElementById("moose-equation-eee09b61-eea3-4bc2-bc2e-ee5c2ce83cf8");katex.render("{{X}_{\\beta }}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-e47b0044-bd05-49a4-ac5d-995f33f255e0");katex.render("{{X}_{\\gamma }}", element, {displayMode:false,throwOnError:false});$ , the Soret term coefficient is $(3)var element = document.getElementById("moose-equation-80edcdbd-3215-427a-a6e4-8588f6212c97");katex.render(" \\begin{aligned} S(x,T) & = {{f}_{\\beta}}{{D}_{\\beta }}({{X}_{\\beta}}(T),T)\\left( \\frac{{{X}_{\\beta}}(T)(1-{{x}_{\\text{Pu}}}-{{X}_{\\beta}}(T))}{1-{{x}_{\\text{Pu}}}} \\right)\\left( \\frac{\\Delta {{H}_{\\beta}}+Q_{\\beta}^{*}}{R{{T}^{2}}} \\right) \\\\ & +(1-{{f}_{\\beta }}){{D}_{\\gamma}}({{X}_{\\gamma}}(T),T)\\left( \\frac{{{X}_{\\gamma }}(T)(1-{{x}_{\\text{Pu}}}-{{X}_{\\gamma }}(T))}{1-{{x}_{\\text{Pu}}}} \\right)\\left( \\frac{\\Delta {{H}_{\\gamma }}+Q_{\\gamma}^{*}}{R{{T}^{2}}} \\right) \\end{aligned}", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-078e7414-dc79-4cc3-8a9f-6e71a7e23b57");katex.render("{{f}_{\\beta }}=({{X}_{\\gamma }}(T)-x)/({{X}_{\\gamma }}(T)-{{X}_{\\beta }}(T))", element, {displayMode:false,throwOnError:false});$ is the phase fraction of $var element = document.getElementById("moose-equation-73a0da39-e530-4da5-a4e6-136529ae8da3");katex.render("\\beta", element, {displayMode:false,throwOnError:false});$ according to the lever rule. In a 2-phase region one ought to take $var element = document.getElementById("moose-equation-2c33f000-a575-4194-b27f-5214648e9041");katex.render("D=0", element, {displayMode:false,throwOnError:false});$ , however this would result in Eq. (1) being purely advective with a transport velocity proportional to the temperature gradient. Mathematically this leads to jump discontinuities in the Zr atom fraction profile at domain boundaries and boundaries between single and 2-phase regions. Moreover, it is well-known that standard centered finite difference or Galerkin finite element schemes are unstable for pure advection, leading to spurious oscillations. Thus to stabilize our Galerkin finite element implementation we add some artificial diffusion and take, $(4)var element = document.getElementById("moose-equation-b1e75de3-9ec5-416f-9d03-d7a858364ec3");katex.render(" D(x,T)={{f}_{\\beta }}{{D}_{\\beta }}({{X}_{\\beta }}(T),T){{p}_{\\beta }}+(1-{{f}_{\\beta }}){{D}_{\\gamma }}({{X}_{\\gamma }}(T),T){{p}_{\\gamma }},", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-50736bcc-861c-4936-952a-4c3e5b75d247");katex.render("{{p}_{\\beta }}", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-bba81f42-2e28-4246-89ac-86a52369d17a");katex.render("{{p}_{\\gamma }}", element, {displayMode:false,throwOnError:false});$ are dimensionless numerical parameters taken as small as possible while maintaining stability. In future work we plan to replace this simple stabilization approach with the SUPG (streamline upwind Petrov-Galerkin) finite element method, which can be also implemented within BISON. The $var element = document.getElementById("moose-equation-9d01cef7-9f21-496e-b03f-d22d966c5bb8");katex.render("D", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-1ce0a141-5b3a-4cf6-b03e-3f3cc05199db");katex.render("S", element, {displayMode:false,throwOnError:false});$ coefficients in the 2-phase $var element = document.getElementById("moose-equation-1e5b505c-97bf-47f8-95c7-5d08032a58da");katex.render("\\alpha + \\delta", element, {displayMode:false,throwOnError:false});$ region are defined in an analogous way.

Smoothing $var element = document.getElementById("moose-equation-723987c5-b6e2-47f6-8ffc-6c8c32f54bc2");katex.render("D", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-6478affd-1208-4d8a-bcab-990c7d806732");katex.render("S", element, {displayMode:false,throwOnError:false});$

As defined above, $var element = document.getElementById("moose-equation-ce3db8c1-7e12-4085-8ec7-f5491b39ec9c");katex.render("D", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-2758d6d6-5539-496c-9cc8-3f8949272c9b");katex.render("S", element, {displayMode:false,throwOnError:false});$ do not vary smoothly (or even continuously) with $var element = document.getElementById("moose-equation-5e3891d4-bd1b-4460-a5e0-ccf4fdfa73da");katex.render("x", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-2a93ed68-6805-4a9b-bee2-c7e6957fd84e");katex.render("T", element, {displayMode:false,throwOnError:false});$ when crossing phase diagram boundaries, and consequently the discretized system will be much more difficult to solve than it would otherwise be. We have ameliorated this effect by smoothing the relative phase contributions at each location in the phase diagram, discussed in detail in the documentation for PhaseUPuZr.

Example Input Syntax

[Kernels<<<{"href": "../../syntax/Kernels/index.html"}>>>]
  [ZrDiffusion]
    type = ZirconiumDiffusion<<<{"description": "Computes the amount of zirconium that is transported across the mesh.", "href": "ZirconiumDiffusion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = zirconium
    temperature<<<{"description": "Coupled temperature"}>>> = T
  []
[]

Input Parameters

temperatureCoupled temperature
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:Coupled temperature
variableThe name of the variable that this residual object operates on
C++ Type:NonlinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The name of the variable that this residual object operates on

Required Parameters

Zr_SoretZr_SoretZr_Soret
Default:Zr_Soret
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Zr_Soret
Zr_Soret_dTZr_Soret_dTMaterial property name for zirconium Soret diffusivity derivative w.r.t. temperature, used for non-AD applications.
Default:Zr_Soret_dT
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Material property name for zirconium Soret diffusivity derivative w.r.t. temperature, used for non-AD applications.
Zr_Soret_dZrZr_Soret_dZrMaterial property name for zirconium Soret diffusivity derivative w.r.t. zirconium, used for non-AD applications.
Default:Zr_Soret_dZr
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Material property name for zirconium Soret diffusivity derivative w.r.t. zirconium, used for non-AD applications.
Zr_diffusivityZr_diffusivityZr_diffusivity
Default:Zr_diffusivity
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Zr_diffusivity
Zr_diffusivity_dTZr_diffusivity_dTMaterial property name for zirconium diffusivity derivative w.r.t. temperature, used for non-AD applications.
Default:Zr_diffusivity_dT
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Material property name for zirconium diffusivity derivative w.r.t. temperature, used for non-AD applications.
Zr_diffusivity_dZrZr_diffusivity_dZrMaterial property name for zirconium diffusivity derivative w.r.t. zirconium, used for non-AD applications.
Default:Zr_diffusivity_dZr
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:Material property name for zirconium diffusivity derivative w.r.t. zirconium, used for non-AD applications.
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
displacementsThe displacements
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The displacements

Optional Parameters

absolute_value_vector_tagsThe tags for the vectors this residual object should fill with the absolute value of the residual contribution
C++ Type:std::vector<TagName>
Controllable:No
Description:The tags for the vectors this residual object should fill with the absolute value of the residual contribution
extra_matrix_tagsThe extra tags for the matrices this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the matrices this Kernel should fill
extra_vector_tagsThe extra tags for the vectors this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the vectors this Kernel should fill
matrix_tagssystemThe tag for the matrices this Kernel should fill
Default:system
C++ Type:MultiMooseEnum
Options:nontime, system
Controllable:No
Description:The tag for the matrices this Kernel should fill
vector_tagsnontimeThe tag for the vectors this Kernel should fill
Default:nontime
C++ Type:MultiMooseEnum
Options:nontime, time
Controllable:No
Description:The tag for the vectors this Kernel should fill

Contribution To Tagged Field Data 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.
diag_save_inThe name of auxiliary variables to save this Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The name of auxiliary variables to save this Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
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
save_inThe name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
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

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/zirconium_diffusion/ad_exact.i)
(examples/metal_fuel/uzr/pin.i)
(test/tests/zirconium_diffusion/exact.i)
(test/tests/zirconium_diffusion/ad_test.i)
(test/tests/zirconium_diffusion/test.i)

References

Yeon Soo Kim, S L Hayes, G L Hofman, and A M Yacout. Modeling of constituent redistribution in U–Pu–Zr metallic fuel. Journal of Nuclear Materials, 359(1-2):17–28, December 2006.

@article{Kim:2006hj,
    author = "Kim, Yeon Soo and Hayes, S L and Hofman, G L and Yacout, A M",
    title = "{Modeling of constituent redistribution in U{\textendash}Pu{\textendash}Zr metallic fuel}",
    journal = "Journal of Nuclear Materials",
    year = "2006",
    volume = "359",
    number = "1-2",
    pages = "17-28",
    month = "December"
}

Yeon Soo Kim, G L Hofman, S L Hayes, and Y H Sohn. Constituent redistribution in U–Pu–Zr fuel during irradiation. Journal of Nuclear Materials, 327(1):27–36, April 2004.

@article{Kim:2004ei,
    author = "Kim, Yeon Soo and Hofman, G L and Hayes, S L and Sohn, Y H",
    title = "{Constituent redistribution in U{\textendash}Pu{\textendash}Zr fuel during irradiation}",
    journal = "Journal of Nuclear Materials",
    year = "2004",
    volume = "327",
    number = "1",
    pages = "27-36",
    month = "April"
}

(test/tests/zirconium_diffusion/test.i)

# This is an example of an axial slice from the T179 EBR-II X419 experiment,
# which had 61U-16Pu-23Zr atom%, 71U-19Pu-10Zr wt%.
#
# Most importantly, this test shows that the new calculation method for zirconium
# redistribution matches the old method, to a maximum difference of 3e-3.

[Mesh]
  coord_type = RZ
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 104
    ny = 1
    xmax = 0.00216
    ymax = 2.16e-05
  []
[]

[Variables]
  [T]
  []
  [zirconium]
    initial_condition = 0.22 #initial mole fraction of Zr in the system
    scaling = 1e17
  []
[]

[ICs]
  [T]
    type = FunctionIC
    variable = T
    function = outer2
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = T
  []
  [ie]
    type = HeatConductionTimeDerivative
    variable = T
  []
  [ZrDiffusion]
    type = ZirconiumDiffusion
    variable = zirconium
    temperature = T
  []
  [diffusion_ie]
    type = TimeDerivative
    variable = zirconium
  []
  [heat_source]
    type = CoupledForce
    variable = T
    v = power
  []
[]

[AuxVariables]
  [fission_rate]
  []
  [power]
  []
[]

[AuxKernels]
  [fission_rate_aux]
    type = ConstantAux
    variable = fission_rate
    value = 8.74117559524e+19
  []
  [power_aux]
    type = ParsedAux
    coupled_variables = 'zirconium fission_rate'
    constant_names = 'A B epf'
    constant_expressions = '-1.12884007 1.24834482 3.2e-11'
    expression = 'fission_rate * epf * (A*zirconium + B)'
    variable = power
  []
[]

[Functions]
  [outer2]
    type = ParsedFunction
    expression = '918.15-3.2e7*(x*x)' #r^2 initial temperature profile
  []
[]

[BCs]
  [temp_rightside]
    type = DirichletBC
    variable = T
    boundary = right
    value = 818.15
  []
[]

[Materials]
  [porosity]
    type = GenericConstantMaterial
    block = 0
    prop_names = porosity
    prop_values = 0.0
  []
  [phase]
    type = PhaseUPuZr
    block = 0
    temperature = T
    X_Pu = 0.16
    X_Zr = zirconium
    AB_temp = 936.15
    CD_temp = 968.15
  []
  [zr_diff]
    type = ZrDiffusivityUPuZr
    block = 0
    temperature = T
    X_Pu = 0.16
    X_Zr = zirconium
    p_alpha = ${fparse 0.35 / 100}
    p_delta = ${fparse 0.35 / 2}
    p_beta = ${fparse 0.35 / 55}
    p_gamma = ${fparse 0.35 / 2}
    D0_scale_alpha = 100.0
    D0_scale_delta = 2.0
    D0_scale_beta = 55.0
    D0_scale_gamma = 2.0
  []
  [thermal_upuzr]
    type = UPuZrThermal
    block = 0
    temperature = T
    X_Pu = 0.16
    X_Zr = zirconium
    thcond_model = galloway
    spheat_model = savage
    A_Pu = .244
    porosity = porosity
  []
  [fuel_density]
    type = ParsedMaterial
    block = 0
    property_name = density
    expression = 10431.0
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-ksp_gmres_restart -pc_type    -sub_pc_type -pc_hypre_type'
  petsc_options_value = ' 1001               hypre         hypre        boomeramg'
  line_search = 'none'

  l_tol = 8e-3
  nl_max_its = 15
  nl_rel_tol = 1e-9
  nl_abs_tol = 1e-18

  l_max_its = 1000
  dtmax = 1e5
  dtmin = 1e-2
  end_time = 1e5

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1.0e3
    optimal_iterations = 25
    iteration_window = 3
    linear_iteration_ratio = 35
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/zirconium_diffusion/ad_exact.i)

# This input file mirrors the method in moose/modules/misc/tests/kernels/thermo_diffusion.i

# Steady-state test for the ZirconiumDiffusion kernel.
#
# This test applies a constant temperature gradient to drive thermo-diffusion
# in the variable u. At steady state, the thermo-diffusion is balanced by
# diffusion due to Fick's Law, so the total flux is
#
#   J = -D grad(u) - S grad(T) = 0
#
# If there are no fluxes at the boundaries, then there is no background flux and
# these two terms must balance each other everywhere:
#
#   grad(u) = -(S / D) grad(T)
#
# Using S = D * u(1-u) * Qstar/R/T^2, and setting Qstar=R the dx can be eliminated to give
#
#   du = -u(1-u) dT/T^2
#
# This can be solved to give the profile for u as a function of temperature:
#
#   u = exp(T^-1 + C) / ( 1 + exp(T^-1 + C ) )
#
# Here, we are using simple heat conduction with Dirichlet boundaries on 0 <= x <= 1
# to give a linear profile for temperature: T = x + 1:
#
#   u = exp(1/(x+1) + C) / ( 1 + exp(1/(x+1) + C ) )
#
# Using the BC u(x=0) = 0.5, then C = -1.
# This analytical result is tracked by the aux variable "correct_u".
#
#
# Postprocessor Values:
# +----------------+----------------+----------------+----------------+
# | time           | correct_u      | temp           | u              |
# +----------------+----------------+----------------+----------------+
# |   0.000000e+00 |   4.243783e-01 |   1.000000e+00 |   0.000000e+00 |
# |   1.000000e+00 |   4.243783e-01 |   1.500000e+00 |   4.243658e-01 |
# +----------------+----------------+----------------+----------------+

[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = 1
    ymin = 0
    ymax = 1
    nx = 10
    ny = 1
  []
[]

[Variables]
  [u]
    initial_condition = 0.5
  []
  [temp]
    initial_condition = 1
  []
[]

[Kernels]
  [zr_diff]
    type = ADZirconiumDiffusion
    variable = u
    temperature = temp
  []

  # Heat diffusion gives a linear temperature profile to drive the Soret diffusion.
  [diffT]
    type = ADDiffusion
    variable = temp
  []
[]

[BCs]
  [left]
    type = ADDirichletBC
    variable = u
    boundary = left
    value = 0.5
  []

  [leftt]
    type = ADDirichletBC
    variable = temp
    boundary = left
    value = 1
  []
  [rightt]
    type = ADDirichletBC
    variable = temp
    boundary = right
    value = 2
  []
[]

[Materials]
  [zr_diff]
    type = ADZrDiffusivityUPuZr
    block = 0
    temperature = temp
    X_Pu = 0
    X_Zr = u
    Qstar_alpha = 8.31446
    Q_scale_alpha = 0
  []
  [phase]
    type = ADPhaseUPuZr
    block = 0
    temperature = 800
    X_Pu = 0
    X_Zr = 0.01
  []
[]

[Functions]
  # The correct (analytical) result for the solution of u
  [concentration_profile]
    type = ParsedFunction
    expression = 'exp(1.0/(x+1.0)-1) / ( 1.0 + exp(1.0/(x+1.0)-1) )'
  []
[]

[AuxVariables]
  [correct_u]
  []
[]

[AuxKernels]
  # This kernel is just being used to plot the correct result for u.
  [copy_from_function]
    type = FunctionAux
    variable = correct_u
    function = concentration_profile
    execute_on = 'initial timestep_begin'
  []
[]

[Executioner]
  type = Steady

  nl_abs_tol = 1e-50
  nl_rel_tol = 1e-12
  l_tol = 1e-10
[]

[Postprocessors]
  [u]
    type = ElementAverageValue
    variable = u
    execute_on = 'initial timestep_end'
  []
  [temp]
    type = ElementAverageValue
    variable = temp
    execute_on = 'initial timestep_end'
  []
  [correct_u]
    type = ElementAverageValue
    variable = correct_u
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  exodus = true
[]

(examples/metal_fuel/uzr/pin.i)

initial_fuel_density = 15800

[GlobalParams]
  energy_per_fission = 3.2e-11  # J/fission
  temperature = temp
[]

[Mesh]
  coord_type = RZ
  [mesh]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 300
    xmax = 0.002
  []
[]

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

[Variables]
  [temp]
    initial_condition = 298
  []
  [X_Zr]
    scaling = 1e10
    initial_condition = 0.225
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 100000 3901600 4333600 6148000 7948800  '
    y = '0 40016 40344 42312 41984 41984  '
  []
  [temp_history]
    type = PiecewiseLinear
    x = '0 100000 3901600 4333600 6148000 7948800  '
    y = '298.0 805.5 803.8 811.1 808.5 807.4'
  []
[]

[Kernels]
  [heat]
    type = ConstitutiveHeatConduction
    variable = temp
    thermal_conductivity = 'thermal_conductivity'
    thermal_conductivity_args = 'temp X_Zr'
    thermal_conductivity_derivs = 'thermal_conductivity_dT thermal_conductivity_dZr'
    extra_vector_tags = 'ref'
  []
  [heat_ie]
    type = ConstitutiveHeatConductionTimeDerivative
    variable = temp
    specific_heat = 'specific_heat'
    specific_heat_args = 'temp'
    specific_heat_derivs = 'specific_heat_dT'
    extra_vector_tags = 'ref'
  []
  [heat_source]
    type = FissionRateHeatSource
    variable = temp
    fission_rate_args = 'X_Zr'
    fission_rate_derivs = 'fission_rate_dZr'
    extra_vector_tags = 'ref'
  []
  [ZrDiffusion]
    type = ZirconiumDiffusion
    variable = X_Zr
    temperature = temp
    extra_vector_tags = 'ref'
  []
  [diffusion_ie]
    type = TimeDerivative
    variable = X_Zr
    extra_vector_tags = 'ref'
  []
[]

[BCs]
  [temp_rightside]
    type = FunctionDirichletBC
    variable = temp
    boundary = right
    function = temp_history
  []
[]

[Materials]
  [fission_rate]
    type = UPuZrFissionRate
    rod_linear_power = power_history
    axial_power_profile = 1.0
    pellet_radius = 0.002
    initial_X_Zr = 0.225
    X_Zr = X_Zr
    X_Pu_function = 0.163
  []
  [burnup]
    type = UPuZrBurnup
    initial_X_Zr = 0.225
    initial_X_Pu = 0.163
    density = ${initial_fuel_density}
  []
  [metal_fuel_thermal]
    type = UPuZrThermal
    spheat_model = savage
    thcond_model = lanl
    porosity = porosity
    X_Zr = X_Zr
    X_Pu = 0.163
  []
  [swelling]
    type = UPuZrVolumetricSwellingEigenstrain
    burnup = 0
    hydrostatic_stress = 0
    eigenstrain_name = 'swelling'
    output_properties = 'porosity gaseous_porosity'
  []
  [fuel_density]
    type = GenericConstantMaterial
    prop_names = density
    prop_values = ${initial_fuel_density}
  []
  [phase]
    type = PhaseUPuZr
    X_Zr = X_Zr
    X_Pu = 0.163
    AB_temp = 965.15
    CD_temp = 995.15
  []
  [zr_diff]
    type = ZrDiffusivityUPuZr
    X_Zr = X_Zr
    X_Pu = 0.163
    p_alpha = 0.20
    p_delta = 0.20
    p_beta  = 0.20
    p_gamma = 0.20
    D0_scale_alpha = 15.0
    D0_scale_delta = 1.0
    D0_scale_beta = 2.0
    D0_scale_gamma = 7.0
    outputs = all
  []
[]

[Dampers]
  [bounded_zr]
    type = BoundingValueNodalDamper
    variable = X_Zr
    max_value = 0.83
    min_value = 0
  []
[]

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

[Executioner]
  type = Transient

  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'
  line_search = 'none'

  end_time = 3e6
  dtmin = 1e-2
  dtmax = 1e6
  nl_max_its = 30
  l_max_its = 100
  nl_abs_tol = 1e-9
  nl_rel_tol = 1e-9
  l_tol = 1e-05

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e3
    optimal_iterations = 8
    iteration_window = 2
    growth_factor = 2
    cutback_factor = .5
    timestep_limiting_function = power_history
    force_step_every_function_point = true
  []
[]

[Postprocessors]
  [dt]
    type = TimestepSize
  []
  [num_lin]
    type = NumLinearIterations
  []
  [num_nonlin]
    type = NumNonlinearIterations
  []
  [ave_temp]
    type = ElementAverageValue
    variable = temp
    block = 0
  []
  [X_Zr]
    type = NodalExtremeValue
    variable = X_Zr
    block = 0
    value_type = max
  []
[]

[Outputs]
  perf_graph = true
  #exodus = true

  [csv_output]
    type = CSV
    precision = 8
    file_base = pin_out
    hide = 'dt num_lin num_nonlin'
    execute_on = final
  []
[]

(test/tests/zirconium_diffusion/exact.i)

# This input file mirrors the method in moose/modules/misc/tests/kernels/thermo_diffusion.i

# Steady-state test for the ZirconiumDiffusion kernel.
#
# This test applies a constant temperature gradient to drive thermo-diffusion
# in the variable u. At steady state, the thermo-diffusion is balanced by
# diffusion due to Fick's Law, so the total flux is
#
#   J = -D grad(u) - S grad(T) = 0
#
# If there are no fluxes at the boundaries, then there is no background flux and
# these two terms must balance each other everywhere:
#
#   grad(u) = -(S / D) grad(T)
#
# Using S = D * u(1-u) * Qstar/R/T^2, and setting Qstar=R the dx can be eliminated to give
#
#   du = -u(1-u) dT/T^2
#
# This can be solved to give the profile for u as a function of temperature:
#
#   u = exp(T^-1 + C) / ( 1 + exp(T^-1 + C ) )
#
# Here, we are using simple heat conduction with Dirichlet boundaries on 0 <= x <= 1
# to give a linear profile for temperature: T = x + 1:
#
#   u = exp(1/(x+1) + C) / ( 1 + exp(1/(x+1) + C ) )
#
# Using the BC u(x=0) = 0.5, then C = -1.
# This analytical result is tracked by the aux variable "correct_u".
#
#
# Postprocessor Values:
# +----------------+----------------+----------------+----------------+
# | time           | correct_u      | temp           | u              |
# +----------------+----------------+----------------+----------------+
# |   0.000000e+00 |   4.243783e-01 |   1.000000e+00 |   0.000000e+00 |
# |   1.000000e+00 |   4.243783e-01 |   1.500000e+00 |   4.243658e-01 |
# +----------------+----------------+----------------+----------------+

[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = 1
    ymin = 0
    ymax = 1
    nx = 10
    ny = 1
  []
[]

[Variables]
  [u]
    initial_condition = 0.5
  []
  [temp]
    initial_condition = 1
  []
[]

[Kernels]
  [zr_diff]
    type = ZirconiumDiffusion
    variable = u
    temperature = temp
  []

  # Heat diffusion gives a linear temperature profile to drive the Soret diffusion.
  [diffT]
    type = Diffusion
    variable = temp
  []
[]

[BCs]
  [left]
    type = DirichletBC
    variable = u
    boundary = left
    value = 0.5
  []

  [leftt]
    type = DirichletBC
    variable = temp
    boundary = left
    value = 1
  []
  [rightt]
    type = DirichletBC
    variable = temp
    boundary = right
    value = 2
  []
[]

[Materials]
  [zr_diff]
    type = ZrDiffusivityUPuZr
    block = 0
    temperature = temp
    X_Pu = 0
    X_Zr = u
    Qstar_alpha = 8.31446
    Q_scale_alpha = 0
  []
  [phase]
    type = PhaseUPuZr
    block = 0
    temperature = 800
    X_Pu = 0
    X_Zr = 0.01
  []
[]

[Functions]
  # The correct (analytical) result for the solution of u
  [concentration_profile]
    type = ParsedFunction
    expression = 'exp(1.0/(x+1.0)-1) / ( 1.0 + exp(1.0/(x+1.0)-1) )'
  []
[]

[AuxVariables]
  [correct_u]
  []
[]

[AuxKernels]
  # This kernel is just being used to plot the correct result for u.
  [copy_from_function]
    type = FunctionAux
    variable = correct_u
    function = concentration_profile
    execute_on = 'initial timestep_begin'
  []
[]

[Executioner]
  type = Steady

  nl_abs_tol = 1e-50
  nl_rel_tol = 1e-12
  l_tol = 1e-10
[]

[Postprocessors]
  [u]
    type = ElementAverageValue
    variable = u
    execute_on = 'initial timestep_end'
  []
  [temp]
    type = ElementAverageValue
    variable = temp
    execute_on = 'initial timestep_end'
  []
  [correct_u]
    type = ElementAverageValue
    variable = correct_u
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/zirconium_diffusion/ad_test.i)

# This is an example of an axial slice from the T179 EBR-II X419 experiment,
# which had 61U-16Pu-23Zr atom%, 71U-19Pu-10Zr wt%.
#
# Most importantly, this test shows that the new calculation method for zirconium
# redistribution matches the old method, to a maximum difference of 3e-3.

[Mesh]
  coord_type = RZ
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 104
    ny = 1
    xmax = 0.00216
    ymax = 2.16e-05
  []
[]

[Variables]
  [T]
  []
  [zirconium]
    initial_condition = 0.22 #initial mole fraction of Zr in the system
    scaling = 1e17
  []
[]

[ICs]
  [T]
    type = FunctionIC
    variable = T
    function = outer2
  []
[]

[Kernels]
  [heat]
    type = ADHeatConduction
    variable = T
  []
  [ie]
    type = ADHeatConductionTimeDerivative
    variable = T
  []
  [ZrDiffusion]
    type = ADZirconiumDiffusion
    variable = zirconium
    temperature = T
  []
  [diffusion_ie]
    type = ADTimeDerivative
    variable = zirconium
  []
  [heat_source]
    type = ADCoupledForce
    variable = T
    v = power
  []
[]

[AuxVariables]
  [fission_rate]
  []
  [power]
  []
[]

[AuxKernels]
  [fission_rate_aux]
    type = ConstantAux
    variable = fission_rate
    value = 8.74117559524e+19
  []
  [power_aux]
    type = ParsedAux
    coupled_variables = 'zirconium fission_rate'
    constant_names = 'A B epf'
    constant_expressions = '-1.12884007 1.24834482 3.2e-11'
    expression = 'fission_rate * epf * (A*zirconium + B)'
    variable = power
  []
[]

[Functions]
  [outer2]
    type = ParsedFunction
    expression = '918.15-3.2e7*(x*x)' #r^2 initial temperature profile
  []
[]

[BCs]
  [temp_rightside]
    type = ADDirichletBC
    variable = T
    boundary = right
    value = 818.15
  []
[]

[Materials]
  [porosity]
    type = ADGenericConstantMaterial
    block = 0
    prop_names = porosity
    prop_values = 0.0
  []
  [phase]
    type = ADPhaseUPuZr
    block = 0
    temperature = T
    X_Pu = 0.16
    X_Zr = zirconium
    AB_temp = 936.15
    CD_temp = 968.15
  []
  [zr_diff]
    type = ADZrDiffusivityUPuZr
    block = 0
    temperature = T
    X_Pu = 0.16
    X_Zr = zirconium
    p_alpha = ${fparse 0.35 / 100}
    p_delta = ${fparse 0.35 / 2}
    p_beta = ${fparse 0.35 / 55}
    p_gamma = ${fparse 0.35 / 2}
    D0_scale_alpha = 100.0
    D0_scale_delta = 2.0
    D0_scale_beta = 55.0
    D0_scale_gamma = 2.0
  []
  [thermal_upuzr]
    type = ADUPuZrThermal
    block = 0
    temperature = T
    X_Pu = 0.16
    X_Zr = zirconium
    thcond_model = galloway
    spheat_model = savage
    A_Pu = .244
    porosity = porosity
  []
  [fuel_density]
    type = ADParsedMaterial
    block = 0
    property_name = density
    expression = 10431.0
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-ksp_gmres_restart -pc_type    -sub_pc_type -pc_hypre_type'
  petsc_options_value = ' 1001               hypre         hypre        boomeramg'
  line_search = 'none'

  l_tol = 8e-3
  nl_max_its = 15
  nl_rel_tol = 1e-9
  nl_abs_tol = 1e-18

  l_max_its = 1000
  dtmax = 1e5
  dtmin = 1e-2
  end_time = 1e5

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1.0e3
    optimal_iterations = 25
    iteration_window = 3
    linear_iteration_ratio = 35
  []
[]

[Outputs]
  exodus = true
[]

(test/tests/zirconium_diffusion/test.i)

# This is an example of an axial slice from the T179 EBR-II X419 experiment,
# which had 61U-16Pu-23Zr atom%, 71U-19Pu-10Zr wt%.
#
# Most importantly, this test shows that the new calculation method for zirconium
# redistribution matches the old method, to a maximum difference of 3e-3.

[Mesh]
  coord_type = RZ
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 104
    ny = 1
    xmax = 0.00216
    ymax = 2.16e-05
  []
[]

[Variables]
  [T]
  []
  [zirconium]
    initial_condition = 0.22 #initial mole fraction of Zr in the system
    scaling = 1e17
  []
[]

[ICs]
  [T]
    type = FunctionIC
    variable = T
    function = outer2
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = T
  []
  [ie]
    type = HeatConductionTimeDerivative
    variable = T
  []
  [ZrDiffusion]
    type = ZirconiumDiffusion
    variable = zirconium
    temperature = T
  []
  [diffusion_ie]
    type = TimeDerivative
    variable = zirconium
  []
  [heat_source]
    type = CoupledForce
    variable = T
    v = power
  []
[]

[AuxVariables]
  [fission_rate]
  []
  [power]
  []
[]

[AuxKernels]
  [fission_rate_aux]
    type = ConstantAux
    variable = fission_rate
    value = 8.74117559524e+19
  []
  [power_aux]
    type = ParsedAux
    coupled_variables = 'zirconium fission_rate'
    constant_names = 'A B epf'
    constant_expressions = '-1.12884007 1.24834482 3.2e-11'
    expression = 'fission_rate * epf * (A*zirconium + B)'
    variable = power
  []
[]

[Functions]
  [outer2]
    type = ParsedFunction
    expression = '918.15-3.2e7*(x*x)' #r^2 initial temperature profile
  []
[]

[BCs]
  [temp_rightside]
    type = DirichletBC
    variable = T
    boundary = right
    value = 818.15
  []
[]

[Materials]
  [porosity]
    type = GenericConstantMaterial
    block = 0
    prop_names = porosity
    prop_values = 0.0
  []
  [phase]
    type = PhaseUPuZr
    block = 0
    temperature = T
    X_Pu = 0.16
    X_Zr = zirconium
    AB_temp = 936.15
    CD_temp = 968.15
  []
  [zr_diff]
    type = ZrDiffusivityUPuZr
    block = 0
    temperature = T
    X_Pu = 0.16
    X_Zr = zirconium
    p_alpha = ${fparse 0.35 / 100}
    p_delta = ${fparse 0.35 / 2}
    p_beta = ${fparse 0.35 / 55}
    p_gamma = ${fparse 0.35 / 2}
    D0_scale_alpha = 100.0
    D0_scale_delta = 2.0
    D0_scale_beta = 55.0
    D0_scale_gamma = 2.0
  []
  [thermal_upuzr]
    type = UPuZrThermal
    block = 0
    temperature = T
    X_Pu = 0.16
    X_Zr = zirconium
    thcond_model = galloway
    spheat_model = savage
    A_Pu = .244
    porosity = porosity
  []
  [fuel_density]
    type = ParsedMaterial
    block = 0
    property_name = density
    expression = 10431.0
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options = '-snes_ksp_ew'
  petsc_options_iname = '-ksp_gmres_restart -pc_type    -sub_pc_type -pc_hypre_type'
  petsc_options_value = ' 1001               hypre         hypre        boomeramg'
  line_search = 'none'

  l_tol = 8e-3
  nl_max_its = 15
  nl_rel_tol = 1e-9
  nl_abs_tol = 1e-18

  l_max_its = 1000
  dtmax = 1e5
  dtmin = 1e-2
  end_time = 1e5

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1.0e3
    optimal_iterations = 25
    iteration_window = 3
    linear_iteration_ratio = 35
  []
[]

[Outputs]
  exodus = true
[]

Description
Example Input Syntax
Input Parameters
Input Files
References