Quadratic drag¶

The drag of a body's wake, the acceleration it takes per unit mass from the air it moves through: a = -(rho C_d A / 2 m) |u| u, u = v - w its velocity through the air (w the air's own, zero in still air, so a wind is a quantity, not a special case), rho the air's density, C_d the body's drag coefficient and A its frontal area.

The law of a body whose Reynolds number rho |u| D / mu is past about 10^3, where C_d is nearly constant (a smooth sphere's about 0.47), and below its drag crisis (a smooth sphere's near 2 x 10^5), where C_d falls; lusine.properties.reynolds_number says which side a body is on. Below about 1 the drag is Stokes's, lusine.laws.drag. Gravity and any other force are the caller's to add. Written once over a vector of any length. Source: Kundu & Cohen, Fluid Mechanics, 6th ed., 10.11 (drag on a sphere); Taylor, Classical Mechanics, 2.4 (quadratic air resistance). Checked against a reference in numpy, which has its own tests.

Source

Kundu & Cohen, Fluid Mechanics, 6th ed., 10.11; Taylor, Classical Mechanics, 2.4.

Proof

reference, curve, invariant.

quadratic_drag.trait

namespace lusine.laws;

vector<float<m/s2>, _> quadratic_drag(vector<float<m/s>, _> velocity, vector<float<m/s>, _> wind, float<kg/m3> density,
                                     float drag_coefficient, float<m2> area, float<kg> mass)
{
  let u = velocity - wind;
  let speed = sqrt(dot(u, u));
  return u * (speed * -((density * drag_coefficient * area) / (2.0 * mass)));
}

The drag of a body's wake, -(rho C_d A / 2 m) |u| u with u its velocity through the air, the law of a thrown ball (a boule at 7 m/s is at a Reynolds number of about 3 x 10^4, far past Stokes's drag, lusine.laws.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
import matplotlib.pyplot as plt
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 import matplotlib.pyplot as plt 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/quadratic_drag.trait'))

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

// lusine.laws.quadratic_drag: the drag of a body's wake, the acceleration it takes per unit mass from the air it moves
// through: a = -(rho C_d A / 2 m) |u| u, u = v - w its velocity through the air (w the air's own, zero in still air,
// so a wind is a quantity, not a special case), rho the air's density, C_d the body's drag coefficient and A its
// frontal area. The law of a body whose Reynolds number rho |u| D / mu is past about 10^3, where C_d is nearly
// constant (a smooth sphere's about 0.47), and below its drag crisis (a smooth sphere's near 2 x 10^5), where C_d
// falls; lusine.properties.reynolds_number says which side a body is on. Below about 1 the drag is Stokes's,
// lusine.laws.drag. Gravity and any other force are the caller's to add. Written once over a vector of any length.
// Source: Kundu & Cohen, Fluid Mechanics, 6th ed., 10.11 (drag on a sphere); Taylor, Classical Mechanics, 2.4
// (quadratic air resistance).
// Checked against a reference in numpy, which has its own tests.

namespace lusine.laws;

vector<float<m/s2>, _> quadratic_drag(vector<float<m/s>, _> velocity, vector<float<m/s>, _> wind, float<kg/m3> density,
                                     float drag_coefficient, float<m2> area, float<kg> mass)
{
  let u = velocity - wind;
  let speed = sqrt(dot(u, u));
  return u * (speed * -((density * drag_coefficient * area) / (2.0 * mass)));
}

The program¶

Each section below 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/quadratic_drag.trait', 'quadratic_drag', 3))

print(trait_run.law_program('laws/quadratic_drag.trait', 'quadratic_drag', 3))

namespace probe;

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

