Simulating Fluids

Made by: Joseph Oliai

There are various ways for simulating fluids and they all have different results, the question is how can we make a fluid that looks and behaves like the real deal? By using physics and math we could find an answer to this question.

Table of contents

• 1. Introduction
• 2. Inspiration
• 3. How do fluid simulations work?
  • 3.1 Navier Stokes Equation
  • 3.2 Eulerian Grids
  • 3.3 Lagrangian Particles
• 4. Smoothed Particle Hydrodynamics
  • 4.1 Spatial Hash Grids
  • 4.2 Applying the forces on the particles
  • 4.3 What is a Kernel Function
  • 4.4 Density
  • 4.5 Pressure
  • 4.6 Viscosity
  • 4.7 External Forces
• 5. Rendering
  • 5.1 Marching Squares
  • 5.2 Linear Interpolated Marching Squares
• 6. Interaction
  • 6.1 Collision Check
  • 6.2 User Interaction
• 7. Performance
  • 7.1 Rendering
  • 7.2 Calculation changes
• 8. Results meets expectations
• 9. Conclusion
• 10. Future Reference
• 11. Sources

1. Introduction

This research will be about fluid simulations. The reason this topic has been chosen is because I wanted to research more about fluids and how I can achieve making them look more realistic. In this research the question “How can an interactive fluid be recreated using Smoothed Particle Hydrodynamics?” will be answered.

We will take a look at something called Smoothed Particle Hydrodynamics (SPH) which is a method for simulating fluids. The results of implementing this method is often a fluid that looks and behaves like the real deal.

2. Inspiration

The biggest inspiration as to why this topic has been chosen is because of a YouTube video by Two Minute Papers that showed an example of just how beautiful fluid simulations can be.

3. How do fluid simulations work

Understanding how fluid simulations exactly work in a mathematical way is really complicated, there is a lot of high level math involved. However there are a lot of resources online explaining the basics of fluid simulations. First of all there are different methods for simulating fluids.

3.1 Navier Stokes Equation

The mathematical equation fluid simulations are based on is an equation called the Navier Stokes Equation. The equation is based on newton’s second law and the conservation of mass (What Are the Navier-Stokes Equations?, n.d.). The Navier Stokes Equation’s formula looks like this.

It exists out of a few components:
ρ is the density.
u is the velocity.
p is the pressure.
t is the time.
τ is viscosity.
g are the external forces.

In the case of fluid simulations we usually want to see mass as density (Numberphile, 2019). As one maybe can already tell the left side of the equation basically states newton’s second law namely mass * acceleration. The right side of the equation is the calculation of forces consisting of internal forces and external forces. The internal forces consist out of three forces; density, pressure and viscosity. The external forces are the forces outside applying upon the fluid, usually just gravity. These forces are very important for the recreation of fluids and in this project they are being used to make Smoothed Particle Hydrodynamics.

3.2 Eulerian Grids

One method of simulating fluids is using Eulerian Grids. An image is provided below to visualize what it looks like.

The way it works is that each square on the grid has a density and diffusion factor (West, n.d.). Based on the density of the square the intensity of the color changes. The higher the density of the square the more alpha the square has. The direction of a square is based on the diffusion and density of the surrounding cells.

Making use of Eulerian grids is one of the more traditional ways of recreating fluid simulations as it was the first method of fluid simulations. However in this research a different approach for simulating fluids is used the reason being that making use of Eulerian Grids usually differs in visual results which was not desirable.

3.3 Lagrangian Particles

In this research Lagrangian Particles are used to recreate fluids. Lagrangian Particles is an approach to creating Fluid simulations, small units with mass are tracked using particles (Månsson, 2013). This approach of tracking the particles is called Smoothed Particle Hydrodynamics. The particles all contain certain properties such as velocity and position. The particles will behave in a certain way based on viscosity, density and external forces such as gravity. 

4. Smoothed Particle Hydrodynamics

4.1 Spatial Hash Grids

The way Smoothed Particle Hydrodynamics work is to check the distance between the particles (Månsson, 2013). A few functions called kernel functions will then return values based on the particles attributes such as density and pressure which results in a certain behavior. Doing so is really intensive since if the user uses a lot of particles every particle is going to check if every other particle is in range. For this problem there is a solution called a Spatial Hash Grid. For this project a video (SimonDev, 2020) has been used explaining how a Spatial Hash Grid works and can be programmed.

