SIMPLE

Solves the Navier-Stokes equations using the SIMPLE algorithm and linear finite volume variables.

Overview

This executioner is based on the algorithm proposed by Patankar and Spalding (1983). The algorithm is based on the splitting of operators and successive correction for the momentum and pressure fields. The formulation implemented in MOOSE has been presented in Jasak (1996) and Juretic (2005). See also the examples and derivations in Moukalled et al. (2016). The concept relies on deriving a pressure equation using the discretized form of the momentum equations together with the continuity constraint. Let's take the steady-state incompressible Navier-Stokes equations in the following form:

(ρuu)(μeff(u+uT))=p+G.\nabla \cdot \left(\rho \vec{u} \otimes \vec{u}\right) - \nabla \cdot \left(\mu_\text{eff} \left(\nabla\vec{u} +\nabla \vec{u}^T \right)\right) = -\nabla p + \vec{G}.(1)(ρu)=0.\nabla \cdot \left(\rho \vec{u}\right) = 0.(2)

Where u\vec{u} denotes the velocity, pp the pressure, ρ\rho the density, and μeff\mu_\text{eff} the effective dynamic viscosity which potentially includes the contributions of eddy viscosity derived from turbulence models. Term G\vec{G} expresses a volumetric source term which can be potentially velocity-dependent. As a first step, we assume that we have a guess for the pressure field, therefore the gradient is known. Furthermore, we assume that the advecting velocity field is known from the previous iteration. By explicitly showing the iteration index, Eq. (1) and Eq. (2) become:

(ρun1un)(μeff(un+un,T))=pn1+G(un1,un).\nabla \cdot \left(\rho \vec{u}^{n-1} \otimes \vec{u}^n\right) - \nabla \cdot \left(\mu_\text{eff} \left(\nabla\vec{u}^n +\nabla \vec{u}^{n,T}\right)\right) = -\nabla p^{n-1} + \vec{G}(\vec{u}^{n-1},\vec{u}^{n}).(3)(ρun)=0.\nabla \cdot \left(\rho \vec{u}^n\right) = 0.(4)

At this point, we should note that the finite volume discretization in MOOSE uses a collocated formulation which has an advantage of being flexible for unstructured meshes. However, in certain scenarios it can exhibit numerical pressure checker-boarding due to the discretization of the pressure gradient and continuity terms. A common approach for tackling this issue is the utilization of the Rhie-Chow interpolation method (See Rhie and Chow (1983) and Moukalled et al. (2016) for a detailed explanation). This means that the face velocities (or face fluxes) are determined using pressure corrections. As we will see later, due to this behavior, the iteration between pressure and velocity will in fact be an iteration between pressure and face velocity. Nevertheless, to keep this in mind we add a subscript to the advecting velocity in our formulation:

(ρuRCn1un)(μeff(un+un,T))=pn1+G(un1,un).\nabla \cdot \left(\rho \vec{u}^{n-1}_{RC} \otimes \vec{u}^n\right) - \nabla \cdot \left(\mu_\text{eff} \left(\nabla\vec{u}^n +\nabla \vec{u}^{n,T}\right)\right) = -\nabla p^{n-1} + \vec{G}(\vec{u}^{n-1},\vec{u}^{n}).(5)(ρuRCn)=0.\nabla \cdot \left(\rho \vec{u}^n_{RC}\right) = 0.(6)

Next, we split the operator acting on u\vec{u} in the momentum equation into two components: a component that incorporates effects that result in contributions to the diagonal of a soon-to-be-generated system matrix and another component that contains everything else. With this in mind, we can rewrite the equation the following, semi-discretized way:

A(un1)un+H(un)=pn1,A(\vec{u}^{n-1})\vec{u}^n + H(\vec{u}^{n}) = -\nabla p^{n-1},(7)

where A(un1)A(\vec{u}^{n-1}) is the diagonal contribution, and H(un1)H(\vec{u}^{n-1}) includes the off-diagonal contributions multiplied by the solution together with any additional volumetric source and sink terms (i.e. the discretized forms of G\vec{G}). One can solve this equation to obtain a new guess for the velocity field. This guess, however, will not respect the continuity equation, therefore we need to correct it. For this, a pressure equation is derived from the following formulation:

Aun+H(un)=pn.A\vec{u}^n + H(\vec{u}^{n}) = -\nabla p^n.(8)

By applying the inverse of the diagonal operator (a very cheap process computationally), we arrive to the following expression:

un+A1H(un)=A1pn,\vec{u}^n + A^{-1}H(\vec{u}^{n}) = -A^{-1}\nabla p^n,(9)

By applying the continuity equation onto un\vec{u}^n (which is a constraint) and assuming that the Rhie-Chow interpolation is used for the velocity, we arrive to a Poisson equation for pressure:

