BISON TRISO Workshop

Implicit, parallel, fully-coupled nuclear fuel performance analysis

Computational Mechanics and Materials Department
Idaho National Laboratory

BISON History

BISON Fuel Performance Code

Finite element-based engineering scale fuel performance code
Solves the fully-coupled thermomechanics and species diffusion equations in 1D, 1.5D, 2D axisymmetric or plane-strain, or full 3D
Used for LWR, ATF, TRISO, and metallic fuels
Applicable to both steady and transient operations and includes LOCA and RIA capability for LWR fuel
Readily coupled to lower length scale material models
Designed for efficient use on parallel computers
Development follows NQA-1 process

2D axisymmetric (or 1.5D)

2D plane strain

$~$

BISON Governing Equations

Energy conservation (transient heat conduction with fission source)
- $\rho c_p \frac{\partial T}{\partial t}=\nabla\cdot (k\nabla T)+E_f\dot{F}$
- $E_f$ is energy per fission. $\dot{F}(x, t)$ is fission density rate.
Species conservation (transient oxygen or fission product diffusion with radioactive decay)
- $\frac{\partial C}{\partial t}=\nabla\cdot D(\nabla C-\frac{Q}{R T^2}\nabla T) - \lambda C + S$
- Equation includes Fickian diffusion, Soret diffusion, radioactive decay, and a source term.
Momentum conservation (Cauchy's equation of equilibrium)
- $\mathbf{\nabla\cdot T+\rho f = 0}$
- $\mathbf{T}$ is the stress tensor, and $\mathbf{f}$ is a body force.

Fuel Behavior: Introduction

At beginning of life, a fuel element is quite simple...

Nakajima et al., Nuc. Eng. Des., 148, 41 (1994)

$~$

but irradiation brings about substantial complexity...

\Longrightarrow

Michel et al., Eng. Frac. Mech., 75, 3581 (2008)

Fuel Fracture

Olander, p. 584 (1978)

Multidimensional Contact
and Deformation

Olander, p. 584 (1978)

Fission Gas

Bentejac et al., PCI Seminar (2004)

Stress Corrosion
Cracking Cladding
Failure

Fuel Behavior Modeling: Coupled Multiphysics

Multiphysics

Fully coupled nonlinear thermo-mechanics
Multiple species diffusion
Neutronics
Thermal-hydraulics
Chemistry

Multi-space scale

Important physics at the atomistic and micro-structural levels
Practical engineering simulations require the continuum level

Multi-time scale

Steady operation ( $\Delta t > 1$ week)
Power ramps/accidents
( $\Delta t < 0.1$ s)

MOOSE, BISON, and Marmot

MOOSE, BISON, and Marmot provide an advanced multidimensional, multiphysics, multiscale fuel performance capability

Atomistic/Mesoscale Material Model Development

Predicts microstructure evolution in fuel and cladding
Used with atomistic methods to develop multiscale material models

$~$

Advanced Multidimensional Fuel Performance Code

Models a wide variety of fuel types and geometries at an engineering scale
Applicable for steady, transient and accident conditions

Simulation framework allowing rapid development of FEM-based applications

BISON test/tests Directory

ADMetallicFuelWastage/
Al2O3/
GrainRadiusPorosity_test/
GraphiteMatrixElasticityTensor/
GraphiteMatrixSpecificHeat/
GraphiteMatrixThermalConductivity/
MetallicFuelWastage/
OxideEnergyDeposition/
SiCPdPenetration/
ThermalFuel_error_messages/
ad_arrhenius_material_property/
ad_d9_thermal/
ad_ht9_thermal/
ad_mc_thermal/
ad_ss316_thermal/
ad_upuzr_burnup/
ad_upuzr_fast_neutron_flux/
ad_upuzr_fission_gas_release/
ad_upuzr_fission_rate/
ad_upuzr_sodium_logging/
ad_upuzr_thermal/
anisotropic_swelling/
arrhenius_diffusion_coef/
arrhenius_material_property/
average_axial_position/
axial_relocation/
burnup_action/
carbon_monoxide_production/
check_error/
circular_cross_section_mesh/
constitutive_heat_conduction/
convective_heat_transfer/
coolant_channel_model/
creep_SiC/
creep_U10Mo/
creep_mox/
creep_uo2/
cumulative_damage_index/
decay_heating/
diffusion_limited_reaction/
dislocation_density/
dryCask/
effective_burnup_aux/
element_integral_power/
example_problem_test/
fast_neutron_flux/
fcci_ht9/
fecral_oxidation/
fgr_fraction/
fgr_percent/
fgr_upuzr/
fill_gas_thermal_conductivity/
fission_gas_1d/
sifgrs/uo2/

sifgrs/u3si2/
fission_gas_release_formas/
fission_rate_LWR/
fission_rate_MOX/
fission_rate_axial/
fission_rate_from_power_density/
fission_rate_heat_source/
fuelrodlinevaluesampler/
gamma_heating/
gap_heat_transfer/
gap_heat_transfer_fission/
gap_heat_transfer_htonly/
gap_heat_transfer_mixedgas/
gap_heat_transfer_radiation/
gap_jump_distance/
gap_perfect_transfer/
generic_material_failure/
grain_radius_aux/
hydride/
hydrogen/
ifba_he_production/
irradiation_growth/
irradiation_growth_Zr4/
layered2D/
layered_1D/
mechTests/
mechZry/
mechanical_uo2/
meso_thcond_test/
monolithicSiCThermal/
mox_oxygen_to_metal_ratio/
mox_pore_velocity/
oxygen_aux/
oxygen_transport/
partial_sum_heat_flux/
percolation/
performance_outputs_action/
phase_transition_zircaloy/
phase_upuzr/
plate_mesh/
plenum_pressure/
plenum_temp/
power_peaking_function/
radial_avg_fuel_enthalpy/
radial_crack_counter/
radial_power_factor/
radioactive_decay/
radius_aux/
relocation_UO2/
side_ave_incr_tensor_component/
side_int_var_incr_postprocess/
side_integral_mass_flux/
fuel_pin_mesh/
fuel_pin_mesh_fipd/

fuel_pin_mesh_generator/
fuel_pin_mesh_generator_fipd/
sodium_coolant_channel/
solid_mechanics_deprecated/
species_source/
stan_neumann/
standard_lwr_outputs_action/
submodel_end_bc/
temperature_jump_distance/
solid_mechanics/
thermalChromium/
thermalCompositeSiC/
thermalD9/
thermalFastMOX/
thermalFeCrAl/
thermalFuel_Amaya/
thermalFuel_Duriez/
thermalFuel_FinkLucuta/
thermalFuel_Halden/
thermalFuel_HaldenMOX/
thermalFuel_HaldenUO2/
thermalFuel_NFIR/
thermalFuel_NFImod/
thermalFuel_Ronchi/
thermalFuel_Staicu/
thermalFuel_Toptan/
thermalFuel_rimLayer/
thermalHT9/
thermalMAMOX/
thermalMOX/
thermalNa/
thermalSilicideFuel/
thermalTests/
thermalUO2/
thermalZrO2/
thermalZry/
thermal_accommodation_coeff/
thermirrad_creep_zr42/
thermo_mech_oxygen/
triso/
triso_failure/
un_swelling/
upuzr_burnup/
upuzr_dictra/
upuzr_diffusivity/
upuzr_fast_neutron_flux/
upuzr_fission_rate/
upuzr_phase_lookup/
void_volume/
zirconium_diffusion/
zrdiffusivity_upuzr/
zrh_formation/
zry_plasticity/

BISON has about 1800 regression tests

BISON examples Directory

1.5D_restart/
1.5D_rodlet_10pellets/
2D-RZ_rodlet_10pellets/
2D_plane_strain_rod/
3D_rodlet_3pellets/
TRISO/
accident_tolerant_fuel/

$~$

axial_relocation/
fast_mox_sifgrs/
hydride_rim/
metal_fuel/
mox_fuel/
multiapp/

$~$

non-cylindrical_fuel/
percolation/
pore_migration/
restart/
spent_fuel/
temperature_tables/

BISON examples/TRISO/full_particle/1D Directory

examples
full_particle_1D.i
gold
tests

Example BISON Source File

src/materials/solid_mechanics/UCOVolumetricSwellingEigenstrain.C

/*************************************************/
/*           DO NOT MODIFY THIS HEADER           */
/*                                               */
/*                     BISON                     */
/*                                               */
/*    (c) 2015 Battelle Energy Alliance, LLC     */
/*            ALL RIGHTS RESERVED                */
/*                                               */
/*   Prepared by Battelle Energy Alliance, LLC   */
/*     Under Contract No. DE-AC07-05ID14517      */
/*     With the U. S. Department of Energy       */
/*                                               */
/*     See COPYRIGHT for full restrictions       */
/*************************************************/

#include "UCOVolumetricSwellingEigenstrain.h"

registerMooseObject("BisonApp", UCOVolumetricSwellingEigenstrain);
registerMooseObject("BisonApp", ADUCOVolumetricSwellingEigenstrain);

template <bool is_ad>
InputParameters
UCOVolumetricSwellingEigenstrainTempl<is_ad>::validParams()
{
  InputParameters params = ComputeEigenstrainBaseTempl<is_ad>::validParams();
  params.addClassDescription(
      "Computes fission-induced swelling (percent per percent FIMA) for UCO.");
  params.addParam<Real>("swelling_rate", 2.9, "Swelling rate (%).");
  params.addParam<Real>("swelling_scale_factor", 1.0, "Multiplier for UCO swelling");
  return params;
}

template <bool is_ad>
UCOVolumetricSwellingEigenstrainTempl<is_ad>::UCOVolumetricSwellingEigenstrainTempl(
    const InputParameters & parameters)
  : ComputeEigenstrainBaseTempl<is_ad>(parameters),
    _swelling_rate(parameters.get<Real>("swelling_rate")),
    _swelling_scale_factor(parameters.get<Real>("swelling_scale_factor")),
    _burnup(this->template getGenericMaterialProperty<Real, is_ad>("burnup")),
    _swelling(this->template declareGenericProperty<Real, is_ad>("swelling"))
{
}

template <bool is_ad>
void
UCOVolumetricSwellingEigenstrainTempl<is_ad>::initQpStatefulProperties()
{
  _swelling[_qp] = 0;
  ComputeEigenstrainBaseTempl<is_ad>::initQpStatefulProperties();
}

template <bool is_ad>
void
UCOVolumetricSwellingEigenstrainTempl<is_ad>::computeQpEigenstrain()
{
  GenericReal<is_ad> volumetric_swelling_strain =
      _swelling_scale_factor * _swelling_rate * _burnup[_qp];
  GenericReal<is_ad> strain_component =
      this->computeVolumetricStrainComponent(volumetric_swelling_strain);
  _swelling[_qp] = volumetric_swelling_strain;

  _eigenstrain[_qp].zero();
  _eigenstrain[_qp].addIa(strain_component);
}

template class UCOVolumetricSwellingEigenstrainTempl<false>;
template class UCOVolumetricSwellingEigenstrainTempl<true>;

(src/materials/solid_mechanics/UCOVolumetricSwellingEigenstrain.C)

Example BISON Test Continued

test/tests/triso/UCOVolumetricSwellingEigenstrain/UCOVolumetricSwellingEigenstrain.i

# UCO fission-induced swelling
# The geometry is a unit cube made of UCO subject to burnup-induced swelling.
#
# The swelling is simply 2.9 * burnup.  Burnup is ramped from 0 to 0.125.
# Thus, swelling increases to 0.3625.  The final volume is 1.36256.

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
  volumetric_locking_correction = true
[]

[Mesh]
  coord_type = XYZ
  [mesh]
    type = GeneratedMeshGenerator
    dim = 3
    xmax = 1.0
    ymax = 1.0
    zmax = 1.0
  []
[]

[AuxVariables]
  [swelling]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Functions]
  [burnup]
    type = PiecewiseLinear
    x = '0 10'
    y = '0 0.125'
  []
