「目的」
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)
・空力係数(グラフ読み取り)
CL3D=0.353、CD3D=0.0138 ＠迎角5deg
CL3D=0.719、CD3D=0.0400 ＠迎角10deg
（試験条件：(推定1気圧)、速度87fps(26.5m/s)、レイノルズ数3.3e+6）
（３）2D揚力係数
レポート-2においてNACA-0012矩形翼の空力係数の風洞試験値が示されています。
複数の試験結果からレポート-3のレイノルズ数に近いものを抽出して下図に示します。
（試験条件：マッハ数0.15、レイノルズ数3.94e+6）
上図に3D風洞試験値も入れておきました。
結果CL3D/CL2D=0.719/1.063=68%と32%もの揚力係数の低下が見られます。
そしてCD3D/CD2D=0.04/0.00934=428%と抵抗係数は約4倍の増加となりました。
この辺りを以降で観察していきます。
１．２ 理論値
（１）楕円翼における揚力係数
迎角10degの風洞試験ではCL2D=1.063に対し、CL3D=0.719と2Dに対して3Dは68％となっています。これに対応する計算として、まずは楕円翼の理論計算値を示します(参考文献[1]14.3.3章)。
dCL’/dα=a0/(1+a0/πAR)=4.603：誘導抵抗による揚力係数の変化
a0=CL/10deg=6.091[1/rad]
CL=1.063：揚力係数(レポート-2の試験値)
AR=6：アスペクト比(レポート-3の模型)
CL’=dCL’/dα×10deg=0.803
結果、CL’=0.803で試験値CL3D=0.719とそれなりに近い値となりました。
（２）楕円翼における抵抗係数
迎角10degの風洞試験ではCD2D=0.00934に対し、CD3D=0.04と約4倍と大きく変わっています(これでは2Dと3Dがまるで別物です)。これに対応する計算として、(1)項と同様に楕円翼の理論計算値を示します(参考文献[1]14.3.2章)。
CDi = CL’2/(πAR)=0.0342
CDi：誘導抗力係数
CL’=0.803：(1)項計算値
結果、試験値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’は基準座標に対して傾きます。そしてその分力が揚力と抵抗となります。これをエクセルにて積算計算してみました。
CL’=ΣCL*(10-αi)/10×cosαi=0.7095
CDi=ΣCL*(10-αi)/10×sinαi=0.0308
結果、
・揚力係数
試験値CL3D=0.719とほぼ一致しました。
・抵抗係数
試験値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
翼中央付近(x=1)での想定される境界層厚は以下となります(参考文献[1]13.6項)。
・層流境界層の場合 δ=5.2x/√Rex=5.2*1/√(26.5*1/1.5e-5)=4[mm]
・乱流境界層の場合 δ=0.37x/(Rex)0.2=0.37*1/(26.5*1/1.5e-5)0.2=21[mm]
一方上記モデルの翼周りメッシュ厚さは最小7[mm]で、乱流境界層と見た場合はギリギリ足りているように見えます。但しこの点は「補足１」にて見直してみます。
２．２ 計算結果 計算結果を以下に示します。
Openfoam関数forceCoeffsを用いた空力係数の計算結果を試験値グラフ上で示します。
結果、
・揚力係数
Openfoam計算結果CL=0.853は2D風洞試験値1.063に対し80％の結果となりそれなりの一致となりました。
・抵抗係数
Openfoam計算結果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
翼中央付近(x=1)での想定される境界層厚は2Dと同じ以下となります(参考文献[1]13.6項)。
・層流境界層の場合 δ=5.2x/√Rex=5.2*1/√(26.5*1/1.5e-5)=4[mm]
・乱流境界層の場合 δ=0.37x/(Rex)0.2=0.37*1/(26.5*1/1.5e-5)0.2=21[mm]
一方上記モデルの翼周りメッシュ厚さも2D同様に最小7[mm]で、乱流境界層と見た場合はギリギリ足りているように見えます。
２．２ 計算結果 計算結果を以下に示します。
ここでｙ+最大値648は少し大きいですが、図示してみると翼後縁端に僅かに出ているのみ(下図)で解析上、影響ないと考えました。
Openfoam関数forceCoeffsを用いた空力係数の計算結果を以下にグラフで示します。
不思議なことにほぼ同じレベルのメッシュ粗さと同じソルバー(SimpleFoam+k-ωSST)でありながら、3D計算ではCL、CD共に風洞試験値とほぼ一致しています(2D計算では抵抗係数が試験値と4倍の開きがありました)。
この点を確認するため、試験結果を更に詳しく見ていきます。
まずは翼端渦の様子を可視化してみます。 翼端渦により吹きおろしが発達している様子が確認できます。尚通常の翼表面流れの下降速度を取り除くため、翼後縁から8[m]後方の断面を取り出しました(完全に取り除ける訳ではありません)。
スパン方向の下降速度の分布を机上計算と比較してみます(下図)。 机上計算とOpenFoamによる計算結果が良く一致しています。これより翼端渦による誘導速度(吹きおろし)が良く再現され、それに伴い揚力係数と抵抗係数が机上計算及び風洞試験値と一致したと考えられます。
更に誘導抵抗そのものを確認してみます。3D翼計算結果を5[cm]幅の短冊状に切り取り、その部位の揚力係数と抵抗係数を求めてみました(下図)。
・揚力係数(CL_3D_CFD)は中央部で最も大きく、翼端に行くにしたがって小さくなる。
・抵抗係数(CD_3D_CFD)は中央部で最も小さく、翼端に行くにしたがって大きくなる。但し中央部でも相当な抵抗を発生しており、抵抗は翼全体で生成させている。
・抵抗係数を圧力による係数部(CD_press)と粘性による係数部(CD_viscos)に分けるとCD_viscosはスパン方向であまり変化がない。このことは誘導抵抗が主に圧力に影響を及ぼすことを示しています。
上記結果から3D_CFDにおけるCD0を求めておきます
CD=0.04087＠y=0[m](上右図より)
また誘導速度wi=1.03[m/s]よりαi=tan-1(1.03/26.5)=2.23[deg]
これよりy=0の局所的誘導抵抗は
CDi＝CL*(10-2.23)/10*sin2.23=0.032 (CLは風洞試験値1.063を使用)
∴CD0=CD-CDi＝0.04087-0.0321=0.00877
この値は風洞試験値CD2D=0.00934に十分近い値になっています。
3DのOpenfoamの抵抗係数計算結果は2Dの風洞試験とも辻褄があっています。
更に翼端渦Γと翼中央の循環Γ0の比較を行ってみます。
下図の領域にてQ基準＞0の領域の渦度(vorticity_y)を積分してみました。(Q及びvorticityはOpenFoamを用いた循環Γと渦度の観察にて説明しております。)
エクセルによる積算結果はΓ=12.2[m2/s](翼後縁から8[m]後方)。
これはNACA0012より予想される循環値Γ0(下記)に対して46%(12.2/26.4)になってしまいました。
Γ0=4πURsin10=26.4[m2/s](平板近似)
U=26.5[m/s]
R=6ft/4
これでは小さい過ぎるので、Q基準>0(ひずみよりも回転運動が優勢を示す)という条件を外してみます。
エクセルによる積算結果はΓ=20.2[m2/s]。
これでも循環値Γ0に対して77%(20.2/26.4)と足りません。
風洞試験結果とOpenfoamの計算結果がよく一致している状況では理論値Γ0と計算値Γもおおよそ一致するべきです。
そこで以下の推測をしてみました。
・片翼端渦の吹きおろしの計算式
v1=0.5*Γ0/2π(b/2+y)
この際、翼端渦は半無限のため0.5倍掛けていますが、吹きおろしの影響は前縁だけではなく、翼全体及びある程度後方まで及ぶと考えられます。例えば翼舷長×0.5の距離を影響範囲と考えると0.5倍ではなく、0.5×(1+sin(tan-1(0.5/3))=0.582倍(翼中央で、ビオ・サバール積分範囲の+90～-9degに相当)となります。すなわち0.582/0.5=1.16倍となり、上述77％を89％にまで引き上げられます。Q基準を無視するという無理をしていますが、この値ならば理論値Γ0と計算値Γがおおよそ一致したと言えると考えます。
念のため迎角5degにおいて同じ計算を行ってみました。
結果をグラフに書き足します。
迎角5degにおいても風洞試験値と良い一致をみることが出来ました。
以上3DOpenfoamの計算は教科書で示される誘導抵抗を良く表現出来ており、更に風洞試験値とも良く一致することを確認出来ました。
４．まとめ
・2DのOpenfoam揚力係数計算値は風洞試験の80％とそれなりの値を示しました。
・2DのOpenfoam抵抗係数計算値は風洞試験よりも4倍に増えてしまい、あまり良い一致とは言えない結果となりました(「補足１」にて補足検討実施)。
・一転して3DのOpenfoam計算値と風洞試験は揚力係数、抵抗係数共にと良く一致しています。Openformの計算ロジックは2Dよりも3Dの方が相性がよさそうです。
・付け足し
paraviewでデジタルデータをエクセルにエクスポートするとき、必要なものだけに限定(Toggle column visibility)しないと異常な値をエクスポートすることが分かりました。具体的には17万個というデータの中に元とは違うデータが14個混ざっていました。原因は不明ですが要注意です。
「補足１」
2.2項の計算結果より2Dの計算結果は試験値からかなりずれています(3D計算は試験値と良く一致)。
ここではメッシュ粗さについて再検討してみます。
2.1項のモデルでは想定される乱流境界層21[mm]に対し、メッシュ厚が7[mm]で何とか足りそうに見えます(メッシュ粗さと解析結果の関係に関してはOpenFoamにおける境界層の確認にて検討)。 しかしメッシュの様子(下図)を調べると翼表面のメッシュが一様でなく、境界層を表現するには適していないと感じます。
そこで翼周りの均一なメッシュ層を生成し、計算しなおしてい見ます。
作成したモデルを以下に示します。メッシュ厚は7[ｍｍ]のままですが、より多くの細かい層を100[mm]用意しました。
本モデルを使った抵抗係数Cdの計算結果を以下に示します。
結果は振動してしまい、収束しそうにありません。翼上面の境界層が一定にならず、剥離点が後退してはまた大きく前方に移動するを繰り返しています。実はOpenFoamを用いた2D円柱周りの流体解析detail1における検討でも結果が収束しない経験をしています。境界層剥離は主流と乱流成分が絡み合う非常に繊細な運動であり、それをCFDで扱うこと自体が難しい解析と考えます。
2D円柱回りの流体解析は翼でいうところの迎角０度に対応します。そこで迎角０度の場合を計算してみます。翼周りのメッシュの様子は以下となります。前回同様メッシュ厚7[ｍｍ]の細かい層を100[mm]用意しました。
計算結果を以下に示します。
計算結果は早々に収束し、Cd＝0.00845となりました。この結果を風洞試験値グラフに重ねてみます(緑点)。迎角0[deg]の時の試験値0.00659と大分良く一致しています。
以上よりOpenfoamによる2D試験値の再現は迎角0[deg]では十分一致出来たのですが、迎角が付くと上手くいきませんでした。これは解析対象である翼周り流れが境界層剥離を伴う繊細な流体運動であるためと考えます。
一方でなぜOpenfoamによる3D試験値の再現が十分な一致となったのかの理由として以下を推定しています。
・今回選んだ3D試験値のアスペクト比が6と小さく、翼端渦による吹きおろしの影響で翼上面の境界層は早々剥離が進み、低抵抗を示さなかった。
・この状況とsimpeFoam+k-OmegaSSTの剥離挙動が一致したため。 （Rex=1.8e+6なのに、双方ともに2D円柱でいうとろこのRe<3e+5の領域の挙動となった。)
余談：今回の翼周りメッシュ強化のために2D円柱解析で用いたsnappyHexMeshDictのaddLayer機能により30層足そうとしたのですが、パソコンで24時間計算してもメッシュ作成が終了しませんでした。addLayerは相当な計算資源を消費するようです。
参考文献
[1]航空宇宙工学テキストシリーズ 空気力学入門 日本航空宇宙学会編 丸善出版
[2]OpenFoamによる熱移動と流れの数値解析 第2版 OpenFoamCAE学会編 森北出版
改訂記録
A改訂 2026.10.25 「補足1」追加
