Pronghorn thermal hydraulics steady state for gFHR

Contact: Dr. Mustafa Jaradat, [email protected]

Sponsor: Dr. Steve Bajorek (NRC)

Model summarized, documented, and uploaded by Dr. Mustafa Jaradat and Dr. Samuel Walker

Pronghorn is the thermal hydraulics solver we use to simulate fluid flow and conjugate heat transfer in the gFHR core. In the same input file, we model the heat transfer in the solid phase in the core and heat conduction in the solid components around the core, for a total of five equations: conservation of mass, x- and y-momentum (in RZ coordinates), fluid and solid energy.

We are using the finite volume capabilities in MOOSE to model all these physics. We use an weakly-compressible approximation for the fluid flow. There are various closures used in the pebble bed for the heat transfer coefficients, the drag models, which are detailed in the Pronghorn manual. For a more detailed analyses of the specific Navier-Stokes equations that are solved for FHRs please see the Mark 1 Pronghorn Model.

Problem Parameters

We first define in the input file header physical quantities such as the core geometry, hydraulic diameters of flow paths, porosities, interpolation methods, core power, fluid blocks, heat transfer coefficients, and a Forcheimer friction coefficient.

# ------------------------------------------------------------------------------
# Description:
# gFHR input
# gFHR salt simulation
# Porous Flow Weakly Compressible Media Finite Volume
# sprvel: superficial formulation of the velocities
# ------------------------------------------------------------------------------
# THIS INPUT IS THE FULL-PHYSICS VERSION
# ------------------------------------------------------------------------------

# Set the conductivity type using the Pronghorn model
pebble_bed_fluid_effective_conductivity_type=FunctorLinearPecletKappaFluid

# Include common input shared between full simulation and testing version
!include gFHR_pronghorn_ss_base.i

# Set a parameter specifically required by the FunctorLinearPecletKappaFluid conductivity type
FunctorMaterials/pebble_bed_fluid_effective_conductivity/porosity=porosity_energy
(pbfhr/gFHR/steady/gFHR_pronghorn_ss.i)

Global Parameters & Debug Block

We also define a few global parameters that will be added to every block that may use them. This is done to reduce the length of the input file and improve its readability.

Here we set the pebble diameter, fluid properties, pressure, interpolation methods, density, acceleration, and Rhie-Chow user object.

[GlobalParams]
  pebble_diameter = ${pebble_diameter}
  fp = fluid_properties
  pressure = pressure

  # FV parameters
  two_term_boundary_expansion = true
  advected_interp_method = ${advected_interp_method}
  velocity_interp_method = ${velocity_interp_method}
  rho = 'rho'
  force_define_density = true
  acceleration = '0.0 -9.8 0.0'
  rhie_chow_user_object = 'pins_rhie_chow_interpolator'
[]
(pbfhr/gFHR/steady/gFHR_pronghorn_ss.i)

Additionally, there's a very short Debug block which simply specifies that we want the variable residual norms to be showed in the output as the model runs.

[Debug]
  show_var_residual_norms = true
[]
(pbfhr/gFHR/steady/gFHR_pronghorn_ss.i)

Mesh and geometry

The mesh input is standardized across MOOSE applications. The Mesh block is different from the one we saw for Griffin since the thermal-hydraulic mesh incorporates flow paths that go beyond the core. Similarly, we use the CartesianMeshGenerator to make our mesh. Then we generate side-sets to define boundaries using SideSetsAroundSubdomainGenerator and SideSetsBetweenSubdomainsGenerator.

The thermal-hydraulic mesh is shown in Figure 1.

Figure 1: Pronghorn thermal fluids mesh from (Ortensi et al., 2023).

