时域有限差分法对平面TE波的MATLAB仿真

时域有限差分法对平面TE波的MATLAB仿真姓名：王云璐学号：2011302021（西北工业大学电子信息学院08041103，陕西西安，710072）摘要：本文分析FDTD算法的基本原理及两种典型边界条件的算法特点，通过MATLAB程序对TE波进行了仿真，模拟了高斯制下完全匹配层中磁场分量瞬态分布。得到了相应的磁场幅值效果图。最后得出用Matlab语言对FDTD算法编程的几点结论。关键词：FDTD；MATLAB；PML；平面TE波1引言电磁场的有限差分一般是在时域上进行的，随着计算机技术的发展和应用，近年来，时域计算方法越来越受到重视。时域有限差分法具有简单、结果直观、网格剖分简单等优点。近些年FDTD发展的十分迅速，在许多方面都有重要应用，包括天线设计，微波电路设计，电磁兼容分析，电磁散射计算，光子学应用等等。时域有限差分(FDTD)算法是K.S.Yee于1966年提出的直接对麦克斯韦方程作差分处理，用来解决电磁脉冲在电磁介质中传播和反射问题的算法。基本思想是：FDTD计算域空间节点采用Yee元胞的方法，同时电场和磁场节点空间与时间上都采用交错抽样；把整个计算域划分成包括散射体的总场区以及只有反射波的散射场区，这两个区域是以连接边界相连接，最外边是采用特殊的吸收边界，同时在这两个边界之间有个输出边界，用于近、远场转换；在连接边界上采用连接边界条件加入入射波，从而使得入射波限制在总场区域；在吸收边界上采用吸收边界条件，尽量消除反射波在吸收边界上的非物理性反射波。本文主要简述了FDTD算法的基本原理，解的稳定性以及边界条件特点，并用Matlab语言进行对平面TE波进行了编程计算。2FDTD基本原理时域有限差分法的主要思想是把Maxwell方程在空间、时间上离散化，用差分方程代替一阶偏微分方程，求解差分方程组，从而得出各网格单元的场值。FDTD空间网格单元上电场和磁场各分量的分布如图1所示。图1.Yee氏网格及其电磁场分量分布电场和磁场被交叉放置，电场分量位于网格单元每条棱的中心，磁场分量位于网格单元每个面的中心，每个磁场(电场)分量都有4个电场(磁场)分量环绕。这样不仅保证了介质分界面上切向场分量的连续性条件得到自然满足，而且还允许旋度方程在空间上进行中心差分运算，同时也满足了法拉第电磁感应定律和安培环路积分定律，也可以很恰当地模拟电磁波的实际传播过程。3解的稳定性FDTD方程是一种显式差分方程，在执行时，存在一个重要的问题：即算法的稳定性问题。这种不稳定性表现为在解显式方程时，随着时间步数的继续增加，计算结果也将无限制地增加。Taflove等于1975年对Yee氏差分格式的稳定性进行了讨论，并导出了对时间步长的限制条件。数值解是否稳定主要取决于时间步长t与空间步长x、y、z的关系。对于非均匀媒质构成的计算空间选用如下的稳定性条件：2221111()()()tvxyz(1)(1)式是空间和时间离散之间应当满足的关系，又称为Courant稳定性条件。若采用均匀立方体网格：xyzs3stv(2)其中，v为计算空间中的电磁波的最大速度。4FDTD算法的几个边界条件边界条件是对FDTD算法应用的另外一个关键问题。时域有限差分法是在计算机的数据存储空间中对连续的实际电磁波的传播过程在时间按空进程上进行数字模拟。在计算电磁场的辐射、散射等问题中，边界总是开放的，电磁场从实际的角度看总是占据无限大空间，由于计算机内存总是有限的，只能模拟有限的空间。也就是说，时域有限差分网格将在某处被截断。如何处理截断边界，使之需要考虑的无限空间的差异尽量的小，这是FDTD算法要解决的重要问题。在当今FDTD算法中，涉及到的边界条件主要有Mur吸收边界条件、廖氏吸收边界条件、超吸收边界条件、完全匹配层(PML)，在本文中主要介绍Mur吸收边界条件和完全匹配层(PML)。4.1Mur吸收边界条件其具体推导过程可参考文献[4]，递推公式如下一阶条件下的递推公式11,1,[1,,]nnnnzzzzctEijEijEijEijct(3)二阶条件下Mur的FDTD递推公式1121/21/21/21/2,1,1[1,1,]111,,1,2()22211,2nnnnzzzznnnxxxnxctEijEijEijEijctctHijHijHijctHij(4)这就是Mur所建议的具有二阶近似的、适用于二维问题的近似吸收边界条件。它在FDTD中有着广泛的应用。4.2PML吸收边界条件1994年，Berengert提出了用完全匹配(PML)来吸收外向电磁波，它将电磁场分量在吸收边界分裂，并能分别对各个分裂的场分量赋予不同的损耗。完全匹配层是一种非物理性的电磁波吸收层，有关性质和在具体条件下PML格式的FDTD递推公式见文献[5]，本文采用PML吸收边界条件。5MATLAB的仿真程序及结果5.1MATLAB程序及相应说明%二维FDTDTE波仿真Clearall；%定义常数pi=3.1415；c=3.0e10；%高斯制下光速f=1.0e15；%频率lambda=c/f；%波长nmax=400；%时间步数del_s=lambda/20；%每最小波长20个采样点del_t=0.5*del_s/c；%迭代时间步长n=182；%真空区域网格数np=9；%pml层数N1=n+2*np；%总网格数N=N1+1；%采样点数M=4；%导电率渐变指数sigma_max=(M+1)/1.50/pi/del_s；%最大导电率%TE波的分量初始化tic；figure(1)；axis([0N0N-0.50.5])；Ex=zeros(N1,N)；%x方向为横向，采样点为网格的横向边，故行数+1Ey=zeros(N,N1)；%y方向为纵向，采样点为网格的纵向边，故列数+1Bz=zeros(N1,N1)；%矩阵行为纵向网格数，矩阵列为横向网格数,循环中用j表示行数，i表示列数Bzx=zeros(N1,N1)；Bzy=zeros(N1,N1)；Bzxx=zeros(nmax,2)；%进入电磁场迭代计算fortt=1:nmaxfori=1:N1ifi=np+1&&i=N1-npdi=0；elseifi=npdi=np-i+0.5；elseifi=N1-np+1di=np+i-N1-0.5；end%di是采样点横向距PML内边界的距离sigma_mx=sigma_max*(di/np)^M；forj=1:N1ifj=np+1&&j=N1-npdj=0；elseifj=npdj=np-j+0.5；elseifj=N1-np+1dj=np+j-N1-0.5；end%dj是采样点纵向距PML内边界的距离sigma_my=sigma_max*(dj/np)^M；ifsigma_mx+sigma_my==0%真空区ifj==100&&i==100t=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；%加入脉冲源elsec_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；endelse%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)；endendendfori=2:N1ifi=np+1&&i=N1-npdi=0；elseifi=npdi=np-i+1；elseifi=N1-np+1di=np+i-N1-1；end%di是采样点横向距PML内边界的距离sigma_ex=sigma_max*(di/np)^M；forj=1:N1cam=(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))；endendfori=1:N1forj=2:N1ifj=np+1&&j=N1-npdj=0；elseifj=npdj=np-j+1；elseifj=N1-np+1dj=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))；endendBzxx(tt,1)=Bz(90,50)；%对靠近边界的中央磁场点采样Bzxx(tt,2)=Bz(50,90)；figure(1)；%可视化处理clf；mesh(Bz)；%磁场的幅值axis([0N0N-0.50.5])；xlabel('i')ylabel('j')drawnow；endfigure(2)；%激励信号示意图plot(Bzxx)；figure(3)；%磁场的幅值分布图title('磁场幅值分布图')；surface(Bz)；shadinginterp；axissquare；toc5.2结果图2.正弦激励下3D磁场幅值分布图图3.正弦激励信号图4.正弦激励下磁场幅值分布图图5.高斯脉冲激励下3D磁场幅值分布图图6.高斯脉冲信号图7.高斯脉冲激励下磁场幅值分布图6结束语从以上分析可以看出，FDTD方法在复杂电磁场的数值计算中是一种比较实用可行的方法，可以帮助我们更加形象准确的理解场。FDTD方法是一种计算复杂目标电磁散射数值的有效方法，FDTD法在电磁场数值分析方面有很大的优越性，而Matlab具有强大的数据处理和图形处理功能，可以快速地编出高效高质量的程序。将二者的优势有效地结合起来，可以将算法迅速程序化，并获得很好的数据处理结果，使研究者可以集中精力在FDTD方法和研究对象本身上，而只需花费少量的时间在程序的实现上。参考文献[1]喻志远.MATLAB模拟的电磁学时域有限差分法[M].国防工业出版社,2012-8.[2]张亮,寇晓艳.FDTD算法的Matlab语言实现[J].延安大学学报,2009,28(2).[3]A.TaylorandM.E.Brodwin.Numericalsolutionofsteady-stateEMscatteringproblemsusingthetime-Maxwell’sequations.[J]IEEETrans.MicrowaveTheroyTech.,Aug.1975，MTT-