火炮内弹道求解与计算
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
火炮内弹道求解与计算
摘要:本文结合火炮内弹道基本方程,得出压力、速度与行程、时间的关系式。并利用了MA TLAB 的程序对该火炮系统的内弹道过程进行求解。
关键词:内弹道基本方程;MA TLAB ;
1.火炮内弹道诸元
火炮内弹道诸元数据如下表所示:
炮膛断面积S
药室容积V 0
弹丸全行程I g
弹丸质量m
装药质量ω
dm 2 dm 3 dm kg kg 0.818
7.92
47.48
15.6
5.5
火药参数如下表所示:
F α
p ρ
燃气比热比k
管状火药长2a
管状火药厚δ2
kJ/kg dm 3/kg kg/dm 3 1 mm mm 960
1
1.6
1.2
260
1.7
协调常量如下表所示:
B ϕ
Ik 挤进压力P0
1 1 kPa ·s MPa 1.602
1.276
1601.9
30
其他所需的参数计算:
1b
==
δα;30
1054.6a
-⨯==
δβ;
01.21=++=βαχ;50.01--=++++=β
ααβ
βαλ;
2.内弹道基本方程组及其解析解法
⎪⎪⎪⎪⎪⎪
⎩⎪⎪
⎪⎪⎪
⎪⎨
⎧
=-=---+===+=v dt
dl mv f V V p Sp dl dv v Sp dt dv
m I p dt dZ Z Z p k 2
02))1((m )1(ϕθωψαωψψρωϕϕλχψ及 方程组建立如上,则考虑三个时期分别求解:
①前期:考虑为定容燃烧过程,则有条件: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.
4.Matlab代码
代码:
function ndd
%100mm加农炮
S=0.818; %枪(炮)膛横断面积dm^2
M=15.6; %弹重kg
V0=7.92; %药室容积dm^3
I_g=47.48; %身管行程dm
P_0 =30000; %起动压力kpa
fai1=1.02; %次要功系数
theta =0.2; %火药热力系数
%=========================================
f=960000; %火药力kg*dm/kg
alpha=1; %余容dm^3/kg
delta=1.6; %火药密度ρkg/dm^3
%==================================
ome=5.5; %装药量kg
u1=1.6184*10^-5; %第一种装药烧速系数dm^3/(s*kg)
n1=1; %装药压力指数n1
lambda=-0.5; %装药形状特征量λ
lambda_s=0; %装药分裂点形状特征量λs
chi=2.01; %装药形状特征量χ
chi_s=0; %装药分裂点形状特征量χs
mu=0; %装药形状特征量μ
et1=1.7*10^-2; %装药药厚δ0
d1=1.7*10^-2; %装药火药内径d
B=1.602;
%=========================================
%常数与初值计算-----------------------------------------------------------------
l_0=V0/S;
Delta=ome/V0;
phi=1.276;
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: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;