内弹道计算

clc; clear;
%%%%%%%%%火炮构造及弹丸诸元
S=2.66e-3;               %枪（炮）膛横断面积S  m^2
V0=1.51e-3;               %药室容积V0  m^3
l_g=3.624;             %身管行程lg  m
m=2.8;                 %弹重m  kg
%%%%%%%%装药条件
f=950e3;                %火药力f  J/kg 
omega=1.16;             %装药量ω  kg
alpha=0.001;            %余容α  m^3/kg
rou_p=1600;             %ρp   kg/m^3
theta =0.25;            %火药热力系数
u1=5.127e-8;            %燃速系数  dm^3/（s*kg）
n=0.83;                 %燃速指数 n
e1=5.5e-4;              %装药药厚e1 dm
d=5.5e-4;               % 
chi=0.75;               %装药形状特征量 χ
lambda=0.12;            %装药形状特征量 λ
mu=0.0;                 %装药形状特征量 μ
chi_s=1.696;            %装药分裂点形状特征量 χs
lambda_s = -0.4104;     %装药分裂点形状特征量 λs
%%%%%%%%%%起始条件
p_0 =3e7;                    %起动压力  Pa
%%%%%%%%%%计算常数       用于计算φ，本程序中已提供该参数，所以以下两个参数未用到
phi1=1.02;                   
lambda2=0;
%K=1.03;                      %运动阻力系数K
%%%%%%%%%%%%常量计算
%phi=phi1+lambda2*omega/m;
%phi=K+omega/3/m;                  
phi=1.168;                        %次要功系数φ
Delta=omega/V0;                   %装填密度
l_0 = V0/S;                       %
vj=sqrt（2*f*omega/（theta*phi*m））; %
B = S^2*e1^2*（f*Delta）^（2-2*n）/（f*omega*phi*m*u1^2）;
lg1=l_g/l_0;
psi_s=chi*（1+lambda+mu）;
%es=（1+chi_s-sqrt（（1+chi_s）^2-4*chi_s*psi_s））/（2*chi_s）;
%rou=e1/es-e1;
%rou=0.2956*（d/2+e1）;
%Zk=（e1+rou）/e1;
Zk=（chi_s+sqrt（chi_s^2+4*psi_s*（1-chi_s）））/（2*psi_s）;  %%%%需确认的部分
rou=Zk*e1-e1;                                          %%%%%%%%%%应先求ρ，进而确定Zk
%%%%%%%%%%%%%%%%初值计算
psi_0=（1/Delta-1/rou_p）/（f/p_0