Accident Tolerant Fuel and Cladding Tutorial

[GlobalParams]
  # Set initial fuel density, other global parameters
  density = 11590.0
  disp_x = disp_x
  disp_y = disp_y
  order = SECOND
  family = LAGRANGE
  energy_per_fission = 3.2e-11  # J/fission
  volumetric_locking_correction = true
  displacements = 'disp_x disp_y'
[]

[Problem]
  coord_type = RZ
  type = ReferenceResidualProblem
  solution_variables = 'disp_x disp_y temp'
  reference_residual_variables = 'saved_x saved_y saved_t'
[]

[Mesh]
  # Import mesh file
  file = u3si2_zircaloy_smeared.e
  patch_size = 10 # For contact algorithm
  patch_update_strategy = auto
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[Variables]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./temp]
    initial_condition = 293.0
  [../]
[]

[UserObjects]
  [./pin_geometry]
    type = FuelPinGeometry
    clad_inner_wall = 5
    clad_outer_wall = 2
    clad_top = 3
    clad_bottom = 1
    pellet_exteriors = 8
  [../]
[]

[AuxVariables]
  # Define auxilary variables
  [./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
  [../]
  [./hoop_stress]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./fract_beta_phase]
    order  = CONSTANT
    family = MONOMIAL
  [../]
  [./coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./oxide_thickness]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./total_hoop_strain]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./saved_x]
  [../]
  [./saved_y]
  [../]
  [./saved_t]
  [../]
  [./creep_rate]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./densification]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./solid_swell]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./gaseous_swell]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[Functions]
  [./power_history]
    type = PiecewiseLinear
    data_file = powerhistory.csv
    scale_factor = 1
  [../]
  [./axial_peaking_factors]
    type = ParsedFunction
    value = 1
  [../]
  [./pressure_ramp]
    type = PiecewiseLinear
    x = '-200 0 1e8'
    y = '6.537e-3 1 1'
    scale_factor = 15.5e6
  [../]
  [./mass_flux_func]
    type = PiecewiseLinear
    x = '-200 0 1e8'
    y = '3800 3800 3800'
  [../]
  [./q]
    type = CompositeFunction
    functions = 'power_history axial_peaking_factors'
  [../]
[]

