齿轮啮合非线性动力学

Seraphena

各位高手，本人最近在做齿轮非线性动力学分析，参照硕士论文《斜齿轮多间隙非线性耦合系统动力学研究》中的10自由度弯-扭-轴-摆耦合动力学模型编写如下程序，但运行结果与原文有出入，而且扭摆振动不收敛，还望各位高手赐教，指出程序的问题所在

function dx=fun(tao,x)

global w

% 齿轮系统参数
mn=6; %法向模数，mm
B=70; %齿宽，mm
belt=17; %螺旋角，°
alpha=25; %法向压力角，°
z_1=34; %小齿轮齿数
z_2=84; %大齿轮齿数
m_1=15.47; %小齿轮质量，kg
m_2=91.24; %大齿轮质量，kg
J_1=74650; %小齿轮转动惯量，kg*mm^2
J_2=2596800; %大齿轮转动惯量，kg*mm^2
% b1=0.107/2; %主动轮轴承径向游隙，mm
% b2=0.075/2; %主动轮轴承轴箱游隙，mm
% b3=0.148/2; %被动轮轴承径向游隙，mm
% b4=0;       %被动轮轴承径向游隙，mm
b_5=0.025/2; %齿面侧隙，mm
e_a=108/1e3;   %齿轮传动误差幅值，mm
k_m=1.19e6; %平均啮合刚度，N/mm
s=0.1;    %刚度变化范围0.9*k_m~1.1*k_m
n=2000;   %主动轮转速，rad/min

d_1=mn*z_1/cosd(belt); %分度圆直径,mm
d_2=mn*z_2/cosd(belt);
r_1=d_1*cosd(alpha)/2; %基圆半径,mm
r_2=d_2*cosd(alpha)/2;

m_e=J_1*J_2/(J_1*r_2^2+J_2*r_1^2); %等效质量,kg
ksi_g=0.08; %齿轮啮合阻尼比，0.03~0.17
c_m=2*ksi_g*sqrt(k_m*m_e); %啮合阻尼
w_n=sqrt(k_m/m_e); %固有频率，rad/s
w_eh=2*pi*z_1*n/60; %齿轮啮合频率,rad/s
% T_eh=2*pi/w_eh; %齿轮啮合周期
w=w_eh/w_n; %无量纲频率

l1=74; %小齿轮支撑轴承距齿轮质心距离
l2=78.5; %大齿轮支撑轴承距齿轮质心距离

T1=10611; %输入轴转矩，N*m
T2=60350; %输出轴转矩，N*m

% 支撑轴承刚度及阻尼参数
k_1x=5.968e6; %支撑轴承径向刚度，N/mm
k_1y=5.968e6; %支撑轴承横向刚度，N/mm
k_1z=5.968e6; %支撑轴承轴向刚度，N/mm
k_2x=5.968e6;
k_2y=5.968e6;
k_2z=5.968e6;
c_1x=0; %%小齿轮轴支撑轴承径向阻尼，设阻尼比为0，则阻尼为0
c_1y=0;
c_1z=0;
c_2x=0;
c_2y=0;
c_2z=0;

w_1x=sqrt(k_1x/m_1);
w_1y=sqrt(k_1y/m_1);
w_1z=sqrt(k_1z/m_1);
w_2x=sqrt(k_2x/m_2);
w_2y=sqrt(k_2y/m_2);
w_2z=sqrt(k_2z/m_2);

% 微分方程参数
a11=c_1x/(m_1*w_n);
a12=(w_1x/w_n)^2;
a13=m_e*sind(alpha)*(1+s*cosd(w*tao))/m_1;
a14=c_m*sind(alpha)/(m_1*w_n);
a21=c_1y/(m_1*w_n);
a22=(w_1y/w_n)^2;
a23=m_e*cosd(alpha)*cosd(belt)*(1+s*cosd(w*tao))/m_1;
a24=c_m*cosd(alpha)*cosd(belt)/(m_1*w_n);
a31=c_1z/(m_1*w_n);
a32=(w_1z/w_n)^2;
a33=m_e*cosd(alpha)*sind(belt)*(1+s*cosd(w*tao))/m_1;
a34=c_m*cosd(alpha)*sind(belt)/(m_1*w_n);
a41=2*(l1+l1)/r_1*a11;
a42=2*(l1+l1)/r_1*a12;
a43=4*a33;
a44=4*a34;
a51=2*a23;
a52=2*a24;

