弹道计算大作业.pdf

《弹道计算大作业.pdf》由会员分享，可在线阅读，更多相关《弹道计算大作业.pdf（14页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、弹道计算大作业目录一、初始条件和要求.21.1 初始条件.21.2 仿真要求.2二、模型的建立.22.1 升力和阻力模型.22.2 大气和重力加速度模型.32.3 无控飞行.32.4 平衡滑翔.42.5 最大升阻比滑翔飞行弹道.4三、仿真结果.53.1 无控飞行弹道仿真.53.2 平衡滑翔弹道仿真.73.3 最大升阻比滑翔弹道仿真.8附录.9一、初始条件和要求1.1 初始条件已知给定的初始条件如下：表 1 初始条件名称参考面积展弦比效率因子质量重力加速度零升阻力系数密度符号给定值1.70.860.9115通过模型计算可得0.02通过模型计算可得单位m2/kgm/s2/radkg/m3SemgC

2、x01.2 仿真要求请使用 Simulink 或 Buildfly 完成以下仿真任务：（1）请完成该导弹的无控飞行弹道仿真；（2）请完成该导弹的平衡滑翔方案飞行弹道仿真；（3）请完成该导弹的最大升阻比滑翔飞行弹道仿真；二、模型的建立2.1 升力和阻力模型已知展弦比为的飞行器的升力线斜率为：Cy1 122(1)根据飞行力学相关知识，飞行器的升力系数和阻力系数为：Cy Cy2CxCx0Cy(2)其中，升力线斜率由（1）式可得；为效率系数：由升力系数和阻力系数，得到导弹的升力和阻力为：12X Cv Sx2Y C1v2Sy21。e(3)2.2 大气和重力加速度模型在计算过程中，大气密度采用如下模型：0

3、.0065H=01T04.25588(4)其中，01.225kg/m3为海平面的大气密度；T0 288.15K。重力加速度采用如下模型：Rdg g0Rd H2(5)其中，g0 9.8，Rd 6371000m为地球半径；H为飞行器距离地面的高度。2.3 无控飞行假设导弹的运动始终在铅垂平面，根据飞行力学知识，得到导弹无控飞行时的运动学和动力学方程为：XdV gsindtmdYgcosdtmVVdxV cosdtdydtV sin(6)在上述模型中，假设俯仰角为 0。2.4 平衡滑翔d 0。因此得到dt所谓的“平衡”可以理解为垂直于速度方向受力平衡，即平衡滑翔时的导弹运动学和动力学方程：XdV g

4、sindtmdYgcosdtmVVdxV cosdtdyV sindtd 0dtY12Cyv S2(7)由于弹道倾角的变化率为常数，方程组中的第二个方程等于 0。这个方程可以用来求攻角。2.5 最大升阻比滑翔飞行弹道联立（1）式、（2）式可得升阻比的表达式为：21 12Cy2Cx2Cx01 122(8)从上式可以看出，由于展弦比、零升阻力系数Cx0为常数，因此升阻比只和攻角有关，是关于攻角的函数。因此要使升阻比达到最大，须使dCyCxd 0得到2Cx0e1 12Cx0eCx0aaCyCy因此，以最大升阻比滑翔时导弹运动学和动力学方程为：dV X gsindtmdYgcosdtmVVdxdtV

5、cosdydtV sinCx0Cay三、仿真结果3.1 无控飞行弹道仿真根据无控弹道模型，写出 s 函数，搭建的仿真模块如下图所示：(9)(10)图 1 无控飞行仿真模块由于初始条件给定，因此模块没有输入；输出有六个，分别为导弹的射程变化、高度变化、速度变化、弹道倾角变化、攻角变化以及密度变化。模块的仿真时间由高度变化决定，当高度降为 0（导弹落到地面上）时仿真结束。导出数据后画图如下：飞 行 器 空 间 运 行 轨 迹20001800951600140090100速 度 随 时 间 变 化 图 像垂直位移y1200100080060040085速度v02000400060008000水 平

6、位 移 x100001200014000807570200065020406080100时 间 t120140160180弹 道 倾 角 随 时 间 变 化 图 像-0.08-0.09-0.1-0.110.180.170.160.150.14攻 角 随 时 间 变 化 图 像弹道倾角-0.12攻角020406080100时 间 t120140160180-0.13-0.14-0.15-0.16-0.17-0.180.130.120.110.10.090.08020406080100时 间 t120140160180图 2 无控飞行时各参数变化3.2 平衡滑翔弹道仿真平衡滑翔弹道仿真模块如下图所

7、示：图 3 平衡滑翔模块取仿真时间为 150s，无输入，输出分别为：导弹的射程变化、高度变化、速度变化、弹道倾角变化、攻角变化以及密度变化。得到各参量时间变化图如下：飞 行 器 空 间 运 行 轨 迹2000180016001400200220速 度 随 时 间 变 化 图 像180垂直位移y12001000800600400200000.51水 平 位 移 x弹 道 倾 角 随 时 间 变 化 图 像10.110.10.50.090.081.522.5x 10速度v160140120100050时 间 t100150攻 角 随 时 间 变 化 图 像弹道倾角00.07攻角-0.5-1-1.5

