有限元法求转子临界转速的MATLAB实现

一、转子临界转速与有限元原理

1.1 临界转速定义

转子系统旋转时，当转速达到某一特定值，系统会因共振发生剧烈振动，此转速称为临界转速。其本质是转子-支承系统的固有频率与旋转频率重合时的现象（无阻尼假设下）。

1.2 有限元建模思路

将转子离散为梁单元（Euler-Bernoulli梁或Timoshenko梁），考虑以下关键效应：

• 弯曲振动：转子横向弯曲变形（主导振动形式）；

• 陀螺效应：旋转时离心力引起的陀螺力矩（影响固有频率分裂）；

• 支承刚度：轴承简化为弹簧（提供径向支撑）。

通过求解广义特征值问题得到系统固有频率，进而确定临界转速（临界转速=固有频率，单位：rad/s 或 rpm）。

二、有限元模型与矩阵推导

2.1 梁单元动力学方程

对于旋转梁单元，考虑横向位移 w(x,t)和转角 θ(x,t)，其无阻尼自由振动方程为：

$(M+G)\ddot{q}+K_q=0$

其中：

• $$q=[w_1,θ_1,w_2,θ_2]$$T为单元节点位移向量（2节点，每节点2自由度：横向位移+转角）；
• M为质量矩阵（含平动和转动惯量）；
• G为陀螺矩阵（描述陀螺效应）；
• K为刚度矩阵（含弯曲刚度和支承刚度）。

2.2 单元矩阵公式（Euler-Bernoulli梁）

（1）刚度矩阵 $K^e$（4×4）

• E：弹性模量（Pa）；I：截面惯性矩（m⁴）；l：单元长度（m）。

（2）质量矩阵 Me（4×4，一致质量矩阵）

• ρ：材料密度（kg/m³）；A：截面面积（m²）。

（3）陀螺矩阵 Ge（4×4）

• Ω：转子旋转角速度（rad/s），注意：陀螺矩阵与旋转速度相关，求解临界转速时需迭代或近似处理（见下文）。

三、MATLAB实现步骤

3.1 主程序框架

function rotor_critical_speed()
    % 有限元法求转子临界转速

    % 1. 转子参数设置
    params = setup_rotor_parameters();

    % 2. 离散化转子（梁单元划分）
    [nodes, elements] = discretize_rotor(params);

    % 3. 组装整体刚度矩阵K、质量矩阵M、陀螺矩阵G
    [K, M, G] = assemble_matrices(nodes, elements, params);

    % 4. 施加边界条件（如简支、固支）
    [K_reduced, M_reduced, G_reduced] = apply_boundary_conditions(K, M, G, params);

    % 5. 求解广义特征值问题（考虑陀螺效应）
    [eigenvalues, eigenvectors] = solve_eigenvalue_problem(M_reduced, K_reduced, G_reduced, params);

    % 6. 计算临界转速并可视化
    critical_speeds = eigenvalues / (2*pi) * 60;  % rad/s → rpm
    visualize_results(nodes, eigenvectors, critical_speeds, params);
end

3.2 参数设置模块

function params = setup_rotor_parameters()
    % 转子几何与材料参数
    params.E = 210e9;          % 弹性模量 (Pa)，钢
    params.rho = 7800;         % 密度 (kg/m³)，钢
    params.D = 0.05;           % 转子直径 (m)
    params.L = 1.0;            % 转子总长 (m)
    params.A = pi*(params.D/2)^2;  % 截面面积 (m²)
    params.I = pi*params.D^4/64;   % 截面惯性矩 (m⁴)

    % 离散化参数
    params.ne = 10;            % 单元数（将转子分为10段）
    params.node_num = params.ne + 1;  % 节点数
    params.dof_per_node = 2;   % 每节点自由度（横向位移w, 转角θ）
    params.total_dof = params.node_num * params.dof_per_node;  % 总自由度

    % 边界条件（简支：节点1和节点ne+1的w=0, θ自由）
    params.support_nodes = [1, params.node_num];  % 支承节点编号
    params.support_type = 'simply_supported';     % 简支

    % 旋转速度（初始猜测，用于陀螺矩阵）
    params.Omega = 1000;        % 初始旋转速度 (rad/s)
