8. PIMPLE Predictor–Corrector Method#

PIMPLE combines a SIMPLE-family nonlinear outer iteration with PISO-style pressure–momentum correctors. The outer iteration refreshes the nonlinear transport coefficients and coupled physical state. The corrector stages then reduce the remaining continuity and momentum-splitting error without repeating every implicit transport solve. This arrangement is especially useful for transient calculations in which one pressure correction per nonlinear iteration is unnecessarily expensive, while a purely PISO-style sequence is not sufficiently robust for the selected time step or mesh.

The Stream formulation was developed from the operator-splitting construction used by PISO [Iss86, TW04], combined with a time-consistent momentum-interpolation treatment [Pas11]. Its original development and validation are documented in the project final report [TW21].

This chapter uses the finite-volume and pressure-correction machinery already developed in:

Chapter 4 for assembled transport rows;
Chapter 5 for the cell-to-face momentum response; and
Chapter 6 for the compressible pressure-correction equation.

PIMPLE stages these components as an implicit predictor followed by repeated pressure–momentum correctors inside a nonlinear outer iteration.

8.1. Iteration Hierarchy#

Let \(n\) denote the accepted physical-time level, \(k\) the nonlinear outer iteration, and \(m\) the PISO corrector. An iteration \(\ell\) performed by an individual linear solver acts only within one equation solve and does not advance the PIMPLE state. Thus,

Index	Role
\(n\)	Accepted physical-time level
\(k\)	SIMPLE-family nonlinear outer iteration
\(m\)	PISO pressure–momentum corrector
\(\ell\)	Iteration internal to one linear-system solve

chapter_2/media/pimple_module_drawing.svg — Fig. 8.1.1 Predictor and corrector nesting for PISO and PIMPLE. PISO advances through the predictor and correctors at each time step. PIMPLE repeats that sequence through the outer nonlinear iteration \(k\) before accepting the new time level.#

For one physical time step, the state sequence is

(8.1.1)#\[\boldsymbol{Q}^{n} \longrightarrow \boldsymbol{Q}^{n,k,0} \longrightarrow \boldsymbol{Q}^{n,k,1} \longrightarrow\cdots\longrightarrow \boldsymbol{Q}^{n,k,N_c} \longrightarrow \boldsymbol{Q}^{n,k+1,0},\]

where \(\boldsymbol{Q}^{n,k,0}\) is the state after the implicit predictor and its first pressure correction, and \(N_c\) is the number of PISO correctors. Once the outer convergence or iteration criterion is satisfied, \(\boldsymbol{Q}^{n+1}\) is accepted from the last outer state.

Three Different Iterations

An outer iteration \(k\) repeats the implicit predictor equations and refreshes their nonlinear coefficients. A corrector \(m\) updates pressure, velocity, face flux, and selected auxiliary equations without another implicit momentum solve. A linear-solver iteration \(\ell\) reduces the algebraic residual of one equation without advancing either of the enclosing algorithmic loops. These iterations solve different problems and are not interchangeable.

8.2. Implicit Predictor#

At outer iteration \(k\), freeze the current nonlinear coefficients and write the time-discrete momentum row for cell \(P\) as

(8.2.1)#\[{}^{(u)}\!a_P^{k}\vec{u}_P^{*,k} = \sum_{N\in\mathcal{N}(P)} {}^{(u)}\!a_{PN}^{k}\vec{u}_N^{*,k} +{}^{(u)}\!\vec{b}_P^{k} +\vec{b}_{P,t}^{n} -\mathcal{G}_P\!\left(p^{k}\right).\]

Here \({}^{(u)}\!a_P^{k}\) contains the spatial and implicit temporal diagonal, \({}^{(u)}\!\vec{b}_P^{k}\) contains explicit spatial and physical-source terms, and \(\vec{b}_{P,t}^{n}\) contains known time-history terms. The discrete pressure force is

(8.2.2)#\[\mathcal{G}_P(p) = \sum_{f\in\mathcal{F}(P)} p_f\vec{A}_f \simeq \Omega_P\nabla p_P.\]