How a Spatial Hash Grid works is that the worldspace is divided into several cells, like so:

Every cell gets an ID and a list of particles that are within the cell, this gets updated every update. The cells are given a Hash Value based on the x position on the grid + y multiplied by the amounts of columns. After doing so we will get a 1 dimensional grid which looks like this:

This grid has a potential to reduce the amount of calculations by a significant amount. Instead of checking whether every particle is in range of every other particle the program can now just check the neighboring cells and not calculate the rest, an example is shown in the figure above.

4.2 Applying the forces on the particles

The forces that are being used are based on the Navier Stokes Equation, density, pressure, viscosity and external forces. The way these forces are partially calculated is making use of something called kernel functions. The calculations and formulas for calculating the forces on the particles were found in a video (Xiangyu Hu, 2021) which is about Smoothed Particle Hydrodynamics.

4.3 What is a Kernel Function

A kernel function is generally used for Machine Learning purposes. However a kernel function can be used for everything that has a x and y component in this case particle positions, but what exactly is a kernel function? In a nutshell a kernel function is a method for faster calculations which otherwise results in a much more difficult calculation (Afonja, 2018).

The Kernel functions I used we’re found in a research paper (Ertekin, 2021) about Smoothed Particle Hydrodynamics. The way the kernel functions mathematically work exactly and why certain values have a specific value was really hard to find and in this research paper there will be no explanation about how they mathematically work.

But how are kernel functions applied to SPH? Let’s say there is a particle that needs to calculate it’s forces. The way the forces are calculated is to check for surrounding particles that are within a certain range, this is called neighborhood. The closer the surrounding particle is to the target particle the higher the weight value of the force. After this a summation of all surrounding particles weight are summarized and an average weight is calculated. The question is how does one decide what the weight value should be multiplied with? For this we use kernel functions.

4.4 Density

When thinking of density we can think of it as a measurement for mass per unit. Density is often used to express grams per cubic centimeter. In other words how much mass is in a certain volume?

For the density a kernel function called a Poly6 (sixth degree polynomial) kernel function is used. The way the density is calculated is with the following formula

The density of particle i is equal to the mass of particle j multiplied by the kernel value based on distance. This value called the density value is then added to the density of particle i. The image below shows the way this formula has been implemented in code.

To show what density may look like a setup is made where every particle is drawn as a ball. Each ball has a color, green being that the particle is really dense and blue being that the particle has a really low density. As one can see in the video below particles that are lower in the fluid have a really high density. The fact that the particles distances at the bottom of the fluid are relatively short in comparison to the top part of the fluid also is an indication that the particles at the bottom are more dense than the top part of the fluid.

There is one issue when using the Poly6 kernel namely that when two particles get really near each other the particles have little to no repulsive force meaning that particles can cluster which is not ideal. To solve this problem there is another kernel function called Spiky Kernel which is used for the calculation of pressure.

4.5 Pressure

Pressure is a force that is applied to an object’s surface area which is then distributed throughout the whole surface area. But how is this applied in this project? When particles get near each other pressure occurs. The amount of pressure is dependent on how close the particles are from each other and what pressure and density the surrounding particles have. For a part of the calculation of pressure a kernel function called Spiky Kernel is used. The following formula is to calculate the pressure.

The way the formula has been implemented in code looks like the following:

In this formula it is important that mass is conserved as this is one of the requirements of the Navier Stokes Equation. This basically means that every particle has a velocity and a reversed velocity and if these are substracted the answer should be 0. This is why there is a minus included in front of the equation. The reason we use a Spiky Gradient Kernel is because when only the Poly6 kernel is used to calculate the forces the particles can cluster. With the Spiky Gradient we resolve this problem and particles have repulsive forces.

To demonstrate the pressure on each particle a setup is made where every particle is drawn as a ball. Each ball has a color, red being that the particle is under a lot of pressure and blue being that the particle is not or little under pressure.

As one can see the particles on the bottom of the fluid are under high pressure whereas the particles on top are not under a lot of pressure at all. This is because as a particle gets lower in a fluid the more mass is stacked on top of the particle resulting in more pressure.

4.6 Viscosity

