波音737-800建模大作业.doc

波音737800飞机飞行模型建立实验 学院：航空自动化 专业：导航制导与控制 1 实验目的 根据飞机所提供的QAR数据，把飞机的飞行过程分为几个阶段，通过受力分析计算得出飞机在各阶段的各个时刻的地速以及飞机当时所处的地球经纬度。这之后，再把计算出来的这些数据与QAR里面的相对应的数据进行比较，得出数据误差。使我们对飞机各阶段的机体受力分析得到验证，最后确定飞机的整个飞行过程的模型。2 实验内容 分析所得的QAR数据，根据QAR数据对飞机的飞行过程进行分阶段处理。然后查找相关资料，对飞机在飞行各阶段过程中进行受力分析。进而用MATLAB软件编写程序，计算出飞机各个阶段的地速和地球经纬度。最后把计算出来的数据和QAR里相应的数据作比较，用MATLAB画出比较曲线图，得出计算误差，建立起飞机的飞行过程模型。在整个实验过程中要修学的课程有：大气数据应用分析 、导航原理与系统 、飞机的飞行性能 、惯性导航原理 、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/end9真空速Computed airspeed飞机相对于空气的运动速度，根据空速可计算地速，从而确定已飞距离和待飞时间。空速表0.51.010马赫数MACH真空速与飞机所在高度的音速之比，当飞机的M数超过临界M数时，飞机的空气动力特性要发生显著变化，飞机的安全性，操纵性出现一系列变化，飞行员此时根据指示空速表不能判断飞机所受空气动力的情况必须测M数。马赫数表0.51.011CAPT显示地速CAPT display groundspeed飞机相对于地面的运动速度。验证自己算的地速用到。=地速表12风速飞行管理计算机Wind speed FMC测量风速，大气数据大气数据计算机13风向飞行管理计算机Wind directionFMC测量风向，大气数据大气数据计算机14航偏角飞行管理计算机DRIFT ANGLE FMC飞行偏离目标航线的角15惯性垂直速度Inertial vertical speed垂直方向的速度升降速度表16垂直加速度Vertical Acceleration垂直方向加速度，Z轴方向加速度17横向（纬）加速度Lateral AccelerationY轴方向加速度，北向18纵向（经）加速度Longitudinal AccelerationX轴方向加速度，东向19左转速表LEFT ENG N1 TACHOM测量左侧和右侧发动机的工作速度，与进气压力表配合可以反映发动机的功率。计算发动机推力的时候用到.转速计20右转速表RIGHT ENG N1 TACHOMETER21飞机总重量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，则由牛顿第二定律，有 (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 地面滑跑阶段NLFnGf1f2图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;%清空命令区域clear all;%清空工作区与全局变量close all;%关闭所有窗口%机体坐标系（b系）axb=0;%X和Y坐标轴上无加速度azb=0; data=xlsread(DH_data);%将EXCEL数据（QAR数据）导入MATLABmach=data(1:42,8);%马赫，真空速与飞机所在高度的音速之比%ias=data(:,7);%空速，飞机相对于空气的运动速度，单位节knot=1海里/小时=1852m/3600sgs=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轴方向速度为0for 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受力分析 此阶段我不考虑风向的影响，此时飞机往上爬升，受力分析可得：FnLf1G图3.4.1 爬升阶段机体受力右视图N图3.4.2 爬升阶段机体受力俯视图图3.4.3 爬升阶段机体受力后视图GrL3.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;%清空命令区域clear all;%清空工作区与全局变量close all;%关闭所有窗口 data=xlsread(DH_data);%将EXCEL数据（QAR数据）导入MATLABmach=data(43:283,8);%马赫，真空速与飞机所在高度的音速之比%ias=data(:,7);%空速，飞机相对于空气的运动速度，单位节knot=1海里/小时=1852m/3600sgs=data(43:283,9);%地速，飞机相对于地面的运动速度，单位节head_mag=data(43:283,4);%磁航向，飞机纵轴在地平面上的投影，与磁子午线的夹角（磁北顺时针转的夹角），单位度degpitch=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);for k=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; grid on;box on;plot(t,vg,t,gs,r,LineWidth,2);hold;title(地速比较图);%图片标题xlabel(飞行时间t:单位s,color,b);%坐标轴注释ylabel(飞行速度:单位m/s,color,b);legend(计算得出的地速,QAR数据给出的地速,2);%添加图例，1表示放在图形内右上角位置 pause();plotyy(t,vg,t,gs); grid on;box on;title(地速比较图);%图片标题xlabel(飞行时间t:单位s,color,b);%坐标轴注释ylabel(飞行速度:单位m/s,color,b);legend(计算得出的地速,QAR数据给出的地速,2);%添加图例，1表示放在图形内右上角位置 程序运行所得出的对地速度结果如下：图3.4.4 地速计算结果比较图图3.4.5 地速计算结果比较图3.4.3 计算结果分析 由图3.4.4和图3.4.5可知，计算的结果和QAR里面的数据还是有很大的误差的。为何出现如此大的误差，首先我们通过受力分析得出的力有很多是随着时间不短变化的，它们不是一个定值。其次，我们还没有把风力对飞机的受力考虑进去。而且误差会随着时间的增加而积累。3.5巡航阶段巡航阶段可以仿照爬升阶段来做，受力分析和数据都是一样的，巡航阶段飞机也并不是一直处在一个高度不变的，这个从图一就能清晰地看出来。我编写的计算地速的MATLAB程序如下：%巡航阶段，F=ma方法 clc;%清空命令区域clear all;%清空工作区与全局变量close all;%关闭所有窗口 data=xlsread(DH_data);%将EXCEL数据（QAR数据）导入MATLABmm,nn=size(data);%计算数据的行mm数和列数nnmach=data(283:mm,8);%马赫，真空速与飞机所在高度的音速之比%ias=data(:,7);%空速，飞机相对于空气的运动速度，单位节knot=1海里/小时=1852m/3600sgs=data(283:mm,9);%地速，飞机相对于地面的运动速度，单位节head_mag=data(283:mm,4);%磁航向，飞机纵轴在地平面上的投影，与磁子午线的夹角（磁北顺时针转的夹角），单位度degpitch=data(283:mm,23);%俯仰角，机体坐标系下用到roll=data(283:mm,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(r