Solving an Impeller by Holding It Still — Rotating Frames and MRF

A spinning impeller is solved while it sits still. The answer still comes out right. That is the usual trick for the steady performance of a pump or a fan. These notes follow the momentum equation in a rotating frame, then one free particle, until the sign of the source terms is checked.

Flip the frame toggle from absolute to relative (MRF) and the arrows inside the rotating zone fall from a rigid rotation at $\boldsymbol{\Omega}\times\mathbf{r}$ toward zero.

The outer pink field stays put. Only the absolute-frame arrows grow toward the interface. The relative frame has already moved that drag into the equation.

Drop the absolute frame and the windmill stops#

A camera bolted to the ground sees the impeller blades whirling. The geometry keeps changing in time. So the solve has to be unsteady.

Mount the camera on a blade and the blade sits still. Only the surrounding fluid flows past it. With the geometry frozen, a steady solve is possible. Turning a rotating problem into a stationary one is the whole idea of the rotating frame.

The price is not free. The frame itself accelerates, so two terms that were never there show up in the momentum equation.

Two fictitious forces sneak into the momentum equation#

Let the angular velocity be $\boldsymbol{\Omega}$, the relative velocity $\mathbf{u}_r$, and position $\mathbf{r}$. The link to absolute velocity is:

where $\mathbf{u}_a$ is the absolute (inertial) velocity and $\boldsymbol{\Omega}\times\mathbf{r}$ is the drag velocity from rotation.

Substituting into Navier–Stokes gives the steady momentum equation for the relative velocity:

The next term, $2\,\boldsymbol{\Omega}\times\mathbf{u}_r$, is the Coriolis acceleration. It bends the direction of motion. The term after that, $\boldsymbol{\Omega}\times(\boldsymbol{\Omega}\times\mathbf{r})$, is the centrifugal acceleration. It pushes outward from the axis. Neither is a real force. They appear only because the frame accelerates.

The Coriolis term enters as a source in the momentum equation. That is what OpenFOAM's MRFSource adds, cell by cell. The work done by Coriolis and the energy the impeller adds are not the same term. That split is written in rotating-frame energy and rothalpy.

MRF — freezing rotation into a frozen rotor#

MRF splits the domain in two: a rotating zone wrapping the rotor, and a stationary zone outside it.

• Rotating zone: solve the relative-velocity equation above. The Coriolis and centrifugal sources are on.
• Stationary zone: solve the ordinary absolute-frame equations.

The blade never actually moves. The mesh is fixed too. The rotation effect is only injected as a source. That is why the method is called a frozen rotor: the relative position of blade and fluid is frozen at one instant while the steady solution is found.

The toggle above is that freeze. In the absolute frame the inner fluid rotates as a solid body at $\boldsymbol{\Omega}\times\mathbf{r}$. The arrows grow toward the interface. Switch to the relative frame and the inner arrows collapse near zero, because the rotating-zone flow now looks stationary. The outer pink field is unchanged.

Stitching the velocity back together at the interface#

The two zones are solved in different frames. Velocities do not match at the boundary. The rotating zone stores $\mathbf{u}_r$, the stationary zone stores $\mathbf{u}_a$.

At the interface the transform is applied directly. Crossing from rotating to stationary, add $\mathbf{u}_a = \mathbf{u}_r + \boldsymbol{\Omega}\times\mathbf{r}$. The other way, subtract it. Scalars such as pressure and turbulence quantities pass straight through. Only vectors get the drag-velocity correction.

Skip this transform and mass and momentum start leaking at the interface. Residuals stall, or an unphysical jet crosses the interface. The common MRF mistake is placing the rotating-zone boundary too close to the blade. Put the interface where the flow is nearly axisymmetric, so the transform error stays small.

MRF versus sliding mesh — what, and when#

MRF is a steady approximation. It cannot capture the changing relative position between the blade and stationary parts such as a volute (rotor–stator interaction), because it freezes one instant.

When the unsteady effect is actually required, move to a sliding mesh. The mesh rotates and the interface is re-interpolated every step. Expensive, but the relative motion stays in the solution.

The usual path is short. Pin the operating point with MRF, then use that field to start a sliding-mesh run for the unsteady detail.

Python — the trajectory of a free particle in a rotating frame#