Define the non-pressure momentum source

(8.2.3)#\[\vec{H}_P^{k}(\vec{u}) = \sum_N {}^{(u)}\!a_{PN}^{k}\vec{u}_N +{}^{(u)}\!\vec{b}_P^{k} +\vec{b}_{P,t}^{n}.\]

The predictor row then exposes the two pieces needed by pressure correction:

(8.2.4)#\[\vec{u}_P^{*,k} = \frac{\vec{H}_P^{k}(\vec{u}^{*,k})}{{}^{(u)}\!a_P^{k}} - \frac{\mathcal{G}_P(p^k)}{{}^{(u)}\!a_P^{k}}.\]

The momentum equation is solved implicitly at this stage. The selected thermal, turbulence, species, phase, and other active predictor equations are likewise advanced through their normal outer-iteration paths before the predictor state is accepted.

8.2.1. Temporal Diagonal Split#

The face flux must retain the same temporal response as the cell momentum equation. Split the momentum diagonal into non-temporal and leading temporal parts,

(8.2.1.1)#\[{}^{(u)}\!a_P^{k} = {}^{(u)}\!\widetilde a_P^{k}+c_{0,P}^{k} = {}^{(u)}\!\widetilde a_P^{k}\left(1+\delta_P^k\right), \qquad \delta_P^k=\frac{c_{0,P}^{k}}{{}^{(u)}\!\widetilde a_P^{k}}.\]

This is the same split derived in Eqs. (5.5.2)–(5.5.5). It prevents the converged face flux from acquiring an artificial dependence on the time step through an unmatched temporal diagonal.

8.2.2. Predictor Face Flux#

For face \(f\) between cells \(P\) and \(N\), define the volumetric and mass fluxes

(8.2.2.1)#\[\mathcal{V}_f=\vec{u}_f\cdot\vec{A}_f, \qquad \dot m_f=\rho_f\left(\mathcal{V}_f-\mathcal{V}_{g,f}\right),\]

where \(\mathcal{V}_{g,f}\) is the grid volumetric flux. Momentum-weighted interpolation converts Eq. (8.2.4) into a pressure-sensitive face flux. Its predictor form can be written

(8.2.2.2)#\[\dot m_f^{*,k} = \widehat{\dot m}_f^{k} -d_f^{k}\left(p_N^k-p_P^k\right), \qquad d_f^{k}>0,\]

where \(\widehat{\dot m}_f^{k}\) contains the interpolated non-pressure momentum response, temporal history, grid motion, and explicitly treated geometry terms. The coefficient \(d_f^{k}\) is the derivative of the face mass flux with respect to the adjacent-cell pressure difference. The same coefficient must be used when the pressure equation corrects that flux.

8.2.3. Predictor Pressure Correction#

The predicted flux generally does not satisfy continuity. Define its cell residual by

(8.2.3.1)#\[R_P^{*,k} = \mathcal{D}_t\!\left(\rho_P^{k}\Omega_P\right) +\sum_{f\in\mathcal{F}(P)}\dot m_f^{*,k},\]

where \(\mathcal{D}_t\) denotes the selected discrete time derivative, including the accepted time-level history. Introduce the predictor pressure correction \(p_P^{\prime,*,k}\). Its direct face response is

(8.2.3.2)#\[\dot m_f^{\prime,*,k} = -d_f^{k}\left(p_N^{\prime,*,k}-p_P^{\prime,*,k}\right) +\dot m_{f,\rho}^{\prime,*,k}.\]

For compressible flow, the acoustic tangent supplies

(8.2.3.3)#\[\rho_P^{\prime,*,k} \approx \left(\frac{\partial\rho}{\partial p}\right)_{s,Y}^{k} p_P^{\prime,*,k} = \frac{p_P^{\prime,*,k}}{a_{s,P}^{2}},\]

and \(\dot m_{f,\rho}^{\prime,*,k}\) carries the upwind face density correction with the provisional relative volumetric flux. Substitution into continuity gives the standard row

