Sumar un término en v y otro en v² acertó en los extremos y falló un 42 % en el medio — la ley de arrastre de dos términos de Prony y Re = 60

Sumar arrastre lineal y cuadrático funciona para Re ≪ 60 y Re ≫ 60, pero entre Re 10 y 300 omite el arrastre de la estela y sobrestima la velocidad terminal en un 30–40 %.

La misma agua, el mismo vidrio, distinta ley de arrastre#

Se dejan caer dos esferas de vidrio en agua. Una mide 0,05 mm de diámetro y la otra 10 mm. La pequeña sufre un arrastre proporcional a su velocidad. La grande, uno proporcional al cuadrado de la velocidad. Misma agua, mismo material. Solo cambia el tamaño.

A finales del siglo XVIII, científicos franceses encontraron estas dos leyes por separado. Luego las fusionaron en una sola fórmula, $F = a v^2 + b v$ . Este artículo comprueba con números dónde se sostiene esa suma y dónde falla. Al final quedan claros el número de Reynolds (la razón entre fuerzas inerciales y viscosas) en el que los dos términos se relevan y el intervalo donde la suma simple se equivoca más.

El recorrido se divide en cinco preguntas.

P1. Entre 1775 y 1800, ¿el arrastre era v o v²?#

En 1752 d'Alembert chocó con su paradoja: un fluido ideal no ejerce arrastre sobre un cuerpo. Con la teoría estancada, los experimentos abrieron camino. En 1775 la Academia Francesa formó un equipo para estudiar la resistencia de los fluidos. Bossut dirigió el trabajo.

Bossut remolcó modelos de barcos en un tanque. Las mediciones mostraron un arrastre que crecía con el cuadrado de la velocidad. Mientras tanto, Du Buat midió el flujo dentro de tuberías. La fricción que observó crecía simplemente en proporción a la velocidad. Dos experimentos, dos respuestas.

Prony realizó experimentos más precisos y los unificó. La resistencia podía escribirse como $a v^2 + b v$ . Coulomb la organizó en dos tipos. La parte que crece con $v^2$ es el arrastre de forma (form drag, la diferencia de presión que impone la geometría del cuerpo). La que crece con $v$ es el arrastre de fricción (friction drag). Esa división alimentó después el trabajo de Navier y Stokes sobre flujo viscoso.

¿Quién tenía razón? Ambos. Los modelos de barco eran grandes y rápidos. Los flujos en tuberías eran lentos y estrechos. Los dos experimentos vieron extremos opuestos de la misma física. El número que separa esos extremos es la siguiente pregunta.

P2. ¿De dónde salen los dos términos?#

Escrito con dos términos, el arrastre sobre una esfera es

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

$\mu$ es la viscosidad, $\rho$ la densidad del fluido, $d$ el diámetro de la esfera y $C_\infty$ el coeficiente de arrastre a alta velocidad (cerca de 0,4 para una esfera).

El primer término es el arrastre de Stokes. Es la solución exacta del flujo lento sin inercia alguna (flujo reptante). El segundo es el arrastre de Newton. Se obtiene si el momento del fluido que golpea la cara frontal se convierte íntegramente en presión.

Una salvedad. El arrastre de Stokes $3\pi\mu d v$ no es fricción pura de superficie. Dos tercios son fricción y un tercio es presión. Por eso "dominado por viscosidad frente a dominado por inercia" es una etiqueta más precisa que la de Coulomb, "fricción frente a forma".

Conviene manipular la simulación de abajo. Al soltar la esfera, ambos arrastres crecen con la velocidad. El punto donde su suma iguala al peso es la velocidad terminal (terminal velocity, la velocidad a la que cesa la aceleración).

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%

Basta mover el control de diámetro de 0,02 mm a 20 mm. En esferas pequeñas la recta azul (b·v) corta la línea del peso; en las grandes lo hace la parábola naranja (a·v²). Al cambiar el fluido a glycerol, incluso una esfera de 10 mm pasa al lado de Du Buat.

P3. ¿Cuándo cambia el control de manos?#

Hay que hallar la velocidad a la que ambos términos se igualan.

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$ es el número de Reynolds basado en el diámetro de la esfera. Con $C_\infty = 0.4$ , $Re = 60$ .

Este resultado reconcilia a Bossut y Du Buat. La ley de arrastre no la fija el tipo de fluido ni la velocidad por sí sola. La fija $Re$ . El flujo lento en tuberías tenía $Re$ pequeño. Los modelos de barco operaban con $Re$ de cientos de miles. Cambiar agua por glycerol en la simulación anterior no era más que mover $Re$ .

El mismo problema de "unir dos límites" aparece en otros lugares. En el calentamiento de partículas en el régimen de transición de Knudsen, los límites continuo y de molécula libre se unían con una media armónica. En la velocidad del sonido de Wood, las compresibilidades se sumaban por fracción de volumen. La forma de mezclar los límites cambia el resultado. Prony eligió la suma simple.

P4. Entonces, ¿basta con sumar?#

Se puede comprobar. Se sueltan desde el reposo seis esferas de vidrio ( $\rho_p = 2500$ kg/m³) en agua. La ecuación de movimiento se integra con arrastre de dos términos.

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

