mesh_stiffness_ssssdd_齿轮_行星齿轮_动力学_

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%石川法求齿轮外啮合刚度%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear; clc; z1=27;%太阳轮 z2=33;%行星轮 m=18; b1=500; b2=500; ha=m; c=0.25*m; d1=m*z1;d2=m*z2; r1=m*z1/2;r2=m*z2/2;%分度圆半径 hf1=1.25*m; hf2=1.25*m; alpha=20*pi/180;%分度圆压力角 alpha_x=20*pi/180;%啮合角 invalpha=tan(alpha)-alpha; db1=d1*cos(alpha); db2=d2*cos(alpha); rb1=db1/2; rb2=db2/2;%基圆半径 rf1=r1-hf1; rf2=r2-hf2;%齿根圆半径 da1=d1+2*ha; da2=d2+2*ha;%齿顶圆 ra1=da1/2; ra2=da2/2;%齿顶圆半径 alpha_a1=acos(rb1/ra1); alpha_a2=acos(rb2/ra2);%齿顶圆压力角 alpha_f1=acos(rb1/rf1); alpha_f2=acos(rb2/rf2);%齿根圆压力角 s=pi*m/2;%分度圆弧齿厚 e=s;%分度圆齿槽宽 sk1=ra1*(s/r1+2*((tan(alpha)-alpha)-(tan(alpha_a1)-alpha_a1))); sk2=ra2*(s/r2+2*((tan(alpha)-alpha)-(tan(alpha_a2)-alpha_a2)));%齿顶圆齿厚 PB1=r1*cos(alpha)*(tan(alpha_a1)-tan(alpha)); PB2=r2*cos(alpha)*(tan(alpha_a2)-tan(alpha)); B1B2=PB1+PB2;%P分度圆相切点 Pb=pi*m*cos(alpha);%基圆齿距 Epsilona=B1B2/Pb; %N1N2为理论啮合线，B1B2为理论啮合线 N1B1=sqrt(ra1^2-rb1^2);%B1：齿顶极限啮合点 N1B2=N1B1-B1B2; % N1B2=-N1B1+B1B2; N2B2=sqrt(ra2^2-rb2^2);%从齿轮圆心做啮合线垂线，与啮合线交点为N2 N2B1=N2B2-B1B2; % N2B1=-N2B2+B1B2; rF2=sqrt(rb2^2+N2B1^2); rF1=sqrt(rb1^2+N1B2^2);%有效齿根圆半径,有效齿根圆半径就是齿轮啮合时渐开线刚进入啮合的那点对应的园半径. alpha_F1=acos(rb1/rF1); alpha_F2=acos(rb2/rF2); sf1=2*rF1*sin(pi/2/z1+invalpha-tan(alpha_F1)+alpha_F1);%齿根圆齿厚 sf2=2*rF2*sin(pi/2/z2+invalpha-tan(alpha_F2)+alpha_F2); h1=sqrt(ra1^2-(sk1/2)^2)-sqrt(rf1^2-(sf1/2)^2);%齿全高 h2=sqrt(ra2^2-(sk2/2)^2)-sqrt(rf2^2-(sf2/2)^2); hr1=sqrt(rF1^2-(sf1/2)^2)-sqrt(rf1^2-(sf1/2)^2);%矩形高 hr2=sqrt(rF2^2-(sf2/2)^2)-sqrt(rf2^2-(sf2/2)^2); hi1=(h1*sf1-hr1*sk1)/(sf1-sk1);%辅助尺寸 hi2=(h2*sf2-hr2*sk2)/(sf2-sk2); %B2CDB1为实际啮合取区 N2C=N2B1+Pb; B2C=B1B2-Pb;%双齿啮合区 CD=Pb-B2C;%单齿啮合区 N1C=N1B1-Pb; Fn=670000/(r1*cos(20/180*pi));%外力3333.33 E=2.05e+008;%弹性模量 v=0.3;%泊松比 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%将双齿啮合区等分100份，累加啮合线位移迭代计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%一个啮合周期在啮合线位移上是一个齿距Pb,Pb由B2C和CD组成%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% n=100; step=B2C/n;%双齿啮合区位移等分，B2C为双齿啮合区 nz1=3000;%齿轮1转速r/min Tz=60/z1/nz1;%循环的周期,单个齿从啮合进入到啮出所用时间（线位移上是一个齿距Pb）。 t1=B2C/Pb*Tz;%双齿啮合区所占时间 t2=CD/Pb*Tz;%单齿啮合区所占时间 step2=t1/n; step4=t2/n;%时间等分 %双齿啮合区设一啮合为i点，一啮合点为j点。 for i=1:n x(i)=i*step;%啮合线位移量 tt(i)=i*step2;%双齿啮合时，时间量变化 xx(i)=Pb+i*step;%双齿啮合时另一个齿的啮合位移量，一个齿距的差值 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%双齿啮合区i啮合点公式中具体参数的计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% N1Bi(i)=N1B2+i*step;%理论啮合线上的位移 O1Bi(i)=sqrt(N1Bi(i)*N1Bi(i)+rb1^2); ai1(i)=acos(rb1/O1Bi(i)); %invalpha=tan(alpha)-alpha;%渐开线函数 gamai1(i)=pi/2/z1+tan(alpha)-alpha-tan(ai1(i))+ai1(i); miui1(i)=ai1(i)-gamai1(i);%载荷作用角 rxi1(i)=O1Bi(i);%载荷作用点与齿轮中心点的距离，啮合点到齿轮中心的距离 hxi1(i)=rxi1(i)*cos(alpha_x-miui1(i))-sqrt(rf1^2-(sf1/2)^2);%c齿轮1节圆高 N2Bi(i)=N2B2-i*step; O2Bi(i)=sqrt(N2Bi(i)^2+rb2^2); ai2(i)=acos(rb2/O2Bi(i)); gamai2(i)=pi/2/z2+tan(alpha)-alpha-tan(ai2(i))+ai2(i); miui2(i)=ai2(i)-gamai2(i); rxi2(i)=O2Bi(i); hxi2(i)=rxi2(i)*cos(alpha_x-miui2(i))-sqrt(rf2^2-(sf2/2)^2);%c齿轮2节圆高 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%双齿啮合区j啮合点公式中具体参数的计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% N1Bj(i)=N1Bi(i)+Pb; O1Bj(i)=sqrt(N1Bj(i)^2+rb1^2); aj1(i)=acos(rb1/O1Bj(i)); gamaj1(i)=pi/2/z1+tan(alpha)-alpha-tan(aj1(i))+aj1(i); miuj1(i)=aj1(i)-gamaj1(i); rxj1(i)=O1Bj(i); hxj1(i)=rxj1(i)*cos(gamaj1(i))-sqrt(rf1^2-(sf1/2)^2); N2Bj(i)=N2Bi(i)-Pb; O2Bj(i)=sqrt(N2Bj(i)^2+rb2^2); aj2(i)=acos(rb2/O2Bj(i)); gamaj2(i)=pi/2/z2+tan(alpha)-alpha-tan(aj2(i))+aj2(i); miuj2(i)=aj2(i)-gamaj2(i); rxj2(i)=O2Bj(i); hxj2(i)=rxj2(i)*cos(gamaj2(i))-sqrt(rf2^2-(sf2/2)^2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%变形部分计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% sigmabri1(i)=12*Fn*cos(miui1(i))^2*(hxi1(i)*hr1*(hxi1(i)-hr1)+hxi1(i)^3/3)/b1/E/sf1^3;%长方形部分变形 sigmabri2(i)=12*Fn*cos(miui2(i))^2*(hxi2(i)*hr2*(hxi2(i)-hr2)+hxi2(i)^3/3)/b2/E/sf2^3; sigmabti1(i)=6*Fn*cos(miui1(i))^2*((hi1-hxi1(i))/(hi1-hr1)*(4-(hi1-hxi1(i))/(hi1-hr1))-2*log((hi1-hxi1(i))/(hi1-hr1))-3)*(hi1-hr1)^3/b1/E/sf1^3;%梯形部分变形量 sigmabti2(i)=6*Fn*cos(miui2(i))^2*((hi2-hxi2(i))/(hi2-hr2)*(4-(hi2-hxi2(i))/(hi2-hr2))-2*log((hi2-hxi2(i))/(hi2-hr2))-3)*(hi2-hr2)^3/b2/E/sf2^3; sigmasi1(i)=2*(1+v)*Fn*cos(miui1(i))^2*(hr1+(hi1-hr1)*log((hi1-hr1)/(hi1-hxi1(i))))/b1/E/sf1;%剪力产生变形量 sigmasi2(i)=2*(1+v)*Fn*cos(miui2(i))^2*(hr2+(hi2-hr2)*log((hi2-hr2)/(hi2-hxi2(i))))/b2/E/sf2; sigmagi1(i)=24*Fn*hxi1(i)*cos(miui1(i))^2/pi/b1/E/sf1^2;%基础部分倾斜产生的变形量 sigmagi2(i)=24*Fn*hxi2(i)*cos(miui2(i))^2/pi/b2/E/sf2^2; sigmap=4*Fn*(1-v^2)/pi/b2/E;%齿面接触部分变形量 sigmabrj1(i)=12*Fn*cos(miuj1(i))^2*(hxj1(i)*hr1*(hxj1(i)-hr1)+hxj1(i)^3/3)/b1/E/sf1^3; sigmabrj2(i)=12*Fn*cos(miuj2(i))^2*(hxj2(i)*hr2*(hxj2(i)-hr2)+hxj2(i)^3/3)/b2/E/sf2^3; sigmabtj1(i)=6*Fn*cos(miuj1(i))^2*((hi1-hxj1(i))/(hi1-hr1)*(4-(hi1-hxj1(i))/(hi1-hr1))-2*log((hi1-hxj1(i))/(hi1-hr1))-3)*(hi1-hr1)^3/b1/E/sf1^3; sigmabtj2(i)=6*Fn*cos(miuj2(i))^2*((hi2-hxj2(i))/(hi2-hr2)*(4-(hi2-hxj2(i))/(hi2-hr2))-2*log((hi2-hxj2(i))/(hi2-hr2))-3)*(hi2-hr2)^3/b2/E/sf2^3; sigmasj1(i)=2*(1+v)*Fn*cos(miuj1(i))^2*(hr1+(hi1-hr1)*log((hi1-hr1)/(hi1-hxj1(i))))/b1/E/sf1; sigmasj2(i)=2*(1+v)*Fn*cos(miuj2(i))^2*(hr2+(hi2-hr2)*log((hi2-hr2)/(hi2-hxj2(i))))/b2/E/sf2; sigmagj1(i)=24*Fn*hxj1(i)*cos(miuj1(i))^2/pi/b1/E/sf1^2; sigmagj2(i)=24*Fn*hxj2(i)*cos(miuj2(i))^2/pi/b2/E/sf2^2; ki1(i)=Fn/(sigmabri1(i)+sigmabri2(i)+sigmabti1(i)+sigmabti2(i)+sigmasi1(i)+sigmasi2(i)+sigmagi1(i)+sigmagi2(i)+sigmap);%双齿啮合时i点刚度 kj(i)=Fn/(sigmabrj1(i)+sigmabrj2(i)+sigmabtj1(i)+sigmabtj2(i)+sigmasj1(i)+sigmasj2(i)+sigmagj1(i)+sigmagj2(i)+sigmap);%双齿啮合时j点刚度 k(i)=ki1(i)+kj(i);%双齿啮合时，综合时变啮合刚度 %k(i)=ki1(i)*kj(i)/(ki1(i)+kj(i)); kkk(1,i)=k(i); TTT(1,i)=tt(i); end step3=CD/n;%单齿啮合区位移等分，CD为单齿啮合区 for i=1:n xxx(i)=B2C+i*step3;%单齿啮合时啮合线位移变化 ttt(i)=t1+i*step4;%单齿啮合时，时间量变化 N1Bi(i)=N1C+i*step3; O1Bi(i)=sqrt(N1Bi(i)^2+rb1^2); ai1(i)=acos(rb1/O1Bi(i)); gamai1(i)=pi/2/z1+invalpha-tan(ai1(i))+ai1(i); miui1(i)=ai1(i)-gamai1(i); rxi1(i)=O1Bi(i); hxi1(i)=rxi1(i)*cos(gamai1(i))-sqrt(rf1^2-(sf1/2)^2); N2Bi(i)=N2C-i*step3; O2Bi(i)=sqrt(N2Bi(i)^2+rb2^2); ai2(i)=acos(rb2/O2Bi(i)); gamai2(i)=pi/2/z2+invalpha-tan(ai2(i))+ai2(i); miui2(i)=ai2(i)-gamai2(i); rxi2(i)=O2Bi(i); hxi2(i)=rxi2(i)*cos(gamai2(i))-sqrt(rf2^2-(sf2/2)^2); sigmabri1(i)=12*Fn*cos(miui1(i))^2*(hxi1(i)*hr1*(hxi1(i)-hr1)+hxi1(i)^3/3)/b1/E/sf1^3; sigmabri2(i)=12*Fn*cos(miui2(i))^2*(hxi2(i)*hr2*(hxi2(i)-hr2)+hxi2(i)^3/3)/b2/E/sf2^3; sigmabti1(i)=6*Fn*cos(miui1(i))^2*((hi1-hxi1(i))/(hi1-hr1)*(4-(hi1-hxi1(i))/(hi1-hr1))-2*log((hi1-hxi1(i))/(hi1-hr1))-3)*(hi1-hr1)^3/b1/E/sf1^3; sigmabti2(i)=6*Fn*cos(miui2(i))^2*((hi2-hxi2(i))/(hi2-hr2)*(4-(hi2-hxi2(i))/(hi2-hr2))-2*log((hi2-hxi2(i))/(hi2-hr2))-3)*(hi2-hr2)^3/b2/E/sf2^3; sigmasi1(i)=2*(1+v)*Fn*cos(miui1(i))^2*(hr1+(hi1-hr1)*log((hi1-hr1)/(hi1-hxi1(i))))/b1/E/sf1; sigmasi2(i)=2*(1+v)*Fn*cos(miui2(i))^2*(hr2+(hi2-hr2)*log((hi2-hr2)/(hi2-hxi2(i))))/b2/E/sf2; sigmagi1(i)=24*Fn*hxi1(i)*cos(miui1(i))^2/pi/b1/E/sf1^2; sigmagi2(i)=24*Fn*hxi2(i)*cos(miui2(i))^2/pi/b2/E/sf2^2; sigmap=4*Fn*(1-v^2)/pi/b2/E; ki2(i)=Fn/(sigmabri1(i)+sigmabri2(i)+sigmabti1(i)+sigmabti