基于势能原理的圆柱齿轮时变啮合刚度计算程序

计算圆柱齿轮的时变啮合刚度。该程序基于势能原理（包括弯曲能、剪切能和轴向压缩能）和赫兹接触理论。

1. 主程序：齿轮啮合刚度计算

%% 基于势能原理的圆柱齿轮时变啮合刚度计算程序
% 作者: MATLAB助手
% 功能: 计算圆柱齿轮副的时变啮合刚度
% 原理: 基于势能法（弯曲、剪切、轴向压缩能）和赫兹接触理论
% 参考: Tian X., 齿轮系统动力学 - 振动、冲击、噪声
clear; clc; close all;

%% 1. 齿轮参数输入
fprintf('========== 圆柱齿轮啮合刚度计算程序 ==========\n');

% 齿轮1参数（小齿轮）
gear1.modulus = 2;        % 模数 m (mm)
gear1.teeth = 30;         % 齿数 z1
gear1.pressure_angle = 20; % 压力角 α (度)
gear1.face_width = 20;    % 齿宽 b (mm)
gear1.youngs_modulus = 2.06e5; % 弹性模量 E (MPa)
gear1.poisson_ratio = 0.3;    % 泊松比 ν
gear1.addendum_coefficient = 1; % 齿顶高系数 ha*
gear1.clearance_coefficient = 0.25; % 顶隙系数 c*
gear1.shift_coefficient = 0;    % 变位系数 x

% 齿轮2参数（大齿轮）
gear2.modulus = 2;        % 模数 m (mm)
gear2.teeth = 50;         % 齿数 z2
gear2.pressure_angle = 20; % 压力角 α (度)
gear2.face_width = 20;    % 齿宽 b (mm)
gear2.youngs_modulus = 2.06e5; % 弹性模量 E (MPa)
gear2.poisson_ratio = 0.3;    % 泊松比 ν
gear2.addendum_coefficient = 1; % 齿顶高系数 ha*
gear2.clearance_coefficient = 0.25; % 顶隙系数 c*
gear2.shift_coefficient = 0;    % 变位系数 x

% 材料参数（假设两齿轮材料相同）
material.E = gear1.youngs_modulus;  % MPa
material.v = gear1.poisson_ratio;   % 泊松比

fprintf('齿轮参数:\n');
fprintf('  模数 m = %.1f mm\n', gear1.modulus);
fprintf('  齿数: z1 = %d, z2 = %d\n', gear1.teeth, gear2.teeth);
fprintf('  压力角 α = %.1f°\n', gear1.pressure_angle);
fprintf('  齿宽 b = %.1f mm\n', gear1.face_width);
fprintf('  弹性模量 E = %.2e MPa\n', material.E);
fprintf('  泊松比 ν = %.2f\n', material.v);

%% 2. 计算齿轮基本几何参数
fprintf('\n2. 计算齿轮几何参数...\n');

% 转换为弧度
alpha = deg2rad(gear1.pressure_angle); % 压力角(弧度)

% 分度圆直径
d1 = gear1.modulus * gear1.teeth; % 齿轮1分度圆直径(mm)
d2 = gear2.modulus * gear2.teeth; % 齿轮2分度圆直径(mm)

% 基圆直径
db1 = d1 * cos(alpha); % 齿轮1基圆直径(mm)
db2 = d2 * cos(alpha); % 齿轮2基圆直径(mm)

% 齿顶圆直径
da1 = d1 + 2 * gear1.modulus * (gear1.addendum_coefficient + gear1.shift_coefficient); % 齿轮1齿顶圆直径(mm)
da2 = d2 + 2 * gear2.modulus * (gear2.addendum_coefficient + gear2.shift_coefficient); % 齿轮2齿顶圆直径(mm)

% 齿根圆直径
df1 = d1 - 2 * gear1.modulus * (gear1.addendum_coefficient + gear1.clearance_coefficient - gear1.shift_coefficient); % 齿轮1齿根圆直径(mm)
df2 = d2 - 2 * gear2.modulus * (gear2.addendum_coefficient + gear2.clearance_coefficient - gear2.shift_coefficient); % 齿轮2齿根圆直径(mm)

% 中心距
a = (d1 + d2) / 2; % 标准中心距(mm)

% 啮合线长度
g_alpha = sqrt(da1^2 - db1^2) / (2 * cos(alpha)) + sqrt(da2^2 - db2^2) / (2 * cos(alpha)) - a * sin(alpha);
pb = pi * gear1.modulus * cos(alpha); % 基圆齿距

% 重合度
epsilon_alpha = g_alpha / (pi * gear1.modulus * cos(alpha));

