生物の飛行と遊泳における計算バイオミメティクス

4759

昆虫規範型微小羽ばたき翼の流体構造連成設計
石原 大輔

1

緒言

本稿では、昆虫規範型 FWPAV（Flapping Wing Pico Air
Vehicle）開発 [1] の一環として取り組んでいる昆虫規範型
微小羽ばたき翼の流体構造連成設計について述べる。
その大筋は文献 [2] に基づく。
1.1 昆虫規範型 FWPAV
昆虫は飛行の利点により地球全体に生息域を拡大させた。その飛行能力は自然淘汰によって、数億年という永い時をかけて洗練されてきた。そのメカニズムを理解し、模倣することによって、高性能な飛行体の開発が期待される [3]。本稿では昆虫規範型 FWPAV を対象とする。FWPAV は翼、アクチュエータ、トランスミッ
ション、支持フレームなどから構成され、その大きさは5cm 以下である [4]。FWNAV における関心の1つは微小化にある [5]。その潜在的な最小寸法は非常に小さな昆虫の寸法であり、進化の帰結からは1mm 程度以下となる。
1.2 昆虫飛行における流体構造連成
昆虫の飛行能力はピッチング、キャンバリング、翼先端軌跡（8の字など）といった羽ばたき翼の特徴的運動 [6] による。特に、ピッチング運動は昆虫飛行の非定
常流体力学 [7] の運動学的基礎をなす。昆虫翼とその基部は柔軟である [8]。ゆえに大きなストローク角を伴う
羽ばたき運動は周囲の空気流からの流体力による翼の変形を引き起こし、更にそれが周囲の空気流の変化を引き起こす（流体構造連成）。このトピックの最近の進歩に対して、計算力学的方法、特に数値流体構造連成
解析方法 [9, 10] の寄与は顕著である。

則、デザインウィンドウ探索）を示し、それを昆虫規範型羽ばたき翼に適用することで、FWPAV の小型化の可能性を示す。

2

流体構造連成設計に対する計算力学的アプ
ローチ

2.1 昆虫翼の形状簡易化モデル
昆虫翼は多数の複合的な構造要素から構成される [17]。
翼膜を支持する翅脈は、飛行中の荷重に対して十分な強さを保ちつつ、ひねりやキャンバリングといった特徴的な変形を許容して、空力性能を改善する [18]。その基部は多数の骨片、筋肉、腱から構成されており、それらが一体となって働くことで、飛翔筋による駆動力を伝達して翼を羽ばたかせるとともに、翼の運動を巧みに制御するトランスミッションとして働く [19]。これらの構造を全てモデル化することは困難であるため、目的に応じて適切に簡易化が行われる [10, 17]。ここでは、計算コストの低減および模倣の目的のために、図1に示す集中化柔軟性 [2, 9, 14] と翅脈 - 膜相互作用 [15] による形状簡易化モデルを用いる。それは複雑な構造の詳細をバネ、梁、シェルとして隠匿し、弾性変形メカニズムを大局的構成関係により記述する。ねじりバネはピッチング運動を生むねじり柔軟性に対応し [9]、昇降バネは翼先端軌跡を生む昇降
柔軟性に対応し [14]、羽ばたきバネはトランスミッションが有する羽ばたきに関する柔軟性に対応する [2]。さらに、翅脈はキャンバー形成の役割に対応する3つの梁によって表される [15]。このモデル翼は昆虫翼の特徴的運動を模擬することができる。昆虫翼の複雑な構造がバ
Insect Body

1.3 流体構造連成設計
著者らは上述の特徴的な運動が流体構造連成によって受動的に生じ得ることを示してきた [9, 11-16]。この力
学的な翼運動の自動化は昆虫の微小化に寄与していると考えられ、昆虫規範型 FWPAV においても利用されている [4]。これを一般化して、人工物の機能を生み出す
ために積極的に流体構造連成を利用する設計（流体構造連成設計）が考えられる。この設計方法では、流体構造
連成現象の複雑さのため、数値流体構造連成解析を中心とする計算力学的アプローチが有効である。ここでは著者らの流体構造連成設計に対する計算力学的アプ
ローチ（翼モデル、直接数値モデリング、動力学的相似

