简介一份基于Timoshenko梁理论的声子晶体梁频散特性计算MATLAB程序面向结构动力学、声学超材料及周期性结构方向的研究者和高年级本科生。代码采用12×12矩阵形式的传递矩阵法将梁离散为若干子段逐段构造考虑剪切效应与转动刚度的局部传递矩阵再组装为全局矩阵进而求解声波在周期梁中的传播特性与带隙范围。资源包仅含1个m文件大小约2KB为独立脚本无附加数据文档使用前需具备MATLAB基础及传递矩阵基本概念。目前已有366人学习下载可作为声子晶体梁建模与传递矩阵编程的入门参考。程序内通常涵盖参数定义、单元矩阵计算、频率扫描及可视化模块用户只需调整几何尺寸、材料属性和单元数量即可复现频响曲线与带隙结构为声隔离、声过滤设计提供理论验证手段。1. 为什么声子晶体梁的带隙计算绕不开 12×12 传递矩阵声子晶体梁的带隙计算第一反应多是 COMSOL 扫频或 ANSYS 的 Floquet 周期边界但参数扫描动辄上千次、每次重画网格的场景下解析/半解析的传递矩阵法才是能反复改参数的工具。这套 12×12 的 Timoshenko 梁传递矩阵把周期单元的 3 个平动、3 个转动自由度和对应 6 个力/力矩合成状态向量单元两端状态用矩阵指数直接映射。对声学超材料、周期隔振、宽梁/厚梁振动分析的从业者它比 4×4 的 Euler 梁模型多考虑了剪切变形与转动惯量高频段和粗短单元上带隙边界偏差可达 10% 以上。2. Timoshenko 梁状态向量与 12×12 系统矩阵的构造2.1 为什么要 12 个状态量而不是 4 个声子晶体梁的每一段都是三维结构轴向纵波、扭转波、两个平面内的弯曲波同时存在。如果只取某个平面内的 w、φ、M、V 四个量虽然能算单一弯曲带隙但遇到截面非对称、需要同时考察纵波与弯曲波耦合或者要核对实验里多个模态叠加的透射峰时4×4 就明显不够用。12 个状态量按「位移组在前、力组在后」排列序号状态量含义序号状态量含义1u轴向位移7N轴力2vy 向横向位移8Vyy 向剪力3wz 向横向位移9Vzz 向剪力4θx扭转角10Mx扭矩5θy绕 y 轴转角11My绕 y 轴弯矩6θz绕 z 轴转角12Mz绕 z 轴弯矩这个排序不是唯一的但一旦定下来A 矩阵每行每列都要严格对应。位移组与力组的对偶排列会让 12×12 传递矩阵具有辛对称性特征值两两互为倒数这是判断程序是否写错的第一道检查点。需要说明的是12×12 里的弯曲、轴向、扭转三个子块在均匀截面下互不耦合带隙结构可以看成三个波型独立叠加。真正让 12×12 比多个 4×4 更有价值的地方在于它保留了完整的状态向量结构一旦单元内出现几何偏心、斜支承或者材料各向异性各子块之间会产生耦合项这时 4×4 必须重写而 12×12 只需要在 A 矩阵对应位置补上非零元素。这也是我建议直接维护 12×12 而不是拆成三个独立求解器的原因。2.2 各物理场的微分方程与矩阵填充对简谐振动时间因子 e^{iωt}Timoshenko 梁在一个平面内的状态方程是四个一阶常微分方程dV/dx -ρAω²vdM/dx V - ρIω²φdφ/dx M/(EI)dv/dx -φ V/(κGA)。第一项是横向力平衡第二项是力矩平衡第三项是弯矩-曲率关系第四项含剪切柔度 1/(κGA)是与 Euler 梁的本质区别。写成 4×4 系数矩阵[0, -1, 0, 1/(κGA)] [0, 0, 1/(EI), 0] [0, -ρIω², 0, 1] [-ρAω², 0, 0, 0]轴向波是 2×2 块 du/dx N/(EA)、dN/dx -ρAω²u扭转波是 dθx/dx Mx/(GJ)、dMx/dx -ρIpω²θx。把两个 2×2 和两个 4×4 按状态向量顺序拼进 12×12 系统矩阵 A单元长度为 L 时传递矩阵即 T exp(A·L)。用矩阵指数而不是手解特征方程拼三角函数是因为 expm 对含混合波场、非对称系数的矩阵能一次性处理。下面这段代码对应 12×12 的 timoshenko 梁.m 里最核心的单段传递矩阵函数function T tm_transfer(L, mat, omega) % L: 单元长度(m); mat: 材料与截面参数; omega: 角频率(rad/s) % 返回 12x12 传递矩阵, state_right T * state_left A zeros(12, 12); % 轴向: u - N A(1, 7) 1 / (mat.E * mat.A); A(7, 1) -mat.rho * mat.A * omega^2; % 扭转: theta_x - Mx A(4, 10) 1 / (mat.G * mat.J); A(10, 4) -mat.rho * mat.Ip * omega^2; % x-y 平面弯曲: v - theta_z - Vy - Mz A(2, 6) -1; A(2, 8) 1 / (mat.kappa * mat.G * mat.A); A(6, 12) 1 / (mat.E * mat.Iz); A(12, 6) -mat.rho * mat.Iz * omega^2; A(12, 8) 1; A(8, 2) -mat.rho * mat.A * omega^2; % x-z 平面弯曲: w - theta_y - Vz - My, 块结构与 xy 面相同 A(3, 5) -1; A(3, 9) 1 / (mat.kappa * mat.G * mat.A); A(5, 11) 1 / (mat.E * mat.Iy); A(11, 5) -mat.rho * mat.Iy * omega^2; A(11, 9) 1; A(9, 3) -mat.rho * mat.A * omega^2; T expm(A * L); end逻辑说明A 的每一行对应一个状态量的空间导数比如第 2 行是 dv/dx -θz Vy/(κGA)所以 A(2,6) -1、A(2,8) 1/(κGA)。这个 4×4 块的特征方程展开后是 k⁴ - ρω²(1/E 1/(κG))k² ρ²ω⁴/(EκG) - ρAω²/(EI) 0正是 Timoshenko 色散关系的标准形式低频极限下回到 Euler 梁的 k ≈ (ρAω²/EI)^{1/4}。写代码时我习惯先用符号计算验证特征多项式再进频率扫描避免某个 4×4 块符号写反而带隙整体漂移。参数方面kappa 是剪切修正系数矩形截面取 5/6、圆截面取 9/10G 是剪切模量Ip 是极惯性矩J 是扭转常数实心圆截面两者相等薄壁开口截面必须区分开。2.3 矩阵指数的数值边界expm 不是万能的。当 ω 很高或 L 很长时A·L 的范数变大矩阵指数的辛对称性会被舍入误差破坏典型表现是特征值不再严格两两成对。我一般用单元长细比做判据单元长度小于目标频段最短波长的 1/10 时直接 expm超过就把单元等分成 k 段每段算完再连乘。声子晶体梁的高阶带隙恰好落在高频区这个分段处理要提前写好。注意expm 的输入 A*L 必须无量纲。omega 用 rad/s、L 用来、E/G 用 Pa 时乘积自然无量纲如果把频率写成 Hz 直接代入带隙会全部错位。3. 周期单胞组装与 Bloch 能带扫描的 MATLAB 实现3.1 声子晶体单胞的传递矩阵组装声子晶体梁的周期单胞通常是两种材料交替比如钢-环氧、铝-橡胶这类大阻抗比配对。设材料 a 长度 La、材料 b 长度 Lb波先经过 a 段再经过 b 段单胞传递矩阵就是 T_cell T_b · T_a乘法顺序不能换。下面给组装函数function Tcell tm_unit_cell(La, Lb, matA, matB, omega) Ta tm_transfer(La, matA, omega); Tb tm_transfer(Lb, matB, omega); Tcell Tb * Ta; % 波先过 a 再过 b, 所以 b 左乘 end这段逻辑可以推广到多段单元也能处理集中质量、弹性支承等不连续点每个不连续点对应一个点传递矩阵按位置连乘进总矩阵即可。点矩阵通常是稀疏的只在约束自由度对应的行非零能复用同一个组装框架。3.2 Bloch 条件与传播常数提取无限周期结构里Bloch 定理要求相邻单胞对应位置的状态向量只差相位因子 e^{ikΛ}Λ La Lb 是单胞长度k 是复波数。代入传递关系 T_cell·Ψ_left e^{ikΛ}·Ψ_lefte^{ikΛ} 正是 T_cell 的特征值。所以求能带结构不需要拼全局矩阵只需对每个频率求 T_cell 的 12 个特征值 λi再取传播常数 μi ln(λi)。μ 为纯虚数对应无衰减传播的 Bloch 波落在通带μ 实部非零对应指数衰减或增长落在禁带。禁带成立的充要条件是该频段内没有任何一个 |λ| 在单位圆上。波数 k 的实部决定相位积累虚部决定幅值衰减。频率从零上升时各分支沿色散曲线爬行一旦进入禁带k 的虚部从零跳变为非零对应 Bloch 波在单胞界面反复反射干涉、无法远距传播。传递矩阵法的优势是单次特征值分解覆盖所有波型不需要像有限元那样在每个 k 点单独求解广义特征值问题。Λ 同时出现在 Bloch 相位因子、态密度估计和隔振设计频段三个地方改变 Λ 相当于整体平移带隙配合 4.3 节的标度律可以快速估算目标频段对应的单胞尺寸。3.3 能带扫描代码与参数说明om 2 * pi * linspace(10, 8000, 800); % 10 Hz ~ 8 kHz, 共 800 个频点 mu_real zeros(12, numel(om)); for i 1:numel(om) Tcell tm_unit_cell(La, Lb, matA, matB, om(i)); lam eig(Tcell); mu_real(:, i) sort(real(log(lam))); % 禁带判定只看传播常数实部 end逻辑说明这里不求虚部因为带隙判据只关心传播常数实部是否为 0。sort 按实部升序排列方便后续 imagesc 画色温图时通带与禁带分明。800 个频点对工程判断足够定位带隙边缘时再对边缘频段做二分细化。eig 默认按模排序传递矩阵特征值两两互为倒数排序后衰减曲线成对对称正好框出禁带上下边界。matA、matB 是包含 E、G、rho、A、Iy、Iz、Ip、J、kappa 的结构体常用值如下参数钢环氧说明E (GPa)2103.5弹性模量G (GPa)801.3剪切模量rho (kg/m³)78501200密度kappa5/65/6矩形截面剪切修正系数提示如果扫描结果里完全没有禁带优先检查两段材料的阻抗比是否足够大阻抗比小于 3 时第一带隙会窄到难以辨认。4. 带隙图谱判读与剪切修正、填充率的敏感性4.1 从特征值分布理解通带和禁带T_cell 是辛矩阵12 个特征值两两互为倒数。无阻尼时要么成对在单位圆上要么成对在单位圆内外对称位置。判据「|λ| 偏离 1 即禁带」可简化成设定模最大偏差阈值。阈值要考虑 expm 的数值误差我通常取 1e-6低于这个阈值的微小波动是舍入噪声不应当成带隙边缘。画能带图时横轴是频率纵轴是归一化波数 Re(k)·Λ/π通带分支是连续斜线禁带处没有任何分支穿过衰减常数单独画一张色温图亮带位置与能带图的空隙完全对应两张图交叉检查是排查程序错误最快的路径。另外建议把 mu_real 的第 1 行和第 12 行单独画线这两条曲线分别是衰减最小的两支代表波衰减最慢的通道实验里透射峰往往出现在它们所在频段判断带隙质量时比看平均衰减更可靠。4.2 剪切修正系数与转动惯量的作用Timoshenko 和 Euler 模型的差异集中在 κ 与 ρI 两项。把 κ 从 5/6 改成 1高频带隙上边缘可能移动 3% 到 8%5000 Hz 以上的误差会超过一个频点间距。反过来要复现实验测得的带隙位置κ 常需按实际截面的翘曲情况微调这是 Timoshenko 模型比 Euler 模型更贴近实验的原因。ρIω² 项对低阶弯曲带隙影响小对高阶光学支和纵弯耦合模态不可忽略。另一个工程经验是梁截面高度与单胞长度比值超过 1/5 时剪切变形对一阶弯曲带隙下边缘的影响就会超过 5%这时候不要试图用 Euler 梁替代。4.3 材料配对与填充率的调参带隙中心频率主要由单胞总长与等效波速决定存在标准标度律单胞总长减半带隙中心频率近似翻倍这是快速正确性检查的利器。带隙宽度主要由阻抗比和填充率决定。钢与橡胶配对时阻抗比接近 50第一带隙可以做得很宽但橡胶的阻尼会引入复波数处理时要保留损耗因子把 E、G 换成复模量。填充率 η La/Λ 的扫描经验弯曲主导的带隙在 η≈0.5 附近最宽轴向带隙对 η 不敏感。扫描 η 只需在循环里改 La不需要改矩阵填充逻辑这也是传递矩阵法相对有限元最大的工程优势。4.4 与有限元或实验的交叉验证拿到代码先别信带隙结果做两步验证。第一把 ω 压到接近 0T_cell 特征值的虚部应接近各波型的静态波数低频弯曲波数与 κGA、EI 的关系一致。第二在 COMSOL 里建 12×1 的单胞加 Floquet 周期边界对比第一带隙上下边缘两个方法偏差在 2% 以内说明矩阵符号和材料参数都正确。偏差大时按 2.2 节的状态向量顺序逐行核对重点看力组与位移组的排列是否错位以及 Iy、Iz 是否对调。5. 有限周期梁透射率计算与 expm 数值稳定性处理5.1 级联传递矩阵与矩阵幂的溢出实际隔振梁是有限周期结构N 个单胞的传递矩阵是 T_cell 的 N 次幂。直接 mpower(Tcell, N) 在 N 超过 50 时会因大特征值指数放大而溢出。稳定做法是先特征分解 T_cell V·Λ·V⁻¹对对角特征值逐项求 Λ^N 再重组既得到 T_total 又拿到传播常数。5.2 高频段 expm 的条件数控制高 ω 区间 A·L 的范数可能到几十expm 的误差直接反映在衰减常数上。我常用两个替代一是分段推进每一步只算 A·(L/k) 的矩阵指数再连乘 k 次工程上 L/k 折算成几何段长取 20 到 50 段已经足够二是直接走 3.2 节的特征值对数法计算传播常数彻底避开 expm 的大范数区间。代价是计算量线性上升对 800 个频点完全可接受。5.3 状态向量初值选取与调试技巧计算透射率时左端自由端的状态向量不能全设为零否则只有平凡零解。常见做法是令 Ψ_left 分别为单位矩阵的每一列对每个初始列算右端响应组合成柔度矩阵后按边界条件提取透射系数。这个技巧在调试阶段尤其好用如果 12 列响应出现非物理的对称性破坏说明状态向量排序或符号约定在某一行写反逐列对比能快速定位到具体是哪个自由度出了问题。本文还有配套的精品资源点击获取
