La fuerza de Coriolis hizo cero trabajo y el impulsor sumó 1800 J/kg — energía en un marco rotante

Quien definió el trabajo tiene una fuerza que no lo hace#

Gaspard-Gustave de Coriolis fue quien dejó "trabajo" y "energía cinética" en la forma que todavía se usa. Su libro de 1829 se titulaba Du calcul de l'effet des machines, sobre el cálculo del efecto de las máquinas. Lo escribió un ingeniero que contaba cuánto entraba y cuánto salía de una rueda hidráulica.

La fuerza que lleva su nombre no hace trabajo alguno. Su magnitud da igual. Este artículo comprueba esa afirmación con un solo impulsor centrífugo. La aceleración de Coriolis llega a 290 g, el libro contable sigue marcando $10^{-16}$ J/kg, y la entalpía total del fluido sube 1800 J/kg. De dónde salen esos 1800, y qué desaparece cuando un solver en marco rotante omite el término correspondiente, es la conclusión.

1832: de calcular ruedas hidráulicas a las coordenadas rotantes#

La serie coreana Mecánica de fluidos en la historia describe el París de esos años. La Revolución de Julio de 1830 expulsó al realista Cauchy, y el republicano Navier ocupó su puesto en la École Polytechnique. En 1832, Navier empezó a investigar junto con Coriolis.

El tema que Coriolis venía masticando eran las ruedas hidráulicas. ¿Cómo se cuentan la energía y el trabajo dentro de una máquina de fluidos que gira? Para hacer esa cuenta hay que sentarse sobre el eje que gira. De ahí nació su estudio de las coordenadas rotantes, y el artículo de 1835 dio a conocer la fuerza que hoy lleva su nombre.

El orden importa. La fuerza de Coriolis no salió de la rotación terrestre ni de la meteorología. Salió del intento de cuadrar el libro energético de una máquina rotante. En la simulación de abajo se puede manipular ese libro directamente.

Con el botón frame se cambia la cámara entre el marco relativo y el absoluto. La trayectoria naranja cambia de forma por completo, mientras la línea verde de la derecha se niega a inclinarse. Al subir omega solo se empina la curva naranja de $h_0$.

Una fuerza perpendicular no llega al libro contable#

En un marco que gira con $\vec{\Omega}$, una partícula de fluido con velocidad relativa $\vec{w}$ obedece

El segundo término del lado derecho es el de Coriolis y el tercero el centrífugo. Ambos son términos inerciales, el precio de haber puesto a girar el sistema de coordenadas. Cómo entran en la ecuación de cantidad de movimiento se trató en marcos de referencia rotantes y MRF. Aquí interesa solo el lado energético.

Para ver energía se hace el producto escalar de cada fuerza con $\vec{w}$. El término de Coriolis muere en el acto.

El producto vectorial es perpendicular a sus dos factores. Es una identidad: sin aproximación y sin condiciones. Vale para cualquier $\Omega$ y cualquier dirección de $\vec{w}$.

El término centrífugo se comporta distinto. Vale $\Omega^2 r\,\hat{r}$, así que cualquier velocidad radial mantiene vivo el producto escalar. A cambio, tiene potencial.

Tener potencial significa que puede pasarse al lado izquierdo y plegarse dentro de una constante.

Por eso lo conservado es $I$, no $h_0$#

Integrar esa ecuación a lo largo de una línea de corriente en condiciones estacionarias y no viscosas da la rotalpía (rothalpy, de rotational + enthalpy).

$h$ es la entalpía estática, $w$ el módulo de la velocidad relativa y $U=\Omega r$ la velocidad periférica del álabe. El signo menos del último término es el potencial centrífugo. El término de Coriolis no deja ni rastro en esta expresión, porque su producto escalar era cero desde el principio.

Su vínculo con la entalpía total absoluta $h_0 = h + c^2/2$ sale de $\vec{c} = \vec{w} + \vec{U}$:

donde $c_\theta$ es la componente de giro de la velocidad absoluta. Es decir, $h_0$ puede subir mientras $I$ se queda quieta, y sube exactamente $\Delta(U c_\theta)$: la ecuación de Euler para turbomáquinas.

Tres libros contables separados en Python#

Se toma un impulsor centrífugo que va de 50 mm a 150 mm de radio, con la altura del canal estrechándose de 20 mm a 8 mm. Con $\Omega = 300$ rad/s y 30 kg/s de agua. Dos ángulos de curvatura hacia atrás, $\beta = 0°$ y $30°$. La continuidad fija la velocidad relativa radial, la rotalpía constante fija la presión, y después se integran por separado en el tiempo los trabajos de la fuerza de Coriolis, la centrífuga y la reacción del álabe.

import numpy as np
 
