MATLAB内弹道程序 - 毕设专用!!!

内弹道及枪膛合力Matlab程序

clear;

close all;

format long

d=0.0127;

S=0.82*0.0127^2;

V0=2.04e-5;

l_0=V0/S;

lg=0.924;

f=1000000;

alpha=0.001;

w=0.017;

rou=1600;

theta=0.2;

phi=1.45;

chi=0.79825;

lamda=0.1387;

mu=-0.043956;

e1=0.00052/2;

u1=7.5991e-10;

Is=e1/u1;

chi_s=1.2645;

lamda_s=-0.31322;

zk=1.4434;

%Ik=447000;

m=0.048;

p0=30e6;

delta=800;

psi0=(1/delta-1/rou)/(f/p0+alpha-1/rou); sigma0=sqrt(1+4*lamda*psi0/chi);

z0=2*psi0/chi/(sigma0+1);%(sigma0-1)/2/lamda; %====赋予初值====%

v(1)=0;

l(1)=0;

p(1)=p0;

z(1)=z0;

psi(1)=psi0;

lpsi(1)=l_0*(1-delta/rou-(alpha-1/rou)*delta*psi(1));

t(1)=0;

h=0.000001;

for i=1:100000

z1=p(i)/Is;

v1=S*p(i)/m/phi;

l1=v(i);

psi1=(z(i)>=0&z(i)<=1).*(chi+2*chi*z(i)*lamda+3*chi*mu*z(i)^2)*z1 +(z(i)>1&z(i)zk).*0;

lpsi1=-l_0*(alpha-1/rou)*delta*psi1;

p1=((f*w/S+p(i)*l_0*delta*(alpha-1/rou))*psi1-theta*phi*m*v1*v(i)/S-p(i )*l1)/(l(i)+lpsi(i));

z2=(p(i)+h*p1/2)/Is;

v2=S*(p(i)+h*p1/2)/m/phi;

l2=v(i)+h*v1/2;

psi2=(z(i)>=0&z(i)<=1).*(chi*z2+2*chi*(z(i)+h*z1/2)*lamda*z2+3*chi*mu*z 2*(z(i)+h*z1/2)^2)+(z(i)>1&z(i)<=zk).*(chi_s*z2+2*chi_s*z2*lamda_s*(z(i )+h*z1/2))+(z(i)>zk).*0;

lpsi2=-l_0*(alpha-1/rou)*delta*psi2;

p2=((f*w/S+(p(i)+h*p1/2)*l_0*delta*(alpha-1/rou))*psi2-theta*phi*m*v2*( v(i)+h*v1/2)/S-(p(i)+h*p1/2)*l2)/((l(i)+h*l1/2)+(lpsi(i)+h*lpsi1/2));

z3=(p(i)+h*p2/2)/Is;

v3=S*(p(i)+h*p2/2)/m/phi;

l3=v(i)+h*v2/2;

psi3=(z(i)>=0&z(i)<=1).*(chi*z3+2*chi*(z(i)+h*z2/2)*lamda*z3+3*chi*mu*z 3*(z(i)+h*z2/2)^2)+(z(i)>1&z(i)zk).*0;

lpsi3=-l_0*(alpha-1/rou)*delta*psi3;

p3=((f*w/S+(p(i)+h*p2/2)*l_0*delta*(alpha-1/rou))*psi3-theta*phi*m*v3*( v(i)+h*v2/2)/S-(p(i)+h*p2/2)*l3)/((l(i)+h*l2/2)+(lpsi(i)+h*lpsi2/2));

z4=(p(i)+h*p3)/Is;

l4=v(i)+h*v3;

v4=S*(p(i)+h*p3)/m/phi;

psi4=(z(i)>=0&z(i)<=1).*(chi*z4+2*chi*(z(i)+h*z3)*lamda*z4+3*chi*mu*z4* (z(i)+h*z3)^2)+(z(i)>1&z(i)=zk).*0;

lpsi4=-l_0*(alpha-1/rou)*delta*psi4;

p4=((f*w/S+(p(i)+h*p3)*l_0*delta*(alpha-1/rou))*psi1-theta*phi*m*v1*v(i )/S-(p(i)+h*p3)*l4)/((l(i)+h*l3)+(lpsi(i)+h*lpsi3));

z(i+1)=z(i)+h*(z1+2*z2+2*z3+z4)/6;

l(i+1)=l(i)+h*(l1+2*l2+2*l3+l4)/6;

v(i+1)=v(i)+h*(v1+2*v2+2*v3+v4)/6;

psi(i+1)=psi(i)+h*(psi1+2*psi2+2*psi3+psi4)/6;

lpsi(i+1)=lpsi(i)+h*(lpsi1+2*lpsi2+2*lpsi3+lpsi4)/6;

p(i+1)=p(i)+h*(p1+2*p2+2*p3+p4)/6;

t(i+1)=t(i)+h;

if l(i)>=lg;

T=i;

break;

end

end

figure(1)

plot(t*10^3,p/1e6,'linewidth',1.5);

hold on

grid on

title('\fontsize{12}\bf(t-p曲线)');

xlabel('\fontsize{12}\bf时间 t 单位:ms');

ylabel('\fontsize{12}\bf膛内压力 p 单位:Mpa');

figure(2)

plot(t*1e3,v,'linewidth',1.5);

grid on

title('\fontsize{12}\bf(v-t曲线)');

xlabel('\fontsize{12}bf时间 t 单位:ms');

ylabel('\fontsize{12}\bf弹丸速度 v 单位:m/s');

figure(3)

plot(l*1e1,v,'linewidth',1.5);

grid on

title('\fontsize{12}\bf(l-v曲线)');

xlabel('\fontsize{12}\bf弹丸行程 l 单位:mm');

ylabel('\fontsize{12}\bf弹丸速度 v 单位:m/s');

figure(4)

plot(l*10,p/1e6,'linewidth',1.5);

grid on

title('\fontsize{12}\bf(l-p曲线)');