Donde la presión total cayó 1,080 Pa la disipación era cero — cuando la constante de Bernoulli cruza las líneas de corriente

En 1738 un padre adelantó seis años el libro de su hijo#

Daniel Bernoulli publicó Hydrodynamica en 1738. En la portada se presentó como "hijo de Johann", un gesto para reconciliarse con su padre. Johann respondió publicando casi el mismo contenido por separado, bajo el título Hydraulica, y presionó al impresor para que estampara el año 1732. Sobre el papel, seis años antes que su hijo.

Siglo y medio después llegó el teorema de la mano de Horace Lamb: al integrar las ecuaciones de Euler aparece el resultado de Bernoulli. Los dos libros por los que pelearon padre e hijo eran dos caras de una misma expresión. Pero ese "al integrar" arrastra una condición, y esa condición es justamente lo que hace que los contornos de presión total se lean mal en el posproceso de CFD.

Este artículo mide la condición con un solo vórtice de Rankine. La conclusión primero: dentro del núcleo en rotación la presión total cae 1,080 Pa y la disipación viscosa es exactamente cero. Y afuera, donde la presión total es perfectamente plana, la disipación no es cero. Las dos salen al revés.

En la simulación de abajo conviene mover la sonda hacia dentro y hacia fuera del núcleo.

Al desplazar probe r, hay que observar dónde se aplana la curva verde de $p_0$ del gráfico derecho. Y en paralelo, mirar el cuadrado del lado izquierdo del lienzo: cuándo conserva su forma y cuándo se deforma. El punto donde ambas observaciones se contradicen es el tema de este artículo.

Integrar las ecuaciones de Euler da Bernoulli, ¿pero en qué dirección?#

Para un flujo estacionario no viscoso, con densidad constante y sin fuerzas de volumen, la ecuación de cantidad de movimiento queda así.

$\mathbf{u}$ es la velocidad, $p$ la presión y $\rho$ la densidad. Al miembro izquierdo se le aplica una identidad vectorial.

$\boldsymbol{\omega} = \nabla\times\mathbf{u}$ es la vorticidad. Al juntar ambas expresiones y agrupar la presión total $p_0 = p + \tfrac{1}{2}\rho|\mathbf{u}|^2$ queda una sola línea.

Esta es la forma de Crocco de la ecuación de Euler (condición incompresible e isentrópica). Si el miembro derecho no es cero, la presión total varía en el espacio. Y lo decisivo es en qué dirección varía.

$\mathbf{u}\times\boldsymbol{\omega}$ es perpendicular a $\mathbf{u}$. Por eso, al multiplicar escalarmente ambos miembros por $\mathbf{u}$, el lado derecho desaparece.

A lo largo de una línea de corriente $p_0$ siempre es constante. Da igual que haya vorticidad o no. Ese es el alcance exacto del teorema de Bernoulli. En cambio, para que $p_0$ sea igual en todas partes hace falta que $\mathbf{u}\times\boldsymbol{\omega} = 0$, lo que en la práctica significa flujo irrotacional. La diferencia entre los dos libros en disputa nunca estuvo aquí: la condición estaba en la nota al pie que agregó Lamb.

Vórtice de Rankine: un núcleo que gira y un exterior que no#

El vórtice de Rankine es el ejemplo donde esta condición se enciende y se apaga dentro de un mismo flujo. Adentro del radio $a$ el fluido gira como un sólido rígido; afuera es un vórtice libre.

$\Omega$ es la velocidad angular del núcleo y $u_\theta$ la velocidad circunferencial. La presión se obtiene integrando el balance radial de cantidad de movimiento $\mathrm{d}p/\mathrm{d}r = \rho u_\theta^2 / r$. Afuera, lo que baja la presión y lo que sube la presión dinámica se cancelan exactamente, y $p_0$ queda fijada en $p_\infty$. Dentro del núcleo no hay cancelación.

En $r = 0$ el déficit vale $\rho\Omega^2 a^2$, es decir, el cuadrado de la velocidad en el borde del núcleo multiplicado por $\rho$. Con $\rho = 1.2$, $\Omega = 60\,\mathrm{s^{-1}}$ y $a = 0.5\,\mathrm{m}$ resulta $u_{\max} = 30$ m/s y el déficit es de 1,080 Pa.

Los dos miembros medidos en el mismo radio con Python#