fprintf('  分度圆直径: d1 = %.2f mm, d2 = %.2f mm\n', d1, d2);
fprintf('  基圆直径: db1 = %.2f mm, db2 = %.2f mm\n', db1, db2);
fprintf('  齿顶圆直径: da1 = %.2f mm, da2 = %.2f mm\n', da1, da2);
fprintf('  齿根圆直径: df1 = %.2f mm, df2 = %.2f mm\n', df1, df2);
fprintf('  中心距: a = %.2f mm\n', a);
fprintf('  啮合线长度: g_α = %.4f mm\n', g_alpha);
fprintf('  基圆齿距: pb = %.4f mm\n', pb);
fprintf('  端面重合度: ε_α = %.4f\n', epsilon_alpha);

%% 3. 定义计算参数
fprintf('\n3. 设置计算参数...\n');

% 啮合点采样
n_points = 100; % 沿啮合线的采样点数
theta_range = linspace(0, 2*pi, n_points); % 一个啮合周期的角度采样

% 轮齿离散化参数
n_slices = 50; % 每个齿的离散切片数

fprintf('  啮合周期采样点数: %d\n', n_points);
fprintf('  轮齿离散切片数: %d\n', n_slices);

%% 4. 计算单个轮齿的刚度（基于势能原理）
fprintf('\n4. 计算单个轮齿刚度（势能法）...\n');

% 计算齿轮1的单个齿刚度
fprintf('  计算齿轮1单齿刚度...\n');
single_tooth_stiffness_1 = zeros(n_points, 1);
for i = 1:n_points
    % 当前啮合点位置（从齿根到齿顶）
    contact_position = (i-1)/(n_points-1); % 0:齿根, 1:齿顶
    
    % 计算齿轮1在当前接触点的单齿刚度
    single_tooth_stiffness_1(i) = calculate_single_tooth_stiffness(...
        gear1, material, contact_position, n_slices);
end

% 计算齿轮2的单个齿刚度
fprintf('  计算齿轮2单齿刚度...\n');
single_tooth_stiffness_2 = zeros(n_points, 1);
for i = 1:n_points
    % 当前啮合点位置（从齿根到齿顶）
    contact_position = (i-1)/(n_points-1); % 0:齿根, 1:齿顶
    
    % 计算齿轮2在当前接触点的单齿刚度
    single_tooth_stiffness_2(i) = calculate_single_tooth_stiffness(...
        gear2, material, contact_position, n_slices);
end

fprintf('  齿轮1单齿刚度范围: %.2e ~ %.2e N/mm\n', ...
    min(single_tooth_stiffness_1), max(single_tooth_stiffness_1));
fprintf('  齿轮2单齿刚度范围: %.2e ~ %.2e N/mm\n', ...
    min(single_tooth_stiffness_2), max(single_tooth_stiffness_2));

%% 5. 计算赫兹接触刚度
fprintf('\n5. 计算赫兹接触刚度...\n');

% 计算综合曲率半径
rho1 = db1/2 * tan(alpha); % 齿轮1在啮合点的近似曲率半径
rho2 = db2/2 * tan(alpha); % 齿轮2在啮合点的近似曲率半径
rho_eq = (rho1 * rho2) / (rho1 + rho2); % 等效曲率半径

% 赫兹接触刚度公式: k_h = (π * E * b) / [4 * (1 - ν^2)]
hertz_stiffness = (pi * material.E * gear1.face_width) / (4 * (1 - material.v^2));

fprintf('  齿轮1曲率半径: ρ1 = %.4f mm\n', rho1);
fprintf('  齿轮2曲率半径: ρ2 = %.4f mm\n', rho2);
fprintf('  等效曲率半径: ρ_eq = %.4f mm\n', rho_eq);
fprintf('  赫兹接触刚度: k_h = %.4e N/mm\n', hertz_stiffness);

%% 6. 计算单齿对啮合刚度
fprintf('\n6. 计算单齿对啮合刚度...\n');

% 串联刚度：1/k_mesh = 1/k_gear1 + 1/k_gear2 + 1/k_hertz
single_pair_stiffness = zeros(n_points, 1);
for i = 1:n_points
    k1 = single_tooth_stiffness_1(i);
    k2 = single_tooth_stiffness_2(i);
    kh = hertz_stiffness;
    
    % 串联刚度计算
    single_pair_stiffness(i) = 1 / (1/k1 + 1/k2 + 1/kh);
end

fprintf('  单齿对啮合刚度范围: %.2e ~ %.2e N/mm\n', ...
    min(single_pair_stiffness), max(single_pair_stiffness));

