Facebook
From Rory Anderson, 1 Week ago, written in Plain Text.
Embed
Download Paste or View Raw
Hits: 62
  1. function Aeroelastic_SA_AL
  2.  
  3. % This code formulates the Structural and Aerodynamic Model for a wing
  4. % represented using an Euler-Bernoulli beam of n elements using the Finite
  5. % Element Method. It develops an aeroelastic model which is used for
  6. % flutter analysis to determine the flutter speed and flutter frequency.
  7.  
  8. % The electric propulsion system is represented as nodal lumped masses at
  9. % chosen nodes, with a thrust applied at each node.
  10.  
  11.  
  12. %%% EXPECTED OUTPUTS
  13. % Vf = flutter speed
  14. % wf = flutter frequency
  15.  
  16. %%% REQUIRED INPUTS
  17. % The code requires specific data to be inputted which relates to the model
  18. % The inputs will be split into sections - 1) Structural Model Inputs
  19. % 2) Aerodynamic Model Inputs 3) Electric Propulsion Inputs
  20.  
  21. n = 32;                         % number of elements
  22. % selected based on results of mesh convergence study
  23.  
  24. % 1) STRUCTURAL MODEL INPUTS
  25. L=7.5;                              % Length of beam [m]
  26. Lel=L/n;                            % Length of each element
  27. h = 0.24;                           % Height of beam [m]
  28. b = 2;                              % Width of beam [m]
  29.  
  30. rho_beam = 2770;                % Density [kg/m^3]
  31. E = 73.1*10^9;                  % Young's modulus [Pa]
  32. G = 28*10^9;                    % Shear Modulus
  33.  
  34. aa = b/2;
  35. bb = h/2;
  36. J = aa*(bb^3)*((16/3)-3.36*(bb/aa)*(1-(bb^4)/(12*aa^4))); % Torsion constant
  37.  
  38. %Damping Parameters
  39. alpha = 1.658525;
  40. beta = 0.000229;
  41.  
  42. % Calculations
  43. I_x = h^3*b/12;                 %Area moment of inertia - X direction
  44. I_y = b^3*h/12;                 %Area moment of inertia - Y direction
  45.  
  46. m0=b*h*rho_beam;                % Mass per unit length
  47. Xt = (m0*J)/(b*h);
  48.  
  49.  
  50. % 2) AERODYNAMIC MODEL INPUTS
  51. c = b;                  % chord
  52. dz = Lel;                 % width of strip
  53. a = 2*vpa(pi);                  % wing lift curve slope
  54. e = 0.25;                  % wing eccentricity
  55. M_thetadot = -0.3;         % Pitch damping term
  56. rho = 0.412;            % air density at 10000m
  57. %rho = 1.225;           % air density at sea level
  58.  
  59. % Range of Airspeeds
  60. V_lb = 0;
  61. V_ub = 500;
  62. V_int = 5;
  63.  
  64. % Calculations
  65. syms V
  66. q = 0.5*rho*V^2;
  67. cza = c*dz*a;
  68. ecza = e*c^2*dz*a;
  69. czM = c^3*dz*M_thetadot;
  70.  
  71.  
  72. %%% 3) ELECTRIC PROPULSION INPUTS
  73.  
  74. % Propellor properties - lumped mass
  75. massnodes = [9 17 25 33];                 % enter nodes where lumped masses exist
  76.  
  77. m = 500;            % mass of propellor [kg]
  78. ypos = 0.12;
  79. xpos = 1;
  80.  
  81. % Propellor properties - thrust
  82. Ft = 1000;
  83.  
  84.  
  85. Ix = m*ypos^2;             % moment of inertia - x
  86. Iy = m*xpos^2;             % moment of inertia - y
  87. Jz = (ypos^2+xpos^2)*m;             % polar moment of inertia
  88.  
  89.  
  90. %%% STRUCTURAL MODEL
  91. % Declare stiffness and mass matrix variables
  92. M = zeros(10);
  93. K = zeros(10);
  94.  
  95. % bernoulli beam - x direction
  96. % Thrust
  97. M_x = [(13/35)*m0*Lel, (11/210)*m0*Lel^2, (9/70)*m0*Lel, (-13/420)*m0*Lel^2;
  98.       (11/210)*m0*Lel^2, (1/105)*m0*Lel^3, (13/420)*m0*Lel^2, (-1/140)*m0*Lel^3;
  99.       (9/70)*m0*Lel, (13/420)*m0*Lel^2, (13/35)*m0*Lel, (-11/210)*m0*Lel^2;
  100.       (-13/420)*m0*Lel^2, (-1/140)*m0*Lel^3, (-11/210)*m0*Lel^2, (1/105)*m0*Lel^3];
  101.  
  102. K_x = [12*(E*I_x/Lel^3), 6*(E*I_x/Lel^2), -12*(E*I_x/Lel^3), 6*(E*I_x/Lel^2);
  103.        6*(E*I_x/Lel^2), 4*(E*I_x/Lel), -6*(E*I_x/Lel^2), 2*(E*I_x/Lel);
  104.        -12*(E*I_x/Lel^3), -6*(E*I_x/Lel^2), 12*(E*I_x/Lel^3), -6*(E*I_x/Lel^2);
  105.        6*(E*I_x/Lel^2), 2*(E*I_x/Lel), -6*(E*I_x/Lel^2), 4*(E*I_x/Lel)];  
  106.  
  107. % bernoulli beam - y direction
  108. % Lift
  109. M_y = M_x;
  110.  
  111. K_y = [12*(E*I_y/Lel^3), 6*(E*I_y/Lel^2), -12*(E*I_y/Lel^3), 6*(E*I_y/Lel^2);
  112.        6*(E*I_y/Lel^2), 4*(E*I_y/Lel), -6*(E*I_y/Lel^2), 2*(E*I_y/Lel);
  113.        -12*(E*I_y/Lel^3), -6*(E*I_y/Lel^2), 12*(E*I_y/Lel^3), -6*(E*I_y/Lel^2);
  114.        6*(E*I_y/Lel^2), 2*(E*I_y/Lel), -6*(E*I_y/Lel^2), 4*(E*I_y/Lel)];  
  115.  
  116. % torsional elements
  117. M_t = [(Xt*Lel)/3, (Xt*Lel)/6; (Xt*Lel)/6, (Xt*Lel)/3];
  118.  
  119. K_t = [(G*J)/Lel, (-G*J)/Lel; (-G*J)/Lel, (G*J)/Lel];
  120.  
  121.  
  122. % assembly of structural mass and stiffness matrices for 2-node beam element
  123. M(1:2, 1:2) = M_x(1:2,1:2);
  124. M(1:2, 6:7) = M_x(1:2, 3:4);
  125. M(3:4, 3:4) = M(1:2, 1:2);
  126. M(3:4, 8:9) = M(1:2, 6:7);
  127. M(5, 5) = M_t(1, 1);
  128. M(5, 10) = M_t(1, 2);
  129. M(6:7, 1:2) = M_x(3:4, 1:2);
  130. M(6:7, 6:7) = M_x(3:4, 3:4);
  131. M(8:9, 3:4) = M(6:7, 1:2);
  132. M(8:9, 8:9) = M(6:7, 6:7);
  133. M(10, 5) = M_t(2, 1);
  134. M(10, 10) = M_t(2, 2);
  135.  
  136. K(1:2, 1:2) = K_x(1:2,1:2);
  137. K(1:2, 6:7) = K_x(1:2, 3:4);
  138. K(3:4, 3:4) = K_y(1:2, 1:2);
  139. K(3:4, 8:9) = K_y(1:2, 3:4);
  140. K(5, 5) = K_t(1, 1);
  141. K(5, 10) = K_t(1, 2);
  142. K(6:7, 1:2) = K_x(3:4, 1:2);
  143. K(6:7, 6:7) = K_x(3:4, 3:4);
  144. K(8:9, 3:4) = K_y(3:4, 1:2);
  145. K(8:9, 8:9) = K_y(3:4, 3:4);
  146. K(10, 5) = K_t(2, 1);
  147. K(10, 10) = K_t(2, 2);
  148.  
  149. % assembly of structural mass and stiffness matrices for any number of elements
  150. if n == 1
  151.     M = M;
  152.     K = K;
  153. else
  154.     n_el = (n+1)*5;             % number of elements in matrix
  155.     M_new = zeros(n_el);        % Working variable
  156.     K_new = zeros(n_el);        % Working variable
  157.    
  158.     M_new(1:10, 1:10) = M;
  159.     K_new(1:10, 1:10) = K;
  160.    
  161.     for i=1:n-1
  162.        
  163.         M_new(i*5+1:i*5+10, i*5+1:i*5+10) = M_new(i*5+1:i*5+10, i*5+1:i*5+10) + M;
  164.         K_new(i*5+1:i*5+10, i*5+1:i*5+10) = K_new(i*5+1:i*5+10, i*5+1:i*5+10) + K;
  165.        
  166.     end
  167.    
  168.    
  169.     % Mass and Stifness matrices
  170.     M = M_new;
  171.     K = K_new;
  172.    
  173. end
  174.  
  175. % assembly of structural damping matrix using Rayleigh damping
  176. C=alpha*M+beta*K;
  177.                
  178. % QUASI-STEADY AERODYNAMIC MODEL
  179.  
  180. % assembly of aerodynamic stiffness matrix for 2-node beam element
  181. kAM = zeros(10);
  182. kAM(3,5) = cza;
  183. kAM(5,5) = ecza;
  184. kAM(6:10, 6:10) = kAM(1:5, 1:5);
  185.  
  186. %assemble matrix
  187. if n == 1;
  188.     kAM = q*kAM;
  189. else
  190.     %if number of elements > 1
  191.     kAM_new = zeros(5*(n+1));
  192.     kAM_new(1:10, 1:10) = kAM(1:10,1:10);
  193.    
  194.     for j = 1:n-1
  195.         kAM_new(i*5+1:i*5+10, i*5+1:i*5+10) = kAM_new(i*5+1:i*5+10, i*5+1:i*5+10) + kAM;
  196.    
  197.     end
  198.    
  199.     kAM = q*kAM_new;
  200.      
  201. end
  202.  
  203. % assembly of aerodynamic damping matrix for 2-node beam element
  204. cAM = zeros(10);
  205. cAM(3,3) = cza;
  206. cAM(3,5) = ecza;
  207. cAM(5,5) = czM;
  208. cAM(6:10,6:10) = cAM(1:5,1:5);
  209.  
  210. % assemble matrix
  211. if n == 1;
  212.     cAM = (q/V)*cAM;
  213. else
  214.     %if number of elements > 1
  215.     cAM_new = zeros(5*(n+1));
  216.     cAM_new(1:10, 1:10) = cAM(1:10,1:10);
  217.    
  218.     for j = 1:n-1
  219.         cAM_new(i*5+1:i*5+10, i*5+1:i*5+10) = cAM_new(i*5+1:i*5+10, i*5+1:i*5+10) + cAM;
  220.    
  221.     end
  222.    
  223.     cAM = (q/V)*cAM_new;
  224.    
  225. end
  226.  
  227.  
  228. %%% Adding propellor mass and thrust
  229. m_lumped = [ ...
  230.     m 0 0 0 0; ...
  231.     0 Ix 0 0 0; ...
  232.     0 0 m 0 0; ...
  233.     0 0 0 Iy 0; ...
  234.     0 0 0 0 Jz ];
  235.  
  236. Kt = [ ...
  237.     0 0 0 0 0; ...
  238.     0 0 0 0 0; ...
  239.     0 0 0 0 Ft; ...
  240.     0 0 0 0 0; ...
  241.     0 0 0 0 0.5*b*Ft ];
  242.  
  243. % adding propellor mass and thrust at specified nodes
  244. if isempty(massnodes) == 0  
  245.     for k = 1:length(massnodes)
  246.         mpos = massnodes(k);
  247.         M(5*mpos-4:5*mpos,5*mpos-4:5*mpos) = M(5*mpos-4:5*mpos,5*mpos-4:5*mpos) + m_lumped(1:5,1:5);
  248.         K(5*mpos-4:5*mpos,5*mpos-4:5*mpos) = K(5*mpos-4:5*mpos,5*mpos-4:5*mpos) - Kt(1:5,1:5);
  249.     end
  250. end  
  251.  
  252.  
  253. % EQUATIONS OF MOTION - assemble matrices A, B, D
  254. % EoM of the form: A*q_dotdot + B*qdot + D*q = 0
  255. % A = structural intertia (mass matrix, including lumped mass)
  256. % B = combined structural and aerodynamic damping
  257. % D = combined structural and aerodynamic stiffness, and thrust matrix
  258.  
  259. % Matrix A
  260. A = M;
  261.  
  262. % Matrix B
  263. B = C - cAM;
  264.  
  265. % Matrix D
  266. D = K - kAM;
  267.  
  268. % Apply boundary conditions
  269. A(:,1:5) = []; % = empty []
  270. B(:,1:5) = [];
  271. D(:,1:5) = [];
  272.  
  273. A(1:5,:) = [];
  274. B(1:5,:) = [];
  275. D(1:5,:) = [];
  276.  
  277. %create system matrix
  278. Q = [zeros(length(A)),eye(length(A)); -inv(A)*D, -inv(A)*B];
  279.  
  280. f_stable = 1;
  281. d_stable = 1;
  282.  
  283. % for range of airspeeds V_air
  284. for V_air = V_lb:V_int:V_ub
  285.     %V_air;
  286.    
  287.     %sub V into system matrix Q
  288.     syms V
  289.     Q_new = double(vpa(subs(Q,V,V_air)));
  290.    
  291.     % obtain eigenvalues of Q
  292.     eigval = eig(Q_new);
  293.    
  294.     % for each eigenvalue
  295.     for j = 1:length(eigval)
  296.         % if real part of complex eigenvalue is positive
  297.         if isreal(eigval(j)) == 0 && real(eigval(j)) >= 0 && f_stable == 1
  298.             % system is unstable
  299.             f_stable = 0;
  300.             % flutter speed = current airspeed
  301.             Vf = V_air;
  302.             wf=abs(imag(eigval(j)));   % flutter frequency
  303.         % if real eigenvalue is positive
  304.         elseif isreal(eigval(j)) == 1 && eigval(j) >= 0 && d_stable == 1
  305.             % system is unstable
  306.             d_stable = 0;
  307.             % divergence speed = current airspeed
  308.             Vd = V_air;
  309.         end
  310.     end
  311. end
  312.  
  313.  
  314. %%% DISPLAY RESULTS
  315. if f_stable == 1
  316.     disp("System remains dynamically stable within range of airspeeds")
  317. elseif f_stable == 0
  318.     dispVf = ['The flutter speed is ',num2str(Vf),' m/s.'];
  319.     dispwf = ['The flutter frequency is ',num2str(wf),' rad/s.'];
  320.     disp(dispVf)
  321.     disp(dispwf)
  322. end
  323.  
  324. if d_stable == 1
  325.     disp('')
  326.     disp("System remains statically stable within range of airspeeds")
  327. elseif d_stable == 0
  328.     disp('')
  329.     dispVd = ['The divergence speed is ',num2str(Vd),' m/s.'];
  330.     disp(dispVd)
  331. end
  332.  
  333. end
  334.