弹道计算大作业

弹道计算大作业目录一、初始条件和要求 )三、仿真结果3.1无控飞行弹道仿真根据无控弹道模型，写出s函数，搭建的仿真模块如下图所示：图1无控飞行仿真模块由于初始条件给定，因此模块没有输入；输出有六个，分别为导弹的射程变化、高度变化、速度变化、弹道倾角变化、攻角变化以及密度变化。模块的仿真时间由高度变化决定，当高度降为0（导弹落到地面上）时仿真结束。导出数据后画图如下： 图2无控飞行时各参数变化3.2平衡滑翔弹道仿真平衡滑翔弹道仿真模块如下图所示：图3平衡滑翔模块取仿真时间为150s，无输入，输出分别为：导弹的射程变化、高度变化、速度变化、弹道倾角变化、攻角变化以及密度变化。得到各参量时间变化图如下：图4平衡滑翔飞行时各参数变化3.3最大升阻比滑翔弹道仿真按最大升阻比飞行时弹道仿真模块如下图所示：图5最大升阻比飞行模块取仿真时间为180s，无输入，输出分别为：导弹的射程变化、高度变化、速度变化、弹道倾角变化、攻角变化以及密度变化。得到各参量时间变化图如下：图4最大升阻比飞行时各参数变化附录附表1无控弹道飞行时完整的s函数无控弹道function[sys,x0,str,ts,simStateCompliance]=trace2(t,x,u,flag)switchflag,case0,[sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes;case1,sys=mdlDerivatives(t,x,u);case2,sys=mdlUpdate(t,x,u);case3,sys=mdlOutputs(t,x,u);case4,sys=mdlGetTimeOfNextVarHit(t,x,u);case9,sys=mdlTerminate(t,x,u);otherwiseDAStudio.error('Simulink:blocks:unhandledFlag',num2str(flag));endfunction[sys,x0,str,ts,simStateCompliance]=mdlInitializeSizessizes=simsizes;sizes.NumContStates=4;sizes.NumDiscStates=0;sizes.NumOutputs=5;sizes.NumInputs=0;sizes.DirFeedthrough=0;sizes.NumSampleTimes=1;sys=simsizes(sizes);x0=[0;2000;100;-5/180*pi];str=[];ts=[00];simStateCompliance='UnknownSimState';functionsys=mdlDerivatives(t,x,u)S=1.7;%参考面积,m^2AR=0.86;%展弦比e=0.9;%效率因子;m=115;%质量,kgg0=9.8;%海平面重力加速度,m/s^2Rd=6371000;%地球半径r=Rd/(Rd+x(2));g=g0*r^2;%飞行器所在高度的重力加速度rho0=1.225;%海平面大气密度，kg/m^3T0=288.15;rho=rho0*(1-0.0065*x(2)/T0)^4.2288;%飞行器所在高度的大气密度alpha=-x(4);%无控飞行时CLa=3.141592*AR/(1+sqrt(1+(AR/2)^2));%升力线斜率,/radCDo=0.02;%零升阻力系数epsilon=1/(pi*e*AR);%诱导阻力因子CL=CLa*alpha;%升力系数CD=CDo+epsilon*CL^2;%阻力系数X=CD*1/2*rho*x(3)^2*S;Y=CL*1/2*rho*x(3)^2*S;%以下为飞行器在铅垂平面的运动方程dx=x(3)*cos(x(4));dy=x(3)*sin(x(4));dv=-X/m-g*sin(x(4));dtheta=Y/(m*x(3))-g*cos(x(4))/x(3);sys=[dx;dy;dv;dtheta];functionsys=mdlUpdate(t,x,u)sys=[];functionsys=mdlOutputs(t,x,u)y1=x(1);y2=x(2);y3=x(3);rho0=1.225;T0=288.15;rho=rho0*(1-0.0065*x(2)/T0)^4.2288;sys=[x(1)x(2)x(3)x(4)rho];functionsys=mdlGetTimeOfNextVarHit(t,x,u)sampleTime=1;%Example,setthenexthittobeonesecondlater.sys=t+sampleTime;functionsys=mdlTerminate(t,x,u)sys=[];附表2平衡滑翔飞行部分代码平衡滑翔飞行functionsys=mdlDerivatives(t,x,u)S=1.7;%参考面积,m^2AR=0.86;%展弦比e=0.9;%效率因子;m=115;%质量,kgg0=9.8;%海平面重力加速度,m/s^2Rd=6371000;%地球半径r=Rd/(Rd+x(2));g=g0*r^2;%飞行器所在高度的重力加速度rho0=1.225;%海平面大气密度，kg/m^3T0=288.15;rho=rho0*(1-0.0065*x(2)/T0)^4.2288;%飞行器所在高度的大气密度CLa=3.141592*AR/(1+sqrt(1+(AR/2)^2));%升力线斜率,/radCDo=0.02;%零升阻力系数epsilon=1/(pi*e*AR);%诱导阻力因子alpha=2*m*g*cos(x(4))/(rho*x(3)^2*S*CLa);CL=CLa*alpha;%升力系数CD=CDo+epsilon*CL^2;%阻力系数X=CD*1/2*rho*sqrt(x(3)^2)*S;Y=CL*1/2*rho*sqrt(x(3)^2)*S;dx=x(3)*cos(x(4));dy=x(3)*sin(x(4));dv=-X/m-g*sin(x(4));dtheta=0;sys=[dx;dy;dv;dtheta];functionsys=mdlOutputs(t,x,u)rho0=1.225;T0=288.15;rho=rho0*(1-0.0065*x(2)/T0)^4.2288;S=1.7;%参考面积,m^2AR=0.86;%展弦比m=115;%质量,kgg0=9.8;%海平面重力加速度,m/s^2Rd=6371000;%地球半径r=Rd/(Rd+x(2));g=g0*r^2;%飞行器所在高度的重力加速度CLa=pi*AR/(1+sqrt(1+(AR/2)^2));%升力线斜率,/radalpha=2*m*g*cos(x(4))/(rho*x(3)^2*S*CLa);y(1)=x(1);y(2)=x(2);y(3)=x(3);y(4)=x(4);y(5)=alpha;y(6)=rho;sys=[y];附表3最大升阻比飞行部分代码最大升阻比飞行functionsys=mdlDerivatives(t,x,u)S=1.7;%参考面积,m^2AR=0.86;%展弦比e=0.9;%效率因子;m=115;%质量,kgg0=9.8;%海平面重力加速度,m/s^2Rd=6371000;%地球半径r=Rd/(Rd+x(2)); g=g0*r^2;%飞行器所在高度的重力加速度rho0=1.225;%海平面大气密度，kg/m^3T0=288.15;rho=rho0*(1-0.0065*x(2)/T0)^4.2288;%飞行器所在高度的大气密度CLa=3.141592*AR/(1+sqrt(1+(AR/2)^2));%升力线斜率,/radCDo=0.02;%零升阻力系数epsilon=1/(pi*e*AR);%诱导阻力因子alpha=sqrt(CDo*pi*AR*e)/CLa;CL=CLa*alpha;%升力系数CD=CDo+epsilon*CL^2;%阻力系数X=CD*1/2*rho*x(3)^2*S;Y=CL*1/2*rho*x(3)^2*S;%以下为飞行器在铅垂平面的运动方程dx=x(3)*cos(x(4));dy=x(3)*sin(x(4));dv=-X/m-g*sin(x(4));dtheta=Y/(m*x(3))-g*cos(x(4))/x(3);sys=[dx;dy;dv;dtheta];functionsys=mdlOutputs(t,x,u)rho0=1.225;T0=288.15