数理解析研究所講究録
第 1808 巻 2012 年 127-132

127

羽ばたき飛翔する蝶のつくる流れと飛行安定性
京大工広大理阪府大生命環境科学

横山直人 (Naoto Yokoyama), 泉田啓 (Kei Senda)
Faculty of Engineering, Kyoto Univ.
飯間信 (Makoto Iima)
Faculty of Science, Hiroshima Univ.
平井規央 (Norio Hirai)

Grad. School of Life & Environmental Sciences, Osaka Pref. Univ.

1

はじめに

人類が製作する航空機の多くは，固定翼あるいは回転翼を持ち，非定常に動作する部分はわずかである．一方で，飛翔する動物の多くは，翅を非定常に羽ばたかせることによって飛翔する．その中でも，ヒラヒラと飛翔する蝶の飛翔は，その美しさから興味深いものである．
昆虫の羽ばたき飛翔は，羽ばたき運動によって生成する渦構造によって特徴づけられる．渦構造と昆虫に作用する力の関係は，実験的・数値的・理論的研究が現在精力的になされている [1,3,5,6].
本研究では，アサギマダラを，胸部，腹部，両翅の 4 リンクでモデル化し，その飛翔を数値シミュレーションした．風洞実験から得られた自由飛翔する生体のアサギマダラの胸部軌道と関節角の運動をモデルの蝶に与えるシミュレーションでは，アサギマダラの飛翔が生成する流れ場を可視化し，その渦構造と蝶に作用する力トルクの関係を明らかにする．また，関節角は風洞実験から得た運動を用いるものの，胸部軌道は運動方程式を解くシミュレーションでは，飛翔の安定性を数値的に調べる．

2

数値計算法
流れ場は，非圧縮粘性流体に対する Navier-Stokes 方程式$\frac{\partial u}{\partial t}+(u\cdot\nabla)u=-\frac{\nabla p}{\rho_{air}}+\nu\nabla^{2}u+f_{IB}$$\nabla\cdot u=0$

(la)

,(lb)

,によって記述される．ここで， は，流れと構造としての蝶の間の相互作用を埋め込み境界法 [4]での値 $\rho_{air}=1.184kg/m^{3}$ を用いる．動粘性係によって求めた体積力を表す．空気の密度は数には数値振動を低減するために，実際より 2 倍大きい動粘性係数の値 $v=2v_{air}$ を用いる．ここ$v_{air}=1.54\cross 10^{-5}m^{2}/s$ は，での空気の動粘性係数である．
で，空気は，流入境界において一様速度 $u=u_{0}e_{x}(u_{0}=1.6m/s)$ で流入する．流出境界条件にはSommerfeld 条件を用いる．また，流入及び流出境界における圧力の境界条件には，Neumann 条件を用いる．主流方向 の計算領域の長さは $5\cross 10^{-1}m$ とした．主流と垂直な方向 (鉛直方向$fi_{B}$$25^{o}C$$25^{o}C$$e_{y}$$e_{x}$

および水平方向$e_{z}$) の境界条件には周期境界条件を用い，その長さは$2.5\cross 10^{-1}m$とした．

計算領域は，主流方向に 4096 点，主流と垂直な方向に 5122 点の等間隔格子によって離散化される．主流方向の微分は 4 次精度中心差分によって得られ，主流と垂直な方向の微分は Fourier 級数によって評価される．

蝶はアサギマダラをモデル化し，楕円体の胸部，同じく楕円体の腹部，左右の両翅の 4 リンク系で近似される．実際のアサギマダラは前翅と後翅を持つがこれらを 1 枚の翅で近似する．蝶の状

128
態は，一般化座標$\theta=(x_{t}, \phi_{t}, \phi_{a}, \phi_{1}, \phi_{r})\in \mathbb{R}^{15}$とその時間微分$\dot{\theta}$