(8.2.3.4)#\[{}^{(p')}\!a_P^{*,k}p_P^{\prime,*,k} = \sum_N {}^{(p')}\!a_{PN}^{*,k}p_N^{\prime,*,k} +{}^{(p')}\!b_P^{*,k}.\]

The detailed acoustic, boundary, and face-coefficient contributions are those derived in Chapter 6.

After solving Eq. (8.2.3.4), the predictor state that enters the PISO loop is

(8.2.3.5)#\[\begin{split}p_P^{k,0} &=p_P^k+\alpha_p p_P^{\prime,*,k},\\ \vec{u}_P^{k,0} &=\vec{u}_P^{*,k}+\vec{u}_P^{\prime,*,k},\\ \dot m_f^{k,0} &=\dot m_f^{*,k}+\dot m_f^{\prime,*,k}.\end{split}\]

The outer predictor may use the SIMPLE or SIMPLEC response described in Chapter 7. Each subsequent PISO corrector uses the unrelaxed momentum response. It updates the state produced by the predictor without repeating the under-relaxed outer iteration.

Predictor-Stage Result

The implicit transport solves, pressure correction, flux correction, and EOS closure together produce the complete state

(8.2.3.6)#\[\boldsymbol{Q}^{n,k,0} = \left\{ p^{k,0},\vec{u}^{\,k,0},\mathcal{V}_f^{k,0}, \Theta^{k,0},\boldsymbol{Y}^{k,0},\rho^{k,0},\ldots \right\}.\]

This state supplies the incoming fields and coefficients for the first PISO corrector.

8.3. PISO Corrector Derivation#

The first corrector is initialized from the complete predictor-stage state:

(8.3.1)#\[p^0=p^{k,0}, \qquad \vec{u}^{\,0}=\vec{u}^{\,k,0}, \qquad \mathcal{V}_f^0=\mathcal{V}_f^{k,0}, \qquad \boldsymbol{Q}^0=\boldsymbol{Q}^{n,k,0}.\]

Consider corrector \(m\), with the incoming fields denoted by superscript \(m\). The latest face flux and thermodynamic state determine a newly assembled momentum diagonal, while the incoming velocity supplies the neighbor values in the corresponding non-pressure momentum contribution. Denote saved incoming quantities by superscript \(o\) and the newly assembled corrector quantities by superscript \(c\):

(8.3.2)#\[\begin{split}{}^{(u)}\!a_P^{o}\vec{u}_P^{m} &= \vec{H}_P^{o}-\mathcal{G}_P(p^m),\\ {}^{(u)}\!a_P^{c}\vec{u}_P^{m+1} &= \vec{H}_P^{c}-\mathcal{G}_P(p^{m+1}).\end{split}\]

The corrector quantities have the schematic dependencies

(8.3.3)#\[\begin{split}{}^{(u)}\!a_P^{c} &= \mathcal{A}_P\!\left(\mathcal{V}^{m},\boldsymbol{Q}^{m}\right),\\ \vec{H}_P^{c} &= \sum_N {}^{(u)}\!a_{PN}^{c}\vec{u}_N^{m} +{}^{(u)}\!\vec{b}_P^{c} +\vec{b}_{P,t}^{n}.\end{split}\]

Here \(\mathcal{A}_P\) denotes assembly of the momentum diagonal from the latest face flux, properties, and source linearizations. Thus, the corrector refreshes the momentum response using the latest state but remains explicit in the velocity update.

With

(8.3.4)#\[p^{m+1}=p^m+p^{\prime,m},\]

subtract the first row in Eq. (8.3.2) from the second:

(8.3.5)#\[{}^{(u)}\!a_P^{c}\vec{u}_P^{m+1} -{}^{(u)}\!a_P^{o}\vec{u}_P^{m} = \left(\vec{H}_P^{c}-\vec{H}_P^{o}\right) -\mathcal{G}_P(p^{\prime,m}).\]

Solving for the corrected velocity gives

