## 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:

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

where 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:

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:

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

with

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:

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] = π 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:

January 5, 2011 at 9:17 pm

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?

January 5, 2011 at 10:45 pm

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.

January 6, 2011 at 2:23 am

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?

January 6, 2011 at 11:04 am

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.

January 5, 2011 at 9:43 pm

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.

January 8, 2011 at 12:40 pm

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?

February 8, 2011 at 11:03 am

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?

February 8, 2011 at 11:13 am

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.

February 8, 2011 at 1:13 pm

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.

February 8, 2011 at 1:36 pm

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.

February 12, 2011 at 8:11 am

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?

February 12, 2011 at 10:46 am

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.

June 21, 2011 at 8:21 am

Hello, Tom!

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

June 21, 2011 at 4:53 pm

Hi Mikhall, I’m afraid I don’t have Skype (or much free time these days). I’ve dropped you a mail.

June 24, 2011 at 8:08 am

Thanks! I was sent some questions to you via e-mail, hope you’ll react as soon as you can.

April 3, 2012 at 9:10 am

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

April 4, 2012 at 6:47 pm

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.

June 5, 2013 at 4:07 pm

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?

June 6, 2013 at 8:21 am

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.

July 7, 2013 at 5:23 am

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?

July 8, 2013 at 10:42 am

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.

September 16, 2013 at 6:53 pm

[…] 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 […]

September 17, 2013 at 8:13 pm

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

September 17, 2013 at 9:13 pm

No worries, I’m glad you found it some help.

January 10, 2014 at 3:20 am

Pretty cool stuff!

I just wanted to ask it your simulation models the incompressiblility of the liquids.

Regards.

January 10, 2014 at 9:10 am

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.