CFD Online Logo CFD Online URL
Home > Forums > Main CFD Forum

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

Register Blogs Members List Search Today's Posts Mark Forums Read

LinkBack Thread Tools Display Modes
Old   March 11, 2014, 15:02
Default Basic and simple SPH-code questions (smoothed particle hydrodynamics)
New Member
Join Date: Mar 2014
Posts: 3
Rep Power: 5
Peacekeeper is on a distinguished road
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 im 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).

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! )
- Collision, viscosity and not using euler integration is scheduled for later…
- Thanks to numpy a fast nearest neighbor algorithm is already integrated

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:

With a spiky kernel (q = r/h for normalization)

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:

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:

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

Peacekeeper is offline   Reply With Quote

Old   March 12, 2014, 16:16
Default density
New Member
Join Date: Apr 2011
Location: Virginia
Posts: 8
Rep Power: 8
lbrieda is on a distinguished road
Send a message via Skype™ to lbrieda
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:
lbrieda is offline   Reply With Quote

Old   March 12, 2014, 17:29
New Member
Join Date: Mar 2014
Posts: 3
Rep Power: 5
Peacekeeper is on a distinguished road
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:

now i dont really get why there is this previous term "15/pi*h^6".

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 ?

Last edited by Peacekeeper; March 12, 2014 at 22:08.
Peacekeeper is offline   Reply With Quote


Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On

Similar Threads
Thread Thread Starter Forum Replies Last Post
[ICEM] Meshing adjacent wall geometry and simple ICEM questions everdimension ANSYS Meshing & Geometry 25 June 20, 2012 04:25
LES basic and simple questions. dshawul Main CFD Forum 0 October 4, 2010 10:56
Basic Blade Design Code apoorv Main CFD Forum 0 April 5, 2010 06:37
Progressing to write my own simple code Jonny6001 Main CFD Forum 5 October 19, 2009 15:09
A simple CFD code for teaching basic CFD? Christoph Lund Main CFD Forum 13 September 14, 2005 04:36

All times are GMT -4. The time now is 08:19.