What is viscosity? Simply said viscosity is a measurement of friction between layers of fluid particles. When two layers of particles slide past each other there is friction, the higher amount of viscosity the more friction occurs thus resulting in a more viscous fluid.

The problem with the Spiky Kernel is that particles that are near each other can create a negative velocity field which causes particles to increase their relative velocity, this causes particles to drive to much apart. For this we use a Laplacian kernel function which ensures that this doesn’t happen. The solution that is used to simulate viscosity is the following:

An image below is provided to show the translated formula in code.

4.7 External Forces

The external forces in this project are the forces that occur besides the calculation of the viscosity, pressure and density such as gravity, walls and user interaction.

5. Rendering

5.1 Marching Squares

To render the particles we need to have an algorithm that can create meshes with various shapes. In this project a YouTube video (Sebastian Lague, 2015) and an article (Wong, 2014) has been used to implement the algorithm. The algorithm used in this project is called Marching Squares. So how does it work?

Let’s suppose the world space is divided in squares. Each square exists out of the corners which we will call control nodes and in between these nodes you have center points, these we will call center nodes.

Every control node holds a binary value. When the corner is active the value is added to the total binary value of the square. In total there are 16 (2^4) possibilities of outcome. When the binary value of a square is assigned a mesh can be generated based on that value with a triangulation table which is shown in the image below.

Suppose a binary value of 7 is used, the square that will be generated will be based on the triangulation table so in the case of a binary value of 7 a line between the left and top center node is drawn.

So how is this applied to SPH? The first step is to check for every square if a particle is in a certain range. If a particle happens to be in range we will make the corner value that is within range active, by doing so a binary value is received and can be fed to the triangulation table which then can be translated to vertices and triangles. To demonstrate how this looks a gif is provided.

In this case the rendered mesh look really blocky, one way to make it more pretty is to reduce the square size which results in more squares drawn in the world and thus creating a higher resolution.

The problem with this solution is that it requires a lot of processing power. Suppose the world is 10 units wide and high and we have a square size of 0.1, this means we have a grid of 100 by 100 squares. And now suppose we have 100 particles. This means we now have to check 100 * 100 * 100 = 1000000 times every update which is really heavy. One way to increase smoothness of the particles other than upping the resolution is called Linear Interpolated Marching Squares.

5.2 Linear Interpolated Marching Squares

Linear Interpolated Marching Squares has the same concept as regular Marching Squares there is just a slight difference. In regular Marching Squares the Control Nodes are either on or off (1 or 0), and lines are being drawn between Center Nodes based on whether the control nodes are on or off. In the Linear Interpolated version, the Control Nodes hold a float value between 0 and 1. With this in mind we can then lerp between those values to get different results other than drawing a line between center nodes. An example is shown in the image below. Because of this the resolution does not has to be as high and the result will be a much smoother mesh which saves processing power.

In this project regular marching squares has been used which is not ideal but because of difficulty and time the linear interpolated marching squares algorithm has not been implemented. The solution for saving processing power was in this case rather easy and the solution was to shorten the world space and lowering the marching squares resolution.

6. Interaction

6.1 Collision Check

Because the particles used in this projects are no prefabs and don’t inherit from monobehaviour a way of collision checking has to be made. In the case of this project the only collision check that has been made is a collision check between cubes and circles (particles).

First of all an object manager is made that stores the values of the cubes being instantiated. When this is done the object manager can communicate the transforms of the objects in the list to the fluid manager. The fluid manager can then check for collisions between the moving particles and objects in the object manager. The procedure is visualized in the image below.

With the help of an article (Day, 2018) that explained how to do collision checks between circles and rectangles I managed to make a method that detects collisions. The first thing was to calculate the distance between the center of the rectangle and the center of the circle. In the image below we can see that everything outside of RH+CR or CR+RH is out of bounds.

After this step we can check if the circles are inside of the rectangle. We can say that if the distance between the center of the circle and the center of the rectangle is larger than RW the circle is colliding, but as one can see in the image below when a circle is in a corner it will not be detected.

To fix this there is a solution. First the Euclidian distance between the center of the circle and the center of the rectangle is calculated. After that we check if the Euclidian distance is equal or larger than the squared radius of the circle. After implementing this we now have a fully functional collision detection system between circles and rectangles which can be used to apply forces to the particles when colliding. This force is equal to the negative normal vector of the particle multiplied by a velocity damping.

6.2 User Interaction

