Interactive 2D Fluid Simulation
Henrik Peegel, Jaakob-Jaan Avvo, Ander Pavlov
Repository link: https://github.com/vapsik/Interactive_Fluid_Simulations
Introduction
The goal of this project was to create a gallery of interactive 2D fluid simulations that enact the Navier-Stokes' equations for incompressible fluids. The algorithm used to simulate the fluids is Stable Fluids based on the Jos Stam's seminal work, a 2001 article in the field of real-time fluid simulation that has been extensively used in video games since its inception. In this algorithm, a simple incompressible (constant density) fluid is represented in a rectangular grid as a field of coupled velocity vectors and pressure scalars that will dictate the movement of the dye field representing the color of the fluid or smoke being simulated. The method relies on semi-Lagrangian advection and, which all allows for stable, large time-step simulation without numerical divergence, combined with a projection step that enforces incompressibility, which makes the velocity field divergence-free by iteratively solving the Poisson's equation. Diffusion, external forces, and boundary interactions are incorporated in a computationally efficient time evolution scheme such that larger time steps can be used without losing physical accuracy or stability.
Navier-Stokes' Equations
The N-S equations are the following {$$\begin{align} & \frac{\partial \vec{u}}{\partial t} = - \frac{\nabla p}{\rho} - (\vec{u} \cdot \nabla)\vec{u} + \nu \nabla^{2} \vec{u} + \frac{\vec{F}_{\text{external}}}{\rho}, \\&\nabla \cdot \vec{u} = 0, \end{align}$$} where the terms are the following:
- {$\vec{x}$} - let that be the position vector which in 2D is just {$\vec{x} = (x,y)$}.
- {$\rho$} - constant density of the fluid,
- {$\nu$} - the viscosity constant,
- {$\vec{u} = \vec{u}(\vec{x}, t)$} - the vector field of velocities,
- {$p = p(\vec{x},t)$} - the scalar-valued pressure field (implicit in the solver),
- {$\vec{F}_{\text{external}}$} - vector field of external forces such as gravity, bouyant forces or other forces.
Because solving for viscosity is computationally expensive and it can be cheaply imitated by applying a time dependent decay on the dye color, this project will simulate a 0 viscosity fluid and apply exponential decay on the dye colours.
The zero-divergence requirement {$\nabla \cdot \vec{u} = 0$} means that in the final solution of each iteration step the sum of fluid components entering each discrete cell will be equal to the components exiting the same cell. This kind of requirement will produce the swirliness of incompressible fluids in the simulation
Discrete Representation of Fluids
Let's say we want to discretely evaluate (solve) in a rectangular spatial domain {$R$} and time {$t>0$} some type of function {$f(x,y,t)$} which is governed by some partial differential equation. Let {$R$} be defined as {$$\begin{align} R = \{(x,y) \in \mathbb{R}| \ 0<x<a, 0<y<a\}, \quad a = \text{const.} >0. \end{align}$$} In finite-difference methods for solving PDE-s in 2D, the spatial domain is usually discretized by some fixed smallest spatial step for all coordinates as {$\Delta x = h_{x} >0, \Delta y = h_{y} >0$}. Time is also discretized to consist of smallest possible timesteps {$\Delta t >0$}. To assure the integer number of smallest spaces in a domain {$R$}, the smallest steps are defined by the integer numbers of cells in each direction. Let that be an integer constant {$N>0$} for both axes. Then the following can be stated about the spatial discretization:
{$$\begin{align} &\Delta x = \Delta y = \frac{a}{N} \end{align}$$} and the position vector {$\vec{x} =(x,y)$} coordinates in the domain can be expressed for some integer {$0\leq i\leq N$} as {$$\begin{align} &x_{i} = i\Delta x, \\ &y_{i} = i\Delta y, \end{align}$$} therefore there are {$N+1$} nodes from {$n=0$} to {$n=N$} where a scalar field {$f(x,y) =f(x_{i}, y_{j})$} can be evaluated therefore all the {$(N+1)(N+1)$} values of {$f(x_{i},y_{i})$} can be represented as a matrix {$$\begin{align} f_{ij}=f(x_{j},y_{i}) = \begin{pmatrix} f(x_{0},y_{0}) & f(x_{1},y_{0}) & \dots&f(x_{n},y_{0}) \\ f(x_{0},y_{1}) & f(x_{1},y_{1}) & \dots & f(x_{n},y_{0}) \\ \vdots&\vdots&\ddots&\vdots \\ f(x_{0},y_{n}) & f(x_{1},y_{n}) & \dots & f(x_{n},y_{n}) \end{pmatrix}. \end{align}$$}
Staggered Discretization of Velocities
In computational fluid dynamics, the pressure is (implicitly) discretized. However, to ensure the numerical ability to create 0-divergence, a staggered grid will be used to map velocity field components on the edges of cells that surround some pressure values. Therefore, if pressure for example takes on a matrix {$(p)_{i,j}$} with dimensions {$(N_{y}) \times (N_{x})$} in the form {$$(p)_{i,j} = \begin{pmatrix} p_{00}&p_{01}&\dots&p_{0,N_{x}-2}&p_{0,N_{x}-1} \\ p_{10}&p_{11}&\dots&p_{1,N_{x}-2}&p_{1,N_{x}-1} \\ \dots&\dots&\ddots&\vdots \\ p_{N_{y},0}&p_{N_{y},1}&\dots&p_{N_{y}-2,N_{x}-2}&p_{N_{y}-2,N_{x}-1} \\ p_{N_{y}-1,0}&p_{N_{y}-1,1}&\dots&p_{N_{y}-1,N_{x}-2}&p_{N_{y}-1,N_{x}-1} \end{pmatrix} \in \mathbb{R}^{N_{y}\times N_{x}},$$} then the corresponding velocity fields {$(u_{x})_{k,l}$} and {$(u_{y})_{m,n}$} will be created such that each pressure value {$p_{ij} = (p)_{i,j}$} will have two vertical neighbors that correspond to the vertical velocities {$(u_{y})_{i,j}$} its cell and {$(u_{y})_{i+1,j}$} below, and two horizontal neighbors that correspond to the horizontal velocities {$(u_{x})_{i,j}$} from the left and {$(u_{x})_{i,j+1}$} from the right.
Therefore {$$\begin{align} &(u_{x})_{k,l} \in \mathbb{R}^{N_{y}\times (N_{x}+1)}, \\ &(u_{y})_{m,n} \in \mathbb{R}^{(N_{y}+1)\times N_{x}}. \end{align}$$}
Computation Implementation
Since the update loop of the simulation includes computations up to time complexity {$\mathcal{O}(kN^2)$}, where {$N$} is the dimension of the grid and {$k \approx 10^1...10^2$} the number of Poisson iterations, the implementation will hold discrete fields info in data textures. The data textures will be passed into compute shaders each frame, where fast parallelized computation will take place. HLSL compute shaders are used alongside with C# scripts for the fluid interactive simulations in Unity game engine.
Update Loop
The simulation will start with initial conditions {$\vec{u}(\vec{x},0) = \vec{u}_{0}$}. Let the solution of the algorithm after some time {$t$} be {$\vec{u}(\vec{x}, t) = \vec{w}_{0}(\vec{x})$}.Then the update loop of the Stable Fluids simulation algorithm proceeds as follows:
{$$\vec{w}_{0}(\vec{x}) \overbrace{ \to }^{ \text{add force} } \vec{w}_{1}(\vec{x}) \overbrace{ \to }^{ \text{advect} } \vec{w}_{2}(\vec{x})\overbrace{ \to }^{ \text{remove divergence}}\vec{w}_{3}(\vec{x}).$$}