— 12 —SENAC Vol. 49, No. 2（2016. 4）［共同研究成果］

低レイノルズ数流れにおける非定常空力特性の数値解析
児島有哉 1)，飯岡大樹 1)，佐々木大輔 2)
1) 金沢工業大学大学院 工学研究科 機械工学専攻2) 金沢工業大学 工学部 航空システム工学科

本研究では，直交格子法に基づく Building-Cube Method を用いて薄翼周りの低レイノルズ数流れの 3 次元非定常解析を行う．本研究では翼近傍の格子を細かくし境界層内まで直接解析を行っている．解析対象はレイノルズ数 5,000 以下の平板翼，NACA0012 および円弧翼である．解析の結果平板翼と NACA0012 では実験結果と概ね一致した結果が得られたが，曲率変化の大きい円弧翼においては実験結果と乖離した．曲率が大きく変化する形状において，流れ場の変化を直交格子で捉える手法の研究が今後必要である．

1.

緒言
近年，Micro Air Vehicle (MAV)と呼ばれる超小型航空機が災害調査や撮影といった様々な場所で活用されている．MAV まわりの流れ場や空力特性は低レイノルズ数領域（Re≤104）のため，大型旅客機などとは異なっている．低レイノルズ数域における翼周りの流れは，翼前縁付近の流れは層流であるが，剥離泡によって流れは乱流へと遷移し，後方で再付着が生じる等，複雑な流れ場を形成する[1]．このような複雑な流れは翼の空力特性に大きな影響を与えることから，MAV に用いる翼型には空力的な課題が残されている．そのため，複雑な流れ場を理解し，空力性能の改善を行うことで工学的な応用が期待できる．
現在の航空機設計では風洞実験と併用し，数値流体力学（Computational Fluid Dynamics, CFD）が使用されている．計算機の高速化とともに数値解析法の研究が盛んに行われており，現在では大規模な計算を行うため並列計算が一般的になっている．並列計算を効率的に行う CFD 手法として直交格子法に基づく BCM（Building-Cube Method）が中橋らにより提案されている[2]．BCM を用いることで構造格子と同様に高次精度解法を採用でき，また効率的に格子を配置することができる．そこで，本研究ではレイノルズ数 5,000 以下における平板翼，NACA0012 及び円弧翼まわりの CFD 解析を直交格子により行い，実験における空力係数との比較を行う．低レイノルズ数では境界層が厚く，境界層内部まで格子を解像することが可能であるため，直交格子法により低レイノルズ数域における複雑な流れ場の解明を目指す．
2.

Building-Cube Method
解析には非圧縮性流体解析ソルバーのBCM［3, 4］を用いる．BCM とはブロック型直交格子に分類される CFD 手法である．BCMの特徴として，複雑形状に対しても格子生成が容易であること，計算領域を”Cube”と呼ばれる領域に分割し，Cube 内を等間隔・同数の”Cell”で分割することにより，均等な負荷分散を達成できること，また空間高次精度解法を採用できること等が挙げられる．Fig.1 に物体近傍の Cube および Cell を示す．
Fig.1 Mesh allocation of NACA0012
本研究で用いたソルバーの支配方程式は非圧縮airfoil (Cube and Cell definition)
性 Navier-Stokes 方程式と連続の式である．ここで，時間積分法は 2 次精度 Adams-Bashforth 法，空間離散化手法は 1 次精度風上差分法（対流項）および 2 次精度中心差分法（拡散項）を使用する．
物体形状は階段状表現として解析を行う．

— 13 —

低レイノルズ数流れにおける非定常空力特性の数値解析

3.

計算対象および計算条件
黒田らの実験[5, 6]と空力係数値を比較検証するために用いた計算対象は，Fig.2 に示すように翼弦長 c を代表長さとした翼厚 0.02c の平板翼，NACA0012 および翼厚 6%c の円弧翼である．計算領域は翼弦長 c を基準に x，y 方向に 48 倍とする．流入条件は層流とし，3 次元解析を行う．計算領域の両端には周期境界条件を適用し，2 次元性を確保している．計算格子数は約 1,200 万～2,000 万セルである．Fig.3 に使用する物体近傍の計算格子（Cube 分割域），Table 1 に計算格子について示す．

(a) 2% Flat plate(b) NACA0012
Fig.2 Airfoil shape(c) 6% Circular arc

Table1 Computational conditions
Flat plate
NACA0012 Circular arc
3Reynolds number
5.0x102.8x1032.8x103
Min. Cell Size
2.93x10-3
2.93x10-3
5.86x10-3
Cells in a Cube
163
163
163
Number of Cubes
3,112
5,112
5,112
Number of Cells
12,746,752
20,938,752 20,938,752(a) Flat plate(b) Flat plate (near theleading edge)(c) NACA0012(d) NACA0012 (near the leadingedge)
Fig.3 Computational mesh(e) Circular arc(f) Circular arc (near theleading edge)

