「目的」
2次元の翼の揚力係数よりも3次元の揚力係数は相当下がる。そこで抵抗係数がどうなるかを確認する。
「方針」
OpenFoamにて2次元翼と3次元翼の空力係数を計算し、そのメカニズムを計算結果から推定してみようと思います。
尚ソルバーとしてOpenfoamを用いた2D翼の抵抗係数確認での検討結果を受けてSimpleFoam-RAS(k-ωSST)を使用した。
１．理論
１．１ 試験値
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抵抗係数
出典：NASA-TM-4074
「Effects of Independent Variation of Mach and Reynolds Numbers on the Low-Speed Aerodynamic Characteristics of the NACA 0012 Airfoil Section」
入手先：NTRS-NASA Technical Report Server
・レポート-3：3D抵抗係数：
出典：NACA-TR-647
「TESTS OF NACA 0009, 0012, AND 0018 AIRFOILS IN THE FULL-SCALE TUNNEL」
入手先：NTRS-NASA Technical Report Server
（１）翼形状
レポート-1にて定義される翼断面形状のうち、今回も揚力係数の時と同じNACA-0012を選定しました。
式は下式で定義され、形状は下図(半断面)となります。
±y=0.6×(0.29690√x-0.12600x-0.35160x^2+0.28430x^3-0.10150x^4)
（２）３D揚力係数
レポート-3においてNACA-0012矩形翼(下図)の抵抗係数の試験値が示されています。
その中から迎角10degと5degの抵抗係数を代表として抽出します。
・形状 NACA-0012 コード長6ft(1.83m)×スパン36ft(10.97m)
・抵抗係数 0.04＠迎角10deg、0.0138＠迎角5deg
（試験条件：(推定1気圧)、速度87fps(26.5m/s)、レイノルズ数3.3e+6）
（３）2D揚力係数
レポート-2においてNACA-0012矩形翼の空力係数の風洞試験値が示されています。
複数の試験結果からレポート-3のレイノルズ数に近いものを抽出して下図に示します。
（試験条件：マッハ数0.15、レイノルズ数3.94e+6）
１．２ 理論値
風洞試験結果(迎角10deg)においてCD2D=0.00934に対し、CD3D=0.04と約4倍と大きく変わっています(これでは2Dと3Dがまるで別物です)。これに対応する計算として、まずは楕円翼の理論計算値を示します(参考文献[1]14.3.2章)。
CDi = CL2/(πAR)=0.0342
CDi：誘導抗力係数
CL’=dCL’/dα×10deg=0.803
AR=6：アスペクト比(レポート-3の模型)
dCL’/dα=a0/(1+a0/πAR)=4.603：誘導抵抗による揚力係数の変化
a0=CL/10deg=6.091
CL=1.063：揚力係数(レポート-2の試験値)
結果、試験値CDi=CD3D-CD2D=0.03とほぼ一致しました。
別物に見えたのは誘導抗力が大きな割合を占めているからのようです(但し本試験は矩形翼)。
続いて同じ議論を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：26.5[m/s] R：6ft/4(ジェーコスフキ変換より:平板近似)
エクエルによりグラフにすると以下となります。
中央で1.66deg、縁では無限発散してしまいます。そこで縁は9.69degで切り捨てしました(平均3.32deg)。この誘導迎角に伴い局所揚力係数CLがグラフのように減少します(試験値CL2D=1.063を基準に計算)。下図に示す通り、減少した揚力L’は基準座標に対して傾きます。そしてその分力が抵抗となります。これをエクセルに適用すると
CDi=ΣCL*(10-αi)/10×sinαi=0.0308
結果、試験値CDi=CD3D-CD2D=0.03とほぼ一致しました。
よって抵抗係数は2D試験値から良く3D試験値を予言出来そうです。
この辺りをOpenfoamを用いて検証していきます。
尚このエクセル計算と楕円翼公式(最小抵抗)の結果がほぼ同じになってしまっているのは、エクセル計算が翼端で適当な切り捨てを行っている点が影響していると考えます。
２．2D翼揚力係数
２．１ 解析モデル＆OpenFoam設定
3D検討の前にまず2D翼についてopenfoamによりCD2D=0.00934を確認します。
・モデル
NACA-0012 2D翼 コード長6[ft](1.83[m])
迎角α：10deg
・空間メッシュ
28(長さ)×14(高さ)×0.05(幅)[m]の直方体を40×30×1の合計1200マスに分割。
翼回りを2倍、4倍、8倍、16倍、64倍と段々に細分化。
・速度分布 左端から一定流26.5[m/s](レポート-3風洞試験に合わせた)
・動粘度：空気 1.5e-5[m2/s]
・解析ソルバー：SimpleFoam(非圧縮・定常乱流解析ソルバー)
・乱流モデル：RAS(k-ωSST)。
・k、nut及びomegaファイルの設定は適当(速度26.5[m/s]で参考文献[2]4.5.5項に従う。)
・fvSchemes及びfvSolutionは以下を使用
\tutorials\incompressible\simpleFoam\motorBike
２．２ 計算結果 計算結果を以下に示します。
Openfoam関数forceCoeffsを用いた空力係数の計算結果を試験値グラフ上で示します。
結果、収束したCD＝0.0388は2D風洞試験値0.00934に対し、約4倍でOpenfoamを用いた2D翼の抵抗係数確認の検討とおよそ同じ結果となりました(計算方法が同じなので当然と言えます)。 2D計算では試験値とかなり離れた結果となりました。
３．3D翼揚力係数
３．１ 解析モデル＆OpenFoam設定
Openfoamにより3D風洞試験値CD=0.04を確認します。
・モデル
NACA-0012 3D翼 コード長6[ft](1.83[m])、スパン長36[ft](10.97[m])
迎角α：10deg
・空間メッシュ
28(長さ)×14(高さ)×20(幅)[m]の直方体を40×30×14の合計16800マスに分割。
翼回りを2倍、4倍、8倍、16倍、64倍と段々に細分化。
・速度分布 左端から一定流26.5[m/s](レポート-3風洞試験に合わせた)
・動粘度：空気 1.5e-5[m2/s]
・解析ソルバー SimpleFoam(非圧縮・定常乱流解析ソルバー)
・乱流モデル：RAS(k-ωSST)。
・k、nut及びomegaファイルの設定は適当(速度26.5[m/s]で参考文献[2]4.5.5項に従う。)
・fvSchemes及びfvSolutionは以下を使用
\tutorials\incompressible\simpleFoam\motorBike
２．２ 計算結果 計算結果を以下に示します。
ここでｙ+最大値648は少し大きいですが、図示してみると翼後縁端に僅かに出ているのみ(下図)で解析上、影響ないと考えました。
Openfoam関数forceCoeffsを用いた空力係数の計算結果を以下にグラフで示します。
不思議なことにほぼ同じレベルのメッシュ粗さと同じソルバー(SimpleFoam+k-ωSST)でありながら、3D計算ではCL、CD共に風洞試験値とほぼ一致しています(2D計算では抵抗係数が試験値と4倍の開きがありました)。
この点を確認するため、試験結果を更に詳しく見ていきます。
まずは翼端渦の様子を可視化してみます。 翼端渦により吹きおろしが発達している様子が確認できます。尚通常の翼表面流れの下降速度を取り除くため、翼後縁から8[m]後方の断面を取り出しました(完全に取り除ける訳ではありません)。
スパン方向の下降速度の分布を机上計算と比較してみます(下図)。 机上計算とOpenFoamによる計算結果が良く一致しています。これより翼端渦による誘導速度(吹きおろし)が良く再現され、それに伴い揚力係数と抵抗係数が机上計算及び風洞試験値と一致したと考えられます。
更に誘導抵抗そのものを確認してみます。3D翼計算結果を5[cm]幅の短冊状に切り取り、その部位の抵抗係数を求めてみました(下図)。
・抵抗係数(CD_3D_CFD)は中央部で最も小さく、翼端に行くにしたがって大きくなる。但し中央部でも相当な抵抗を発生しており、抵抗は翼全体で生成させている。
・抵抗係数を圧力による係数部(CD_press)と粘性による係数部(CD_viscos)に分けるとCD_viscosはスパン方向であまり変化がない。このことは誘導抵抗が主に圧力に影響を及ぼすことを示しています。
念のため迎角5degにおいて同じ計算を行ってみました。
結果をグラフに書き足します。
迎角5degにおいても風洞試験値と良い一致をみることが出来ました。
以上3DOpenfoamの計算は教科書で示される誘導抵抗を良く表現出来ており、更に風洞試験値とも良く一致することを確認出来ました。
４．まとめ
・2DのOpenfoam計算値は風洞試験よりも4倍に増えてしまい、あまり良い一致とは言えない結果となりました(但しそもそも議論している抵抗値が小さいため、差が大きく見えているという観点があります)。
・一転して3DのOpenfoam計算値と風洞試験と良く一致しています。しかも同時に揚力係数も良く一致しました。Openformの計算ロジックは2Dよりも3Dの方が相性がよさそうです。
参考文献
[1]航空宇宙工学テキストシリーズ 空気力学入門 日本航空宇宙学会編 丸善出版
[2]OpenFoamによる熱移動と流れの数値解析 第2版 OpenFoamCAE学会編 森北出版
HA-technodefy 