%% 7. 计算多齿对啮合刚度（考虑重合度）
fprintf('\n7. 计算时变啮合刚度（考虑重合度）...\n');

% 计算啮合刚度随时间（或转角）的变化
time_mesh_stiffness = zeros(n_points, 1);
contact_ratio = zeros(n_points, 1);

% 计算每个位置的接触齿对数
for i = 1:n_points
    % 当前啮合点位置（0到1，表示从啮入到啮出）
    pos_ratio = (i-1)/(n_points-1);
    
    % 根据重合度确定接触齿对数
    if epsilon_alpha < 2
        % 单双齿交替啮合区
        if pos_ratio < (epsilon_alpha - 1)
            % 双齿啮合区
            contact_ratio(i) = 2;
            % 找到相邻齿的刚度（相位差一个基节）
            j = mod(i + round(n_points * pb/g_alpha) - 1, n_points) + 1;
            time_mesh_stiffness(i) = single_pair_stiffness(i) + single_pair_stiffness(j);
        elseif pos_ratio > 1 - (epsilon_alpha - 1)
            % 双齿啮合区
            contact_ratio(i) = 2;
            j = mod(i - round(n_points * pb/g_alpha) - 1, n_points) + 1;
            time_mesh_stiffness(i) = single_pair_stiffness(i) + single_pair_stiffness(j);
        else
            % 单齿啮合区
            contact_ratio(i) = 1;
            time_mesh_stiffness(i) = single_pair_stiffness(i);
        end
    else
        % 重合度大于2的情况（较少见）
        contact_ratio(i) = floor(epsilon_alpha);
        time_mesh_stiffness(i) = contact_ratio(i) * single_pair_stiffness(i);
    end
end

fprintf('  时变啮合刚度范围: %.2e ~ %.2e N/mm\n', ...
    min(time_mesh_stiffness), max(time_mesh_stiffness));
fprintf('  平均啮合刚度: %.2e N/mm\n', mean(time_mesh_stiffness));
fprintf('  刚度波动幅度: %.2f%%\n', ...
    100*(max(time_mesh_stiffness)-min(time_mesh_stiffness))/mean(time_mesh_stiffness));

%% 8. 计算平均啮合刚度（解析法验证）
fprintf('\n8. 计算平均啮合刚度（ISO标准方法）...\n');

% ISO 6336-1:2006方法计算平均啮合刚度
% 刚度修正系数
c_gamma = 0.8; % 对于直齿轮
c_beta = 1;    % 对于直齿轮

% 理论刚度
c_th = c_gamma * c_beta / (0.04723 + 0.15551/gear1.teeth + 0.25791/gear2.teeth - 0.00635*gear1.pressure_angle - 0.11654*gear1.shift_coefficient/gear1.teeth - 0.00193*gear2.shift_coefficient - 0.24188*gear2.shift_coefficient/gear2.teeth + 0.00529*gear1.shift_coefficient^2 + 0.00182*gear2.shift_coefficient^2);

% 考虑齿宽和模数的刚度
k_iso = material.E * gear1.face_width * c_th; % N/mm

fprintf('  ISO标准平均啮合刚度: %.2e N/mm\n', k_iso);
fprintf('  势能法平均啮合刚度: %.2e N/mm\n', mean(time_mesh_stiffness));
fprintf('  两种方法差异: %.2f%%\n', ...
    100*abs(mean(time_mesh_stiffness)-k_iso)/k_iso);

%% 9. 可视化结果
fprintf('\n9. 生成可视化图表...\n');

figure('Position', [100, 100, 1200, 800]);

% 9.1 单齿刚度对比
subplot(2, 3, 1);
x_ratio = linspace(0, 1, n_points);
plot(x_ratio, single_tooth_stiffness_1/1e3, 'b-', 'LineWidth', 2);
hold on;
plot(x_ratio, single_tooth_stiffness_2/1e3, 'r--', 'LineWidth', 2);
xlabel('啮合点位置 (0:齿根, 1:齿顶)');
ylabel('单齿刚度 (N/μm)');
title('单齿刚度变化');
legend(['齿轮1 (z=' num2str(gear1.teeth) ')'], ['齿轮2 (z=' num2str(gear2.teeth) ')'], 'Location', 'best');
grid on;

% 9.2 单齿对啮合刚度
subplot(2, 3, 2);
plot(x_ratio, single_pair_stiffness/1e3, 'g-', 'LineWidth', 2);
xlabel('啮合点位置 (0:齿根, 1:齿顶)');
ylabel('单齿对啮合刚度 (N/μm)');
title('单齿对啮合刚度变化');
grid on;

