Pronghorn thermal hydraulics steady state for Mark 1

Pronghorn is the thermal hydraulics solver we use to simulate fluid flow and conjugate heat transfer in the FHR 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), fluid and solid energy.

We are using the finite volume capabilities in MOOSE to model all these physics. We use an incompressible approximation for the fluid flow, with a Boussinesq approximation to model buoyancy. There are various closures used in the pebble bed for the heat transfer coefficients, the drag models, which are detailed in the Pronghorn manual.

We first define in the input file header physical quantities such as the pebble geometry, the material compositions as well as some constant material properties.

blocks_fluid = '3 4 5 6'
blocks_solid = '1 2 6 7 8 9 10'

# Material compositions
UO2_phase_fraction = 1.20427291e-01
buffer_phase_fraction = 2.86014816e-01
ipyc_phase_fraction = 1.59496539e-01
sic_phase_fraction = 1.96561801e-01
opyc_phase_fraction = 2.37499553e-01
TRISO_phase_fraction = 3.09266232e-01
core_phase_fraction = 5.12000000e-01
fuel_matrix_phase_fraction = 3.01037037e-01
shell_phase_fraction = 1.86962963e-01

# FLiBe properties #TODO: Rely only on PronghornFluidProps
fluid_mu = 7.5e-3 # Pa.s at 900K [3]
# k_fluid =  1.1     # suggested constant [3]
# cp_fluid = 2385    # suggested constant [3]
rho_fluid = 1970.0 # kg/m3 at 900K [3]
alpha_b = 2e-4 # /K from [4]

# Graphite properties
heat_capacity_multiplier = 1e0 # >1 gets faster to steady state
solid_rho = 1780.0
# solid_k = 26.0
solid_cp = '${fparse 1697.0*heat_capacity_multiplier}'

# Core geometry
pebble_diameter = 0.03
bed_porosity = 0.4
IR_porosity = 0
OR_porosity = 0.1123
plenum_porosity = 0.5
model_inlet_rin = 0.45
model_inlet_rout = 0.8574
model_vol = 10.4
model_inlet_area = '${fparse 3.14159265 * (model_inlet_rout * model_inlet_rout - model_inlet_rin * model_inlet_rin)}'

# The convention for friction factors changed
darcy_conversion_plenum = '${fparse rho_fluid / plenum_porosity / fluid_mu}'
# The porosity simplfies out from the previous definition of the coefficients
darcy_conversion_or = '${fparse rho_fluid / fluid_mu}'
forchheimer_partial_conversion = 2 # divide further by speed

# Outer reflector drag parameters
# TODO: tune using CFD
# TODO: verify current values (imported from [1])
Ah = '${fparse 1337.76 * darcy_conversion_or}'
Bh = '${fparse 2.58 * forchheimer_partial_conversion}'
Av = '${fparse 599.30 * darcy_conversion_or}'
Bv = '${fparse 0.95 * forchheimer_partial_conversion}'
Dh = 0.02582
Dv = 0.02006

# Plenum drag parameters. The core hot legs cannot be represented in 2D RZ
# accurately, so this should be tuned to obtain the desired mass flow rates
plenum_friction = '${fparse 0.5 * darcy_conversion_plenum}'

# Operating parameters
mfr = 976.0 # kg/s, from [2]
total_power = 236.0e6 # W, from [2]
inlet_T_fluid = 873.15 # K, from [2]
inlet_vel_y = '${fparse mfr / model_inlet_area / rho_fluid}' # superficial
(pbfhr/mark1/steady/ss1_combined.i)

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.

[GlobalParams]
  rho = ${rho_fluid}
  porosity = porosity_viz
  characteristic_length = 0.01 #${pebble_diameter}
  pebble_diameter = ${pebble_diameter}
  speed = 'speed'

  fp = fp
  T_solid = T_solid
  T_fluid = T_fluid
  pressure = pressure

  vel_x = 'superficial_vel_x'
  vel_y = 'superficial_vel_y'

  fv = true
  rhie_chow_user_object = 'pins_rhie_chow_interpolator'
[]
(pbfhr/mark1/steady/ss1_combined.i)

Mesh and geometry

The mesh input is standardized across MOOSE applications. The Mesh block is similar to the one we saw for Griffin. We load a different mesh file for this application. The mesh for a CFD simulation should be as aligned as possible to the flow lines to avoid false diffusion. The CFD mesh should also be as orthogonal as possible as there is no skewness correction implemented.

[Mesh]
  coord_type = RZ
  # Mesh should be fairly orthogonal for finite volume fluid flow
  # If you are running this input file for the first time, run core_with_reflectors.py
  # in pbfhr/meshes using Cubit to generate the mesh
  # Modify the parameters (mesh size, refinement areas) for each application
  # neutronics, thermal hydraulics and fuel performance
  uniform_refine = 1
  [fmg]
    type = FileMeshGenerator
    file = 'meshes/core_pronghorn.e'
  []
  [barrel]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = 6
    paired_block = 7
    input = fmg
    new_boundary = 'barrel_wall'
  []
  [OR_inlet]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'abs(y) < 1e-10'
    new_sideset_name = 'OR_horizontal_bottom'
    included_subdomains = '6'
    input = barrel
  []
  [OR_outlet]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'abs(y - 5.3125) < 1e-10'
    new_sideset_name = 'OR_horizontal_top'
    included_subdomains = '6'
    input = OR_inlet
  []
  [add_bed_right]
    type = SideSetsBetweenSubdomainsGenerator
    input = OR_outlet
    new_boundary = 'bed_right'
    primary_block = '6'
    paired_block = '4 5'
  []
  [remove_bad_sideset]
    type = BoundaryDeletionGenerator
    input = add_bed_right
    boundary_names = bed_left
  []
  [add_sideset_again]
    type = SideSetsBetweenSubdomainsGenerator
    input = remove_bad_sideset
    new_boundary = 'bed_left'
    primary_block = '3'
    paired_block = '1 2'
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

The geometry is also specified in the WallDistance model, which is used to account for the presence of the wall in material closures. For the Mk1-FHR we use a WallDistanceAngledCylindricalBed which is specific to the shape of this core. This is defined in the UserObjects block.

[UserObjects]
  [wall_dist]
    type = WallDistanceAngledCylindricalBed
    outer_radius_x = '0.8574 0.8574 1.25 1.25 0.89 0.89'
    outer_radius_y = '0.0 0.709 1.389 3.889 4.5125 5.3125'
    inner_radius_x = '0.45 0.45 0.35 0.35 0.71 0.71'
    inner_radius_y = '0.0 0.859 1.0322 3.889 4.5125 5.3125'
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

Variables

The fluid variables are defined automatically by the NavierStokesFV action: the two velocity components, pressure and fluid temperature. They are defined over the active region and the pebble reflector, the two fluid flow regions. The model may evolve in the future to model flow in the inner and outer reflectors.

