第 25 回数値流体力学シンポジウム
A04-1

平均揚力最大時の羽ばたき翼の三次元効果の解析
Analysis of three-dimensional effects on flapping wing aerodynamics associated with atime-averaged lift maximized wing motion○ 宇賀神誠也, 東理大院,
東京都千代田区九段北 1-14-6, E-mail:ugajin@flab.isas.jaxa.jp
青野光,ISAS/JAXA, 神奈川県相模原市由野台 3-1-1, E-mail:aono@flab.isas.jaxa.jp
野々村拓,ISAS/JAXA, 神奈川県相模原市由野台 3-1-1, E-mail:nonomura@flab.isas.jaxa.jp
大山聖,ISAS/JAXA, 神奈川県相模原市由野台 3-1-1, E-mail:oyama@flab.isas.jaxa.jp
藤井孝藏,
ISAS/JAXA, 神奈川県相模原市由野台 3-1-1, E-mail:fujii@flab.isas.jaxa.jp
山本誠,
東理大,
東京都千代田区九段北 1-14-6, E-mail:yamamoto@rs.kagu.tus.ac.jp
Seiya UGAJIN, Tokyo Univ. of Science, 1-14-6 Kudankita, Chiyoda-ku, Tokyo, JAPAN
Hikaru AONO, ISAS/JAXA, 3-1-1 Yoshinodai, Sagamihara, Kanagawa, JAPANTaku NONOMURA, ISAS/JAXA, 3-1-1 Yoshinodai, Sagamihara, Kanagawa, JAPANAkira OYAMA, ISAS/JAXA, 3-1-1 Yoshinodai, Sagamihara, Kanagawa, JAPAN
Kozo FUJII, ISAS/JAXA, 3-1-1 Yoshinodai, Sagamihara, Kanagawa, JAPANMakoto YAMAMOTO, Tokyo Univ. of Science, 1-14-6 Kudankita, Chiyoda-ku, Tokyo, JAPANAn aerodynamic design problem of a flapping airfoil/wing in low Reynolds number (Re=103) is explored with themulti-objective design exploration framework coupled with a Navier-Stokes solver. It is found that a tradeoff
relationship among lift maximization, thrust maximization, and required power minimization for forward flight
condition. Qualitative comparison with three-dimensional and two-dimensional flapping wing is conducted tounderstand the effect of dimensionality of the flow field and wing motion. It can be seen that formation anddevelopment of leading-edge vortex (LEV) is the different between three-dimennsional wing and two-dimensional
wing as well as the LEV on three-dimensional wing generates more lift than that on two-dimensional wing. On theother hand, wing tip vortex largely decrease the lift.
１．序論
近年，ローバーやランダーのような着陸機ではなく，飛行タイプの火星探査機（火星航空機）が様々な国で提案され[1-3]，注目を集めている（Fig.1）．しかし，火星大気の密度は地球に比べ非常に薄く，想定される火星航空機の飛行条件ではレイノルズ数（Re数）が非常に小さくなる．一般的な固定翼タイプや回転翼タイプの航空機翼の空力特性は Re 数の低下と共に悪化することが知られており，想定される火星航空機の Re 数（103～104 程度）では，流れは容易に翼面からの剥離を起こしやすく，高い揚力を出すことが難しい．従って，火星航空機には高い揚力を確保するための何らかの工夫が必要である．
現在，火星航空機システムとして考えられているのが，固定翼
機・回転翼機・羽ばたき機・気球などである．それぞれの飛行システムには長所・短所があり，この内，羽ばたき機は以下の利点から期待されている．
・ 想定される Re 数が 103 程度で，地球上での昆虫のそれとおおよそ一致しており，昆虫などが低 Re 数下で羽ばたき飛行を行う事実から，火星のような低 Re 数でも高い揚力の発生や優れた飛行性能を発揮できる可能性がある・ 火星という未開の荒地において，比較的狭い領域での垂直離着陸が可能である
・ 固定・回転翼に比べて突風や横風に強く，嵐などが多いとされる火星環境下でも対応できる
しかしながら，羽ばたき運動の空力メカニズムが十分に明らかになっていないため，自律制御が難しいという問題点がある．
羽ばたき翼の空力性能向上を図るためには，想定される火星羽
ばたき航空機の Re 数に於いて，特徴的な羽ばたき運動である，・ ホバリング時の運動
・ 巡航時の運動の 2 つの空力メカニズムの理解が不可欠であると考えられる．

