Motion with drag and variable mass

Introductory problems usually ignore air resistance and assume that the mass of a moving body never changes. Real projectiles feel drag, and the drag force depends on speed in a way that makes the equations nonlinear. Rockets throw away most of their mass as they fly, so Newton's second law in the form \(F = ma\) does not even apply directly. This lesson solves both kinds of problem, exactly where possible and numerically otherwise. The results explain why a baseball travels about half as far as the vacuum formula predicts, why the best launch angle is less than \(45^\circ\), and why every orbital rocket is built in stages.

Linear and quadratic drag

A body moving through a fluid feels a drag force opposite to its velocity. Two regimes matter. At very low speeds, or for very small bodies, the fluid flows smoothly around the object and the drag is proportional to speed. For a sphere of radius \(R\) in a fluid of viscosity \(\eta\), George Stokes found in 1851 that \(F = 6\pi\eta Rv\). At higher speeds the fluid breaks into a turbulent wake, and the drag is proportional to the square of the speed,

where \(\rho\) is the fluid density, \(A\) the frontal area and \(C_D\) a drag coefficient of order one (about 0.47 for a smooth sphere over a wide range of speeds, about 0.35 for a baseball, about 1 for an upright person). The dividing line is the Reynolds number \(\mathrm{Re} = \rho vD/\eta\), the ratio of the two drag mechanisms: linear drag dominates for \(\mathrm{Re} \lesssim 1\) and quadratic drag for \(\mathrm{Re} \gtrsim 1000\). For a baseball at 40 m/s in air (\(\eta = 1.8\times10^{-5}\ \mathrm{Pa\,s}\)), \(\mathrm{Re} \approx 2\times10^5\), firmly quadratic. For a 10 μm droplet of fog settling at 3 mm/s, \(\mathrm{Re} \approx 0.002\), firmly linear.

Linear drag: an exact solution

With linear drag \(\mathbf F = -b\mathbf v\), the equation \(m\dot{\mathbf v} = m\mathbf g - b\mathbf v\) is linear and its components separate. Define the time constant \(\tau = m/b\). The horizontal velocity decays exponentially, \(v_x = v_{x0}e^{-t/\tau}\), so the horizontal distance approaches the finite limit \(v_{x0}\tau\) however long the flight lasts. The vertical equation \(\dot v_y = -g - v_y/\tau\) has the solution

The velocity relaxes exponentially to the terminal velocity \(-v_t\), the speed at which drag balances weight. Integrating once more gives the trajectory. A useful check: expanding \(e^{-t/\tau} \approx 1 - t/\tau + t^2/2\tau^2\) for \(t \ll \tau\) recovers \(y = v_{y0}t - \tfrac12 gt^2\), the vacuum result. Linear drag describes a sediment grain settling in a lake, a droplet in a cloud, and Robert Millikan's oil drops of 1909.

Quadratic drag

With quadratic drag \(\mathbf F = -mk|\mathbf v|\mathbf v\) (where \(k = C_D\rho A/2m\)), the components are coupled through \(|\mathbf v| = \sqrt{v_x^2 + v_y^2}\) and there is no closed-form solution for a general trajectory. Purely vertical motion can still be solved. For a body falling from rest, with \(y\) measured downward, \(\dot v = g - kv^2\). The terminal speed is \(v_t = \sqrt{g/k}\), and separating variables,

The speed reaches 76% of terminal speed after one time \(v_t/g\) and 96% after two. Writing \(v\,dv/dy = g - kv^2\) and integrating gives the speed as a function of distance fallen, \(v^2 = v_t^2(1 - e^{-2ky})\), so the characteristic distance is \(1/2k\). A skydiver in a spread position, with \(v_t \approx 55\ \mathrm{m/s}\), reaches about 96% of terminal speed after falling roughly 400 m.

Projectiles with drag

For a projectile with quadratic drag we integrate numerically. The program below follows a baseball launched at 45 m/s (about 100 mph), using RK4 from lesson 1, and searches for the launch angle that gives the longest range. It ignores the ball's spin, which in a real hit adds lift through the Magnus effect and can extend the range by 10 to 20 m.

import math

def rk4(f, y, dt):
    k1 = f(y); k2 = f([a + dt / 2 * b for a, b in zip(y, k1)])
    k3 = f([a + dt / 2 * b for a, b in zip(y, k2)]); k4 = f([a + dt * b for a, b in zip(y, k3)])
    return [a + dt / 6 * (b1 + 2 * b2 + 2 * b3 + b4) for a, b1, b2, b3, b4 in zip(y, k1, k2, k3, k4)]

# A baseball (m = 0.145 kg, d = 74 mm) with quadratic drag F = (1/2) rho Cd A v^2
g, rho, Cd, m, d = 9.81, 1.2, 0.35, 0.145, 0.074
k = 0.5 * rho * Cd * math.pi * d * d / 4 / m          # drag constant per unit mass (1/m)
vt = math.sqrt(g / k)
print(f"terminal speed = {vt:.1f} m/s")

def flight(v0, deg, drag=True):
    th = math.radians(deg)
    y = [0.0, 0.0, v0 * math.cos(th), v0 * math.sin(th)]
    def f(s):
        v = math.hypot(s[2], s[3]); c = k * v if drag else 0.0
        return [s[2], s[3], -c * s[2], -g - c * s[3]]
    dt = 0.001
    while True:
        new = rk4(f, y, dt)
        if new[1] < 0:                       # interpolate the landing point
            frac = y[1] / (y[1] - new[1])
            return y[0] + frac * (new[0] - y[0])
        y = new

