octaveで6自由度飛行シミュレーションを行う(1)

本編．前回の記事で示したスクリプトの中身の解説です．

概要

Matlab互換のフリーソフト octaveを使って，簡単な6自由度フライトシミュレーションを行う．6自由度というのは，xyzの三次元軌道とピッチロールヨーの三軸姿勢運動の自由度を合わせたもの．要するに「どんな飛行経路」を「どんな姿勢」で飛んで行くのかを三次元空間でシミュレートする．簡単に済ませるのが目的なので，有名なテキスト「航空機力学入門」に従って，定常飛行周りで線形化された系を計算する．

飛行機の6自由度運動の定式

テキストの1章2章には運動学の記述や力学の運動方程式など，ちょっと頭が痛くなるような話が載っている．ここは原理レベルの話なのでとても大切な話なのだけど，思い切ってスルー．実際に使うのは3章だけなので．微小擾乱近似（要するに線形化）の枠組みで満足できなくなったら，1章や2章を読んで非線形力学系を解けば良いと思う．
でも大抵はフライトダイナミクスを線形化できるということで，運動方程式を次の形に書こうというわけ．

$\dot{\mathbf{x}} = \mathbf{Ax + Bu}$

$\mathbf{x}$ は飛行機の状態ベクトル， $\mathbf{u}$ は操作量ベクトル（舵やエンジンの動き）である．とりあえず簡単のために模型飛行機を飛ばす要領で $\mathbf{u}$ を無視する．Let it beな飛行機になる．

$\dot{\mathbf{x}} = \mathbf{Ax}$

これは大学の応用数学でやる線形連立常微分方程式で，実は手計算でも解ける．
しかし僕は頭も手も使いたくないので、octaveのODEソルバを使うことにする．

問題は行列Aの中身である．
線形化することで，飛行機の運動は「縦」「横/方向」の2系統に分離される．
つまり，定常飛行からちょっと外れただけなら，ピッチング方向の運動とロール/ヨー方向の運動は連成しない．逆に言えば，線形化してもロールとヨーは連成する（連成現象の例：ダッチロールとスパイラル）．

縦の運動方程式

縦運動に関しては，テキストの(3.32)式で記述されている．
(3.32)式を行列の形で書き直せば次のようになる．

$\displaystyle \frac{d}{dt} \left( \begin{array}{c} u \\ \alpha \\ q \\ \theta \end{array} \right) = \left( \begin{array}{cccc} X_u & X_a & -W_0 & -g\cos{\theta_0} \\ Z_u/U_0 & Z_a/U_0 & (U_0+Z_q)/U_0 & -g\sin{\theta_0}/U_0\\ M_u & M_a & M_q & 0\\ 0 & 0 & 1 & 0 \end{array} \right) \left( \begin{array}{c} u \\ \alpha \\ q \\ \theta \end{array} \right)$

$M_{\dot{\alpha}}=0$ と近似すること，および $\dot{\theta} = q$ の関係を使って式変形することがポイント．
これは定係数なので，もう話は終わったようなものだ．コードは次のようになる．
行列の係数に出てくる数値はすべてテキストに載っているので，該当ページを探して各自打ち込んでください．

縦運動のコード

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);

横・方向の運動方程式

次に，横・方向運動についても同じようにする．
数式がテキストに書かれているので，粛々と行列表記に翻訳する．

$\displaystyle \frac{d}{dt} \left( \begin{array}{c} \beta \\ p \\ r \\ \phi \\ \psi \end{array} \right) = \left( \begin{array}{ccccc} Y_{\beta} & (W_0 + Y_p) & -(U_0 Y_r) & -g\cos{\theta_0}& 0 \\ L'_{\beta} & L'_p & L'_r & 0 & 0 \\ N'_{\beta} & N'_p & N'_r & 0 & 0\\ 0 & 1 & \tan{\theta_0} & 0 & 0 \\ 0 & 0 & \sec{\theta_0} & 0 & 0 \end{array} \right) \left( \begin{array}{c} \beta \\ p \\ r \\ \phi \\ \psi \end{array} \right)$

これも同じように解く．繰り返しになるけど，行列の係数の数値はテキストに全部書かれているから，それを打ち込んでください．
（前回の記事で示したスクリプトを参考に）．

横/方向運動のコード

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);

これでフライトシミュレータのコアは出来上がり．得られた $\mathbf{x}_{lon}$ と $\mathbf{x}_{lat}$ がフライトログというわけ．

一つの飛行機の動きを2系統に分解して独立に解いてしまうのは不思議な感じがするけれど，これで十分精度が出るというのだから，飛行機はそういう力学系なんだと受け入れることにする．

余談

独立な２つのシステムなので，1つにくっつけて書いても全く問題ない．
$\displaystyle \frac{d}{dt} \left( \begin{array}{c} \mathbf{x}_{lat}\\ \hdashline \mathbf{x}_{lon} \end{array} \right) = \left( \begin{array}{c:c} \mathbf{A}_{lat}&0\\ \hdashline 0&\mathbf{A}_{lon} \end{array} \right) \left( \begin{array}{c} \mathbf{x}_{lat}\\ \hdashline \mathbf{x}_{lon} \end{array} \right)$
行列を対角ブロックにして並べるにはblkdiag()関数を使えば良い．また，初期値ベクトルもvertcat()でくっ付ける．そして常微分方程式を解く(いつものlsode関数)．
2つのシステムは同居しているだけであって，お互いに全く干渉しない．つまり別々に解いた時と結果は同じ．
前回の記事で示したスクリプトは後のことを考えて，２つのシステムを一つに束ねて解いている．
得られた $\mathbf{x}$ の中身をplot()関数で描画してやると，こんな感じのデータが入っている（この図を書いたときにどんな初期値だったかは忘れたけど・・・）．
f:id:aeroelasticity:20141117023423p:plain

おわりに

ポイントなのは，微小擾乱の仮定が入っているということ．
$\dot{\mathbf{x}} = \mathbf{Ax + Bu}$ という形で運動方程式が書けて，さらにこれが定係数になっているのは釣り合い状態近傍で線形化しているからだ．しかし例えば，迎角 $\alpha$ が大きくなり過ぎると失速現象が起こることは有名．この現象は数理的には係数行列を変化させるはずなのだけど，そういう仕掛けはこのスクリプトには書いていない．その点は心に留めておく必要がある．失速を含むような激しいマニューバを計算すると精度が悪化する.