Fig. 1 Examples of proposed Mars explorer
ホバリング時の運動については，これまでに地球上の昆虫をモデルにした多くの研究がなされており，昆虫はホバリング時に流れを翼前縁で剥離させ，生じる強い前縁剥離渦の負圧を利用することで高い揚力を得ていることが実験や数値計算を用いた研究により明らかになっている[4-6]．
一方巡航時の運動については， MAV（Micro Air Vehicle）開発に向けた羽ばたきの研究の中で行われている．これらの研究での
Re 数は 104～105 であり，高い推進効率を出すような羽ばたき運動
は，ホバリング時にみられた前縁剥離渦を発生させるといった大規模剥離を伴わない運動であるという報告がなされている[7,8]．しかし，火星航空機にこれらの知見を適用するには以下の問題点がある．
・ 火星航空機の想定される 103 程度の Re 数での巡航時に於ける運動についての研究が十分されていないため，この領域での揚力・推力の発生メカニズムの理解が不十分である．
・ 多くの MAV 設計に向けた 3 次元羽ばたき運動の流れ場に関する研究は，主な運動の定義として、翼の上下運動のみを考慮しているが，鳥や昆虫が利用している実際の羽ばたき運動には，翼の迎角を時間的に変化させるような複雑な運動であり，このような運動に対しての議論が不十分である．
1

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
A04-1
岡部ら[9]は火星航空機を想定したレイノルズ数領域中で2次元的な羽ばたき運動をする翼について数値解析を行った．また，最適化手法[10]と組み合わせ，数値計算結果から得られる翼の空力性能に優劣をつけることにより，揚力を増大させる羽ばたき運動や，推力を増大させる羽ばたき運動についてそれぞれの流れ場に対して議論を行った．岡部らの研究結果から，103 程度の Re 数での 2
次元的な羽ばたき運動では，揚力を大きくしたい場合はダウンス
トローク（打ち下ろし運動）時に強い前縁剥離渦を発生させる必要があり，推力を大きくしたい場合はダウンストローク・アップ
ストローク（打ち上げ運動）が対称の運動にすることによる渦の発生が垂直力の推力に大きく寄与することが分かっている．
これらの研究で得られた知見は非常に有用であるが，3 次元的な羽ばたき運動に関しても当てはまるという確証はなく，翼性能に寄与している流体現象の理解にむけて議論した研究が不足している．そこで本研究では，3 次元羽ばたき運動を上下運動とピッチ運動を組み合せた運動を用いて定義し，過去の研究[9]と同様に最適化手法と数値計算を組みわせた解析を行う．本解析により効率的に揚力や推力を増大させる羽ばたき運動を分類し，複雑になり得る翼運動に於いて発生する 3 次元的な翼周りの渦の構造や発達の翼性能への影響を議論する．さらに、流れ場の三次元性の影
響の詳細な議論を行うために 2 次元羽ばたき翼との定量的な比較を行う．

Plunging motion




・



2




  02

yx

Pitching motion
Leading edge is fixed

1

z

0

x

Fig. 2 Schematic of flapping motion

k

2 fc
U

(4)

である．これらの式の中の変数 k, h,  0 ,  1 , を適当に変化させる
ことにより羽ばたき運動を実現する． flap はフラッピング角を表しており，f が最大のプランジ角である．k は式(4)で表される羽ばたきの周波数 f を無次元化したものであり，h はコード長 c 無次元化された翼端での上下運動振幅の大きさである．0 はオフセット角であり，1 は翼前縁を中心に回転させるピッチ運動の変化を表すピッチ角である．は上下運動とピッチ運動の位相差を示している．U∞は一様流の速度であり，t は一様流速度とコード長で無次元化された時間である．下記に設計変数/設計空間を示す．
【設計変数】
・ 羽ばたきの無次元周波数， k・ 上下運動（プランジ運動）の振幅， h・ ピッチ運動の振幅，1・ ピッチ運動のオフセット，0
・ ピッチ運動とプランジ運動の位相差，
【設計変数空間】
・ k : 0.2 ~ 1.0
・ h : 0.5 ~ 2.5
・ 1 : 15°~ 50°・ 1 : 0° ~ 40°
・  : 30°~ 270°
最後に，数値解析に用いた流れ場条件を示す．本研究では，扱う翼の寸法を Michelson らが NASA に提案した火星探査ミッショ
ン[3]を参考にした．翼のコード長を 0.1 [m]，スパン長を 0.5 [m]とし，火星を約 10 [m/s]の速度で巡航している状態を想定している．解析に用いた Re 数は 103 であり，翼のコード長を代表長さとしている．翼型には，低 Re 数に於いても空力特性が良いといわれている[11]薄翼 NACA0002 を用いた．翼平面は矩形となっている．

(2)



z
root

Airfoil shape: NACA 0002

ピッチ運動

 pitch  1 sin kt   

tip

lspan

(1)

 flap   f sin kt fh

２．問題設定(1) 設計問題の定式化
本研究の目的は，火星探査用の羽ばたき機に適した羽ばたき運
動を探すことである．火星羽ばたき機に適した性能を出す運動とは，高い揚力係数を出す，必要パワーを小さく抑える，また推力係数を大きくする運動のことである．ここでいう揚力係数，推力
係数は羽ばたき運動 1 周期分（ダウンストロークとアップストロークを 1 回ずつの運動）を時間平均した値であり，推力は翼の抵
抗分を差し引いた正味の推力のことである．そこでこの運動 1 周
期分を時間平均した揚力係数，推力係数，および必要パワーを無次元化したもの（必要パワー係数）の 3 つを目的関数として最適化を行った．以下に多目的最適化問題の定義を述べる．
【目的関数】
・ 時間平均揚力係数 CL,ave の最大化
・ 時間平均推力係数 CT,ave の最大化
・ 時間平均必要パワー係数 CPower,ave の最小化
【制約条件】・ CL,ave > 0
・ CT,ave > 0
また以下のような羽ばたき運動を定義し，設計変数/設計変数空間を定義した．まず，上下運動とピッチ運動を組み合わせた 3 次元
羽ばたき運動（Fig. 2）は，双方の運動を次に示す sin 関数によって定義する．以下にこの 2 つの運動を表現した式を示す．
・ プランジ運動（上下運動）h f  sin 1   ls 

Symmetry Plane

(3)

ここで，
2

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
A04-1
また，翼は剛体翼を仮定している．
(2) 数値計算手法
基礎方程式は 3 次元圧縮性 Navier-Stokes 方程式を用いた．対流
項の離散化には 3 次精度 MUSCL[12]法を用い，粘性項の離散化に
は 2 次精度中心差分を用いた．時間積分には 1 次精度 ADI-SGS
陰解法[13]を用いた．本研究では低 Re 数であることに加え，乱流に伴って発生する細かな渦を解く必要は無く，剥離により生ずる
比較的大きな渦構造を捉えればよいと考えたため，乱流モデルは使用せず層流計算とした．
本研究で用いた計算格子を Fig. 3，4 に示す．計算格子には Fig.
3 に示すような H 型格子を用いた．格子点数は 201×100×100 の合計約 200 万点である．翼の前後上下の外部境界はコード長の 20倍とした．本研究では，3 次元羽ばたき運動特有の現象であるス
パン方向流れが重要な渦構造の 1 つであると考え，翼スパン方向の外部境界は，コード長の 50 倍をとっている．また，翼周りの渦を詳細に捉える必要があることから，Fig. 4 に示すように翼面付
近に格子を集中させ，更には翼端部で密となるような格子分布とし，十分な格子解像度を保った．

c

３．結果及び考察(1) 最適化結果
初めに目的関数間の傾向を考察するため，最適化によって得られた個体に関する情報を基にデータマイニングを行った．本研究では世代数 30，総個体数 496（初期個体を含む）のうち非劣解となるものが 198 であった．ここで多目的最適化によって得られた
3 つの目的関数のうち非劣解となる個体の関係を示した散布図をFig. 5 に示す．全体に以下の傾向があった．
① 高揚力係数を得る羽ばたき運動では推力はほとんど発生しない② ①とは反対に高い推力係数を得る羽ばたき運動の場合，ほとんど揚力が発生しない③ 低いパワーでの運動では揚力も推力もほとんど発生しない④ 最大揚力よりも最大推力を発生させる運動の方がよりパワーが必要である
以上のことから，3 つの目的関数にはトレードオフの関係があるといえる．本研究では最大揚力係数の得られた運動に着目し，解析を進めていく．
CLmax,ave solution

