有大神用有限元法求解过转子的谐响应吗？
有大神用有限元法求解过转子的谐响应吗？ 本帖最后由 1713573225 于 2015-12-9 16:57 编辑{:{39}:}% 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等微分方程求解方法对系统进行求解，得到系统的响应曲线
页:
[1]