[Modules/TensorMechanics/Master]
  [./pellets]
    block = pellet_type_1
    add_variables = true
    strain = FINITE
    eigenstrain_names = 'fuel_thermal_strain fuel_volumetric_strain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
    save_in = 'saved_x saved_y'
  [../]
  [./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'
    save_in = 'saved_x saved_y'
  [../]
[]

[Kernels]
  [./gravity]
    type = Gravity
    variable = disp_y
    value = -9.81
    save_in = saved_y
  [../]
  [./heat]
    type = HeatConduction
    variable = temp
    save_in = saved_t
  [../]
  [./heat_ie]
    type = HeatConductionTimeDerivative
    variable = temp
    save_in = saved_t
  [../]
  [./heat_source]
     type = NeutronHeatSource
     variable = temp
     block = pellet_type_1
     burnup_function = burnup
     save_in = saved_t
  [../]
[]

[Burnup]
  [./burnup]
    block = pellet_type_1
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    num_radial = 81
    num_axial = 11
    fuel_pin_geometry = pin_geometry
    fuel_volume_ratio = 1.0
    RPF = RPF
  [../]
[]

[AuxKernels]
  [./fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    block = clad
    rod_ave_lin_pow = power_history
    axial_power_profile = axial_peaking_factors
    factor = 3e13
    execute_on = timestep_begin
  [../]
  [./fast_neutron_fluence]
    type = FastNeutronFluenceAux
    variable = fast_neutron_fluence
    block = clad
    fast_neutron_flux = fast_neutron_flux
    execute_on = timestep_begin
  [../]
  [./grain_radius]
    type = GrainRadiusAux
    block = pellet_type_1
    variable = grain_radius
    temp = temp
    execute_on = linear
  [../]
  [./hoop_stress]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = hoop_stress
    scalar_type = HoopStress
    execute_on = timestep_end
  [../]
  [./total_hoop_strain]
    type = RankTwoScalarAux
    rank_two_tensor = total_strain
    variable = total_hoop_strain
    scalar_type = HoopStress
    execute_on = timestep_end
  [../]
  [./conductance]
    type = MaterialRealAux
    property = gap_conductance
    variable = gap_cond
    boundary = 10
  [../]
  [./coolant_htc]
    type = MaterialRealAux
    property = coolant_channel_htc
    variable = coolant_htc
    boundary = '1 2 3'
  [../]
  [./oxide]
    type = MaterialRealAux
    variable = oxide_thickness
    property = scale_thickness
    boundary = 2
  [../]
  [./creep_rate]
    type = MaterialRealAux
    variable = creep_rate
    property = creep_rate
    execute_on = timestep_end
    block = clad
  [../]
  [./fract_bphase]
    type = MaterialRealAux
    variable = fract_beta_phase
    property = fract_beta_phase
    block = clad
  [../]
  [./densfication]
    type = MaterialRealAux
    property = densification
    variable = densification
    block = pellet_type_1
  [../]
  [./solid_swell]
    type = MaterialRealAux
    property = solid_swelling
    variable = solid_swell
    block = pellet_type_1
  [../]
  [./gaseous_swell]
    type = MaterialRealAux
    property = gaseous_swelling
    variable = gaseous_swell
    block = pellet_type_1
  [../]
[]

[Contact]
  [./pellet_clad_mechanical]
    master = 5
    slave = 10
    system = constraint
    normal_smoothing_distance = 0.1
    disp_x = disp_x
    disp_y = disp_y
    penalty = 1e7
  [../]
[]

[ThermalContact]
  [./thermal_contact]
    type = GapHeatTransferLWR
    variable = temp
    master = 5
    slave = 10
    initial_moles = initial_moles
    gas_released = fis_gas_released
    contact_pressure = contact_pressure
    quadrature = true
  [../]
[]

[BCs]
  [./no_x_all] # pin pellets and clad along axis of symmetry (y)
    type = DirichletBC
    variable = disp_x
    boundary = 12
    value = 0.0
  [../]
  [./no_y_clad_bottom] # pin clad bottom in the axial direction (y)
    type = DirichletBC
    variable = disp_y
    boundary = '1'
    value = 0.0
  [../]
  [./no_y_fuel_bottom] # pin fuel bottom in the axial direction (y)
    type = DirichletBC
    variable = disp_y
    boundary = '1020'
    value = 0.0
  [../]
  [./Pressure]
    [./coolantPressure]
      boundary = '1 2 3'
      function = pressure_ramp
    [../]
  [../]
  [./PlenumPressure]
    [./plenumPressure]
      boundary = 9
      initial_pressure = 2.0e6
      startup_time = 0
      R = 8.3143
      output_initial_moles = initial_moles
      temperature = ave_temp_interior
      volume = gas_volume
      material_input = fis_gas_released
      output = plenum_pressure
    [../]
  [../]
[]

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

[Materials]
  [./fuel_thermal]
    type = ThermalSilicideFuel
    block = pellet_type_1
    thermal_conductivity_model = WHITE
    silicon_mole_fraction = 0.4
    temp = temp
  [../]
  [./fuel_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1.4e11
    poissons_ratio = 0.17
    block = pellet_type_1
  [../]
  [./fuel_elastic_stress]
    type = ComputeFiniteStrainElasticStress
    block = pellet_type_1
  [../]
  [./fuel_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = pellet_type_1
    temperature = temp
    stress_free_temperature = 295.0
    thermal_expansion_coeff = 15.0e-6
    eigenstrain_name = fuel_thermal_strain
  [../]
  [./fuel_volumetric_swelling]
    type = U3Si2VolumetricSwellingEigenstrain
    block = pellet_type_1
    temperature = temp
    burnup_function = burnup
    save_densification = true
    save_solid_swelling = true
    eigenstrain_name = fuel_volumetric_strain
  [../]

  [./clad_thermal]
    type = ThermalZry
    temp = temp
    block = clad
  [../]
  [./clad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 7.5e10
    poissons_ratio = 0.3
    block = clad
  [../]
  [./clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = elastic
    inelastic_models = 'clad_zrycreep clad_plasticity'
    relative_tolerance = 1e-5
    block = clad
  [../]
  [./phase]
    type = ZrPhase
    block = clad
    temperature = temp
    numerical_method = 2
  [../]
  [./clad_zrycreep]
    type = ZryCreepLOCAErbacherLimbackHoppeUpdate
    block = clad
    temperature = temp
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    model_irradiation_creep = true
    model_primary_creep = true
    relative_tolerance = 1e-5
    max_inelastic_increment = 1e-4
    temperature_standard_thermal_creep_end = 700
    temperature_loca_creep_begin = 900
    zircaloy_material_type = stress_relief_annealed
  [../]
  [./thermal_expansion]
    type = ZryThermalExpansionMatproEigenstrain
    block = clad
    temperature = temp
    stress_free_temperature = 295.0
    eigenstrain_name = clad_thermal_eigenstrain
  [../]
  [./irradiation_swelling]
    type = ZryIrradiationGrowthEigenstrain
    block = clad
    fast_neutron_fluence = fast_neutron_fluence
    zircaloy_material_model_type = stress_relief_annealed
    eigenstrain_name = clad_irradiation_strain
  [../]
  [./clad_plasticity]
    type = ZryPlasticityUpdate
    block = clad
    temperature = temp
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    relative_tolerance = 1e-5
    cold_work_factor = 0.5
    plasticity_model_type = PNNL
    zircaloy_alloy_type = 4
  [../]
  [./fission_gas_release]
    type = Sifgrs
    block = pellet_type_1
    temp = temp
    burnup_function = burnup
    grain_radius = grain_radius
    gbs_model = true
  [../]
  [./clad_density]
    type = Density
    block = clad
    density = 6511.0
    disp_r = disp_x
    disp_z = disp_y
  [../]
  [./fuel_density]
    type = Density
    block = pellet_type_1
    disp_r = disp_x
    disp_z = disp_y
  [../]
  [./oxidationcladding]
    type = OxidationCladding
    boundary = 2
    clad_inner_radius = 4.1783e-3
    clad_outer_radius = 4.7498e-3
    normtemp_model = 0
    hightemp_model = 0
    temperature = temp
    fast_neutron_flux = fast_neutron_flux
    use_coolant_channel = true
    oxygen_wtfract_initial = 0.0012
    oxidation_scalef = 1.0
  [../]
[]

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

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

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

  line_search = 'none'

  l_max_its = 100
  l_tol = 8e-3

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

  start_time = -200
  n_startup_steps = 1
  end_time = 1e8

  dtmax = 1e6
  dtmin = 1e-3

  [./TimeStepper]
    type = IterationAdaptiveDT
    dt = 2.0e2
    force_step_every_function_point = true
    timestep_limiting_function = power_history
    max_function_change = 3e20
    optimal_iterations = 10
    iteration_window = 2
    linear_iteration_ratio = 100
    postprocessor_dtlim = material_timestep
  [../]
  [./Quadrature]
    order = FIFTH
    side_order = SEVENTH
  [../]
[]

[Postprocessors]
  [./avg_fuel_surface]
    type = SideAverageValue
    boundary = 10
    variable = temp
  [../]
  [./ave_temp_interior]
    type = SideAverageValue
    boundary = 9
    variable = temp
    execute_on = 'initial linear'
  [../]
  [./clad_inner_vol]
    type = InternalVolume
    boundary = 7
  [../]
  [./pellet_volume]
    type = InternalVolume
    boundary = 8
  [../]
  [./avg_clad_temp]
    type = SideAverageValue
    boundary = 7
    variable = temp
  [../]
  [./fis_gas_produced]
    type = ElementIntegralFisGasGeneratedSifgrs
    variable = temp
    block = pellet_type_1
  [../]
  [./fis_gas_released]
    type = ElementIntegralFisGasReleasedSifgrs
    variable = temp
    block = pellet_type_1
  [../]
  [./fis_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    variable = temp
    block = pellet_type_1
  [../]
  [./fis_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    variable = temp
    block = pellet_type_1
  [../]
  [./gas_volume]
    type = InternalVolume
    boundary = 9
    execute_on = 'initial linear'
  [../]
  [./flux_from_clad]
    type = SideFluxIntegral
    variable = temp
    boundary = 5
    diffusivity = thermal_conductivity
  [../]
  [./flux_from_fuel]
    type = SideFluxIntegral
    variable = temp
    boundary = 10
    diffusivity = thermal_conductivity
  [../]
  [./_dt]
    type = TimestepSize
  [../]
  [./num_lin_it]
    type = NumLinearIterations
  [../]
  [./num_nonlin_it]
    type = NumNonlinearIterations
  [../]
  [./tot_lin_it]
    type = CumulativeValuePostprocessor
    postprocessor = num_lin_it
  [../]
  [./tot_nonlin_it]
    type = CumulativeValuePostprocessor
    postprocessor = num_nonlin_it
  [../]
  [./alive_time]
    type = PerformanceData
    event = ALIVE
  [../]
  [./rod_total_power]
    type = ElementIntegralPower
    variable = temp
    burnup_function = burnup
    block = pellet_type_1
  [../]
  [./rod_input_power]
    type = FunctionValuePostprocessor
    function = power_history
    scale_factor = 0.1186
  [../]
  [./average_burnup]
    type = ElementAverageValue
    block = pellet_type_1
    variable = burnup
  [../]
  [./oxide_thickness]
    type = ElementExtremeValue
    block = clad
    variable = oxide_thickness
  [../]
  [./fis_gas_percent]
    type = FGRPercent
    value1 = fis_gas_released
    value2 = fis_gas_produced
  [../]
  [./material_timestep]
    type = MaterialTimeStepPostprocessor
    block = clad
  [../]
[]

[Outputs]
  # Define output file(s)
  interval = 1
  exodus = true
  color = false
  csv = true
  print_linear_residuals = true
  print_perf_log = true
  [./console]
    type = Console
    max_rows = 25
  [../]
[]
(examples/tensor_mechanics/accident_tolerant_fuel/u3si2_zircaloy.i)