《毕业设计论文时域有限差分法对平面TE波的MATLAB仿真.doc》由会员分享，可在线阅读，更多相关《毕业设计论文时域有限差分法对平面TE波的MATLAB仿真.doc（34页珍藏版）》请在三一办公上搜索。
1、时域有限差分法对平面TE波的MATLAB仿真摘 要 时域有限差分法是由有限差分法发展出来的数值计算方法。自1966年Yee在其论文中首次提出时域有限差分以来，时域有限差分法在电磁研究领域得到了广泛的应用。主要有分析辐射条线、微波器件和导行波结构的研究、散射和雷达截面计算、分析周期结构、电子封装和电磁兼容的分析、核电磁脉冲的传播和散射以及在地面的反射及对电缆传输线的干扰、微光学元器件中光的传播和衍射特性等等。由于电磁场是以场的形态存在的物质，具有独特的研究方法，采取重叠的研究方法是其重要的特点，即只有理论分析、测量、计算机模拟的结果相互佐证，才可以认为是获得了正确可信的结论。时域有限差分法就是实
2、现直接对电磁工程问题进行计算机模拟的基本方法。在近年的研究电磁问题中，许多学者对时域脉冲源的传播和响应进行了大量的研究，主要是描述物体在瞬态电磁源作用下的理论。另外，对于物体的电特性，理论上具有几乎所有的频率成分，但实际上，只有有限的频带内的频率成分在区主要作用。 文中主要谈到了关于高斯制下完全匹配层的差分公式的问题，通过MATLAB程序对TE波进行了仿真，模拟了高斯制下完全匹配层中磁场分量瞬态分布。得到了相应的磁场幅值效果图。关键词：时域有限差分 完全匹配层 MATLAB 磁场幅值效果图目 录摘 要1目 录2第一章 绪 论31.1 课题背景与意义31.2 时域有限差分法的发展与应用32.1
3、Maxwell方程和Yee氏算法62.2 FDTD的基本差分方程82.3 时域有限差分法相关技术102.3.1 数值稳定性问题102.3.2 数值色散112.3.3 离散网格的确定122.4 吸收边界条件122.4.1 一阶和二阶近似吸收边界条件132.4.2 二维棱边及角顶点的处理162.4.3 完全匹配层182.5 FDTD计算所需时间步的估计22第三章 MATLAB的仿真的程序及模拟243.1 MATLAB程序及相应说明243.2 出图及结果273.2.1程序部分273.2.2 所出的效果图28第四章 结 论30参考文献31第一章 绪 论1.1 课题背景与意义20世纪60年代以来，随着计
4、算机技术的发展，一些电磁场的数值计算方法逐步发展起来，并得到广泛应用，其中主要有：属于频域技术的有限元法（FEM）、矩量法(MM)和单矩法等；属于时域技术方面的时域有限差分法(FDTD)、传输线矩阵法(TLM)和时域积分方程法等。此外，还有属于高频技术的几何衍射理论(GTD)和衍射物理理论(PLD)等。各种方法都具有自己的特点和局限性，在实际中经常把它们相互配合而形成各种混合方法12。其中FDTD是一种已经获得广泛应用并且有很大发展前景的时域数值计算方法。时域有限差分（FDTD）方法于1966年由K.S.Yee3 提出并迅速发展，且获得广泛应用。K.S.Yee用后来被称作Yee氏网格的空间离散
5、方式，把含时间变量的Maxwell旋度方程转化为差分方程，并成功地模拟了电磁脉冲与理想导体作用的时域响应。但是由于当时理论的不成熟和计算机软硬件条件的限制，该方法并未得到相应的发展。20世纪80年代中期以后，随着上述两个条件限制的逐步解除，FDTD便凭借其特有的优势得以迅速发展。它能方便、精确地预测实际工程中的大量复杂电磁问题，应用范围几乎涉及所有电磁领域，成为电磁工程界和理论界研究的一个热点。目前，FDTD日趋成熟，并成为分析大部分实际电磁问题的首选方法。另外，利用矩量法求解电磁场问题时，要用到并失Green函数。对于某些问题，可以找到其解析形式的并失Green函数；而对于复杂的问题，很难找
6、到其解析形式的并失Green函数，这样就使得问题无法解决。作为时域分析中的一个重要数值方法，FDTD不存在这样的问题。1.2 时域有限差分法的发展与应用经过四十多年的发展，FDTD已发展成为一种成熟的数值计算方法。在发展过程中，几乎都是围绕几个重要问题展开的，即数值稳定性、计算精度、数值色散、激励源技术以及开域电磁问题的吸收边界条件等。数值稳定和计算精度对任何一种数值计算方法都是至关重要的。A.Taylor和M.E.Brodwin4利用本征值方法给出了直角坐标系下FDTD的空间步长与时间步长之间的关系。X.Min等5研究了存在边界条件时FDTD的稳定性问题。对于数值色散，与实际的物理色散不同，
7、它是由电磁场量在空间和时间上的对波动方程作差分近似处理造成的。这种色散引起的误差造成在计算区域内传播的电磁波逐渐畸变67。K. L. Shlager 等8比较了二维和三维空间中几种正交网格算法的色散误差。当采用其他变形或非正交网格时，必须重新分析其数值稳定性和色散特性911，P.Monk 和 E.Suli12分析了不均匀长方体网格算法的稳定性。激励源的设计和引入也是FDTD的一个重要任务。目前，应用最广泛的激励源引入技术是总场/散射场体系12。对于散射问题，通常在FDTD计算空间中引入连接边界，它将整个计算空间划分为内部的总场区和外部的散射场区，如图1-1。利用Huygens原理，可以在连接边
8、界处引入入射场，使入射场的加入变得简单易行。图1-1开域电磁问题中，为了在有限的计算空间内模拟无限空间中的电磁问题，必须在计算空间的截断边界处设置吸收边界条件。吸收边界条件从开始简单的插值边界，已经发展了多种吸收边界条件。在早期得到广泛应用的是G.Mur13的一阶和二阶吸收边界条件，它是基于B.Engquist和A.Majda14的单向波方程而提出的差分格式，在FDTD仿真区域外边界具有0.5%到5%的反射系数。目前应用最广泛的是J.P.Berenger15-17的分裂式完全匹配层，以及Z.S.Sacks等18和S.D.Gedney20的各向异性介质的完全匹配层，它们可使FDTD模拟的最大动态
9、范围达到80dB。另一方面，为了更好的拟合研究对象的形状，克服台阶逼近带来的误差，D.E.Merewether19提出了柱坐标系下的网格剖分方法，R.Holland20提出了球坐标系下的网格剖分方法，P.Monk和E.Suli12提出了变网格步长方法，S.S.Zivanovic等21和P.Thoma等22提出了亚网格技术（即在一般区域采用粗网格，在电磁场快变区域采用精细网格）。利用这些技术，可以更精确地模拟各种复杂的结构，适应各种复杂的介质，提高了复杂介质中数值计算的精度。时域模拟一般获得的是近场电磁信息，为了得到诸如天线方向图或散射体雷达散射截面之类的远场信息，必须获得计算区域以外的频域场或
10、瞬态场。多位学者在这方面做了许多工作，发展了一种高效的时域近远场变换方法23-26。借助这种方法，可以实现由计算区域内近场数据到计算区域外远场数据的外推。目前，粗糙面散射的FDTD，传递函数在FDTD中的应用，周期介质、各向异性介质、色散介质和含有集中元件的FDTD，以及网络并行FDTD技术等方面也取得了很大进展。FDTD在迅速发展的同时，也获得了非常广泛的应用。目前，它几乎被应用到了电磁场工程中的各个方面，例如：电磁散射、生物电磁计量学、辐射天线的分析、微波器件和导行波结构的研究、散射和雷达截面的计算、周期结构的分析、电子封装和电磁兼容的分析、核电磁脉冲传播和散射的分析、以及微光学元器件中光
11、的传播和衍射特性的分析等。随着新技术的不断提出，其应用范围和成效正在迅速地扩大和提高。第二章 时域有限差分法的基本原理Maxwell方程是描述宏观电磁现象的一组基本方程。这组方程即可以写成微分形式，又可以写成积分形式。FDTD方法由Maxwell旋度方程的微分形式出发，利用二阶精度的中心差分近似，直接将微分运算转换为差分运算，这样达到了在一定体积内和一段时间上对连续电磁场数据的抽样压缩。2.1 Maxwell方程和Yee氏算法 根据27中电磁场基本方程组的微分形式，若在无源空间，其空间中的媒质是各向同性、线性和均匀的，即媒质的参数不随时间变化且各向同性，则Maxwell旋度方程可写成： （2-
12、1a） （2-1b）式中，是电场强度，单位为伏/米（V/m）；是磁场强度，单位为安/米（A/m）；表示介质介电系数，单位为法拉/米（F/m）； 表示磁导系数，单位为亨利/米（H/m）；表示介质电导率，单位为西门子/米（S/m）；表示导磁率，单位为欧姆/米（）。在直角坐标系中，（2-1）式可化为如下六个标量方程： （2-2） （2-3）这六个偏微分方程是FDTD算法的基础。 K.S.Yee3在1966年建立了如图2-1所示的空间网格，这就是著名的Yee氏元胞网格。图2-1 Yee氏网格及其电磁场分量分布并引入如下的差分近似方法对（2-2）、(2-3)式中的六个偏微分方程进行了差分离散。令代表或在
13、直角坐标系中某一分量，在时间和空间域中的离散可记为 (2-4)式中，、和分别是长方体网格沿x、y、z方向的空间步长，是时间步长，i、j、k分别是沿x、y、z方向的网格编号，n是时间步数。对关于时间和空间的一阶偏导数取中心差分近似，具有二阶精度，即 (2-5a) (2-5b) (2-5c) (2-5d) 在FDTD中，空间上连续分布的电磁场物理量离散的空间排布如图2-1所示。由图可见，电场和磁场分量在空间交叉放置，使得在每个坐标平面上每个电场分量被磁场环绕，每个磁场分量也被电场环绕。这种电磁场的空间结构与电磁感应和电磁波传播的规律相符，在每一个网格单元都能满足法拉第感应定律和安培环流定律。各分量
14、的空间相对位置也适合于Maxwell方程的差分计算，能够恰当地描述电磁场的传播特性。同时，电场和磁场在时间上交替抽样，抽样时间间隔相差半个时间步，使Maxwell旋度方程离散以后构成显式差分方程，从而可以在时间上迭代求解，而不需要进行矩阵求逆运算。因此，由给定相应电磁问题的初始条件，FDTD就可以逐步推进地求得以后各个时刻空间电磁场的分布。2.2 FDTD的基本差分方程根据上述原则，可将（2-2）、（2-3）式离散为如下的差分方程形式： （2-6a） (2-6b) (2-6c) (2-6d) (2-6e) (2-6f) 式中， (2-7a)， (2-7b)(2-6)式就是FDTD的基本差分方程
15、组。从式中可以看出，方程组中含有半个空间步和半个时间步，为了便于编程，可将（2-6）式改写成如下形式28：(2-8a)(2-8b)(2-8c)(2-8d)(2-8e)(2-8f)根据上述FDTD差分方程组可得出计算电磁场的时域推进计算方法，如图2-2所示。已知t1=t0= 0时刻空间各处的电磁场初始值计算t2=t1+ 时刻空间各处的磁场值计算t1=t2+ 时刻空间各处的电场值循环n次 图2-2 FDTD在时域的交叉半步逐步推进计算式(2-8a)(2-8c)的等号左边的电场值是第n次循环的电场值，等号右边的电场值是第n-1次循环存储在内存中的电场值，磁场值是本次循环计算得到的磁场值；式(2-8d
16、)(2-8f) 等号左边的磁场值是第n次循环的磁场值，等号右边的磁场值和电场值都是第n-1次循环存储在内存中的场值。这样，就解决了半个时间步在程序中无法表示的问题，而且也没有破坏电磁场在时间上逐步推进的逻辑关系。2.3 时域有限差分法相关技术2.3.1 数值稳定性问题上述FDTD方程是一种显式差分方程，在执行时，存在一个重要的问题：即算法的稳定性问题。这种不稳定性表现为在解显式方程时，随着时间步数的继续增加，计算结果也将无限制地增加。Taflove等4于1975年对Yee氏差分格式的稳定性进行了讨论，并导出了对时间步长的限制条件。数值解是否稳定主要取决于时间步长与空间步长、的关系。对于非均匀媒
17、质构成的计算空间选用如下的稳定性条件： (2-9)(2-9)式是空间和时间离散之间应当满足的关系，又称为Courant稳定性条件。若采用均匀立方体网格： (2-10)其中，为计算空间中的电磁波的最大速度。2.3.2 数值色散 FDTD方程组是对Maxwell旋度方程进行差分近似，在进行数值计算时，将会在计算网格中引起数字波模的色散，即在FDTD网格中，电磁波的相速与频率有关，电磁波的相速度随波长、传播方向及变量离散化的情况不同而改变。这种关系由非物理因素引起，且色散将导致非物理因素引起的脉冲波形畸变、人为的各向异性和虚假折射等现象。显然，色散与空间、时间的离散间隔有关，如下式所示： (2-11
18、)式(2-11)是三维情况下在FDTD方法中的单色平面波数值色散关系的一般形式，它表明FDTD计算中波的传播速度与传播方向有关。式中、分别是波矢量沿、方向的分量，是角频率，是被模拟的均匀介质中的光速。与数值色散关系相对应，在无耗介质中的单色平面波，色散解析关系是： (2-12)由式(2-11)可知，当式(2-11)中的、均趋于零时，它就趋于式(2-12)。也就是说数值色散是由于用近似差分替代连续微分而引起的，而且在理论上可以减小到任意程度，只要此时时间步长和空间步长都足够小，但这将大大增加所需的计算机存储空间和计算时间，并使累积误差增加。因此，在实际计算中要根据问题的性质和计算机的软硬件条件来
19、选择合适的时间步长和空间步长。为获得理想的色散关系，问题空间分割应按照小于正常网格的原则进行。一般选取的最大空间步长为，为所研究范围内电磁波的最小波长。由上分析说明，数值色散在用FDTD法分析电磁场传播中的影响是不可能避免的，但我们可以尽可能的减小数值色散的影响。2.3.3 离散网格的确定 无论是简单目标还是复杂目标，在进行FDTD离散时网格尺寸的确定，除了受计算资源的限制不可能取得很小外，还需要考虑以下几个因素：1目标离散精确度的要求。网格应当足够小以便能精确模拟目标几何形状和电磁参数。2FDTD方法本身的要求。主要是考虑色散误差的影响。设网格为立方体，所关心频段的频率上限为，对应波长为，则
20、考虑FDTD的数值色散要求 （2-13）通常。上式是根据已知所关心频率上限情况下来确定FDTD网格尺寸的；反之，若给定，则FDTD计算结果可用的上限频率也随之确定。3入射波的要求。入射波的上限截止频率应包含所关心频率范围，即。2.4 吸收边界条件 由时域有限差分法的基本原理可知，在利用时域有限差分法研究电磁场时，需在全部问题空间建立Yee氏网格空间，并存储每个单元网格上任一时间步的六个场分量用于下一时间步的计算。而在对于辐射、散射这类开放系统的实际研究中，不可能有无限大的存储空间。因此，必须在某处将网格空间截断，且在截断边界网格点处运用特殊的场分量计算方法，使得向边界面行进的波在边界处保持“外
21、向行进”特性、无明显的反射现象，并且不会使内部空间的场产生畸变，从而用有限网格空间模拟电磁波在无界空间中传播的情况。具有这种功能的边界条件称之为吸收边界条件，或辐射边界条件，或网格截断条件2931，如图2-3所示。图 2-3 附加截断边界使计算区域变为有限域从FDTD的基本差分方程组可以看出，在截断边界面上切向场分量的计算需要利用计算空间以外的电磁场分量，因此FDTD基本差分方程对这些截断边界面上的场分量失效。如何处理截断边界上的场分量，使之与需要考虑的无限空间有尽量小的差异，是FDTD中必须很好解决的一个重要问题。实际上，这是要求在误差可容忍的范围内，计算空间中的外向波能够顺利通过截断边界面
22、而不引起波的明显反射，使有限计算空间的数值模拟与实际情况趋于一致，对外向波而言，就像在无限大空间中传播一样。所以，需要一种截断边界网格处的特殊计算方法，它不仅要保证边界场计算的必要精度，而且还要大大消除非物理因素引起的波反射，使得用有限的网格空间就能模拟电磁波在无限空间中的传播。但是如果处理不当，截断边界面可能造成较大反射，构成数值模拟误差的一部分，甚至可能造成算法不稳定。加于截断边界场分量符合上述要求的算法就称为吸收边界条件（Absorbing Boundary Conditions）。2.4.1 一阶和二阶近似吸收边界条件 在截断边界附近通常没有激励源。考虑齐次波动方程 (2-14)式中，
23、表示直角坐标系下任意电磁场分量。 B.Engquist和A.Majda15利用偏微分算子对式(2-13)作因式分解，并分别取其Taylor级数展开式中的第一项和前两项近似，导出了适合直角坐标系下FDTD吸收边界条件的单向波动方程，这就是Engquist-Majda吸收边界条件。设三维长方体FDTD区域0xa，0yb，0z=np+1&i=N1-np di=0; elseif i=N1-np+1 di=np+i-N1-0.5; end %di是采样点横向距PML内边界的距离 sigma_mx=sigma_max*(di/np)M; for j=1:N1 if j=np+1&j=N1-np dj=0
24、; elseif j=N1-np+1 dj=np+j-N1-0.5; end %dj是采样点纵向距PML内边界的距离 sigma_my=sigma_max*(dj/np)M; if sigma_mx+sigma_my=0 %真空区 if j=100&i=100 t=30; %可选择高斯脉冲 term=(tt-t); pulse=exp(-4*pi*term2/202); 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(
25、i,j); Bz(i,j)=Bz(i,j)+Eterm1-Eterm2+pulse;%加入脉冲源 else 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; end else %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
26、)-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); end end end for i=2:N1 if i=np+1&i=N1-np di=0; elseif i=N1-np+1 di=np+i-N1-1; end %di是采样点横向距PML内边界的距离 sigma_
27、ex=sigma_max*(di/np)M; for j=1:N1 cam=(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); end end for i=1:N1 for j=2:N1 if j=np+1&j=N1-np dj=0; elseif j=N1-np+1 dj=np+j-N1-1; end %dj是采样点纵向距PML内边界的距离 sigma_ey=sigma_max*(
28、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); end end Bzxx(tt,1)=Bz(90,50); %对靠近边界的中央磁场点采样 Bzxx(tt,2)=Bz(50,90);3.2 出图及结果3.2.1程序部分 figure(1); %可视化处理 clf; mesh(Bz); %磁场的幅值 axis(0 N 0 N -0.5 0.5); xlabel(i) ylabel(j) drawnow;endfigure(2);plot(Bzxx);figure(3);title(磁场幅值分布图);surface(Bz);shading interp;axis square;toc3.2.2 所出的效果图正弦时谐场源辐射效果图高斯脉冲辐射效果图
