REBEKADirichletBC

Creates a temperature boundary condition specific to the REBEKA LOCA simulations.

Description

The REBEKADirichletBC boundary condition is specific for application to the REBEKA LOCA simulations. It provides a Dirichlet temperature boundary condition with a parabolic profile along the axial direction that is based on the conditions of the REBEKA experimental tests (Erbacher et al., 1982; Markiewicz and Erbacher, 1988; Erbacher et al., 1990). The temperature profile is based on calculations of axial heat conduction within the fuel rod simulator and convection to the outer steam atmosphere, based on available information on materials and experimental conditions. In addition, an azimuthal temperature variation can be specified if the boundary condition is applied in 3D simulation.

Example Input Syntax

For 2D axisymmetric simulations:

[BCs<<<{"href": "../../syntax/NuclearMaterials/BCs/index.html"}>>>]
  [inner_temperature]
    type = REBEKADirichletBC<<<{"description": "Creates a temperature boundary condition specific to the REBEKA LOCA simulations.", "href": "REBEKADirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
    function_tempm<<<{"description": "Function for the peak temperature vs time."}>>> = temperature_func
    minimum_temperature<<<{"description": "Minimum temperature considered (K)."}>>> = 573.
    translation<<<{"description": "Translation along the y-axis (m)."}>>> = 0.1625
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 4
  []
[]

For 3D with an azimuthal variation:

[BCs<<<{"href": "../../syntax/NuclearMaterials/BCs/index.html"}>>>]
  [inner_temperature]
    type = REBEKADirichletBC<<<{"description": "Creates a temperature boundary condition specific to the REBEKA LOCA simulations.", "href": "REBEKADirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
    function_tempm<<<{"description": "Function for the peak temperature vs time."}>>> = temperature_func
    translation<<<{"description": "Translation along the y-axis (m)."}>>> = 0.1625
    vary_azimuthally<<<{"description": "For 3D problems if the bc is to vary azimuthally"}>>> = true
    azimuthal_variation_pi2<<<{"description": "Azimuthal variation factor across a pi/2 angle."}>>> = -30.
    minimum_temperature<<<{"description": "Minimum temperature considered (K)."}>>> = 573.0
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 4
  []
[]

Input Parameters

boundaryThe list of boundary IDs from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundary IDs from the mesh where this object applies
function_tempmFunction for the peak temperature vs time.
C++ Type:FunctionName
Unit:(no unit assumed)
Controllable:No
Description:Function for the peak temperature vs time.
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

azimuthal_variation_pi20Azimuthal variation factor across a pi/2 angle.
Default:0
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Azimuthal variation factor across a pi/2 angle.
diag_save_inThe name of auxiliary variables to save this BC'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 BC's diagonal jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
displacementsThe displacements
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The displacements
minimum_temperature0Minimum temperature considered (K).
Default:0
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Minimum temperature considered (K).
save_inThe name of auxiliary variables to save this BC'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 BC's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
translation0Translation along the y-axis (m).
Default:0
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Translation along the y-axis (m).
vary_azimuthallyFalseFor 3D problems if the bc is to vary azimuthally
Default:False
C++ Type:bool
Controllable:No
Description:For 3D problems if the bc is to vary azimuthally

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

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.
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

matrix_tagssystem timeThe tag for the matrices this Kernel should fill
Default:system time
C++ Type:MultiMooseEnum
Options:nontime, system, time
Controllable:No
Description:The tag for the matrices this Kernel should fill
vector_tagsresidualThe tag for the vectors this Kernel should fill
Default:residual
C++ Type:MultiMooseEnum
Options:nontime, time, residual
Controllable:No
Description:The tag for the vectors this Kernel should fill

Tagging Parameters

Input Files

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_08MPa/rebeka_singlerod_2d_08MPa.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_14MPa/rebeka_singlerod_2d_14MPa.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_10MPa/rebeka_singlerod_2d_10MPa.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_10MPa/rebeka_singlerod_2d_10MPa_aniso_hunt.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_06MPa/rebeka_singlerod_2d_06MPa.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_12MPa/rebeka_singlerod_2d_12MPa_aniso.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_10MPa/rebeka_singlerod_2d_10MPa_1pt5_aniso.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_04MPa/rebeka_singlerod_2d_04MPa.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_04MPa/rebeka_singlerod_2d_04MPa_aniso.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_06MPa/rebeka_singlerod_2d_06MPa_aniso.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_01MPa/rebeka_singlerod_2d_01MPa.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_02MPa/rebeka_singlerod_2d_02MPa_aniso.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_12MPa/rebeka_singlerod_2d_12MPa.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_10MPa/rebeka_singlerod_2d_10MPa_aniso.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_14MPa/rebeka_singlerod_2d_14MPa_aniso.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_3d_10MPa_01Kpers_90deg/rebeka_singlerod_3d_10MPa_90deg.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_10MPa/rebeka_singlerod_2d_10MPa_1pt5.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_08MPa/rebeka_singlerod_2d_08MPa_aniso.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_02MPa/rebeka_singlerod_2d_02MPa.i)
(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_01MPa/rebeka_singlerod_2d_01MPa_aniso.i)

References

F. J. Erbacher, H. J. Neitzel, H. Rosinger, H. Schmidt, and K. Wiehr. Burst criterion of Zircaloy fuel claddings in a loss-of-coolant accident. In Zirconium in the Nuclear Industry, Fifth Conference, ASTM STP 754, D.G. Franklin Ed., 271–283. American Society for Testing and Materials, 1982.

@INPROCEEDINGS{erbacher_ea_1982,
    author = "Erbacher, F. J. and Neitzel, H. J. and Rosinger, H. and Schmidt, H. and Wiehr, K.",
    title = "Burst criterion of {Z}ircaloy fuel claddings in a loss-of-coolant accident",
    booktitle = "Zirconium in the Nuclear Industry, Fifth Conference, ASTM STP 754, D.G. Franklin Ed.",
    organization = "American Society for Testing and Materials",
    year = "1982",
    pages = "271-283"
}

F.J. Erbacher, H.J. Neitzel, and K. Wiehr. Technical Report KfK 4781, Kernforschungszentrum, Karlsruhe, 1990.

@techreport{erbacher_et_al_1990,
    author = "Erbacher, F.J. and Neitzel, H.J. and Wiehr, K.",
    title = "",
    year = "1990",
    number = "{KfK} 4781",
    institution = "Kernforschungszentrum, Karlsruhe"
}

M. E. Markiewicz and F.J. Erbacher. Experiments on ballooning in pressurized and transiently heated Zircaloy-4 tubes. Technical Report KfK 4343, Kernforschungszentrum Karlsruhe GmbH (Germany, Kernforschungszentrum Karlsruhe, Germany, 1988.

@TECHREPORT{markiewicz_erbacher_1988,
    author = "Markiewicz, M. E. and Erbacher, F.J.",
    title = "Experiments on ballooning in pressurized and transiently heated {Z}ircaloy-4 tubes",
    year = "1988",
    number = "KfK 4343",
    address = "Kernforschungszentrum Karlsruhe, Germany",
    institution = "Kernforschungszentrum Karlsruhe GmbH (Germany"
}

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_10MPa/rebeka_singlerod_2d_10MPa.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

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

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

[Variables]
  [temperature]
    initial_condition = 573.0
  []
[]

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+07  1.e+07' # 100 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = true
    strain = FINITE
    decomposition_method = EigenSolution
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = MaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = MaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = MaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = MaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.
    translation = 0.1625
    boundary = 4
  []
  [no_y_midsection]
    type = DirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[Constraints]
  [bottom_section_plane]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    primary = 101 # node on boundary
    secondary = 1 # boundary
    penalty = 1.e+10
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    oxide_thickness = scale_thickness
  []
[]

[Materials]
  [thermal]
    type = ZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ZryElasticityTensor
    block = cladding
    matpro_youngs_modulus = false
    matpro_poissons_ratio = false
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [clad_zrycreep]
    type = ZryCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501
    temperature_standard_thermal_creep_end = 500
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = StrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[Executioner]
  type = Transient
  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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-06
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp]
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [min_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = min
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = cladding
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = cladding
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geo
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
  [top_disp_r_clad] # this is mid height
    type = NodalVariableValue
    variable = disp_x
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [top_disp_r_clad_slice] # this is mid height matched to the 1.5d
    type = NodalVariableValue
    variable = disp_x
    nodeid = 3 #coords (0.005375, 0.1625)
  []
  [top_disp_z_clad]
    type = NodalVariableValue
    variable = disp_y
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [stress_xx] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_xx
    elementid = 0
  []
  [stress_yy] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_yy
    elementid = 0
  []
  [stress_zz] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_zz
    elementid = 0
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_10MPa_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  [fuel_pin_geo]
    type = FuelPinGeometry
    clad_outer_wall = '2'
    clad_inner_wall = '4'
    include_fuel = false
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_3d_10MPa_01Kpers_90deg/rebeka_singlerod_3d_10MPa_90deg.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
  order = SECOND
  family = LAGRANGE
[]

[Problem]
[]

[Mesh]
  [mesh]
    type = FileMeshGenerator
    file = 3d_rebeka_singlerod.e
  []
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
  [temperature]
    initial_condition = 573.0
  []
[]

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+07  1.e+07' # 100 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = false
    strain = FINITE
    decomposition_method = EigenSolution
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase] # Fraction of beta phase in Zry
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness] # ZrO2 scale thickness (m)
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total] # Total oxigen weight fraction (oxide+metal) (/)
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress] # Hoop stress at cladding burst
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = MaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = MaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = MaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = MaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    translation = 0.1625
    vary_azimuthally = true
    azimuthal_variation_pi2 = -30.
    minimum_temperature = 573.0
    boundary = 4
  []
  [no_y_midsection]
    type = DirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [no_z_frontside]
    type = DirichletBC
    variable = disp_z
    boundary = 5
    value = 0.
  []
  [no_x_backside]
    type = DirichletBC
    variable = disp_x
    boundary = 6
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    number_axial_zone = 15
    oxide_thickness = scale_thickness
  []
