Fluid Flow Simulation with the Navier-Stokes Equation

Tap and drag to add turbulence to the simulation.

Wikipeida says that the Navier-Stokes equation is

ρ \frac{D 𝐮}{D t} = - \nabla p + \nabla \cdot 𝛕 + ρ 𝐠

where

$ρ$ is density
$𝐮$ is the velocity vector
$p$ is the pressure
$𝛕$ is the stress tensor
$g$ is any other acceleration

\nabla

is the gradient, equal to

(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z}, . . .)

for however many dimensions. This simulation is only 2 dimensional so in this case it is just

(\frac{\partial}{\partial x}, \frac{\partial}{\partial y})

\nabla \cdot

is the divergence and is the dot product of the gradient with the vector so it is equal to to

(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z} . . .) \cdot (F_{x}, F_{y}, F_{z}, . . .) = \frac{\partial F_{x}}{\partial x} + \frac{\partial F_{y}}{\partial y} + \frac{\partial F_{z}}{\partial z}, . . .

for however many dimensions. In 2 dimentions it's

\frac{\partial F_{x}}{\partial x} + \frac{\partial F_{y}}{\partial y}

\frac{D}{D t}

is the material derivative which is

\frac{\partial}{\partial t} + 𝐮 \cdot \nabla

, which measures how much material will change measured at a point moving with the fluid.
Here it is taking the matieral derivative of the velocity, or saying how much a particles velocity will change as it flows through the fluid.

\frac{D 𝐮}{D t} = \frac{\partial 𝐮}{\partial t} + 𝐮 \cdot \nabla 𝐮

so the origional equation becomes

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

\frac{\partial 𝐮}{\partial t} + 𝐮 \cdot \nabla 𝐮 = - \frac{1}{ρ} \nabla p + \frac{1}{ρ} \nabla \cdot 𝛕 + 𝐠

𝛕

doesn't matter here since it is the stress tensor used for viscosity and when the shape of the box is changing size and volume, but both of those are ignored, and we can also don't have any external forces like gravity or electromagnetism so we can ignore

𝐠

\frac{\partial 𝐮}{\partial t} + 𝐮 \cdot \nabla 𝐮 = - \frac{1}{ρ} \nabla p

\frac{\partial 𝐮}{\partial t} = - 𝐮 \cdot \nabla 𝐮 - \frac{1}{ρ} \nabla p

\frac{\partial 𝐮}{\partial t} = - 𝐮 \cdot (\frac{\partial 𝐮}{\partial x}, \frac{\partial 𝐮}{\partial y}) - \frac{1}{ρ} \nabla p

\frac{\partial 𝐮}{\partial t} = - 𝐮_{x} \frac{\partial 𝐮}{\partial x} - 𝐮_{y} \frac{\partial 𝐮}{\partial y} - \frac{1}{ρ} \nabla p

Then expanding the vector

[\begin{array}{c} \frac{\partial 𝐮_{x}}{\partial t} \\ \frac{\partial 𝐮_{y}}{\partial t} \end{array}] = - [\begin{array}{c} 𝐮_{x} \frac{\partial 𝐮_{x}}{\partial x} \\ 𝐮_{x} \frac{\partial 𝐮_{y}}{\partial x} \end{array}] - [\begin{array}{c} 𝐮_{y} \frac{\partial 𝐮_{x}}{\partial y} \\ 𝐮_{y} \frac{\partial 𝐮_{y}}{\partial y} \end{array}] - \frac{1}{ρ} [\begin{array}{c} \frac{\partial p}{\partial x} \\ \frac{\partial p}{\partial y} \end{array}]

[\begin{array}{c} \frac{\partial 𝐮_{x}}{\partial t} \\ \frac{\partial 𝐮_{y}}{\partial t} \end{array}] = - [\begin{array}{c} 𝐮_{x} \frac{\partial 𝐮_{x}}{\partial x} + 𝐮_{y} \frac{\partial 𝐮_{x}}{\partial y} + \frac{1}{ρ} \frac{\partial p}{\partial x} \\ 𝐮_{x} \frac{\partial 𝐮_{y}}{\partial x} + 𝐮_{y} \frac{\partial 𝐮_{y}}{\partial y} + \frac{1}{ρ} \frac{\partial p}{\partial y} \end{array}]

Applying the finite distance method:

\frac{\partial f (x)}{\partial x} \approx \frac{f (x + Δ x) - f (x)}{Δ x}

