Realtime Fluid Simulation in Java

The Basics

The most common approach to simulating fluids is to use SPH (Smoothed Particle Hydrodynamics) This approach uses lots of particles to represent the fluid and they interact with each other by exerting a range of different forces on each other, all of which can be derived from the original navier-stokes equations. The forces are pressure, viscosity, external (boundaries, gravity or use interactions) and surface tension.

In every new frame these forces are calculated for each particle based on formulas and a smoothing kernel, which is a function that scales the force according to how close the particles are, so the forces stay local. The forces are then summed and divided by the mass to calculate acceleration (F=ma) which is integrated to calculate the velocity and then new position.

Each particle is defined to have a set mass, position, velocity. As well as universal constants for the entire fluid such as viscosity, pressure stiffness, smoothing width.

 

Implementation

The implementation I used worked in java libgdx, and simulated a 2d fluid of any parameters, the defaults are chosen to mimic water.

The algorithm used is as follows:

For each particle:

Compute the list of neighbours (in range of the smoothing width) to that particle (1)

For each particle:

Calculate the density of the particle (2)

Calculate the pressure of the particle (3)

For each particle:

Calculate the forces of the particle (4)

For each particle:

Calculate the acceleration of the particle and integrate (5)

(1) Neighbour Search

Illustration:

Capture

The code divides the world into a grid of width h (smoothing width of the kernel) and adds the particles in each grid block to a separate list and then adds each list into a hash map, therefore linking the list and the coordinates of that section together.

particlesGrid.clear();
 for(int i = 0 ; i < particles.size();i++)
 {
     p = particles.get(i);
      x = (int)(Math.floor((double)(p.pos.x)/gridwidth));
      y = (int)(Math.floor((double)(p.pos.y)/gridwidth));  //Generating the grid coordinates
     if(particlesGrid.containsKey(x+","+y))
     {
         particlesGrid.get(x+","+y).add(p);  //Adding to the list
     }else
     {
         a = new ArrayList();
         a.add(p);
         particlesGrid.put(x + "," + y, a); //Creating a new section in the grid
     }
     p.particleneighbour.clear();
 }

Then, starting from the most bottom left section of the grid it starts calculating the neighbours, it iterates

double checkingDistance = hwidth*hwidth;
 for(int i = 0 ; i < particles.size();i++)
 {
     p = particles.get(i);
     getNeighbours(p, checkingDistance); //Generate the list of neighbours for each particle
     p.density = 0;
 }
private void getNeighbours(Particle p, double checkingDistance) {
    int x = (int)(Math.floor((double)(p.pos.x)/ gridwidth));
    int y = (int)(Math.floor((double)(p.pos.y)/ gridwidth));

        parTemp.clear();
        parTemp.addAll(particlesGrid.get((x) + "," + (y)));
        Particle p2;
        for (int j = 0; j < parTemp.size(); j++) {
            p2 = parTemp.get(j);
            if ((p.num != p2.num) && (p.pos.dst2(p2.pos) <= checkingDistance)) {
                p.particleneighbour.add(p2);            //Add particles in its own grid
            }
        }

        parTemp.clear();
        if (particlesGrid.containsKey((x + 1) + "," + (y))) {   //Check the grid sections "ahead" of it
            parTemp.addAll(particlesGrid.get((x + 1) + "," + (y)));
        }

        if (particlesGrid.containsKey((x) + "," + (y + 1))) {
            parTemp.addAll(particlesGrid.get((x) + "," + (y + 1)));
        }

        if (particlesGrid.containsKey((x + 1) + "," + (y + 1))) {
            parTemp.addAll(particlesGrid.get((x + 1) + "," + (y + 1)));
        }
        if (particlesGrid.containsKey((x - 1) + "," + (y + 1))) {
            parTemp.addAll(particlesGrid.get((x - 1) + "," + (y + 1)));
        }

        for (int j = 0; j < parTemp.size(); j++) {
            p2 = parTemp.get(j);
            if ((p.pos.dst2(p2.pos) <= checkingDistance)) {
                p.particleneighbour.add(p2);
                p2.particleneighbour.add(p);  //Add them to each other lists (symmetry)
            }
        }
}

(2+3) Density and Pressure 

Looping through each particle it calculates the density of the particle and the effect of it on the other particles (symmetry again), as pressure is just a function of the density it can be readily calculated in this loop.

The formula for density is: (A list of kernels is provided at the end of the post)

Capture2The formula for pressure used is taits equation:pair_sph_tait

