stability_lobe.zip_stability_lobe_动力学建模_铣_铣削动力学_铣削稳定性

clear % parameters N = 2; % number of teeth Kt = 6e8; % tangential cutting force coef?cient(N/m2) Kn = 2e8; % normal cutting force coef?cient(N/m2) w0 = 922*2*pi; % angular natural frequency (rad/s) zeta = 0.011; % relative damping (1) m_t = 0.03993; % mass (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; % starting depth of cut (m) w_fi = 10e-3; % ?nal depth of cut (m) o_st = 5e3; % starting spindle speed (rpm) o_fi = 25e3; % ?nal spindle speed (rpm) % computational parameters k = 40; % number of discretization interval over one period intk = 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(m + 2, m + 2); % matrix D d = ones(m + 1, 1); d(1 : 2) = 0; D = D + diag(d, -1); D(3,1) = 1; % numerical integration of speci?c cutting force coef?cient according to Equation (37) for i = 1 : k dtr = 2*pi/N/k; %t, if = 2 /N h_i(i) = 0; for j = 1 : N % loop for tooth j for h = 1 : intk % loop for numerical integration of h fi(h) = i*dtr + ( j - 1)*2*pi/N + h*dtr/intk; if(fi(h) >= fist)*(fi(h) <= fiex) g(h) = 1; % tooth is in the cut else g(h) = 0; % tooth is out of cut end end h_i(i) = h_i(i) + sum(g.*(Kt.*cos(fi) + Kn.*sin(fi)).*sin(fi))/intk; 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; % depth of cut % construct transition matrix Fi Fi = eye(m + 2, m + 2); for i = 1 : m A = zeros(2, 2); % matrix Ai A(1,2) = 1; A(2,1) = -w0^2 - h_i(i)*w/m_t; A(2,2) = -2*zeta*w0; B = zeros(2, 2); % matrix Bi B(2,1) = h_i(i)*w/m_t; P = expm(A*dt); % matrix Pi R = (expm(A*dt) - eye(2))*inv(A)*B; % matrix Ri D(1 : 2, 1 : 2) = P; D(1 : 2, m + 1) = wa*R(1 : 2, 1 : 1); D(1 : 2, m + 2) = wb*R(1 : 2, 1 : 1); 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')