The solid temperature is defined in the entire domain, as the solid phase temperature in the pebble bed, and as the regular solid temperature in the rest of the core. It is defined as shown in the block below. The NavierStokesFV syntax may in the future include the definition of the solid variables as porous flow simulation capabilities increase.

[Variables]
  [T_solid]
    type = INSFVEnergyVariable
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

The fluid system includes conservation equations for fluid mass, momentum, and energy. The conservation of energy for the solid / solid phase is solved simultaneously.

commentnote

We could have used another MultiApp setup to compute the solid temperature, with the justification that it evolves on a longer timescale so it takes longer to reach steady state. However, the cost of the more expensive full-core solves and of the additional fluid flow solves, while the solid temperature is converging, is sufficiently offset by not needing to iterate the coupling between two applications. The workflow is also generally simplified.

commentnote

Legacy kernel syntax for the mass/momentum/energy equations, that does not use the NavierStokesFV action, is show here.

Conservation of fluid mass

The conservation of mass is,

with the density of the fluid, the porosity of the homogenized medium and the interstitial velocity. We reformulate this equation in terms of the superficial or Darcy velocity .

Here the system will be simplified by modeling the flow as incompressible. (The effect of buoyancy will be re-introduced later with the Boussinesq approximation.) The simplified conservation of mass is then given by,

This conservation is added by the NavierStokesFV action by simply specifying the compressibility of the fluid and the porous media treatment:


[Modules]
  [NavierStokesFV]
    # General parameters
    compressibility = 'incompressible'
    porous_medium_treatment = true
    ...
  []
[]

Conservation of fluid momentum

This system also includes the conservation of momentum in the -direction. We use a transient simulation to reach steady state, and the porous media Navier Stokes equation can be written as:

In this model, gravity will point in the negative -direction so the quantity is zero for , the x-direction velocity,

The effective viscosity is the Brinkman viscosity. There is insufficient agreement in the literature about whether this term should be considered and its magnitude. The study on which this model is based neglected it (Novak et al., 2021). We simply used the regular fluid viscosity as a placeholder for future refinements.

The drag term represents the interphase drag. It accounts for both viscous and inertial effects. It can be anisotropic in solid components such as the reflector and isotropic in the pebble bed. It's a superposition of a linear and quadratic in velocity terms, which are dominant in low/high velocity regions respectively. A common correlation for the pebble bed is to use an Ergun drag coefficient, detailed in the Pronghorn manual. To model the outer reflector using a porous medium, we would need to tune the anisotropic drag coefficient using CFD simulations.

Finally, we must collect all of the terms on one side of the equation and express the equation in terms of the superficial velocity. This gives the form that is implemented for the FHR model, (1)

This equation is added automatically by specifying the compressibility and the porous media option to the NavierStokesFV action. The flow equations are restricted to the fluid domain using the block parameter, and the fluid energy equation is added similarly.


[Modules]
  [NavierStokesFV]
    # General parameters
    compressibility = 'incompressible'
    porous_medium_treatment = true
    block = ${blocks_fluid}
    ...

The friction terms in the porous media treatment are specified in the action as below:


[Modules]
  [NavierStokesFV]
    ...
    # Friction in porous media
    friction_types = 'darcy forchheimer'
    friction_coeffs = 'Darcy_coefficient Forchheimer_coefficient'

The conservation of momentum in the -direction is analogous, but it also includes the Boussinesq approximation in order to capture the effect of buoyancy. Note that this extra term is needed because of the approximation that the fluid density is uniform and constant,

For each kernel describing the -momentum equation, there is a corresponding kernel for the -momentum equation. The additional Boussinesq and gravity kernels for this equation are added by specifying:


[Modules]
  [NavierStokesFV]
    ...
    boussinesq_approximation = true

Conservation of fluid energy

The conservation of the fluid energy can be expressed as, (2) where is the fluid specific enthalpy, is the fluid effective thermal conductivity, is the fluid temperature, the solid temperature, and is the heat transfer coefficient between the two phases. The heat generation term is featured in the solid phase energy equation, and heat is transfered by convection between the two phases.

We neglect the work terms of gravity, friction and the buoyancy force compared to the heat generation. It is also expected that the energy released from nuclear reactions will be very large compared to pressure work terms. Consequently, we will use the simplified form,

The effective thermal diffusivity represents both the molecular diffusion and thermal dispersion in the fluid. We use a closure with a linear Peclet number dependence valid at low Reynolds number (<800).

The interphase heat transfer coefficient represents convective heat transfer between the solid and fluid phase. We use the Wakao correlation (Wakao et al., 1979), detailed in the Pronghorn manual and commonly used in pebble-bed reactor analysis.

The fluid energy equation is added automatically through the add_energy_equation parameter, and the volumetric heat convection with the solid phase is specified as below:


[NavierStokesFV]
  ...
  add_energy_equation = true
  ...

  # Porous flow parameters
  ambient_convection_blocks = ${blocks_pebbles}
  ambient_convection_alpha = 'alpha alpha'
  ambient_temperature = 'T_solid T_solid'

Conservation of solid energy

The conservation of energy in the solid phase may be written in terms of the temperature as: (3)

with the effective solid thermal conductivity in the pebble bed, and the solid thermal conductivity in the other solid parts and the heat source.

Effective solid diffusion accounts for

  • conduction inside pebbles and radiation between pebbles across a transparent fluid

  • conduction inside pebbles and conduction in the fluid between pebbles

  • conduction inside pebbles and conduction between pebbles at contact areas

Each effect is represented by a separate correlation, lumped together in the effective thermal diffusivity.

The first term of Eq. (3)—the energy time derivative—is captured by the kernel,

[FVKernels]
  [temp_solid_time]
    type = PINSFVEnergyTimeDerivative
    variable = T_solid
    cp = 'cp_s'
    rho = 'rho_s'
    porosity = 0
    is_solid = true
    block = ${blocks_solid}
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

The second term—the effective diffusion and diffusion of heat—corresponds to the kernels,

[FVKernels]
  [temp_solid_conduction]
    type = FVDiffusion
    variable = T_solid
    coeff = 'k_s'
    block = ${blocks_solid}
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

The third term—the fluid-solid heat convection—is added by the kernel,

[FVKernels]
  [temp_fluid_to_solid]
    type = PINSFVEnergyAmbientConvection
    variable = T_solid
    is_solid = true
    h_solid_fluid = 'alpha'
    block = ${blocks_fluid}
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

The final term—the heat source—is added by the kernel,

[FVKernels]
  [temp_solid_source]
    type = FVCoupledForce
    variable = T_solid
    v = power_distribution
    block = '3'
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

Auxiliary system and initialization

We define two auxiliary variables in this simulation: porosity and the power distribution. While porosity is intrinsically more of a material property for the homogenized medium, we may need to compute its gradient on cell faces in the future versions of the models, something which is currently not possible with material properties. Power distribution is an auxiliary field that is passed on from the neutronics simulation, and used to compute the energy conservation in the solid phase.