What is nice about the scripts used in this project is that everything is really easy to setup. The user only has to add a prefab to the scene. After that the user can tweak the values in the inspector such as mass, viscosity and how many particles are spawned when right clicking.

The user can left click to thresh the particles around and the user can use “C” to spawn cubes into the scene. The size and spawn cooldown of these cubes can be tweaked in the object manager component.

7 Performance

7.1 Rendering

So all these calculations of distances and forces can be really heavy and we want a way to enhance the performance. Let’s start with looking what the performance is at the moment and which classes and methods take up a lot of memory. I noticed that world space makes a lot of difference even though the amounts of particles are the same? So how can that be? Even though only 6 particles are used in a space of 10 by 10, the FPS drops to 4 with the deep profiler on (The deep profiler in unity can specifically see which methods take in a lot of memory, though this is nice it drops the performance of the simulation significantly.)

In the image above one can see which methods take a lot of memory. Subtraction takes in 16% of the memory, this is a part in my code where the squared distance is substracted from a particle’s position. For the mesh generator a lot of distance calculation is done to decide which marching squares should be active. So let’s see what happens when the mesh generator code is commented out. Will it help the performance?

Nice! The fps went from 4 particles at 4 fps to 100 particles with 4 fps (With the deep profiler on) and the weird math calculations are out of the list in the deep profiler. That’s a huge difference. So how can we improve the performance? Sadly the marching squares algorithm is just not competent enough to render in real time. How the code works at the moment is that each particle is distance checking with each grid element in the marching squares algorithm, and ideally we want to use something like the hash grid system we made for the particles forces calculations, or a shader that can render these particles.

To achieve better results we have to scrap the marching squares algorithm and use a more simple solution. There is no real name for this solution so I named it alpha cutout rendering and the inspiration for this idea comes from a research paper (Månsson, 2013). The main idea behind this type of rendering is that each particle becomes a somewhat large circle. After we transformed the particles in circles the circles get a shader that makes them blurry where the edges of the circles have an opacity of 0 and the center an opacity of 1. When this step has been applied we can assign a threshold that says e.g. that everything above an alpha of 0.3 should be rendered sharply and the rest should be cut out of the image. Particles that are near each other have overlapping edges which causes the edges to have a bit higher alpha which creates the illusion that the particles are one homogamous blob. For the rendering method I used a preexisting project (Vue code, 2018) which does exactly what I described.

So when we apply this method how does it look like in comparison to the marching squares algorithm? And how does it perform?

I can proudly say that first of all it looks better shapewise and second of all the performance is much better! The FPS is consistent and the world space can be much bigger since the world is not divided in squares anymore and the particles don’t have to check their distances with every square.

Let’s compare our data. 60 particles with the new rendering method has a whopping 400 FPS! (without deep profiler). Whereas the old marching squares algorithm has 15 fps with 60 particles (without deep profiler), that’s an improvement of almost 27 times. Problem solved right?

7.2 Calculation changes

When observing the fluid I noticed that the fluid was not behaving correctly. The viscosity factor did not contribute a whole lot and overall the movement of the fluid was a bit off. So I dived into the code and found a few things. First of all there was a equation error in the code. The viscosity Gradient was never added to the viscosity force which means that the whole time viscosity was calculated incorrectly.

The second thing I did was looking up more research papers about SPH and I found that the Kernel functions are different for 3D and 2D (VIJAYKUMAR, 2012). During this project the kernel functions that had been used are meant for 3D, which meant I had to change the functions to a 2D version. So what did my fluid look like before and how did it look after? For the comparison I tested two scenarios. Throwing the particles in the air and adding particles.

Before results; So as you can see the fluid is really bouncy and it seems really tense. This is due to viscosity not playing a role in the fluid.

After results;

Not only does it look better but the performance is even better. The simulation can handle more particles. This is due to the fact that before the new calculations were implemented a lot of particles accumulated on the bottom. This resulted into having to distance check for more particles on the bottom grids. For every particle that is added in a grid the amount of calculations that have to be done grows exponentially (Reminder a particle can only distance check with particles that are within the same grid). Further more the particles are making little groups of globs which it did way less before the new implementation.

8 Results meets expecations

