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. Un solo impulsor centrífugo basta para rehacer el libro. 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, lo dicen los números.

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. Al pasar el botón frame de relativo a absoluto solo cambia de forma la trayectoria naranja. La línea verde de la derecha se queda horizontal.

La trayectoria naranja se dobla con la cámara. La línea verde no. Al subir omega solo se empina la curva 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

Al arrastrar backsweep de 0° a 30° solo baja el trabajo morado del álabe, de 1800.00 a 1737.98 J/kg.

La flecha naranja de Coriolis es la marca más larga de la pantalla y su libro sigue en cero. Al subir omega esa barra tampoco crece.

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.

La velocidad periférica de entrada es 15.0 m/s y la de salida 45.0 m/s. Sustituidas en $(U_2^2 - U_1^2)/2$, el cuadrado de 45 es 2025, el de 15 es 225, y la mitad de la diferencia es 900. El trabajo centrífugo impreso, 900.00 J/kg, coincide con esa expresión hasta el centésimo. Con curvatura nula el trabajo del par del álabe es 1800.00, exactamente el doble del centrífugo. La potencia en el marco absoluto es par por $\Omega$, así que los 900 del libro relativo se cuentan dos veces. A 30° el trabajo del álabe baja a 1737.98 J/kg y el término de Euler $U c_\theta$ baja al mismo 1737.98. El gap imprime 0.0e+00. El trabajo centrífugo se queda en 900.00. Solo se movió el giro. $w_t$ apunta contra la rotación, de modo que $c_\theta$ de salida es menor que $U$ y $\Delta(U c_\theta)$ queda por debajo de 1800. La deriva de la rotalpía es 1.1e-13 J/kg a 0° y 7.1e-14 a 30°. Es $I$ sostenida como constante dentro de la doble precisión, no una pérdida nueva. La altura de 183.5 m es 1800.00 dividido entre 9.81. La de 177.2 m es 1737.98 dividido entre el mismo 9.81. Al apagar la fuente, la altura imprime 0.0 m con cualquier curvatura. El residuo de Coriolis a 30°, 2.2e-16 J/kg, no mueve ese 0.0 m. Los 290 g del arranque son una aceleración, no una potencia. $(-2\vec{\Omega}\times\vec{w})\cdot\vec{w} = 0$ sigue siendo una identidad con $\Omega = 300$ rad/s, aunque la aceleración llegue a 290 veces la gravedad. Un caudal de 30 kg/s por un canal que pasa de 20 mm a 8 mm solo aumenta $w_r$. La flecha se alarga. El trabajo permanece en 0.0e+00. En 50 mm, $U = \Omega r = 15.0$ m/s; en 150 mm es 45.0 m/s. El potencial centrífugo se cierra con esos dos radios.

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.

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. Paga el eje.

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. Si falta un término en el cambio de coordenadas, el libro se desvía. 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. A lo largo de una línea de corriente la constante es $I$, y $h_0 = I + U c_\theta$ no lo es. El pariente sin rotación, donde la presión total falla entre líneas de corriente, está en la constante de Bernoulli al cruzar una línea de corriente.

Una altura impresa de 0.0 m significa que a $h_0$ le falta la fuente#

Dentro de la zona rotante la constante es $I$. Al pasar a la zona estacionaria vuelve a ser $h_0$. Llevar $h_0$ con la velocidad relativa exige que vaya pegado $\Omega\,\partial(rc_\theta)/\partial s$. Al apagar esa fuente, 1800.00 J/kg pasa a 0.00 y la altura cae de 183.5 m a 0.0 m. Los residuos se ven limpios. Convergieron a cero. Refinar la malla deja el cero donde está.

Sumar Coriolis como fuente de energía es el arreglo equivocado. Su potencia imprime 0.0e+00, y a 30° sigue en $2.2\times10^{-16}$ J/kg. El término que corresponde ahí es la fuente de giro, no Coriolis. Él llegó a las coordenadas rotantes contando la eficiencia de ruedas hidráulicas. En el libro que armó, la fuerza con su nombre siempre es cero. Ciento noventa años después, dejar ese cero en cero es lo que mantiene el trabajo del eje en 1800 J/kg.