[]

[Materials]
  [thermal]
    type = ZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = cladding
    youngs_modulus = 1.0e11
    poissons_ratio = 0.3
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [clad_zrycreep]
    type = ZryCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501 # see Erbacher et al., 1982
    temperature_standard_thermal_creep_end = 500
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = StrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[Executioner]
  type = Transient
  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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-04
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp]
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_singlerod_3d_10MPa_90deg_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_08MPa/rebeka_singlerod_2d_08MPa.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

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

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

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

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '8.e+06  8.e+06' # 80 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = false
    strain = FINITE
    decomposition_method = EigenSolution
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    eigenstrain_names = 'clad_thermal_eigenstrain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = MaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = MaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = MaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = MaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.
    translation = 0.1625
    boundary = 4
  []
  [no_y_midsection]
    type = DirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[Constraints]
  [bottom_section_plane]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    primary = 101 # node on boundary
    secondary = 1 # boundary
    penalty = 1.e+10
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    number_axial_zone = 15
    oxide_thickness = scale_thickness
  []
[]

[Materials]
  [thermal]
    type = ZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = cladding
    youngs_modulus = 1.0e11
    poissons_ratio = 0.3
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [clad_zrycreep]
    type = ZryCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501 # see Erbacher et al., 1982
    temperature_standard_thermal_creep_end = 500
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = StrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[Executioner]
  type = Transient
  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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-04
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp] # average temperature of cladding exterior
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = cladding
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = cladding
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geo
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_08MPa_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  [fuel_pin_geo]
    type = FuelPinGeometry
    clad_outer_wall = '2'
    clad_inner_wall = '4'
    include_fuel = false
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_14MPa/rebeka_singlerod_2d_14MPa.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

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

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

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

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.4e+07  1.4e+07' # 140 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = false
    strain = FINITE
    decomposition_method = EigenSolution
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = MaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = MaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = MaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = MaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.0
    translation = 0.1625
    boundary = 4
  []
  [no_y_midsection]
    type = DirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[Constraints]
  [bottom_section_plane]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    primary = 101 # node on boundary
    secondary = 1 # boundary
    penalty = 1.e+10
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    number_axial_zone = 15
    oxide_thickness = scale_thickness
  []
[]

[Materials]
  [thermal]
    type = ZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = cladding
    youngs_modulus = 1.0e11
    poissons_ratio = 0.3
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [clad_zrycreep]
    type = ZryCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501 # see Erbacher et al., 1982
    temperature_standard_thermal_creep_end = 500
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = StrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[Executioner]
  type = Transient
  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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-04
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp] # average temperature of cladding exterior
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []

  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = cladding
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = cladding
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geo
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_14MPa_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  [fuel_pin_geo]
    type = FuelPinGeometry
    clad_outer_wall = '2'
    clad_inner_wall = '4'
    include_fuel = false
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_10MPa/rebeka_singlerod_2d_10MPa.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

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

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

[Variables]
  [temperature]
    initial_condition = 573.0
  []
[]

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+07  1.e+07' # 100 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = true
    strain = FINITE
    decomposition_method = EigenSolution
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = MaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = MaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = MaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = MaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.
    translation = 0.1625
    boundary = 4
  []
  [no_y_midsection]
    type = DirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[Constraints]
  [bottom_section_plane]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    primary = 101 # node on boundary
    secondary = 1 # boundary
    penalty = 1.e+10
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    oxide_thickness = scale_thickness
  []
[]

