The PPE formulation

$ \renewcommand{\vec}[1]{ {\bf #1} } \newcommand{\bnabla}{ \vec{\nabla} } \newcommand{\Rey}{Re} \def\vechat#1{ \hat{ \vec{#1} } } \def\mat#1{#1} $ The incompressibility condition may be 'replaced' by specifying an equation for the pressure. Taking the divergence of the Navier--Stokes equation leads to the system

$\label{eq:NSeqs}

        \left\{\begin{array}{rcl}
    \partial_t \vec{u} & = & L \,\vec{u} + \vec{N} - \bnabla p \, , \\
                \nabla^2 p & = & \bnabla\cdot\vec{N} ,
        \end{array}\right.$

where $L$ is a linear operator and $\vec{N}$ represents nonlinear terms. The boundary conditions are $\vec{u}=\vec{0}$, and we retain that $\bnabla\cdot\vec{u}=0$ must be satisfied everywhere, i.e. also on the boundary. There is no boundary condition explicitly on the pressure.

PPE-formulation with correct boundary conditions

Consider the case where the spatial discretisation splits the Navier--Stokes equations into a set one-dimensional problems, here into problems for radially-dependent Fourier modes.

Let $\vec{u}$ denote the velocity, $\vec{N}$ denote nonlinear terms, and $\mat{X}$ and $\mat{Y}$ be matrices associated with implicit timestepping of the viscous terms for a particular mode. The time-discretised Navier–Stokes equations for this mode may be written in the form

$\label{eq:NSdisc}

        \left\{\begin{array}{rcl}
    \mat{X}\, \vec{u}^{q+1} & = & \mat{Y}\, \vec{u}^q 
                + \vec{N}^{q+\frac{1}{2}} - \bnabla p \, , \\
                \nabla^2 p & = & \bnabla\cdot(\mat{Y}\, \vec{u}^q 
                + \vec{N}^{q+\frac{1}{2}}) ,
        \end{array}\right.$

where $q$ denotes time $t_q$, which is sixth order in $r$ for $\vec{u}^{q+1}$ and second order for $p$, where the solenoidal condition is not explicitly imposed. Symmetry conditions provide the conditions at the axis. The difficulty is in imposing the remaining four — this system should be inverted, in principle, simultaneously for $p$ and $\vec{u}^{q+1}$ with boundary conditions $\vec{u}^{q+1}=\vec{0}$ and $\bnabla\cdot\vec{u}^{q+1}=0$ on $r=R$ (Rempfer 2006). In practice it would be preferable to invert for $p$ first then for $\vec{u}^{q+1}$, but the boundary conditions to not involve $p$ directly.

Note that the $\mat{Y}\,\vec{u}^q$ term has been included in the right-hand side of the pressure-Poisson equation, the divergence of which should be small. Assume that pressure boundary condition is known: the right-hand side of the Navier–Stokes equation is then projected onto the space of solenoidal functions though $p$ and hence after inversion, $\vec{u}^{q+1}$ will be solenoidal.

Consider the ‘bulk’ solution, $\{\bar{\vec{u}},\bar{p}\}$, obtained from solution of the following:

$\label{eq:NSbulk}

        \left\{\begin{array}{rcl}
    \mat{X}\, \bar{\vec{u}} & = & \mat{Y}\, \vec{u}^q 
                + \vec{N}^{q+\frac{1}{2}} - \bnabla \bar{p} \, , \\
                \nabla^2 \bar{p} & = & \bnabla\cdot(\mat{Y}\, \vec{u}^q 
                + \vec{N}^{q+\frac{1}{2}}) ,
        \end{array}\right.$

with boundary conditions $\bar{\vec{u}}=\vec{0}$ and $\partial_{r}\bar{p}=0$. Introduce the following systems:

$\label{eq:NSp0}

        \left\{\begin{array}{rcl}
         \vec{u}' & = & -\bnabla p' \, , \\
                \nabla^2 p' & = & 0,
        \end{array}\right.$

with boundary condition $\partial_{r}p'=1$ on $r=R$, and

$\label{eq:NSu0}

        \left\{\begin{array}{rcl}
    \mat{X}\, \vec{u}' & = & \vec{0},
        \end{array}\right.$

with boundary conditions $u'_+=1$, $u'_-=1$, $u'_z=\mathrm{i}$ on $r=R$ (see Decoupling_the_equations). The system \ref{eq:NSp0} provides a linearly independent function $\vec{u}'_1$ that may be added to $\bar{\vec{u}}$ without affecting the right-hand side in \ref{eq:NSbulk}, but altering (to correct) the boundary condition applied. ($\mat{X}$ is dropped on the left-hand side of \ref{eq:NSp0} since $\bnabla^2(\bnabla p')=\bnabla(\nabla^2 p')-\bnabla \wedge\bnabla \wedge\bnabla p'=0$.) Similarly the system \ref{eq:NSu0} provides a further three functions $\vec{u}'_j$ (where $j=2,3,4$; the +,- and $z$ components may be considered separately). The superposition

$\label{eq:usuperpos}

  \vec{u}^{q+1} = \bar{\vec{u}} + \sum_{j=1}^4 a_j\, \vec{u}'_j \,$

may be formed in order to satisfy the four boundary conditions, $\vec{u}^{q+1}=\vec{0}$ and $\bnabla\cdot\vec{u}^{q+1}=0$ on $r=R$. Substituting \ref{eq:usuperpos} into the boundary conditions, they may be written

$\label{eq:velBCs}

  \mat{A}\,\vec{a} = -\vec{g}(\bar{\vec{u}}) ,$

where $\mat{A}=\mat{A}(\vec{g}(\vec{u}'))$ is a 4$\times$4 matrix. The appropriate coefficients required to satisfy the boundary conditions are thereby recovered from solution of this small system for $\vec{a}$. If needed, the pressure may be calculated using the same $a_j$:

$\label{eq:pressure}

  p^{q+1} = \bar{p} + \sum_{j=1}^4 a_j\, p'_j \, .$

The error in the boundary conditions $g_j(\vec{u}^{q+1})$ using the influence-matrix technique is at the level of the machine epsilon, typically order 1e-18. The functions $u'_j(r)$, the matrix $\mat{A}$ and its inverse may all be precomputed. The boundary conditions for $\vec{u}'$ have been chosen so that that $u'_\pm$ are pure real, $u'_z$ is pure imaginary, and $\mat{A}$ is real. For each timestep, this application of the influence matrix technique requires only evaluation of the deviation from the boundary condition, multiplication by a 4$\times$4 real matrix, and the addition of only two functions to each component of $\vec{u}$, each either pure real or pure imaginary. Compared to the evaluation of nonlinear terms, the computational overhead is negligible.