TRISO Fuel Tutorial

# This example is 2D-RZ analysis of a TRISO fuel particle. Fully coupled
# heat transfer and solid mechanics, plus diffusion of the fission product
# species cesium (Cs) are simulated.  The mesh includes contact surfaces
# between the buffer and IPyC layers to facilitate a gap opening between
# these layers. These surfaces are initially in mechanical contact but
# are assumed to have no strength in tension. A coarse mesh is used to
# provide a short run time.

# The calculation simulates fuel-life in three steps. The first step is an
# irradiation period, where constant power and a fixed particle surface
# temperature (1500 K) are assumed over a lifetime of 76 Ms (2.4 yrs).
# For the second step, fuel removal and storage are simulated by setting
# the reactor power and Cs source terms to zero, reducing the particle
# surface temperature to ambient (300 K), and then holding it at ambient
# temperature for 100 days. A third and final step simulates accident
# behavior by increasing the particle surface temperature from ambient
# to 2073 K over 2 hrs, and then holding it at this elevated temperature
# for an additional 200 hrs. At the particle outer boundary, the Cs
# concentration is held at zero and the pressure at ambient during the
# entire simulation. The particle is assumed to be stress-free at an
# initial temperature of 1500 K.
#
# Details about this simulation are given in Section 4 of the following
# article: J. D. Hales, R. L. Williamson, S. R. Novascone, D. M. Perez,
# B. W. Spencer and G. Pastore, "Multidimensional multiphysics simulation
# of TRISO particle fuel", Journal on Nuclear Materials, Vol. 443, p. 531,
# 2013.

[GlobalParams]
  density = 11000.0          # kg/m^3
  order = SECOND
  family = LAGRANGE
  disp_x = disp_x
  disp_y = disp_y
  displacements = 'disp_x disp_y'
[]

[Mesh]
  file = triso2Dmed.e
[]

[Problem]
  coord_type = RZ
[]

[Variables]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./temp]
    initial_condition = 1500.0
  [../]
  [./conc]
    initial_condition = 0.0
    scaling = 1e18
  [../]
[]

[AuxVariables]
  [./fluence]
  [../]
  [./fast_neutron_flux]
  [../]
  [./fission_rate]
    block = fuel
  [../]
  [./burnup]
    block = fuel
  [../]
  [./stress_xx]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_yy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_zz]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_xy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_yz]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./stress_xz]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./hydrostatic_stress]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./gap_condSlave]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./creep_xx]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./creep_yy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./creep_zz]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[Functions]
  [./power_history]
    type = PiecewiseLinear
    x = '0 76e6 76.001e6'
    y = '1 1     0'
  [../]
  [./temp_bc]
    type = PiecewiseLinear
    x = '0    76e6  76.001e6 84.641e6 84.6482e6'
    y = '1500 1500  300  300     2073'
  [../]
  [./k_function]
    type = PiecewiseLinear
    x = '0     200e6'
    y = '4e-37 4e-37'
  [../]
  [./d1_function]
    type = ParsedFunction
    value = 'exp(t/4.5e25)'
  [../]
  [./d_gap]
    type = PiecewiseLinear
    x = '1500  2100'
    y = '1e-14 1e-12'
  [../]
[]

[SolidMechanics]
  [./solid]
    temp = temp
    disp_r = disp_x
    disp_z = disp_y
  [../]
[]

[Kernels]
  [./heat_ie]
    type = HeatConductionTimeDerivative
    variable = temp
  [../]
  [./heat]
    type = HeatConduction
    variable = temp
  [../]
  [./heat_source]
     type = NeutronHeatSource
     variable = temp
     block = fuel
     energy_per_fission = 3.2e-11        # units of J/fission
     fission_rate = fission_rate
  [../]
  [./mass_ie]
    type = TimeDerivative
    variable = conc
  [../]
  [./mass]
    type = ArrheniusDiffusion
    variable = conc
  [../]
  [./mass_source]
    type = BodyForce
    variable = conc
    function = power_history
    value = 1.22e-5     # units of moles/m**3-s
    block = fuel
  [../]
  [./mass_decay]
    type = Decay
    variable = conc
  [../]
[]

