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

二、差分算法x h ∆= ()()f f x h f x ∆=+-0()()()()lim ()()()()=2x df f x f x f x h f x dx x x hf x f x h hf x h f x h h∆→∆∆+-=≈=∆∆--=+--前向差分 后向差分中心差分 222()1()()()2!df x d f x f x h f x h h dx dx +=+++222()1()()()2!df x d f x f x h f x h h dx dx -=-++333()2()()()23!df x d f x f x h f x h h h dx dx +--=++三、时域有限差分算法步骤(1)采用一定的网格划分方式离散化场域；(2)对场内的偏微分方程及各种边界条件进行差分离散化处理，建立差分格式，得到差分方程组；(3)结合选定的代数方程组的解法，编制程序，求边值问题的数值解。

四、吸收边界条件1、问题的提出在电磁场的辐射和散射问题中，边界总是开放的，电磁场占据无限大空间，而计算机内存是有限的，所以只能模拟有限空间。

即：时域有限差分网格将在某处被截断。

这要求在网格截断处不能引起波的明显反射，因而对向外传播的波而言，就像在无限大的空间传播一样，一种行之有效的方法是在截断处设置一种吸收边界条件。

使传播到截断出的波被边界吸收而不产生反射。

2、良匹配层基本原理1994年Beernger 提出一种新颖的由吸收媒质构成的良匹配层(PML)的概念,这种人工设计的良匹配层由有耗,导电,导磁媒质组成,可吸收任意入射角,任意频率,任意偏振态的入射电磁波,且透射波以原来波速传播,且有相同的特征阻抗,在垂直入射方向衰减。

根据推导，只要这种媒质的结构参数满足一定条件,这种媒质就可以起到完全吸收入射波的作用。

五、计算所用到的电磁场公式根据麦克斯韦方程组和时域差分法，得到以下公式：()()()()1/21/21/21/2nn n n y y x x H k H k E k E k t x +-+---=∆∆()()()()11/21/201/21/211n n n n y y x x H k H k E k E k txμ++++--+-=-⋅∆∆12=1/21/21/2()()[(1/2)(1/2)]n n n nx x y y rE k E k H k H k ε++=++--1/21/21(1/2)=(1/2)-[(1)()]2n nn n y y x x H k H k E k E k +++++-六、仿真结果timeE -f i e l di n p u tfrequencyp o i n t s /w a v e l e n g t h-4time [ns]R e f l e c t e d f i e l dfrequency [GHz]R e f l e c t i o n i n d Btime [ns]T r a n s m i t t e d f i e l d-3frequency [GHz]T r a n s m i s s i o n i n d B七、源代码%*********************************************************************** % 时域有限差分法仿真一维TE波在分裂场完全匹配层吸收边界条件下的传输%*********************************************************************** %%% Eref Ein Erans% PEC PML 介质PML PEC% +-----+---+------+---------------------+-----+------+-->% z=0 z=Lz% r=1 Nref Nin Ntrans r=Nz+1%%*********************************************************************** % 物理常量mu0 = 4e-7 * pi; % 磁导率c0 = 299792458; % 光速eps0 = 1/c0^2/mu0; % 电导率eta0 = sqrt(mu0/eps0); % 波导率GHz = 10^9;mm = 10^(-3);%*********************************************************************** % 参数初始化%*********************************************************************** Lz = 1; % 长度单位：米Tmax = 4*Lz/c0; % 时间最大值Dz = 1*mm; % 空间网格尺寸R = 1/sqrt(3); %每一单位时间计算的网格单元数freq = 5*GHz; % 激励源的中心频率fmax = 9*GHz; % 画图的最大频率fmin = 1*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 % 中心频率处每一段波长的采样点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之间的介质10.624 0.628;0.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:objectsepsr(obj_ind(n,1):obj_ind(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 = linspace(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)%***********************************************************************% 波激励%***********************************************************************Hdelay = Dt/2-Dz/2/c0;% 从左侧加入射脉冲tau = 80.0e-12;delay = 2.1*tau;Ein = sin(omega*(t-delay)).*exp(-((t-delay).^2/tau^2));Hin = 1/eta0*sin(omega*(t-delay+Hdelay)).*exp(-((t-delay+Hdelay).^2/tau^2));% 谐波从左侧入射的时间Einmax = max(abs(Ein));%***********************************************************************% 时间步进%***********************************************************************for n = 1:Nt;% 从-Dt/2到Dt/2更新磁场Hy = Hy - (Dt/mu0/Dz) * diff(Ex,1); % 主网格Hy(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更新EEx(2:Nz) = (Ex(2:Nz).*(eps0*epsr/Dt-sig/2) - diff(Hy,1)/Dz)./(eps0*epsr/Dt+sig/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)-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) == 0figure(1);subplot(2,1,1);plot(z,Ex,'LineWidth',1);title(['time is ',num2str(n*Dt*10^9),'ns'])axis tight;ax=axis;axis([ax(1:2) -1.2*Einmax 1.2*Einmax]);ax=axis;ylabel('E-field');for p=1:objectsline([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);endsubplot(2,1,2);plot([E1abc(:)' Ex' E2abc(:)'],'LineWidth',1)title(['time is ',num2str(n*Dt*10^9),'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:Nabc,1:2));pause(0.01);endEref(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)/10^9;t = t*10^9; % 纳米秒（ns）figure(2);% 入射波clf;subplot(2,3,1);plot(t,Ein,'r','LineWidth',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/10^9);warning on;set(A(2),'yLim',[0 40],'yTick',[0:5:40],'xlim',[fmin fmax]/10^9,'xGrid','on','yGrid','on','LineWidth',linew); set(A(1),'xlim',[fmin fmax]/10^9,'LineWidth',linew);set(get(A(2),'Ylabel'),'String','points/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/10^9 fmax/10^9 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*log10(abs(Ehtrans(freqN1:freqN)./Ehin(freqN1:freqN))),'LineWidth',linew); axis tight;ax=axis;axis([fmin/10^9 fmax/10^9 ax(3:4)]);xlabel('frequency [GHz]');ylabel('Transmission in dB');tr=interp1(f,abs(Ehtrans(freqN1:freqN)./Ehin(freqN1:freqN)),5) re=interp1(f,abs(Ehref(freqN1:freqN)./Ehin(freqN1:freqN)),5) tr^2+re^2。