科里奥利力做的功是零，叶轮却加了1800 J/kg — 旋转坐标系里的能量与转焓

旋转坐标系中科里奥利力的功率严格为零。扬程仍然出现，是因为那个力正是叶片必须顶住的对象。

定义了"功"的人，他的力却不做功

把力学中的"功"与"动能"整理成今天形式的人是加斯帕尔·科里奥利。他1829年著作的标题是《论机器效果的计算》。那是一位工程师面对水车、清点投入与产出而写下的书。

然而以他命名的力并不做功。无论量级多大都一样。本文用一台离心叶轮来核对这件事。科里奥利加速度达到290 g，账本上却只写着 $10^{-16}$ J/kg，而流体的总焓上升了1800 J/kg。这1800从哪里来，以及旋转坐标系求解器漏掉对应项时会失去什么，就是全文的结论。

1832年，从水车计算走向旋转坐标系#

韩国的连载《历史中的流体力学》这样记述当年的巴黎。1830年七月革命把保王派的柯西赶了出去，共和派的纳维叶被任命到巴黎综合理工学院。随后在1832年，纳维叶开始与科里奥利合作研究。

科里奥利此前一直啃的题目是水车。在旋转的流体机械里，能量与功该怎么清点？要做这笔账，人就得坐到旋转的轴上去。旋转坐标系的研究由此而来，1835年的论文让今天的科里奥利力为世人所知。

顺序很重要。科里奥利力不是从地球自转或气象学里长出来的。它是为了对上旋转流体机械的能量账本才出现的。在下面的模拟中直接操作这本账。

omega [rad/s]300backsweep [deg]0

Flip the frame switch: the orange trail changes shape completely, yet the green rothalpy line never tilts. At the current radius 50 mm, I = 0.0 J/kg against h0 = 0.0 J/kg. Across the whole channel h0 climbs by 1800.0 J/kg — exactly U2·c_theta2 − U1·c_theta1 = 1800.0. Raise omega and only the orange line responds.

用 frame 按钮在相对系与绝对系之间切换视角。橙色轨迹会完全换一个形状，而右侧的绿线始终不倾斜。把 omega 滑块拉高，只有橙色的 $h_0$ 变陡。

垂直的力上不了账

在以角速度 $\vec{\Omega}$ 旋转的坐标系中，相对速度为 $\vec{w}$ 的流体质点满足

\frac{D\vec{w}}{Dt} = -\frac{\nabla p}{\rho} - 2\vec{\Omega}\times\vec{w} - \vec{\Omega}\times(\vec{\Omega}\times\vec{r})

右端第二项是科里奥利项，第三项是离心力项。两者都是把坐标系转起来所付的惯性代价。这两项如何进入动量方程，在旋转参考系与MRF中讲过。这里只看能量一侧。

要看能量，就把各个力与 $\vec{w}$ 点乘。科里奥利项当场死掉。

(-2\vec{\Omega}\times\vec{w})\cdot\vec{w} = 0

因为叉乘的结果与两个因子都垂直。这是恒等式，没有近似也没有条件。 $\Omega$ 取多少、 $\vec{w}$ 朝哪个方向，都是零。

离心力不同。它等于 $\Omega^2 r\,\hat{r}$ ，只要有径向速度，点乘就活着。作为交换，这一项拥有势函数。

\Omega^2 r\,\hat{r} = -\nabla\left(-\frac{U^2}{2}\right), \qquad U = \Omega r

有势函数意味着它可以挪到方程左端，折进常数里面。

所以守恒的是 $I$ ，不是 $h_0$ #

在定常、无粘条件下沿流线积分上式，得到的就是转焓(rothalpy，rotational + enthalpy)。

I = h + \frac{w^2}{2} - \frac{U^2}{2} = \text{const}

$h$ 是静焓， $w$ 是相对速度大小， $U=\Omega r$ 是叶片圆周速度。最后一项的负号就是离心势。科里奥利项因为点乘本来就是零，在这个式子里连痕迹都不留。

它与绝对系总焓 $h_0 = h + c^2/2$ 通过 $\vec{c} = \vec{w} + \vec{U}$ 相连：

h_0 = I + U c_\theta

其中 $c_\theta$ 是绝对速度的旋绕分量。也就是说 $I$ 保持常数的同时 $h_0$ 可以上升，上升量恰好是 $\Delta(U c_\theta)$ ，即欧拉透平机械方程。

用Python分别记三本账#

取一台半径从50 mm到150 mm、流道高度由20 mm收窄到8 mm的离心叶轮。 $\Omega = 300$ rad/s，水30 kg/s。后弯角 $\beta$ 看0°和30°两种。连续方程定下径向相对速度，转焓守恒定下压力。然后把科里奥利力、离心力、叶片反力所做的功分别对时间积分。

import numpy as np
 
OMEGA = 300.0            # 角速度 [rad/s]
RHO = 1000.0             # 密度 [kg/m^3]
MDOT = 30.0              # 质量流量 [kg/s]
R1, R2 = 0.05, 0.15      # 进出口半径 [m]
B1, B2 = 0.020, 0.008    # 进出口流道高度 [m]
GRID = np.linspace(R1, R2, 4001)
 
 
def channel_state(r, beta_deg):
    """半径 r 处的相对速度分量 (w_r, w_t)、叶片速度 U、绝对旋绕速度 c_t。"""
    b = B1 + (B2 - B1) * (r - R1) / (R2 - R1)
    w_r = MDOT / (RHO * 2.0 * np.pi * r * b)          # 连续方程决定的径向分量
    w_t = -w_r * np.tan(np.radians(beta_deg))         # 后弯叶片 → 与旋转反向
    u = OMEGA * r
    return w_r, w_t, u, u + w_t
 
 
