📄 pml3d.m
字号:
clear
c=3e8;
muz=4*pi*1e-7;
epsz=1/(c*c*muz);
f=1e+9;
lambda=c/f;
omega=2*pi*f;
ds=lambda/12;
dt=ds/(2*c);
muz=4*pi*1e-7;
R0=1e-5;
%************************************************
ie=50; %x向空间步
je=50; %y向空间步
ke=50; %z向空间步
ib=ie+1;
jb=je+1;
kb=ke+1;
is=25; %激励源x向位置
js=25; %激励源y向位置
ks=25; %激励源z向位置
iebc=8; %左、右PML层层数
jebc=8; %上、下PML层层数
kebc=8; %前、后PML层层数
ibbc=iebc+1;
jbbc=jebc+1;
kbbc=kebc+1;
nmax=100; %时间步
dlta=3; %观测点距激励源步数
%************************************************
Exy=zeros(ie,jb,kb);
Exz=zeros(ie,jb,kb);
Eyx=zeros(ib,je,kb);
Eyz=zeros(ib,je,kb);
Ezx=zeros(ib,jb,ke);
Ezy=zeros(ib,jb,ke);
Hxy=zeros(ie,jb,kb);
Hxz=zeros(ie,jb,kb);
Hyx=zeros(ib,je,kb);
Hyz=zeros(ib,je,kb);
Hzx=zeros(ib,jb,ke);
Hzy=zeros(ib,jb,ke);
fid1=fopen('m1.txt','wt');
fid2=fopen('m2.txt','wt');
%************************************************
sigxmax=-log(R0)/(2*iebc*ds)*3*epsz*c; %x向最大电导率
sigymax=-log(R0)/(2*jebc*ds)*3*epsz*c; %y向最大电导率
sigzmax=-log(R0)/(2*kebc*ds)*3*epsz*c;
% ca=linspace(ds,ie*ds,ie);
% cb=linspace(ds,ie*ds,ie);
% cp=linspace(ds,ie*ds,ie);
% cq=linspace(ds,ie*ds,ie);
% rou=linspace(ds,ie*ds,ie);
for i=1:ie %x向维电导率分布
if i<=iebc
sigx(i)=sigxmax*((iebc+1-i)/iebc)^2; %左侧PML层
elseif i>=ie-iebc+1 %右侧PML层
sigx(i)=sigxmax*((i-ie+iebc)/iebc)^2;
else
sigx(i)=0; %中间自由空间
end
roux(i)=muz*sigx(i)/epsz; %磁导率
cax(i)=(1-(dt*sigx(i))/(2*epsz))/(1+(dt*sigx(i))/(2*epsz));
cbx(i)=(dt/epsz/ds)/(1+(dt*sigx(i))/(2*epsz));
cpx(i)=(1-(dt*roux(i))/(2*muz))/(1+(dt*roux(i))/(2*muz));
cqx(i)=(dt/muz/ds)/(1+(dt*roux(i))/(2*muz));
end
for j=1:je %y向维电导率分布
if j<=jebc
sigy(j)=sigymax*((jebc+1-j)/jebc)^2; %下侧PML层
elseif j>=je-jebc+1 %上侧PML层
sigy(j)=sigymax*((j-je+jebc)/jebc)^2;
else
sigy(j)=0; %中间自由空间
end
rouy(j)=muz*sigy(j)/epsz; %磁导率
cay(j)=(1-(dt*sigy(j))/(2*epsz))/(1+(dt*sigy(j))/(2*epsz));
cby(j)=(dt/epsz/ds)/(1+(dt*sigy(j))/(2*epsz));
cpy(j)=(1-(dt*rouy(j))/(2*muz))/(1+(dt*rouy(j))/(2*muz));
cqy(j)=(dt/muz/ds)/(1+(dt*rouy(j))/(2*muz));
end
for k=1:ke %z向维电导率分布
if k<=kebc
sigz(k)=sigzmax*((kebc+1-k)/kebc)^2; %前侧PML层
elseif k>=ke-kebc+1 %后侧PML层
sigz(k)=sigzmax*((k-ke+kebc)/kebc)^2;
else
sigz(k)=0; %中间自由空间
end
rouz(k)=muz*sigz(k)/epsz; %磁导率
caz(k)=(1-(dt*sigz(k))/(2*epsz))/(1+(dt*sigz(k))/(2*epsz));
cbz(k)=(dt/epsz/ds)/(1+(dt*sigz(k))/(2*epsz));
cpz(k)=(1-(dt*rouz(k))/(2*muz))/(1+(dt*rouz(k))/(2*muz));
cqz(k)=(dt/muz/ds)/(1+(dt*rouz(k))/(2*muz));
end
% Ex=Exy+Exz;
% Ey=Eyx+Eyz;
% Ez=Ezy+Ezx;
%************************************************
tic
figure(1);clf;
set(gcf,'Renderer','OpenGL');
for n=1:nmax
for i=2:ie
for j=1:je
for k=1:ke
if(i==is&&j==js&&k==ks)
Hxy(i,j,k)=sin(omega*n*dt); %源激励
else
Hxy(i,j,k)=cpy(j).*Hxy(i,j,k)+cqy(j).*(Ezx(i,j,k)+Ezy(i,j,k)-Ezx(i,j+1,k)-Ezy(i,j+1,k));
end
end
end
end
for i=1:ie %
for j=1:je
for k=1:ke
Hxz(i,j,k)=cpz(k).*Hxz(i,j,k)+cqz(k).*(Eyx(i,j,k+1)+Eyz(i,j,k+1)-Eyx(i,j,k)-Eyz(i,j,k));
end
end
end
% Hx=Hxy+Hxz;
for i=1:ie
for j=1:je
for k=1:ke
Hyz(i,j,k)=cpz(k).*Hyz(i,j,k)+cqz(k).*(Exy(i,j,k)+Exz(i,j,k)-Exy(i,j,k+1)-Exz(i,j,k+1));
end
end
end
for i=1:ie
for j=1:je
for k=1:ke
Hyx(i,j,k)=cpx(i).*Hyx(i,j,k)+cqx(i).*(Ezx(i+1,j,k)+Ezy(i+1,j,k)-Ezx(i,j,k)-Ezy(i,j,k));
end
end
end
% Hy=Hyx+Hyz;
for i=1:ie
for j=1:je
for k=1:ke
Hzx(i,j,k)=cpx(i).*Hzx(i,j,k)+cqx(i).*(Eyx(i,j,k)+Eyz(i,j,k)-Eyx(i+1,j,k)-Eyz(i+1,j,k));
end
end
end
for i=1:ie
for j=1:je
for k=1:ke
Hzy(i,j,k)=cpy(j).*Hzy(i,j,k)+cqy(j).*(Exy(i,j+1,k)+Exz(i,j+1,k)-Exy(i,j,k)-Exz(i,j,k));
end
end
end
% Hz=Hzx+Hzy;
for i=1:ie
for j=2:je
for k=1:ke
Exy(i,j,k)=cpy(j).*Exy(i,j,k)+cqy(j).*(Hzx(i,j,k)+Hzy(i,j,k)-Hzx(i,j-1,k)-Hzy(i,j-1,k));
end
end
end
for i=1:ie
for j=1:je
for k=2:ke
Exz(i,j,k)=cpz(k).*Exz(i,j,k)+cqz(k).*(Hyx(i,j,k-1)-Hyx(i,j,k)+Hyz(i,j,k-1)-Hyz(i,j,k));
end
end
end
for i=1:ie
for j=1:je
for k=2:ke
Eyz(i,j,k)=cpz(k).*Eyz(i,j,k)+cqz(k).*(Hxy(i,j,k)-Hxy(i,j,k-1)+Hxz(i,j,k)-Hxz(i,j,k-1));
end
end
end
for i=2:ie
for j=1:je
for k=1:ke
Eyx(i,j,k)=cpx(i).*Eyx(i,j,k)+cqx(i).*(Hzx(i-1,j,k)-Hzx(i,j,k)+Hzy(i-1,j,k)-Hzy(i,j,k));
end
end
end
for i=2:ie
for j=1:je
for k=1:ke
Ezx(i,j,k)=cpx(i).*Ezx(i,j,k)+cqx(i).*(Hyx(i,j,k)-Hyx(i-1,j,k)+Hyz(i,j,k)-Hyz(i-1,j,k));
end
end
end
for i=1:ie
for j=2:je
for k=1:ke
Ezy(i,j,k)=cpy(j).*Ezy(i,j,k)+cqy(j).*(Hxy(i,j-1,k)-Hxy(i,j,k)+Hxz(i,j-1,k)-Hxz(i,j,k));
end
end
end
%************************************************
m1(n)=Hyx(is-dlta,js,ks)+Hyz(is-dlta,js,ks);
m2(n)=Hyx(is+dlta,js,ks)+Hyz(is+dlta,js,ks);
end
toc
fprintf(fid1,'%f\n',m1);
fprintf(fid2,'%f\n',m2);
fclose(fid1);
fclose(fid2);
plot(m1)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -