Why my fluids don’t flow

December 14, 2010

I have an unopened copy of Digital Color Management sitting on my desk. It’s staring at me accusingly.

In order to keep myself distracted from its dirty looks, I’ve been tinkering around with fluid simulation. Miles Macklin has done some great work with Eulerian (grid based) solvers, so in an effort to distance myself from the competition, I’m sticking to 2D Lagrangian (particle based) simulation.

Until recently, I’d always thought that particle based fluid simulation was complicated and involved heavy maths. This wasn’t helped by the fact that most of the papers on the subject have serious sounding names like Particle-based Viscoelastic Fluid Simulation, Weakly compressible SPH for free surface flows, or even Smoothed Particle Hydrodynamics and Magnetohydrodynamics.

It wasn’t until I finally took the plunge and tried writing my own Smoothed Particle Hydrodynamics simulation that I found that it can be quite easy, provided you work from the right papers. SPH has a couple of advantages over grid based methods: it is trivial to ensure that mass is exactly conserved, and free-surfaces (the boundary between fluid and non-fluid) come naturally. Unfortunately, SPH simulations have a tendency to explode if the time step is too large and getting satisfactory results is heavily dependent on finding “good” functions with which to model the inter-particle forces.

I had originally intended to write an introduction to SPH, but soon realised that it would make this post intolerably long, so instead I’ll refer to the papers that I used when writing my own sim. Pretty much every SPH paper comes with an introduction to the subject, invariably in section 2. Related Work.

The first paper I tried implementing was Particle-Based Fluid Simulation for Interactive Applications by Müller et. al. It serves as a great introduction to SPH with a very good discussion of kernel weighting functions, but I had real difficulty getting decent results. In the paper pressure, viscosity and surface tension forces are modeled using following equations:

\textbf{f}_i^{pressure} = -\sum_{j}m_j\frac{p_j}{\rho_j}\nabla{W(\textbf{r}_i-\textbf{r}_j, h)}

\textbf{f}_i^{viscosity} = \mu\sum_{j}m_j\frac{\textbf{v}_j-\textbf{v}_i}{\rho_j}\nabla^2W(\textbf{r}_i-\textbf{r}_j, h)

c_S(\textbf{r}) = \sum_jm_j\frac{1}{\rho_j}W(\textbf{r}-\textbf{r}_j, h)

\textbf{f}_i^{surface} = -\sigma\nabla^2c_S\frac{\textbf{n}}{|\textbf{n}|}

The pressure for each particle is calculated from its density using:

P_i = k(\rho_i - \rho_0) where \rho_0 is the some non-zero rest density.

The first problem I encountered was with the pressure model; it only acts as a repulsive force if the particle density is greater than the rest density. If a particle has only a small number of neighbours, the pressure force will attract them to form a cluster of particles all sharing the same space. In my experiments, I often found large numbers of clusters of three or four particles all in the same position. It took me a while to figure out what was going on because Müller states that the value of the rest density “mathematically has no effect on pressure forces”, which is only true given a fairly uniform density of particles far from the boundary.

The second problem I found was with the surface tension force. It was originally developed for multiphase fluid situations with no free surfaces and doesn’t behave well near the surface boundary; in fact it can actually pull the fluid into concave shapes. Additionally, because it’s based on a Laplacian, it’s very sensitive to fluctuations in the particle density, which are the norm at the surface boundary.

After a week or so of trying, this was my best result:

From the outset, you can see the surface tension force is doing weird things. Even worse, once the fluid starts to settle the particles tended to stack on top of each and form a very un-fluid blob.

On the up side, I did create possibly my best ever bug when implementing the surface tension model; I ended up with something resembling microscopic life floating around under the microscope:

The next paper I tried was Particle-based Viscoelastic Fluid Simulation by Clavet et al. I actually had a lot of success with their paper and had a working implementation of their basic model up and running in less than two hours. Albeit minus the viscoelasticity. In addition to the pressure force described in Müller’s paper, they model “near” density and pressure, which are similar to their regular counterparts but with a zero rest density and different kernel functions:

