铅垂平面内无控飞行导弹的
弹道仿真及分析报告
姓名：韩宏伟 学号：1120100202
班级：01811001 学院：宇航学院
日期：2013年6月
铅垂平面内无控飞行导弹弹道的仿真及分析报告
（一）摘要
本次实验所做的是无控飞行导弹的弹道仿真和分析。本次报告首先从无控导弹的数
学模型进行分析，确定所需参数和求解所需的积分方法，之后对在迭代过程需要用到的参数数据进行插值处理。整个插值过程采用MATLAB进行编程插值，本次实验中采取边迭代边插值，所以迭代程序仍由MATLAB编制。最后把得出的结论用曲线的形式表示出来，以直观分析无控飞行弹道。最后分别分析每一个参数的变化规律，从而得出导弹飞行过程中所有参数变化的规律。
关键词：弹道 插值 龙格库塔法 MATLAB
Abstract
The aim of this experience is to make simulation and analysis for ballistic of non-controlled missile. And firstly, test report gives out analysis about mathematical model of missile,then, needed parameters and integral method have been chosen, moreover, parameters have been interpolated through using MATLAB. When using MATLAB program to iterate, computer also goes on last step. After that, some variables’ curves will be displayed, which can be used to analyze specific ballistic of the missile. At last, every parameter’s change law has been given out, and then whole laws of missile is apparent.
Keywords: ballistic, interpolate, Rung-Kutta method, MATLAB
（二）实验内容及要求
根据描述飞行器在铅垂平面内运动的数学模型，编制某导弹的铅垂平面无控飞行弹道仿真程序，利用计算机解算初始段无控飞行弹道，对初始段弹道参数的变化规律进行分析。
（三）实验具体仿真及分析过程
1.数学模型 1.1 原理分析
导弹在空间中的飞行可以用严格的数学模型表示出来，本次实验以铅垂面内飞行导弹的弹道数学模型为基础仿真并分析其具体飞行过程。
铅垂面内无控导弹数学模型如下：
（1）
其中（空气动力及空气动力矩）：
1?V2S21 Y?cy?V2S2X?cx（2）
1?z??2zMz?Mz??M???(m??m?)?VSLzzzzz2
以上数学模型基于导弹仅在铅垂面的运动，也即只考虑导弹在铅垂面内运动时的参数变化。
1.2 模型分析方法
在分析仿真过程中主要需要求解微分方程组，直接得到解析解的情况和可能性极小，所以采用微分方程的数值解法。龙格-库塔法就是一种常用的数值求解工具，由于其拥有
精度高和易于编程的优点，现在微分方程多采用此法，本次试验也选用龙格-库塔四阶法作为仿真迭代过程的积分方法，且步长选为h=0.005。
2.实验工具
本次实验选用Matlab作为仿真软件。
3.实验初值及原始数据 3.1 准备原始数据
m 求解导弹运动方程组，必须给出所需的原始数据，它们一般来源于总体初步设计、
估算和实验结果。这些原始数据可能是以曲线或表格形式给出的，也可以用拟合表达式给定。要求解方程组（1），需要如下数据：
（1）标准大气参数，包括大气密度ρ，声速C和重力加速度g。 （2）导弹气动力及气动力矩相关数据。
（3）推力P，燃料秒流量mc，质心位置xG和转动惯量Jz。 （4）导弹外形尺寸，特征面积及特征长度。 （5）积分初值。
3.2 原始数据
3.2.1 初始数据
x=0(m) y=20.0(m) ?=18? ?=18? v=20(m/s) m=52.38(kg) ?z=0(rad/s)
其中：x，y是导弹位置坐标；?是导弹俯仰角；?是导弹弹道倾角；v是导弹飞行速度；m是导弹质量；?z 是导弹俯仰角速度。
3.2.2 攻角和马赫数范围 攻角ɑ 0～10o; 马赫数Ma 0～0.9 注：给出的范围仅用于插值计算。 3.2.3 基本常数表
特征面积(m2) 0.0227
特征长度(m) 1.8 毛翼展(m) 0.5 音速(m/s) 343.13 大气密度(kg/m3) 1.225 重力加速度 （m/s2） 9.8 3.2.4 阻力系数表
马赫数 攻角ɑ(?) Ma 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
0 .4177 .3858 .3779 .3785 .3787 .3829 .3855 .4082 .4947 2 .4404 .4086 .4007 .4015 .4018 .4062 .4091 .4321 .5192 4 .5219 .4903 .4827 .4838 .4846 .4897 .4934 .5175 .6073 6 .6603 .6290 .6218 .6234 .6249 .6310 .6363 .6621 .7571 8 .8534 .8226 .8160 .8184 .8209 .8284 .8358 .8641 .9672 10 1.1023 1.0723 1.0666 1.0700 1.0738 1.0835 1.0938 1.1254 1.2392 3.2.5 升力系数表
马赫数 Ma 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 2 攻角ɑ(?) 4 1.4758 1.4807 1.4858 1.4923 1.5007 1.5134 1.5304 6 2.2870 2.2942 2.3014 2.3107 2.3227 2.3409 2.3661 8 3.0713 3.0814 3.0915 3.1039 3.1197 3.1436 3.1775 10 3.8463 3.8598 3.8731 3.8891 3.9092 3.9401 3.9835 .0000 .6430 .0000 .6454 .0000 .6480 .0000 .6512 .0000 .6554 .0000 .6617 .0000 .6698 0.8 0.9
.0000 .6792 .0000 .6933 1.5501 1.5935 2.3950 2.4706 3.2162 3.3273 4.0323 4.1790 3.2.6 推力数据
t(s) .000 .15 .49 2.11 2.27 3.53 8.78 25.45 42.80 43.68 44.08 P(kgf) 331.2 614.3 505.4 607.8 48.65 43.97 42.01 41.00 40.80 40.79 2.22 注： 1kgf=1gN (其中kgf是单位质量力，g为当地重力加速度)。
第一级工作结束时间：2.1126s， 第二级工作结束时间：44.0832s。
3.2.7 发动机质量秒流量
t(s) 秒流量(kg/s) 3.2.8 导弹转动惯量
0. 2.1 2.105 44.1 44.105 100 0. 2.362 2.362 0.21059 0.21059 0. t(s) Jz(kg*m2)
.0 2.0 2.4 6.4 10.4 14.4 18.4 22.4 26.4 30.4 34.0 38.4 42.4 44.0 8.35 7.88 7.86 7.81 7.78 7.75 7.73 7.71 7.70 7.70 7.69 7.69 7.69 7.69 3.2.9 导弹重心（起自头部）
t(s) .0 2.0 2.4 10.0 18.0 26.0 32.0 38.0 42.0 44.0 XG(m) .9381 .9095 .9091 .9026 .8969 .8928 .8907 .8896 .8895 .8896
?m3.2.10 静稳定力矩系数z0?Xg?Xg0
马赫数 Ma 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 2 4 攻角ɑ(?) 6 8 10 0.0000 -0.0104 -0.0341 -0.0564 -0.0771 -0.0985 0.0000 -0.0104 -0.0341 -0.0564 -0.0770 -0.0983 0.0000 -0.0104 -0.0341 -0.0564 -0.0769 -0.0982 0.0000 -0.0105 -0.0342 -0.0564 -0.0768 -0.0979 0.0000 -0.0104 -0.0339 -0.0560 -0.0761 -0.0969 0.0000 -0.0093 -0.0314 -0.0521 -0.0708 -0.0903 0.0000 -0.0080 -0.0286 -0.0477 -0.0650 -0.0829 0.0000 -0.0065 -0.0252 -0.0425 -0.0578 -0.0739 0.0000 -0.0053 -0.0229 -0.0391 -0.0538 -0.0693 注：当导弹重心变化时的修正公式： ααmzα?m1z0α?cy(Xg?Xg0)/L
3.2.11 阻尼力矩系数导数
当Xg=.9381时 马赫数 Ma 0.1 0.2 0.3 0 2 攻角ɑ(?) 4 6 8 10 -0.4686 -0.4829 -0.4982 -0.5130 -0.5272 -0.5409 -0.4707 -0.4850 -0.5003 -0.5150 -0.5292 -0.5429 -0.4744 -0.4886 -0.5039 -0.5186 -0.5327 -0.5464 0.4 0.5 0.6 0.7 0.8 0.9 -0.4797 -0.4939 -0.5090 -0.5237 -0.5378 -0.5514 -0.4882 -0.5022 -0.5173 -0.5318 -0.5458 -0.5593 -0.5089 -0.5227 -0.5376 -0.5520 -0.5658 -0.5791 -0.5366 -0.5502 -0.5649 -0.5790 -0.5927 -0.6058 -0.5738 -0.5871 -0.6014 -0.6153 -0.6287 -0.6415 -0.6272 -0.6407 -0.6553 -0.6694 -0.6830 -0.6960
当Xg=.8896时 马赫数 Ma 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
3.2.12 插值算法
攻角ɑ(?) 0 2 4 6 8 10 -0.6179 -0.6384 -0.6600 -0.6805 -0.6999 -0.7182 -0.6207 -0.6410 -0.6626 -0.6830 -0.7024 -0.7207 -0.6253 -0.6455 -0.6670 -0.6874 -0.7067 -0.7249 -0.6319 -0.6521 -0.6734 -0.6937 -0.7129 -0.7310 -0.6424 -0.6624 -0.6835 -0.7036 -0.7226 -0.7406 -0.6669 -0.6866 -0.7074 -0.7272 -0.7459 -0.7636 -0.6997 -0.7190 -0.7395 -0.7589 -0.7774 -0.7948 -0.7435 -0.7624 -0.7824 -0.8014 -0.8194 -0.8365 -0.8069 -0.8266 -0.8474 -0.8672 -0.8859 -0.9035 气动数据插值——等距双变元样条插值；
推力、重心、转动惯量等——不等距一元线性插值。
4.弹道计算的积分算法（即龙格库塔法）
龙格库塔四阶法的本质是利用微分方程右端子函数在若干点上的函数值的线性组
合来构造近似公式，构造时要求近似公式在任一点处的Taylor展开式与解处的Taylor展开式的前几项重合，从而使近似公式达到所需的阶数，和更高的精度。
四阶龙格库塔法的计算公式如下：
?dy?f?x,y??对于形如?dx??y(x0)?y0计算公式为：
x??a,b?y?Rn 的微分方程组
K1?f?xi,yi?
hh??
K2?f?xi?,yi?K1?22??hh??K3?f?xi?,yi?K2?22 ??
K4?f?xi?h,yi?hK3?yi?1?yi?h?K1?2K2?2K3?K4? 65.数据仿真
5.1程序编制
程序编制的基本思路如下： main函数 实参 Rungkuta子函数 插值数据 实参 主函数 结果 Shujuchuli子函数，主要对已有数据进行插值，并返回所需点的插值结果数据。
注：（1）程序中所用变量如X，Y和方程式中导弹位置对应； （2）程序中Wz为俯仰角速度；
（3）程序中Big_theta为导弹俯仰角；
（4）程序中Theta指弹道倾角；
（5）程序中alpha1指每一个迭代时刻的导弹攻角； （6）程序中M指导弹质量。
（7）程序中其它参数如cx,cy,L,S等完全与方程（1），（2）一一对应。 （8）剩余为程序中的中间变量，可在程序中轻松辨识。 5.2 仿真结果及图线
运行程序得到如下导弹飞行过程中的各变化参数的变化曲线：
导弹飞行弹道曲线150100Y/m500050010001500X/m200025003000
导弹飞行速度曲线300250200V/(m/s)1501005000246t/s81012
导弹俯仰角变化曲线2015105?/:o0-5-10-150246t/s81012
导弹弹道倾角变化曲线2015105?/:o0-5-10-150246t/s81012
导弹攻角变化曲线3.532.521.5a/rad10.50-0.5-1-1.50246t/s81012
导弹俯仰角速度变化曲线0.10.050wz/(rad/s)-0.05-0.1-0.15-0.2-0.250246t/s81012
导弹质量变化曲线53525150m/kg49484746450246t/s81012
5.3 结果分析
5.3.1 飞行弹道变化规律
飞行弹道近似为一个抛物线，起飞阶段在发动机的推动下先爬升，爬升过程中随着时间推移（本质是弹道倾角缓慢较小），爬升率下降；当t=6s，导弹爬升到最高点，即当x=1234.68m，y达到最高点，ymax=145.5m；然后开始下落，此时弹道倾角和俯仰角都继续减小为负值，重力沿切向的分量增大，并有阻力变为动力，使导弹加速下落，在最终的射程为2679.1m.
导弹飞行水平距离曲线30002000X/m1000002468t/s导弹飞行高度曲线1012150100Y/m5000246t/s81012
5.3.2 飞行速度变化规律
从导弹飞行速度曲线可以看出，速度曲线可以近似为两段折线。在0-2.1126s，即第一级发动机工作期间，导弹在较高的切向加速度加速飞行，在第一级发动机停止工作，第二级发动机结束工作之前，导弹维持加速飞行，但加速度非常小，近似为0.5左右. 原因是推力是影响导弹切向加速度和速度的决定性因素，由于导弹飞行始终保持在较小的攻角，在所有切向力中，推力大小远远大于气动阻力和重力沿速度方向分量。因此，在导弹第一级发动机提供的大推力下先以很高的加速度加速飞行来达到一定速度，再在第二级发动机工作时维持或进一步提高速度，以保证预定的射程。
5.3.3 攻角变化规律
攻角和俯仰角与弹道偏角之间存在简单的几何关系：α=?-θ。因此，由初值可知，在导弹发射时，攻角α=0°，由于速度矢量下偏的角速度大于弹体偏转的角速度，导致弹体产生了不断增加的攻角，但导弹是俯仰静稳定的，因此产生的俯仰静稳定力矩使弹体攻角又不断减小，但这种调节过程存在一定滞后，因此攻角也会产生较大超调和较长调节时间的的振荡，又弹体受到阻尼力矩，所以攻角的变化收敛于平衡攻角。导弹是无控飞行的，相当于舵偏角?z?0o，所以平衡攻角α也必然为0°。
5.3.4 俯仰角速度变化规律
俯仰角速度的变化规律和攻角的变化有类似之处。因为根据微分方程，俯仰角速度的俯仰静稳定力矩和阻尼力矩二者决定，导弹静稳定力矩又取决与攻角的变化，由导弹攻角变化规律可知攻角的变化是先振荡后在收敛为零值。因此经过积分后，俯仰角速度也是
振荡的并且在攻角为零后收敛于常值。但其振荡曲线的相位落后于攻角曲线。
5.3.5 弹道倾角变化规律
由弹道倾角曲线可以看出，在导弹整个飞行过程中始终不断减小。导弹发射后的初始段，速度较小，攻角也很小，因此在速度的法向，升力和推力分量也很小，重力分量G 很大，合力为负，于是弹道倾角下降较快；随着攻角的增加，推力在速度法向的分量不断增大，升力也由于速度的增大而增大，这都使弹道倾角的减小变缓。因为初始段推力的扰动和升力存在振荡，因而弹道倾角减小的同时伴随轻微的振荡。2.1126s后，第二级发动机工作，推力减小，升力逐渐稳定，因此弹道倾角在近似恒定的变化率下下降。
5.3.5 俯仰角变化规律
从俯仰角变化曲线可以看出，其与弹道倾角变化趋势相似，在整个飞行过程中，
o
俯仰角整体趋势是不断减小的。导弹发射时两者同为18，由于随时间推移，弹体产生了攻角并不断加大，由于导弹是俯仰静稳定的，所以攻角又会出现收敛性振荡，结果是使导弹弹体跟随速度矢量偏转，从而产生了俯仰角的变化。然而随着弹体攻角收敛于0o,俯仰角的变化曲线也收敛于弹道倾角变化的曲线。
5.3.6 导弹质量变化规律
导弹质量变化规律较为简单，因为在两级发动机分别工作时间段，其秒流量几乎恒为常量，所以由发动机质量微分方程可知，其变化规律和最终曲线一样，是两段下降直线，即分段线性变化规律。
（四）结语
本次实验仅针对无控飞行弹道的仿真和分析，但其意义却很大，为以后有控弹道设计和仿真提供了必要的经验和数据。在现在导弹弹道设计中常常以理论弹道作为参考，再加之经验校正及优化便可得到实际弹道，所以它是很有辅助课程学习价值的。这次实验不仅加强了作者对基本概念的理解，也提升了对数学方法和软件编程应用的能力。所以这次设计对以后课程学习及未来可能的工作都有极大的帮助。
在本次实验中，岳振江同学，赵建博同学等都提出了很多建设性的意见和建议，冀四梅老师辛勤的耕耘是学生获取只是最大的源泉，特此一并表示感谢！
（五）参考文献
[1].钱杏芳，林瑞雄，赵亚男.导弹飞行力学[M].北京：北京理工大学出版社，2012 [2].丁丽娟，程杞元.数值计算方法[M].北京：高等教育出版社，2011
[3].刘浩.MATLAB R2012完全自学一本通[M].北京：电子工业出版社，2013 [4].史峰，邓森，陈冰.MATLAB函数速查手册[M].北京：中国铁道出版社，2011