[]

[Physics/SolidMechanics/QuasiStatic]
  add_variables = true
  [stuff]
    block = '0'
    strain = FINITE
    eigenstrain_names = 'UCO_swelling_eigenstrain'
    generate_output = 'vonmises_stress stress_xx stress_yy stress_zz'
  []
[]

[AuxKernels]
  [swelling]
    type = MaterialRealAux
    variable = swelling
    property = swelling
    block = '0'
    execute_on = linear
  []
[]

[BCs]
  [no_x]
    type = DirichletBC
    variable = disp_x
    boundary = 'left'
    value = 0
  []

  [no_y]
    type = DirichletBC
    variable = disp_y
    boundary = 'bottom'
    value = 0
  []

  [no_z]
    type = DirichletBC
    variable = disp_z
    boundary = 'back'
    value = 0
  []
[]

[Materials]

  [burnup]
    type = GenericFunctionMaterial
    prop_names = burnup
    prop_values = burnup
  []

  [UCO_VolumetricSwellingEigenstrain]
    type = UCOVolumetricSwellingEigenstrain
    eigenstrain_name = UCO_swelling_eigenstrain
  []

  [stress]
    type = ComputeFiniteStrainElasticStress
  []

  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 10
    poissons_ratio = 0.4
  []

  [UCO_density]
    type = Density
    block = '0'
    density = 11250.0
  []

