计算电磁学之FDTD算法的MATLAB语言实现

1、South China Normal University课程设计实验报告课程名称： 计算电磁学 指导老师： 专业班级： 2014级 电路与系统姓 名： 学 号： FDTD算法的MATLAB语言实现摘要：时域有限差分(FDTD)算法是K.S.Yee于1966年提出的直接对麦克斯韦方程作差分处理,用来解决电磁脉冲在电磁介质中传播和反射问题的算法。其基本思想是:FDTD计算域空间节点采用Yee元胞的方法,同时电场和磁场节点空间与时间上都采用交错抽样；把整个计算域划分成包括散射体的总场区以及只有反射波的散射场区,这两个区域是以连接边界相连接,最外边是采用特殊的吸收边界,同时在这两个边界之间有个输出边

2、界,用于近、远场转换;在连接边界上采用连接边界条件加入入射波,从而使得入射波限制在总场区域；在吸收边界上采用吸收边界条件,尽量消除反射波在吸收边界上的非物理性反射波。本文主要结合FDTD算法边界条件特点,在特定的参数设置下，用MATLAB语言进行编程，在二维自由空间TEz网格中，实现脉冲平面波。关键词：FDTD；MATLAB；算法1 绪论1.1 课程设计背景与意义 20世纪60年代以来，随着计算机技术的发展，一些电磁场的数值计算方法逐步发展起来，并得到广泛应用，其中主要有：属于频域技术的有限元法（FEM）、矩量法(MM)和单矩法等；属于时域技术方面的时域有限差分法(FDTD)、传输线矩阵法(T

3、LM)和时域积分方程法等。其中FDTD是一种已经获得广泛应用并且有很大发展前景的时域数值计算方法。时域有限差分（FDTD）方法于1966年由K.S.Yee提出并迅速发展，且获得广泛应用。K.S.Yee用后来被称作Yee氏网格的空间离散方式，把含时间变量的Maxwell旋度方程转化为差分方程，并成功地模拟了电磁脉冲与理想导体作用的时域响应。但是由于当时理论的不成熟和计算机软硬件条件的限制，该方法并未得到相应的发展。20世纪80年代中期以后，随着上述两个条件限制的逐步解除，FDTD便凭借其特有的优势得以迅速发展。它能方便、精确地预测实际工程中的大量复杂电磁问题，应用范围几乎涉及所有电磁领域，成为电

4、磁工程界和理论界研究的一个热点。目前，FDTD日趋成熟，并成为分析大部分实际电磁问题的首选方法。1.2 时域有限差分法的发展与应用经过四十多年的发展，FDTD已发展成为一种成熟的数值计算方法。在发展过程中，几乎都是围绕几个重要问题展开的，即数值稳定性、计算精度、数值色散、激励源技术以及开域电磁问题的吸收边界条件等。数值稳定和计算精度对任何一种数值计算方法都是至关重要的。其中激励源的设计和引入也是FDTD的一个重要任务。目前，应用最广泛的激励源引入技术是总场/散射场体系。对于散射问题，通常在FDTD计算空间中引入连接边界，它将整个计算空间划分为内部的总场区和外部的散射场区，如图1.1。利用Huy

5、gens原理，可以在连接边界处引入入射场，使入射场的加入变得简单易行。图1.1总场/散射场区2 FDTD法基本原理 2.1 Maxwell方程和Yee氏算法根据电磁场基本方程组的微分形式，若在无源空间，其空间中的煤质是各向同性、线性和均匀的，即煤质的参数不随时间变化且各向同性，则Maxwell旋度方程可写成： 式中，E是电场强度，单位为伏/米（V/m）；H是磁场强度，单位为安/米（A/m）；表示介质介电系数，单位为法拉/米（F/m）；表示磁导系数，单位为亨利/米（H/m）；表示介质电导率，单位为西门子/米（S/m）；m表示导磁率，单位为欧姆/米（/m）。  K.S.Yee在1966年