Elevation Motion
StrokePlane

Wing

CamberFlapping Motion(Down Stroke)

Wing
Chord

Pitching Motion(a) 昆虫羽ばたき翼
Flapping axis
Spring for flappingyTorsional
axis

O

Spring for torsion
Spring for elevation

φ

x
Wing

Φ
Stroke
plane

χ

z

筆者紹介

Elevationaxisθ(b) 形状簡易化モデル
図1

いしはら だいすけ（著者紹介は P.2）

Vol.29, No.2 2024

Leading Edge

（7）

昆虫羽ばたき翼の形状簡易化モデル

f I,Ef MR
exchangeableα

2.2 直接的数値モデリングによる流体構造連成解析
本流体構造連成設計は、薄く柔軟な羽ばたき翼の動的大変形と周囲の流動の連成問題を正確かつ安定に解くことを前提とする。この連成問題の支配方程式は
∂σsji
dvs
ρs i =dt∂x j
ρf

σfi j nfj + σsi j nsj = 0,

vfi = vsi

∂vfi
=0
∂xi

f
f I,C

Re
f Vf

f

f
P

Ca

f

図2

(1)

∂σfji
∂vfi
∂vf
+ ρf (vfj − vmj ) i =
,
∂t
∂x j
∂x j

f I,LssE

Fluid

ネ、梁、シェルで表されているので、それらの実装により、昆虫翼の工学解 [20] を導くことができる。

f fI, E : inertial force due to local accelerationf fI, C : inertial force due to advection
f fV : viscos force
f fP : pressure

Structure

特集

4760

f fI, L : inertial force
f fE : elastic force

流体構造連成の動力学的相似則

2.4 デザインウィンドウ探索
MEMS の構造設計 [23] と同様に、FWNAV の構造設計
では、作製可能性など定式化し難い設計制約を考慮して、設計問題を徐々に修正しながら、最終的な設計解に至るので、唯一の最適解を探索するよりも設計要求を満足する解の集合（デザインウィンドウ）を探索する方が有効な場合がある [24, 25]。

(2a, b)(3a, b)である。ここで式 (1)、(2) はそれぞれ弾性体と非圧縮性
粘性流体の支配方程式、式 (3) は流体構造境界面に課される連成条件である。翼に作用する重力は、その他の項に比べて十分小さいので考えない。これらの有限要素離散式から得られる一体型方程式を、系全体の運動
方程式と流体圧力ポアソン方程式に分割して解くプロ
ジェクション解法 [21, 22] を用いる。流体圧力ポアソン方
程式は、流体の非圧縮性拘束と等価である。ここではSchur 補行列を回避するため、反復法に基づく分割 [22]を用いる。流体構造境界面は ALE 法により直接的に記述され、連成条件は一体型定式化により直接的に解かれる。構造大変形は Total Lagrange 法により定式化される。分割後の方程式は良条件化されているため、メッ
シュ分割による対角スケーリング前処理付き並列 CG 法を用いる。以上の詳細は文献 [16] に譲る。

3

動力学的相似則に関する数値計算例

3.1 問題設定
図1の集中化ねじり柔軟性のみを考慮し、代表量 L、T、V、C をそれぞれ平均翼弦長 cw、平均羽ばたき速さ
Vm = 2r2ΦLwfφ、羽ばたき周波数 fφ の逆数、ねじり柔軟性 Cθ によって与える。ここで r2は翼面積の断面2次モー
メントの無次元半径、Φ はストローク角、Lw は翼長である。このとき、α = 1/(r2rAΦ) となる。但し、rA は2Lw/
cw で定義される翼アスペクト比である。集中化ねじり
柔軟性を表すねじりバネを細長い板バネで与えた実装モデルを、図3に示す。簡単のため、翼形状を長方形とする。長方形の r2は1/30.5であり、多くの双翅目の翼のr2はこの値に近い。周囲の流体領域は十分な大きさの直方体で与え、全ての壁面上に滑り無し条件を課す。
翼前縁に課される強制変位は、正弦波状の時刻歴の羽ばたき角変位によって与える。ここでその1階微分の羽ばたき角速度は、実際の昆虫の観測に基づき、各半ス
トロークの中盤で一定となる台形波によって与える。
振幅は1周期かけて正負のピーク間の大きさが Φ となるまで線形的に増加させる。翼に MITC シェル要素、流体に P1P1要素を用いて離散化する。

