MATLABstability_lobe_颤振稳定性分析_铣_动力学_铣削_叶瓣图_源码

clear all close all clc tic %计时开始，检验程序性能 % 定义x方向参数 kx = 2.1425e6; % N/mx方向刚度 wnx = 4035*2*pi; % rad/s,w=2*pi*f 转速大 频率大 阻尼比（模态参数） zetax = 0.01; %x方向阻尼比 % 定义y方向参数 ky = 2.1425e6; % N/m y方向刚度 wny =4035*2*pi; % rad/s y方向阻尼比 zetay = 0.01; %y方向阻尼比 phis = 0*pi/180; % rad 切入角 （逆铣） phie = 180*pi/180; % rad 切出角 （逆铣） % 切削力系数 % Ks = 750e6; % N/m^2 % beta = 68; % deg 螺旋角（以圆柱螺旋立铣刀为背景的） Kn = 1550/2550; Kt = 1.482e9; % N/m^2 切向切削力系数 kr=0.889e9; %径向切削力系数 Nt = 2; %齿数 alphaxx = 0.5*((cos(2*phie)-2*Kn*phie+Kn*sin(2*phie))-(cos(2*phis)-2*Kn*phis+Kn*sin(2*phis))); alphaxy = 0.5*((-sin(2*phie)-2*phie+Kn*cos(2*phie))-(-sin(2*phis)-2*phis+Kn*cos(2*phis))); alphayx = 0.5*((-sin(2*phie)+2*phie+Kn*cos(2*phie))-(-sin(2*phis)+2*phis+Kn*cos(2*phis))); alphayy = 0.5*((-cos(2*phie)-2*Kn*phie-Kn*sin(2*phie))-(-cos(2*phis)-2*Kn*phis-Kn*sin(2*phis))); wnmax = max([wnx wny]); w = (0:1:2*wnmax/2/pi)'*2*pi; % frequency, rad/s前面不是换成弧度了吗？那不是弧度，那是角频率 FRFxx = (wnx^2/kx)./(wnx^2 - w.^2 + i*2*zetax*wnx.*w); % m/N FRFyy = (wny^2/ky)./(wny^2 - w.^2 + i*2*zetay*wny.*w); for cnt = 1:length(w)%取所有的频率 % Oriented FRF FRF_or = [alphaxx*FRFxx(cnt) alphaxy*FRFyy(cnt); alphayx*FRFxx(cnt) alphayy*FRFyy(cnt)]; % m/N这里就是求频响函数特征多项式的特征根的矩阵运算方程就是altitas中的定向传函矩阵，是一系列两行两列的矩阵 % Calculate two eigenvalues E = eig(FRF_or);%求特征方程对应矩阵的特征根 temp = E(1);%只是将特征根的一个根暂时存储在这里 lambda1(cnt) = temp;%求出这个特征根对应的频 temp = E(2);% lambda2(cnt) = temp; if (cnt > 1) dot_prod1 = real(lambda2(cnt))*real(lambda2(cnt-1)) + imag(lambda2(cnt))*imag(lambda2(cnt-1)); dot_prod2 = real(lambda2(cnt))*real(lambda1(cnt-1)) + imag(lambda2(cnt))*imag(lambda1(cnt-1)); if (dot_prod2 > dot_prod1)% 判断特征根2小于特征根1 temp = lambda2(cnt); lambda2(cnt) = lambda1(cnt);%将大值赋给特征根2 lambda1(cnt) = temp;%将小值赋给特征根1 end end end lambda1 = lambda1'; blim1 = (2*pi/Nt/Kt)./((real(lambda1)).^2 + (imag(lambda1)).^2) .* (real(lambda1) .* (1 + (imag(lambda1)./real(lambda1)).^2)); % m分别求出小特征值和大特征值所对应的极限切深论文中公式23 [index1] = find(blim1 > 0);%找出极限切深大于零的值存储的是blim1 所对应的位置号就是频率 blim1 = blim1(index1); blim1 = blim1*1e6; % um w1 = w(index1);%找出大于零的极限切深所对应的频率， psi1 = atan2(imag(lambda1), real(lambda1));%求正切对应的弧度值用它来计算极限切深所对应的转速 psi1 = psi1(index1);%极限切深大于零求正切所对应的弧度值 epsilon1 = pi - 2*psi1;%切入切出角所对应的相位差 % omega11 = (60/Nt)*w1./(epsilon1 + 2*0*pi); % rpm这里的12345代表公式中的K即叶瓣数 % omega12 = (60/Nt)*w1./(epsilon1 + 2*1*pi); % omega13 = (60/Nt)*w1./(epsilon1 + 2*2*pi); % omega14 = (60/Nt)*w1./(epsilon1 + 2*3*pi); % omega15 = (60/Nt)*w1./(epsilon1 + 2*4*pi); % omega16 = (60/Nt)*w1./(epsilon1 + 2*5*pi); lobes=6; w=w1; Omega=zeros(lobes,length(w)); %计算转速，一个epsilon 对应多个转速 for cnt=1:lobes Omega(cnt,:)=w/Nt./((cnt-1)*2*pi + epsilon1)*60; %omega（cnt,:）是取出第cnt行的数据 end %% %循环去交线 x=cell(lobes,1);%cell为取最小整数 y=cell(lobes,1); z=cell(lobes,1); % 对每条叶瓣进行步距为1的插值操作 CZ=max(Omega); for cnt=1:lobes %x{cnt}=ceil(Omega(cnt,1)):1: floor(Omega(cnt,length(Omega(cnt,:))));%floor为圆整 x{cnt}= 1:10:CZ; x{cnt}=x{cnt}';%x{cnt}控制的是插值的横坐标， 横坐标为转速 y{cnt}=interp1(Omega(cnt,:),blim1,x{cnt}); %interp1() 一维插值函数调用转速与极限切深的插值y{cnt}是极限切深 z{cnt}=interp1(Omega(cnt,:),w,x{cnt}); %颤振频率的插值，转速与颤振频率的插值z{cnt}是转速 end %去相交曲线准备，B保存不重复转速 A=[]; for cnt=1:lobes A=cat(1,A,x{cnt});%cat就是将两个矩阵链接起来 A和x{cnt}按列连接起来 end B=unique(A);%就是将矩阵中重复的元素剔除这里是算出不重复的转速 C=zeros(length(B),3); C(:,1)=B; %取出第一列数据 C(:,2)=1/eps; %机器的浮点运算误差限 C(:,3)=1/eps; %eps为最小值，取倒数为最大值 %从左至右，判断切宽，取小值作为最终值，C(:,1)转速，C(：，2）切宽，C（：，3）颤振频率 for cnt=lobes:-1:1 for i=1:length(x{cnt}) N=find(C(:,1)==x{cnt}(i,1));%从每个叶瓣中寻找对应的转速 if y{cnt}(i,1)<C(N,2) %比较切深 C(N,2)=y{cnt}(i,1);%就让转速对应的切深赋值 C(N,3)=z{cnt}(i,1);%Z为转速和频率表的插值 end end end % 绘制SLD图 %[AX,H1,H2]=plotyy(C(:,1),C(:,2),C(:,1),C(:,3)); plot(C(:,1),C(:,2),'r-') hold on axis([10000 90000 0 90])%设置横坐标主轴范围 %axis([500 6000 0 1])%设置横坐标主轴范围 set(gca,'FontSize', 10)%设置标注显示字体大小 xlabel('主轴转速(rpm)') ylabel('纵向切削深度 (mm)') toc% 计时结束