Resolver un rodete dejándolo quieto — Marcos rotativos y MRF

Si no se resta la velocidad de arrastre 1.2, la x absoluta se aleja de 1. Omitir el mismo término en la interfaz deja escapar masa.

Se resuelve un rodete que gira dejándolo quieto. Y la respuesta sale bien. Es el truco más usado para el rendimiento estacionario de una bomba o un ventilador. Estas notas siguen la ecuación de cantidad de movimiento en un marco rotativo, y una partícula libre, hasta comprobar el signo de los términos fuente.

Al pasar el control de marco de absoluto a relativo (MRF), las flechas dentro de la zona rotativa caen de una rotación de sólido rígido a $\boldsymbol{\Omega}\times\mathbf{r}$ hacia cero.

omega1.00

Switch to the relative frame: the inner cyan arrows collapse toward zero because the rotating zone becomes steady, while the pink outer field is unchanged. The dashed circle is where the two velocity descriptions must be transformed into each other.

El campo exterior rosa no se mueve. Solo crecen, hacia la interfaz, las flechas del marco absoluto. El marco relativo ya metió ese arrastre dentro de la ecuación.

Si se abandona el marco absoluto, el molino se detiene#

Una cámara fijada al suelo ve las palas del rodete girar. La geometría cambia con el tiempo sin parar. Por eso hace falta un cálculo no estacionario (unsteady).

Con la cámara sobre una pala, la pala queda quieta. Solo el fluido alrededor circula de forma relativa. Con la geometría congelada, un cálculo estacionario (steady) se vuelve posible. Convertir un problema rotativo en uno estático es toda la idea del marco rotativo.

El precio no es gratis. Como el propio marco acelera, en la ecuación de cantidad de movimiento aparecen dos términos que antes no estaban.

Dos fuerzas ficticias se cuelan en la ecuación de cantidad de movimiento#

Sea la velocidad angular $\boldsymbol{\Omega}$ , la velocidad relativa $\mathbf{u}_r$ y la posición $\mathbf{r}$ . La relación con la velocidad absoluta es:

\mathbf{u}_a = \mathbf{u}_r + \boldsymbol{\Omega}\times\mathbf{r}

donde $\mathbf{u}_a$ es la velocidad absoluta (inercial) y $\boldsymbol{\Omega}\times\mathbf{r}$ es la velocidad de arrastre por la rotación.

Al sustituir en Navier–Stokes se obtiene la ecuación de cantidad de movimiento estacionaria para la velocidad relativa:

(\mathbf{u}_r\cdot\nabla)\mathbf{u}_r + 2\,\boldsymbol{\Omega}\times\mathbf{u}_r + \boldsymbol{\Omega}\times(\boldsymbol{\Omega}\times\mathbf{r}) = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u}_r

El término que sigue, $2\,\boldsymbol{\Omega}\times\mathbf{u}_r$ , es la aceleración de Coriolis. Desvía la dirección del movimiento. El término posterior, $\boldsymbol{\Omega}\times(\boldsymbol{\Omega}\times\mathbf{r})$ , es la aceleración centrífuga. Empuja hacia afuera desde el eje. Ninguna es una fuerza real. Aparecen solo porque el marco acelera.

El término de Coriolis entra como fuente en la ecuación de cantidad de movimiento. Eso es lo que el MRFSource de OpenFOAM añade, celda a celda. El trabajo de Coriolis y la energía que sube el rodete no son el mismo término. Ese corte está en energía del marco rotativo y rotalpía.

MRF — congelar la rotación en un frozen rotor#

MRF parte el dominio en dos: una zona rotativa que envuelve al rotor y una zona estacionaria afuera.

Zona rotativa: se resuelve la ecuación de velocidad relativa de arriba. Las fuentes de Coriolis y centrífuga quedan activas.
Zona estacionaria: se resuelven las ecuaciones ordinarias del marco absoluto.

La pala nunca se mueve de verdad. La malla también está fija. Solo se inyecta el efecto de rotación como una fuente. Por eso el método se llama frozen rotor (rotor congelado): la posición relativa entre pala y fluido queda congelada en un instante mientras se busca la solución estacionaria.

El control de arriba es ese congelamiento. En el marco absoluto el fluido interior gira como sólido rígido a $\boldsymbol{\Omega}\times\mathbf{r}$ . Las flechas crecen hacia la interfaz. Al cambiar al marco relativo, las flechas interiores se desploman cerca de cero, porque el flujo de la zona rotativa parece quieto. El campo exterior rosa no cambia.

