简介针对板结构动力学特性分析需求这份资源提供基于MATLAB自编有限元程序的完整解决方案面向机械、土木、航空航天等领域的学生、工程师与科研人员用于求解一块板的固有频率和频率响应。通过定义板的几何参数、材料属性、边界条件与加载情况可构建有限元模型并求解动力学方程为评估结构稳定性、避免共振及优化设计提供依据。包内共3个文件包含ban.m主程序与2个txt说明文件分别对应板响应求解和板固有频率求解文件总大小仅3KB整体轻量、结构清晰便于快速理解与二次开发。目前已有180人学习下载适合正在学习有限元编程或需要验证板类结构动态特性的读者。借助该程序使用者可掌握板单元建模、刚度矩阵与质量矩阵组装、特征值求解及频率响应的基本流程同时txt文档中还提供了关键公式和求解步骤说明可帮助减少自行推导与调试时间直接应用于简单板结构的动力学分析或教学演示。1. MATLAB 有限元算板固有频率与频率响应先避开两个坑拿到一块板要算前几阶固有频率和稳态频率响应最省事的办法是在 MATLAB 里做一次有限元特征值分析。但工程的坑往往不在求解器本身而在单元选型、边界约束和阻尼处理。不少人先把板当成梁单元去算得出一组频率后才发现高阶模态完全对不上也有人固有频率算对了频率响应的幅值单位又搞错曲线画出来和实验差好几个数量级。下面把这块板的完整计算链路过一遍控制方程、板单元离散、K/M 矩阵组装、eigs 求广义特征值、模态叠加法算频响函数每一段都会给出可以直接复制的 MATLAB 代码和参数表。2. 板振动控制方程与有限元离散从偏微分方程到特征值问题2.1 一块板与一根梁的本质差异控制方程上的第一处分歧梁的控制方程是四阶常微分方程板则升级成二维双调和方程这是第一个必须转变认知的地方。各向同性薄板在横向分布载荷 p(x,y,t) 下的 Kirchhoff 控制方程为D∇⁴w ρh ∂²w/∂t² p(x,y,t)其中 w 是中面挠度h 为板厚ρ 为材料密度D 是弯曲刚度 D Eh³/(12(1-ν²))。二维 Laplacian 的四次方意味着边界上需要两个条件而不是梁的两个条件简单平移。更关键的是板的多边形边界和曲率耦合使得同一阶模态形状在不同方向上变形特征完全不同这决定了后续有限元每个节点至少要有挠度和两个转角三个自由度。2.2 Kirchhoff 板模型与 Mindlin 板模型怎么选看 h/L选错板理论是算完结果才发现不准的最常见原因。判断标准只有一个厚度与特征尺寸的比值 h/L。以一块边长 1 m、厚 10 mm 的钢板为例h/L 0.01Kirchhoff 薄板理论完全适用但如果换成 20 mm 厚、同样边长的板剪切变形和转动惯量开始贡献误差再套 Kirchhoff 公式基频会偏高几个百分点。板类型h/L 范围适用模型MATLAB 离散时需要注意薄板小于 1/10Kirchhoff 板每节点 w、θx、θy 三个自由度选取经典板单元中厚板1/10 到 1/5Mindlin-Reissner 板需考虑横向剪切应变避免剪切锁死厚板大于 1/5三维实体或高阶壳板假设本身不再成立建议直接做实体单元工程上绝大多数平面板件落在薄板区间本文的代码也以 Kirchhoff 板为准。Mindlin 板虽然可以做但如果网格不够密剪切锁死会让固有频率大幅偏高这是初学 MATLAB 有限元时经常踩的坑。2.3 离散成 K、M 之后的广义特征值问题板的偏微分方程没有解析通解只有少数规则边界有封闭解。四边简支矩形板是经典反例之外的少数好运情况其固有频率为ω_mn π²[(m²/a²) (n²/b²)]·sqrt(D/(ρh))其中 a、b 是矩形边长m、n 是半波数。这个公式有两个用途第一作为后续有限元代码的验证基准第二说明频率不是简单按模态阶数线性排列的板的长宽比会改变模态顺序这也是为什么扫频图的峰值顺序经常出乎意料。有限元离散后振动方程变成 2 阶常微分方程组M·d²u/dt² K·u F(t)令 u U·e^(jωt) 并忽略阻尼就得到广义特征值问题 K·U ω²M·U。这里的 ω 是圆频率固有频率 f ω/(2π)。对于无限自由度连续体离散后 K 和 M 的规模决定你能算多少阶模态也决定高阶频率的精度。下一章的板单元组装目标就是生成这两个矩阵。3. 用 MATLAB 手写矩形板单元组装 K、M 并提取前几阶固有频率3.1 板单元自由度与 ACM 元的形函数写过梁单元的读者对刚度矩阵组装并不陌生板单元只是把自由度从每节点 2 个变成 3 个挠度 w、绕 x 轴转角 θx ∂w/∂y、绕 y 轴转角 θy ∂w/∂x。常用的 4 节点矩形板单元是 ACM 元每个单元 12 个自由度形函数基于不完全 4 次多项式1、ξ、η、ξ²、ξη、η²、ξ³、ξ²η、ξη²、η³、ξ³η、ξη³ACM 元不是严格协调元单元边界的法向转角并不完全连续但它通过 patch test在固有频率计算中收敛性良好是经典板单元中最容易在 MATLAB 里手写的一种。实现时不以显式形函数公式为入口而是建立 12×12 的节点自由度插值矩阵再对它求逆这样代码更短且不容易抄错。3.2 生成网格和单元矩阵组装的 MATLAB 代码以一块 1 m × 1.2 m、厚 10 mm 的钢板为例材料 E 210 GPa、ν 0.3、ρ 7850 kg/m³用 20 × 16 网格做划分。网格生成与单元矩阵组装代码如下% 板尺寸与材料参数 LX 1.0; LY 1.2; h 0.01; E 210e9; nu 0.3; rho 7850; % 网格nx, ny 为单元数 nx 20; ny 16; x linspace(0, LX, nx1); y linspace(0, LY, ny1); [X, Y] meshgrid(x, y); % 节点编号与坐标 nodeId reshape(1:numel(X), size(X)); nodeX X(:); nodeY Y(:); % 单元矩阵每行4个节点逆时针排列 elem zeros(nx*ny, 4); for ix 1:nx for iy 1:ny e (iy-1)*nx ix; n1 nodeId(iy, ix); n2 nodeId(iy, ix1); n3 nodeId(iy1, ix1); n4 nodeId(iy1, ix); elem(e,:) [n1 n2 n3 n4]; end end % 组装全局 K、M ndof length(nodeX) * 3; K sparse(ndof, ndof); M sparse(ndof, ndof); for e 1:size(elem,1) n elem(e,:); a max(nodeX(n)) - min(nodeX(n)); b max(nodeY(n)) - min(nodeY(n)); [Ke, Me] acm_rect_ke(a, b, E, nu, rho, h); % 节点自由度映射每个节点连续3个自由度 dofMap reshape([3*(n-1)1; 3*(n-1)2; 3*(n-1)3], 1, 12); K(dofMap, dofMap) K(dofMap, dofMap) Ke; M(dofMap, dofMap) M(dofMap, dofMap) Me; endmeshgrid 生成的 X、Y 展开后列方向是 y 轴变化这在 nodeId 和 elem 的矩阵填充里要保持同步。单元边长 a、b 直接从节点坐标差计算是为了避免从 IX 和 IY 索引重复换算导致的越界错误。ACM 元的单元矩阵需要单独实现。function [Ke, Me] acm_rect_ke(a, b, E, nu, rho, h) % ACM 矩形板单元每节点自由度 [w; dw/dx; dw/dy] % 单元坐标系原点在矩形中心x 方向 [-a/2, a/2] D0 E*h^3 / (12*(1-nu^2)); C D0 * [1 nu 0; nu 1 0; 0 0 (1-nu)/2]; % 12 个多项式基底 P (xi,eta) [1 xi eta xi^2 xi*eta eta^2 ... xi^3 xi^2*eta xi*eta^2 eta^3 xi^3*eta xi*eta^3]; dPdxi (xi,eta) [0 1 0 2*xi eta 0 3*xi^2 2*xi*eta eta^2 0 3*xi^2*eta eta^3]; dPdeta (xi,eta) [0 0 1 0 xi 2*eta 0 xi^2 2*xi*eta 3*eta^2 xi^3 3*xi*eta^2]; ddPddxi (xi,eta) [0 0 0 2 0 0 6*xi 2*eta 0 0 6*xi*eta 0]; ddPddeta (xi,eta) [0 0 0 0 0 2 0 0 2*xi 6*eta 0 6*xi*eta]; ddPddxidy (xi,eta) [0 0 0 0 1 0 0 2*xi 2*eta 0 3*xi^2 3*eta^2]; % 4 个节点的自然坐标 xiN [-1 1 1 -1]; etaN [-1 -1 1 1]; % 构造自由度插值矩阵 A12x12 A zeros(12,12); row 1; for i 1:4 xi xiN(i); eta etaN(i); A(row, :) P(xi, eta); row row 1; A(row, :) (2/a) * dPdxi(xi, eta); row row 1; A(row, :) (2/b) * dPdeta(xi, eta); row row 1; end invA inv(A); % 2x2 高斯积分 g 1/sqrt(3); gp1 [-g g g -g]; gp2 [-g -g g g]; J a*b / 4; Ke zeros(12,12); Me zeros(12,12); for q 1:4 xi gp1(q); eta gp2(q); Nv P(xi, eta) * invA; % 形函数值 1x12 B [ - (4/a^2) * ddPddxi(xi, eta); - (4/b^2) * ddPddeta(xi, eta); - (8/(a*b)) * ddPddxidy(xi, eta) ] * invA; Ke Ke B * C * B * J; Me Me Nv * Nv * J; end Me rho * h * Me; end代码里 A 矩阵的每一行对应一个自由度的插值条件。前三行是节点 1 的 w、dw/dx、dw/dy 在 12 个基底上的取值后续节点同理。求逆后乘基底向量就得到形函数在积分点的值B 矩阵则由二阶导数基底乘 invA 得到。高斯积分选 2×2 点对于含 3 次和 4 次项的单元刚度积分正好够用换更高阶高斯点不会改变结果只会增加计算量。3.3 求解广义特征值与振型归位组装完成后用 eigs 求广义特征值问题的前若干阶低频模态% 先施加边界条件见第 4 章这里假设 freeDof 已定义 [Vf, Df] eigs(K(freeDof, freeDof), M(freeDof, freeDof), 8, smallestabs); freq sqrt(diag(Df)) / (2*pi); % 特征向量映射回全局自由度 V zeros(ndof, 8); V(freeDof, :) Vf;eigs 求解广义特征值问题时smallestabs 返回模最小的 8 个特征值也就是最低频段。老版本 MATLAB 用 sm两者语义相同。如果模型没有约束或约束不足K 半正定eigs 会提示矩阵奇异这时先检查 fixedDof 的长度。特征向量 V 的每一列是模态振型但符号没有物理意义正负翻转只影响模态形状的显示不影响频率和频响计算。4. 不同边界条件下板的固有频率对比与网格收敛性检查4.1 简支、固支、悬臂板的约束自由度怎么给边界条件在板单元里直接表现为对节点自由度的约束。四边简支约束四个边全部节点的 w 自由度转角 θx、θy 不约束因为简支边允许截面转动但阻止挠度四边固支则把边界节点的 w、θx、θy 全部约束为 0。悬臂板是一条边固支、其余边自由约束方式与固支边相同。% 四边简支约束所有边界节点的 w bottom nodeId(1,:); top nodeId(end,:); left nodeId(:,1); right nodeId(:,end); edgeNode unique([bottom, top, left, right]); % 每个节点的第 1 个自由度是 w fixedDof 3*(edgeNode-1) 1; % 四边固支约束边界节点全部 3 个自由度 fixedDof []; for i 1:length(edgeNode) fixedDof [fixedDof, 3*(edgeNode(i)-1) (1:3)]; end freeDof setdiff(1:ndof, unique(fixedDof));简支边只约束 w 是工程上最常用的处理方式。严格从板理论讲简支边还需要弯矩为零的条件但 ACM 元中弯矩为零是自然边界条件有限元解会自动满足因此不需要额外施加。固支边则必须约束转动自由度漏掉 θy 会让算出来的频率偏低。4.2 四边简支板的收敛性网格从 4×4 加密到 24×24用第 2 章的解析解 49.2 Hz 做基准1 m × 1 m、厚 10 mm 钢板可以观察到板单元随网格加密的收敛趋势。网格规模自由度总数基频计算值与解析解误差量级说明4×4455% 到 10%网格太稀高阶模态误差更大8×82431% 到 3%前 3 阶可用高阶仍偏大16×168670.2% 到 1%工程计算够用24×2418750.1% 左右适合作为收敛基准ACM 元按 h² 收敛网格每加密一倍误差大致缩小到原来的四分之一。所以从 8×8 加密到 16×16 时误差会明显跳变这是判断计算结果是否可信的一个快速手段。多算几阶模态时前 5 阶比基频收敛慢如果只需要前 3 阶16×16 网格完全够用如果要取 10 阶以上建议加密到 32×32。4.3 边界条件对固有频率的倍数关系同一块板的边界条件不同固有频率能差出近一倍。四边固支板基频约是四边简支板的 1.8 到 2.0 倍悬臂板则比四边简支低得多因为自由边降低了整体刚度。实际计算时可以先跑一次 8×8 网格快速估一遍频率范围再加密网格做正式计算这样比直接上 32×32 省很多时间也能避免边界约束写错导致结果偏差过大而事后返工。5. 板的频率响应模态叠加法计算幅频特性5.1 频率响应函数公式和质量归一化稳态频率响应通常用模态叠加法计算。设激励频率为 ω激励自由度编号为 p响应自由度编号为 q则位移频响函数为H_pq(ω) Σᵣ φ_pᵣ·φ_qᵣ / (ω_r² − ω² 2jζ_rω_rω)其中 φ_r 必须是质量归一化振型满足 φ_rᵀ M φ_r 1。eigs 返回的特征向量虽然关于 M 正交但幅值未归一化直接用会导致频响幅值整体错误。% V 是全局特征向量矩阵M 是全局质量矩阵 nmode size(V,2); Vn zeros(size(V)); for r 1:nmode mr V(:,r) * M * V(:,r); Vn(:,r) V(:,r) / sqrt(mr); end5.2 用 MATLAB 计算点激励下的幅频曲线% 激励点作用在节点 p响应取节点 q p 100; q 205; % 按实际网格节点编号修改 zeta 0.01 * ones(nmode,1); % 模态阻尼比 fscan linspace(20, 120, 2001); % 频率扫描范围单位 Hz omega 2*pi * fscan; Hresp zeros(size(omega)); for k 1:length(omega) w omega(k); s 0; for r 1:nmode w_r 2*pi * freq(r); s s Vn(p,r)*Vn(q,r) / (w_r^2 - w^2 2i*zeta(r)*w_r*w); end Hresp(k) s; end figure; semilogy(fscan, abs(Hresp)); xlabel(频率 (Hz)); ylabel(位移频响幅值 (m/N));频率扫描的间隔要足够密半功率带宽通常只有固有频率的百分之几2001 个点在 20 到 120 Hz 内大约 0.05 Hz 分辨率足够分辨峰值。如果曲线毛刺多首先加密 fscan不要急着改阻尼。频响峰值对应的频率应该与第 4 章算出的固有频率一致这是对全部分析链路最直接的验证。5.3 用半功率带宽从频率响应中验证阻尼比响应峰值下降 3 dB 处的两个频率之差除以中心频率就是阻尼比的近似值ζ (ω₂ − ω₁) / (2ω_r)如果频响曲线在固有频率附近算出一个 ζ 0.015而输入用的是 ζ 0.01多半是扫描频率分辨率不够导致峰值偏低、带宽偏大。把频率点数提高一倍再验证一次这个技巧能同时检验频率步长和模态叠加的截断阶数是否足够。之后把激励从点力换成都布压力或者基础激励公式里只改激励自由度的分布向量模态叠加框架不用动。本文还有配套的精品资源点击获取 MATLAB有限元分析板结构固有频率与频响曲线实战
简介针对板结构动力学特性分析需求这份资源提供基于MATLAB自编有限元程序的完整解决方案面向机械、土木、航空航天等领域的学生、工程师与科研人员用于求解一块板的固有频率和频率响应。通过定义板的几何参数、材料属性、边界条件与加载情况可构建有限元模型并求解动力学方程为评估结构稳定性、避免共振及优化设计提供依据。包内共3个文件包含ban.m主程序与2个txt说明文件分别对应板响应求解和板固有频率求解文件总大小仅3KB整体轻量、结构清晰便于快速理解与二次开发。目前已有180人学习下载适合正在学习有限元编程或需要验证板类结构动态特性的读者。借助该程序使用者可掌握板单元建模、刚度矩阵与质量矩阵组装、特征值求解及频率响应的基本流程同时txt文档中还提供了关键公式和求解步骤说明可帮助减少自行推导与调试时间直接应用于简单板结构的动力学分析或教学演示。1. MATLAB 有限元算板固有频率与频率响应先避开两个坑拿到一块板要算前几阶固有频率和稳态频率响应最省事的办法是在 MATLAB 里做一次有限元特征值分析。但工程的坑往往不在求解器本身而在单元选型、边界约束和阻尼处理。不少人先把板当成梁单元去算得出一组频率后才发现高阶模态完全对不上也有人固有频率算对了频率响应的幅值单位又搞错曲线画出来和实验差好几个数量级。下面把这块板的完整计算链路过一遍控制方程、板单元离散、K/M 矩阵组装、eigs 求广义特征值、模态叠加法算频响函数每一段都会给出可以直接复制的 MATLAB 代码和参数表。2. 板振动控制方程与有限元离散从偏微分方程到特征值问题2.1 一块板与一根梁的本质差异控制方程上的第一处分歧梁的控制方程是四阶常微分方程板则升级成二维双调和方程这是第一个必须转变认知的地方。各向同性薄板在横向分布载荷 p(x,y,t) 下的 Kirchhoff 控制方程为D∇⁴w ρh ∂²w/∂t² p(x,y,t)其中 w 是中面挠度h 为板厚ρ 为材料密度D 是弯曲刚度 D Eh³/(12(1-ν²))。二维 Laplacian 的四次方意味着边界上需要两个条件而不是梁的两个条件简单平移。更关键的是板的多边形边界和曲率耦合使得同一阶模态形状在不同方向上变形特征完全不同这决定了后续有限元每个节点至少要有挠度和两个转角三个自由度。2.2 Kirchhoff 板模型与 Mindlin 板模型怎么选看 h/L选错板理论是算完结果才发现不准的最常见原因。判断标准只有一个厚度与特征尺寸的比值 h/L。以一块边长 1 m、厚 10 mm 的钢板为例h/L 0.01Kirchhoff 薄板理论完全适用但如果换成 20 mm 厚、同样边长的板剪切变形和转动惯量开始贡献误差再套 Kirchhoff 公式基频会偏高几个百分点。板类型h/L 范围适用模型MATLAB 离散时需要注意薄板小于 1/10Kirchhoff 板每节点 w、θx、θy 三个自由度选取经典板单元中厚板1/10 到 1/5Mindlin-Reissner 板需考虑横向剪切应变避免剪切锁死厚板大于 1/5三维实体或高阶壳板假设本身不再成立建议直接做实体单元工程上绝大多数平面板件落在薄板区间本文的代码也以 Kirchhoff 板为准。Mindlin 板虽然可以做但如果网格不够密剪切锁死会让固有频率大幅偏高这是初学 MATLAB 有限元时经常踩的坑。2.3 离散成 K、M 之后的广义特征值问题板的偏微分方程没有解析通解只有少数规则边界有封闭解。四边简支矩形板是经典反例之外的少数好运情况其固有频率为ω_mn π²[(m²/a²) (n²/b²)]·sqrt(D/(ρh))其中 a、b 是矩形边长m、n 是半波数。这个公式有两个用途第一作为后续有限元代码的验证基准第二说明频率不是简单按模态阶数线性排列的板的长宽比会改变模态顺序这也是为什么扫频图的峰值顺序经常出乎意料。有限元离散后振动方程变成 2 阶常微分方程组M·d²u/dt² K·u F(t)令 u U·e^(jωt) 并忽略阻尼就得到广义特征值问题 K·U ω²M·U。这里的 ω 是圆频率固有频率 f ω/(2π)。对于无限自由度连续体离散后 K 和 M 的规模决定你能算多少阶模态也决定高阶频率的精度。下一章的板单元组装目标就是生成这两个矩阵。3. 用 MATLAB 手写矩形板单元组装 K、M 并提取前几阶固有频率3.1 板单元自由度与 ACM 元的形函数写过梁单元的读者对刚度矩阵组装并不陌生板单元只是把自由度从每节点 2 个变成 3 个挠度 w、绕 x 轴转角 θx ∂w/∂y、绕 y 轴转角 θy ∂w/∂x。常用的 4 节点矩形板单元是 ACM 元每个单元 12 个自由度形函数基于不完全 4 次多项式1、ξ、η、ξ²、ξη、η²、ξ³、ξ²η、ξη²、η³、ξ³η、ξη³ACM 元不是严格协调元单元边界的法向转角并不完全连续但它通过 patch test在固有频率计算中收敛性良好是经典板单元中最容易在 MATLAB 里手写的一种。实现时不以显式形函数公式为入口而是建立 12×12 的节点自由度插值矩阵再对它求逆这样代码更短且不容易抄错。3.2 生成网格和单元矩阵组装的 MATLAB 代码以一块 1 m × 1.2 m、厚 10 mm 的钢板为例材料 E 210 GPa、ν 0.3、ρ 7850 kg/m³用 20 × 16 网格做划分。网格生成与单元矩阵组装代码如下% 板尺寸与材料参数 LX 1.0; LY 1.2; h 0.01; E 210e9; nu 0.3; rho 7850; % 网格nx, ny 为单元数 nx 20; ny 16; x linspace(0, LX, nx1); y linspace(0, LY, ny1); [X, Y] meshgrid(x, y); % 节点编号与坐标 nodeId reshape(1:numel(X), size(X)); nodeX X(:); nodeY Y(:); % 单元矩阵每行4个节点逆时针排列 elem zeros(nx*ny, 4); for ix 1:nx for iy 1:ny e (iy-1)*nx ix; n1 nodeId(iy, ix); n2 nodeId(iy, ix1); n3 nodeId(iy1, ix1); n4 nodeId(iy1, ix); elem(e,:) [n1 n2 n3 n4]; end end % 组装全局 K、M ndof length(nodeX) * 3; K sparse(ndof, ndof); M sparse(ndof, ndof); for e 1:size(elem,1) n elem(e,:); a max(nodeX(n)) - min(nodeX(n)); b max(nodeY(n)) - min(nodeY(n)); [Ke, Me] acm_rect_ke(a, b, E, nu, rho, h); % 节点自由度映射每个节点连续3个自由度 dofMap reshape([3*(n-1)1; 3*(n-1)2; 3*(n-1)3], 1, 12); K(dofMap, dofMap) K(dofMap, dofMap) Ke; M(dofMap, dofMap) M(dofMap, dofMap) Me; endmeshgrid 生成的 X、Y 展开后列方向是 y 轴变化这在 nodeId 和 elem 的矩阵填充里要保持同步。单元边长 a、b 直接从节点坐标差计算是为了避免从 IX 和 IY 索引重复换算导致的越界错误。ACM 元的单元矩阵需要单独实现。function [Ke, Me] acm_rect_ke(a, b, E, nu, rho, h) % ACM 矩形板单元每节点自由度 [w; dw/dx; dw/dy] % 单元坐标系原点在矩形中心x 方向 [-a/2, a/2] D0 E*h^3 / (12*(1-nu^2)); C D0 * [1 nu 0; nu 1 0; 0 0 (1-nu)/2]; % 12 个多项式基底 P (xi,eta) [1 xi eta xi^2 xi*eta eta^2 ... xi^3 xi^2*eta xi*eta^2 eta^3 xi^3*eta xi*eta^3]; dPdxi (xi,eta) [0 1 0 2*xi eta 0 3*xi^2 2*xi*eta eta^2 0 3*xi^2*eta eta^3]; dPdeta (xi,eta) [0 0 1 0 xi 2*eta 0 xi^2 2*xi*eta 3*eta^2 xi^3 3*xi*eta^2]; ddPddxi (xi,eta) [0 0 0 2 0 0 6*xi 2*eta 0 0 6*xi*eta 0]; ddPddeta (xi,eta) [0 0 0 0 0 2 0 0 2*xi 6*eta 0 6*xi*eta]; ddPddxidy (xi,eta) [0 0 0 0 1 0 0 2*xi 2*eta 0 3*xi^2 3*eta^2]; % 4 个节点的自然坐标 xiN [-1 1 1 -1]; etaN [-1 -1 1 1]; % 构造自由度插值矩阵 A12x12 A zeros(12,12); row 1; for i 1:4 xi xiN(i); eta etaN(i); A(row, :) P(xi, eta); row row 1; A(row, :) (2/a) * dPdxi(xi, eta); row row 1; A(row, :) (2/b) * dPdeta(xi, eta); row row 1; end invA inv(A); % 2x2 高斯积分 g 1/sqrt(3); gp1 [-g g g -g]; gp2 [-g -g g g]; J a*b / 4; Ke zeros(12,12); Me zeros(12,12); for q 1:4 xi gp1(q); eta gp2(q); Nv P(xi, eta) * invA; % 形函数值 1x12 B [ - (4/a^2) * ddPddxi(xi, eta); - (4/b^2) * ddPddeta(xi, eta); - (8/(a*b)) * ddPddxidy(xi, eta) ] * invA; Ke Ke B * C * B * J; Me Me Nv * Nv * J; end Me rho * h * Me; end代码里 A 矩阵的每一行对应一个自由度的插值条件。前三行是节点 1 的 w、dw/dx、dw/dy 在 12 个基底上的取值后续节点同理。求逆后乘基底向量就得到形函数在积分点的值B 矩阵则由二阶导数基底乘 invA 得到。高斯积分选 2×2 点对于含 3 次和 4 次项的单元刚度积分正好够用换更高阶高斯点不会改变结果只会增加计算量。3.3 求解广义特征值与振型归位组装完成后用 eigs 求广义特征值问题的前若干阶低频模态% 先施加边界条件见第 4 章这里假设 freeDof 已定义 [Vf, Df] eigs(K(freeDof, freeDof), M(freeDof, freeDof), 8, smallestabs); freq sqrt(diag(Df)) / (2*pi); % 特征向量映射回全局自由度 V zeros(ndof, 8); V(freeDof, :) Vf;eigs 求解广义特征值问题时smallestabs 返回模最小的 8 个特征值也就是最低频段。老版本 MATLAB 用 sm两者语义相同。如果模型没有约束或约束不足K 半正定eigs 会提示矩阵奇异这时先检查 fixedDof 的长度。特征向量 V 的每一列是模态振型但符号没有物理意义正负翻转只影响模态形状的显示不影响频率和频响计算。4. 不同边界条件下板的固有频率对比与网格收敛性检查4.1 简支、固支、悬臂板的约束自由度怎么给边界条件在板单元里直接表现为对节点自由度的约束。四边简支约束四个边全部节点的 w 自由度转角 θx、θy 不约束因为简支边允许截面转动但阻止挠度四边固支则把边界节点的 w、θx、θy 全部约束为 0。悬臂板是一条边固支、其余边自由约束方式与固支边相同。% 四边简支约束所有边界节点的 w bottom nodeId(1,:); top nodeId(end,:); left nodeId(:,1); right nodeId(:,end); edgeNode unique([bottom, top, left, right]); % 每个节点的第 1 个自由度是 w fixedDof 3*(edgeNode-1) 1; % 四边固支约束边界节点全部 3 个自由度 fixedDof []; for i 1:length(edgeNode) fixedDof [fixedDof, 3*(edgeNode(i)-1) (1:3)]; end freeDof setdiff(1:ndof, unique(fixedDof));简支边只约束 w 是工程上最常用的处理方式。严格从板理论讲简支边还需要弯矩为零的条件但 ACM 元中弯矩为零是自然边界条件有限元解会自动满足因此不需要额外施加。固支边则必须约束转动自由度漏掉 θy 会让算出来的频率偏低。4.2 四边简支板的收敛性网格从 4×4 加密到 24×24用第 2 章的解析解 49.2 Hz 做基准1 m × 1 m、厚 10 mm 钢板可以观察到板单元随网格加密的收敛趋势。网格规模自由度总数基频计算值与解析解误差量级说明4×4455% 到 10%网格太稀高阶模态误差更大8×82431% 到 3%前 3 阶可用高阶仍偏大16×168670.2% 到 1%工程计算够用24×2418750.1% 左右适合作为收敛基准ACM 元按 h² 收敛网格每加密一倍误差大致缩小到原来的四分之一。所以从 8×8 加密到 16×16 时误差会明显跳变这是判断计算结果是否可信的一个快速手段。多算几阶模态时前 5 阶比基频收敛慢如果只需要前 3 阶16×16 网格完全够用如果要取 10 阶以上建议加密到 32×32。4.3 边界条件对固有频率的倍数关系同一块板的边界条件不同固有频率能差出近一倍。四边固支板基频约是四边简支板的 1.8 到 2.0 倍悬臂板则比四边简支低得多因为自由边降低了整体刚度。实际计算时可以先跑一次 8×8 网格快速估一遍频率范围再加密网格做正式计算这样比直接上 32×32 省很多时间也能避免边界约束写错导致结果偏差过大而事后返工。5. 板的频率响应模态叠加法计算幅频特性5.1 频率响应函数公式和质量归一化稳态频率响应通常用模态叠加法计算。设激励频率为 ω激励自由度编号为 p响应自由度编号为 q则位移频响函数为H_pq(ω) Σᵣ φ_pᵣ·φ_qᵣ / (ω_r² − ω² 2jζ_rω_rω)其中 φ_r 必须是质量归一化振型满足 φ_rᵀ M φ_r 1。eigs 返回的特征向量虽然关于 M 正交但幅值未归一化直接用会导致频响幅值整体错误。% V 是全局特征向量矩阵M 是全局质量矩阵 nmode size(V,2); Vn zeros(size(V)); for r 1:nmode mr V(:,r) * M * V(:,r); Vn(:,r) V(:,r) / sqrt(mr); end5.2 用 MATLAB 计算点激励下的幅频曲线% 激励点作用在节点 p响应取节点 q p 100; q 205; % 按实际网格节点编号修改 zeta 0.01 * ones(nmode,1); % 模态阻尼比 fscan linspace(20, 120, 2001); % 频率扫描范围单位 Hz omega 2*pi * fscan; Hresp zeros(size(omega)); for k 1:length(omega) w omega(k); s 0; for r 1:nmode w_r 2*pi * freq(r); s s Vn(p,r)*Vn(q,r) / (w_r^2 - w^2 2i*zeta(r)*w_r*w); end Hresp(k) s; end figure; semilogy(fscan, abs(Hresp)); xlabel(频率 (Hz)); ylabel(位移频响幅值 (m/N));频率扫描的间隔要足够密半功率带宽通常只有固有频率的百分之几2001 个点在 20 到 120 Hz 内大约 0.05 Hz 分辨率足够分辨峰值。如果曲线毛刺多首先加密 fscan不要急着改阻尼。频响峰值对应的频率应该与第 4 章算出的固有频率一致这是对全部分析链路最直接的验证。5.3 用半功率带宽从频率响应中验证阻尼比响应峰值下降 3 dB 处的两个频率之差除以中心频率就是阻尼比的近似值ζ (ω₂ − ω₁) / (2ω_r)如果频响曲线在固有频率附近算出一个 ζ 0.015而输入用的是 ζ 0.01多半是扫描频率分辨率不够导致峰值偏低、带宽偏大。把频率点数提高一倍再验证一次这个技巧能同时检验频率步长和模态叠加的截断阶数是否足够。之后把激励从点力换成都布压力或者基础激励公式里只改激励自由度的分布向量模态叠加框架不用动。本文还有配套的精品资源点击获取