北航飞行力学理论与应用课程大作业2014第4组

1

飞行力学理论及应用

大作业

院（系） 自动化科学与电气工程学院 专 业

控制科学与工程 朱付涛 BY1403158

学生姓名 李琛 SY1403520

左京兴 BY1303165

2014年1月4日

飞行力学理论大作业

大作业分工情况：

李琛：飞机运动方程建模 左京兴：小扰动线性化

朱付涛：求解飞机时域响应和附录

第一部分 飞机运动建模

在平面地球假设下，飞机飞行在平面大地上空，不考虑地球曲率和地球自转。 首先列出气动力和气动力矩方程：

?X?qSCx??Y?qSCy ??Z?qSCz?L?qSblC??M?qScmC?N?qSbCn?

其中，X,Y,Z为气动力在飞机机体系的分量，坐标系选为美式坐标系。L,M,N为飞机气动力矩在机体上的分量。q为动压，S为飞机的机翼面积，b为翼展，c为平均气动弦长。

1.1 飞机质心动力学方程

???MrM??2?MrM???M?MrM?和式(5.4.9)f?maC出发，由式(5.1.7)aM?aoM?rM推导飞机的质心动力学方

程。本文推导了两种情况下的质心动力学方程：风轴系中的质心动力学方程和体轴系中的质心动力学方程。

(1) 风轴系中的质心动力学方程

式（5.4.9）在风轴系中写为

fw?macw

(1.1)

式（5.1.7）是一般坐标系中惯性加速度的表示，将其写为地轴系中的表达，此时r?就是质心相对地球的速度，用VE来表示。假定地轴固定于惯性空间，则??0。FE原点的加速度a0就是与地球转动有关的向心加速度，它与g相比可以略去。其中的向心加速度项??r?通常也可略去。另外令r??VE,哥氏加速度为2?EVE，则有

EEaCE?VEE?2?EVE

(1.2)

由于FW相对于FE的角速度为(?W??E)，利用

Lbava?vb??bvb

(1.3)

将式(1.2)转换到风轴系得到

EEaCW?LWEaCE?LWE(VEE?2?EVE)?V?(???)WV?2LWE?VEWWEEWEEEE (1.4)

其中最后一项可以化简为：

EEEEELWE?EVE?(LWE?ELEW)(LWEVEE)??WVW

(1.5)

于是得到下式：

EEEEEEaCW?VW?(?W??E)WVW?2?EVE?VW?(?W??E)WVW

(1.6)

其中

?V?VE????Wxw??pW??pEW?

W??0???W?yw?,?WW??q?,?E??qE? ??0????Wzw??W???W?W??rW????rEW??利用矩阵表示式

?0??za?ya?

??a???za0???xa ??

???ya?xa0??可得到?EW

W和?W的表示为：

??0?rWqW???rEWqEW?

?W?0W??rW0pE?EW0pE?W? ???q?,?W??rWWpW0?????qEWpEW0??将式(1.7)和式(1.9)代入到式(1.6)，并假设大气静止，即W=0，则可得：

??aC?xw??

??a??VE?Cyw????V(rW?rW) ??aC?? ?V(qEzw???W?qW)??将式(1.10)代入到式(1.1)中可得到标量方程为：

??fxw?mV

?fmV(rEyw?W?rW) ??fzw??V(qEW?qW)大气飞行的力矢量由两部分组成：气动反作用力A和重力mg，即

f?A?mg

在风轴系FW中，A的分量可以写成下面的形式： ?XW?

A??D??Txw?w???YW?????C?????T?yw

??ZW??????L???Tzw??上式中，D是阻力，C是侧力，L是升力。

重力的表达如下式所示：

?0?

g???0?V

???g??因此利用风轴系的欧拉角，可得

(1.7)

(1.8)

(1.9)

(1.10)

(1.11)

(1.12)

(1.13)

(1.14)

mgW?mLMVgV

cos?wcos?wcos?wsin?w??m??sin?wsin?wcos?w?cos?wsin?wsin?wsin?wsin?w?cos?wcos?w??cos?wsin?wcos?w?sin?wsin?wcos?wsin?wsin?w?sin?wcos?w??0??0?(1.15)

sin?wcos?w????cos?wcos?w????g???sin?w??sin?w??mg??sin??wcos?w?cos???wcos?w??将式(1.11)、式(1.13)和式(1.16)。代入到式(1.12)中可得如下方程组：

??Txw?D?mgsin?w?mV

?T?C?mgcos?Eywwsin?w?mV(rW?rW) ??Tcos?Ezw?L?mgcos?ww??mV(qW?qW)当忽略地球转动时，rE和qEWW项为零，则有

??Txw?D?mgsin?w?mV

?Tyw?C?mgcos?wsin?w?mVrW ??Tzw?L?mgcos?wcos?w??mVqW由此得到了风轴系下飞机的质心动力学方程(1.18)。 (2) 体轴系中的质心动力学方程

式（5.4.9）在体轴系中可写为

fB?maCB

设FBEB相对于FE的角速度为(???)，则根据式Lbava?vb??bvb及式(1.2)可得：