[]

[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 = 50
  l_tol = 1e-2

  nl_max_its = 150
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-10

  start_time = 0.0
  end_time = 10.0
  dt = 1.0

[]

[Postprocessors]
  [_dt]
    type = TimestepSize
  []

  [volume_UCO]
    type = VolumePostprocessor
    use_displaced_mesh = true
  []

  [swelling_UCO_max]
    type = ElementExtremeValue
    value_type = 'max'
    variable = swelling
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  csv = true
[]

(test/tests/triso/UCOVolumetricSwellingEigenstrain/UCOVolumetricSwellingEigenstrain.i)

Case	Geometry	Description
1	SiC layer	Elastic only
2	IPyC layer	Elastic only
3	IPyC/SiC layer	Elastic with no fluence
4a	IPyC/SiC layer	Swelling and no creep
4b	IPyC/SiC layer	Creep and no swelling
4c	IPyC/SiC layer	Creep and swelling
4d	IPyC/SiC layer	Creep- and fluence-dependent swelling
5	TRISO	350 $\mu$ m kernel, real conditions
6	TRISO	500 $\mu$ m kernel, real conditions
7	TRISO	Same as 6 with high BAF pyC
8	TRISO	Same as 6 with cyclic temperature
10	HFR-K3	10% FIMA, $5.3\times10^{-25}$ n/m $^2$ fluence
11	HFR-P4	10% FIMA, $7.2\times10^{-25}$ n/m $^2$ fluence

TRISO Thermal Models

Thermal modeling for TRISO fuel follows thet same pattern as for any other fuel:

\rho c_p \frac{\partial T}{\partial t} = \nabla \cdot \left(k \nabla T \right) + E_f \dot{F}

That is, we need to
1. Define density, specific heat, and thermal conductivity, which may be functions of temperature or other parameters.
2. Ef is the energy released in a single fission event and F dot is the volumetric fission rate.
3. Invoke the heat conduction and heat conduction time derivative kernels.

\rho c_p \frac{\partial T}{\partial t} = \nabla \cdot \left(k \nabla T \right) + E_f \dot{F}

$~$

[heat_ie]
  type = HeatConductionTimeDerivative
  variable = temperature
  extra_vector_tags = 'ref'
[]
[heat]
  type = HeatConduction
  variable = temperature
  extra_vector_tags = 'ref'
[]
[heat_source]
  type = NeutronHeatSource
  variable = temperature
  block = fuel
  fission_rate = fission_rate
  extra_vector_tags = 'ref'
[]

(examples/TRISO/parfume/parfume.i)

TRISO Mechanical Models

Mechanical modeling for TRISO fuel follows the same pattern as for any other fuel:

\nabla \cdot \bf{T} + \rho \bf{f} = 0

This is more involved than thermal modeling.
1. Define constitutive response (define, e.g., an elasticity tensor and an object to convert strain to stress). 2. Define so-called eigenstrains (thermal strain, irradiation strain)
The code provides input shortcuts to invoke the computation of strain and the divergence of stress. This shortcut may or may not be used.

Creep + irradiation Strain

$~$

[Materials]
  [IPyC_elasticity_tensor]
    type = PyCElasticityTensor
    block = IPyC
    temperature = temperature
    initial_BAF = 1.045
  []

  [IPyC_stress]
    type = PyCCEGACreep
    block = IPyC
    temperature = temperature
  []

  [IPyC_IIDC_strain]
    type = PyCCEGAIrradiationEigenstrain
    block = IPyC
    eigenstrain_name = IPyC_IIDC_strain
    temperature = temperature
  []
[]

(examples/TRISO/parfume/parfume.i)

Added Moles of Gas

If using Sifgrs:

[Postprocessors]
  [fis_gas_released]
    type = ElementIntegralFisGasReleasedSifgrs
    block = fuel
  []
[]

(assessment/TRISO/benchmark/IAEA_CRP-6/fuel_performance/case_11/case_11_1D.i)

For CO production with UO $_2$ fuel:

[Postprocessors]
  [co_production]
    type = CarbonMonoxideProduction
    total_fissions = total_fissions
    initial_enrichment = 0.14029
    execute_on = 'initial timestep_end'
  []
[]

(test/tests/carbon_monoxide_production/carbon_monoxide_production_test.i)

TRISO single layer failure probability

In the Weibull theory, the failure probability is

P_f = 1 - \exp \left(-\int_V \left(\sigma/\sigma_0\right)^m dV \right)

The stress integration above is performed using the principle of independent action (PIA) model as follows:

\int_V \sigma^m dV = \int_V \left( \sigma_1^m + \sigma_2^m + \sigma_3^m \right)dV

[Postprocessors]
  [Weibull_failure_probability_IPyC]
    type = WeibullFailureProbability
    block = IPyC
    weibull_modulus = 6
    characteristic_strength = characteristic_strength_IPyC
  []
[]

(test/tests/triso_failure/triso_1d_weibull_probability.i)

[Materials]
  [characteristic_strength_PyC]
    type = PyCCharacteristicStrength
    temperature = temp
    X = 1.02
    flux_conversion_factor = 0.85
    block = 'IPyC OPyC'
  []
[]

(test/tests/triso_failure/triso_ipyc_characteristic_strength.i)

Effective mean strength for the layer

The effective mean strength for the layer is defined to be:

\sigma_{ms} = \frac{\sigma_0}{I_n^{\frac{1}{m}}}, I_n = \frac{\int_V \left( \sigma_1^m + \sigma_2^m + \sigma_3^m \right)dV}{\sigma_c^m}

$~$

[Postprocessors]
  [strength_SiC]
    type = WeibullEffectiveMeanStrength
    block = SiC
    weibull_modulus = 6
  []
[]

(test/tests/triso_failure/triso_1d_asphericity_failure.i)

[Materials]
  [characteristic_strength]
    type = GenericConstantMaterial
    prop_values = '9640000'
    prop_names = 'characteristic_strength'
    block = SiC
  []
[]

(test/tests/triso_failure/triso_1d_asphericity_failure.i)

TRISO failure determination

Whether or not particle failure occurs is determined by comparing the maximum stress and a strength that is sampled from a Weibull distribution

\sigma_0 > \sigma_{ms} ?

[Postprocessors]
  [failure_indicator_SiC]
    type = WeibullFailureOutputUsingCorrelation
    block = SiC
    weibull_modulus = 6
    stress_name = max_principal_stress
    effective_mean_strength = strength_SiC
  []
[]

(test/tests/triso_failure/triso_ipyc_characteristic_strength.i)

[Postprocessors]
  [triso_failure]
    type = TRISOFailureEvaluation
    SiC_failure = failure_indicator_SiC
  []
[]

(test/tests/triso_failure/triso_1d_asphericity_failure.i)

High-fidelity analysis and stress correlation

High-Fidelity Analysis

Perform high-fidelity analysis on cracked particles using parameters at their mean values.

$~$

$\rightarrow$

Stress Correlation

Perform similar one-dim analysis with same parameters.
Produce correlation functions to be used for calculating maximum stress in SiC layer.

$~$

\sigma_c = \frac{\bar{\sigma}_{highfidelity}}{\bar{\sigma}_{onedim}}\sigma_{onedim}

$~$

I_n = \frac{\int{\left(\sigma_1^m + \sigma_2^m + \sigma_3^m)\right) dV}}{\sigma_c^m}