で与えられる．ここで， は$x_{t}$は，胸部の重心座標， は慣性座標系から 1-2-3 オイラー角で変換される胸部の姿勢，それぞれ，胸部の重心に張った座標から 1-3-2 オイラー角で変換される腹部，左翅，右翅の関節角で$\phi_{a},$$\phi_{t}$$\phi_{1},$$\phi_{r}$ある．

本研究では，左右対称な羽ばたき運動を行うと仮定する．その結果，流れ場も矢状面に対して,ほぼ対象である．このとき，一般化座標は，$x_{t}=(x_{t}, y_{t}, 0),$の 7 自由度のみをもつ．ここで， は胸部重心の主流方$\phi_{t}=(0, \theta_{t}, 0),$$\phi_{1}=(\beta+\pi, -\eta, -\theta),$$\phi_{r}=(-\beta, -\eta, \theta)$$\phi_{a}=(0,0, -\theta_{a})$$x_{t}$

向の位置， は鉛直方向位置， は胸部姿勢， は腹部関節角， は概フラッピング角， は概りードラグ角， は概フェザリング角である．
一般化座標の支配方程式は，Euler-Langrange 方程式$\theta_{a}$$\theta_{t}$$y_{t}$$\beta$$\eta$$\theta$(2)$M \ddot{\theta}+\dot{M}\dot{\theta}-\frac{1}{2}\frac{\partial}{\partial\theta}(\dot{\theta}M\dot{\theta})+\frac{\partial V}{\partial\theta}=\tau_{air}+\tau_{cont}$で与えられる．ここで，は質量行列， は重力によるポテンシャルである．また，は埋め込み境界法によって得られた流れと蝶の相互作用による力 $fiB$ を，一般化座標 に対応す
は胸部に作用する制御トルクや関節が発生する制御トルる一般化力に変換したものであり，$M\in \mathbb{R}^{15\cross 15}$$V$$\tau_{air}$$\theta$$\tau_{cont}$

クである．

本研究では，2 種類の数値計算を行う．1 つの計算では，生体の蝶を風洞内で自由飛翔させた実
験映像から蝶の胸部飛翔軌道や関節角を読み，一般化座標 を得る．このとき，一般化座標 の時$\theta$系列を得ることで， と$\dot{\theta}$によって$\tau_{air}$$\ddot{\theta}$$\theta$も得られ，式 (2) の左辺が得られる．また，数値計算内で埋め込み境界法が得られるので: 結果的に制御トルク$\tau_{c}$$nt$。

が得られる．この数値計算を規程運動シミュレーションと呼ぶことにする．

もう 1 つの数値計算では，生体の蝶の自由飛翔映像から蝶の関節角を読み，これを用いる．すなと関節座標に形式的に分解すると，
わち，一般化座標 を胸部座標$\theta$$\theta_{2}$$\theta_{1}=(x_{t}, \theta_{t})$$\theta_{2}=(\phi_{a}, \phi_{1}, \phi_{r})$に風洞実験で得られた関節角の値を与える．質量行列や一般化力も同様に分解し，胸部座標に関する部分を添字 1 で表し，関節座標に関する部分を添字 2 で表すことにする．このとき式 (2) は
形式的に

(3)

,

$(\begin{array}{ll}M_{11} M_{12}M_{21} M_{22}\end{array})(\begin{array}{l}.\theta_{1}.\cdot\theta_{2}\end{array})+(\begin{array}{l}h_{1}h_{2}\end{array})=(\begin{array}{l}\tau_{cont1}\tau_{cont2}\end{array})$と書き換えることができる．ここで，$h=(h_{1}, h_{2})= \dot{M}\dot{\theta}-\frac{1}{2}\frac{\partial}{\partial\theta}(\dot{\theta}M\dot{\theta})+\frac{\partial V}{\partial\theta}-\tau_{air}.$である．一般的に，蝶は関節の運動を制御することは可能であるが，体幹である胸部の位置及び姿$\tau_{cont1}=0$ と仮定する．このとき，勢を直接制御することは難しい．したがって，$\ddot{\theta}_{1}=-M_{11}^{-1}(M_{12}\ddot{\theta}_{2}+h_{1})$を得る．同時に，$\tau_{cont2}=h_{2}-M_{21}M_{11}^{-1}(M_{12}\ddot{\theta}_{2}+h_{1})$.

