clear all
close all
clc
%% ==================== 公共参数定义 ====================
% 结构动力学参数
kx = 2.1425e6; % x方向刚度 (N/m)
wnx = 4035*2*pi; % x方向固有角频率 (rad/s)
zetax = 0.01; % x方向阻尼比
ky = 2.1425e6; % y方向刚度 (N/m)
wny = 4035*2*pi; % y方向固有角频率 (rad/s)
zetay = 0.01; % y方向阻尼比
% 计算质量 (添加缺失的质量参数)
m_t = kx / wnx^2; % 质量 (kg)
% 切削参数
Kt = 1.482e9; % 切向切削力系数 (N/m²)
Kn = 1550/2550; % 法向力系数比例
Nt = 4; % 刀具齿数
aD = 1; % 径向切深比 (100%)
up_or_down = -1; % 铣削方式: -1=顺铣, 1=逆铣
% 计算范围设置
stx = 200; % 主轴转速分点数 (半离散法)
sty = 100; % 轴向切深分点数 (半离散法)
w_st = 0e-3; % 起始轴向切深(m)
w_fi = 90e-3; % 终止轴向切深(m) → 90mm
o_st = 10e3; % 起始主轴转速(rpm)
o_fi = 90e3; % 终止主轴转速(rpm) → 90,000rpm
% 根据铣削类型计算切入切出角
if up_or_down == 1 % 逆铣
fist = 0; % 切入角
fiex = acos(1 - 2*aD); % 切出角
elseif up_or_down == -1 % 顺铣
fist = acos(2*aD - 1); % 切入角
fiex = pi; % 切出角
end
%% ==================== 零阶近似法(ZOA)计算 ====================
fprintf('开始零阶近似法(ZOA)计算...\n');
tic
% 方向系数计算 (修正公式)
alphaxx = (1/(2*pi)) * ((cos(2*fiex) - 2*Kn*fiex + Kn*sin(2*fiex)) - ...
(cos(2*fist) - 2*Kn*fist + Kn*sin(2*fist)));
alphaxy = (1/(2*pi)) * ((-sin(2*fiex) - 2*fiex + Kn*cos(2*fiex)) - ...
(-sin(2*fist) - 2*fist + Kn*cos(2*fist)));
alphayx = (1/(2*pi)) * ((-sin(2*fiex) + 2*fiex + Kn*cos(2*fiex)) - ...
(-sin(2*fist) + 2*fist + Kn*cos(2*fist)));
alphayy = (1/(2*pi)) * ((-cos(2*fiex) - 2*Kn*fiex - Kn*sin(2*fiex)) - ...
(-cos(2*fist) - 2*Kn*fist - Kn*sin(2*fist)));
% 频率范围设置
wnmax = max([wnx wny]);
w = (0:10:2*wnmax)' * pi; % 更密集的频率向量(rad/s)
% 频响函数计算
FRFxx = (1/kx) ./ (1 - (w/wnx).^2 + 1i*2*zetax*(w/wnx)); % m/N (标准形式)
FRFyy = (1/ky) ./ (1 - (w/wny).^2 + 1i*2*zetay*(w/wny)); % m/N
% 特征值计算
lambda = zeros(length(w), 2);
for cnt = 1:length(w)
% 定向频响矩阵
FRF_or = [alphaxx*FRFxx(cnt), alphaxy*FRFyy(cnt);
alphayx*FRFxx(cnt), alphayy*FRFyy(cnt)];
% 计算特征值
E = eig(FRF_or);
lambda(cnt, :) = E;
end
% 稳定性极限计算 (修正公式)
blim = zeros(length(w), 2);
for i = 1:2
ReL = real(lambda(:, i));
ImL = imag(lambda(:, i));
blim(:, i) = -2*pi./(Nt*Kt*ReL); % 正确的稳定性极限公式
end
% 筛选正切深值
valid_idx = blim(:, 1) > 0 | blim(:, 2) > 0;
blim = blim(valid_idx, :);
w_valid = w(valid_idx);
% 提取最小临界切深
blim_min = min(blim, [], 2);
% 相位计算
psi = atan2(imag(lambda(valid_idx, 1)), real(lambda(valid_idx, 1)));
epsilon = pi - 2*psi; % 相位差
% 叶瓣生成
lobes = 6; % 生成6个叶瓣
Omega_zoa = zeros(lobes, length(w_valid));
for cnt = 1:lobes
Omega_zoa(cnt,:) = w_valid/(cnt*2*pi + epsilon)' * (60/(2*pi*Nt));
end
% 合并并筛选最小切深
all_points = [];
for lobe = 1:lobes
[sorted_omega, idx] = sort(Omega_zoa(lobe, :));
valid_idx = find(~isnan(sorted_omega) & sorted_omega >= o_st & sorted_omega <= o_fi);
all_points = [all_points; sorted_omega(valid_idx)', blim_min(idx(valid_idx))];
end
% 按转速排序并提取最小值
[~, idx] = sort(all_points(:, 1));
C_zoa = all_points(idx, :);
min_blim = [];
current_omega = unique(round(C_zoa(:, 1)));
for omega = current_omega'
idx = find(round(C_zoa(:, 1)) == omega);
min_blim = [min_blim; omega, min(C_zoa(idx, 2))];
end
t_zoa = toc;
fprintf('零阶近似法计算完成, 耗时 %.2f 秒\n', t_zoa);
%% ==================== 半离散化法计算 ====================
fprintf('开始半离散化法计算...\n');
tic
% 数值计算参数
k = 40; % 一个周期内的离散点数
m = k; % 时滞周期离散点数
% 方向系数计算
dtr = 2*pi/Nt/k; % 角度步长
h_i = zeros(k, 1);
for i = 1:k
angle = (i-1)*dtr; % 当前角度
for h = 1:k % 积分简化
fi_val = angle + (h-1)*dtr;
% 判断刀齿是否参与切削
if fi_val >= fist && fi_val <= fiex
g = 1;
else
g = 0;
end
% 切削力系数 (简化计算)
h_i(i) = h_i(i) + g * Kt * sin(fi_val) * (cos(fi_val) + Kn*sin(fi_val));
end
h_i(i) = h_i(i) / k; % 平均
end
% 初始化结果矩阵
ss = zeros(stx+1, sty+1);
dc = zeros(stx+1, sty+1);
ei = zeros(stx+1, sty+1);
% 主轴转速循环
for x = 1:stx+1
o = o_st + (x-1)*(o_fi - o_st)/stx; % 当前转速
tau = 60/o/Nt; % 时滞(一个齿的通过时间)
dt = tau/m; % 时间步长
% 轴向切深循环
for y = 1:sty+1
w_val = w_st + (y-1)*(w_fi - w_st)/sty; % 当前切深
% 构建状态转移矩阵Φ
Fi = eye(2); % 初始为单位矩阵
for i = 1:m
% 构建系统矩阵A
A = [0, 1;
-wnx^2 - h_i(i)*w_val/m_t, -2*zetax*wnx];
% 矩阵指数计算 (使用Pade近似提高精度)
A_dt = A*dt;
Fi_step = expm(A_dt); % 使用矩阵指数
% 累积状态转移矩阵
Fi = Fi_step * Fi;
end
% 存储结果
ss(x,y) = o;
dc(x,y) = w_val;
ei(x,y) = max(abs(eig(Fi))); % 最大特征值模
end
% 显示进度
if mod(x, 20) == 0
fprintf('半离散化法进度: %.1f%%\n', x/(stx+1)*100);
end
end
% 提取稳定性边界
contour_matrix = contourc(ss(:,1), dc(1,:), ei', [1, 1]);
if ~isempty(contour_matrix)
start_idx = 1;
boundaries = {};
while start_idx < size(contour_matrix, 2)
level = contour_matrix(1, start_idx);
n_points = contour_matrix(2, start_idx);
x_vals = contour_matrix(1, start_idx+1:start_idx+n_points);
y_vals = contour_matrix(2, start_idx+1:start_idx+n_points);
boundaries{end+1} = [x_vals; y_vals];
start_idx = start_idx + n_points + 1;
end
% 合并边界
x_semi = [];
y_semi = [];
for i = 1:length(boundaries)
x_semi = [x_semi, boundaries{i}(1, :)];
y_semi = [y_semi, boundaries{i}(2, :)];
end
y_semi = y_semi * 1e3; % 转换为mm
else
x_semi = [];
y_semi = [];
end
% 过滤无效点
valid_idx = x_semi >= o_st & x_semi <= o_fi & y_semi >= w_st*1e3 & y_semi <= w_fi*1e3;
x_semi = x_semi(valid_idx);
y_semi = y_semi(valid_idx);
% 排序
[x_semi, sort_idx] = sort(x_semi);
y_semi = y_semi(sort_idx);
t_semi = toc;
fprintf('半离散化法计算完成, 耗时 %.2f 秒\n', t_semi);
%% ==================== 结果可视化 ====================
figure('Position', [100, 100, 900, 600], 'Color', 'w');
% 绘制ZOA结果
plot(min_blim(:, 1), min_blim(:, 2)/1e3, 'r-', 'LineWidth', 2.5, 'DisplayName', '零阶近似法(ZOA)');
hold on;
% 绘制半离散化法结果
if ~isempty(x_semi)
plot(x_semi, y_semi, 'b--', 'LineWidth', 2.5, 'DisplayName', '半离散化法');
end
% 设置图形属性
title('铣削稳定性叶瓣图对比', 'FontSize', 16, 'FontWeight', 'bold');
xlabel('主轴转速 (rpm)', 'FontSize', 14);
ylabel('极限轴向切深 (mm)', 'FontSize', 14);
legend('Location', 'northeast', 'FontSize', 12);
grid on;
box on;
axis([o_st, o_fi, 0, w_fi*1e3]);
set(gca, 'FontSize', 12, 'LineWidth', 1.5);
% 添加计算时间标注
text(0.02, 0.95, sprintf('ZOA计算时间: %.2f秒', t_zoa), 'Units', 'normalized', ...
'FontSize', 10, 'BackgroundColor', [1,1,1,0.7]);
text(0.02, 0.90, sprintf('半离散化法计算时间: %.2f秒', t_semi), 'Units', 'normalized', ...
'FontSize', 10, 'BackgroundColor', [1,1,1,0.7]);
% 添加参数表
param_str = {
sprintf('刚度 (kx, ky): %.2e N/m', kx),
sprintf('固有频率 (fnx, fny): %.0f Hz', wnx/(2*pi)),
sprintf('阻尼比 (ζx, ζy): %.3f', zetax),
sprintf('质量: %.2f kg', m_t),
sprintf('切向力系数 (Kt): %.2e N/m²', Kt),
sprintf('齿数 (Nt): %d', Nt),
sprintf('径向切深比 (aD): %.1f', aD),
sprintf('铣削方式: %s', char(up_or_down==1)*'逆铣' + char(up_or_down==-1)*'顺铣'))
};
annotation('textbox', [0.15, 0.15, 0.25, 0.2], 'String', param_str, ...
'FitBoxToText', 'on', 'BackgroundColor', [1,1,1,0.7], ...
'FontSize', 9, 'EdgeColor', 'k');
%% ==================== 计算性能分析 ====================
fprintf('\n============ 计算性能分析 ============\n');
fprintf('零阶近似法(ZOA):\n');
fprintf(' 计算点数: %d\n', length(w));
fprintf(' 计算时间: %.2f 秒\n', t_zoa);
fprintf(' 速度: %.0f 点/秒\n', length(w)/t_zoa);
fprintf('\n半离散化法:\n');
fprintf(' 主轴转速点数: %d\n', stx);
fprintf(' 轴向切深点数: %d\n', sty);
fprintf(' 总计算点: %d\n', (stx+1)*(sty+1));
fprintf(' 计算时间: %.2f 秒\n', t_semi);
fprintf(' 速度: %.1f 点/秒\n', ((stx+1)*(sty+1))/t_semi);
fprintf('\n半离散化法比ZOA慢 %.1f 倍\n', t_semi/t_zoa);
最新发布
本文提供了一套针对Internet Explorer 6浏览器的兼容性样式解决方案，包括高度、定位及背景抖动等问题的修复方法。通过使用特定的CSS表达式，确保网页在IE6下也能正常显示。
扫一扫
730
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
