彈道計(jì)算大作業(yè)

彈道計(jì)算大作業(yè)目錄一、初始條件和要求 )三、仿真結(jié)果3.1無(wú)控飛行彈道仿真根據(jù)無(wú)控彈道模型，寫(xiě)出s函數(shù)，搭建的仿真模塊如下圖所示：圖1無(wú)控飛行仿真模塊由于初始條件給定，因此模塊沒(méi)有輸入；輸出有六個(gè)，分別為導(dǎo)彈的射程變化、高度變化、速度變化、彈道傾角變化、攻角變化以及密度變化。模塊的仿真時(shí)間由高度變化決定，當(dāng)高度降為0（導(dǎo)彈落到地面上）時(shí)仿真結(jié)束。導(dǎo)出數(shù)據(jù)后畫(huà)圖如下： 圖2無(wú)控飛行時(shí)各參數(shù)變化3.2平衡滑翔彈道仿真平衡滑翔彈道仿真模塊如下圖所示：圖3平衡滑翔模塊取仿真時(shí)間為150s，無(wú)輸入，輸出分別為：導(dǎo)彈的射程變化、高度變化、速度變化、彈道傾角變化、攻角變化以及密度變化。得到各參量時(shí)間變化圖如下：圖4平衡滑翔飛行時(shí)各參數(shù)變化3.3最大升阻比滑翔彈道仿真按最大升阻比飛行時(shí)彈道仿真模塊如下圖所示：圖5最大升阻比飛行模塊取仿真時(shí)間為180s，無(wú)輸入，輸出分別為：導(dǎo)彈的射程變化、高度變化、速度變化、彈道傾角變化、攻角變化以及密度變化。得到各參量時(shí)間變化圖如下：圖4最大升阻比飛行時(shí)各參數(shù)變化附錄附表1無(wú)控彈道飛行時(shí)完整的s函數(shù)無(wú)控彈道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;%質(zhì)量,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);%無(wú)控飛行時(shí)CLa=3.141592*AR/(1+sqrt(1+(AR/2)^2));%升力線斜率,/radCDo=0.02;%零升阻力系數(shù)epsilon=1/(pi*e*AR);%誘導(dǎo)阻力因子CL=CLa*alpha;%升力系數(shù)CD=CDo+epsilon*CL^2;%阻力系數(shù)X=CD*1/2*rho*x(3)^2*S;Y=CL*1/2*rho*x(3)^2*S;%以下為飛行器在鉛垂平面的運(yùn)動(dòng)方程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;%質(zhì)量,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;%零升阻力系數(shù)epsilon=1/(pi*e*AR);%誘導(dǎo)阻力因子alpha=2*m*g*cos(x(4))/(rho*x(3)^2*S*CLa);CL=CLa*alpha;%升力系數(shù)CD=CDo+epsilon*CL^2;%阻力系數(shù)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;%質(zhì)量,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;%質(zhì)量,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;%零升阻力系數(shù)epsilon=1/(pi*e*AR);%誘導(dǎo)阻力因子alpha=sqrt(CDo*pi*AR*e)/CLa;CL=CLa*alpha;%升力系數(shù)CD=CDo+epsilon*CL^2;%阻力系數(shù)X=CD*1/2*rho*x(3)^2*S;Y=CL*1/2*rho*x(3)^2*S;%以下為飛行器在鉛垂平面的運(yùn)動(dòng)方程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