简介面向机械工程与齿轮传动系统动力学分析的一份 MATLAB 源码包主要解决直齿轮副非线性动力学建模与求解问题。资源以单自由度简化模型为基础基于牛顿第二定律建立包含质量、加速度、力与扭矩的运动方程并通过相图刻划位置与速度的状态轨迹帮助判断系统稳定性与周期性利用傅里叶变换分析振动信号识别由啮合冲击、制造误差或载荷波动引起的频谱特征服务于齿轮故障诊断与结构优化。包内为1个可运行的 .m 脚本压缩后仅 1KB体积小巧代码结构清晰便于读者快速修改参数并重新运行。已有 678 人学习适合需要入门齿轮非线性动力学或希望借助编程完成建模、求解与可视化验证的工程技术人员也可在此基础上扩展含间隙、时变啮合刚度等复杂因素进一步研究分岔与混沌现象。1. 单自由度动力学方程为什么值得用MATLAB跑一遍一个弹簧、一个质量、一个阻尼器单自由度动力学方程看起来比 Navier-Stokes 简单得多但真用手推解简谐激励的瞬态响应时初始位移和初始速度决定的那两项常数经常让人在积分结果上栽跟头。这也是为什么我始终建议先用 MATLAB 把 m x c x k x F(t) 完整跑一遍不是为了得到一条曲线而是为了验证你对这个系统的直觉——共振峰值并不是恰好出现在固有频率处而是在略偏左的地方。这类模型在机械振动、结构减震、动态测试和控制系统里都是被反复调用的基础单元新手可以用它掌握 ode45 和状态空间老手可以用它快速评估算法参数。2. 把单自由度动力学方程改写成MATLAB可执行的ode45状态空间在写任何求解脚本之前先把方程标准化这一步常被跳过。2.1 运动方程标准化先除以质量还是保留原参数原始方程是m x(t) c x(t) k x(t) F(t)我一般习惯先除以 m写成x 2 ζ ω_n x ω_n² x f(t)其中 ω_n sqrt(k/m)ζ c / (2 sqrt(m k))。这样做的直接好处是二阶项系数为 1ode45 内部的相对误差控制对各项量级更敏感如果 m 是 1e-3 而 k 是 1e6直接用原始系数写回调数值梯度会跨好几个量级容易触发不必要的步长加密甚至让结果看起来像高频振荡。更重要的是ζ 和 ω_n 这两个参数直接决定响应的形态看到 ζ 小于 1就能预期响应是衰减振荡而不是单调趋近这对排查求解结果是很有用的。也要注意除以 m 只适合单自由度方程本身。如果后面扩展到多自由度系统质量矩阵存在耦合不能简单逐行除而应使用状态空间下的质量矩阵求逆。单自由度时不需要计较这一步但思路要清楚。2.2 状态空间两个一阶方程替换一个二阶方程ode45 只能处理一阶常微分方程组而动力学方程是二阶的所以要做降阶。令 y1 xy2 x则由 y1 y2 和原方程可得y1 y2 y2 (F(t) - c y2 - k y1) / m这是一个标准的二维状态空间。外激励 F(t) 用函数句柄表示这样自由振动、简谐激励、随机载荷都不需要改导数函数只换句柄即可。如果 F(t) 是线性系统的外部输入这个写法也方便以后改成时变参数比如 c(t) 或 k(t)。2.3 最小可运行代码sdof函数加ode45调用把状态空间写成 MATLAB 函数建议单独保存为 sdofODE.m方便多个脚本共用function dydt sdofODE(t, y, m, c, k, Ffun) % 单自由度动力学方程的状态空间形式 % y(1) 位移, y(2) 速度 dydt zeros(2,1); dydt(1) y(2); dydt(2) (Ffun(t) - c*y(2) - k*y(1)) / m; end然后调用 ode45 并绘图m 1.2; % 质量 kg c 0.4; % 阻尼系数 N*s/m k 30; % 刚度 N/m Ffun (t) 0; % 自由振动 tspan [0 10]; y0 [0.02; 0]; opts odeset(RelTol, 1e-8, AbsTol, 1e-10); [t, y] ode45((t,y) sdofODE(t, y, m, c, k, Ffun), tspan, y0, opts); plot(t, y(:,1), b-, LineWidth, 1.2); xlabel(时间 t (s)); ylabel(位移 x (m)); grid on;(t,y) sdofODE(t,y,m,c,k,Ffun)是匿名函数它把当前时刻 t 和状态 y 传给 sdofODE同时把外层变量 m、c、k、Ffun 捕获到回调里这样主程序里改参数无需改函数签名。opts 里的 RelTol 从默认 1e-6 收紧到 1e-8因为小阻尼系统的峰值持续很多个周期累计相位误差会放大AbsTol 给到 1e-10 是为了让接近零的位移阶段不被绝对误差限制。输出的 y 是两列第一列是位移第二列是速度。注意如果 MATLAB 脚本里直接写局部函数R2016b 起支持把 sdofODE 放在脚本末尾为了兼容旧版本建议单独保存为 sdofODE.m。参数与建议取值见表。参数含义建议取值 / 影响m质量kg和 k 一起决定固有频率c阻尼系数N·s/mc2√(mk) 时为临界阻尼k刚度N/m决定回复力大小tspan积分时间区间至少覆盖 3 个衰减时间常数y0初始位移和初始速度初速不为零时瞬态幅值会明显改变opts误差控制结构高频率或小阻尼时收紧 AbsTol常见错误是只设置 RelTol 不设置 AbsTol。对位移量级为毫米的模型AbsTol 如果仍是默认 1e-6会把 1e-7 以下的真实微振动直接当零处理。另一个常见错误是搞混 y 的列顺序导致把速度画成了位移。2.4 用固有频率和衰减周期验证方程写没写错求解之前先计算理论值数值解才谈得上可信。用下面几行wn sqrt(k/m); fn wn/(2*pi); zeta c/(2*sqrt(m*k)); wd wn*sqrt(1-zeta^2); Td 2*pi/wd; fprintf(固有频率 fn %.3f Hz, 阻尼比 zeta %.4f, 衰减周期 Td %.4f s\n, ... fn, zeta, Td);在欠阻尼状态下自由振动的相邻正峰值时间差应接近 Td。如果你从 plot 里量出来的间距和 Td 差超过 2%要检查参数单位比如位移用 mm力用 N刚度应当按 N/mm 换算而不是混用 m 和 mm。若 zeta 大于 1响应不再振荡上面 wd 的公式变成虚数此时就不要用 Td 验证改看数值解是否单调趋近 F/k。写完这四步单自由度动力学方程在 MATLAB 里就跑通了。下一步是处理各种外部激励这比看自由衰减曲线更有工程意义。3. 单自由度动力学方程在MATLAB中的激励与响应三种典型输入怎么设工程中很少只有初始位移作用更多是持续的力输入。3.1 激励函数句柄的统一写法因为 sdofODE 里的外激励是 Ffun(t)所以不同载荷只改一行。自由振动、简谐激励和阶跃载荷的句柄如下% 自由振动 Ffun_free (t) 0; % 简谐激励幅值 2 N圆频率 8 rad/s F0 2; w 8; Ffun_harmonic (t) F0 * sin(w*t); % 阶跃载荷t0 时持续 2 N F0 2; Ffun_step (t) F0 * (t 0);简谐激励是最常用的输入注意 MATLAB 的 sin 默认角度单位是弧度w 的单位是 rad/s 时sin(w*t) 的周期 T 2π/w。阶跃载荷的写法在数学上没问题但 t0 处是不连续的ode45 会在断点附近自动加密步长。如果发现这段积分很慢可以把阶跃看作“0 到 t0 无激励t0 之后有激励”两个阶段前一段结束时的状态作为后一段初值继续求解这样能避免事件检测带来的额外开销。3.2 单次求解和稳态峰值提取把上一章的求解过程封装成函数方便后面扫描。函数返回时间、状态和稳态峰值function [t, y, Amp] run_sdof(m, c, k, Ffun, y0, tspan) opts odeset(RelTol, 1e-8, AbsTol, 1e-10); [t, y] ode45((t,y) sdofODE(t, y, m, c, k, Ffun), tspan, y0, opts); n length(t); idx floor(0.7*n):n; % 后 30% 作为稳态区间 Amp max(abs(y(idx,1))); end这里用后 30% 的数据点计算稳态振幅。对一般激励频率不是特别接近固有频率的情况足够但接近共振时瞬态衰减很慢假设 5 倍时间常数后进入稳态更好。可以把 0.7 换成基于衰减时间常数 τ1/(ζω_n) 的自动计算在调用前根据 c、k、m 预估。稳态区间太短会把瞬态峰值误当成稳态峰值导致幅频曲线峰值偏大。使用自定义函数时确保 sdofODE 也已经放在 MATLAB 搜索路径中或者和 run_sdof 写在同一个脚本里。3.3 幅频响应扫描的完整代码用 run_sdof 做频率扫描一次性得到不同激励频率下的响应幅值m 1; k 100; c 0.2; wn sqrt(k/m); F0 1; y0 [0 0]; freq linspace(0.4*wn, 1.6*wn, 80); Amp zeros(size(freq)); for i 1:length(freq) w freq(i); Ffun (t) F0 * sin(w*t); [~, ~, A] run_sdof(m, c, k, Ffun, y0, [0 100]); Amp(i) A; end plot(freq/wn, Amp, ro-); xlabel(频率比 r \omega/\omega_n); ylabel(稳态振幅 (m)); grid on;linspace(0.4*wn, 1.6*wn, 80)把激励圆频率从 0.4 倍固有频率平滑扫到 1.6 倍80 个点一般足够看到共振峰。每个频率点都重新调用 ode45所以总共 80 次数值积分数据量不夸张。需要更快时用 parfor 把 for 换成 parfor但要注意 Ffun 在并行池中序列化参数 m、c、k 会作为广播变量传给每个 worker如果内存紧张反而串行更快。频率比是动力学里的标准无量纲量。从绘制出的曲线可以看到峰值出现在 r 略小于 1 的位置而不是 r1。这个偏移大小由阻尼比决定下一章会专门对比不同阻尼比曲线。如果你的扫描结果峰值特别尖锐可能是阻尼比很小且扫描点不够密这时应在 wn 附近加密频率点例如在共振区单独用更小的步长。3.4 不同激励类型该怎么选仿真时长仿真时长不能用同一个经验值。下表是三类激励下 tspan 的建议。激励类型tspan 建议需要观察的量自由振动3~5 倍衰减时间常数衰减率是否和理论一致简谐激励至少 20 个激励周期稳态幅值是否稳定阶跃载荷到响应超过稳态值的 98% 为止峰值时间和过冲量阶跃响应的峰值时间大约在 0.5 到 1 个固有周期附近过冲量只由阻尼比控制和激励幅值大小无关。这些量在结构设计中常被直接用作评价指标所以很值得在 MATLAB 里精确求解。接下来我们把阻尼比和频率比同时作为变量看参数扫描怎么做。4. 单自由度动力学方程的MATLAB参数扫描阻尼比和频率比对振动峰值的影响单自由度系统最有价值的输出之一是幅频响应曲线族。4.1 阻尼比扫描一批曲线代替多次试算保持质量和刚度不变阻尼比从 0.02 到 0.2 取四个值每个值对应不同 c。MATLAB 代码如下m 1; k 200; wn sqrt(k/m); zetaVec [0.02 0.05 0.1 0.2]; freqRatio linspace(0.4, 1.6, 100); F0 1; y0 [0 0]; figure; hold on; for iZ 1:length(zetaVec) zeta zetaVec(iZ); c 2*zeta*sqrt(m*k); AmpVec zeros(size(freqRatio)); for iR 1:length(freqRatio) w freqRatio(iR)*wn; Ffun (t) F0*sin(w*t); tspan [0 max(60/w, 50/(zeta*wn))]; [~, ~, AmpVec(iR)] run_sdof(m, c, k, Ffun, y0, tspan); end plot(freqRatio, AmpVec/(F0/k), LineWidth, 1.5); end xlabel(频率比 \omega/\omega_n); ylabel(位移放大因子 x/(F0/k)); legend(arrayfun((z) sprintf(\\zeta%.2f, z), zetaVec, uni, 0)); hold off; grid on;tspan 里用max(60/w, 50/(zeta*wn))作为仿真终点两个表达式分别保证激励周期数和衰减时间常数都足够。小阻尼ζ0.02时时间常数很长50/(ζω_n) 会覆盖整个瞬态衰减频率比 r 达到 1.0 附近时瞬态衰减慢需要额外关注。AmpVec/(F0/k)把振幅除以静位移得到无量纲放大因子这是动力学中比绝对位移更通用的指标。四条曲线放在同一张图上可以看到峰值随阻尼比增大而显著变矮。4.2 有阻尼共振峰位置最大响应不在固有频率处幅频响应峰值对应的频率其实不是 ω_n而是 ω_n sqrt(1-2ζ²)。验证一下zeta 0.1; r_peak sqrt(1 - 2*zeta^2); fprintf(峰值频率比 r_peak %.4f\n, r_peak);当 ζ0.1 时r_peak 约为 0.98995和固有频率只差 1%当 ζ0.2 时r_peak 约为 0.9798差别仍不大。工程上常把峰值频率直接称为共振频率但精确辨识阻尼比时要按这个偏移做修正。数值扫描得到的峰值位置还会受频率分辨率影响如果扫描步长是 0.01那么观测到的峰值频率误差至少为半个步长。峰值半功率法会利用幅值下降 3dB 对应的频率差所以峰值位置偏移对阻尼辨识影响不大下面给出方法。4.3 半功率带宽法在MATLAB里辨识阻尼比假设已经扫描得到频率数组 freq 和幅值数组 Amp。峰值频率附近按半功率点计算阻尼比的代码[peak, idx_peak] max(Amp); half peak/sqrt(2); above Amp half; idx_up find(diff(above)1); idx_down find(diff(above)-1); if ~isempty(idx_up) ~isempty(idx_down) f1 interp1(Amp(idx_up:idx_up1), freq(idx_up:idx_up1), half); f2 interp1(Amp(idx_down:idx_down1), freq(idx_down:idx_down1), half); zeta_est (f2 - f1) / (2 * freq(idx_peak)); else zeta_est NaN; endhalf peak/sqrt(2)对应幅值下降 3 dB 的点。above标记哪些频率点高于半功率阈值find(diff(above)1)找到第一次由低到高的转折find(diff(above)-1)找到由高到低的转折。interp1做一维插值补足扫描点之间的空白。如果扫描分辨率太粗导致阈值附近没有采样点idx_up或idx_down为空结果就是 NaN此时必须在峰值附近加密频率点。阻尼比很小时半功率带宽很窄半功率点可能落在扫描步长内部所以这一步实际也在检验你的频率扫描够不够细。对数值仿真来说直接在频响曲线上做半功率带宽是可行的对实验数据需要先排除噪声和泄漏造成的毛刺。4.4 阻尼比、质量和刚度对响应的不同影响下表总结了三个基础参数改变时单自由度系统响应的变化方向。参数变化固有频率峰值幅值峰值频率质量 m 增大降低不变或降低需看阻尼比向低偏移刚度 k 增大升高静位移 F0/k 减小放大因子不变向高偏移阻尼比 ζ 增大几乎不变明显降低略向低偏移表格里需要说明若 c 固定m 增大会使 ζ 减小峰值反而可能升高若 ζ 固定m 对放大因子影响不大。这就是为什么在做参数设计时不能单独说“增加质量就能减小振动”要同时考虑阻尼比的变化。这一章的扫描脚本已经能处理绝大多数单自由度参数分析任务。如果只需要理论幅频响应直接用解析式 1/sqrt((1-r²)² (2ζr)²) 会比 80 次 ode45 快得多数值解的意义在于验证瞬态过程和辨识实验参数。后面我们讨论如何验证这些数值结果避免被错误的单位或积分设置带偏。5. 单自由度动力学方程的MATLAB验证与排错解析解对照与误差来源数值解不是可信的充分条件。我每次写完新的动力学求解器都会先用解析解对照一次再做参数扫描。5.1 用解析解验证ode45的瞬态响应欠阻尼自由振动的解析解是x(t) e^(-ζω_n t) [A cos(ω_d t) B sin(ω_d t)]其中 ω_d ω_n sqrt(1-ζ²)A x0B (v0 ζω_n x0)/ω_d。用代码生成理论曲线和 ode45 输出叠加m 1; k 100; c 0.2; wn sqrt(k/m); zeta c/(2*sqrt(m*k)); wd wn*sqrt(1-zeta^2); x0 0.02; v0 0; A x0; B (v0 zeta*wn*x0)/wd; x_analytic exp(-zeta*wn*t) .* (A*cos(wd*t) B*sin(wd*t)); plot(t, y(:,1), b-, t, x_analytic, r--); legend(ode45, 解析解);这里直接用了 t 作为理论曲线的横坐标所以必须确保 t 与 y 是同一组时间点。如果曲线完全重合说明状态空间和参数传递正确。如果只在小振幅处有偏差优先检查 AbsTol如果后半段误差越来越大优先检查 RelTol 和单位转换。5.2 一个能快速定位问题的排错顺序排错时可以按三个来源定位。现象可能原因检查点峰值比理论低稳态区间选太短idx 起始点后移频率偏移扫描步长太大在共振区加密频率点相位持续漂移刚度或质量单位错误复核 k/m 量纲最后一个技巧是用 odeset 里的 Stats 查看积分步数。小阻尼系统会提交大量步数这是正常的如果步数少得异常且结果不像解析解说明误差控制可能被放宽了或者状态方程写成了稳定但错误的系统。本文还有配套的精品资源点击获取
MATLAB单自由度振动仿真：ode45状态空间与幅频响应分析
简介面向机械工程与齿轮传动系统动力学分析的一份 MATLAB 源码包主要解决直齿轮副非线性动力学建模与求解问题。资源以单自由度简化模型为基础基于牛顿第二定律建立包含质量、加速度、力与扭矩的运动方程并通过相图刻划位置与速度的状态轨迹帮助判断系统稳定性与周期性利用傅里叶变换分析振动信号识别由啮合冲击、制造误差或载荷波动引起的频谱特征服务于齿轮故障诊断与结构优化。包内为1个可运行的 .m 脚本压缩后仅 1KB体积小巧代码结构清晰便于读者快速修改参数并重新运行。已有 678 人学习适合需要入门齿轮非线性动力学或希望借助编程完成建模、求解与可视化验证的工程技术人员也可在此基础上扩展含间隙、时变啮合刚度等复杂因素进一步研究分岔与混沌现象。1. 单自由度动力学方程为什么值得用MATLAB跑一遍一个弹簧、一个质量、一个阻尼器单自由度动力学方程看起来比 Navier-Stokes 简单得多但真用手推解简谐激励的瞬态响应时初始位移和初始速度决定的那两项常数经常让人在积分结果上栽跟头。这也是为什么我始终建议先用 MATLAB 把 m x c x k x F(t) 完整跑一遍不是为了得到一条曲线而是为了验证你对这个系统的直觉——共振峰值并不是恰好出现在固有频率处而是在略偏左的地方。这类模型在机械振动、结构减震、动态测试和控制系统里都是被反复调用的基础单元新手可以用它掌握 ode45 和状态空间老手可以用它快速评估算法参数。2. 把单自由度动力学方程改写成MATLAB可执行的ode45状态空间在写任何求解脚本之前先把方程标准化这一步常被跳过。2.1 运动方程标准化先除以质量还是保留原参数原始方程是m x(t) c x(t) k x(t) F(t)我一般习惯先除以 m写成x 2 ζ ω_n x ω_n² x f(t)其中 ω_n sqrt(k/m)ζ c / (2 sqrt(m k))。这样做的直接好处是二阶项系数为 1ode45 内部的相对误差控制对各项量级更敏感如果 m 是 1e-3 而 k 是 1e6直接用原始系数写回调数值梯度会跨好几个量级容易触发不必要的步长加密甚至让结果看起来像高频振荡。更重要的是ζ 和 ω_n 这两个参数直接决定响应的形态看到 ζ 小于 1就能预期响应是衰减振荡而不是单调趋近这对排查求解结果是很有用的。也要注意除以 m 只适合单自由度方程本身。如果后面扩展到多自由度系统质量矩阵存在耦合不能简单逐行除而应使用状态空间下的质量矩阵求逆。单自由度时不需要计较这一步但思路要清楚。2.2 状态空间两个一阶方程替换一个二阶方程ode45 只能处理一阶常微分方程组而动力学方程是二阶的所以要做降阶。令 y1 xy2 x则由 y1 y2 和原方程可得y1 y2 y2 (F(t) - c y2 - k y1) / m这是一个标准的二维状态空间。外激励 F(t) 用函数句柄表示这样自由振动、简谐激励、随机载荷都不需要改导数函数只换句柄即可。如果 F(t) 是线性系统的外部输入这个写法也方便以后改成时变参数比如 c(t) 或 k(t)。2.3 最小可运行代码sdof函数加ode45调用把状态空间写成 MATLAB 函数建议单独保存为 sdofODE.m方便多个脚本共用function dydt sdofODE(t, y, m, c, k, Ffun) % 单自由度动力学方程的状态空间形式 % y(1) 位移, y(2) 速度 dydt zeros(2,1); dydt(1) y(2); dydt(2) (Ffun(t) - c*y(2) - k*y(1)) / m; end然后调用 ode45 并绘图m 1.2; % 质量 kg c 0.4; % 阻尼系数 N*s/m k 30; % 刚度 N/m Ffun (t) 0; % 自由振动 tspan [0 10]; y0 [0.02; 0]; opts odeset(RelTol, 1e-8, AbsTol, 1e-10); [t, y] ode45((t,y) sdofODE(t, y, m, c, k, Ffun), tspan, y0, opts); plot(t, y(:,1), b-, LineWidth, 1.2); xlabel(时间 t (s)); ylabel(位移 x (m)); grid on;(t,y) sdofODE(t,y,m,c,k,Ffun)是匿名函数它把当前时刻 t 和状态 y 传给 sdofODE同时把外层变量 m、c、k、Ffun 捕获到回调里这样主程序里改参数无需改函数签名。opts 里的 RelTol 从默认 1e-6 收紧到 1e-8因为小阻尼系统的峰值持续很多个周期累计相位误差会放大AbsTol 给到 1e-10 是为了让接近零的位移阶段不被绝对误差限制。输出的 y 是两列第一列是位移第二列是速度。注意如果 MATLAB 脚本里直接写局部函数R2016b 起支持把 sdofODE 放在脚本末尾为了兼容旧版本建议单独保存为 sdofODE.m。参数与建议取值见表。参数含义建议取值 / 影响m质量kg和 k 一起决定固有频率c阻尼系数N·s/mc2√(mk) 时为临界阻尼k刚度N/m决定回复力大小tspan积分时间区间至少覆盖 3 个衰减时间常数y0初始位移和初始速度初速不为零时瞬态幅值会明显改变opts误差控制结构高频率或小阻尼时收紧 AbsTol常见错误是只设置 RelTol 不设置 AbsTol。对位移量级为毫米的模型AbsTol 如果仍是默认 1e-6会把 1e-7 以下的真实微振动直接当零处理。另一个常见错误是搞混 y 的列顺序导致把速度画成了位移。2.4 用固有频率和衰减周期验证方程写没写错求解之前先计算理论值数值解才谈得上可信。用下面几行wn sqrt(k/m); fn wn/(2*pi); zeta c/(2*sqrt(m*k)); wd wn*sqrt(1-zeta^2); Td 2*pi/wd; fprintf(固有频率 fn %.3f Hz, 阻尼比 zeta %.4f, 衰减周期 Td %.4f s\n, ... fn, zeta, Td);在欠阻尼状态下自由振动的相邻正峰值时间差应接近 Td。如果你从 plot 里量出来的间距和 Td 差超过 2%要检查参数单位比如位移用 mm力用 N刚度应当按 N/mm 换算而不是混用 m 和 mm。若 zeta 大于 1响应不再振荡上面 wd 的公式变成虚数此时就不要用 Td 验证改看数值解是否单调趋近 F/k。写完这四步单自由度动力学方程在 MATLAB 里就跑通了。下一步是处理各种外部激励这比看自由衰减曲线更有工程意义。3. 单自由度动力学方程在MATLAB中的激励与响应三种典型输入怎么设工程中很少只有初始位移作用更多是持续的力输入。3.1 激励函数句柄的统一写法因为 sdofODE 里的外激励是 Ffun(t)所以不同载荷只改一行。自由振动、简谐激励和阶跃载荷的句柄如下% 自由振动 Ffun_free (t) 0; % 简谐激励幅值 2 N圆频率 8 rad/s F0 2; w 8; Ffun_harmonic (t) F0 * sin(w*t); % 阶跃载荷t0 时持续 2 N F0 2; Ffun_step (t) F0 * (t 0);简谐激励是最常用的输入注意 MATLAB 的 sin 默认角度单位是弧度w 的单位是 rad/s 时sin(w*t) 的周期 T 2π/w。阶跃载荷的写法在数学上没问题但 t0 处是不连续的ode45 会在断点附近自动加密步长。如果发现这段积分很慢可以把阶跃看作“0 到 t0 无激励t0 之后有激励”两个阶段前一段结束时的状态作为后一段初值继续求解这样能避免事件检测带来的额外开销。3.2 单次求解和稳态峰值提取把上一章的求解过程封装成函数方便后面扫描。函数返回时间、状态和稳态峰值function [t, y, Amp] run_sdof(m, c, k, Ffun, y0, tspan) opts odeset(RelTol, 1e-8, AbsTol, 1e-10); [t, y] ode45((t,y) sdofODE(t, y, m, c, k, Ffun), tspan, y0, opts); n length(t); idx floor(0.7*n):n; % 后 30% 作为稳态区间 Amp max(abs(y(idx,1))); end这里用后 30% 的数据点计算稳态振幅。对一般激励频率不是特别接近固有频率的情况足够但接近共振时瞬态衰减很慢假设 5 倍时间常数后进入稳态更好。可以把 0.7 换成基于衰减时间常数 τ1/(ζω_n) 的自动计算在调用前根据 c、k、m 预估。稳态区间太短会把瞬态峰值误当成稳态峰值导致幅频曲线峰值偏大。使用自定义函数时确保 sdofODE 也已经放在 MATLAB 搜索路径中或者和 run_sdof 写在同一个脚本里。3.3 幅频响应扫描的完整代码用 run_sdof 做频率扫描一次性得到不同激励频率下的响应幅值m 1; k 100; c 0.2; wn sqrt(k/m); F0 1; y0 [0 0]; freq linspace(0.4*wn, 1.6*wn, 80); Amp zeros(size(freq)); for i 1:length(freq) w freq(i); Ffun (t) F0 * sin(w*t); [~, ~, A] run_sdof(m, c, k, Ffun, y0, [0 100]); Amp(i) A; end plot(freq/wn, Amp, ro-); xlabel(频率比 r \omega/\omega_n); ylabel(稳态振幅 (m)); grid on;linspace(0.4*wn, 1.6*wn, 80)把激励圆频率从 0.4 倍固有频率平滑扫到 1.6 倍80 个点一般足够看到共振峰。每个频率点都重新调用 ode45所以总共 80 次数值积分数据量不夸张。需要更快时用 parfor 把 for 换成 parfor但要注意 Ffun 在并行池中序列化参数 m、c、k 会作为广播变量传给每个 worker如果内存紧张反而串行更快。频率比是动力学里的标准无量纲量。从绘制出的曲线可以看到峰值出现在 r 略小于 1 的位置而不是 r1。这个偏移大小由阻尼比决定下一章会专门对比不同阻尼比曲线。如果你的扫描结果峰值特别尖锐可能是阻尼比很小且扫描点不够密这时应在 wn 附近加密频率点例如在共振区单独用更小的步长。3.4 不同激励类型该怎么选仿真时长仿真时长不能用同一个经验值。下表是三类激励下 tspan 的建议。激励类型tspan 建议需要观察的量自由振动3~5 倍衰减时间常数衰减率是否和理论一致简谐激励至少 20 个激励周期稳态幅值是否稳定阶跃载荷到响应超过稳态值的 98% 为止峰值时间和过冲量阶跃响应的峰值时间大约在 0.5 到 1 个固有周期附近过冲量只由阻尼比控制和激励幅值大小无关。这些量在结构设计中常被直接用作评价指标所以很值得在 MATLAB 里精确求解。接下来我们把阻尼比和频率比同时作为变量看参数扫描怎么做。4. 单自由度动力学方程的MATLAB参数扫描阻尼比和频率比对振动峰值的影响单自由度系统最有价值的输出之一是幅频响应曲线族。4.1 阻尼比扫描一批曲线代替多次试算保持质量和刚度不变阻尼比从 0.02 到 0.2 取四个值每个值对应不同 c。MATLAB 代码如下m 1; k 200; wn sqrt(k/m); zetaVec [0.02 0.05 0.1 0.2]; freqRatio linspace(0.4, 1.6, 100); F0 1; y0 [0 0]; figure; hold on; for iZ 1:length(zetaVec) zeta zetaVec(iZ); c 2*zeta*sqrt(m*k); AmpVec zeros(size(freqRatio)); for iR 1:length(freqRatio) w freqRatio(iR)*wn; Ffun (t) F0*sin(w*t); tspan [0 max(60/w, 50/(zeta*wn))]; [~, ~, AmpVec(iR)] run_sdof(m, c, k, Ffun, y0, tspan); end plot(freqRatio, AmpVec/(F0/k), LineWidth, 1.5); end xlabel(频率比 \omega/\omega_n); ylabel(位移放大因子 x/(F0/k)); legend(arrayfun((z) sprintf(\\zeta%.2f, z), zetaVec, uni, 0)); hold off; grid on;tspan 里用max(60/w, 50/(zeta*wn))作为仿真终点两个表达式分别保证激励周期数和衰减时间常数都足够。小阻尼ζ0.02时时间常数很长50/(ζω_n) 会覆盖整个瞬态衰减频率比 r 达到 1.0 附近时瞬态衰减慢需要额外关注。AmpVec/(F0/k)把振幅除以静位移得到无量纲放大因子这是动力学中比绝对位移更通用的指标。四条曲线放在同一张图上可以看到峰值随阻尼比增大而显著变矮。4.2 有阻尼共振峰位置最大响应不在固有频率处幅频响应峰值对应的频率其实不是 ω_n而是 ω_n sqrt(1-2ζ²)。验证一下zeta 0.1; r_peak sqrt(1 - 2*zeta^2); fprintf(峰值频率比 r_peak %.4f\n, r_peak);当 ζ0.1 时r_peak 约为 0.98995和固有频率只差 1%当 ζ0.2 时r_peak 约为 0.9798差别仍不大。工程上常把峰值频率直接称为共振频率但精确辨识阻尼比时要按这个偏移做修正。数值扫描得到的峰值位置还会受频率分辨率影响如果扫描步长是 0.01那么观测到的峰值频率误差至少为半个步长。峰值半功率法会利用幅值下降 3dB 对应的频率差所以峰值位置偏移对阻尼辨识影响不大下面给出方法。4.3 半功率带宽法在MATLAB里辨识阻尼比假设已经扫描得到频率数组 freq 和幅值数组 Amp。峰值频率附近按半功率点计算阻尼比的代码[peak, idx_peak] max(Amp); half peak/sqrt(2); above Amp half; idx_up find(diff(above)1); idx_down find(diff(above)-1); if ~isempty(idx_up) ~isempty(idx_down) f1 interp1(Amp(idx_up:idx_up1), freq(idx_up:idx_up1), half); f2 interp1(Amp(idx_down:idx_down1), freq(idx_down:idx_down1), half); zeta_est (f2 - f1) / (2 * freq(idx_peak)); else zeta_est NaN; endhalf peak/sqrt(2)对应幅值下降 3 dB 的点。above标记哪些频率点高于半功率阈值find(diff(above)1)找到第一次由低到高的转折find(diff(above)-1)找到由高到低的转折。interp1做一维插值补足扫描点之间的空白。如果扫描分辨率太粗导致阈值附近没有采样点idx_up或idx_down为空结果就是 NaN此时必须在峰值附近加密频率点。阻尼比很小时半功率带宽很窄半功率点可能落在扫描步长内部所以这一步实际也在检验你的频率扫描够不够细。对数值仿真来说直接在频响曲线上做半功率带宽是可行的对实验数据需要先排除噪声和泄漏造成的毛刺。4.4 阻尼比、质量和刚度对响应的不同影响下表总结了三个基础参数改变时单自由度系统响应的变化方向。参数变化固有频率峰值幅值峰值频率质量 m 增大降低不变或降低需看阻尼比向低偏移刚度 k 增大升高静位移 F0/k 减小放大因子不变向高偏移阻尼比 ζ 增大几乎不变明显降低略向低偏移表格里需要说明若 c 固定m 增大会使 ζ 减小峰值反而可能升高若 ζ 固定m 对放大因子影响不大。这就是为什么在做参数设计时不能单独说“增加质量就能减小振动”要同时考虑阻尼比的变化。这一章的扫描脚本已经能处理绝大多数单自由度参数分析任务。如果只需要理论幅频响应直接用解析式 1/sqrt((1-r²)² (2ζr)²) 会比 80 次 ode45 快得多数值解的意义在于验证瞬态过程和辨识实验参数。后面我们讨论如何验证这些数值结果避免被错误的单位或积分设置带偏。5. 单自由度动力学方程的MATLAB验证与排错解析解对照与误差来源数值解不是可信的充分条件。我每次写完新的动力学求解器都会先用解析解对照一次再做参数扫描。5.1 用解析解验证ode45的瞬态响应欠阻尼自由振动的解析解是x(t) e^(-ζω_n t) [A cos(ω_d t) B sin(ω_d t)]其中 ω_d ω_n sqrt(1-ζ²)A x0B (v0 ζω_n x0)/ω_d。用代码生成理论曲线和 ode45 输出叠加m 1; k 100; c 0.2; wn sqrt(k/m); zeta c/(2*sqrt(m*k)); wd wn*sqrt(1-zeta^2); x0 0.02; v0 0; A x0; B (v0 zeta*wn*x0)/wd; x_analytic exp(-zeta*wn*t) .* (A*cos(wd*t) B*sin(wd*t)); plot(t, y(:,1), b-, t, x_analytic, r--); legend(ode45, 解析解);这里直接用了 t 作为理论曲线的横坐标所以必须确保 t 与 y 是同一组时间点。如果曲线完全重合说明状态空间和参数传递正确。如果只在小振幅处有偏差优先检查 AbsTol如果后半段误差越来越大优先检查 RelTol 和单位转换。5.2 一个能快速定位问题的排错顺序排错时可以按三个来源定位。现象可能原因检查点峰值比理论低稳态区间选太短idx 起始点后移频率偏移扫描步长太大在共振区加密频率点相位持续漂移刚度或质量单位错误复核 k/m 量纲最后一个技巧是用 odeset 里的 Stats 查看积分步数。小阻尼系统会提交大量步数这是正常的如果步数少得异常且结果不像解析解说明误差控制可能被放宽了或者状态方程写成了稳定但错误的系统。本文还有配套的精品资源点击获取