z
x

20 c

2.5

y

0.156

CPower

0.8

1.0

1.046

20 c

5c

2.0
1.5

50 c

c

CL
1.0
0.5

Fig. 3 Computational grid

0.0
0.0

0.2

0.4

0.6

1.2

1.4

CT

TE

CTmax,ave solution
Fig. 5 Scatter Plots

75 point
LE
Tip

50 point

Table 1 に最大時間平均揚力係数の解の設計変数を示す．以下，
本論ではこの設計変数を用いた 2 次元羽ばたき翼の解析を行い，3次元の空力特性と比較する．

Root

Fig. 4 Computational grid (top view)(3) 最適化計算手法
パラメトリックな解析から羽ばたき運動の性質を把握するのは困難であると考え，ここでは多目的遺伝的アルゴリズムを用いて最適化を行った．多目的遺伝的アルゴリズムには Deb によって提案された，Non-dominated solution genetic algorithm-II (NSGA-II)[14]を用いた．NSGA-II では，個体間の優越関係[15]と混雑距離を用いたバイナリトーナメント選択により親個体ペアが選択される．選択された親個体ペアに実数値交叉手法である SBX[16]と多項式突然変異手法である Polynomial mutation[17]が適用され，子個体が生成される．親個体と子個体を組み合わせた母集団から非優越ソー
トと混雑距離を用いて次世代の親個体が選ばれる．
本研究では，個体数 16，世代数 30，交差率 0.7，突然変異率 0.1として最適化をおこなった．

Table 1 Wing kinematics for lift maximization of a flapping wing
k [-]
h [-]
0 [deg.] 1 [deg.] deg.
1.00
22.0
45.7
102
2.48

(2) 羽ばたき運動中の翼性能の推移
Fig. 6 に本解析で得られた 3 次元羽ばたき翼と 2 次元羽ばたき
翼の 1 周期あたりの揚力係数の推移を示す．ここで，3 次元羽ばたき翼は翼面の平均揚力であり，2 次元羽ばたき翼の運動は 3 次
元羽ばたき翼の翼端と同じ運動をする．Fig. 6 より，2 次元と 3 次
元共に，
主にダウンストローク中に揚力を発生させる運動である．
無次元時間 0 から 0.35 までの間および，0.9 から 1.0 までの間は 2
次元と 3 次元共に似たような揚力係数の時間変化を示しているが、無次元時間約 0.35 から 0.9 の間に 2 次元と 3 次元の揚力係数が最大になる羽ばたきのタイミングが異なっており，また 2 次元羽ばたき翼の無次元時間 0.8 において 3 次元羽ばたき翼にはみられな
い揚力係数の急激な低下がみられた．
3

Copyright © 2011 by JSFM

22.0

3D CL
pitch

18.0

14.0

2D CL
plunge

110
1

90

2

70
50

6.0

3

30

2.0

10

-2.0

-10

-6.0

Angle [deg.]

10.0
CL