$~$

$\rightarrow$

Failure Probability Determination

Perform 1D simulation to calculate stresses in Monte Carlo sampling
Use Weibull statistics to determine failure probability for particles.

IPyC cracking

$~$

[Mesh]
  coord_type = RZ
  [mesh]
    type = TRISO2DMeshGenerator
    elem_type = quad4
    coordinates = '0 ${coordinates1} ${coordinates2} ${coordinates2} ${coordinates3} ${coordinates4} ${coordinates5}'
    mesh_density = '20 8 0 4 4 4'
    block_names = 'fuel buffer IPyC SiC OPyC'
    num_sectors = 60
    aspect_ratio = ${aspect_ratio}
    all_bottom_left = True
  []
[]
[BCs]
  [no_disp_x]
    type = DirichletBC
    variable = disp_x
    boundary = xzero
    value = 0.0
  []

  [no_disp_y]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0.0
  []
[]

(examples/TRISO/correlation_function/h_asphericity/triso_asphericity.i)

$~$

IPyC cracking

Stress correlation function: $\sigma_c = - \frac{\bar{\sigma}_{2D}}{\bar{\sigma}_{1D}}\sigma_{1D}$

[Materials]
  [characteristic_strength_SiC]
    type = GenericConstantMaterial
    prop_values = '9640000'
    block = SiC
    prop_names = 'characteristic_strength'
  []

  [characteristic_strength_PyC]
    type = PyCCharacteristicStrength
    temperature = temperature
    X = 1.02
    block = 'IPyC OPyC'
  []
