失速特性を考慮した翼型設計ocXgag その2 公開と使い方

さっそく、公開および使い方についてみていこう。

公開

PyInstallerを用いて、Windows 64bit版のPCを用いてbuildした。
twitterのフォロワーさんの協力によってWindows 32bit版でも動作することを確認している。

いつも通りDropBoxの公開リンクを用いた公開です。
https://dl.dropboxusercontent.com/u/36164244/ocXgag_win64build.zip

ocXgag_win64buildには、前述のとおりbuildしたocXgag.exeおよび解析を行うxfoil.exe、および翼型のサンプルとソースコードが入っている。
pythonをいじれる人は、ソースコードから起動したほうが早いかもしれない。
今回のbuildはすべてexeにまとめたので、xfoil.exeとocXgag.exeが同じディレクトリにあれば、どの場所にあろうと動作する。
また、consoleは出さないようにしたかったのだけれども、なぜか--noconsoleオプションをつけてbuildするとXFOILが起動しない。
申し訳ないけれどもconsoleがいちいちうるさいと感じる人はウィンドウの最小化をしてほしい。一応出しておけば、ちゃんと動作しているかどうかの確認になる。

mac版,Linux版は作成できる環境にないので、私の手では作れなさそう。どなたか協力してくれる方がいるとありがたい。

XFOILを使っている関係上、このプログラムのライセンスはGPLなのでソースコード公開の義務がある。
今回のbuildのソースコードはzipファイルに同梱した他、githubにてアップデートを行っていく予定なので、そちらで追ってもらえるとよいかと。
(現状githubへのpublishにエラーが出ているので、でき次第この記事を更新します)

使い方

このプログラムは、失速特性を維持しつつ、抗力最小を狙って翼型を「修正」するものだ。
したがって揚力係数や翼厚の大幅な変更は難しいし、そういった翼型の「作成」はXGAGを使ってほしい。
http://d.hatena.ne.jp/yuukivel/20140420/1397998798

exeから起動すると、pythonやら様々なモジュールを解凍する作業が入るため、起動にそこそこ時間がかかる。
起動すると次のような画面が現れる。基本的な手順は図の中に示した通り。

エラーが出る場合は大体(というかほとんど)XFOILがうまく回っていないことに起因する。
XFOIL用の一時翼型ファイルを書き出し⇒XFOILに読ませる⇒そのファイルを削除
といった感じでXFOILを回しているのだけれど、時々Permission deniedによってXFOILがうまく回る解析条件でもエラーが出ることがある。
エラーのダイアログが出たら、とりあえずもう一回翼型のオープンや計算実行をしてもらうか、それでもダメな場合は解析条件を変えてほしい。
もちろん解析する翼型が存在しなかったりしたら、エラーダイアログが出るのは当たり前だが。

続いて、Open Foilをクリックして翼型を開くとXFOILの解析が同時に行われ、各種性能と最適化の設定値のとりあえずの値がセットされる。
また、Setボタンの横にあるNorm and Dero from fileというチェックボックスについて。
ここにチェックがあると、翼型を開く際に、XFOILによる翼型点数の変更(n=200)と、翼型の翼弦長を1とするNormおよび前縁を原点にしつつ翼弦線をx軸と平行にするDeroが行われる。
初めて開く際はやるといいのだけれど、これをやるかやらないかによって微妙にCL等が変わるのでお好みで。

ここで最適化の設定値について説明する。

• CL target：揚力係数CLの目標値。下のtolと合わせて(1-tol)CL target<=CL<=(1+tol)CLという制約条件。
• Thn(%)target：翼厚(%)の目標値。この設定値以上の翼厚という制約条件。
• Thn Calc Point(%)：翼厚の計算位置(%)。どこで翼厚の値を計算するかということ
• TEang(deg)：後縁開き角(deg)の目標値。この値以上になるようにという制約条件。後縁開き角とは。上面および下面の後縁角の差のこと。
• xt upper(%), xt lower(%)：修正PARSECにおける最大翼厚位置。詳しくは後ほど
• dXtr tolerance(%)：Hの制約値を作る際に、元となるHの分布をどれだけ左右に動かすかというパラメータ。小さくすればそれだけHの分布は似てくるが、なかなか収束解を見つけるのが難しくなる。
• tol(% to CL tgt)：上で説明したとおり。この値はCLの許容誤差がCL targetに対して何%かということを表している。
• max Iter：解修正動作の最大繰り返し数。

xt upper,xt lowerについて
このプログラムは翼型のY座標を以下のようにして表し、繰り返し修正する。

ここでからを決定するために、

に対して以下のような条件を用いる。
1.翼型の前縁を通る条件
2.翼型の後縁を通る条件
3.xtにてあるY座標をとおる条件
4.xtにてY座標の微分がゼロとなる条件
5.xtにてY座標の2階微分値がある値となる条件
6.前縁半径がある値となる条件
7.後縁角がある値となる条件
これらのパラメータを範囲を絞って最適化していくことで、修正の方向を逐次的に得ているのだけれども、xtの値だけはこちらで指定してあげる必要がある。
なぜならxtだけ動かしてほかの値をゼロとして評価関数の勾配を取ろうとしても、修正される翼型に変化はないので、必ず勾配が0になってしまうから。
無理やり勾配を取ろうとしたけれども、結局収束性が悪化したので、固定値とした。
また、Local Aerodynamic Approximationにおける勾配を取る際にXFOILの出力桁数が問題となってしまうため、勾配を取る際のパラメータの幅は今後の課題となっている。

