コリオリ力の仕事はゼロなのに羽根車は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}$ がどの向きでもゼロです。

遠心力は違います。 $-\vec{\Omega}\times(\vec{\Omega}\times\vec{r}) = \Omega^2 r\,\hat{r}$ であり、半径方向の速度があれば内積は生き残ります。その代わりこの項はポテンシャルをもちます。

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

ポテンシャルがあるということは、方程式の左辺へ移して定数のなかに畳み込めるということです。

だから保存されるのは $h_0$ ではなく $I$ です#

定常・非粘性の条件で上の式を流線に沿って積分すると出てくるのがロータルピー(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$ です。

コリオリは水車の効率を数えているうちに回転座標系にたどり着きました。彼が作った帳簿で、彼の名の力はいつも0円です。回転するものを計算するときにそのゼロをゼロのまま置くことが、190年後に私たちがその帳簿を引き継ぐやり方です。