Skip to content

Instantly share code, notes, and snippets.

@mbebar

mbebar/Stability_Control.m

Created November 15, 2017 07:49
Show Gist options
  • Star 3 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save mbebar/844e4186f452d403f4ab714fbe5c43d6 to your computer and use it in GitHub Desktop.
Save mbebar/844e4186f452d403f4ab714fbe5c43d6 to your computer and use it in GitHub Desktop.
Download ZIP
MATLAB - Aircraft Stability & Control
Raw
Stability_Control.m
%**************************************************************************
% AENG 4400 - Stability & Control
% Final Project - Fall 2016
%**************************************************************************
%
%
%
% "Will Engineer for coffee"
%
% ...not kidding. Please send coffee ASAP (Espresso Preferred)
%
%**************************************************************************
% Boeing 747-200, Cessna T-37, McDonnell Douglas F-4 Phantom
%**************************************************************************
clc,clear
%--------------------------(Select an Aircraft)---------------------------%
aircraft = selectAircraft;
if strcmp(aircraft,'Boeing 747-200')
[m,g,b,S,qbar,MAC,VP1,IxxB,IyyB,IzzB,IxzB,CD1,CDu,CDalpha,CTx1, ...
CTxu,CL0,CL1,CLu,CLalpha,CL_alphaDot,CLq,CM0,CM1,CMu,CMalpha, ...
CM_alphaDot,CMq,CMT1,CMTu,CMTalpha,CD_deltaE,CL_deltaE,CM_deltaE, ...
Clbeta,Clp,Clr,CYbeta,CYp,CYr,Cnbeta,CnTbeta,Cnp,Cnr,Cl_deltaA, ...
Cl_deltaR,Cy_deltaA,Cy_deltaR,Cn_deltaA,Cn_deltaR] = B747;
elseif strcmp(aircraft,'Cessna T-37')
[m,g,b,S,SH,qbar,MAC,VP1,IxxB,IyyB,IzzB,IxzB,CD1,CDu,CDalpha,CTx1, ...
CTxu,CL0,CL1,CLu,CLalpha,CL_alphaDot,CLq,CM0,CM1,CMu,CMalpha, ...
CM_alphaDot,CMq,CMT1,CMTu,CMTalpha,CD_deltaE,CL_deltaE,CM_deltaE, ...
Clbeta,Clp,Clr,CYbeta,CYp,CYr,Cnbeta,CnTbeta,Cnp,Cnr,Cl_deltaA, ...
Cl_deltaR,Cy_deltaA,Cy_deltaR,Cn_deltaA,Cn_deltaR] = T37;
end
%******************************[Problem 1]*******************************%
%---------------(Non-Dimensional Aerodynamic Coefficients)---------------%
%Longitudinal Aerodynamic Coefficients-----------------------------------%
% Lift Coefficient Sensitivity due to Angle of Attack
%CL_alpha = CL_alphaW + CL_alphaH*(etaH)*(S_H/S)*(1-downwash);
% Moment Coefficient Sensitivity due to Angle of Attack
%CM_alpha = (CL_alphaW*(xbarCG - xbarAC)) ...
%- (CL_alphaH)*(etaH)*(S_H/S)*(1-downwash)*(xbarAC_H - xbarCG);
% Moment Coefficient Sensitivity due to Angle of Attack Rate of Change
%CM_alphaDot = (-2*CL_alphaH*etaH)*(S_H/S) ...
%* ((xbarAC_H - xbarCG)^2)*(downwash);
% Moment Coefficient Sensitivity due to Pitch Rate
%CMq = (-2*CL_alphaH*etaH)*(S_H/S)*((xbarAC_H - xbarCG)^2);
%Lateral Aerodynamic Coefficients----------------------------------------%
%Cl_beta = Cl_beta; % -> Contribution to Rolling Moment due to Side Force
%Clp = Clp; % -> Contribution to Rolling Moment due to Roll Rate
%Cn_beta = Cn_beta; % -> Contribution to Yawing Moment due to Side Force
%Cnr = Cnr; % -> Contribution to Yawing Moment due to Yaw Rate
%******************************[Problem 2]*******************************%
%---------------------(Trim Angle of Attack @ Cruise)--------------------%
% Determine Trim Angle of Attack
AOA1_num = [(CL1-CL0) CL_deltaE; (-CM0) CM_deltaE];
AOA1_den = [CLalpha CL_deltaE; CMalpha CM_deltaE];
AOA1 = (det(AOA1_num))/(det(AOA1_den)); % Steady State AOA (rad)
% Transfer Product and Moments of Inertia from Body Axis to Stability Axis
A = [(cos(AOA1))^2, (sin(AOA1))^2, - sin(2*AOA1) ;
(sin(AOA1) )^2, (cos(AOA1))^2, sin(2*AOA1);
(0.5*sin(2*AOA1)), ( - 0.5*sin(2*AOA1)), cos(2*AOA1)];
ib = [IxxB, IzzB, IxzB]';
c = A*ib;
% Product and Moments of Inertia Transferred to Stability Axis
Ixx = c(1,1);
Izz = c(2,1);
Ixz = c(3,1);
I1 = (Ixz/Ixx);
I2 = (Ixz/Izz);
%******************************[Problem 3]*******************************%
%--------(Longitudinal & Lateral Dimensional Derivatives @ Cruise)-------%
%************************************************%
% Longitudinal Dimensional Derivatives
%------------------------------------------------%
Xu = ((-qbar*S)*(CDu + (2*CD1)))/(m*VP1); % (sec^-1)
Xalpha = ((-qbar*S)*(CDalpha - CL1))/(m); % (ft/sec^2)
Zu = ((-qbar*S)*(CLu + 2*CL1))/(m*VP1); % (sec^-1)
Z_alphaDot = -(qbar*S*MAC*CL_alphaDot)/(2*m*VP1); % (ft/sec)
Mu = ((qbar*S*MAC)*(CMu + 2*CM1))/(VP1*IyyB); % (1/ft*sec)
Malpha = (qbar*S*MAC*CMalpha)/(IyyB); % (1/sec^2)
M_alphaDot = ((qbar*S*MAC*CM_alphaDot)/(IyyB))*(MAC/(2*VP1)); % (sec^-1)
XTu = ((qbar*S)*(CTxu + 2*CTx1))/(m*VP1); % (sec^-1)
Zalpha = ((-qbar*S)*(CLalpha + CD1))/(m); % (ft/sec^2)
Zq = -(qbar*S*MAC*CLq)/(2*m*VP1); % (ft/sec)
MTu = ((qbar*S*MAC)*(CMTu + 2*CMT1))/(VP1*IyyB); % (1/ft*sec)
MTalpha = (qbar*S*MAC*CMTalpha)/(IyyB); % (sec^-2)
Mq = ((qbar*S*MAC*CMq)/(IyyB))*(MAC/(2*VP1)); % (sec^-1)
X_deltaE = -(qbar*S*CD_deltaE)/(m); % (ft/sec^2)
Z_deltaE = -(qbar*S*CL_deltaE)/(m); % (ft/sec^2)
M_deltaE = (qbar*S*MAC*CM_deltaE)/(IyyB); % (sec^-2)
%------------------------------------------------------------------------%
% Export to Function File
% Longitudinal Dimensional Derivatives -> "longitudinalTFCoefficients.m"
[A1,B1,C1,D1,E1,Au,Bu,Cu,Du,A_alpha,B_alpha,C_alpha,D_alpha, ...
A_theta,B_theta,C_theta] = longitudinalTFCoefficients(g,AOA1,VP1, ...
Xu,Xalpha,Zu,Z_alphaDot,Mu,Malpha,M_alphaDot,XTu,Zalpha,Zq, ...
MTu,MTalpha,Mq,X_deltaE,Z_deltaE,M_deltaE);
%************************************************%
% Lateral Dimensional Derivatives
%------------------------------------------------%
Ybeta = (qbar*S*CYbeta)/(m); % (ft/sec^2)
Yp = ((qbar*S*CYp)/(m))*(b/(2*VP1)); % (ft/sec^2)
Yr = ((qbar*S*CYr)/(m))*(b/(2*VP1)); % (ft/sec^1)
Lbeta = (qbar*S*b*Clbeta)/(Ixx); % (sec^-2)
Lp = (qbar*S*(b^2)*Clp)/(2*Ixx*VP1); % (sec^-1)
Lr = (qbar*S*(b^2)*Clr)/(2*Ixx*VP1); % (sec^-1)
Nbeta = (qbar*S*b*Cnbeta)/(Izz); % (1/sec^2)
NTbeta = (qbar*S*b*CnTbeta)/(Izz); % (sec^-2)
Np = (qbar*S*(b^2)*Cnp)/(2*Izz*VP1); % (sec^-1)
Nr = (qbar*S*(b^2)*Cnr)/(2*Izz*VP1); % (sec^-1)
Y_deltaA = (qbar*S*Cy_deltaA)/(m); % (ft/sec^2)
Y_deltaR = (qbar*S*Cy_deltaR)/(m); % (ft/sec^2)
L_deltaA = (qbar*S*b*Cl_deltaA)/(Ixx); % (sec^-2)
L_deltaR = (qbar*S*b*Cl_deltaR)/(Ixx); % (sec^-2)
N_deltaA = (qbar*S*b*Cn_deltaA)/(Izz); % (sec^-2)
N_deltaR = (qbar*S*b*Cn_deltaR)/(Izz); % (sec^-2)
%------------------------------------------------------------------------%
% Export to Function File
% Lateral Dimensional Derivatives -> "lateralTFCoefficients.m"
[A2,B2,C2,D2,E2,A_betaA,B_betaA,C_betaA,D_betaA,A_phiA,B_phiA, ...
C_phiA,A_psiA,B_psiA,C_psiA,D_psiA, A_betaR,B_betaR,C_betaR, ...
D_betaR,A_phiR,B_phiR,C_phiR,A_psiR,B_psiR,C_psiR,D_psiR] = ...
lateralTFCoefficients(g,VP1,I1,I2,Ybeta,Yp,Yr,Lbeta,Lp,Lr, ...
Nbeta,Np,Nr,Y_deltaA,Y_deltaR,L_deltaA,L_deltaR,N_deltaA,N_deltaR);
%******************************[Problem 4]*******************************%
%----------(Longitudinal & Lateral Stability of Flight Dynamics)---------%
%------------------------(Longitudinal Stability)------------------------%
% Longitudinal Characteristic Equation
Den1 = [A1 B1 C1 D1 E1];
% Velocity with respect to Elevator Deflection
%Num_u = Au*(s^3) + Bu*(s^2) + Cu*(s) + Du;
Num_u = [Au Bu Cu Du];
% Angle of Attack with respect to Elevator Deflection
Num_AOA = [A_alpha B_alpha C_alpha D_alpha];
% Pitch Angle with respect to Elevator Deflection
Num_Pitch = [A_theta B_theta C_theta];
%******************************************************%
%------------------------------------------------------%
% Longitudinal Transfer Functions: Elevator Deflection
u_E = tf(Num_u,Den1);
AOA_E = tf(Num_AOA,Den1);
Pitch_E = tf(Num_Pitch,Den1);
%------------------------------------------------------%
%******************************************************%
% Determine the Longitudinal Dynamic Stability of the Aircraft
Routh = D1*((B1*C1)-(A1*D1))-((B1^2)*B1);
if A1 <= 0
disp('Longitudinal dynamic stability NOT satisfied!')
elseif B1 <= 0
disp('Longitudinal dynamic stability NOT satisfied!')
elseif C1 <= 0
disp('Longitudinal dynamic stability NOT satisfied!')
elseif D1 <= 0
disp('Longitudinal dynamic stability NOT satisfied!')
elseif B1 <= 0
disp('Longitudinal dynamic stability NOT satisfied!')
elseif Routh <= 0
disp('Longitudinal dynamic stability NOT satisfied!')
end
%---------------------------(Lateral Stability)--------------------------%
% Lateral Characteristic Equation
Den2 = [A2 B2 C2 D2 E2 0];
% Angle of Sideslip with respect to Aileron Deflection
Num_betaA = [A_betaA B_betaA C_betaA D_betaA 0];
% Roll Angle with respect to Aileron Deflection
Num_phiA = [A_phiA B_phiA C_phiA 0];
% Yaw Angle with respect to Aileron Deflection
Num_psiA = [A_psiA B_psiA C_psiA D_psiA];
% Angle of Sideslip with respect to Rudder Deflection
Num_betaR = [A_betaR B_betaR C_betaR D_betaR 0];
% Roll Angle with respect to Rudder Deflection
Num_phiR = [A_phiR B_phiR C_phiR 0];
% Yaw Angle with respect to Rudder Deflection
Num_psiR = [A_psiR B_psiR C_psiR D_psiR];
%***********************************************%
%-----------------------------------------------%
% Lateral Transfer Functions: Aileron Deflection
Sideslip_A = tf(Num_betaA,Den2);
Roll_A = tf(Num_phiA,Den2);
Yaw_A = tf(Num_psiA,Den2);
%-----------------------------------------------%
% Lateral Transfer Functions: Rudder Deflection
Sideslip_R = tf(Num_betaR,Den2);
Roll_R = tf(Num_phiR,Den2);
Yaw_R = tf(Num_psiR,Den2);
%-----------------------------------------------%
%***********************************************%
% Determine the Lateral Dynamic Stability of the Aircraft
Routh_Lateral = D2*((B2*C2)-(A2*D2))-((B2^2)*B2);
if A2 <= 0
disp('Lateral directional dynamic stability NOT satisfied!')
elseif B2 <= 0
disp('Lateral-directional dynamic stability NOT satisfied!')
elseif C2 <= 0
disp('Lateral-directional dynamic stability NOT satisfied!')
elseif D2 <= 0
disp('Lateral-directional dynamic stability NOT satisfied!')
elseif B2 <= 0
disp('Lateral-directional dynamic stability NOT satisfied!')
elseif Routh_Lateral <= 0
disp('Lateral-directional dynamic stability NOT satisfied!')
end
%******************************[Problem 5]*******************************%
%------(Short Period, Phugoid, Dutch Roll, and Roll Characteristics)-----%
%************************************************%
% Solve and Plot Poles for Longitudinal Dynamics
%------------------------------------------------%
poles_1 = roots(Den1);
subplot(121)
plot(poles_1,'*'); grid on
title('Longitudinal Poles (S-Domain)')
xlabel('Real Axis')
ylabel('Imaginary Axis')
% Natural Frequency for Longitudinal Short Period Mode
omega_sp = sqrt(abs(poles_1(1,1)^2));
% Natural Frequency for Longitudinal Phugoid Mode
omega_ph = sqrt(abs(poles_1(3,1)^2));
% Damping Ratio for Longitudinal Short Period Mode
zeta_sp = abs(real(poles_1(1,1)))/(omega_sp);
% Damping Ratio for Longitudinal Phugoid Mode
zeta_ph = abs(real(poles_1(3,1)))/(omega_ph);
%******************************************************%
% Solve and Plot Poles for Lateral-Directional Dynamics
%------------------------------------------------------%
poles_2 = roots(Den2);
subplot(122)
plot(poles_2,'*'); grid on
title('Lateral-Directional Poles (S-Domain)')
xlabel('Real Axis')
ylabel('Imaginary Axis')
% Natural Frequency for Second-Order Lateral Dynamics (Dutch Roll Mode)
omega_dr = sqrt(abs(poles_2(2,1)^2));
% Damping Ratio for Second-Order Lateral Dynamics (Dutch Roll)
zeta_dr = abs(real(poles_2(2,1)))/(omega_dr);
% Rise Time for First-Order Lateral Dynamics (Rolling Mode)
Tr_rolling = (-1)/(poles_2(4,1));
% Rise Time for First-Order Lateral Dynamics (Spiral Mode)
Tr_spiral = (-1)/(poles_2(5,1));
%******************************[Problem 6]*******************************%
%-----------------(Simulation of Aircraft Flight Dynamics)---------------%
% Time column vector for simulation <12001 time steps)
t = [0:0.025:300]';
%******************************%
% Elevator Maneuvers
%------------------------------%
% Basic Elevator Maneuvers
Library = menu('Elevator Inputs', 'Single Doublet Impulse',...
'Multiple Doublet Impulses', 'Single Doublet', ...
'Multiple Doublets');
% Single Doublet Impulse (Elevator)
if Library == 1
for i = 1:200
de(i,1) = 0.0/57.3;
end;
for i = 201:205
de(i,1) = -4/57.3;
end;
for i = 206:1200
de(i,1) = 0.0/57.3;
end;
for i = 1201: 1205
de(i,1) = 4/57.3;
end;
for i = 1206:12001
de(i,1) = 0.0/57.3;
end;
end;
% Multiple Doublets Impulse (Elevator)
if Library == 2
echo off;
for i = 1:200
de(i,1) = 0.0/57.3;
end;
for i = 201:215
de(i,1) = -4/57.3;
end;
for i = 216:400
de(i,1) = 0.0/57.3;
end;
for i = 401:415
de(i,1) = 4.0/57.3;
end;
for i = 416:800
de(1,1) = 0.0/57.3;
end;
for i = 801:815
de( 1,1) = -2.0/57.3;
end;
for i = 816:1000
de(1,1) = 0.0/57.3;
end;
for i = 1001:1015
de(i,1) = 2.0/57.3;
end;
for i = 1016:12001
de(i,1) = 0.0/57.3;
end;
end;
% Single Doublet (Elevator)
if Library == 3
for i = 1:200
de(i,1) = 0.0/57.3;
end;
for i = 201:400
de(i,1) = -4/57.3;
end;
for i = 401:1200
de(i,1) = 0.0/57.3;
end;
for i =1201:1400
de(i,1) = 4/57.3;
end;
for i =1401:12001
de(i,1) = 0.0/57.3;
end;
end;
% Multiple Doublets (Elevator)
if Library == 4
for i = 1:200
de(i,1) = 0.0/57.3;
end;
for i = 201: 300
de(i,1) = -2/57.3;
end;
for i = 301:1600
de(i,1) = 0.0/57.3;
end;
for i = 1601: 1700
de(i,1) = 2/57.3;
end;
for i = 1701:3000
de(i,1) = 0.0/57.3;
end;
for i = 3001: 3100
de(i,1) = -1/57.3;
end;
for i = 3101:4400
de(i,1) = 0.0/57.3;
end;
for i = 4401: 4500
de(i,1) = 1/57.3;
end;
for i = 4501:12001
de(i,1) = 0.0/57.3;
end;
end;
% Simulation for Elevator Maneuvers
Esys_1 = AOA_E;
alphaE = lsim(Esys_1,de,t);
Esys_2 = u_E;
uE = lsim(Esys_2,de,t);
Esys_3 = Pitch_E;
pitchE = lsim(Esys_3,de,t);
% Plotting for Elevator Maneuvers
ans_E = menu('Responses to Elevator Maneuver', 'Angle of Attack', 'Velocity', 'Pitch', 'Elevator', 'All');
if ans_E == 1
close all;
alphaE_deg = alphaE*(180/pi);
plot(t,alphaE_deg, 'r');
title ('Angle of Attack vs. Time for Elevator Maneuver');
xlabel ('Time (s)');
ylabel ('Alpha (deg)');
grid on;
end;
if ans_E == 2
close all;
plot(t,uE, 'k');
title('Small Peturbation Velocity vs. Time for Elevator Maneuver ');
xlabel('Time (s)');
ylabel('u (ft/sec)');
grid on;
end;
if ans_E == 3
close all;
pitchE_deg = pitchE*(180/pi);
plot(t,pitchE_deg, 'b');
title('Small Peturbation Pitch Angle vs. Time for Elevator Maneuver');
xlabel ('Time (s)');
ylabel ('Theta (deg)');
grid on;
end;
if ans_E == 4
close all;
de_deg = de*(180/pi);
plot (t,de_deg, 'm');
title ('Elevator Deflection vs. Time');
xlabel ('Time (s)');
ylabel ('Elevator Deflection (deg) ');
grid on;
end;
if ans_E == 5
close all;
alphaE_deg = alphaE * (180/pi);
pitchE_deg = pitchE * (180/pi);
de_deg = de * (180/pi);
subplot(221),plot(t,alphaE_deg,'r');
title ('Angle of Attack vs. Time for Elevator Maneuver');
xlabel ('Time (s)');
ylabel ('Alpha (deg)');
grid on;
subplot(222),plot(t,uE,'k');
title ('Velocity vs. Time for Elevator Maneuver');
xlabel ('Time (s)');
ylabel ('u (ft/s)');
grid on;
subplot(223),plot(t,pitchE_deg,'b');
title ('Pitch Angle vs. Time for Elevator Maneuver');
xlabel ('Time (s)');
ylabel ('Theta (deg) ');
grid on;
subplot(224),plot(t,de_deg,'m');
title ('Elevator Deflection vs. Time');
xlabel ('Time (s)');
ylabel ('Elevator Deflection (deg) ');
grid on;
end;
%******************************%
% Aileron and Rudder Maneuvers
%------------------------------%
% Basic Ailerons/Rudder Maneuvers
Library = menu ('AILERON/RUDDER INPUTS', 'Single Doublet Impulse Aileron',...
'Multiple Doublet Impulses Aileron', 'Single Doublet Aileron', ...
'Multiple Doublets Aileron', 'Single Doublet Impulse Rudder', ...
'Multiple Doublet Impulses Rudder', 'Single Doublet Rudder',...
'Multiple Doublets Rudder');
% Single Impulse Aileron
if Library == 1
for i = 1:200
da(i,1) =0.0/57.3;
end;
for i = 201:205
da(i,1) =-8/57.3;
end;
for i = 206:3200
da(i,1) =0.0/57.3;
end;
for i = 3201: 3205
da(i,1) =5/57.3;
end;
for i = 3206:12001
da(i,1) =0.0/57.3;
end;
end;
%Multiple Impulse Aileron
if Library == 2
echo off;
for i = 1:200
da(i,1) =0.0/57.3;
end;
for i = 201:205
da(i,1) =-5/57.3;
end;
for i = 206:900
da(i,1) =0.0/57.3;
end;
for i = 901: 905
da(i,1) =5/57.3;
end;
for i = 906: 1800
da(1,1)=0.0/57.3;
end;
for i = 1801:1805
da( 1,1) =-21/57.3;
end;
for i = 1806: 3000
da(1,1)=0.0/57.3;
end;
for i = 3001:3005
da(i,1)=2/57.3;
end;
for i = 3006: 12001
da(i,1)=0.0/57.3;
end;
end;
%Single Step Aileron
if Library == 3
for i = 1:380
da(i,1)=0.0/57.3;
end;
for i = 381: 400
da(i,1)=-4/57.3;
end;
for i = 401: 3280
da(i,1)=0.0/57.3;
end;
for i =328:3300
da(i,1)=4/57.3;
end;
for i =3301: 12001
da(i,1)=0.0/57.3;
end;
end;
% Multiple Steps Aileron
if Library == 4
for i = 1:380
da(i,1)=0.0/57.3;
end;
for i = 381: 400
da(i,1)=-4/57.3;
end;
for i = 401:1680
da(i,1)=0.0/57.3;
end;
for i = 1681: 1700
da(i,1)=4/57.3;
end;
for i = 1701:3080
da(i,1)=0.0/57.3;
end;
for i = 3081: 3100
da(i,1)=-2/57.3;
end;
for i = 3101:4080
da(i,1)=0.0/57.3;
end;
for i = 4081: 4100
da(i,1)=2/57.3;
end;
for i = 4101:12001
da(i,1)=0.0/57.3;
end;
end;
% Single Impulse Doublet Rudder
if Library == 5
for i = 1:200
dr(i,1)= 0.0/57.3;
end;
for i = 201:205
dr(i,1) =-4/57.3;
end;
for i = 206: 1200
dr(i,1)=0.0/57.3;
end;
for i = 1201:1205
dr(i,1)=4/57.3;
end;
for i = 1206:12001
dr(i,1)=0.0/57.3;
end;
end;
%Multiple Impulse Doublet Rudder
if Library == 6
for i = 1:200
dr(i,1)=0.0/57.3;
end;
for i = 201:215
dr(i,1) =-4/57.3;
end;
for i = 216:400
dr(i,1)=0.0/57.3;
end;
for i = 401:415
dr(i,1)=4/57.3;
end;
for i = 416:800
dr(i,1)=0.0/57.3;
end;
for i = 801: 815
dr(i,1) =-4/57.3;
end;
for i = 816:1000
dr(i,1)=0.0/57.3;
end;
for i = 1001:1015
dr(i,1)=4/57.3;
end;
for i = 1016:12001
dr(i,1)=0.0/57.3;
end;
end;
% Single Doublet Rudder
if Library == 7
for i = 1:360
dr(i,1)=0.0/57.3;
end;
for i = 361:400
dr(i,1) =-5/57.3;
end;
for i = 401:1360
dr(i,1)=0.0/57.3;
end;
for i = 1361:1400
dr(i,1)=5/57.3;
end;
for i = 1401:12001
dr(i,1)=0.0/57.3;
end;
end;
% Multiple Doublets Rudder
if Library == 8
for i = 1:300
dr(i,1) =0.0/57.3;
end;
for i = 361:400
dr(i,1) =-5/57.3;
end;
for i = 401:760
dr(i,1) =0.0/57.3;
end;
for i = 761:800
dr(i,1) =5/57.3;
end;
for i = 801:1160
dr(i,1) =0.0/57.3;
end;
for i = 1161:1200
dr(i,1) =-2.5/57.3;
end;
for i = 1201:1860
dr(i,1) =0.0/57.3;
end;
for i = 1861:1900
dr(i,1) =2.5/57.3;
end;
for i = 1901:12001
dr(i,1) =0.0/57.3;
end;
end;
% Simulation
if Library<= 4
sys_1 = Sideslip_A;
betaA = lsim(sys_1,da,t);
sys_2 = Roll_A;
phiA = lsim(sys_2,da,t);
sys_3 = Yaw_A;
psiA = lsim(sys_3,da,t);
end;
if Library > 4
sys_4 = Sideslip_R;
betaR = lsim(sys_4,dr,t);
sys_5 = Roll_R;
phiR = lsim(sys_5,dr,t);
sys_6 = Yaw_R;
psiR = lsim(sys_6,dr,t);
end;
% Plotting
if Library <= 4
ans_a = menu('Responses to Aileron Maneuver', 'Beta', 'Phi', 'Psi', 'Aileron', 'All');
if ans_a == 1
close all;
betaA_deg = betaA*(180/pi);
plot(t,betaA_deg, 'r');
title ('Beta vs. Time for Aileron Maneuver');
xlabel ('Time (s)');
ylabel ('Beta (deg)');
grid on;
end;
if ans_a == 2
close all;
phiA_deg = phiA*(180/pi);
plot(t,phiA_deg, 'k');
title('Phi vs. Time for Aileron Maneuver');
xlabel('Time (s)');
ylabel('Phi (deg)');
grid on;
end;
if ans_a == 3
close all;
psiAdeg = psiA*(180/pi);
plot(t,psiAdeg, 'b');
title('Psi vs. Time for Aileron Maneuver');
xlabel ('Time (s)');
ylabel ('Psi (deg)');
grid on;
end;
if ans_a == 4
close all;
da_deg = da*(180/pi);
plot (t,da_deg, 'm');
title ('Aileron Deflection vs. Time');
xlabel ('Time (s)');
ylabel ('Aileron Deflection (deg) ');
grid on;
end;
if ans_a == 5
close all;
betaA_deg = betaA * (180/pi);
phiA_deg = phiA * (180/pi);
psiA_deg = psiA * (180/pi);
da_deg = da * (180/pi);
subplot(221),plot(t,betaA_deg,'r');
title ('Beta vs. Time for Aileron Maneuver');
xlabel ('Time (s)');
ylabel ('Beta (deg)');
grid on;
subplot(222),plot(t,phiA_deg,'k');
title ('Phi vs. Time for Aileron Maneuver');
xlabel ('Time (s)');
ylabel ('Phi (deg)');
grid on;
subplot(223),plot(t,psiA_deg,'b');
title ('Psi vs. Time for Aileron Maneuver');
xlabel ('Time (s)');
ylabel ('Psi (deg) ');
grid on;
subplot(224),plot(t,da_deg,'m');
title ('Aileron Deflection vs. Time');
xlabel ('Time (s)');
ylabel ('Aileron Deflection (deg) ');
grid on;
end;
end;
if Library >4
ans_r = menu('Responses to Rudder Maneuver', 'Beta', 'Phi', 'Psi', 'Rudder', 'All');
if ans_r == 1
close all;
betaR_deg = betaR*(180/pi);
plot(t,betaR_deg, 'r');
title ('Beta vs. Time for Rudder Manuever');
xlabel ('Time (s)');
ylabel ('Beta (deg)');
grid on;
end;
if ans_r == 2
close all;
phiR_deg = phiR*(180/pi);
plot(t,phiR_deg, 'k');
title ('Phi vs. Time for Rudder Maneuver');
xlabel ('Time (s) ');
ylabel ('Phi (deg)');
grid on;
end;
if ans_r == 3
close all;
psiR_deg = psiR*(180/pi);
plot(t,psiR_deg, 'b');
title('Psi vs. Time for Rudder Maneuver');
xlabel('Time (s)');
ylabel('Psi (deg)');
grid on;
end;
if ans_r == 4
close all;
dr_deg = dr*(180/pi);
plot (t,dr_deg,'m');
title('Rudder Deflection vs. Time');
xlabel ('Time (s)');
ylabel ('Rudder Deflection (deg)');
grid on;
end;
if ans_r == 5
close all;
betaR_deg = betaR * (180/pi);
phiR_deg = phiR * (180/pi);
psiR_deg = psiR * (180/pi);
dr_deg = dr * (180/pi);
subplot(221),plot(t,betaR_deg,'r');
title('Beta vs. Time for Rudder Maneuver');
xlabel('Time (s)');
ylabel('Bet a (deg)');
grid on;
subplot(222),plot(t,phiR_deg,'k');
title('Phi vs. Time for Rudder Maneuver');
xlabel('Time (s)');
ylabel('Phi (deg) ');
grid on;
subplot(223),plot(t,psiR_deg,'b');
title ('Psi vs. Time for Rudder Maneuver');
xlabel ('Time (s)');
ylabel('Psi (deg)');
grid on;
subplot(224),plot(t,dr_deg,'m');
title('Rudder Deflection vs. Time');
xlabel('Time (s)');
ylabel('Rudder Deflection (deg)');
grid on;
end;
end;
%******************************[Problem 7]*******************************%
%-------------------------(Sensitivity Analysis)-------------------------%
% short period natural frequency and damping ratio equations
S_H_S = SH/S;
%CM_alphadot = -2*CL_alphaH*NH*(S_H_S)*(xbar_ACH-xbarCG)^2*de_dalpha;
%CMq = -2.2*CL_alphaH*NH*(S_H_S)*(xbar_ACH-xbarCG)^2;
%CLalpha = CL_alphaW+CL_alphaW*NH*SH_S*(1-de_dalpha);
%CMalpha = CL_alphaW*(xbarCG-xbarACWB)-CL_alphaH*NH*SH_S*(1-de_dalpha)*(xbarACH-xbarCG);
%M_alphadot = (CM_alphadot*qbar*S*MAC)/(2*Iyy*VP1);
%Mq = (CMq*qbar*S*MAC^2);
%M_alpha = (CMalpha*qbar*S*MAC)/Iyy;
%Z_alpha = (-qbar*S*(CLalpha+CD1))/m;
%wn_SP = sqrt(((Zalpha*Mq)/VP1)-Malpha);
%Z_SP = -(M_alphadot+Mq+(Zalpha/VP1))/(2*wn_SP);
% Dutch Roll Natural Frequency and Damping Ratio and the rolling and spiral
% time constants
%Lr = (qbar*S*Clr/Ixx)*(b/(2*U));
%Lp = (qbar*S*Clp/Ixx)*(b/(2*U));
%Ybeta = (qbar*S*Cybeta/m);
%Yp = (qbar*S*Cyp/m)*(b/(2*U));
%Yr = (qbar*S*Cyr/m)*(b/(2*U));
%Lbeta = (qbar*S*Clbeta*b/Ixx);
%Nbeta = (qbar*S*Cnbeta*b/Izz);
%Np = (qbar*S*b*Cnp/Izz)*(b/(2*U));
%Nr = (qbar*S*b*Cnr/Izz)*(b/(2*U));
%I1 = Ixz/Ixx;
%I2 = Ixz/Izz;
%A2 = VP1*(1-(I1*I2));
%B2 = -Ybeta*(1-(I1*I2))- VP1*(Lp+Nr+(I1*Np)+(I2*Lr));
%C2 = VP1*((Lp*Nr)-(Lr*Np))+ Ybeta*(Nr+Lp+(I1*Np)+(I2*Lr)) ...
% - Yp*(Lbeta+(Nbeta*I1)) + VP1*((Lbeta*I2)+Nbeta) ...
% - Yr*((Lbeta*I2)+Nbeta);
%D2 = -Ybeta*((Lp*Nr)-(Np*Lr)) + Yp*((Lbeta*Nr)-(Lr*Nbeta))...
% - g*(Lbeta+(I1*Nbeta)) + VP1*((Lbeta*Np)-(Lp*Nbeta))...
% - Yr*((Lbeta*Np)-(Lp*Nbeta));
%E2 = g*((Lbeta*Nr)-(Lr*Nbeta));
% Lateral Characteristic Equation
%Den2 = [A2 B2 C2 D2 E2 0];
% Determine the Lateral Dynamic Stability of the Aircraft
%Routh_Lateral = D2*((B2*C2)-(A2*D2))-((B2^2)*B2);
%if A2 <= 0
%disp('Lateral directional dynamic stability NOT satisfied!')
%elseif B2 <= 0
%disp('Lateral-directional dynamic stability NOT satisfied!')
%elseif C2 <= 0
%disp('Lateral-directional dynamic stability NOT satisfied!')
%elseif D2 <= 0
%disp('Lateral-directional dynamic stability NOT satisfied!')
%elseif B2 <= 0
%disp('Lateral-directional dynamic stability NOT satisfied!')
%elseif Routh_Lateral <= 0
%disp('Lateral-directional dynamic stability NOT satisfied!')
%end
%poles_2 = roots(Den2);
% Natural Frequency for Second-Order Lateral Dynamics (Dutch Roll Mode)
%omega_dr = sqrt(abs(poles_2(2,1)^2));
% Damping Ratio for Second-Order Lateral Dynamics (Dutch Roll)
%zeta_dr = abs(real(poles_2(2,1)))/(omega_dr);
% Rise Time for First-Order Lateral Dynamics (Rolling Mode)
%Tr_rolling = (-1)/(poles_2(4,1));
% Rise Time for First-Order Lateral Dynamics (Spiral Mode)
%Tr_spiral = (-1)/(poles_2(5,1));
xbarCG = 0.25*(MAC);
figure
subplot(221),plot(t,betaA_deg,'r');
title ('Beta vs. Time for Aileron Maneuver');
xlabel ('Time (s)');
ylabel ('Beta (deg)');
grid on;
subplot(221),plot(zeta_sp,xbarCG,'r');
title('Damping SP vs. CG');
xlabel('Center of Gravity');
ylabel('Damping Ratio');
grid on;
subplot(222),plot(omega_sp,xbarCG)
title('Natural Frequency SP vs. CG');
xlabel('Center of Gravity');
ylabel('Natural Frequency');
grid on;
subplot(223),plot(zeta_sp,S_H_S)
title('Damping SP vs. SH/S');
xlabel('Center of Gravity');
ylabel('Damping Ratio');
grid on;
subplot(224),plot(omega_sp,S_H_S)
title('Damping SP vs. CG');
xlabel('Center of Gravity');
ylabel('Damping Ratio');
grid on;
%plot(zeta_dr,Clbeta)
%plot(Tr_rolling,Clbeta)
%plot(Tr_spiral,Clbeta)
%plot(zeta_dr,Cnbeta)
%plot(Tr_rolling,Cnbeta)
%plot(Tr_spiral,Cnbeta)
%plot(zeta_dr,Cnr)
%plot(Tr_rolling,Cnr)
%plot(Tr_spiral,Cnr)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment