Adding a v Term to a v² Term Got Both Ends Right and Missed the Middle by 42%

Summing linear and quadratic drag works for Re ≪ 60 and Re ≫ 60, but for Re between 10 and 300 it leaves out the wake drag and overpredicts terminal velocity by 30–40%.

Same water, same glass, different drag law#

Drop two glass beads into water. One is 0.05 mm across, the other 10 mm. The small one feels drag proportional to its speed. The large one feels drag proportional to the square of its speed. Same water, same material. Only the size differs.

French scientists in the late 18th century found these two laws separately. Then they merged them into one formula, $F = a v^2 + b v$ . This post checks, with numbers, where that sum holds and where it fails. You will come away with the Reynolds number (the ratio of inertial to viscous forces) at which the two terms trade places, and the range where the simple sum misses the most.

We walk through it as five questions.

Q1. Between 1775 and 1800, was drag v or v²?#

In 1752 d'Alembert hit his paradox: an ideal fluid exerts zero drag on a body. With theory stuck, experiments took over. In 1775 the French Academy set up a team to study fluid resistance. Bossut ran the work.

Bossut towed model ships through a tank. Careful measurements showed drag growing with the square of speed. Meanwhile Du Buat measured flow inside pipes. The friction he saw grew simply in proportion to speed. Two experiments, two answers.

Prony ran more precise experiments and tied them together. Resistance could be written as $a v^2 + b v$ . Coulomb organized this into two kinds of drag. The part growing with $v^2$ is form drag (the pressure difference set by body shape). The part growing with $v$ is friction drag. That split later fed into Navier's and Stokes's work on viscous flow.

So who was right? Both. The model ships were large and fast. The pipe flows were slow and narrow. The two experiments saw opposite ends of the same physics. The number that separates those ends is the next question.

Q2. Where do the two terms come from?#

Written as two terms, the drag on a sphere is

F = \underbrace{3\pi\mu d\, v}_{b\,v} + \underbrace{C_\infty \frac{\pi}{8}\rho d^2 v^2}_{a\,v^2}

$\mu$ is viscosity, $\rho$ fluid density, $d$ sphere diameter, and $C_\infty$ the high-speed drag coefficient (about 0.4 for a sphere).

The first term is Stokes drag. It is the exact solution for slow flow with no inertia at all (creeping flow). The second is Newton drag. It follows if the momentum of the fluid hitting the front face turns straight into pressure.

One caveat. Stokes drag $3\pi\mu d v$ is not pure skin friction. Two-thirds of it is friction and one-third is pressure. So "viscous-dominated vs inertia-dominated" is a more accurate label than Coulomb's "friction vs form."

Try it yourself in the simulation below. Drop a bead and watch both drag forces grow as it speeds up. Where their sum equals the weight is the terminal velocity (the speed at which acceleration stops).

d = 0.30 mm fluid:

v = 0.00 mm/s · v_t = 1000 mm/s · Re = 0.00 · linear share 100.0% · quadratic 0.0%

Sweep the diameter slider from 0.02 mm to 20 mm. For small beads the blue line (b·v) meets the weight line; for large ones the orange parabola (a·v²) does. Switch the fluid to glycerol and even a 10 mm bead moves over to Du Buat's side.

Q3. When does control change hands?#

Find the speed at which the two terms are equal.

3\pi\mu d\, v = C_\infty \frac{\pi}{8}\rho d^2 v^2 \quad\Longrightarrow\quad Re \equiv \frac{\rho v d}{\mu} = \frac{24}{C_\infty}

$Re$ is the Reynolds number based on sphere diameter. With $C_\infty = 0.4$ , $Re = 60$ .

This result reconciles Bossut and Du Buat. The drag law is set neither by the kind of fluid nor by speed alone. It is set by $Re$ . Slow pipe flow had small $Re$ . Model ships ran at $Re$ in the hundreds of thousands. Swapping water for glycerol in the simulation above was just moving $Re$ .

The same "bridge two limits" problem shows up elsewhere. For particle heating in the Knudsen transition regime, the continuum and free-molecular limits were joined by a harmonic mean. For Wood's sound speed, compressibilities were added by volume fraction. How you blend the limits changes the answer. Prony chose plain addition.

Q4. So is adding them enough?#

