In this chapter we will consider the pressure term. That is the under-lined term in the Navier-Stokes equations below,

$\frac{\partial \mathbf{u}}{\partial t} = -(\mathbf{u} \cdot \nabla)\mathbf{u} - \underline{\nabla p} + \nu \nabla^2 \mathbf{u}.$

This term is perhaps the most involved and least intuitive of these equations. This is perhaps because the pressure in a fluid is a somewhat illusive concept. As already mentioned, we will find an implicit problem that needs to be solved in an intermediate step before actually computing the pressure gradient. To add insult to injury, we cannot really conceptualize this term for a helper scalar field, so we will do the full vector calculus.

Before we start, we should get some common misconceptions out of the way. If you have learned about forced flows (e.g. pipe flow) in your fluid dynamics classes, your intuition may tell you that the pressure gradient in externally applied and prescribed by some sort of pumping mechanism. Although a pressure drop can be treated as a boundary condition, a more relevant view (from the perspective of our solver) is demonstrated by the Bernoulli principle. If a steady pipe flow contains some constriction (e.g. a Venturi), the approaching flow accelerates towards the narrowing. This accelerating force comes from the pressure gradient as Bernoulli realized this region of high velocities is associated with a lower pressure. It thus appears that the pressure field is directly related to the flow itself. But how?

Another confusing conception is the thermodynamic meaning of pressure. Especially for those who are versed incompressible flows it can be hard to go without some equation of state. Indeed, for in compressible flows in simple (i.e. at least rigid) geometries, pressure loses its thermodynamic meaning. The pressure can be seen to adjust such that the flow remains incompressible. This entails that we do not need to time integrate some evolution equation for the pressure as we can derive it directly as a function of the flow.

## A projection method

The fundamental theorem of vector calculus states that any vector field $\mathbf{v}$ can be decomposed as,

$\mathbf{v} = \mathbf{R} + \nabla \phi,$

where,

$\nabla \cdot \mathbf{R} = 0.$

such that the divergence in $\mathbf{v}$ is captured by the scalar-gradient term. We already have the tools to make this so-called Helmholtz decomposition numerically. If we take the divergence of the equation,

$\nabla \cdot \mathbf{v} = \nabla \cdot \mathbf{R} + \nabla \cdot \nabla \phi,$ $\nabla \cdot \mathbf{v} = 0 + \nabla \cdot \nabla \phi.$

We see the structure of a Poisson equation where the source term is the divergence of the vector field $\mathbf{v}$. If we have solved for $\phi$, we can also find $\mathbf{R}$ by,

$\mathbf{R} = \mathbf{v} - \nabla \phi,$

completing the decomposition. At the risk of sounding cryptic, we could say that $\mathbf{R}$ is the projection of $\mathbf{v}$ onto the space of divergence-free vector fields1. Lets code it up! In ns.h,

...

void diffusion (...) {
...
}

void projection (double vx[N][N], double vy[N][N], double phi[N][N]) {

}
...

Note that it will prove important to have the phi field as input, altough we could in principle declare and initialize it in the function’s scope. The function takes the steps described above, overwriting the input vector field.

#include "common.h"
#include "poisson.h"
...
void projection (double vx[N][N], double vy[N][N], double phi[N][N]) {
// Compute the divergence
double div[N][N];
foreach() {
val(div, 0, 0) = (val(vx, 1, 0) - val(vx, -1, 0))/(2*Delta);
val(div, 0, 0) += (val(vy, 0, 1) - val(vy, 0, -1))/(2*Delta);
}
// Compute the helper scalar field
poisson (phi, div);
// Reject the divergent part
foreach() {
val (vx, 0, 0) -= (val(phi, 1, 0) - val(phi, -1, 0))/(2*Delta);
val (vy, 0, 0) -= (val(phi, 0, 1) - val(phi, 0, -1))/(2*Delta);
}
}
...

For now, I propose to only test the syntax of this function (using gcc), and not do a proper test, as we will already change it in the next section and we will test it soon enough in a proper flow setting.

## Chorin’s projection method

