High Temperature Engineering Test Reactor (HTTR) Steady State Multiphysics Model Description

Mesh

The 'Mesh' block loads a homogenized mesh from an Exodus file.

[Mesh]
  [fmg]
    type = FileMeshGenerator
    file = '../mesh/full_core/thermal_mesh_in.e'
  []
[]

Figure 2: Radial cut of the 3D homogenized coarse mesh loaded in full_core_ht_steady.i. The reactor pressure vessel exchanges heat with the core via radiation and conduction.

Boundary Conditions (BCs)

The 'BCs' block imposes a boundary conditions that are all applied to the temperature variable.

'RPV_in_BC' imposes a forced convective boundary condition on the inner side of the RPV.

[BCs]
  [RPV_in_BC]
    # FIXME: this BC adds energy to the RPV that is not removed from the fluid! (but only ~40kW for 9MW so about 0.6K error on the inlet fluid temperature)
    type = ConvectiveFluxFunction # (Robin BC)
    variable = Tsolid
    boundary = '200' # inner RPV
    coefficient = 30.7 # calculated value for 9MW forced convection # old value: 1e2 W/K/m^2 - arbitrary value
    T_infinity = T_init
  []
[]

'RPV_out_BC' imposes a natural convective boundary condition on the outer side of the RPV.

[BCs]
  [RPV_out_BC]
    # k \nabla T = h (T- T_inf) at RPV outer boundary
    type = ConvectiveFluxFunction # (Robin BC)
    variable = Tsolid
    boundary = '100' # outer RPV
    coefficient = 2.58 # calculated value for natural convection following NRC report method
    T_infinity = T_outside # matches Figure 13 value concrete shield
  []
[]

'radiative_BC' imposes a radiative boundary condition on the outer side of the RPV.

[BCs]
  [radiative_BC]
    # radiation from RPV to bio shield
    type = InfiniteCylinderRadiativeBC
    boundary = '100' # outer RPV
    variable = Tsolid
    boundary_radius = ${rpv_outer_radius}
    boundary_emissivity = ${rpv_emissivity}
    cylinder_radius = ${vcs_radius}
    cylinder_emissivity = ${vcs_emissivity} # matching NRC value, 0.79 for inner RPV
    Tinfinity = T_outside # matches Figure 13 value concrete shield
  []
[]

Kernels

The 'heat_conduction' block defines the heat conduction in the reactor pressure vessel.

[Kernels]
  [heat_conduction]
    type = HeatConduction
    variable = Tsolid
    block = '${rpv_blocks}'
  []
[]

The 'aniso_heat_conduction' block defines the anisotropic heat conduction inside the core. The anisotropy accounts for holes in the radial direction, within fuel and CR blocks.

[Kernels]
  [aniso_heat_conduction]
    type = AnisoHeatConduction
    variable = Tsolid
    thermal_conductivity = aniso_thermal_conductivity
    block = '${full_fuel_blocks} ${full_cr_blocks} ${full_pr_blocks}' # all blocks but RPV
  []
[]

The time derivative of the energy equation is added using 'HeatConductionTimeDerivative'.

[Kernels]
  [heat_ie]
    type = HeatConductionTimeDerivative
    variable = Tsolid
    block = '${rpv_blocks}'
  []
[]

[Kernels]
  [aniso_heat_ie]
    type = HeatConductionTimeDerivative
    variable = Tsolid
    block = '${full_fuel_blocks} ${full_cr_blocks} ${full_pr_blocks}' # all blocks but RPV
    specific_heat = aniso_specific_heat
  []
[]

A term that accounts for radiative and conductive heat transfers through the gap between the fuel pins and the graphite blocks is added.

[Kernels]
  [heat_loss_conductance_active_fuel]
    type = Removal
    variable = Tsolid
    block = '${fuel_blocks_31pin} ${fuel_blocks_33pin}'
    sigr = gap_conductance_mat
  []
[]

[Kernels]
  [heat_gain_conductance_active_fuel]
    type = MatCoupledForce
    variable = Tsolid
    v = inner_Twall
    block = '${fuel_blocks_31pin} ${fuel_blocks_33pin}'
    material_properties = gap_conductance_mat
  []
[]

A term that accounts for the heat extracted by fluid convection in the fuel and CR cooling channels is represented with those two kernels:

[Kernels]
  [heat_loss_convection_outer]
    type = Removal
    variable = Tsolid
    block = '${full_fuel_blocks} ${full_cr_blocks}'
    sigr = Hw_outer_homo_mat
  []
[]

[Kernels]
  [heat_gain_convection_outer]
    type = MatCoupledForce
    variable = Tsolid
    v = Tfluid
    block = '${full_fuel_blocks} ${full_cr_blocks}'
    material_properties = Hw_outer_homo_mat
  []
[]

Thermal Contact

