2014-09-04

人力飛行機の「安定」について

こんにちは。yuukivelです。

某琵琶湖での大会の放映が終わりましたね。昔出ていたTT部門は変わらず放送時間TT部門だったわけですが（笑）

さて、今回は人力飛行機の安定性について少し記事にするとしよう。

2010年に、teamFというメチャメチャ速かったチームの機体が、「上反角を抑える」「エルロンを積極的に用いる」etc...といった設計思想だったことから、その後の学生チームはこぞって「上反角を抑えて不安定にする」設計思想を採用した。確かに感覚的に不安定だとロール角をつけやすいと感じるのかもしれない。

ただし、それは、エルロンがきちんと効果を発揮して、それでいてパイロットが上反角がなくてずるっずるに滑る機体を制御できる場合に限る。

teamFに対して、今年優勝したチームは上手くて、エルロンがきちんと効いていてかつ上反角も十分で、ものすごく扱いやすい機体に仕上がっていたらしい。

式とか難しいお話に移る前に、ラジコン飛行機において、上反角が不足している場合の飛び方を動画で見てもらおう。（先輩ごめんなさい！）

この機体が無尾翼であるということはちょっと置いておいて、横（ロールとヨー）の運動に関して。動画を見てもらえれば分かるとおり、完全なスパイラル不安定だ。この機体はラダー機ではなくてエレボンで操舵するのだけれど、スパイラルダイブに陥ってしまうと抜け出すのがとても難しい。

ここからちょっと式を出して小難しいお話。
参考文献は言わずと知れた名著「航空機力学入門」

航空機力学入門

作者: 加藤寛一郎
出版社/メーカー: 東京大学出版会
発売日: 1982/11
メディア: 単行本
クリック: 15回
この商品を含むブログ (4件) を見る

今回のお話は二つ

スパイラル不安定について
ラダー機の上反角効果（航空機力学では $L_\beta$ の大小で表す)について

ここで、ラダー機とはラダーのみでロールを制御する機体。
先に僕の考えを書いておくとラダー機でもエルロン付きでも、きちんと上反角を着けた方が操縦しやすくて良いと思う。ラジコンでも上反角がつけばつくほど操縦しやすくなるし、エルロンによって揚力差がつけば、翼の剛性如何に関わらず、自ずと機体のロールは始まる。スパイラルダイブに入ったら無理だけれど。

用いる運動方程式

$-L_\beta+(\frac{d}{dt}-L_p)p-((I_{xy}/I_{xx})\frac{d}{dt}+L_r)r=L_{\delta a}\delta_{a}+L_{\delta r}\delta_{r}$
$-N_\beta-((I_{xz}/I_{zz})\frac{d}{dt}+N_p)p+(\frac{d}{dt}-N_r)r=N_{\delta a}\delta_{a}+N_{\delta r}\delta_{r}$

これは、ロール角速度とヨー角速度に関する線形化された運動方程式である。本当は $\beta$ （横滑り角)に関する式もあって、この方程式が実は、「速く曲がりたいなら機首を下げろ」の根拠になるのだけれど、今回はラダーによってある横滑り角になったとして話を進めるので、略。

//////////2014/9/6 訂正////////
「早く曲がりたいなら機首を下げろ」について $\beta$ の微分方程式のみ影響するような書き方になっていますが、実際は動座標系における見かけの力の影響によるもので、この式だけの影響ではありません。申し訳ありません。
///////////////////////////////

Iは慣性モーメント。添え字は軸。xzとかは連成項。
Lは有次元ローリングモーメント係数。添え字は何で微分したかを表す。（もちろん線形運動方程式なのでこれらは定数）例えば $L_\beta$ はローリングモーメントを横滑り角で微分した値に比例する。何度も書いている通り、上反角効果だ。
Nは有次元ヨーイングモーメント係数。添え字についてはLと同様。
また、 $\delta$ について。これは操舵機器の舵角による微分値を表す。今回はa:エルロン r:ラダーだ。だから $N_{\delta r}$ はヨーイングモーメントをラダー舵角で微分した値に比例する。これはラダーの効きを表している。

有次元係数は、方程式が見やすいように、慣性モーメントで除されている。だから慣性モーメントが小さければ有次元係数は大きくなる

さて、これらの方程式から、仮定によって項を落としていこう。こちらの式はラダー機の上反角効果を考えるときに用いる
まず、簡単のために、 $I_{xz}$ と $I_{xy}$ を0としよう。
そして、エルロンのヨーモーメント（アドバースヨーの原因)とラダーのローリングモーメントを無視しよう。

そして、ラダー機なのでエルロンの効果を外そう。

簡単化した結果は以下の通り
$-L_\beta\beta+(\frac{d}{dt}-L_p)p-L_rr=0$
$-N_\beta\beta-N_pp+(\frac{d}{dt}-N_r)r=N_{\delta r}\delta_{r}$

スパイラル不安定について

スパイラル不安定＝スパイラルモード不安定とは、機体が旋回するにつれてロールが大きくなっていく運動を指す。そして、そのまま回復不可能になって墜落するのがスパイラルダイブだ。
人力飛行機界隈ではスパイラルダイブを極端に恐れていて、それもあって短スパン化が進むのだろう。（僕はそれよりも短スパン化による出力増加の方が圧倒的に問題だと思っているが)

スパイラル不安定によって旋回中にロール角とロール角速度がある値以上になると、エルロン、ラダーの操舵力をオーバーして復帰不可能になる。だからこのスパイラル不安定を小さくすることが、タイムトライアル部門では有効だ。もちろんバカでかいエルロンとラダーで、不安定が強い機体でも操舵力を大きくして旋回を可能にするという解もある（それがまさに冒頭のteamFの設計思想だ）

航空機力学入門によると
スパイラルモードの根は以下のように近似される。（根とは制御用語で、複素平面における根の位置によって、運動の減衰や振動数が決定される)
$-\lambda_s\simeq-E/D$
根は虚軸の左側にあればそれは安定な根だから
$E/D>0$ であれば安定である。
$D$ は一般的に正であるらしい。従って $E$ の正負がスパイラルモード安定に寄与する。正であれば安定で、大きければ大きいほどスパイラルモード安定は強い。逆に小さければ小さいほどスパイラル不安定は強く、0から負に遠く離れてしまう事は避けなければならない。

$E$ は以下のように表される。
$E=(g/U_0)(L_{\beta}N_r-N_{\beta}L_r)$
これは本当は安定軸での記述なのだけれど、人力飛行機の設計であれば普通機体軸と安定軸は一致するのが普通。飛んでいるときのテールパイプは水平になるようにするため。

これより、 $L_\beta$ と $N_r$ が大きければ安定方向になる。上反角効果とラダー減衰（動ファクタ比と関連がある）が大きければ安定方向にむかう。
逆に $N_\beta$ と $L_r$ が小さければ安定方向となる。すなわちでっかい垂直尾翼と長いスパンは不安定の原因となる。

人力飛行機はこの $L_r$ が大きくて（スパンが長いから)スパイラルモードを安定にすることは難しい。不安定をいかに小さくするかが勝負だ。

そこで上反角効果が重要になる。上反角を大きくすることで、スパイラルモードは安定に向かう。

私が昔設計した機体の肝はそこだった。 $L_r$ はもともと大きいし、今更スパンを伸ばした所で大きくは変わらない事が解析してみて分かった。だからスパンを伸ばすことで旋回時必要な出力を小さくしつつ、上反角を十分に確保することで旋回を可能にしようという設計思想だった。

スパイラルモードについてよく誤解をもたれるのが、「スパイラルダイブというのは、旋回中の内側と外側の速度差によって起こる」という理解だ。これは一部正しい。そのまま $L_r$ の説明だからだ。しかしそれだけではなく、スパイラルダイブは上反角や垂直尾翼の減衰等の関連によって起こる。だからスパンが長いなら、ちゃんと上反角をつける、というのが操縦しやすい機体の定石となる。

ラダー機の上反角効果について

「不安定な機体は望む方向のロール角速度を得やすい」本当だろうか。
ここで不安定な機体とは、上反角効果の小さい機体を意味する。

先ほど出した簡易化した運動方程式を持ってこよう
$-L_\beta\beta+(\frac{d}{dt}-L_p)p-L_rr=0$
$-N_\beta\beta-N_pp+(\frac{d}{dt}-N_r)r=N_{\delta r}\delta_{r}$

そして、ある一定のロール角速度、ヨー角速度のとき、ある一定の横滑り角をラダーによって与えた時の挙動を考えるので、下の式は不要となり、旋回にてロールを変化させるため運動を考えるのに必要な方程式は
$\frac{d}{dt}p=L_\beta\beta+L_rr+L_pp$
まで単純化される。

ここまで来れば簡単だろう。ロール角速度を素早く変化させる（ロール角加速度を大きくする）ためには、上反角効果を大きくすることが重要となる。
右辺のうち、 $L_rr+L_pp$ は運動を阻害しようとする項だ。何度も言っているとおり、スパンが長いほどこれらの項は大きくなる。

私はよく「ラダー機において、上反角が大きい事はエルロンが大きい事に等しい」と言っている。その根拠がこれだ。
信じられない人は、是非いちど、上反角のついていないラジコン機を作って飛ばしてみるといい。すぐに上反角を付けたくなること請け合いだ。ラダー機だと操縦不能だから。

ここでこの章の結論を出そう。「不安定な機体が曲がりやすいわけではない」
$L_\beta$ を大きくする方法としては、上反角を大きくする事の他に慣性モーメントを小さくするという解もある。しかし運動を阻害する項も全部慣性モーメントで除されているので、あまり効果は上がらない。
安定な機体だからといって、決して曲がりにくい訳ではないことがきちんと理解されて欲しいと思う今日この頃。

終わりに

いかがだっただろうか。この手の話はかなりの頻度でされているので、今回は式を出すこともいとわずにきちんと議論してみた。
仮定が雑なところもあるが、大学の低学年の人も読むかもしれないことを考えてできるだけ制御論を出さなくても分かるようにしたのでご容赦願いたい。

機械、特に飛行機は、だいたい感覚と合っているのだけれど、時々クリティカルなところで感覚と合わない。これは仕方ないことだ。
だから数式をきちんと理解、せめて傾向だけでも把握して、実際に飛ばした結果やラジコンとかその他飛行機の感覚も勘定して、直観だけでなく、きちんとした理論的な根拠を与えて設計するのがいいと思うのだが。。。

まぁそんなこと言いまくっていると、「あいつは理論派」とか言われてしまうわけで、、、感覚と理論の擦り合わせ、そして理論屋さんとして絶対に譲れない所ははっきり主張していく、といったバランスが大切なんだろうな

2014-08-29

翼型の”美しさ”について&機械学習

翼型について

ご無沙汰しておりました。yuukivelです。

長かった大学院入試がやっと終わりました（面接まだあるけれども）(実はこれを書いているのは面接の待ち時間だったりする)いやはや大変だった。
とにかくやりたいことできず、つらいわでほんと、大変だったとしか言えないです。
やりきった感はあるけれども。合格発表までしばしやりたいことやるとします。

さて、本題に入ろう。

今回のタイトルはいかがわしい（？）感満載だが、「翼型は美しい！」って事が言いたいのではく、ちゃんと工学的なお話なので安心を。。。

早速だが、下に示す翼型のうち、どちらが”美しい”と言えるだろうか。

たぶん大半の方が、「上の方が美しい」というのではないだろうか。

そこでもう一つ質問しよう。

なにを基準に美しいと判断したのか。

この質問は難しい、本当に難しい。

そこで、今回の記事のテーマは「翼型の美醜をいかに数値的に表すか」というテーマだ。ちょっとArtisticかな？中二病だからしかたないね！

主に3DCADプログラミングの界隈で、曲線の美しさをいかに定義するかということに関して多くの先行研究がある。
歪みエネルギー最小化曲線網に基づく曲面の生成
といった論文を紹介しておこう。

ここの論文にあるとおり、曲線の滑らかさについては、曲線の歪みエネルギー
$R=\frac{x''y'-x'y''}{((y')^2+(x')^2)^{3/2}}$

$E=\int{R^2dx}$

で求められる。ここで翼型の場合、前縁付近で
$\frac{dy}{dx}=\infty$
となるので、精度は落ちるけれども、x,yについての差分表現とした。

