时域有限差分法仿真一维TE波在分裂场完全匹配层【含源码】

1、时域有限差分法仿真一维TE波在分裂场完全匹配层吸收边界条件下的传输时域有限差分法仿真一维TE波在分裂场完全匹配层吸收边界条件下的传输一、时域有限差分法 (FDTD, Finite-Difference Time-Domain)FDTD是1966年K.S.Yee发表在AP上的一篇论文建立起来的，后被称为Yee网格空间离散方式核心思想是把带时间变量的Maxwell旋度方程转化为差分形式，模拟出电子脉冲和理想导体作用的时域响应号称目前计算电磁学界最受关注，最时髦的算法，但还在发展完善之中国外已有多种基于FDTD算法的电磁场计算的软件:XFDTD等。二、差分算法 三、时域有限差分算法步骤(1)采用一定

2、的网格划分方式离散化场域；(2)对场内的偏微分方程及各种边界条件进行差分离散化处理，建立差分格式，得到差分方程组；(3)结合选定的代数方程组的解法，编制程序，求边值问题的数值解。四、吸收边界条件1、问题的提出在电磁场的辐射和散射问题中，边界总是开放的，电磁场占据无限大空间，而计算机内存是有限的，所以只能模拟有限空间。即：时域有限差分网格将在某处被截断。这要求在网格截断处不能引起波的明显反射，因而对向外传播的波而言，就像在无限大的空间传播一样，一种行之有效的方法是在截断处设置一种吸收边界条件。使传播到截断出的波被边界吸收而不产生反射。2、良匹配层基本原理1994年Beernger提出一种新颖的由

3、吸收媒质构成的良匹配层(PML)的概念,这种人工设计的良匹配层由有耗,导电,导磁媒质组成,可吸收任意入射角,任意频率,任意偏振态的入射电磁波,且透射波以原来波速传播,且有相同的特征阻抗,在垂直入射方向衰减。根据推导，只要这种媒质的结构参数满足一定条件,这种媒质就可以起到完全吸收入射波的作用。五、计算所用到的电磁场公式根据麦克斯韦方程组和时域差分法，得到以下公式：六、仿真结果七、源代码%*% 时域有限差分法仿真一维TE波在分裂场完全匹配层吸收边界条件下的传输%*% Eref Ein Erans% PEC PML 介质 PML PEC% +-+-+-+-+-+-+-% z=0 z=Lz% r=1

4、Nref Nin Ntrans r=Nz+1%*% 物理常量mu0 = 4e-7 * pi; % 磁导率c0 = 299792458; % 光速eps0 = 1/c02/mu0; % 电导率eta0 = sqrt(mu0/eps0); % 波导率 GHz = 109;mm = 10(-3);%*% 参数初始化%*Lz = 1; % 长度单位：米Tmax = 4*Lz/c0; % 时间最大值Dz = 1*mm; % 空间网格尺寸R = 1/sqrt(3); %每一单位时间计算的网格单元数freq = 5*GHz; % 激励源的中心频率fmax = 9*GHz; % 画图的最大频率fmin = 1

