「目的」
2次元の翼の揚力係数よりも3次元の揚力係数は相当下がる。このメカニズムを調べる。
「方針」 OpenFoamにて2次元翼と3次元翼の空力係数を計算し、そのメカニズムを計算結果から推定してみようと思います。
１．理論
１．１ 試験値
OpenFoamの計算に先立ち、比較対象とすべき試験値を定めておきます。
今回は入手し易さから以下を選定しました。
・レポート-1：翼形状
出典：NACA-TR-460「THE CHARACTERISTICS OF 78 RELATED AIRFOIL SECTIONS FROM TESTS IN THE VARIABLE-DENSITY WIND TUNNEL」
入手先：NTRS-NASA Technical Report Server
・レポート-2：2D揚力係数
出典：NACA-TM-100019「A Critical Assessment of Wind Tunnel Results for the NACA 0012 Airfoil」
入手先：NTRS-NASA Technical Report Server
・レポート-3：3D揚力係数：
出典：NACA-WR-L-406「SOME EFFECTS OF REYNOLDS AND MACH NUMBERS ON THE LIFT OF AN NACA 0012 RECTANGULAR WING IN THE NACA 19-FOOT PRESSURE TUNNEL」 入手先：NTRS-NASA Technical Report Server
（１）翼形状
レポート-1にて定義される翼断面形状のうち、今回はNACA-0012を選定しました(適当)。
式は下式で定義され、形状は下図(半断面)となります。
±y=0.6×(0.29690√x-0.12600x-0.35160x2+0.28430x3-0.10150x4)
（２）３D揚力係数
レポート-3においてNACA-0012矩形翼(下図)の揚力係数の試験値が示されています。
・形状 NACA-0012 コード長24inch(0.610[m])×スパン144inch(3.66[m])
・揚力係数 0.316＠迎角5deg
（試験条件：１気圧、マッハ数0.073、レイノルズ数1.07e+6）
（３）2D揚力係数
レポート-２にてNACA-0012翼形状の2D揚力係数の試験値が示されています。
・揚力微係数 βdCL/dα=0.105@レイノルズ数1e+6
迎角5deg及びβ=√(1-M2)＝√(1-0.0732)＝0.997より
揚力係数CL=0.527
ちなみに薄板の理論値はCL=2πα=0.548となり十分近い値になっています(参考文献[1]13.5.2章)。
１．２ 理論値
試験値からはCL2D=0.527に対し、CL3D=0.316と60％となっています。
一方当事業所で調べた限りでは2D翼と3D矩形翼の揚力係数の変化に関する理論計算値は見当たりませんでした。その代わり以下に示す楕円翼における理論計算値を示します(参考文献[1]14.3.3章)。
dCL/d(α-αi) = a0
CL：楕円翼揚力係数
α：迎角
αi：=CL/π(AR)、誘導迎え角(翼端渦Γによる吹きおろし)
a0：2D翼での揚力微係数(揚力係数の迎角微分)
AR：アスペクト比
上式を積分＆微分することで下式が得られます。
dCL/dα = a0/(1+a0/πAR）
この式を拡大解釈するとAR⇒∞ではdCL/dα=a0となります。これは2D翼揚力微係数と見て取れます。この式に試験値dCL2D/dα=0.527/5[deg]＝6.039を代入してみると(レポート-3よりAR＝6)
dCL/dα = 6.039/(1+6.039/6π）= 4.574
割合にすると4.574/6.039＝76％となり楕円翼の理論を利用するとこの段階で24％減の説明がつくことになりますが、16%足りません。この点は楕円翼理論式のためと考えます。
そこで同じ議論を3D矩形翼に適応してみます(当事業所オリジナル:他機関検証無)。
3D矩形翼では翼中央で発生する循環Γ0は一定で翼端で翼端渦に移行します。この場合、翼で発生する誘導速度は下式で計算出来ると考えます(※１)。
v1=0.5*Γ0/2π(b/2+y)
v2=0.5*Γ0/2π(b/2-y) （Γ0の前の0.5倍は翼端渦が半無限のため）
∴v=v1+v2=Γ0*b/π(b2-4y2)
誘導迎角αi=v/Uより
αi=v/U=Γ0*b/πU(b2-4y2)
Γ0=4πURsinθ
U：25[m/s]
R：24inch/4(ジェーコスフキ変換より)
グラフにすると以下となります。
中央で0.83[deg]、縁では無限発散してしまいます。そこで縁は6.1[deg]で切り捨てしました。 結果、誘導迎角の平均は1.927[deg]となりました。
ここで楕円翼で用いられる下式をここでも代用します。
αi=CL/πAR
∴πAR＝CL/αi=15.67：πARを変数と仮置き(オズワルドの効率係数を含んだ形)。
CL=0.527：CL2Dを使用
αi=1.927[deg]
本式をdCL/dα = a0/(1+a0/πAR）に代入すると
dCL/dα = 6.039/(1+6.039/15.67）= 4.359
割合にすると4.359/6.039＝72％となり楕円翼の76%よりも4%悪化した結果となりました。
つまり2D翼での揚力係数は3D翼になる段階で翼端渦による吹きおろしのために少なくとも28％は減ってしまうことが見通せました(依然12%の違いが残っていますが)。
この辺りの現象をopenfoamを用いて視覚的に確認していこうと考えます。
※１：循環Γによる誘導速度はv=Γ/2πr(ビオ・サバールの法則)で計算しています。これはOpenFoamを用いた循環Γと渦度の観察で示したストークスの定理とは状況が異なります。ストークスの定理では円周は渦度の縁でしたが、誘導速度では縁という訳ではありません。ビオ・サバールの法則は複素ポテンシャル空間で速度場を表現するために導入された数学上の式なのに対し、ストークスの定理は実空間での渦度の性質を表現しており、似た式ではありますが使い訳が必要と考えます。
２． 2D翼揚力係数
２．１ 解析モデル＆OpenFoam設定
3D検討の前にまず2D翼についてopenfoamにより試験値0.527を確認します。
・モデル
NACA-0012 2D翼 コード長24inch
迎角α：5[deg]
・空間メッシュ
10(長さ)×5(高さ)×0.05(幅)[ｍ]の直方体を40×20×1の合計800マスに分割。
翼回りを2倍、4倍、8倍、16倍、32倍と段々に細分化。
・積分幅：adjustTimeStep=yesとしてクーラン数による自動設定。
・クーラン数 Co=0.3
・速度分布 左端から一定流25[m/s](M=0.073相当:風洞試験に合わせた)
・動粘度：空気 1.5e-5[m2/s]
・解析ソルバー PimpleFoam(非圧縮・非定常乱流解析ソルバー)
・乱流モデルはRAS(kEpsilon)。
・k及びnutファイルの設定は適当(速度25[m/s]で参考文献[2]4.5.5項に従う。)
・fvSchemes及びfvSolutionは以下を使用
\tutorials\incompressible\pimpleFoam\RAS\pitzDaily
２．２ 計算結果
計算結果を以下に示します。
openfoam関数forceCoeffsを用いた空力係数の計算結果を以下にグラフで示します。
結果、収束したCL＝0.443は2D試験値0.527に対し、やや小さい近い値となってしまいました。
メッシュの改善を行うとCL=0.497位まで改善するのですが、次項3D計算と同じメッシュ粗さにするため、このままとします。当社経験上、レイノルズ数が小さい間は試験値とCFD結果は良好な一致を示すのですが、実用領域である2e+5を超えてくると試験値とCFD結果に顕著なずれが出てくると思います。
尚この計算値はメッシュ粗さとクーラン数に強く依存しており、CL=0.3～0.5位まで簡単に変動してしまいます。こういう点からもCFD計算結果の数値の絶対値を議論する際は試験とのキャリブレーションが必須と考えます。
３．3D翼揚力係数
３．１ 解析モデル＆OpenFoam設定
3Dopenfoamにより試験値CL=0.316を確認します。
・モデル
NACA-0012 3D翼 コード長24inch、スパン長144inch
迎角α：5[deg]
・空間メッシュ
10(長さ)×5(高さ)×8(幅)[ｍ]の直方体を40×20×8の合計6400マスに分割。
翼回りを2倍、4倍、8倍、16倍、32倍と段々に細分化。
・積分幅：adjustTimeStep=yesとしてクーラン数による自動設定。
・クーラン数 Co=0.3
・速度分布 左端から一定流25[m/s](M=0.073相当:風洞試験に合わせた)
・動粘度：空気 1.5e-5[m2/s]
・解析ソルバー PimpleFoam(非圧縮・非定常乱流解析ソルバー)
・乱流モデルはRAS(kEpsilon)。
・k及びnutファイルの設定は適当(速度25[m/s]で参考文献[2]4.5.5項に従う。)
・fvSchemes及びfvSolutionは以下を使用
\tutorials\incompressible\pimpleFoam\RAS\pitzDaily
２．２ 計算結果
計算結果を以下に示します。
openfoam関数forceCoeffsを用いた空力係数の計算結果を以下にグラフで示します。
CFD結果はCL＝0.362となり3D試験値0.316に対しそれなりに近い値なのですが、2D_CFD結果(CL=0.443)との割合は0.362/0.443=82%となり、当事業所で実施した机上推定値72%(1.2項)と10%ずれてしまいました。このずれの原因を以下に考えてみます。
原因を考えるにあたり、まず吹きおろしの様子をビジュアル化してみます。上述エクセルでは中央部でも吹きおろしが0.36[m/s]程度発生しているはずですが3D計算では通常の翼表面流れ(25[m/s]で5[deg]の向き)と相まって吹きおろしを直接見て取れません。ただ吹きおろしの様子がきれいな円形ではなく、解析空間の狭さから歪んでいることが見て取れます。
そこで直接吹きおろしを見出す代わりに翼端渦Γと翼中央の循環Γ0の比較を行ってみます。
下図の領域にてQ基準＞0の領域の渦度(vorticity_y)を積分してみました。
(Q及びvorticityはOpenFoamを用いた循環Γと渦度の観察にて説明しております。)
結果はΓ=2.56(y=4[ｍ]で極大)
これはNACA0012より予想される循環値Γ0(下記)に対して76%(2.56/3.376)になっていました。
Γ0=L/(ρU)=3.376[m2/s]
ρ：密度、ここでは「１」
U=25[m/s]
L=0.5ρU2*CL*S=0.5*252*0.443*24″=84.39[N]
すなわち3D計算では翼中央で発生している循環Γ0が翼端渦に移行しそこなっている分、誘導速度(吹きおろし)が減っていると予想します(この計算では24%程度移行出来なかったようです)。そうすると前述のCFD計算値82％と机上計算値72%の差をだいたい説明出来ると考えます(18÷0.76＝24%≒28%)。
なお原因としては以下を考えています。
・翼端メッシュの粗さ：循環渦をきれいに90度向きを変えれていない
・解析空間狭さ：吹きおろしがきれいに円形を描けない
結果、CFD計算は机上計算と一致する方向だと分かりました。そうだとすると風洞試験結果であるCL2D=0.527とCL3D=0.316の割合60％と、上述机上計算値72%との差12％の原因が掴めません。
これ以上は風洞試験結果をより詳細に検討する必要があり、本稿ではここまでが限界と考えました。あくまで推定ですが以下を要因の一つと感じます。
・以下の式で片翼端渦の吹きおろしを計算しています。
v1=0.5*Γ0/2π(b/2+y)
この際、翼端渦は半無限のため0.5倍掛けていますが、吹きおろしの影響は前縁だけではなく、翼全体及びある程度後方まで及ぶと考えられます。例えば翼舷長×2の距離を影響範囲と考えると0.5倍ではなく、0.5×(1+sin(tan-1(2/6))=0.658倍(ビオ・サバール積分範囲の+90～-18degに相当)となり、試験値を十分説明出来る状態になります。
実務上はCFD計算による誘導速度は試験値よりも少な目になってしまうことを留意すべきと考えます。
４．まとめ
・他論文より3D矩形翼の揚力係数は2D翼に対して60％に減じています。
・OpenFoamと机上計算により同程度の減少(72%)を確認出来ましたが、試験値との違い12％分を説明出来切れていません。
・翼端渦により生成される誘導速度の様子をビジュアル化出来ました。しかし主たる流れとの分離は困難です(主流が25[m/s]に対し、誘導速度はせいぜい1[m/s]程度しかないため)。
・OpenFoam計算では翼中央の循環が翼端で翼端渦に漏れなく移行出来るかが重要でした。
参考文献
[1]航空宇宙工学テキストシリーズ 空気力学入門 日本航空宇宙学会編 丸善出版
[2]OpenFoamによる熱移動と流れの数値解析 第2版 OpenFoamCAE学会編 森北出版
HA-technodefy 