飛行器運動模型的建立

作者：由 Mond 發表于遊戲時間：2022-10-14

掌握了上一篇文章中的各項基礎知識，就可以順暢的完成飛行器運動方程的推導和模型的建立。基本思路就是利用上一篇文章中的動力學方程以及各座標系之間的變換關係來推導方程組。從而確定狀態量和狀態量微分之間的關係。有了該微分方程，就可以搭建出運動模型。

在本篇文章中，我們假設分析的物件飛行器是面對稱的，且飛行速度和高度有一定限制，足以我們把大地當作是平面且不旋轉的，那麼大地上的座標系就是慣性系。

座標系定義

值得一提的是，目前主流的兩種座標系定義方式，分別為美製座標系和蘇制座標系，同樣對應有兩種符號系統。在目前的論文和高校中，大家使用的多是美製座標系，也是現在的國標，主要是為了學術交流。但鑑於原國標——蘇制座標系的各種優點，在研發單位工程應用中，蘇制座標系和符號系統使用更為普遍。筆者做畢設時採用的美製座標系和符號系統，為了把兩種都熟悉一下，在本系列文章中我會使用蘇制座標系和符號系統。對於美製座標系和符號系統，讀者們可參閱吳森堂教授的《飛行控制系統》。

使用到的座標系定義如表

上面的座標系定義中，提到了空速和地速，這裡有必要做個區分：空速是相對未受飛行器擾動的氣流的速度，即相對空氣的速度；地速是飛行器相對地面的速度。當無風時，空氣相對地面靜止，此時空速與地速相等。但

當有風時，空速與地速是不相等的

。很多教科書和論文上對空速和地速不做區分，這是不合適的，因為這會使得模型無法加入風擾動。當有風時，假設風速為

$\boldsymbol{w}$

，則空速

$\boldsymbol{u}$

和地速

$\boldsymbol{v}$

的關係為

$\boldsymbol{v}=\boldsymbol{u}+\boldsymbol{w}\\$

另外，與地面座標系各軸平行的，原點位於飛行器質心的座標系，稱之為與飛行器牽連的地面座標系

，當不關係地面座標系位置時，可以用此座標系代替。

為了實現各座標系之間的轉換，我們引入下面各角度的定義。下面各角度的方向都是符合右手定則的，這也是蘇制座標系的一大好處。

機體座標系和地面座標系

飛行器的姿態角由機體座標系與地面座標系之間的轉換關係確定。

偏航角

$\psi$

：機體軸

在

平面內的投影與地面軸

的夾角，機頭向左偏為正；

俯仰角

$\vartheta$

：機體軸

與

平面的角度，機頭上抬為正；

滾轉角

$\gamma$

：機體軸

與包含機體軸

的鉛垂平面之間的角度，右翼下傾、左翼上傾為正。

將地面座標系按照

的順序分別轉動

$\psi$

、

$\vartheta$

和

$\gamma$

的角度，即可得到機體座標系，因此二者的轉換關係可確定為：

$\begin{equation} \begin{gathered}\boldsymbol{\mathrm{B}}_{d}^t=\boldsymbol{\mathrm{B}}_x(\gamma)\boldsymbol{\mathrm{B}}_z(\vartheta)\boldsymbol{\mathrm{B}}_y(\psi)\left[ \begin{array}{ccc} \cos\vartheta\cos\psi & \sin\vartheta & -\cos\vartheta\sin\psi\\ \sin\gamma\sin\psi-\cos\gamma\cos\psi\sin\vartheta & \cos\gamma\cos\vartheta & \cos\psi\sin\gamma+\cos\gamma\sin\psi\sin\vartheta \\ \cos\gamma\sin\psi+\sin\gamma\cos\psi\sin\vartheta & -\sin\gamma\cos\vartheta & \cos\gamma\cos\psi-\sin\gamma\sin\psi\sin\vartheta \end{array} \right] \end{gathered} \end{equation}\\$

氣流座標系和機體座標系

氣流角由氣流座標系和機體座標系之間的關係確定。

迎角

$\alpha$

：氣流系的

在飛行器對稱平面上的投影

與機體軸

之間的角度，

位於

下方為正；

側滑角

$\beta$

：氣流系的

與飛行器對稱平面之間的夾角，

在飛行器對稱平面右側為正。

由氣流系繞

軸轉過

$\beta$

角得到穩定（半機體）系：

$\begin{equation} \begin{gathered} \boldsymbol{\mathrm{B}}_q^b=\begin{bmatrix} \cos\beta&0&-\sin\beta\\ 0&1&0\\ \sin\beta&0&\cos\beta \end{bmatrix} \end{gathered} \end{equation}\\$

由穩定座標系繞

軸轉過

$\alpha$

角得到機體系：

