Butterfly_Effect( )

地面と別れる方法

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

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

概要

Matlab互換のフリーソフトoctaveを使って、簡単な6自由度フライトシミュレーションを行う。6自由度というのは,xyzの3軸並進運動とピッチ/ロール/ヨーの3軸姿勢運動の自由度を合わせたもの。要するに「どんな飛行経路」を「どんな姿勢」で飛んで行くのかを3次元空間でシミュレートするフレームワークだ。1次近似として、定常飛行周りで線形化された系を計算する.

航空機の6自由度運動の定式化

テキストの1/2章には運動学や運動方程式など、取っつきにくい解説が並ぶ。ここは原理の話なので重要なのだけれど、実は後回しにしても良い。線形化された系(微小擾乱近似)で不足を感じたら、1/2章に立ち戻って非線形力学系を理解すればよいからだ。不足を感じるケースとは、高高度で高機動するために迎角60°で飛ぶことを考える場合等だ。そんな特殊な事情でもない限りは線形化の範疇なので、運動方程式を次の形に書くことができる.

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

 \mathbf{x}は飛行機の状態ベクトル \mathbf{u}は操作量ベクトル(舵やエンジンの動き)である.なお簡単のために制御項 \mathbf{u}を無視する。

 \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系統に分解して独立に解いてしまうのは不思議な感じがするけれど、飛行機はそういう性質をもった力学系なのである。

余談

独立な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つのシステムは同居しているだけであって,お互いに干渉しない.つまり別々に解いた時と結果は同じ。
前回の記事で示したスクリプトは後のことを考えて、2つのシステムを一つに束ねて解いている。
得られた \mathbf{x}の中身をplot()関数で描画してやると,こんな感じのデータが入っている。

おわりに

繰り返しになるが、微小擾乱の仮定のもとで線形化されていることは常に頭の片隅に置いておく必要がある。
 \dot{\mathbf{x}} = \mathbf{Ax + Bu}という形で運動方程式が書けて,さらにこれが定係数になっているのは、釣合状態からの微小擾乱の範疇での話。例えば,迎角 \alphaが大きくなり過ぎると失速現象が生じ、係数行列を変化させるが,そういう仕掛けはこのスクリプトには含まれていない。

次回はこの計算結果を利用して飛行軌道を可視化する。