[Materials]
  [thermal]
    type = ZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ZryElasticityTensor
    block = cladding
    matpro_youngs_modulus = false
    matpro_poissons_ratio = false
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [clad_zrycreep]
    type = ZryCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501
    temperature_standard_thermal_creep_end = 500
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = StrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[Executioner]
  type = Transient
  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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-06
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp]
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [min_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = min
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = cladding
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = cladding
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geo
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
  [top_disp_r_clad] # this is mid height
    type = NodalVariableValue
    variable = disp_x
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [top_disp_r_clad_slice] # this is mid height matched to the 1.5d
    type = NodalVariableValue
    variable = disp_x
    nodeid = 3 #coords (0.005375, 0.1625)
  []
  [top_disp_z_clad]
    type = NodalVariableValue
    variable = disp_y
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [stress_xx] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_xx
    elementid = 0
  []
  [stress_yy] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_yy
    elementid = 0
  []
  [stress_zz] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_zz
    elementid = 0
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_10MPa_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  [fuel_pin_geo]
    type = FuelPinGeometry
    clad_outer_wall = '2'
    clad_inner_wall = '4'
    include_fuel = false
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_10MPa/rebeka_singlerod_2d_10MPa_aniso_hunt.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

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

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

[Variables]
  [temperature]
    initial_condition = 573.0
  []
[]

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+07  1.e+07' # 100 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
  # Circumferential steady-state creep test and analysis of Zircaloy-4 fuelcladding
  # Nuclear Engineering and Design, 53, 2312-2322 (2021)
  [F]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.738 0.57 0.57'
  []
  [G]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.174 0.45 0.45'
  []
  [H]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.588 0.48 0.48'
  []
  [L]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [M]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [N]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = true
    strain = FINITE
    decomposition_method = TaylorExpansion
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
    use_automatic_differentiation = true
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = ADMaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = ADMaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = ADMaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = ADMaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = ADMaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = ADMaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = ADMaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = ADMaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.
    translation = 0.1625
    boundary = 4
  []
  [no_y_midsection]
    type = ADDirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[Constraints]
  [bottom_section_plane]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    primary = 101 # node on boundary
    secondary = 1 # boundary
    penalty = 1.e+10
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    oxide_thickness = scale_thickness
    use_ad = true
  []
[]

[Materials]
  [thermal]
    type = ADZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ADZryElasticityTensor
    block = cladding
    matpro_youngs_modulus = false
    matpro_poissons_ratio = false
  []
  [clad_stress]
    type = ADComputeMultipleInelasticStress
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [hill_constants]
    type = ADHillConstants
    hill_constants = "0.956 0.240 0.304 1.0 1.0 1.0"
    # function_names = "F G H L M N"
    # temperature = temperature
  []
  [clad_zrycreep]
    type = ADZryAnisoCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501
    temperature_standard_thermal_creep_end = 500
    fract_beta_phase_name = 'fract_beta_phase'
  []
  [thermal_expansion]
    type = ADZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = ADStrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ADZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ADZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ADZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-06
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp]
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [min_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = min
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = cladding
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = cladding
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geo
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
  [top_disp_r_clad] # this is mid height
    type = NodalVariableValue
    variable = disp_x
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [top_disp_r_clad_slice] # this is mid height matched to the 1.5d
    type = NodalVariableValue
    variable = disp_x
    nodeid = 3 #coords (0.005375, 0.1625)
  []
  [top_disp_z_clad]
    type = NodalVariableValue
    variable = disp_y
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [stress_xx] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_xx
    elementid = 0
  []
  [stress_yy] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_yy
    elementid = 0
  []
  [stress_zz] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_zz
    elementid = 0
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_10MPa_aniso_hunt_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  [fuel_pin_geo]
    type = FuelPinGeometry
    clad_outer_wall = '2'
    clad_inner_wall = '4'
    include_fuel = false
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_06MPa/rebeka_singlerod_2d_06MPa.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

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

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

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

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '6.e+06  6.e+06' # 60 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = false
    strain = FINITE
    decomposition_method = EigenSolution
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = MaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = MaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = MaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = MaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.
    translation = 0.1625
    boundary = 4
  []
  [no_y_midsection]
    type = DirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[Constraints]
  [bottom_section_plane]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    primary = 101 # node on boundary
    secondary = 1 # boundary
    penalty = 1.e+10
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    number_axial_zone = 15
    oxide_thickness = scale_thickness
  []
[]

[Materials]
  [thermal]
    type = ZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = cladding
    youngs_modulus = 1.0e11
    poissons_ratio = 0.3
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [clad_zrycreep]
    type = ZryCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501 # see Erbacher et al., 1982
    temperature_standard_thermal_creep_end = 500
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = StrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[Executioner]
  type = Transient
  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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-04
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp] # average temperature of cladding exterior
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = cladding
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = cladding
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geo
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_06MPa_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  [fuel_pin_geo]
    type = FuelPinGeometry
    clad_outer_wall = '2'
    clad_inner_wall = '4'
    include_fuel = false
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_12MPa/rebeka_singlerod_2d_12MPa_aniso.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

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

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

[Variables]
  [temperature]
    initial_condition = 573.0
  []
[]

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '12.e+06  12.e+06' # 120 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
  # Circumferential steady-state creep test and analysis of Zircaloy-4 fuel cladding
  # Nuclear Engineering and Design, 53, 2312-2322 (2021)
  [F]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.738 0.738 0.57'
  []
  [G]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.174 0.174 0.45'
  []
  [H]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.588 0.588 0.48'
  []
  [L]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [M]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [N]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = true
    strain = FINITE
    decomposition_method = TaylorExpansion
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
    use_automatic_differentiation = true
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = ADMaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = ADMaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = ADMaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = ADMaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = ADMaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = MaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.
    translation = 0.1625
    boundary = 4
  []
  [no_y_midsection]
    type = ADDirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[Constraints]
  [bottom_section_plane]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    primary = 101 # node on boundary
    secondary = 1 # boundary
    penalty = 1.e+10
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    oxide_thickness = scale_thickness
    number_axial_zone = 15
    use_ad = true
  []
[]