$\begin{equation} \begin{gathered} \boldsymbol{\mathrm{B}}_b^t=\left[ \begin{array}{ccc} \cos\alpha & \sin\alpha & 0 \\ -\sin\alpha & \cos\alpha & 0 \\ 0 & 0 & 1 \end{array} \right] \end{gathered} \end{equation}\\$

由機體座標系到氣流座標系的旋轉矩陣為：

$\begin{equation} \begin{gathered} \boldsymbol{\mathrm{B}}_{q}^t=\left[ \begin{array}{ccc} \cos\alpha\cos\beta & \sin\alpha & -\cos\alpha\sin\beta \\ -\sin\alpha\cos\beta & \cos\alpha & \sin\alpha\sin\beta \\ \sin\beta & 0 & \cos\beta \end{array} \right] \end{gathered} \end{equation}\\$

寫到這裡，忽然想起，一般的教科書上這一部分都是會配一個座標軸的圖的，但平面的圖我覺得很難表現出立體的三軸旋轉，很是看不懂，所以我就不放圖了，大家只要知道這個旋轉的過程就好了。

航跡座標系和地面座標系

航跡角可由地面座標系和航跡座標系之間的關係定義。

航跡偏轉角（航向角）

$\psi_s$

：航跡系的

（也就是地速向量

$\boldsymbol{v}$

）在地平面上的投影與地面座標系的

的夾角，投影在

軸左側為正；

航跡傾斜角

$\theta$

：軸

（也就是地速向量

$\boldsymbol{v}$

）與地平面的夾角，向上飛為正。

地面座標系轉到航跡座標系可先繞

軸轉

$\psi_s$

，再繞

軸轉

$\theta$

。旋轉矩陣為

$\boldsymbol{\mathrm{B}}_d^h=\boldsymbol{\mathrm{B}}_z(\theta)\boldsymbol{\mathrm{B}}_y(\psi_s)=\begin{bmatrix} \cos\theta\cos\psi_s & \sin\theta & -\cos\theta\sin\psi_s\\ -\sin\theta\cos\psi_s & \cos\theta & \sin\theta\sin\psi_s\\ \sin\psi_s & 0 & \cos\psi_s \end{bmatrix}\\$

航跡座標系和氣流座標系

航跡座標系與氣流座標系的關係在有風時很複雜，但在無風時，空速與地速重合，也就是

與

重合，兩個座標系之間只有一個角度：繞速度向量的滾轉角，簡稱速度滾轉角。

速度滾轉角

$\gamma_s$

：氣流座標系

軸與鉛垂的

平面的夾角，

軸位於右面為正。

此時航跡座標系轉到氣流座標系的旋轉矩陣為

$\boldsymbol{\mathrm{B}}_h^q=\begin{bmatrix} 1 & 0 & 0\\ 0 & \cos\gamma_s & \sin\gamma_s\\ 0 & -\sin\gamma_s & \cos\gamma_s \end{bmatrix}\\$

運動方程的建立

作用在飛行器上的力和力矩

飛行器在空中所受力為發動機推力

$\boldsymbol{\mathrm{P}}$

，空氣動力

$\boldsymbol{\mathrm{R}}$

和重力

$\boldsymbol{\mathrm{G}}$

。其中，

$\boldsymbol{\mathrm{R}}$

又包含在

氣流座標系下

分解的三個方向的力

$\left[-Q,\;Y,\;Z\right]^T$

，分別為阻力、升力和側力；在

機體座標系下

分解的三個方向的力矩

$\left[M_x,\;M_y,\;M_z\right]^T$

，一般情況下發動機推力位於對稱平面內，設其與

夾角為

$\varphi_p$

，與質心的偏心距為

。（如果是二元推力向量的話，可能就需要加入另一對夾角和偏心距，按同樣的方法分析計算即可。）

飛行器所受

合力

在機體座標系下可表示為

$\begin{equation} \begin{gathered} \left[ \begin{array}{c} F_x \\ F_y \\ F_z \end{array} \right]=\left[ \begin{array}{c} P\cos\varphi_p \\ P\sin\varphi_p \\ 0 \end{array} \right]+\boldsymbol{\mathrm{B}}_q^t\left[ \begin{array}{c} -Q \\ Y \\ Z \end{array} \right]+\boldsymbol{\mathrm{B}}_d^t\left[ \begin{array}{c} 0 \\ -mg \\ 0 \end{array} \right] \end{gathered} \end{equation}\\$

飛行器所受

合力矩

在機體座標系下可表示為

$\begin{equation} \begin{gathered} \left[ \begin{array}{c} \sum M_x \\ \sum M_y \\ \sum M_z \end{array} \right]=\left[ \begin{array}{c} M_x \\ M_y \\ M_z \end{array} \right]+\left[ \begin{array}{c} 0 \\ 0 \\ -Pe_p \end{array} \right] \end{gathered} \end{equation}\\$

上面的力和力矩計算方式如下。

重力

