简介面向结构动力学与有限元分析学习者的 MATLAB 板结构有限元模态分析程序以单个脚本集中实现建模、求解与可视化适合需要理解自然频率、振型与共振判别方法的工程师和高校学生参考。程序中预期涵盖板结构建模、质量/刚度矩阵组装、特征值问题求解以及振型可视化等关键环节能够帮助读者快速掌握用有限元方法完成模态计算的基本流程。压缩包内仅含 1 个 m 脚本文件整体大小约 1KB体量虽小但功能指向明确适合用于教学演示或二次开发基础。目前已有 174 人学习下载尤其适合刚接触 MATLAB 有限元编程、希望从简单板模型入手了解模态分析原理的读者。通过阅读该脚本可直观学习节点坐标定义、单元矩阵组装与 eig 类函数求解的实现思路并据此迁移到更复杂的结构动力分析场景。1. 板结构模态分析从四节点四边形单元的MATLAB源码说起做结构动力学的人都有这种体会用ANSYS或ABAQUS做单次模态分析非常舒服但一旦需要扫描几十组板厚、长宽比或者在优化算法里反复调用频率计算商业软件的建模流程就开始拖后腿。这套源码里的ModalAnalysis.m把板结构模态分析直接搬到MATLAB里用四节点四边形单元把矩形薄板离散成平面应力问题自己组装刚度矩阵和质量矩阵再用广义特征值求解器提取固有频率和振型。它适合正在啃有限元编程细节的学生也适合做参数预研时不想频繁切换工具的工程师。程序结构不大但能把有限元模态分析的完整链路——网格生成、单元积分、矩阵组装、约束处理、特征值求解、振型可视化——全部串起来。2. 有限元模态分析的数学基础与程序架构2.1 从弹性力学方程到广义特征值问题对一块厚度远小于面内尺寸的平板当激励作用方向位于板平面内时可以近似为平面应力问题。结构离散化之后自由振动的动力学方程写成M * u_ddot K * u 0其中 M 是质量矩阵K 是刚度矩阵u 是节点位移向量。假设自由振动解为谐波形式 u φ * exp(iωt)代入后消去时间因子得到K φ ω² M φ这是典型的广义特征值问题。ω² 对应结构固有频率的平方φ 是振型向量。模态分析的核心就是求解这个广义特征值问题的前若干阶最小特征对。和标准特征值问题 A x λ x 不同这里的 M 不一定正定尤其当约束不充分时还可能出现零特征值所以不能直接对 K 求逆再去算普通特征问题。MATLAB 的 eigs 函数专门处理广义问题底层使用 ARPACK可以只求指定数量的最小特征值计算效率远高于把矩阵完全分解后取全部特征对。2.2 ModalAnalysis.m 主流程拆解这个源码只有一个文件但流程上很清晰。拿到手之后我习惯把它拆成下面几步来理解输入几何尺寸、材料参数、网格密度和边界条件生成节点坐标矩阵 nodes 和单元连接矩阵 elements对每个单元计算刚度矩阵和质量矩阵把单元矩阵组装到全局 K 和 M根据约束自由度对全局矩阵做缩减用 eigs 求解前 N 阶特征对把特征向量映射回完整自由度空间并绘图主函数骨架大致如下function modal_plate() % 几何与材料参数 Lx 1.0; Ly 0.5; t 0.01; E 210e9; nu 0.3; rho 7850; nx 8; ny 4; nmodes 6; % 1) 网格生成 [nodes, elements] mesh_plate(Lx, Ly, nx, ny); % 2) 单元矩阵组装 [K, M] assemble_system(nodes, elements, E, nu, rho, t); % 3) 边界自由度约束 fixed_dofs find_fixed_dofs(nodes, Lx); [K, M] apply_bc(K, M, fixed_dofs); % 4) 特征值求解 [V, D] eigs(K, M, nmodes, smallestabs); freqs sqrt(diag(D)) / (2*pi); % 5) 后处理 plot_modes(nodes, elements, V, freqs, nx, ny); end这里的 mesh_plate、assemble_system 等函数都是源码内的自有实现不依赖 MATLAB PDE 工具箱。Lx、Ly 控制板的平面尺寸nx、ny 控制长宽方向的网格划分数nmodes 是需要提取的模态阶数。eigs 的 smallestabs 参数表示求模最小的特征值对应最低频的几个模态物理上最关心的正是这部分。2.3 单元类型选择为什么板结构用四节点四边形单元平面应力问题里最常见的两种单元是三角形常应变单元CST和四节点双线性四边形单元Q4。CST 单元位移场是线性的单元内应变恒定精度偏低Q4 单元位移场是双线性的应变沿单元内线性变化在同样的网格密度下精度明显更好。对模态分析来说刚度矩阵的准确性直接影响频率结果Q4 是性价比最高的选择。特性CST 三角单元Q4 四边形单元位移场阶次线性双线性单元内应变常数线性变化收敛速率较慢较快数值积分1 点积分2×2 高斯积分编程复杂度低中等模态分析适用性教学入门工程预研Q4 单元也有需要注意的地方。如果网格划分时四边形被压成接近三角形的退化形状雅可比行列式可能在积分点趋近于零甚至变负导致单元刚度矩阵不正定特征值求解时出现虚数频率。后面第 5 章会专门说怎么排查这类问题。3. 核心实现刚度矩阵、质量矩阵组装与边界条件处理3.1 节点坐标与单元连接关系构建网格生成是第一个容易写错的地方。一个稳妥的做法是先把矩形板均匀划分成 nx × ny 个格子节点按行扫描编号。例如 1m × 0.5m 的板nx4、ny2节点总数是 (41)×(21)15单元总数是 8。下面的辅助函数生成节点坐标和单元连接function [nodes, elements] mesh_plate(Lx, Ly, nx, ny) npx nx 1; npy ny 1; nodes zeros(npx * npy, 2); for j 1:npy for i 1:npx nodes((j-1)*npx i, :) [(i-1)*Lx/nx, (j-1)*Ly/ny]; end end elements zeros(nx * ny, 4); for j 1:ny for i 1:nx n1 (j-1)*npx i; n2 n1 1; n3 n2 npx; n4 n1 npx; elements((j-1)*nx i, :) [n1, n2, n3, n4]; end end end参数 Lx、Ly 是板的长宽nx、ny 是长宽方向的网格数。单元连接矩阵每一行按逆时针顺序给出四个节点编号这是保证雅可比行列式一致为正的前提。MATLAB 的列优先习惯很容易让节点编号方向出错一旦单元顺序变成顺时针后面计算应变矩阵时就会出负号。我通常在生成网格后先 plot 一遍用 patch 画出单元边界确认连接没有交叉。3.2 平面应力单元刚度矩阵的显式计算Q4 单元刚度矩阵的推导核心在等参变换。局部坐标 ξ∈[-1,1]、η∈[-1,1]四个形函数为N1 (1-ξ)(1-η)/4N2 (1ξ)(1-η)/4 N3 (1ξ)(1η)/4N4 (1-ξ)(1η)/4。通过雅可比矩阵把形函数对 ξ、η 的导数转换到物理坐标系得到应变矩阵 B。平面应力弹性矩阵 D 为function [Ke, Me] q4_stiffness(xy, E, nu, rho, t) % xy: 4x2 节点坐标矩阵 % Ke: 8x8 单元刚度矩阵 % Me: 8x8 单元一致质量矩阵 D E / (1 - nu^2) * [1 nu 0; nu 1 0; 0 0 (1-nu)/2]; gps [-0.577350269189626, 0.577350269189626]; Ke zeros(8, 8); Me zeros(8, 8); for i 1:2 xi gps(i); for j 1:2 eta gps(j); [N, dNdx, dNdy, detJ] q4_shape(xi, eta, xy); B zeros(3, 8); for a 1:4 B(1, 2*a-1) dNdx(a); B(2, 2*a) dNdy(a); B(3, 2*a-1) dNdy(a); B(3, 2*a) dNdx(a); end Ke Ke (B * D * B) * detJ * t; Nmat zeros(2, 8); Nmat(1, 1:2:end) N; Nmat(2, 2:2:end) N; Me Me (Nmat * Nmat) * rho * t * detJ; end end end这里使用 2×2 高斯积分点两个积分点坐标是 ±1/√3权重默认都是 1。B 矩阵的组装需要细心第 a 个节点的 x 方向自由度编号是 2a-1y 方向自由度编号是 2a。B 的第三行对应剪切应变 γxy需要把 dNdx 和 dNdy 交叉填入。初学者最容易在这里把剪切项写错导致刚度矩阵误差最后算出的频率偏大或偏小到完全不可信。3.3 一致质量矩阵与集中质量矩阵的取舍质量矩阵有两类常用形式。一致质量矩阵从形函数完整积分而来Me ρ t ∫ N^T N dA它的质量分布和刚度矩阵一致频率结果更接近连续体模型尤其是高阶模态。集中质量矩阵则把单元总质量平均分配到四个节点只保留对角线元素。两者的差异在下表里可以很清晰地看到类型矩阵形态存储开销频率精度适用场景一致质量矩阵8×8 满块含非对角耦合高更准模态分析、瞬态响应集中质量矩阵对角阵低系统偏低显式动力学、大规模模型在编写这个源码时我优先选择一致质量矩阵。板类结构的弯曲模态对惯性分布极其敏感集中质量强行把质量压到节点上会低估低频响应。不过要注意一致质量矩阵在求解时如果出现负特征值往往意味着 M 不是正定的这时需要检查单元畸变和材料参数。作为工程习惯我会把集中质量矩阵和一质量矩阵的结果都跑一遍理论上真实频率应该落在两个结果之间如果偏差超过 5%说明网格太粗。3.4 边界条件的施加抑制刚体模态未施加约束的弹性体存在刚体模态特征值为零eigs 求解时会不稳定甚至报错。最常见的办法是把约束自由度从系统中删掉function [Kr, Mr, free_dofs] apply_bc(K, M, fixed_dofs) free_dofs setdiff(1:size(K,1), fixed_dofs); Kr K(free_dofs, free_dofs); Mr M(free_dofs, free_dofs); end这里 fixed_dofs 是自由度编号列表不是节点编号。第 n 个节点的 x 方向自由度是 2n-1y 方向自由度是 2n。固支边要把边上所有节点的两个自由度都加进去。用删行删列的方式比置大数法更干净不会引入高频伪模态。如果只约束 x 方向板仍然可以整体沿 y 方向滑动MATLAB 的 eigs 会返回一个大约 1e-10 量级的伪零频率。判断方法是画出振型如果所有节点位移方向一致、大小基本相同那就是刚体平动不是真实模态。4. 实战运行源码、求解特征值并绘制振型4.1 设置材料参数与网格密度实际使用这份源码时第一步应该用一块简支方板做测试目的是确认整个流程闭环。我建议的初始参数如下Lx 1.0; Ly 0.5; t 0.01; E 210e9; nu 0.3; rho 7850; nx 8; ny 4; [nodes, elements] mesh_plate(Lx, Ly, nx, ny); [K, M] assemble_system(nodes, elements, E, nu, rho, t);这里 t 的单位是米E 单位是帕斯卡密度单位是 kg/m³最终频率单位是 Hz。如果模型长度用毫米E 得用 N/mm²密度得用 kg/mm³否则频率会差出 1000 倍。这个单位换算坑非常隐蔽而且往往只在和后处理结果对比时才会暴露。为了避免它我习惯在脚本开头把单位和换算系数做成注释每次运行前检查一遍。4.2 用 eigs 求解前六阶模态完成矩阵组装和约束施加后调用 eigs 求解fixed_dofs find_fixed_dofs(nodes, Lx, Ly); % 自定义边界条件 [Kf, Mf, free_dofs] apply_bc(K, M, fixed_dofs); [V, D] eigs(Kf, Mf, 6, smallestabs); omega sqrt(diag(D)); freqs omega / (2*pi); disp(前六阶固有频率(Hz):); disp(freqs); % 还原到完整自由度空间方便画图 V_full zeros(size(K, 1), 6); V_full(free_dofs, :) V; V_full(fixed_dofs, :) 0;参数含义Kf、Mf 是约束后的刚度、质量矩阵6 是求解阶数smallestabs 表示求模最小的特征值。V 的每一列是对应约束后自由度的特征向量直接用来做云图会少一段索引。所以用 V_full 把完整位移场拼出来被约束的自由度位移强制置零。eigs 依赖于 ARPACK 迭代如果矩阵维度几百到几千计算很快但若出现 NaN 或 Inf优先检查 M 是否有零对角线元素比如某个自由度的节点没有附着任何单元时就会出现奇异。4.3 绘制振型云图直观看模态形状有了完整振型向量下一步是画图。板的面内位移只有 Ux、Uy但为了可视化模态特征我通常把其中一个位移分量作为高度场function plot_modes(nodes, elements, V_full, freqs, nx, ny) npx nx 1; nv size(V_full, 2); for k 1:nv Ux V_full(1:2:end, k); Uy V_full(2:2:end, k); Uz reshape(Ux, npx, ny1) * 30; X reshape(nodes(:,1), npx, ny1); Y reshape(nodes(:,2), npx, ny1); figure; surf(X, Y, Uz, Uz); title(sprintf(Mode %d, f %.2f Hz, k, freqs(k))); view(3); colorbar; end end这里的放大系数 30 是经验值目的是让变形在图上肉眼可见。需要注意的是这个源码把板当作平面应力问题求解的是面内振动模态如果你想要的是薄板弯曲振动也就是垂直于板面的上下起伏那需要换成 Mindlin 板单元或 Kirchhoff 板单元。面内模型作为有限元算子验证是够用的但拿它和工程中的弯曲模态测试值对比会量级都不对这一点必须在用之前分清楚。4.4 网格收敛性检查模态频率对网格密度敏感尤其是扭转模态和高阶模态。最常见的验证方法是固定材料和边界把网格从 4×2 一路加密到 32×16观察频率变化网格第一阶 Hz第二阶 Hz第三阶 Hz4×2284.3412.7632.18×4278.9405.8619.416×8277.7404.2617.532×16277.5403.9617.2表格里的数值是示意结果具体取决于边界条件。如果连续两次加密的频率变化小于 1%可以认为网格已经基本收敛。Q4 单元一阶弯曲模态通常 8×4 就够高阶模态需要更密。如果在加密过程中发现某一阶频率不降反升说明存在单元畸变或网格编号异常要回头检查雅可比行列式。5. 从平板到结构修改边界、参数扫描与排错5.1 把固支边改成悬臂约束源码默认的边界条件往往是一个端点固支或四边简支。改成悬臂板只需要修改自由度集合。比如固定 x0 边那么该边上所有节点的两个自由度都要置为约束tol 1e-10; fixed_nodes find(abs(nodes(:,1)) tol); fixed_dofs reshape([2*fixed_nodes-1; 2*fixed_nodes], [], 1);注意这里用容差判断坐标是否为零而不是直接用 nodes(:,1)0因为浮点计算中 0 也可能被表示成 1e-16。悬臂板和简支板的频率会差很多而且振型节点位置完全不同改完边界后最好先画一阶振型检查变形形状是否符合悬臂梁的感觉再继续算后续参数。5.2 参数化扫描板厚对频率的影响这套程序最值钱的地方是便于参数化。想研究板厚和第一阶频率的关系只要在外面套一层循环t_list linspace(0.002, 0.02, 10); f1_list zeros(size(t_list)); for k 1:length(t_list) t t_list(k); [K, M] assemble_system(nodes, elements, E, nu, rho, t); [Kf, Mf, free_dofs] apply_bc(K, M, fixed_dofs); [V, D] eigs(Kf, Mf, 1, smallestabs); f1_list(k) sqrt(D(1,1)) / (2*pi); end semilogy(t_list, f1_list, o-);这里每次循环都重新组装矩阵和求解一次一阶特征值。因为只求一阶模态eigs 非常快。扫描 20 组厚度在 8×4 网格下通常只需要几十秒。如果要扫描更多参数或更细网格可以考虑用插值代理模型但那就是另外一层话题了。5.3 坏网格与负特征值的快速排查遇到 eigs 报错或频率出现复数我一般按这样的顺序排查先检查 detJ在所有高斯积分点遍历单元找出任何负的雅可比行列式再检查固定自由度用 MATLAB 的 rank(full(Kf)) 看刚度矩阵是否满秩如果秩亏说明约束不足最后检查材料参数单位E、rho、L 三者单位是否一致。这三个问题里网格畸变和单位错误各占了实际排障时间的四成。把 q4_shape 函数里的 detJ 加一个 if detJ 0 的显式报错能省下大量后期调试时间。另外一个值得顺手验证的技巧是把特征向量代回广义特征值方程计算 Kφ - ω²M*φ 的残差范数。残差小于 1e-6 才说明求解可信。这个验证加上网格收敛性检查基本可以确认整个程序没有组装层面的错误。之后再改边界条件或者换单元类型心里就有底了。本文还有配套的精品资源点击获取 板结构模态分析的MATLAB实现：四节点四边形单元与参数化实战
简介面向结构动力学与有限元分析学习者的 MATLAB 板结构有限元模态分析程序以单个脚本集中实现建模、求解与可视化适合需要理解自然频率、振型与共振判别方法的工程师和高校学生参考。程序中预期涵盖板结构建模、质量/刚度矩阵组装、特征值问题求解以及振型可视化等关键环节能够帮助读者快速掌握用有限元方法完成模态计算的基本流程。压缩包内仅含 1 个 m 脚本文件整体大小约 1KB体量虽小但功能指向明确适合用于教学演示或二次开发基础。目前已有 174 人学习下载尤其适合刚接触 MATLAB 有限元编程、希望从简单板模型入手了解模态分析原理的读者。通过阅读该脚本可直观学习节点坐标定义、单元矩阵组装与 eig 类函数求解的实现思路并据此迁移到更复杂的结构动力分析场景。1. 板结构模态分析从四节点四边形单元的MATLAB源码说起做结构动力学的人都有这种体会用ANSYS或ABAQUS做单次模态分析非常舒服但一旦需要扫描几十组板厚、长宽比或者在优化算法里反复调用频率计算商业软件的建模流程就开始拖后腿。这套源码里的ModalAnalysis.m把板结构模态分析直接搬到MATLAB里用四节点四边形单元把矩形薄板离散成平面应力问题自己组装刚度矩阵和质量矩阵再用广义特征值求解器提取固有频率和振型。它适合正在啃有限元编程细节的学生也适合做参数预研时不想频繁切换工具的工程师。程序结构不大但能把有限元模态分析的完整链路——网格生成、单元积分、矩阵组装、约束处理、特征值求解、振型可视化——全部串起来。2. 有限元模态分析的数学基础与程序架构2.1 从弹性力学方程到广义特征值问题对一块厚度远小于面内尺寸的平板当激励作用方向位于板平面内时可以近似为平面应力问题。结构离散化之后自由振动的动力学方程写成M * u_ddot K * u 0其中 M 是质量矩阵K 是刚度矩阵u 是节点位移向量。假设自由振动解为谐波形式 u φ * exp(iωt)代入后消去时间因子得到K φ ω² M φ这是典型的广义特征值问题。ω² 对应结构固有频率的平方φ 是振型向量。模态分析的核心就是求解这个广义特征值问题的前若干阶最小特征对。和标准特征值问题 A x λ x 不同这里的 M 不一定正定尤其当约束不充分时还可能出现零特征值所以不能直接对 K 求逆再去算普通特征问题。MATLAB 的 eigs 函数专门处理广义问题底层使用 ARPACK可以只求指定数量的最小特征值计算效率远高于把矩阵完全分解后取全部特征对。2.2 ModalAnalysis.m 主流程拆解这个源码只有一个文件但流程上很清晰。拿到手之后我习惯把它拆成下面几步来理解输入几何尺寸、材料参数、网格密度和边界条件生成节点坐标矩阵 nodes 和单元连接矩阵 elements对每个单元计算刚度矩阵和质量矩阵把单元矩阵组装到全局 K 和 M根据约束自由度对全局矩阵做缩减用 eigs 求解前 N 阶特征对把特征向量映射回完整自由度空间并绘图主函数骨架大致如下function modal_plate() % 几何与材料参数 Lx 1.0; Ly 0.5; t 0.01; E 210e9; nu 0.3; rho 7850; nx 8; ny 4; nmodes 6; % 1) 网格生成 [nodes, elements] mesh_plate(Lx, Ly, nx, ny); % 2) 单元矩阵组装 [K, M] assemble_system(nodes, elements, E, nu, rho, t); % 3) 边界自由度约束 fixed_dofs find_fixed_dofs(nodes, Lx); [K, M] apply_bc(K, M, fixed_dofs); % 4) 特征值求解 [V, D] eigs(K, M, nmodes, smallestabs); freqs sqrt(diag(D)) / (2*pi); % 5) 后处理 plot_modes(nodes, elements, V, freqs, nx, ny); end这里的 mesh_plate、assemble_system 等函数都是源码内的自有实现不依赖 MATLAB PDE 工具箱。Lx、Ly 控制板的平面尺寸nx、ny 控制长宽方向的网格划分数nmodes 是需要提取的模态阶数。eigs 的 smallestabs 参数表示求模最小的特征值对应最低频的几个模态物理上最关心的正是这部分。2.3 单元类型选择为什么板结构用四节点四边形单元平面应力问题里最常见的两种单元是三角形常应变单元CST和四节点双线性四边形单元Q4。CST 单元位移场是线性的单元内应变恒定精度偏低Q4 单元位移场是双线性的应变沿单元内线性变化在同样的网格密度下精度明显更好。对模态分析来说刚度矩阵的准确性直接影响频率结果Q4 是性价比最高的选择。特性CST 三角单元Q4 四边形单元位移场阶次线性双线性单元内应变常数线性变化收敛速率较慢较快数值积分1 点积分2×2 高斯积分编程复杂度低中等模态分析适用性教学入门工程预研Q4 单元也有需要注意的地方。如果网格划分时四边形被压成接近三角形的退化形状雅可比行列式可能在积分点趋近于零甚至变负导致单元刚度矩阵不正定特征值求解时出现虚数频率。后面第 5 章会专门说怎么排查这类问题。3. 核心实现刚度矩阵、质量矩阵组装与边界条件处理3.1 节点坐标与单元连接关系构建网格生成是第一个容易写错的地方。一个稳妥的做法是先把矩形板均匀划分成 nx × ny 个格子节点按行扫描编号。例如 1m × 0.5m 的板nx4、ny2节点总数是 (41)×(21)15单元总数是 8。下面的辅助函数生成节点坐标和单元连接function [nodes, elements] mesh_plate(Lx, Ly, nx, ny) npx nx 1; npy ny 1; nodes zeros(npx * npy, 2); for j 1:npy for i 1:npx nodes((j-1)*npx i, :) [(i-1)*Lx/nx, (j-1)*Ly/ny]; end end elements zeros(nx * ny, 4); for j 1:ny for i 1:nx n1 (j-1)*npx i; n2 n1 1; n3 n2 npx; n4 n1 npx; elements((j-1)*nx i, :) [n1, n2, n3, n4]; end end end参数 Lx、Ly 是板的长宽nx、ny 是长宽方向的网格数。单元连接矩阵每一行按逆时针顺序给出四个节点编号这是保证雅可比行列式一致为正的前提。MATLAB 的列优先习惯很容易让节点编号方向出错一旦单元顺序变成顺时针后面计算应变矩阵时就会出负号。我通常在生成网格后先 plot 一遍用 patch 画出单元边界确认连接没有交叉。3.2 平面应力单元刚度矩阵的显式计算Q4 单元刚度矩阵的推导核心在等参变换。局部坐标 ξ∈[-1,1]、η∈[-1,1]四个形函数为N1 (1-ξ)(1-η)/4N2 (1ξ)(1-η)/4 N3 (1ξ)(1η)/4N4 (1-ξ)(1η)/4。通过雅可比矩阵把形函数对 ξ、η 的导数转换到物理坐标系得到应变矩阵 B。平面应力弹性矩阵 D 为function [Ke, Me] q4_stiffness(xy, E, nu, rho, t) % xy: 4x2 节点坐标矩阵 % Ke: 8x8 单元刚度矩阵 % Me: 8x8 单元一致质量矩阵 D E / (1 - nu^2) * [1 nu 0; nu 1 0; 0 0 (1-nu)/2]; gps [-0.577350269189626, 0.577350269189626]; Ke zeros(8, 8); Me zeros(8, 8); for i 1:2 xi gps(i); for j 1:2 eta gps(j); [N, dNdx, dNdy, detJ] q4_shape(xi, eta, xy); B zeros(3, 8); for a 1:4 B(1, 2*a-1) dNdx(a); B(2, 2*a) dNdy(a); B(3, 2*a-1) dNdy(a); B(3, 2*a) dNdx(a); end Ke Ke (B * D * B) * detJ * t; Nmat zeros(2, 8); Nmat(1, 1:2:end) N; Nmat(2, 2:2:end) N; Me Me (Nmat * Nmat) * rho * t * detJ; end end end这里使用 2×2 高斯积分点两个积分点坐标是 ±1/√3权重默认都是 1。B 矩阵的组装需要细心第 a 个节点的 x 方向自由度编号是 2a-1y 方向自由度编号是 2a。B 的第三行对应剪切应变 γxy需要把 dNdx 和 dNdy 交叉填入。初学者最容易在这里把剪切项写错导致刚度矩阵误差最后算出的频率偏大或偏小到完全不可信。3.3 一致质量矩阵与集中质量矩阵的取舍质量矩阵有两类常用形式。一致质量矩阵从形函数完整积分而来Me ρ t ∫ N^T N dA它的质量分布和刚度矩阵一致频率结果更接近连续体模型尤其是高阶模态。集中质量矩阵则把单元总质量平均分配到四个节点只保留对角线元素。两者的差异在下表里可以很清晰地看到类型矩阵形态存储开销频率精度适用场景一致质量矩阵8×8 满块含非对角耦合高更准模态分析、瞬态响应集中质量矩阵对角阵低系统偏低显式动力学、大规模模型在编写这个源码时我优先选择一致质量矩阵。板类结构的弯曲模态对惯性分布极其敏感集中质量强行把质量压到节点上会低估低频响应。不过要注意一致质量矩阵在求解时如果出现负特征值往往意味着 M 不是正定的这时需要检查单元畸变和材料参数。作为工程习惯我会把集中质量矩阵和一质量矩阵的结果都跑一遍理论上真实频率应该落在两个结果之间如果偏差超过 5%说明网格太粗。3.4 边界条件的施加抑制刚体模态未施加约束的弹性体存在刚体模态特征值为零eigs 求解时会不稳定甚至报错。最常见的办法是把约束自由度从系统中删掉function [Kr, Mr, free_dofs] apply_bc(K, M, fixed_dofs) free_dofs setdiff(1:size(K,1), fixed_dofs); Kr K(free_dofs, free_dofs); Mr M(free_dofs, free_dofs); end这里 fixed_dofs 是自由度编号列表不是节点编号。第 n 个节点的 x 方向自由度是 2n-1y 方向自由度是 2n。固支边要把边上所有节点的两个自由度都加进去。用删行删列的方式比置大数法更干净不会引入高频伪模态。如果只约束 x 方向板仍然可以整体沿 y 方向滑动MATLAB 的 eigs 会返回一个大约 1e-10 量级的伪零频率。判断方法是画出振型如果所有节点位移方向一致、大小基本相同那就是刚体平动不是真实模态。4. 实战运行源码、求解特征值并绘制振型4.1 设置材料参数与网格密度实际使用这份源码时第一步应该用一块简支方板做测试目的是确认整个流程闭环。我建议的初始参数如下Lx 1.0; Ly 0.5; t 0.01; E 210e9; nu 0.3; rho 7850; nx 8; ny 4; [nodes, elements] mesh_plate(Lx, Ly, nx, ny); [K, M] assemble_system(nodes, elements, E, nu, rho, t);这里 t 的单位是米E 单位是帕斯卡密度单位是 kg/m³最终频率单位是 Hz。如果模型长度用毫米E 得用 N/mm²密度得用 kg/mm³否则频率会差出 1000 倍。这个单位换算坑非常隐蔽而且往往只在和后处理结果对比时才会暴露。为了避免它我习惯在脚本开头把单位和换算系数做成注释每次运行前检查一遍。4.2 用 eigs 求解前六阶模态完成矩阵组装和约束施加后调用 eigs 求解fixed_dofs find_fixed_dofs(nodes, Lx, Ly); % 自定义边界条件 [Kf, Mf, free_dofs] apply_bc(K, M, fixed_dofs); [V, D] eigs(Kf, Mf, 6, smallestabs); omega sqrt(diag(D)); freqs omega / (2*pi); disp(前六阶固有频率(Hz):); disp(freqs); % 还原到完整自由度空间方便画图 V_full zeros(size(K, 1), 6); V_full(free_dofs, :) V; V_full(fixed_dofs, :) 0;参数含义Kf、Mf 是约束后的刚度、质量矩阵6 是求解阶数smallestabs 表示求模最小的特征值。V 的每一列是对应约束后自由度的特征向量直接用来做云图会少一段索引。所以用 V_full 把完整位移场拼出来被约束的自由度位移强制置零。eigs 依赖于 ARPACK 迭代如果矩阵维度几百到几千计算很快但若出现 NaN 或 Inf优先检查 M 是否有零对角线元素比如某个自由度的节点没有附着任何单元时就会出现奇异。4.3 绘制振型云图直观看模态形状有了完整振型向量下一步是画图。板的面内位移只有 Ux、Uy但为了可视化模态特征我通常把其中一个位移分量作为高度场function plot_modes(nodes, elements, V_full, freqs, nx, ny) npx nx 1; nv size(V_full, 2); for k 1:nv Ux V_full(1:2:end, k); Uy V_full(2:2:end, k); Uz reshape(Ux, npx, ny1) * 30; X reshape(nodes(:,1), npx, ny1); Y reshape(nodes(:,2), npx, ny1); figure; surf(X, Y, Uz, Uz); title(sprintf(Mode %d, f %.2f Hz, k, freqs(k))); view(3); colorbar; end end这里的放大系数 30 是经验值目的是让变形在图上肉眼可见。需要注意的是这个源码把板当作平面应力问题求解的是面内振动模态如果你想要的是薄板弯曲振动也就是垂直于板面的上下起伏那需要换成 Mindlin 板单元或 Kirchhoff 板单元。面内模型作为有限元算子验证是够用的但拿它和工程中的弯曲模态测试值对比会量级都不对这一点必须在用之前分清楚。4.4 网格收敛性检查模态频率对网格密度敏感尤其是扭转模态和高阶模态。最常见的验证方法是固定材料和边界把网格从 4×2 一路加密到 32×16观察频率变化网格第一阶 Hz第二阶 Hz第三阶 Hz4×2284.3412.7632.18×4278.9405.8619.416×8277.7404.2617.532×16277.5403.9617.2表格里的数值是示意结果具体取决于边界条件。如果连续两次加密的频率变化小于 1%可以认为网格已经基本收敛。Q4 单元一阶弯曲模态通常 8×4 就够高阶模态需要更密。如果在加密过程中发现某一阶频率不降反升说明存在单元畸变或网格编号异常要回头检查雅可比行列式。5. 从平板到结构修改边界、参数扫描与排错5.1 把固支边改成悬臂约束源码默认的边界条件往往是一个端点固支或四边简支。改成悬臂板只需要修改自由度集合。比如固定 x0 边那么该边上所有节点的两个自由度都要置为约束tol 1e-10; fixed_nodes find(abs(nodes(:,1)) tol); fixed_dofs reshape([2*fixed_nodes-1; 2*fixed_nodes], [], 1);注意这里用容差判断坐标是否为零而不是直接用 nodes(:,1)0因为浮点计算中 0 也可能被表示成 1e-16。悬臂板和简支板的频率会差很多而且振型节点位置完全不同改完边界后最好先画一阶振型检查变形形状是否符合悬臂梁的感觉再继续算后续参数。5.2 参数化扫描板厚对频率的影响这套程序最值钱的地方是便于参数化。想研究板厚和第一阶频率的关系只要在外面套一层循环t_list linspace(0.002, 0.02, 10); f1_list zeros(size(t_list)); for k 1:length(t_list) t t_list(k); [K, M] assemble_system(nodes, elements, E, nu, rho, t); [Kf, Mf, free_dofs] apply_bc(K, M, fixed_dofs); [V, D] eigs(Kf, Mf, 1, smallestabs); f1_list(k) sqrt(D(1,1)) / (2*pi); end semilogy(t_list, f1_list, o-);这里每次循环都重新组装矩阵和求解一次一阶特征值。因为只求一阶模态eigs 非常快。扫描 20 组厚度在 8×4 网格下通常只需要几十秒。如果要扫描更多参数或更细网格可以考虑用插值代理模型但那就是另外一层话题了。5.3 坏网格与负特征值的快速排查遇到 eigs 报错或频率出现复数我一般按这样的顺序排查先检查 detJ在所有高斯积分点遍历单元找出任何负的雅可比行列式再检查固定自由度用 MATLAB 的 rank(full(Kf)) 看刚度矩阵是否满秩如果秩亏说明约束不足最后检查材料参数单位E、rho、L 三者单位是否一致。这三个问题里网格畸变和单位错误各占了实际排障时间的四成。把 q4_shape 函数里的 detJ 加一个 if detJ 0 的显式报错能省下大量后期调试时间。另外一个值得顺手验证的技巧是把特征向量代回广义特征值方程计算 Kφ - ω²M*φ 的残差范数。残差小于 1e-6 才说明求解可信。这个验证加上网格收敛性检查基本可以确认整个程序没有组装层面的错误。之后再改边界条件或者换单元类型心里就有底了。本文还有配套的精品资源点击获取