飛行力学入門-1

( Updated on  )
physics math aerodynamics control-theory

今回からは飛行力学について解説していきます。

最初は、航空機の運動方程式を導出します。（今回の内容）

次回は微小擾乱理論に基づき、4つの連立方程式を線形化していきます。

そうすることで、伝達関数や状態方程式を用いた運動解析が可能となります。

1 航空機の非線形運動方程式

1.0. 座標変換

静止座標系OI−XIYIZIO_I-X_IY_IZ_Iと動座標系O−XYZO-XYZを考える。 まず、静止座標系における速度、角速度を以下のように成分表示することにする。

v⃗c=[i⃗j⃗k⃗][U V W]\mathbf{\vec{v}_c}=\begin{bmatrix}\vec i&\vec j&\vec k\end{bmatrix} \begin{bmatrix}U\\\ V\\\ W \end{bmatrix} ω⃗=[i⃗j⃗k⃗][P Q R]\mathbf{\vec \omega}=\begin{bmatrix}\vec i&\vec j&\vec k\end{bmatrix} \begin{bmatrix}P\\\ Q\\\ R \end{bmatrix}

このとき、U=[i⃗j⃗k⃗]\mathbf U=\begin{bmatrix}\vec i&\vec j&\vec k\end{bmatrix}を基準基底といい、各座標軸方向の単位ベクトルを並べたものである。 また、v⃗c\mathbf{\vec v_c}、ω⃗\mathbf{\vec \omega}を幾何ベクトル、[UVW]T\begin{bmatrix}U&V&W\end{bmatrix}^Tや[UVW]T\begin{bmatrix}U&V&W\end{bmatrix}^Tを座標成分ベクトルと区別する。(矢印の有無で書き分ける。) このとき任意のベクトルr\mathbf rについて

ddtr⃗=r˙+ω⃗×r⃗\frac d {dt} {\vec r}=\dot {\mathbf r}+\vec \omega \times \vec r

が成り立つ。ただし、

r˙=dXdti⃗+dYdtj⃗+dZdtk⃗\dot {\mathbf r}=\frac {dX}{dt}\vec i+\frac {dY}{dt}\vec j+\frac {dZ}{dt}\vec k

は動座標系から見たv⃗c\vec v_cの時間変化率である。 更に、ω⃗\vec \omegaの外積をとることは

Ω=[0−R−Q −R0−P −Q−P0]\Omega=\begin{bmatrix}0& -R& -Q \\\ -R&0&-P \\\ -Q& -P&0\end{bmatrix}

の行列積を取ることと同義であり、このΩ\Omegaを外積行列という。

1.1. 動座標系で記述した運動方程式

Newtonの第二法則を動座標系で記述することを考える。

1.1.1. 機体重心の運動方程式

質点運動に関するNewtonの第二法則は運動量保存則

mdv⃗cdt=∑F⃗m\frac{d \vec v_c}{dt}=\sum \vec F

で記述される。 機体に固定された動座標系で表現すると、

m(Ωv⃗c+v⃗c˙)=U⃗∑Fm(\Omega\vec v_c +\dot {\vec v_c})=\vec{U}\sum \mathbf F

成分で書き下すと、以下の式を得る。

[U˙ V˙ W˙]=[0−R−Q −R0−P −Q−P0][U V W]+1m[FX FY FZ]\begin{bmatrix}\dot U \\\ \dot V \\\ \dot W\end{bmatrix} =\begin{bmatrix}0&-R&-Q \\\ -R&0&-P\\\ -Q&-P&0\end{bmatrix} \begin{bmatrix}U\\\ V\\\ W\end{bmatrix} +\frac 1 m \begin{bmatrix}F_X \\\ F_Y\\\ F_Z\end{bmatrix}

(Ω\Omegaは交代行列であることに注意)

1.1.2. 姿勢運動の方程式

回転運動に関するNewtonの第二法則は角運動量保存則

ddt∫Vr⃗×V⃗dm=M⃗\frac d {dt}\int_V\vec r\times \vec V dm=\vec M

で記述される。 機体に固定された動座標系で表現すると、

U⃗Jω˙+U⃗ΩJω=U⃗M\vec U \mathbf J \mathbf{ \dot \omega}+ \vec U \mathbf \Omega \mathbf J\mathbf \omega = \vec U \mathbf M

ただし、J\mathbf Jは慣性テンソルである。 U⃗−1\vec U^{-1}を両辺にかけて以下の式を得る。

ω˙=−J−1ΩJω+J−1M\mathbf{ \dot \omega}= -\mathbf J^{-1} \mathbf \Omega \mathbf J\mathbf \omega + \mathbf J^{-1} \mathbf M

