7. Pressure–Velocity Coupling Algorithms

Chapter 4, Chapter 5, and Chapter 6 developed the common components of a collocated, segregated pressure-based method: finite-volume transport assembly, momentum-weighted face fluxes, and the pressure-correction equation. This chapter now asks the narrower comparative question: which momentum-correction terms are approximated, which pressure equation is solved, and how often is the correction repeated in SIMPLE, SIMPLEC, SIMPLER, PISO, and PIMPLE?

The discussion begins with one common correction identity so that the named methods can be compared by their mathematical choices rather than by separate lists of procedural steps. General method theory is kept distinct from the short implementation note at the end of the chapter.

7.1. Common Momentum-Correction Identity

The five algorithms compared here begin from the same linearized momentum equation. With the coefficients and non-pressure source terms frozen at outer iteration \(k\), write the provisional momentum row for cell \(P\) as

(7.1.1)\[a_P^{\mathrm{pred}}\vec{u}_P^* = \sum_{N\in\mathcal{N}(P)} a_{P,N}\vec{u}_N^* + \vec{b}_{P,\mathrm{np}}^{\,\mathrm{pred}} - \Omega_P \nabla p_P^k .\]

Only momentum coefficients occur in this chapter, so the equation-family prescript \({}^{(u)}\) is suppressed. The neighbor coefficient \(a_{P,N}\) belongs to the row for \(P\); it is not the neighboring cell’s diagonal. The subscript \(\mathrm{np}\) identifies the nonpressure part of the RHS, with pressure displayed separately.

Here \(a_P^{\mathrm{pred}}\) is the diagonal actually used in the momentum predictor, including any algebraic velocity under-relaxation. If \(a_P\) denotes the unrelaxed diagonal and \(\alpha_v\) is the velocity-relaxation factor, then

(7.1.2)\[a_P^{\mathrm{pred}}=\frac{a_P}{\alpha_v},\]

and \(\vec{b}_{P,\mathrm{np}}^{\,\mathrm{pred}}\) includes the corresponding deferred relaxation source. This distinction matters because the same relaxed momentum operator must be used when the pressure response is derived.

The provisional velocity is a momentum solution at the current pressure; it is not expected to satisfy discrete continuity. For the purpose of deriving the unrelaxed momentum response, introduce

(7.1.3)\[p_P^{\mathrm{corr}}=p_P^k+p_P', \qquad \vec{u}_P^{\mathrm{corr}}=\vec{u}_P^*+\vec{u}_P'.\]

If the corrected fields were to satisfy the same frozen momentum row exactly, then