[Mesh]
  block_id = '  1     2     3     4     5    6    7    8    9     10           12     13      14      20      21  ' #30 is a buffer layer to protect the porosity smoothing and ensure conservation
  block_name = 'flow1 flow2 flow3 flow4 flow5 dio1 dio2 dio3 flowp reflector   barrel vessel1 vessel2 inlet outlet'
  coord_type = RZ
  rz_coord_axis = Y

  [cartesian_mesh]
    type = CartesianMeshGenerator
    dim = 2

    # keep mesh in the pebble bed roughly 3xpebble diameter ~12 cm for gFHR
    dx = '0.8 0.1 0.3 0.2 0.2 0.2 0.02 0.05 0.04'
    ix = '  6   1   3   2   2   2    1    2    1'

    dy = '0.04 0.05  0.3 0.3  0.154735 2.78523 0.154735  0.6  0.5  0.05 0.1 0.2 0.4'
    iy = '   1    2    2   2         1      18        1    4    4     2   1   2   2'

    subdomain_id = '
   14 13 13 13 13 13 13 13 13
    2  1  1  1  1  1  1  1 13
    2  2 10 10 10 10 12  1 13
    2  2  9 10 10 10 12  1 13
    3  3  3 10 10 10 12  1 13
    3  3  3 10 10 10 12  1 13
    3  3  3 10 10 10 12  1 13
    4  9  5 10 10 10 12  1 13
    4 10  5 10 10 10 12  1 13
    4 10  5  6  7  8  8  1 13
    4 10 10 10 10 10 12  1 13
    4 10 10 10 10 10 12  1 20
   21 10 10 10 10 10 12 12 13'
  []

  [CL_SS]
    type = SideSetsAroundSubdomainGenerator
    input = cartesian_mesh
    fixed_normal = true
    normal = '-1 0 0'
    block = '2 3 4 14 21'
    new_boundary = 'centerline'
  []
  [Inlet]
    type = SideSetsAroundSubdomainGenerator
    input = CL_SS
    fixed_normal = true
    normal = '1 0 0'
    block = '20'
    new_boundary = 'inlet'
  []
  [Outlet]
    type = SideSetsAroundSubdomainGenerator
    input = Inlet
    fixed_normal = true
    normal = '0 1 0'
    block = '21'
    new_boundary = 'outlet'
  []
  [vertical_walls_01]
    type = SideSetsBetweenSubdomainsGenerator
    input = Outlet
    primary_block = '1 20'
    paired_block = 13
    new_boundary = wall1
  []
  [vertical_walls_02]
    type = SideSetsBetweenSubdomainsGenerator
    input = vertical_walls_01
    primary_block = 2
    paired_block = 14
    new_boundary = wall1
  []
  [vertical_walls_03]
    type = SideSetsBetweenSubdomainsGenerator
    input = vertical_walls_02
    primary_block = '1 8'
    paired_block = 12
    new_boundary = wall1
  []
  [vertical_walls_04]
    type = SideSetsBetweenSubdomainsGenerator
    input = vertical_walls_03
    primary_block = '1 2 3 4 5 6 7 8 9 21'
    paired_block = 10
    new_boundary = wall1
  []
  [RCSS_Boundary]
    type = SideSetsAroundSubdomainGenerator
    input = vertical_walls_04
    fixed_normal = true
    include_only_external_sides = true
    block = '13 14'
    normal = '1 0 0'
    new_boundary = RCSS
  []
  [Mass_Flow_Sensor_1]
    type = SideSetsBetweenSubdomainsGenerator
    input = RCSS_Boundary
    primary_block = '6'
    paired_block = 7
    new_boundary = MF_1
  []
  [Mass_Flow_Sensor_2]
    type = SideSetsBetweenSubdomainsGenerator
    input = Mass_Flow_Sensor_1
    primary_block = '7'
    paired_block = 8
    new_boundary = MF_2
  []
  [Mass_Flow_Sensor_3]
    type = SideSetsBetweenSubdomainsGenerator
    input = Mass_Flow_Sensor_2
    primary_block = '1'
    paired_block = 2
    new_boundary = MF_3
  []
  uniform_refine = 0
[]
(pbfhr/gFHR/steady/gFHR_pronghorn_ss.i)

Physics - Navier Stokes Finite Volume Physics

Next we introduce a simplified approach to solving the thermal-hydraulic equations using the Navier Stokes Finite Volume Physics

The Physics directly set up the numerical form of the thermal-hydraulic equations that will be solved rather than making the user define every single kernel explicitly. In this case the action uses a weakly-compressible formulation of the Navier-Stokes equations and adds in porous_medium_treatment to model flow through the graphite pebbles.

Boundary conditions are set here, with no slip wall_boundaries, inlet_boundaries, and outlet_boundaries where in fluid dynamics fashion, the inlet is a fixed-velocity, and the outlet is a fixed-pressure.

The heat convection from blocks to the coolant is defined in the ambient_convection section. Friction factors using the darcy and forchheimer friction factors are also read into the equations here. Lastly, the remaining advection interpolation methods for momentum and mass are defined as well.

A Physics is set up for every equation that is being solved. The Physics communicate with each other to verify their parameters to some extent, attempting to minimize the number of errors. The equations being solved here are:

  • the mass and momentum equations by the Flow Physics

