octaveで6自由度飛行シミュレーションを行う(1)
概要
Matlab互換のフリーソフトoctaveを使って、簡単な6自由度フライトシミュレーションを行う。6自由度というのは,xyzの3軸並進運動とピッチ/ロール/ヨーの3軸姿勢運動の自由度を合わせたもの。要するに「どんな飛行経路」を「どんな姿勢」で飛んで行くのかを3次元空間でシミュレートするフレームワークだ。1次近似として、定常飛行周りで線形化された系を計算する.
航空機の6自由度運動の定式化
テキストの1/2章には運動学や運動方程式など、取っつきにくい解説が並ぶ。ここは原理の話なので重要なのだけれど、実は後回しにしても良い。線形化された系(微小擾乱近似)で不足を感じたら、1/2章に立ち戻って非線形力学系を理解すればよいからだ。不足を感じるケースとは、高高度で高機動するために迎角60°で飛ぶことを考える場合等だ。そんな特殊な事情でもない限りは線形化の範疇なので、運動方程式を次の形に書くことができる.
は飛行機の状態ベクトル,は操作量ベクトル(舵やエンジンの動き)である.なお簡単のために制御項を無視する。
これは線形連立常微分方程式で,大学の応用数学を学んでいれば手計算でも解ける.
しかしもっと楽な方法の一つとして、ここではoctaveのODEソルバを使うことにする.
微分方程式を解くことよりも大事なこと、それは行列Aの中身である.
飛行機の運動は「縦」「横/方向」の2系統に分離されていることが分かる.
つまり,線形化の範疇ではピッチ運動とロール/ヨー運動は連成しない。そしてロール運動とヨー運動は連成している(ダッチロールとスパイラル)。
これは航空機特有の性質である。
縦の運動方程式
縦運動に関しては,テキストの(3.32)式で記述されている.
(3.32)式を行列の形で書き直せば次のようになる.
と近似すること,およびの関係を使って式変形することがポイントだ。
これは定係数なので、解くことは簡単。スクリプトは次のようになる。
縦運動のコード
clear % 縦のシステム global A_lat; A_lat =[Xu, Xa, -W0, -g*cos(theta0); Zu/U0, Za/U0, (U0+Zq)/U0, -g*sin(theta0)/U0; Mu, Ma, Mq, 0; 0, 0, 1, 0]; % 運動方程式の定義 % x = [u,α,q,θ]' function dx = lat_system(x,t) global A_lat; dx = A_lat * x; end % 10secを100回刻みにして微分方程式を解く t = linspace(0,10,100); x0 = [-5; 0.1; 0.8; 0.1];%テキトーな初期値 x_lat = lsode(@lat_system,x0,t);
横・方向の運動方程式
次に,横・方向運動についても同じようにする。
数式がテキストに書かれているので,行列表記に翻訳するだけで良い。
これも同じように解く。なお行列の係数は航空機力学入門の当該頁に記載されている。
(前回の記事で示したスクリプトを参考に).
横/方向運動のコード
clear % 横/方向のシステム global A_lon; A_lon = [Yb, (W0+Yp), -(U0-Yr), g*cos(theta0), 0; Lb_, Lp_, Lr_, 0, 0; Nb_, Np_, Nr_, 0, 0; 0, 1, tan(theta0), 0, 0; 0, 0, sec(theta0), 0, 0]; % 運動方程式の定義 % x = [β,p,r,φ,ψ]' function dx = lon_system(x,t) global A_lon; dx = A_lon * x; end % 10secを100回刻みにして微分方程式を解く t = linspace(0,10,100); x0 = [0.5;1.0;1.0;0.0;0.0]; %テキトーな初期値 x_lon = lsode(@lon_system,x0,t);
これでフライトシミュレータのコアとなる部分は出来上がり。得られたとがフライトログである。
一つの飛行機の動きを2系統に分解して独立に解いてしまうのは不思議な感じがするけれど、飛行機はそういう性質をもった力学系なのである。