JP2005263073A - 浮上移動システム及び浮上移動装置 - Google Patents

Abstract

【課題】 流体との相互作用による複雑な羽の変形や、地面などの近隣障害物および周囲流速の影響を考慮した制御が可能で、狭所での作業に適した羽ばたき浮上移動システムおよび浮上移動装置を提供する。
【解決手段】 制御装置は、流体・構造連成解析を用いたfree-flightシミュレーションを用いて経路予測する。それにより、周囲障害物および自己の加速度からの地面効果その他の流体の影響を正確に勘案し、浮上移動装置の羽部の羽ばたき駆動手法を決定する。その後、制御装置は、その駆動手法をアクチュエータに行なわせる処理を実行する。
【選択図】 図２６

Description

本発明は、移動システムに関するものであり、特に、羽を羽ばたかせることによって所定の力を生み出すことで浮上を可能にする浮上移動システムおよび浮上移動装置に関するものである。

従来の航空機に対して、機動性において優れている羽ばたき飛行の工学的実現を目指した研究が、近年盛んになっている。

Georgia Tech Research Institute のRobert Michelson らは、“Beyond Biologically-Inspired Insect Flight,”（von Karman Institute for Fluid Dynamics RTO/AVT Lecture Series on Low Reynolds Number Aerodynamics on Aircraft Including Applications in Emerging UAV Technology, Brussels Belgium, 24-28 November 2003）の中で、羽ばたき飛行を用いた惑星探査を提案している。この他に、災害救助、高所作業等、人間の手の届かない領域での作業への浮上移動システムの応用が期待されている。
"Beyond Biologically-Inspired Insect Flight," （von Karman Institute for Fluid Dynamics RTO/AVT Lecture Series on Low Reynolds Number Aerodynamics on Aircraft Including Applications in Emerging UAV Technology, Brussels Belgium, 24-28 November 2003）

しかしながら、接地移動とちがい、浮上移動については、浮力を用いるものを除き、常に重力に対抗する力を発生させている必要があり、それも、特に狭所での作業では、地面効果によりこの力の変動も大きいため、これらを正確に予測、算出する必要があるが、従来の手法ではこれらを正確に求めることは不可能であった。以下、具体的に例を挙げて説明する。

上述のRobertらの論文では、羽ばたきに飛行における流体力の発生は、羽を剛体と仮定した第一近似的なものであり、かつ、障害物などの影響は考慮されていない。

ところが、実際には羽は変形せず飛行するというのは不可能であり、また、狭所では、地面効果により数倍程度の浮上力変動が起きることが考えられる。また、上記狭所内に元来存在する流れ（以後、周囲流速と称する）による浮上力変動も考えられる。これらの浮上力変動により羽の変形も当然異なったものとなり、両者の相乗作用により、飛行形態は当初の予測とは全く異なったものになってしまう。

以上より、従来提案されてきた手法では、複雑な羽の変形や、地面などの近隣障害物および周囲流速の影響を考慮した制御は不可能であり、最も活用が求められる狭所での作業に適さないと言う問題があった。

本発明は、上記従来の問題点を解決するためになされたものであり、狭所でも作業可能な浮上移動システムおよび浮上移動装置を提供することである。

本発明における浮上移動システムは、羽部とこの羽部を駆動する駆動部、およびこの駆動部を制御する制御部とを有し、この駆動部を用いて、当該羽部を能動的に運動させることで得られる周囲流体からの反力によって移動を行なう停空可能な浮上移動装置を用いた浮上移動システムであり、上記制御部において流体・構造連成解析を用いたfree-flightシミュレーションを行ない、この結果を用いて、制御を行なうことを特徴とする。

この構成により、羽の変形、およびその相乗作用による浮上力変化に正確に対応した浮上移動システムが実現される。

また、本発明における浮上移動システムは、当該浮上移動装置の構造モデルを予め有することを特徴とする。

これにより、制御時における計算を削減することができる。なお、ここで言う構造モデルとは、羽や胴体、およびその周囲の、羽や胴体と連結されている、流出しない部位を指す。

また、本発明における浮上移動システムは、上記駆動部の出力を、上記制御部からの入力信号に関連づけて予め有することを特徴とする。

また、本発明における浮上移動システムは、上記関連づけについて、上記駆動部に対する負荷の影響が考慮されていても良い。

一般に、駆動部のシミュレーションは、電磁場構造連成解析や燃焼解析などの複雑なシミュレーションを必要とするが、これをブラックボックスとして、予め数値的関連づけとして有することにより、計算を簡略化することができる。

また、本発明における浮上移動システムは、上記駆動部の出力に基づき、前記構造モデルの節点における境界条件を与えることを特徴とする。

これにより、制御部からの制御信号により駆動部より出力された力または変位が、構造の駆動として当該浮上装置を運動させるモデルが簡便に実現できる。

また、本発明における浮上移動システムは、上記駆動部の変位出力に基づき、前記構造モデルの節点における変位を境界条件として与えることを特徴とする。

駆動部における駆動源が超音波モータやステッピングモータなどの変位出力であるアクチュエータである場合は、上記の手法により最も簡便に駆動部の挙動が構造モデルに反映される。

また、本発明における浮上移動システムは、上記駆動部の力出力に基づき、前記構造モデルの節点における外力を境界条件として与えることを特徴とする。

駆動部における駆動源がボイスコイルモータや内燃機関などの力出力であるアクチュエータである場合は、上記の手法により最も簡便に駆動部の挙動が構造モデルに反映される。

ひるがえって、本発明における浮上移動システムは、上記free-flightシミュレーション内でのモデルに置いて、羽部と胴体部の接合部は、モデル上は節点を共有せず、計算上、上記接合すべき節点間で、節点外力もしくは節点変位を相互に与えることを特徴とする。

これにより、計算上の取り扱いが非常に厄介であるが空力的には殆ど影響のない、胴体と羽の接合部分をモデルに含むことなく、かつ、力学的に矛盾無く取り扱うことができる。

また、本発明における浮上移動システムは、胴体のモデルを流体・構造連成解析中に含まず、構造についての計算のみ行なうことを特徴とする。

胴体の空力的影響が、構造に基づく力（慣性力、粘性力、弾性力等）に比べて非常に小さい場合は、胴体のモデルを流体・構造連成解析から独立させ、構造についてのみの計算を行なうことで計算量を減らすことができる。また逆に、計算機能力を、胴体の挙動を綿密に解析することに振り向けることができるため、影響の大きい、構造に基づく力の計算を正確に行なうことができるため、より正確にfree-flight解析を行なうことが可能になる。

また、本発明における浮上移動システムは、上記流体・構造連成計算に置いて、羽の位置が胴体に対する相対位置で与えられることを特徴とする。

これにより、最も簡便に駆動装置の変位を羽の運動に反映させることができる。

また、本発明における浮上移動システムは、上記流体・構造連成計算に置いて、胴体の運動が、羽の節点における節点力によって定められることを特徴とする。

これにより、最も簡便に羽の力を胴体に作用させることができる。

また、本発明における浮上移動システムは、上記流体・構造連成計算に置いて、羽の節点における外力を、当該羽の駆動装置から発生する力に基づき与えることを特徴とする。

これにより、最も簡便に駆動装置の発生力を羽の運動に反映させることができる。

また、本発明における浮上移動システムは、上記流体・構造連成計算に置いて、当該計算における流体モデル全体および胴体モデルに対して、オフセットを与えて（相殺速度を加えて）行なうことを特徴とする。

また、本発明における浮上移動システムは、上記流体・構造連成計算に置いて、胴体の運動ベクトルのマイナス１倍の運動ベクトルを、当該計算における流体モデル全体および胴体モデルにオフセットとして与えることを特徴とする。（前述の相殺速度として胴体の速度ベクトルの−１倍のベクトルが流体モデル全体の対地速度の速度ベクトルに加えられる。）
これにより、胴体の解析領域に対する移動が最小限にとどめられるため、計算の効率を高めることができる。すなわち、ある誤差の範囲内では、解析領域を最も小さくとることができるため、最も計算量を小さくできる。このオフセットは、胴体の運動ベクトルのマイナス１倍であるときに、この効果が最大になる。

また、本発明における浮上移動システムは、速度実測手段を有し、当該実測手段によって得られた速度と、上記流体・構造連成計算によって算出された速度との差を、周囲流速とすることを特徴とする。

これにより、最も簡便に周囲の気流状態を得ることができる。

ひるがえって、本発明における浮上移動装置は、羽部とこの羽部を駆動する駆動部、およびこの駆動部を制御する制御部とを有し、この駆動部を用いて、当該羽部を能動的に運動させることで得られる周囲流体からの反力によって移動を行なう停空可能な浮上移動装置であり、上記制御部において流体・構造連成解析を用いたfree-flightシミュレーションを行ない、この結果を用いて、制御を行なうことを特徴とする。

この構成により、羽の変形、およびその相乗作用による浮上力変化に正確に対応して飛行可能な浮上移動装置が実現される。

また、本発明における浮上移動装置は、通信手段を有し、上記流体・構造連成解析を、物理的に独立した計算部にて行なうことを特徴とする。

これにより、流体・構造連成解析に必要なメモリ等のハードウェアの質量を浮上移動させる必要が無くなるため、浮上移動が容易になる。

また、本発明における浮上移動装置は、形状計測手段を有し、これにより計測された周囲障害物の形状を、上記流体・構造連成解析における境界条件として与えることを特徴とする。

これにより、周囲物体からの地面効果による流体力の変化を反映して正確に飛行可能な浮上移動装置が実現される。

以下、本発明の他の局面の浮上移動システムについて説明する。なお、他の局面の浮上移動システムは、前述の浮上移動システムまたは浮上移動装置のそれぞれと実質的に同様の構成を有していれば、前述の浮上移動システムまたは浮上移動装置のそれぞれによって得られる効果と同様の効果を得ることができる。

本発明の他の局面の浮上移動システムは、図１、図２および図２５に示すように、胴体部（支持構造９）と、胴体部に取り付けられた羽部（羽４）と、羽部を駆動する駆動部（アクチュエータ２、ベアリング３１,３２、伝達軸３３、ストッパ３５１,３５２、突起３４、ストッパ移動機構３６１,３６２）と、駆動を制御する制御部（制御装置５）と、周囲環境の形状を測定する環境センサ（近傍形状計測センサ５２）と、胴体部、羽部および駆動部のうち少なくともいずれか１つの運動に関する物理量を測定する運動センサ（６軸加速度センサ５１）とを備えた浮上移動システムである。

制御部は、図２６のフローチャートに示すように、羽部および胴体部のそれぞれの構造に関する数値モデル、周囲環境の流体に関する数値モデル、および、駆動部が羽部を駆動する駆動態様と駆動部がその駆動態様で駆動しているときに胴体部にかかる力との関係を示す関係データ（Ｆ２とｄの関係：図２０参照）とを有する記憶部（図２６のＳＳ３のＲＯＭ５３）を有している。

また、制御部は、図２６のフローチャートに示すように、周囲環境の形状のデータ（近傍形状計測センサ５２により取得されたデータ）、運動に関する物理量のデータ（６軸加速度センサ５１により取得されたデータ）、構造に関する数値モデル（羽モデルおよび胴体モデル、または構造モデル）、流体に関する数値モデル（流体モデル）、ならびに、関係データ（ホバリング時のＦ２とｄの関係）を用いて、流体の挙動を無視して、駆動部の最適駆動方法を決定する手段（ＳＳ４）を含んでいる。

また、制御部は、図２６のフローチャートに示すように、最適駆動方法のデータ（ＳＳ４のｄＬ、ｄＲ）を用いて、流体の挙動を考慮して、流体・構造連成有限要素法による羽ばたき飛行のfree-flight解析の結果に基づいて、最適駆動方法のデータを修正し、修正された駆動方法のデータを取得する手段（ＳＳ５、ＳＳ６、ＳＳ７）と、修正された駆動方法のデータを用いて、駆動部を駆動する手段（ＳＳ８、ＳＳ９）とを含んでいる。

また、流体・構造連成有限要素法によるfree-flight解析における流体・構造連成解析においては、胴体および駆動部の構造に関する数値モデルを用いることなく、羽部の構造に関する数値モデルのみが、構造に関する数値モデルとして用いられる。

また、浮上移動システムは、胴体、羽部、駆動部、制御部、環境センサ、および運動センサの全てが１つのロボットに搭載されたスタンドアロン型の羽ばたき飛行ロボットであってもよい。

一方、浮上移動システムは、胴体、羽部、駆動部、制御部、環境センサおよび運動センサを含む羽ばたきロボットと、羽ばたきロボットとは別個独立して外部に設けられた制御部とを有し、羽ばたきロボットが、制御部と通信することによって、羽ばたき飛行を実行するものであってもよい。

また、形状センサは、羽ばたきロボットに設けられた複数のＣＣＤ（Charge Coupled Device）であり、運動センサは、羽ばたきロボット自身の加速度を測定するセンサであることが望ましい。

また、前述の関係データは、ホバリング時の羽ばたき方を基本データとし、基本データに対して変更を加えることによって他の羽ばたき方が特定される形式で記憶部に格納されていることが望ましい（図２０参照）。

また、制御部は、周囲環境の形状のデータ、運動に関する物理量のデータ、構造に関する数値モデル、流体に関する数値モデル、および、関係データを用いて、羽ばたきロボットにとっての最適経路を算出し（図２６のＳＳ２のパスプランニング）、最適経路のデータに基づいて、最適駆動方法のデータを算出する（図２６のＳＳ５〜ＳＳ７）ことが望ましい。

また、流体・構造連成有限要素法による羽ばたき飛行のfree-flight解析においては、修正された駆動手法によって駆動部が駆動されたときの羽ばたきロボットのシミュレーション上の飛行経路が１回または複数回算出され、シミュレーション上の飛行経路と最適飛行経路との差が所定値以下であれば、制御部が、修正された駆動手法を用いて、駆動部を制御することが望ましい（図２６のＳＳ５〜ＳＳ９）。

前述の運動に関する物理量には、羽ばたきロボット自身の移動の速度ベクトルが含まれており、制御部は、羽ばたきロボット自身の移動の速度ベクトルとは逆向きで大きさが同じ速度ベクトルを羽ばたきロボットの周囲流体の対地速度ベクトルに加えて、流体・構造連成有限要素法による羽ばたき飛行のfree-flight解析を実行することが望ましい。

また、流体・構造連成有限要素法による羽ばたき飛行のfree-flight解析においては、駆動部の力の出力態様または変位の出力態様のデータを用いて、構造に関する数値モデルの節点における境界条件として外力または変位のデータが作成されることが望ましい。

また、流体・構造連成有限要素法による羽ばたき飛行のfree-flight解析においては、胴体部の１つの節点と胴体部の１つの節点と接合される羽部の１つの節点とに、同一節点変位、または、方向が逆向きで互いに大きさが等しい節点外力が与えられるものとして、解析が行なわれることが望ましい。

より具体的には、前述の修正された駆動方法のデータを取得する手段は、図１４に示すように、次のような（１）〜（９）に示す手段を有している。

（１） 運動に関する物理量のデータを用いて、第１所定時刻における前記胴体部の第１絶対位置を算出する手段（図１４のＳ１の胴体の位置Ｇ）。

（２） 運動に関する物理量のデータを用いて、第１所定時刻における胴体部の第１姿勢を算出する手段（図１４の姿勢Ｅ）。

（３） 最適駆動方法のデータを用いて、仮に、第１所定時刻と第１所定時刻よりも後の第２所定時刻との間の胴体部に対する羽部の相対変位を算出する手段（図１４のＳ１の羽相対位置ａ）。

（４） 第１絶対位置のデータ、第１姿勢のデータ、および羽部の相対変位のデータを用いて、第２所定時刻における胴体部に対する羽部の絶対位置を算出する手段（図１４のＳ１）。

（５） 羽部の絶対位置のデータを用いて、胴体部にかかる並進力および胴体部にかかるトルクを算出する手段（図１４のＳ２）。

（６） 並進力のデータおよびトルクのデータを用いて、第２所定時刻における胴体部の第２絶対位置および第２所定時刻における胴体部の第２姿勢を算出する手段（図１４のＳ３）。

（７） 第１絶対位置のデータおよび第１姿勢のデータと第２絶対位置のデータおよび前記第２姿勢のデータとの差を算出する手段（図１４のＳ４）。

（８） 差が所定値以下であるか否かを判別する手段（図１４のＳ４）。

（９） 差が所定値以下である場合に、羽部の相対変位のデータを修正された駆動方法のデータとして用いて、駆動部を制御する手段（図２５の超音波モータドライバ）。

本発明における浮上移動システムによれば、流体・構造連成解析を用いたfree-flightシミュレーションを用いて経路予測することにより、周囲障害物からの地面効果その他の流体の影響を正確に勘案し、羽ばたき駆動手法を決定することができるため、障害物周辺からの非定常気流などによる衝突などの問題を回避することができる。このため、災害救助などの狭所作業に適した羽ばたき飛行ロボットが実現できる。

本発明の一実施の形態の羽ばたき浮上移動装置について、図１〜図１２を用いて説明する。

（全体の構成）
まず、本実施の形態の浮上移動装置の全体の主要な構成について、図１および図２を用いて説明する。なお、本実施の形態の浮上移動装置では、特に断らない限り、説明の簡便のため、羽を駆動する主要な部分は左右で面対称の一対であるものとする。したがって、以後においては、浮上移動装置の左半分についてのみ説明を行ない、左半分に対して面対称である構成要素が右半分について配されているものとする。但しこれは説明の簡便のためであり、本発明の浮上移動装置においては、左右対称であることは、必要な条件ではない。例えば、当該構成要素が３組以上配されている構成や、１組のみ用いる構成も可である。

本実施の形態の羽ばたき浮上移動装置１は、図１および図２に示すように、支持構造９上に配されたアクチュエータ２、このアクチュエータ２によりベアリング３１，３２を介して駆動される伝達軸３３、および伝達軸３３により駆動される羽４を主要な構成要素とする。

伝達軸３３には、伝達軸３３の並進運動（アクチュエータ２の回動：伝達軸３３が扇型状の軌跡を描く運動）を用いて伝達軸３３に伝達軸３３が延びる方向の軸まわりの回動（羽４を伝達軸３３が延びる方向に見たときに羽４が扇型状の軌跡を描く運動）を与える機構として突起３４が固定されている。突起３４は、アクチュエータ２の回動運動の両端付近において、ストッパ３５１,３５２に接触する。突起３４はストッパ３５１,３５２の接触面上を滑りながら移動することにより、伝達軸３３に回転運動を発生させる。突起３４は各々ストッパ３５１，３５２に応じて２つの独立した接触部位を有し、図２では単一の構成要素として描かれているが、２つの突起を配したものと等価である。

ストッパ３５１，３５２は支持構造９に対する相対的な位置を、支持構造９に配されたストッパ移動機構３６１，３６２によって変更することができる。ストッパ移動機構３６１，３６２は、突起３４がストッパ３５１，３５２に接触およびこの表面上を移動する際の反力を上回る、ストッパ位置保持力を有する。

また、ベアリング３１，３２は、羽４がその翼弦長方向を図１における鉛直下向きに向けた状態を挟んで±３０度の回転制限角の範囲のみに回転可能である、伝達軸回転制限機構を兼ねている。

なお、この回転制限角は説明の簡便のための一例であり、この数値そのものは本発明の範囲を制限するものではない。例えば羽４が剛体であるなら、回転制限角は±８０度程度であることが望ましい。

