What repulsion direction do we use for Smoothed Particle Hydrodynamics when radius is 0?

Question

When doing SPH, the paper by Kelagar recommends using a particular kernel for pressure induced forces between particles. The kernel it recommends is the following when the radius is within the kernel radius:

(15/(pi*h^9)) * (h - r)^3

where h is the kernel radius, and r is the radius we are interested in calculating the value of a function at.

The paper then states that the gradient of this function is

(-45/(pi*h^9))*((r_vec)/r)*(h-r)^2

where r_vec is now the vector from the center of the kernel to the point we are interested in. As length of r_vec goes to 0 from the positive direction, the paper states that this gradient approaches:

(-45/(pi*h^6))

But this is a scalar, not a vector. In order for there to be a repulsion between the two points we're interested in, there needs to be a direction to repel in.

What direction should we use for when two particles are right next to each other?

Answer 1

I assume that first expression is meant to be a potential. The negative gradient (derivative with respect to r) is then the force. This gradient is a vector, always pointing toward or away from the center. This appears correct for the second expression.

r_vec is, according to what you say, a vector pointing away from the origin to a point at some distance r away. (r_vec/r) is then a unit vector to specify direction. This works at every point except the origin itself, where it can be declared undefined, or declared to be zero. Zero is the average value of (r_vec/r) over all "nearby" points. This means zero force.

Normally in particle simulations with pair-wise forces, we ignore forces of a particle on itself, and of two particle at the same exact position. What about two particle very close, and you have a force law that goes like 1/r, 1/(r^2), or similar? Nobody wants a divide by zero fault. Usually there's a small radius below which the potential is constant matching the given potential formula at the boundary of that radius. Particles too close together have zero force, just so that the simulation won't crash. It may seem unphysical for the force to suddenly cease just inside that boundary when it is fiercely strong just outside it. But we strive to avoid such situations. Keep count of such incidences, and if there are too many, the simulation has gone bad. Maybe a smaller time step is needed.

Luckily you don't have a 1/r type of force, but still you have that nasty r_vec/r whose direction can swing wildly. The same technique of making force zero below a certain tiny radius will help.

But that third expression bothers me. If it's force at r=0, then starting with the force law in the second expression, I'm not sure how the third expression comes about. The problem of it looking scalar while, if it is supposed to be force, expecting vector could be resolved by understanding that it is meant to be the radial component of a force vector. Just multiply the expression by (r_vec/r), the familiar unit-magnitude vector. OTOH, it has no defined direction, so it is nonsense.

Better overall solution: start with a new potential function, one that smoothly levels off and is flat right at r=0, like exp(-r^2) or 1/(1+r^2). The given potential peaks sharply. You want something more like Instead of declaring force zero inside some small zone, the force would just naturally be zero at r=0. Find a flat-at-origin potential that approximates the given one well outside some small radius.