铅垂平面飞行弹道仿真改进版讲解
过重力补偿比例导引弹道仿真研究
过重力补偿比例导引弹道仿真研究蒋 明, 宋卫东(军械工程学院兵器测试中心,河北 石家庄 050003)[摘 要] 激光半主动末制导炮弹能实现对目标的“攻顶”效应,分析了“攻顶”技术及其理论实质,即过重力补偿比例导引律的应用,以弹丸和目标的相对运动方程为基础,通过M atlab仿真计算的方法,分析了比例导引系数、弹丸与目标的速度比、过重力补偿值等主要参数的变化对弹道性能的影响,提出了导引律实现中注意的几个问题及改进制导系统的建议.[关键词] 导引规律; 比例导引法; 过重力补偿; 弹道; 仿真[中图分类号] TJ765.4 [文献标识码] AResearch on the Proporti ona l Nav i ga ti on w ith Co m pen s a ti on overGrav ity Tra jectory S im ul a ti onJ iang M ing, Song W eidong(W eapons Test Centre,O rdnance Engineering College,Shijiazhuang 050003,China)Abstract:Laser Se m i2active Ter m inal Hom ing Pr ojectile can actualize attack on the t op of the target.It is analyzed that the technique of t op2attack and its theoretical matter is the app licati on of the Pr oporti onal Navigati on with Compensati on over Gravity.Based on the relative move ment equati on gr oup of the p r ojec2 tile and the target,using Matlab t o make the si m ulati on and calculati on,it is analyzed that the affecti on of the maj or para meter’change t o the traject ory,such as p r oporti onal navigati on quotiety,the vel ocity ra2 ti o of p r ojectile and target,the compensati on over gravity,and s o on.Some p r oble m s in the actualizati on of guidance and s ome suggesti on on the i m p r ove ment of the guidance and contr ol syste m are offered in the paper.Keywords:guidance la w; p r oporti onal navigati on; compensati on over gravity; traject ory; si m ula2 ti on1 引 言比例导引(P N)是一种有效的寻的制导规律,广泛应用于自寻的战术导弹和制导炮弹的的末段制导[1].由于在实际应用中还要考虑到目标机动、许用过载、末段落角约束等条件的制约,因此须对比例导引加以改进以满足作战要求.以攻击工事和装甲目标的顶部(简称“攻顶”)为主要战术性能的激光半主动末制导炮弹在其弹道末制导段采用了一种论实质,以弹丸和目标的相对运动方程为基础,计算分析主要弹道参数对弹道性能的影响,提出制导系统设计和导引律实现中应注意的几个问题.2 末制导炮弹的“攻顶”技术激光半主动末制导炮弹的“攻顶”技术主要基于比例导引律的实现.2.1 比例导引在末制导中的实现制导炮弹的末制导段采用比例导引来提高命中精度,即弹丸速度矢量V 的转动角速度σ·正比于目标视线的转动角速度q ·,其数学表达式为:σ·=k ·q·(1) 弹丸头部的导引头是弹载核心部件,为一陀螺跟踪装置,它能测量出弹丸目标之间连线即目标视线的角速度,通过制导系统的运算和信号变换,形成与目标视线角速度成正比例的控制信号.在该控制信号的作用下导引头修正线圈产生修正力矩,使导引头陀螺向目标方向进动,驱使导引头光轴与目标视线间的误差角趋近于零.在制导过程中,该控制信号同时送入控制舵机,利用气动舵控制弹丸的飞行姿态从而实现弹丸的导引飞行,保证弹丸稳定跟踪目标.但是在比例导引的实现过程中还要考虑重力的影响,需要利用制导系统为弹丸提供重力补偿.2.2 “攻顶”技术的过程描述所谓的“攻顶”技术就是指在弹道的末制导段上给导引头引入重力补偿控制信号,使导引头产生附加进动,形成控制指令使弹丸以较大的落角攻击目标顶部的技术.该控制信号取自自动驾驶仪惯性陀螺产生的重力补偿信号,它给导引头施加一个附加力矩,驱使导引头的陀螺转子在铅垂平面内向下误进动,引起导引头产生反向调整的控制指令,与目标视线角速度形成的控制信号一起送入控制舵机完成比例导引.加入重力补偿信号之后其弹道特征为弹道在末制导段的前半段小幅上抬,后半段以较大的角度进行俯冲攻击,使弹丸的落角增大,实现“攻顶”[4].2.3 “攻顶”技术的理论实质由以上描述可知导引头测得的角速度信号包括两部分:一部分是目标视线角速度,一部分是重力补偿引起的附加进动角速度.测量值是这两个角速度的代数和.所以导引律(1)式应变为:σ·=k ·(q ·+w (t ))(2) 若令q ·g =q ·+w (t ),则式中q ·g 为导引头进动角速度,是测量值;q ·为目标视线角速度;w (t )为导引头附加进动角速度,本文仿真选取w (t )为一小值常量w g .考察弹丸的重力补偿情况可以发现,弹丸的惯导系统已经为弹丸提供了所需的重力补偿,而导引头附加进动角速度的产生相当于额外增加了一部分重力补偿值,故而重力补偿值比实际需要值大.而导引弹道的导引全过程又遵循比例导引,故可以称(2)式描述的导引律为过重力补偿比例导引律,w (t )视为过重力补偿信号.“攻顶”技术的理论实质即为过重力补偿比例导引的实践应用.3 弹道仿真的实现弹道仿真的实现分三个步骤[5]:建立弹道模型;建立仿真模型;改变参数,观察结果.(1)建立弹道模型d rd t=V T co sηT -V cos ηrd qd t=-V T sinηT +V sin ηq =ηT +σT q =η+σ(3) 由于假定弹丸处于理想控制状态,只研究其运动学特征,故选择弹丸和目标的相对运动方程(3)以及导引规律方程(2)为弹道模型.式(3)中各参数意义参见文献[3].(2)建立仿真模型根据方程(3),在Matlab 中建立仿真模型.(3)改变参数,观察结果改变各弹道参数,如初始弹道角、比例系数、弹丸目标速度比、过重力补偿值等,观察仿真结果,并进行相关分析.·27·战术导弹技术 TacticalM issile Technol ogy March,2008,(2)2012-05-18########################2012-05-18########################2012-05-18图2 c 变化时的弹道曲线图1 k 变化时的弹道曲线4 主要参数对弹道性能的影响本文主要分析了比例导引系数k 、弹丸与目标的速度比c 、过重力补偿值w g 等三个参数对弹道性能的影响.假定弹丸做等速飞行,目标做等速直线运动.4.1 比例导引系数k 对弹道性能的影响弹丸的末段速度相对较小,假定在200m /s 至400m /s 之间,而弹丸与目标的速度比c 较大,假定c =20.由图1可知:当k =1时末段弹道弯曲比较大,随着k 的增加,整个弹道的弯曲程度比较均匀,并且很明显制导初始段有一个小幅度上抬,后半段为俯冲弹道.通过分析可以得到以下结论:(1)弹道先上抬后俯冲的主要原因是由于引入了过重力补偿信号,导引头产生了附加进动角速度w g ,随着k 的增加,俯冲程度变大.(2)随着k 的增加,弹道初始段过载变大,弹末段过载趋近于某一常值.(3)弹道中间某点处将出现一个过载零点,这是因为视线角速度q ·是变化的,当q ·与w g 这两个角速度的代数和为零时,出现弹道过载零点,零点两侧过载方向相反.k 越大零点出现的时刻越早.4.2 弹丸与目标的速度比c 对弹道性能的影响由图2可以看出,弹丸与目标的速度比对弹道影响较大,c 越小弹道越弯曲;c 越大,弹道越光滑.c 对弹道的影响主要表现在:(1)c >15,弹道的弯曲程度变得均匀,表现为过载变化幅值不大;c <5,弹道初始段上抬程度变大,末段俯冲程度也变大,表现为弹道过载的急剧变化,例如图2中c =1.5时的情况.(2)c 越大,导引飞行的时间越短,如c =30时,t =10.5s ;c 越小,导引飞行的时间越长,如c=1.5时,t =31s .4.3 过重力补偿值w g 对弹道性能的影响过重力补偿值w g 为一小量,并且具有角速度的物理意义,所以在选取时取w g 为0(°)/s ～5(°)/s ,w g =0时即为重力补偿等于实际重力,其他即为过重力补偿的情况.由图3可以得出如下结论:(1)w g 越小,弹道越平直;w g 越大,弹道越弯曲,弹道的上抬和俯冲越明显,对目标的“攻顶”效果越好.(2)w g ≤2(°)/s ,弹道过载较小,当w g =0时,弹道末段过载趋于0;w g >2(°)/s ,弹道过载变大,弹道中间有过载零点.5 结 论通过以上分析和仿真计算可知,比例导引系数k,弹丸和目标的速度比c 以及过重力补偿值w g 是影响弹道性能的三个主要参数.因此在导引律的实·37·战术导弹技术 TacticalM issile Technol ogy March,2008,(2)2012-05-18########################2012-05-18########################2012-05-18图3 w g 变化时的弹道曲线现和设计弹丸的制导控制系统时,需要注意以下几个方面:(1)弹上的制导系统参数的选择要考虑其对弹道的综合影响,例如k 增大使得弹道中段和末段过载变化均匀,但加大了初始过载值;w g 的增大使得“攻顶”效果显著同时加大了全弹道上的过载.(2)由于弹丸末制导段速度变化不大,故c 对弹道性能的影响反应了弹丸对目标运动状态的选择和要求.因此在设计制导控制系统时要考虑其对目标的适用范围. (3)本文的计算主要基于k 和w g 在全弹道上为定量,但是制导控制系统的自适应和智能化的发展要求实现最优控制,故可以进一步改进导引律,使比例导引系数k 和过重力补偿值w g 随着弹道需求的变化而变化.(4)过重力补偿比例导引和传统比例导引一样过分依靠目标视线角速度q ·这一个观测量来控制,而在目标机动性加大或者出现测量误差时,难以实现精确控制而导致脱靶,因此可以发展新的比例导引律,使其能基于更多的观测量实现精确制导.[参 考 文 献][1] 雷虎民,刘兴堂,陈新海.比例导引的指令加速度[J ].电光与控制,1999(1):19～24.[2] 陈宏,雷鸣.激光寻的制导导引头技术[J ].光电子技术,2002,27(1):53～57.[3] 钱杏芳,张鸿端,林端雄.导弹飞行力学[M ].北京工业学院出版社.1987,12.[4] 翁彦,张天桥.激光半主动寻的制导在滚转导弹上的应用[J ].弹箭与制导学报,2002(1).[5] 张中南,童幼堂,张卫峰.比例导引法导引弹道仿真研究[J ].战术导弹技术,2005(2):56～59.(上接第48页)5.2 复合打击战法研究的内容对复合目标进行打击时,哪种型号的弹先打,哪种型号的弹后打,其总体毁伤效果是不一样的.复合式打法研究的内容还包括:导弹武器组合优化分析,即根据目标类型特性、战斗部毁伤效应和毁伤指标,科学组合各类型导弹对目标实施最有效的打击,以获取最大的军事效益;武器打击顺序研究,即根据目标价值、作战部队目前状态,结合作战任务,确定最佳的武器打击顺序;此外还包括多型号武器打击组合效应研究、耗弹量分析与优化及最佳复合打击方案的确定等.[参 考 文 献][1] 戴开政,王志平.伊拉克战争述评[M ].桂林:桂林陆军学院图书馆,2003.[2] 杨世荣.目标分析[M ].第二炮兵工程学院,2002,6:67.[3] 李新其等.美军在近期几场实战中目标选择的特点及对二炮的启示.二炮军事学术,2005(1):69～71.[4] 刘怡昕等.论对复合目标的复合式打法及其射击效率.射击学报,2004(2):4～6.[5] 章克凌,汪新红,张伟.精确制导武器对阵地工程损伤评估研究现状.科技研究,2002(5):4～6.·47·战术导弹技术 TacticalM issile Technol ogy March,2008,(2)2012-05-18########################2012-05-18########################2012-05-18file:///C|/Users/Administrator/Desktop/新建文本文档.txt通信/电子电脑、 杂志、 会议、 劳动合同、 生活休闲、 考试、 股票。
铅垂平面飞行弹道仿真改进版
数值实验作业1 实验背景实验名称：铅垂平面飞行弹道仿真及分析实验内容与要求：根据描述飞行器在铅垂平面内运动的数学模型，编制某导弹的铅垂平面无控飞行弹道仿真程序,利用计算机解算初始段无控飞行弹道,对初始段弹道参数的变化规律进行分析.2 建立数学模型：3计算方法研究确定数值积分方法和积分步长使用算法：四阶龙格—库塔法。
使用步长：h=0。
0054 原始数据：1)初值=0 x=0（m) y=20。
0（m）ϑ=18︒θ=18︒v=20(m/s) ωz （rad/s）m=52。
38（kg)2)攻角与马赫数范围(仅用于插值计算)攻角=0～10︒马赫数=0.1～0。
95）推力数据第一级工作结束时间：2。
1126s，第二级工作结束时间：44。
0832s 6)发动机质量秒流量7）转动惯量8）导弹重心（起自头部）9）静稳定力矩系数00Xg Xg z m =αα当导弹重心变化时的修正公式：L Xg Xg c m m y z z /)(001-+=αααα10）阻尼力矩导数z z m ϖ使用的插值算法:气动数据插值——等距双变元抛物线插值; 推力、重心、转动惯量等-—不等距一元线性插值。
5 空气动力和空气动力矩表达式22212121()2zz xy z z z z z z zX c V S Y c V SM M M m m V SL ϖϖααρραϖαϖρ===+=+6 编制计算程序计算机算法采用了matlab 实现,源程序见附件.该程序有八个函数组成，各函数之间的调用关系如下图所示.1) 子函数initl 的功能是输入求解导弹运动方程组所需的原始数据.2) 子函数rk_4是四阶龙格-库塔法积分算法子函数，其中调用了子函数dery 。
3) 子函数dery 功能是计算微分方程组的右端函数,其中调用了子函数interp.4)插值子函数interp 功能是计算所有需要插值的参数,其中调用了子函数interp11和interp33。
某型飞机垂尾工作平台仿真计算及优化
民用飞机在航线运营时，在一定飞行小时或一定循环架次后，需进行各项检查、维修和维护任务。
如何提高飞机的维修效率、降低飞机的维修时间是国内外各家航空维修公司的目标。
民用飞机大型工作平台不仅可以方便飞机某些位置的检查和维修，也可以保证飞机几个部位同时进行维修，这就可以大大降低飞机维修时间和成本，提高经济效益。
由于大型工作平台多用于飞机局部难以达到的区域，而飞机某些部位的高度需要达到5 m甚至更高的高度，这就使得某些大型工作平台高度、重量等参数比一般设备大。
这些设备在使用的时候会存在人员和飞机方面的安全问题。
设备的庞大，造成对设备整体结构强度、移动等方面要求较高。
对于所有大型工作平台来说，在结构设计时，平台的设DOI：10.16660/ k i.1674-098X.2016.21.010某型飞机垂尾工作平台仿真计算及优化汪邦中（上海飞机设计研究院四性与产品支援设计研究部 上海 200235）摘 要：该文介绍了民用飞机大型工作梯类地面设备在飞机维修时可以提高维修效率、降低维修成本和保证飞机某些部位的维修可达等功能，并提出了工作梯仿真计算的必要性。
该文以某型飞机垂尾工作平台为例，利用有限元方法对大型工作平台整体建立有限元模型，通过有限元仿真计算方法对其整体结构强度进行分析。
根据分析结果，对工作平台的局部结构给出优化建议。
关键词：地面设备 工作平台 仿真计算中图分类号：V21文献标识码：A文章编号：1674-098X(2016)07(c)-0010-04图1 某型飞机垂尾工作平台三维图材料牌号密度(g/cm 3)弹性模量E (GPa)泊松比屈服极限(MPa)R P0.2强度极限(MPa)R m 6061-T4 2.73680.33110180Q2357.852000.3235>375表1 材料属性. All Rights Reserved.计刚度应达到要求，保证平台在使用时不碰到飞机结构造成结构损伤；抽拉板的刚度达到要求，保证维修人员安全；对平台结构进行整体强度校核，保证工作平台整体安全。
飞行器姿态控制方法综述
飞行器姿态控制方法综述一.引言经过一个世纪的发展，各种飞行器如雨后春笋般出现，从飞机、导弹到火箭、卫星，从宇宙飞船、航天飞机、空间站到月球探测器、火星探测器。
这些飞行器能在空中按预定的轨迹运动总离不开它的姿态控制系统，飞行器在空间的运动是十分复杂的。
为使问题简单化，总是将一飞行器的空间运动分解为铅锤平面的纵向运动和水平面内的侧向运动，将飞行器在空间的角运动分解成俯仰、偏航和滚动三个角运动。
由于角运.动使飞行器的姿态发生变化，所以对角运动的控制就是对飞行器姿态的控制。
对于飞行器姿态的控制，不同的飞行器需要不同的策略，本文主要就飞行器姿态控制方法的应用与发展作一一论述。
二.姿态控制的数学模型要控制飞行器的姿态，就是要控制使飞行器三个姿态角发生变化的力矩大小。
飞行器的姿态模型可以认为是一类不确定MIMO 仿射非线性系统，如式(1)所示：()//()//()//(c o s s i n )/c o s c o s s i ns i n t a n c o s t a n x y z y x x x x x z x x x y y y x x y x y z z z x x x z x y z I I I M I I I I M I I I I M I ωωωωωωωωωψωθωθϕϕωθωθθωθϕωωθϕ=-+⎧⎪=-+⎪⎪=-+⎪⎨=-⎪⎪=+⎪=+-⎪⎩ （1） 式中，x 、y 、z 下标表示空间飞行器的三个主轴方向；I 表示相对于飞行器质心的惯量矩，设飞行器是主轴对称的，则惯量积可以忽略；ω表示飞行器相对于惯性空间的角速度;M 表示控制力矩；,,ψϕθ分别是飞行器的欧拉角。
控制了M 的大小，就可以控制飞行器按我们期望的轨迹运动。
M 由飞行器上的执行机构产生，常见的有空气舵、推力矢量发动机、反作用飞轮、喷气执行机构或由其它环境力执行机构。
三.飞行器姿态控制方法3.1空气动力控制根据运动的相对性原理和气体流动时的基本定律，当飞行器在大气中以一定的速度飞行时，飞行器都会受到空气动力的作用。
(518618873) 典型导引规律三维弹道仿真分析
总第164 期200 8年第2期舰船电子工程Vo l.28No. 211 0Sh i p E lec tron ic En gin ee ringt3典型导引规律三维弹道仿真分析周纪元 1 ) 童幼堂2) 张磊1) 王亚1)(海军大连舰艇学院研究生1队1 )大连11 601 8)( 海军大连舰艇学院舰载武器系2 )大连11 601 8)摘要介绍几种典型导弹导引规律及其约束方程,运用MA T L AB 语言对纯追踪法、平行接近法和比例导引法的理想弹道进行了三维仿真。
绘制导弹速度和目标速度变化时,其对应的三维理想弹道,并进行比较分析。
关键词导引规律;弹道;仿真; MA T LAB中图分类号TP391. 91 引言导弹在攻击目标时,为了使导弹稳定的跟踪并命中目标,导引系统必须及时准确地提供导引信号,而导引信号又必须严格遵循某种规律,这种规律就是导引规律。
导引规律一般分为两大类,即两点法和三点法。
两点法又包括纯追踪法、比例导引法和平行接近法及固定前置角法等[ 1 ] 。
目前,大多数文献给出的两点法的运动学弹道多为导弹和目标在同一平面内的二维弹道,本文使用 MA T LAB 语言,对两点法前三种导引规律的理想弹道进行了三维仿真。
绘制出了三维仿真弹道,分析比较了不同导引规律对应的导弹和目标运动参数变化时的弹道变化特性。
2 两点法导引规律式中: r = d r/ d t:距离变化率;θ: V 与目标视线的夹角, 即目标运动前置角;β: Vm与目标视线的夹角;即导弹运动前置角;φ: 目标视线与x 轴夹角, 即视线角; r:目标与导弹之间的距离。
如图1所示。
图1导弹与目标相对关系图由图1可得各夹角之间的相互关系为:φ= η+θ (3)φ= β+γ(4)将( 3) 、( 4 )式代入(1)、( 2 ) 式可得:r = V t co s(φ- η) - V m co s(φ- γ) (5)两点法是确定导弹与目标两点在空间相对位置的方法。
铅垂面内弹道设计、解算与分析
铅垂面内弹道设计、解算与分析檀望春，李禹志，李坤猷，陈国勇北京理工大学飞行器设计与工程专业摘要：基于MATLAB的对铅垂面内弹道设计的数值模拟仿真。
在给定的参数条件下设计合理的飞行方案并研究导弹在不同飞行状态下（爬升段、平飞段、俯冲段）攻角、速度等的变化及其飞行轨迹。
关键词：数值模拟、飞行轨迹、铅垂面。
The Design, Analysis and Solver of Vertical Plane TrajectoryTAN Wang-chun, LI Yu-zhi, LI Kun-you, CHEN Guo-yongAbstract: Design the vertical plane trajectory which based on MATLAB simulation. Under the conditions of the given parameters design the missile flight plan and study with the angle of attack, speed of change and its flight path at different flight conditions (climb segment, level flight segment, dive segment) .Keywords: Numerical simulation, Flight path, Vertical plane.导弹作为现代战争中最重要的武器之一一直备受重视。
导弹以它的精确制导打击目标和强大的破坏性能改变整个战场的格局给现代战争带来了深刻的影响。
导弹的运动以及控制是十分复杂的。
研究导弹的空间飞行状况首先应当研究导弹的弹道学方面的问题。
将导弹的运动看成质点的运动并且假定导弹的制导系统是完全理想的。
在导弹飞行的各个状态下导弹所受外力的和力矩为零。
因此，在质点上作用重力、推力、空气动力等，研究这些力和导弹的运动之间的关系，加上制导系统理想工作的约束关系式，就可以求出导弹质心的运动轨迹——弹道、飞行速度和过载等的参数。
飞行器控制课程设计报告书
航天学院导弹制导控制系统设计与Matlab仿真班级：02020902小组成员：导弹的制导系统设计与仿真敌方坦克正在自东向西作水平匀速直线运动，VT=15m/s，我方反坦克步兵使用反坦克导弹进行反击，假设导弹速度飞行中近似不变，VM=250m/s，开始引导瞬间发射点距离目标RT0=2000m，导弹距离发射点RM0=20m，且qM0=qT0=80○，试利用比例导引法对该导弹的攻击过程进行仿真，并绘制弹道曲线和导弹的法向过载曲线。
1）第一步：选取坐标系和攻击平面通常步兵使用反坦克导弹攻击坦克时，发射点和坦克基本在同一水平面，并且采用比例导引法的导弹整个飞行过程也几乎在这个水平面，因此攻击平面主要考虑为发射点和坦克所在的水平面，而导弹在纵向铅垂平面的运动可以近似认为是等高度水平飞行。
另一方面，步兵使用的反坦克导弹射程一般不超过 5km ，因此可以认为地面是水平大地且不考虑地球自转等影响。
因此可以选取地面坐标系Axyz ，原点 A 与制导站（发射点）重合，Az 轴平行与目标运动方向且指向正东，Az 轴指向正北， Ay指向天向，攻击平面为 Axz平面（如下图所示）。
图表1反坦克导弹攻击坦克时的坐标系和攻击平面选取2）第二步：导弹的受力分析和矢量分解现在分析导弹的受力情况，通常导弹会受到重力、气动力（升力、侧向力和阻力）、发动机推力的作用。
对于本例题来说，由于反坦克导弹攻击坦克主要在攻击水平面，因此可以近似认为在铅垂方向上，重力和升力相互平衡，因此不再考虑；而对于阻力，通常其与速度方向相反，会减小导弹的速度，在本例题中为了简单期间，暂时不考虑阻力的影响；对于发动机推力，由于反坦克导弹的发动机通常在发射初期就已经燃烧完毕，因此仿真阶段为无动力的自由飞行阶段，发动机推力为零。
由于导弹主要在水平面飞行，因此可以认为导弹的速度也在水平面，所以重力不会影响到导弹的速度。
另一方面前面假设导弹的阻力和推力为零时，因此可以认为导弹在攻击过程中的速度恒定不变。
水下飞行器的弹道仿真
深度 要求 的提 高，飞行 器在 水 中段 的运 动 时间增 长 ， 对于 静不稳 定 的飞行器 来 说 ，由于入水 攻角 的影 响 ， 其姿态 是不 断发 散的 ； 同样地 ，提 高发 射深度 后 ， 由
于水 中段行 程加长 ， 飞行器 出水速 度 降低 ，使 得 飞行 器 出水 时间加长 ，导致 其 出水姿 态进一 步加 大 。
２０ ０ ９年 第 ６ 期 总第 ３ ４ ０ 期
文 章编号 ：１ ０ —１ ２２ ０ ） ６０ ０ －４ ０ ４７ ８ （０ ９ ０ ・０ ｌ ０
导 弹 与 航 天 运 载 技 术
ＭＩ Ｓ ＬＥ Ｓ Ｉ Ｓ ＡＮＤ Ｐ Ｓ ＡＣＥ ＶＥＨＩ ＥＳ ＣＬ
Ｎ Ｏ６ ２ ９ ． ００ Ｓ ｕｍ Ｏ．０ Ｎ ３４
１ 铅 垂 面 内无控 弹 道 的仿 真
１ 无控 弹道模 型的建 立 ． １ 根 据 飞行器 的受力 情况 ，水 弹道 可分为 筒 中段 、 水 中段 、 空 中段 。本文 研 究飞行器 在铅 垂面 内的运 动 情况 ，具 体分 为水 中段 和 出水段两 部 分 ，建 立水下 飞
行器 的弹道模 型：
真 ，根 据飞 行器无 控运 动时 的合外 力矩 ，增加 了控 制
力矩 ，建立 了水下 飞行 器有控 弹道 的数 学模 型 ，并对
的性 能 。 冷发射 状态 下的 飞行器 在水 中运动 时 ， 由于 发射
＝ ［ 一 ＋ ｇ（）］ ｉ （ ） 互 ＋ Ｆ ｓ 南 。 十 ＋Ｇ —ｉ ｎ
（ ＋
（－ Ｉ
Ｉ ［ 也（＾ 一 ＋ ＋ 一 ： ＋ｔ 五 一 ） 一 Ｃ ｍ ｇ
ｃｇ ＧｃｓＪ ＋ ） ￣ 匠 一 ｏ
．
十米 ，但 由于环 境和 起始 条件都 很苛 刻 ，涉 及 的 问题 比较 复杂 ，它 的设计 、研 究和试 验在 水下 飞行器 研制
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值实验作业1 实验背景实验名称：铅垂平面飞行弹道仿真及分析实验内容与要求：根据描述飞行器在铅垂平面内运动的数学模型，编制某导弹的铅垂平面无控飞行弹道仿真程序，利用计算机解算初始段无控飞行弹道，对初始段弹道参数的变化规律进行分析。
2 建立数学模型：d3计算方法研究确定数值积分方法和积分步长使用算法：四阶龙格—库塔法。
使用步长：h=0.0054 原始数据：1）初值=0(rad/s)x=0(m) y=20.0(m) ϑ=18︒θ=18︒v=20(m/s) ωzm=52.38(kg)2）攻角与马赫数范围(仅用于插值计算)攻角=0~10︒马赫数=0.1~ 0.935）推力数据第一级工作结束时间：2.1126s，第二级工作结束时间：44.0832s 6）发动机质量秒流量7）转动惯量8）导弹重心（起自头部）9）静稳定力矩系数00Xg Xg z m =αα当导弹重心变化时的修正公式：L Xg Xg c m m y z z /)(001-+=αααα10）阻尼力矩导数z z m ϖ使用的插值算法：气动数据插值——等距双变元抛物线插值； 推力、重心、转动惯量等——不等距一元线性插值。
5 空气动力和空气动力矩表达式22212121()2zz xy z z z z z z zX c V S Y c V SM M M m m V SL ϖϖααρραϖαϖρ===+=+6 编制计算程序计算机算法采用了matlab 实现,源程序见附件.该程序有八个函数组成,各函数之间的调用关系如下图所示.1) 子函数initl 的功能是输入求解导弹运动方程组所需的原始数据.2) 子函数rk_4是四阶龙格-库塔法积分算法子函数,其中调用了子函数dery. 3) 子函数dery 功能是计算微分方程组的右端函数,其中调用了子函数interp.4)插值子函数interp 功能是计算所有需要插值的参数,其中调用了子函数interp11和interp33.5) 子函数interp11是不等距单变元线性插值函数,主要用于转动惯量Jz, 质心位置xg, 推力p 等单变量的插值.6) 子函数interp33是等距双变元抛物线-线性插值函数,主要用于气动力系数cx, cy, 气动力矩系数z m α,z z m ϖ,等双变元参数的插值,该函数调用了子函数interp31.7) 子函数interp31是等距单变元抛物线插值函数,被interp33调用,完成子函数interp33的等距双变元抛物线-线性插值功能.关于程序中出现的数组和变量名,作如下说明: acx:阻力系数数组; acy: 升力系数数组;amzaf:某一质心位置下的静稳定性导数数组; amzwz:阻尼力矩系数导数数组;axg:质心位置Xg 随时间的变化规律; ajz:转动惯量Jz 随时间的变化规律; ap:起飞,续航发动机的推力值;amc: 起飞,续航发动机的燃料质量的秒流量值; agc:质心位置变化的始末值;andm:气动数据插值所需的Ma 数的最小,最大值; andaf: 气动数据插值所需的攻角α的最小,最大值;y:存放积分结果的数组,该数组在程序开始时存放积分初值; dy:存放右端函数数值的数组;b:存放三个时间值的数组,其中b(2)存放起飞发动机工作结束时间;b(3)存放续航发动机工作结束时间;程序中一些主要的变量名有: L —特征长度; S —特征面积; SONIC —声速C; RHO —大气密度ρ;h —积分步长;程序中其他变量都是存放中间结果的变量.7 计算结果运行程序得到弹道曲线, 速度曲线, 攻角曲线如下:通过运用龙格库塔法和一些插值方法求解导弹运动方程组,获得了导弹各运动参数的变化规律.通过该算法获得的变化规律比较接近真实情况.附件:源程序function main()clear allclear globalglobal y;global ii;ii=0;h=0.005;initl();while y(7)>=0ii=ii+1;result(ii);rk_4(8,h);endsavedata(ii);drawing();%原始数据初始化function initl()global acx acy ajz amzaf amzwz axg ap amc agc andm andaf b L S SONIC RHO;global y;y=[0 20. 18 0 18 0 20 52.38];%马赫数ma%攻角alpha%三个时间b(1)为导弹离轨时间b(2)为起飞发动机工作结束时间b(3)为续航发动机工作结束时间b=[0 2.1126 44.0832];%系数表维数n1=9;n2=6;%andm最小,最大值andm(1)=0.1;andm(2)=0.9;%andaf最小,最大值andaf(1)=0;andaf(2)=10;%阻力系数acx=[.4177 .4404 .5219 .6603 .8534 1.1023;.3858 .4086 .4903 .6290 .8226 1.0723;.3779 .4007 .4827 .6218 .8160 1.0666;.3785 .4015 .4838 .6234 .8184 1.07;.3787 .4018 .4846 .6249 .8209 1.0738;.3829 .4062 .4897 .6310 .8248 1.0835;.3855 .4091 .4934 .6363 .8358 1.0938;.4082 .4321 .5175 .6621 .8641 1.1254;.4947 .5192 .6073 .7571 .9672 1.2392];%升力系数acy=[.0000 .6430 1.4758 2.2870 3.0713 3.8463;.0000 .6454 1.4807 2.2942 3.0915 3.8731;.0000 .6480 1.4858 2.3014 3.0915 3.8731;.0000 .6512 1.4923 2.3107 3.1039 3.8891;.0000 .6554 1.5007 2.3227 3.1197 3.9092;.0000 .6617 1.5134 2.3409 3.1426 3.9401;.0000 .6698 1.5304 2.3661 3.1775 3.9835;.0000 .6792 1.5501 2.3950 3.2162 4.0323;.0000 .6933 1.5935 2.4706 3.3273 4.1790];%推力ap=[.000 .15 .49 2.11 2.27 3.53 8.78 25.45 42.80 43.68 44.08;3245.76 6020.14 4952.92 5956.44 476.77 430.906 411.698 401.8 399.84 399.742 21.756]; %发动机质量秒流量amc=[0. 2.1 2.105 44.1 44.105 100;2.362 2.362 0.21059 0.21059 0.0000 0.0000];%转动惯量ajz=[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];%导弹重心axg=[.0 2.0 2.4 10.0 18.0 26.0 32.0 38.0 42.0 44.0;.9381 .9095 .9091 .9026 .8969 .8928 08907 .8896 .8895 .8896];agc=[.9381 .8896];%静稳定力矩系数amzaf=[0.0000 -0.0104 -0.0341 -0.0564 -0.0771 -0.0985;0.0000 -0.0104 -0.0341 -0.0564 -0.0770 -0.09830.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];%阻尼力矩导数%当xg=.9381时amzwz=[-0.4686 -0.4829 -0.4928 -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.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时%其他参数S=0.0227; L=1.8; SONIC=343.13; RHO=1.225;%四阶龙格-库塔法子函数function rk_4(n,h)global y dy;dy=zeros(n,1);old_y=zeros(1,n);y1=zeros(1,n);a(1)=h/2;a(2)=h/2;a(3)=h;a(4)=h;dery(y);for i=1:nold_y(i)=y(i);endfor j=1:3for i=1:ny1(i)=old_y(i)+a(j)*dy(i);y(i)=y(i)+a(j+1)*dy(i)/3;enddery(y1);endfor i=1:ny(i)=y(i)+a(1)*dy(i)/3;end%右端子函数function dery(y)global dy ;global L S SONIC RHO;global ma abs_alpha cx cy mzaf mzwz jz alpha p mc;aa=zeros(1,4);q=RHO*y(2)*y(2)/2;ma=y(2)/SONIC;alpha=y(5)-y(3);abs_alpha=abs(alpha);interp();aa(1)=sin(y(3)*pi/180);aa(2)=cos(y(3)*pi/180);aa(3)=sin(alpha*pi/180);aa(4)=cos(alpha*pi/180);if alpha<0cy=-cy;endxf=cx*q*S;yf=cy*q*S;wzt=y(4)*L/y(2);dy(1)=1;dy(2)=(p*aa(4)-xf-9.81*y(8)*aa(1))/y(8);dy(3)=(p*aa(3)+yf-9.81*y(8)*aa(2))/(y(2)*y(8))/pi*180;%转换成度dy(4)=(mzaf*alpha+mzwz*wzt)*q*S*L/jz;dy(5)=y(4)/pi*180;dy(6)=y(2)*aa(2);dy(7)=y(2)*aa(1);dy(8)=-mc;%插值子函数function interp()global acx acy ajz amzaf amzwz axg ap amc agc andm andaf b L ; global ma abs_alpha cx cy mzaf mzwz xg jz p mc;global y;if y(1)<b(3)%a=interp11(tt,3,y(1));xg=interp11(axg,10,y(1));jz=interp11(ajz,14,y(1));p=interp11(ap,11,y(1));mc=interp11(amc,6,y(1));else xg=axg(2,10);jz=ajz(2,14);p=ap(2,11);mc=amc(2,6);endcx=interp33(ma,abs_alpha,9,6,andm,andaf,acx);cy=interp33(ma,abs_alpha,9,6,andm,andaf,acy);mzaf0=interp33(ma,abs_alpha,9,6,andm,andaf,amzaf);if abs_alpha~=0.mzaf=mzaf0+cy*(xg-agc(1))/(abs_alpha*L);%去掉了*RAD else mzaf=mzaf0;endmzwz=interp33(ma,abs_alpha,9,6,andm,andaf,amzwz);%不等距单变元线性插值子函数function res=interp11(yy,n,x)for j=1:(n-1)if x<=yy(1,j+1)i=j;break;else i=n-1;endendres=yy(2,i)+(yy(2,i+1)-yy(2,i))*(x-yy(1,i))/(yy(1,i+1)-yy(1,i));%等距双变元抛物线线性插值子函数function res=interp33(x,qq,n1,n2,a,bb,yy)h=(bb(2)-bb(1))/(n2-1);i=fix((qq-bb(1))/h+1);if (i-1)<0i=1;elseif (i-n2)>=0i=n2-1;endyy1=interp31(x,n1,i,a,yy);yy2=interp31(x,n1,i+1,a,yy);res=((qq-(i-1)*h-bb(1))*yy2-(qq-i*h-bb(1))*yy1)/h;%等距单变元抛物线插值子函数function res=interp31(x,n1,i,a,yy)h=(a(2)-a(1))/(n1-1);i1=fix((x-a(1))/h+1);if (i1-1)<=0i1=1;else if (i1-n1)>=0i1=n1-2;endendx0=a(1)+(i1-1)*h;x1=x0+h;x2=x0+2*h;c0=0.5*(x-x1)*(x-x2)/(h^2);c1=-(x-x0)*(x-x2)/(h^2);c2=0.5*(x-x0)*(x-x1)/(h^2);res=c0*yy(i1,i)+c1*yy(i1+1,i)+c2*yy(i1+2,i);%保存节点处数据function result(k)global R1 R2 R3 R4 R5 R6 R7 R8 R9;global y;R1(k)=y(1);R2(k)=y(2);R3(k)=y(3);R4(k)=y(4);R5(k)=y(5);R6(k)=y(6);R7(k)=y(7);R8(k)=y(8);R9(k)=y(5)-y(3);%输出数据文件子函数function savedata(ii)global R1 R2 R3 R4 R5 R6 R7 R8 R9;fp=fopen('c:\\data.dat','w');fprintf(fp,'T V SIGOMA WZ E X Y M ALPHA\n');for i=1:iifprintf(fp,'%f %f %f %f %f %f %f %f %f\n',R1(i),R2(i),... R3(i),R4(i),R5(i),R6(i),R7(i),R8(i),R9(i));endfclose(fp);%绘图子函数function drawing()global R1 R2 R3 R4 R5 R6 R7 R8 R9;figure('numbertitle','on','name','弹道曲线')plot(R6,R7);title('弹道曲线')xlabel('x/m');ylabel('y/m');figure('numbertitle','on','name','速度曲线')plot(R1,R2);title('速度曲线')xlabel('T/s');ylabel('V/(m/s)');figure('numbertitle','on','name','攻角曲线')plot(R1,R9);title('攻角曲线')xlabel('T/s');ylabel('alpha/度');。