以上の構成により、アクチュエータ２が、アクチュエータ２の回動中心軸回りに往復回動運動をすると、その両端付近において、突起３４がストッパ３５１，３５２に接触することにより、伝達軸３３に、伝達軸３３が延びる方向の軸である回動中心軸回りに回動する運動が与えられる。その後、往復回動運動の方向が逆転し、突起３４がストッパ３５１，３５２から離れ、伝達軸３３の軸回り回転角は、伝達軸回転制限機構により固定されるため、羽４は周囲流体に対して一定の迎え角のまま移動する。以上により、羽４は、アクチュエータ２の往復回動運動中心近傍では羽４にかかる、一般的な航空機の翼と同様の揚力によって、また、当該往復回動運動両端付近では羽４の回転運動によって生み出される回転揚力によって、浮上移動装置１に、鉛直上向きの力を与え、これにより、浮上移動装置１は鉛直上向きに浮上する。

また、詳細については後述するが、本実施の形態の浮上移動装置１は、ストッパ移動機構３６１，３６２を用いて、左右それぞれの羽４のアクチュエータ２の回動の振幅の中心を、別個独立して進行方向前側または進行方向後側にずらすことにより、進行方向を左右のいずれかに変更すること、ならびに、前進および後退のいずれかの状態に変更することが可能である。なお、浮上移動装置１の前後方向に延びる軸回りの回転動作は、左右の羽４の、アクチュエータ２における往復回動運動の振幅の大きさおよび周波数を異ならせることにより、実現されるものとする。

アクチュエータ２およびストッパ移動機構３６１，３６２それぞれは制御装置５により制御が決定される。

制御装置５には、６軸加速度センサ５１より浮上移動装置１の加速度および角加速度が、近傍形状計測センサ５２より浮上移動装置１の近傍の障害物の配置が与えられる。なお、浮上移動装置１の加速度は、後述する流体・構造連成解析によるfree-flightシミュレーションに用いるに適した計測方法で計測された値が望ましい。つまり、本実施の形態においては重心で計測することが望ましいが、部品のレイアウト等の都合上、当該重心位置が６軸加速度センサ５１を配置するに適さない位置である場合がある。しかし、支持構造９に対して、大きく形状もしくは位置の変化する部位は、軽量な羽部のみであり、それ以外は支持構造９に固定されているため、質量の殆どは支持構造９に固定されている。そこで、６軸加速度センサ５１は支持構造９に固定しておき、この固定位置と浮上移動装置１の羽部を除いた部位の重心との位置関係は固定されているので、６軸加速度センサ５１で計測された加速度を上記位置関係に基づき変換することで、浮上移動装置１の加速度を得ることができる。

制御装置５にはＲＯＭ５３より制御テーブルや後述する流体・構造連成解析によるfree-flightシミュレーションの基礎モデル情報、例えば初期の解析空間メッシュの形状や、羽のメッシュ情報、物性値などが与えられる。

制御装置５には、後述する流体・構造連成解析によるfree-flightシミュレーション演算部５５２および羽駆動演算部５５１が内包されており、前者の結果を得て後者が最適な羽駆動方法の算出を行なう。また、後者による羽駆動算出結果の良否を前者シミュレーションにより判断することができる。さらに、これらプランニングと判断を反復して行なうことにより最適経路が求められる（パスプランニング）。

これらの計算結果はＲＡＭ５４に格納される。また、ＲＡＭ５４より得られたこれまでの計算結果等により、制御装置５は次のステップでのパスプランニングを行なうことができる。

制御装置５により算出された結果に基づき、制御装置５は、超音波モータドライバ５５に制御信号を送る。これは超音波モータドライバ５５により駆動信号に変換され、アクチュエータ２に与えられ、羽が駆動される。例えば、制御装置５はＰＷＭ周波数を求め、これに対応したパルスを超音波モータドライバ５５に与え、これに従い上記超音波モータをＰＷＭ駆動する等の手法が考えられる。

また、アクチュエータ２およびストッパ移動機構３６１，３６２、および制御装置５は、電源６より供給される電力によって稼働する。

（支持構造）
次に、支持構造９について説明する。支持構造９は、支持構造としての機能を損なわない範囲内の質量であって、より軽量であることが望ましい。この範囲は、従来のエンジニアリングの技術で求めることができる。たとえば、軽量化のあまり剛性が低下し過ぎてしまい、羽４を駆動させるべきエネルギが支持構造９の振動となって散逸することがないこと等が、その範囲を決定するための基準となる。本実施の形態の浮上移動装置１においては、軽量なカーボングラファイトを支持構造９に用いる。

（羽駆動部）
次に、羽の駆動部について、図２〜図７を用いて説明する。

（アクチュエータ）
まず、アクチュエータ２について、図２を用いて説明する。

本実施の形態の浮上移動装置１においては、アクチュエータ２として、超音波モータと一般に呼ばれているモータを用いる。本実施の形態においてアクチュエータ２として用いる超音波モータは、図２に示されるように、支持構造９に固定されているステータ部２１に対してロータ部２２を回動させることができる。

（ベアリングおよびその回転制限機構）
ベアリング３１,３２は、その外周部３１ａ，３２ａがロータ部２２に固定されている。ベアリング３１,３２は、外周部３１ａ，３２ａに対して内周部３１ｂ，３２ｂが円滑に回転できるようになっており、内周部３１ｂ，３２ｂには伝達軸３３が固定されている。したがって、アクチュエータ２は、伝達軸３３が延びる方向の回動軸を含む面であって、アクチュエータ２の主表面と平行な面において、伝達軸３３をアクチュエータ２の回動中心点回りに回動させることが可能である。

また、このベアリングの可動回転範囲は約６０度に制限されている。これは機械的な当たり面を外周部３１ａ，３２ａおよび内周部３１ｂ，３２ｂに設けることで簡便に実現できる。なお、この実現方法は、本浮上移動機構の性能を損なわない限り本発明の機能を損なうものではなく、上記の方法に限るものではない。そこで、図２および図３では説明の簡便のためこれを黒円弧状のマークで抽象的に表わしている。

（ストッパおよび突起による伝達軸回転機構）
次に、突起３４およびストッパ３５１，３５２について図２〜図５を用いて説明する。

説明の簡便のため、図２に示される、アクチュエータ２の円盤と相似形の円を底面に有する円柱の周面における展開図を図３に示す。

突起３４は、ベアリング３１，３２に設けられた伝達軸回転制限機構により、±３０度の範囲で、伝達軸３３の中心軸を回転中心にして自由に回転運動できる。図３に示されたのは回転が−３０度、つまり最も時計回りに位置する状態である。

この図３において、伝達軸３３は右方向に運動するとする。この際、羽４に対しては、図３において左上方向に流体力がかかっているため、伝達軸３３に対して羽４が下方にあることから、伝達軸３３には時計回り方向のトルクが羽４より与えられる。しかし、伝達軸回転制限機構により、伝達軸３３は−３０度より時計方向に回転しないため、図３に示す状態では、羽４にかかる流体力は、伝達軸３３を介して、浮上移動装置１に、図３において左上方向への力として伝達される。

さらに図３の状態より伝達軸３３が右に移動すると、突起３４はストッパ３５１に接触した後、更なる伝達軸３４の右方向への運動によって、ストッパ３５１の面上を滑りながら移動。なお、この接触の態様はアクチュエータ２の駆動態様、ストッパ３５１と突起３４の図３における上下および左右方向の位置関係、ストッパ３５１の移動の有無、およびストッパ３５１の形状により決定することができるが、本実施の形態においては、接触時の衝撃を和らげ、かつ、伝達軸３３の回転が緩慢に始まるよう、突起３４に対してほぼ０の接触角を有する、２次以上の曲線であるとする。

さらに、伝達軸３３そのものとストッパ３５１との衝突を防ぐため、突起３４は、図３に示されるように、伝達軸３３より離して、上方左右に配されている。

伝達軸の運動が右端に達し、伝達軸３３は左方向に運動を始めると、羽４にかかる流体力は右上方向となり、伝達軸３３には反時計回りのトルクがかかるため、伝達軸回転制限機構により、伝達軸回転角は＋３０度に至り、これを維持したまま、図３における左方向への運動を行なう。その後、今度はストッパ３５２と接触し、上記のストッパ３５１との接触により起きた運動の、図３における左右対称の運動を行なう。

本実施の形態においては、図４における羽の運動を意図したが、実際には、さらに羽４の変形が起きるため、図４に等価な羽４の運動を実現するには、伝達軸回転制限機構の回転制限範囲を、羽４が剛体である場合より小さく設定するのが最も簡便である。結果的に、羽４の運動は図５に示される運動となる。

なお、本実施の形態の浮上移動装置１における突起３４およびストッパ３５１，３５２は、突起３４の機能を果たせるのであれば、必ずしも図３に示した形状、位置である必要はない。例えば伝達軸３３をベアリング３１よりさらに羽４の反対方向に延長した先に、伝達軸３３の軸に対して、図２に示す方向と逆方向に突起３４が配され、これに対応した位置にストッパ３５１，３５２が配されていたとしても、同様の羽４の運動が実現できる。

なお、本実施の形態の浮上移動装置においては、アクチュエータ２の回動角は±４５°である。また、±４５°のアクチュエータ２の回動角によって生じる伝達軸３３が延びる方向の回動中心軸回りの回動角は±３０°と設定されているものとする。

（ストッパ移動機構）
次に、ストッパ移動機構３６１，３６２について、図３および図７を用いて説明する。

ストッパ移動機構３６１，３６２は、支持構造９に設置されており、ストッパ３５１，３５２がアクチュエータ２の外周上の点の回動の軌跡と相似形の円弧状の軌跡を描くようにストッパ３５１，３５２の位置を移動させる機能を有する。このストッパ３５１，３５２の位置の移動は、図３においては、ストッパ３５１，３５２が水平方向（左右方向）に移動することに相当する。

このストッパ３５１，３５２の位置を浮上移動装置１の前方または後方にずらし、かつ、羽４のアクチュエータ２の回動領域を浮上移動装置１の前方または後方側にずらすことで、羽ばたき運動の振幅を変更したり、回転中心を変更することができる。これにより、後述する飛行方向の転換などの制御を簡単に行なうことができる。

また、ストッパ移動機構３６１，３６２によって、伝達軸３３の伝達軸回り回転角を制御することが可能である。たとえば、突起３４がストッパ３５１，３５２に接触している間に、伝達軸３３の移動方向と同一方向にストッパを動かせば、伝達軸３３の伝達軸回り回転は、ストッパの移動が無い場合に比べ緩慢になる。

ストッパ移動機構３６１，３６２における、ストッパ推進力発生機構は、必要なトルクや質量および消費エネルギなどの条件を満たすなら特にその構成に制約はない。そこで、本実施の形態の浮上移動装置のストッパ移動機構３６１，３６２においては、図７に示すように、応答性に優れた超音波リニアアクチュエータ３７を用いる。図７は、図３の展開図において、紙面に垂直にストッパ移動機構３６１を切断したときのストッパ移動機構３６１およびストッパ３５１の断面図である。なお、ストッパ３５２およびストッパ移動機構３６２についてはこの鏡像の構成となる。なお、軽量化のため、ストッパ移動機構３６１とストッパ３５１との接触部であって超音波リニアアクチュエータが存在しない部分には、潤滑性に優れたテフロン（Ｒ）ベアリング３８が用いられている。また、空気から羽４へ加えられる反力に起因して超音波リニアアクチュエータ３７に加えられる力の影響を低減するため、超音波リニアアクチュエータ３７は伝達軸３３が延びる方向をその法線に有する面上にその面に沿って平行に延びるように設けられている。

また、ストッパ３５１，３５２は、突起３４とストッパ移動機構３６１，３６２の干渉を避けるため、ストッパ移動機構３６１，３６２より突起３４側に配されている。

（羽部）
続いて、羽４について図８および図１２を用いて説明する。

［羽の主要構成］
図１２（ａ）に示されるように、羽４は前縁４１、内側枝部４２１、中央枝部４２２、羽面縁取り部４４１、連結部４３、羽面部膜４４２、枝部膜４５からなり、ステンレスからなる梁状部材である前縁部４１、内側枝部４２１、中央枝部４２２、羽面縁取り部４４１、連結部４３、に対して、フィルム状部材である羽面部膜４４２、枝部膜４５が張られてなる構成である。また、連結部４３は静止状態に置いて、上方にたわみを有している。

羽面部膜４４、枝部膜４５は室内用軽飛行機などに用いられているポリエチレン薄膜や、マイクロフィルムなどの軽量なものを用いることが望ましい。しかし、この膜の性質の詳細は本発明の実施要件に関係しないため、本実施の形態では説明の簡便のため、この膜の性質の詳細は略する。すなわちこの膜の剛性はステンレスに比べ遙かに小さいが、極端に変形したり、破れたりして空力特性に影響を与えるのでないとする。

［前縁部］
前縁部４１は、厚さ２０μｍであり、長さ３０ｍｍ、幅は２.５ｍｍである。その断面図である図１２（ｂ）から分かるように、山と谷をそれぞれ１つずつ有する浪板状の構造を有している。山と谷の間隔は１ｍｍ、山の高さと谷の深さはそれぞれ０．５ｍｍである。このような浪板構造を、昆虫の羽における類似の形状の呼称を用いて、コルゲーションと称する。また、コルゲーションの、羽駆動部に近い一端は、並進自由度、回転自由度共に伝達軸３３に拘束されている。以後、このように並進および回転量の自由度は、共に拘束されている状態を、単に拘束されていると称する。このコルゲーションにより、前縁部４１は、その長手方向に高い剛性を有する。また、前後方向については通常の平板の曲げと同じ剛性である。

［内側枝部］
内側枝部４２１は厚さ４０μｍ、幅約０．２ｍｍ、長さ約１０ｍｍで、その前縁に近い側の一端は、伝達軸３３に対して拘束されている。また、もう一端は、連結部４３に拘束されている。

上記の構成により、内側枝部４２１は伝達軸３３にその運動を拘束され、また、前縁部４１を除く、羽４の他の外周部分に比べて厚さが２倍あるため変形が小さく、連結部４３を介して羽面縁取り部４４１の一部を拘束することで、羽面部膜４４２の迎え角を拘束する役割を有する。

［中央枝部］
中央枝部４２２は厚さ２０μｍ、幅約０．２ｍｍ、長さ約１０ｍｍで、その前縁に近い側の一端は、前縁部４１に対して拘束されている。また、もう一端は、連結部４３の内側枝部に拘束されていないもう一端および、羽面縁取り部４４１の前縁部4１に拘束されていないもう一端に拘束されている。

上記の構成により羽面部膜４４２を保持している。

［連結部］
連結部４３は厚さ２０μｍ、幅約０．２ｍｍ、長さ約１５ｍｍで、内側枝部に拘束されている。また、連結部４３は図１２（ａ）に示されるように、上方に約２ｍｍの膨らみを有している。

［羽面縁取り部］
羽面縁取り部４４１は、厚さ２０μｍ、幅約０．２ｍｍ、長さ約２０ｍｍで、中央枝部４２２と前縁部４１にそれぞれ一端を拘束されている。

［羽の変形の様子］
羽４の変形については、以下に挙げる現象が支配的である。先ず、前縁部４１がその長手方向の曲げに対して高い剛性を有し前後方向のみ変形する。また、本実施の形態において羽４は図６に示されるような運動をするため、羽の変形は殆ど前後方向への曲げとして生じる。また本実施の形態においては羽4の運動速度が大きい駆動部より遠い部分の変形が大きい。内側枝部４２１は伝達軸３３に拘束され、変形も小さい。連結部４３は膨らみの曲率が変化する曲げ変形については剛性が小さく、その全長が変化する引っ張りについては剛性が高い。

以上の変形を示すのが、図８に示される加重と変位の関係である。図８は前縁部４１の羽駆動部に近い側の一端を固定した際の、図１２（ａ）における点Ｅの荷重変位関係である。連結部４３が曲げ変形を生じ、主に膨らみの曲率が変わる範囲の荷重では非常に変形が大きく、この曲率が小さくなり、連結部４３が引っ張り変形に転じると、殆ど変形しなくなることが分かる。

羽４が左羽であるとした場合の、打ち下ろしにおける羽４の形状を図９に示す。

［羽の変形の効果］
図９に示されるように、羽４は、その移動速度が小さい内側部分は変形が小さいため迎え角が大きく、逆に移動速度が大きい外側部分は変形が大きいため迎え角が小さい、という状態を、羽の変形によって、打ち下ろし、打ち上げの両方においてとることができる。

このため、打ち下ろし、打ち上げ共に、図１１に示される羽ばたき方のように、伝達軸３３の捻り角のみを制御することで、羽の各部位における迎え角を適切な値にすることができる。本実施の形態では、羽の移動速度が小さい内側では迎え角を大きくとって揚力を増し、羽の移動速度が大きい外側では迎え角を小さくとって失速を防止すると共に、流体の抵抗を受けにくくしている。

また、打ち上げと打ち下ろしの間で、羽の変形が反転する際は、浮上力が低下してしまうため、この時間はなるべく短いことがのぞましい。このためには、羽を前後方向に変形しやすくしておくことが有効である。ところが、羽を前後方向に変形しやすくすると、上記迎え角が適切に取れないと言う問題がある。そこで本実施の形態においては、図８に示すように、連結部４３の曲げ剛性と引っ張り剛性の違いを利用して、非線型な荷重変位関係を実現することで、上記迎え角に達すれば殆ど変形しないが、それまでは大きく変形する羽によって、この相反する問題を解決している。

（羽形状の他のバリエーション）
本実施の形態においては、作成可能な範囲でなるべく効率の良い羽を用いたが、効率の悪化を許容できるのであれば、本実施の形態に示された羽と一部の機能を同じくする、簡略化された羽を用いることも可である。この場合、羽の作成においてコストダウンが可能であるというメリットがある。

（飛行制御）
続いて、飛行制御手法について図６〜図２４を用いて説明する。

（基礎的飛行制御）
先ず、飛行の制御方法について図６、図１０、および図１１を用いて説明する。なお、ここに示した飛行の制御方法は第一義近似的な例示であり、後述する流体・構造連成解析を用いた羽駆動手法算出の反復計算における初期値であると解釈すべきである。

また、ある羽ばたき飛行の状態を生み出す羽ばたき方を変更する手法は、以下に記述された手法以外にも存在するが、そういった羽ばたき方を採用した際にも本実施の形態における構成およびアルゴリズムを用いることができる。

羽ばたき飛行により、３次元空間をくまなく移動するには、前進後退、左右への旋回、上昇および下降の３つの運動要素が実現される必要がある。以下、各運動要素およびその基本動作となるホバリングを実現する手法について述べる。

（ホバリング）
本実施の形態の浮上移動装置においては、図１０に示す羽ばたき方でホバリングが可能である。突起３４およびストッパ３５１，３５２それぞれの形状は、ストッパ移動機構３６１，３６２によりストッパ３５１，３５２が移動しない状態で、アクチュエータ２が伝達軸３３に±４５°の回動をさせる際に、この図１０に示す運動が実現されるよう設計されている。

このため、図６に示されるように、左右のアクチュエータ２を回動の振幅の中心が浮上移動装置１のほぼ真横方向になる状態で、振幅の中心に対して前後に±４５°の振幅で羽４を回動させることでホバリングが実現される。つまり、左のアクチュエータ２の回動の振幅中心と左の羽部４の伝達軸３３が延びる方向の先端を結ぶ線と、右のアクチュエータ２の回動の振幅中心と右の羽部４の伝達軸３３が延びる方向の先端を結ぶ線とが、ほぼ一直線上になるように、左右の羽４それぞれが前後方向に往復運動すれば、浮上移動装置１はホバリングを行なう。

説明の簡略のため、以後この羽ばたき方を、前後略対称な羽ばたき方、もしくは、オフセットゼロの羽ばたき方と称する。

（前進および後退）
上述の羽４の前後方向の±４５°の回動により、羽４には、鉛直下向きの流れ以外に、伝達軸３３が延びる方向において、アクチュエータ２側から羽４の先端側へ向かうように、図１１に示すような流れが発生する。この流れを利用することで羽ばたき浮上移動装置１に前後方向の移動をさせることができる。

