diff --git a/RoccarasoFlight/postprocessing/RPS/cleanData/ROC2/engineFlightData.mat b/RoccarasoFlight/postprocessing/RPS/cleanData/ROC2/engineFlightData.mat index cfd6deb231917c80950f74c2dbfb8ad51dc084f0..aef81bab5ef6ee0f4df703b96c10110556d95a46 100644 Binary files a/RoccarasoFlight/postprocessing/RPS/cleanData/ROC2/engineFlightData.mat and b/RoccarasoFlight/postprocessing/RPS/cleanData/ROC2/engineFlightData.mat differ diff --git a/RoccarasoFlight/postprocessing/RPS/configDataAnalyzer.m b/RoccarasoFlight/postprocessing/RPS/configDataAnalyzer.m index 5e45f0fb52fdfdf25262ec17e7b29681e7982524..441dcc97598c2c3f464aa411ec2b11bf716684e7 100644 --- a/RoccarasoFlight/postprocessing/RPS/configDataAnalyzer.m +++ b/RoccarasoFlight/postprocessing/RPS/configDataAnalyzer.m @@ -1,10 +1,9 @@ %% FLIGHT and ACQUISITION SETTINGS -settings.launchSite = 'Roccaraso'; % [-] Launch site name settings.testName = 'ROC2'; % [-] Test name settings.frPT = 1000; % [Hz] PTs sampling rate settings.frSrv = 10; % [Hz] SRVs sampling rate -settings.fileCSVpath = strcat('../../../', settings.launchSite, '/data/MOTORE/'); +settings.fileCSVpath = '../../data/MOTORE/'; % Load mission parameters run(strcat('engineData/missions/', settings.testName)); @@ -16,6 +15,20 @@ settings.engineName = 'Furia'; % Retrieve engine data run(strcat('engineData/engines/', settings.engineName, '/', settings.engineName)); +% Load CEA database +load('engineData/engines/CEA_db.mat'); +engine.CEA_db = CEA_db; +clear CEA_db + +%% SIMULATION SETTINGS +settings.interpMthd = 1; + +%% CONSTANTS +settings.g0 = 9.806; % [m/s2] Acceleration of Gravity, constant +settings.R = 8314.51; % [kJ/mol-K] +settings.z = 1410; % [m] Altitude above-sea-level +[~, ~, settings.Pamb, ~] = atmosisa(settings.z); % [Pa] Ambient initial pressure + %% ODE SETTINGS % !! Do not change unless you know what you are doing !! settings.odeEngineSettings = odeset('RelTol', 1e-5, 'AbsTol', 1e-5, 'Events', @eventBurnOut); @@ -28,5 +41,4 @@ settings.bisection.xTol = 1e-3; settings.bisection.absTol = 1e-3; %% PLOT SETTINGS -settings.plots = true; % [-] True if you want to see the plots - +settings.plots = true; % [-] True if you want to see the plots \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m b/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m index b3f9c00e32f4fa8641288d81ef355cb38d87864e..ca73d482882adb77bf26974dfaa359276ef2373b 100644 --- a/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m +++ b/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m @@ -30,9 +30,20 @@ end %% DATA PROCESSING AND ANALYSIS engine.acquisition = flightData; -engine = tankEvacuation(engine, settings); +% Computing the mass flow rate of oxidizer +engine = tankEvacuation(engine); + +% Simulating the engine +engineSim = combustionSimulation(engine, settings); + +% Computing actual performance +engineFly = performanceReconstruction(engine, engineSim); %% PLOTTING +fprintf('DATA ANALYSIS CONCLUDED: \n\n') +fprintf('\t - Total Impulse simulated: %.2f\n', engineSim.performances.Itot) +fprintf('\t - Total Impulse computed: %.2f\n', engineFly.Itot) + if settings.plots f1 = figure; yyaxis left @@ -45,6 +56,21 @@ if settings.plots stairs(flightData.timeSrvMain, flightData.srvMain) xlabel('Time [s]'); ylabel('SRV Opening [-]'); title('Tank pressure values') legend('Top tank', 'Bottom tank', 'Combustion chamber') + axis padded + + f2 = figure; + plot(engineSim.time, engineSim.mfr.mp) + grid minor + ylabel('Mass Flow [kg/s]'); xlabel('Time [s]'); title('Propellant MFR') + axis padded + + f3 = figure; + plot(engineSim.time, engineSim.performances.T) + hold on; grid minor + plot(engineFly.time, engineFly.T) + ylabel('Thrust [N]'); xlabel('Time [s]'); title('Thrust compare') + legend('Simulated', 'Computed') + axis padded end diff --git a/RoccarasoFlight/postprocessing/RPS/engineData/engines/CEA_db.mat b/RoccarasoFlight/postprocessing/RPS/engineData/engines/CEA_db.mat index 2712efdac18d9d0c55b4166d865a2d70d69286e8..92de427b9ecb77f7f42a48eb7eda47fdef275435 100644 Binary files a/RoccarasoFlight/postprocessing/RPS/engineData/engines/CEA_db.mat and b/RoccarasoFlight/postprocessing/RPS/engineData/engines/CEA_db.mat differ diff --git a/RoccarasoFlight/postprocessing/RPS/engineData/engines/Furia/Furia.m b/RoccarasoFlight/postprocessing/RPS/engineData/engines/Furia/Furia.m index 09f094e2b292763bb2615892ff7fae2645ea4120..7190b74c04cf1b46d10ba6f40da6562379955075 100644 --- a/RoccarasoFlight/postprocessing/RPS/engineData/engines/Furia/Furia.m +++ b/RoccarasoFlight/postprocessing/RPS/engineData/engines/Furia/Furia.m @@ -44,10 +44,12 @@ engine.data.d_pipe = 1.25e-2; % [m] Diameter of the feedline pipes engine.data.dInj = 1.5e-3; % [m] Diameter of injection holes engine.data.Aox = engine.data.Ninj* ... pi/4*engine.data.dInj.^2; % [m^2] Injection area - +engine.data.alphaFit = [-11.7668580891917, 21.5544619916252, -8.85332436329970]; % [-] Scaling factor between top tank pressure and injector pressure +engine.data.alphaDeg = length(engine.data.alphaFit) - 1; + %%% EFFICIENCIES -engine.data.etaCT = 0.96; -engine.data.etaCs = 0.85; +engine.data.etacT = 0.97; +engine.data.etaCs = 0.84; %%% COMBUSTION CHAMBER engine.data.rt = 0.0149; % [m] Throat radius @@ -56,6 +58,17 @@ engine.data.lambda = 0.985; % [-] Loss due to nozzle divergence engine.data.Rc_CM = 0.5; engine.data.R_pre = 1.15; % [-] Ratio between pre-chamber length and port radius engine.data.R_post = 2.35; % [-] ratio between post-chamber length and port radius +engine.data.Pc0 = 10; % [bar] Initial pressure value + +%%% PRE CC +engine.data.Dpre = 66e-3; +engine.data.Lpre = 29.3e-3; +engine.data.Vpre = pi/4*engine.data.Dpre^2*engine.data.Lpre; + +%%% POST CC +engine.data.Dpost = 62e-3; +engine.data.Lpost = 60e-3; +engine.data.Vpost = pi/4*engine.data.Dpost^2*engine.data.Lpost; %%% NOZZLE engine.data.zopt = 0; % [m] Altitude of optimal expansion diff --git a/RoccarasoFlight/postprocessing/RPS/src/Ap2RegressionCoor.m b/RoccarasoFlight/postprocessing/RPS/src/Ap2RegressionCoor.m new file mode 100644 index 0000000000000000000000000000000000000000..4324600f74fb39d8559fb2f0c19dff1ca67cd459 --- /dev/null +++ b/RoccarasoFlight/postprocessing/RPS/src/Ap2RegressionCoor.m @@ -0,0 +1,30 @@ +function [r, dr, Ab] = Ap2RegressionCoor(Ap, dAp, engine) +%{ +Ap2RegressionCoor - this function takes the port area of the grain and +gives back the regression coordinate and its derivative. This is done +accordin to the grain type specified in the structure engine. + +CALLED SCRIPTS: + +CALLED FUNCTIONS: + +CALLED DATA FILES: + +REVISIONS: +- 0, 12/09/2022, Release, Alberto Boffi, Paolo Gnata, PRP +Department + +Copyright © 2022, Skyward Experimental Rocketry, PRP department +All rights reserved + +SPDX-License-Identifier: GPL-3.0-or-later + +%} + +%% INPUT +L = engine.data.Lgrain; + +%% CALCULATIONS +r = sqrt(Ap./pi); +dr = dAp./(2*pi.*r); +Ab = L*2*pi.*r; \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/RPS/src/combustionSimulation.m b/RoccarasoFlight/postprocessing/RPS/src/combustionSimulation.m new file mode 100644 index 0000000000000000000000000000000000000000..19cad597d0dd54d74ed290e02fa2abc746045092 --- /dev/null +++ b/RoccarasoFlight/postprocessing/RPS/src/combustionSimulation.m @@ -0,0 +1,136 @@ +function engineSim = combustionSimulation(engine, settings) +%{ +combustionCompare - this function simulates combustion chamber behavior +using estimated/measured mox curve and compares it with those resulting +from combustion with most reliable a, n coefficients retrieved until now + +CALLED SCRIPTS: / + +CALLED FUNCTIONS: / + +CALLED DATA FILES: / + +REVISIONS: +-1, 20/04/2023, Release, Alberto Boffi, PRP Department + +Copyright © 2023, Skyward Experimental Rocketry, PRP Department +All rights reserved + +SPDX-License-Identifier: GPL-3.0-or-later + +%} + +%% INPUT + +%%% GRAIN +L = engine.data.Lgrain; + +%%% ENGINE +epsilon = engine.data.epsilon; +At = pi * engine.data.rt^2; +lambda = engine.data.lambda; +etacT = engine.data.etacT; + +%%% PROPELLANT +rhoFu = engine.data.fu.rhof; + +%%% ENGINE GUESSES +Ap0 = pi*engine.data.r0^2; +Pc0 = engine.data.Pc0*1e5; + +%%% SETTINGS +fsolveOptions = settings.fsolveSettings; +tSpan = engine.acquisition.cutData.timePTop - engine.acquisition.cutData.timePTop(1); +Pamb = settings.Pamb * 1e-5; +switch settings.interpMthd + case 1 + interpFun = 'interp1'; + case 3 + interpFun = 'spline'; + otherwise + error(['\nThe chosen interpolation method cannot be used in this code.\n ...' ... + 'Look at configTestSimulator.m for the allowed options']) +end + +%%% CONSTANTS +g0 = settings.g0; + + +%% COMBUSTION SIMULATION +% Defining the time +engineSim.time = engine.acquisition.cutData.timePTop - engine.acquisition.cutData.timePTop(1); +engineSim.mfr.mox = engine.performances.mox; + +% Performing Marxman simulation +[solution] = ode45(@(t, Ap) marxmanRegressionOde(t, Ap, engine, engineSim, interpFun), tSpan, Ap0); + +% Evaluating solution derivative +[Ap, dAp] = deval(tSpan, solution); +engineSim.cc.Ap = Ap'; +engineSim.cc.dAp = dAp'; +clear solution + +% Converting into radial coordinate +[~, dr, ~] = Ap2RegressionCoor(Ap, dAp, engine); + + +%% CHAMBER MASS BALANCE & PERFORMANCES +% Computing fuel mass flow rate +engineSim.mfr.mf = L*dAp'*rhoFu; + +% Computing oxidizer-to-fuel ratio +OF = engineSim.mfr.mox./engineSim.mfr.mf; +OF(1) = OF(2); +engineSim.performances.OF = OF; + +% Integration of mass balance equation to retrieve pressure +[solution] = ode45(@(t, Pc) massBalance(t, Pc, engine, engineSim, settings, interpFun), tSpan, Pc0); + +% Evaluating solution derivative +[Pc, ~] = deval(tSpan, solution); + +% Conversion from Pa to bar +Pc = Pc'*1e-5; +clear solution + + +%% PERFORMANCES +% Gamma +y = interp2(engine.CEA_db.OF_input, engine.CEA_db.P_input, engine.CEA_db.gamma', OF, Pc, 'spline'); + +% Exit pressure +PeFun = @(Pe, Pc, y) -1./epsilon + ((y + 1)./2).^(1./(y - 1)) .* ((Pe/Pc).^(1./y)) .* sqrt((y + 1)./(y - 1).*(1 - (Pe./Pc).^((y - 1)./y))); +Pe = []; +for ii = 1:length(Pc) + jj = Pc(ii); + kk = y(ii); + Pe = [Pe fsolve(@(Pe) PeFun(Pe, jj, kk), Pamb, fsolveOptions)]; +end + +% Turning the array +Pe = Pe'; + +% Thrust coefficient +cT = sqrt(2*y.^2./(y - 1).*(2./(y + 1)).^((y + 1)./(y - 1))).*sqrt(1 - (Pe./Pc).^((y - 1)./y)) * lambda + (Pe - Pamb)./Pc.*epsilon; + +% Thrust +T = cT.*Pc*1e5*At*etacT; + +% Specific impulse +Isp = T./engineSim.mfr.mox/g0; + +% Total impulse +Itot = trapz(engineSim.time, T); + + +%% OUTPUT +engineSim.cc.dr = dr; +engineSim.cc.Pc = Pc; +engineSim.performances.Pe = Pe; +engineSim.performances.T = T; +engineSim.performances.cT = cT; +engineSim.performances.Isp = Isp; +engineSim.performances.Itot = Itot; +engineSim.mfr.mp = engineSim.mfr.mox + engineSim.mfr.mf; + + diff --git a/RoccarasoFlight/postprocessing/RPS/src/dataLoader.m b/RoccarasoFlight/postprocessing/RPS/src/dataLoader.m index e884a7dd4414bcc1dd3180d3f5ab78be24f728a8..253274bc6057b1794967114771f7c2cc552b8f9e 100644 --- a/RoccarasoFlight/postprocessing/RPS/src/dataLoader.m +++ b/RoccarasoFlight/postprocessing/RPS/src/dataLoader.m @@ -53,7 +53,7 @@ if strcmp(settings.testName, 'ROC1') elseif strcmp(settings.testName, 'ROC2') % Defining the burning time - tb = 3.3; + tb = 4.0; % Defining the first index ip = 8481 * settings.frPT; @@ -75,6 +75,17 @@ elseif strcmp(settings.testName, 'ROC2') flightData.PEng = rawData.PEng(ip:jp); flightData.srvMain = rawData.srvMain(is:js); flightData.srvVent = rawData.srvVent(is:js); + + + % Defining manual cutting indeces + [~, ic] = min(abs(0.980 - flightData.timePTop)); + [~, jc] = min(abs(0.980 + tb - flightData.timePTop)); + + flightData.cutData.timePTop = flightData.timePTop(ic:jc); + flightData.cutData.timePEng = flightData.timePEng(ic:jc); + + flightData.cutData.PTankTop = flightData.PTankTop(ic:jc); + flightData.cutData.PEng = flightData.PEng(ic:jc); end % Fixing pressures diff --git a/RoccarasoFlight/postprocessing/RPS/src/marxmanRegressionOde.m b/RoccarasoFlight/postprocessing/RPS/src/marxmanRegressionOde.m new file mode 100644 index 0000000000000000000000000000000000000000..e39c70dcf2181a2c8460f30cc181f439448c20a3 --- /dev/null +++ b/RoccarasoFlight/postprocessing/RPS/src/marxmanRegressionOde.m @@ -0,0 +1,11 @@ +function dAp = marxmanRegressionOde(t, Ap, engine, engineSim, interpFun) + +%%% INPUT +a = engine.data.fu.a; +n = engine.data.fu.n; + +%%% INTERPOLATION +mox = eval(strcat(interpFun, '(engineSim.time, engineSim.mfr.mox, t)')); + +% Port area derivative +dAp = 2*a*sqrt(pi.*Ap).*(mox./Ap).^n; \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/RPS/src/massBalance.m b/RoccarasoFlight/postprocessing/RPS/src/massBalance.m new file mode 100644 index 0000000000000000000000000000000000000000..cccae0b178926426dd03cfc3d989501114729b20 --- /dev/null +++ b/RoccarasoFlight/postprocessing/RPS/src/massBalance.m @@ -0,0 +1,46 @@ +function dPc = massBalance(t, Pc, engine, engineSim, settings, interpFun) + +%% INPUT + +%%% GRAIN +Vpre = engine.data.Vpre; +Vpost = engine.data.Vpost; +L = engine.data.Lgrain; +etacS = engine.data.etaCs; + +%%% NOZZLE +At = pi * engine.data.rt^2; + +%%% CONSTANTS +R = settings.R; + + +%% INTERPOLATION +mox = eval(strcat(interpFun, '(engineSim.time, engineSim.mfr.mox, t)')); +mf = eval(strcat(interpFun, '(engineSim.time, engineSim.mfr.mf, t)')); +dAp = eval(strcat(interpFun, '(engineSim.time, engineSim.cc.dAp, t)')); +Ap = eval(strcat(interpFun, '(engineSim.time, engineSim.cc.Ap, t)')); + +if t == 0 + OF = engineSim.performances.OF(2); +else + OF = mox/mf; +end + +% CEA thermodynamics parameters +% CC temperature +Tc = interp2(engine.CEA_db.OF_input, engine.CEA_db.P_input, engine.CEA_db.T_comb', OF, Pc*1e-5, 'spline'); + +% Gamma +y = interp2(engine.CEA_db.OF_input, engine.CEA_db.P_input, engine.CEA_db.gamma', OF, Pc*1e-5, 'spline'); + +% Molecular mass +mwc = interp2(engine.CEA_db.OF_input, engine.CEA_db.P_input, engine.CEA_db.mw_comb', OF, Pc*1e-5, 'spline'); + +% Characteristic velocity +cS = sqrt(R./mwc.*Tc)./(sqrt(y .* (2./(y + 1)).^((y + 1)./(y - 1)))); + + +%% EQUATION +% chamber pressure +dPc = (mox + mf - Pc*At/etacS/cS - Pc*dAp*L/R/mwc/Tc/etacS^2) * (R/mwc*Tc*etacS^2/(Ap*L+ Vpre + Vpost)); \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/RPS/src/performanceReconstruction.m b/RoccarasoFlight/postprocessing/RPS/src/performanceReconstruction.m new file mode 100644 index 0000000000000000000000000000000000000000..3d1ceffe7a85824028d66424a96855dc00091642 --- /dev/null +++ b/RoccarasoFlight/postprocessing/RPS/src/performanceReconstruction.m @@ -0,0 +1,52 @@ +function engineFly = performanceReconstruction(engine, engineSim) +%{ +combustionCompare - this function simulates combustion chamber behavior +using estimated/measured mox curve and compares it with those resulting +from combustion with most reliable a, n coefficients retrieved until now + +CALLED SCRIPTS: / + +CALLED FUNCTIONS: / + +CALLED DATA FILES: / + +REVISIONS: +-1, 24/09/2023, Release, Paolo Gnata, PRP Department + +Copyright © 2023, Skyward Experimental Rocketry, PRP Department +All rights reserved + +SPDX-License-Identifier: GPL-3.0-or-later + +%} + +%% INPUT + +%%% DATA +time = engine.acquisition.cutData.timePEng - engine.acquisition.cutData.timePEng(1); +Pc = engine.acquisition.cutData.PEng; + +%%% PARAMETERS +cT = engineSim.performances.cT; +etaCT = engine.data.etacT; +lambda = engine.data.lambda; +At = pi * engine.data.rt^2; + + +%% RESHAPING +% Interpolating the thrust coefficient +cT = interp1(engineSim.time, cT, time); + + +%% PERFORMANCES +% Computing the thrust +T = Pc .* cT * At * etaCT * lambda * 1e5; + +% Computing the total impulse +Itot = trapz(time, T); + + +%% OUTPUT +engineFly.time = time; +engineFly.T = T; +engineFly.Itot = Itot; \ No newline at end of file diff --git a/RoccarasoFlight/postprocessing/RPS/src/tankEvacuation.m b/RoccarasoFlight/postprocessing/RPS/src/tankEvacuation.m index 8e90e499b58e0fb049d3800c89be0a395b724464..89e9769dacf413a1737cc0b58d072067c8a0c694 100644 --- a/RoccarasoFlight/postprocessing/RPS/src/tankEvacuation.m +++ b/RoccarasoFlight/postprocessing/RPS/src/tankEvacuation.m @@ -1,4 +1,4 @@ -function engine = tankEvacuation(engine, settings) +function engine = tankEvacuation(engine) %{ tankEvacuation - this function outputs oxydizer mass flow rate usinìg NHNE model @@ -23,11 +23,12 @@ SPDX-License-Identifier: GPL-3.0-or-later d_pipe = engine.data.d_pipe; N2O = engine.data.ox; -pTank = engine.acquisition.PTankTop; +pTank = engine.acquisition.cutData.PTankTop; Aox = engine.data.Aox; cd = engine.data.cd_inj; -alpha = 1; -time = engine.acquisition.timePTop; +alphaFit = engine.data.alphaFit; + +alphaFun = @(x) sum(alphaFit .* x.^(engine.data.alphaDeg:-1:0), 2); %% CALCULATIONS @@ -40,7 +41,9 @@ v = 10; for i = 1 : length(p1) rho1(i) = interp1(N2O.P, N2O.rhoL, p1(i)); - p1(i) = pTank(i) - (alpha/(2*Aox.^2*rho1(i)).*mox(i).^2)*1e-5; + if i > 1 + p1(i) = pTank(i) - alphaFun(mox(i - 1)) .* (0.5*(mox(i - 1)/Aox).^2)./rho1(i) * 1e-5; + end s1(i) = interp1(N2O.P, N2O.sL, p1(i)); h1(i) = interp1(N2O.P, N2O.hL, p1(i)); p2(i) = p1(i) * N2O.eta_cr(p1(i)); @@ -55,15 +58,15 @@ for i = 1 : length(p1) rhoV2(i) = interp1(N2O.P, N2O.rhoV, p2(i)); rho2(i) = 1/(x2(i) * 1/rhoV2(i) + (1-x2(i))*1/rhoL2(i)); k(i) = sqrt((p1(i)-p2(i))./(pTank(i)-p2(i))); - mox(i) = Aox * cd * (k(i)*sqrt(2*abs(p1(i)-p2(i))*1e5*rho1(i)) + ... + mox(i + 1) = Aox * cd * (k(i)*sqrt(2*abs(p1(i)-p2(i))*1e5*rho1(i)) + ... rho2(i)*sqrt(2*abs(h1(i)-h2(i))))/(1+k(i)); - v(i+1) = mox(i)./rho1(i)/(pi*d_pipe^2/4); + v(i + 1) = mox(i)./rho1(i)/(pi*d_pipe^2/4); end %% OUTPUT % ozydizer mass flow rate -engine.performances.mox = mox; +engine.performances.mox = mox(1:end - 1);