M=[MNL]T\mathbf M=\begin{bmatrix}M&N&L\end{bmatrix}^Tとおけば、

[U˙ V˙ W˙]=[Ixx−Ixy−Ixz −IyxIyy−Iyz −Izx−IzyIzz]−1(−[0−R−Q −R0−P −Q−P0][Ixx−Ixy−Ixz −IyxIyy−Iyz −Izx−IzyIzz][P Q R]+[M N L])\begin{bmatrix}\dot U\\\ \dot V\\\ \dot W\end{bmatrix} =\begin{bmatrix}I_{xx}&-I_{xy}&-I_{xz}\\\ -I_{yx}&I_{yy}&-I_{yz}\\\ -I_{zx}&-I_{zy}&I_{zz}\end{bmatrix}^{-1} \Big( -\begin{bmatrix}0& -R & -Q \\\ -R & 0 & -P \\\ -Q & -P & 0 \end{bmatrix} \begin{bmatrix}I_{xx}&-I_{xy}&-I_{xz} \\\ -I_{yx}&I_{yy}&-I_{yz} \\\ -I_{zx} & -I_{zy} & I_{zz} \end{bmatrix} \begin{bmatrix}P \\\ Q \\\ R\end{bmatrix} +\begin{bmatrix}M \\\ N \\\ L\end{bmatrix} \Big)

1.2. 座標系の種類と座標成分の変換

1.2.1. 地面固定座標系

上で導いた二つの方程式を表現するために、一つ固定座標系を用意する必要がある。 そこでZZ軸を下向きにとる地上に固定された座標系XEYEZEX_EY_EZ_Eを導入する。 ここでは、XEYEZEX_EY_EZ_Eの回転は無視できるとし、慣性系の一つとして考えることができるものとする。

1.2.1. 機体軸系

座標原点を機体重心、XBX_B軸を機体の幾何学的基準線におき、ZBZ_B軸をXBX_Bと垂直な方向に機体の下方へとり、YBY_Bを右手直交系をなすように取った座標系を機体軸系という。 また、XBX_B、YBY_B、ZBZ_Bの3軸のことをそれぞれroll軸、pitch軸、yaw軸という。

1.2.2. 安定軸系

機体軸系をYBY_B軸まわりに迎え角−α-\alphaだけ回転させ、飛行速度を対称面に射影した方向にXSX_S軸を取る。このXSYSZSX_SY_SZ_S座標系を安定軸系という。 安定軸系では、常にZSZ_S成分の速度ベクトルは0となる。

1.2.3. 風軸系

対気速度ベクトルが対称面からβ\betaだけずれているとする。 このとき、安定軸系をZSZ_S軸まわりに横滑り角β\betaだけ回転させて、対気速度ベクトルの方向にXWX_W軸を取った座標系XWYWZWX_WY_WZ_Wを風軸系という。 α\alpha、β\betaに関する次の式を用いて、計算をできるだけ簡単にしていく。

tan⁡α=WU, tan⁡β=VU2+W2\tan\alpha=\frac W U,\ \tan\beta=\frac{V}{\sqrt{U^2+W^2}}

1.2.4. 慣性主軸系

慣性テンソルJ\mathbf Jは正対称行列だから、適当な座標軸を取ることにより対角行列にすることができる。 このときの座標系を慣性主軸系という。

1.3. 速度・角速度の座標変換

回転操作は線形変換であるため、ある行列R\mathbf Rの行列積をとることに置き換えることができる。 このR\mathbf Rを回転行列という。 ここでは、Euler角における姿勢表現から回転行列を求め、角速度成分の関係を表すキネマティック方程式を導く。

1.3.1. Euler角による姿勢表現

Euler角で姿勢を表現するときは 座標軸の回転順序が重要であり、航空機の場合は321系、つまりZ1Z_1軸→Y2Y_2→Z3Z_3の順に、それぞれΨ\Psi、Θ\Theta、Φ\Phiだけ回転させることで表現する。

1.3.2. 速度ベクトルの座標変換

地上固定座標系での速度をve\mathbf v_e、機体軸系での速度をvc=[UVW]T\mathbf v_c=\begin{bmatrix} U & V & W \end{bmatrix}^T すると、1.2.1.から地面固定座標系から機体軸系への速度の座標変換は

