clc;
clear;
%带有完全匹配层PML
dt=0.002; %时间采样间隔，单位ms
t=dt:dt:3; %时间总长，单位ms
f0=20; %声源频率，单位kHz
ft=(1-2*pi^2*f0^2.*(t-1/f0).^2).*exp(-pi^2*f0^2.*(t-1/f0).^2); %声源
dr=0.01;
dz=0.01; %空间网格步长，单位m
r=-2:dr:2;
z=-2:dz:4; %计算空间大小，单位m
r0=0;
z0=0; %声源位置
l=length(t);
m=length(r);
n=length(z);
vp=zeros(n,m); % 介质纵波速度
vs=zeros(n,m); % 介质横波速度
p=zeros(n,m); % 介质密度
L=zeros(n,m);
N=zeros(n,m); % 拉梅系数
m0=find(r==r0);
n0=find(z==z0); %声源在网格中的位置
%模型
for i=1:n
for j=1:m
if abs(r(j))<0.1 %井径，单位m
%井中参数
vp(i,j)=1.5; %声波速度，单位m/ms
vs(i,j)=0;
p(i,j)=1000; %密度，单位kg/m3
else vp(i,j)=3.5;
vs(i,j)=2.0;
p(i,j)=2500;
end
N(i,j)=p(i,j)*vs(i,j)^2;
L(i,j)=p(i,j)*(vp(i,j)^2-2*vs(i,j)^2);
% if abs(r(j))==0.1
% N(i+1,j)=2*(1/N(i,j)+1/N(i+2,j))^(-1);
% p(i+1,j)=0.5*(p(i,j)+p(i+2,j)); %井壁处流固界面条件
% end
end
end
% for i=1:n
% for j=425:435
% if i+100>=sqrt(3)/2*j&&i+100<=sqrt(3)/2*j+2
% vp(i,j)=2;
% end
% end
% end
%************************稳定性条件**************************
maxvp=max(max(vp));
maxvs=max(max(vs));
if sqrt(2)*maxvp*dt/dr>=1
error('参数设置有问题，差分不稳定');
end
P=zeros(n,m);
Q=P;S=P;U=P;W=P;
P1=P;Q1=Q;S1=S;U1=U;W1=W; %应力和速度的物理量
U1x=P;U1z=P;Ux=P;Uz=P;
W1x=P;W1z=P;Wx=P;Wz=P;
P1x=P;P1z=P;Px=P;Pz=P;
Q1x=P;Q1z=P;Qx=P;Qz=P;
S1x=P;S1z=P;Sx=P;Sz=P; %分裂式PML吸收条件
l0=round(t/dt)+1;
pml=20; %吸收层网格
DP=pml*dr; %吸收层厚度
%衰减函数
dxj=zeros(n,m);
dxj2=dxj;
dzi=dxj;
dzi2=dxj;
xoleft=DP;
xoright=(m-pml)*dr;
RC=0.0001;
for i=1:n
for j=1:m
x=j*dr;
if j>=1 && j<pml
do= 3*vp(i,j)/(2*DP)*log10(1/RC);
dxj(i,j)=do*((xoleft-x)/DP)^4;
% dxj2(i,j)=do*((xoleft-x-0.5*dr)/DP)^4;
elseif j>m-pml && j<=m
do= 3*vp(i,j)/(2*DP)*log10(1/RC);
dxj(i,j)=do*((x-xoright)/DP)^4;
% dxj2(i,j)=do*((x+0.5*dr-xoright)/DP)^4;
else
dxj(i,j)=0;dxj2(i,j)=0;
end
z=i*dz;
if i>=1 && i<pml
do= 3*vp(i,j)/(2*DP)*log10(1/RC);
dzi(i,j)=do*((xoleft-z)/DP)^4;
% dzi2(i,j)=do*((xoleft-z-0.5*dz)/DP)^4;
elseif i>n-pml && i<=n
do= 3*vp(i,j)/(2*DP)*log10(1/RC);
dzi(i,j)=do*((z-xoright)/DP)^4;
% dzi2(i,j)=do*((z+0.5*dz-xoright)/DP)^4;
else
dzi(i,j)=0;dzi2(i,j)=0;
end
end
end
for k=1:800
k
if k<=l
P1x(n0,m0)=P1x(n0,m0)+ft(k);
% P1z(n0,m0)=ft(k);
Q1x(n0,m0)=Q1x(n0,m0)+ft(k);
end
for i=2:n-1
for j=2:m-1
px=0.5*(p(i,j)+p(i,j+1));
pz=0.5*(p(i,j)+p(i+1,j));
if j~=m0
y1=(dt/px)*((P1x(i,j+1)-P1x(i,j))/dr+(P1z(i,j+1)-P1z(i,j))/dr+(P1x(i,j+1)+P1x(i,j)+P1z(i,j+1)+P1z(i,j))/(2*(abs(j-m0)+1)*dr));
y2=(dt/px)*((S1x(i,j)-S1x(i-1,j))/dz+(S1z(i,j)-S1z(i-1,j))/dz);
y3=(dt/pz)*((S1x(i,j)-S1x(i,j-1))/dr+(S1z(i,j)-S1z(i,j-1))/dr+(S1x(i,j)+S1x(i,j-1)+S1z(i,j)+S1z(i,j-1))/(2*(abs(j-m0)+1)*dr));
y4=(dt/pz)*((Q1x(i+1,j)-Q1x(i,j))/dz+(Q1z(i+1,j)-Q1z(i,j))/dz);
else
y1=(dt/px)*((P1x(i,j+1)-P1x(i,j))/dr+(P1z(i,j+1)-P1z(i,j))/dr+(P1x(i,j+1)+P1x(i,j)+P1z(i,j+1)+P1z(i,j))/(2*dr));
y2=(dt/px)*((S1x(i,j)-S1x(i-1,j))/dz+(S1z(i,j)-S1z(i-1,j))/dz);
y3=(dt/pz)*((S1x(i,j)-S1x(i,j-1))/dr+(S1z(i,j)-S1z(i,j-1))/dr+(S1x(i,j)+S1x(i,j-1)+S1z(i,j)+S1z(i,j-1))/(2*dr));
y4=(dt/pz)*((Q1x(i+1,j)-Q1x(i,j))/dz+(Q1z(i+1,j)-Q1z(i,j))/dz);
end
Ux(i,j)=(1/(1+0.5*dt*dxj(i,j)))*((1-0.5*dt*dxj(i,j))*U1x(i,j)+y1);
Uz(i,j)=(1/(1+0.5*dt*dzi(i,j)))*((1-0.5*dt*dzi(i,j))*U1z(i,j)+y2);
Wx(i,j)=(1/(1+0.5*dt*dxj(i,j)))*((1-0.5*dt*dxj(i,j))*W1x(i,j)+y3);
Wz(i,j)=(1/(1+0.5*dt*dzi(i,j)))*((1-0.5*dt*dzi(i,j))*W1z(i,j)+y4);
end
end
U=Ux+Uz;W=Wx+Wz;
U1x=Ux;W1x=Wx;
U1z=Uz;W1z=Wz;
for i=2:n-1
for j=2:m-1
Nxz=1/(0.5*(1/N(i,j)+1/N(i,j+1)));
if j~=m0
y5=dt*(L(i,j)+2*N(i,j))*((Ux(i,j)-Ux(i,j-1))/dr+(Uz(i,j)-Uz(i,j-1))/dr+(Ux(i,j)+Ux(i,j-1)+Uz(i,j)+Uz(i,j-1))/(2*(abs(j-m0)+1)*dr));
y6=dt*L(i,j)*((Wx(i,j)-Wx(i-1,j))/dz+(Wz(i,j)-Wz(i-1,j))/dz);
y7=dt*L(i,j)*((Ux(i,j)-Ux(i,j-1))/dr+(Uz(i,j)-Uz(i,j-1))/dr+(Ux(i,j)+Ux(i,j-1)+Uz(i,j)+Uz(i,j-1))/(2*(abs(j-m0)+1)*dr));
y8=dt*(L(i,j)+2*N(i,j))*((Wx(i,j)-Wx(i-1,j))/dz+(Wz(i,j)-Wz(i-1,j))/dz);
y9=dt*Nxz*((Wx(i,j+1)-Wx(i,j))/dr+(Wz(i,j+1)-Wz(i,j))/dr);
y0=dt*Nxz*((Ux(i+1,j)-Ux(i,j))/dz+(Uz(i+1,j)-Uz(i,j))/dz);
else
y5=dt*(L(i,j)+2*N(i,j))*((Ux(i,j)-Ux(i,j-1))/dr+(Uz(i,j)-Uz(i,j-1))/dr+(Ux(i,j)+Ux(i,j-1)+Uz(i,j)+Uz(i,j-1))/(2*dr));
y6=dt*L(i,j)*((Wx(i,j)-Wx(i-1,j))/dz+(Wz(i,j)-Wz(i-1,j))/dz);
y7=dt*L(i,j)*((Ux(i,j)-Ux(i,j-1))/dr+(Uz(i,j)-Uz(i,j-1))/dr+(Ux(i,j)+Ux(i,j-1)+Uz(i,j)+Uz(i,j-1))/(2*dr));
y8=dt*(L(i,j)+2*N(i,j))*((Wx(i,j)-Wx(i-1,j))/dz+(Wz(i,j)-Wz(i-1,j))/dz);
y9=dt*Nxz*((Wx(i,j+1)-Wx(i,j))/dr+(Wz(i,j+1)-Wz(i,j))/dr);
y0=dt*Nxz*((Ux(i+1,j)-Ux(i,j))/dz+(Uz(i+1,j)-Uz(i,j))/dz);
end
Px(i,j)=(1/(1+0.5*dt*dxj(i,j)))*((1-0.5*dt*dxj(i,j))*P1x(i,j)+y5);
Pz(i,j)=(1/(1+0.5*dt*dzi(i,j)))*((1-0.5*dt*dzi(i,j))*P1z(i,j)+y6);
Qx(i,j)=(1/(1+0.5*dt*dxj(i,j)))*((1-0.5*dt*dxj(i,j))*Q1x(i,j)+y7);
Qz(i,j)=(1/(1+0.5*dt*dzi(i,j)))*((1-0.5*dt*dzi(i,j))*Q1z(i,j)+y8);
Sx(i,j)=(1/(1+0.5*dt*dxj(i,j)))*((1-0.5*dt*dxj(i,j))*S1x(i,j)+y9);
Sz(i,j)=(1/(1+0.5*dt*dzi(i,j)))*((1-0.5*dt*dzi(i,j))*S1z(i,j)+y0);
end
end
P=Px+Pz;Q=Qx+Qz;S=Sx+Sz;
P1x=Px;Q1x=Qx;S1x=Sx;
P1z=Pz;Q1z=Qz;S1z=Sz;
for ii=1:m0-1
P(:,m0-ii)=P(:,m0+ii);
end
for ii=1:m0-1
S(:,m0-ii)=S(:,m0+ii);
end
for ii=1:m0-1
Q(:,m0-ii)=Q(:,m0+ii);
end
for ii=1:m0-1
W(:,m0-ii)=W(:,m0+ii);
end
for ii=1:m0-1
U(:,m0-ii)=-U(:,m0-1+ii);
end
PP0(k)=P(n0+200,m0);
% PP1(k)=P(n0+440,m0);
% PP2(k)=P(n0+480,m0);
% PP3(k)=P(n0+520,m0);
% PP4(k)=P(n0+560,m0);
% PP5(k)=P(n0+600,m0);
end