P_i = k(\rho_i - \rho_0)

P_i^{near} = k^{near} \rho_i^{near}

\rho_i = \sum_j (1 - \frac{|\textbf{r}_i - \textbf{r}_j|}{h}) ^ 2

\rho_i^{near} = \sum_j (1 - \frac{|\textbf{r}_i - \textbf{r}_j|}{h}) ^ 3

This near pressure ensures a minimum spacing and as an added bonus performs a decent job of modelling surface tension too. This is the first simulation I ran using their pressure and viscosity forces:

Although initial results were promising, I struggled when tweaking the parameters to find a good balance between a fluid that was too compressible and one that was too viscous. Also, what I really wanted was to do multiphase fluid simulation. This wasn’t covered in the viscoelastic paper, so my next port of call was Weakly compressible SPH for free surface flows by Becker et al. In this paper, surface tension is modeled as:

\textbf{f}_i^{surface} = -\frac{\kappa}{m_i} \sum_j m_j W(\textbf{r}_i-\textbf{r}_j) (\textbf{r}_i - \textbf{r}_j)

They also discuss using Tait’s equation for the pressure force, rather than one based on the ideal gas law:

P_i = B((\frac{\rho_i}{\rho_0})^\gamma - 1) with \gamma = 7

I gave that a shot, but the large exponent caused the simulation to explode unless I used a really small time step. Instead, I found that modifying the pressure forces from the viscoelastic paper slightly gave a much less compressible fluid without the requirement for a tiny time step:

\rho_i = \sum_j (1 - \frac{|\textbf{r}_i - \textbf{r}_j|}{h}) ^ 3

\rho_i^{near} = \sum_j (1 - \frac{|\textbf{r}_i - \textbf{r}_j|}{h}) ^ 4

Here’s one of my more successful runs:

And here is a slightly simplified version of the code behind it. Be warned, it’s quite messy; I’m rather enjoying hacking code together these days:

#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <memory.h>
#include <glut.h>


#define kScreenWidth 640
#define kScreenHeight 480
#define kViewWidth 10.0f
#define kViewHeight (kScreenHeight*kViewWidth/kScreenWidth)
#define kPi 3.1415926535f
#define kParticleCount 3000

#define kRestDensity 82.0f
#define kStiffness 0.08f
#define kNearStiffness 0.1f
#define kSurfaceTension 0.0004f
#define kLinearViscocity 0.5f
#define kQuadraticViscocity 1.0f

#define kParticleRadius 0.05f
#define kH (6*kParticleRadius)
#define kFrameRate 20
#define kSubSteps 7

#define kDt ((1.0f/kFrameRate) / kSubSteps)
#define kDt2 (kDt*kDt)
#define kNorm (20/(2*kPi*kH*kH))
#define kNearNorm (30/(2*kPi*kH*kH))


#define kEpsilon 0.0000001f
#define kEpsilon2 (kEpsilon*kEpsilon)


struct Particle
{
    float x;
    float y;

    float u;
    float v;

    float P;
    float nearP;

    float m;

    float density;
    float nearDensity;
    Particle* next;
};

struct Vector2
{
    Vector2() { }
    Vector2(float x, float y) : x(x) , y(y) { }
    float x;
    float y;
};

struct Wall
{
    Wall() { }
    Wall(float _nx, float _ny, float _c) : nx(_nx), ny(_ny), c(_c) { }
    float nx;
    float ny;
    float c;
};

struct Rgba
{
    Rgba() { }
    Rgba(float r, float g, float b, float a) : r(r), g(g), b(b), a(a) { }
    float r, g, b, a;
};

struct Material
{
    Material() { }
    Material(const Rgba& colour, float mass, float scale, float bias) : colour(colour) , mass(mass) , scale(scale) , bias(bias) { }
    Rgba colour;
    float mass;
    float scale;
    float bias;
};

