3D Light Water Reactor Fuel Rod Tutorial

In capturing phenomena that break the axisymmetry of typical LWR problems, three-dimensional models can be employed. Examples of these applications encompass pellets with missing surfaces or arbitrary misalignment. Two approaches to thermomechanical contact can be used in three-dimensions. Node (or point) to surface (NTS) enforcement and a surface to surface (mortar) enforcement.

Node-to-surface input files

Below is an input file containing a setup for an NTS thermomechanical enforcement of pellets whose geometry is assumed to be smeared. I.e., no discrete geometry effects are considered in the finite element mesh. These examples use frictionless contact.

initial_fuel_density = 10431.0

[GlobalParams<<<{"href": "../syntax/GlobalParams/index.html"}>>>]
  density = ${initial_fuel_density} # initial fuel density 95.0% of theoretical (10980 kg/m3)
  displacements = 'disp_x disp_y disp_z'
  order = SECOND
  family = LAGRANGE
  energy_per_fission = 3.2e-11 # J/fission
[]

[Mesh<<<{"href": "../syntax/Mesh/index.html"}>>>]
  patch_size = 20
  patch_update_strategy = iteration
  partitioner = centroid
  centroid_partitioner_direction = y
  [mesh]
    type = FileMeshGenerator<<<{"description": "Read a mesh from a file.", "href": "../source/meshgenerators/FileMeshGenerator.html"}>>>
    file<<<{"description": "The filename to read."}>>> = smearedTest3.e
  []
[]

[Problem<<<{"href": "../syntax/Problem/index.html"}>>>]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
  group_variables = 'disp_x disp_y disp_z'
[]

[Variables<<<{"href": "../syntax/Variables/index.html"}>>>]
  [temp]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 580
  []
[]

[AuxVariables<<<{"href": "../syntax/AuxVariables/index.html"}>>>]
  [fast_neutron_flux]
    block = clad
  []
  [fast_neutron_fluence]
    block = clad
  []
  [grain_radius]
    block = pellet_type_1
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 10e-6
  []
  [gap_cond]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
  []
  [coolant_htc]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
  []
  [hoop_inelastic_strain]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    block = clad
  []
[]

[Functions<<<{"href": "../syntax/Functions/index.html"}>>>]
  [power_history]
    type = PiecewiseLinear<<<{"description": "Linearly interpolates between pairs of x-y data", "href": "../source/functions/PiecewiseLinear.html"}>>>
    data_file<<<{"description": "File holding CSV data"}>>> = powerhistory.csv
    scale_factor<<<{"description": "Scale factor to be applied to the ordinate values"}>>> = 1
  []
  [axial_peaking_factors]
    type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../source/functions/MooseParsedFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = 1
  []
  [pressure_ramp]
    type = PiecewiseLinear<<<{"description": "Linearly interpolates between pairs of x-y data", "href": "../source/functions/PiecewiseLinear.html"}>>>
    x<<<{"description": "The abscissa values"}>>> = '-200 0'
    y<<<{"description": "The ordinate values"}>>> = '0 1'
  []
[]

