Un 9% menos de coeficiente de arrastre dejó la gota un 35% más lenta — dónde se enganchan las correcciones de deformación y vaporización

La vaporización recorta $C_D$, pero la desaceleración se monta sobre $C_D/d$. En cuanto el diámetro empieza a encogerse, el coeficiente ahorrado regresa.

En 1983, el coeficiente de arrastre no cuadraba frente a una gota en llamas#

Renksizbulut y Yuen midieron el arrastre de gotas que se vaporizaban dentro de una corriente de gas caliente. El arrastre medido resultó sistemáticamente menor que el previsto por la correlación esférica. Y seguía siendo menor aunque las gotas conservaran la forma esférica. El vapor que salía de la superficie estaba empujando la capa límite hacia afuera.

Por esa misma época se reportaba el desajuste contrario. Una gota expuesta a una corriente rápida se achata, y achatada arrastra más que una esfera. Esa es la rama que O'Rourke y Amsden empaquetaron como modelo TAB en 1987.

Por eso el coeficiente de arrastre de gota en un código de spray suele tener tres capas. Una correlación esférica multiplicada por una corrección de deformación y una corrección de vaporización. Este artículo separa las tres capas y mide cuánto desplaza cada una la trayectoria en condiciones idénticas. La conclusión por delante: el modelo que más recorta el coeficiente es el que se detiene primero.

Al subir juntos el control de velocidad relativa y el diámetro, el aplastamiento de la gota y la barra de arrastre se mueven en el mismo sentido.

relative speed25 m/sdiameter40 umgas density5.0

y 0.000y_ss 0.000C_D x1.00We_r 0.00Re 0

Raise the speed and watch the drop flatten into a disc while the amber bar climbs past the grey sphere value. The oscillation overshoots the dashed steady line by almost 2x, so the peak drag is roughly twice what an equilibrium estimate would give. Turn liquid damping off and the ring never settles. Push far enough and y hits 1 — TAB stops and calls it breakup.

El instante en que la barra supera el valor gris de la esfera es la corrección de deformación. Conviene mirar también cuánto se pasa la deformación real de la línea punteada estacionaria. Subir solo el diámetro apenas mueve la barra. La velocidad relativa es lo que la hace saltar.

Tres correlaciones corrigen el mismo punto de formas distintas#

El punto de partida es el coeficiente de arrastre de una gota esférica. Se usa una correlación del tipo Schiller–Naumann.

$C_{D,\mathrm{sph}} = \frac{24}{Re}\left(1 + \tfrac{1}{6} Re^{2/3}\right), \qquad Re < 1000$

$Re = \rho_g u_{\mathrm{rel}} d / \mu_g$ es el número de Reynolds construido con el diámetro de la gota y la velocidad relativa. El término $24/Re$ es el arrastre de Stokes y el paréntesis es la corrección inercial.

Las dos correcciones que se acoplan aquí tienen caracteres opuestos.

Corrección	Factor multiplicativo	Dirección	Fundamento
Ninguna (esfera)	$1$	—	Esfera rígida, sin transferencia de masa en superficie
Deformación TAB	$1 + 2.632\,y$	arrastre sube	Crece el área frontal
Vaporización (R–Y)	$(1+B_M)^{-0.2}$	arrastre baja	El vapor empuja la capa límite

$y$ es la deformación adimensional del modelo TAB. $y = 0$ es una esfera y $y = 1$ es el criterio de ruptura. $B_M$ es el número de transferencia de masa de Spalding (la diferencia de fracción másica de vapor entre superficie y campo lejano, adimensionalizada), el mismo número tratado en el término fuente de evaporación y la ley $d^2$ .

La deformación se monta sobre el coeficiente, no sobre el área frontal#

TAB trata la gota como un oscilador armónico amortiguado. La presión dinámica del gas empuja, la tensión superficial restituye y la viscosidad del líquido amortigua.

\frac{d^2 y}{dt^2} = \frac{C_F}{C_b}\frac{\rho_g u_{\mathrm{rel}}^2}{\rho_l r^2} - \frac{C_k \sigma}{\rho_l r^3} y - \frac{C_d \mu_l}{\rho_l r^2}\frac{dy}{dt}

$r$ es el radio de la gota, $\sigma$ la tensión superficial y $\mu_l$ la viscosidad del líquido. Las constantes empleadas son $C_F = 1/3$ , $C_k = 8$ , $C_d = 5$ , $C_b = 1/2$ . El proceso por el cual esta ecuación llega hasta la ruptura se trató en el número de Weber y la ruptura TAB.