#define kMaxNeighbourCount 64
struct Neighbours
{
    const Particle* particles[kMaxNeighbourCount];
    float r[kMaxNeighbourCount];
    size_t count;
};

size_t particleCount = 0;
Particle particles[kParticleCount];
Neighbours neighbours[kParticleCount];
Vector2 prevPos[kParticleCount];
Vector2 relaxedPos[kParticleCount];
Material particleMaterials[kParticleCount];
Rgba shadedParticleColours[kParticleCount];

#define kWallCount 4
Wall walls[kWallCount] =
{
    Wall( 1,  0, 0),
    Wall( 0,  1, 0),
    Wall(-1,  0, -kViewWidth),
    Wall( 0, -1, -kViewHeight)
};

#define kCellSize kH
const size_t kGridWidth = (size_t)(kViewWidth / kCellSize);
const size_t kGridHeight = (size_t)(kViewHeight / kCellSize);
const size_t kGridCellCount = kGridWidth * kGridHeight;
Particle* grid[kGridCellCount];
size_t gridCoords[kParticleCount*2];


struct Emitter
{
    Emitter(const Material& material, const Vector2& position, const Vector2& direction, float size, float speed, float delay)
        : material(material), position(position), direction(direction), size(size), speed(speed), delay(delay), count(0)
    {
        float len = sqrt(direction.x*direction.x + direction.y*direction.y);
        this->direction.x /= len;
        this->direction.y /= len;
    }
    Material material;
    Vector2 position;
    Vector2 direction;
    float size;
    float speed;
    float delay;
    size_t count;
};

#define kEmitterCount 2
Emitter emitters[kEmitterCount] =
{
    Emitter(
        Material(Rgba(0.6f, 0.7f, 0.9f, 1), 1.0f, 0.08f, 0.9f),
        Vector2(0.05f*kViewWidth, 0.8f*kViewHeight), Vector2(4, 1), 0.2f, 5, 0),
    Emitter(
        Material(Rgba(0.1f, 0.05f, 0.3f, 1), 1.4f, 0.075f, 1.5f),
        Vector2(0.05f*kViewWidth, 0.9f*kViewHeight), Vector2(4, 1), 0.2f, 5, 6),
};


float Random01() { return (float)rand() / (float)(RAND_MAX-1); }
float Random(float a, float b) { return a + (b-a)*Random01(); }


void UpdateGrid()
{
    // Clear grid
    memset(grid, 0, kGridCellCount*sizeof(Particle*));

    // Add particles to grid
    for (size_t i=0; i<particleCount; ++i)
    {
        Particle& pi = particles[i];
        int x = pi.x / kCellSize;
        int y = pi.y / kCellSize;

        if (x < 1)
            x = 1;
        else if (x > kGridWidth-2)
            x = kGridWidth-2;

        if (y < 1)
            y = 1;
        else if (y > kGridHeight-2)
            y = kGridHeight-2;

        pi.next = grid[x+y*kGridWidth];
        grid[x+y*kGridWidth] = &pi;

        gridCoords[i*2] = x;
        gridCoords[i*2+1] = y;
    }
}


void ApplyBodyForces()
{
    for (size_t i=0; i<particleCount; ++i)
    {
        Particle& pi = particles[i];
        pi.v -= 9.8f*kDt;
    }
}


void Advance()
{
    for (size_t i=0; i<particleCount; ++i)
    {
        Particle& pi = particles[i];

        // preserve current position
        prevPos[i].x = pi.x;
        prevPos[i].y = pi.y;

        pi.x += kDt * pi.u;
        pi.y += kDt * pi.v;
    }
}


