Putting the Navier–Stokes Equations Together: The Pressure Poisson Equation

Navier–Stokes, part 4: joining acceleration, forces and the no-pile-up rule, and the equation that finds the pressure at every time step

October 15th, 2026

navier-stokes
incompressible flow
pressure poisson equation
jacobi iteration
fluid dynamics
pde
physics simulation

Posts 10, 11 and 12 each built one piece of the Navier–Stokes equations: the acceleration of a fluid parcel, the pushes that act on it, and the rule that water cannot pile up. This post joins the pieces, then asks how to step the result forward in time, as Posts 2–6 did for the heat and wave equations. The answer leads to the pressure, which has no update rule of its own, and to the equation that determines it at every time step: the pressure Poisson equation.

Three loose pieces

Posts 10–12 left three separate results: (A) the parcel's acceleration, (B) the pushes on it per kilogram, and (C) the constraint that water cannot pile up. On their own, they don't talk to each other, and none of them can predict anything. (A) is bookkeeping: it computes an acceleration from a given velocity field, but doesn't say what the acceleration will be. (B) lists pushes, but says nothing about what they do. (C) only forbids certain fields.

$\begin{aligned} \text{(A)}\quad & \frac{D\mathbf{u}}{Dt} = \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} \\ \text{(B)}\quad & -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} + \mathbf{g} \\ \text{(C)}\quad & \nabla\cdot\mathbf{u} = 0 \end{aligned}$

Newton's law $F = ma$ is the glue between (A) and (B), where $F$ denotes a force (N), $m$ a mass (kg) and $a$ an acceleration (m/s²). It is the same move as on the string in Post 3, whose acceleration and tension force only became the wave equation once $F = ma$ joined them. (C) is the one genuinely new ingredient: the string had no rule that forbade certain states. The next three sections recap each piece on its own.

Piece (A): the parcel's acceleration

$F = ma$ is about one piece of matter, but our data is Eulerian: a table of velocities at fixed places, through which water flows (Post 10). The material derivative computes from it the acceleration of the parcel passing through. In words: a parcel's velocity changes for two independent reasons, and its total acceleration is the sum of both.

$\underbrace{\dfrac{D\mathbf{u}}{Dt}}_{\text{the answer}} = \underbrace{\dfrac{\partial \mathbf{u}}{\partial t}}_{\substack{\text{reason 1:} \\ \text{the river changes}}} + \underbrace{(\mathbf{u}\cdot\nabla)\mathbf{u}}_{\substack{\text{reason 2: drifting} \\ \text{into different water}}}$(1)

Where $\mathbf{u} = (u_x, u_y, u_z)$ denotes the velocity field (m/s), $D\mathbf{u}/Dt$ the acceleration of the parcel currently at position $\mathbf{x} = (x, y, z)$ (m/s²), $\partial \mathbf{u}/\partial t$ the change of the velocity at that fixed place, which is what a sensor on a bridge pillar reports, and $(\mathbf{u}\cdot\nabla)\mathbf{u}$ the change because the parcel drifts into water that moves differently, with $(\mathbf{u}\cdot\nabla) := u_x\,\partial/\partial x + u_y\,\partial/\partial y + u_z\,\partial/\partial z$ (Post 10). Equation (1) contains no force and no mass, so no physics.

Piece (B): the pushes per kilogram

A parcel has no natural size: a bigger parcel feels more force, but it is heavier by the same factor. Dividing by the mass $m = \rho\,\Delta V$ cancels the volume $\Delta V$ and leaves N/kg = m/s², an acceleration that belongs to the place, not to our arbitrary parcel (Post 12):

$\underbrace{-\frac{1}{\rho}\nabla p}_{\text{pressure}} + \underbrace{\nu\nabla^2\mathbf{u}}_{\text{friction}} + \underbrace{\mathbf{g}}_{\text{gravity}}$(2)

Where $p$ denotes the pressure (Pa), $\rho$ the constant density (kg/m³), $\nu$ the kinematic viscosity (m²/s), $\nabla^2\mathbf{u}$ the Laplacian of the velocity (1/(m·s)), and $\mathbf{g} = (0, 0, -9.81)$ m/s² gravity. Pressure: only differences push; $\nabla p$ points uphill, the minus sign pushes from high to low pressure, and $1/\rho$ makes light fluids accelerate more.

Friction pulls each parcel's velocity toward the average of its neighbours, the heat equation's smoothing applied to velocity (Post 12):