可以透過飛行器的質量和目前所處位置的重力加速度得到。質量主要與飛行器的自重、掛載和燃油相關，可透過實驗測量得出。而重力加速度主要與所處經緯度和高度相關，現在已經可以透過模型可以計算得出。

推力

大小與飛行高度

，空速大小

，油門開度

$\delta_p$

相關。

$P = P(u,h,\delta_p)\\$

空氣動力和動力矩

分量可透過下式計算：

$\begin{equation} \begin{gathered} \begin{cases} Q=c_x \dfrac{1}{2}\rho u^{2} S \\ Y=c_y \dfrac{1}{2}\rho u^{2} S \\ Z=c_z \dfrac{1}{2}\rho u^{2} S \\ M_x= m_x \dfrac{1}{2}\rho u^2 Sb\\ M_y= m_y \dfrac{1}{2}\rho u^2 Sb\\ M_z= m_z \dfrac{1}{2}\rho u^2 Sc_A \end{cases} \end{gathered} \end{equation}\\$

式中

為飛行器

參考面積

，一般由結構設計相關負責人員給出，飛控設計人員不需要關心是怎麼算的；

$\rho$

為當前高度上的大氣密度；

為空速；

為

機翼展長

；

為

平均氣動弦長

。

機翼幾何引數定義

式中的

$1/2\rho u^2$

被稱作

動壓

，在飛機上又被稱作指示空速或錶速。飛機能否有足夠的升力安全飛行很大程度上就取決於錶速，因此是很重要的一個引數。

另外，上式中的

、

分別稱為

阻力系數

、

升力係數

和

側力系數

；

、

分別稱為

滾轉力矩係數

、

偏航力矩係數

、

俯仰力矩係數

。這些係數合稱氣動引數，由空氣動力學實驗得出。同樣，搞飛控的同學不需要關心是怎麼算的，但需要知道如何利用這些引數。

一般來講，上述的氣動引數和飛行器的狀態量相關。

$\begin{equation} \begin{gathered} \begin{cases} c_x=c_x(Ma,Re,\alpha,\delta_z;\beta,\delta_x,\delta_y)\\ c_y=c_y(Ma,\alpha,\delta_z; \dot{\alpha},\omega_z)\\ c_z=c_z(Ma,\beta,\delta_y;\omega_y, \dot{\beta})\\ m_x=m_x(Ma,\beta,\delta_x,\delta_y,\omega_x,\omega_y,\alpha)\\ m_y=m_y(Ma,\beta,\delta_y,\delta_x,\omega_y,\omega_x,\alpha; \dot{\beta})\\ m_z=m_z(Ma,\alpha,\delta_z,\omega_z, \dot{\alpha}, \dot{\delta_e}) \end{cases} \end{gathered} \end{equation}\\$

式中，

為馬赫數，

為雷諾數，

$\delta_x$

、

$\delta_y$

、

$\delta_z$

分別為副翼舵、方向舵和升降舵的偏轉角度。

有時給出的氣動引數是各個引數關於狀態量的導數，此時稱為氣動導數，就需要用氣動導數乘以該狀態量，然後做加和得到最終的力或力矩係數。如

$c_y=c_{y0}+c_{y\alpha}\alpha+c_{y\delta_z}\delta_z\\$

有時，給出的直接是在一系列狀態點下的不同的力和力矩係數，此時就需要透過插值或擬合方程的方式得到任一狀態點處對應的力和力矩係數。

事實上，狀態點不同，即使是氣動導數也會不同，還是需要透過插值運算或擬合方程來得到當前狀態點處的導數值。

另外需要提到的一點是，副翼舵左上偏右下偏為正；方向舵右偏為正；升降舵下偏為正。各舵偏正向定義均與右手定則四個手指的指向相吻合，這也是蘇制座標系的一個優點。而對於一般佈局（非鴨式佈局，方向舵和升降舵均置於重心之後）的飛機而言，

正向舵偏會產生負向力矩

。所以舵機模型均有一個

環節，必不可少。

動力學方程

對於飛行器按照

剛體的動力學

分析方法進行分析，為了方便查閱，把一般剛體的動力學方程寫在這裡。

剛體質心平移動力學方程

$m\left(\frac{\mathrm{d}}{\mathrm{d}t}\begin{bmatrix} v_x\\v_y\\v_z \end{bmatrix}+\begin{bmatrix} 0&-\omega_z&\omega_y\\ \omega_z&0&-\omega_x\\ -\omega_y&\omega_x&0 \end{bmatrix}\begin{bmatrix} v_x\\v_y\\v_z \end{bmatrix} \right)=\begin{bmatrix} F_x\\F_y\\F_z \end{bmatrix}\\$

式中各分量均為在某一旋轉座標系中的分解，該座標系角速度在其自身中的分解為

$[\omega_x,\;\omega_y,\;\omega_z]^\mathrm{T}$

