OSSローターダイナミクス解析ソフト【ROSS】の紹介

Contents

はじめに

今回はオープンソースのローターダイナミクス解析ソフト “ROSS”について紹介します。

ターボ機械シリーズ第三回、RDの続きで公式チュートリアルの和訳解説になります。使い方を覚えれば何かと便利なので、是非トライしてみてください！

ローターダイナミクスの基本については前回記事を読んでみて下さいね！

この記事はこういうことを知りたい人向け

ローターダイナミクスへの理解度をもう少し深めたい

ローターダイナミクス解析を体験してみたい

今回の記事のレーティング：

ローターダイナミクス解析

商用ソフト

RD解析は手計算でやるには大きなボリュームになってしまうため、解析ソフトを使うのが一般的です。基本的には有償の商用ソフトで、FEAのような感じのものをイメージしてもらえればと思います。

専用の計算ソフトとしては、こういうのがあります。多くは一般的な設計ツールのアドインとして提供されていたりするから、中々個人では手が出せないのが悩みですね。

RD解析ソフト

Dyrobes: 独立系専用ソフトで米国ターボ業界で広く利用されている。比較的歴史が新しく、積極的に開発がされている。
XLrotor: 同様の主要商用ソフト
COMSOL Rotordynamics: COMSOL用アドインで汎用性が高い
DYNRot：Matlab用ソフトで機能性が高い
AMrotor:Matlab用アドインでベーシックな機能に留まる。今回紹介のROSSに近い
MyRot:日本RD研究の始祖、防大松下教授が開発した純国産初にしてほぼ唯一のVBAアプリケーション。歴史が長い。
XLTRC2: テキサスA&M大学が開発しているVBAアプリケーション
Ansys RD: Ansys用アドイン、三次元計算がメイン
MSC RD: Nastran用アドイン、三次元計算がメイン

計算に当たって必要な試験

計算ソフトへの入力諸元としては、軸やローターの慣性、剛性などがありますが、実測値が特に重要となってくるのは軸受やシールの特性です。
軸受の剛性が既知でない場合、単体試験でデータを取得しソフトへの入力値として使います。計算を回す前に必要なステップは以上です。

他の要素に関しも同じです。

磁気浮上軸受や動的な制御を行う軸受の場合は更に試験が複雑化してきてしまい、従来ではまかり通ってしまった経験則や”感じニアリング”では通用しなくなってきているのが現状です。
また一般的なボールベアリングでも”ガタ”影響などにより軸受は非線形な特性を示すことがほとんどです。

ROSSの使い方

それでは早速ソフトを導入してサンプルの計算を回してみましょう。詳細については公式HPに書かれています:

https://ross.readthedocs.io/en/stable/

インストール

動作環境として、python 3.7以上が必要です。ライブラリとして提供されているため、下記よりインストールが可能です。その他、実行に必要な各ライブラリ(numpyやpyplotなど)は逐次入れてください。

pip install ross-rotordynamics

手動で入れたい場合は、下記から各verがdl可能です。

ROSS is a library written in Python for rotordynamic analysis.
https://github.com/petrobras/ross
105 forks.
136 stars.
25 open issues.

Recent commits:

Release v2.0.0rc2, raphaeltimbo
Merge pull request #1132 from jguarato/fix/bearing-elementReturn kwargs in arguments of BearingElements, GitHub
Return kwargs in init arguments in BearingElements, jguarato
Release v2.0.0rc1, raphaeltimbo
Merge pull request #1127 from jguarato/remove/4dof-modelRemove 4 DoF (leave only 6 DoF) model, GitHub

ROSSには公式チュートリアルも提供されており、下記から履修可能です。

https://ross.readthedocs.io/en/stable/tutorials/tutorial_part_1.html

全部英語なので分かりづらい人もいるかも。。。と思ってこの記事を書きました！中身は大体同じです。

計算の流れ

ROSSは基本原則として各要素のコードブロックを記述し、pythonの実行順序に従って上から順に計算が進みます。大まかな計算の流れは下記の通りです。

チュートリアル

公式チュートリアルは全文英語、かつ実行するにあたって誤解を生みやすいような点もあるため、今回は管理人が用意したカスタマイズ版を使用してください。下記リポジトリを落としたら、ROSSの中身を使っていきます。

Turbo Compressor Design Tools
https://github.com/Gertrud-Violett/TurboTools
1 forks.
3 stars.
0 open issues.

