冲压发动机流推力计算程序matlab

psi=7;%循环静温比
f=0.0675;%燃料空气当量比
fh_pr=2818000;%f*hpr[J/kg]
h_pr=fh_pr/f;%反应热[J/kg]
T0=222;%参考温度K
R=289.3;%热力学常数
C_pc=1090;%压缩过程定压比热容[J/（kg*K）]
C_pb=1510;%热附加过程定压比热容[J/（kg*K）]
gamma_c=1.362;%压缩过程比热比
gamma_b=1.238;%热附加过程比热比
a0=sqrt（1.4*R*T0）;%入口声速[m/s]
F=zeros（413）;%未装机比推力[N*s/kg]
eta_0=zeros（413）;%总效率
eta_c=zeros（413）;%绝热压缩效率
eta_e=zeros（413）;%绝热膨胀效率
I_sp=zeros（413）;%比冲[s]
for i=1:41
    for j=1:3
    Ma0=4+0.2*（i-1）;    
    V0=Ma0*a0;
    Sa0=V0*（1+R*T0/V0^2）;
    T3=psi*T0;
    V3=sqrt（V0^2-2*C_pc*T0*（psi-1））;
    %Sa3=V3*（1+R*T3/V3^2）;
    c=0.005+j*0.005;%取值0.01、0.015、0.02
    psi1=1+（（gamma_c-1）/2）*Ma0^2*（1-（1+c）^2*（1+（2-2*psi）/（（gamma_c-1）*Ma0^2）））;
    eta_c_=（1-0.05*j）*（（1-1/psi1）/（1-1/psi））;%eta1取值0.95、0.9、0.85
    eta_c（ij）=eta_c_;
    ratio_P_3_0=（psi/（psi*（1-eta_c_）+eta_c_））^（C_pc/R）;%压缩部件压强比
    ratio_A_3_0=psi*（V3/V0）/ratio_P_3_0;%压缩部件面积比
    b=0.05+0.05*j;%取值0.1、0.15、0.2
    V4=V3*（（1+f*0.5）/（1+f）-b/（1+f））;