[AuxVariables]
  [power_distribution]
    type = MooseVariableFVReal
    block = '3'
  []
  [porosity_viz]
    type = MooseVariableFVReal
    block = ${blocks_fluid}
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

There are numerous options to initialize a Variable or an AuxVariable in MOOSE. The initial_condition can be set directly in the relevant Variables block. It can also be set using an initial condition, ICs, block. The velocity variables have to be initialized to a non-zero value, to avoid numerical issues with advection at 0 velocity. Initializations are specified in the NavierStokesFV action:


[Modules]
  [NavierStokesFV]
    # Initial conditions
    initial_velocity = '1e-6 ${inlet_vel_y} 0'
    initial_pressure = 2e5
    initial_temperature = 800.0

We initialize in the [ICs] block the power distribution, which will be overriden by the power distribution provided by Griffin, and the solid temperature. The solid temperature can take a long time to relax to its steady state value, so a closer reasonable flat guess can save computation time.

[ICs]
  [pow_init1]
    type = FunctionIC
    variable = power_distribution
    function = ${power_density}
    block = '3'
  []
  [core]
    type = FunctionIC
    variable = T_solid
    function = 900
    block = '1 2 3 4 5 6 7 8'
  []
  [bricks]
    type = FunctionIC
    variable = T_solid
    function = 350
    block = '9 10'
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

Another important part of the simulation initialization is the viscosity ramp-down. It is very difficult to initialize a high-Reynolds number fluid simulation from an initial guess and let it relax to the steady state simulation. To avoid those numerical difficulties, we initialize the simulation with a very high viscosity, making the flow very slow and easy to solve for. We then ramp-down viscosity to its value from (Lopez et al., 2013) over a few seconds. This is not particularly expensive considering the length of the relaxation pseudo-transient. The ramp down is performed using a piecewise linear function and a functionalized material property.

[Functions]
  [mu_func]
    type = ParsedFunction
    expression = 'if(num_fixed_point>2, 1, 1000*exp(-3*t) + 1)'
    symbol_names = 'num_fixed_point'
    symbol_values = 'num_fixed_point'
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

Boundary conditions

commentnote

Legacy syntax for the boundary conditions of the mass/momentum/energy equations, that does not use the NavierStokesFV action, is shown here.

We first define the inlet of the core. We specify the velocity of the fluid at the inlet and its temperature. In this simplified model of the Mk1-FHR, there is no flow coming from the inner reflector. The velocity and temperature are currently set to a constant value, but will be coupled in the future to a SAM simulation of the primary loop.


[Modules]
  [NavierStokesFV]
    # Inlet boundary conditions
    inlet_boundaries = 'bed_horizontal_bottom OR_horizontal_bottom'
    momentum_inlet_types = 'fixed-velocity fixed-velocity'
    momentum_inlet_function = '0 ${inlet_vel_y} 0 0'
    energy_inlet_types = 'fixed-temperature heatflux'
    energy_inlet_function = '${inlet_T_fluid} 0'
    # so the flux BCs have to be used consistently across all equations

Since the model does not include flow in the inner and outer reflector, we define wall boundary conditions on these surfaces. We chose free-slip boundary conditions as the friction on the fluid is heavily dominated by the friction with the pebbles, so wall friction may be neglected. The heat transfer at the pebble bed wall is assumed to take place in the solid phase only, so it is set to 0 heat flux.


[Modules]
  [NavierStokesFV]
    ...
    # Wall boundary conditions
    wall_boundaries = 'bed_left barrel_wall'
    momentum_wall_types = 'slip slip'
    energy_wall_types = 'heatflux heatflux'
    energy_wall_function = '0 0'

The outflow boundary condition is a pressure boundary condition. Since this model does not include flow in the outer reflector, all the flow goes through the defueling chute. The velocity is very high in this region, causing a large pressure drop. This is a known result of the simplified model. This boundary condition is a fully developed flow boundary condition. It may only be used sufficiently far from modifications of the flow path.


[Modules]
  [NavierStokesFV]
    ...
    # Outlet boundary conditions
    outlet_boundaries = 'bed_horizontal_top plenum_top OR_horizontal_top'
    momentum_outlet_types = 'fixed-pressure fixed-pressure fixed-pressure'
    pressure_function = '2e5 2e5 2e5'

Finally, we define the boundary conditions for the solid temperature. They are implicitly defined to be zero flux boundary conditions on the center, top and bottom surfaces. For the center, this represents the axi-symmetry of the model. For the top and bottom surfaces, this is equivalent to perfect insulation. Future evolutions of the model will distinguish between the pebble bed solid phase and top and bottom surfaces of the core regions. Finally, on the outer cylindrical boundary of the model, we define a fixed temperature boundary condition, as an approximation of the outer environment.

[FVBCs]
  [outer]
    type = FVDirichletBC
    variable = T_solid
    boundary = 'brick_surface'
    value = '${fparse 35 + 273.15}'
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

Material properties

There are four main types of materials that may be found in Pronghorn input files:

  • solid material properties

  • closure relations (, , )

  • fluid properties

Solid material properties

The properties of the firebrick around the core are specified directly, using an ADGenericConstantMaterial, neglecting their temperature dependence.

[FunctorMaterials]
  [firebrick_properties]
    type = ADGenericFunctorMaterial
    prop_names = 'rho_s         cp_s        k_s'
    prop_values = '${solid_rho} ${solid_cp} 0.26'
    block = '10'
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

The properties of the other solid phases are obtained through volumetric mixing. The base material properties of each type of graphite, UO2 and other materials are defined first:

[UserObjects]
  [graphite]
    type = FunctionSolidProperties
    rho_s = 1780
    cp_s = '${fparse 1800.0 * heat_capacity_multiplier}'
    k_s = 26.0
  []

  [pebble_graphite]
    type = FunctionSolidProperties
    rho_s = 1600.0
    cp_s = 1800.0
    k_s = 15.0
  []

  [pebble_core]
    type = FunctionSolidProperties
    rho_s = 1450.0
    cp_s = 1800.0
    k_s = 15.0
  []

  [UO2]
    type = FunctionSolidProperties
    rho_s = 11000.0
    cp_s = 400.0
    k_s = 3.5
  []

  [pyc]
    type = PyroliticGraphite # (constant)
  []

  [buffer]
    type = PorousGraphite # (constant)
  []

  [SiC]
    type = FunctionSolidProperties
    rho_s = 3180.0
    cp_s = 1300.0
    k_s = 13.9
  []

  [solid_flibe]
    type = FunctionSolidProperties
    rho_s = 1986.62668
    cp_s = 2416.0
    k_s = 1.0665
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

