riccati传递矩阵法计算转子不平衡响应.rar_Riccati_不平衡_不平衡响应_临界转速_临街转速、传递矩阵法

% transfermatrix_calculate.m % 计算传递矩阵 function EE=transfermatrix_calculate(N,wi,K,m_k,Jp_k,Jd_k,l,EI,rr,U) % ========================================================================= % 初值设为0 % ========================================================================= for i=1:1:N+1 for j=1:1:2 for k=1:1:2 ud11(j,k,i)=0; ud12(j,k,i)=0; ud21(j,k,i)=0; ud22(j,k,i)=0; end end end for i=1:1:N for j=1:1:2 Ff(j,1,i)=0; % 2*1*N维矩阵 Fe(j,1,i)=0; for k=1:1:2 us11(j,k,i)=0; us12(j,k,i)=0; us21(j,k,i)=0; us22(j,k,i)=0; end end end for i=1:1:N for j=1:1:2 for k=1:1:2 u11(j,k,i)=0; u12(j,k,i)=0; u21(j,k,i)=0; u22(j,k,i)=0; end end end % ========================================================================= % 计算质点上的传递矩阵——点矩阵的一部分！ % ========================================================================= for i=1:1:N+1 ud11(1,1,i)=1; ud11(1,2,i)=0; ud11(2,1,i)=0; ud11(2,2,i)=1; ud21(1,1,i)=0; ud21(1,2,i)=0; ud21(2,1,i)=0; ud21(2,2,i)=0; ud22(1,1,i)=1; ud22(1,2,i)=0; ud22(2,1,i)=0; ud22(2,2,i)=1; end % ========================================================================= % 计算质点上的传递矩阵——点矩阵的一部分！ % ========================================================================= for i=1:1:N+1 ud12(1,1,i)=0; ud12(2,1,i)=m_k(i)*wi^2-K(i); ud12(1,2,i)=(Jp_k(i)-Jd_k(i))*wi^2; % 考虑陀螺力矩 ud12(2,2,i)=0; end % ========================================================================= % 以下为计算无质量梁上的传递矩阵——场矩阵 % ========================================================================= for i=1:1:N us11(1,1,i)=1; us11(1,2,i)=l(i); us11(2,1,i)=0; us11(2,2,i)=1; us12(1,1,i)=0; us12(1,2,i)=0; us12(2,1,i)=0; us12(2,2,i)=0; us21(1,1,i)=l(i)^2/(2*EI(i)); us21(1,2,i)=(l(i)^3*(1-rr(i)))/(6*EI(i)); us21(2,1,i)=l(i)/EI(i); us21(2,2,i)=l(i)^2/(2*EI(i)); us22(1,1,i)=1; us22(1,2,i)=l(i); us22(2,1,i)=0; us22(2,2,i)=1; end % ========================================================================= % 计算质点上离心力产生的传递矩阵 (前N节点，质量点和轴段相结合的传递矩阵) % ========================================================================= for i=1:1:N Fe(1,1,i)=U(i)*wi^2*(l(i)^3*(1-rr(i)))/(6*EI(i)); Fe(2,1,i)=U(i)*wi^2*l(i)^2/(6*EI(i)); Ff(1,1,i)=U(i)*wi^2*l(i); Ff(2,1,i)=U(i)*wi^2; end Fe(1,1,N+1)=0; Fe(2,1,N+1)=0; % 第N+1个质量点的离心力矩阵 Ff(1,1,N+1)=0; Ff(2,1,N+1)=U(N+1)*wi^2; % ========================================================================== % 此处全为计算中间量 % ========================================================================= for i=1:1:N+2 Su(1,1,i)=0; Su(1,2,i)=0; Su(2,1,i)=0; Su(2,2,i)=0; Sn(1,1,i)=0; Sn(1,2,i)=0; Sn(2,1,i)=0; Sn(2,2,i)=0; SS(1,1,i)=0; SS(1,2,i)=0; SS(2,1,i)=0; SS(2,2,i)=0; end for i=1:1:2 P(i,1)=0; % 初始端面的边界条件 P为2*1*N维矩阵 for j=1:1:2 SS1(i,j)=0; Ud11(i,j)=0; Ud12(i,j)=0; Ud21(i,j)=0; Ud22(i,j)=0; Us11(i,j)=0; Us12(i,j)=0; Us21(i,j)=0; Us22(i,j)=0; end end % ========================================================================= % (e)形成最终传递矩阵 % ========================================================================= % Ud11,Ud12,Ud21,Ud22 为最终参与计算的传递矩阵 for i=1:1:N u11(:,:,i)=us11(:,:,i)*ud11(:,:,i)+us12(:,:,i)*ud21(:,:,i); u12(:,:,i)=us11(:,:,i)*ud12(:,:,i)+us12(:,:,i)*ud22(:,:,i); u21(:,:,i)=us21(:,:,i)*ud11(:,:,i)+us22(:,:,i)*ud21(:,:,i); u22(:,:,i)=us21(:,:,i)*ud12(:,:,i)+us22(:,:,i)*ud22(:,:,i); end u11(:,:,N+1)=ud11(:,:,N+1); u12(:,:,N+1)=ud12(:,:,N+1); u21(:,:,N+1)=ud21(:,:,N+1); u22(:,:,N+1)=ud22(:,:,N+1); for i=1:1:N+1 Ud11=u11(:,:,i); Ud12=u12(:,:,i); Ud21=u21(:,:,i); Ud22=u22(:,:,i); FF=Ff(:,:,i); FE=Fe(:,:,i); SS(:,:,i+1)=(Ud11*SS1+Ud12)*inv(Ud21*SS1+Ud22); PP(:,:,i+1)=(Ud11*P+FF)-SS(:,:,i+1)*(Ud21*P+FE); Su(:,:,i)=Ud21*SS1+Ud22; Sn(:,:,i)=inv(Ud21*SS1+Ud22); % 计算振型时用到 UF(:,:,i)=(Ud21*P+FE); SS1=SS(:,:,i+1); P=PP(:,:,i+1); end % ==========================(2)计算不平衡响应EE============================= EE(:,:,N+2)=-inv(SS(:,:,N+2))*PP(:,:,N+2); for i=N+1:-1:1 EE(:,:,i)=Sn(:,:,i)*EE(:,:,i+1)-Sn(:,:,i)*UF(:,:,i); % EE为2*1*N+1维数组 end