波音737-800建模大作业

波音737—800飞机飞行模型建立实验

学院：航空自动化

专业：导航制导与控制

1 实验目的

根据飞机所提供的QAR数据，把飞机的飞行过程分为几个阶段，通过受力分析计算得出飞机在各阶段的各个时刻的地速以及飞机当时所处的地球经纬度。这之后，再把计算出来的这些数据与QAR里面的相对应的数据进行比较，得出数据误差。使我们对飞机各阶段的机体受力分析得到验证，最后确定飞机的整个飞行过程的模型。

2 实验内容

分析所得的QAR数据，根据QAR数据对飞机的飞行过程进行分阶段处理。然后查找相关资料，对飞机在飞行各阶段过程中进行受力分析。进而用MATLAB软件编写程序，计算出飞机各个阶段的地速和地球经纬度。最后把计算出来的数据和QAR里相应的数据作比较，用MA TLAB画出比较曲线图，得出计算误差，建立起飞机的飞行过程模型。在整个实验过程中要修学的课程有：《大气数据应用分析》、《导航原理与系统》、《飞机的飞行性能》、《惯性导航原理》、《MATLAB应用与编程》等等。

3 实验步骤

3.1 QAR数据分析

QAR数据分析

数据英文数据意义和用途所用仪表备注

1 东经Present

Position

Longitude

由0°本初子午线向东、西递增到180°导航仪

2 北纬Present

Position

Latitude

赤道向北递增到90°导航仪

3 磁航向Heading

Magnetic

飞机纵轴在地平面上的投影，与磁子午线的

夹角（磁北顺时针转的夹角）。磁偏角：地

球表面任一点的磁子午圈同地理子午圈的夹

角。

磁罗盘上

有罗差修

正器，已经

抵消罗差，

所以磁罗

盘测的基

本就是磁

航向。

4 标准气压高度ALTITUDE 飞机到标准气压平面的垂直距离气压式高度表

5 左无线电高度RADIO HEIGHT

Left

飞机到地面的垂直距离

无线电高

度表

6 机场标高AIR/GROUND 机场与海平面的垂直高度

7 左主起落架Left main

gear air/end

起落架用于在地面停放及滑行时支撑飞机并

使飞机在地面上灵活运动，并吸收飞机运动

时产生的撞击载荷。主要用来判断飞机是否

起飞。

8 右主起落架Right main gear air/end

9 真空速Computed

airspeed

飞机相对于空气的运动速度，根据空速可计

算地速，从而确定已飞距离和待飞时间。

空速表0.5~1.0

10 马赫数MACH 真空速与飞机所在高度的音速之比，当飞机

的M数超过临界M数时，飞机的空气动力特

马赫数表0.5~1.0

性要发生显著变化，飞机的安全性，操纵性出现一系列变化，飞行员此时根据指示空速表不能判断飞机所受空气动力的情况必须测M数。

11 CAPT显示地速CAPT display

groundspeed

飞机相对于地面的运动速度。验证自己算的

地速用到。=

地速表

12 风速飞行管理

计算机

Wind speed

FMC

测量风速，大气数据

大气数据

计算机

13 风向飞行管理

计算机

Wind

direction

FMC

测量风向，大气数据

大气数据

计算机

14 航偏角飞行管

理计算机

DRIFT ANGLE

FMC

飞行偏离目标航线的角

15 惯性垂直速度Inertial

vertical

speed

垂直方向的速度

升降速度

表

16 垂直加速度Vertical

Acceleration

垂直方向加速度，Z轴方向加速度

17 横向（纬）加

速度

Lateral

Acceleration

Y轴方向加速度，北向

18 纵向（经）加

速度

Longitudinal

Acceleration

X轴方向加速度，东向

19 左转速表LEFT ENG N1

TACHOM

测量左侧和右侧发动机的工作速度，与进气

压力表配合可以反映发动机的功率。计算发