a61=c_2x/(m_2*w_n);
a62=(w_2x/w_n)^2;
a63=-m_e*sind(alpha)*(1+s*cosd(w*tao))/m_2;
a64=-c_m*sind(alpha)/(m_2*w_n);
a71=c_2y/(m_2*w_n);
a72=(w_2y/w_n)^2;
a73=-m_e*cosd(alpha)*cosd(belt)*(1+s*cosd(w*tao))/m_2;
a74=-c_m*cosd(alpha)*cosd(belt)/(m_2*w_n);
a81=c_2z/(m_2*w_n);
a82=(w_2z/w_n)^2;
a83=-m_e*cosd(alpha)*sind(belt)*(1+s*cosd(w*tao))/m_2;
a84=-c_m*cosd(alpha)*sind(belt)/(m_2*w_n);
a91=2*(l2+l2)/r_2*a61;
a92=2*(l2+l2)/r_2*a62;
a93=-4*a83;
a94=-4*a84;
a101=-2*a73;
a102=-2*a74;

g1=-2*T1*1e3/(b_5*m_1*r_1*w_n^2);
g2=2*T2*1e3/(b_5*m_2*r_2*w_n^2);

%间隙非线性函数（用3次函数拟合间隙）
fy1=0.0107*x(1)^3-0.1173*x(1);         %f(p1)
fy2=0.0107*x(3)^3-0.1173*x(3);        %f(p2)
fy3=0.0031*x(5)^3-0.0568*x(5);        %f(p3)
fy6=0.0215*x(11)^3-0.1143*x(11);    %f(p6)
fy7=0.0215*x(13)^3-0.1143*x(13);    %f(p7)
fy8=x(15); %f(p8)
fy11=0.0638*(x(1)-x(11)-e_a*sind(alpha)*sind(w*tao)/b_5)^3+...
    0.1737*(x(1)-x(11)-e_a*sind(alpha)*sind(w*tao)/b_5);     %f(p11)
fy12=0.0638*(x(3)-x(13)-e_a*cosd(alpha)*cosd(belt)*sind(w*tao)/b_5)^3+...
    0.1737*(x(3)-x(13)-e_a*cosd(alpha)*cosd(belt)*sind(w*tao)/b_5);     %f(p12)
fy13=0.0638*(x(5)-x(15)-e_a*cosd(alpha)*sind(belt)*sind(w*tao)/b_5)^3+...
    0.1737*(x(5)-x(15)-e_a*cosd(alpha)*sind(belt)*sind(w*tao)/b_5);      %f(p13)

fz1=x(2)-x(12)-w*e_a*sind(alpha)*cosd(w*tao)/b_5; %p11导数
fz2=x(4)-x(14)-w*e_a*cosd(alpha)*cosd(belt)*cosd(w*tao)/b_5; %p12导数
fz3=x(6)-x(16)-w*e_a*cosd(alpha)*sind(belt)*cosd(w*tao)/b_5; %p13导数


% df=zeros(20,1);%初始化列向量
dx(1)=x(2);
dx(2)=-a11*x(2)-a12*fy1-a13.*fy11-a14*fz1;
dx(3)=x(4);
dx(4)=-a21*x(4)-a22*fy2-a23.*fy12-a24*fz2;
dx(5)=x(6);
dx(6)=-a31*x(6)-a32*fy3-a33.*fy13-a34*fz3;
dx(7)=x(8);
dx(8)=-a41*x(8)-a42*x(7)-a43.*fy13-a44*fz3;
dx(9)=x(10);
dx(10)=-a51.*fy12-a52*fz2;
dx(11)=x(12);
dx(12)=-a61*x(12)-a62*fy6-a63.*fy11-a64*fz1;
dx(13)=x(14);
dx(14)=-a71*x(14)-a72*fy7-a73.*fy12-a74*fz2;
dx(15)=x(16);
dx(16)=-a81*x(16)-a82*fy8-a83.*fy13-a84*fz3;
dx(17)=x(18);
dx(18)=-a91*x(18)-a92*x(17)-a93.*fy13-a94*fz3;
dx(19)=x(20);
dx(20)=-a101.*fy12-a102*fz2;
dx=[dx(1);dx(2);dx(3);dx(4);dx(5);dx(6);dx(7);dx(8);dx(9);dx(10);dx(11);dx(12);dx(13);dx(14);dx(15);dx(16);dx(17);dx(18);dx(19);dx(20)];


