如何用ode45求10自由度齿轮扭振模型？

zhlei566

各位高高手，本人初学齿轮系统动力学，建立了一个考虑齿轮时变刚度的扭振分析模型，我用四五阶龙格库塔法（ode45）求解时，发现运算的时间很长，不知是什么问题，望各位高手不吝赐教。程序如下：
    扭振模型：

1. function doty=GearTorqueVibratory(t,y)
   
2. %转动惯量kg.m2
   
3. doty=zeros(10,2);
   
4. Jin=4.217823;
   
5. J1=0.008952085;
   
6. J2=0.190956155;
   
7. J3=0.020671561;
   
8. J4=0.766038806;
   
9. J5=0.090060;
   
10. J6=4.24048334;
    
11. J7=0.6044254;
    
12. J8=27.9071658;
    
13. Jout=58.764562;
    
14. %齿轮基圆半径m
    
15. R1=(63.171e-3)/2;
    
16. R2=(234.749e-3)/2;
    
17. R3=(97.732e-3)/2;
    
18. R4=(321.941e-3)/2;
    
19. R5=(137.755e-3)/2;
    
20. R6=(451.530e-3)/2;
    
21. R7=(213.959e-3)/2;
    
22. R8=(601.758e-3)/2;
    
23. %轴扭转刚度Nm.rad
    
24. k1=252851;
    
25. k2=1936800;
    
26. k3=6186712;
    
27. k4=33459824;
    
28. k5=10190124;
    
29. Tin=819.4;%输入扭矩N.m
    
30. Tout=91675;%输出轴扭矩N.m
    
31. % 传动轴阻尼计算
    
32. lanmda=0.01;
    
33. c1=2*lanmda*sqrt(k1/(1/Jin+1/J1));
    
34. c2=2*lanmda*sqrt(k2/(1/J2+1/J3));
    
35. c3=2*lanmda*sqrt(k3/(1/J4+1/J5));
    
36. c4=2*lanmda*sqrt(k4/(1/J6+1/J7));
    
37. c5=2*lanmda*sqrt(k5/(1/J8+1/Jout));
    
38. %扭振模型基本参数矩阵
    
39. J=[Jin 0 0 0 0 0 0 0 0 0;
    
40.     0 J1 0 0 0 0 0 0 0 0;
    
41.     0 0 J2 0 0 0 0 0 0 0;
    
42.     0 0 0 J3 0 0 0 0 0 0;
    
43.     0 0 0 0 J4 0 0 0 0 0;
    
44.     0 0 0 0 0 J5 0 0 0 0;
    
45.     0 0 0 0 0 0 J6 0 0 0;
    
46.     0 0 0 0 0 0 0 J7 0 0;
    
47.     0 0 0 0 0 0 0 0 J8 0;
    
48.     0 0 0 0 0 0 0 0 0 Jout];
    
49. T=[Tin 0 0 0 0 0 0 0 0 Tout]';
    
50. %齿轮时变啮合刚度N.B/m
    
51. k12=GT1(t)/R1^2*1e6;
    
52. k34=GT2(t)/R2^2*1e6;
    
53. k56=GT3(t)/R3^2*1e6;
    
54. k78=GT4(t)/R4^2*1e6;   
    
55. %齿轮副阻尼计算
    
56. lanmdag=0.05;
    
57. c12=2*lanmdag*sqrt(k12*R1^2*R2^2*J1*J2/(R1^2*J1+R2^2*J2));
    
58. c34=2*lanmdag*sqrt(k34*R3^2*R4^2*J3*J4/(R3^2*J3+R4^2*J4));
    
59. c56=2*lanmdag*sqrt(k56*R5^2*R6^2*J5*J6/(R5^2*J5+R6^2*J6));
    
60. c78=2*lanmdag*sqrt(k78*R7^2*R8^2*J7*J8/(R7^2*J7+R8^2*J8));
    
61. %刚度矩阵   
    
