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?

I don’t really understand what they are for and no kernel that is supplied online does actually equal to “1”.

How are you computing the volume? Because the volume of the kernel is its integral in the range [-h,h], and it has to be 1.

About 2: what you can do (at least that's what I did) was to set the mass of your particles so that their SPH rest density is your desired rest density. it's fairly easy, basically you do:

• create 1 particle P with a full neighborhood filled of particles around it
• set mass of all particles to one.
• compute the density of P, call it invMass.
• Since the density is proportional to the mass of your particles, set the mass of your particles is then 1/invMass * rest_density.

That way, your rest density is going to be the desired one.

PS: About going boom. SPH is highly unstable, really. Well, at least it is more unstable the more you want to conserve volume. In my (little) experience you'll need a well-tuned viscous force to compensate for your pressure forces to avoid explosions. I used Monaghan's.

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

Concerning the integral, I have attached another excel page so that you can verify the Poly6 kernel in 2d.

As for how to use physical quantities in SPH well, I don't understand your question...

I told you how to adjust the mass of the particles so that SPH tries to conserve the correct density of the fluid via the state equation (and the forces that arise from there...) In that regard, the SPH is going to be correct.

The mass of the particles indeed is going to depend on the kernel. However, keep in mind that your particles are ultimately samples, not "physical balls". So for example, if you want to know the mass of a box of fluid of 1x1x1 meters a SPH, you don't actually add up the masses of the particles that are inside¹. You get the particles that are inside, you get their densities, you average them (weighted by how much of their smoothing circle/sphere is inside of the box, if you want accuracy), and you multiply by your volume (1m^3). Since you are approximately conserving the density, your box will weight approximately 1000Kg.

Does that explanation help?

As for viscosity, that's a trickier one. The viscosity I used in the past was Monaghan's artificial viscosity, and it is called artificial for a reason (it has no physical meaning whatsoever). I found the proper values by tuning...

¹ I had seen in some lecture notes that you actually add up the masses of the particles of SPH to find the mass inside a volume. That is correct (if not really accurate unless you take the smoothing volume into account) with some methods for initializing the particle's mass. However, with those methods what you don't get is a physical rest density value. Keep in mind that there is no problem with that. The simulation is going to be the same as long as the ratio between your rest density and mass is the same. You are actually free to choose how you want to proceed.

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.

Glad to be of help!

About the rest density. Ok, what your SPH solver will do is to enforce your rest density via pressure forces. Indeed, whatever your state equation, it will give you pressures that depend on your current density and your rest density.

Now, as said, SPH is unstable, and its unstable the higher the forces it generates (thus, the higher the timestep or the higher your current situation deviates from your target situation, in which the fluid is not compressed). So what you want, is to avoid compression as much as possible at all times. This begins with your initial state.

You want to initialize your simulation in a correct state. So you have two options: for a given mass and rest density, put the particles in the right positions (which can be hell to do) or given some particle positions and rest density, compute which has to be the mass of the particles. How to do this?

Before the simulation (we will discard the particles we are going to create):

- Create a particle (let's call it P) and set its mass to 1 (or whatever value, for that matter).

- Create its neighboring particles. Use whatever algorithm you are going to use in your real simulation. If for example you create particles every N meters in X and Y uniformly, just do that. The closer P resembles any particle inside your final fluid volume the better.

- Now, using the same algorithms and kernels you are going to use in your simulation, compute the density of P.

Now you know, for your setup, whats the ratio between your density and your particle mass. Now, do you want to have a particular rest density (eg. water's)? Easy! The mass you have to use is currentMass * watersDensity / currentDensity.

As for your timestep, there are formulas to compute it automatically. Basically you want a timesteps so that your particles don't move in one timestep more than a fraction of the smoothing length. Although it is a little more involved than that. I am sorry I don't recall the formula right now. I'll look for it if you are interested.

PS: You're not going slow, don't worry. SPH is a simple idea but it is filled with little details that make life... interesting.

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?