3D Light Water Reactor Fuel Rod Tutorial

under construction

The current example uses the soon to be deprecated Solid Mechanics module, and will be updated to use a the newer Tensor Mechanics module.

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

[Mesh]
  file = smearedTest3.e
  patch_size = 5 #40
  patch_update_strategy = auto
  partitioner = centroid
  centroid_partitioner_direction = y
[]

[Problem]
  type = ReferenceResidualProblem
  solution_variables = 'disp_x disp_y disp_z temp'
  reference_residual_variables = 'saved_x saved_y saved_z saved_t'
[]

[Variables]
  # Define dependent variables and initial conditions
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./disp_z]
  [../]
  [./temp]
    initial_condition = 580     # set to coolant inlet temp
  [../]
[]

[AuxVariables]
  # Define auxilary variables
  [./saved_x]
  [../]
  [./saved_y]
  [../]
  [./saved_z]
  [../]
  [./saved_t]
  [../]
  [./fast_neutron_flux]
    block = clad
  [../]
  [./fast_neutron_fluence]
    block = clad
  [../]
  [./grain_radius]
    block = pellet_type_1
    initial_condition = 10e-6
  [../]
  [./stress_xx]      # stress aux variables are defined for output; this is a way to get integration point variables to the output file
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_yy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_zz]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./vonmises]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./creep_strain_hoop]
    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 = ParsedFunction
     value = '1'
#    type = PiecewiseBilinear
#    data_file = peakingfactors5.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'
  [../]
[]

[SolidMechanics]
  # Specify that we need solid mechanics (divergence of stress)
  [./solid]
    temp = temp
    save_in_disp_x = saved_x
    save_in_disp_y = saved_y
    save_in_disp_z = saved_z
  [../]
[]