$\nu\nabla^2\mathbf{u} = \nu\begin{bmatrix} \dfrac{\partial^2 u_x}{\partial x^2} + \dfrac{\partial^2 u_x}{\partial y^2} + \dfrac{\partial^2 u_x}{\partial z^2} \\[2mm] \dfrac{\partial^2 u_y}{\partial x^2} + \dfrac{\partial^2 u_y}{\partial y^2} + \dfrac{\partial^2 u_y}{\partial z^2} \\[2mm] \dfrac{\partial^2 u_z}{\partial x^2} + \dfrac{\partial^2 u_z}{\partial y^2} + \dfrac{\partial^2 u_z}{\partial z^2} \end{bmatrix}$(3)

Where $\partial^2 u_x/\partial y^2$ denotes the curvature of $u_x$ along $y$ (1/(m·s)), and likewise for every other pairing. Each row is one velocity component, and only that component appears in it: the x-velocity is smoothed only with the neighbours' x-velocities. Within a row, each term is one pair of neighbours (left/right, below/above, behind/in front), and the pairs add. In 2D, the third row and the z-terms drop out, and on a grid, each second slope is Post 2's three-point stencil:

$\begin{aligned} \nabla^2 u_x &\approx \frac{u_{x,\text{right}} + u_{x,\text{left}} + u_{x,\text{up}} + u_{x,\text{down}} - 4\,u_{x,\text{here}}}{h^2} \\ &= \frac{4}{h^2}\Bigg(\underbrace{\frac{u_{x,\text{right}} + u_{x,\text{left}} + u_{x,\text{up}} + u_{x,\text{down}}}{4}}_{\text{neighbours' average}} - u_{x,\text{here}}\Bigg) \end{aligned}$(4)

The bracket is "neighbours minus me", so the sign is already right: positive means I am slower than my neighbours, and friction accelerates me toward $+x$; negative means I am faster, and friction holds me back.

So $\nabla^2\mathbf{u} = (2, -1)$ does not mean "the difference along x and along y": +2 says my x-velocity is below my neighbours' average x-velocity (pushed toward $+x$), and −1 says my y-velocity is above theirs (pushed toward $-y$). Its unit is 1/(m·s), not m/s, so "2" is not "2 m/s faster"; only $\nu$ turns it into m/s².

Gravity gives every kilogram the same $\mathbf{g}$, and at rest the hydrostatic pressure balances it (Post 12).

Piece (C): the no-pile-up rule

Nothing in (A) or (B) knows that water cannot be squeezed. (C) says: every tiny box sends out exactly as much water as it receives.

$\frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} + \frac{\partial u_z}{\partial z} = 0$(5)

Where $\partial u_x/\partial x$ denotes how fast $u_x$ changes along $x$ (1/s), likewise for $y$ and $z$; the sum is $\nabla\cdot\mathbf{u}$, the net fraction of a tiny box's volume flowing out per second (Post 11). Unlike in the friction vector, each component is differentiated only along its own direction, because only flow through a face counts, and the result is one number. It is Post 11's continuity equation with constant density (squeeze one way, bulge the other), and with no time derivative, it is a constraint, not an evolution equation.

Why (A) and (B) cannot stop a pile-up. Together they are Newton's law applied to each parcel on its own, with pushes only from its surroundings; no term checks whether, after all parcels have moved, some box receives more than it sends. The brake is in (B): crowding raises the pressure, which pushes the fluid back out. But nothing says what value $p$ takes. In a gas, crowding raises the density and with it the pressure, so the rule comes for free. In water the density never changes, the chain is cut, and something else must fix $p$: (C). (B) provides the brake; (C) sets how hard to press it.

Assembling the equations

$F = ma$ divided by $m$ gives $a = F/m$. (B) already is force per kilogram, so (A) becomes the left side and (B) the right side:

$\underbrace{\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u}}_{\text{(A): acceleration}} = \underbrace{-\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} + \mathbf{g}}_{\text{(B): pushes per kg}}$(6)

Where the left side is Equation (1) and the right side Equation (2), all in m/s²: at every place, the acceleration of the parcel passing through equals the sum of the pushes per kilogram on it. It is a vector equation, so three equations.

With Post 12's split $p = p_{\text{rest}} + p_{\text{stir}}$ into the hydrostatic part and the extra from the motion, $-\frac{1}{\rho}\nabla p_{\text{rest}} + \mathbf{g} = \mathbf{0}$, so papers drop $\mathbf{g}$ and rename $p_{\text{stir}}$ to $p$ (valid for constant density and no free surface). A driven flow gets a body force per kilogram $\mathbf{f}$ (m/s²) instead, lowercase because a capital $F$ would be in newtons. From here on it is left out.

(C) cannot be added into this vector equation in m/s², because it is one number in 1/s. The two are written side by side, joined by "and", and both must hold at the same time:

$\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} \qquad \text{and} \qquad \nabla\cdot\mathbf{u} = 0$(7)

