MSFR Griffin-Pronghorn Model

Contact: Mauricio Tano, mauricio.tanoretamales.at.inl.gov

Model link: Griffin-Pronghorn Steady-State Model

The MultiApp system is used to separate the neutronics and the fluid dynamics problems. The fluid system is solved by the subapp and it uses the run_ns.i input files. (Here "ns" is an abbreviation for Navier-Stokes.)

The fluid system includes conservation equations for fluid mass, momentum, and energy as well as the conservation of delayed neutron precursors.

Conservation of fluid mass

The conservation of mass is,

where is the fluid density and is the velocity vector.

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 equation is automatically added by the NavierStokesFV action when selecting the incompressible model.


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

Legacy syntax for the mass equation is included here.

Conservation of fluid momentum

This system also includes the conservation of momentum in the -direction. A fairly general form of the steady-state condition is,

where is the component of the velocity, is the pressure, is the component of the viscous friction force, and is the gravity vector.

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

Practical simulations require modifications to the momentum equations in order to model the effects of turbulence without explicitly resolving the turbulent structures. Here, we will apply the Reynolds-averaging procedure and the Boussinesq hypothesis so that the effect of turbulent momentum transfer is modeled with a term analogous to viscous shear,

where is the turbulent viscosity.

Here, a zero-equation model based on the mixing length model is used. In this model the turbulent viscosity is defined as: and

The standard Prandtl's mixing length model dictates that has a linear dependence on the distance to the nearest wall. However, for this simulation we implement a capped mixing length model (Escudier, 1966) that defines the mixing length as

where is the Von Karman constant, as in Escudier's model and has length units and represents the thickness of the velocity boundary layer.

The viscous friction is treated with two different models for different regions of the reactor. The bulk of the friction is expected to occur in the heat exchanger region where there is a large surface area in contact with the salt to facilitate heat transfer. Rather than explicitly modeling the heat exchanger geometry, a simple, tunable model will be used, where is a tunable volumetric friction coefficient, which is selected to match the desired pressure drop, in conjunction with the pump head value.

A typical viscous friction model could be used for the regions outside of the heat exchanger. However, the viscous effects will be negligible due to the large, uniform eddy viscosity model used here. Consequently, the viscous force will be neglected outside of the heat exchanger,

By convention, we must collect all of the terms on one side of the equation. This gives the form that is implemented for the MSFR model,

(1)

The first few terms in Eq. (1), including the advection of momentum, the pressure gradient and the laminar diffusion are automatically handled by the NavierStokesFV action.

The mixing length turbulence model is added by specifying additional parameters to the action. This basic model is currently the most limiting issue in terms of accuracy.


[Modules]
  [NavierStokesFV]
    ...
    # Turbulence parameters
    turbulence_handling = 'mixing-length'
    turbulent_prandtl = ${Pr_t}
    von_karman_const = ${von_karman_const}
    mixing_length_delta = 0.1
    mixing_length_walls = 'shield_wall reflector_wall'
    mixing_length_aux_execute_on = 'initial'
    ...

The friction in the heat exchanger is specified by passing the following parameters to the NavierStokesFV action. A quadratic friction model with a coefficient tuned to obtain the desired mass flow rate is selected.


[Modules]
  [NavierStokesFV]
    ...
    # Heat exchanger
    friction_blocks = 'hx'
    friction_types = 'FORCHHEIMER'
    friction_coeffs = ${friction}
    ...

This model does not include a friction force for the other blocks so no kernel needs to be specified for those blocks.

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

where is the pump head driving the flow, is the expansion coefficient, is the fluid temperature, and is a reference temperature value. The pump head is tuned such that the imposed mass flow rate is ~18500 kg/s.

For each kernel describing the -momentum equation, there is a corresponding kernel for the -momentum equation. They are similarly added by the NavierStokesFV. An additional kernel is added to add a volumetric pump component to the system:

