diff --git a/Data Analysis/DataPlotsOlympus_03_12_22_Hotfire.m b/Data Analysis/DataPlotsOlympus_03_12_22_Hotfire.m new file mode 100644 index 0000000..e8c61be --- /dev/null +++ b/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_timet_test_start & test_time= high_idx + savedata.(this_fn) = savedata.(this_fn)(low_idx:high_idx); + end + end + savedata.t_offset = max(0, savedata.t_offset - savedata.test_time(1)); + savedata.test_time = savedata.test_time - savedata.test_time(1); + savedata.tick = savedata.tick - savedata.tick(1); + hold off; + end + + if skip % check if user pressed skip + break + end + if reset_data % if user chose to reset, reset data to originally parsed data and re-start loop w/o saving + savedata = logdata; + continue + end + % Save struct to mat file + + if i>0 + name_mod = strcat('_',int2str(i)); + else + name_mod = ''; + end + i = i +1; + + % Add PropSim's preferred names to struct and filter + filter_arr = ones(window_size,1)/window_size; + savedata.pft = filter(filter_arr, 1,savedata.(pft_id)); + savedata.pot = filter(filter_arr, 1,savedata.(pot_id)); + savedata.pom = filter(filter_arr, 1,savedata.(pom_id)); + savedata.pcc = filter(filter_arr, 1,savedata.(pcc_id)); + savedata.ft = -filter(filter_arr, 1,savedata.(ft_id)); + savedata.ft_offset = mean(savedata.ft(end-10:end)); + save(fullfile(p,strcat(n,'_trimmed',name_mod,'.mat')), '-struct','savedata'); + + if new_selection % If user wants to do a new selection in addition to this one, repeat w/ original datas + savedata = logdata; + else + break; + end + end + close(h); +end + +function processLogFile(filename, required_sensors) + % filename: path/name of SD card txt file + % ft_id: fuel tank pressure id + % pcc_id: combustion chamber pressure id + % pom_id: ox manifold pressure id + % pot_id: ox tank ullage pressure id + num_solenoids = 8; + num_ematches = 8; + data = readlines(filename); + data_struct = struct(); + firsttime = true; + datalen = length(data); + delete_idx = []; + for i = 1:datalen + fline = data(i); + try + json_struct = jsondecode(fline); + fn = fieldnames(json_struct); + for k=1:numel(fn) + this_fn = fn{k}; + if isstruct(json_struct.(this_fn)) + subfn = fieldnames(json_struct.(this_fn)); + for j=1:numel(subfn) + this_subfn = subfn{j}; + if ~isfield(data_struct,this_subfn) + data_struct.(this_subfn) = zeros(datalen,1); + end + data_struct.(this_subfn)(i) = json_struct.(this_fn).(this_subfn); + end + elseif ischar(json_struct.(this_fn)) + if ~isfield(data_struct,this_fn) + data_struct.(this_fn) = cell(datalen,1); + end + data_struct.(this_fn){i} = {json_struct.(this_fn)}; + else + if ~isfield(data_struct,this_fn) + data_struct.(this_fn) = zeros(datalen,1); + end + data_struct.(this_fn)(i) = json_struct.(this_fn); + end + end + if firsttime + for l = 1:length(required_sensors) + if ~isfield(data_struct, required_sensors{l}) + return + end + end + end + catch + delete_idx = [delete_idx; i]; + end + end + fn = fieldnames(data_struct); + for k=1:numel(fn) + data_struct.(fn{k})(delete_idx) = []; + end + try + data_struct.test_time = (data_struct.tick)./1000; % convert to seconds from msec + svstates = int2bit((data_struct.SV).',num_solenoids, false).'; + emstates = int2bit((data_struct.EM).',num_ematches, false).'; + for i = 1:num_solenoids + data_struct.(strcat('SV',int2str(i))) = svstates(:,i); + end + for i = 1:num_ematches + data_struct.(strcat('EM',int2str(i))) = emstates(:,i); + end + + [p,n,~] = fileparts(filename); + save(fullfile(p,strcat(n,'.mat')), '-struct','data_struct'); + catch ME + error("Please download the communications toolbox!"); + end +end \ No newline at end of file diff --git a/Data Analysis/QuailDataPlotting.m b/Data Analysis/QuailDataPlotting.m new file mode 100644 index 0000000..b9e2472 --- /dev/null +++ b/Data Analysis/QuailDataPlotting.m @@ -0,0 +1,249 @@ +%% Quail Test Data Plotting +clearvars; +close all; +clc; + +[oxdatafile, oxpath] = uigetfile('*.mat','Select the ox data to be plotted.'); +oxdata = load(fullfile(oxpath,oxdatafile)); +[fueldatafile, fuelpath] = uigetfile('*.mat','Select the fuel data to be plotted.',oxpath); +fueldata = load(fullfile(fuelpath, fueldatafile)); +%% Constants and Channel Mapping +psiToPa = 6894.75729; % pressure conversion +KtoF = @(K) 1.8*(K-273) + 32; % temp conversion +Ntolbf = 0.224809; % force conversion +LTom3 = 1E-3; % volume conversion + +Pamb = 12.5*psiToPa; % Pa, ambient pressure +gamma_GN2 = 1.4; % adiabatic constant of nitrogen +rho_fuel = 1000; % kg/m3, density of fuel +V_fuel_loaded = 1.9*LTom3; % m3, volume of kerosene loaded +V_fuel_tank = 3.95*LTom3; % m3, volume of fuel tank +P_fuel_initial = 553*psiToPa; % Pa, set point of fuel (get from plots) +fuel_blowdown_time = 17.5; %sec, end of fuel blowdown + +fuel_tank_ch_id = 'PT2'; +ox_actuators = {{'Ox Fill', 'SV5'}; + {'Ox Bleed', 'SV6'}; + {'Ox Vent', 'SV8'}; + {'Ox Main Valve', 'SV2'}; + }; +ox_pressures = {{'Ox Pre-Fill', 'PT3'}; + {'Ox Post-Fill', 'PT4'}; + {'Ox Ullage', 'PT6'}; + {'Ox Manifold','PT5'}; + }; +fuel_actuators = {{'Fuel Fill', 'SV3'}; + {'Fuel Bleed', 'SV4'}; + {'Fuel Abort', 'SV7'}; + {'Fuel Main Valve', 'SV1'}; + }; +fuel_pressures = {{'Fuel Sump', 'PT2'}; + {'Fuel Source', 'PT1'}; + }; +cc_actuators = {{'Fuel Main Valve', 'SV1'}; + {'Ox Main Valve', 'SV2'}; + {'Igniter 1', 'EM1'}; + {'Igniter 2', 'EM2'}; + }; +cc_pressures = {{'CC', 'PT7'}}; +engine_pressures = {{'Fuel Sump', 'PT2'}; + {'Ox Ullage', 'PT6'}; + {'CC', 'PT7'}; + {'Ox Manifold','PT5'}; + }; +loads = {{'Right LC', 'LC1'}; + {'Left LC', 'LC2'}; + }; + +temps = {{'CC Wall', 'TC1'}; + {'Ambient', 'TC2'}; + }; + +%% Plotting + +% Fuel Stuff +fuelt = fueldata.test_time - fueldata.test_time(1); +if isfield(fueldata,'t_offset') + fuelt = fuelt - fueldata.t_offset; +end +figure('Name','Fuel Summary'); hold on; +title('Fuel Summary') +yyaxis left +xlabel('time, s'); +ylabel('pressure, psi'); +for i = 1:length(fuel_pressures) + this_map = fuel_pressures{i}; + plot(fuelt,fueldata.(this_map{2})/psiToPa,'DisplayName',strcat(this_map{2},' - ',this_map{1})); +end +yyaxis right +xlabel('time, s'); +ylabel('State'); +svs = zeros(length(fuel_actuators),1); +for i = 1:length(fuel_actuators) + this_map = fuel_actuators{i}; + svs(i) = plot(fuelt,fueldata.(this_map{2}),'DisplayName',strcat(this_map{2},' - ',this_map{1})); +end +legend; +set(gca,'ButtonDownFcn',@(src,~)togglevisibility(svs)) + +% Ox Stuff +oxt = oxdata.test_time - oxdata.test_time(1); +if isfield(oxdata,'t_offset') + oxt = oxt - oxdata.t_offset; +end + +figure('Name','Ox Summary'); hold on; +title('Ox Summary') +yyaxis left +xlabel('time, s'); +ylabel('pressure, psi'); +for i = 1:length(ox_pressures) + this_map = ox_pressures{i}; + plot(oxt,oxdata.(this_map{2})/psiToPa,'DisplayName',strcat(this_map{2},' - ',this_map{1})); +end +yyaxis right +xlabel('time, s'); +ylabel('State'); +svs = zeros(length(ox_actuators),1); +for i = 1:length(ox_actuators) + this_map = ox_actuators{i}; + svs(i) = plot(oxt,oxdata.(this_map{2}),'DisplayName',strcat(this_map{2},' - ',this_map{1})); +end +legend; +set(gca,'ButtonDownFcn',@(src,~)togglevisibility(svs)) + +% CC Stuff +oxt = oxdata.test_time - oxdata.test_time(1); +if isfield(oxdata,'t_offset') + oxt = oxt - oxdata.t_offset; +end + +figure('Name','CC Summary'); hold on; +title('CC Summary') +yyaxis left +xlabel('time, s'); +ylabel('pressure, psi'); +for i = 1:length(cc_pressures) + this_map = cc_pressures{i}; + plot(oxt,oxdata.(this_map{2})/psiToPa,'DisplayName',strcat(this_map{2},' - ',this_map{1})); +end +yyaxis right +xlabel('time, s'); +ylabel('State'); +svs = zeros(length(cc_actuators),1); +for i = 1:length(cc_actuators) + this_map = cc_actuators{i}; + svs(i) = plot(oxt,oxdata.(this_map{2}),'DisplayName',strcat(this_map{2},' - ',this_map{1})); +end +legend; +set(gca,'ButtonDownFcn',@(src,~)togglevisibility(svs)) + +% Engine Summry +oxt = oxdata.test_time - oxdata.test_time(1); +if isfield(oxdata,'t_offset') + oxt = oxt - oxdata.t_offset; +end + +figure('Name','Engine Summary'); hold on; +title('Engine Summary') +xlabel('time, s'); +ylabel('pressure, psi'); +for i = 1:length(engine_pressures) + this_map = engine_pressures{i}; + plot(oxt,oxdata.(this_map{2})/psiToPa,'DisplayName',strcat(this_map{2},' - ',this_map{1})); +end +legend; + +% Fuel CdA +time_idx = fuelt >= 0 & fuelt <= fuel_blowdown_time; +v = sqrt(2*(fueldata.(fuel_tank_ch_id)(time_idx)-Pamb)/rho_fuel); +V_ullage_initial = V_fuel_tank - V_fuel_loaded; +q = diff(V_ullage_initial*((P_fuel_initial./fueldata.(fuel_tank_ch_id)(time_idx)).^(1/gamma_GN2)-1))./diff(fuelt(time_idx)); +q = [q; q(end)]; +calced_fuelvol = trapz(fuelt(time_idx),q); +scaling_factor = V_fuel_loaded/calced_fuelvol; +window_size = 5; +fil = ones(window_size,1)./window_size; +v = filter(fil,1,v); +q = filter(fil,1,q); +cda = scaling_factor*q./v; +mean_cda = mean(cda); + +fprintf("Expected fuel volume : %0.2f L, Calculated fuel volume : %0.2f L \n", V_fuel_loaded/LTom3, calced_fuelvol/LTom3) +fprintf("Avg. Fuel CdA: %0.2f mm^2 \n", 1E6*mean_cda); + +figure('Name','Fuel CdA'); hold on; +title('Fuel CdA'); +xlabel('time, s'); +ylabel('CdA, mm^2'); +plot([0 fuel_blowdown_time],mean_cda*1E6*[1 1],'--') +plot(fuelt(time_idx),cda*1E6); + +% Loads +figure('Name','Loads'); hold on; +title('Loads'); +xlabel('time, s'); +ylabel('load, lbf'); +fil = ones(10,1)/10; +oxhold = 0; +fuelhold = 0; +for i = 1:length(loads) + this_map = loads{i}; + ox_filtered = filter(fil,1,oxdata.(this_map{2}))*Ntolbf; + oxhold = oxhold + ox_filtered; + fuel_filtered = filter(fil,1,fueldata.(this_map{2}))*Ntolbf; + fuelhold = fuelhold + fuel_filtered; + plot(oxt,ox_filtered,'DisplayName',strcat('Ox ', this_map{2},' - ',this_map{1})); + plot(fuelt,fuel_filtered,'DisplayName',strcat('Fuel ', this_map{2},' - ',this_map{1})); +end +plot(oxt,-(oxhold-mean(oxhold(end-5:end))),'DisplayName','Ox Total'); +plot(fuelt,fuelhold,'DisplayName','Fuel Total'); +legend; + +% Temps +figure('Name','Temps'); hold on; +title('Temps'); +xlabel('time, s'); +ylabel('temp, F'); +fil = ones(10,1)/10; +oxhold = 0; +fuelhold = 0; +for i = 1:length(temps) + this_map = temps{i}; + ox_filtered = filter(fil,1,KtoF(oxdata.(this_map{2}))); + oxhold = oxhold + ox_filtered; + fuel_filtered = filter(fil,1,KtoF(fueldata.(this_map{2}))); + plot(oxt,ox_filtered,'DisplayName',strcat('Ox ', this_map{2},' - ',this_map{1})); + plot(fuelt,fuel_filtered,'DisplayName',strcat('Fuel ', this_map{2},' - ',this_map{1})); +end +legend; + +% Battery +figure('Name','Battery'); hold on; +title('Battery Data'); +yyaxis left; +xlabel('time, s'); +ylabel('voltage, V'); +plot(oxt,oxdata.V_batt,'DisplayName','Ox Data Voltage'); +plot(fuelt, fueldata.V_batt,'DisplayName', 'Fuel Data Voltage'); + +yyaxis right; +ylabel('current, A'); +plot(oxt,oxdata.I_batt,'DisplayName','Ox Data Current'); +plot(fuelt, fueldata.I_batt,'DisplayName', 'Fuel Data Current'); +legend; + + +%% Prompt for merged save file. + + +function togglevisibility(vislist,~,~) + for i = 1:length(vislist) + h = vislist(i); + if strcmp(get(h,'visible'),'on') + set(h,'visible','off') + else + set(h,'visible','on') + end + end +end \ No newline at end of file diff --git a/Data Analysis/SimulateLiquid_ColdFlow_3_5_22.m b/Data Analysis/SimulateLiquid_ColdFlow_3_5_22.m new file mode 100644 index 0000000..df655a8 --- /dev/null +++ b/Data Analysis/SimulateLiquid_ColdFlow_3_5_22.m @@ -0,0 +1,156 @@ +%% Run Performance Code +cd('..') +addpath(fullfile(pwd, 'Supporting Functions')) + +clear +close all + +% Unit Conversion +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 + + +%% Options +test_data.test_plots_on = 1; % Import tests data and plot against simulation data +test_data.test_data_file = "03_04_22_ColdflowFuel"; % File from which to import test data +test_data.t_offset = -6.85; % Time offset of test data wrt simulation data [s] + +% 1: simulate combustion (hot fire), 0: no combustion (cold flow) +mode.combustion_on = 0; +% 1: simulate flight conditions (i.e. acceleration head), 0: ground test +% conditions +mode.flight_on = 0; +% 'hybrid' for solid fuel, liquid oxidizer, 'liquid' for liquid fuel and +% oxidizer +mode.type = 'liquid'; + +%% Input Parameters + +inputs.CombustionData = fullfile('.\Combustion Data', 'CombustionData_T1_N2O.mat'); + +%-------Gases----------------------- + +helium = Gas(); +helium.c_v = 3.12; % J/kg*K +helium.molecular_mass = 4.0026e-3; % kg/mol + +nitrogen = Gas(); +nitrogen.c_v = 0.743e3; % J/kg*K +nitrogen.molecular_mass = 2*14.0067e-3; % kg/mol + +%-------Injector Properties---------- + +%Injector Exit Area +inputs.ox.injector_area = 28.5*0.95*mm_to_m^2; +inputs.fuel.injector_area = 2.63*mm_to_m^2; + +% Ball Valve Time to Injector Area (s) +inputs.dt_valve_open = 0.01; + +%Discharge Coefficient +inputs.ox.Cd_injector = 1; +inputs.fuel.Cd_injector = 1; + +%-------Rocket Properties-------- +%Rocket Dry Mass +inputs.mass_dry_rocket = 50*lbm_to_kg; + +%-------Oxidizer Properties-------- +%Tank Volume +inputs.ox.V_tank = 13.4*L_to_m3; + +%Nitrous Volume +inputs.ox.V_l = 12.5*L_to_m3; + +%Tank Inner Diameter +inputs.ox.tank_id = 4.75*in_to_m; + +%Distance from Bottom of Tank to Injector +inputs.ox.h_offset_tank = 13*in_to_m; + +%Main Flow Line Diameter +inputs.ox.d_flowline = .4*in_to_m; + +%Tank Temperature (K) +inputs.ox.T_tank = N2O_FindT(538*psi_to_Pa); + +%-------Oxidizer Pressurant Properties-------- + +inputs.ox_pressurant = Pressurant('oxidizer'); +inputs.ox_pressurant.gas_properties = helium; +inputs.ox_pressurant.set_pressure = 800*psi_to_Pa; +inputs.ox_pressurant.storage_initial_pressure = 4500*psi_to_Pa; +inputs.ox_pressurant.tank_volume = 3.5*L_to_m3; +inputs.ox_pressurant.flow_CdA = 8*mm_to_m^2; + +%Are You Supercharging? (0 for 'NO' 1 for 'YES') +inputs.ox_pressurant.active = 0; + +%-------Fuel Properties-------- + +%Tank Volume +inputs.fuel.V_tank = 3.95*L_to_m3; + +%Nitrous Volume +inputs.fuel.V_l = 1.75*L_to_m3; + +%Tank Inner Diameter +inputs.fuel.tank_id = 4.75*in_to_m; + +%Distance from Bottom of Tank to Injector +inputs.fuel.h_offset_tank = 2*in_to_m; + +%Main Flow Line Diameter(in) +inputs.fuel.d_flowline = .25*in_to_m; + +inputs.fuel.rho = 1000; %Kg/m^3 + +%-------Fuel Pressurant Properties-------- + +inputs.fuel_pressurant = Pressurant('fuel'); +inputs.fuel_pressurant.gas_properties = nitrogen; +inputs.fuel_pressurant.set_pressure = 500*psi_to_Pa; +inputs.fuel_pressurant.storage_initial_pressure = 4500*psi_to_Pa; +inputs.fuel_pressurant.tank_volume = 0.0*L_to_m3; +inputs.fuel_pressurant.flow_CdA = 8*mm_to_m^2; + +%Are You Supercharging? (0 for 'NO' 1 for 'YES') +inputs.fuel_pressurant.active = 1; + +%-------Other Properties-------- + +%Combustion chamber dimensions +inputs.length_cc = 10*in_to_m; +inputs.d_cc = 3.75*in_to_m; + +%Estimated nozzle efficiency +inputs.nozzle_efficiency = 0.95; +inputs.nozzle_correction_factor = 0.9830; + +% Estimated combustion efficiency +inputs.c_star_efficiency = 0.87; + +%Nozzle Throat diameter +inputs.d_throat = 2.56e-2; + +%Expansion Ratio +inputs.exp_ratio = 3.5; + +%Ambient Temperature +inputs.T_amb = convtemp(60,'F','K'); + +%Ambient Pressure +inputs.p_amb = 12.5*psi_to_Pa; + +% Load Combustion Data +inputs.comb_data = load(inputs.CombustionData); +inputs.comb_data = inputs.comb_data.CombData; + +%% Run Performance Code +PerformanceCode(inputs, mode, test_data); +cd('Data Analysis') diff --git a/Data Analysis/SimulateLiquid_HotFire_6_8_19.m b/Data Analysis/SimulateLiquid_HotFire_6_8_19.m index 41c2362..c8874b1 100644 --- a/Data Analysis/SimulateLiquid_HotFire_6_8_19.m +++ b/Data Analysis/SimulateLiquid_HotFire_6_8_19.m @@ -98,7 +98,7 @@ inputs.fuel.V_tank = 1.88*L_to_m3; %Nitrous Volume -inputs.fuel.V_l = 0.415*L_to_m3; +inputs.fuel.V_l = 0.815*L_to_m3; %Tank Inner Diameter inputs.fuel.tank_id = 3.75*in_to_m; @@ -115,7 +115,7 @@ inputs.fuel_pressurant = Pressurant('fuel'); inputs.fuel_pressurant.gas_properties = nitrogen; -inputs.fuel_pressurant.set_pressure = 715*psi_to_Pa; +inputs.fuel_pressurant.set_pressure = 750*psi_to_Pa; inputs.fuel_pressurant.storage_initial_pressure = 4500*psi_to_Pa; inputs.fuel_pressurant.tank_volume = 0.0*L_to_m3; inputs.fuel_pressurant.flow_CdA = 8*mm_to_m^2; diff --git a/Data Analysis/SimulateLiquid_Hotfire_3_12_22.m b/Data Analysis/SimulateLiquid_Hotfire_3_12_22.m new file mode 100644 index 0000000..914a32d --- /dev/null +++ b/Data Analysis/SimulateLiquid_Hotfire_3_12_22.m @@ -0,0 +1,157 @@ +%% Run Performance Code +cd('..') +addpath(genpath(fullfile(pwd, 'Supporting Functions'))) + +clear +close all + +% Unit Conversion +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 + + +%% Options +test_data.test_plots_on = 1; % Import tests data and plot against simulation data +test_data.test_data_file = "03_12_22_Hotfire"; % File from which to import test data +test_data.t_offset = -0.2; % Time offset of test data wrt simulation data [s] +test_data.ft_offset = 26*lbf_to_N; + +% 1: simulate combustion (hot fire), 0: no combustion (cold flow) +mode.combustion_on = 1; +% 1: simulate flight conditions (i.e. acceleration head), 0: ground test +% conditions +mode.flight_on = 0; +% 'hybrid' for solid fuel, liquid oxidizer, 'liquid' for liquid fuel and +% oxidizer +mode.type = 'liquid'; + +%% Input Parameters + +inputs.CombustionData = fullfile('.\Combustion Data', 'CombustionData_T1_N2O.mat'); + +%-------Gases----------------------- + +helium = Gas(); +helium.c_v = 3.12; % J/kg*K +helium.molecular_mass = 4.0026e-3; % kg/mol + +nitrogen = Gas(); +nitrogen.c_v = 0.743e3; % J/kg*K +nitrogen.molecular_mass = 2*14.0067e-3; % kg/mol + +%-------Injector Properties---------- + +%Injector Exit Area +inputs.ox.injector_area = 28.5*0.95*mm_to_m^2; +inputs.fuel.injector_area = 1.85*mm_to_m^2; + +% Ball Valve Time to Injector Area (s) +inputs.dt_valve_open = 0.01; + +%Discharge Coefficient +inputs.ox.Cd_injector = 1; +inputs.fuel.Cd_injector = 1; + +%-------Rocket Properties-------- +%Rocket Dry Mass +inputs.mass_dry_rocket = 50*lbm_to_kg; + +%-------Oxidizer Properties-------- +%Tank Volume +inputs.ox.V_tank = 13.4*L_to_m3; + +%Nitrous Volume +inputs.ox.V_l = 8.75*L_to_m3; + +%Tank Inner Diameter +inputs.ox.tank_id = 4.75*in_to_m; + +%Distance from Bottom of Tank to Injector +inputs.ox.h_offset_tank = 13*in_to_m; + +%Main Flow Line Diameter +inputs.ox.d_flowline = .4*in_to_m; + +%Tank Temperature (K) +inputs.ox.T_tank = N2O_FindT(585*psi_to_Pa); + +%-------Oxidizer Pressurant Properties-------- + +inputs.ox_pressurant = Pressurant('oxidizer'); +inputs.ox_pressurant.gas_properties = helium; +inputs.ox_pressurant.set_pressure = 800*psi_to_Pa; +inputs.ox_pressurant.storage_initial_pressure = 4500*psi_to_Pa; +inputs.ox_pressurant.tank_volume = 3.5*L_to_m3; +inputs.ox_pressurant.flow_CdA = 8*mm_to_m^2; + +%Are You Supercharging? (0 for 'NO' 1 for 'YES') +inputs.ox_pressurant.active = 0; + +%-------Fuel Properties-------- + +%Tank Volume +inputs.fuel.V_tank = 3.95*L_to_m3; + +%Nitrous Volume +inputs.fuel.V_l = 1.8*L_to_m3; + +%Tank Inner Diameter +inputs.fuel.tank_id = 4.75*in_to_m; + +%Distance from Bottom of Tank to Injector +inputs.fuel.h_offset_tank = 2*in_to_m; + +%Main Flow Line Diameter(in) +inputs.fuel.d_flowline = .25*in_to_m; + +inputs.fuel.rho = 795; %Kg/m^3 + +%-------Fuel Pressurant Properties-------- + +inputs.fuel_pressurant = Pressurant('fuel'); +inputs.fuel_pressurant.gas_properties = nitrogen; +inputs.fuel_pressurant.set_pressure = 567*psi_to_Pa; +inputs.fuel_pressurant.storage_initial_pressure = 4500*psi_to_Pa; +inputs.fuel_pressurant.tank_volume = 0.0*L_to_m3; +inputs.fuel_pressurant.flow_CdA = 8*mm_to_m^2; + +%Are You Supercharging? (0 for 'NO' 1 for 'YES') +inputs.fuel_pressurant.active = 1; + +%-------Other Properties-------- + +%Combustion chamber dimensions +inputs.length_cc = 10*in_to_m; +inputs.d_cc = 3.75*in_to_m; + +%Estimated nozzle efficiency +inputs.nozzle_efficiency = 0.95; +inputs.nozzle_correction_factor = 0.9830; + +% Estimated combustion efficiency +inputs.c_star_efficiency = 0.87; + +%Nozzle Throat diameter +inputs.d_throat = 2.56e-2; + +%Expansion Ratio +inputs.exp_ratio = 3.5; + +%Ambient Temperature +inputs.T_amb = convtemp(60,'F','K'); + +%Ambient Pressure +inputs.p_amb = 12.5*psi_to_Pa; + +% Load Combustion Data +inputs.comb_data = load(inputs.CombustionData); +inputs.comb_data = inputs.comb_data.CombData; + +%% Run Performance Code +PerformanceCode(inputs, mode, test_data); +cd('Data Analysis') diff --git a/Supporting Functions/Combustion Chamber/CombustionCalc.m b/Supporting Functions/Combustion Chamber/CombustionCalc.m new file mode 100644 index 0000000..32a623c --- /dev/null +++ b/Supporting Functions/Combustion Chamber/CombustionCalc.m @@ -0,0 +1,55 @@ +function [gamma_ex, iMW, R_ex, iTc, c_star] = ... + CombustionCalc(OF,Pc,inputs,CombData) +%CombustionCalc Calculate combustion properties. +% Inputs: +% OF: oxidizer to fuel mass flow ratio +% Pc: combustion chamber static pressure +% inputs: motor input structure +% CombDataFit: fit object with combustion properties over range of OF +% and chamber pressure +% Outputs: +% gamma_ex: ratio of specific heats of combustion gases +% iMW: molecular mass of combustion gases +% R_ex: exhaust specific gas constant +% iTc: stagnation temperature of combustion gases +% c_star: characteristic velocity of combustion + +%% Constants +R_u = 8.3144621; % Universal gas constant [J/mol*K] + +%% Calculate Combustion Properties +% Correct OF ratio or chamber pressure if over boundaries +if OF < CombData.OF_range(1) + OF = CombData.OF_range(1); +elseif OF > CombData.OF_range(2) + OF = CombData.OF_range(2); +end +if Pc < CombData.Pc_range(1) + Pc = CombData.Pc_range(1); +elseif Pc > CombData.Pc_range(2) + Pc = CombData.Pc_range(2); +end + +% Use empirical formula to calculate chamber temperature, exhaust +% properties +iTc = FastInterp2(CombData.OF, CombData.Pc, CombData.Tc, OF, Pc); +iMW = FastInterp2(CombData.OF, CombData.Pc, CombData.M, OF, Pc); +R_ex = R_u/iMW; +gamma_ex = FastInterp2(CombData.OF, CombData.Pc, CombData.gamma, OF, Pc); +c_star = FastInterp2(CombData.OF, CombData.Pc, CombData.c_star, OF, Pc); + +% Apply c-star efficiency +c_star = c_star*inputs.c_star_efficiency; +iTc = iTc*inputs.c_star_efficiency^2; + +% Check Tc, MW, gamma_ex for invalid values +if (iTc < 0) || ~isreal(iTc) || isnan(iTc) + error('Invalid Tc calculated: Tc = %.3f', iTc); +end +if (iMW < 0) || ~isreal(iMW) || isnan(iMW) + error('Invalid MW calculated: MW = %.3f', iMW); +end +if (gamma_ex < 0) || ~isreal(gamma_ex) || isnan(gamma_ex) + error('Invalid gamma_ex calculated: gamma_ex = %.3f', gamma_ex); +end +end \ No newline at end of file diff --git a/Supporting Functions/Integration.m b/Supporting Functions/Integration.m index 0a4e8b4..4148793 100644 --- a/Supporting Functions/Integration.m +++ b/Supporting Functions/Integration.m @@ -691,62 +691,6 @@ function LiquidRecord(time, state, state_dot, F_thrust, p_crit, ... T_dot_drain = dW/(state.m_fueltank_press*inputs.fuel_pressurant.gas_properties.c_v); end -function [gamma_ex, iMW, R_ex, iTc, c_star] = ... - CombustionCalc(OF,Pc,inputs,CombData) -%CombustionCalc Calculate combustion properties. -% Inputs: -% OF: oxidizer to fuel mass flow ratio -% Pc: combustion chamber static pressure -% inputs: motor input structure -% CombDataFit: fit object with combustion properties over range of OF -% and chamber pressure -% Outputs: -% gamma_ex: ratio of specific heats of combustion gases -% iMW: molecular mass of combustion gases -% R_ex: exhaust specific gas constant -% iTc: stagnation temperature of combustion gases -% c_star: characteristic velocity of combustion - -%% Constants -R_u = 8.3144621; % Universal gas constant [J/mol*K] - -%% Calculate Combustion Properties -% Correct OF ratio or chamber pressure if over boundaries -if OF < CombData.OF_range(1) - OF = CombData.OF_range(1); -elseif OF > CombData.OF_range(2) - OF = CombData.OF_range(2); -end -if Pc < CombData.Pc_range(1) - Pc = CombData.Pc_range(1); -elseif Pc > CombData.Pc_range(2) - Pc = CombData.Pc_range(2); -end - -% Use empirical formula to calculate chamber temperature, exhaust -% properties -iTc = FastInterp2(CombData.OF, CombData.Pc, CombData.Tc, OF, Pc); -iMW = FastInterp2(CombData.OF, CombData.Pc, CombData.M, OF, Pc); -R_ex = R_u/iMW; -gamma_ex = FastInterp2(CombData.OF, CombData.Pc, CombData.gamma, OF, Pc); -c_star = FastInterp2(CombData.OF, CombData.Pc, CombData.c_star, OF, Pc); - -% Apply c-star efficiency -c_star = c_star*inputs.c_star_efficiency; -iTc = iTc*inputs.c_star_efficiency^2; - -% Check Tc, MW, gamma_ex for invalid values -if (iTc < 0) || ~isreal(iTc) || isnan(iTc) - error('Invalid Tc calculated: Tc = %.3f', iTc); -end -if (iMW < 0) || ~isreal(iMW) || isnan(iMW) - error('Invalid MW calculated: MW = %.3f', iMW); -end -if (gamma_ex < 0) || ~isreal(gamma_ex) || isnan(gamma_ex) - error('Invalid gamma_ex calculated: gamma_ex = %.3f', gamma_ex); -end -end - function [T_dot_press, m_dot_press] = PressurantFlow(p_presstank, ... T_presstank, p_downtank, T_downtank, m_cv_downtank, pressurant) %PressurantFlow Models flow of pressurant gas from source tank to diff --git a/Supporting Functions/PerformanceCode.m b/Supporting Functions/PerformanceCode.m index 68fc3a0..bbb9033 100644 --- a/Supporting Functions/PerformanceCode.m +++ b/Supporting Functions/PerformanceCode.m @@ -141,6 +141,29 @@ if test_data.test_plots_on % Load Imported Data for comparison [test_time, pft, pom, pot, we, ft, pcc] = LoadDataVars(test_data.test_data_file, test_data.t_offset); + if isfield(test_data,'ft_offset') + ft = ft + test_data.ft_offset; + end + if ~isempty(pcc) % calc thrust! + ft_calc = zeros(size(test_time)); + for j = 1:length(test_time) + [~, t_ind] = min(abs(test_time(j)-time)); % find index of sim time to use sim OF val + area_throat = pi*(inputs.d_throat)^2/4; + d_exit = inputs.d_throat*sqrt(inputs.exp_ratio); + exit_area = area_throat*inputs.exp_ratio; + thisOF = OF(t_ind); + if ~isnan(thisOF) && pcc(j) > inputs.p_amb + [gamma_ex, iMW, ~, iTc, c_star] = ... + CombustionCalc(thisOF,pcc(j),inputs,inputs.comb_data); + [ ~, p_exit, p_shock, u_exit, ~, M_e ] = NozzleCalc( inputs.d_throat, ... + d_exit, iTc, pcc(j), gamma_ex, iMW, inputs.p_amb); + m_dot_ex = pcc(j)*area_throat/c_star; + C_T = inputs.nozzle_correction_factor*(m_dot_ex*u_exit + ... + (p_exit-inputs.p_amb)*exit_area)/(area_throat*pcc(j)); + ft_calc(j) = c_star*C_T*m_dot_ex; + end + end + end else test_time = []; pft = []; @@ -148,6 +171,7 @@ pot = []; we = []; ft = []; + ft_calc = []; pcc = []; end @@ -165,9 +189,14 @@ if mode.combustion_on axes('parent',tab1) plot(time, F_thrust./lbf_to_N,... - test_time,ft./lbf_to_N) + test_time,ft_calc./lbf_to_N,test_time,ft./lbf_to_N) xlabel('Time (s)') ylabel('Thrust Force (lbf)') + if isempty(test_time) + legend('Simulation') + else + legend('Simulation', 'Calculated','Measured') + end end axes('parent',tab2) @@ -377,22 +406,47 @@ end function [test_time, pft, pom, pot, we, ft, pcc] = LoadDataVars(filename, t_offset) - test_time = 0; - pft = 0; - pom = 0; - pot = 0; - we = 0; - ft = 0; - pcc = 0; + pft = []; + pom = []; + pot = []; + we = []; + ft = []; + pcc = []; if contains(filename, '/') || contains(filename, '\') - load(filename) + s = load(filename); else - load(['.\Test Data\' filename]) + s = load(strcat('.\Test Data\',filename)); end - test_time = test_time + t_offset; - if length(test_time) == 1 + try + if isfield(s,'t_offset') + test_time = s.test_time - s.test_time(1) - s.t_offset + t_offset; + else + test_time = s.test_time- s.test_time(1) + t_offset; + end + catch error('Variable ''test_time'' not present in imported data.\n'); end + if isfield(s,'pft') + pft = s.pft; + end + if isfield(s,'pom') + pom = s.pom; + end + if isfield(s,'pot') + pot = s.pot; + end + if isfield(s,'we') + we = s.we; + end + if isfield(s,'ft') + ft = s.ft; + end + if isfield(s,'ft_offset') + ft = ft - s.ft_offset; + end + if isfield(s,'pcc') + pcc = s.pcc; + end % Create emp_oxtanky vectors for missing data % n_indices = length(test_time); diff --git a/Test Data/03_04_22_ColdflowFuel.mat b/Test Data/03_04_22_ColdflowFuel.mat new file mode 100644 index 0000000..f5e701e Binary files /dev/null and b/Test Data/03_04_22_ColdflowFuel.mat differ diff --git a/Test Data/03_06_22_ColdflowOx.mat b/Test Data/03_06_22_ColdflowOx.mat new file mode 100644 index 0000000..2143967 Binary files /dev/null and b/Test Data/03_06_22_ColdflowOx.mat differ diff --git a/Test Data/03_12_22_Hotfire.mat b/Test Data/03_12_22_Hotfire.mat new file mode 100644 index 0000000..76d1943 Binary files /dev/null and b/Test Data/03_12_22_Hotfire.mat differ