主函数：
T=2*pi/w;
y0=zeros(20,1); %初始值
options=odeset('relTol',1e-3,'Maxstep',0.1);
m=100;  %每周期采样点数
[tt,x]=ode45(@fun,[0:T/m:400*T],y0,options);

运行结果原文运行结果

求各位高手不吝赐教

MVH

从程序本身上看不出什么问题
1. 详细检查相关的公式有没有搞错
2. 比对一下参数是否完全一致
3. 追踪一下啮合刚度等参数，看是否合理

Seraphena

MVH 发表于 2015-11-30 10:16
从程序本身上看不出什么问题
1. 详细检查相关的公式有没有搞错
2. 比对一下参数是否完全一致

我看到有的论文将扭转方向的振动方程进行合并，然后引入新的自由度，不知道这对扭转方向的收敛有什么影响。现在单独解算扭转方向的2自由度微分方程，没有进行合并，结果和前述一样仍不收敛

Seraphena

MVH 发表于 2015-11-30 10:16
从程序本身上看不出什么问题
1. 详细检查相关的公式有没有搞错
2. 比对一下参数是否完全一致

您还在吗，请问我用ode45进行求解，除位移输出外，怎样追踪刚度参数？@MVH

MVH

Seraphena 发表于 2015-12-1 15:47
我看到有的论文将扭转方向的振动方程进行合并，然后引入新的自由度，不知道这对扭转方向的收敛有什么影响 ...

国内的文章,对于其分析结果的正确性有的时候可以持一定的保留意见,另外我不知道原文阻尼是否明确给出来了,这个经常是文章中可以调整的参数,而且也和收敛性有很大关系

MVH

Seraphena 发表于 2015-12-2 10:55
您还在吗，请问我用ode45进行求解，除位移输出外，怎样追踪刚度参数？@MVH

可以考虑在odefun中增加一个全局变量作为输出

Seraphena

MVH 发表于 2015-12-3 08:28
国内的文章,对于其分析结果的正确性有的时候可以持一定的保留意见,另外我不知道原文阻尼是否明确给出来了 ...

原文给出了阻尼比和阻尼公式c_m=2*ksi_g*sqrt(k_m*r_1^2*r_2^2*J_1*J_2/(J_1*r_1^2+J_2*r_2^2))，而有的文章中阻尼c_m=2*ksi_g*sqrt(k_m*m_e)，请问应该选择哪一个，这两个公式算出的阻尼可能会相差2个数量级。如果从微分方程归一化来考虑，采用第一个公式算系数a14=c_m*cos(belt)/(m_1*w_n)的结果是无量纲的吗？

kiu

楼主你好，我最近也在做这方面的研究，请问你的是几自由度的啊？

wzx1993

我运行了怎么都不会出图的  ，求助啊各位
??? function dx=fun(tao,x)
    |
Error: Function definitions are not permitted in this context.

frogfish

wzx1993 发表于 2016-4-5 15:37
我运行了怎么都不会出图的  ，求助啊各位
??? function dx=fun(tao,x)
    |

function dx=fun(tao,x)
这部分的内容要单独定义成一个文件

建议找本matlab的基础方面的书籍，先了解一下matlab的基础知识吧

wzx1993

frogfish 发表于 2016-4-6 15:17
function dx=fun(tao,x)
这部分的内容要单独定义成一个文件

谢谢，有什么简单入门的书可以推荐吗？刚刚接触matlab不久。

wzx1993

??? Input argument "x" is undefined.

Error in ==> fun0 at 109
fy1=0.0107*x(1)^3-0.1173*x(1);  需要怎么定义x

wzx1993

请问谁有[非线性振动] 齿轮啮合非线性动力学这个论文的全部程序？急求帮助，谢谢大家！

desolate

wzx1993 发表于 2016-4-10 19:41
请问谁有[非线性振动] 齿轮啮合非线性动力学这个论文的全部程序？急求帮助，谢谢大家！

fun函数是个子程序，x是需要传递进来的变量

desolate

wzx1993 发表于 2016-4-10 19:41
请问谁有[非线性振动] 齿轮啮合非线性动力学这个论文的全部程序？急求帮助，谢谢大家！

自己想办法联系原文作者吧。这种求助一般都不会有结果的