第 25 回数値流体力学シンポジウム
A04-1
全面での平均揚力は正に保たれている要因であると考えられる．
(3) スパン方向の揚力分布
・ フェーズ 1 とフェーズ 2 の比較
フェーズ 1（2D flapping CL max）およびフェーズ 2（3D flapping
CL max）における 2 次元羽ばたき翼と 3 次元羽ばたき翼の揚力係数が最大値になる羽ばたきのタイミングの違いに関して，翼周りの 3 次元性の渦が揚力の発生にどのように寄与しているかをより詳細に議論するために，Fig. 11 に 3 次元翼のスパン方向の揚力係数の分布を示す．Fig. 11 から全体として以下の２つのことがあげ
られる．
I) 75%スパン付近から翼端までの各断面の揚力係数はフェー
ズ 1 の方がフェーズ 2 よりも高いII) 逆に翼根から 75%スパン付近にかけては，フェーズ 1 よりもフェーズ 2 の方がより揚力係数が高い
3 次元羽ばたき翼の場合は翼全面での平均をとるため，翼全面での揚力はフェーズ 1 よりも 2 の方が高くなり，2 次元羽ばたき翼の最大揚力係数のタイミングと異なることがわかる．
・ 翼端付近の 3 次元性の検証
フェーズ 2 における 3 次元羽ばたき翼の 3 次元性とみられる特性を議論するために，3 次元羽ばたき翼の翼根から翼端までの10%スパン刻みでの断面と同じ運動をする 2 次元羽ばたき翼を解析し，比較する．Fig. 12 にフェーズ 2 の時のそれぞれのスパン方向の断面の 2 次元羽ばたき翼と 3 次元羽ばたき翼の揚力係数の分布を示す．
翼根から 60%スパン位置にかけては，3 次元翼の各スパン位置
での揚力係数がその位置の翼断面と同じ運動をする 2 次元羽ばたき翼の揚力係数よりも高い値であるのに対して，70%スパン位置
から 2 次元羽ばたき翼の方が 3 次元羽ばたき翼よりも高い値を示している．特に，80%スパンで急激な揚力の低下が発生しており，90%スパン付近から翼端の手前においては翼端渦の影響から少し揚力は回復するが 2 次元羽ばたきの揚力係数よりは下回っている．
ここで，Fig. 13～16 にフェーズ 2 における 3 次元羽ばたき翼の60%スパンから翼端までの４つの翼断面の翼周りの圧力係数分布
コンターとその断面と同様の運動をする 2 次元羽ばたき翼周りの圧力係数分布コンター及び，各断面における翼面上の圧力係数の分布を示す．
① 90%スパン位置
90%スパンにおいて，Fig. 13から2次元と3次元の違いは主に，翼下面の高圧領域と翼上面の後縁付近の低圧領域である．90%ス
パンは翼端に近いため，翼端から発生する渦の影響を受けて，翼下面の高圧領域が上面に移動する．これにより下面の高圧領域が2 次元羽ばたき翼よりも狭くなっていると考えられる．次に，Fig.
13 から2 次元羽ばたき翼は後縁渦による低圧領域が存在していることがわかるが，3 次元の場合，Fig. 9 でみられるような翼端から発生する渦が存在しており，後縁渦の効果がほとんど得られないと考えられる．これらのことから，90%スパンの位置では 3 次元
羽ばたき翼よりも 2 次元羽ばたき翼の方がより揚力が高い．
② 80%スパン位
80%スパンの位置における 2 次元羽ばたき翼も 90%スパン位置と同様に Fig. 14 から後縁に存在する渦によって揚力を得ている
ことがわかる．その一方で，3 次元羽ばたき翼の場合は，翼上面で大規模な剥離を起こすことが Fig. 14 からわかる．これにより翼上面に低圧領域が存在せず，Fig. 14 のような圧力分布となり，揚力の低下を招くことがわかる．

-30

-10.0

Up stoke

Down stroke

Up stroke

-14.0

-50

-70
0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
t* and wing motion for one
Fig. 6 Time histories of lift coefficientmotion cycle
ここで 2 次元羽ばたき翼の無次元時間 0.50，0.55，0.60，0.75，
0.80，0.85.の翼周りの流れ場をみるために各時間における圧力係
数分布のコンターを Fig. 7 に示す．また Fig. 6 でのフェーズ 1 ～3 における3 次元羽ばたき翼の流れ場を議論するために翼周りの圧力係数分布のコンターおよび翼表面の速度勾配テンソル第 2 不変量の等値面（Q=1）を用いて渦構造を表したもの（ｙ方向渦度で面塗り）を，Fig. 8～10 にそれぞれ示す．
・ T= 0.50 (フェーズ 1 )フェーズ 1 は 2 次元羽ばたき翼の揚力係数が最大になる時間である．Fig. 7 から，2 次元羽ばたき翼は上下運動と迎角の関係により下面に大きな高圧領域が存在している．また，翼上面には後縁付近に渦が存在する状態であり，これらが揚力を最大にする要因と考えられる．その一方で，Fig. 8 より 3 次元羽ばたき翼上面の渦は全体的に 2 次元羽ばたき翼の場合と比べて弱いので 2 次元よりも揚力が発生していないといえる．この渦が弱くなる現象は，3次元羽ばたきの特徴的な翼端から発生する渦と 3 次元的な前縁渦の発達が影響していると考えられる．
・ T= 0.55 (フェーズ 2)
フェーズ 2 は 3 次元羽ばたき翼の揚力係数が最大になる時間である．2 次元の場合は Fig. 7 から翼上面に存在していた渦が翼の後縁に移動して翼から離れている．これがフェーズ 1 と比べて揚力が低下した原因と考えられる．一方で，Fig. 9 をみると 3 次元
羽ばたき翼の翼端付近では 2 次元と同様に翼上面から渦が完全に離れている様子がみられるが，翼根から 40%スパンにかけてはまだ強い前縁渦が翼上面に存在している．このことから，この時間が最も 3 次元羽ばたきの揚力を稼いでいると考えられる．Fig. 7の 2 次元羽ばたき翼の流れ場の特徴は各スパンのどの断面とも一致しないため，詳細な流れ場の議論は次の章で行うことにする．
・ T= 0.80 (フェーズ 3)
フェーズ 3 は 2 次元羽ばたき翼の揚力係数が最小になる時間である．この谷は 3 次元羽ばたき翼では見られない．Fig. 7 から 2
次元羽ばたき翼の場合，翼上面に高圧領域が存在していることがわかる．これにより 2 次元羽ばたき翼では揚力係数が負になる．
その一方で，3 次元羽ばたき翼では Fig. 10 より翼端からの渦の影
響により翼上面の高圧領域が存在していない．また，40%スパンから翼根において未だ翼上面に渦が存在している．つまり，翼根から翼の半分付近までは前縁渦により正の揚力が発生しており，
翼の半分から翼端では高圧領域が存在しないため 2 次元羽ばたき
翼ほどの著しい揚力の低下が見られないことが 3 次元羽ばたき翼