実際に翼型の歪みエネルギーを計算してみると結構面白くて、事実、上に表示した翼型の歪みエネルギーはそれぞれ
上：0.1351
下：0.3768
となり、上の翼型が圧倒的に翼型歪みエネルギーが小さくなる。

他にも、曲率の連続性が曲線の美しさに影響するとあったり、様々なようだ。
STUDIOの曲率連続性評価 - Autodesk Knowledge Network

ここで、とりあえず美しさの一つのパラメータとして翼型歪みエネルギーがあることが分かったが、これだけではないと思われる。
キャンバーや翼厚の分布のさせ方やその値自体、最大キャンバや最大翼厚の位置といったことも影響してくるのではないだろうか。

そこで、翼型に関する私自身の好みをはっきりさせ、翼型設計に活かしていくために、これらのパラメータを収集し、最近はやりの
「機械学習」
を用いて、自分の好みの翼型と好みじゃない翼型の判別（点数付け）を行ってみた。

なぜこれをしたかというと、このブログでさんざん述べてきた「遺伝的アルゴリズムによる翼型設計」と深く関係がある。

現在、距離競技機の設計を行っているけれども、そこで現在公開しているXGAGの改良として、
「レイノルズ数（翼弦長）および有次元翼厚の最適化を含めた遺伝的アルゴリズムによる翼型設計コード
『Tubdotgag』」
を開発した。
（開発者専用試作設計コードとかいえばかっこいいでしょ？えへへ)
（設計が落ち着いたらこちらも公開しようと思います。ただし最適化が強すぎてちょっと使いにくいのは確か。またXGAGにもましてピーキーな翼型がでることがあるので、そうした翼型を弾くためにも翼型の速度分布や圧力分布、境界層パラメータを見る事ができる知識が必要です）

遺伝的アルゴリズムによる翼型設計は、最終的には翼型設計者が形状を見て、その翼型を採用するかを決めることが多い。つまりはいかに性能が良くてもその翼型が”汚”ければ採用されない。
そこで、GA内の評価関数に翼型の美醜を組み込む事によって、もう最初から汚いのはいらない、としてしまえば良いのではないかと思ったのが事の始まりだ。
確かに探索範囲が狭くなってしまって、たとえば汚いと判断された翼型の中にも美しいものがあって、素晴らしい性能を持つものがあるかもしれない。GAのコンセプトとして系の多様性をできるだけ大きくするのが重要だが、少しそれには反するのかなと。
しかしこの効果は大きくて、実際に組み込んでみると、本当に好みの翼型ばかり出るようになった。

さて、実際にどのように翼型の美醜を組み込んでいるのか見てみよう。

まず、収集するデータだ

翼型歪みエネルギー
翼型キャンバ分布歪みエネルギー
翼型翼厚分布歪みエネルギー
前縁曲率
後縁角度
最大キャンバ
最大キャンバ位置
最大翼厚
最大翼厚位置
美醜どちらかの（人の）判定結果

の10個だ。何が影響するかは分からないので、とにかく雑多にデータを集めてみた。曲率の連続性等も今後とれればいいなと。

収集したデータはフリーの機械学習ソフト「Weka」で読み込むことができるフォーマットになっている。
http://www.cs.waikato.ac.nz/ml/weka/
データをWekaで可視化するとこんな感じ

また、Wekaについて参考にした書籍は以下の通り。

フリーソフトではじめる機械学習入門

作者: 荒木雅弘
出版社/メーカー: 森北出版
発売日: 2014/03/29
メディア: 単行本（ソフトカバー）
この商品を含むブログ (4件) を見る

評価関数に組み込むためには、美醜を数値化しなければならない。その数値化には「マハラノビス距離」を用いた。
マハラノビス距離とは感覚的に、「多次元正規分布において、その分布の山の中で何合目あたりにいるか」の指標だ（すみません僕自身初心者なのであまりよく分かっていないのです。。。）
つまりは、ある値、例えば翼型歪みエネルギーのデータをとって、その平均を求める。データが平均値付近に集まっていれば、平均値から大きく外れた個体のマハラノビス距離はとても大きくなる。
逆にデータが平均から散らばって存在していれば、普通の距離では前のデータと同じくらい離れていたデータでも、マハラノビス距離は小さくなる。

詳しくはこのWebサイトに詳しい。
http://www.geisya.or.jp/~mwm48961/statistics/hanbetu1.htm

マハラノビス距離は以下のように定義される
$C$ :共分散行列
$\bold{x}-\mu$ :値を計算するデータと平均値の差のベクトル
$Maharanobis Distance = \sqrt{(\bold{x}-\mu)^TC(\bold{x}-\mu)}$

共分散行列はpythonのnumpyでもoctaveでもビルトインファンクションで計算してくれるのでかなり楽だ。

Tubdotgagではマハラノビス距離が閾値以上になった場合にその値に応じて罰則を付加している。

機械学習を最適化に組み込むのはとても面白いと感じた。
そして遺伝的アルゴリズムの拡張はいろいろやりやすくて楽しい。

あと某琵琶湖での大会にて、「曲線の美しさとか、生物模倣とか、凄く興味あるんです」と言ってくれた設計の人もいて、凄く嬉しく思った事もかいておかなければならない。

2014-04-20

XGAG ver2 公開！

翼型について

ご無沙汰しておりました。yuukivelです。いかがお過ごしでしょうか。
私自身、研究室に配属され、もう音速とか越えまくってますけどって感じで過ごしています。

さて、遺伝的アルゴリズムによる翼型最適化GUIアプリケーション
「XGAG」

のアップデートができましたので、ここで配布します。

/////////////(4月21日追記ここから)///////////
既定翼型のリンクエラーのバグを改善しました。~~明日32bit版も2.02に更新します。~~※都合により遅れています。コメントにてご指摘下さった事象の改善も含め、でき次第公開します。
64bitOS向け
XGAG_ver2.02 64bitOS

32bitOS向け
XGAG_ver2.01 32bitOS

32bitOSにて上のもので動かなかった場合、こちらを使用してみて下さい。（コンソールが表示されますが最小化して使用して下さい）
XGAG_ver2.01 32bitOS console
/////////////(4月21日追記ここまで)///////////