Anulando todas las derivadas se obtiene la deformación estacionaria.

$y_{ss} = \frac{C_F}{C_b C_k}We_r = \frac{We_r}{12}, \qquad We_r = \frac{\rho_g u_{\mathrm{rel}}^2 r}{\sigma}$

Aquí hay un punto en el que los practicantes se equivocan seguido. Es común que un código meta $y_{ss}$ directamente en la corrección de arrastre. Pero este oscilador casi no está amortiguado. La razón de amortiguamiento de una gota de n-heptano de 40 μm es del orden de $10^{-2}$ . Partiendo del reposo, $y$ rebasa $y_{ss}$ y sube hasta casi el doble.

Es decir, el arrastre calculado con el valor estacionario apenas alcanza la mitad de lo que entrega la primera oscilación. El problema no es el coeficiente sino el instante en que se lee ese coeficiente.

La vaporización empuja la capa límite y recorta el coeficiente#

La correlación de Renksizbulut–Yuen trata la vaporización como una corrección multiplicativa.

$C_D = \frac{24}{Re}\left(1 + 0.2\,Re^{0.63}\right)(1+B_M)^{-0.2}$

El exponente $-0.2$ es pequeño. Con $B_M = 0.6$ el factor vale $0.91$ ; con $B_M = 2$ vale $0.80$ . Incluso una gota que arde con violencia pierde apenas un 20% del coeficiente.

El principio es el mismo de la teoría de película (film theory). El flujo de vapor que sale de la superficie engrosa la capa límite. El gradiente de velocidad se suaviza, el esfuerzo cortante en la pared baja y el arrastre disminuye. El número de Nusselt de la transferencia de calor recibe una corrección de la misma forma.

Cuatro gotas disparadas desde la misma tobera en Python#

Se dispara una gota en gas en reposo y se la sigue durante 1 ms. La ecuación de gobierno es unidimensional.

$\frac{du_d}{dt} = -\frac{3}{4}\frac{\rho_g}{\rho_l}\frac{C_D}{d}\,u_d|u_d|$

Los cuatro casos solo difieren en el coeficiente de arrastre. A es la esfera, B es la esfera multiplicada por la deformación TAB, C es solo la corrección de vaporización y D suma a esa corrección la contracción por la ley $d^2$ .

import math
 
RHO_G, MU_G = 4.98, 3.3e-5      # aire a 700 K, 10 bar
RHO_L, SIG, MU_L = 620.0, 0.0145, 3.0e-4   # n-heptano
D0, U0 = 40e-6, 25.0            # gota tras la ruptura, velocidad de inyeccion
K_EVAP = 1.0e-6                 # constante de la ley d^2 [m^2/s]
CF, CK, CD_T, CB = 1.0 / 3.0, 8.0, 5.0, 0.5   # constantes TAB
 
def cd_sphere(re):
    if re < 1e-8:
        return 0.0
    return 24.0 / re * (1.0 + re ** (2.0 / 3.0) / 6.0) if re < 1000.0 else 0.424
 
def cd_distorted(re, y):
    return cd_sphere(re) * (1.0 + 2.632 * max(0.0, min(1.0, y)))
 
def cd_blowing(re, bm):
    if re < 1e-8:
        return 0.0
    return 24.0 / re * (1.0 + 0.2 * re ** 0.63) * (1.0 + bm) ** -0.2
 
def tab_oscillate(y, ydot, u, d, dt):
    r = 0.5 * d
    acc = (CF / CB) * RHO_G * u * u / (RHO_L * r * r) \
        - (CK * SIG / (RHO_L * r ** 3)) * y \
        - (CD_T * MU_L / (RHO_L * r * r)) * ydot
    ydot += acc * dt
    y += ydot * dt
    return y, ydot
 