OMEGA = 300.0            # velocidad angular [rad/s]
RHO = 1000.0             # densidad [kg/m^3]
MDOT = 30.0              # caudal másico [kg/s]
R1, R2 = 0.05, 0.15      # radio de entrada / salida [m]
B1, B2 = 0.020, 0.008    # altura del canal en entrada / salida [m]
GRID = np.linspace(R1, R2, 4001)
 
 
def channel_state(r, beta_deg):
    """En el radio r: componentes relativas (w_r, w_t), velocidad de álabe U, giro absoluto c_t."""
    b = B1 + (B2 - B1) * (r - R1) / (R2 - R1)
    w_r = MDOT / (RHO * 2.0 * np.pi * r * b)          # componente radial que fija la continuidad
    w_t = -w_r * np.tan(np.radians(beta_deg))         # álabe curvado atrás -> contra la rotación
    u = OMEGA * r
    return w_r, w_t, u, u + w_t
 
 
def march_channel(beta_deg):
    """Con la rotalpía constante, se levantan p/rho y la entalpía total absoluta h0 en el canal."""
    w_r, w_t, u, _ = channel_state(R1, beta_deg)
    i_const = 0.5 * (w_r**2 + w_t**2) - 0.5 * u**2    # tomando p1/rho = 0 como referencia
    cols = []
    for r in GRID:
        w_r, w_t, u, c_t = channel_state(r, beta_deg)
        p = i_const - 0.5 * (w_r**2 + w_t**2) + 0.5 * u**2
        cols.append((u, c_t, p + 0.5 * (w_r**2 + c_t**2),
                     p + 0.5 * (w_r**2 + w_t**2) - 0.5 * u**2))
    return np.array(cols)                              # U, c_t, h0, I
 
 
def power_ledger(beta_deg):
    """Integra en el tiempo, por separado, la potencia de Coriolis, centrífuga y de reacción."""
    om = np.array([0.0, 0.0, OMEGA])
    st = np.array([channel_state(r, beta_deg) for r in GRID])
    w_r, w_t = st[:, 0], st[:, 1]
    dwt_dr = np.gradient(w_t, GRID)
    p_cor = np.zeros_like(GRID)
    p_cen = np.zeros_like(GRID)
    p_bld = np.zeros_like(GRID)
    for k, r in enumerate(GRID):
        w = np.array([w_r[k], w_t[k], 0.0])            # base ortonormal local (r, theta, z)
        f_cor = -2.0 * np.cross(om, w)
        f_cen = -np.cross(om, np.cross(om, np.array([r, 0.0, 0.0])))
        acc_t = w_r[k] * dwt_dr[k] + w_r[k] * w_t[k] / r   # con el término de curvatura cilíndrica
        f_bld_t = acc_t - f_cor[1]                     # la aceleración de giro restante es del álabe
        p_cor[k] = f_cor @ w
        p_cen[k] = f_cen @ w
        p_bld[k] = OMEGA * r * f_bld_t                 # par x velocidad angular = potencia absoluta
    dt = 1.0 / w_r                                     # dt = dr / w_r
    ints = [np.trapezoid(p * dt, GRID) for p in (p_cor, p_cen, p_bld)]
    return ints + [np.max(np.abs(p_cor))]
 
 
def transported_h0(beta_deg, with_source):
    """Integra directamente el transporte de h0. El término fuente es Omega * d(r c_theta)/dr."""
    st = np.array([channel_state(r, beta_deg) for r in GRID])
    src = OMEGA * np.gradient(GRID * st[:, 3], GRID) if with_source else np.zeros_like(GRID)
    return np.trapezoid(src, GRID)
 
 
for beta in (0.0, 30.0):
    tab = march_channel(beta)
    d_h0 = tab[-1, 2] - tab[0, 2]
    euler = tab[-1, 0] * tab[-1, 1] - tab[0, 0] * tab[0, 1]
    drift = np.max(np.abs(tab[:, 3] - tab[0, 3]))
    cor, cen, bld, cor_peak = power_ledger(beta)
    ok = transported_h0(beta, True)
    bad = transported_h0(beta, False)
    print(f"beta = {beta:4.1f} deg   U1 = {tab[0,0]:4.1f} m/s   U2 = {tab[-1,0]:4.1f} m/s")
    print(f"  delta h0 from rothalpy = {d_h0:9.2f} J/kg")
    print(f"  U*c_theta (Euler)      = {euler:9.2f} J/kg   gap {abs(d_h0-euler):.1e}")
    print(f"  rothalpy max drift     = {drift:9.1e} J/kg")
    print(f"  work by Coriolis       = {cor:9.1e} J/kg   (peak power {cor_peak:.1e} W/kg)")
    print(f"  work by centrifugal    = {cen:9.2f} J/kg")
    print(f"  work by blade torque   = {bld:9.2f} J/kg")
    print(f"  h0 transport, source   = {ok:9.2f} J/kg   head {ok/9.81:5.1f} m")
    print(f"  h0 transport, dropped  = {bad:9.2f} J/kg   head {bad/9.81:5.1f} m")

beta =  0.0 deg   U1 = 15.0 m/s   U2 = 45.0 m/s
  delta h0 from rothalpy =   1800.00 J/kg
  U*c_theta (Euler)      =   1800.00 J/kg   gap 0.0e+00
  rothalpy max drift     =   1.1e-13 J/kg
  work by Coriolis       =   0.0e+00 J/kg   (peak power 0.0e+00 W/kg)
  work by centrifugal    =    900.00 J/kg
  work by blade torque   =   1800.00 J/kg
  h0 transport, source   =   1800.00 J/kg   head 183.5 m
  h0 transport, dropped  =      0.00 J/kg   head   0.0 m
