back to listing index

### Fluid Simulation (with WebGL demo) - Zero Wind :: Jamie Wong

[web search]
Original source (jamie-wong.com)
Clipped on: 2016-08-27

# Fluid Simulation (with WebGL demo)

Posted on August 5, 2016

Click and drag to change the fluid flow. Double click to reset.

Note: The demos in this post rely on WebGL features that might not be implemented in mobile browsers.

About a year and a half ago, I had a passing interest in trying to figure out how to make a fluid simulation. At the time, it felt just a bit out of my reach, requiring knowledge of shaders, vector calculus, and numerical computation that were all just a little bit past my grasp. At the time, I was working through the Fluid Simulation Course Notes from SIGGRAPH 2007, and was struggling with the math. Now armed with a bit more knowledge and a lot more time, and with the help of other less dense resources like GPU Gems Chapter 38. Fast Fluid Dynamics Simulation on the GPU, I was finally able to figure out enough to get something working. I am still a beginner at simulations like this, and I’m going to brazenly ignore things like numerical stability, but hopefully I can help leapfrog you past a few places I got stuck.

We’re going to work with the simplest 2D fluid simulation, where the entire area is full of fluid, and we’re going to ignore viscosity.

# The Velocity Field

As compared to a rigid, unrotating solid, where every bit of the thing has to be moving in the same direction at the same speed, each bit of a fluid might be moving differently. One way to model this is to use a vector field representing velocity. For any given (x,y) (x, y) (x,y) coordinate, this field will tell you the velocity of the fluid at that point.

u(x,y)=(ux,uy) \vec u(x, y) = (u_x, u_y) u(x,y)=(ux,uy)

A nice way to get an intuition about what a given field looks like is to sample the function in a grid of points, then draw arrows starting at each grid point whose size and orientation are dictated by the value of the function at that point. For the purposes of this post, we’re always going to be working over the domain x[1,1] x \in [-1, 1] x[1,1], y[1,1] y \in [-1, 1] y[1,1].

For instance, here’s a very simple field u(x,y)=(1,0) \vec u(x, y) = (1, 0) u(x,y)=(1,0) representing everything moving at a constant speed to the right.

And here’s a more interesting one u(x,y)=(x,y) \vec u(x, y) = (x, y) u(x,y)=(x,y) where things move away from the origin, increasing in speed the farther away from the origin they are.

We’re going to play with this one, u(x,y)=(sin(2πy),sin(2πx)) \vec u(x, y) = \left( \sin (2 \pi y), \sin (2 \pi x) \right) u(x,y)=(sin(2πy),sin(2πx)), since it creates some interesting visual results once we start making the fluid move accordingly.

For a more thorough introduction to vector fields, check out the Introduction to Vector Fields video on Khan Academy. The rest of the videos on multivariate calculus might prove useful for understanding concepts in fluid flow too.

Now then, let’s get things moving.

Advection is the transfer of a property from one place to another due to the motion of the fluid. If you’ve got some black dye in some water, and the water is moving to the right, then surprise surprise, the black dye moves right.

If the fluid is moving in a more complex manner, that black dye will get pulled through the liquid in a more complex manner.

Before we dive into how advection works, we need to talk a bit about the format of the data underlying these simulations.

The simulation consist of two fields: color and velocity. Each field is represented by a two dimensional grid. For simplicity, we use the same dimensions as the output pixel grid.

Previously, I described the velocity field as an analytical function u(x,y) \vec u(x, y) u(x,y). In practice, that analytical function is only used to initialize the grid values.

The simulation runs by stepping forward bit-by-bit in time, with the state of the color and velocity grids depending only on the state of the color and velocity grids from the previous time step. We’ll use u(p,t) \vec u(\vec p, t) u(p,t) to represent the velocity grid at 2d position p \vec p p and time t t t, and c(p,t) \vec c(\vec p, t) c(p,t) to represent the color in the same manner.

So how do we move forward in time? Let’s just talk about how the color field changes for now. If we consider each grid point as a little particle in the fluid, then one approach is to update the color of the fluid where that particle will be, one time step in the future.

c(p+u(p,t)Δt,t+Δt):=c(p,t) \vec c(\vec p + \vec u(\vec p, t) \Delta t, t + \Delta t) := \vec c(\vec p, t) c(p+u(p,t)Δt,t+Δt):=c(p,t)