end

3.3 转子离散化

function [nodes, elements] = discretize_rotor(params)
    % 生成节点坐标和单元连接关系
    nodes = zeros(params.node_num, 1);  % 节点坐标（x轴，m）
    for i = 1:params.node_num
        nodes(i) = (i-1)*params.L/params.ne;  % 均匀分布
    end

    elements = zeros(params.ne, 2);  % 单元连接（节点i, 节点j）
    for i = 1:params.ne
        elements(i, :) = [i, i+1];
    end
end

3.4 矩阵组装模块

function [K, M, G] = assemble_matrices(nodes, elements, params)
    % 组装整体刚度矩阵K、质量矩阵M、陀螺矩阵G
    total_dof = params.total_dof;
    K = zeros(total_dof, total_dof);
    M = zeros(total_dof, total_dof);
    G = zeros(total_dof, total_dof);

    % 遍历所有单元
    for e = 1:size(elements, 1)
        node_i = elements(e, 1);  % 单元起点节点
        node_j = elements(e, 2);  % 单元终点节点
        l = nodes(node_j) - nodes(node_i);  % 单元长度 (m)

        % 1. 单元刚度矩阵K^e
        Ke = beam_stiffness_matrix(params.E, params.I, l);

        % 2. 单元质量矩阵M^e（一致质量矩阵）
        Me = beam_mass_matrix(params.rho, params.A, params.I, l);

        % 3. 单元陀螺矩阵G^e（与旋转速度Ω相关）
        Ge = beam_gyroscopic_matrix(params.rho, params.I, l, params.Omega);

        % 4. 自由度映射（节点i: dof 2i-1=w, 2i=θ；节点j: dof 2j-1=w, 2j=θ）
        dofs = [2*node_i-1, 2*node_i, 2*node_j-1, 2*node_j];  % 单元自由度索引

        % 5. 组装到整体矩阵
        K(dofs, dofs) = K(dofs, dofs) + Ke;
        M(dofs, dofs) = M(dofs, dofs) + Me;
        G(dofs, dofs) = G(dofs, dofs) + Ge;
    end
end

% 梁单元刚度矩阵（4×4）
function Ke = beam_stiffness_matrix(E, I, l)
    Ke = (E*I/l^3) * [12, 6*l, -12, 6*l;
                     6*l, 4*l^2, -6*l, 2*l^2;
                     -12, -6*l, 12, -6*l;
                     6*l, 2*l^2, -6*l, 4*l^2];
end

% 梁单元质量矩阵（4×4，一致质量）
function Me = beam_mass_matrix(rho, A, I, l)
    term1 = rho*A*l/420 * [156, 22*l, 54, -13*l;
                          22*l, 4*l^2, 13*l, -3*l^2;
                          54, 13*l, 156, -22*l;
                          -13*l, -3*l^2, -22*l, 4*l^2];
    term2 = rho*I*l/30 * [36, 3*l, -36, 3*l;
                          3*l, 4*l^2, -3*l, -l^2;
                          -36, -3*l, 36, -3*l;
                          3*l, -l^2, -3*l, 4*l^2];
    Me = term1 + term2;
end

% 梁单元陀螺矩阵（4×4）
function Ge = beam_gyroscopic_matrix(rho, I, l, Omega)
    Ge = (rho*I*Omega/(30*l)) * [36, 3*l, -36, 3*l;
                                3*l, 4*l^2, -3*l, -l^2;
                                -36, -3*l, 36, -3*l;
                                3*l, -l^2, -3*l, 4*l^2];
end

3.5 边界条件施加

以简支边界为例（节点1和节点n的横向位移w=0，转角θ自由）：

function [K_red, M_red, G_red] = apply_boundary_conditions(K, M, G, params)
    % 简支边界：节点1和节点n的w=0（自由度1和2n-1）
    fixed_dofs = [1, 2*params.node_num-1];  % 固定自由度索引（w1=0, wn=0）
    all_dofs = 1:params.total_dof;
    free_dofs = setdiff(all_dofs, fixed_dofs);  % 释放自由度

    % 缩减矩阵（仅保留自由自由度）
    K_red = K(free_dofs, free_dofs);
    M_red = M(free_dofs, free_dofs);
    G_red = G(free_dofs, free_dofs);