The 'ThermalContact' block is used to model radiation and conduction between the permanent reflectors and the reactor pressure vessel. The gap is filled with cooled helium gas. The two boundaries involved are '200' (boundary id) on the reactor pressure vessel and '1' (boundary id) on the permanent reflectors.

[ThermalContact]
  [RPV_gap]
    type = GapHeatTransfer
    gap_geometry_type = 'CYLINDER'
    emissivity_primary = 0.8
    emissivity_secondary = 0.8 # varies from 0.6 to 1.0 in RPV paper, 0.79 in NRC paper
    variable = Tsolid # graphite -> rpv gap
    primary = '200'
    secondary = '1'
    gap_conductivity = 0.2091 # for He at 453K, 2.8 MPa
    quadrature = true
    #paired_boundary = '1'
    cylinder_axis_point_1 = '0 0 0'
    cylinder_axis_point_2 = '5.22 0 0'
  []
[]

Materials

The material solid thermal properties are defined using 'HeatConductionMaterial' and 'AnisoHeatConductionMaterial'; this includes thermal conductivity and specific heat.

[Materials]
  [graphite_moderator]
    type = AnisoHeatConductionMaterial
    block = '${full_fuel_blocks} ${full_cr_blocks} ${full_pr_blocks}' # all blocks but RPV
    temperature = Tsolid
    base_name = 'aniso'
    reference_temperature = 0.0
    thermal_conductivity_temperature_coefficient_function = IG110_k
    thermal_conductivity = '1    0    0
                            0 0.54    0
                            0    0 0.54' # isotropic for testing only, otherwise one in the x-direction
    specific_heat = IG110_cp
  []
[]

[Materials]
  [ss304_rpv]
    type = HeatConductionMaterial
    block = '${rpv_blocks}'
    thermal_conductivity_temperature_function = ssteel_304_k
    specific_heat = 500 # J/kg/K
    temp = Tsolid
  []
[]

The thermal conductivity is defined as a function of temperature defined in the 'Functions' block. It is divided by the temperature to get the desired behavior with AnisoHeatConductionMaterial.

[Functions]
  [IG110_k]
    # thermal conductivity divided by the temperature (dirty trick to make AnisoHeatConductionMaterial define a temperature-independent anistropic thermal conductivity)
    type = ParsedFunction
    symbol_values = '6.632e+01 -4.994e-02 1.712e-05' # [W/m/K]
    symbol_names = 'a0        a1         a2       '
    expression = '((a0 + a1 * t + a2 * t * t) - 1) / t'
  []
[]

Additionally the densities of the various regions are defined in the 'Materials' block.

[Materials]
  [graphite_density_fuel33]
    type = Density
    block = '${fuel_blocks_33pin}'
    density = '${fparse 1770 * vol_frac_fuel33}'
  []
[]

[Materials]
  [graphite_density_fuel31]
    type = Density
    block = '${fuel_blocks_31pin}'
    density = '${fparse 1770 * vol_frac_fuel31}'
  []
[]

[Materials]
  [graphite_density_rr_fuel33]
    type = Density
    block = '${rr_fuel_blocks_33pin}'
    density = '${fparse 1760 * vol_frac_fuel33}'
  []
[]

[Materials]
  [graphite_density_rr_fuel31]
    type = Density
    block = '${rr_fuel_blocks_31pin}'
    density = '${fparse 1760 * vol_frac_fuel31}'
  []
[]

[Materials]
  [graphite_density_cr]
    type = Density
    block = '${full_cr_blocks}'
    density = '${fparse 1770 * vol_frac_cr}'
  []
[]

[Materials]
  [graphite_density_rest]
    type = Density
    block = '${full_pr_blocks}'
    density = '${fparse 1760}' # In reality, the actual PR blocks are made of PGX graphite of density 1.732 g/cc but here we assume IG110 graphite for now.
  []
[]

[Materials]
  [ss304_density]
    type = Density
    block = '${rpv_blocks}'
    density = 800 #1/10th of the typical value (8000) for ss304 in kg/m^3 to accelerate steady-state solution
  []
[]

Mesh

The 'Mesh' block loads a mesh to simulate the 5 fuel pins within the fuel column.

[Mesh]
  [fmg]
    type = FileMeshGenerator
    file = '../mesh/fuel_element/HTTR_fuel_pin_2D_refined_m_5pins_axialref.e'
  []
  uniform_refine = 0 # do not modify if using ThermalContact
  coord_type = RZ
  rz_coord_axis = X
[]

Figure 3: The 2D heterogeneous mesh in RZ geometry loaded in fuel_elem_steady.i. Shown is one of the five fuel pins from the mesh; composed of a fuel block, in grey, a graphite sleeve, in red, and a graphite moderator block, in green. Areas 1, 2, and 3 contain the thermal contact boundaries and will be described in greater detail in the next section.

Thermal Contact

