Axisymmetric (2D-RZ) Light Water Reactor Fuel Rod Tutorial

# This model is a linear element, 10 discrete fuel pellet stack (pellet_type_1) with a fine mesh.

[GlobalParams]
  # Set initial fuel density, other global parameters
  density = 10431.0
  initial_porosity = 0.05
  energy_per_fission = 3.2e-11  # J/fission
  volumetric_locking_correction = true
  displacements = 'disp_x disp_y'
[]

[Problem]
  # Specify coordinate system type
  coord_type = RZ
  type = ReferenceResidualProblem
  solution_variables = 'disp_x disp_y temp'
  reference_residual_variables = 'saved_x saved_x saved_t'
[]

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

[Variables]
  # Define dependent variables and initial conditions
  [./temp]
    initial_condition = 580.0     # set initial temp to coolant inlet
  [../]
[]

[AuxVariables]
  # Define auxilary variables
  [./saved_x]
  [../]
  [./saved_y]
  [../]
  [./saved_t]
  [../]
  [./fast_neutron_flux]
    block = clad
  [../]
  [./fast_neutron_fluence]
    block = clad
  [../]
  [./grain_radius]
    block = pellet_type_1
    initial_condition = 10e-6
  [../]
  [./creep_strain_rate]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./gap_cond]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./coolant_htc]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[Functions]
  # Define functions to control power and boundary conditions
  [./power_history]
    type = PiecewiseLinear   # reads and interpolates an input file containing rod average linear power vs time
    data_file = powerhistory.csv
    scale_factor = 1
  [../]
  [./axial_peaking_factors]      # reads and interpolates an input file containing the axial power profile vs time
    type = PiecewiseBilinear
    data_file = peakingfactors.csv
    scale_factor = 1
    axis = 1 # (0,1,2) => (x,y,z)
  [../]
  [./pressure_ramp]              # reads and interpolates input data defining amplitude curve for fill gas pressure
    type = PiecewiseLinear
    x = '-200 0'
    y = '0 1'
  [../]
  [./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_relocation_strain 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_strain clad_irradiation_strain'
    generate_output = 'vonmises_stress effective_creep_strain stress_xx stress_yy stress_zz'
    save_in = 'saved_x saved_y'
  [../]
[]

[Kernels]
  [./gravity]       # body force term in stress equilibrium equation
    type = Gravity
    variable = disp_y
    value = -9.81
  [../]
  [./heat]         # gradient term in heat conduction equation
    type = HeatConduction
    variable = temp
    save_in = saved_t
  [../]
  [./heat_ie]       # time term in heat conduction equation
    type = HeatConductionTimeDerivative
    variable = temp
    save_in = saved_t
  [../]
  [./heat_source]  # source term in heat conduction equation
    type = NeutronHeatSource
    variable = temp
    save_in = saved_t
    block = pellet_type_1 # fission rate applied to the fuel (block 2) only
    burnup_function = burnup
  [../]
[]

[Burnup]
  [./burnup]
    block = pellet_type_1
    rod_ave_lin_pow = power_history          # using the power function defined above
    axial_power_profile = axial_peaking_factors     # using the axial power profile function defined above
    num_radial = 80
    num_axial = 11
    a_lower = 0.00324    # mesh dependent!
    a_upper = 0.12184    # mesh dependent!
    fuel_inner_radius = 0
    fuel_outer_radius = .0041
    fuel_volume_ratio = 0.987775 # for use with dished pellets (ratio of actual volume to cylinder volume)
    order = CONSTANT
    family = MONOMIAL
    RPF = RPF
    i_enrich = '0.05 0.95 0 0 0 0'

    #N235 = N235 # Activate to write N235 concentration to output file
    #N238 = N238 # Activate to write N238 concentration to output file
    #N239 = N239 # Activate to write N239 concentration to output file
    #N240 = N240 # Activate to write N240 concentration to output file
    #N241 = N241 # Activate to write N241 concentration to output file
    #N242 = N242 # Activate to write N242 concentration to output file
  [../]
[]

[AuxKernels]
  # Define auxilliary kernels for each of the aux variables
  [./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
  [../]
  [./creep_strain_rate]
    type = MaterialRealAux
    property = creep_rate
    variable = creep_strain_rate
    block = clad
    execute_on = timestep_end
  [../]
  [./conductance]
    type = MaterialRealAux
    property = gap_conductance
    variable = gap_cond
    boundary = 10
    execute_on = 'linear'
  [../]
  [./coolant_htc]
    type = MaterialRealAux
    property = coolant_channel_htc
    variable = coolant_htc
    boundary = '1 2 3'
    execute_on = 'linear'
  [../]
[]

[Contact]
  # Define mechanical contact between the fuel (sideset=10) and the clad (sideset=5)
  [./pellet_clad_mechanical]
    master = 5
    slave = 10
    system = Constraint
    formulation = kinematic
    model = frictionless
    penalty = 1e7
  [../]
[]

[ThermalContact]
  # Define thermal contact between the fuel (sideset=10) and the clad (sideset=5)
  [./thermal_contact]
    type = GapHeatTransferLWR
    variable = temp
    master = 5
    slave = 10
    initial_moles = initial_moles       # coupling to a postprocessor which supplies the initial plenum/gap gas mass
    gas_released = fis_gas_released     # coupling to a postprocessor which supplies the fission gas addition
    contact_pressure = contact_pressure
    quadrature = true
  [../]
[]

[BCs]
# Define boundary conditions
  [./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] # apply coolant pressure on clad outer walls
    [./coolantPressure]
      boundary = '1 2 3'
      factor = 15.5e6
      function = pressure_ramp   # use the pressure_ramp function defined above
    [../]
  [../]
  [./PlenumPressure] # apply plenum pressure on clad inner walls and pellet surfaces
    [./plenumPressure]
      boundary = 9
      initial_pressure = 2.0e6
      startup_time = 0
      R = 8.3143
      output_initial_moles = initial_moles       # coupling to post processor to get initial fill gas mass
      temperature = ave_temp_interior            # coupling to post processor to get gas temperature approximation
      volume = gas_volume                        # coupling to post processor to get gas volume
      material_input = fis_gas_released          # coupling to post processor to get fission gas added
      output = plenum_pressure                   # coupling to post processor to output plenum/gap pressure
    [../]
  [../]
[]

[CoolantChannel]
  [./convective_clad_surface] # apply convective boundary to clad outer 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
  [../]
[]

[NuclearMaterials]
  [./UO2]
    [./fuel]
      block = pellet_type_1
      temperature = temp
      stress_free_temperature = 580.0
      diameter = 0.0082
      burnup_relocation_stop = 0.03
      thermal_model = NIFR
      grain_radius = grain_radius
      gbs_model = true
    [../]
  [../]
  [./ZirconiumAlloy]
    [./clad]
      block = clad
      temperature = temp
      stress_free_temperature = 295.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 = '-snes_ksp_ew'
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
  petsc_options_value = 'lu       superlu_dist'

  line_search = 'none'

  l_max_its = 50
  l_tol = 8e-3
  nl_max_its = 15
  nl_rel_tol = 1e-4
  nl_abs_tol = 1e-10

  start_time = -200
  n_startup_steps = 1
  end_time = 8.0e7

  dtmax = 2e6
  dtmin = 1

  [./TimeStepper]
    type = IterationAdaptiveDT
    dt = 2e2
    optimal_iterations = 8
    iteration_window = 2
    linear_iteration_ratio = 100
    growth_factor = 2
    cutback_factor = .5
  [../]
  [./Quadrature]
    order = THIRD
    side_order = FIFTH
  [../]
[]

[Postprocessors]
  # Define postprocessors (some are required as specified above; others are optional; many others are available)
  [./ave_temp_interior]            # average temperature of the cladding interior and all pellet exteriors
    type = SideAverageValue
    boundary = 9
    variable = temp
    execute_on = 'initial linear'
  [../]
  #[./centerline_temp]
  #  type = SideAverageValue
  #  boundary = 12
  #  variable = temp
  #  execute_on = linear
  #[../]
  [./clad_inner_vol]              # volume inside of cladding
    type = InternalVolume
    boundary = 7
    #outputs = exodus
    execute_on = 'initial timestep_end'
  [../]
  [./pellet_volume]               # fuel pellet total volume
    type = InternalVolume
    boundary = 8
    #outputs = exodus
    execute_on = 'initial timestep_end'
  [../]
  [./avg_clad_temp]               # average temperature of cladding interior
    type = SideAverageValue
    boundary = 7
    variable = temp
    execute_on = 'initial linear'
  [../]
  [./ave_fuel_temp]
    type = ElementAverageValue
    block = pellet_type_1
    variable = temp
    execute_on = 'initial linear'
  [../]
  [./fis_gas_produced]           # fission gas produced (moles)
    type = ElementIntegralFisGasGeneratedSifgrs
    variable = temp
    block = pellet_type_1
    execute_on = 'linear'
  [../]
  [./fis_gas_released]           # fission gas released to plenum (moles)
    type = ElementIntegralFisGasReleasedSifgrs
    variable = temp
    block = pellet_type_1
    execute_on = 'linear'
  [../]
  [./fis_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    variable = temp
    block = pellet_type_1
    outputs = exodus
    execute_on = 'linear'
  [../]
  [./fis_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    variable = temp
    block = pellet_type_1
    outputs = exodus
    execute_on = 'linear'
  [../]
  [./fisson_gas_release]
    type = FGRPercent
    value1 = fis_gas_released
    value2 = fis_gas_produced
    execute_on = 'linear'
  [../]

  [./gas_volume]
    type = InternalVolume
    boundary = 9
    execute_on = 'initial linear'
  [../]
  [./flux_from_clad]           # area integrated heat flux from the cladding
    type = SideFluxIntegral
    variable = temp
    boundary = 5
    diffusivity = thermal_conductivity
  [../]
  [./flux_from_fuel]          # area integrated heat flux from the fuel
    type = SideFluxIntegral
    variable = temp
    boundary = 10
    diffusivity = thermal_conductivity
  [../]

  [./_dt]                     # time step
    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 # rod height
  [../]

  [./mid_penetration]
    type = NodalVariableValue
    nodeid = 3781 #!!Mesh dependent!!
    variable = penetration
  [../]
  [./central_fuel_temp]
    type = NodalVariableValue
    variable = temp
    nodeid = 3781 # !! Mesh dependent
  [../]
  [./max_fuel_temp]
    type = NodalExtremeValue
    block = pellet_type_1
    value_type = max
    variable = temp
  [../]
  [./max_clad_temp]
    type = NodalExtremeValue
    block = clad
    value_type = max
    variable = temp
  [../]
  [./average_vonMises_fuel]
    type = ElementAverageValue
    variable = vonmises_stress
    block = pellet_type_1
  [../]
  [./average_vonMises_clad]
    type = ElementAverageValue
    variable = vonmises_stress
    block = clad
  [../]

  [./effective_creep_strain]
    type = ElementAverageValue
    block = clad
    variable = effective_creep_strain
  [../]
  [./effective_creep_strain_rate]
    type = ElementAverageValue
    block = clad
    variable = creep_strain_rate
  [../]
[]

[VectorPostprocessors]
  [./clad_dia]
    type = NodalValueSampler
    variable = disp_x
    boundary = 2
    sort_by = y
    outputs = 'outfile_clad_radial_displacement'
  [../]
  [./pellet_dia]
    type = NodalValueSampler
    variable = disp_x
    boundary = 10
    sort_by = y
    outputs = 'outfile_fuel_radial_displacement'
  [../]
[]

[Outputs]
  # Define output file(s)
  exodus = true
  color = false
  csv = true
  [./console]
    type = Console
    perf_log = true
    max_rows = 25
  [../]
  [./outfile_clad_radial_displacement]
    type = CSV
    execute_on = 'FINAL'
  [../]
  [./outfile_fuel_radial_displacement]
    type = CSV
    execute_on = 'FINAL'
  [../]
[]
(examples/tensor_mechanics/2D-RZ_rodlet_10pellets/2Dtm_discrete_finiteStrain_action.i)