齿轮变刚度计算出现错误
本帖最后由 VibInfo 于 2015-11-2 11:01 编辑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)=;
du(:,1)=;
ddu(:,1)=;
%
%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)=;
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=
%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 上面是我花了两天时间编的纯扭转模型，各位帮看看对不对？我看论文上的u值和动载荷系数R都是>0的，我的不是。另外还想请问接下来该如何分析，刚入门，也不知道该做些什么是有用的，望各位指点。 什莫叫做“齿轮变刚度”？楼主介绍一下，或许能有更多的人参与进来讨论 请问你那求齿轮变刚度是参考什么编的，有什么依据吗
回复 板凳 vib 的帖子
齿轮变刚度这里面的意思就，齿轮系统中的时变啮合刚度，就是接触过程刚度是随着位置（或者说是时间）变化的 这里我使用矩形波作为齿轮的变刚度没考虑间隙误差，初始位移，速度，角速度都为0，初始转矩90，负载转矩0，sun 轮无支撑，ring轮有rus=8000的支撑。结果出来挺奇怪的这个程序更全面一些
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)=;
du(:,1)=;
ddu(:,1)=;
%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)=;
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 结果如下：我的结果好像不稳定，也不像是振动 无水老师给指点一下吧
回复 9楼 redstonegu 的帖子
我也在调试，但是我就是找不到原因哈，还有方程的求解为什么不用ode呢？ 谢谢无水老师帮忙。
我这是老师指定让用newmark方法的, 我刚开始做行星轮动力学研究，今天让老师看了一下，老师说我的初始条件不对，Ks*X=F,求一下初始条件x0,再带入到方程中。也就是先求出静态解。
我今天的程序，又调了调，
附件中的（1）图是当 x0=[;这个初始条件是我随便给的.
附件中的（2）图是当x0=;
附件中（3）是程序。
在调试的过程中发现：转速ns，初始x0 ，和dt 改变时，结果相差很大。有时完全不能理解。
还有就是结果应是周期的我的不是。 有可能是 初始的 齿轮系统参数有不匹配吧.
请问楼主 你的初始参数 是鉴戒 哪个参考文献的啊？
我想找一组 匹配的完整的齿轮系统的 初始参数却一直没找到。。 讨论的人好少啊 sjxu100 发表于 2015-11-2 10:18
讨论的人好少啊
要碰上对的人才能讨论起来 VibInfo 发表于 2015-11-2 11:02
要碰上对的人才能讨论起来
恩，刚开始学齿轮，这个帖子对我帮助比较大，但是有些问题不太懂，看这个发帖时间估计作者也不会再看了
页:
[1]
2