5、*GHz; Nz = round(Lz/Dz); % z轴上的单元格数量（round 四舍五入）n_pict = round(Nz/20); % 每n_pict步画一次图linew = 2; % 画图的线宽Dt = Dz*R/c0; % 时间步长Nt = round(Tmax/Dt); % 时间步长数量Nfft = max(2(round(log2(Nt)+1),2*1024); %做傅里叶变换点的数量z = linspace(0,Lz,Nz+1); % Z轴坐标t = Dt*0:Nt-1; % 时间lambda = c0/freq; % 中心波长的激励源ppw = lambda/Dz % 中

6、心频率处每一段波长的采样点omega = 2.0*pi*freq; % 角频率Nabc = 1*10; % PML边界左边和右边的层数Nin = 4; % 表面散射场Nref = 2; % 由近到远场表面Ntrans = Nz-2;Labc = Nabc*Dz; % PML边界长度%*% 材料参数%*Nmedia=3; % 材料数量media_par = 1.0 0.0; % 真空：电容率和电导率 3.4 1*10(-12); % 3.0 0; objects = 0; % 介质的数量obj_z = Dz/2+0.500 0.504; % Z1和Z2之间的介质1 0.624 0.628; 0.

7、24 0.26;obj_m = 2 2 1; % 介质的材料obj_c = rrb; % 介质的颜色epsr = ones(Nz-1,1);sig = zeros(Nz-1,1);obj_ind = ;obj_ind(:,1) = ceil(obj_z(:,1)/Dz);obj_ind(:,2) = floor(obj_z(:,2)/Dz);obj_frac(:,1) = obj_ind(:,1)-obj_z(:,1)/Dz;obj_frac(:,2) = -obj_ind(:,2)+obj_z(:,2)/Dz;for n=1:objects epsr(obj_ind(n,1):obj_ind

8、(n,2) = media_par(obj_m(n),1); sig(obj_ind(n,1):obj_ind(n,2) = media_par(obj_m(n),2);end%*% 分配场矢量%*Ex = zeros(Nz+1 ,1);Hy = zeros(Nz,1);Eref = zeros(Nt,1);Etrans = zeros(Nt,1);% 吸收边界条件E1abc = zeros(Nabc+1,1); % 左端E2abc = zeros(Nabc+1,1); % 右端H1abc = zeros(Nabc,1);H2abc = zeros(Nabc,1);zpmlH = linspa

9、ce(0,Labc,Nabc)/Labc;zpmlE = linspace(Dz/2,Labc-Dz/2,Nabc-1)/Labc;abc_pow = 3;sigma_max = -(abc_pow+1)*log(10(-4)/(2*eta0*Labc);sigmaH2 = sigma_max*zpmlH.abc_pow*mu0/eps0;sigmaE2 = sigma_max*zpmlE.abc_pow;sigmaE1 = sigmaE2(Nabc-1):-1:1);sigmaH1 = sigmaH2(Nabc:-1:1);figure(4)plot(sigmaH2)%*% 波激励%*Hde

10、lay = Dt/2-Dz/2/c0;% 从左侧加入射脉冲tau = 80.0e-12;delay = 2.1*tau;Ein = sin(omega*(t-delay).*exp(-(t-delay).2/tau2);Hin = 1/eta0*sin(omega*(t-delay+Hdelay).*exp(-(t-delay+Hdelay).2/tau2);% 谐波从左侧入射的时间Einmax = max(abs(Ein);%*% 时间步进%*for n = 1:Nt; % 从-Dt/2到Dt/2更新磁场 Hy = Hy - (Dt/mu0/Dz) * diff(Ex,1); % 主网格 H

11、y(Nin) = Hy(Nin) - (Dt/mu0/Dz) * Ein(n); H1abc = (H1abc.*(mu0/Dt-sigmaH1/2) - diff(E1abc)/Dz) ./ (mu0/Dt+sigmaH1/2); % 从ABC向左 H2abc = (H2abc.*(mu0/Dt-sigmaH2/2) - diff(E2abc)/Dz) ./ (mu0/Dt+sigmaH2/2); % 从ABC向右 % 从0到Dt更新E Ex(2:Nz) = (Ex(2:Nz).*(eps0*epsr/Dt-sig/2) - diff(Hy,1)/Dz)./(eps0*epsr/Dt+sig

12、/2); % 除去边界的主网格 Ex(Nin) = Ex(Nin) - Hin(n)/(Dz*eps0/Dt); Ex(1) = Ex(1) - (Dt /eps0) * (Hy(1)-H1abc(Nabc)/Dz; %左边ABC和主网格之间的边界 E1abc(Nabc+1) = Ex(1); % 补充一点，以简化方案 E1abc(2:Nabc) = (E1abc(2:Nabc).*(eps0/Dt-sigmaE1/2) - diff(H1abc,1)/Dz) ./ (eps0/Dt+sigmaE1/2); Ex(Nz+1) = Ex(Nz+1) - (Dt /eps0) * (H2abc(1

13、)-Hy(Nz)/Dz; E2abc(1) = Ex(Nz+1); E2abc(2:Nabc) = (E2abc(2:Nabc).*(eps0/Dt-sigmaE2/2) - diff(H2abc,1)/Dz) ./ (eps0/Dt+sigmaE2/2); % 在选定点采样电场 if mod(n,n_pict) = 0 figure(1); subplot(2,1,1); plot(z,Ex,LineWidth,1); title(time is ,num2str(n*Dt*109),ns) axis tight; ax=axis; axis(ax(1:2) -1.2*Einmax 1.2*

14、Einmax); ax=axis; ylabel(E-field); for p=1:objects line(obj_z(p,1) obj_z(p,1), ax(3:4),color,obj_c(p); line(obj_z(p,2) obj_z(p,2), ax(3:4),color,obj_c(p); line(obj_z(p,1) obj_z(p,2), ax(3) ax(3),color,obj_c(p),LineWidth,4); line(obj_z(p,1) obj_z(p,2), ax(4) ax(4),color,obj_c(p),LineWidth,4); end sub

15、plot(2,1,2); plot(E1abc(:) Ex E2abc(:),LineWidth,1) title(time is ,num2str(n*Dt*109),ns) ax = axis; line(Nabc Nabc, ax(3:4),color,r); line(Nin+Nabc Nabc, ax(3:4),color,b); line(Ntrans+2+Nabc Nabc, ax(3:4),color,b); line(Nz+2+Nabc Nz+2+Nabc, ax(3:4),color,r); axis tight; ylabel(E-field); %plot(Eabc(2

16、:Nabc,1:2); pause(0.01); end Eref(n) = Ex(Nref); Etrans(n) = Ex(Ntrans);end%*% 后处理%*Ehin = fft(Ein,Nfft)/Nt*2;Df = 1/(Dt*Nfft);freqN = round(fmax/Df);freqN1 = round(fmin/Df);f = linspace(Df*freqN1,Df*freqN,freqN-freqN1+1)/109;t = t*109; % 纳米秒（ns）figure(2);% 入射波clf;subplot(2,3,1);plot(t,Ein,r,LineWid

17、th,linew)axis tight;xlabel(time);ylabel(E-field);subplot(2,3,4);warning off;A=plotyy(f,abs(Ehin(freqN1:freqN),f,c0./f/Dz/109); warning on;set(A(2),yLim,0 40,yTick,0:5:40,xlim,fmin fmax/109,xGrid,on,yGrid,on,LineWidth,linew);set(A(1),xlim,fmin fmax/109,LineWidth,linew);set(get(A(2),Ylabel),String,poi

18、nts/wavelength);set(get(A(1),Ylabel),String,input);set(get(A(1),Xlabel),String,frequency);% 反射subplot(2,3,2);plot(t,Eref,LineWidth,linew); axis tight;xlabel(time ns);ylabel(Reflected field);subplot(2,3,5);Ehref = fft(Eref,Nfft)/Nt*2;plot(f,20*log10(abs(Ehref(freqN1:freqN)./Ehin(freqN1:freqN),LineWidth,linew); axis tight;ax=axis;axis(fmin/109 fmax/109 ax(3:4);axis tight;xlabel(frequency GHz);ylabel(Reflection in dB);% 传播subplot(2,3,3);plot(t,Etrans,LineWidth,linew); axis tight;xlabel(time ns);ylabel(Transmitted field);subplot(2,3,6);Ehtrans = fft(Etrans,Nfft)/Nt*2;plot(f,20