Where $p$ now denotes the pressure caused by the motion, and the other symbols are as in Equation (6). These are the incompressible Navier–Stokes equations; the left one is called the momentum equation. A solution is a pair of fields $(\mathbf{u}, p)$ that satisfies both at every place and every instant. Note who appears where: $\mathbf{u}$ in both equations, $p$ only in the momentum equation, while (C) involves $\mathbf{u}$ alone.

What solving these equations means

The goal is the same as for the heat rod in Post 2. The state at $t = 0$ is given; nobody solves for it. The PDE is a rule that computes the rate of change from the current snapshot, an Euler step does new = old + $\Delta t$ × rate, where $\Delta t$ denotes the time step (s), and this repeats. "Solving the system" means: does it give us a rule to step everything forward? Rearranging the momentum equation so that the rate stands alone gives

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

Where the symbols are as in Equation (7). Everything on the right can be computed from the current snapshot of $\mathbf{u}$, from its slopes and curvatures. Everything except $\nabla p$.

Note that $\partial \mathbf{u}/\partial t$ is the change at a fixed place (the bridge sensor of Post 10), not a parcel's acceleration. The reading there changes for two reasons. The pushes $-\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u}$ accelerate the water at that place. And the water at that place keeps being replaced by water from upstream, which brings its own velocity: this is $-(\mathbf{u}\cdot\nabla)\mathbf{u}$, the delivery. It is not a force, just the convective term of Post 10 seen from the sensor instead of from the leaf. In short: change at a fixed place = delivery + pushes.

$ change at a fixed place = delivery + pushes

// eq. (8)

● water that was here● water delivered from upstream● pressure push● friction push● new readingΔt = 1 s · illustrative numbers

delivery −(u·∇)u

+1.0 m/s²

pressure −(1/ρ)∇p

−0.5 m/s²

friction ν∇²u

−0.3 m/s²

at the sensor ∂u/∂t

+0.2 m/s²

the parcel itself slows down (pushes: −0.8 m/s²), yet the sensor reading rises (+0.2 m/s²): faster water was delivered.

The missing piece

To evaluate $\nabla p$, we need $p$, and we don't have it. Putting $p$ into the initial condition fails twice. It only helps once: after the first step, the pressure at $t = \Delta t$ is needed, and there is no equation $\partial p/\partial t = \dots$ to update it. And $p$ isn't free to choose: with the wrong $p$, the new velocity field violates (C) and fluid piles up; only one pressure field keeps (C) satisfied. So the initial condition is only $\mathbf{u}(\mathbf{x}, 0)$, and pressure is not part of the state: it is recomputed from the current state at every step. The question for the rest of the post is: given the current snapshot of $\mathbf{u}$, how do we compute the one pressure field that keeps $\nabla\cdot\mathbf{u} = 0$ true after the step?

The Pressure Term

Take Equation (8) with the current snapshot of $\mathbf{u}$ in hand. Every term on its right side can be computed from that snapshot: the delivery from the slopes of $\mathbf{u}$, the friction from its curvatures. The only thing we do not have is $p$. If we had it, we could compute $\partial \mathbf{u}/\partial t$ at every point and take the Euler step, exactly as for the heat rod. So the whole problem comes down to one question: where does $p$ come from?

To look for it, write out everything the system offers: the momentum equation, one line per velocity component, and (C):

$\begin{aligned} \frac{\partial u_x}{\partial t} + (\mathbf{u}\cdot\nabla)\,u_x &= -\frac{1}{\rho}\frac{\partial p}{\partial x} + \nu\nabla^2 u_x && \text{(x-component)} \\ \frac{\partial u_y}{\partial t} + (\mathbf{u}\cdot\nabla)\,u_y &= -\frac{1}{\rho}\frac{\partial p}{\partial y} + \nu\nabla^2 u_y && \text{(y-component)} \\ \frac{\partial u_z}{\partial t} + (\mathbf{u}\cdot\nabla)\,u_z &= -\frac{1}{\rho}\frac{\partial p}{\partial z} + \nu\nabla^2 u_z && \text{(z-component)} \\ \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} + \frac{\partial u_z}{\partial z} &= 0 && \text{(C)} \end{aligned}$(9)

The three velocity lines contain $p$, but each also contains a rate of change ($\partial u_x/\partial t$, $\partial u_y/\partial t$ or $\partial u_z/\partial t$), which is exactly what we want to compute and do not know either. One equation with two unknowns gives us neither.
The (C) line contains no $p$ at all. It only talks about the velocity.
No line contains $\partial p/\partial t$, so there is no rule to step the pressure forward the way we step the velocity.

So there is no equation of the form "$p$ = something computed from $\mathbf{u}$". We have to build one, and since (C) is the only line not already used to step the velocity, that is where it has to come from.

An equation for the pressure

If $p$ can be computed from $\mathbf{u}$, a loop still works: given $\mathbf{u}^n$, compute $p^n$ from it, compute the rate with Equation (8), Euler-step to $\mathbf{u}^{n+1}$, and repeat. The path to the core problem:

$p$ cannot be given in the initial condition, because one time step later it would be missing again: there is no $\partial p/\partial t$.
So $p$ must be solved for at every time step.
That needs an equation "something with $p$ = something with $\mathbf{u}$", and we have none: the momentum equation contains a second unknown, $\partial \mathbf{u}/\partial t$, and (C) tells us nothing about $p$. We have to build one by combining the two.

The constraint at every instant

The left side: the pile-up meter must not move

If the divergence is zero now and must still be zero a moment later, it must not change:

$\begin{aligned} \frac{\partial}{\partial t}\left(\nabla\cdot\mathbf{u}\right) &= \frac{\partial}{\partial t}\left(\frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} + \frac{\partial u_z}{\partial z}\right) \\ &= \frac{\partial}{\partial t}\left(\frac{\partial u_x}{\partial x}\right) + \frac{\partial}{\partial t}\left(\frac{\partial u_y}{\partial y}\right) + \frac{\partial}{\partial t}\left(\frac{\partial u_z}{\partial z}\right) = 0 \end{aligned}$(10)

Where $\frac{\partial}{\partial t}\left(\frac{\partial u_x}{\partial x}\right)$ denotes how fast the slope $\partial u_x/\partial x$ changes at a fixed place (1/s²), likewise for the other terms.

Equation (10) says that at every fixed place, the net outflow rate does not change in time: the pile-up meter must not move. It starts at zero, because the initial $\mathbf{u}$ satisfies (C), so it stays at zero.

The update must not pile anything up

$\partial \mathbf{u}/\partial t = (\partial u_x/\partial t, \partial u_y/\partial t, \partial u_z/\partial t)$ is itself a vector field, an update arrow at every place: exactly the update in the Euler step

$\mathbf{u}^{n+1} = \mathbf{u}^n + \Delta t\,\frac{\partial \mathbf{u}}{\partial t}$(11)

Where $\mathbf{u}^n$ denotes the velocity field at step $n$, $\mathbf{u}^{n+1}$ the field one step later, and $\partial \mathbf{u}/\partial t$ the rate computed from $\mathbf{u}^n$ with Equation (8). Its divergence, $\nabla\cdot(\partial \mathbf{u}/\partial t) = \frac{\partial}{\partial x}\left(\frac{\partial u_x}{\partial t}\right) + \frac{\partial}{\partial y}\left(\frac{\partial u_y}{\partial t}\right) + \frac{\partial}{\partial z}\left(\frac{\partial u_z}{\partial t}\right)$, asks: would the update alone pile fluid up or drain it?

Take the divergence of the Euler step. The slope of a sum is the sum of the slopes and the constant $\Delta t$ pulls out. Adding the three parts:

$\nabla\cdot\mathbf{u}^{n+1} = \underbrace{\nabla\cdot\mathbf{u}^n}_{=\,0} + \Delta t\;\nabla\cdot\left(\frac{\partial \mathbf{u}}{\partial t}\right)$(12)

Where $\nabla\cdot\mathbf{u}^n = 0$ because the current field satisfies (C). So the new field passes (C) only if the update's divergence is zero. Adding a change that piles nothing up to a field that piles nothing up gives a field that piles nothing up.

Why both sides say the same thing

$\frac{\partial}{\partial t}\left(\nabla\cdot\mathbf{u}\right) = 0 \quad\iff\quad \nabla\cdot\left(\frac{\partial \mathbf{u}}{\partial t}\right) = 0$(13)

Where $\iff$ denotes "is equivalent to". Don't read $\partial u_x/\partial t$ only as "the rate of change of the velocity at a point". At a face of a tiny box, $u_x$ times the face's area is the volume flowing through that face per second (Post 11), so $\partial u_x/\partial t$ there is the rate of change of the inflow through that face. The left side says that the rate of change of the box's net outflow must stay zero, so if the net outflow was zero everywhere initially, it stays that way. The right side says the same through the faces: the rates of change of the inflows and outflows at the box's faces must balance too. In 1D, the change of the inflow at the left face equals the change of the outflow at the right face. In 3D, a single pair of faces may be unbalanced; only the sum over all six faces must stay zero.

Eliminating the update

We were stuck with one equation and two unknowns, the update $\partial \mathbf{u}/\partial t$ and the pressure $p$. Equation (13) is a second fact about the update, though not the update itself: its divergence is zero. We substitute: Equation (8) says what the update $\partial \mathbf{u}/\partial t$ equals, so take its divergence and set it to zero:

$\nabla\cdot\left(-(\mathbf{u}\cdot\nabla)\mathbf{u} - \frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u}\right) = 0$(14)

