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