简介针对声子晶体声传播特性分析的一维传递矩阵法MATLAB程序包面向声学、凝聚态物理及结构振动领域的研究人员与学习者用于快速模拟周期性复合结构的声波传播特性。压缩包仅包含单个m文件体积约7KB结构精简便于直接调用与二次修改。程序基于一维传递矩阵法将晶体逐段分解构建总体传递矩阵并结合边界条件求解不同频率下的声波传递率据描述其计算正确率可达98%可作为科研与工程设计的辅助工具。已有246人学习下载使用者可通过修改单元长度、材料密度、声速和周期数等参数获得频率-传递率曲线直观分析禁带位置与传输衰减特性为声屏障、噪声控制及新型声学器件设计提供数据支撑。程序实现包含结构定义、传递矩阵串联、边界条件设置与结果可视化等关键环节适合希望掌握计算原理并快速上手的入门与中级用户。1. 拿到 huisang_v13.zip 先别急着跑传递矩阵法算声子晶体到底在算什么如果你像我一样收到过huisang_v13.zip这种命名接近随机数的压缩包大概率是同行或学生丢来的声子晶体计算程序。这里面的传递矩阵法是声子晶体能带和传递率计算里最常用的一类半解析方法把每一层介质的状态向量沿厚度方向乘起来再用 Bloch 周期条件提取色散关系要算有限周期结构的透射率就把单胞矩阵连乘 N 次。它能解决什么不需要建模、不需要划网格几毫秒扫完一条 0–20 kHz 的传递率曲线特别适合参数扫描找禁带。适合谁做超材料隔声、周期杆/周期板减振、声学超表面的人。第一次打开这个压缩包先别急着点 run花十分钟搞清楚它用的是物理状态基还是幅值基——这两个基的矩阵长得像、符号习惯差得远后面很多翻车都是从这里开始的。2. 传递矩阵法核心表达式从位移连续到特征矩阵的推导与参数选择TMM 在一维分层结构里的地位相当于查表算法频率定下来材料定下来每一层就是一个 2×2 矩阵从头乘到尾边界条件就自动满足了。先别急着往压缩包里找代码把这两条路径在脑子里过一遍后面改参数才不会瞎。2.1 为什么传递矩阵法能在一维声子晶体里横着走一维声子晶体最常见的模型是周期叠层A 层和 B 层交替堆叠纵向波垂直入射。对这种结构每一层内部的波动解就是两个方向的行波叠加不需要像有限元那样把每层再切成一大堆单元。常见做法是用声阻抗 Zρc 和相位 φωd/c 概括整层层与层之间只要满足位移连续和应力连续两个条件就够了。有限元能算任意 2D/3D 结构但做参数扫描时成本感人一个频率一个模型扫 2000 个频点就得解 2000 次非稀疏矩阵。而传递矩阵法扫完一整条传递率曲线只需要一个循环每个频点只做几个 2×2 矩阵乘起来再加一次特征值分解这个成本差异就是它在一维周期结构里“横着走”的原因。局限也很明确它只能吃分层结构和平面波垂直入射这盘菜要算二维周期格子或弹性波斜入射这包就只能当参考了。2.2 层传递矩阵的两种写法物理状态基与幅值基先说物理状态基。取位移 u 和正应力 σ 组成状态向量单层介质内部解是 u(x)Ae^{-jkx}Be^{jkx}应力和位移通过弹性模量 E 关联。定义声阻抗 Zρc、相位 φkd单层从左到右的传递矩阵是[[cos φ, sin φ / Z], [-Z sin φ, cos φ]]这个矩阵行列式为 1两个特征值是 e^{±jφ}分别对应右行波和左行波穿过这一层的相位积累。注意左下角的负号时间因子取 e^{jωt} 还是 e^{-iωt}负号位置会互换。你拿到的老代码十有八九最多只按其中一种约定写成直接照抄很容易把透射率的相位弄反。再说幅值基。每层内只保留前进波和后退波振幅 A、B。穿过一个厚度为 d 的均匀层振幅变化只有相位项传播矩阵写成对角阵P diag(e^{-jφ}, e^{jφ})穿过一个阻抗从 Z1 变到 Z2 的界面由位移连续和应力连续可以解出界面矩阵I(1→2) [[(1r)/2, (1-r)/2], [(1-r)/2, (1r)/2]]其中rZ1/Z2当两测阻抗相等时I 退化成单位阵物理意义正好是“没有反射没有透射损失”。用幅值基算透射率最后只需要取全局矩阵 M 的 M22 的倒数这是最干净的一条路。符号含义常用单位ρ材料密度kg/m3c纵波声速m/sd单层厚度mZρc特征阻抗kg/(m2·s)ω角频率rad/sφkd单层相位厚度radN周期数个2.3 Bloch 条件与能带映射trace(Q) 的临界值就是禁带无限周期结构的核心是 Floquet 定理在周期边界处状态向量只差一个复数倍 λe^{jqa}a 是单胞厚度q 是 Bloch 波数。把幅值基下的单胞矩阵 Q 定义为Q I(B→A) P(B) I(A→B) P(A)这里的坐标取在 A 层起点也就是每个胞刚进入 A 层的位置。Q 把当前胞 A 层起点处的振幅向量映射到下一个胞 A 层起点。因为 Q 是 2×2 且行列式为 1特征值是 λ 和 1/λ于是trace(Q) λ 1/λ 2cos(qa)当|trace(Q)| 2时q 变成复数Bloch 波在传播方向上指数衰减这就是禁带。画色散图时只要把每个频率的 trace(Q)/2 画在纵轴上越过 ±1 的区间就是禁带画传递率曲线时同一频率范围对应透射低谷。这一条能直接验证你手上压缩包里的能带图和传递率曲线对不对。3. 用 huisang_v13.zip 里的核心函数跑通单胞传递矩阵最小复现流程与代码解读理论立住之后动手第一步是认文件第二步是找一个能跑的最小例子第三步是验证输出。别一上来就换成自己的材料参数先用钢/环氧这种经典组合跑一遍出了问题容易定位。3.1 先认文件与输入v13 常见目录结构和主函数入口压缩包解压进去常见做法是先看有没有 README 或注释头。我收到的这类包通常分成三类文件主扫频脚本、单胞矩阵构造函数、画图脚本。文件名不一定叫transfer_matrix还是tmm_main但输入基本逃不出这五个量输入描述常见量纲材料密度数组每层 ρkg/m3声速数组每层 c纵波还是剪切波看题目设定m/s厚度数组每层 dm周期数 N有限结构胞数个扫频范围起止频率和频点数Hz如果压缩包里没有说明文档先找主函数入口里最外层那个 for 循环循环里每轮算一个频率点内部调用的就是单胞矩阵函数。用 v13 这类编号很高的包还有个额外好处量大说明迭代多核心公式大概率是被验证过的坏处是变量名可能改得和论文对不上所以不要只看函数名要看矩阵乘法顺序。3.2 用钢/环氧双层声子晶体跑通传递矩阵核心代码与参数表下面我用 Python 和 numpy 复现最核心的幅值基流程。宿主介质取钢单胞由钢层和环氧层组成N8。代码里保留注释方便对照前面的公式。import numpy as np # 材料与几何输入 rho np.array([7800.0, 1180.0]) # 钢, 环氧树脂 (kg/m3) c np.array([5200.0, 2200.0]) # 纵波声速 (m/s) d np.array([0.02, 0.02]) # 单层厚度 (m) N 8 # 周期数 Z rho * c # 各层特征阻抗 Z0 Z[0] # 宿主介质用钢的阻抗 a d[0] d[1] # 单胞总厚度 # 幅值基基础矩阵 def interface(Z_left, Z_right): r Z_left / Z_right return np.array([[(1 r) / 2, (1 - r) / 2], [(1 - r) / 2, (1 r) / 2]], dtypecomplex) def propagation(phi): return np.diag([np.exp(-1j * phi), np.exp(1j * phi)]) def cell_transfer(omega): I_host_A interface(Z0, Z[0]) # 宿主 - 钢 I_A_B interface(Z[0], Z[1]) # 钢 - 环氧 I_B_A interface(Z[1], Z[0]) # 环氧 - 下一胞的钢 I_B_host interface(Z[1], Z0) # 环氧 - 宿主 P_A propagation(omega * d[0] / c[0]) P_B propagation(omega * d[1] / c[1]) Q I_B_A P_B I_A_B P_A # 单胞 Bloch 矩阵 M_N I_B_host P_B I_A_B P_A np.linalg.matrix_power(Q, N-1) I_host_A return Q, M_N # 扫频 freqs np.linspace(100, 20000, 2000) omega 2 * np.pi * freqs trans_rate np.zeros_like(freqs) gap_flag np.zeros_like(freqs) for i, om in enumerate(omega): Q, M cell_transfer(om) t 1.0 / M[1, 1] # 透射系数 trans_rate[i] np.abs(t) ** 2 # 能量透射率 gap_flag[i] 1.0 if np.abs(np.trace(Q)) 2.0 else 0.0逻辑说明interface按左介质到右介质的顺序拼矩阵这正是前面推导的 I(1→2)。propagation里的-1j对应 e^{jωt} 约定下的右行波相位积累如果你手上的项目用的是 e^{-iωt}把两个指数符号整体互换即可计算|t|^2的结果不变。cell_transfer里 Q 的乘法顺序举足轻重必须先 P_A 再界面再 P_B方向反了 trace 不变但 M_N 会错。最后一个关键点在M[1, 1]M 是把左侧宿主里的振幅向量映射到右侧宿主里右侧没有背向波所以是[A_tr; 0] M [A_in; B_ref]解出来透射系数就是 M22 的倒数。3.3 输出解读色散曲线、禁带区间与传递率低谷对照跑完这段代码把trans_rate画成 dB 曲线把gap_flag画成阴影区间叠在同一个图里。你会看到两个结果在|trace(Q)| 2的频率区间传递率曲线往下掉得又快又深在通带内传递率曲线会来回波动这是有限周期结构两端界面反射造成的驻波波纹不代表有禁带。若在某个小区间gap_flag是 1 但传递率还在 -10 dB 附近先怀疑透射公式取错矩阵元素再怀疑两侧宿主阻抗不一致还硬用 1/M22。判断禁带边界还有一个工程习惯以传递率第一次穿过 -20 dB 的交点为准不要直接把 trace 过 ±2 的点当成绝对边界。有限 N 的禁带边缘本来就有过渡区N 越大过渡区越窄谱上的边缘点向 trace 临界点靠近。所以在 v13 这类包里如果看到能带图和传递率图边界差了几个百分点优先怀疑 N 设太小。4. 参数怎么设把传递率曲线调到想要的位置和深度跑通只是第一步。实际做隔声或减振设计时要回答的问题是“禁带能不能压到我的工作频段”这就要学会调参数。4.1 三个必调参数密度比/声速比、厚度比、周期数周期总厚度 ad1d2 是第一个要动的量。禁带中心频率大体上反比于 aa 缩小一半整段禁带向高频平移一倍左右。常见做法是先定 a 再定频率别单独缩放某一层厚度那样会同时改变周期和占空比曲线变化不好归因。厚度比 d1:d2 是第二个要动的量。保持 a 不变把两层的相位厚度往接近的方向调禁带宽度会变大当两层的光学厚度 ωd/c 接近相等时第一禁带一般做到最宽。这也是为什么很多声子晶体算例都选钢和环氧各一半不是为了材料便宜是为了把禁带调宽。阻抗比 Z1/Z2 是第三个要动的量。密度比和声速比共同决定阻抗对比度对比度越大同一厚度比下的禁带越宽。钢/环氧的阻抗比大约在十几倍出来的第一禁带已经能覆盖接近一个倍频程想更宽就换橡胶或加截面变化但阻抗比过大会让状态基矩阵的条件数变得很难看这时更要坚持用幅值基。参数调整方式主要影响周期厚度 a等比缩放 d1、d2禁带整体平移厚度比 d1:d2固定 a 改占空比禁带宽度、通带波纹阻抗比 Z1/Z2换材料、改截面禁带宽度对比度越大越宽周期数 N增大 N传递率谷深、边界更陡4.2 带阻尼材料怎么写复声速和复模量的落点真实结构里树脂、橡胶这类材料必须有阻尼。常见做法是给阻尼层一个复声速或复弹性模量比如环氧层写成eta 0.02 c[1] c[1] * (1 1j * eta)在 e^{jωt} 约定下c 的虚部取正意味着波数 k 有负虚部右行波 e^{-jkx} 沿 x 衰减这是物理正确的阻尼方向。如果算完发现不加阻尼的系统透射率在某些频段算出了大于 1 的值先查这一行虚部符号而不是去调周期数。加了阻尼后禁带判断就不能再看|trace(Q)| 2了因为 trace 变成复数临界条件模糊掉。更稳的判断是直接看 Q 的特征值模长禁带内最大特征值模长明显大于 1通带内约等于 1。这样也顺便能提取衰减常数留到后面最终验证用。4.3 用扫频结果反推有效参数传递率低谷与 Bloch 衰减线对照一个实用的自检是在扫频循环里顺手记录 Q 的最大特征值模长lam_max np.max(np.abs(np.linalg.eigvals(Q))) alpha np.log(lam_max) / a # 单位: Np/m把 alpha 曲线叠在传递率曲线下面你会发现传递率的深谷和 alpha 的高原严格对齐。若某段 alpha 明明大于 0 而传递率还在 -20 dB 左右晃荡说明有限结构的边界反射把禁带里的衰减“短路”了常见原因是两侧宿主介质选得太软能量从表面漏过去了。这时该改的是宿主介质或入射端边界条件而不是继续加周期数。5. 声子晶体传递率计算的 5 个高频踩坑现象、原因与解决下面这五条是我在类似压缩包代码里反复见到的问题每条都按现象、原因、解决三步写排查起来可以直接照着走。5.1 频率单位混乱色散曲线整体平移 6.28 倍现象是能带图里第一个禁带的中心频率比自己手算的值大了约 2π 倍或者反向光谱出现不存在的“鬼禁带”。原因是扫频循环里用的是 Hz但传递矩阵里的 k 需要的是波数波数里的 ω 必须用角频率。直接把 f 塞进 φ 里等效把声速放大了 6.28 倍。解决统一在循环入口写omega 2*np.pi*f检查单层相位用phi omega*d/c不是f*d/c。5.2 阻抗比悬殊时矩阵条件数爆炸结果出现 NaN 或透射率大于 1现象钢和橡胶这类阻抗差两三数量级的组合状态基写法算到 N10 附近开始出现 NaN或者明明是无源结构透射率却出现几十甚至上百。原因是物理状态基的层矩阵里含有 Z 和 1/Z 作为离对角项相乘几次后动态范围直接冲出浮点精度。解决改用幅值基也就是前面代码里的interfacepropagation写法如果必须保留状态基先把状态向量做阻抗归一化比如换成[√Z·u, σ/√Z]再去做连乘。5.3 直接把含宿主界面的 M1 做 N 次幂透射率曲线出现假波纹现象传递率曲线在通带内多出一组 1/N 周期的细波纹禁带深度反而变浅。原因是有人图省事把单个胞含左右宿主界面的矩阵整体做幂次而实际上相邻胞之间没有宿主层中间只应该重复 Q 的内部界面。解决按标准写法先搭内部单胞矩阵再做matrix_power(Q, N-1)最后只在两端各乘一个宿主界面矩阵。代码里我把这个顺序写死了照着抄基本不会触发这条。5.4 取错 M 矩阵元素禁带位置对但透射深度全错现象gap_flag 标出来的禁带位置和传递率低谷大体能对上但谷底只到 -6 dB 而不是预期的 -40 dB或者反射率曲线长得像透射率。原因是透射系数在幅值基下的公式是t 1/M22不是t M22也不是M11不少人把反射系数和透射系数拿反了。另一个隐藏原因两侧宿主介质阻抗不一致M 已经不满足“右侧无背向波时 M22 的倒数就是透射系数”这个前提。解决先拿一块均质板做 sanity check令 A、B 材料参数完全相同此时透射率应该是|exp(-jφ)|^2如果对不上说明 M 的界面向或取行列式有问题。5.5 阻尼虚部符号反了无源系统算出了能量增益现象加了复声速后高频段透射率反而涨到 1 以上尤其接近周期谐振频率时出现尖峰。原因是时间约定和复模量虚部符号不匹配。e^{jωt} 下 c 的虚部取正代表衰减取负代表介质主动供能e^{-iωt} 下则反过来。解决先在一个单层均质板里试令入射和出射介质完全一样透射率应严格小于等于 1如果大于 1把11j*eta改成1-1j*eta重试。6. 进阶用有限周期透射率曲线反推 Bloch 衰减常数多一条验证手段参数扫完、坑也踩过一轮我习惯再加一道验证把 Q 矩阵的特征值拉开画一条衰减常数曲线用来判断有限周期透射率低谷是“真禁带”还是“两端边界阻抗不匹配的假象”。这条曲线不依赖 N只依赖单胞材料和几何所以它比直接看传递率更深一层。实现上在扫频循环里补三行eigs np.abs(np.linalg.eigvals(Q)) lam_max np.max(eigs) alpha[i] np.log(lam_max) / (d[0] d[1]) # Np/m通带内lam_max约等于 1alpha 接近 0禁带内 lam_max 明显大于 1alpha 变成正值。把 alpha 大于 0 的区间和gap_flag叠在一起两者应该重合。重合不上的时候问题多半出在 Q 的构造顺序而不是材料参数——调材料是白费功夫。这个 alpha 还能用来估算周期数对隔声量的贡献。每增加一个胞禁带中心频率附近的隔离度增量约为20*log10(lam_max)dB。比如用钢/环氧 20mm/20mmN8 和 N16 在禁带中心的传递率差应该约等于 8 倍的这个值。如果实际算出来差太多去检查 M 的朝向而不是怀疑声速数据少了个零。注意这个估算只在离开带边界足够远的频率才准带边界附近有限结构的透射深谷还受两端界面反射影响衰减常数反而不如表面阻抗条件重要。我拿到这类压缩包第一件事就是把 trace(Q) 的临界边界、alpha 曲线和有限周期传递率的 -20 dB 交点叠在一张图上三线对不上就先查矩阵朝向而不是急着调厚度。希望帮到你。本文还有配套的精品资源点击获取