固支边界有限大平板隔声理论及Matlab实现

机械工人杨师傅

已于 2022-11-18 20:52:55 修改

阅读量975

点赞数 1

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

于 2022-11-17 20:27:13 首次发布

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

版权

隔声同时被 2 个专栏收录

7 篇文章 13 订阅

订阅专栏

论文复现

7 篇文章 4 订阅

订阅专栏

1.理论基础(模态叠加法)

1.1基础参数定义

首先定义板的基本参数： $L_x \ L_y$ 为板的长度和宽度； $\rho_s$ 为板密度； $h$ 为板厚度； $D$ 为板的弯曲刚度； $E$ 为板杨氏模量； $v$ 为板泊松比； $\eta$ 为板损耗因子；

$D=\frac{E(1+j\eta)h^3}{12(1-v^2)}$

其次定义声学参数： $\rho_0$ 为空气密度； $c_0$ 为空气声速； $P_{in}$ 为入射声波幅值； $\phi$ 为入射俯仰角； $\theta$ 为入射方位角； $\omega$ 为角频率； $k_0=\frac{\omega}{c_0}$ 为波数； $k_x \ k_y \ k_z$ 为入射声波在坐标轴三个方向的波数分量；

$k_x=k_0sin\phi cos\theta$

$k_y=k_0sin\phi sin\theta$

$k_z=k_0cos\phi$

1.2动力学方程

在声场作用下板的动力学响应方程为

$D\nabla^4w+\rho_sh\frac{\partial^2w}{\partial t^2}=jw\rho_0(\Phi_{in}-\Phi_{t})$

梯度算子:

$\nabla^4=(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2})^2$

板的横向弯曲位移：

$w(x,y,t)=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}\Psi_{mn}(x,y)p_{mn}(t)$

入射声场、辐射声场速度势函数：

$\Phi_{in}=\Phi_{inc}+\Phi_{ref}=\sum_{m}^{\infty}\sum_{n}^{\infty}I_{mn}\Psi_{mn}(x,y)e^{-j(k_z-\omega t)}+\sum_{m}^{\infty}\sum_{n}^{\infty}\beta_{mn}\Psi_{mn}(x,y)e^{-j(k_z-\omega t)}$

$\Phi_{t}=\sum_{m}^{\infty}\sum_{n}^{\infty}\xi_{mn}\Psi_{mn}(x,y)e^{-j(k_z-\omega t)}$

固支边界模态型函数：

$\Psi_{mn}(x,y)=(1-cos(\frac{2\pi mx}{L_x}))(1-cos(\frac{2\pi ny}{L_y}))$

简支边界模态型函数：

$\Psi_{mn}(x,y)=sin(\frac{\pi mx}{L_x})sin(\frac{\pi ny}{L_y})$

以 $\alpha_{mn}$ 为模态位移系数的模态位移：

$p_{mn}=\alpha_{mn}e^{j\omega t}$

空气-板耦合界面上速度连续，入射声波 $I_{mn}$ 、反射声波 $\beta_{mn}$ 、透射声波 $\xi_{mn}$ 的速度势幅值系数关系：

$\beta_{mn}=I_{mn}-\frac{\omega \alpha_{mn}}{k_z}$

$\xi_{mn}=\frac{\omega \alpha_{mn}}{k_z}$

$I_{mn}=\int_{0}^{L_x}\int_{0}^{L_y}P_{in}e^{-j(xk_x+yk_y)}\Psi_{mn}(x,y)dxdy$

1.3模态叠加

运用模态叠加法，将动力学方程等式做模态型函数的加权积分：两边同时乘以模态型函数并对xy求积分

$\int_{0}^{L_x}\int_{0}^{L_y}(D\nabla^4w+\rho_sh\frac{\partial^2w}{\partial t^2}-jw\rho_0(\Phi_{in}-\Phi_{t}))\Psi_{rs}(x,y)dxdy$

利用1.2节中的公式进行展开

