プロペラ・風車の設計解析ソフトXROTORの使い方

[[]]
人力飛行機などのプロペラや風力発電用の風車の解析や設計をやりたい人向けの記事です。
ブレードが回転するものの設計ならかなり汎用的に使えるはずです。
CUIソフトなので少し取っ付きにくいですが、慣れれば短時間で細かく設計・解析ができます。
解説を書いていますが、関連論文を全部精読しているわけではないので間違っている所が多くあるかもしれません。

プロペラ解析・設計ソフトの概要

ソフトの製作者はMITのMark Drela教授。人力飛行機の長距離飛行の世界記録を持っているDaedalusの機体設計のために開発されたようです。
公式HPのユーザーガイドによると以下の機能があります。

  • 誘導損失最小の回転翼の設計(プロペラや風車)
  • 任意の回転翼の形状の入力
  • 回転翼の形のインタラクティブな変更
  • 任意の回転翼の誘導損失最小となるねじり角の最適化
  • 回転翼の解析
  • プロペラに入って来る流れの影響の考慮
  • 構造解析と荷重がかかった状態のねじり角の解析
  • 音響解析をした騒音範囲予想
  • 幾何データや空力パラメータ、性能マップのプロット

回転翼の空力的な設計と解析に加えて、構造解析や音響解析もできる機能があります。
プロペラや風車の適応範囲はとても広いです。飛行機や飛行船、ホバークラフトのプロペラ、扇風機の羽、PCファン、揚力型の風力発電のタービンの設計まで可能です。ヘリコプターや水中プロペラもできるようですが、元々空気中かつ低Reの航空機のプロペラの設計ソフトとして開発されているので精度は不明です。
空気中であっても解析できないものとしては2重反転プロペラや回転面から傾斜角がついたプロペラがあります。

設計手法

プロペラの設計手法として大きく分けて二種類あります。

  • Betzの条件を用いて適切な近似を行うことで設計を行う手法

この中でよく使われる方法はLarrabeeの設計法というものです。このLarrabeeからもう少し近似を行ったAdkins・Liebeckの設計法という方法は以前に人力飛行機のプロペラの設計解析できるようなエクセルファイルを作ったりしましたが、近似式の変形がアヤシく、条件によっておおきく誤差が生まれる理論になっています。

  • 渦法と呼ばれる手法

プロペラが回ることで後流が誘導されるのをビオ・サバールの法則を用いて計算して性能計算を行い、プロペラ『効率』の最適化問題を解くことで設計する方法です。前述のBetzの条件の手法ではプロペラ『損失』最小の最適化問題を解くので何を最適にするかの考え方が違います。

XROTORでは改善されたLarrabeeの設計法を用いています。元となるLarrabeeの設計法の弱点としてレイノルズ数が小さい領域(Re<300,000)では誘導エネルギ損失だけではなく形状抗力の損失が大きくなるのに対応していない問題があります。これをXROTORではレイノルズ数による形状抗力の変化のモデルを取り入れることによって改善しています。

プロペラ設計・解析ソフト

ダウンロード・インストール

XROTOR Download Page
からソースコードが配布されているのでこれをWindowsLinuxでビルドすると使えます。
Linux、特にUbuntuで使うためには過去の記事を参考にしてください。非常に簡単です。
XROTORのUbuntuへの簡単インストール方法
XROTORがフリーソフトになってたのでUbuntu上で動かしてみた

実行

Windows環境で.exeを作った場合は別だが、LinuxのターミナルやCygwinなどから起動する方法。

xrotor

もしくは

 xrotor save_file case_file

で起動する。ここでsave_fileは以前に作ったプロペラの形状のファイル。case_fileは解析ファイル。これらの引数を書くと現在のディレクトリ(フォルダ)と同じ場所に保存されてたファイルの中身が読み込まれる。save_fileを書くとLOADコマンドをしたのと同じ。case_fileはOPERメニューでCGETコマンドしたのと同じ。

XROTORの全体の構成としてはトップレベルメニューがあって、そこからサブメニューにわかれている。(図)
一般的な使い方としては

 XROTOR   c>   desi
.DESI   c>  <Return>
 XROTOR   c> 

とするとサブメニューに飛んでからトップレベルメニューに戻ってくる。
?コマンドを入れるとそのメニュー中で使用できるコマンドの一覧が表示される。

使い方


起動するとトップレベルメニューに入れるのでそこから各サブメニューに移動して設計・解析を行う。
DESIサブメニューで設計してOPERサブメニューで解析を行なってPLOTコマンドで解析結果を見るのが基本的な流れになる。
最後に各サブメニューのコマンド内容の日本語訳を載せておきます。

