齿轮变刚度计算出现错误

redstonegu

clc;
N=3;
ha=1;ns=5000;
z1=40;z2=33;z3=106;m=3;a=22*pi/180; pb=pi*m*cos(a);
ra1=m*z1/2+ha*m;ra2=m*z2/2+ha*m;ra3=m*z3/2-ha*m;
rb1=m*z1*cos(a)/2;rb2=m*z2*cos(a)/2;rb3=m*z3*cos(a)/2;
%质量
ms=4.75;m1=2.88;mr=16.92;
%j1=m1*rb1^2;j2=m2*rb2^2;j3=m3*rb3^2;
M=[ms,0,0;
   0,m1,0;
   0,0,mr];
%%%刚度
l1=sqrt(ra1^2-rb1^2)+sqrt(ra2^2-rb2^2)-(m*z1/2+m*z2/2)*sin(a);
l2=sqrt(ra2^2-rb2^2)-(sqrt(ra3^2-rb3^2)-(m*z3/2-m*z2/2)*sin(a));
e1=l1/pb;e2=l2/pb;
w=2*pi*ns*z1/60;
T=2*pi/w;
Ks1=1.0*(1e+9);Ks2=6.0*(1e+8);Asp1=(Ks1+Ks2)/2;Asp2=(Ks1-Ks2)/2;
Kr1=1.2*(1e+9);Kr2=8.0*(1e+8);Arp1=(Kr1+Kr2)/2;Arp2=(Kr1-Kr2)/2;
%?????1/d??
d=100000;
totaltime=0.008;
tt=totaltime*d+2;
%totaltime?????????tt????????????????d???
u=zeros(N,tt);
du=zeros(N,tt);
ddu=zeros(N,tt);
f1=zeros(tt);
f2=zeros(tt);
Gs1=zeros(tt);
Gr1=zeros(tt);
u(:,1)=[10;0;0];
du(:,1)=[0;0;0];
ddu(:,1)=[0;0;0];
%
%nuwmark参数
b=0.25;   r=0.5;  dt=1/d;
a0=1/(b*(dt)^2); a1=r/(b*dt);
a2=1/(b*dt);     a3=1/(2*b)-1;
a4=r/b-1;        a5=0.5*dt*((r/b)-2);
a6=dt*(1-r);     a7=r*dt;
f=zeros(N,tt);
Ef=zeros(N,tt);
t1=1;
%?t1??????????????
TD=90;
TL=306;
for t=0:1/d:totaltime
      %??detat=0.1???t=0:0.1:totalt
    f(:,t1)=[TD;0;TL];
    t1=t1+1;
end
t1=1;
for t=0:1/d:totaltime   
ks1=Asp1+Asp2*square(w*t,(e1-1)*100);
%ks2=Asp1+Asp2*square(w*(t+z1/3*T),(e1-1)*100);
%ks3=Asp1+Asp2*square(w*(t+z1*2/3*T),(e1-1)*100);
kr1=Arp1+Arp2*square(w*t,(e2-1)*100);
%kr2=Arp1+Arp2*square(w*(t-z3/3*T),(e2-1)*100);
%kr3=Arp1+Arp2*square(w*(t-z3*2/3*T),(e2-1)*100);
kss=0.0008;   krr=0.0008;
K=[kss+ks1,ks1,0;ks1,(ks1+kr1),-kr1;0,-kr1,kr1+krr]
%EK?????????
    EK=K+a0*M;
        %??t+detat????????
    Ef(:,t1+1)=f(:,t1+1)+M*(a0*u(:,t1)+a2*du(:,t1)+a3*ddu(:,t1));
    %??t+detat??????
    u(:,t1+1)=inv(EK)*Ef(:,t1+1);
    %??t+detat??????????
    ddu(:,t1+1)=a0*(u(:,t1+1)-u(:,t1))-a2*du(:,t1)-a3*ddu(:,t1);
    du(:,t1+1)=du(:,t1)+a6*ddu(:,t1)+a7*ddu(:,t1+1);
   f1(t1)=ks1*(u(1,t1)+u(2,t1));f2(t1)=kr1*(u(2,t1)-u(3,t1));
   Gs1(t1)=rb1*ks1*f1(t1)/TD;  Gr1(t1)=rb3*kr1*f2(t1)/TL;
    t1=t1+1;
end
t=0:1/d:totaltime+1/d;
subplot(2,2,1),plot(t,u(1,:)),xlabel('t'),ylabel('u');
%subplot(2,2,2),plot(t,du),xlabel('t'),ylabel('du');
%subplot(2,2,3),plot(t,ddu),xlabel('t'),ylabel('ddu');
subplot(2,2,4),plot(t,Gs1),xlabel('t'),ylabel('R1');
%plot(t,Ef(1,1:tt));hold on
%plot(t,Gs1); hold on
%plot(t,f1); hold on
%plot(t,u(1,:)); hold on
%plot(t,du(1,:)); hold on
%plot(t,ddu(1,:));
grid on
clear

