The MOOSE Navier-Stokes Module

July 2025

Overview

Main capabilities:

  • Fluid types:

    • Incompressible

    • Weakly-compressible

    • Compressible

    • Two-phase mixture model

  • Flow regimes:

    • Laminar

    • Turbulent

  • Flow types:

    • Free flows

    • Flows in porous media describing homogenized structures

Utilizes the following techniques:

  • Stabilized Continuous Finite Element Discretization (maintained, but not intensely developed) (Peterson et al. (2018))

  • Finite Volume Discretization with nonlinear solver (Lindsay et al. (2023))

  • Finite Volume Discretization with linear solver

Applicability of each solver

Nonlinear (Newton method) finite volume

  • 2D problems

  • only solver with support for porous medium

  • turbulence models (mixing-length, k-)

  • fully compressible flow problems

Linear (SIMPLE/PIMPLE) finite volume

  • 3D problems

You can check that the combination of models you want is supported in Table 2 on this page The nonlinear (SIMPLE/PIMPLE) finite volume is being replaced by the linear implementation. It is currently still only used for two-phase flow studies.

The Porous Navier-Stokes equations (Vafai (2015), Radman et al. (2021))

  • Conservation of mass:

(2)
  • Conservation of (linear) momentum:

(3)

  • - porosity

  • - real (or interstitial) velocity

  • - superficial velocity

  • - density

  • - effective dynamic viscosity (laminar + turbulent)

  • - pressure

  • - friction tensor (using correlations)

  • If then and we get the original Navier-Stokes equations back.

The Porous Navier-Stokes equations (contd.)

  • These are supplemented by the conservation of energy:

    Assumptions: Weakly-compressible fluid, kinetic energy of the fluid can be neglected

(1)
  • - specific heat

  • - fluid temperature

  • - effective thermal conductivity

  • - volumetric heat-transfer coefficient (using correlations)

  • - solid temperature

  • - volumetric source term

The Porous Navier-Stokes equations (contd.)

The equations are extended with initial and boundary conditions:

  • Few examples:

    • Velocity Inlet:

    • Pressure Outlet:

    • No-slip, insulated walls:

    • Heat flux on the walls:

    • Free slip on walls:

    • ...

Useful Acronyms

Good to know before navigating source code and documentation:

  • NS - Navier-Stokes

  • FV - Finite Volume

  • I - Incompressible

  • WC - Weakly-Compressible

  • C - Compressible

  • P - Porous-Medium

Few Examples:

  • PINSFV: Porous-Medium Incompressible Navier-Stokes Finite Volume

  • WCNSFV: Weakly-Compressible Navier-Stokes Finite Volume

  • WCNSLinearFV: Weakly-Compressible Navier-Stokes Linear Finite Volume

Building Input Files

The building blocks in MOOSE for terms in the PDEs are Kernels for FE, FVKernels for FV, and LinearFVKernels for SIMPLE-FV:

The building blocks in MOOSE for boundary conditions are BCs for FE, FVBCs for Newton-FV, and LinearFVBCs for SIMPLE-FV:

  • Inlet Velocity:

  • Inlet Temperature:

  • Outlet pressure:

  • Heatflux:

  • Freeslip:

  • ...

A Simple Example: Laminar Free Flow in a Channel

Let us consider the following (fictional) material properties:

  • (we are using incompressible formulation in this case)

Detailed Input File (nonlinear/Newton FV)

mu = 1.1
rho = 1.1
l = 2
U = 1
advected_interp_method = 'average'
velocity_interp_method = 'rc'

[GlobalParams]
  rhie_chow_user_object = 'rc'
[]

[UserObjects]
  [rc]
    type = INSFVRhieChowInterpolator
    u = vel_x
    v = vel_y
    pressure = pressure
  []
[]

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = 10
    ymin = ${fparse -l / 2}
    ymax = ${fparse l / 2}
    nx = 100
    ny = 20
  []
  uniform_refine = 0
[]

[Variables]
  [vel_x]
    type = INSFVVelocityVariable
    initial_condition = 1
  []
  [vel_y]
    type = INSFVVelocityVariable
    initial_condition = 1
  []
  [pressure]
    type = INSFVPressureVariable
  []
