「目的」
Opnefoamの回転チュートリアルを用いてヘリコプター・ブレードの基礎的流体解析を行う。
(ブレード周りの流れの様子を視覚で確認することを目的とします。)
「方針」
遠心ポンプの流体解析では、SRFPimpleFoam(非圧縮・非定常乱流解析ソルバー)を用いていましたが、SRFでは外周に行くに従いコリオリ力の誤差が積もり重なって現実では起こらないであろう動きを示していました。この難点によりヘリコプター・ブレード解析では計算(下図)が破綻してしまいました(底面から気流が吸い上げられています)。この難点は回避出来ないと判断して、別のソルバーであるMRF(Multiple Reference Frame)を用いることとしました。
１．解析モデル
１．１ モデル方針
ラジコン・ヘリの500クラスをイメージして以下の形状を検討対象とします。
・ブレード直径 1050[ｍｍ]
・ブレード長 450[mm]
・ブレード厚み 6[mm](メッシュ粗さから適当に厚く設定)
・回転座標系直径 1160[mm](メッシュ粗さからぼこぼこの円になっています)
・解析空間 1600×1600×180[mm]
・ブレード回転数 2800[rpm]程度
・機体重量 2[kg]程度(仮定)
１．２ OpenFoam設定
openfoamの設定を以下に示します。
・メッシュ
1800mm×1800×180[mm]を26x26x12の合計8112マスに分割。
ブレード周りを4倍、8倍、16倍と段々に細分化。
回転系定義面はtopoSetDictにより設定。
・積分幅：adjustTimeStep=yesとしてクーラン数による自動設定。
発散防止のため1巡目のみ1e-8刻み固定で5e-5秒まで計算。
・クーラン数 Co=0.5(巻末「補足1」にて記述)
・回転速度2800[rpm](反時計回り)。
・速度境界：ブレード表面はnonslip、その他はinletOutlet uniform (0 0 0)。
・初期速度：下降気流1[m/s]
・動粘度：空気 1e-5[m2/s]
・解析ソルバーはMRF+PimpleFoam(非圧縮・非定常乱流解析ソルバー)：オリジナルSimpleFoamから変更
・乱流モデルはRAS(kEpsilon)。
・k及びnutファイルの設定は適当(初期速度1[m/s]で参考文献[1]4.5.5項に従う。)
上記以外のパラメータは以下のOpenFoamチュートリアルの物をそのまま使用。
tutorials\incompressible\simpleFoam\mixerVessel2D
本ソルバーはマルチプル・レファレンス・フレーム(MRF)機能を持ったpimpleFoamです。遠心ポンプの流体解析で示した回転系解析空間を静止系解析空間で囲むモデルとなっており、その特徴は以下となります。
・コリオリ力と遠心力を考慮(SRFと同じ)
・回転系と静止系の間で速度及び流量が回転を考慮して引き渡される。
・回転系にある個体と静止系にある個体はその位置を変えることなく解析が進む。
したがって実際の回転系で起きる局所的時間変化は計算されず定常状態へ向かって計算を行うそうです。
２．計算結果
２．１ 下降速度1m/sの場合
以下に回転座標系での速度Urelのz成分とブレード表面圧力pの計算結果を示します。
[y+：Min27.9 Max2310 Ave385]
綺麗に螺旋状のブレード回り流れを視覚で確認することが出来ました。
OpenFoamを用いた遠心ポンプの流体解析「補足４」と同様にopenfoamの汎用関数「forceCoeffs」を用いて推力を計算しました。forceCoeffsの計算結果を以下に図で示します。
上図はまだ完全に収束していませんがCL＝66.9、CmYaw＝-17.1を用います。
forceCoeffsの計算結果より推力Fと反トルクTが以下の式で求められます。
F=1/2*CL=33.5[N](3.4[kgf])
T=1/2*CmYaw=-8.53[Nm]
推力に関しては仮定自重2㎏の約1.7倍となり、500クラスヘリとして妥当な値と考えられます(ピッチ角5度のため更に推力を上げられる状態です)。この値に関しては次項にて検証していきます。
尚forceCoeffsに頼らず、以下の式を計算することでも直接推力が計算出来ます(それなりの手間がかかりますが)。
F=ΣAi×nzi×Pi+Ai×τzi
Ai:i要素面積、nzi:i要素法線ｚ成分、Ｐi:i要素での圧力、τzi=壁面せん断応力ｚ成分
続いて下降流速度を0.1、3、5[m/s]の場合を計算し、結果を前述1[m/s]と一緒にエクセルに示します。
図がわかる通り下降流が大きくなるほど推力Fが減っていき、最終的には推力がマイナスになると予想されます。これは下降流が大きくなっていくとブレードの迎え角が小さくなっていくことが原因と考えられます。
尚ヘリコプターではホバリングから前進運動に移行する過程で揚力が増加する現象が起きます（転移揚力(トランスレーショナル・リフト)(参考文献[2]2章12項))。この現象は高迎角におけるブレード上面剥離が下降流によって緩和されることが原因の一つと考えます。しかし巻末「補足１」にて示しますが、本解析では翼面上の発生する乱流を扱うことが困難なため、この転移揚力を本モデルで解析することは残念ながら困難と考えます。
３．検証
３．１ 速度分布による推力推定
１項にてopenfoamを用いた計算結果を机上計算で検証していきたいと考えます。まずは教科書に出てくる推力推論式で確認していきます。
下面でのZ方向速度分布(ブレード面より85[mm]下方)を下図に示します。下記推力の基本理論式(参考文献[2]12章2項)では円状の一定下降流を仮定していますが、MRFの特質上、上下流のむらがあります。円としての適用は無理ですが、微小面積毎に下記式に当てはめて積分した結果はF＝43.9[N]※1となりました。これはopenfoamの計算結果33.5[N]にそれなりの一致をしています。
F=２×ρ×A×ｖ2
ρ：密度（CFD内部計算では1）
Ａ：ローター面積
v：誘導速度(ブレードによる上下方向の増速分)
※１：計算結果UrelZ(Urel:回転系に回転速度を足す)から下降速度-1[m/s]を引き、かつ上昇気流部位は省いた条件で
メッシュ毎に２×ρ×A×(UrelZ+1)2を計算し、合計した(四角い面全面)。尚ρ=1とした。
３．２ 2D計算値を用いた推定
机上計算のみでMRF計算値を推算出来ると検算として大変便利です。
そのためにまずブレードの2D断面形状を用いたCFDを実施、揚抗比CL/CDを求めました。
〇モデル
・解析空間 2100×180×20[mm]
・ブレード断面形状 50×6[mm]、端部R3[mm]
〇OpenFoam設定
・メッシュ
2100×180×20[mm]を67x24x1の合計1608マスの2Dに分割。
ブレード周りを2倍、4倍、8倍と段々に細分化。
・積分幅：adjustTimeStep=yesとしてクーラン数による自動設定。
・クーラン数 Co=0.5
・速度境界：入口146.6[m/s:]、出口はinletOutlet。ブレード表面はnoSlip、その他はslip。
・動粘度：空気 1e-5[m2/s]
・解析ソルバーはPimpleFoam(非圧縮・非定常乱流解析ソルバー)
・乱流モデルはRAS(kEpsilon)。
・k及びnutファイルの設定は適当(初期速度146.6[m/s]×1%で参考文献[1]4.5.5項に従う。)
・fvSchemesは「OpenFoamを用いた円柱周りの流体解析detail1」と、fvSolutionはMRFと同一とした。
結果が以下となります。
CL=0.528 / CD=0.124 / L/D=4.26
[y+：Min142 Max2960 Ave709]
尚高レイノルズ数2D計算の注意事項を巻末「補足2」に示します。
ブレードを2D形状の連続と仮定し、推力Fを以下の式で表現してみました。
F=∫1/2ρ(rω)2・CL・hdr=1/6*ρCLhω2R3
T=∫1/2ρ(rω)2・CL・r・hdr=1/8*ρCLhω2R4
ρ：密度(ここでは1)
ｒ：ブレード半径 75～525[mm]
ω：回転数2800[rpm]
CL：揚力係数 上述計算より0.528
CD：抵抗係数 上述計算より0.124
h：ブレード幅 50[mm]
∴ F=1/6*ρCLhω2(0.5253-0.0753)*2(左右分)=109[N]
T=1/8*ρCDhω2(0.5254-0.0754)*2(左右分)=10.1[Nm]
推力Fはopenfoamの計算結果33.5[N]の約3倍となっています。
一方反トルクTはopenfoamの計算結果8.53[N]と大体同じ値となっています。
これ以降の検討の方針はなんとかして2D計算から得られるCL/CDを用いて3D計算を表現することを目標とします。
まずは翼の特性をよく表現している揚抗比から検討していきます。
MRF計算結果：CL/CD=33.5÷(8.53/.394)=1.55
代表長さとして6/8*R=6/8*.525=0.394[m]を使用
2D計算結果 ：CL/CD=.528/.124=4.26
よって2D計算値4.26に対し3D計算値1.55は64％減という結果で、ヘリブレードでは単体2D翼に対し揚抗比がかなり悪化しています。
流体解析では「鶏と卵」の関係で因果関係を明確に分けることが難しいですが、あえて主要因を探っていこうと考えます。
悪化原因として以下の４つを想定し、次項から検討してみます。
①2Dと3Dによる影響
②ブレード半径方向での速度分布の影響
③他方ブレードの後流(ウエイク)の影響
④遠心力＆コリオリ力の影響
３．３ 2Dと3Dによる影響
MRF計算による揚抗比が2D計算に対して悪化している原因を推定していきます。
まずは今まで検討してきたブレードを３次元一様流の中にいれた状態をpimpeFoamにて計算してみました。
(速度は2800[rpm]@アーム長500[mm]に対応する146.6[m/s]。OpenFoaam設定はMRFと同じとし、MRF機能のみ停止しました)。
[y+：Min138 Max1920 Ave638]
両翼端で渦が発生している様子が確認出来ます。大きな渦は剥離とセットになっており、剥離は流管の圧力低下を疎外するため揚力低下につながります。計算結果は以下となります。
CL=0.331 / CD=0.121 / L/D=2.73
結果、2D計算値4.26から36%減という状況で、前項MRF計算結果の64%減の大半をこの3D化により説明出来そうです。
念のため速度を半分にした場合も計算しておきました。結果が以下でほぼ一定でした。
CL=0.330 / CD=0.122 / L/D=2.71
この様子を教科書の式を用いてもう少し確認しておきます。
揚力線理論より迎角減少及び抵抗増加に関して(参考文書[3]14.3.2項)
αi = CL/(πAR)＝0.0187[rad]=1.1[deg]
CDi = CL2/(πeAR)=0.0123
αi：誘導迎え角
CDi：誘導抵抗抗力係数
CL：揚力係数 2D計算値0.528(本来は3D計測値となりますがここは2Dを用います)
e：オズワルドの効率係数。揚力が楕円分布で１。今回は適当に0.8を使用
AR：アスペクト比 9
以上より迎角減少分1.1degより揚力は1.1/5=22％減少、抵抗は0.0123増加していると計算されます。この数値から推定される3D推定値は以下となりますが揚力係数CLについてはOpenFoam計算結果とはかなりずれています。
推定CL 0.528(2DCL)×0.78(迎角22%減)＝0.412 vs 0.331(3Dブレード単体)
推定CD 0.124(2DCD)×0.78(迎角22%減)+0.0123(誘導抵抗)=0.109 vs 0.121(3Dブレード単体)
結果CDに関しては大体一致していますが、CLに関しては理論推定よりも20%低い結果となりました。しかしあくまで揚力線理論に対する比較であり3D計算値が安全側の結果のため、このまま使用していきます(理論推定値よりもCLが下がるのは当然の結果とも感じます)。
３．４ ブレード半径方向での速度分布の影響
下図にMRF計算結果から得られるブレード後端の流束の様子を示します。2D計算では渦は垂直方向の平面で表現されますが、下図の流線は水平方向にも移動しています。このため水平方向にも剥離＆渦が発生していると予想されます。抵抗はこの剥離により流管の圧力回復が疎外されることが直接の起因となります。MRF計算では水平方向にも渦が発生し、これが背圧の回復を妨げることで揚抗比が悪化した一因と推測します。
そこでこの点を確認するため以下のような傾斜のある速度分布を持った３DモデルをpimpleFoamで計算してみました(遠心力及びコリオリ力は含まれていません)。
結果が以下となります。
CL=0.186 / CD=0.0701 / L/D=2.66
(速度は半径500ｍｍでの速度146.6で代表計算。よってL/Dのみ他と比較可能な数値となります) 結果、2D計算値4.26から38%減という状況でした。前項の誘導抗力の影響が36％でしたので傾斜速度による影響はわずか2%といったところになりました（見た感じは影響が大きそうに見えたのですが）。
[y+：Min51.4 Max1900 Ave467]
３．５ 他方ブレードの後流(ウエイク)の影響
下図にMRF計算結果から得られるブレード後端の流束の様子を示します。前方ブレードの後流が後方ブレードに到達している様子が見て取れます。後流の中は乱流動粘性係数nutが大きく成長しており、剥離が起き易くなります。すると背面での圧力回復が不十分になるため抵抗の原因となります
[y+：Min133 Max1580 Ave743]
そこでこの点を確認するため以下のような2D空間でブレード2枚を配置し、pimpleFoamで計算してみました。OpenFoam計算設定は3.2項と同じです。
後方ブレードの位置を以下のように修正しました。
・MRF計算においてはブレード先端での双方の円周距離はπR(525[mm])＝1.65[m]。
・下降流速度1[m/s]のため半回転のうちに降下する距離10.7[mm]となす角は0.37°。
・解析空間の観点からブレード間距離を500mmとしたため、500tan0.37=3[mm]高い位置に配置。
結果が以下となります。
前方ブレード CL=0.448 / CD=0.103 / L/D=4.35
後方ブレード CL=0.466 / CD=0.114 / L/D=4.09
結果、前後比較で6%減という状況でした。こちらもある程度の原因を担っているようですが、本筋といった感じはありません。
尚CL計算結果が3.2項2D計算結果0.528と大分違いますが、この点は巻末「補足２」にて確認しています。
３．６ 遠心力＆コリオリ力の影響
最後に遠心力とコリオリ力の影響を検討していくのですが、これらを単独で抽出して検討することは難しいと考えます。
一方前項までの影響検討より3D化で36%、傾斜速度で2%、ウエイクで6%、合計44%となり、目標の64%の7割弱に達しています。
そこでここでは回転数2800[rpm]を変化させることでこの64%という値の変動を見ていき、遠心力とコリオリ力の影響を間接的に推定しようと思います。
コリオリ力と遠心力は以下の式となります(質量m=一定と仮定)。注意すべきはコリオリ力と遠心力の向きが逆さになる点です。符号からするとコリオリ力が遠心力をカバーし、今回の揚抗比を改善する側に働きそうという点です。
・コリオリ力 F=-2ｍvω=-2mrω2
・遠心力 F=mrω2
2800[rpm]、2000[rpm]、1400[rpm]における揚抗比の対2D比率の計算結果を以下に示します。
予想通り回転数が上がるにつれて対2D比率が改善(右肩下り)しています。近似式より回転数0ではその比率が67.5％になっています。よって2800[rpm]でのコリオリ力＆遠心力は3.5％(=67.5-64)となります。但しこの3.5％は揚抗比改善側に働いています。
以上2項をまとめますと以下となります。
揚抗比悪化原因の割合(あくまで目安)
①2Dと3Dによる影響 ：36%
②ブレード半径方向での速度分布の影響 ： 2%
③他方ブレードの後流(ウエイク)の影響 ： 6%
④遠心力＆コリオリ力の影響 : -4%
⑤不明分 ：24%
合計 ：64％(2D/3D比)
完全に明確化出来たわけではありませんが、ある程度の指標になると考えます。
また抵抗の絶対値も確認しておきます。3.4項傾斜速度の計算結果よりDrag=16.95[N]。これに代表長さ0.395mと左右分を掛けると13.4[Nm]。MRF結果の8.6[Nm]に対して1.6倍で揚力ほどの違いはないようです。
よってブレード回転による性能低下は揚力が主で、それにつられて揚抗比も悪化しているという全体像が見えてきました。
3．まとめ
MRF+PimplFoamによるラジコンヘリのブレード回りの流体解析より以下のことが分かりました。
・ブレードにより生成される流束は綺麗な螺旋形状を形成している。
・MRF特有の解法により下降流にはむらが見られる(実際はそのむらが回転し、平均化される)。
・下降流の初期速度が大きくなるほど、ブレードへの迎角が減り、結果として推力が減っていく。
・2D翼性能からヘリブレード性能の概略推算の目途を得ることが出来た。
一方でMRFの特性上以下の点は不明のままとなりました。
・MRFでは一方のブレードに発生する後流を時間的に再現しておらず、他方のブレードとの干渉は定常解として表現される(干渉は見れるが時間的にダイナミックな解ではない)。
「補足１」クーラン数と計算結果
「OpenFoamを用いた円柱周りの流体解析detail1」にてクーラン数に注目しました。ここでクーラン数0.1とクーラン数0.01の場合の円柱CD計算値を試験データと比較しています。 結果、全く同じモデルでありながらCDに約3倍の違いがでてしまいました。
[Drag coefficient vs Reynolds number for all tests from Wieselsberger (1921)]
この時の計算ではクーラン数を下げていくと局所的にepsilonが急増し、それに伴い下記式に従いnutが減少しました。nutが減少すると流体は比較的自由に運動することが可能となり、流束が円柱から剥離するのを抑え込む方に作用したと推察しています。 しかしこのようにクーラン数で結果が大きく変わるレイノルズ数領域(臨界領域)を正確に解析することは実験値との比較なしでは難しいと考えます。
今回の計算のレイノルズは以下で臨界レイノルズ数(2e+5)の領域に入っています。
Re＝V*L/nu=7.7e+5
Ｖ=ω*Ｒ=154[m/s]
ω：2800[rpm]
R：ブレード半径525[mm]
L：代表長さ翼舷長50[mm]
nu：空気の動粘度1e-5[m2/s]
このためクーラン数によって計算値が大きく異なる心配があります。そこで本文で計算したクーラン数Co=0.5に対し、同様の計算をクーラン数Co=0.1で計算しなおしてみました(Co＝0.05は計算時間が膨大になるため断念しました)。
結果を以下の図に示します。
結果としてクーラン数別の揚力CLとCmYawモーメントの計算結果は殆ど一致しています。例え臨界レイノルズ数領域でも迎角5度では上面剥離が少ないため、クーラン数による計算値の差が出なかったと考えています。
「補足２」高レイノルズ数2D計算時の注意事項
3.2項の2Dブレード解析にて不思議な現象に出会いましたので、ここにまとめておきます。
「事象」
2D解析を行う上で後流解析空間をどれくらいにするかでパラメータスタディーを実施しました。
結果が以下となります。解析空間全長が下から110、160、210、260[cm]の４つとなります。
CL/CDの様子をグラフにすると以下になります。解析空間を増やしていっても収束するどころか260cmでは急減しています。これではなにを使っていいか分かりません。
ちなみにブレード長は5cm(=c)、解析空間全長110[cm]の場合でも前方6ｃ、後方16cで一般的には十分な空間と考えます。
結論から申し上げますと前方空間が足りませんでした。このため乱流動粘性係数nutが十分拡散する前にブレードに到達し、CLが不安定になっていたようです。前方空間を50cm伸ばした計算結果を以下に示します（下段：未伸長、上段：50[cm]伸長）。通常の計算ではこれほどまで前方空間を取る必要はないのですが、この計算ではレイノルズ数が7e+5(=146.6*0.05/1e-5)で、臨界レイノルズ数2e+5の領域に達しているため、このようなケアが必要なのだと考えます。
[y+：Min120 Max1570 Ave677@110+50cm]
結果、安定したCL、CDが得られるようになりました。板の理論値CL=2πsin5[deg]=0.548とも十分一致しています。但し空間を必要以上に広げていくと不安定になっています。原因はつかめ切れていませんが、CFDの注意すべき点だと感じます(試験値とのコリレーションが望まれます)。
参考文献
[1]OpenFoamによる熱移動と流れの数値解析 第2版 OpenFoamCAE学会編 森北出版
[2]ヘリコプターの工学と操縦 筒井 善直著 酣燈社
[3]空気力学入門 日本航空宇宙学会編 丸善出版
HA-technodefy 