In Chorin’s seminal paper2 on numerical flow-problem solving, the above projection method is used to compute the pressure term. His special insight was on how to compute it in conjunction with the advection and diffusion terms. As such, Chorin’s projection method may also referred to as Chorin’s operator splitting, describing how and in what sequence the various terms should be handled for time advancement. In the first stage of our time step, a provisional velocity field ($\mathrm{u}^*$) is computed at $t_{n + 1} = t_n + \mathrm{d}t$, which considers every term except the pressure term.

$\mathrm{u}^* = \mathbf{u}_n + \mathrm{d}t \left( -(\mathbf{u}_n \cdot \nabla)\mathbf{u}_n + \nu \nabla^2 \mathbf{u}_n\right).$

We can find $\mathbf{u}_{n + 1}$ by projecting $\mathbf{u}^*$ on to the space of divergence-free vector fields. But what about the pressure? phi can not be it as is does not have the right units. Well, it is the pressure, but its off by a factor $\mathrm{d}t$ . See,

$\mathbf{u}_{n + 1} = \mathbf{u}^* - \mathrm{d}t \nabla p,$

taking the divergence,

$\nabla \cdot \mathbf{u}_{n + 1} = \nabla \cdot \mathbf{u}^* - \nabla \cdot \mathrm{d}t \nabla p,$ $0 = \nabla \cdot \mathbf{u}^* - \nabla^2 \mathrm{d}tp,$

This is the previous Poisson equation with $\phi \leftarrow \mathrm{d}tp$.

$\nabla \cdot \mathbf{u}^* = \nabla^2 \mathrm{d}tp,$

Or even better, we can divide by the time-step size ($\mathrm{d}t$)

$\frac{\nabla \cdot \mathbf{u}^*}{\mathrm{d}t} = \nabla^2 p.$

In view of this, I propose to alter the projection() function a little,

void projection (double vx[N][N], double vy[N][N],
double phi[N][N], double dt) {
// Compute the divergence devided by dt
double div[N][N];
foreach() {
val(div, 0, 0) = (val(vx, 1, 0) - val(vx, -1, 0))/(2*Delta);
val(div, 0, 0) += (val(vy, 0, 1) - val(vy, 0, -1))/(2*Delta);
val(div, 0, 0) /= dt;
}
// Compute the pressure field
poisson (phi, div);
// Reject the divergent part
foreach() {
val (vx, 0, 0) -= dt*(val(phi, 1, 0) - val(phi, -1, 0))/(2*Delta);
val (vy, 0, 0) -= dt*(val(phi, 0, 1) - val(phi, 0, -1))/(2*Delta);
}
}

Using $p$ as the helper scalar field for projection instead of any other is important as it means that for slowly evolving flows, the pressure field from the previous time step will serve as a good initial guess. This would not be the case for $\phi$, which could change a lot with varying time-step sizes. Further, scaling the divergence with dt introduces a good moment to consider the meaning of our tolerance on the maximum residual (tolerance). If we wish the tolerance to represent the maximum divergence in our solution we should scale it with the time step size.

void projection (double vx[N][N], double vy[N][N],
double phi[N][N], double dt) {
// store desired tolerance and set tolerance on divergence/dt
double tol = tolerance;
tolerance /= dt;
...
// Reset tolerance
tolerance = tol;
}

We cant really diagnose this low-tolerance divergence as the gradients in our Poisson solver were defined in a staggered fashion. Any divergence we can approximate will be polluted with an additonal discretization error. Further, our projection method is only approximately as we did not opt for the 5-point stencil version.

Notice that with Chorin’s pressure-term treatment, we are not doing a regular forward Euler step, but rather a backward step for $p$, as it is evaluated for the solution at $t_{n+1}$.

At this moment, we can combine our functions in to a Navier-Stokes equation time-stepping function. We will do this in the next chapter.

1. It seems more like a vector rejection to me.↩︎

2. Chorin, Alexandre Joel. “Numerical solution of the Navier-Stokes equations.” Mathematics of computation 22.104 (1968): 745-762.↩︎

The marvelous design of this website is taken from Suckless.org