aEEE?E)EECB?LBEaCE?LBE(VE?2?EVE)?VEB?(?B?BVEB?2LBE?EVE

其中最后一项根据教材中式（4.6.8）?a?Lab?bLba可化简为:

LEEELEEEBE?EVE?(LBE?EEB)(LBELE)??BVB 则式(1.20)可改为：

aCB?VEB?(?B??E)BVEB

式(1.22)矢量的标量分量为

?u?? VE????Wx??p?pEB?????W?,???q?y?B??,?EB???qE?B?vB

??w????Wz????r????E?rB??根据式(1.8)可得

???0?rq?0?rEqEBB?

?B?B??B??r0?p?,?EB???rEB0pE?B? ???qp0?????qEBpEB0??

(1.17)

(1.18)

(1.19)

(1.20)

(1.21)

(1.22)

(1.23)

(1.24)

将式(1.23)和式(1.24)代入到式(1.22)中，并考虑大气静止的情况，即W=0，得到如下方程组：

??aEE ?CBx?u?(q?qB)w?(r?rB)v?av?(r?rE?(p?pECBy?B)uB)w ???aCBz?w?(p?pEEB)v?(q?qB)u类似的可得：

mgB?mLBVgV?cos?cos?cos?sin??sin???

?m??sin?sin?cos??cos?sin?sin?sin?sin??cos?cos?sin?cos???0?0?

??cos?sin?cos??sin?sin?cos?sin?sin??sin?cos?cos?cos????????g????sin???mg??cos?sin??????cos?cos??

考虑到体轴系中的气动力表达式为：

?X? A?B???Y

???Z??并利用式(1.12)可得：

fB?AB?mgB

即

??fBx?X?mgsin?

?fBy?Y?mgcos?sin? ??fBz?Z?mgcos?cos?将式(1.29)和式(1.25)代入到式(1.19)中得到体轴系下的质心动力学方程为：

? ?X?mgsin??m[u?(q?qEB)w?(r?rEB)v]?Y?mgcos?sin??m[v?(r?rEEB)u?(p?pB)w] ??Z?mgcos?cos??m[w?(p?pEEB)v?(q?qB)u]当忽略地球转动时，pEEB,qB,rEB为零，则有：

??X?mgsin??m(u?qw?rv)

?Y?mgcos?sin??m(v?ru?pw) ??Z?mgcos?cos??m(w?pv?qu)1.2飞机转动动力学方程

由式（5.4.12）G?h出发，推导飞机的转动动力学方程（考虑发动机转子产生的角动量）。 仅在体轴系下考虑：G?h，式中

(1.25)

(1.26)

(1.27)

(1.28)

(1.29)

(1.30)

(1.31)

h??R?Vdm

(1.32)

将式(1.32)转换成矩阵符号为：

hI??RIVIdm??RIRIdm

(1.33)

根据式RI?LIB(RB??BRB)则可得h在FB中的分量为

hB?LBIhI??LBIRILIBRBdm??LBIRILIB?BRBdm

(1.34)

根据教材中式（4.7.4），有RB?LBIRILIB，则上式可以化简为：

hB??RBRBdm??RB?BRBdm

(1.35)

利用?R??R?可将上式中第二项改写为

?R?BBRBdm???RBRB?Bdm?JB?B

(1.36)

则式(1.35)可重新写为：

hB??RBRBdm?JB?B

(1.37)

在以上分析的基础上，设两个体轴系FB1和FB2，并考虑转子产生的角动量，FB2相对于FB1有一角速度??，则对这两个参考系的角动量为：

ri?hB?RBRBdm?JB?B??hB?111111?i ?ri?hB2??RB2RB2dm?JB2(?B1???)??hB2i? (1.38)

令RB2RB2dm?hB，则可将其转换成下列形式：

2?eehB??TB2B1RB1TB1B2RB2dm?TB2B1?RB1(RB1???B1RB1)dm2

??TB2B1???RB1RB1dm??RB1RB1??B1dm??TB2B1??RB1RB1dm?JB1??B1??? (1.39)

当hB为零，即RB1RB1dm?JB1??B1?0，则得到了平均轴系下的关系式：

2e?

rihB?JB?B??hB

i (1.40)

下面再回到(5.4.12)式，它在参考系FI中的形式为：

GI?hI

(1.41)

那么转换到体轴系中为：

GB?LBIGI?hB??BhB

(1.42)

其中

GB为：

?L??

GB??M????N?? (1.43)

当采用平均轴系时，将式(1.40)代入到式(1.42)中得到：

ririGB?JB?B?JB?B??BJB?B??hB???BhB

ii (1.44)

对于刚体，则上式中可以忽略JB，并进行标量展开为：

?L??Ix?M????I???xy??N?????Izx?IxyIy?Iyz?Izx??p??0?rq??Ix??????I?Iyz??q?r0?p?????xyIz??r?????qp0???????Izx?IxyIy?Iyz??ri?ri?hhx???x??Izx??p???ii???????riri?Iyz??q????hy???B??hy? (1.45)

?i??i???Iz?rr????hri???hi??zz?????i??i?将上式写成方程组的形式为：

?ririri22L?Ip?I(q?r)?I(r?pq)?I(q?rp)?(I?I)qr?rh?qh?h???xyzzxxyyzyzx?iii??rrr22?M?Iyq?Izx(r?p)?Ixy(p?qr)?Iyz(r?pq)?(Iz?Ix)rp?r?hxi?p?hzi??hyiiii??N?Ir?I(p2?q2)?I(q?rp)?I(p?qr)?(I?I)pq?qhri?phri?hri???zxyyzzxxyxyz?iii?

通常Cxz是飞行器的对称面，此时Ixy?Iyz?0,则上式可以简化为

(1.46)

?riririL?Ip?I(r?pq)?(I?I)qr?rh?qh?h???xzxyzyzx?iii??rrr22?M?Iyq?Izx(r?p)?(Iz?Ix)rp?r?hxi?p?hzi??hyi

iii??N?Ir?I(p?qr)?(I?I)pq?qhri?phri?hri???zzxxyxyz?iii?列出整个方程的标量展开式没有什么好处，对于忽略?的且没有转子项的局限情况，即对于刚体，式（1.46）标量方程为

?L?Ixp?Izx(r?pq)?(Iy?Iz)qr?22?M?Iyq?Izx(r?p)?(Iz?Ix)rp ??N?Izr?Izx(p?qr)?(Ix?Iy)pq

(1.47)

1.3 飞机质心运动学方程

由速度矢量三角形VE?V+W出发，推导飞机的质心运动学方程。飞机的质心运动学方程在风轴系和体轴系展开。

(1) 风轴系的质心运动学方程

E由VE?V+W可写出风轴系中的方程VW?VW?WW,进而可以得到牵连垂直坐标系中的分量为：

VVE?TVW(VW?WW)

(1.48)

式中

?V??Wxw??,W??W?

VW??0??W?yw????0???Wzw??E?Vxv??cos?Wcos?W?E???Vyv???cos?Wsin?WE??Vzv??sin?W??? (1.49)

将式(1.48)展开成标量形式为：

sin?Wsin?Wcos?W?cos?Wsin?Wsin?Wsin?Wsin?W?cos?Wcos?Wsin?Wcos?Wcos?Wsin?Wcos?W?sin?Wsin?W??V?Wxw??W? (1.50) cos?Wsin?Wsin?W?sin?Wcos?W?yw?????cos?Wcos?W??Wzw??当忽略地球曲率时（?V?0），即可认为FV平行于FE，并且假设大气相对于地球静止，即W=0，而质心的位置坐标

?xE,yE,zE?可写成下面的形式：

?xE?Vcos?Wcos?W??yE?Vcos?Wsin?W ?z??Vsin?W?E (1.51)

上式便是风轴系中飞机质心运动学方程。

(2) 体轴系中运动学方程

由VE?V+W可写出体轴系中的方程VBE?VB?WB,进而得到牵连垂直坐标系中的分量：

VVE?TVB(VB?WB)

(1.52)

式中

?u??Wx??,W??W?

VB??v??B?y????w???Wz?? (1.53)

当忽略地球曲率时（?V?0），即可认为FV平行于FE，并且假设大气相对于地球静止，即W=0，而质心的位置坐标

?xE,yE,zE?可写成下面的形式：

?xE??u??y??T?,?,??v?

????E?VB????zE???w?? (1.54)

上式便是体轴系中的飞机质心运动学方程。

1.4飞机转动运动学方程

BV由(5.2.5)式????i??j????k??出发，推导飞机的转动运动学方程。

(1) 体轴系转动运动学方程 由式(5.2.5)可得

?1??0???T(?)??1??T(?)T(?)????V?i??j3??k2????0y??x??x???0???0???0??0??? ??1?? (1.55)

所以

?p?????V???q????sin??

?????cos???cos?sin???

??r?????cos?cos???sin???将上式变换得到欧拉角的变化率的表示为：

???p? ?qsin?tan??rcos?tan????qcos??rsin? ????(qsin??rcos?)sec?上式就是体轴系中的飞机转动运动学方程。 (2) 风轴系转动运动学方程

与体轴系中的推导过程完全相同，最终得到的转动运动学方程为：

???W?PW?QWsin?Wtan?W?RWcos?Wtan?W

??W?QWcos?W?RWsin? ?W??W?(QWsin?W?RWcos?W)sec?W1.5六自由度全量方程组

当无风时( W = 0 )，对于具有对称面的刚体飞机，建立飞机的6自由度全量运动方程组 飞机的质心动力学方程：

??X?mgsin??m(u?qw?rv)

?Y?mgcos?sin??m(v?ru?pw) ??Z?mgcos?cos??m(w?pv?qu)对于具有对称面的刚体飞机的转动动力学方程为：

??L?Ixp?Izx(r?pq)?(Iy?Iz)qr

?M?I2?p2yq?Izx(r)?(Iz?Ix)rp ??N?Izr?Izx(p?qr)?(Ix?Iy)pq飞机的质心运动学方程为：

? ?xE??u??y?E?T?,?,??? z?VB???v??

?E????w??飞机的转动运动学方程为：

????p?qsin?tan??rcos?tan?

???qcos??rsin? ????(qsin??rcos?)sec?第二部分 小扰动线性化

设基准运动为对称定常直线水平飞行，对飞机全量运动方程组进行小扰动线性化处理。

(1.56)

(1.57)

(1.58)

(1.59)

(1.60)

(1.61)

(1.62)

2.1飞机全量运动方程组

飞机的全量运动方程组如下所示，在此基础上进行小扰动线性化。

?X?mgsin??m(u?qw?rv)??Y?mgcos?sin??m(v?ru?pw)?Z?mgcos?cos??m(w?pv?qu)??L?Ixp?Izx(r?pq)?(Iy?Iz)qr?22?M?Iyq?Izx(r?p)?(Iz?Ix)rp??N?Izr?Izx(p?qr)?(Ix?Iy)pq ???p?qsin?tan??rcos?tan?????qcos??rsin?????(qsin??rcos?)sec??xE??u??y??T?,?,??v?????E?VB????zE???w??

由所给的F-16气动参数表，可以得到气动力和气动力矩的表达式为：

1?V2S(Cx?Cxq?q)?T21Y??V2S(Cy?Cyr?r?Cyp?p)21Z??V2S(Cz?Czq?q)2

1L??V2Sc(Cl?Clr?r?Clp?p?Cl?a??a?Cl?r??r)21M??V2Sc(Cm?Cmq?q)21N??V2Sc(Cn?Cnr?r?Cnp?p?Cn?a??a?Cn?r??r)2X?气动力系数需要按如下方法通过插值求出：

Cx?Cx(?,?e)Cxq?Cxq(?)Cyr?Cyr(?)Czq?Czq(?)Cl?a?Cl?a(?,?)Cn?a?Cn?a(?,?)Cl?r?Cl?r(?,?)Cn?r?Cn?r(?,?)Cyp?Cyp(?)

Cy??0.02??0.021(?a/20)?0.086(?r/30)Cz?A(?)(1?(beta/57.3)^2)?0.19(?e/25)Cl?Cl(?,?)Cm?Cm(?,?e)Cn?Cn(?,?)

Clr?Clr(?)Cmq?Cmq(?)Cnr?Cnr(?)Clp?Clp(?)Cnp?Cnp(?)以上为飞机的全量运动方程组。

2.2小扰动线性化方程组

设基准运动为对称定常直线水平飞行，对刚体飞机全量运动方程组进行小扰动线性化处理。

对如下纵向非线性方程

X?mgsin??m(u?qw?rv)Z?mgcos?cos??m(w?pv?qu)M?Iyq?Izx(r2?p2)?(Iz?Ix)rp (2.1)

??qcos??rsin?进行小扰动线性化，有

X0??X?mg(sin?0???cos?0)?m?u

Z0??Z?mgcos(??)cos(?0???)?m(w?qu)M0??M?Iyq (2.2)

?0????qcos??r??q其中上述四个方程的基准运动方程为：

X0?mgsin?0?0

Z0?mgcos?0?0M0?0 (2.3)

?0?0用式(2.2)减去式(2.3)并化简得到：

?X??u??g??cos?0?m??w??Z?g??sin??uq00?m ??M?q??Iy?????q? (2.4)

然后将所给的力和力矩的导数形式

??X?Xu??u?Xw?w?Xq?q?X??T???T?X??e???e?? ??Z?Zu??u?Zw?w?Zq?q?Z??T???T?Z??e???e????M?Mu??u?Mw?w?Mq?q?M??T???T?M??e???e (2.5)

代入到式(2.4)中得到：

X??eX??TXq?XuXw?u??u?w?q?gcos?????????e?0Tmmmmm?Z??eZ??T??Zq?ZuZww??u?w??uq?gsin?????????e??0?0Tmmmm?m? ??ZZMMM?q?u?u?ww?qq???T??T???e??eIyIyIyIyIy??????q (2.6)

T然后以x?[?uwq??]为状态变量，以c?[??T??e]T为控制变量，可将式(2.6)线化成x?Ax?Bc的

形式，即

?Xu?m??u???w??Zu????m?q?????Mu?????I?y??0XwmZwmMwIy0XqmZqmMqIy1?u0?X??T??gcos?0?????u??m?Z??T??gsin?0??w??????q???m????M??T0??????I?y???0??0X??e??m?Z??e??????m??T?

??M??e??e?Iy??0?? (2.7)

以上是纵向方程组的线化形式，同理可得到横侧向方程组（其中状态变量x?[vpr?]T，控制变量

c?[?a?r]T）的线化形式为：

?Yv?m?v???p??IzLv?IzxNv2????IxIz?Izx?r?????IzxLv?IxNv????II?I2?xzzx?0?YpmIzLp?IzxNp2IxIz?IzxYr?u0mIzLr?IzxNr2IxIz?IzxIzxLr?IxNr2IxIz?Izxtan?0IzxLp?IxNp2IxIz?Izx1?Y?a?gcos?0??m??v??IL?IN?0??p??z?azx2?a????IxIz?Izx??r??????IzxL?a?IxN?a0?????2II?Ixzzx??0?0???Y?r??m?IzL?r?IzxN?r?2???a?(2.8) IxIz?Izx???r??IzxL?r?IxN?r??2?IxIz?Izx?0??

2.3求气动导数

在第二部分中已经建立了分立的小扰动线性化方程，式子中涉及到了无因次量的导数，需要求出这些导数的计算式，计算过程如下，此过程中w?u?/r2d,v?u?/r2d。要建立小扰动线性化方程，在纵向，需要求出气动力和气动力矩关于[uwq导数。

(1) 轴向力X的导数

导数是在基准状态下求得的，所以有p?q?r??a??r?0 轴向力X??T?e]的导数，在横侧向，需要求出气动力和气动力矩关于[vpr?a?r]的

1?V2SCx?T的导数有Xu、Xw、Xq、X??e、X??T： 2??X???X???V0???X???Cx???V0?1??Cx?1Xu?????????uSC??uS??????????????u0SCwsin?0??u0SCxu0x?0?2 ??u0???V0???u0???Cx???u0???u0?2??u0?其中，

X?mgsin?0?Cx?X另外有：

11?u02S?Cwsin?0?Cw?mg?u02S22?0??Cxu?0

1?u0SCx??r2d2?Cx1Cxq??Xq??u0ScCxq?(q/(2u0/c))4w?u?/r2d?Xw?X??eX??T12??u0SCx?e2?T???T

(2) 法向力Z的导数

轴向力Z?1?V2SCz的导数有Zu、Zw、Zq、Z??e、Z??T： 21Zu???Su0Cwcos?0??u0SCzuCzu?021w?u?/r2d?Zw??u0SCz??r2d2?Cz1Czq??Zq??u0ScCzq

?(q/(2u0/c))4Z??e?Z??T12?u0SCz?e2?0Cz?e??0.1925(3) 力矩M的导数

力矩M?1?V2ScCm的导数有Mu、Mw、Mq、M??e、M??T： 21Mu??u0ScCm??u0ScCmuCmu?021w?u?/r2d?Mw??u0ScCm??r2d2?Cm1Cmq??Mq??u0Sc2Cmq

?(q/(2u0/c))4M??e?M??T12?u0ScCm?e2?01?V2SCy的导数有Yv、Yp、Yr、Y??、Y?r： 2 (4) 侧向力Y的导数

侧向力Y?1?u0SCy??r2d2?Cy1Cyp??Yp??u0SbCyp?(p/(2u0/b))4v?u?/r2d?Yv?Cyr?Y????Cy?(r/(2u0/b))?Yr?1?u0SbCyr 412?u0SCy??212Y?r??u0SCy?r2 (5) 力矩L的导数

力矩L?Cy???Cy?r0.021200.086?301?V2SbCl的导数有Lv、Lp、Lr、L??、L?r： 21v?u?/r2d?Lv??u0SbCl??r2d2?Cl1Clp??Lp??u0Sb2Clp?(p/(2u0/b))4Clr?L????Cl1?Lr??u0Sb2Clr

?(r/(2u0/b))412?u0SbCl??212L?r??u0SbCl?r2(6) 力矩N的导数

力矩N?1?V2SbCn的导数有Nv、Np、Nr、N??、N?r： 21?u0SbCn??r2d21Np??u0Sb2Cnp41Nr??u0Sb2Cnr

412N????u0SbCn??212N?r??u0SbCn?r2Nv?第三部分 飞机时域响应求解

已知某飞机数据，建立描述飞机运动特性的全量运动方程。飞机的惯性数据如下表所示：

表1 飞机的惯性数据

重量(lbs) W 25000 重量(kg) W Ix Ix 9496 转动惯量(slug-ft2) Iy 55814 Iy 75673.6077 Iz 63100 Iz Ixz 982 Ixz 转动惯量(kg-m2) 85552.0953 1331.4130 11339.8093 12874.8446 飞机机翼几何尺寸如下表所示：

表2 飞机的几何尺寸 翼展b(ft) 30 9.144 机翼面积S(ft2) 300 27.8709 平均气动弦c(ft) 11.32 3.4503 翼展b(m) 机翼面积S(m2) 平均气动弦c(m) 由于不同的大气密度计算公式算出的大气密度差别很大，所以采用常用的大气模型，1976年美国标准大气

模型USSA76，该模型以地理纬度45°32′33″地区海平面为基准，全年实际大气参数的统计平均值为标准大气参数。大气基准值为：温度T0?288.15K，密度?0?1.225kg/m3，声速a0?340.294m/s。USSA适用高度范围为0km1000km，并分成低层和高层部分。低层高度为0km86km，在这里只介绍低层大气。

低层大气空气分子量为常数，大气是完全气体，处在静平衡状态，其位势高度为hp?1g0?gdh，重力加速

02h度为g?g0?R(R?h)?，h为海拔，g0?9.80665ms/2，与g0对应的地球半径为R?6356766m。于是，

hp?Rh(R?h)?h(1?h/R)。

大气参数分段计算公式： (1). 0km?h?11.0191km

?W?1?hp44.3308? ?T?288.15W????W4.25590?(2). 11.0191km?h?20.0631km

?W?exp?(14.9647?hp)6.3416??? ?T?216.650???1.5898?10?1?W0??h?20.0631km的大气参数在此不列出。声速a?20.0468T。

3.1飞机全量运动方程组

?X?mgsin??m(u?qw?rv)??Y?mgcos?sin??m(v?ru?pw)?Z?mgcos?cos??m(w?pv?qu)??L?Ixp?Izx(r?pq)?(Iy?Iz)qr?22?M?Iyq?Izx(r?p)?(Iz?Ix)rp??N?Izr?Izx(p?qr)?(Ix?Iy)pq ???p?qsin?tan??rcos?tan?????qcos??rsin?????(qsin??rcos?)sec??xE??u??y??T?,?,??v?????E?VB????zE???w??

其中，气动力和气动力矩为：

1?V2S(Cx?Cxq?q)?T21Y??V2S(Cy?Cyr?r?Cyp?p)21Z??V2S(Cz?Czq?q)2

1L??V2Sc(Cl?Clr?r?Clp?p?Cl?a??a?Cl?r??r)21M??V2Sc(Cm?Cmq?q)21N??V2Sc(Cn?Cnr?r?Cnp?p?Cn?a??a?Cn?r??r)2X?气动力系数需要按如下方法通过插值求出：

Cx?Cx(?,?e)Cxq?Cxq(?)Cyr?Cyr(?)Czq?Czq(?)Cl?a?Cl?a(?,?)Cn?a?Cn?a(?,?)Cl?r?Cl?r(?,?)Cn?r?Cn?r(?,?)Cyp?Cyp(?)

Cy??0.02??0.021(?a/20)?0.086(?r/30)Cz?A(?)(1?(beta/57.3)^2)?0.19(?e/25)Cl?Cl(?,?)Cm?Cm(?,?e)Cn?Cn(?,?)Clr?Clr(?)Cmq?Cmq(?)Cnr?Cnr(?)Cnp?Cnp(?)Clp?Clp(?)下面以定直平飞状态[120m/s,1000m]作为基准运动状态，完成纵向力和力矩的配平，得到配平迎角、配平升降舵偏角和配平油门的大小。

3.2纵向力和力矩的配平

地面的重力加速度g=9.8066 m/s2、大气密度?=1.225kgm。M?MXi?Myj?Mzk，定直平飞状态下，

3MX?0，My?0，M?Mz。

纵向方向配平需要用到的力和力矩方程为：

?T?X?mgsin??m?u?qw?rv??0??Z?mgcos?cos??m?w?pv?qu??0 ?2?M?CmqSL?0.5Cm?u0Sc?0 (2.9)

纵向气动力和气动力矩的表达式为

Z?111?V2S(Cz?Czq?q)X??V2S(Cx?Cxq?q)?TM??V2Sc(Cm?Cmq?q) 222

又有??0、???和p?q?r?0，可将上式进一步化简为

由气动参数表查表可知，Cz=A*(1-(beta/57.3)^2)-0.19*(elevator/25.0)，在定直平飞状态中，beta为零，

?T?0.5Cx?V2S?mgsin??0?2?0.5Cz?VS?mgcos??0?C?0?m

(2.10)

Cz?A(?)?0.19?e25，于是配平方程进一步写为：

?T(V,H,?T)?0.5Cx(?,?e)?V2S?mgsin??0??e?2 ?0.5(A(?)?0.19)?VS?mgcos??025???Cm(?,?e)?0在给定飞行状态[V,H,?,S,m,g]时，上式为关于[?,?e,?T]的非线性方程组。用matlab的fsolve函数可以先解出[?,?e]，再解出[?T]。

每一种飞行状态对应的速度和高度如下表所示：

表3 飞机的飞行状态

飞行状态 速度(m/s) 高度(m) 飞行状态 速度(m/s) 高度(m) 1 2 3 4 5 6 7 8 9 10 11 12 50 60 70 80 90 100 110 120 130 140 150 160 0 0 0 0 0 0 1000 1000 1000 1000 5000 5000 13 14 15 16 17 18 19 20 21 22 23 24 170 180 190 200 210 220 230 240 250 260 270 280 5000 10000 10000 10000 10000 11000 11000 11000 12000 13000 14000 15000 配平攻角、升降舵偏角及油门的值为：

403550.90.80-5配平攻角（°）252015105050配平发动机油门开度（0-1）100150200飞机速度（米/秒）250300300.70.60.50.40.30.2配平舵偏角elevator（°）-10-15-20100150200飞机速度（米/秒）250300-25500.150100150200飞机速度（米/秒）250300

图1 对应所有飞行状态的配平结果

我们这一组（第四组）所取飞行状态的速度为120m/s，高度为1000m，其配平攻角、舵偏角和油门开度为：6.3296°、-0.5714°和0.169。

3.3小扰动线性化方程组

(1) 纵向小扰动方程

在第二部分中已经建立了分立的小扰动线性化方程和气动导数的求解公式，根据3.2节的配平结果，可以求得每一个气动导数，进一步求得纵向和横侧向的小扰动线性化方程。其中，纵向的小扰动线性化方程为：

??u?? -0.0180 0.0543 0.4402 -9.7438???u?? 5.7520 -0.0008??w?? -0.1624 -0.5900 111.0733 -1.0808??w?? 0 -0.1495??????????????T?

??q?? -0.0000 -0.0044 -0.7730 0??q?? 0 -0.1053?????e??????????? 0 0 1.0000 0?? 0 0????????其中，纵向系统矩阵的特征值和对应的特征向量为：

? -0.6858 + 0.6992i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i?? 0.0000 + 0.0000i -0.6858 - 0.6992i 0.0000 + 0.0000i 0.0000 + 0.0000i??? ? 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0046 + 0.0846i 0.0000 + 0.0000i??? 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.0046 - 0.0846i??和

? 0.0271 - 0.0323i 0.0271 + 0.0323i -0.9914 + 0.0000i -0.9914 + 0.0000i?? 0.9991 + 0.0000i 0.9991 + 0.0000i 0.1305 - 0.0093i 0.1305 + 0.0093i??? ? -0.0008 + 0.0062i -0.0008 - 0.0062i -0.0007 + 0.0001i -0.0007 - 0.0001i??? 0.0051 - 0.0039i 0.0051 + 0.0039i 0.0021 + 0.0086i 0.0021 - 0.0086i??(2) 横侧向小扰动方程

横侧向的小扰动线性化方程为：

?v??-0.1879 0.1243 -119.2771 9.7438??v??0.0207 0.0564??p??-0.1308 -2.4626 0.7854 0??p??-8.1387 2.0932?????????????a?

??r??0.2469 -0.0539 -0.3323 0??r??-0.3356 -1.0339????r?????????? 0 1.0000 0.1109 0? 0 0????????横侧向系统矩阵的特征值和对应的特征向量为：

? -0.2377 + 5.4470i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i?? 0.0000 + 0.0000i -0.2377 - 5.4470i 0.0000 + 0.0000i 0.0000 + 0.0000i??? ? 0.0000 + 0.0000i 0.0000 + 0.0000i -2.5357 + 0.0000i 0.0000 + 0.0000i??? 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0283 + 0.0000i??和

? 0.9987 + 0.0000i 0.9987 + 0.0000i -0.3654 + 0.0000i 0.1219 + 0.0000i?? -0.0140 + 0.0182i -0.0140 - 0.0182i -0.8661 + 0.0000i 0.0190 + 0.0000i?? ?? 0.0006 - 0.0454i 0.0006 + 0.0454i 0.0197 + 0.0000i 0.0806 + 0.0000i??? 0.0025 + 0.0024i 0.0025 - 0.0024i 0.3407 + 0.0000i 0.9891 + 0.0000i??

至此，求出了小扰动线性化方程的具体形式。

3.4计算时域响应

基于全量运动方程和小扰动线化方程，计算飞机各状态变量的时域响应。飞机的状态变量为：体轴系速度

[u,v,w]；角速度[p,q,r]；姿态角[?,?,?]；位移[x,y,z]。输入为舵偏角[?a,?e,?r]；发动机油门开度[?T]。

配平攻角、升降舵和油门为6.3296°、-0.5714和0.169。

(1). 状态变量初值为零，外界输入为零

对于全量方程：

初值为零：[u,v,w,p,q,r? 0 120*sin(6.3296/57.3) 0 0 0 0 ,?,?,,x,y,z=[120*cos(6.3296/57.3) ]6.3296/57.3 0 0 0 0]。

外界输入为零：[?a,?e,?r,?T]=[0 -0.5714°0 0.169]。 对于小扰动方程：

初值为零：[?u,w,q,??]?[0,0,0,0]和[v,p,r,?]?[0,0,0,0]。 外界输入为零：[?a,?e,?r,?T]?[0,0,0,0]。 (2). 纵向模态激励下，外界输入为零

纵向特征值为[?0.00?46对于全量方程：

对于长周期模态，模态激励为[u??u,v,w,p,?q,?r,?,???0 ,x,=[120*cos(6.3296/57.3) yz120*sin(6.3296/57.3) 0 0 0 0 6.3296/57.3 0 0 0 0]+[-0.9914,0.1305,-0.0007,0.0021]，外界输入为[?a,?e,?r,?T]=[0 -0.5714°0 0.169]；对于短周期模态，模态激励为[u??u,v,w,p,q,r,?,????,?,x,y,z]=[120*cos(6.3296/57.3) 0 120*sin(6.3296/57.3) 0 0 0 0 6.3296/57.3 0 0 0 0]+ [0.0271,0.9991,-0.0008,0.0051]，外界输入为[?a,?e,?r,?T]=[0 -0.5714° 0 0.169]。 对于小扰动方程：

对于长周期模态，模态激励为[?u,w,q,??]?[-0.9914,0.1305,-0.0007,0.0021]，外界输入为

0.i0和8[?0.6858?0.6992i]。对应的特征向量的实部分别为

[-0.9914,0.1305,-0.0007,0.0021]和[0.0271,0.9991,-0.0008,0.0051]。

[?a,?e,?r,?T]?[0,0,0,0]；对于短周期模态，模态激励为[?u,w,q,??]?[0.0271,0.9991,-0.0008,0.0051]，

外界输入为[?a,?e,?r,?T]?[0,0,0,0]。

(3). 侧向模态激励下，外界输入为零

侧向特征值为[?0.2377?5.4470i]、[?2.5357]和[0.0283]，现在取一个稳定模态[?2.5357]和一个不稳定模态[0.0283]。对应的特征向量分别为[-0.3654,-0.8661,0.0197,0.3407]和[0.1219,0.0190,0.0806,0.9891]。 对于全量方程：

对于长周期模态，模态激励为[u,v,w,p,q,r,?,?,?,x,y,z]=[120*cos(6.3296/57.3) 0 120*sin(6.3296/57.3) 0 0 0

0 6.3296/57.3 0 0 0 0]+ [0.1219,0.0190,0.0806,0.9891]，外界输入为[?a,?e,?r,?T]=[0 -0.5714°0 0.169]；对于短周期模态，模态激励为[u,v,w,p,q,r,?,?,?,x,y,z]=[120*cos(6.3296/57.3) 0 120*sin(6.3296/57.3) 0 0 0 0 6.3296/57.3 0 0 0 0]+ [-0.3654,-0.8661,0.0197,0.3407]，外界输入为[?a,?e,?r,?T]=[0 -0.5714° 0 0.169]。

对于小扰动方程：

对于长周期模态，模态激励为[v,p,r,?]?[0.1219,0.0190,0.0806,0.9891]，外界输入为

[?a,?e,?r,?T]?[0,0,0,0]；对于短周期模态，模态激励为[v,p,r,?]?[-0.3654,-0.8661,0.0197,0.3407]，外界

输入为[?a,?e,?r,?T]?[0,0,0,0]。

(4). 升降舵单位阶跃输入（1o≈0.01745 rad）

对于全量方程：

初值为零：[u,v,w,p,q,r,?,?,?,x,y,z]=[120*cos(6.3296/57.3) 0 120*sin(6.3296/57.3) 6.3296/57.3 0 0 0 0]。

外界输入为[?a,?e,?r,?T]=[0 -0.5714° 0 0.169]+[ 0 step 0 0]。 对于小扰动方程：

初值为零：[?u,w,q,??]?[0,0,0,0]和[v,p,r,?]?[0,0,0,0]。 外界输入为[?a,?e,?r,?T]?[0,0,0,0]?[0,step,0,0]。 (5). 副翼单位阶跃输入（1o≈0.01745 rad）

对于全量方程：

初值为零：[u,v,w,p,q,r,?,?,?,x,y,z]=[120*cos(6.3296/57.3) 0 120*sin(6.3296/57.3) 6.3296/57.3 0 0 0 0]。

外界输入为[?a,?e,?r,?T]=[0 -0.5714° 0 0.169]+[ step 0 0 0]。 对于小扰动方程：

初值为零：[?u,w,q,??]?[0,0,0,0]和[v,p,r,?]?[0,0,0,0]。 外界输入为[?a,?e,?r,?T]?[0,0,0,0]?[step,0,0,0]。

3.4.1初值为零输入为零

(1) 全量方程模型

全量方程模型在MATLAB中的实现如下图所示：

图2 零状态零输入的全量方程模型

全量方程在初值为零无输入情况下的响应结果如下图所示：

0 0 0 0 0 0 0 0

119.2710.8119.2680.6119.2660.413.2305飞机轴向速度（米/秒）飞机法向速度（米/秒）119.2640.20-0.2-0.4-0.6飞机侧向速度（米/秒）0100020003000仿真时间（秒）4000500013.23119.26213.2295119.2613.229119.258-0.8119.256-113.22850100020003000仿真时间（秒）400050000100020003000仿真时间（秒）40005000

图3 初值为零输入为零全量方程的三轴速度

10.80.6654x 10-610.80.6飞机滚转角速度（弧度/秒）飞机俯仰角速度（弧度/秒）0.40.20-0.2-0.4-0.6-0.8-10100020003000仿真时间（秒）400050003210-1-2-3-40100020003000仿真时间（秒）40005000飞机偏航角速度（弧度/秒）0.40.20-0.2-0.4-0.6-0.8-10100020003000仿真时间（秒）40005000

图4 初值为零输入为零全量方程的三轴角速度

10.80.60.11050.40.40.11050.110510.80.6飞机滚转角（弧度）飞机俯仰角（弧度）0.20-0.2-0.4-0.6-0.8-10100020003000仿真时间（秒）400050000.11050.11050.11050.1104飞机偏航角（弧度）0100020003000仿真时间（秒）400050000.20-0.2-0.4-0.60.11040.1104-0.8-10100020003000仿真时间（秒）40005000

图5 初值为零输入为零全量方程的三轴姿态角

(2) 小扰动线性化方程模型

小扰动线性化模型在MATLAB中的实现如下图所示：

图6 零状态零输入的小扰动线性化模型

小扰动线性化模型在初值为零无输入的情况下的响应如下图所示：

120.510.81200.61414.5飞机轴向速度（米/秒）飞机法向速度（米/秒）119.50.20-0.2-0.4-0.6-0.8飞机侧向速度（米/秒）05101520仿真时间（秒）25300.413.511913118.512.511805101520仿真时间（秒）2530-11205101520仿真时间（秒）2530

图7 初值为零输入为零小扰动方程的三轴姿速度

10.80.610.80.610.80.6飞机滚转角速度（弧度/秒）飞机俯仰角速度（弧度/秒）0.40.20-0.2-0.4-0.6-0.8-105101520仿真时间（秒）25300.40.20-0.2-0.4-0.6-0.8-105101520仿真时间（秒）2530飞机偏航角速度（弧度/秒）0.40.20-0.2-0.4-0.6-0.8-105101520仿真时间（秒）2530

图8 初值为零输入为零小扰动方程的三轴姿角速度

10.80.60.411.5飞机滚转角（弧度）0.20-0.2-0.4-0.6-0.8-105101520仿真时间（秒）2530飞机俯仰角（弧度）0.50-0.5-105101520仿真时间（秒）2530

图9 初值为零输入为零小扰动方程的姿态角

结果的分析：从图3-5和图7-9可以看出，零状态零输入情况下全量方程和小扰动方程的响应结果是一致的。小扰动方程的响应一直为零，这与小扰动方程的特性是一致的。全量方程的响应与小扰动方程的响应结果一致，在开始的地方有振荡调整过程，其纵向状态的最终值有微小的振荡。

3.4.2纵向模态激励下

全量方程在纵向模态激励下的响应，分为2种情况：以对应最大模特征值[?0.6858?0.6992i]的特征向量的实部为状态初始值；以对应最小模特征值[?0.0046?0.0846i]的特征向量的实部为状态初始值。 (1) 全量方程模型 全量方程模型的MATLAB实现如下图所示：

图10纵向模态激励全量方程模型

以对应最大模特征值[?0.6858?0.6992i]的特征向量的实部为状态初始值，纵向模态激励全量模型的响应如下图所示：

119.35114.20.80.614.114飞机轴向速度（米/秒）飞机法向速度（米/秒）119.30.20-0.2-0.4-0.6飞机侧向速度（米/秒）05001000仿真时间（秒）150020000.413.913.813.713.613.513.413.3119.25-0.8119.205001000仿真时间（秒）15002000-113.20510仿真时间（秒）1520

图11 初值为零大特征值纵向模态激励全量方程的三轴速度

10.80.60.5x 10-310.800.6飞机滚转角速度（弧度/秒）飞机俯仰角速度（弧度/秒）0.40.20-0.2-0.4-0.6-0.8-105001000仿真时间（秒）15002000飞机偏航角速度（弧度/秒）05001000仿真时间（秒）15002000-0.50.40.20-0.2-0.4-0.6-0.8-1-1.5-2-2.5-105001000仿真时间（秒）15002000

图12 初值为零大特征值纵向模态激励全量方程的三轴角速度

10.80.60.1150.40.40.1170.11610.80.6飞机滚转角（弧度）飞机俯仰角（弧度）0.20-0.2-0.4-0.6-0.8-105001000仿真时间（秒）150020000.1140.1130.1120.111飞机偏航角（弧度）05001000仿真时间（秒）150020000.20-0.2-0.4-0.60.110.109-0.8-105001000仿真时间（秒）15002000

图13 初值为零大特征值纵向模态激励全量方程的三轴姿态角

以对应最小模特征值[?0.0046?0.0846i]的特征向量的实部为状态初始值：

120.2120119.810.80.613.3613.3413.32飞机轴向速度（米/秒）飞机法向速度（米/秒）飞机侧向速度（米/秒）05001000仿真时间（秒）15002000119.6119.4119.2119118.8118.6118.4118.205001000仿真时间（秒）150020000.40.20-0.2-0.4-0.6-0.8-113.313.2813.2613.2413.2213.213.1805001000仿真时间（秒）15002000

图14 初值为零小特征值纵向模态激励全量方程的三轴速度

10.80.6864x 10-410.80.6飞机滚转角速度（弧度/秒）飞机俯仰角速度（弧度/秒）0.40.20-0.2-0.4-0.6-0.8-105001000仿真时间（秒）15002000飞机偏航角速度（弧度/秒）05001000仿真时间（秒）1500200020-2-4-6-8-100.40.20-0.2-0.4-0.6-0.8-105001000仿真时间（秒）15002000

图15 初值为零小特征值纵向模态激励全量方程的三轴角速度

10.80.60.1180.1160.11410.80.60.40.4飞机滚转角（弧度）飞机俯仰角（弧度）0.20-0.2-0.4-0.6-0.8-105001000仿真时间（秒）150020000.1120.110.1080.106飞机偏航角（弧度）05001000仿真时间（秒）150020000.20-0.2-0.4-0.60.1040.102-0.8-105001000仿真时间（秒）15002000

图16 初值为零小特征值纵向模态激励全量方程的三轴姿态角

(2) 小扰动线性化方程模型

小扰动线性化模型的MATLAB实现如下图所示：

图17纵向模态激励小扰动方程模型

小扰动方程的响应，大特征值[?0.6858?0.6992i]对应的模态激励下的响应：

119.310.8119.2950.614.5飞机轴向速度（米/秒）飞机法向速度（米/秒）飞机侧向速度（米/秒）0246仿真时间（秒）810119.290.40.20-0.2-0.4-0.614119.285119.2813.5119.275119.27-0.8119.2650246仿真时间（秒）810-1130246仿真时间（秒）810

图18 初值为零大特征值纵向模态激励小扰动方程的三轴速度

10.80.60.5x 10-310.800.6飞机滚转角速度（弧度/秒）飞机俯仰角速度（弧度/秒）0.40.20-0.2-0.4-0.6-0.8-10246仿真时间（秒）810飞机偏航角速度（弧度/秒）0246仿真时间（秒）810-0.50.40.20-0.2-0.4-0.6-0.8-1-1.5-2-2.5-10246仿真时间（秒）810

图19 初值为零大特征值纵向模态激励小扰动方程的三轴角速度

10.80.1160.60.40.1150.117飞机滚转角（弧度）0.20-0.2-0.4-0.6飞机俯仰角（弧度）0246仿真时间（秒）8100.1140.1130.1120.111-0.8-10.110246仿真时间（秒）810

图20 初值为零大特征值纵向模态激励小扰动方程的姿态角

小扰动方程的响应，小特征值[?0.0046?0.0846i]对应的模态激励下的响应：

120.2120119.810.813.40.613.45飞机轴向速度（米/秒）飞机法向速度（米/秒）119.4119.2119118.8118.6118.4118.205001000仿真时间（秒）150020000.20-0.2-0.4-0.6飞机侧向速度（米/秒）05001000仿真时间（秒）15002000119.60.413.3513.313.2513.213.15-0.8-113.105001000仿真时间（秒）15002000

图21 初值为零小特征值纵向模态激励小扰动方程的三轴速度

10.80.686420-2-4-6-8x 10-410.80.6飞机滚转角速度（弧度/秒）飞机俯仰角速度（弧度/秒）0.40.20-0.2-0.4-0.6-0.8-105001000仿真时间（秒）15002000飞机偏航角速度（弧度/秒）05001000仿真时间（秒）150020000.40.20-0.2-0.4-0.6-0.8-105001000仿真时间（秒）15002000

图22 初值为零小特征值纵向模态激励小扰动方程的三轴角速度

10.80.60.1180.1160.1140.4飞机滚转角（弧度）飞机俯仰角（弧度）05001000仿真时间（秒）150020000.20-0.2-0.4-0.6-0.8-10.1120.110.1080.1060.1040.10205001000仿真时间（秒）15002000

图23 初值为零小特征值纵向模态激励小扰动方程的姿态角

分析：对比图11-13和图14-16可知，全量方程模型在两种纵向激励作用下的响应差别不大。对比图18-20和图21-23可知，小扰动方程模型在长周期模态激励和短周期激励作用下的响应差别明显，短周期响应在4秒时趋于稳定，而长周期响应需要近1000秒时间才能趋于稳定。

从图11-23可以看出，小扰动线性化方程模型便于分析飞机的运动特性，能够准确分析飞机的长周期模态和短周期模态。而全量方程模型虽然能够体现出一定的短周期特性，但整体上分析不出飞机的运动特性。

3.4.3横侧向模态激励下

飞机横侧向特征值为[?0.2377?5.4470i]、[?2.5357]和[0.0283]，现在取一个稳定模态[?2.5357]和一个不稳定模态[0.0283]来求取响应。对应的特征向量分别为[-0.3654,-0.8661,0.0197,0.3407]和

[0.1219,0.0190,0.0806,0.9891]。

(1) 全量方程模型

全量方程的MATLAB实现如下图所示：

图24横侧向模态激励全量方程模型

以大特征值[?2.5357]对应的特征向量为初始状态求得响应如下图所示：

119.8119.7119.621.5113.413.313.2飞机轴向速度（米/秒）飞机法向速度（米/秒）飞机侧向速度（米/秒）05001000仿真时间（秒）150020000.50-0.5-1-1.5-213.11312.912.812.712.612.5119.5119.4119.3119.2119.1119-2.505001000仿真时间（秒）15002000-305001000仿真时间（秒）15002000

图25 初值为零大特征值横侧向模态激励全量方程的速度

0.10-0.100.021x 10-30.03飞机滚转角速度（弧度/秒）飞机俯仰角速度（弧度/秒）飞机偏航角速度（弧度/秒）05001000仿真时间（秒）15002000-0.2-0.3-0.4-0.5-0.6-0.7-0.8-0.905001000仿真时间（秒）150020000.01-10-2-0.01-3-0.02-4-0.03-5-0.0405001000仿真时间（秒）15002000

图26 初值为零大特征值横侧向模态激励全量方程的角速度

0.350.30.250.1140.1130.150.1120.2飞机滚转角（弧度）飞机俯仰角（弧度）0.20.150.10.050-0.050.1110.110.1090.108飞机偏航角（弧度）05001000仿真时间（秒）150020000.10.0500.1070.106-0.0505001000仿真时间（秒）1500200005001000仿真时间（秒）15002000

图27 初值为零大特征值横侧向模态激励全量方程的姿态角

以小特征值[0.0283]对应的特征向量为初始状态求得响应如下图所示：

1801.41.211517010飞机轴向速度（米/秒）飞机法向速度（米/秒）0.80.60.40.20-0.2150飞机侧向速度（米/秒）0510仿真时间（秒）152016050140-5130-10120-151100510仿真时间（秒）1520-200510仿真时间（秒）1520

图28 初值为零小特征值横侧向模态激励全量方程的速度

0.030.020.010-0.01-0.02-0.03-0.04-0.050.050.0900.08飞机滚转角速度（弧度/秒）飞机俯仰角速度（弧度/秒）-0.05飞机偏航角速度（弧度/秒）0510仿真时间（秒）15200.07-0.10.06-0.150.05-0.20.04-0.250.030510仿真时间（秒）1520-0.30.020510仿真时间（秒）1520

图29 初值为零小特征值横侧向模态激励全量方程的角速度

1.0510.950.20.10.60-0.10.50.7飞机滚转角（弧度）飞机俯仰角（弧度）-0.2-0.3-0.4-0.5-0.6飞机偏航角（弧度）0510仿真时间（秒）15200.90.850.80.750.70.650.40.30.20.1-0.70510仿真时间（秒）1520-0.800510仿真时间（秒）1520

图30 初值为零小特征值横侧向模态激励全量方程的姿态角

(2) 小扰动线性化方程模型

小扰动线性化模型在MATLAB中的实现如下图所示：

图31横侧向模态激励小扰动方程模型

小扰动线性化模型在横侧向模态激励下的响应，大实根的特征向量为状态初始值：

120.50.05014.5120-0.0514飞机轴向速度（米/秒）飞机法向速度（米/秒）-0.1-0.15-0.2-0.25-0.3-0.35119.5飞机侧向速度（米/秒）0123仿真时间（秒）4513.511913118.512.51180123仿真时间（秒）45-0.4120123仿真时间（秒）45

图32 初值为零大特征值横侧向模态激励小扰动方程的速度

0-0.1-0.210.80.61520x 10-3飞机滚转角速度（弧度/秒）飞机俯仰角速度（弧度/秒）-0.3-0.4-0.5-0.6-0.7-0.8-0.90.40.20-0.2-0.4-0.6-0.8飞机偏航角速度（弧度/秒）0123仿真时间（秒）4510500123仿真时间（秒）45-1-50123仿真时间（秒）45

图33 初值为零大特征值横侧向模态激励小扰动方程的角速度

0.350.310.251.5飞机滚转角（弧度）0.20.150.10.05飞机俯仰角（弧度）0123仿真时间（秒）450.50-0.50-0.05-10123仿真时间（秒）45

图34 初值为零大特征值横侧向模态激励小扰动方程的姿态角

小实根的特征向量为状态初始值：

120.52.514.5120214飞机轴向速度（米/秒）飞机法向速度（米/秒）119.51.5飞机侧向速度（米/秒）0204060仿真时间（秒）8010013.5119113118.50.512.51180204060仿真时间（秒）801000120204060仿真时间（秒）80100

图35 初值为零小特征值横侧向模态激励小扰动方程的速度

0.3510.81.40.30.61.2飞机滚转角速度（弧度/秒）飞机俯仰角速度（弧度/秒）0.25飞机偏航角速度（弧度/秒）0204060仿真时间（秒）801000.40.20-0.2-0.4-0.6-0.810.20.80.150.60.10.40.050.200204060仿真时间（秒）80100-100204060仿真时间（秒）80100

图36 初值为零小特征值横侧向模态激励小扰动方程的角速度

18161411.5飞机滚转角（弧度）飞机俯仰角（弧度）0204060仿真时间（秒）801001210864200.50-0.5-10204060仿真时间（秒）80100

图37 初值为零小特征值横侧向模态激励小扰动方程的姿态角

分析：横侧向长周期模态对应的特征值为稳定特征点，而短周期模态对应的特征值是不稳定的，这从图25-30和图32-37也能看出，无论是全量方程模型还是小扰动线性化模型，短周期模态激励下的状态响应都是不稳定的。对于长周期模态激励下的响应，小扰动线性化模型的结果能够很快的收敛；而全量方程模型的结果收敛速度较慢。

3.4.4升降舵单位阶跃输入

全量方程模型在MATLAB中的实现如下图所示。其中，配平升降舵偏角为-0.5714，配平油门开度为0.169，升降舵阶跃输入的幅值为0.05。 (1) 全量方程模型

图38 零状态升降舵阶跃输入全量方程模型

全量模型在升降舵阶跃输入下的响应结果如下图所示：

12912812710.813.613.40.6飞机轴向速度（米/秒）飞机法向速度（米/秒）12512412312212112011905001000仿真时间（秒）150020000.20-0.2-0.4-0.6飞机侧向速度（米/秒）05001000仿真时间（秒）150020001260.413.21312.812.612.4-0.8-112.205001000仿真时间（秒）15002000

图39 零状态升降舵阶跃输入全量方程的速度

10.80.63210-1-2-3-4-5x 10-310.80.6飞机滚转角速度（弧度/秒）飞机俯仰角速度（弧度/秒）0.40.20-0.2-0.4-0.6-0.8-105001000仿真时间（秒）15002000飞机偏航角速度（弧度/秒）05001000仿真时间（秒）150020000.40.20-0.2-0.4-0.6-0.8-105001000仿真时间（秒）15002000

图40 零状态升降舵阶跃输入全量方程的角速度

10.80.60.40.150.140.1310.80.60.4飞机滚转角（弧度）飞机俯仰角（弧度）0.20-0.2-0.4-0.6-0.8-105001000仿真时间（秒）15002000飞机偏航角（弧度）05001000仿真时间（秒）150020000.120.110.10.090.080.070.060.20-0.2-0.4-0.6-0.8-105001000仿真时间（秒）15002000

图41 零状态升降舵阶跃输入全量方程的姿态角

(2) 小扰动线性化方程模型

小扰动线性化方程模型在MATLAB中的实现如下图所示：

图42 零状态升降舵阶跃输入小扰动方程模型

小扰动线性化模型在升降舵单位阶跃输入下的响应结果如下图所示：

30028026010.80.651510飞机轴向速度（米/秒）飞机法向速度（米/秒）飞机侧向速度（米/秒）05001000仿真时间（秒）150024022020018016014012010005001000仿真时间（秒）15000.40.20-0.2-0.4-0.6-0.8-10-5-10-15-20-2505001000仿真时间（秒）1500

图43 零状态升降舵阶跃输入小扰动方程的速度

10.80.60.060.040.020-0.02-0.04-0.06-0.08-0.110.80.6飞机滚转角速度（弧度/秒）飞机俯仰角速度（弧度/秒）0.40.20-0.2-0.4-0.6-0.8-105001000仿真时间（秒）1500飞机偏航角速度（弧度/秒）05001000仿真时间（秒）15000.40.20-0.2-0.4-0.6-0.8-105001000仿真时间（秒）1500

图44 零状态升降舵阶跃输入小扰动方程的角速度

10.80.60.400.5飞机滚转角（弧度）0.20-0.2-0.4-0.6-0.8-105001000仿真时间（秒）1500飞机俯仰角（弧度）-0.5-105001000仿真时间（秒）1500

图45 零状态升降舵阶跃输入小扰动方程的姿态角

分析：从图39-45可以看出，飞机状态量最终趋于稳定，说明飞机纵向动态稳定；对比图39-41和图43-45可以看出全量方程模型和小扰动模型的动态调整时间都较长。飞机的纵向动态稳定和2对复根的实部都小于0相对应。

3.4.5副翼单位阶跃输入

副翼阶跃输入下全量方程在MATLAB中的实现如下图所示，副翼阶跃输入幅值为0.01。 (1) 全量方程模型

图46 零状态副翼阶跃输入全量方程模型

全量模型在副翼阶跃输入下的响应如下图所示：

1800.115170010飞机轴向速度（米/秒）飞机法向速度（米/秒）飞机侧向速度（米/秒）0102030仿真时间（秒）4050160-0.15150-0.20140-0.3-5130-0.4120-0.5-101100102030仿真时间（秒）4050-0.6-150102030仿真时间（秒）4050

图47 零状态副翼阶跃输入全量方程的速度

0-0.005-0.010.020-0.020-0.005-0.01飞机滚转角速度（弧度/秒）飞机俯仰角速度（弧度/秒）飞机偏航角速度（弧度/秒）0102030仿真时间（秒）4050-0.015-0.02-0.025-0.03-0.035-0.04-0.045-0.04-0.06-0.08-0.1-0.12-0.14-0.16-0.015-0.02-0.025-0.03-0.035-0.04-0.045-0.050102030仿真时间（秒）40500102030仿真时间（秒）4050

图48 零状态副翼阶跃输入全量方程的角速度

0-0.1-0.20-0.3-0.40.20.10-0.2飞机滚转角（弧度）飞机俯仰角（弧度）-0.4-0.5-0.6-0.7-0.8-0.9-10102030仿真时间（秒）4050-0.1-0.2-0.3-0.4-0.5-0.6飞机偏航角（弧度）0102030仿真时间（秒）4050-0.6-0.8-1-1.2-1.40102030仿真时间（秒）4050

图49 零状态副翼阶跃输入全量方程的姿态角

(2) 小扰动线性化方程模型

小扰动线性化方程在MATLAB中的实现如下图所示：

图50 零状态副翼阶跃输入小扰动方程模型

小扰动线性化模型在副翼单位阶跃输入下的响应如下图所示：

120.5514.5012014飞机轴向速度（米/秒）飞机法向速度（米/秒）119.5飞机侧向速度（米/秒）0510152025仿真时间（秒）303540-5-1013.5119-1513-20118.5-2512.51180510152025仿真时间（秒）303540-30120510152025仿真时间（秒）303540

图51 零状态副翼阶跃输入小扰动方程的速度

0-1-2-3-4-5-6-7-810.80.60-2-4飞机滚转角速度（弧度/秒）飞机俯仰角速度（弧度/秒）0.40.20-0.2-0.4-0.6-0.8飞机偏航角速度（弧度/秒）0510152025仿真时间（秒）303540-6-8-10-12-14-16-18-200510152025仿真时间（秒）3035400510152025仿真时间（秒）303540-1

图52 零状态副翼阶跃输入全量方程的角速度

01.5-501飞机滚转角（弧度）-100飞机俯仰角（弧度）0510152025仿真时间（秒）3035400.5-1500-200-0.5-250-10510152025仿真时间（秒）303540

图53 零状态副翼阶跃输入全量方程的姿态角

%---------飞机飞行状态4的配平攻角、升降舵偏角、油门开度---------% Result_alpha = resultalfaelevator(8,1); Result_elevator = resultalfaelevator(8,2); Result_Throtte = resultalfaelevator(8,3);

r2dalpha = Result_alpha; %配平攻角 alpha=r2dalpha/r2d;

delta_e = Result_elevator; %配平的舵偏角 Throtte = Result_Throtte; %配平的油门值

%---------配平时的侧滑角为零，副翼和方向舵均为零，姿态角速度为零---------% beta = 0; r2dbeta=0;

delta_a=0; delta_r=0; p=0; q=0; r=0;

%%--------------计算配平基础上的小扰动的气动导数-------------%% %---------计算阻尼导数----------%

Cx_q = interp1(given_Alpha,given_Cxq,r2dalpha,'spline'); Cy_r = interp1(given_Alpha,given_Cyr,r2dalpha,'spline'); Cy_p = interp1(given_Alpha,given_Cyp,r2dalpha,'spline'); Cz_q = interp1(given_Alpha,given_Czq,r2dalpha,'spline'); Cl_r = interp1(given_Alpha,given_Clr,r2dalpha,'spline'); Cl_p = interp1(given_Alpha,given_Clp,r2dalpha,'spline'); Cm_q = interp1(given_Alpha,given_Cmq,r2dalpha,'spline'); Cn_r = interp1(given_Alpha,given_Cnr,r2dalpha,'spline'); Cn_p = interp1(given_Alpha,given_Cnp,r2dalpha,'spline');

%---------计算对舵偏角的导数----------%

Cx_delta_e=(interp2(given_Elevator,given_Alpha,given_Cx,delta_e+0.0001,r2dalpha,'spline')-interp2(given_Elevator,given_Alpha,given_Cx,delta_e-0.0001,r2dalpha,'spline'))/0.0002;

Cm_delta_e=(interp2(given_Elevator,given_Alpha,given_Cm,delta_e+0.0001,r2dalpha,'spline')-interp2(given_Elevator,given_Alpha,given_Cm,delta_e-0.0001,r2dalpha,'spline'))/0.0002;

Cl_delta_a=interp2(given_delta_beta,given_Alpha,given_Cl_delta_a,r2dbeta,r2dalpha,'spline'); Cl_delta_r=interp2(given_delta_beta,given_Alpha,given_Cl_delta_r,r2dbeta,r2dalpha,'spline'); Cn_delta_a=interp2(given_delta_beta,given_Alpha,given_Cn_delta_a,r2dbeta,r2dalpha,'spline'); Cn_delta_r=interp2(given_delta_beta,given_Alpha,given_Cn_delta_r,r2dbeta,r2dalpha,'spline');

%---------计算对攻角和侧滑角的导数----------%

Cx0=interp2(given_Elevator,given_Alpha,given_Cx,delta_e,r2dalpha,'spline');

Cz0=interp1(given_Alpha,given_Cz,r2dalpha,'spline')*(1-beta^2)-0.19*(delta_e/25.0); Cm0=interp2(given_Elevator,given_Alpha,given_Cm,delta_e,r2dalpha,'spline');

Cx_a=(interp2(given_Elevator,given_Alpha,given_Cx,delta_e,r2dalpha+0.0001,'spline')-interp2(given_Elevator,given_Alpha,given_Cx,delta_e,r2dalpha-0.0001,'spline'))/0.0002;

Cz_a=(interp1(given_Alpha,given_Cz,r2dalpha+0.0001,'spline')-interp1(given_Alpha,given_Cz,r2dalpha-0.0001,'spline'))/0.0002;

Cm_a=(interp2(given_Elevator,given_Alpha,given_Cm,delta_e,r2dalpha+0.0001,'spline')-interp2(given_Elevator,given_Alpha,given_Cm,delta_e,r2dalpha-0.0001,'spline'))/0.0002; Cl_b = (sign(r2dbeta+0.0001)*interp2(given_beta,given_Alpha,given_Cl,abs(r2dbeta+0.0001),r2dalpha,'spline')-sign(r2dbeta-0.0001)*interp2(given_beta,given_Alpha,given_Cl,abs(r2dbeta-0.0001),r2dalpha,'spline'))/0.0002; Cn_b = (sign(r2dbeta+0.001)*interp2(given_beta,given_Alpha,given_Cn,abs(r2dbeta+0.001),r2dalpha,'spline')-(sign(r2dbeta-0.001)*interp2(given_beta,given_Alpha,given_Cn,abs(r2dbeta-0.001),r2dalpha,'spline')))/0.0002;

%---------计算轴向力对油门开度的导数----------%

Interp_throtte_000=interp2(given_height,given_mah,given_throtte_000,H/0.3048,Mah,'spline')*(14.5939*0.3048); Interp_throtte_077=interp2(given_height,given_mah,given_throtte_077,H/0.3048,Mah,'spline')*(14.5939*0.3048); Interp_throtte_100=interp2(given_height,given_mah,given_throtte_100,H/0.3048,Mah,'spline')*(14.5939*0.3048); matrix_interp=[Interp_throtte_000 Interp_throtte_077 Interp_throtte_100]; throtte_value=[0 0.77 1];

Cx_delta_T = ((interp1(throtte_value,matrix_interp,Throtte+0.0001))-(interp1(throtte_value,matrix_interp,Throtte-0.0001)))/0.0002;

%------------计算轴向气动力X的气动导数-------------------% Cw = m*g/(0.5*rou*V1^2*S); Xu = -rou*V1*S*Cw*sin(alpha); Xw = 0.5*rou*V1*S*Cx_a*r2d; Xq = 0.25*rou*V1*S*c_*Cx_q;

XdeltaE = 0.5*rou*V1^2*S*Cx_delta_e; XdeltaT = Cx_delta_T;

%------------计算法向气动力Z的气动导数-------------------% Zu = -rou*S*V1*Cw*cos(alpha); Zw = 0.5*rou*V1*S*Cz_a*r2d; Zq = 0.25*rou*V1*S*c_*Cz_q; ZdeltaT = 0;

ZdeltaE = 0.5*rou*V1^2*S*(-0.19/25);

%------------计算俯仰力矩M的气动导数-------------------% Mu = rou*V1*S*c_*Cm0;

Mw = 0.5*rou*V1*S*c_*Cm_a*r2d; Mq = 0.25*rou*V1*S*c_^2*Cm_q; MdeltaT = 0;

MdeltaE = 0.5*rou*V1^2*S*c_*Cm_delta_e;

%------------计算侧向力Y的气动导数-------------------% Yv = 0.5*rou*V1*S*(-0.02)*r2d; Yp = 0.25*rou*V1*S*b*Cy_p; Yr = 0.25*rou*V1*S*b*Cy_r;

YdeltaA = 0.5*rou*V1^2*S*0.021/20; YdeltaR = 0.5*rou*V1^2*S*0.086/30;

%------------计算滚转力矩L的气动导数-------------------% Lv = 0.5*rou*V1*S*b*Cl_b*r2d; Lp = 0.25*rou*V1*S*b^2*Cl_p; Lr = 0.25*rou*V1*S*b^2*Cl_r;

LdeltaA = 0.5*rou*V1^2*S*b*Cl_delta_a; LdeltaR = 0.5*rou*V1^2*S*b*Cl_delta_r;

%------------计算偏航力矩N的气动导数-------------------% Nv = 0.5*rou*V1*S*b*Cn_b*r2d; Np = 0.25*rou*V1*S*b^2*Cn_p; Nr = 0.25*rou*V1*S*b^2*Cn_r;

NdeltaA = 0.5*rou*V1^2*S*b*Cn_delta_a; NdeltaR = 0.5*rou*V1^2*S*b*Cn_delta_r;

5.2.2 求小扰动线性化方程

%%---------计算纵向小扰动线性化方程矩阵的各元素----------%%

ZA11 = Xu/m; ZA12 = Xw/m; ZA13 = Xq/m; ZA14 = -g*cos(alpha); ZA21 = Zu/m; ZA22 = Zw/m; ZA23 = Zq/m+V1; ZA24 = -g*sin(alpha); ZA31 = Mu/Iy; ZA32 = Mw/Iy; ZA33 = Mq/Iy; ZA34 = 0; ZA41 = 0; ZA42 = 0; ZA43 = 1; ZA44 = 0;

ZB11 = XdeltaT/m; ZB12 = XdeltaE/m; ZB21 = ZdeltaT/m; ZB22 = ZdeltaE/m; ZB31 = MdeltaT/Iy; ZB32 = MdeltaE/Iy; ZB41 = 0; ZB42 = 0;

ZA =[ ZA11 ZA12 ZA13 ZA14; ZA21 ZA22 ZA23 ZA24; ZA31 ZA32 ZA33 ZA34; ZA41 ZA42 ZA43 ZA44]

ZB=[ZB11 ZB12; ZB21 ZB22; ZB31 ZB32; ZB41 ZB42]

%%---------计算横侧向小扰动线性化方程矩阵的各元素----------%%

HA11 = Yv/m; HA12 = Yp/m; HA13 = Yr/m-V1; HA21 = (Iz*Lv+Ixz*Nv)/(Ix*Iz-Ixz^2); HA22 = (Iz*Lp+Ixz*Np)/(Ix*Iz-Ixz^2); HA23 = (Iz*Lr+Ixz*Nr)/(Ix*Iz-Ixz^2); HA24 = 0; HA31 = (Ixz*Lv+Ix*Nv)/(Ix*Iz-Ixz^2); HA32 = (Ixz*Lp+Ix*Np)/(Ix*Iz-Ixz^2); HA33 = (Ixz*Lr+Ix*Nr)/(Ix*Iz-Ixz^2); HA34 = 0;

HA41 = 0; HA42 = 1; HA43 = tan(alpha); HA14 = g*cos(alpha); HA44 = 0;

HB11 = YdeltaA/m; HB12 = YdeltaR/m;

HB21 = (Iz*LdeltaA+Ixz*NdeltaA)/(Ix*Iz-Ixz^2); HB22 = (Iz*LdeltaR+Ixz*NdeltaR)/(Ix*Iz-Ixz^2); HB31 = (Ixz*LdeltaA+Ix*NdeltaA)/(Ix*Iz-Ixz^2); HB32 = (Ixz*LdeltaR+Ix*NdeltaR)/(Ix*Iz-Ixz^2); HB41 = 0; HB42 = 0;

HA=[HA11 HA12 HA13 HA14; HA21 HA22 HA23 HA24; HA31 HA32 HA33 HA34; HA41 HA42 HA43 HA44] HB=[HB11 HB12; HB21 HB22; HB31 HB32; HB41 HB42]

%%---------计算小扰动线性化方程的特征值和特征矩阵----------%% [z1,lambdaz1]=eig(ZA) [h1,lambdah1]=eig(HA)

5.3 全量方程求时域响应

全量方程示意图如下所示：

5.3.1 计算力和力矩子函数

function [X, Y, Z, L, M, N] = fcn(delta_e, delta_r, delta_a, delta_T, state)

you = state(1); v = state(2); w = state(3); p = state(4); q = state(5); r = state(6); fai = state(7); theta = state(8); psi = state(9); ze = state(12);

H = 1000 - ze;

V = sqrt(you^2 + v^2 + w^2);

%地球半径% R = 6356766;

%-------飞机的常数（国际单位制）---------%

c_= 3.4503;S = 27.8709;b = 9.144;r2d=57.29577951;

Vv = sqrt(you^2+v^2+w^2);

alpha = atan(w/you);

beta = asin(v/Vv);

r2dalpha=r2d*alpha; r2dbeta=r2d*beta;

%-------飞机的惯性数据（国际单位制）---------%

Ix = 12874.8446; Iy = 75673.6077; Iz = 85552.0953; Ixz = 1331.4130; Isum=Ix*Iz-Ixz^2;

%气动参数

given_Alpha = [-10 -5 0 5 10 15 20 25 30 35 40 45]; given_Elevator = [-24 -12 0 12 24];

given_Cz =[0.7700 0.2410 -0.1000 -0.4160 -0.7310 -1.0530 -1.3660 -1.6460 -1.9170 -2.1200 -2.2480 -2.2290]; given_Cx = [-0.099 -0.048 -0.022 -0.04 -0.083 -0.081 -0.038 -0.02 -0.038 -0.073 -0.081 -0.04 -0.021 -0.039 -0.076 -0.063 -0.021 -0.004 -0.025 -0.072 -0.025 0.016 0.032 0.006 -0.046 0.044 0.083 0.094 0.062 0.012 0.097 0.127 0.128 0.087 0.024 0.113 0.137 0.13 0.085 0.025 0.145 0.162 0.154 0.1 0.043 0.167 0.177 0.161 0.11 0.053 0.174 0.179 0.155 0.104 0.047 0.166 0.167 0.138 0.091 0.04]; given_Cm = [0.205 0.081 -0.046 -0.174 -0.259 0.168 0.077 -0.02 -0.145 -0.202 0.186 0.107 -0.009 -0.121 -0.184 0.196 0.11 -0.005 -0.127 -0.193 0.213 0.11 -0.006 -0.129 -0.199 0.251 0.141 0.01 -0.102 -0.15 0.245 0.127 0.006 -0.097 -0.16 0.238 0.119 -0.001 -0.113 -0.167 0.252 0.133 0.014 -0.087 -0.104 0.231 0.108 0 -0.084 -0.076 0.198 0.081 -0.013 -0.069 -0.041 0.192 0.093 0.032 -0.006 -0.005]; given_Cl = [0 -0.001 -0.003 -0.001 0 0.07 0.009 0 -0.004 -0.009 -0.01 -0.01 -0.01 -0.011 0 -0.008 -0.017 -0.02 -0.022 -0.023 -0.023 0 -0.012 -0.024 -0.03 -0.034 -0.034 -0.037 0 -0.016 -0.03 -0.039 -0.047 -0.049 -0.05 0 -0.019 -0.034 -0.044 -0.046 -0.046 -0.047 0 -0.02 -0.04 -0.05 -0.059 -0.068 -0.074 0 -0.02 -0.037 -0.049 -0.061 -0.071 -0.079 0 -0.015 -0.016 -0.023 -0.033 -0.06 -0.091 0 -0.008 -0.002 -0.006 -0.036 -0.058 -0.076 0 -0.013 -0.1 -0.014 -0.035 -0.062 -0.077 0 -0.015 -0.19 -0.027 -0.035 -0.059 -0.076]; given_Cn = [0 0.018 0.038 0.056 0.064 0.074 0.079 0 0.019 0.042 0.057 0.077 0.086 0.09 0 0.018 0.042 0.059 0.076 0.093 0.106 0 0.019 0.042 0.058 0.074 0.089 0.106 0 0.019 0.043 0.058 0.073 0.08 0.096 0 0.018 0.039 0.053 0.057 0.062 0.08 0 0.013 0.03 0.032 0.029 0.049 0.068 0 0.007 0.017 0.012 0.007 0.022 0.03 0 0.004 0.004 0.002 0.012 0.028 0.064 0 -0.014 -0.035 -0.046 -0.034 -0.012 0.015 0 -0.017 -0.047 -0.071 -0.065 -0.002 0.011 0 -0.033 -0.057 -0.073 -0.041 -0.013 -0.001];

given_Cl_delta_a = [-0.041 -0.041 -0.042 -0.04 -0.043 -0.044 -0.043 -0.052 -0.053 -0.053 -0.052 -0.049 -0.048 -0.049 -0.053 -0.053 -0.052 -0.051 -0.048 -0.048 -0.047 -0.056 -0.053 -0.051 -0.052 -0.049 -0.047 -0.045 -0.05 -0.05 -0.049 -0.048 -0.043 -0.042 -0.042 -0.056 -0.051 -0.049 -0.048 -0.042 -0.041 -0.037 -0.082 -0.066 -0.043 -0.042 -0.042 -0.02 -0.003 -0.059 -0.043 -0.035 -0.037 -0.036 -0.028 -0.013 -0.042 -0.038 -0.026 -0.031 -0.025 -0.013 -0.01 -0.038 -0.027 -0.016 -0.026 -0.021 -0.014 -0.003 -0.027 -0.023 -0.018 -0.017 -0.016 -0.011 -0.007

-0.017 -0.016 -0.014 -0.012 -0.011 -0.01 -0.008];

given_Cl_delta_r = [0.005 0.007 0.013 0.018 0.015 0.021 0.023 0.017 0.016 0.013 0.015 0.014 0.011 0.01 0.014 0.014 0.011 0.015 0.013 0.01 0.011 0.01 0.014 0.012 0.014 0.013 0.011 0.011 -0.005 0.013 0.011 0.014 0.012 0.01 0.011 0.009 0.009 0.009 0.014 0.011 0.009 0.01 0.019 0.012 0.008 0.014 0.011 0.008 0.008 0.005 0.005 0.005 0.015 0.01 0.01 0.01 0 0 -0.002 0.013 0.008 0.006 0.006 -0.005 0.004 0.005 0.011 0.008 0.005 0.014 -0.011 0.009 0.003 0.006 0.007 0 0.02 0.008 0.007 0.005 0.001 0.003 0.001 0 ]; given_Cn_delta_a = [0.001 0.002 -0.006 -0.011 -0.015 -0.024 -0.022 -0.027 -0.014 -0.008 -0.011 -0.015 -0.01 0.002 -0.017 -0.016 -0.006 -0.01 -0.014 -0.004 -0.003 -0.013 -0.016 -0.006 -0.009 -0.012 -0.002 -0.005 -0.012 -0.014 -0.005 -0.008 -0.011 -0.001 -0.003 -0.016 -0.019 -0.008 -0.006 -0.008 0.003 -0.001 0.001 -0.021 -0.005 0 -0.002 0.014 -0.009 0.017 0.002 0.007 0.004 0.002 0.006 -0.009 0.011 0.012 0.004 0.007 0.006 -0.001 -0.001 0.017 0.015 0.007 0.1 0.012 0.004 0.003 0.008 0.015 0.006 0.004 0.011 0.004 -0.002 0.016 0.011 0.006 0.1 0.011 0.006 0.001];

given_Cn_delta_r = [-0.018 -0.028 -0.037 -0.048 -0.043 -0.052 -0.062 -0.052 -0.051 -0.041 -0.045 -0.044 -0.034 -0.034 -0.052 -0.043 -0.038 -0.045 -0.041 -0.036 -0.027 -0.052 -0.046 -0.04 -0.045 -0.041 -0.036 -0.028 -0.054 -0.045 -0.04 -0.044 -0.04 -0.035 -0.027 -0.049 -0.049 -0.038 -0.045 -0.038 -0.028 -0.027 -0.059 -0.057 -0.037 -0.047 -0.034 -0.024 -0.023 -0.051 -0.052 -0.03 -0.048 -0.035 -0.023 -0.023 -0.03 -0.03 -0.027 -0.049 -0.035 -0.02 -0.019 -0.037 -0.033 -0.024 -0.045 -0.029 -0.016 -0.009 -0.026 -0.03 -0.019 -0.033 -0.022 -0.01 -0.025 -0.013 -0.008 -0.013 -0.016 -0.009 -0.014 -0.01];

given_throtte_000=[1060 670 880 1140 1500 1860 635 425 690 1010 1330 1700 60 25 345 755 1130 1525 -1020 -710 -300 350 910 1360 -2700 -1900 -1300 -247 600 1100 -3600 -1400 -595 -342 -200 700]; given_throtte_077=[12680 9150 6200 3950 2450 1400 12680 9150 6313 4040 2470 1400 12610 9312 6610 4290 2600 1560 12640 9839 7090 4660 2840 1660 12390 10176 7750 5320 3250 1930 11680 9848 8050 6100 3800 2310];

given_throtte_100=[20000 15000 10800 7000 4000 2500 21420 15700 11225 7323 4435 2600 22700 16860 12250 8154 5000 2835 24240 18910 13760 9285 5700 3215 26070 21075 15975 11115 6860 3950 28886 23319 18300 13484 8642 5057];

given_Cxq = [-0.267 -0.11 0.308 1.34 2.08 2.91 2.76 2.05 1.5 1.49 1.83 1.21];

given_Cyr = [0.882 0.852 0.876 0.958 0.962 0.974 0.819 0.483 0.59 1.21 -0.493 -1.04]; given_Cyp = [-0.108 -0.108 -0.188 0.11 0.258 0.226 -0.344 0.362 0.611 0.529 0.298 -2.27]; given_Czq = [-8.8 -25.8 -28.9 -31.4 -31.2 -30.7 -27.7 -28.2 -29 -29.8 -38.3 -35.3]; given_Clr = [-0.126 -0.126 0.063 0.113 0.208 0.23 0.319 0.437 0.68 0.1 0.447 -0.33]; given_Clp = [-0.36 -0.359 -0.443 -0.42 -0.383 -0.375 -0.329 -0.294 -0.23 -0.21 -0.12 -0.1]; given_Cmq = [-7.21 -0.54 -5.23 -5.26 -6.11 -6.64 -5.69 -6 -6.2 -6.4 -6.6 -6];

given_Cnr = [-0.38 -0.363 -0.378 -0.386 -0.37 -0.453 -0.55 -0.582 -0.595 -0.637 -1.02 -0.804]; given_Cnp = [0.061 0.052 0.052 -0.012 -0.013 -0.024 0.05 0.15 0.13 0.158 0.24 0.15];

given_delta_beta = [-30 -20 -10 0 10 20 30]; given_beta = [0 5 10 15 20 25 30];

%%---------发动机推力插值，（英美单位制）-------------%% given_height=[0 10000 20000 30000 40000 50000]; given_mah=[0 0.2 0.4 0.6 0.8 1]';

本文来源：https://www.bwwdw.com/article/dmea.html