薄膜声学超材料模态叠加法Matlab实现(Zhang Yuguang论文复现)

机械工人杨师傅

已于 2022-11-27 16:51:01 修改

阅读量1.1k

点赞数 5

分类专栏：隔声论文复现文章标签： matlab

于 2022-11-26 19:43:20 首次发布

本文链接：https://blog.csdn.net/YGZ11113/article/details/128048546

版权

隔声同时被 2 个专栏收录

7 篇文章 10 订阅

订阅专栏

论文复现

7 篇文章 2 订阅

订阅专栏

1.理论分析

设超材料的一个周期单元的长度为Lx、Ly；薄膜密度为 $\rho_s$ ，薄膜张紧力为T；振子尺寸为lx、ly，振子密度为 $\rho_{mass}$ ，振子位置为 $x_0,y_0$ ；板的横向位移为w

$w(x,y,t)=\sum_{m=1,n=1}^{MN}W_{mn}(x,y)q_{mn}(t)$

定义m、n分别为x、y方向的模态数，模态型函数为：

$W_{mn}(x,y)=sin(\frac{m\pi x}{Lx})sin(\frac{m\pi y}{Ly})$

列出声波入射在板上引起的振动微分方程

$\rho_{s}\frac{\partial ^{2}w}{\partial t^{2}}+\rho_{mass}h(x,y,x_0,y_0,l_x,l_y)\frac{\partial ^{2}w}{\partial t^{2}}+2\rho_{0}c_{0}\frac{\partial w}{\partial t}-T\nabla^2w=2P_{in}e^{jwt}$

其中 $\rho_{0},c_{0}$ 为空气密度， $P_{in}$ 为入射声波幅值；

根据振动微分方程求解得到以下等式

$-\omega^{2}\left \{ \left [\ \mathbf{M} \right ] + \left [\ \mathbf{Q} \right ] \right \}\left \{ \mathbf{\widetilde{q}} \right \}+j\omega \left [\ \mathbf{Q} \right ] \left \{ \mathbf{\widetilde{q}} \right \}+ \left [\ \mathbf{K} \right ] \left \{ \mathbf{\widetilde{q}} \right \}=2P_{in}\left \{ \mathbf{H} \right \}$

各矩阵的求解方式如下：

