火炮内弹道求解与计算

火炮内弹道求解与计算文档编制序号：[KKIDT-LLE0828-LLETD298-POI08]

火炮内弹道求解与计算

摘要：本文结合火炮内弹道基本方程，得出压力、速度与行程、时间的关系式。并利用了MATLAB 的程序对该火炮系统的内弹道过程进行求解。 关键词：内弹道基本方程；MATLAB ；

1.火炮内弹道诸元

火炮内弹道诸元数据如下表所示：

炮膛断面积S

药室容积V 0

弹丸全行程I g

弹丸质量m

装药质量ω

dm 2 dm 3 dm kg kg

火药参数如下表所示：

F

燃气比热比k

管状火药长2a

管状火药厚δ2

kJ/kg dm 3/kg kg/dm 3 1 mm mm 960

1

260

协调常量如下表所示：

B

Ik 挤进压力P0

1 1 kPa ·s MPa

30

其他所需的参数计算：

1b

==

δα；30

1054.6a

-⨯==

δβ；

01.21=++=βαχ；50.01-

-=++++=β

ααβ

βαλ；

2.内弹道基本方程组及其解析解法

方程组建立如上，则考虑三个时期分别求解：

①前期：考虑为定容燃烧过程，则有条件：MPa p p V V v x 30,0,0,000======

则有025.011

V 0

00

0=-+-

=

ραρω

ψp f ，013.021

4100=-+=

λ

ψχλ

Z

令99.04100

=+=ψχ

λ

σ ②第一时期：将前期的参量计算得出之后，代入方程组，解算第一时期的v 、p 值。

考虑ψV 平均法，利用2

0ψ

ψψψV V V V +==

若设x=Z-Z 0

则可得x x m

SI v k 3.658==ϕ，ψψθψωθψωl l x

B S f V V x B f p +-=+-

=2

222

③第二时期：考虑第二时期无火药燃烧，则有： 设极限速度66.162812=-=

m

k f v j ϕω

）（

)1()(122

111j k

k k j v v l l l l v v -++-=-，l l v v S f P j +-⋅=12

2

1ω

利用①~③可得各个时期的p-l ，v-l 曲线。

3.使用MATLAB 对内弹道进行求解

由于解析解方法较为繁琐，并且需要相当多的简化才能进行计算，因此考虑使用MATLAB 对内弹道方程进行求解与仿真，描绘p-t 、p-l 、v-t 、v-l 曲线，如下图所示。最大膛压约为800MPa ，出膛速度大约为1000m/s.

代码

代码： function ndd

%100mm加农炮

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

M=; %弹重 kg

V0=; %药室容积 dm^3

I_g=; %身管行程 dm

P_0 =30000; %起动压力 kpa

fai1=; %次要功系数

theta =; %火药热力系数

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

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

delta=; %火药密度ρ kg/dm^3

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

ome=; %装药量 kg

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

n1=1; %装药压力指数n1

lambda=; %装药形状特征量λ

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

chi=; %装药形状特征量χ

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

mu=0; %装药形状特征量μ

et1=*10^-2; %装药药厚δ0

d1=*10^-2; %装药火药内径d

B=;

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

%常数与初值计算-----------------------------------------------------------------

l_0=V0/S;

Delta=ome/V0;

phi=;

v_j=196*f*ome/(phi*theta*M);

v_j=sqrt(v_j);

Z_s=1;

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

fl = min(fl);

[tt,y] = ode45(@ndd_fun,0::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>=I_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/M);

pt=pd*(1+ome/2/fai1/M);

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