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).*(chi_s*z1+2*chi_s*z1*lamda_s*z(i))+(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).*(chi_s*z3+2*chi_s*z3*lamda_s*(z(i) +h*z2/2))+(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).*(chi_s*z4+2*chi_s*z4*lamda_s*(z(i)+h*z 3))+(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曲线)');
xlabel('\fontsize{12}\bf弹丸行程 l 单位:mm');
ylabel('\fontsize{12}\bf膛内压力 p 单位：Mpa');
%弹丸膛内时期膛底合力 F_pt 单位:N
F_pt=1.34*10^(-4)*p;
%后效期膛底合力 F_pt1 单位:N
t1=0.002288:0.000001:0.006537;
F_pt1=12231.10*exp(-(t1-0.002288)/0.00107);
figure(5)
plot(t*10^3,F_pt,t1*10^3,F_pt1,':','linewidth',1.5);
legend('F_pt','F_pt1');
grid on
title('\fontsize{12}\bf(t-F_pt曲线)');
xlabel('\fontsize{12}\bf时间 t 单位:ms');
ylabel('\fontsize{12}\bf膛底合力 F_pt 单位:N');
tspan=length(t)/30;
tspan=1:ceil(tspan):length(t);
tspan(end)=length(t);
fprintf('t(ms) p(Mpa) v(m/s) l(dm) z psi ');
format short g;
Result=[t(tspan)' (p(tspan)/1e6)' v(tspan)' l(tspan)' z(tspan)' psi(tspan)']。