Microscale simulation of the TRISO fuel matrix

Geometry

Similar to the pebbles, we represent the extended TRISO particle with a 1D spherical mesh. The CartesianMeshGenerator allows us to align mesh boundaries with the geometrical boundaries of each material in the TRISO. The graphite matrix around the TRISO particles is represented by an additional layer of graphite, sized to preserve the total packing fraction of graphite.

[Mesh]
  coord_type = 'RSPHERICAL'
  type = MeshGeneratorMesh
  block_id = '1 2 3 4 5 6'
  block_name = 'uo2 buffer ipyc sic opyc graphite_matrix'

  [gen_mesh]
    type = CartesianMeshGenerator
    dim = 1

    dx = '${rUO2} ${fparse rbuffer - rUO2} ${fparse ripyc - rbuffer}
          ${fparse rSiC - ripyc} ${fparse ropyc - rSiC}
          ${fparse stainsby_sphere_radius - ropyc}'
    ix = '20 10 4 4 4 20'
    subdomain_id = '1 2 3 4 5 6'
  []
[]

Defining the equation: variables and kernels

The temperature variable in this simulation is the unshifted temperature, or the temperature solution of the microscale equation, not the physical temperature. We use the auxiliary system, an AuxKernel and an AuxVariable to compute the physical temperature, which would be used to compute the material properties.

[Variables]
  [T_unshifted]
  []
[]

[AuxVariables]
  [fuel_matrix_heat_source]
  []
  [T]
  []
[]

[AuxKernels]
  # Shift variable to average 0
  [T]
    type = ShiftedAux
    variable = T
    unshifted_variable = T_unshifted
    shift_postprocessor = negative_average_T_unshifted
    execute_on = 'TIMESTEP_END'
  []
[]

The microscale equation is decomposed in four kernels:

• the time derivative kernel. It's inactive by default as we neglect the time dependence of the temperature at the microscale. The time scale for the evolution of the microscale is much smaller. We instead solve a steady state problem.

[Kernels]
  [time]
    type = ADHeatConductionTimeDerivative
    variable = T_unshifted
    specific_heat = 'cp_s'
    density_name = 'rho_s'
    block = '1 2 3 4 5 6'
  []
[]

• the heat conduction kernel. We only need to specify the name of the thermal conductivity material property and the name of the temperature variable.

[Kernels]
  [diffusion]
    type = ADHeatConduction
    variable = T_unshifted
    thermal_conductivity = 'k_s'
    block = '1 2 3 4 5 6'
  []
[]

• the heat source. This term is only present in the center of the TRISO particle, in the fuel. It is assumed uniform, as we neglect spatial self-shielding over the particle for the heat deposition. We rescale it based on the actual UO2 volume, as the heat source was averaged over the solid phase / fuel matrix at various stages of the simulation.

[Kernels]
  [heat_source]
    type = HeatSrc
    variable = T_unshifted
    heat_source = fuel_matrix_heat_source
    scaling_factor = '${fparse 1.0/fuel_matrix_phase_fraction/TRISO_phase_fraction/UO2_phase_fraction}'
    block = '1'
  []
[]

• the heat sink. This term makes the total heat source have a zero average, so that it represents the fluctuation term in our multiscale approach. It is distributed over the entire microscale fuel matrix.

[Kernels]
  [heat_sink_for_zero_average_over_matrix]
    type = HeatSrc
    variable = T_unshifted
    heat_source = fuel_matrix_heat_source
    scaling_factor = '${fparse - 1.0/fuel_matrix_phase_fraction}'
    block = '1 2 3 4 5 6'
  []
[]

We add a boundary condition to this heat diffusion problem. This boundary condition, null temperature on the right for the unshifted temperature, is arbitrary, but it is required for the diffusion problem. The default zero gradient boundary conditions on both sides would allow an infinity of solutions, all offset by a constant. This boundary condition is superseded by the normalization condition imposed later on the shifted temperature.

[BCs]
  [right]
    type = DirichletBC
    variable = T_unshifted
    value = 0.0
    boundary = 'right'
  []
[]

Outputs

We measure using Postprocessors the maximum and average temperatures of each phase, which will be transferred for visualization and analysis purposes to a core-scale application.

[Postprocessors]
  [average_T_unshifted]
    type = ElementAverageValue
    variable = T_unshifted
    execute_on = 'LINEAR'
  []
  [negative_average_T_unshifted]
    type = ScalePostprocessor
    value = average_T_unshifted
    scaling_factor = -1.0
    execute_on = 'TIMESTEP_END'
  []
  [surface_T]
    type = PointValue
    variable = T
    point = '${stainsby_sphere_radius} 0.0 0.0'
    execute_on = 'TIMESTEP_END'
  []

  # for visualization only
  [max_T_UO2]
    type = ElementExtremeValue
    variable = T
    value_type = max
    block = '1'
    execute_on = 'TIMESTEP_END'
  []
  [average_T_UO2]
    type = ElementAverageValue
    variable = T
    block = '1'
    execute_on = 'TIMESTEP_END'
  []
  [average_T_matrix]
    type = ElementAverageValue
    variable = T
    block = '6'
    execute_on = 'TIMESTEP_END'
  []
[]