% 9.3 时变啮合刚度
subplot(2, 3, 3);
yyaxis left;
plot(x_ratio, time_mesh_stiffness/1e3, 'k-', 'LineWidth', 2);
ylabel('总啮合刚度 (N/μm)');
yyaxis right;
plot(x_ratio, contact_ratio, 'm:', 'LineWidth', 1.5);
ylabel('接触齿对数');
xlabel('啮合点位置 (0:齿根, 1:齿顶)');
title('时变啮合刚度与接触齿对数');
legend('总啮合刚度', '接触齿对数', 'Location', 'best');
grid on;

% 9.4 刚度分量对比
subplot(2, 3, 4);
stiffness_components = zeros(n_points, 3);
for i = 1:n_points
    k1 = single_tooth_stiffness_1(i);
    k2 = single_tooth_stiffness_2(i);
    kh = hertz_stiffness;
    k_total = 1/(1/k1 + 1/k2 + 1/kh);
    
    stiffness_components(i, 1) = 1/(1/k1) * (k_total/(1/k1 + 1/k2 + 1/kh));
    stiffness_components(i, 2) = 1/(1/k2) * (k_total/(1/k1 + 1/k2 + 1/kh));
    stiffness_components(i, 3) = 1/(1/kh) * (k_total/(1/k1 + 1/k2 + 1/kh));
end

area(x_ratio, stiffness_components/1e3);
xlabel('啮合点位置 (0:齿根, 1:齿顶)');
ylabel('刚度分量 (N/μm)');
title('啮合刚度分量分解');
legend('齿轮1变形', '齿轮2变形', '赫兹接触', 'Location', 'best');
grid on;

% 9.5 轮齿几何和受力分析
subplot(2, 3, 5);
% 绘制齿廓示意图
theta = linspace(0, 2*pi/gear1.teeth, 100);
r_base = db1/2;
r_tip = da1/2;
r_root = df1/2;

% 渐开线齿廓（简化）
inv_angle = @(phi) tan(phi) - phi; % 渐开线函数
phi_range = linspace(0, acos(r_base/(r_tip)), 50);
x_tooth = r_base * (cos(phi_range) + phi_range .* sin(phi_range));
y_tooth = r_base * (sin(phi_range) - phi_range .* cos(phi_range));

% 旋转到正确位置
rotation_angle = pi/(2*gear1.teeth) - inv_angle(alpha);
R = [cos(rotation_angle), -sin(rotation_angle); 
     sin(rotation_angle), cos(rotation_angle)];
tooth_points = R * [x_tooth; y_tooth];

plot(tooth_points(1,:), tooth_points(2,:), 'b-', 'LineWidth', 2);
hold on;
% 绘制基圆、分度圆、齿顶圆
theta_circle = linspace(0, 2*pi, 100);
plot(r_base*cos(theta_circle), r_base*sin(theta_circle), 'k--');
plot((d1/2)*cos(theta_circle), (d1/2)*sin(theta_circle), 'g--');
plot(r_tip*cos(theta_circle), r_tip*sin(theta_circle), 'r--');

axis equal;
xlabel('x (mm)');
ylabel('y (mm)');
title('轮齿几何示意图');
legend('齿廓', '基圆', '分度圆', '齿顶圆', 'Location', 'best');
grid on;

% 9.6 刚度统计分布
subplot(2, 3, 6);
histogram(time_mesh_stiffness/1e3, 20, 'FaceColor', 'c', 'EdgeColor', 'k');
xlabel('啮合刚度 (N/μm)');
ylabel('频率');
title('啮合刚度分布直方图');
grid on;

% 添加统计信息
avg_stiffness = mean(time_mesh_stiffness/1e3);
std_stiffness = std(time_mesh_stiffness/1e3);
text(0.05, 0.95, sprintf('平均值: %.1f', avg_stiffness), 'Units', 'normalized', 'FontSize', 10);
text(0.05, 0.90, sprintf('标准差: %.1f', std_stiffness), 'Units', 'normalized', 'FontSize', 10);
text(0.05, 0.85, sprintf('最小值: %.1f', min(time_mesh_stiffness/1e3)), 'Units', 'normalized', 'FontSize', 10);
text(0.05, 0.80, sprintf('最大值: %.1f', max(time_mesh_stiffness/1e3)), 'Units', 'normalized', 'FontSize', 10);

