PANAIRとPANROCKETとLaWGS - A Plane On The Sand

ご無沙汰しております。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座標系で対称かどうか。