Sodium Cooled Metal Fuel Reactor Tutorial

[GlobalParams]
  density = 15800.0
  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_x saved_t'
[]

[Mesh]
  # file = x447.e     # Uncomment if using Exodus mesh file and comment lines from type to elem_type
  # rod specific parameters
  type = SmearedPelletMesh
  clad_thickness = 0.381e-03
  pellet_outer_radius = 2.192e-03
  pellet_height = 343.0e-3
  clad_top_gap_height = 271.3e-3
  clad_gap_width = 0.348e-03
  top_bot_clad_height = 2.24e-3 # arbitrary
  clad_bot_gap_height = 0.31e-3 # arbitrary
  # meshing parameters
  clad_mesh_density = customize
  pellet_mesh_density = customize
  nx_p = 6
  ny_p = 260
  nx_c = 4
  ny_c = 260
  ny_cu = 3
  ny_cl = 3
  pellet_quantity = 1
  elem_type = QUAD8
  # mesh options
  patch_size = 50
  patch_update_strategy = auto
  partitioner = centroid
  centroid_partitioner_direction = y
[]

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

[AuxVariables]
  [./saved_x]
  [../]
  [./saved_y]
  [../]
  [./saved_t]
  [../]
  [./fission_rate]
    block = pellet
  [../]
  [./burnup]
    block = pellet
  [../]
  [./fast_neutron_flux]
    block = clad
  [../]
  [./fast_neutron_fluence]
    block = clad
  [../]

  # Aux variables for output
  [./stress_xx]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_yy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_zz]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./vonmises]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./hydrostatic_stress]
    block = pellet
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./creep_strain_mag]
    order = CONSTANT
    family = MONOMIAL
  [../]

  [./gap_cond]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./cumulative_damage_index]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./element_failed]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./volumetric_strain]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[Functions]
  [./power_history]
    type = PiecewiseLinear
    x = '0 1e5'
    y = '0 29274'
  [../]
  [./axial_peaking_factors]
    type = PowerPeakingFunction
    fit = EBR-II_low
    pellet_length = 343.0e-3
    pellet_y_start = 2.55e-3
  [../]
[]

