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

今回と次回で，計算結果のポスト処理について解説する．

軌道データを生成する

これまでに，定常飛行状態からの微小擾乱の時間発展を計算してきた．
そこから飛行軌道を生成する方法は，テキスト5章の最初の方にサラッと書かれている．1章に基づいて軌道計算すればよいらしい．
1章では，慣性座標系と機体に埋め込まれた動座標系が存在することが説明されていた．
また，お互いの座標変換についても書かれている．オイラー角による3軸回転を用いる方法だ．
この辺り，実は力学(運動学)の話として奥が深くて，あまり詳細に立ち入ると紙面を食う．
余談だが3次元回転について一番真摯な態度で望んでいるのは衛星屋さんだと思われ，ちゃんと理解したい人には「人工衛星の力学と制御ハンドブック」がオススメ．ちなみに物理屋/ゲーム屋/航空屋/衛星屋のそれぞれの3次元回転の扱い方を比較すると，カルチャーの違いが出ていて面白い．

人工衛星の力学と制御ハンドブック―基礎理論から応用技術まで

さて，そろそろ軌道計算していこう．まず，この座標変換式が基本になる.

は慣性座標系での速度ベクトル，は動座標系上での速度ベクトル．3次元回転を扱っているから，回転行列Rは3×3．

三次元の回転行列を定義する．

関数を定義するため, "Rotation_Matrix.m"というファイル名で以下のスクリプトを保存する．
保存場所は，フライトシミュレーションのスクリプトと同じディレクトリが良い．

1; %各軸周りの回転行列の定義

%x軸回転
function R_x = Rotation_X(Psi)
	R_x = [ cos(Psi), sin(Psi), 0;
               -sin(Psi), cos(Psi), 0;
                       0,        0, 1 ];
end

%y軸回転
function R_y = Rotation_Y(Theta)
	R_y = [ cos(Theta), 0, -sin(Theta);
                         0, 1,           0;
                sin(Theta), 0,   cos(Theta) ];
end

%z軸回転
function R_z = Rotation_Z(Phi)
	R_z = [1,         0,        0;
               0,  cos(Phi), sin(Phi);
               0, -sin(Phi), cos(Phi) ];
end

慣性系から見た機体速度ベクトルの積分計算を行う

まず，の右辺のをはっきりさせておく．これは動座標系から見た飛行機の進行速度UVWだ．テキストでは元々を仮定しているので，としている．しかし特に近似する必要も無く，そのまま書けば

となる．これは迎角と横滑り角の定義から明らか．
これで，慣性系から見た機体速度ベクトルを計算するのに必要な情報はすべて揃った．
ここは横着して，これまでに書いたスクリプトを改築して飛行軌道も同時に解いてしまおう．

フライトシミュレーションをアップデートする

フライトシミュレーションのスクリプトを次のようにアップデートする．ファイル名は"aircraft_dynamics.m"とした．

%aircraft_dynamics.m
clear;cla
global A; global U0;

source("Rotation_Matrix.m"); %回転行列定義を読み込む

%有次元安定微係数
Xu = -0.01; Zu = -0.1; Mu = 0.001;
Xa = 30;    Za = -200; Ma = -4;
Xq = 0.3;   Zq = -5;   Mq = -1;
Yb = -45;   Lb_= -2;   Nb_= 1;
Yp = 0.5;   Lp_= -1;   Np_= -0.1;
Yr = 3;     Lr_= 0.2;  Nr_=-0.2;

%その他のパラメタ
W0 = 0;     U0 = 100;  theta0 = 0.05;
g  = 9.8; %重力加速度

%縦のシステム
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];

%横・方向のシステム
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];

%対角ブロックとしてシステムを結合する
A = blkdiag(A_lat,A_lon);

%さらに飛行軌道分のスペースを確保しておく
A = blkdiag(A,zeros(3));

function dx = dynamical_system(x,t)
% x = [u,alpha,q,theta,beta,p,r,phi,psi,x,y,z]
	global A; global U0;
	dx = A*x;

	u = x(1)+U0;%速度
	UVW = [u;u*x(5);u*x(2)];%速度ベクトル[U,V,W]
	dX = (Rotation_X(-x(9)) * Rotation_Y(-x(4)) * Rotation_Z(-x(8))) * UVW;
	dx(10) = dX(1);
	dx(11) = dX(2);
	dx(12) = dX(3);
end

%計算条件の設定
endurance	= 100;%飛行時間[sec]
step		= 10;%1.0[sec]あたりの時間ステップ数
t = linspace(0,endurance,endurance*step);

%初期値 x0 = [u,alpha,q,theta,	beta,p,r,phi,psi,	x,y,z]
x0_lat		= [10,0.1,0.4,0.2]';%縦の初期値
x0_lon		= [0.0,0.6,0.4,0.2,0.2]';%横・方向の初期値
x0_pos		= [0,0,-1000]';%飛行機の初期位置
x0 = vertcat(x0_lat,x0_lon,x0_pos);

%運動方程式を解く
x = lsode(@dynamical_system,x0,t);

printf('run successfully.\n')

実行してみる

実行して，新たに付け加えた状態量x(10),x(11),x(12)の中身を確認する．

octave:1>aircraft_dynamics
run successfully.
octave:2>plot3(x(:,10),x(:,11),x(:,12))

こんなのが表示される．

それらしい結果が出てきた．

次回（最終回）は，今回得られた軌道と姿勢の時系列データから飛行機のアニメーションを作る．

航空機力学入門

Fundamentals of Aerodynamics, 5E (Mcgraw Hill Series in Aeronautical and Aerospace Engineering)