%==========================================================================
%                       W H I R L   F L U T T E R 2                       %
%==========================================================================
%
close all
clear all
clc
nazev='Whirl flutter 2';
%
%==========================================================================
%                           VSTUPNÍ   PARAMETRY                           %
%==========================================================================
%
% Výbìr modelu tlumící síly
% 
viskozni_tlumeni=1;          % NASTAV HODNOTU 1 PRO VÝPOÈET S VISKÓZNÍM 
%                              TLUMENÍM (pro hodnotu 0 probìhne výpoèet bez
%                              vnitøního tlumení)
%
% Interval vlastních frekvencí izolovaných módù
%
o=1;
oo=0.01:0.01:15;              % interval vlastních frekvencí izolovaného 
%                               módu ve smìru zatáèení (horizont. rovina)
for o=1:numel(oo)
    f_psi=oo(o);              % !!! OBA INTEVRALY MUSÍ BÝT STEJNÉ !!!

n=1;
nn=0.01:0.01:15;              % interval vlastních frekvencí izolovaného 
%                               módu ve smìru klopení (vertikální rovina)
%                             
for n=1:numel(nn)
    f_theta=nn(n);
%
% Parametry proudu vzduchu 
%
H=0;                        % výška letu MSA [m]
T=15-0.0065*H;              % teplota vzduchu MSA [°C]
ro_0=1.225;                 % hustota vzduchu pøi hladinì moøe MSA [kg/m3]
ro=ro_0*(1-H/44308)^4.2553; % hustota vzduchu MSA [kg/m3]
az=(331.57+0.607*T);        % rychlost zvuku v atmosféøe v závislosti na 
%                             teplotì (aproximace pro suchý vzduch v 
%                             rozsahu teplot -100°C až 100°C) [m/s]
v_EAS=137.5;                % rychlost letu na hladinì moøe MSA [m/s]
v=v_EAS*sqrt(ro_0/ro);      % true air speed [m/s]
Mach=v/az;                  % Machovo èíslo [1]
q=0.5*ro*v^2;               % dynamický tlak [kg/m*s2]
%
% Parametry vrtule
%
Nb=4;                       % poèet listù vrtule [1]
R=1.150;                    % polomìr vrtule [m]
Dp=2*R;                     % prùmìr vrtule [m]
Fp=pi*R^2;                  % plocha vrtulového disku [m2]
r0=0.2*R;                   % radius krytu/zaèátek profilu vrtule [m]
n_max=2080;                 % max. otáèky vrtule [1/min]
Big_omega=2*pi*n_max/60;    % úhlová rychlost vrtule [rad/s]
my=v/(Big_omega*R);         % rychlostní pomìr (blade advance ratio) [1]
cr=0.1665;                  % referenèní tìtiva profilu (0,75R) [m]
Mrmax=sqrt(v^2+Big_omega*...% maximální lokální machovo èíslo [1]
      R)/az;
c_y_alpha_max=10;           % maximální stoupání nestlaèitelné vztlakové 
%                             èáry vrtulového listu
aM=c_y_alpha_max/...        % maxiální stoupání vztlakové èáry vrtulového
   sqrt(1-Mrmax^2);%          listu pøi transsonických rychlostech [1/rad]                           
%
r01=r0;                     % složky vektoru vzdálenosti øezù vrtulového 
r02=0.25*R;%                  listu mìøené od støedu vrtule poèínaje 
r03=0.30*R;%                  vzdáleností r0 za pøedním krytem až po vnìjší 
r04=0.35*R;%                  polomìr R(pro øešení aerodynamiky pomocí 
r05=0.40*R;%                  Strip Theory) [m]
r06=0.45*R;
r07=0.50*R;
r08=0.55*R;
r09=0.60*R;
r10=0.65*R;
r11=0.70*R;
r12=0.75*R;
r13=0.80*R;
r14=0.85*R;
r15=0.90*R;
r16=0.95*R;
r17=1.00*R;
%
c01=0.1357;                 % složky vektoru tìtiv jednotlivých øezù 
c02=0.1488;%                  vrtulového listu [m]
c03=0.1580;
c04=0.1679;
c05=0.1732;
c06=0.1741;
c07=0.1732;
c08=0.1725;
c09=0.1709;
c10=0.1697;
c11=0.1686;
c12=cr;
c13=0.1617;
c14=0.1571;
c15=0.1502;
c16=0.1380;
c17=0.1212;
%
a01=4.55;                   % složky vektoru stoupání vztlakové èáry 
a02=5.35;%                    profilu v jednotlivých øezech listu [1/rad]
a03=6.10;
a04=6.52;
a05=6.60;
a06=6.62;
a07=6.70;
a08=6.90;
a09=7.00;
a10=7.10;
a11=7.32;
a12=7.60;
a13=8.10;
a14=8.55;
a15=8.85;
a16=9.35;
a17=10.00;
%
%                           % pøíslušné vektory 
%
r_vektor=[r01 r02 r03 r04 r05 r06 r07 r08 r09 r10 r11 r12 r13 r14 r15...
          r16 r17];
%
a_vektor=[a01 a02 a03 a04 a05 a06 a07 a08 a09 a10 a11 a12 a13 a14 a15...
          a16 a17];
%
c_vektor=[c01 c02 c03 c04 c05 c06 c07 c08 c09 c10 c11 c12 c13 c14 c15...
          c16 c17];
%
pomerny_c_vektor=c_vektor/cr;  % vektor pomìrných (bezrozmìrných) tìtiv [1]
%
eta_vektor=r_vektor/R;      % vektor pomìrných (bezrozmìrných) polomìrù [1]
%
% Parametry motoru a jeho uložení 
%
Jx=20.61;                   % polární moment setrvaènosti samotné vrtule k 
%                             jejímu støedu (podélná osa x) [kg*m2]
Jy=126.03;                  % moment setrvaènosti celé soustavy k ose y 
%                             (vodorovná osa) [kg*m2]
Jz=126.12;                  % moment setrvaènosti celé soustavy k ose z 
%                             (svislá osa) [kg*m2]
Omega_theta=2*pi*f_theta;   % vlastní úhlová frekvence klopivého módu [Hz]
Omega_psi=2*pi*f_psi;       % vlastní úhlová frekvence zatáèivého módu [Hz]
%
K_theta=Jy*Omega_theta^2;   % Tuhost klopivé pružiny [N*m/rad]
K_psi=Jz*Omega_psi^2;       % Tuhost zatáèivé pružiny [N*m/rad]
%
gamma_theta=0.02;           % Konstrukèní tlumení ve smìru klopení [1]
gamma_psi=0.02;             % Konstrukèní tlumení ve smìru zatáèení [1]
%
ksi_theta=gamma_theta/2;    % Pomìrný útlum ve smìru klopení (viskózní) [1]
ksi_psi=gamma_psi/2;        % Pomìrný útlum ve smìru zatáèení (viskózní)[1]
%
b_theta=2*ksi_theta*...     % koeficient viskózního tlumení ve smìru
        Jy*Omega_theta;%      klopení [N*m*s/rad]
b_psi=2*ksi_psi*Jz*...      % koeficient viskózního tlumení ve smìru
      Omega_psi;%             zatáèení [N*m*s/rad]
%
a=1.045;                    % délka vyložení (od støedu vrtule k uchycení 
%                             motoru) [m]
%
%==========================================================================
%                   VÝPOÈET AERODYNAMICKÝCH KOEFICIENTÙ                   %
%==========================================================================
%
% vektor lokálních redukovaných frekvencí jednotlivých øezù listu [1]
%
kp_vektor=c_vektor./(2*R*sqrt(my^2+eta_vektor.^2));
%
% vektor složek Theodorsenovy funkce reprezentujících èást "ve fázi" 
% (opìt pro každý øez)
%
F_vektor=(besselj(1,kp_vektor).*(besselj(1,kp_vektor)+...
          bessely(0,kp_vektor))+bessely(1,kp_vektor).*...
          (bessely(1,kp_vektor)-besselj(0,kp_vektor)))./...
          (((besselj(1,kp_vektor)+bessely(0,kp_vektor)).^2)+...
          ((bessely(1,kp_vektor)-besselj(0,kp_vektor)).^2));
