clc;clear;close
a=0.00005; %波速
dx=0.002; dy=0.002; %空间步长c
dt=0.01; %时间步长
Lx=0.7; Ly=1; %求解域
Ti=150; %传播时间
U0=0.5; %声强
x=0:dx:Lx;y=0:dy:Ly;t=0:dt:Ti;
[X,Y]=meshgrid(x,y);
nx=length(x); ny=length(y);
U=zeros(length(x),length(y),2); %存放结果
i=1;j=1;k=2; %初始化变量
U(:,1,1)=U0;U(:,1,2)=U0;
beta=0.005;K=1/(1+beta*dt);%衰减系数设置
%找到左边的弧面==============================
Ucx=Lx/2; Ucy=Ly/2;
M=(((X-Ucx).^2+(Y-Ucy).^2)<(Ly/2-dx)^2);
M=M+M([end,1:end-1],:)
[Xbw,Ybw]=find((Y<Ly/2).*(M==1));
%主循环==============================
for it=dt:dt:Ti
i=1;k=k+1; %下标变量
for ix=dx:dx:Lx-dx
i=i+1; j=1; %下标变量
for mk=1:length(Ybw); U(mk,1:Xbw(mk),k)=0;end %左边界
U(i,ny,k)=U(i,ny-1,k); %右边界
for iy=dy:dy:Ly-dy
j=j+1;
U(1,j,k)=0*U(3,j,k); U(nx,j,k)=0*U(nx-2,j,k); %上下边界
U(i,j,k)=K*(a*(U(i,j-1,k-1)-2*U(i,j,k-1)+U(i,j+1,k-1))/dx^2*dt^2+…
K*(a*(U(i-1,j,k-1)-2*U(i,j,k-1)+U(i+1,j,k-1))/dy^2*dt^2)+…
+(2+beta*dt)*U(i,j,k-1)-U(i,j,k-2));
end
for mk=1:length(Ybw) %波动源
U(Ybw(mk),Xbw(mk),k)=sin(it);
end
end
% 后处理==============================
if mod(k,50)==0;Ut=U(:,:,k);imagesc(Ut,[-5,5]);
colormap jet; axis equal; axis off;%后处理
drawnow;k
end
end