# 9. The pressure-gradient force

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 fields^{1}.
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(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}
// Compute the helper scalar field
(phi, div);
poisson // Reject the divergent part
() {
foreach(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);
val }
}
...
```

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 paper^{2} 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(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;
val}
// Compute the pressure field
(phi, div);
poisson // Reject the divergent part
() {
foreach(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);
val }
}
```

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;
/= dt;
tolerance ...
// Reset tolerance
= tol;
tolerance }
```

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.