Then these properties are mixed using CompositeSolidProperties. We specify the volumetric fraction of each material and how the properties should be combined. The default option is simply volume averaging. The Chiew correlation is used for mixing thermal conductivities in the fuel matrix in the pebble. Note that this step is obtaining macroscopic phase solid properties, the temperature inside the pellet is resolved later on using a MultiApp. The effective thermal conductivity is determined by a closure relation, further detailed below.

[UserObjects]
  [TRISO]
    type = CompositeSolidProperties
    materials = 'UO2 buffer pyc SiC pyc'
    fractions = '${UO2_phase_fraction} ${buffer_phase_fraction} ${ipyc_phase_fraction} ${sic_phase_fraction} ${opyc_phase_fraction}'
  []

  [fuel_matrix]
    type = CompositeSolidProperties
    materials = 'TRISO pebble_graphite'
    fractions = '${TRISO_phase_fraction} ${fparse 1.0 - TRISO_phase_fraction}'
    k_mixing = 'chiew'
  []

  [pebble]
    type = CompositeSolidProperties
    materials = 'pebble_core fuel_matrix pebble_graphite'
    fractions = '${core_phase_fraction} ${fuel_matrix_phase_fraction} ${shell_phase_fraction}'
  []

  [inner_reflector]
    type = CompositeSolidProperties
    materials = 'solid_flibe graphite'
    fractions = '${IR_porosity} ${fparse 1.0 - IR_porosity}'
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

Finally, these materials, defined as user objects, are placed in each subdomain of the mesh in the Materials block.

[FunctorMaterials]
  [solid_fuel_pebbles]
    type = PronghornSolidFunctorMaterialPT
    solid = pebble
    block = '3'
  []

  [solid_blanket_pebbles]
    type = PronghornSolidFunctorMaterialPT
    solid = graphite
    block = '4'
  []

  [plenum_and_OR]
    type = PronghornSolidFunctorMaterialPT
    solid = graphite
    block = '5 6 8'
  []

  [IR]
    type = PronghornSolidFunctorMaterialPT
    solid = inner_reflector
    block = '1 2'
  []

  [barrel_and_vessel]
    type = PronghornSolidFunctorMaterialPT
    solid = stainless_steel
    block = '7 9'
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

Closure relations

The closure relations used for friction and effective thermal diffusivities are documented in the Pronghorn manual.

[FunctorMaterials]
  [alpha]
    type = FunctorWakaoPebbleBedHTC
    block = ${blocks_pebbles}
  []

  [drag]
    type = FunctorKTADragCoefficients
    block = ${blocks_pebbles}
  []

  [kappa]
    type = FunctorLinearPecletKappaFluid
    block = ${blocks_pebbles}
  []

  [kappa_s]
    type = FunctorPebbleBedKappaSolid
    emissivity = 0.8
    Youngs_modulus = 9e9
    Poisson_ratio = 0.136
    wall_distance = wall_dist
    block = ${blocks_pebbles}
    T_solid = T_solid
    acceleration = '0 -9.81 0'
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

Fluid properties

The fluid properties are similar to variable materials but are specific to Pronghorn. They act as an aggregator of both materials properties and variables and allow the computation of closure inputs such as the Reynolds and Prandtl number.

[FluidProperties]
  [fp]
    type = FlibeFluidProperties
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)
[FunctorMaterials]
  [fluidprops]
    type = GeneralFunctorFluidProps
    block = ${blocks_fluid}
    mu_rampdown = mu_func
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

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. The petsc_options_iname/ivalue specify a reasonably scaling pre-conditioner. The factor shift avoids numerical issues when there are zeros in the diagonal of the Jacobian, for example when using Lagrange multipliers to enforce constraints, which some iterations of this model used. MOOSE defaults to GMRES for the linear solver, and we increase the size of the base used as using too small a Krylov vector base is a common issue.

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 since the bricks around the core are very slow to heat up. We use an adaptive time-stepper to increase the time step over the course of the simulation. Because we are using an implicit scheme, we can take very large time steps and still have a stable simulation. An alternative strategy is to artificially decrease the specific heat capacity of the solid materials for the steady state solve.

[Executioner]
  type = Transient

  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -ksp_gmres_restart'
  petsc_options_value = 'asm      lu           NONZERO                   200'
  line_search = 'none'

  # Iterations parameters
  l_max_its = 500
  l_tol = 1e-8
  nl_max_its = 25
  nl_rel_tol = 5e-7
  nl_abs_tol = 5e-7

  # Automatic scaling
  automatic_scaling = true

  # Problem time parameters
  dtmin = 0.1
  dtmax = 2e4
  end_time = 1e6

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.15
    cutback_factor = 0.5
    growth_factor = 2.0
  []

  # Time integration scheme
  scheme = 'implicit-euler'

  # Steady state detection.
  steady_state_detection = true
  steady_state_tolerance = 1e-8
  steady_state_start_time = 200000
[]
(pbfhr/mark1/steady/ss1_combined.i)

MultiApp Coupling

We couple the thermal-hydraulics simulation with simulation of the pebbles spread throughout the core. The solid phase temperature is then the boundary condition to a heat conduction problem that allows us to examine the temperature profile in fueled and reflector pebbles. The sub-application is defined using the MultiApps block.

[MultiApps]
  [coarse_mesh]
    type = TransientMultiApp
    execute_on = 'TIMESTEP_END'
    input_files = 'ss3_coarse_pebble_mesh.i'
    cli_args = 'Outputs/console=false'
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

Pronghorn is running as a sub-application of Griffin. As such it does not need to handle transfers with Griffin, as they are defined in Griffin's input file. We transfer the power distribution, previously received from Griffin, and the solid phase temperature in the pebble bed.

[Transfers]
  [fuel_matrix_heat_source]
    type = MultiAppProjectionTransfer
    to_multi_app = coarse_mesh
    source_variable = power_distribution
    variable = power_distribution
  []
  [pebble_surface_temp]
    type = MultiAppProjectionTransfer
    to_multi_app = coarse_mesh
    source_variable = T_solid
    variable = temp_solid
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

Outputs

We first define a few postprocessors to:

  • help quantify the convergence of our mesh.

[Postprocessors]
  [max_Tf]
    type = ElementExtremeValue
    variable = T_fluid
    block = ${blocks_fluid}
  []

  [max_vy]
    type = ElementExtremeValue
    variable = superficial_vel_y
    block = ${blocks_fluid}
  []

  [pressure_in]
    type = SideAverageValue
    boundary = 'bed_horizontal_bottom'
    variable = pressure
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

  • pass inlet and outlet conditions to a future SAM model of the primary loop

[Postprocessors]
  [pressure_in]
    type = SideAverageValue
    boundary = 'bed_horizontal_bottom'
    variable = pressure
    execute_on = 'INITIAL TIMESTEP_END'
  []

  [T_flow_out]
    type = SideAverageValue
    boundary = 'bed_horizontal_top plenum_top OR_horizontal_top'
    variable = T_fluid
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

  • examine conservation of mass

