内弹道计算程序

Ro1=0; %药型系数α1
%=========================================
%常数与初值计算-----------------------------------------------------------------
l_0=W0/A;
内弹道计算程序.txt始终相信，这世间，相爱的原因有很多，但分开的理由只有一个--爱的还不够。人生有四个存折：健康 情感 事业和金钱。如果健康消失了，其他的存折都会过期。 %59nian130
A=0.87; %枪(炮)膛横断面积A dm^2
G=19;%33.4; %弹重 kg
aa=max(px);
M=find(px==aa);
Pm=[tt(M)*l_0*1000/v_j lx(M)*l_0 vx(M)*v_j/10 px(M)*f*Delta/100 pt(M) pd(M) psi(M) Z(M)];
%ll=length(tt);
ran=find(Z>=1);
options = odeset('outputfcn','odeplot');
[tt,y] = ode45(@ndd_fun,0:100,[Z_0;0;0],options,C);
l = y(:,2);
l = l*l_0;
fl = find(l>=l_g);
fl = min(fl);
[tt,y] = ode45(@ndd_fun,0:0.005:fl,[Z_0;0;0],options,C);
W0=2.04; %药室容积 dm^3
l_g=25.0; %身管行程 dm
P_0 =30000; %起动压力 kpa
fai1=1.02; %次要功系数
K=1.03; %运动阻力系数φ1
theta =0.2; %火药热力系数
Result = [t(tspan) p(tspan) v(tspan) l(tspan)]
format;
ylabel('\fontsize{8}\bfp (kg/cm^{2})');
title('\fontsize{8}\bft-p曲线');
subplot(2,2,2)
plot(t,v,'linewidth',2);
grid on;
xlabel('\fontsize{8}\bft (ms)');
psij=[tt(jie)*l_0*1000/v_j lx(jie)*l_0 vx(jie)*v_j/10 px(jie)*f*Delta/100 pt(jie) pd(jie) psi(jie) Z(jie)];
pg=[tt(end)*l_0*1000/v_j lx(end)*l_0 vx(end)*v_j/10 px(end)*f*Delta/100 pt(end) pd(end) psi(end) Z(end)];
chi=1.00716; %第一种装药形状特征量χ1
chi_s=0; %第一种装药分裂点形状特征量χ1s
mu=0; %第一种装药形状特征量μ1
et1=1.14*10^-2; %第一种装药药厚δ01
d1=2.5*10^-2; %第一种装药火药内径d1
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];
%=========================================
f=950000; %火药力 kg*dm/kg
alpha=1; %余容 dm
%==================================
ome=2.2;%12.9; %第一种装药量 kg
u1=5.0024*10^-5; %第一种装药烧速系数 dm^3/(s*kg)
n1=0.82; %第一种装药的压力指数n1
lambda=-0.0071; %第一种装药形状特征量λ1
lambda_s=0; %第一种装药分裂点形状特征量λ1s
title('\fontsize{8}\bfl-p曲线');
subplot(2,2,4)
plot(l,v,'linewidth',2);
grid on;
xlabel('\fontsize{8}\bfl (dm)');
ylabel('\fontsize{8}\bfv (m/s)');
title('\fontsize{8}\bfl-v曲线');
Z = y(:,1);lx = y(:,2); vx = y(:,3);
psi = (Z>=0&Z<1).*( chi*Z.*(1 + lambda*Z + mu*Z) ) +...%%%%%%%%%
(Z>=1&Z<Z_s).*( chi_s*Z.*(1 + lambda_s*Z) ) +...
(Z>=Z_s)*1;
l_psi = 1 - (Delta/delta)*(1-psi) - alpha*Delta*psi;
px = ( psi - vx.*vx )./( lx + l_psi );
p = px*f*Delta/100;
v = vx*v_j/10;
l = lx*l_0;
tspan = length(t)/20;
tspan = 1:ceil(tspan):length(t);
tspan(end) = length(t);
fprintf(' t(ms) p(kg/cm^2) v(m/s) l(dm)');
format short g;
Ry1=[Zf;psij;pg;Pm];
Ry2=[tt*l_0*1000/v_j lx*l_0 vx*v_j/10 px*f*Delta/100 pt pd psi Z];
subplot(2,2,1);
plot(t,p,'linewidth',2);
grid on;
xlabel('\fontsize{8}\bft (ms)');
ran=min(ran);
Zf=[tt(ran)*l_0*1000/v_j lx(ran)*l_0 vx(ran)*v_j/10 px(ran)*f*Delta/100 pt(ran) pd(ran) psi(ran) Z(ran)];
jie=find(psi>=1);
jie=min(jie);
p_0=P_0/(f*Delta);
psi_0=(1/Delta - 1/delta)/(f/P_0 + alpha - 1/delta);
Z_0=(sqrt(1+4*psi_0*lambda/chi) - 1)/(2*lambda);
%解算子-----------------------------------------------------------------------
Delta=ome/W0;
phi=K + ome/(3*G);
v_j=196*f*ome/(phi*theta*G);
v_j=sqrt(v_j);
B = 98*(et1*A)^2/( u1*u1*f*ome*phi*G );
B=B*(f*Delta)^(2-2*n1);
Z_s=1+Ro1*(d1/2+et1)/et1;
t = tt*l_0*1000/v_j;
fl = find(l>=l_g);
fl = min(fl)+1;
p(fl:end)=[];v(fl:end)=[];l(fl:end)=[];t(fl:end)=[];
pd=px*f*Delta/100/(1+ome/3/fai1/G);
pt=pd*(1+ome/2/fai1/G);
ylabel('\fontsize{8}\bfv (m/s)');
title('\fontsize{8}\bft-v曲线');
subplot(2,2,3)
plot(l,p,'linewidth',2);
grid on;
xlabel('\fontsize{8}\bfl (dm)');
ylabel('\fontsize{8}\bfp (kg/cm^{2})');