Table of Contents Show
Attractive forces
The simulation already looks very good at this stage and there is also a whole range of setting options and freely definable parameters. However, the interaction between the particles has not yet been dealt with. I have therefore been looking for a way to model this behaviour mathematically. The equations I found were quite complex at the time. I was looking for a formula that was sufficiently accurate but also easy to implement.
After some (what felt like an eternity) time, I came across the Lennard-Jones potential. This principle has its origins in chemistry. There it is used to calculate the interactions between molecules or atoms, e.g. to visualise the so-called Van der Waals forces. This is done in a very simple but surprisingly efficient way. The term “simple” is of course very subjective. However, if you bear in mind that theoretical chemistry is ultimately quantum physics, then the Lennard-Jones potential is astonishingly simple.
The potential calculates repulsion, attraction and equilibrium, i.e. precisely the three states that are crucial for particle simulation. To avoid a brute force approach that takes all particles into account, the neighbour search is always a good option. Only particles within a certain radius of the particle currently under consideration are taken into account. The Lennard-Jones equation itself is
V(d) = 4 * epsilon * [(sigma / d)12 - (sigma / d)6]V(d) is the potential energy as a function of the distance d between two particles. In the script, this corresponds to the force acting on neighbouring particles. The factor epsilon corresponds to the force of attraction of a particle. sigma is the distance at which the force becomes 0. In VEX the formula looks like this:
lj_potential = 4 * epsilon * (pow(sigma / d, 12) - pow(sigma / d, 6));Since I was ultimately unable to find an empirical or experimental value for sigma, I set it to 1, which also simplifies the formula, because (1 / sigma)12 can also be written as sigma-12.
With equations of this type, it is also often the case that the forces become extreme at very small distances. To prevent this, the values in VEX can be limited using the clamp() function. To do this, you define a minimum and a maximum – this ultimately corresponds to the epsilon (here: force_limit) in the above equation. I have already replaced sigma with distance here.
float force_magnitude = clamp(4 * (pow(distance, -12) - pow(distance, -6)), -force_limit, force_limit);Here’s to good neighbours!
I now had to apply this equation to each neighbouring particle of the particle currently under consideration. The search for neighbours is very nicely solved in VEX and there are several functions for this purpose. I decided to use nearpoints(). However, it took me some time to get this function to work at all. This was due to a small, easily overlooked peculiarity in the POP system.
The input tab of the POP Wrangle has a drop-down menu for each virtual input from 1 to 4. Only by selecting the correct option is it possible to use the neighbouring functions of VEX in the wrangle. If the wrong option is selected, the number of neighbour particles is always zero! The only option that works in this context is “Myself” at Input 0.
To call the neighbouring particles in sequence, create a loop. There are also several options for loops in VEX, but the following variant is the most suitable here.
foreach (int neighbour; nearpoints(geoself(), @P, interaction_radius)) {
//Calculate the force and apply it to the particle
}In the loop, a counter called neighbor is initialised and since there are no “broken” particles, neighbor is an integer. nearpoints() captures the neighbouring particles at the current position @P in the radius of interaction_radius. The geoself() function refers to the above-mentioned input 0.
What is still required is the distance between the current particle and the neighbours as well as the direction in which Lennard-Jones should act. All these calculations take place for each neighbour within the foreach loop.
float interaction_radius = 0.75; // Search radius for neighbour particles
float force_limit = 0.025; // Min and max force
foreach (int neighbour; nearpoints(geoself(), @P, interaction_radius)) {
vector neighbour_pos = point(geoself(), "P", neighbor);
vector direction = normalize(neighbor_pos - @P);
float distance = length(neighbour_pos - @P);
// Lennard-Jones implementation
float force_magnitude = clamp(4 * (pow(distance, -12) - pow(distance, -6)), -force_limit, force_limit);
vector neighbour_force = force_magnitude * direction;
v@force += neighbour_force;
}In the loop, the position of the neighbouring particle currently under consideration is determined and finally normalised. This function transforms a vector so that its length is always 1 – the result is a so-called unit vector that represents a direction, not a magnitude. And the distance is calculated from the difference between the position of the particle and a neighbour. I was then able to use this distance in the Lennard-Jones equation. Since attraction and repulsion are directed forces, I still had to include the direction vector in the determination of the neighbour_force.
The summation of the force (v@force = neighbour_force;) also takes place in the loop, because only the neighbours are treated here. However, the force total_force, which acts on the current particle, is still calculated outside the loop.
// Calculate total force
vector total_force = (stokes_force + buoyancy_force + wave_force + turbulence_force) * damping + v@force;Here, v@force also reappears from the loop for the neighbouring particles, as the current particle is of course also affected by Lennard-Jones. However, as the potential does not depend on the depth, the force is added outside the bracket.
The value for interaction_radius depends on the number of particles in the simulation domain. A radius of 1 or more is appropriate for a small number of particles. If there are a large number of particles, the value can be set correspondingly smaller.
