インペラを止めたまま解く — 回転基準系と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 の圧力–速度連成と同じ場所を疑います。