企业官网建设流程全解析

做动力学分析的人对“悬臂梁冲击”这个题应该都不陌生。机械臂在意外碰撞时根部应力怎么算，无人机起落架着陆瞬间的冲击响应怎么估，分离机构解锁后的瞬态位移怎么看，这些工程问题落到最简单的模型上，往往就是一个悬臂梁自由端受冲击载荷。问题看起来边界条件简单，但真要把解析解和有限元解都跑通，中间牵扯到模态叠加、振型正交、数值积分时间步长、单元刚度矩阵组装这些环节，每一步都有坑。

我把自己做过的一个小项目完整梳理一遍：用解析方法，也就是模态叠加法，求解悬臂梁自由端受半正弦冲击载荷时的动力响应；再用自编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 悬臂梁的固有频率与振型函数

悬臂梁的边界条件是固定端 x=0 处位移和转角为零，自由端 x=L 处弯矩和剪力为零。把这四个边界条件代进去，经过一番推导，会得到一个关于 β 的特征方程：

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)，因为力只作用在 x=L 这个点上。

我这次算例用的是无阻尼模型，ζₙ=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 和 e+1，对应局部自由度 [w_e, θ_e, w_{e+1}, θ_{e+1}]，全局编号是 [2e-1, 2e, 2e+1, 2e+2]，然后把单元矩阵累加到总矩阵对应位置。

边界条件的处理用的是“划行划列法”，思路很直接：把固定端节点对应的全局自由度从求解集合里移除。例如第一个节点的自由度 1 和 2 被约束，那我只需要求解自由度 3 到最后的子矩阵。这种方法的好处是缩小了求解规模，坏处是如果你后面需要输出所有节点的位移，还需要把约束自由度补零放回去。代码里我用 dof_free = 3:2*nNode 提取自由度集，并把自由端载荷加到对应位置。

3.4 Newmark-β方法时间积分

动力学方程 M ü + C u̇ + K u = F(t) 是二阶常微分方程组，需要时间积分方法逐步推进。我用了 Newmark-β 方法，这是一种广泛使用的隐式时间积分法，核心是两个参数 β 和 γ 控制精度和稳定性。

对于线性问题，取 β=0.25、γ=0.5 时就是平均加速度法，它的特点是无条件稳定：无论时间步长取多大，结果不会因为数值原因发散。但无条件稳定不代表时间步可以随便取，因为时间步太大时高频响应会被严重过滤掉，导致峰值偏小。

Newmark 方法的实现可以概括为三步：

1. 计算等效刚度矩阵 K_hat = K + a0·M + a1·C。
2. 对每个时间步，基于当前位移、速度、加速度，计算等效载荷向量。
3. 解线性方程组 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*e+1, 2*e+2]; 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, nSteps+1); theta = zeros(nNode, nSteps+1); 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, nSteps+1); v = zeros(ndof, nSteps+1); a = zeros(ndof, nSteps+1); 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(i+1)) ... + M*(a0*u(:,i) + a2*v(:,i) + a3*a(:,i)) ... + C*(a1*u(:,i) + a4*v(:,i) + a5*a(:,i)); u(:,i+1) = Keff \ Fhat; v(:,i+1) = a1*(u(:,i+1)-u(:,i)) - a4*v(:,i) - a5*a(:,i); a(:,i+1) = a0*(u(:,i+1)-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, nSteps+1); % 求各阶特征根 betaL betaL = zeros(1, Nmode