声子晶体一维透射率计算：传递矩阵法与MATLAB实现

简介：面向声子晶体与振动传递特性研究需求，提供一款MATLAB程序，通过一维传递矩阵法求解周期性复合结构的声波传输特性，重点输出声子晶体的传递率与透射谱，适用于科研人员、研究生及相关工程技术人员。压缩包体积仅7KB，内含1个m脚本文件，单文件实现模型定义、矩阵装配与结果展示，轻量易用。已有246人浏览学习。程序允许自定义单元长度、材料密度与声速、晶体周期数等参数，自动组装总传递矩阵并计算不同频率下的传递率，据描述正确率可达98%，可为声学器件设计、噪声控制及热声制冷等场景提供快速仿真支持。借助该脚本，使用者无需自行推导波动方程和矩阵级联过程，只需调整输入参数即可得到频率-传递率曲线，并据此观察禁带位置、评估隔声性能，为后续实验对比与参数优化提供便利。

1. 声子晶体与传递矩阵法：一维透射率计算，从模型到可跑代码

设计隔振结构时经常遇到这样的场景：钢板厚度加到 20 毫米，某个频段的振动还是稳稳地穿过去。问题往往不在刚度，而在结构本身就是一条畅通的声通道。把两种声阻抗差异足够大的材料按周期交替排列，声波在界面上来回反射，特定频段内会出现类似半导体能带的禁带——这种周期性复合结构就是声子晶体。要预判禁带位置和透射衰减幅度，最省事的工具是一维传递矩阵法。这份 huisang_v13.zip 里的 MATLAB 程序 huisang_v13.m 做的就是这件事：输入材料密度、声速、层厚和周期数，输出频率-传递率曲线，资源说明称正确率可达 98%。适合正在做周期结构隔振、声学超材料设计、结构声学仿真的人拿来做快速验证。

2. 传递矩阵法的数学底子：单层四极子矩阵与总矩阵的拼装

2.1 状态向量与单层矩阵：声波在一层均匀材料里的传播

一维声子晶体可以简化成若干层均匀介质沿 x 方向堆叠，每一层都有自己的密度ρ、声速 c 和厚度 t。声阻抗 Z = ρc，波数 k = ω/c。层内压力场的通解是：

p(x) = A e^{-ikx} + B e^{ikx}

这里 e^{-ikx} 对应沿 +x 方向传播的波，e^{+ikx} 对应沿 -x 方向传播的波。质点振速与压力通过阻抗联系：

u(x) = (A e^{-ikx} - B e^{ikx}) / Z

把 x = 0 处定义为入射面，x = t 处定义为出射面，那么入射面的压力 p(0) = A + B，振速 u(0) = (A - B)/Z，出射面同理。消去 A 和 B，就得到把出射面状态映射回入射面的单层传递矩阵：

function M = tl_matrix(k, t, Z) % 单层均匀介质的传递矩阵 % 输入: k 波数(rad/m), t 层厚(m), Z 声阻抗(Pa*s/m) % 输出: M 把出射面状态[p; u]映射到入射面状态[p; u] M = [ cos(k*t), -1i*Z*sin(k*t); -1i/Z*sin(k*t), cos(k*t) ]; end

矩阵对角项都是 cos(kt)，对应波在层内传播的相位积累；副对角项带虚数单位 i 和阻抗 Z，体现压力和振速两个物理量之间的耦合。这段代码就是整份资源的核心函数，其他所有逻辑都在围绕它做乘法。注意这里的状态向量用的是 [p; u]，也就是压力和质点振速。如果改成 [p; Zu]，矩阵可以写成更漂亮的对称形式，但两种写法结果等价，关键是全程序必须统一约定，否则后面计算透射系数时会出符号错误。

2.2 总矩阵组装与禁带来源：为什么周期性能挡波

有了单层矩阵，剩下的工作就是乘法。若入射端状态向量为 V_in，出射端 V_out，多层结构满足 V_in = M_total · V_out，那么总矩阵由各层矩阵按从入射侧到出射侧的次序左乘得到：

M_total = M_1 · M_2 · … · M_N

对 A/B 交替、周期数为 N 的结构，假设首尾都是 A 层，总矩阵就是：