たとえば、ストッパ移動機構３６１，３６２を各々前方もしくは後方に移動させ、アクチュエータ２の回動の振幅の中心を前側または後側にずらすことで、羽４の前後方向の往復運動である羽ばたき運動を前側または後側に偏らせる。それにより、上述の伝達軸３３が延びる方向において、アクチュエータ２側から羽４の先端側へ向かう流れを浮上移動装置１の後方または前方へ向けることができる。そのため、浮上移動装置１は推進する力または後退する力を得ることができる。左右の羽４について等しく羽４の根元から先端へ向かう流れを等しく後方へ向ける羽ばたき方に変更すれば、羽ばたき浮上機構１は前進する。また、伝達軸３３が延びる方向においてアクチュエータ２側から羽４の先端側へ向かう流れを、左右の羽４について等しく、浮上移動装置１の前方へ向ける羽ばたき方に変更すれば、浮上移動装置は後退することが可能となる。この様子を図１１に示す。
以後これを、説明の簡便のため、上記アクチュエータ２の回動の振幅の中心を前方向にｘ度ずらしたものを、オフセットｘ度の羽ばたき方と称する。

（左右への旋回）
上記の前進及び後退を、左右の羽について別々に行なえば、左右への旋回を実現できる。
例えば左羽を前進させれば右側に、右羽を前進させれば左側に旋回を行なうことができる。

またに、上述の羽４の表面において、伝達軸３３が延びる方向においてアクチュエータ２側から羽４の先端側へ向かう流れを用いることで、浮上移動装置１を左右のいずれかへ旋回させることが可能である。

たとえば、左の羽４のみ上述の伝達軸３３が延びる方向においてアクチュエータ２側から羽４の先端側へ向かう流れを後方に向ければ、羽ばたき浮上機構１は右へ旋回しながら前進する。逆に、右の羽４のみ上述の伝達軸３３が延びる方向においてアクチュエータ２側から羽４の先端側へ向かう流れを後方に向ければ、羽ばたき浮上機構１は左へ旋回しながら前進する。

また、左の羽４のみ上述の伝達軸３３が延びる方向においてアクチュエータ２側から羽４の先端側へ向かう流れを前方に向け、かつ、右の羽４のみ上述の伝達軸３３が延びる方向においてアクチュエータ２側から羽４の先端側へ向かう流れを後方に向ければ、羽ばたき浮上装置１は、その上方から見て反時計回りにその場で回転することができる。
以後これを、上記オフセットを左右別々に扱い、右羽オフセットｘ度、左羽オフセットｙ度の羽ばたき方、のように称する。

（浮上および降下）
浮上および降下を行なう最も単純な方法は、アクチュエータ２の駆動を変更し、羽ばたき周波数を上昇または低下させる手法が挙げられる。

また、前述の手法以外の浮上および降下を行なう手法としては、羽ばたきのストローク（振幅）を大きくする手法、すなわちアクチュエータ２による伝達軸３３の前後方向の回動の振幅を大きくする手法が挙げられる。

すなわち、アクチュエータ２の回動の振幅を大きくして、かつ、ストッパ移動機構３６１，３６２によって、ストッパ３５１，３５２を、伝達軸３３の移動方向と同方向に、伝達軸３３の回動の振幅に相当するだけ移動させる。たとえば、アクチュエータ２の回動を±４５°から±５５°に変更したとすると、ストッパ移動機構３６１，３６２によってストッパ３５１，３５２を、各々＋１０度、−１０度移動させる手法が挙げられる。

これらを総合して、以後、右羽オフセットｘ１度、右羽振幅拡大ｘ２度、左羽オフセットｙ１度、左羽振幅拡大ｙ２度、の羽ばたき方、というように称する。

（姿勢変化への対処）
上述の運動変更により、浮上移動装置１の姿勢は変化する。姿勢の変化が望ましくない場合には、浮上移動装置１の重心を、伝達軸３３の位置より低くしておけばよい。これによって、ある程度の時間ホバリングすることによって、左右の伝達軸３３の中心直下に浮上移動装置１の重心が位置する姿勢に安定する。

（流体・構造連成解析手法）
（連成解析手法）
本実施の形態における流体・構造連成解析には、文献“Analysis of fluid-structure interaction problems with structural buckling and large domain changes by ALE finite element method.” （Qun Zhang, Toshiaki Hisada，Comput. methods Appl. Mech. Engrg. 190, 6341-6357 （2001））に示されている、ＡＬＥ法 （Arbitrary Lagrangian-Eulerian method）に基づく流体・構造連成有限要素法解析を用いる。以後これを、fluid-structure interaction と finite element methodの頭文字をとって、ALE-FSI-FEMと称する。

ＡＬＥ法は、任意に運動する点から観測した流体の挙動を定式化することが可能であり、これと有限要素法を組合せることにより、任意に変形する有限要素（メッシュ）を用いて流体挙動の解析を実現している。さらに、流体・構造境界では両者の運動が一致していることおよび、当該境界に働く力については線形和が成り立つことから、上記ＡＬＥ法による流体挙動に関するメッシュ方程式と、構造挙動に関するメッシュ方程式は、当該境界における共有節点を介して直接マージ可能であり、これを解くことにより流体・構造連成問題の解を得ることができる。

さらに、本実施の形態においては、空間のメッシュについては、弾性体の構成則を仮定したアルゴリズムに基づき、羽と、解析領域最外周に与えられた境界条件に従う静的問題解くことで求められる。これは、解析空間が仮想の弾性体で満たされており、弾性体内部に固定された節点が、羽と解析領域外周の運動に応じて起きる弾性体の変形に従って移動する状態に等しい。これはあくまでメッシュ形状を求めるための手法であり、上記弾性体の性質は、メッシュの破綻が起こらない限りは、流体及び構造の挙動算出には影響を与えない。このため、同手法では多くの場合、解析精度が望まれる度合いに応じて上記の仮想弾性体の剛性を高くし、メッシュ形状が角に歪むことで計算精度が低くなることを防いでいる。

（free-flight解析手法）
本実施の形態では、文献「流体・構造連成有限要素法による羽ばたき飛行のfree-flight解析」（濱本、太田、原、久田、日本ロボット学会第21回学術講演会予稿集、2003）によるfree-flight解析手法を用いる。その概略は以下の通りである。

（羽と胴体の相互作用の扱い）
最も忠実な手法においては、アクチュエータの挙動も含め、制御信号から機械的運動までの全ての現象について物理的シミュレーションを行なうべきであるが、これは、計算コストがかかる上に、アクチュエータの挙動は人工的に取り扱いが容易であるべく設計されている場合が多いので、本実施の形態ではこのアクチュエータの部分は、ある制御入力に対してある定められたテーブルに基づく変位もしくは力を出力するものとして、ブラックボックス化して扱う。

このため、羽４とそれ以外の構造部位（以後、胴体９４と称する）は、このブラックボックスを介して力ないしは変位をやりとりするため、羽４と胴体９４は、モデル上では不連続であり、計算上、変位もしくは節点力をそれぞれの節点に加えることで、アクチュエータにより羽４と胴体９４の双方が移動し、相互に力を及ぼしあう状態を実現する。よって、アクチュエータの出力は、解析モデルにおけるある（複数可能）節点に対する変位境界条件もしくは節点力として、上記ALE-FSI-FEMに与えられる。この様子を図１３に示す。

本実施の形態においては、超音波モータの特性上、制御信号からアクチュエータの変位が直接得られるために、アクチュエータの位置を基準とする相対変位にて、上記節点の変位境界条件を与える。

また、本実施の形態においては、胴体９４が周囲流体の挙動に与える影響は、羽４の周囲流体の挙動に与える影響に比べ極端に小さいため、羽４のみについて、上記ALE-FSI-FEMを行ない、胴体９４については構造解析のみを用いることによりその挙動を算出する。さらに、支持構造９は軽量であり、その変形も小さいことから、胴体９４を、ある質量とモーメントを有する剛体として、計算をより簡略化する。

当然、これは合理的な計算を行なうための手法であり、他の手法も可である。例えば、胴体９４をモデルに含めてALE-FSI-FEMを行なう手法も可である。また、胴体９４の周囲流体に与える影響は小さいが、変形もしくは内部での質量移動等が飛行に対して無視できないほどの影響を与えると考えられる場合には、胴体９４はやはり流体・構造連成解析には含まないものの、その構造解析は緻密に行なうことで必要な精度を合理的に得る手法も可である。

（計算フロー）
以上の想定に基づき、ある与えられた制御信号から得られる浮上移動装置の挙動を、ALE-FSI-FEMを用いて計算のように算出する。なお、時刻ｔにおける運動状態、すなわち位置、速度および加速度の６成分は既知であるものとする。

ステップ１：時刻ｔの運動状態および胴体９４に加わる力から、時刻ｔ＋Δｔにおける胴体９４の位置を算出。

ステップ２：予め定められた制御信号より、時刻ｔ→ｔ＋Δｔにおけるアクチュエータ２の挙動を算出。

ステップ３：アクチュエータ２の挙動に基づき、羽４の胴体９４に対する相対変位を算出。

ステップ４：上記相対変位を元に、羽４の、空間座標系での運動を算出。

ステップ５：上記運動を用いて、流体・構造連成解析を行ない、時刻ｔ＋Δｔにおける胴体９４および羽４および周囲流体の状態を仮に算出。

ステップ６：上記計算により求まった、羽４よりアクチュエータ２を介して胴体９４に与えられる反力を算出。

ステップ７：時刻ｔより時刻ｔ＋Δｔにおいて、胴体９４に加わる力が線形に変化するものとして、時刻ｔ＋Δｔにおける胴体９４の位置を再度算出。すなわち、ステップ１で求めた時刻ｔ＋Δｔでの胴体の位置を修正する。

ステップ８：上記修正量があるスレッショルド以上であるなら、ステップ２に戻り、再度時刻ｔ→ｔ＋Δｔについての流体・構造連成解析を行なう。このように反復計算を行なうことにより、矛盾のない、時刻ｔ＋Δｔでの胴体９４と羽４、および周囲流体の状態が求まる。すなわち、上記修正量があるスレッショルド以下であるなら、時刻ｔ→ｔ＋Δｔの挙動は解析が完了したものとして、引き続き、ｔをｔ＋Δｔと読み替えて、引き続いて次の時刻での挙動を求める。

当然、これは、本実施の形態での想定における一例であり、胴体と羽の相互作用の扱いによってこれらの手法は異なる。

例えば、羽の駆動を胴体に対する強制変位でなく節点力で行なう手法も可であり、その場合、上記強制変位を当該節点に与えられる外力に読み替えればよい。

（本実施の形態でのアルゴリズム）
ここでは、実際に本発明者らが行なった解析をより具体的に示す。時刻ｔ→ｔ＋Δｔにおけるアクチュエータ２の挙動が決定しているので、羽４におけるある節点について、胴体に対する相対変位ａ_iの時刻歴は既に与えられている。時刻ｔでの空間座標系における上記節点の座標ｘについては、Ｇを胴体９４の重心位置、Ｅ_iを胴体９４に固定された直交規定ベクトルとして、

の式で与えられる。このａ_iがgiven parameterであり、他の値は以下のようにして求める。

先に説明したとおり、今回のシミュレーションでは、胴体９４の空力における作用は殆ど影響がなく、流体・構造連成解析を行うのは羽４のモデルのみで、胴体９４は羽４と力学的にインタラクションするのみであるとした。つまり胴体９４の運動方程式は剛体の力学に基づき、その質量をＭ、重力加速度をｇとして、

である。但し、θは方向余弦が(λ，μ，ν)である軸ｘ周りの角加速度で、この軸周りの胴体の慣性モーメントＩｘは、

であるとする。但し、各Ｉ_i(ｉ：ｘ、ｙ、ｚ）およびＩ_ij（ｉｊ：ｘｙ、ｙｚ、ｚｘ)については、ｄｍを微小質量として、

で与えられる。

今回の胴体９４のモデルについては各パラメータは表１のようになる。

また、羽はその根元の部位の節点における強制変位によって駆動されているので、このi番目の節点に働く節点力をｆ_i、重心に対するこの節点の位置ベクトルをｒ_iとすると、胴体が受ける並進力ＦとトルクＴは、

によって決まる。但し、総和記号は、羽４上の上記強制変位により駆動されている全ての節点についての総和を意味する。そして、以下の手順で胴体と羽を反復計算する。

まず、解析初期の過渡状態を避けるため、約２周期のtetheredシミュレーションを行う。ここまでは当然、胴体の速度及び加速度、角速度及び角加速度は０である。

上記準備解析完了時点で、胴体に掛かる並進力Ｆ及びトルクＴを算出する。

これにより時刻ｔと時刻ｔ＋Δｔでの加速度が求まるので、この間は線形に加速度が変化するものとすると、時刻ｔではその重心位置がＧ_tであった胴体の、時刻ｔ0＋Δｔでの胴体の位置が加速度の２階積分の結果、下記の式のごとく求まる。

姿勢についてはＧをθに読み替えればよい。

但し、初回の反復では時刻ｔ0＋Δtでの節点力は分かっていないので、加速度ｄ²Ｇ_t＋Δt／ｄｔ²は時刻ｔ0での値と変化しないとする。

こうして仮に求まった時刻ｔ0＋Δtでの胴体の位置と姿勢から、さらに時刻ｔ0＋Δｔでの仮の羽の位置が、式（１）により定まる。
これで、時刻ｔ0よりｔ0＋Δｔに至る羽の動作が決定されたので、この羽の運動について流体・構造連成解析を行なう。これにより、時刻ｔ0＋Δｔにおける節点力ＦとトルクＴが式（５）より算出される。

さらに、これにより、時刻ｔ0＋Δtにおけるｄ²Ｇ／ｄｔ²が式（２）から算出される。

さらに、これにより、胴体９４の位置と姿勢が式（６）により更新される。

以後同様に、式（１）により、時刻ｔ0＋Δｔにおける羽の位置が決定され、流体・構造連成解析により時刻ｔ0＋Δｔにおける節点力がさらに更新され、これより式（２）、式（５）、および式（６）を経て胴体の位置と姿勢が更新される。ここで、更新された胴体の位置と姿勢の変化が、あるスレッショルド以下であるなら収束したと見なし、次タイムステップでの解析に進む。

以上のように、胴体の位置算出と流体構造連成解析を交互に反復収束させることで、重力下で胴体と羽が力学的に拘束された状態での運動の追跡、すなわちfree-flight状態のシミュレーションが可能になる。

以上の概念を図１４に示す。

（モデリング）
以下、図１５〜図１９を用いて、浮上移動装置１およびその周囲流体の、シミュレーションにおける数値表現の作成方法、つまりモデリング手法について説明する。

本実施の形態では、左右アクチュエータの回転中心を結ぶ線の中点を原点とし、左右アクチュエータを結ぶ線を、右アクチュエータから左アクチュエータに向かう方向を正としてｘ軸をとり、後方を正としてｙ軸をとり、鉛直方向上方にｚ軸をとるものとする。

（羽のモデリング）
本実施の形態では、予めＲＯＭ５３に基本となる羽４のメッシュ情報を格納しておく。なお、本実施の形態では、羽全体がメッシュのサイズに比較して薄く、曲げ変形が主要な変形であると考えられるので、計算の合理化のためにシェル要素を用いる。図１５に、本発明者らが用いた羽４のメッシュを示す。これらメッシュを表現するデータを、nodeposition_wing.datおよび、nodelinkage_wing.datという名前のデータファイルに納める。但しこれはデータの認識の便宜のためであり、以下のいずれの説明においても、上記ファイル名も含め、データファイルという概念でデータが存在することは本発明の必須の要件ではない。例えば、ＲＯＭ５３における、上記データセットの先頭のアドレスにて認識するなどの手法も可である。

nodeposition_wing.datは、上記羽４を構成する各節点についての、節点番号及びｘ，ｙ，ｚ座標を記したもので、番号と節点は１対１対応している。すなわち、「節点番号、ｘ座標値、ｙ座標値、ｚ座標値」と表わされるデータが、上記羽４を構成する節点の数だけ並んでいるデータである。

nodelinkage_wing.datは、上記羽４を構成する各三角形（以後、「要素」とも称する）についての、上から見た反時計回りの節点の上記節点番号を並べたものである。すなわち、「要素番号、節点番号１、節点番号２、節点番号３」と表わされるデータが、上記羽４を構成する三角形の数だけ並んでいるデータである。

すなわち、この２つのデータから、羽４の、三角形で近似された空間形状が表現される。

さらに、parameters_wing.datなるファイルには、「要素番号、ヤング率、ポアソン比、厚み」と表わされるデータが、各要素に１対１対応して記される。

これにより、羽４の物理的形状、構成が表現される。

（羽ばたき方のモデリング）
羽４はアクチュエータ２によって駆動される。ホバリング時の羽ばたき位相φと、図６に示される羽の姿勢を表わすα、βは、図１０のように表わされる。但し、０≦φ＜πが打ち上げ、π≦φ＜２πが打ち下ろしで、周期運動であるのでφ＝２πはφ＝０と同一の状態である。この一周期のデータがphi_angle.datとしてＲＯＭ５３に格納されている。また、アクチュエータ２の角速度は、駆動ＰＷＭ制御端子２５２におけるDuty比により比例的に制御される。図１６にDuty比Ｄと羽ばたき周期T、すなわち、φが１度だけ０から２πまで到達するのに経過する時間との関係を示す。これをduty_T.datと称する。

アクチュエータ２は、図示されない駆動方向制御端子２５１と、駆動ＰＷＭ制御端子２５２により制御される。駆動方向制御端子２５１のactive/inactiveにより、αが正/負の方向に回転する。すなわち、制御装置５は駆動方向制御端子２５１を、羽ばたき位相が０≦φ＜π、すなわち打ち上げの際には、activeに、π≦φ＜２π、すなわち打ち下ろしの際にはinactiveに設定し、さらに、Duty比Ｄに相当するＰＷＭ信号を駆動ＰＷＭ制御端子２５２に入力することで、アクチュエータ２の回転速度と回転方向を制御する。

さらに羽４の駆動は図６におけるα、βより、図１７に示される、羽の根元付近の３つの節点の位置に変換されて表される。図１７に示すように、最も羽の根元の点の、アクチュエータ２の回転中心に対する位置ベクトルＰ０および、これより羽の長軸方向に沿う側の節点への位置ベクトルＰ１、同じく羽面内にあるＰ１に直交する方向の節点への位置ベクトルＰ２および、胴体の重心Ｇに対するアクチュエータ２の回転中心位置を表す位置ベクトルＱＧを用いて、各Ｐ０、Ｐ１、Ｐ２にまつわる節点について、重心に対する、胴体座標系における相対位置ベクトルの３つの成分は、

と表わされる。これが、式（１）におけるａとしてALE-FSI-FEMを用いたfree-flight解析に与えられる。これを便宜上、wingposition.datと称する。

（解析空間のモデリング）
周囲に全く障害物のない状態での解析領域を半径Ｒの球体としてとる。本発明者らの実験によれば、障害物の存在が、羽４の流体力に影響を与える距離は概ね２０ｃｍまでであるので、ここでは解析領域を半径２０ｃｍの球体とし、その略中央に浮上移動装置が配されているものとする。この状態で、この解析領域について三角錐メッシュを作成したものが、ＲＯＭ５３に格納されている。形状の把握が容易であるよう、表面のメッシュのみ表示したものを図１８に示す。実際には、解析空間内部には三角錐のメッシュが作成されている。

羽４と同様、解析空間における節点の座標を表わすnodeposition_air.datおよび、各要素を構成する節点のつながりを表わしたnodelinkage_air.datなるデータにより、上記解析空間は表現され、これがＲＯＭ５３に格納されている。

さらに、parameters_air.datなるファイルには、「要素番号、質量密度、粘性係数、体積弾性率」と表わされるデータが、各要素に１対１対応して記される。

