Higher order terms of Range Kutta for SPH?

[muffinman123]

I am using smooth-particle hydrodymamics(SPH) method to model a deformable object. It is based on the procedure found in Markus Gross's point based animation book, the algorithm is attached in the image below.

This algorithm uses Euler for numerical integration as seen in step 21 and 22, and I am trying to see if I can implement Runge Kutta method. However, when I look at the procedure I couldn't see how I can calculate the higher order terms.

in Range Kutta
y' = f(t,y)
k1 = euler's method
k2 = hf(t + h/2, y + k1/2)
k3 = hf(t + h/2, y + k2/2)
k4 = hf(t + h, y + k3)

in SPH
u'(t+1) = u'(t) + dt * F(t)/m

I see that f(t,y) should be the acceleration, but in SPH acceleration is obtained from forces summation from neighboring elements.

in this page
http://gafferongames.com/game-physics/integration-basics/

it describes the f(t,y) to be mass spring approximation of acceleration

a = -k * x - b*v

but is this suffice to use this replacement to calculate k2, k3, and k4?

what other ways can be used to approximate the acceleration?

 

[AlephZero]

The way to do this using RK is:

First find the values of k1 for each the particle in the model.

Them to find k2 for each particle, you assume each particle is at position y + k1 / 2, sum up the forces, and divide by the mass to get the acceleration. You can interpret this "physically" as finding how the forces have changed half way through the Euler step.

Then to find the values of k3 you assume each particle is at position y + k2/2, and so on. That is not so easy to explain "phyiscally", you just have to accept that mathematically it's another correction to the forces to improve the accuracy of the method

And similarly to find all the values of k4.

Finally, you find the position of each particle at the end of the step from the RK formula involving k1, k2, k3, and k4.

Presumably, in this simulation nothing is explicitly dependent on time, so the fact that you calculate k2 k3 and k4 at particular times (k2 and k3 half way through the step and k4 at the end) doesn't make any difference to the calculations. But if you were applying some external forces that varied with time, you would need to use the correct times for each "sub-step" of the RK calculation.

 

[muffinman123]

what other numerical methods out there that doesn't require looking ahead in this manner?

because my simulation is running on real time with random external forces applied at different times.

I tried leapfrog but the displacements could not reach equilibrium and the model blew itself up so I had to settle for Euler.

 

[AlephZero]

You could try Verlet. That seems to be popular in the particle simulation world.

I don't really see why you are describing R-K as "look-ahead". It is just doing 4 steps of calcualtion before you decide what the answer is.

No numerical method is going to be accurate for changes in loads that happen at a higher frequency than the length of the time steps. If the loads are changing with frequencies up to F Hz, in practice you probably need at least 4 or 5 steps to cover one cycle of the loading, i.e. time steps of about 1 / 4F or 1 / 5F seconds.

It may not make much practical difference if you assume your "random loads" remain constant for the length of each step.

 

[muffinman123]

isn't leapfrog a version of verlet? I already know leapfrog is unstable in my application, how can I adjust it to forcefully stabilize it?

 

[AlephZero]

Unless I'm missing something, Verlet seems to be exactly the same as what used to be called the central difference algorithm (and "leapfrog" would be make sense as a description of the same method).

If that is true, it is only conditionally stable, so reducing the time step might fix the problem. If that doesn't work, you could introduce some artificial damping, for example see http://en.wikipedia.org/wiki/Verlet_integration#Applications, but of course that would take energy out of the model which may not be what you want.

On the other hand, if the cause of the instability is your "random" applied loads, you might need to reduce the amount of "randomness" and take out some of the high frequency content of the loading. I would start with the sub-problem of getting the model running with the particles interacting, but without the "random" input.

Euler is also conditionally stable and only first order accurate, so if it is "working", the reason may be because the low accuracy is (maybe by good luck rather than design!) filtering out whatever is causing the problem when you try to get a more accurate solution.

 

[muffinman123]

is damping usually introduced in velocity or displacement?
if displacement reaches stable condition first, wouldn't velocity still be fluctuating because velocity calculation does not involve displacement at all?

or should the damping go into the force calculations in SPH?

 

[AlephZero]

You can estimate the velocity at time t from (x(t+h) - x(t-h))/2h and include a force term proportional to that. When you rearrange the equation, you get x(t+h) = a linear function of x(t) and x(t-h), similar to Verlet without damping. That's where the equation in the wiki page comes from.