6、建立了如图2.1所示的空间网格，这就是著名的Yee氏元胞网格。图2.1 Yee氏网格及其电磁场分量分布 电场和磁场被交叉放置,电场分量位于网格单元每条棱的中心,磁场分量位于网格单元每个面的中心,每个磁场(电场)分量都有4个电场(磁场)分量环绕。这样不仅保证了介质分界面上切向场分量的连续性条件得到自然满足,而且还允许旋度方程在空间上进行中心差分运算,同时也满足了法拉第电磁感应定律和安培环路积分定律,也可以很恰当地模拟电磁波的实际传播过程。2.2 时域有限差分法相关技术2.2.1数值稳定性问题FDTD方程是一种显式差分方程，在执行时，存在一个重要的问题：即算法的稳定性问题。这种不稳定性表现为在解显

7、式方程时，随着时间步数的继续增加，计算结果也将无限制地增加。Taflove于1975年对Yee氏差分格式的稳定性进行了讨论，并导出了对时间步长的限制条件。数值解是否稳定主要取决于时间步长t与空间步长x、y、z的关系。对于非均匀媒质构成的计算空间选用如下的稳定性条件：上式是空间和时间离散之间应当满足的关系，又称为Courant稳定性条件。若采用均匀立方体网格：x=y=z=sD其中，v为计算空间中的电磁波的最大速度。2.2.2数值色散 FDTD方程组是对Maxwell旋度方程进行差分近似，在进行数值计算时，将会在计算网格中引起数字波模的色散，即在FDTD网格中，电磁波的相速与频率有关，电磁波的相速

8、度随波长、传播方向及变量离散化的情况不同而改变。这种关系由非物理因素引起，且色散将导致非物理因素引起的脉冲波形畸变、人为的各向异性和虚假折射等现象。为获得理想的色散关系，问题空间分割应按照小于正常网格的原则进行。一般选取的最大空间步长为max=min/20，min为所研究范围内电磁波的最小波长。由上分析说明，数值色散在用FDTD法分析电磁场传播中的影响是不可能避免的，但我们可以尽可能的减小数值色散的影响。2.2.3离散网格的确定 无论是简单目标还是复杂目标，在进行FDTD离散时网格尺寸的确定，除了受计算资源的限制不可能取得很小外，还需要考虑以下几个因素： 1目标离散精确度的要求。网格

9、应当足够小以便能精确模拟目标几何形状和电磁参数。 2FDTD方法本身的要求。主要是考虑色散误差的影响。设网格为立方体x=y=z=，所关心频段的频率上限为f max，对应波长为f minl，则考虑FDTD的数值色散要求min /N通常N10。上式是根据已知所关心频率上限情况下来确定FDTD网格尺寸d的；反之，若给定，则FDTD计算结果可用的上限频率也随之确定。2.3 吸收边界条件由时域有限差分法的基本原理可知，在利用时域有限差分法研究电磁场时，需在全部问题空间建立Yee氏网格空间，并存储每个单元网格上任一时间步的六个场分量用于下一时间步的计算。而在对于辐射、散射这类开放系统的实际研究中

10、，不可能有无限大的存储空间。因此，必须在某处将网格空间截断，且在截断边界网格点处运用特殊的场分量计算方法，使得向边界面行进的波在边界处保持“外向行进”特性、无明显的反射现象，并且不会使内部空间的场产生畸变，从而用有限网格空间模拟电磁波在无界空间中传播的情况。具有这种功能的边界条件称之为吸收边界条件，或辐射边界条件，或网格截断条件，如图2.2所示图 2.2 附加截断边界使计算区域变为有限域从FDTD的基本差分方程组可以看出，在截断边界面上切向场分量的计算需要利用计算空间以外的电磁场分量，因此FDTD基本差分方程对这些截断边界面上的场分量失效。如何处理截断边界上的场分量，使之与