ve=RX3(Φ)(RY3(Θ)(RZ3(Ψ)vc)) =[cos⁡Ψsin⁡Ψ0 −sin⁡Ψcos⁡Ψ0 001][cos⁡Θ0−sin⁡Θ 010 sin⁡Θ0cos⁡Θ][100 0cos⁡Φsin⁡Φ 0−sin⁡Φcos⁡Φ][U V W] =[cos⁡(Ψ)cos⁡(Θ)sin⁡(Φ)sin⁡(Θ)cos⁡(Ψ)−sin⁡(Ψ)cos⁡(Φ)sin⁡(Φ)sin⁡(Ψ)+sin⁡(Θ)cos⁡(Φ)cos⁡(Ψ)sin⁡(Ψ)cos⁡(Θ)sin⁡(Φ)sin⁡(Ψ)sin⁡(Θ)+cos⁡(Φ)cos⁡(Ψ)−sin⁡(Φ)cos⁡(Ψ)+sin⁡(Ψ)sin⁡(Θ)cos⁡(Φ) −sin⁡(Θ)sin⁡(Φ)cos⁡(Θ)cos⁡(Φ)cos⁡(Θ)][U V W]\begin{aligned} \mathbf v_e&=R_{X_3}(\Phi)(R_{Y_3}(\Theta)(R_{Z_3}(\Psi)v_c))\\\ &=\begin{bmatrix}\cos\Psi&\sin\Psi&0 \\\ -\sin\Psi&\cos\Psi&0 \\\ 0&0&1\end{bmatrix} \begin{bmatrix}\cos\Theta&0&-\sin\Theta \\\ 0&1&0 \\\ \sin\Theta&0&\cos\Theta\end{bmatrix} \begin{bmatrix}1&0&0 \\\ 0&\cos\Phi&\sin\Phi \\\ 0&-\sin\Phi&\cos\Phi\end{bmatrix} \begin{bmatrix} U \\\ V \\\ W \end{bmatrix}\\\ &=\left[\begin{matrix}\cos{\left(\Psi \right)} \cos{\left(\Theta \right)} & \sin{\left(\Phi \right)} \sin{\left(\Theta \right)} \cos{\left(\Psi \right)} - \sin{\left(\Psi \right)} \cos{\left(\Phi \right)} & \sin{\left(\Phi \right)} \sin{\left(\Psi \right)} + \sin{\left(\Theta \right)} \cos{\left(\Phi \right)} \cos{\left(\Psi \right)}\\\sin{\left(\Psi \right)} \cos{\left(\Theta \right)} & \sin{\left(\Phi \right)} \sin{\left(\Psi \right)} \sin{\left(\Theta \right)} + \cos{\left(\Phi \right)} \cos{\left(\Psi \right)} & - \sin{\left(\Phi \right)} \cos{\left(\Psi \right)} + \sin{\left(\Psi \right)} \sin{\left(\Theta \right)} \cos{\left(\Phi \right)}\\\ - \sin{\left(\Theta \right)} & \sin{\left(\Phi \right)} \cos{\left(\Theta \right)} & \cos{\left(\Phi \right)} \cos{\left(\Theta \right)}\end{matrix}\right]\begin{bmatrix} U \\\ V \\\ W \end{bmatrix} \end{aligned}

と表現される。これを航法方程式という。

1.3.3. 角速度ベクトルの座標変換

Z1Z_1軸周りの角速度をもつ角速度ベクトルΨ˙\mathbf {\dot\Psi}を機体軸系の成分ベクトルに変換すると、1.3.1.よりRZ1(Ψ)R_{Z_1}(\Psi)→RY2(Θ)R_{Y_2}(\Theta)→RX3(Φ)R_{X_3}(\Phi)の順に作用させたから、

Ψ˙3=RX3(Φ)(RY3(Θ)(RZ3(Ψ)Ψ˙))\mathbf {\dot\Psi_3}=R_{X_3}(\Phi)(R_{Y_3}(\Theta)(R_{Z_3}(\Psi)\mathbf {\dot\Psi}))

Y2Y_2軸周りの角速度をもつ角速度ベクトルΘ˙\mathbf {\dot\Theta}を機体軸系の成分ベクトルに変換すると、1.3.1.よりRY2(Θ)R_{Y_2}(\Theta)→RX3(Φ)R_{X_3}(\Phi)の順に作用させたから、

Θ˙3=RX3(Φ)(RY3(Θ)Θ˙)\mathbf {\dot\Theta_3}=R_{X_3}(\Phi)(R_{Y_3}(\Theta)\mathbf {\dot\Theta})

X3X_3軸周りの角速度をもつ角速度ベクトルΦ˙\mathbf {\dot\Phi}を機体軸系の成分ベクトルに変換すると、1.3.1.よりRX3(Φ)R_{X_3}(\Phi)を作用させたから、