(A1H(un))=(A1pn).\nabla \cdot \left(A^{-1}H(\vec{u}^{n})\right) = -\nabla \cdot \left(A^{-1}\nabla p^n\right).(10)

This equation is solved for a pressure which can be used to correct the face velocities in a sense that they respect the continuity equation. This correction already involves a Rhie-Chow interpolation, considering that the A1HA^{-1}H and A1A^{-1} fields are interpolated to the faces in a discretized form:

uRC,fn+1=(A1H(un))f(A1)fpfn.\vec{u}^{n+1}_{RC,f} = - \left(A^{-1}H(\vec{u}^{n})\right)_f - (A^{-1})_f \nabla p^n_f.(11)

This correction applies the continuity constraint in an iterative manner, while ensuring the lack of numerical pressure checker-boarding phenomena.

The next guess for the velocity, however, does not necessarily respect the momentum equation. Therefore, the momentum prediction and pressure correction steps need to be repeated until both the momentum and continuity equations are satisfied.

commentnote

The iterative process above is not stable if the full update is applied every time. This means that the variables need to be relaxed. Specifically, it is a common practice to relax the pressure when plugging it back to the gradient term in the momentum predictor:

pn+1,=pn+λp(pn+1pn),p^{n+1,*} = p^n + \lambda_p (p^{n+1}-p^n),(12)

where pn+1,p^{n+1,*} is the relaxed field and λp(0,1]\lambda_p \in (0,1] is the corresponding relaxation parameter.

commentnote

To help the solution process of the linear solver, we add options to ensure diagonal dominance through the relaxation of equations. This is done using the method mentioned in Juretic (2005), meaning that a numerical correction is added to the diagonal of the system matrix and the right hand side. This is especially useful for advection-dominated systems.

commentnote

Currently, this solver only respects the following execute_on flags: INITAL, TIMESTEP_BEGIN, and FINAL, other flags are ignored. MultiApps and the corresponding MultiappTransfers are executed at FINAL only.

Example Input Syntax

The setup of a problem with the segregated solver in MOOSE is slightly different compared to conventional monolithic solvers. In this section, we highlight the main differences. For setting up a 2D simulation with the SIMPLE algorithm, we need three linear systems in MOOSE: one for each momentum component and another for the pressure. The different systems can be created within the Problem block:

[Problem]
  linear_sys_names = 'u_system v_system pressure_system'
  previous_nl_solution_required = true
[]
(contrib/moose/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/2d-velocity-pressure.i)

It is visible that we requested that MOOSE keeps previous solution iterates as well. This is necessary to facilitate the relaxation processes mentioned in the overview. Next, we create linear FV variables and assign them to the given systems.