さて、計算を実行すると次のような画面になる。

Res to Orgのボタンで計算結果をOriginalFoilにセットしなおし、計算が続けられるようになっている。
まだ計算のインタラプトは実装していないので、maxIterの値を大きくしすぎると、解が振動したり、XFOILの解析の限界にぶち当たってうまく回らなくなったりすることがあるので注意。
個人的にはmaxIterを10回くらいにして様子を見ながら計算を進めていくのがいいなと思っている。

Local Aerodynamic Approximationは研究段階であり、まだそこまで精度がよいわけではないので、各繰り返しによって制約条件を満たしていたり満たしていなかったりする。
また収束(制約条件を満たしつつ設計変数の修正量ベクトルのノルムが10e-4以下)になることもほとんどない。
Feasibleとは、最適化における専門用語の１つで、制約条件を満たしている状態(実行可能)
Infeasibleとは、制約条件を満たしていない状態(実行不可能)
したがって、実際には収束していなくても結果を見て、Feasibleでかつこれいいなと思ったものをエクスポートして使ってもらう感じになる。

計算実行後は次の図の通り。

Update Hconst at Res to Org　というのは、計算結果の翼型のHの分布を使って、Hの制約条件をアップデートしますかということ。
チェックしなければ最初の翼型を基にしたHの制約条件のままとなり、基本的にはmaxIterを増やしただけでそのまま計算を続けることと同義になる。
チェックすればHconstの値がアップデートされるので、収束解を得やすくなる。
特にCLを増やす要求だったり、翼厚を大きくする要求のときはHの値が大きくなりやすいので、多少のHの分布を許容するのであれば、これにチェックをつけるとよい。
また、ResultFoilの横にNorm of Direction = 0と出たときは制約条件を満たす解が存在する見込みがないといわれていてどんづまり状態であるので、これにチェックをつけると計算が続けられる可能性がある。

おわりに

なにか不明な点があれば、コメント欄に記入もしくは私のtwitterやらfacebookやらに伝えていただけると助かります。

さて、失速特性という翼型設計にとってなかなか難しい問題をなんとか自分なりに定式化して解いてみようという試みなのだが、いかがだろうか。
この手法をもちいて設計した機体が、良い機体だと評価されればソフトを作った身としてはうれしいなぁ。
もちろん、自分で設計した飛行機がよい評価を受けなければならないのはそうなんだけど、やっぱり新しい手法を用いて飛行機をつくれるのは楽しいね。

失速特性を考慮した翼型設計ocXgag その1 背景

とてもとてもご無沙汰しておりました。Yuukivelです。

私事ですが、今いる大学の博士課程への進学を決めました。
これから3年間、研究に(趣味の)研究に精進して行けたらと思います。

また、2016年の琵琶湖の某大会において、私たちのチームFlightWorksは1029.75mと少々残念な結果に終わってしまいました。

ザッパーン

で、本記事はその結果の考察をうけて作った、新たな翼型設計プログラムの紹介と公開記事となります。

FWフライト考察の抜粋

さて、FWのフライトにおいて、大きく二つの原因によって飛行時の必要出力が増大し、短距離にて着水したと考える。

1.機体の速度安定が非常に小さいもしくは不安定だったため、飛行の大部分をバックサイド側にて行っていた。
2.翼端の翼型の失速特性が悪く、速度不安定によって減速しがちな機体において、頻繁に層流剥離を起こしていた。

以上の2点。

1はこの記事の本題ではないので詳しく述べないが、速度不安定というのは、加速をすれば抗力が減ってさらに加速、減速をすれば抗力が増えてさらに減速するという状態。
水平飛行の飛行機において、誘導抗力と形状抗力が等しい速度で飛べば抗力は最小となる。
以前の記事に書いた通り、この機体はスパンを固定して最適設計を行ったわけだから、基本的にはこの底の抗力を目指して設計されることになる。

すなわちどういうことか。機体の抗力が小さすぎるということだ。より正確に言えば形状抗力が小さすぎる、翼が薄すぎる、ということ。
まさに狙った通りの設計だったことが証明されたわけだけれども、これによってバックサイド側の飛行出力とフロントサイド側の飛行出力が接近し、圧縮機でいうところのサージングのような現象が起こったというわけだ。

2つめ。
以下のグラフを見てほしい。

