Step 3
In this step, we make the flow problem heterogeneous. Above the pebble bed in a pebble-bed reactor is usually a cavity that is only filled with fluid (i.e., it has a porosity of 1). The interface between the cavity and the bed is modeled as a sharp discontinuity in both porosity and correlations (i.e., momentum and energy transfer; e.g., very small pressure drop in the cavity and much larger pressure drop in the bed).
Modifying the Geometry
We distinguish the two regions, cavity and bed, by assigning the elements in each to different blocks. This is accomplished in the mesh generation process.
Instead of using GeneratedMeshGenerator
, we use the CartesianMeshGenerator
.
[Mesh]
[gen]
type = CartesianMeshGenerator
dim = 2
dx = '${bed_radius}'
ix = '6'
dy = '${bed_height} ${cavity_height}'
iy = '40 2'
subdomain_id = '1 2'
[]
[]
(htgr/generic-pbr-tutorial/step3.i)Note that we added entries in subdomain_id
. The order of entries in the subdomain_id
parameter is:
(0,0) (1,0), (2,0), ... , (0,1), (1,1), ... (nx,ny)
where each (ix,iy)
is one entry that gives the subdomain (or synonymously block id) of zone located at grid coordinates ix
and iy
. The indexing goes from small to large coordinate values, so we fill the subdomain_id
from the left bottom to the right top.
For convenience, we label block 1 as bed and block 2 as cavity using the RenameBlockGenerator
:
[Mesh]
[rename_blocks]
type = RenameBlockGenerator
old_block = '1 2'
new_block = 'bed cavity'
input = gen
[]
[]
(htgr/generic-pbr-tutorial/step3.i)That allows us to use block names rather than ids in the block restrictions of the input file.
The geometry for Step 3 is depicted in Figure 1. The mesh can be generated by:
./pronghorn-opt -i step3.i --mesh-only
creating the file step3_in.e
that can be visualized using paraview.