[Materials]
  [converter]
    type = MaterialADConverter
    reg_props_in = 'fract_beta_phase'
    ad_props_out = 'ad_fract_beta_phase'
  []
  [thermal]
    type = ADZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    block = cladding
    youngs_modulus = 1.0e11
    poissons_ratio = 0.3
  []
  [clad_stress]
    type = ADComputeMultipleInelasticStress
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [hill_constants]
    type = ADHillConstants
    hill_constants = "0.738 0.174 0.588 1.0 1.0 1.0"
    function_names = "F G H L M N"
    temperature = temperature
  []
  [clad_zrycreep]
    type = ADZryAnisoCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501
    temperature_standard_thermal_creep_end = 500
    fract_beta_phase_name = 'ad_fract_beta_phase'
  []
  [thermal_expansion]
    type = ADZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = ADStrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ADZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[Executioner]
  type = Transient
  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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-06
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp]
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [min_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = min
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = cladding
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = cladding
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geo
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
  [top_disp_r_clad] # this is mid height
    type = NodalVariableValue
    variable = disp_x
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [top_disp_r_clad_slice] # this is mid height matched to the 1.5d
    type = NodalVariableValue
    variable = disp_x
    nodeid = 3 #coords (0.005375, 0.1625)
  []
  [top_disp_z_clad]
    type = NodalVariableValue
    variable = disp_y
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [stress_xx] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_xx
    elementid = 0
  []
  [stress_yy] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_yy
    elementid = 0
  []
  [stress_zz] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_zz
    elementid = 0
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_12MPa_aniso_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  [fuel_pin_geo]
    type = FuelPinGeometry
    clad_outer_wall = '2'
    clad_inner_wall = '4'
    include_fuel = false
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_10MPa/rebeka_singlerod_2d_10MPa_1pt5_aniso.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.
[GlobalParams]
  order = SECOND
  family = LAGRANGE
  volumetric_locking_correction = false
  displacements = 'disp_x'
[]

[Mesh]
  coord_type = RZ
  [layered1D_mesh]
    type = Layered1DMeshGenerator
    clad_mesh_density = customize
    include_fuel = false
    include_clad = true
    include_plenum = false
    pellet_outer_radius = 0.00465
    clad_gap_width = 0
    clad_thickness = 0.000725
    fuel_height = 0.1625
    slices_per_block = 10
    pellet_bottom_coor = 0.0
    nx_c = 1
  []
  patch_update_strategy = auto
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[UserObjects]
  [pin_geometry]
    type = Layered1DFuelPinGeometry
    include_fuel = false
    mesh_generator = layered1D_mesh
  []
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
[]

[Variables]
  [temperature]
    initial_condition = 573.0
  []
[]

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+07  1.e+07' # 100 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
  [clad_axial_pressure]
    type = CladdingAxialPressureFunction
    plenum_pressure = inner_pressure_func
    coolant_pressure = outer_pressure_func
    coolant_pressure_scaling_factor = 1
    fuel_pin_geometry = pin_geometry
  []
  # Circumferential steady-state creep test and analysis of Zircaloy-4 fuelcladding
  # Nuclear Engineering and Design, 53, 2312-2322 (2021)
  [F]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.738 0.57 0.57'
  []
  [G]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.174 0.45 0.45'
  []
  [H]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.588 0.48 0.48'
  []
  [L]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [M]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [N]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
[]

[Physics/SolidMechanics/Layered1D]
  [clad]
    block = clad
    add_variables = true
    strain = FINITE
    add_scalar_variables = true
    out_of_plane_strain_name = strain_yy
    fuel_pin_geometry = pin_geometry
    out_of_plane_pressure_function = clad_axial_pressure
    decomposition_method = TaylorExpansion
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    eigenstrain_names = 'clad_thermal_eigenstrain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz
      creep_strain_xx creep_strain_yy creep_strain_xy creep_strain_zz strain_xx
      strain_yy strain_zz hoop_stress'
    group_scalar_vars_in_reference_residual = false
    mesh_generator = layered1D_mesh
    use_automatic_differentiation = true
  []
[]

[AuxVariables]
  [disp_y] ## Required for easier visualization in Paraview
  []
  [disp_z] ## Required for easier visualization in Paraview
  []
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [effective_creep_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = ADMaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [effective_creep_strain]
    type = ADMaterialRealAux
    property = effective_creep_strain
    variable = effective_creep_strain
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = ADMaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = ADMaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = ADMaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = MaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  # [inner_temperature]
  #   type = REBEKADirichletBC
  #   variable = temperature
  #   function_tempm = temperature_func
  #   minimum_temperature = 573.0
  #   translation = 0.1625
  #   boundary = 5
  # []

  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 5
      function = inner_pressure_func
    []
  []
[]

# [CoolantChannel]
#  [convective_clad_surface]
#     boundary = 2
#     variable = temperature
#     inlet_temperature = 473.
#     inlet_pressure    = 1.e+05
#     inlet_massflux    = 1. # kg/m^2-sec # almost stagnant steam
#     rod_diameter      = 10.75e-03
#     rod_pitch         = 1.26e-02  # default
#     oxide_thickness = scale_thickness
#   []
# []

[Materials]
  [converter]
    type = MaterialADConverter
    reg_props_in = 'fract_beta_phase'
    ad_props_out = 'ad_fract_beta_phase'
  []
  [thermal]
    type = ADZryThermal
    block = clad
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ADZryElasticityTensor
    block = clad
    matpro_youngs_modulus = false
    matpro_poissons_ratio = false
  []
  [clad_stress]
    type = ADComputeMultipleInelasticStress
    inelastic_models = 'clad_zrycreep'
    block = clad
  []
  [hill_constants]
    type = ADHillConstants
    hill_constants = "0.738 0.174 0.588 1.0 1.0 1.0"
    function_names = "F G H L M N"
    temperature = temperature
  []
  [clad_zrycreep]
    type = ADZryAnisoCreepLOCAUpdate
    block = clad
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501
    temperature_standard_thermal_creep_end = 500
    fract_beta_phase_name = 'ad_fract_beta_phase'
  []
  [thermal_expansion]
    type = ADZryThermalExpansionMATPROEigenstrain
    block = clad
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = StrainAdjustedDensity
    block = clad
    strain_free_density = 6550
  []
  [phase]
    type = ZrPhase
    block = clad
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ADZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    # use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = 'ad_fract_beta_phase'
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[Executioner]
  type = Transient
  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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-06
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp]
    type = LayeredSideAverageValuePostprocessor
    boundary = 2
    variable = temperature
    fuel_pin_geometry = pin_geometry
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = clad
    value_type = max
    variable = temperature
  []
  [min_clad_temp]
    type = NodalExtremeValue
    block = clad
    value_type = min
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = clad
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = clad
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = clad
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = clad
    value_type = max
    variable = creep_rate_aux
  []
  [max_eff_creep_strain]
    type = ElementExtremeValue
    block = clad
    value_type = max
    variable = effective_creep_strain
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = clad
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = clad
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = clad
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = clad
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = clad
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = clad
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = pin_geometry
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
  [top_disp_r_clad_slice] # this is mid height
    type = NodalVariableValue
    variable = disp_x
    nodeid = 28 # 1 line element
  []
  [inner_pressure_func]
    type = FunctionValuePostprocessor
    function = inner_pressure_func
    execute_on = 'initial linear'
  []
  [stress_xx] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_xx
    elementid = 9 #single line element
  []
  [stress_yy] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_yy
    elementid = 9 #single line element
  []
  [stress_zz] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_zz
    elementid = 9 #single line element
  []
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_10MPa_1pt5_out_chkfile
    show = 'ave_clad_exterior_temp max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_04MPa/rebeka_singlerod_2d_04MPa.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

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

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

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

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '4.e+06  4.e+06' # 40 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = false
    strain = FINITE
    decomposition_method = EigenSolution
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = MaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = MaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = MaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = MaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.0
    translation = 0.1625
    boundary = 4
  []
  [no_y_midsection]
    type = DirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[Constraints]
  [bottom_section_plane]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    primary = 101 # node on boundary
    secondary = 1 # boundary
    penalty = 1.e+10
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    number_axial_zone = 15
    oxide_thickness = scale_thickness
  []
[]

[Materials]
  [thermal]
    type = ZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = cladding
    youngs_modulus = 1.0e11
    poissons_ratio = 0.3
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [clad_zrycreep]
    type = ZryCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501 # see Erbacher et al., 1982
    temperature_standard_thermal_creep_end = 500
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = StrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[Executioner]
  type = Transient
  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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-04
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp] # average temperature of cladding exterior
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = cladding
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = cladding
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geo
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_04MPa_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  [fuel_pin_geo]
    type = FuelPinGeometry
    clad_outer_wall = '2'
    clad_inner_wall = '4'
    include_fuel = false
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_04MPa/rebeka_singlerod_2d_04MPa_aniso.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

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

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

[Variables]
  [temperature]
    initial_condition = 573.0
  []
[]

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '4.e+06  4.e+06' # 10 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
  # Circumferential steady-state creep test and analysis of Zircaloy-4 fuel cladding
  # Nuclear Engineering and Design, 53, 2312-2322 (2021)
  [F]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.738 0.738 0.57'
  []
  [G]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.174 0.174 0.45'
  []
  [H]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.588 0.588 0.48'
  []
  [L]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [M]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [N]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = true
    strain = FINITE
    decomposition_method = TaylorExpansion
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
    use_automatic_differentiation = true
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = ADMaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = ADMaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = ADMaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = ADMaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = ADMaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = ADMaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = ADMaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = ADMaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.
    translation = 0.1625
    boundary = 4
  []
  [no_y_midsection]
    type = ADDirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[Constraints]
  [bottom_section_plane]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    primary = 101 # node on boundary
    secondary = 1 # boundary
    penalty = 1.e+10
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    oxide_thickness = scale_thickness
    number_axial_zone = 15
    use_ad = true
  []
[]