このグラフは、この機体の翼端の翼型を、迎角0.5°、レイノルズ数25万で解析した際の、形状係数H(Shape Factor)の値である。
形状係数はこの後詳しく解説するが、この値が層流において3.55(西山哲男「翼型学」より)、乱流においては2.2から2.4(http://adg.stanford.edu/aa200b/blayers/turbseparation.html)を超えると剥離が起こるといわれている。
形状係数が大きければ、基本的には摩擦抗力係数は小さくなるが、その分、剥離に近い状態になる。
この翼型はHの履歴を見ると、上面、下面ともに早い段階でHの値を引き上げてさらに増加させることで、抵抗を小さくしようとしていると考えられる。
しかしながら早い段階で層流剥離の基準値である3.55を超え、さらに下面も高い形状係数値になっている。
以下にAG24の同条件での形状係数を示す。

AG24の形状係数分布は比較的一般的な翼型のものと一致する。
翼型は遷移点付近で急激にHの値を引き上げて、できる限り剥離しないように形状係数を設計することが多いようだ。

すなわち何が言いたいかというと、FWの機体の翼端の翼型は、性能を追求するあまり非常に失速特性が悪い。
速度不安定によって機体が急激に減速すると、失速特性の悪い翼端が層流剥離を起こし、ストンとロールが始まる。

FWの記録は、ほとんどこれが原因だろう。
翼端の翼型の形状が普通だった、そして翼根の翼型の特異性を気にかけすぎたことによって、設計においてこの点を見逃していた。

そこで、失速特性を考慮した翼型設計になるわけですよ。

ocXgag --Optimal Control of Viscous/Invicid Interaction for General Airfoil Generation

公開と使い方についてはその2に譲るとして、基本的な考え方を紹介しよう。

名前の由来は前に作ったXGAGにoptimal Controlを付け足した感じ。Viscous/Invicid InteractionがXFOILにかかっているので、X。

まず翼型設計コンセプトは以下。
「境界層成長に影響を与える形状係数Hの値を失速特性のわかっている翼型と同様に保ちつつ、揚力・翼厚の調整および抗力の最小化を行う」

境界層方程式は、境界層運動量厚Θを状態量に、以下の微分方程式で記述される。

ここでは翼型の形状(と境界層厚さ)によって一意に決まってしまうため、多少変わるのは仕方ない。
他方、の値は(翼型学の本によると)とHの値によって決まるため、境界層の成長の仕方は形状係数Hに大きく依存する。

ここで形状係数HまたはShape Factor(Shape ParameterとかShape Functionとか言ったりもする)について。
Shape Factor は境界層排除厚と境界層運動量厚の比である。



この値は境界層内速度分布が立つ、すなわち剥離に近くなれば大きくなり、壁付近にて速度分布の傾きが大きい剥離しにくい速度分布であれば小さくなる。
また先述したように境界層方程式の運動量厚の成長に大きな影響を持つ。

そして「境界層が似た状態にあれば失速特性も似てくるだろう」と安直な考えのもと最適化問題を定式化すると以下のようになる

傾向として後縁開き角が小さくなりがちだったので、これも制約条件に加えた。

Hconstの値の設定の仕方については、オリジナル翼型のHの値をdXtrだけ左右に動かしたものの最大値を制約条件の値とした。
これによって、遷移点の前後によるHの変化を吸収することができる。

この非常に制約条件の多い最適化問題を効率的に解くため、遺伝的アルゴリズムのようなメタヒューリスティクスではなく、勾配情報を利用するSLSQPという方法を用いている。
また、勾配評価ごとにXFOILを回していたのでは有限時間で計算が終わらないため、ある点周りの境界層情報を近似して設計変数の変化に対する各パネルの圧力分布や形状係数値を表し、その積分によって揚力係数や抗力係数を求める「Local Aerodynamic Approximation」法を用いている。この方法はまだ研究段階のため、収束に難があるけれども、計算に時間のかかる流体と関数評価回数の多い最適制御問題を橋渡しするものとして注目している。

翼型の表現については、ベースとなる翼型に修正PARSEC法による微小振動を加える形で局所的に翼型を修正し、その繰り返しで解を探索していく。
正直Optimal Controlか？といわれればそうなのだけれども、一応シューティング法っぽい手法を使っているということで許してほしい(笑)
本当であれば、XFOILを使わず境界層計算も自前で実装して境界層の状態量と形状を同時に解くDirect Collocation法で実装したかった。

実際にDAE31と、この手法を使ってCLを上げ(CL=1.12)、翼厚を薄くした翼型(9%以上)の揚力係数-迎角のグラフを比較してみる。

二段階失速の開始点はほぼ変わらず、全体的にオリジナルDAE31を上にシフトしたような翼型となった。
解析点の抗力係数については、オリジナルDAEが80.4カウントなのに対し、修正後が76.6カウントとなり、薄翼化による抗力係数の削減が実現できている。
すなわち、失速特性を維持しつつ、桁や指定の循環値に合わせた翼型の調整が可能となる。

さて、ではその2にて、配布および実際の使い方を見ていく。

"紅鮭"三面図とXFLR5

やっと少し涼しくなってまいりました。いかがお過ごしでしょうか。

さて、今回は前回の記事でまだまだ公開しないといっていた三面図とXFLR5で作成したモデルを公開しようと思う。
XFLR5はver6.10からバイナリファイルの拡張子が.wpaから.xflに変更されたので注意。

まずはXFLR5から。

benishake.xfl

続いて三面図

FlightWorks_benishake.pdf

安定性設計とXFLR5

今回、尾翼を設計するのにXFLR5をフルに活用した。
XFLR5の良さは、モデルを作ってしまえばそこそこの精度で「高速に」結果を出してくれることにある。
特に安定性解析は優秀で、推算の面倒な有次元安定微係数をそれこそパッと出してくれて、固有振動数や減衰率を出してくれる。

たとえば慣性能率や有次元安定微係数を見たいなら、[Stability Analysis]を実行した後、[View Log File]でログファイルを開けばそこに入っている。
こうしたログファイルをtxtとして保存すれば、プログラムから参照してさまざまな計算にもちいることができる。

今回尾翼設計の指針にしたのは、おなじみの尾翼容積の値(動ファクタは正直見てもわからない)と、縦短周期モードの挙動、そして横ダッチロールモードとスパイラルモードの挙動。
またunityを使ってStability AnalysisのログファイルとCLCDCM-alphaのグラフを読み込んでシミュレータを作成し、実際の操縦感覚を見ながら設計を行った。
後者については、設計⇒シミュレーションのサイクルをかなり高速に実行できるようになったので相当便利。まさにXFLR5様様。

XFLR5の慣性能率(慣性モーメント)については、設定でpropeller等の質点を追加していくことで、そこそこ推算できる。一応3DCADのモデルと照らしてみたけれども、オーダーはあっている。その他制度検証についてはXFLR5のドキュメンテーションに乗っている。ということで今回はこの結果を信頼するということで、、、

水平尾翼の大きさは
①容積比からざっと計算
②解析をして縦短周期固有振動数と減衰率を見る(ダイアログに出るFNが固有振動数,ζが減衰率)
③航空機力学入門にのっているグラフ(加藤寛一郎,航空機力学入門,p183,図8.2)の領域の中でよさそうなところに落ち着いているか確認
の手順で行った。その結果ほとんどピッチがオーバーシュートしない、びっちりピッチが決まる機体となった。
ただし、テールの長さがどんどん伸びてしまい、Distの中でも最大級の大きさに、、、(某二人乗りは除く)

垂直尾翼については
①それっぽい形を作る
②安定解析
③スパイラルモード基準を満たす(振幅倍増時間が20秒以上)まで面積を調整する
の順で大まかな形をつくり、シミュレーターを使って微調整した。
その結果、人力飛行機に珍しくスパイラルモードが安定となり、ダッチロールも振動的ではなくなった。
これはDist機にしては少し大きい主翼上反量と、長いテールが影響しているのかなと考えている。

また両尾翼ともテーパーのきつい形状になったのは、揚力を発生した際に少しでも楕円循環に近づけて抵抗を減らすため。またテールが長いのも面積を減らして抵抗を減らすため。
飛行機設計の王道は尾翼に強いテーパーをつけないで翼端失速を防ぐのだけれども、今回は十分にレイノルズ数が大きくまず失速しないので抵抗減に振った。

とりあえず安定性設計についてはこんなところ。作ったシミュレータ、かなり遊べるので、これでパイロットをイジメつつ腕を上げてもらおう(笑)
機会があったらプロペラのことも記事にする、かな？

人力飛行機空力設計と粒子群最適化

ご無沙汰しておりました。Yuukivelです。
直近の前刷り原稿締め切りでバタバタして、お盆に入って実家でダラダラしております。
一人暮らしだとなかなか好物の梨とか食べられないし、やはり食の充実という意味で実家はいいですね！

さて、、

今回は以前記事にした人力飛行機設計について、OBチームにて実際に作ることになった機体の空力設計について解説する。
前記事はこちら→　http://d.hatena.ne.jp/yuukivel/20140101/1388592540

先に三面図を公開しておきましょう。ちょうど良い時期ですし、様々に参考にしていただければ。

一応PDFとXFLR5ファイルについては、まだ変更の可能性があるので今後公開するかも、としておく。

最適化について

本題の空力設計の話に入る前に、前提の知識として、「最適化」について解説しておく。
TLを見ていると、最適化を研究している者としてかなり気になる用法でこの言葉が使われていることがある。
あまり詳しく立ち入ると記事が長くなってしまうので、ここでは簡単に。

最適化とは、「ある評価関数があって、その評価関数を実行可能な範囲で最大・最小化すること」だ。
一応、最大化でも最小化でも、評価関数の正負を逆にすれば最大化にも最小化にもできるので、ここでは全て最小化、とする。

実行可能な範囲全てをくまなく調べて、本当にここより小さな解は存在しない、という解を大域的最適解という。
また、実行可能な範囲で、その解の周りだけみると最小となっている解を局所的最適解という。

ある解が大域的最適解だと言い切るためには、評価関数が以下である必要がある。
・実行可能な範囲で、関数の形が分かっている
・実行可能な範囲で、凸関数であることが証明できる

凸関数とは、例えば二次関数のように、谷が1つしか無いような関数のこと。

上のグラフはを表している。
このうち赤い○で囲ったのが「xが-3から3の範囲で」大域的最適解といえる。ちなみにこの関数であれば-3から0の範囲で凸関数。2階微分がこの範囲は正だからだ。
しかし、もしかしたらxが遙か彼方で谷を持つかもしれない。つまりは真の意味での大域的最適解は「神のみぞ知る」解といえる。

人間は厳密な意味での大域的最適解を得るのは不可能であり、ましてや関数の形が分からないし凸関数であることも証明できない問題には局所最適解で甘んじるしかないのだ。
局所最適解であることの必要十分条件としてKKT条件(Karush-Kuhn-Tucker condition)がある。これを満たせば、一応局所最適解ですよということは示される。
KKT条件について詳しくは立ち入らないが、実行可能かどうか、その点でのラグランジアンの傾きが0かどうか、といった条件をひとまとめにしたもの。

言いたかったことはなにかというと
・「最適化」とは評価関数があって、実行可能な範囲で局所最適解を求めること、ということ。それ「改善」じゃない？と言いたくなることが多々。
・最適性の条件を満たしていればそれはその評価関数において最適解であるということ。あとこの部分性能悪いから、それ最適解じゃないよとかいうのは、ちょっとこまる。だいたい評価関数の選び方はコンセプトで決まってしまうのだから。

最適解を求める方法には大きく分けて二種類あって
・勾配法(ニュートン法とか逐次二次計画法(SQP)とか)　
関数の傾きの情報を使って(局所)最適解を探索する。関数が微分可能でないと使えない。速いと言われるけどヤコビアンやヘッシアンが分かっていないとそれらを評価するのに時間がかかる。
・メタヒューリスティクス(遺伝的アルゴリズムとか粒子群最適化とか)
ランダムサーチの強化版みたいなもの。理論的根拠は希薄だが上手くいく。厳密な局所最適解は得られないが、関数が連続ですらなくても実行可能でロバスト。ちなみにXFOILを用いた最適化はFortranの小数点の問題(丸め誤差)のせいでメタヒューリスティクスでしか動かせません。

これらのうち長所短所をみながら問題に合わせて使っていく。

で、今回はメタヒューリスティクスの中でも粒子群最適化(Particle Swarm Optimizer:PSO)をつかって人力飛行機を設計しましたよというお話

空力設計(コンセプト)

この機体の空力設計におけるコンセプトは、以前の記事で作ったものと変わっていなくて、
・構造優先設計
・抵抗の徹底的な削減(薄翼化・高いテーパー比の尾翼等)
ということ。抵抗の削減であって出力の最小化ではないことも注意。取り回しの関係で翼幅は33.75mに固定。また設計速度は風に負けないように少々早めの7.5m/s。こうした制約の中で、そして存在する桁を使って、できるかぎり抵抗の小さい機体を設計しましょう、というコンセプト。

前者については以前の記事の通りで、桁を(上手い人に)設計してもらって、実際に作って重量を実測した上で、循環分布を設計した。
桁の剛性は重さの割に結構高く、たわみ2.5mを実現する循環分布を求めると楕円循環分布よりほんの少し揚力を真ん中によせた形になった。
そうこうして、各翼の端での循環の値が求まったとして、次の「α制約分散PSOによる主翼自動設計」の項に行きましょう。

空力設計(α制約分散PSOによる主翼自動設計)

各翼での循環値が定まったら、その循環値と指定した翼厚をみたす翼型の中から最小抵抗の翼型を設計する。
最適化問題風に記述すると以下のようになる。

翼型については4つの翼型の線形和と追加キャンバーで表現する。

1.使用する4翼型の翼型座標ファイルを用意しx座標の点の位置をそろえておく
2.係数に沿って翼型座標ファイルのy座標を単純に足し算する(引き算も可)
3.係数に従ってNACA4digit系列のキャンバーラインを作成する。
4.線形和で得た翼型のy座標に3.で得たキャンバーを足し算する。

以上の手順で翼型をパラメーターで表現する。以前GAでやっていたころは追加キャンバーをPARSECの応用でやっていたが、今回はNACA4digit系のキャンバーを利用した。
こうすることで得られる翼型がぐんと滑らかになった。

で、こうしてパラメータ表現された翼型に対しPSOを適用するのだが、そもそもPSOとは、、
Wikipediaの記事がシンプルでよいのでリンクを張っておきます。 粒子群最適化
PSOは魚群や鳥群の餌に群がる際の情報伝達をモデルとしたメタヒューリスティクスであり、同じメタヒューリスティクスである遺伝的アルゴリズム(GA)に比べて、実数問題や多峰性問題に向いていると言われる。
今回はPSOを並列実行できるように拡張し、octaveを7プロセス並列で一世代あたり490個体で計算させた。
PSOは各世代の位置と速度を、以下の式に従って更新する。

このときlbestは各プロセスでの過去から現在までの最良個体値、gbestは全プロセスを通した最良個体値とする。
ωやC,σは適当な定数。挙動をみながら調整する必要有り。

またgbestやlbestを決める上で最良個体を選ばなければならないが、今回は制約条件が存在するため、αレベル比較を実装した。
http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/1457-30.pdf
αレベル比較とは、制約条件の「満足度」を定義し、比較する両個体の満足度がともにα以上であれば目的関数の比較を、どちらかでもα以下であれば満足度の大小で比較を行う方法。
詳しい実装方法は論文に載っているので省略します。

このようにPSOにαレベル比較を実装することで、「欲しい循環を満たしながら、抵抗が最小となる翼型」を探索することが可能となった。
従って、
1.循環分布を各セクションに対して求める。
2.循環γ=0.5UCCLを満たす翼型を順次求める
3.各翼の中身を線形補間で翼型を求める。
のループを自動で回すことで、可能な限り抵抗の小さい主翼を求めることができる。
もちろんメタヒューリスティクスは厳密な最適解を求めることはできないから、実行毎に解が変わる可能性がある。
そこは自分が満足するのが出てきたらおわり、ということで(笑)
主翼最適設計とせず、主翼自動設計と言っているのは、こうしたところに理由がある。
1⇒2⇒3のループを繰り返し実行して、解を求めるようにすれば主翼最適設計といっても口がむずがゆくならない。

実際に設計を行った様子を動画にあげているので、もしよろしければ。

(見返したら翼自動最適設計といっていてむずがゆい、、確かに循環最適⇒翼型1最適⇒翼型2最適・・・となっているはずなので最適設計ではあるのだけれど）

おわりに

尾翼やペラの設計についても話したかったのだけれど、書く力が尽きてしまったので、いったんここで記事を締めくくろうと思う。
実際XFLR5の安定性解析結果をぱっとUnityで作ったフライトシミュに読み込んで安定性を議論・・・とかやっててそっちの話もしたいのだけれど、また今度ということで。。。
もうとにかく最適設計、最適設計、最適設計とやってきた機体なので実際にどう飛ぶか非常に楽しみ。

机上の空論となるか、計算は人間を追い越すか。

みたいな煽り文句をつければ、話題のタネになるかな(笑)

追記)三面図の厨二版あるんだけど、あまりに評判が悪いので悲しんでいます。上手く表現したつもりだったんだけどなぁ。。

PANAIRとPANROCKETとLaWGS

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

大学を移って早二十日。やっと航空を専門とすることができ、充実した生活を送っております。
今回の記事は、研究でもお世話になっているPDAS(Public Domain Aeronautical Software)について、と、このPDASの目玉プログラムである（であろう）「PANAIR」について紹介します。
PDASといっても名古屋の某航空宇宙ベンチャーじゃないですよ（笑）

それでは本題。

まえがきでも述べたとおり、PDAS(Public Domain Aeronautical Software)というコンテンツがある。
リンクは下記。
http://www.pdas.com/index.html
このページにははNASAを中心に1970年代から80年代ごろにかけて使用されていたプログラムが無償で公開されている。
中身はたとえば、このブログでもたびたび取り上げている「XFOIL」の前に用いられていた「EPPLER CODE」だったり、大気の状態を計算する「ATOMSPHERE」だったりかなり古かったりする。
（EPPLER CODEはXFOILにあるような境界層外縁速度を求めるイテレーションがなく、パネル法で求めた速度を境界層外縁速度としている。大体私が昔作って公開した翼型解析ソフトと似たようなもの）

ところが、古いといっても現在でも通用するプログラムがいくつかある。
たとえば「Hypersonic Arbitrary Body Program」
http://www.pdas.com/hyper.html
は、LaWGSフォーマットで記述された任意形状の物体の、超音速における空力係数をたとえば「Newtonian法」や「Tangent Wedge/Cone法」等を用いて求める。
実は私は研究でこのプログラムにとてもお世話になっている。現在においてもごく超音速域の空力係数を求めるのに上記の方法を用いるのは一般的だ。

そしてこのPDASの目玉プログラムのうちの一つであろうプログラムがタイトルにものせた
「PANAIR」
である。（http://www.pdas.com/panair.html）

PANAIRは名前が格好いいのもさることながら、その機能も素晴らしい。
その名前に含まれる通り、任意形状の機体周りの流れを一次元圧縮性修正の非粘性三次元パネル法で求める、というものだ。
一次元の圧縮性修正があるため（一応）局所マッハ数が4以下の流れを計算してくれる。ただしマニュアルにもある通り、遷音速では著しく精度が悪化する。

たとえば、私の研究対象であるOsculating Conical Waveriderのマッハ0.5における圧力係数分布を求めると下の画像のようになる。（もう古いので出してもいいでしょう）

航空機だけでなくこれまであまり計算されてこなかった（であろう）モデルロケットの圧力分布と抗力係数等を求めるのにも適している。
下の画像はA8モーターの搭載を想定したモデルロケットの迎角0度における圧力係数分布である。

PANAIRは基本ポテンシャル流であるので、一応高迎角の空力係数の推算も可能である。
ところがこれはプログラムが回るという意味であって、WAKE（後流）の設定がユーザー任せであることと粘性の影響がないことから精度についてはお察しである。

たとえば上記のロケットのCL・CD・CMcgを-180度から180度の範囲で計算すると以下のようになる。
実際には0度から90度まで計算してそのほかはひっくり返したり何なりして求めている。

各空力係数は動圧と胴体断面積で無次元化し、モーメント係数の基準長さは機体全長としている。
また重心は全長の50%位置にあるとして計算した。

各係数の補間には、RBF（Radial Based Function:放射基底関数）ネットワークを用いて補間した。スプラインや多項式補間よりよさそうな印象。
これについて立ち入ると記事ひとつできてしまうので略。

で、空力係数まで出してしまったら軌道計算するしかないということで、以前から作っているロケットシミュレータを組み込んだ「PANROCKET」というプログラムを作り、公開することとしよう。
本当はPANROCKETはpythonで書いてXGAGみたいなpyQtによるGUIを実装したいのだが、時間が取れるかどうか。。。
せっかくpyQtを使うのであればQopenGLを使って打ち上げのアニメーションを作りたいと思っている。

PANROCKETの使い方

コードはZipの中にまとめて入っています。

PANROCKET

PANROCKETはOctaveで動作し、「PANROCKET.m」が本体である。PANROCKET.mのなかで機体形状の定義と各種スクリプトの実行を行う。
軌道解析の設定は「PR_motion_init.m」で設定する。設定項目についてはmq_rocketと基本は同じだ。
PR_graphでdt*step秒ごとのロケット位置と姿勢の画像を出力し、アニメーションを作れる機能があるが、Ghostscriptのpathをoctaveに通していないと実行できなかったり、姿勢が描画されてなかったり、ご機嫌をとるのが難しい。

また、コードの文字エンコードは今回からshift-jisに戻した。以前までUTF-8を用いていたので、「文字化けして読めない！」といった声も多かった。

一応Zipの中には、panin.exeとpanair.exeを入れてあるが、cygwinでコンパイルしたので、環境が違うとたぶん動かない。
動かないという人は、前述のPDASのページからPANAIR.f90とPANIN.f90をダウンロードしてコンパイルして使ってほしい。
PANAIRというフォルダにpanair.exeとpanin.exeさえあれば、必要なものはすべてコードで用意して実行してくれる。

Macの人はがんばってPANAIR.mを改変してください。これがPANAIRの実行スクリプトです。
まぁ、octaveでコンパイルしてoctfile作ってもいいんだけど、それは最近存在を知ったので今回はパス。

そのうち機体重量と機体慣性モーメントを機体の形状から推算できるようにしようと思っているが、これもそのうち。

mq_rocketと比較して弾道頂点での軌道がかなりリアルになっている印象をうける。
また、αⅢロケットというモデルロケットを意識した感じで機体をインプットしたら80mぐらい上昇とシミュレートされ、ほぼ文献と一致した。
こうしたシミュレーションやっているものとしては、どれだけ現実と一致しているのか気になるところで、何とか精度検証したいものだ。
ただ本業の航空とちがって、モデルロケットは知識も経験もあまりないので難しいところ。

LaWGS(Langray Wireframe Geometry Standard)フォーマット

上の記事にもところどころ名前が出ているが、PDASのプログラムを使うに当たってLaWGSが扱えるようになると、格段に便利になる。
PANAIRの形状や計算条件のインプットには専用のインプットファイルが必要で、これの記述がまた難しい。
そこでPANAIRにはPANIN (http://www.pdas.com/panin.html) というプリプロセッサが用意されていて、WGSファイルと設定を読み込んでPANAIRのインプットファイルを作成してくれる。
（実はPANINあってもPANAIRをまわすのは大変で、パネルの向きや境界条件によって簡単にまわったりまわらなかったりする。パネルのネットワーク（郡）が複数あるともうカオスで、計算されたCpを見ながらパネルの向き等を調整するといった地獄の作業が始まる。ただし翼だけといったひとつのネットワークで記述できるものなら比較的簡単にまわすことができる）

前述のhyperもWGSさえつくればすぐに計算してくれる。

特筆して便利なのが、機体形状の可視化だ。
たとえばhlp(http://www.pdas.com/hlp.html)というプログラムは、PostScriptやGnuplotで任意の位置から見た機体形状を可視化してくれる。

またwgs2wrl(http://www.pdas.com/wgs2wrl.html)はwgsを直接vrml(たぶん1.0)に変換する。
たとえばcortona3d(http://www.cortona3d.com/cortona3d-viewers)とかを使えば、機体形状をワンクリックでぐりぐり見ることが可能だ。

私自身はWAKEもWGSで記述してしまうので、Viewerで見るときにはWAKEもくっついている（笑）

このようにPDASのソフトウェアを使うにあたってLaWGSは非常に強力な武器となる。
中身自体はいたって簡単で、たとえば上のモデルロケットのWGSの冒頭は以下のようになっている。

Created by Yuukivel from wgs_make
'NOSEup'
1 5 5 0   0 0 0   0 0 0    1 1 1  0 
50.0000000000 0.0000000000 12.5000000000 
37.5000000000 0.0000000000 12.1030729569 
25.0000000000 0.0000000000 10.8253175473 
12.5000000000 0.0000000000 8.2679728471 
0.0000000000 0.0000000000 0.0000000000 
50.0000000000 4.7835429046 11.5484941564 
37.5000000000 4.6316455013 11.1817813854 
25.0000000000 4.1426696754 10.0012893149 
12.5000000000 3.1640162278 7.6386108888 
0.0000000000 0.0000000000 0.0000000000 
50.0000000000 8.8388347648 8.8388347648 
37.5000000000 8.5581649610 8.5581649610 
25.0000000000 7.6546554462 7.6546554462 
12.5000000000 5.8463396668 5.8463396668 
0.0000000000 0.0000000000 0.0000000000 
50.0000000000 11.5484941564 4.7835429046 
37.5000000000 11.1817813854 4.6316455013 
25.0000000000 10.0012893149 4.1426696754 
12.5000000000 7.6386108888 3.1640162278 
0.0000000000 0.0000000000 0.0000000000 
50.0000000000 12.5000000000 0.0000000000 
37.5000000000 12.1030729569 0.0000000000 
25.0000000000 10.8253175473 0.0000000000 
12.5000000000 8.2679728471 0.0000000000 
0.0000000000 0.0000000000 0.0000000000 
'BODYup'
1 5 5 0   0 0 0   0 0 0    1 1 1  0 
287.0000000000 0.0000000000 12.5000000000 
227.7500000000 0.0000000000 12.5000000000 
168.5000000000 0.0000000000 12.5000000000 
109.2500000000 0.0000000000 12.5000000000 
50.0000000000 0.0000000000 12.5000000000 
287.0000000000 4.7835429046 11.5484941564 
227.7500000000 4.7835429046 11.5484941564 
168.5000000000 4.7835429046 11.5484941564 
109.2500000000 4.7835429046 11.5484941564 
50.0000000000 4.7835429046 11.5484941564 
287.0000000000 8.8388347648 8.8388347648 
227.7500000000 8.8388347648 8.8388347648 
168.5000000000 8.8388347648 8.8388347648 
109.2500000000 8.8388347648 8.8388347648 
50.0000000000 8.8388347648 8.8388347648 
287.0000000000 11.5484941564 4.7835429046 
227.7500000000 11.5484941564 4.7835429046 
168.5000000000 11.5484941564 4.7835429046 
109.2500000000 11.5484941564 4.7835429046 
50.0000000000 11.5484941564 4.7835429046 
287.0000000000 12.5000000000 0.0000000000 
227.7500000000 12.5000000000 0.0000000000 
168.5000000000 12.5000000000 0.0000000000 
109.2500000000 12.5000000000 0.0000000000 
50.0000000000 12.5000000000 0.0000000000 
(続く)

中身を見ていく。
1行目：タイトル

Created by Yuukivel from wgs_make

2行目:ネットワーク名。ネットワークは各パネル郡をまとめたものと思ってもらえれば。

'NOSEup'

3行目：ネットワークの情報
1つ目の１はネットワークナンバー（ほとんど関係ない）
2つ目の5は記述するライン（点の集まり）のネットワーク内における総数
3つ目の5は記述するラインの点の個数
この二つと4行目からの点の集まりがWGSの本質なので後述で詳しく説明します。
4つめはlocal座標である平面で対称かどうか（対称なしなら0、xz平面対称なら1、xy平面対称なら2、zy平面対称なら3）
その後の3つは局所座標の原点に対する回転。
その後の3つは局所座標の原点に対する移動。
その後の3つの1は各軸に対する拡大縮小。以上のこれら3つはまず使わない。
その後の0はglobal座標系で対称かどうか。

1 5 5 0   0 0 0   0 0 0    1 1 1  0 

そして4行め以降は (x y z)の座標点の繰り返しだ。

これらの記述の仕方はWGSのマニュアルに詳しく書いてある。(http://www.pdas.com/refs/tm85767.pdf)
マニュアルから抜粋した画像を載せる。これにあるとおり、各コンター（ライン）の並びで記述される。

このようにLaWGSを記述してやれば、PDASのソフトウェアを用いて様々なことを行う事ができる。
PANAIRは例えば人力飛行機のフェアリング設計や胴体と翼の干渉などの調査に活かすことができる。
hyperはSSTOやTSTOの設計に活かすことができる。

PDASのソフトはせっかく無償で公開されているのに、日本では殆ど存在が知られておらず、日本語ページの記述もほとんどない。
これらを有効活用すれば、個人の航空宇宙活動のレベルが高まるだろうと期待している。

Windmizeアップデート

時間が取れたので連投します。

直そう直そうと思っていたのですが時間が取れず、今やっとアップデートできました。
「Windmize」のアップデートです。

Windmize ver1.10
(64bit版 Windows7のみ動作確認しました。ご注意下さい)

ソースコードはこちら
https://github.com/NaotoMorita/Circulation_design

今回のアップデートは、UIは変えず、バグフィックスを行いました。

たわみ計算をするときに桁にかかる重力をkgfからNに変換してない。。。
気づいたときには青ざめたものでした。本当に申し訳なく思います。

堅くなる方にミスっていたのでまだ良かったですが、こういうミスはなくさなければなりませんね。

GPS「NMEA」フォーマットとロケットシミュ

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

卒論の検印もすみ、卒論発表も終わらせ、卒業まで一直線となりました。
２月末に論文2本の締め切りがあるわけですが、それさえ済ませればとうとう晴れて航空宇宙工学を生業とする学生となるわけです。

さて本題です。

最近訳あって、GPSロガーに手を出している。
飛行ロボットなりやっていくにあたってGPSの利用は必須なので、必要な知識というわけだ。

今回はGPSのフォーマット「NMEA」について。
そしてせっかくだから昔作ったロケットシミュ「mq_rocket」というソフトの結果をNMEAではけるようにした。

NMEAのフォーマットの説明はググればいくらでも出てくる。
例えば
http://www.hiramine.com/physicalcomputing/general/gps_nmeaformat.html
とか。

NMEAはさまざまなセンテンスの集まりであり、その情報には重複も含まれる。
疑似NMEAを作成するに当たって、今回は$GPGGAのセンテンスのみ利用した。３次元座標の表現にはこれで十分。
逆にmbedの解説記事とかで目にする$GPRMCだけにすると高度の情報が不足する。

疑似NMEAを作るにあたって注意すべきはチェックサム。これが間違っているとソフトはNMEAを読んでくれない。
チェックサムはの作り方について。
$以降のセンテンスを全て数値コードに変換する。
そののちそれらを前から順次排他的論理和xorをとり、16進数に変換する。

ここについては既存のコードを参考にした。

さて、NMEAを作るには軌道の緯度経度標高とジオイド高のデータが必要だ。
シミュレーションで出てくるのはXYZのデータであり、これを緯度経度に変換しなければならない。
そのために地球中心・地球固定直交座標系ECEF座標系を経由して求める必要がある。

これについての参考資料は以下のリンクより
理解するためのGPS測位計算プログラム入門(その1)

計算の手順としては、

射点の緯度経度標高ジオイド高
↓
射点のECEF座標
↓
シミュレーションの結果（地平座標）を回転行列と射点ECEF座標によりECEF座標に変換
↓
得られたECEF座標を緯度経度楕円体高に変換
↓
楕円体高を標高に変換。（標高は楕円体高からジオイド高を引けばよい）

という流れ。

その結果がこちら

これは
NMEA to KMZ file converter というソフトを使ってkmvに変換して表示させている。
http://homepage2.nifty.com/k8/gps/indexj.htm#003
3Dにチェックを入れるのを忘れずに。
本来であればgoogle earthでNMEAは直接読めるはずだが、なぜかできないのでkmzを経由した。

あと、久しぶりにこのプログラムを弄ったついでに、動画出力機能を付けた。
windowsメディアメーカーで作成できる。
出力の選択をすると、figフォルダ内に画像を大量に出力するので、これらをムービーメーカーで再生時間を合わせるだけ。
再生時間の設定をちゃんと同期するように気を付けて下さい。

その結果がこれ。上の軌道のやつです。

あと、いつも通り作ったプログラムを公開します。
mq_rocket

NMEAについてはよく分かったので、ロガーの製作をやりたい。
新作飛行ロボットの設計も終わったし、春休みはやりたいことがたくさんだ。