毕业设计（论文）-基于FDTD的E和H型源激励的矩形波导的截止频率和动态特性研究.doc

摘 要 时域有限差分法是由有限差分法发展出来的数值计算方法。自1966年yee在其论文中首次提出时域有限差分以来，时域有限差分法在电磁研究领域得到了广泛的应用。主要有分析辐射条线、微波器件和导行波结构的研究、散射和雷达截面计算、分析周期结构、电子封装和电磁兼容的分析、核电磁脉冲的传播和散射以及在地面的反射及对电缆传输线的干扰、微光学元器件中光的传播和衍射特性等等。由于电磁场是以场的形态存在的物质，具有独特的研究方法，采取重叠的研究方法是其重要的特点，即只有理论分析、测量、计算机模拟的结果相互佐证，才可以认为是获得了正确可信的结论。时域有限差分法就是实现直接对电磁工程问题进行计算机模拟的基本方法。在近年的研究电磁问题中，许多学者对时域脉冲源的传播和响应进行了大量的研究，主要是描述物体在瞬态电磁源作用下的理论。另外，对于物体的电特性，理论上具有几乎所有的频率成分，但实际上，只有有限的频带内的频率成分在区主要作用。 文中主要分析了源处不同性质的激励源对波导截止频率特性的影响，得到了在二维情况下对于矩形波导的场分量瞬态图以及截止频率特性图。详细模拟了te和tm模式场结构分布以及脊波导中的电磁波传播的场分量瞬态分布。关键词：时域有限差分 波导 截止频率 场结构特性目录摘 要4目录5第一章 绪论61.1课题背景与意义61.2时域有限差分法的发展与应用71.3时域有限差分法的特点91.4本课题的主要研究内容12第二章 时域有限差分法的基本原理132.1 maxwell方程和yee氏算法132.2 fdtd的基本差分方程152.3时域有限差分法相关技术172.3.1数值稳定性问题172.3.2数值色散182.3.3 fdtd离散网格的确定192.4吸收边界条件192.4.1一阶和二阶近似吸收边界条件202.4.2 二维棱边及角顶点的处理232.4.3三维棱边及角顶点的处理262.4.3 berenger完全匹配层262.5激励源设计302.5.1激励源类型302.5.2激励源的设置312.6 fdtd计算所需时间步的估计32第三章fdtd法分析矩形波导的截止频率343.1 波导理论343.2 傅立叶变换算法的特点343.3 矩形波导的截止频率35第四章波导的场结构特性434.1波导的场结构特性434.2脊波导的场结构特性48第五章结论50参考文献51攻读学位期间的研究成果55第一章 绪论1.1 课题背景与意义1873年maxwell建立电磁场基本方程以来，电磁场理论的应用和发展已经有一百多年的历史了。现代技术的许多方面都与电磁场，尤其是与高频电磁场有关，复杂的高频电磁场系统的分析和综合，以及高频电磁场与复杂目标相互作用的分析和计算，都成为现代技术发展的重要课题。应当说，通过理论分析计算，若能获得精确的解析解，对分析电磁结构的特性具有重要的指导意义。然而，只有一些典型几何形状和结构相对简单的问题才有可能求得严格的解析解。当代电磁场工程中高频电磁场问题的主要特点是电磁系统的高度复杂性。虽然对很多典型问题的解析分析仍然能帮助人们加深对电磁规律的认识，但作为工程问题的解决，解析方法往往无能为力。能够较广泛发挥作用的，唯有各种数值计算方法。数值计算方法在电磁场研究中是一种非常有效的方法，它可以结合研究对象及其周围环境的电磁参数，赋以相应的边界条件，在频域或者时域求解maxwell方程初值问题。因此，数值计算方法是理论分析计算的重要工具，可以为设计提供很有价值的参考。20世纪60年代以来，随着计算机技术的发展，一些电磁场的数值计算方法逐步发展起来，并得到广泛应用，其中主要有：属于频域技术的有限元法（fem）、矩量法(mm)和单矩法等；属于时域技术方面的时域有限差分法(fdtd)、传输线矩阵法(tlm)和时域积分方程法等。此外，还有属于高频技术的几何衍射理论(gtd)和衍射物理理论(pld)等。各种方法都具有自己的特点和局限性，在实际中经常把它们相互配合而形成各种混合方法12。其中fdtd是一种已经获得广泛应用并且有很大发展前景的时域数值计算方法。时域有限差分（fdtd）方法于1966年由k.s.yee3 提出并迅速发展，且获得广泛应用。k.s.yee用后来被称作yee氏网格的空间离散方式，把含时间变量的maxwell旋度方程转化为差分方程，并成功地模拟了电磁脉冲与理想导体作用的时域响应。但是由于当时理论的不成熟和计算机软硬件条件的限制，该方法并未得到相应的发展。20世纪80年代中期以后，随着上述两个条件限制的逐步解除，fdtd便凭借其特有的优势得以迅速发展。它能方便、精确地预测实际工程中的大量复杂电磁问题，应用范围几乎涉及所有电磁领域，成为电磁工程界和理论界研究的一个热点。目前，fdtd日趋成熟，并成为分析大部分实际电磁问题的首选方法。对于波导结构，其稳态特性的分析主要应用有限元法或矩量法等频域分析法。假定波导中传播的电磁波按简谐规律变化，则将maxwell方程中有关的运算简化为乘以的代数运算即可，由此把一个四维电磁场问题简化成三维空间的场问题。因此，这种频域分析法对分析波导结构的稳态特性是非常有效的。然而，随着微波与毫米波技术的发展，仅仅研究波导结构的稳态特性是远远不够的，还必须研究其瞬态特性。若在频域中研究波导结构的瞬态特性，由于要对每一个频率点进行单独计算，如果入射波具有较宽的频谱，那么计算工作量会很大，将四维电磁场问题化为三维场问题的频域分析法的优越性就不明显。而且，频域分析一般假设在单模传输的条件下，如果遇到波导不均匀问题，在不均匀区附近往往会激励起高次模，单模传输不能反映瞬态电磁波传输的真实过程。此外，频域分析以线性叠加原理为基础，原则上只对线性波导系统有效，对于非线性结构，频域分析法就显得十分困难。而在时域中分析波导结构的瞬态特性，尽管遇到的是四维电磁场问题，用计算机进行数值模拟需要较大的内存与较多的计算时间，但这对于计算机技术高速发展的今天来说已不是大问题。而且，时域分析只需一次计算就能得到所研究波导结构在整个频谱上的传输特性。如果入射波具有较宽的频谱，时域分析的优势就十分明显。并且，时域分析不以叠加原理为基础，处理非线性波导系统更方便，它更能反映不均匀区的瞬态电磁波的真实传输特性。另外，利用矩量法求解电磁场问题时，要用到并失green函数。对于某些问题，可以找到其解析形式的并失green函数；而对于复杂的问题，很难找到其解析形式的并失green函数，这样就使得问题无法解决。作为时域分析中的一个重要数值方法，fdtd不存在这样的问题。1.2 时域有限差分法的发展与应用经过四十多年的发展，fdtd已发展成为一种成熟的数值计算方法。在发展过程中，几乎都是围绕几个重要问题展开的，即数值稳定性、计算精度、数值色散、激励源技术以及开域电磁问题的吸收边界条件等。数值稳定和计算精度对任何一种数值计算方法都是至关重要的。a.taylor和m.e.brodwin4利用本征值方法给出了直角坐标系下fdtd的空间步长与时间步长之间的关系。x.min等5研究了存在边界条件时fdtd的稳定性问题。对于数值色散，与实际的物理色散不同，它是由电磁场量在空间和时间上的对波动方程作差分近似处理造成的。这种色散引起的误差造成在计算区域内传播的电磁波逐渐畸变67。k. l. shlager 等8比较了二维和三维空间中几种正交网格算法的色散误差。当采用其他变形或非正交网格时，必须重新分析其数值稳定性和色散特性911，p.monk 和 e.suli12分析了不均匀长方体网格算法的稳定性。激励源的设计和引入也是fdtd的一个重要任务。目前，应用最广泛的激励源引入技术是总场/散射场体系13。对于散射问题，通常在fdtd计算空间中引入连接边界，它将整个计算空间划分为内部的总场区和外部的散射场区，如图1-1。利用huygens原理，可以在连接边界处引入入射场，使入射场的加入变得简单易行。图1-1开域电磁问题中，为了在有限的计算空间内模拟无限空间中的电磁问题，必须在计算空间的截断边界处设置吸收边界条件。吸收边界条件从开始简单的插值边界，已经发展了多种吸收边界条件。在早期得到广泛应用的是g.mur14的一阶和二阶吸收边界条件，它是基于b.engquist和a.majda15的单向波方程而提出的差分格式，在fdtd仿真区域外边界具有0.5%到5%的反射系数。目前应用最广泛的是j.p.berenger16-18的分裂式完全匹配层，以及z.s.sacks等19和s.d.gedney20的各向异性介质的完全匹配层，它们可使fdtd模拟的最大动态范围达到80db。另一方面，为了更好的拟合研究对象的形状，克服台阶逼近带来的误差，d.e.merewether21提出了柱坐标系下的网格剖分方法，r.holland22提出了球坐标系下的网格剖分方法，p.monk和e.suli12提出了变网格步长方法，s.s.zivanovic等23和p.thoma等24提出了亚网格技术（即在一般区域采用粗网格，在电磁场快变区域采用精细网格）。利用这些技术，可以更精确地模拟各种复杂的结构，适应各种复杂的介质，提高了复杂介质中数值计算的精度。时域模拟一般获得的是近场电磁信息，为了得到诸如天线方向图或散射体雷达散射截面之类的远场信息，必须获得计算区域以外的频域场或瞬态场。多位学者在这方面做了许多工作，发展了一种高效的时域近远场变换方法25-28。借助这种方法，可以实现由计算区域内近场数据到计算区域外远场数据的外推。目前，粗糙面散射的fdtd，传递函数在fdtd中的应用，周期介质、各向异性介质、色散介质和含有集中元件的fdtd，以及网络并行fdtd技术等方面也取得了很大进展。fdtd在迅速发展的同时，也获得了非常广泛的应用。目前，它几乎被应用到了电磁场工程中的各个方面，例如：电磁散射、生物电磁计量学、辐射天线的分析、微波器件和导行波结构的研究、散射和雷达截面的计算、周期结构的分析、电子封装和电磁兼容的分析、核电磁脉冲传播和散射的分析、以及微光学元器件中光的传播和衍射特性的分析等。随着新技术的不断提出，其应用范围和成效正在迅速地扩大和提高。1.3 时域有限差分法的特点fdtd法对电磁场e、h分量在空间和时间上采取交替抽样的离散方式，每一个e(或h)场分量周围有四个h(或e)场分量环绕，应用这种离散方式将含时间变量的麦克斯韦旋度方程转化为一组差分方程，并在时间轴上逐步推进地求解空间电磁场。yee提出的这种抽样方式后来被称为yee元胞。fdtd方法是求解麦克斯韦微分方程的直接时域方法。在计算中将空间某一样本点的电场(或磁场)与周围格点的磁场(或电场)直接相关联，且介质参数已赋值给空间每一个元胞，因此这一方法可以处理复杂目标和非均匀介质物体的电磁场问题。作为一种电磁场的数值计算方法，时域有限差分法具有一些非常突出的特点，最主要包括以下几点：首先，fdtd是一种求解maxwell方程的直接时域方法。在直角坐标系中，时域有限差分法直接把含时间变量的麦克斯韦旋度方程在yee时网格空间中转换为差分方程。在这种差分格式中，每个网格点上的电场（磁场）分量仅与它相邻的磁场（或点场）分量及上一时间步该点的场值有关。在每一时间步计算网格空间格点的电场和磁场分量，随着时间步的推进即能直接模拟电磁波的传播及其与物体的相互作用过程。时域有限差分法把各类问题都作为初值问题来处理，使电磁波的时域特性被直接反映出来。这一特点使它能直接给出非常丰富的电磁场问题的时域信息，为复杂的物理过程描绘出清晰的物理图像。如果需要频域信息，则只需对时域信息进行傅里叶变换。为获得宽频带的信息，只需在宽频谱的脉冲激励下进行一次计算。这样不仅保证了faraday电磁感应定律和ampere环路积分定律，以及介质分界面上切向场分量的连续性条件得到自然满足，而且可以方便地将maxwell旋度方程在空间和时间上进行二阶精度的中心差分运算，从而转化成一组有限差分方程，然后就可以在相应的边界条件和初始条件下求解这组方程。于是，随着时间的逐步推进，便能直接模拟电磁波的传播及其与物体相互作用的过程，从而可以得到整个计算空间随时间变化的电磁场。由此可见，fdtd可以很恰当地模拟电磁波的实际传播过程，直接给出非常丰富的电磁场问题的时域信息，给复杂的物理过程描绘出清晰的物理图像。如果需要频域信息，只需在宽频谱的脉冲激励下进行一次计算，然后对获得的时域信息做fourier变换即可。其次，时域有限差分法具有广泛的适用性。由于时域有限差分法的直接出发点时概括了电磁场普遍规律的麦克斯韦方程，这就预示着该方法应具有最广泛的适用性。近几年的发展完全证实了这点。从具体的算法看，在时域有限差分法的差分格式中，被模拟空间电磁性质的参量是按空间网格给出的，因此，只需为相应空间点设定适当的参数，就可模拟各种复杂的电磁结构。媒质的非均匀性、各向异性、色散特性和非线性等均能很容易得进行精确模拟。由于在网格空间中电场和磁场各分量是被交叉放置的，而且计算中用差分代替了微商，使得介质交接面上的边界条件能自然得到满足，这就为模拟复杂的结构提供了极大的方便。对于散射、辐射、传播、透入、或吸收问题，只要能对波源和被研究物体的结构进行正确模拟，时域有限差分法就能给出正确的答案。第三，时域有限差分法能够节约存储空间和计算空间，并且十分适合并行计算。在时域有限差分中，每个网格电场和磁场的六个分量及其上一时间步的值是必须存储的，此外还有描述各网格电磁性质的参数以及吸收边界体条件和连接边界条件的有关参量，它们一般是空间网格总数n的数倍。所以，时域有限差分法所需要的存储空间直接由设定的网格空间决定，与网格总数n成正比。在计算时，每个网格的电磁场都按同样的差分格式计算，所以，所需的主要计算时间也与网格总数n成正比。相比之下，若离散单元也是n，则矩量法所需的存储空间与成正比，而所需的cpu时间则与至成正比。当n较大时，二者的差别是很明显的。所以，当n很大时，时域有限差分法往往是更合适的方法。很多复杂的电磁场问题不能计算往往不是没有可选用的方法，而是由于计算条件的限制。当代计算机发展方向是利用并行处理技术进一步提高计算速度。并行计算机的发展推动了数值计算中并行处理的研究，适合并行计算的方法将更多的发挥作用。如前面所指出的时域有限差分法的计算特点是，每一网格上的电场（或磁场）值与其周围相邻网格点处的磁场（或电场）及其上一时间步的场值有关，这使得它特别适合并行计算。实行并行计算可使时域有限差分法所需的存储空间和计算时间减少为只与成正比。以直角坐标系中的立方体网格空间为例，若每个坐标方向的网格数为n，这计算网格空间的网格总数。若用个处理器，则每一处理器只需记忆和处理一行中一个场分量的有关信息，行可同时处理。这样，对于一个确定的时间步，全部运行时间就正比于完成一行处理所需的时间，这时间又正比于一行中一个场分量的个数n，即。由此看来，时域有限差分法将随着并行计算机的发展而越来越显得重要。第四，时域有限差分法的计算程序具有很好的通用性。由于麦克斯韦方程是时域有限差分法计算任何问题的数学模型，因而它的基本差分方程对广泛的问题是不变的。此外，吸收边界条件和连接条件对很多问题是可以通用的，而计算对象的模拟是通过给网格赋予参数来实现，和以上各部分没有直接联系，可以独立进行。因此一个基础的时域有限差分法计算程序，对广泛的电磁场问题具有通用性，对不同的问题或不同的计算对象只需修改有关部分，而大部分是共同的。第五，简单、直观、容易掌握。由于时域有限差分直接从maxwell方程出发，不需要任何导出方程，这样就避免了使用更多的数学工具，使得它成为所有电磁场的计算方法中最简单的一种。其次，由于它能直接在时域中模拟电磁波的传播及其与物体作用的物理过程，所以它又是非常直观的一种方法。这样，这一方法很容易得到推广，并在广泛的领域发挥作用。1.4 本课题的主要研究内容根据以上介绍，确定了本文的研究方向，即时域有限差分研究波导中电磁波的传播特性，以验证fdtd中对于在源处加入不同性质的激励源对分析波导传输特性的有效性；在此基础上，分析了波导的场结构分布特性。本文共分五章，各章的安排如下：第一章阐述了数值计算法在电磁理论中的重要性，然后介绍了时域有限差分法的发展与应用，最后介绍了时域有限差分法的特点，明确了研究方向。第二章结合波导的特点，介绍了时域有限差分法的基本原理，主要包括：（1）yee氏网格和fdtd基本差分方程的建立，以及fdtd的数值稳定性和数值色散；（2）几种适合波导结构计算的激励源，以及波导结构中的总场/散射场体系；（3）g. mur的一阶和二阶边界条件，以及j. p. berenger的分裂式完全匹配层第三章将fdtd应用到波导结构的具体算例中分析波导的截止频率特性，验证了用matlab所编写的fdtd程序的正确性。第四章用fdtd法分析波导的传输特性。（1）介绍了傅立叶变换的基本理论，给出了时域与频域变换的基本形式。（2）对波导激励源处加入性质不同的源即电场源和磁场源，得到了相关的截止特性，并对结果作了分析和比较。第五章对本文的研究工作做了全面的总结。第二章 时域有限差分法的基本原理maxwell方程是描述宏观电磁现象的一组基本方程。这组方程即可以写成微分形式，又可以写成积分形式。fdtd方法由maxwell旋度方程的微分形式出发，利用二阶精度的中心差分近似，直接将微分运算转换为差分运算，这样达到了在一定体积内和一段时间上对连续电磁场数据的抽样压缩。2.1 maxwell方程和yee氏算法 根据29中电磁场基本方程组的微分形式，若在无源空间，其空间中的媒质是各向同性、线性和均匀的，即媒质的参数不随时间变化且各向同性，则maxwell旋度方程可写成： （2-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)式中的六个偏微分方程进行了差分离散。令代表或在直角坐标系中某一分量，在时间和空间域中的离散可记为 (2-4)式中，、和分别是长方体网格沿x、y、z方向的空间步长，是时间步长，i、j、k分别是沿x、y、z方向的网格编号，n是时间步数。对关于时间和空间的一阶偏导数取中心差分近似，具有二阶精度，即 (2-5a) (2-5b) (2-5c) (2-5d) 在fdtd中，空间上连续分布的电磁场物理量离散的空间排布如图2-1所示。由图可见，电场和磁场分量在空间交叉放置，使得在每个坐标平面上每个电场分量被磁场环绕，每个磁场分量也被电场环绕。这种电磁场的空间结构与电磁感应和电磁波传播的规律相符，在每一个网格单元都能满足法拉第感应定律和安培环流定律。各分量的空间相对位置也适合于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的基本差分方程组。从式中可以看出，方程组中含有半个空间步和半个时间步，为了便于编程，可将（2-6）式改写成如下形式30：(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)(2-8f) 等号左边的磁场值是第n次循环的磁场值，等号右边的磁场值和电场值都是第n-1次循环存储在内存中的场值。这样，就解决了半个时间步在程序中无法表示的问题，而且也没有破坏电磁场在时间上逐步推进的逻辑关系。2.3时域有限差分法相关技术2.3.1数值稳定性问题上述fdtd方程是一种显式差分方程，在执行时，存在一个重要的问题：即算法的稳定性问题。这种不稳定性表现为在解显式方程时，随着时间步数的继续增加，计算结果也将无限制地增加。taflove等4于1975年对yee氏差分格式的稳定性进行了讨论，并导出了对时间步长的限制条件。数值解是否稳定主要取决于时间步长与空间步长、的关系。对于非均匀媒质构成的计算空间选用如下的稳定性条件： (2-9)(2-9)式是空间和时间离散之间应当满足的关系，又称为courant稳定性条件。若采用均匀立方体网格： (2-10)其中，为计算空间中的电磁波的最大速度。2.3.2数值色散 fdtd方程组是对maxwell旋度方程进行差分近似，在进行数值计算时，将会在计算网格中引起数字波模的色散，即在fdtd网格中，电磁波的相速与频率有关，电磁波的相速度随波长、传播方向及变量离散化的情况不同而改变。这种关系由非物理因素引起，且色散将导致非物理因素引起的脉冲波形畸变、人为的各向异性和虚假折射等现象。显然，色散与空间、时间的离散间隔有关，如下式所示： (2-11)式(2-11)是三维情况下在fdtd方法中的单色平面波数值色散关系的一般形式，它表明fdtd计算中波的传播速度与传播方向有关。式中、分别是波矢量沿、方向的分量，是角频率，是被模拟的均匀介质中的光速。与数值色散关系相对应，在无耗介质中的单色平面波，色散解析关系是： (2-12)由式(2-11)可知，当式(2-11)中的、均趋于零时，它就趋于式(2-12)。也就是说数值色散是由于用近似差分替代连续微分而引起的，而且在理论上可以减小到任意程度，只要此时时间步长和空间步长都足够小，但这将大大增加所需的计算机存储空间和计算时间，并使累积误差增加。因此，在实际计算中要根据问题的性质和计算机的软硬件条件来选择合适的时间步长和空间步长。为获得理想的色散关系，问题空间分割应按照小于正常网格的原则进行。一般选取的最大空间步长为，为所研究范围内电磁波的最小波长。由上分析说明，数值色散在用fdtd法分析电磁场传播中的影响是不可能避免的，但我们可以尽可能的减小数值色散的影响。2.3.3 fdtd离散网格的确定 无论是简单目标还是复杂目标，在进行fdtd离散时网格尺寸的确定，除了受计算资源的限制不可能取得很小外，还需要考虑以下几个因素：1目标离散精确度的要求。网格应当足够小以便能精确模拟目标几何形状和电磁参数。2fdtd方法本身的要求。主要是考虑色散误差的影响。设网格为立方体，所关心频段的频率上限为，对应波长为，则考虑fdtd的数值色散要求 （2-13）通常。上式是根据已知所关心频率上限情况下来确定fdtd网格尺寸的；反之，若给定，则fdtd计算结果可用的上限频率也随之确定。 3入射波的要求。入射波的上限截止频率应包含所关心频率范围，即。2.4吸收边界条件 由时域有限差分法的基本原理可知，在利用时域有限差分法研究电磁场时，需在全部问题空间建立yee氏网格空间，并存储每个单元网格上任一时间步的六个场分量用于下一时间步的计算。而在对于辐射、散射这类开放系统的实际研究中，不可能有无限大的存储空间。因此，必须在某处将网格空间截断，且在截断边界网格点处运用特殊的场分量计算方法，使得向边界面行进的波在边界处保持“外向行进”特性、无明显的反射现象，并且不会使内部空间的场产生畸变，从而用有限网格空间模拟电磁波在无界空间中传播的情况。具有这种功能的边界条件称之为吸收边界条件，或辐射边界条件，或网格截断条件3133，如图2-3所示。图 2-3 附加截断边界使计算区域变为有限域从fdtd的基本差分方程组可以看出，在截断边界面上切向场分量的计算需要利用计算空间以外的电磁场分量，因此fdtd基本差分方程对这些截断边界面上的场分量失效。如何处理截断边界上的场分量，使之与需要考虑的无限空间有尽量小的差异，是fdtd中必须很好解决的一个重要问题。实际上，这是要求在误差可容忍的范围内，计算空间中的外向波能够顺利通过截断边界面而不引起波的明显反射，使有限计算空间的数值模拟与实际情况趋于一致，对外向波而言，就像在无限大空间中传播一样。所以，需要一种截断边界网格处的特殊计算方法，它不仅要保证边界场计算的必要精度，而且还要大大消除非物理因素引起的波反射，使得用有限的网格空间就能模拟电磁波在无限空间中的传播。但是如果处理不当，截断边界面可能造成较大反射，构成数值模拟误差的一部分，甚至可能造成算法不稳定。加于截断边界场分量符合上述要求的算法就称为吸收边界条件（absorbing boundary conditions）。2.4.1一阶和二阶近似吸收边界条件 在截断边界附近通常没有激励源。考虑齐次波动方程 (2-14)式中，表示直角坐标系下任意电磁场分量。 b.engquist和a.majda15利用偏微分算子对式(2-13)作因式分解，并分别取其taylor级数展开式中的第一项和前两项近似，导出了适合直角坐标系下fdtd吸收边界条件的单向波动方程，这就是engquist-majda吸收边界条件。设三维长方体fdtd区域0 xa，0yb，0zd，这时有六个截断边界，其一阶和二阶解析吸收边界条件的具体形式见表2-1，其中代表任一直角场分量。g. mur14对表中的吸收边界条件引入了一种简单有效的差分数值算法，即对时间和空间的偏微分取二阶中心差分近似，将单向波方程离散化，便形成了著名的g.mur的一阶和二阶吸收边界条件，其总体虚假反射在1%5%。图2-5给出了反射系数（反射波与入射波）与入射角的关系，由图2-4可看出，仅当入射角较小时其反射系数较小。表2-1 三维长方体fdtd区域的一阶和二阶吸收边界条件一阶近似二阶近似图2-4 近似吸收边界条件作用后残留的反射波与入射波之比根据yee氏元胞网格的特点，在fdtd截断边界面上只有电场切向分量和磁场法向分量。以界面为例，此界面仅有，节点。由于fdtd中的计算式不涉及区域，即不涉及截断边界界面外的节点。所以，吸收边界条件将不考虑，而只考虑电场切向分量和。以为例，对一阶和二阶吸收边界条件分别有 (2-15a) (2-15b)将式(2-15b)分别在距离边界半个空间步长的辅助网格点处及时刻离散，此时的位置在，并对各项做差分近似，同时对辅助网格点处的场值应用线性插值，假设，可得到一阶吸收边界条件在三维情况的形式： (2-16)式中代表截断边界上切向场分量。将式(2-16)所示二阶吸收边界推广到长方体元胞各边、和不相等情形。设边界为，对于分量有 (2-17)同理，可以得到所有截断边界面上切向场分量的一阶和二阶差分式。需要注意的是，用二阶近似条件计算界面上与棱边相邻的一列节点时会涉及棱边上的场值。因此，要想避免用到棱边上的场值，只需对截断边界面上切向场分量的计算按以下两种情况区别对待：（1）截断边界面上与棱边相邻的一列场分量采用一阶差分式；（2）截断边界面上其它场分量采用二阶差分式。这样，就完全不必考虑棱边上的场分量，避免了计算棱边上场分量所带来的误差。实际计算表明，这样做提高了吸收边界条件的精度和fdtd计算的稳定性。但是，就目前fdtd的发展来看，g.mur的一阶和二阶吸收边界条件已不能满足目前高精度计算的要求。2.4.2 二维棱边及角顶点的处理对于二维电磁场问题，在用fdtd计算边界处的元胞时，将涉及到截断边界外侧的节点。如对于二维tm情况，电磁场分量有，由图2-5可见，在用fdtd计算边界处的tm元胞和时并不涉及截断边界以外或的节点。只有涉及截断边界外侧的节点。因此，只需给出边界处切向场分量的吸收边界条件。同样，对于te波只需给出边界处切向场分量的吸收边界条件。表2-2为矩形截断边界面上节点的二阶mur吸收边界条件。图2-5 二维tm左截断边界元胞表2-2 矩形截断边界四边上的二阶mur吸收边界条件截断边界位置二阶mur吸收边界条件在二维矩形计算区域的角点，吸收边界条件的离散式需特殊考虑。假设角点附近只有向外传播的行波，且传播方向沿角点处元胞的对角线，如图2-6所示的矩形区域左下角点处的tm元胞为例，导出适用于角点的吸收边界条件离散式。根据courant稳定条件式，有。在点与对角点之间取一点，对于沿对角线的外行波，有以下等式： （2-18）由于传播距离很小，上式中略去了振幅的衰减。在利用线性插值 （2-19）由式（2-18）解出后代入（2-17）式得 （2-20）或（2-21）对于矩形域的其它角点可作类似的处理。对于二维te波情况，将式（2-19），（2-21）式中的换为即可。 图2-6 矩形域四个角点2.4.3三维棱边及角顶点的处理 在三维fdtd长方体计算区域有6个截断边界面和12条棱边，如图2-7所示。与x轴平行的棱边上仅有场分量，与y轴平行的棱边上仅有场分量，与z轴平行的棱边上仅有场分量。图2-7 三维fdtd长方体计算区域实际上，可以将二维角点的吸收边界条件式（2-21）应用于三维棱边。2.4.3 berenger完全匹配层完全匹配层1618（perfectly matched layer, pml）是1994年由j.p.berenger首先提出，并将其设置在fdtd计算区域截断边界处，用来吸收外向电磁波。berenger假设将电磁场分量在pml介质中分裂，并分别对各个分裂的场分量赋以不同的损耗。这就相当于在fdtd区域截断边界外设置了一种特殊的非物理的吸收介质层，该层介质的波阻抗与相邻介质的波阻抗完全匹配，因而外向波将无反射地穿过分界面进入pml。同时，由于pml为有耗介质，而且不依赖于外向波的入射角和波阻抗，即使为有限厚度，外向波在其中也会迅速衰减。在实际计算中，pml是目前一种很常用的吸收边界条件，有很好的吸收效果，其总的网格噪声能量是使用普通吸收边界条件时的1/107，可使fdtd模拟的最大动态范围达到80db。在pml介质中，其网格剖分方式与常规fdtd网格完全一致，每个场分量在网格中的位置也不变，只是都被分裂为两个子分量。这样，式(2-1)可写成 (2-22a) (2-22b) (2-22c) (2-22d) (2-22e) (2-22f) (2-22g) (2-22h) (2-22i) (2-22j) (2-17k) (2-22l)式中，为电导率和磁导率，描述了pml介质的各向异性。当且时， pml介质退化为普通有耗介质;当且时，退化为自由空间中的maxwell方程。同时，要满足pml介质的重要基本条件即阻抗匹配条件，如下式所示： (2-23)式中和分别为真空的介电常数和真空的磁导率。对式(2-18)做差分处理，可得到pml介质中的fdtd差分方程表达式：(2-23a)(2-23b)(2-23c)(2-23d)(2-23e) (2-23f) (2-23g) (2-23h) (2-23i) (2-23j) (2-23k) (2-23l)式中， (2-24)， (2-25)其中，； 在实际计算中，在fdtd计算中pml的设置不可能延伸到半无限空间，只能是有限厚度，因此pml的外侧边界需要特殊处理，通常情况下采用理想导体边界截断。这样，透入pml中的外向波到达理想导体边界处会反射回来，重新进入计算区域，pml的反射系数不再等于零。pml介质内沿方向的电导率分布通常采取以下函数形式 （2-26）式中，为pml层的厚度，为pml层靠近fdtd分界面的距离，是电导率分布阶数，表示pml中电导率变化的程度，取整数。式（2-26）说明了实际计算中pml介质层厚度必须取若干个空间步长，使电导率从fdtd-pml分界面的0渐变到pml最外面的，避免电导率跃变太大，尽量消除数值反射。这样，如果相对于fdtd-pml分界面定义的外向波入射角为，则pml内侧表面反射系数为 （2-27）使用pml吸收边界条件时，首先要选定三个参数：pml层数，电导率分布阶数，pml表面反射系数。大量的数值实验表明：（1）当层数n固定时，减小，即增加pml的衰减，可以使局部及总体误差都单调地减小。然而，当减小到一定程度后，这种现象不再出现，原因是存在由空间网格引起的固有误差。（2）增加pml层数，可以使局部及总体误差都单调地减小，但是n过大又会使计算量剧增，需折衷考虑吸收效果和计算量。（3）电导率分布阶数的改变不会影响计算量，却会影响pml的吸收效果，一般情况下，越大，吸收效果越好。因此，在实际计算中参数的选择随所处理问题的不同而不同，需综合考虑，并由数值实验寻找最佳值。2.5激励源设计实际的电磁场问题总是要在合适的激励源激励下求解，因此，恰当地将激励源引入到fdtd网格之中，对于正确地模拟电磁场问题是至关重要的。在激励源的引入过程中，为了尽量减小计算机内存占用和计算时间、提高整个程序的效率，通常要求激励源的实现尽可能地紧凑，即在fdtd网格中只用很少的几个电（磁）场分量就可实现对激励源的恰当模拟。就是说，应将被研究媒质在真实源激励下这一完整条件在数值计算中尽可能近似的“复现”出来，或者说对所研究的主要参数来讲，二者是等效的34。2.5.1激励源类型在关于时域有限差分求解的问题中，可以设置各种类型的源，这些源可以按不同的形式和性质来分类。从空间分布来看，有面源、线源、点源等。典型的面源有常见的平面波源。在研究柱面波的问题时，常设置线源。从频谱特性来看，有工作于一个频率上的周期变化的连续波源，亦可以是占有较宽频谱的波源。从源的时变特点看，主要有二大类：一类是随时间周期变化的时谐源；另一类是对时间呈冲击函数形式的波源。后者时变源包括矩形脉冲、高斯脉冲、上升余弦脉冲和核电脉冲等形式。从源处场的性质看，有可分为型和型源。即在源面上仅赋值与电场量或磁场量，一般常用的是型源。应当指出，源的时变特性和其频谱特性是相关的，二者由傅立叶变换相联系的，但源的时域频域特性是与其空间分布状态不相关的。任一种空间分布均可采用不同时变特性的源，例如平面型波源，其时间变化规律既可以是时谐型的，亦可以是任一冲击函数型的。下面简要介绍几种时变源。1. 时谐型 （2-28）式中是频率，是时间步长，为整数，代表递增的时间步数。2. 高斯脉冲型 （2-29）式中为常数，决定了高斯脉冲的宽度。脉冲峰值出现在时刻，应合适的选择 与，以保证源的初始条件，即。3. 升余弦脉冲 （2-30）为脉冲底座宽度。4. 微分高斯脉冲 （2-31）5. 调制高斯脉冲型 （2-32）上式右边第一项为基波表达式，中心频率为；第二项为高斯函数形式，通常取基波的个周期，即。2.5.2激励源的设置fdtd仿真中，激励源的正确设置是一个非常重要的部分，它必须反映场的分布，而且能够覆盖所关心的频段。就散射问题而言35，我们所关心的主要问题是散射场的性质。但在源场存在的空间中，得到的是总场(入射场+散射场)，而不是散射场。另一方面不论是散射场还是入射场还是总场均应满足麦克斯韦方程式，也就是说，fdtd算法对任一种的求解均是适用的。为了在fdtd运算过程中“提取”出散射场，如何设置入射波源是一个值得研究的问题。一种解决的办法是将fdtd网格空间划分成总场和散射场区。如图2-8所示，将fdtd的计算空间分成二个区域，其分界面为假想的连接面。 图 2-8 总场区和散射场区的划分区域1称为fdtd网格空间的内区，标作总场区，在此区内“装入”所研究的散射体。在总场区内存在模拟的平面波源，该平面波照射到散射体上产生散射场。在该区域内，麦克斯韦方程作用于总场矢的各分量。区域2为网格空间的外区，标为散射场区。该区中不存在入射波。麦克斯韦方程仅作用于散射场矢的各分量。该区的外表面周界为吸收边界，它能近似无反射的吸收外向行进的散射波。它的存在模拟将场的空间扩展到无限远处。连接边界为网格空间的内区和外区的公共界面，在此界面设置平面波源以在区域1内产生平面波。该平面波具有预定的时间波形、入射角和极化角该连接边界的存在几乎完全将入射波“禁锢”在区域1内，并且能有效地将外向散射波转移进入外区2内。除了平面波源外，还有其他类型的激励源。这些初始源的设置大多模拟实际源存在条件，或是激发出所需的场分布。2.6 fdtd计算所需时间步的估计 为了使计算达到稳定,通常计算所需要时间步按照电磁波往返穿越fdtd计算区对角线35次来估计。若fdtd计算区总元胞数为，则对角线上元胞约为 （三维)。按照courant稳定条件，设计算中心区，即穿越对角线一次需要时间步为。总计算时间步约需步。对于二维情况则约为。或者说，计算时间步大约等于fdtd计算区对角线上元胞数目的1220倍。实际上，计算所需时间步还与散射体具体形状、结构有关。图2-9给出了应用fdtd分析电磁场问题时的程序流程图图2-9 fdtd 程序流程图第三章 fdtd法分析矩形波导的截止频率本章首先介绍了波导的基本理论，然后介绍了傅立叶变换的特点，最后将fdtd应用到具体的波导结构中，求解矩形波导的截止频率，验证了用matlab所编写的fdtd程序的正确性，为场结构图的分析奠定基础。3.1 波导理论波导是用来引导电磁波在有限空间中传播，使电磁波不至于扩散到漫无边际的空间中去的结构的总称。在各类电子设备中，电磁能量的传递是通过传输线进行的。电磁能量是以电磁波的形式在传输线中传播的。也就是说，传输线是用来引导电磁波做定向传播的一种导波结构，所以，传输线又可以成为波导36 。在传输线引导下定向传播的电磁波被称为导行电磁波，简称为导波。研究导行电磁波在波导中的传播，就是求解电磁波在波导横截面上的分布规律及其沿轴向的传播特性。波导的种类很多，在工程中常见的有平行双线传输线、同轴线、圆波导、矩形波导、微带、介质波导、脊波导等等，如图所示。 图3-1a 平行双线传输线 图3-1b 同轴线 图3-1c 圆波导 图3-1d 矩形波导 图3-1e 微带线3.2 傅立叶变换算法的特点快速fourier变换算法3738 (fast fouriertransform algorithm,fft)使许多科学领域在内插三角多项式的使用上有了革命性的进展。此方法须将数据点的数目能够加以分解成为某一乘幂的形式,通常是以2的乘幂为数目,即逐次分半。由于fft算法的核心思想就是尽量减少乘法次数，在逐次分半法中恒取，并把一次点变换的过程分解为步。每步计算一个简化的傅氏变换，工作量为n个复运算，因而总工作量为。由于分点数n的选取在多数场合下人们有主动权，因此取这一特殊形式并不是严重的限制。现用公式计算全部，其中（正变换），是已知复数数列。表面上看有 个的项，而这个项中实际上只有个不同的值，把中各项先按归类，然后把同类项中先加起来再和相乘，这样就能大量地减少乘法运算的次数。3.3 矩形波导的截止频率 在波导和传输线本征值问题中，导波结构沿纵向是不变的，波导中的电磁场沿z方向的变化可忽略不计，即，则三维问题便转化为二维问题。因此，根据式(2-6)可得横磁模式（tm模式）、横电模式（te模式）的电场和磁场纵向分量的差分表达式分别为te模：(3-1a) (3-1b) (3-1c)tm模： (3-2a) (3-2b) (3-2c)其中， ， , (3-3)以二维波导3940的设计条件与傅立叶变换的特点，采用fdtd法，场分量可离散为 （3-4）式中，n为时间步数。为所研究模式的空间分布函数，为其截止频率，是一个复常数。对（3-4）式进行离散傅立叶变换，可以得到频率的离散函数 (3-5)将式（3-4）代入式（3-5），由离散傅立叶变换和连续傅立叶变换的关系，可得 （3-6）从上式可看出,取极大值时所对应的频率值近似等于模式s的截止频率,因此只要由式(3-6)求出的各极大值,通过时域迭代后,在各极大值点上进行傅立叶变换即可获得截止频率和相应的场分布。在波导内任意一点场量做傅立叶变换。以高斯脉冲(2-25)为例，得高斯脉冲的频域表达式为 式中为高斯脉冲的频宽，其它参数同(2-25)式中的意义一样。根据29中，对于矩形波导的截止频率可根据式（3-7）求得， （3-7）其中，a,b为波导尺寸，m,n取不同时为零整数。相应的截止波长为 （3-8）由式（3-8）可知，如果波导的尺寸决定，那么，m,n的值越小，截止波长越长。模和模分别是te波和tm波中有最长截止波长的模式，即最低模式，而模的截止波长又比模的截止波长长，也就是说在矩形波导中具有最长的截止波长，通常将模作为主模或基模。由（3-7）、（3-8）可见截止频率和波长都与工作频率无关，仅与波导的尺寸和模式相关。为了比较不同类型源的激励特性，在波导中心线上以，极化后的线源激励的波导截止频率计算结果如图3-2所示。由图是可知，计划线源激励的、模的模场分布均在中心线上分量的最大值。 图3-2 线源ey激励的te模以的矩形波导为例，在归一化频率为05范围内进行计算，为主模的截止频率。根据解析解可知在这个频率范围内的最高模为，根据29可知模与模简并。在归一化频率范围内共有14个模。7个模。网格为6030，源参数取t=15dt(s)，t0=140dt(s)，时间步长取dt=2.710-13s。取点源位置m+n-1=4，在20000个时间步后对波导截面上一点(也应避开波节点)对hz和ez分别作dft，如图3-4中的谐振峰对应了所有导模的截止频率。表3-1和3-2对低阶te模和tm模截止频率的理论值、参考文献和计算值进行了比较。（a）te模 （b）tm 模 图3-3文献37的点源激励的模式结果 （a）te模 （b）tm 模 图3-4 本文的点源激励结果 表3-1 矩形波导te模截止频率（ghz）模式te10te20 te01te11te21te30te31te40 te02te12te41 te22te32 te 50解析解21.08242.16347.14059.62863.24576.01084.32686.92194.279105.41文献121.07642.15247.14059.64563.15876.01486.23486.83394.210105.38本文21.0942.3546.9359.5262.9576.1284.7086.9994.43105.3表3-2矩形波导tm模截止频率（ghz）模式tm11tm21tm31tm12tm22 tm41tm32解析解47.14059.62876.01086.92194.279105.410文献147.14058.94376.01486.83394.842105.380本文46.81259.44176.27886.79894.165105.203根据fdtd法在分析电磁场问题时，激励源的模拟时是常重要的。根据fdtd的基本原理与激励源