把叶轮固定住来求解 — 旋转参考系与MRF

不减去1.2的牵连速度，绝对轨迹的x就会离开1。界面上漏掉同一项，质量会漏出去。

把旋转的叶轮固定住再求解。答案仍然是对的。求泵或风机的稳态性能时，最常用的办法就是这个。旋转参考系里动量方程怎么变，MRF 怎么把旋转冻成稳态，这里一直追到一个自由粒子。

把坐标系开关从绝对拨到相对（MRF），旋转区内部的箭头会从 $\boldsymbol{\Omega}\times\mathbf{r}$ 的刚体旋转掉到接近零。

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.

外面的粉色场不动。只有绝对坐标系的箭头越靠近界面越大。相对坐标系已经把那份牵连放进方程里了。

抛弃绝对坐标系，风车就停了

固定在地面的相机看到叶轮叶片快速旋转。几何形状随时间不断变化。所以需要非稳态（unsteady）求解。

把相机装到叶片上，叶片是静止的。只有周围的流体相对流动。几何形状不变，稳态（steady）求解就成为可能。把旋转体问题变成静止问题，就是旋转参考系的全部想法。

代价不是免费的。坐标系本身在加速，动量方程里多出了原本没有的两项。

混进动量方程的两个假想力

设旋转角速度为 $\boldsymbol{\Omega}$ ，相对速度为 $\mathbf{u}_r$ ，位置为 $\mathbf{r}$ 。它与绝对速度的关系为：

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

其中 $\mathbf{u}_a$ 是绝对（惯性）坐标系速度， $\boldsymbol{\Omega}\times\mathbf{r}$ 是旋转带来的牵连速度。

代入 Navier–Stokes，得到关于相对速度的稳态动量方程：