In order to run these simulations in real-time at high resolution, we want to implement them on the GPU. It turns out that this method of updating the value at the new location of the particle is difficult to implement on the GPU.

First, the position we want to write, p+u(p,t)Δt \vec p + \vec u(\vec p, t) \Delta t p+u(p,t)Δt might not lie on a grid point, so we’d have to distribute the impact of the write across the surrounding grid points. Second, many of our imaginary particles might end up in the same place, meaning we need to analyze the entire grid before we decide what the new values of each grid point might be.

So, instead of figuring out where our imaginary particles at the grid points go to, we’ll figure out where they came from in order to calculate the next time step.

c(p,t+Δt):=c(pu(p)Δt,t) \vec c(\vec p, t + \Delta t) := \vec c(\vec p - \vec u(\vec p) \Delta t, t) c(p,t+Δt):=c(pu(p)Δt,t)

With this scheme, we only need to write to a single grid point, and we don’t need to consider the contributions of imaginary particles coming from multiple different places.

The last teensy hurdle is figuring out the value of c(pu(p,t)Δt,t) \vec c(\vec p - \vec u(\vec p, t) \Delta t, t ) c(pu(p,t)Δt,t), since pu(p,t) \vec p - \vec u(\vec p, t) pu(p,t) might not be at a grid point. We can hop this hurdle using bilinear interpolation on the surrounding 4 grid points (the ones linked by the dashed grey rectangle above).

# Advecting the Velocity Field

Barring a bizarre sequence of perfectly aligned fans underneath the liquid, there’s no reason why the velocity field wouldn’t change over time. Just as black ink would move through the fluid, so too will the velocity field itself! Just as we can advect c \vec c c through u \vec u u, we can also advect u \vec u u through itself!

Intuitively you can think of it this way: a particle moving in a certain direction will continue moving in that direction, even after it’s moved.

Since we’re storing velocity in a grid just like we did with color, we can use the exact same routine to advect velocity through itself. Below, watch the velocity change over time, with an initial velocity field of u=(1,sin(2πy)) \vec u = (1, \sin(2 \pi y)) u=(1,sin(2πy)).

See how the changes you make by dragging propagate through space via advection.

If you tried playing around with this, and saw a bunch of weird hard edges and might’ve thought to yourself “I don’t think fluids work like that…”, you’d be right. We’re missing an important ingredient, but before we look at the solution, let’s take a closer look at the problem.

# Divergent Fields

Something about the velocity field below makes this intuitively not feel like a fluid. Fluids just don’t behave like this.

Same problem with this one…

If you look at where the arrows are pointing in each of the above 2 simulations, you’ll see that there are spots where the all the arrows point away from that spot, and others where all the arrows point toward that spot. Assuming the volume of the liquid is staying constant, the density of the fluid has to be changing for such a velocity field to be possible.

Water is roughly incompressible. That means that at every spot, you have to have the same amount of fluid entering that spot as leaving it.

Mathematically, we can represent this fact by saying a field is divergence-free. The divergence of a velocity field u \vec u u, indicated with div(u) div(\vec u) div(u) or u \nabla \cdot \vec u u, is a measure of how much net stuff is entering or leaving a given spot in the field. For our 2D velocity field, it’s defined like this:

u=(x,y)(ux,uy)=uxx+uyy\begin{aligned} \nabla \cdot \vec u &= \left( \frac{\partial}{\partial x}, \frac{\partial}{\partial y} \right) \cdot \left( u_x, u_y \right) \\ &= \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} \end{aligned}u=(x,y)(ux,uy)=xux+yuy

The first of the two not-very-fluidy fields above has an equation u(x,y)=(x,y) \vec u(x, y) = (x, y) u(x,y)=(x,y). Taking the divergence, we find:

u=x(x)+y(y)=1+1=2\begin{aligned} \nabla \cdot \vec u &= \frac{\partial}{\partial x}(x) + \frac{\partial}{\partial y}(y) \\ &= 1 + 1 \\ &= 2 \end{aligned}u=x(x)+y(y)=1+1=2

This positive value tells us that, in all places, more stuff is leaving that point than entering it. In physical terms, this means that the density is decreasing uniformly everywhere.

The other not-very-fluidy field has an equation u(x,y)=sin(2πx),0) \vec u(x, y) = \sin(2 \pi x), 0) u(x,y)=sin(2πx),0). If we look at its divergence, we see:

u=x(sin(2πx))+y(0)=2πcos(2πx)\begin{aligned} \nabla \cdot \vec u &= \frac{\partial}{\partial x}(\sin(2 \pi x)) + \frac{\partial}{\partial y}(0) \\ &= 2 \pi \cos (2 \pi x) \end{aligned}u=x(sin(2πx))+y(0)=2πcos(2πx)

Which tells us that in some places, density is increasing (where u<0 \nabla \cdot \vec u < 0 u<0), and in others, density is decreasing (where u>0 \nabla \cdot \vec u > 0 u>0).

Doing the same operation on the more fluidy looking swirly velocity field u=(sin(2πy),sin(2πx) \vec u = (\sin ( 2 \pi y), \sin ( 2 \pi x ) u=(sin(2πy),sin(2πx) that you saw in the section about advection, we discover u=0 \nabla \cdot \vec u = 0 u=0.

An incompressible fluid will have a divergence of zero everywhere. So, if we want our simulated fluid to look kind of like a real fluid, we better make sure it’s divergence-free.

Since our velocity field undergoes advection and can be influenced by clicking and dragging around the fluid, having an initially divergence-free velocity field isn’t enough to guarantee that the field will continue to be divergence-free. For example, if we take our swirly simulation and start advecting the velocity field through itself, we end up with something divergent:

So we need a way of taking a divergent field and making it divergence-free. To understand what force makes that happen in the real world, we need to talk about some honest-to-goodness physics.

# Navier-Stokes

The Navier-Stokes equations describe the motion of fluids. Here are the Navier-Stokes equations for incompressible fluid flow:

ut=uu1ρp+ν2u+Fu=0\begin{aligned} & \frac{\partial \vec u}{\partial t} = -\vec u \cdot \nabla \vec u -\frac{1}{\rho} \nabla p + \nu \nabla^2 \vec u + \vec F \\ \\ & \nabla \cdot \vec u = 0 \end{aligned}tu=uuρ1p+ν2u+Fu=0

Where u \vec u u is the velocity field, ρ \rho ρ is density, p p p is pressure, ν \nu ν is the kinematic viscosity, and F \vec F F is external forces acting upon the fluid.

Since we’re pretending the viscosity of our fluid is zero, we can drop the ν \nu ν term in the first equation. In our simple simulation, external forces are only applied by dragging the mouse, so we’ll ignore that term for now, opting to allow it to influence the velocity field directly.

Dropping those terms, we’re left with the following:

ut=uu1ρp \frac{\partial \vec u}{\partial t} = -\vec u \cdot \nabla \vec u -\frac{1}{\rho} \nabla p tu=uuρ1p

We can expand this to its partial derivative form, expanding vector components to leave us with only scalar variables.

[uxtuyt]=[uxxuxyuyxuyy][uxuy]1ρ[pxpy][uxtuyt]=[uxuxxuyuxy1ρpxuxuyxuyuyy1ρpy]\begin{aligned} \begin{bmatrix} \frac{\partial u_x}{\partial t} \\ \\ \frac{\partial u_y}{\partial t} \end{bmatrix} &= - \begin{bmatrix} \frac{\partial u_x}{\partial x} & \frac{\partial u_x}{\partial y} \\ \\ \frac{\partial u_y}{\partial x} & \frac{\partial u_y}{\partial y} \end{bmatrix} \begin{bmatrix} u_x \\ \\ u_y \end{bmatrix} - \frac{1}{\rho} \begin{bmatrix} \frac{\partial p}{\partial x} \\ \\ \frac{\partial p}{\partial y} \end{bmatrix} \\ \begin{bmatrix} \frac{\partial u_x}{\partial t} \\ \\ \frac{\partial u_y}{\partial t} \end{bmatrix} &= \begin{bmatrix} - u_x \frac{\partial u_x}{\partial x} - u_y \frac{\partial u_x}{\partial y} - \frac{1}{\rho} \frac{\partial p}{\partial x} \\ \\ - u_x \frac{\partial u_y}{\partial x} - u_y \frac{\partial u_y}{\partial y} - \frac{1}{\rho} \frac{\partial p}{\partial y} \end{bmatrix} \end{aligned}tuxtuytuxtuy=xuxxuyyuxyuyuxuyρ1xpyp=uxxuxuyyuxρ1xpuxxuyuyyuyρ1yp

Remembering that these fields are all functions on (x,y,t) (x, y, t) (x,y,t), we can approximate the partial derivatives with finite differences. For instance, we can approximate the partial derivative of ux u_x ux with respect to t t t like so:

uxtux(x,y,t+Δt)ux(x,y,t)Δt \frac{\partial u_x}{\partial t} \approx \frac{u_x(x, y, t + \Delta t) - u_x(x, y, t)}{\Delta t} tuxΔtux(x,y,t+Δt)ux(x,y,t)

Because the procedure ends up being the same for both components, we’ll focus on only the x x x component here. Applying finite differences to all of the partial derivatives, we have this:

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{aligned} \frac{u_x(x, y, t + \Delta t) - u_x(x, y, t)}{\Delta t} =& -u_x(x, y, t) \frac{u_x(x + \epsilon, y, t) - u_x(x - \epsilon, y, t)}{2 \epsilon} \\ & -u_y(x, y, t) \frac{u_x(x, y + \epsilon, t) - 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{aligned} Δtux(x,y,t+Δt)ux(x,y,t)=ux(x,y,t)2ϵux(x+ϵ,y,t)ux(xϵ,y,t)uy(x,y,t)2ϵux(x,y+ϵ,t)ux(x,yϵ,t)ρ12ϵp(x+ϵ,y,t)p(xϵ,y,t)

Ultimately what we want is u(x,y,t+Δt) \vec u(x, y, t + \Delta t) u(x,y,t+Δt), which will tell us, for a given point, what the velocity will be at the next time step. So let’s solve for that by rearranging the big long formula above:

ux(x,y,t+Δt)=ux(x,y,t)ux(x,y,t)Δtux(x+ϵ,y,t)ux(xϵ,y,t)2ϵuy(x,y,t)Δtux(x,y+ϵ,t)ux(x,yϵ,t)2ϵ1ρΔtp(x+ϵ,y,t)p(xϵ,y,t)2ϵ \begin{aligned} u_x(x, y, t + \Delta t) = & u_x(x, y, t) \\ & - u_x(x, y, t) \Delta t \frac{u_x(x + \epsilon, y, t) - u_x(x - \epsilon, y, t)}{2 \epsilon} \\ & -u_y(x, y, t) \Delta t \frac{u_x(x, y + \epsilon, t) - u_x(x, y - \epsilon, t)}{2 \epsilon} \\ & -\frac{1}{\rho} \Delta t \frac{p(x + \epsilon, y, t) - p(x - \epsilon, y, t)}{2 \epsilon} \end{aligned} ux(x,y,t+Δt)=ux(x,y,t)ux(x,y,t)Δt2ϵux(x+ϵ,y,t)ux(xϵ,y,t)uy(x,y,t)Δt2ϵux(x,y+ϵ,t)ux(x,yϵ,t)ρ1Δt2ϵp(x+ϵ,y,t)p(xϵ,y,t)

If you look at the first three terms in this expression, what does it look like they conceptually represent? It looks like they represent the next velocity after we’ve taken into account changes due to the motion of the fluid itself. That sounds an awful like advection as discussed earlier. In fact, it will work quite well if we substitute the velocity field after it’s undergone advection. We’ll call the advected velocity field ua \vec u ^ a ua. So now we have:

ux(x,y,t+Δt)=uxa(x,y,t)1ρΔtp(x+ϵ,y,t)p(xϵ,y,t)2ϵ \begin{aligned} u_x(x, y, t + \Delta t) = u^a_x(x, y, t) -\frac{1}{\rho} \Delta t \frac{p(x + \epsilon, y, t) - p(x - \epsilon, y, t)}{2 \epsilon} \end{aligned} ux(x,y,t+Δt)=uxa(x,y,t)ρ1Δt2ϵp(x+ϵ,y,t)p(xϵ,y,t)

So after all of that, we have an equation that relates the velocity field at the next time tick to the current velocity field after it’s undergone advection, followed by application of pressure.

We know that a divergence-free field that undergoes advection isn’t necessarily still divergence-free, and yet we know that the Navier-Stokes equations for impressible flow represent divergence-free velocity fields, so therefore we have our answer about what in nature prevents the velocity field from becoming divergent: pressure!

# Solving for Pressure

Now that we have an equation that relates u \vec u u to p p p. This is where the math gets messy. We start from the second Navier-Stokes equation for incompressible flow, applied at time t+Δt t + \Delta t t+Δt, and apply finite differences again:

u=0uxx+uyy=0ux(x+ϵ,y,t+Δt)ux(xϵ,y,t+Δt)2ϵ+uy(x,y+ϵ,t+Δt)uy(x,yϵ,t+Δt)2ϵ=0 \begin{aligned} \nabla \cdot \vec u & = 0 \\ \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} & = 0 \\ \frac{ u_x(x + \epsilon, y, t + \Delta t) - u_x(x - \epsilon, y, t + \Delta t) }{ 2 \epsilon } + \\ \\ \frac{ u_y(x, y + \epsilon, t + \Delta t) - u_y(x, y - \epsilon, t + \Delta t) }{ 2 \epsilon } & = 0 \end{aligned} uxux+yuy2ϵux(x+ϵ,y,t+Δt)ux(xϵ,y,t+Δt)+2ϵuy(x,y+ϵ,t+Δt)uy(x,yϵ,t+Δt)=0=0=0

Here, we can substitute our equations for u \vec u u expressed in terms of ua \vec u ^ a ua and p p p to get this monster:

0=12ϵ((uxa(x+ϵ,y,t)1ρΔtp(x+2ϵ,y,t)p(x,y,t)2ϵ)(uxa(xϵ,y,t)1ρΔtp(x,y,t)p(x2ϵ,y,t)2ϵ)+(uya(x,y+ϵ,t)1ρΔtp(x,y+2ϵ,t)p(x,y,t)2ϵ)(uya(x,yϵ,t)1ρΔtp(x,y,t)p(x,y2ϵ,t)2ϵ))\begin{aligned} 0 = \frac{1}{2 \epsilon} \left( \left( u^a_x(x + \epsilon, y, t) -\frac{1}{\rho} \Delta t \frac{p(x + 2\epsilon, y, t) - p(x, y, t)}{2 \epsilon} \right) \right. \\ \left. - \left( u^a_x(x - \epsilon, y, t) -\frac{1}{\rho} \Delta t \frac{p(x, y, t) - p(x - 2\epsilon, y, t)}{2 \epsilon} \right) \right. \\ \left. + \left( u^a_y(x, y + \epsilon, t) -\frac{1}{\rho} \Delta t \frac{p(x, y + 2\epsilon, t) - p(x, y, t)}{2 \epsilon} \right) \right. \\ \left. - \left( u^a_y(x, y - \epsilon, t) -\frac{1}{\rho} \Delta t \frac{p(x, y, t) - p(x, y - 2\epsilon, t)}{2 \epsilon} \right) \right) \end{aligned}0=2ϵ1((uxa(x+ϵ,y,t)ρ1Δt2ϵp(x+2ϵ,y,t)p(x,y,t))(uxa(xϵ,y,t)ρ1Δt2ϵp(x,y,t)p(x2ϵ,y,t))+(uya(x,y+ϵ,t)ρ1Δt2ϵp(x,y+2ϵ,t)p(x,y,t))(uya(x,yϵ,t)ρ1Δt2ϵp(x,y,t)p(x,y2ϵ,t)))

Rearranging to have all of the p p p terms on the left and all the ua \vec u ^ a ua terms on the right, and multiplying both sides by 2ϵ 2 \epsilon 2ϵ, we have:

Δt2ϵρ(4p(x,y,t)p(x+2ϵ,y,t)p(x2ϵ,y,t)p(x,y+2ϵ,t)p(x,y2ϵ,t))=uxa(x+ϵ,y,t)uxa(xϵ,y,t)+uya(x,y+ϵ,t)uya(x,yϵ,t) -\frac{\Delta t}{2 \epsilon \rho} \left( \begin{matrix} 4 p(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{matrix} \right) = \begin{matrix} u^a_x(x + \epsilon, y, t) \\ -u^a_x(x - \epsilon, y, t) \\ +u^a_y(x, y + \epsilon, t) \\ -u^a_y(x, y - \epsilon, t) \end{matrix} 2ϵρΔt4p(x,y,t)p(x+2ϵ,y,t)p(x2ϵ,y,t)p(x,y+2ϵ,t)p(x,y2ϵ,t)=uxa(x+ϵ,y,t)u