一、获取代码方式

获取代码方式1：
完整代码已上传我的资源：【力学】基于matlab立铣刀力模拟仿真【含Matlab源码 193期】

获取代码方式2：
通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。

备注：
订阅紫极神光博客付费专栏，可免费获得1份代码（有效期为订阅日起，三天内有效）；

二、部分源代码

% m file to simulate helical end milling forces %
clear;

% Angle conversions
d2r = pi/180; % degree to radian
r2d = 180/pi; % radian to degree

% Workpeice data %
% AL 6061 %
tau_s = 209.24; % shear strength [N/mm2]
phi_c = 27.67*d2r; % shear angle [rad]
B_a = 32.27*d2r; % Friction angle [rad]
kte = 23.66; % tangential edge force coefficent [N/mm]
kfe = 5.57; % normal edge force coefficent [N/mm]

% 4 flute HSS cylindrical endmill %
D = 25.4; % diameter [mm]
r = D/2; % radius [mm]
B = d2r*30.0; % helix angle [rad]
a_n = d2r*5.0; % normal rake angle [rad]
Nf = 2; % number of flutes []

% Machining conditions %
n = 1200; % spindle speed [rpm]
a = 4.0; % axis depth of cut [mm]
f = 1200; % feedrate [mm/min]
calib_factor = 200; % force calibration factor [N/V]

% Assumptions %
i = B; % oblique angle = helix angle [rad]
eta = i; % chip flow angle = oblique angle [rad] [Stabler's chip flow rule]
B_n=atan(tan(B_a)*cos(eta)); %[rad] Friction Coefficient Ba and Shear Stress, Ts, are the same in both orthogonal & Oblique Cutting
phi_n = phi_c; %[rad] Orthogonal Shear Angle = Normal Shear Angle in Oblique
a_r = a_n; %[rad] Orthogonal Normal Rake Angle = Rake Angle in Oblique 
kre = kfe; %[N/mm] Normal Edge Force Coefficient in Orthogonal Cutting is Equivalent to Radial Edge Force in Milling
% Feed Per Tooth
c = f/(Nf*n); % feed per tooth [mm]

% Cutting coefficents
den = sqrt(cos(phi_n+B_n-a_n)^2 + tan(eta)^2*sin(B_n)^2); % Denominator
ktc = (tau_s/sin(phi_c))*(cos(B_n-a_n)+tan(i)*tan(eta)*sin(B_n))/(den);
krc = (tau_s/sin(phi_c))*(cos(B_n-a_n)*tan(i)-tan(eta)*sin(B_n))/(den);
kfc = (tau_s/(sin(phi_c)*cos(i)))*((sin(B_n-a_n))/(den));

% Discritization %
Ts = 0.0001; % Sampling period [sec]
Trev = (n/60); % period of one revolution [sec]
rot_number = 4; % number of revolutions to consider []
dphi = Trev*2*pi*Ts; % angular increment
phi = (0:dphi:rot_number*2*pi); % rotation array for bottom of flute 1 [rad]
Nt = length(phi); % number of time (or rotation) samples to consider
t = (0:Ts:rot_number/Trev); % time array
%t_meas = (0:ts:(N_meas-1)*Ts)' % time array for measured forces [sec]

phi_p = (2*pi)/Nf; % Angle between teeth [rad]

% Vertical Integration %
dz = 0.1; % vertical integration increment [mm]
Nz = a/dz; % number of vertical segments
z = (dz:dz:a); % [mm] vertical integration array

% Immersion %
im = [0.25,0.5,0.75,1]*D; %[mm] Downmilling 1/4,1/2,3/4,1 

% allocate array; columns are each immersion%
Fx = zeros(length(im),Nt); % x axis force history [N]
Fy = zeros(length(im),Nt); % y axis force history [N]
Ft = zeros(length(im),Nt); % t axis force history [N]
F = zeros(length(im),Nt); % resultant force history (in x-y plane) [N] 
T = zeros(length(im),Nt); % torque history [N*m]
P = zeros(length(im),Nt); % Power history

% Main Loop for Simulation %
for count=1:length(im)
    phi_st = pi - (acos((r-im(count))/r)); % [rad] start angle
    phi_ex = pi; %[rad] exit angle
    for kt=1:Nt % Rotation counter
        for kz=1:Nz
            for kf=1:Nf
                phi_bottom = phi(kt) + phi_p*(kz-1); % Immersion angle for tooth
                si = ((2*tan(B)*z(kf))/D); % Correction factor for helx
                phi_cur = phi_bottom-si; % Current phi

                % Ensure current angle between 0 and 2pi
                if phi_cur>=0
                    phi_cur2=phi_cur-2*pi*floor(phi_cur/(2*pi));
                else
                    phi_cur2=phi_cur+2*pi*floor(-phi_cur/(2*pi)+1);
                end

                if (phi_cur2>=phi_st) && (phi_cur2 <=phi_ex)
                    h = c*sin(phi_cur2);

                    dFt = dz * (ktc * h + kte); % tangential force contribution [N]
                    dFr = dz * (krc * h + kre); % radial force contribution [N]
                    dFx = -dFt*cos(phi_cur2) - dFr*sin(phi_cur2); % x axis force contribution [N]
                    dFy = dFt*sin(phi_cur2) - dFr*cos(phi_cur2); % y axis force contribution [N]
                    % (dF dFR dFx dFy)

                    Ft(count,kt) = Ft(count,kt) + dFt; % Intregrate tangential force contribution
                    Fx(count,kt) = Fx(count,kt) + dFx; % Intregrate x axis force contribution
                    Fy(count,kt) = Fy(count,kt) + dFy; % Intregrate y axis force contribution
                end
            end
        end
        F(count,kt) = sqrt(Fx(count,kt)^2 + Fy(count,kt)^2); % resultant Cutting forces history [N]
        T(count,kt) = 1e-3*(D/2)*Ft(count,kt); % cutting torque history [N*m]
        P(count,kt) = T(count,kt)*n*2*pi/60; % Power [W]
    end
end

三、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 门云阁.MATLAB物理计算与可视化[M].清华大学出版社，2013.