Figure 1: Geometry for Step 3.
(htgr/generic-pbr-tutorial/step3.i)
# ==============================================================================
# Model description:
# Step3 - Step2 with cavity on top.
# ------------------------------------------------------------------------------
# Idaho Falls, INL, August 15, 2023 04:03 PM
# Author(s): Joseph R. Brennan, Dr. Sebastian Schunert, Dr. Mustafa K. Jaradat
# and Dr. Paolo Balestra.
# ==============================================================================
bed_height = 10.0
bed_radius = 1.2
cavity_height = 0.5
bed_porosity = 0.39
outlet_pressure = 5.5e6
T_fluid = 300
density = 8.6545
pebble_diameter = 0.06
mass_flow_rate = 60.0
flow_area = '${fparse pi * bed_radius * bed_radius}'
flow_vel = '${fparse mass_flow_rate / flow_area / density}'
[Mesh]
[gen]
type = CartesianMeshGenerator
dim = 2
dx = '${bed_radius}'
ix = '6'
dy = '${bed_height} ${cavity_height}'
iy = '40 2'
subdomain_id = '1 2'
[]
[rename_blocks]
type = RenameBlockGenerator
old_block = '1 2'
new_block = 'bed cavity'
input = gen
[]
coord_type = RZ
[]
[FluidProperties]
[fluid_properties_obj]
type = HeliumFluidProperties
[]
[]
[Modules]
[NavierStokesFV]
# general control parameters
compressibility = 'weakly-compressible'
porous_medium_treatment = true
# material property parameters
density = rho
dynamic_viscosity = mu
# porous medium treatment parameters
porosity = porosity
porosity_interface_pressure_treatment = 'bernoulli'
# initial conditions
initial_velocity = '1e-6 1e-6 0'
initial_pressure = 5.4e6
# boundary conditions
inlet_boundaries = top
momentum_inlet_types = fixed-velocity
momentum_inlet_function = '0 -${flow_vel}'
wall_boundaries = 'left right'
momentum_wall_types = 'slip slip'
outlet_boundaries = bottom
momentum_outlet_types = fixed-pressure
pressure_function = ${outlet_pressure}
# friction control parameters
friction_types = 'darcy forchheimer'
friction_coeffs = 'Darcy_coefficient Forchheimer_coefficient'
[]
[]
[FunctorMaterials]
[fluid_props_to_mat_props]
type = GeneralFunctorFluidProps
fp = fluid_properties_obj
porosity = porosity
pressure = pressure
T_fluid = ${T_fluid}
speed = speed
characteristic_length = ${pebble_diameter}
[]
[drag_pebble_bed]
type = FunctorKTADragCoefficients
fp = fluid_properties_obj
pebble_diameter = ${pebble_diameter}
porosity = porosity
T_fluid = ${T_fluid}
T_solid = ${T_fluid}
block = bed
[]
[drag_cavity]
type = ADGenericVectorFunctorMaterial
prop_names = 'Darcy_coefficient Forchheimer_coefficient'
prop_values = '0 0 0 0 0 0'
block = cavity
[]
[porosity_material]
type = ADPiecewiseByBlockFunctorMaterial
prop_name = porosity
subdomain_to_prop_value = 'bed ${bed_porosity}
cavity 1'
[]
[]
[Executioner]
type = Transient
end_time = 100
[TimeStepper]
type = IterationAdaptiveDT
iteration_window = 2
optimal_iterations = 8
cutback_factor = 0.8
growth_factor = 2
dt = 1e-3
[]
dtmax = 5
line_search = l2
solve_type = 'NEWTON'
petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu NONZERO superlu_dist'
nl_rel_tol = 1e-8
nl_abs_tol = 1e-5
nl_max_its = 15
[]
[Postprocessors]
[inlet_mfr]
type = VolumetricFlowRate
advected_quantity = rho
vel_x = 'superficial_vel_x'
vel_y = 'superficial_vel_y'
boundary = 'top'
rhie_chow_user_object = pins_rhie_chow_interpolator
[]
[outlet_mfr]
type = VolumetricFlowRate
advected_quantity = rho
vel_x = 'superficial_vel_x'
vel_y = 'superficial_vel_y'
boundary = 'bottom'
rhie_chow_user_object = pins_rhie_chow_interpolator
[]
[inlet_pressure]
type = SideAverageValue
variable = pressure
boundary = top
outputs = none
[]
[outlet_pressure]
type = SideAverageValue
variable = pressure
boundary = bottom
outputs = none
[]
[pressure_drop]
type = ParsedPostprocessor
pp_names = 'inlet_pressure outlet_pressure'
expression = 'inlet_pressure - outlet_pressure'
[]
[integral_density]
type = ADElementIntegralFunctorPostprocessor
functor = rho
outputs = none
[]
[average_density]
type = ParsedPostprocessor
pp_names = 'volume integral_density'
expression = 'integral_density / volume'
[]
[integral_mu]
type = ADElementIntegralFunctorPostprocessor
functor = mu
outputs = none
[]
[average_mu]
type = ParsedPostprocessor
pp_names = 'volume integral_mu'
expression = 'integral_mu / volume'
[]
[area]
type = AreaPostprocessor
boundary = bottom
outputs = none
[]
[volume]
type = VolumePostprocessor
[]
[]
[Outputs]
exodus = true
[]
(htgr/generic-pbr-tutorial/step3.i)
# ==============================================================================
# Model description:
# Step3 - Step2 with cavity on top.
# ------------------------------------------------------------------------------
# Idaho Falls, INL, August 15, 2023 04:03 PM
# Author(s): Joseph R. Brennan, Dr. Sebastian Schunert, Dr. Mustafa K. Jaradat
# and Dr. Paolo Balestra.
# ==============================================================================
bed_height = 10.0
bed_radius = 1.2
cavity_height = 0.5
bed_porosity = 0.39
outlet_pressure = 5.5e6
T_fluid = 300
density = 8.6545
pebble_diameter = 0.06
mass_flow_rate = 60.0
flow_area = '${fparse pi * bed_radius * bed_radius}'
flow_vel = '${fparse mass_flow_rate / flow_area / density}'
[Mesh]
[gen]
type = CartesianMeshGenerator
dim = 2
dx = '${bed_radius}'
ix = '6'
dy = '${bed_height} ${cavity_height}'
iy = '40 2'
subdomain_id = '1 2'
[]
[rename_blocks]
type = RenameBlockGenerator
old_block = '1 2'
new_block = 'bed cavity'
input = gen
[]
coord_type = RZ
[]
[FluidProperties]
[fluid_properties_obj]
type = HeliumFluidProperties
[]
[]
[Modules]
[NavierStokesFV]
# general control parameters
compressibility = 'weakly-compressible'
porous_medium_treatment = true
# material property parameters
density = rho
dynamic_viscosity = mu
# porous medium treatment parameters
porosity = porosity
porosity_interface_pressure_treatment = 'bernoulli'
# initial conditions
initial_velocity = '1e-6 1e-6 0'
initial_pressure = 5.4e6
# boundary conditions
inlet_boundaries = top
momentum_inlet_types = fixed-velocity
momentum_inlet_function = '0 -${flow_vel}'
wall_boundaries = 'left right'
momentum_wall_types = 'slip slip'
outlet_boundaries = bottom
momentum_outlet_types = fixed-pressure
pressure_function = ${outlet_pressure}
# friction control parameters
friction_types = 'darcy forchheimer'
friction_coeffs = 'Darcy_coefficient Forchheimer_coefficient'
[]
[]
[FunctorMaterials]
[fluid_props_to_mat_props]
type = GeneralFunctorFluidProps
fp = fluid_properties_obj
porosity = porosity
pressure = pressure
T_fluid = ${T_fluid}
speed = speed
characteristic_length = ${pebble_diameter}
[]
[drag_pebble_bed]
type = FunctorKTADragCoefficients
fp = fluid_properties_obj
pebble_diameter = ${pebble_diameter}
porosity = porosity
T_fluid = ${T_fluid}
T_solid = ${T_fluid}
block = bed
[]
[drag_cavity]
type = ADGenericVectorFunctorMaterial
prop_names = 'Darcy_coefficient Forchheimer_coefficient'
prop_values = '0 0 0 0 0 0'
block = cavity
[]
[porosity_material]
type = ADPiecewiseByBlockFunctorMaterial
prop_name = porosity
subdomain_to_prop_value = 'bed ${bed_porosity}
cavity 1'
[]
[]
[Executioner]
type = Transient
end_time = 100
[TimeStepper]
type = IterationAdaptiveDT
iteration_window = 2
optimal_iterations = 8
cutback_factor = 0.8
growth_factor = 2
dt = 1e-3
[]
dtmax = 5
line_search = l2
solve_type = 'NEWTON'
petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu NONZERO superlu_dist'
nl_rel_tol = 1e-8
nl_abs_tol = 1e-5
nl_max_its = 15
[]
[Postprocessors]
[inlet_mfr]
type = VolumetricFlowRate
advected_quantity = rho
vel_x = 'superficial_vel_x'
vel_y = 'superficial_vel_y'
boundary = 'top'
rhie_chow_user_object = pins_rhie_chow_interpolator
[]
[outlet_mfr]
type = VolumetricFlowRate
advected_quantity = rho
vel_x = 'superficial_vel_x'
vel_y = 'superficial_vel_y'
boundary = 'bottom'
rhie_chow_user_object = pins_rhie_chow_interpolator
[]
[inlet_pressure]
type = SideAverageValue
variable = pressure
boundary = top
outputs = none
[]
[outlet_pressure]
type = SideAverageValue
variable = pressure
boundary = bottom
outputs = none
[]
[pressure_drop]
type = ParsedPostprocessor
pp_names = 'inlet_pressure outlet_pressure'
expression = 'inlet_pressure - outlet_pressure'
[]
[integral_density]
type = ADElementIntegralFunctorPostprocessor
functor = rho
outputs = none
[]
[average_density]
type = ParsedPostprocessor
pp_names = 'volume integral_density'
expression = 'integral_density / volume'
[]
[integral_mu]
type = ADElementIntegralFunctorPostprocessor
functor = mu
outputs = none
[]
[average_mu]
type = ParsedPostprocessor
pp_names = 'volume integral_mu'
expression = 'integral_mu / volume'
[]
[area]
type = AreaPostprocessor
boundary = bottom
outputs = none
[]
[volume]
type = VolumePostprocessor
[]
[]
[Outputs]
exodus = true
[]