[Physics]
  [NavierStokes]
    [Flow]
      [salt]
        # General parameters
        compressibility = 'weakly-compressible'
        porous_medium_treatment = true
        block = ${blocks_fluid}
        gravity = '0 -9.81 0'

        # Initial conditions
        initial_pressure = 1e5
        initial_velocity = '1e-10 1e-10 0'

        # Variables
        velocity_variable = 'superficial_vel_x superficial_vel_y'
        pressure_variable = 'pressure'

        # Wall boundary conditions
        wall_boundaries = 'wall1'
        momentum_wall_types = 'slip'

        # Inlet boundary conditions
        inlet_boundaries = 'inlet'
        momentum_inlet_types = 'fixed-velocity'
        momentum_inlet_functors = '-0.242984118 0'

        # Outlet boundary conditions
        outlet_boundaries = 'outlet'
        momentum_outlet_types = 'fixed-pressure'
        pressure_functors = 'Outlet_Average_Pressure_Function'

        # Friction in porous media
        friction_types = 'darcy forchheimer'
        friction_coeffs = 'combined_linear combined_quadratic'
        use_friction_correction = True
        consistent_scaling = 200
        # porosity_smoothing_layers = 8

        # Numerical scheme
        momentum_advection_interpolation = 'upwind'
        mass_advection_interpolation = 'upwind'
      []
    []
  []
[]
(pbfhr/gFHR/steady/gFHR_pronghorn_ss.i)

  • the fluid energy conservation equation by the FluidHeatTransfer Physics

[Physics]
  [NavierStokes]
    [FluidHeatTransfer]
      [salt]
        fluid_temperature_variable = 'T_fluid'
        initial_temperature = 823.14
        block = ${blocks_fluid}

        energy_inlet_types = 'fixed-temperature'
        energy_inlet_functors = 'Inlet_Temperature_BC'

        # Wall convection is modeled with a volumetric term
        energy_wall_types = 'heatflux'
        energy_wall_functors = '0'

        thermal_conductivity = 'kappa'
        effective_conductivity = false

        # Porous flow parameters
        ambient_convection_blocks = ${blocks_fluid}
        ambient_convection_alpha = 'alpha'
        ambient_temperature = 'T_solid'
      []
    []
  []
[]
(pbfhr/gFHR/steady/gFHR_pronghorn_ss.i)

  • the solid energy conservation equation by the SolidHeatTransfer Physics

[Physics]
  [NavierStokes]
    [SolidHeatTransfer]
      [solid]
        # On all blocks
        block = '${blocks_bed} ${blocks_non_bed}'
        initial_temperature = 823.14

        thermal_conductivity_solid = 'eff_solid_conductivity'
        cp_solid = 'cp_s'
        rho_solid = 'rho_s'

        # Heat source
        external_heat_source = power_density
        external_heat_source_blocks = 'flow3'

        # Porous flow parameters
        ambient_convection_blocks = ${blocks_fluid}
        ambient_convection_alpha = 'alpha'
      []
    []
  []
[]
(pbfhr/gFHR/steady/gFHR_pronghorn_ss.i)

Auxiliary Variables

Next the AuxVariables are listed which can be used to define fields/variables that are not solved for. They may be used for coupling, for computing material properties, for outputting quantities of interests, and many other uses. In this input file, we define auxiliary variables for items like the Courant-Friedrichs-Lewy (CFL) condition CFL_Now, density rho_var, or heat capacity cp_var. Additionally, we have auxiliary variables that are used in Multi-app Transfers such as the power_density from the neutronics solve, and rho_var from the thermal hydraulics solve. Lastly we have auxiliary variables used for heat transfer or friction factor correlations.

[AuxVariables]
  [Volume_of_Cell]
    order = CONSTANT
    family = MONOMIAL
    block = ${blocks_fluid}
    initial_condition = 0.00012499797344389663
  []
  [CFL_Now]
    type = MooseVariableFVReal
    block = ${blocks_fluid}
  []
  [Dt]
    type = MooseVariableFVReal
    block = ${blocks_fluid}
  []
  [porosity_1]
    type = MooseVariableFVReal
    block = ${blocks_fluid}
  []
  [real_vel_x]
    type = MooseVariableFVReal
    initial_condition = 0
  []
  [real_vel_y]
    type = MooseVariableFVReal
    initial_condition = 0.0
  []
  [power_density]
    type = MooseVariableFVReal
    initial_condition = '${fparse core_power_density}'
    block = 'flow3'
  []
  [rho_var]
    type = MooseVariableFVReal
  []
  [cp_var]
    type = MooseVariableFVReal
  []
  [T_fluid_WD]
    type = MooseVariableFVReal
  []
  [Wall_Distance]
    type = MooseVariableFVReal
  []
  [htc_material]
    type = MooseVariableFVReal
  []
  [Adjusted_HTC_Material]
    type = MooseVariableFVReal
  []
  [kappa_fluid]
    type = MooseVariableFVReal
  []
  [Darcy_Value]
    type = MooseVariableFVReal
  []
  [Forchheimer_Value]
    type = MooseVariableFVReal
  []