void CalculatePressure()
{
    for (size_t i=0; i<particleCount; ++i)
    {
        Particle& pi = particles[i];
        size_t gi = gridCoords[i*2];
        size_t gj = gridCoords[i*2+1]*kGridWidth;

        neighbours[i].count = 0;

        float density = 0;
        float nearDensity = 0;
        for (size_t ni=gi-1; ni<=gi+1; ++ni)
        {
            for (size_t nj=gj-kGridWidth; nj<=gj+kGridWidth; nj+=kGridWidth)
            {
                for (Particle* ppj=grid[ni+nj]; NULL!=ppj; ppj=ppj->next)
                {
                    const Particle& pj = *ppj;

                    float dx = pj.x - pi.x;
                    float dy = pj.y - pi.y;
                    float r2 = dx*dx + dy*dy;
                    if (r2 < kEpsilon2 || r2 > kH*kH)
                        continue;

                    float r = sqrt(r2);
                    float a = 1 - r/kH;
                    density += pj.m * a*a*a * kNorm;
                    nearDensity += pj.m * a*a*a*a * kNearNorm;

                    if (neighbours[i].count < kMaxNeighbourCount)
                    {
                        neighbours[i].particles[neighbours[i].count] = &pj;
                        neighbours[i].r[neighbours[i].count] = r;
                        ++neighbours[i].count;
                    }
                }
            }
        }

        pi.density = density;
        pi.nearDensity = nearDensity;
        pi.P = kStiffness * (density - pi.m*kRestDensity);
        pi.nearP = kNearStiffness * nearDensity;
    }
}


void CalculateRelaxedPositions()
{
    for (size_t i=0; i<particleCount; ++i)
    {
        const Particle& pi = particles[i];

        float x = pi.x;
        float y = pi.y;

        for (size_t j=0; j<neighbours[i].count; ++j)
        {
            const Particle& pj = *neighbours[i].particles[j];
            float r = neighbours[i].r[j];
            float dx = pj.x - pi.x;
            float dy = pj.y - pi.y;

            float a = 1 - r/kH;

            float d = kDt2 * ((pi.nearP+pj.nearP)*a*a*a*kNearNorm + (pi.P+pj.P)*a*a*kNorm) / 2;

            // relax
            x -= d * dx / (r*pi.m);
            y -= d * dy / (r*pi.m);

            // surface tension
            if (pi.m == pj.m)
            {
                x += (kSurfaceTension/pi.m) * pj.m*a*a*kNorm * dx;
                y += (kSurfaceTension/pi.m) * pj.m*a*a*kNorm * dy;
            }

            // viscocity
            float du = pi.u - pj.u;
            float dv = pi.v - pj.v;
            float u = du*dx + dv*dy;
            if (u > 0)
            {
                u /= r;

                float a = 1 - r/kH;
                float I = 0.5f * kDt * a * (kLinearViscocity*u + kQuadraticViscocity*u*u);

                x -= I * dx * kDt;
                y -= I * dy * kDt;
            }

        }

        relaxedPos[i].x = x;
        relaxedPos[i].y = y;
    }
}


void MoveToRelaxedPositions()
{
    for (size_t i=0; i<particleCount; ++i)
    {
        Particle& pi = particles[i];
        pi.x = relaxedPos[i].x;
        pi.y = relaxedPos[i].y;
        pi.u = (pi.x - prevPos[i].x) / kDt;
        pi.v = (pi.y - prevPos[i].y) / kDt;
    }
}


void ResolveCollisions()
{
    for (size_t i=0; i<particleCount; ++i)
    {
        Particle& pi = particles[i];

        for (size_t j=0; j<kWallCount; ++j)
        {
            const Wall& wall = walls[j];
            float dis = wall.nx*pi.x + wall.ny*pi.y - wall.c;
            if (dis < kParticleRadius)
            {
                float d = pi.u*wall.nx + pi.v*wall.ny;
                if (dis < 0)
                    dis = 0;
                pi.u += (kParticleRadius - dis) * wall.nx / kDt;
                pi.v += (kParticleRadius - dis) * wall.ny / kDt;
            }
        }
    }
}


