Pythonでフライトシミュレーションをする第1回

最近やったことから書いてみようということで，Pythonを用いた航空機の飛行シミュレーションについて書きます．

1.モチベーション

つい数週間前に大学の夏学期が終わり，無為に時間を過ごしていた感が否めなかったので，復習がてら航空機力学，数値計算を扱える題材はないかなと思っていたところ，『航空機力学入門　(加藤寛一郎・大屋昭男・柄沢研治著　東京大学出版会)』の著者である加藤寛一郎さんの『壊れた尾翼　日航ジャンボ機墜落の真実 (講談社＋α文庫)』に出会いました．

航空機力学入門

作者:加藤寛一郎
メディア: 単行本

壊れた尾翼　日航ジャンボ機墜落の真実 (講談社＋α文庫)

作者:加藤寛一郎
発売日: 2013/09/27
メディア: Kindle版

1985年の有名な日航機墜落事故について解説する本(かなり本格的)なのですが，その中で出てくる事故機のブラックボックス(コックピットでの会話，飛行機の速度，姿勢などのデータを格納した，車で言うドラレコみたいなやつ)の飛行データによる飛行再現シミュレーションがかなり面白そうでした．

フライトシミュレーション，実装でググったところ，ちょうど「航空機力学入門の実装」みたいな記事

butterfly-effect.hatenablog.com

を見つけたのでやってやろう!と思いました．ただOctaveとかいう数値計算に特化した言語で書かれていて意味不明だったので，Python用に書き直す必要がありました．が，まあ結果的には連立微分方程式をただ脳死で打ち込むだけなので楽でした．

2.必要な方程式

航空機は3次元運動なので重心の運動方程式で3成分，剛体の運動方程式で3成分の計6本で記述されます．なお，以下の方程式は機体固定座標系での記述なので，飛行経路を計算する際に地面固定座標系に直すにはオイラー角の座標変換が必要です(それは追々)．

$m(\dot{U}+QW-RV)=-mg\sin{\varTheta}+X_a$

$m(\dot{V}+RU-PW)=mg\cos{\varTheta}\sin{\varPhi}+Y_a$

$m(\dot{W}+PV-QU)=mg\cos{\varTheta}\cos{\varPhi}+Z_a$

$\dfrac{d^{*}\mathbf{H}_c}{dt}＋\mathbf{\omega}\times\mathbf{H}_c=\mathbf{M}$

剛体運動方程式はちょっと打ち込むのが怠すぎたのでベクトル表記で省略．

$\mathbf{H}_c$ は機体重心周りの角運動量,

$\mathbf{\omega}=\begin{bmatrix}P&Q&R\end{bmatrix}^{T}$ は重心周りの角速度ベクトル,

$\mathbf{M}$ は重心周りの空気力，推力などによるモーメントです．詳しくはテキストを参照．

ご覧の通り，未知変数は機体飛行速度 $U,V,W$ ，そして角速度 $P,Q,R$ の6個の1階非線形微分方程式です．

このままの形だとさすがに扱いにくいので，線形化をします．

3.線形化された縦の運動方程式

航空機の運動が釣り合い状態からの微小擾乱(例えば $U \rightarrow U_{0} + u$ )だと仮定して近似します．(どの程度の擾乱まで許されるのかわからない．)

縦の運動

縦の状態量

$\mathbf{x}(t)=\begin{bmatrix} u&\alpha& q& \theta \\ \end{bmatrix}^{T}$

制御量

$\mathbf{u}(t)=\begin{bmatrix}\delta_e & \delta_t \end{bmatrix}^{T}$

とすると，以下のように書けます．

$\mathbf{A}=\begin{bmatrix}X_{u} & X_{\alpha} & 0 & -g \\ Z_{u}/U_{0} & Z_{\alpha}/U_{0} & 1+Z_{q}/U_{0} & 0 \\ M_{u}+M_{\dot{\alpha}}(Z_{u}/U_{0}) & M_{\alpha}+ M_{\dot{\alpha}}(Z_{\alpha}/U_{0}) & M_{q}+M_{\dot{\alpha}}(1+Z_{q}/U_{0}) & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix}$

$\mathbf{B}= \begin{bmatrix} 0 & X_{\delta_{t}} \\ Z_{\delta_{t}}/U_{0} & Z_{\delta_{t}}/U_{0} \\ M_{\delta_{e}}+M_{\dot{\alpha}} (Z_{\delta_{e}}/U_{0}) & M_{\delta_{t}}+M_{\dot{\alpha}} (Z_{\delta_{t}}/U_{0}) \\ 0 & 0 \end{bmatrix}$

に対し $\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u}$

行列の成分に出てくる $X_u$ とかは有次元安定微係数というやつで機体固有の定数です．注意が必要なのが安定軸で記述しているので

$\varTheta_{0}=W_{0}=0$

です．この方程式は線形なので実装しやすい形になりました．

これを4次のルンゲクッタで数値計算させました．(安定微係数は教科書に載っていたロッキード P2V-7機の値です．) Scipyのodeintモジュールを使えば良いのですが，復習のため封印しました．

import numpy as np
import matplotlib.pyplot as plt
from math import *

#以下，テキストの数値例を写経
#縦の有次元安定微係数
Xu=-0.0215
Zu=-0.227
Mu=0
Xa=14.7
Za=-236
Ma=-3.78
Madot=-0.28
Xq=0
Zq=-5.76
Mq=-0.992

X_deltat=0
Z_deltae=-12.9
Z_deltat=0
M_deltae=-2.48
M_deltat=0

#安定軸での水平釣り合い飛行条件での飛行速度などは以下の通り
U0=293.8
W0=0
g=9.80
theta0=0.0 

#講義プリントの固有運動のところを参考に，4次のルンゲクッタを実行．
#x=[u,a,theta,q] 運動変数
#縦の運動パラメータの微分計算
#u1=[delta_e,delta_t]　は操縦項
def deriv_pitch(x,u1):
    du=Xu*x[0]+Xa*x[1]-g*np.cos(theta0)*x[2]-W0*x[3]+X_deltat*u1[1]
    da=Zu/U0*x[0]+Za/U0*x[1]-g*np.sin(theta0)/U0*x[2]+(U0+Zq)/U0*x[3]+Z_deltae/U0*u1[0]+Z_deltat/U0*u1[1]
    dtheta=x[3]
    dq=(Madot*Zu/U0+Mu)*x[0]+(Ma+Madot*Za/U0)*x[1]-Madot*g*np.sin(theta0)/U0*x[2]+(Madot*(U0+Za)+Mq*U0)/U0*x[3]+(M_deltae+Madot*Z_deltae/U0)*u1[0]+(M_deltat+Madot*Z_deltat/U0)*u1[1]
    dx=np.array([du,da,dtheta,dq])
    return dx
#時間の刻み幅
dt=0.01
#時間の計算範囲
T= np.arange(0.0,300.0,dt)
N = len(T)
X = np.zeros((N,4))
#微小擾乱の初期値
u0=0
a0=0
q0=0
theta_ini=0
#操縦項を定数(操縦桿固定)
u1=[0.01,0]
#運動変数を格納する行列
X[0,:]=np.array([u0,a0,theta_ini,q0])

for i in range(1,N):
    k1 = deriv_pitch(X[i-1,:],u1)
    k2 = deriv_pitch(X[i-1,:]+0.5*dt*k1,u1)
    k3 = deriv_pitch(X[i-1,:]+0.5*dt*k2,u1)
    k4 = deriv_pitch(X[i-1,:]+dt*k3,u1)
    X[i,:] = X[i-1,:] + dt/6 * (k1+2*k2+2*k3+k4)
    

plt.plot(T,X[:,0])
plt.xlabel('time [s]')
plt.ylabel('u [ft/s]')
plt.grid()
plt.show()

plt.plot(T,X[:,1])
plt.xlabel('time [s]')
plt.ylabel('alpha(AoA) [rad]')
plt.grid()
plt.show()

plt.plot(T,X[:,2])
plt.xlabel('time [s]')
plt.ylabel('theta [rad]')
plt.grid()
plt.show()

plt.plot(T,X[:,3])
plt.xlabel('time [s]')
plt.ylabel('q [rad/s]')
plt.grid()
plt.show()