当前位置：首页 > news >正文

MATLAB 实现一维声子晶体相关计算：从传递矩阵法到能带图、响应图与弥散关系

news 2026/5/11 13:14:10

MATLAB实现传递矩阵法计算一维声子晶体能带图，响应图，弥散关系计算程序传递矩阵法一维声子晶体超材料声子晶体能带图计算

最近在研究超材料领域中一维声子晶体相关内容，发现通过传递矩阵法结合 MATLAB 能高效地计算其能带图、响应图以及弥散关系，今天就来和大家分享一下实现过程。

传递矩阵法原理简介

传递矩阵法是求解一维声子晶体问题的常用方法。对于由不同材料交替组成的一维声子晶体结构，我们可以将其看成是多个具有不同声学特性的单元层的组合。在每一个单元层内，波动方程具有特定的解形式，通过在不同层之间的边界条件匹配，可以得到从一层到另一层的传递关系，这种关系就可以用传递矩阵来描述。

假设一维声子晶体由两种材料 A 和 B 交替排列组成，对于沿 x 方向传播的弹性波，在第 j 层的位移场和应力场可以用一个 2x2 的传递矩阵 $Mj$ 来联系相邻界面的状态。整个结构的总传递矩阵 $M$ 就是各层传递矩阵的连乘，即 $M = \prod{j = 1}^{N} M_j$，其中 N 是结构的总层数。

MATLAB 实现能带图计算

代码部分

% 参数设置 a = 1; % 晶格常数 fmin = 0; % 最小频率 fmax = 10; % 最大频率 npoints = 100; % 频率点数 kvec = linspace(0, pi/a, 100); % 波矢范围 freq = linspace(fmin, fmax, npoints); % 材料参数 rho1 = 1; % 材料1密度 rho2 = 2; % 材料2密度 c1 = 1; % 材料1声速 c2 = 1.5; % 材料2声速 l1 = 0.4 * a; % 材料1层厚度 l2 = 0.6 * a; % 材料2层厚度 band = zeros(length(kvec), length(freq)); for ik = 1:length(kvec) k = kvec(ik); for ifreq = 1:length(freq) omega = 2 * pi * freq(ifreq); M1 = [cos(k * l1) 1i / (rho1 * omega * c1) * sin(k * l1); -1i * rho1 * omega * c1 * sin(k * l1) cos(k * l1)]; M2 = [cos(k * l2) 1i / (rho2 * omega * c2) * sin(k * l2); -1i * rho2 * omega * c2 * sin(k * l2) cos(k * l2)]; M = M1 * M2; % 这里简单假设两层为一个周期 eigval = eig(M); cond = abs(eigval(1)) <= 1 && abs(eigval(2)) <= 1; if cond band(ik, ifreq) = 1; end end end figure; [X, Y] = meshgrid(kvec * a / pi, freq); contourf(X, Y, band, 'LineStyle', 'none'); xlabel('k*a/\pi'); ylabel('Frequency (Hz)'); title('Band Structure of 1D Phononic Crystal'); colorbar;

代码分析

参数设置部分：首先定义了晶格常数a，频率范围fmin到fmax以及频率和波矢的采样点数npoints。同时设定了两种材料的密度rho1、rho2，声速c1、c2以及它们各自层的厚度l1、l2。这些参数是后续计算的基础，不同的材料特性和结构参数会导致不同的能带图结果。
初始化和循环部分：初始化一个二维矩阵band来存储能带信息。通过两层循环，外层遍历波矢kvec，内层遍历频率freq。在每一步循环中，计算当前波矢和频率下的传递矩阵M1和M2，这里根据波动方程在不同材料层中的解来构建传递矩阵。
判断条件和绘图部分：计算总传递矩阵M的特征值eigval，根据特征值的模是否都小于等于 1 来判断该频率和波矢组合是否处于能带内（这是基于传递矩阵法的能带判断条件）。如果满足条件，就在band矩阵对应位置标记为 1。最后使用meshgrid和contourf函数绘制能带图，横坐标为归一化波矢，纵坐标为频率。

响应图计算

代码部分

% 继续使用上述材料和结构参数 omega = 2 * pi * linspace(fmin, fmax, npoints); response = zeros(length(omega), 1); for i = 1:length(omega) M1 = [cos(omega(i) / c1 * l1) 1i / (rho1 * omega(i)) * sin(omega(i) / c1 * l1); -1i * rho1 * omega(i) * sin(omega(i) / c1 * l1) cos(omega(i) / c1 * l1)]; M2 = [cos(omega(i) / c2 * l2) 1i / (rho2 * omega(i)) * sin(omega(i) / c2 * l2); -1i * rho2 * omega(i) * sin(omega(i) / c2 * l2) cos(omega(i) / c2 * l2)]; M = M1 * M2; response(i) = abs(M(1, 1)); end figure; plot(linspace(fmin, fmax, npoints), response); xlabel('Frequency (Hz)'); ylabel('Response Magnitude'); title('Response Diagram of 1D Phononic Crystal');

代码分析

参数和初始化：继续使用前面定义的材料和结构参数，定义频率范围omega并初始化响应向量response。
循环计算：在频率循环中，同样构建每一层的传递矩阵M1和M2并得到总传递矩阵M。这里取传递矩阵M的 (1, 1) 元素的模作为该频率下的响应幅值，这是一种常见的对结构响应的度量方式。不同的响应定义可能会根据具体研究目的而改变。
绘图：最后使用plot函数绘制频率与响应幅值的关系图，即响应图，直观展示结构在不同频率下的响应特性。