(8.3.6)#\[\vec{u}_P^{m+1} = \frac{{}^{(u)}\!a_P^{o}}{{}^{(u)}\!a_P^{c}}\vec{u}_P^{m} +\frac{\vec{H}_P^{c}-\vec{H}_P^{o}}{{}^{(u)}\!a_P^{c}} -\frac{\mathcal{G}_P(p^{\prime,m})}{{}^{(u)}\!a_P^{c}}.\]

Equation (8.3.6) separates three effects: the change in the diagonal, the change in the explicit momentum source, and the new pressure correction.

8.3.1. Timestep-Consistent Hatted Velocity#

Apply the temporal split of Eq. (8.2.1.1) to both corrector states and define

(8.3.1.1)#\[\widehat{\vec{u}}_P^{m} = \frac{\vec{H}_P^{c}-\vec{H}_P^{o}} {{}^{(u)}\!\widetilde a_P^{c}}.\]

Equation (8.3.6) becomes

(8.3.1.2)#\[\vec{u}_P^{m+1} = \frac{1+\delta_P^{o}}{1+\delta_P^{c}} \frac{{}^{(u)}\!\widetilde a_P^{o}}{{}^{(u)}\!\widetilde a_P^{c}} \vec{u}_P^{m} +\frac{\widehat{\vec{u}}_P^{m}}{1+\delta_P^{c}} -\frac{\mathcal{G}_P(p^{\prime,m})} {{}^{(u)}\!\widetilde a_P^{c}(1+\delta_P^{c})}.\]

The temporal factor is now explicit and can be carried consistently to the face. For an interior face between \(P\) and \(N\), let \(w_P+w_N=1\) denote the cell-to-face interpolation weights. Interpolate the old and current temporal ratios, the non-temporal diagonal ratio, and the hatted cell velocity as

(8.3.1.3)#\[\begin{split}\delta_f^{o}=w_P\delta_P^o+w_N\delta_N^o, \qquad \delta_f^{c}=w_P\delta_P^c+w_N\delta_N^c, \\ r_{a,f} = w_P\frac{{}^{(u)}\!\widetilde a_P^{o}}{{}^{(u)}\!\widetilde a_P^{c}} +w_N\frac{{}^{(u)}\!\widetilde a_N^{o}}{{}^{(u)}\!\widetilde a_N^{c}}, \qquad \widehat{\vec{u}}_f^{m} = w_P\widehat{\vec{u}}_P^{m} +w_N\widehat{\vec{u}}_N^{m}.\end{split}\]

The hatted corrector volume flux is then

(8.3.1.4)#\[\widehat{\mathcal{V}}_f^{m} = \frac{(1+\delta_f^o)r_{a,f}}{1+\delta_f^c} \mathcal{V}_f^{m} + \frac{\widehat{\vec{u}}_f^{m}\cdot\vec{A}_f} {1+\delta_f^c}.\]

This is the principal timestep-consistent PIMPLE flux relation. It advances the old face flux with the same diagonal-ratio and temporal factors used by the cell velocity, then adds the interpolated change in the non-pressure momentum source. The pressure response is applied next.

8.3.2. Corrector Pressure Response#

For the time-consistent form, the interior-face mass-flux coefficient is

(8.3.2.1)#\[d_f^{m} = \rho_f A_f \frac{ \mathcal{I}_f\!\left(\dfrac{\Omega}{{}^{(u)}\!\widetilde a^{c}}\right) }{ (1+\delta_f^c)\ell_f },\]

where \(\ell_f\) is the projected cell-center spacing defined in Eq. (5.3.3.4). No momentum under-relaxation factor appears in Eq. (8.3.2.1): the PISO corrector uses the raw momentum response.

On a non-orthogonal mesh, write the signed lagged geometric correction from the preceding corrector as \(\mathcal{E}_f(p^{\prime,m-1})\). The corrector volume-flux increment has the generic form

(8.3.2.2)#\[\mathcal{V}_f^{\prime,m} = -\frac{d_f^{m}}{\rho_f} \left(p_N^{\prime,m}-p_P^{\prime,m}\right) +\frac{d_f^{m}}{\rho_f} \mathcal{E}_f\!\left(p^{\prime,m-1}\right).\]