[Materials]
  [thermal]
    type = ADZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    block = cladding
    youngs_modulus = 1.0e11
    poissons_ratio = 0.3
  []
  [clad_stress]
    type = ADComputeMultipleInelasticStress
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [hill_constants]
    type = ADHillConstants
    hill_constants = "0.738 0.174 0.588 1.0 1.0 1.0"
    function_names = "F G H L M N"
    temperature = temperature
  []
  [clad_zrycreep]
    type = ADZryAnisoCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501
    temperature_standard_thermal_creep_end = 500
    fract_beta_phase_name = 'fract_beta_phase'
  []
  [thermal_expansion]
    type = ADZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = ADStrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ADZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ADZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ADZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-06
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp]
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [min_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = min
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = cladding
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = cladding
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geo
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
  [top_disp_r_clad] # this is mid height
    type = NodalVariableValue
    variable = disp_x
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [top_disp_r_clad_slice] # this is mid height matched to the 1.5d
    type = NodalVariableValue
    variable = disp_x
    nodeid = 3 #coords (0.005375, 0.1625)
  []
  [top_disp_z_clad]
    type = NodalVariableValue
    variable = disp_y
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [stress_xx] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_xx
    elementid = 0
  []
  [stress_yy] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_yy
    elementid = 0
  []
  [stress_zz] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_zz
    elementid = 0
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_04MPa_aniso_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  [fuel_pin_geo]
    type = FuelPinGeometry
    clad_outer_wall = '2'
    clad_inner_wall = '4'
    include_fuel = false
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_06MPa/rebeka_singlerod_2d_06MPa_aniso.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

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

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

[Variables]
  [temperature]
    initial_condition = 573.0
  []
[]

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '6.e+06  6.e+06' # 10 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
  # Circumferential steady-state creep test and analysis of Zircaloy-4 fuel cladding
  # Nuclear Engineering and Design, 53, 2312-2322 (2021)
  [F]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.738 0.738 0.57'
  []
  [G]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.174 0.174 0.45'
  []
  [H]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.588 0.588 0.48'
  []
  [L]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [M]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [N]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = true
    strain = FINITE
    decomposition_method = TaylorExpansion
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
    use_automatic_differentiation = true
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = ADMaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = ADMaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = ADMaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = ADMaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = ADMaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = MaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.
    translation = 0.1625
    boundary = 4
  []
  [no_y_midsection]
    type = ADDirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[Constraints]
  [bottom_section_plane]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    primary = 101 # node on boundary
    secondary = 1 # boundary
    penalty = 1.e+10
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    oxide_thickness = scale_thickness
    number_axial_zone = 15
    use_ad = true
  []
[]

[Materials]
  [converter]
    type = MaterialADConverter
    reg_props_in = 'fract_beta_phase'
    ad_props_out = 'ad_fract_beta_phase'
  []
  [thermal]
    type = ADZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    block = cladding
    youngs_modulus = 1.0e11
    poissons_ratio = 0.3
  []
  [clad_stress]
    type = ADComputeMultipleInelasticStress
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [hill_constants]
    type = ADHillConstants
    hill_constants = "0.738 0.174 0.588 1.0 1.0 1.0"
    function_names = "F G H L M N"
    temperature = temperature
  []
  [clad_zrycreep]
    type = ADZryAnisoCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501
    temperature_standard_thermal_creep_end = 500
    fract_beta_phase_name = 'ad_fract_beta_phase'
  []
  [thermal_expansion]
    type = ADZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = ADStrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ADZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[Executioner]
  type = Transient
  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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-06
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp]
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [min_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = min
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = cladding
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = cladding
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geo
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
  [top_disp_r_clad] # this is mid height
    type = NodalVariableValue
    variable = disp_x
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [top_disp_r_clad_slice] # this is mid height matched to the 1.5d
    type = NodalVariableValue
    variable = disp_x
    nodeid = 3 #coords (0.005375, 0.1625)
  []
  [top_disp_z_clad]
    type = NodalVariableValue
    variable = disp_y
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [stress_xx] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_xx
    elementid = 0
  []
  [stress_yy] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_yy
    elementid = 0
  []
  [stress_zz] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_zz
    elementid = 0
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_06MPa_aniso_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  [fuel_pin_geo]
    type = FuelPinGeometry
    clad_outer_wall = '2'
    clad_inner_wall = '4'
    include_fuel = false
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_01MPa/rebeka_singlerod_2d_01MPa.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

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

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

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

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+06  1.e+06' # 10 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = false
    strain = FINITE
    decomposition_method = EigenSolution
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = MaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = MaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = MaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = MaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.
    translation = 0.1625
    boundary = 4
  []
  [no_y_midsection]
    type = DirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[Constraints]
  [bottom_section_plane]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    primary = 101 # node on boundary
    secondary = 1 # boundary
    penalty = 1.e+10
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02
    number_axial_zone = 15
    oxide_thickness = scale_thickness
  []
[]

[Materials]
  [thermal]
    type = ZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = cladding
    youngs_modulus = 1.0e11
    poissons_ratio = 0.3
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [clad_zrycreep]
    type = ZryCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501 # see Erbacher et al., 1982
    temperature_standard_thermal_creep_end = 500
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = StrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[Executioner]
  type = Transient
  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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-04
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp] # average temperature of cladding exterior
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = cladding
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = cladding
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geo
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_01MPa_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  [fuel_pin_geo]
    type = FuelPinGeometry
    clad_outer_wall = '2'
    clad_inner_wall = '4'
    include_fuel = false
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_02MPa/rebeka_singlerod_2d_02MPa_aniso.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

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

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

[Variables]
  [temperature]
    initial_condition = 573.0
  []
[]

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '2.e+06  2.e+06' # 10 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
  # Circumferential steady-state creep test and analysis of Zircaloy-4 fuel cladding
  # Nuclear Engineering and Design, 53, 2312-2322 (2021)
  [F]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.738 0.738 0.57'
  []
  [G]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.174 0.174 0.45'
  []
  [H]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.588 0.588 0.48'
  []
  [L]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [M]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [N]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = true
    strain = FINITE
    decomposition_method = TaylorExpansion
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
    use_automatic_differentiation = true
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = ADMaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = ADMaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = ADMaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = ADMaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = ADMaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = ADMaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = ADMaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = ADMaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.
    translation = 0.1625
    boundary = 4
  []
  [no_y_midsection]
    type = ADDirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[Constraints]
  [bottom_section_plane]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    primary = 101 # node on boundary
    secondary = 1 # boundary
    penalty = 1.e+10
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    oxide_thickness = scale_thickness
    number_axial_zone = 15
    use_ad = true
  []
[]

[Materials]
  [thermal]
    type = ADZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    block = cladding
    youngs_modulus = 1.0e11
    poissons_ratio = 0.3
  []
  [clad_stress]
    type = ADComputeMultipleInelasticStress
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [hill_constants]
    type = ADHillConstants
    hill_constants = "0.738 0.174 0.588 1.0 1.0 1.0"
    function_names = "F G H L M N"
    temperature = temperature
  []
  [clad_zrycreep]
    type = ADZryAnisoCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501
    temperature_standard_thermal_creep_end = 500
    fract_beta_phase_name = 'fract_beta_phase'
  []
  [thermal_expansion]
    type = ADZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = ADStrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ADZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ADZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ADZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-06
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp]
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [min_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = min
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = cladding
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = cladding
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geo
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
  [top_disp_r_clad] # this is mid height
    type = NodalVariableValue
    variable = disp_x
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [top_disp_r_clad_slice] # this is mid height matched to the 1.5d
    type = NodalVariableValue
    variable = disp_x
    nodeid = 3 #coords (0.005375, 0.1625)
  []
  [top_disp_z_clad]
    type = NodalVariableValue
    variable = disp_y
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [stress_xx] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_xx
    elementid = 0
  []
  [stress_yy] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_yy
    elementid = 0
  []
  [stress_zz] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_zz
    elementid = 0
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_02MPa_aniso_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  [fuel_pin_geo]
    type = FuelPinGeometry
    clad_outer_wall = '2'
    clad_inner_wall = '4'
    include_fuel = false
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_12MPa/rebeka_singlerod_2d_12MPa.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

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

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

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

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.2e+07  1.2e+07' # 120 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = false
    strain = FINITE
    decomposition_method = EigenSolution
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = MaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = MaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = MaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = MaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.0
    translation = 0.1625
    boundary = 4
  []
  [no_y_midsection]
    type = DirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[Constraints]
  [bottom_section_plane]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    primary = 101 # node on boundary
    secondary = 1 # boundary
    penalty = 1.e+10
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    number_axial_zone = 15
    oxide_thickness = scale_thickness
  []
[]

[Materials]
  [thermal]
    type = ZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = cladding
    youngs_modulus = 1.0e11
    poissons_ratio = 0.3
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [clad_zrycreep]
    type = ZryCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501 # see Erbacher et al., 1982
    temperature_standard_thermal_creep_end = 500
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = StrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[Executioner]
  type = Transient
  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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-04
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp] # average temperature of cladding exterior
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = cladding
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = cladding
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geo
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_12MPa_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  [fuel_pin_geo]
    type = FuelPinGeometry
    clad_outer_wall = '2'
    clad_inner_wall = '4'
    include_fuel = false
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_10MPa/rebeka_singlerod_2d_10MPa_aniso.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

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

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

[Variables]
  [temperature]
    initial_condition = 573.0
  []
[]

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+07  1.e+07' # 100 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
  # Circumferential steady-state creep test and analysis of Zircaloy-4 fuel cladding
  # Nuclear Engineering and Design, 53, 2312-2322 (2021)
  [F]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.738 0.738 0.57'
  []
  [G]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.174 0.174 0.45'
  []
  [H]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.588 0.588 0.48'
  []
  [L]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [M]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [N]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = true
    strain = FINITE
    decomposition_method = TaylorExpansion
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
    use_automatic_differentiation = true
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = ADMaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = ADMaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = ADMaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = ADMaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = ADMaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = ADMaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = ADMaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = ADMaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.
    translation = 0.1625
    boundary = 4
  []
  [no_y_midsection]
    type = ADDirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[Constraints]
  [bottom_section_plane]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    primary = 101 # node on boundary
    secondary = 1 # boundary
    penalty = 1.e+10
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    oxide_thickness = scale_thickness
    use_ad = true
  []
[]

