Thermal Fluid and Thermal Conduction Model

Thermal Fluids Model

The thermal-hydraulics model for Pronghorn is based on ongoing NEAMS work (Schunert et al., 2022) (Schunert et al., 2021) with some additional improvements. The Pronghorn model uses the weakly compressible finite volume formulation for discretizing the fluid mass, fluid momentum, fluid energy, and solid energy conservation equations. The model includes a riser and bypass flow channels, and the fueling chute as an open flow region.

The cold fluid from the circulators enters the core via the vertical risers in the reflector region. The flow then enters the cold plenum, where the flow is diverted into the cavity, upper reflector, and control and shutdown system bypass channels. From the upper cavity, the fluid enters the active core region, then the lower reflector, and finally the outlet plenum.

The radiative heat transfer at the outer boundary of the reactor vessel has a small impact on steady-state calculations and a larger impact during the loss of forced cooling (DLOFC) transients. Also, during the loss of flow transient, fluid inflow and outflow boundary conditions are not changed to wall boundary conditions, which leads to a significant change in the helium leaving the core due to thermal expansion, which tends to increase temperature estimates.

The model for the thermal fluid calculations can be found in this input. A depiction of the geometry and materials in the Pronghorn thermal-hydraulic model and the fluid flow path is shown in Figure 1.

Figure 1: Depiction of the geometry and masterials in the Pronghorn thermal-hydraulic model including the fluid flow path.

The boundary condition for the Pronghorn model of the HTR-PM include: * Inlet helium flow rate of 96.0 kg/s. * The inlet fluid temperature is set to 523.15 K. * The outlet fluid pressure is set to Pa. * Inlet velocity, outlet pressure, slip-wall, and symmetry boundary conditions are used for the fluid mass, momentum, and energy equations. * All walls are assumed to be adiabatic for the fluid energy equation, and conjugate heat transfer is treated as a volumetric phenomenon. * The solid energy equations have adiabatic boundary conditions except for the outside of the pressure vessel * Radiative and convective boundary conditions were applied between the pressure vessel and isothermal cylindrical reactor cavity cooling system (RCCS) panel with an inner diameter of 4 m, a temperature of K, a heat transfer coefficient of 5 W/m.K, and a surface emissivities are assumed to be 0.8.

The model input starts by defining geometric parameters for the problem. Then, global parameters are defined in [GlobalParams] and the mesh is defined in the [Mesh] block

[GlobalParams]
  acceleration = '0.0 -9.81 0.0' # Gravity acceleration (m/s2).
  fp = fluid_properties_obj
  porosity = 'porosity'
  pebble_diameter = ${pebble_diameter}
  T_solid = T_solid
  rhie_chow_user_object = pins_rhie_chow_interpolator
[]