Volver a coser la velocidad en la interfaz#

Las dos zonas se resuelven en marcos distintos. Las velocidades no coinciden en la frontera: la zona rotativa guarda $\mathbf{u}_r$ y la estacionaria guarda $\mathbf{u}_a$ .

En la interfaz se aplica la transformación directamente. Al cruzar de la zona rotativa a la estacionaria, se suma $\mathbf{u}_a = \mathbf{u}_r + \boldsymbol{\Omega}\times\mathbf{r}$ ; en sentido contrario, se resta. Los escalares como la presión y las cantidades de turbulencia pasan sin cambio. Solo los vectores reciben la corrección de la velocidad de arrastre.

Si falta esta transformación, la masa y la cantidad de movimiento empiezan a fugarse en la interfaz. Los residuos se estancan, o aparece un chorro no físico que cruza la interfaz. El error más común al configurar MRF es colocar la frontera de la zona rotativa demasiado cerca de la pala. La interfaz debe ir donde el flujo es casi axisimétrico, para que el error de transformación quede pequeño.

MRF frente a sliding mesh — qué y cuándo#

MRF es una aproximación de estado estacionario. No puede captar el cambio de posición relativa entre la pala y las partes estáticas como una voluta (interacción rotor–estator), porque congela un instante.

Cuando de verdad se necesitan efectos no estacionarios, se pasa a una sliding mesh (malla deslizante). La malla gira y la interfaz se vuelve a interpolar en cada paso. Costosa, pero el movimiento relativo queda en la solución.

Aspecto	MRF (frozen rotor)	Sliding mesh
Movimiento de malla	ninguno (fija)	rotación real
Tiempo	estacionario	no estacionario
Costo	bajo	alto
Interacción rotor–estator	no aproximable	captada
Uso	verificación de diseño, curvas de rendimiento	ruido, pulsación de presión

El camino habitual es corto. Se fija el punto de operación con MRF y esa solución arranca una corrida de sliding mesh para ver el detalle no estacionario.

Python — la trayectoria de una partícula libre en un marco rotativo#

Una partícula libre, sin fuerza alguna, comprueba los términos fuente. En el marco absoluto debe volar en línea recta. Si al integrar solo Coriolis y centrífugo en el marco rotativo, y transformar de vuelta, se recupera esa recta, los signos están bien.

import numpy as np
 
OMEGA = np.array([0.0, 0.0, 1.2])   # vector de velocidad angular (rad/s), eje z
 
def rotating_frame_acceleration(r, v_rel):
    """Aceleracion aparente de una particula libre (fuerza real cero) en el
    marco rotativo. Solo sobreviven los terminos de Coriolis y centrifugo."""
    coriolis = -2.0 * np.cross(OMEGA, v_rel)              # -2 Omega x u_r
    centrifugal = -np.cross(OMEGA, np.cross(OMEGA, r))    # -Omega x (Omega x r)
    return coriolis + centrifugal
 
def integrate_rotating_particle(r0, v0, dt, steps):
    """Integra la trayectoria del marco relativo con una variante velocity-Verlet."""
    r = np.array(r0, float); v = np.array(v0, float)
    traj = [r.copy()]
    a = rotating_frame_acceleration(r, v)
    for _ in range(steps):
        r = r + v * dt + 0.5 * a * dt * dt
        a_half = rotating_frame_acceleration(r, v)        # depende de la velocidad -> una correccion
        v = v + 0.5 * (a + a_half) * dt
        a = rotating_frame_acceleration(r, v)
        traj.append(r.copy())
    return np.array(traj)
 
def relative_to_absolute(traj_rel, dt):
    """Rota la trayectoria relativa de vuelta al marco absoluto: en Omega*t cada instante."""
    w = OMEGA[2]
    out = []
    for k, r in enumerate(traj_rel):
        t = k * dt
        c, s = np.cos(w * t), np.sin(w * t)
        out.append([c * r[0] - s * r[1], s * r[0] + c * r[1], r[2]])
    return np.array(out)
 
# Se eligen condiciones iniciales para que el movimiento sea recto en el marco absoluto
V = 0.8
r0 = np.array([1.0, 0.0, 0.0])
v_abs0 = np.array([0.0, V, 0.0])                # velocidad absoluta: recta en +y
v_rel0 = v_abs0 - np.cross(OMEGA, r0)           # u_r = u_a - Omega x r
 