struct item {
  host vector<float<m/s>, 3> velocity;
  host vector<float<m/s>, 3> wind;
  host float<kg/m3> density;
  host float drag_coefficient;
  host float<m2> area;
  host float<kg> mass;
  vector<float<m/s2>, 3> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { items.result_0 = lusine.laws.quadratic_drag(items.velocity, items.wind, items.density, items.drag_coefficient, items.area, items.mass); }

In [4]:

from lusine.kernels.laws.quadratic_drag import quadratic_drag as reference

from lusine.kernels.laws.quadratic_drag import quadratic_drag as reference

Reference: velocities and winds at random¶

The kernel's acceleration within binary32's roundings of the float64 reference's, for a boule in air (rho 1.204 kg/m^3, C_d 0.47, D 74 mm, 0.7 kg).

In [5]:

rho, cd, r, m = 1.204, 0.47, 0.037, 0.7
area = float(np.pi * r * r)
v = F32(rng.normal(0.0, 8.0, (500, 3))).astype(float)
w = F32(rng.normal(0.0, 3.0, (500, 3))).astype(float)
got = trait_run.law('laws/quadratic_drag.trait', 'quadratic_drag', vector_length=3, velocity=v, wind=w, density=rho,
                    drag_coefficient=cd, area=area, mass=m)
want = reference(v, w, float(F32(rho)), float(F32(cd)), float(F32(area)), float(F32(m)))
err = np.max(np.abs(got - want) / np.maximum(np.abs(want), 1e-6))
assert err <= 64 * U, err
print(f'{len(v)} cases within {err / U:.1f} unit roundoffs of the reference')

rho, cd, r, m = 1.204, 0.47, 0.037, 0.7 area = float(np.pi * r * r) v = F32(rng.normal(0.0, 8.0, (500, 3))).astype(float) w = F32(rng.normal(0.0, 3.0, (500, 3))).astype(float) got = trait_run.law('laws/quadratic_drag.trait', 'quadratic_drag', vector_length=3, velocity=v, wind=w, density=rho, drag_coefficient=cd, area=area, mass=m) want = reference(v, w, float(F32(rho)), float(F32(cd)), float(F32(area)), float(F32(m))) err = np.max(np.abs(got - want) / np.maximum(np.abs(want), 1e-6)) assert err <= 64 * U, err print(f'{len(v)} cases within {err / U:.1f} unit roundoffs of the reference')

500 cases within 4.2 unit roundoffs of the reference

In [6]:

FLIGHT = '''namespace flight;

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

struct body { vector<float<m>, 3> position; vector<float<m/s>, 3> velocity; vector<float<m/s2>, 3> acceleration; };
lattice<body, 1> bodies(spacing = 1 m, boundary = wall);
host struct {
  vector<float<m/s2>, 3> gravity; vector<float<m/s>, 3> wind; vector<float<rad/s>, 3> spin;
  float<kg/m3> density; float drag_coefficient; float<m2> area; float<m> radius; float<kg> mass; float<s> dt;
} settings;

// the body's acceleration: gravity, the wake's drag and the lift of its spin
vector<float<m/s2>, 3> pulled(vector<float<m/s>, 3> velocity, vector<float<m/s2>, 3> gravity, vector<float<m/s>, 3> wind,
                              vector<float<rad/s>, 3> spin, float<kg/m3> density, float drag_coefficient,
                              float<m2> area, float<m> radius, float<kg> mass)
{
  let lift = lusine.laws.lift_coefficient(spin_ratio = lusine.laws.spin_ratio(velocity = velocity, wind = wind,
                                                                               angular_velocity = spin, radius = radius));
  return (gravity + lusine.laws.quadratic_drag(velocity = velocity, wind = wind, density = density,
                                               drag_coefficient = drag_coefficient, area = area, mass = mass))
      + lusine.laws.magnus(velocity = velocity, wind = wind, angular_velocity = spin, density = density,
                           lift_coefficient = lift, area = area, mass = mass);
}

void start()
{
  bodies.acceleration = pulled(bodies.velocity, settings.gravity, settings.wind, settings.spin, settings.density,
                               settings.drag_coefficient, settings.area, settings.radius, settings.mass);
}

void main()
{
  (bodies.position, bodies.velocity) = lusine.integrator.kick_drift(bodies.position, bodies.velocity,
                                                                    bodies.acceleration, settings.dt);
  bodies.acceleration = pulled(bodies.velocity, settings.gravity, settings.wind, settings.spin, settings.density,
                               settings.drag_coefficient, settings.area, settings.radius, settings.mass);
  bodies.velocity = lusine.integrator.kick(bodies.velocity, bodies.acceleration, settings.dt);
}
'''
LIBRARY = ['laws/quadratic_drag.trait', 'laws/magnus.trait', 'laws/spin_ratio.trait', 'laws/lift_coefficient.trait',
           'integrator/kick.trait']
G = [0.0, -9.80665, 0.0]
BOULE = dict(density=1.204, drag_coefficient=0.47, area=float(np.pi * 0.037 ** 2), radius=0.037, mass=0.7)

def fly(velocity, steps, dt=0.001, spin=(0.0, 0.0, 0.0), wind=(0.0, 0.0, 0.0), position=(0.0, 0.0, 0.0), **over):
    run = trait_run.Program(LIBRARY, FLIGHT, {'flight.bodies': [1]})
    settings = dict(BOULE, gravity=G, wind=list(wind), spin=list(spin), dt=dt)
    settings.update(over)
    for name, value in settings.items():
        run.store('flight.settings', name, value)
    run.store('flight.bodies', 'position', [list(position)])
    run.store('flight.bodies', 'velocity', [list(velocity)])
    run.step('flight.start')
    p, v = run.steps('flight.main', steps, [('flight.bodies', 'position'), ('flight.bodies', 'velocity')])
    return p.reshape(steps, 3), v.reshape(steps, 3)

FLIGHT = '''namespace flight; lusine.version version = { { number = 1 } }; struct body { vector<float<m>, 3> position; vector<float<m/s>, 3> velocity; vector<float<m/s2>, 3> acceleration; }; lattice<body, 1> bodies(spacing = 1 m, boundary = wall); host struct { vector<float<m/s2>, 3> gravity; vector<float<m/s>, 3> wind; vector<float<rad/s>, 3> spin; float<kg/m3> density; float drag_coefficient; float<m2> area; float<m> radius; float<kg> mass; float<s> dt; } settings; // the body's acceleration: gravity, the wake's drag and the lift of its spin vector<float<m/s2>, 3> pulled(vector<float<m/s>, 3> velocity, vector<float<m/s2>, 3> gravity, vector<float<m/s>, 3> wind, vector<float<rad/s>, 3> spin, float<kg/m3> density, float drag_coefficient, float<m2> area, float<m> radius, float<kg> mass) { let lift = lusine.laws.lift_coefficient(spin_ratio = lusine.laws.spin_ratio(velocity = velocity, wind = wind, angular_velocity = spin, radius = radius)); return (gravity + lusine.laws.quadratic_drag(velocity = velocity, wind = wind, density = density, drag_coefficient = drag_coefficient, area = area, mass = mass)) + lusine.laws.magnus(velocity = velocity, wind = wind, angular_velocity = spin, density = density, lift_coefficient = lift, area = area, mass = mass); } void start() { bodies.acceleration = pulled(bodies.velocity, settings.gravity, settings.wind, settings.spin, settings.density, settings.drag_coefficient, settings.area, settings.radius, settings.mass); } void main() { (bodies.position, bodies.velocity) = lusine.integrator.kick_drift(bodies.position, bodies.velocity, bodies.acceleration, settings.dt); bodies.acceleration = pulled(bodies.velocity, settings.gravity, settings.wind, settings.spin, settings.density, settings.drag_coefficient, settings.area, settings.radius, settings.mass); bodies.velocity = lusine.integrator.kick(bodies.velocity, bodies.acceleration, settings.dt); } ''' LIBRARY = ['laws/quadratic_drag.trait', 'laws/magnus.trait', 'laws/spin_ratio.trait', 'laws/lift_coefficient.trait', 'integrator/kick.trait'] G = [0.0, -9.80665, 0.0] BOULE = dict(density=1.204, drag_coefficient=0.47, area=float(np.pi * 0.037 ** 2), radius=0.037, mass=0.7) def fly(velocity, steps, dt=0.001, spin=(0.0, 0.0, 0.0), wind=(0.0, 0.0, 0.0), position=(0.0, 0.0, 0.0), **over): run = trait_run.Program(LIBRARY, FLIGHT, {'flight.bodies': [1]}) settings = dict(BOULE, gravity=G, wind=list(wind), spin=list(spin), dt=dt) settings.update(over) for name, value in settings.items(): run.store('flight.settings', name, value) run.store('flight.bodies', 'position', [list(position)]) run.store('flight.bodies', 'velocity', [list(velocity)]) run.step('flight.start') p, v = run.steps('flight.main', steps, [('flight.bodies', 'position'), ('flight.bodies', 'velocity')]) return p.reshape(steps, 3), v.reshape(steps, 3)

Curve: the fall to the terminal velocity¶

A boule dropped from rest, gravity and the wake's drag stepped by the library's velocity Verlet: its speed reaches sqrt(2 m g / (rho C_d A)) and stays there.

In [7]:

vt = np.sqrt(2 * 0.7 * 9.80665 / (1.204 * 0.47 * np.pi * 0.037 ** 2))
p, v = fly([0.0, 0.0, 0.0], 30000, dt=0.002)
ts = 0.002 * np.arange(1, 30001)
speed = -v[:, 1]
want = vt * np.tanh(9.80665 * ts / vt)          # the quadratic law's fall from rest, in closed form
err = np.max(np.abs(speed - want))
assert abs(speed[-1] - vt) < 1e-3 * vt and err < 2e-3 * vt, (speed[-1], vt, err)
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(ts, speed, label='stepped by the emitted code')
ax.plot(ts, want, '--', label='v_t tanh(g t / v_t)')
ax.set_xlabel('t (s)'); ax.set_ylabel('speed (m/s)'); ax.legend()
plt.show()
print(f'terminal {speed[-1]:.2f} m/s against sqrt(2 m g / (rho C_d A)) = {vt:.2f}; worst {err:.3f} m/s off the closed form')

vt = np.sqrt(2 * 0.7 * 9.80665 / (1.204 * 0.47 * np.pi * 0.037 ** 2)) p, v = fly([0.0, 0.0, 0.0], 30000, dt=0.002) ts = 0.002 * np.arange(1, 30001) speed = -v[:, 1] want = vt * np.tanh(9.80665 * ts / vt) # the quadratic law's fall from rest, in closed form err = np.max(np.abs(speed - want)) assert abs(speed[-1] - vt) < 1e-3 * vt and err < 2e-3 * vt, (speed[-1], vt, err) fig, ax = plt.subplots(figsize=(7, 3)) ax.plot(ts, speed, label='stepped by the emitted code') ax.plot(ts, want, '--', label='v_t tanh(g t / v_t)') ax.set_xlabel('t (s)'); ax.set_ylabel('speed (m/s)'); ax.legend() plt.show() print(f'terminal {speed[-1]:.2f} m/s against sqrt(2 m g / (rho C_d A)) = {vt:.2f}; worst {err:.3f} m/s off the closed form')

terminal 75.09 m/s against sqrt(2 m g / (rho C_d A)) = 75.11; worst 0.015 m/s off the closed form

Curve: a coast¶

With no gravity a body slows as v(t) = v0 / (1 + k v0 t), k = rho C_d A / 2 m: a boule launched at 7 m/s.

In [8]:

k = 1.204 * 0.47 * np.pi * 0.037 ** 2 / (2 * 0.7)
p, v = fly([7.0, 0.0, 0.0], 5000, dt=0.002, gravity=[0.0, 0.0, 0.0])
ts = 0.002 * np.arange(1, 5001)
want = 7.0 / (1.0 + k * 7.0 * ts)
err = np.max(np.abs(v[:, 0] - want))
assert err < 1e-4, err
print(f'after {ts[-1]:.0f} s: {v[-1, 0]:.4f} m/s, the closed form {want[-1]:.4f}; worst {err:.2g} m/s')

k = 1.204 * 0.47 * np.pi * 0.037 ** 2 / (2 * 0.7) p, v = fly([7.0, 0.0, 0.0], 5000, dt=0.002, gravity=[0.0, 0.0, 0.0]) ts = 0.002 * np.arange(1, 5001) want = 7.0 / (1.0 + k * 7.0 * ts) err = np.max(np.abs(v[:, 0] - want)) assert err < 1e-4, err print(f'after {ts[-1]:.0f} s: {v[-1, 0]:.4f} m/s, the closed form {want[-1]:.4f}; worst {err:.2g} m/s')

after 10 s: 6.2406 m/s, the closed form 6.2406; worst 2.3e-05 m/s

Invariant: the air's frame¶

A body moving with the wind feels no drag; a wind on a body at rest pushes it as moving through still air holds it back.

In [9]:

def drag(v, w):
    return trait_run.law('laws/quadratic_drag.trait', 'quadratic_drag', vector_length=3, velocity=v, wind=w,
                         density=1.204, drag_coefficient=0.47, area=float(np.pi * 0.037 ** 2), mass=0.7)
v = np.array([[3.0, -1.0, 2.0]])
assert np.all(drag(v, v) == 0.0)
assert np.array_equal(drag(np.zeros((1, 3)), v), -drag(v, np.zeros((1, 3))))
print('with the wind, none; a wind is a motion through the air')

def drag(v, w): return trait_run.law('laws/quadratic_drag.trait', 'quadratic_drag', vector_length=3, velocity=v, wind=w, density=1.204, drag_coefficient=0.47, area=float(np.pi * 0.037 ** 2), mass=0.7) v = np.array([[3.0, -1.0, 2.0]]) assert np.all(drag(v, v) == 0.0) assert np.array_equal(drag(np.zeros((1, 3)), v), -drag(v, np.zeros((1, 3)))) print('with the wind, none; a wind is a motion through the air')

with the wind, none; a wind is a motion through the air