[Postprocessors]
  [mass_flow_out]
    type = VolumetricFlowRate
    boundary = 'bed_horizontal_top plenum_top OR_horizontal_top'
    advected_variable = ${rho_fluid}
    advected_quantity = ${rho_fluid}
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]
(pbfhr/mark1/steady/ss1_combined.i)

  • examine conservation of energy (missing diffusion term, under development)

# ==============================================================================
# Model description
# Application : Pronghorn
# ------------------------------------------------------------------------------
# Idaho Falls, INL, August 10, 2020
# Author(s): Dr. Guillaume Giudicelli, Dr. Paolo Balestra, Dr. April Novak
# If using or referring to this model, please cite as explained on
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
# ==============================================================================
# - Coupled fluid-solid thermal hydraulics model of the Mk1-FHR
# ==============================================================================
# - The Model has been built based on [1-2].
# ------------------------------------------------------------------------------
# [1] Multiscale Core Thermal Hydraulics Analysis of the Pebble Bed Fluoride
#     Salt Cooled High Temperature Reactor (PB-FHR), A. Novak et al.
# [2] Technical Description of the “Mark 1” Pebble-Bed Fluoride-Salt-Cooled
#     High-Temperature Reactor (PB-FHR) Power Plant, UC Berkeley report 14-002
# [3] Molten salts database for energy applications, Serrano-Lopez et al.
#     https://arxiv.org/pdf/1307.7343.pdf
# [4] Heat Transfer Salts for Nuclear Reactor Systems - chemistry control,
#     corrosion mitigation and modeling, CFP-10-100, Anderson et al.
#     https://neup.inl.gov/SiteAssets/Final%20%20Reports/FY%202010/
#     10-905%20NEUP%20Final%20Report.pdf
# ==============================================================================
# MODEL PARAMETERS
# ==============================================================================
# Problem Parameters -----------------------------------------------------------

blocks_pebbles = '3 4'
blocks_fluid = '3 4 5 6'
blocks_solid = '1 2 6 7 8 9 10'

# Material compositions
UO2_phase_fraction = 1.20427291e-01
buffer_phase_fraction = 2.86014816e-01
ipyc_phase_fraction = 1.59496539e-01
sic_phase_fraction = 1.96561801e-01
opyc_phase_fraction = 2.37499553e-01
TRISO_phase_fraction = 3.09266232e-01
core_phase_fraction = 5.12000000e-01
fuel_matrix_phase_fraction = 3.01037037e-01
shell_phase_fraction = 1.86962963e-01

# FLiBe properties #TODO: Rely only on PronghornFluidProps
fluid_mu = 7.5e-3 # Pa.s at 900K [3]
# k_fluid =  1.1     # suggested constant [3]
# cp_fluid = 2385    # suggested constant [3]
rho_fluid = 1970.0 # kg/m3 at 900K [3]
alpha_b = 2e-4 # /K from [4]

# Graphite properties
heat_capacity_multiplier = 1e0 # >1 gets faster to steady state
solid_rho = 1780.0
# solid_k = 26.0
solid_cp = '${fparse 1697.0*heat_capacity_multiplier}'

# Core geometry
pebble_diameter = 0.03
bed_porosity = 0.4
IR_porosity = 0
OR_porosity = 0.1123
plenum_porosity = 0.5
model_inlet_rin = 0.45
model_inlet_rout = 0.8574
model_vol = 10.4
model_inlet_area = '${fparse 3.14159265 * (model_inlet_rout * model_inlet_rout - model_inlet_rin * model_inlet_rin)}'

# The convention for friction factors changed
darcy_conversion_plenum = '${fparse rho_fluid / plenum_porosity / fluid_mu}'
# The porosity simplfies out from the previous definition of the coefficients
darcy_conversion_or = '${fparse rho_fluid / fluid_mu}'
forchheimer_partial_conversion = 2 # divide further by speed

# Outer reflector drag parameters
# TODO: tune using CFD
# TODO: verify current values (imported from [1])
Ah = '${fparse 1337.76 * darcy_conversion_or}'
Bh = '${fparse 2.58 * forchheimer_partial_conversion}'
Av = '${fparse 599.30 * darcy_conversion_or}'
Bv = '${fparse 0.95 * forchheimer_partial_conversion}'
Dh = 0.02582
Dv = 0.02006

# Plenum drag parameters. The core hot legs cannot be represented in 2D RZ
# accurately, so this should be tuned to obtain the desired mass flow rates
plenum_friction = '${fparse 0.5 * darcy_conversion_plenum}'

# Operating parameters
mfr = 976.0 # kg/s, from [2]
total_power = 236.0e6 # W, from [2]
inlet_T_fluid = 873.15 # K, from [2]
inlet_vel_y = '${fparse mfr / model_inlet_area / rho_fluid}' # superficial
power_density = '${fparse total_power / model_vol / 258 * 236}' # adjusted using power pp

# ==============================================================================
# GEOMETRY AND MESH
# ==============================================================================

[Mesh]
  coord_type = RZ
  # Mesh should be fairly orthogonal for finite volume fluid flow
  # If you are running this input file for the first time, run core_with_reflectors.py
  # in pbfhr/meshes using Cubit to generate the mesh
  # Modify the parameters (mesh size, refinement areas) for each application
  # neutronics, thermal hydraulics and fuel performance
  uniform_refine = 1
  [fmg]
    type = FileMeshGenerator
    file = 'meshes/core_pronghorn.e'
  []
  [barrel]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = 6
    paired_block = 7
    input = fmg
    new_boundary = 'barrel_wall'
  []
  [OR_inlet]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'abs(y) < 1e-10'
    new_sideset_name = 'OR_horizontal_bottom'
    included_subdomains = '6'
    input = barrel
  []
  [OR_outlet]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'abs(y - 5.3125) < 1e-10'
    new_sideset_name = 'OR_horizontal_top'
    included_subdomains = '6'
    input = OR_inlet
  []
  [add_bed_right]
    type = SideSetsBetweenSubdomainsGenerator
    input = OR_outlet
    new_boundary = 'bed_right'
    primary_block = '6'
    paired_block = '4 5'
  []
  [remove_bad_sideset]
    type = BoundaryDeletionGenerator
    input = add_bed_right
    boundary_names = bed_left
  []
  [add_sideset_again]
    type = SideSetsBetweenSubdomainsGenerator
    input = remove_bad_sideset
    new_boundary = 'bed_left'
    primary_block = '3'
    paired_block = '1 2'
  []
[]

[GlobalParams]
  rho = ${rho_fluid}
  porosity = porosity_viz
  characteristic_length = 0.01 #${pebble_diameter}
  pebble_diameter = ${pebble_diameter}
  speed = 'speed'

  fp = fp
  T_solid = T_solid
  T_fluid = T_fluid
  pressure = pressure

  vel_x = 'superficial_vel_x'
  vel_y = 'superficial_vel_y'

  fv = true
  rhie_chow_user_object = 'pins_rhie_chow_interpolator'
