Hi!

Iam trying to develop a very basic 2d SPH program in python and it seems that in this forum there are actually some people who know more about this topic^^

I have read a ton of literature but as English is not my main language im running into a few problems here. Right now I created a basic SPH program that calculates the movement based purely on the pressure gradients. I think my main problem is the understanding of the basic kernel interpolation.

I see the same formulas everywhere but sometimes I wish someone supplied a basic example with values to compare with.

to calculate the density of a specific particle i need to use this:

1. 2D Kernels:

I don’t really understand what they normalization factors (15/pih^6) are for and no kernel that is supplied online does actually equal to “1”. Some are even called the same but the volume created by the kurve (for 2D) is not even close to one. Also i found different formulas for the same kernel. the Spiky Kernel can be written as “factor * (h-r)^3” or “factor * (1-(r/h))^3”.

“factor * (h-r)^3” with the factor being “10/(PI()*I3^5)” does actually provide the same volume even if “h” is changed. So I assumed it should be correct. The value is actually ~ 4.6 (not 1!)

I tried it in excel for 2d! i hope i did it the right way. maybe my whole understanding of how it works is wrong...

As I use different Kernels for calculating the density (poly6) and pressure (spiky) I think the totally different values do fu** up my calculations. Using poly6 my density estimation is around “455” while using the spiky kernel leads to densitys in the range of millions! This brings me to my second question

2. People talk about using the rest density ( for “p = (p – p0)” ) like a real density! They use 1000 kg/m³ for it. But isn’t the density basically a pseudo density? How can I compute realistic densities and values out of a summation of weighted mass points? Like I said – my rest density is highly influenced by the kernel and is between 100 and 10000000

attachments:

- simple excel sheet in rar. maybe someone can edit it to be correct?

it goes all boom:

no viscosity so kinetic energy is preserved?

hm?

thanks for the input! the idea with the "invMass" makes sense to me and i will try it. but it still doesnt answer how people can archieve physical values ?

regarding the volume of the kernel. i hacked it together in excel like this:

the idea was to multiply the surface below one side of the function with "2* PI * r" - (Rotational Symmetry)

1. i discretized my values from min (e.g. "0") to max (e.g. "1") into 50 points.

2. calculated W(r,h) with the formula at all points

3. simple little rectangles for the integral ( r0 * W(r0,h) + (r1-r0) * W(r1,h) + (r2-r1) * W(r1,h) + (r3-r2) .....)

4. multiplied the sum with the 2 * PI * r with r beeing equal to h.

i agree that the excel sheet is a bit unclear ^^ i attached a slightly improved version and here is a gif so you dont have to download it.

when you look at the red box for the volume it stays the same around ~4.6

i really appreciate your detailed help. thanks alot for fixing the ecxel. i hope i dont annoy you with my frequent questions ;) i found my mistake in the excel sheet. obviously the functions where alright but my calcualation of the integrals wrong.

now its starting to flow! iam quite exited ^^

note:

- no viscosity yet - instead iam simply dampening the velocity by -0.1 in every timestep

- walls just move the particle back to their surface if crossed

- very low stiffnes of the fluid. so its compressible

ok i think the next step is a closer look at what you wrote about calculating the rest density tomorrow. right now i choosed working random values. sorry that iam a bit slow here. can you give an example for a 0.1m x 0.1m space with 100 particles and water ? your method seems to need an additional first loop for initialising the simulation?

edit: i think thats indeed what lets my simulation require such a tiny timestep (t=0.0003) or it explodes. as initially the particles are too close or to far away from each other HUGE accelerations are calculated. after it slowed down i quess i could change the timestep to be quicker.

Hi!

i was rather busy the last days so my progress isnt that impressive, but i got your method to work today! works like charm.

its still rather slow due to python. in the future when everything is setup i might transfer everything to cython to speed it up. but thats not a problem right now. next comes the artificial viscosity.

i have a few general questions left here:

- all papers simulate water. but as sph is actually compressible what about an ideal gas? is that possible? actually ideal gas doesnt even have viscosity but probably that´d make the whole simulation quite unstable.

- i noticed that i use the spiky kernel for the pressure calculation AND the density calculation. as you can see above it works quite good. using the poly kernel for the density like its done in the literature results in clumping of the particles. yes, they clump together altough the spiky(-pressure) kernel should avoid that. i have no idea how this can happen.

Hi there Peacekeeper. To avoid clumping, it might be a good idea to apply a repulsive-only kernel between each interacting particle pair, where force is proportional to either distance or a local pressure value ony dependent on those two particles. As far as I know (I'm no expert but have some experience with sph) it's the only way to avoid clumping. In my experience it works best if the particle pair kernel is spikier than the pressure kernel, which means that repulsive force is usually weaker than pressure force, except when particles are very near each other.

Cheers Mike

yep and that should be the case here. as i can use the poly 6 kernel for calculating the density and pressure at each particle but for the pressure force (gradient) i use the spiky kernel. so if particles come close to each other there should still be a repulsive force

Well okay, might just be a matter of fine-tuning then. What happens if you turn off pressure kernel, so the only forces influencing the particles are from the repulsive-only kernel?