(4)も得られる．式 (4) を 2 次 Adams-Bashforth

法によって時間発展させるこのシミュレーションを，自由飛翔シミュレーションと呼ぶことにする．

3

数値計算結果3.1

規程運動シミュレーション

生体の蝶を風洞内で自由飛翔させた実験映像から得られた一般化座標$\theta$

によって飛翔運動を再現する規程運動シミュレーションは，系が概ね周期的な状態での結果を示す．ここで，系が周期的

129Time $[s]$Time $[s]$図 1:(左) 蝶に作用する揚力，抗力，胸部重心まわりのピッチングモーメント．(右) 実験で得られた飛翔軌道を再現するための補償揚力，補償抗力，補償トルク．になった後で，翅の打ち下ろしを開始する時刻，すなわち， が最大となる時刻の 1 つを時刻の原点 $t=0$ に選んだ．規程運動シミュレーションにおける 1 周期 $(T=0.112s)$ 間の，流体による揚に力，抗力，胸部重心まわりのピッチングモーメントを図 1 に示す．翅の打ち下ろしの間$\beta$$t\lessapprox T/2$揚力が得られ，打ち上げの間 $t>T\sim/2$ に負の抗力，すなわち，推力が得られている．また，打ち下ろしの間に負 (頭を下げる向き) のピッチングモーメントが作用し，打ち上げの間に正 (頭を上げる向き) のピッチングモーメントが作用している．これは，翅に作用する力の平均的な位置 (いわゆる空力中心) が，胸部重心よりも下流側に存在するためである．
また，実験で得られた胸部軌道$\theta_{1}$を再現するために，胸部に作用する補償力$\tau_{c}$$ntl\neq 0$。

が必要である．これらを図 1 右に示す．これらが小さいほど自由飛翔実験が再現されたことに対応するが，残念ながら，補償揚力や補償抗力は流体による揚力や抗力と同程度の大きさ，ピツチングモー
メントに関しては，数倍の大きさとなっている．補償力では，基準振動数に対して 3 倍あるいは 4倍振動が顕著である．これは，胸部加速度 において 3 倍あるいは 4 倍振動が顕著であることに$\ddot{\theta}_{1}$由来し，時間・空間的に高精度な実験的観測が必要であることを示唆している．
流れ場の構造を $Q$ 値，すなわち，速度勾配テンソルの 2 次不変量 $Q=\nabla^{2}p/(2\rho_{air})>0$ を用いて可視化する．ここで $t=T/4$ における $Q=1.5\cross 10^{4}s^{-2}$ の等値面を図 2 に示す．翼端付近に，数値的な振動が見えるが，これらが生成する力は無視できるほど小さく，移流もされないので，これらの影響はないものと考える．また， 値により大きな値を用いるとこの振動は現れない．流れ場は非常に複雑であるが，その左右対称性は高いことがわかる．
複雑な渦構造の中にいくつかの強い渦構造が見られる．打ち下ろし中の前縁にはシート状の前$Q$

縁渦 (LEVd) が形成されている．また，両翼端から細い渦が螺旋状に絡みあった翼端渦 (WTVd)が伸びている．この細い渦は右翅から出発したものは主流方向の渦度を持ち，左翅から出発したものは主流と逆行する方向の渦度を持つ．このため，左右翅から出発した WTVd に挟まれる領域では，鉛直下向きの流れが誘起され，その反力として蝶に揚力が作用する．また，粘性によって細い渦は後流で結合し，1 本の渦を構成する．さらに，打ち上げの際に生成される翼端渦 (WTVu)