Let's check. Release six glass beads ( $\rho_p = 2500$ kg/m³) from rest in water. Integrate the equation of motion with two-term drag.

m \frac{dv}{dt} = (\rho_p - \rho) V g - b v - a v^2

$m$ is the sphere mass, $V$ its volume, and $g$ gravity. Added mass (the inertia of fluid dragged along with the sphere) is ignored.

Compare terminal velocity against a reference: the Schiller–Naumann correlation $C_D = \frac{24}{Re}(1 + 0.15 Re^{0.687})$ . It is known to match experiments within a few percent for $Re < 1000$ .

import math
 
# Glass beads settling in water: Du Buat (linear) vs Bossut (quadratic) drag
RHO_P, RHO, MU, G = 2500.0, 1000.0, 1.0e-3, 9.81
C_INF = 0.4          # quadratic (form) drag coefficient of the two-term law
 
def weight(d):
    return (RHO_P - RHO) * math.pi * d**3 / 6 * G
 
def prony_drag(d, v):
    lin = 3 * math.pi * MU * d * v                 # Du Buat / Stokes:  b v
    quad = C_INF * math.pi / 8 * RHO * d**2 * v**2  # Bossut / Newton:  a v^2
    return lin, quad
 
def vt_prony(d):
    a = C_INF * math.pi / 8 * RHO * d**2
    b = 3 * math.pi * MU * d
    return (-b + math.sqrt(b * b + 4 * a * weight(d))) / (2 * a)
 
def cd_schiller(re):
    return 24 / re * (1 + 0.15 * re**0.687) if re < 1000 else 0.44
 
def vt_reference(d):
    lo, hi = 1e-9, 10.0
    for _ in range(200):
        v = 0.5 * (lo + hi)
        re = RHO * v * d / MU
        f = cd_schiller(re) * 0.5 * RHO * v * v * math.pi * d * d / 4
        lo, hi = (v, hi) if f < weight(d) else (lo, v)
    return 0.5 * (lo + hi)
 
def settle(d, steps_per_tau=400):
    m = RHO_P * math.pi * d**3 / 6
    vt = vt_prony(d)
    tau = vt / ((1 - RHO / RHO_P) * G)       # time scale of acceleration
    dt = tau / steps_per_tau
    v, t = 0.0, 0.0
    while v < 0.99 * vt:
        lin, quad = prony_drag(d, v)
        v += dt * (weight(d) - lin - quad) / m
        t += dt
    return t, tau, vt
 
print(" d[mm]   Re_t    vt[mm/s]  ref[mm/s]  err%   lin%   t99/tau  Stokes/vt  Newton/vt")
for dmm in [0.05, 0.1, 0.3, 1.0, 3.0, 10.0]:
    d = dmm * 1e-3
    t99, tau, vt = settle(d)
    lin, quad = prony_drag(d, vt)
    ref = vt_reference(d)
    v_st = weight(d) / (3 * math.pi * MU * d)
    v_nw = math.sqrt(weight(d) / (C_INF * math.pi / 8 * RHO * d * d))
    re = RHO * vt * d / MU
    print(f"{dmm:6.2f} {re:7.2f} {vt*1e3:10.2f} {ref*1e3:10.2f} {100*(vt-ref)/ref:6.1f} "
          f"{100*lin/(lin+quad):6.1f} {t99/tau:8.2f} {v_st/vt:10.2f} {v_nw/vt:10.2f}")
re_x = 24 / C_INF  # 3 pi mu d v = C_inf pi/8 rho d^2 v^2
print(f"crossover Re where lin = quad: 24/C_inf = {re_x:.1f}")

 d[mm]   Re_t    vt[mm/s]  ref[mm/s]  err%   lin%   t99/tau  Stokes/vt  Newton/vt
  0.05    0.10       2.04       1.98    2.9   99.8     4.59       1.00      24.27
  0.10    0.81       8.07       7.29   10.6   98.7     4.55       1.01       8.68
  0.30   17.16      57.21      40.21   42.3   77.8     3.92       1.29       2.12
  1.00  193.49     193.49     145.95   32.6   23.7     2.93       4.22       1.14
  3.00 1121.20     373.73     365.75    2.2    5.1     2.70      19.69       1.03
 10.00 6973.63     697.36     667.76    4.4    0.9     2.65     117.23       1.00