いつもの事ながらソースコードをGitHubにて公開しています。
https://github.com/NaotoMorita/NFoilDesign

さてさて、今回のアップデートで何が変わったのか紹介します。

端的にまとめますと、
実用レベルまで完成度を高め、GA自体のカスタマイズもできるようにした。
ということです。

具体的に変更点を示します。
まずは遺伝的アルゴリズム関連から
1.遺伝子情報を「Binary表現」から「Gray表現」に変更し、収束性が向上した。
2.MOGA(多目的遺伝的アルゴリズム）によく用いられる「シェアリング手法」を導入し、世代が進んでも幅広く探索できるようになった。
3.毎世代投入するエリート遺伝子の評価関数値を毎世代更新することで、最適化計算実行中も評価関数の値を変えたりできるようになった（あまりオススメはしません）

その他UIや設定に関して
4.翼型の入出力の際にXFOILの補間、翼弦長ノーマライズ、翼型迎角リセットを実行することによって、出力される翼型の速度分布がなめらかになり、Ver1.00で推奨された速度分布平滑化が不必要になった。
5.最適化する各係数の値を表示し、設計に活かせるようになった。
6.最適化する各係数の値の範囲を、後から設定で変えられるようになった。
7.グラフの表示が見やすくなった。
8.追加キャンバーの形状がグラフに表示されるようになった。

大きくは以上8点です。他にも細かいdebugや表示の変更を行っています。

では、具体的に1から順に解説していきますね

1.Gray Cording

Gray Cording（グレイコーディング)ってなんぞやって人が大半だと思います。簡単に説明します。
XGAG ver1では遺伝子表現にbinary(2進数)表現を用いていました。この遺伝子を遺伝的アルゴリズムによって組み換えていきます。
例えば、遺伝的アルゴリズムによってある変数値を10進数にして11から12に変えたいという要請があったとしましょう。
このとき2進数表現では、1011から1100の変更となり、1と0を2箇所入れ替えなければなりません。逆に考えて、ある2進数にて0と1を1箇所入れ替えると10進数が2とか3とか変わることがあります。
これが遺伝的アルゴリズムの収束に悪影響を与えます。これを防ぐのに用いられるのが「Gray Cording」です。
Gray Codeにて0と1を1箇所入れ替えたとき、10進数も1しか変わらないようになっています。
「グレイコード」とかでググるともっと詳しい解説が出てきます。

この実装により、収束性が大分改善されました。
Gray Cording実装自体もpythonの「binstr」というライブラリ
https://pypi.python.org/pypi/binstr/1.3
を見つけてこの中の関数をちょちょっと使うだけなので凄く簡単でした。

2.シェアリング手法

多目的遺伝的アルゴリズム(MOGA)では、パレート解(ある評価関数の値を改善するのに他の評価関数の値を下げなければいけない解)を効率よく探索するために、個体が密集している地域の評価関数の値を罰則として下げる「シェアリング手法」がよく用いられます。
XGAGが採用するGAはMOGAではありませんが、シェアリング手法を導入することで、局所解に陥りにくくできるのではないかと考えました。
実際、単一評価関数遺伝的アルゴリズムにシェアリング手法を導入した例を示しているものもあります。
適応的免疫アルゴリズムを用いた多峰性関数最適化に関する研究
XGAGで用いているシェアリングは表現空間シェアリング手法であり、近い形状をもつ翼型がたくさん集まっている所に罰則を科します。
各最適化係数のユークリッド距離を計算し、シェアリング半径という定数内に入ったものが、その距離に応じた罰則を受けます。
ここで難しい数式を出したくないので、詳しく知りたい方は以下の論文を参照して下さい
多目的分散遺伝的アルゴリズムにおけるシェアリング，収束判定，及び解の評価手法の検討

シェアリングを導入することによって、気持ち一つか二つの解に到達しやすくなったかな？という感じがしています。シェアリングを導入する前は実行するたびに出てくる解が違ったりするので、良い形状を探すのに人力を使っていることもありました（今のバージョンでもそれが必要な場合はありますが）ただし、シェアリング半径は設定できるようにしましたが、この値を大きな値にすると全然収束しなくなるので注意！シェアリング半径を0にするとシェアリングを行わなくできます。

「進化！」と出ているのに評価関数の値が下がったりすることがありますが、それはだいたいシェアリングのせいです。もしくは後に述べる評価関数再計算のせいとかですね。バグや退化ではありません。

3.評価関数再計算

これは（説明しなくても）いいですね。その通りです。あまり推奨されることではありませんが、例えば、翼厚が欲しい値にならないなぁ…とか感じたら、アルゴリズム実行中に評価関数の係数を変えて強制的に解を流すことができるようになりました。
レイノルズ数とか迎角もいじれますが、ちょっと高難度かなと。

4.翼型インプットに対するXFOILの積極的な利用

XFOILには、翼型の曲率が大きい所ほどたくさん点を入れてくれるといった便利な翼型補間機能(PANE)があったり、翼弦長が1じゃない翼型も相似拡大縮小して1にしてくれたり(NORM)、翼型ファイルの段階で迎え角（厳密には表現が違いますが）ついちゃったりしてる翼型の迎え角を0にしてくれたり(DERO)、いろいろ便利な機能があります。これらの機能を積極的に使うことによって、きちんと整理されていない翼型ファイルをXGAGに突っ込んでも大丈夫になりました。（もちろんきちんとしていた方が良いですが）また、この処置を行うことによって、ver1.00では必須な、難易度の高い速度分布調整を行う必要がなくなりました。つまり、XGAGで翼型を作って、そのままXFLR5とかで解析、実際に使用、といった流れが可能になったわけです。
計算の揚力係数が1.200だとか、翼厚11.00%とかの翼型が生成できたら、調整することなくそのまま使うことができるということです！（確認はしてくださいね）

以下はXGAGで作った翼型を速度分布調整することなくXFLR5解析に掛けた結果。圧力分布に不自然な振動がないのがわかります。XGAGとXFLR5ではともにXFOILを使っているので当たりまえっちゃ当たりまえだが、揚力係数の値が1.200（設計目標値）であることにも注目

5.各最適化係数の表示

以下の画像のとおり、最適化する8つのパラメータを表示しました。

