Added test data from Winter 2022

Data Analysis/DataPlotsOlympus_03_12_22_Hotfire.m

@@ -0,0 +1,200 @@
+%% Data Plots - 03/12/2022 Hotfire
+clear
+close all
+
+
+% Unit Conversion
+g = 9.8066; %m/s2
+psi_to_Pa = 6894.75729; % 1 psi in Pa
+in_to_m = 0.0254; % 1 in in m
+mm_to_m = 1e-3; % 1 mm in m
+lbf_to_N = 4.44822162; % 1 lbf in N
+lbm_to_kg = 0.453592; % 1 lbm in kg
+atm_to_Pa = 101325; % 1 atm in Pa
+L_to_m3 = 1e-3; % 1 L in m^3
+
+%% Import Data
+load('../Test Data/03_12_22_Hotfire.mat')
+test_time = test_time - t_offset; % offset time
+dp_fuel = pft - pcc;
+dp_ox = pot - pcc;
+
+%% Inputs - Initial Values, start/stop times
+t_burn_start = 0.2;
+t_burn_stop = 11.89;
+t_test_start = -2;
+t_test_stop = 25;
+
+m_ox = 8.0; % kg, loaded mass of oxidizer (liq + gas)
+rho_fuel = 795; % kg/m3, density of fuel (kerosene)
+V_fueltank = 3.95E-3; % m3, volume of fuel tank
+V_fuel = 2.0e-3; % m3, volume of fuel loaded into tank
+m_fuel_remaining = 0.5616; %kg, mass of fuel remaining at end of burn (estimated from propsim)
+% ^^ this is only needed because ox ran out waaaayyy before fuel
+gamma = 1.4; % unitless, N2 adiabatic constant
+A_star = pi/4*(2.545e-2)^2; % m^2, area of nozzle throat
+
+% constants from PropSim - idealized performance from comparison run 
+c_star_theo = 1358.1; % m/s, C* value
+C_f_theo = 1.3531; % unitless, Cf value
+Isp_theo = c_star_theo*C_f_theo; % m/s^2, specific impulse (defined as C* times Cf)
+
+%% Calculated Constants - Values calced from inputs
+burn_ind = test_time>t_burn_start & test_time<t_burn_stop;
+t_burn_start_ind = find(burn_ind,1,'first');
+t_burn_stop_ind = find(burn_ind,1,'last');
+burn_time = test_time(t_burn_stop_ind) - test_time(t_burn_start_ind);
+test_ind = test_time>t_test_start & test_time<t_test_stop;
+
+V_fuelullage = V_fueltank-V_fuel; %m3, volume of fuel pressurant
+dm_fuel = rho_fuel*V_fuel - m_fuel_remaining; %kg,  fuel mass expended during burn
+dm_tot = dm_fuel+m_ox;
+pft_start = pft(t_burn_start_ind); %Pa, initial fuel tank pressure
+pft_stop = pft(t_burn_stop_ind); %Pa, final fuel tank pressure
+
+%% Corrections - load cell shift + tare
+ft_original = ft;
+ft_correction = zeros(size(ft));
+% Add zero shift from static friction (difference between mean + weight before and mean after)
+ft_initial_zero = dm_tot*g + mean(ft(t_burn_start_ind-40:t_burn_start_ind));
+ft_final_zero = mean(ft(t_burn_stop_ind+60:t_burn_stop_ind+60+40));
+ft_delta =   ft_initial_zero - ft_final_zero;
+ft_correction(t_burn_start_ind:end) = ft_correction(t_burn_start_ind:end) + ft_delta;
+% Remove weight, approximating roughly constant blowdown mdot
+ft_correction(1:t_burn_start_ind) = ft_correction(1:t_burn_start_ind) + dm_tot*g;
+ft_correction(burn_ind) = ft_correction(burn_ind) + linspace(dm_tot*g,0,t_burn_stop_ind-t_burn_start_ind+1)';
+% Tare
+ft_correction = ft_correction - ft_initial_zero;
+% Adjust ft
+ft_original = ft_original - ft_final_zero;
+ft = ft+ft_correction;
+
+%% Calculate fuel CdA from data
+m_fuel = rho_fuel * V_fuelullage*((pft_start./pft).^(1/gamma)-1);
+m_dot_fuel = [0; diff(m_fuel)./diff(test_time)];
+dm_fuel_test = rho_fuel * V_fuelullage*((pft_start/pft_stop)^(1/gamma)-1);
+fuel_scaling_factor = dm_fuel/dm_fuel_test;
+m_dot_fuel = fuel_scaling_factor*m_dot_fuel;
+fuel_cda = m_dot_fuel./(sqrt(2*rho_fuel)*sqrt(dp_fuel));
+fuel_cda_avg = dm_fuel./(sqrt(2*rho_fuel)*...
+    trapz(test_time(burn_ind),sqrt(dp_fuel(burn_ind))));
+
+%% Calculate performance characteristics (in time and avg)
+impulse = trapz(test_time(burn_ind),ft(burn_ind));
+pcc_int = trapz(test_time(burn_ind),pcc(burn_ind));
+pcc_diff = SimpleFilter(test_time,pcc,1.0,'high'); % high pass filter on pcc
+pcc_rms_int = sqrt(trapz(test_time(burn_ind),pcc_diff(burn_ind).^2));
+C_f = ft./(pcc*A_star);
+C_f_avg = impulse/(pcc_int*A_star);
+c_star_avg = (pcc_int*A_star)/(dm_tot);
+Isp_avg = impulse/(dm_tot);
+OF_avg = m_ox/dm_fuel;
+ft_avg = impulse/burn_time;
+pcc_avg = pcc_int/burn_time;
+
+%% Print out avg. performance 
+fprintf('Calculations:\n')
+fprintf('Impulse: %.1f kN*s (%.0f lbf*s)\n', impulse/1e3, impulse/lbf_to_N)
+fprintf('Burn Time: %.1f s\n', burn_time)
+fprintf('Avg. Thrust: %.2f kN (%.0f lbf)\n', ft_avg/1e3, ft_avg/lbf_to_N)
+fprintf('Avg. Pcc: %.2f MPa (%.0f psi)\n', pcc_avg/1e6, pcc_avg/psi_to_Pa)
+fprintf('Avg. Pcc RMS: %.2f %%\n', pcc_rms_int/pcc_int*1e2)
+fprintf('Isp: %.0f m/s (%.0f s)\n', Isp_avg, Isp_avg/9.81)
+fprintf('Isp efficiency: %.1f %%\n', Isp_avg/Isp_theo*1e2)
+fprintf('C_f: %.2f \n', C_f_avg)
+fprintf('C_f efficiency: %.1f %%\n', C_f_avg/C_f_theo*1e2)
+fprintf('C*: %.0f m/s (%.1f s)\n', c_star_avg, c_star_avg/9.81)
+fprintf('C* efficiency: %.1f %%\n', c_star_avg/c_star_theo*1e2)
+fprintf('OF (avg): %.1f\n', OF_avg)
+fprintf('Fuel Mass Spent: %.2f kg (%.2f lbm)\n', dm_fuel, dm_fuel/lbm_to_kg)
+fprintf('Ox. Mass Spent: %.2f kg (%.2f lbm)\n', m_ox, m_ox/lbm_to_kg)
+fprintf('Fuel CdA: %.2f mm^2\n', fuel_cda_avg*1e6)
+
+%% Pretty pretty plots of time data
+time = test_time;
+figure('Name', 'Pressures')
+plot((time(test_ind)), pot(test_ind)/psi_to_Pa,...
+    (time(test_ind)), pom(test_ind)/psi_to_Pa,...
+    (time(test_ind)), pft(test_ind)/psi_to_Pa,...
+    (time(test_ind)), pcc(test_ind)/psi_to_Pa)
+xlabel('Time (s)');
+ylabel('Pressures (psi)');
+legend({'Ox. Tank','Ox. Manifold','Fuel Tank','Combustion Chamber'})
+title('Engine Pressure Overview')
+
+
+figure('Name','CC Pressure')
+plot((time(test_ind)), pom(test_ind)/psi_to_Pa,...
+    (time(test_ind)), pft(test_ind)/psi_to_Pa,...
+    (time(test_ind)), pcc(test_ind)/psi_to_Pa)
+xlabel('Time (s)');
+ylabel('Pressures (psi)');
+legend({'Ox. Manifold','Fuel Tank','Combustion Chamber'})
+title('CC Pressure Overview')
+
+figure('Name','Pressure Drops')
+plot((time(burn_ind)),1e2*dp_fuel(burn_ind)./pft(burn_ind),...
+    (time(burn_ind)), 1e2*dp_ox(burn_ind)./pot(burn_ind))
+xlabel('Time (s)');
+ylabel('Pressure Drop (% of tank pressure)');
+legend({'Fuel','Oxidizer'})
+title('Pressure Drop over Injector')
+
+figure('Name','Thrust')
+plot(time(test_ind), ft_original(test_ind)/lbf_to_N, time(test_ind), ft(test_ind)/lbf_to_N)
+xlabel('Time (s)');
+ylabel('Thrust (lbf)');
+legend('Original','Corrected');
+title('Thrust')
+
+% Performance Parameters
+figure('Name','C_f')
+plot(time(burn_ind), C_f(burn_ind))
+xlabel('Time (s)');
+ylabel('Thrust Coefficient ()');
+title('C_f, Thrust Coefficient');
+
+% CdA Plots
+figure('Name','PFT')
+plot(time(burn_ind),pft(burn_ind)/psi_to_Pa)
+xlabel('Time (s)');
+ylabel('Pressure Fuel Tank  (psi)');
+title('Fuel Tank Pressure');
+
+figure('Name','Fuel Mass Consumed')
+plot(time(burn_ind),m_fuel(burn_ind))
+xlabel('Time (s)');
+ylabel('Fuel Mass Consumed (kg)');
+
+figure('Name', 'Fuel Mass Flow')
+plot(time(burn_ind),m_dot_fuel(burn_ind))
+xlabel('Time (s)');
+ylabel('Fuel Mass Flow Rate (kg/s)');
+title('Fuel Mass Flow Rate');
+
+figure('Name','Fuel CdA')
+plot(time(burn_ind),fuel_cda(burn_ind)*1e6)
+% ylim([0 2*median(fuel_cda(burn_ind))*1e6])
+xlabel('Time (s)');
+ylabel('Fuel CdA (mm^2)');
+title('Fuel CdA');
+
+FolderName = '../Outputs';   % Your destination folder
+FigList = findobj(allchild(0), 'flat', 'Type', 'figure');
+for iFig = 1:length(FigList)
+  FigHandle = FigList(iFig);
+  FigName   = get(FigHandle, 'Name');
+  saveas(FigHandle, fullfile(FolderName, strcat(FigName, '.png')));
+end
+
+
+%% Additional functions
+function smoothed = SimpleFilter(time,values,tau,type)
+    dt = median(diff(time));
+    Wn = dt/(tau); % cut-off frequency
+    values(isnan(values)) = 0;
+    [b, a] = butter(1,Wn,type); % butter-worth filter of order 1 of type specified (high/low or band)
+    smoothed = filter(b, a, values); % apply filter
+    shift_ind = round(1/(4*Wn));
+    smoothed(1:(end-shift_ind)) = smoothed((1+shift_ind):end);
+end