%
% vektor složek Theodorsenovy funkce reprezentujících èást "mimo fázi" 
% (opìt pro každý øez)
%
G_vektor=-((bessely(1,kp_vektor).*bessely(0,kp_vektor)+...
          besselj(1,kp_vektor).*besselj(0,kp_vektor))./...
          (((besselj(1,kp_vektor)+bessely(0,kp_vektor)).^2)+...
          ((bessely(1,kp_vektor)-besselj(0,kp_vektor)).^2)));
%
% integrál funkce pomìrné šíøky [1]
%
pomocny_integral=trapz(eta_vektor,pomerny_c_vektor);
%
% štíhlost vrtule [1]
%
Ar=(R*(1-r0/R)^2)/(cr*pomocny_integral);
%
% pomocná funkce pro korekci nadzvukového obtékání 
% (cut-off value of the compresible lift curve slope) [1]
%
M_fce=(Mach^2)*(1+(eta_vektor.^2)/(my^2));
%
a_fce=(1-(a_vektor/aM).^2);
%
for u=1:length(M_fce)
    
    if M_fce(1,u) >= a_fce(1,u)
        
        Mfce_final(1,u)=a_fce(1,u);
        
    else Mfce_final(1,u)=M_fce(1,u);
        
    end
end
%
% vektor hodnot funkce II1 v jednotlivých øezech 
% (k integraci po polomìru vrtule) 
%
II1_vektor=(a_vektor.*c_vektor.*F_vektor)./((((my^2)+(eta_vektor.^2)).^...
            (0.5)).*(2+Ar*((1-Mfce_final).^(0.5))));
%
% vektor hodnot funkce JJ1 v jednotlivých øezech 
% (k integraci po polomìru vrtule) 
%
JJ1_vektor=(a_vektor.*c_vektor.*G_vektor)./((((my^2)+(eta_vektor.^2)).^...
           (0.5)).*(2+Ar*((1-Mfce_final).^(0.5))));
%
% vektor hodnot funkce II2 v jednotlivých øezech 
% (k integraci po polomìru vrtule) 
%
II2_vektor=((eta_vektor.^2).*a_vektor.*c_vektor.*F_vektor)./...
           ((((my^2)+(eta_vektor.^2)).^(0.5)).*(2+Ar*((1-Mfce_final).^...
           (0.5))));
%
% vektor hodnot funkce JJ2 v jednotlivých øezech 
% (k integraci po polomìru vrtule) 
%
JJ2_vektor=((eta_vektor.^2).*a_vektor.*c_vektor.*G_vektor)./...
           ((((my^2)+(eta_vektor.^2)).^(0.5)).*(2+Ar*((1-Mfce_final).^...
           (0.5))));
%
% vektor hodnot funkce II3 v jednotlivých øezech 
% (k integraci po polomìru vrtule) 
%
II3_vektor=((eta_vektor.^4).*a_vektor.*c_vektor.*F_vektor)./...
           ((((my^2)+(eta_vektor.^2)).^(0.5)).*(2+Ar*((1-Mfce_final).^...
           (0.5))));
%
% vektor hodnot funkce JJ3 v jednotlivých øezech 
% (k integraci po polomìru vrtule) 
%
JJ3_vektor=((eta_vektor.^4).*a_vektor.*c_vektor.*G_vektor)./...
           ((((my^2)+(eta_vektor.^2)).^(0.5)).*(2+Ar*((1-Mfce_final).^...
           (0.5))));
%
% integrály hodnot funkcí II1 až JJ3 (metoda trapz) 
%
II1_integral=trapz(eta_vektor,II1_vektor);
JJ1_integral=trapz(eta_vektor,JJ1_vektor);
II2_integral=trapz(eta_vektor,II2_vektor);
JJ2_integral=trapz(eta_vektor,JJ2_vektor);
II3_integral=trapz(eta_vektor,II3_vektor);
JJ3_integral=trapz(eta_vektor,JJ3_vektor);
%
% integrály vrtulových listù I1 až J3 
%
I1=(Nb/4)*(1/(2*pi))*(my^2)*Ar*II1_integral/cr;
J1=(Nb/4)*(1/(2*pi))*(my^2)*Ar*JJ1_integral/cr;
I2=(Nb/4)*(1/(2*pi))*my*Ar*II2_integral/cr;
J2=(Nb/4)*(1/(2*pi))*my*Ar*JJ2_integral/cr;
I3=(Nb/4)*(1/(2*pi))*Ar*II3_integral/cr;
J3=(Nb/4)*(1/(2*pi))*Ar*JJ3_integral/cr;
%
% hodnoty aerodynamických koeficientù v závislosti na rychlosti letu
%
c_ztheta=-I1*4*Big_omega*cr/v;
c_zq=J2*4*Big_omega*cr/v;
c_mtheta=-J2*2*Big_omega*cr/v;
c_mq=-I3*2*Big_omega*cr/v;
c_ytheta=-J1*4*Big_omega*cr/v;
c_yq=-I2*4*Big_omega*cr/v;
c_ntheta=-I2*2*Big_omega*cr/v;
c_nq=-J3*2*Big_omega*cr/v;
%
%==========================================================================
%                 VÝPOÈET MATIC A ØEŠENÍ ROVNICE FLUTTERU                 %
%==========================================================================
%
m11=Jy;                             % prvky matice hmotnosti [kg*m2]           
m12=0;
m21=0;
m22=Jz;
%
M=[m11 m12                          % matice hmotnosti 
   m21 m22];
%
g11=0;                              % prvky gyroskopické matice [kg*m2/s]
g12=Jx*Big_omega;
g21=-g12;
g22=0;
%
G=[g11 g12                          % gyroskopická matice
   g21 g22];
%
k11=K_theta;                        % prvky matice tuhosti [N*m/rad]
k12=0;
k21=0;
k22=K_psi;
%
K=[k11 k12
   k21 k22];                        % matice tuhosti 
%
%                                   % matice tuhosti s implementací 
%                                     konstrukèního tlumení (komplexní 
%                                     forma tuhosti)
%
K_s_tlumenim=[(1+1i*gamma_theta)*k11 k12
               k21 (1+1i*gamma_psi)*k22]; 
%
d11=b_theta;                        % prvky matice viskózního tlumení 
d12=0;%                               [N*m*s/rad]
d21=0;
d22=b_psi;
%
D=[d11 d12
   d21 d22];                        % matice viskózního tlumení
%
%                                   % prvky matice aerodynamického 
%                                     tlumení [1]
%
da11=(a/Dp)*c_mtheta-0.5*c_mq-((a^2)/(Dp^2))*c_ztheta+(a/(2*Dp))*c_zq;
da12=-(a/Dp)*c_ntheta+0.5*c_nq-((a^2)/(Dp^2))*c_ytheta+0.5*(a/Dp)*c_yq;
da21=-da12;
da22=da11;
%
DA=[da11 da12                       % matice aerodynamického tlumení
    da21 da22];
%
%                                   % prvky matice aerodynamické tuhosti[1]
%
ka11=(a/Dp)*c_ztheta-c_mtheta;
ka12=c_ntheta+(a/Dp)*c_ytheta;
ka21=-ka12;
ka22=ka11;
%
KA=[ka11 ka12                       % matice aerodynamické tuhosti
    ka21 ka22];