（周囲物体の形状および運動反映）
近傍形状計測センサ５２より得られた周囲形状に基づき、解析空間を以下のように決定する。但し、説明の簡便のため、周囲の障害物は、それ単体で位置を変化させるものではないとする。すなわち、周囲の障害物の運動は、浮上移動装置の移動により起きる相対的なものであるとする。もし上記想定外の場合においても、対応する節点にそれに応じた流速を与えれば良い。

再記するが、初期状態で、解析領域は半径２０ｃｍの球体とし、その略中央に浮上移動装置が配されているものとされている。

ここで、上記解析領域最外面における節点と、上記原点を結ぶ線分と、周囲形状表面の交点を求め、当該節点がこの交点に一致するような強制変位を求める。これを全ての最外周面における節点に対して行なう。上記解析領域を弾性体と仮想して、上述の強制変位に従い各節点を移動すれば、周囲形状に一致した形状の解析領域が得られる。また、上記強制変位を与えた節点の速度は、当該節点の流速でもあるので、この流速を境界条件としてALE-FSI-FEMに与えれば、周囲の物体の運動による流れへの影響を導入することができる。 すなわち、本実施の形態における周囲物体の運動の反映は、上記節点における流速データをvelocity_boundary.datとして有することに他ならない。このようにして、新たにメッシュを生成することなく、初期にＲＯＭ５３に格納されたメッシュのみで解析領域の形状変化に対応することができる。

なお、本実施の形態では、ALE-FSI-FEMにおけるメッシュ制御に、仮想弾性体の変形を用いているので、これら解析空間の形状を反映したメッシュの形状算出は独立して行なう必要はなく、各節点における強制変位を与えれば、ALE-FSI-FEMの演算により上記処理は行なわれる。このため、本実施の形態における解析領域形状の反映は、上記強制変位を与える節点の位置を、BCdata.datなるデータファイルとして有することに他ならない。

（胴体のモデリング）
胴体は、先述のように、重心位置を（０，０，−０．００４）とする、表１に示す、物理的パラメータを有する剛体であるとする。故に、当該胴体の運動にかかるパラメータは並進及び回転に関する慣性であり、前者は質量であり、後者は回転慣性モーメントである。これらのパラメータと重心位置を、bodyparams.datというデータファイルに格納して有する。

（羽駆動決定手法）
（制御手法の概要）
先ず、制御手法の概要について、図１９を用いて説明する。

例として、図１９に示すように、幅２０ｃｍの途中で折れる、くの字型の通路を通過する際の制御手法を考える。この移動経路を目標経路として、idealmotion.datと称する。

この際の移動経路およびそれを実現するための羽の駆動方法は複数考えられるが、本実施の形態では説明の簡便のため、後述する基礎的飛行制御を用い、上記通路の略中心を通る経路を実現する羽の駆動方法を、流体・構造連成解析を用いて求める手法について解説する。

（具体的な駆動決定手法の説明）
以下、具体的なデータを元に、図１９〜図２４を用いて、本実施の形態における羽駆動決定手法を例示する。なお、以下の説明では、あるデータ処理のまとまりを一つの構成要素として扱い、図番も独立して与えているが、これは説明の便宜上のことであって、実際にハードウェアとして独立している必要はない。また、プログラムとして独立している必要もない。当然、これらを統合したプログラムを格納した単体のＣＰＵにて全ての処理が行なわれることがメモリやハードウェア節約の観点からは、容易に想定される。

（羽駆動決定プロセスの概要）
本実施の形態では、説明の簡便のために、上述のモデリングの項で説明した、打ち上げと打ち下ろしの速度比率ｄの変更のみにて制御を行なうものとする。すなわち、本実施の形態における羽駆動の算出は、左右各羽の、上記比率ｄの時刻歴を決定することに換言される。

（既知のデータ）
上述のモデリングの項で説明したデータは、変更の必要がないものなので、便宜上、予めＲＯＭ５３に格納されている。すなわち、nodeposition_wing.dat、nodeposition_air.dat、nodelinkage_wing.dat、nodelinkage_air.dat、parameters_wing.dat、parameters_air.dat、parameters_body.dat、phi_angle.dat、duty_T.datが予めＲＯＭ５３に格納されている。

また、上記基本制御による、基本制御パラメータがＲＯＭ５３に格納されている。便宜上これを、motiontable.datと称する。これは、羽駆動演算部５５１が羽駆動の導出演算を行なう初期値として用いるもので、ホバリング状態からの上記基本制御による羽ばたき方変更により生ずる運動態様の変化が格納されている。本実施の形態においては説明の簡便のため、羽ばたき周波数一定である場合の、打ち下ろし・打ち上げ動作における速度の比率ｄに対する、前後方向への、胴体への作用点に与えられる前後方向の力Ｆ２が、後方を正として、図２０のごとく与えられているものとする。

当然、これらは一例示であり、データ格納位置、内容等は、その目的とするアプリケーションにより変化する。また、ＲＡＭはＲＯＭの代用が可能であるので、これらのデータがＲＡＭ５４に格納されている構成も可である。

例えば、羽の一部が飛行中に欠落するようなアプリケーションの場合、nodelinkage_wing.datをＲＡＭ５４に格納しておくことで、上記のような形状変化に対応できる。

（近傍形状把握）
先ず、近傍形状計測センサ５２から得られた近傍形状データが、境界演算部５５３により、境界における座標、およびその運動速度に変換され、この値がBCdata.datとしてＲＡＭ５４に格納される。

（胴体挙動計測）
胴体運動演算部５５４により、６軸加速度センサ５１より計測された６自由度の加速度データが、浮上移動装置１の重心における加速度および角加速度に変換される。また、この値を一階ないし二階積分することにより、浮上移動装置１の重心の速度、角速度、位置、姿勢が求まる。これらの値がbodymotion.datとして、ＲＡＭ５４に格納される。

説明の簡便のため、ホバリング状態からスタートするものとする。すなわちｔ＝０における浮上移動装置１の速度、加速度、角速度、角加速度は全て０である。

（反復初回の羽駆動データ算出）
羽駆動演算部５５１は、上記目標経路に対して、これを実現する上記比率ｄの時刻歴を、上記motiontable.datに基づいて、左右それぞれの羽について求める。以後、これをｄＬ、ｄＲとする。

胴体を剛体と仮定して、motiontable.datのデータより、ｄＬ、ｄＲの時刻歴が、例えば以下の手法により求まる。
胴体を質量Ｍ₀、慣性モーメントＩ₀の左右対称な剛体として、上記胴体重心に対する、胴体と左右の羽の力学的作用点の位置ベクトルをｌＬ、ｌＲ、上記力学的作用点にかかる力をそれぞれＦＬ、ＦＲとすると、重心の加速度ｄ²Ｇ₀／ｄｔ²と、重心を通る方向余弦Ｅの軸周りの角加速度ｄ²θ₀／ｄt²はそれぞれ、
ｄ²Ｇ₀／ｄｔ²＝（ＦＬ＋ＦＲ）／Ｍ、
ｄ²θ₀／ｄｔ²＝（ｌＬ×ＦＬ＋ｌＲ×ＦＲ）・Ｅ／Ｉ₀
で表される。

説明の簡便のためにｘｙ２次元平面内での運動に限定すると、各ベクトルの成分を添え字で表すとして、この面内での加速度は、
ｄ²Ｇx₀／ｄｔ²＝（ＦＬx＋ＦＲx）／Ｍ、
ｄ²Ｇy₀／ｄｔ²＝（ＦＬy＋ＦＲy）／Ｍ
である。

この面内での回転加速度は、
ｄ²θxy₀／ｄｔ²＝（ＦＲy・ｌＲx−ＦＲx・ｌＲy＋ＦＬy・ｌＬx−ＦＬx・ｌＬy)／Ｉ₀
で表される。胴体の姿勢は既知であり、ｌＬ、ｌＲも既知であるので、これはＦＬ、ＦＲについて解くことができ、ＦＬ、ＦＲの値が得られる。

ＦＬ、ＦＲとｄＬ、ｄＲの関係はmotiontable.datにより与えられているので、これにより、ｄＬ、ｄＲの時刻歴が求められる。

図２１に、これにより求まった、右羽と左羽の上記比率ｄの時刻歴、ｄＬ、ｄＲを示す。

以上のように、ＲＯＭ５３に格納されたmotiontable.datより、比率ｄの時刻歴を逆算することができる。なお、本実施の形態に示された比率ｄの時刻歴算出手法はあくまで一例であり、本発明の権利範囲をこれにより制限するものではない。

（Free-flightシミュレーションへの受け渡し）
このようにして求まった比率ｄによって、duty_T.datを用いてある時刻における羽ばたきの位相φが求まり、このφから、phi_angle.datを用いてα、βの値が求まる。すなわち、α、βの時刻歴が、例えば図２１における時刻ｔ1の左羽について、図２２のように求められる。すなわち、ここでは打ち上げが打ち下ろしより速いため、前進方向に力が得られ、左側がより速く前に進むことで右折することがわかる。さらに、式（７）により、羽４のモデルの根元の節点、Ｐ０、Ｐ１、Ｐ２の座標の時刻歴が求まり、各節点の初期位置を差し引けば、上記節点に対する強制変位の時刻歴が求まる。同様に、右羽についても当該データが求まるので、羽の駆動態様が求まる。これを便宜上、wingposition.datと称する。

但し、浮上移動装置１の移動速度が大きい場合などのように、解析領域の範囲を超えて移動してしまうことや、実際に求めたい、移動方向の空間に対して広い解析領域をとれない等の問題が生じる場合は、ある速度で移動している観測枠から観測した状態でのシミュレーションを行なうことにより、この移動を相殺することが望ましい。つまり、解析空間全体にオフセットを与える。例えば浮上移動装置１が、１０ｍ／ｓで前進している場合、これに並行に１０ｍ／ｓで移動する観測枠から観測した解析領域を定めれば、常に浮上移動装置１の解析領域に対する位置は変化しないため、最大限に解析領域を活用した解析が行なえる。この場合、観測枠の運動ベクトルのマイナス１倍を、流速および胴体の運動に与えれば矛盾のない解析が行なえる。

当然、このオフセット量は、胴体の運動そのものであることが、解析領域を最大限活用できて、上記の観点からは有利である。

（Free-flightシミュレーション）
上記ＲＯＭ５３に格納されたデータにより浮上移動装置及びその周囲流体のモデルが与えられ、さらに、上記反復初回の羽駆動データとして算出されたwingposition.datより羽の運動が与えられており、また、ＲＡＭ５４に格納されたBCdata.datより境界の形状と移動速度が与えられているので、現時点での胴体の運動態様bodymotion.datを初期値とするfree-flightシミュレーションが可能になる。

これらのデータより、free-flightシミュレーション演算部５５２は羽ばたき浮上装置１の運動態様を予め算出し、predictmotion.datとして、羽駆動演算部５５１に渡す。

（２回目以降の反復）
初回反復は周囲障害物などの地面効果等を考慮せずに決定された羽駆動手法であるので、当然、図２３に示されるように、predictmotion.datと、目標経路idealmotion.datとは異なる。

このpredictmotion.datと、idealmotion.datの差異から、羽駆動演算部５５１は、上記羽ばたきの比率ｄＬ、ｄＲの修正を行なう。ここでは、図２４に示すように、同一時刻における、進行方向の差異を表わすθideal、θpreを用いて比例的に、
ｄ＝１＋（ｄ−１）×θideal／θpre
として修正を行なうことで、より目標経路に近い経路が実現できる駆動方法が算出される。

なお、本来これは一種のパスプランニング問題であるが、本発明における目的は、パスプランニング手法の開発ではないため、ここに挙げた手法は一例示であると考えるべきであり、他の効率の良いパスプランニング手法を採用することは、本発明の権利範囲を損ねるものではない。例えばＲＡＭ５４に格納されている、前回の反復によるpredictmotion.datをＲＡＭ５４に格納しておき、ｄＬ、ｄＲの変更による飛行態様の変化の感度を求め、これによりｄＬ、ｄＲの変更を行なう手法も可である。これによれば、計算はやや複雑になるが、理想状態での計測データを用いる前者の手法よりも、実際の影響を反映した正確な羽駆動算出が行なえる等のメリットがある。

このようにして求まった新たなｄＬ、ｄＲを用いて、新たなwingposition.datを算出し、改めて、free-flightシミュレーション演算部５５２により、predictmotion.datを求める。

このようにして求まった、新たなpredictmotion.datと、idealmotion.datを比較し、その差異が、許容範囲より大きい場合は反復計算を継続する。逆にこれが、許容範囲以内であるならば、反復計算を打ちきり、この時点でのｄＬ、ｄＲを、羽駆動手法として採用する。これをduty_T.datにより、Duty比の時刻歴に変換し、左右アクチュエータに出力することで、流体・構造連成解析によるfree-flightシミュレーションを用いた制御が実現する。

（データフロー）
図２５に本実施の形態における羽駆動決定のデータフローを示す。

なお、説明の簡便のため、本実施の形態ではスタンドアロンタイプを想定したが、流体・構造連成計算等における計算量が大きく、羽ばたき浮上移動装置に搭載できる回路規模では現実的な計算時間で計算を完了できないなどの場合は、当該計算部やその他の処理部分を、浮上移動装置外に独立して設け、通信機能を浮上移動装置に付加し、これと通信を行なうことで同様の機能を実現する手法も当然可である。

（フローチャート）
図２６に本実施の形態における羽駆動決定のフローチャートを示す。このフローチャートは、本発明の他の局面の浮上移動装置の基本概念を表わすものである。

（浮上可能要件）
本発明者らの実験によれば、図２４に示す羽４の運動により発生する浮上力の最大値は、羽１枚あたり約０．１３ｇｆである。また、この浮上力を得る際に必要なアクチュエータ２の駆動トルクは最大約１ｇｆ・ｃｍである。

羽４の質量は約５ｍｇである。主軸３３の質量は約３ｍｇである。突起３４の質量は約２ｍｇである。ストッパ３５１，３５２の質量は約６ｍｇである。アクチュエータ２の質量は約８０ｍｇである。

支持構造９、制御システム５および電源６の重量の合計が６８ｍｇ以内となるように構成すれば、浮上移動装置１は浮上が可能である。但し、電源６を無線で供給し、無線送信される電源の変換部分および制御システム５をワンチップに集積して、そのワンチップを支持構造９上にパッケージングすれば、支持構造９、制御システム５の合計が５〜１０ｍｇ以内となるように構成することができる。

（駆動エネルギ供給方法）
本実施の形態の浮上移動装置１においては、駆動エネルギ供給方法は、浮上機能を損なわない限りその手法に制約はない。

たとえば、浮上移動装置１内に充電池または燃料電池等を内蔵する手法、または、外部のエネルギ供給源から電波を用いて、羽４に設けられたアンテナに電力を送電する手法などが考えられる。

（その他補足）
上述したストッパ、突起、およびストッパ移動機構は一例であり、浮上移動に支障無く、本実施の形態に示す羽ばたき方を実現することができるのであれば、浮上移動装置の駆動機構は上述した機構以外のものであってもよい。

本実施の形態の浮上移動装置１においては、現段階で浮上移動を実現することが可能である機構を提示するために、アクチュエータ２として超音波モータを用いたが、浮上移動装置１が浮上移動できるのであれば、超音波モータ以外の駆動源を用いても良い。また、アクチュエータ２の運動は、最終的に伝達軸３３に前述した運動を行なわせるものであるならば特に限定が必要なものではない。たとえば、軽量化を狙って高分子材料を用いたリニアアクチュエータ等を使用し、リンク機構によって伝達軸３３を運動させる手法が用いられてもよい。また、２サイクルエンジンのような往復運動が容易な内燃機関を用いる手法も考えられる。

また、伝達軸３３および羽４の運動も、本実施の形態に示したものに限らない。例えば伝達軸３３を完全な剛体とすることは難しく、伝達軸３３がしなる運動をすることも考えられ、この場合は、伝達軸３３の先端が横に８の字を描く運動をすることが考えられるが、極端に柔らかく、浮上力を損なうような伝達軸３３で無い限りこれは特に本実施の形態の必要要件を逸脱するものではない。またさらに、アクチュエータ２の回転半径を無限に大きいと考えると、伝達軸３３は、単純に往復運動を行なうが、この場合も本実施の形態の要件は満足される。

羽４の形状は、浮上移動装置が効率的に浮上移動可能なものの一例であり、浮上移動を実現することが可能であるならば、羽４の構成または材料などは、他の構成または材料であってもよい。

また、回転角制限機構については、その制限範囲は固定である必要はない。図５に示されるように、空気力によって羽の変形は異なるので、最適な迎え角を維持するためにこれを制御することは浮上の効率化に繋がる。

以下、本発明において用いる流体構造連成解析の一例を説明する。本発明に用いられている流体構造連成解析の手法および用語の定義は以下の手法と同様のものが用いられている。たとえば、羽モデルまたは構造モデルという用語は、以下の構造に関する数値モデルを示すものであり、流体モデルという用語は、以下の流体に関する数値モデルである。

（実施の形態１）
まず、表２〜表７および図２７〜図４６を用いて、実施の形態１の流体・構造連成数値モデルの作成方法およびそれを用いた羽ばたき飛行ロボットの製造方法を説明する。

本実施の形態の流体・構造連成数値モデルの作成方法は、昆虫が空気中で羽ばたき飛行する場合の、昆虫の羽根の構造および昆虫の羽ばたき飛行態様を解析して得られる、流体としての空気に関する数値モデルおよび羽根の構造に関する数値モデルを作成するためのものである。

また、本実施の形態の羽ばたき飛行ロボットの製造方法は、前述の流体・構造連成数値モデルの作成方法により作成された数値モデルを用いて、昆虫の羽根の構造および昆虫の羽ばたき飛行態様が模倣された羽ばたき飛行ロボットを製造するためのものである。

具体的には、流体に関する数値モデルとは、流体に関する流速および圧力の数値モデルのことを指す。また、構造に関する数値モデルとは、主に、流体に接する構造物の移動および変形などの運動の態様の数値モデル、ならびに、構造物に作用する内部応力などの力の数値モデルを指す。

本実施の形態の流体・構造連成数値モデルの作成方法においては、実際のトンボの飛翔に関する数値モデル化について説明する。より具体的に言うと、飛翔するトンボの羽根の構造に関する数値モデルおよび羽根の周囲の空気に関する数値モデルを作成する手法について説明する。すなわち、羽ばたくトンボの羽根の周囲の空気の流速および圧力の数値モデルと、トンボの羽根の構造が周囲の空気から受ける圧力、トンボの羽根を駆動する外力（駆動力）、トンボの羽根構造の垂直応力（圧力または張力）および剪断力などの力、トンボの羽の並進移動または回転移動、ならびに、トンボの羽根の変形などの数値モデルとを、実物のトンボの空気中での羽ばたき飛行において計測された計測データに流体・構造連成解析を用いることによって求める手法を述べる。

なお、本実施の形態においては、力学的要素についてのみ数値モデル作成を行なうが、筋肉の駆動モデル、筋肉を駆動する神経系のモデル、脳における情報処理のモデル、および、情報処理のもとになる感覚器官におけるセンシングのモデルなどの数値モデルについても、同様の手法を適用することが可能である。

（流体・構造連成数値モデルの作成の手順の概要）
まず、流体・構造連成数値モデルの作成の手順の概要について図２６を用いて説明する。

本実施の形態における流体・構造連成数値モデル（単に「数値モデル」という場合もある。）の作成の手順は、１．基準構造の計測、２．基準構造のモデリング、３．実構造物の計測、４．基準構造の補間による、実構造等価数値数値モデルの作成、５．羽ばたき代表点の位置の時系列データの計測、６．本体（羽以外の構造部分）のモデリング、７．流体・構造連成解析の各プロセスに分けられる。