sgtitle(['圆柱齿轮啮合刚度分析 (z1=' num2str(gear1.teeth) ', z2=' num2str(gear2.teeth) ...
         ', m=' num2str(gear1.modulus) 'mm, ε_α=' sprintf('%.3f', epsilon_alpha) ')'], ...
         'FontSize', 14, 'FontWeight', 'bold');

%% 10. 保存结果
fprintf('\n10. 保存计算结果...\n');

% 保存刚度数据
stiffness_results.epsilon_alpha = epsilon_alpha;
stiffness_results.single_tooth_stiffness_1 = single_tooth_stiffness_1;
stiffness_results.single_tooth_stiffness_2 = single_tooth_stiffness_2;
stiffness_results.hertz_stiffness = hertz_stiffness;
stiffness_results.single_pair_stiffness = single_pair_stiffness;
stiffness_results.time_mesh_stiffness = time_mesh_stiffness;
stiffness_results.contact_ratio = contact_ratio;
stiffness_results.avg_stiffness = mean(time_mesh_stiffness);
stiffness_results.iso_stiffness = k_iso;

% 保存齿轮参数
gear_params.gear1 = gear1;
gear_params.gear2 = gear2;
gear_params.material = material;
gear_params.geometry.d1 = d1;
gear_params.geometry.d2 = d2;
gear_params.geometry.db1 = db1;
gear_params.geometry.db2 = db2;
gear_params.geometry.da1 = da1;
gear_params.geometry.da2 = da2;
gear_params.geometry.df1 = df1;
gear_params.geometry.df2 = df2;
gear_params.geometry.a = a;
gear_params.geometry.g_alpha = g_alpha;
gear_params.geometry.pb = pb;

save('gear_mesh_stiffness_results.mat', 'stiffness_results', 'gear_params');
fprintf('  结果已保存到 gear_mesh_stiffness_results.mat\n');

% 生成报告文件
generate_report(gear_params, stiffness_results);

fprintf('\n========== 计算完成 ==========\n');

2. 关键函数：单齿刚度计算（势能法）

function k_tooth = calculate_single_tooth_stiffness(gear, material, contact_position, n_slices)
% 基于势能法计算单个轮齿的刚度
% 输入:
%   gear - 齿轮结构体
%   material - 材料属性结构体
%   contact_position - 接触点位置 (0:齿根, 1:齿顶)
%   n_slices - 离散切片数
% 输出:
%   k_tooth - 单齿刚度 (N/mm)

% 齿轮参数
m = gear.modulus;           % 模数 (mm)
z = gear.teeth;             % 齿数
alpha = deg2rad(gear.pressure_angle); % 压力角 (弧度)
b = gear.face_width;        % 齿宽 (mm)
E = material.E;             % 弹性模量 (MPa)
v = material.v;             % 泊松比

% 计算齿轮几何
d = m * z;                  % 分度圆直径 (mm)
db = d * cos(alpha);        % 基圆直径 (mm)
da = d + 2 * m;             % 齿顶圆直径 (简化) (mm)
df = d - 2 * 1.25 * m;      % 齿根圆直径 (简化) (mm)

% 计算接触点到轮齿对称线的距离
% 接触点半径
r_contact = df/2 + contact_position * (da/2 - df/2);

% 计算接触点处的压力角
alpha_contact = acos(db / (2 * r_contact));

% 计算力臂 (接触点到齿根危险截面的距离)
% 齿根危险截面位置 (30°切线法简化)
% 实际应用中可以使用更精确的方法
h_f = 1.25 * m; % 齿根高
r_root = df/2;  % 齿根圆半径

% 离散化轮齿截面
y_positions = linspace(r_root, r_contact, n_slices); % 从齿根到接触点
dy = (r_contact - r_root) / (n_slices - 1); % 切片厚度

% 初始化能量积分
U_bending = 0;   % 弯曲能
U_shear = 0;     % 剪切能
U_axial = 0;     % 轴向压缩能

% 单位载荷
F = 1; % 单位力 (N)

