diff --git a/RoccarasoFlight/postprocessing/AFD/descent_para_plot.mat b/RoccarasoFlight/postprocessing/AFD/descent_para_plot.mat
new file mode 100644
index 0000000000000000000000000000000000000000..8be38e73c75f0ab0db3f7ff80ff1c408ec242cab
Binary files /dev/null and b/RoccarasoFlight/postprocessing/AFD/descent_para_plot.mat differ
diff --git a/RoccarasoFlight/postprocessing/AFD/estimationCA.m b/RoccarasoFlight/postprocessing/AFD/estimationCA.m
index d7f08131d7d65be148bd90a8ed5d58513e4c6e31..6f8df89b604f07b9498a9b3ff32eb1a3a1c5099c 100644
--- a/RoccarasoFlight/postprocessing/AFD/estimationCA.m
+++ b/RoccarasoFlight/postprocessing/AFD/estimationCA.m
@@ -4,11 +4,6 @@ function results = estimationCA(main, alt0)
NPOINTS = 1000;
MSTRUCT = 25.2885; % [kg] Gemini structural mass
AREF = 0.15^2 * pi * 0.25;
- PROP_FINAL_MASS = 0.7799 + 0.4445; % [kg] initial propellant mass
- tIGN = 0.97500; % Ignition time in the thrust curve
- tCO = 5.32702; % Cutoff time in the thrsut curve
- mDot = 0.9660; % [kg/s] Averaged propellant mass flow rate
- DELTA_TIME = 0.4988;
load("CA_alpha0.mat");
% acceleration raw
@@ -18,49 +13,25 @@ function results = estimationCA(main, alt0)
accXB = movmean(accXB_raw, 30);
% Motor data
- load("../RPS/cleanData/ROC2/engineFlightData.mat");
- timeThrust = flightData.timeThrust;
- indexTimeThrust = and(timeThrust > tIGN, timeThrust < tCO);
- timeThrust = timeThrust(indexTimeThrust);
- timeThrust = timeThrust - timeThrust(1) - DELTA_TIME;
-
- Thrust = flightData.Thrust(indexTimeThrust);
-
- tNodal = [timeThrust(1), -0.2248, 3.33322, timeThrust(end)];
- indexTransIGN = and(timeThrust>=tNodal(1), timeThrust<tNodal(2));
- indexTransCO = and(timeThrust>tNodal(3), timeThrust<=tNodal(4));
-
- thrustTransIGN = Thrust(indexTransIGN);
- thrustTransCO = Thrust(indexTransCO);
-
- offsetIGN = thrustTransIGN(1);
- offsetCO = thrustTransCO(end);
-
- thrustTransIGN = thrustTransIGN - offsetIGN;
- thrustTransCO = thrustTransCO - offsetCO;
-
- thrustTransIGN = thrustTransIGN * (1 + offsetIGN/thrustTransIGN(end));
- thrustTransCO = thrustTransCO * (1 + offsetCO/thrustTransCO(1));
-
- Thrust = [thrustTransIGN; Thrust(and(timeThrust>=tNodal(2), timeThrust<=tNodal(3))); thrustTransCO];
- mMotor = linspace(PROP_FINAL_MASS + mDot*(tCO - tIGN), PROP_FINAL_MASS, sum(indexTimeThrust));
- timeThrust = timeThrust - timeThrust(1);
-
+ motor = load("roccEngine.mat");
+ timeThrust = motor.time;
+ Thrust = motor.Thrust;
+ mMotor = motor.mass;
+
% Velocity
- tNas = main.NASData(1:main.NASApogeeIndex, 1);
- altNas = -main.NASData(1:main.NASApogeeIndex, 4);
- velNas = vecnorm(main.NASData(1:main.NASApogeeIndex, 5:7), 2, 2);
-
- load sensorPitot.mat
- indNasPit = sensorTot_WithPitot.nas.time < main.tApogee;
+ tNasFlight = main.NASData(1:main.NASApogeeIndex, 1);
+ altNasFlight = -main.NASData(1:main.NASApogeeIndex, 4);
+ velNasFlight = vecnorm(main.NASData(1:main.NASApogeeIndex, 5:7), 2, 2);
- tNasPit = sensorTot_WithPitot.nas.time(indNasPit);
- velNasPit = vecnorm(sensorTot_WithPitot.nas.states(indNasPit, 4:6), 2,2);
- altNasPit = -sensorTot_WithPitot.nas.states(indNasPit, 3);
+ load nasSim.mat
+ indNasSim = sensorTot_NoPitot.nas.time < main.tApogee;
+ tNasSim = sensorTot_NoPitot.nas.time(indNasSim);
+ velNasSim = vecnorm(sensorTot_NoPitot.nas.states(indNasSim, 4:6), 2,2);
+ altNasSim = -sensorTot_NoPitot.nas.states(indNasSim, 3);
% Normalize data
- timeRef = linspace(max([tNas(1), tAcc(1), tNasPit(1)]), min([tNas(end), tAcc(end), tNasPit(end)]) , NPOINTS);
+ timeRef = linspace(max([tNasFlight(1), tAcc(1), tNasSim(1)]), min([tNasFlight(end), tAcc(end), tNasSim(end)]) , NPOINTS);
indexEndThrust = find(timeRef > timeThrust(end), 1, 'first') - 1;
accXBRef = interp1(tAcc, accXB, timeRef);
@@ -71,68 +42,67 @@ function results = estimationCA(main, alt0)
mRocket = (mMotor(end) + MSTRUCT) * ones(1, NPOINTS);
mRocket(1:indexEndThrust) = interp1(timeThrust, mMotor, timeRef(1:indexEndThrust)) + MSTRUCT;
- velRef = interp1(tNas, movmean(velNas, 10), timeRef);
+ velNasFlightNorm = interp1(tNasFlight, movmean(velNasFlight, 10), timeRef);
+ velNasSimNorm = interp1(tNasSim, velNasSim, timeRef);
- velNasPitRef = interp1(tNasPit, velNasPit, timeRef);
-
- altRef = interp1(tNasPit, altNasPit, timeRef);
+ velVertNasFlightNorm = interp1(tNasFlight, movmean(-main.NASData(1:main.NASApogeeIndex, 7), 10), timeRef);
+ velVertNasSimNorm = interp1(tNasSim, -sensorTot_NoPitot.nas.states(indNasSim, 6), timeRef);
+
+ altFlightNorm = interp1(tNasFlight, altNasFlight, timeRef);
+ altSimNorm = interp1(tNasSim, altNasSim, timeRef);
% Air density and sonic velocity
- [~, vSonRef, ~, rhoRef] = atmosisa(altRef + alt0);
- machRef = velRef./vSonRef;
- machPitRef = velNasPitRef./vSonRef;
+ [~, vSonSimNorm, ~, rhoSimNorm] = atmosisa(altSimNorm + alt0);
+ [~, vSonFlightNorm, ~, rhoFlightNorm] = atmosisa(altFlightNorm + alt0);
+
+ machFlightNorm = velNasFlightNorm./vSonFlightNorm;
+ machSimNorm = velNasSimNorm./vSonSimNorm;
% index opening
indexOpeningARB = find(timeRef > 5.12, 1, 'first');
% Compute CA
- CA = (ThrustRef - mRocket.*accXBRef) ./ (0.5 .* rhoRef .* velRef.^2 * AREF);
+ CA_NasFlight = (ThrustRef - mRocket.*accXBRef) ./ (0.5 .* rhoFlightNorm .* velNasFlightNorm.^2 * AREF);
- CA_NasPitot = (ThrustRef - mRocket.*accXBRef) ./ (0.5 .* rhoRef .* velNasPitRef.^2 * AREF);
+ CA_NasSim = (ThrustRef - mRocket.*accXBRef) ./ (0.5 .* rhoSimNorm .* velNasSimNorm.^2 * AREF);
- tCO = tCO -1.8519;
indexCO = find(ThrustRef == 0, 1, 'first');
- mRocketNoPwr = mRocket(indexCO:end);
- CANoPwr = CA(indexCO:end);
CA_100ARB = zeros(1, length(timeRef));
CA_0ARB = zeros(1, length(timeRef));
for i = 1:length(timeRef)
- [~, indexMach] = nearestValSorted(MachRef, machPitRef(i));
- [~, indexAlt] = nearestValSorted(AltRef, altRef(i));
+ [~, indexMach] = nearestValSorted(MachRef, machSimNorm(i));
+ [~, indexAlt] = nearestValSorted(AltRef, altSimNorm(i));
CA_100ARB(i) = CA_MD(indexMach, indexAlt, 3);
CA_0ARB(i) = CA_MD(indexMach, indexAlt, 1);
end
- ThrustComp = 0.5 .* rhoRef(1:indexCO) .* velRef(1:indexCO).^2 .* AREF .* CA_0ARB(1:indexCO) + mRocket(1:indexCO).*accXBRef(1:indexCO);
+ ThrustComp = 0.5 .* rhoSimNorm(1:indexCO) .* velNasFlightNorm(1:indexCO).^2 .* AREF .* CA_0ARB(1:indexCO) + mRocket(1:indexCO).*accXBRef(1:indexCO);
- figure
- plot(timeRef(indexCO:end),CANoPwr);
- hold on
- plot(timeRef(indexCO:end), CA_NasPitot(indexCO:end));
- plot(timeRef(indexCO:end), [CA_0ARB(indexCO:indexOpeningARB) CA_100ARB(indexOpeningARB+1:end)]);
- xlabel('Time [s]');
- ylabel('CA [-]');
- title('CA from telemetry');
- grid on
- legend('NAS', 'NAS with PITOT', 'MD', 'Location','southoutside', 'Orientation', 'horizontal');
-
%% SAVE RESULTS
results.time = timeRef;
results.Thrust = ThrustRef;
results.ThrustComputed = ThrustComp;
- results.velocity.vel_pitot = velNasPitRef;
- results.velocity.val_noPitot = velRef;
-
- results.alt = altRef;
+ results.velocity.vel_NasSim = velNasSimNorm;
+ results.velocity.vel_NasFlight = velNasFlightNorm;
+ results.velocity.velVert_NasSim = velVertNasSimNorm;
+ results.velocity.velVert_NasFlight = velVertNasFlightNorm;
+
+ results.altitude.alt_NasSim = altSimNorm;
+ results.altitude.alt_NasFlight = altFlightNorm;
+
+ results.mach.mach_NasSim = machSimNorm;
+ results.mach.mach_NasFlight = machFlightNorm;
+
+
results.CA.CA100 = CA_100ARB;
results.CA.CA0 = CA_0ARB;
- results.CA.CA_pitot = CA_NasPitot;
- results.CA.CA_noPitot = CA;
+ results.CA.CA_NasSim = CA_NasSim;
+ results.CA.CA_NasFlight = CA_NasFlight;
results.mass.propMass = mMotor;
results.mass.rocketMass = mRocket;
diff --git a/RoccarasoFlight/postprocessing/AFD/mainPostprocess.m b/RoccarasoFlight/postprocessing/AFD/mainPostprocess.m
index 3616d5a4580deab3333a7d604fcc6a5dae1556b5..c7b1fde450c038776b8e966e81df3b186d692a5f 100644
--- a/RoccarasoFlight/postprocessing/AFD/mainPostprocess.m
+++ b/RoccarasoFlight/postprocessing/AFD/mainPostprocess.m
@@ -3,13 +3,14 @@ clc;
close all;
%% Select options
-opt.flagTraj = false;
+opt.flagTraj = true;
opt.flagCA = true;
alt0 = 1414;
%% adding path and recall data
addpath("..\commonFunctions\");
-addpath(genpath("..\..\..\..\msa-toolkit\commonFunctions\"));
+addpath(genpath("..\..\..\..\msa-toolkit\commonFunctions\"));
+addpath(genpath("..\..\..\..\msa-toolkit\simulator\src\"))
main = importData('main');
%% Trajectory
@@ -21,18 +22,18 @@ end
%% CA Estimation
if opt.flagCA
+ clc
res = estimationCA(main,alt0);
-
+ indexCO = res.events.shutdonw;
+ indexOpeningARB = res.events.ARB_opening;
end
%% TEST Simulation
if opt.flagCA
-
- indexCO = res.events.shutdonw;
- indexOpeningARB = res.events.ARB_opening;
configSimPostp;
+ settings.Local = [1414, 25 + 273, 88369];
% CONTINUE
Thrust = res.Thrust(1:indexCO);
@@ -60,8 +61,61 @@ if opt.flagCA
settings.motor.expThrust = Thrust;
settings.tb = settings.motor.expTime(end);
settings.tControl = settings.tb;
+
+ settings.z0 = settings.z0 - 100;
[Tf, Yf, Ta, Ya, bound_value] = stdRun(settings);
+ load ascent_plot.mat;
+
+ %%
+ close all;
+
+ figure
+ plot(res.time, res.altitude.alt_NasFlight);
+ hold on
+ plot(res.time, res.altitude.alt_NasSim);
+ plot(Ta, -Ya(:, 3))
+ xlabel('Time [s]');
+ ylabel('Altitude [m]');
+ title('Altitude');
+ grid on
+ legend('NAS - flight', 'NAS - simulated', 'Simulation', 'Location','southoutside', 'Orientation', 'horizontal');
+
+ figure
+ plot(res.time, res.velocity.vel_NasFlight);
+ hold on
+ plot(res.time, res.velocity.vel_NasSim);
+ plot(Ta, vecnorm(Ya(:, 4:6), 2, 2));
+ xlabel('Time [s]');
+ ylabel('Velocity [m/s]');
+ title('Total Velocity');
+ grid on
+ legend('NAS - flight', 'NAS - simulated', 'Simulation', 'Location','southoutside', 'Orientation', 'horizontal');
+
+ figure
+ plot(res.time, res.velocity.velVert_NasFlight);
+ hold on
+ plot(res.time, res.velocity.velVert_NasSim);
+ plot(Ta, -data_ascent.velocities(3, :));
+ xlabel('Time [s]');
+ ylabel('Velocity [m/s]');
+ title('Vertical Velocity');
+ grid on
+ legend('NAS - flight', 'NAS - simulated', 'Simulation', 'Location','southoutside', 'Orientation', 'horizontal');
+
+
+ figure
+ plot(res.time(indexCO:end),res.CA.CA_NasFlight(indexCO:end));
+ hold on
+ plot(res.time(indexCO:end),res.CA.CA_NasSim(indexCO:end));
+ plot(res.time(indexCO:end), [res.CA.CA0(indexCO:indexOpeningARB) res.CA.CA100(indexOpeningARB+1:end)]);
+ plot(data_ascent.integration.t, data_ascent.coeff.CA);
+ xlabel('Time [s]');
+ ylabel('CA [-]');
+ title('CA from telemetry');
+ grid on
+ legend('NAS - flight', 'NAS - simulated', 'MD', 'Simulation', 'Location','southoutside', 'Orientation', 'horizontal');
+
end
diff --git a/RoccarasoFlight/postprocessing/AFD/nasSim.mat b/RoccarasoFlight/postprocessing/AFD/nasSim.mat
new file mode 100644
index 0000000000000000000000000000000000000000..60c73ce495cf34f198d334dbb3e23d0d9b30c4ab
Binary files /dev/null and b/RoccarasoFlight/postprocessing/AFD/nasSim.mat differ
diff --git a/RoccarasoFlight/postprocessing/AFD/roccEngine.mat b/RoccarasoFlight/postprocessing/AFD/roccEngine.mat
new file mode 100644
index 0000000000000000000000000000000000000000..d5c3231a657a2d4c57092071be9c4047465c4795
Binary files /dev/null and b/RoccarasoFlight/postprocessing/AFD/roccEngine.mat differ
diff --git a/RoccarasoFlight/postprocessing/AFD/src/roccEngine.mat b/RoccarasoFlight/postprocessing/AFD/src/roccEngine.mat
new file mode 100644
index 0000000000000000000000000000000000000000..d5c3231a657a2d4c57092071be9c4047465c4795
Binary files /dev/null and b/RoccarasoFlight/postprocessing/AFD/src/roccEngine.mat differ
diff --git a/RoccarasoFlight/postprocessing/AFD/src/thrustBricolage.m b/RoccarasoFlight/postprocessing/AFD/src/thrustBricolage.m
new file mode 100644
index 0000000000000000000000000000000000000000..59bde01708f551575b9a785f4d942ca1065a4f64
--- /dev/null
+++ b/RoccarasoFlight/postprocessing/AFD/src/thrustBricolage.m
@@ -0,0 +1,56 @@
+clear;
+clc;
+close all;
+
+%% Load data
+
+load("..\..\RPS\msaEngineData.mat"); % curva nuova
+load("..\..\RPS\cleanData\ROC2\engineFlightDataOld.mat"); % curva originale
+addpath("..\..\..\..\..\msa-toolkit\commonFunctions\miscellaneous\");
+
+%% Isolate transients
+tNodal = [0.982008, 1.25, 4.79702, 5.10202];
+timeOld = flightData.timeThrust;
+ThrustOld = flightData.Thrust;
+
+indIGN = and(timeOld > tNodal(1), timeOld < tNodal(2));
+indCO = and(timeOld > tNodal(3), timeOld < tNodal(4));
+
+timeIGN = timeOld(indIGN);
+ThrustIGN = ThrustOld(indIGN);
+
+timeCO = timeOld(indCO);
+ThrustCO = ThrustOld(indCO);
+
+%% Add transient to computed thrust
+
+timeComp = msaEngineData.timeThrustReal;
+ThrustComp = msaEngineData.ThrustReal;
+
+timeIGN = timeIGN - timeIGN(1);
+timeComp = timeComp - timeComp(1) + timeIGN(end);
+timeCO = timeCO - timeCO(1) + timeComp(end);
+
+%% Reshape transient
+
+T0 = ThrustComp(1);
+TEnd = ThrustComp(end);
+
+ThrustIGN = ThrustIGN - ThrustIGN(1); % parte da 0
+ThrustCO = ThrustCO - ThrustCO(end); % finisce a 0
+
+ThrustIGN = ThrustIGN * T0/ThrustIGN(end);
+ThrustCO = ThrustCO * TEnd/ThrustCO(1);
+
+%% Plot results
+
+time = [timeIGN(1:end-1); timeComp; timeCO(2:end)];
+Thrust = [ThrustIGN(1:end-1); ThrustComp; ThrustCO(2:end)];
+mass = reshapeVector(time, msaEngineData.time, msaEngineData.Mp);
+
+plot(time, Thrust);
+
+save('roccEngine.mat', 'time', 'Thrust', 'mass');
+
+
+
diff --git a/RoccarasoFlight/postprocessing/GNC/postProcess.m b/RoccarasoFlight/postprocessing/GNC/postProcess.m
index 88f9b13b84944283fd5da13246e8a6e1604dac91..373d4d6be0bce280c9728c0e394570935d4156e0 100644
--- a/RoccarasoFlight/postprocessing/GNC/postProcess.m
+++ b/RoccarasoFlight/postprocessing/GNC/postProcess.m
@@ -505,6 +505,8 @@ title('Pressure')
legend
xlabel('Time [s]')
ylabel('Pressure [Pa]')
+xline(motor.mainServoData(find(motor.mainServoData(:,2)>0.965,1),1),'k--','OPEN MAIN','HandleVisibility','off')
+xline(motor.mainServoData(find(motor.mainServoData(:,2)<0.001,1),1),'k--','CLOSE MAIN','HandleVisibility','off')
% apogee
subplot(3,1,3)
plot(main.MEAControllerStatus(main.MEAControllerStatus(:,1)<t_end,1),main.MEAControllerStatus(main.MEAControllerStatus(:,1)<t_end,3),'DisplayName','OBSW')
@@ -514,6 +516,8 @@ plot(sensorTot_NoPitot.mea.time,sensorTot_NoPitot.mea.predictedApogee,'DisplayNa
plot(sensorTot_WithPitot.mea.time,sensorTot_WithPitot.mea.predictedApogee,'DisplayName','Simulator WithPitot' )
plot(main.NASData(main.NASData(:,1)<t_end,1),-main.NASData(main.NASData(:,1)<t_end,4),'DisplayName','NAS Altitude')
yline(main.MEATarget,'--','SHUTDOWN TARGET','HandleVisibility','off')
+xline(motor.mainServoData(find(motor.mainServoData(:,2)>0.965,1),1),'k--','OPEN MAIN','HandleVisibility','off')
+xline(motor.mainServoData(find(motor.mainServoData(:,2)<0.001,1),1),'k--','CLOSE MAIN','HandleVisibility','off')
legend
title('Apogee prediction')
xlabel('Time [s]')
diff --git a/RoccarasoFlight/postprocessing/RPS/cleanData/ROC2/engineFlightData.mat b/RoccarasoFlight/postprocessing/RPS/cleanData/ROC2/engineFlightData.mat
index aef81bab5ef6ee0f4df703b96c10110556d95a46..580fd2205e9eaad10b14e2b4a6085a66b9a27bdc 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 441dcc97598c2c3f464aa411ec2b11bf716684e7..d8ff216dafbd62b3d4b5de2232a2d0ad686884f6 100644
--- a/RoccarasoFlight/postprocessing/RPS/configDataAnalyzer.m
+++ b/RoccarasoFlight/postprocessing/RPS/configDataAnalyzer.m
@@ -22,6 +22,7 @@ clear CEA_db
%% SIMULATION SETTINGS
settings.interpMthd = 1;
+settings.saveResults = 1;
%% CONSTANTS
settings.g0 = 9.806; % [m/s2] Acceleration of Gravity, constant
diff --git a/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m b/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m
index ca73d482882adb77bf26974dfaa359276ef2373b..c86f3356f311c2a972f5860ec61b769af0cbc109 100644
--- a/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m
+++ b/RoccarasoFlight/postprocessing/RPS/dataAnalyzer.m
@@ -36,14 +36,31 @@ engine = tankEvacuation(engine);
% Simulating the engine
engineSim = combustionSimulation(engine, settings);
-% Computing actual performance
+% Computing actual performances
engineFly = performanceReconstruction(engine, engineSim);
+% Loading MSA simulated motor for the flight
+load("HREmotors.mat");
+engineMSA = HREmotors(59);
+clear HREmotors;
+
%% 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.saveResults
+ msaEngineData.Mp = engineSim.mass.Mp;
+ msaEngineData.Mox = engineSim.mass.Mox;
+ msaEngineData.Mfu = engineSim.mass.Mfu;
+ msaEngineData.mp = engineSim.mfr.mp;
+ msaEngineData.time = engineSim.time;
+ msaEngineData.Thrust = engineSim.performances.T;
+ msaEngineData.ThrustReal = engineFly.T;
+ msaEngineData.timeThrustReal = engineFly.time;
+ save('msaEngineData', 'msaEngineData');
+end
+
if settings.plots
f1 = figure;
yyaxis left
@@ -68,8 +85,9 @@ if settings.plots
plot(engineSim.time, engineSim.performances.T)
hold on; grid minor
plot(engineFly.time, engineFly.T)
+ plot(engineMSA.t, engineMSA.T)
ylabel('Thrust [N]'); xlabel('Time [s]'); title('Thrust compare')
- legend('Simulated', 'Computed')
+ legend('Simulated', 'Computed', 'MSA-SFT06')
axis padded
end
diff --git a/RoccarasoFlight/postprocessing/RPS/engineData/engines/CEA_db.mat b/RoccarasoFlight/postprocessing/RPS/engineData/engines/CEA_db.mat
index 92de427b9ecb77f7f42a48eb7eda47fdef275435..e1e3861b8edc88853c5be0794f2974b912720d30 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 7190b74c04cf1b46d10ba6f40da6562379955075..343a6d76c4a146f4025f6851f0575307881f0054 100644
--- a/RoccarasoFlight/postprocessing/RPS/engineData/engines/Furia/Furia.m
+++ b/RoccarasoFlight/postprocessing/RPS/engineData/engines/Furia/Furia.m
@@ -58,7 +58,7 @@ 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
+engine.data.Pc0 = 2; % [bar] Initial pressure value
%%% PRE CC
engine.data.Dpre = 66e-3;
@@ -74,3 +74,8 @@ engine.data.Vpost = pi/4*engine.data.Dpost^2*engine.data.Lpost;
engine.data.zopt = 0; % [m] Altitude of optimal expansion
engine.data.RAO_Lratio = 0.8; % [-] Rao length ratio
+
+%% ENGINE TEST PARAMETERS
+engine.postFire.Mox = 0.4445; % [kg] Oxidizer residual mass @ X = 1
+engine.postFire.Mfu = 0.7799; % [kg] Fuel residual mass
+
diff --git a/RoccarasoFlight/postprocessing/RPS/msaData/HREmotors.mat b/RoccarasoFlight/postprocessing/RPS/msaData/HREmotors.mat
new file mode 100644
index 0000000000000000000000000000000000000000..b7b3c10914d77aee38c5c66ed75e6f72bdf9926f
Binary files /dev/null and b/RoccarasoFlight/postprocessing/RPS/msaData/HREmotors.mat differ
diff --git a/RoccarasoFlight/postprocessing/RPS/msaEngineData.mat b/RoccarasoFlight/postprocessing/RPS/msaEngineData.mat
new file mode 100644
index 0000000000000000000000000000000000000000..7ad64beb50dd65246eaf943eb67c0b7183694f8f
Binary files /dev/null and b/RoccarasoFlight/postprocessing/RPS/msaEngineData.mat differ
diff --git a/RoccarasoFlight/postprocessing/RPS/src/combustionSimulation.m b/RoccarasoFlight/postprocessing/RPS/src/combustionSimulation.m
index 19cad597d0dd54d74ed290e02fa2abc746045092..1bfaeaf8c5f67537f3d6fdc6068dfdb133e9878b 100644
--- a/RoccarasoFlight/postprocessing/RPS/src/combustionSimulation.m
+++ b/RoccarasoFlight/postprocessing/RPS/src/combustionSimulation.m
@@ -59,6 +59,8 @@ g0 = settings.g0;
%% COMBUSTION SIMULATION
% Defining the time
engineSim.time = engine.acquisition.cutData.timePTop - engine.acquisition.cutData.timePTop(1);
+
+% Saving the oxidizer mass flow in the out-struct
engineSim.mfr.mox = engine.performances.mox;
% Performing Marxman simulation
@@ -83,8 +85,12 @@ OF = engineSim.mfr.mox./engineSim.mfr.mf;
OF(1) = OF(2);
engineSim.performances.OF = OF;
+% Fixing convergence problems in the first transitory phase (performances are not relevant in this phase)
+ii = find(OF < 2);
+OF(ii) = 2;
+
% Integration of mass balance equation to retrieve pressure
-[solution] = ode45(@(t, Pc) massBalance(t, Pc, engine, engineSim, settings, interpFun), tSpan, Pc0);
+[solution] = ode23s(@(t, Pc) massBalance(t, Pc, engine, engineSim, settings, interpFun), tSpan, Pc0);
% Evaluating solution derivative
[Pc, ~] = deval(tSpan, solution);
@@ -108,10 +114,10 @@ for ii = 1:length(Pc)
end
% Turning the array
-Pe = Pe';
+Pe = real(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;
+cT = sqrt(2*y.^2./(y - 1).*(2./(y + 1)).^((y + 1)./(y - 1))).*sqrt(1 - (Pe./Pc).^((y - 1)./y)) * lambda;
% Thrust
T = cT.*Pc*1e5*At*etacT;
@@ -122,6 +128,19 @@ Isp = T./engineSim.mfr.mox/g0;
% Total impulse
Itot = trapz(engineSim.time, T);
+% Computing the oxidizer mass evolution
+Mox = @(t, y) interp1(engineSim.time, engineSim.mfr.mox, t);
+[~, Mox] = ode45(Mox, engineSim.time, engine.postFire.Mox);
+Mox = flip(Mox);
+
+% Computing the fuel mass evolution
+Mfu = @(t, y) interp1(engineSim.time, engineSim.mfr.mf, t);
+[~, Mfu] = ode45(Mfu, engineSim.time, engine.postFire.Mfu);
+Mfu = flip(Mfu);
+
+% Computing the propellant mass evolution
+Mp = Mfu + Mox;
+
%% OUTPUT
engineSim.cc.dr = dr;
@@ -132,5 +151,8 @@ engineSim.performances.cT = cT;
engineSim.performances.Isp = Isp;
engineSim.performances.Itot = Itot;
engineSim.mfr.mp = engineSim.mfr.mox + engineSim.mfr.mf;
+engineSim.mass.Mox = Mox;
+engineSim.mass.Mfu = Mfu;
+engineSim.mass.Mp = Mp;
diff --git a/RoccarasoFlight/postprocessing/RPS/src/dataLoader.m b/RoccarasoFlight/postprocessing/RPS/src/dataLoader.m
index 253274bc6057b1794967114771f7c2cc552b8f9e..e84abd52a4b7f233f30f1a538ef4b17990015490 100644
--- a/RoccarasoFlight/postprocessing/RPS/src/dataLoader.m
+++ b/RoccarasoFlight/postprocessing/RPS/src/dataLoader.m
@@ -52,8 +52,8 @@ if strcmp(settings.testName, 'ROC1')
flightData.srvVent = rawData.srvVent(is:js);
elseif strcmp(settings.testName, 'ROC2')
- % Defining the burning time
- tb = 4.0;
+ % Defining the burning time - without ending transitory
+ tb = 3.5;
% Defining the first index
ip = 8481 * settings.frPT;
@@ -78,8 +78,8 @@ elseif strcmp(settings.testName, 'ROC2')
% Defining manual cutting indeces
- [~, ic] = min(abs(0.980 - flightData.timePTop));
- [~, jc] = min(abs(0.980 + tb - flightData.timePTop));
+ [~, ic] = min(abs(1.3 - flightData.timePTop));
+ [~, jc] = min(abs(1.3 + tb - flightData.timePTop));
flightData.cutData.timePTop = flightData.timePTop(ic:jc);
flightData.cutData.timePEng = flightData.timePEng(ic:jc);
diff --git a/RoccarasoFlight/postprocessing/RPS/src/massBalance.m b/RoccarasoFlight/postprocessing/RPS/src/massBalance.m
index cccae0b178926426dd03cfc3d989501114729b20..d5cd381a54c737dc82d0bb0db54e8f7f8678873d 100644
--- a/RoccarasoFlight/postprocessing/RPS/src/massBalance.m
+++ b/RoccarasoFlight/postprocessing/RPS/src/massBalance.m
@@ -21,8 +21,8 @@ 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);
+if mox/mf < 2 || isnan(mox/mf)
+ OF = 2;
else
OF = mox/mf;
end
diff --git a/RoccarasoFlight/postprocessing/RPS/src/tankEvacuation.m b/RoccarasoFlight/postprocessing/RPS/src/tankEvacuation.m
index 89e9769dacf413a1737cc0b58d072067c8a0c694..ff911f16855c7160d2f54a11064c65605452b6bb 100644
--- a/RoccarasoFlight/postprocessing/RPS/src/tankEvacuation.m
+++ b/RoccarasoFlight/postprocessing/RPS/src/tankEvacuation.m
@@ -63,6 +63,14 @@ for i = 1 : length(p1)
v(i + 1) = mox(i)./rho1(i)/(pi*d_pipe^2/4);
end
+% Adding transitory to injection mox
+% tTrans = 0.1;
+% [~, j] = min(abs(engine.acquisition.cutData.timePTop - engine.acquisition.cutData.timePTop(1) - tTrans));
+% moxTrans = linspace(mox(1), mox(2), j)';
+% mox = [moxTrans ; mox(3:end)];
+% time = [linspace(0, engine.acquisition.cutData.timePTop(j) - engine.acquisition.cutData.timePTop(1), j)' ; engine.acquisition.cutData.timePTop(3:end) - 2*engine.acquisition.cutData.timePTop(1) + engine.acquisition.cutData.timePTop(j)];
+% mox = interp1(time, mox(1:end-1), engine.acquisition.cutData.timePTop - engine.acquisition.cutData.timePTop(1));
+
%% OUTPUT
% ozydizer mass flow rate
diff --git a/RoccarasoFlight/postprocessing/commonFunctions/atmosphereAlt.m b/RoccarasoFlight/postprocessing/commonFunctions/atmosphereAlt.m
new file mode 100644
index 0000000000000000000000000000000000000000..52764d6cb71e24875a543be72113a95d559ba391
--- /dev/null
+++ b/RoccarasoFlight/postprocessing/commonFunctions/atmosphereAlt.m
@@ -0,0 +1,35 @@
+function absoluteAltitude = atmosphereAlt(pressure, local)
+
+
+
+atmCostants = [0.0065, 401.87434, 1.86584515, 1.225, ...
+ 101325, 288.15, 287.05]; % atmosisa constants: [a, gamma*R, R*a, rho0, p0, T0, R]
+g = 9.81;
+P0 = atmCostants(5);
+T0 = atmCostants(6);
+
+if length(local) == 4 % if T, P, rho (at local altitude) and local altitude are given as inputs, T0, P0 and rho0 are computed in order to have a more accurate atmosphere model
+ altLocal = local(1);
+ TLocal = local(2);
+ PLocal = local(3);
+ T0 = TLocal + atmCostants(1)*altLocal;
+ theta0 = TLocal/T0;
+ P0 = PLocal/(theta0^(g/atmCostants(3)));
+elseif length(local) == 3 % if T, P (at local altitude) and local altitude are given as inputs, T0, P0 and rho0 are computed in order to have a more accurate atmosphere model
+ altLocal = local(1);
+ TLocal = local(2);
+ PLocal = local(3);
+ T0 = TLocal + atmCostants(1)*altLocal;
+ theta0 = TLocal/T0;
+ P0 = PLocal/(theta0^(g/atmCostants(3)));
+elseif length(local) == 2 % if T (at a local altitude) and local altitude are given as inputs, T0, P0 and rho0 are computed in order to have a more accurate atmosphere model
+ altLocal = local(1);
+ TLocal = local(2);
+ T0 = TLocal + atmCostants(1)*altLocal;
+end
+
+theta = (pressure./P0).^(atmCostants(3)./g);
+T = T0.*theta;
+absoluteAltitude = (T0 - T)./atmCostants(1);
+
+end