これによって、「FX76MP120の割合高いなぁ…」とか「もっとDAE41使えばいいのに…」とか考えられるだけでなく、「この翼型前のと違うぞ！」とか議論することができます。

6.各最適化係数の範囲の設定

以前から要望のあった、最適化係数の範囲を変更できるようにしました。
「設定」タブの「遺伝子係数設定」で変更することができます。

この画面で先に説明したシェアリング半径の値を変えることができます。繰り返しますが、シェアリングを行わない場合はこの値を0として下さい。0.1とか1とかに設定すると、収束性が悪くなります。個人的にはデフォルト値ぐらいで良いんじゃないかなと。研究不足で何とも言えませんが。（実装したのが昨日だった）

7.グラフの表示の改良

具体的には結果の履歴タブで、点を線でつなぎました。これできるようにするのも大変だったんですよ！（汗

8.追加キャンバー表示

先に述べたとおりです。追加キャンバーは、翼型の微調整を行うためのものです。だからもちろんまっすぐだったりして変化してないように見えることもありますが、ちゃんと最適化されています！

いかがでしたでしょうか。いろいろ改善しました…大変でした…（遠い目

XGAGを用いて、もっと凄い飛行機がつくれたら、それは私の本望です。良いお話が聞けることを楽しみにしていますね！

2014-03-11

構造制約を考慮した循環分布最適化GUIアプリケーション「Windmize」

人力飛行機設計

ご無沙汰しておりました。yuukivelです。

今回、また新作ができましたので、この場で公開させてもらおうと思います。
そうです。以前から作る作る言っていた構造制約を考慮した循環分布最適化GUIアプリケーション、その名も
Windmize

ダウンロードはこちら
Windmize (for Windows)
※SmartScreenによって弾かれることがあります。ご注意下さい。
※このソフトを使用した事によるいかなる損害も責任は負いかねます。実際に設計に用いるときは確認の計算を行うことを強く推奨します。

また、Githubにてソースコードを公開しています。
https://github.com/NaotoMorita/Circulation_design

GUIの構築にはQtのpython向けwraperであるpyQt4を用いています。
この記事の後、このアプリケーションを題材にして、pyQt4の紹介及び解説を行うつもりです。
この記事はとりあえず、Windmizeの使い方について記述しようと思います。

また、このアプリケーションについて2014年3月9日に電気通信大学で行われたHPA交流会における発表で宣伝させていただきました。
このときに用いたスライドをslideshareにアップしました。

2014春HPA交流会発表(公開用)

Windmizeでできる事

Windmizeは、このブログで何度も触れてきた、構造制約を考慮した循環分布最適化、通称「TR-797の方法」を人力飛行機向けに改造したTR797改手法を用いています。
それはすなわち、翼端でのたわみの値（最大たわみ）を構造制約条件として、そのたわみを実現する循環分布の中で、誘導抵抗が最小となる循環分布を探索します。
求める循環分布は、インプットされた桁の区切りを基準に、これまた以前記事にした多角形化行列を用いて、多角形の状態で計算されます。これによって、テーパー比を決めるような単純な問題にも、このアプリケーションを使うことができます。
また、構造制約条件を外せば、多角形循環での最小誘導抵抗（すなわち楕円循環分布のようなもの）を求めることができますので、設計として何もできていない状態でのスパンや機速を決定する初期設計ツールとしても、バンバン使うことができます。
計算結果はCSVで出力できるようにしていますので、今まで用いられてきた設計ツールや図面出力ツールに値を入れることも可能だと思います。
その証拠として、エクセルで出力した構造制約を考慮した循環分布としない循環分布のグラフを載せておきますね。

では実際に使い方の方を見ていきましょう。

Windmizeの使い方

Windmize.exeはそれ単体では動作しません（入っているフォルダから出して使用することはできない）ので、デスクトップ等から起動したいときはショートカットを作って下さい。もしくはフォルダ「Windmize」を丸ごと移動して下さい。
Windmize.exeをクリックして、アプリを起動すると下図のような画面になります。デフォルトでセットしてある種々の値は、このブログで触れてきた「紅鮭-肥」の設計数値（実際には微妙に異なる）となります。

とりあえず「計算」を押してみて下さい。計算中を示すプログレスバーが動き、値が出力されれば、アプリケーションは使えています。
時間が許すのであれば、デフォルト設定を少し変えて、どのように値が動くのか遊んでみると良いと思います。ある程度設計に関わっている人であれば、何をやっているのかすぐにわかると思います。

では、詳細な設定の仕方を説明しますね。

上のグラフとその下の行は設計結果の表示のための領域です。計算が終わったらここに結果が表示されます。
構造制約係数とか揚力制約係数の値を見て設計できる人とか神様だと思うけれども、一応表示しておきました。

その下が計算の設計変数を設定する部分です。まずは最適化変数の設定の部分の説明をしますね。

最適化変数の設定

揚力[kgf]

そのままです。必要な揚力をkgfで入力します。

速度[m/s]

飛行機が飛行する速度を入力して下さい。Dist機であれば6.5〜8.00[m/s]、TT機であれば8.5〜16[m/s]くらいだと思います。

最大たわみ[mm]

最外翼の終端が翼根と比較してどれだけ上に行っているか（最大たわみ値）をmmで入力します。このプログラムは「ある桁があって、その桁で指定した最大たわみ値を実現する循環分布をもとめる」プログラムですので、構造制約にとって最も重要な値となります。

ワイヤー取付位置

ワイヤーに代表される「翼を下に引っ張るもの」を取り付けるスパン方向の位置です。ワイヤーじゃなくても例えば双発機のプロペラとかもここで考慮できます。

ワイヤー下向引張

ワイヤーに代表される「翼を下に引っ張るもの」の下向きに引っ張る力を単位[N]で入力して下さい。注意すべきは「ワイヤー張力」ではないということです。この値を0にすると、一般的な片持ち機が設計できます。

dy[mm]