11、需要考虑的无限空间有尽量小的差异，是FDTD中必须很好解决的一个重要问题。实际上，这是要求在误差可容忍的范围内，计算空间中的外向波能够顺利通过截断边界面而不引起波的明显反射，使有限计算空间的数值模拟与实际情况趋于一致，对外向波而言，就像在无限大空间中传播一样。所以，需要一种截断边界网格处的特殊计算方法，它不仅要保证边界场计算的必要精度，而且还要大大消除非物理因素引起的波反射，使得用有限的网格空间就能模拟电磁波在无限空间中的传播。但是如果处理不当，截断边界面可能造成较大反射，构成数值模拟误差的一部分，甚至可能造成算法不稳定。加于截断边界场分量符合上述要求的算法就称为吸收边界条件（Absorbin

12、g Boundary Conditions）。2.4 FDTD计算所需时间步的估计 为了使计算达到稳定,通常计算所需要时间步按照电磁波往返穿越FDTD计算区对角线35次来估计。若FDTD计算区总元胞数为N，则对角线上元胞约为（三维)。按照Courant稳定条件，设计算中心区t=/（2c），即穿越对角线一次需要时间步为。总计算时间步约为 需() 。对于二维情况则约为 ´ 。或者说，计算时间步大约等于FDTD计算区对角线上元胞数目的1220倍。实际上，计算所需时间步还与散射体具体形状、结构有关。图2.3给出了应用FDTD分析电磁场问题时的程序流程图图2.3 

13、FDTD 程序流程图3 课程设计要求及MATLAB语言实现3.1课程设计要求Ø 采用总场/散射场技术，在二维自由空间TEz网格中，实现脉冲平面波。Ø 总场区域100*100个网格。Ø 周边散射场区域20个网格。Ø 散射场区域外采用PML吸收边界条件。Ø 平面波0度角入射。使用高斯脉冲，中心频率为1GHz,带宽为1/e。Ø 采用正方形Yee网格，离散步长=1.5cm，时间离散步长为=/2c.3.2 MATLAB程序% 此程序实现自由空间的平面波 % 此程序采用总场/散射场技术 % 此程序使用PML设置吸收边界条件 %KE=12

14、8;%真空区域网格数IE=KE;JE=KE;T=0; NSTEP=300;%迭代次数e=2.71828;% 初始化 %dz=zeros(KE,KE);ez=zeros(KE,KE);hx=zeros(KE,KE);hy=zeros(KE,KE);ihx=zeros(KE,KE);ihy=zeros(KE,KE);ga=ones(KE,KE);ez_inc=zeros(1,JE);hx_inc=zeros(1,JE);% PML 参数 %gi2=ones(1,IE);gi3=gi2;fi1=zeros(1,IE);fi2=ones(1,IE);fi3=fi2;gj2=ones(1,JE);gj3

15、=gj2;fj1=zeros(1,JE);fj2=ones(1,JE);fj3=fj2;npml=8;for i=1:8 %设置PML层中的参数 xnum=npml-i; xd=npml; xxn=xnum/xd; xn=0.33*xxn3; gi2(i)=1.0/(1.0+xn); gi2(IE-1-i)=1.0/(1.0+xn); gi3(i)=(1.0-xn)/(1.0+xn); gi3(IE-1-i)=(1.0-xn)/(1.0+xn); xxn=(xnum-0.5)/xd; xn=0.25*xxn3; fi1(i)=xn; fi1(IE-2-i)=xn; fi2(i)=1.0/(1.

16、0+xn); fi2(IE-2-i)=1.0/(1.0+xn); fi3(i)=(1.0-xn)/(1.0+xn); fi3(IE-2-i)=(1.0-xn)/(1.0+xn); endfor j=1:8 xnum=npml-j; xd=npml; xxn=xnum/xd; xn=0.33*xxn3; gj2(j)=1.0/(1.0+xn); gj2(JE-1-j)=1.0/(1.0+xn); gj3(j)=(1.0-xn)/(1.0+xn); gj3(JE-1-j)=(1.0-xn)/(1.0+xn); xxn=(xnum-0.5)/xd; xn=0.25*xxn3; fj1(j)=xn;