[]
(pbfhr/gFHR/steady/gFHR_pronghorn_ss.i)

Functions

The Functions block is defined next. Here functions are defined that can be used elsewhere in the input file. Most of these function define porosity for specific locations in the model using the ParsedFunction which is a versatile way of setting up various functions.

[Functions]
  [dt_function]
    type = ParsedFunction
    symbol_values = 'dt'
    symbol_names = 'dt'
    expression = 'dt'
  []
  [porosity_core]
    type = ParsedFunction
    expression = 0.4011505932
  []
  [porosity_cones]
    type = ParsedFunction
    expression = 0.75
  []
  [porosity_plena]
    type = ParsedFunction
    expression = 0.9999
  []
  [porosity_pipe]
    type = ParsedFunction
    expression = 0.9999
  []
  [porosity_diode]
    type = ParsedFunction
    expression = 0.9999
  []
  [Outlet_Average_Pressure_Function]
    type = ParsedFunction
    symbol_names = 'OAP'
    expression = OAP
    symbol_values = 'Outlet_Pressure_BC'
  []
  [mu_rampdown]
    type = ParsedFunction
    expression = 1
    symbol_names = 'NP'
    symbol_values = 'Num_Picards'
  []
[]
(pbfhr/gFHR/steady/gFHR_pronghorn_ss.i)

Auxiliary Kernels

Similarly we can have Auxiliary Kernels which operate on the Auxiliary Variables. Similarly - the ParsedAux is a powerful tool in performing various actions that operate on the Auxiliary Variables. Specifically - heat transfer, turbulence, and friction factor correlations can be implemented here.

[AuxKernels]
  [volume_Calc]
    type = VolumeAux
    variable = Volume_of_Cell
    block = ${blocks_fluid}
  []
  [dt_kernel]
    type = FunctionAux
    variable = Dt
    function = dt_function
  []
  [CFL_Calc]
    type = ParsedAux
    variable = CFL_Now
    coupled_variables = 'superficial_vel_x superficial_vel_y Volume_of_Cell Dt'
    expression = '(abs(superficial_vel_x) + abs(superficial_vel_y))*Dt/(Volume_of_Cell)^(0.5)'
    block = ${blocks_fluid}
  []
  [real_vel_x_calc]
    type = ParsedAux
    variable = real_vel_x
    coupled_variables = 'superficial_vel_x porosity_1'
    expression = 'superficial_vel_x / porosity_1'
    block = ${blocks_fluid}
  []
  [real_vel_y_calc]
    type = ParsedAux
    variable = real_vel_y
    coupled_variables = 'superficial_vel_y porosity_1'
    expression = 'superficial_vel_y / porosity_1'
    block = ${blocks_fluid}
  []
  [rho_out]
    type = FunctorAux
    functor = 'rho'
    variable = 'rho_var'
    execute_on = 'timestep_begin timestep_end Final'
    block = ${blocks_fluid}
  []
  [cp_out]
    type = FunctorAux
    functor = 'cp'
    variable = 'cp_var'
    execute_on = 'timestep_begin'
    block = ${blocks_fluid}
  []
  [Whole_Domain_Fluid_Temp]
    type = ParsedAux
    variable = T_fluid_WD
    coupled_variables = 'T_fluid'
    expression = 'T_fluid'
    block = ${blocks_fluid}
  []
  [Whole_Domain_Fluid_Temp_2]
    type = ConstantAux
    variable = T_fluid_WD
    value = 0
    block = 'reflector barrel vessel1 vessel2'
  []
  [porosity_out]
    type = FunctorAux
    functor = 'porosity'
    variable = 'porosity_1'
    block = ${blocks_fluid}
    execute_on = 'initial'
  []
  [kappa_out]
    type = ADFunctorVectorElementalAux
    functor = 'kappa'
    variable = 'kappa_fluid'
    component = '0'
    execute_on = 'timestep_end final'
    block = ${blocks_fluid}
  []
  [Wall_Distance_Aux_Kernel]
    type = WallDistanceMixingLengthAux
    variable = Wall_Distance
    von_karman_const = 1
    von_karman_const_0 = 1
    walls = 'wall1'
  []
  [Adjust_HTC]
    type = ParsedAux
    variable = Adjusted_HTC_Material
    coupled_variables = 'htc_material Wall_Distance'
    expression = 'htc_material*Wall_Distance'
    block = ${blocks_fluid}
  []
  [d_var_aux]
    type = ADFunctorVectorElementalAux
    variable = Darcy_Value
    functor = combined_linear
    block = ${blocks_fluid}
    component = '0'
  []
  [f_var_aux]
    type = ADFunctorVectorElementalAux
    variable = Forchheimer_Value
    functor = combined_quadratic
    block = ${blocks_fluid}
    component = '0'
  []