8、050时 间 t1001500.060.050.040.030.020.01050时 间 t100150图 4 平衡滑翔飞行时各参数变化3.3 最大升阻比滑翔弹道仿真按最大升阻比飞行时弹道仿真模块如下图所示：图 5 最大升阻比飞行模块取仿真时间为 180s，无输入，输出分别为：导弹的射程变化、高度变化、速度变化、弹道倾角变化、攻角变化以及密度变化。得到各参量时间变化图如下：飞 行 器 空 间 运 行 轨 迹250010095200090150085速 度 随 时 间 变 化 图 像垂直位移y1000速度v02000400060008000水 平 位 移 x1000012000140008075

9、50070065-50060020406080100时 间 t120140160180弹 道 倾 角 随 时 间 变 化 图 像0.31.5攻 角 随 时 间 变 化 图 像0.210.1弹道倾角-0.1攻角0-0.5020406080100时 间 t120140160180-1000.5-0.2-0.3-0.420406080100时 间 t120140160180图 4 最大升阻比飞行时各参数变化附录附表 1 无控弹道飞行时完整的s 函数无控弹道function sys,x0,str,ts,simStateCompliance=trace2(t,x,u,flag)switch flag,c

10、ase 0,sys,x0,str,ts,simStateCompliance=mdlInitializeSizes;case 1,sys=mdlDerivatives(t,x,u);case 2,sys=mdlUpdate(t,x,u);case 3,sys=mdlOutputs(t,x,u);case 4,sys=mdlGetTimeOfNextVarHit(t,x,u);case 9,sys=mdlTerminate(t,x,u);otherwise DAStudio.error(Simulink:blocks:unhandledFlag,num2str(flag);endfunction

11、 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 =0 0;simStateCompliance=UnknownSimState;function

12、sys=mdlDerivatives(t,x,u)S=1.7;%参考面积,m2 AR=0.86;%展弦比 e=0.9;%效率因子;m=115;%质量,kg g0=9.8;%海平面重力加速度,m/s2 Rd=6371000;%地球半径 r=Rd/(Rd+x(2);g=g0*r2;%飞行器所在高度的重力加速度 rho0=1.225;%海平面大气密度，kg/m3 T0=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);%升力线斜率,/

13、rad CDo=0.02;%零升阻力系数 epsilon=1/(pi*e*AR);%诱导阻力因子 CL=CLa*alpha;%升力系数 CD=CDo+epsilon*CL2;%阻力系数 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;function sys=mdlUpdate(t,x,u)sys=;functio

14、n sys=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;function sys=mdlGetTimeOfNextVarHit(t,x,u)sampleTime=1;%Example,set the next hit to be one second later.sys=t+sampleTime;function sys=mdlTerminate(t,x,u)sys=;附表 2 平衡滑翔飞行部分代码平

15、衡滑翔飞行function sys=mdlDerivatives(t,x,u)S=1.7;%参考面积,m2 AR=0.86;%展弦比 e=0.9;%效率因子;m=115;%质量,kg g0=9.8;%海平面重力加速度,m/s2 Rd=6371000;%地球半径 r=Rd/(Rd+x(2);g=g0*r2;%飞行器所在高度的重力加速度 rho0=1.225;%海平面大气密度，kg/m3 T0=288.15;rho=rho0*(1-0.0065*x(2)/T0)4.2288;%飞行器所在高度的大气密度 CLa=3.141592*AR/(1+sqrt(1+(AR/2)2);%升力线斜率,/rad C

16、Do=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*CL2;%阻力系数 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;function sys=mdlOutputs(t,x,u)rho0=1.

17、225;T0=288.15;rho=rho0*(1-0.0065*x(2)/T0)4.2288;S=1.7;%参考面积,m2AR=0.86;%展弦比m=115;%质量,kgg0=9.8;%海平面重力加速度,m/s2Rd=6371000;%地球半径r=Rd/(Rd+x(2);g=g0*r2;%飞行器所在高度的重力加速度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;sy

18、s=y;附表 3 最大升阻比飞行部分代码最大升阻比飞行function sys=mdlDerivatives(t,x,u)S=1.7;%参考面积,m2 AR=0.86;%展弦比 e=0.9;%效率因子;m=115;%质量,kg g0=9.8;%海平面重力加速度,m/s2 Rd=6371000;%地球半径 r=Rd/(Rd+x(2);g=g0*r2;%飞行器所在高度的重力加速度 rho0=1.225;%海平面大气密度，kg/m3 T0=288.15;rho=rho0*(1-0.0065*x(2)/T0)4.2288;%飞行器所在高度的大气密度 CLa=3.141592*AR/(1+sqrt(1+

19、(AR/2)2);%升力线斜率,/rad CDo=0.02;%零升阻力系数 epsilon=1/(pi*e*AR);%诱导阻力因子 alpha=sqrt(CDo*pi*AR*e)/CLa;CL=CLa*alpha;%升力系数 CD=CDo+epsilon*CL2;%阻力系数 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;function sys=mdlOutputs(t,x,u)rho0=1.225;T0=288.15;rho=rho0*(1-0.0065*x(2)/T0)4.2288;AR=0.86;%展弦比 e=0.9;%效率因子;CLa=3.141592*AR/(1+sqrt(1+(AR/2)2);%升力线斜率,/rad CDo=0.02;%零升阻力系数 alpha=sqrt(CDo*pi*AR*e)/CLa;sys=x(1)x(2)x(3)x(4)alpha rho;