。

剛體繞質心轉動動力學方程

$\frac{\mathrm{d}}{\mathrm{d}t}\begin{bmatrix} H_x\\H_y\\H_z \end{bmatrix}+\begin{bmatrix} 0&-\omega_z&\omega_y\\ \omega_z&0&-\omega_x\\ -\omega_y&\omega_x&0 \end{bmatrix}\begin{bmatrix} H_x\\H_y\\H_z \end{bmatrix} =\begin{bmatrix} M_x\\M_y\\M_z \end{bmatrix}\\$

該式各分量是在與剛體固連的座標系中的分解，也就是說剛體的角速度和該座標系的角速度相同，在該座標系中的分解為

$[\omega_x,\;\omega_y,\;\omega_z]^\mathrm{T}$

。而動量矩計算方式如下。

$\begin{bmatrix} H_x\\H_y\\H_z \end{bmatrix}=\begin{bmatrix} I_{x}&-I_{xy}&-I_{zx}\\ -I_{xy}&I_y&-I_{yz}\\ -I_{zx}&-I_{yz}&I_z \end{bmatrix}\begin{bmatrix} \omega_x\\\omega_y\\\omega_z \end{bmatrix}\\$

由此我們可以知道，飛行器動力學方程能夠分成兩部分：

質心平移的動力學方程，該方程可以分別在

機體座標系

、

航跡座標系

和

氣流座標系

等中進行描述；

繞質心轉動的動力學方程，該方程需要在與飛行器固連的座標系中進行描述，也就是

機體座標系

。

下面將對其進行推導，請讀者和我一起，一定要自己推一遍。可以先自己推導一下，再和我的結果對照一下來驗證，這樣加深理解。

Mathematica

和

MATLAB

是兩個不錯的計算軟體，在推導時可以靈活應用。

質心平移的動力學方程

機體座標系下

地速在機體座標系下的分解寫作

$[v_{xt},v_{yt},v_{zt}]^\mathrm{T}$

，而機體座標系的旋轉角速度就是飛機本身的旋轉角速度，在機體座標系下的分解寫作

$[\omega_x,\omega_y,\omega_z]^\mathrm{T}$

。飛行器所受合力在機體座標系中的表示上面已給出。那麼代入一般剛體的質心平移動力學方程就為下式。

$m\left(\frac{\mathrm{d}}{\mathrm{d}t}\begin{bmatrix} v_{xt}\\v_{yt}\\v_{zt} \end{bmatrix}+\begin{bmatrix} 0&-\omega_z&\omega_y\\ \omega_z&0&-\omega_x\\ -\omega_y&\omega_x&0 \end{bmatrix}\begin{bmatrix} v_{xt}\\v_{yt}\\v_{zt} \end{bmatrix} \right)=\left[ \begin{array}{c} P\cos\varphi_p \\ P\sin\varphi_p \\ 0 \end{array} \right]+\boldsymbol{\mathrm{B}}_q^t\left[ \begin{array}{c} -Q \\ Y \\ Z \end{array} \right]+\boldsymbol{\mathrm{B}}_d^t\left[ \begin{array}{c} 0 \\ -mg \\ 0 \end{array} \right]$

一般將微分項寫在等式的左邊，其餘項放在右邊，做移項後就是下式。

$\frac{\mathrm{d}}{\mathrm{d}t}\begin{bmatrix} v_{xt}\\v_{yt}\\v_{zt} \end{bmatrix}=-\begin{bmatrix} 0&-\omega_z&\omega_y\\ \omega_z&0&-\omega_x\\ -\omega_y&\omega_x&0 \end{bmatrix}\begin{bmatrix} v_{xt}\\v_{yt}\\v_{zt} \end{bmatrix}+ \frac{1}{m}\left(\left[ \begin{array}{c} P\cos\varphi_p \\ P\sin\varphi_p \\ 0 \end{array} \right]+\boldsymbol{\mathrm{B}}_q^t\left[ \begin{array}{c} -Q \\ Y \\ Z \end{array} \right]+\boldsymbol{\mathrm{B}}_d^t\left[ \begin{array}{c} 0 \\ -mg \\ 0 \end{array} \right]\right)\\$

展開後得到的結果如下。

$\begin{cases} \dot{v}_{xt}=\omega_zv_{yt}-\omega_yv_{zt}+\left( P\cos\varphi_p -Q\cos\alpha\cos\beta+Y\sin\alpha-Z\cos\alpha\sin\beta -mg\sin\vartheta \right)/m\\ \dot{v}_{yt}=\omega_xv_{zt}-\omega_zv_{xt}+\left( P\sin\varphi_p +Q\sin\alpha\cos\beta+Y\cos\alpha+Z\sin\alpha\sin\beta -mg\cos\gamma\cos\vartheta \right)/m\\ \dot{v}_{zt}=\omega_yv_{xt}-\omega_xv_{yt}+\left( -Q\sin\beta+Z\cos\beta +mg\sin\gamma\cos\vartheta \right)/m\\ \end{cases}\\$

