做动力学分析的人对“悬臂梁冲击”这个题应该都不陌生。机械臂在意外碰撞时根部应力怎么算无人机起落架着陆瞬间的冲击响应怎么估分离机构解锁后的瞬态位移怎么看这些工程问题落到最简单的模型上往往就是一个悬臂梁自由端受冲击载荷。问题看起来边界条件简单但真要把解析解和有限元解都跑通中间牵扯到模态叠加、振型正交、数值积分时间步长、单元刚度矩阵组装这些环节每一步都有坑。我把自己做过的一个小项目完整梳理一遍用解析方法也就是模态叠加法求解悬臂梁自由端受半正弦冲击载荷时的动力响应再用自编Matlab有限元程序做同一个算例对比两条位移时程曲线。文章会给出完整的Matlab代码解释每个关键函数的物理含义和数值实现细节。适合正在学结构动力学、有限元入门或者需要做瞬态响应分析的同学参考。代码不依赖任何商业工具箱用MATLAB基础包就能跑起来。1. 悬臂梁冲击问题的工程背景与求解思路1.1 冲击问题为什么比静力问题麻烦静力问题求解的是位移和应力在某个固定载荷下的平衡状态方程是 K u F解一次就够了。但冲击问题不一样载荷随时间剧烈变化而且作用时间往往很短结构内部的惯性力不能忽略。比如一个质量块以一定速度撞到悬臂梁自由端力的峰值可能很大但持续时间只有几十毫秒结构根本来不及达到静力平衡状态。这时候控制方程变成了 M ü C u̇ K u F(t)是一组关于时间的常微分方程组。很多人刚接触动力学时容易犯一个错误把峰值载荷直接当静载荷加到结构上然后校核最大应力。这种做法在载荷变化非常缓慢的时候是保守的但冲击载荷频谱很宽会激发结构的高阶模态导致实际响应可能比静力解大好几倍。所以冲击分析必须走时程分析路线把每一时刻的位移、速度、加速度都算出来再从中提取峰值响应。1.2 为什么拿悬臂梁做基准题悬臂梁是整个结构动力学里面最经典的连续体模型之一。固定端位移和转角为零自由端自由边界条件写起来干净。更关键的是欧拉-伯努利梁的振型和固有频率有解析表达式可以拿来当“标准答案”验证有限元程序的正确性。工程上悬臂梁结构也到处都是支架、悬臂吊、天线桅杆、机械臂大臂甚至电路板上的引脚都可以简化成悬臂梁。所以拿它做冲击响应基准题既有理论代表性又有工程适用性。如果连悬臂梁的动力学都算不对那算复杂结构的结果可信度就很低。1.3 解析解与有限元解怎么分工这个项目里我同时用了两种方法目的不是比谁更厉害而是互相校验。解析解基于模态叠加法思路是把连续体的响应分解成一系列固有振型的线性组合每个振型对应一个单自由度方程。只要材料线弹性、几何小变形这个方法是精确的误差主要来自模态截断和数值积分。有限元法则先把连续体离散成若干个梁单元每个节点有挠度和转角两个自由度然后用Newmark-β方法做时间积分。有限元的优点是可以推广到变截面、复杂边界、非线性和多体结构缺点是需要仔细处理单元数量和时间步长。把两者摆在一起对比既能验证有限元程序写没写对又能反过来评估模态叠加法在截断高阶模态后的误差到底有多大。2. 解析解推导从梁振动方程到模态叠加2.1 控制方程和分离变量欧拉-伯努利梁的自由振动控制方程是ρA ∂²w/∂t² EI ∂⁴w/∂x⁴ f(x, t)其中 w 是梁的横向挠度EI 是抗弯刚度ρA 是单位长度质量f(x,t) 是分布力。这个方程的物理含义很直白惯性力加上弹性恢复力等于外力。对于自由振动设 f(x,t)0令 w(x,t)φ(x)q(t)代入后可以把时间和空间变量分开得到两个方程。空间部分满足d⁴φ/dx⁴ - β⁴ φ 0其中 β⁴ ρA ω² / EI。这个四阶常微分方程的通解可以写成三角函数和双曲函数的组合具体形式由边界条件决定。2.2 悬臂梁的固有频率与振型函数悬臂梁的边界条件是固定端 x0 处位移和转角为零自由端 xL 处弯矩和剪力为零。把这四个边界条件代进去经过一番推导会得到一个关于 β 的特征方程cosh(βL) · cos(βL) -1这个方程没有闭式解只能用数值方法求根。前五阶 βL 的值大概是 1.8751、4.6941、7.8548、10.9955、14.1372。从第二阶开始每一阶都比上一阶略小于 (2n-1)π/2这个规律在后面写求根代码时非常有用。振型函数可以写成φ(x) cosh(βx) - cos(βx) α [sinh(βx) - sin(βx)]系数 α 由自由端弯矩为零的条件确定α -(cosh(βL) cos(βL)) / (sinh(βL) sin(βL))有了振型函数固有频率就可以通过 β 算出来ω β² √(EI / (ρA))注意频率和 β 是平方关系所以高阶模态的频率上升很快。这也是冲击响应分析里高频模态不可忽略的原因——冲击载荷虽然持续时间短但它能在瞬间把能量注入到高阶模态里。2.3 模态叠加法求解冲击响应有了振型和固有频率就可以把实际受迫振动的位移展开成振型的叠加w(x,t) Σ φₙ(x) qₙ(t)代入受迫振动方程利用振型关于质量矩阵和刚度矩阵的正交性每个模态坐标 qₙ(t) 满足一个独立的单自由度方程qₙ 2ζₙωₙ qₙ ωₙ² qₙ Fₙ(t) / Mₙ其中Mₙ ∫ ρA φₙ² dx 是第 n 阶模态的广义质量Fₙ(t) ∫ φₙ(x) f(x,t) dx 是广义力ζₙ 是模态阻尼比对于自由端集中力 F(t)广义力简化为 Fₙ(t) F(t) · φₙ(L)因为力只作用在 xL 这个点上。我这次算例用的是无阻尼模型ζₙ0那么杜哈梅积分给出qₙ(t) (1 / (Mₙ ωₙ)) ∫₀ᵗ Fₙ(τ) sin[ωₙ(t-τ)] dτ这个积分可以用数值积分来求也可以用解析方法求。我为了代码通用性直接用 trapz 做数值积分虽然慢一点但换载荷函数时不用改公式。2.4 解析解的两个“隐形假设”模态叠加法虽然精确但有两个前提条件很容易被忽略。第一材料必须线弹性变形必须是小变形。冲击力大到引起塑性变形或者大挠度时模态叠加法不再适用因为振型本身已经变了。第二欧拉-伯努利梁理论忽略了剪切变形和转动惯量所以长细比越大精度越高。对于短粗梁比如长细比小于 10 的结构应该改用 Timoshenko 梁理论。我算例里梁长 1 米截面 50mm × 50mm长细比 20满足欧拉-伯努利梁的适用条件。代码实现时还有一个小坑振型的归一化方式会影响 φₙ(L) 和 Mₙ但最终响应 w 不变因为两个量同时缩放。这相当于数学上的“比例不变性”写程序时不用担心归一化方式选得不对。3. 有限元程序架构从单元矩阵到Newmark积分3.1 为什么自己写Matlab有限元而不是直接上商业软件有人可能会问现在 Ansys、Abaqus 这么好用为什么还要自己用 Matlab 写有限元程序我的看法是商业软件适合算大型复杂模型但不适合用来理解算法。对于悬臂梁冲击这个尺度的题商业软件的建模、网格划分、求解设置反而更繁琐。自己写程序最大的好处是每一步都透明单元矩阵长什么样、怎么组装、边界条件怎么施加、时间积分怎么推进全部可以打印出来逐行检查。这也是有限元教学里一直保留“手写程序”这个环节的原因。另一个实际原因是参数化研究方便。我想看单元数量从 5 个变成 40 个时结果怎么变化商业软件里要么改网格重新求解要么写脚本而小程序里一个 for 循环就搞定了。3.2 欧拉-伯努利梁单元的刚度矩阵与质量矩阵每个梁单元有两个节点每个节点有挠度 w 和转角 θ 两个自由度。单元长度 le是弹性模量 E截面惯性矩 I单元刚度矩阵是经典的四阶方阵k EI/le³ * [ 12 6le -12 6le; 6le 4le² -6le 2le²; -12 -6le 12 -6le; 6le 2le² -6le 4le² ]这个矩阵大家可能见过无数次但要注意坐标约定自由度顺序是 [w1, θ1, w2, θ2]θ 以逆时针为正。组装时一定要保证局部自由度和全局自由度一一对应否则矩阵位置放错结果会非常离谱。质量矩阵我选了一致质量矩阵不是集中质量矩阵m ρA·le/420 * [ 156 22le 54 -13le; 22le 4le² 13le -3le²; 54 13le 156 -22le; -13le -3le² -22le 4le² ]一致质量矩阵由单元形函数积分得到能更准确地描述质量在单元内的连续分布。集中质量矩阵是假设质量集中在节点上求固有频率会偏低而且时程分析中对高频振型的精度明显不如一致质量矩阵。对于冲击这种宽带激励问题建议优先用一致质量矩阵。3.3 系统组装与边界条件处理我采用的自由度编号规则是第 i 个节点的全局自由度为 2i-1挠度和 2i转角自由度总数为 2×节点数。组装循环里单元 e 连接节点 e 和 e1对应局部自由度 [w_e, θ_e, w_{e1}, θ_{e1}]全局编号是 [2e-1, 2e, 2e1, 2e2]然后把单元矩阵累加到总矩阵对应位置。边界条件的处理用的是“划行划列法”思路很直接把固定端节点对应的全局自由度从求解集合里移除。例如第一个节点的自由度 1 和 2 被约束那我只需要求解自由度 3 到最后的子矩阵。这种方法的好处是缩小了求解规模坏处是如果你后面需要输出所有节点的位移还需要把约束自由度补零放回去。代码里我用 dof_free 3:2*nNode 提取自由度集并把自由端载荷加到对应位置。3.4 Newmark-β方法时间积分动力学方程 M ü C u̇ K u F(t) 是二阶常微分方程组需要时间积分方法逐步推进。我用了 Newmark-β 方法这是一种广泛使用的隐式时间积分法核心是两个参数 β 和 γ 控制精度和稳定性。对于线性问题取 β0.25、γ0.5 时就是平均加速度法它的特点是无条件稳定无论时间步长取多大结果不会因为数值原因发散。但无条件稳定不代表时间步可以随便取因为时间步太大时高频响应会被严重过滤掉导致峰值偏小。Newmark 方法的实现可以概括为三步计算等效刚度矩阵 K_hat K a0·M a1·C。对每个时间步基于当前位移、速度、加速度计算等效载荷向量。解线性方程组 K_hat · u_new F_hat然后更新速度和加速度。具体系数我在后续代码里给出这里先记住一个原则初始加速度必须用 M(F(0) - K·u0 - C·v0) 精确求出不能直接设为 0否则从第一步开始就会引入误差。4. Matlab代码实现分模块讲解4.1 主脚本参数设置与整体流程主脚本负责定义所有参数然后调用有限元程序和模态叠加程序。我用半正弦脉冲模拟冲击载荷F(t) F0 · sin(πt/Td)当 0 ≤ t ≤ Td否则为 0选用半正弦而不是矩形脉冲是因为矩形脉冲起点和终点是阶跃突变频谱尾部衰减慢需要很多模态才能准确捕捉而半正弦脉冲的频谱衰减相对快一些对比时收敛性更好。% cantilever_impact_main.m clc; clear; close all; % 几何与材料参数 L 1.0; % 梁长 [m] b 0.05; h 0.05; % 矩形截面宽高 [m] A b*h; % 截面面积 [m^2] I b*h^3/12; % 截面惯性矩 [m^4] E 210e9; % 弹性模量 [Pa] rho 7850; % 密度 [kg/m^3] % 冲击载荷参数 F0 1000; % 冲击力峰值 [N] Td 0.01; % 冲击持续时间 [s] t_end 0.10; % 总计算时间 [s] % 数值离散参数 nElem 20; % 单元数量 dt 1e-4; % 时间步长 [s] nSteps round(t_end/dt); % 载荷函数句柄 force_func (t) F0 * sin(pi * min(t, Td) / Td) .* (t Td); % 有限元求解 [t_fem, w_fem, theta_fem] fem_cantilever_impact(L, A, I, E, rho, nElem, dt, nSteps, force_func); % 解析解模态叠加法 Nmode 8; % 模态截断数 [t_ana, w_ana] modal_impact_solution(L, A, I, E, rho, Nmode, dt, nSteps, force_func); % 对比自由端挠度时程 figure(Color,white); plot(t_fem, w_fem(end,:), b-, LineWidth, 1.5); hold on; plot(t_ana, w_ana, r--, LineWidth, 1.5); xlabel(时间 t [s]); ylabel(自由端挠度 w [m]); legend(有限元解 (20单元), 解析解 (模态叠加)); title(悬臂梁自由端受冲击载荷的响应对比); grid on;4.2 单元矩阵函数与系统组装beam_k_e 函数根据单元长度返回刚度矩阵beam_m_e 返回一致质量矩阵。这里注意质量矩阵需要用到单位长度质量 rhoA我把它作为参数传入。function ke beam_k_e(E, I, le) ke E*I / 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]; end function me beam_m_e(rhoA, le) 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]; end组装和约束都在 fem_cantilever_impact 函数里完成。函数首先初始化全零的总体矩阵然后循环所有单元把单元矩阵累加到对应自由度位置。最后提取自由度为 3 到 2*nNode 的子矩阵也就是去掉固定端的挠度和转角自由度。function [t, w, theta] fem_cantilever_impact(L, A, I, E, rho, nElem, dt, nSteps, force_func) nNode nElem 1; ndof 2*nNode; le L / nElem; K zeros(ndof, ndof); M zeros(ndof, ndof); for e 1:nElem idx [2*e-1, 2*e, 2*e1, 2*e2]; K(idx, idx) K(idx, idx) beam_k_e(E, I, le); M(idx, idx) M(idx, idx) beam_m_e(rho*A, le); end % 固定端约束节点1的挠度和转角均设为0 dof_free 3:ndof; Kf K(dof_free, dof_free); Mf M(dof_free, dof_free); Cf zeros(size(Kf)); % 无阻尼 % 载荷作用位置自由端挠度自由度全局编号 ndof-1 fdof find(dof_free ndof-1); u0 zeros(length(dof_free), 1); v0 zeros(length(dof_free), 1); % Newmark 时间积分 [u, v, a, t] newmark_beta(Mf, Cf, Kf, u0, v0, dt, nSteps, ... (tt) load_vector(tt, fdof, force_func, length(dof_free))); % 还原完整节点位移 w zeros(nNode, nSteps1); theta zeros(nNode, nSteps1); w(2:end, :) u(1:2:end, :); % 节点2开始为自由挠度 theta(2:end, :) u(2:2:end, :); end function F load_vector(t, fdof, force_func, nfree) F zeros(nfree, 1); F(fdof) force_func(t); end组装这里有个值得注意的细节循环变量 idx 用的是节点自由度映射而不是局部单元自由度。很多初学者容易把 idx 写成 [1 2 3 4]那样所有单元矩阵都堆到同一个位置结果肯定是错的。检查组装是否正确的一个简单方法是打印组装后 K 矩阵的带宽如果带宽不对说明自由度编号映射出了问题。4.3 Newmark-β求解器代码newmark_beta 函数是动力学求解的核心。它接收质量、阻尼、刚度矩阵和初始条件以及载荷函数返回每一步的位移、速度和加速度。function [u, v, a, t] newmark_beta(M, C, K, u0, v0, dt, nSteps, Ffun) beta 0.25; gamma 0.5; % 平均加速度法 ndof length(u0); u zeros(ndof, nSteps1); v zeros(ndof, nSteps1); a zeros(ndof, nSteps1); t (0:nSteps) * dt; % 初始加速度 a(:,1) M \ (Ffun(0) - K*u0 - C*v0); u(:,1) u0; v(:,1) v0; % 预处理系数 a0 1 / (beta*dt^2); a1 gamma / (beta*dt); a2 1 / (beta*dt); a3 1/(2*beta) - 1; a4 gamma/beta - 1; a5 dt/2 * (gamma/beta - 2); Keff K a0*M a1*C; for i 1:nSteps Fhat Ffun(t(i1)) ... M*(a0*u(:,i) a2*v(:,i) a3*a(:,i)) ... C*(a1*u(:,i) a4*v(:,i) a5*a(:,i)); u(:,i1) Keff \ Fhat; v(:,i1) a1*(u(:,i1)-u(:,i)) - a4*v(:,i) - a5*a(:,i); a(:,i1) a0*(u(:,i1)-u(:,i)) - a2*v(:,i) - a3*a(:,i); end end这个求解器有几个关键点第一等效刚度矩阵 Keff 只需要计算一次不用在每个时间步重新组装这是隐式方法的最大优势。第二载荷函数 Ffun 在每一步只调用一次但要注意它返回的是列向量。第三Newmark 方法是无条件稳定的所以即使时间步大也不会发散但时间步太大会导致响应峰值偏小这个在后面的收敛性分析里会说到。4.4 模态叠加解析解代码模态叠加部分需要求解固有频率和振型。频率根用到 fzero这是 MATLAB 基础包里的函数不需要额外工具箱。求根初始值我选在 (2n-1)π/2 附近因为当 βL 增大时悬臂梁频率根越来越接近奇数的 π/2 倍。function [t, w] modal_impact_solution(L, A, I, E, rho, Nmode, dt, nSteps, force_func) t (0:nSteps) * dt; x_tip L; w zeros(1, nSteps1); % 求各阶特征根 betaL betaL zeros(1, Nmode