动机推力的时候用到.

转速计

20 右转速表RIGHT ENG N1 TACHOMETER

21 飞机总重量GROSS WEIGHT 受力分析计算重力用到

22 迎角Angle of

attack

对于固定翼飞机，机翼的前进方向（相当于

气流的方向）和翼弦（与机身轴线不同）的

夹角叫迎角，也称为攻角，它是确定机翼在

气流中姿态的基准。计算发动机功率P时用

迎角传感

器

23 俯仰角CAPT DISPLAY

PITCH ATT

飞机机身轴线或者机翼弦线和水平线的夹

角。起飞后，计算阻力在水平和垂直方向上

的分力用到。

陀螺地平

仪

24 滚转角Capt Disp

Roll Att

飞机机翼所在x轴相对于z轴所转动的角度。

陀螺侧滑

仪

根据上述数据分析我得到整个飞行过程的几个图形用以直观感受：

图一

无线电测高的极限为5500英尺，所以红色曲线在下方。

图二表明的是飞机从始发点到终点的整个经纬度和高度的实时曲线，由于飞机是从北京飞往重庆的，所以曲线从里到外，经纬度在不断减小，高度由低到高再由高到低。

图二

图三是为了验证地速=空速+风速，红色曲线代表的是空速+风速的值，蓝色曲线代表的是QAR中地速的值，由于没有考虑风向，再加上QAR数据本身不准确，所以两天曲线重合度并不高。

图三

3.2 理论分析

(1)比力

加速度计是用来测量载体相对惯性空间运动加速度的传感器，加速度计测量的并不是载体的运动加速度，而是载体相对惯性空间的绝对加速度和引力加速度之差，这叫做比力。

比力的定义:

设一质点P ，质量为m ，在惯性坐标系中的位置矢量为R ，则由牛顿第二定律，有

R ..

22m ==dt

R

d m F (1)

引外F F F += (2) 这里 mG F =引 (3)

其中。错误！未找到引用源。为作用在P 点上的引力，G 为引力加速度。由此得 R m

mG F ..

=+引 (4)

或

G F R

-=..m 外

错误！未找到引用源。(5)

定义比力f 为

f m

=外

F (6)

则

(7)

即比力是作用在单位质量上的外力，比力也称作非引力加速度。

(2)比力方程

载体相对地球运动，地球又相对惯性空间运动，因此，对地球表面的运动载体而言，由加速度计输出的比力表示了载体相对惯性系的非引力加速度；而对于在地球表面导航的载体，需要知道载体相对地球系的加速度，比力方程表示了载体、地球系、惯性系这三者之间的运动关系。设载体在地心惯性坐标系中的位置矢量为R ，则利用矢量的相对导数和绝对导数之间的关系。载体的位置矢量R 在地心惯性坐标系中的导数可表达为

R w ie e

i

dt

dR dt

dR ?+=

(8)

式中错误！未找到引用源。为载体相对地球的速度；错误！未找到引用源。为地球自转的角速度，错误！未找到引用源。为地球自转产生的牵连加速度。

用错误！未找到引用源。代表载体相对地球的运动速度，即错误！未找到引用源。，则有

R w V ie ep i

dt

dR ?+= （9)

将上式两边在惯性系中求导，得

)(R 2

2R w dt

d

dt dVep dt d ie i

i

i

?+=

（10)

考虑

ie w 为常值，则

dt

dR dt

dV dt d i

ie

i

i

w ep ?+=

2

2R

(11)

因Vep 的各分量是沿平台坐标系（理论上沿导航坐标系，一般选取地理坐标系）的，故平台以坐标系P 作为动坐标系。则

ep ip p

i

v w ep ep dt

dV dt

dV ?+=

(12)