は，右翅から出発したものは鉛直上向き，左翅から出発したものは鉛直下向きの渦度を持つ．したがって，WTVu の間では，主流方向の流れが加速され，その反力が蝶に作用する推進力である．
WTVd の上流に，打ち下ろしから打ち上げへの間に生成される後縁渦 (TEV) がリング状の見ら付近の揚力れる．TEV は鉛直下向きの流れを誘起し，蝶に揚力を生成し，これが，図 1 のの平坦な領域の起源である．また，後流において WTVu と TEV の間に，非常に乱れた構造が存在する．これは，打ち下ろしの際にできた LEVd の下に形成される 2 次流れ (SFV) が後流に流れた$t\approx T/2$

130

$Q=1.5\cross 10^{4}s^{-2}$ . (上) 上方から見た図．(下) 左
図 2: 等 $Q$ 値面によって可視化された流れ場．

から見た図．

ものである．このように，流れ場の構造は，Brodsky[2] が提案した “chain of vortices” よりはるかに多彩である．

流れ場の構造と蝶に作用する力の関係をより詳細に見るために，$t=T/4$ における，蝶に作用する力の分布と値による渦構造を図 3 に示す．図 2 で用いた $Q=1.5\cross 10^{4}s^{-2}$ よりも大きな$Q=1\cross 10^{5}s^{-2}$を用いたため，打ち下ろし中の LEVd が前縁付近のみで見えるようになり，その下の構造が見えている．LEVd の下の SFV の中にも比較的強い渦管 ( $VT$ ) が存在することがわか$Q$る．翅に作用する力は，LEVd の下，特に $VT$ が存在する部分で大きいことが見て取れる．$VT$ の分布の蝶に作用する揚力や抗力に対する寄与は大きくないが，ピッチングモーメントに対する影$VT$ のよう響は大きい．したがって，蝶の姿勢変化を伴う飛翔をシミュレーションするためには，な微細な構造を解像することが必要である．

3.2

自由飛翔シミュレーション

生体の蝶を風洞内で自由飛翔させた実験映像から得られた関節角
胸部軌道$\theta_{2}$を与え，式 (4) に基づいてを時間発展させる自由飛翔シミュレーションは，規程運動シミュレーションで十分周期的になった後，概フラッピング角が最大となる時刻以降 $ntl=0$ とするシミュレーションによってなされる．自由飛翔シュミレーションでの 1 周期間の蝶に作用する力の時間変化と胸部軌道を，図 4 に示す．揚力は十分に発生しているため，1 周期後の鉛直方向の位置は，シミュレーション開始時とほぼ同じである．一方で，推力の発生は充分でなく，流入する空気に対して後退していることがわかる．ここで，流入する空気とともに移動する座標系から見れば，1. $lm/s$ 程度で前進していることに注意しておく．しかしながら，姿勢は周期的でなく，打ち下ろし時の蝶の頭が十分$\theta_{1}$$\tau_{c}$。

131$Q=1\cross 10^{5}s^{-2}$ . 翅に作用する
図 3: 蝶に作用する力と等 $Q$ 値面によって可視化された渦構造．
力の大きさは色の濃さで表現される．

41’$y_{t}x_{\theta_{t}^{t}}$$\equiv$$\prime^{\prime’}\prime\cdots\cdots\prime\prime^{\prime’}../’/\cdot\cdot-\prime^{\prime’}$$N\frac{\overline{}\vee oE\cross}{o\fcircle\geq},$$23$$\frac{}{o}\frac{o}{o,\infty}$

1$\prime\prime^{\vee^{\nearrow}}\prime.\cdot\prime^{\prime’}\prime.\prime.\cdot\prime.\cdot\cdot\prime’.\cdot/^{\prime’}$

0.

$8$

0.

$6$

0.

$4$

$\overline{\frac{\underline{\varpi\varpi}}{\triangleleft\frac {}{}c\infty\omega}}$$\mathring{\overline{L}\equiv\vee co}$$\vdash\sum_{arrow,\circ} 0 arrow--.--.-.-.\vee\prime\vee^{-}$$\underline{\frac{o}{\varpi 0}}$

$1L\simeq\subset 0\underline{o}\omega -1 \vdash\Sigma 0$
$-2$
$0$

0.028

0.056

0.084

0.2
0.112