プロセス１〜プロセス４により、羽根の形状と剛性とを数値で表現する、羽根の実構造等価モデルの数値モデル（以後、「実構造等価数値モデル」と称する。）を作成し、プロセス５により、羽根の羽ばたき飛行態様の数値モデル、すなわち、羽根の駆動力モデルを作成する（以後、「羽ばたき運動数値モデル」と称する。）。

なお、プロセス６の本体のモデルは、説明の簡便のため、形状および質量分布のみが数値化されているものとする。つまり、本体のモデルは、流体に対する境界条件を与える役割、および、羽根の支点において並進移動や回転のための慣性を与える役割のみに用いられているものとする。

但し、実際には、本体の数値モデルも、その位置、姿勢および形状を変化させることが考えられるが、たとえ、これらを考慮する場合であっても、位置、姿勢および形状の変化についても、数値モデルとしては、羽根の流体・構造連成数値モデルの作成方法と同様の方法を適用して作成することができる。

より具体的な数値モデルの作成の手順は、次に述べるようなものである。

まず、必要に応じて、あるトンボ（以後、「個体Ａ」と称する。）の羽根を生体から切離すなどして、羽根の構造に関する物理量を精密に測定する。それにより、その測定された物理量を用いて基準構造数値モデルを作成する。また、他の個体（以後、「個体Ｂ」と称する。）の羽ばたき飛行態様の物理量を測定して、その測定された物理量を用いて羽ばたき運動数値モデルを作成する。

次に、固体Ｂに対しては、個体Ｂを損傷しない状態で計測可能な羽根の部位のうちの代表点の構造に関する物理量を測定する。この代表点の物理量のデータを用いて、個体Ａの構造に関する数値モデル、すなわち、基準構造数値モデルに何らかの変換を加えることで、個体Ｂの羽根の構造の数値モデルと等価とみなすことができる数値モデル、すなわち、実構造等価数値モデルを作成する。

そして、個体Ｂの羽の実構造等価数値モデルを、個体Ｂより直接測定した羽ばたき方によって駆動する際の流体の挙動および構造の挙動を、流体・構造連成解析によって算出する。これにより、個体Ｂについて、周囲の流体からの影響を含めた羽ばたき運動時の流体・構造連成数値モデルが作成される。

次に、前述の手法により得られた流体・構造連成数値モデルをベースとして、羽ばたき飛行ロボットを製造する手順について述べる。

以下、各手順について、図２７〜図４０を用いて詳細に説明する。

（基準構造数値モデルの構造に関する物理量の計測）
まず、基準構造数値モデルを作成するための、個体Ａの構造に関する物理量の計測を行なう。一般に、構造物の運動方程式は、バネとダンパ、つまり、変位に対する弾性力と速度に対する減衰比とを用いて、外力と加速度との式で表わされる。

ところで、一般に、構造における内部応力の減衰は、構造内部において運動エネルギが熱エネルギへ変換されることにより行なわれる場合に起きるが、これは構造の変化または破壊が起きていることに相当する。ところが、たとえば、羽ばたき周波数３０Ｈｚのトンボが羽根を自律的要因で破壊することなく１週間以上飛び続けることを考えると、数周期での羽ばたき運動では、上記羽根の構造の変化（塑性変形）は極めて小さいものと考えられる。よって、実施の形態においては、構造における内部応力の減衰はゼロであるものとする。すなわち、羽根に関しては、弾性変形はするが、塑性変形はしないものとして取り扱うこととする。

これにより、構造における運動方程式は、変位に対する弾性、すなわち、羽根の剛性、羽根の質量、および、外力によって表わされることになる。

すなわち、羽の構造モデル作成に必要なパラメータは、羽の形状、剛性および質量の３つである。以後、これらの計測方法について図２８〜図４３を用いて述べる。

（基準構造数値モデルの形状計測）
まず、羽根の形状計測に関して、図２８および図２９を用いて説明する。

図２８に示すように、本実施の形態にて用いるトンボの羽根２は、枝状に分布する梁構造２１に膜構造２２が貼られている構成となっている。後に述べる剛性モデリングのために、この梁構造２１と膜構造２２との位置関係を把握しておく必要がある。

このために、本発明者らが形状計測に用いたのは以下の手法である。

まず、３次元での羽根の形状の把握のため、一般的に市販されているＸ−Ｙステージ３３と、このＸ−Ｙステージ３３の主表面にほぼ垂直な方向の距離を計測可能な状態に配されたレーザ距離計３２とを用いる。

図２９に示すように、レーザ距離計３２をＸ−Ｙステージ３３の主表面内のある位置に定位させた状態で、羽根２までの距離を、レーザ距離計３２によって計測する。これにより、羽根２全体の形状が把握される。すなわち、Ｘ−Ｙステージ３３の主表面上のレーザ距離計３２の位置を特定可能な位置指定パラメータｘ（ｉ），ｙ（ｉ）と、この際のレーザ距離計３２の計測値ｚ（ｉ）とを用いて、パラメータ（ｘ（ｉ），ｙ（ｉ），ｚ（ｉ））の集合で表わされる形状が、測定された羽根２の形状である。

また、図３に示すように、画像スキャナまたはデジタルカメラなどの画像取得手段により、Ｘ−Ｙステージ３３の主表面にほぼ平行な面内における羽根２の画像を取得し、これにより梁構造２１のＸ，Ｙ方向の位置を得ることができる。羽根２をＸ−Ｙステージ３３の主表面のほぼ平行な位置に配置すれば、羽根２における梁構造２１は、Ｚ方向に重複して存在することはないので、梁構造２１は前述の画像より計測されたｘ（ｉ），ｙ（ｉ）に対応するレーザ距離計３２により計測されたＺ方向の距離、すなわち、計測値ｚ（ｉ）の位置にのみ存在することになる。

これによって、すべての梁構造２１のｘ、ｙ、ｚを座標軸とするときの配置が決定される。なお、羽根２の形状を把握する手法は、３次元で梁構造２１の配置を含めて羽根２の形状を把握することができる手法であれば他の手法であってもよい。

たとえば、市販のレーザ距離計３２の多くは、レーザ光が十分に戻ってこない場合に、計測精度が低下する。このレーザ光が十分に戻ってこない状態を検出するために、市販のレーザ距離計３２の多くには、レーザの戻り光量をモニタするためのモニタ機能が付加されており、このモニタ機能によって、レーザ距離計３２の測定時における反射率を算出することができる。なお、トンボの羽の場合、膜構造２２はほとんど透明であるので、散乱光の大小を、羽根２の位置に対してマッピングすることによって、羽根２の梁構造２１を検出する手法を用いてもよい。

（基準構造数値モデルの剛性計測）
羽根の剛性計測は、基本的には一般的な手法がそのまま適用可能である。より詳細言うと、次のような手法である。まず、前述したように、羽根２は梁構造２１と膜構造２２との複合構造であるので、羽根２の各代表的な部位については、梁構造２１と膜構造２２とを分離する。この分離された梁構造２１および膜構造２２それぞれについて引張り強度試験、曲げ強度試験などを行なうことにより、剛性を決定する。さらに、膜構造２２の一部に一定形状の孔を開ける。その孔形状の変形から、膜構造２２にかかる張力を逆算する。

しかしながら、この手法は煩雑であるという問題がある。また、この手法では、梁構造と膜構造の相互作用の結果、曲げ剛性が決定されるので、両者の誤差を含んでしまうという問題もある。

ところで、羽根の変形は、曲げによる変形と、引張りによる変形とに分けられるが、昆虫の羽根は非常に薄いため、羽根の変形を大きく支配するものは曲げによる変形である。すなわち、昆虫の羽根は、歪が微小でかつ変形が大きい。

そこで、本実施の形態においては、羽根２をシェル構造の集合体として近似する。すなわち、１つの膜構造２２および梁構造２１の複合要素を１つのシェル構造として取扱う。これにより、羽根の剛性については、羽根全体をシェル構造の複合体と考えて、梁構造２１と膜構造２２とを複合した状態のまま測定すればよい。このように、羽根全体をシェル構造の複合体と考えて、羽根および流体の挙動を主に決定する羽根の曲げ剛性を、昆虫の羽根から直接測定することで、実質上、実物の昆虫の羽根の曲げ剛性に近い値の曲げ剛性を有する羽根の基準構造数値モデルを作成することができる。また、梁構造と膜構造との複合要素に比べて梁構造を除外した分だけ、基準構造数値モデル作成のための計算量を低減することができ、作業の効率化を図ることができる。

なお、羽根２の各部位について、個別に計測を行なうため、羽根２が損壊（分割）した状態で測定せざるを得なくなることが十分考えられる。そのため、後述する実構造等価数値モデル作成の観点から、これらの剛性測定は、羽根２が均一に乾燥した状態で行なわれることが望ましい。

（基準構造数値モデルの質量分布の計測）
剛性の計測の項で説明したように、羽根２の構造を梁構造２１と膜構造２２とに分けるのであれば、分けた状態での羽根の代表的な部位の梁構造２１と膜構造２２との各々の質量を計測すればよい。しかしながら、羽根２をシェル構造として取扱うのであれば、羽根２の各部位の梁構造２１および膜構造２２の複合構造の質量を計測すればよい。

また、実際には、羽根２は生体より切離されると乾燥を始め、質量が軽減されてしまうので、まず、生体より切離した時点で羽根２全体の通常状態の質量を計測しておく。その後、羽根２を十分乾燥させた後に再び羽根２全体の乾燥状態の質量を計測する。乾燥した羽根２を測定するための各部位ごとに分割して、分割された各々の部位について質量を測定する。そして、通常状態の羽根２全体の質量に対する乾燥状態の羽根２全体の質量の比で、先の乾燥状態での分割された羽根２の各部位の計測された質量を除算することで、通常状態での羽根２の分割された各部位の質量に変換することが望ましい。

これらより得られた質量を、質量を測定した断片のサイズにより除算しておき、その除算された値を質量分布とする。たとえば、シェル構造の場合、断片の面積で断片の質量を除算することにより、単位面積当りの質量が求まる。

以上により、羽根２の形状、および、羽根２の代表部位での剛性と質量とが測定される。

（基準構造のモデリング）
続いて、基準構造となる個体Ａの羽根２の構造のモデリングについて、図２８〜図３０を用いて説明する。

一般に、構造の解析に多用される手法は、羽根２を単位構造の集合として分割し、かつ、分割された単位構造に対して種々の物性値を付与する手法である。この単位構造を一般的にメッシュと呼ぶ。

本実施の形態では、羽根の単位構造をシェル構造として取扱う。すなわち、羽根２は、上記の形状、質量分布、および、剛性からなる物性値を有するシェル構造の集合体によって表現される。

（形状のモデリング）
形状に関しては、各メッシュを構成するノード（メッシュを示す線の交点）の位置を付与する。また、解析手法によっては、各ノードの有する姿勢（複数の点位置データ）を付与する。

この際、メッシュ分割は、梁構造２１の方向を基準に行なわれることが望ましい。特に、コルゲーションと呼ばれる凹凸が梁構造２１を稜線として走っており、このコルゲーションを境として、羽根２の物性値も極端に変化するので、１つのメッシュは、この梁構造２１を横切らないことが望ましい（これが梁構造２１の配置の把握が重要な理由である）。これらの点以外は、メッシュ作成方法は従来から用いられてきた手法と比べ、特に変更を要する部分はない。

（剛性のモデリング）
次に、各メッシュに剛性を付与する。

あるメッシュにおける剛性、すなわち、ある外力に対する変形を決定する基本的なパラメータは、ヤング率、ポアソン比、およびメッシュの厚みである。

前述のように本実施の形態における羽は、歪みが微小でかつ変形が大きいものであるので、これらシェル構造の変形は、ヤング率と断面２次モーメントとの積によりほぼ決定される。なお、ポアソン比は変形にほとんど影響を与えないので、一般的な値とされている０．３を用いる。また、後述する実構造等価数値モデル作成の観点から、ヤング率は、全代表位置での計測の結果得られた値の平均など、ある程度ラフな値を用いて構わない。そして、各代表位置において計測された曲げ剛性を用いて、理論解または数値計算結果からメッシュの厚みを逆算する。

たとえば、長さｌ、高さｈ、幅ｂの梁の一端を固定し、他端に高さ方向に荷重ｗを与えた際の高さ方向の梁の変位ｘは、ｘがｌに比べて微小であるなら、ｘ＝ｗ×ｌ³／（Ｅ×ｂ×ｈ³／１２）と表わされ、シェルの幅ｂは測定により求まるため、ヤング率Ｅが既知であれば、梁の高さｈを算出することができる。ここでは、Ｅ＝１．０×１０⁹とする。

これにより、各代表位置でのメッシュの厚みが決定されるので、これを補間した値を各メッシュに付与する。なお、各代表位置での厚みは、梁構造２１と膜構造２２との複合構造として曲げ剛性を表現するための厚みであり、梁構造２１そのものの厚みおよび膜構造２２そのものの厚みのいずれとも異なる。例として、ＭＩＴＣ（Mixed Interpolation of Tensolial Components）４節点シェル要素を用いたモデルの厚み分布を図３０に示す。

（質量のモデリング）
質量は、単純に、各代表位置での質量計測結果から算出される単位面積当りの質量を補間して求まる。メッシュの位置での単位面積当りの質量に、メッシュの面積を乗じた値を付与すればよい。

以上により、羽根２の形状、剛性、および、質量分布を表現した構造に関する物理量の数値モデルが作成される。以後これを、基準構造数値モデルと称する。

（実構造等価数値モデルの作成）
ここでは、図３１〜図３３を用いて、生体である個体Ｂを損壊せずに固体Ｂの代表的な部位の形状と代表的な部位の剛性とを計測する手法について説明する。

（代表形状の計測）
代表形状を計測する場合は、まず、羽根２０を撮影した画像などから、代表的な構造、たとえば、羽根２０の外形線や、梁構造２１の特徴的な分岐点からの距離を計測する。これにより、基準構造数値モデルに対する形状の比率α、すなわち、ある共通の基準点から測定点までの羽根２０の基準構造数値モデルに対する距離の比率を算出する。たとえば、外形線を計測し、羽根２０の長手方向の長さが、基準構造数値モデルｌ₁、ここで計測した個体Ｂすなわち羽根２０の外形線がｌ₂であるなら、形状の比率α＝ｌ₂／ｌ₁となる。

（代表剛性の計測）
ここでは、生体である固体Ｂを破壊せずに、羽根２０の代表的な部位の剛性を計測する手法を説明する。

一般に、剛性測定は、ある物体の一端を固定し、他の部位にある荷重をかけ、これに対する変位量を計測することで行なわれる。このため、上述の固定は、他の部位にかけられる荷重によって変化しないだけの強固なものであることが望ましい。ところが、生体を生理的に損傷することなく測定するには、荷重が羽根２０に与えるダメージが小さいことが要求される。

そこで、生体の羽根２０を損壊せずに、代表的な剛性を正確に測定するためには、以下の２通りの手法が考えられる。

１つは、羽根２０の凹凸形状に合わせた固定手段を用いる手法である。しかしながら、昆虫の羽の個別の形状の差異によっては、それに応じた固定手段を用意する必要があるため、固定手段を用意する手間が煩雑である。

そこで、固定による形状の差異の少ない部分に対して上記手法を適用することが望ましい。たとえば、図３１に示す羽根２０の先端後部はほぼ平面であるので、この領域を図３２に示す如く線上に固定し、他の部位に荷重をかけ、変位を測定する。

または、丈夫な羽根２０の梁構造２１部分、特に、羽根２０の梁構造２１の根元部分を固定する手法も有効である。たとえば、図３３に示すように、梁構造２１の根元を固定しておき、羽根２０の先端部分において荷重変位関係を計測する。

（基準構造数値モデルの変換による実構造等価数値モデルの作成）
ここでは、基準構造数値モデルを、固体Ｂすなわち羽根２０の生体の代表値計測によって得られた代表的な測定値を利用して変換し、実構造等価数値モデルを作成する手法について説明する。

まず、代表形状の計測によって求まった形状の比率αによって、基準構造数値モデルの拡大縮小を行なう。この拡大縮小されたモデルを中間モデルと称する。

すなわち、各メッシュのノードの座標とメッシュの厚みとをそれぞれＰ、Ｔ、Ｐ、Ｔそれぞれを代表形状によって求まった形状の比率αによって変換された後のノードの座標をＰ′、メッシュの厚みをＴ′とすると、Ｐ′＝Ｐ×αおよびＴ′＝Ｔ×αの関係が成立する。中間モデルにおいては、形状と質量が個体Ｂのそれをほぼ反映したものとなっている。

続いて、この中間モデルに、上記個体Ｂの代表剛性の計測に用いた条件と同様の力学的条件、すなわち、固定条件と荷重条件とを数値モデル上で与え、数値解析により構造の変形量を求める。そして、生体の代表値剛性と同等の荷重における変位がβであるならば、各メッシュのヤング率Ｅを、Ｅ′＝Ｅ×βによって求まるＥ′に変更する。

すなわち、中間モデルにおける、ある力学的条件に対する変位の値が、生体より計測された同条件での変位の値の０．５倍であるならば、中間モデルのヤング率を０．５倍にすることで、生体と同等の荷重変位関係を有するように中間モデルを変換する。

これにより、中間モデルのヤング率を変換したモデルは、形状、質量および剛性それぞれが、個体Ｂにおける羽根２０とほぼ一致する。これ以後、実構造等価数値モデルと称する。以上のプロセスにより、実構造等価数値モデルが作成される。

（羽ばたき代表点位置計測）
続いて、個体Ｂにおける羽ばたき代表点の位置計測により、羽根２０を駆動する姿勢（位置データの数値モデル）、すなわち、羽ばたき方のモデルを作成する手法について説明する。羽根２０の姿勢（位置データの数値モデル）は、第一義的にはアクチュエータの駆動態様そのもの、すなわち、羽根の支点の位置と姿勢である。ところが、羽根の支点は固体Ｂの体内にあり、実質上測定は不可能である。また、実質上有用であると思われる、高速度カメラ等の映像からの計測では、点の位置は計測できても姿勢は計測できない。

そこで、羽根の支点と実質上等しい位置と姿勢を有する、高速度カメラにて撮影可能な３点の位置を計測し、この３点がなす平面の姿勢をこの３点それぞれの姿勢とみなすことで、羽根２０の姿勢を近似的に決定する。すなわち、羽根２０の支点を含む、変形量の微小な領域は、羽根２０の支点と同じ姿勢である。つまり、この領域内における３点のなす平面の姿勢は、羽根２０の支点の姿勢を表わしている。以上より、羽根２０の姿勢は、この領域内に含まれる３点の位置データで表わすことができる。

より具体的には、羽根２０の上のなるべく変形の小さい位置に３つのポイントマーカを配し、高速度カメラなどで撮影した羽根２０の画像からポイントマーカの位置を計測する手法が考えられる。姿勢の決定の精度から、３つのポイントマーカによって形成される角度は９０°に近いことが望ましい。また、より望ましくは、３つのポイントマーカが直角二等辺三角形をしていることが望ましい。これらの位置計測を、個体Ｂの飛翔を高速度撮影した画像から、ポイントマーカの位置を計測することによって行なう。

この点の位置計測手法としては、従来から用いられている羽根の所定の点の位置計測手法を用いる。本実施の形態の流体・構造連成数値モデルの作成方法においては、上記ポイントマーカの位置を、２方向から高速度カメラにより撮影した画像上の位置を用いて算出した。但し、高速度カメラの画像は、量子化されている。したがって、本実施の形態の流体・構造連成数値モデルの作成方法においては、量子化誤差に起因した所定の点の位置ずれを緩和するために、スムージングを行なっている。流体から及ぼされる力は速度の関数であるため、スムージングは、上記ポイントマーカの位置の時刻歴を時刻で微分した、速度の時刻歴について行なうことが望ましい。なお、スムージングが実行可能であるためには、前述の速度の時刻歴は、連続（連続関数）である必要がある。また、直接直交座標計での位置の時刻歴の時間微分値のスムージングを行なったのでは、各ポイントマーカ間の距離が変化してしまう可能性がある。そのため、このスムージングは、３点により構成される平面の位置と姿勢とを分離した状態で、すなわち、この３点により構成される平面の並進速度と角加速度について行なわれることがより望ましい。なお、点の位置計測手法は、羽根２０のある３点の位置を決定することができる手法、すなわち、羽根の所定の領域の姿勢を決定することができる手法であれば、前述の手法に限定されるものではない。