[]
(pbfhr/gFHR/steady/gFHR_pronghorn_ss.i)

Fluid Properties

Very briefly we introduce the Fluid Properties which imports thermophysical properties at evaluates this data at the temperature and pressure that is calculated by the thermal hydraulic solution. Here FLiBe molten salt thermophysical properties are loaded in.

[FluidProperties]
  [fluid_properties]
    type = FlibeFluidProperties
  []
[]
(pbfhr/gFHR/steady/gFHR_pronghorn_ss.i)

Functor Materials

Next we define Functor Materials which define specific materials which are functors and can be used and incorporated within the Kernels or Auxiliary Kernels to perform specific calculations.

A few examples here are the GeneralFunctorFluidProps which performs the evaluation of thermophysical properties of the FliBe salt from the FluidProperties and the thermal hydraulic variables that are being solved for (i.e., T_fluid, pressure, speed). Additionally, the characteristic length or the hydraulic diameter, porosity, and friction factors of flow regions in the model can also be defined here and used in the Action or FVKernels.

[FunctorMaterials]
  [GFFP]
    type = GeneralFunctorFluidProps
    T_fluid = 'T_fluid'
    characteristic_length = characteristic_length
    fp = fluid_properties
    porosity = 'porosity'
    pressure = 'pressure'
    speed = speed
    mu_rampdown = 'mu_rampdown'
    block = ${blocks_fluid}
  []

  ## Characteristic length
  [characteristic_length]
    type = PiecewiseByBlockFunctorMaterial
    prop_name = characteristic_length
    block = ${blocks_fluid}
    subdomain_to_prop_value = '1  ${D_H_downcomer}
                               2  ${D_H_flow2}
                               3  ${pebble_diameter}
                               4  ${pebble_diameter}
                               5  ${D_H_flow5}
                               6  ${D_H_diode}
                               7  ${D_H_diode}
                               8  ${D_H_diode}
                               9  ${pebble_diameter}
                               20 ${D_H_downcomer}
                               21 ${D_H_outlet}'
  []

  ## Porosity
  [porosity_energy]
    type = ADPiecewiseByBlockFunctorMaterial
    prop_name = porosity_energy
    subdomain_to_prop_value = '1 ${porosity_free_flow}
                               2 ${porosity_free_flow}
                               3 ${porosity_bed}
                               4 ${porosity_cone}
                               5 ${porosity_free_flow}
                               6 ${porosity_free_flow}
                               7 ${porosity_free_flow}
                               8 ${porosity_free_flow}
                               9 ${porosity_cone}
                               20 ${porosity_free_flow}
                               21 ${porosity_free_flow}
                               10 ${porosity_solid}
                               12 ${porosity_solid}
                               13 ${porosity_solid}
                               14 ${porosity_solid}'
  []

  [porosity]
    type = ADPiecewiseByBlockFunctorMaterial
    prop_name = porosity
    subdomain_to_prop_value = '1 ${porosity_free_flow}
                               2 ${porosity_free_flow}
                               3 ${porosity_bed}
                               4 ${porosity_cone}
                               5 ${porosity_free_flow}
                               6 ${porosity_free_flow}
                               7 ${porosity_free_flow}
                               8 ${porosity_free_flow}
                               9 ${porosity_cone}
                               20 ${porosity_free_flow}
                               21 ${porosity_free_flow}
                               10 ${porosity_solid}
                               12 ${porosity_solid}
                               13 ${porosity_solid}
                               14 ${porosity_solid}'
  []

  [diode]
    type = NSFVFrictionFlowDiodeFunctorMaterial
    direction = '-1 0 0'
    additional_linear_resistance = '1000 0 0'
    additional_quadratic_resistance = '100 0 0'
    base_linear_friction_coefs = 'Darcy_coefficient'
    base_quadratic_friction_coefs = 'Forchheimer_coefficient'
    sum_linear_friction_name = 'diode_linear'
    sum_quadratic_friction_name = 'diode_quad'
    block = '6 7 8'
    turn_on_diode = true
  []

  [combine_linear_friction]
    type = ADPiecewiseByBlockVectorFunctorMaterial
    prop_name = 'combined_linear'
    subdomain_to_prop_value = '1 Darcy_coefficient
                               2 Darcy_coefficient
                               3 Darcy_coefficient
                               4 Darcy_coefficient
                               5 Darcy_coefficient
                               6 diode_linear
                               7 diode_linear
                               8 diode_linear
                               9 Darcy_coefficient
                               20 Darcy_coefficient
                               21 Darcy_coefficient'
  []

  [combine_quadratic_friction]
    type = ADPiecewiseByBlockVectorFunctorMaterial
    prop_name = 'combined_quadratic'
    subdomain_to_prop_value = '1 Forchheimer_coefficient
                               2 Forchheimer_coefficient
                               3 Forchheimer_coefficient
                               4 Forchheimer_coefficient
                               5 Forchheimer_coefficient
                               6 diode_linear
                               7 diode_linear
                               8 diode_linear
                               9 Forchheimer_coefficient
                               20 Forchheimer_coefficient
                               21 Forchheimer_coefficient'
  []

  ## Thermal conductivity of the fluid
  [pebble_bed_fluid_effective_conductivity]
    type = ${pebble_bed_fluid_effective_conductivity_type}
    block = ${blocks_bed}
    porosity = porosity_energy
  []

  [non_pebble_bed_kappa]
    type = ADGenericVectorFunctorMaterial
    prop_names = 'kappa'
    prop_values = 'k k k'
    block = ${blocks_non_bed}
  []

  ## Drag coefficients
  [pebble_bed_drag_coefficient]
    type = FunctorKTADragCoefficients
    T_solid = T_solid
    T_fluid = T_fluid
    porosity = porosity
    block = ${blocks_bed}
  []

  [churchill_drag_coefficient]
    type = FunctorChurchillDragCoefficients
    block = 'flow1 inlet dio1 dio2 dio3'
  []

  [drag_open_flow_regions]
    type = FunctorIsotropicFunctorDragCoefficients
    Darcy_coefficient_sub = 0
    Forchheimer_coefficient_sub = ${Forch_open}
    block = 'flow2 flow5 outlet'
  []

  ## Heat transfer coefficients
  [pebble_bed_alpha]
    type = FunctorPetrovicPebbleBedHTC
    block = ${blocks_bed}
    porosity = porosity_energy
  []

  [equivalent_alpha_flow1]
    type = FunctorDittusBoelterWallHTC
    C = '${fparse C_DB * ApV_flow_1}'
    block = flow1
  []

  [equivalent_alpha_diodes]
    type = FunctorDittusBoelterWallHTC
    C = '${fparse C_DB * ApV_diode}'
    block = 'dio1 dio2 dio3'
  []

  [convert_wall_htc_to_alpha]
    type = ADGenericFunctorMaterial
    prop_names = 'alpha'
    prop_values = 'wall_htc'
    block = 'dio1 dio2 dio3 flow1'
  []

  # in some regions we do not care about volumetric
  # heat transfer, but we might have to do sanity checks here
  [alpha_null]
    type = ADGenericFunctorMaterial
    prop_names = 'alpha'
    prop_values = '0'
    block = 'flow2 flow5 outlet  inlet'
  []

  ## Solid phase properties
  [pebble_bed_kappa_s]
    type = FunctorPebbleBedKappaSolid
    radiation = BreitbachBarthels
    solid_conduction = ChanTien
    fluid_conduction = ZBS
    emissivity = ${pebble_emissivity}
    Youngs_modulus = 9e9
    Poisson_ratio = 0.136
    wall_distance = bed_geometry
    T_solid = T_solid
    block = ${blocks_bed}
    porosity = porosity_energy
  []

  [pebble_bed_solid_properties]
    type = PronghornSolidFunctorMaterialPT
    solid = graphite
    T_solid = T_solid
    block = ${blocks_bed}
  []

  [steel]
    type = ADGenericFunctorMaterial
    block = 'barrel vessel1 vessel2'
    prop_names = 'rho_s cp_s k_s'
    prop_values = '8000.0 ${fparse 500.0} 21.5'
  []

  [reflector_full_density_graphite]
    type = ADGenericFunctorMaterial
    block = 'reflector'
    prop_names = 'rho_s cp_s k_s'
    prop_values = '1740.0 ${fparse 1697.0} 26.0'
  []

  [time_derivatives_constants]
    type = ADGenericFunctorMaterial
    prop_names = 'drho_s_dt'
    prop_values = '0'
    block = '${blocks_bed}'
  []

  # Modeling solid everywhere. Open flow regions are modeled
  # as very low density/conductivity regions. Heat transfer
  # through them will be facilitated by conjugate heat transfer
  # from the adjacent walls that is converted to a volumetric term
  [dummy_open_flow_region_solid_props]
    type = ADGenericFunctorMaterial
    prop_names = 'rho_s cp_s k_s'
    prop_values = '1 1 1e-4'
    block = '${blocks_free_flow}'
  []

  # write k_s into kappa_s so we can use the same
  # kernel everywhere
  [non_pebble_bed_eff_solid_conductivity]
    type = ADGenericVectorFunctorMaterial
    prop_names = 'eff_solid_conductivity'
    prop_values = 'k_s k_s k_s'
    block = '${blocks_free_flow} reflector barrel vessel1 vessel2'
  []

  [pebble_bed_eff_solid_conductivity]
    type = ADGenericVectorFunctorMaterial
    prop_names = 'eff_solid_conductivity'
    prop_values = 'kappa_s kappa_s kappa_s'
    block = '${blocks_bed}'
  []