計算に用いられる翼素の(y軸に投影した）幅です。この値を大きくすれば計算は速くなりますし、小さくすれば精度は良くなります。しかし、精度が良くなると言っても限界はありますので、一般の人力機のスパンであれば50[mm]でいいと思います。なんでこの値を変更できるようにしたのかというと、ラジコン機の設計に使うときに50mmだと大きすぎるからです。
※dyを大きくしすぎると落ちます。あたりまえっちゃあたりまえですが

ここまでが設計変数を設定する部分です。続いて桁に関する設定に移ります。

桁に関する設定

桁に関する設定でできる事は3つ。

各翼の終端を入力すること（つまり中央翼はここまで、内翼はここまで、みたいな入力の仕方をする)

翼の数を増やしたいときは「列追加」減らしたいときは「列削除」を押して下さい。ただし、これによって挿入された桁の剛性・線密度は必ず固定の値が初期値として入りますので注意して下さい。（例えば第5翼を削除して、列を増やして第5翼を復活させても、設定されていた剛性・線密度の値は復活しない）
また、もちろん内側の翼より外側の翼が翼根側に来てても上手く計算できません。現在のバージョンでは、そのまま落ちます。

詳細に桁剛性・線密度を入力すること。

これは桁詳細設定内でできます。桁詳細設定を開くと下のようなダイアログが出てきます。

このダイアログでは、桁を4分割して、それぞれに対して桁剛性と線密度を設定できます。区切りの考え方は、その桁の翼根側を0[mm]とした座標系で考えます。これで、俗に言う「階段積層」が（一部制限はありますが）考慮できます。各翼に対する入力は上部のタブで管理されています。
桁の剛性や線密度をいちいち入力するのは面倒くさいですね。このダイアログではコピー&ペーストに対応していますので、メモ帳などにメモしてお使い下さい。え？セーブ機能を実装しろって？それは次くらいのアップデートで実装します（汗

調整係数によって、各翼の剛性と線密度を大まかに調整すること

先に書いたようにいちいち剛性や線密度を入力することは、設計の仕上げにやることとして非常に重要なのですが、設計の最中に細かい設定をするのは非常に面倒な作業になります。従って、この調整係数を用いて、各翼の剛性と線密度を調整できます。単純に、剛性と線密度にこの値を掛けているだけです。もちろん剛性と線密度は比例しなかったり、積層するような構造部材では実際にはとびとびの値になりますので、調整係数で大まかに値の見当をつけて、いいものができたと思ったらその剛性値を狙って桁を設計し、できたら実際の剛性・線密度をWindmizeに入力してやってみる、といった流れが良いでしょう。

ここまでが設定の説明です。

使い方のコツ

このアプリケーションは言うなれば「揚力を中央に寄せるためのアプリケーション」です。従って、

このような明らかに楕円循環（緑点線）より中央に揚力が寄っている循環分布なら良いのですが、

このように,外側によってしまった循環は「もったいない循環分布設計」となります。この計算手法ではたわみ角による揚力の損失も考慮した循環分布を求めますので、「たわみによる揚力の減少を考慮して、すこし外側に揚力を寄せる」必要はありません。

また、このような、翼端での循環が負になるような循環分布も不自然です。

またWindmizeに内包されている材料計算は、空力最適化のためのそこまで厳密ではない材料計算です。それなりに実際と計算のズレが出ることはご了承下さい。

値のCSV出力について

CSV出力ボタンを押せば計算結果がCSVに出力されます。
これは普通ですが、中身に少し注意してもらいたい点があります。
だーっと並んでいる循環の値やたわみの値についてです。
この手法では、各翼素の値は翼素の中心の値で計算しています。従って、各翼終端での値は、だーっと並んでいる、各翼中心の値には含まれていません。
CSVの中身の方にも注意がありますが、循環値だけは各翼終端での値がCSVに記されています。
それだけ、注意して下さい（CSV開いてみれば何のことかすぐわかると思います）

最後に

構造制約を考慮した循環分布最適化によって、循環分布を設計する時代がやってきたなと感じます。私が人力飛行機の設計を始めた頃には楕円循環分布を用いる、もしくは感覚に頼るしかなく、ある意味盲目的な所もありました。この「Windmize」によって、楕円循環分布を用いるにしても、きちんとした理由が存在する、そういった時代になってほしいなと思っています。
また、このツールをネタに用いたpyQt4の解説も少し行います。pyQt4は日本語の解説が少なくて、コードを書くときに苦労しました。Windmizeからアプリケーション作成に興味を持ってくれて、新しいことをアプリにまとめて、共有する人がもっとたくさん出て欲しいです。

それでは、みなさん「Windmize」をよろしくお願いします。

2014-01-16

ロケットシミュレータとQuaternion演算

飛翔体運動

こんにちは。yuukivelです。今回は飛行機の話題からちょっと離れて、飛翔体（ロケットとか）のシミュレータを作った（作っている）ので、この公開と、この中で用いた座標変換処理におけるQuaternion（クオータニオン：四元体）について解説する。

まあとりあえず先にダウンロードリンクを載せておこう。Octave スクリプトで書いているが、MATLABでも動くし、PythonのGUI版も作りかけだったりする。
（使用して損害が生じたとしても責任は負いかねますのでご了承下さい。）
mq_rocket zipファイルダウンロード
 githubはこちら
https://github.com/NaotoMorita/mq_rocket/

パラメータをいじる場合はrocket_parameter.m内の変数の値を変える。

本当はPythonでGUIを積んだアプリケーションの形で公開するつもりだった。そうでないのは、言い訳をすると、

Numpyでの演算で生じる謎の微小な周期的数値振動
グラフ表示のライブラリであるMatplotlibのスパゲッティーコードとイケてNASA。３次元プロットの軸まわりの処理が残念。

の二つ。克服したら公開します。最終的にはOpenGLやDirectX等を用いて３次元のモデル表示まで持っていく予定。

さて、今回はこのスクリプトの中で用いた四元数、すなわちQuaternionについて紹介する。

まずあなたが三軸のジャイロセンサーを手に入れて、値も綺麗に取れたとしよう。そうしたら何をするだろうか？
当然、角度を求めようと積分するだろう。

ではその積分はどうする？X軸、Y軸、Z軸で積分するだろうか。当たり前だと思う人もいるだろう。
否、それは大間違いなのである。

ジャイロが存在する座標系は「動座標系」であり、私たちが知りたいのは、ある地面から固定された座標から見たジャイロセンサーの傾きの角度、すなわち「オイラー角」なのである。

