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

はじめに

日本語で書かれた飛行力学のテキストと言えば，「航空機力学入門」が有名．教科書の中ではやや難しい部類に入るので，一部の飛行機ファンの間では鬼門扱いされてるとかされてないとか．
ただ，コアの部分（3章）は恐ろしく単純なので，それについて具体的に示そうと思う．
1~2章はただの前置きだし，4章以降は3章で構築した力学系を様々な観点から考察してる内容なので，3章を理解することが1番大切.

octaveで実装してみた

ともかく実際に飛行力学を数値計算してみる．
テキストの中核部分である3章をoctaveで書いたらこうなる．
制御項をオミットしてやればこんなにも単純（制御項入れても単純だけど）．

clear;cla
global A;

%運動方程式
function dx = dynamical_system(x,t)
	% x = [u;α;q;θ; β;p;r;φ;ψ]
	global A;
	dx = A*x;
end

%有次元安定微係数
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);

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

%初期値 x0 = [u,alpha,q,theta, beta,p,r,phi,psi]
x0_lat		= [10;0.1;0.4;0.2]; %縦の初期値
x0_lon		= [0.0;0.6;0.4;0.2;0.2]; %横・方向の初期値
x0 = vertcat(x0_lat,x0_lon);

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

printf('run successfully.\n')

このスクリプトをoctave上で実行して，「run successfully.」と画面に表示されたら正常に動いてる．得られた変数ベクトルxの中には、いわゆるフライトログが格納されている．

おわりに

特にカッコイイ図が出てくるわけでもなく，計算も一瞬で終わってしまって，大変あっけなかったと思う．でも，このスクリプトで求めた計算結果をポスト処理して可視化してやれば，こんな感じの絵を書くことができる（可視化の方が大変だけど）．物理演算としてそれだけの内容をやっているということを押さえておくと、モチベーションにもなる.
f:id:aeroelasticity:20141117011527p:plain
載せたスクリプトの有次元安定微係数は適当に打ち込んだダミー値なので，ちゃんとした数字が欲しい場合はテキストを買って該当ページの数値を打ち込んでください（実験機の数値が載ってる）．例えばF22戦闘機の解析がしたいときはF22の安定微係数を何とか揃えて打ち込んでください．見つからないときは航空機力学入門に載ってる微係数推算式かDATCOMで一つ一つ求めるか，自力で空力解析&慣性モーメント計算しよう(?）．
今回はイントロということで，あえて解説をゴッソリ削って，スクリプトを示すに留めた．
やってることの解説や可視化方法については，次回以降ということで．