Se calcula $\nabla p_0$ por diferencias centradas y se compara con $\rho\,u_\theta\omega_z$. La función de disipación $\phi = \mu(2S_{r\theta})^2$ se mide en el mismo paso.

import math
 
RHO, MU = 1.2, 1.8e-5      # kg/m^3, Pa*s
OMEGA, A = 60.0, 0.5       # 1/s, m   -> u_max = 30 m/s
P_INF = 101325.0           # Pa
 
def rankine_velocity(r):
    return OMEGA * r if r < A else OMEGA * A * A / r
 
def rankine_spin(r):
    return 2.0 * OMEGA if r < A else 0.0
 
def rankine_pressure(r):
    if r >= A:
        return P_INF - 0.5 * RHO * (OMEGA * A * A / r) ** 2
    return P_INF - RHO * OMEGA**2 * A**2 + 0.5 * RHO * OMEGA**2 * r**2
 
def total_head(r):
    return rankine_pressure(r) + 0.5 * RHO * rankine_velocity(r) ** 2
 
def shear_dissipation(r):
    s = 0.0 if r < A else -OMEGA * A * A / r**2   # 2*S_rtheta
    return MU * s * s
 
def crocco_residual(r, h=1e-6):
    lhs = (total_head(r + h) - total_head(r - h)) / (2 * h)   # d(p0)/dr
    rhs = RHO * rankine_velocity(r) * rankine_spin(r)         # rho*(u x omega)_r
    return lhs, rhs
 
print("  r[m]   u[m/s]   p[Pa]      p0[Pa]     dp0/dr    rho*u*w    phi[W/m^3]")
for r in (0.10, 0.25, 0.40, 0.60, 1.00):
    lhs, rhs = crocco_residual(r)
    print("%6.2f %8.2f %10.1f %10.1f %9.1f %10.1f %11.4f"
          % (r, rankine_velocity(r), rankine_pressure(r), total_head(r),
             lhs, rhs, shear_dissipation(r)))
 
def bernoulli_gap(r1, r2):
    return total_head(r2) - total_head(r1)
 
print()
print("core  p0(0.00) - p0(%.2f) = %8.1f Pa" % (A, total_head(0.0) - total_head(A)))
print("outer p0(%.2f) - p0(1.00) = %8.1f Pa" % (A, total_head(A) - total_head(1.0)))
print("cross-streamline gap r=0.10 -> 0.40 : %8.1f Pa" % bernoulli_gap(0.10, 0.40))
print("along-streamline  gap r=0.40 -> 0.40 : %8.1f Pa" % bernoulli_gap(0.40, 0.40))
print("dissipation at r=0.25 (core)  : %.4f W/m^3" % shear_dissipation(0.25))
print("dissipation at r=0.60 (outer) : %.4f W/m^3" % shear_dissipation(0.60))

  r[m]   u[m/s]   p[Pa]      p0[Pa]     dp0/dr    rho*u*w    phi[W/m^3]
  0.10     6.00   100266.6   100288.2     864.0      864.0      0.0000
  0.25    15.00   100380.0   100515.0    2160.0     2160.0      0.0000
  0.40    24.00   100590.6   100936.2    3456.0     3456.0      0.0000
  0.60    25.00   100950.0   101325.0       0.0        0.0      0.0313
  1.00    15.00   101190.0   101325.0       0.0        0.0      0.0040
 
core  p0(0.00) - p0(0.50) =  -1080.0 Pa
outer p0(0.50) - p0(1.00) =      0.0 Pa
cross-streamline gap r=0.10 -> 0.40 :    648.0 Pa
along-streamline  gap r=0.40 -> 0.40 :      0.0 Pa
dissipation at r=0.25 (core)  : 0.0000 W/m^3
dissipation at r=0.60 (outer) : 0.0313 W/m^3

La cuarta y la quinta columna coinciden en todos los radios. La relación de Crocco se cumple hasta los decimales. Y dentro del núcleo la diferencia de presión total entre $r = 0.10$ y $r = 0.40$ es de 648 Pa. Entre dos ángulos del mismo radio es cero. Aplicar Bernoulli a lo largo de una línea de corriente funciona; aplicarlo cruzándolas deja escapar 648 Pa.

Donde la presión total era plana la disipación no era cero#

La última columna está invertida. La disipación es cero dentro del núcleo, donde la presión total cae 1,080 Pa, y no es cero afuera, donde la presión total es exactamente plana.