[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
  []
[]
(contrib/moose/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/2d-velocity-pressure.i)

The kernels are then created within the LinearFVKernels block. The fundamental terms that contribute to the face fluxes in the momentum equation (stress and advection terms) are lumped into one kernel. Furthermore, instead of adding contribution from the continuity equation, we build an anisotropic diffusion (Poisson) equation for pressure:

[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
  []
[]
(contrib/moose/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/2d-velocity-pressure.i)

By default, the coupling fields corresponding to A1HA^{-1}H and A1A^{-1} are handled by functor called HbyA and Ainv, respectively. These fields are generated by RhieChowMassFlux under the hood. This means that we need to add the user object responsible for generating these fields:

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
  []
[]
(contrib/moose/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/2d-velocity-pressure.i)

As a last step, we add the SIMPLE executioner:

[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 = 50
  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
[]
(contrib/moose/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/2d-velocity-pressure.i)

Input Parameters

  • momentum_systemsThe solver system(s) for the momentum equation(s).

    C++ Type:std::vector<SolverSystemName>

    Unit:(no unit assumed)

    Controllable:No

    Description:The solver system(s) for the momentum equation(s).

  • pressure_systemThe solver system for the pressure equation.

    C++ Type:SolverSystemName

    Unit:(no unit assumed)

    Controllable:No

    Description:The solver system for the pressure equation.

  • rhie_chow_user_objectThe rhie-chow user-object

    C++ Type:UserObjectName

    Unit:(no unit assumed)

    Controllable:No

    Description:The rhie-chow user-object

Required Parameters

  • print_fieldsFalseUse this to print the coupling and solution fields and matrices throughout the iteration.

    Default:False

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:Use this to print the coupling and solution fields and matrices throughout the iteration.

  • turbulence_field_min_limitThe lower limit imposed on turbulent quantities. The recommended value for robustness is 1e-8.

    C++ Type:std::vector<double>

    Unit:(no unit assumed)

    Controllable:No

    Description:The lower limit imposed on turbulent quantities. The recommended value for robustness is 1e-8.

  • verboseFalseSet to true to print additional information

    Default:False

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:Set to true to print additional information

Optional Parameters

  • control_tagsAdds user-defined labels for accessing object parameters via control logic.

    C++ Type:std::vector<std::string>

    Unit:(no unit assumed)

    Controllable:No

    Description:Adds user-defined labels for accessing object parameters via control logic.

  • enableTrueSet the enabled status of the MooseObject.

    Default:True

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:Set the enabled status of the MooseObject.

  • outputsVector of output names where you would like to restrict the output of variables(s) associated with this object

    C++ Type:std::vector<OutputName>

    Unit:(no unit assumed)

    Controllable:No

    Description:Vector of output names where you would like to restrict the output of variables(s) associated with this object

Advanced Parameters

  • momentum_absolute_tolerance1e-05The absolute tolerance on the normalized residual of the momentum equation.

    Default:1e-05

    C++ Type:double

    Unit:(no unit assumed)

    Controllable:No

    Description:The absolute tolerance on the normalized residual of the momentum equation.

  • num_iterations1000The number of momentum-pressure-(other fields) iterations needed.

    Default:1000

    C++ Type:unsigned int

    Unit:(no unit assumed)

    Controllable:No

    Description:The number of momentum-pressure-(other fields) iterations needed.

  • pressure_absolute_tolerance1e-05The absolute tolerance on the normalized residual of the pressure equation.

    Default:1e-05

    C++ Type:double

    Unit:(no unit assumed)

    Controllable:No

    Description:The absolute tolerance on the normalized residual of the pressure equation.

Iteration Control Parameters

  • momentum_equation_relaxation1The relaxation which should be used for the momentum equation. (=1 for no relaxation, diagonal dominance will still be enforced)

    Default:1

    C++ Type:double

    Unit:(no unit assumed)

    Controllable:No

    Description:The relaxation which should be used for the momentum equation. (=1 for no relaxation, diagonal dominance will still be enforced)

  • pressure_variable_relaxation1The relaxation which should be used for the pressure variable (=1 for no relaxation).

    Default:1

    C++ Type:double

    Unit:(no unit assumed)

    Controllable:No

    Description:The relaxation which should be used for the pressure variable (=1 for no relaxation).

Relaxation Parameters

  • momentum_l_abs_tol1e-50The absolute tolerance on the normalized residual in the linear solver of the momentum equation.

    Default:1e-50

    C++ Type:double

    Unit:(no unit assumed)

    Controllable:No

    Description:The absolute tolerance on the normalized residual in the linear solver of the momentum equation.

  • momentum_l_max_its10000The maximum allowed iterations in the linear solver of the momentum equation.

    Default:10000

    C++ Type:unsigned int

    Unit:(no unit assumed)

    Controllable:No

    Description:The maximum allowed iterations in the linear solver of the momentum equation.

  • momentum_l_tol1e-05The relative tolerance on the normalized residual in the linear solver of the momentum equation.

    Default:1e-05

    C++ Type:double

    Unit:(no unit assumed)

    Controllable:No

    Description:The relative tolerance on the normalized residual in the linear solver of the momentum equation.

  • pressure_l_abs_tol1e-10The absolute tolerance on the normalized residual in the linear solver of the pressure equation.

    Default:1e-10

    C++ Type:double

    Unit:(no unit assumed)

    Controllable:No

    Description:The absolute tolerance on the normalized residual in the linear solver of the pressure equation.

  • pressure_l_max_its10000The maximum allowed iterations in the linear solver of the pressure equation.

    Default:10000

    C++ Type:unsigned int

    Unit:(no unit assumed)

    Controllable:No

    Description:The maximum allowed iterations in the linear solver of the pressure equation.

  • pressure_l_tol1e-05The relative tolerance on the normalized residual in the linear solver of the pressure equation.

    Default:1e-05

    C++ Type:double

    Unit:(no unit assumed)

    Controllable:No

    Description:The relative tolerance on the normalized residual in the linear solver of the pressure equation.

Linear Iteration Control Parameters

  • momentum_petsc_optionsSingleton PETSc options for the momentum equation

    C++ Type:MultiMooseEnum

    Unit:(no unit assumed)

    Options:-dm_moose_print_embedding, -dm_view, -ksp_converged_reason, -ksp_gmres_modifiedgramschmidt, -ksp_monitor, -ksp_monitor_snes_lg-snes_ksp_ew, -ksp_snes_ew, -snes_converged_reason, -snes_ksp, -snes_ksp_ew, -snes_linesearch_monitor, -snes_mf, -snes_mf_operator, -snes_monitor, -snes_test_display, -snes_view

    Controllable:No

    Description:Singleton PETSc options for the momentum equation

  • momentum_petsc_options_inameNames of PETSc name/value pairs for the momentum equation

    C++ Type:MultiMooseEnum

    Unit:(no unit assumed)

    Options:-ksp_atol, -ksp_gmres_restart, -ksp_max_it, -ksp_pc_side, -ksp_rtol, -ksp_type, -mat_fd_coloring_err, -mat_fd_type, -mat_mffd_type, -pc_asm_overlap, -pc_factor_levels, -pc_factor_mat_ordering_type, -pc_hypre_boomeramg_grid_sweeps_all, -pc_hypre_boomeramg_max_iter, -pc_hypre_boomeramg_strong_threshold, -pc_hypre_type, -pc_type, -snes_atol, -snes_linesearch_type, -snes_ls, -snes_max_it, -snes_rtol, -snes_divergence_tolerance, -snes_type, -sub_ksp_type, -sub_pc_type

    Controllable:No

    Description:Names of PETSc name/value pairs for the momentum equation

  • momentum_petsc_options_valueValues of PETSc name/value pairs (must correspond with "petsc_options_iname" for the momentum equation

    C++ Type:std::vector<std::string>

    Unit:(no unit assumed)

    Controllable:No

    Description:Values of PETSc name/value pairs (must correspond with "petsc_options_iname" for the momentum equation

  • pressure_petsc_optionsSingleton PETSc options for the pressure equation

    C++ Type:MultiMooseEnum

    Unit:(no unit assumed)

    Options:-dm_moose_print_embedding, -dm_view, -ksp_converged_reason, -ksp_gmres_modifiedgramschmidt, -ksp_monitor, -ksp_monitor_snes_lg-snes_ksp_ew, -ksp_snes_ew, -snes_converged_reason, -snes_ksp, -snes_ksp_ew, -snes_linesearch_monitor, -snes_mf, -snes_mf_operator, -snes_monitor, -snes_test_display, -snes_view

    Controllable:No

    Description:Singleton PETSc options for the pressure equation

  • pressure_petsc_options_inameNames of PETSc name/value pairs for the pressure equation

    C++ Type:MultiMooseEnum

    Unit:(no unit assumed)

    Options:-ksp_atol, -ksp_gmres_restart, -ksp_max_it, -ksp_pc_side, -ksp_rtol, -ksp_type, -mat_fd_coloring_err, -mat_fd_type, -mat_mffd_type, -pc_asm_overlap, -pc_factor_levels, -pc_factor_mat_ordering_type, -pc_hypre_boomeramg_grid_sweeps_all, -pc_hypre_boomeramg_max_iter, -pc_hypre_boomeramg_strong_threshold, -pc_hypre_type, -pc_type, -snes_atol, -snes_linesearch_type, -snes_ls, -snes_max_it, -snes_rtol, -snes_divergence_tolerance, -snes_type, -sub_ksp_type, -sub_pc_type

    Controllable:No

    Description:Names of PETSc name/value pairs for the pressure equation

  • pressure_petsc_options_valueValues of PETSc name/value pairs (must correspond with "petsc_options_iname" for the pressure equation

    C++ Type:std::vector<std::string>

    Unit:(no unit assumed)

    Controllable:No

    Description:Values of PETSc name/value pairs (must correspond with "petsc_options_iname" for the pressure equation

Petsc Control Parameters

  • pin_pressureFalseIf the pressure field needs to be pinned at a point.

    Default:False

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:If the pressure field needs to be pinned at a point.

  • pressure_pin_pointThe point where the pressure needs to be pinned.

    C++ Type:libMesh::Point

    Unit:(no unit assumed)

    Controllable:No

    Description:The point where the pressure needs to be pinned.

  • pressure_pin_value0The value which needs to be enforced for the pressure.

    Default:0

    C++ Type:double

    Unit:(no unit assumed)

    Controllable:No

    Description:The value which needs to be enforced for the pressure.

Pressure Pin Parameters

    Restart Parameters

    References

    1. Hrvoje Jasak. Error analysis and estimation for the finite volume method with applications to fluid flows. PhD thesis, Imperial College London (University of London), 1996.[BibTeX]
    2. Franjo Juretic. Error analysis in finite volume CFD. PhD thesis, Imperial College London (University of London), 2005.[BibTeX]
    3. Fadl Moukalled, L Mangani, Marwan Darwish, and others. The finite volume method in computational fluid dynamics. Volume 6. Springer, 2016.[BibTeX]
    4. Suhas V Patankar and D Brian Spalding. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. In Numerical prediction of flow, heat transfer, turbulence and combustion, pages 54–73. Elsevier, 1983.[BibTeX]
    5. Chae M Rhie and Wei-Liang Chow. Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA journal, 21(11):1525–1532, 1983.[BibTeX]