Access wn_i_all using the indexing variable i. Refer the code below.
%% Linear Aeroelastic model of a flexible wing using assumed mode shapes approach
clc;
close all;
global s c x_f m e EI GJ rho M_theta a_w A_inv C E B D
% Setting Parameters
s = 7.5; %semi span (m)
c = 2; %chord (m)
x_f = (0.48*c); %Flexural axis (m)
m = 100; %mass per unit area (kgm^-2)
e = 0.23; %Eccentricity Ratio
EI = 3675000; %Flexural rigitidy (Nm^2)
GJ = 1890000 ;%Torsional rigitidy(Nm^2)
rho = 1.225; %Air Density (Kgm^-3)
M_theta = -1.2; %Non-Dimensional pitch damping derivative
a_w = 2*pi(); %Lift Curve slove
%Inertia, Structural, Aerodynamic Damping , Aerodynamic Stiffness matrices
A_i =[(s^5) (s^4/4)*((c^2/2)-(c*x_f));
(s^4/4)*((c^2/2)-(c*x_f)) (s^3/3)*(c^3/3)-c^2*x_f+c*c*x_f^2];%Inertria metrix of system
E = [4*EI*s 0;
0 GJ*s];%Structural Matrix
B = [(c*a_w*s^5)/10 0;
(-c^2*e*a_w*s^4)/8 (-c^3*s^3*M_theta)/24];%Aerodynamic damping matrix
C = [0 (c*s^4*a_w);
0 (-e*c^2*s^3*a_w)/6];% Aerodynamic Stiffness matrix
D = [0 0;
0 0];%Structural Damping matrix
A_inv = inv(A_i);
Q = zeros(2,2);
I = [1 0; 0 1];
%Frequency Domain
V_air = linspace(0, 140, 1000);
wn_i_all = zeros(4, length(V_air));
zet_i = zeros(4, length(V_air));
for i = 1:length(V_air)
C_t = (rho*V_air(i).*(B+D));
K_t = rho*V_air(i).*(C+E);
A = [Q I ; A_inv*C_t A_inv*K_t];
[eig_v, eig_D] = eig(A);
d_eig_D_i = diag(eig_D);
d_eig_D_i = sort( d_eig_D_i );
wn_i = sqrt( real( d_eig_D_i ).^2 + imag( d_eig_D_i ).^2 );
zet_i = -real( d_eig_D_i ) ./ wn_i;
wn_i_all(:, i) = wn_i;
end
plot(wn_i_all)
%Calling ODE45 Solver
%t_range=[0 40];
%Y0=[1 2];
%[t,x]= ode45(@LAMFWM,t_range,Y0)
