diff --git a/classes/components/WindMatlab.m b/classes/components/WindMatlab.m index 613fd819a01bcbbc10edb343ff51791cbe982310..b47c45d6b4cf8f56fb1c72a65b41beb5796d75cc 100644 --- a/classes/components/WindMatlab.m +++ b/classes/components/WindMatlab.m @@ -40,7 +40,7 @@ classdef WindMatlab < Component end function [uw, vw, ww] = getVels(obj, z, t, Hour, Day) - h = -z + obj.environment.z0; + h = z + obj.environment.z0; if h < 0 h = 0; end diff --git a/functions/odeFunctions/ballistic.m b/functions/odeFunctions/ballistic.m index 3baa2f2e25737032d33d04041bb9b719ffcbe667..516465a30894d75455d8be2529443fef21474c3b 100644 --- a/functions/odeFunctions/ballistic.m +++ b/functions/odeFunctions/ballistic.m @@ -62,9 +62,8 @@ q3 = Y(13); opening = Y(14); dOpening = 0; -if -z < 0 % z is directed as the gravity vector - z = 0; -end +altitude = -z; % z is directed as the gravity vector (NED reference) +% if altitude < 0, altitude = 0; end persistent isFlipped; if isempty(isFlipped) @@ -74,7 +73,7 @@ end %% CONSTANTS S = rocket.crossSection; % [m^2] cross surface C = rocket.diameter; % [m] caliber -g = environment.g0/(1 + (-z*1e-3/6371))^2; % [N/kg] module of gravitational field +g = environment.g0/(1 + (altitude*1e-3/6371))^2; % [N/kg] module of gravitational field tb = rocket.motor.cutoffTime; % [s] Burning Time local = [environment.z0, environment.temperature, ... % vector containing inputs for atmosphereData environment.pressure, environment.rho]; @@ -86,9 +85,9 @@ Q = Q/norm(Q); %% ADDING WIND (supposed to be added in NED axes); if isa(wind, 'WindMatlab') - [uw, vw, ww] = wind.getVels(z, t); + [uw, vw, ww] = wind.getVels(altitude, t); else - [uw, vw, ww] = wind.getVels(-z); + [uw, vw, ww] = wind.getVels(altitude); end dcm = quatToDcm(Q); @@ -100,19 +99,18 @@ vr = v - windVels(2); wr = w - windVels(3); % Body to Inertial velocities -Vels = dcm'*[u; v; w]; -V_norm = norm([ur vr wr]); +vels = dcm'*[u; v; w]; +velsNorm = norm([ur vr wr]); %% ATMOSPHERE DATA -absoluteAltitude = -z + environment.z0; +absoluteAltitude = altitude + environment.z0; [~, a, P, rho] = atmosphereData(absoluteAltitude, g, local); -M = V_norm/a; -M_value = M; +mach = velsNorm/a; %% TIME-DEPENDENTS VARIABLES if t < tb - m = interpLinear(rocket.motor.time, rocket.motor.mass, t); + m = interpLinear(rocket.motor.time, rocket.mass, t); T = interpLinear(rocket.motor.time, rocket.motor.thrust, t); Pe = interpLinear(rocket.motor.time, rocket.motor.pe, t); T = T + rocket.motor.ae*(Pe - P); @@ -131,7 +129,7 @@ Ixx = I(1); Iyy = I(2); Izz = I(3); Ixxdot = Idot(1); Iyydot = Idot(2); Izzdot = Idot(3); %% AERODYNAMICS ANGLES -if not(abs(ur) < 1e-9 || V_norm < 1e-9) +if not(abs(ur) < 1e-9 || velsNorm < 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 @@ -158,7 +156,7 @@ if abs(alpha)>25*pi/180 || abs(beta)>25*pi/180 coeffsValues = interpN( rocket.coefficientsHighAOA.total,... {rocket.coefficientsHighAOA.state.alphas, rocket.coefficientsHighAOA.state.machs, ... rocket.coefficientsHighAOA.state.betas, rocket.coefficientsHighAOA.state.altitudes}, ... - [alpha, M, beta, absoluteAltitude]); + [alpha, mach, beta, absoluteAltitude]); angle0 = [alpha beta]; isFlipped = true; else @@ -167,10 +165,10 @@ else coeffsValues = interpN( rocket.coefficients.total(:, :, :, :, :, 1, end),... {rocket.coefficients.state.alphas, rocket.coefficients.state.machs, ... rocket.coefficients.state.betas, rocket.coefficients.state.altitudes}, ... - [alpha, M, beta, absoluteAltitude]); + [alpha, mach, beta, absoluteAltitude]); angle0 = [alpha beta]; else - [coeffsValues, angle0] = rocket.interpCoeffs(t, alpha, M, beta, absoluteAltitude, opening); + [coeffsValues, angle0] = rocket.interpCoeffs(t, alpha, mach, beta, absoluteAltitude, opening); end end @@ -211,9 +209,9 @@ end % end %% RAMP / FREE FLIGHT -if -z < environment.effectiveRampLength*sin(omega) % No torque on the launchpad +if altitude < environment.effectiveRampLength*sin(omega) % No torque on the launchpad Fg = m*g*sin(omega); % [N] force due to the gravity - fX = 0.5*rho*V_norm^2*S*CA; + fX = 0.5*rho*velsNorm^2*S*CA; F = -Fg +T -fX; du = F/m; @@ -238,8 +236,8 @@ if -z < environment.effectiveRampLength*sin(omega) % No torque on the launc else %% FORCES % first computed in the body-frame reference system - qdyn = 0.5*rho*V_norm^2; % [Pa] dynamics pressure - qdynL_V = 0.5*rho*V_norm*S*C; + qdyn = 0.5*rho*velsNorm^2; % [Pa] dynamics pressure + qdynL_V = 0.5*rho*velsNorm*S*C; fX = qdyn*S*CA; % [N] x-body component of the aerodynamics force fY = qdyn*S*CY; % [N] y-body component of the aerodynamics force @@ -258,10 +256,10 @@ else dw = F(3)/m - p*v + q*u; % Rotation - dp = (Iyy - Izz)/Ixx*q*r + qdynL_V/Ixx*(V_norm*Cl+Clp*p*C/2) - Ixxdot*p/Ixx; - dq = (Izz - Ixx)/Iyy*p*r + qdynL_V/Iyy*(V_norm*Cm + (Cmad+Cmq)*q*C/2)... + dp = (Iyy - Izz)/Ixx*q*r + qdynL_V/Ixx*(velsNorm*Cl+Clp*p*C/2) - Ixxdot*p/Ixx; + dq = (Izz - Ixx)/Iyy*p*r + qdynL_V/Iyy*(velsNorm*Cm + (Cmad+Cmq)*q*C/2)... - Iyydot*q/Iyy; - dr = (Ixx - Iyy)/Izz*p*q + qdynL_V/Izz*(V_norm*Cn + (Cnr*r+Cnp*p)*C/2)... + dr = (Ixx - Iyy)/Izz*p*q + qdynL_V/Izz*(velsNorm*Cn + (Cnr*r+Cnp*p)*C/2)... - Izzdot*r/Izz; % Compute the aerodynamici roll angle @@ -286,7 +284,7 @@ OM = [ 0 -p -q -r ; dQQ = 1/2*OM*Q'; %% FINAL DERIVATIVE STATE ASSEMBLING -dY(1:3) = Vels; +dY(1:3) = vels; dY(4:6) = [du; dv; dw]; dY(7:9) = [dp; dq; dr]; dY(10:13) = dQQ; @@ -298,10 +296,10 @@ dY = dY'; if nargout == 2 parout.integration.t = t; - parout.interp.M = M_value; + parout.interp.M = mach; parout.interp.alpha = alpha_value; parout.interp.beta = beta_value; - parout.interp.alt = -z; + parout.interp.alt = altitude; parout.interp.mass = m; parout.interp.inertias = [Ixx, Iyy, Izz]; @@ -310,7 +308,7 @@ if nargout == 2 parout.rotations.dcm = dcm; - parout.velocities = Vels; + parout.velocities = vels; parout.forces.AeroDyn_Forces = [fX, fY, fZ]; parout.forces.T = T; diff --git a/missions/2024_Lyra_Portugal_October/config/environmentConfig.m b/missions/2024_Lyra_Portugal_October/config/environmentConfig.m index 8f66bc9fad44f1801980ef370106a7bab7a2390f..44865c27c91f14758f50da7cb947fa15192be549 100644 --- a/missions/2024_Lyra_Portugal_October/config/environmentConfig.m +++ b/missions/2024_Lyra_Portugal_October/config/environmentConfig.m @@ -4,13 +4,15 @@ environment = Environment(); -environment.lat0 = 39.388727; % [deg] Launchpad latitude -environment.lon0 = -8.287842; % [deg] Launchpad longitude -environment.z0 = 160; % [m] Launchpad Altitude +environment.lat0 = 39.388727; % [deg] Launchpad latitude +environment.lon0 = -8.287842; % [deg] Launchpad longitude +environment.z0 = 160; % [m] Launchpad Altitude +environment.omega = 85; % [deg] Launchpad elevation +environment.phi = 0; % [deg] Launchpad azimuth environment.pin1Length = 0.5603; % [m] Distance from the upper pin to the upper tank cap environment.pin2Length = 0.2055; % [m] Distance from the lower pin to the lower tank cap -environment.rampLength = 12; % [m] Total launchpad length +environment.rampLength = 12; % [m] Total launchpad length -environment.temperature = []; % [deg] Ground temperature correction -environment.pressure = []; % [Pa] Ground pressure correction -environment.rho = []; % [Kg/m^3] Gorund air density correction \ No newline at end of file +environment.temperature = []; % [deg] Ground temperature correction +environment.pressure = []; % [Pa] Ground pressure correction +environment.rho = []; % [Kg/m^3] Gorund air density correction \ 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 d64866c0fdf4fd9a3a1f12574df265237048da6b..0b91679d1f75e83d41d0fa08f0bc091c1b4bfb16 100644 --- a/missions/2024_Lyra_Portugal_October/config/rocketConfig.m +++ b/missions/2024_Lyra_Portugal_October/config/rocketConfig.m @@ -52,9 +52,6 @@ airbrakes.xCg = []; % [m] Cg relative to bay u airbrakes.multipleAB = false; % If true, multiple and smooth airbrakes opening will be simulated airbrakes.opening = [0 1 0.5]; % aerobrakes, 1-2-3 for 0%, 50% or 100% opened airbrakes.deltaTime = [5 5 5]; % aerobrakes, configurations usage time -airbrakes.multipleAB = false; % If true, multiple and smooth airbrakes opening will be simulated -airbrakes.opening = [1]; % aerobrakes, 1-2-3 for 0%, 50% or 100% opened -airbrakes.deltaTime = []; % aerobrakes, configurations usage time airbrakes.n = 3; % [-] number of brakes airbrakes.height = linspace(0, 0.0363, 3); % [m] Block airbrakes opening coordinate ( First entry must be 0! ) diff --git a/missions/2024_Lyra_Portugal_October/config/windConfig.m b/missions/2024_Lyra_Portugal_October/config/windConfig.m index 332c74eea10b207f02afd7a9d0d019b509c8ceb1..d85fceb252b9962007ee000ade640a2ac337e6f8 100644 --- a/missions/2024_Lyra_Portugal_October/config/windConfig.m +++ b/missions/2024_Lyra_Portugal_October/config/windConfig.m @@ -12,9 +12,9 @@ windCustom = WindCustom(); windCustom.altitudes = [0 200 2000]; % [m] Altitudes at which a distribution change occurs -windCustom.magnitudeDistribution = ["g", "u", "u"]; % [-] Distribution type: "u" - uniform, "g" - gaussian +windCustom.magnitudeDistribution = ["u", "u", "u"]; % [-] Distribution type: "u" - uniform, "g" - gaussian windCustom.magnitudeParameters = [7 2 10; % [m/s] Distribution parameters: "u" - [min; max], "g" - [mu; sigma] - 0.5 9 20]; + 7 2 10]; windCustom.azimuthDistribution = ["u", "u", "u"]; % [-] Distribution type: "u" - uniform, "g" - gaussian windCustom.azimuthParameters = 0*pi/180 * ones(2,3); % [deg] Distribution parameters: "u" - [min; max], "g" - [mu; sigma] diff --git a/settings/odeConfig.m b/settings/odeConfig.m index 81014b6dc2e11a7062985cf62ed9a8490a32024e..351594795377c1a8b8c8947be24c5d918e4bbd2d 100644 --- a/settings/odeConfig.m +++ b/settings/odeConfig.m @@ -4,7 +4,7 @@ ode.finalTime = 2000; % [s] Final integration time ode.optAscent = odeset('Events', @eventApogee, 'InitialStep', 1); % ODE options for ascend ode.optAscentDelayPara= odeset('InitialStep', 1); % ODE options for due to the opening delay of the parachute -ode.optAscentMultipleAB = odeset('Events', @eventAB, 'InitialStep', 1); % ODE options for opening of the airbrakes +ode.optAscentMultipleAB = odeset('Events', @eventApogee, 'InitialStep', 1); % ODE options for opening of the airbrakes ode.optParachute = odeset('Events', @eventParaCut); % ODE options for the parachutes ode.optDescent = odeset('Events', @eventLanding,... 'RelTol', 1e-3, 'AbsTol', 1e-3); % ODE options for ballistic descent