MATLAB脚本实现了一个转子系统的参数扫描和分岔分析

% 参数扫描范围
clc;
clear;
close all;S_values = 500:200:20000; % 转速范围% 定义系统参数
N = 5; % 质量点数量
num_nodes = N; % 节点数
num_dofs_per_node = 4; % 每个节点的自由度数
num_elements = num_nodes-1; % 单元数
total_dofs = num_nodes * num_dofs_per_node; % 总自由度数% 示例参数
m = [2.3182 2.5249 3.0181 4.9458 15.6646]*1; % 质量 (kg)Jd = [5.5671e-4 1.3564e-4 1.1411e-3 6.3456e-4 4.8254e-2]*1; % 转动惯量 (kg*m^2)
Jp = [2.7835e-4 6.7821e-5 5.7059e-4 3.1728e-4 2.4127e-2]*1; % 极转动惯量 (kg*m^2)e = [50 50]*1e-6; % 偏心距 (m) 5
% S=10000;% r/min% phi = [0, pi/4, pi/2, 3*pi/4, pi]; % 相位角 (rad)
l = [0.246 0.1045 0.118 0.1795]*1; % 长度 (m)
nu = [0.3, 0.3, 0.3, 0.3, 0.3]; % 泊松比% f=S/60;
% omega=2*pi*f;
vi=0.3;
% 定义材料和截面属性（示例值）
% E = 2e11; % 弹性模量
% R0=40e-3;
% I=pi*R0^4/4; %直径惯量
% rho = 7860; % 材料密度
%  质量矩阵、陀螺矩阵、刚度矩阵
% -------------质量矩阵
z0=zeros(total_dofs/2, total_dofs/2);
mj=[];
J1=[];
for i=1:Nmj=[mj, m(i) Jd(i)];J1=[J1,0,Jp(i)];% 惯性矩阵
end
Mx = diag(mj);
My = Mx;
M=[Mx z0;z0 My];
M_global=M;
% 陀螺矩阵
J2=diag(J1);% ---转子刚度矩阵----------
% 参数初始化
K =z0; % 初始化刚度矩阵% 计算b参数
b1 = zeros(N-1, 1);
b2 = zeros(N-1, 1);
b3 = zeros(N-1, 1);EI=[2.4627e5 4.5537e5 5.1703e5 1.3473e7];
EA=[1.2070e9 1.9007e9 2.0866e9 1.584e10];for i = 1:N-1b1(i) = 12 * EI(i) / (l(i)^3);b2(i) = 0.5 * l(i) * b1(i);b3(i) = 1/6 * l(i)^2 * b1(i);
end% 填充刚度矩阵
for i = 1:Nif i == 1K(1, 1) = b1(1);K(1, 2) = b2(1);K(1, 3) = -b1(1);K(1, 4) = b2(1);K(2, 2) = l(1)*b2(1)-b3(1);K(2, 3) = -b2(1);K(2, 4) = b3(1);elseif i<=N-1K(2*i-1, 2*i-1) = b1(i-1) + b1(i);K(2*i-1, 2*i) = -b2(i-1) + b2(i);K(2*i-1, 2*i+1) = -b1(i);K(2*i-1, 2*i+2) = b2(i);K(2*i, 2*i) = sum(l(i-1:i)*b2(i-1:i)) - sum(b3(i-1:i));K(2*i, 2*i+1) = -b2(i);K(2*i, 2*i+2) = b3(i);elseK(2*i-1, 2*i-1) = b1(i-1);K(2*i-1, 2*i) = b2(i-1);K(2*i, 2*i) = l(i-1)*b2(i-1)-b3(i-1);end
end% 对称填充
K= K + K';% 转子总刚度矩阵
Kr=[K z0;z0 K];% 系统总刚度矩阵K_global=Kr;% 计算阻尼矩阵（假设为比例阻尼）
aa=298;
bb=3.3482e-6;S = S_values(1);
f = S / 60;
omega = 2 * pi * f;
T=60/S;
dt = T/256; % 时间步长 1e-5
t0 = 0:dt:200*T; % 时间向量
% 初始化位移、速度、加速度向量
u = zeros(total_dofs, length(t0))+1e-6;
v = zeros(total_dofs, length(t0))+1e-6;
a = zeros(total_dofs, length(t0))+1e-6;
F = zeros(total_dofs, length(t0))+1e-6;
Fe=F;%  %Newmark-β方法参数
beta=0.25;                  % β
delta=0.5;                  % δ
a0=1/(beta*dt^2);           % 1/(βΔt^2)
a1=delta/(beta*dt);         % δ/(βΔt)
a2=1/(beta*dt);             % 1/(βΔt)
a3=1/(2*beta)-1;            % 1/(2β)-1
a4=delta/beta-1;            % δ/β-1
a5=dt/2*(delta/beta-2);     % Δt/2*(δ/β-2)
a6=dt*(1-delta);            % Δt(1-δ)
a7=dt*delta;                % Δtδbifurcation_data1 = [];bifurcation_data3 = [];bifurcation_data5 = [];
bifurcation_data7 = [];bifurcation_data9 = [];bifurcation_data = cell(5,1);
poincare_data = cell(5,1);  % 新增初始化
for k = 1:5bifurcation_data{k} = [];poincare_data{k} = [];  % 每个质量点单独存储庞加莱数据
end
bifurcation_data0 = cell(5,1);% 对每个转速值进行模拟
tic
for idx = 1:length(S_values)S = S_values(idx);f = S / 60;omega = 2 * pi * f;T=60/S;dt = T/256; % 时间步长 1e-5
t0 = 0:dt:200*T; % 时间向量% 动力学计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 定义时间参数t_transient = 20*T;steady_index = find(t0 >= t_transient, 1);G=omega*[z0 J2;-J2' z0;];C_global = aa* M_global+bb*K_global-G;% %偏心% %pianxin=[9,19];me=20;% pian=[20 20.0]*1e-6;fai1=0;fai2=0;fai3=0;%开始计算num=0;% bearing=[1,9,7,15];% bearing=[3,7,13,17];bearing=[3,4,7,8,13,14,17,18];% zhongli=[9,11,13,15];zhongli=2*N+1:2:4*N-1;Nb=15;Kb=13.34e9;R=63.9e-3;r=40.1e-3;c=0*1e-6;%轴承参数g=10;anglea=0:2*pi/Nb:(2*pi-2*pi/Nb);%参数准备anglea=anglea';xishu=[65/(120+65) 75/(130+75)];% xishu=r/(R+r);jilu=zeros(length(m),length(t0));jilu_KL=[];al=[];Fb1=zeros(length(bearing),length(t0));% 将最后的结果作为下一次积分的初值u(:,1) = u(:,end);v(:,1) = v(:,end);a(:,1) = a(:,end);F = zeros(total_dofs, length(t0))+1e-10;Fe=F;% 迭代计算for i = 2:length(t0)t=(i-1)*dt;% % 保持架的转速omega_cage = omega * xishu(1);% [Fx1,Fy1,Fz1,Mx1,My1]=zhouchengli(Nb,omega_cage,t);[Fx1,Fy1,Fz1,Mx1,My1]=zhouchengli2(u(3,i-1),u(13,i-1),0,u(14,i-1),u(4,i-1),Nb,omega_cage,t);% % 保持架的转速omega_cage = omega * xishu(2);% [Fx2,Fy2,Fz2,Mx2,My2]=zhouchengli(Nb,omega_cage,t);[Fx2,Fy2,Fz2,Mx2,My2]=zhouchengli2(u(7,i-1),u(17,i-1),0,u(18,i-1),u(8,i-1),Nb,omega_cage,t);Fb1(:,i-1)=[Fx1;Mx1;Fx2;Mx2;Fy1;My1;Fy2;My2];% 轴承力F(bearing,i-1)=F(bearing,i-1)-2*[Fx1,Mx1,Fx2,Mx2,Fy1,My1,Fy2,My2]';% 重力F(zhongli,i-1)=F(zhongli,i-1)-m(1:N)'*g;% 偏心力F(pianxin,i-1)=F(pianxin,i-1)+[me*e(1)*omega^2*cos(omega*t+fai1);me*e(1)*omega^2*sin(omega*t+fai1);];% m(2)*e(1)*omega^2*cos(omega*t+fai1);% m(3)*e(2)*omega^2*cos(omega*t+fai1);Ke=K_global+a0*M_global+a1*C_global;         % 等效刚度矩阵Fe(:,i-1)=F(:,(i-1))+M_global*(a0*u(:,i-1)+a2*v(:,i-1)+a3*a(:,i-1))+C_global*(a1*u(:,i-1)+a4*v(:,i-1)+a5*a(:,i-1));   % 求解t+Δt时刻的等效载荷u(:,i)=Ke\Fe(:,i-1);                                                                      % 求解t+Δt时刻的位移a(:,i)=a0*(u(:,i)-u(:,i-1))-a2*v(:,i-1)-a3*a(:,i-1);                                    % 求解t+Δt时刻的加速度v(:,i)=v(:,i-1)+a6*a(:,i-1)+a7*a(:,i);                                                  % 求解t+Δt时刻的速度endend% 动力学计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%t_transient = 20*T;% 提取稳态阶段数据（t >= t_transient）steady_indices = find(t0 >= t_transient);if isempty(steady_indices)continue;end% x_steady = y(steady_indices, 1);for j3=4x_steady = u(2*j3-1,steady_indices:end);% 寻找位移的局部极大值% [pks, ~] = findpeaks(x_steady(j3,:));for ff=1:20pks(ff)=max(x_steady(1+round(ff*T/dt):round((ff+1)*T/dt)));end% 保存分岔数据if ~isempty(pks)bifurcation_data0{j3} = [bifurcation_data0{j3};omega*ones(length(pks), 1), pks(:)];endendelapsed_time = toc; % 已耗时fprintf('已完成 %.2f %%，耗时 %.1f 秒 \n', idx * 100 / length(S_values), elapsed_time);
end% 绘制分岔图
j2=4;
figure;
scatter(bifurcation_data0{j2}(:,1)*60/2/pi, bifurcation_data0{j2}(:,2));
xlabel('\omega');
ylabel('Peak Displacement');
title('Rotor System Bifurcation Diagram');
grid on;figure;
plot(bifurcation_data0{j2}(:,1)*60/2/pi, bifurcation_data0{j2}(:,2),'k.');
xlabel('\omega');
ylabel('Peak Displacement');
title('Rotor System Bifurcation Diagram');
grid on;figure
plot(u(2*j2-1),v(2*j2-1),'b');
title([' 集中质量', num2str(j2), '相图 ']);

该MATLAB脚本实现了一个转子系统的参数扫描和分岔分析，具体功能如下：

---

1. 转子系统建模

• 系统参数定义：
  • 质量点数量 N = 5，每个节点自由度数为4（横向位移、转角等）。
  • 定义了质量 m、转动惯量 Jd 和 Jp、偏心距 e、单元长度 l、弹性模量 EI 和 EA 等参数。
• 矩阵构建：
  • 生成质量矩阵 M、刚度矩阵 K 和陀螺矩阵 G，并通过比例阻尼模型构建阻尼矩阵 C_global。

---

2. 参数扫描分析

• 转速范围扫描：
  • 对转速 S_values = 500:200:20000 rpm 进行循环，模拟不同转速下的系统响应。
• 动力学计算：
  • 使用 Newmark-β方法 求解转子系统的瞬态动力学响应，考虑以下激励：
    • 轴承力：通过函数 zhouchengli2 计算非线性轴承力（需自定义）。
    • 重力：作用在垂直方向。
    • 偏心力：由质量偏心引起的周期性激励（频率与转速相关）。

---

3. 分岔分析

• 稳态数据提取：
  • 舍弃前20个周期的瞬态响应，提取稳态阶段的位移数据。
  • 对每个转速下的稳态位移，记录局部极大值（峰值）。
• 分岔图绘制：
  • 绘制第4个质量点的分岔图，横轴为转速（ω），纵轴为位移峰值，展示系统随转速变化的非线性行为（如周期、拟周期或混沌运动）。

---

4. 结果可视化

• 分岔图：
  • 使用散点图（scatter）和点线图（plot）展示分岔数据。
• 相图：
  • 绘制第4个质量点的位移-速度相图，用于观察系统的动力学状态（如极限环、混沌吸引子等）。

---

关键功能总结

1. 参数化建模：通过矩阵组装实现多自由度转子系统的动力学建模。
2. 非线性响应分析：结合轴承力、偏心力等非线性因素，模拟转子在不同转速下的振动特性。
3. 分岔检测：通过峰值提取和绘图，识别系统的临界转速和非线性现象（如倍周期分岔）。

---

注意事项

• 需补充自定义函数 zhouchengli2 的实现（轴承力计算）。
• 分岔分析仅针对第4个质量点，可通过修改 j2 或循环分析其他节点。
• 计算耗时较长（因参数扫描范围大），可通过并行计算优化。

该脚本适用于转子动力学研究，特别是临界转速和非线性振动行为的分析。