如何获得飞机运动方程

tuchooy

一、什么是惯性坐标系二、飞机主要使用的坐标系三、非惯性坐标系下的牛顿第二定律四、施加在飞机上的力和力矩
前面一篇文章J Pan：飞机是怎么飞起来的 引起了很多人的兴趣，在私信和留言里面大家讨论了很多其实是关于飞机控制的问题，说明大家对航空知识的热情还是很高的，这自然是好事情，因为这个事情在美国，NASA每年是需要花费大量资金来做的。
说到飞机控制，一般称之为飞行控制，主要解决飞机的稳定性和操纵性问题，而这两者往往是矛盾的，是需要权衡的，比如军机就要机动性好，也就是操纵性强，而民机则要求稳定性好。要进行飞行控制，首要解决的是飞机飞行时的建模问题，这一部分相对比较繁琐，也是很多人望而却步的一部分，然后又是一切工作的基础，没有模型，谈何控制呢？
飞机建模，要从坐标系说起，就是要解决飞机在三维空间中的运动是如何描述的，本质就是如何在非惯性坐标系下建立牛顿第二定律。
<hr/>一、什么是惯性坐标系

三维空间中，忽略飞机的弹性变形，飞机可以看成是一个刚体，其运动可以用两个概念来表示：平移和旋转。
平移，可以用向量来表示，因为它不仅有大小，还有方向。为了描述向量，应该先确定一个具体的坐标系，明确该坐标系的线性基 [\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3] 后，才能够确定一个向量  \mathbf{\vec{a}} 在该坐标系下的坐标：
\mathbf{\vec{a}}=[\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3] \begin{bmatrix} a_1 \\ a_2 \\ a_3 \end{bmatrix}  
\begin{bmatrix} a_1 & a_2 & a_3 \end{bmatrix} ' 称之为坐标值，注意坐标值和向量并不是完全等价的，这是我们常常忽略的部分，坐标值和坐标系合在一起才能描述向量的完整信息。话句话说，向量是客观存在的，它是不随坐标系变化的，但是它可以通过某个坐标系及其在该坐标系下的坐标值进行描述。注意坐标系的选取可以是任意的，比如我们也可以选择坐标系 [\mathbf{e}'_1, \mathbf{e}'_2, \mathbf{e}'_3] ，此时：
[\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3] \begin{bmatrix} a_1 \\ a_2 \\ a_3 \end{bmatrix} = [\mathbf{e}_1', \mathbf{e}_2', \mathbf{e}_3'] \begin{bmatrix} a_1' \\ a_2' \\ a_3' \end{bmatrix}
显然，同一个向量，在不同的坐标系下，其坐标值是不同的，向量坐标的具体取值，和向量本身和选取的坐标系相关。那是不是就意味着我们可以随意选择坐标系呢？——答案是否定的，因为有些坐标系，有一些特殊的性质，我们可以加以利用，最常用的坐标系是——惯性坐标系。
简单点来说，如果有一个坐标系，不受外力时，一切物体总保持与参考系的匀速直线运动状态或相对静止状态，那它就是惯性坐标系，反之则不是。这是牛顿第一定律确定的事。
比如对于我们碰到的大多数问题而言，地球就可以看成是一个近似理想的惯性坐标系，相对于地球静止或者做匀速运动的参考系，也可以认为是惯性坐标系。在惯性坐标系，牛顿第二定律可以表示为：
F=ma
如果一个坐标系相对于惯性坐标系有加速度，则该坐标系不是惯性坐标系。比如，你在高铁的桌面上放一个小球，假设桌面是完全光滑的，高铁加速时，高铁上的人会发现小球运动。很容易发现小球运动状态的改变并非受力的原因。所以，这时的高铁就是非惯性系。


这个时候牛顿第二定律就不是 F=ma 了，因为没有力，也能产生加速度——因为我们采用的坐标系（加速的高铁）相对于惯性坐标系产生了加速度，这样在非惯性系里面，需要对牛顿第二定律 F=ma 的形式进行修正。
如果能选择惯性坐标系作为参考坐标系，那自然是最好的，我们可以直接利用牛二定理的简洁形式，但是对于飞机而言，事情要更复杂一点，请看下图：


飞机的运动姿态可以加速的高铁复杂多了，因此坐在飞机上的飞行员，天然会以自己的视角作为参考系，飞行员可以看成和飞机固连在一起，因此此时飞机的机体作为了一个参考坐标系，而这个参考系，很显然是非惯性坐标系。
当然不仅仅于此了，我们知道，飞机的升力来自于空气，也就是说，研究气动力的时候，最好以气流的方向作为参考系，这也是一个非惯性参考系。
实际上还有其他坐标系需要用，我们需要详细了解其中几个比较重要的。
<hr/>二、飞机主要使用的坐标系

最常用的坐标系当然是地面坐标系（Ground axis），固定于大地表面，原点为海平面或地面某点，轴OZ铅锤向下，轴OX和OY在水平面内，只要满足右手直角坐标系要求，在大地上的方向可以任意选取。
另一个重要的坐标系是机体坐标系（Body axis），与飞机固连在一起，原点O位于飞机的质心，轴OXb平行于机身轴段或机翼平均气动弦，指向机头方向，OYb沿飞机的对称方向，OZb满足右手定则，示意图如下：


机体坐标系和地面坐标系的关系如下：

https://zippy.gfycat.com/
https://www.zhihu.com/video/1128451143017402368
其中：
\Psi ——偏航角（angle of yaw）：
\Theta ——俯仰角（angle of pitch）；
\Phi ——滚转角（angle of roll）：
也就是，地面坐标系通过 \Psi,\Theta,\Phi 三个角度的旋转（顺序关系不能乱哦）就可以变换到机体坐标系。
有了机体坐标系，我们还要衍生出另外两个常用且重要的坐标系，即气流坐标系（Wind axis）和稳定性坐标系（Stability axis）。
先来说一下稳定性坐标系，如下图所示：


它是由机体坐标系沿Y轴旋转 \alpha 角度（攻角或迎角）得到的，当不考虑侧风的时候，即空气沿飞机机体的X轴线运动，气动力在该坐标系下很容易的到。
对于有侧风的情况，要略微再复杂一点，即将稳定性坐标系沿该坐标系的Z轴旋转 \beta 角度（侧滑角）得到，示意图如下：


其中：
\alpha ——攻角或迎角（angle of attack or angle of incidence）；
\beta ——侧滑角（angle of side slip ）；
至此，最重要的几个角度定义基本齐活了。
<hr/>说到这，我们就可以稍微总结一下：我们要对飞机建模，就要有坐标系；坐标系有惯性系和非惯性系两种，而牛顿第二定律 F=ma 只在惯性系成立；飞机机体是一个非惯性坐标系，那在非惯性坐标系下，牛顿第二定律长什么样子呢？
三、非惯性坐标系下的牛顿第二定律

我们知道，在惯性坐标系力里面，牛顿第二定律为：
F=ma
这是高中知识，到了大学之后，我们一般写成更通用的形式：
\color{deeppink}{\frac{d(m\vec{\mathbf{v}})}{dt}=F}
\color{deeppink}{\frac{d(\vec{\mathbf{H}})}{dt}=M}
即动量的变化率是力，角动量的变化率是力矩。现在我们就要看一下这两个公式在非惯性坐标系下是什么样的。
3.1 非惯性坐标系下的力平衡方程
注意，以下讨论都是基于机体坐标系进行的，假定机体坐标系的线性基在 x,y,z 三个方向的线性基为[\mathbf{i}, \mathbf{j}, \mathbf{k}] （该坐标系在绕地面坐标系旋转），则任意一个在惯性坐标系下的向量  \mathbf{\vec{a}} 在机体坐标系下可表示为：
\mathbf{\vec{a}} =[\mathbf{i}, \mathbf{j}, \mathbf{k}] \begin{bmatrix} a_x \\ a_y \\ a_z \end{bmatrix} =(\mathbf{i}a_x+\mathbf{j}a_y+ \mathbf{k}a_z )
则向量  \mathbf{\vec{a}} 随时间的变化为：
\frac{d\mathbf{\vec{a}}}{dt}=\frac{d(\mathbf{i}a_x+\mathbf{j}a_y+ \mathbf{k}a_z) }{dt}
对上式的右半部分进行微分，注意 a_x,a_y,a_z 以及 \mathbf{i}, \mathbf{j}, \mathbf{k} 都是变量，可以得到：
\begin{equation} \begin{aligned}  \frac{d\mathbf{\vec{a}}}{dt} &= \underbrace{(\frac{\partial{a_x}}{\partial t}\mathbf{i}+ \frac{\partial{a_y}}{\partial t}\mathbf{j}+ \frac{\partial{a_z}}{\partial t}\mathbf{k})}_{\frac{\partial{\mathbf{\dot{\vec a}}}}{\partial t}}\\ &+\underbrace{(a_x\frac{\partial{\mathbf{i}}}{\partial t}+a_y\frac{\partial{\mathbf{j}}}{\partial t} +a_z\frac{\partial{\mathbf{k}}}{\partial t})}_{\bm{\omega}\times \mathbf{\vec a}}\\ &=\color{deeppink}{\frac{\partial{\mathbf{\dot{\vec a}}}}{\partial t}+\bm{\omega}\times \mathbf{\vec a}} \end{aligned} \end{equation}   
其中
\bm{\omega}=P\mathbf{i}+Q\mathbf{j}+R\mathbf{k}
\bm{\omega} 表示机体（body）坐标系绕惯性（inertial）坐标系旋转矢量， P,Q,R 分别为机体坐标系绕地面坐标系的滚转角速率，俯仰角速率以及偏航角速率。至于
(a_x\frac{\partial{\mathbf{i}}}{\partial t}+a_y\frac{\partial{\mathbf{j}}}{\partial t} +a_z\frac{\partial{\mathbf{k}}}{\partial t})=\bm{\omega}\times \mathbf{\vec a}
是怎么得到的，大家感兴趣可以翻阅一下刚体力学的教材，可以通过基本定义得到，这里不详细推导。将矢量形式写成矩阵形式为：
\bm{\omega}\times \mathbf{\vec a}=\begin{Vmatrix} \mathbf{i} & \mathbf{j} & \mathbf{k} \\ P & Q &R \\ a_x & a_y & a_z \\ \end{Vmatrix}
有了前面的铺垫，后面就方便多了，假设飞机的速度矢量（在惯性坐标系）为 \mathbf{\vec v} ，
\mathbf{\vec{v}} =[\mathbf{i}, \mathbf{j}, \mathbf{k}] \begin{bmatrix} U \\ V \\ W \end{bmatrix}
U,V,W 为速度矢量在机体坐标系下的坐标值。则速度矢量的可在机体坐标系表示为：
\color{deeppink}{\frac{d \mathbf{\vec v}}{dt}=\frac{\partial{\mathbf{\vec v}}}{\partial t}+\bm{\omega}\times \mathbf{\vec v}}  
所以：
\frac{d \mathbf{\vec v}}{dt}=\begin{bmatrix} \dot U \\ \dot V\\ \dot W \end{bmatrix}+\begin{Vmatrix} \mathbf{i} & \mathbf{j} & \mathbf{k} \\ P & Q &R \\ U & V & W \\ \end{Vmatrix}  
在惯性坐标系，满足牛顿第三定律：
m\frac{d \mathbf{\vec v}}{dt}=\mathbf{F}
同时，可以将在惯性坐标系的里转换到机体坐标系：
\mathbf{F}=[\mathbf{i}, \mathbf{j}, \mathbf{k}] \begin{bmatrix} F_x \\ F_y \\ F_z \end{bmatrix}
综合以上各式：
m\begin{bmatrix} \dot U \\ \dot V\\ \dot W \end{bmatrix}+m\begin{Vmatrix} \mathbf{i} & \mathbf{j} & \mathbf{k} \\ P & Q &R\\ U & V & W\\ \end{Vmatrix} =\begin{bmatrix} F_x \\ F_y \\ F_z \end{bmatrix}
展开即可以获得：
\color{deeppink}{\left\{ \begin{array}{lcl}    m(\dot U+QW-RV)=F_x  \\   m(\dot V+RU-PW)=F_y\\ m(\dot W+PV-QU)=F_z  \end{array} \right.}

3.2 非惯性坐标系下的力矩平衡方程
有了前面的铺垫就简单多了，在惯性坐标系下：
\color{deeppink}{\frac{d(\vec{\mathbf{H}})}{dt}=M}
其中角动量的定义为：
\vec{\mathbf{H}}=\int\vec{\mathbf{r}}\times \vec{\mathbf{v}}dm
\mathbf{\vec{r}}  为飞机中任意一微元在地面参考系中的矢量，其可在机体坐标系中表示为：
\mathbf{\vec{r}} =[\mathbf{i}, \mathbf{j}, \mathbf{k}] \begin{bmatrix} x \\ y \\ z \end{bmatrix}
则速度矢量可表示为：
\vec{\mathbf{v}}=\frac{d \vec{\mathbf{r}}}{dt}=\frac{\partial{\mathbf{\vec r}}}{\partial t}+\bm{\omega}\times \mathbf{\vec r}
于是可以的到地面参考系中的角动量方程：
\vec{\mathbf{H}}=\begin{bmatrix}  I_{xx}& -I_{xy} & -I_{xz} \\ - I_{xy}& I_{yy} & -I_{yz} \\  -I_{xz}& -I_{yz} & I_{zz}  \end{bmatrix} \begin{bmatrix} P \\ Q \\R \end{bmatrix}
其中 I_{xy} 等为飞机在机体坐标系惯性矩，飞机各平面的定义方式如下：


由于飞机沿 y 轴的对称性，有
I_{xy}=\int xydm=0
I_{yz}=\int yzdm=0
所以可以简化为
\vec{\mathbf{H}}=\begin{bmatrix}  I_{xx}& 0 & -I_{xz} \\ 0& I_{yy} & 0 \\  -I_{xz}& 0 & I_{zz}  \end{bmatrix} \begin{bmatrix} P \\ Q \\R \end{bmatrix}
由机体坐标系和地面坐标系的转换关系，可以得到：
\color{deeppink}{\frac{d \vec{\mathbf{H}}}{dt}=\frac{\partial{\mathbf{\vec H}}}{\partial t}+\bm{\omega}\times \mathbf{\vec H}}
展开后为：
\color{deeppink}{\dot P I_{xx}+QR(I_{zz}-I_{yy})-(\dot R+PQ)I_{xz}=M_x  }   
\color{deeppink}{\dot Q I_{yy}-PR(I_{zz}-I_{xx})+(P^2-R^2)I_{xz}=M_y }  
\color{deeppink}{\underbrace{\dot R I_{zz}}_{\text{accelaration}}+\underbrace{PQ(I_{yy}-I_{xx})}_{\text{precession}}+\underbrace{(QR-\dot P)I_{xz}}_{\text{coupling}}=M_z }  
其中 L,M,N 分别为滚转力矩、俯仰力矩和航向力矩。注意，由于历史原因， L 有时候表示升力（lift），又是表示滚转力矩，可通过上下文区分。
可见，力矩平衡方程主要组成包括三部分，即角加速度项（angular acceleration terms），角速率进动项（gyro precession terms）和耦合项（coupling terms）。下面我们分别来解释一下。
角加速度项是最容易理解的，以滚转力矩方程为例，忽略进动项和耦合项，有：
\dot P I_{xx}=L
啥意思呢？——如果有一个滚转力矩加在飞机上，比如副翼有个偏转，加速度项表明飞机回产生一个滚转角加速度，转动惯量 I_{xx} 越大，获得的角加速度就越小，这是我们非常熟悉的概念。
再来看看进动项，为简单起见，我们还是先忽略加速度项和耦合项，同时假定 I_{zz} 非常小，也可以忽略，这时候滚转力矩平衡方程就简化为：
-I_{yy}QR=L
这样还是看不出来什么，再稍微变个形：
-R\underbrace{(I_{yy}Q)}_{\text{momentum(pitch)}}=L
式中 I_{yy}Q 表示俯仰轴的角动量（转动惯量乘以角速率），假如飞行员给副翼一个偏转指令，产生滚转力矩 L ，从这个方程可以看出，飞机在出现翻滚的同时，还会出现一个航向的运动 R ，这就是我们常说的荷兰滚。
耦合项描述了飞机的惯性耦合倾向。说着有点拗口，还是拿公式分析一下，忽进动项，假定零偏航速率（R），施加的俯仰力矩（M）也为零，有：
P^2I_{xz}=-I_{yy} \dot Q
这个方程表示啥意思呢？——飞行员在做翻滚运动时（P），飞机会出现一个俯仰角加速度，出现非期望的抬头，如果不进行压杆控制，就可能会出现飞行事故。在上世纪40-50年代，人们对于耦合项的理解还不是很深入时，就出现过几起这样的事故。
<hr/>四、施加在飞机上的力和力矩

前面我们得到了飞机的六自由度运动方程，这些方程之间是相互耦合的，给我们分析飞机的行为带来了很大的困难，通过对这六个方程深入研究，我们发现，可以将其分成两组，这两组耦合性相对较小——一组我们称之为纵向方程组，一组称之为横航向方程组。
纵向方程组为一个旋转，两个平动：
\color{deeppink}{\left\{ \begin{array}{lcl}    m(\dot U+QW-RV)=F_x  \\   \dot Q I_{yy}-PR(I_{zz}-I_{xx})+(P^2-R^2)I_{xz}=M_y \\ m(\dot W+PV-QU)=F_z  \end{array} \right.}
横航向为两个旋转，一个平动：
\color{deeppink}{\left\{ \begin{array}{lcl}    \dot P I_{xx}+QR(I_{zz}-I_{yy})-(\dot R+PQ)I_{xz}=M_x  \\   m(\dot V+RU-PW)=F_y \\ \dot R I_{zz}+PQ(I_{yy}-I_{xx})+(QR-\dot P)I_{xz}=M_z  \end{array} \right.}
好了进行到现在，我们还是只解决了飞机运动方程的左半部分，而右半部分的力及机具还完全没有分析。
我们知道，飞机飞行时，受到的的力主要有三种：重力、气动力以及发动机推力，其中气动力又可以分为升力和阻力。产生的力矩主要有两种：气动力矩和发动机偏转力矩（矢量发动机），先来看看力。
重力显然在地面坐标系描述比较简单；气动力在气体坐标系下描述比较简单；当不考虑侧风时，在稳定性坐标系下描述简单；发动机推力在机体坐标系下描述比较简单。最通用的做法就是把他们都统一到机体坐标系（可通过坐标变换实现），这时飞机受到的力以及力矩为：
纵向：
\color{deeppink}{\left\{ \begin{array}{lcl}     F_x=-mgsin\Theta + (-Dcos\alpha+Lsin\alpha)+Tcos\Phi_T  \\    M_y=M_A+M_T \\ F_z= mgcos \Phi cos\Theta + (-Dsin\alpha-Lcos\alpha)-Tsin\Phi_T \end{array} \right.}
横航向：
\color{deeppink}{\left\{ \begin{array}{lcl}    M_x=M_{Ax}+M_{Tx}  \\    F_y=mgsin\Phi cos\Theta+F_{A_y}+F_{T_y} \\  M_z=M_{Az}+M_{Tz}  \end{array} \right.}
力矩的产生比较复杂，一般通过飞机的吹风或者CFD计算确定，具体处理方式后续线性化的时候再说。
参考文献：

• Thomas R. Yechout. Introduction to Aircraft Flight Mechanics Performance, Static Stability, Dynamic Stability, and Classical Feedback Control. AIAA
  

生锈的手指、

大佬6

wqt0806

期待大佬的后续

kimicheng

这本飞行力学书是我也最喜欢的，可以我影印的上下册，被人拿走一册，我的笔记都没了。

Li—_™

坐标变换那一块怎么看都是欧拉角啊

A...周东

膜！对着书想象角度变换真的头疼[捂脸]

ask终结

大佬终于更新了[赞同]

没有驾照的司机

侧滑角跟滚转角有什么区别啊？[好奇][好奇]

yangyang110

跟我们水下航行器是一样的理论基础。但想问楼主一下，像我们这种在不可压缩流体中缓和运动的航行器，是把水动力展开成惯性力、粘性力的线性、二阶三阶和耦合系数，而航空里是怎么处理空气动力的呢？不可能光是航速攻角和升阻力矩那么简单吧？

jasonhsu

还真是那么粗暴，大多把力展成线性项，最多也就再增加个攻角的二次项。

咏蘅誓唁

好吧[捂脸]可能空气密度比水小得多线性够用了，航速又太快随便一个回转半径都几十千米，而我们水下就是在原地旋转跳跃闭着眼的[小情绪]

wqt0806

硬

仓颉

这是航理吗？

PLA

算是吧

a840177949

哇…以后看来要学

SOPHIL❤Sophie

我想问下关于上面力学的知识，可以推荐相关的书籍学习吗？最近在学四旋翼动力学建模

kccsong

对飞行超感兴趣 谢谢楼主

异乡人

大大 你好，除了加速度项，其他两项是否可以理解为施加力矩后，飞机朝着总体能量最小的方向变化，而耦合项 进动项就是能量最小方向

嗯?

想问地面坐标系变换到机体坐标系时三个角度旋转的顺序关系不能乱呢？
[疑惑]

freejaja

大佬，我有个疑问，就是在非惯性坐标系下的力矩平衡方程这一小节中关于角动量的计算，我试着算了下，发现每个基向量下除了有转动惯量和惯性积部分，还有r的xyz的参量与参量微分的乘积部分，而这部分内容在角动量最终矩阵式子中没有体现出来