M_total = (M_A · M_B)^N · M_A

这里有个方向问题：先经过的层，矩阵反而乘在左边。我在 5.1 节专门展开这个坑。禁带出现的数学判据可以用 Bloch 定理解释。对单个周期胞，定义它的传递矩阵 M_cell，则周期性结构中波的传播波数 K 满足：

cos(Kd) = (M_cell(1,1) + M_cell(2,2)) / 2

其中 d = d_a + d_b 是晶格常数。当上式右边落在 [-1, 1] 区间外时，K 变成复数，虚部对应波的指数衰减，这个频带就是禁带。这个判据非常实用，因为它在扫频之前就能预估禁带位置。传递矩阵法比分块有限元省计算量的原因也在这里——扫一个频率点只需要几次 2×2 矩阵乘法，结构再复杂也只是矩阵链变长，不会像有限元那样在网格收敛和稀疏求解上花时间。

2.3 透射率的定义与边界条件：入射波、反射波与透射波的衔接

传递矩阵算出来只是第一步，最终要的是透射率。把入射端声场写成幅值为 1 的入射波加反射波：p_in = 1 + r，u_in = (1 - r)/Z_in；出射端只有沿 +x 方向的透射波：p_out = t，u_out = t/Z_out。代入 V_in = M_total · V_out，得到两个方程：

1 + r = M11·t + M12·t/Z_out
(1 - r)/Z_in = M21·t + M22·t/Z_out

化简成关于反射系数 r 和透射系数 t 的二元线性方程组。这里的透射系数是压力幅值比，功率透射率 τ = |t|²。工程上经常用分贝表示，20·log10(|t|) 和 10·log10(τ) 数值相同，但要想清楚用的是振幅还是功率，严格讲这两种叫法不能混用。入射端和出射端的介质不一定是层材料本身，最常见的情况是两侧都是空气，Z_in = Z_out = 412 rayl。这个边界条件的差异对结果影响很大，我在 5.5 节再讲。

3. huisang_v13.m 核心逻辑拆解：从参数段到频率-透射率谱

3.1 结构定义段：材料表与几何参数先写清楚

拿到这份资源，第一步是找到程序的参数定义段，通常在脚本最顶部。以最常见的实现为例，核心输入是六个量：A 层密度、A 层声速、B 层密度、B 层声速、两层厚度和周期数。代码大致长这样：

% 1) 结构参数段 rho_A = 7800; % A层密度 kg/m^3 (钢) c_A = 5100; % A层声速 m/s (纵波声速) rho_B = 1100; % B层密度 kg/m^3 (环氧树脂) c_B = 1300; % B层声速 m/s d_a = 0.01; % A层厚度 m d_b = 0.005; % B层厚度 m N = 8; % 周期数

声阻抗由密度和声速直接导出：Z_A = rho_A * c_A，Z_B = rho_B * c_B。改代码时不要只想着改厚度，材料换掉时这六个参数要一起换。版本迭代类的程序（比如名字带 v13）经常会把参数段越改越复杂，甚至有从外部 Excel 读参数的版本，但核心还是这六个量。如果看到 MATLAB 脚本里用 input() 提示输入，输入顺序多半也是密度、声速、厚度、周期数。

3.2 单层矩阵函数：别把传递矩阵当黑匣子

上节的 tl_matrix 函数是整段代码的心脏。它的物理含义很直白：对角项是波数乘以厚度的相位余弦，副对角项是阻抗和相位正弦的耦合。程序中可能有名为 transfer_matrix、layer_matrix 之类的函数，内部基本都是这个四极子形式。如果看到指数形式：

M = [exp(1i*k*t) 0; 0 exp(-1i*k*t)]; % 错误示范

这种写法在低频和周期数小时还能用，一旦 k*t 超过几百，exp 的虚部会震荡，数值上不如 cos/sin 稳定。正确写法是 2.1 节那种 cos/sin 形式。若程序中混用了两种形式，结果大概率在高频段出现毛刺或发散。

3.3 扫频循环与总矩阵组装：把每一层拼起来

组装总矩阵的循环是所有代码里最容易写错的地方。我习惯从出射侧往入射侧乘：

