内弹道及枪膛合力Matlab程序clear;close all; format longd=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*z2*(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*z3*(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; endend figure(1)plot(t*10^3,p/1e6,'linewidth',1.5);hold ongrid ontitle('\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 ontitle('\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 ontitle('\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 ontitle('\fontsize{12}\bf(l-p曲线)');xlabel('\fontsize{12}\bf弹丸行程 l 单位:mm');ylabel('\fontsize{12}\bf膛内压力 p 单位：Mpa'); %弹丸膛内时期膛底合力 F_pt 单位:NF_pt=1.34*10^(-4)*p;%后效期膛底合力 F_pt1 单位:Nt1=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 ontitle('\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)']。