将上面（9）（12）代入（11）

)()2(2

2R R w w v w w ep ie ie ep ep ie p

i

dt

dv dt d ??+?++=

(13)

令

.

ep p

V ep dt

dV =则有

)()2(2

2.

R R w w v w w v ie ie ep ep ie ep i

dt d ??+?++= (14)

由G f R +=..

)()2(.

R w w v w w v G f ie ie ep ep ie ep ??+?++=+ (15)

考虑到地球的重力场是地球引力和地球自转产生的离心力的矢量和，即

)(R w w G g ie ie ??-= (16)

则

g v w w v f ep ep ie ep -?++=)2(.

(17)

上式即为比力方程，它是惯性导航中的系统中的一个基本方程，比力方程说明了加速度计输出的比力中所包含的物理量，其中.

ep v 是导航所需要的，其他的量都需要实时扣除掉。

(3) 加速度信息的提取 式(17)改写为

B ep ep ie ep a f g v w w f v -=-?+-=])2[(.

（18） 式中通常称B a 为有害加速度，因为它对计算所需要的.

ep v 带来了麻烦，必须要从测得的比力f 中补偿掉，才能提取出载体的运动加速度.

ep v 。B a 包含两部分，一部分是重力加速度g ，另一部分是哥氏加速度和法向加速度。

由于加速度是在机体坐标系下得到，故需要将其转换到地理坐标系上，即为

????

??????=??

????????bz by bx n b tz ty tx a a a C a a a ，

n b C 为由机体坐标系转移到地理坐标系的转移矩阵。将在地理坐标系下的加速度代入

（18）中，即可得到所需要的加速度。经过积分即可得到各轴的速度。

那么经纬度

00??+=?dt R v t

t y

t y

00c o s λ?λ+=?dt R v t

t

x t

x

0?为初始纬度，0λ为初始位移，?为纬度。

)sin 1(1?e R R e t

x -=

)s i n 321(12?e e R R e

t

y -+= e 为椭球度，e R 为赤道半径。

几个主要的坐标系：

1.地心惯性坐标系（i 系）。地球绕太阳公转，其公转速度为29.8km/h,地心和日心距离810496.1?km,公转周期为365.2422d ，向心加速度为4

1005.6-?g ，公转角速度为

。

041.0/h ，由于地球公转角速度很小，在研究地球表面的导航时，对导航的精度影响一般可

以不考虑，可将惯性坐标系的原点取在地心，且原点随地球移动，z 轴是沿地球的自转轴，x 轴，y 轴在赤道平面内，指向太阳系外的任意恒星。

2.地球坐标系（e 系）。随地球一起转动，其原点在地心，z 轴沿地球自转的方向，x 轴在赤道平面内，与零度子午线相交，y 轴与x 轴，z 轴构成右手直角坐标系，

3.地理坐标系（g 系）。原点位于载体质心，其中一坐标轴沿当地地理垂线的方向，另外两个轴在当地水平面内分别沿当地经线和纬线的切线方向。本次采用x 轴指向东，y 轴指向北，z 轴垂直与当地水平面，沿当地垂线向上。

4.载体坐标系（b 系）。原点与载体质心重合，对于飞机等巡航载体，x 轴沿载体横轴向右，y 轴沿载体纵轴向前，z 轴沿载体竖轴向下。

5.平台坐标系（p 系）。描述平台是惯导系统中平台指向的坐标系，它与平台固连。

6.导航坐标系（n 系）。它是惯性导航系统在求解导航参数时所采用的坐标系，通常，它与导航系统所在的位置有关，对平台式惯导系统来说，理想的平台坐标系就是导航坐标系，一般选取地理坐标系，对捷联式惯导系统来说，导航参数并不在载体坐标系内求解，它必须在加速度计信号分解到某个求解导航参数较为方便的坐标系，再进行导航计算，这个方便求解导航参数的坐标系就是导航坐标系，一般也是取地理坐标系。

几种坐标变换：

1.惯性坐标系--地球坐标系

由惯性坐标系i i i Z Y X 到地球坐标系e e e Z Y X ，只需经过一次平面旋转：

[]????

?

?????-==1000cos sin 0sin cos t w t w t w t w t w C ie ie ie ie z ie e i i

2.惯性坐标系--地理坐标系

需要两次坐标变换： 第一次：

[]

????

????

?

?-=10

00cos sin 0sin cos /

////

λλλλλi

z 第二次：

????

??????---=??????

?

???????++-++=???

???+-L L L L L L L L L i

Y sin 0cos 010cos 0sin )2π

cos(0)2πsin(010)2πsin(0)2πcos()2π(/ 于是，由于地心惯性坐标系到地理坐标系的坐标变换矩阵为

[]

????

?

?????------=????

?

?????-??????????---=???

???+-=L L L L L L L L L L L C i

i z Y n i sin cos sin cos cos 0cos sin cos sin sin sin cos 10

0cos sin 0sin cos sin 0cos 010cos 0sin )2π(///////

/

///

/

λλλλλλλλλλ

λ 式中，L--地理纬度，t w l l ie +-=0/λ--黄经 3.地球坐标系--地理坐标系 第一次：

[]????

?

?????-=1000cos sin 0sin cos λλλλλe

z

第二次：

????

??????---=??????

?

???????++-++=??????+-L L L L L L L L L e

Y sin 0cos 010cos 0sin )2π

cos(0)2πsin(010)2πsin(0)2πcos()2π(/ 于是，由于地球坐标系到地理坐标系的坐标变换矩阵为：

[]

??

??

?

?????------=?????

?????-??????????---=??????+-=L L L L L L L L L L L C e

e z Y n e sin cos sin cos cos 0cos sin cos sin sin sin cos 1000cos sin 0sin cos sin 0cos 010cos 0sin )2π(/

/

λλλλλλλλλλλ

式中，t w l l ie -=-=/0λλ，为本地地理纬度经度与原点地理经度之差。 4.地理坐标系--运载体坐标系

航向角（?）——运载纵轴与北向轴之间的夹角，在水平面测量，顺时针为正； 俯仰角（θ）——运载体纵轴与水平面之间的夹角，在垂直面中测量，抬头为正； 横滚角（φ）——运载体横轴与水平面之间的夹角，在横截面测量，左边抬起为正。 第一次：

[]????

??????-=10

00cos sin 0sin cos ?

????D

第二次：

[]????

??????-=θθ

θθθcos 0sin 010sin 0cos /

b

y

第三次：

[]????

?

?????-=φφ

φφφcos sin 0sin cos 0001

//

b

x

因此，由地理坐标系到运载体坐标系的坐标变换矩阵为

[][][]????

?

?????

+-+++--=??

??

?

?????-??????????-??????????-==φθφθ?φ?φθ?φ?φθφθ?φ?φθ?φ?θφ?θ?????θθ

θθ

φφ

φφ?θφcos cos cos sin sin sin cos cos sin cos sin sin sin cos sin sin sin cos cos sin sin cos cos sin sin cos sin cos cos 1000cos sin 0sin cos cos 0sin 010sin 0cos cos sin 0sin cos 0001///D

y x b n b b C

3.3 地面滑跑阶段 3.3.1受力分析

F n 为飞机发动机的推力，f 1为空气阻力，f 2为地面摩擦力，L 为飞机升力，N 为地面支持力，

G 为飞机所受的当地重力。

3.3.2地速计算

在机体坐标系下，由牛顿第二定律可列出以下等式：

/)02

1

(=-

-

==a m f

f

F a a c z

n

c y c

x 飞机 (1)

N

L

F n G

f 1

f 2

图3.3.1滑跑阶段机体受力右视图

P

P

N F

F

M C

S C

f

s

e

n

L

w D

a

L

G N L 0

12

2

2

1

/

1****4/1****5.0**

*5.0εερρυυ+=

≈-=== (2)

N f

*2

μ=

F e

为飞机的额定推力，

N

1

为飞机发动机的低压转子转速，

P

s

为当地大气压强，ε为修

正系数，

P

为标准大气压强；ρ为当地空气密度，υ为飞机空速，

S

w

为飞机机翼总面积，

C

D

为空气阻力系数；μ为跑道地面摩擦系数；

C

L

为飞机升力系数；

M

a

为飞机的马赫

数。

上述数据近似取为

F

e

=27300磅,

C

D

=0.02,

C

L

=1.8,

S

w

=1171平方英尺，

M

a

=0.23，

P

=101325pa,在水泥地板上取μ=0.03。根据以上数据，再结合QAR 数据表，

我用MATLAB 编写的一个计算此阶段各个时刻的地速的程序：

%滑跑阶段，42秒，F=ma 方法

clc;%清空命令区域

clear all ;%清空工作区与全局变量 close all ;%关闭所有窗口 %机体坐标系（b 系）

axb=0;%X 和Y 坐标轴上无加速度 azb=0;

data=xlsread('DH_data');%将EXCEL 数据（QAR 数据）导入MATLAB mach=data(1:42,8);%马赫，真空速与飞机所在高度的音速之比

%ias=data(:,7);%空速，飞机相对于空气的运动速度，单位节knot=1海里/小时=1852m/3600s

gs=data(1:42,9);%地速，飞机相对于地面的运动速度，单位节

N1=91/60;%N1转子取平均值，单位转/秒

gs=gs*1852/3600;%地速，单位从节换算为米/秒 ias=gs;%空速，QAR 数据不对，故取近似取地速

Fe=2*27300*4.448;%27300为单台发动机的额定推力，2表示双发动机，4.448是把单位换

算为牛顿的系数

m=167500*0.454;%飞机质量，将单位从磅换算为千克

g=9.8;%重力加速度，单位米/二次方秒

miu=0.03;%飞机在水泥地上受到的摩擦力系数

Cd=0.2;%空气阻力系数

Cl=1.8;%飞机升力系数

rou=1.23;%当地空气密度，单位千克/立方米

Sw=1171*0.093;%飞机机翼总面积，平方英尺换算为平方米

P0=101325;%标准大气压强

Ps=77541;%飞机所处大气压强

%以下参数为向量

kexi=1/4*mach.^2;%修正系数

G=m*g;

Fn=Fe*N1*sqrt(Ps*(1+kexi)/P0);%发动机推力

F1=0.5*rou*ias.^2*Sw*Cd;%空气阻力

L=0.5*rou*ias.^2*Sw*Cl;%飞机升力

N=G-L;%地面支持力

F2=miu*N;%地面摩擦力

ayb=(Fn-F1-F2)/m;%飞机Y轴方向的加速度

vyb=1:42;

vyb(1)=0;%初始时刻Y轴方向速度为0

for k=2:42

vyb(k)=vyb(k-1)+ayb(k);

end

hold on;

grid on;

plot(vyb,'*');

plot(gs,'r*');

title('滑跑阶段地速比较图');%图片标题

xlabel('飞行时间t:单位s','color','b');%坐标轴注释

ylabel('飞行速度:单位m/s','color','b');

legend('计算得出的地速','QAR数据给出的地速',2);%添加图例，2表示放在图形内左上角位置

此程序执行的结果为：

图3.3.2 地速计算结果比较图

3.3.3 计算结果分析

由图3.3.2可知，滑跑阶段总共42秒，误差还是非常大的。我们分析原因在于，一方面计算公式都采用近似公式，另一方面给出的初始值并不准确，同时忽略了风向，计算的累计误差越来越大。

α

3.4 爬升阶段 3.

4.1受力分析

此阶段我不考虑风向的影响，此时飞机往上爬升，受力分析可得：

F n

L

f 1

G

图3.4.1 爬升阶段机体受力右视图

r

L

N

图3.4.2 爬升阶段机体受力俯视图

3.4.2 地速计算

如图3.4.1所示，G 为飞机所受的重力，Fn 为飞机发动机的实时推力，f 1为飞机所受的空气阻力，L 为飞机所受到的升力。此时飞机是无翻转角的往上爬升。θ为飞机的俯仰角，α为飞机的航向角。

m

a m f F a a G L G c n

c y

c

x 飞机

飞机

/)cos *(/)sin *0

z

1

(θθ-=--== (3)

P

P N F F C f C s

e

n

D

L

L 0

1

2

1

2

/

1*

*****5.0***5.0ερρυυ+=== (4)

但是在图3.4.3，r 为飞机的滚转角。当飞机有滚转角时，在机体坐标系下，飞机的受力为：

m a

m a r G L r G c

z c x 飞机

飞机/cos */sin *-== （5）

在爬升阶段，由于飞机的飞行高度是不断上升的过程，在此过程中大气的各项数据也随着飞机的海拔高度的变化而变化。在公式（3）和（4）中，空气密度ρ和P

s

可以根据海拔高度

的变化去查询大气数据表而得出。升力系数

C

L

取为2，,阻力系数

C

D

取为3.5。

我编写计算此阶段的对地速度MATLAB 程序如下：

%爬升阶段（241秒），F=ma 方法

clc;%清空命令区域

clear all ;%清空工作区与全局变量 close all ;%关闭所有窗口

data=xlsread('DH_data');%将EXCEL 数据（QAR 数据）导入MATLAB mach=data(43:283,8);%马赫，真空速与飞机所在高度的音速之比

%ias=data(:,7);%空速，飞机相对于空气的运动速度，单位节knot=1海里/小时

=1852m/3600s

gs=data(43:283,9);%地速，飞机相对于地面的运动速度，单位节

head_mag=data(43:283,4);%磁航向，飞机纵轴在地平面上的投影，与磁子午线的夹角（磁北顺时针转的夹角），单位度deg

pitch=data(43:283,23);%俯仰角，机体坐标系下用到

roll=data(43:283,24);%滚转角

head_mag=head_mag*pi/180;%角度换算

pitch=pitch*pi/180;

roll=roll*pi/180;

N1=91/60;%N1转子取平均值，单位转/秒

gs=gs*1852/3600;%地速，单位从节换算为米/秒

ias=gs;%空速，QAR数据不对，故近似取地速

Fe=2*27300*4.448;%27300为单台发动机的额定推力，2表示双发动机，4.448是把单位换算为牛顿的系数

m=167500*0.454;%飞机质量，将单位从磅换算为千克

g=9;%重力加速度，单位米/二次方秒

Cd=3.5;%空气阻力系数

Cl=2;%飞机升力系数

rou=1.23;%当地空气密度，单位千克/立方米

Sw=1171*0.093;%飞机机翼总面积，平方英尺换算为平方米

P0=101325;%标准大气压强

Ps=77541;%飞机所处大气压强

%以下参数为向量

kexi=1/4*mach.^2;%修正系数

G=m*g;

Fn=Fe*N1*sqrt(Ps*(1+kexi)/P0);%发动机推力

F1=0.5*rou*ias.^2*Sw*Cd;%空气阻力

L=0.5*rou*ias.^2*Sw*Cl;%飞机升力

%机体坐标系（b系）

axb=(G*sin(roll))/m;

ayb=(Fn-F1-G*sin(pitch))/m;%飞机Y轴方向的加速度

azb=(L-G*cos(roll)-G*cos(pitch))/m;

%地理坐标系（g系）

x1=cos(head_mag).*cos(roll);

x2=sin(head_mag).*sin(pitch).*sin(roll);

x3=cos(pitch).*sin(head_mag);

x4=cos(head_mag).*sin(roll);

x5=sin(head_mag).*sin(pitch).*cos(roll);

axg=axb.*x1-x2+ayb.*x3+azb.*x4+x5;