2.3 流体構造連成の動力学的相似則
支配方程式 (1)〜(3) の次元解析から求まる流体構造連成系の動力学的相似則は、局所加速度と移流による慣性 力 の 比 α、 移 流 に よ る 慣 性 力 と 粘 性 力 の 比 で あ る
Reynolds 数 Re、 流 体 動 圧 と 構 造 弾 性 力 の 比 で あ る
Cauchy 数 Ca、構造と流体の質量比 MR から構成される：
α = L/(TV), Re = ρf LV/µf , Ca = ρf V 2 L3 C, MR = mw /(ρf L3 )y(4a, b, c, d)

Strokeangle O

こ こ で L、T、V は そ れ ぞ れ 寸 法、 時 間、 速 度 の 代 表
量、μf、C、mw はそれぞれ流体粘性係数、構造コンプ
ライアンス、構造質量を表す。図2に示すように、これらの4つの無次元数は、支配方程式 (1)〜(3) 中の全ての力学項を相互に関連付けるので、幾何学的相似の下で、異なる流体構造連成系を互いに力学的相似にできる。最初に、昆虫データを図1の計算モデルに正確に導入するために用いられ [9]、次に、スケール拡大モデル実験に用いられた [11, 12]。本流体構造連成解析手法とその実装プログラムが支配方程式 (1)〜(3) と整合していることの検証にも用いられた [2, 16]。

x Strokeplane
LeadingedgeΦPlate springPitch angle
Stroke angular z
displacement
Wing plate
図3

集中化ねじり柔軟性翼モデルの実装

3.2 2つの異なる流体構造連成系の動力学的相似
双翅目の昆虫とそのスケール拡大モデルを考える。
前者に対して後者の翼長を約20倍、羽ばたき周波数を約1/200倍、流体を空気からシリコン油に変更する。幾何 学 的 相 似 の 無 次 元 パ ラ メ ー タ は r2 = 1/30.5、rA =
（8）

計算工学

生物の飛行と遊泳における計算バイオミメティクス

7.267、Φ = 107.5°であり、動力学的相似の無次元パラメータは α = 0.127、Re = 824、Ca = 0.244、MR = 0.111である。2つの異なる流体構造連成系でピッチ角と無次元化された揚力の時刻歴（図4）は非常によく一致した。
ゆえに数値流体構造連成解法とその実装プログラム
は、支配方程式 (1)〜(3) と整合している。
Lift coefficient

4

2

0.0

0.5

1.0
Time [cycle]

Insect

図4

1.5

2.0

Dynamically scaled

揚力係数の時刻歴

3.3 受動的翼運動における羽ばたき周波数の効果
動力学的相似則に基づき、羽ばたき翼の性能を系統的に評価する。ガガンボの無次元パラメータ α = 0.073、
Re = 254、Ca = 0.19、MR = 16、Φ = 123°、rA = 10.7に対して、fφ を変えることで、無次元周波数 Fr = fn/f（
φ fn は翼の 固 有 周 波 数）を3.32〜5.80と 変 え な が ら、Φ を80〜140°と調整することで、Re、Ca、MR を一定に保ちつつ、α のみ0.112〜0.0637と変わるようにした。この設定は、翼慣性力だけをその他の力学項に対して相対的に変化させるものである。この結果、図5に示すように、ピッチ角がストロークの切り返しの直後で顕著に変化する一方で、各半ストローク中盤ではあまり変化しないように調整された。図6にストロークの切り返

Pitch [deg]

90
45
0
-45

6

6.25

-90

6.5

7

Time [cycle]

Fr = 3.3

Fr = 4.1

図5

Fr = 5.1

Fr = 5.8

ピッチ角の時刻歴(a) Φ = 123°、Fr = 5.1、
α = 0.073 (cranefly)
図6

6.75

(b) Φ = 80°、Fr = 3.3、
α = 0.11

6.65 cycle における流れ場。コンターは z 方向の渦度、ベクトルは流速を表す。

Vol.29, No.2 2024

2.5次元構造を有する微小羽ばたき翼の流体構
造連成設計

図7に2.5次元構造を有する微小羽ばたき翼の流体構
造連成設計の概念を示す。流体構造連成による受動的
ピッチング運動に加えて、微小な基部励振が共振によって十分大きな羽ばたき角に増幅される。前章で示したように、羽ばたき周波数は受動的運動で重要なパラメータなので、それを設計パラメータとする。この翼の最も顕著な特徴は2.5次元構造であり、MEMS マイ
クロマシニングによって作製可能である [2, 20]。

1
0

し直後にピッチ角がピークをとる時点付近での流れ場を示す。同図に示すように、いずれの場合にも昆虫羽
ばたき飛行の非定常流体力学理論で重要な前縁渦が明確に捉えられている。また Fr が小さい方が下向きに強い流れが生じており、その時点のピッチ角と揚力がよ
り大きいことと対応している。

4

3

4761

（9）

y
Base excitationof micro actuator

PlatespringAmplification
using resonanceFlapping motion
Aerodynamicforce
FSI-cause ofpitching motion

Gravity
Lift forcex

Leading edge
beam
z

Wing
membrane

Pitching motion

図7 2.5次元構造を有する微小羽ばたき翼の流体構造連成設計

翼膜の材質をポリイミド、翼長2.5mm、翼弦長0.8mm、膜厚 tw = 1.6μm とした。これらの寸法は小さなハエと同程 度 で あ る。 翼 前 縁 の 材 質 を シ リ コ ン 単 結 晶、 幅
0.1mm、厚さ0.05mm とした。但し、このような複合材の微小翼の作製は困難であるため [26]、現在は全てポリイミドで作製している [5, 20, 27, 28]。板バネは翼膜と同じ層に形成されるとし、その長さを0.05mm とする。周囲の流体は空気とする。翼基部の励振は、利用可能な圧電アクチュエータから、面外方向振幅0.08mm、羽ばたき周波数 fφ = 100〜1000Hz と仮定する。揚力 FL が小さなハエの自重と同程度以上である場合を満足解とする。図8に fφ = 428Hz のときの翼基部と翼先端の面外方向変位の時刻歴を示す。共振によって後者は前者の19倍に増幅された。このとき、準定常状態において、ス
トローク角とピッチ角はそれぞれ実際の昆虫と同程度の73°と62°に達した。このとき両翼による揚力（鉛直方向の流体力の総和）の平均は実際の昆虫の自重を超える。このときの翼運動と周囲の流れを図9に示す。同図に示すように、ピッチング運動が適切に生じていること、前縁渦が各半ストローク中盤で安定的に付着していることが確認できる。揚力の周波数応答関数から、
デザインウィンドウは fφ = 410〜460Hz であり、満足解の存在が示された。

特集

4762

■参考文献

Tip displacement in the
direction of x-axis [mm]

2

[1]

1
0

0

5

10

[2]
[3]

15

-1

[4]

-2
図8

Time [cycle][5]

翼基部と翼先端の面外変位の時刻歴

[6]
[7]
[8]
[9]

(a) 8.48 cycle

[10]
[11]

(b) 8.72 cycle

[12]
[13]
[14]
[15]

(c) 8.96 cycle
図9

[16]

(d) 9.20 cycle[17]共振駆動による翼運動と周囲の流れ（https://www.mdpi.com/2311-5521/7/1/26の動画参照）

5

結言

流体構造連成による翼運動の自動化は、人工的な羽ばたき飛行体の電気機械的複雑さの低減に寄与する。
この流体構造連成設計の複雑さを取り扱うために、本稿ではその計算力学的アプローチについて述べ、昆虫
規範型微小羽ばたき翼への適用を通じて、FWNAV の
小型化の可能性を示した。

