火炮内弹道求解与计算 This model paper was revised by the Standardization Office on December 10, 2020
火炮内弹道求解与计算
摘要：本文结合火炮内弹道基本方程，得出压力、速度与行程、时间的关系式。并利用了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()(12
2
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)]; %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); grid on; xlabel('\fontsize{8}\bft (ms)'); 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)'); 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})'); 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曲线'); 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; Result = [t(tspan) p(tspan) v(tspan) l(tspan)] format; %-------------------------------------------------------------------------- function dy = ndd_fun(t,y,C) chi=C(1);lambda=C(2);lambda_s=C(3);chi_s=C(4);Z_s=C(5);mu=C(12); theta=C(6);B=C(7);V=C(8);Delta=C(9);delta=C(10);alpha=C(11); Z = y(1); l = y(2); v = 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; p = ( psi - v*v )/( l + l_psi ); dy(1) = sqrt(theta/(2*B))*(p^V)*(Z>=0&Z<=Z_s); dy(2) = v; dy(3) = theta*p/2; dy = [dy(1);dy(2);dy(3)]; 参考文献 [1]郭新鹏,赵军利.基于MATLAB的枪炮内弹道程序设计及仿真[J].高校理科研究. [2]吴晶,刘金元.局域MATLAB的舰炮内弹道计算模块的GUI设计[J].舰船电子工 程.2014,34,6:94-98. Wu Jing,Liu design of calculation module of the interior ballistic for ship gun based on matlab[J].Ship electronic Engineering,2014,34,6:94-98. [3]钱林方.火炮弹道学[M].2009:116-183.
相关主题
文本预览