AEROサブメニュー

Larrabeeの設計法の弱点をなくすためにRe数と翼型の性能を指定する必要がある。それがAEROサブメニュー。
もし一つのaero sectionしか存在しないと全てがその条件ということになる。もし複数あれば線形補間されることになる。




f = -0.1 to -0.2 ;高Re乱流 Re > 2M
f = -0.5 to -1.5 ;低Re領域 Re ~ 200K..800K
f = -0.3 to -0.5 ;ほぼ層流 Re < 100K

低Re領域においてはfが強く影響している。
失速前と失速後の2つのモデルを簡単に合わせたものを翼型特性のモデルとしている。
この扱いは地上の2次元翼ではそこそこ正しいが、実際のローターは強い遠心力とコリオリ力によって境界層流れを剥がそうとするのでかなり異なる。
Mcritは圧縮性の効果に効いてくるファクター。
\delta C_D = K (Mach - M_{crit})^2
圧縮性が効いてくるときは試験されたデータをもってモデルを確認する必要がある。

技術的なこと

解析における手法の違い

FORMコマンドによって以下に述べる誘導速度と誘導損失の計算手法を選択できる。
両手法はローターのブレードを揚力線として扱う。共に円盤荷重が比較的小さく、したがって後流の収縮や自己変形が小さいと仮定している。

  • Granded Momentum Formulation

プロペラにおいて古典的手法でE.E Larabeeが復興さえたもの。Betz-Prandtの先端損失誤差に頼っている。これはローターが低進行率もしくは多くのブレードがある場合を仮定している。結果としてこの手法は進行率0.5以上では適応できないとある。主なアドバンテージは極めて計算量が少ないこと。

  • Potential Formulation

これは少し現代的なアプローチである。リジッドならせん状の後流について螺旋対称のポテンシャル流を解くものであり、これはブレードの数や進行率が任意で適応できる。これはGoldsteinの2つと4つのブレードの計算をすべてのブレード数と任意の円盤荷重に拡張したものである。デメリットとしてはこのアプローチはGraded Momentumより計算量が多くなる。しかし現在の計算機では気にするほどではない。この手法は微小な半径のローターハブと取り付けられたナセルも取り扱うことができる。これはBetz-Prandtlの誤差では適切に取り扱えないものである。この手法はまたローターの先端の境界条件(ダクトの有無)を変えることが出来る。これはTOP LEVELでDUCTコマンドで有効無効にできる。

  • Vortex Formulation

渦法。これは誘導速度計算による離散化された渦後流を使うもの。離散化された直線渦(無限遠下流側の揚力線から得られたリジッドな螺旋状の後流)によってプロペラ揚力線上の速度が計算される。このオプションは傾斜や後退角の付いたブレードのプロペラのためのもの。この渦法はポテンシャル流の手法より少し計算時間がかかるより色々なプロペラ形状を扱うことが出来る。後流が自己変形しないことに注意。渦法はOPERかDESIメニューのVRTXコマンドで有効になる。理論的には非平面のプロペラも解析可能だが、XROTORではそのような幾何形状の入力をサポートしていない。

JMAPコマンド

このコマンドで作られるデータはCp vs Jの等ブレード角と等効率の等高線としての2次元プロットである。プロペラのパフォーマンスの要約であり、アメリカのプロペラメーカHamilton Standardなどによって使われている。このコマンドはポテンシャル流の手法を使っているときならいつでも実行することができる。ギザギザの等高線の生データはかなり早く生成できるがより高品質なプロットを求めるために実行時間が必要である。JMAPコマンドで生成されるデータはテキストエディタなどで読むことができない。Ubuntuなどにはデフォルトで入っている"JPlot"というソフトでグラフにすることが出来る。

jplot
Enter dump filename s> (dumpファイル名)
go
go
y

のコマンドでグラフ化できる。

単位と無次元化

XROTORではSI単位系が使われている。しかし内部の計算は全て無次元化してから計算しているので全てを別の単位系でインプットすればその単位系で出力される。例えばmm単位でインプットして計算すると出力はm単位のように出力されるが、実際はmm単位の結果になっている。
T_c = \frac{T}{\frac{1}{2} \rho V^2 \pi R^2} 、 P_c = \frac{P}{\frac{1}{2} \rho V^3 \pi R^2}
C_t = \frac{T}{\rho D^4 n^2} 、 C_p = \frac{P}{\rho D^5 n^3}
J = \frac{V}{n D} = \pi \frac{V}{W R} 、:W = w\pi n[rad/s]
\eta = \frac{T_c}{P_c} = J \frac{Ct}{Cp} 、 G = B \frac{\frac{\Gamma}{VR}}{2\pi L_w}
\eta _i L_w = \frac{V}{WR}
nは回転数/秒[rps]、Gammaは循環、Lwは後流進行率。eta_iは非粘性の効率、V\WRは幾何的な進行率。
Gは理論的なペラ効率の上限にどれだけ近づいているかの指標で1.000に近いほど良い。