[]

[FVKernels]
  [mass]
    type = INSFVMassAdvection
    variable = pressure
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = ${rho}
  []

  [u_advection]
    type = INSFVMomentumAdvection
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = ${rho}
    momentum_component = 'x'
  []
  [u_viscosity]
    type = INSFVMomentumDiffusion
    variable = vel_x
    mu = ${mu}
    momentum_component = 'x'
  []
  [u_pressure]
    type = INSFVMomentumPressure
    variable = vel_x
    momentum_component = 'x'
    pressure = pressure
  []

  [v_advection]
    type = INSFVMomentumAdvection
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    velocity_interp_method = ${velocity_interp_method}
    rho = ${rho}
    momentum_component = 'y'
  []
  [v_viscosity]
    type = INSFVMomentumDiffusion
    variable = vel_y
    mu = ${mu}
    momentum_component = 'y'
  []
  [v_pressure]
    type = INSFVMomentumPressure
    variable = vel_y
    momentum_component = 'y'
    pressure = pressure
  []
[]

[FVBCs]
  [inlet-u]
    type = INSFVInletVelocityBC
    boundary = 'left'
    variable = vel_x
    function = '${U}'
  []
  [inlet-v]
    type = INSFVInletVelocityBC
    boundary = 'left'
    variable = vel_y
    function = '0'
  []
  [walls-u]
    type = INSFVNoSlipWallBC
    boundary = 'top bottom'
    variable = vel_x
    function = 0
  []
  [walls-v]
    type = INSFVNoSlipWallBC
    boundary = 'top bottom'
    variable = vel_y
    function = 0
  []
  [outlet_p]
    type = INSFVOutletPressureBC
    boundary = 'right'
    variable = pressure
    function = '0'
  []
[]

[Executioner]
  type = Steady
  solve_type = 'NEWTON'
  nl_rel_tol = 1e-12
[]

[Preconditioning]
  active = FSP
  [FSP]
    type = FSP
    # It is the starting point of splitting
    topsplit = 'up' # 'up' should match the following block name
    [up]
      splitting = 'u p' # 'u' and 'p' are the names of subsolvers
      splitting_type  = schur
      # Splitting type is set as schur, because the pressure part of Stokes-like systems
      # is not diagonally dominant. CAN NOT use additive, multiplicative and etc.
      #
      # Original system:
      #
      # | Auu Aup | | u | = | f_u |
      # | Apu 0   | | p |   | f_p |
      #
      # is factorized into
      #
      # |I             0 | | Auu  0|  | I  Auu^{-1}*Aup | | u | = | f_u |
      # |Apu*Auu^{-1}  I | | 0   -S|  | 0  I            | | p |   | f_p |
      #
      # where
      #
      # S = Apu*Auu^{-1}*Aup
      #
      # The preconditioning is accomplished via the following steps
      #
      # (1) p* = f_p - Apu*Auu^{-1}f_u,
      # (2) p = (-S)^{-1} p*
      # (3) u = Auu^{-1}(f_u-Aup*p)

      petsc_options_iname = '-pc_fieldsplit_schur_fact_type  -pc_fieldsplit_schur_precondition -ksp_gmres_restart -ksp_rtol -ksp_type'
      petsc_options_value = 'full                            selfp                             300                1e-4      fgmres'
    []
    [u]
      vars = 'vel_x vel_y'
      petsc_options_iname = '-pc_type -pc_hypre_type -ksp_type -ksp_rtol -ksp_gmres_restart -ksp_pc_side'
      petsc_options_value = 'hypre    boomeramg      gmres    5e-1      300                 right'
    []
    [p]
      vars = 'pressure'
      petsc_options_iname = '-ksp_type -ksp_gmres_restart -ksp_rtol -pc_type -ksp_pc_side'
      petsc_options_value = 'gmres    300                5e-1      jacobi    right'
    []
  []
  [SMP]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_shift_type'
    petsc_options_value = 'lu       NONZERO'
  []
[]

[Outputs]
  print_linear_residuals = true
  print_nonlinear_residuals = true
  [out]
    type = Exodus
    hide = 'Re lin cum_lin'
  []
  [perf]
    type = PerfGraphOutput
  []