航跡座標系下

在航跡座標系下的地速可寫作

$[v,0,0]^\mathrm{T}$

，根據從地面繫到航跡系的轉換過程，可以知道航跡座標系相對慣性系（地面系）的轉動角速度可分解為沿

軸的

$\dot\psi_s$

和沿

軸的

$\dot\theta$

。因此航跡座標系的轉動角速度在航跡座標系下的表示為

$\begin{bmatrix} \omega_{h,xh} \\ \omega_{h,yh} \\ \omega_{h,zh} \end{bmatrix}=\boldsymbol{\mathrm{B}}_d^h \begin{bmatrix} 0 \\ \dot\psi_s \\ 0 \end{bmatrix} + \begin{bmatrix} 0 \\ 0 \\ \dot\theta \end{bmatrix} = \begin{bmatrix} \dot\psi_s\sin\theta \\ \dot\psi_s\cos\theta \\ \dot\theta \end{bmatrix}\\$

則動力學方程可表示為

$m\left(\frac{\mathrm{d}}{\mathrm{d}t}\begin{bmatrix} v\\0\\0 \end{bmatrix} + \begin{bmatrix} 0&-\dot\theta&\dot\psi_s\cos\theta\\ \dot\theta&0&-\dot\psi_s\sin\theta\\ -\dot\psi_s\cos\theta&\dot\psi_s\sin\theta&0 \end{bmatrix}\begin{bmatrix} v\\0\\0 \end{bmatrix}\right) = \boldsymbol{\mathrm{B}}_q^h\boldsymbol{\mathrm{B}}_t^q\left[ \begin{array}{c} P\cos\varphi_p \\ P\sin\varphi_p \\ 0 \end{array} \right]+\boldsymbol{\mathrm{B}}_q^h\left[ \begin{array}{c} -Q \\ Y \\ Z \end{array} \right]+\boldsymbol{\mathrm{B}}_d^h\left[ \begin{array}{c} 0 \\ -mg \\ 0 \end{array} \right]\\$

在有風的情況下，從氣流繫到航跡系的變換矩陣將十分複雜，我們這裡不作計算，也不用它來建模。但在

無風

時，變換矩陣就會變得很簡略。此時將上式展開如下

$\begin{cases} m\dot v = P\cos(\alpha+\varphi_p)\cos\beta -Q -mg\sin\theta\\ mv\dot\theta = P[\cos(\alpha+\varphi_p)\sin\beta\sin\gamma_s+\sin(\alpha+\varphi_p)\cos\gamma_s] +Y\cos\gamma_s-Z\sin\gamma_s -mg\cos\theta\\ -mv\cos\theta\dot\psi_s = P[-\cos(\alpha+\varphi_p)\sin\beta\cos\gamma_s+\sin(\alpha+\varphi_p)\sin\gamma_s] +Y\sin\gamma_s+Z\cos\gamma_s \end{cases}\\$

氣流座標系下

首先計算重力在氣流座標系下的表示。

$\begin{bmatrix} G_{xq} \\ G_{yq} \\ G_{zq} \end{bmatrix} = \boldsymbol{\mathrm{B}}_t^q\boldsymbol{\mathrm{B}}_d^t\left[ \begin{array}{c} 0 \\ -mg \\ 0 \end{array} \right]\\$

展開可得

$\begin{cases} G_{xq}=mg(\sin\alpha\cos\beta\cos\gamma\cos\vartheta-\cos\alpha\cos\beta\sin\vartheta+\sin\beta\sin\gamma\cos\vartheta)\\ G_{yq}=-mg(\cos\alpha\cos\gamma\cos\vartheta+\sin\alpha\sin\vartheta)\\ G_{zq}=mg(-\sin\alpha\sin\beta\cos\gamma\cos\vartheta+\cos\alpha\sin\beta\sin\vartheta+\cos\beta\sin\gamma\cos\vartheta) \end{cases}\\$

設氣流座標系相對機體座標系的角速度在氣流座標系下的分解為

$[\omega_{q\rightarrow t,xq},\,\omega_{q\rightarrow t,yq},\,\omega_{q\rightarrow t,zq}]$

。根據座標變換來分析，氣流系的轉動角速度由繞

軸的

$-\dot{\alpha}$

和繞

軸的

$-\dot{\beta}$

。

$\begin{bmatrix} \omega_{q\rightarrow t,xq} \\ \omega_{q\rightarrow t,yq} \\ \omega_{q\rightarrow t,zq} \end{bmatrix} = \boldsymbol{\mathrm{B}}_t^q\begin{bmatrix} 0 \\ 0 \\ -\dot\alpha \end{bmatrix}+\begin{bmatrix} 0 \\ -\dot\beta \\ 0 \end{bmatrix} = \begin{bmatrix} -\dot\alpha\sin\beta \\ -\dot\beta \\ -\dot\alpha\cos\beta \end{bmatrix}\\$

