铅垂平面飞行弹道仿真改进版

数值实验作业1实验背景实验名称：铅垂平面飞行弹道仿真及分析实验内容与要求：根据描述飞行器在铅垂平面内运动的数学模型，编制某导弹的铅垂平面无控飞行弹道仿真程序，利用计算机解算初始段无控飞行弹道，对初始段弹道参数的变化规律进行分析。2建立数学模型：cossinsincoscossinzzzzzzzzzdVmPXGdtdmVPYGdtdJMMdtddtdxVdtdyVdtdmmdt3计算方法研究确定数值积分方法和积分步长使用算法：四阶龙格—库塔法。使用步长：h=0.0054原始数据：1）初值x=0(m)y=20.0(m)=18=18v=20(m/s)z=0(rad/s)m=52.38(kg)2）攻角与马赫数范围(仅用于插值计算)攻角=0~10马赫数=0.1~0.93）阻力系数表马赫数攻角()02468100.1.4177.4404.5219.6603.85341.10230.2.3858.4086.4903.6290.82261.07230.3.3779.4007.4827.6218.81601.06660.4.3785.4015.4838.6234.81841.07000.5.3787.4018.4846.6249.82091.07380.6.3829.4062.4897.6310.82841.08350.7.3855.4091.4934.6363.83581.09380.8.4082.4321.5175.6621.86411.12540.9.4947.5192.6073.7571.96721.23924）升力系数表马赫数攻角()02468100.1.0000.64301.47582.28703.07133.84630.2.0000.64541.48072.29423.08143.85980.3.0000.64801.48582.30143.09153.87310.4.0000.65121.49232.31073.10393.88910.5.0000.65541.50072.32273.11973.90920.6.0000.66171.51342.34093.14363.94010.7.0000.66981.53042.36613.17753.98350.8.0000.67921.55012.39503.21624.03230.9.0000.69331.59352.47063.32734.17905）推力数据t(s).000.15.492.112.273.538.7825.4542.8043.6844.08P(kgf)331.2614.3505.4607.848.6543.9742.0141.0040.8040.792.22第一级工作结束时间：2.1126s，第二级工作结束时间：44.0832s6）发动机质量秒流量t(s)0.2.12.10544.144.105100秒流量(kg/s)2.3622.3620.210590.210590.0.7）转动惯量t(s).02.02.46.410.414.418.422.426.430.434.038.442.444.0Jz(kgms)8.357.887.867.817.787.757.737.717.707.707.697.697.697.698）导弹重心（起自头部）t(s).02.02.410.018.026.032.038.042.044.0XG(m).9381.9095.9091.9026.8969.8928.8907.8896.8895.88969）静稳定力矩系数00XgXgzm马赫数攻角()02468100.10.0000-0.0104-0.0341-0.0564-0.0771-0.09850.20.0000-0.0104-0.0341-0.0564-0.0770-0.09830.30.0000-0.0104-0.0341-0.0564-0.0769-0.09820.40.0000-0.0105-0.0342-0.0564-0.0768-0.09790.50.0000-0.0104-0.0339-0.0560-0.0761-0.09690.60.0000-0.0093-0.0314-0.0521-0.0708-0.09030.70.0000-0.0080-0.0286-0.0477-0.0650-0.08290.80.0000-0.0065-0.0252-0.0425-0.0578-0.07390.90.0000-0.0053-0.0229-0.0391-0.0538-0.0693当导弹重心变化时的修正公式：LXgXgcmmyzz/)(00110）阻尼力矩导数zzm当Xg=.9381时马赫数攻角()02468100.1-0.4686-0.4829-0.4982-0.5130-0.5272-0.54090.2-0.4707-0.4850-0.5003-0.5150-0.5292-0.54290.3-0.4744-0.4886-0.5039-0.5186-0.5327-0.54640.4-0.4797-0.4939-0.5090-0.5237-0.5378-0.55140.5-0.4882-0.5022-0.5173-0.5318-0.5458-0.55930.6-0.5089-0.5227-0.5376-0.5520-0.5658-0.57910.7-0.5366-0.5502-0.5649-0.5790-0.5927-0.60580.8-0.5738-0.5871-0.6014-0.6153-0.6287-0.64150.9-0.6272-0.6407-0.6553-0.6694-0.6830-0.6960当Xg=.8896时马赫数攻角()02468100.1-0.6179-0.6384-0.6600-0.6805-0.6999-0.71820.2-0.6207-0.6410-0.6626-0.6830-0.7024-0.72070.3-0.6253-0.6455-0.6670-0.6874-0.7067-0.72490.4-0.6319-0.6521-0.6734-0.6937-0.7129-0.73100.5-0.6424-0.6624-0.6835-0.7036-0.7226-0.74060.6-0.6669-0.6866-0.7074-0.7272-0.7459-0.76360.7-0.6997-0.7190-0.7395-0.7589-0.7774-0.79480.8-0.7435-0.7624-0.7824-0.8014-0.8194-0.83650.9-0.8069-0.8266-0.8474-0.8672-0.8859-0.903511）其它参数特征面积S(m2)特征长度L(m)毛翼展(m)音速SONIC(m/s)大气密度(kg/m3)0.02271.80.5343.131.225使用的插值算法：气动数据插值——等距双变元抛物线插值；推力、重心、转动惯量等——不等距一元线性插值。5空气动力和空气动力矩表达式22212121()2zzxyzzzzzzzXcVSYcVSMMMmmVSL6编制计算程序计算机算法采用了matlab实现,源程序见附件.该程序有八个函数组成,各函数之间的调用关系如下图所示.子函数initl主函数main子函数rk_4子函数result子函数savedata子函数drawing子函数dery子函数interp子函数interp33子函数interp31子函数interp111)子函数initl的功能是输入求解导弹运动方程组所需的原始数据.2)子函数rk_4是四阶龙格-库塔法积分算法子函数,其中调用了子函数dery.3)子函数dery功能是计算微分方程组的右端函数,其中调用了子函数interp.4)插值子函数interp功能是计算所有需要插值的参数,其中调用了子函数interp11和interp33.5)子函数interp11是不等距单变元线性插值函数,主要用于转动惯量Jz,质心位置xg,推力p等单变量的插值.6)子函数interp33是等距双变元抛物线-线性插值函数,主要用于气动力系数cx,cy,气动力矩系数zm,zzm,等双变元参数的插值,该函数调用了子函数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计算结果运行程序得到弹道曲线,速度曲线,攻角曲线如下:通过运用龙格库塔法和一些插值方法求解导弹运动方程组,获得了导弹各运动参数的变化规律.通过该算法获得的变化规律比较接近真实情况.附件:源程序functionmain()clearallclearglobalglobaly;globalii;ii=0;h=0.005;initl();whiley(7)=0ii=ii+1;result(ii);rk_4(8,h);endsavedata(ii);drawing();%原始数据初始化functioninitl()globalacxacyajzamzafamzwzaxgapamcagcandmandafbLSSONICRHO;globaly;y=[020.1801802052.38];%马赫数ma%攻角alpha%三个时间b(1)为导弹离轨时间b(2)为起飞发动机工作结束时间b(3)为续航发动机工作结束时间b=[02.112644.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.85341.1023;.3858.4086.4903.6290.82261.0723;.3779.4007.4827.6218.81601.0666;.3785.4015.4838.6234.81841.07;.3787.4018.4846.6249.82091.0738;.3829.4062.4897.6310.82481.0835;.3855.4091.4934.6363.83581.0938;.4082.4321.5175.6621.86411.1254;.4947.5192.6073.7571.96721.2392];%升力系数acy=[.0000.64301.47582.28703.07133.8463;.0000.64541.48072.29423.09153.8731;.0000.64801.48582.30143.09153.8731;.0000.65121.49232.31073.10393.8891;.0000.65541.50072.32273.11973.9092;.0000.66171.51342.34093.1