以下、本実施の形態の流体・構造連成解析で用いる実構造等価数値モデルを、図４１〜４５を用いて具体的に示す。

本実施の形態においては、羽の根元の要素についてこの３点の位置を与える。図４１は、図３０における羽根の根元の要素、つまり左下角の要素のみ抽出したものに、本実施の形態での解析で用いたノード番号と要素番号を付記したものである。ここでは、ホバリング状態を想定し、ノード７１は原点すなわち（０，０，０）に固定とする。そして、この羽根の姿勢は、図４２に示すように、残りの２点、すなわちノード８３とノード３３７の位置すなわちｘ、ｙ、ｚ座標の時刻歴によって数値モデル化される。

筆者らが、前述の東らの文献より算出した羽根の挙動を、図４２に示した座標系において、羽根の打ち下ろしについて時刻と共に示したのが図４３である。また、図４４にノード８３の、図４５にノード３３７のｘ、ｙ、ｚ座標値の時刻歴を示す。

（流体・構造連成解析）
ここで、流体・構造連成解析によって、数値モデルを作成する手法について、図３０〜図３７、図４１、図４６を用いて説明する。

本発明者らが用いた解析手法は、張 群（１９９９年度東京大学学位論文「構造座屈と領域変形を行なう構造・流体連成問題のＡＬＥ（Arbitrary Lagrangian-Eulerian Method）有限要素解析」）により示されたＡＬＥ有限要素法解析による、流体・構造の強連成解析方法である。以下にその適用手法を述べる。なお、ＡＬＥとは、参照座標系を用いてオイラー表記（流体）とラグランジェ表記（構造）を統一的に扱う手法である。また、有限要素法は、解析領域を有限要素に分割して、その要素内で近似を用いて統合した方程式を解く手法のことである。この手法は、差分法と異なり、自由形状が扱え、数学的に収束性が証明されている。

まず、数値モデルを作成するための仮想空間を設定する必要がある。本発明者らの計算によると、羽根の長さ４ｃｍ、羽ばたき周波数３０Ｈｚである場合、概ね２０ｃｍ以遠において空気の流れはほぼ定常流となる。このため、数値モデルを作成するための構造物を含む半径２０ｃｍの球体内の空間を、流体・構造連成解析の対象とすることで、計算量を減らすことができる。

ここでは説明の簡便のため、１辺２０ｃｍの立方体のケースのほぼ中央のトンボが位置し、停空飛翔すなわちホバリングを行なっている際の流体・構造連成数値モデルを作成するものとする。

ＡＬＥ有限要素法においては、次のＡ〜Ｄの４つパラメータを流体および構造の両者において設定する必要がある。なお、以下、Ａは各節点の座標、Ｂは各節点のコネクティビティ、Ｃは各節点の境界条件、Ｄは各要素の物性値とする。

（構造要素）
Ａ、Ｂ、Ｄは実構造等価数値モデルを作成した時点で決定される。

Ｃについては、羽ばたき運動数値モデルとして計測された羽根２０に記したマーカの時系列的に表現された位置データを、実構造等価数値モデルのメッシュのノードの強制変位として与える。

（流体要素）
流体要素はメッシュが用意されていないので、たとえば、市販のメッシュ作成ソフトなどを用いて作成する。本発明者らが用いた流体・構造連成解析手法においては、流体のメッシュにおいても構造のメッシュと共通のノードを有する必要があるため、まず、羽根を含む平面のメッシュを四辺形要素にて作成し、これを上下に掃引（スイープ）することによって六面体メッシュを作成する。また、羽根以外の四辺形メッシュについては、六面体メッシュの作成後、削除した。これを図３４および図３５に示す。

これにより、流体要素についても各節点の座標、コネクティビティが作成される。また、境界条件として、羽根のノードと、立方体のケースの外壁を構成するノードについては、固着境界条件を付与するとともに、空気の質量密度、粘性係数および体積弾性率を、各流体要素の物性値として与えた。

そして、定常的にホバリングしている状態に収束した運動の態様に関する数値モデルを得るために、１周期の羽ばたきを、複数回反復して行なう状態を解析することにした。

［実際のデータ］
解析で実際に用いたデータを表２〜表７を用いて説明する。

各節点の座標については、流体・構造両者に共通である。よって、これは、表２に示す通り、各ノードに番号を振ってNode1、Node2・・・とし、各ノードのｘ、ｙ、ｚ座標と共に列記したものとなる。これを表２では、NodeCoords.datと名付けている。

続いて、羽根の構造を指定するため、羽根を構成するノードのコネクティビティを指定する必要がある。これは、図４１に示される要素１について、左下角の点から反時計回りに、Node71、Node83、Node149、Node337、と指定することで、この４節点シェル要素を表現する。すなわち、表３に示すとおり、羽根を構成する各シェルに番号を振って、ShellElement1、ShellElement2、・・・とし、各シェルを構成するノード番号を、先に説明した順に列挙したものを列記したものとなる。これを表３ではShellMesh.datと名付けている。

同様に、流体領域に作成したメッシュについても、コネクティビティを指定する。但し、流体領域は６面体要素で表わされている。さらに、圧力というスカラー量を表現するため、各流体のメッシュについてさらに１点、圧力節点と呼ばれる節点を付加する。コネクティビティを指定する順は、６面体要素下面の４点を反時計回りに、続いて上面の４点を同じく反時計回りに、最後に圧力節点を、それぞれノード番号を列記することによって、合計９点のコネクティビティの集合が、表４のように表わされる。表４においてはこれを、FluidMesh.datと名付けている。

また、流体、構造それぞれに物性値を付与する必要がある。表５は構造の物性値を付与するものであり、表６は流体の物性値を付与するものであり、それぞれ、ShellMaterial.dat、FluidMaterial.datと名付けられている。ShellMaterial.datにおいては、羽の各構造要素について、ヤング率、ポアソン比、比重、厚さが列記されている。また、FluidMaterial.datにおいては、各流体要素について、粘性、比重、体積弾性率が列記されている。

また、表７には、ShellMotion.datにおいて、羽根の節点におけるｘ、ｙ、ｚ座標の時刻歴が表わされている。

なお、以上のデータは一例であり、データ形式、およびその数値はこれに限るものではない。

［計算結果］
以下、上記の手法により得られる数値モデルについての一例を示す。なお、この数値モデルは上記データにより得られたものであり、この数値モデルに限るものではない。

図３６に、定常状態における羽根２０の周囲での流体の挙動を数値モデルとして算出した結果を示す。矢印が羽根２０のまわりの流速分布を表わしている。なお、計算時間の都合上、図３６では、本体構造を略し、一方の羽根のみが動作するハーフモデルで解析した場合の流速分布を示している。また、表示の煩雑さを避ける都合上、羽根２０に垂直な断面についてのみ流速分布を表示してある。

また、この手法により算出される、ノード７１，８３，３３７におけるｙ軸方向の節点力の合計を示したのが図３７である。強制変位を与えているのはこの３点なので、この点における節点力の合計が、本体に加わる力である。最初不規則であった節点力が、最終的には周期的な挙動に収束している。すなわちこの時点で流体の挙動も構造の挙動も周期的となり、これはホバリング時の流体及び構造の挙動と等価である。但し、図４３に示すように、浮上力の発生する方向は、ｙ軸負方向である。筆者らの計算によると、重力加速度９．８ｍ／ｓｅｃ²の条件においては、図３７から分かるように、羽根１枚につき０．１ｇ程度の質量を浮上させることができる。

また、この際の羽根にかかる駆動トルクを図４６に示す。なお、θｙはｙ軸正の方向への右ねじ回転方向へのトルクを、θｚは、ｚ軸正の方向への右ねじ回転方向へのトルクを示す。

これは、アクチュエータの駆動トルクそのものであるので、以上のデータより、トルク１．５ｇｆ・ｃｍ〜３．５ｇｆ・ｃｍのアクチュエータを用いて、羽根１枚につき０．１ｇの質量を地上より浮上させることができると言える。

これにより、定常的なホバリングの際の流れの数値モデルを作成することができ、浮上可能重量や、これを実現するアクチュエータのトルクを求めることができる。

（解析結果のロボット作製への応用）
ここでは、流体を含めた羽ばたきの数値モデルのロボット制御方法への適用等について、図２７、図３７〜図４０を用いて説明する。

上述の手法により得られた数値モデルは、後述するＩ.アクチュエータ制御方法にて説明するが、そのままでも羽ばたき飛行ロボットの制御に用いることが可能である。また、上述の流体・構造連成解析により得られた数値モデルそのものから、昆虫の利用している空力を解明し、解明された空力を利用して羽ばたき飛行ロボットを製造することもできる。また、実際の羽ばたき飛行ロボットの製造に関して、感度解析を用いてこの流体・構造連成数値モデルを変更し、この感度解析により変更された数値モデルに基づいて羽ばたき飛行ロボットの羽根の構造や羽ばたき飛行の態様を決定することは工業上非常に有用である。

羽ばたきロボッのト制御方法の作成に関して、特に有用である、Ｉ.アクチュエータの制御方法の作成手法、および、ＩＩ.昆虫の羽根をベースとした昆虫の羽根と異なる羽根の制御方法の作成手法を以下説明する。

（アクチュエータ制御方法作成手法）
図３８は、羽ばたき飛行ロボット制御系を示した模式図である。

説明の簡便のため、アクチュエータ１２２とその制御系については、羽根１２１の打ち上げ、打ち下ろし方向に駆動する、自由度が１のアクチュエータである。また、本来複数枚あるべき羽根について、１枚の羽根だけの制御を例にして説明することにする。この羽１２１を昆虫の羽根と同様の運動をさせる手法を考える。なお、この羽ばたき飛行ロボットにおける羽ばたき飛行の態様等は羽根１２１の水平面に対してなす角度θの時系列データで表現することにする。

一般的に、アクチュエータの駆動は、アクチュエータのパワー（駆動力の数値）の指定によって行なわれている。したがって、アクチュエータの駆動力モデルがない場合、羽根の変位を測定して、羽根の変位を駆動力モデルに変換する機構が必要になる。なぜなら、アクチュエータの駆動力（パワー）が一定であっても、羽根にかかる負荷により羽根の変位は異なるからである。たとえば、図３８のアクチュエータ１２２の場合、羽根１２１の速度によって、流体より羽根１２１が受ける力が変化する等の理由により、羽根１２１をある姿勢に定位させる駆動力（パワー）は異なる。

また、羽根の変位、すなわち、羽根の所定の位置の時系列データをフィードバックすることにより制御するタイプのアクチュエータは、その性質上、羽根が目標位置に到達する時刻以前に、羽根が目標位置に到達したときに必要とする駆動力（パワー）を得なければならないため、余分なエネルギを必要とする。その結果、アクチュエータそのものの大きさが大きくなってしまう。

また、アクチュエータ１２２の制御に必要なパワーは、流体と構造との連成によって決まるので、この方法では、羽根を駆動するのにどれだけのパワーが必要かが、羽根構造の実物を製造してから初めて判明する。そのため、アクチュエータのパワーにかなりのマージンを見込んでアクチュエータ１２２の仕様を決定する必要がある。このため、たとえば、高速度カメラから撮影した画像より求めた羽根１２１の運動を実現する駆動を、直接アクチュエータに行なわせると、アクチュエータのパワーにかなりのマージンを持たせる必要が生じる。

逆に、アクチュエータ１２２の制御が、アクチュエータ１２２の発するパワー制御のみで可能ならば、必要かつ十分なパワーを有するアクチュエータを、必要かつ十分なエネルギで駆動することができる。また、この場合、羽根にかかる負荷を含めて羽根を制御しているので、羽根にかかる負荷による制御のズレを調整するための制御を必要とない。

そこで、上述の流体・構造連成数値モデルを利用すれば、直接アクチュエータ１２２を駆動すべきトルクが与えられるので、容易にアクチュエータ１２２の制御方法を作成することができる。

この手法によれば、より軽量なアクチュエータおよびエネルギ供給源によって飛行することが可能になるため、この手法を適用することは羽ばたき飛行ロボットには極めて有用である。

但し、羽ばたき飛行ロボットに、羽根の位置を検出する装置を付加し、その装置から羽根の位置の時系列データをフィードバックさせることは、たとえば、駆動誤差の吸収に用いるために有用であるので、必ずしも羽ばたき飛行ロボットに羽根の位置を検出する装置を設けなくてもよいということではない。

図３８におけるアクチュエータ１２２にかかる電圧とトルクとの関係を、アクチュエータ１２２に関して計測すれば、直接、制御に用いる電圧を決定することができる。そのため、アクチュエータドライブ回路１２３に向けて発生させる信号を、制御装置１２４に内包しておき、制御の際には、この電圧を出力するように、制御装置１２４がアクチュエータドライブ回路１２３に信号を送信すれば、所望の羽ばたき制御が実現される。

（数値モデルをベースとした羽ばたき飛行ロボットの製造）
多くの場合、人間が製造したアクチュエータと、昆虫の筋肉には特性の違いがあり、昆虫の筋肉よる駆動と同一の駆動が最適な駆動であるとは限らない。

また、羽根２０と等価な構造物を製造することは、実物の羽根の構造の複雑さを考慮すると、効率的ではない場合もある。

逆に、構造モデルや羽ばたき運動数値モデルを変更した場合に、駆動力モデルに与える影響が把握できるならば、トンボの羽根にトンボの羽ばたき方をさせる以外の羽ばたき飛行の態様の制御を導出することが可能になる。そのため、羽ばたき飛行ロボットのデザインや動作に多様な変更を加えることができる。

ここでは、感度解析と一般的に呼ばれる手法を応用して、上述のプロセスにより得られた個体Ｂの羽ばたき飛行の流体・構造連成数値モデルをベースとし、このベースの数値モデルを変形した羽根を有する羽ばたき飛行ロボットを製造する手法について説明する。

（各種感度解析について）
まず、羽根１３１の形状における感度解析について説明する。感度解析とは、羽１３１の形状に微妙な変更を加え、これを上記一連の流体・構造連成解析を行なった結果の流体に関する数値モデルおよび構造に関する数値モデルと、形状に変更加える前の流体に関する数値モデルおよび構造に関する数値モデルとを比較することにより、感度すなわち羽根１３１の形状の変更に対するあるターゲットパラメータの変化の度合を求めることである。言いかえれば、形状の数値モデルの変更により生じた所定の数値モデルの変化を求めることである。

たとえば、図３９における節点Ｐｔのｚ座標をＰｔからＰｔ×（１＋δｚ）と変更した数値モデルを用い上記流体・構成連成解析を行なったところ、図１３における羽１３１の支点反力の重力加速度方向成分Ｆｚが、Ｆｚ×（１＋δｆ）に変化したとすると、節点Ｐｔのｚ座標の変更に対する支点反力の重力加速度成分の変化、すなわち、感度はδｆ／δｚと定義される。

続いて、羽ばたき方、すなわち、羽ばたき飛行態様の感度解析について説明する。羽ばたき方の数値モデルは、時系列に変化する羽根の姿勢、すなわち、時系列的な羽根の位置データに感度解析のための微妙な変更を加える。それにより、上記一連の流体・構造連成解析を行なった結果の羽ばたき方での流体に関する数値モデルおよび構造に関する数値モデルと、変更前の羽ばたき方での流体に関する数値モデルおよび構造に関する数値モデルとを比較する。それにより、感度、すなわち、羽ばたき方の数値モデルの変更に対する、あるターゲットパラメータの変化を求める。言いかえれば、羽ばたき方の数値モデルの変更により生じた特定の数値モデルの変化を求める。

たとえば、図３８における羽根１２１の振幅θｗを、θｗ×（１＋δθ）と変更した数値モデルを用い、上記流体・構造連成解析を行なったところ、図３８における羽根１２１の支点反力の重力加速度方向成分Ｆｚが、Ｆｚ×（１＋δｇ）に変化したとすると、振幅θｗの変更に対する支点反力の重力加速度成分の変化、すなわち、感度はδｇ／δθと定義される。

これ以後、特に断らない限り、上記ターゲットパラメータとして、羽ばたき飛行にとって重要である、羽ばたき浮上力を採用することにする。より具体的にいうと、図３８における羽根１２１の羽ばたき１周期あたりの羽根１２１の支点にかかる重力加速度と反対の方向の力の平均値を、本実施の形態では特定ターゲットパラメータとする。この特定ターゲットパラメータの値が、正であれば、羽ばたき飛行ロボットは上昇し、負であれば、羽ばたき飛行ロボットは下降する。また、特定ターゲットパラメータの値がほぼゼロであれば、羽ばたき飛行ロボットはホバリングを行なう状態となる。そして、以後は、羽根の構造や羽ばたき方の変更前の羽根の実構造等価数値モデルを用いて作成された羽根の支点にかかる重力加速度と反対の方向の力の平均値の数値モデルと、羽根の構造や羽ばたき方の変更を行なった後の羽根の支点にかかる重力加速度と反対の方向の力の平均値の数値モデルとが等しくなる羽ばたき方を導出する手法について述べる。

（異なる形状の羽に対する羽ばたき方の導出）
まず、羽の形状を変更した際の羽ばたき方を求める手法について説明する。より具体的には、図４０に示すように、流体・構造連成の数値モデルの作成に用いた実構造等価数値モデルと異なる人工の羽根１４３の数値モデルに対して、実構造等価数値モデルにおける羽ばたき浮上力と等しい羽ばたき浮上力を得ることのできる羽ばたき方を導出する手法について述べる。

上記の感度解析の説明においては、節点を個別に付加したが、あるパラメータに基づき、まとめて羽根の構造に関する数値モデルを変化させることで、このまとまった羽根の構造に関する数値モデルの変更に対する特定のパラメータの感度を求めることができる。例として、羽根のサイズを拡大、縮小した場合の羽ばたき態様の変化、すなわち、羽根の形状変化に対する羽ばたき浮上力の感度を求める場合が考えられる。たとえば、補間比率をδ１として、羽根のサイズを（１＋δ１）倍した場合には、羽ばたき飛行ロボットの羽ばたき浮上力は（１＋δ２）倍になったとする。

また、羽ばたきの振幅θｗ（１＋δ３）倍にした際には、羽ばたき飛行ロボットの羽ばたき浮上力が（１＋δ４）倍になったとする。この結果より、羽根のサイズを（１＋δ１）倍した場合には、羽ばたきの振幅θｗを、（１−δ３×δ２／δ４）倍にすれば、羽根のサイズを変更する前と同等の羽ばたき浮上力を得ることができる。このように、羽根の形状の変化を１つのパラメータで置き換えて、この１つのパラメータの変更に対する浮上力の変化を求める方が、解析手法を簡略化することができる。

この手法をさらに発展させ、羽根の実構造等価数値モデル１４１の構造を表わすパラメータ全体を、１つのパラメータによって変更する場合を考える。図４０に示すように、人工の羽根１４３の構造（形状等）に関する数値モデルのパラメータを表わすパラメータ群Ｍｍ（ｉ）と、先のトンボの実構造等価数値モデルの構造に関する数値モデルのパラメータを表わすパラメータ群Ｍｎ（ｉ）とを定義する。

