火炮内弹道求解与计算

火炮内弹道求解与计算摘要：本文结合火炮内弹道基本方程，得出压力、速度与行程、时间的关系式。

并利用了MATLAB 的程序对该火炮系统的内弹道过程进行求解。

关键词：内弹道基本方程；MATLAB ；1.火炮内弹道诸元火炮内弹道诸元数据如下表所示：炮膛断面积S药室容积V 0弹丸全行程I g弹丸质量m装药质量ωdm 2 dm 3 dm kg kg 0.8187.9247.4815.65.5火药参数如下表所示：F燃气比热比k管状火药长2a管状火药厚δ2kJ/kg dm 3/kgkg/dm 3 1 mm mm 96011.61.22601.7协调常量如下表所示：BIk 挤进压力P01 1 kPa ·s MPa 1.6021.2761601.930其他所需的参数计算：1b==δα；301054.6a-⨯==δβ；01.21=++=βαχ；50.01--=++++=βααββαλ；2.内弹道基本方程组及其解析解法方程组建立如上，则考虑三个时期分别求解：①前期：考虑为定容燃烧过程，则有条件：MPa p p V V v x 30,0,0,000======则有025.011V 0000=-+-=ραρωψp f ，013.0214100=-+=λψχλZ令99.04100=+=ψχλσ ②第一时期：将前期的参量计算得出之后，代入方程组，解算第一时期的v 、p 值。

考虑ψV 平均法，利用20ψψψψV V V V +==若设x=Z-Z 0则可得x x mSI v k 3.658==ϕ，ψψθψωθψωl l xB S f V V x B f p +-=+-=2222③第二时期：考虑第二时期无火药燃烧，则有： 设极限速度66.162812=-=mk f v j ϕω）（)1()(122111j k k k j v vl l l l v v -++-=-，l l v v S f P j +-⋅=1221ω 利用①~③可得各个时期的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^2M=15.6; %弹重kgV0=7.92; %药室容积dm^3I_g=47.48; %身管行程dmP_0 =30000; %起动压力kpafai1=1.02; %次要功系数theta =0.2; %火药热力系数%=========================================f=960000; %火药力kg*dm/kgalpha=1; %余容dm^3/kgdelta=1.6; %火药密度ρkg/dm^3%==================================ome=5.5; %装药量kgu1=1.6184*10^-5; %第一种装药烧速系数dm^3/(s*kg) n1=1; %装药压力指数n1lambda=-0.5; %装药形状特征量λlambda_s=0; %装药分裂点形状特征量λschi=2.01; %装药形状特征量χchi_s=0; %装药分裂点形状特征量χsmu=0; %装药形状特征量μet1=1.7*10^-2; %装药药厚δ0d1=1.7*10^-2; %装药火药内径dB=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_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;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_s).*( chi_s*Z.*(1 + lambda_s*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 Jinyuan.GUI 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.。