「目的」
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・CD・r・hdr=1/8*ρCDhω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を用いて3D計算を表現することを目標とします(CLは試験値とCFD計算及び理論値が比較的良い一致を示しますが、CDは物体表面の質、CFDメッシュ粗さに強く依存するため今回は検討対象外としました)。
MRF計算結果 ：F(揚力)=33.5[N]
2Dを用いた机上単純計算：F(揚力)=109[N]
よって2D机上計算値109[N]に対しMRF計算値33.5[N]は31％(=33.5/109)と非常に小さい値になってしまいました(本来ならば一致して欲しい値です)。
流体解析では「鶏と卵」の関係で因果関係を明確に分けることが難しいですが、あえて主要因を探っていこうと考えます。
悪化原因として以下の3つを想定し、次項から検討してみます。
①誘導迎角の考慮
②他方ブレードの後流(ウエイク)の影響
③遠心力＆コリオリ力の影響
３．３ 誘導迎角の考慮
MRFによるCL計算値が2Dを用いた机上計算に対して悪化している原因を推定していきます。
ここでは以下の2点を検討します。
・3D化に伴う誘導迎角の考慮
(揚力線理論(参考文献[3]14.3項)：Openfoamを用いた2D翼と3D翼の空力係数検討参照の事)
・半径方向の速度変化の考慮
この２点を確認するため以下のような傾斜のある速度分布を持った３DモデルをpimpleFoamで計算してみました(遠心力及びコリオリ力は含まれていません)。
そこでこの点を確認するため以下のような傾斜のある速度分布を持った３DモデルをpimpleFoamで計算してみました(遠心力及びコリオリ力は含まれていません)。
結果が以下となります。
[y+：Min51.4 Max1900 Ave467]
F(揚力)=45.1[N]
結果、2Dを用いた単純机上計算値109[N]に対し41%という状況で、MRF計算値33.5[N]に対してかなり近い値となりました（主要因と見なせます）。
この様子を教科書の式を用いてもう少し確認しておきます。
Openfoamを用いた2D翼と3D翼の空力係数検討にて示したのですが、NACA0012という翼形状では風洞試験から2Dから3Dに移行すると誘導速度により揚力係数CLは60%(40%減)になってしまいます。
減少の主要因は吹きおろし(誘導速度)による実効迎角の減少でした。ここでも誘導速度の机上計算より実効迎角の減少を計算してみます。
ビオ・サバールの法則より翼上循環の微分より局所的誘導速度が求められます(参考文献[3]14.3.1項)。
dw=(1/4π)*[1/(y-y0)]*(dΓ(y)/dy)*dy
固定翼機の矩形翼ではΓ(y)=一定で
dΓ(y)/dy=0
となりますが、回転翼では
dΓ(y)/dy=4πωRsinα
Γ=4πURsinα(薄板公式を代用)
U=yω
R=C/4
C：ブレード幅50[mm]
ω：2800[rpm]=293.2[rad/s]
α：5[deg]
このdΓ(y)/dyを用いてdwをy=Ri～Roまで積分することで座標y0での吹きおろし速度w(y0)を計算します。
w(y0)=Rω/4*sinα*[LN(Ro-y0)-LN(y0-Ri)]
Ro：525[mm]
Ri：75[mm]
エクセルにて計算した誘導速度(吹きおろし)と誘導迎角を以下に示します。
(翼両端は翼対数関数上、発散してしまうため、誘導速度を両隣の値で代用しました。)
結果、興味深い様子がうかがえます。ブレード先端では迎角0.5度相当の吹きおろしが発生しているのに対し、ブレード付け根では迎角3度相当の吹き上げが発生しています。この点をOpenFoam計算結果で確認してみました。
確かにブレード付け根(図左側)では吹きおろし速度(Uz)が急減していることを確認出来ました(もともと5度の迎角があるため完全に吹きおろしがなくなるわけではありません)。
このエクセル誘導迎角と2D計算値CL=0.528を用いてブレード全体の揚力を計算してみました。
結果が以下です。
F(揚力)=54.6[N]
これはopenfoam計算結果45.1[N]と十分近い値となっています。つまり誘導迎角と速度分布による影響を揚力線理論により机上計算した結果とopenfoam計算結果がよく一致したことになります。
３．４ 他方ブレードの後流(ウエイク)の影響
下図に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
結果、CLに関しては後方ブレードの方が良くなってしまいました。これは巻末補足1及び巻末補足2で触れていますが、揚力係数CLはnutに大きく影響を受けるためと考えます。
以上、後流(ウエイク)の影響は本計算では影響を確認出来ませんでした(その程度と考えることも出来ます)。
尚CL計算結果が3.2項2D計算結果0.528と大分違いますが、この点は巻末「補足２」にて確認しています。
３．５ 遠心力＆コリオリ力の影響
最後に遠心力とコリオリ力の影響を検討していくのですが、これらを単独で抽出してCFD計算することはまさしくMRF解析であり別の方法で計算することは難しいと考えます。
そこでここでは回転数2800[rpm]を変化させることでこの31%という値の変動を見ていき、遠心力とコリオリ力の影響を間接的に推定しようと思います。
コリオリ力と遠心力は以下の式となります(質量m=一定と仮定)。注意すべきはコリオリ力と遠心力の向きが逆さになる点です。符号からするとコリオリ力が遠心力をカバーし、今回のCLを改善する側に働きそうという点です。
・コリオリ力 F=-2ｍvω=-2mrω2
・遠心力 F=mrω2
2800[rpm]、2000[rpm]、1400[rpm]における揚抗比の対2D比率の計算結果を以下に示します。
予想通り回転数が上がるにつれて対2D比率が改善(右肩上がり)しています。近似式より回転数0ではその比率が26％になっています。よって2800[rpm]でのコリオリ力＆遠心力の影響は5％(=31-26)となります。但しこの5％はCL比改善側に働いています。
尚本推定は回転数2乗の効果を遠心力とコリオリ力のみとしたかなり強引な仮定であり、あくまで参考となります。
以上2項をまとめますと以下となります。
揚抗比悪化原因の割合(あくまで目安)
①誘導迎角及び速度分布による影響 ：59%減(41%)
②他方ブレードの後流(ウエイク)の影響 ：0%(不明)
③遠心力＆コリオリ力の影響 :-5%(増分)
④不明分 ：15%減
合計 69％減(31%)
完全に明確化出来たわけではありませんが、ある程度の指標になると考えます。
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 