利用matlab编程实现的声波的时域差分方程

一、 声波方程的数值计算

网格离散

声波基本方程如下
$$\begin{gathered} \frac{\partial p}{\partial t}=-\rho c^2\left(\frac{\partial v_x}{\partial x}+\frac{\partial v_y}{\partial y}+\frac{\partial v_z}{\partial z}\right) \\ \frac{\partial v_z}{\partial t}=-\frac{1}{\rho }\frac{\partial p}{\partial z} \end{gathered}$$
离散网格采用cartesian网格划分法，对于声压P、质点振速v，进行离散化得到：
$$p^n(i,j,k)=p[i\Delta x,j\Delta y,k\Delta z,n\Delta t]$$
速度离散
$$\begin{gathered} v_x^{n+1/2}(i-1/2,j,k)=v_x[(i-1/2)\Delta x,j\Delta y,k\Delta z,(n+1/2)\Delta t] \\ {v}_y^{n+1/2}(i,j-1/2,k)=v_x[i\Delta x,(j-1/2)\Delta y,k\Delta z,(n+1/2)\Delta t] \\ {v}_z^{n+1/2}(i,j,k-1/2)=v_x[i\Delta x,j\Delta y,(k-1/2)\Delta z,(n+1/2)\Delta t] \end{gathered}$$
将振动速度和声压进行离散，然后把质点振动速度和声压进行在时间上相差半步进行抽样，进而实现空间，时间的四维模型中实现交替抽样。
将离散后的方程带入声波方程如下：
$$\begin{gathered} p^{n+1}(i,j,k)=p^n(i,j,k)-\rho c^2(i,j,k)\cdot \Delta t\cdot \{\frac{1}{\Delta x}\left[v_x^{n+1/2}\left(i+\frac{1}{2},j,k\right)-v_x^{n+1/2}\left(i-\frac{1}{2},j,k\right)\right] \\ +\frac{1}{\Delta y}\left[v_y^{n+1/2}\left(i,j+\frac{1}{2},k\right)-v_y^{n+1/2}\left(i,j-\frac{1}{2},k\right)\right]+\frac{1}{\Delta z}\left[v_z^{n+1/2}\left(i,j+\frac{1}{2},k\right)-v_z^{n+1/2}\left(i,j-\frac{1}{2},k\right)\right]\}\} \\ {v}_x^{n+1/2}\left(i-\frac{1}{2},j,k\right)=v_x^{n-1/2}\left(i-\frac{1}{2},j,k\right)-\frac{\Delta t}{\rho \cdot \Delta x}\left[p^n(i,j,k)-p^n(i-1,j,k)\right] \\ {v}_y^{n+1/2}\left(i,j-\frac{1}{2},k\right)=v_y^{n-1/2}\left(i,j-\frac{1}{2},k\right)-\frac{\Delta t}{\rho \cdot \Delta y}\left[p^n(i,j,k)-p^n(i,j-1,k)\right] \\ {v}_z^{n+1/2}\left(i,j,k-\frac{1}{2}\right)=v_z^{n-1/2}\left(i,j,k-\frac{1}{2}\right)-\frac{\Delta t}{\rho \cdot \Delta z}\left[p^n(i,j,k)-p^n(i,j,k-1)\right] \end{gathered}$$

时域差分方程是一种显式差分方程，在求解过程中，需要考虑算法的稳定性问题。在求解过程中的不稳定性表现为，随着时间步数的增加，计算结果也会无限制的增加。之前研究者们，对时间步长的限制条件进行了讨论，得出数值解是否稳定取决于时间步长和空间步长的关系。在非均匀媒介构成的空间中，需要进行合理的选择，得到时域有限差分的时间空间关系如下：
$$\Delta t⩽\frac{1}{c_{\max }\sqrt{{\left(\frac{1}{\Delta x}\right)}^2+{\left(\frac{1}{\Delta y}\right)}^2+{\left(\frac{1}{\Delta z}\right)}^2}}$$
其中必须确定频率和步长的关系如下：
$$\Delta s<=\frac{{\lambda }_{\min }}{N}$$

设定非均匀区域为中间十字形状（简单设置非均匀区域），如果需要可以更加复杂非均匀区域的可以再设定一个密度或者振速矩阵，或者改成for循环，利用条件语句设定，如果需要联系我
代码如下：

clc;
clear;
global Nx Ny Nt dx dy dt
Nx=300; Ny=300; Nt=900;%设置采样点数，采样时间点数
dx=1e-4;    %x方向和Y方向的步长
dy=1e-4;     %空间步长
dt=1e-8;   %时间步长
f=1e6;      %激励频率
p=p_(1000,1000,1000,1000);   %均匀 pl（声速1，声速2，密度1，密度2）
p1=p_(1000,1000,1000,1500);  %非均匀
dp=p-p1;                     %二者差值 
for k=1:4:Nt
figure (1);
subplot(1,3,2);   
pcolor(p(:,:,k));
shading interp;
colormap('bone');
axis equal;axis([0,200,0,200]);title('圆形密度均匀30*30')

subplot(1,3,1);
pcolor(p1(:,:,k));
shading interp;
colormap('bone');
axis equal;axis([0,200,0,200]);title('非均匀')

subplot(1,3,3);
pcolor(dp(:,:,k));
shading interp;
colormap('bone');
axis equal;axis([0,200,0,200]);title('均匀')
end
function p=p_(c1,c2,q1,q2)
global Nx Ny Nt dx dy dt
Nx=300; Ny=300; Nt=900;%设置采样点数，采样时间点数
dx=1e-4;    %x方向和Y方向的步长
dy=1e-4;     %空间步长
dt=1e-8;   %时间步长
f=1e6;      %激励频率
%基础条件设定完成
ux=zeros(Nx,Ny,Nt);
uy=zeros(Nx,Ny,Nt);
p=zeros(Nx,Ny,Nt);
cx=zeros(Nx,Ny,Nt);
cy=zeros(Nx,Ny,Nt);
D=@(v)[0;(v(2:end-1)-v(1:end-2));0];
D1=@(v)[0,(v(2:end-1)-v(1:end-2)),0];
D2=@(v)[(v(2:end)-v(1:end-1)); 0];
D3=@(v)[(v(2:end)-v(1:end-1)),0];
for k=1:Nt-1
    p(100,100,k)=sin(2*pi*f*k*dt);  %中心声源点
    for j=1:Nx
        if (j-100)^2<=900
                    c=c1;   q=q1;              
          
            else 
                    c=c2; q=q2;                %介质
        end      
      ux(:,j,k+1)=ux(:,j,k)-dt/(dx*q)*D(p(:,j,k));
      uy(j,:,k+1)=uy(j,:,k)-dt/(dy*q)*D1(p(j,:,k));
      cx(:,j,k+1)= D2(ux(:,j,k+1));
      cy(j,:,k+1)= D3(uy(j,:,k+1));
    end
     for i=1:Nx
        for j=1:Nx
            
    p(i,j,k+1)=p(i,j,k)-(q*c^2*dt*(cx(i,j,k+1)/dx+cy(i,j,k+1)/dy));
        end
     end
end                  
end

图像如下：

当然为了直观也可以画二维数值图像
代码如下：

figure (2)
for k=1:4:Nt
plot((p(100,:,k)));axis([0,200,-1,1]);
drawnow;
end

这是声波方程的时域差分的数值计算，通过在循环中改变密度速度来改变传播介质。此代码没有施加吸收边界，因此瞬态声波传导到边界后会发生反射。邮箱dynyanhao@163.com