The lagged correction is zero for the first corrector. Each subsequent corrector recovers the geometric contribution left explicit by the preceding stage. The corresponding linearized mass-flux response is

(8.3.2.3)#\[\dot m_f^{\prime,m} = \rho_f\mathcal{V}_f^{\prime,m} +\left(\widehat{\mathcal{V}}_f^{m}-\mathcal{V}_{g,f}\right) \rho_{\mathrm{up}}^{\prime,m}.\]

The second term is present for compressible flow and represents the upwind acoustic density response. The acoustic term enters the mass-balance linearization, whereas the volume-flux correction contains the pressure-gradient and lagged geometric responses in Eq. (8.3.2.2).

Substituting Eq. (8.3.2.3) and the acoustic unsteady-density response into continuity produces the corrector pressure row.

Final Corrector Pressure Row

For cell \(P\) at corrector \(m\),

(8.3.2.4)#\[{}^{(p')}\!a_P^{m}p_P^{\prime,m} = \sum_N {}^{(p')}\!a_{PN}^{m}p_N^{\prime,m} +{}^{(p')}\!b_P^{m},\]

with the coefficient groups

(8.3.2.5)#\[\begin{split}{}^{(p')}\!a_P^{m} &= {}^{(p')}\!a_{P,t} +\sum_f d_f^{m} +{}^{(p')}\!a_{P,\rho} +{}^{(p')}\!a_{P,\mathrm{bc}},\\ {}^{(p')}\!a_{PN}^{m} &= d_f^{m}+{}^{(p')}\!a_{PN,\rho},\\ {}^{(p')}\!b_P^{m} &= -\widehat R_P^{m} +{}^{(p')}\!b_{P,t} +{}^{(p')}\!b_{P,\mathrm{nonorth}}^{m-1} +{}^{(p')}\!b_{P,\mathrm{bc}}.\end{split}\]

Here \(\widehat R_P^{m}\) is the continuity residual formed from the hatted corrector flux, \({}^{(p')}\!a_{P,t}\) is the implicit acoustic time contribution, and \({}^{(p')}\!a_{P,\rho}\) and \({}^{(p')}\!a_{PN,\rho}\) are the upwind compressible face contributions. Interior pressure-response faces contribute \(d_f^{m}\) to the diagonal and to the corresponding neighbor coefficient. A boundary whose normal flux responds to pressure contributes the corresponding one-sided diagonal term; a prescribed-flux boundary does not acquire that response.

8.3.3. Completing One Corrector#

Solving Eq. (8.3.2.4) completes corrector \(m\):

(8.3.3.1)#\[\begin{split}p_P^{m+1} &=p_P^{m}+p_P^{\prime,m},\\ \vec{u}_P^{m+1} &= \frac{1+\delta_P^{o}}{1+\delta_P^{c}} \frac{{}^{(u)}\!\widetilde a_P^{o}}{{}^{(u)}\!\widetilde a_P^{c}} \vec{u}_P^{m} +\frac{\widehat{\vec{u}}_P^{m}}{1+\delta_P^{c}} -\frac{1}{{}^{(u)}\!a_P^{c}} \mathcal{G}_P\!\left(p^{\prime,m}\right),\\ \mathcal{V}_f^{m+1} &=\widehat{\mathcal{V}}_f^{m} +\mathcal{V}_f^{\prime,m}.\end{split}\]

The first two velocity terms are the hatted corrector update derived in Eq. (8.3.1.2); the final term is its pressure response. Pressure, cell velocity, and face volume flux therefore advance as one corrected state from \(m\) to \(m+1\).

Conservation Pairing

The pressure row and the face-flux update must use the same \(d_f^{m}\), acoustic density response, orientation, and boundary participation. The cell velocity correction is a pressure-force response, while the face-flux correction is the quantity that enforces discrete continuity. Updating one without the matching other two breaks the pressure-based coupling.

8.4. Thermal State and EOS Closure#

After pressure, velocity, and face flux have been corrected, an explicit thermal corrector updates the transported thermal variable \(\Theta\) by evaluating its assembled scalar row with the incoming neighbor values:

(8.4.1)#\[\Theta_P^{m+1} = \frac{ {}^{(\Theta)}\!b_P^{m} +\sum_N {}^{(\Theta)}\!a_{PN}^{m}\Theta_N^m }{{}^{(\Theta)}\!a_P^{m}}.\]

The coefficients and source use the corrected mass flux together with the current pressure, velocity, diffusion, source, and time-history data. When a thermal correction is not part of the corrector stage, the incoming thermal state is retained until the next outer predictor.

For compressible flow, the complete density is then reconstructed from the EOS:

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

When total energy or total enthalpy is transported, the corrected velocity is used to remove the kinetic-energy contribution before the static thermal state is passed to the EOS. When temperature is transported, the corrected pressure and temperature close the EOS directly.

The acoustic estimate \(\rho'=p'/a_s^2\) belongs only to the linear pressure equation. Equation (8.4.2) supplies the accepted nonlinear density for the next corrector or outer iteration.

Corrector-Stage Result

At the end of corrector \(m\), pressure, cell velocity, face flux, transported auxiliary fields, and the thermodynamic state all carry the same stage index \(m+1\). The next corrector therefore begins from a complete state rather than from pressure and velocity alone.

8.5. Completing the Outer PIMPLE Iteration#

The full outer iteration can now be summarized as follows:

Assemble and implicitly solve the active predictor equations using the current outer state.
Construct the momentum-consistent predictor face flux.
Assemble and solve the predictor pressure-correction equation.
Correct pressure, cell velocity, and face flux; update the thermal state and rebuild the EOS state.
For \(m=0,\ldots,N_c-1\):
1. reassemble the corrector momentum diagonal and non-pressure source;
2. form the timestep-consistent hatted cell velocity and face flux;
3. solve the corrector pressure equation, including the lagged geometric contribution;
4. update selected explicit corrector equations; and
5. rebuild the thermodynamic state.
Promote the final corrector fields to outer iteration \(k+1\).
Repeat the outer loop until its convergence criterion or maximum iteration count is reached, then accept the new physical-time state.

The PISO correctors reduce the pressure–velocity splitting error inside one outer iteration. The outer loop still performs the nonlinear work: it refreshes transport properties, advective coefficients, source linearizations, and any other model couplings that the explicit correctors do not fully recover.

Stability and Time-Step Interpretation

PIMPLE can tolerate a larger time step than a single-pass PISO treatment when the additional outer reassembly is sufficient to control nonlinear and mesh coupling errors. This is not an unconditional Courant-number guarantee. Accuracy remains controlled by temporal resolution, and stability remains dependent on the mesh, physical models, equation linearizations, corrector count, and linear-solver convergence.

PIMPLE Algorithm Summary

Predictor: solve the nonlinear transport systems, construct a momentum-consistent face flux, solve one pressure-correction row, and rebuild the thermodynamic state.

Corrector: update the momentum response explicitly, solve the corrector pressure row, apply the same correction to pressure, velocity, and face flux, then update selected auxiliary equations and the EOS.

Outer advance: promote the complete corrected state and repeat the predictor only when another nonlinear iteration is required.

8.6. References#

[Iss86]

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

[Pas11]

A. Pascau. Cell face velocity alternatives in a structured colocated grid for the unsteady navier–stokes equations. International Journal for Numerical Methods in Fluids, 65:812–833, 2011. doi:10.1002/fld.2215.

[TW04]

S. Thakur and J. Wright. A multiblock operator-splitting algorithm for unsteady flows at all speeds in complex geometries. International Journal for Numerical Methods in Fluids, 46:383–413, 2004.

[TW21]

S. Thakur and J. Wright. A fast and robust pimple-based algebraic vof method for two-phase compressible flows. Technical Report NASA SBIR Phase I Final Report, Contract 80NSSC20C0536, Streamline Numerics, Inc., March 2021.