[Mesh]
  type = MeshGeneratorMesh
  block_id = '1 2 3 4 5 6 7 8 10 12 61 71 9 11'
  block_name = 'pebble_bed
                top_reflector
                bottom_reflector
                top_cavity
                hot_plenum
                cold_plenum
                side_reflector
                carbon_brick
                core_barrel
                rpv
                riser
                bypass
                refl_barrel_gap
                barrel_rpv_gap'
  uniform_refine = 1

  [cartesian_mesh]
    type = CartesianMeshGenerator
    dim = 2

    dx = ' 0.250 0.250 0.250 0.250 0.250 0.250
           0.010 0.050
           0.130
           0.080 0.080 0.080 0.200 0.120 0.010 0.240
           0.150 0.040 0.160 0.150 '
    ix = ' 1 1 1 1 1 1
           1 1
           1
           1 1 1 2 1 1 1
           1 1 1 1 '

    dy = ' 0.400 0.400 0.100 0.100
           0.800 0.300 0.200 0.300 0.216 0.412
           0.550 0.550 0.550 0.550 0.550 0.550 0.550 0.550 0.550 0.550
           0.550 0.550 0.550 0.550 0.550 0.550 0.550 0.550 0.550 0.550
           0.760 0.712 0.300
           0.400 0.400 '
    iy = ' 1 1 1 1
           2 1 1 1 1 1
           1 1 1 1 1 1 1 1 1 1
           1 1 1 1 1 1 1 1 1 1
           2 2 1
           1 1 '

    subdomain_id = ' 8  8 8 8 8 8  8 8 8  8 8 8 8  8 8 8  9  10  11  12
                     8  8 8 8 8 8  8 8 8  8 8 8 8  8 8 8  9  10  11  12
                     7  7 7 7 7 7  7 7 7  7 7 7 7  7 8 8  9  10  11  12
                     7  7 7 7 7 7  7 7 7  7 7 7 7  7 8 8  9  10  11  12
                     5  5 5 5 5 5  5 5 5  7 7 7 61 7 8 8  9  10  11  12
                     3  3 3 3 3 3  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     3  3 3 3 3 3  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     3  3 3 3 3 3  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     3  3 3 3 3 3  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     3  3 3 3 3 3  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     1  1 1 1 1 1  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     4  4 4 4 4 4  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     4  2 2 2 2 2  7 7 71 7 7 7 61 7 8 8  9  10  11  12
                     6  6 6 6 6 6  6 6 6  6 6 6 6  7 8 8  9  10  11  12
                     7  7 7 7 7 7  7 7 7  7 7 7 7  7 8 8  9  10  11  12
                     8  8 8 8 8 8  8 8 8  8 8 8 8  8 8 8  9  10  11  12 '
  []

  # Side sets for gap conductance model.
  [reflector_barrel_gap_inner]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = 8
    paired_block = 9
    input = cartesian_mesh
    new_boundary = reflector_barrel_gap_inner
  []
  [reflector_barrel_gap_outer]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = 10
    paired_block = 9
    input = reflector_barrel_gap_inner
    new_boundary = reflector_barrel_gap_outer
  []
  [barrel_rpv_gap_inner]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = 10
    paired_block = 11
    input = reflector_barrel_gap_outer
    new_boundary = barrel_rpv_gap_inner
  []
  [barrel_rpv_gap_outer]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = 12
    paired_block = 11
    input = barrel_rpv_gap_inner
    new_boundary = barrel_rpv_gap_outer
  []

  # Side sets for inflow and outflow conditions.
  [reactor_inlet]
    type = ParsedGenerateSideset
    included_subdomains = '61'
    combinatorial_geometry = 'abs(y-1) < 1e-3'
    fixed_normal = true
    normal = '0 -1 0'
    input = barrel_rpv_gap_outer
    new_sideset_name = reactor_inlet
  []
  [reactor_outlet]
    type = SideSetsAroundSubdomainGenerator
    block = '5'
    fixed_normal = true
    normal = '1 0 0'
    input = reactor_inlet
    new_boundary = reactor_outlet
  []
  [riser_walls]
    type = ParsedGenerateSideset
    included_subdomains = '61'
    included_neighbors = '7'
    combinatorial_geometry = 'y > 1 + 1e-3'
    input = reactor_outlet
    new_sideset_name = riser_walls
  []

  # Side sets for wall boundaries.
  [cold_plenum_walls]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '6'
    paired_block = '7'
    input = riser_walls
    new_boundary = cold_plenum_walls
  []
  [hot_plenum_walls]
    type = ParsedGenerateSideset
    included_subdomains = 5
    included_neighbors = 7
    combinatorial_geometry = 'x < 1.69'
    input = cold_plenum_walls
    new_sideset_name = hot_plenum_walls
  []
  [bypass_walls]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '71'
    paired_block = '7'
    input = hot_plenum_walls
    new_boundary = 'bypass_wall'
  []
  [pbed_inner]
    type = ParsedGenerateSideset
    included_subdomains = '1 2 3 4 5 6'
    combinatorial_geometry = '( abs(x - 0.000) < ${geometric_tolerance} &
                              y > ${fparse 1.000 - geometric_tolerance} &
                              y < ${fparse 16.00 + geometric_tolerance} )'
    input = bypass_walls
    new_sideset_name = pbed_inner
  []
  [pbed_outer]
    type = ParsedGenerateSideset
    included_subdomains = '1 2 3 4'
    combinatorial_geometry = '( abs(x - ${pbed_r}) < ${geometric_tolerance} &
                              y > ${fparse 1.800 - geometric_tolerance} &
                              y < ${fparse 15.70 + geometric_tolerance} ) '
    input = pbed_inner
    new_sideset_name = pbed_outer
  []
  [bypass_hot_plenum_interface]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = '71'
    paired_block = '5'
    new_boundary = 'bypass_hot_plenum_interface'
    input = pbed_outer
  []

  coord_type = RZ
[]

The definition of the pronghorn Navier-Stokes equation for the solution of the weakly-compressible fluid, and the properties of the fluid go into Physics block.

[Physics]
  [NavierStokes]
    [Flow]
      [all]
        # basic settings
        block = ${fluid_blocks}
        compressibility = 'weakly-compressible'
        gravity = '0.0 -9.81 0.0'

        # Porous treatement
        porous_medium_treatment = true
        friction_types = 'darcy forchheimer'
        friction_coeffs = 'Darcy_coefficient Forchheimer_coefficient'
        consistent_scaling = ${scaling}
        porosity_smoothing_layers = 0
        use_friction_correction = true

        # fluid properties
        density = 'rho'
        dynamic_viscosity = 'mu'

        # initial conditions
        initial_velocity = '1e-6 1e-6 0'
        initial_pressure = '${p_outlet}'

        # boundary conditions
        inlet_boundaries = 'reactor_inlet'
        momentum_inlet_types = 'flux-mass'
        flux_inlet_pps = 'set_inlet_mfr'
        flux_inlet_directions = '0 1 0'

        outlet_boundaries = 'reactor_outlet'
        momentum_outlet_types = 'fixed-pressure'
        pressure_functors = '${p_outlet}'

        wall_boundaries = 'pbed_inner pbed_outer hot_plenum_walls cold_plenum_walls riser_walls bypass_wall'
        momentum_wall_types = 'symmetry   slip       slip             slip              slip        slip'

        # numerical scheme
        pressure_face_interpolation = average
        momentum_advection_interpolation = upwind
        mass_advection_interpolation = upwind
        # prevents solution jump on future restarts
        time_derivative_contributes_to_RC_coefficients = false
      []
    []
    [FluidHeatTransfer]
      [all]
        block = ${fluid_blocks}

        # numerical scheme
        energy_advection_interpolation = upwind
        system_names = 'nl0'

        # convective heat transfer
        ambient_convection_blocks = '1 2 3 5 6 61 71'
        ambient_convection_alpha = 'alpha'
        ambient_temperature = 'T_solid'

        # fluid properties
        thermal_conductivity = 'kappa'
        specific_heat = 'cp'

        # initial conditions
        initial_temperature = '${T_inlet}'

        # boundary conditions
        # see Flow physics for list of boundaries
        energy_inlet_types = 'flux-mass'
        energy_inlet_functors = '${T_inlet}'

        energy_wall_types = 'heatflux   heatflux   heatflux         heatflux          heatflux    heatflux'
        energy_wall_functors = '0          0          0                0                 0           0'
      []
    []
  []
[]

Material properties are defined in the [Materials] block

[Materials]
  ## natural htc for BC
  [natural_htc_mat]
    type = ADGenericConstantMaterial
    prop_names = 'natural_htc'
    prop_values = '5'
  []
[]

The characteristics of the execution of the problem are defined in [Executioner] block

[Executioner]
  type = Transient

  # solver parameters
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -ksp_gmres_restart -pc_factor_shift_type -mat_mumps_icntl_20'
  petsc_options_value = 'lu       mumps                         100                 NONZERO              0'
  automatic_scaling = true
  nl_abs_tol = 1e-5
  line_search = l2
  nl_max_its = 50

  # time stepping
  end_time = 1e6
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.5
    optimal_iterations = 9
    iteration_window = 2
    growth_factor = 2
    cutback_factor = 0.5
  []
[]

And control of output is defined in [Outputs] block