void Render()
{
    glClearColor(0.02f, 0.01f, 0.01f, 1);
    glClear(GL_COLOR_BUFFER_BIT);

    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    glOrtho(0, kViewWidth, 0, kViewHeight, 0, 1);

    glEnable(GL_POINT_SMOOTH);
    glEnable(GL_BLEND);
    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

    for (size_t i=0; i<particleCount; ++i)
    {
        const Particle& pi = particles[i];
        const Material& material = particleMaterials[i];

        Rgba& rgba = shadedParticleColours[i];
        rgba = material.colour;
        rgba.r *= material.bias + material.scale*pi.P;
        rgba.g *= material.bias + material.scale*pi.P;
        rgba.b *= material.bias + material.scale*pi.P;
    }

    glEnableClientState(GL_VERTEX_ARRAY);
    glEnableClientState(GL_COLOR_ARRAY);

    glPointSize(2.5f*kParticleRadius*kScreenWidth/kViewWidth);

    glColorPointer(4, GL_FLOAT, sizeof(Rgba), shadedParticleColours);
    glVertexPointer(2, GL_FLOAT, sizeof(Particle), particles);
    glDrawArrays(GL_POINTS, 0, particleCount);

    glDisableClientState(GL_COLOR_ARRAY);
    glDisableClientState(GL_VERTEX_ARRAY);

    glutSwapBuffers();
}


void EmitParticles()
{
    if (particleCount == kParticleCount)
        return;

    static int emitDelay = 0;
    if (++emitDelay < 3)
        return;

    for (size_t emitterIdx=0; emitterIdx<kEmitterCount; ++emitterIdx)
    {
        Emitter& emitter = emitters[emitterIdx];
        if (emitter.count >= kParticleCount/kEmitterCount)
            continue;

        emitter.delay -= kDt*emitDelay;
        if (emitter.delay > 0)
            continue;

        size_t steps = emitter.size / (2*kParticleRadius);

        for (size_t i=0; i<=steps && particleCount<kParticleCount; ++i)
        {
            Particle& pi = particles[particleCount];
            Material& material = particleMaterials[particleCount];
            ++particleCount;

            ++emitter.count;

            float ofs = (float)i / (float)steps - 0.5f;

            ofs *= emitter.size;
            pi.x = emitter.position.x - ofs*emitter.direction.y;
            pi.y = emitter.position.y + ofs*emitter.direction.x;
            pi.u = emitter.speed * emitter.direction.x*Random(0.9f, 1.1f);
            pi.v = emitter.speed * emitter.direction.y*Random(0.9f, 1.1f);
            pi.m = emitter.material.mass;

            material = emitter.material;
        }
    }

    emitDelay = 0;
}


void Update()
{
    for (size_t step=0; step<kSubSteps; ++step)
    {
        EmitParticles();

        ApplyBodyForces();
        Advance();
        UpdateGrid();
        CalculatePressure();
        CalculateRelaxedPositions();
        MoveToRelaxedPositions();
        UpdateGrid();
        ResolveCollisions();
    }

    glutPostRedisplay();
}


int main (int argc, char** argv)
{
    glutInitWindowSize(kScreenWidth, kScreenHeight);
    glutInit(&argc, argv);
    glutInitDisplayString("samples stencil>=3 rgb double depth");
    glutCreateWindow("SPH");
    glutDisplayFunc(Render);
    glutIdleFunc(Update);
    
    memset(particles, 0, kParticleCount*sizeof(Particle));
    UpdateGrid();

    glutMainLoop();
    
    return 0;
}

I’m pretty happy with the results, even if at three seconds per frame for the video above, my implementation isn’t exactly fast. Here are a few other videos from various stages of development:

About these ads