$m$ es la masa de la esfera, $V$ su volumen y $g$ la gravedad. Se desprecia la masa añadida (added mass, la inercia del fluido que la esfera arrastra consigo).

La velocidad terminal se compara con una referencia: la correlación de Schiller–Naumann $C_D = \frac{24}{Re}(1 + 0.15 Re^{0.687})$ . Se sabe que coincide con los experimentos dentro de unos pocos por ciento para $Re < 1000$ .

import math
 
# Esferas de vidrio sedimentando en agua: arrastre de Du Buat (lineal) vs Bossut (cuadrático)
RHO_P, RHO, MU, G = 2500.0, 1000.0, 1.0e-3, 9.81
C_INF = 0.4          # coeficiente de arrastre cuadrático (de forma) de la ley de dos términos
 
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)       # escala de tiempo de la aceleración
    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

El código anterior puede ejecutarse directamente en el navegador. Dos esferas idénticas se hunden lado a lado. Una sufre el arrastre de dos términos de Prony; la otra, el de referencia de Schiller–Naumann.

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

Con el valor por defecto d = 0,30 mm se obtiene err 42,3 %, lin 77,8 % y t99/τ 3,92, igual que en el bloque de salida. La curva de la derecha es el error en todo el rango de diámetros. Al subir C∞ a 0,6 el error de las esferas grandes se vuelve negativo; al bajar ρp, Re_t disminuye para un mismo diámetro y el pico se desplaza a la derecha.

La salida dice tres cosas. Primero, ambos extremos ajustan bien: 2,9 % en 0,05 mm y 2,2 % en 3 mm. Segundo, usar un solo término colapsa en el extremo opuesto. Solo Stokes hace que la esfera de 10 mm caiga 117 veces más rápido. Solo Newton hace que la de 0,05 mm caiga 24 veces más rápido. Tercero, el problema está en el medio: 42 % de error en 0,3 mm ( $Re_t = 17$ ) y 33 % en 1 mm ( $Re_t = 193$ ).

El tiempo de aceleración también sigue la ley. Con arrastre solo lineal, la velocidad se acerca de forma exponencial, así que $t_{99}/\tau = \ln 100 = 4.61$ . Con arrastre solo cuadrático, se acerca como $\tanh$ , así que $t_{99}/\tau = \operatorname{artanh} 0.99 = 2.65$ . Los valores 4,59 y 2,65 de la salida son esos dos límites; los diámetros intermedios pasan entre ellos.

P5. ¿Por qué falla en el medio?#

Al reescribir la ley de dos términos como coeficiente de arrastre, la causa aparece.

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)

Ambas tienden a $24/Re$ cuando $Re \to 0$ . La diferencia está en el segundo término dentro del paréntesis.

Con $Re = 17$ , Prony da $1.41 + 0.40 = 1.81$ . Schiller–Naumann da $1.41 \times 2.05 = 2.90$ . El arrastre queda un 38 % por debajo de la referencia. El peso es el mismo, así que la esfera de Prony debe caer más rápido para equilibrarse. El resultado es una velocidad terminal un 42 % mayor.

¿De dónde sale el arrastre que falta? Cuando $Re$ supera 1, la simetría delante-detrás del flujo se rompe. Detrás de la esfera se forma una estela lenta y, cerca de $Re \approx 20$ , el flujo superficial se separa en una zona de recirculación. Ese arrastre de presión de la estela nace de efectos viscosos e inerciales mezclados. No pertenece limpiamente ni a $v$ ni a $v^2$ . La suma de Prony no tiene lugar para ese término mixto.

Por eso la suma falla más cuando ambos términos tienen tamaño parecido, justo alrededor de $Re \approx 60$ . También por eso los dos experimentos históricos acertaban cada uno a su manera. Ambos quedaban fuera de esa ventana.

Dónde sigue viva la ley de dos términos#

Pese a esa ventana de fallo, la forma $a v^2 + b v$ nunca desapareció. La ecuación de Ergun y la ley de Darcy–Forchheimer para medios porosos tienen exactamente esa forma. La fórmula de Prony para pérdidas en tuberías comparte la misma raíz. En medios porosos la suma de dos términos funciona razonablemente bien incluso en un rango amplio de $Re$ de poro, porque los coeficientes se ajustan a experimentos.

En los códigos de partículas Euler–Lagrange la historia es otra. El $Re$ de una partícula individual suele caer justo en este intervalo problemático (1–1000). Gotas de aerosol, partículas de lecho fluidizado y sedimentos lo hacen. Usar ahí $24/Re + 0.4$ desvía la velocidad terminal un 30–40 %. Los tiempos de sedimentación y de residencia se desplazan en la misma proporción.

Al elegir un modelo de arrastre para un código, conviene revisar en este orden:

Extraer primero el rango de $Re$ de las partículas. Por debajo de 1 basta Stokes; por encima de 1000, un $C_D$ constante.
Si gran parte de la masa cae entre 1 y 1000, usar una correlación con término intermedio, como Schiller–Naumann.
Correr una prueba de sedimentación de una sola esfera y compararla con las filas de 0,3 mm y 1 mm de la tabla anterior. Un error del 30 % delata una suma de dos términos escondida.

Los barcos de Bossut y las tuberías de Du Buat eran la misma agua vista a distinto $Re$ . El hueco entre ambos no se llena solo con una suma.