Bzy=zeros(N1,N1); Bzxx=zeros(nmax,2); %进入电磁场迭代计算 for tt=1:nmax for i=1:N1
if i>=np+1&&i<=N1-np di=0; elseif i<=np di=np-i+0.5; elseif i>=N1-np+1 di=np+i-N1-0.5;
end %di是采样点横向距PML内边界的距离 sigma_mx=sigma_max*(di/np)^M; for j=1:N1
if j>=np+1&&j<=N1-np dj=0; elseif j<=np dj=np-j+0.5; elseif j>=N1-np+1 dj=np+j-N1-0.5;
end %dj是采样点纵向距PML内边界的距离 sigma_my=sigma_max*(dj/np)^M;
if sigma_mx+sigma_my==0 %真空区 if j==100&&i==100
t=30; %可选择高斯脉冲 term=(tt-t);
% pulse=exp(-4*pi*term^2/20^2); pulse=sin(2*pi*tt/40); %可选正弦时谐源 c_miu=c*del_t/del_s;
Eterm1=c_miu*(Ex(i,j+1)-Ex(i,j));
Eterm2=c_miu*(Ey(i+1,j)-Ey(i,j));
Bz(i,j)=Bz(i,j)+Eterm1-Eterm2+pulse;%加入脉冲源 else
c_miu=c*del_t/del_s;
Eterm1=c_miu*(Ex(i,j+1)-Ex(i,j)); Eterm2=c_miu*(Ey(i+1,j)-Ey(i,j)); Bz(i,j)=Bz(i,j)+Eterm1-Eterm2; end
else %PML区
cpm=(1-2*c*sigma_mx*del_t)/(1+2*c*sigma_mx*del_t); cqm=c*del_t/(1+2*c*sigma_mx*del_t)/del_s; Bzx(i,j)=cpm*Bzx(i,j)-cqm*(Ey(i+1,j)-Ey(i,j));
cpm=(1-2*c*sigma_my*del_t)/(1+2*c*sigma_my*del_t); cqm=c*del_t/(1+2*c*sigma_my*del_t)/del_s; Bzy(i,j)=cpm*Bzy(i,j)+cqm*(Ex(i,j+1)-Ex(i,j)); Bz(i,j)=Bzx(i,j)+Bzy(i,j); end end end for i=2:N1
if i>=np+1&&i<=N1-np di=0; elseif i<=np di=np-i+1; elseif i>=N1-np+1 di=np+i-N1-1;
end %di是采样点横向距PML内边界的距离 sigma_ex=sigma_max*(di/np)^M; for j=1:N1
cam=(1-2*c*sigma_ex*del_t)/(1+2*c*sigma_ex*del_t);
cbm=c*del_t/(1+2*c*sigma_ex*del_t)/del_s; Ey(i,j)=cam*Ey(i,j)-cbm*(Bz(i,j)-Bz(i-1,j)); end end for i=1:N1 for j=2:N1
if j>=np+1&&j<=N1-np dj=0; elseif j<=np dj=np-j+1; elseif j>=N1-np+1 dj=np+j-N1-1;
end %dj是采样点纵向距PML内边界的距离 sigma_ey=sigma_max*(dj/np)^M;
cam=(1-2*c*sigma_ey*del_t)/(1+2*c*sigma_ey*del_t); cbm=c*del_t/(1+2*c*sigma_ey*del_t)/del_s; Ex(i,j)=cam*Ex(i,j)+cbm*(Bz(i,j)-Bz(i,j-1)); end end
Bzxx(tt,1)=Bz(90,50); %对靠近边界的中央磁场点采样 Bzxx(tt,2)=Bz(50,90);
3.2 出图及结果 3.2.1程序部分
figure(1); %可视化处理 clf;
mesh(Bz); %磁场的幅值 axis([0 N 0 N -0.5 0.5]); xlabel('i') ylabel('j')
drawnow; end figure(2); plot(Bzxx); figure(3);
title('磁场幅值分布图'); surface(Bz); shading interp; axis square; toc
3.2.2 所出的效果图
正弦时谐场源辐射效果图
高斯脉冲辐射效果图