% 2) 扫频与总矩阵组装 f = linspace(100, 8000, 4000); % 频率范围与采样点数 tau = zeros(size(f)); Z_A = rho_A * c_A; Z_B = rho_B * c_B; for fi = 1:numel(f) kA = 2*pi*f(fi)/c_A; kB = 2*pi*f(fi)/c_B; MA = tl_matrix(kA, d_a, Z_A); MB = tl_matrix(kB, d_b, Z_B); M_total = eye(2); for n = N:-1:1 % 注意从出射侧往入射侧乘 M_total = M_total * MA * MB; end M_total = M_total * MA; % 若结构以A层收尾，补上最后一层 tau(fi) = transmittance(M_total, Z_A, Z_B); end

循环里 M_total = M_total * MA * MB 这条语句，如果直接写成 for n = 1:N，从入射侧开始乘，方向就会反。对称结构（A/B 交替且首尾同层）恰好掩盖了这个问题，换成非对称结构立刻现原形。所以不管你看到的原始代码怎么写，先核对矩阵乘法次序再跑。频率数组用 linspace 均匀取点，如果禁带容而窄，后面会遇到频点太疏的问题，5.2 节有解法。

3.4 透射系数求解与绘图：二元方程组的直接解法

透射系数的求解函数可以直接按 2.3 节的方程组写：

function tau = transmittance(M, Z_in, Z_out) % 由总传递矩阵和两端介质阻抗求功率透射率 M11 = M(1,1); M12 = M(1,2); M21 = M(2,1); M22 = M(2,2); A = [1, -(M11 + M12/Z_out); 1, Z_in*(M21 + M22/Z_out)]; b = [-1; 1]; x = A \ b; % x(1)=反射系数 r, x(2)=透射系数 t t = x(2); tau = abs(t)^2; end

最后画图用 semilogx 看对数频率轴：

semilogx(f, 20*log10(tau + eps), 'LineWidth', 1.2); xlabel('Frequency (Hz)'); ylabel('Transmission Loss (dB)'); grid on;

这里 tau 是功率比，20·log10(|t|) 和 10·log10(τ) 数值相同，注意在坐标轴标签写明单位是 dB。如果画出来是一条几乎没有起伏的平线，先检查频率上限是不是太低，或者阻抗比设置太小，别急着怀疑矩阵写错。程序里若有附加绘图，通常还会画带隙边界辅助线，那部分主要用于论文排版，日常快速验证用基础曲线就够。

4. 把参数填进去跑起来：材料组合、周期数与透射率曲线的复现

4.1 选材料组合：阻抗比几乎决定一切

用传递矩阵法扫频之前，先做两个判断：声阻抗比够不够大，晶格常数与目标频段匹不匹配。阻抗比 ρ1c1 / ρ2c2 越大，禁带衰减越深，这是物理本质决定的。常见的材料和阻抗参考值如下：

钢和环氧的阻抗比约 28，第一带隙内衰减能做到 -30 dB 以上；钢和铝的阻抗比只有 2.9，带隙很浅，做隔振结构时基本可以放弃这个组合。选材料组合时先查这张表，阻抗比小于 5 的组合，后面周期数堆到很大也补不回来。这是做声子晶体设计时第一眼该看的数。

4.2 周期数与厚度：带隙位置和深度怎么调

晶格常数 d = d_a + d_b 决定带隙中心频率。第一带隙中心大致落在 Bragg 频率附近：

f_B ≈ c_eff / (2d)

其中 c_eff 可以用串联平均粗略估计：2/c_eff = 1/c_A + 1/c_B。周期数 N 决定的是衰减深度，基本不改变带隙宽度。经验关系：N 每翻一倍，带隙内最大衰减大约增加 6 dB。按这个规律可以列一张快速估算表：

厚度调整优先动晶格常数，周期数只做深度微调。例如想把 4 kHz 附近的振动隔掉，先算 d ≈ c_eff / (2f_B)，再分配两层厚度比例。

4.3 三个可直接跑的输入场景

