Kernel syntax (legacy) for the 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 expressed by the FVKernels/mass kernel:

[FVKernels]
  [mass]
    type = PINSFVMassAdvection
    variable = pressure
  []
[]

(Note that this kernel uses variable = pressure even though pressure does not appear in the mass conservation equation. This is an artifact of the way systems of equations are described in MOOSE. Since the number of equations must correspond to the number of non-linear variables, we are essentially using the variables to "name" each equation. The momentum equations will be named by the corresponding velocity variable. That just leaves the pressure variable for naming the mass equation, even though pressure does not appear in the mass equation.)

Kernel syntax (legacy) for the 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)

The first term in Eq. (1)—the momentum time derivative—is input with the time derivative kernel:

[FVKernels]
  [vel_x_time]
    type = PINSFVMomentumTimeDerivative
    variable = vel_x
    momentum_component = x
  []
[]

The second term—the advection of momentum—is handled by a PINSFVMomentumAdvection kernel

[FVKernels]
  [vel_x_advection]
    type = PINSFVMomentumAdvection
    variable = vel_x
    momentum_component = x
  []
[]

The third term—the pressure gradient—is handled by a PINSFVMomentumPressure kernel,

[FVKernels]
  [u_pressure]
    type = PINSFVMomentumPressure
    variable = vel_x
    momentum_component = 'x'
  []
[]

The third term—the effective diffusion—with a PINSFVMomentumDiffusion kernel,

[FVKernels]
  [vel_x_viscosity]
    type = PINSFVMomentumDiffusion
    variable = vel_x
    momentum_component = x
  []
[]

And the fourth term–the friction term—with a PNSFVMomentumFriction kernel,

[FVKernels]
  [u_friction]
    type = PINSFVMomentumFriction
    variable = vel_x
    Darcy_name = 'Darcy_coefficient'
    Forchheimer_name = 'Forchheimer_coefficient'
    momentum_component = 'x'
  []
[]

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,

[FVKernels]
  [buoyancy_boussinesq]
    type = PINSFVMomentumBoussinesq
    variable = vel_y
    ref_temperature = ${inlet_T_fluid}
    momentum_component = 'y'
    alpha_name = 'alpha_b'
  []

  [gravity]
    type = PINSFVMomentumGravity
    variable = vel_y
    momentum_component = 'y'
  []
[]

Kernel syntax (legacy) for the 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 first term of Eq. (2)—the energy time derivative—is captured by the kernel,

[FVKernels]
  [temp_fluid_time]
    type = PINSFVEnergyTimeDerivative
    variable = temp_fluid
    cp = 'cp'
    is_solid = false
  []
[]

The second term of Eq. (2)—energy advection—is expressed by,

[FVKernels]
  [temp_fluid_advection]
    type = PINSFVEnergyAdvection
    variable = temp_fluid
    advected_quantity = 'rho_cp_temp'
  []
[]

The third term—the effective diffusion of heat—corresponds to the kernel,

[FVKernels]
  [temp_fluid_conduction]
    type = PINSFVEnergyAnisotropicDiffusion
    variable = temp_fluid
    effective_diffusivity = false
    kappa = 'kappa'
  []
[]

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

[FVKernels]
  [temp_solid_to_fluid]
    type = PINSFVEnergyAmbientConvection
    variable = temp_fluid
    is_solid = false
    h_solid_fluid = 'alpha'
  []
[]

Legacy explicit boundary conditions syntax

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.

[FVBCs]
  [inlet_vel_y]
    type = INSFVInletVelocityBC
    variable = vel_y
    function = ${inlet_vel_y}
    boundary = 'bed_horizontal_bottom'
  []

  [inlet_temp_fluid]
    type = FVDirichletBC
    variable = temp_fluid
    value = '${fparse inlet_T_fluid}'
    boundary = 'bed_horizontal_bottom'
  []
[]

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.

[FVBCs]
  [free-slip-wall-x]
    type = INSFVNaturalFreeSlipBC
    boundary = 'bed_left barrel_wall'
    variable = vel_x
    momentum_component = x
    type = INSFVNaturalFreeSlipBC
    boundary = 'bed_left barrel_wall'
    variable = vel_x
    momentum_component = x
  []
[]

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.

[FVBCs]
  [outlet_p]
    type = INSFVOutletPressureBC
    variable = pressure
    function = 2e5 # not too far from atm for matprop evaluations
    boundary = 'bed_horizontal_top OR_horizontal_top plenum_top'
  []
[]

Material properties

[FunctorMaterials]
  [ins_fv]
    type = INSFVPrimitiveSuperficialVarMaterial
    T_fluid = 'temp_fluid'
    block = ${blocks_fluid}
    p_constant = 1e5
    T_constant = 900
  []
[]

References