螺旋锥齿轮动力学求助
本帖最后由 wzy123 于 2014-5-14 21:50 编辑按照王立华博士论文上的螺旋锥齿轮振动模型编写了一个动力学方程，采用的是newmark算法，但是算出来的结果不是收敛的，newmark算法没有问题，测试过四个自由度的计算，结果是正确的。后来又改成了算8个自由度的模型，结果依然不正确。
现在对程序怀疑有两点：
1.刚度阻尼的取值的数量级有问题，因此推倒出了的刚度阻尼矩阵有问题
2.动力学模型有问题（8自由度的模型问题概率应该比较小）
8自由度的MATLAB程序如
%%%%%%齿轮基本参数%%%%%
m=7.5;z1=27;z2=31;b=51;
rou1=7.8e-9;%齿轮密度
rou2=7.8e-9;
alpha=20/180*pi;%压力角
beta1=35/180*pi;%中点螺旋角
beta2=35/180*pi;
sigma=90/180*pi;%轴交角
%%%%%%%基本参数计算%%%%%%
delta2=atan(z2*sin(sigma)/(z1+z2*cos(sigma)));
delta1=sigma-delta2;%节锥角
d1=m*z1;d2=m*z2;%节圆直径
Re=d2/(2*sin(delta2));%外锥距
hk=1.7*m;%工作齿高
ht=1.888*m;%全齿高
hae2=0.46*m+0.39*(z1*cos(delta2)/(z2*cos(delta1)))*m;%齿轮2齿顶高
hae1=hk-hae2;
hfe1=ht-hae1;%齿根高
hfe2=ht-hae2;
c=ht-hk;%顶隙
theta_f1=atan(hfe1/Re);%齿根角
theta_f2=atan(hfe2/Re);
delta_f1=delta1-theta_f1;%根锥角
delta_f2=delta2-theta_f2;
theta_a1=theta_f2;
theta_a2=theta_f1;
delta_a1=delta1+theta_a1;%面角度
delta_a2=delta2+theta_a2;
de1=d1+2*hae1*cos(delta1);%外径
de2=d2+2*hae1*cos(delta2);
X_e1=Re*cos(delta1)-hae1*sin(delta1);
X_e2=Re*cos(delta2)-hae1*sin(delta2);
p=pi*m;%节距
r1=m*z1/2;r2=m*z2/2;
rp1=m*z1/2;rp2=m*z2/2;
%%%%%%质量系数%%%%%%
m1=rou1*pi*(d1/2)^2*b*cos(delta1);
m2=rou2*pi*(d2/2)^2*b*cos(delta2);
I1=1/2*m1*(d1/2)^2;
I2=1/2*m2*(d2/2)^2;
%%%%%%%阻尼刚度系数%%%%%%%%%
cx1=0.1;cy1=0.1;cz1=0.1;
cx2=0.1;cy2=0.1;cz2=0.1;
kx1=1e6;ky1=1e6;kz1=1e6;
kx2=1e6;ky2=1e6;kz2=1e6;
%%%%%%%%化简用系数%%%%
a1=cos(delta1)*sin(alpha);
a2=cos(delta2)*sin(alpha);
a3=cos(delta1)*cos(alpha)*sin(beta1)+sin(delta1)*sin(alpha);
a4=cos(delta2)*cos(alpha)*sin(beta2)+sin(delta2)*sin(alpha);
a5=sin(delta1)*cos(alpha)*sin(beta1)+cos(delta1)*sin(alpha);
a6=sin(delta2)*cos(alpha)*sin(beta2)+cos(delta2)*sin(alpha);
a7=sin(delta2)*cos(alpha)*sin(beta2)+cos(delta2)*sin(alpha);
a8=sin(delta2)*sin(alpha)-cos(delta2)*cos(alpha)*sin(beta2);
a9=cos(alpha)*cos(beta2);
T1=100000;
T2=T1*z1/z2;
km=1e6;
cm=0.1;
n=1000;
wn=n*z1/60;
T=1/wn;
dt=1/50*T;
tend=4*T;
V=;
M=diag(V);
K=[kx1+a1*a7*km-a3*a7*km-a5*a7*km-a5*a7*r1*km-a2*a7*kma4*a7*kma6*a7*km a6*a7*r2*km;
-a1*a8*km ky1+a3*a8*kma5*a8*kma5*a8*r1*km a2*a8*km-a4*a8*km-a6*a8*km-a6*a8*r2*km;
-a1*a9*km a3*a9*km kz1+a5*a9*kma5*a9*r1*kma2*a9*km-a4*a9*km-a6*a9*km-a6*a9*r2*km;
a1*a9*r1*km-a3*a9*r1*km-a5*a9*r1*km-a5*a9*r1*r1*km-a2*a9*r1*kma4*a9*r1*kma6*a9*r1*kma6*a9*r1*r2*km;
-a1*a7*km a3*a7*km a5*a7*kma5*a7*r1*kmkx2+a2*a7*km-a4*a7*km-a6*a7*km-a6*a7*r2*km;
a1*a8*km-a3*a8*km-a5*a8*km-a5*a8*r1*km -a2*a8*kmky2+a4*a8*km a6*a8*km a6*a8*r2*km;
a1*a9*km-a3*a9*km-a5*a9*km-a5*a9*r1*km-a2*a9*kma4*a9*km kz2+a6*a9*kma6*a9*r2*km;
-a1*a9*r2*kma3*a9*r2*kma5*a9*r2*kma5*a9*r1*r2*kma2*a9*r2*km-a4*a9*r2*km-a6*a9*r2*km-a6*a9*r2*r2*km];
C=[cx1+a1*a7*cm-a3*a7*cm-a5*a7*cm-a5*a7*r1*cm-a2*a7*cma4*a7*cma6*a7*cm a6*a7*r2*cm;
-a1*a8*cm cy1+a3*a8*cma5*a8*cma5*a8*r1*cm a2*a8*cm-a4*a8*cm-a6*a8*cm-a6*a8*r2*cm;
-a1*a9*cm a3*a9*cm cz1+a5*a9*cma5*a9*r1*cma2*a9*cm-a4*a9*cm-a6*a9*cm-a6*a9*r2*cm;
a1*a9*r1*cm-a3*a9*r1*cm-a5*a9*r1*cm-a5*a9*r1*r1*cm-a2*a9*r1*cma4*a9*r1*cma6*a9*r1*cma6*a9*r1*r2*cm;
-a1*a7*cm a3*a7*cm a5*a7*cma5*a7*r1*cmcx2+a2*a7*cm-a4*a7*cm-a6*a7*cm-a6*a7*r2*cm;
a1*a8*cm-a3*a8*cm-a5*a8*cm-a5*a8*r1*cm -a2*a8*cmcy2+a4*a8*cm a6*a8*cm a6*a8*r2*cm;
a1*a9*cm-a3*a9*cm-a5*a9*cm-a5*a9*r1*cm-a2*a9*cma4*a9*cm cz2+a6*a9*cma6*a9*r2*cm;
-a1*a9*r2*cma3*a9*r2*cma5*a9*r2*cma5*a9*r1*r2*cma2*a9*r2*cm-a4*a9*r2*cm-a6*a9*r2*cm-a6*a9*r2*r2*cm];
%%%%%%积分系数%%%%%%%
gama=0.5;
beta=0.25;
b0=1/beta/dt^2;
b1=gama/beta/dt;
b2=1/beta/dt;
b3=1/2/beta-1;
b4=gama/beta-1;
b5=dt/2*(gama/beta-2);
b6=dt*(1-gama);
b7=gama*dt;
%%%%%%%%%%矩阵规模%%%%%%%%
K1=K+b0*M+b1*C;
Nd=zeros(8,floor(tend/dt)+1);
Nv=zeros(8,floor(tend/dt)+1);
Na=zeros(8,floor(tend/dt)+1);
f=zeros(8,floor(tend/dt)+1);
%%%%%%%初值%%%%%%%%%%
en0=1e-3;
en2=(1e-3)*wn;
f(:,1)=[a7*km*en0+a7*cm*en2;-a8*km*en0-a8*cm*en2;-a9*km*en0-a9*cm*en2;T1+r1*a9*km*en0+r1*a9*cm*en2;
-a7*km*en0-a7*cm*en2;a8*km*en0+a8*cm*en2;a9*km*en0+a9*cm*en2;-T2-r2*a9*km*en0-r2*a9*cm*en2];%初始载荷
Nd(:,1)=0;%初始位移
Nv(:,1)=0;%初始速度
Na(:,1)=inv(M)*(f(:,1)-K*Nd(:,1)-C*Nv(:,1));%初始加速度
%%%%%动力学响应求解%%%%%%%
t=0:dt:tend;
for i=2:1:length(t)
en=(1e-3)*sin(wn*t(i));
en1=(1e-3)*wn*cos(wn*t(i));
f1=[a7*km*en+a7*cm*en1;-a8*km*en-a8*cm*en1;-a9*km*en-a9*cm*en1;T1+r1*a9*km*en+r1*a9*cm*en1;
-a7*km*en-a7*cm*en1;a8*km*en+a8*cm*en1;a9*km*en+a9*cm*en1;-T2-r2*a9*km*en-r2*a9*cm*en1];
f2=f1+M*(b0*Nd(:,i-1)+b2*Nv(:,i-1)+b3*Na(:,i-1))+C*(b1*Nd(:,i-1)+b4*Nv(:,i-1)+b5*Na(:,i-1));
Nd(:,i)=inv(K1)*f2;
Na(:,i)=b0*(Nd(:,i)-Nd(:,i-1))-b2*Nv(:,i-1)-b3*Na(:,i-1);
Nv(:,i)=Nv(:,i-1)+b6*Na(:,i-1)+b7*Na(:,i);
end
求各位指导下，毕设课题有点着急啊。。。
从程序上来看拟采用的是显式的newmark法
一般而言,显式方法具有计算过程简捷,所占计算内存小,计算效率高等优点,但时间积分步长的选取须受到稳定性条件的限制。
隐式法无论对线性还是非线性问题,其数值稳定性较好,有助于选择时间积分步长,但由于每步积分都需求解高阶线性代数方程组,特别是对于非线性问题还需重新计算〔C〕,〔K〕矩阵,重新进行三角分解,这对大型工程问题的分析带来很大困难,其计算工作量是十分惊人的,计算费用也是相当可观的。
所以建议调整步长看看吧，实在不行还是应该改成隐式的吧
犟牛 发表于 2014-5-23 23:34
从程序上来看拟采用的是显式的newmark法
一般而言,显式方法具有计算过程简捷,所占计算内存小,计算效率高 ...
不是说β 、γ 为按精度和稳定性要求调节公式特性的参数，一般 β 取值范围为0 ≤ β ≤0.5。当 γ ≥0.5β>0.25(0.5+γ^2)时 Newmark-β法是无条件稳定的。那么时间步长对稳定性收敛该如何调整
哥们你的问题解决了没有啊，我也遇到了这个问题
页:
[1]