[]

# ==============================================================================
# VARIABLES AND KERNELS
# ==============================================================================

[Problem]
  kernel_coverage_check = false
[]

[Modules]
  [NavierStokesFV]
    # General parameters
    compressibility = 'incompressible'
    porous_medium_treatment = true
    add_energy_equation = true
    boussinesq_approximation = true
    block = ${blocks_fluid}

    # Material properties
    # density should be explicitly defined as constant
    density = ${rho_fluid}
    dynamic_viscosity = 'mu'
    thermal_conductivity = '0' #kappa'
    specific_heat = 'cp'
    thermal_expansion = 'alpha_b'
    porosity = 'porosity'

    # Boussinesq parameters
    gravity = '0 -9.81 0'
    ref_temperature = ${inlet_T_fluid}

    # Initial conditions
    initial_velocity = '1e-6 ${inlet_vel_y} 0'
    initial_pressure = 2e5
    initial_temperature = 800.0

    # Wall boundary conditions
    wall_boundaries = 'bed_left barrel_wall'
    momentum_wall_types = 'slip slip'
    energy_wall_types = 'heatflux heatflux'
    energy_wall_function = '0 0'

    # Inlet boundary conditions
    inlet_boundaries = 'bed_horizontal_bottom OR_horizontal_bottom'
    momentum_inlet_types = 'fixed-velocity fixed-velocity'
    momentum_inlet_function = '0 ${inlet_vel_y}; 0 0'
    energy_inlet_types = 'fixed-temperature heatflux'
    energy_inlet_function = '${inlet_T_fluid} 0'
    # so the flux BCs have to be used consistently across all equations

    # Outlet boundary conditions
    outlet_boundaries = 'bed_horizontal_top plenum_top OR_horizontal_top'
    momentum_outlet_types = 'fixed-pressure fixed-pressure fixed-pressure'
    pressure_function = '2e5 2e5 2e5'

    # Porous flow parameters
    ambient_convection_blocks = ${blocks_pebbles}
    ambient_convection_alpha = 'alpha'
    ambient_temperature = 'T_solid'

    # Friction in porous media
    friction_types = 'darcy forchheimer'
    friction_coeffs = 'Darcy_coefficient Forchheimer_coefficient'

    # Numerical scheme
    momentum_advection_interpolation = 'upwind'
    mass_advection_interpolation = 'upwind'
    energy_advection_interpolation = 'upwind'
  []
[]

[Variables]
  [T_solid]
    type = INSFVEnergyVariable
  []
[]

[FVKernels]
  # Fluid heat diffusion
  [temp_fluid_conduction]
    type = PINSFVEnergyAnisotropicDiffusion
    variable = T_fluid
    effective_diffusivity = false
    kappa = 'kappa'
  []

  # Solid Energy equation.
  [temp_solid_time_core]
    type = PINSFVEnergyTimeDerivative
    variable = T_solid
    cp = 'cp_s'
    rho = ${solid_rho}
    is_solid = true
    block = ${blocks_fluid}
  []
  [temp_solid_time]
    type = PINSFVEnergyTimeDerivative
    variable = T_solid
    cp = 'cp_s'
    rho = 'rho_s'
    porosity = 0
    is_solid = true
    block = ${blocks_solid}
  []
  [temp_solid_conduction_core]
    type = FVDiffusion
    variable = T_solid
    coeff = 'kappa_s'
    block = ${blocks_fluid}
  []
  [temp_solid_conduction]
    type = FVDiffusion
    variable = T_solid
    coeff = 'k_s'
    block = ${blocks_solid}
  []
  [temp_solid_source]
    type = FVCoupledForce
    variable = T_solid
    v = power_distribution
    block = '3'
  []
  [temp_fluid_to_solid]
    type = PINSFVEnergyAmbientConvection
    variable = T_solid
    is_solid = true
    h_solid_fluid = 'alpha'
    block = ${blocks_fluid}
  []
[]

[FVInterfaceKernels]
  [diffusion_interface]
    type = FVDiffusionInterface
    boundary = 'bed_left'
    subdomain1 = '3 4 5'
    subdomain2 = '1 2 6'
    coeff1 = 'kappa_s'
    coeff2 = 'k_s'
    variable1 = 'T_solid'
    variable2 = 'T_solid'
  []
[]

# ==============================================================================
# AUXVARIABLES AND AUXKERNELS
# ==============================================================================
[AuxVariables]
  [power_distribution]
    type = MooseVariableFVReal
    block = '3'
  []
  [porosity_viz]
    type = MooseVariableFVReal
    block = ${blocks_fluid}
  []
[]

[AuxKernels]
  [eps]
    type = FunctorAux
    variable = porosity_viz
    functor = porosity
    block = ${blocks_fluid}
    execute_on = 'INITIAL'
  []
[]

# ==============================================================================
# INITIAL CONDITIONS AND FUNCTIONS
# ==============================================================================
[ICs]
  [pow_init1]
    type = FunctionIC
    variable = power_distribution
    function = ${power_density}
    block = '3'
  []
  [core]
    type = FunctionIC
    variable = T_solid
    function = 900
    block = '1 2 3 4 5 6 7 8'
  []
  [bricks]
    type = FunctionIC
    variable = T_solid
    function = 350
    block = '9 10'
  []
[]

[Functions]
  # Only perform a viscosity rampdown for the first coupling iteration
  [mu_func]
    type = ParsedFunction
    expression = 'if(num_fixed_point>2, 1, 1000*exp(-3*t) + 1)'
    symbol_names = 'num_fixed_point'
    symbol_values = 'num_fixed_point'
  []
[]

# ==============================================================================
# BOUNDARY CONDITIONS
# ==============================================================================
[FVBCs]
  [outer]
    type = FVDirichletBC
    variable = T_solid
    boundary = 'brick_surface'
    value = '${fparse 35 + 273.15}'
  []
[]

# ==============================================================================
# FLUID PROPERTIES, MATERIALS AND USER OBJECTS
# ==============================================================================
[FluidProperties]
  [fp]
    type = FlibeFluidProperties
  []
[]

