内弹道计算

59-130加农炮内弹道计算

function ndd

%59-130

A=1.394; %枪(炮)膛横断面积A dm^2

G=33.4; %弹重kg

W0=18.56; %药室容积dm^3

l_g=59.52; %身管行程dm

P_0 =30000; %起动压力kpa

fai1=1.02; %次要功系数

K=1.03; %运动阻力系数φ1

theta =0.2; %火药热力系数

%=========================================

f=950000; %火药力kg*dm/kg

alpha=1; %余容dm^3/kg

delta=1.6; %火药重度γ

%==================================

ome=12.9; %第一种装药量kg

u1=5.0024*10^-5; %第一种装药烧速系数dm^3/(s*kg)

n1=0.82; %第一种装药的压力指数n1

lambda=-0.0071; %第一种装药形状特征量λ1

lambda_s=0; %第一种装药分裂点形状特征量λ1s

chi=1.00716; %第一种装药形状特征量χ1

chi_s=0; %第一种装药分裂点形状特征量χ1s

mu=0; %第一种装药形状特征量μ1

et1=1.14*10^-2; %第一种装药药厚δ01

d1=2.5*10^-2; %第一种装药火药内径d1

Ro1=0; %药型系数α1

%=========================================

%常数与初值计算----------------------------------------------------------------- l_0=W0/A;

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;

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);

%解算子-----------------------------------------------------------------------

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);

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);

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>=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;

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);

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);

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);

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)];

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);