Φ˙3=RX3(Φ)Ψ˙\mathbf {\dot\Phi_3}=R_{X_3}(\Phi)\mathbf {\dot\Psi}

以上から

ω=RX3(Φ)(Ψ˙+RY3(Θ)(Θ˙+RZ3(Ψ)Ψ˙))\mathbf \omega =R_{X_3}(\Phi)(\mathbf {\dot\Psi}+R_{Y_3}(\Theta)(\mathbf {\dot\Theta}+R_{Z_3}(\Psi)\mathbf {\dot\Psi}))

が得られる。 しかし、実際にセンサーから得ることができるのはω\mathbf \omegaであるので、ω\mathbf\omegaから[Φ˙Θ˙Ψ˙]T\left[\begin{matrix}\dot\Phi&\dot\Theta&\dot\Psi\end{matrix}\right]^Tを得る式を計算する。 上の式を成分表示すると、

[P Q R]=[100 0cos⁡Φsin⁡Φ 0−sin⁡Φcos⁡Φ]([Φ˙ 0 0]+[cos⁡Θ0−sin⁡Θ 010 sin⁡Θ0cos⁡Θ]([0 Θ˙ 0]+[cos⁡Ψsin⁡Ψ0 −sin⁡Ψcos⁡Ψ0 001][0 0 Ψ˙])) =[10−sin⁡Θ 0cos⁡Φsin⁡Φcos⁡Θ 0−sin⁡Φcos⁡Φcos⁡Θ][Φ˙ Θ˙ Ψ˙]\begin{aligned} \begin{bmatrix}P\\\ Q\\\ R\end{bmatrix}&= \begin{bmatrix}1&0&0 \\\ 0&\cos\Phi&\sin\Phi \\\ 0&-\sin\Phi&\cos\Phi\end{bmatrix}(\begin{bmatrix}\dot\Phi\\\ 0\\\ 0\end{bmatrix}+\begin{bmatrix}\cos\Theta&0&-\sin\Theta \\\ 0&1&0 \\\ \sin\Theta&0&\cos\Theta\end{bmatrix}(\begin{bmatrix}0 \\\ \dot\Theta\\\ 0\end{bmatrix}+\begin{bmatrix}\cos\Psi&\sin\Psi&0 \\\ -\sin\Psi&\cos\Psi&0 \\\ 0&0&1\end{bmatrix}\begin{bmatrix}0\\\ 0\\\ \dot\Psi\end{bmatrix}))\\\ &=\begin{bmatrix}1&0&-\sin\Theta \\\ 0&\cos\Phi&\sin\Phi\cos\Theta \\\ 0&-\sin\Phi&\cos\Phi\cos\Theta\end{bmatrix} \begin{bmatrix}\dot\Phi\\\ \dot\Theta\\\ \dot\Psi\end{bmatrix} \end{aligned}

よって逆行列をかけてやれば

[Φ˙ Θ˙ Ψ˙]=[10−sin⁡Θ 0cos⁡Φsin⁡Φcos⁡Θ 0−sin⁡Φcos⁡Φcos⁡Θ]−1[PQR] =[1sin⁡Φtan⁡Θcos⁡Φtan⁡Θ 0cos⁡Φ−sin⁡Φ 0sin⁡Φ/cos⁡Θcos⁡Φ/cos⁡Θ][PQR]\begin{aligned} \begin{bmatrix}\dot\Phi\\\ \dot\Theta\\\ \dot\Psi\end{bmatrix}&= \begin{bmatrix}1&0&-\sin\Theta \\\ 0&\cos\Phi&\sin\Phi\cos\Theta \\\ 0&-\sin\Phi&\cos\Phi\cos\Theta\end{bmatrix}^{-1} \begin{bmatrix}P\\ Q\\ R\end{bmatrix}\\\ &=\begin{bmatrix}1&\sin\Phi\tan\Theta&\cos\Phi\tan\Theta \\\ 0&\cos\Phi&-\sin\Phi \\\ 0&\sin\Phi/\cos\Theta&\cos\Phi/\cos\Theta\end{bmatrix} \begin{bmatrix}P\\ Q\\ R\end{bmatrix} \end{aligned}

これはEulerのキネマティックス方程式と呼ばれる。

1.4. まとめ

剛体の非線形運動方程式は以下の４つの連立方程式として表される。