オイラー角は固定座標系、すなわち「慣性座標系」からみた動座標系の傾きを示すものであり、それ自体が回転を定義する。
このため、オイラー角から、３次元座標の回転を表すことができる。これをDCM(Direction Cosine Matrix)という。

さて、このオイラー角と回転行列を用いて飛翔体の運動方程式を記述した際の長所と短所をまとめておこう。

オイラー角で運動方程式を記述した場合

長所

わかりやすい。定義も自然なので、おかしければすぐに分かる。後処理も簡単。

短所

三角関数、逆三角関数を多用するので、計算が重くなりやすい。これはマシンパワーが無いマイコンとかで積分するときには致命的
そもそも見かけの力を定義する式が煩雑でわかりにくい。
ジンバルロックを起こす。

DCMで運動方程式を記述した場合

長所

計算が楽

短所

姿勢更新を行うことができないので、DCM単体では「そもそも運動方程式が記述できない」（！）
単なる行列なので、見ても何が何だか分からない。

ここで登場するのが、QuaternionとDCMを組み合わせた運動方程式の記述なのだ。

QuaternionとDCMを用いて運動方程式を記述する

Quaternionを用いると以下の事ができる。

代数演算だけで回転行列が記述できる。
角速度を用いたQuaternion更新の公式がある。
逆Quaternionを用いれば逆回転行列を求めるのも簡単

こうした特徴の御陰で、軽く、精度の良いシミュレータを実現できる。飛行機の運動方程式も、XFOIL.mとの組み合わせてやれば、微小擾乱といった近似を外して軌跡を求めたりできる。（というより、最適化理論と組み合わせて設計に使えるものをいま作ろうかなとか思ってる。話変わるど、グラディウスを遺伝的アルゴリズムで解く（？）とかすごいよね）

以下にクォータニオンを用いてロケットの運動方程式を書く際のフローチャートを載せる。そのまんまmq_rocketのフローチャートである。

この中で出てくる速度座標系とは機体座標上での速度のX方向にX軸を合わせて設定する座標系。揚力や抗力はこの座標系の軸方向に働くので、拘束で高精度な力見積もりが可能だ。（ただしCLの求め方とかざっくりだけど。）

CLって言葉を使って思い出した。僕がコード内で使っている言葉は基本的に飛行機のやつを使っています。早くなれないと…

詳しいQuaternionの用い方は以下の論文を参考にして欲しい
クォータニオン計算便利ノート
またオイラー角とクォータニオンの比較を行っている論文として
TM-636 クオータニオンとオイラー角によるキネマティックス表現の比較について
この二つをオススメする。

他にも航空宇宙工学便覧等も参考になる。

航空宇宙工学便覧

作者: 日本航空宇宙学会
出版社/メーカー: 丸善
発売日: 2005/12/01
メディア: 大型本
この商品を含むブログを見る

Quaternionが扱えるようになると、運動方程式の記述がグッと楽になる。Quaternionは3Dモデルの座標変換にも使われているようだし、これからステップアップしていくために必須のツールだと思う。

さて、テスト乗り越えたら、またイロイロ手を出すぞー（あ、テスト終わる前もいろいろやってるか（笑））

2014-01-03

夕暮れ時の街を飛行ロボットで空撮してみた。

こんにちは。yuukivelです。

今回は人力飛行機を離れて、ラジコン飛行機の話題をば。

手っ取り早く言うと、ラジコン飛行機を用いて空撮してみた。たんなる動画の宣伝記事です、はい。

まず最初に動画をどうぞ。

こちらはアクションカムについてるGPSで機体の速度や位置を表示している。

元動画はこちら

メチャクチャいいっすね。BGMかけたりしながら見てると全然飽きないッす。
昔遊んだ公園とか、上空から見られるようになるとは…感動。

夕焼けも綺麗だ…

用いた機体

飛行ロボット「Daven」中ペラ複葉無尾翼の機体。

この機体について話し始めると長くなるので略。
とある大会に出した機体を改造した。

飛ばしてる様子はこちら。

この機体、機体形状の割に安定性が高くて、かつ風にも結構強い。足に問題点があったりするのだけれど、大会の時には自動操縦を入れたり、なかなかブンブン飛び回ってくれる機体だ。
設計速度は3m/s、200gなので、かなりゆっくり飛ぶ。ペイロードも結構ある。

用いたカメラ

SONY　HDR-AS30V
http://www.sony.jp/actioncam/products/HDR-AS30V/
ソニーのアクションカム

年始のヨドバシの福袋、（たぶん)GoProを無視して、このHDR-AS30Vの約90gという軽さ、GPS搭載ということで買った。
最初に紹介した動画みたいに、簡単にGPSの情報が表示できてすごくいい。

私の趣味だったりするスキーの時にも重宝しそうだ（笑）付属ケースに入れれば防水だし、耐衝撃、防塵仕様。
どんな動画が撮れてるのかWi-Fiで飛ばしてスマホで見れたりする。スマホがリモコン代わりになる。
画質もいい。いいもの買ったなぁとしみじみ。

あと、人力飛行機とかでも、動画撮りながらGPS取れるし、まさに夢ひろがりんぐ。

2014-01-02

人力飛行機を割と本気で設計してみた。第二報

人力飛行機設計

さて、第二報です。今回は翼型設計、というより自分で自分のアプリケーション
「XGAG」
XGAG　遺伝的アルゴリズムによる翼型設計GUIアプリケーション
をつかって翼型を設計してみた結果およびノウハウ等を書きたいと思う。

実際に設計した翼型は第一報の方のXFLR5に入ってます。加えて翼形状の下にダウンロードリンクを張っておきました。

まずは翼型設計コンセプト。二つあげられる。

薄翼にして、とにかく低抵抗であること
揚力係数が循環分布設計によって求められたものとできる限り一致すること

である。
最近はどこかの機体の影響もあってか、翼根の翼厚が13%を越えるようなぶっとい翼型が流行っているけれども、今回はワイヤー機ということで、できるだけ翼厚を小さくして、形状抗力を小さくしようという考えが根底にある。また、せっかく循環分布きちんと設計したんだから、できるだけそれに合わせようと頑張っている。

紅鮭・肥の主翼はセクション毎に翼型を設計して、それを線形に混合するという形を取っている。こうするとプランクとか線形に変化していくことになって作りやすいし再現性上がるのかなとか考えている。

