波音737-800建模大作业.doc
《波音737-800建模大作业.doc》由会员分享，可在线阅读，更多相关《波音737-800建模大作业.doc（27页珍藏版）》请在冰点文库上搜索。
波音737—800飞机飞行模型建立实验
学院：
航空自动化
专业：
导航制导与控制
1实验目的
根据飞机所提供的QAR数据，把飞机的飞行过程分为几个阶段，通过受力分析计算得出飞机在各阶段的各个时刻的地速以及飞机当时所处的地球经纬度。
这之后，再把计算出来的这些数据与QAR里面的相对应的数据进行比较，得出数据误差。
使我们对飞机各阶段的机体受力分析得到验证，最后确定飞机的整个飞行过程的模型。
2实验内容
分析所得的QAR数据，根据QAR数据对飞机的飞行过程进行分阶段处理。
然后查找相关资料，对飞机在飞行各阶段过程中进行受力分析。
进而用MATLAB软件编写程序，计算出飞机各个阶段的地速和地球经纬度。
最后把计算出来的数据和QAR里相应的数据作比较，用MATLAB画出比较曲线图，得出计算误差，建立起飞机的飞行过程模型。
在整个实验过程中要修学的课程有：
《大气数据应用分析》、《导航原理与系统》、《飞机的飞行性能》、《惯性导航原理》、《MATLAB应用与编程》等等。
3实验步骤
3.1QAR数据分析
QAR数据分析
数据
英文
数据意义和用途
所用仪表
备注
1
东经
PresentPositionLongitude
由0°本初子午线向东、西递增到180°
导航仪
2
北纬
PresentPositionLatitude
赤道向北递增到90°
导航仪
3
磁航向
HeadingMagnetic
飞机纵轴在地平面上的投影，与磁子午线的夹角（磁北顺时针转的夹角）。
磁偏角：
地球表面任一点的磁子午圈同地理子午圈的夹角。
磁罗盘上有罗差修正器，已经抵消罗差，所以磁罗盘测的基本就是磁航向。
4
标准气压高度
ALTITUDE
飞机到标准气压平面的垂直距离
气压式高度表
5
左无线电高度
RADIOHEIGHTLeft
飞机到地面的垂直距离
无线电高度表
6
机场标高
AIR/GROUND
机场与海平面的垂直高度
7
左主起落架
Leftmaingearair/end
起落架用于在地面停放及滑行时支撑飞机并使飞机在地面上灵活运动，并吸收飞机运动时产生的撞击载荷。
主要用来判断飞机是否起飞。
8
右主起落架
Rightmaingearair/end
9
真空速
Computed
airspeed
飞机相对于空气的运动速度，根据空速可计算地速，从而确定已飞距离和待飞时间。
空速表
0.5~1.0
10
马赫数
MACH
真空速与飞机所在高度的音速之比，当飞机的M数超过临界M数时，飞机的空气动力特性要发生显著变化，飞机的安全性，操纵性出现一系列变化，飞行员此时根据指示空速表不能判断飞机所受空气动力的情况必须测M数。
马赫数表
0.5~1.0
11
CAPT显示地速
CAPTdisplay
groundspeed
飞机相对于地面的运动速度。
验证自己算的地速用到。
=
地速表
12
风速飞行管理计算机
WindspeedFMC
测量风速，大气数据
大气数据计算机
13
风向飞行管理计算机
Winddirection
FMC
测量风向，大气数据
大气数据计算机
14
航偏角飞行管理计算机
DRIFTANGLEFMC
飞行偏离目标航线的角
15
惯性垂直速度
Inertialverticalspeed
垂直方向的速度
升降速度表
16
垂直加速度
VerticalAcceleration
垂直方向加速度，Z轴方向加速度
17
横向（纬）加速度
LateralAcceleration
Y轴方向加速度，北向
18
纵向（经）加速度
LongitudinalAcceleration
X轴方向加速度，东向
19
左转速表
LEFTENGN1TACHOM
测量左侧和右侧发动机的工作速度，与进气压力表配合可以反映发动机的功率。
计算发动机推力的时候用到.
转速计
20
右转速表
RIGHTENGN1TACHOMETER
21
飞机总重量
GROSSWEIGHT
受力分析计算重力用到
22
迎角
Angleofattack
对于固定翼飞机，机翼的前进方向（相当于气流的方向）和翼弦（与机身轴线不同）的夹角叫迎角，也称为攻角，它是确定机翼在气流中姿态的基准。
计算发动机功率P时用
迎角传感器
23
俯仰角
CAPTDISPLAYPITCHATT
飞机机身轴线或者机翼弦线和水平线的夹角。
起飞后，计算阻力在水平和垂直方向上的分力用到。
陀螺地平仪
24
滚转角
CaptDispRollAtt
飞机机翼所在x轴相对于z轴所转动的角度。
陀螺侧滑仪
根据上述数据分析我得到整个飞行过程的几个图形用以直观感受：
图一
无线电测高的极限为5500英尺，所以红色曲线在下方。
图二表明的是飞机从始发点到终点的整个经纬度和高度的实时曲线，由于飞机是从北京飞往重庆的，所以曲线从里到外，经纬度在不断减小，高度由低到高再由高到低。
图二
图三是为了验证地速=空速+风速，红色曲线代表的是空速+风速的值，蓝色曲线代表的是QAR中地速的值，由于没有考虑风向，再加上QAR数据本身不准确，所以两天曲线重合度并不高。
图三
3.2理论分析
（1）比力
加速度计是用来测量载体相对惯性空间运动加速度的传感器，加速度计测量的并不是载体的运动加速度，而是载体相对惯性空间的绝对加速度和引力加速度之差，这叫做比力。
比力的定义:
设一质点P，质量为m，在惯性坐标系中的位置矢量为R，则由牛顿第二定律，有
（1）
（2）
这里（3）
其中。
为作用在P点上的引力，G为引力加速度。
由此得
（4）
或（5）
定义比力f为
（6）
则
（7）
即比力是作用在单位质量上的外力，比力也称作非引力加速度。
（2）比力方程
载体相对地球运动，地球又相对惯性空间运动，因此，对地球表面的运动载体而言，由加速度计输出的比力表示了载体相对惯性系的非引力加速度；而对于在地球表面导航的载体，需要知道载体相对地球系的加速度，比力方程表示了载体、地球系、惯性系这三者之间的运动关系。
设载体在地心惯性坐标系中的位置矢量为R，则利用矢量的相对导数和绝对导数之间的关系。
载体的位置矢量R在地心惯性坐标系中的导数可表达为
（8）
式中为载体相对地球的速度；为地球自转的角速度，为地球自转产生的牵连加速度。
用代表载体相对地球的运动速度，即，则有
（9）
将上式两边在惯性系中求导，得
（10）
考虑为常值，则
（11）
因Vep的各分量是沿平台坐标系（理论上沿导航坐标系，一般选取地理坐标系）的，故平台以坐标系P作为动坐标系。
则
（12）
将上面（9）（12）代入（11）
（13）
令则有
（14）
由
得
（15）
考虑到地球的重力场是地球引力和地球自转产生的离心力的矢量和，即
（16）
则
（17）
上式即为比力方程，它是惯性导航中的系统中的一个基本方程，比力方程说明了加速度计输出的比力中所包含的物理量，其中是导航所需要的，其他的量都需要实时扣除掉。
（3）加速度信息的提取
式（17）改写为
（18）
式中通常称为有害加速度，因为它对计算所需要的带来了麻烦，必须要从测得的比力f中补偿掉，才能提取出载体的运动加速度。
包含两部分，一部分是重力加速度g，另一部分是哥氏加速度和法向加速度。
由于加速度是在机体坐标系下得到，故需要将其转换到地理坐标系上，即为
，
为由机体坐标系转移到地理坐标系的转移矩阵。
将在地理坐标系下的加速度代入（18）中，即可得到所需要的加速度。
经过积分即可得到各轴的速度。
那么经纬度
为初始纬度，为初始位移，为纬度。
e为椭球度，为赤道半径。
几个主要的坐标系：
1.地心惯性坐标系（i系）。
地球绕太阳公转，其公转速度为29.8km/h,地心和日心距离km,公转周期为365.2422d，向心加速度为g，公转角速度为/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.惯性坐标系--地球坐标系
由惯性坐标系到地球坐标系，只需经过一次平面旋转：
2.惯性坐标系--地理坐标系
需要两次坐标变换：
第一次：
第二次：
于是，由于地心惯性坐标系到地理坐标系的坐标变换矩阵为
式中，L--地理纬度，--黄经
3.地球坐标系--地理坐标系
第一次：
第二次：
于是，由于地球坐标系到地理坐标系的坐标变换矩阵为：
式中，，为本地地理纬度经度与原点地理经度之差。
4.地理坐标系--运载体坐标系
航向角（）——运载纵轴与北向轴之间的夹角，在水平面测量，顺时针为正；
俯仰角（）——运载体纵轴与水平面之间的夹角，在垂直面中测量，抬头为正；
横滚角（）——运载体横轴与水平面之间的夹角，在横截面测量，左边抬起为正。
第一次：
第二次：
第三次：
因此，由地理坐标系到运载体坐标系的坐标变换矩阵为
3.3地面滑跑阶段
N
L
Fn
G
f1
f2
图3.3.1滑跑阶段机体受力右视图
3.3.1受力分析
Fn为飞机发动机的推力，f1为空气阻力，f2为地面摩擦力，L为飞机升力，N为地面支持力，G为飞机所受的当地重力。
3.3.2地速计算
在机体坐标系下，由牛顿第二定律可列出以下等式：
（1）
（2）
为飞机的额定推力，为飞机发动机的低压转子转速，为当地大气压强，为修正系数，为标准大气压强；为当地空气密度，为飞机空速，为飞机机翼总面积，为空气阻力系数；为跑道地面摩擦系数；为飞机升力系数；为飞机的马赫数。
上述数据近似取为=27300磅,=0.02,=1.8,=1171平方英尺，=0.23，=101325pa,在水泥地板上取=0.03。
根据以上数据，再结合QAR数据表，我用MATLAB编写的一个计算此阶段各个时刻的地速的程序：
%滑跑阶段，42秒，F=ma方法
clc;%清空命令区域
clearall;%清空工作区与全局变量
closeall;%关闭所有窗口
%机体坐标系（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
fork=2:
42
vyb（k）=vyb（k-1）+ayb（k）;
end
holdon;
gridon;
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受力分析
此阶段我不考虑风向的影响，此时飞机往上爬升，受力分析可得：
Fn
L
f1
G
图3.4.1爬升阶段机体受力右视图
N
图3.4.2爬升阶段机体受力俯视图
α
图3.4.3爬升阶段机体受力后视图
G
r
L
3.4.2地速计算
如图3.4.1所示，G为飞机所受的重力，Fn为飞机发动机的实时推力，f1为飞机所受的空气阻力，L为飞机所受到的升力。
此时飞机是无翻转角的往上爬升。
为飞机的俯仰角，α为飞机的航向角。
（3）
（4）
但是在图3.4.3，r为飞机的滚转角。
当飞机有滚转角时，在机体坐标系下，飞机的受力为：
（5）
在爬升阶段，由于飞机的飞行高度是不断上升的过程，在此过程中大气的各项数据也随着飞机的海拔高度的变化而变化。
在公式（3）和（4）中，空气密度和可以根据海拔高度的变化去查询大气数据表而得出。
升力系数取为2，,阻力系数取为3.5。
我编写计算此阶段的对地速度MATLAB程序如下：
%爬升阶段（241秒），F=ma方法
clc;%清空命令区域
clearall;%清空工作区与全局变量
closeall;%关闭所有窗口
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;
y1=sin（head_mag）.*cos（roll）;
y2=cos（head_mag）.*sin（pitch）.*sin（roll）;
y3=cos（pitch）.*cos（head_mag）;
y4=sin（head_mag）.*sin（roll）;
y5=cos（head_mag）.*sin（pitch）.*cos（roll）;
ayg=axb.*y1+y2+ayb.*y3+azb.*y4-y5;
z1=-sin（roll）.*cos（head_mag）;
z2=sin（head_mag）;
z3=cos（head_mag）.*cos（roll）;
azg=axb.*z1+ayb.*z2+azb.*z3;
vxg=1:
241;
vyg=1:
241;
vg=1:
241;
vxg
（1）=0;
vyg
（1）=gs
（1）;
vg
（1）=gs
（1）;
fork=2:
241
vxg（k）=vxg（k-1）+axg（k）;
vyg（k）=vyg（k-1）+ayg（k）;
vg（k）=sqrt（vxg（k）^2+vyg（k）^2）;
end
t=1:
241;
hold;
gridon;
boxon;
p
升级会员 