场景一：钢-环氧，d_a = 10 mm，d_b = 5 mm，N = 8。这是最标准的声子晶体配置，阻抗比约 28，带隙中心在 2～5 kHz 范围内，适合验证程序是否跑通。场景二：铝-空气，铝板与空气层交替，阻抗比接近 3 万，带隙极深，但要记住空气层内没有刚性声传播，一维模型只适用于法向入射的平面波，低频段结果会偏离物理实际。场景三：铅-橡胶，低波速材料组合，带隙下移到低频，适合隔振设计，但铅的密度大，实际结构重量要专门评估。

每组配置跑完输出一条透射率曲线，先用 -20 dB 作为判据把禁带边界读出来，再决定要不要做参数微调。

4.4 改代码的三个入口：把脚本改造成可复用函数

拿到手先搜 rho 或 density 找材料参数段，搜 linspace 或 logspace 找频率段，搜 N 找周期数。如果这份资源的脚本里没有参数输入提示，最快的方式是把它包成函数：

function tau = run_v13(f, rho_A, c_A, d_a, rho_B, c_B, d_b, N, Z_in, Z_out) % 一键复用的传递率计算函数 % f 为频率向量(Hz), Z_in/Z_out 为入射端和出射端介质阻抗 % 返回 tau 为功率透射率(与 f 等长) Z_A = rho_A * c_A; Z_B = rho_B * c_B; tau = zeros(size(f)); for fi = 1:numel(f) kA = 2*pi*f(fi)/c_A; kB = 2*pi*f(fi)/c_B; MA = tl_matrix(kA, d_a, Z_A); MB = tl_matrix(kB, d_b, Z_B); M_total = eye(2); for n = N:-1:1 M_total = M_total * MA * MB; end M_total = M_total * MA; tau(fi) = transmittance(M_total, Z_in, Z_out); end end

这个封装把输入端和出射端阻抗单独拎出来，方便做空气耦合修正，后文提到的坑都能从参数入口解决。

5. 避坑手册：矩阵方向、频点疏密与数值稳定性的五个常见翻车点

5.1 矩阵串联顺序反了：透射谱整体失效

现象：程序跑出的曲线像模像样，但把结构改成非对称层序后，结果与参考解相差很大，带隙位置明显偏移。

原因：传递矩阵满足 V_in = M_total · V_out，先经过的层要乘在左侧。写循环时用 for n = 1:N 直接累乘，从入射侧开始编号，结果矩阵方向反了。对称结构恰好掩盖了这个错误——首尾同层时正反方向得到同一个总矩阵。

解决：循环从 N 倒着写到 1，或者把各层矩阵放进 cell 数组后用 fliplr 倒序相乘。我在 3.3 节已经用了这个写法，抄代码时重点核对这一行：

for n = N:-1:1 M_total = M_total * MA * MB; end

5.2 频点太疏：禁带边界被采样漏掉

现象：图上看到窄窄的深谷，边界毛刺明显，把采样点数减半后谷底只剩一个孤点。

原因：linspace 均匀取点，每个频率步长可能跨几百赫兹。第一禁带宽度往往只有几百赫兹甚至更窄，边界点落在两格之间就被跳过了。

解决：改用 logspace，低频既覆盖又加密，或者对目标频段做两段加密采样：

f = logspace(log10(100), log10(8000), 4000);

这样低频段也不会因为步长过大把窄带结构漏掉。跑完先看曲线是否光滑，出现明显折线就先加密频点再谈结果。

5.3 高频段出现 NaN：数值稳定性崩了

现象：周期数 N 较大或频率超过 20 kHz 时，透射率返回 NaN 或 Inf，曲线断裂。

原因：传递矩阵的虚部随 kt 增长震荡，虽然 sin/cos 有界，但当阻抗比极大或两端介质阻抗差异悬殊时，2×2 线性方程组 A\b 的条件数变差，复矩阵求解出现数值发散。改用 e 指数形式的层矩阵会更快崩溃。

解决：先退回双层平均思路降低动态范围，把状态向量从 [p; u] 换成 [p; Zu]，矩阵元素量级接近，条件数大幅改善；还不行就用最小二乘 x = pinv(A)*b 替代左除。另外周期数 N 超过 32 时先把透射率降权处理，或者改用散射矩阵法做级联，哪种都不建议直接硬算。