% 对每个切片计算能量
for i = 1:n_slices
    % 当前切片位置
    y = y_positions(i);
    
    % 计算当前截面的几何参数
    % 1. 计算齿厚 (简化方法)
    % 在半径为y处的圆周齿厚
    s_y = calculate_tooth_thickness(y, gear);
    
    % 2. 计算截面面积和惯性矩
    A = b * s_y; % 截面面积 (mm²)
    
    % 惯性矩 (矩形截面)
    I = (b * s_y^3) / 12; % 截面惯性矩 (mm⁴)
    
    % 3. 计算弯矩、剪力和轴力
    % 力臂: 从当前截面到接触点的距离
    moment_arm = r_contact - y;
    
    % 载荷分解:
    % 法向力Fn分解为切向力Ft和径向力Fr
    % 在接触点处，载荷分解为:
    F_t = F * cos(alpha_contact); % 切向力
    F_r = F * sin(alpha_contact); % 径向力
    
    % 当前截面处的弯矩
    M = F_t * moment_arm;
    
    % 当前截面处的剪力 (切向分量)
    V = F_t;
    
    % 当前截面处的轴力 (径向分量)
    N = F_r;
    
    % 4. 计算能量密度并积分
    % 弯曲能密度
    if I > 0
        dU_bending = (M^2) / (2 * E * I) * dy;
        U_bending = U_bending + dU_bending;
    end
    
    % 剪切能密度 (考虑剪切形状系数)
    k_shear = 1.2; % 矩形截面的剪切形状系数
    dU_shear = (k_shear * V^2) / (2 * G(E, v) * A) * dy;
    U_shear = U_shear + dU_shear;
    
    % 轴向压缩能密度
    dU_axial = (N^2) / (2 * E * A) * dy;
    U_axial = U_axial + dU_axial;
end

% 总势能
U_total = U_bending + U_shear + U_axial;

% 单齿刚度 (基于能量法: k = F²/(2U)，其中F=1)
if U_total > 0
    k_tooth = 1 / (2 * U_total);
else
    k_tooth = 1e10; % 极大值，避免除零
end

% 确保刚度为正
k_tooth = max(k_tooth, 1e-6);

end

function s = calculate_tooth_thickness(r, gear)
% 计算在半径r处的齿厚
% 简化方法: 线性插值

m = gear.modulus;
z = gear.teeth;
alpha = deg2rad(gear.pressure_angle);

% 分度圆齿厚
s_pitch = pi * m / 2;

% 分度圆半径
r_pitch = m * z / 2;

% 基圆半径
r_base = r_pitch * cos(alpha);

% 齿顶圆半径
r_tip = r_pitch + m;

% 齿根圆半径
r_root = r_pitch - 1.25 * m;

% 如果半径小于基圆半径，使用基圆处齿厚
if r < r_base
    s = s_pitch * r / r_pitch - 2 * r * (inv(alpha) - inv(acos(r_base/r)));
else
    % 渐开线部分: 使用渐开线函数计算
    alpha_r = acos(r_base / r);
    s = 2 * r * (s_pitch/(2*r_pitch) + inv(alpha) - inv(alpha_r));
end

% 确保齿厚不为负
s = max(s, 0.1 * m);

end

function G = G(E, v)
% 计算剪切模量
G = E / (2 * (1 + v));
end

function inv_phi = inv(phi)
% 渐开线函数 inv(φ) = tan(φ) - φ
inv_phi = tan(phi) - phi;
end

3. 辅助函数：生成报告

function generate_report(gear_params, stiffness_results)
% 生成计算结果报告

fprintf('\n========== 计算结果报告 ==========\n');
fprintf('齿轮参数:\n');
fprintf('  齿轮1: z1=%d, m=%.2fmm, α=%.1f°, b=%.1fmm\n', ...
    gear_params.gear1.teeth, gear_params.gear1.modulus, ...
    gear_params.gear1.pressure_angle, gear_params.gear1.face_width);
fprintf('  齿轮2: z2=%d, m=%.2fmm, α=%.1f°, b=%.1fmm\n', ...
    gear_params.gear2.teeth, gear_params.gear2.modulus, ...
    gear_params.gear2.pressure_angle, gear_params.gear2.face_width);

fprintf('\n几何参数:\n');
fprintf('  分度圆直径: d1=%.2fmm, d2=%.2fmm\n', ...
    gear_params.geometry.d1, gear_params.geometry.d2);
fprintf('  基圆直径: db1=%.2fmm, db2=%.2fmm\n', ...
    gear_params.geometry.db1, gear_params.geometry.db2);
fprintf('  齿顶圆直径: da1=%.2fmm, da2=%.2fmm\n', ...
    gear_params.geometry.da1, gear_params.geometry.da2);
fprintf('  中心距: a=%.2fmm\n', gear_params.geometry.a);
fprintf('  啮合线长度: g_α=%.4fmm\n', gear_params.geometry.g_alpha);
fprintf('  基圆齿距: pb=%.4fmm\n', gear_params.geometry.pb);
fprintf('  端面重合度: ε_α=%.4f\n', stiffness_results.epsilon_alpha);