4. 結果および考察
平板翼，NACA0012 及び円弧翼まわりの 3 次元非定常流体解析を行い，時間平均の空力係数を求めた．平板では迎角 0deg から 20deg まで 2.5deg ずつ変化させ，NACA0012 および円弧翼では0deg から 12deg まで 2deg ずつ変化させ，流れが収束するまでそれぞれ独立して解析を行った．
4.1 平板翼の空力特性
Fig.4 に示すのは，平板翼における実験結果と解析結果の空力係数の比較である．図中の Theoryとは平板翼の揚力係数理論解である．実験結果における平板翼の揚力係数 CL は非線形性を示しているが，解析結果でも同様の傾向を示している．一方，抗力係数 CD は解析結果のほうが小さい値を示したが，傾向はほぼ同じである．しかし，実験結果における最大揚力係数は 12.6deg であるのに対し，解析結果では 20deg まで揚力係数が伸び続けている．この理由として層流解析を行っていることが挙げられる．実験では高迎角で翼の前縁で剥離した流れがすぐに乱流に遷移していると考えられる．一方で，解析においては翼後方においても層流のまま流れているため，揚力— 14 —

SENAC Vol. 49, No. 2（2016. 4）の低下が起きず失速が緩やかになると考えられる．また，本解析では対流項に 1 次精度風上差分を用いているため，粘性が強く剥離が遅れているのではないかと考えられる．
Fig.5 に時間平均の平板翼周りの流線と主流方向の速度分布を示す．平板翼の流れは低迎角から前縁での剥離が見られる．2.5deg における流れ場は翼周りに渦が見られない．5deg になると翼上面に剥離泡が生じ，非定常性が表れている．10deg になると翼上面の渦は大きく成長し，強い非定常性が表れている．また，コード長 50%位置に 2 次剥離渦と呼ばれる反時計回りに流れる渦も確認できる．これは実験においても確認されている．
4.2 NACA0012 の空力特性
Fig.6 に，NACA0012 における実験結果と解析結果の空力係数の比較を示す．実験結果におけるNACA0012 の揚力係数 CL は，やや非線形性を示しているが，解析結果はほぼ線形を示している．
一方で抗力係数 CD は低迎角においては実験結果をやや下回り，高迎角になるほど値が近づいている．
Fig.7 は NACA0012 周りの流線と主流方向の速度分布である．2deg での流れ場は剥離点が翼中央からやや翼後縁よりにあり小さな剥離渦を形成している．6deg になると剥離点はやや前縁方向に移動し後縁に見られる剥離渦が成長している．10deg になると剥離点が翼前縁に移動し翼上面に大きな剥離泡が形成されている．高迎角付近になると揚力傾斜が小さくなることから，この剥離点の移動が揚力係数 CL の非線形性に影響を与えているものと考えられる．NACA0012 周りの流れは剥離した後再付着することのない流れとなっている．
4.3 円弧翼の空力特性
Fig.8 は，円弧翼における実験結果と解析結果の空力係数の比較を示している．また，Fig.9 は，円弧翼周りの流線と主流方向の速度分布を表している．Fig.8 より，実験結果では円弧翼の揚力係数 CL はほぼ線形の傾向を示しており，解析結果においても低迎角付近でやや値が大きいが，実験結果と同じような傾向を示している．しかし，高迎角になると実験値との乖離が見られる．この理由としては Fig.9 に示すように，2deg および 6deg では最大キャンバー高さ付近で流れが剥離するのに対し，10deg では流れが前縁から大きく剥離しており，翼上面に現れる大きな渦の影響であると考えられる．しかし，Fig.3 (c)のように翼上面の格子は翼近傍と比較して粗くなっているた
め流れが正確にとらえられていない可能性がある．解析結果の抗力係数 CD は低迎角付近では実験結果とほぼ一致している．しかし，同様に高迎角になると実験結果の傾向からそれている．
Fig.9 より，2deg における流れ場は強い非定常性が表れている．しかし，迎角が大きくなると非定常性は弱くなる．低迎角時には後縁付近に渦が形成されているが，高迎角になると剥離点が前縁に移動し渦は翼上面全体に形成されている．このことから高迎角になると急激に揚力係数が増加しているのではないかと考えられる．
4.4 翼型の空力特性比較実験において平板翼と NACA0012 は揚力係数に非線形性を示し，円弧翼は 2 つの翼型に比べ線形性を示す．また，厚翼である NACA0012 は薄翼である平板翼や円弧翼に対して，揚力係数は低い値を示す．解析結果では，平板翼では 12deg まではその特徴を捉えることができた．NACA0012に対しては 6deg までは実験の傾向を捉えることができているが，それ以降は傾向からそれている．
円弧翼は全体的に実験の傾向を捉えることができていない．NACA0012 と円弧翼において傾向を捉えることができなかった理由は翼近傍では格子が密であるが，そのまわりの格子は粗くなっているため，迎角が増すと成長する渦を正確に捉えることができていない可能性がある．また，これら 2 つの翼型は曲面形状であり，直交格子を用いた階段状表現では，物体形状を正確に表現できていないこともその理由として考えられる．
Fig.10 に示すのは NACA0012 と円弧翼における z 軸方向の等渦度面の上の速度分布である．2degでの NACA0012 における z 方向の速度変化はあまり見られていない．しかし，10deg になると 3

