Basic and simple SPH-code questions (smoothed particle hydrodynamics)
 

Hi!
This is my first post in this forum, so I really appreciate any help. This problem is driving me crazy in the last days. And sorry for my English …
As a hobby i´m trying to improve my python (3.3) skills and my knowledge in fluid dynamics. Recently I discovered some interesting paper concerning SPH (Smoothed Particle Hydrodynamics).

http://software.intel.com/en-us/arti...-games-part-15
http://obswww.unige.ch/lastro/confer...f/lecture4.pdf

The Internet is full of it, so I decided to create my own very basic version of it. My goals are:

- Acceleration only based on the pressure gradient for now (tough enough! :D )
- Collision, viscosity and not using euler integration is scheduled for later…
- Thanks to numpy a fast nearest neighbor algorithm is already integrated :D

But after getting the first (rather unstable) results I begin to wonder if I really understand the math behind it. I will begin with the two most basic things to do in sph - get the density and pressure with the kernel.
Regarding the density I used this equatation:
http://110.imagebam.com/download/yJG...13634895/1.jpg
With a spiky kernel (q = r/h for normalization)
http://109.imagebam.com/download/_Da...13634896/2.jpg
After blindly programming everything the results worry me:

How can I compute density by rather random weighted functions? How can I compute density by adding weighted point masses together? In the end I still have a mass. If I increase the range h I get information from more and more particles so the “density”-value will rise. Some kernels return “1” for zero “q” (= maximum value of W) some return random values like “3.92”. So I get totally random values for my particle density. In my head as an (unprofessional) example it looks like this:
http://104.imagebam.com/download/6Ly...13637142/3.jpg
Shouldn’t I somehow tie the accumulated masses to the observed volume (or surface) to get the density? But I have not seen that in the codes I looked at:

"dist = particle.Position - ((FluidParticle) particles[nIdx]).Position;
particle.Density += particle.Mass * this.SKGeneral.Calculate(ref dist);”

Most authors use this to calculate the pressure afterwards:
http://111.imagebam.com/download/ToJ...13634898/4.jpg
P_0 ist called the rest density and usually 998 kg/m³ for water. Then it really stops making any sense to me:

Pi=k * (0.0217 kg(?) -998 kg/m³ ) .. äh what?

Result is obviously a totally wrong pressure.
My next question is how to use a gradient when I have no scalar-field but a simple function W(r,h). I think I did it the right way but before I go into that in detail I hope my most noobish questions can be answered ;)

Greetings
JP

density
 

JP, what you describe is what is in fact done in the Particle In Cell method (i.e. my article), where the particles are described by mass and velocity, and you obtain the density by dividing the weighted sum by the cell volume.

How did you implement your kernel function? The second link you posted shows how the volume cancels out, by the dV term in the integral replaced by m/rho.

Here is another good intro article on the kernel interpolation: http://www.icpsr.umich.edu/CrimeStat...tChapter.8.pdf

probably my understanding of the way how to use a kernel is truly wrong.

i made a very basic approach and used the functions provided in the little box:
Spiky Kernel: (1-q)^3 with q beeing r(distance)/h(smoothing length)
i treated all particle positions as 2d vectors so for example p1(1,1) and p2(3,2). so i used cartesien coordinates and calculated the distance between those - placed the result into the simple function and thats it - got the value for that particular particle - e.g. (1- sqrt(2^2+1^2) / h ))^3 . after adding all results of the particles in range together i got my "density" (each multiplied by its mass like posted above)

thats all rather simple and i think iam missing something mathematically.
most authors actually dont even supply the kernel function in such a simple manner as posted above. usually they are provided like this:
http://109.imagebam.com/download/Cx1...13923868/1.JPG
now i dont really get why there is this previous term "15/pi*h^6".

edit:
ok the kernels i used before missed the normalization term (15/pi*h^6) so the integral didnt satisfie the unity condition.
now i wonder why there are different normalization terms for 2d and 3d. isnt the kernel basically always 1D ?
http://imageupper.com/s02/1/4/K1394658893236966_1.jpg