謝辞

[18]
[19]
[20]
[21]
[22]
[23]
[24]

東レ株式会社の支援および JSPS 科研費の助成に謝意
を表する。

[25]
[26]
[27]
[28]

（10）

昆虫の自由飛行と知能を包括したマルチフィジクス連成
解析による真の計算 Bio-Mimetics、基盤研究 (A)、課題番
号23H00475.
Ishihara, D. Fluids 2022, 7, 26.
Liu, H.; Ravi, S.; Kolomenskiy, D.; Tanaka, H. Philos. Trans. R.
Soc. 2016, B371, 20150390.
Wood, R.J.; Finio, B.; Karpelson, M.; Ma, K.; Pérez-Arancibia,
N.O.; Sreetharan, P.S.; Tanaka, H.; Whitney, J.P. Int. J. Robot.
Res. 2012, 31 (11), pp. 1292-1302.
Ishihara, D.; Murakami, S.; Ohira, N.; Ueo, J.; Takagi, M.;
Urakawa, K.; Horie, T. Proceedings of the 15th Annual IEEE
International Conference on Nano/Micro Engineered andMolecular Systems 2020, pp. 27-30.
Ellington, C.P. Philos. Trans. R. Soc. Lond. 1984, B305, 41-78.
Dickinson, M.H.; Lehmann, F.-O.; Sane, P.S. Science 1999,
284, 1954-1960.
Ennos, A.R. J. Exp. Biol. 1988, 140, 137-160.
Ishihara, D.; Horie, T.; Denda, M. J. Exp. Biol. 2009, 212,
1-10.
Nakata, T.; Liu, H. J. Comput. Phys. 2012, 231, 1822-1847.
Ishihara, D.; Yamashita, Y.; Horie, T.; Yoshida, S.; Niho, T. J.
Exp. Biol. 2009, 212, 3882-3891.
Ishihara, D.; Horie, T.; Niho, T. Bioinspiration Biomim. 2014,
9, 046009.
Ishihara, D.; Horie, T. Bioinspiration Biomim. 2017, 12,016008.
Ishihara, D. Phys. Rev. E 2019, 98, 032411.
Ishihara, D.; Onishi, M.; Sugikawa, K. Biomimetics 2023, 8,571.
Ishihara, D.; Onishi, M. Int. J. Numer. Methods Fluids 2024,
96 (4), 435-481.
Wootton, R.J.; Herbert, R.C.; Young, P.G.; Evans, K.E. Philos.
Trans. R. Soc. Lond. 2003, B358, 1577-1587.
Du, G.; Sun, M.J. Exp. Biol. 2010, 213, 2273-2283.
Beatus, T.; Cohen, I. Phys. Rev. E 2015, 92, 022712.
Shankar, V.; Shirakawa, N.; Ishihara, D. Biomimetics 2024, 9,133.
Ishihara, D.; Yoshimura, S. Int. J. Numer. Methods Eng. 2005,
64, 167-203.
Ishihara, D.; Horie, T. Comput. Modeling Eng. Sci. 2014, 101,
421-440.
Ishihara, D.; Jeong, M.J.; Yoshimura, S.; Yagawa, G. Comput.
Struct. 2002, 80, 2469-2481.
Rashmikant, Ishihara, D. Int. J. Mech. Mater. Des. 2023, 19,407-429.
Rashmikant; Onishi, M.; Sugikawa, K.; Ishihara, D. Eng. Appl.
Comput. Fluid Mech. 2023, 17 (1), 2271053.
Murakami, S; Ishihara, D.; Araki, M.; Ohira, N.; Ito, T.; Horie,
T. Micro Nano Lett. 2017, 12 (11), 913-915.
Rashmikant; Ishihara, D.; Suetsugu, R.; Ramegowda, P. C.
Eng. Res. Express 2021, 3, 045006.
Rashmikant; Suetsugu, R.; Onishi, M.; Ishihara, D.
Proceedings of 2023 IEEE 36th International Conference on
Micro Electro Mechanical Systems 2023, 534-537.

計算工学