fprintf('\n刚度计算结果:\n');
fprintf('  赫兹接触刚度: %.4e N/mm\n', stiffness_results.hertz_stiffness);
fprintf('  齿轮1单齿刚度范围: %.2e ~ %.2e N/mm\n', ...
    min(stiffness_results.single_tooth_stiffness_1), ...
    max(stiffness_results.single_tooth_stiffness_1));
fprintf('  齿轮2单齿刚度范围: %.2e ~ %.2e N/mm\n', ...
    min(stiffness_results.single_tooth_stiffness_2), ...
    max(stiffness_results.single_tooth_stiffness_2));
fprintf('  单齿对啮合刚度范围: %.2e ~ %.2e N/mm\n', ...
    min(stiffness_results.single_pair_stiffness), ...
    max(stiffness_results.single_pair_stiffness));
fprintf('  时变啮合刚度范围: %.2e ~ %.2e N/mm\n', ...
    min(stiffness_results.time_mesh_stiffness), ...
    max(stiffness_results.time_mesh_stiffness));
fprintf('  平均啮合刚度: %.2e N/mm\n', stiffness_results.avg_stiffness);
fprintf('  ISO标准平均刚度: %.2e N/mm\n', stiffness_results.iso_stiffness);
fprintf('  刚度波动幅度: %.2f%%\n', ...
    100*(max(stiffness_results.time_mesh_stiffness) - ...
    min(stiffness_results.time_mesh_stiffness)) / stiffness_results.avg_stiffness);

fprintf('\n啮合特性:\n');
avg_contact_ratio = mean(stiffness_results.contact_ratio);
fprintf('  平均接触齿对数: %.2f\n', avg_contact_ratio);
fprintf('  单齿啮合比例: %.1f%%\n', ...
    100*sum(stiffness_results.contact_ratio == 1)/length(stiffness_results.contact_ratio));
fprintf('  双齿啮合比例: %.1f%%\n', ...
    100*sum(stiffness_results.contact_ratio == 2)/length(stiffness_results.contact_ratio));

% 保存报告到文件
report_filename = 'gear_stiffness_report.txt';
fid = fopen(report_filename, 'w');
if fid ~= -1
    fprintf(fid, '圆柱齿轮啮合刚度分析报告\n');
    fprintf(fid, '生成时间: %s\n\n', datestr(now));
    
    fprintf(fid, '1. 齿轮参数\n');
    fprintf(fid, '   齿轮1: z1=%d, m=%.2fmm, α=%.1f°, b=%.1fmm\n', ...
        gear_params.gear1.teeth, gear_params.gear1.modulus, ...
        gear_params.gear1.pressure_angle, gear_params.gear1.face_width);
    fprintf(fid, '   齿轮2: z2=%d, m=%.2fmm, α=%.1f°, b=%.1fmm\n\n', ...
        gear_params.gear2.teeth, gear_params.gear2.modulus, ...
        gear_params.gear2.pressure_angle, gear_params.gear2.face_width);
    
    fprintf(fid, '2. 几何参数\n');
    fprintf(fid, '   分度圆直径: d1=%.2fmm, d2=%.2fmm\n', ...
        gear_params.geometry.d1, gear_params.geometry.d2);
    fprintf(fid, '   基圆直径: db1=%.2fmm, db2=%.2fmm\n', ...
        gear_params.geometry.db1, gear_params.geometry.db2);
    fprintf(fid, '   端面重合度: ε_α=%.4f\n\n', stiffness_results.epsilon_alpha);
    
    fprintf(fid, '3. 刚度计算结果\n');
    fprintf(fid, '   赫兹接触刚度: %.4e N/mm\n', stiffness_results.hertz_stiffness);
    fprintf(fid, '   平均啮合刚度: %.4e N/mm\n', stiffness_results.avg_stiffness);
    fprintf(fid, '   刚度波动幅度: %.2f%%\n\n', ...
        100*(max(stiffness_results.time_mesh_stiffness) - ...
        min(stiffness_results.time_mesh_stiffness)) / stiffness_results.avg_stiffness);
    
    fprintf(fid, '4. 啮合特性\n');
    fprintf(fid, '   平均接触齿对数: %.2f\n', avg_contact_ratio);
    
    fclose(fid);
    fprintf('  报告已保存到 %s\n', report_filename);
else
    fprintf('  无法创建报告文件\n');
end

end

4. 简化版本（快速计算）

%% 圆柱齿轮啮合刚度简化计算程序
% 快速计算平均啮合刚度和时变刚度

clear; clc; close all;

% 1. 输入基本参数
m = 2;           % 模数 (mm)
z1 = 30;         % 小齿轮齿数
z2 = 50;         % 大齿轮齿数
alpha = 20;      % 压力角 (度)
b = 20;          % 齿宽 (mm)
E = 2.06e5;      % 弹性模量 (MPa)
v = 0.3;         % 泊松比