風車の場合

XROTORの入力と出力はプロペラがデフォルトになっているが、風車設計と解析も可能である。風車で覚えておく必要のあることは負の(CL、推力、トルク、パワー、レシプロ効率)になることである。複雑になる部分は CLの値が負になるプロペラに対して、風車のブレードの翼型は常にプロペラと上下反対になる。揚力曲線などは原点に対して対称である。

a_0 → -a_0 、 CL_0 → -CL_0
CL_min → -CL_min 、 CL_max → -CL_max 

AEROメニューのREFLコマンドはこの反対にするのを自動でやってくれる。別の方法としてXROTOR.DEFファイルを変えても出来る。風車を設計するときは負のCLによって推力かパワーが負になる。

出力されたPostscriptファイル

Linuxの場合imagemagickなどのソフトを使うと

convert plot.ps <変換後のファイル(例えばplot.jpg)>

を端末上でコマンドを打つと変換できる。画像が粗い場合はconvert -density 250 plot.ps <変換後のファイル.jpg>にすると良い。
横長のpsファイルはうまく変換出来無いことが多いのでその場合はInkscapeIllustratorなどで変化すると良い。

トップレベルメニュー

   QUIT   プログラムを終了
  .AERO   翼型特性サブメニュー
   NAME s 名前をつける

   ATMO r 標準大気から空気のプロパティを設定(rには高度(km)を入力)
   VSOU r 音速の変更(m/s)
   DENS r 空気密度の変更(kg/m^3)
   VISC r 空気の粘度の変更(kg/m/s)

   NACE f ナセルの幾何ファイルの設定
   DUCT r ダクトか自由端かのスイッチ
   VRAT r ダクトの速度比の変更

   ARBI   任意のローターの幾何条件の入力
  .DESI   デザインサブメニュー
  .MODI   修正サブメニュー
  .OPER   設計点以外での特性計算サブメニュー

  .BEND   構造荷重とたわみの計算サブメニュー
  .NOIS   音響解析サブメニュー
   JMAP   Cp vs J operating map

   INTE   指定の半径に補間
   SAVE f ローターをリスタートファイルに保存
   SAVO f ローターを古いタイプのリスタートファイルに保存
   LOAD f ローターをリスタートファイルから読み込み
   WDEF f 現在の設定をxrotor.defファイルに書き込み

   VPUT f 後流の速度データの保存
   VGET f 後流の速度データの読み込み
   VCLR   後流の速度データの消去

   HARD   現在のプロットをPostscript形式で保存
   WIND   風車/プロペラのプロットモード変更
   PLOP   プロットオプション
   DISP   現在の設計点の表示

ARBIコマンドで任意の幾何形状のローターを入力できる。
r/R, chord/R, blade angle.を書いたファイルを読み込む。点と点はスプラインによって結ばれる

.DESIサブメニュー

   INPU    設計パラメータの入力... design rotor
  .EDIT    設計パラメータの編集... design rotor

  .FORM    後流と速度の公式の選択
   ATMO r  標準大気から空気のプロパティを設定(高度km)

   DISP    現在の設計点の表示
   TERS    簡潔出力(terse)か冗長出力(verbose)かの選択
   ITER i  ニュートン法の最大反復回数の変更
   N    i  点の個数の変更
   NAME s  名前の変更

   PLOT i  様々なローターのパラメータのプロット(1〜12))
   ANNO    現在のプロットの注釈
   HARD    現在のプロットをPostscript形式で保存
   SIZE r  プロットの大きさの変更

INPUコマンドで設計を行う。以下の2つはどちらかしか選択できない。
{ advanced ratio , rpm} , {thrust , power}

.MODIサブメニュー

   BLAD i  ブレード数の変更
   MODC    コード長の分布の修正
   MODB    ブレードのねじり角の分布の修正
   SCAL rr 現在のコード長から大きさ変更(AとB) (Chord)new = (Chord)old * [A + B(r/R)]
   TLIN r  線形にねじり角を増やす(deg)
   RTIP r  先端半径の変更(コード長/Rは一定)
   RHUB r  ハブ半径の変更
   RWAK r  ハブの後流のbody半径からの変位の変更
   RAKE r  ブレードの傾斜角の変更
   XPAX r  ピッチ-軸x/cの変更

   CLPI r  Clip current prop at inner radius (preserve chord)
   CLPO r  Clip current prop at outer radius (preserve chord)

   OPTI    現在の平面形でのブレードのねじり角の最適化

   PLOT i  様々なパラメータでのプロット(1〜12))
   ANNO    現在のプロットの注釈
   HARD    現在のプロットをPostscript形式で保存
   SIZE r  プロットの大きさの変更