[FunctorMaterials]
  # solid material properties
  [solid_fuel_pebbles]
    type = PronghornSolidFunctorMaterialPT
    solid = pebble
    block = '3'
  []
  [solid_blanket_pebbles]
    type = PronghornSolidFunctorMaterialPT
    solid = graphite
    block = '4'
  []
  [plenum_and_OR]
    type = PronghornSolidFunctorMaterialPT
    solid = graphite
    block = '5 6 8'
  []
  [IR]
    type = PronghornSolidFunctorMaterialPT
    solid = inner_reflector
    block = '1 2'
  []
  [barrel_and_vessel]
    type = PronghornSolidFunctorMaterialPT
    solid = stainless_steel
    block = '7 9'
  []
  [firebrick_properties]
    type = ADGenericFunctorMaterial
    prop_names = 'rho_s         cp_s        k_s'
    prop_values = '${solid_rho} ${solid_cp} 0.26'
    block = '10'
  []

  # FLUID
  [alpha_boussinesq]
    type = ADGenericFunctorMaterial
    prop_names = 'alpha_b rho'
    prop_values = '${alpha_b} ${rho_fluid}'
    block = ${blocks_fluid}
    define_dot_functors = false
  []
  [fluidprops]
    type = GeneralFunctorFluidProps
    block = ${blocks_fluid}
    mu_rampdown = mu_func
  []

  # closures in the pebble bed
  [alpha]
    type = FunctorWakaoPebbleBedHTC
    block = ${blocks_pebbles}
  []
  [drag]
    type = FunctorKTADragCoefficients
    block = ${blocks_pebbles}
  []
  [kappa]
    type = FunctorLinearPecletKappaFluid
    block = ${blocks_pebbles}
  []
  [kappa_s]
    type = FunctorPebbleBedKappaSolid
    emissivity = 0.8
    Youngs_modulus = 9e9
    Poisson_ratio = 0.136
    wall_distance = wall_dist
    block = ${blocks_pebbles}
    T_solid = T_solid
    acceleration = '0 -9.81 0'
  []

  # closures in the outer reflector and the plenum
  [Fh]
    type = ADParsedFunctorMaterial
    property_name = Fh
    expression = 'Bh / Dh / speed'
    functor_symbols = 'Bh Dh speed'
    functor_names = '${Bh} ${Dh} speed'
  []
  [Fv]
    type = ADParsedFunctorMaterial
    property_name = Fv
    expression = 'Bv / Dv / speed'
    functor_symbols = 'Bv Dv speed'
    functor_names = '${Bv} ${Dv} speed'
  []
  [drag_OR]
    type = FunctorAnisotropicFunctorDragCoefficients
    Darcy_coefficient = '${fparse Ah / Dh / Dh} ${fparse Av / Dv / Dv} ${fparse Av / Dv / Dv}'
    Forchheimer_coefficient = 'Fh Fv Fv'
    block = '6'
  []
  [drag_plenum]
    type = ADGenericVectorFunctorMaterial
    prop_names = 'Darcy_coefficient Forchheimer_coefficient'
    prop_values = '${plenum_friction} ${plenum_friction} ${plenum_friction}
                   0 0 0'
    block = '5'
  []
  [alpha_OR_plenum]
    type = ADGenericFunctorMaterial
    prop_names = 'alpha'
    prop_values = '0.0'
    block = '5 6'
  []
  [kappa_OR_plenum]
    type = FunctorKappaFluid
    block = '5 6'
  []
  [kappa_s_OR_plenum]
    type = FunctorVolumeAverageKappaSolid
    block = '5 6'
  []

  # porosity
  [porosity]
    type = ADPiecewiseByBlockFunctorMaterial
    prop_name = 'porosity'
    subdomain_to_prop_value = '3 ${bed_porosity}
                               4 ${bed_porosity}
                               5 ${plenum_porosity}
                               6 ${OR_porosity}
                               7 ${OR_porosity}' # !!!
  []
[]

[UserObjects]
  [graphite]
    type = FunctionSolidProperties
    rho_s = 1780
    cp_s = '${fparse 1800.0 * heat_capacity_multiplier}'
    k_s = 26.0
  []
  [pebble_graphite]
    type = FunctionSolidProperties
    rho_s = 1600.0
    cp_s = 1800.0
    k_s = 15.0
  []
  [pebble_core]
    type = FunctionSolidProperties
    rho_s = 1450.0
    cp_s = 1800.0
    k_s = 15.0
  []
  [UO2]
    type = FunctionSolidProperties
    rho_s = 11000.0
    cp_s = 400.0
    k_s = 3.5
  []
  [pyc]
    type = PyroliticGraphite # (constant)
  []
  [buffer]
    type = PorousGraphite # (constant)
  []
  [SiC]
    type = FunctionSolidProperties
    rho_s = 3180.0
    cp_s = 1300.0
    k_s = 13.9
  []
  [TRISO]
    type = CompositeSolidProperties
    materials = 'UO2 buffer pyc SiC pyc'
    fractions = '${UO2_phase_fraction} ${buffer_phase_fraction} ${ipyc_phase_fraction} ${sic_phase_fraction} ${opyc_phase_fraction}'
  []
  [fuel_matrix]
    type = CompositeSolidProperties
    materials = 'TRISO pebble_graphite'
    fractions = '${TRISO_phase_fraction} ${fparse 1.0 - TRISO_phase_fraction}'
    k_mixing = 'chiew'
  []
  [pebble]
    type = CompositeSolidProperties
    materials = 'pebble_core fuel_matrix pebble_graphite'
    fractions = '${core_phase_fraction} ${fuel_matrix_phase_fraction} ${shell_phase_fraction}'
  []
  [stainless_steel]
    type = StainlessSteel
  []
  [solid_flibe]
    type = FunctionSolidProperties
    rho_s = 1986.62668
    cp_s = 2416.0
    k_s = 1.0665
  []
  [inner_reflector]
    type = CompositeSolidProperties
    materials = 'solid_flibe graphite'
    fractions = '${IR_porosity} ${fparse 1.0 - IR_porosity}'
  []
  [wall_dist]
    type = WallDistanceAngledCylindricalBed
    outer_radius_x = '0.8574 0.8574 1.25 1.25 0.89 0.89'
    outer_radius_y = '0.0 0.709 1.389 3.889 4.5125 5.3125'
    inner_radius_x = '0.45 0.45 0.35 0.35 0.71 0.71'
    inner_radius_y = '0.0 0.859 1.0322 3.889 4.5125 5.3125'
  []
[]

# ==============================================================================
# EXECUTION PARAMETERS
# ==============================================================================
[Executioner]
  type = Transient

  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -ksp_gmres_restart'
  petsc_options_value = 'asm      lu           NONZERO                   200'
  line_search = 'none'

  # Iterations parameters
  l_max_its = 500
  l_tol = 1e-8
  nl_max_its = 25
  nl_rel_tol = 5e-7
  nl_abs_tol = 5e-7

  # Automatic scaling
  automatic_scaling = true

  # Problem time parameters
  dtmin = 0.1
  dtmax = 2e4
  end_time = 1e6

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.15
    cutback_factor = 0.5
    growth_factor = 2.0
  []

  # Time integration scheme
  scheme = 'implicit-euler'

  # Steady state detection.
  steady_state_detection = true
  steady_state_tolerance = 1e-8
  steady_state_start_time = 200000
[]