[]

[Postprocessors]
  [Re]
    type = ParsedPostprocessor
    expression = '${rho} * ${l} * ${U}'
  []
  [lin]
    type = NumLinearIterations
  []
  [cum_lin]
    type = CumulativeValuePostprocessor
    postprocessor = lin
  []
[]
(modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/2d-rc-no-slip.i)

Detailed Input File (linear/SIMPLE FV)

mu = 2.6
rho = 1.0
advected_interp_method = 'average'

[Mesh]
  [mesh]
    type = CartesianMeshGenerator
    dim = 2
    dx = '0.3'
    dy = '0.3'
    ix = '3'
    iy = '3'
  []
[]

[Problem]
  linear_sys_names = 'u_system v_system pressure_system'
  previous_nl_solution_required = true
[]

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    initial_condition = 0.5
    solver_sys = u_system
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0.0
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0.2
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'x'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
  []
  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'y'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []
  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []
[]

[LinearFVBCs]
  [inlet-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_x
    functor = '1.1'
  []
  [inlet-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_y
    functor = '0.0'
  []
  [walls-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom'
    variable = vel_x
    functor = 0.0
  []
  [walls-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom'
    variable = vel_y
    functor = 0.0
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'right'
    variable = pressure
    functor = 1.4
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = vel_x
    use_two_term_expansion = false
    boundary = right
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = vel_y
    use_two_term_expansion = false
    boundary = right
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-10
  pressure_l_abs_tol = 1e-10
  momentum_l_tol = 0
  pressure_l_tol = 0
  rhie_chow_user_object = 'rc'
  momentum_systems = 'u_system v_system'
  pressure_system = 'pressure_system'
  momentum_equation_relaxation = 0.8
  pressure_variable_relaxation = 0.3
  num_iterations = 100
  pressure_absolute_tolerance = 1e-10
  momentum_absolute_tolerance = 1e-10
  momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
  momentum_petsc_options_value = 'hypre boomeramg'
  pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
  pressure_petsc_options_value = 'hypre boomeramg'
  print_fields = false
[]

[Outputs]
  exodus = true
[]
(modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/2d-velocity-pressure.i)

Rhie-Chow Interpolation

  • When using linear interpolation for the advecting velocity and the pressure gradient, we encounter a checker-boarding effect for the pressure field. Neighboring cell center solution values are decoupled, which is undesirable.

  • Solution for this issue: Rhie-Chow interpolation for the advecting velocity. For more information on the issue and the derivation of the interpolation method, see Moukalled et al. (2016).

The Navier-Stokes Physics

A simplified syntax has been developed, it relies on the Physics system in MOOSE. Each Physics is in charge of creating one equation. The input below shows the flow and energy conservation equations. Only the finite volume discretization, for the Newton and PIMPLE solvers, support this syntax.

mu = 1
rho = 1
k = 1e-3
cp = 1
u_inlet = 1
T_inlet = 200
h_cv = 1.0

[Mesh]
  [mesh]
    type = CartesianMeshGenerator
    dim = 2
    dx = '5 5'
    dy = '1.0'
    ix = '50 50'
    iy = '20'
    subdomain_id = '1 2'
  []
[]

[Variables]
  [T_solid]
    type = MooseVariableFVReal
  []
[]

[AuxVariables]
  [porosity]
    type = MooseVariableFVReal
    initial_condition = 0.5
  []
[]

[Physics]
  [NavierStokes]
    [Flow]
      [flow]
        compressibility = 'incompressible'
        porous_medium_treatment = true

        density = ${rho}
        dynamic_viscosity = ${mu}
        porosity = 'porosity'

        initial_velocity = '${u_inlet} 1e-6 0'
        initial_pressure = 0.0

        inlet_boundaries = 'left'
        momentum_inlet_types = 'fixed-velocity'
        momentum_inlet_function = '${u_inlet} 0'

        wall_boundaries = 'top bottom'
        momentum_wall_types = 'noslip symmetry'

        outlet_boundaries = 'right'
        momentum_outlet_types = 'fixed-pressure'
        pressure_function = '0.1'

        mass_advection_interpolation = 'average'
        momentum_advection_interpolation = 'average'
      []
    []
    [FluidHeatTransfer]
      [heat]
        thermal_conductivity = ${k}
        specific_heat = ${cp}

        # Reference file sets effective_conductivity by default that way
        # so the conductivity is multiplied by the porosity in the kernel
        effective_conductivity = false

        initial_temperature = 0.0

        energy_inlet_types = 'heatflux'
        energy_inlet_function = '${fparse u_inlet * rho * cp * T_inlet}'

        energy_wall_types = 'heatflux heatflux'
        energy_wall_function = '0 0'

        ambient_convection_alpha = ${h_cv}
        ambient_temperature = 'T_solid'

        energy_advection_interpolation = 'average'
      []
    []
  []
[]

[FVKernels]
  [solid_energy_diffusion]
    type = FVDiffusion
    coeff = ${k}
    variable = T_solid
  []
  [solid_energy_convection]
    type = PINSFVEnergyAmbientConvection
    variable = 'T_solid'
    is_solid = true
    T_fluid = 'T_fluid'
    T_solid = 'T_solid'
    h_solid_fluid = ${h_cv}
  []
[]

[FVBCs]
  [heated-side]
    type = FVDirichletBC
    boundary = 'top'
    variable = 'T_solid'
    value = 150
  []
[]

[Executioner]
  type = Steady
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu NONZERO'
  nl_rel_tol = 1e-14
[]

# Some basic Postprocessors to examine the solution
[Postprocessors]
  [inlet-p]
    type = SideAverageValue
    variable = pressure
    boundary = 'left'
  []
  [outlet-u]
    type = SideAverageValue
    variable = superficial_vel_x
    boundary = 'right'
  []
  [outlet-temp]
    type = SideAverageValue
    variable = T_fluid
    boundary = 'right'
  []
  [solid-temp]
    type = ElementAverageValue
    variable = T_solid
  []
[]

[Outputs]
  exodus = true
  csv = false
[]
(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/heated/2d-rc-heated-physics.i)

Recommendations for Building Input Files

  • Check that the combination of models you want is supported in Table 2 there

  • Start from an example. If no example, find a test that looks similar.

  • If possible, use the Navier Stokes finite volume Physics syntax

  • Use Rhie-Chow interpolation for the advecting velocity

    • Other interpolation techniques may lead to checkerboarding/instability

  • Start with first-order advected-interpolation schemes (e.g. upwind)

  • Make sure that the pressure is pinned for incompressible/weakly-compressible simulations in closed systems

  • For monolithic solvers (the default at the moment) use a variant of LU preconditioner

  • For complex monolithic systems monitor the residuals of every variable

  • Try to keep the number of elements relatively low (limited scaling with LU preconditioner). Larger problems can be solved in a reasonable wall-time as segregated-solve techniques are introduced

  • Try to utilize porosity_smoothing_layers or a higher value for the consistent_scaling if you encounter oscillatory behavior in case of simulations using porous medium

Because the complexity of the LU preconditioner is in the worst case scenario, its performance in general largely depends on the matrix sparsity pattern, and it is much more expensive compared with an iterative approach.

Validation studies and additional examples

Turbulence modeling for Molten Salt Reactors

Open-Pronghorn flow test bed (to be released soon)

References

  1. Alexander Lindsay, Guillaume Giudicelli, Peter German, John Peterson, Yaqi Wang, Ramiro Freile, David Andrs, Paolo Balestra, Mauricio Tano, Rui Hu, and others. Moose navier–stokes module. SoftwareX, 23:101503, 2023.
  2. Fadl Moukalled, L Mangani, Marwan Darwish, and others. The finite volume method in computational fluid dynamics. Volume 6. Springer, 2016.
  3. John W. Peterson, Alexander D. Lindsay, and Fande Kong. Overview of the incompressible navier–stokes simulation capabilities in the moose framework. Advances in Engineering Software, 119:68–92, 2018.
  4. Stefan Radman, Carlo Fiorina, and Andreas Pautz. Development of a novel two-phase flow solver for nuclear reactor analysis: algorithms, verification and implementation in openfoam. Nuclear Engineering and Design, 379:111178, 2021.
  5. Kambiz Vafai. Handbook of porous media. Crc Press, 2015.