[FVKernels]
  [pump]
    type = INSFVBodyForce
    variable = vel_y
    functor = ${pump_force}
    block = 'pump'
    momentum_component = 'y'
  []
[]
(msr/msfr/steady/run_ns.i)

The Boussinesq buoyancy approximation is added to the model using this parameter:


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

Boundary conditions include standard velocity wall functions at the walls to account for the non-linearity of the velocity in the boundary layer given the coarse mesh and symmetry at the center axis of the MSFR. They are added by the action based on user-input parameters:


[Modules]
  [NavierStokesFV]
    ...
    # Boundary conditions
    wall_boundaries = 'shield_wall reflector_wall fluid_symmetry'
    momentum_wall_types = 'wallfunction wallfunction symmetry'
    ...

Auxkernels are used to compute the wall shear stress obtained by the standard wall function model, the dimensionless wall distance and the value for the eddy viscosity. These are used for analysis purposes.

[AuxKernels]
  [mixing_len]
    type = WallDistanceMixingLengthAux
    walls = 'shield_wall reflector_wall'
    variable = mixing_len
    execute_on = 'initial'
    block = 'fuel'
    von_karman_const = ${von_karman_const}
    delta = 0.1 # capping parameter (m)
  []
  [wall_shear_stress]
    type = WallFunctionWallShearStressAux
    variable = wall_shear_stress
    walls = 'shield_wall reflector_wall'
    block = 'fuel'
    mu = 'total_viscosity'
  []
  [wall_yplus]
    type = WallFunctionYPlusAux
    variable = wall_yplus
    walls = 'shield_wall reflector_wall'
    block = 'fuel'
    mu = 'total_viscosity'
  []
  [turbulent_viscosity]
    type = INSFVMixingLengthTurbulentViscosityAux
    variable = eddy_viscosity
    block = 'fuel pump hx'
    execute_on = 'initial'
  []
[]
(msr/msfr/steady/legacy/run_ns.i)

Legacy syntax for the momentum equation is included here.

Conservation of fluid energy

The steady-state conservation of energy can be expressed as, where is the fluid specific enthalpy, is the thermal conductivity, and is the volumetric heat generation rate.

Here it is expected that the energy released from nuclear reactions will be very large compared to pressure work terms. Consequently, we will use the simplified form,

As with the momentum equations, a practical solver requires a model for the turbulent transport of energy, where is the eddy diffusivity for heat. It is related to the eddy diffusivity for momentum (i.e. the kinematic eddy viscosity) as where is the turbulent Prandtl number.

As with the molecular viscosity in the momentum equations, the simple turbulence model used here will overwhelm the thermal conductivity. Neglecting that term and moving the heat generation to the left-hand-side of the equation gives,

The heat generation due to nuclear reactions is computed by the neutronics solver, and this distribution will be used directly in the term. The effect of the heat exchanger will also be included in the term as a volumetric heat loss per the model, where is a coefficient (equal to the surface area density times the heat transfer coefficient) and is the temperature of the coolant on the secondary side of the heat exchanger.

Note that all material properties, including and , are assumed constant. It will therefore be convenient to move the factors outside of the divergence operators and divide the entire equation by that factor, (2)

Eq. (2) is the final equation that is implemented in this model. It is added to the equation system using a parameter in the NavierStokesFV action:


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

The first term—energy advection— is automatically added by the action syntax. The second term—turbulent diffusion of heat— is added based on the momentum turbulent term, with the Prandtl number also to be specified to adjust the diffusion coefficient.

The third term—heat source and loss—is covered by two kernels. First, the nuclear heating computed by the neutronics solver is included with,


[Modules]
  [NavierStokesFV]
    ...
    # Heat source
    external_heat_source = power_density
    ...

And second, the heat loss through the heat exchanger is implemented with the kernel,


[Modules]
  [NavierStokesFV]
    ...
    # Heat exchanger
    ambient_convection_blocks = 'hx'
    ambient_convection_alpha = ${fparse 600 * 20e3} # HX specifications
    ambient_temperature = ${T_HX}
    ...