[Kernels]
  # Define kernels for the various terms in the PDE system
  [./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
    block = pellet_type_1        # fission rate applied to the fuel (block 2) only
    fission_rate = fission_rate  # coupling to the fission_rate aux variable
    save_in = saved_t
  [../]
[]

[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 = 2.49e-3 #2.26e-3
    a_upper = 2.621e-2 #1.2086e-1
    fuel_inner_radius = 0
    fuel_outer_radius = .0041
    #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
    RPF = RPF
  [../]
[]

[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
  [../]
  [./stress_xx]               # computes stress components for output
    type = MaterialTensorAux
    tensor = stress
    variable = stress_xx
    index = 0
    execute_on = timestep_end     # for efficiency, only compute at the end of a timestep
  [../]
  [./stress_yy]
    type = MaterialTensorAux
    tensor = stress
    variable = stress_yy
    index = 1
    execute_on = timestep_end
  [../]
  [./stress_zz]
    type = MaterialTensorAux
    tensor = stress
    variable = stress_zz
    index = 2
    execute_on = timestep_end
  [../]
  [./vonmises]
    type = MaterialTensorAux
    tensor = stress
    variable = vonmises
    quantity = vonmises
    execute_on = timestep_end
  [../]
  [./creep_strain_hoop]
    type = MaterialTensorAux
    block = clad
    tensor = creep_strain
    variable = creep_strain_hoop
    quantity = Hoop
    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'
  [../]
[]

[Contact]
  # Define mechanical contact between the fuel (sideset=10) and the clad (sideset=5)
  [./pellet_clad_mechanical]
    master = 5
    slave = 10
#    formulation = penalty
    penalty = 1e8
    tangential_tolerance = 1e-4
    model = frictionless
    normal_smoothing_distance = 0.1
    system = constraint
  [../]
[]

[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
    tangential_tolerance = 1e-4
    contact_pressure = contact_pressure
    quadrature = true
    normal_smoothing_distance = 0.1
  [../]
[]

[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_z_all] # pin pellets and clad along axis of symmetry (y)
    type = DirichletBC
    variable = disp_z
    boundary = 13
    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
      displacements = 'disp_x disp_y disp_z'
    [../]
  [../]
[]

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

[Materials]
  # Define material behavior models and input material property data
  [./fuel_thermal]                       # temperature and burnup dependent thermal properties of UO2 (bison kernel)
    type = ThermalFuel
    block = pellet_type_1
    model = 2 # Fink-Lucuta
    temp = temp
    burnup = burnup
  [../]
  [./fuel_solid_mechanics_swelling]      # free expansion strains (swelling and densification) for UO2 (bison kernel)
    type = VSwellingUO2
    gas_swelling_type = MATPRO
    block = pellet_type_1
    temp = temp
    burnup = burnup
  [../]
  [./fuel_creep]                         # thermal and irradiation creep for UO2 (bison kernel)
    type = Elastic #CreepUO2
    block = pellet_type_1
    temp = temp
#    fission_rate = fission_rate
    youngs_modulus = 2.e11
    poissons_ratio = .345
    thermal_expansion = 10e-6
#    grain_radius = 10.0e-6
#    oxy_to_metal_ratio = 2.0
  [../]
  [./fuel_relocation]
    type = RelocationUO2
    block = pellet_type_1
    burnup = burnup
    diameter = 0.0082
    q = q
    gap = 160e-6 # diametral gap
#    burnup_relocation_stop = 0.0136
  [../]
  [./clad_thermal]                       # general thermal property input (elk kernel)
    type = HeatConductionMaterial
    block = clad
    thermal_conductivity = 16.0
    specific_heat = 330.0
  [../]
  [./clad_solid_mechanics]               # thermoelasticity and thermal and irradiation creep for Zr4 (bison kernel)
    type = MechZry #ThermalIrradiationCreepZr4
    block = clad
    temp = temp
    fast_neutron_flux = fast_neutron_flux
    fast_neutron_fluence = fast_neutron_fluence
    youngs_modulus = 7.5e10
    poissons_ratio = 0.3
    thermal_expansion = 5.0e-6
    model_thermal_expansion = false
    stress_free_temperature = 295.0
  [../]
#  [./clad_irrgrowth]
#    type = IrradiationGrowthZr4
#    block = clad
#    fast_neutron_fluence = fast_neutron_fluence
#  [../]
  [./fission_gas_release]              # Forsberg-Massih fission gas release mode
    type = Sifgrs
    block = pellet_type_1
    temp = temp
    fission_rate = fission_rate        # coupling to fission_rate aux variable
    grain_radius = grain_radius
    gbs_model = true
  [../]
  [./clad_density]
    type = Density
    block = clad
    density = 6551.0
  [../]
  [./fuel_density]
    type = Density
    block = pellet_type_1
  [../]
[]

[Dampers]
   [./limitT]
     type = MaxIncrement
     variable = temp
     max_increment = 50
   [../]
[]

[Executioner]

  # PETSC options:
  #   petsc_options
  #   petsc_options_iname
  #   petsc_options_value
  #
  # controls for linear iterations
  #   l_max_its
  #   l_tol
  #
  # controls for nonlinear iterations
  #   nl_max_its
  #   nl_rel_tol
  #   nl_abs_tol
  #
  # time control
  #   start_time
  #   dt
  #   optimal_iterations
  #   iteration_window
  #   linear_iteration_ratio

  type = Transient
  solve_type = 'PJFNK'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -ksp_gmres_restart -snes_ksp_ew_rtol0 -snes_ksp_ew_rtolmax -snes_ksp_ew_gamma -snes_ksp_ew_alpha -snes_ksp_ew_alpha2 -snes_ksp_ew_threshold'
  petsc_options_value = ' lu       superlu_dist                  51                              0.5                  0.9                  1                  2                   2                    0.1'

  line_search = 'none'

  l_max_its = 100
  l_tol = 8e-3

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

  start_time = -200
  end_time = 8.0e7

  dtmax = 2e6
  dtmin = 1

  [./TimeStepper]
    type = IterationAdaptiveDT
    dt = 2e2
    optimal_iterations = 8
    linear_iteration_ratio = 100
  [../]

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

[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'
  [../]
  [./clad_inner_vol]              # volume inside of cladding
    type = InternalVolume
    boundary = 7
    outputs = exodus
  [../]
  [./pellet_volume]               # fuel pellet total volume
    type = InternalVolume
    boundary = 8
    outputs = exodus
  [../]
  [./avg_clad_temp]               # average temperature of cladding interior
    type = SideAverageValue
    boundary = 7
    variable = temp
  [../]
  [./fis_gas_produced]           # fission gas produced (moles)
    type = ElementIntegralFisGasGeneratedSifgrs
    variable = temp
    block = pellet_type_1
  [../]
  [./fis_gas_released]           # fission gas released to plenum (moles)
    type = ElementIntegralFisGasReleasedSifgrs
    variable = temp
    block = pellet_type_1
  [../]
  [./fis_gas_grain]
    type = ElementIntegralFisGasGrainSifgrs
    variable = temp
    block = pellet_type_1
    outputs = exodus
  [../]
  [./fis_gas_boundary]
    type = ElementIntegralFisGasBoundarySifgrs
    variable = temp
    block = pellet_type_1
    outputs = exodus
  [../]
  [./gas_volume]                # 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
  [../]
  [./alive_time]
    type = PerformanceData
    event = ALIVE
  [../]
  [./tot_nonlin_it]
    type = CumulativeValuePostprocessor
    postprocessor = num_nonlin_it
  [../]
  [./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 # rod height (2 pellets)
  [../]
  [./average_fission_rate]
    type = AverageFissionRate
    rod_ave_lin_pow = power_history
  [../]
[]

[Outputs]
  # Define output file(s)
  exodus = true
  [./console]
    type = Console
    perf_log = true
    max_rows = 25
  [../]
[]
(examples/3D_rodlet_3pellets/smearedTest3D.i)