$4D\pi^4L_xL_y \left \{ \left[ 3(\frac{r}{Lx})^4+3(\frac{s}{Ly})^4+2(\frac{r}{Lx})^2(\frac{s}{Ly})^2 \right ]\alpha_{rs}+\sum_{n}^{\infty}2(\frac{r}{L_x})^4\alpha_{rn}+\sum_{m}^{\infty}2(\frac{s}{L_y})^4\alpha_{ms} \right \}+\frac{9L_xL_y}{4}Q_{rs}+\sum_{n}^{\infty}\frac{3L_xL_y}{2}Q_{rn}+\sum_{m}^{\infty}\frac{3L_xL_y}{2}Q_{ms}+\sum_{m}^{\infty}\sum_{n}^{\infty}L_xL_yQ_{mn}=2j\omega\rho_0I_{rs}$

其中

$Q_{mn}=-\rho_sh\omega^2\alpha_{mn}+\frac{i2\omega^2\rho_0}{k_z}\alpha_{mn}$

方程经过 $1\leqslant r\leqslant M ,1\leqslant s\leqslant N$ 模态截断，可以表示成关于模态位移系数 $\alpha_{rs}$ 的MN个方程组，写成矩阵形式

$[\mathbf{M}_{rs}]_{MN\times MN}[\mathbf{\alpha}_{rs}]_{MN\times 1}=[\mathbf{F}_{rs}]_{MN\times 1}$

其中

$[\mathbf{\alpha}_{rs}]_{MN\times 1}=[\alpha_{11}\cdots \alpha_{M1}\alpha_{12}\cdots \alpha_{M2}\cdots\alpha_{MN}]_{MN\times 1}^T$

$[\mathbf{F}_{rs}]_{MN\times 1}=2j\omega\rho_0[\mathbf{I}_{11}\cdots \mathbf{I}_{M1}\mathbf{I}_{12}\cdots\mathbf{I}_{M2}\cdots\mathbf{I}_{MN}]_{MN\times 1}^T$

$[\mathbf{M}_{rs}]_{MN\times MN}=4D\pi^4L_xL_y(\Delta_1^{*1}+\Delta_1^{*2}+\Delta_1^{*3})-(\rho_sh\omega^2+j\omega\rho_0\frac{2\omega}{k_z})(\Delta_2^{*1}+\Delta_2^{*2}+\Delta_2^{*3}+\Delta_2^{*4})$

$\lambda_{rs}^{*1}=3(\frac{r}{Lx})^4+3(\frac{s}{Ly})^4+2(\frac{r}{Lx})^2(\frac{s}{Ly})^2$

$\Delta_1^{*1} = \mathbf{diag}[\lambda_{11}^{*1}\cdots \lambda_{M1}^{*1}\lambda_{12}^{*1}\cdots \lambda_{M2}^{*1}\cdots\lambda_{MN}^{*1}]_{MN\times MN}^T$

$\lambda_{1s}^{*2}=\frac{2s^4}{L_y^4}\begin{bmatrix} 0 \; 1\;1 \;\cdots \;1\\ 1\;0 \;1 \; \cdots\;1\\ 1 \; 1\;0\;\cdots\;1\\ \vdots \; \vdots \;\vdots \;\vdots \; \\ 1\;1\;1\;\cdots 0\; \end{bmatrix}_{M\times N}$

$\Delta_1^{*2} = \mathbf{diag}[\lambda_{11}^{*2} \; \lambda_{12}^{*2}\; \lambda_{13}^{*2}\; \cdots \; \lambda_{1N}^{*2}]_{MN\times MN}^T$

$\lambda_{1}^{*3}=\frac{2}{L_x^4}\mathbf{diag}[1^4\; 2^4 \; \cdots M^4]_{M\times N}$

$\Delta_{1}^{*3}=\begin{bmatrix} 0 \; \lambda_{1}^{*3}\;\lambda_{1}^{*3} \;\cdots \;\lambda_{1}^{*3}\\ \lambda_{1}^{*3}\;0 \;\lambda_{1}^{*3} \; \cdots\;\lambda_{1}^{*3}\\ \lambda_{1}^{*3} \; \lambda_{1}^{*3}\;0\;\cdots\;\lambda_{1}^{*3}\\ \vdots \; \vdots \;\vdots \;\vdots \; \\ \lambda_{1}^{*3}\;\lambda_{1}^{*3}\;\lambda_{1}^{*3}\;\cdots 0\; \end{bmatrix}_{M\times N}$

