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

ρDuDt=−∇p+∇⋅τ+ρg\rho\dfrac{D\mathbf{u}}{Dt} = -\nabla p + \nabla \cdotp \mathbf{ \tau } + \rho \mathbf{g}

where


ρ\rho is density
u\mathbf{u} is the velocity vector
pp is the pressure
τ\mathbf{\tau} is the stress tensor
gg is any other acceleration

∇\nabla is the gradient, equal to (∂∂x,∂∂y,∂∂z,...)({\partial \over \partial x}, {\partial \over \partial y}, {\partial \over \partial z}, ...) for however many dimensions. This simulation is only 2 dimensional so in this case it is just (∂∂x,∂∂y)({\partial \over \partial x}, {\partial \over \partial y}).
∇⋅{\nabla \cdotp} is the divergence and is the dot product of the gradient with the vector so it is equal to to(∂∂x,∂∂y,∂∂z...)⋅(Fx,Fy,Fz,...)=∂Fx∂x+∂Fy∂y+∂Fz∂z,...({\partial \over \partial x}, {\partial \over \partial y}, {\partial \over \partial z} ...) \cdotp (F_x, F_y, F_z, ...) = {\partial F_x\over \partial x} + {\partial F_y\over \partial y} + {\partial F_z\over \partial z}, ... for however many dimensions. In 2 dimentions it's∂Fx∂x+∂Fy∂y{\partial F_x\over \partial x} + {\partial F_y\over \partial y}.
DDtD\mathbf{} \over Dt is the material derivative which is ∂∂t+u⋅∇{\partial \over \partial t} + \mathbf{u} \cdotp \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.
DuDt=∂u∂t+u⋅∇u{D\mathbf{u} \over Dt} = {\partial\mathbf{u} \over \partial t} + \mathbf{u} \cdotp \nabla \mathbf{u}
so the origional equation becomes
ρ(∂u∂t+u⋅∇u)=−∇p+∇⋅τ+ρg\rho({\partial\mathbf{u} \over \partial t} + \mathbf{u} \cdotp \nabla \mathbf{u}) = -\nabla p + \nabla \cdotp \mathbf{ \tau } + \rho \mathbf{g}
∂u∂t+u⋅∇u=−1ρ∇p+1ρ∇⋅τ+g{\partial\mathbf{u} \over \partial t} + \mathbf{u} \cdotp \nabla \mathbf{u} = - \frac{1}{\rho} \nabla p + \frac{1}{\rho} \nabla \cdotp \mathbf{ \tau } + \mathbf{g}
τ\mathbf{\tau} 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 g\mathbf{g}.
∂u∂t+u⋅∇u=−1ρ∇p{\partial\mathbf{u} \over \partial t} + \mathbf{u} \cdotp \nabla \mathbf{u} = - \frac{1}{\rho} \nabla p
∂u∂t=−u⋅∇u−1ρ∇p{\partial\mathbf{u} \over \partial t} = - \mathbf{u} \cdotp \nabla \mathbf{u} - \frac{1}{\rho} \nabla p
∂u∂t=−u⋅(∂u∂x,∂u∂y)−1ρ∇p{\partial\mathbf{u} \over \partial t} = - \mathbf{u} \cdotp ({\partial \mathbf{u} \over \partial x}, {\partial \mathbf{u} \over \partial y}) - \frac{1}{\rho} \nabla p
∂u∂t=−ux∂u∂x−uy∂u∂y−1ρ∇p{\partial\mathbf{u} \over \partial t} = - \mathbf{u}_x {\partial \mathbf{u} \over \partial x} - \mathbf{u}_y {\partial \mathbf{u} \over \partial y} - \frac{1}{\rho} \nabla p
Then expanding the vector
[∂ux∂t∂uy∂t]=−[ux∂ux∂xux∂uy∂x]−[uy∂ux∂yuy∂uy∂y]−1ρ[∂p∂x∂p∂y]\begin{bmatrix}{\partial\mathbf{u}_x \over \partial t} \\{\partial\mathbf{u}_y \over \partial t}\end{bmatrix} = -\begin{bmatrix}\mathbf{u}_x{\partial\mathbf{u}_x \over \partial x} \\\mathbf{u}_x{\partial\mathbf{u}_y \over \partial x}\end{bmatrix}-\begin{bmatrix}\mathbf{u}_y{\partial\mathbf{u}_x \over \partial y} \\\mathbf{u}_y{\partial\mathbf{u}_y \over \partial y}\end{bmatrix}-\frac{1}{\rho}\begin{bmatrix}{\partial p \over \partial x} \\{\partial p \over \partial y}\end{bmatrix}
[∂ux∂t∂uy∂t]=−[ux∂ux∂x+uy∂ux∂y+1ρ∂p∂xux∂uy∂x+uy∂uy∂y+1ρ∂p∂y]\begin{bmatrix}{\partial\mathbf{u}_x \over \partial t} \\{\partial\mathbf{u}_y \over \partial t}\end{bmatrix} = -\begin{bmatrix}\mathbf{u}_x{\partial\mathbf{u}_x \over \partial x} + \mathbf{u}_y{\partial\mathbf{u}_x \over \partial y} + \frac{1}{\rho} {\partial p \over \partial x}\\\mathbf{u}_x{\partial\mathbf{u}_y \over \partial x} + \mathbf{u}_y{\partial\mathbf{u}_y \over \partial y} + \frac{1}{\rho} {\partial p \over \partial y}\end{bmatrix}
Applying the finite distance method:
∂f(x)∂x≈f(x+Δx)−f(x)Δx{\partial f(x) \over \partial x} \approx \frac{f(x + \Delta x) - f(x)}{\Delta x}
onto ∂ux∂t{\partial \mathbf{u}_x \over \partial t} makes:
ux(x,y,t+Δt)−ux(x,y,t)Δt=−ux(x,y,t)ux(x+ϵ,y,t)−ux(x−ϵ,y,t)2ϵ−uy(x,y,t)ux(x,y+ϵ,t)−ux(x,y−ϵ,t)2ϵ−1ρp(x+ϵ,y,t)−p(x−ϵ,y,t)2ϵ\begin{split}\frac{\mathbf{u}_x(x, y, t + \Delta t) - \mathbf{u}_x(x, y, t)}{\Delta t} = &-\mathbf{u}_x(x, y, t) \frac{\mathbf{u}_x(x + \epsilon, y, t) - \mathbf{u}_x(x - \epsilon, y, t)}{2 \epsilon} \\&-\mathbf{u}_y(x, y, t) \frac{\mathbf{u}_x(x, y + \epsilon, t) - \mathbf{u}_x(x, y - \epsilon, t)}{2 \epsilon} \\&-\frac{1}{\rho} \frac{p(x + \epsilon, y, t) - p(x - \epsilon, y, t)}{2 \epsilon}\end{split}