A free particle under no force checks the source terms. In the absolute frame it must fly in a straight line. Integrate only Coriolis and centrifugal terms in the rotating frame, transform back, and a restored straight line means the signs are right.

import numpy as np
 
OMEGA = np.array([0.0, 0.0, 1.2])   # angular velocity vector (rad/s), about z
 
def rotating_frame_acceleration(r, v_rel):
    """Apparent acceleration of a free particle (zero real force) in the
    rotating frame. Only the Coriolis and centrifugal terms survive."""
    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):
    """Integrate the relative-frame trajectory with a velocity-Verlet variant."""
    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)        # velocity-dependent -> one correction
        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):
    """Rotate the relative trajectory back to the absolute frame: by Omega*t each instant."""
    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)
 
# Choose initial conditions so the motion is a straight line in the absolute frame
V = 0.8
r0 = np.array([1.0, 0.0, 0.0])
v_abs0 = np.array([0.0, V, 0.0])                # absolute velocity: straight along +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()      # for a straight line, x stays ~ 1.0
print(f"max x deviation in absolute frame: {x_dev:.4e}")
print(f"absolute end point: {traj_abs[-1].round(3)}")

The initial data is hard-coded. $\boldsymbol{\Omega}=(0,0,1.2)$ rad/s and $\mathbf{r}_0=(1,0,0)$. The drag speed $|\boldsymbol{\Omega}\times\mathbf{r}|$ is $\Omega r=1.2$. Absolute speed is $V=0.8$ along $+y$, so $\mathbf{u}_a=(0,0.8,0)$. Relative velocity is the difference. Only the $y$ component changes: $0.8-1.2=-0.4$. The line v_rel0 is that subtraction.

At $t=0$ the Coriolis acceleration read from the equation has size $|2\boldsymbol{\Omega}\times\mathbf{u}_r|=2\times 1.2\times 0.4=0.96$, pointing in $-x$. Centrifugal acceleration has size $\Omega^2 r=1.44$, and the code's sign points outward, $+x$. Moving the two left-hand terms to the right-hand side matches the minus signs in the code.

The step is $dt=0.002$ and there are 1500 of them. The run covers 3 seconds. In that time the turntable rotates $\Omega t=1.2\times 3=3.6$ rad. If absolute $x$ stays near 1, the signs are right. Flip one sign and the particle flings outward. The $x$ departure grows.

The integrator is a velocity-Verlet variant. The position update carries $0.5\,a\,dt^2$, and because acceleration depends on velocity the velocity is corrected once. The transform-back angle at each instant is $wt$, with $w$ the $z$ component of $\Omega$, 1.2. One step is $1.2\times 0.002=0.0024$ rad. Fifteen hundred steps fill 3.6 rad. How far absolute $x$ leaves 1 is the maximum after that rotation is undone.

Raise the $\Omega$ slider and the gray line stays straight. Only the cyan path curls into a spiral.

At $\Omega=0$ the two overlap. At the code's 1.2 rad/s, cyan is what an observer sees after 3 seconds and 3.6 rad. Whether gray stays straight is the sign check.

Pitfalls to check before the code moves into a case#

Cross-check the axis and the direction of $\boldsymbol{\Omega}$ against the drawing. A flipped sign sends thrust the wrong way. The solve does not diverge. It returns a plausible wrong answer, which is harder to catch.

Place the interface where the flow is smooth. Slicing through the middle of a blade-tip vortex inflates the transform error. Make the rotating zone comfortably larger than the blade.

Do not leave an MRF field as the last answer. It is a steady approximation. If rotor–stator pulsation matters, check once more with a sliding mesh.

A flipped sign kills the straight line first#

The free particle has no force. Absolute $x$ has to stay at 1. Skip $\boldsymbol{\Omega}\times\mathbf{r}$ in v_rel0, or flip one of the Coriolis and centrifugal signs, and that line breaks first. If the $x$ departure grows over 3 seconds and 3.6 rad, suspect the source terms. The thrust curve goes wrong after that.

The same order carries into the impeller. Axis direction, then $\mathbf{u}_a=\mathbf{u}_r+\boldsymbol{\Omega}\times\mathbf{r}$ on the interface, then the performance curve. A sign that failed to restore the straight line, copied into the cell source, can drop the residual and still reverse the thrust. If that residual remains as a pressure checkerboard, the same spot shows up in SIMPLE pressure–velocity coupling.