The boundary conditions are added similarly, by simply specifying 0 heat flux for symmetry and adiabatic boundaries.


[Modules]
  [NavierStokesFV]
    ...
    energy_wall_types = 'heatflux heatflux heatflux'
    energy_wall_function = '0 0 0'

Legacy syntax for the energy equation is included here.

Conservation of delayed neutron precursors

The steady-state conservation of delayed neutron precursors (DNP) can be expressed as, (3) where , and are the concentration, the decay constant and the production fraction per neutron generated from fission of the -th group of DNP respectively. is the total number of groups of DNP. is the fission neutron production rate. is the eddy diffusivity for DNP, where is the turbulnt viscosity and is the turbulent Schmidt number. It is noted the DNP term in the neutron transport equation is scaled with -effective as the prompt fission term for steady-state eigenvalue calculations, where -effective (the eigenvalue) is part of the solution.

It is added to the equation system using a parameter in the NavierStokesFV action:


[Modules]
  [NavierStokesFV]
    ...
    add_scalar_equation = true
    ...

The terms of Eq. (2) are added with


[Modules]
  [NavierStokesFV]
    ...
    # Precursor advection, diffusion and source term
    passive_scalar_names = 'c1 c2 c3 c4 c5 c6'
    Sc_t = '${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t} ${Sc_t}'
    passive_scalar_coupled_source = 'fission_source c1; fission_source c2; fission_source c3;
                                     fission_source c4; fission_source c5; fission_source c6;'
    passive_scalar_coupled_source_coeff = '${beta1} ${lambda1_m}; ${beta2} ${lambda2_m};
                                           ${beta3} ${lambda3_m}; ${beta4} ${lambda4_m};
                                           ${beta5} ${lambda5_m}; ${beta6} ${lambda6_m}'

Neutronics

With Griffin, the process of converting the basic conservation equations into MOOSE variables and kernels is automated with the TransportSystems block:

[TransportSystems]
  particle = neutron
  equation_type = eigenvalue

  G = 6

  VacuumBoundary = 'outer'

  [diff]
    scheme = CFEM-Diffusion
    n_delay_groups = 6
    external_dnp_variable = 'dnp'
    family = LAGRANGE
    order = FIRST
    fission_source_aux = true

    # For PJFNKMO
    assemble_scattering_jacobian = true
    assemble_fission_jacobian = true
  []
[]
(msr/msfr/steady/run_neutronics.i)

Details about neutron transport equations can be found in Griffin theory manual.

Here we are specifying an eigenvalue neutronics problem using 6 energy groups (G = 6) solved via the diffusion approximation with a continuous finite element discretization scheme (scheme = CFEM-Diffusion).

Note the external_dnp_variable = 'dnp' parameter. This is a special option needed for liquid-fueled MSRs which signals that the conservation equations for DNP will be handled "externally" from the default Griffin implementation which assumes that the precursors do not have the turbulent treatment. This parameter is referencing the dnp array auxiliary variable defined by the MSRNeutronicsFluidCoupling class, which in this case is an array auxiliary variable with 6 components, corresponding to the 6 DNP groups used here. Support within the Framework for array variables is somewhat limited. For example, not all of the multiapp transfers work with array variables, and the Navier-Stokes module does not include the kernels that are needed to advect an array variable. For this reason, there is also a separate auxiliary variable for each DNP group, which are also created within MSRNeutronicsFluidCoupling.

The MSRNeutronicsFluidCoupling block,

[MSRNeutronicsFluidCoupling]
  fluid_blocks = 'fuel pump hx'
  solid_blocks = 'reflector shield'
  n_dnp = 6
  use_transient_multi_app = false
  fluid_app_name = ns
  fluid_input_file = run_ns.i
  initial_fluid_temperature = 873.15
  fluid_temperature_name_neutronics_app = tfuel
  fluid_temperature_name_solid_suffix = _constant
  fluid_temperature_name_fluid_app = T_fluid
  dnp_name_prefix = c
  dnp_array_name = dnp