62. K=[k1 -k1 0 0 0 0 0 0 0 0;
    
63.     -k1 k1+R1^2*k12 -R1*R2*k12 0 0 0 0 0 0 0;
    
64.     0 -R1*R2*k12 k2+R2^2*k12 -k2 0 0 0 0 0 0;
    
65.     0 0 -k2 k2+R3^2*k34 -R3*R4*k34 0 0 0 0 0;
    
66.     0 0 0 -R3*R4*k34 k3+R4^2*k34 -k3 0 0 0 0;
    
67.     0 0 0 0 -k3 k3+R5^2*k56 -R5*R6*k56 0 0 0;
    
68.     0 0 0 0 0 -R5*R6*k56 k4+R6^2*k56 -k4 0 0;
    
69.     0 0 0 0 0 0 -k4 k4+R7^2*k78 -R7*R8*k78 0;
    
70.     0 0 0 0 0 0 0 -R7*R8*k78 k5+R8^2*k78 -k5;
    
71.     0 0 0 0 0 0 0 0 -k5 k5];
    
72. %阻尼矩阵
    
73. C=[c1 -c1 0 0 0 0 0 0 0 0;
    
74.     -c1 c1+R1^2*c12 -R1*R2*c12 0 0 0 0 0 0 0;
    
75.     0 -R1*R2*c12 c2+R2^2*c12 -c2 0 0 0 0 0 0;
    
76.     0 0 -c2 c2+R3^2*c34 -R3*R4*c34 0 0 0 0 0;
    
77.     0 0 0 -R3*R4*c34 c3+R4^2*c34 -c3 0 0 0 0;
    
78.     0 0 0 0 -c3 c3+R5^2*c56 -R5*R6*c56 0 0 0;
    
79.     0 0 0 0 0 -R5*R6*c56 c4+R6^2*c56 -c4 0 0;
    
80.     0 0 0 0 0 0 -c4 c4+R7^2*c78 -R7*R8*c78 0;
    
81.     0 0 0 0 0 0 0 -R7*R8*c78 c5+R8^2*c78 -c5;
    
82.     0 0 0 0 0 0 0 0 -c5 c5];   
    
83. %初始值
    
84. %激励力
    
85. P=[Tin 0 0 0 0 0 0 0 0 Tout]';
    
86. invJ=inv(J);
    
87. M=invJ*P
    
88. N=invJ*C
    
89. L=invJ*K
    
90. doty=[y(2)
    
91.      M(1,1)-N(1,1)*y(2)-N(1,2)*y(4)-L(1,1)*y(1)-L(1,2)*y(3)
    
92.      y(4)
    
93.      M(2,1)-N(2,1)*y(2)-N(2,2)*y(4)-N(2,3)*y(6)-L(2,1)*y(1)-L(2,2)*y(3)-L(2,3)*y(5)
    
94.      y(6)
    
95.      M(3,1)-N(3,2)*y(4)-N(3,3)*y(6)-N(3,4)*y(8)-L(3,2)*y(3)-L(3,3)*y(5)-L(3,4)*y(7)
    
96.      y(8)
    
97.      M(4,1)-N(4,3)*y(6)-N(4,4)*y(8)-N(4,5)*y(10)-L(4,3)*y(5)-L(4,4)*y(7)-L(4,5)*y(9)
    
98.      y(10)
    
99.      M(5,1)-N(5,4)*y(8)-N(5,5)*y(10)-N(5,6)*y(12)-L(5,4)*y(7)-L(5,5)*y(9)-L(5,6)*y(11)
    
100.      y(12)
     
101.      M(6,1)-N(6,5)*y(10)-N(6,6)*y(12)-N(6,7)*y(14)-L(6,5)*y(9)-L(6,6)*y(11)-L(6,7)*y(13)
     
102.      y(14)
     
103.      M(7,1)-N(7,6)*y(12)-N(7,7)*y(14)-N(7,8)*y(16)-L(7,6)*y(11)-L(7,7)*y(13)-L(7,8)*y(15)
     
104.      y(16)
     
105.      M(8,1)-N(8,7)*y(14)-N(8,8)*y(16)-N(8,9)*y(18)-L(8,7)*y(13)-L(8,8)*y(15)-L(8,9)*y(17)
     
106.      y(18)
     
107.      M(9,1)-N(9,8)*y(16)-N(9,9)*y(18)-N(9,10)*y(20)-L(9,8)*y(15)-L(9,9)*y(17)-L(9,10)*y(19)
     
108.      y(20)
     
109.      M(10,1)-N(10,9)*y(18)-N(10,10)*y(20)-L(10,9)*y(17)-L(10,10)*y(19)]
     
110. end

复制代码

