CU Spaceflight Wiki

I successfully (at least I believe it works) made simple MatLab script which calculates key conical nozzle parameters with input of:

• Combustion chamber temperature and pressure
• Gas properties (Molecular mass and cp/cv coefficient)
• Ambient pressure
• Half angle of nozzle
• Combustion products flow rate.

All these values are calculated under following assumptions:

• Working gas is homogeneous
• All combustion products are gaseous
• Working gas is a perfect gas
• Flow is adiabatic
• Friction and boundary layer effects are neglected
• There are no shock waves or discontinuities in the nozzle flow
• Flow is steady.
• All gases at the exhaust have axially directed velocity
• Gas pressure, velocity, temperature and density are uniform across any section normal to the nozzle axis.
• Flow is frozen i.e. there is no change in gas composition during nozzle expansion

(All theory is taken from a lovely book “Rocket Propulsion Elements” by G.P. Sutton and O. Biblarz)

Here is a code:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Nozzle and gas properties definitions (You need to fill this section, User!)%%
T1 = 2000; % Combustion chamber temperature [K]
p1 = 2 * 10^6; % Combustion chamber pressure [Pa]
k = 1.3; % cp/cv of the gas 
miu = 16; % Molar mass of gas [kg/kmol]
p2 = 10^5; % Ambient pressure [Pa]
mdot = 0.1028; % Mass flow rate [kg/s]
alpha = 15 * pi/180; % half angle of the conical nozzle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Key preformance parameters equations (Just lots of equations) %%

R = 8314.3 / miu;
disp ('Supersonic throat condition')
condition = p1/p2 - ((k+1)/2)^(k/(k-1)) 
disp ('All calculations are valid only if above number is positive') 

disp ('Throat area') 
At = mdot * (p1 * k * sqrt( (2/(k+1)) ^ ((k+1)/(k-1)) / (k * R * T1)))^(-1)

disp (' Throat temperature [K]') 
Tt = 2* T1 /(k+1)

disp('Ratio of exit to throat area')
epsilon = (((k+1)/2)^(1/(k-1)) * (p2/p1)^(1/k) * sqrt( (k+1)/(k-1) * (1-(p2/p1)^((k-1)/k))) )^(-1)

disp ('Exit area [m^2]') 
A2= epsilon * At % Exit area [m^2]

disp ('Nozzle lenght [m]')
L = cot(alpha) * sqrt(At/pi) * (sqrt(epsilon)-1) % Nozzle lenght [m]

disp ('Exit velocity [m/s]')
v2 = sqrt ( 2*k /(k-1) * R * T1 * (1 - (p2/p1)^((k-1)/k) ) ) 

disp ('Total Theoretical Thrust [N]')
F = mdot * v2 

disp ('Theoretical Specific Impulse') 
Is = v2/9.81
%% End of the script %% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%