17、fj1(JE-2-j)=xn; fj2(j)=1.0/(1.0+xn); fj2(JE-2-j)=1.0/(1.0+xn); fj3(j)=(1-xn)/(1.0+xn); fj3(JE-2-j)=(1.0-xn)/(1.0+xn); end% 总场/散射场边界 %ia=28;ib=IE-ia-1;ja=28;jb=JE-ja-1;ez_inc_low_m1=0;ez_inc_low_m2=0;ez_inc_high_m1=0;ez_inc_high_m2=0;% 信号源 %t0=80;%脉冲高度spread=12;%脉冲宽度% 网格 %ddx=0.015; % 离散步长 %dt=ddx/(2

18、*3e8); % 计算时间离散步长 %epsz=8.851*1e-12; % 自由空间的介电常数 % 迭代求解电场与磁场 % for n=1:NSTEP a=1; T=T+1; for j=2:JE ez_inc(j)= ez_inc(j)+0.5*( hx_inc(j-1)- hx_inc(j); end % 入射波缓冲区 % ez_inc(1) =ez_inc_low_m2; ez_inc_low_m2=ez_inc_low_m1; ez_inc_low_m1=ez_inc(2); ez_inc(JE-1) =ez_inc_high_m2; ez_inc_high_m2=ez_inc_hi

19、gh_m1; ez_inc_high_m1=ez_inc(JE-2); % 计算 dx 区域 % for i=2:KE for j=2:KE dz(i,j)=gi3(i)*gj3(j)*dz(i,j)+gi2(i)*gj2(j)*0.5*(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1); end end % 输入高斯脉冲信号源 % pulse=exp(-0.5*(t0-T)2/spread2); % ez_inc(3)=pulse; % 入射波 DZ 参数 % for i=ia:ib dz(i,ja)= dz(i,ja)+0.5*hx_inc(ja-1); dz(i,j

20、b)= dz(i,jb)-0.5*hx_inc(jb); end % 计算电场dx区域 % for i=1:KE for j=1:KE ez(i,j)=ga(i,j)*dz(i,j); end end % 设置电场 EZ 边缘为0, 作为PML的参数% for j=1:JE-1 ez(1,j)=0; ez(IE,j)=0; end for i=1:IE-1 ez(i,1)=0; ez(i,JE)=0; end for j=1:JE-1 hx_inc(j)=hx_inc(j)+0.5*( ez_inc(j)-ez_inc(j+1); end % 计算磁场hx区域 % for i=1:KE for

21、 j=1:KE-1 curl_e=ez(i,j)-ez(i,j+1); ihx(i,j)=ihx(i,j)+fi1(i)*curl_e; hx(i,j)=fj3(j)*hx(i,j)+fj2(j)*0.5* (curl_e+ihx(i,j); end end % 入射磁场Hx区域参数 % for i=ia:ib hx(i,ja-1)= hx(i,ja-1)+0.5*ez_inc(ja); hx(i,jb)= hx(i,jb)-0.5*ez_inc(jb); end % 计算磁场hy区域% for i=1:KE-1 for j=1:KE curl_e=ez(i+1,j)-ez(i,j); ihy

22、(i,j)=ihy(i,j)+fj1(j)*curl_e; hy(i,j)=fi3(i)*hy(i,j)+fi2(i)*0.5* (curl_e+ihy(i,j); end end % 入射磁场hy区域参数 % for j=ja:jb hy(ia-1,j)= hy(ia-1,j)-0.5*ez_inc(j); hy(ib,j)= hy(ib,j)+0.5*ez_inc(j); end % 图形展示 % m=moviein(500); x=1:KE; y=1:KE; X,Y=meshgrid(x,y); surf(X,Y,ez); axis(0 KE 0 KE -1 1 ); set(gcf,'Color', 'white', 'Number', 'off', 'Name', sprintf('Simulation FDTD 2D, Iteration = %i', n); title( sprintf('t = %.3f nsec',n*dt*1e9) xlabel('x'); ylabel('y'); zlabel('Ez'); m(a)=getframe(gcf); a=a+1; end 下图为该