onto $\frac{\partial 𝐮_{x}}{\partial t}$ makes:

\begin{aligned} \frac{𝐮_{x} (x, y, t + Δ t) - 𝐮_{x} (x, y, t)}{Δ t} = & - 𝐮_{x} (x, y, t) \frac{𝐮_{x} (x + ϵ, y, t) - 𝐮_{x} (x - ϵ, y, t)}{2 ϵ} \\ - 𝐮_{y} (x, y, t) \frac{𝐮_{x} (x, y + ϵ, t) - 𝐮_{x} (x, y - ϵ, t)}{2 ϵ} \\ - \frac{1}{ρ} \frac{p (x + ϵ, y, t) - p (x - ϵ, y, t)}{2 ϵ} \end{aligned}

Then isolate $𝐮_{x} (x, y, t + Δ t)$ :

\begin{aligned} 𝐮_{x} (x, y, t + Δ t) = & 𝐮_{x} (x, y, t) \\ - Δ t 𝐮_{x} (x, y, t) \frac{𝐮_{x} (x + ϵ, y, t) - 𝐮_{x} (x - ϵ, y, t)}{2 ϵ} \\ - Δ t 𝐮_{y} (x, y, t) \frac{𝐮_{x} (x, y + ϵ, t) - 𝐮_{x} (x, y - ϵ, t)}{2 ϵ} \\ - Δ t \frac{1}{ρ} \frac{p (x + ϵ, y, t) - p (x - ϵ, y, t)}{2 ϵ} \end{aligned}

The first 3 terms are just the advected velocity, because its equal to

𝐮_{x} - Δ t 𝐮 \cdot \nabla 𝐮_{x}

. Make a new variable

a (x, y, t)

for the first 3 terms, the advected velocity.

\begin{aligned} 𝐮_{x} (x, y, t + Δ t) = & a_{x} (x, y, t) \\ - Δ t \frac{1}{ρ} \frac{p (x + ϵ, y, t) - p (x - ϵ, y, t)}{2 ϵ} \end{aligned}

In incompressable fluids, density is constant so

ρ

is a constant. If

ρ

is constant then

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

so the conservation of mass equation:

\frac{d ρ}{d t} + \nabla \cdot (ρ 𝐮) = 0

becomes

\nabla \cdot (ρ 𝐮) = 0

and since $ρ$ is a constant it can be divided out

\nabla \cdot 𝐮 = 0

One point can't have 2 neighboring points both pulling from it or else that means it is going in 2 directions and creating new mass, or if a cell has 2 neighbors going into it mass is deleted since density can't change, so the divergence of the velocity at every point must equal zero.

So the solution must also have zero divergence

\frac{\partial 𝐮_{x}}{\partial x} + \frac{\partial 𝐮_{y}}{\partial y} = 0

Then also apply the finite difference method

\begin{aligned} 0 = & \frac{𝐮_{x} (x + ϵ, y, t + Δ t) - 𝐮_{x} (x - ϵ, y, t + Δ t)}{2 ϵ} \\ \frac{𝐮_{y} (x, y + ϵ, t + Δ t) - 𝐮_{y} (x, y - ϵ, t + Δ t)}{2 ϵ} \end{aligned}

We can substitute the previous equation for the 2 derivatives.

\begin{aligned} 0 = \frac{1}{2 ϵ} ( & (a_{x} (x + ϵ, y, t) - Δ t \frac{1}{ρ} \frac{p (x + 2 ϵ, y, t) - p (x, y, t)}{2 ϵ}) \\ - & (a_{x} (x - ϵ, y, t) - Δ t \frac{1}{ρ} \frac{p (x, y, t) - p (x - 2 ϵ, y, t)}{2 ϵ}) \\ + & (a_{y} (x, y + ϵ, t) - Δ t \frac{1}{ρ} \frac{p (x, y + 2 ϵ, t) - p (x, y, t)}{2 ϵ}) \\ - & (a_{y} (x, y - ϵ, t) - Δ t \frac{1}{ρ} \frac{p (x, y, t) - p (x, y - 2 ϵ, t)}{2 ϵ})) \end{aligned}

Then move the knows and unknowns to the left and right sides.

\begin{aligned} a_{x} (x + ϵ, y, t) - a_{x} (x - ϵ, y, t) \\ + & a_{y} (x, y + ϵ, t) - a_{y} (x, y - ϵ, t) \end{aligned} = \begin{aligned} - (- & Δ t \frac{1}{ρ} \frac{p (x + 2 ϵ, y, t) - p (x, y, t)}{2 ϵ}) \\ + & Δ t \frac{1}{ρ} \frac{p (x, y, t) - p (x - 2 ϵ, y, t)}{2 ϵ}) \\ - & Δ t \frac{1}{ρ} \frac{p (x, y + 2 ϵ, t) - p (x, y, t)}{2 ϵ}) \\ + & Δ t \frac{1}{ρ} \frac{p (x, y, t) - p (x, y - 2 ϵ, t)}{2 ϵ})) \end{aligned}

Then simplify

\begin{aligned} a_{x} (x + ϵ, y, t) - a_{x} (x - ϵ, y, t) \\ + & a_{y} (x, y + ϵ, t) - a_{y} (x, y - ϵ, t) \end{aligned} = - \frac{Δ t}{2 ϵ ρ} (\begin{array}{c} 4 p (x, y, t) \\ - p (x + 2 ϵ, y, t) - p (x - 2 ϵ, y, t)) \\ - p (x, y + 2 ϵ, t) - p (x, y - 2 ϵ, t)) \end{array})

In the code this systems of equations could be solved with a massive matrix with each voxel as a variable and have very perfect calculations, or just do it iteratevly and approximate with an iterative solver. This code uses the Gauss Seidel method. Set

- \frac{Δ t}{2 ϵ ρ} b (x, y, t)

to the left hand side of the equation.

p (x, y) = \frac{p (x + ϵ, y) + p (x - ϵ, y) + p (x, y + ϵ) + p (x, y - ϵ) - b (x, y)}{4}

And with that system of equations iterate by setting the initial guess for every point $p^{0} (x, y) = 0$ . Then iterate with

p_{n + 1} (x, y) = \frac{p_{n} (x + ϵ, y) + p_{n} (x - ϵ, y) + p_{n} (x, y + ϵ) + p_{n} (x, y - ϵ) - b (x, y)}{4}

After the divergence is cleared, the velocities just need to move both the color, themselves, and any other values in the direction they are pointing.

In code this is optimised by pulling instead from behind instead of pushing out in front because the end point of the velicity will be inbetween 4 squares and writing values is almost always more expensive than reading values, so instead it will read the 4 values of the squares if the velocity vector was rotated 180 degrees or muliplied by -1, linear interpolate the expected value at the sub-voxel position, then write it into only one square saving time and making the simulation faster at about the same accuracy as the intuitive method.

Now the hardest part is done, but there is still diffusion where quares take the average values of their neighbors over time this is modeled by

c_{n + 1} (x, y) = c_{n} (x, y) + k * (\frac{a_{n} (x - ϵ, y) + c_{n + 1} (x + ϵ, y) + a_{n + 1} (x, y - ϵ) + a_{n + 1} (x, y + ϵ)}{4} - c_{n} (x, y))

between the new and the current attrubute,

a

and the percent it goes towards it in that time tick,

k

this is again a system of equations that can easily be solved with the Gauss-Seidel method

For boundaries and walls, the velocity in that voxel must be both 0 in the y and x direction, so when divergence clearing is calculated it will make the velocity near only be able to flow parrallel to obsticles and the boundary. Since the divergence clearing isn't perfect it doesn't make velocity exacly parralel so some fluid leaks off the edge or into walls, so a script should iterate through every voxel and make the x velocity 0 if there are obstacles or boundaries to the left or right, and y velocity 0 if there are up or down, exept a more computationally efficient way is to just set the voxels which are considered boundaries to have a x and y velocity of 0, then the divergence clearing will automatically make flow parralel to them

For the numerical implementation of derivatives it is just

\frac{d f (x)}{d x} = l i m_{h \to 0} \frac{f (x + h) - f (x)}{h}

which can also be

\frac{d f (x)}{d x} = l i m_{h \to 0} \frac{f (x + h) - f (x - 0)}{h - (- 0)}

then to

\frac{d f (x)}{d x} = l i m_{h \to 0} \frac{f (x + h) - f (x - h)}{h - (- h)}

\frac{d f (x)}{d x} = l i m_{h \to 0} \frac{f (x + h) - f (x - h)}{2 h}

and in code form with the voxels where the position must be integers the most accurate would be would be

\frac{d f (x)}{d x} = \frac{f (x + ϵ) - f (x - ϵ)}{2}

used sources:
https://jamie-wong.com/2016/08/05/webgl-fluid-simulation/