31 Responses to “Why my fluids don’t flow”

  1. Niall Says:

    Hey, I’m the guy that just left a comment on one of your videos.
    Excuse me if I go on a bit of a ramble, you’ve been experiencing the same sort of issues & discoveries I have recently!

    I fully implemented the viscoelastics paper a few weeks ago in flash – http://niallmeister.deviantart.com/art/Viscoelastic-Fluid-Sim-174243339 – albeit with a few small changes. I’ve since redone it quite a bit, switched to Leapfrog integration and completely changed the spring algorithm to suit (hey, that whole bit isn’t even slightly physically accurate anyway, as long as it’s stable right?).

    The particle clumping effect you’re describing is pretty well known, referred to as a tensile instability. There have been various efforts to reduce / eliminate the effect, a good paper on it is “SPH Without a Tensile Instability” by J.J. Monaghan. I agree though, the results you get from double density relaxation are far nicer looking and more stable in general.

    There’s a common, very annoying issue specifically with the “Müller-esque” formulation of SPH. In 2 fluids of a high density contrast a gap tends to form between them, it’s detailed and a fix is proposed in the paper Density Contrast SPH Interfaces http://www.ifi.uzh.ch/arvo/vmml/admin/upload/Solenthaler_sca08.pdf
    I’ve tried doing multiphase sims with double density relaxation, things like changing the pressure force to massi*massj*(pressi+pressj)*weight+(pressneari+pressnearj)*weight*weight and so on but the gap starts to appear there too (and it’s not very physically viable, I don’t particularly mind though :P). Your implementation is pretty clever, I’m just wondering how you arrived at the kNorm and kNearNorm variables?

    • Tom Madams Says:

      Hi Niall, it’s great to hear from people working on the same stuff. I’d definitely noticed both the problems you mentioned, but hadn’t come across either of those papers, to thanks for the links. I’m looking forward to checking them out.

      Congrats on the Flash simulations, I’m impressed that the performance is so good. You’ve motivated me to try and implement the full viscoelastic sim.

      The kNear and kNearNorm variables are the normalization constants for the two kernel functions. They’re the scaling values required to ensure that the kernels integrate to 1.

      If I remember correctly, the normalization constants are folded directly into the stiffness and rest pressure constants in the viscoelastic paper. I guess that’s in keeping with the tone of that paper (keeping things simple), but I felt like it was hiding important details about the SPH method in general, so reinserted them explicitly back into my code.

      When analytically integrating the kernel functions in 2D, you have to add an extra 2*pi*r term to account for the fact that the functions are specified in a circular coordinate system.

      Given the kernel function:
      (1-r/h)^3

      Integrating from 0 to h in circular coordinates becomes:
      Integrate[((1-r/h)^3) * r*2*pi, {r, 0, h}]
      Here’s the result courtesy of wolframalpha.com

      Which evaluates to:
      (2*pi*h^2)/20

      Dividing by this constant ensures the kernel is normalized.

      kNearNorm is similarly derived.

      Since most of the papers talk about SPH in 3D, they’re working in a spherical coordinate system, so the integration is different. That caught me out for the longest time.

      • Niall Says:

        Thanks for the speedy elaboration! I’m only 13 so the maths side of this stuff is still pretty new to me and I haven’t fully dived into integrals yet, I somewhat understand though.
        One thing I’m wondering about, it doesn’t seem like when integrating the equations of motion for the particles you don’t take into account their mass, for acceleration. Am I overlooking it or has it been omitted?

      • Tom Madams Says:

        You’re right, I apparently forgot to account for mass in the viscosity force calculation.

        As for surface tension and pressure relaxation, I actually calculate the acceleration directly (rather than the force), which is where the division by pi.m comes from in lines 275, 276, 281 and 282.

        That being said, it’s entirely possible that I’ve made a mistake somewhere as my maths isn’t terribly good either :/

        By the way, there were a couple of mistakes in my previous reply, which I’ve now corrected. Most notably, the extra integration term should be 2*pi*r, not pi*r.

  2. Niall Says:

    Sorry, wrong link, that’s an even older version, here’s the correct one http://niallmeister.deviantart.com/art/Viscoelastic-Jelly-183686746

    Also sorry for inadvertently plugging my stuff, haha.

  3. Niall Says:

    Hey, I’ve implemented your method of doing things myself in flash. I haven’t experimented with it much as I’m posting this comment about 10 minutes after getting a working build because I’ve noticed an issue. The mass/density contrast, aargh!

    Here’s a demo with a timestep of .1, with 2 fluids of mass 1 and 30 respectively – http://www.fileize.com/view/c82523d7-915/ – mouse to interact, sorry that it’s a bit like watching paint dry at the moment.
    I also have yet to try varying rest densities, that could possibly contribute to the effect as well.

    Do you think I might’ve coded this wrong or have you experienced the same thing?

  4. kaoD Says:

    Beautiful videos! I got all scared of SPH when I faced it the first time, I might give it a try.

    I tried your code to see if I could make a video generator out of it, for the sake of entertainmient, but particle radius below 0.04f makes it crash. How do you make the videos then, if it crashes for small radius? With a huge rendering area?

    • Tom Madams Says:

      I never got any hard crashes when running it locally. Then again, the code I posted is a much simplified version of what I actually use at home, so it’s possible I introduced a bug somewhere.

      In the posted videos, the particle radius is unchanged at 0.05 and kViewWidth is increased to something like 80 or more if I remember correctly, which zooms the camera out.

      I hope this helps; I’d love to see your videos if you come up with anything cool.

      • kaoD Says:

        With crash I mean unconsistent behaviour. All particles explode and go crazy.

        Anyways I think you might’ve introduced a bug, everything tends to splash a lot, even sudden splashes out of nowhere at the bottom of the field. Particles tend to flow like marbles or sand too. I’ll play with constants and see if I can fix it.

        I guess the part you simplified is the drawing part, right? I can clearly see the spaces between adjacent particles, so I guess the real code does some smoothing.

        • Tom Madams Says:

          Ah that’s not a bug, it’s what I meant by “Unfortunately, SPH simulations have a tendency to explode if the time step is too large” – it’s caused by forces in the simulation being too large for the integration method.

          You could replace the simple forward Euler with some more complicated integration scheme, but the simplest thing to do is to increase the number of substeps per frame. For some of the videos, I was using around 20 substeps.

          Alternatively, if you don’t mind your fluid being a little uh, compressible, you could also reduce the stiffness constants.

          • kaoD Says:

            Nice, got it to work and had some fun tweaking constants, but it takes so long per frame after a large number of particles is spawned, and re-running it from the beggining after each tweak is a bit time consuming.

            I’m planning on doing this realtime or at least a bit faster, implementing a quadtree for the neighbor check and Leapfrog integration to reduce substeps. If I succeed in doing so I will let you know :)

            By the way, was I right when I guessed that the simplified part is the drawing one? A big number of particles would make it look continous like in your videos, or it had some kind of near-deformation to make it look liquid?

          • Tom Madams Says:

            I’m looking forward to seeing what you come up with :)

            The rendering code is actually more or less unchanged – you could always increase the size of the constant in glPointSize if you feel the points are too small.

            The major simplification in the code I posted is the removal of multi-threading.

            The code as posted is set up to run on multiple threads; if you’re only running on one thread, you could rearrange it and probably speed it up by around 30%. For example, applying forces to both particles involved in an interaction at the same time.

  5. Mikhail Says:

    Hello, Tom!
    I would like to talk with you about SPH fluids, could you drop me your skype?

  6. yanioaioan Says:

    Hi there,

    i was lookin gat your code above and i can’t seem to find where do u implement the Tait’s equation of state to calculate pressure (as it was done in 2007 paper by Becker et. al ). I am having an sph algorithm where i firstly calculate my densities based on neighbour search and then i try to use Tait’s EOS so as to calculate pressure but it’s explodes. i see you use the simple gas law as i did.

    Any ideas to share would be great!

    Cheers

    • Tom Madams Says:

      My code is based mostly on the Particle-based Viscoelastic Fluid Simulation paper by Clavet et al, which IIRC doesn’t use Tait’s equation. Instead they opted for more a approximate model that tends to explode less.

  7. Jumes Says:

    Hi,
    I run your code, then I find there is large space between each particle.It seems that it is just a circle.Your demo video looks better.Building surface tringle mesh? How did you render the water surface?

    • Tom Madams Says:

      For the most part, the videos are rendered using many more particles than in the demo.

      To improve upon simple circles for rendering the particles, one common technique is to render gaussian splats for each particle into a off-screen accumulation buffer and then perform a per-pixel thresholding operation when rendering the accumulation buffer to the back buffer. Another technique that was used PixelJunk Shooter to hide some of the large gaps (especially when the fluid was falling) was to decay the accumulation buffer by some amount each frame, rather than clearing it.

      I don’t have any references to hand for either of these techniques, but hopefully that should give you pointers.

  8. Bakatare Says:

    This may sound awfully retarded, but what exactly is the gradient of the kernel function in the very first equation here (e. g. nabla W)? I mean, W is not a scalar field, just a simple function of one argument (the distance). Must be some gap in my math, have read several papers on SPH, it’s used everywhere and I totally don’t get it. Could you, please, explain how to compute that?

    • Tom Madams Says:

      Oh dear, I’m afraid it’s been a bit too long since I wrote this post and I’ve forgotten all the details :(
      I didn’t actually use the pressure model that the first equation describes, (I use the model from Particle-based Viscoelastic Fluid Simulation) so that gradient doesn’t appear in the source code I provided, which is why that might be confusing.

      The the thing to remember is that although W might be defined in terms of distance from particle i, it still forms a scalar field. Think about a particle i in 2D space: W(r_i, r_j, h) forms a non-zero scalar field in the region of space described by a circle of radius h centered at r_i. I hope this helps at least a little.


  9. […] but getting the code to run stable without explosions and crashes is far from simple). Then I found this blog post by Tom Madams that had sample source code. With that I was finally able to get SPH working. It can […]

  10. Me Says:

    Hi Tom,

    I just wanted to say thanks for sharing your code. It helped me a lot getting my own SPH working.

    http://softologyblog.wordpress.com/2013/09/17/multiphase-smoothed-particle-hydrodynamics/

    Regards,
    Jason

  11. Rhody Lugo Says:

    Pretty cool stuff!
    I just wanted to ask it your simulation models the incompressiblility of the liquids.
    Regards.

    • Tom Madams Says:

      No, the simulation makes several approximations in order to provide more stability at larger time steps. This is why you can see some compressibility in the videos.

  12. JPK Says:

    recently im trying to develop my own version of this in python. somehow this code doesnt make use of the standard formulas used in most of the papers. i have a few questions – maybe someone who is interested in sph can contact me so we can solve it together. probably the author of this isnt working on sph for a few years now :D

    1. The Kernels – why do they need to be normalized. and why are there often normalization factors for 2d and 3d given? isnt the Kernel basically always a 1d function? some even say they need to be symmetric W(r,h)=W(-r,h) – but as r is the Euclidian distance this doesnt make any sense to me.

    2.How can you convert pseudo densities and pressures to real-life quantities? I wonder how you can compute a real density by adding weighted masses together??

    3. Someone above postet a similar question regarding the gradient. as i only have a simple function – is it correct that i just calculated the direction vector between particle i and j. and then i multiplied the vector with the result of the W-function (1. derivative). so i have the direction pointing away from the calculated particle and the length is the value of the function. later i added all vectors together to get the final gradient?

    4. are there papers regarding compressible fluids like air? ive only seen water :(

  13. crococode Says:

    Hey Tom! Thanks a lot for your code posted here. I used your basic layout to get my implementation started. Still fighting with some surface tension issues for varying particle sizes, but the first results were really good:

    http://crococode.wordpress.com/2014/04/11/smoothed-particle-hydrodynamics-even-more/

  14. hassan Says:

    Hi
    Have you used the kernel gradient correction term in your code? please

  15. Golzar Says:

    Hi Tom
    For interaction between two phase (different particles with different density ) what did you do it?
    Regards


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

Follow

Get every new post delivered to your Inbox.

%d bloggers like this: