Quadratic drag, live¶
The drag of a body's wake, the acceleration it takes per unit mass from the air it moves through: a = −(ρ C_d A / 2m) |u| u, u = v − w its velocity through the air, ρ the air's density, C_d the body's drag coefficient and A its frontal area. It is the law of a thrown ball, far past the slow flow where drag is linear. Below, two boules (a smooth sphere's C_d of 0.47, 3.7 cm in radius, 0.7 kg) are stepped by the library's velocity Verlet (lusine.integrator.kick_drift, the drag lusine.laws.quadratic_drag, lusine.integrator.kick), a hundredth of a second a step and two hundred steps a second, on your GPU as you read: one dropped from rest under gravity, the other thrown along at 7 m/s with no gravity at all. Beside them the same program computes the closed forms, the fall v_t tanh(g t / v_t) toward the terminal speed v_t = √(2 m g / (ρ C_d A)), and the coast v₀ / (1 + k v₀ t) with k = ρ C_d A / 2m, and checks each against its boule at every step.
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 fall: the chart draws the falling boule's speed and the closed form against its time, a point a step: the speed rising to the terminal speed, about 75 m/s, and staying there. Every thirty seconds both boules start again. Below it, the time, the speed and the terminal speed.
- The coast: the thrown boule's speed and its closed form, slowing from 7 m/s with nothing to keep it going.
- The checks: the chart draws each boule's miss of its closed form at each step as a share of what its check allows, on a log scale, so a check holds while its line stays below 1. falls allows 0.02% of the terminal speed, coasts 0.001% of the throw. the air's frame holds while a ball moving with a wind of (3, −1, 2) m/s has felt no drag at all, and that wind has pushed a ball at rest exactly as still air holds back a ball moving as fast.
Every number here is computed by the program on the GPU; the page only reads the values back and draws them. Gravity, the air's density, the boule's coefficient, area and mass, the throw and the step are the program's settings, which the host may write; this page offers no control to change them yet.
The program¶
The boules, the closed forms and the checks, the panels the page draws, and the one procedure, main, run each step: one velocity Verlet step for each boule, the drag taken at the step's start and at its middle, and the checks after it.
// Quadratic drag's live page: a boule dropped from rest, and another thrown along with no gravity, each stepped by the
// library's velocity Verlet under the drag of its wake, beside the closed forms of the fall and the coast, with the
// page's checks computed here on the GPU and shown by the panels the program declares.
namespace quadratic_drag;
lusine.version version = { { number = 1 } };
// The reader's settings: gravity, the air's density, the boule's drag coefficient, frontal area (3.7 cm in radius) and
// mass, the coasting boule's speed at the start, and the step.
host struct {
float<m/s2> gravity = 9.80665 m/s2; float<kg/m3> density = 1.204 kg/m3; float drag_coefficient = 0.47;
float<m2> area = 0.0043008403 m2; float<kg> mass = 0.7 kg; float<m/s> throw = 7.0 m/s; float<s> dt = 0.01 s;
} settings;
// The falling boule, its velocity up and its height, and the coasting one, its velocity along and how far it has
// gone; each one's acceleration; their time the sum of the steps they have taken.
struct {
float<m> height; float<m/s> fall; float<m/s2> fall_acceleration; float<m> distance; float<m/s> coast;
float<m/s2> coast_acceleration; float<s> time;
} boules;
// The closed forms at their time: the fall from rest, v_t tanh(g t / v_t) downward, v_t = sqrt(2 m g / (rho C_d A))
// the terminal speed; the coast, v0 / (1 + k v0 t), k = rho C_d A / (2 m); and the drag in the air's frame: on a ball
// moving with a wind, and on a ball at rest in that wind and on one moving through still air as fast.
struct {
float<m/s> terminal; float<m/s> speed; float<m/s> fall; float<m/s> coast; float fall_share; float coast_share;
vector<float<m/s2>, 3> with_the_wind; vector<float<m/s2>, 3> in_the_wind; vector<float<m/s2>, 3> through_still_air;
} closed;
// The checks: the fall within 0.02% of the terminal speed of its closed form at every step; the coast within 0.001% of
// the throw of its closed form at every step; and the air's frame, no drag moving with the wind, and a wind pushing a
// ball at rest exactly as still air holds back a ball moving as fast.
struct { bool falls = true; bool coasts = true; bool air_frame = 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<m/s> terminal_speed(float<kg> mass, float<m/s2> gravity, float<kg/m3> density, float drag_coefficient,
float<m2> area)
{
return sqrt(2.0 * mass * gravity / (density * drag_coefficient * area));
}
float<m/s> falling(float<m/s> terminal, float<m/s2> gravity, float<s> time)
{
let x = gravity * time / terminal;
let e = exp(-2.0 * x);
return -terminal * (1.0 - e) / (1.0 + e);
}
float<m/s> coasting(float<m/s> throw, float<kg/m3> density, float drag_coefficient, float<m2> area, float<kg> mass,
float<s> time)
{
return throw / (1.0 + density * drag_coefficient * area / (2.0 * mass) * throw * time);
}
vector<float<m/s>, 3> breeze() { return vector<float<m/s>, 3>(3.0 m/s, -1.0 m/s, 2.0 m/s); }
vector<float<m/s>, 3> still() { return vector<float<m/s>, 3>(0.0 m/s, 0.0 m/s, 0.0 m/s); }
bool nothing(vector<float<m/s2>, 3> a) { return a == vector<float<m/s2>, 3>(0.0 m/s2, 0.0 m/s2, 0.0 m/s2); }
bool reversed(vector<float<m/s2>, 3> a, vector<float<m/s2>, 3> b) { return a == -b; }
output table<lusine.widget> panels = {
lusine.series { { panel = "the fall", chart = "speed", quantity = closed.speed, steps = 3000, against = boules.time,
label = "stepped" } },
lusine.series { { panel = "the fall", chart = "speed", quantity = closed.fall, steps = 3000, against = boules.time,
label = "v_t tanh(g t / v_t)" } },
lusine.value { { panel = "the fall", quantity = boules.time, label = "time" } },
lusine.value { { panel = "the fall", quantity = closed.speed, label = "stepped speed" } },
lusine.value { { panel = "the fall", quantity = closed.terminal, label = "terminal speed" } },
lusine.series { { panel = "the coast", chart = "speed along", quantity = boules.coast, steps = 3000,
against = boules.time, label = "stepped" } },
lusine.series { { panel = "the coast", chart = "speed along", quantity = closed.coast, steps = 3000,
against = boules.time, label = "v0 / (1 + k v0 t)" } },
lusine.value { { panel = "the coast", quantity = boules.coast, label = "stepped speed" } },
lusine.value { { panel = "the coast", quantity = closed.coast, label = "closed form" } },
lusine.series { { panel = "the checks", chart = "share of the tolerance", quantity = closed.fall_share, steps = 3000,
log = true, label = "fall" } },
lusine.series { { panel = "the checks", chart = "share of the tolerance", quantity = closed.coast_share,
steps = 3000, log = true, label = "coast" } },
lusine.value { { panel = "the checks", quantity = checks.falls, label = "falls" } },
lusine.value { { panel = "the checks", quantity = checks.coasts, label = "coasts" } },
lusine.value { { panel = "the checks", quantity = checks.air_frame, label = "the air's frame" } },
};
void main()
{
// Every thirty seconds both boules start again: the one dropped from rest, the other thrown.
boules.height = boules.time >= 30.0 s ? 0.0 m : boules.height;
boules.fall = boules.time >= 30.0 s ? 0.0 m/s : boules.fall;
boules.distance = boules.time >= 30.0 s ? 0.0 m : boules.distance;
boules.coast = boules.time == 0.0 s || boules.time >= 30.0 s ? settings.throw : boules.coast;
boules.time = boules.time >= 30.0 s ? 0.0 s : boules.time;
boules.fall_acceleration = lusine.laws.quadratic_drag(velocity = boules.fall, wind = 0.0 m/s, settings.density,
settings.drag_coefficient, settings.area, settings.mass)
- settings.gravity;
boules.coast_acceleration = lusine.laws.quadratic_drag(velocity = boules.coast, wind = 0.0 m/s, settings.density,
settings.drag_coefficient, settings.area, settings.mass);
(boules.height, boules.fall) = lusine.integrator.kick_drift(position = boules.height, velocity = boules.fall,
acceleration = boules.fall_acceleration, settings.dt);
(boules.distance, boules.coast) = lusine.integrator.kick_drift(position = boules.distance, velocity = boules.coast,
acceleration = boules.coast_acceleration,
settings.dt);
boules.fall_acceleration = lusine.laws.quadratic_drag(velocity = boules.fall, wind = 0.0 m/s, settings.density,
settings.drag_coefficient, settings.area, settings.mass)
- settings.gravity;
boules.coast_acceleration = lusine.laws.quadratic_drag(velocity = boules.coast, wind = 0.0 m/s, settings.density,
settings.drag_coefficient, settings.area, settings.mass);
boules.fall = lusine.integrator.kick(velocity = boules.fall, acceleration = boules.fall_acceleration, settings.dt);
boules.coast = lusine.integrator.kick(velocity = boules.coast, acceleration = boules.coast_acceleration,
settings.dt);
boules.time = boules.time + settings.dt;
closed.terminal = terminal_speed(settings.mass, settings.gravity, settings.density, settings.drag_coefficient,
settings.area);
closed.speed = -boules.fall;
closed.fall = -falling(closed.terminal, settings.gravity, boules.time);
closed.coast = coasting(settings.throw, settings.density, settings.drag_coefficient, settings.area, settings.mass,
boules.time);
closed.fall_share = abs(closed.speed - closed.fall) / (0.0002 * closed.terminal);
closed.coast_share = abs(boules.coast - closed.coast) / (0.00001 * settings.throw);
checks.falls = checks.falls && closed.fall_share <= 1.0;
checks.coasts = checks.coasts && closed.coast_share <= 1.0;
closed.with_the_wind = lusine.laws.quadratic_drag(velocity = breeze(), wind = breeze(), settings.density,
settings.drag_coefficient, settings.area, settings.mass);
closed.in_the_wind = lusine.laws.quadratic_drag(velocity = still(), wind = breeze(), settings.density,
settings.drag_coefficient, settings.area, settings.mass);
closed.through_still_air = lusine.laws.quadratic_drag(velocity = breeze(), wind = still(), settings.density,
settings.drag_coefficient, settings.area, settings.mass);
checks.air_frame = checks.air_frame && nothing(a = closed.with_the_wind)
&& reversed(a = closed.in_the_wind, b = closed.through_still_air);
lamps.colour = lamp_colour(holds = checks.falls && checks.coasts && checks.air_frame);
}
The law's other proof (the drag against a reference over velocities and winds at random) is on its proof page.