the 2 dof milling analysis.m

在加工过程中

clear 
% parameters
N=2;        % number of teeth
Kt=6e8;     % tangential cutting force coefficient(N/m^2)
Kn=2e8;     % normal cutting force coefficient(N/m^2)
w0x=922*2*pi; %angular natural frequency x(rad/s)
zetax=0.011;  % relative damping x(1)
m_tx=0.03993; % mass x(kg)
w0y=922*2*pi; %angular natural frequency y(rad/s)
zetay=0.011;  % relative damping y(1)
m_ty=0.03993; % mass y(kg)
aD=0.05;      % radial depth of cut
up_or_down=-1; % 1: up-milling, -1: down-milling
if up_or_down==1  % up-milling
    fist=0;       % start angle
    fiex=acos(1-2*aD); % exit angle
elseif up_or_down==-1  % down-milling
    fist=acos(2*aD-1);           % start angle
    fiex= pi;               % exit angle
end 
stx=400;           % steps of spindle speed
sty=200;           % steps of depth of cut
w_st=0e-3;         % strating depth of cut(m)
w_fi=10e-3;        % final depth of cut (m)
o_st=5000;          % starting spindle speed (rpm)
o_fi=25000;         % final spindle speed (rpm)
%% computational parameters     
k=40;              % number of discretization interval over one period
fix k=20;          % number of numerical integration steps for equation(37)
m=k;               % since time delay=time period
wa=1/2;            % since time delay=time period 
wb=1/2;            % since time delay=time period
D=zeros(2*m+4,2*m+4); % matrix D
d=ones(2*m+2,1);
d(1:4)=0;
D=D+diag(d,-2);
D(5,1)=1;
D(6,2)=1;
%% numerical integration of specific cutting force coefficient according to
%% Equations (40)-(43)
for  i=1:k    % loop for tooth j
    dtr=2*pi/N/k;          
    hxx(i)=0;
    hxy(i)=0;
    hyx(i)=0;
    hyy(i)=0;
    for j=1:N        % loop for tooth j
        for h=1:k   % loop for numerical integration of hi
            fi(h)=i*dtr+(j-1)*2*pi/N+h*dtr/k;
            if (fi(h)>=fist)*(fi(h)<=fiex)
                g(h)=1;
            else
                g(h)=0;
            end
        end
        hxx(i)=hxx(i)+sum(g.*(Kt.*cos(fi)+Kn.*sin(fi)).*sin(fi))/k;
        hxy(i)=hxy(i)+sum(g.*(Kt.*cos(fi)+Kn.*sin(fi)).*cos(fi))/k;
        hyx(i)=hyx(i)+sum(g.*(-Kt.*sin(fi)+Kn.*cos(fi)).*sin(fi))/k;
        hyy(i)=hyy(i)+sum(g.*(-Kt.*sin(fi)+Kn.*cos(fi)).*cos(fi))/k;
    end 
end
%% start of computation
for x=1:stx+1         % loop for spindle speeds
    o=o_st+(x-1)*(o_fi-o_st)/stx;   % spindle speed 
    tau=60/o/N;                     % time delay
    dt=tau/(m);                     % time step
    for y=1:sty+1                   % loop for depth of cuts
    w=w_st+(y-1)*(w_fi-w_st)/sty;
%% construct transition matrix Fi
    Fi=eye(2*m+4,2*m+4);
    for i=1:m
        A=zeros(4,4);
        A(1,3)=1;
        A(2,4)=1;
        A(3,1)=-w0x^2-hxx(i)*w/m_tx;
        A(3,2)=-hxy(i)*w/m_tx;
        A(3,3)=-2*zetax*w0x;
        A(4,1)=-hyx(i)*w/m_ty;
        A(4,2)=-w0y^2-hyy(i)*w/m_ty;
        A(4,4)=-2*zetay*w0y;
        B=zeros(4,4);
        B(3,1)=hxx(i)*w/m_tx;
        B(3,2)=hxy(i)*w/m_tx;
        B(4,1)=hyx(i)*w/m_ty;
        B(4,2)=hyy(i)*w/m_ty;
        P=expm(A*dt);          % matrix Pi
        R=(expm(A*dt)-eye(4))*inv(A)*B;    %matrix Ri
        D(1:4,1:4)=P;
        D(1:4,(2*m+1):(2*m+2))=wa*R(1:4,1:2);
        D(1:4,(2*m+3):(2*m+4))=wb*R(1:4,1:2);
        Fi=D*Fi;          % transition matrix 
    end 
    ss(x,y)=o;      % matrix of spindle speeds
    dc(x,y)=w;      % matrix of depth of cuts 
    ei(x,y)=max(abs(eig(Fi))); % matrix of eigenvalues
    end
    stx+1-x
end
figure
contour(ss,dc,ei,[1,1],'k')