%BurnSim V0.0.1
%固体火箭发动机燃烧过程仿真器

clear all;
%参数名                     单位      含义
d_nozzle_begin     = 6      % mm       喷嘴初始直径
nozzle_burn_speed  = 0      % mm/s     喷嘴烧蚀速率
d_hole(1)        =  6       % mm       药柱内孔径
D_outer        = 24             % mm       药柱外径
N             = 2              %          药柱能够燃烧的横断面数
L(1)        = 70               % mm       药柱可燃部分长度
prop         = 'KNSU'               %          推进剂名称

%初始值                                             单位
pi = 3.1415926;                                 %圆周率  
d_nozzle(1) = d_nozzle_begin;                   %   mm
r(1) = burnsim_propellant( 0.1 ,prop,'burnspeed');  %   mm/s
V_gram(1) =(D_outer^2-d_hole(1)^2)*pi/4*L(1)/1000; %   cm^3
c = burnsim_propellant(0,prop,'vc');            %   m/s
de = burnsim_propellant(0,prop,'density');      %   g/cm^3
t=0;                                            % s
dt          = 0.001                % s        分析时间片长度,数值小精度高,数值大需要内存少,修改后注意要更改图像的X坐标名称

i = 1;%循环计数
while( V_gram(i) > 0 )
%计算公式                                                       得数含义              单位
Ae(i) = d_nozzle(i)^2 * pi / 4;                                %  当前喷管直径        mm^2
S_burn(i) =( D_outer^2 - d_hole(i)^2 )*pi/4*N + d_hole(i)*pi*L(i); %  当前燃烧面积        mm^2
Kn(i) = S_burn(i) / Ae(i) ;                                    %  喷燃比
P(i) = Kn(i) * r(i) * c * de                                  %  压力                MPa
r(i+1) = burnsim_propellant(P(i),prop,'burnspeed');            %  当前燃速            mm/s
d_hole(i+1) = d_hole(i) + 2*r(i)*dt;                           %  下一时刻的药柱内孔径 mm
L(i+1) = L(i) - N * dt * r(i);                                 %  下一时刻的药柱长度  mm
V_gram(i+1) = (D_outer^2 - d_hole(i)^2)* pi /4 *L(i+1) /1000;    %  下一时刻的药柱体积  cm^3
d_nozzle(i+1) = d_nozzle(i) + 2*nozzle_burn_speed*dt;           %下一时刻的喷管内径 mm

    disp('T = ');      disp(t);
    disp('Vgram = ');  disp(V_gram(i));
    disp('[s:9] = ');      disp(P(i));
    disp('Kn = ');     disp(Kn(i));
    
i = i+1;
t = t + dt;
    
end

公式的单位有问题？我计算出来压力比自己估计的大10倍左右。