简介：一份基于Timoshenko梁理论的12×12传递矩阵法声子晶体梁MATLAB计算程序，面向结构动力学、声子晶体及波动控制方向的研究者与学习者。程序将Timoshenko梁的剪切变形与转动惯量纳入12×12矩阵模型，通过周期单元局部传递矩阵的级联构建全局传递矩阵，进而扫描频率并求解声波传播特性，可用于分析带隙位置、频率响应及模态形态，探索几何尺寸、材料参数对声隔离与声过滤效果的影响。压缩包内含1个m文件，整体仅2KB，代码结构简明，适合作为理论验证与二次开发的基础脚本。已有367人学习下载，对于希望快速理解传递矩阵法在声子晶体梁中应用、或需要可运行示例辅助课程设计与科研计算的读者，具有直接参考价值。

1. 用12×12的Timoshenko梁传递矩阵算声子晶体带隙：快、稳、能解释

做周期复合梁的带隙计算时，最常用的是有限元扫频和传递矩阵两条路。有限元能处理复杂截面，但每次调几何参数都要重新分网、求解、提取模态，算一条频散曲线往往要几十分钟起步，调参体验很差。我一般先把周期单胞写成“12×12的Timoshenko梁传递矩阵”，通过Bloch周期条件求特征值，几分钟就能把带隙边界扫清楚，还能顺手把每个频段对应的模态衰减趋势看明白。这个方案对声子晶体梁、压电分流梁、三明治周期梁都适用，适合结构工程师和做减振降噪的研发人员当一个“极速预筛工具”，跑完再决定要不要上有限元复核。

2. 状态向量怎么选：Timoshenko梁的4个分项与12×12扩展矩阵的构造

2.1 Euler梁和Timoshenko梁的边界：深梁和高频谁会翻车

Euler-Bernoulli梁理论假设截面在变形后仍垂直于中性轴，忽略剪切变形和转动惯量。这个假设在“细长梁 + 低频”时误差很小，但在声子晶体梁里往往不成立：为了在低频段打开带隙，单胞厚度通常做得比较大；为了算高阶频带，扫频范围又常常覆盖到数倍于第一带隙的频率。这两个条件叠加后，Euler梁算出的色散关系会有明显偏差，带隙边界可能整体偏移10%以上。我在对比过几组深梁算例后就不再对Euler梁结果做“修正系数”了，直接换Timoshenko梁模型，把横向剪切变形和截面转动惯量都放进控制方程，高频段曲线踏实很多。

Timoshenko梁有两组材料参数是需要显式给定的：弯曲刚度EI和剪切刚度GA。注意这里的G不是弹性模量E，而是剪切模量，并且要乘一个截面剪切系数κ。矩形截面取5/6，圆形截面取9/10。很多初算代码习惯性写成GA=E·A，等于把κ漏掉了，后果是高频带隙位置整体上移，后面第4章会专门讲这个坑。

2.2 场传递矩阵：波数解、剪切系数和材料参数的落点

在频率域里，均匀Timoshenko梁的振动方程可以写成一组关于横向位移w和截面转角φ的常微分方程组。设解为w=a·e^{ikx}、φ=b·e^{ikx}，代入后得到一个关于k²的二次特征方程，解出两个波数k₁²和k₂²，再取正负根就得到四个传播常数。这四个波数对应两对传播/衰减模态，是组装场传递矩阵的基础。

这部分我直接用数值方法处理：给定频率ω、材料参数和截面几何，先用roots解出k²，再代回方程求出每个波数对应的幅值比β=φ/w，然后组装状态向量通解。状态向量我固定取4维顺序：

[w, φ, M, Q]^T

其中M为弯矩，Q为剪力。这里有个符号约定问题：剪力方向如果不一致，后面周期传递矩阵的特征值图像会被“镜像”，看起来像带隙和通带互换，实际是状态向量定义不同导致的。后文避坑章节会再强调。

有了通解形式，把梁段左端x=0和右端x=L的4个状态量分别写出来，就能得到一组线性关系：

q_R = F(L, ω) · q_L

F就是该均匀梁段的4×4场传递矩阵。这个矩阵既不含有限元离散误差，又比直接解析展开简洁，是整条频散曲线的计算基石。

2.3 三层单胞的12×12扩展矩阵组装与内部自由度消元

声子晶体梁的周期单胞通常由两种或三种材料层叠组成。以AB A三层单胞为例，下层和上层是材料A，中间是材料B。按传统做法，直接把三层场矩阵连乘得到4×4单胞传递矩阵：

T_cell = F_A · F_B · F_A

这个式子写法很干净，但工程上有两个别扭的地方：一是如果界面存在粘结层、脱粘或局部刚度削弱，需要在每层之间插入4×4的点传递矩阵，连乘式子会变得难以维护；二是三层结构的中间两个截面是内部自由度，直接连乘时数值误差会在双曲函数项里快速放大，厚单胞或高频段容易算出NaN。

我采用的做法是构造12×12扩展传递矩阵。把单胞左端截面、第一个内部界面、第二个内部界面的状态向量并成一个12维列向量：

v = [q_L; q_int1; q_int2]

三个子段各自满足场传递关系：

q_int1 = F_A · q_L q_int2 = F_B · q_int1 q_R = F_A · q_int2

把这三个关系写成矩阵形式，就得到一个12×12的扩展传递矩阵T_ext。写成矩阵后，可以先用分块高斯消元把内部自由度q_int1和q_int2消掉，得到只含左右端状态的4×4矩阵T_per，再进入带隙特征值计算。这个“先扩展再消元”的做法比直接连乘多写几行代码，但换来了两类好处：中间需要修改界面模型时，直接在12×12的对应分块里插入点矩阵即可；消元过程天然避免了大波数双曲项先相乘后抵消的数值灾难，单胞厚度较大时数值稳定性明显更好。

2.4 一个用于下文代码演示的单胞参数表

为了让后面MATLAB片段可直接复现，我把三层单胞的参数固定如下。这个例子取铝-橡胶-铝，橡胶层提供强阻抗失配，容易在中低频段打开较宽的弯曲波带隙。

梁宽取20 mm，梁高按三层总厚20 mm计算。Timoshenko梁需要给出截面惯性矩I=A·t²/12，其中t是该层厚度；铝层剪切系数取κ=5/6，橡胶层因截面变形较强，也可先按5/6处理，精细建模时可用有限元标定等效κ。

3. 从12×12矩阵到频散曲线：MATLAB实现与参数设置

3.1 周期边界条件与Bloch特征值判定

单胞的传播特性由Bloch周期条件连接左右端状态：

q_R = μ · q_L

μ是传播常数，写成μ=e^{-iqΛ}，其中q是波数，Λ是单胞总长度。当|μ|=1时，波在该频率下可以无衰减通过周期结构，对应通带；当|μ|≠1时，波幅沿单胞方向指数衰减，对应带隙。

具体操作上，对扫频范围内的每个频率点：先算各层场矩阵，组装12×12扩展矩阵，消元得到4×4周期传递矩阵T_per，再解其特征值。特征值的模如果偏离1超过设定阈值，就判定该频率落在带隙内。我自己常用的阈值是1e-6，但更高阶频带里特征值退化，有时需要把阈值放宽到1e-4再做平滑。小于阈值的特征值模、对应的衰减常数以及带隙边界频率，就是最终要输出的三件套。

3.2 MATLAB代码：Timoshenko梁场矩阵与12×12消元

下面这段代码是整套算法的核心。第一个函数负责生成均匀Timoshenko梁段的4×4场传递矩阵，第二个函数负责组装三层单胞的12×12扩展矩阵并消元得到T_per。代码只保留算子部分，便于直接理解逻辑和修改参数。