[Materials]
  [thermal]
    type = ADZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ADZryElasticityTensor
    block = cladding
    matpro_youngs_modulus = false
    matpro_poissons_ratio = false
  []
  [clad_stress]
    type = ADComputeMultipleInelasticStress
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [hill_constants]
    type = ADHillConstants
    hill_constants = "0.738 0.174 0.588 1.0 1.0 1.0"
    function_names = "F G H L M N"
    temperature = temperature
  []
  [clad_zrycreep]
    type = ADZryAnisoCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501
    temperature_standard_thermal_creep_end = 500
    fract_beta_phase_name = 'fract_beta_phase'
  []
  [thermal_expansion]
    type = ADZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = ADStrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ADZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ADZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ADZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-06
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp]
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [min_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = min
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = cladding
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = cladding
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geo
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
  [top_disp_r_clad] # this is mid height
    type = NodalVariableValue
    variable = disp_x
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [top_disp_r_clad_slice] # this is mid height matched to the 1.5d
    type = NodalVariableValue
    variable = disp_x
    nodeid = 3 #coords (0.005375, 0.1625)
  []
  [top_disp_z_clad]
    type = NodalVariableValue
    variable = disp_y
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [stress_xx] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_xx
    elementid = 0
  []
  [stress_yy] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_yy
    elementid = 0
  []
  [stress_zz] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_zz
    elementid = 0
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_10MPa_aniso_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  [fuel_pin_geo]
    type = FuelPinGeometry
    clad_outer_wall = '2'
    clad_inner_wall = '4'
    include_fuel = false
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_14MPa/rebeka_singlerod_2d_14MPa_aniso.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

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

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

[Variables]
  [temperature]
    initial_condition = 573.0
  []
[]

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '14.e+06  14.e+06' # 140 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
  # Circumferential steady-state creep test and analysis of Zircaloy-4 fuel cladding
  # Nuclear Engineering and Design, 53, 2312-2322 (2021)
  [F]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.738 0.738 0.57'
  []
  [G]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.174 0.174 0.45'
  []
  [H]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.588 0.588 0.48'
  []
  [L]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [M]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [N]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = true
    strain = FINITE
    decomposition_method = TaylorExpansion
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
    use_automatic_differentiation = true
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = ADMaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = ADMaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = ADMaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = ADMaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = ADMaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = MaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.
    translation = 0.1625
    boundary = 4
  []
  [no_y_midsection]
    type = ADDirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[Constraints]
  [bottom_section_plane]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    primary = 101 # node on boundary
    secondary = 1 # boundary
    penalty = 1.e+10
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    oxide_thickness = scale_thickness
    number_axial_zone = 15
    use_ad = true
  []
[]

[Materials]
  [converter]
    type = MaterialADConverter
    reg_props_in = 'fract_beta_phase'
    ad_props_out = 'ad_fract_beta_phase'
  []
  [thermal]
    type = ADZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    block = cladding
    youngs_modulus = 1.0e11
    poissons_ratio = 0.3
  []
  [clad_stress]
    type = ADComputeMultipleInelasticStress
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [hill_constants]
    type = ADHillConstants
    hill_constants = "0.738 0.174 0.588 1.0 1.0 1.0"
    function_names = "F G H L M N"
    temperature = temperature
  []
  [clad_zrycreep]
    type = ADZryAnisoCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501
    temperature_standard_thermal_creep_end = 500
    fract_beta_phase_name = 'ad_fract_beta_phase'
  []
  [thermal_expansion]
    type = ADZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = ADStrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ADZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[Executioner]
  type = Transient
  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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-06
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp]
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [min_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = min
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = cladding
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = cladding
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geo
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
  [top_disp_r_clad] # this is mid height
    type = NodalVariableValue
    variable = disp_x
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [top_disp_r_clad_slice] # this is mid height matched to the 1.5d
    type = NodalVariableValue
    variable = disp_x
    nodeid = 3 #coords (0.005375, 0.1625)
  []
  [top_disp_z_clad]
    type = NodalVariableValue
    variable = disp_y
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [stress_xx] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_xx
    elementid = 0
  []
  [stress_yy] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_yy
    elementid = 0
  []
  [stress_zz] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_zz
    elementid = 0
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_14MPa_aniso_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  [fuel_pin_geo]
    type = FuelPinGeometry
    clad_outer_wall = '2'
    clad_inner_wall = '4'
    include_fuel = false
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_3d_10MPa_01Kpers_90deg/rebeka_singlerod_3d_10MPa_90deg.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
  order = SECOND
  family = LAGRANGE
[]

[Problem]
[]

[Mesh]
  [mesh]
    type = FileMeshGenerator
    file = 3d_rebeka_singlerod.e
  []
[]

[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
  [temperature]
    initial_condition = 573.0
  []
[]

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+07  1.e+07' # 100 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = false
    strain = FINITE
    decomposition_method = EigenSolution
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase] # Fraction of beta phase in Zry
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness] # ZrO2 scale thickness (m)
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total] # Total oxigen weight fraction (oxide+metal) (/)
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress] # Hoop stress at cladding burst
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = MaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = MaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = MaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = MaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    translation = 0.1625
    vary_azimuthally = true
    azimuthal_variation_pi2 = -30.
    minimum_temperature = 573.0
    boundary = 4
  []
  [no_y_midsection]
    type = DirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [no_z_frontside]
    type = DirichletBC
    variable = disp_z
    boundary = 5
    value = 0.
  []
  [no_x_backside]
    type = DirichletBC
    variable = disp_x
    boundary = 6
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    number_axial_zone = 15
    oxide_thickness = scale_thickness
  []
[]

[Materials]
  [thermal]
    type = ZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = cladding
    youngs_modulus = 1.0e11
    poissons_ratio = 0.3
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [clad_zrycreep]
    type = ZryCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501 # see Erbacher et al., 1982
    temperature_standard_thermal_creep_end = 500
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = StrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[Executioner]
  type = Transient
  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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-04
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp]
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_singlerod_3d_10MPa_90deg_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_10MPa/rebeka_singlerod_2d_10MPa_1pt5.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.
[GlobalParams]
  order = SECOND
  family = LAGRANGE
  volumetric_locking_correction = false
  displacements = 'disp_x'
[]

[Mesh]
  coord_type = RZ
  [layered1D_mesh]
    type = Layered1DMeshGenerator
    clad_mesh_density = customize
    include_fuel = false
    include_clad = true
    include_plenum = false
    pellet_outer_radius = 0.00465
    clad_gap_width = 0
    clad_thickness = 0.000725
    fuel_height = 0.1625
    slices_per_block = 10
    pellet_bottom_coor = 0.0
    nx_c = 1
  []
  patch_update_strategy = auto
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[UserObjects]
  [pin_geometry]
    type = Layered1DFuelPinGeometry
    include_fuel = false
    mesh_generator = layered1D_mesh
  []
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
[]

[Variables]
  [temperature]
    initial_condition = 573.0
  []
[]

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+07  1.e+07' # 100 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
  [clad_axial_pressure]
    type = CladdingAxialPressureFunction
    plenum_pressure = inner_pressure_func
    coolant_pressure = outer_pressure_func
    coolant_pressure_scaling_factor = 1
    fuel_pin_geometry = pin_geometry
  []
[]

[Physics/SolidMechanics/Layered1D]
  [clad]
    block = clad
    add_variables = true
    strain = FINITE
    add_scalar_variables = true
    out_of_plane_strain_name = strain_yy
    fuel_pin_geometry = pin_geometry
    out_of_plane_pressure_function = clad_axial_pressure
    decomposition_method = EigenSolution
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    eigenstrain_names = 'clad_thermal_eigenstrain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz
      creep_strain_xx creep_strain_yy creep_strain_xy creep_strain_zz strain_xx
      strain_yy strain_zz hoop_stress'
    group_scalar_vars_in_reference_residual = false
    mesh_generator = layered1D_mesh
  []
[]

[AuxVariables]
  [disp_y] ## Required for easier visualization in Paraview
  []
  [disp_z] ## Required for easier visualization in Paraview
  []
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [effective_creep_strain]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = MaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [effective_creep_strain]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = effective_creep_strain
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = MaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = MaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = MaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.0
    translation = 0.1625
    boundary = 5
  []

  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 5
      function = inner_pressure_func
    []
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1. # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    oxide_thickness = scale_thickness
  []
[]

[Materials]
  [thermal]
    type = ZryThermal
    block = clad
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ZryElasticityTensor
    block = clad
    matpro_youngs_modulus = false
    matpro_poissons_ratio = false
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep'
    block = clad
  []
  [clad_zrycreep]
    type = ZryCreepLOCAUpdate
    block = clad
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501
    temperature_standard_thermal_creep_end = 500
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = clad
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = StrainAdjustedDensity
    block = clad
    strain_free_density = 6550
  []
  [phase]
    type = ZrPhase
    block = clad
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[Executioner]
  type = Transient
  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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-06
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp]
    type = LayeredSideAverageValuePostprocessor
    boundary = 2
    variable = temperature
    fuel_pin_geometry = pin_geometry
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = clad
    value_type = max
    variable = temperature
  []
  [min_clad_temp]
    type = NodalExtremeValue
    block = clad
    value_type = min
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = clad
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = clad
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = clad
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = clad
    value_type = max
    variable = creep_rate_aux
  []
  [max_eff_creep_strain]
    type = ElementExtremeValue
    block = clad
    value_type = max
    variable = effective_creep_strain
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = clad
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = clad
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = clad
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = clad
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = clad
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = clad
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = pin_geometry
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
  [top_disp_r_clad_slice] # this is mid height
    type = NodalVariableValue
    variable = disp_x
    nodeid = 28 # 1 line element
  []
  [inner_pressure_func]
    type = FunctionValuePostprocessor
    function = inner_pressure_func
    execute_on = 'initial linear'
  []
  [stress_xx] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_xx
    elementid = 9 #single line element
  []
  [stress_yy] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_yy
    elementid = 9 #single line element
  []
  [stress_zz] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_zz
    elementid = 9 #single line element
  []
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_10MPa_1pt5_out_chkfile
    show = 'ave_clad_exterior_temp max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_08MPa/rebeka_singlerod_2d_08MPa_aniso.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

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

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

[Variables]
  [temperature]
    initial_condition = 573.0
  []
[]

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '8.e+06  8.e+06' # 80 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
  # Circumferential steady-state creep test and analysis of Zircaloy-4 fuel cladding
  # Nuclear Engineering and Design, 53, 2312-2322 (2021)
  [F]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.738 0.738 0.57'
  []
  [G]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.174 0.174 0.45'
  []
  [H]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.588 0.588 0.48'
  []
  [L]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [M]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [N]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = true
    strain = FINITE
    decomposition_method = TaylorExpansion
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
    use_automatic_differentiation = true
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = ADMaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = ADMaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = ADMaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = ADMaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = ADMaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = MaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.
    translation = 0.1625
    boundary = 4
  []
  [no_y_midsection]
    type = ADDirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[Constraints]
  [bottom_section_plane]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    primary = 101 # node on boundary
    secondary = 1 # boundary
    penalty = 1.e+10
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    oxide_thickness = scale_thickness
    number_axial_zone = 15
    use_ad = true
  []
[]

[Materials]
  [converter]
    type = MaterialADConverter
    reg_props_in = 'fract_beta_phase'
    ad_props_out = 'ad_fract_beta_phase'
  []
  [thermal]
    type = ADZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    block = cladding
    youngs_modulus = 1.0e11
    poissons_ratio = 0.3
  []
  [clad_stress]
    type = ADComputeMultipleInelasticStress
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [hill_constants]
    type = ADHillConstants
    hill_constants = "0.738 0.174 0.588 1.0 1.0 1.0"
    function_names = "F G H L M N"
    temperature = temperature
  []
  [clad_zrycreep]
    type = ADZryAnisoCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501
    temperature_standard_thermal_creep_end = 500
    fract_beta_phase_name = 'ad_fract_beta_phase'
  []
  [thermal_expansion]
    type = ADZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = ADStrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ADZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[Executioner]
  type = Transient
  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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-06
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp]
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [min_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = min
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = cladding
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = cladding
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geo
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
  [top_disp_r_clad] # this is mid height
    type = NodalVariableValue
    variable = disp_x
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [top_disp_r_clad_slice] # this is mid height matched to the 1.5d
    type = NodalVariableValue
    variable = disp_x
    nodeid = 3 #coords (0.005375, 0.1625)
  []
  [top_disp_z_clad]
    type = NodalVariableValue
    variable = disp_y
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [stress_xx] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_xx
    elementid = 0
  []
  [stress_yy] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_yy
    elementid = 0
  []
  [stress_zz] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_zz
    elementid = 0
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_08MPa_aniso_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  [fuel_pin_geo]
    type = FuelPinGeometry
    clad_outer_wall = '2'
    clad_inner_wall = '4'
    include_fuel = false
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_02MPa/rebeka_singlerod_2d_02MPa.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

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

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

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

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '2.e+06  2.e+06' # 20 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = false
    strain = FINITE
    decomposition_method = EigenSolution
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = MaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = MaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = MaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = MaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = MaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = MaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = MaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = MaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.
    translation = 0.1625
    boundary = 4
  []
  [no_y_midsection]
    type = DirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []

  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[Constraints]
  [bottom_section_plane]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    primary = 101 # node on boundary
    secondary = 1 # boundary
    penalty = 1.e+10
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    number_axial_zone = 15
    oxide_thickness = scale_thickness
  []
[]

[Materials]
  [thermal]
    type = ZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = cladding
    youngs_modulus = 1.0e11
    poissons_ratio = 0.3
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [clad_zrycreep]
    type = ZryCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501 # see Erbacher et al., 1982
    temperature_standard_thermal_creep_end = 500
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = StrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[Executioner]
  type = Transient
  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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-04
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp] # average temperature of cladding exterior
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = cladding
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = cladding
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geo
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_02MPa_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  [fuel_pin_geo]
    type = FuelPinGeometry
    clad_outer_wall = '2'
    clad_inner_wall = '4'
    include_fuel = false
  []
[]

(assessment/LWR/validation/LOCA_REBEKA_cladding_burst_tests/analysis/rebeka_2d_01MPa/rebeka_singlerod_2d_01MPa_aniso.i)

# Simulation of REBEKA single-rod, cladding-only LOCA tests.

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

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

[Variables]
  [temperature]
    initial_condition = 573.0
  []
[]

