A particle settling through a fluid¶

The speed a small particle sinks at through a still fluid.

Stokes's law for a sphere at a low Reynolds number, its weight less its buoyancy balanced by the drag: w = (rho_p - rho_f) g d^2 / (18 mu), downward where the particle is denser than the fluid and upward (negative) where it is lighter. Past a Reynolds number rho_f w d / mu of about 1 the drag grows faster and the law overstates the speed; a game's grains of sand or silt in water are within it. Source: Stokes, Trans. Camb. Phil. Soc. 9, 8 (1851); Batchelor, An Introduction to Fluid Dynamics, section 4.9. Checked against a reference in numpy, which has its own tests.

Source

Stokes, Trans. Camb. Phil. Soc. 9, 8 (1851); Batchelor, An Introduction to Fluid Dynamics, section 4.9.

Proof

reference, invariant, limit.

settling.trait

namespace lusine.laws;

float<m/s> settling(float<kg/m3> particle, float<kg/m3> fluid, float<m/s2> gravity, float<m> diameter,
                    float<Pa*s> viscosity)
{
  return (particle - fluid) * gravity * diameter * diameter / (18.0 * viscosity);
}

Stokes's law, the speed a small sphere sinks at, its weight less its buoyancy against the drag.

In [1]:

import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
    root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
from lusine import trait_run
F32 = np.float32
U = float(np.finfo(F32).eps) / 2        # binary32's unit roundoff, 2^-24
rng = np.random.default_rng(2026)

import os, sys root = os.getcwd() while not os.path.isfile(os.path.join(root, 'docsite.json')): root = os.path.dirname(root) sys.path.insert(0, os.path.join(root, 'tools')) import numpy as np from lusine import trait_run F32 = np.float32 U = float(np.finfo(F32).eps) / 2 # binary32's unit roundoff, 2^-24 rng = np.random.default_rng(2026)

The kernel¶

The law as the library states it.

In [2]:

print(trait_run.source('laws/settling.trait'))

print(trait_run.source('laws/settling.trait'))

// lusine.laws.settling: the speed a small particle sinks at through a still fluid. Stokes's law for a sphere at a low
// Reynolds number, its weight less its buoyancy balanced by the drag:
//   w = (rho_p - rho_f) g d^2 / (18 mu),
// downward where the particle is denser than the fluid and upward (negative) where it is lighter. Past a Reynolds
// number rho_f w d / mu of about 1 the drag grows faster and the law overstates the speed; a game's grains of sand or
// silt in water are within it.
// Source: Stokes, Trans. Camb. Phil. Soc. 9, 8 (1851); Batchelor, An Introduction to Fluid Dynamics, section 4.9.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.laws;

float<m/s> settling(float<kg/m3> particle, float<kg/m3> fluid, float<m/s2> gravity, float<m> diameter,
                    float<Pa*s> viscosity)
{
  return (particle - fluid) * gravity * diameter * diameter / (18.0 * viscosity);
}

The program¶

Each section runs the kernel over a lattice of items, one a case: the host sets every parameter, and main() binds the kernel over the items.

In [3]:

print(trait_run.law_program('laws/settling.trait', 'settling'))
from lusine.kernels.laws.settling import settling as reference
def law(**kw):
    return trait_run.law('laws/settling.trait', 'settling', **kw)

print(trait_run.law_program('laws/settling.trait', 'settling')) from lusine.kernels.laws.settling import settling as reference def law(**kw): return trait_run.law('laws/settling.trait', 'settling', **kw)

namespace probe;

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

struct item {
  host float<kg/m3> particle;
  host float<kg/m3> fluid;
  host float<m/s2> gravity;
  host float<m> diameter;
  host float<Pa*s> viscosity;
  float<m/s> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { items.result_0 = lusine.laws.settling(items.particle, items.fluid, items.gravity, items.diameter, items.viscosity); }

Reference: sand and silt in water and air¶

Grains from 2 to 200 micrometres, quartz to clay, in water and in air: the kernel within binary32's roundings of the reference.

In [4]:

d = rng.uniform(2e-6, 2e-4, 400)
p = rng.uniform(1800.0, 2700.0, 400)
fluid = np.where(rng.random(400) < 0.5, 1000.0, 1.2)
mu = np.where(fluid > 2.0, 1.0e-3, 1.8e-5)
got, want = law(particle=p, fluid=fluid, gravity=9.81, diameter=d, viscosity=mu), reference(p, fluid, 9.81, d, mu)
err = np.max(np.abs(got - want) / np.abs(want))
assert err <= 16 * U, err / U
print(f'within {err / U:.1f} x 2^-24 over {len(d)} grains')

d = rng.uniform(2e-6, 2e-4, 400) p = rng.uniform(1800.0, 2700.0, 400) fluid = np.where(rng.random(400) < 0.5, 1000.0, 1.2) mu = np.where(fluid > 2.0, 1.0e-3, 1.8e-5) got, want = law(particle=p, fluid=fluid, gravity=9.81, diameter=d, viscosity=mu), reference(p, fluid, 9.81, d, mu) err = np.max(np.abs(got - want) / np.abs(want)) assert err <= 16 * U, err / U print(f'within {err / U:.1f} x 2^-24 over {len(d)} grains')

within 4.0 x 2^-24 over 400 grains

Invariant: the force balance¶

Stokes's drag on the sinking sphere, 3 pi mu d w, equals its weight less its buoyancy, (rho_p - rho_f) g pi d^3 / 6.

In [5]:

w = law(particle=p, fluid=fluid, gravity=9.81, diameter=d, viscosity=mu)
drag = 3 * np.pi * mu * d * w
net = (p - fluid) * 9.81 * np.pi * d ** 3 / 6
assert np.allclose(drag, net, rtol=1e-5)
print('the drag balances the weight less the buoyancy for every grain')

w = law(particle=p, fluid=fluid, gravity=9.81, diameter=d, viscosity=mu) drag = 3 * np.pi * mu * d * w net = (p - fluid) * 9.81 * np.pi * d ** 3 / 6 assert np.allclose(drag, net, rtol=1e-5) print('the drag balances the weight less the buoyancy for every grain')

the drag balances the weight less the buoyancy for every grain

Limit: neutral and light¶

A particle as dense as the fluid hangs in it; a lighter one rises.

In [6]:

w = law(particle=np.array([1000.0, 900.0]), fluid=1000.0, gravity=9.81, diameter=1e-4, viscosity=1e-3)
assert w[0] == 0.0 and w[1] < 0.0
print('neutral hangs; lighter rises')

w = law(particle=np.array([1000.0, 900.0]), fluid=1000.0, gravity=9.81, diameter=1e-4, viscosity=1e-3) assert w[0] == 0.0 and w[1] < 0.0 print('neutral hangs; lighter rises')

neutral hangs; lighter rises