有大神用有限元法求解过转子的谐响应吗？

1713573225

有大神用有限元法求解过转子的谐响应吗？

1713573225

% 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





这是传递矩阵法 有限元法怎么弄呀  求大神指点！

缱绻

所谓有限元法就是用有限元的理论建立相关的运动微分方程，得到其质量矩阵，刚度矩阵以及陀螺力矩矩阵等，然后在方程的右端加上相应的激振力，形成非齐次运动微分方程
最后在此基础上用龙格库塔法或者newmark等微分方程求解方法对系统进行求解，得到系统的响应曲线