[U˙ V˙ W˙]=[0−R−Q −R0−P −Q−P0][U V W]+1m[FX FY FZ] [P˙ Q˙ R˙]=[Ixx−Ixy−Ixz −IyxIyy−Iyz −Izx−IzyIzz]−1(−[0−R−Q −R0−P −Q−P0][Ixx−Ixy−Ixz −IyxIyy−Iyz −Izx−IzyIzz][P Q R]+[M N L]) [Ue Ve We]=[cos⁡(Ψ)cos⁡(Θ)sin⁡(Φ)sin⁡(Θ)cos⁡(Ψ)−sin⁡(Ψ)cos⁡(Φ)sin⁡(Φ)sin⁡(Ψ)+sin⁡(Θ)cos⁡(Φ)cos⁡(Ψ)sin⁡(Ψ)cos⁡(Θ)sin⁡(Φ)sin⁡(Ψ)sin⁡(Θ)+cos⁡(Φ)cos⁡(Ψ)−sin⁡(Φ)cos⁡(Ψ)+sin⁡(Ψ)sin⁡(Θ)cos⁡(Φ) −sin⁡(Θ)sin⁡(Φ)cos⁡(Θ)cos⁡(Φ)cos⁡(Θ)][U V W] [Φ˙ Θ˙ Ψ˙]=[1sin⁡Φtan⁡Θcos⁡Φtan⁡Θ 0cos⁡Φ−sin⁡Φ 0sin⁡Φ/cos⁡Θcos⁡Φ/cos⁡Θ][P Q R]\begin{aligned} \begin{bmatrix}\dot U\\\ \dot V\\\ \dot W\end{bmatrix} &=\begin{bmatrix}0&-R&-Q\\\ -R&0&-P\\\ -Q&-P&0\end{bmatrix} \begin{bmatrix}U\\\ V\\\ W\end{bmatrix} +\frac 1 m \begin{bmatrix}F_X\\\ F_Y\\\ F_Z\end{bmatrix} \\\ \begin{bmatrix}\dot P\\\ \dot Q\\\ \dot R\end{bmatrix} &=\begin{bmatrix}I_{xx}&-I_{xy}&-I_{xz}\\\ -I_{yx}&I_{yy}&-I_{yz}\\\ -I_{zx}&-I_{zy}&I_{zz}\end{bmatrix}^{-1} \Big( -\begin{bmatrix}0&-R&-Q\\\ -R&0&-P\\\ -Q&-P&0\end{bmatrix} \begin{bmatrix}I_{xx}&-I_{xy}&-I_{xz}\\\ -I_{yx}&I_{yy}&-I_{yz}\\\ -I_{zx}&-I_{zy}&I_{zz}\end{bmatrix} \begin{bmatrix}P\\\ Q\\\ R\end{bmatrix} +\begin{bmatrix}M\\\ N\\\ L\end{bmatrix} \Big) \\\ \begin{bmatrix} U_e \\\ V_e \\\ W_e \end{bmatrix} &=\left[\begin{matrix}\cos{\left(\Psi \right)} \cos{\left(\Theta \right)} & \sin{\left(\Phi \right)} \sin{\left(\Theta \right)} \cos{\left(\Psi \right)} - \sin{\left(\Psi \right)} \cos{\left(\Phi \right)} & \sin{\left(\Phi \right)} \sin{\left(\Psi \right)} + \sin{\left(\Theta \right)} \cos{\left(\Phi \right)} \cos{\left(\Psi \right)}\\\sin{\left(\Psi \right)} \cos{\left(\Theta \right)} & \sin{\left(\Phi \right)} \sin{\left(\Psi \right)} \sin{\left(\Theta \right)} + \cos{\left(\Phi \right)} \cos{\left(\Psi \right)} & - \sin{\left(\Phi \right)} \cos{\left(\Psi \right)} + \sin{\left(\Psi \right)} \sin{\left(\Theta \right)} \cos{\left(\Phi \right)}\\\ -\sin{\left(\Theta \right)} & \sin{\left(\Phi \right)} \cos{\left(\Theta \right)} & \cos{\left(\Phi \right)} \cos{\left(\Theta \right)}\end{matrix}\right]\begin{bmatrix} U \\\ V \\\ W \end{bmatrix} \\\ \begin{bmatrix}\dot\Phi\\\ \dot\Theta\\\ \dot\Psi\end{bmatrix}&= \begin{bmatrix}1&\sin\Phi\tan\Theta&\cos\Phi\tan\Theta \\\ 0&\cos\Phi&-\sin\Phi \\\ 0&\sin\Phi/\cos\Theta&\cos\Phi/\cos\Theta\end{bmatrix} \begin{bmatrix}P\\\ Q\\\ R\end{bmatrix} \end{aligned}