Skip to content
Snippets Groups Projects
Commit e796b866 authored by Marco Luigi Gaibotti's avatar Marco Luigi Gaibotti
Browse files

[Refactoring-ode][classes] Bugfix

parent fce2292c
Branches
No related tags found
1 merge request!3Refactoring ODE functions
......@@ -13,28 +13,33 @@ classdef Parafoil < Para
properties
name = ''
surface double % [m^2] Surface
deltaSMax double % [-] aerodynamic control coefficient - symmetric, max value
deltaSMax double % [-] Aerodynamic control coefficient - symmetric, max value
mass double % [kg] Parachute Mass
mass double % [kg] Parafoil 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
finalAltitude double % [m] Final altitude of the parafoil
semiWingSpan double % [m] Parafoil semiwingspan
MAC double % [m] Parafoil mean aero chord
cd0 double % [-] aerodynamic control coefficient - asymetric, CD0
cdAlpha2 double % [-] aerodynamic control coefficient - asymetric, CDAlpha2
cdDeltaA double % [-] aerodynamic control coefficient - asymetric, CDDeltaA
cdAlpha double % [-] aerodynamic control coefficient - asymetric, CDAlpha
cdSurface double % [-] aerodynamic control coefficient - asymetric, CDSurface
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
clSurface double % [-] aerodynamic control coefficient - asymetric, CLSurface
cLP double % [-] aerodynamic control coefficient - asymetric, Clp
cLPhi double % [-] aerodynamic control coefficient - asymetric, ClPhi
cLSurface
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
cNSurface double % [-] aerodynamic control coefficient - asymetric, CnDeltaA
end
properties(SetAccess = private)
......@@ -52,7 +57,13 @@ classdef Parafoil < Para
obj = obj(:, :, 1);
end
function set.inertia(obj, inertia)
obj.inertia = inertia;
obj.updateAll();
end
function updateAll(obj)
if isempty(obj.inertia), return; end
obj.inverseInertia = obj.inertia\eye(3);
end
end
......
......@@ -48,6 +48,5 @@ value = altitude - rocket.parachutes(para, stage).finalAltitude;
isterminal = 1;
direction = 0;
end
function [dY,parout] = descentParafoil(t, Y, rocket, environment, wind, descentData, wrapper)
function [dY,parout] = descentParafoil(t, Y, rocket, environment, wind, ~, descentData, wrapper)
arguments
t
Y
rocket Rocket
environment Environment
wind % WindCustom {mustBeA(wind, {'WindCustom', 'WindMatlab'})}
~
descentData struct
wrapper = [] %DataWrapper = DataWrapper.empty
end
......@@ -75,17 +76,18 @@ inertia = parafoil.inertia; % 3x3 inertia matr
inverseInertia = parafoil.inverseInertia; % 3x3 inverse inertia matrix
% aerodynamic coefficients
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;
cd0 = parafoil.cd0; cdAlpha = parafoil.cdAlpha;
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 = parafoil.cnDeltaA;
CDDeltaA = parafoil.cdDeltaA;
CLDeltaA = parafoil.clDeltaA;
ClDeltaA = parafoil.clDeltaA; % !!!
cNSurface = parafoil.cNSurface;
cdSurface = parafoil.cdSurface;
cLSurface = parafoil.cLSurface;
clSurface = parafoil.clSurface;
% aerodynamic control coefficients - symmetric
deltaSMax = parafoil.deltaSMax; % max value
......@@ -138,8 +140,8 @@ deltaANormalized = deltaA / deltaSMax;
%% FORCES
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 * cLSurface) * [wr; 0; -ur];
drag = qFactor * (cd0 + alpha^2 * cdAlpha + deltaANormalized * cdSurface) * [ur; vr; wr];
forceAero = (lift - drag)*surface;
......@@ -150,9 +152,9 @@ Fg = dcm*[0; 0; mass*g]; % [N] force due to the gravity in body frame
F = Fg + forceAero; % [N] total forces vector
%% MOMENT
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
Mx = (eul(1) * ClPhi *nNorm + 0.5 * semiWingSpan* cLP * p + deltaANormalized * clSurface *nNorm) * semiWingSpan;
My = (alpha * cMAlpha *nNorm+ 0.5 * MAC * cMQ * q + cM0*nNorm) * MAC;
Mz = (r * 0.5 * semiWingSpan* cNR + deltaANormalized * cNSurface*nNorm) * semiWingSpan; % bizzarre
momentAero = [Mx; My; Mz] * qFactor * surface;
%% derivatives computations
......@@ -170,7 +172,7 @@ OM = [ 0 -p -q -r ;
dQQ = 1/2*OM*Q';
%% angular acceleration
angAcc = inverseInertia*(momentAero - cross(omega, inertia .* omega));
angAcc = inverseInertia*(momentAero - cross(omega, inertia * omega));
% %% actuator dynamics
%
......@@ -205,9 +207,6 @@ dY(1:3) = vNED;
dY(4:6) = bodyAcc;
dY(7:9) = angAcc;
dY(10:13) = dQQ;
% dY(14) = ddeltaA;
dY = dY';
if nargout == 2 || ~isempty(wrapper)
parout.state = [];
......
......@@ -76,16 +76,24 @@ para(2, 2).surface = 2.04; % [m^2] parafoil surface - vel
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).finalAltitude = 0;
para(2, 2).cd0 = 0.25;
para(2, 2).cdAlpha2 = 0.12;
para(2, 2).cdDeltaA = 0.01;
para(2, 2).cdAlpha = 0.12;
para(2, 2).cdSurface = 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
para(2, 2).clSurface = -0.0035;
para(2, 2).cLP = -0.84;
para(2, 2).cLPhi = -0.1;
para(2, 2).cLSurface = -0.0035;
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).cNSurface = 0.0115;
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment