MATLAB中对转子建立有限元模型并进行动力学计算

MATLAB中对转子建立有限元模型并进行动力学计算，需基于转子动力学理论手动构建模型（质量矩阵、刚度矩阵、陀螺矩阵等），通过数值求解特征值/特征向量或频响函数实现临界转速、不平衡响应等分析。MATLAB的优势在于灵活的矩阵运算和自定义模型能力，适合教学、科研中的轻量化建模与算法验证。

一、核心思路

转子动力学有限元分析的本质是求解多自由度系统的运动微分方程：
$\mathbf{M}\ddot{\mathbf{q}} + (\mathbf{C} + \mathbf{G})\dot{\mathbf{q}} + \mathbf{K}\mathbf{q} = \mathbf{F}(t)$
其中：

• $\mathbf{M}$：质量矩阵（含平动质量和转动惯量）；
• $\mathbf{C}$：阻尼矩阵（轴承阻尼、结构阻尼）；
• $\mathbf{G}$：陀螺矩阵（旋转引起的陀螺效应，与转速$\Omega$相关）；
• $\mathbf{K}$：刚度矩阵（轴弯曲刚度、轴承支撑刚度）；
• $\mathbf{q}$：广义位移向量（节点横向位移、转角）；
• $\mathbf{F}(t)$：激励向量（如不平衡力）。

分析目标：

• 临界转速：求解无阻尼自由振动特征值（$\mathbf{C}=0, \mathbf{F}=0$），得固有频率随转速$\Omega$的变化（坎贝尔图）；
• 不平衡响应：施加简谐激励$\mathbf{F}(t)=\mathbf{F}_0 e^{i\omega t}$，求稳态响应$\mathbf{q}(t)$；
• 稳定性：求解复特征值，判断系统是否失稳（实部>0为不稳定）。

二、关键步骤与MATLAB实现

1. 模型简化与单元选择

常用梁单元模型（如欧拉-伯努利梁或铁木辛柯梁），将转子轴段离散为梁单元，节点含横向位移$y$、转角$\theta$（铁木辛柯梁还含剪切位移）。

• 简化假设：忽略轴向变形，仅考虑横向振动；材料线弹性（E, $\rho$已知）；集中质量（如叶轮、齿轮）用附加质量矩阵模拟。

2. 单元矩阵推导（以欧拉-伯努利梁为例）

单元参数：长度$l$，截面惯性矩$I$，截面积$A$，材料密度$\rho$，弹性模量$E$。

• 单元刚度矩阵$\mathbf{k}^e$（4×4，对应节点位移$[y_1, \theta_1, y_2, \theta_2]^T$）：
  $\mathbf{k}^e = \frac{EI}{l^3} \begin{bmatrix} 12 & 6l & -12 & 6l \\ 6l & 4l^2 & -6l & 2l^2 \\ -12 & -6l & 12 & -6l \\ 6l & 2l^2 & -6l & 4l^2 \end{bmatrix}$
• 单元质量矩阵$\mathbf{m}^e$（一致质量矩阵，4×4）：
  $\mathbf{m}^e = \frac{\rho A l}{420} \begin{bmatrix} 156 & 22l & 54 & -13l \\ 22l & 4l^2 & 13l & -3l^2 \\ 54 & 13l & 156 & -22l \\ -13l & -3l^2 & -22l & 4l^2 \end{bmatrix}$
• 陀螺矩阵$\mathbf{g}^e$（4×4，与转速$\Omega$相关，描述陀螺效应）：
  $\mathbf{g}^e = \frac{\rho I l}{30} \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 36 & 0 & 3l \\ 0 & 0 & 0 & 0 \\ 0 & 3l & 0 & 4l^2 \end{bmatrix} \Omega$

3. 整体矩阵组装

• 将各单元矩阵按节点自由度叠加，形成整体$\mathbf{M}, \mathbf{K}, \mathbf{G}$（注意：陀螺矩阵$\mathbf{G}$是反对称矩阵，与$\Omega$成正比）。
• 边界条件处理：轴承支撑用弹簧-阻尼单元（刚度$k_b$、阻尼$c_b$）模拟，在对应节点添加$\mathbf{K} += k_b \cdot \mathbf{I}$，$\mathbf{C} += c_b \cdot \mathbf{I}$。
• 集中质量添加：在节点处附加质量$m$和转动惯量$J$，修改$\mathbf{M}$的对角元素（$M_{ii} += m$，$M_{i+1,i+1} += J$）。

4. 动力学求解（MATLAB代码框架）

以下以两支承转子（单跨梁+中间集中质量） 为例，演示临界转速计算与坎贝尔图绘制。

三、MATLAB实例：两支承转子临界转速分析