[AuxKernels]
  [./fast_neutron_flux]
    type = FastNeutronFluxAux
    variable = fast_neutron_flux
    function = power_history
    factor = 5e17
    execute_on = timestep_begin
  [../]
  [./fluence]
    type = FastNeutronFluenceAux
    variable = fluence
    fast_neutron_flux = fast_neutron_flux
    execute_on = timestep_begin
  [../]
  [./fission_rate]
    type = FissionRateAux
    variable = fission_rate
    block = fuel
    function = power_history
    value = 3.89e19
    execute_on = timestep_begin
  [../]
  [./burnup]
    type = BurnupAux
    variable = burnup
    block = fuel
    fission_rate = fission_rate
    molecular_weight = 0.270  # units of kg/mole
    execute_on = timestep_begin
  [../]
  [./stress_xx]
    type = MaterialTensorAux
    tensor = stress
    variable = stress_xx
    index = 0
    execute_on = timestep_end
  [../]
  [./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
  [../]
  [./stress_xy]
    type = MaterialTensorAux
    tensor = stress
    variable = stress_xy
    index = 3
    execute_on = timestep_end
  [../]
  [./stress_yz]
    type = MaterialTensorAux
    tensor = stress
    variable = stress_yz
    index = 4
    execute_on = timestep_end
  [../]
  [./stress_xz]
    type = MaterialTensorAux
    tensor = stress
    variable = stress_xz
    index = 5
    execute_on = timestep_end
  [../]
  [./hydrostatic_stress]
    type = MaterialTensorAux
    tensor = stress
    variable = hydrostatic_stress
    quantity = hydrostatic
    execute_on = timestep_end
  [../]
  [./creep_xx]
    type = MaterialTensorAux
    tensor = creep_strain
    variable = creep_xx
    index = 0
    block = 'buffer IPyC SiC OPyC'
    execute_on = timestep_end
  [../]
  [./creep_yy]
    type = MaterialTensorAux
    tensor = creep_strain
    variable = creep_yy
    index = 1
    block = 'buffer IPyC SiC OPyC'
    execute_on = timestep_end
  [../]
  [./creep_zz]
    type = MaterialTensorAux
    tensor = creep_strain
    variable = creep_zz
    index = 2
    block = 'buffer IPyC SiC OPyC'
    execute_on = timestep_end
  [../]
  [./conductanceSlave]
    type = MaterialRealAux
    property = gap_conductance
    variable = gap_condSlave
    boundary = BufferGapBndry
    execute_on = linear
  [../]
[]

[Contact]
  [./pellet_clad_mechanical]
    master = 15
    slave = 17
    penalty = 1e5
    model = frictionless
    formulation = penalty
  [../]
[]

[ThermalContact]
  [./thermal_contact]
    type = GapHeatTransferLWR
    variable = temp
    master = 15
    slave = 17
    initial_moles = initial_moles                      # coupling to a postprocessor which supplies the initial plenum/gap gas mass
    gas_released = 'fis_gas_released co_production'    # coupling to postprocessors which supply the fission gas addition, co addition
    gas_released_fractions = '0 0 0.153 0.847 0 0 0 0 0 0
                              0 0 0     0     0 0 0 1 0 0'
    tangential_tolerance = 1e-6
    #quadrature = true
  [../]
  [./cesium_contact]
    type = GapHeatTransfer
    variable = conc
    master = 15
    slave = 17
    tangential_tolerance = 1e-6
    gap_conductivity_function = d_gap
    gap_conductivity_function_variable = temp
    appended_property_name = _conc
    #quadrature = true
  [../]
[]

[BCs]
# pin particle along symmetry planes
  [./no_disp_x]
    type = DirichletBC
    variable = disp_x
    boundary = xzero
    value = 0.0
  [../]
  [./no_disp_y]
    type = DirichletBC
    variable = disp_y
    boundary = yzero #1
    value = 0.0
  [../]
# fix temperature on free surface
  [./freesurf_temp]
    type = FunctionDirichletBC
    variable = temp
    boundary = exterior
    function = temp_bc
  [../]
# fix concentration on free surface
  [./freesurf_conc]
    type = PresetBC
    variable = conc
    boundary = exterior
    value = 0.0
  [../]
  [./PlenumPressure] #  apply plenum pressure on clad inner walls and pellet surfaces
    [./plenumPressure]
      boundary = BufferGapVol
      initial_pressure = 0
      startup_time = 1.0e4
      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 = volumeGas                                 # coupling to post processor to get gas volume
      material_input = 'fis_gas_released co_production'  # coupling to post processor to get fission gas added, co added
      output = plenum_pressure                           # coupling to post processor to output plenum/gap pressure
    [../]
  [../]
[]

[Materials]
  [./radio_active_decay_constant]
    type = RadioActiveDecayConstant
    block = 'fuel buffer IPyC SiC OPyC'
    lambda = 7.297e-10   # units:(1/sec)  The constant for Cesium
  [../]
  [./fission_gas_release]              # Sifgrs fission gas release mode
    type = Sifgrs
    block = fuel
    temp = temp
    fission_rate = fission_rate        # coupling to fission_rate aux variable
    grain_radius_const = 5.0e-6
  [../]
  [./fuel_thermal]                       # temperature and burnup dependent thermal properties of UO2 (bison kernel)
    type = ThermalFuel
    model = 2
    block = fuel
    temp = temp
    burnup = burnup
    initial_porosity = 0.0
  [../]
  [./fuel_solid_mechanics_swelling]      # free expansion strains (swelling and densification) for UO2 (bison kernel)
    type = VSwellingUO2
    gas_swelling_type = MATPRO
    block = fuel
    temp = temp
    burnup = burnup
  [../]
  [./fuel_disp]                         # thermal and irradiation creep for UO2 (bison kernel)
    type = Elastic
    block = fuel
    disp_r = disp_x
    disp_z = disp_y
    temp = temp
    youngs_modulus = 2.2e11
    poissons_ratio = .345
    thermal_expansion = 10e-6
  [../]
  [./fuel_den]
    type = Density
    block = fuel
  [../]
  [./fuel_conc]
    type = ArrheniusDiffusionCoef
    block = fuel
    d1 = 5.6e-8               # m^2/s
    q1 = 209.0e+3             # J/mol
    d2 = 5.2e-4               # m^2/s
    q2 = 362.0e+3             # J/mol
    gas_constant = 8.3143     # J/K-mol
    temp = temp
  [../]
  [./PyCBuffer]
    type = PyCIrradiationStrain
    block = buffer
    fluence = fluence
    pyc_type = buffer
  [../]
  [./buffer_disp]
    type = CreepPyC
    block = buffer
    temp = temp
    youngs_modulus = 2.0e10
    poissons_ratio = 0.23
    thermal_expansion = 5.65e-6
    disp_r = disp_x
    disp_z = disp_y
    flux = fast_neutron_flux
    density = 1000.0                   #kg/m^3
  [../]
  [./buffer_temp]
    type = HeatConductionMaterial
    block = buffer
    thermal_conductivity = 0.5        # J/m-s-K
    specific_heat = 720.0             # J/kg-K
  [../]
  [./buffer_den]
    type = Density
    density = 1000.0                   #kg/m^3
    block = buffer
  [../]
  [./buffer_conc]
    type = ArrheniusDiffusionCoef
    block = buffer
    d1 = 1.0e-12                        # m^2/s
    q1 = 0.0
    d2 = 0.0
    q2 = 0.0
    gas_constant = 8.3143     # J/K-mol
    temp = temp
  [../]
  [./IPyC_densify]
    type = PyCIrradiationStrain
    block = IPyC
    fluence = fluence
    pyc_type = dense
  [../]
  [./IPyC_disp]
    type = CreepPyC
    block = IPyC
    temp = temp
    youngs_modulus = 4.74e10
    poissons_ratio = 0.23
    thermal_expansion = 5.65e-6
    disp_r = disp_x
    disp_z = disp_y
    flux = fast_neutron_flux
    density = 1900.0                        # kg/m^3
  [../]
  [./IPyC_temp]
    type = HeatConductionMaterial
    block = IPyC
    thermal_conductivity = 4.0          # J/m-s-K
    specific_heat = 720.0               # J/kg-K
  [../]
  [./IPyC_den]
    type = Density
    density = 1900.0                        # kg/m^3
    block = IPyC
  [../]
  [./IPyC_conc]
    type = ArrheniusDiffusionCoef
    block = IPyC
    d1 = 6.3e-8                   # m^2/s
    q1 = 222.0e+3                    # J/mol
    d2 = 0.0
    q2 = 0.0
    gas_constant = 8.3143     # J/K-mol
    temp = temp
  [../]
  [./SiC_disp]
    type = CreepSiC
    block = SiC
    disp_r = disp_x
    disp_z = disp_y
    temp = temp
    k_function = k_function
    fast_neutron_flux = fast_neutron_flux
    youngs_modulus = 3.4e11
    poissons_ratio = .13
    thermal_expansion = 4.9e-6
  [../]
  [./SiC_temp]
    type = HeatConductionMaterial
    block = SiC
    thermal_conductivity = 13.9          # J/m-s-K
    specific_heat = 620.0                # J/kg-K
  [../]
  [./SiC_den]
    type = Density
    density = 3180.0                     # kg/m^3
    block = SiC
  [../]
  [./SiC_conc]
    type = ArrheniusDiffusionCoef
    block = SiC
    d1 = 5.5e-14                  # m^2/s
    d1_function = d1_function
    d1_function_variable = fluence
    q1 = 125.0e+3                 # J/mol
    d2 = 1.6e-2                   # m^2/s
    q2 = 514.0e+3                 # J/mol
    gas_constant = 8.3143         # J/K-mol
    temp = temp
  [../]
  [./OPyC_densify]
    type = PyCIrradiationStrain
    block = OPyC
    fluence = fluence
    pyc_type = dense
  [../]
  [./OPyC_disp]
    type = CreepPyC
    block = OPyC
    temp = temp
    youngs_modulus = 4.74e10
    poissons_ratio = 0.23
    thermal_expansion = 5.65e-6
    disp_r = disp_x
    disp_z = disp_y
    flux = fast_neutron_flux
    density = 1900.0                         # kg/m^3
  [../]
  [./OPyC_temp]
    type = HeatConductionMaterial
    block = OPyC
    thermal_conductivity = 4.0                # J/m-s-K
    specific_heat = 720.0                     # J/kg-K
  [../]
  [./OPyC_den]
    type = Density
    density = 1900.0                         # kg/m^3
    block = OPyC
  [../]
  [./OPyC_conc]
    type = ArrheniusDiffusionCoef
    block = OPyC
    d1 = 6.3e-8                             # m^2/s
    q1 = 222.0e+3                           # J/mol
    d2 = 0.0
    q2 = 0.0
    gas_constant = 8.3143                   # J/K-mol
    temp = temp
  [../]
[]

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

[Debug]
  show_var_residual_norms = true
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options_iname = '-ksp_gmres_restart -pc_type -pc_hypre_type -pc_hypre_boomeramg_max_iter'
  petsc_options_value = '201                hypre    boomeramg      4'

  line_search = 'none'

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

  l_tol = 1e-3
  l_max_its = 50

  start_time = 0.0
  end_time = 85.3682e6
  dt = 20

  dtmax = 2e6
  dtmin = 1

  [./TimeStepper]
    type = IterationAdaptiveDT
    dt = 1.0e2
    optimal_iterations = 6
    growth_factor = 1.5
    linear_iteration_ratio = 100
    time_t = '0    76e6  76.001e6 84.641e6 84.6482e6'
    time_dt = '20  20     20      20       20'
  [../]

  [./Quadrature]
    order = THIRD
  [../]
[]

[Postprocessors]
  [./Cs_release]
    type = SideIntegralMassFlux
    variable = conc
    boundary = exterior
    execute_on = timestep_end
  [../]
  [./dt]
    type = TimestepSize
    execute_on = timestep_end
  [../]
  [./fis_gas_produced]           # fission gas produced (moles)
    type = ElementIntegralFisGasGeneratedSifgrs
    variable = temp
    block = fuel
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  [../]
  [./fis_gas_released]           # fission gas released to plenum (moles)
    type = ElementIntegralFisGasReleasedSifgrs
    variable = temp
    block = fuel
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  [../]
  [./volumeTotal]
    type = InternalVolume
    boundary = exterior
    execute_on = 'initial timestep_end'
  [../]
  [./volumeFuel]
    type = InternalVolume
    boundary = fuel
    execute_on = 'initial timestep_end'
  [../]
  [./volumeGas]
    type = InternalVolume
    boundary = BufferGapVol
    # ro = 3.125e-4
    # ri = 2.125e-4
    # vb = 4/3*pi*(ro^3-ri^3) = 8.76e-11
    # buffer density = 1000
    # PyC density    = 1900
    # fill ratio = 10/19
    # vb*10/19 = 4.6e-11
    # Must remove 4.6e-11 m^3 from the volume
    addition = -4.6e-11
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  [../]
  [./volumeBufferShell]
    type = InternalVolume
    boundary = BufferGapVol
    execute_on = 'initial timestep_end'
  [../]
  [./ave_temp_interior]
    type = SideAverageValue
    boundary = BufferGapVol
    variable = temp
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  [../]

# Postprocessors for CO production

  [./total_fission_rate]
    type = ElementIntegralPower
    variable = temp
    fission_rate = fission_rate
    block = fuel
    energy_per_fission = 1.0
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  [../]
  [./total_fissions]
    type = TotalVariableValue
    value = total_fission_rate
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  [../]
  [./avg_surface_temp]
    type = SideAverageValue
    variable = temp
    boundary = exterior
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  [../]
  [./time_int_surf_temp]
    type = TotalVariableValue
    value = avg_surface_temp
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  [../]
  [./co_production]
    type = CarbonMonoxideProduction
    total_fissions = total_fissions
    time_int_surf_temp = time_int_surf_temp
    execute_on = 'initial linear nonlinear timestep_begin timestep_end'
  [../]
  [./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
  [../]
[]

[Outputs]
  exodus = true
  [./console]
    type = Console
    perf_log = true
    max_rows = 25
  [../]
[]
(examples/triso_particle_2D-RZ/triso2D_accident.i)