Fluid Flow Simulation with the Navier-Stokes Equation

Click and drag to add velocity to the simulation. The color is the pressure.

The Navier-Stokes equation is below

ρ \frac{D \vec{u}}{D t} = - \nabla p + \nabla \cdot 𝛕 + ρ \vec{g}

$ρ$ is density.
$\vec{u}$ is the velocity vector.
$p$ is the pressure.
$𝛕$ is the stress tensor.
$\vec{g}$ is any body acceleration on the particles.

For this simulation the only body acceleration $𝐠$ is the mouse dragging. We can just ignore it.

ρ \frac{D \vec{u}}{D t} = - \nabla p + \nabla \cdot 𝛕

The material derivative for a velocity field is defined below.

\frac{D}{D t} = \frac{\partial}{\partial t} + \vec{u} \cdot \nabla

So original equation expanded is below.

ρ (\frac{\partial \vec{u}}{\partial t} + \vec{u} \cdot \nabla \vec{u}) = - \nabla p + \nabla \cdot 𝛕

The stress tensor is defined as

τ = μ (\nabla \vec{u} + (\nabla \vec{u})^{⊤}) + λ (\nabla \cdot \vec{u}) I

where $I$ is the identity matrix, and $μ, λ$ are viscosity constants. Simple algebreic manipulation yields a solution for the change in time.

\frac{\partial \vec{u}}{\partial t} = - \frac{1}{ρ} \nabla p + \frac{1}{ρ} \nabla \cdot 𝛕 - \vec{u} \cdot \nabla \vec{u}

Incompressable constraint

But we still need to define what $\nabla p$ is. To do this we add the additional constraint that the fluid is incompressable, meaning the density $ρ$ is constant over space and time.

\frac{\partial ρ}{\partial t} = 0

\nabla ρ = \vec{0}

The change in density follows from the continuity equation in a velocity field.

\frac{\partial ρ}{\partial t} + \nabla \cdot (ρ \vec{u}) = 0

The time derivative is zero from our constraint, and we apply the chain rule on the right.

ρ \nabla \cdot \vec{u} + \vec{u} \cdot \nabla ρ = 0