% 2. 计算重合度
d1 = m * z1;
d2 = m * z2;
db1 = d1 * cosd(alpha);
db2 = d2 * cosd(alpha);
da1 = d1 + 2 * m;
da2 = d2 + 2 * m;
a = (d1 + d2) / 2;
g_alpha = sqrt(da1^2 - db1^2)/(2*cosd(alpha)) + sqrt(da2^2 - db2^2)/(2*cosd(alpha)) - a*sind(alpha);
pb = pi * m * cosd(alpha);
epsilon_alpha = g_alpha / pb;

% 3. 计算赫兹接触刚度
k_h = (pi * E * b) / (4 * (1 - v^2));

% 4. 计算单齿刚度（简化公式）
% 经验公式: 单齿刚度与模数、齿宽、材料有关
k_single1 = 0.8 * E * b;  % 齿轮1单齿刚度 (N/mm)
k_single2 = 0.8 * E * b;  % 齿轮2单齿刚度 (N/mm)

% 5. 计算单齿对啮合刚度
k_pair = 1 / (1/k_single1 + 1/k_single2 + 1/k_h);

% 6. 计算时变刚度
n_points = 100;
pos = linspace(0, 1, n_points);
k_mesh = zeros(n_points, 1);

for i = 1:n_points
    if epsilon_alpha < 2
        if pos(i) < (epsilon_alpha - 1) || pos(i) > (2 - epsilon_alpha)
            k_mesh(i) = 2 * k_pair;  % 双齿啮合
        else
            k_mesh(i) = k_pair;      % 单齿啮合
        end
    else
        k_mesh(i) = floor(epsilon_alpha) * k_pair;
    end
end

% 7. 显示结果
fprintf('齿轮参数: m=%dmm, z1=%d, z2=%d, α=%d°, b=%dmm\n', m, z1, z2, alpha, b);
fprintf('重合度: ε_α = %.3f\n', epsilon_alpha);
fprintf('赫兹接触刚度: k_h = %.2e N/mm\n', k_h);
fprintf('单齿对啮合刚度: k_pair = %.2e N/mm\n', k_pair);
fprintf('平均啮合刚度: k_avg = %.2e N/mm\n', mean(k_mesh));
fprintf('刚度波动: %.1f%%\n', 100*(max(k_mesh)-min(k_mesh))/mean(k_mesh));

% 8. 绘图
figure;
plot(pos, k_mesh/1e3, 'b-', 'LineWidth', 2);
xlabel('啮合位置 (0:啮入, 1:啮出)');
ylabel('啮合刚度 (N/μm)');
title(sprintf('圆柱齿轮时变啮合刚度 (z1=%d, z2=%d, ε_α=%.3f)', z1, z2, epsilon_alpha));
grid on;

参考代码 基于势能原理计算圆柱齿轮啮合刚度的程序 www.3dddown.com/cnb/50789.html

5. 使用说明

5.1 程序特点

1. 基于势能原理：考虑弯曲、剪切和轴向压缩变形能
2. 考虑赫兹接触：计算齿面接触变形
3. 时变特性：考虑重合度影响，计算单双齿交替啮合的刚度变化
4. 参数化设计：可方便修改齿轮参数

5.2 主要计算步骤

1. 输入齿轮几何参数和材料属性
2. 计算齿轮基本几何尺寸
3. 基于势能法计算单齿刚度
4. 计算赫兹接触刚度
5. 计算单齿对啮合刚度（串联模型）
6. 考虑重合度，计算多齿对啮合刚度（并联模型）
7. 输出时变啮合刚度曲线

5.3 关键参数影响

• 模数增大 → 齿厚增加 → 刚度增大
• 齿宽增大 → 截面面积增大 → 刚度线性增大
• 齿数增加 → 齿形更瘦长 → 刚度可能减小
• 压力角增大 → 齿根更厚 → 刚度增大

5.4 注意事项

1. 程序中使用了一些简化假设，实际齿轮的齿形可能更复杂
2. 势能法计算中，轮齿离散切片数影响计算精度
3. 赫兹接触刚度公式适用于线接触情况
4. 实际应用中还需考虑齿轮体变形、制造误差等因素

5.5 扩展应用

1. 斜齿轮：修改为螺旋齿参数，考虑轴向重合度
2. 变位齿轮：考虑变位系数对齿形的影响
3. 修形齿轮：考虑齿廓修形对接触的影响
4. 动态分析：将刚度用于齿轮系统动力学分析