4

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
A04-1

U∞
t*= 0.50(フェーズ 1)U∞t*= 0.60

t*= 0.55(フェーズ 2)-20

t*= 0.75

Cp

t*= 0.85

t*= 0.80(フェーズ 3)15

Fig. 7 Snapshots of pressure distribution around a 2D flapping airfoil
t*=0.50
6
Wing Tip Vortex

Leading Edge Vortex
Root

20% span

40% span
Cp

U∞
U∞

Isosurface is colored by y-vorticity
Isosurface of 2nd invariant of the velocity gradient tensor

60% span

80% span-12

Tip

Fig. 8 Snapshots of vortical structures and pressure distribution of a 3D flapping wing (phase1)
t*=0.55
6Root

20% span

40% spanCp

U∞

Isosurface is colored by y-vorticitynd

Isosurface of 2 invariant of the velocity gradient tensor

60% span

80% span

Tip-12Fig. 9 Snapshots of vortical structures and pressure distribution of a 3D flapping wing (phase2)
t*=0.80
6Root

20% span

40% spanCp

U∞

Isosurface is colored by y-vorticity
Isosurface of 2nd invariant of the velocity gradient tensor

60% span

80% span

Tip-12Fig. 10 Snapshots of vortical structures and pressure distribution of a 3D flapping wing (phase3)

5

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
A04-1③ 70%スパン位置
3 次元羽ばたき翼の70%スパンの位置ではFig. 15 より80%位置で完全に剥離していた前縁渦が存在しており，上面低圧領域が大きくなることがわかる．さらに Fig. 15 から翼下面の圧力分布が 2
次元羽ばたき翼の下面の圧力分布とほぼ同等になることがわかる．
翼端から離れるにつれて下面の高圧領域が 2 次元と同等になっていることから，90%スパンや 80%スパンでみられるような 2 次元
羽ばたき翼より小さく見積もられる翼下面の高圧領域の原因は翼端から発生する渦の影響であると考えられる．
④ 60%スパン位置
スパン60%の位置では逆に2次元羽ばたき翼で渦が翼面から離れており，翼上面の低圧領域が小さくなっている．一方で，3 次
元羽ばたき翼では前縁渦が翼上面に存在しており大きな低圧領域を形成することが Fig. 16 からわかる．この 3 次元羽ばたき翼上面の渦の挙動は 2 次元羽ばたき翼と異なることから，翼端から発生する渦やスパン方向に流れを持つ 3 次元的な前縁渦が翼根から60%スパンでは揚力を増加させるように促したと考えられるが，現在検討中である．

12.0

CL

8.0

2.0
0.0
0 10 20 30 40 50 60 70 80 90 100(Root)(Tip)
Distance of span direction [%]Fig. 11 Spanwise lift distribution

12.0

CL