The only unknown left is $p$; everything else comes from the current snapshot of $\mathbf{u}$. Choose $p$ so that the resulting update piles up fluid nowhere. The result will not be a literal "$p$ = formula", but "(slopes of $p$) = (things computed from $\mathbf{u}$)", which is still enough to find $p$.

Three shares of the pile-up

The divergence is a sum of slopes, and the slope of a sum is the sum of the slopes. For any vector fields $\mathbf{a}$, $\mathbf{b}$, $\mathbf{c}$, the x-part is $\frac{\partial}{\partial x}(a_x + b_x + c_x) = \frac{\partial a_x}{\partial x} + \frac{\partial b_x}{\partial x} + \frac{\partial c_x}{\partial x}$, where $a_x$, $b_x$, $c_x$ denote their x-components, so $\nabla\cdot(\mathbf{a} + \mathbf{b} + \mathbf{c}) = \nabla\cdot\mathbf{a} + \nabla\cdot\mathbf{b} + \nabla\cdot\mathbf{c}$. Constants such as $-1$, $1/\rho$ and $\nu$ pull out. Equation (14) splits into three shares:

$\underbrace{-\nabla\cdot\big[(\mathbf{u}\cdot\nabla)\mathbf{u}\big]}_{\text{delivery's share}} \; \underbrace{-\;\frac{1}{\rho}\,\nabla\cdot\left(\nabla p\right)}_{\text{pressure's share}} \; + \underbrace{\nu\,\nabla\cdot\left(\nabla^2\mathbf{u}\right)}_{\text{friction's share}} = 0$(15)

Where each share is one number per place in 1/s²: the pile-up that this part of the update would cause. The pile-up caused by delivery, by pressure and by friction must add up to zero at every place.

In 2D, $p(x, y)$ is one number per place with no direction: a landscape, or a grayscale image, with hills of high and valleys of low pressure. You can ask it two different questions.

1. Which way is uphill, and how steep? The answer is the gradient $\nabla p = (\partial p/\partial x, \partial p/\partial y)$, an arrow at every place. It is what pushes: $-\frac{1}{\rho}\nabla p$ points downhill.

2. Am I a peak, a dip, or neither? The answer is the Laplacian, one number per place; on a grid, the stencil of Equation (4) applied to $p$:

$\nabla^2 p = \frac{\partial^2 p}{\partial x^2} + \frac{\partial^2 p}{\partial y^2} \approx \frac{p_{\text{right}} + p_{\text{left}} + p_{\text{up}} + p_{\text{down}} - 4\,p_{\text{here}}}{h^2}$(16)

Where $p_{\text{here}}$ denotes the pressure at a grid point and $p_{\text{right}}$, $p_{\text{left}}$, $p_{\text{up}}$, $p_{\text{down}}$ the pressures at its four neighbours (Pa). On a hilltop, $\nabla^2 p$ is negative; in a bowl, positive; on flat ground or an evenly tilted slope, zero. The tilted slope is the key case: the gradient is not zero, so there is a push, but it pushes all the water equally hard in the same direction, which moves water but piles it up nowhere. So the Laplacian measures exactly the "peak or dip" part of the push, the part that sends water away from a point or toward it. Pressure's share contains exactly this quantity:

$\nabla\cdot\left(\nabla p\right) = \frac{\partial}{\partial x}\left(\frac{\partial p}{\partial x}\right) + \frac{\partial}{\partial y}\left(\frac{\partial p}{\partial y}\right) + \frac{\partial}{\partial z}\left(\frac{\partial p}{\partial z}\right) = \nabla^2 p$(17)

Where the z-term is added for 3D. The gradient tells in which direction the pressure rises, and so which way a parcel accelerates: downhill, because of the minus sign, and per kilogram, because of $1/\rho$. The divergence of the gradient (not "the gradient of the gradient", which would be a 3×3 table like $\nabla\mathbf{u}$; the divergence keeps only its diagonal, summed) tells whether $p$ here is above, below or equal to its neighbours' average. Around a pressure peak, all push arrows point away from the peak: a positive net outflow. So "$p$ is above its neighbours' average" and "the pressure push drives fluid out of this box" are the same statement, once in curvatures and once in outflow. This is the brake from piece (C), now as a formula.

Drag the box onto the hilltop: the push on the parcel right on top is zero, yet every face of the box pushes water out.

$ hill: drag the box onto the peak

// push vs pile-up

lighter = higher p→ push −(1/ρ)∇p→ out / → in across a facep = p₀ − 500 (x² + y²) Paρ = 1000 kg/m³

push at centre −(1/ρ)∇p

(+1.60, +0.85) m/s²

∇²p

−2000 Pa/m²

net outflow −(1/ρ)∇²p

+2.00 /s²

the box

drains

box at (1.60, 0.85) m

$\nu\nabla^2\mathbf{u}$ is the vector field that tells in which direction friction accelerates the water at each place ($\nabla^2\mathbf{u}$ alone points the same way, without the strength $\nu$). Its divergence checks how these accelerations change around a point and sums the changes: it tells whether, because of acceleration differences around the point, the box there would start to gain or lose fluid. Only own-direction changes count (the x-part along $x$, and so on), because only pushes through a face move fluid in or out. And "would start to" because it concerns the update, exactly what would break (C) one time step later. In one line: $\nabla\cdot(\nu\nabla^2\mathbf{u})$ asks, is friction pushing more fluid out of this tiny box than into it?

$ drag the box: each lane is braked or pushed as a whole

// friction's share

~ flow u→ friction's push ν∇²uthrough a face: → more out · → more inwater, ν = 1 mm²/s · box 3 × 3 mm · arrows scaled

friction at the centre

0.64 mm/s²

left + right faces

0.00 /s²

top + bottom faces

0.00 /s²

pile-up change ∇·(ν∇²u)

0.000 /s²

The answer is always no. Treat the velocity field on a 2D grid as an image with two channels, $u_x$ and $u_y$. Friction is a blur-like filter: Equation (4) computes "neighbours minus me" with fixed weights, the same kernel on both channels:

$K_{\text{fric}} = \frac{1}{h^2}\begin{bmatrix} 0 & 1 & 0 \\ 1 & -4 & 1 \\ 0 & 1 & 0 \end{bmatrix}$(18)

Where $K_{\text{fric}}$ denotes the friction kernel, laid out like the grid around a pixel. Divergence is an edge-detector-like filter: $(u_{x,\text{right}} - u_{x,\text{left}})/(2h) + (u_{y,\text{up}} - u_{y,\text{down}})/(2h)$, take fixed-weight differences and sum both channels into one number per pixel, the pile-up map.

Both are filters in a CNN sense, the same small kernel slid over the whole grid, and for such filters the order doesn't matter: blurring then edge-detecting equals edge-detecting then blurring. Either way, each output pixel is a weighted sum of nearby input pixels, with each weight a product of one blur weight and one edge-detector weight; both orders give the same products at the same offsets. This is the commutativity of convolution. So

$\nabla\cdot\left(\nabla^2\mathbf{u}\right) = \nabla^2\big(\underbrace{\nabla\cdot\mathbf{u}}_{=\,0}\big) = 0$(19)

Where the left side blurs first, the middle term takes the divergence first. The pile-up map $\nabla\cdot\mathbf{u}$ is all zeros by (C), and blurring an all-zero map gives all zeros. Physically, friction only mixes existing velocities, and mixing velocities that pile up nothing cannot create a pile-up, so pressure never has to fight friction.

$-\nabla\cdot\big[(\mathbf{u}\cdot\nabla)\mathbf{u}\big]$ is the pile-up that the delivery part of the update would cause if it acted alone: positive drains the box, negative overfills it. It is built only from $\mathbf{u}$ and its slopes, and therefore is entirely computable from the current snapshot.

Water can't pile up, so shouldn't this be zero everywhere? No: $\nabla\cdot\mathbf{u}$, the divergence of the velocity, is always zero, but this is the divergence of only one part of the update, a hypothetical. Delivery never acts alone; only the total update must pile nothing up:

$\underbrace{-\nabla\cdot\big[(\mathbf{u}\cdot\nabla)\mathbf{u}\big]}_{\text{delivery's share}} + \underbrace{\left(-\frac{1}{\rho}\nabla^2 p\right)}_{\text{pressure's share}} + \underbrace{0}_{\text{friction's share}} = 0$(20)

Where pressure's share uses Equation (17) and friction's share is Equation (19). It works like a bank account that must stay at zero: income and expenses can each be large, as long as they cancel.

Worked example A channel carries water along $+x$ straight into a wall, where it splits into an upper and a lower arm. The pressure equation makes the pressure bulge up in front of the wall exactly as much as needed, and the gradient of that pressure pushes the water back and sideways, into the arms instead of piling up. The water also can't simply slow down along $x$, which alone would pile it up: as $u_x$ drops toward the wall, $u_y$ grows into the arms at exactly the matching rate. The pressure peak causes that turn, pushing back toward $-x$ and sideways toward $\pm y$.

$ drag the blue marker, press release; hover to probe

// t-junction

■ with pressure · ▢ without (delivery only)→ velocity · → pressure push −(1/ρ)∇p█ pressure, brighter = higher · contours every 0.04 Pawater at 2 cm/s, frictionless · computing…

pressure at the marker

+0.000 Pa

u_x · u_y at the marker

2.00 · 0.00 cm/s

area · with pressure

—

area · without pressure

—

The pressure Poisson equation

Equation (20) is $-\nabla\cdot\big[(\mathbf{u}\cdot\nabla)\mathbf{u}\big] - \frac{1}{\rho}\nabla^2 p = 0$. Bring delivery's share to the right side:

$-\frac{1}{\rho}\nabla^2 p = \nabla\cdot\big[(\mathbf{u}\cdot\nabla)\mathbf{u}\big]$(21)

Multiply both sides by $-\rho$:

$\boxed{\nabla^2 p = -\rho\,\nabla\cdot\big[(\mathbf{u}\cdot\nabla)\mathbf{u}\big]}$(22)

Where $\nabla^2 p$ denotes the Laplacian of the unknown pressure (Pa/m²), and $\nabla\cdot\big[(\mathbf{u}\cdot\nabla)\mathbf{u}\big]$ the divergence of the convective term, computed from the current velocity (1/s²). This is the pressure Poisson equation. At every place, the pressure must bulge up or dip down relative to its neighbours by exactly the amount needed to cancel the pile-up that delivery would cause. Where delivery would overfill a box, $p$ sits above its neighbours' average and pushes water away; where delivery would drain a box, $p$ sits below and draws water in. It isn't literally "$p$ = …", because the left side only says how $p$ compares with its neighbours, and those are unknown too.

From the equation to numbers

One equation per grid point

With the stencil of Equation (16), the pressure Poisson equation at each grid point reads

$\frac{p_{\text{right}} + p_{\text{left}} + p_{\text{up}} + p_{\text{down}} - 4\,p_{\text{here}}}{h^2} = -\rho\,\nabla\cdot\big[(\mathbf{u}\cdot\nabla)\mathbf{u}\big] =: s$(23)

Where $s$ denotes the source: one known number per grid point (Pa/m²), computed from the current $\mathbf{u}$. The pressure at each place must rise or fall with the pile-up, exactly enough that the water can't pile up. In a tiny 1D example with grid points 0 to 4, it gives the equation at point $i$:

$\frac{p_{i+1} + p_{i-1} - 2\,p_i}{h^2} = s_i$(24)

Where $p_i$ denotes the pressure at grid point $i$ and $s_i$ the known source there.

The two end points sit at walls, and we simply take their pressures $p_0$ and $p_4$ as given. That leaves three unknowns, $p_1$, $p_2$ and $p_3$. Writing Equation (24) at each of these points, multiplied by $h^2$, gives one equation per unknown:

$\begin{aligned} p_2 + p_0 - 2\,p_1 &= h^2 s_1 \\ p_3 + p_1 - 2\,p_2 &= h^2 s_2 \\ p_4 + p_2 - 2\,p_3 &= h^2 s_3 \end{aligned}$(25)

Every equation contains more than one unknown, so none of them can be solved on its own. But these three equations for three unknowns form a linear system: solved together, they pin down all three pressures. In 2D it works the same way, just with more points.

What each equation says about its own point

Multiplying Equation (24) by $h^2$, adding $2p_i$ and subtracting $h^2 s_i$ on both sides, and dividing by 2 gives

$p_i = \underbrace{\dfrac{p_{i+1} + p_{i-1}}{2}}_{\substack{\text{(1) be like} \\ \text{your neighbours}}} \; \underbrace{-\;\dfrac{h^2}{2}\,s_i}_{\substack{\text{(2) but} \\ \text{stick out}}}$(26)

Where the symbols are as in Equation (24). Part (1) says "be like your neighbours": take the average of the two neighbouring pressures. Part (2) says "but stick out by a certain amount", with the size and direction decided by $s_i$:

$s_i$ negative: water would pile up at point $i$, so the pressure must rise above its neighbours. That is why the minus sign is there: a negative $s_i$ makes $p_i$ higher than its neighbours, and the water is pushed away.
$s_i$ positive: water would drain from point $i$, so the pressure must dip below its neighbours to draw water in.
$s_i$ zero: delivery neither piles up nor drains anything there, so $p_i$ is exactly its neighbours' average, with no peak and no dip. In 1D, the pressure runs in a straight line through that point.

The circle of conditions

Writing Equation (26) at each of the three unknown points gives three recipes:

$\begin{aligned} p_1 &= \frac{p_2 + p_0}{2} - \frac{h^2}{2}\,s_1 \\ p_2 &= \frac{p_3 + p_1}{2} - \frac{h^2}{2}\,s_2 \\ p_3 &= \frac{p_4 + p_2}{2} - \frac{h^2}{2}\,s_3 \end{aligned}$(27)

Where $p_0$ and $p_4$ are the given wall values. Each line looks like a recipe for one pressure. But look at the ingredients:

To compute $p_1$, you need $p_2$.
To compute $p_2$, you need $p_1$ and $p_3$.
To compute $p_3$, you need $p_2$.

Every recipe needs the result of another one. Whichever line you start with, it asks for a value you don't have yet. So these lines are not calculations you can do one after the other. They are conditions that the correct pressures must satisfy at the same time.

Guess, then improve

The way out: since we don't know the true values, we guess them. We plug the guesses into the right-hand sides and get new values. These are better than the guesses but not yet exact, so we repeat with the new values, again and again.

With numbers: write the stick-out part as $b_i := -\frac{h^2}{2}s_i$, how far point $i$ must stick out above its neighbours' average, so each condition reads $p_i = (p_{i+1} + p_{i-1})/2 + b_i$. Take walls $p_0 = p_4 = 0$, and let water pile up only at the middle point: $b_2 = 1$, $b_1 = b_3 = 0$. Start with all guesses at zero. Each round computes all three new values from the previous round:

Round	$p_1 = \frac{p_2 + 0}{2}$	$p_2 = \frac{p_1 + p_3}{2} + 1$	$p_3 = \frac{0 + p_2}{2}$
0 (guess)	0	0	0
1	0	1	0
2	0.5	1	0.5
3	0.5	1.5	0.5
4	0.75	1.5	0.75
5	0.75	1.75	0.75
…	→ 1	→ 2	→ 1

The values creep toward $(1, 2, 1)$. Plugging this in, all three conditions hold at once: $1 = (2 + 0)/2$, $2 = (1 + 1)/2 + 1$, $1 = (0 + 2)/2$. The result is a peak in the middle, where the water would pile up, sloping down to the walls. This procedure is called the Jacobi iteration. Real solvers are faster, but they solve the same system.

Why the neighbours rise too

The source sits only at point 2, yet $p_1$ and $p_3$ rose as well. Point 1 has no stick-out, so its condition says $p_1 = (p_2 + p_0)/2$: be exactly the average of your neighbours. One of those neighbours is raised, so $p_1$ is dragged up halfway. The same holds for $p_3$.

On a grid of 100 points with the source still only in the middle, every other point is the average of its neighbours, so the pressure runs in straight lines: a tent, with its peak in the middle and straight slopes all the way to both walls. One local pile-up changes the pressure everywhere.

In the table, the change spread by one point per round. But those rounds are only the solver's guesses. No physical time passes; all of them happen within one single instant, before the next time step. Physically, the pressure adjusts everywhere instantly.

So pressure and heat both even out. The Jacobi iteration is literally a heat equation with a source, run until it stops changing. But we don't step a field forward in time. We let the pressure field settle into the one shape that satisfies every condition at once: its Laplacian is prescribed to equal the known $s$, and we adjust $p$ until it matches.

$ step the solver to find the pressure that stops the pile-up

// pressure poisson

● pressure guess p◆ neighbours' average▲ water would pile up (▼ drain), stick-out b○ next round┊ leftover = gap to nextledger: ▮ delivery + ▮ pressure = ● leftover (up = piles up)walls p = 0 · arbitrary units

solver round

0 · time t frozen

largest leftover

1.00

update of p₂

p₂ ← (0 + 0)/2 + 1.00 = 1.00 · now 0 · leftover +1.00

speed

points

The full loop

The "missing" equation for $p$ was hidden in (C) all along. (C) does not contain $p$, but it forces $p$ to take exactly these values. One time step, starting from a velocity field $\mathbf{u}^n$ that satisfies (C):

$\begin{aligned} &\text{1. compute the source:}\quad s = -\rho\,\nabla\cdot\big[(\mathbf{u}^n\cdot\nabla)\mathbf{u}^n\big] \\ &\text{2. solve for all points at once:}\quad \nabla^2 p^n = s \\ &\text{3. compute the rate:}\quad \frac{\partial \mathbf{u}}{\partial t} = -(\mathbf{u}^n\cdot\nabla)\mathbf{u}^n - \frac{1}{\rho}\nabla p^n + \nu\nabla^2\mathbf{u}^n \\ &\text{4. Euler step:}\quad \mathbf{u}^{n+1} = \mathbf{u}^n + \Delta t\,\frac{\partial \mathbf{u}}{\partial t} \end{aligned}$(28)

Where $\mathbf{u}^n$ denotes the velocity field at step $n$, $s$ the source, $p^n$ the pressure at step $n$, $\partial \mathbf{u}/\partial t$ the rate at each fixed grid point and $\Delta t$ the time step; then repeat with $\mathbf{u}^{n+1}$. On a real grid, small numerical errors make $\mathbf{u}$ drift away from (C) over many steps, so practical solvers use a refined version of this loop, the projection method, covered in a later post.

The next post looks at why the convective term makes Navier–Stokes nonlinear, unlike the heat and wave equations, and at the Reynolds number, which scores the tug of war between convection and viscosity.