A particle settling through a fluid, live¶

Stokes's law for a small sphere at a low Reynolds number, the speed it sinks at through a still fluid, its weight less its buoyancy balanced by the drag: w = (ρₚ − ρf) g d² / (18 μ), downward where the particle is denser than the fluid and upward where it is lighter. Below, a grain of fine sand's size, a tenth of a millimetre across, in water, is swept from 500 kg/m³ lighter than the water to 2000 kg/m³ denser, 25 kg/m³ a step, twenty steps a second, and at each the library's law (lusine.laws.settling) gives its settling speed on your GPU as you read. Beside it the same program computes the two forces the speed must balance: Stokes's drag on a sphere moving at it, 3π μ d w, and the grain's weight less its buoyancy, (ρₚ − ρf) g π d³ / 6.

This page runs the program on your GPU through WebGPU, so it needs a browser with WebGPU (a recent Chrome, Edge, Safari or Firefox). Without it the program does not run and nothing is drawn below.

What it shows¶

The strip above the controls is green while every check holds, and turns red for good once one fails.
The grain: the first chart draws the settling speed against how much denser than the water the grain is, a point a step over one whole sweep: up where it is lighter, down where it is denser. The second draws the drag at that speed and the weight less the buoyancy, one over the other. Below them, the difference, the speed and the two forces.
The checks, each held at every step: the forces balance, the drag within 0.001% of the weight less the buoyancy; denser sinks, the speed down exactly where the grain is denser than the water and up exactly where it is lighter; neutral hangs, a grain as dense as the water not moving at all.

Every number here is computed by the program on the GPU; the page only reads the values back and draws them. The water's density and viscosity, gravity, the grain's diameter and the sweep are the program's settings, which the host may write; this page offers no control to change them yet.

The program¶

The sweep, the law's answers and the checks, the panels the page draws, and the one procedure, main, run each step: the sweep moves on, the law is asked at the grain's density and at the water's, and the checks are taken.

// Stokes's settling law's live page: a grain's density swept from lighter than the water it is in to denser, the
// library's settling speed at each, beside the force balance the speed must strike, with the page's checks computed
// here on the GPU and shown by the panels the program declares.

namespace settling;

lusine.version version = { { number = 1 } };

// The reader's settings: the water's density and viscosity, gravity, the grain's diameter (fine sand's), and the
// sweep: how much denser than the water the grain starts and how far each step moves it.
host struct {
  float<kg/m3> fluid = 1000.0 kg/m3; float<Pa*s> viscosity = 0.001 Pa*s; float<m/s2> gravity = 9.81 m/s2;
  float<m> diameter = 0.0001 m; float<kg/m3> lowest = -500.0 kg/m3; float<kg/m3> stride = 25.0 kg/m3;
} settings;

// The sweep: the step it is at, from 0 to 100 and round again, starting where the grain is as dense as the water, and
// how much denser than the water it is there.
struct { integer step = 19; float<kg/m3> difference; } sweep;

// The settling speed at the sweep, downward positive; the drag Stokes's law puts on a sphere moving at it,
// 3 pi mu d w, and the grain's weight less its buoyancy, (rho_p - rho_f) g pi d^3 / 6.
struct { float<m/s> speed; float<N> drag; float<N> net_weight; } law;

// The checks: the drag balancing the weight less the buoyancy within 0.001% at every density; the speed's sign the
// density's, a denser grain sinking and a lighter one rising; and a grain as dense as the water hanging in it.
struct { bool balance = true; bool sign = true; bool neutral = true; } checks;

// A strip of lamps the page draws, lit green while every check holds and red once one fails.
struct lamp { vector<float, 4> colour; };
lattice<lamp, 2> lamps(spacing = 1 m, boundary = wall);

vector<float, 4> lamp_colour(bool holds)
{
  return holds ? vector<float, 4>(0.15, 0.6, 0.3, 1.0) : vector<float, 4>(0.75, 0.15, 0.15, 1.0);
}

float<N> stokes_drag(float<Pa*s> viscosity, float<m> diameter, float<m/s> speed)
{
  return 3.0 * 3.14159265 * viscosity * diameter * speed;
}

float<N> net_weight(float<kg/m3> difference, float<m/s2> gravity, float<m> diameter)
{
  return difference * gravity * 3.14159265 * diameter * diameter * diameter / 6.0;
}

output table<lusine.widget> panels = {
  lusine.series { { panel = "the grain", chart = "settling speed", quantity = law.speed, steps = 101,
                    against = sweep.difference, label = "w" } },
  lusine.series { { panel = "the grain", chart = "forces", quantity = law.drag, steps = 101,
                    against = sweep.difference, label = "drag at w" } },
  lusine.series { { panel = "the grain", chart = "forces", quantity = law.net_weight, steps = 101,
                    against = sweep.difference, label = "weight less buoyancy" } },
  lusine.value { { panel = "the grain", quantity = sweep.difference, label = "denser than water by" } },
  lusine.value { { panel = "the grain", quantity = law.speed, label = "settling speed" } },
  lusine.value { { panel = "the grain", quantity = law.drag, label = "drag" } },
  lusine.value { { panel = "the grain", quantity = law.net_weight, label = "weight less buoyancy" } },
  lusine.value { { panel = "the checks", quantity = checks.balance, label = "the forces balance" } },
  lusine.value { { panel = "the checks", quantity = checks.sign, label = "denser sinks" } },
  lusine.value { { panel = "the checks", quantity = checks.neutral, label = "neutral hangs" } },
};

void main()
{
  sweep.step = sweep.step >= 100 ? 0 : sweep.step + 1;
  sweep.difference = settings.lowest + float(sweep.step) * settings.stride;
  law.speed = lusine.laws.settling(particle = settings.fluid + sweep.difference, settings.fluid, settings.gravity,
                                   settings.diameter, settings.viscosity);
  law.drag = stokes_drag(settings.viscosity, settings.diameter, law.speed);
  law.net_weight = net_weight(sweep.difference, settings.gravity, settings.diameter);
  checks.balance = checks.balance && abs(law.drag - law.net_weight) <= 0.00001 * abs(law.net_weight);
  checks.sign = checks.sign && (law.speed > 0.0 m/s) == (sweep.difference > 0.0 kg/m3)
                && (law.speed < 0.0 m/s) == (sweep.difference < 0.0 kg/m3);
  checks.neutral = checks.neutral
                   && lusine.laws.settling(particle = settings.fluid, settings.fluid, settings.gravity,
                                           settings.diameter, settings.viscosity) == 0.0 m/s;
  lamps.colour = lamp_colour(holds = checks.balance && checks.sign && checks.neutral);
}

The law's other proof (the speed against a reference over sand and silt in water and air) is on its proof page.