翼根翼型「BENIFOIL_1」

BENIFOIL_1.dat
本当はもっと薄くして形状抵抗を減らすつもりだった。けれども9%程度の翼厚では、設計点を外れたとたんに性能ががた落ちしてしまう上に、見た目もガタガタのものが多く、さらに下面のキャンバーがきつすぎて変な形になってしまうものが多かった。変な形って何だよって感じだが、違和感と言えばいいのかな？やっぱし綺麗な形使いたいからね。

翼厚を11%程度にすることで安定した性能と、作りやすそうな形になった。前縁はDAE31、後縁はFX76MP120ぽくなった。性能は下のポーラーカーブ参照
Re:440000 0-10° XFOIL解析値

とにかく翼根で高いCLが要求されたことから、着目すべき点は「失速特性」および「高CLでの抵抗」である。翼根でCm気にするのは野暮ってもんです。そのためには上面の圧力係数を前縁にて急激に引き下げる事が有効だ。この手法が採られているのがDAE31やFX76MPシリーズであるからこれをベース翼型に使ってやる。この条件設定によって相対的にれらの翼型が際立ってくるので、それを修正するような翼型を他のベースにすることで、既存翼型の改善を狙う。。。という感じか。ちょっと良さそうなのができたらそれを保存してまたその翼型をベースに。。。といった形で欲しいCLを実現した。速度分布修正はXFLR5で、速度分布の高周波成分を取り除くhanning filterをかけた。この翼型は翼根向けなので、取付角を微調整して、厳密にCLを合わせる感じ。

性能としては低CL域ではDAEにほんの少し劣るものの、使用CL域で抵抗はDAE31より小さい。失速特性もDAEより良いと言えるだろう。欠点は後縁形状で、少し細く作りにくいかなと思うがFX76MP120より全然マシ。

この翼型を中央翼、内翼に混合無しで用いる。内翼は捻り下げであり、CL1.125付近まで捻り下げる。

外翼翼型その1「BENIFOIL_2」

BENIFOIL_2.dat
薄くして形状抵抗を落とした翼型。設計CLを1.012とすこし低く設定した。そこでこのCL域でのの性能確保（低抵抗）を一番の目的とした。

うーん、この翼型はあんまり言うことないなぁ。。。（笑）
ポーラーカーブはこちら。高CL域ではDAE31の圧勝だが、低CLでは低い抗力である。御陰で最大揚抗比はDAE31を凌駕していたりする。
Re:390000 0-10°XFOIL解析値

着目すべきは下面の失速による性能低下がBENIFOIL_2の方が小さいという点だ。高CLでおおきいキャンバをもつ翼型は下面が乱流遷移する事がある。低CLを実現するときは、やはり大きいキャンバの翼型を捻り下げて使うより、小さいキャンバの翼型を設計してしまった方がいいなと強く感じる今日この頃。

外翼翼型その2「BENIFOIL_3」

BENIFOIL_3.dat
BENIFOIL_2とほとんど同じ形状だが、微妙に太い。設計点は迎角1.48°、CL0.832。BENIFOIL_2との形状の類似は偶然。（ベース翼型にBENIFOIL_2が入っているので完全に偶然とは言えない）
しかしこれは
外翼その2を取り外して飛行することができるこの設計
にとってかなりプラスである。

いつもどおりポーラカーブを載せておく。
Re:390000 0-10°XFOIL解析値

ちょうどDAE31だと下面の乱流遷移が始まるところだし、DAE41にすると使用CLが高すぎて性能がでない。
また、BENIFOIL_2と比べると少し失速特性が穏やかになっていて、レイノルズ数が急激に下がり始める最外翼での使用に適していると思われる。

XGAGで作るときはBENIFOIL_2をベースに修正を加えていった感じだった。

最外翼翼型「BENIFOIL_4」

BENIFOIL_4.dat
構造制約を考慮した最適循環分布計算では、その狙いの通り、翼端で発生する揚力は小さくなる。揚力は翼弦長×揚力係数CLに比例するから、翼弦長を大きくしすぎると設計CLが低くなって取付角マイナスもしくはキャンバーが下にそるとか残念な事が起こる。逆に翼弦長を小さくしすぎると、レイノルズ数が小さくなるうえに揚力係数が高くなるからどんどん失速に近づいて。。。といったジレンマが発生しやすくなる。BENIFOIL4はその中間を狙った感じ。キャンバーはほとんどなし。翼厚も太いが太すぎず。。。な値を選定した。

ポーラーカーブでDAE41と比較してみよう。
Re:220000 -4-4°XFOIL解析値

DAE41はまさに楕円循環分布の翼端のために作られた翼型のようなもので、楕円循環だとちょうどいいCLだし抵抗は小さいしで、その存在感の薄さの割にいい翼型。
この翼型に対し、CL0.5付近から抗力の増大が見られるものの、それ以下での抵抗はDAE41より小さくなっている。

ただしこの翼型は、1.48°の設計迎え角で使用すると、下面に乱流遷移が見られる。これが翼端で剛性の低い桁と相まってちょっとだけ不安。

まとめとおまけ

第二報では翼型に焦点をあてて、その設計例やノウハウを紹介してみた。
総じて言えることは（自分で言うのも何だけど）
マジXGAG作って良かった！
ということ。今までみたいに翼型に縛られた桁設計をしなくていいし、なにより、だいぶ主翼の性能を上げることができる。
また、外翼その2を取り外して翼幅30m機にするというのも、XGAG無ければ設計難しいだろうなーと思う。

話は外れるけれども、やはり翼幅30mと34mでどれだけ違うのか、解析では分かるけれども実際どうなのか、やっぱり気になりませんか？
34mでの運用はパイロットが慣れるまで結構大変だと思うし、30mの方が風が多少強くても飛ばせそうだ。

こんな風に、今までなかった機体の、時代を半歩ぐらい進めてくれる機体の実現をめざして、ちまちま活動をはじめているところだ。

私の大好きな曲の歌詞をもじってとりあえず第一報、第二報の締めくくりとすることにしよう。

世界分かつ空に耳を澄ませば、招く声が響く。
旧知の軛と袂を分かち、そして生まれる
新世界