RunDesignHybrid.m

clc
clear all
close all

%% Code for Hybrid System Design
% Based on G. Zilliac's slides on "Hybrid Propulsion System Design" 
% Not yet validated on real designs. 
% This code is for PMMA - need to upgrade for propellant selection and
% built in CEA. 

%% Conversion factors

PA_TO_PSI = 0.000145038;
M_TO_IN = 39.3701;
KG_TO_LBS = 2.20462;

%% INPUTS

% chamber and atmospheric pressure
p_c = 350 / PA_TO_PSI; % Chamber pressure [Pa]
p_o = 101325; % Atmospheric pressure [Pa]

% results from CEA obtained by maximizing C* using PMMA
of_ratio = 4.5; % Optimal O/F ratio 
% ISP = 257.13; % [s]
gamma = 1.131659; % gamma in the chamber
c_star = 1661.77; % [m/s]
ep = 2.426; % nozzle area ratio for full expansion at 1 atm

ep = 4.54;

% Regression Coefficients 
a = 0.00146; 
n = 0.35; 

% efficiencies
zeta_d = 1.07; % Discharge correction factor (ratio of actual to ideal mass flow rate)
zeta_v = 0.928; % Velocity correction factor (sqrt of energy conversion efficiency)
eff_c_star = 0.85; % c star efficiency

% other data
N = 1; % Number of ports
rho_f = 1180; % Fuel density [kg/m^3]
t_b = 10; % Burn time [s]
F = 1250; % Average thrust [N] taken as input
R_f = 3/M_TO_IN / 2; % Grain outer radius [m]
L_grain = 20 / M_TO_IN; % grain length

%% CALCULATIONS

% % calculates nozzle exit pressure
% syms ep
% temp1 = ((gamma + 1)/2)^(1/(gamma -1));
% temp2 = (p_o/p_c)^(1/gamma);
% temp3 = (gamma + 1)/(gamma - 1);
% temp4 = 1 - (p_o/p_c)^((gamma -1)/gamma);
% eqn = 0 == temp1*temp2*sqrt(temp3*temp4) - 1/ep;
% [sol] = vpasolve(eqn,ep)

%% Calculations
% calculates nozzle exit pressure
syms p_e
temp1 = ((gamma + 1)/2)^(1/(gamma -1));
temp2 = (p_e/p_c)^(1/gamma);
temp3 = (gamma + 1)/(gamma - 1);
temp4 = 1 - (p_e/p_c)^((gamma -1)/gamma);
eqn = 0 == temp1*temp2*sqrt(temp3*temp4) - 1/ep;
[sol] = vpasolve(eqn,p_e);

for j = 1: size(sol,1)
    if isreal(sol(j)) && sol(j)>0
        p_e = double(vpa(sol(j))); % Nozzle exit pressure [Pa]
        break
    end
end

% check that the result is consistent with CEA data up to 2% error
assert(p_e/p_o-1.0 < 0.02,...
    'Difference between exit pressure and 1 atm exceeds 2%');

% calculates initial thrust coefficient 
temp1 = (2*gamma^2)/(gamma-1);
temp2 = (2/(gamma + 1))^((gamma+1)/(gamma-1));
temp3 = 1 - (p_e/p_c)^((gamma-1)/gamma);
temp4 = (p_e - p_o)*ep/p_c;
C_f = sqrt(temp1 * temp2 * temp3) - temp4; % Initial thrust coefficient

% efficiencies
eff_nozz = zeta_d * zeta_v; % Nozzle efficiency

% calculates initial throat area
A_t = F / (eff_nozz * C_f * p_c); % [m^2]
R_t = sqrt(A_t/pi);

% calculates initial propellant mass flow rate 
mdot = p_c*A_t / (eff_c_star * c_star); % [kg/s]
mdot_f = mdot / (of_ratio + 1); % [kg/s]
mdot_ox = of_ratio * mdot_f; % [kg/s]

% R_f = .171/2*39.37;
% a = .223;
% n = .5;
% mdot_ox = 1.511*2.204;

% calculates maximum grain inner radius (such that no fuel is left)
syms R_i
temp1 = (R_f) ^(2*n + 1) - (R_i) ^ (2*n + 1);
temp2  = a * (2*n + 1) * (mdot_ox/pi) ^ n;
eqn = t_b == temp1 / temp2;
sol = vpasolve(eqn, R_i);

for j = 1: size(sol,1)
    if isreal(sol(j)) && sol(j)>0
        R_i_max = double(vpa(sol(j))); % Grain inner radius [m]
        break
    end
end

% fix grain length and compute inner radius from it
% L = L_grain;
% R_i = ( mdot_f/(a*(mdot_ox/pi)^n) .* ...
%         1./(rho_f*2*pi*L)) .^ (1/(1-2*n));
% % check that we are below R_i_max
% assert(R_i<R_i_max,...
%     ['An inner radius greater than maximum allowable inner radius ',...
%     'fixed by burn time was given for desired thrust and grain length.']); 
% thickness = R_f - R_i;

% alternatively, fix R_i and compute L
R_i = 0.5 / M_TO_IN;
G_ox = mdot_ox/(pi*(R_i^2)); %[kg/(m^2*s)]
rdot = (a * (G_ox)^n); %[m/s]
L = mdot_f / (2 * pi * N * R_i * rho_f * rdot); %[m]
        
% calculates required fuel and oxidizer mass
m_f = rho_f * L*pi*(R_f^2 - R_i^2);
m_ox = mdot_ox * t_b; % [kg]

%% For plotting

% instantaneous O/F ratio
t = (0:.1:t_b);
inst_of_ratio = zeros(1,size(t,2));
for i = 1:size(t,2)
    temp1 = 1/(2*rho_f*L*a);
    temp2 = (mdot_ox/(pi*N))^(1-n);
    temp3 = (a * (2*n + 1) * (mdot_ox/(pi*N))^n * t(i) + R_i^(2*n + 1))^((2*n-1)/(2*n+1));
    inst_of_ratio(i) = temp1*temp2*temp3; % Instantaneous O/R ratio 
end

% instantaneous inner radius assuming constant mdot_ox 
t = (0:.1:t_b); R_i_inst = zeros(size(t));
for i = 1:size(t,2)
    R_i_inst(i) = (a*(2*n+1)*(mdot_ox/pi)^n*t(i) + R_i^(2*n+1))^(1/(2*n+1));
end

% sensitivity of R_i with L
L_range = (10:.5:20) / M_TO_IN;
R_i_range = ( mdot_f/(a*(mdot_ox/pi)^n) .* ...
            1./(rho_f*2*pi*L_range)) .^ (1/(1-2*n));

%% Plots
close all

% plot of R_i(t)
figure(1)
hold on
plot(t, R_i_inst * M_TO_IN)
plot(t, ones(length(t),1) * R_f * M_TO_IN, 'color', 'black')
xlabel('$t$ [s]'); ylabel('$R_i(t)$ [in]'); 

% plot of R_i(L)
figure(2)
hold on 
plot(L_range * M_TO_IN, R_i_range * M_TO_IN)
% annotations
X = [18, 18]; Y = [0, R_i*M_TO_IN];
plot(X, Y, 'color', 'k', 'linestyle', '--')
X = [0, 18]; Y = [R_i*M_TO_IN, R_i*M_TO_IN];
plot(X, Y, 'color', 'k', 'linestyle', '--')
xlabel('$L$ [in]'); ylabel('$R_i$ [in]'); title('Required $R_i$ given $T$ and $L$');
xlim([L_range(1)*M_TO_IN; L_range(end)*M_TO_IN])