低レイノルズ数流れにおける非定常空力特性の数値解析— 15 —次元的な渦が発生し，z 方向の速度に変化が表れているのがわかる．それに対し，円弧翼では 2degと 10deg のどちらでも z 方向の速度変化はあまり見られない．NACA0012 と円弧翼では翼の高さは同じであるがキャンバー位置が NACA0012 のほうが前方にある．そのため剥離点が前方にあり，剥離の範囲も大きく強い非定常な流れとなっていると考えられる．

5. 結言
本研究では，BCM を用いて低レイノルズ数流れにおける平板翼，NACA0012 および円弧翼周りの 3 次元非定常流体解析を行い風洞実験値と比較した．今回，平板翼においては高迎角時を除き，実験値と同様の傾向を捉えることができた．しかし，NACA0012 では 6deg 程度まで実験値の特徴を捉えることができる一方，それ以降の迎角では実験値からやや逸れる結果となった．また，円弧翼では全体的に実験値との差があり，その特徴を捉えることができなかった．その原因として翼の近傍では格子が細かいのに対し，周辺では粗くなっていることが挙げられる．また，翼近傍においても階段状表現であるため物体近傍の流れが正確に捉えられていない可能性がある．
低レイノルズ数域における複雑流れに対して直交格子法の優位性を示すためには，今後，翼近傍全体の流れの変化が大きい領域に対して細分化格子を適用する等，より大規模な格子による大規模流体解析が必要となる．将来的に，直交格子法に基づく CFD 解析手法の信頼性が実証できると，現在注目されている超小型航空機やマルチコプタの運用時間向上に適用可能である．

— 16 —

SENAC Vol. 49, No. 2（2016. 4）

0.7

Theory0.6

CL simu.
CD simu.

0.5

CL exp.

CL,CD0.4

CD exp.

0.3
0.2
0.1
0-0.1

Fig.4 Aerodynamic characteristics of flat plate

-113
5
7
9Angle of Attack, α, deg1113

Fig.6 Aerodynamic characteristics of NACA0012(a) 2deg(a) 2.5deg(b) 6deg(b) 5deg(c) 10deg(c) 10deg

Fig.5 x-direction velocity distributions andstreamlines NACA0012

Fig.7 x-direction velocity distributions andstreamlines around flat plate

— 17 —

低レイノルズ数流れにおける非定常空力特性の数値解析(a) 2deg (NACA0012)

Fig.8 Aerodynamic characteristics of circular arc(b) 10deg (NACA0012)(a) 2deg(c) 2deg (Circular arc)(b) 6deg(d) 10deg (Circular arc)
Fig.10 z-direction velocity distributions on
z-direction vorticity isosurface(c) 10deg

Fig.9 x-direction velocity distributions andstreamlines around circular arc

— 18 —

SENAC Vol. 49, No. 2（2016. 4）

謝辞
本研究は，東北大学サイバーサイエンスセンターのスーパーコンピュータを利用することで実現することができた．また，研究にあたっては同センター関係各位に有益なご指導とご協力をいただいた．
参考文献
1] 近藤勝俊，“低レイノルズ数で優れた空力性能を持つ翼型設計に向けた CFD を用いたパラメトリック解析”，第 52 回飛行機シンポジウム講演論文集，長崎，(2014)．
[2] K. Nakahashi, and Kim, L. S.,“Building-Cube Method for Large-Scale, High Resolution FlowComputations,” AIAA Paper 2004-0423, (2004).
[3] R．Sakai, “Study of Practical Large-Scale Turbulent Flow Simulation Method Using Cartesian Mesh,”
Ph.D Thesis，Tohoku University, (2013).
[4] 坂井玲太郎，大林茂，松尾裕一，中橋和博，“Building-Cube Method を用いた実用的な大規模乱流解析手法の構築”， 第 45 回流体力学講演会/航空宇宙数値シミュレーション技術シンポジウ
ム 2013，船堀，(2013)．
[5] 黒田達哉，岡本正人，“低レイノルズ数翼型フェザリング応答特性”，第 44 回日本航空宇宙学会年会講演会講演集，東京，(2013)．
[6] T. Kuroda and M. Okamoto，”Unsteady Aerodynamic Forces Measurements on Oscillating Airfoils
with Heaving and Feathering Motions at Very Low Reynolds Number”, Proceedings of APISAT,
Takamatsu, (2013).