[]
(pbfhr/gFHR/steady/gFHR_pronghorn_ss.i)

Postprocessors

Nearing the end of the model, the Postprocessors block specifies calculations that should be made regarding the solution that may be used during picard iterations or for data analysis after convergence. There are many different types used in this model including ElementExtremeValue, ParsedPostprocessor which also gives the user a wide berth of possibilities to calculate specific values of interest, ElementIntegralVariablePostprocessor, and SideAverageValue to name a few. Here both postprocessors for the thermal-hydraulic solve are listed which operate on the pressure, velocity, and temperatures solved for in the numerical simulation, as well as values that were transferred in from Griffin like the power_density.

[Postprocessors]
  [dt]
    type = TimestepSize
  []
  [Max_CFL]
    type = ElementExtremeValue
    variable = CFL_Now
    block = ${blocks_fluid}
  []
  [New_Dt]
    type = ParsedPostprocessor
    pp_names = 'Max_CFL dt'
    #Relaxed and limited growth timestep so that we dont end up in a cyclical problem
    expression = 'min((10000*t/Max_CFL * dt/2 + dt/2) , dt*2)'
    execute_on = 'timestep_end'
    use_t = true
  []
  [Num_Picards]
    type = Receiver
    default = 0
  []

  [Inlet_Temperature_BC]
    type = Receiver
    default = 823.15
    execute_on = 'INITIAL TIMESTEP_END TRANSFER'
  []
  [Outlet_Pressure_BC]
    type = Receiver
    default = 2e5
    execute_on = 'INITIAL TIMESTEP_END TRANSFER'
  []
  [Outlet_Velocity_y]
    type = SideAverageValue
    boundary = 'outlet'
    variable = superficial_vel_y
    execute_on = 'INITIAL TIMESTEP_END TRANSFER'
  []

  ### Griffin values
  [total_power]
    type = ElementIntegralVariablePostprocessor
    block = 'flow3'
    variable = power_density
    execute_on = 'INITIAL TIMESTEP_END TRANSFER'
  []

  ### Pronghorn Values ###
  [Outlet_Velocity_PH]
    type = ParsedPostprocessor
    expression = 'mass_flow_out / area_SAM / density'
    pp_names = 'mass_flow_out'
    constant_names = 'area_SAM density'
    constant_expressions = '1.327511 2201'
    execute_on = 'INITIAL TIMESTEP_END TRANSFER'
  []
  [Outlet_Pressure_PH]
    type = SideAverageValue
    boundary = 'outlet'
    variable = pressure
    execute_on = 'INITIAL TIMESTEP_END TRANSFER'
  []
  [Outlet_Temperature_PH]
    type = SideAverageValue
    boundary = 'outlet'
    variable = T_fluid
    execute_on = 'INITIAL TIMESTEP_END TRANSFER'
  []
  [Heatflux_RCCS]
    type = SideIntegralVariablePostprocessor
    boundary = 'RCSS'
    variable = T_solid
  []
  [Tfluid_core]
    type = ElementAverageValue
    block = 'flow3'
    variable = T_fluid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tfluid_bottom]
    type = ElementAverageValue
    block = 'flow2'
    variable = T_fluid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tfluid_upper]
    type = ElementAverageValue
    block = 'flow4 flow5'
    variable = T_fluid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tfluid_diode]
    type = ElementAverageValue
    block = 'dio1 dio2 dio3'
    variable = T_fluid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Tfluid_downcomer]
    type = ElementAverageValue
    block = 'flow1'
    variable = T_fluid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Reflector_Temp]
    type = ElementAverageValue
    block = 'reflector'
    variable = T_solid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Barrel_Temp]
    type = ElementAverageValue
    block = 'barrel'
    variable = T_solid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [Vessel_Temp]
    type = ElementAverageValue
    block = 'vessel1 vessel2'
    variable = T_solid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  ### ###
  [inlet_mdot]
    type = Receiver
    default = 1173.0
    execute_on = 'INITIAL TIMESTEP_END TRANSFER'
  []
  [inlet_area]
    type = AreaPostprocessor
    boundary = 'inlet'
    execute_on = 'INITIAL TIMESTEP_END TRANSFER'
  []
  [mass_flow_out]
    type = VolumetricFlowRate
    boundary = 'outlet'
    advected_quantity = 'rho_var'
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    execute_on = 'INITIAL TIMESTEP_END TRANSFER'
  []
  [mass_flow_in]
    type = VolumetricFlowRate
    boundary = 'inlet'
    advected_quantity = 'rho_var'
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    execute_on = 'INITIAL TIMESTEP_END TRANSFER'
  []
  [mass_flow_1]
    type = VolumetricFlowRate
    boundary = 'MF_1'
    advected_quantity = 'rho_var'
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    execute_on = 'TIMESTEP_END TRANSFER'
  []
  [mass_flow_2]
    type = VolumetricFlowRate
    boundary = 'MF_2'
    advected_quantity = 'rho_var'
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    execute_on = 'TIMESTEP_END TRANSFER'
  []
  [mass_flow_3]
    type = VolumetricFlowRate
    boundary = 'MF_3'
    advected_quantity = 'rho_var'
    vel_x = superficial_vel_x
    vel_y = superficial_vel_y
    execute_on = 'TIMESTEP_END TRANSFER'
  []