[Outputs]
  csv = true
  exodus = true
  checkpoint = true
  [console]
    type = Console
    hide = 'area_pp_reactor_inlet set_inlet_mfr'
  []
  print_linear_converged_reason = false
  print_linear_residuals = false
  print_nonlinear_converged_reason = false
[]

Pebble Conduction Model:

The pebble heat conduction model solves the 1D spherical conduction problem assuming a thermal equilibrium approximation in pebble simulations during transient calculations. Several sources of heat transfer nonuniformity around the pebble were ignored: the coolant flow orientation, pebble-to-pebble contact, pebble-to-reflector contact, and radiation. A Dirichlet boundary condition is set at the TRISO surface to obtain the fuel temperature, and a Neumann boundary condition would improve energy conservation.

The model of the pebble conduction is found in /htgr/htr-pm-2/pebble_triso.i. The model set up is visually depicted in Figure 2.

Figure 2: Depiction of the pebble conduction model.

In this model, the mesh is defined as follows

[Mesh]
  block_id = '1 2 3 4 5 6 7'
  block_name = 'core
                 shell
                 kernel
                 buffer
                 ipyc
                 sic
                 opyc'
  dim = 1
  [pebble_mesh]
    type = CartesianMeshGenerator
    dim = 1
    dx = '2.50e-02 5.00e-03'
    ix = '15 3'
    subdomain_id = '1 2'
  []
  [triso_mesh]
    type = CartesianMeshGenerator
    dim = 1
    dx = '2.50e-04 9.00e-05 4.00e-05 3.50e-05 4.00e-05'
    ix = '21 8 3 3 3'
    subdomain_id = '3 4 5 6 7'
  []
  [mesh_combine]
    type = CombinerGenerator
    inputs = 'pebble_mesh triso_mesh'
  []
  [pebble_surface]
    type = SideSetsAroundSubdomainGenerator
    block = '2'
    fixed_normal = 1
    normal = '1 0 0'
    input = mesh_combine
    new_boundary = pebble_surface
  []
  [triso_surface]
    type = SideSetsAroundSubdomainGenerator
    block = '7'
    fixed_normal = 1
    normal = '1 0 0'
    input = pebble_surface
    new_boundary = triso_surface
  []
  coord_type = 'RSPHERICAL'
[]

The model also defines pebble and TRISO temperatures as follows

[Variables]
  [T_pebble]
    block = '1 2'
    # initial_condition = ${initial_temperature}
  []
  [T_triso]
    block = '3 4 5 6 7'
    # initial_condition = ${initial_temperature}
  []
[]

Materials are defined in [Materials] block

[Materials]
  [pebble_core]
    type = PronghornSteadyStateSolidMaterial
    T_solid = T_pebble
    solid = pebble_core
    block = '1'
  []
  [pebble_shell]
    type = PronghornSteadyStateSolidMaterial
    T_solid = T_pebble
    solid = gmatrix
    block = '2'
  []

  # TRISO.
  [kernel]
    type = PronghornSteadyStateSolidMaterial
    T_solid = T_triso
    solid = kernel
    block = '3'
  []
  [buffer]
    type = PronghornSteadyStateSolidMaterial
    T_solid = T_triso
    solid = buffer
    block = '4'
  []
  [ipyc]
    type = PronghornSteadyStateSolidMaterial
    T_solid = T_triso
    solid = ipyc
    block = '5'
  []
  [sic]
    type = PronghornSteadyStateSolidMaterial
    T_solid = T_triso
    solid = sic
    block = '6'
  []
  [opyc]
    type = PronghornSteadyStateSolidMaterial
    T_solid = T_triso
    solid = opyc
    block = '7'
  []
[]

materials thermal conductivity and burnup are defined in [Functions] block:

[Functions]
  [uo2_k]
    type = ParsedFunction
    expression = 'if(bnp < 1e-10, (115.8/(7.5408+17.692*(t/1000)+3.6142*(t/1000)^2)+7410.5*(t/1000)^(-5./2.)*exp(-16.35/(t/1000))),
           (115.8/(7.5408+17.692*(t/1000)+3.6142*(t/1000)^2)+7410.5*(t/1000)^(-5./2.)*exp(-16.35/(t/1000)))*
           (1.09/bnp^(3.265)+0.0643*sqrt(t/bnp))*atan(1./(1.09/bnp^(3.265)+0.0643*sqrt(t/bnp)))*
           (1+0.019*bnp/(3.-0.019*bnp)*(1.+exp(-(t-1200)/100))^(-1))*
           (1.-0.2/(1+exp((t-900.)/80.))) )'
    symbol_names = 'bnp'
    symbol_values = 'fima'
  []
  [buffer_k]
    type = ParsedFunction
    expression = 244.3/2*t^(-0.574)*(970/(2.2*(1930.-970)+970))*(1.-0.336*(1.-exp(-1.005*gam))-3.50e-2*gam)
    symbol_names = 'gam'
    symbol_values = 'fluence'
  []
  [pyc_k]
    type = ParsedFunction
    expression = 244.3*t^(-0.574)*(1900/(2.2*(1930.-1900)+1900))*(1.-0.336*(1.-exp(-1.005*gam))-3.50e-2*gam)
    symbol_names = 'gam'
    symbol_values = 'fluence'
  []
  [sic_k]
    type = ParsedFunction
    expression = (17885/t+2.)*exp(-0.1277*gam)
    symbol_names = 'gam'
    symbol_values = 'fluence'
  []
  [gmatrix_k]
    type = ParsedFunction
    expression = 47.4*(1-9.7556E-4*(t-373.15)*exp(-6.036E-4*(t-273.15)))*(1740/(2.2*(1700.-1740)+1740))*(1.-0.336*(1.-exp(-1.005*gam))-3.50e-2*gam)
    symbol_names = 'gam'
    symbol_values = 'fluence'
  []
  [fluence]
    type = ParsedFunction
    expression = 7.41611E-06*bnp*bnp*bnp-5.36979E-06*bnp*bnp+1.37527E-02*bnp-4.48921E-02
    symbol_names = 'bnp'
    symbol_values = 'burnup'
  []
  [fima]
    type = ParsedFunction
    expression = -2.022642E-06*bnp*bnp+1.053601E-03*bnp
    symbol_names = 'bnp'
    symbol_values = 'burnup'
  []
[]

These properties are linked to a type of object in the [UserObject] block as follows

[UserObjects]
  [kernel]
    type = FunctionSolidProperties
    k_s = uo2_k
  []
  [buffer]
    type = FunctionSolidProperties
    k_s = buffer_k
  []
  [ipyc]
    type = FunctionSolidProperties
    k_s = pyc_k
  []
  [sic]
    type = FunctionSolidProperties
    k_s = sic_k
  []
  [opyc]
    type = FunctionSolidProperties
    k_s = pyc_k
  []
  [gmatrix]
    type = FunctionSolidProperties
    k_s = gmatrix_k
  []

  # Mixtures.
  [triso]
    type = CompositeSolidProperties
    materials = 'kernel buffer ipyc sic opyc'
    fractions = '0.1659  0.2514  0.1653  0.1762  0.2412' # volume fractions.
    k_mixing = 'series'
  []
  [pebble_core]
    type = CompositeSolidProperties
    materials = 'triso gmatrix'
    fractions = '0.090484107  0.909515893' # volume fractions.
    k_mixing = 'chiew'
  []
[]

The boundary conditions are applied in the [BCs] block

[BCs]
  [pebble_surface_temp]
    type = PostprocessorDirichletBC
    variable = T_pebble
    postprocessor = T_surface
    boundary = 'pebble_surface'
  []
  [triso_surface_temp]
    type = PostprocessorDirichletBC
    variable = T_triso
    postprocessor = pebble_core_average_temp
    boundary = 'triso_surface'
  []
[]

Finally, the characteristics of the execution process are provided in

[Executioner]
  type = Steady
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
  line_search = 'l2'
  # Linear/nonlinear iterations.
  nl_abs_tol = 1e-8
[]

References