[Kernels]
  # Define kernels for the various terms in the PDE system
  [./TensorMechanics] #continuum mechanics stress divergence
    use_displaced_mesh = true #Incremental formulation
    save_in = 'saved_x saved_y'
  [../]
  [./gravity]
    type = Gravity
    variable = disp_y
    value = -9.81
    save_in = 'saved_x'
  [../]
  [./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
    fission_rate = fission_rate
    save_in = 'saved_t'
  [../]
[]

[AuxKernels]
  [./fission_rate]
    type = FissionRateAuxLWR
    variable = fission_rate
    block = pellet
    axial_power_profile = axial_peaking_factors
    rod_ave_lin_pow = power_history
    pellet_diameter = 4.384e-3
    execute_on = timestep_begin
  [../]
  [./burnup]
    type = BurnupMetalAux
    block = pellet
    fission_rate = fission_rate
    variable = burnup
    execute_on = timestep_begin
  [../]
  [./fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    block = clad
    factor = 0.6e19
    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
  [../]

  [./stress_xx]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_xx
    index_j = 0
    index_i = 0
    execute_on = timestep_end
  [../]
  [./stress_yy]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_yy
    index_j = 1
    index_i = 1
    execute_on = timestep_end
  [../]
  [./stress_zz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zz
    index_j = 2
    index_i = 2
    execute_on = timestep_end
  [../]
  [./vonmises]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = vonmises
    scalar_type = VonMisesStress
    execute_on = timestep_end
  [../]
  [./hydrostatic_stress]
    type = RankTwoScalarAux
    rank_two_tensor = stress
    variable = hydrostatic_stress
    scalar_type = Hydrostatic
    execute_on = timestep_end
  [../]
  [./creep_strain_mag]
    type = RankTwoScalarAux
    block = clad
    rank_two_tensor = creep_strain
    variable = creep_strain_mag
    scalar_type = EffectiveStrain
    execute_on = timestep_end
  [../]

  [./conductance]
    type = MaterialRealAux
    property = gap_conductance
    variable = gap_cond
    boundary = 10
  [../]

  [./cdf_amount]
    boundary = '1 2 3'
    type = MaterialRealAux
    property = cdf_failure
    variable = cumulative_damage_index
  [../]
  [./failed_element]
    boundary = '1 2 3'
    type = MaterialRealAux
    property = failed
    variable = element_failed
  [../]
  [./volumetric_strain]
    type = RankTwoScalarAux
    rank_two_tensor = total_strain
    variable = volumetric_strain
    scalar_type = VolumetricStrain
    execute_on = timestep_end
  [../]
[]

[Contact]
  [./pellet_clad_mechanical]
    master = 5
    slave = 10
    penalty = 1e14
    model = coulomb
    formulation = tangential_penalty
    friction_coefficient = 0.4
    system = constraint
    normalize_penalty = true
    tangential_tolerance = 1e-3
    normal_smoothing_distance = 0.1
  [../]
[]

[ThermalContact]
  [./thermal_contact]
    type = GapHeatTransfer
    variable = temp
    master = 5
    slave = 10
    quadrature = true
    gap_conductivity = 61.0
    min_gap = 1e-6
  [../]
[]

[BCs]
  [./no_x_all]
    type = DirichletBC
    variable = disp_x
    boundary = 12
    value = 0.0
  [../]
  [./no_y_fuel]
    type = PenaltyDirichletBC
    penalty = 1e10
    variable = disp_y
    boundary = 20
    value = 0.0
  [../]
  [./no_y_clad]
    type = DirichletBC
    variable = disp_y
    boundary = 1
    value = 0.0
  [../]

  [./Pressure]
    [./coolantPressure]
      boundary = '1 2 3'
      factor = 0.151e6 # Pa
    [../]
  [../]
  [./PlenumPressure]
    [./plenumPressure]
      boundary = 9
      initial_pressure = 0.084e6 # Pa
      startup_time = 0
      R = 8.3143
      temperature = ave_temp_interior
      volume = gas_volume
      output = plenum_pressure
      material_input = fis_gas_released
    [../]
  [../]
[]

[CoolantChannel]
  [./convective_clad_surface]
    boundary = '1 2 3'
    variable = temp
    inlet_temperature   = 648.0    # K
    inlet_pressure      = 0.151e6  # Pa
    inlet_massflux      = 2437.8   # kg/m^2-sec
    coolant_material    = sodium
    rod_diameter        = 5.84e-3  # m
    rod_pitch           = 6.912e-3 # m
    linear_heat_rate    = power_history
    axial_power_profile = axial_peaking_factors
    subchannel_geometry = triangular
  [../]
[]

[Materials]
  [./fuel_elasticity_tensor]
    type = UPuZrElasticityTensor
    block = pellet
    temperature = temp
  [../]
  [./fuel_elastic_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = nonlinear
    inelastic_models = 'fuel_upuzrcreep'
    block = pellet
  [../]
  [./fuel_upuzrcreep]
    type = UPuZrCreepUpdate
    block = pellet
    fission_rate = fission_rate
    temperature = temp
    porosity = porosity
  [../]
  [./fuel_strain]
    type = ComputeAxisymmetricRZFiniteStrain
    block = pellet
    eigenstrain_names = 'fuel_thermal_strain fuel_volumetric_strain'
  [../]
  [./fuel_thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = pellet
    thermal_expansion_coeff = 1.18e-5
    temperature = temp
    stress_free_temperature = 295.0
    eigenstrain_name = fuel_thermal_strain
  [../]
  [./fuel_volumetric_swelling]
    type = UPuZrVolumetricSwellingEigenstrain
    block = pellet
    temp = temp
    fission_rate = fission_rate
    burnup = burnup
    hydrostatic_stress = hydrostatic_stress
    eigenstrain_name = fuel_volumetric_strain
    outputs = all
    output_properties = 'porosity gaseous_porosity gas_swelling solid_swelling'
  [../]
  [./metal_fuel_thermal]
    type = ThermalUPuZr
    block = pellet
    X_Zr = 0.225
    X_Pu = 0
    spheat_model = savage
    thcond_model = lanl
    porosity = porosity
    temp = temp
  [../]
  [./fuel_density]
    type = Density
    block = pellet
  [../]
  [./Fission_Gas_Release]
    type = FgrUPuZr
    block = pellet
    fission_rate = fission_rate
  [../]

  [./clad_elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 1.88e11
    poissons_ratio = 0.236
    block = clad
  [../]
  [./clad_stress]
    type = ComputeMultipleInelasticStress
    tangent_operator = nonlinear
    inelastic_models = 'clad_ht9creep'
    block = clad
  [../]
  [./clad_strain]
    type = ComputeAxisymmetricRZFiniteStrain
    block = clad
    eigenstrain_names = 'clad_thermal_eigenstrain'
  [../]
  [./clad_ht9creep]
    type = HT9CreepUpdate
    block = clad
    temperature = temp
    fast_neutron_flux = fast_neutron_flux
  [../]
  [./thermal_expansion]
    type = ComputeThermalExpansionEigenstrain
    block = clad
    thermal_expansion_coeff = 1.2e-5
    temperature = temp
    stress_free_temperature = 293.0
    eigenstrain_name = clad_thermal_eigenstrain
  [../]
  [./clad_thermal]
    type = ThermalHT9
    block = clad
    temp = temp
  [../]
  [./clad_density]
    type = Density
    block = clad
    density = 7874.0
  [../]
  [./longHT9_failure]
    type = FailureCladHT9
    boundary = '1 2 3'
    method = cdf_long
    temperature = temp
    hoop_stress = stress_zz # Since 2D-RZ
  [../]
[]

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

[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 = 100
  l_tol = 1e-3

  nl_max_its = 50
  nl_rel_tol = 1e-4
  nl_abs_tol = 1e-8

  end_time = 5.685e+07
  dtmin = 1
  dtmax = 1e5

  [./Quadrature]
    order = fifth
    side_order = seventh
  [../]
  [./TimeStepper]
    type = IterationAdaptiveDT
    postprocessor_dtlim = creep_timestep
    dt = 1e2
    iteration_window = 4
    optimal_iterations = 10
  [../]
[]

[Postprocessors]
  [./_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
  [../]

  [./ave_temp_interior]
    type = SideAverageValue
    boundary = 9
    variable = temp
    execute_on = 'initial linear'
    outputs = exodus
  [../]
  [./avg_clad_temp]
    type = ElementAverageValue
    variable = temp
    outputs = exodus
    block = clad
  [../]
  [./peak_clad_temp]
    type = ElementExtremeValue
    variable = temp
    value_type = max
    block = clad
  [../]
  [./peak_fuel_temp]
    type = ElementExtremeValue
    variable = temp
    value_type = max
    block = pellet
  [../]

  [./max_hydro]
    type = ElementExtremeValue
    variable = hydrostatic_stress
    value_type = max
    block = pellet
  [../]
  [./min_hydro]
    type = ElementExtremeValue
    variable = hydrostatic_stress
    value_type = min
    block = pellet
  [../]

  [./peak_porosity]
    type = ElementExtremeValue
    variable = porosity
    value_type = max
    block = pellet
  [../]
  [./clad_inner_vol]
    type = InternalVolume
    boundary = 7
    outputs = exodus
  [../]
  [./pellet_volume]
    type = InternalVolume
    boundary = 8
    outputs = exodus
  [../]
  [./gas_volume]
    type = InternalVolume
    boundary = 9
    execute_on = 'initial linear'
    outputs = exodus
  [../]
  [./clad_fuel_gap]
    type = NodalMaxValue
    variable = penetration
    boundary = 10
  [../]

  [./flux_from_clad]
    type = SideFluxIntegral
    variable = temp
    boundary = 5
    diffusivity = thermal_conductivity
    outputs = exodus
  [../]
  [./flux_from_fuel]
    type = SideFluxIntegral
    variable = temp
    boundary = 10
    diffusivity = thermal_conductivity
    outputs = exodus
  [../]

  [./rod_total_power]
    type = ElementIntegralPower
    variable = temp
    fission_rate = fission_rate
    block = pellet
    outputs = exodus
  [../]
  [./LHGR_W_per_cm]
    type = FunctionValuePostprocessor
    function = power_history
    scale_factor = 0.01
  [../]
  [./average_burnup]
    type = ElementAverageValue
    block = pellet
    variable = burnup
  [../]
  [./max_cdf]
    type = ElementExtremeValue
    value_type = max
    variable = cumulative_damage_index
  [../]
  [./fis_gas_produced]
    type = ElementIntegralFisGasProduce
    variable = temp # dummy variable
    block = pellet
  [../]
  [./fis_gas_released]
    type = ElementIntegralFisGasRelease
    variable = temp # dummy variable
    block = pellet
    execute_on = 'initial timestep_end'
  [../]

  [./creep_timestep]
    type = MaterialTimeStepPostprocessor
    block = pellet
  [../]
  [./disp_x_9076]
    type = NodalVariableValue
    nodeid = 9075
    variable = disp_x
  [../]
  [./disp_y_9076]
    type = NodalVariableValue
    nodeid = 9075
    variable = disp_y
  [../]
  [./hydrostatic_stress]
    type = ElementAverageValue
    variable = hydrostatic_stress
    execute_on = 'initial timestep_end'
    block = pellet
  [../]

  [./solid_swelling]
    type = ElementAverageValue
    variable = solid_swelling
    block = pellet
  [../]
  [./gas_swelling]
    type = ElementAverageValue
    variable = gas_swelling
    block = pellet
  [../]
  [./volumetric_strain]
    type = ElementAverageValue
    variable = volumetric_strain
    block = pellet
  [../]
  [./fission_rate]
    type = ElementAverageValue
    variable = fission_rate
    block = pellet
  [../]
  [./porosity]
    type = ElementAverageValue
    variable = porosity
    block = pellet
  [../]
[]

[VectorPostprocessors]
  [./clad_x_disp]
    type = NodalValueSampler
    variable = disp_x
    boundary = 2
    sort_by = y
  [../]
  [./fuel_cl_temp]
    type = NodalValueSampler
    variable = temp
    boundary = 12
    sort_by = y
  [../]
  [./fuel_surf_temp]
    type = NodalValueSampler
    variable = temp
    boundary = 10
    sort_by = y
  [../]
  [./clad_inn_temp]
    type = NodalValueSampler
    variable = temp
    boundary = 5
    sort_by = y
  [../]
  [./clad_out_temp]
    type = NodalValueSampler
    variable = temp
    boundary = 2
    sort_by = y
  [../]
[]

[Outputs]
  interval = 10
  color = false
  exodus = true
  print_perf_log = true
  csv = true
  sync_times = '1e3 5e3 1e4 5e4 1e5 5e6 1e6 5e6 10e6 15e6 20e6 25e6 30e6 35e6 40e6 45e6 50e6'
  [./console]
    type = Console
    max_rows = 25
    interval = 1
    output_linear = true
  [../]
[]
(examples/metal_fuel/uzr/x447_tm.i)