Contents

clear;
clc;

this script will take inputs from design decisions and CEA and will return an M_dot

givens

max pressure: ~650psi from stand+expected loss Thrust: we decide Mm from combustion reaction and chem math

% cea run gives:
% gamma
% chamber temp
% R bar

%{
Sutton Equation:
V_ex = sqrt(((2*gammat*R*T0)/(gammat-1))*(1-(14.7/500)^((gammat-1)/gammat)));

Thrust Equation:
T=M_dot*V_ex
%}

%solving for molecular mass of products
%{
reaction:
9N2O + C2H6O = 3CO + 9N2 + 4H2O
%}
Wg = 3*(12.011+15.999)+9*(2*(14.007))+4*(2*(1.0078)+15.999); %weight in grams of product
Mm = Wg/(3+9+4); % grams/ mole of product
Mm = Mm/1000; % to get grams/mole to kilograms/mole



% CEA results
gammat=1.2714;
T0=2196.6;
Mm_CEA=20.195;
Mm_CEA=Mm_CEA/1000;
R_bar_CEA=8.3144598/Mm_CEA;

% calculated deet
R_bar_calc=8.3144598/Mm;         %universal gas constant made specific by deviding by specific molar mass
P_chamber=.75*(1000-100);   %[psi]
P_exit=14.7;                %[psi]

V_ex = sqrt(((2*gammat*R_bar_CEA*T0)/(gammat-1))*(1-(P_exit/P_chamber)^((gammat-1)/gammat)));
% M_dot=Thrust/V_ex;

%{
ok, at this point i understand how which inputs to take from what

next is to create a script that will take in:
    Mm,Gamma,T_chamber
for a range of o/f ratio, then output V_exit to calc m_dot for given thrust

to plot m_dot vs o/f ratio for a given chamber pressure & thrust
%}

graphic function

%inputs from CEA RUN @ 650psi 2-3 o/f ratio in increments of 0.1

Gamma_Vect = [1.2917 1.2871 1.2826 1.2784 1.2743 1.2704 1.2666 1.2628 1.2592 1.2556 1.2521];
M_bar_vect = [18.531 18.882 19.224 19.557 19.880 20.195 20.502 20.800 21.090 21.373 21.648];
T_chamber = [1892.30 1957.80 2020.30 2080.10 2137.30 2192.00 2244.30 2294.50 2342.60 2388.60 2432.80];


P_chamber=650;
R_bar_CEA=8.3144598./(M_bar_vect/1000);


%make this a for loop to iterate through each item in list duh


for i = 1:length(Gamma_Vect)
    V_ex_vect(i)=exitVelocityCalc(Gamma_Vect(i),R_bar_CEA(i),T_chamber(i),P_exit,P_chamber);
end

fprintf('Exit Velocity [m/s]\n')
fprintf('%10.2f', V_ex_vect)
fprintf('\n')

of_vect=linspace(2,3,11);

%Thrust=m_dot*v_ex
thrust=1000; %[lbf]
thrust=thrust*4.44822; %[N]
m_dot_vect=thrust./V_ex_vect;

fprintf('Mass Flow Rate [kg/s]\n')
fprintf('%10.2f', m_dot_vect)
fprintf('\n')

yyaxis left
plot(of_vect, V_ex_vect, 'b-', 'LineWidth', 1.5);
ylabel('Exit Velocity [m/s]')

yyaxis right
plot(of_vect, m_dot_vect, 'r-', 'LineWidth', 1.5);
ylabel('Mass flow Rate [kg/s]')

xlabel('O/F Ratio')
title(sprintf('O/F Ratio vs Exit Velocity - given P_c=%d PSI and T=%d lbf', P_chamber, thrust/4.44822))
grid on

Exit Velocity [m/s]
   2079.35   2100.01   2118.91   2136.13   2152.07   2166.65   2180.10   2192.78   2204.47   2215.39   2225.64
Mass Flow Rate [kg/s]
      2.14      2.12      2.10      2.08      2.07      2.05      2.04      2.03      2.02      2.01      2.00

Published with MATLAB® R2024b