Recent commits:

Updated shaft, NHI_Owner
added general plot MD, Trude
update comp map plot, Trude
fix range bug. added levels to Comp Map Plot, Trude
Released v0.1 Turbine Mapping Tool, Trude

記事内はコードの抜粋のみのため、直接dlしたファイルを参照するのがおすすめです。

公式を読みながら自分で進めたい方は一気に次の章に進んじゃって大丈夫です。

ライブラリの導入

実行に必要なライブラリをインポートします、インストールされていない場合はpipからinstallをしてください。

import os
from pathlib import Path
import ross as rs
import numpy as np
import ross.bearing_seal_element
import plotly.graph_objects as go

動作確認

ROSSをインストールしたら、TurboTools内部のROSS/Tutorialに動作に必要なすべてのコード一式が入っています。まずはこれで動くか、確認をしてみましょう。なお、わかりやすさのため要素毎に.pyを分けていますが、一つのファイルに書いてしまっても問題なく動作します。

ターミナルから下記で動作確認をしてください。(IDE上でなく、ターミナル直打ち推奨です。)

> python Materials.py
> python Rotor.py

この段階でoutputフォルダの中にrotor_pic.htmlという図が保存されているはずです。html iframe画像のためズームイン・アウトや各種数値の確認ができます。

次に本体の実行確認をします。計算の実行に数分程度かかるため、気長に待ちましょう。計算が成功すれば諸元の出力と共に、各種画像がoutputフォルダに生成されているはずです。

> python MAIN.py
Undamped natural frequencies[rad/s, rpm]:
 [ 144.72997813  205.33566432  386.64342244  444.34294291 1098.8108737
 1559.12921705], [ 1382.06948593  1960.81115821  3692.1727137   4243.16254755
 10492.87092432 14888.58730879]

Damped natural frequencies[rad/s, rpm]:
 [ 144.24521631  204.6256489   386.35509925  444.01366549 1098.81073232
 1559.12896433], [ 1377.44035159  1954.03101033  3689.41943006  4240.01817978
 10492.86957426 14888.58489548]

Damping ratio for each mode:
 [0.08177791 0.0830885  0.03861166 0.03849075 0.00050728 0.00056937]

Log Decrement for each mode:
 [0.51555254 0.52387193 0.24278527 0.24202386 0.00318731 0.00357746]

Whirl Direction for each mode:
 [1.  0.5 1.  0.  1.  0. ]
Extrapolating bearing coefficients. Be careful when post-processing the results.

最後のエラー”Extrapolating bearing coefficients”は、軸受剛性の定義が甘いのでプログラム側で適当に補完しちゃうよ、という警告です。動作には問題ありませんが結果がおかしいときは定義の見直しをしみましょう。

ローターの構築

材料の定義

軸系に使用される材料の密度、縦弾性係数、横弾性係数あるいはポアソン比を設定します。ここでは軟鋼での定義の例3通り、クロモリ、アルミ合金といった形で指定しています。

Material.py
#Sample Steel
# from E and G_s
steel = rs.Material(name="Steel", rho=7810, E=211e9, G_s=81.2e9)
steel.save_material()
# from E and Poisson
steel2 = rs.Material(name="Steel", rho=7810, E=211e9, Poisson=0.3)
# from G_s and Poisson
steel3 = rs.Material(name="Steel", rho=7810, G_s=81.2e9, Poisson=0.3)

ChroMoly1 = rs.Material(name="CrMo4340", rho=7900, E=212e9, G_s=82.2e9)
ChroMoly1.save_material()


#Aluminum Alloy
AC1 = rs.Material(name="A2618", rho=2780, E=78e9, Poisson=0.27)
AC1.save_material()

注意点！

注意点!：材料を定義すると、Materials.tomlというデータベースファイルがrossインストールディレクトリに保存されます。誤った材料データを使わないために、必ず下記のことを心がけてください！

“マテリアルの名前を確認する”
“MAIN計算実行前にMaterial.pyを再実行する (古いデータが上書きされる)”

他の場所に変更したい方はross本体のmaterialメソッドを書き換えてください

軸の定義

次に軸の定義です。軸は各ノードにおいて、長さ、内径、外径、材質及び各種オプションを定義します。