$\Delta_{2}^{*1}=\frac{9L_xL_y}{4}\mathbf{diag}[1\; 1 \; \cdots 1]_{M\times N}$

$\lambda_{2}^{*2}=\frac{3L_xL_y}{2}\begin{bmatrix} 0 \; 1\;1 \;\cdots \;1\\ 1\;0 \;1 \; \cdots\;1\\ 1 \; 1\;0\;\cdots\;1\\ \vdots \; \vdots \;\vdots \;\vdots \; \\ 1\;1\;1\;\cdots 0\; \end{bmatrix}_{M\times N}$

$\Delta_2^{*2} = \mathbf{diag}[\lambda_{2}^{*2} \; \lambda_{2}^{*2}\; \lambda_{2}^{*2}\; \cdots \; \lambda_{2}^{*2}]_{MN\times MN}^T$

$\lambda_{2}^{*3}=\frac{3L_xL_y}{2}\mathbf{diag}[1\; 1 \; \cdots 1]_{M\times N}$

$\Delta_{2}^{*3}=\begin{bmatrix} 0 \; \lambda_{2}^{*3}\;\lambda_{2}^{*3} \;\cdots \;\lambda_{2}^{*3}\\ \lambda_{2}^{*3}\;0 \;\lambda_{2}^{*3} \; \cdots\;\lambda_{2}^{*3}\\ \lambda_{2}^{*3} \; \lambda_{2}^{*3}\;0\;\cdots\;\lambda_{2}^{*3}\\ \vdots \; \vdots \;\vdots \;\vdots \; \\ \lambda_{2}^{*3}\;\lambda_{2}^{*3}\;\lambda_{2}^{*3}\;\cdots 0\; \end{bmatrix}_{M\times N}$

$\lambda_{2}^{*4}=L_xL_y\begin{bmatrix} 0 \; 1\;1 \;\cdots \;1\\ 1\;0 \;1 \; \cdots\;1\\ 1 \; 1\;0\;\cdots\;1\\ \vdots \; \vdots \;\vdots \;\vdots \; \\ 1\;1\;1\;\cdots 0\; \end{bmatrix}_{M\times N}$

$\Delta_{2}^{*4}=\begin{bmatrix} 0 \; \lambda_{2}^{*4}\;\lambda_{2}^{*4} \;\cdots \;\lambda_{2}^{*4}\\ \lambda_{2}^{*4}\;0 \;\lambda_{2}^{*4} \; \cdots\;\lambda_{2}^{*4}\\ \lambda_{2}^{*4} \; \lambda_{2}^{*4}\;0\;\cdots\;\lambda_{2}^{*4}\\ \vdots \; \vdots \;\vdots \;\vdots \; \\ \lambda_{2}^{*4}\;\lambda_{2}^{*4}\;\lambda_{2}^{*4}\;\cdots 0\; \end{bmatrix}_{M\times N}$

$[\mathbf{F}_{rs}]_{MN\times 1}$ 运算过程可以简写为

$[\mathbf{F}_{rs}]_{MN\times 1}=2j\omega\rho_0P_{in}[f_{11}\cdots f_{M1}f_{12}\cdots f_{M2}\cdots f_{MN}]_{MN\times 1}^T$