El motivo es que la disipación va unida a la deformación, no a la rotación. La rotación de sólido rígido solo hace girar al elemento de fluido, sin deformarlo. Como el tensor de velocidad de deformación es cero, $\phi = 0$. El vórtice libre es lo contrario: la vorticidad es cero, pero al ser $u_\theta \propto 1/r$ el lado interno y el externo pasan a velocidades distintas. El elemento se cizalla sin parar.

El cuadrado de la simulación anterior lo muestra tal cual. Dentro del núcleo (rojo) el cuadrado gira conservando su forma, y afuera (violeta) se derrumba en un paralelogramo. Qué se conserva y qué no en un sistema en rotación ya apareció en fuerza de Coriolis y rotalpía, y aquí el eje es el mismo. Girar y hacer trabajo se anotan en libros contables distintos.

Cinco flujos en una misma tabla#

La vorticidad y la disipación son independientes. Las cuatro combinaciones existen de verdad.

Lo importante es que en la tabla hay una sola celda en negrita. Lo único que merece llamarse pérdida es que $p_0$ disminuya a lo largo de una línea de corriente. En las otras cuatro filas, que $p_0$ varíe espacialmente es pura geometría de la vorticidad, no energía que se haya esfumado.

Cuando dos contornos de presión total se parecen por motivos distintos#

En la práctica esta distinción se rompe en el contorno de presión total del plano de salida. Al imponer como condición de entrada una capa límite atmosférica o un perfil de tubería desarrollado, $p_0$ varía en la dirección $y$ desde la primera celda. La pérdida todavía es cero. El déficit de presión total que genera una estela se dibuja con el mismo contorno. Y ese sí es pérdida real.

Abajo se pueden correr los dos canales en paralelo con la misma escala de color.

Al pulsar match outlet spread, los rangos de $p_0$ en la salida de ambos canales se igualan. Es decir, con el contorno de salida solo ya no se pueden distinguir. En su lugar hay que fijarse en si el color de las partículas cambia mientras avanzan. Las del canal superior mantienen su color, y solo las del canal inferior cambian de color al pasar el cuerpo.

Por eso, si el coeficiente de pérdida se calcula tomando como referencia el promedio de $p_0$ en la sección, la entrada cizallada arroja un valor distinto de cero aunque no haya pérdida alguna. La referencia tiene que ser el $p_0$ de entrada de esa misma línea de corriente. Y aun usando promedios ponderados por caudal másico, entrada y salida deben promediarse del mismo modo para que la diferencia sea solo pérdida.

Numéricamente se agrega algo más. Con una malla gruesa, $p_0$ se aplana artificialmente en las zonas de rotación. Si la difusión numérica difumina la vorticidad, el miembro derecho de $\nabla p_0 = \rho\,\mathbf{u}\times\boldsymbol{\omega}$ se reduce y el déficit del núcleo sale menos profundo de lo que es. En una malla que atraviesa el núcleo de un vórtice, la profundidad del déficit de presión total sirve como indicador de la resolución de vorticidad. Igual que en la validación en tubería, donde el error de diámetro entraba al caudal en cuarta potencia, también aquí el error entra por una magnitud que cuesta ver.

Con lo mismo se topó d'Alembert en 1752#

Johann nunca terminó de aceptar la teoría de la viscosidad de Newton, y su hijo Daniel y su discípulo Euler tampoco. Si el flujo alrededor de un cuerpo se resuelve solo con teoría no viscosa, la resistencia sale cero. Por eso la comunidad entró en pánico cuando d'Alembert lo publicó en 1752. El problema con el que d'Alembert chocó en esa misma época en la ecuación de onda tiene una única raíz común: hasta dónde se admite el conjunto de soluciones que la ecuación permite.

Escrito con la notación que se usa hoy, queda así. En flujo irrotacional y no viscoso $p_0$ es una sola constante en todo el dominio, con lo cual la presión delante y detrás del cuerpo resulta simétrica y la integral da cero. Para que haya resistencia, $p_0$ tiene que caer en algún lado a lo largo de una línea de corriente. Quien crea ese lugar es la viscosidad y la vorticidad que produce en la pared.

Así que, al abrir un contorno de presión total, la pregunta no es "cuánto cayó". Es "¿cayó a lo largo de la línea de corriente o cruzándola?". Lo primero es pérdida; lo segundo, vorticidad.