pk.zip_Chatter_PK法_颤振pk_颤振分析pk法_颤振计算pk法

%%%%%%%%%%%%% 赵永辉《气动弹性力学与控制》中P-K法 %%%%%%%%%%%%%%%%% clear clc % %%%%%%%%%%%%参数设置 mu=20; a=-0.4; % b=0.127; xapha=0.03; rapha=sqrt(0.25); % wh=55.9; % wapha=64.1; % Rw=wh/wapha; Rw=0.6; V=2:0.01:4; n=length(V); for j=1:n Vnon=V(j); % Vnon=1.0; % [A]{q}=z[B]{q};p^2M{q}=(V^2*B-K){q} K=[Rw^2 0;0 rapha^2]; M=[1 xapha;xapha rapha^2]; omiga=sqrt(eig(K,M)); % omiga=sort(omiga); % k10=omiga(1)/Vnon; % k20=omiga(2)/Vnon; if omiga(1)<omiga(2) k10=omiga(1)/Vnon; k20=omiga(2)/Vnon; else k10=omiga(2)/Vnon; k20=omiga(1)/Vnon; end k11=1; k=0; while abs(k11-k)>1e-6 k=k10; if k<0.5 Fk=1-0.165/(1+0.00207025/k^2)-0.335/(1+0.09/k^2); Gk=-1/k*(0.0075075/(1+0.00207025/k^2)+0.1005/(1+0.09/k^2)); else Fk=1-0.165/(1+0.001681/k^2)-0.335/(1+0.1024/k^2); Gk=-1/k*(0.006765/(1+0.001681/k^2)+0.1072/(1+0.1024/k^2)); end Lh=1-1i*2/k*(Fk+1i*Gk); Lapha=1/2-1i*1/k*(1+2*Fk+1i*Gk)-2/k^2*(Fk+1i*Gk); Mh=1/2; Mapha=3/8-1i*1/k; B=k^2/mu*[Lh Lapha-(0.5+a)*Lh;Mh-(0.5+a)*Lh Mapha-(0.5+a)*(Lapha+Mh)+(0.5+a)^2*Lh]; z=eig((Vnon)^2*B-K,M); p=-sqrt(z); if imag(p(1))<0 p(1)=-p(1); end if imag(p(2))<0 p(2)=-p(2); end % if abs(imag(p(1))-omiga(1))<abs(imag(p(2))-omiga(1)) % k11=imag(p(1))/Vnon; % else % k11=imag(p(2))/Vnon; % end % [P,I]=sort(imag(p)); % PP=[p(I(1));p(I(2))]; % % k11=imag(PP(1))/Vnon; k11=min(imag(p(1)),imag(p(2)))/Vnon; k10=k11; end if imag(p(1))<imag(p(2)) gama1(j)=real(p(1))/imag(p(1)); g1(j)=2*gama1(j); w1(j)=imag(p(1)); else gama1(j)=real(p(2))/imag(p(2)); g1(j)=2*gama1(j); w1(j)=imag(p(2)); end % % gama1(j)=real(PP(1))/imag(PP(1)); % g1(j)=2*gama1(j); % w1(j)=imag(PP(1)); %%%%%%%%%%%%%%%%% k21=1; k=0; while abs(k21-k)>1e-6 k=k20; if k<0.5 Fk=1-0.165/(1+0.00207025/k^2)-0.335/(1+0.09/k^2); Gk=-1/k*(0.0075075/(1+0.00207025/k^2)+0.1005/(1+0.09/k^2)); else Fk=1-0.165/(1+0.001681/k^2)-0.335/(1+0.1024/k^2); Gk=-1/k*(0.006765/(1+0.001681/k^2)+0.1072/(1+0.1024/k^2)); end Lh=1-1i*2/k*(Fk+1i*Gk); Lapha=1/2-1i*1/k*(1+2*Fk+1i*Gk)-2/k^2*(Fk+1i*Gk); Mh=1/2; Mapha=3/8-1i*1/k; B=k^2/mu*[Lh Lapha-(0.5+a)*Lh;Mh-(0.5+a)*Lh Mapha-(0.5+a)*(Lapha+Mh)+(0.5+a)^2*Lh]; z=eig((Vnon)^2*B-K,M); p=-sqrt(z); % if imag(p(1))*imag(p(2))<0 if imag(p(1))<0 p(1)=-p(1); end if imag(p(2))<0 p(2)=-p(2); end % if abs(imag(p(1))-omiga(2))<abs(imag(p(2))-omiga(2)) % k21=imag(p(1))/Vnon; % else % k21=imag(p(2))/Vnon; % end % [P,I]=sort(imag(p)); % PP=[p(I(1));p(I(2))]; % k21=imag(PP(2))/Vnon; k21=max(imag(p(1)),imag(p(2)))/Vnon; k20=k21; end if imag(p(1))<imag(p(2)) gama2(j)=real(p(2))/imag(p(2)); g2(j)=2*gama2(j); w2(j)=imag(p(2)); else gama2(j)=real(p(1))/imag(p(1)); g2(j)=2*gama2(j); w2(j)=imag(p(1)); end % end subplot(2,1,1) plot(V/sqrt(mu),g1,'o-') hold on plot(V/sqrt(mu),g2,'ro-') grid on subplot(2,1,2) plot(V/sqrt(mu),w1,'o-'); hold on plot(V/sqrt(mu),w2,'ro-') grid on % plot(V*wapha*b,g2,'ro-')