則氣流座標系相對慣性系（地面系）的角速度就是相對機體系的角速度和機體系相對地面系的角速度的疊加。在氣流座標系下表示為

$\begin{bmatrix} \omega_{qxq} \\ \omega_{qyq} \\ \omega_{qzq} \end{bmatrix} = \begin{bmatrix} -\dot\alpha\sin\beta \\ -\dot\beta \\ -\dot\alpha\cos\beta \end{bmatrix} + \boldsymbol{\mathrm{B}}_t^q\begin{bmatrix} \omega_x \\ \omega_y \\ \omega_z \end{bmatrix} = \begin{bmatrix} -\dot\alpha\sin\beta+\omega_x\cos\alpha\cos\beta+\omega_y\sin\alpha\cos\beta+\omega_z\sin\beta\\ -\dot\beta+\omega_x\sin\alpha+\omega_y\cos\alpha\\ -\dot\alpha\cos\beta-\omega_x\cos\alpha\sin\beta+\omega_y\sin\alpha\sin\beta+\omega_z\cos\beta \end{bmatrix}\\$

在有風的情況下，地速在氣流座標系下的分解會稍複雜。且在飛行過程中，更加關心空速，控制物件一般也為飛行器的空速而不是地速。但動力學方程中的速度是地速，這樣就涉及到空速和地速之間的關係，在氣流座標系下會十分複雜。因此我們只考慮無風條件。在

無風

時，地速可用空速代替，氣流座標系下寫作

$[u,\;0,\;0]^\mathrm{T}$

。動力學方程為

$m\left(\frac{\mathrm{d}}{\mathrm{d}t}\begin{bmatrix} u\\0\\0 \end{bmatrix} + \begin{bmatrix} 0&-\omega_{qzq}&\omega_{qyq}\\ \omega_{qzq}&0&-\omega_{qxq}\\ -\omega_{qyq}&\omega_{qxq}&0 \end{bmatrix}\begin{bmatrix} u\\0\\0 \end{bmatrix}\right) = \boldsymbol{\mathrm{B}}_t^q\left[ \begin{array}{c} P\cos\varphi_p \\ P\sin\varphi_p \\ 0 \end{array} \right]+\left[ \begin{array}{c} -Q \\ Y \\ Z \end{array} \right]+\begin{bmatrix} G_{xq} \\ G_{yq} \\ G_{zq} \end{bmatrix}\\$

將上式展開得

$\begin{cases} m\dot u=P\cos(\alpha+\varphi_p)\cos\beta - Q+G_{xq}\\ mu\dot\beta=P\cos(\alpha+\varphi_p)\sin\beta - Z+mu(\omega_x\sin\alpha+\omega_y\cos\alpha)-G_{zq}\\ mu\cos\beta\dot\alpha=-P\sin(\alpha+\varphi_p)-Y+mu(-\omega_x\cos\alpha\sin\beta+\omega_y\sin\alpha\sin\beta+\omega_z\cos\beta)-G_{yq} \end{cases}\\$

繞質心轉動的動力學方程

轉動動力學方程在機體座標系下進行描述，飛行器的轉動角速度在機體座標系下的分解表示為

$[\omega_x,\;\omega_y,\;\omega_z]^\mathrm{T}$

。代入一般剛體的轉動動力學方程中可得

對於面對稱飛行器來說，

$I_{yz}=I_{zx}=0$

，因此動量矩可寫作

$\begin{bmatrix} H_x\\H_y\\H_z \end{bmatrix}=\begin{bmatrix} I_{x}&-I_{xy}&0\\ -I_{xy}&I_y&0\\ 0&0&I_z \end{bmatrix}\begin{bmatrix} \omega_x\\\omega_y\\\omega_z \end{bmatrix}=\begin{bmatrix} \omega_xI_x-\omega_yI_{xy}\\ -\omega_xI_{xy}+\omega_yI_y\\ \omega_zI_z \end{bmatrix}\\$

將動力學方程展開可得

$\begin{equation} \begin{gathered} I_x\dot\omega_x-I_{xy}\dot\omega_y+\omega_x\omega_zI_{xy}-\omega_y\omega_zI_y+\omega_y\omega_zI_z=M_x\\ -I_{xy}\dot\omega_x+I_y\dot\omega_y+\omega_x\omega_zI_x-\omega_y\omega_zI_{xy}-\omega_x\omega_zI_z=M_y\\ I_z\dot\omega_z-\omega_x\omega_yI_x+{\omega_y}^2I_{xy}-{\omega_x}^2I_{xy}+\omega_x\omega_yI_y=M_z-Pe_p \end{gathered} \end{equation}\\$

