octaveで6自由度飛行シミュレーションを行う(0)
はじめに
日本の飛行力学の教科書と言えば「航空機力学入門」。
一般向けの航空工学本と比べると内容が難しいので、学生の間では鬼門扱いされるとか。
反面、勘所(第3章)さえ押さえれば、一気に取っつき易くなる本でもある。
1~2章は前置き，3章が力学系の構築、4章以降は3章を多角的に考察している構成だからだ。
3章をoctaveで実装してみた
ともかく3章を実際に数値計算してみることが最短ルート。
テキストの中核部分である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の中には飛行状態履歴、いわゆるフライトログが格納されている。
