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.

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$ .

A continuación se puede hacer crecer la deformación de forma manual.

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.

Al subir la velocidad relativa y el diámetro, la gota se aplasta hasta formar un disco y la barra del coeficiente supera el valor esférico gris. Lo que conviene observar es cuánto sobrepasa la $y$ real al valor estacionario $y_{ss}$ marcado con línea punteada.

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

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.

A continuación se disparan las cuatro gotas desde la misma tobera.

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.

Al subir $B_M$ , el carril C adelanta a la esfera gris. Desde ese estado, al subir $K$ solo se rezaga el carril D, que usa el mismo coeficiente. Volviendo a $K = 0$ , D se superpone exactamente sobre C, de modo que los dos efectos quedan separados de una sola vez.

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. Cerca de una pared puede tratarse de un caso como la crisis de arrastre y la separación de la capa límite, donde el propio $Re$ se salió del rango de validez de la correlación.

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.