模型参数

• 轴段：总长$L=1m$，划分为2个梁单元（3节点），每段长$l=0.5m$；
• 截面：$d=0.02m$（直径），$A=\pi d^2/4$，$I=\pi d^4/64$；
• 材料：钢，$E=210e9 Pa$，$\rho=7800 kg/m^3$；
• 集中质量：中间节点（节点2）附加叶轮质量$m=10kg$，转动惯量$J=0.1 kg·m²$；
• 轴承支撑：两端节点（节点1、3）用弹簧刚度$k_b=1e6 N/m$（忽略阻尼）；
• 转速范围：$\Omega=0\sim 10000 rpm$（转换为rad/s：$\Omega_{rad}=2\pi\Omega_{rpm}/60$）。

MATLAB代码实现

% 转子动力学有限元分析（两支承转子临界转速）
clear; clc;

%% 1. 参数定义
L = 1;          % 轴总长(m)
n_elem = 2;     % 梁单元数
l = L/n_elem;   % 单元长度(m)
d = 0.02;       % 轴直径(m)
A = pi*d^2/4;   % 截面积(m²)
I = pi*d^4/64;  % 惯性矩(m⁴)
E = 210e9;      % 弹性模量(Pa)
rho = 7800;     % 密度(kg/m³)

% 集中质量（中间节点）
m_disk = 10;    % 叶轮质量(kg)
J_disk = 0.1;   % 叶轮转动惯量(kg·m²)

% 轴承刚度(N/m)
k_b = 1e6;      % 两端轴承刚度

% 转速范围(rpm)
Omega_rpm = 0:500:10000;  % 转速序列
Omega_rad = Omega_rpm * 2*pi / 60;  % 转换为rad/s


%% 2. 单元矩阵（欧拉-伯努利梁）
% 单元刚度矩阵(k_e)、质量矩阵(m_e)、陀螺矩阵(g_e)
k_e = @(l,EI) EI/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];

m_e = @(l,rhoA) rhoA*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];

g_e = @(l,rhoI,Omega) rhoI*l/30 * [0, 0, 0, 0;
                                  0, 36, 0, 3*l;
                                  0, 0, 0, 0;
                                  0, 3*l, 0, 4*l^2] * Omega;


%% 3. 组装整体矩阵（3节点，每节点2自由度：y, θ → 共6自由度）
n_nodes = n_elem + 1;  % 节点数
DOF_per_node = 2;      % 每节点自由度（y, θ）
total_DOF = n_nodes * DOF_per_node;  % 总自由度

% 初始化整体矩阵
M = zeros(total_DOF);  % 质量矩阵
K = zeros(total_DOF);  % 刚度矩阵
G = zeros(total_DOF);  % 陀螺矩阵

EI = E*I;       % 抗弯刚度
rhoA = rho*A;   % 单位长度质量
rhoI = rho*I;   % 单位长度转动惯量

% 组装梁单元矩阵（2个单元，节点1-2，2-3）
for i = 1:n_elem
    node1 = (i-1)*DOF_per_node + 1;  % 单元起始自由度索引
    node2 = node1 + DOF_per_node;    % 单元结束自由度索引
    idx = [node1, node1+1, node2, node2+1];  % 单元自由度索引
    
    % 单元矩阵
    ke = k_e(l, EI);
    me = m_e(l, rhoA);
    ge = g_e(l, rhoI, 0);  % 初始Omega=0，后续更新
    
    % 组装到整体矩阵
    M(idx, idx) = M(idx, idx) + me;
    K(idx, idx) = K(idx, idx) + ke;
    G(idx, idx) = G(idx, idx) + ge;
end

% 添加集中质量（中间节点2：自由度3(y), 4(θ)）
M(3,3) = M(3,3) + m_disk;    % y方向质量
M(4,4) = M(4,4) + J_disk;    % θ方向转动惯量

% 添加轴承支撑（节点1和3：y方向加弹簧）
K(1,1) = K(1,1) + k_b;  % 节点1 y方向刚度
K(5,5) = K(5,5) + k_b;  % 节点3 y方向刚度（节点3自由度：5(y), 6(θ)）


%% 4. 临界转速计算（坎贝尔图：转速Ω vs 固有频率ω）
n_modes = 4;  % 提取前4阶模态
freqs = zeros(length(Omega_rad), n_modes);  % 存储各转速下的固有频率