The ThermalContact block solves for the heat conduction and radiation in the gap between the graphite sleeve and the top of the fuel compact, the gap between the side of graphite sleeve and side of the fuel compact, and the gap between the side of the graphite sleeve and the side of the moderator block.

[ThermalContact]
  [gap_top]
    type = GapHeatTransfer
    gap_geometry_type = 'PLATE'
    emissivity_primary = 0.8 # MHTGR: 0.85, Shi PhD: 0.8
    emissivity_secondary = 0.8 # MHTGR: 0.85, Shi PhD: 0.8
    # coarser side should be primary
    # (in theory using quadrature = true should make it independent but other errors appear and not exactly conservative)
    primary = Sleeve_top
    secondary = Fuel_top
    gap_conductivity = 0.2735 # He at 668K, 2.8 MPa - Original: 1e-5 W/m/K, Shi PhD: 0.35 (but no gap on the inside)
    variable = temp
    quadrature = true
  []
[]

Figure 4: A close up of the fuel pin mesh, Area 1, showing the gap between the 'Sleeve_top' (the section of the graphite sleeve to the right, in red) and the 'Fuel_top' (the section of the fuel pin to the left, in grey), where heat transfer takes place; the boundaries are shown as faint pink lines.

[ThermalContact]
  [gap_side]
    type = GapHeatTransfer
    gap_geometry_type = 'CYLINDER'
    emissivity_primary = 0.8 # MHTGR: 0.85, Shi PhD: 0.8
    emissivity_secondary = 0.8 # MHTGR: 0.85, Shi PhD: 0.8
    primary = Sleeve_side
    secondary = Fuel_side
    gap_conductivity = 0.2735 # He at 668K, 2.8 MPa - Original: 1e-5 W/m/K, Shi PhD: 0.35 (but no gap on the inside)
    variable = temp
    quadrature = true
  []
[]

Figure 5: A close up of the fuel pin mesh, Area 2, showing the 'Sleeve_side' (the section of the graphite sleeve at the top, in red) and the 'Fuel_side' (the section of the fuel pin at the bottom, in grey), which contains a very slight gap in which heat transfer takes place. The boundaries are shown as faint blue and yellow lines.

[ThermalContact]
  [cooling_channel]
    type = GapHeatTransfer
    gap_geometry_type = 'CYLINDER'
    emissivity_primary = 0.8 # MHTGR: 0.85, Shi PhD: 0.8
    emissivity_secondary = 0.8 # MHTGR: 0.85, Shi PhD: 0.8
    primary = Inner_outer_ring
    secondary = Relap7_Tw
    gap_conductivity = 0.2309 # He at 523K, 2.8MPa
    variable = temp
    quadrature = true
  []
[]

Figure 6: A close up of the fuel pin mesh, Area 3, showing the gap between the 'Inner_outer_ring' (the outer wall of the graphite sleeve, in red) and the 'Relap7_Tw' (the inner wall of the graphite block, in green) where heat transfer takes place. The boundaries are shown as faint teal and black lines.

Boundary Conditions (BCs)

We impose a zero heat flux using a Neumann boundary condition on the bottom of the graphite sleeve.

[BCs]
  [Neumann]
    type = NeumannBC
    variable = temp
    value = 0.0
    boundary = 'Neumann Sleeve_bot'
  []
[]

A pseudo-Dirichlet boundary is also defined; it weakly imposes 'T = Tmod' on the outer most radial wall of the graphite block.

[BCs]
  [outside_bc]
    type = CoupledConvectiveHeatFluxBC
    variable = temp
    boundary = 'Outer_outer_ring'
    T_infinity = Tmod
    htc = 1e5
  []
[]

The convective heat transfer with the helium channel is modeled using 'CoupledConvectiveHeatFluxBC'.

The 'convective_inner' block specifies a heat transfer coefficient for the inner radius of the cooling channel in W/K/.

[BCs]
  [convective_inner]
    type = CoupledConvectiveHeatFluxBC
    boundary = 'Relap7_Tw'
    variable = temp
    T_infinity = Tfluid
    htc = Hw_inner
  []
[]

The 'convective_outer' block specifies a heat transfer coefficient for the outer radius of the cooling channel in W/K/.

[BCs]
  [convective_outer]
    type = CoupledConvectiveHeatFluxBC
    boundary = 'Inner_outer_ring'
    variable = temp
    T_infinity = Tfluid
    htc = Hw_outer
  []
[]

Kernels

The 'Kernels' block defines a steady-state heat conduction equation. A source term is defined using the 'CoupledForce' kernel to model how much power is deposited in the fuel compact.

[Kernels]
  [heat]
    type = HeatConduction
    variable = temp
  []
  # [heat_ie] # omit to speed up steady state convergence
  #   type = HeatConductionTimeDerivative
  #   variable = temp
  # []
  [heat_source]
    type = CoupledForce
    variable = temp
    block = '${fuel_blocks}'
    v = scaled_power_density
  []
[]