我列了一对齿的动力学方程，不知道哪错了，大家帮忙看看
function dy=donglixue(t,y)beta=18; %螺旋角
p=560; %输入功率，kw
n1=4100; %输入转速，r/min
alfa=20; %齿形角
d1=0.2081895; %齿轮1分度圆直径，m
d2=0.5425545; %齿轮2分度圆直径，m
g=0.08; %啮合阻尼比
k=7850; %铁的密度，kg/m3
km=18.553*10^9; %平均刚度
i=2.6061; %传动比
z1=33; %齿轮1的齿数
b=0.065; %齿宽，m
b3=2.5*10^(-4); %侧隙
ea=1.08*10^(-4); %误差幅值
f=pi/2; %相位角
s=0.2; %啮合刚度幅值ka/km
k1=0.5*10^6; %齿轮1的y方向支承刚度
k2=0.5*10^6; %齿轮1的z方向支承刚度
k3=0.5*10^6; %齿轮2的y方向支承刚度
k4=0.5*10^6; %齿轮2的z方向支承刚度
c1=0; %齿轮1的y方向支承阻尼
c2=0; %齿轮1的z方向支承阻尼
c3=0; %齿轮2的y方向支承阻尼
c4=0; %齿轮2的z方向支承阻尼
f1=0; %齿轮1支承系统所受的横向外载荷
f2=0; %齿轮2支承系统所受的横向外载荷
f3=0; %齿轮1支承系统所受的径向外载荷
f4=0; %齿轮2支承系统所受的径向外载荷
at=atan(tan(alfa*pi/180)/cos(beta*pi/180));%端面齿形角
r1=d1*cos(at)/2; %齿轮1基圆半径，m
r2=d2*cos(at)/2; %齿轮2基圆半径，m
m1=pi*(d1/2)^2*b*k; %齿轮1的质量，kg
m2=pi*(d2/2)^2*b*k; %齿轮2的质量，kg
j1=m1*r1^2/2; %齿轮1的转动惯量
j2=m2*r2^2/2; %齿轮2的转动惯量
me=j1*j2/(j1*r2^2+j2*r1^2);%质量当量
t1=9.55*10^6*p/n1; %转矩，N.mm
n2=n1/i; %输出转速，r/min
t2=9.55*10^6*p/n2; %转矩，N.mm
wn=sqrt(km/me); %固有频率
wh=(n1/60)*pi*z1/30;%啮合频率
w=wh/wn;%激励频率
ch=2*g*sqrt(km*r1^2*r2^2*j1*j2/(r1^2*j1+r2^2*j2));%啮合阻尼
w1=sqrt(k1/m1);
w2=sqrt(k2/m1);
w3=sqrt(k3/m2);
w4=sqrt(k4/m2);
c11=c1/(2*m1*wn);
c13=ch*cos(beta*pi/180)/(2*m1*wn);
c22=c3/(2*m2*wn);
c23=ch*cos(beta*pi/180)/(2*m2*wn);
c34=c2/(2*m1*wn);
c36=ch*sin(beta*pi/180)/(2*m1*wn);
c45=c4/(2*m2*wn);
c46=ch*sin(beta*pi/180)/(2*m2*wn);
c53=c13;
c63=c23;
k11=w1^2/wn^2;
k13=(1+s*sin(w*t+f))*cos(beta*pi/180)*me/m1;
k22=w3^2/wn^2;
k23=(1+s*sin(w*t+f))*cos(beta*pi/180)*me/m2;
k34=w2^2/wn^2;
k36=(1+s*sin(w*t+f))*sin(beta*pi/180)*me/m1;
k45=w4^2/wn^2;
k46=(1+s*sin(w*t+f))*sin(beta*pi/180)*me/m2;
k53=k13;
k63=k23;
g1=f1/(m1*b3*wn^2);
g2=f2/(m2*b3*wn^2);
g3=f3/(m1*b3*wn^2);
g4=f4/(m2*b3*wn^2);
g5=t1*10^(-3)/(m1*r1*b3*wn^2);
g6=t2*10^(-3)/(m2*r2*b3*wn^2);
e=ea*sin(w*t+f)/b3; %无量纲时变啮合误差
ey=e*cos(beta*pi/180);%y方向的无量纲时变啮合误差
ez=e*sin(beta*pi/180);%z方向的无量纲时变啮合误差
y1=y(1)+y(9)-y(3)+y(11)-ey;
y2=y(5)-y(7)+(-y(1)+y(2)-y(9)-y(11))*tan(beta*pi/180)-ez;
y11=y(2)+y(10)-y(4)+y(12)-w*ea*cos(beta*pi/180)*cos(w*t+f)/(b3*wn);
y22=y(6)-y(8)+(-y(2)+y(4)-y(10)-y(12))*tan(beta*pi/180)-w*ea*sin(beta*pi/180)*cos(w*t+f)/(b3*wn);
if y1>1
fy1=y1-1;
else if abs(y1)<=1
fy1=0;
else
fy1=y1+1;
end
end
if y2>1
fy2=y2-1;
else if abs(y2)<=1
fy2=0;
else
fy2=y2+1;
end
end
dy=[y(2);g1-2*c11*y(2)-k11*y(1)-k13*fy1-2*c13*y11;
y(4);-g2-2*c22*y(4)-k22*y(3)+k23*fy1+2*c23*y11;
y(6);g3-2*c34*y(6)-k34*y(5)-k36*fy2-2*c36*y22;
y(8);-g4-2*c45*y(8)-k45*y(7)+k46*fy2+2*c46*y22;
y(10);g5-k53*fy1-2*c53*y11;
y(12);-g6+k63*fy1+2*c63*y11]; 程序表达的原公式最好也能贴出来，大家对比者比较容易看，直接看程序太费劲了 最好多点说明。
建议楼主!!
建议先把你的想法和实现的意图说明一下!要不然大家看起来费劲!!
页:
[1]