dt = 0.002
traj_rel = integrate_rotating_particle(r0, v_rel0, dt, 1500)
traj_abs = relative_to_absolute(traj_rel, dt)
 
x_dev = np.abs(traj_abs[:, 0] - 1.0).max()      # si es recta, x se mantiene ~ 1.0
print(f"desviacion maxima de x en el marco absoluto: {x_dev:.4e}")
print(f"punto final absoluto: {traj_abs[-1].round(3)}")

Los datos iniciales están fijos en el código. $\boldsymbol{\Omega}=(0,0,1.2)$ rad/s y $\mathbf{r}_0=(1,0,0)$ . La rapidez de arrastre $|\boldsymbol{\Omega}\times\mathbf{r}|$ vale $\Omega r=1.2$ . La velocidad absoluta es $V=0.8$ en $+y$ , así que $\mathbf{u}_a=(0,0.8,0)$ . La velocidad relativa es la diferencia. Solo cambia la componente $y$ : $0.8-1.2=-0.4$ . La línea v_rel0 es esa resta.

En $t=0$ la aceleración de Coriolis leída de la ecuación tiene tamaño $|2\boldsymbol{\Omega}\times\mathbf{u}_r|=2\times 1.2\times 0.4=0.96$ , hacia $-x$ . La centrífuga tiene tamaño $\Omega^2 r=1.44$ , y el signo del código apunta hacia afuera, $+x$ . Pasar los dos términos del lado izquierdo al derecho coincide con los signos menos del código.

El paso es $dt=0.002$ y hay 1500. La integración cubre 3 segundos. En ese tiempo la plataforma gira $\Omega t=1.2\times 3=3.6$ rad. Si la $x$ absoluta se queda cerca de 1, los signos están bien. Invertir un signo lanza la partícula hacia afuera y la desviación de $x$ crece.

El integrador es una variante de velocity-Verlet. La actualización de posición lleva $0.5\,a\,dt^2$ , y como la aceleración depende de la velocidad, la velocidad se corrige una vez. El ángulo de vuelta en cada instante es $wt$ , con $w$ la componente $z$ de $\Omega$ , 1.2. Un paso vale $1.2\times 0.002=0.0024$ rad. Mil quinientos pasos llenan 3.6 rad. Cuánto se aleja la $x$ absoluta de 1 es el máximo después de deshacer esa rotación.

Al subir el control de $\Omega$ , la recta gris se queda quieta y solo la trayectoria cian se enrosca.

omega0.90

At omega = 0 both paths coincide. Increase omega and the cyan trail curls into a spiral — the Coriolis deflection that an observer riding the turntable mistakes for a force.

Con $\Omega=0$ ambas coinciden. En los 1.2 rad/s del código, el cian es lo que ve un observador tras 3 segundos y 3.6 rad. Que el gris siga recto es la prueba del signo.

Trampas antes de pasar el término al caso#

Conviene contrastar el eje y la dirección de $\boldsymbol{\Omega}$ con el plano una vez más. Un signo invertido manda el empuje al revés. El cálculo no diverge: da una respuesta verosímil pero equivocada, y por eso cuesta verla.

La interfaz va donde el flujo sea suave. Cortar por el centro de un vórtice de punta de pala infla el error de transformación. La zona rotativa tiene que quedar holgadamente más grande que la pala.

Una solución MRF no se deja como respuesta final. Es una aproximación estacionaria. Si importa la pulsación rotor–estator, se verifica una vez más con una sliding mesh.

Un signo invertido mata primero la recta#

La partícula libre no tiene fuerza. La $x$ absoluta tiene que quedarse en 1. Si en v_rel0 no se resta $\boldsymbol{\Omega}\times\mathbf{r}$ , o se invierte un signo entre Coriolis y centrífugo, esa recta se rompe primero. Si la desviación de $x$ crece a lo largo de 3 segundos y 3.6 rad, el sospechoso es el término fuente. La curva de empuje se equivoca después.

Al pasar al rodete el orden es el mismo. Dirección del eje, luego $\mathbf{u}_a=\mathbf{u}_r+\boldsymbol{\Omega}\times\mathbf{r}$ en la interfaz, y después la curva de rendimiento. Un signo que no restauró la recta, copiado a la fuente de la celda, puede bajar el residuo y aun así invertir el empuje. Si ese residuo queda como un tablero de presión, el mismo sitio aparece en el acoplamiento presión–velocidad de SIMPLE.