function F = timoshenko_field(L, rho, E, G, kappa, A, I, omega) % 返回均匀Timoshenko梁的4x4场传递矩阵 % 状态向量顺序: [w; phi; M; Q] % 单位统一: N-m-kg-s lambda2 = rho*A*omega^2 / (kappa*G*A); % 与平动有关的项 beta2 = rho*I*omega^2 / (E*I); % 与转动惯量有关的项 % 波数满足: k^4 - (lambda2 + beta2)*k^2 + lambda2*beta2 - lambda2*k^2 = 0 % 这里直接用幅值比构造，具体展开见推导 p = [1, -(lambda2 + beta2), lambda2*beta2]; r2 = roots([1,0,-p(2),0,p(3)]); r2 = r2(imag(r2)>=0); % 取两个物理相关的k^2 kk = sqrt(r2); % 两个传播常数，另一个方向取负号 k = [kk(1); -kk(1); kk(2); -kk(2)]; betaCoef = zeros(4,1); for j = 1:4 kk0 = k(j); betaCoef(j) = 1i*kk0 / (E*I*kk0^2 + kappa*G*A - rho*I*omega^2) ... * kappa*G*A; % beta = phi/w 的幅值比 end M = @(x) [exp(1i*k*x).'; betaCoef.'.*exp(1i*k*x).']; % M(x) = [w(x), phi(x), M(x), Q(x)] 的基函数矩阵，这里省略M和Q的显式写法 % 实际代码需补全 M(x) 第三、四行，见下方说明 A_in = M(0); A_out = M(L); F = A_out / A_in; % q_R = F * q_L end

上面代码里betaCoef的推导来自Timoshenko梁第二式，把w和φ幅值比约束住；M(x)的第3、4行分别按弯矩M=EI·dφ/dx、剪力Q=κGA·(φ−dw/dx)写出即可，注意符号约定要与状态向量一致。这段代码用矩阵求逆/来得到场矩阵，单段长度不大时没问题；如果单段长度超过数十毫米且频率很高，双曲项增长会让A_in条件数变差，建议改用消元求解线性方程组而不是显式求逆。

下一个代码块展示12×12扩展矩阵的组装与消元过程。

function T_per = assemble_periodic3(F_A, F_B, F_A2) % 三层单胞: A-B-A % 输入: 三个4x4场传递矩阵 % 输出: 4x4单胞周期传递矩阵 T_per T_ext = zeros(12,12); % 块1: q_int1 = F_A * q_L T_ext(1:4, 1:4) = F_A; T_ext(1:4, 5:8) = -eye(4); % 块2: q_int2 = F_B * q_int1 T_ext(5:8, 5:8) = F_B; T_ext(5:8, 9:12) = -eye(4); % 块3: q_R = F_A * q_int2 T_ext(9:12,9:12) = F_A2; % 消去内部自由度: q_int1, q_int2 % 由第1,2块回代，得到 q_R 与 q_L 的直接关系 T_per = F_A2 / T_ext(9:12,9:12) ... * T_ext(9:12,9:12); % 该行仅示意结构 % 实际消元: 先解得 q_int2，再代入 RHS = zeros(4,4); RHS(1:4,1:4)=eye(4); % q_L 单位输入 % 用分块回代计算 T_per，具体代码省略 end

第二个函数里最后几行只是展示“消元”这一关键动作的意图。实际工程代码一般不会对12×12矩阵做完整求逆，而是按三层串联做两次4×4回代：先由q_L得到q_int1，再由q_int1得到q_int2，最后得到q_R，本质上等价于F_A*F_B*F_A。之所以保留12×12的结构写法，是为了在中间任意界面插入点传递矩阵时，只需要在T_ext对应分块上加矩阵而非重排整体逻辑，这个可维护性优势会体现在修改界面条件时。

omega = linspace(10*2*pi, 3000*2*pi, 800); % 10 Hz ~ 3 kHz，用 rad/s muAbs = zeros(size(omega)); kVal = zeros(size(omega)); for i = 1:numel(omega) w = omega(i); F_Al = timoshenko_field(0.005, 2700, 70e9, 70e9/(2*(1+0.33)), 5/6, ...); F_Rub = timoshenko_field(0.010, 1100, 0.1e9, 0.1e9/(2*(1+0.47)), 5/6, ...); T_per = assemble_periodic3(F_Al, F_Rub, F_Al); mu = eig(T_per); [~, idx] = min(abs(abs(mu) - 1)); muAbs(i) = abs(mu(idx)); kVal(i) = angle(mu(idx)) / 0.020; % 单胞长0.02 m end

扫频段里，F_Al和F_Rub的输入参数要按第2.4节表格补全：铝层的截面面积A=0.005×0.02=1e-4 m²，惯性矩I≈8.33e-11 m⁴；橡胶层A=0.01×0.02=2e-4 m²，I≈1.67e-10 m⁴。剪切模量G由E和ν换算。每一步的eig(T_per)特征值模如果取到距离1最近的那个，说明该频率在带隙内投影最小；真正要输出带隙边界，需要把abs(mu)整体画出来，看哪段频率全部大于1。

3.3 扫频参数：范围、步长和带隙判定阈值

扫频范围的经验值是从第一阶弯曲共振的1/10开始，到目标带隙最高边界的1.5倍。比如想设计一个120 Hz附近隔振的周期梁，扫到1500 Hz就足够覆盖前两三阶带隙。步长不是越小越好：步长太密，8000个频率点跑一轮也要几分钟；步长太粗又会漏掉那些只有几赫兹宽的窄带隙。我一般先按对数间隔扫400个点看轮廓，锁定带隙边界附近后再用2 Hz甚至0.5 Hz步长加密，总体计算量受控。

判定阈值方面，|μ|偏离1的量级要做到1e-2以下才算严格带隙。数值误差在全频段都会有轻微影响，需要把“因为算法误差导致的假带隙”排除。一个简单办法是连续观察相邻频率点：真实带隙的|μ|曲线呈平滑抛物线上凸，数值假带隙通常是单点脉冲。后处理的代码里加一个五点滑动平均即可，不必用复杂滤波。

4. 计算声子晶体梁带隙的5个翻车现场：症状、原因与修复

4.1 数值溢出：双曲项爆炸导致矩阵NaN

现象：单胞厚度增大或扫到高频时，传递矩阵元素出现NaN或Inf，频散曲线在某个频率处突然断裂。 原因：Timoshenko场矩阵里包含e^{kL}和e^{-kL}项，k为实数（衰减模态）时，kL较大，双曲函数值暴涨。直接组装后做矩阵减法，大数吃小数，信息丢失，再往后就是NaN。 解决：优先用12×12扩展矩阵分块消元而非显式连乘；长度单位尽量用mm而不是m，降低kL量级；对大波数模态做归一化处理，即每算一段后按最大值缩放到1附近再传播。这属于数值卫生问题，改完通常能撑到更高频段。

4.2 状态向量符号约定不一致，带隙图像被“镜像”

现象：同一组参数，自己代码画出的带隙和文献对不上，通带位置像左右镜像，或者带隙边界频率一模一样但带内衰减特性完全不同。 原因：剪力Q的符号方向、弯矩M的符号方向在不同文献里定义不同。状态向量latex公式看着都一样，代码里正负号差一个负号，最终特征值虚部符号反转，图像就像被翻了一面。 解决：先写一个30元素的单胞算0 Hz附近的第一通带斜率，和Euler梁理论对比。如果斜率一致说明符号正确；如果差一个负号，把Q的定义整体取反。我把这个自检放在开工第一步，比对着文献调半天快得多。

4.3 剪切系数取值错误，带隙边界整体偏移

现象：矩形截面梁的带隙边界比有限元结果偏高约5%~8%，数值上很接近但就是对不上。 原因：剪切刚度写成GA而不是κGA，κ被漏掉。矩形截面κ=5/6，漏掉等效于把剪切刚度放大了20%，高频段剪切效应被低估，色散关系自然偏硬。 解决：检查所有kappa*G*A的传参位置，尤其是橡胶层这类低剪切模量材料，κ的影响会更加显著。另外注意，复合梁如果用等效单层模型，κ不能随意取，需要通过截面剪应力分布反算，否则不如直接分三层建模。

4.4 扫频步长太粗，窄带隙被漏检

现象：粗扫时发现某频段|μ|全部小于1，判断为无带隙；加密后却发现这里有一条仅5 Hz宽的窄带隙，衰减还不小。 原因：带隙窄时，|μ|曲线只在很窄的频率区间内超过阈值，400个点的对数扫频根本踩不到峰值区间。 解决：两段式扫描策略先粗后加密。粗扫用对数间隔400点锁定大轮廓，然后对|μ|>0.98的区域做1 Hz步长直线扫频。这个习惯能让窄带隙不漏检，同时总计算量不失控。

4.5 界面刚度被当成刚性连接，去耦带隙失真

现象：理论预测的低频带隙很宽，实验和有限元却测不到，只有中频段的一个窄带隙勉强对得上。 原因：层间界面不是理想粘结，存在厚度很薄的胶层或微小脱粘。刚性连接假设把界面近似成完美连续，实际结构里上下层的剪切力无法完全传递，低频去耦模态被理论吞掉了。 解决：在12×12扩展矩阵对应位置插入界面点传递矩阵，用界面弹簧模型描述层间刚度。界面刚度可以先取胶层剪切模量除胶厚得到分布刚度，再换算成集中参数。插点矩阵的方法在2.3节的组装结构里改起来非常顺手，这也是我保留12×12扩展形式的一个实际原因。

5. 用有限元和传递率曲线把带隙“验”回来

5.1 有限元Floquet周期边界校验的三步

传递矩阵结果始终是解析模型，截面形状复杂或材料非线性时必须用有限元交叉验证。我常用COMSOL或ANSYS的特征频率分析配合Floquet周期边界：取一个单胞，在左右端面加周期性位移/转角约束，设置扫描参数为波数q沿单胞长度方向扫过第一布里渊区，求解特征频率。三步操作分别是：第一步，几何建模时保证左右端面网格节点一一对应；第二步，在周期边界设置里指定波数q的实部和虚部映射方式；第三步，扫q从0到π/Λ，画出频散曲线并叠加传递矩阵的结果。两者带隙边界对得上，说明解析模型边界条件没设错。

5.2 我给自己定的三条校验习惯

我不太信任不经过交叉验证的“纯解析带隙图”，所以固定给自己三条习惯：第一，新代码写完先跑均匀铝梁，把频散曲线和瑞利-里兹解对照，检查斜率、拐点和截止频率；第二，每次改材料参数前，把旧参数单胞的有限元结果缓存下来，方便对比回归；第三，如果实验室有条件，做一根8到12个单胞的周期梁，一端压电片激励、另一端测响应，传递率曲线在带隙频段的跌落深度一般能超过20 dB，这个实验数据是最终说服自己的证据。

前两条习惯花的时间不超过半天，但能省下后面整整一周的调参焦虑。声子晶体梁的计算本质上是用矩阵把波传播规律理清楚，12×12扩展传递矩阵只是让这个过程更可控、更好改、更不容易数值翻车。工程里没有一劳永逸的公式，只有一套值得长期保留的验证流程，希望这套流程也能帮到你。

本文还有配套的精品资源，点击获取