弹道模型matlab代码

弹道模型matlab代码

matlab内弹道程序

1、火炮内弹道求解与计算摘要：本文结合火炮内弹道基本方程，得出压力、速度与行程、时间的关系式。并利用了MATLAB的程序对该火炮系统的内弹道过程进行求解。关键词：内弹道基本方程；MATLAB；1.火炮内弹道诸元火炮内弹道诸元数据如下表所示：炮膛断面积S药室容积V0弹丸全行程Ig弹丸质量m装药质量dm2dm3dmkgkg0.8187.9247.4815.65.5火药参数如下表所示：F燃气比

热比k管状火药长2a管状火药厚

kJ/kgdm3/kgkg/dm31mmmm96011.61.22601.7协调常量如下表所示：BIk挤进压力P011kPasMPa1.6021.2761601.930其他所需的参数计。

2、算：；2.内弹道基本方程组及其解析解法方程组建立如上，则考虑三个时

期分别求解：前期：考虑为定容燃烧过程，则有条件：则有，令第一时期：将前期的参量计算得出之后，代入方程组，解算第一时期的v、p值。考虑平均法，利用若设x=Z-Z0则可得，第二时期：考虑第二时期无火药燃烧，则有：设极限速度，利用可得各个时期的p-l，v-l曲线。3.使用MATLAB对内弹道进行求解由于解析

解方法较为繁琐，并且需要相当多的简化才能进行计算，因此考虑使用MATLAB

对内弹道方程进行求解与仿真，描绘p-t、p-l、v-t、v-l曲线，如下图所示。最大膛压约为800MPa，出膛速度大约为1000m/s.4.Matlab代码代。

3、码：functionndd%100mm加农炮S=0.818;%枪(炮)膛横断面积

dm2M=15.6;%弹重kgV0=7.92;%药室容积dm3I_g=47.48;%身管行程

dmP_0=30000;%起动压力kpafai1=1.02;%次要功系数theta=0.2;%火药热力系数%=f=;%火药力kg*dm/kgalpha=1;%余容dm3/kgdelta=1.6;%火药密度

kg/dm3%=ome=5.5;%装药量kgu1=1.6184*10-5;%第一种装药烧速系数

dm3/(s*kg)n1=1;%装药压力指数n1lambda=-0.5;%。

4、装药形状特征量lambda_s=0;%装药分裂点形状特征量schi=2.01;%装药形状特征量chi_s=0;%装药分裂点形状特征量smu=0;%装药形状特征量

et1=1.7*10-2;%装药药厚0d1=1.7*10-2;%装药火药内径dB=1.602;%=%常数

与初值计算-

l_0=V0/S;Delta=ome/V0;phi=1.276;v_j=196*f*ome/(phi*theta*M);v_j=sqrt(

v_j);Z_s=1;p_0=P_0/(f*Delta);psi_0=(1/Delta-1/delta)/(f/P_0+alpha-

1/delta);Z_0=(s。

5、qrt(1+4*psi_0*lambda/chi)-1)/(2*lambda);%解算子-

C=zeros(1,12);C(1)=chi;C(2)=lambda;C(3)=lambda_s;C(4)=chi_s;C(5)=Z_s;% C(6)=theta;C(7)=B;C(8)=n1;C(9)=Delta;C(10)=delta;C(11)=alpha;C(12)=mu; C;y0=Z_0;0;0;psi_0;options=odeset(outputfcn,odeplot);tt,y=ode45(ndd_fun, 0:100,Z_0;0;0,options,C);l=y(:,2)。

6、;l=l*l_0;fl=find(l=I_g);fl=min(fl);tt,y=ode45(ndd_fun,0:0.005:fl,Z_0;0;0 ,options,C);Z=y(:,1);lx=y(:,2);vx=y(:,3);psi=(Z=0&Z=1&Z=Z_s)*1;l_psi=1-(Delta/delta)*(1-psi)-alpha*Delta*psi;px=(psi-

vx.*vx)https:///tags/(lx+l_psi);p=px*f*Delta/100;v=vx*v_j/10;l =lx*l_0;t=tt*l_0*10。

7、

00/v_j;fl=find(l=I_g);fl=min(fl)+1;p(fl:end)=;v(fl:end)=;l(fl:end)=;t(fl:end)=;p d=px*f*Delta/100/(1+ome/3/fai1/M);pt=pd*(1+ome/2/fai1/M);aa=max(px );M=find(px=aa);Pm=tt(M)*l_0*1000/v_jlx(M)*l_0vx(M)*v_j/10px(M)*f*Delta /100pt(M)pd(M)psi(M)Z(M);%ll=length(tt);ran=find(Z=1);ran=min(ran);Zf=t t。

8、

(ran)*l_0*1000/v_jlx(ran)*l_0vx(ran)*v_j/10px(ran)*f*Delta/100pt(ran)pd(ran) psi(ran)Z(ran);jie=find(psi=1);jie=min(jie);psij=tt(jie)*l_0*1000/v_jlx(jie)*l_0v x(jie)*v_j/10px(jie)*f*Delta/100pt(jie)pd(jie)psi(jie)Z(jie);pg=tt(end)*l_0*1000 /v_jlx(end)*l_0vx(end)*v_j/10px(end)*f*Delta/1。

9、

00pt(end)pd(end)psi(end)Z(end);Ry1=Zf;psij;pg;Pm;Ry2=tt*l_0*1000/v_jlx*l_ 0vx*v_j/10px*f*Delta/100ptpdpsiZ;subplot(2,2,1);plot(t,p,linewidth,2);grido n;xlabel(fontsize8bft(ms);ylabel(fontsize8bfp(kg/cm2);title(fontsize8bft-p曲线);subplot(2,2,2)plot(t,v,linewidth,2);gridon;xlabel(fontsize8。

10、bft(ms);ylabel(fontsize8bfv(m/s);title(fontsize8bft-v曲

线);subplot(2,2,3)plot(l,p,linewidth,2);gridon;xlabel(fontsize8bfl(dm);ylabel(f ontsize8bfp(kg/cm2);title(fontsize8bfl-p曲

线);subplot(2,2,4)plot(l,v,linewidth,2);gridon;xlabel(fontsize8bfl(dm);ylabel(f ontsize8bfv(m/s);title(fontsize8bfl-v曲线);。