$f_{mn}(k_x,k_y)=\left\{\begin{matrix} L_xL_y \; ; \; k_x=0\&k_y=0\\ \frac{4jn^2\pi^2L_x(1-e^{-jk_yL_y})}{k_y(k_y^2L_y^2-4n^2\pi^2)} \; ; \; k_x=0\&k_y\neq0\\ \frac{4jm^2\pi^2L_y(1-e^{-jk_xL_x})}{k_x(k_x^2L_x^2-4m^2\pi^2)} \; ; \; k_x\neq0\&k_y=0\\ \frac{-16m^2n^2\pi^4(1-e^{-jk_xL_x})(1-e^{-jk_yL_y})}{k_xk_y(k_x^2L_x^2-4m^2\pi^2)(k_y^2L_y^2-4n^2\pi^2)} \; ; \; k_x\neq0\&k_y\neq0\\ \end{matrix}\right.$

通过矩阵运算求得模态位移系数矩阵

$[\mathbf{\alpha}_{rs}]_{MN\times 1}=[\mathbf{M}_{rs}]_{MN\times MN}^{-1}[\mathbf{F}_{rs}]_{MN\times 1}$

1.4求解隔声量

根据入射声压 $p_{in}$ ，共轭入射声速 $v_{in}^*$ ，入射声能量可以写为

$\Pi_{in} = 1/2Re\int\int_Ap_{in}v_{in}^*dA$

根据

$p_{in} =j\omega\rho_0\Phi_{in}= j\omega\rho_0 \left[2Ie^{-j(k_xx+k_yy)}-\sum_{m,n}\frac{\omega}{k_z}\alpha_{1,mn}\times (1-cos(\frac{2\pi mx}{L_x}))(1-cos(\frac{2\pi ny}{L_y})) \right]$

入射声能量可以改写为

$\Pi_{in} =\frac{\rho_0\omega^2}{2c_0}\int\int_A\left | \Phi_{in} \right |dA=\frac{\rho_0\omega^2}{2c_0}\int\int_A\left |2Ie^{-j(k_xx+k_yy)}-\sum_{m,n}\frac{\omega}{k_z}\alpha_{1,mn}\times (1-cos\frac{2\pi mx}{L_x})(1-cos\frac{2\pi ny}{L_y}) \right |^2dA=\frac{\rho_0\omega^2}{2c_0}\left | 4I^2\int\int_Ae^{-2j(k_xx+k_yy)}dA-4I\frac{\omega}{k_z}\sum_{m,n}\alpha_{1,mn}\int\int_A e^{-j(k_xx+k_yy)}(1-cos\frac{2\pi mx}{L_x})(1-cos\frac{2\pi ny}{L_y})dA+(\frac{\omega}{k_z})^2\sum_{mn,kl}\alpha_{1,mn}\alpha_{1,kl}\int\int_A(1-cos\frac{2\pi mx}{L_x})(1-cos\frac{2\pi ny}{L_y})(1-cos\frac{2\pi kx}{L_x})(1-cos\frac{2\pi ly}{L_y})dA \right |$

同理投射声能量可以写为

$\Pi_{t} =\frac{\rho_0\omega^2}{2c_0}\int\int_A\left | \Phi_{t} \right |dA=\frac{\rho_0\omega^2}{2c_0}\int\int_A\left |\sum_{m,n}\frac{\omega}{k_z}\alpha_{1,mn}\times (1-cos\frac{2\pi mx}{L_x})(1-cos\frac{2\pi ny}{L_y}) \right |^2dA=\frac{\rho_0\omega^2}{2c_0}\left | (\frac{\omega}{k_z})^2\sum_{mn,kl}\alpha_{1,mn}\alpha_{1,kl}\int\int_A(1-cos\frac{2\pi mx}{L_x})(1-cos\frac{2\pi ny}{L_y})(1-cos\frac{2\pi kx}{L_x})(1-cos\frac{2\pi ly}{L_y})dA \right |$

2.Matlab实现

2.1构建 $\Delta$ 矩阵

function [delta11,delta12,delta13,delta21,delta22,delta23,delta24] = GetDelta(Lx,Ly,M,N)
%% delta11,delta12求解    
    delta11 = [];
    delta12 = [];
    diagone = ones(M,N)-diag(ones(M,1),0);
    for n = 1:1:N
        for m = 1:1:M
            lamdars1 = 3*((m/Lx)^4)+3*((n/Ly)^4)+2*((m/Lx)^2)*((n/Ly)^2);%lamdars1计算出来是数值
            delta11 = [delta11,lamdars1];
        end
        lamdars2 = (2*(n^4)/(Ly^4)).*diagone;%lamdars2计算出来是一个M*N的矩阵
        delta12 = blkdiag(delta12,lamdars2); 
    end
    delta11 = diag(delta11,0);    
%% delta13求解
    lamda13 = (zeros(M,N)+diag(ones(M,1),0));
    for i = 1:1:M
        lamda13(i,i) = (2/(Lx^4))*(i^4);
    end
    delta13 = GenerateDiag(lamda13,M,N,1);
%% delta21求解
    delta21 = (9*Lx*Ly/4)*diag(ones(M*N,1),0);
%% delta22求解
    lamda22 = (3*Lx*Ly/2).*diagone; 
    delta22 = GenerateDiag(lamda22,M,N,2);
%% delta23求解
    lamda23 = (3*Lx*Ly/2).*diag(ones(M,1),0);
    delta23 = GenerateDiag(lamda23,M,N,1);
%% delta24求解
    lamda24 = (Lx*Ly).*diagone;
    delta24 = GenerateDiag(lamda24,M,N,1);
end

构建过程中需要利用矩阵堆叠出对角矩阵，此功能写成子函数GenerateDiag(lamda,M,N,flag)

function delta = GenerateDiag(lamda,M,N,flag)
%% 构建以为同维度0矩阵对角元素、以目标矩阵为其他元素的对角阵
if flag==1
    tempx = [];
    tempy = [];
    tempz = [];
    %构建出全为lamda的矩阵
        for i = 1:1:M
            tempx = [tempx lamda];
        end
        for j = 1:1:N
            tempy = [tempy;tempx];
        end
    %构建出以lamda为对角的矩阵
        for i = 1:1:M
            tempz = blkdiag(tempz,lamda);
        end
    %做差得到结果
        delta = tempy-tempz;
%% 构建以相同目标矩阵为对角元素、以同维度0矩阵为其他元素的对角阵（其中输入lamda为目标矩阵）.0
elseif flag==2
    delta = [];    
    for i = 1:1:M
        delta = blkdiag(delta,lamda);
    end
%% flag标志位异常，抛出错误
else
    error('flag标志位异常！')
end
end

2.2构建 $[\mathbf{F}_{rs}]_{MN\times 1}$ 矩阵

function [F,I] = GetF(M,N,Pin,kx,ky,Lx,Ly,w,rho0,j)
    F = zeros(M*N,1);
    I = zeros(M*N,1);
    k = 1;
    for n = 1:1:N
        for m = 1:1:M
            I(k,1) = Pin*GetFx(kx,ky,Lx,Ly,j,m,n);
            F(k,1) = 2*j*w*rho0*I(k,1);
            k = k+1;
        end
    end
end

Fx分类讨论子函数

function fx = GetFx(kx,ky,a,b,j,m,n)
if kx==0&&ky==0
    fx = a*b;
elseif kx==0&&ky~=0
    fx = (4*j*n^2*pi^2*a*(1-exp(-j*ky*b)))/(ky*(ky^2*b^2-4*n^2*pi^2));
elseif kx~=0&&ky==0
    fx = (4*j*m^2*pi^2*b*(1-exp(-j*kx*a)))/(kx*(kx^2*a^2-4*m^2*pi^2));
elseif kx~=0&&ky~=0
    fx = -(16*m^2*n^2*pi^4*(1-exp(-j*ky*b))*(1-exp(-j*kx*a)))/(kx*ky*(ky^2*b^2-4*n^2*pi^2)*(kx^2*a^2-4*m^2*pi^2));
else
    error('kx、ky错误');
end

2.3求解隔声量

function TL = GetTL(alpha1,rho0,w,c0,I,kx,ky,kz,a,b,M,N,j)
%alpha由MN*1重排为M*N；为了（m,n）能直接取到对应元素
alpha1 = reshape(alpha1,[M,N]); 
K = M;
L = N;
Sum11 = 0;%算子11
Sum12 = 0;%算子12
Sum13 = 0;%算子13
for m = 1:1:M
    for n = 1:1:N
        if kx==0&&ky==0
             Sum11 = a*b;
             Sum12 = Sum12+alpha1(m,n)*a*b;
        elseif kx==0&&ky~=0
             Sum11 = j*a*(-1+exp(-2*j*ky*b))/(2*ky);
             Sum12 = Sum12+alpha1(m,n)*4*j*n^2*pi^2*a*(1-exp(-j*ky*b))/(ky*(ky^2*b^2-4*n^2*pi^2));
        elseif kx~=0&&ky==0
             Sum11 = j*b*(-1+exp(-2*j*kx*a))/(2*kx);
             Sum12 = Sum12+alpha1(m,n)*4*j*m^2*pi^2*b*(1-exp(-j*kx*a))/(kx*(kx^2*a^2-4*m^2*pi^2));
        elseif kx~=0&&ky~=0
             Sum11 = -(1-exp(-2*j*kx*a)-exp(-2*j*ky*b)+exp(-2*j*(kx*a+ky*b)))/(4*kx*ky);
             Sum12 = Sum12+alpha1(m,n)*16*(-1)*m^2*n^2*pi^4*(1-exp(-j*ky*b)*(1-exp(-j*kx*a)))/(kx*ky*(kx^2*a^2-4*m^2*pi^2)*(ky^2*b^2-4*n^2*pi^2));
         end   
        for k = 1:1:K
            for i = 1:1:L
                if m==k&&n==i
                    Sum13 = Sum13+alpha1(m,n)*alpha1(k,i)*(9*a*b/4);
                elseif m==k&&n~=i
                    Sum13 = Sum13+alpha1(m,n)*alpha1(k,i)*(3*a*b/2);
                elseif m~=k&&n==i
                    Sum13 = Sum13+alpha1(m,n)*alpha1(k,i)*(3*a*b/2);
                elseif m~=k&&n~=i
                    Sum13 = Sum13+alpha1(m,n)*alpha1(k,i)*(a*b);
                end   
            end
        end
    end
end
Win =(rho0*w^2/(2*c0))*abs(4*I^2*Sum11-(4*I*w/kz)*Sum12+(w/kz)^2*Sum13); 
Wr = (rho0*w^2/(2*c0))*abs((w/kz)^2*Sum13);
TL = 10*log10(Win/Wr);
end

2.4主子函数

function RigidPlate()
    j = sqrt(-1);
%% 预分配内存
    F=zeros(10000,1);
    TLbare = zeros(10000,1);
    bar= waitbar(0,'Simulation in process');%显示运算进度
%% 空气参数定义
    rho0 = 1.21;%空气密度
    c0 = 343;%空气声速
    phi = 0;%入射角
    theta = 0;%方位角
    Pin = 1;%入射波幅值
%% 板参数定义
    rho = 2700;%板密度
    E = 70e9;%板杨氏模量
    v = 0.3;%板泊松比
    h = 0.0005;%板厚度
    Lx = 0.06;%板x方向长度
    Ly = 0.06;%板y方向长度
%% 平板参数定义
    D = E*(1+j*0.005)*h^3/(12*(1-v^2));                              %板的弯曲刚度
    %模态截断
    M = 20;                                                          
    N = 20;
%% 求解delta算子
    [delta11,delta12,delta13,delta21,delta22,delta23,delta24] = GetDelta(Lx,Ly,M,N);
%% 隔声系数计算
    fmin = 1;%计算下限
    fmax = 10000;%计算上限
    df = 1;
    for f = fmin:df:fmax
        s=['Simulation in process:' num2str(ceil((f-fmin)*100/(fmax-fmin))) '%'];
        waitbar((f-fmin)/(fmax-fmin),bar,s);%进度显示
        w = 2*pi*f;
        k0 = w/c0;
        kx = k0*sin(phi)*cos(theta);
        ky = k0*sin(phi)*sin(theta);
        kz = k0*cos(phi);
        %求解Mrs、Frs
        Mrs = (4*D*(pi^4)*Lx*Ly).*(delta11+delta12+delta13)-(rho*h*(w^2)+j*w*rho0*2*w/kz).*(delta21+delta22+delta23+delta24);
        [Frs,I] = GetF(M,N,Pin,kx,ky,Lx,Ly,w,rho0,j);
        abare = Mrs\Frs;
        F(f) = f;
        TLbare(f) = GetTL(abare,rho0,w,c0,Pin,kx,ky,kz,Lx,Ly,M,N,j);
    end
    semilogx(F(fmin:df:fmax),TLbare(fmin:df:fmax),'k','linewidth',2);
end

仿真时调用主子函数即可。