多孔材料吸声系数计算解析解：基于穿孔板阻抗和层间递推

本代码使用解析解方法计算多孔材料的吸声系数，基于穿孔板阻抗和层间递推，并分析了吸声系数随频率的变化趋势。

clc
clear all
f=10:5:1000;
%% 几何参数设置
D = 0.05; % 单胞直径
Tp = 0.01; % 穿孔板厚度
dp = 0.01; % 穿孔板直径
dt = 0.03; % 多孔材料腔体顶部直径
db = 0.02; % 多孔材料腔体底部直径
H = 0.05; % 腔体单胞高度
N = 5; % 分层数量
faipp = dp^2 / D^2; % 穿孔刚性板的穿孔系数
xi = 0.425 * dp * (1 - 1.41 * sqrt(faipp) + 0.34 * sqrt(faipp^3) + 0.07 * sqrt(faipp^5)); % the end correction length

%%% 材料参数设置
epsilonm = 0.89; % 孔隙率
Rf0 = 305767; % 流阻
tau0 = 3.7; % 曲折因子
Lv0 = 2.16e-5; % 粘滞特征长度
Lth0 = 5.36e-5; % 热特征长度
rho0 = 1.2; % air density
c0 = 343;
yeta = 1.84e-5; % dynamic viscosity of air
gamma = 1.4; % 比热率
P0 = 101325; % 大气压
k0 = 0.026; % 空气的导热系数
cp = 1000; % 比热容
Pr = yeta * cp / k0; % 普朗特常数
Z0 = rho0 * c0; % 空气的特征阻抗

%% 分层
for m = 1:N % 层数循环
    di(m) = (dt - db) * (m - 0.5) / N + db; % 每层腔体直径
    Vi(m) = pi * (D^2 - di(m)^2) * H / (4 * N); % 每层多孔材料体积
    faipi(m) = di(m) / D; % 每层的穿孔系数
    hi(m) = H / N; % 每层的高度
end

%% 分层求多孔与空气之间的界面面积
for m = 1:N - 1
    Gi(m) = pi * di(m) * H / N + pi * (di((m+1))^2 - di(m)^2) / 4;
end
Gi(N) = pi * di(N) * H / N;
Lamd_i = 2 * Vi ./ Gi; % 多孔材料的扩散特性长度