# ==============================================================================
# MULTIAPPS FOR PEBBLE MODEL
# ==============================================================================
[MultiApps]
  [coarse_mesh]
    type = TransientMultiApp
    execute_on = 'TIMESTEP_END'
    input_files = 'ss3_coarse_pebble_mesh.i'
    cli_args = 'Outputs/console=false'
  []
[]

[Transfers]
  [fuel_matrix_heat_source]
    type = MultiAppProjectionTransfer
    to_multi_app = coarse_mesh
    source_variable = power_distribution
    variable = power_distribution
  []
  [pebble_surface_temp]
    type = MultiAppProjectionTransfer
    to_multi_app = coarse_mesh
    source_variable = T_solid
    variable = temp_solid
  []
[]

# ==============================================================================
# POSTPROCESSORS DEBUG AND OUTPUTS
# ==============================================================================
[Postprocessors]
  # For future SAM coupling
  # [inlet_vel_y]
  #   type = Receiver
  #   default = ${inlet_vel_y}
  # []
  [inlet_temp_fluid]
    type = Receiver
    default = ${inlet_T_fluid}
  []
  [inlet_mdot]
    type = Receiver
    default = ${mfr}
  []
  # [outlet_pressure]
  #   type = Receiver
  #   default = ${outlet_pressure}
  # []
  [max_Tf]
    type = ElementExtremeValue
    variable = T_fluid
    block = ${blocks_fluid}
  []
  [max_vy]
    type = ElementExtremeValue
    variable = superficial_vel_y
    block = ${blocks_fluid}
  []
  [power]
    type = ElementIntegralVariablePostprocessor
    variable = power_distribution
    block = '3'
    execute_on = 'INITIAL TIMESTEP_BEGIN TRANSFER TIMESTEP_END'
  []
  [mass_flow_out]
    type = VolumetricFlowRate
    boundary = 'bed_horizontal_top plenum_top OR_horizontal_top'
    advected_variable = ${rho_fluid}
    advected_quantity = ${rho_fluid}
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [T_flow_out]
    type = SideAverageValue
    boundary = 'bed_horizontal_top plenum_top OR_horizontal_top'
    variable = T_fluid
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [pressure_in]
    type = SideAverageValue
    boundary = 'bed_horizontal_bottom'
    variable = pressure
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [pressure_out]
    type = SideAverageValue
    boundary = 'bed_horizontal_top plenum_top OR_horizontal_top'
    variable = pressure
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [pressure_drop]
    type = DifferencePostprocessor
    value1 = pressure_out
    value2 = pressure_in
  []

  # Energy balance
  # Energy balance will be shown once #18119 #18123 are merged in MOOSE
  # [outer_heat_loss]
  #   type = ADSideDiffusiveFluxIntegral
  #   boundary = 'brick_surface'
  #   variable = T_solid
  #   diffusivity = 'k_s'
  #   execute_on = 'INITIAL TIMESTEP_END'
  # []
  [flow_in_m]
    type = VolumetricFlowRate
    boundary = 'bed_horizontal_bottom OR_horizontal_bottom'
    advected_quantity = 'rho_cp_temp'
  []
  # [diffusion_in]
  #   type = ADSideVectorDiffusivityFluxIntegral
  #   variable = T_fluid
  #   boundary = 'bed_horizontal_bottom OR_horizontal_bottom'
  #   diffusivity = 'kappa'
  # []
  # diffusion at the top is 0 because of the fully developped flow assumption
  [flow_out]
    type = VolumetricFlowRate
    boundary = 'bed_horizontal_top plenum_top OR_horizontal_top'
    advected_quantity = 'rho_cp_temp'
  []
  [core_balance]
    type = ParsedPostprocessor
    pp_names = 'power flow_in_m flow_out' #diffusion_in  outer_heat_loss'
    expression = 'power - flow_in_m - flow_out' # + diffusion_in + outer_heat_loss'
  []

  # Bypass
  [mass_flow_OR]
    type = VolumetricFlowRate
    boundary = 'OR_horizontal_top'
    advected_quantity = ${rho_fluid}
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [mass_flow_plenum]
    type = VolumetricFlowRate
    boundary = 'plenum_top'
    advected_quantity = ${rho_fluid}
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [bypass_fraction]
    type = ParsedPostprocessor
    pp_names = 'mass_flow_OR mass_flow_out'
    expression = 'mass_flow_OR / mass_flow_out'
  []
  [plenum_fraction]
    type = ParsedPostprocessor
    pp_names = 'mass_flow_plenum mass_flow_out'
    expression = 'mass_flow_plenum / mass_flow_out'
  []

  # Miscellaneous
  [h]
    type = AverageElementSize
    outputs = 'console csv'
    execute_on = 'timestep_end'
  []
  [mu_factor]
    type = FunctionValuePostprocessor
    function = 'mu_func'
  []
  [num_fixed_point]
    type = Receiver
    default = 0
  []
[]

[Outputs]
  csv = true
  [console]
    type = Console
    hide = 'pressure_in pressure_out mass_flow_OR mass_flow_plenum max_vy flow_in_m flow_out h
            num_fixed_point'
  []
  [exodus]
    type = Exodus
  []
  [checkpoint]
    type = Checkpoint
    num_files = 2
    execute_on = 'FINAL'
  []
  # Reduce base output
  print_linear_converged_reason = false
  print_linear_residuals = false
  print_nonlinear_converged_reason = false
[]
(pbfhr/mark1/steady/ss1_combined.i)

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.

[Outputs]
  csv = true
  [console]
    type = Console
    hide = 'pressure_in pressure_out mass_flow_OR mass_flow_plenum max_vy flow_in_m flow_out h
            num_fixed_point'
  []
  [exodus]
    type = Exodus
  []
  [checkpoint]
    type = Checkpoint
    num_files = 2
    execute_on = 'FINAL'
  []
  # Reduce base output
  print_linear_converged_reason = false
  print_linear_residuals = false
  print_nonlinear_converged_reason = false
[]
(pbfhr/mark1/steady/ss1_combined.i)

References

  1. Serrano Lopez, J. Fradera, and S. Cuesta-Lopez. Molten salts database for energy applications. https://arxiv.org/pdf/1307.7343.pdf, 2013.[BibTeX]
  2. A.J. Novak, S. Schunert, R.W. Carlsen, P. Balestra, R.N. Slaybaugh, and R.C. Martineau. Multiscale thermal-hydraulic modeling of the pebble bed fluoride-salt-cooled high-temperature reactor. Annals of Nuclear Energy, 2021. doi:10.1016/j.anucene.2020.107968.[BibTeX]
  3. N. Wakao, S. Kaguei, and T. Funazkri. Effect of Fluid Dispersion Coefficients on Particle-to-Fluid Heat Transfer Coefficients in Packed Beds: Correlation of Nusselt Number. Chemical Engineering Science, 34:325–336, 1979. doi:10.1016/0009-2509(79)85064-2.[BibTeX]