寫作容易解出的形式如下

$\begin{cases} \dot\omega_x=\dfrac{1}{{I_xI_y-I_{xy}}^2}\left[ I_yM_x+I_{xy}M_y-I_{xy}(I_x+I_y-I_z)\omega_x\omega_z+({I_y}^2+{I_{xy}}^2-I_yI_z)\omega_y\omega_z \right]\\ \dot\omega_y=\dfrac{1}{I_xI_y-{I_{xy}}^2}\left[ I_xM_y+I_{xy}M_x-I_{xy}(-I_x-I_y+I_z)\omega_y\omega_z+(-{I_x}^2-{I_{xy}}^2+I_xI_z)\omega_x\omega_z \right]\\ \dot\omega_z=\dfrac{1}{I_z}\left[ M_z-Pe_p-(-I_x+I_y)\omega_x\omega_y+I_{xy}({\omega_x}^2-{\omega_y}^2) \right] \end{cases}\\$

運動學方程

運動學方程主要由速度和角速度從不同座標系變換到地面座標系推導而出。在地面系下，飛行器質心所處位置表示為

$\begin{bmatrix} x_d \\ y_d \\ z_d \end{bmatrix}\\$

其中，

代表了飛行器所處的高度，又可寫作

。

那麼飛行器質心的運動速度，也就是地速可寫作

$[x_d,\;y_d,\;z_d]^\mathrm{T}$

。

機體座標系下

把速度分量由機體座標系變換到地面座標系，使用下式：

$\begin{cases} \dot{x}_d=v_{xt}\cos\vartheta\cos\psi+v_{yt}(\sin\gamma\sin\psi-\cos\gamma\sin\vartheta\cos\psi)+v_{zt}(\cos\gamma\sin\psi+\sin\gamma\sin\vartheta\cos\psi) \\ \dot{y}_d=v_{xt}\sin\vartheta+v_{yt}\cos\gamma\cos\vartheta-v_{zt}\sin\gamma\cos\vartheta \\ \dot{z}_d=-v_{xt}\cos\vartheta\sin\psi+v_{yt}(\sin\gamma\cos\psi+\cos\gamma\sin\vartheta\sin\psi)+v_{zt}(\cos\gamma\cos\psi-\sin\gamma\sin\vartheta\sin\psi) \end{cases}\\$

飛機繞質心轉動的角速度向量在機體座標系下分解為

$[\omega_x,\;\omega_y,\;\omega_z]^\mathrm{T}$

，根據地面座標系轉到機體座標系的過程，其又可分解為沿

的

$\dot{\psi}$

、沿

的

$\dot{\vartheta}$

和沿

的

$\dot{\gamma}$

。而沿

的

$\dot{\vartheta}$

又可分解為沿

的

$\dot{\vartheta}\cos\psi$

和沿

的

$\dot{\vartheta}\sin\psi$

。則角速度間的關係如下

$\begin{bmatrix} \omega_x \\ \omega_y \\ \omega_z \end{bmatrix} = \boldsymbol{\mathrm{B}}_d^t\begin{bmatrix} \dot\vartheta\sin\psi \\ \dot\psi \\ \dot\vartheta\cos\psi \end{bmatrix}+ \begin{bmatrix} \dot\gamma \\ 0 \\ 0 \end{bmatrix}\\$

展開解得

$\begin{cases} \dot\gamma=\omega_x-\tan\vartheta(\omega_y\cos\gamma-\omega_z\sin\gamma)\\ \dot\psi=\dfrac{1}{\cos\vartheta}(\omega_y\cos\gamma-\omega_z\sin\gamma)\\ \dot\vartheta=\omega_y\sin\gamma+\omega_z\cos\gamma \end{cases}\\$

航跡座標系下

把速度分量由航跡座標系變換到地面座標系下，使用下式

$\begin{bmatrix} \dot x_d \\ \dot y_d \\ \dot z_d \end{bmatrix}=\boldsymbol{\mathrm{B}}_h^d\begin{bmatrix} v \\ 0\\0 \end{bmatrix}\\$

展開解得

$\begin{cases} \dot x_d = v\cos\theta\cos\psi_s\\ \dot y_d = v\sin\theta \\ \dot z_d = -v\cos\theta\sin\psi_s \end{cases}\\$

其他關係式

設風速

$\boldsymbol{w}$

在地面座標系下的分解為

$[w_{xd},\;w_{yd},\;w_{zd}]^\mathrm{T}$

。則空速可表示為

$\begin{bmatrix} u_{xt} \\ u_{yt} \\ u_{zt} \end{bmatrix}=\begin{bmatrix} v_{xt} \\ v_{yt} \\ v_{zt} \end{bmatrix}-\boldsymbol{\mathrm{B}}_d^t\begin{bmatrix} w_{xg} \\ w_{yg} \\ w_{zg} \end{bmatrix}\\$

