基于Johnson-Champoux-Allard模型的多孔材料声学性能计算

该代码实现了基于Johnson-Champoux-Allard模型的多孔材料声学性能计算。其中，该模型考虑了穿孔板和多孔材料本身的效应，通过递推求解多层多孔材料的阻抗，进而计算反射系数和吸声系数。具体实现过程如下：

首先，定义了一些几何参数和材料参数，例如单胞直径、穿孔板厚度、穿孔板直径、多孔材料腔体顶部直径、多孔材料腔体底部直径、腔体单胞高度、分层数量、孔隙率、流阻、曲折因子、粘滞特征长度、热特征长度、空气密度、声速等等。然后，将频率进行循环，计算穿孔常数、多孔材料基体密度、每层穿孔腔体的等效密度、多孔材料体积模量、每层穿孔腔体的等效体积模量、每层的阻抗和波数、每层间的界面面积等参数。接着，通过阻抗递推算法，求解多层多孔材料的阻抗，并计算反射系数和吸声系数。最后，使用plot函数画出吸声系数随频率变化的曲线。

需要注意的是，该代码中的一些参数需要根据具体情况进行修改，例如多孔材料的几何形状、孔隙率、材料特性等等。此外，该模型是一种经验模型，仅适用于一定范围内的频率和孔隙率，无法用于预测材料的高频或低频性能。

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