float kernel;
Iterator iter;
for(int i = 0 ; i < particles.size();i++)
    {
         p = particles.get(i);
        iter =p.particleneighbourden.iterator();
        while(iter.hasNext())
        {
            p2 = (Particle) iter.next();
            kernel = densityKernel(dis(p, p2), hwidth);
            p.density += p2.getMass()*kernel;
            p2.density += p.getMass()*kernel;
            p2.particleneighbourden.remove(p);
        }
       p.density += p.getMass()*densityKernel(dis(p,p),hwidth);
        p.pressure = (float)((kstiffness*pressurefactor*(1/restDensity))*(Math.pow(p.density/restDensity,pressurefactor) - 1));
    }
private float densityKernel(Vector2 r, float h) {
    float rmagsqr = r.len2();
    float hsqr = h*h;
    if(rmagsqr>= 0 && rmagsqr<= hsqr)
    {
        return (float)Math.pow(hsqr-rmagsqr,3)*calcDensityConstant;
//calcDensityConstant = (float)((315d)/(64d*Math.PI*Math.pow(hwidth,9)));
}else { return 0f; } }


(4) Forces

Looping through each particle it calculates the force of the particle (no symmetry this time). The forces are:

Force of pressure:  Capture3

Force of viscosity:Capture4

Force of surface tension (only added if the last condition holds):Capture9Capture6Capture7Capture8Capture10

Force of gravity (although added in last step 5):Capture5

Vector2 pressureForce ;
Particle p2;
Vector2 viscoForce ;
Vector2 surfaceForcen;
for(int i = 0 ; i < particles.size();i++)
    {
       p = particles.get(i);
       pressureForce = new Vector2(0,0);
       viscoForce = new Vector2(0,0);
       surfaceForcen = new Vector2(0,0);
       float surfaceForcedn = 0;

        for(int j = 0 ; j < p.particleneighbour.size();j++)
        {
            p2 = p.particleneighbour.get(j);

            pressureForce.sub(pressureKernel(dis(p, p2), hwidth).scl((p2.getMass() * (p.pressure + p2.pressure) * (1f / (2f * p2.density)))));

            viscoForce.add((p2.getVel().sub(p.vel)).scl(p2.getMass() * (1 / p2.density) * viscoKernel(dis(p, p2), hwidth)));

            surfaceForcen.add(dpoly6(dis(p, p2), hwidth).scl(p2.getMass() / p2.density));
            surfaceForcedn += (ddpol6(dis(p, p2), hwidth))*(p2.getMass() / p2.density);
        }

        float surfacel = surfaceForcen.len2();
        if(surfacel >= thresholdsqr)
        {

            p.Fsurface = surfaceForcen.scl(surfaceTension*(-surfaceForcedn/(float)Math.sqrt(surfacel)));
            p.redness = 0.5f;
        }else
        {
            p.redness = 0;
            p.Fsurface = Vector2.Zero;
        }
        p.Fsurface = surfaceForcen;
        viscoForce.scl(viscosity);
        p.Fviscosity = viscoForce;
        p.Fpressure = pressureForce;
    }

(5) Intergration

Looping through each particle it sums all the forces, divides by density to find the acceleration, then uses the leap frog time integration scheme to calculate the new velocity and position (leap frog has built-in average damping which is important for a smooth sph simulation). It then uses a simple boundary model that reflects the dampened velocity and position if the particle is at the boundary of the container.

Force summation:Capture11

Leap frog integration:

Capture12Capture13

Capture14

Capture15

for(int i = 0 ; i < particles.size();i++)
    {
        p = particles.get(i);

        if(!p.solid) {

            Vector2 totalForce = new Vector2(0,0);


           totalForce.add(p.Fpressure);
            totalForce.add(p.Fviscosity);
            totalForce.add(p.Fsurface);
            totalForce.scl(1 / p.density);
            totalForce.add(gravity);



            p.aftervel = p.beforevel.cpy().add(totalForce.scl(deltaTime));
            p.vel = (p.beforevel.cpy().add(p.aftervel)).scl(0.5f);
            p.beforevel = p.aftervel;
            p.pos.add(p.aftervel.cpy().scl(deltaTime));

            if(p.pos.y<=0)
                {

                    if(p.aftervel.y<0) {
                        p.aftervel.scl(cr, -cr);
                    }
                    p.vel = p.aftervel;
                     p.beforevel = p.aftervel;
                    p.pos.set(p.pos.x,-p.pos.y);
                }
            if(p.pos.x<=0)
                {

                    if(p.aftervel.x<0) {
                        p.aftervel.scl(-cr, cr);
                    }
                    p.vel = p.aftervel;
                    p.beforevel = p.aftervel;
                    p.pos.set(-p.pos.x, p.pos.y);
                }
            if(p.pos.x>=2)
                {

                    if(p.aftervel.x>0) {
                        p.aftervel.scl(-cr, cr);
                    }
                    p.vel = p.aftervel;
                    p.beforevel = p.aftervel;
                    p.pos.set(2-(p.pos.x-2),p.pos.y);
                }
        }
    }