[]
(pbfhr/gFHR/steady/gFHR_pronghorn_ss.i)

Preconditioning

Next we declare the Preconditioning block which specifies what kind of preconditioning matrix the user wants to implement. Here the Single Matrix Preconditioner (SMP) is used with specific PETSc options selected.

[Preconditioning]
  [SMP]
    type = SMP
    full = true
    petsc_options_iname = ' -pc_type -pc_factor_mat_solver_type
                            -ksp_gmres_restart -pc_factor_shift_type'
    petsc_options_value = ' lu superlu_dist 50 NONZERO'
  []
[]
(pbfhr/gFHR/steady/gFHR_pronghorn_ss.i)

Executioner - Solving the Equations

The Executioner block defines how the non-linear system will be solved. Since this is a transient problem, we are using a Transient executioner. Pronghorn makes use of automatic differentiation, more information here, to compute an exact Jacobian, so we make use of the Newton method to solve the non-linear system.

We then specify the solver tolerances. These may be loosened to obtain a faster solve. They should initially be set fairly small to obtain a tight convergence, then relaxed to improve performance.

Finally we set the time parameters for the simulation. We choose a very large runtime. We use an adaptive time-stepper to increase the time step over the course of the simulation which uses a postprocessor to define what the new time step should be. Because we are using an implicit scheme, we can take very large time steps and still have a stable simulation.