[Functions]
  [temperature_func]
    type = PiecewiseLinear
    x = '0.   700. '
    y = '573. 1273.' # From 300 to 1000 C at 1 K/s
  []
  [inner_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+06  1.e+06' # 10 bar
  []
  [outer_pressure_func]
    type = PiecewiseLinear
    x = '0.      700.  '
    y = '1.e+05  1.e+05' # atmospheric pressure
  []
  # Circumferential steady-state creep test and analysis of Zircaloy-4 fuel cladding
  # Nuclear Engineering and Design, 53, 2312-2322 (2021)
  [F]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.738 0.738 0.57'
  []
  [G]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.174 0.174 0.45'
  []
  [H]
    type = PiecewiseLinear
    x = '293 1073 1400'
    y = '0.588 0.588 0.48'
  []
  [L]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [M]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
  [N]
    type = PiecewiseLinear
    x = '-1000 10000'
    y = '1.0 1.0'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [clad]
    block = cladding
    add_variables = true
    strain = FINITE
    decomposition_method = TaylorExpansion
    eigenstrain_names = 'clad_thermal_eigenstrain'
    cylindrical_axis_point1 = '0 0 0'
    cylindrical_axis_point2 = '0 1 0'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz strain_xx
      strain_yy strain_zz hoop_stress'
    use_automatic_differentiation = true
  []
[]

[AuxVariables]
  [creep_rate_aux]
    order = CONSTANT
    family = MONOMIAL
  []
  [creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [fract_beta_phase]
    order = CONSTANT
    family = MONOMIAL
  []
  [scale_thickness]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfract_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [oxywtfgain_total]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst_stress]
    order = CONSTANT
    family = MONOMIAL
  []
  [burst]
    order = CONSTANT
    family = MONOMIAL
  []
[]

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

[AuxKernels]
  [creep_rate_aux]
    type = ADMaterialRealAux
    variable = creep_rate_aux
    property = creep_rate
    execute_on = timestep_end
  []
  [creep_strain_mag]
    type = ADMaterialRealAux
    property = effective_creep_strain
    variable = creep_strain_mag
    execute_on = timestep_end
  []
  [fract_bphase]
    type = ADMaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
  []
  [scl_thickness]
    type = ADMaterialRealAux
    variable = scale_thickness
    property = oxide_scale_thickness
    boundary = 2
  []
  [ofract_total]
    type = ADMaterialRealAux
    variable = oxywtfract_total
    property = current_oxygen_weight_frac_total
    boundary = 2
  []
  [ofgain_total]
    type = ADMaterialRealAux
    boundary = 2
    variable = oxywtfgain_total
    property = oxygen_weight_frac_gained_total
  []
  [sigmaburst]
    type = ADMaterialRealAux
    variable = burst_stress
    property = burst_stress
    boundary = 2
  []
  [hasburst]
    type = ADMaterialRealAux
    variable = burst
    property = failed
    boundary = 2
    execute_on = timestep_end
  []
[]

[BCs]
  [inner_temperature]
    type = REBEKADirichletBC
    variable = temperature
    function_tempm = temperature_func
    minimum_temperature = 573.
    translation = 0.1625
    boundary = 4
  []
  [no_y_midsection]
    type = ADDirichletBC
    variable = disp_y
    boundary = 3
    value = 0.
  []
  [Pressure]
    [outer_pressure] #  apply steam pressure on clad outer wall
      boundary = 2
      function = outer_pressure_func
    []
    [inner_pressure] #  apply He pressure on clad inner wall
      boundary = 4
      function = inner_pressure_func
    []
  []
[]

[Constraints]
  [bottom_section_plane]
    type = EqualValueBoundaryConstraint
    variable = disp_y
    primary = 101 # node on boundary
    secondary = 1 # boundary
    penalty = 1.e+10
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temperature
    inlet_temperature = 473.
    inlet_pressure = 1.e+05
    inlet_massflux = 1.0 # kg/m^2-sec # almost stagnant steam
    rod_diameter = 10.75e-03
    rod_pitch = 1.26e-02 # default
    oxide_thickness = scale_thickness
    number_axial_zone = 15
    use_ad = true
  []
[]

[Materials]
  [thermal]
    type = ADZryThermal
    block = cladding
    temperature = temperature
  []
  [clad_elasticity_tensor]
    type = ADComputeIsotropicElasticityTensor
    block = cladding
    youngs_modulus = 1.0e11
    poissons_ratio = 0.3
  []
  [clad_stress]
    type = ADComputeMultipleInelasticStress
    inelastic_models = 'clad_zrycreep'
    block = cladding
  []
  [hill_constants]
    type = ADHillConstants
    hill_constants = "0.738 0.174 0.588 1.0 1.0 1.0"
    function_names = "F G H L M N"
    temperature = temperature
  []
  [clad_zrycreep]
    type = ADZryAnisoCreepLOCAUpdate
    block = cladding
    temperature = temperature
    model_irradiation_creep = false
    model_primary_creep = false
    model_thermal_creep = true
    temperature_loca_creep_begin = 501
    temperature_standard_thermal_creep_end = 500
    fract_beta_phase_name = 'fract_beta_phase'
  []
  [thermal_expansion]
    type = ADZryThermalExpansionMATPROEigenstrain
    block = cladding
    temperature = temperature
    stress_free_temperature = 573.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [density]
    type = ADStrainAdjustedDensity
    block = cladding
    strain_free_density = 6550
  []
  [phase]
    type = ADZrPhase
    block = cladding
    temperature = temperature
    numerical_method = 2
  []
  [oxidation]
    type = ADZryOxidation
    boundary = 2
    temperature = temperature
    clad_inner_radius = 0.00465
    clad_outer_radius = 0.005375
    normal_operating_temperature_model = epri_kwu_ce
    high_temperature_model = leistikow
    use_coolant_channel = true
  []
  [clad_failure_criterion]
    type = ADZryCladdingFailure
    boundary = 2
    failure_criterion = combined_overstress_and_plastic_instability
    hoop_stress = hoop_stress
    effective_strain_rate_creep = creep_rate_aux
    temperature = temperature
    fraction_beta_phase = fract_beta_phase
    fraction_oxygen_gain = oxywtfgain_total
  []
[]

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

[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 = 100
  l_tol = 1e-2
  nl_max_its = 15
  nl_rel_tol = 1.e-06
  nl_abs_tol = 1.e-08
  start_time = 0
  n_startup_steps = 1
  end_time = 700.
  dtmax = 1.
  dtmin = 0.00000001

  [TimeStepper]
    type = PostprocessorDT
    postprocessor = material_timestep
    dt = 1.
  []
[]

[Postprocessors]
  [ave_clad_exterior_temp]
    type = SideAverageValue
    boundary = 2
    variable = temperature
  []
  [max_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = max
    variable = temperature
  []
  [min_clad_temp]
    type = NodalExtremeValue
    block = cladding
    value_type = min
    variable = temperature
  []
  [material_timestep]
    type = MaterialTimeStepPostprocessor
    block = cladding
  []
  [max_oxygen_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = oxywtfract_total
  []
  [max_betaph_fract]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = fract_beta_phase
  []
  [max_creep_rate]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_rate_aux
  []
  [max_creep_strain_mag]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = creep_strain_mag
  []
  [max_hoop_strain]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = strain_zz
  []
  [max_hoop_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = hoop_stress
  []
  [vonmises_stress_clad]
    type = ElementAverageValue
    block = cladding
    variable = vonmises_stress
  []
  [min_burst_stress]
    type = ElementExtremeValue
    block = cladding
    value_type = min
    variable = burst_stress
  []
  [burst]
    type = ElementExtremeValue
    block = cladding
    value_type = max
    variable = burst
  []
  [peak_hoop_strain]
    type = ElementExtremeValue
    value_type = max
    variable = strain_zz
    block = cladding
  []
  [zry_burst_opening_area]
    type = ZryBurstOpening
    fuel_pin_geometry = fuel_pin_geo
    peak_hoop_strain = peak_hoop_strain
    estimate = limiting
    opening_shape = rectangle
    output = area
  []
  [top_disp_r_clad] # this is mid height
    type = NodalVariableValue
    variable = disp_x
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [top_disp_r_clad_slice] # this is mid height matched to the 1.5d
    type = NodalVariableValue
    variable = disp_x
    nodeid = 3 #coords (0.005375, 0.1625)
  []
  [top_disp_z_clad]
    type = NodalVariableValue
    variable = disp_y
    nodeid = 0 #coords (0.005375, 0.1625)
  []
  [stress_xx] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_xx
    elementid = 0
  []
  [stress_yy] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_yy
    elementid = 0
  []
  [stress_zz] # stess in the top Element
    type = ElementalVariableValue
    variable = stress_zz
    elementid = 0
  []
[]

[StandardLWRFuelRodOutputs]
  rod_component = clad
  cladding_blocks = cladding
  plenum_boundary_name = 4
  external_clad_boundary_name = 2
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  csv = true
  [console]
    type = Console
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    file_base = rebeka_2d_01MPa_aniso_out_chkfile
    show = 'average_interior_clad_temperature max_hoop_strain max_hoop_stress'
    execute_on = 'FINAL'
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = 'burst > 0'
  []
  [fuel_pin_geo]
    type = FuelPinGeometry
    clad_outer_wall = '2'
    clad_inner_wall = '4'
    include_fuel = false
  []
[]

Description
Example Input Syntax
Input Parameters
Input Files
References