diff --git a/classes/bays/Payload.m b/classes/bays/Payload.m index 4b54423f198fbde8563a00e836f102a734d39f74..0f9f7edb3f5528648e3142d3cb85cfe508438e49 100644 --- a/classes/bays/Payload.m +++ b/classes/bays/Payload.m @@ -24,30 +24,7 @@ classdef Payload < Bay length % [m] Total bay length mass % [kg] Total bay mass inertia % [kg*m^2] Total bay inertia (Body reference) - inertiaMatrix % [kg*m^2] 3x3 Inertia Matrix (the diagonal is made of the inertia vector values) xCg % [m] Cg relative to bay upper side - - paraPrev double % Reference to the previous parachute in the multistage descent - semiWingSpan double % [m] semiwingspan - MAC double % [m] mean aero chord - surface double % [m^2] payload surface - - deltaSMax double % [-] aerodynamic control coefficient - symmetric, max value - - CD0 double % [-] aerodynamic control coefficient - asymetric, CD0 - CDAlpha2 double % [-] aerodynamic control coefficient - asymetric, CDAlpha2 - CL0 double % [-] aerodynamic control coefficient - asymetric, CL0 - CLAlpha double % [-] aerodynamic control coefficient - asymetric, CLAlpha - Cm0 double % [-] aerodynamic control coefficient - asymetric, Cm0 - CmAlpha double % [-] aerodynamic control coefficient - asymetric, CmAlpha - Cmq double % [-] aerodynamic control coefficient - asymetric, Cmq - CLDeltaA double % [-] aerodynamic control coefficient - asymetric, CLDeltaA - Cnr double % [-] aerodynamic control coefficient - asymetric, Cnr - CnDeltaA double % [-] aerodynamic control coefficient - asymetric, CnDeltaA - CDDeltaA double % [-] aerodynamic control coefficient - asymetric, CDDeltaA - Clp double % [-] aerodynamic control coefficient - asymetric, Clp - ClPhi double % [-] aerodynamic control coefficient - asymetric, ClPhi - ClDeltaA double % [-] aerodynamic control coefficient - asymetric, ClDeltaA end properties (Dependent) diff --git a/classes/components/Environment.m b/classes/components/Environment.m index 4c6806100c2b9d00047b6260938f8d2b0cefdb89..e9b47cc6aa304805dea5d0c326ec95379bfcbc75 100644 --- a/classes/components/Environment.m +++ b/classes/components/Environment.m @@ -84,6 +84,9 @@ classdef Environment < Component end function obj = updateLocal(obj) + if isempty(obj.temperature), obj.temperature = 288.15; end + if isempty(obj.gamma), obj.gamma = 1.4; end + obj.local = [obj.z0, obj.temperature, ... % vector containing inputs for atmosphereData obj.pressure, obj.rho]; end diff --git a/classes/components/Parafoil.m b/classes/components/Parafoil.m index 364256ba1714c63dbfd3c44f86b720026c371d2b..c79ac4ac1ff56f5804463748f37b1ea68acf7f5a 100644 --- a/classes/components/Parafoil.m +++ b/classes/components/Parafoil.m @@ -11,22 +11,34 @@ classdef Parafoil < Para % file properties - name = '' - surface double % [m^2] Surface - mass double % [kg] Parachute Mass - openingDelay double % [s] drogue opening delay - finalAltitude double % [m] Final altitude of the parachute - chordLength double % [m] Shock Chord Length - chordK double % [N/m^2] Shock Chord Elastic Constant - chordC double % [Ns/m] Shock Chord Dynamic Coefficient - expulsionSpeed double - cx double % [/] Parachute Longitudinal Drag Coefficient - cd double % [/] Parachute Drag Coefficient - cl double % [/] Parachute Lift Coefficient - m double % [m^2/s] Coefficient of the surface vs. time opening model - nf double % [/] Adimensional Opening Time - deployDuration double % [s] Time to get the parachute to full aperture - forceCoefficient double % [-] Empirical coefficient to obtain correct peak force at deployment + name = '' + surface double % [m^2] Surface + deltaSMax double % [-] aerodynamic control coefficient - symmetric, max value + + mass double % [kg] Parachute Mass + + inertia % [kg*m^2] 3x3 Inertia Matrix (the diagonal is made of the inertia vector values) + + semiWingSpan double % [m] semiwingspan + MAC double % [m] mean aero chord + + cd0 double % [-] aerodynamic control coefficient - asymetric, CD0 + cdAlpha2 double % [-] aerodynamic control coefficient - asymetric, CDAlpha2 + cdDeltaA double % [-] aerodynamic control coefficient - asymetric, CDDeltaA + cl0 double % [-] aerodynamic control coefficient - asymetric, CL0 + clAlpha double % [-] aerodynamic control coefficient - asymetric, CLAlpha + clDeltaA double % [-] aerodynamic control coefficient - asymetric, CLDeltaA + clP double % [-] aerodynamic control coefficient - asymetric, Clp + clPhi double % [-] aerodynamic control coefficient - asymetric, ClPhi + cm0 double % [-] aerodynamic control coefficient - asymetric, Cm0 + cmAlpha double % [-] aerodynamic control coefficient - asymetric, CmAlpha + cmQ double % [-] aerodynamic control coefficient - asymetric, Cmq + cnR double % [-] aerodynamic control coefficient - asymetric, Cnr + cnDeltaA double % [-] aerodynamic control coefficient - asymetric, CnDeltaA + end + + properties(SetAccess = private) + inverseInertia end methods @@ -39,5 +51,9 @@ classdef Parafoil < Para obj@Para(mission, varIn); obj = obj(:, :, 1); end + + function updateAll(obj) + obj.inverseInertia = obj.inertia\eye(3); + end end end \ No newline at end of file diff --git a/functions/odeFunctions/descentParafoil.m b/functions/odeFunctions/descentParafoil.m index c6590ad22673c2aa0fa3a91a44f05b5df1a05b72..1a328c55fa82631a0cf1e562d9d60b9f89d9b374 100644 --- a/functions/odeFunctions/descentParafoil.m +++ b/functions/odeFunctions/descentParafoil.m @@ -51,11 +51,11 @@ parafoil = rocket.parachutes(descentData.para, descentData.stage); % saturation on servo angle, needed if the actuation is faster than the % integration step -% if deltaA > contSettings.rocket.payload.uMax -% deltaA = contSettings.rocket.payload.uMax; +% if deltaA > contSettings.parafoil.uMax +% deltaA = contSettings.parafoil.uMax; % flagAngleSaturation = true; -% elseif deltaA < contSettings.rocket.payload.uMin -% deltaA = contSettings.rocket.payload.uMin; +% elseif deltaA < contSettings.parafoil.uMin +% deltaA = contSettings.parafoil.uMin; % flagAngleSaturation = true; % else % flagAngleSaturation = false; @@ -64,31 +64,30 @@ parafoil = rocket.parachutes(descentData.para, descentData.stage); %% CONSTANTS % environment g = environment.g0/(1 + (altitude*1e-3/6371))^2; % [N/kg] module of gravitational field -local = [environment.z0, environment.temperature, ... % vector containing inputs for atmosphereData - environment.pressure, environment.rho]; +local = environment.local; % vector containing inputs for atmosphereData % geometry -semiWingSpan = rocket.payload.semiWingSpan; % [m] wingspan -MAC = rocket.payload.MAC; % [m] mean aero chord -surface = rocket.payload.surface; % [m^2] payload surface -mass = rocket.payload.mass; % [kg] -inertia = rocket.payload.inertia; % 3x3 inertia matrix -inverseInertia = rocket.payload.inverseInertia; % 3x3 inverse inertia matrix +semiWingSpan = parafoil.semiWingSpan; % [m] wingspan +MAC = parafoil.MAC; % [m] mean aero chord +surface = parafoil.surface; % [m^2] payload surface +mass = parafoil.mass; % [kg] +inertia = parafoil.inertia; % 3x3 inertia matrix +inverseInertia = parafoil.inverseInertia; % 3x3 inverse inertia matrix % aerodynamic coefficients -CD0 = rocket.payload.CD0; CDAlpha2 = rocket.payload.CDAlpha2; -CL0 = rocket.payload.CL0; CLAlpha = rocket.payload.CLAlpha; -Cm0 = rocket.payload.Cm0; CmAlpha = rocket.payload.CmAlpha; Cmq = rocket.payload.Cmq; -Cnr = rocket.payload.Cnr; -Clp = rocket.payload.Clp; ClPhi = rocket.payload.ClPhi; +CD0 = parafoil.cd0; CDAlpha2 = parafoil.cdAlpha2; +CL0 = parafoil.cl0; CLAlpha = parafoil.clAlpha; +Cm0 = parafoil.cm0; CmAlpha = parafoil.cmAlpha; Cmq = parafoil.cmQ; +Cnr = parafoil.cnR; +Clp = parafoil.clP; ClPhi = parafoil.clPhi; % aerodynamic control coefficients - asymmetric -CnDeltaA = rocket.payload.CnDeltaA; -CDDeltaA = rocket.payload.CDDeltaA; -CLDeltaA = rocket.payload.CLDeltaA; -ClDeltaA = rocket.payload.ClDeltaA; +CnDeltaA = parafoil.cnDeltaA; +CDDeltaA = parafoil.cdDeltaA; +CLDeltaA = parafoil.clDeltaA; +ClDeltaA = parafoil.clDeltaA; % !!! % aerodynamic control coefficients - symmetric -deltaSMax = rocket.payload.deltaSMax; % max value +deltaSMax = parafoil.deltaSMax; % max value %% ROTATIONS @@ -97,7 +96,7 @@ Q = Q/norm(Q); %% ADDING WIND (supposed to be added in NED axes); if isa(wind, 'WindMatlab') - [uw, vw, ww] = wind.getVels(-altitude, t); + [uw, vw, ww] = wind.getVels(altitude, t); else [uw, vw, ww] = wind.getVels(altitude); end @@ -113,19 +112,18 @@ ur = abs(u - windVels(1)); % abs to evade problem of negative vx body in forces vr = v - windVels(2); wr = w - windVels(3); - % Body to Inertial velocities -Vels_NED = dcm'*[u; v; w]; +vNED = dcm'*[u; v; w]; % relative velocity norm -V_norm = norm([ur vr wr]); +nNorm = norm([ur vr wr]); %% ATMOSPHERE DATA absoluteAltitude = altitude + environment.z0; [~, ~, ~, rho] = atmosphereData(absoluteAltitude, g, local); %% AERODYNAMICS ANGLES -if not(abs(ur) < 1e-9 || V_norm < 1e-9) +if not(abs(ur) < 1e-9 || nNorm < 1e-9) alpha = atan(wr/ur); % beta = atan(vr/ur); % beta = asin(vr/V_norm) is the classical notation, Datcom uses this one though. % alpha_tot = atan(sqrt(wr^2 + vr^2)/ur); % datcom 97' definition @@ -138,29 +136,27 @@ end deltaANormalized = deltaA / deltaSMax; %% FORCES -qFactor = 0.5*rho*V_norm; % [Pa * s / m] dynamic pressure/vNorm +qFactor = 0.5*rho*nNorm; % [Pa * s / m] dynamic pressure/vNorm -LIFT = qFactor * (CL0 + alpha * CLAlpha + deltaANormalized * CLDeltaA) * [wr; 0; -ur]; -DRAG = qFactor * (CD0 + alpha^2 * CDAlpha2 + deltaANormalized * CDDeltaA) * [ur; vr; wr]; +lift = qFactor * (CL0 + alpha * CLAlpha + deltaANormalized * CLDeltaA) * [wr; 0; -ur]; +drag = qFactor * (CD0 + alpha^2 * CDAlpha2 + deltaANormalized * CDDeltaA) * [ur; vr; wr]; -F_AERO = (LIFT - DRAG)*surface; +forceAero = (lift - drag)*surface; % Inertial to body gravity force (in body frame): Fg = dcm*[0; 0; mass*g]; % [N] force due to the gravity in body frame % total force assembly -F = Fg + F_AERO; % [N] total forces vector +F = Fg + forceAero; % [N] total forces vector %% MOMENT -Mx = (eul(1) * ClPhi *V_norm + 0.5 * semiWingSpan* Clp * p + deltaANormalized * ClDeltaA *V_norm) * semiWingSpan; -My = (alpha * CmAlpha *V_norm+ 0.5 * MAC * Cmq * q + Cm0*V_norm) * MAC; -Mz = (r * 0.5 * semiWingSpan* Cnr + deltaANormalized * CnDeltaA*V_norm) * semiWingSpan; % bizzarre - -M_AERO = [Mx; My; Mz] * qFactor * surface; -M = M_AERO; % no inertia considerations in this model +Mx = (eul(1) * ClPhi *nNorm + 0.5 * semiWingSpan* Clp * p + deltaANormalized * ClDeltaA *nNorm) * semiWingSpan; +My = (alpha * CmAlpha *nNorm+ 0.5 * MAC * Cmq * q + Cm0*nNorm) * MAC; +Mz = (r * 0.5 * semiWingSpan* Cnr + deltaANormalized * CnDeltaA*nNorm) * semiWingSpan; % bizzarre +momentAero = [Mx; My; Mz] * qFactor * surface; %% derivatives computations -omega = [p;q;r]; +omega = [p; q; r]; % acceleration bodyAcc = F/mass - cross(omega,[u;v;w]); @@ -174,7 +170,7 @@ OM = [ 0 -p -q -r ; dQQ = 1/2*OM*Q'; %% angular acceleration -angAcc = inverseInertia*(M - cross(omega, inertia .* omega)); +angAcc = inverseInertia*(momentAero - cross(omega, inertia .* omega)); % %% actuator dynamics % @@ -195,9 +191,9 @@ angAcc = inverseInertia*(M - cross(omega, inertia .* omega)); % deltaA_ref = deltaA_ref_vec(idx_deltaA,2); % end % -% ddeltaA = (deltaA_ref-deltaA)/contSettings.rocket.payload.deltaA_tau; -% if abs(ddeltaA) >contSettings.rocket.payload.deltaA_maxSpeed -% ddeltaA = sign(deltaA_ref-deltaA)*contSettings.rocket.payload.deltaA_maxSpeed; % concettualmente sta roba è sbagliata perchè dipende dal passo di integrazione, fixare +% ddeltaA = (deltaA_ref-deltaA)/contSettings.parafoil.deltaA_tau; +% if abs(ddeltaA) >contSettings.parafoil.deltaA_maxSpeed +% ddeltaA = sign(deltaA_ref-deltaA)*contSettings.parafoil.deltaA_maxSpeed; % concettualmente sta roba è sbagliata perchè dipende dal passo di integrazione, fixare % end % % if flagAngleSaturation @@ -205,7 +201,7 @@ angAcc = inverseInertia*(M - cross(omega, inertia .* omega)); % end %% FINAL DERIVATIVE STATE ASSEMBLING -dY(1:3) = Vels_NED; +dY(1:3) = vNED; dY(4:6) = bodyAcc; dY(7:9) = angAcc; dY(10:13) = dQQ; @@ -224,7 +220,7 @@ if nargout == 2 || ~isempty(wrapper) parout.velocities = [u; v; w]; - parout.forces.aero = F_AERO; + parout.forces.aero = forceAero; parout.forces.gravity = Fg; parout.air = []; diff --git a/missions/2024_Lyra_Portugal_October/config/paraConfig.m b/missions/2024_Lyra_Portugal_October/config/paraConfig.m index d8e1d5dbacd4b82de6b2b66c81b0365fbcddd22b..24ff5c0a1a827f43fc34bab0787f04d4eda068aa 100644 --- a/missions/2024_Lyra_Portugal_October/config/paraConfig.m +++ b/missions/2024_Lyra_Portugal_October/config/paraConfig.m @@ -2,22 +2,22 @@ % parachute 1 para(1, 1) = Parachute(); -para(1,1).name = 'DROGUE chute'; -para(1,1).surface = 0.6; % [m^2] Surface -para(1,1).mass = 0.15; % [kg] Parachute Mass -para(1,1).cd = 0.75; % [/] Parachute Drag Coefficient -para(1,1).cl = 0; % [/] Parachute Lift Coefficient -para(1,1).openingDelay = 1; % [s] drogue opening delay -para(1,1).finalAltitude = 350; % [m] Final altitude of the parachute -para(1,1).cx = 1.4; % [/] Parachute Longitudinal Drag Coefficient -para(1,1).chordLength = 1.5; % [m] Shock Chord Length -para(1,1).chordK = 7200; % [N/m^2] Shock Chord Elastic Constant -para(1,1).chordC = 0; % [Ns/m] Shock Chord Dynamic Coefficient -para(1,1).m = 1; % [m^2/s] Coefficient of the surface vs. time opening model -para(1,1).nf = 12; % [/] Adimensional Opening Time -para(1,1).expulsionSpeed = 5; % [m/s] Expulsion Speed -para(1,1).forceCoefficient = 1.8; % [-] Empirical coefficient to obtain correct peak force at deployment -para(1,1).deployDuration = 0.1; % [s] Time to get the parachute to full aperture +para(1, 1).name = 'DROGUE chute'; +para(1, 1).surface = 0.6; % [m^2] Surface +para(1, 1).mass = 0.15; % [kg] Parachute Mass +para(1, 1).cd = 0.75; % [/] Parachute Drag Coefficient +para(1, 1).cl = 0; % [/] Parachute Lift Coefficient +para(1, 1).openingDelay = 1; % [s] drogue opening delay +para(1, 1).finalAltitude = 350; % [m] Final altitude of the parachute +para(1, 1).cx = 1.4; % [/] Parachute Longitudinal Drag Coefficient +para(1, 1).chordLength = 1.5; % [m] Shock Chord Length +para(1, 1).chordK = 7200; % [N/m^2] Shock Chord Elastic Constant +para(1, 1).chordC = 0; % [Ns/m] Shock Chord Dynamic Coefficient +para(1, 1).m = 1; % [m^2/s] Coefficient of the surface vs. time opening model +para(1, 1).nf = 12; % [/] Adimensional Opening Time +para(1, 1).expulsionSpeed = 5; % [m/s] Expulsion Speed +para(1, 1).forceCoefficient = 1.8; % [-] Empirical coefficient to obtain correct peak force at deployment +para(1, 1).deployDuration = 0.1; % [s] Time to get the parachute to full aperture % parachute 2 para(2, 1) = Parachute(); @@ -27,6 +27,7 @@ para(2, 1).surface = 14; % [m^2] Surface para(2, 1).mass = 1.05; % [kg] Parachute Mass para(2, 1).cd = 0.6; % [/] Parachute Drag Coefficient para(2, 1).cl = 0; % [/] Parachute Lift Coefficient +para(2, 1).openingDelay = 0; % [s] drogue opening delay para(2, 1).finalAltitude = 0; % [m] Final altitude of the parachute para(2, 1).cx = 1.15; % [/] Parachute Longitudinal Drag Coefficient para(2, 1).chordLength = 6; % [m] Shock Chord Length @@ -34,6 +35,7 @@ para(2, 1).chordK = 3000; % [N/m^2] Shock Chord Elastic Constant para(2, 1).chordC = 0; % [Ns/m] Shock Chord Dynamic Coefficient para(2, 1).m = 1; % [m^2/s] Coefficient of the surface vs. time opening model para(2, 1).nf = 8.7; % [/] Adimensional Opening Time +para(2 ,1).expulsionSpeed = 0; % [m/s] Expulsion Speed para(2, 1).forceCoefficient = 2.2; % [-] Empirical coefficient to obtain correct peak force at deployment para(2, 1).deployDuration= 0.9; % [s] Time to get the parachute to full aperture @@ -55,21 +57,35 @@ para(1, 2).chordC = 0; % [Ns/m] Shock Chord Dynamic Coefficie para(1, 2).m = 1; % [m^2/s] Coefficient of the surface vs. time opening model para(1, 2).nf = 12; % [/] Adimensional Opening Time para(1, 2).expulsionSpeed = 10; % [m/s] Expulsion Speed - +para(1, 2).forceCoefficient = 0; % [-] Empirical coefficient to obtain correct peak force at deployment +para(1, 2).deployDuration = 0; % [s] Time to get the parachute to full aperture + + % parachute 2 para(2, 2) = Parafoil(); para(2, 2).name = "Payload AIRFOIL"; -para(2, 2).mass = 0.50; % [kg] Parachute Mass +para(2, 2).mass = 0.50; % [kg] Parafoil Mass + para(2, 2).surface = 0.11; % [m^2] Surface -para(2, 2).cd = 0; % [/] Parachute Drag Coefficient -para(2, 2).cl = 0; % [/] Parachute Lift Coefficient -para(2, 2).openingDelay = 1; % [s] drogue opening delay -para(2, 2).finalAltitude = 0; % [m] Final altitude of the parachute -para(2, 2).cx = 1.4; % [/] Parachute Longitudinal Drag Coefficient -para(2, 2).chordLength = 1.5; % [m] Shock Chord Length -para(2, 2).chordK = 7200; % [N/m^2] Shock Chord Elastic Constant -para(2, 2).chordC = 0; % [Ns/m] Shock Chord Dynamic Coefficient -para(2, 2).m = 1; % [m^2/s] Coefficient of the surface vs. time opening model -para(2, 2).nf = 12; % [/] Adimensional Opening Time -para(2, 2).expulsionSpeed = 10; % [m/s] Expulsion Speed +para(2, 2).deltaSMax = 0.1; % max value + +para(2, 2).semiWingSpan = 2.55/2; % [m] settings.para(2, 2).b: semiwingspan - vela nuova: 2.55/2; - vela vecchia: 2.06/2; +para(2, 2).MAC = 0.8; % [m] mean aero chord +para(2, 2).surface = 2.04; % [m^2] parafoil surface - vela nuova 2.04; - vela vecchia: 1.64; +para(2, 2).inertia = [0.42, 0, 0.03; + 0, 0.4, 0; + 0.03, 0, 0.053]; % [kg m^2] [3x3] inertia matrix parafoil +para(2, 2).cd0 = 0.25; +para(2, 2).cdAlpha2 = 0.12; +para(2, 2).cdDeltaA = 0.01; +para(2, 2).cl0 = 0.091; +para(2, 2).clAlpha = 0.9; +para(2, 2).clDeltaA = -0.0035; +para(2, 2).clP = -0.84; +para(2, 2).clPhi = -0.1; +para(2, 2).cm0 = 0.35; +para(2, 2).cmAlpha = -0.72; +para(2, 2).cmQ = -1.49; +para(2, 2).cnR = -0.27; +para(2, 2).cnDeltaA = 0.0115; \ No newline at end of file diff --git a/missions/2024_Lyra_Portugal_October/config/rocketConfig.m b/missions/2024_Lyra_Portugal_October/config/rocketConfig.m index 5f3420171f1601851930a7eb159b39f834419b9d..ce5f8a07193d58592bc3d918c0e7462af0d3788f 100644 --- a/missions/2024_Lyra_Portugal_October/config/rocketConfig.m +++ b/missions/2024_Lyra_Portugal_October/config/rocketConfig.m @@ -29,33 +29,6 @@ payload.nosePower = 3/4; % [-] No payload.nosePMod = 1.291586e+00; % [-] P coefficient for modified nosecone shapes payload.noseCMod = 9.272342e-03; % [-] C coefficient for modified nosecone shapes -payload.paraPrev = [2, 1]; % Reference to the previous parachute in the multistage descent - % (in this case, the payload drogue: 2nd stage, 1st parachute) -payload.semiWingSpan = 2.55/2; % [m] settings.payload.b: semiwingspan - vela nuova: 2.55/2; - vela vecchia: 2.06/2; -payload.MAC = 0.8; % [m] mean aero chord -payload.surface = 2.04; % [m^2] payload surface - vela nuova 2.04; - vela vecchia: 1.64; -payload.inertiaMatrix = [0.42, 0, 0.03; - 0, 0.4, 0; - 0.03, 0, 0.053]; % [kg m^2] [3x3] inertia matrix payload+payload -%%% Aerodynamic constants -% aerodynamic control coefficients - asymmetric -payload.CD0 = 0.25; -payload.CDAlpha2 = 0.12; -payload.CL0 = 0.091; -payload.CLAlpha = 0.9; -payload.Cm0 = 0.35; -payload.CmAlpha = -0.72; -payload.Cmq = -1.49; -payload.CLDeltaA = -0.0035; -payload.Cnr = -0.27; -payload.CnDeltaA = 0.0115; -payload.CDDeltaA = 0.01; -payload.Clp = -0.84; -payload.ClPhi = -0.1; -payload.ClDeltaA = 0.01; -% aerodynamic control coefficients - symmetric -payload.deltaSMax = 0.1; % max value - %% RCS recovery = Recovery();