[Executioner]
  type = Transient # Pseudo transient to reach steady state.
  solve_type = 'NEWTON'
  petsc_options = '-snes_converged_reason'
  line_search = 'none'

  # Problem time parameters.
  end_time = 1e+6
  dtmin = 1e-6
  dtmax = 10000

  # Iterations parameters.
  l_max_its = 50
  l_tol = 1e-6

  nl_max_its = 25
  nl_rel_tol = 1e-10
  nl_abs_tol = 1e-11
  nl_forced_its = 1

  # Steady state detection.
  steady_state_detection = true
  steady_state_tolerance = 1e-12
  steady_state_start_time = 0

  auto_advance = true

  # Adaptive time step scheme.
  [TimeStepper]
    type = PostprocessorDT
    dt = 1e-1
    postprocessor = New_Dt
  []

  # Automatic scaling
  automatic_scaling = true
  off_diagonals_in_auto_scaling = true
  compute_scaling_once = false
[]
(pbfhr/gFHR/steady/gFHR_pronghorn_ss.i)

Outputs

Lastly we define the output forms for the model which will be filled with the information listed in the Postprocessor block and the Variables and Auxiliary Variables blocks.

We then define an Exodus output. This will have the multi-dimensional distributions of the quantities Pronghorn solved for: the fluid velocities, the pressure and the temperature fields. We also include the material properties to help us understand the behavior of the core.

We also define a CSV output which will list all of the postprocessors for each time step.

[Outputs]
  print_linear_converged_reason = false
  print_linear_residuals = false
  print_nonlinear_converged_reason = false
  perf_graph = true
  csv = true
  execute_on = 'INITIAL FINAL'
  [console]
    type = Console
    outlier_variable_norms = false
    all_variable_norms = false
  []
  [exo]
    type = Exodus
    execute_on = 'FINAL'
  []
[]
(pbfhr/gFHR/steady/gFHR_pronghorn_ss.i)

References

  1. Javier Ortensi, Cole M. Mueller, Stefano Terlizzi, Guillaume Giudicelli, and Sebastian Schunert. Fluoride-cooled high-temperature pebble-bed reactor reference plant model. Technical Report INL/RPT-23-72727, Idaho National Laboratory, 2023.[BibTeX]