悬臂梁连续体振动这是结构动力学里最经典的入门题。不过说它“入门”不代表简单——很多做有限元仿真的人第一次用Matlab算模态就是从悬臂梁开始的。一端固定、一端自由边界条件清晰质量沿长度连续分布刚度沿全长连续变化这样的结构没法用单自由度弹簧-质量块去糊弄只能老老实实把梁当成一个连续体来处理。也正是这个原因悬臂梁成为验证理论公式、数值方法、收敛性的最佳试验田。这篇文章我会把整个研究过程拆开讲清楚连续体模型是怎么建立的、控制方程怎么解、Matlab里解析解和有限元解怎么做、振型怎么可视化、以及实操中我踩过的坑。这套内容适合正在学结构动力学或有限元的学生也适合做机械设备振动分析、模态测试的工程师。如果你已经会用Matlab跑简单的循环和矩阵操作那这篇文章的代码你直接照着敲就能出结果。不需要额外工具箱纯脚本就能跑通。1. 悬臂梁连续体振动从工程直觉到数学描述1.1 悬臂梁到底是谁悬臂梁在工程里太常见了。机械臂的臂身、风力发电机叶片的简化模型、机翼的梁式简化结构、MEMS里的微型谐振梁再往大了说高架桥的悬臂施工段在浇筑时也近似一根超长的悬臂结构。它们共同的特征是一端固定约束另一端自由承受横向弯曲变形。研究悬臂梁振动核心目标不是看它怎么弯而是搞清楚它“自己会怎么抖动”——也就是在没有任何外界持续激励的情况下结构按什么频率、什么形态做自由振动。这个频率叫固有频率这个形态叫振型。工程上最怕的是外界激励频率恰好和固有频率重合引发共振所以算准固有频率是振动分析的第一道门槛。连续体这三个字是重点。梁上每一个截面都是一个质量点而不是某几个离散的质量块。理论上它有无限多个固有频率和振型不像单自由度系统只有一个固有频率。这个“无限多”的特征正好是有限元方法要逼近的核心目标。1.2 连续体模型和单自由度系统的分水岭如果你以前只接触过单自由度振子可能很难理解“一根梁同时有好多个固有频率”到底怎么回事。可以这么类比一把吉他弦拨动时会同时发出基音和高次泛音弦这个连续体天然就有多个振动形态。再想远一点把一根梁想象成很多个小质量块用弹簧串起来的链块数越多系统自由度越多固有频率的数目也越多。当块数趋于无穷弹簧和质量块“糊”成一根连续体所有频率和振型信息都浓缩到了一个偏微分方程里。单自由度系统的运动方程是常微分方程质量、刚度都是标量连续体梁的运动方程是偏微分方程因为位移既是时间的函数又是位置坐标的函数。这意味着自由度是无限的但实际工程中我们只关心前几阶低频模态因为高频模态要么难以激发要么在能耗散下很快衰减。这也是为什么解析解和有限元解都只算前几阶就够了。1.3 频率和振型你要解决的核心问题整个Matlab代码的核心任务就两件事算出固有频率画出振型。固有频率决定了结构在什么激励频率下会共振工程上用来做避振设计振型决定了结构在某阶共振时“长什么样”——哪里变形大、哪里是节点这直接关系到传感器的安装位置、激励点的选择也关系到后续模态叠加法的响应计算。举个例子如果某一阶振型的节点恰好在你想要布置加速度计的位置那这个测点在这阶模态下几乎测不到信号典型的“白装了”。所以代码里我会同时做两套算法一套基于欧拉-伯努利梁理论推导解析解直接用频率方程求根另一套是有限元离散组装刚度矩阵和质量矩阵后解广义特征值问题。两套结果互相验证这才是研究振动模型最靠谱的闭环。2. 理论基础欧拉-伯努利梁方程的建立与求解2.1 控制方程里面每个符号都是物理意义讲到欧拉-伯努利梁公式长这样EI * ∂⁴w(x,t)/∂x⁴ ρA * ∂²w(x,t)/∂t² 0这里w是梁的横向挠度x是沿梁长方向的位置t是时间。EI是抗弯刚度E是弹性模量I是截面惯性矩ρA是线密度ρ是密度A是横截面积。这个方程的物理含义很直白梁在某个截面一侧内力弯矩带来的弹性恢复力和该截面的惯性力相平衡。推导它的时候默认了两个假设平截面假设变形后截面仍保持平面且垂直于中性轴以及小变形假设。换句话说只考虑弯曲变形、忽略剪切变形和转动惯量这也是欧拉-伯努利梁和铁木辛柯梁的根本区别。对一根长细比很大的梁长度远大于截面尺寸欧拉-伯努利模型精度足够。因为剪切变形和转动惯量对低频模态的影响和长细比的平方成反比梁越细长误差越小。这就像细面条甩起来靠弯曲弹性回弹而粗弹簧棒就要考虑剪切效应了。2.2 分离变量把偏微分方程拆成两个常微分方程求解这个偏微分方程最经典的手法叫分离变量法。先把解写成空间函数和时间函数的乘积w(x,t) W(x) * T(t)代入控制方程整理后能得到一个关于时间的方程和一个关于空间的方程。时间部分解出来是一个简谐函数T(t) Acos(ωt) Bsin(ωt)其中ω就是待求的固有频率。空间部分写成更紧凑的形式W(x) - β⁴ * W(x) 0其中β⁴ ω² * ρA / EI。这个参数β非常关键它同时包含了几何参数、材料参数和频率后面求解频率方程时主要就靠它。你看原本一个含有时间项的空间偏微分方程一旦假设是简谐振动就和时间无关了只剩一个四阶常微分方程这就是特征值问题的数学本质。这个常微分方程的通解长这样W(x) C₁cosh(βx) C₂sinh(βx) C₃cos(βx) C₄sin(βx)四个待定系数需要四个边界条件来定。对悬臂梁来说这四个条件清清楚楚写在结构两端。2.3 四个边界条件把频率方程逼出来悬臂梁左端固定位移为零、转角为零右端自由弯矩为零、剪力为零。写成数学语言固定端x0W(0) 0W(0) 0 自由端xLW(L) 0W(L) 0把通解代进去会得到一个四个方程组成的齐次线性方程组。齐次方程要有非零解系数行列式必须等于零。这个“行列式等于零”的条件展开之后就是著名的悬臂梁频率方程cosh(βL) * cos(βL) 1 0这个方程里没有EI也没有ρA只有βL这个整体变量这就是频率参数的无量纲化。前几阶的根是阶次βL11.87510424.69409137.854757410.995541514.137168每多一阶βL近似多π所以高阶根可以用(n-0.5)*π快速估算。有了βL固有频率直接套公式ω (βL)² * sqrt(EI / (ρA * L⁴))这里特别注意量纲。βL是无量纲数但β本身是1/米所以频率公式里必须除一个L⁴才能凑出rad/s的量纲。我第一次手算的时候就栽在这忘了L⁴导致频率大了几个数量级检查半天才发现是单位没对齐。2.4 什么情况不能再用欧拉-伯努利模型很多教程讲到这里就结束了但实际工程里有一个问题必须留意欧拉-伯努利模型什么时候失效。当你处理的梁长细比小于10或者分析的是高频模态时剪切变形和转动惯量的影响会显著变大。这时候要换成铁木辛柯梁模型控制方程从四阶偏微分变成耦合方程组频率不再是βL²这样的简单形式而要和截面形状、剪切系数一起解一个更复杂的超越方程。还有复合材料梁、层合梁需要考虑截面翘曲和横向剪切这些都不能贪图方便直接用欧拉-伯努利。所以在Matlab实现之前先确认你研究的梁确实是细长梁、频率范围在中低阶否则整个模型地基就打歪了。这也是我说“连续体振动模型研究”的第一步不是写代码而是弄清楚你的模型边界在哪。3. Matlab代码实现解析解和有限元双线并行3.1 代码组织一个主脚本加两个函数Matlab写这类程序我习惯不把所有内容堆在一个脚本里太乱。我的组织方式是一个主脚本控制流程两个函数分别做解析解和有限元辅助计算。整个项目三个文件main_beam_vibration.m参数定义、有限元装配、特征值求解、绘图betaL_roots.m用符号变化扫描频率方程返回前N阶βLanalytical_mode.m给定βL返回解析振型曲线这样的好处是逻辑隔离。后面如果你想把梁改成变截面只需要动有限元部分的单元矩阵解析解和绘图代码都不用大改。主脚本头几步长这样clear; clc; close all; % 几何与材料参数 L 1.0; % 梁长单位 m EI 2100; % 抗弯刚度单位 N·m^2 rhoA 7.8; % 线密度单位 kg/m numElem 40; % 有限元网格单元数 nModes 5; % 需要计算的模态数这里的EI和rhoA是我随便给的一组量级合理的值。实际用的时候换成自己结构的等效刚度和等效线密度即可。3.2 解析解fzero在超越方程里挖根频率方程cosh(βL)*cos(βL)10是一个超越方程没法用初等函数求根。Matlab里最朴素的方案是用fzero配合区间扫描。核心思路是用一个很小的步长在βL轴上扫凡是相邻两个点函数值异号说明中间必然跨了一个根直接调fzero锁定它。function betaL betaL_roots(N, bmax) % 返回悬臂梁频率方程的前N个根 betaL(1:N) if nargin 2 bmax (N 1) * pi; % 第N根约在 (N-0.5)pi 附近留出余量 end fun (b) cosh(b) .* cos(b) 1; step 0.01; betaL []; for b step : step : bmax if fun(b - step) * fun(b) 0 % 异号说明区间内有根 r fzero(fun, [b - step, b]); if isempty(betaL) || min(abs(betaL - r)) 1e-6 betaL(end1) r; %#okAGROW end end end betaL betaL(1:N); % 如果长度不足把bmax加大即可 end步长选择有讲究。0.01的步长对前五阶绰绰有余因为相邻根的间距接近π再稀疏的扫描也不会漏根。步长太小会让循环次数过多步长太大则可能漏根万一函数值从正到正还有一个尖峰扫描就错过了。fzero调用时区间两个端点必须函数值异号否则会直接报错所以外面先判断变号是必须的。算出来的betaL直接代入频率公式% 解析固有频率 betaL_val betaL_roots(nModes, 20); omegaExact (betaL_val / L).^2 .* sqrt(EI / rhoA); fExact omegaExact / (2 * pi);这里注意频率单位的换算。omegaExact是圆频率单位rad/s工程上习惯看Hz所以除以2π。我见过好几个新手把这两个弄混画频谱图时横轴频率全差6.28倍。3.3 有限元装配Hermite梁单元的分步拆解解析解只能算等截面、均质、规则边界的悬臂梁。想处理变截面、复杂边界、带附加质量块的梁必须走有限元。悬臂梁的有限元最常用的是二维梁单元Euler-Bernoulli beam element。每个节点有两个自由度横向位移w和转角θ。所以每个单元有4个自由度左端w、左端θ、右端w、右端θ。单元的刚度矩阵和质量矩阵是现成的直接抄标准形式。单元长度记为LeLe L / numElem; % 单元刚度矩阵 Ke EI / Le^3 * [12 6*Le -12 6*Le; 6*Le 4*Le^2 -6*Le 2*Le^2; -12 -6*Le 12 -6*Le; 6*Le 2*Le^2 -6*Le 4*Le^2]; % 单元一致质量矩阵 Me rhoA * Le / 420 * [156 22*Le 54 -13*Le; 22*Le 4*Le^2 13*Le -3*Le^2; 54 13*Le 156 -22*Le; -13*Le -3*Le^2 -22*Le 4*Le^2];这里有两个选择值得解释一下。第一刚度矩阵用的是精确单元这来自于梁控制方程的解而质量矩阵我用了“一致质量矩阵”consistent mass matrix不是把质量直接平均分配到节点上的集中质量矩阵。一致质量矩阵和高阶单元配合收敛速度更快低频结果通常比集中质量更准。缺点是矩阵带宽更大、非对角元多但对小规模问题完全无所谓。第二单元长度Le要动态计算不要直接写死。后面做收敛性分析时numElem一变Le必须跟着变。很多人复制代码时不注意单元数改了但Le还是1结果全错。接下来组装全局矩阵。用一个循环遍历每个单元把单元矩阵按自由度编号累加进全局矩阵nn numElem 1; % 节点数 ndof 2 * nn; % 总自由度 K zeros(ndof, ndof); M zeros(ndof, ndof); for e 1:numElem % 单元两端节点对应的全局自由度 dofs [2*e-1, 2*e, 2*e1, 2*e2]; K(dofs, dofs) K(dofs, dofs) Ke; M(dofs, dofs) M(dofs, dofs) Me; end自由度编号的规则是节点1是w和θ占自由度1和2节点2占3和4以此类推。这种编号方式让单元矩阵在全局矩阵里落在连续的对角块上方便阅读和调试。如果你习惯按“所有w在一起、所有θ在另一块”的方式编号也可以但边界条件处理的时候要跟着变别搞混。3.4 边界条件处理和广义特征值求解悬臂梁左端固定意味着第一个节点的w和θ必须等于0。在有限元里最直接的做法是删除全局矩阵中对应的自由度行和列。因为编号第1、2个自由度恰好是固定端所以代码非常简单% 固定端边界删除自由度1和2 K(1:2, :) []; K(:, 1:2) []; M(1:2, :) []; M(:, 1:2) [];删除自由度法只适合约束不多的小模型好处是矩阵变小、没有引入近似误差。缺点是当有多个约束、或约束自由度在矩阵中间位置时索引容易错。更通用的做法是划行划列法或者乘大数法罚函数法但这些对初学者理解成本更高。悬臂梁这种只有单一固定端的情况删除法是最清晰的。接下来解广义特征值问题[V, D] eig(K, M); ev diag(D); omegaFEM sqrt(max(real(ev), 0)); % 特征值是 ω²开方得到圆频率 % 排序并过滤数值零模态 [omegaFEM, idx] sort(omegaFEM); V V(:, idx); tol max(omegaFEM) * 1e-6; keep omegaFEM tol; omegaFEM omegaFEM(keep); V V(:, keep);这里有一个很容易踩的坑。eig(K,M)返回的特征值对角阵D对角线上的值是λ ω²不是频率本身。所以一定要开方。另外数值求解得到的极小特征值可能是负数比如-1e-12开方会得到虚数。所以先用real取实部、和0取最大值确保不会得到NaN或者虚数再排序、再过滤。过滤的tol设置很有讲究。如果是完全刚体模态特征值几乎是0和第一阶频率相差至少几百倍过滤掉完全无压力。但如果你把tol设得太大可能会误删低阶柔性模态那后面的振型就错位了。我一般用最大频率乘以1e-6作为阈值实测比较稳。4. 振型可视化让模态从矩阵里走出来4.1 解析振型和有限元振型的绘图代码特征值解出来后V矩阵里的每一列就是一个振型向量。不过振型向量里同时包含w和θ直接拿来画图不行得先把每个节点的w分量挑出来。因为我们用的自由度编号是“w、θ交替”而固定端的w被删掉了所以剩下振型向量v的顺序是节点1的w、节点1的θ、节点2的w、节点2的θ……因此挠度节点值用 v(1:2:end) 取奇数位转角节点值用 v(2:2:end) 取偶数位。再把固定端的0补回去wNode [0; v(1:2:end)]; thNode [0; v(2:2:end)]; xNode linspace(0, L, numElem 1);解析振型更方便。前面已经推导过模态函数W(x) cosh(βx) - cos(βx) - σ * (sinh(βx) - sin(βx))其中σ (cosh(βL) cos(βL)) / (sinh(βL) sin(βL))。我把这个封装成函数function [W, xq] analytical_mode(betaL_val, L, npts) % 悬臂梁第n阶解析振型返回最大挠度归一化结果 if nargin 3 npts 200; end b betaL_val / L; xq linspace(0, L, npts); sigma (cosh(betaL_val) cos(betaL_val)) / ... (sinh(betaL_val) sin(betaL_val)); W cosh(b .* xq) - cos(b .* xq) - ... sigma .* (sinh(b .* xq) - sin(b .* xq)); W W / max(abs(W)); % 归一化到最大位移1 end画图的时候用subplot分开展示前几阶对比解析和有限元结果figure(Color, w, Position, [100 100 900 700]); for j 1:3 subplot(2, 2, j); hold on; % 解析解 [W, xq] analytical_mode(betaL_val(j), L, 300); plot(xq, W, k-, LineWidth, 1.5, DisplayName, 解析解); % 有限元节点解 v V(:, j); wNode [0; v(1:2:end)] / max(abs(v(1:2:end))); thNode [0; v(2:2:end)]; plot(xNode, wNode, ro, MarkerSize, 4, DisplayName, FEM节点); grid on; legend; title(sprintf(第 %d 阶模态, f %.3f Hz, j, omegaExact(j) / (2*pi))); xlabel(x (m)); ylabel(归一化振型); end振型归一化很重要。有限元和解析解算出来的振型绝对幅值没有比较意义因为振型本来就是一个相对形状差一个常数倍很正常。所以统一除以最大绝对值大家画在一起才有可比性。4.2 有限元单元内插值节点位移变连续曲线如果你只画节点值单元数少的时候振型会看起来像折线不够平滑。要得到光滑曲线可以在每个单元内部用Hermite插值恢复位移曲线因为梁单元的形函数本来就是三次Hermite插值插出来和理论振型吻合得很好。function [xg, ug] hermite_interp(xNode, wNode, thNode) % 沿梁长度在每个单元内做三次Hermite插值 xg []; ug []; for e 1:numel(xNode) - 1 Le xNode(e1) - xNode(e); xLocal linspace(xNode(e), xNode(e1), 15); t (xLocal - xNode(e)) / Le; N1 1 - 3*t.^2 2*t.^3; N2 Le * (t - 2*t.^2 t.^3); N3 3*t.^2 - 2*t.^3; N4 Le * (-t.^2 t.^3); u N1 * wNode(e) N2 * thNode(e) ... N3 * wNode(e1) N4 * thNode(e1); xg [xg, xLocal]; %#okAGROW ug [ug, u]; %#okAGROW end end注意N2和N4这两个形函数代表转角自由度的影响所以乘以单元长度Le来对齐量纲。如果不乘Le转角对位移的贡献会直接少一个长度量纲插值结果完全不对。插值之后把曲线叠到解析解的图上你会看到两组曲线几乎重叠。这种“眼见为实”的验证比看误差数据更让人安心。4.3 结果对比与网格收敛性以EI2100、rhoA7.8、L1为例解析解的前五阶频率如下表。右侧是我用40个单元跑出来的有限元结果阶次βL解析频率 f (Hz)FEM频率 f (Hz)相对误差11.8751049.199.190.01%24.69409157.5757.580.02%37.854757161.1161.20.06%410.995541316.0316.60.18%514.137168522.0523.90.36%你会发现单元数固定时阶数越高误差越大。这个规律很正常因为高阶振型波长更短需要更细的网格才能分辨。解决方法是要算准第N阶模态网格数量至少要到N的五倍以上。如果是复杂结构或者高阶频率还得加倍。收敛性测试也很简单把numElem分别设成5、10、20、40跑一遍观察第一阶频率的数值变化。第一次跑可能只差了百分之几越往后越小最后稳定在解析解附近。这个过程不仅能验证代码还能让你直观感受到“网格无关解”是怎么回事。5. 实操中反复踩到的问题与排查顺序5.1 特征值出现负值或复数先别慌第一次跑eig(K,M)看到D的对角线上有负值第一反应可能是“完了矩阵组装错了”。不一定。对悬臂梁来说删除固定端自由度后已经不存在刚体模态正常情况特征值全是正数。但数值计算过程中如果某个自由度约束不彻底或者矩阵条件数很差最小特征值可能算出轻微的负值比如-1e-13这种量级。正确做法不是修改程序而是在后处理时过滤掉那些远小于主频率的伪模态。具体就用我前面写的tol过滤法把接近零的特征值直接删掉。如果你发现过滤之后频率依然有几个低于物理常识的值那才需要回头检查边界条件是否删对了行列。另外要留意K矩阵是否对称。eig(K,M)要求K对称、M对称且正定如果组装时单元矩阵累加顺序写错导致K不对称特征值就会大量出现复数。这时候最有效的排查方法是用短暂测试将K打印出来看对角线两边的元素是否严格对称一眼就能看出问题。5.2 有限元频率偏高还是偏低先检查自由度索引有限元频率相对解析解通常偏高把结构变“硬”了但随着网格细化逐步逼近解析解。如果你算出来的第一阶频率严重偏小那多半不是网格问题而是矩阵组装时单元刚度贡献被漏掉了或者边界条件没有正确施加梁在根部的约束等于虚设相当于一根更长的结构在振动频率自然低得离谱。有一个很实用的检查技巧画出第一阶振型看形态对不对。悬臂梁第一阶振型是“固定端位移为0自由端最大曲线没有节点”如果画出来固定端位移不为0别急着调特征值先去看自由度删除的索引顺序。打印K和M减完前后的尺寸确认少了2行2列每个矩阵的总自由度从2*(numElem1)变成了2*numElem这样就对上了。5.3 eig和eigs怎么选小规模模型用eig一次算出所有特征对完全没问题。40个单元的悬臂梁只有80个自由度eig几乎是瞬间完成。但当你把网格加密到上千个单元自由度上到几千几万直接用eig会非常慢甚至有内存溢出风险。这时候改用eigs只求前几阶[V, D] eigs(K, M, nModes, smallestabs);eigs基于迭代法只算你指定数量的特征对效率高得多。但要注意两点第一eigs要求矩阵是稀疏矩阵最好用sparse构造K和M第二不同版本的Matlab对eigs的调用方式有细微差异老版本有的不认识smallestabs这个选项需要换成sm遇到报错先去看版本说明。5.4 振型符号不一致导致对比图“打架”你以为解析解和有限元解是对齐的可画出来发现一个向上一个向下第一反应可能是代码有虫。实际不是振型乘以-1仍然是同一个模态因为求解特征值问题天然带有符号任意性。数值求得的特征向量列向量整体乘一个负号对应的频率完全相同。对比图里要处理这个问题最简单的方法是强制归一化符号一致比如约定振型最大绝对值处的符号为正refIdx find(abs(W) max(abs(W)), 1); if W(refIdx) 0 W -W; end有限元那边同样处理。这样画出来两组曲线必定方向一致。这个小细节我第一次对比时没注意白折腾了一晚上最后发现只是符号问题。实际操作中的几点延伸建议做完这套悬臂梁连续体振动模型最基本的频率和振型计算就跑通了。我个人在实际操作中最受益的一点是永远先跑解析解验证数值框架。哪怕你后面要做变截面梁、带集中质量的梁也先留一个最简单的悬臂梁案例当基准任何改动之后都能立刻发现是哪里出了问题。这套代码往后扩展的空间很大。把单元质量矩阵换成集中质量矩阵可以对比两种离散方式的差异在方程里加阻尼项就能做受迫振动的稳态响应把EI和rhoA从常数改成关于x的函数就变成变截面梁。还有一条路是给梁末端加一个质量块或弹簧地基边界条件从固定-自由变成固定-弹性约束这就在朝工程实际结构走了。悬臂梁就像结构动力学里的一把尺子方法对了它永远是平的、准的。把这把尺子握牢后面做板、壳、三维结构底层逻辑都是一样的。