Application - Pressure Poisson Equation

The momentum equation for the velocity field $\vec{u}$ in a fluid is

$$\begin{equation} \frac{\partial \vec{u}}{\partial t}+(\vec{u}\cdot\nabla)\vec{u}=-\frac{1}{\rho}\nabla p + \nu \nabla^2\vec{u} \end{equation}$$

where $p$ is the pressure, $\nu$ is the fluid viscosity and $\rho$ is the fluid density. With three velocity components, plus the pressure, we have four unknowns but only three equations. For compressible fluids, we have an equation of state to complete the system. In the incompressible case, we don't have an equation of state and we need an additional constraint from somewhere else.

This is what we do: take the divergence of the momentum equation, apply the incompressibility constraint to cancel some terms, and get an equation for the pressure. It's a pretty cool trick.

Conservation of mass for an incompressible fluid requires that the divergence of $\vec{u}$ must be zero:

$$\nabla \cdot \vec{u} = 0$$

Writing out the momentum equation in $x$ and $y$ components (for two-dimensional flow), we get

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial x}+\nu \left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right)$$

$$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial y}+\nu\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right)$$

We take the divergence of the momentum equation and then apply the incompressibility constraint. After some wrangling and cancellations, this leaves us with the pressure Poisson equation:

$$\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} = -\rho\left(\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial v}{\partial y} \right)$$

Which is an equation of the form

$$\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = b$$

Imagine we discretize a domain using a uniform mesh of points in each spatial direction, as in the figure below:

Then the left-hand side of the Poisson equation, i.e., the Laplacian differential operator applied to $p$, is discretized using 2nd-order central differences as follows

$$\frac{p^n_{i+1, j} - 2p^n_{i,j} + p^n_{i - 1, j}}{\Delta x ^2} + \frac{p^n_{i, j+1} - 2p^n_{i,j} + p^n_{i, j-1}}{\Delta y ^2}$$

where subscripts $i,j$ denote the spatial location on a Cartesian coordinate system and superscripts $n$ denote a point in time.

We apply an appropriate finite-difference discretization to the momentum equation (forward-time, backward-space for the 1st-order terms) and also assume a uniform mesh, so $\Delta x = \Delta y$.

Using this discretized form in the Poisson equation, we will leave only the $p_{i,j}$ terms in the left-hand side, and move the other terms to the right. Then we say that we can update all the values of $p_{i,j}$ using the values at the neighboring points for both $p$ and $u, v$. This update, repeated many times, happens to converge to the solution of Poisson's equation.

Solution procedure

Initial velocity field

We start with a velocity field in $u$ and $v$ at some timestep $n$.

Calculate pressure

Then, we iteratively solve the Poisson equation for pressure, as described above. Starting with an initial guess, the values $p_{i,j}$ are updated using the neighboring values of $p$, $u$ and $v$ at $(i+1,j)$ and $(i,j+1)$. The updates can be written as follows, where the $k$ superscript denotes an iteration in 'pseudo-time':

$$\begin{align} p_{i,j}^{k+1} &= \frac{1}{4}\left(p_{i+1,j}^{k}+p_{i-1,j}^{k}+p_{i,j+1}^{k}+p_{i,j-1}^{k}\right) \\ &-\frac{\rho \Delta x}{16} \left( \frac{2}{\Delta t} \left(u_{i+1,j} - u_{i-1,j} + v_{i,j+1} - v_{i,j-1}\right) \right . \\ &-\frac{2}{\Delta x}\left(u_{i,j+1} - u_{i,j-1} \right) \left(v_{i+1,j} - v_{i-1,j} \right) \\ &- \left . \frac{\left(u_{i+1,j} - u_{i-1,j} \right)^2}{\Delta x} - \frac{ \left(v_{i,j+1} - v_{i,j-1} \right)^2 }{\Delta x} \right) \\ \end{align}$$

In other words, we repeatedly apply the Poisson equation until the pressure reaches a quasi-steady state.

Update the velocity

Once the pressure field reaches its quasi-steady state via the Poisson equation, we use that field for the current time step, $p^n$, to solve for the velocity components $u$ and $v$ at the next timestep, $n+1$.

The momentum equation in the $u$ direction:

$$\begin{align} u_{i,j}^{n+1} = u_{i,j}^{n} &- \frac{\Delta t}{\Delta x} \left( u_{i,j}^{n}(u_{i,j}^{n}-u_{i-1,j}^{n}) + v_{i,j}^{n} (u_{i,j}^{n}-u_{i,j-1}^{n}) + \frac{1}{2 \rho}(p_{i+1,j}^{n}-p_{i-1,j}^{n}) \right) \\ &+\frac{\nu \Delta t}{\Delta x^2}\left(u_{i+1,j}^{n} + u_{i-1,j}^{n} + u_{i,j+1}^{n} + u_{i,j-1}^{n} -4u_{i,j}^{n}\right) \end{align}$$

The momentum equation in the $v$ direction:

$$\begin{align} v_{i,j}^{n+1} = v_{i,j}^{n} &- \frac{\Delta t}{\Delta x} \left( u_{i,j}^{n}(v_{i,j}^{n}-v_{i-1,j}^{n}) + v_{i,j}^{n} (v_{i,j}^{n}-v_{i,j-1}^{n}) + \frac{1}{2 \rho}(p_{i,j+1}^{n}-p_{i,j-1}^{n}) \right) \\ &+\frac{\nu \Delta t}{\Delta x^2}\left(v_{i+1,j}^{n} + v_{i-1,j}^{n} + v_{i,j+1}^{n} + v_{i,j-1}^{n} -4v_{i,j}^{n}\right) \end{align}$$

Then, rinse and repeat.

What we left out

There are various subtleties that we left out here, to get quickly to the equations we need in the code. First, there are some variations on the form of the pressure Poisson equation, depending on what terms involving $\nabla\cdot\vec{u}$ one chooses to cancel (this has caused long arguments in the literature!). Second, we say nothing about the boundary conditions, which can cause some trouble (and more arguments!). And third, we show only the simplest iterative method for solving the Poisson equation, which also happens to be the slowest to converge. This is just meant to be a pedagogical example and discussion of these subtleties would be part of a full-fledged CFD course.