五一数学建模联赛寻找黑匣子.docx
《五一数学建模联赛寻找黑匣子.docx》由会员分享，可在线阅读，更多相关《五一数学建模联赛寻找黑匣子.docx（20页珍藏版）》请在冰点文库上搜索。
五一数学建模联赛寻找黑匣子
2014第十一届五一数学建模联赛
承诺书
我们仔细阅读了五一数学建模联赛的竞赛规则．
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与本队以外的任何人（包括指导教师）研究、讨论与赛题有关的问题．
我们知道，抄袭别人的成果是违反竞赛规则的,如果引用别人的成果或其它公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出．
我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性．如有违反竞赛规则的行为，我们愿意承担由此引起的一切后果．
我们授权五一数学建模联赛赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）．
我们参赛选择的题号为（从A/B/C中选择一项填写）：
A
我们的参赛报名号为：
1122参赛组别（研究生或本科或专科）：
所属学校（请填写完整的全名）
参赛队员（打印并签名）：
1.
2.
3.
日期：
2014年5月4日
获奖证书邮寄地址：
邮政编码
2014第十一届五一数学建模联赛
编号专用页
竞赛评阅编号（由竞赛评委会评阅前进行编号）：
评阅记录
评
阅
人
评
分
备
注
裁剪线裁剪线裁剪线
竞赛评阅编号（由竞赛评委会评阅前进行编号）：
参赛队伍的参赛号码：
（请各参赛队提前填写好）：
2014年第十一届五一数学建模联赛
题目寻找黑匣子
摘要
飞机在航行中可能遇到机械故障、恶劣天气、电磁波干扰等突发事件导致飞机坠毁，快速准确地找到黑匣子是破解飞机失事原因的关键．本文针对飞机失去动力后的坠落轨迹问题求解，分别建立了黑匣子在水中和空中的力学几何模型，并在MATLAB中求解出黑匣子轨迹的微分方程的数值解，并画出了坠落轨迹图。
针对问题一，本文选取力学几何模型中的微分方程，根据题中给出的飞机具体位置，找到与方程有关的相关变量，运用空气动力学公式和牛顿第二定律公式得到与阻力系数相关的变量，然后对运行轨迹和运动加速度进行分解，对水平方向和竖直方向分别进行数值微分，求得飞机水平运动距离为5701.4287m，最后将其距离值转化为经纬度值，得到落水方位为
。
针对问题二，本文与问题一同样选取力学几何模型的微分方程，只是黑匣子在水的相关阻力系数发生了改变，在竖直方向多了一个浮力．查阅资料得到相关数据，本问题假设黑匣子落到海底以下6000m，用微分方程解得黑匣子在海中运行的水平距离为109.0649m,根据海底第一区域的长度比较，得知黑匣子落在第Ⅰ区域内。
针对问题三，本题的模型建立在模型二的基础上，只是在水平方向上多了一个水平方向的力．首先确定洋流方向与飞机航行的方向，再根据飞机处的地理位置可以求出微分方程的相关阻力系数，进而得到黑匣子的在下落深度分别为1000米、2000米、3000米时其水平飞行距离为:
349.9474m、 881.3538m、1280.868m。
利用这个模型可以减少寻找黑匣子的时间，节省物力，进而可以分析空难的原因，找到减少空难的方法。
关键词力学几何模型；动力学公式；微分方程；阻力系数；MATLAB
一、问题重述
1.1问题背景
飞机是远距离航行的交通方式之一，其主要特点是速度快，安全性高．据统计，飞机是汽车、火车、轮船等几种交通方式中事故率最低的交通方式，但是飞机一旦发生事故，乘客的生还几率非常小．黑匣子是飞机专用的电子记录设备之一，它能记录各种飞行参数，供事故分析使用．黑匣子记录的参数包括：
飞机停止工作或失事坠毁前半小时的语音对话和两小时的飞行高度、速度、航向、爬升率、下降率、加速情况、耗油量、起落架放收、格林尼治时间、飞机系统工作状况和发动机工作参数等．
假设有一架飞机在高空中飞行时突然发生事故，此时飞行高度为10000米飞行速度是800公里/小时，航向东北方向45°，飞机在地面的投影位置为南纬22.0度，东经88.0度．
1.2问题提出
1.假定飞机在发生事故时突然失去动力，考虑飞机在降落过程中受到空气气流的影响，建立数学模型，描述飞机坠落轨迹并推测黑匣子的落水点．
2.假设黑匣子落水之后，不考虑洋流流动对黑匣子沉降过程的影响，建立模型描述黑匣子在水中沉降过程轨迹．如图1所示，假设黑匣子落水点所对应的海底位置为1，落水时沿着图1中指定的虚线方向沉海，给出黑匣子沉在海底的位置，并指出在图形中的哪个区域范围．
图1黑匣子在水下沉降过程中的海底剖面图
3.考虑洋流流动对黑匣子在水中沉降的影响，建立模型描述在有洋流流动的情况下黑匣子沉降轨迹方程，并求解出黑匣子沉入水下1000m，2000m和3000m时离落水点的方位．
二、问题分析
2.1问题一的分析
对于问题一，要求考虑在气流的影响下飞机的运动轨迹和落地点．当飞机未发生事故时，飞机是处于平衡状态，发生事故后，飞机将失去动力后在气流的影响下坠落．而飞机在向下坠落过程中会受到气流的干扰因素不稳定，所以我们把一切阻碍飞机向下坠落的力看成阻力，在联系飞机所受到的重力，利用空气动力学公式，我们对阻力进行分解写出在x轴和y轴的加速度的表达式．通过查阅相关数据，最后利用MATLAB软件，对时间进行分段并以高度为上限，通过不断循环求出时间和累加求出水平路程并画出轨迹，最后利用经纬度每一度之间的相差的距离关系解决此问题，求出黑匣子的落水点．
2.2问题二的分析
对于问题二，在不考虑洋流的情况下，要求建立黑匣子在水中的运动轨迹模型并指出黑匣子的沉海区域．忽略黑匣子在接触水面瞬间能量的损失，即物体速度的大小和方向在接触水面的前一时刻和后一时刻保持不变．由于物体在水中的运动过程中，还会受到水的阻力和浮力的影响，根据空气动力学方程求出阻力系数．与问题一类似，再根据物体在水中的受力情况，分别求出水平方向和竖直方向的加速度，建立微分方程，并利用MATLAB求其数值解得出黑匣子的落水区域和沉降轨迹．
2.3问题三的分析
对于问题三，求解飞机在洋流流动的影响下，分别求解黑匣子沉入水下1000m，2000m和3000m时离落水点的方位．我们可以在问题二的基础上加上水平方向的速度，竖直方向同样受到一个竖直向下的重力加速度、向上的浮力加速度、向上的阻力加速度的和加速度，这样就可以建立一个关于竖直方向加速度和水平方向加速度的微分方程．通过MATLAB来求解建立的微分方程，在MATLAB所建立的坐标轴中分别取1000m、2000m、3000m平行于x轴的直线，通过时间这一中间变量建立水平方向距离和竖直方向距离的关系，就可以得到不同海底深度黑匣子的位置．
三、模型假设
.
[1]飞机失事后，黑匣子一直保持和飞机的运动轨迹相同.
[2]地球表面假设为一个平面.
[3]坠落过程在一个竖直平面内发生.
[4]气流和洋流都不分层，介质是均匀的.
[5]气流的大小和方向是恒定不变的，且为一个常值.
[6]重力加速度一直保持不变.
[7]不考虑风对海洋表层海水的影响.
[8]飞机在与海水接触的瞬间，水平速度和竖直方向的速度大小及方向不变.
四、符号说明
符号
意义
空气的阻力系数
海水的阻力
空气的密度
海水的密度
黑匣子的质量
黑匣子只受阻力时，水平方向产生的加速度
黑匣子只受阻力时，竖直方向产生的加速度
水平方向的分速度
竖直方向的分速度
黑匣子在空气中的速度
水平方向的合加速度
竖直方向的合加速度
重力加速度
飞机失事时的高度
黑匣子的长度
黑匣子的宽度
黑匣子的高度
黑匣子的体积
黑匣子的表面积
黑匣子接触介质的面积
黑匣子在地面上的水平路程
黑匣子飞行方位的角度
飞机失事时黑匣子所处的纬度
黑匣子的重力
飞机在空中的所受的浮力
飞机在空中的动力
空中风对飞机产生的阻力
介质的阻力系数
介质的密度
物体在运动方向接触介质的横截面积
物体在介质中运动的速度
物体的加速度
物体的合力
物体的质量
五、模型的建立与求解
5.1模型的建立与求解
对于问题，本文中对于黑匣子坠落有关的数据，没有直接给出，并且还要根据这些数据和实际情况进行综合考虑，对此本文从以下几个步骤进行回答：
在飞机失事前，飞机一直在10km的高空保持平衡状态，受力分析[1]如图一．
由图一可得在平衡时竖直方向和水平方向分别有
，
当飞机突然失事后，失去了动力，受到重力和风的影响一直向下坠落．由于飞机在高空下落过程中受到的风的影响较复杂，根据运动学中的阻力是阻碍物体运动的力，所以我们把飞机受到风的阻力的方向确定为运动的反方向[2]．
由空气动力学公式得，飞机所受阻力f为
由牛顿第二定律得，飞机水平方向和竖直方向所受阻力f为
由
（1）和
（2）可得
对失事后飞机运动过程中阻力产生的加速度a进行水平和竖直方向进行分解如图2所示．
由图可得阻力产生的水平方向的加速度ax1为
竖直方向的加速度ay1为
由于飞机在水平方向不受其他力的影响，首先令
将分解后的x轴与y轴方向的相关变量带入到微分方程，记
为随着时间变化的速度变化的值，即
，然后在对速度进行按x轴和y轴分解，飞机速度方向与水平方向的夹角为
，由简单的几何关系得
所以在水平方向的加速度ax为
而在竖直方向还受重力的影响，所以竖直方向的加速度
为
根据题中给出的飞机具体位置，在网上查找出空气密度、大气压强、飞机重量、飞机横截面积等．在网上查得资料
，
，
，
，
代入（6）可得
．
利用MATLAB[3]，设置飞机高度10000m的值为上限，设置时间为0.1的断点，通过循环求出飞机降落所花时间，然后对每一个时间段内飞机所运动的水平距离进行加总求和得出黑匣子的水平距离约为33532.2米（程序和结果见附件一）．
将水平运动距离，用公式转化成经纬度．
为黑匣子航行方向与纬度的夹角，
为经度，
为黑匣子在地面上投影的水平距离，经查阅资料可知，经纬度与距离的换算公式为：
（所求经度，所求纬度）=（原经度+
），原纬度+
）．将
，
带入公式得所求黑匣子的落水点为:
。
5.2问题二模型的建立与求解
对于问题二，要求求出黑匣子在水中不受洋流影响时，其在水中的运动轨迹和黑匣子的沉底位置．
由于黑匣子在水中的受力与问题一类似，可以得到水平方向加速度
为
竖直方向的加速度为
为
由于在水中受到浮力和水的阻力的影响，根据空气动力学公式首先对
进行求解．
由问题一可得
已知空气对物体的阻力是水对物体阻力的800倍[5]，可得
黑匣子的表面积为
与问题一类似利用MATLAB软件，设置黑匣子高度为6000m的值为上限，设置时间为0.1的断点，通过循环求出黑匣子降落所花时间，然后求出每一极短时间段内黑匣子的水平运动距离，然后对其运动路程进行加总求和，得出黑匣子轨迹如图所示（程序和结果见附件二）．
因为假设黑匣子可以坠落到6000m时，黑匣子水平移动的距离还小于5km，所以黑匣子落在第Ⅰ区域．
5.3问题三模型的建立与求解
对于问题三，求解黑匣子在洋流流动的影响下，分别求解黑匣子沉入水下1000m，2000m和3000m时离落水点的方位．本文在问题二的基础上增添了洋流流动这一变量，然后通过查阅资料知道，飞机所处位置为南印度洋，此地的洋流与黑匣子航行的方向都是东北方向，再通过查阅资料可以得到此地的平均洋流速度，借助问题二的思路，同样可以建立关于水平距离和竖直距离与时间的关系，可以得到几何力学模型的微分方程[6]．
在黑匣子刚好进入水后有一个初始速度
和洋流速度
．可以得到黑匣子在水平方向的加速度为
在竖直方向黑匣子受到一个竖直向下的一个重力加速度
；在竖直方向上受到一个浮力
，产生的加速度
；在竖直方向还有一个阻力所产生的加速度
．这样就可以得到竖直方向的加速度
通过MATLAB在离海平面不同深度设置几个平行于x轴的常数方程，然后结合黑匣子沉降过程，最后对其方程进行求解，得到黑匣子的沉降轨迹方程和在不同深度的x轴方向的距离（程序和结果见附件三）。
1000米时水平飞行距离为:
349.9474m
2000米时水平飞行距离为:
881.3538m
3000米时水平飞行距离为:
1280.868m
六、模型评价及推广
6.1模型的优点
（1）对于模型得到的结果，并能联系全文不同模型所得结果，合理地分析，反复推测，最后验证模型的可行性．
（2）模型可以做到对事故地段的逼真还原，再现当时的情景，从而可以进一步飞机坠落轨迹研究．同时，在以后的事故做出科学合理地预测，从而较快的找出黑匣子并分析事发缘由．
（3）在模型求解过程中充分利用了MATLAB等软件．较好的解决了问题，得到了较理想的结果．充分用了题目中的各种信息，并较好的结合了实际情况对模型进行检验．
6.2模型的缺点
（1）在空气中时，只对问题进行了初步分析，由于在空气中风速和风向的不确定性，很难得到较理想的结果．
（2）在飞机接触水面时，飞机的能量会损失一部分，造成速度方向和大小发生改变．
七、参考文献
[1]朱慧，3种空气阻力模型及运动规律研究[J]，常德师范学院学报（自然科学版）：
2001年01期
[2]王子明；阻力和地球自转影响下的抛物体运动学方程[J]，商丘师范学院学院学报：
2003年02期
[3]苏金明，阮沈勇.MATLAB.1实用指南.北京：
电子工业出版社，2002.337-343
[4]狄龙，张启迪，存在空气阻力情况下的平抛运动轨迹研究[J]，物理通报，2014
（1）：
56-57.
[5]熊小峰，游森勇，物块落水后运动过程的动力学分析[J],教学的实践与认识，第41卷第14期：
2011，91-98.
[6]陈文义，张伟.流体力学[M],天津：
天津大学出版社，2004
附录一
附件一：
飞机飞行轨迹图程序及结果
函数文件：
functiondx=ci611（t,v）
ax=-0.0004.*（（（v
（1）.^2+v
（2）.^2）.^（1/2））.*v
（1））;
ay=9.8-0.0004.*（（v
（1）.^2+v
（2）.^2）.^（1/2））.*v
（2）;
dx=[ax;ay];
程序文件：
clear
clc
t_end=1;%假设从1秒开始計算
sum_vy=0;
whilesum_vy<10000%只在<10000m執行
[tv]=ode45（'ci611',[0:
0.1:
t_end],[2000/90]）;
%固定累加時間每秒0.1
%[T,Y]=ode45（odefun,tspan,y0）
%odefun是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名
%tspan是区间[t0tf]或者一系列散点[t0,t1,...,tf]
%y0是初始值向量《Simulink与信号处理》
%T返回列向量的时间点
%Y返回对应T的求解列向量
vy=v（:
2）;%取其y轴速度
product_vy=vy.*0.1;%下落高度
sum_vy=sum（product_vy）;%下落总高度之和
t_end=t_end+0.1;%时间每0.1变化
end%循环，直到高度大于或等于10000
disp（['坠落时间为:
'num2str（t_end）'s']）%其坠落时间
[tv]=ode45（'ci611',[0,t_end],[2000/90]）;
vx=v（:
1）;%取其x轴速度
vy=v（:
2）;%取其y轴速度
time=t;
s1=time
（1）.*vy
（1）;
product_vx1=vx.*time;
c1=cumsum（product_vx1）;
x
（1）=time
（1）;
b
（1）=s1;
fori=1:
length（time）-1%求下落高度
x（i+1）=time（i+1）-time（i）;
product_vy=vy（i）.*x（i）;
s1=s1+product_vy;
b（i+1）=s1;
end
plot（c1,10000-b）;
title（'飞行轨迹'）;xlabel（'水平路程（m）'）;ylabel（'竖直路程（m）'）
disp（['水平飞行距离为:
'num2str（max（c1））,'m']）
disp（['落水时水平飞行速度:
'num2str（vx（length（vx）））,'m/s']）
disp（['落水时竖直飞行速度:
'num2str（max（vy（length（vy））））,'m/s']）
运行结果如下
水平飞行距离为:
33307.2505m
落水时水平飞行速度:
1.9971m/s
落水时竖直飞行速度:
156.4475m/s
附件二
黑匣子沉降轨迹程序及结果
函数文件：
functiondx=ci6111（t,v）
ax=-0.0012.*（（（v
（1）.^2+v
（2）.^2）.^（1/2））.*v
（1））;
ay=9.8-0.0012.*（（v
（1）.^2+v
（2）.^2）.^（1/2））.*v
（2）-9.8*1000*0.01/20000;
dx=[ax;ay];
程序文件：
clear
clc
t_end=0.2;%假设由0.2秒开始計算
sum_vy=0;
whilesum_vy<6000%只在<6000m執行
[tv]=ode45（'ci6111',[0:
0.1:
t_end],[1.9971156.4475]）;
%固定累加時間每秒0.1
vy=v（:
2）;%取其y轴速度
vx=v（:
1）;
product_vy=vy.*0.1;%下落高度
sum_vy=sum（product_vy）;%下落总高度之和
t_end=t_end+0.1;%时间每0.1变化
end%循环，直到高度大于或等于6000
disp（['坠落时间为:
'num2str（t_end）'s']）
[tv]=ode45（'ci6111',[0,t_end],[1.9971156.4475]）;
vx=v（:
1）;%取其x轴速度
vy=v（:
2）;%取其y轴速度
time=t;
s1=time
（1）.*vy
（1）;
product_vx1=vx.*time;
c1=cumsum（product_vx1）;%求其水平运动路程
x
（1）=time
（1）;
b
（1）=s1;
fori=1:
length（time）-1
x（i+1）=time（i+1）-time（i）;
product_vy=vy（i）.*x（i）;%下落高度
s1=s1+product_vy;
b（i+1）=s1;
end
plot（c1,-b）;title（'沉降轨迹'）;
xlabel（'水平路程（m）'）;ylabel（'竖直路程（m）'）
disp（['水平飞行距离为:
'num2str（max（c1））,'m']）
运行结果如下
水平飞行距离为:
109.0649m
附件三
黑匣子在海水中的沉降轨迹程序及结果（有洋流）
函数文件：
functiondx=ci61111（t,v）
ax=-0.0012.*（（v
（1）+0.25）.^
（2））;
ay=9.8-0.0012.*（（v
（1）.^2+v
（2）.^2）.^（1/2））.*v
（2）-9.8*1000*0.01/20000;
dx=[ax;ay];
程序文件：
clear
clc
t_end=0.2;%假设由0.2秒开始計算
sum_vy=0;
whilesum_vy<6000%只在<6000m執行
[tv]=ode45（'ci61111',[0:
0.1:
t_end],[1.9971156.4475]）;
%固定累加時間每秒0.1，
vy=v（:
2）;%取其y轴速度
vx=v（:
1）;
product_vy=vy.*0.1;%下落高度
sum_vy=sum（product_vy）;%下落总高度之和
t_end=t_end+0.1;%时间每0.1变化
end%循环，直到高度大于或等于6000
[tv]=ode45（'ci61111',[0,t_end],[1.9971156.4475]）;
vx=v（:
1）;%取其x轴速度
vy=v（:
2）;%取其y轴速度
time=t;
s1=time
（1）.*vy
（1）;
product_vx=vx.*time;
c1=cumsum（product_vx）;%求其水平运动路程
x
（1）=time
（1）;
b
（1）=s1;
fori=1:
length（time）-1
x（i+1）=time（i+1）-time（i）;
product_vy=vy（i）.*x（i）;%下落高度
s1=s1+product_vy;
b（i+1）=s1;
end
cc2=ones（1,length（c1））.*（-1000）;
plot（c1,-b,c1,cc2,'-',c1,cc2.*2,'-',c1,cc2.*3,'-'）;
title（'沉降轨迹'）;xlabel（'水平路程（m）'）;ylabel（'竖直路程（m）'）;
t_end1=0.2;
t_end2=0.2;
t_end3=0.2;
sum_vy1=0;
sum_vy2=0;
sum_vy3=0;
whilesum_vy1<1000%只在<1000m執行
[t1v1]=ode45（'ci61111',[0:
0.1:
t_end1],[1.9971156.4475]）;
%固定累加時間每秒0.1，
vy1=v1（:
2）;%取其y轴速度
product_vy1=vy1.*0.1;%下落高度
sum_vy1=sum（product_vy1）;%下落总高度之和
t_end1=t_end1+0.1;%时间每0.1变化
end
运行结果如下
[t1v1]=ode45（'ci61111',[0,t_end1],[1.9971156.4475]）;
vx1=v1（:
1）;%取其x轴速度
time1=t1;
product_vx1=vx1.*time1;
c2=cumsum（product_vx1）;%求其水平运动路程
whilesum_vy2<2000%只在<2000m執行
[t2v2]=ode45（'ci61111',[0:
0.1:
t_end2],[1.9971156.4475]）;
%固定累加
升级会员 