齿轮副时变刚度见附件（GT1、GT2、GT3、GT4）

1. clc;
   
2. clear;
   
3. t_span=linspace(0,3,100);
   
4. Y0=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
   
5. [t,y]=ode45('GearTorqueVibratory',t_span,Y0);
   
6. plot(t,y)

复制代码

xiaoyaoyu3700

楼主你好，我是刚刚接触齿轮故障检测的新人，对于齿轮副的齿数比，是让其等于1好还是不为1好呢？各会有什么优缺点吗，希望您能指点一二，谢谢。

zhlei566

xiaoyaoyu3700 发表于 2014-5-10 20:40
楼主你好，我是刚刚接触齿轮故障检测的新人，对于齿轮副的齿数比，是让其等于1好还是不为1好呢？各会有什么 ...

常规设计，齿数比一般不为整数，这样的好处是齿轮在运转过程中每一对齿都有相互啮合的机会，进而磨损较为均匀，也不会产生啮合损伤的累积，这也是目前主流设计观点；随着硬齿面齿轮的广泛应用，齿面渗碳淬火后不易磨损，齿数互质也就失去了存在的意义，如SEW减速器现在开始用整数速比了，整数速比的好处是可以根据每个轮齿的齿距偏差进行选配，使固定的几对齿啮合，提高啮合的平稳性。

xiaoyaoyu3700

zhlei566 发表于 2014-5-11 10:10
常规设计，齿数比一般不为整数，这样的好处是齿轮在运转过程中每一对齿都有相互啮合的机会，进而磨损较为 ...

嗯嗯，楼主说的很专业，我现在就是在纠结搭建两个齿轮啮合的试验台时，是让两个齿轮的齿数相同好呢？还是不同好呢？如果相同的话那齿轮的振动信号中就会只有一个转频及其倍频分量了，貌似分析起来也比较简单，但是好像和现实应用不太相符，因为感觉现实中相互啮合的齿轮齿数比一般是不相同的，也就是一大一小的。那您认为就齿轮故障诊断这方面来说，试验台的齿数比有什么限制和要求吗？谢谢

zhlei566

xiaoyaoyu3700 发表于 2014-5-11 20:00
嗯嗯，楼主说的很专业，我现在就是在纠结搭建两个齿轮啮合的试验台时，是让两个齿轮的齿数相同好呢？还是 ...

既然要搭建试验台，考虑代表性和常规性，建议齿数互质。为便于加载同时降低测功机或磁粉加载器的采购成本，建议速比不易过大，结合你的需要定义即可。

xiaoyaoyu3700

zhlei566 发表于 2014-5-11 21:09
既然要搭建试验台，考虑代表性和常规性，建议齿数互质。为便于加载同时降低测功机或磁粉加载器的采购成本 ...

嗯嗯，非常感谢您的建议，希望以后能对多多跟您交流

acai1238

怎么没有人回答

meiyongyuandeze

如果系统不存在强非线性项，只是个线性系统，如果长时间不收敛的话可能是系统出现发散运动，可能是系统参数问题！如果存在非线性环节，可能是积分步长或者选择的积分方法不对！

luotao1014

正在学习

Zeng_0902

学习了

sjxu100

lz问题解决了么

Edinburgh

sjxu100 发表于 2015-11-9 21:21
lz问题解决了么

个人理解楼主的问题其实不是很明确
长时间运行不知道是否能够得到结果
很可能是由于模型的收敛性差造成的

另外这类问题建议也可以尝试一下其他算法
比如newmark算法等

sjxu100

Edinburgh 发表于 2015-11-10 09:15
个人理解楼主的问题其实不是很明确
长时间运行不知道是否能够得到结果
很可能是由于模型的收敛性差造成 ...

做齿轮的好多都在用newmark算法，比ode算法的好在哪里呢？能解决高自由度数值算法依赖初值，不容易收敛的问题么

kkkttt

sjxu100 发表于 2015-11-10 19:48
做齿轮的好多都在用newmark算法，比ode算法的好在哪里呢？能解决高自由度数值算法依赖初值，不容易收敛的 ...

相对于rk算法，newmark算法虽然精度略低，但是计算效率更好
更重要的是newmark算法稳定性较好，在一定条件下可以做到无条件稳定（稳定性和步长没关系）

Seraphena

正在学习中，希望能有帮助