end

3.6 特征值求解（考虑陀螺效应）

陀螺效应导致系统矩阵变为非对称，需求解二次特征值问题：

$det(λ^2M+λG+K)=0$

其中 λ=jω（ω为固有频率）。采用多项式特征值求解器（如MATLAB的 polyeig）：

function [eigenvalues, eigenvectors] = solve_eigenvalue_problem(M, K, G, params)
    % 求解二次特征值问题：λ²M + λG + K = 0
    A0 = K;          % 常数项矩阵
    A1 = G;          % 一次项矩阵
    A2 = M;          % 二次项矩阵

    % 调用多项式特征值求解器
    [eigenvectors, eigenvalues] = polyeig(A0, A1, A2);

    % 提取正实部特征值（物理意义：固有频率）
    valid_idx = real(eigenvalues) > 0 & imag(eigenvalues) == 0;
    eigenvalues = eigenvalues(valid_idx);
    eigenvectors = eigenvectors(:, valid_idx);

    % 按升序排序
    [eigenvalues, idx] = sort(eigenvalues);
    eigenvectors = eigenvectors(:, idx);
end

3.7 结果可视化

function visualize_results(nodes, eigenvectors, critical_speeds, params)
    % 可视化振型与临界转速
    figure('Position', [100, 100, 1200, 600]);

    % 1. 临界转速表格
    subplot(1,2,1);
    column_names = {
   '阶次', '固有频率 (rad/s)', '临界转速 (rpm)'};
    row_names = arrayfun(@(x) sprintf('%d', x), 1:length(critical_speeds), 'UniformOutput', false);
    table_data = [row_names', num2cell([eigenvalues, critical_speeds])];
    uitable('Data', table_data, 'ColumnName', column_names, 'RowName', []);
    title('转子临界转速计算结果');

    % 2. 一阶振型（横向位移）
    subplot(1,2,2);
    first_mode = eigenvectors(:, 1);  % 一阶振型
    w = first_mode(1:2:end);  % 提取横向位移（w1, w2, ..., wn）
    plot(nodes, w, 'b-o', 'LineWidth', 2);
    xlabel('轴向位置 (m)'); ylabel('横向位移 (m)');
    title(sprintf('一阶振型 (临界转速: %.0f rpm)', critical_speeds(1)));
    grid on;
end

四、示例与结果分析

4.1 示例参数

• 转子：钢质细长转子，直径50mm，长度1m，两端简支；

• 离散化：10个梁单元（11个节点）；

• 材料：E=210GPa，ρ=7800kg/m³。

4.2 计算结果

4.3 振型分析

• 一阶振型：转子中部振幅最大（类似简支梁一阶弯曲振型）；

• 二阶振型：转子呈“S”形，两处反弯点；

• 三阶振型：转子呈“M”形，三处反弯点。

参考代码 有限元求转子临界转速 www.youwenfan.com/contentalh/60046.html

五、关键问题与优化

5.1 陀螺效应的影响

高速旋转时，陀螺效应会导致正向涡动（与旋转同向）和反向涡动（与旋转反向）频率分裂，需在特征值问题中保留陀螺矩阵 G。

5.2 阻尼的影响

实际转子存在阻尼（如轴承油膜阻尼），需在方程中加入阻尼矩阵 C，求解复特征值问题（阻尼固有频率）。

5.3 优化方向

• 高阶单元：采用Timoshenko梁单元（考虑剪切变形）；

• 非线性支承：轴承刚度随转速变化（如油膜轴承）；

• 并行计算：大规模转子系统（如汽轮机转子）需用稀疏矩阵求解器（eigs）。

六、总结

本MATLAB程序通过有限元法建立了转子-梁单元的动力学模型，考虑了弯曲刚度、陀螺效应和简支边界条件，求解广义特征值问题得到临界转速。代码模块化设计，可扩展至复杂转子系统（如多盘转子、弹性支承），为旋转机械的振动分析与平衡设计提供理论工具。