Then isolate ux(x,y,t+Δt)\mathbf{u}_x(x, y, t + \Delta t):

ux(x,y,t+Δt)=ux(x,y,t)−Δtux(x,y,t)ux(x+ϵ,y,t)−ux(x−ϵ,y,t)2ϵ−Δtuy(x,y,t)ux(x,y+ϵ,t)−ux(x,y−ϵ,t)2ϵ−Δt1ρp(x+ϵ,y,t)−p(x−ϵ,y,t)2ϵ\begin{split}\mathbf{u}_x(x, y, t + \Delta t) = &\mathbf{u}_x(x, y, t)\\&-\Delta t \mathbf{u}_x(x, y, t) \frac{\mathbf{u}_x(x + \epsilon, y, t) - \mathbf{u}_x(x - \epsilon, y, t)}{2 \epsilon} \\&-\Delta t \mathbf{u}_y(x, y, t) \frac{\mathbf{u}_x(x, y + \epsilon, t) - \mathbf{u}_x(x, y - \epsilon, t)}{2 \epsilon} \\&-\Delta t \frac{1}{\rho} \frac{p(x + \epsilon, y, t) - p(x - \epsilon, y, t)}{2 \epsilon}\end{split}
The first 3 terms are just the advected velocity, because its equal to ux−Δtu⋅∇ux\mathbf{u}_x - \Delta t \mathbf{u} \cdotp \nabla\mathbf{u}_x. Make a new variable a(x,y,t)a(x, y, t) for the first 3 terms, the advected velocity.ux(x,y,t+Δt)=ax(x,y,t)−Δt1ρp(x+ϵ,y,t)−p(x−ϵ,y,t)2ϵ\begin{split}\mathbf{u}_x(x, y, t + \Delta t) = &a_x(x, y, t)\\&-\Delta t \frac{1}{\rho} \frac{p(x + \epsilon, y, t) - p(x - \epsilon, y, t)}{2 \epsilon}\end{split}In incompressable fluids, density is constant so ρ\rho is a constant. If ρ\rho is constant then dρdt=0{d\rho \over dt} = 0 so the conservation of mass equation:
dρdt+∇⋅(ρu)=0{d\rho \over dt} + \nabla \cdotp (\rho \mathbf{u}) = 0
becomes
∇⋅(ρu)=0\nabla \cdotp (\rho \mathbf{u}) = 0
and since ρ\rho is a constant it can be divided out
∇⋅u=0\nabla \cdotp \mathbf{u} = 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
∂ux∂x+∂uy∂y=0{\partial \mathbf{u}_x \over \partial x} + {\partial \mathbf{u}_y \over \partial y} = 0
Then also apply the finite difference method0=ux(x+ϵ,y,t+Δt)−ux(x−ϵ,y,t+Δt)2ϵuy(x,y+ϵ,t+Δt)−uy(x,y−ϵ,t+Δt)2ϵ\begin{split}0 = &{\mathbf{u}_x(x + \epsilon, y, t + \Delta t) - \mathbf{u}_x(x - \epsilon, y, t + \Delta t) \over 2 \epsilon} \\ &{\mathbf{u}_y(x, y + \epsilon, t + \Delta t) - \mathbf{u}_y(x, y - \epsilon, t + \Delta t) \over 2 \epsilon}\end{split}We can substitute the previous equation for the 2 derivatives.
0=12ϵ((ax(x+ϵ,y,t)−Δt1ρp(x+2ϵ,y,t)−p(x,y,t)2ϵ)−(ax(x−ϵ,y,t)−Δt1ρp(x,y,t)−p(x−2ϵ,y,t)2ϵ)+(ay(x,y+ϵ,t)−Δt1ρp(x,y+2ϵ,t)−p(x,y,t)2ϵ)−(ay(x,y−ϵ,t)−Δt1ρp(x,y,t)−p(x,y−2ϵ,t)2ϵ))\begin{split}0 = \frac{1}{2\epsilon}\Bigg( & \bigg( a_x(x + \epsilon, y, t) - \Delta t \frac{1}{\rho} \frac{p(x + 2\epsilon, y, t) - p(x, y, t)}{2 \epsilon} \bigg) \\ - & \bigg( a_x(x - \epsilon, y, t) - \Delta t \frac{1}{\rho} \frac{p(x, y, t) - p(x - 2\epsilon, y, t)}{2 \epsilon} \bigg) \\ + & \bigg( a_y(x, y + \epsilon, t) - \Delta t \frac{1}{\rho} \frac{p(x, y + 2\epsilon, t) - p(x, y, t)}{2 \epsilon} \bigg) \\ - & \bigg( a_y(x, y - \epsilon, t) - \Delta t \frac{1}{\rho} \frac{p(x, y, t) - p(x, y - 2\epsilon, t)}{2 \epsilon} \bigg) \Bigg)\end{split}Then move the knows and unknowns to the left and right sides.ax(x+ϵ,y,t)−ax(x−ϵ,y,t)+ay(x,y+ϵ,t)−ay(x,y−ϵ,t)=−(−Δt1ρp(x+2ϵ,y,t)−p(x,y,t)2ϵ)+Δt1ρp(x,y,t)−p(x−2ϵ,y,t)2ϵ)−Δt1ρp(x,y+2ϵ,t)−p(x,y,t)2ϵ)+Δt1ρp(x,y,t)−p(x,y−2ϵ,t)2ϵ))\begin{split} & a_x(x + \epsilon, y, t) - a_x(x - \epsilon, y, t) \\+ & a_y(x, y + \epsilon, t) - a_y(x, y - \epsilon, t)\end{split} =\begin{split}-\Bigg( - & \Delta t \frac{1}{\rho} \frac{p(x + 2\epsilon, y, t) - p(x, y, t)}{2 \epsilon} \bigg) \\ + & \Delta t \frac{1}{\rho} \frac{p(x, y, t) - p(x - 2\epsilon, y, t)}{2 \epsilon} \bigg) \\ - & \Delta t \frac{1}{\rho} \frac{p(x, y + 2\epsilon, t) - p(x, y, t)}{2 \epsilon} \bigg) \\ + & \Delta t \frac{1}{\rho} \frac{p(x, y, t) - p(x, y - 2\epsilon, t)}{2 \epsilon} \bigg) \Bigg)\end{split}Then simplifyax(x+ϵ,y,t)−ax(x−ϵ,y,t)+ay(x,y+ϵ,t)−ay(x,y−ϵ,t)=−Δt2ϵρ(4p(x,y,t)−p(x+2ϵ,y,t)−p(x−2ϵ,y,t))−p(x,y+2ϵ,t)−p(x,y−2ϵ,t)))\begin{split} & a_x(x + \epsilon, y, t) - a_x(x - \epsilon, y, t) \\+ & a_y(x, y + \epsilon, t) - a_y(x, y - \epsilon, t)\end{split} =-\frac{\Delta t}{2\epsilon \rho}\begin{pmatrix} 4p(x, y, t) \\- p(x + 2\epsilon, y, t) - p(x - 2\epsilon, y, t)) \\- p(x, y + 2\epsilon, t) - p(x, y - 2\epsilon, t))\end{pmatrix}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 −Δt2ϵρb(x,y,t)-\frac{\Delta t}{2\epsilon \rho} b(x, y, t) to the left hand side of the equation.
p(x,y)=p(x+ϵ,y)+p(x−ϵ,y)+p(x,y+ϵ)+p(x,y−ϵ)−b(x,y)4p(x, y) = \dfrac{p(x+\epsilon, y) + p(x-\epsilon, y) + p(x, y+\epsilon) + p(x, y-\epsilon) - b(x, y)}{4}