redstonegu

上面是我花了两天时间编的纯扭转模型，各位帮看看对不对？我看论文上的u值和动载荷系数R都是>0的，我的不是。另外还想请问接下来该如何分析，刚入门，也不知道该做些什么是有用的，望各位指点。

vib

什莫叫做“齿轮变刚度”？楼主介绍一下，或许能有更多的人参与进来讨论

讨论求知

请问你那求齿轮变刚度是参考什么编的，有什么依据吗

无水1324

齿轮变刚度这里面的意思就，齿轮系统中的时变啮合刚度，就是接触过程刚度是随着位置（或者说是时间）变化的

redstonegu

这里我使用矩形波作为齿轮的变刚度没考虑间隙误差，初始位移，速度，角速度都为0，初始转矩90，负载转矩0，sun 轮无支撑，ring轮有rus=8000的支撑。结果出来挺奇怪的

redstonegu

clc;
N=5;
ha=1;ns=20;
z1=40;z2=33;z3=106;m=0.003;a=22*pi/180; pb=pi*m*cos(a);
ra1=m*z1/2+ha*m;ra2=m*z2/2+ha*m;ra3=m*z3/2-ha*m;
rb1=m*z1*cos(a)/2; rb2=m*z2*cos(a)/2; rb3=m*z3*cos(a)/2;
%mass
ms=4.75;m1=2.88;m2=2.88;m3=2.88;mr=16.92;
%j1=m1*rb1^2;j2=m2*rb2^2;j3=m3*rb3^2;
M=[ms,0,0,0,0;
   0,m1,0,0,0;
   0,0,m2,0,0;
   0,0,0,m3,0;
   0,0,0,0,mr];
%%%刚度
l1=sqrt(ra1^2-rb1^2)+sqrt(ra2^2-rb2^2)-(m*z1/2+m*z2/2)*sin(a);
l2=sqrt(ra2^2-rb2^2)-(sqrt(ra3^2-rb3^2)-(m*z3/2-m*z2/2)*sin(a));
e1=l1/pb;e2=l2/pb;
w=2*pi*ns*z1;T=2*pi/w;
Ks1=1.0e8;Ks2=6.0e7;Asp1=(Ks1+Ks2)/2;Asp2=(Ks1-Ks2)/2;
AKs=Ks1*(e1-1)+Ks2*(2-e1);
Kr1=1.2e8;Kr2=8e7;Arp1=(Kr1+Kr2)/2;Arp2=(Kr1-Kr2)/2;
AKr=Kr1*(e2-1)+Kr2*(2-e2);
%C阻尼
c0=0.07
cs=2*c0*square(AKs/(1/ms+1/m1));
cr=2*c0*square(AKr/(1/mr+1/m1));
cs1=cs;cs2=cs;cs3=cs;
cr1=cr;cr2=cr;cr3=cr;
C=[cs1+cs2+cs3,     cs1   ,     cs2   ,   cs3    ,      0;
      cs1     ,  (cs1+cr1),      0    ,     0    ,   -cr1;
      cs2     ,      0    , (cs1+cr2) ,     0    ,   -cr2;
      cs3     ,      0    ,      0    ,(cs1+cr3) ,   -cr3;
       0      ,     -cr1  ,   -cr2    ,   -cr3   ,cr1+cr2+cr3]
%step
d=40000;
totaltime=0.01;
tt=totaltime*d+2;
%totaltime?????????tt????????????????d???
u=zeros(N,tt);du=zeros(N,tt);ddu=zeros(N,tt);
fs1=zeros(1,tt);fs2=zeros(1,tt);fs3=zeros(1,tt);
fr1=zeros(1,tt);fr2=zeros(1,tt);fr3=zeros(1,tt);
Gs1=zeros(1,tt);Gr1=zeros(1,tt);
u(:,1)=[0;0;0;0;0];
du(:,1)=[0;0;0;0;0];
ddu(:,1)=[0;0;0;0;0];
%nuwmark参
b=0.25;   r=0.5;  dt=1/d;
a0=1/(b*(dt)^2); a1=r/(b*dt);
a2=1/(b*dt);     a3=1/(2*b)-1;
a4=r/b-1;        a5=0.5*dt*((r/b)-2);
a6=dt*(1-r);     a7=r*dt;
f=zeros(N,tt);
Ef=zeros(N,tt);
t1=1;
%force
TD=90;
TL=0;
for t=0:1/d:totaltime
      %??detat=0.1???t=0:0.1:totalt
    f(:,t1)=[TD;0;0;0;TL];
    t1=t1+1;