ここで、補間により新たに作成した微小変化数値モデル、すなわち、補間構造数値モデル１４２のパラメータ群Ｍｂ＿ｎｅｗ（ｉ）について、
Ｍｂ＿ｏｌｄ（ｉ）＝Ｍｎ（ｉ）
Ｍｂ＿ｎｅｗ（ｉ）＝Ｍｂ＿ｏｌｄ（ｉ）＋δ５×（Ｍｍ（ｉ）−Ｍｎ（ｉ））…（１）
とすれば、パラメータδ５を、補間比率と考えて、先のδ１と同様に扱うことができる。すなわち、Ｍｂ＿ｎｅｗ（ｉ）に対して、Ｍｎ（ｉ）と同じ羽ばたき浮上力を得る羽ばたきの振幅θｗを求めることができる。

さらに、解析の結果得られたＭｂ＿ｎｅｗ（ｉ）を新たにＭｂ＿ｏｌｄ（ｉ）として、再び、上記（１）式に戻って、各種感度解析とこれに基づいた羽ばたき振幅θｗの変更を繰返す。

これを繰返して、Ｍｂ＿ｎｅｗ（ｉ）を更新し、それに応じて羽ばたきの振幅θｗを修正していくので、新たなモデルＭｍ（ｉ）における、元の実構造等価数値モデルＭｎ（ｉ）と等しい羽ばたき浮上の態様を実現する羽ばたきの振幅θｗが求まる。

さらに具体的には、たとえば、図４０に示すように、人工の羽根１４３の構造数値モデルにおけるメッシュの相対的な配置は実構造等価数値モデル１４１と同一として、ノードの位置と各メッシュの厚さ、およびヤング率のみを変更していく。それにより、実構造等価数値モデルを人工の羽根１４３の数値モデルに近づくように修正を加えていく。たとえば、Ｍｍ（ｉ）、Ｍｎ（ｉ）、Ｍｂ＿ｏｌｄ（ｉ）、Ｍｂ＿ｎｅｗ（ｉ）に含まれるパラメータは、共通の配列構造であって、ノードの座標、ノードの姿勢、メッシュの厚さ、および、メッシュのヤング率を並べたものであってもよい。なお、Ｍｍ（ｉ）等のパラメータ群についてはここで示したものが一例であり、これに限るものではない。

たとえば、精度は落ちるがメッシュの相対的位置関係が異なる場合においても、何らかの補間を用いることで、上記の手法を適用することが可能である。

また、Ｍｍ（ｉ）などを、羽根の形状や剛性分布を表わす関数の係数群とすることも可能である。なお、ここに挙げた感度解析は、各パラメータの変化が線形であると想定できる、十分小さいδに対して行なわれることが前提となる。

また、上記の説明においては説明の簡便のため、羽ばたき方の変更について羽ばたきの振幅θｗのみの説明をしたが、実際には羽ばたきの振幅θｗを変更することによって、羽ばたき浮上力以外の浮上を表わすパラメータも変化する。

よって、実際には、羽ばたき方を表わす各種のパラメータの、羽ばたき力（羽ばたき浮上力以外の羽ばたきに関する力：たとえば、羽ばたき浮上力とともに発生する推進力など）に対する感度を算出し、補間された構造の数値モデルに対応する羽ばたき方の数値モデルとして、前述の羽ばたき力の変化を相殺するような羽ばたき方を表わすパラメータを求める。すなわち、構造の数値モデルの補間によって生じた羽ばたき力の変化を打消すうような羽ばたきを表わすパラメータの組を、アクチュエータ１２２の特性も考慮に入れて、線形計画法などを用いて決定することになる。具体的には、羽ばたき方を表わす各種のパラメータとは、羽根の根元近傍のノードに与えられる時系列で表わされた位置データであり、羽ばたき浮上力は、本体構造における並進加速度のデータおよび回転加速度データで表わされる。

前述の事項をまとめると次のようになる。

前述の羽ばたき飛行ロボットの製造方法においては、まず、人工的に作成された人工羽根１４３の構造に関する数値モデルを作成する。次に、補間の基準となる基準構造数値モデルＭｂ＿ｏｌｄ（０）＝Ｍｎ（０）を準備する。その後、基準構造数値モデルＭｎ（０）の運動態様に対応する基準運動数値モデルＮ（０）を準備する。この基準運動数値モデルＮ（０）は、予め人工的に作成された羽根の構造の数値モデルに対応する運動数値モデルであってもよいが、前述の昆虫の羽根の実物の測定により作成された実構造等価数値モデルに対応する運動数値モデルを用いれば、昆虫の羽根を基準として人工羽根１４３を設計することができるため、設計を効率的に行なうことができる。

次に、基準構造数数値モデルＭｎ（０）と基準運動数値モデルＮ（０）とを用いて、流体・構造連成解析を行なうことにより、基準構造数値モデルＭｎ（０）の構造に関する数値モデルＳ（０）および流体に関する数値モデルＦ（０）を算出する。

その後、人工羽根１４３の構造に関する数値モデルと、基準構造数値モデルＭｎ（０）とを、所定の補間比率にて補間した第１補間構造数値モデルＭｂ＿ｎｅｗ（１）＝Ｍｂ＿ｏｌｄ（０）＋δ５×（Ｍｍ（１）−Ｍｎ（０））を作成する。

また、構造に関する数値モデルＳ（０）および流体に関する数値モデルＦ（０）から、第１補間構造数値モデルＭｂ＿ｎｅｗ（１）と基準運動数値モデルＮ（０）とを用いて、流体・構造連成解析を行なったときの、構造に関する数値モデルＳ（１，０）および流体に関する数値モデルＦ（１，０）へ、変化した場合を考える。その変化により生じた特定の数値モデル（たとえば、浮上力）の変化の度合よりも、その特定の数値モデルの変化の度合が小さくなるように、基準運動数値モデルＮ（０）を変更する。それにより、第１補間構造数値モデルＭｂ＿ｎｅｗ（１）に対応する第１補間構造運動数値モデルＮ（１）を作成する。

その後、第１補間構造数値モデルＭｂ ｎｅｗ（１）と第１補間構造運動数値モデルＮ（１）とを用いて、流体・構造連成解析を行なうことにより、第１補間構造数値モデルＭｂ＿ｎｅｗ（１）の構造に関する数値モデルＳ（１）および流体に関する数値モデルＦ（１）を算出する。

次に、人工羽根１４３の構造に関する数値モデルと、第１補間構造数値モデルＭｂ＿ｎｅｗ（１）とを、所定の補間比率にて補間した第２補間構造数値モデルＭｂ＿ｎｅｗ（２）＝Ｍｂ＿ｏｌｄ（１）＋δ５×（Ｍｍ（２）−Ｍｎ（１））を作成する。

また、構造に関する数値モデルＳ（１）および流体に関する数値モデルＦ（１）から、第２補間構造数値モデルＭｂ＿ｎｅｗ（２）と第１補間構造運動数値モデルＮ（１）とを用いて、流体・構造連成解析を行なったときの、構造に関する数値モデルＳ（２，１）および流体に関する数値モデルＦ（２，１）へ、変化した場合を考える。その変化により生じた特定の数値モデル（たとえば、浮上力）の変化の度合よりも、その特定の数値モデルの変化の度合が小さくなるように、第１補間構造運動数値モデルＮ（１）を変更する。それにより、第２補間構造数値モデルＭｂ＿ｎｅｗ（２）に対応する第２補間構造運動数値モデルＮ（２）を作成する。

次に、第２補間構造数値モデルＭｂ＿ｎｅｗ（２）と第２補間構造運動数値モデルＮ（２）とを用いて、流体・構造連成解析を行なうことにより、第２補間構造数値モデルＭｂ＿ｎｅｗ（２）の構造に関する数値モデルＳ（２）および流体に関する数値モデルＦ（２）を算出する。

以後においては、ｉを順次繰上げていくことにより、補間構造数値モデルＭｂ＿ｎｅｗ（ｉ）が人工羽根１４３の構造に関する数値モデルと一致するかまたは近似されるまで前述の各数値モデル（補間構造運動数値モデルＮ（ｉ）、構造に関する数値モデルＳ（ｉ）および流体に関する数値モデルＦ（ｉ））を順次更新する。

それにより、補間構造数値モデルＭｂ＿ｎｅｗ（ｉ）を作成するステップ、補間構造運動数値モデルＮ（ｉ）を作成するステップおよび流体構造・連成解析を行なうことにより構造に関する数値モデルＳ（ｉ）および流体に関する数値モデルＦ（ｉ）を作成するステップを累積的に繰返す。

その結果、人工羽根１４３の構造に関する数値モデルと一致するかまたは近似された補間構造数値モデルＭｂ＿ｎｅｗ（一致または近似）に対応する補間構造運動数値モデルＭｂ＿ｎｅｗ（一致または近似）を得る。そして、補間構造運動数値モデルＭｂ＿ｎｅｗ（一致または近似）を用いて、人工的に作成された人工羽根１４３を駆動するアクチュエータドライブ回路１２３を制御する制御装置１２４の制御方法を決定する。

なお、補間構造運動数値モデルＮ（ｉ）は、補間構造数値モデルＭｂ＿ｎｅｗ（ｉ）に対応する運動数値モデルである。また、構造に関する数値モデルＳ（ｉ）および流体に関する数値モデルＦ（ｉ）は、補間構造数値モデルＭｂ＿ｎｅｗ（ｉ）の羽根を補間構造運動数値モデルＮ（ｉ）の羽ばたき方で羽ばたかせる条件で行なう流体・構造連成解析の結果得られる数値モデルである。

このような方法によれば、補間構造数値モデルＭｂ＿ｎｅｗ（ｉ）を人工羽根１４３の構造に関する数値モデルに近づけていくように補間していくことで、特定の数値モデル（たとえば、浮上力）に関しては人工羽根１４３の構造に関する数値モデルに、基準構造数値モデルＭｂ＿ｏｌｄ（０）（たとえば、実構造等価モデル）に対応する運動数値モデルＮ（０）での羽ばたき飛行態様に近い羽ばたき飛行態様を行なわせることが可能となる。その結果、人工羽根１４３に関しても、基準構造Ｍｂ＿ｏｌｄ（０）（たとえば、昆虫の実物）の羽ばたき飛行態様に近い羽ばたき飛行態様となるように、人工羽根１４３を駆動するアクチュエータドライブ回路１２３を制御する制御装置１２４の制御方法を決定することができる。

（羽構造の最適化）
上述の項目では、人工の羽根１４３構造を目指して羽根の構造を変更し、これに伴い羽ばたき方を変更することで、人工の羽１４３構造における羽ばたき方を導出する手法について説明した。しかしながら、感度解析を用いることにより、ある仕様を満たすように実構造等価数値モデルに変更を加えていくことで、新たな羽根構造を導出して作成する手法も可能である。

つまり、ある羽ばたき飛行ロボットに求められている仕様に適した羽根の構造を作成することができる。たとえば、羽根の形状についてのみ説明すれば、羽根の各節点について、ｘ、ｙ、ｚのそれぞれの方向の移動に対する、羽ばたき飛行の態様のパラメータの感度を算出し、このパラメータが最適な仕様に近づくよう、羽根の各節点を羽ばたき飛行の態様の感度に応じて移動させることを繰返す。

より具体的な例として、羽ばたきによる浮上力を最大にしたい場合は、羽根の各節点の移動に対する羽ばたき浮上力の感度を求め、その感度が増加するように各節点をその感度に応じて微小量移動させる。たとえば、感度が負である節点については、負の方向に移動させる。なお、感度は非線形に変化していくので、形状の変更のために感度解析を改めて行なう必要がある。また、上述のように羽根の特性を縮尺したパラメータを用いて感度解析を行なう手法も可能である。たとえば、ある羽ばたき周波数においてある羽ばたき浮上力を得ることのできる羽根の元の羽根に対する拡大比率を算出し、羽根構造を作成する場合等が考えられる。

（羽ばたき方の最適化）
上述の感度解析は、羽ばたき飛行ロボットの羽ばたき飛行態様の制御方法を最適化するために単独に用いることができる。

たとえば、図３７に示すように、昆虫の羽ばたきにおける浮上力には変動が大きい場合がある。これを、感度解析に基づき、羽根の姿勢、すなわち、羽根の位置の時系列のデータを変更をすることにより、羽ばたき浮上力を維持したままで、より変動の少ない羽ばたき方を導出することが可能である。このようにして、人工物に対応する要求仕様に基づいた制御方法を作成することが可能である。

さらに、羽ばたきロボットに所定の羽ばたき運動態様をもたらす羽ばたき駆動態様の数値モデルを作成する手法も考えられる。これは、羽根の駆動態様を表わすパラメータＷ（ｉ）について、Ｗ（ｉ）の変更に対する、羽ばたき飛行ロボットの運動態様の変化についての感度解析を行ない、解析結果に基づいて、羽ばたき駆動態様を変更していくことで行なわれる。

このように、羽ばたき駆動態様の数値モデルを作成する手法であれば、浮上する条件を満たしていない羽ばたき駆動態様の数値モデルであっても、容易に、羽ばたき駆動態様の数値モデルの性質を調べることができる。たとえば、羽ばたき飛行ロボットが右旋回をする場合、従来の、実験による手法では、浮上させている状態の中で、右旋回の可能な羽ばたき駆動態様を調べていく必要があった。しかしながら、本実施の形態の流体・構造連成解析により作成された羽ばたき駆動態様の数値モデルを用いた手法では、まず、右旋回する羽ばたき駆動態様の数値モデルを作成し、その後で、右旋回する羽ばたき駆動態様を行なうことに起因した浮上力の変化を補償するように、この羽ばたき駆動態様の数値モデルを変更していくという手法を採用することができる。

（アクチュエータ駆動の最適化）
上述の議論と同様に、アクチュエータの駆動についても変更を行なうことができる。たとえば、アクチュエータ１２２の応答が高周波数帯域で悪くなる場合など、ある一定以上の周波数でのアクチュエータ１２２の駆動を行ないたくない場合に、アクチュエータ１２２の駆動に要求される駆動周波数を、羽根の構造の変更により低くする手法を例として説明する。

高速度カメラにおける１秒当りの画像取得枚数、すなわち、サンプリング周波数を２×ｆｃとおくと、羽根１２１の駆動に関しては最も高い制御周波数はｆｃである。ここで、まず、羽ばたき方に関してｆｃ×（１−δｃ）になる周波数にてローパスフィルタリングを行なう。このように定まった羽ばたき方に基づき上述の流体・構造連成解析を行なう。さらに、フィルタリング周波数の変更に対する、羽ばたき浮上力の感度を算出する。また、上述の羽１２１の形状の変更による、羽ばたき浮上力の感度を算出する。

これらの感度により、フィルタリング周波数の変更によって起きる羽ばたき浮上力の変化を補うことができる、羽根１２１の変更された形状を導出することができる。

（羽ばたき飛行ロボットの製造）
上述の感度解析を用いる手法を組合せることにより、羽１２１の形状が実構造等価数値モデルと異なる羽ばたき飛行ロボットについても、前述の数値モデルをベースとした、同様の羽ばたき浮上力を得られる羽ばたき方を導出して作成することが可能である。

また、アクチュエータの駆動特性を考慮した羽根構造とその制御方法を作成することも可能である。なお、ここに示した流体・構造連成数値モデルのみならず、神経系モデルや筋肉の駆動モデル、昆虫の情報処理モデルについても、同様の手法を用いて、昆虫の当該数値モデルをベースとして人工物の仕様に基づいた新たな数値モデルを作成して、この数値モデルに基づいた羽ばたき飛行ロボットを製造することが可能である。これらの羽ばたき飛行ロボットの製造の手順をまとめて図２７に示す。

（その他）
（モデリング手順について）
トンボなどの羽根は、種が同じであれば、個体差はあるものの、概ね同一形状をしている。また、羽根は非常に複雑な形状をしており、それゆえに構造や剛性の緻密なモデリングには非常に時間がかかると考えられる。

なお、本実施の形態の流体・構造連性解析においては、羽根の構造や運動態様の最適化手法として、最も明示的である感度解析という手段を用いたが、最適化手法としての別の手法を用いる手法も可能である。他の有効な手法として、ニューラルネットワークを用いた学習を用いる手法や、遺伝アルゴリズムを用いた最適化を行なう手法が挙げられる。

たとえば、遺伝アルゴリズムを用いて羽根の設計を最適化する手法においては、まず、前述の各羽根の構造の数値モデルを表わすパラメータ群Ｍ（ｉ）の各要素を遺伝子としてコーディングを行なう。これにより表わされる羽根の構造の数値モデルが、空気中で所定の運動をする場合の流体・構造連性解析を行なう。これにより得られた結果を、たとえば、浮上力の大小といった、何らかの評価関数によってその良否を評価する。その評価に基づいて、よい構造の遺伝子をもつ数値モデル同士をかけ合わせる（組合せる）というプロセスを繰返すことで、最適な構造に近い羽根の構造の数値モデルが作成される。

また、この遺伝アルゴリズムを用いる手法は羽ばたき運動態様の数値モデルの作成にも応用することができる。すなわち、その一例として、羽根の運動を表わすパラメータをＷ（ｉ）とおき、このＷ（ｉ）に関して、前述のＭ（ｉ）と同様に遺伝アルゴリズムを適用する手法がある。その手法においては、まず、たとえば、羽根の水平面に対する角度を表わす関数ｗ、羽ばたきの１周期にかかる時間をＴ、ω＝２×π／Ｔ、最大サンプリング周波数をｆc、ｎを０≦ｎ≦Ｔ／ｆcとなる整数とする。このとき、羽根の運動を表わすパラメータをＷ（ｉ）を、周波数成分ごとに、ｗ＝Ａ（ｎ）×ｓｉｎ（ｎ×ω×ｔ）＋Ｂ（ｎ）×ｃｏｓ（ｎ×ω×ｔ）という式に分解した際の、各運動を表わすパラメータｗに対するＡ（ｎ）、Ｂ（ｎ）の集合を羽根の運動態様の数値モデルとする手法が挙げられる。

ところが、羽根は生体より切離されると急速に水分を失い、剛性が変化してしまう。また、同様の理由から、計測を行なっている期間中は、トンボの生理状態を維持する必要がある。

また、環境保護の観点から、データを採取したトンボは、計測終了後、自然に帰すことが望ましいとともに、計測中は昆虫の生理状態を維持する必要があり、計測に時間がかかる場合や、えさや水を与えることや、気温や湿度の管理などかなりの設備と手間が必要になる。

これらの理由により、羽根の剛性測定はトンボを損壊せずに行なうことが望ましいが、羽根を分解しないと、緻密な形状や剛性の計測を行なうことが困難となる。そこで、本実施の形態においては、羽根の緻密な計測に用いる個体と、羽根の駆動の計測に用いる個体とを別にすることでこれを解決している。なお、上記条件を完全に満たすことを必要としなければ、他の手法を用いて羽根と羽ばたき方のモデリングを行なうことは可能である。

たとえば、これらの測定が昆虫の測定上の寿命に対して十分短時間に行なえるなら、各個体に対して個別に羽根の実構造等価数値モデルを作成することが可能なので、これらの一連の羽根の基準構造数値モデルの作成のためのプロセスは必須ではない。

（流体・構造連成解析の手法について）
流体・構造連成解析には複数の種類の提案がなされている。

最も単純な手法として、高速度カメラにより撮影された羽根の映像より、流体のみについて移動境界問題を解くことによって流速を決定する手法が考えられる。また、羽根の構造は単独で変形を与えて解析することができる。

しかしながら、この手法においては、羽根のすべての部位（点）の移動を、すべての解析ステップにおいて計測する必要があり、データ量は膨大なものとなる。また、羽ばたき飛行の結果の解析に留まるため、上述の羽ばたき飛行ロボットなどの応用に用いることはできない。