def run_droplet(law, evaporating, bm=0.6, t_end=1.0e-3, dt=5e-8):
    u, x, d, t = U0, 0.0, D0, 0.0
    y, ydot, ymax = 0.0, 0.0, 0.0
    cd_sum, cd0, n = 0.0, None, 0
    while t < t_end and u > 1e-6:
        re = RHO_G * u * d / MU_G
        if law == 'sphere':
            cd = cd_sphere(re)
        elif law == 'tab':
            y, ydot = tab_oscillate(y, ydot, u, d, dt)
            ymax = max(ymax, y)
            cd = cd_distorted(re, y)
        else:
            cd = cd_blowing(re, bm)
        if cd0 is None:
            cd0 = cd
        cd_sum += cd
        n += 1
        u += -0.75 * RHO_G * cd * u * u / (RHO_L * d) * dt
        x += u * dt
        if evaporating:
            d = math.sqrt(max(1e-18, d * d - K_EVAP * dt))
        t += dt
    return dict(x_mm=x * 1e3, u=u, d_um=d * 1e6, cd0=cd0,
                cd_mean=cd_sum / n, ymax=ymax)
 
CASES = [
    ('A sphere            ', 'sphere', False),
    ('B sphere x TAB      ', 'tab', False),
    ('C blowing (B_M=0.6) ', 'blow', False),
    ('D blowing + d2-law  ', 'blow', True),
]
 
print('Re0 = %.0f   We_r = %.2f   t = 1.0 ms' %
      (RHO_G * U0 * D0 / MU_G, RHO_G * U0 * U0 * 0.5 * D0 / SIG))
print('case                  Cd(t=0)  <Cd>    u(1ms)   x(1ms)   d(1ms)')
base = None
for name, law, ev in CASES:
    r = run_droplet(law, ev)
    if base is None:
        base = r
    print('%s  %6.3f  %6.3f  %6.2f   %6.3f   %5.1f  (%+.1f%% x)' %
          (name, r['cd0'], r['cd_mean'], r['u'], r['x_mm'], r['d_um'],
           100.0 * (r['x_mm'] / base['x_mm'] - 1.0)))
 
tab = run_droplet('tab', False)
print('TAB peak distortion y_max = %.3f  ->  Cd multiplier %.2fx' %
      (tab['ymax'], 1.0 + 2.632 * tab['ymax']))

Re0 = 151   We_r = 4.29   t = 1.0 ms
case                  Cd(t=0)  <Cd>    u(1ms)   x(1ms)   d(1ms)
A sphere               0.910   1.708    3.36    9.261    40.0  (+0.0% x)
B sphere x TAB         0.910   2.233    2.66    7.490    40.0  (-19.1% x)
C blowing (B_M=0.6)    0.828   1.537    3.68    9.732    40.0  (+5.1% x)
D blowing + d2-law     0.828   2.096    2.18    8.868    24.5  (-4.3% x)
TAB peak distortion y_max = 0.632  ->  Cd multiplier 2.66x

Al subir el control de $B_M$ y la constante de evaporación $K$ , solo se separa el alcance de los carriles C y D.

injection u025 m/sd040 umB_M0.60K1.00 mm2/s

A 0.00 mm (--)B 0.00 mm (-100.0%)C 0.00 mm (-100.0%)D 0.00 mm (-100.0%)

Push B_M up: lane C's drag coefficient drops and it pulls ahead of the grey sphere. Now push K up and lane D uses that same reduced coefficient yet falls behind — the shrinking diameter multiplies the deceleration faster than blowing removes it. Setting K to 0 collapses D onto C, which is the cleanest way to see which of the two effects is doing the work.

Subir $B_M$ hace que C adelante a la esfera gris. Subir $K$ desde ahí deja atrás solo a D, que usa el mismo coeficiente. Al volver $K$ a 0, D cae encima de C. Ahí se separan el +5.1% y el -4.3% de la salida impresa.

La corrección de deformación reduce la penetración un 19.1%, porque en la primera oscilación el coeficiente sube hasta 2.66 veces su valor base. Aplicando solo la corrección de vaporización ocurre lo contrario y la penetración crece un 5.1%. Hasta aquí todo sigue el signo de cada correlación.

Por qué la gota con el coeficiente más bajo se detuvo primero#

El problema es D. Usa exactamente la misma expresión del coeficiente de arrastre que C. Su coeficiente inicial también vale 0.828. Sin embargo, tras 1 ms su velocidad es de 2.18 m/s, la más baja de las cuatro. Un 35% más lenta que los 3.36 m/s del caso A.

Releer la ecuación de desaceleración da la respuesta. En el lado derecho aparece $C_D/d$ . Mientras el coeficiente de arrastre se recortaba un 9%, el diámetro cayó de 40 μm a 24.5 μm, un 39% menos. El divisor se movió mucho más que el numerador.

