1、基于 matlab 质点弹道数值计算1. 无控弹质点弹道模型(不考虑姿态运动影响)向量形式* MERGEFORMAT (1.1)()2DxyzSCdVIJKVWgt m其中：、 、 分别为地面坐标系 、 、 轴上的单位向量。地面坐标系是与地球表面固连的坐标系。IJKz坐标原点 定义在炮口断面中心， 轴沿水平线指向射击方向， 轴垂直向上， 轴按右手法则确定。oxyz为速度向量， 为风速向量， 、 、 分别为速度向量在地面坐标系三轴上投影向量的幅值，VWVyz风速向量的分量定义类同。 222()()()xyyzV为空气密度为特征面积，对于轴对称弹取 ， 为弹径。S24dS为弹箭质量m为重力加速度向
2、量，一般计算时取 。g 29.8065/gms为阻力系数(作业中忽略诱导阻力) DC分量形式表示* MERGEFORMAT (1.2)()2DxxSCVW* MERGEFORMAT (1.3)yygm* MERGEFORMAT (1.4)()Dzz质点位置(射程方向) * MERGEFORMAT (1.5)0txXVd(飞行高度) * MERGEFORMAT (1.6)tyY(侧偏) * MERGEFORMAT (1.7)0tzZ2. 仿真模型参数空气密度 ：按照美国标准大气推广委员会(COESA)提出的模型进行计算(随高度变化)，建议采用SIMULINK 中的 COESA Atmospher
3、e Model 模块。当地音速： 按照美国标准大气推广委员会(COESA)提出的模型进行计算(随高度变化)。弹径 105mm质量 14.97kg阻力系数随飞行马赫数的变化如下表：马赫数 0 0.875 0.925 0.965 0.99 1.025 1.085 1.19 1.35 1.8 2 2.5阻力系数 0.124 0.124 0.15 0.2 0.35 0.375 0.415 0.415 .385 0.335 0.318 0.276仿真时采用线性插值的方式得到对应飞行状态的阻力系数，建议采用 SIMULINK 中的 Lookup Table。模块。仿真条件：炮口初速 ，发射点海拔高度为 0
4、，落点海拔为 0。493Vm/s3. 作业内容分别在射角 、 、 的情况下绘制弹丸飞行速度 、飞行高度随时间变化曲线，并求出全飞3056 V行时间、弹道顶点高度、落点侧偏、射程(发射点到落点的水平距离)、落速(落点飞行速度)、落角(落点弹道倾角) 。采用 ODE45 算法求解有关微分方程。4. 作业提交方式交作业时以 WORD 文档形式打印提交，必须包含图形数据(以图形形式)和计算结果(以表格形式)。源程序以文档附件形式提交：应包含 Simulink 仿真程序截图及有关的.m 脚本文件，保证作业结果可以复现。空军标准气象条件 3228.34,1.25/,9.806/onononKkgms()6
5、0yy54.83/(.31)onHy()20.4)Cyy21onRg637Rm子函数:当地声速 C.mfunction C=C(y);T0=288.34;Ty=T0-0.00586.*y;Cy=20.046.*sqrt(Ty);C=Cy;End重力加速度 g.mfunction g=g(y);g0=9.806;R=6371000;gy=g0.*(1-2.*y/R);g=gy;end密度 rou.mfunction rou=rou(y);rou0=1.225;rouy=rou0.*(1-2.0323e-5.*y)4.83;rou=rouy;end解微分方程组主函数 maguo.m functio
6、n d=maguo(t,x);m=14.97;d=0.105;S=pi/4.*d2;Ma=0 0.875 0.925 0.965 0.99 1.025 1.085 1.19 1.35 1.8 2 2.5;Cx0=0.124 0.124 0.15 0.2 0.35 0.375 0.415 0.415 0.385 0.335 0.318 0.276;p1=polyfit(Ma,Cx0,5);d= -(polyval(p1,x(7).*S.*rou(x(5).*x(8).*(x(1)-15)./m./2;-(polyval(p1,x(7).*S.*rou(x(5).*x(8).*(x(2)-0)./
7、m./2-g(x(5);-(polyval(p1,x(7).*S.*rou(x(5).*x(8).*(x(3)-10)./m./2; x(1);x(2);x(3); x(7)-x(8)./C(x(5);x(8)-sqrt(x(1)-15).2+(x(2)-0).2+(x(3)-10).2); %x(9)-sqrt(x(1).2+(x(2).2+(x(3).2); %end调用：clearclcfor k=1:3sita0=30,45,60;t=37.848,54.465,70.81;v0=493;vx0=v0*cos(sita0(k)/180*pi);vy0=v0*sin(sita0(k)/1
8、80*pi);vz0=0;x0=0;y0=0;z0=0;Ma0=v0/C(0); %Cvm0=sqrt(vx0-15).2+(vy0-0).2+(vz0-10).2); %tspan=0,t(k);Xx0= vx0 vy0 vz0 x0 y0 z0 Ma0 vm0 v0; %M=eye(7);M(7,7)=0;M(8,8)=0;M(9,9)=0;options=odeset(Mass,M);t,x=ode15s(maguo,tspan,Xx0,options); %maguo.mfiguresubplot(2,3,1);plot(t,x(:,9);grid on subplot(2,3,2)p
9、lot(t,x(:,5);grid on%subplot(2,3,3); plot(t,x(:,4);grid on subplot(2,3,4);plot(t,x(:,5);grid on%subplot(2,3,5);plot(t,x(:,6);grid onsubplot(2,3,5);plot(t,x(:,8);grid onshegao=max(x(:,5)luosu=x(end,9)shec=sqrt(x(end,4).2+x(end,5).2+x(end,6).2)luojiao=acos(x(end,1)./(sqrt(x(end,1).2)+x(end,2).2)./pi.*180cepian=x(end,6)figureplot3(x(:,6),x(:,4),x(:,5);grid onclearend射角 全飞行时间（S） 弹道顶点高度 落角 射程(KM) 落速 落点侧偏30 37.848 1830.6 89.5353 11.247 313.1140 119.263645 54.465 3610.3 89.7682 14.251 321.4250 141.989760 70.81 5713.0 89.8593 15.426 327.1971 87.6354