%
% Pøevod na náhradní soustavu + výpoèet matice vlastních èísel/vlastních
% vektorù pro všechny rychlosti letu
%
M_vlnka=M;                          % celková matice hmotnosti
%                                    
%
DA_roznasobena=q*Fp*((Dp^2)/v)*DA;
%
KA_roznasobena=q*Fp*Dp*KA;
%
if viskozni_tlumeni==0
    %
    D_vlnka=G+DA_roznasobena;       % celková matice tlumení                           
    %
    K_vlnka=K+KA_roznasobena;       % celková matice tuhosti 
    %
elseif viskozni_tlumeni==1
    %
    D_vlnka=D+G+DA_roznasobena;     % celková matice tlumení
    %
    K_vlnka=K+KA_roznasobena;       % celková matice tuhosti 
    %
else
    error('Nastavení modelu tlumení na zaèátku skriptu je nesprávné.')
end
%
%                                   % náhradní soustava - matice flutteru
%
J=[zeros(2,2) eye(2,2)
   M_vlnka\-K_vlnka M_vlnka\-D_vlnka];
%
%                                   % Matice vlastních èísel pro hodnoty  
%                                     frekvencí ve vertikálním smìru
%
EIG_matrix(n,:)=eig(J);
%
n=n+1;
    % 
end  
%
EIG_matrix2=EIG_matrix;             % tøídìní vlastních èísel - vynulování 
EIG_matrix2(imag(EIG_matrix2)<=0)=0;% reálných hodnot pøíslušných statické
cols_with_all_zeros=...             % divergenci a fyzikálnì nemožným
    find(all(EIG_matrix2==0));      % módùm)                                  
%                                     
%                                   % odstranìní sloupcù po fyzikálnì
%                                   nevýznamných módech a seøazení
%                                   dopøedného a zpìtného módu)
%
EIG_matrix2(:,cols_with_all_zeros)=[];
EIG_matrix_final=sort(EIG_matrix2,2,'descend');
%
%                                   % (vymazání stabilního - dopøedného -
%                                      módu)
%
cols_with_all_real_negative=find(all(real(EIG_matrix_final)<0));
EIG_matrix_final(:,cols_with_all_real_negative)=[];
%
%                                   % Celková matice vlastních èísel
%                                     odpovídajících zpìtnému módu. Sloupce
%                                     matice odpovídají vlastním
%                                     frekvencím z intervalu psí
%                                     øádky pak z intervalu theta.                                
%                     
EIG_matrix_all_psi_theta(:,o)=EIG_matrix_final;
% 
o=o+1;
%  
end
%
%==========================================================================
%                        VYŠETØENÍ STABILITNÍ MEZE                        %
%==========================================================================
%
EIG_real=real(EIG_matrix_all_psi_theta);
%
reseni_psi=[];
reseni_theta=[];
%
for j=1:size(EIG_real,2)
    for i=2:size(EIG_real,1)
      if  EIG_real(i,j) < 0 && EIG_real(i-1,j) > 0
          reseni_theta = [reseni_theta,nn(i)];
          reseni_psi = [reseni_psi,oo(j)];
      end
        
    end
end
%
reseni_psi2=[];
reseni_theta2=[];
%
for i=1:size(EIG_real,1)
    for j=2:size(EIG_real,2)
      if  EIG_real(i,j) < 0 && EIG_real(i,j-1) > 0
          reseni_theta2 = [reseni_theta2,nn(i)];
          reseni_psi2 = [reseni_psi2,oo(j)];
      end
        
    end
end
%
%==========================================================================
%                           GRAF STABILITNÍ MEZE                          %
%==========================================================================
%
figure(1)
hold on
plot(reseni_psi,reseni_theta,'k')
plot(reseni_psi2,reseni_theta2,'k')
hold off
axis equal
xlabel('f_{\Psi} [Hz]')
ylabel('f_{\Theta} [Hz]')
xlim([0,15])
ylim([0,15])
txt = {'FLUTTER'};
text(2,4,txt);
txt = {'STABILITY'};
text(7.5,9,txt);
%
% save('PomocnySoubor.mat','reseni_psi','reseni_theta','reseni_psi2',...
%      'reseni_theta2','-append')