PISO

More information about Issa's PISO (Pressure-Implicit with Splitting of Operators) algorithm is given in References ,,,,. Rather than solve all of the coupled equations in a coupled or iterative sequential fashion, PISO splits the operators into an implicit predictor and multiple explicit corrector steps. This scheme is not thought of as iterative, and very few corrector steps are necessary to obtain desired accuracy. At each time step with buoyantBoussinesqPisoFoam, velocity and temperature are predicted, and then pressure and velocity are corrected. Temperature is not corrected for an unknown reason; however, Oliveira and Issa recommend temperature correction for incompressible flow using the Boussinesq approximation. It is possible that since this is not really compressible flow because the Boussinesq approximation is used, the maximum difference in temperature in a problem suitable for this solver is small and a corrector is not necessary.

The PISO algorithm can be understood, for simplicity, by considering a one-dimensional, inviscid flow along the $x$-direction in which gravity acts in that direction as well. Here $u$ could be replaced by $\bar{u}$ in RANS or LES with no consequence to understanding the PISO algorithm. The momentum equation simplifies to

$\frac{\partial u}{\partial t} + \frac{\partial}{\partial x}\left( u u \right) = -\frac{\partial p}{\partial x} + \rho_k g.$

(19)

1 Predictor

As an example, if Euler implicit time stepping is used with linear interpolation of values to the cell faces and linearization of the convective term by taking the convective velocity to be from the old time step $n$ (this treatment of the convective term is the same in as OpenFOAM according to Nilsson, then the discretized implicit velocity predictor form of Equation 19 is

$\left[ \frac{1}{\Delta t}+\left( \frac{u_{i+\frac{1}{2}}^n-u_{i-\frac{1}{2}}^n}{2 \Delta x} \right) \right] \Delta V u_i^{*} + \left( \frac{u_{i+\frac{1}{2}}^n}{2 \Delta x} \right) \Delta V u_{i+1}^{*} - \left( \frac{u_{i-\frac{1}{2}}^n}{2 \Delta x} \right) \Delta V u_{i-1}^{*} = \frac{u_i^n}{\Delta t} \Delta V - \left. \left( \frac{\partial p}{\partial x} \right) \right|_i^n \Delta V + \left. \left( \rho_k g \right) \right|_i^n \Delta V,$

(20)

where the predicted values are denoted by $*$. Notice that pressure is taken from the old time level $n$ since it is yet unknown, which makes the predicted velocity non-divergence free. Equation 20 is what is actually solved in UEqn.H as described in Section 3.1. However, the cell volumes $\Delta V$ must be divided out as follows to get the correct coefficient matrices and vectors that are used in the corrector step

$\left[ \frac{1}{\Delta t}+\left( \frac{u_{i+\frac{1}{2}}^n-u_{i-\frac{1}{2}}^n}{2 \Delta x} \right) \right] u_i^{*} + \left( \frac{u_{i+\frac{1}{2}}^n}{2 \Delta x} \right) u_{i+1}^{*} - \left( \frac{u_{i-\frac{1}{2}}^n}{2 \Delta x} \right) u_{i-1}^{*} = \frac{u_i^n}{\Delta t} - \left. \left( \frac{\partial p}{\partial x} \right) \right|_i^n + \left. \left( \rho_k g \right) \right|_i^n.$

(21)

In solution vector form, this becomes (where vector refers to the vector of the solution to $u_i^{*}$ at all points $i$) as

$\mathbf{C} \mathbf{u}^{*} = \mathbf{r}-\mathbf{\nabla p}^n + \mathbf{\rho_k g}^n,$

(22)

where $\mathbf{C}$ is the coefficient array multiplying the solution $\mathbf{u}^{*}$ vector and $\mathbf{r}$ are the right-hand side explicit source terms excluding the pressure gradient. One can see that the inclusion of the viscous and turbulent stress terms would modify the coefficient matrix $\mathbf{C}$ and not change the general form of the matrix-vector equation. This equation can be changed to

$\mathbf{A} \mathbf{u}^{*} + \mathbf{H'} \mathbf{u}^{*} = \mathbf{r}-\mathbf{\nabla p}^n + \mathbf{\rho_k g}^n,$

(23)

where $\mathbf{A}$ is the diagonal matrix of $\mathbf{C}$ and $\mathbf{H'}$ is the off-diagonal matrix as $\mathbf{A}$ (in other words, $\mathbf{A}+\mathbf{H'} = \mathbf{C}$). Using a matrix solver, Equation 23 can be readily solved for the predicted velocity $\mathbf{u}^*$. Equation 23 is a generalization of what is solved during the velocity predictor step of buoyantBoussinessqPisoFoam as discussed in Section 3.1.

2 Corrector

Next, the discretized explicit velocity corrector is written as

$\left[ \frac{1}{\Delta t}+\left( \frac{u_{i+\frac{1}{2}}^n-u_{i-\frac{1}{2}}^n}{2 \Delta x} \right) \right] u_i^{**} + \left( \frac{u_{i+\frac{1}{2}}^n}{2 \Delta x} \right) u_{i+1}^{*} - \left( \frac{u_{i-\frac{1}{2}}^n}{2 \Delta x} \right) u_{i-1}^{*} = \frac{u_i^n}{\Delta t} - \left. \left( \frac{\partial p}{\partial x} \right) \right|_i^* + \left. \left( \rho_k g \right) \right|_i^n.$

(24)

Here, the first corrected velocity $\mathbf{u}^{**}$ is being solved from the predicted velocity $\mathbf{u}^{*}$, old velocity $\mathbf{u}^{n}$, and the first corrected pressure $\mathbf{p}^{*}$. The problem is that the corrected pressure is yet unknown --- all that is known is the old pressure. Like in Equation 23, Equation 24 can be expressed in matrix-vector form as

$\mathbf{A} \mathbf{u}^{**} + \mathbf{H'} \mathbf{u}^{*} = \mathbf{r}-\mathbf{\nabla p}^* + \mathbf{\rho_k g}^n.$

(25)

Introducing $\mathbf{H} = \mathbf{r} - \mathbf{H'} \mathbf{u}^*$ and inverting $\mathbf{A}$ (which is easy since it is diagonal), Equation 25 becomes

$\mathbf{u}^{**} = \mathbf{A}^{-1} \mathbf{H} - \mathbf{A}^{-1} \mathbf{\nabla p}^* + \mathbf{A}^{-1} \mathbf{\rho_k g}^n.$

(26)

The point of the corrector step is to make the corrected velocity field divergence free so that it adheres to the continuity equation. Therefore, taking the divergence of Equation 26 and recognizing that $\mathbf{\nabla u}^{**} = 0$ due to continuity (Equation 1) yields a Poisson equation for the first corrected pressure

$\mathbf{\nabla}^2 \left( \mathbf{A}^{-1} \mathbf{p}^{*} \right) = \mathbf{\nabla} \cdot \left( \mathbf{A}^{-1} \mathbf{H} + \mathbf{A}^{-1} \mathbf{\rho_k g}^n \right).$

(27)

With the first corrected pressure $\mathbf{p}^{*}$, Equation 26 can be solved for the first corrected velocity $\mathbf{u}^{**}$.

Further correction steps can be applied using the same $\mathbf{A}$ matrix and $\mathbf{H}$ vector (this is convenient computationally since they can be stored in the computer's memory once and recalled as necessary). For example the second correction step would consist of the equations.

$\mathbf{\nabla}^2 \left( \mathbf{A}^{-1} \mathbf{p}^{**} \right) = \mathbf{\nabla} \cdot \left( \mathbf{A}^{-1} \mathbf{H} + \mathbf{A}^{-1} \mathbf{\rho_k g}^n \right)$

(28)

and

$\mathbf{u}^{***} = \mathbf{A}^{-1} \mathbf{H} - \mathbf{A}^{-1} \mathbf{\nabla p}^{**} + \mathbf{A}^{-1} \mathbf{\rho_k g}^n,$

(29)

where $\mathbf{p}^{**}$ and $\mathbf{u}^{***}$ are the second corrected pressure and velocity, respectively. Equations 26 and 27 are a generalization of what is solved during the corrector step in buoyantBoussinessqPisoFoam as discussed in Section 3.3.

In this example, first-order Euler implicit time stepping with second-order linear interpolation has been used. This method is equally applicable with other forms of implicit time stepping, such as Crank-Nicholson or second-order backward, and interpolation schemes. The same general results will follow but with modifications to the $\mathbf{A}$ matrix and $\mathbf{H}$ vector. Issa states that if a second-order accurate time stepping scheme is used, then three corrector steps should be used to reduce the discretization error due to the PISO algorithm to second-order. For flow in which a temperature equation (or other equations) are necessary, Oliveira and Issa also includes a corrector for those equations in the corrector loop. buoyantBoussinesqPisoFoam does not do this --- temperature is simply predicted from the last time level and never corrected. Therefore, it may need to be modified to provide temperature correction.