那麼由空速在氣流座標系和機體座標系下的關係可得

$\begin{bmatrix} u_{xt} \\ u_{yt} \\ u_{zt} \end{bmatrix} = \boldsymbol{\mathrm{B}}_q^t\begin{bmatrix} u \\ 0 \\ 0 \end{bmatrix} = \begin{bmatrix} u\cos\alpha\cos\beta \\ -u\sin\alpha\cos\beta \\ u\sin\beta \end{bmatrix}\\$

由此可得，由機體座標系中的速度分量求氣流角和空速大小的表示式為

$\begin{cases} \alpha=\arctan\left(\dfrac{-u_{yt}}{u_{xt}}\right)\\ \beta=\arcsin\left(\dfrac{u_{zt}}{u}\right)\\ u=\sqrt{{u_{xt}}^2+{u_{yt}}^2+{u_{zt}}^2} \end{cases}\\$

上式左右對

求導可得

$\begin{equation} \begin{aligned} & \frac{1}{\cos^2\alpha}\dot{\alpha}=\frac{\dot{u}_{yt}u_{xt}-u_{yt}\dot{u}_{xt}}{{u_{xt}}^2} \\ & \cos\beta\dot{\beta}=\frac{\dot{u}_{zt}u-u_{zt}\dot{u}}{u^{2}} \\ & u\dot{u}=u_{xt}\dot{u}_{xt}+u_{yt}\dot{u}_{yt}+u_{zt}\dot{u}_{zt} \end{aligned} \end{equation}\\$

則氣流角和空速的導數可以透過下式求得。

$\begin{cases} \dot u=\dfrac{u_{xt}\dot{u}_{xt}+u_{yt}\dot{u}_{yt}+u_{zt}\dot{u}_{zt}}{u} \\ \dot\alpha=\cos^2\alpha\dfrac{\dot{u}_{yt}u_{xt}-u_{yt}\dot{u}_{xt}}{{u_{xt}}^2}\\ \dot\beta=\dfrac{\dot{u}_{zt}u-u_{zt}\dot{u}}{\cos\beta u^{2}} \end{cases}\\$

由航跡系變換到地面系再變換到機體系，和由航跡系變換到氣流系再變換到機體系的變換矩陣應該相同，在

無風

的情況下就可以得出以下關係式。

$\boldsymbol{\mathrm{B}}_q^t\boldsymbol{\mathrm{B}}_h^q=\boldsymbol{\mathrm{B}}_d^t\boldsymbol{\mathrm{B}}_h^d\\$

令兩邊的第 3 行第 1 列元素相等，可得

$\sin\beta=[\sin\gamma\sin\vartheta\cos(\psi-\psi_s)+\cos\gamma\sin(\psi-\psi_s)]\cos\theta-\sin\gamma\cos\vartheta\sin\theta\\$

令第 2 行第 1 列元素相等，可得

$\sin\alpha=\{[\cos\gamma\sin\vartheta\cos(\psi-\psi_s)-\sin\gamma\sin(\psi-\psi_s)]\cos\theta-\cos\gamma\cos\vartheta\sin\theta\}/\cos\beta\\$

同理，由

$\boldsymbol{\mathrm{B}}_h^d\boldsymbol{\mathrm{B}}_q^h=\boldsymbol{\mathrm{B}}_t^d\boldsymbol{\mathrm{B}}_q^t$

，令第 2 行第 3 列元素相等，可得

$\sin\gamma_s=[\cos\alpha\sin\beta\sin\vartheta-\sin\alpha\sin\beta\cos\gamma\cos\vartheta+\cos\beta\sin\gamma\cos\vartheta]/\cos\theta\\$

結語

至此，我們後面建模要用到的方程就推導完畢了。一般的飛行器模型包含 12 個狀態量，根據所使用的方程不同，具體是哪 12 個量作為狀態量也存有差別。

值得一提的是，作為狀態量，一定是要先求出其微分，然後對微分做積分得到狀態量的值。Simulink 中對系統進行分析時會自動把積分模組的輸入當做狀態量，這一點在建模時一定要注意。比如利用

其他關係式

直接求出氣流角和空速值，就不能將其作為狀態量，因為不是透過積分得到的。

在《飛行器運動方程》一書中，介紹了目前使用的兩種主流運動模型結構，分別是 H-T 和 T-T。第一個字母代表質心平移動力學方程是在哪個座標系下建立的。第二個字母則代表繞質心轉動的動力學方程在哪個座標系下建立。而在實際運用中，T-T 在有風的情況下形式更簡單，而 H-T 則十分繁雜。因此我們後續只使用 T-T 來建立我們的模型。其對應的 12 個狀態量為

$[v_{xt},v_{yt},v_{zt},\gamma,\psi,\vartheta,\omega_{x},\omega_{y},\omega_{z},x_d,y_d,z_d]$