%% 频率循环
for n = 1:length(f)
    w = 2 * pi * f(n);
    y = dp * sqrt(rho0 * w / (4 * yeta)); % 穿孔常数

    % 多孔材料基体密度
    rhom1 = tau0 * rho0 / epsilonm;
    rhom2 = epsilonm * Rf0 / (j * w * tau0 * rho0);
    rhom3 = sqrt(1 + j * 4 * w * yeta * rho0 * tau0^2 / (epsilonm^2 * tau0^2 * Lv0^2));
    rhom = rhom1 * (1 + rhom2 * rhom3);

    % 穿孔腔体每层的等效密度
    s = di .* sqrt(w .* rho0 ./ yeta) ./ 2; % 穿孔直径与粘性边界层之比
    rhopi1 = rho0 ./ faipi;
    rhopi2 = (2 ./ (s .* sqrt(-j))) .* ((besselj(1, s .* sqrt(-j))) ./ (besselj(0, s .* sqrt(-j))));
    rhopi = rhopi1 ./ (1 - rhopi2);

    % 多孔材料体积模量
    Km1 = gamma * P0 / epsilonm;
    Km2 = (j * (8 * k0) / (w * Lv0^2 * cp * rho0));
    Km3 = sqrt(1 + j * (w * Lv0^2 * cp * rho0) / 16 / k0);
    Km4 = (gamma - 1) / (1 - Km2 * Km3);
    Km5 = 1 / (gamma - Km4);
    Km = Km1 * Km5;

    % 穿孔腔体每层的等效体积模量
    Kpi1 = gamma * P0 ./ faipi;
    Kpi2 = (gamma - 1) .* (2 ./ (s .* sqrt(-j * Pr))) .* ((besselj(1, s .* sqrt(-j * Pr))) ./ (besselj(0, s .* sqrt(-j * Pr))));
    Kpi = Kpi1 ./ (1 + Kpi2);

    % 计算每层的阻抗和波数
    Ei = D^2 .* (log(1 ./ faipi) - 3/2 + 2 .* faipi - faipi.^2 ./ 2) ./ 16; % 穿孔腔体每层的等效杨氏模量
    wdi = (1 - faipi) .* P0 ./ (epsilonm .* Rf0 .* Ei); % 压力扩散效应的特征频率
    Qi = j .* w .* P0 ./ (epsilonm .* Km .* wdi); % 压力扩散效应的特征波数

    % 这是空气和穿孔孔洞之间界面的静态热渗透率
    Md_i = 8 * Ei ./ (Lamd_i.^2 .* (1 - faipi)); % 空间几何形状的一个量化参数
    Fi = 1 - Qi ./ (Qi + sqrt(1 + Qi .* Md_i ./ 2)); % 在高磁导率对比情况下的频率相关函数
    rhoeffi = 1 ./ (1 ./ rhopi + (1 - faipi) ./ rhom); % 每层的密度
    Keffi = 1 ./ (1 ./ Kpi + Fi .* (1 - faipi) ./ Km); % 每层的体积模量

    %% 计算穿孔板阻抗
    ki = w .* sqrt(rhoeffi ./ Keffi); % 每层的波数
    Zci = sqrt(rhoeffi .* Keffi); % 每层的特征阻抗
    Zpp1 = j .* w .* rho0 * Tp ./ faipp;
    Zpp2 = 1 - (2 * besselj(1, y * sqrt(-j)) / (y * sqrt(-j) * besselj(0, y * sqrt(-j))));
    Zpp3 = 2 * sqrt(2 * yeta * w * rho0) / faipp;
    Zpp4 = 2 * (j .* w .* rho0 .* xi) ./ faipp;
    Zpp = Zpp1 .* Zpp2 + Zpp3 + Zpp4; % 穿孔板阻抗

    %% 阻抗递推
    Zs(1) = -1j * Zci(1) * cot(ki(1) * hi(1));
    for i = 2:N
        Zs(i) = Zci(i) * (Zci(i) - 1j * Zs(i-1) * cot(ki(i) * hi(i))) / (Zs(i-1) - 1j * Zci(i) * cot(ki(i) * hi(i)));
    end
    Zs = Zs(N) + Zpp;

    % 反射系数
    R(n) = (Zs - Z0) / (Zs + Z0);
    % 吸声系数
    alpha(n) = 1 - abs(R(n))^2;
end

%% 画图
plot(f, alpha)
hold on
% 加载数据，画图对比
% load('WPPporousHelmholtz.txt')
% plot(WPPporousHelmholtz(:,1), WPPporousHelmholtz(:,2), 'ro')
% legend('解析解', '数值解') 这段代码debug

运行结果

不加数据对比，结果如下：

吸声系数随频率变化曲线

可以看到，在低频时，吸声系数随着频率的增加而逐渐降低，这是因为在低频时，波长比多孔材料的孔隙尺寸大，声波能够穿透多孔材料，而在高频时，波长比多孔材料的孔隙尺寸小，声波被多孔材料所吸收，因此吸声系数随频率增加而增加。

代码解读

代码主要分为以下几个部分：

参数设置：设置多孔材料的几何参数、材料参数以及空气参数。
分层计算：将多孔材料分为N层，计算每层的几何参数。
频率循环：遍历不同的频率，计算每层材料的等效密度、体积模量、阻抗和波数。
穿孔板阻抗计算：计算穿孔板的阻抗。
阻抗递推：利用递推方法计算多层材料的总阻抗。
吸声系数计算：根据总阻抗和空气的特征阻抗计算吸声系数。
画图：绘制吸声系数随频率的变化曲线。

注意事项

代码中使用的参数需要根据实际情况进行调整。
代码使用解析解方法计算吸声系数，可能会存在一定误差。
代码中使用的公式来自声学理论，需要具备一定的声学基础才能理解。

总结

本代码使用解析解方法计算多孔材料的吸声系数，并分析了吸声系数随频率的变化趋势。该方法可以帮助理解多孔材料的吸声机理，并为多孔材料的设计提供参考。