crossover Re where lin = quad: 24/C_inf = 60.0

Now run that code in your browser. Two identical beads sink side by side. One feels Prony's two-term drag, the other the Schiller–Naumann reference.

d = 0.30 mm C∞ = 0.40 ρp = 2500 kg/m³ steps/frame 8

step 0 · Prony v = 0.00 · ref v = 0.00 mm/s · t99/τ = …
Re_t 17.16 · vt 57.21 · ref 40.21 · err 42.3% · lin 77.8% · Stokes/vt 1.29 · Newton/vt 2.12

At the default d = 0.30 mm you get err 42.3%, lin 77.8%, and t99/τ 3.92, matching the output block. The right-hand curve is the error across all diameters. Raise C∞ to 0.6 and the large-bead error turns negative; lower ρp and Re_t drops at a given diameter, pushing the peak to the right.

The output says three things. First, both ends fit well: 2.9% at 0.05 mm, 2.2% at 3 mm. Second, using only one term collapses at the opposite end. Stokes alone makes the 10 mm bead 117 times too fast. Newton alone makes the 0.05 mm bead 24 times too fast. Third, the middle is the problem: 42% off at 0.3 mm ( $Re_t = 17$ ) and 33% off at 1 mm ( $Re_t = 193$ ).

Acceleration time follows the law too. With linear drag only, speed approaches its limit exponentially, so $t_{99}/\tau = \ln 100 = 4.61$ . With quadratic drag only, it approaches as $\tanh$ , so $t_{99}/\tau = \operatorname{artanh} 0.99 = 2.65$ . The 4.59 and 2.65 in the output are those two limits; intermediate diameters pass between them.

Q5. Why does it miss in the middle?#

Rewrite the two-term law as a drag coefficient and the cause appears.

C_D^{\text{Prony}} = \frac{24}{Re} + C_\infty, \qquad C_D^{\text{S–N}} = \frac{24}{Re}\left(1 + 0.15\, Re^{0.687}\right)

Both go to $24/Re$ as $Re \to 0$ . The difference is the second term inside the parentheses.

Plug in $Re = 17$ . Prony gives $1.41 + 0.40 = 1.81$ . Schiller–Naumann gives $1.41 \times 2.05 = 2.90$ . Drag comes out 38% below the reference. The weight is the same, so the Prony bead must fall faster to balance. The result is a terminal velocity 42% too high.

Where does the missing drag come from? Once $Re$ passes 1, the fore–aft symmetric flow breaks. A slow wake forms behind the sphere, and around $Re \approx 20$ the surface flow separates into a recirculation zone. That wake pressure drag comes from viscous and inertial effects mixed together. It belongs cleanly to neither $v$ nor $v^2$ . Prony's sum has no slot for this mixed term.

That is why the sum misses most when the two terms are similar in size, right around $Re \approx 60$ . It is also why both historical experiments were right on their own terms. Both sat outside this window.

Where the two-term law still lives#

Despite the failure window, the $a v^2 + b v$ form never went away. The Ergun equation and the Darcy–Forchheimer law for porous media have exactly this shape. Prony's pipe-loss formula shares the same root. In porous media the two-term sum works reasonably well even across a wide range of pore-scale $Re$ , because the coefficients are fitted to experiments.

Euler–Lagrange particle codes are a different story. Single-particle $Re$ often falls right in this troublesome range (1–1000). Spray droplets, fluidized-bed particles, and settling sediment all do. Using $24/Re + 0.4$ there puts terminal velocity off by 30–40%. Settling and residence times shift by the same amount.

When choosing a drag model for your code, check in this order:

Extract the range of particle $Re$ first. Below 1, Stokes is enough; above 1000, a constant $C_D$ is enough.
If much of the mass sits between 1 and 1000, use a correlation with an intermediate term, such as Schiller–Naumann.
Run one single-bead settling test and compare with the 0.3 mm and 1 mm rows in the table above. A 30% error signals a hidden two-term sum.

Bossut's ships and Du Buat's pipes were the same water seen at different $Re$ . The gap between them cannot be filled by addition alone.

Adding a v Term to a v² Term Got Both Ends Right and Missed the Middle by 42% — Prony's Two-Term Drag Law and Re = 60