时域有限差分法论文.doc

时域有限差分法 1 选题背景在多种可用的数值方法中，时域有限差分法(fdtd)是一种新近发展起来的可选方法。1966年，k.s.yee首次提出电磁场数值计算的新方法时域有限差分法(finite difference- time domain，简称fdtd)。经历了二十年的发展fdtd法才逐渐走向成熟。上世纪80年代后期以来fdtd法进入了一个新的发展阶段，即由成熟转为被广泛接受和应用的阶段。fdtd法是解决复杂问题的有效方法之一，是一种直接基于时域电磁场微分方程的数值算法，它直接在时域将maxwell旋度方程用二阶精度的中心差分近似，从而将时域微分方程的求解转换为差分方程的迭代求解。是电磁场和电磁波运动规律和运动过程的计算机模拟。原则上可以求解任意形式的电磁场和电磁波的技术和工程问题，并且对计算机内存容量要求较低、计算速度较快、尤其适用于并行算法。现在fdtd法己被广泛应用于天线的分析与设计、目标电磁散射、电磁兼容、微波电路和光路时域分析、生物电磁剂量学、瞬态电磁场研究等多个领域1。经过了近四十年的发展，fdtd法在计算方法和应用上取得了大量成果。近几年来，讨论fdtd法的深入发展和实际应用的文章几乎按指数增长，目前fdtd法的许多重要问题得到了很好的解决，已经发展成为一种成熟的数值计算方法。随着计算机数据处理性能的快速提高和计算机价格的下降，使得fdtd法的应用范围越来越广，而fdtd法本身在应用中又有新的发展.2 原理分析2.1 fdtd的yee元胞e,h场分量取样节点在空间和时间上采取交替排布，利用电生磁，磁生电的原理 图1 yee模型如图1所示，yee单元有以下特点2：1）e与h分量在空间交叉放置，相互垂直；每一坐标平面上的e分量四周由h分量环绕，h分量的四周由e分量环绕；场分量均与坐标轴方向一致。2）每一个yee元胞有8个节点，12条棱边，6个面。棱边上电场分量近似相等，用棱边的中心节点表示，平面上的磁场分量近似相等，用面的中心节点表示。3）每一场分量自身相距一个空间步长，e和h相距半个空间步长4）每一场分量自身相距一个时间步长，e和h相距半个时间步长，电场取n时刻的值，磁场取n+0.5时刻的值；即：电场n时刻的值由n-1时刻的值得到，磁场n+0.5时刻的值由n-0.5时刻的值得到；电场n时刻的旋度对应n+0.5时刻的磁场值，磁场n+0.5时刻的旋度对应(n+0.5)+0.5时刻的电场值，逐步外推。5）3个空间方向上的时间步长相等， 以保证均匀介质中场量的空间变量与时间变量完全对称。应用这种离散方式，将含时间变量的maxwell方程转化为一组差分方程，并在时间轴上逐步推进地求解空间电磁场。由电磁问题的初值和边界条件，就可以逐步推进地求解以后各时刻空间电磁场分布。2.2 maxwell方程fdtd的差分格式麦克斯韦第一、二方程 （1）式中，时电流密度，反映电损耗，是磁流密度，单位，反映磁损耗。主要与上式对应。各向同性介质中的本构关系： （2）其中是磁阻率，计算磁损耗的。以为变量，在直角坐标中，展开麦克斯韦第一、二方程，分别为 （3） （4）令代表在直角坐标中的任何一个分量，离散符号取为 （5）关于时间和空间的一阶偏导数取中心差分近似为 （6）可以看出，每一节点上沿某一方向场分量的一阶偏微分可以用在该方向上相邻两点的一阶中心差商来描述，将式(1)用一阶中心差商方程取代，整理后便得到一阶差分方程，它具有二阶精度3。 yee元胞如图1所示，规定为 1）剖分节点与场分量所在棱边中点不同，场分量的位置，即节点是yee元胞节点的相对位置，不需要单独编码；2）当空间存在媒质分界面时，场量自动满足场的连续性条件，电磁分量的取样方式不仅符合法拉第电磁感应定律和安培环路定律的自然结构，也符合麦克斯韦方程的差分计算。其次，时间步长可以取为电磁波传播一个空间步长所需时间的一半，因此与在时间顺序上交替抽样，时间间隔相差半个时间步长。2.3 一维问题 均匀平面波（tem波）是一维问题，设电磁波沿z轴方向传播，则，场量和介质参数均与x，y无关，即，麦克斯韦方程为 （7）和 （8）旋转坐标轴后可以只保留一组公式4，设保留（7）yee元胞如图2所示 图2 一维yee元胞差分格式为 （9） （10） 如果介质无损耗，则2.4 二维问题三维通常是散射问题，二维是te、tm波问题，一维是tem波问题。在二维场中，所有物理量与z坐标无关，既。于是在te和tm波的表达式分别为te波（） （11）tm波（） （12）图3分别给出了tm波和te波的yee元胞图 图3 tm波的yee元胞 图4 te波的yee元胞对于te波，只要令，在上，不随z变化，m中去掉k即可得到：式中： （13）式中： （14）式中， （15） 对tm波，只要令，在上，不随z变化，m中去掉k，即可得到：式中， （16）式中， （17）式中： （18） 为了编写统一的te和tm 波二维fdtd程序，可将描述te波差分公式(13)(15)中相应的标号整体移动1/2，即坐标（x,y）分别沿x和y轴方向移动半个网格，并将离散时间也移动半个时间步长，式(13)(15)可以重新写为 式中： （19）式中： （20）式中， （21）可以看出，te波的fdtd公式(19)(21)与tm波的fdtd公式(16)(18)形式相同，给编程带来极大方便。注意te波和tm波之间的对偶关系5，即这样就可以编写统一的计算程序了。2.5 三维问题（直角坐标系）2.5.1电场时间推进差分格式节点的3个电场分量分别用、位置上的表示，以式（3）中第一个公式为例： 在时间步，对节点的离散公式为：上式中的第二项用平均值来替代是因为离散方程中电场的时间取样是整数n，磁场的时间取样是n+1/2，所以只能取n及n+1时电场的平均值。实际也证明这个平均值使fdtd算法具有数值稳定性。整理后，将作为未知数，其余作为迭代计算的已知数 （22） 同理，式（3）中其它两个公式的离散形式为 （23） （24）以上三式是电场的时间推进计算公式。2.5.2磁场时间推进差分格式节点的3个磁场分量分别用、位置上的表示，同样，讨论式（4）中第一个公式，设观察点为的节点，即在时刻，对节点的离散公式为： （25） 同理，式（4）中其它两个公式的离散形式为 （26） （27）以上三式是磁场的时间推进计算公式。 时域推进计算框图（交叉半步逐步推进）若已知时空间各节点处的电场值（赋初值）计算时空间各节点处的磁场值，式（25）（27）计算时空间各节点处的电场值，式（22）（24）在编程中，为了使电场和磁场有相同的数量级（为减小误差），可对h或e进行“归一化”处理，即：用取代，用取代，式中是自由空间波阻抗。计算结果再分别除以和乘以即可。可以看出，这种离散方法电场和磁场在时间顺序上交替抽样，抽样间隔相差半个时间步长，使麦克斯韦方程离散后成为显示差分方程，从而可以在时间上迭代求解，不需矩阵求逆。给定初值后，可以逐步推进，求得以后各个时刻点的空间电磁分布。这是fdtd法的最大特点。2.6 解的稳定性在fdtd中，时间增量和空间增量、之间不是相互独立的，它们的取值必须满足一定的关系，以避免数值结果的不稳定，表现为随着时间步数的增加，计算结果发散。造成解不稳定的因素有多种： 误差因素：计算机在计算过程中，原始数据可能有误差，如系数阵建立过程中产生的误差，而每次运算由于只能保留有限位数而又产生误差，误差的积累有可能淹没真正解，使计算结果不可靠，即不稳定； 计算方法不合适； 、离散间隔不当等。为了确定数值解稳定的条件，有许多推导方式，结论相同6。2.6.1时间步长稳定性要求一般情况下 （28）2.6.2时间步长与空间步长的关系三维 （29）在非均匀区域，v取最大值。真空中v=c(光速)。若是正方体yee元胞，那么 （30）若是正方形yee元胞（二维），那么 （31）若是线段等分yee元胞（一维），那么 （32）2.7 数值色散当波传播的速度是频率的函数，即速度与频率有关时，称其波为色散波。色散的原因有多种：由于媒质是金属或各向异性等；由于载波体形状（也称几何色散）；由于数值计算方法等（也称数值色散）原因。用差分法计算时，会在计算网格中引起模拟波模的色散，即在时域有限差分网格中，数值波模的传播速度随频率变化，即速度在网格中随数值波模在网格中的传播方向以及离散化情况的不同而改变。这种色散将导致非物理因素引起的脉冲波形畸变，人为的各向异性及虚假折射现象，都会给计算带来误差，因此要尽量减小数值色散7。为了减小数值色散，除了满足式（31）外，可取比式（30）更高的要求 （33）并令电磁波沿网格的对角线方向传播，这时有，可以得到理想的色散关系。3 程序设计及模拟结果3.1一维fdtd算法及模拟结果%*% 1-d fdtd code with simple radiation boundary conditions%*% this matlab m-file implements the finite-difference time-domain% solution of maxwells curl equations over a one-dimensional space% lattice comprised of uniform grid cells.% to illustrate the algorithm, a sinusoidal wave (1ghz) propagating % in a nonpermeable lossy medium (epsr=1.0, sigma=5.0e-3 s/m) is % modeled. the simplified finite difference system for nonpermeable% media (discussed in section 3.6.6 of the text) is implemented.% the grid resolution (dx = 1.5 cm) is chosen to provide 20% samples per wavelength. the courant factor s=c*dt/dx is set to% the stability limit: s=1. in 1-d, this is the magic time step.% the computational domain is truncated using the simplest radiation% boundary condition for wave propagation in free space: % ez(imax,n+1) = ez(imax-1,n)% to execute this m-file, type fdtd1d at the matlab prompt.% this m-file displays the fdtd-computed ez and hy fields at every % time step, and records those frames in a movie matrix, m, which is% played at the end of the simulation using the movie command.%*clear%*% fundamental constants%*cc=2.99792458e8; %speed of light in free spacemuz=4.0*pi*1.0e-7; %permeability of free spaceepsz=1.0/(cc*cc*muz); %permittivity of free spacefreq=1.0e+9; %frequency of source excitationlambda=cc/freq; %wavelength of source excitationomega=2.0*pi*freq;%*% grid parameters%* ie=200; %number of grid cells in x-direction ib=ie+1; dx=lambda/20.0; %space increment of 1-d latticedt=dx/cc; %time stepomegadt=omega*dt;nmax=round(12.0e-9/dt); %total number of time steps%*% material parameters%*eps=1.0;sig=5.0e-3;%*% updating coefficients for space region with nonpermeable media%* scfact=dt/muz/dx; ca=(1.0-(dt*sig)/(2.0*epsz*eps)/(1.0+(dt*sig)/(2.0*epsz*eps);cb=scfact*(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps); %*% field arrays%*ez(1:ib)=0.0;hy(1:ie)=0.0;%*% movie initialization%* x=linspace(dx,ie*dx,ie); subplot(2,1,1),plot(x,ez(1:ie)/scfact,r),axis(0 3 -1 1);ylabel(ez); subplot(2,1,2),plot(x,hy,b),axis(0 3 -3.0e-3 3.0e-3);xlabel(x (meters);ylabel(hy); rect=get(gcf,position);rect(1:2)=0 0; m=moviein(nmax/2,gcf,rect); %*% begin time-stepping loop%* for n=1:nmax %*% update electric fields%* ez(1)=scfact*sin(omegadt*n); rbc=ez(ie);ez(2:ie)=ca*ez(2:ie)+cb*(hy(2:ie)-hy(1:ie-1);ez(ib)=rbc; %*% update magnetic fields%* hy(1:ie)=hy(1:ie)+ez(2:ib)-ez(1:ie); %*% visualize fields%* if mod(n,2)=0; rtime=num2str(round(n*dt/1.0e-9); subplot(2,1,1),plot(x,ez(1:ie)/scfact,r),axis(0 3 -1 1);title(time = ,rtime, ns);ylabel(ez); subplot(2,1,2),plot(x,hy,b),axis(0 3 -3.0e-3 3.0e-3);title(time = ,rtime, ns);xlabel(x (meters);ylabel(hy); m(:,n/2)=getframe(gcf,rect); end %*% end time-stepping loop%* end movie(gcf,m,0,10,rect);3.2二维fdtd算法及模拟结果%*% 2-d fdtd te code with pml absorbing boundary conditions%*% this matlab m-file implements the finite-difference time-domain% solution of maxwells curl equations over a two-dimensional% cartesian space lattice comprised of uniform square grid cells.% to illustrate the algorithm, a 6-cm-diameter metal cylindrical % scatterer in free space is modeled. the source excitation is % a gaussian pulse with a carrier frequency of 5 ghz.% the grid resolution (dx = 3 mm) was chosen to provide 20 samples% per wavelength at the center frequency of the pulse (which in turn% provides approximately 10 samples per wavelength at the high end% of the excitation spectrum, around 10 ghz).% the computational domain is truncated using the perfectly matched% layer (pml) absorbing boundary conditions. the formulation used % in this code is based on the original split-field berenger pml. the% pml regions are labeled as shown in the following diagram: % -% | | back pml | |% -% |l | /| r|% |e | (ib,jb) | i|% |f | | g|% |t | | h|% | | main grid | t|% |p | | |% |m | | p|% |l | (1,1) | m|% | |/ | l|% -% | | front pml | |% -% to execute this m-file, type fdtd2d at the matlab prompt.% this m-file displays the fdtd-computed ex, ey, and hz fields at % every 4th time step, and records those frames in a movie matrix, % m, which is played at the end of the simulation using the movie % command.%*clear %*% fundamental constants%* cc=2.99792458e8; %speed of light in free spacemuz=4.0*pi*1.0e-7; %permeability of free spaceepsz=1.0/(cc*cc*muz); %permittivity of free space freq=5.0e+9; %center frequency of source excitationlambda=cc/freq; %center wavelength of source excitationomega=2.0*pi*freq; %*% grid parameters%* ie=100; %number of grid cells in x-directionje=50; %number of grid cells in y-direction ib=ie+1;jb=je+1; is=15; %location of z-directed hard sourcejs=je/2; %location of z-directed hard source dx=3.0e-3; %space increment of square latticedt=dx/(2.0*cc); %time step nmax=300; %total number of time steps iebc=8; %thickness of left and right pml regionjebc=8; %thickness of front and back pml regionrmax=0.00001;orderbc=2;ibbc=iebc+1;jbbc=jebc+1;iefbc=ie+2*iebc;jefbc=je+2*jebc;ibfbc=iefbc+1;jbfbc=jefbc+1; %*% material parameters%* media=2; eps=1.0 1.0;sig=0.0 1.0e+7;mur=1.0 1.0;sim=0.0 0.0; %*% wave excitation%* rtau=160.0e-12;tau=rtau/dt;delay=3*tau; source=zeros(1,nmax);for n=1:7.0*tau source(n)=sin(omega*(n-delay)*dt)*exp(-(n-delay)2/tau2);end %*% field arrays%* ex=zeros(ie,jb); %fields in main grid ey=zeros(ib,je);hz=zeros(ie,je); exbcf=zeros(iefbc,jebc); %fields in front pml regioneybcf=zeros(ibfbc,jebc);hzxbcf=zeros(iefbc,jebc);hzybcf=zeros(iefbc,jebc); exbcb=zeros(iefbc,jbbc); %fields in back pml regioneybcb=zeros(ibfbc,jebc);hzxbcb=zeros(iefbc,jebc);hzybcb=zeros(iefbc,jebc); exbcl=zeros(iebc,jb); %fields in left pml regioneybcl=zeros(iebc,je);hzxbcl=zeros(iebc,je);hzybcl=zeros(iebc,je); exbcr=zeros(iebc,jb); %fields in right pml regioneybcr=zeros(ibbc,je);hzxbcr=zeros(iebc,je);hzybcr=zeros(iebc,je); %*% updating coefficients%* for i=1:media eaf =dt*sig(i)/(2.0*epsz*eps(i); ca(i)=(1.0-eaf)/(1.0+eaf); cb(i)=dt/epsz/eps(i)/dx/(1.0+eaf); haf =dt*sim(i)/(2.0*muz*mur(i); da(i)=(1.0-haf)/(1.0+haf); db(i)=dt/muz/mur(i)/dx/(1.0+haf);end%*% geometry specification (main grid)%*% initialize entire main grid to free space caex(1:ie,1:jb)=ca(1); cbex(1:ie,1:jb)=cb(1); caey(1:ib,1:je)=ca(1);cbey(1:ib,1:je)=cb(1); dahz(1:ie,1:je)=da(1);dbhz(1:ie,1:je)=db(1); % add metal cylinder diam=20; % diameter of cylinder: 6 cmrad=diam/2.0; % radius of cylinder: 3 cmicenter=4*ie/5; % i-coordinate of cylinders centerjcenter=je/2; % j-coordinate of cylinders center for i=1:iefor j=1:je dist2=(i+0.5-icenter)2 + (j-jcenter)2; if dist2 = rad2 caex(i,j)=ca(2); cbex(i,j)=cb(2); end dist2=(i-icenter)2 + (j+0.5-jcenter)2; if dist2 = rad2 caey(i,j)=ca(2); cbey(i,j)=cb(2); endendend %*% fill the pml regions%* delbc=iebc*dx;sigmam=-log(rmax/100.0)*epsz*cc*(orderbc+1)/(2*delbc);bcfactor=eps(1)*sigmam/(dx*(delbcorderbc)*(orderbc+1); % front region caexbcf(1:iefbc,1)=1.0;cbexbcf(1:iefbc,1)=0.0;for j=2:jebc y1=(jebc-j+1.5)*dx; y2=(jebc-j+0.5)*dx; sigmay=bcfactor*(y1(orderbc+1)-y2(orderbc+1); ca1=exp(-sigmay*dt/(epsz*eps(1); cb1=(1.0-ca1)/(sigmay*dx); caexbcf(1:iefbc,j)=ca1; cbexbcf(1:iefbc,j)=cb1;endsigmay = bcfactor*(0.5*dx)(orderbc+1);ca1=exp(-sigmay*dt/(epsz*eps(1);cb1=(1-ca1)/(sigmay*dx);caex(1:ie,1)=ca1;cbex(1:ie,1)=cb1;caexbcl(1:iebc,1)=ca1;cbexbcl(1:iebc,1)=cb1;caexbcr(1:iebc,1)=ca1;cbexbcr(1:iebc,1)=cb1; for j=1:jebc y1=(jebc-j+1)*dx; y2=(jebc-j)*dx; sigmay=bcfactor*(y1(orderbc+1)-y2(orderbc+1); sigmays=sigmay*(muz/(epsz*eps(1); da1=exp(-sigmays*dt/muz); db1=(1-da1)/(sigmays*dx); dahzybcf(1:iefbc,j)=da1; dbhzybcf(1:iefbc,j)=db1; caeybcf(1:ibfbc,j)=ca(1); cbeybcf(1:ibfbc,j)=cb(1); dahzxbcf(1:iefbc,j)=da(1); dbhzxbcf(1:iefbc,j)=db(1);end % back region caexbcb(1:iefbc,jbbc)=1.0;cbexbcb(1:iefbc,jbbc)=0.0;for j=2:jebc y1=(j-0.5)*dx; y2=(j-1.5)*dx; sigmay=bcfactor*(y1(orderbc+1)-y2(orderbc+1); ca1=exp(-sigmay*dt/(epsz*eps(1); cb1