彈道計算大作業(yè)

彈道計算大作業(yè)目錄一、初始條件和要求 )三、仿真結果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