beta = 30.0 deg   U1 = 15.0 m/s   U2 = 45.0 m/s
  delta h0 from rothalpy =   1737.98 J/kg
  U*c_theta (Euler)      =   1737.98 J/kg   gap 0.0e+00
  rothalpy max drift     =   7.1e-14 J/kg
  work by Coriolis       =   2.2e-16 J/kg   (peak power 1.4e-12 W/kg)
  work by centrifugal    =    900.00 J/kg
  work by blade torque   =   1737.98 J/kg
  h0 transport, source   =   1737.98 J/kg   head 177.2 m
  h0 transport, dropped  =      0.00 J/kg   head   0.0 m

La deriva de la rotalpía es de $10^{-13}$ J/kg, es decir, redondeo de doble precisión. La cuenta de Coriolis marca exactamente 0.0 con curvatura nula y $2.2\times10^{-16}$ a 30°. Ese segundo número es residuo de coma flotante, no física.

900 y 1800 — quién pagó la otra mitad#

Dos números llaman la atención. La fuerza centrífuga hizo 900 J/kg de trabajo, pero la entalpía total subió 1800 J/kg. Exactamente el doble.

Esos 900 son dinero que solo circula dentro del marco relativo. Se reparten entre presión y energía cinética relativa. Es una cuenta distinta de los 1800 que el fluido recibió de verdad en el marco absoluto.

Los otros 900 los pagó el álabe. En un álabe radial la velocidad relativa no tiene componente de giro. Para que siga así, algo debe cancelar exactamente la fuerza de Coriolis $2\Omega w_r$, y ese algo es la cara de presión del álabe. Por reacción el fluido recibe una fuerza de igual magnitud, cuyo par es $2\Omega r w_r$.

Aquí se bifurca. En el marco relativo esa fuerza es perpendicular a $\vec{w}$, así que no hace trabajo. En el marco absoluto el álabe gira a $\Omega$, de modo que par por velocidad angular es potencia. Al integrar sale $2\times 900 = 1800$ J/kg: la línea work by blade torque de la salida.

Dicho de forma directa, la fuerza de Coriolis no le entrega nada al fluido. Lo obliga a empujar el álabe, y el trabajo del eje entra por esa vía de reacción. La fuerza solo elige la dirección; quien paga es el eje.

La flecha naranja de Coriolis es la más larga de la pantalla y su cuenta sigue clavada en cero. Al arrastrar el deslizador backsweep solo cambia la cuenta morada del álabe. Al subir omega, la barra naranja sigue sin crecer.

Cuando un solver en marco rotante transporta $h_0$ tal cual#

Hasta aquí la física. En el código, el accidente ocurre en un sitio previsible.

Se trata de qué variable se transporta en la ecuación de energía sobre una zona rotante. Lo que la velocidad relativa transporta y se conserva sin fuente es $I$. Si se quiere llevar $h_0$ con la velocidad relativa, tiene que venir acompañada de un término fuente.

con $s$ la coordenada a lo largo de la línea de corriente. Al omitirlo se obtienen las dos últimas líneas de la salida. Con fuente: 1800 J/kg, altura de 183.5 m. Sin ella: 0.00 J/kg, altura de 0 m. El fluido atraviesa el impulsor y no le pasa nada.

El síntoma confunde porque los residuos se ven bien. La ecuación converge limpiamente, solo que converge a la respuesta con altura cero. Refinar la malla no ayuda: cero sigue siendo cero.

En los solvers de la familia OpenFOAM, esta forma aparece cuando la definición de energía total del lado rhothermo y la corrección MRF quedan desalineadas. Como término que falta en una transformación de coordenadas, comparte raíz con el error visto en los términos de curvatura del FVM polar. Por eso mismo la línea acc_t de arriba lleva el término de curvatura $w_r w_\theta / r$. Al quitarlo, la cuenta del álabe con 30° de curvatura se desvía a 1802.06.

Tres preguntas al encontrarse con una zona rotante#

Primera: en este caso, ¿qué se mantiene constante, $h_0$ o $I$? Dentro de la zona rotante es $I$. Al pasar a la zona estacionaria vuelve a ser $h_0$. Conviene revisar qué se está empalmando en la interfaz.

Segunda: ¿el término fuente de la ecuación de energía está realmente conectado? Si la altura sale extrañamente baja, o exactamente cero, ahí es donde hay que mirar primero.

Tercera: ¿se está "añadiendo" el término de Coriolis como fuente de energía? Su potencia es idénticamente cero. Si da la impresión de que ahí falta algo, lo que hace falta no es Coriolis sino $\Omega\,\partial(rc_\theta)/\partial s$.

Coriolis llegó a las coordenadas rotantes contando la eficiencia de ruedas hidráulicas. En el libro contable que él inventó, la fuerza con su nombre siempre marca cero. Mantener ese cero en cero al calcular máquinas rotantes es la forma en que heredamos su contabilidad, 190 años después.