内弹道计算程序

%内弹道计算程序
function ndd
clc;clear;
%构造与装填条件------------------------------------
S =1.681; %枪(炮)膛横断面积 dm^2
q = 10.24; %弹重 kg
W_0 = 10.35; %药室容积 dm^3
l_g = 62.25; %身管行程 dm
%进程的进行条件-------------------------------------
P_0 = 45000; %起动压力 kpa
theta = 0.25; %火药热力系数
K = 1.222; %次要功系数
%装填条件--------------------------------------------
delta =1.6; %火药密度 kg/dm^3
f = 948000; %火药力 kg*dm/kg
alpha =1 %余容 dm^3/kg
omega = 10.35; %药量 kg
V = 0.9627; %烧速指数
u =0.0000088; %烧速系数 dm^3/(s*kg)
%火药特征(仅适用于多孔火药)-----------------------------
e_1 = 0.0088; %厚度(肉厚)/2 dm
d_0 = 0.022; %孔道直径 dm
D_0 = 0.1364; %药粒直径(梅花形不用管此项) dm
c = 1.3 %药粒长度/2 dm
n = 7; %孔数
flag = 1; %1:圆柱形多孔火药 2:梅花形七孔火药 3:梅花形十四孔火药 4:梅花形十九孔火药
%形状函数 ------
Flag = [ 1 n 0 D_0 0 0.2956;
2 8 65.2968 d_0+4*e_1 d_0+2*e_1 0.1547;
8/3 47/3 141.4764 d_0+4*e_1 d_0+2*e_1 0.1547;
3 21 195.8903 d_0+4*e_1 d_0+2*e_1 0.1547];
F1 = Flag(flag,1);
F2 = Flag(flag,2);
F3 = Flag(flag,3);
F4 = Flag(flag,4);
F5 = Flag(flag,5);
F6 = Flag(flag,6);
Pi_1 = ( F1*F4 + F2*d_0 )/(2*c);
Q_1 = ( F3*F5*F5 + F1*F4*F4 - F2*d_0^2 )/( 4*c^2 );
beta = e_1/c;
chi = ( Q_1 + 2*Pi_1 )*beta/Q_1;
lambda = ( n - 1 - 2*Pi_1 )*beta/( Q_1 + 2*Pi_1 );
mu = ( 1 - n )*beta^2/( Q_1 + 2*Pi_1 );
psi_s = chi*( 1 + lambda + mu );
rho = F6*( d_0/2 + e_1 );
Z_s = 1 + rho/e_1;
chi_s = ( psi_s*Z_s^2 - 1 )/( Z_s^2 - Z_s );
lambda_s = psi_s/chi_s - 1;
%常数与初值计算-------------------------------------------------------
l_0 = W_0/S;
Delta = omega/W_0;
phi = K + omega/(3*q);
v_j = 196*f*omega/(phi*theta*q);
v_j = sqrt(v_j);
B = 98*(e_1*S)^2/( u*u*f*omega*phi*q );
B = B*(f*Delta)^(2-2*V);
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,11);
C(1)=chi;C(2)=lambda;C(3)=lambda_s;C(4)=chi_s;C(5)=Z_s;C(12)=mu;
C(6)=theta;C(7)=B;C(8)=V;C(9)=Delta;C(10)=delta;C(11)=alpha;
%解算子-------------------------------------------------------------------
options = odeset('outputfcn','odeplot');
[t,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);
[t,y] = ode45(@ndd_fun,0:0.005:fl,[Z_0;0;0],options,C);
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 );
p = p*f*Delta/98.0665;
v = v*v_j/10;
l = l*l_0;
t = t*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)=[];
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}\bfp-t曲线');
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}\bfv-t曲线');
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}\bfp-l曲线');
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}\bfv-l曲线');
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;