简介本资源是一套面向计算机、电子信息工程及数学等专业本科生的Matlab动力学与振动分析实践代码集专为课程设计、期末大作业及毕业设计场景打造帮助学习者快速掌握系统建模、运动方程求解与频响分析等核心能力。压缩包共36个文件含32个功能完备的.m主程序文件实现参数化建模、数值仿真与结果可视化、3个.t模板文件用于数据格式适配与接口扩展及1张说明性PNG图整体仅156KB轻量易用。已有76人下载学习代码兼容Matlab 2014/2019a/2024a多版本所有脚本均采用清晰模块化结构关键变量与算法步骤均配有中文注释附带可直接运行的案例数据省去数据准备环节显著提升仿真实验效率。1. 用 MATLAB 做动力学与振动分析不是调用几个函数就完事——它要你真正理解系统建模、数值求解与物理响应之间的闭环关系很多刚接触《振动力学》或《机械系统动力学》课程的学生拿到“matlab代码.rar”压缩包后第一反应是解压、运行、看图——结果报错、曲线不对、相位反了、频谱毛刺一堆。问题不在代码本身而在于没把“单自由度阻尼振动方程”和ode45的输入格式对齐没意识到fft默认零频在首项而工程频谱要求中心对齐更没注意bode绘图时采样率不足会引发混叠失真。这套工作流面向的是需要复现教材案例如倪振华《振动力学》第3章受迫振动、验证实验数据电机振动信号频谱特征提取、或搭建仿真原型机器人关节柔性动力学的工程师与高年级本科生。它不依赖 Simulink 图形界面全部基于脚本化建模与可追溯计算核心能力是从微分方程出发生成时域响应转换为频域特征再通过参数敏感性分析反推结构刚度/阻尼。下面我们就从最基础的单自由度系统开始一层层拆解真实项目中必须跨过的三道坎建模规范、求解器配置、结果可信度验证。2. 把物理方程写成 ode45 能解的形式状态变量定义、雅可比矩阵显式化与初始条件物理意义校验2.1 单自由度有阻尼受迫振动的标准建模流程动力学建模的第一步不是敲代码而是明确系统自由度、建立牛顿第二定律或拉格朗日方程。以质量-弹簧-阻尼器串联系统为例其运动微分方程为$$ m\ddot{x} c\dot{x} kx F_0 \cos(\omega t) $$该二阶常微分方程不能直接传给ode45必须降阶为一阶方程组。标准做法是定义状态变量$ x_1 x $位移$ x_2 \dot{x} $速度则导数关系为$ \dot{x}_1 x_2 $$ \dot{x}_2 \frac{1}{m} \left[ -c x_2 - k x_1 F_0 \cos(\omega t) \right] $这个转换不是数学技巧而是物理约束ode45求解器内部采用自适应步长若状态变量物理量纲混乱如用加速度而非速度作状态会导致误差估计失效步长剧烈震荡甚至发散。2.2 编写可调试的 odefun 函数避免隐式依赖与全局变量陷阱常见错误是把参数m,c,k写成全局变量或硬编码在函数内。正确做法是用匿名函数闭包传递参数保证函数纯度与可复用性% 参数定义实际项目中应从 config.mat 或 JSON 加载 m 2.5; % kg c 8.0; % N·s/m k 120; % N/m F0 15; % N omega 6.0; % rad/s % 构造状态方程函数句柄 —— 注意t 必须是第一个输入x 是第二个 odefun (t, x) [x(2); ... (1/m) * (-c*x(2) - k*x(1) F0*cos(omega*t))]; % 初始条件x(0)0.02 m, v(0)0 m/s → 物理意义明确 x0 [0.02; 0]; % 时间跨度需覆盖至少 5 个激励周期且满足奈奎斯特采样 tspan [0, 5*(2*pi/omega)]; % 5 个完整周期 % 调用 ode45 —— 不指定相对误差容限时默认为 1e-3对振动问题常不够 [t, x] ode45(odefun, tspan, x0, odeset(RelTol, 1e-6, AbsTol, 1e-9));提示odeset中RelTol和AbsTol必须同时设置。仅调RelTol会导致小位移阶段绝对误差超标仅调AbsTol会在大振幅阶段相对误差失控。振动问题典型容限组合是RelTol1e-6AbsTol1e-9对应毫米级位移下亚微米精度。2.3 雅可比矩阵显式提供加速求解并提升刚性系统稳定性当系统阻尼极小如 $ c 0.1\sqrt{km} $或存在高频模态耦合时ODE 系统呈现刚性特征。此时ode45可能步长过小、耗时剧增。显式提供雅可比矩阵可显著改善性能% 雅可比矩阵 J d(odefun)/d(x)按列排列J(:,1)∂f/∂x1, J(:,2)∂f/∂x2 jacfun (t,x) [0, 1; ... -k/m, -c/m]; % 将雅可比嵌入选项 opts odeset(RelTol,1e-6, AbsTol,1e-9, Jacobian, jacfun); [t, x] ode45(odefun, tspan, x0, opts);雅可比矩阵在此例中为常数矩阵但若系统含非线性弹簧如 $ f_s kx \alpha x^3 $则雅可比需实时计算此时Jacobian应设为函数句柄而非数值矩阵。3. 从时域响应到频域特征FFT 参数设置、窗函数选择与功率谱密度物理标定3.1 FFT 前必须做的三件事去直流、补零、选窗直接对x(:,1)位移序列做fft得到的频谱必然失真。正确预处理流程如下% 提取位移信号稳态段剔除初始瞬态 N_transient round(0.2*length(t)); % 剔除前 20% 瞬态响应 x_disp x(N_transient:end, 1); % 1. 去直流分量 —— 否则零频幅值淹没有效频谱 x_disp x_disp - mean(x_disp); % 2. 选择汉宁窗Hanning抑制频谱泄漏长度与信号一致 win hanning(length(x_disp)); x_win x_disp .* win; % 3. 补零至 2 的整数次幂加速 FFT但不增加频率分辨率 N_fft 2^nextpow2(length(x_win)); X_fft fft(x_win, N_fft); % 计算单边幅值谱物理意义各频率分量的位移幅值 amp_spec (2/N_fft) * abs(X_fft(1:N_fft/21)); freq_vec (0:N_fft/2) * (1/(t(2)-t(1))) / N_fft; % Hz 单位注意fft输出是双边谱abs()后需乘 2 并取前半除 DC 和 Nyquist 点外才能得到真实幅值。t(2)-t(1)是实际采样间隔绝不能用t(end)/length(t)近似——因ode45输出时间点非均匀。3.2 功率谱密度PSD的工程标定从pwelch到 g²/Hz 单位转换实验室振动台输出常以加速度单位m/s²给出而 PSD 标准单位是 (m/s²)²/Hz。若原始信号是位移m需先微分得加速度% 对位移信号二次微分得加速度使用差分法注意边界处理 dt mean(diff(t)); % 平均采样间隔 acc gradient(gradient(x_disp, dt), dt); % 使用 pwelch 计算 PSD —— 自动分段、加窗、平均抗噪能力强 [pxx, f_pxx] pwelch(acc, hanning(1024), 512, 1024, 1/dt); % 单位转换若 acc 单位为 m/s²则 pxx 单位为 (m/s²)²/Hz % 若需转换为 g²/Hz1 g 9.80665 m/s²则 pxx_g2 pxx / (9.80665^2);pwelch的三个关键参数窗长1024、重叠点数512、FFT 点数1024决定了频谱平滑度与频率分辨率的平衡。窗长越长频率分辨率越高Δf fs/N_window但时间局部性越差重叠越多平均次数越多方差越小。3.3 共振峰识别与阻尼比提取半功率带宽法的 MATLAB 实现共振频率处的 PSD 峰值对应系统固有频率其宽度反映阻尼大小。半功率带宽法3dB 带宽是工程常用方法% 找到 PSD 主峰排除 DC 附近低频干扰 [~, idx_peak] max(pxx_g2(5:end)); % 跳过前 4 点0~4 Hz 常为噪声 idx_peak idx_peak 4; f_res f_pxx(idx_peak); % 计算半功率点峰值功率的一半 p_half pxx_g2(idx_peak) / 2; % 向左找第一个低于 p_half 的点 idx_left find(pxx_g2(1:idx_peak) p_half, 1, last); % 向右找第一个低于 p_half 的点 idx_right find(pxx_g2(idx_peak:end) p_half, 1, first) idx_peak - 1; % 3dB 带宽 Δf_3dB f_right - f_left delta_f f_pxx(idx_right) - f_pxx(idx_left); % 阻尼比 ζ Δf_3dB / (2 * f_res) zeta_est delta_f / (2 * f_res); fprintf(估算阻尼比 ζ %.4f\n, zeta_est);该方法要求 PSD 峰形对称且信噪比 20 dB。若峰形畸变应改用拟合 SDOF 系统频率响应函数FRF的方法。4. 多自由度系统建模实战从质量-刚度矩阵组装到模态叠加法验证4.1 用物理参数直接构建 M、C、K 矩阵避免手算耦合项错误两自由度系统如车辆悬架模型需显式写出质量、阻尼、刚度矩阵。以车体质量 $ m_1 $、轮胎质量 $ m_2 $、悬架刚度 $ k_1 $、轮胎刚度 $ k_2 $ 为例% 物理参数 m1 1200; k1 25000; c1 1200; % 车体-悬架 m2 50; k2 200000; c2 0; % 轮胎 % 组装全局质量矩阵 M对角阵 M diag([m1, m2]); % 刚度矩阵 KK11k1k2, K12K21-k2, K22k2 K [k1k2, -k2; ... -k2, k2]; % 阻尼矩阵 C同理 C [c1c2, -c2; ... -c2, c2]; % 外部激励路面不平度作为基础激励转化为等效作用力 % 假设路面位移 u(t) A*sin(ωt)则等效力向量 F_ext [k2*u; k2*u]矩阵组装必须符合力学约定对角线元素为自作用项非对角线为耦合作用项符号由相对位移方向决定。手算易错建议用符号计算工具Symbolic Math Toolbox辅助验证。4.2 将 MCK 系统降阶为状态空间确保维度匹配与物理一致性多自由度系统需将 $ M\ddot{x} C\dot{x} Kx F $ 降阶为 $ \dot{z} A z B u $ 形式。标准变换为$ z [x; \dot{x}] $则 $ \dot{z} [\dot{x}; \ddot{x}] $$ \ddot{x} M^{-1}(F - C\dot{x} - Kx) $MATLAB 实现n size(M,1); A [zeros(n), eye(n); ... -M\K, -M\C]; B [zeros(n); M\eye(n)]; C_out [eye(n), zeros(n)]; % 输出位移 D_out zeros(n); % 构建状态空间模型用于后续 bode、step 分析 sys ss(A, B, C_out, D_out); % 验证计算模态频率eig(K,M) 应与 bode 峰值一致 [V,D] eig(K,M); omega_n sqrt(diag(D)); % rad/s f_n omega_n/(2*pi); % Hz disp(理论固有频率 (Hz):); disp(f_n);提示eig(K,M)返回广义特征值其平方根即固有圆频率。若f_n与bode(sys)中的谐振峰频率偏差 1%说明矩阵组装有误如刚度符号、质量位置错位。4.3 模态叠加法验证用前 2 阶模态重构时域响应模态叠加法是验证数值解精度的黄金标准。对无阻尼系统响应可表示为$$ x(t) \sum_{i1}^{r} q_i(t) \phi_i $$其中 $ \phi_i $ 为第 i 阶模态向量$ q_i(t) $ 为广义坐标。MATLAB 实现% 提取前 2 阶模态已归一化 Phi V(:,1:2); % 2×2 模态矩阵 % 计算模态质量、刚度矩阵 M_phi Phi * M * Phi; % 对角阵 K_phi Phi * K * Phi; % 对角阵 % 每阶模态独立求解q_i ω_i² q_i φ_i^T F(t) / m_i q0 Phi * x0(1:2); % 初始广义位移 dq0 Phi * x0(3:4); % 初始广义速度 % 对每阶构造 ode 函数并求解 q_sol zeros(length(t), 2); for i 1:2 omega_i sqrt(K_phi(i,i)/M_phi(i,i)); f_i (t) Phi(:,i) * F_func(t) / M_phi(i,i); % 广义力 odefun_q (t,q) [q(2); -omega_i^2*q(1) f_i(t)]; [~, q_i] ode45(odefun_q, t, [q0(i); dq0(i)]); q_sol(:,i) q_i(:,1); end % 重构物理位移 x_modal q_sol * Phi;将x_modal与ode45直接求解结果对比若 RMS 误差 0.5%说明数值解可靠否则需检查ode45容限或矩阵组装。5. 工程级振动分析技巧从电机振动信号数据集加载到故障特征频率标记5.1 加载公开振动数据集如 CWRU 轴承数据并重采样对齐许多用户搜索“声音振动信号电机数据集”实际指 Case Western Reserve UniversityCWRU轴承故障数据。其原始采样率为 12 kHz但常需重采样以匹配模型采样率% 下载并加载 CWRU 数据假设已存为 .mat load(12kDriveEnd_B014.mat); % 内含 signal 变量 fs_orig 12000; % 原始采样率 fs_target 5000; % 目标采样率需 2×最高关注频率 % 重采样使用 resample 避免混叠自动设计抗混叠滤波器 signal_rs resample(signal, fs_target, fs_orig); % 时间向量 t_data (0:length(signal_rs)-1) / fs_target; % 提取稳态段跳过启停瞬态 N_start round(0.1*length(t_data)); N_end round(0.9*length(t_data)); signal_steady signal_rs(N_start:N_end); t_steady t_data(N_start:N_end);resample内置抗混叠滤波器比decimate更适合保留故障特征频率如轴承外圈故障频率 BPFO。若需更高保真度可用designMultirateFIR手动设计 FIR 滤波器。5.2 标记电机故障特征频率BPFO、BPFI、BSF、FTF 的 MATLAB 计算轴承故障频率取决于几何参数与转速。CWRU 数据标注了驱动端转速 RPM据此计算RPM 1772; % 示例转速实际从文件名或元数据读取 n RPM/60; % 转速Hz % CWRU 轴承参数6205-2RS 轴承 d 0.0159; % 滚子直径 (m) D 0.072; % 节径 (m) N 9; % 滚子数 alpha 0; % 接触角深沟球轴承≈0 % 计算故障频率单位Hz BPFO N*n/2 * (1 - d/D*cos(alpha)); % 外圈故障 BPFI N*n/2 * (1 d/D*cos(alpha)); % 内圈故障 BSF D*n/(2*d) * (1 - (d/D*cos(alpha))^2); % 滚子故障 FTF n/2 * (1 - d/D*cos(alpha)); % 保持架故障 fprintf(BPFO%.1f Hz, BPFI%.1f Hz, BSF%.1f Hz, FTF%.1f Hz\n, ... BPFO, BPFI, BSF, FTF);这些频率应在 PSD 图上用垂直线标记便于人工判读。例如figure; plot(f_pxx, 10*log10(pxx_g2)); hold on; xline(BPFO, --r, BPFO); xline(BPFI, --g, BPFI); xlabel(Frequency (Hz)); ylabel(PSD (g^2/Hz)); legend(PSD, BPFO, BPFI);5.3 振动信号包络谱分析提取早期微弱冲击特征轴承早期故障在时域呈微弱冲击在频域被基频和谐波淹没。包络谱Envelope Spectrum可增强冲击特征% 1. 带通滤波聚焦于共振频带如 3–5 kHz [b,a] butter(4, [3000,5000]/(fs_target/2), bandpass); signal_bp filtfilt(b,a,signal_steady); % 2. 解调取绝对值 低通滤波截止频率 1 kHz signal_abs abs(signal_bp); [b_lp,a_lp] butter(4, 1000/(fs_target/2), low); envelope filtfilt(b_lp,a_lp,signal_abs); % 3. 对包络信号做 FFT N_env 2^nextpow2(length(envelope)); env_fft fft(envelope, N_env); env_amp (2/N_env) * abs(env_fft(1:N_env/21)); f_env (0:N_env/2) * (fs_target/N_env); % 4. 在包络谱中标记故障频率此时单位为 Hz对应冲击重复率 figure; plot(f_env, env_amp); xline(BPFO, --r); xline(BPFI, --g); xlabel(Envelope Frequency (Hz)); ylabel(Amplitude); title(Envelope Spectrum);包络谱峰值若出现在 BPFO 或其倍频即为外圈故障确证。此方法比原始 PSD 更敏感是工业预测性维护的核心技术。振动分析不是 MATLAB 函数的堆砌而是物理建模、数值方法、信号处理三者的严密闭环。每一次ode45的收敛、每一处pwelch的峰值、每一个xline标出的 BPFO都在验证你对系统本质的理解是否到位。当电机振动信号的包络谱清晰显示 BPFO 倍频族时那不是代码的胜利是你把课本公式、实验数据与工程直觉真正焊在了一起。本文还有配套的精品资源点击获取