流体と構造の相互依存を持った問題を解く手法として、流体と構造の支配方程式を交互に計算する弱連成法や流体と構造を含めた系全体の方程式（連成系方程式）を一挙（同時）に計算する強連成法が提案されている。前述した張群らによれば、本実施の形態のような流体・構造相互の依存性が強い問題に関しては、強連成法を用いることが最も効率的で最適であるとしている。

なお、前述した張群らの流体・構造連成構造解析は、生体以外の構造物についてなされたものであり、本実施の形態の流体・構造連成構造解析は、空気中での昆虫の羽ばたき運動を例とした生体の挙動についてなされたものである。このように生体に流体・構造連成構造解析を適用することにより、その生体を模倣したロボットを製造することが容易になることが本発明の特徴である。

（実施の形態２）
次に、表８および図４７〜図５２を用いて、実施の形態２の流体・構造連成数値モデルの作成方法を説明する。

本実施の形態の流体・構造連成数値モデルの作成方法において用いられる構造の数値モデルは、図３０に示す実施の形態１の流体・構造連成数値モデルの作成方法において用いられる構造の数値モデルの羽根の厚さが０．３５ｍｍ、０．１８ｍｍ、０．１５ｍｍ、０．１２ｍｍ、０．１ｍｍ、０．０５ｍｍである部分が、それぞれ、０．０４５ｍｍ、０．０１２ｍｍ、０．０１０ｍｍ、０．００８ｍｍ、０．００６ｍｍ、０．００４ｍｍに置き換えられている。

また、実施の形態１の流体・構造連成数値モデルの作成方法においては、羽ばたきの運動が略鉛直方向であったため、浮上力の発生する方向は、図４３において、ｙ軸正方向であったが、本実施の形態の流体・構造連成数値モデルの作成方法においては、羽ばたきの方向が略水平方向であるように変更されたので、浮上力の発生する方向はｚ軸正方向である。よって、浮上力はこの向きを正とする。

また、本実施の形態の流体・構造連成数値モデルの作成方法においては、実施の形態１の流体・構造連成数値モデルの作成方法において使用した図３４および図３５のメッシュ構造が、図４７および図４８に示すメッシュ構造に置き換えられている。

また、本実施の形態の流体・構造連成数値モデルの作成方法においては、実施の形態１の流体・構造連成数値モデルの作成方法において得られた図３７に示す浮上力と時刻との関係が、図４９に示すような関係になっている。

また、本実施の形態の流体・構造連成数値モデルの作成方法においては、実施の形態１の流体・構造連成数値モデルの作成方法において得られた図４４および図４５に示す位置と時間との関係が、図５０および図５１に示す関係になっている。

また、本実施の形態の流体・構造連成数値モデルの作成方法においては、実施の形態１の流体・構造連成数値モデルの作成方法において得られた図４６に示すトルクと時刻との関係が、図５２に示すような関係になっている。但し、図５２において、Ｔθはθ方向の駆動を行なうためのトルクであり、Ｔβはβ方向の駆動を行なうためのトルクである。

本実施の形態の流体・構造連成数値モデルの作成方法においては、実施の形態１の流体・構造連成数値モデルの作成方法で使用した表２のNode83のｘ座標値の０．００１が０．０００８８９に置き換えられており、ｙ座標値の０．００１が０．０００８８９に置き換えられている。

また、本実施の形態の流体・構造連成数値モデルの作成方法においては、実施の形態１の流体・構造連成数値モデルの作成方法で使用した表５のヤング率１．００Ｅ+０９が０．５Ｅ＋０９に置き換えられており、比重１．２Ｅ＋０３が０．７Ｅ＋０３に置き換えられており、厚さ０．３５Ｅ−０３が４．５Ｅ−０５に置き換えられている。

また、ShellMotion.datにおいて、羽根の節点におけるｘ、ｙ、ｚ座標の時刻歴は、表６の代わりに表８に示すようになる。

また、本実施の形態の流体・構造連成数値モデルの作成方法において用いられる構成および手法は、その他の図および表においては、実施の形態１の流体・構造連成数値モデルの作成方法と同様である。

なお、今回開示された実施の形態はすべての点で例示であって制限的なものではないと考えられるべきである。本発明の範囲は上記した説明ではなく特許請求の範囲によって示され、特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれることが意図される。

本発明の一実施の形態の浮上移動システムに用いられる浮上移動装置の構成を示す概略図である。 本発明の一実施の形態の浮上移動システムに用いられる浮上移動装置の駆動機構の主要部の構成を示す概略図である。 本発明の一実施の形態の浮上移動システムに用いられる浮上移動装置の突起およびストッパの機能を説明するための展開図である。 浮上移動に用いる羽の運動態様の一形態を示す概念図である。 本発明の一実施の形態の浮上移動システムに用いられる浮上移動装置の羽の動作を示す概念図である。 本発明の一実施の形態の浮上移動システムに用いられる浮上移動装置の羽の挙動を表わす概念図である。 本発明の一実施の形態の浮上移動システムに用いられる浮上移動装置のストッパ移動機構の構成を説明する断面図である。 本発明の一実施の形態の浮上移動システムに用いられる浮上移動装置の羽の剛性を、荷重に対する変位量を用いて示すグラフである。 本発明の一実施の形態の浮上移動システムに用いられる浮上移動装置の、羽の変形の様子を表わす概略図である。 本発明の一実施の形態の浮上移動システムに用いられる浮上移動装置の羽の角度と時間との関係とを表わすグラフである。 本発明の一実施の形態の浮上移動システムに用いられる浮上移動装置の方向制御の方法を示す概念図である。 本発明の一実施の形態の浮上移動システムに用いられる浮上移動装置の羽の形状を表わす図である。 本発明の一実施の形態の浮上移動システムの、羽の位置表現方法を示す概念図である。 本発明の一実施の形態の浮上移動システムに用いられる、free-flight解析のフローチャートである。 本発明の一実施の形態の浮上移動システムに用いる浮上移動装置の、羽の厚み分布を示す概念図である。 本発明の一実施の形態の浮上移動システムに用いる浮上移動装置の、羽ばたき駆動信号におけるduty比と駆動周波数の関係を表わすグラフである。 本発明の一実施の形態の浮上移動システムの、羽の位置を表現する節点の位置表現方法を説明する概略図である。 本発明の一実施の形態の浮上移動システムのシミュレーションに用いられるメッシュ分割を表わす図である。 本発明の一実施の形態の浮上移動システムの動作説明のための問題設定を表わす概略図である。 本発明の一実施の形態の浮上移動システムに用いる浮上移動装置の、羽ばたき駆動における、打ち上げ周波数に対する打ち下ろし周波数と、前後方向の流体力の関係を表わすグラフである。 本発明の一実施の形態の浮上移動システムの、ある目標経路を実現するための羽ばたき駆動手法を例示するグラフである。 本発明の一実施の形態の浮上移動システムに用いる浮上移動装置の、羽の駆動方法を表わすグラフである。 本発明の一実施の形態の浮上移動システムに用いる浮上移動装置の、シミュレーションにより求まった経路と目標経路の関係を示す概略図である。 本発明の一実施の形態の浮上移動システムに用いる浮上移動装置の、羽ばたき駆動方法修正のための角度差を示す概略図である。 本発明の一実施の形態の浮上移動システムのデーターフロー図である。 本発明の一実施の形態の浮上移動システムの処理のフローチャートである。 実施の形態１に係る、流体・構造連成数値モデルの作成方法および流体・構造連成数値モデルに感度解析による変更を加えた羽ばたき飛行ロボット製造プロセスを表わす図である。 実施の形態１に係る、実際の昆虫の羽根における梁構造および膜構造の状態を示す図である。 実施の形態１に係る、羽根の形状の計測手法を説明するための図である。 実施の形態１に係る、羽根の厚みの分布を表わす図である。 実施の形態１に係る、羽根におけるほぼ平坦な領域を表わす図である。 実施の形態１に係る、羽根における代表点での剛性を測定する１つの手法を説明するための図である。 実施の形態１に係る、羽根における代表点での剛性を測定するさらに他の手法を説明するための図である。 実施の形態１に係る、流体のメッシュの作成方法を説明するための図である。 実施の形態１に係る、流体のメッシュの外形線を表わす図である。 実施の形態１に係る、流体・構造連成数値モデルのうち流体の挙動を示す図である。 実施の形態１に係る、流体・構造連成数値モデルが羽ばたいた場合の羽根の支点にかかる支点反力の推移を示す図である。 実施の形態１に係る、羽ばたき飛行ロボットの概略構成を示す図である。 実施の形態１に係る、感度解析の手法を説明するための図である。 実施の形態１に係る、羽根の形状を補間する手法を説明するための図である。 実施の形態１に係る、実構造等価数値モデルの具体例を説明するための図である。 実施の形態１に係る、実構造等価数値モデルの具体例を説明するための図である。 実施の形態１に係る、実構造等価数値モデルの具体例を説明するための図である。 実施の形態１に係る、実構造等価数値モデルの具体例を説明するための図である。 実施の形態１に係る、実構造等価数値モデルの具体例を説明するための図である。 実施の形態１に係る、実構造等価数値モデルにかかる駆動トルクを説明するための図である。 実施の形態２に係る、流体のメッシュの作成方法を説明するための図である。 実施の形態２に係る、流体のメッシュの外形線を表わす図である。 実施の形態２に係る、流体・構造連成数値モデルが羽ばたいた場合の羽根の支点にかかる支点反力の推移を示す図である。 実施の形態２に係る、実構造等価数値モデルの具体例を説明するための図である。 実施の形態２に係る、実構造等価数値モデルの具体例を説明するための図である。 実施の形態２に係る、実構造等価数値モデルにかかる駆動トルクを説明するための図である。

符号の説明

１ 浮上移動装置、２ アクチュエータ、４ 羽、５ 制御装置、６ 電源、９ 支持構造、２１ ステータ、２２ ロータ、３１,３２ ベアリング、３１ａ，３２ａ 外周部、３１ｂ，３２ｂ 内周部、３３ 伝達軸、３４ 突起、３５１，３５２ ストッパ、３６１，３６２ ストッパ移動機構、３７ 超音波リニアアクチュエータ、３８ テフロン（Ｒ）ベアリング、５１ ６軸加速度センサ、 ５２ 近傍形状計測センサ、 ５３ ＲＯＭ、５４ ＲＡＭ、５５１ 羽駆動演算部、５５２ Free-flightシミュレーション演算部。

Claims (29)

1. 羽部とこの羽部を駆動する駆動部、およびこの駆動部を制御する制御部とを有し、この駆動部を用いて、当該羽部を能動的に運動させることで得られる周囲流体からの反力によって移動を行なう停空可能な浮上移動装置を用い、前記制御部において流体・構造連成解析を用いたfree-flightシミュレーションを行ない、この結果を用いて、制御を行なうことを特徴とする浮上移動システム。
2. 前記浮上移動装置の構造における数値モデルを予め有することを特徴とする請求項１に記載の浮上移動システム。
3. 前記駆動部の出力を、前記制御部からの入力信号に関連づけて予め有することを特徴とする請求項１に記載の浮上移動システム。
4. 前記関連づけについて、前記駆動部に対する負荷の影響が考慮されていることを特徴とする請求項３に記載の浮上移動システム。
5. 前記駆動部の出力に基づき、前記構造モデルの節点における境界条件を与えることを特徴とする請求項２〜４のいずれかに記載の浮上移動システム。
6. 前記駆動部の変位出力に基づき、前記構造モデルの節点における変位を境界条件として与えることを特徴とする請求項５に記載の浮上移動システム。
7. 前記駆動部の力出力に基づき、前記構造モデルの節点における外力を境界条件として与えることを特徴とする請求項５に記載の浮上移動システム。
8. 前記free-flightシミュレーション内でのモデルに置いて、羽部と胴体部の接合部は、モデル上は節点を共有せず、計算上、前記接合すべき節点間で、節点外力もしくは節点変位を相互に与えることを特徴とする請求項１に記載の浮上移動システム。
9. 胴体のモデルを流体・構造連成解析中に含まず、構造についての計算のみ行なうことを特徴とする、請求項８に記載の浮上移動システム。
10. 前記流体・構造連成計算に置いて、羽の位置が胴体に対する相対位置で与えられることを特徴とする請求項１に記載の浮上移動システム。
11. 前記流体・構造連成計算に置いて、胴体の運動が、羽の節点における節点力によって定められることを特徴とする請求項１に記載の浮上移動システム
12. 前記流体・構造連成計算に置いて、羽の節点における外力を、当該羽の駆動装置から発生する力に基づき与えることを特徴とする請求項１に記載の浮上移動システム。
13. 前記流体・構造連成計算に置いて、当該計算における流体モデル全体および胴体モデルに対して、オフセットを与えて行なうことを特徴とする請求項１に記載の浮上移動システム。
14. 前記流体・構造連成計算に置いて、胴体の運動ベクトルのマイナス１倍の運動ベクトルを、当該計算における流体モデル全体および胴体モデルにオフセットとして与えることを特徴とする請求項１３に記載の浮上移動システム。
15. 羽部とこの羽部を駆動する駆動部、およびこの駆動部を制御する制御部とを有し、この駆動部を用いて、当該羽部を能動的に運動させることで得られる周囲流体からの反力によって移動を行なう停空可能な浮上移動装置であり、前記制御部において流体・構造連成解析を用いたfree-flightシミュレーションを行ない、この結果を用いて、制御を行なうことを特徴とする浮上移動装置。
16. 通信手段を有し、前記流体・構造連成解析を、物理的に独立した計算部にて行なうことを特徴とする、請求項１５に記載の浮上移動装置。
17. 形状計測手段を有し、これにより計測された周囲障害物の形状を、前記流体・構造連成解析における境界条件として与えることを特徴とする、請求項１５に記載の浮上移動装置。
18. 胴体部と、
    前記胴体部に取り付けられた羽部と、
    前記羽部を駆動する駆動部と、
    前記駆動部を制御する制御部と、
    周囲環境の形状を測定する環境センサと、
    前記胴体部、前記羽部および前記駆動部のうち少なくともいずれか１つの運動に関する物理量を測定する運動センサとを備えた浮上移動システムであって、
    前記制御部は、
    前記羽部および前記胴体部のそれぞれの構造に関する数値モデル、前記周囲環境の流体に関する数値モデル、および、前記駆動部が前記羽部を駆動する駆動態様と前記駆動部が前記駆動態様で駆動しているときに前記胴体部にかかる力との関係を示す関係データとを有する記憶部と、
    前記周囲環境の形状のデータ、前記運動に関する物理量のデータ、前記構造に関する数値モデル、前記流体に関する数値モデル、ならびに、前記関係データを用いて、流体の挙動を無視して、前記駆動部の最適駆動方法を決定する手段と、
    前記最適駆動方法のデータを用いて、流体の挙動を考慮して、流体・構造連成有限要素法による羽ばたき飛行のfree-flight解析の結果に基づいて、前記最適駆動方法のデータを修正し、修正された駆動方法のデータを取得する手段と、
    前記修正された駆動方法のデータを用いて、前記駆動部を駆動する手段とを含む、浮上移動システム。
19. 前記流体・構造連成有限要素法によるfree-flight解析における流体・構造連成解析においては、前記胴体および前記駆動部の構造に関する数値モデルが用いられることなく、前記羽部の構造に関する数値モデルのみが、前記構造に関する数値モデルとして用いてられる、請求項１８に記載の浮上移動システム。
20. 前記浮上移動システムは、前記胴体、前記羽部、前記駆動部、前記制御部、前記環境センサ、および前記運動センサの全てが１つのロボットに搭載されたスタンドアロン型の羽ばたき飛行ロボットである、請求項１８に記載の浮上移動システム。
21. 前記浮上移動システムは、
    前記胴体、前記羽部、前記駆動部、前記制御部、前記環境センサおよび前記運動センサを含む羽ばたきロボットと、
    前記羽ばたきロボットとは別個独立して外部に設けられた前記制御部とを有し、
    前記羽ばたきロボットは、前記制御部と通信することによって、羽ばたき飛行を実行する、請求項１８に記載の浮上移動システム。
22. 前記形状センサは、前記羽ばたきロボットに設けられた複数のＣＣＤ（Charge Coupled Device）であり、前記運動センサは、前記羽ばたきロボット自身の加速度を測定するセンサである、請求項２０または２１に記載の浮上移動システム。
23. 前記関係データは、ホバリング時の羽ばたき方を基本データとし、前記基本データに対して変更を加えることによって他の羽ばたき方が特定される形式で前記記憶部に格納されている、請求項１８に記載の浮上移動システム。
24. 前記制御部は、
    前記周囲環境の形状のデータ、前記運動に関する物理量のデータ、前記構造に関する数値モデル、前記流体に関する数値モデル、および、前記関係データを用いて、羽ばたきロボットにとっての最適経路を算出し、
    前記最適経路のデータに基づいて、前記最適駆動方法のデータを算出する、請求項１８に記載の浮上移動システム。
25. 前記流体・構造連成有限要素法による羽ばたき飛行のfree-flight解析においては、
    前記修正された駆動方法によって前記駆動部が駆動されたときの前記羽ばたきロボットのシミュレーション上の飛行経路が１回または複数回算出され、
    前記シミュレーション上の飛行経路と前記最適飛行経路との差が所定値以下であれば、前記制御部は、前記修正された駆動方法を用いて、前記駆動部を制御する、請求項１８に記載の浮上移動システム。
26. 前記運動に関する物理量には、前記羽ばたきロボット自身の移動の速度ベクトルが含まれ、
    前記制御部は、前記羽ばたきロボット自身の移動の速度ベクトルとは逆向きで大きさが同じ速度ベクトルを前記羽ばたきロボットの周囲の流体の対地速度ベクトルに加えて、前記流体・構造連成有限要素法によるfree-flight解析を実行する、請求項２０または２１に記載の浮上移動システム。
27. 前記流体・構造連成有限要素法による羽ばたき飛行のfree-flight解析においては、
    前記駆動部の力の出力態様または変位の出力態様のデータを用いて、前記構造に関する数値モデルの節点における境界条件として外力または変位のデータが作成される、請求項１８に記載の浮上移動システム。
28. 前記流体・構造連成有限要素法による羽ばたき飛行のfree-flight解析においては、
    前記胴体部の１つの節点と該胴体部の１つの節点と接合される前記羽部の１つの節点とに、同一節点変位、または、方向が逆向きで互いに大きさが等しい節点外力が与えられるものとして、解析が行なわれる、請求項１８に記載の浮上移動システム。
29. 前記修正された駆動方法のデータを取得する手段は、
    前記運動に関する物理量のデータを用いて、第１所定時刻における前記胴体部の第１絶対位置を算出する手段と、
    前記運動に関する物理量のデータを用いて、前記第１所定時刻における前記胴体部の第１姿勢を算出する手段と、
    前記最適駆動方法のデータを用いて、仮に、前記第１所定時刻と前記第１所定時刻よりも後の第２所定時刻との間の前記胴体部に対する前記羽部の相対変位を算出する手段と、
    前記第１絶対位置のデータ、前記第１姿勢のデータ、および前記羽部の相対変位のデータを用いて、前記第２所定時刻における前記胴体部に対する前記羽部の絶対位置を算出する手段と、
    前記羽部の絶対位置のデータを用いて、前記胴体部にかかる並進力および前記胴体部にかかるトルクを算出する手段と、
    前記並進力のデータおよび前記トルクのデータを用いて、前記第２所定時刻における前記胴体部の第２絶対位置および前記第２所定時刻における前記胴体部の第２姿勢を算出する手段と、
    前記第１絶対位置のデータおよび前記第１姿勢のデータと前記第２絶対位置のデータおよび前記第２姿勢のデータとの差を算出する手段と、
    前記差が所定値以下であるか否かを判別する手段と、
    前記差が所定値以下である場合に、前記羽部の相対変位のデータを前記修正された駆動方法のデータとして用いて、前記駆動部を制御する手段とを有する、請求項１８に記載の浮上移動システム。

Images