[]
(msr/msfr/steady/run_neutronics.i)

creates a variety of objects needed on the neutronics side of the neutronics-fluid coupling for MSRs:

  • aux variables for the aforementioned DNP concentration variables and additionally, fluid temperature,

  • the aux kernel to form the DNP array from its components,

  • initial conditions for the fluid temperature if not restarting the simulation,

  • a MultiApp (TransientMultiApp for transients, FullSolveMultiApp otherwise) to run the fluid application, and

  • transfers for the DNP concentrations, power density, fission source, and fluid temperature.

The run_ns.i subapp is responsible for computing the precursor distributions, and the distributions are transferred from the subapp to the main app by MSRNeutronicsFluidCoupling and then internally copied from the individual DNP group variables into the DNP array variable using an aux kernel.

Also note that solving the neutronics problem requires a set of multigroup cross sections. Generating cross sections is a topic that is left outside the scope of this example. A set has been generated for the MSFR problem and stored in the repository using Griffin's ISOXML format. These cross sections are included by the blocks,

[Materials]
  [fuel]
    type = CoupledFeedbackNeutronicsMaterial
    library_name = 'msfr_xs'
    library_file = '../mgxs/msfr_xs.xml'
    grid_names = 'tfuel'
    grid_variables = 'tfuel'
    plus = true
    isotopes = 'pseudo'
    densities = '1.0'
    is_meter = true
    material_id = 2
    block = 'fuel pump hx'
  []
  [shield]
    type = CoupledFeedbackNeutronicsMaterial
    library_name = 'msfr_xs'
    library_file = '../mgxs/msfr_xs.xml'
    grid_names = 'tfuel'
    grid_variables = 'tfuel_constant'
    isotopes = 'pseudo'
    densities = '1.0'
    is_meter = true
    material_id = 5
    block = 'shield'
  []
  [reflector]
    type = CoupledFeedbackNeutronicsMaterial
    library_name = 'msfr_xs'
    library_file = '../mgxs/msfr_xs.xml'
    grid_names = 'tfuel'
    grid_variables = 'tfuel_constant'
    isotopes = 'pseudo'
    densities = '1.0'
    is_meter = true
    material_id = 1
    block = 'reflector'
  []
[]
(msr/msfr/steady/run_neutronics.i)

CoupledFeedbackNeutronicsMaterial is able to use the temperature transferred from the fluid system for evaluating multigroup cross sections based on a table lookup on element quadrature points to bring in the feedback effect. It also has the capability of adjusting cross sections based on fluid density. In this model, the fluid density change is negligible.

The neutronics solution is normalized to the rated power 3000 MW with the PowerDensity input block:

[PowerDensity]
  power = 3e9
  power_density_variable = power_density
  power_scaling_postprocessor = power_scaling
  family = MONOMIAL
  order = CONSTANT
[]
(msr/msfr/steady/run_neutronics.i)

The power density is evaluated with the normalized neutronics solution. It provides the source of the fluid energy equation. Because the fluid energy equation is discretized with FV, we evaluate the power density variable with constant monomial.

The calculation is driven by Eigenvalue, an executioner available in the MOOSE framework. The PJFNKMO option for the solve_type parameter is able to drive the eigenvalue calculation with the contribution of DNP to the neutron transport equation as an external source scaled with -effective.

[Executioner]
  type = Eigenvalue
  solve_type = PJFNKMO

  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart '
  petsc_options_value = 'hypre boomeramg 50'
  l_max_its = 50

  free_power_iterations = 4 # important to obtain fundamental mode eigenvalue

  nl_abs_tol = 1e-9
  fixed_point_abs_tol = 1e-9
  fixed_point_max_its = 20
[]
(msr/msfr/steady/run_neutronics_base.i)

References

  1. M.P Escudier. The Distribution of Mixing-Length in Turbulent Flows Near Walls. PhD thesis, Imperial College, 1966.[BibTeX]