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 - the Pressure-Poisson Equation (PPE). Taking the divergence of the Navier--Stokes equation leads to the system

$\label{eq:NSeqs} (1) \qquad \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 no-slip 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. (See Boundary_conditions.)

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 Fourier mode. The time-discretised Navier–Stokes equations for this mode may be written in the form

$\label{eq:NSdisc} (2) \qquad \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 provides the conditions at the axis. The difficulty is in imposing the remaining four boundary conditions — 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} (3) \qquad \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} (4) \qquad \left\{\begin{array}{rcl} \mat{X}\,\vec{u}' & = & -\bnabla p' \, , \\ \nabla^2 p' & = & 0, \end{array}\right.$

with boundary conditions $\vec{u}'=\vec{0}$ and $\partial_{r}p'=1$ on $r=R$, and

$\label{eq:NSu0} (5) \qquad \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 $(4)$ provides a linearly independent function $\vec{u}'_1$ that may be added to $\bar{\vec{u}}$ without affecting the right-hand side in $(3)$, but altering (to correct) the boundary condition applied. Similarly the system $(5)$ 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} (6) \qquad \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$. Let $\vec{g}(\vec{u})$ be a 4-vector composed of these boundary conditions, such that they are satisfied when $\vec{g}(\vec{u})=\vec{0}$. Substituting the superposition $(6)$ into the boundary conditions, they may be written

$\label{eq:velBCs} (7) \qquad \mat{A}\,\vec{a} = -\vec{g}(\bar{\vec{u}}) ,$

where $\mat{A}=\mat{A}(\vec{g}(\vec{u}'))$ is a 4$\times$4 matrix and the 4-vector $\vec{a}$ is composed of the $a_j$. Thus, the appropriate coefficients required to satisfy the boundary conditions are recovered from solution of this small system for $\vec{a}$.

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-16. 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.

If needed, the pressure may be calculated using the same $a_1$ as in $(6)$ using the $p'$ from $(4)$:

$\label{eq:pressure} (8a) \qquad \tilde{p} = \bar{p} + a_1\, p' \, \qquad$ (for each Fourier mode, then)

$ (8b) \qquad p = \tilde{p} - \frac{1}{2}\vec{u}\cdot\vec{u}\, \qquad$ (in physical space).

The adjustment in $(8b)$ arises when the Navier-Stokes equations are in rotational form (see Equations_and_parameters).