建立转子系统有限元模型 newmark算法得到结果超大 怎么回事

0502dongfang

两根轴  联轴器联接 划分了20个单元
clc;clear;close;

E=200e6;
I=3.97e-8;
L=0.061;
k1=BeamElement(E,I,L);
I=1.0235e-007;
L=0.09;
k2=BeamElement(E,I,L);
I=1.0235e-007;
L=0.112;
k3=BeamElement(E,I,L);
I=1.2566e-007;
L=0.116;
k4=BeamElement(E,I,L);
I=1.2566e-007;
L=0.095;
k5=BeamElement(E,I,L);
I=1.2566e-007;
L=0.061;
k6=BeamElement(E,I,L);
I=1.2566e-007;
L=0.2;
k7=BeamElement(E,I,L);
I=1.2566e-007;
L=0.185;
k8=BeamElement(E,I,L);
I=1.2566e-007;
L=0.05;
k9=BeamElement(E,I,L);
I=1.0235e-007;
L=0.106;
k10=BeamElement(E,I,L);
I=1.0235e-007;
L=0.02;
k11=BeamElement(E,I,L);
I=1.0235e-007;
L=0.151;
k12=BeamElement(E,I,L);
I=1.0235e-007;
L=0.13;
k13=BeamElement(E,I,L);
I=1.2566e-007;
L=0.104;
k14=BeamElement(E,I,L);
I=1.2566e-007;
L=0.183;
k15=BeamElement(E,I,L);
I=1.2566e-007;
L=0.2;
k16=BeamElement(E,I,L);
I=1.2566e-007;
L=0.2;
k17=BeamElement(E,I,L);
I=1.2566e-007;
L=0.171;
k18=BeamElement(E,I,L);
I=1.0235e-007;
L=0.139;
k19=BeamElement(E,I,L);
I=1.0235e-007;
L=0.06;
k20=BeamElement(E,I,L);
K=zeros(42,42);
K=BeamAssemble(K,k1,1,2);
K=BeamAssemble(K,k2,2,3);
K=BeamAssemble(K,k3,3,4);
K=BeamAssemble(K,k4,4,5);
K=BeamAssemble(K,k5,5,6);
K=BeamAssemble(K,k6,6,7);
K=BeamAssemble(K,k7,7,8);
K=BeamAssemble(K,k8,8,9);
K=BeamAssemble(K,k9,9,10);
K=BeamAssemble(K,k10,10,11);
K=BeamAssemble(K,k11,11,12);
K=BeamAssemble(K,k12,12,13);
K=BeamAssemble(K,k13,13,14);
K=BeamAssemble(K,k14,14,15);
K=BeamAssemble(K,k15,15,16);
K=BeamAssemble(K,k16,16,17);
K=BeamAssemble(K,k17,17,18);
K=BeamAssemble(K,k18,18,19);
K=BeamAssemble(K,k19,19,20);
K=BeamAssemble(K,k20,20,21);%整体刚度矩阵


L=0.061;
U=0.3393/L;
m1=MBeamElement(U,L);
L=0.09;
U=0.8008/L;
m2=MBeamElement(U,L);
L=0.112;
U=0.9966/L;
m3=MBeamElement(U,L);
L=0.116;
U=1.1437/L;
m4=MBeamElement(U,L);
L=0.095;
U=0.9367/L;
m5=MBeamElement(U,L);
L=0.061;
U=0.6014/L;
m6=MBeamElement(U,L);
L=0.2;
U=1.9719/L;
m7=MBeamElement(U,L);
L=0.185;
U=1.824/L;
m8=MBeamElement(U,L);
L=0.05;
U=0.493/L;
m9=MBeamElement(U,L);
L=0.106;
U=0.9432/L;
m10=MBeamElement(U,L);
L=0.02;
U=0.178/L;
m11=MBeamElement(U,L);
L=0.151;
U=1.3436/L;
m12=MBeamElement(U,L);
L=0.13;
U=1.1568/L;
m13=MBeamElement(U,L);
L=0.104;
U=1.0254/L;
m14=MBeamElement(U,L);
L=0.183;
U=1.8043/L;
m15=MBeamElement(U,L);
L=0.2;
U=1.9719/L;
m16=MBeamElement(U,L);
L=0.2;
U=1.9719/L;
m17=MBeamElement(U,L);
L=0.171;
U=1.686/L;
m18=MBeamElement(U,L);
L=0.139;
U=1.2369/L;
m19=MBeamElement(U,L);
L=0.06;
U=0.5339/L;
m20=MBeamElement(U,L);
M=zeros(42,42);
M=MBeamAssemble(M,m1,1,2);
M=MBeamAssemble(M,m2,2,3);
M=MBeamAssemble(M,m3,3,4);
M=MBeamAssemble(M,m4,4,5);
M=MBeamAssemble(M,m5,5,6);
M=MBeamAssemble(M,m6,6,7);
M=MBeamAssemble(M,m7,7,8);
M=MBeamAssemble(M,m8,8,9);
M=MBeamAssemble(M,m9,9,10);
M=MBeamAssemble(M,m10,10,11);
M=MBeamAssemble(M,m11,11,12);
M=MBeamAssemble(M,m12,12,13);
M=MBeamAssemble(M,m13,13,14);
M=MBeamAssemble(M,m14,14,15);
M=MBeamAssemble(M,m15,15,16);
M=MBeamAssemble(M,m16,16,17);
M=MBeamAssemble(M,m17,17,18);
M=MBeamAssemble(M,m18,18,19);
M=MBeamAssemble(M,m19,19,20);
M=MBeamAssemble(M,m20,20,21);   %整体质量矩阵