8.0
6.0
4.0
2.0
0.0
0 10 20 30 40 50 60 70 80 90 100(Root)(Tip)
Distance of span direction [%]Fig. 12 Comparison of spanwise lift distribution between
two-and three-dimensional simulation at phase 2-18.0

3D Cp (90% span)
2D Cp (90% span)

3D Cp (80% span)
2D Cp (80% span)

-14.0

-10.0

-10.0

-6.0
Upper surface

-2.0

2D

Cp

Cp

3D CL
2D CL

10.0

-18.0

-6.0

2D

-2.0

2.0

2.0

Lower surface

6.0

6.0
3D

10.0

3D
0.2 0.4 0.6 0.8
1
-14
Cp
7
Distance in chord direction
Fig. 15 Comparison of two and three dimensional wing: pressure
coefficient(left); pressure contours(right). Black arrow indicatesdirection of wing motion. (70% span)

10.0

0

0.2 0.4 0.6 0.8
1
-14
Cp
7
Distance in chord direction
Fig. 13 Comparison of two and three dimensional wing: pressure
coefficient(left); pressure contours(right). Black arrow indicatesdirection of wing motion. (90% span)

0

-18.0

-18.0
3D Cp (70% span)
2D Cp (70% span)-14.0

3D Cp (60% span)
2D Cp (60% span)

-14.0
-10.0

-10.0
-6.0

2D

Cp

Cp

6.0
4.0

以上のことから，
羽ばたき運動のダウンストローク中に際して，3 次元性の前縁渦や翼端から発生する渦の流れが 3 次元羽ばたき
翼の 70%スパン付近から翼端における揚力を 2 次元羽ばたき翼と比べて低下させ，翼根から 70%スパン付近までの揚力を 2 次元羽
ばたき翼と比べて大きく増加させることに寄与していることが分
かった．

-14.0

3D CL (フェーズ 1)
3D CL (フェーズ 2)

10.0

-6.0

-2.0

-2.0

2.0

2.0

6.0

6.0

2D

3D
0.2 0.4 0.6 0.8
1
-14
Cp
7
Distance in chord direction
Fig. 16 Comparison of two and three dimensional wing: pressure
coefficient(left); pressure contours(right). Black arrow indicatesdirection of wing motion. (60% span)

3D
0.2 0.4 0.6 0.8
1
-14
Cp
7
Distance in chord direction
Fig. 14 Comparison of two and three dimensional wing: pressure
coefficient(left); pressure contours(right). Black arrow
indicates direction of wing motion. (80% span)

10.0

10.0

0

0

6

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
A04-1
４．結言
上下運動とピッチ運動を組み合わせた羽ばたき運動における空力設計を 3 次元 CFD 解析と最適化手法を組み合わせて行った．
また，揚力を最大にする羽ばたき運動の渦の挙動とその翼性能への影響について調査を 2 次元羽ばたき運動と比較することにより行った．本研究で得られた知見は以下の通りである．
・ ダウンストローク中において，2 次元羽ばたき翼と３次元羽
ばたき翼では最大揚力が発生するタイミングが異なる
・ 2 次元羽ばたき翼は打ち上げ初めに大きく負の揚力を発生す
るが，3 次元羽ばたき翼では発生しない
・ 2 次元羽ばたき翼が最大揚力係数を得ている時間の 3 次元羽
ばたき翼について 3 次元羽ばたき翼は翼根から 70%スパンにかけての前縁渦がまだ完全に発達していないので，揚力係数が最大値をとっていない・ 3 次元羽ばたき翼が最大揚力係数を得ている時間の 3 次元羽
ばたき翼について 3 次元羽ばたき翼の 70%スパンから翼端までは揚力が低下しているが，翼根から 70%スパンまでは前縁渦が発達することで揚力が増加しているので揚力係数が最大値となっている 80%スパンにおいて，翼端から発生する渦の影響による翼下面の高圧領域の縮小と翼上面の大規模な剥離により 2
次元羽ばたき翼に比べて大幅な揚力低下が起こっている 翼根から 70%スパン付近において，同条件で計算した 2
次元羽ばたき翼では渦が剥離しているのに対し，3 次元羽
ばたき翼では前縁渦が翼面上に存在しているため，2 次元
羽ばたき運動よりも高い揚力を得ている
以上の様に本解析では，2 次元羽ばたき運動には見られない 3 次
元羽ばたき運動の流体的特性について示した．これにより，スパン方向の流れは揚力発生と流れ場の 3 次元性に大きく寄与していることが分かった．本研究では最適化によって得られた揚力を最大にする羽ばたかせ方の結果についてだけの議論であったので，今後は推力を最大にする運動や必要パワーを最小にする運動を比較し，それぞれの特徴に合わせ，流れ場と翼性能の関係を調査する予定である．