5.4 声速取错：带隙中心频率漂移 30%

现象：计算带隙中心与实验或参考论文相差两到三成，带宽看着宽，但位置对不上。

原因：声子晶体一维杆模型的纵波声速 c_L = sqrt(E(1-ν)/(ρ(1+ν)(1-2ν)))，很多人顺手填了杆声速 sqrt(E/ρ)，两个值差得远。另外各层必须用各自材料的纵波速度，不能直接用等效声速反过来算单层相位。

解决：先用两种材料的纵波声速分别代入，算一遍带隙中心，再用 Bragg 条件 f_B ≈ c_eff/(2d) 交叉验证。如果程序输出与这个估计差超过 20%，先回材料参数表，而不是怀疑程序本身。

5.5 空气耦合被忽略：隔声量高估成幻觉

现象：算出来带隙内衰减 -60 dB，实测只有 -25 dB，差距大到没法收场。

原因：程序默认入射端和出射端介质就是结构材料本身，Z_in = Z_A、Z_out = Z_B。实际测试时声波从空气进入结构，再回空气，两端介质阻抗是空气的 412 rayl，与固体差两到三个数量级，边界反射把透射能量压掉了很大一部分。模型没改，结果自然成了理想化的结构内部传递率，不是实测意义上的隔声量。

解决：把两端阻抗改成空气：

Z_air = 1.2 * 343; tau = run_v13(f, rho_A, c_A, d_a, rho_B, c_B, d_b, N, Z_air, Z_air);

改完之后再看带隙深度，通常比原来浅 20～30 dB，这才是和实验对得上的数。这是这份资源里最容易被忽略、但对工程判断最有价值的一个参数。

6. 进阶用法：批量扫描填充率与阻抗比，五分钟锁定带隙区间

6.1 批量扫描脚本：把单点计算变成设计扫盲

传递矩阵法最大的好处是算得够快，单次扫频千个点只要秒级，所以可以直接做参数扫描，把填充率 r = d_a/(d_a+d_b) 和阻抗比 Z_ratio 两个设计变量扫一遍，统计每个组合下的带隙起止频率。核心做法是在 4.4 节封装的函数外面再套一层循环：

fractions = 0.1:0.05:0.9; Z_ratios = [10, 20, 40, 80]; bandgap_lo = zeros(numel(fractions), numel(Z_ratios)); bandgap_hi = zeros(numel(fractions), numel(Z_ratios)); for i = 1:numel(fractions) for j = 1:numel(Z_ratios) rho_B = rho_A / sqrt(Z_ratios(j)); % 保持阻抗比 c_B = c_A * sqrt(Z_ratios(j)); d_b = d_a * (1 - fractions(i)) / fractions(i); tau = run_v13(f, rho_A, c_A, d_a, rho_B, c_B, d_b, N, Z_air, Z_air); idx = find(20*log10(tau+eps) < -20); % 以-20dB为带隙阈值 if ~isempty(idx) bandgap_lo(i,j) = f(idx(1)); bandgap_hi(i,j) = f(idx(end)); end end end

这段扫描在解析法下几分钟跑完，得到的是填充率-阻抗比-带隙范围的映射，可以直接画成等值线图 or 色块图，定位最优设计区间。注意阈值 -20 dB 可以按实际隔声需求改，不是固定值。

6.2 带隙边界的验证判据：用 Bragg 估算交叉检查

盲目的参数扫描容易陷入局部最优，先用 Bragg 条件粗筛。第一带隙中心约在 c_eff/(2d)，带隙上下沿大致对称地分布在中心两侧。扫描结果如果带隙中心偏离这个估算超过 20%，大概率是材料参数或矩阵方向出了问题。反过来，代码输出与 Bragg 估算吻合时，再继续做有限元交叉验证才有效率。

6.3 与有限元交叉验证的省时顺序

解析法把 -20 dB 频段圈出来之后，再用 COMSOL 或 ANSYS 在这个窄区间内布点验证，避免开盲盒式仿真。从那以后，我每次换材料组合都会先让传递矩阵跑一遍把频段圈出来，五分钟内看到趋势再上重型软件，既不浪费算力，也不至于漏掉关键禁带。希望帮到你。

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