%========================================================================== % 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')