求助：newmark-beta法算齿轮扭振为啥速度和位移不收敛？

zhlei566

各位高手，我采用newmark-beta法算一个齿轮系统扭振时，出现速度和位移随着时间逐渐增大的现象，不收敛，不知是什么原因，如图所示，请高手指点一二。
程序：
clc;
clear;
close;
N=10;%自由度
dt=0.002;%步长
T=3;%总时间
%转动惯量
Jin=4.217823;
J1=0.008952085;
J2=0.190956155;
J3=0.020671561;
J4=0.766038806;
J5=0.090060;
J6=4.24048334;
J7=0.6044254;
J8=27.9071658;
Jout=58.764562;
%齿轮基圆直径m
R1=(63.171e-3)/2;
R2=(234.749e-3)/2;
R3=(97.732e-3)/2;
R4=(321.941e-3)/2;
R5=(137.755e-3)/2;
R6=(451.530e-3)/2;
R7=(213.959e-3)/2;
R8=(601.758e-3)/2;
%齿轮及轴扭转刚度Nm.rad
k1=252851;
k2=1936800;
k3=6186712;
k4=33459824;
k5=10190124;
%有限元法计算啮合刚度
k12=9.852e8;%齿轮啮合刚度N.B/m
k34=1.466e9;
k56=3.31e9;
k78=3.146e9;
Tin=819.4;%输入扭矩N.m
Tout=-91675;%输出轴扭矩N.m
% 传动轴阻尼计算
lanmda=0.01;
c1=2*lanmda*sqrt(k1/(1/Jin+1/J1));
c2=2*lanmda*sqrt(k2/(1/J2+1/J3));
c3=2*lanmda*sqrt(k3/(1/J4+1/J5));
c4=2*lanmda*sqrt(k4/(1/J6+1/J7));
c5=2*lanmda*sqrt(k5/(1/J8+1/Jout));
%齿轮副阻尼计算
lanmdag=0.05;
c12=2*lanmdag*sqrt(k12*R1^2*R2^2*J1*J2/(R1^2*J1+R2^2*J2));
c34=2*lanmdag*sqrt(k34*R3^2*R4^2*J3*J4/(R3^2*J3+R4^2*J4));
c56=2*lanmdag*sqrt(k56*R5^2*R6^2*J5*J6/(R5^2*J5+R6^2*J6));
c78=2*lanmdag*sqrt(k78*R7^2*R8^2*J7*J8/(R7^2*J7+R8^2*J8));
%扭振模型基本参数矩阵
M=[Jin 0 0 0 0 0 0 0 0 0;
    0 J1 0 0 0 0 0 0 0 0;
    0 0 J2 0 0 0 0 0 0 0;
    0 0 0 J3 0 0 0 0 0 0;
    0 0 0 0 J4 0 0 0 0 0;
    0 0 0 0 0 J5 0 0 0 0;
    0 0 0 0 0 0 J6 0 0 0;
    0 0 0 0 0 0 0 J7 0 0;
    0 0 0 0 0 0 0 0 J8 0;
    0 0 0 0 0 0 0 0 0 Jout];
K=[k1 -k1 0 0 0 0 0 0 0 0;
    -k1 k1+R1^2*k12 -R1*R2*k12 0 0 0 0 0 0 0;
    0 -R1*R2*k12 k2+R2^2*k12 -k2 0 0 0 0 0 0;
    0 0 -k2 k2+R3^2*k34 -R3*R4*k34 0 0 0 0 0;
    0 0 0 -R3*R4*k34 k3+R4^2*k34 -k3 0 0 0 0;
    0 0 0 0 -k3 k3+R5^2*k56 -R5*R6*k56 0 0 0;
    0 0 0 0 0 -R5*R6*k56 k4+R6^2*k56 -k4 0 0;
    0 0 0 0 0 0 -k4 k4+R7^2*k78 -R7*R8*k78 0;
    0 0 0 0 0 0 0 -R7*R8*k78 k5+R8^2*k78 -k5;
    0 0 0 0 0 0 0 0 -k5 k5];
C=[c1 -c1 0 0 0 0 0 0 0 0;
    -c1 c1+R1^2*c12 -R1*R2*c12 0 0 0 0 0 0 0;
    0 -R1*R2*c12 c2+R2^2*c12 -c2 0 0 0 0 0 0;
    0 0 -c2 c2+R3^2*c34 -R3*R4*c34 0 0 0 0 0;
    0 0 0 -R3*R4*c34 c3+R4^2*c34 -c3 0 0 0 0;
    0 0 0 0 -c3 c3+R5^2*c56 -R5*R6*c56 0 0 0;
    0 0 0 0 0 -R5*R6*c56 c4+R6^2*c56 -c4 0 0;
    0 0 0 0 0 0 -c4 c4+R7^2*c78 -R7*R8*c78 0;
    0 0 0 0 0 0 0 -R7*R8*c78 c5+R8^2*c78 -c5;
    0 0 0 0 0 0 0 0 -c5 c5];
tt=0:dt:T;
RecordLength=length(tt);
x=zeros(N,RecordLength);
v=zeros(N,RecordLength);
a=zeros(N,RecordLength);
%初值
x(:,1)=0;
v(:,1)=0;
a(:,1)=0;
deta=0.50;
alpha=0.25;
a0=1/alpha/dt/dt;
a1=deta/alpha/dt;
a2=1/alpha/dt;
a3=1/2/alpha-1;
a4=deta/alpha-1;
a5=dt*(deta/alpha-2)/2;
a6=dt*(1-deta);
a7=deta*dt;
for i=2:RecordLength
    t(i)=i*dt;
    P=[Tin 0 0 0 0 0 0 0 0 Tout]';%平稳运行阶段载荷激励
    K_=K+a0*M+a1*C;%考虑阻尼和不考虑阻尼
    iK=inv(K_);
    Pe=P+M*(a0*x(:,i-1)+a2*v(:,i-1)+a3*a(:,i-1))+C*(a1*x(:,i-1)+a4*v(:,i-1)+a5*a(:,i-1));%考虑阻尼
    x(:,i)=iK*Pe;
    a(:,i)=a0*(x(:,i)-x(:,i-1))-a2*v(:,i-1)-a3*a(:,i-1);
    v(:,i)=v(:,i-1)+a6*a(:,i-1)+a7*a(:,i);
end

计算结果如图片所示

dgyezw007

先检查一下刚度矩阵是否奇异，半正定系统允许刚体位移，随着时间推移位移肯定是无穷大的，但速度不应该

zhlei566

dgyezw007 发表于 2014-5-12 21:46
先检查一下刚度矩阵是否奇异，半正定系统允许刚体位移，随着时间推移位移肯定是无穷大的，但速度不应该

我检查了下，刚度矩阵为非奇异矩阵啊，应该有唯一解，纠结！位移随时间逐渐增大可以理解，可我用ode45算，速度也是随着时间趋于无穷大。

wzy123

同学我最近也遇到这个问题了，算的是螺旋锥齿轮的动力学，算出来的结果页不收敛，请问你解决了么。能加下qq767315580讨论下么

华电机械

童鞋，请问一下，你的这个问题解决了么？求指教~~~

Kevin_HIT

学习了学习了

hnczf739892267

楼主，你的问题解决了吗