(7.1.4)\[a_P^{\mathrm{pred}} \left(\vec{u}_P^*+\vec{u}_P'\right) = \sum_N a_{P,N}\left(\vec{u}_N^*+\vec{u}_N'\right) + \vec{b}_{P,\mathrm{np}}^{\,\mathrm{pred}} - \Omega_P\nabla\left(p_P^k+p_P'\right).\]

Subtracting Eq. (7.1.1) gives the exact correction identity for that linearization:

(7.1.5)\[a_P^{\mathrm{pred}}\vec{u}_P' = \sum_N a_{P,N}\vec{u}_N' - \Omega_P\nabla p_P'.\]

The raw pressure correction is used in this response derivation. A SIMPLE-family outer iteration may subsequently under-relax the primary pressure update, \(p^{k+1}=p^k+\alpha_p p'\), without changing the velocity/flux response required by the pressure-correction equation.

This identity exposes the central approximation made by SIMPLE-family methods. The correction at \(P\) depends on the corrections at neighboring cells, so using the exact inverse would require another coupled vector solve inside the pressure-correction construction. The named methods differ in how they approximate or recover that missing coupling.

Assumptions in the Common Correction Equation

• The momentum coefficients and explicit sources are frozen during one pressure-correction stage.

• The pressure correction represents the pressure change relative to the pressure used in the momentum predictor.

• The velocity and face-flux response uses an approximation to the inverse of the frozen momentum operator.

• For compressible flow, the pressure equation retains an acoustic density response, while the complete density is rebuilt from the EOS after the corrected thermodynamic state is available.

7.1.1. SIMPLE Local-Diagonal Response

SIMPLE neglects the neighbor-correction sum in Eq. (7.1.5). The resulting local response is

(7.1.1.1)\[\vec{u}_P' \approx - D_P^{\mathrm{S}}\nabla p_P', \qquad D_P^{\mathrm{S}} = \frac{\Omega_P}{a_P^{\mathrm{pred}}} = \frac{\alpha_v\Omega_P}{a_P}.\]

The approximation is inexpensive and robust, but the omitted neighbor corrections must be recovered gradually through outer iterations. This is one reason SIMPLE commonly uses under-relaxation.

7.1.2. SIMPLEC Consistent Response

SIMPLEC begins from the same exact identity but approximates neighboring velocity corrections by the local correction over the compact stencil,

(7.1.2.1)\[\vec{u}_N'\approx\vec{u}_P'.\]

Equation (7.1.5) then becomes

(7.1.2.2)\[\left(a_P^{\mathrm{pred}}-\sum_N a_{P,N}\right)\vec{u}_P' \approx - \Omega_P\nabla p_P',\]

and hence

(7.1.2.3)\[\vec{u}_P' \approx - D_P^{\mathrm{C}}\nabla p_P', \qquad D_P^{\mathrm{C}} = \frac{\Omega_P} {a_P^{\mathrm{pred}}-\sum_N a_{P,N}} = \frac{\Omega_P} {a_P/\alpha_v-\sum_N a_{P,N}}.\]

For an unrelaxed momentum predictor, \(\alpha_v=1\), this reduces to the familiar denominator \(a_P-\sum_N a_{P,N}\). The neighbor-coefficient sum, rather than an unexplained empirical multiplier, is the formal SIMPLEC distinction. Transient terms, source linearization, non-orthogonal treatment, and the face interpolation used in an actual unstructured implementation modify how this effective response is evaluated, but not why it differs from SIMPLE.

7.1.3. From Cell Correction to Face-Flux Correction

Continuity is a balance of face fluxes, not cell velocities. The approximate cell response must therefore be converted into a pressure-sensitive face response. On a collocated mesh this is done by the momentum-weighted interpolation derived in Chapter 5. For an interior face it supplies the generic pressure-difference correction

(7.1.3.1)\[\dot m_f' = -d_f\left(p_N'-p_P'\right) + \dot m_{f,\rho}',\]

where \(d_f\) contains the selected SIMPLE or SIMPLEC momentum response and \(\dot m_{f,\rho}'\) is the additional pressure-induced density-flux term for compressible flow. Substitution into discrete continuity produces the pressure-correction system

(7.1.3.2)\[\mathbf{A}_{p'}p'=\mathbf{b}_{p'},\]

as derived in Chapter 6. Momentum interpolation and the named coupling algorithm therefore play different roles: the former constructs a compatible face flux on a collocated mesh; the latter decides how the pressure–velocity correction is approximated and repeated.

7.1.4. Corrected State and EOS Reconstruction

For SIMPLE or SIMPLEC, one outer correction has the schematic form

(7.1.4.1)\[\begin{split}p_P^{k+1} &= p_P^k+\alpha_p p_P', \\ \vec{u}_P^{k+1} &= \vec{u}_P^*-D_P\nabla p_P', \\ \dot m_f^{k+1} &= \hat m_f^*+\dot m_f',\end{split}\]

where \(D_P\) and the corresponding face coefficient \(d_f\) use the selected momentum-response approximation. Pressure relaxation has been shown explicitly; the precise relaxation conventions for velocity and flux depend on how the predictor coefficients were formed.

For compressible flow, the relation

(7.1.4.2)\[\rho'\approx\frac{p'}{a_s^2}\]

is used only inside the linear pressure equation. It defines the hatted pressure-equation estimate \(\widehat{\rho}^{\,k+1}=\rho^k+\rho'\), not the complete density update. After the corrected pressure and the updated thermal and composition variables are available, the next nonlinear density is rebuilt from the EOS:

(7.1.4.3)\[\rho_P^{k+1} = \mathcal{R}\!\left( p_P^{k+1}, \Theta_P^{k+1}, \boldsymbol{Y}_P^{k+1}, \ldots \right).\]

Corrected velocity also enters this reconstruction when total energy or total enthalpy must first be converted to its static form. This staging is common to the compressible variants even though their pressure update and number of correctors differ.

Common-Machinery Quick Reference

SIMPLE drops the neighbor-correction sum. SIMPLEC approximates it and replaces \(a_P^{\mathrm{pred}}\) by \(a_P^{\mathrm{pred}}-\sum_N a_{P,N}\) in the local response. Momentum-weighted interpolation turns that response into a face coefficient \(d_f\). Discrete continuity then determines \(p'\). The acoustic density correction belongs to that linear equation; the EOS supplies the complete next-iteration density.

7.2. Algorithm Families

The methods most commonly associated with pressure-based CFD differ mainly in three places:

1. how the momentum predictor is formed,

2. whether one solves for a pressure-like variable before the correction equation,

3. how many times the pressure-correction step is repeated within a time step.

At a high level, all of them reuse the same building blocks:

(7.2.1)\[\text{momentum predictor} \;\rightarrow\; \text{pressure or pressure-correction solve} \;\rightarrow\; \text{velocity/flux correction}.\]

The distinction lies in how tightly those blocks are coupled.

7.2.1. SIMPLE

SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) [Pat80] uses one momentum predictor and one pressure-correction equation per outer nonlinear iteration. In its standard form, the velocity-correction relation neglects the neighbor-correction summation, so the method takes the form

(7.2.1.1)\[A_u^k \vec{u}^{*} = \mathcal{H}_u^k - \Omega \nabla p^k\]

(7.2.1.2)\[A_p^k p' = b_p^k(\dot{m}^{*})\]

(7.2.1.3)\[p^{k+1} = p^k + \alpha_p p'\]

(7.2.1.4)\[\vec{u}^{\,k+1} = \vec{u}^{*} - D_u^k \nabla p'\]

Here \(D_u^k\) is the diagonal approximation to the velocity-correction operator; in a cell-local reading one may think of \(D_{u,P}^k \approx \Omega_P/a_P^{\mathrm{pred}} =\alpha_v\Omega_P/a_P\). The \(()\) notation in \(b_p^k(\dot{m}^{*})\) is only indicating dependence: the pressure-correction right-hand side is assembled from the provisional face-flux field \(\dot{m}_f^{*}\), i.e. from the current mass-balance residual of that provisional state. It is not a multiplication by \(\dot{m}^{*}\).

What makes SIMPLE distinctive is its simplicity and robustness: the pressure field enters through a single correction equation per outer iteration, but the method usually relies more strongly on under-relaxation than the other variants.

7.2.2. SIMPLEC

SIMPLEC (SIMPLE-Consistent) [VDR84] follows the same outer structure as SIMPLE, but uses a more consistent approximation for the neglected neighbor-velocity-correction terms. In this compact iteration-level view, the momentum predictor and pressure-correction solve retain the same form:

(7.2.2.1)\[A_u^k \vec{u}^{*} = \mathcal{H}_u^k - \Omega \nabla p^k\]

(7.2.2.2)\[A_p^k p' = b_p^k(\dot{m}^{*})\]

while the corrected fields are written as

(7.2.2.3)\[p^{k+1} = p^k + \alpha_p p'\]

(7.2.2.4)\[\vec{u}^{\,k+1} = \vec{u}^{*} - D_u^{\mathrm{C},k}\nabla p', \qquad D_{u,P}^{\mathrm{C},k} = \frac{\Omega_P} {a_P^{\mathrm{pred}}-\sum_N a_{P,N}}\]

Thus SIMPLEC is not a different predictor/corrector architecture. It uses the same sequence with the consistent neighbor-correction approximation derived in Eq. (7.1.2.3). The corresponding face coefficient in the pressure equation must use the same response. Because it recovers more of the omitted momentum coupling within one correction, SIMPLEC commonly requires less restrictive relaxation than SIMPLE, although the admissible choice remains problem dependent.

7.2.3. SIMPLER

SIMPLER (SIMPLE-Revised) [Pat80] introduces an additional pressure equation before the momentum predictor. The key object is the pressure-free or hatted velocity formed from the non-pressure part of the momentum solve,

(7.2.3.1)\[\hat{\vec{u}}^k = \left(A_u^k\right)^{-1}\mathcal{H}_u^k,\]

This inverse notation is only shorthand for “solve the momentum linear system with the non-pressure right-hand side.” No explicit dense matrix inverse is formed in the implementation. From this one forms a pressure equation for an improved pressure estimate \(\tilde{p}\):

(7.2.3.2)\[A_{\tilde p}^k \tilde{p} = b_{\tilde p}^k(\hat{\vec{u}}^k)\]

Here \(A_{\tilde p}^k\) is the operator for the preliminary SIMPLER pressure equation, and \(b_{\tilde p}^k(\hat{\vec{u}}^k)\) is the corresponding continuity-based right-hand side formed from the hatted velocity field.

The preliminary pressure is relaxed, when required, to give the pressure used by the momentum predictor:

(7.2.3.3)\[p_{\mathrm{pre}}^{k+1} = (1-\alpha_{\tilde p})p^k + \alpha_{\tilde p}\tilde p .\]

The momentum equations are then solved with this improved pressure,

(7.2.3.4)\[A_u^k \vec{u}^{*} = \mathcal{H}_u^k - \Omega \nabla p_{\mathrm{pre}}^{k+1}.\]

The remaining pressure-correction equation is used only to remove the residual continuity error:

(7.2.3.5)\[A_p^k p' = b_p^k(\dot{m}^{*})\]

(7.2.3.6)\[\vec{u}^{\,k+1} = \vec{u}^{*} - D_u^k \nabla p'\]

The face flux is corrected by the same \(p'\) response used in continuity. In the defining SIMPLER organization, however, \(p'\) is not added to the primary pressure as it is in SIMPLE. The primary pressure is the relaxed preliminary pressure \(p_{\mathrm{pre}}^{k+1}\); the subsequent pressure-correction solve is used to correct velocity and face flux. This is the exact answer to the common question “How is SIMPLER different from SIMPLE?”: SIMPLER obtains its pressure from an additional pressure equation before the momentum predictor, whereas SIMPLE obtains its next pressure by adding a correction after the predictor.

7.2.4. PISO

PISO (Pressure-Implicit with Splitting of Operators) is aimed primarily at transient calculations[Iss86]. Instead of relying on many outer under-relaxed iterations, it performs one momentum predictor and then applies multiple pressure-correction stages within the same time step:

(7.2.4.1)\[A_u^n \vec{u}^{(0)} = \mathcal{H}_u^n - \Omega \nabla p^n,\]

with the stage index initialized by \(p^{(0)} \equiv p^n\). For each corrector stage \(m=1,\ldots,N_c\), one then writes

(7.2.4.2)\[A_p^{(m)} p'^{(m)} = b_p^{(m)}\!\left(\dot{m}^{(m-1)}\right)\]

(7.2.4.3)\[p^{(m)} = p^{(m-1)} + p'^{(m)}\]

(7.2.4.4)\[\vec{u}^{(m)} = \vec{u}^{(m-1)} - D_u^{(m)} \nabla p'^{(m)}\]

The special feature of PISO is the repeated corrector loop inside the time step. Each new corrector uses the already corrected fluxes from the previous stage and includes deferred correction information left unresolved by the preceding approximate momentum response. Pressure, velocity, and face flux are therefore progressively corrected without reassembling every nonlinear transport coefficient as a full SIMPLE outer iteration would. In a compressible formulation, the density and thermodynamic state must also remain consistent with the corrected stage state; the exact placement of that reconstruction is part of the particular PISO implementation.

7.2.5. PIMPLE

PIMPLE combines a SIMPLE-like outer loop with PISO-like inner correctors. In a compact iteration-level summary, one outer iteration begins with the momentum predictor

(7.2.5.1)\[A_u^k \vec{u}^{*} = \mathcal{H}_u^k - \Omega \nabla p^k,\]

The predictor pressure correction then produces the state that initializes the inner-corrector loop,

(7.2.5.2)\[p^{k,0} = p^k + \alpha_p p'^{k,*}, \qquad \vec{u}^{\,k,0} = \vec{u}^{*}+\vec{u}'^{k,*}.\]

For each corrector stage \(m=0,\ldots,N_c-1\), one then writes

(7.2.5.3)\[A_p^{k,m} p'^{\,k,m} = b_p^{k,m}\!\left(\dot{m}^{\,k,m}\right)\]

(7.2.5.4)\[p^{k,m+1} = p^{k,m} + p'^{\,k,m}\]

(7.2.5.5)\[\vec{u}^{\,k,m+1} = \widehat{\vec{u}}^{\,k,m} - D_u^{k,m} \nabla p'^{\,k,m},\]

where \(\widehat{\vec{u}}^{\,k,m}\) includes the explicit change in the momentum row between correctors. The matching pressure-free face flux is used to form the pressure-equation residual.

After the inner correctors are completed, the outer iteration is updated by

(7.2.5.6)\[p^{k+1} \equiv p^{k,N_c}, \qquad \vec{u}^{\,k+1} \equiv \vec{u}^{\,k,N_c}.\]

The outer loop plays the role of SIMPLE by updating the nonlinear state, while the inner corrector loop plays the role of PISO by tightening continuity within the time step. This is why PIMPLE is often described as a hybrid of SIMPLE and PISO: it is designed to retain the robustness of outer iterations while reducing the number of such iterations required for transient problems.

The predictor, timestep-consistent corrector flux, compressible pressure row, thermodynamic closure, and Stream iteration hierarchy are derived in Chapter 8.

7.3. Comparison of the Methods

The common algebra makes the distinctions compact:

Method	Approximate momentum response	Source of primary pressure	Corrector structure	Typical relaxation role
SIMPLE	Local diagonal, \(D_P^{\mathrm{S}}=\Omega_P/a_P^{\mathrm{pred}}\)	\(p^{k+1}=p^k+\alpha_p p'\)	One pressure correction per outer iteration	Usually important for steady nonlinear convergence
SIMPLEC	Consistent local response, \(D_P^{\mathrm{C}}=\Omega_P/(a_P^{\mathrm{pred}}-\sum_N a_{P,N})\)	Same additive pressure update as SIMPLE	One pressure correction per outer iteration	Often less restrictive because the correction is stronger
SIMPLER	SIMPLE- or SIMPLEC-type response after a preliminary pressure solve	Relaxed solution of the preliminary pressure equation	Preliminary pressure solve, momentum predictor, then one pressure correction for velocity and flux	Applied to the preliminary pressure and, as needed, other blocks
PISO	Successive corrections recover deferred momentum coupling	Sum of staged pressure corrections	Multiple correctors within a time step	Designed to reduce dependence on outer under-relaxation
PIMPLE	SIMPLE-like outer response plus PISO-like staged corrections	Updated through the inner correctors inside each outer iteration	Outer nonlinear loop containing multiple correctors	Outer relaxation may remain useful for difficult nonlinear coupling

Algorithm-Family Quick Reference

SIMPLE and SIMPLEC differ in the inverse momentum response used by the correction. SIMPLER differs in where the primary pressure comes from. PISO differs in repeating pressure corrections before a full nonlinear reassembly. PIMPLE nests that repeated-corrector idea inside an outer SIMPLE-like nonlinear loop.

7.4. Relation to the Stream Source Organization

The main Stream pressure-based selector establishes the SIMPLE, SIMPLEC, and SIMPLER rule families. The SIMPLER path contains a separate preliminary pressure system and advances the relaxed result of that system as the primary pressure, while its pressure-correction result is used by the velocity and flux corrections. Repeated-corrector formulations are organized in the separate src/piso/ and src/pimple/ modules; the latter exposes an explicit corrector iteration level and a configurable number of corrector iterations.

Stream evaluates the same SIMPLEC response in a transient-aware multiplier form. In the representative decomposition \(a_P=a_{\mathrm{sp}}+a_t\), the quantity \(\delta\simeq a_t/a_{\mathrm{sp}}\) measures the temporal contribution relative to the non-temporal momentum diagonal; its exact value follows the selected time-integration weights. With \(\alpha_v\) denoting the momentum-relaxation factor, the multiplier relative to the SIMPLE response is

(7.4.1)\[\kappa = \frac{\alpha_v}{1-\alpha_v+\delta}, \qquad D_P^{\mathrm{C}} = (1+\kappa)D_P^{\mathrm{S}}.\]

For the simple decomposition \(\delta=a_t/a_{\mathrm{sp}}\), when the source-free interior continuity-subtracted row gives \(a_{\mathrm{sp}}\simeq\sum_N a_{P,N}\), this is algebraically equivalent to \(D_P^{\mathrm{C}}=\Omega_P/(a_P/\alpha_v-\sum_N a_{P,N})\). The cell correction uses the local ratio, while the compact face closure uses its face ratio. The resulting face coefficient must be the same in the pressure equation and the subsequent face-flux correction. The momentum-interpolation chapter derives this additional face approximation in Eq. (5.9.2.2); it is not generally an exact interpolation of the cell response.

These source landmarks identify the intended algorithmic roles; they do not replace case-specific schedule inspection when establishing which optional module and rule path a particular run assembled. Detailed fact and matrix-store names are collected in the solver framework reference.

7.5. References

[Iss86]

R.I. Issa. Solution of the implicitly discretised fluid flow equations by operator splitting. Journal of Computational Physics, 62:40–65, 1986.

[Pat80] (1,2)

S.V. Patankar. Numerical Heat Transfer and Fluid Flow. Hemisphere, 1980.

[VDR84]

J.P. Van Doormaal and G.D. Raithby. Enhancements of the simple method for predicting incompressible fluid flows. Numerical Heat Transfer, 7(2):147–163, 1984. URL: https://www.tandfonline.com/doi/abs/10.1080/01495728408961817, doi:10.1080/01495728408961817.