TCA齿轮故障分类.zip

syms PhiG ThetaG SGG ETA1 PhiP ThetaP SGP ETA2=0.2*pi/180; %安装误差 J=0;%小轮在2坐标系下沿x方向位移 V=0;%小轮在2坐标系下沿y方向位移 H=0;%小轮在2坐标系下沿z方向位移 Gamma=90*pi/180; %-----------轴交角 %--------大轮convex参数 iG=0; %-------------刀倾角 jG=0; %-------------刀转角 ThetaIG=59.2342*pi/180; %-------------初始摇台角 SRG=(2.85995-0.01)*25.4; %-------------半径设定 DeltaAG=0.0; %-------------机床中心到后部距离 DeltaBG=0.0; %-------------滑动底座 Gammam1G=(67.68333+0.02)*pi/180; %-------------机床安装角 Em1G=0; %-------------机床偏置距 RatioG=1.051674; %-------------滚比 RadiusG=(2.9675-0.01)*25.4; %-------------凸面刀盘半径 AlphaG=(22-0.01)*pi/180; %-------------刀片角 %--------小轮concave参数 iP=0; %-------------刀倾角 jP=0; %-------------刀转角 ThetaIP=63.94203*pi/180; %-------------初始摇台角 SRP=2.947802*25.4; %-------------半径设定 DeltaAP=-0.04023062*25.4; %-------------机床中心到后部距离 DeltaBP=0.01167273*25.4; %-------------滑动底座 Gammam1P=16.86667*pi/180; %-------------机床安装角 Em1P=0.1545759*25.4; %-------------机床偏置距 RatioP=3.242698536; %-------------滚比 RadiusP=2.965621*25.4; %-------------凹面刀盘半径 AlphaP=18.04567*pi/180; %-------------刀片角 %----------大轮齿面方程 QG=ThetaIG-PhiG/RatioG; %----刀具向量---- RgV=[(RadiusG+SGG*sin(AlphaG))*cos(ThetaG); (RadiusG+SGG*sin(AlphaG))*sin(ThetaG); -SGG*cos(AlphaG); 1]; %----刀具法向量---- NgV=[cos(AlphaG)*cos(ThetaG); cos(AlphaG)*sin(ThetaG); sin(AlphaG)]; %----大轮转换矩阵---- MbtG=[cos(iG),0,sin(iG),0; 0,1,0,0; -sin(iG),0,cos(iG),0; 0,0,0,1]; McbG=[-sin(jG),-cos(jG),0,SRG; cos(jG),-sin(jG),0,0; 0,0,1,0; 0,0,0,1]; Mm1cG=[cos(QG),-sin(QG),0,0; sin(QG),cos(QG),0,0; 0,0,1,0; 0,0,0,1]; Mnm1G=[1,0,0,0; 0,1,0,Em1G; 0,0,1,-DeltaBG; 0,0,0,1]; MqnG=[cos(Gammam1G),0,sin(Gammam1G),-DeltaAG; 0,1,0,0; -sin(Gammam1G),0,cos(Gammam1G),0; 0,0,0,1]; M1qG=[1,0,0,0; 0,cos(PhiG),-sin(PhiG),0; 0,sin(PhiG),cos(PhiG),0; 0,0,0,1]; %----刀具法向量使用的转换矩阵---- LbtG=MbtG(1:3,1:3); LcbG=McbG(1:3,1:3); Lm1cG=Mm1cG(1:3,1:3); Lnm1G=Mnm1G(1:3,1:3); LqnG=MqnG(1:3,1:3); L1qG=M1qG(1:3,1:3); R1=M1qG*MqnG*Mnm1G*Mm1cG*McbG*MbtG*RgV; N1=L1qG*LqnG*Lnm1G*Lm1cG*LcbG*LbtG*NgV; %-------大轮啮合方程------ Rm1V=Mm1cG*McbG*MbtG*RgV; %----由于转换后向量元素是四个，无法进行计算，把最后一个无关的元素去掉-- Rm1V(4,:)=[]; Nm1V=Lm1cG*LcbG*LbtG*NgV; Omega1m1G=[cos(Gammam1G); 0; sin(Gammam1G)]; Omegacm1G=[0; 0; 1/RatioG]; Om1AG=[0; -Em1G; DeltaBG]; Vrm1V=cross((Omegacm1G-Omega1m1G),Rm1V)-cross(Om1AG,Omega1m1G); %----啮合方程---- % fG=dot(Nm1V,Vrm1V); %----------小轮齿面方程 QP=ThetaIP+PhiP/RatioP; %-------------摇台角 % %----凹面刀刃上的点在其固连坐标系下的向量---- RpC=[(RadiusP+SGP*sin(AlphaP))*cos(ThetaP); (RadiusP+SGP*sin(AlphaP))*sin(ThetaP); -SGP*cos(AlphaP); 1]; %----凹面刀具法向量---- NpC=[cos(AlphaP)*cos(ThetaP); cos(AlphaP)*sin(ThetaP); sin(AlphaP)]; %----转换矩阵---- MbtP=[cos(iP),0,sin(iP),0; 0,1,0,0; -sin(iP),0,cos(iP),0; 0,0,0,1]; McbP=[-sin(jP),-cos(jP),0,SRP; cos(jP),-sin(jP),0,0; 0,0,1,0; 0,0,0,1]; Mm1cP=[cos(QP),sin(QP),0,0; -sin(QP),cos(QP),0,0; 0,0,1,0; 0,0,0,1]; Mnm1P=[1,0,0,0; 0,1,0,Em1P; 0,0,1,-DeltaBP; 0,0,0,1]; MqnP=[cos(Gammam1P),0,sin(Gammam1P),-DeltaAP; 0,1,0,0; -sin(Gammam1P),0,cos(Gammam1P),0; 0,0,0,1]; M1qP=[1,0,0,0; 0,cos(PhiP),-sin(PhiP),0; 0,sin(PhiP),cos(PhiP),0; 0,0,0,1]; %----刀具法向量使用的转换矩阵---- LbtP=MbtP(1:3,1:3); LcbP=McbP(1:3,1:3); Lm1cP=Mm1cP(1:3,1:3); Lnm1P=Mnm1P(1:3,1:3); LqnP=MqnP(1:3,1:3); L1qP=M1qP(1:3,1:3); R2=M1qP*MqnP*Mnm1P*Mm1cP*McbP*MbtP*RpC; N2=L1qP*LqnP*Lnm1P*Lm1cP*LcbP*LbtP*NpC; %------小轮啮合方程----- %----由于啮合方程在m1坐标系下求解，所以把刀具向量转换到m1坐标系--- Rm1C=Mm1cP*McbP*MbtP*RpC; %----由于转换后向量元素是四个，无法进行计算，把最后一个无关的元素去掉-- Rm1C(4,:)=[]; %----把刀具法向量转换到m1坐标系下--- Nm1C=Lm1cP*LcbP*LbtP*NpC; %----齿坯在m1坐标系下转动角速度---- Omega1m1P=[cos(Gammam1P); 0; sin(Gammam1P)]; %----摇台在m1坐标系下转动角速度---- Omegacm1P=[0; 0; 1/RatioP]; %----在m1坐标系下刀具轴线的偏置---- Om1AP=[0; -Em1P; DeltaBP]; %----摇台和齿坯相对速度---- Vrm1C=cross((Omegacm1P-Omega1m1P),Rm1C)-cross(Om1AP,Omega1m1P); %----啮合方程---- % fP=dot(Nm1C,Vrm1C); %------大小轮之间转换矩阵 Mf11=[1,0,0,0; 0,cos(ETA1),sin(ETA1),0; 0,-sin(ETA1),cos(ETA1),0; 0,0,0,1]; Mf2f1=[cos(Gamma),0,sin(Gamma),0; 0,1,0,0; -sin(Gamma),0,cos(Gamma),0; 0,0,0,1]; M2f2=[1,0,0,-J; 0,cos(pi+ETA2),-sin(pi+ETA2),-V; 0,sin(pi+ETA2),cos(pi+ETA2),-H; 0,0,0,1]; Lf11=Mf11(1:3,1:3); Lf2f1=Mf2f1(1:3,1:3); L2f2=M2f2(1:3,1:3); %都转到f1坐标系下 Rf1G=Mf11*R1; Rf1P=Mf2f1^(-1)*M2f2^(-1)*R2; %大小轮齿面法向量 NG=Lf11*N1; NP=Lf2f1^(-1)*L2f2^(-1)*N2; %大小轮齿面向量 RG=Rf1G(1:3,:); RP=Rf1P(1:3,:); eq1=dot(Nm1V,Vrm1V); eq2=dot(Nm1C,Vrm1C); eq3=RG(1,1)-RP(1,1); eq4=RG(2,1)-RP(2,1); eq5=RG(3,1)-RP(3,1); eq6=NG(1,1)+NP(1,1); eq7=NG(2,1)+NP(2,1); F0=[eq1;eq2;eq3;eq4;eq5;eq6;eq7]; % syms PhiG ThetaG SGG ETA1 PhiP ThetaP SGP A0=[ diff(eq1,PhiG),diff(eq1,ThetaG),diff(eq1,SGG),diff(eq1,ETA1),diff(eq1,PhiP),diff(eq1,ThetaP),diff(eq1,SGP); diff(eq2,PhiG),diff(eq2,ThetaG),diff(eq2,SGG),diff(eq2,ETA1),diff(eq2,PhiP),diff(eq2,ThetaP),diff(eq2,SGP); diff(eq3,PhiG),diff(eq3,ThetaG),diff(eq3,SGG),diff(eq3,ETA1),diff(eq3,PhiP),diff(eq3,ThetaP),diff(eq3,SGP); diff(eq4,PhiG),diff(eq4,ThetaG),diff(eq4,SGG),diff(eq4,ETA1),diff(eq4,PhiP),diff(eq4,ThetaP),diff(eq4,SGP); diff(eq5,PhiG),diff(eq5,ThetaG),diff(eq5,SGG),diff(eq5,ETA1),diff(eq5,PhiP),diff(eq5,ThetaP),diff(eq5,SGP); diff(eq6,PhiG),diff(eq6,ThetaG),diff(eq6,SGG),diff(eq6,ETA1),diff(eq6,PhiP),diff(eq6,ThetaP),diff(eq6,SGP); diff(eq7,PhiG),diff(eq7,ThetaG),diff(eq7,SGG),diff(eq7,ETA1),diff(eq7,PhiP),diff(eq7,ThetaP),diff(eq7,SGP); ]; x=[0; 2.7186; 3.81; 0; 0; 0.3689; 3.81 ]; PhiG=x(1,1); ThetaG=x(2,1); SGG=x(3,1); ETA1=x(4,1); PhiP=x(5,1); ThetaP=x(6,1); SGP=x(7,1); F0=eval(F0); A0=eval(A0); H=inv(A0); x0=x; F=F0; i=1; while 1 x=x-H*F; deltax=H*F; i=i+1; F=[eq1;eq2;eq3;eq4;eq5;eq6;eq7]; PhiG=x(1,1); ThetaG=x(2,1); SGG=x(3,1); ETA1=x(4,1); PhiP=x(5,1); ThetaP=x(6,1); SGP=x(7,1); F=eval(F); yy=F-F0; F0=F; s=x-x0 x0=x; H=H+(s-H*yy)*(s-H*yy)'/((s-H*yy)'*yy); i=i+1; if(abs(deltax(1,1))<0.5e-5) break end end % for zz=1:1:100 % %----------------- % PhiG=x(1,zz); % ThetaG=x(2,zz); % SGG=x(3,zz); % ThetaP=x(4,zz); % SGP=x(5,zz); % %----------- % F(:,zz)=eval(F); % A=[A,eval(A)]; % Z=A(:,zz); % H=Z^(-1); % % x(:,zz+1)=x(:,zz)-H*F(:,zz); % PhiG=x(1,zz+1); % ThetaG=x(2,zz+1); % SGG=x(3,zz+1); % ThetaP=x(4,zz+1); % SGP=x(5,zz+1); % F(:,zz+1)=eval(F); % A(:,zz+1)=(F(:,zz+1)-F(:,zz))/(x(:,zz+1)-x(:,zz)) % error=abs(x(1,zz+1)-x(1,zz)) % % if error<0.001 % break % else % x=x; % end % end % for zz=0:1:10 % x=qi-A^(-1)*F; % if (x-qi)<0.00000001 % break % else % qi=x; % end % end