当前位置：首页 > news >正文

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

news 2026/6/13 5:19:46

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_elemnode1 = (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_elemnode1 = (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;endG = G_temp;% 求解复特征值问题：(K - ω²M + iΩG)Φ = 0 → 简化为实特征值（忽略阻尼）% 注：严格来说需求解复特征值，此处简化为无阻尼情况（G仅影响频率偏移）[V, D] = eig(K, M);  % 广义特征值问题 KΦ = ω²MΦ → ω²=Domega_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_modesf_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*GC = 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