def march_channel(beta_deg):
    """令转焓保持常数，沿流道建立 p/rho 与绝对总焓 h0。"""
    w_r, w_t, u, _ = channel_state(R1, beta_deg)
    i_const = 0.5 * (w_r**2 + w_t**2) - 0.5 * u**2    # 以 p1/rho = 0 为基准
    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):
    """沿轨迹分别对科里奥利力、离心力、叶片反力的功率作时间积分。"""
    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])            # 局部 (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   # 含柱坐标曲率项
        f_bld_t = acc_t - f_cor[1]                     # 剩下的旋绕加速归叶片压力面
        p_cor[k] = f_cor @ w
        p_cen[k] = f_cen @ w
        p_bld[k] = OMEGA * r * f_bld_t                 # 力矩 × 角速度 = 绝对系功率
    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):
    """直接积分 h0 的输运方程。源项为 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

转焓的漂移是 $10^{-13}$ J/kg，也就是双精度舍入的量级。科里奥利这本账在后弯角为0°时正好是 0.0，在30°时是 $2.2\times10^{-16}$ 。后者是浮点残渣，不是物理。

900与1800 — 另一半是谁付的#

有两个数字很显眼。离心力做的功是900 J/kg，可总焓上升了1800 J/kg。正好两倍。

\int \Omega^2 r\,w_r\,dt = \int_{r_1}^{r_2} \Omega^2 r\,dr = \frac{U_2^2 - U_1^2}{2}

这900是只在相对系内部流转的钱，分头进入压力与相对动能。它和绝对系中流体真正收到的1800是两个科目。

剩下的900由叶片支付。在径向叶片上，相对速度没有旋绕分量。要保持这一点，就必须有东西精确抵消科里奥利力 $2\Omega w_r$ ，那个东西就是叶片的压力面。作为反作用，流体受到同样大小的力，其力矩为 $2\Omega r w_r$ 。

分岔就在这里。在相对系中，那个力垂直于 $\vec{w}$ ，所以不做功。在绝对系中，叶片以 $\Omega$ 旋转，力矩乘角速度就是功率。积分得到 $2\times 900 = 1800$ J/kg，也就是输出中 work by blade torque 那一行。

归纳起来：科里奥利力一分钱也不给流体。它只是让流体去推叶片，轴功沿着那条反作用路径进来。力只决定方向，付账的是轴。

omega [rad/s]300backsweep [deg]0

The orange Coriolis arrow is the longest one on screen — up to 292 g — and its account still reads 0.0e+0 J/kg. Meanwhile the blade account fills to 0 of 1800 J/kg. Drag backsweep and only the purple bar moves: the perpendicular force never pays, it only decides where the blade has to push.

橙色的科里奥利箭头是画面上最长的，账本却牢牢贴在零上。拖动 backsweep 滑块，只有紫色的叶片科目在变。把 omega 调高，橙色柱条依旧不长。

当旋转坐标系求解器照搬输运 $h_0$ #

以上是物理。代码里出事的位置是固定的。

问题在于旋转网格区域内的能量方程用什么当输运变量。被相对速度携带、且无源守恒的是 $I$ 。若要用相对速度携带 $h_0$ ，就必须带上源项。

\vec{w}\cdot\nabla h_0 = \Omega\,\frac{\partial (r c_\theta)}{\partial s}

$s$ 是沿流线的坐标。漏掉这一项就得到输出的最后两行。有源项时1800 J/kg，扬程183.5 m；没有时0.00 J/kg，扬程0 m。流体穿过了叶轮，却什么都没发生。

症状之所以迷惑人，是因为残差看起来很正常。方程本身收敛得很干净，只是收敛到扬程为零的那个解。加密网格也没用，零还是零。

在OpenFOAM系求解器中，把MRF与可压缩能量方程一起用时，如果 rhothermo 一侧的总能量定义与MRF修正对不上，就是这个样子。就"坐标变换中漏掉一项"这一点而言，它与极坐标FVM的曲率项中看到的错误同根。上面代码在 acc_t 一行加入曲率项 $w_r w_\theta / r$ 也是同样的理由。把它去掉，30°后弯角下叶片账目就会偏到1802.06。

遇到旋转网格先问的三件事

第一，这个算例里保持常数的是 $h_0$ 还是 $I$ ？在旋转区域内部是 $I$ ，越过静止区域又变回 $h_0$ 。检查界面上到底衔接了什么。

第二，能量方程的源项是否真的接上了？扬程异常偏低、或者正好为零时，先看这一项。

第三，是不是把科里奥利项当能量源"加"了进去？那一项的功率严格为零。如果觉得那里似乎该放点什么，需要的不是科里奥利，而是 $\Omega\,\partial(rc_\theta)/\partial s$ 。

科里奥利是在清点水车效率的过程中走到旋转坐标系的。在他造出来的这本账里，冠他名字的力永远是零元。计算旋转的东西时把那个零保持为零，就是190年后我们接手这本账的方式。