Shaft.py
n = 10
shaft_elem = [
    rs.ShaftElement(
        L=0.05,
        idl=0.0,
        odl=0.1,
        material=steel,
        shear_effects=True,
        rotary_inertia=True,
        gyroscopic=True,
    )
    for _ in range(n)
]

n = 10
shaft_elem2 = [
    rs.ShaftElement(
        L=0.05,
        idl=0.0,
        odl=0.12,
        material=steel,
        shear_effects=True,
        rotary_inertia=True,
        gyroscopic=True,
    )
    for _ in range(n)
]

#Collar
shaft_elem3 = rs.ShaftElement(
        L=0.05,
        idl=0.12,
        odl=0.16,
        material=steel,
        shear_effects=True,
        rotary_inertia=True,
        gyroscopic=True,
        n=18
    )

shaft_elem4 = rs.ShaftElement(
        L=0.05,
        idl=0.12,
        odl=0.16,
        material=steel,
        shear_effects=True,
        rotary_inertia=True,
        gyroscopic=True,
        n=19
    )

shafts=[shaft_elem,shaft_elem2,shaft_elem3,shaft_elem4]

ここでは全く同じ径の軸要素が10ノード連続で並んでいるため、for文で同じものを複数生成しました。生成した要素は最後に順番を揃えてひとつのlistに格納する必要があります。

サンプルコードにはいくつかの定義のためのオプションがコメントアウトされているので、色々と試してみてどのような結果になるか見てみるとわかりやすいかと思います。

Diskの定義

インペラやタービンなど、円盤状の回転部品をDiskを呼んでいます。Shaft同様に内径、外径、長さ等を指定します。直接慣性モーメントや質量として定義することも可能で、後述します。

n= はShaftのどの位置に配置されるかを示したもので、n=0からみた累計ノード数です。

Disk.py
#DISK=======================================================================
disk0 = rs.DiskElement.from_geometry(
    n=1, material=steel, width=0.20, i_d=0.10, o_d=0.23
)
disk1 = rs.DiskElement.from_geometry(
    n=20, material=steel, width=0.05, i_d=0.12, o_d=0.28
)
disk2 = rs.DiskElement.from_geometry(
    n=10, material=steel, width=0.3, i_d=0.12, o_d=0.20
)

#Create Disk Array==========================================================
disks = [disk0, disk1, disk2]

ROSSの基本的な使い方は以上です。このような感じで他の要素も定義をして行きます。

軸受の定義

Bearing_Seal.py
#BEARING DEFINITIONS========================================================
stfx = 1e6
stfy = 0.8e6
#bearing0 = rs.BearingElement(n=3, kxx=stfx, kyy=stfy, cxx=1e2, n_link=None)

bearing0 = rs.BearingElement(
    n=3, 
    kxx=[1.0e6, 2.0e6, 3.0e6],
    kyy=[0.8e6, 1.8e6, 2.8e6],
    cxx=[1.0e3, 1.5e3, 2.0e3],
    frequency=[0, 1000, 2000],
    n_link=None
)

bearing1 = rs.BearingElement(
    n=16, 
    kxx=[5.0e6, 10.0e6, 15.0e6],
    kyy=[5.0e6, 10.0e6, 15.0e6],
    cxx=[1.0e3, 1.5e3, 2.0e3],
    frequency=[0, 1000, 2000],
)
# Constant Coefficient Seals================================================
stfx = 1e2
stfy = 0.8e2
seal1 = rs.SealElement(n=2, kxx=stfx, kyy=stfy, cxx=1e3, cyy=0.8e3)
seal2 = rs.SealElement(n=17, kxx=stfx, kyy=stfy, cxx=1e3, cyy=0.8e3)
seal3 = rs.SealElement(n=9, kxx=stfx, kyy=stfy, cxx=1e1, cyy=0.8e1)


bearings = [bearing0, bearing1, seal1, seal2, seal3]

軸受及びシールは剛性kxx, kyy及び減衰係数 cxx (cyy) 等を指定します。ここでは回転速度毎に定義していますが、他にもいくつかのモードがあります。

n_link = はラジアル方向の錬成要素かどうか、の指定オプションで、多軸回転機械やダンパ等の定義に用います。のちど説明をします。

軸受特性の確認

入力した軸受の剛性プロットを手早く出力することができます。