Time $[s]$Time $[s]$図 4: 自由飛翔シミュレーションにおける (左) 蝶に作用する揚力，抗力，胸部重心まわりのピッチ

ングモーメントと (右) 胸部重心位置および胸部姿勢．

下がらず，1 周期後の胸部ピッチング角が lrad 程度まで大きくなり，その後仰向けになって墜落する．これは，打ち下ろし中の負のピッチングモーメントの絶対値が，腹部に作用する重力が胸部重心まわりに作る正のトルク5Nm より小さいためである．
$3\cross 10^{-}$

打ち下ろし中の負のピッチングモーメントの小ささを解消するために，羽ばたきの平均概りードラグ角をわずかに変えた自由飛翔シミュレーションを行った．観測された平均概リードラグ角はrad であったが，平均概リードラグ角をに変えた．その他の羽ばたき運動$\overline{\eta}\approx-9.1\cross 10^{-2}$$7\overline{\eta}/4$のパラメータは，前の自由飛翔シミュレーションと同じものを用いた．このときの 3 周期間の蝶に作用する力の時間変化と胸部軌道を図 5 に示す．平均概リードラグ角の変化はわずかなため，揚力や抗力や抗力，ピッチングモーメントの変化は大きくない．しかし，胸部の姿勢は不安定ではある
2rad だ
ものの，3 周期程度は適切な範囲内に入っている．これは，平均概リードラグ角をけ小さく，すなわち，翅をわずかに後方に引くことによって，打ち下ろし時の負のピツチングモー$6.8\cross 10^{-}$

メントをわずかに大きくしたことによる．このような羽ばたき運動のわずかな変化が飛翔軌道の大きな変化を作り出している．

132

$0$

0.056 0.112 0.168 0224

0.280.336$0$Tlme $[s]$0.056

0.112 0.1680224

028

0.336

Time $[s]$図 5: 概リードラグ角を変化させた自由飛翔シミュレーションにおける (左) 蝶に作用する揚カ，抗力，胸部重心まわりのピッチングモーメントと (右) 胸部重心位置および胸部姿勢．

4

まとめ
本研究では，生体の蝶の羽ばたき運動に基づいて，胸部軌道と関節角のすべての羽ばたき運動を与える規程運動シミュレーションと，関節角のみを与え胸部軌道は Euler-Lagrange 方程式に基づいて数値的に解く自由飛翔シミュレーションを行った．規程運動シミュレーションでは，羽ばたき運動が生成する渦構造を可視化し，渦構造と蝶に作用する力の関係を明らかにした．ここで，これまで考えられてきたものより，流れ場は複雑であることがわかった．その中でも，強い渦構造が存在し，揚力や推力の発生源であることがわかった．一方で，自由飛翔シミュレーションでは，わずかな羽ばたき運動の差が飛翔軌道の大きな変化をもたらすことを示した．

本研究の数値計算は京都大学基礎物理学研究所の計算設備にて行った．

参考文献[1] H. Aono, F. Liang, and H. Liu. Near- and far-field aerodynamics in insect hovering flight:$Exp$ . Biol., 211:239-257, 2008.
an integrated computational study.
[2] A. K. Brodsky. The evolution of insect flight. Oxford Univ. $Pr.$ , 1994.
[3] M. H. Dickinson and K. G. G\"otz. The wake dynamics and flight forces of the fruit fly$Exp$ . Biol., 199:2085-2104, 1996.
drosophila melanogaster.
[4] R. Mittal and G. Iaccarino. Immersed boundary methods. Annu. Rev. Fluid Mech., 37:239261, 2005.
[5] R. Ramamurti and W. C. Sandberg.
three-dimensional computational study of the$Exp$ . Biol., 205:1507-1518, 2002.
aerodynamic mechanisms of insect flight.
[6] X. X. Wang and Z. N. Wu. Stroke-averaged lift forces due to vortex rings and their mutualinteractions for a flapping flight model. J. Fluid Mech., 654:453-472, 2010.
$J.$

$J.$

$A$

$J.$