Timoshenko梁12×12传递矩阵与声子晶体带隙计算方法

简介面向声子晶体梁波动特性研究的MATLAB程序包基于Timoshenko梁理论构建12×12传递矩阵。相比欧拉-伯努利梁该理论计入剪切变形与转动刚度适合宽梁、薄壁梁等精细分析场景。程序采用传递矩阵法处理周期结构边界条件可计算频散曲线、带隙与模式形状适合结构声学、周期材料领域的师生和工程师使用。压缩包内共1个.m文件大小仅2KB以纯代码形式给出未附带数据文件便于阅读核心算法并直接运行调试已有367人浏览学习说明该资源在周期结构波动研究中具有一定参考价值。借助这份代码可直观理解12×12传递矩阵的推导与数值实现灵活调整梁的几何、材料参数与边界条件探索不同设计对声波传播带隙的影响为声隔离、声过滤等实际工程问题提供理论分析工具。1. 12×12的Timoshenko梁传递矩阵到底在算什么用12×12的Timoshenko梁传递矩阵去算声子晶体梁的带隙是结构减振工程里性价比最高的一条路子。和单段梁常见的4×4弯曲传递矩阵不同这个12×12矩阵把轴向伸缩、扭转、两个平面内的弯曲全部装进同一个状态向量周期结构两边一拼再套Bloch周期条件就能直接扫出禁带范围。它特别适合两类人一类是要快速评估材料组合和几何尺寸的工程师另一类是想先拿解析模型定好带隙边界、再交给有限元复核的研究者。下面从状态矩阵构造讲起落到能直接跑的周期单元代码最后聊几处我实际踩过的坑。2. 从运动方程到12×12状态矩阵先解决“列什么状态”再写矩阵指数2.1 为什么是12个分量完整梁单元的广义位移与广义力一条轴线沿x方向的等截面直梁左端截面有6个广义位移量轴向位移u、横向位移v和w、绕x轴扭转角θx、绕y轴和z轴的转角θy与θz。与之对应的右端截面有6个广义力/力矩分量轴力N、剪力Vy和Vz、扭矩T、弯矩My和Mz。左端6个位移加6个力恰好12个状态量右端同样12个。所谓传递矩阵就是把这12个左端状态映射到12个右端状态的线性算子。我做程序时状态向量固定排成[ u, v, w, θx, θy, θz, N, Vy, Vz, T, My, Mz ]前6个是位移/转角后6个是力/力矩。两个平面内的弯曲在双轴对称截面下并不耦合所以12×12矩阵实际是四个块对角子矩阵2×2的轴向块、2×2的扭转块、4×4的xy平面弯曲块、4×4的xz平面弯曲块。但写成完整12×12的好处很明显以后截面做成单轴对称或者带偏心质量块弯轴耦合直接往矩阵里填就行不需要改框架。2.2 用状态空间组装Timoshenko梁的动力控制方程Timoshenko梁和Euler-Bernoulli梁最大的区别是它同时考虑剪切变形和截面转动惯量。对谐波振动位移和力都取 e^{iωt} 的时间因子一段均匀梁的控制方程可以写成一阶状态空间形式d z / dx A(ω) · z其中A是12×12的常数矩阵它的非零元素全部与频率ω有关。这样做的好处是不用去背那种四行四列的弯曲传递矩阵闭式解只要把A拼对用矩阵指数 expm 一算就有场传递矩阵程序逻辑非常统一。轴向和扭转各自是独立的两块。xy平面内的弯曲v、θz、Vy、Mz和xz平面内的弯曲w、θy、Vz、My形式对称只差惯性矩对应的方向。这里最容易出错的是转角符号我采用的约定是 θz dv/dx 的弹性转角贡献而 θy -dw/dx 的弹性转角贡献截面剪力方向按右手坐标系补上对应的剪切角。2.3 场传递矩阵的数值实现与量纲检查下面这段代码直接构造12×12的状态矩阵A。材料参数用字典传进去截面几何参数也一并处理import numpy as np from scipy.linalg import expm def state_matrix(omega, p): 构造Timoshenko梁的12x12状态矩阵A(omega) 状态向量顺序 [u, v, w, thx, thy, thz, N, Vy, Vz, T, My, Mz] p 需要字段 rho, E, G, kappa, A, Iy, Iz, Ip rho p[rho] E p[E] G p[G] kappa p[kappa] A_cross p[A] Iy p[Iy] Iz p[Iz] Ip p[Ip] A np.zeros((12, 12)) # 轴向u N/(EA)N -rho*A*w^2*u A[0, 6] 1.0 / (E * A_cross) A[6, 0] -rho * A_cross * omega**2 # 扭转thx T/(G*Ip)T -rho*Ip*w^2*thx A[3, 9] 1.0 / (G * Ip) A[9, 3] -rho * Ip * omega**2 # xy平面弯曲v, thz, Vy, Mz A[1, 5] 1.0 # v thz Vy/(kappa*G*A) A[1, 7] 1.0 / (kappa * G * A_cross) A[5, 11] 1.0 / (E * Iz) # thz Mz/(E*Iz) A[7, 1] -rho * A_cross * omega**2 # Vy -rho*A*w^2*v A[11, 5] -rho * Iz * omega**2 # Mz Vy - rho*Iz*w^2*thz A[11, 7] 1.0 # xz平面弯曲w, thy, Vz, My A[2, 4] -1.0 # w -thy Vz/(kappa*G*A) A[2, 8] 1.0 / (kappa * G * A_cross) A[4, 10] 1.0 / (E * Iy) # thy My/(E*Iy) A[8, 2] -rho * A_cross * omega**2 # Vz -rho*A*w^2*w A[10, 4] -rho * Iy * omega**2 # My Vz - rho*Iy*w^2*thy A[10, 8] 1.0 return A def field_matrix(omega, p, L): 场传递矩阵z_R F * z_L return expm(state_matrix(omega, p) * L)关于参数有一点要特别提醒这里的A是截面积别和状态矩阵A搞混Iy和Iz分别是绕y轴和z轴的截面惯性矩Ip是极惯性矩。矩形截面b×h若b沿y向、h沿z向则 Iz bh³/12Iy b³h/12Ip Iy Iz。量纲上A矩阵乘以长度L之后是纯数expm才能给出正确的无量纲化矩阵指数。矩阵指数在数值上等于把一段均匀梁的精确闭式解全部展开所以它天然包含了轴向、扭转、弯曲三个波族在高频下的所有信息。使用它之前建议做一个自检对单段均匀梁令周期单元长度为L直接把F的特征值取对数得到三个方向的波数应当和解析色散关系一致。我通常拿这个当冒烟测试能通过就说明A矩阵的符号和位置基本没写错。3. 拼出声子晶体梁的周期单元Bloch条件与带隙计算3.1 两种材料拼接是声子晶体梁的最小可复现模型声子晶体梁能出现弯曲波带隙最常见的是布拉格散射机制周期单元里至少有两种不同波阻抗的材料段弹性波在界面反复反射某些频率段内没有实数波数解波就传不过去。工程上最常见也最好复现的模型就是一段材料A、一段材料B交替排列。周期单元从第j段左端传到第j段右端总传递矩阵是两段场传递矩阵按波的传播方向相乘。注意矩阵乘法的顺序先过A段再过B段那么单元左端状态先被F_A作用再被F_B作用所以T F_B F_A很多新手在这里把顺序写反结果算出来的带隙位置差得离谱而且是那种“看起来挺合理但和有限元对不上”的差。def unit_cell_matrix(omega, pA, pB, LA, LB): 周期单元传递矩阵先A段后B段 FA field_matrix(omega, pA, LA) FB field_matrix(omega, pB, LB) return FB FA这里LA和LB是两个材料段的长度单元总长Lcell LA LB。两段的截面几何可以相同只换材料参数这是最常见的设计也可以同时换截面尺寸那就变成几何型声子晶体。代码层面没有区别A矩阵里对应参数填进去即可。3.2 特征值扫描传播常数怎么从传递矩阵里提出来周期结构满足Bloch条件单元右端状态等于左端状态乘一个复常数λλ e^{iqLcell}其中q是波数。代入 z_R T · z_L得到特征值问题( T - λ·I ) · z_L 0于是问题变成对每个频率求T的特征值。无阻尼情况下传递矩阵是辛矩阵特征值成对出现λ和1/λ。如果某个特征值落在单位圆上也就是|λ|1说明该频率有纯实数波数波能通过对应通带如果出现成对的实数特征值一个|λ|1、一个|λ|1说明波在周期单元内指数衰减对应带隙。实际扫频时我一般不用严格判断|λ|是否等于1而是看对数波数的虚部。把λ取复对数再除以iLcell得到q -i·log(λ)/Lcell。通带里q接近实数带隙里q有实部也有虚部虚部就是每单位长度的衰减系数。def scan_band(freqs, pA, pB, LA, LB): 扫频并返回通带频率-波数数据与带隙衰减数据 freqs: 一维数组单位Hz pA/pB: 材料参数字典 LA/LB: 两段长度 Lcell LA LB bands [] atten [] for f in freqs: omega 2.0 * np.pi * f T unit_cell_matrix(omega, pA, pB, LA, LB) lam np.linalg.eigvals(T) for lam_i in lam: qL np.log(lam_i) / 1j # 无量纲波数 q*Lcell if abs(qL.imag) 1e-5: # 实波数落在通带 q_real np.angle(lam_i) # 限制在[-pi, pi] bands.append((f, q_real)) else: atten.append((f, abs(qL.real))) return bands, atten对每个频率12×12矩阵给出12个特征值对应6个正向传播模态和6个反向传播模态。所以画出来的色散图不是一条曲线而是好几条同时存在轴向波、扭转波、两个平面内的弯曲波各占一枝。扫频时建议先把频率上限设到目标带隙频率的1.5倍左右步长取小一点先粗扫一遍观察趋势再对感兴趣频段细扫。粗扫结果里带隙最直观的表现是色散曲线在某段频率内完全断开同时衰减数据里出现明显大于零的平台。我的习惯是先看atten里有没有连续凸起有凸起再去bands里核对断带位置两段数据对应上才算带隙可靠。4. 带隙对参数的敏感度剪切系数、长度占比与材料失配4.1 Timoshenko梁在多高频率必须被认真对待很多教材里的声子晶体梁例子都用Euler-Bernoulli梁的4×4传递矩阵因为推导简单写起来也快。但工程实际里的梁尤其是用于减振的短粗梁、深梁剪切变形影响并不小。到底什么时候Timoshenko模型不能省我自己的判断标准是看截面高度与弯曲波长的比值。当梁高超过波长的十分之一Euler-Bernoulli理论就会明显高估弯曲波速带隙边界算出来偏高越高频偏得越厉害。声子晶体梁要抑制的往往正是中高频振动恰好落在Timoshenko效应不可忽略的区间。所以在12×12矩阵里剪切修正系数κ、截面惯性矩Iy和Iz都必须认真给这直接决定弯曲弯带的带隙边界。轴向和扭转模态不受κ影响但那些模式往往不是你最关心的。4.2 用参数扫描把带隙调向目标频段带隙工程本质上就是调两个东西材料失配程度和单元内两段长度的占比。材料失配决定带隙宽度。密度比和弹性模量比越悬殊界面反射越强带隙越宽。工程上最常见的组合是金属加聚合物或金属加橡胶。金属段提供刚度聚合物段提供质量对比。比如铝-环氧组合E和ρ都有好几倍差距做一维梁周期结构很容易在几千赫兹附近得到可观的弯曲波带隙。长度占比决定带隙中心频率和带宽的权衡。设r LA / (LALB)当r接近0.5时两段几何上等长布拉格反射最强第一带隙通常最宽偏离0.5后中心频率会移动带宽通常收窄。如果目标频段偏低就整体加大单元长度如果目标频段偏高整体缩小单元。这个规律在周期性梁和杆上都成立属于布拉格散射的基本特征。实际做参数扫描时我一般把r从0.1扫到0.9步长0.05材料不变对每个r跑一遍扫频记录第一带隙起止频率。然后画一条“带隙宽度随r变化”的曲线曲线的峰值位置就是几何上最优的占比。计算很快几分钟就能出结果比直接上有限元试错强太多。下表是一组典型低碳钢-环氧组合的示例参数实际算的时候照着填就行。参数钢段环氧段密度 ρ7850 kg/m³1180 kg/m³弹性模量 E206 GPa3.0 GPa泊松比 ν0.300.38剪切模量 GE/(2(1ν))E/(2(1ν))截面宽 b0.02 m0.02 m截面高 h0.02 m0.02 m剪切修正系数 κ5/65/6需要提醒的是钢和环氧的泊松比差别虽然不大但G的差异会直接影响Timoshenko梁的剪切刚度导致弯曲模态的高频色散曲线形态不同。如果只改E和ρ不改G带隙宽度会算不准这是参数扫描里经常踩的坑。4.3 剪切修正系数的选取与量化影响剪切修正系数κ不是无量纲的“经验凑数”它取决于截面剪应力分布。矩形实心截面取5/6圆形实心截面取6/7这是最常用的两个值。但如果是薄壁矩形管、工字钢这类截面κ就不能随便套得按Cowper公式算或者用截面有限元做一次静力剪切分析标定。κ对带隙的影响随频率增大而放大。低频段第一带隙内弯曲波长长剪切变形占比小κ给5/6还是1.0带隙边界可能只差百分之几到高阶带隙或者高频窄带κ的影响能到百分之十几这时候用错κ算出来的带隙和实验对不上基本就找不到原因了。我吃过这个亏一开始图省事把κ全部设成1.0结果高频段色散曲线和有限元差了将近15%后来查了一圈才发现是剪切修正系数没按截面取。另外如果梁是层合结构或者泡沫芯夹层等效剪切刚度不能用单一材料G简单计算。这时候最稳妥的做法是先做等效单层参数标定再进传递矩阵。等效参数不准后面一切都是白算。5. 避坑12×12传递矩阵最常见的五个翻车点5.1 矩阵乘错顺序导致带隙左右颠倒现象扫出来的带隙区段在频率轴上像是被镜像了本该低频禁带的地方出现通带本该通带的地方出现禁带或者带隙宽度看起来差不多位置完全不对。原因周期单元总传递矩阵的乘法顺序写反。如果写成F_A F_B而不是F_B F_A等于把单元内的材料顺序颠倒了。对于两种材料组成的单元顺序颠倒后表面看起来还是同一组材料但界面反射相位关系变掉带隙位置自然偏移。解决先明确传播方向状态向量从单元左端传到右端先经过哪段就先乘哪段。我的检查方法是把LA和LB设成相等、材料也设成相同此时T应该等于单段长度为LALB的场传递矩阵。如果不是顺序肯定错了。5.2 高频段矩阵指数数值病态现象频率扫到几十千赫兹以上特征值里出现1e15量级的大数和1e-15量级的小数通带判断完全失真衰减数据杂乱无章。原因expm(A·L)在L较大或ω较高时矩阵指数里的双曲函数项指数增长矩阵条件数急剧恶化特征值分解的数值误差被放大。解决一是在构造F之前把单元内每段再细分成若干子段用F_sub F_sub ... 代替整段F这相当于数值上给矩阵指数“降火”二是改用阻抗矩阵或散射矩阵格式做周期单元装配辛结构更稳定三是扫描频率上限不要超过实际关心的范围高频段不是声子晶体梁的典型工作区没必要硬算。5.3 剪切修正系数选错高频带隙偏掉百分之十几现象低频第一带隙和有限元对得上第二、第三带隙边界明显偏高而且差距随频率增大。原因κ取值没有按截面形状选或者用了各向同性材料的矩形截面公式去套复合材料截面。解决矩形实心截面取5/6实心圆取6/7薄壁截面按Cowper公式计算。如果截面内有多种材料先做等效剪切刚度标定再反算κ。高频带隙越重要这一步越不能省。5.4 Bloch条件符号取反导致衰减曲线缺失现象通带色散曲线正常但带隙区域看不到明显的衰减凸起衰减数据全是零或者全是乱值。原因Bloch条件里的传播方向符号写反。如果用了 z_L e^{iqL} z_R 而不是 z_R e^{iqL} z_L特征值互为倒数但取对数后实部和虚部符号会互换衰减信息被错误归零。解决固定约定 z_R T z_L特征值λ对应左端态到右端态。用单段均匀梁验证均匀梁没有带隙所有λ都应该在单位圆上如果出现|λ|明显偏离1符号约定就有问题。5.5 只用4×4弯曲矩阵漏掉轴向和扭转通带现象设计目标是抑制弯曲振动结果实验里某个频率下端部振动还是大。一查传递矩阵只算了弯曲波轴向波在某段频率里其实是通带能量走了轴通路。原因声子晶体的带隙往往是频带特定的。4×4弯曲传递矩阵只能给出弯曲波带隙轴向波和扭转波带隙需要单独算。如果结构存在弯轴耦合路径局部激振能量会通过轴向模态泄漏出去。解决直接用12×12传递矩阵把所有波族同时算出来。判断带隙时要看波族弯曲带隙、轴向带隙、扭转带隙可能在不同频段出现。想衰减哪个模态就盯着对应那一枝色散曲线设计单元别拿弯曲带隙去覆盖轴向问题。6. 用辛性条件自检矩阵再拿有限元对一次带隙6.1 辛性条件不用实验也能快速验证传递矩阵写没写错无阻尼情况下传递矩阵是辛矩阵满足 T^T · J · T J其中J是12×12的分块反对称矩阵。这个条件非常严格状态向量排序、符号约定、A矩阵任何一处写错辛性检查都会失败。我每次写完代码第一件事就是跑这个自检比对着解析解找错快得多。def check_symplectic(T, tol1e-8): 检查传递矩阵是否满足辛性条件 J np.block([[np.zeros((6, 6)), np.eye(6)], [-np.eye(6), np.zeros((6, 6))]]) residual T.T J T - J return np.max(np.abs(residual)) tol如果残差大于1e-6先检查状态向量排序是不是前6个位移、后6个力再检查A矩阵里v、w两个弯曲块的符号对称性。因为辛性条件对符号极其敏感它能把“看起来合理但就是不对”的矩阵直接揪出来省掉大量排错时间。6.2 和有限元对拢的标准流程传递矩阵算完至少要和有限元对一个频段再往下走。标准流程是只建一个周期单元的几何模型两端截面加Floquet周期边界条件波数q从0扫到π/Lcell求解特征频率。然后把各q下的频率点画成曲线和传递矩阵扫出来的色散曲线叠在一起。两条线在前三阶模态内误差小于3%到5%基本可以确认建模正确。这里有个细节有限元梁单元如果用Euler-Bernoulli类型高频段会天然和Timoshenko结果偏离那不是你声子晶体模型错了是有限元单元理论选错了。比较时有限元模型要么用实体单元要么用Timoshenko梁单元否则两者差异说不清楚。实体单元要注意网格在截面厚度方向至少划分两层才能反映出剪切变形。我自己的习惯是先跑辛性自检再跑单段均匀梁的解析色散对照最后才和有限元对周期单元。前两步都过了有限元对拢基本一次成功。12×12传递矩阵这套东西最难的不是算是把符号约定和材料参数搞对这两点把住了声子晶体梁的带隙设计就是一条很顺的流水线。希望帮到你。本文还有配套的精品资源点击获取