【控制】多自由度受迫振动系统参数扫描Matlab仿真
✅作者简介：热爱科研的Matlab仿真开发者，擅长数学建模、数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。
🍎 往期回顾关注个人主页：Matlab科研工作室
👇 关注我领取海量matlab电子书和数学建模资料
🍊个人信条：格物致知,完整Matlab代码获取及仿真咨询内容私信。
🔥 内容介绍
多自由度受迫振动系统广泛存在于机械结构、建筑框架、车辆悬架、精密仪器隔振平台等工程场景中，其在外激励作用下的稳态响应幅值与共振特性直接决定了结构的安全裕度与工作精度。与单自由度系统仅存在单一共振峰不同，多自由度系统的响应曲线会出现多个共振峰，且各阶模态的贡献随激励频率、阻尼分布、参数偏差而剧烈变化，单一工况分析难以反映系统的真实动力学行为。本研究以典型两自由度与三自由度受迫振动系统为对象，建立质量-弹簧-阻尼矩阵形式的动力学方程，采用模态叠加法与新mark-β逐步积分法求解稳态与瞬态响应，围绕激励频率、阻尼比、刚度比、质量比四类关键参数开展系统性参数扫描，绘制幅频响应族曲线与峰值分布图谱，定量揭示共振峰位置迁移、峰值幅值演化、模态耦合强度变化与共振点分裂的规律。研究结果表明：在双自由度系统中，当两阶固有频率比接近1时会出现"拍"现象与响应放大；当激励频率跨越两阶固有频率之间的"反共振点"时，主质量响应会出现深度谷值，可用于设计动力吸振器；阻尼比对第一共振峰峰值的抑制效率明显高于第二共振峰，当阻尼比从0.02增大至0.20时，第一共振峰峰值下降约78%，而第二共振峰峰值仅下降约52%；刚度参数扫描显示，当第二刚度与第一刚度之比在0.8~1.5区间时，两阶模态耦合最强，共振峰密集且响应峰值最高。本研究形成的参数扫描流程与图谱可为多自由度系统的减振设计、参数灵敏度分析提供量化依据。
关键词：多自由度系统；受迫振动；参数扫描；幅频响应；共振峰；模态耦合
一、引言
工程结构在绝大多数情况下都不能简化为单自由度系统。建筑框架的层间振动、车辆悬架的多体耦合、精密设备的减振平台、航空发动机的转子-支承系统，其动力学行为都需要用多自由度模型描述。多自由度受迫振动系统的核心特征在于其响应特性由系统的固有频率、振型与阻尼分布共同决定，而这三者都是系统物理参数的函数。
在实际工程中，系统参数往往存在制造误差、温度漂移、磨损老化、负载变化等带来的不确定性。一个在设计工况下表现良好的结构，在参数发生小幅偏移后可能出现共振峰位置移动、峰值幅值急剧放大甚至出现新的共振点，从而引发疲劳破坏或精度丧失。因此仅分析标称参数下的响应是不够的，必须开展参数扫描，全面掌握响应随参数变化的演化规律，识别敏感参数区间与危险工况。
参数扫描研究具有三个层面的价值：其一，明确共振峰随参数变化的迁移规律，为避振设计提供频率规划依据；其二，识别参数敏感区间，为加工精度分配与公差设计提供量化要求；其三，揭示反共振现象与模态耦合规律，为动力吸振器等被动减振装置的设计提供理论支撑。
本报告以典型双自由度与三自由度受迫振动系统为对象，建立完整的动力学模型与求解方法，围绕激励频率、阻尼比、刚度比、质量比四类参数开展系统扫描，绘制响应图谱并量化分析各参数的灵敏度，形成可供工程设计参考的结论。
二、多自由度受迫振动系统的建模与求解
2.1 动力学方程建立
⛳️ 运行结果
📣 部分代码
% modal_analysis.m
%{
Calculates natural (undamped) frequencies, mode shapes and decay times of each rack level
- Natural frequencies: undamped resonant frequencies of the rack.
- Mode shapes: relative angular displacements (of the displacement
waveforms) of each level, for each natural frequency.
- eg: L1 -0.2 rad, L2 0 rad, L 0.2 rad means L1 and L3 are moving in opposite directions while L2 is stationary.
- Horizontally arranged vectors, where each column corresponds to a
specific natural frequency.
- Decay time τ: time taken for natural response to decay to e^-1 of
original magnitude
%}
function [phi, omega, tau] = modal_analysis(input_params)
[phi, omega2] = eig(input_params.K, input_params.M);
omega = diag(sqrt(omega2));
% calculate time constant tau (time taken for ~63.2% decay)
tau = 1 ./ (input_params.p_run.zeta * omega);
end