参考文献(1)
Guynn, M. D., Croom, M. A., Smith, S. C., Parks, R. W. andGelhausen, P. A., ”Evolution of a Mars Airplane Concept for the
ARES Mars Scout Mission,” AIAA Paper 2003-6578, 2003.
(2)
田中義輝，岡部能幸，鈴木大晴，中村久美子，久保大輔，
徳弘雅世，李家賢一，地質・地形探査用火星航空機の概念設計について，
日本航空宇宙学会第 36 期年会講演会講演集（2005）, 61-64.
(3)
Robert C. Michelson, Messan A. naqvi, Extraterrestrial Flight(Entmopter-Based Mars surveyor ), Low Re Aerodynamics on
Aircraft Including Applications in Emerging UAV Technology
RTO-AVT von Karman Institute for Fluid Dynamics Lecture
Series Nov. 24-28,2003.
(4)
Weis-Fgh, T., “Quick Estimate of Flight Fitness in Hovering
Animals, Including Novel Mechanism for Lift Production,”
Journal of Experimental Biology, Vol. 59, 1973, pp. 169-230.
(5)
Ellington, C. P., Usherwood, J. R., “Lift and Drag Characteristicsof Rotary and Flapping Wings,” edited by Thomas J. Mueller,
American Institute of Aeronautics and Astronautics, Inc., 2001.
(6)
Liu, H., Ellington, C. P., Kawachi, K., Van den Berg, C., andWillmott, A. P., “A Computational Pfuid Dynamic Study of
Hawkmoth Hovering,” Journal of Experimental Biology, Vol. 201,
No. 4, 1998, pp. 461-477.
(7)
Isogai, K., Shinmoto, Y., “Effects of dynamic stall phenomena on
propulsive efficiency and thrust of a flapping airfoil,” AIAA Paper
97-1926, 1997.
(8)
Ramamurti, R., Sandberg, W., “Simulation of Flow About
Flapping Airfoils Using Finite Element Incompressible Flow
Solver,” AIAA Journal, Vol. 39, No. 2, February 2001.
(9)
岡部能幸，下山幸治，藤井孝藏，“羽ばたき型火星航空機の羽ばたき運動に関する考察”，第 19 回数値流体力学シンポ
ジウム講演論文集，B2-2, 2005(10) Deb, K., Multi-Objective Optimization using Evolutionary
Algorithms, John Wiley & Sons, Ltd., 2001.
(11) Peter J. Kunz, Analysis and Design of Airfoils for Use at
Ultra-Low Reynolds numbers, Fixed and Flapping Wing
Aerodynamics for Micro Air Vehicle Applications (Progress in
Astronautics and Aeronautics Volume 195) Chapter 3, p35-60,
2001(12) Van Leer, B.: Towards the Ultimate Conservation Difference
Scheme. IV. A New Approach to Numerical Convection, Journalof Computational Physics, 1977, Vol.23(3), 276-299(13) Nishida, H. and Nonomura, T.: ADI-SGS Scheme on Ideal
Magnetohydrodynamics, Journal of Computational Physics, 2009,Vol.228, 3182-3188(14) Deb, K., Pratap, A., Agarwal, S. and Meyarivan, T.: A fast andelitist multiobjective genetic algorithm: NSGA-II, volutionary
Computation, IEEE Transactions, 2002, Vol.6(2), pp.182-197.
(15) Goldberg, D.E.: Genetic Algorithms in Search, Optimization andmachine Learning., Addison Wesley, 1989.
(16) Deb, K. and Agrawal, R.B.: Simulated Binary Crossover for
Continuous Search Space, Somplex Systems, 1994, Vol.50(9),pp.115-148.
(17) Deb, K. and Goyal, M.: A Combined Genetic Adaptive Search(geneAS) for Engineering Design, Computer Science andInformatics, 1996, Vol.26(4), pp.30-45.

7

Copyright © 2011 by JSFM