And with that system of equations iterate by setting the initial guess for every point p0(x,y)=0p^0(x, y) = 0. Then iterate with
pn+1(x,y)=pn(x+ϵ,y)+pn(x−ϵ,y)+pn(x,y+ϵ)+pn(x,y−ϵ)−b(x,y)4p_{n+1}(x, y) = \dfrac{p_n(x+\epsilon, y) + p_n(x-\epsilon, y) + p_n(x, y+\epsilon) + p_n(x, y-\epsilon) - 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
cn+1(x,y)=cn(x,y)+k∗(an(x−ϵ,y)+cn+1(x+ϵ,y)+an+1(x,y−ϵ)+an+1(x,y+ϵ)4−cn(x,y))c_{n+1}(x, y) = c_n(x, y) + k*\big(\dfrac{a_n(x-\epsilon, y) + c_{n+1}(x+\epsilon, y) + a_{n+1}(x, y-\epsilon) + a_{n+1}(x, y+\epsilon)}{4} - c_n(x, y)\big)
between the new and the current attrubute,aa and the percent it goes towards it in that time tick,kk
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
df(x)dx=limh→0f(x+h)−f(x)h{df(x) \over dx} = lim_{h \to 0}{f(x+h)-f(x) \over h}

which can also be

df(x)dx=limh→0f(x+h)−f(x−0)h−(−0){df(x) \over dx} = lim_{h \to 0}{f(x+h)-f(x-0) \over h-(-0)}

then to

df(x)dx=limh→0f(x+h)−f(x−h)h−(−h){df(x) \over dx} = lim_{h \to 0}{f(x+h)-f(x-h) \over h-(-h)}

df(x)dx=limh→0f(x+h)−f(x−h)2h{df(x) \over dx} = lim_{h \to 0}{f(x+h)-f(x-h) \over 2h}

and in code form with the voxels where the position must be integers the most accurate would be would be
df(x)dx=f(x+ϵ)−f(x−ϵ)2{df(x) \over dx} = {f(x+\epsilon)-f(x-\epsilon) \over 2}


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