二维局域共振声子晶体超材料带隙计算matlab程序

二维局域共振声子晶体超材料的带隙计算可以采用有限元或者传输矩阵方法，以下是MATLAB代码实现传输矩阵方法的带隙计算。

首先定义声波在晶体中的传播方向，假设声波沿着x方向传播，则波矢量可以表示为kx，ky=0。

接下来定义传输矩阵的参数，包括声速c，密度rho，晶格常数a，单元大小d，和板厚度h。

c = 343; % m/s
rho = 1.225; % kg/m^3
a = 0.05; % m
d = 0.02; % m
h = 0.01; % m

然后定义晶体的布局，包括晶格的大小和形状，以及晶体中单元的数量。

n = 10; % number of unit cells
N = n*2+1; % number of nodes
L = a*n*2; % length of crystal

接下来计算晶体中每个单元的声速和密度，以及每个单元的传输矩阵。

for i = 1:n
    % acoustic properties and transfer matrix for unit cell
    c1 = c*(1+0.1*(i-n-1)); % speed of sound
    rho1 = rho*(1+0.1*(i-n-1))^(-1); % density
    k1 = 2*pi*c1/a; % wavevector
    Z1 = rho1*c1; % acoustic impedance
    M1 = [1, 1/Z1; Z1, 1]; % transfer matrix
    T1 = expm(1i*k1*d*M1*h); % transfer matrix for unit cell
    T(:,:,i) = T1;
end

接下来计算整个晶体的传输矩阵，以及计算带隙和色散关系。

% transfer matrix for entire crystal
Ttotal = T(:,:,n);
for i = n-1:-1:1
    Ttotal = T(:,:,i)*Ttotal;
end
for i = n+1:N-1
    Ttotal = T(:,:,i-n)*Ttotal;
end

% band structure
kmax = 2*pi/L;
k = linspace(0, kmax, 500);
w = zeros(size(k));
for i = 1:length(k)
    % frequency and wavevector
    M = Ttotal(1:2,1:2);
    w(i) = c*k(i)*sqrt(det(M));
end
plot(k*a,w/2/pi);
xlabel('ka');
ylabel('f (Hz)');

以上是MATLAB程序的基本框架，可以根据具体的问题进行修改和调整。