diff --git a/classes/Rocket.m b/classes/Rocket.m
index 7a58ba51d933ee7cedfc9e54ec3593e46c031973..8329a643282511be908adb5b0a29a4c0cea9a7cb 100644
--- a/classes/Rocket.m
+++ b/classes/Rocket.m
@@ -32,10 +32,14 @@ classdef Rocket < Component
end
properties(SetAccess = private)
+ cutoffInertia double % [kg*m^2] Total inertia at motor cutoff
+ cutoffMass double % [m] Total mass at motor cutoff
coefficients % [-] Aerodynamic coefficients
coefficientsHighAOA % [-] Aerodynamic coefficients for high angle of attacks
absolutePositions % [m] Bay start positions - 0 is set at nose base
absolutePositionsNoMotor % [m] Bay start positions, no Motor - 0 is set at nose base
+ crossSection % [m^2] Rocket cross sectional area
+
bays (1, :) Bay % [Bay] All bays
baysNoMotor (1, :) Bay % [Bay] All bays, except Motor
end
@@ -66,12 +70,17 @@ classdef Rocket < Component
obj.absolutePositionsNoMotor = [absPositions(1:4), absPositions(6:end)];
end
- function updateLength(obj)
+ function updateGeometry(obj)
if ~isempty(obj.length)
return;
end
% Measures rear - tip
obj.length = (obj.absolutePositions(end) + obj.bays(end).length) - obj.absolutePositions(1);
+
+ if ~isempty(obj.crossSection)
+ return;
+ end
+ obj.crossSection = 0.25*pi*obj.diameter^2;
end
function updateMassNoMotor(obj)
@@ -133,16 +142,22 @@ classdef Rocket < Component
motorInertia + baysInertia];
end
+ function updateCutoff(obj)
+ obj.cutoffInertia = interpLinear(obj.motor.time, obj.inertia, obj.motor.cutoffTime);
+ obj.cutoffMass = interpLinear(obj.motor.time, obj.mass, obj.motor.cutoffTime);
+ end
+
function updateAll(obj)
% Note: properties without motor must be updated first
obj.updateAbsolutePositions;
- obj.updateLength;
+ obj.updateGeometry;
obj.updateMassNoMotor;
obj.updateMass;
obj.updateXCgNoMotor;
obj.updateXCg;
obj.updateInertiaNoMotor;
obj.updateInertia;
+ obj.updateCutoff;
end
end
diff --git a/classes/bays/Motor.m b/classes/bays/Motor.m
index 41229b1f4329ab808b2798b8109cccb24be91a96..5bf438c832a61c1e9c2759115ba36cad075e53b7 100644
--- a/classes/bays/Motor.m
+++ b/classes/bays/Motor.m
@@ -90,6 +90,10 @@ classdef Motor < Bay
obj.pe = chosenMotor.Pe;
obj.ae = chosenMotor.Ae;
obj.fuselageMass = chosenMotor.mFus;
+
+ if isempty(obj.cutoffTime) || obj.cutoffTime > obj.time(end)
+ obj.cutoffTime = obj.time(end);
+ end
end
end
end
\ No newline at end of file
diff --git a/functions/odeFunctions/ballistic.m b/functions/odeFunctions/ballistic.m
new file mode 100644
index 0000000000000000000000000000000000000000..77ddb2ac87ea9de975c84785378a28f01d9570d9
--- /dev/null
+++ b/functions/odeFunctions/ballistic.m
@@ -0,0 +1,338 @@
+function [dY, parout] = ballistic(t, Y, rocket, environment, wind, newSettings, settings)
+arguments
+ t
+ Y
+ rocket Rocket
+ environment Environment
+ wind WindCustom {mustBeA(wind, {'windCustom', 'windMatlab'})}
+ newSettings Settings
+ settings
+end
+%{
+ballistic - ode function of the 6DOF Rigid Rocket Model
+
+INPUTS:
+- t, double [1, 1] integration time [s];
+- Y, double [13, 1] state vector [ x y z | u v w | p q r | q0 q1 q2 q3]:
+
+ * (x y z), NED{north, east, down} horizontal frame;
+ * (u v w), body frame velocities;
+ * (p q r), body frame angular rates;
+ * (q0 q1 q2 q3), attitude unit quaternion;
+- settings, struct(motor, CoeffsE, CoeffsF, para, ode, stoch, prob, wind), rocket data structure;
+
+OUTPUTS:
+- dY, double [16, 1] state derivatives
+- parout, struct, interesting fligth quantities structure (aerodyn coefficients, forces and so on..)
+
+
+CALLED FUNCTIONS: windMatlabGenerator, windInputGenerator, quatToDcm, interpCoeffs
+
+NOTE: To get the NED velocities the body-frame must be multiplied by the
+conjugated of the current attitude quaternion
+E.G. quatrotate(quatconj(Y(:,10:13)),Y(:,4:6))
+
+REVISIONS:
+-#0 31/12/2014, Release, Ruben Di Battista
+
+-#1 16/04/2016, Second version, Francesco Colombi
+
+-#2 01/01/2021, Third version, Adriano Filippo Inno
+
+Copyright © 2021, Skyward Experimental Rocketry, AFD department
+All rights reserved
+
+SPDX-License-Identifier: GPL-3.0-or-later
+
+%}
+
+
+% recalling the states
+% x = Y(1);
+% y = Y(2);
+z = Y(3);
+u = Y(4);
+v = Y(5);
+w = Y(6);
+p = Y(7);
+q = Y(8);
+r = Y(9);
+q0 = Y(10);
+q1 = Y(11);
+q2 = Y(12);
+q3 = Y(13);
+
+%% 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
+tb = rocket.motor.cutoffTime; % [s] Burning Time
+local = [environment.z0, environment.temperature, ... % vector containing inputs for atmosphereData
+ environment.pressure, environment.rho];
+omega = environment.omega;
+
+%% INERTIAS
+if t < tb
+ I = interpLinear(settings.motor.expTime, settings.I, t);
+ Idot = interpLinear(settings.motor.expTime, settings.Idot, t);
+else
+ I = rocket.cutoffInertia;
+ Idot = zeros(3, 1);
+end
+Ixx = I(1); Iyy = I(2); Izz = I(3);
+Ixxdot = Idot(1); Iyydot = Idot(2); Izzdot = Idot(3);
+
+%% QUATERION ATTITUDE
+Q = [q0 q1 q2 q3];
+Q = Q/norm(Q);
+
+%% ADDING WIND (supposed to be added in NED axes);
+if isa(wind, 'windMatlab')
+ [uw, vw, ww] = wind.getVels(z, t);
+else
+ [uw, vw, ww] = wind.getVels(-z);
+end
+
+dcm = quatToDcm(Q);
+wind = dcm*[uw; vw; ww];
+
+% Relative velocities (plus wind);
+ur = u - wind(1);
+vr = v - wind(2);
+wr = w - wind(3);
+
+% Body to Inertial velocities
+Vels = dcm'*[u; v; w];
+V_norm = norm([ur vr wr]);
+
+%% ATMOSPHERE DATA
+if -z < 0 % z is directed as the gravity vector
+ z = 0;
+end
+
+absoluteAltitude = -z + environment.z0;
+[~, a, P, rho] = atmosphereData(absoluteAltitude, g, local);
+
+M = V_norm/a;
+M_value = M;
+
+%% TIME-DEPENDENTS VARIABLES
+if t < tb
+ m = interpLinear(settings.motor.expTime, settings.mTotalTime, t);
+ T = interpLinear(settings.motor.expTime, settings.motor.expThrust, t);
+ Pe = interpLinear(settings.motor.expTime, settings.motor.Pe, t);
+ T = T + settings.motor.Ae*(Pe - P);
+else % for t >= tb the flight condition is the empty one (no interpolation needed)
+ m = rocket.cutoffMass;
+ T = 0;
+end
+
+%% AERODYNAMICS ANGLES
+if not(abs(ur) < 1e-9 || V_norm < 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
+else
+ alpha = 0;
+ beta = 0;
+end
+
+alpha_value = alpha;
+beta_value = beta;
+
+%% CHOSING THE EMPTY CONDITION VALUE
+% interpolation of the coefficients with the value in the nearest condition of the Coeffs matrix
+
+if t >= settings.tControl && M <= settings.MachControl
+ c = settings.control;
+else
+ c = 1;
+end
+
+%% INTERPOLATE AERODYNAMIC COEFFICIENTS:
+
+[coeffsValues, angle0] = interpCoeffs(t, alpha, M, beta, absoluteAltitude,...
+ c, settings);
+
+% Retrieve Coefficients
+
+
+CA = coeffsValues(1); CYB = coeffsValues(2); CY0 = coeffsValues(3);
+CNA = coeffsValues(4); CN0 = coeffsValues(5); Cl = coeffsValues(6);
+Clp = coeffsValues(7); Cma = coeffsValues(8); Cm0 = coeffsValues(9);
+Cmad = coeffsValues(10); Cmq = coeffsValues(11); Cnb = coeffsValues(12);
+Cn0 = coeffsValues(13); Cnr = coeffsValues(14); Cnp = coeffsValues(15);
+%--------------------------------------------
+%Clb = coeffsValues(16);
+%--------------------------------------------
+
+% XCP_value = coeffsValues(16);
+
+% compute CN,CY,Cm,Cn (linearized with respect to alpha and beta):
+alpha0 = angle0(1); beta0 = angle0(2);
+
+CN = (CN0 + CNA*(alpha - alpha0));
+CY = (CY0 + CYB*(beta - beta0));
+Cm = (Cm0 + Cma*(alpha - alpha0));
+Cn = (Cn0 + Cnb*(beta - beta0));
+
+% XCPlon = Cm/CN;
+
+if abs(alpha) <= pi/180
+ XCPlon = Cma/CNA;
+else
+ XCPlon = Cm/CN;
+end
+
+% XCPlat = Cn/CY;
+
+
+if abs(beta) <= pi/180
+ XCPlat = Cnb/CYB;
+else
+ XCPlat = Cn/CY;
+end
+
+% if Cn == 0 && CY == 0
+% XCPlat = -5;
+% end
+
+
+
+%%
+if -z < settings.lrampa*sin(omega) % No torque on the launchpad
+
+ Fg = m*g*sin(omega); % [N] force due to the gravity
+ X = 0.5*rho*V_norm^2*S*CA;
+ F = -Fg +T -X;
+ du = F/m;
+
+ dv = 0;
+ dw = 0;
+ dp = 0;
+ dq = 0;
+ dr = 0;
+
+ alpha_value = NaN;
+ beta_value = NaN;
+ Y = 0;
+ Z = 0;
+ XCPlon = NaN;
+ XCPlat = NaN;
+
+ if T < Fg % No velocity untill T = Fg
+ du = 0;
+ end
+
+ XCPtot = NaN;
+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;
+
+ X = qdyn*S*CA; % [N] x-body component of the aerodynamics force
+ Y = qdyn*S*CY; % [N] y-body component of the aerodynamics force
+ Z = qdyn*S*CN; % [N] z-body component of the aerodynamics force
+ Fg = dcm*[0; 0; m*g]; % [N] force due to the gravity in body frame
+ F = Fg + [-X+T, Y, -Z]'; % [N] total forces vector
+
+ %-----------------------------------------------------
+ %F = Fg + [-X+T*cos(chi), Y+T*sin(chi), -Z]'; % [N] total forces vector
+ %-----------------------------------------------------
+
+%% STATE DERIVATIVES
+ % velocity
+ du = F(1)/m - q*w + r*v;
+ dv = F(2)/m - r*u + p*w;
+ 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)...
+ - Iyydot*q/Iyy;
+ dr = (Ixx - Iyy)/Izz*p*q + qdynL_V/Izz*(V_norm*Cn + (Cnr*r+Cnp*p)*C/2)...
+ - Izzdot*r/Izz;
+
+ % Compute the aerodynamici roll angle
+ [~, phi] = getAlphaPhi(alpha, beta);
+
+ % Aerodynamic-force coefficient in the alpha-total plane
+ CFaTot = sin(phi)*CY + cos(phi)*(-CN);
+ % Aerodynanic-moment coefficient in the alpha-total plane
+ CMaTot = cos(phi)*Cm - sin(phi)*Cn;
+ if CFaTot ~= 0
+ XCPtot = CMaTot/CFaTot;
+ else
+ XCPtot = XCPlon;
+ end
+end
+% Quaternions
+OM = [ 0 -p -q -r ;
+ p 0 r -q ;
+ q -r 0 p ;
+ r q -p 0 ];
+
+dQQ = 1/2*OM*Q';
+
+%% FINAL DERIVATIVE STATE ASSEMBLING
+dY(1:3) = Vels;
+dY(4) = du;
+dY(5) = dv;
+dY(6) = dw;
+dY(7) = dp;
+dY(8) = dq;
+dY(9) = dr;
+dY(10:13) = dQQ;
+dY = dY';
+
+%% SAVING THE QUANTITIES FOR THE PLOTS
+
+if nargout == 2
+ parout.integration.t = t;
+
+ parout.interp.M = M_value;
+ parout.interp.alpha = alpha_value;
+ parout.interp.beta = beta_value;
+ parout.interp.alt = -z;
+ parout.interp.mass = m;
+ parout.interp.inertias = [Ixx, Iyy, Izz];
+
+ parout.wind.NED_wind = [uw, vw, ww];
+ parout.wind.body_wind = wind;
+
+ parout.rotations.dcm = dcm;
+
+ parout.velocities = Vels;
+
+ parout.forces.AeroDyn_Forces = [X, Y, Z];
+ parout.forces.T = T;
+
+ parout.air.rho = rho;
+ parout.air.P = P;
+
+ parout.accelerations.body_acc = [du, dv, dw];
+ parout.accelerations.ang_acc = [dp, dq, dr];
+
+ parout.coeff.CA = CA;
+ parout.coeff.CYB = CYB;
+ parout.coeff.CNA = CNA;
+ parout.coeff.Cl = Cl;
+ parout.coeff.Clp = Clp;
+ %--------------------------------------------
+ %parout.coeff.Clb = Clb;
+ %--------------------------------------------
+ parout.coeff.Cma = Cma;
+ parout.coeff.Cmad = Cmad;
+ parout.coeff.Cmq = Cmq;
+ parout.coeff.Cnb = Cnb;
+ parout.coeff.Cnr = Cnr;
+ parout.coeff.Cnp = Cnp;
+ parout.coeff.XCPlon = XCPlon;
+ parout.coeff.XCPlat = XCPlat;
+
+
+ parout.coeff.XCPtot = XCPtot;
+
+
+end
\ No newline at end of file
diff --git a/missions/2024_Lyra_Roccaraso_September/config/rocketConfig.m b/missions/2024_Lyra_Roccaraso_September/config/rocketConfig.m
index cb93a0f901521b5bd7cd4f22d15f2bda20051bc4..6851a934aa0667627d38b07e681024668e787a2d 100644
--- a/missions/2024_Lyra_Roccaraso_September/config/rocketConfig.m
+++ b/missions/2024_Lyra_Roccaraso_September/config/rocketConfig.m
@@ -15,7 +15,7 @@ rocket.xCgNoMotor = []; % [m] OVERRIDE
%% PLD - Includes Payload + Nose
payload = Payload();
-payload.length = 0.6; % [m] Total bay length
+payload.length = 0.6; % [m] Total bay length
payload.mass = 3.38512; % [kg] Total bay mass
payload.inertia = ...
[1032892397; 5461775539; 5450863094]*1e-9; % [kg*m^2] Total bay inertia (Body reference)
@@ -73,7 +73,7 @@ airbrakes.servoOmega = 150*pi/180; % [rad/s] Servo-motor
motor = Motor();
motor.name = 'HRE_FURIA-Rv2-T04T03'; % [-] Motor name
-motor.cutoffTime = inf; % [s] Cutoff time
+motor.cutoffTime = []; % [s] OVERRIDE Cutoff time
motor.ignitionTransient = 0.4; % [s] Ignition transient
motor.cutoffTransient = 0.3; % [s] Cut-off transient