v0 = 45.0
print(f"vacuum range at 45 deg: {flight(v0, 45, drag=False):.1f} m (formula {v0 ** 2 / g:.1f} m)")
best = max(range(20, 56), key=lambda a: flight(v0, a))
for a in (30, 35, best, 45):
    print(f"with drag at {a} deg: range {flight(v0, a):.1f} m")
print(f"best launch angle with drag: {best} deg")

# Rocket equation: dv = ve ln(m0/mf). One stage versus two stages carrying the same propellant.
ve, payload = 3000.0, 1000.0
def stage_dv(m0, mprop):
    return ve * math.log(m0 / (m0 - mprop))
# single stage: 90 t of propellant, 10 t of structure
dv1 = stage_dv(payload + 100e3, 90e3)
# two stages: each carries 45 t of propellant with 5 t of structure; the empty first stage is dropped
m_upper = payload + 50e3
dv2 = stage_dv(m_upper + 50e3, 45e3) + stage_dv(m_upper, 45e3)
print(f"single stage delta-v = {dv1 / 1000:.2f} km/s, two stages = {dv2 / 1000:.2f} km/s")
print(f"mass ratio needed for 9.4 km/s to orbit in one stage: {math.exp(9400 / ve):.1f}")

terminal speed = 39.7 m/s
vacuum range at 45 deg: 206.4 m (formula 206.4 m)
with drag at 30 deg: range 105.5 m
with drag at 35 deg: range 109.7 m
with drag at 40 deg: range 111.1 m
with drag at 45 deg: range 110.0 m
best launch angle with drag: 40 deg
single stage delta-v = 6.65 km/s, two stages = 8.19 km/s
mass ratio needed for 9.4 km/s to orbit in one stage: 23.0

In a vacuum the ball would travel 206 m. With drag it travels only about 110 m at \(45^\circ\), and the best angle drops to about \(40^\circ\). The reason is that drag robs the ball of speed most heavily when it is moving fastest, early in the flight. A lower trajectory spends less time in the air, so it loses less speed, and the horizontal component is a larger share of what remains. For a golf ball or a javelin the optimum is lower still. The terminal speed of 39.7 m/s is also revealing: the ball is launched faster than terminal speed, so on the way down it actually slows toward \(v_t\), and the descending branch of the path is steeper than the ascending one.

Bodies with changing mass

Newton wrote his second law as \(\mathbf F = d\mathbf p/dt\), but applying it to a body whose mass changes requires care, because the mass that leaves or arrives carries momentum of its own. The safe method is to apply momentum conservation to a closed system. Consider a rocket of mass \(m\) moving at velocity \(v\). In a short time \(dt\) it ejects mass \(|dm|\) of exhaust at speed \(v_e\) relative to the rocket, so the exhaust moves at \(v - v_e\). With an external force \(F\) (gravity or drag), the momentum balance between \(t\) and \(t + dt\) is

Expanding and dropping the second-order term \(|dm|\,dv\) gives \(m\,dv - v_e|dm| = F\,dt\). Writing \(|dm| = -dm\), since the rocket's mass decreases,

The term \(-v_e\,dm/dt\) is the thrust: exhaust mass flow rate times exhaust speed. With no external force, separating variables gives Konstantin Tsiolkovsky's rocket equation of 1903,

where \(m_0\) is the starting mass and \(m_f\) the mass when the propellant is gone. The velocity gained depends only on the exhaust speed and the mass ratio, not on how fast the propellant is burned. The logarithm is the tyranny of rocketry: doubling the mass ratio adds only \(v_e\ln 2 \approx 0.69v_e\). Chemical engines reach exhaust speeds of about 3 km/s (kerosene and oxygen) to 4.4 km/s (hydrogen and oxygen in vacuum), while reaching low Earth orbit requires a velocity change of about 9.4 km/s once gravity and drag losses are counted.

A mass ratio of 23 means the rocket must be about 96% propellant by mass at liftoff, leaving 4% for tanks, engines and payload. Nobody can build tanks that light, which is why the program finds that a single stage carrying 90 t of propellant and 10 t of structure reaches only 6.65 km/s. Staging solves the problem. Split the same propellant and structure into two stages, and drop the first stage's empty tanks before lighting the second. The second stage no longer has to accelerate dead weight, and the same hardware reaches 8.19 km/s. Every orbital launcher since Sputnik in 1957 has used at least two stages; reusable boosters recover the first stage but still drop it.

Gravity loss

Near the ground, gravity subtracts \(g\) from the rocket's acceleration for every second it spends thrusting upward. Integrating \(\dot v = -v_e\dot m/m - g\) for a vertical burn of duration \(t_b\) gives \(\Delta v = v_e\ln(m_0/m_f) - gt_b\). The penalty \(gt_b\) is the gravity loss, typically 1 to 1.5 km/s for a real launch. It favors high thrust and short burns. Real rockets turn over toward horizontal soon after clearing the dense lower atmosphere, so that much of the thrust builds orbital speed rather than fighting gravity, and drag loss adds a further 0.1 to 0.3 km/s.

Other variable-mass systems

The same momentum bookkeeping applies whenever mass is gained or lost. A raindrop falling through a cloud sweeps up droplets; a conveyor belt receives sand from a hopper; a chain falls through a hole in a table. In each case write the momentum of the whole system, including the mass being added or removed, at two nearby times. Attempting to use \(F = d(mv)/dt\) for the body alone, without accounting for the velocity of the incoming or outgoing mass, gives wrong answers whenever that velocity differs from zero in your frame.