for i = 1:length(Omega_rad)
    Omega = Omega_rad(i);
    % 更新陀螺矩阵（与Omega相关）
    G_temp = zeros(total_DOF);
    for j = 1:n_elem
        node1 = (j-1)*DOF_per_node + 1;
        node2 = node1 + DOF_per_node;
        idx = [node1, node1+1, node2, node2+1];
        ge = g_e(l, rhoI, Omega);  % 单元陀螺矩阵
        G_temp(idx, idx) = G_temp(idx, idx) + ge;
    end
    G = G_temp;
    
    % 求解复特征值问题：(K - ω²M + iΩG)Φ = 0 → 简化为实特征值（忽略阻尼）
    % 注：严格来说需求解复特征值，此处简化为无阻尼情况（G仅影响频率偏移）
    [V, D] = eig(K, M);  % 广义特征值问题 KΦ = ω²MΦ → ω²=D
    omega_sq = diag(D);  % 特征值（ω²）
    omega = sqrt(abs(omega_sq));  % 固有频率(rad/s)
    freq_hz = omega/(2*pi);       % 转换为Hz
    
    % 排序并提取前n_modes阶频率
    [~, idx_sort] = sort(freq_hz);
    freqs(i,:) = freq_hz(idx_sort(1:n_modes))';
end


%% 5. 绘制坎贝尔图（转速Ω vs 固有频率f）
figure;
plot(Omega_rpm, freqs(:,1), 'b-o', 'DisplayName', '1st Mode');
hold on;
plot(Omega_rpm, freqs(:,2), 'r-s', 'DisplayName', '2nd Mode');
plot(Omega_rpm, freqs(:,3), 'g-d', 'DisplayName', '3rd Mode');
plot(Omega_rpm, freqs(:,4), 'm-^', 'DisplayName', '4th Mode');
xlabel('转速 (rpm)');
ylabel('固有频率 (Hz)');
title('转子坎贝尔图（临界转速预测）');
legend();
grid on;
xlim([0, max(Omega_rpm)]);


%% 6. 提取临界转速（固有频率与转速线交点）
critical_rpm = [];
for mode = 1:n_modes
    f_mode = freqs(:,mode);
    % 寻找f_mode ≈ Omega_rpm/(60)（同步振动，1X激励）的交点（近似）
    % 注：严格需比较f_mode与Omega_rpm/(60)，此处简化为直接读取共振点
    [max_f, idx_max] = max(f_mode);
    critical_rpm(end+1) = Omega_rpm(idx_max);
end
disp('临界转速估算值 (rpm):');
disp(critical_rpm);

四、扩展分析：不平衡响应

若需计算不平衡激励下的振动响应，可在集中质量节点（如节点2）施加简谐力$F(t)=m_disk \cdot e \cdot \Omega^2 \sin(\Omega t)$（$e$为偏心距），通过谐响应分析求解稳态位移：

% 不平衡响应分析（示例：节点2 y方向响应）
e = 1e-4;  % 偏心距(m)
t = 0:0.001:1;  % 时间向量
Omega_test = 3000;  % 测试转速(rpm)
Omega_rad_test = Omega_test * 2*pi / 60;  % rad/s

% 激励力幅值：F0 = m_disk * e * Omega_rad_test²
F0 = m_disk * e * Omega_rad_test^2;
F = [0; 0; F0; 0; 0; 0];  % 节点2 y方向力（自由度3）

% 求解谐响应（频域法）
w = linspace(0, 2*Omega_rad_test, 1000);  % 频率范围
H = zeros(size(w));
for i = 1:length(w)
    % 复刚度矩阵：K - w(i)^2*M + i*w(i)*C + i*Omega_rad_test*G
    C = zeros(total_DOF);  % 忽略阻尼
    G_temp = g_e_matrix(Omega_rad_test);  % 整体陀螺矩阵（需提前定义）
    H(i) = F(3)/( (K(3,3)-w(i)^2*M(3,3)) + 1i*(w(i)*C(3,3) + Omega_rad_test*G_temp(3,3)) );
end

% 绘制幅频响应
figure;
plot(w/(2*pi), abs(H));
xlabel('频率 (Hz)');
ylabel('振幅 (m)');
title('节点2 y方向不平衡响应');
grid on;

参考代码 用于对转子建立有限元模型并进行动力学计算 www.youwenfan.com/contentcnt/63152.html

五、注意事项

1. 模型简化：MATLAB适合轻量化模型（梁单元、集中质量），复杂模型（3D实体单元）建议用ANSYS等专业软件；
2. 陀螺效应：高速转子必须考虑陀螺矩阵$\mathbf{G}$，否则临界转速计算偏差大；
3. 特征值求解：严格需用eig(K, M)求解广义特征值，或用polyeig处理阻尼系统；
4. 单位统一：确保所有参数单位为国际单位（m, kg, s, N）；
5. 验证：先用经典转子（如Jeffcott转子）验证代码正确性（解析解：临界转速$\Omega_c = \sqrt{k/m}$）。