So what is the end result of this project. Well the initial idea was to make a 3D version of fluid simulations which performs well, but after some research it seemed to be one of the most difficult topics I have ever learned about and I changed my scope to make a simulation in 2D. The first iteration with marching squares looked cool but performed really badly and honestly I really wanted to improve on that so I did.
The second iteration of the fluid actually resulted in something more stable and prettier in my opinion. Not only did it look like the countless videos I’ve watched on the internet, but it also gave me a lot of knowledge about several topics like Spatial Hash Grids which can be used for games with a lot of entities, or Marching Cubes/Squares which can be used for terrain generation. So it is save to say that I am really happy with the results and the journey I had in making a fluid.

9 Conclusion

The results of this project is a fluid that pretty much acts and behaves in the way a realistic fluid would do. This is due to the Navier stokes equation. The performance is not great when a user wants to use a lot of particles and there is definitely room for improvement in terms of calculations of the forces and calculations of rendering.

10 Future reference

If I were to work more on this project I definitely would improve neighbor searching since there is definitely room for improvement. I chose not to do that in my project reason being that I had to rewrite almost 80% of my code if I did. Also it would be nice if the fluid could interact with poly shapes instead of just squares. Furthermore it would be a very interesting topic to make this project into a 3D simulation rather than a 2D simulation.

11 Sources

Afonja, T. (2018, 13 juli). Kernel Functions – Towards Data Science. Medium. Geraadpleegd op 8 november 2021, van https://towardsdatascience.com/kernel-function-6f1d2be6091

Day, J. (2018, 23 juli). Collision Detection – Circles, Rectangles and Polygons. Game Development Blog. Geraadpleegd op 1 november 2021, van https://www.gamedevelopment.blog/collision-detection-circles-rectangles-and-polygons/

Ertekin, B. (2021, augustus). Fluid Simulation using Smoothed Particle Hydrodynamics. Burak Ertekin. https://nccastaff.bournemouth.ac.uk/jmacey/MastersProject/MSc15/06Burak/BurakErtekinMScThesis.pdf

khskarl. (2016, 26 december). sph-unity. Github. Geraadpleegd op 1 oktober 2021, van https://github.com/khskarl/sph-unity/tree/master/Assets

Månsson, D. (2013). Interactive 2D Particle-based Fluid Simulation for Mobile Devices. Daniel Månsson. https://www.diva-portal.org/smash/get/diva2:676516/FULLTEXT01.pdf

Numberphile. (2019, 27 augustus). Navier-Stokes Equations – Numberphile [Video]. YouTube. https://www.youtube.com/watch?v=ERBVFcutl3M&t=343s&ab_channel

Sebastian Lague. (2015, 1 maart). [Unity] Procedural Cave Generation (E02. Marching Squares) [Video]. YouTube. https://www.youtube.com/watch?v=yOgIncKp0BE&t=1s&ab_channel=SebastianLague

SimonDev. (2020, 30 november). Spatial Hash Grids & Tales from Game Development [Video]. YouTube. https://www.youtube.com/watch?v=sx4IIQL0x7c&ab_channel=SimonDev

West, M. (z.d.). Practical Fluid Dynamics: Part 1. Gamasutra. Geraadpleegd op 28 september 2021, van https://www.gamasutra.com/view/feature/1549/practical_fluid_dynamics_part_1.php?print=1

What Are the Navier-Stokes Equations? (z.d.). Simscale. Geraadpleegd op 1 oktober 2021, van https://www.simscale.com/docs/simwiki/numerics-background/what-are-the-navier-stokes-equations/

Wong, J. (2014, 19 augustus). Metaballs and Marching Squares. Jamie Wong. Geraadpleegd op 1 november 2021, van http://jamie-wong.com/2014/08/19/metaballs-and-marching-squares/

Xiangyu Hu. (2021, 9 februari). Introduction on SPH method [Video]. YouTube. https://www.youtube.com/watch?v=2RoITtjpE-w&ab_channel=XiangyuHu

VIJAYKUMAR, A. (2012). Smoothed Particle Hydrodynamics Simulation for Continuous Casting. ADITHYA VIJAYKUMAR. https://www.diva-portal.org/smash/get/diva2:573583/FULLTEXT01.pdf

Vue Code. (2018, 31 juli). Unity Liquid Simulation Tutorial [Video]. YouTube. https://www.youtube.com/watch?v=3r6o8CfRlA4&ab_channel=VueCode