Marine-Snow-Simulation in Houdini

On the wave of success

What was still missing was an upward and downward movement that represents the wave motion on the surface. A sine function, which represents a simple wave, is actually always suitable for this type of movement. To define the wave, frequency and amplitude are required, i.e. the number of waves per time interval and the wave height in metres. The time dependency can be read off directly with @Time. The entire formula is therefore

float frequency = 0.2; // Number of waves per time interval
float amplitude = 1.1; // Wave height in metres
float oscillation = amplitude * sin(2 * M_PI * frequency * @Time);
vector wave_force = set(0, oscillation, 0);

The last line converts the oscillation into a vector that only acts along the Y-axis. In physics, forces add up and therefore the total force is

vector total_force = gravity_force + bouyancy_force + wave_force;

To allow this force to act on the particles, it is converted into a vector attribute v@force. The forces are totalled with each time step:

v@force += total_force;

With the correct values, the result of the simulation actually shows the desired behaviour, but is far too undifferentiated overall, i.e. ultimately totally boring. If you check @force in the geometry spreadsheet, you can see that the force changes over time, but all particles experience exactly the same force. The VEX script with the base forces now looks like this:

// Initial variables
float g = 9.81; // Gravitational acceleration in m/s^2
float radius_p = 0.005; // Particle radius in m
float rho_p = 1050.0; // Particle density in kg/m^3
float rho_f = 1030.0; // Fluid density in kg/m^3, e.g. water
float frequency = 0.2; // Number of waves per time interval
float amplitude = 1.1; // Wave height in m

// Intitial calculations
float volume_p = (4.0 / 3.0) * M_PI * pow(radius_p, 3);
float mass_p = rho_p * volume_p;
float oscillation = amplitude * sin(2 * M_PI * frequency * @Time);

// Calculate force components
vector gravity_force = mass_p * set(0, -g, 0);
vector bouyancy_force = set(0, (rho_f * volume_p - mass_p) * g, 0);
vector wave_force = set(0, oscillation, 0);
vector total_force = gravity_force + bouyancy_force + wave_force;
v@force += total_force;

Is it all just coincidence?

In nature, things are rarely uniform or uniform and there are always variations. Two parameters can be randomised here: radius_p and rho_p, i.e. particle radius and density. Both values influence the particle mass and therefore the force acting on the particles. In order to randomise values in VEX, a specific syntax is usually used.

float radius_p_var = 0.002; // Particle variance in m
float rho_p_var = 30.0; // Particle density variance in kg/m^3

// Generate random values for radius and particle density

float radius_p_rnd = fit01(rand(@id), radius_p - radius_p_var, radius_p radius_p_var);
float rho_p_rnd = fit01(rand(@id), rho_p - rho_p_var, rho_p rho_p_var);

The rand(@id) function generates values between 0 and 1, where @id is used as a seed to generate a different random number for each particle. In many cases, the point number @ptnum is used as the seed. However, as the number of particles changes constantly during the simulation, @ptnum does not provide reliable values. The fit01() function transforms the random numbers so that the values are between the two specified variables. For the density, I would then have random values between

rho_p - rho_p_var: 1050.0 - 30.0 = 1020.0
rho_p + rho_p_var: 1050.0 + 30.0 = 1080.0

This means that there is at least a small variation, but it is often the subtle effects that make the difference. A look at the Geometry Spreadsheet shows the changes.

Stokes to the rescue!

The realisation from my first attempts was that my model was too simple. After some searching, I was able to find a very promising concept: Stokes’ Law. The “inventor” of this law is George Stokes, who together with Claude Navier described the famous Navier-Stokes equations. These equations are the centrepiece of every fluid simulation. Stokes wanted to achieve a more realistic result, as the equation describes the behaviour of particles in a liquid medium. Stokes’ equation looks like this in its VEX formulation and will replace gravity_force as gravity_force assumes a vacuum.

float eta = 0.001; // Viscosity of water in Pa*s
float sink_velocity = (2.0 / 9.0) * (radius_p_rnd * radius_p_rnd) * (rho_p_rnd - rho_f) * g / eta;

The variable eta describes the viscosity of the surrounding medium in [Pa*s]. For water, the value is 0.001 Pa*s.

However, the problem now arises that the result of Stokes’ law is a velocity, while buoyancy is a force. The two cannot simply be combined. In order to convert the sinking speed into a force, it must again be dependent on the particle mass, because in principle F = m * a still applies!

vector stokes_force = mass_p * set(0, -sink_velocity, 0); // Stokes-based sinking force

The total force is now

vector total_force = stokes_force + bouyancy_force + wave_force;

As the velocities of the particles have increased significantly as a result of the new calculation, I have introduced damping. The damping is a factor and a value of 1 does not change the total force. This is therefore like a switch and the value can then be set < 1 if required.

float damping = 0.5;
vector total_force = (stokes_force + bouyancy_force + wave_force) * damping;

Strange currents

The differences between the versions may not be so obvious, but the simulation with Stokes appears more appealing, because the particles are now moving in a fluid. Now it was time to add an ocean current. After all, the particles of Marine Snow don’t just move up and down in the ocean! VEX provides several noise functions and my immediate choice was curlnoise(). To make the noise time-dependent, I added @Time to the position @P. And to convert curl_noise back into a force, a mass dependency is needed again. The following therefore applies

float noise_strength = 0.5; //Strength of the curl noise
vector curl_noise = curlnoise(@P @Time) * noise_strength;
vector turbulence_force = mass_p * curl_noise;

However, the dependence on the mass means that very large values are required for noise_strength in order to see any change at all. Values around 1000 would be appropriate here. In order to be able to calculate with reasonable values, I have abandoned the principle of physical consistency here and used a simplified calculation instead:

vector turbulence_force = curlnoise(@P * 0.1 @Time) * noise_strength;

The factor 0.1 is an empirical value that I inserted after some tests. This made the simulation look more harmonious. A good starting value for noise_strength was 0.5, but this value also changed in the course of the work. This new force is also added to force.

vector total_force = (stokes_force + buoyancy_force + wave_force + turbulence_force) * damping;