$M_{mn}=\rho_s \int_{0}^{L_x}\int_{0}^{L_y}W_{mn}\sum_{r=1,s=1}^{M*N}W_{rs}dxdy=\left\{\begin{matrix} 0,m \neq r || n\neq s\\ \frac{L_xL_y}{4},m=r \And\And n==s \end{matrix}\right.$

$I_{mn,rs}= \int_{x_0}^{x_0+l_x}\int_{y_0}^{y_0+l_y}W_{mn}W_{rs}dxdy$

$C_{mn}=2\rho_0c_0 \int_{0}^{L_x}\int_{0}^{L_y}W_{mn}\sum_{r=1,s=1}^{M*N}W_{rs}dxdy$

$K_{mn}=-T \int_{0}^{L_x}\int_{0}^{L_y}W_{mn}\nabla^2\sum_{r=1,s=1}^{M*N}W_{rs}dxdy=\left\{\begin{matrix} 0,m \neq r|| n\neq s\\ \frac{L_xL_y}{4}\left [ \left ( \frac{r\pi}{L_x} \right )^2+\left ( \frac{s\pi}{L_y} \right )^2 \right ],m=r\And\And n==s \end{matrix}\right.$

$H_{mn}= \int_{0}^{L_x}\int_{0}^{L_y}W_{mn}dxdy$

$\mathbf{M}=\begin{bmatrix} M_{11} &0 & 0 & 0\\ 0& M_{12}& 0 &0 \\ \vdots & \vdots& \vdots &\vdots \\ 0& 0& 0& M_{MN} \end{bmatrix}=\mathbf{diag}\begin{bmatrix} M_{11} & \cdots & M_{M1}&M_{12} & \cdots& M_{M2} & \cdots & M_{MN} \end{bmatrix}$ $\mathbf{K}=\begin{bmatrix} K_{11} &0 & 0 & 0\\ 0& K_{12}& 0 &0 \\ \vdots & \vdots& \vdots &\vdots \\ 0& 0& 0& K_{MN} \end{bmatrix}=\mathbf{diag}\begin{bmatrix} K_{11} & \cdots & K_{M1}&K_{12} & \cdots& K_{M2} & \cdots & K_{MN} \end{bmatrix}$ $\mathbf{C}=\begin{bmatrix} C_{11} &0 & 0 & 0\\ 0& C_{12}& 0 &0 \\ \vdots & \vdots& \vdots &\vdots \\ 0& 0& 0& C_{MN} \end{bmatrix}=\mathbf{diag}\begin{bmatrix} C_{11} & \cdots & C_{M1}&C_{12} & \cdots& C_{M2} & \cdots & C_{MN} \end{bmatrix}$

求解得模态位移系数

$\left \{ \mathbf{\widetilde{q}} \right \}=\frac{P_{in}}{-\omega^2\left [\ \mathbf{M} \right ]+\left [\ \mathbf{Q} \right ]+j\omega\left [\ \mathbf{C} \right ]-\left [\ \mathbf{K} \right ]}\left \{ \mathbf{H} \right \}$

透射系数

$t_{p}=\frac{P_t}{P_{in}}=\frac{2\rho_0c_0\omega}{L_xL_y}\left \{ \mathbf{H} \right \}^T\frac{P_{in}}{-\omega^2\left [\ \mathbf{M} \right ]+\left [\ \mathbf{Q} \right ]+j\omega\left [\ \mathbf{C} \right ]-\left [\ \mathbf{K} \right ]}\left \{ \mathbf{H} \right \}$

隔声量

$STL=20log_{10}(\frac{1}{t_p})$

2.Matlab实现

2.1构建M、C矩阵

M、C矩阵形式类似，故采用同一个子函数计算

function [Mmn,Cmn] = GenerateMC(M,N,Lx,Ly,rhos,rho0,c0)
Mmn = [];
Cmn = [];
TempM = 0;    
TempC = 0;
for m = 1:M
    for n = 1:N
        for r = 1:M
            for s = 1:N
                 fun = @(x,y) sin(m.*pi.*x./Lx).*sin(n.*pi.*y./Ly).*sin(r.*pi.*x./Lx).*sin(s.*pi.*y./Ly);
                 Temp = integral2(fun,0,Lx,0,Ly);
                 TempM = TempM+rhos*Temp;
                 TempC = TempC+2*rho0*c0*Temp;
            end
        end
        Mmn = blkdiag(Mmn,TempM);
        Cmn = blkdiag(Cmn,TempC);
        TempM = 0;    
        TempC = 0;
    end
end
end

2.2构建K矩阵

function Kmn = GenerateK(M,N,Lx,Ly,T)
Kmn = [];
TempK = 0;
for m = 1:M
    for n = 1:N
        for r = 1:M
            for s = 1:N
                 fun = @(x,y) T*((r*pi/Lx)^2+(s*pi/Ly)^2)*sin(m.*pi.*x./Lx).*sin(n.*pi.*y./Ly).*sin(r.*pi.*x./Lx).*sin(s.*pi.*y./Ly);
                 Temp = integral2(fun,0,Lx,0,Ly);
                 TempK = TempK+Temp;
            end
        end
        Kmn = blkdiag(Kmn,TempK);
        TempK = 0;
    end
end
end

2.3构建H矩阵

function Hmn = GenerateH(M,N,Lx,Ly)
Hmn = zeros(M*N,1);
k = 1;
for m = 1:M
    for n = 1:N
        fun = @(x,y) sin(m.*pi.*x./Lx).*sin(n.*pi.*y./Ly);
        Hmn(k,1) = integral2(fun,0,Lx,0,Ly);
        k = k+1;
    end
end
end

2.4构建Q矩阵

function Qmn = GenerateQ(M,N,Lx,Ly,rhomass,x0,y0,lx,ly)
k = 1;
Qmn = zeros(M*N,M*N);
for m = 1:M
    for n = 1:N
        for r = 1:M
            for s = 1:N
                fun = @(x,y) sin(m.*pi.*x./Lx).*sin(n.*pi.*y./Ly).*sin(r.*pi.*x./Lx).*sin(s.*pi.*y./Ly);
                if mod(k,M*N)==0
                    Qmn(floor(k/(M*N)),9) = rhomass*integral2(fun,x0,x0+lx,y0,y0+ly);
                else
                    Qmn(floor(k/(M*N))+1,mod(k,M*N)) = rhomass*integral2(fun,x0,x0+lx,y0,y0+ly);
                end
                k = k+1;
            end
        end
    end
end
end

2.5主子函数

function MembraneMetamaterialPlate()
STL = zeros(10000,1);
F = zeros(10000,1);
%% 参数定义
j = sqrt(-1);
rho0 = 1.21;%空气密度
c0 = 343;%空气声速
rhos = 0.0912;%薄膜面密度
T = 487.4;%薄膜张力
Lx = 0.0274;%薄膜x方向尺寸
Ly = 0.0274;%薄膜y方向尺寸
m = 0.00032;%质量块质量，单位是kg
r = 0.00193;%质量块半径
lx = sqrt(pi*r^2);%质量块x方向尺寸
ly = sqrt(pi*r^2);%质量块y方向尺寸
x0 = (Lx-lx)/2;%质量块左下角位置x方向
y0 = (Ly-ly)/2;%质量块左下角位置y方向
rhomass = m/(lx*ly);%质量块面密度
M = 3;%x方向模态数
N = 3;%y方向模态数
%% 生成M、I、C、K、H矩阵
[Mmn,Cmn] = GenerateMC(M,N,Lx,Ly,rhos,rho0,c0);
Kmn = GenerateK(M,N,Lx,Ly,T);
Hmn = GenerateH(M,N,Lx,Ly);
Qmn = GenerateQ(M,N,Lx,Ly,rhomass,x0,y0,lx,ly);
for f = 1:10000
    w = 2*pi*f;
    tp = abs((2*rho0*c0*w/(Lx*Ly))*Hmn'*(((j*w^2*(Mmn+Qmn))+w*Cmn-j*Kmn)\(Hmn)));
    STL(f) = 20*log10(1/tp);
    F(f) = f;
end
semilogx(F,STL);
xlim([100,7000]);
end