% 二叶桨动特性分析(双桨叶五剖面有限元模型),考虑离心力影响
% 旋翼总直径:0.81534m,桨毂直径:50.2mm,自转转速:3800RPM
clear; clc; close all;
%% 1. 基本参数定义
% 几何参数
rotor_diameter = 0.81534; % 旋翼总直径(m,从一侧桨尖到另一侧)
hub_diameter = 0.0502; % 桨毂直径(m)
hub_radius = hub_diameter / 2; % 桨毂半径(m)
% 单桨叶实际长度 = (旋翼直径 - 桨毂直径) / 2
L_blade = (rotor_diameter - hub_diameter) / 2;
n_elem_per_blade = 5; % 每片桨叶单元数量
n_node_per_blade = n_elem_per_blade + 1; % 每片桨叶节点数量
n_blades = 2; % 二叶桨
elem_len = L_blade / n_elem_per_blade; % 每个单元长度(m)
% 旋转参数
RPM = 3800; % 自转转速 (RPM)
Omega = RPM * 2 * pi / 60; % 角速度 (rad/s)
fprintf('旋转角速度: %.2f rad/s (%.0f RPM)\n', Omega, RPM);
% 节点编号方案(总节点数11):
% - 节点1:公共桨毂中心
% - 桨叶1:节点1, 2, 3, 4, 5, 6
% - 桨叶2:节点1, 7, 8, 9, 10, 11
total_nodes = 1 + (n_node_per_blade - 1) * n_blades;
% 存储两片桨叶的节点位置坐标(从桨毂边缘开始)
node_pos = cell(n_blades, 1);
for b = 1:n_blades
% 节点位置从桨毂边缘(hub_radius)到桨尖(hub_radius + L_blade)
node_pos{b} = linspace(hub_radius, hub_radius + L_blade, n_node_per_blade);
end
% 显示几何参数
disp('=== 几何参数 ===');
fprintf('旋翼总直径: %.5f m\n', rotor_diameter);
fprintf('桨毂直径: %.5f m\n', hub_diameter);
fprintf('单桨叶实际长度: %.5f m\n', L_blade);
fprintf('桨叶起始位置(距中心): %.5f m\n', hub_radius);
fprintf('桨尖位置(距中心): %.5f m\n', hub_radius + L_blade);
fprintf('单元长度: %.5f m\n', elem_len);
%% 2. 导入各剖面质量矩阵(6x6,DYMORE convention)
M_profiles = cell(n_elem_per_blade, 1);
% 单元1(80mm处剖面 - 从桨毂边缘算起)
M_profiles{1} = [
0.289079e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.173399e-13, -0.853452e-03;
0.000000e+00, 0.289079e+00, 0.000000e+00, -0.173399e-13, 0.000000e+00, 0.000000e+00;
0.000000e+00, 0.000000e+00, 0.289079e+00, 0.853452e-03, 0.000000e+00, 0.000000e+00;
0.000000e+00, -0.173399e-13, 0.853452e-03, 0.533859e-04, 0.000000e+00, 0.000000e+00;
0.173399e-13, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.461949e-06, -0.585358e-15;
-0.853452e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, -0.585358e-15, 0.529240e-04;
];
% 单元2(160mm处剖面)
M_profiles{2} = [
0.530206e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -0.390687e-18, -0.211993e-02;
0.000000e+00, 0.530206e+00, 0.000000e+00, 0.390687e-18, 0.000000e+00, 0.000000e+00;
0.000000e+00, 0.000000e+00, 0.530206e+00, 0.211993e-02, 0.000000e+00, 0.000000e+00;
0.000000e+00, 0.390687e-18, 0.211993e-02, 0.179590e-03, 0.000000e+00, 0.000000e+00;
-0.390687e-18, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.155400e-05, 0.222282e-19;
-0.211993e-02, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.222282e-19, 0.178036e-03;
];
% 单元3(240mm处剖面)
M_profiles{3} = [
0.407341e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -0.480828e-13, -0.142755e-02;
0.000000e+00, 0.407341e+00, 0.000000e+00, 0.480828e-13, 0.000000e+00, 0.000000e+00;
0.000000e+00, 0.000000e+00, 0.407341e+00, 0.142755e-02, 0.000000e+00, 0.000000e+00;
0.000000e+00, 0.480828e-13, 0.142755e-02, 0.106001e-03, 0.000000e+00, 0.000000e+00;
-0.480828e-13, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.917227e-06, 0.197552e-14;
-0.142755e-02, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.197552e-14, 0.105084e-03;
];
% 单元4(320mm处剖面)
M_profiles{4} = [
0.216953e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -0.655954e-19, -0.554883e-03;
0.000000e+00, 0.216953e+00, 0.000000e+00, 0.655954e-19, 0.000000e+00, 0.000000e+00;
0.000000e+00, 0.000000e+00, 0.216953e+00, 0.554883e-03, 0.000000e+00, 0.000000e+00;
0.000000e+00, 0.655954e-19, 0.554883e-03, 0.300693e-04, 0.000000e+00, 0.000000e+00;
-0.655954e-19, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.260190e-06, 0.528450e-20;
-0.554883e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.528450e-20, 0.298091e-04;
];
% 单元5(400mm处剖面)
M_profiles{5} = [
0.465444e-01, 0.000000e+00, 0.000000e+00, 0.000000e+00, -0.898566e-20, -0.551386e-04;
0.000000e+00, 0.465444e-01, 0.000000e+00, 0.898566e-20, 0.000000e+00, 0.000000e+00;
0.000000e+00, 0.000000e+00, 0.465444e-01, 0.551386e-04, 0.000000e+00, 0.000000e+00;
0.000000e+00, 0.898566e-20, 0.551386e-04, 0.138398e-05, 0.000000e+00, 0.000000e+00;
-0.898566e-20, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.119756e-07, 0.204937e-21;
-0.551386e-04, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.204937e-21, 0.137200e-05;
];
%% 3. 导入各剖面刚度矩阵(6x6,Timoshenko梁模型)
K_profiles = cell(n_elem_per_blade, 1);
% 单元1(80mm处剖面)
K_profiles{1} = [
0.38543891E+08, -0.66963691E-20, -0.56273690E-23, 0.61213913E-25, -0.31799568E-05, -0.11379358E+06;
-0.66963691E-20, 0.12953125E+08, -0.21620513E+00, 0.94822720E-03, 0.60167247E-25, 0.54036883E-22;
-0.56273690E-23, -0.21620513E+00, 0.24713367E+07, -0.44623520E+04, 0.28654185E-24, -0.14632065E-24;
0.61213913E-25, 0.94822720E-03, -0.44623520E+04, 0.10072960E+03, 0.42614533E-25, -0.24468768E-26;
-0.31799568E-05, 0.60167247E-25, 0.28654185E-24, 0.42614533E-25, 0.61593281E+02, 0.33417005E-05;
-0.11379358E+06, 0.54036883E-22, -0.14632065E-24, -0.24468768E-26, 0.33417005E-05, 0.70565324E+04;
];
% 单元2(160mm处剖面)
K_profiles{2} = [
0.70694289E+08, -0.87920645E-20, -0.33253495E-23, 0.50586537E-25, -0.55380916E-07, -0.28265461E+06;
-0.87920645E-20, 0.23757660E+08, -0.30709474E+00, 0.45802332E-02, -0.39530577E-25, 0.74789743E-22;
-0.33253495E-23, -0.30709474E+00, 0.45327500E+07, -0.11084295E+05, 0.72548387E-24, 0.65879888E-26;
0.50586537E-25, 0.45802332E-02, -0.11084295E+05, 0.33885503E+03, 0.11432801E-24, 0.15012058E-27;
-0.55380916E-07, -0.39530577E-25, 0.72548387E-24, 0.11432801E-24, 0.20720002E+03, 0.73770970E-05;
-0.28265461E+06, 0.74789743E-22, 0.65879888E-26, 0.15012058E-27, 0.73770970E-05, 0.23738259E+05;
];
% 单元3(240mm处剖面)
K_profiles{3} = [
0.54312139E+08, -0.83696149E-20, 0.62221045E-23, 0.63188965E-25, 0.18737318E-04, -0.19033989E+06;
-0.83696149E-20, 0.18252228E+08, 0.48802966E-02, -0.68435025E-03, 0.14045370E-25, 0.84467105E-22;
0.62221045E-23, 0.48802966E-02, 0.34823399E+07, -0.74640662E+04, 0.48828173E-24, -0.46403406E-24;
0.63188965E-25, -0.68435025E-03, -0.74640662E+04, 0.20000448E+03, 0.74856810E-25, -0.49685200E-26;
0.18737318E-04, 0.14045370E-25, 0.48828173E-24, 0.74856810E-25, 0.12229702E+03, 0.32527444E-05;
-0.19033989E+06, 0.84467105E-22, -0.46403406E-24, -0.49685200E-26, 0.32527444E-05, 0.14011158E+05;
];
% 单元4(320mm处剖面)
K_profiles{4} = [
0.28927020E+08, -0.55019423E-20, -0.19951886E-23, -0.31680102E-25, 0.20243815E-07, -0.73984363E+05;
-0.55019423E-20, 0.97212628E+07, 0.97022513E-01, -0.92168379E-04, 0.13735065E-25, 0.37530444E-22;
-0.19951886E-23, 0.97022513E-01, 0.18547260E+07, -0.29012603E+04, 0.17499720E-24, 0.24725797E-25;
-0.31680102E-25, -0.92168379E-04, -0.29012603E+04, 0.56735317E+02, 0.26440805E-25, 0.51927984E-27;
0.20243815E-07, 0.13735065E-25, 0.17499720E-24, 0.26440805E-25, 0.34692013E+02, -0.31423757E-06;
-0.73984363E+05, 0.37530444E-22, 0.24725797E-25, 0.51927984E-27, -0.31423757E-06, 0.39745468E+04;
];
% 单元5(400mm处剖面)
K_profiles{5} = [
0.62059240E+07, -0.11241305E-20, 0.12818189E-24, -0.17512347E-27, 0.53922014E-07, -0.73518077E+04;
-0.11241305E-20, 0.20855732E+07, -0.13819523E-02, 0.24279922E-06, 0.42246092E-27, 0.33534716E-23;
0.12818189E-24, -0.13819523E-02, 0.39790589E+06, -0.28829546E+03, 0.13609137E-25, 0.71224377E-27;
-0.17512347E-27, 0.24279922E-06, -0.28829546E+03, 0.26113084E+01, 0.12534636E-26, 0.17057581E-29;
0.53922014E-07, 0.42246092E-27, 0.13609137E-25, 0.12534636E-26, 0.15967422E+01, 0.40783515E-08;
-0.73518077E+04, 0.33534716E-23, 0.71224377E-27, 0.17057581E-29, 0.40783515E-08, 0.18293344E+03;
];
%% 4. 计算离心力引起的几何刚度矩阵(初应力刚度矩阵)
K_geo_profiles = cell(n_elem_per_blade, 1);
for i = 1:n_elem_per_blade
% 从质量矩阵提取单元质量信息
mass_elem = M_profiles{i}(1,1); % 轴向质量
Iy_elem = M_profiles{i}(5,5); % 绕y轴转动惯量
Iz_elem = M_profiles{i}(6,6); % 绕z轴转动惯量
% 单元位置(距旋转中心)- 取单元中点位置
nodeA_pos = node_pos{1}(i);
nodeB_pos = node_pos{1}(i+1);
r_mid = (nodeA_pos + nodeB_pos) / 2; % 单元中点半径
% 计算离心力: F = m * r * Ω²
F_centrifugal = mass_elem * r_mid * Omega^2;
% 构建Timoshenko梁单元的几何刚度矩阵(考虑离心力)
% 参考旋转梁理论推导的几何刚度矩阵形式
K_geo_elem = zeros(6,6);
% 轴向刚度增强(由离心力引起)
K_geo_elem(1,1) = F_centrifugal / elem_len;
% 弯曲刚度增强(考虑转动惯量影响)
K_geo_elem(5,5) = (Iy_elem * Omega^2) / elem_len; % 绕y轴
K_geo_elem(6,6) = (Iz_elem * Omega^2) / elem_len; % 绕z轴
% 剪切刚度影响
K_geo_elem(2,2) = 0.1 * F_centrifugal / elem_len; % 经验系数
K_geo_elem(3,3) = 0.1 * F_centrifugal / elem_len; % 经验系数
% 扭转刚度影响
K_geo_elem(4,4) = 0.05 * F_centrifugal * r_mid / elem_len; % 经验系数
K_geo_profiles{i} = K_geo_elem;
end
% 显示刚度矩阵特性
disp('=== 刚度矩阵特性 ===');
for i = 1:n_elem_per_blade
fprintf('单元%d - 结构刚度最大项: %.3e, 几何刚度最大项: %.3e\n', ...
i, max(max(abs(K_profiles{i}))), max(max(abs(K_geo_profiles{i}))));
end
%% 5. 组装整体质量矩阵、刚度矩阵和几何刚度矩阵
DOF_per_node = 6; % 每个节点6个自由度
DOF_total = DOF_per_node * total_nodes; % 总自由度:11节点×6=66
M_global = zeros(DOF_total); % 整体质量矩阵
K_global = zeros(DOF_total); % 整体刚度矩阵
K_geo_global = zeros(DOF_total); % 整体几何刚度矩阵
% 为每片桨叶定义完整节点列表(包含根部)
blade_nodes = {
[1, 2, 3, 4, 5, 6]; % 桨叶1节点
[1, 7, 8, 9, 10, 11]; % 桨叶2节点
};
% 循环组装两片桨叶的单元
for b = 1:n_blades % 桨叶编号:1(左)、2(右)
for i = 1:n_elem_per_blade % 每个桨叶的5个单元
% 从节点列表中获取当前单元的两个节点
nodeA = blade_nodes{b}(i);
nodeB = blade_nodes{b}(i+1);
% 计算自由度索引
dofA = (nodeA - 1)*DOF_per_node + 1 : nodeA*DOF_per_node;
dofB = (nodeB - 1)*DOF_per_node + 1 : nodeB*DOF_per_node;
% 检查索引范围
if any(dofA > DOF_total) || any(dofB > DOF_total)
error('自由度索引超出范围: 桨叶=%d, 单元=%d, nodeA=%d, nodeB=%d', b, i, nodeA, nodeB);
end
% 获取单元矩阵
M_elem = M_profiles{i};
K_elem = K_profiles{i} / elem_len; % 结构刚度矩阵
K_geo_elem = K_geo_profiles{i}; % 几何刚度矩阵
% 组装质量矩阵
M_global(dofA, dofA) = M_global(dofA, dofA) + M_elem;
M_global(dofB, dofB) = M_global(dofB, dofB) + M_elem;
% 组装结构刚度矩阵
K_global(dofA, dofA) = K_global(dofA, dofA) + K_elem;
K_global(dofB, dofB) = K_global(dofB, dofB) + K_elem;
K_global(dofA, dofB) = K_global(dofA, dofB) - K_elem;
K_global(dofB, dofA) = K_global(dofB, dofA) - K_elem;
% 组装几何刚度矩阵
K_geo_global(dofA, dofA) = K_geo_global(dofA, dofA) + K_geo_elem;
K_geo_global(dofB, dofB) = K_geo_global(dofB, dofB) + K_geo_elem;
K_geo_global(dofA, dofB) = K_geo_global(dofA, dofB) - K_geo_elem;
K_geo_global(dofB, dofA) = K_geo_global(dofB, dofA) - K_geo_elem;
end
end
%% 6. 应用边界条件(桨毂中心固定)
fixed_dof = 1:DOF_per_node; % 节点1(桨毂中心)的6个自由度全部固定
free_dof = setdiff(1:DOF_total, fixed_dof); % 自由自由度
% 缩减矩阵
K_red = K_global(free_dof, free_dof);
M_red = M_global(free_dof, free_dof);
K_geo_red = K_geo_global(free_dof, free_dof);
% 组合刚度矩阵(结构刚度 + 几何刚度)
K_combined_red = K_red + K_geo_red;
%% 7. 求解固有频率和振型(考虑与不考虑离心力两种情况)
% 情况1: 不考虑离心力
[V_no_cent, D_no_cent] = eig(K_red, M_red);
lambda_no_cent = diag(D_no_cent);
frequencies_no_cent = sqrt(abs(lambda_no_cent)) / (2*pi);
[frequencies_no_cent, idx_no_cent] = sort(frequencies_no_cent);
V_no_cent = V_no_cent(:, idx_no_cent);
% 情况2: 考虑离心力
[V_with_cent, D_with_cent] = eig(K_combined_red, M_red);
lambda_with_cent = diag(D_with_cent);
frequencies_with_cent = sqrt(abs(lambda_with_cent)) / (2*pi);
[frequencies_with_cent, idx_with_cent] = sort(frequencies_with_cent);
V_with_cent = V_with_cent(:, idx_with_cent);
%% 8. 结果输出与对比
disp(' ');
disp('=== 模态分析结果对比 ===');
disp('(考虑与不考虑3800RPM离心力作用)');
fprintf('旋翼总直径: %.5f m\n', rotor_diameter);
fprintf('单桨叶实际长度: %.5f m\n', L_blade);
disp(' ');
% 创建对比表格
disp('阶数 | 不考虑离心力 (Hz) | 考虑离心力 (Hz) | 变化量 (Hz) | 变化率 (%)');
disp('---------------------------------------------------------------------');
for i = 1:min(10, length(frequencies_no_cent))
freq1 = frequencies_no_cent(i);
freq2 = frequencies_with_cent(i);
delta = freq2 - freq1;
delta_percent = (delta / freq1) * 100;
fprintf('%4d | %19.4f | %18.4f | %11.4f | %10.2f\n', ...
i, freq1, freq2, delta, delta_percent);
end
%% 9. 振型可视化(对比两种情况)
dof_indices = [2, 3, 4, 5, 6]; % 2=剪切1, 3=剪切2, 4=扭转, 5=弯曲1, 6=弯曲2
dof_name = {'剪切方向1', '剪切方向2', '扭转方向', '弯曲方向1', '弯曲方向2'};
for fig = 1:length(dof_indices)
% 创建对比图
figure('Name', ['振型对比 - ', dof_name{fig}], 'Position', [100, 100, 1000, 800]);
for mode = 1:3 % 前3阶模态
% 不考虑离心力的振型
subplot(3, 2, 2*mode-1);
% 提取左桨叶(桨叶1)振型
disp_left = zeros(1, n_node_per_blade);
pos_left = zeros(1, n_node_per_blade);
for node_idx = 1:n_node_per_blade
actual_node = blade_nodes{1}(node_idx);
dof_idx = (actual_node - 1)*DOF_per_node + dof_indices(fig);
pos_left(node_idx) = node_pos{1}(node_idx);
if ismember(dof_idx, free_dof)
red_idx = find(free_dof == dof_idx);
disp_left(node_idx) = V_no_cent(red_idx, mode);
else
disp_left(node_idx) = 0;
end
end
% 提取右桨叶(桨叶2)振型
disp_right = zeros(1, n_node_per_blade);
pos_right = zeros(1, n_node_per_blade);
for node_idx = 1:n_node_per_blade
actual_node = blade_nodes{2}(node_idx);
dof_idx = (actual_node - 1)*DOF_per_node + dof_indices(fig);
pos_right(node_idx) = node_pos{2}(node_idx);
if ismember(dof_idx, free_dof)
red_idx = find(free_dof == dof_idx);
disp_right(node_idx) = V_no_cent(red_idx, mode);
else
disp_right(node_idx) = 0;
end
end
% 归一化
max_amp = max([abs(disp_left), abs(disp_right)]);
if max_amp > 0
disp_left = disp_left / max_amp;
disp_right = disp_right / max_amp;
end
% 绘图
rectangle('Position', [-hub_radius, -1.2, hub_diameter, 2.4], ...
'FaceColor', [0.8 0.8 0.8], 'EdgeColor', 'k', 'LineWidth', 1);
hold on;
plot(pos_left, disp_left, 'o-', 'LineWidth', 1.5, 'Color', [0.2 0.5 0.8], 'DisplayName', '左桨叶');
plot(-pos_right, disp_right, 's-', 'LineWidth', 1.5, 'Color', [0.8 0.3 0.3], 'DisplayName', '右桨叶');
plot(0, 0, 'ko', 'MarkerSize', 8, 'MarkerFaceColor', 'k');
xlabel('距中心距离(m)');
ylabel('归一化振幅');
title(['第', num2str(mode), '阶 (无离心力), f=', num2str(frequencies_no_cent(mode), '%.4f'), 'Hz']);
legend('Location', 'best');
grid on;
x_range = max(hub_radius + L_blade) * 1.1;
xlim([-x_range, x_range]);
ylim([-1.2, 1.2]);
plot([hub_radius, hub_radius], [-1.2, 1.2], 'k--', 'LineWidth', 0.5);
plot([-hub_radius, -hub_radius], [-1.2, 1.2], 'k--', 'LineWidth', 0.5);
box on;
hold off;
% 考虑离心力的振型
subplot(3, 2, 2*mode);
% 提取左桨叶(桨叶1)振型
disp_left_cent = zeros(1, n_node_per_blade);
for node_idx = 1:n_node_per_blade
actual_node = blade_nodes{1}(node_idx);
dof_idx = (actual_node - 1)*DOF_per_node + dof_indices(fig);
if ismember(dof_idx, free_dof)
red_idx = find(free_dof == dof_idx);
disp_left_cent(node_idx) = V_with_cent(red_idx, mode);
else
disp_left_cent(node_idx) = 0;
end
end
% 提取右桨叶(桨叶2)振型
disp_right_cent = zeros(1, n_node_per_blade);
for node_idx = 1:n_node_per_blade
actual_node = blade_nodes{2}(node_idx);
dof_idx = (actual_node - 1)*DOF_per_node + dof_indices(fig);
if ismember(dof_idx, free_dof)
red_idx = find(free_dof == dof_idx);
disp_right_cent(node_idx) = V_with_cent(red_idx, mode);
else
disp_right_cent(node_idx) = 0;
end
end
% 归一化
max_amp_cent = max([abs(disp_left_cent), abs(disp_right_cent)]);
if max_amp_cent > 0
disp_left_cent = disp_left_cent / max_amp_cent;
disp_right_cent = disp_right_cent / max_amp_cent;
end
% 绘图
rectangle('Position', [-hub_radius, -1.2, hub_diameter, 2.4], ...
'FaceColor', [0.8 0.8 0.8], 'EdgeColor', 'k', 'LineWidth', 1);
hold on;
plot(pos_left, disp_left_cent, 'o-', 'LineWidth', 1.5, 'Color', [0.2 0.5 0.8], 'DisplayName', '左桨叶');
plot(-pos_right, disp_right_cent, 's-', 'LineWidth', 1.5, 'Color', [0.8 0.3 0.3], 'DisplayName', '右桨叶');
plot(0, 0, 'ko', 'MarkerSize', 8, 'MarkerFaceColor', 'k');
xlabel('距中心距离(m)');
ylabel('归一化振幅');
title(['第', num2str(mode), '阶 (有离心力), f=', num2str(frequencies_with_cent(mode), '%.4f'), 'Hz']);
legend('Location', 'best');
grid on;
x_range = max(hub_radius + L_blade) * 1.1;
xlim([-x_range, x_range]);
ylim([-1.2, 1.2]);
plot([hub_radius, hub_radius], [-1.2, 1.2], 'k--', 'LineWidth', 0.5);
plot([-hub_radius, -hub_radius], [-1.2, 1.2], 'k--', 'LineWidth', 0.5);
box on;
hold off;
end
sgtitle(['双桨叶振型对比 - ', dof_name{fig}, ' (3800 RPM离心力影响)']);
end
帮我分析一下这个代码
最新发布