And the gradient is zero also zero from our constraint, and we can divide out the density (assuming it isn't zero).

\nabla \cdot \vec{u} = 0

Thus our incompressability condition implies that the velocity field has zero divergence.

Pressure

Now we can try to find solve for pressure given the constraint that divergence of velocity is zero everywhere. First start with the time evolution equation with the stress tensor expanded.

\frac{\partial \vec{u}}{\partial t} = - \frac{1}{ρ} \nabla p + \frac{1}{ρ} μ \nabla \cdot (\nabla \vec{u} + (\nabla \vec{u})^{⊤}) + \frac{1}{ρ} λ \nabla \cdot (\nabla \cdot \vec{u}) I - \vec{u} \cdot \nabla \vec{u}

The $λ$ term is equal to zero since the velocity divergence is zero by the constriant.

\frac{\partial \vec{u}}{\partial t} = - \frac{1}{ρ} \nabla p + \frac{1}{ρ} μ \nabla \cdot (\nabla \vec{u} + (\nabla \vec{u})^{⊤}) - \vec{u} \cdot \nabla \vec{u}

Distribute the divergence operator.

\frac{\partial \vec{u}}{\partial t} = - \frac{1}{ρ} \nabla p + \frac{1}{ρ} μ (\nabla \cdot \nabla \vec{u} + \nabla \cdot (\nabla \vec{u})^{⊤}) + \frac{1}{ρ} λ \nabla \cdot (\nabla \cdot \vec{u}) I - \vec{u} \cdot \nabla \vec{u}

And $\nabla \cdot (\nabla \vec{u})^{⊤} = (\nabla \cdot \vec{u}) \nabla^{⊤} = 0$ because once again the divergence is defined to be zero.

\frac{\partial \vec{u}}{\partial t} = - \frac{1}{ρ} \nabla p + \frac{1}{ρ} μ \nabla \cdot \nabla \vec{u} - \vec{u} \cdot \nabla \vec{u}

Then take the divergence of both sides

\nabla \cdot \frac{\partial \vec{u}}{\partial t} = - \frac{1}{ρ} \nabla \cdot \nabla p + \frac{1}{ρ} μ \nabla \cdot \nabla \cdot \nabla \vec{u} - \nabla \cdot (\vec{u} \cdot \nabla \vec{u})

And derivatives usually commute unless it's some really weird function. So we can move the divergence inside on the left hand side. Also the middle term is zero because $\nabla \cdot \nabla \cdot \nabla \vec{u} = \nabla \cdot \nabla (\nabla \cdot \vec{u}) = 0$ .

\frac{\partial \nabla \cdot \vec{u}}{\partial t} = - \frac{1}{ρ} \nabla \cdot \nabla p - \nabla \cdot (\vec{u} \cdot \nabla \vec{u})

And once again the gradient is always zero so the left is zero.

0 = - \frac{1}{ρ} \nabla \cdot \nabla p - \nabla \cdot (\vec{u} \cdot \nabla \vec{u})

So we get a constraint on the pressure.

\nabla \cdot \nabla p = - ρ \nabla \cdot (\vec{u} \cdot \nabla \vec{u})

Uniqueness of the solution for pressure

We need to show that the formula above should give $p$ uniquelly up to a constant. Consider two solutions

\nabla \cdot \nabla p_{1} = - ρ \nabla \cdot (\vec{u} \cdot \nabla \vec{u})

\nabla \cdot \nabla p_{2} = - ρ \nabla \cdot (\vec{u} \cdot \nabla \vec{u})

One can see that

\nabla \cdot \nabla (p_{1} - p_{2}) = 0

But we can add any function $\nabla \cdot \nabla ϕ = 0$ to a solution and it's still valid. So the solution's aren't unique up to a constant. The fix is we apply a boundary conditions.

With viscosity

If $μ \geq 0$ we enforce the boundary condition by setting $\vec{u} = \vec{0}$ so that it doesn't slip on the wall.

\frac{\partial \vec{u}}{\partial t} = - \frac{1}{ρ} \nabla p + \frac{1}{ρ} μ \nabla \cdot \nabla \vec{u} - \vec{u} \cdot \nabla \vec{u}

Becomes

0 = - \frac{1}{ρ} \nabla p + \frac{1}{ρ} μ \nabla \cdot \nabla \vec{u}

Then

\nabla p = μ \nabla \cdot \nabla \vec{u}

No viscosity

If $μ = 0$ we have the boundary condition by setting $\vec{u} \cdot \vec{n} = \vec{0}$ so that the fluid doesn't flow into or out of the wall. Then dotting both sides by $\vec{n}$ gives.

\frac{\partial \vec{u}}{\partial t} \cdot \vec{n} = - \frac{1}{ρ} \nabla p \cdot \vec{n} - \vec{n} \cdot (\vec{u} \cdot \nabla \vec{u})

And the left is zero so

\nabla p \cdot \vec{n} = - ρ \vec{n} \cdot (\vec{u} \cdot \nabla \vec{u})

Let the difference between two solutions be the below. We want to show that it's a constant everywhere.

ϕ = (p_{1} - p_{2})

And we know

\nabla \cdot \nabla ϕ = 0

. And because the two solutions have the same gradient on the boundary

\nabla ϕ = \vec{0}

on the wall. Thus dotting it with the normal vector it's also zero.

\vec{n} \cdot \nabla ϕ = 0

And there's also the vector calculus identity below.

\nabla \cdot (ϕ \nabla ϕ) = \nabla ϕ \cdot \nabla ϕ + ϕ \nabla \cdot \nabla ϕ

But the right most term is zero.

\nabla \cdot (ϕ \nabla ϕ) = \nabla ϕ \cdot \nabla ϕ

Then if we integrate over the volume

\int_{Ω} d V \nabla \cdot (ϕ \nabla ϕ) = \int_{Ω} d V \nabla ϕ \cdot \nabla ϕ

But we can also apply the divergence rule

\int_{Ω} d V \nabla \cdot (ϕ \nabla ϕ) = \oint_{Ω} d S (ϕ \nabla ϕ) \cdot \vec{n} = \oint_{Ω} d S ϕ (\nabla ϕ \cdot \vec{n})

And we know $\nabla ϕ \cdot \vec{n} = 0$ on the boundary.

\int_{Ω} d V \nabla ϕ \cdot \nabla ϕ = 0

And since magnitude is always positive it must be that:

\nabla ϕ = \vec{0}

Which means it's a constant.

Discreetizing

We have the proof that it works on a vector field but in reality we'll need to discreetize it for the simulation. We will have our velocity field $u$ live on a $N \times N$ grid with spacing $Δ s$ . We then update it with two steps. We first advect it without accounting for pressure into an intermediate velocity $u^{*}$ .

u^{*} (x, y) = u (x, y) - Δ t (u \cdot \nabla u)

Here we define the divergence with a symmetric finite differences

\nabla f (x, y) = ⟨ \frac{f (x + Δ s, y) - f (x - Δ s, y)}{2 Δ s}, \frac{f (x, y + Δ s) - f (x, y - Δ s)}{2 Δ s} ⟩

We define the next velocity to be

u_{next} (x, y) = u^{*} (x, y) - Δ t \frac{1}{ρ} \nabla p (x, y)

As usual we want the divergence of it to be zero so

0 = \nabla \cdot u_{next} (x, y) = \frac{u_{next} (x + Δ s, y) - u_{next} (x - Δ s, y)}{2 Δ s} + \frac{u_{next} (x, y + Δ s) - u_{next} (x, y - Δ s)}{2 Δ s}

And if we substitute it in

0 = \frac{u^{*} (x + Δ s, y) - Δ t \frac{1}{ρ} \nabla p (x + Δ s, y) - (u^{*} (x - Δ s, y) - Δ t \frac{1}{ρ} \nabla p (x - Δ s, y))}{2 Δ s} + \frac{u^{*} (x, y + Δ s) - Δ t \frac{1}{ρ} \nabla p (x, y + Δ s) - (u^{*} (x, y - Δ s) - Δ t \frac{1}{ρ} \nabla p (x, y - Δ s))}{2 Δ s}

Then rearrange

\frac{Δ t}{ρ} \frac{\nabla p (x + Δ s, y) - \nabla p (x - Δ s, y) + \nabla p (x, y + Δ s) - \nabla p (x, y - Δ s)}{2 Δ s} = \frac{u^{*} (x + Δ s, y) - u^{*} (x - Δ s, y)}{2 Δ s} + \frac{u^{*} (x, y + Δ s) - u^{*} (x, y - Δ s)}{2 Δ s}

And the denominators don't matter.

\frac{Δ t}{ρ} \nabla p (x + Δ s, y) - \nabla p (x - Δ s, y) + \nabla p (x, y + Δ s) - \nabla p (x, y - Δ s) = u^{*} (x + Δ s, y) - u^{*} (x - Δ s, y) + u^{*} (x, y + Δ s) - u^{*} (x, y - Δ s)

Then expand even more

\frac{p (x + Δ s, y) - 2 p (x, u) + p (x - Δ s, y)}{2 Δ s} - \frac{p (x, y + Δ s) - 2 p (x, u) + p (x, y - Δ s)}{2 Δ s} = \frac{ρ}{Δ t} (u^{*} (x + Δ s, y) - u^{*} (x - Δ s, y) + u^{*} (x, y + Δ s) - u^{*} (x, y - Δ s))

Note that this is like a discrete version of the $\nabla \cdot \nabla p = - ρ \nabla \cdot (u \cdot \nabla u)$ from earlier, because $\nabla \cdot u^{*} = \nabla \cdot u + \nabla \cdot (u \cdot \nabla u) = \nabla \cdot (u \cdot \nabla u)$ . Now if you solve for pressure of the current cell it gives the below.

p (x, y) = \frac{1}{4} (p (x + Δ s, y) + p (x - Δ s, y) + p (x, y + Δ s) + p (x, y - Δ s) + \frac{2 ρ Δ s}{Δ t} (u^{*} (x + Δ s, y) - u^{*} (x - Δ s, y) + u^{*} (x, y + Δ s) - u^{*} (x, y - Δ s)))

And this is just a linear system of equations you can solve with whatever method. Its technically just a really big matrix.

Boundary

We also have to define how it handles the boundary conditions. We skipped over how the gradient operator should be handled when it goes off the grid. We focus on the case of no viscosity. So we want

u_{next} \cdot \vec{n} = 0

Then try and figure it out

u_{next, x} n_{x} + u_{next, y} n_{y} = 0

Then substitue

(u_{x}^{*} (x, y) - Δ t \frac{1}{ρ} (\nabla p (x, y))_{x}) n_{x} + (u_{y}^{*} (x, y) - Δ t \frac{1}{ρ} (\nabla p (x, y))_{y}) n_{y} = 0

And expand again

(u_{x}^{*} (x, y) - Δ t \frac{1}{ρ} \frac{p (x + Δ s, y) - p (x - Δ s, y)}{2 Δ s}) n_{x} + (u_{y}^{*} (x, y) - Δ t \frac{1}{ρ} \frac{p (x, y + Δ s) - p (x, y - Δ s)}{2 Δ s}) n_{y} = 0

In the simulation we have 4 possible normal vectors for each of the four sides. On the left and right sides

u_{x}^{*} (x, y) - Δ t \frac{1}{ρ} \frac{p (x + Δ s, y) - p (x - Δ s, y)}{2 Δ s} = 0

Then

p (x + Δ s, y) - p (x - Δ s, y) = \frac{2 ρ Δ s}{Δ t} u_{x}^{*} (x, y)

And similarily for the top and bottom sides

p (x, y + Δ s) - p (x, y - Δ s) = \frac{2 ρ Δ s}{Δ t} u_{y}^{*} (x, y)

Thus on the boundary you can solve for the pressure that's outside grid. For example $p (x - Δ x, y)$ when $N = 0$ . And for velocity you can just assume its zero outside the grid because it would be divergence free anyways since the boundary cells are constrained to be orthogonal to the normal vector.

sources

VortexShedding github for a lot of the webgl shader code