第一次在实验室给一根钢尺贴加速度计、拿力锤敲下去频谱仪上冒出第一个尖峰的时候我脑子里冒出一个很自然的想法这个尖峰对应的频率能不能用一根悬臂梁直接算出来答案是能。而且当你把悬臂梁连续体振动模型在 Matlab 里完整跑通以后才会真正理解“连续体”这三个字背后那套偏微分方程的威力。这篇博文就是我实现这个模型的全过程记录从欧拉-伯努利梁方程的推导、特征方程的数值求解到振型函数、模态叠加法算时域响应最后落在结果验证和几个高频踩坑点上。适合正在学结构动力学、做课程设计、或者想找一套基准解来验证自己有限元模型的同学参考。1. 为什么悬臂梁连续体振动模型值得从零实现一遍1.1 悬臂梁在工程里到底代表什么桥墩、吊车臂、机翼的简化模型、机床刀架这些结构都能抽象成一根悬臂梁。甚至你手里那把钢尺在绝大多数情况下也是标准的悬臂梁一端固定一端自由边界条件清清楚楚。悬臂梁在结构动力学里的地位有点像“hello world”在编程里的地位但它不是个玩具。因为它的边界条件明确、理论解存在又不至于简单到让人忽略了振动问题的本质。很多真实工程问题在细长结构假设下都可以先拿悬臂梁估算前几阶固有频率有了这个数量级概念再决定要不要上更复杂的模型。我建议每个学振动的人亲手实现一遍悬臂梁连续体模型。不是因为课程大纲写了而是它能让你把下面这条逻辑链完整走一遍建立运动微分方程分离变量把偏微分方程转成特征值问题引入边界条件得到特征方程数值求解特征根得到固有频率和振型用模态叠加法计算动力响应。这条链几乎覆盖了所有线性振动分析的核心方法。搞懂它之后以后无论做多自由度系统、有限元还是读文献里的频响函数公式思路都不会乱。1.2 连续体、集中质量和有限元三者到底差在哪我经常被问一个问题既然有有限元了为什么还要研究连续体模型这个问题值得掰开揉碎说清楚。集中质量模型是把梁离散成若干个质点和无质量弹簧自由度有限。优点是直观写代码特别快缺点是精度依赖离散点数而且它的高阶模态对应的是虚构的弹簧-质量组合和真实梁的高频行为之间会有系统性偏差。比如你把一根梁分 5 段第 5 阶频率可能和真实值差出百分之十几。连续体模型用偏微分方程描述整个梁的变形自由度是无限的解析解给出完整的频率谱。欧拉-伯努利梁方程就是连续体模型最常见的代表适合细长梁的低频振动分析。它的最大价值是结果可以作为“基准解”用来检验各种数值方法的收敛性。有限元模型本质上是把连续体再离散回有限自由度但它的“聪明”之处在于基函数构造和收敛性分析所以同样自由度下有限元比简单的集中质量模型精确得多。在复杂几何、变截面、复杂边界条件下有限元是首选但如果你做的是一根等截面直梁连续体模型的解析解几乎零成本就能给你一个非常可靠的结果。我的习惯是先用连续体模型估算再上有限元和实验对照。三者不是互斥关系而是互为验证。2. 欧拉-伯努利梁方程从物理到特征值问题2.1 四阶偏微分方程是怎么推出来的很多人一看到振动方程就背公式但我建议这次推导跟着走一遍因为后面代码里的“特征方程”和“振型系数”全靠这一步。取梁上的一段微元 dx左侧有剪力 Q 和弯矩 M右侧是 QdQ 和 MdM梁上作用有分布载荷 p(x,t) 和惯性力 ρA dx · ∂²w/∂t²。对微元列力平衡和力矩平衡再忽略转动惯量的影响可以得到∂Q/∂x -p ρA ∂²w/∂t²∂M/∂x Q再把材料力学里的弯矩-曲率关系 M EI ∂²w/∂x² 代进去消去 Q 和 M就得到了欧拉-伯努利梁的横向振动方程EI ∂⁴w/∂x⁴ ρA ∂²w/∂t² p(x,t)这个方程看着简单背后有三个假设平截面假设变形后截面保持平面并垂直于中性轴、忽略剪切变形、忽略转动惯量。所以它的适用条件很明确——细长梁、低频段。等到了验证章节我会专门说这些假设什么时候开始失效。2.2 边界条件与特征方程的形成过程自由振动时 p0。设 w(x,t) W(x) sin(ωt)代入之后空间部分会变成EI W(x) ω² ρA W(x)令 β⁴ ω² ρA / (EI)方程就变成 W β⁴ W。它的通解是W(x) C1 cosh(βx) C2 sinh(βx) C3 cos(βx) C4 sin(βx)悬臂梁的边界条件分两组。固定端 x0 处位移为零、转角为零即W(0)0W(0)0自由端 xL 处弯矩为零、剪力为零即W(L)0W(L)0把通解代入四个边界条件。先看固定端两个条件可以得到 C3-C1、C4-C2于是 W(x) 简化为W(x) C1(cosh βx − cos βx) C2(sinh βx − sin βx)再看自由端的弯矩和剪力条件。把上式代进去整理成一个关于 C1、C2 的线性方程组。非零解要求系数矩阵行列式为零最终得到悬臂梁的特征方程cos(λ) · cosh(λ) 1 0其中 λ βL。这个方程没有闭式解它的根决定了梁的固有频率所以我们必须用数值方法求解。很多初学的人会在这一步之间直接背诵结论我建议至少在纸上写一遍行列式展开的过程因为自由端弯矩、剪力的正负号在后面写振型函数时真的会害人。2.3 前几阶固有频率与一个算例特征方程的根记为 λn。前六阶数值如下表所示高阶根渐近趋近于 (n−0.5)π。阶数 nλnλn²11.8751043.51601524.69409122.03449137.85475761.697214410.995541120.901918514.137168199.859454617.278760298.555548固有频率公式是ωn (λn/L)² · √(EI/ρA)频率值 fn ωn / (2π)。我做示例时常用的参数是L1m矩形截面宽 b10mm、高 h10mmE210GPa密度 ρ7850 kg/m³。把这些值代进去前六阶频率大约为8.36Hz、52.37Hz、146.6Hz、287.4Hz、475.1Hz、709.6Hz。顺带说一个我记忆里的常识悬臂梁第一阶和第二阶频率之比大约是 1:6.27而简支梁是 1:4。差异这么大完全是边界条件不同造成的。这个比例可以当经验值记着用来快速判断一个动力学问题里悬臂假设是否合理。3. Matlab实现求根、算频率、画振型3.1 特征方程怎么在 Matlab 里稳准狠地求根f(λ) cosλ·coshλ 1 在正半轴上有无穷多个根直接用 fzero 必须提供根附近的区间。我的做法是先大步长扫描符号变化再在变化区间里调用 fzero。根间距大约 π所以扫描步长取 0.05 已经非常安全对前几十阶都不会漏根。function lambda findCantileverRoots(N) % 求解 cos(lambda)*cosh(lambda)1 0 的前 N 个正根 lambda zeros(N, 1); x 0.01; % 从略大于0开始 dx 0.05; % 扫描步长 k 0; while k N f0 cos(x)*cosh(x) 1; x1 x dx; f1 cos(x1)*cosh(x1) 1; if f0*f1 0 k k 1; lambda(k) fzero((t) cos(t)*cosh(t)1, [x, x1]); end x x1; end end这里有个细节fzero 的区间参数要求两端函数值异号上面的符号检查正好保证了这一点。实际写代码时建议先画一张 f(λ) 的图亲眼确认前几阶根的大概位置再开始扫。这个习惯能省掉夜里调试的时间。如果你不喜欢循环加区间扫描也有更简洁的方案第 k 阶特征根的渐近位置在 (k−0.5)π 附近直接用 fzero 单点初值lambda(k) fzero((t) cos(t)*cosh(t)1, (k - 0.5)*pi);这个方案在低阶时初值离根也不远第一阶初值 1.5708真根 1.8751只要对初步结果做一个“单调递增、互不相等”的合理性检查就很好用。3.2 振型函数公式、代码与质量归一化求得 λn 之后振型函数可以写成Wn(x) cosh(βn x) − cos(βn x) − σn [sinh(βn x) − sin(βn x)]其中系数 σn 由自由端边界条件给出σn (cosh λn cos λn) / (sinh λn sin λn)对应 Matlab 函数function W cantileverModeShape(x, L, lambda) % 计算悬臂梁第n阶振型未归一化 beta lambda / L; bx beta .* x; sigma (cosh(lambda) cos(lambda)) / (sinh(lambda) sin(lambda)); W (cosh(bx) - cos(bx)) - sigma .* (sinh(bx) - sin(bx)); end直接用它画图会发现不同阶的振型幅值差异很大做模态叠加时也很不方便。所以我更推荐在得到原始振型后做质量归一化让每一阶都满足∫₀^L ρA · Wn²(x) dx 1归一化的代码是x linspace(0, L, 500).; % 坐标网格 W_raw cantileverModeShape(x, L, lambda(n)); Mn trapz(x, rhoA * W_raw.^2); % 模态质量 W_n W_raw / sqrt(Mn); % 质量归一化后的振型质量归一化之后后续模态方程里的广义质量全部变成 1广义力直接按投影算不需要每次除模态质量响应计算会清爽很多。3.3 一段跑通主流程的骨架代码下面这段主脚本把求根、频率、归一化振型串起来顺便画前四阶振型图整套流程五分钟之内就能跑通% 参数定义SI单位制 L 1.0; b 0.01; h 0.01; E 210e9; rho 7850; A b * h; I b * h^3 / 12; % 求前六阶特征根 lambda findCantileverRoots(6); % 计算固有频率 omega (lambda / L).^2 .* sqrt(E*I / (rho*A)); freq omega / (2*pi); % 坐标网格与振型矩阵 N 6; x linspace(0, L, 500).; W_mat zeros(length(x), N); for n 1:N W_raw cantileverModeShape(x, L, lambda(n)); Mn trapz(x, rho*A * W_raw.^2); W_mat(:, n) W_raw / sqrt(Mn); end % 绘制前四阶归一化振型 figure; for n 1:4 subplot(2,2,n); plot(x, W_mat(:, n), LineWidth, 1.5); title(sprintf(第%d阶振型, fn%.2f Hz, n, freq(n))); grid on; end画出来的前四阶振型应该是第一阶单向弯曲第二阶出现一个反弯点第三阶两个反弯点第四阶三个反弯点。固定端位移为零、斜率为零自由端弯矩为零对应曲率为零这些形状特征一眼就能对上。4. 从模态到响应让悬臂梁动起来4.1 模态叠加法的基本套路有了频率和振型就可以计算悬臂梁在任意动载荷下的响应。核心方法是模态叠加法思路是把真实位移按振型展开w(x,t) Σ qn(t) · Wn(x)把这一项代回偏微分方程利用振型关于质量矩阵和刚度矩阵的正交性方程会解耦成一个个独立的二阶常微分方程q̈n 2ζnωn q̇n ωn² qn Fn(t)其中广义力Fn(t) ∫₀^L f(x,t) · Wn(x) dx每个模态方程都可以单独数值积分最后再把所有 qn(t)·Wn(x) 叠加回来。这就是模态叠加法。需要强调的是无阻尼线性系统的振型一定满足正交性解耦才能成立如果阻尼是非比例阻尼模态方程就会耦合那就得直接解原始方程组了。4.2 集中力激励、阻尼与 ode45 求解实际试验里最常见的是力锤或激振器施加的集中力。设力作用在 xx0表达式写为f(x,t) F0 · sin(2π fd t) · δ(x − x0)利用 δ 函数的筛选性质广义力变成一个简单式Fn(t) F0 · Wn(x0) · sin(2π fd t)不需要积分直接把激励点位置的振型值取出来乘以力幅就行。这个便利是我建议做质量归一化的原因之一。阻尼的处理上工程里最常用的是模态阻尼比 ζn。如果没有实验数据我通常先给每阶 0.02 左右的统一阻尼比如果重点关注共振峰附近再按实测曲线拟合各阶阻尼。下面是一个用 ode45 求解前五阶模态响应并合成观测点位移的代码框架N 5; x0 0.3; % 激励点 x_meas 0.8; % 响应观测点 F0 10; % 激励幅值 fd freq(2); % 激励频率设为第二阶附近 zeta 0.02; % 统一模态阻尼比 t_out 0:0.01:5; % 固定输出时间轴 w_series zeros(size(t_out)); q_series zeros(length(t_out), N); for n 1:N wn omega(n); Wn_x0 cantileverModeShape(x0, L, lambda(n)); Wn_m cantileverModeShape(x_meas, L, lambda(n)); % 广义力这里用归一化振型但只差一个常数倍数影响不大 Fn (t) F0 * Wn_x0 * sin(2*pi*fd*t); % 每个模态方程转成一阶方程组 odefun (t,y) [y(2); Fn(t) - 2*zeta*wn*y(2) - wn^2*y(1)]; [~, y] ode45(odefun, t_out, [0; 0]); q_series(:, n) y(:, 1); w_series w_series y(:, 1) * Wn_m; end plot(t_out, w_series);几个实现上的细节用固定 t_out 作为 ode45 的输出时间点保证每阶模态响应时间轴一致后面叠加才不会错位初始条件给零体现从静止状态开始被激励起来的完整过程如果只关心稳态幅值可以取响应平稳段后评估或者直接用频响函数公式不需要积分。4.3 响应图怎么读以及截断阶数怎么定看响应时程时最开始那一段波形比较乱那是瞬态响应和稳态响应叠加的阶段。激励开始的几个周期里各阶模态都在“起振”互相干涉经过几个衰减时间常数后瞬态基本衰减掉剩下稳态简谐响应。工程上测频响函数时都喜欢加窗或者多次平均本质就是把瞬态影响压掉。截断阶数不是越多越好而要看你关心的频率范围。一般来说N 取到激励频率的 2~3 倍位置就够了。比如激励频率 50Hz梁前五阶频率覆盖到 475Hz那 N 取 5 基本没问题。如果你发现某阶广义坐标 qn 的幅值远小于其他阶说明这一阶模态对当前激励的贡献可以忽略直接砍掉不心疼。5. 结果验证怎么确认代码没算错5.1 振型正交性是最好的自检手段我每次写完振动相关的代码第一件事不是看频率对不对而是检查振型正交性。质量归一化之后应该有∫₀^L ρA · Wm(x) · Wn(x) dx δmn对应代码M_ij zeros(N, N); for i 1:N for j 1:N Wi W_mat(:, i); Wj W_mat(:, j); M_ij(i,j) trapz(x, rhoA * Wi .* Wj); end end disp(M_ij);理想情况下对角元应该是 1非对角元接近 0。用 trapz 做数值积分会引入非常小的误差所以对角元落在 0.9999~1.0001、非对角元在 1e-12 量级就算正常。如果你看到对角元明显不是 1说明归一化算错了如果非对角元有 1e-3 这种量级那振型函数的符号、系数 σn 很可能出了问题。这个检查能一次性暴露很大一类 bug。5.2 频率与文献值及有限元结果的对照频率本身也可以对照。教科书上的悬臂梁频率常数列出来对比任何一阶差出万分之几以上基本可以断定代码里哪一环出错了。如果手头有通用有限元软件用梁单元算一下同尺寸悬臂梁的前几阶频率效果更好。细长梁的欧拉-伯努利解析解和梁单元有限元解在前几阶应该非常接近第一阶往往差不到 0.1%第二阶也在零点几个百分点以内。等到第三、四阶剪切变形的影响开始显现解析解频率会略微偏高。梁越矮胖跨高比越小偏差出现得越早。这个现象不是代码 bug而是模型假设本身的边界。5.3 欧拉-伯努利模型的适用边界在哪什么情况下要小心连续体解析解我给自己定的经验线是跨高比 L/h 20 时前几阶用欧拉-伯努利模型没问题L/h 在 10~20 之间高阶频率误差可能超过 1%需要考虑剪切变形L/h 10建议直接上 Timoshenko 梁或者二维平面模型。另一个失效场景是高频段。即便是细长梁当阶数升高到波长接近截面尺寸时欧拉-伯努利模型给出的频率偏高因为真实梁存在剪切挠度和转动惯量。如果项目关注的是第 20 阶、第 50 阶这种高频动态响应那解析解的参考价值就要打个问号。变截面梁、轴向力作用下的梁、大变形几何非线性问题也都不再适合用这个简单连续体模型。遇到这些情况有限元是更务实的路线但连续体模型的解仍然是验证程序正确性最好的起点。6. 实操经验Matlab实现中的三个高频坑6.1 特征根漏求和 fzero 的 bracket 问题fzero 这个函数如果只给单点初值它有可能在求根过程中跳到别的根上去如果给区间但区间端点函数值同号会直接报错。所以扫描步长不能拍脑袋取太大的数最好先画一张 f(λ) 的曲线看看根的大概分布。根间隔约 π用 0.05 的步长扫描前六阶几乎不可能漏。唯一要小心的是 f(λ) 在 λ 很大的时候cosh(λ) 增长极快计算 f 值本身没问题但 fzero 的收敛判断可能会因为函数值过于陡峭而多花几次迭代无伤大雅。6.2 双曲函数在高阶时的数值行为振型公式里的 cosh 和 sinh 增长非常快。到第 20 阶时 λ ≈ 61cosh(61) 已经是 1e25 量级到第 100 阶时 λ ≈ 313cosh(313) 约 1e135。double 类型没有溢出但 σn 的计算里会隐约出现两个大数相减丢精度的问题表现为高阶振型在自由端附近出现微小波纹噪声。处理办法很简单前 10 阶以内随便用原始公式如果非要算到更高阶σn 直接取 1 的渐近近似或者改用先提出公共大指数因子的写法数值上更稳定。实际上大多数工程场景根本接触不到第 100 阶但如果你做传感器优化布点、高频超声类问题就得把这件事放在心上。6.3 单位制和参数的工程习惯最后一条是最朴实但最烦人的坑单位制。结构振动计算最怕 cm、mm、m 混着用。我推荐全部走 SI长度一律用 m弹性模量用 Pa密度用 kg/m³力用 N时间用 s。频率先算出来是 rad/s要除以 2π 才是 Hz这一步特别容易漏。我自己早年就在这上面吃过亏算出的第一阶频率差了 6.28 倍排查了很久才发现。另一个建议是把几何和材料参数集中在脚本顶部定义成一个结构体p.L 1.0; p.b 0.01; p.h 0.01; p.E 210e9; p.rho 7850;后面所有函数统一从这个结构体取参数改尺寸时就不用满文件搜索替换。这个习惯在参数扫描、优化计算时尤其值钱。最后分享一个我自己的教训。第一次跑这个模型时我把自由端的剪力边界条件符号写反了频率前几阶都很正常唯独某一阶对不上怎么检查代码都看不出问题。后来翻教材回到弯矩、剪力正负号的定义才发现问题出在推导而不是 Matlab 代码。从那以后我养成了新习惯写代码前先手推一遍前两阶的边界矩阵确认符号和方程完全一致再动手写振型函数。写代码本身很快真正花时间的是把物理问题在纸上理清楚。这个习惯推荐给你们。等你把连续体模型跑顺了下一步可以做两件很有价值的事一是把实验频响数据拿来在同一套代码里拟合各阶阻尼比让模型“对标”真实结构二是把梁改成变截面或者加轴向力观察频率和振型怎么变化。那时你就不是单纯抄一个悬臂梁模型而是真正掌握了连续体振动分析的基本功。