MODCとMODBコマンドはコードとβの分布をマウスからいじれるとあるが、うまくいかない。
OPTIコマンドは設計と変更の両方を行えるコマンド。ねじり角(βの分布)をMIL循環でコード分布は一定のまま再配置してくれる。ねじり角の"optimizes"の意味である。
このOPTIコマンドは注意しないといけない。MIL循環は失速について考慮していない。また最適な(optimum)MIL分布は全体においての最適とは限らなく、他の動作点ではわるくなるかもしれない。
他の動作点での性能の良し悪しは緩和係数(relaxation factor)による。緩和係数は1.000より小さいほうがオフデザイン時の性能がよくなる傾向がある。

.OPERサブメニュー

   ADVA r   進行率の指定[-]
   RPM  r   rpmの指定
   THRU r   推力の指定[N]
   TORQ r   トルクの指定[N-m]
   POWE r   パワーの指定[W]

   ASEQ rrr 進行率を変化させたときの計算
   RSEQ rrr rpmを変化させたときの計算
   BSEQ rrr ピッチ角を変化させったときの計算
   VSEQ rrr ピッチ角を固定してスピードを変化させたときの計算
   CLRC     貯めたデータを消去
   ADDC     貯めたデータを追加
   CPUT f   貯めたデータをファイルにエクスポート
   CGET f   データをファイルから呼び出し
   CASE i   場合を選ぶ

   ATMO r   標準大気から空気のプロパティを設定(高度km)
   VELO r   巡航速度の変更
   ANGL r   ピッチ角の変更
   PVAR f   エンジンのrpm/power曲線の入力

   FORM     後流と速度の公式の選択

   NAME s   事例の名前の変更
   WRIT f   ファイルに現在の動作点を書き出し
   DISP     現在の動作状態の表示
   INIT     次の解析事例の初期化
   REIN     Re-initialize prop to known operating state
   TERS     出力を簡易か冗長かスイッチ
   ITER i   ニュートン法繰り返し数の最大値の変更
   N    i   点の個数の変更

   PLOT i  様々なローターのパラメータのプロット(1〜12))
   ANNO    現在のプロットの注釈
   HARD    現在のプロットをPostscript形式で保存
   SIZE r  プロットの大きさの変更

ADDC,ASEQ,RSEQ,VSEQ,BSEQコマンドで解析したものは"case accumulator"に貯められる。

.AEROサブメニュー

   DISP  翼型特性の表示
   NEW   新しいセクションを作る
   DEL   指定のセクションの消去
   EDIT  セクションの空力データの編集

   DISP  翼型特性の表示
   LIFT  揚力特性の変更
   DRAG  抗力特性の変更
   MOVE  セクションの位置r/Rの変更
   REFL  風車(windmill)モードのためにCLとCDカーブの正負反転
   PLOT  翼型特性のプロット

   READ  ファイルから翼型特性の読み込み
   WRIT  ファイルへ翼型特性のエクスポート
   PLOT  翼型特性のプロット
   ANNO  現在のプロットに注釈をつける
   HARD  現在のプロットをPostscript形式で保存

.BENDサブメニュー

荷重とたわみの解析

   READ f ファイルからブレードの構造特性読み込み
   EVAL   構造荷重とたわみの評価
   CLR    全てのたわみの削除

   DEFL   Set new twist  =  static  +  structural twist
   REST   Set new twist  =  static twist
   SETS   Set static twist = current - structural twist

   WRIT f 構造解析をファイルに保存

   PLOT i  様々なローターのパラメータのプロット(1〜12))
   ANNO    現在のプロットの注釈
   HARD    現在のプロットをPostscript形式で保存
   SIZE r  プロットの大きさの変更

.NOISサブメニュー

騒音解析

   P   rrr 観測点 x,y,z での音響p(t)の計算
   FOOT rr dB騒音マップの計算
   NTIM i  サンプリング数の変更
   UNIT    単位系の変更 m,ft

   AOC  r  ブレードの断面積area/c^2 の設定
   AFIL f  ファイルからブレードの断面積area/c^2

   PLOT i  様々なローターのパラメータのプロット(1〜12))
   ANNO    現在のプロットの注釈
   HARD    現在のプロットをPostscript形式で保存
   SIZE r  プロットの大きさの変更