航天飞行动力学大作业(2、3章)

大作业导弹质量0320mkg，转动惯量2315zJkgm，发动机推力2000PN，质量秒消耗量0.46/smkgs，初始速度0250/Vms，初始位置00xm，初始高度07000Hm，初始弹道倾角00，初始俯仰角00，初始攻角00，初始俯仰角速度00/rads，初始速度0250/Vms，参考面积20.45mrefS，参考长度2.5mrefL，升力线斜率0.25/,0.05/zyyCC，阻力系数20.20.03xC，俯仰力矩系数导数0.1/,0.024/zzzmm，阻尼力矩系数导数2zzm，大气密度计算公式：4.25588000001.2495288.150.0065TTTTHT1）飞行方案9100xm*2000cos(0.0003141.1)5000Hx，0.0/smkgs2）飞行方案240009100mxm*3000Hm，0.46/smkgs3）飞行方案24000&&0xmy4）目标位置为30000Txm，采用比例导引法（k=3）。计算：1）纵向理想弹道（采用瞬时平衡假设）。2）飞行器的纵向短周期扰动的动力系数2425342235,,,,aaaaa，并分析其在特征点的自由扰动的稳定性，并计算弹体传递函数(),()ynWsWs（不考虑下洗影响）。二、问题求解⒈采用瞬时平衡假设计算理想弹道根据瞬时平衡假设，得到导弹在铅锤平面内运动的质心运动方程组：cossinsincoscossinzszzzdVmPXmgdtdmVPYmgdtdxVdtdyVdtdmmdtmm①对X,Y有：xyXCqSYCqS②对方案1，有***2000cos(0.0003141.1)5000(-)+(-)zHxkHHkHH&&&③0.0/smkgs④联立①②③④，应用Runge-Kutta法求解，调试得1,0.7kk时，弹道终点误差不超过20m。弹道轨迹如下：速度变化如下：弹道倾角变化如下：对方案2，有**3050(-)+zHmkHHkH⑤0.46/smkgs⑥联立①②⑤⑥，调试得，0.16,0.12kk时，误差在50m以内。弹道轨迹如下：速度变化如下：弹道倾角变化如下：对方案3，目标固定，有**00+()ddqkdtdtkqq⑦00,q可分别由上阶段终止条件和此阶段初始条件求出。sindqrVdt⑧**()()zkk⑨联立①⑦⑧⑨求解，取4,8,6kkk时落点误差在10m以内。弹道轨迹如下：速度变化如下：弹道倾角变化如下：三阶段整体轨迹如下：弹道轨迹飞行速度变化弹道倾角变化2、不考虑气动力下洗影响，计算飞行器的纵向短周期扰动的动力系数，并分析其在各时间点的自由扰动的稳定性，并计算在各个特征秒弹体传递函数(),()ynWsWs⑴0xm,有：02200240025003400350()134.2167()6.7108()-1.6106()0.0514()()0.0053()zzzzzzzzzMaJMaJMaJPYamVYamV由于24223413.60950aaa，所以纵向自由扰动稳定。253424352422340.0087aaaaKaaa2512534243513.6158aTaaaa24223410.2710Taaa223424223418.19572aaaaa1222(1)0.0087(13.61581)()(21)(0.07349.86211)KTssWssTsTssss02220.2219()(21)0.07349.86211ynVKWsgTsTsss⑵9100xm,有：02200240025003400350()97.4253()4.8713()1.1691()0.0644()()0.0056()zzzzzzzzzMaJMaJMaJPYamVYamV由于24223411.14550aaa，所以纵向自由扰动稳定253424352422342512534243524223422342422340.009211.395410.299514.59832aaaaKaaaaTaaaaTaaaaaaaa12220222(1)0.0092(11.39541)()(21)(0.08978.74441)0.1609()(21)0.08978.74441ynKTssWssTsTssssVKWsgTsTsss⑶24000xm,有：02200240025003400350()109.3729()5.4686()-1.3125()0.0807()()0.0075()zzzzzzzzzMaJMaJMaJPYamVYamV由于24223414.29500aaa，所以纵向自由扰动稳定。25342435242234251253424352422342234242234-0.01038.937410.264414.47192aaaaKaaaaTaaaaTaaaaaaaa12220222(1)0.0103(8.93741)()(21)(0.06997.65271)0.1913()(21)0.06997.65271ynKTssWssTsTssssVKWsgTsTsss⑷30000xm,有：02200240025003400350()263.5062()13.1753()-3.1621()0.1057()()0.0143()zzzzzzzzzMaJMaJMaJPYamVYamV由于24223441.02790aaa，所以纵向自由扰动稳定。25342435242234251253424352422342234242234-0.01276.055010.156120.57702aaaaKaaaaTaaaaTaaaaaaaa12220222(1)0.0127(6.05501)()(21)(0.02446.42411)0.3149()(21)0.02446.42411ynKTssWssTsTssssVKWsgTsTsss三、代码1、main.m(第一问主函数，进行方程组求解及画图，可根据需要对注释部分进行取消)clearallop=odeset('Events',@eventHeight);[T1,Y1,T1end,Y1end,evennum]=ode45(@Height,[0,inf],[250,0,7000,0],op);t1=T1end;t1init1=Y1end;Y1endh=2000*cos(0.000314*1.1*Y1end(:,2))+5000;op=odeset('Events',@eventHeightTwo);[T2,Y2,T2end,Y2end,evennum]=ode45(@HeightTwo,[t1,inf],init1,op);Y2endinit2=Y2end;T2endt2=T2end;op=odeset('Events',@eventHeightThree);[T3,Y3,T3end,Y3end,evennum]=ode45(@HeightThree,[t2,inf],init2,op);Y3endT3endfigure(1)plot(Y1(:,2),Y1(:,3),'color','b');holdon%plot(Y1(:,2),h,':','color','r');%legend('实际弹道','理想弹道');%xlabel('x'),ylabel('H/H*'),title('飞行轨迹');%plot(Y1(:,2),Y1(:,1),'color','b');%xlabel('x'),ylabel('V'),title('飞行速度');%plot(Y1(:,2),Y1(:,4),'color','b');%xlabel('x'),ylabel('弹道倾角'),title('弹道倾角');%h=2000*cos(0.000314*1.1*Y1end(:,2))+5000;%hd=-2000*0.000314*1.1*sin(0.000314*1.1*Y1end(:,2))*Y1end(:,1)*cos(Y1end(:,4));%delta=-1*(Y1end(:,3)-h)-0.7*(Y1end(:,1)*sin(Y1end(:,4))-hd);%delta%plot(Y1(:,2),delta,'color','b');%delta=-0.16*(Y2end(:,3)-3050)-0.12*Y2end(:,1)*sin(Y2end(:,4));%deltaplot(Y2(:,2),Y2(:,3),'color','b');holdon%plot(Y2(:,2),3050,':','color','r');%xlim([910024000]);%legend('实际弹道','理想弹道');%xlabel('x'),ylabel('H/H*'),title('弹道轨迹');plot(Y3(:,2),Y3(:,3),'color','b');xlabel('x'),ylabel('H'),title('弹道轨迹');figure(2)plot(Y1(:,2),Y1(:,1));holdonplot(Y2(:,2),Y2(:,1));holdonplot(Y3(:,2),Y3(:,1));xlabel('x'),ylabel('V'),title('飞行速度');figure(3)plot(Y1(:,2),Y1(:,4));holdonplot(Y2(:,2),Y2(:,4));holdonplot(Y3(:,2),Y3(:,4));xlabel('x'),ylabel('弹道倾角'),title('弹道倾角');%xlabel('x'),ylabel('H'),title('飞行轨迹')%plot(Y2(:,2),Y2(:,1),'color','b');%xlabel('x'),ylabel('V'),title('飞行速度');%plot(Y2(:,2),Y2(:,4),'color','b');%xlabel('x'),ylabel('弹道倾角'),title('弹道倾角');%plot(Y3(:,2),Y3(:,1),'color','b');%xlabel('x'),ylabel('V'),title('飞行速度');%plot(Y3(:,2),Y3(:,3),'color','b');%xlabel('x'),ylabel('H'),title('弹道轨迹');%plot(Y3(:,2),Y3(:,4),'color','b');%xlabel('x'),ylabel('弹道倾角'),title('弹道倾角')2、main2.m(第二问函数，使用时需键入h,v,m,)Jz=315;h=;rou=1.2495*(1-0.0065*h/288.15)^4.25588;v=;m=;a22=2/Jz;a24=0.1*0.5*rou*v^2*0.45*2.5/Jz;a25=-0.024*0.5*rou*v^2*0.45*2.5