[]

[Postprocessors]
  [SiC_stress]
    type = ElementalVariableValue
    elementid = 6300
    variable = tangential_stress
  []

  [strength_SiC]
    type = WeibullEffectiveMeanStrength
    block = SiC
    weibull_modulus = 6
  []
[]

(examples/TRISO/correlation_function/h_asphericity/triso_asphericity.i)

Asphericity

Stress correlation function: $\sigma_c = \frac{\bar{\sigma}_{2D}}{\bar{\sigma}_{1D-min}} + \frac{\Delta \bar{\sigma}_{2D}}{\Delta \bar{\sigma}_{1D}} \Delta \sigma_{1D}$

[Mesh]
  coord_type = RZ
  [mesh]
    type = TRISO2DMeshGenerator
    elem_type = quad4
    coordinates = '0 ${coordinates1} ${coordinates2} ${coordinates2} ${coordinates3} ${coordinates4} ${coordinates5}'
    mesh_density = '20 8 0 4 4 4'
    block_names = 'fuel buffer IPyC SiC OPyC'
    num_sectors = 60
    aspect_ratio = ${aspect_ratio}
    all_bottom_left = True
  []
[]

(examples/TRISO/correlation_function/h_asphericity/triso_asphericity.i)

$~$

[Postprocessors]
  [SiC_stress]
    type = ElementalVariableValue
    elementid = 6300
    variable = tangential_stress
  []

  [strength_SiC]
    type = WeibullEffectiveMeanStrength
    block = SiC
    weibull_modulus = 6
  []
[]

(examples/TRISO/correlation_function/h_asphericity/triso_asphericity.i)

BISON TRISO WorkshopImplicit, parallel, fully-coupled nuclear fuel performance analysis Computational Mechanics and Materials Department Idaho National Laboratory