Físicamente ocurre así. El arrastre escala con el área ( $d^2$ ) y la inercia con la masa ( $d^3$ ). La razón va como $1/d$ . Cuanto más se seca la gota, más desaceleración produce el mismo arrastre. Equivale a decir que el tiempo de respuesta $\tau_d \simeq \rho_l d^2 / 18\mu_g$ se reduce con $d^2$ .

$\frac{a}{m} \propto \frac{d^2}{d^3} = \frac{1}{d}$

La vaporización tiene entonces dos efectos simultáneos: baja el coeficiente de arrastre y agranda la desaceleración. El primero queda atado a un exponente de $-0.2$ ; el segundo actúa directamente sobre el diámetro. Basta integrar un poco más de tiempo para que gane el segundo.

La razón entre los coeficientes iniciales 0.910 y 0.828 es 0.91. Con $B_M = 0.6$ , $(1+B_M)^{-0.2}$ es justo ese 0.91, y 0.910 multiplicado por él da 0.828. Con $B_M = 2$ el factor es 0.80, un recorte del 20%. El 9% de este caso es lo que ese exponente débil realmente compra.

La deformación va al revés. $We_r = 4.29$ , así que $y_{ss} = We_r/12$ vale 0.358. El $y_{\max} = 0.632$ impreso es cerca de 1.8 veces ese valor. Al meter 0.632 en $1 + 2.632\,y$ el multiplicador sale 2.66, el de la tabla. Con el valor estacionario 0.358 el multiplicador se queda en 1.94. Esa diferencia separa el coeficiente medio de B, 2.233, del de A, 1.708, y recorta la penetración de 9.261 mm a 7.490 mm, el −19.1%.

La posición de D queda solo un 4.3% por debajo de A, pero la velocidad cae un 35%: 8.868 mm frente a 9.261 mm, 2.18 m/s frente a 3.36 m/s. Al principio el diámetro sigue en 40 μm, así que D viaja junto con C. El frenado se acumula al llegar a 1 ms, cuando el diámetro ya es 24.5 μm. Restar $(24.5\,\mu\mathrm{m})^2$ de $(40\,\mu\mathrm{m})^2$ deja unos 1,000 μm², que es $K = 1.0\times 10^{-6}$ m²/s por 1 ms. El coeficiente medio también sube, de 1.537 en C a 2.096 en D. Un diámetro menor baja $Re$ y sube $24/Re$ . La fórmula que recortó el coeficiente un 9% termina con un promedio temporal más alto.

Qué correlación activar y cuándo#

Situación	Corrección de deformación	Corrección de vaporización	Motivo
$We_r < 0.5$ , baja temperatura	apagada	apagada	$y_{ss} < 0.04$ , corrección menor al 4%
Campo cercano de spray diésel	encendida	encendida	Velocidad relativa máxima, domina la deformación
Campo lejano de spray	apagada	encendida	$We$ cae tras la desaceleración, solo queda la vaporización
Partículas sólidas	apagada	apagada	Ni deformación ni transferencia de masa
Cerca del punto crítico	apagada	con cuidado	$\sigma \to 0$ , TAB sin sentido · $B_M$ diverge

La última fila da problemas reales con frecuencia. Cuando la tensión superficial tiende a cero, el término restitutivo de TAB desaparece y $y$ diverge. En ese régimen hay que apagar la corrección de deformación y pasar a un modelo de interfaz difusa.

Si el arrastre sigue sin cuadrar aunque $We$ sea suficientemente bajo, hay que mirar en otra parte y no en las correcciones. Puede tratarse de un caso como la cavitación por debajo de la presión de vapor, donde la fase que la correlación daba por sentada ya no es esa fase.

Cambiar una correlación mueve el caso de validación entero#

En este cálculo la penetración osciló entre -19.1% y +5.1%. No se tocó la malla, ni el paso de tiempo, ni el modelo de turbulencia. Todo el experimento consistió en encender y apagar dos factores adimensionales multiplicados sobre una sola gota.

Por eso, cuando una validación de spray discrepa del experimento, apretar primero la malla no es el orden correcto. Primero hay que verificar cuál es la expresión del coeficiente de arrastre, con qué valores entran $y$ y $B_M$ en ella, y si $y$ salió del estado estacionario o de resolver la oscilación. Si el diámetro se está reduciendo con el tiempo, conviene imprimir $C_D/d$ en lugar de $C_D$ . Mirando solo el coeficiente nunca se llega a ver por qué la gota se detuvo primero.