[Kernels<<<{"href": "../syntax/Kernels/index.html"}>>>]
  [gravity]
    type = Gravity<<<{"description": "Apply gravity. Value is in units of acceleration.", "href": "../source/kernels/Gravity.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    value<<<{"description": "Value multiplied against the residual, e.g. gravitational acceleration"}>>> = -9.81
  []
  [heat]
    type = HeatConduction<<<{"description": "Diffusive heat conduction term $-\\nabla\\cdot(k\\nabla T)$ of the thermal energy conservation equation", "href": "../source/kernels/HeatConduction.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
    extra_vector_tags<<<{"description": "The extra tags for the vectors this Kernel should fill"}>>> = 'ref'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative<<<{"description": "Time derivative term $\\rho c_p \\frac{\\partial T}{\\partial t}$ of the thermal energy conservation equation.", "href": "../source/kernels/HeatConductionTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
    extra_vector_tags<<<{"description": "The extra tags for the vectors this Kernel should fill"}>>> = 'ref'
  []
  [heat_source]
    type = NeutronHeatSource<<<{"description": "Compute heat generation due to fission.", "href": "../source/kernels/NeutronHeatSource.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temp
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = pellet_type_1
    fission_rate<<<{"description": "Coupled Fission Rate"}>>> = fission_rate
    extra_vector_tags<<<{"description": "The extra tags for the vectors this Kernel should fill"}>>> = 'ref'
  []
[]

[Physics<<<{"href": "../syntax/Physics/index.html"}>>>/SolidMechanics<<<{"href": "../syntax/Physics/SolidMechanics/index.html"}>>>/QuasiStatic<<<{"href": "../syntax/Physics/SolidMechanics/QuasiStatic/index.html"}>>>]
  [pellets]
    block<<<{"description": "The list of ids of the blocks (subdomain) that the stress divergence kernels will be applied to"}>>> = pellet_type_1
    add_variables<<<{"description": "Add the displacement variables"}>>> = true
    strain<<<{"description": "Strain formulation"}>>> = FINITE
    eigenstrain_names<<<{"description": "List of eigenstrains to be applied in this strain calculation"}>>> = 'fuel_relocation_eigenstrain fuel_thermal_strain fuel_volumetric_swelling_eigenstrain'
    generate_output<<<{"description": "Add scalar quantity output for stress and/or strain"}>>> = 'vonmises_stress hydrostatic_stress stress_xx stress_yy stress_zz strain_xx strain_yy strain_zz'
    extra_vector_tags<<<{"description": "The tag names for extra vectors that residual data should be saved into"}>>> = 'ref'
  []
  [clad]
    block<<<{"description": "The list of ids of the blocks (subdomain) that the stress divergence kernels will be applied to"}>>> = clad
    add_variables<<<{"description": "Add the displacement variables"}>>> = true
    strain<<<{"description": "Strain formulation"}>>> = FINITE
    eigenstrain_names<<<{"description": "List of eigenstrains to be applied in this strain calculation"}>>> = 'clad_thermal_eigenstrain clad_irradiation_strain'
    generate_output<<<{"description": "Add scalar quantity output for stress and/or strain"}>>> = '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'
    extra_vector_tags<<<{"description": "The tag names for extra vectors that residual data should be saved into"}>>> = 'ref'
  []
[]

[Burnup<<<{"href": "../syntax/Burnup/index.html"}>>>]
  [burnup]
    block<<<{"description": "The blocks where radial power factor should be computed."}>>> = pellet_type_1
    rod_ave_lin_pow<<<{"description": "Rod average linear power function."}>>> = power_history
    axial_power_profile<<<{"description": "Axial power peaking function."}>>> = axial_peaking_factors
    num_radial<<<{"description": "Number of radial divisions in secondary grid used to compute radial power profile."}>>> = 80
    num_axial<<<{"description": "Number of axial divisions in secondary grid used to compute radial power profile."}>>> = 11
    a_lower<<<{"description": "The lower axial coordinate of the fuel stack. Required if fuel_pin_geometry is not specified."}>>> = 2.49e-3
    a_upper<<<{"description": "The upper axial coordinate of the fuel stack. Required if fuel_pin_geometry is not specified."}>>> = 2.621e-2
    fuel_inner_radius<<<{"description": "The inner radius of the fuel."}>>> = 0
    fuel_outer_radius<<<{"description": "The outer radius of the fuel."}>>> = 0.0041
    fuel_volume_ratio<<<{"description": "Reduction factor for deviation from right circular cylinder fuel.  The ratio of actual volume to right circular cylinder volume."}>>> = 1.0
    RPF<<<{"description": "Specifies that the radial power factor is required."}>>> = RPF
  []
[]

[AuxKernels<<<{"href": "../syntax/AuxKernels/index.html"}>>>]
  [fast_neutron_flux]
    type = FastNeutronFluxAux<<<{"description": "Computes fast neutron flux.", "href": "../source/auxkernels/FastNeutronFluxAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = fast_neutron_flux
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = clad
    rod_ave_lin_pow<<<{"description": "The power history function."}>>> = power_history
    axial_power_profile<<<{"description": "The axial power profile function."}>>> = axial_peaking_factors
    factor<<<{"description": "The fast neutron flux if the function, RALP, and q_variable are not given. A scaling factor if the function, RALP, or q_variable is given.  If RALP or q_variable is given, it is recommended to use a value of 3e13 (n/(m^2s)/(W/m))."}>>> = 3e13
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_begin
  []
  [fast_neutron_fluence]
    type = FastNeutronFluenceAux<<<{"description": "Computes fast neutron fluence from fast_neutron_flux.", "href": "../source/auxkernels/FastNeutronFluenceAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = fast_neutron_fluence
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = clad
    fast_neutron_flux<<<{"description": "Coupled Fast Neutron Flux."}>>> = fast_neutron_flux
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_begin
  []
  [grain_radius]
    type = GrainRadiusAux<<<{"description": "Computes grain evolution using an empirical model.", "href": "../source/auxkernels/GrainRadiusAux.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = pellet_type_1
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = grain_radius
    temperature<<<{"description": "Coupled temperature (K)"}>>> = temp
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = linear
  []
  [conductance]
    type = MaterialRealAux<<<{"description": "Outputs element volume-averaged material properties", "href": "../source/auxkernels/MaterialRealAux.html"}>>>
    property<<<{"description": "The material property name."}>>> = gap_conductance
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = gap_cond
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 10
  []
  [coolant_htc]
    type = MaterialRealAux<<<{"description": "Outputs element volume-averaged material properties", "href": "../source/auxkernels/MaterialRealAux.html"}>>>
    property<<<{"description": "The material property name."}>>> = coolant_channel_htc
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = coolant_htc
    boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = 2
  []
  [hoop_inelastic_strain]
    type = RankTwoScalarAux<<<{"description": "Compute a scalar property of a RankTwoTensor", "href": "../source/auxkernels/RankTwoScalarAux.html"}>>>
    rank_two_tensor<<<{"description": "The rank two material tensor name"}>>> = creep_strain
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = hoop_inelastic_strain
    scalar_type<<<{"description": "Type of scalar output"}>>> = HoopStress
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = timestep_end
  []
[]

[Contact<<<{"href": "../syntax/Contact/index.html"}>>>]
  [pellet_clad_mechanical]
    primary<<<{"description": "The list of boundary IDs referring to primary sidesets"}>>> = 5
    secondary<<<{"description": "The list of boundary IDs referring to secondary sidesets"}>>> = 10
    formulation<<<{"description": "The contact formulation"}>>> = kinematic
    penalty<<<{"description": "The penalty to apply.  This can vary depending on the stiffness of your materials"}>>> = 1e14
    normalize_penalty<<<{"description": "Whether to normalize the penalty parameter with the nodal area."}>>> = true
    model<<<{"description": "The contact model to use"}>>> = frictionless
    normal_smoothing_distance<<<{"description": "Distance from edge in parametric coordinates over which to smooth contact normal"}>>> = 0.1
  []
[]

[ThermalContact<<<{"href": "../syntax/Modules/HeatTransfer/ThermalContact/index.html"}>>>]
  [thermal_contact]
    type = GasGapHeatTransfer
    variable = temp
    primary = 5
    secondary = 10
    initial_moles = initial_moles
    gas_released = fission_gas_released
    tangential_tolerance = 1e-4
    contact_pressure = contact_pressure
    quadrature = true
    normal_smoothing_distance = 0.1
  []
[]

[BCs<<<{"href": "../syntax/NuclearMaterials/BCs/index.html"}>>>]
  [no_x_all]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 12
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [no_z_all]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_z
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 13
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [no_y_clad_bottom]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 1
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [no_y_fuel_bottom]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 1020
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [Pressure<<<{"href": "../syntax/BCs/Pressure/index.html"}>>>]
    [coolantPressure]
      boundary<<<{"description": "The list of boundaries (ids or names) from the mesh where this object applies"}>>> = '1 2 3'
      factor<<<{"description": "The magnitude to use in computing the pressure"}>>> = 15.5e6
      function<<<{"description": "The function that describes the pressure"}>>> = pressure_ramp
    []
  []
  [PlenumPressure<<<{"href": "../syntax/BCs/PlenumPressure/index.html"}>>>]
    [plenumPressure]
      boundary<<<{"description": "The list of boundary IDs from the mesh where the pressure will be applied"}>>> = 9
      initial_pressure<<<{"description": "The initial pressure in the cavity.  If not given, a zero initial pressure will be used."}>>> = 2.0e6
      startup_time<<<{"description": "The amount of time during which the pressure will ramp from zero to its true value."}>>> = 0
      R<<<{"description": "The universal gas constant for the units used."}>>> = 8.3143
      output_initial_moles<<<{"description": "The name to use when reporting the initial moles of gas"}>>> = initial_moles
      temperature<<<{"description": "The name of the average temperature postprocessor value."}>>> = interior_temp
      volume<<<{"description": "The name of the postprocessor(s) that holds the value of the internal volume in the cavity"}>>> = gas_volume
      material_input<<<{"description": "The name of the postprocessor(s) that holds the amount of material injected into the plenum."}>>> = fission_gas_released
      output<<<{"description": "The name to use for the cavity pressure value"}>>> = plenum_pressure
    []
  []
[]

[CoolantChannel<<<{"href": "../syntax/CoolantChannel/index.html"}>>>]
  [convective_clad_surface]
    boundary<<<{"description": "The boundary where the coolant channel calculation will run"}>>> = '1 2 3'
    variable<<<{"description": "The name of the variable representing temperature"}>>> = temp
    inlet_temperature<<<{"description": "Inlet temperature in K "}>>> = 580 # K
    inlet_pressure<<<{"description": "Inlet pressure in Pa"}>>> = 15.5e6 # Pa
    inlet_massflux<<<{"description": "Inlet mass flux in kg/m^2-sec"}>>> = 3800 # kg/m^2-sec
    rod_diameter<<<{"description": "Rod diameter in meter"}>>> = 0.948e-2 # m
    rod_pitch<<<{"description": "Rod pitch in meter"}>>> = 1.26e-2 # m
    linear_heat_rate<<<{"description": "Linear heat generation rate in W/m"}>>> = power_history
    axial_power_profile<<<{"description": "Axial power profile"}>>> = axial_peaking_factors
  []
[]

[Materials<<<{"href": "../syntax/Materials/index.html"}>>>]
  [fuel_thermal]
    type = UO2Thermal<<<{"description": "Computes thermal conductivity and specific heat capacity of uranium dioxide fuel.", "href": "../source/materials/UO2Thermal.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = pellet_type_1
    thermal_conductivity_model<<<{"description": "The thermal conductivity model."}>>> = NFIR
    temperature<<<{"description": "Coupled temperature"}>>> = temp
    burnup_function<<<{"description": "Burnup function"}>>> = burnup
    initial_porosity<<<{"description": "Initial porosity.  Must be between 0.0 and 1.0."}>>> = 0.05
  []
  [fuel_elastic_stress]
    type = ComputeFiniteStrainElasticStress<<<{"description": "Compute stress using elasticity for finite strains", "href": "../source/materials/ComputeFiniteStrainElasticStress.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = pellet_type_1
  []
  [fuel_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor<<<{"description": "Compute a constant isotropic elasticity tensor.", "href": "../source/materials/ComputeIsotropicElasticityTensor.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = pellet_type_1
    youngs_modulus<<<{"description": "Young's modulus of the material."}>>> = 2.0e11
    poissons_ratio<<<{"description": "Poisson's ratio for the material."}>>> = 0.345
  []
  [fuel_relocation]
    type = UO2RelocationEigenstrain<<<{"description": "Accounts for cracking and relocation of fuel pellet fragments in the radial direction and is necessary for accurate modeling of LWR fuel. The linear heat rate must either be a functionor a variable.", "href": "../source/materials/solid_mechanics/UO2RelocationEigenstrain.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = pellet_type_1
    burnup_function<<<{"description": "Burnup function"}>>> = burnup
    diameter<<<{"description": "As fabricated cold diameter of pellet in meters"}>>> = 0.0082
    rod_ave_lin_pow<<<{"description": "linear power function."}>>> = power_history
    axial_power_profile<<<{"description": "axial power peaking function."}>>> = axial_peaking_factors
    diametral_gap<<<{"description": "As fabricated cold diametral pellet-cladding gap in meters."}>>> =160e-6
    relocation_activation1<<<{"description": "First activation linear power in W/m term in relocation model"}>>> = 5000
    burnup_relocation_stop<<<{"description": "Burnup at which relocation strain stops (in FIMA)"}>>> = 0.02
    eigenstrain_name<<<{"description": "Material property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator."}>>> = 'fuel_relocation_eigenstrain'
  []
  [fuel_swelling]
    type = UO2VolumetricSwellingEigenstrain<<<{"description": "Computes and sums the change in fuel pellet volume due to densification and fission product release. This class applies a volumetric strain correction before adding the strain from this class to the diagonal entries of the eigenstrain tensor.", "href": "../source/materials/solid_mechanics/UO2VolumetricSwellingEigenstrain.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = pellet_type_1
    burnup_function<<<{"description": "Burnup function"}>>> = burnup
    temperature<<<{"description": "Coupled Temperature in Kelvin"}>>> = temp
    initial_fuel_density<<<{"description": "Initial fuel density in kg-UO2/m^3"}>>> = 10431.0
    eigenstrain_name<<<{"description": "Material property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator."}>>> = 'fuel_volumetric_swelling_eigenstrain'
  []
  [fuel_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain<<<{"description": "Computes eigenstrain due to thermal expansion with a constant coefficient", "href": "../source/materials/ComputeThermalExpansionEigenstrain.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = pellet_type_1
    thermal_expansion_coeff<<<{"description": "Thermal expansion coefficient"}>>> = 10.0e-6
    temperature<<<{"description": "Coupled temperature"}>>> = temp
    stress_free_temperature<<<{"description": "Reference temperature at which there is no thermal expansion for thermal eigenstrain calculation"}>>> = 580.0
    eigenstrain_name<<<{"description": "Material property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator."}>>> = 'fuel_thermal_strain'
  []
  [fission_gas_release]
    type = UO2Sifgrs<<<{"description": "Recommended fission gas model to account for generation of fission gasses in nuclear fuel", "href": "../source/materials/UO2Sifgrs.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = pellet_type_1
    temperature<<<{"description": "Coupled temperature"}>>> = temp
    fission_rate<<<{"description": "Coupled fission rate variable (fiss/m^3/s)"}>>> = fission_rate
    grain_radius<<<{"description": "Coupled grain Radius"}>>> = grain_radius
    gbs_model<<<{"description": "Activates the grain-boundary sweeping model"}>>> = true
  []
  [fuel_density]
    type = StrainAdjustedDensity<<<{"description": "Creates density material property", "href": "../source/materials/StrainAdjustedDensity.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = pellet_type_1
    strain_free_density<<<{"description": "Material property for strain-free density"}>>> = ${initial_fuel_density}
  []

  [clad_thermal]
    type = HeatConductionMaterial<<<{"description": "General-purpose material model for heat conduction", "href": "../source/materials/HeatConductionMaterial.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = clad
    thermal_conductivity<<<{"description": "The thermal conductivity value"}>>> = 16.0
    specific_heat<<<{"description": "The specific heat value"}>>> = 330.0
  []
  [clad_elasticity_tensor]
    type = ZryElasticityTensor<<<{"description": "Either provides constant elasticity constants for Zircaloy cladding or calculates the Young's modulus and Poisson's ratio for Zircaloy cladding using MATPRO relations as a function of temperature and fast neutron fluence.", "href": "../source/materials/solid_mechanics/ZryElasticityTensor.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = clad
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress<<<{"description": "Compute state (stress and internal parameters such as plastic strains and internal parameters) using an iterative process.  Combinations of creep models and plastic models may be used.", "href": "../source/materials/ComputeMultipleInelasticStress.html"}>>>
    tangent_operator<<<{"description": "Type of tangent operator to return.  'elastic': return the elasticity tensor.  'nonlinear': return the full, general consistent tangent operator."}>>> = elastic
    inelastic_models<<<{"description": "The material objects to use to calculate stress and inelastic strains. Note: specify creep models first and plasticity models second."}>>> = 'clad_zrycreep'
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = clad
  []
  [clad_zrycreep]
    type = ZryCreepLimbackHoppeUpdate<<<{"description": "Computes the Limback-Andersson thermal primary and secondary creep and the Hoppe irradiation creep for Zircaloy cladding.  This material must be run in conjunction with ComputeMultipleInelasticStress.", "href": "../source/materials/solid_mechanics/ZryCreepLimbackHoppeUpdate.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = clad
    temperature<<<{"description": "The coupled temperature (K)"}>>> = temp
    fast_neutron_flux<<<{"description": "The fast neutron flux"}>>> = fast_neutron_flux
    fast_neutron_fluence<<<{"description": "The fast neutron fluence"}>>> = fast_neutron_fluence
    model_irradiation_creep<<<{"description": "Set true to activate irradiation induced creep"}>>> = true
    model_primary_creep<<<{"description": "Set true to activate primary creep"}>>> = true
    model_thermal_creep<<<{"description": "Set true to activate steady state thermal creep"}>>> = true
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain<<<{"description": "Computes eigenstrain due to anisotropic thermal expansion in Zircaloy cladding using Matpro correlations.", "href": "../source/materials/solid_mechanics/ZryThermalExpansionMATPROEigenstrain.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = clad
    temperature<<<{"description": "Coupled temperature"}>>> = temp
    stress_free_temperature<<<{"description": "Reference temperature at which there is no thermal expansion for thermal eigenstrain calculation"}>>> = 580.0
    eigenstrain_name<<<{"description": "Material property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator."}>>> = clad_thermal_eigenstrain
  []
  [irradiation_swelling]
    type = ZryIrradiationGrowthEigenstrain<<<{"description": "Computes eigenstrain from irradiation growth in Zircaloy cladding using either the Franklin or ESCORE models.", "href": "../source/materials/solid_mechanics/ZryIrradiationGrowthEigenstrain.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = clad
    fast_neutron_fluence<<<{"description": "The fast neutron fluence"}>>> = fast_neutron_fluence
    zircaloy_material_type<<<{"description": "Type of irradiation growth volumetric swelling formulation."}>>> = stress_relief_annealed
    eigenstrain_name<<<{"description": "Material property name for the eigenstrain tensor computed by this model. IMPORTANT: The name of this property must also be provided to the strain calculator."}>>> = clad_irradiation_strain
  []
  [clad_density]
    type = StrainAdjustedDensity<<<{"description": "Creates density material property", "href": "../source/materials/StrainAdjustedDensity.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = clad
    strain_free_density<<<{"description": "Material property for strain-free density"}>>> = 6551.0
  []
[]

[Preconditioning<<<{"href": "../syntax/Preconditioning/index.html"}>>>]
  [SMP]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
  []
[]

[Executioner<<<{"href": "../syntax/Executioner/index.html"}>>>]

  type = Transient
  solve_type = 'PJFNK'

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

  line_search = 'none'

  l_max_its = 50
  l_tol = 8e-3

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

  start_time = -200
  end_time = 3.0e7

  dtmax = 2e6
  dtmin = 1

  [TimeStepper<<<{"href": "../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = IterationAdaptiveDT
    dt = 200
    optimal_iterations = 15
    iteration_window = 3
    linear_iteration_ratio = 100
  []

  [Quadrature<<<{"href": "../syntax/Executioner/Quadrature/index.html"}>>>]
    order<<<{"description": "Order of the quadrature"}>>> = fifth
    side_order<<<{"description": "Order of the quadrature for sides"}>>> = seventh
  []
[]

[Postprocessors<<<{"href": "../syntax/Postprocessors/index.html"}>>>]
  [interior_temp]
    type = SideAverageValue<<<{"description": "Computes the average value of a variable on a sideset. Note that this cannot be used on the centerline of an axisymmetric model.", "href": "../source/postprocessors/SideAverageValue.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 7
    variable<<<{"description": "The name of the variable which this postprocessor integrates"}>>> = temp
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'initial linear'
  []
  [clad_inner_vol]
    type = InternalVolume<<<{"description": "Computes the volume of an enclosed area by performing an integral over a user-supplied boundary.", "href": "../source/postprocessors/InternalVolume.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 7
    outputs<<<{"description": "Vector of output names where you would like to restrict the output of variables(s) associated with this object"}>>> = exodus
  []
  [pellet_volume]
    type = InternalVolume<<<{"description": "Computes the volume of an enclosed area by performing an integral over a user-supplied boundary.", "href": "../source/postprocessors/InternalVolume.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 8
    outputs<<<{"description": "Vector of output names where you would like to restrict the output of variables(s) associated with this object"}>>> = exodus
  []
  [avg_clad_temp]
    type = SideAverageValue<<<{"description": "Computes the average value of a variable on a sideset. Note that this cannot be used on the centerline of an axisymmetric model.", "href": "../source/postprocessors/SideAverageValue.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 7
    variable<<<{"description": "The name of the variable which this postprocessor integrates"}>>> = temp
  []
  [fis_gas_produced]
    type = ElementIntegralFisGasGeneratedSifgrs<<<{"description": "Reports the fission gas that is produced in moles.  To be used in combination with the Sifgrs model.", "href": "../source/postprocessors/ElementIntegralFisGasGeneratedSifgrs.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = pellet_type_1
  []
  [fission_gas_released]
    type = ElementIntegralFisGasReleasedSifgrs<<<{"description": "Reports the fission gas that is released to the plenum in moles.  To be used in combination with the Sifgrs model.", "href": "../source/postprocessors/ElementIntegralFisGasReleasedSifgrs.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = pellet_type_1
  []
  [fis_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs<<<{"description": "Reports the fission gas within the grains in moles.  To be used in combination with the Sifgrs model.", "href": "../source/postprocessors/ElementIntegralFisGasGrainSifgrs.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = pellet_type_1
    outputs<<<{"description": "Vector of output names where you would like to restrict the output of variables(s) associated with this object"}>>> = exodus
  []
  [fis_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs<<<{"description": "Reports the fission gas that is on the grain boundary", "href": "../source/postprocessors/ElementIntegralFisGasBoundarySifgrs.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = pellet_type_1
    outputs<<<{"description": "Vector of output names where you would like to restrict the output of variables(s) associated with this object"}>>> = exodus
  []
  [gas_volume]
    type = InternalVolume<<<{"description": "Computes the volume of an enclosed area by performing an integral over a user-supplied boundary.", "href": "../source/postprocessors/InternalVolume.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 9
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'initial linear'
  []
  [flux_from_clad]
    type = SideDiffusiveFluxIntegral<<<{"description": "Computes the integral of the diffusive flux over the specified boundary", "href": "../source/postprocessors/SideDiffusiveFluxIntegral.html"}>>>
    variable<<<{"description": "The name of the variable which this postprocessor integrates"}>>> = temp
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 5
    diffusivity<<<{"description": "The name of the diffusivity material property that will be used in the flux computation. This must be provided if the variable is of finite element type"}>>> = thermal_conductivity
  []
  [flux_from_fuel]
    type = SideDiffusiveFluxIntegral<<<{"description": "Computes the integral of the diffusive flux over the specified boundary", "href": "../source/postprocessors/SideDiffusiveFluxIntegral.html"}>>>
    variable<<<{"description": "The name of the variable which this postprocessor integrates"}>>> = temp
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 10
    diffusivity<<<{"description": "The name of the diffusivity material property that will be used in the flux computation. This must be provided if the variable is of finite element type"}>>> = thermal_conductivity
  []
  [average_fissionrate]
    type = ElementAverageValue<<<{"description": "Computes the volumetric average of a variable", "href": "../source/postprocessors/ElementAverageValue.html"}>>>
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = 3
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = fission_rate
  []
  [rod_total_power] # should be 1/4 of the rod_input_power as we are using in quarter symmetry
    type = ElementIntegralPower<<<{"description": "Computes the power given the fission rate and energy per fission.", "href": "../source/postprocessors/ElementIntegralPower.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = temp
    fission_rate<<<{"description": "Coupled fission rate"}>>> = fission_rate
    block<<<{"description": "The list of blocks (ids or names) that this object will be applied"}>>> = pellet_type_1
  []
  [rod_input_power]
    type = FunctionValuePostprocessor<<<{"description": "Computes the value of a supplied function at a single point (scalable)", "href": "../source/postprocessors/FunctionValuePostprocessor.html"}>>>
    function<<<{"description": "The function which supplies the postprocessor value."}>>> = power_history
    scale_factor<<<{"description": "A scale factor to be applied to the function"}>>> = 0.02372
  []
  [average_fission_rate]
    type = AverageFissionRate<<<{"description": "Computes the average fission rate (average over time and space).", "href": "../source/postprocessors/AverageFissionRate.html"}>>>
    rod_ave_lin_pow<<<{"description": "The power function"}>>> = power_history
  []
[]

[PerformanceMetricOutputs<<<{"href": "../syntax/PerformanceMetricOutputs/index.html"}>>>]
[]

[Outputs<<<{"href": "../syntax/Outputs/index.html"}>>>]
  perf_graph<<<{"description": "Enable printing of the performance graph to the screen (Console)"}>>> = true
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
  [console]
    type = Console<<<{"description": "Object for screen output.", "href": "../source/outputs/Console.html"}>>>
    max_rows<<<{"description": "The maximum number of postprocessor/scalar values displayed on screen during a timestep (set to 0 for unlimited)"}>>> = 25
  []
  [chkfile]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../source/outputs/CSV.html"}>>>
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = FINAL
    show<<<{"description": "A list of the variables and postprocessors that should be output to the Exodus file (may include Variables, ScalarVariables, and Postprocessor names)."}>>> = 'fission_gas_released plenum_pressure interior_temp gas_volume'
  []
[]

0,  1e4,    1e8
0, 2.5e4, 2.5e4

Mortar input files

An alternative and more recent approach to the thermal and mechanical interface problem is that based on the mortar finite element method. In this approach, the two discretized surfaces create an additional intermediate lower-dimensional mesh on which constraints and residuals are resolved. This approach tends to converge more robustly, provides smoother primal variable output on contacting surfaces, and avoids the need for overintegrating finite element surfaces. Steps to convert an NTS input file to the mortar approach include: The use of thermal LWR and contact mechanics actions, the specification of mortar executioner's solver settings, removal of the quadrature block in the executioner, and block restriction of kernels when such block(s) had not been specified in the input file block. Note that the mortar framework creates lower-dimensional blocks on which typical Bison kernels should not act.

Below are excerpts of actions using the mortar approach to the thermomechanical problem.

Contact mechanics action

[Contact<<<{"href": "../syntax/Contact/index.html"}>>>]
  [pellet_clad_mechanical]
    primary<<<{"description": "The list of boundary IDs referring to primary sidesets"}>>> = 5
    secondary<<<{"description": "The list of boundary IDs referring to secondary sidesets"}>>> = 10
    formulation<<<{"description": "The contact formulation"}>>> = mortar
    model<<<{"description": "The contact model to use"}>>> = coulomb
    c_normal<<<{"description": "Parameter for balancing the size of the gap and contact pressure for a mortar formulation. This purely numerical parameter affects convergence behavior and, in general, should be larger for stiffer materials. It is recommended that the user tries out various orders of magnitude for this parameter if the default value generates poor contact convergence."}>>> = 1e+18
    c_tangential<<<{"description": "Numerical parameter for nonlinear mortar frictional constraints"}>>> = 1e+18
    friction_coefficient<<<{"description": "The friction coefficient"}>>> = 0.5
  []
[]

As guidance for users to set up frictional problems, these 3D mortar thermomechanical examples include friction on the pellet-cladding interface. Note that other options in the action, such as Lagrange multiplier scaling and treatment for edge-to-edge contact are also available.

LWR gap heat transfer action

[ThermalContactMortar<<<{"href": "../syntax/ThermalContactMortar/index.html"}>>>]
  [thermal_contact]
    secondary_variable<<<{"description": "The secondary variable"}>>> = temp
    primary_boundary<<<{"description": "The primary surface"}>>> = '5'
    secondary_boundary<<<{"description": "The secondary surface"}>>> = '10'
    gas_released<<<{"description": "The postprocessor(s) that will give the gas released."}>>> = fis_gas_released_model
    initial_moles<<<{"description": "The Postprocessor that will give the initial moles of gas."}>>> = initial_moles
    jump_distance_model<<<{"description": "Jump distance model: DIRECT=specify distances directly;LANNING=computed as a function of gas properties (Lanning and Hahn, 1975)TOPTAN=computed as a function of gas properties (Toptan et al., 2019)"}>>> = LANNING
    plenum_pressure<<<{"description": "The name of the plenum_pressure postprocessor value."}>>> = plenum_pressure
    roughness_coef<<<{"description": "The coefficient for the roughness summation."}>>> = 3.2
    roughness_secondary<<<{"description": "The roughness of the secondary surface (m)."}>>> = 1e-6
    roughness_primary<<<{"description": "The roughness of the primary surface (m)."}>>> = 2e-6
    emissivity_primary<<<{"description": "The emissivity of the fuel surface"}>>> = 0.8
    emissivity_secondary<<<{"description": "The emissivity of the cladding surface"}>>> = 0.8
  []
[]

Note that the 3D-pellet mortar examples have an artificially high fuel thermal expansion coefficient to exhibit mechanical contact early in the simulation. This is done solely to speed up the simulation into mechanical contact, which is a typical form of simulation failure. Users can manually modify the input file to restore a realistic fuel thermal expansion coefficient and extend the simulation end_time.

Figure 1: Pellet-cladding normal contact pressure in a discrete pellet simulation with mortar contact.

(examples/3D_rodlet_3pellets/smeared/smearedTest3D.i)


initial_fuel_density = 10431.0

[GlobalParams]
  density = ${initial_fuel_density} # initial fuel density 95.0% of theoretical (10980 kg/m3)
  displacements = 'disp_x disp_y disp_z'
  order = SECOND
  family = LAGRANGE
  energy_per_fission = 3.2e-11 # J/fission
[]

[Mesh]
  patch_size = 20
  patch_update_strategy = iteration
  partitioner = centroid
  centroid_partitioner_direction = y
  [mesh]
    type = FileMeshGenerator
    file = smearedTest3.e
  []
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
  group_variables = 'disp_x disp_y disp_z'
[]

[Variables]
  [temp]
    initial_condition = 580
  []
[]

[AuxVariables]
  [fast_neutron_flux]
    block = clad
  []
  [fast_neutron_fluence]
    block = clad
  []
  [grain_radius]
    block = pellet_type_1
    initial_condition = 10e-6
  []
  [gap_cond]
    order = CONSTANT
    family = MONOMIAL
  []
  [coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  []
  [hoop_inelastic_strain]
    order = CONSTANT
    family = MONOMIAL
    block = clad
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    data_file = powerhistory.csv
    scale_factor = 1
  []
  [axial_peaking_factors]
    type = ParsedFunction
    expression = 1
  []
  [pressure_ramp]
    type = PiecewiseLinear
    x = '-200 0'
    y = '0 1'
  []
[]

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

[Physics/SolidMechanics/QuasiStatic]
  [pellets]
    block = pellet_type_1
    add_variables = true
    strain = FINITE
    eigenstrain_names = 'fuel_relocation_eigenstrain fuel_thermal_strain fuel_volumetric_swelling_eigenstrain'
    generate_output = 'vonmises_stress hydrostatic_stress stress_xx stress_yy stress_zz strain_xx strain_yy strain_zz'
    extra_vector_tags = 'ref'
  []
  [clad]
    block = clad
    add_variables = true
    strain = FINITE
    eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_strain'
    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'
    extra_vector_tags = 'ref'
  []
[]

[Burnup]
  [burnup]
    block = pellet_type_1
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    num_radial = 80
    num_axial = 11
    a_lower = 2.49e-3
    a_upper = 2.621e-2
    fuel_inner_radius = 0
    fuel_outer_radius = 0.0041
    fuel_volume_ratio = 1.0
    RPF = RPF
  []
[]

[AuxKernels]
  [fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    block = clad
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    factor = 3e13
    execute_on = timestep_begin
  []
  [fast_neutron_fluence]
    type = FastNeutronFluenceAux
    variable = fast_neutron_fluence
    block = clad
    fast_neutron_flux = fast_neutron_flux
    execute_on = timestep_begin
  []
  [grain_radius]
    type = GrainRadiusAux
    block = pellet_type_1
    variable = grain_radius
    temperature = temp
    execute_on = linear
  []
  [conductance]
    type = MaterialRealAux
    property = gap_conductance
    variable = gap_cond
    boundary = 10
  []
  [coolant_htc]
    type = MaterialRealAux
    property = coolant_channel_htc
    variable = coolant_htc
    boundary = 2
  []
  [hoop_inelastic_strain]
    type = RankTwoScalarAux
    rank_two_tensor = creep_strain
    variable = hoop_inelastic_strain
    scalar_type = HoopStress
    execute_on = timestep_end
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    formulation = kinematic
    penalty = 1e14
    normalize_penalty = true
    model = frictionless
    normal_smoothing_distance = 0.1
  []
[]

[ThermalContact]
  [thermal_contact]
    type = GasGapHeatTransfer
    variable = temp
    primary = 5
    secondary = 10
    initial_moles = initial_moles
    gas_released = fission_gas_released
    tangential_tolerance = 1e-4
    contact_pressure = contact_pressure
    quadrature = true
    normal_smoothing_distance = 0.1
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 12
    value = 0.0
  []
  [no_z_all]
    type = DirichletBC
    variable = disp_z
    boundary = 13
    value = 0.0
  []
  [no_y_clad_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 1
    value = 0.0
  []
  [no_y_fuel_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 1020
    value = 0.0
  []
  [Pressure]
    [coolantPressure]
      boundary = '1 2 3'
      factor = 15.5e6
      function = pressure_ramp
    []
  []
  [PlenumPressure]
    [plenumPressure]
      boundary = 9
      initial_pressure = 2.0e6
      startup_time = 0
      R = 8.3143
      output_initial_moles = initial_moles
      temperature = interior_temp
      volume = gas_volume
      material_input = fission_gas_released
      output = plenum_pressure
    []
  []
[]

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

[Materials]
  [fuel_thermal]
    type = UO2Thermal
    block = pellet_type_1
    thermal_conductivity_model = NFIR
    temperature = temp
    burnup_function = burnup
    initial_porosity = 0.05
  []
  [fuel_elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = pellet_type_1
  []
  [fuel_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = pellet_type_1
    youngs_modulus = 2.0e11
    poissons_ratio = 0.345
  []
  [fuel_relocation]
    type = UO2RelocationEigenstrain
    block = pellet_type_1
    burnup_function = burnup
    diameter = 0.0082
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    diametral_gap =160e-6
    relocation_activation1 = 5000
    burnup_relocation_stop = 0.02
    eigenstrain_name = 'fuel_relocation_eigenstrain'
  []
  [fuel_swelling]
    type = UO2VolumetricSwellingEigenstrain
    block = pellet_type_1
    burnup_function = burnup
    temperature = temp
    initial_fuel_density = 10431.0
    eigenstrain_name = 'fuel_volumetric_swelling_eigenstrain'
  []
  [fuel_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = pellet_type_1
    thermal_expansion_coeff = 10.0e-6
    temperature = temp
    stress_free_temperature = 580.0
    eigenstrain_name = 'fuel_thermal_strain'
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = pellet_type_1
    temperature = temp
    fission_rate = fission_rate
    grain_radius = grain_radius
    gbs_model = true
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = pellet_type_1
    strain_free_density = ${initial_fuel_density}
  []

  [clad_thermal]
    type = HeatConductionMaterial
    block = clad
    thermal_conductivity = 16.0
    specific_heat = 330.0
  []
  [clad_elasticity_tensor]
    type = ZryElasticityTensor
    block = clad
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep'
    block = clad
  []
  [clad_zrycreep]
    type = ZryCreepLimbackHoppeUpdate
    block = clad
    temperature = temp
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    model_irradiation_creep = true
    model_primary_creep = true
    model_thermal_creep = true
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = clad
    temperature = temp
    stress_free_temperature = 580.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [irradiation_swelling]
    type = ZryIrradiationGrowthEigenstrain
    block = clad
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = stress_relief_annealed
    eigenstrain_name = clad_irradiation_strain
  []
  [clad_density]
    type = StrainAdjustedDensity
    block = clad
    strain_free_density = 6551.0
  []
[]

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

[Executioner]

  type = Transient
  solve_type = 'PJFNK'

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

  line_search = 'none'

  l_max_its = 50
  l_tol = 8e-3

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

  start_time = -200
  end_time = 3.0e7

  dtmax = 2e6
  dtmin = 1

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 200
    optimal_iterations = 15
    iteration_window = 3
    linear_iteration_ratio = 100
  []

  [Quadrature]
    order = fifth
    side_order = seventh
  []
[]

[Postprocessors]
  [interior_temp]
    type = SideAverageValue
    boundary = 7
    variable = temp
    execute_on = 'initial linear'
  []
  [clad_inner_vol]
    type = InternalVolume
    boundary = 7
    outputs = exodus
  []
  [pellet_volume]
    type = InternalVolume
    boundary = 8
    outputs = exodus
  []
  [avg_clad_temp]
    type = SideAverageValue
    boundary = 7
    variable = temp
  []
  [fis_gas_produced]
    type = ElementIntegralFisGasGeneratedSifgrs
    block = pellet_type_1
  []
  [fission_gas_released]
    type = ElementIntegralFisGasReleasedSifgrs
    block = pellet_type_1
  []
  [fis_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    block = pellet_type_1
    outputs = exodus
  []
  [fis_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    block = pellet_type_1
    outputs = exodus
  []
  [gas_volume]
    type = InternalVolume
    boundary = 9
    execute_on = 'initial linear'
  []
  [flux_from_clad]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 5
    diffusivity = thermal_conductivity
  []
  [flux_from_fuel]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 10
    diffusivity = thermal_conductivity
  []
  [average_fissionrate]
    type = ElementAverageValue
    block = 3
    variable = fission_rate
  []
  [rod_total_power] # should be 1/4 of the rod_input_power as we are using in quarter symmetry
    type = ElementIntegralPower
    variable = temp
    fission_rate = fission_rate
    block = pellet_type_1
  []
  [rod_input_power]
    type = FunctionValuePostprocessor
    function = power_history
    scale_factor = 0.02372
  []
  [average_fission_rate]
    type = AverageFissionRate
    rod_ave_lin_pow = power_history
  []
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  exodus = true
  [console]
    type = Console
    max_rows = 25
  []
  [chkfile]
    type = CSV
    execute_on = FINAL
    show = 'fission_gas_released plenum_pressure interior_temp gas_volume'
  []
[]

(examples/3D_rodlet_3pellets/discrete_full/3d_3pellets_mortar.i)


initial_fuel_density = 10431.0

[GlobalParams]
  density = ${initial_fuel_density} #95% TD (TD = 10980)
  displacements = 'disp_x disp_y disp_z'
  order = FIRST
  family = LAGRANGE
  energy_per_fission = 3.2e-11 # J/fission (205 Mev)
  volumetric_locking_correction = true
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
  group_variables = 'disp_x disp_y disp_z'
  converge_on = 'disp_x disp_y disp_z temp'
[]

[Mesh]
  [file]
    type = FileMeshGenerator
    file = DiscreteThreePellets3D_full_HEX8.e
  []
  partitioner = centroid
  centroid_partitioner_direction = y
  patch_size = 5
  patch_update_strategy = auto
[]

[Variables]
  [temp]
    initial_condition = 580.0
  []
[]

[AuxVariables]
  [fission_rate]
    block = 3
  []
  [fast_neutron_flux]
  []
  [fast_neutron_fluence]
  []
  [grain_radius]
    block = 3
    initial_condition = 5e-6
  []
  [frictional_status]
    family = LAGRANGE
    order = FIRST
  []
  [gas_swell]
    order = CONSTANT
    family = MONOMIAL
    block = 3
  []
  [hoop_inelastic_strain]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [gap_conductance]
    order = CONSTANT
    family = MONOMIAL
  []
  [pid]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 1.0e4 1.0e8'
    y = '0  1.0  1.0'
    scale_factor = 25e3 # 25 kW/m peak power.
  []
  [axial_peaking_factors]
    type = ParsedFunction
    expression = 1
  []
  [coolant_pressure_ramp]
    type = PiecewiseLinear
    x = '-200 0'
    y = '0 1'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [pellets]
    block = 3
    add_variables = true
    strain = FINITE
    eigenstrain_names = 'fuel_relocation_eigenstrain fuel_thermal_strain fuel_volumetric_swelling_eigenstrain'
    generate_output = 'vonmises_stress hydrostatic_stress stress_xx stress_yy stress_zz strain_xx strain_yy strain_zz'
    extra_vector_tags = 'ref'
    temperature = temp
  []
  [clad]
    block = 1
    add_variables = true
    strain = FINITE
    eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_strain'
    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'
    extra_vector_tags = 'ref'
    temperature = temp
  []
[]

[Burnup]
  [burnup]
    block = 3
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    num_radial = 80
    num_axial = 11
    a_lower = 2.49e-3
    a_upper = 2.621e-2
    fuel_inner_radius = 0
    fuel_outer_radius = 0.0041
    fuel_volume_ratio = 1.0
    RPF = RPF
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temp
    extra_vector_tags = 'ref'
    block = '1 3'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temp
    extra_vector_tags = 'ref'
    block = '1 3'
  []
  [heat_source_fuel]
    type = NeutronHeatSource
    variable = temp
    block = 3
    fission_rate = fission_rate
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [fissionrate]
    type = FissionRateGeneral
    fission_rate_formulation = GENERIC
    variable = fission_rate
    block = 3
    value = 5.3548e+14
    fission_rate_function = power_history
  []
  [frictional_state]
    type = MortarFrictionalStateAux
    tangent_one = pellet_clad_mechanical_tangential_lm
    tangent_two = pellet_clad_mechanical_tangential_3d_lm
    boundary = 10
    contact_pressure = pellet_clad_mechanical_normal_lm
    variable = frictional_status
    mu = 0.5
  []
  [fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    block = 1
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    factor = 2.34e+13
    execute_on = timestep_begin
  []
  [fast_neutron_fluence]
    type = FastNeutronFluenceAux
    variable = fast_neutron_fluence
    block = 1
    fast_neutron_flux = fast_neutron_flux
    execute_on = timestep_begin
  []
  [grain_radius]
    type = GrainRadiusAux
    block = 3
    variable = grain_radius
    temperature = temp
    execute_on = linear
  []
  [gas_swell]
    type = MaterialRealAux
    block = 3
    variable = gas_swell
    property = deltav_v0_bubble_GB
    execute_on = timestep_end
  []
  [hoop_inelastic_strain]
    type = RankTwoScalarAux
    rank_two_tensor = creep_strain
    variable = hoop_inelastic_strain
    scalar_type = HoopStress
    execute_on = timestep_end
  []
  [pid]
    type = ProcessorIDAux
    variable = pid
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    formulation = mortar
    model = coulomb
    c_normal = 1e+18
    c_tangential = 1e+18
    friction_coefficient = 0.5
  []
[]

[ThermalContactMortar]
  [thermal_contact]
    secondary_variable = temp
    primary_boundary = '5'
    secondary_boundary = '10'
    gas_released = fis_gas_released_model
    initial_moles = initial_moles
    jump_distance_model = LANNING
    plenum_pressure = plenum_pressure
    roughness_coef = 3.2
    roughness_secondary = 1e-6
    roughness_primary = 2e-6
    emissivity_primary = 0.8
    emissivity_secondary = 0.8
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 1004
    value = 0.0
  []
  [no_z_all]
    type = DirichletBC
    variable = disp_z
    boundary = 1004
    value = 0.0
  []
  [no_y_clad_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 1
    value = 0.0
  []
  [no_x_clad_bottom]
    type = DirichletBC
    variable = disp_x
    boundary = 1
    value = 0.0
  []
  [no_z_clad_bottom]
    type = DirichletBC
    variable = disp_z
    boundary = 1
    value = 0.0
  []
  [no_z_fuel_bottom_point]
    type = DirichletBC
    variable = disp_z
    boundary = 1110
    value = 0.0
  []
  [no_z_fuel_bottom_point_y]
    type = DirichletBC
    variable = disp_y
    boundary = 1110
    value = 0.0
  []
  [no_x_fuel_bottom_point]
    type = DirichletBC
    variable = disp_x
    boundary = 1120
    value = 0.0
  []
  [no_x_fuel_bottom_point_y]
    type = DirichletBC
    variable = disp_y
    boundary = 1120
    value = 0.0
  []
  [Pressure]
    [coolantPressure]
      boundary = 2
      factor = 15.5e6
      function = coolant_pressure_ramp
    []
  []

  [PlenumPressure]
    [plenumPressure]
      boundary = 9
      initial_pressure = 2.0e6
      startup_time = 0.0
      material_input = fis_gas_released_model
      output_initial_moles = initial_moles
      R = 8.3143
      temperature = interior_temp
      volume = gas_volume
      output = plenum_pressure
    []
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temp
    inlet_temperature = 580 # K
    inlet_pressure = 15.5e6 # PA
    inlet_massflux = 3880 # kg/m^2-sec
    rod_diameter = 0.95e-2 # m
    rod_pitch = 1.26e-2 # m
    linear_heat_rate = power_history
    axial_power_profile = axial_peaking_factors
  []
[]

[Materials]
  [fuel_thermal]
    type = UO2Thermal
    block = 3
    thermal_conductivity_model = NFIR
    temperature = temp
    burnup_function = burnup
    initial_porosity = 0.05
  []
  [fuel_elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 3
  []
  [fuel_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = 3
    youngs_modulus = 2.0e11
    poissons_ratio = 0.345
  []
  [fuel_relocation]
    type = UO2RelocationEigenstrain
    block = 3
    burnup_function = burnup
    diameter = 0.00836
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    diametral_gap =50.0e-6
    relocation_activation1 = 5000
    burnup_relocation_stop = 0.02
    eigenstrain_name = 'fuel_relocation_eigenstrain'
  []
  [fuel_swelling]
    type = UO2VolumetricSwellingEigenstrain
    block = 3
    burnup_function = burnup
    temperature = temp
    initial_fuel_density = 10431.0
    eigenstrain_name = 'fuel_volumetric_swelling_eigenstrain'
  []
  [fuel_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = 3
    # thermal_expansion_coeff = 10.0e-6 (reference)
    # We are artificially increasing the fuel expansion to simulate mechanical contact within reasonable 'example' time
    thermal_expansion_coeff = 50.0e-6
    temperature = temp
    stress_free_temperature = 580.0
    eigenstrain_name = 'fuel_thermal_strain'
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = 3
    temperature = temp
    fission_rate = fission_rate
    grain_radius = grain_radius
    gbs_model = true
    burnup = burnup
    diff_coeff_option = TURNBULL_D1_D2
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = 3
    strain_free_density = ${initial_fuel_density}
  []

  [clad_thermal]
    type = HeatConductionMaterial
    block = 1
    thermal_conductivity = 16.0
    specific_heat = 330.0
  []
  [clad_elasticity_tensor]
    type = ZryElasticityTensor
    block = 1
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep'
    block = 1
  []
  [clad_zrycreep]
    type = ZryCreepLimbackHoppeUpdate
    block = 1
    temperature = temp
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    model_irradiation_creep = true
    model_primary_creep = true
    model_thermal_creep = true
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = 1
    temperature = temp
    stress_free_temperature = 580.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [irradiation_swelling]
    type = ZryIrradiationGrowthEigenstrain
    block = 1
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = stress_relief_annealed
    eigenstrain_name = clad_irradiation_strain
  []
  [clad_density]
    type = StrainAdjustedDensity
    block = 1
    strain_free_density = 6551.0
  []
[]

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

[Executioner]
  type = Transient
  solve_type = 'NEWTON'

  petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -mat_mffd_err -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu    superlu_dist   1e-6          NONZERO               1e-13'
  line_search = 'basic'

  l_max_its = 10
  nl_max_its = 40
  nl_rel_tol = 1e-7
  nl_abs_tol = 1e-8

  start_time = -200
  dtmin = 1.0
  end_time = 25200

  nl_div_tol = 1e+40
  # For a regular thermal expansion value, use following end_time
  # end_time = 4.0e7

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 200
    optimal_iterations = 15
    iteration_window = 3
    growth_factor = 2.0
    cutback_factor = 0.5
  []
[]

[Postprocessors]
  [clad_inner_vol]
    type = InternalVolume
    boundary = 7
    scale_factor = 1.0
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [pellet_volume]
    type = InternalVolume
    boundary = 8
    scale_factor = 1.0
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [gas_volume]
    type = InternalVolume
    boundary = 9
    scale_factor = 1.0
    execute_on = 'initial linear'
  []
  [interior_temp]
    type = SideAverageValue
    boundary = 7
    variable = temp
    execute_on = 'initial linear'
  []
  [fis_gas_produced_model]
    type = ElementIntegralFisGasGeneratedSifgrs
    block = 3
    execute_on = 'initial timestep_end'
  []
  [fis_gas_produced]
    type = ScalePostprocessor
    value = fis_gas_produced_model
    scaling_factor = 1.0
    execute_on = 'initial timestep_end'
  []
  [fis_gas_released_model]
    type = ElementIntegralFisGasReleasedSifgrs
    block = 3
    execute_on = 'initial timestep_end'
  []
  [fission_gas_released]
    type = ScalePostprocessor
    value = fis_gas_released_model
    scaling_factor = 1.0
    execute_on = 'initial timestep_end'
  []
  [fis_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    block = 3
    outputs = exodus
  []
  [fis_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    block = 3
    outputs = exodus
  []
  [power_history]
    type = FunctionValuePostprocessor
    function = power_history
    execute_on = 'initial timestep_end'
  []
  [flux_from_clad_model]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 5
    diffusivity = thermal_conductivity
  []
  [flux_from_clad]
    type = ScalePostprocessor
    value = flux_from_clad_model
    scaling_factor = 1.0
    execute_on = 'initial timestep_end'
  []
  [flux_from_fuel_model]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 10
    diffusivity = thermal_conductivity
  []
  [flux_from_fuel]
    type = ScalePostprocessor
    value = flux_from_fuel_model
    scaling_factor = 1.0
    execute_on = 'initial timestep_end'
  []
  [average_burnup]
    type = RodAverageBurnup
    burnup_function = burnup
  []
  [average_fissionrate]
    type = ElementAverageValue
    block = 3
    variable = fission_rate
    execute_on = 'initial timestep_end'
  []
  [rod_total_power_model]
    type = ElementIntegralPower
    variable = temp
    fission_rate = fission_rate
    block = 3
    execute_on = 'initial timestep_end'
  []
  [rod_total_power]
    type = ScalePostprocessor
    value = rod_total_power_model
    scaling_factor = 1.0
    execute_on = 'initial timestep_end'
  []
  [rod_input_power]
    type = FunctionValuePostprocessor
    function = power_history
    scale_factor = 0.03
    execute_on = 'initial timestep_end'
  []
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  time_step_interval =  1
  exodus = true
  [checkpoint]
    type = Checkpoint
    time_step_interval =  1
    file_base = ckpoint_mortar_dup
    num_files = 2
  []
  [console]
    type = Console
    solve_log = true
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    execute_on = FINAL
    show = 'fission_gas_released plenum_pressure interior_temp gas_volume'
  []
[]

(examples/3D_rodlet_3pellets/discrete_full/3d_3pellets_mortar.i)


initial_fuel_density = 10431.0

[GlobalParams]
  density = ${initial_fuel_density} #95% TD (TD = 10980)
  displacements = 'disp_x disp_y disp_z'
  order = FIRST
  family = LAGRANGE
  energy_per_fission = 3.2e-11 # J/fission (205 Mev)
  volumetric_locking_correction = true
[]

[Problem]
  type = ReferenceResidualProblem
  reference_vector = 'ref'
  extra_tag_vectors = 'ref'
  group_variables = 'disp_x disp_y disp_z'
  converge_on = 'disp_x disp_y disp_z temp'
[]

[Mesh]
  [file]
    type = FileMeshGenerator
    file = DiscreteThreePellets3D_full_HEX8.e
  []
  partitioner = centroid
  centroid_partitioner_direction = y
  patch_size = 5
  patch_update_strategy = auto
[]

[Variables]
  [temp]
    initial_condition = 580.0
  []
[]

[AuxVariables]
  [fission_rate]
    block = 3
  []
  [fast_neutron_flux]
  []
  [fast_neutron_fluence]
  []
  [grain_radius]
    block = 3
    initial_condition = 5e-6
  []
  [frictional_status]
    family = LAGRANGE
    order = FIRST
  []
  [gas_swell]
    order = CONSTANT
    family = MONOMIAL
    block = 3
  []
  [hoop_inelastic_strain]
    order = CONSTANT
    family = MONOMIAL
    block = 1
  []
  [gap_conductance]
    order = CONSTANT
    family = MONOMIAL
  []
  [pid]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [power_history]
    type = PiecewiseLinear
    x = '0 1.0e4 1.0e8'
    y = '0  1.0  1.0'
    scale_factor = 25e3 # 25 kW/m peak power.
  []
  [axial_peaking_factors]
    type = ParsedFunction
    expression = 1
  []
  [coolant_pressure_ramp]
    type = PiecewiseLinear
    x = '-200 0'
    y = '0 1'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  [pellets]
    block = 3
    add_variables = true
    strain = FINITE
    eigenstrain_names = 'fuel_relocation_eigenstrain fuel_thermal_strain fuel_volumetric_swelling_eigenstrain'
    generate_output = 'vonmises_stress hydrostatic_stress stress_xx stress_yy stress_zz strain_xx strain_yy strain_zz'
    extra_vector_tags = 'ref'
    temperature = temp
  []
  [clad]
    block = 1
    add_variables = true
    strain = FINITE
    eigenstrain_names = 'clad_thermal_eigenstrain clad_irradiation_strain'
    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'
    extra_vector_tags = 'ref'
    temperature = temp
  []
[]

[Burnup]
  [burnup]
    block = 3
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    num_radial = 80
    num_axial = 11
    a_lower = 2.49e-3
    a_upper = 2.621e-2
    fuel_inner_radius = 0
    fuel_outer_radius = 0.0041
    fuel_volume_ratio = 1.0
    RPF = RPF
  []
[]

[Kernels]
  [heat]
    type = HeatConduction
    variable = temp
    extra_vector_tags = 'ref'
    block = '1 3'
  []
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = temp
    extra_vector_tags = 'ref'
    block = '1 3'
  []
  [heat_source_fuel]
    type = NeutronHeatSource
    variable = temp
    block = 3
    fission_rate = fission_rate
    extra_vector_tags = 'ref'
  []
[]

[AuxKernels]
  [fissionrate]
    type = FissionRateGeneral
    fission_rate_formulation = GENERIC
    variable = fission_rate
    block = 3
    value = 5.3548e+14
    fission_rate_function = power_history
  []
  [frictional_state]
    type = MortarFrictionalStateAux
    tangent_one = pellet_clad_mechanical_tangential_lm
    tangent_two = pellet_clad_mechanical_tangential_3d_lm
    boundary = 10
    contact_pressure = pellet_clad_mechanical_normal_lm
    variable = frictional_status
    mu = 0.5
  []
  [fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    block = 1
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    factor = 2.34e+13
    execute_on = timestep_begin
  []
  [fast_neutron_fluence]
    type = FastNeutronFluenceAux
    variable = fast_neutron_fluence
    block = 1
    fast_neutron_flux = fast_neutron_flux
    execute_on = timestep_begin
  []
  [grain_radius]
    type = GrainRadiusAux
    block = 3
    variable = grain_radius
    temperature = temp
    execute_on = linear
  []
  [gas_swell]
    type = MaterialRealAux
    block = 3
    variable = gas_swell
    property = deltav_v0_bubble_GB
    execute_on = timestep_end
  []
  [hoop_inelastic_strain]
    type = RankTwoScalarAux
    rank_two_tensor = creep_strain
    variable = hoop_inelastic_strain
    scalar_type = HoopStress
    execute_on = timestep_end
  []
  [pid]
    type = ProcessorIDAux
    variable = pid
  []
[]

[Contact]
  [pellet_clad_mechanical]
    primary = 5
    secondary = 10
    formulation = mortar
    model = coulomb
    c_normal = 1e+18
    c_tangential = 1e+18
    friction_coefficient = 0.5
  []
[]

[ThermalContactMortar]
  [thermal_contact]
    secondary_variable = temp
    primary_boundary = '5'
    secondary_boundary = '10'
    gas_released = fis_gas_released_model
    initial_moles = initial_moles
    jump_distance_model = LANNING
    plenum_pressure = plenum_pressure
    roughness_coef = 3.2
    roughness_secondary = 1e-6
    roughness_primary = 2e-6
    emissivity_primary = 0.8
    emissivity_secondary = 0.8
  []
[]

[BCs]
  [no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 1004
    value = 0.0
  []
  [no_z_all]
    type = DirichletBC
    variable = disp_z
    boundary = 1004
    value = 0.0
  []
  [no_y_clad_bottom]
    type = DirichletBC
    variable = disp_y
    boundary = 1
    value = 0.0
  []
  [no_x_clad_bottom]
    type = DirichletBC
    variable = disp_x
    boundary = 1
    value = 0.0
  []
  [no_z_clad_bottom]
    type = DirichletBC
    variable = disp_z
    boundary = 1
    value = 0.0
  []
  [no_z_fuel_bottom_point]
    type = DirichletBC
    variable = disp_z
    boundary = 1110
    value = 0.0
  []
  [no_z_fuel_bottom_point_y]
    type = DirichletBC
    variable = disp_y
    boundary = 1110
    value = 0.0
  []
  [no_x_fuel_bottom_point]
    type = DirichletBC
    variable = disp_x
    boundary = 1120
    value = 0.0
  []
  [no_x_fuel_bottom_point_y]
    type = DirichletBC
    variable = disp_y
    boundary = 1120
    value = 0.0
  []
  [Pressure]
    [coolantPressure]
      boundary = 2
      factor = 15.5e6
      function = coolant_pressure_ramp
    []
  []

  [PlenumPressure]
    [plenumPressure]
      boundary = 9
      initial_pressure = 2.0e6
      startup_time = 0.0
      material_input = fis_gas_released_model
      output_initial_moles = initial_moles
      R = 8.3143
      temperature = interior_temp
      volume = gas_volume
      output = plenum_pressure
    []
  []
[]

[CoolantChannel]
  [convective_clad_surface]
    boundary = 2
    variable = temp
    inlet_temperature = 580 # K
    inlet_pressure = 15.5e6 # PA
    inlet_massflux = 3880 # kg/m^2-sec
    rod_diameter = 0.95e-2 # m
    rod_pitch = 1.26e-2 # m
    linear_heat_rate = power_history
    axial_power_profile = axial_peaking_factors
  []
[]

[Materials]
  [fuel_thermal]
    type = UO2Thermal
    block = 3
    thermal_conductivity_model = NFIR
    temperature = temp
    burnup_function = burnup
    initial_porosity = 0.05
  []
  [fuel_elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = 3
  []
  [fuel_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    block = 3
    youngs_modulus = 2.0e11
    poissons_ratio = 0.345
  []
  [fuel_relocation]
    type = UO2RelocationEigenstrain
    block = 3
    burnup_function = burnup
    diameter = 0.00836
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    diametral_gap =50.0e-6
    relocation_activation1 = 5000
    burnup_relocation_stop = 0.02
    eigenstrain_name = 'fuel_relocation_eigenstrain'
  []
  [fuel_swelling]
    type = UO2VolumetricSwellingEigenstrain
    block = 3
    burnup_function = burnup
    temperature = temp
    initial_fuel_density = 10431.0
    eigenstrain_name = 'fuel_volumetric_swelling_eigenstrain'
  []
  [fuel_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = 3
    # thermal_expansion_coeff = 10.0e-6 (reference)
    # We are artificially increasing the fuel expansion to simulate mechanical contact within reasonable 'example' time
    thermal_expansion_coeff = 50.0e-6
    temperature = temp
    stress_free_temperature = 580.0
    eigenstrain_name = 'fuel_thermal_strain'
  []
  [fission_gas_release]
    type = UO2Sifgrs
    block = 3
    temperature = temp
    fission_rate = fission_rate
    grain_radius = grain_radius
    gbs_model = true
    burnup = burnup
    diff_coeff_option = TURNBULL_D1_D2
  []
  [fuel_density]
    type = StrainAdjustedDensity
    block = 3
    strain_free_density = ${initial_fuel_density}
  []

  [clad_thermal]
    type = HeatConductionMaterial
    block = 1
    thermal_conductivity = 16.0
    specific_heat = 330.0
  []
  [clad_elasticity_tensor]
    type = ZryElasticityTensor
    block = 1
  []
  [clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep'
    block = 1
  []
  [clad_zrycreep]
    type = ZryCreepLimbackHoppeUpdate
    block = 1
    temperature = temp
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    model_irradiation_creep = true
    model_primary_creep = true
    model_thermal_creep = true
  []
  [thermal_expansion]
    type = ZryThermalExpansionMATPROEigenstrain
    block = 1
    temperature = temp
    stress_free_temperature = 580.0
    eigenstrain_name = clad_thermal_eigenstrain
  []
  [irradiation_swelling]
    type = ZryIrradiationGrowthEigenstrain
    block = 1
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_type = stress_relief_annealed
    eigenstrain_name = clad_irradiation_strain
  []
  [clad_density]
    type = StrainAdjustedDensity
    block = 1
    strain_free_density = 6551.0
  []
[]

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

[Executioner]
  type = Transient
  solve_type = 'NEWTON'

  petsc_options = '-snes_converged_reason -ksp_converged_reason -snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -mat_mffd_err -pc_factor_shift_type -pc_factor_shift_amount'
  petsc_options_value = 'lu    superlu_dist   1e-6          NONZERO               1e-13'
  line_search = 'basic'

  l_max_its = 10
  nl_max_its = 40
  nl_rel_tol = 1e-7
  nl_abs_tol = 1e-8

  start_time = -200
  dtmin = 1.0
  end_time = 25200

  nl_div_tol = 1e+40
  # For a regular thermal expansion value, use following end_time
  # end_time = 4.0e7

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 200
    optimal_iterations = 15
    iteration_window = 3
    growth_factor = 2.0
    cutback_factor = 0.5
  []
[]

[Postprocessors]
  [clad_inner_vol]
    type = InternalVolume
    boundary = 7
    scale_factor = 1.0
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [pellet_volume]
    type = InternalVolume
    boundary = 8
    scale_factor = 1.0
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [gas_volume]
    type = InternalVolume
    boundary = 9
    scale_factor = 1.0
    execute_on = 'initial linear'
  []
  [interior_temp]
    type = SideAverageValue
    boundary = 7
    variable = temp
    execute_on = 'initial linear'
  []
  [fis_gas_produced_model]
    type = ElementIntegralFisGasGeneratedSifgrs
    block = 3
    execute_on = 'initial timestep_end'
  []
  [fis_gas_produced]
    type = ScalePostprocessor
    value = fis_gas_produced_model
    scaling_factor = 1.0
    execute_on = 'initial timestep_end'
  []
  [fis_gas_released_model]
    type = ElementIntegralFisGasReleasedSifgrs
    block = 3
    execute_on = 'initial timestep_end'
  []
  [fission_gas_released]
    type = ScalePostprocessor
    value = fis_gas_released_model
    scaling_factor = 1.0
    execute_on = 'initial timestep_end'
  []
  [fis_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    block = 3
    outputs = exodus
  []
  [fis_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    block = 3
    outputs = exodus
  []
  [power_history]
    type = FunctionValuePostprocessor
    function = power_history
    execute_on = 'initial timestep_end'
  []
  [flux_from_clad_model]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 5
    diffusivity = thermal_conductivity
  []
  [flux_from_clad]
    type = ScalePostprocessor
    value = flux_from_clad_model
    scaling_factor = 1.0
    execute_on = 'initial timestep_end'
  []
  [flux_from_fuel_model]
    type = SideDiffusiveFluxIntegral
    variable = temp
    boundary = 10
    diffusivity = thermal_conductivity
  []
  [flux_from_fuel]
    type = ScalePostprocessor
    value = flux_from_fuel_model
    scaling_factor = 1.0
    execute_on = 'initial timestep_end'
  []
  [average_burnup]
    type = RodAverageBurnup
    burnup_function = burnup
  []
  [average_fissionrate]
    type = ElementAverageValue
    block = 3
    variable = fission_rate
    execute_on = 'initial timestep_end'
  []
  [rod_total_power_model]
    type = ElementIntegralPower
    variable = temp
    fission_rate = fission_rate
    block = 3
    execute_on = 'initial timestep_end'
  []
  [rod_total_power]
    type = ScalePostprocessor
    value = rod_total_power_model
    scaling_factor = 1.0
    execute_on = 'initial timestep_end'
  []
  [rod_input_power]
    type = FunctionValuePostprocessor
    function = power_history
    scale_factor = 0.03
    execute_on = 'initial timestep_end'
  []
[]

[PerformanceMetricOutputs]
[]

[Outputs]
  perf_graph = true
  time_step_interval =  1
  exodus = true
  [checkpoint]
    type = Checkpoint
    time_step_interval =  1
    file_base = ckpoint_mortar_dup
    num_files = 2
  []
  [console]
    type = Console
    solve_log = true
    output_linear = true
    max_rows = 25
  []
  [chkfile]
    type = CSV
    execute_on = FINAL
    show = 'fission_gas_released plenum_pressure interior_temp gas_volume'
  []
[]

Contact mechanics action
LWR gap heat transfer action

3D Light Water Reactor Fuel Rod Tutorial

Node-to-surface input files

Mortar input files

Contact mechanics action

LWR gap heat transfer action

(examples/3D_rodlet_3pellets/smeared/smearedTest3D.i)

(examples/3D_rodlet_3pellets/smeared/powerhistory.csv)

(examples/3D_rodlet_3pellets/discrete_full/3d_3pellets_mortar.i)

(examples/3D_rodlet_3pellets/discrete_full/3d_3pellets_mortar.i)