a=2*[0.02/1520+0.02/6080]/[1/6080^2-1/1520^2];%×è?á????
b=(0.02*6080-0.02*1520)/(6080^2-1520^2);
C=a*M+b*K;  %阻尼矩阵


%Newmark  算法
n=length(M);
t(1)=0;
z(:,1)=zeros(n,1);
z1(:,1)=zeros(n,1);
z2(:,1)=zeros(n,1);
gama=0.5;
dt=1;
delta=0.25;
a0=1/(delta*dt^2);
a1=gama/(delta*dt);
a2=1/(delta*dt);
a3=1/(2*delta)-1;
a4=gama/delta-1;
a5=dt*(gama/(2*delta)-1);
a6=dt*(1-gama);
a7=gama*dt;
K1=K+a0*M+a1*C; %等效刚度矩阵
t_max=10;

j=1;
t(1)=0;
q=zeros(42,1);
while t(j)<t_max;
Mc=5.6;        
r=0.0002;        
w=1800;
e1=r*sin(w*t(j));
e2=r*cos(w*t(j));
frx=-Mc*e1*(2*w^2)*cos(w*t(j));     %齿套激振力
fry=-Mc*e1*(2*w^2)*sin(w*t(j));
fRx=-Mc*e2*(2*w^2)*cos(w*t(j));
fRy=Mc*e1*(2*w^2)*cos(w*t(j));
Frx=-(6.56e3)*cos(w*t(j));       %扭矩产生力
Fry=-(6.56e3)*sin(w*t(j));
FRx=-(6.56e3)*sin(w*t(j));
FRy=(6.56e3)*cos(w*t(j));
Fr=([0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 frx+Frx fry+Fry fRx+FRx fRy+FRy 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0])';
%重力
G=([0 -19 0 0 0 0 0 0 0 -51.6 0 -263.2 0 -263.2 0 0 0 0 0 0 0 -13.05 0 -13.05 0 0 0 0 0 -218.7 0 -218.7 0 -218.7 0 -218.7 0 0 0 0 0 0])';
Q1=0.04*0.003*w^2;  
Q2=0.04*0.007*w^2;
FQ1x=-0.04*0.003*w^2*cos(w*t(j)+5*pi/18);%不平衡力
FQ1y=-0.04*0.003*w^2*sin(w*t(j)+5*pi/18);
FQ2x=-0.04*0.007*w^2*cos(w*t(j)+pi/3);
FQ2y=0.04*0.007*w^2*sin(w*t(j)+pi/3);
Q=([0 0 0 0 0 0 0 0 FQ1x FQ1y 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 FQ2x FQ2y 0 0 0 0 0 0 0 0 0 0])';
F=Fr+Q+G; %整体受力  没有考虑轴承油膜力
   q(:,j+1)=F+M*(a0*z(:,j)+a2*z1(:,j)+a3*z2(:,j))+C*(a1*z(:,j)+a4*z(:,j)+a5*z2(:,j));                        
   z(:,j+1)=inv(K1)*q(:,j+1);
   z2(:,j+1)=a0*(z(:,j+1)-z(:,j))-a2*z1(:,j)-a3*z2(:,j);
   z1(:,j+1)=z1(:,j)+a6*z2(:,j)+a7*z2(:,j+1);
   %z1(:,j+1)=a1*(z(:,j+1)-z(:,j))-a4*z1(:,j)-a5*z2(:,j);
      j=j+1;
      t(j)=t(j-1)+dt;
end  
  

1713573225

考虑油膜力应该如何做

xiaowuu

楼主你好，我也是算出来结果超大，你怎么解决的

desolate

有可能是刚度矩阵中刚度值比较大，所以造成矩阵运算过程的截断误差

可以考虑对系统模型进行归一化处理

dollfish000

适当加些阻尼