#Check Bearing specs========================================================
bearing_fig=bearing0.kxx.plot() #Enter any bearing no and attribute
#bearing_fig.show()
bearing_fig.write_html("./output/bearing_pic.html")

PointMass

質点の定義です。体積ゼロの質量m (kg)をノードnに追加します。例えばロックナットやバランスウェイト等を再現するのに使います。

PointMass.py
pm0 = rs.PointMass(n=20, m=0.1)
pointmass = [pm0]

ローター構築のための設定は以上です！

Rotor.py

import os
from pathlib import Path
import ross as rs
import numpy as np
import Materials
import Shaft
import Disk
import Bearing_Seal
import PointMass

shafts=Shaft.shafts
disks=Disk.disks
bearings=Bearing_Seal.bearings
pointmass=PointMass.pointmass

#ROTOR DEFINITIONS==========================================================
rotor = rs.Rotor(shafts, disks, bearings)
#rotor = rs.Rotor(shafts, disks, bearings, pointmass) #for n_link option

# plotting the rotor model==================================================
rotorfig = rotor.plot_rotor()
#rotorfig.show()
rotorfig.write_html("./output/rotor_pic.html")

#Save & Load ROTOR==========================================================
rotor.save('./output/rotor.toml')

設定した各要素をimportし、rotorインスタンスを作ります。このようにすべての要素を設定した上でRotor.pyを実行すると、outputフォルダにこんな感じのrotor.tomlというファイルが生成されます。

[parameters]

["ShaftElement_ShaftElement 0"]
L = 0.05
idl = 0.0
odl = 0.1
idr = 0.0
odr = 0.1
n = 0
axial_force = 0
torque = 0
shear_effects = true
rotary_inertia = true
gyroscopic = true
shear_method_calc = "cowper"
tag = "ShaftElement 0"

["ShaftElement_ShaftElement 1"]
L = 0.05
idl = 0.0
odl = 0.1
idr = 0.0
odr = 0.1
n = 1
axial_force = 0
torque = 0
shear_effects = true
rotary_inertia = true
gyroscopic = true
shear_method_calc = "cowper"
tag = "ShaftElement 1"

これが構築されたローターの中身です。tomlファイルはいちいちpythonで作る必要はなく、csvやjsonから読み込んで変換しても構いません。最終的にここで記述された形式でRotor.pyに読み込みができれば計算できます。

また、それぞれの作った要素を個別にsaveして後で別の計算に利用することも可能です。

csvで直接書いたほうがわかりやすいという方向けに、変換用ツールのcsv2toml.pyも同梱してあるのでご活用ください。ROSSの機能としてエクセルからの読み込みも可能で、後述します。

計算の実行

MAIN.pyから計算します。ターミナルに表示されるそれぞれの出力内容は下記の通りです。

Undamped natural frequencies[rad/s, rpm]: 非減衰固有振動数r (rad/s 及び rpm)
Damped natural frequencies[rad/s, rpm]: 減衰固有振動数
Damping ratio for each mode:各モードの減衰比
Log Decrement for each mode: 対数減衰率δ
Whirl Direction for each mode:　各モードのホワール向き (1:前回り/ 0:後ろ回り)

MAINの各項目についてこのように記述してます：

1.1 Static analysis: 静的梁としての解析結果
1.2 modeplot: 2d及び3dのモードシェイプの描画します。
1.3 campbell: キャンベル線図
1.4 freq_response: 周波数応答線図
1.5 unbalance_response: 不釣り合いを設定した際の応答 (ampsにて設定)。plotする際には、計算済のspeedをrad/sで指定する必要があります。(linspaceの関係上、離散的にしか結果が出力されない)
1.6 time_response: 強制外力を与えられた際の応答→振動等の外乱を想定。非定常状態のいて系が収束するかどうかの確認に使います。
1.7 ucs:非減衰危険速度線図。軸受剛性を設定して描画します。(計算に時間がかかります)

これで基本的な計算は一通り完了しました。
それぞれのプロットがoutputに出力いるので確認してみましょう。(各図は拡大縮小や回転ができます)

ここで挙げたチュートリアルサンプルでは、シール特性の試験機を想定しています。モード4以上は軸が変形していて弾性軸となっている様子が良く分かりますね。プロット結果から軸の動きをイメージしてみましょう！

ROSSが使えるようになりました！次の章ではもう少し実用的なターボ機械を想定して、参考事例を元に各機能を紹介します。

1 2