end
t1=1;
for t=0:1/d:totaltime   
ks1=Asp1+Asp2*square(w*t,(e1-1)*100);
ks2=Asp1+Asp2*square(w*(t+z1/3*T),(e1-1)*100);
ks3=Asp1+Asp2*square(w*(t+z1*2/3*T),(e1-1)*100);
kr1=Arp1+Arp2*square(w*t,(e2-1)*100);
kr2=Arp1+Arp2*square(w*(t-z3/3*T),(e2-1)*100);
kr3=Arp1+Arp2*square(w*(t-z3*2/3*T),(e2-1)*100);
ksu=0;   kru=8000;
K=[ks1+ks2+ks3+ksu,     ks1   ,   ks2    ,  ks3     ,   0;
         ks1      , (ks1+kr1) ,    0     ,   0      ,-kr1;
         ks2      ,     0     ,(ks1+kr2) ,   0      ,-kr2;
         ks3      ,     0     ,    0     , (ks1+kr3),-kr3;
          0       ,    -kr1   ,   -kr2   ,   -kr3   ,kr1+kr2+kr3+kru]
%EK?????????
    EK=K+a0*M+a1*C;
        %??t+detat????????
    Ef(:,t1+1)=f(:,t1+1)+M*(a0*u(:,t1)+a2*du(:,t1)+a3*ddu(:,t1))...
        +C*(a1*u(:,t1)+a4*du(:,t1)+a5*ddu(:,t1));
    %??t+detat??????
    u(:,t1+1)=inv(EK)*Ef(:,t1+1);
    %??t+detat??????????
    ddu(:,t1+1)=a0*(u(:,t1+1)-u(:,t1))-a2*du(:,t1)-a3*ddu(:,t1);
    du(:,t1+1)=du(:,t1)+a6*ddu(:,t1)+a7*ddu(:,t1+1);
    fs1(t1)=ks1*(u(1,t1)+u(2,t1)); fs2(t1)=ks2*(u(1,t1)+u(3,t1));  fs3(t1)=ks1*(u(1,t1)+u(4,t1));
    fr1(t1)=kr1*(u(5,t1)-u(2,t1));fr2(t1)=kr1*(u(5,t1)-u(3,t1));  fr3(t1)=kr1*(u(5,t1)-u(4,t1));
    Gs1(t1)=3*rb1*max(max(fs1(t1),fs2(t1)),fs3(t1))/TD;
   % Gr1(t1)=3*rb3*max(max(fr1(t1),fr2(t1)),fr3(t1))/TL;
    t1=t1+1;
   
end
t=0:1/d:totaltime+1/d;
subplot(2,2,1),plot(t/T,u(1,:)),xlabel('t/T'),ylabel('u');
subplot(2,2,2),plot(t/T,du(1,:)),xlabel('t/T'),ylabel('du');
subplot(2,2,3),plot(t/T,ddu(1,:)),xlabel('t/T'),ylabel('ddu');
subplot(2,2,4),plot(t/T,Gs1),xlabel('t/T'),ylabel('R1');
grid;
clear

redstonegu

结果如下：我的结果好像不稳定，也不像是振动

redstonegu

无水老师给指点一下吧

无水1324

我也在调试，但是我就是找不到原因哈，
还有方程的求解为什么不用ode呢？

redstonegu

谢谢无水老师帮忙。
我这是老师指定让用newmark方法的, 我刚开始做行星轮动力学研究，今天让老师看了一下，老师说我的初始条件不对，Ks*X=F,求一下初始条件x0,再带入到方程中。也就是先求出静态解。
我今天的程序，又调了调，
附件中的（1）图是当 x0=[[0.000006;0;0;0;0];这个初始条件是我随便给的.
附件中的（2）图是当x0=[0;0;0;0;0];
附件中（3）是程序。
在调试的过程中发现：转速ns，初始x0 ，和dt 改变时，结果相差很大。有时完全不能理解。
                                   还有就是结果应是周期的我的不是。

zhj0231984

有可能是 初始的 齿轮系统参数有不匹配吧.
请问楼主 你的初始参数 是鉴戒 哪个参考文献的啊？
我想找一组 匹配的完整的齿轮系统的 初始参数  却一直没找到。。

sjxu100

讨论的人好少啊

VibInfo

sjxu100 发表于 2015-11-2 10:18
讨论的人好少啊

要碰上对的人才能讨论起来

sjxu100

VibInfo 发表于 2015-11-2 11:02
要碰上对的人才能讨论起来

恩，刚开始学齿轮，这个帖子对我帮助比较大，但是有些问题不太懂，看这个发帖时间估计作者也不会再看了