Optimization

A large part of the optimization is due to the fact that the forces between particles are symmetrical, so the force exerted on pairs of particles by each other is the same. Therefore, the forces do not always have to be calculated as they may have already been calculated beforehand.

Computing the neighbours of a given particle can be sped up by using a grid based search. So the world is split into a grid of block (of width h) and then only the neighbouring blocks and the blocks the particle in need to be checked. Furthermore, because the neighbours are also symmetrical (if particle a is particle b’s neighbour then particle b is also particle a’s neighbour), therefore only the blocks “forward” of the particle need to be checked, so any previous block that has already been iterated over isn’t then re-iterated over.

Results

The actual code to simulate the fluid is not long at all, I used libgdx to help render as well as using the built in Vector2 class for vector addition and scaling.

It is able to simulate ~500 particles at 60fps and ~1000 particles at 30fp, with most of the cpu time being used for force calculation, specifically the surface tension.

The parameters I used for defaults are here:

timestep = 0.003f;
viscosity = 1.4f;
gravity = new Vector2(0f,-9.82f);

pressurefactor = 7;
restDensity = 80f;
mass = 0.01f;
cr = 0.95f;

kstiffness = 200f;
thresholdsqr = 55f;
simTime = 0;
realTime = 0;
hwidth = 0.1f;
gridwidth = hwidth;
surfaceTension = 0.05f;

calcDensityConstant = (float)((315d)/(64d*Math.PI*Math.pow(hwidth,9)));
calcPressureConstant = (float)((-45d)/(Math.PI*Math.pow(hwidth,6)));
calcViscConstant = (float)((45d)/(Math.PI*Math.pow(hwidth,6)));
calcSurfConstant = (float)((-945d)/(32d*Math.PI*Math.pow(hwidth,9)));

The kernels I used:CaptureCCaptureBCaptureA

In code form:

private Vector2 dpoly6(Vector2 r, float h) {
    float rmagsqr = r.len2();
    float hsqr = h*h;
    if(rmagsqr> 0 && rmagsqr<= hsqr)
    {
        return (r.scl((float)(calcSurfConstant*Math.pow(hsqr-rmagsqr,2))));
    }else
    {
        return Vector2.Zero;
    }
}

private float ddpol6(Vector2 r, float h) {

    float rmagsqr = r.len2();
    float hsqr = h*h;


    if(rmagsqr> 0 && rmagsqr<= hsqr)
    {
        return (calcSurfConstant*(hsqr-rmagsqr)*(3*hsqr-7*rmagsqr));
    }else
    {
        return 0;
    }
}

private float viscoKernel(Vector2 r, float h) {
    float rmag = r.len();
    if(rmag> 0 && rmag<= h)
    {
        return (h-rmag)*calcViscConstant;
    }else
    {
        return 0f;
    }
}

private float densityKernel(Vector2 r, float h) {
    float rmagsqr = r.len2();
    float hsqr = h*h;
    if(rmagsqr>= 0 && rmagsqr<= hsqr)
    {
        return (float)Math.pow(hsqr-rmagsqr,3)*calcDensityConstant;
    }else
    {
        return 0f;
    }
}

private Vector2 pressureKernel(Vector2 r, float h) {
    float rmag = r.len();
    if(rmag> 0 && rmag<= h)
    {
        return (r.scl((float)(calcPressureConstant*Math.pow(h-rmag,2)*(1f/rmag))));
    }else
    {
        return Vector2.Zero;
    }
}
//hello Roweena!! ;)

Some images and gifs of the simulation for rendering I simply rendered each particle as a blue circle, if they are a surface particle they are drawn more purple.

Dam Break:

The velocity of the particles:

6

The internal forces (pressure, viscosity and surface tension):

Splash test:

 

Splash test of viscous fluid:

 

Droplets in zero gravity coalesce:


Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s