(\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

从左边数的下一项 $2\,\boldsymbol{\Omega}\times\mathbf{u}_r$ 是科里奥利加速度，让运动方向偏转。再下一项 $\boldsymbol{\Omega}\times(\boldsymbol{\Omega}\times\mathbf{r})$ 是离心加速度，从轴向外推。两者都不是真实的力。它们只因坐标系加速而出现。

科里奥利项作为源项进入动量方程。这正是 OpenFOAM 的 MRFSource 逐单元加上去的内容。科里奥利做的功，和叶轮抬高的能量，不是同一个项。这个分叉写在旋转系能量与转焓里。

MRF — 把旋转冻结成 frozen rotor#

MRF 把区域分成两块：包住转子的旋转区（rotating zone），和它外面的静止区（stationary zone）。

旋转区：求解上面的相对速度方程。科里奥利、离心源项被打开。
静止区：求解普通的绝对坐标系方程。

叶片实际上并不移动。网格也是固定的。只是把旋转效应作为源项注入旋转区的方程。所以这种方式叫 frozen rotor（冻结转子）：在某一瞬间把叶片与周围流体的相对位置冻结，然后求稳态解。

上面拨过的开关就是这次冻结。绝对坐标系中，内部流体以 $\boldsymbol{\Omega}\times\mathbf{r}$ 做刚体旋转。箭头越靠近界面越大。切换到相对坐标系，内部箭头塌到接近零，因为旋转区的流动看起来静止了。外部粉色的场保持不变。

在界面上把速度重新缝合

两个区域用不同的坐标系求解。于是边界上的速度对不上。旋转区存的是 $\mathbf{u}_r$ ，静止区存的是 $\mathbf{u}_a$ 。

在界面上直接套用上面的变换。从旋转区跨到静止区时，按 $\mathbf{u}_a = \mathbf{u}_r + \boldsymbol{\Omega}\times\mathbf{r}$ 相加；反方向则相减。压力、湍流量这类标量原样通过。只有矢量接受牵连速度的修正。

漏掉这个变换，界面上的质量和动量就开始泄漏。残差降不下来，或者出现横穿界面的非物理射流。MRF 设置里最常见的错误，是把旋转区边界放得离叶片太近。应把界面放在流动近乎轴对称的位置，变换误差才小。

MRF 与 sliding mesh，何时用哪个#

MRF 是稳态近似。它无法捕捉叶片与蜗壳（volute）等静止结构相对位置随时间变化的效应（rotor–stator 相互作用），因为它把某一瞬间冻结了。

真正需要非稳态效应时，转向 sliding mesh（滑移网格）。网格真正旋转，每一步重新插值界面。代价高，但那份相对运动留在解里。

项目	MRF (frozen rotor)	Sliding mesh
网格运动	无（固定）	真实旋转
时间	稳态	非稳态
成本	低	高
rotor–stator 相互作用	无法近似	能捕捉
用途	设计校核、性能曲线	噪声、压力脉动

实务流程通常是这样：用 MRF 快速锁定工作点，再把该解作为初场跑 sliding mesh，看非稳态细节。

Python — 旋转坐标系中自由粒子的轨迹#

用一个不受任何力的自由粒子来核对源项。在绝对坐标系中它应当沿直线飞行。如果在旋转坐标系中只用科里奥利、离心项积分，再变回绝对坐标系时直线被复原，源项的符号就是对的。

import numpy as np
 
OMEGA = np.array([0.0, 0.0, 1.2])   # 旋转角速度矢量 (rad/s), 绕 z 轴
 
def rotating_frame_acceleration(r, v_rel):
    """自由粒子（真实力为零）在旋转坐标系中受到的视加速度。
    只剩科里奥利项和离心项。"""
    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):
    """用 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)        # 依赖速度 -> 修正一次
        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):
    """把相对坐标系轨迹变回绝对坐标系: 每个时刻反向旋转 Omega*t。"""
    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)
 
# 选取初始条件，使运动在绝对坐标系中是一条直线
V = 0.8
r0 = np.array([1.0, 0.0, 0.0])
v_abs0 = np.array([0.0, V, 0.0])                # 绝对速度: 沿 +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()      # 若为直线, x 保持 ~ 1.0
print(f"绝对坐标系 x 最大偏离: {x_dev:.4e}")
print(f"绝对坐标系终点: {traj_abs[-1].round(3)}")

初值写在代码里。 $\boldsymbol{\Omega}=(0,0,1.2)$ rad/s， $\mathbf{r}_0=(1,0,0)$ 。牵连速度 $|\boldsymbol{\Omega}\times\mathbf{r}|$ 的大小是 $\Omega r=1.2$ 。绝对速度 $V=0.8$ ，沿 $+y$ ，所以 $\mathbf{u}_a=(0,0.8,0)$ 。相对速度是两者之差。只有 $y$ 分量变成 $0.8-1.2=-0.4$ 。v_rel0 这一行就是这次减法。

$t=0$ 时，从方程读出的科里奥利加速度大小是 $|2\boldsymbol{\Omega}\times\mathbf{u}_r|=2\times 1.2\times 0.4=0.96$ ，方向为 $-x$ 。离心加速度大小是 $\Omega^2 r=1.44$ ，代码的符号指向轴外的 $+x$ 。把动量方程左边那两项挪到右边，就和代码里的负号对齐。

时间步长 $dt=0.002$ ，步数 1500。积分时间是 3 秒。这段时间里转盘转过 $\Omega t=1.2\times 3=3.6$ rad。绝对轨迹的 $x$ 留在 1 附近，就说明这个符号是对的。哪怕翻转一个符号，粒子都会被甩向外侧， $x$ 偏离变大。

积分是 velocity-Verlet 的一种变体。位置更新里有 $0.5\,a\,dt^2$ ，加速度依赖速度，所以速度要再修正一次。变回去的转角是每个时刻的 $wt$ ， $w$ 是 $\Omega$ 的 $z$ 分量 1.2。一步的角度是 $1.2\times 0.002=0.0024$ rad。1500 步填满 3.6 rad。绝对坐标的 $x$ 离开 1 多远，看的是把这个旋转转回去之后的最大值。

把 $\Omega$ 滑块往上推，灰色直线不动，只有青色轨迹卷成螺旋。

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.

$\Omega=0$ 时两者重合。滑块来到代码里的 1.2 rad/s，青色就是转了 3 秒、3.6 rad 的观测者看到的弯曲。灰色还直不直，就是源项符号的判定。

落到代码之前要盯住的陷阱

把旋转轴和 $\boldsymbol{\Omega}$ 的方向再跟图纸核对一遍。符号反了，推力就反向。它不发散，给出一个貌似合理却错误的答案，因此很难抓住。

把界面放在流动平滑的地方。从叶尖涡正中切过会放大变换误差。旋转区要比叶片留得足够大。

不要把 MRF 解当成最后的答案。它是稳态近似。若 rotor–stator 脉动很重要，就再用 sliding mesh 核对一次。

符号一反，直线先死

自由粒子没有力。绝对坐标的 $x$ 应该停在 1。v_rel0 里不减 $\boldsymbol{\Omega}\times\mathbf{r}$ ，或者科里奥利和离心的符号翻转一个，这条直线会先垮。3 秒、3.6 rad 里 $x$ 偏离变大，就该怀疑源项。推力曲线是在那之后才错的。

搬到叶轮计算时顺序一样。先是轴的方向，然后是界面上的 $\mathbf{u}_a=\mathbf{u}_r+\boldsymbol{\Omega}\times\mathbf{r}$ ，再是性能曲线。没能把直线还原的符号如果写进单元源项，残差可以下降，推力却是反的。那份残差若留成压力的棋盘格，就要怀疑和 SIMPLE 压力–速度耦合同一处。