第一部分 有限差分时域法

1、电磁场计算机辅助分析,北京交通大学 电子信息工程学院,课程介绍,(1-1),电磁波在实际环境中的传播过程十分复杂。大多数电磁波传播问题 (如复杂目标散射，复杂结构天线的辐射，复杂地形及海面对电磁波传播的影响)，由于边界条件的影响，很难通过偏微分或积分形式的麦克斯韦方程组得到电磁波的解析解。因此发展了一系列的数值方法来得到近似解。目前主要的数值解方法有矩量法(MoM)、有限元法(FEM)、边界元法(BEM)以及时域有限差分法(FDTD)等。,本门课程主要介绍矩量法(MoM)、有限元法(FEM)、以及时域有限差分法(FDTD)的算法原理，以及电磁仿真软件HFSS（基于FEM ）、FEKO（基于Mo

2、M）的基本应用。,(1-1),要求同学们：,理解时域有限差分法(FDTD)算法原理，基于该算法利用计算机编程实现对简单电磁问题的求解; 理解有限元法(FEM)算法原理，能够利用该算法对简单电磁问题进行求解; 理解矩量法(MoM)算法原理， 能够利用该算法对简单电磁问题进行求解; 掌握HFSS仿真软件的应用方法，能够利用该软件实现对简单天线的仿真分析; 掌握FEKO/NEC仿真软件的应用方法，能够利用该软件实现对简单天线的仿真分析。,主要参考文献: 1Kane S.Yee, Numerical solution of initial boundary value problems involvi

3、ng Maxwells equations in isotropic media, IEEE Trans. on Antennas and Propagation Vol. AP-14, No.3, 302-307 (1966). 2王长清, 祝西里. 电磁场计算中的时域有限差分法. 北京: 北京大学出版社.1994. 3葛德彪, 闫玉波. 电磁波时域有限差分方法.西安电子科技大学出版社,2011. 4金建铭（美）.王建国译.电磁场有限元方法.西安电子科技大学出版社.1998.,5哈林登(美).计算电磁场的矩量法。国防工业出版社. 1981. 6何国瑜,卢才成,洪家才,邓晖.电磁散射的计算和测

4、量.北京航空航天大学出版社.2006. 7A.Taflove, M.E.Brodwin, Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Maxwells equations, IEEE Trans. on Microwave Theory Tech., MTT-23, pp623-630, Aug. 1975. 8R.C. Booton, Computational Methods for Electromagnetics and Microwave

5、s, Wiley, 1992.,答疑地点：9教学楼南318；时间：每周五下午,联系方式：;考核方式：平时考勤、作业、期末考试成绩,期末成绩计算方法：平时考勤10%+作业20%+期末考试成绩70%,期末考试时间：暂定第16周周三第二节（即本门课最后一堂课）。,数值方法的发展历史,(1-1),1943年柯朗（Courant）首次提出有限元法(FEM),用于近似求解数理边值问题。20世纪50年代应用于飞机设计，在六七十年代被引入电磁场问题的求解。 结构力学仿真软件ANSYS、电磁仿真软件HFSS的计算方法采用有限元法.,1966年K.S.Yee提出时域有限差分法(FDTD)

6、，采用差分格式模拟电磁场在时间和空间演变过程。 XFDTD、EMPIRE、Fidelity等电磁仿真软件的计算方法采用有限差分时域法.,1968年Harrington提出距量法(MoM)，用于近似求解电磁场的边值问题。 NEC、FEKO、Microwave Office等电磁仿真软件的计算方法采用矩量法.,电磁场数值求解的基本步骤,1、离散化(discretization): 将所研究空间通过合适的方法离散为若干顶点（vertex）,将时间离散为若干个时间点。两个邻近顶点间的距离称为网格尺寸（meshsize）,两个相邻时间点的间隔称为时间步进(timestep)。 2、建立代数方程组：将偏微

7、分或积分方程近似或离散化为代数方程组。代数方程组中的未知数即为离散点的解。 3、计算：利用计算机通过直接或迭代技术解代数方程组。注意，所得代数方程组的确切解只是原始偏微分方程的近似解，因为离散化过程会引入离散误差。,数值求解的有效性,(1-1),只有当离散后代数方程组的解是稳定(stable)及收敛(convergent)，才可以说数值求解是有效的，即可以用代数方程组的解来代替原来电磁场偏微分方程组的解。,1、稳定性是指在数值求解过程中，代数方程组的解始终是有界的。 2、收敛性是指当离散间隔趋于零时，代数方程组的解在空间任意一点和任意时刻都一致趋于原方程的解。通常以计算结果的误差是否小于预设值

8、来判断是否收敛。 *此处的收敛性与HFSS中自适应迭代所要求的收敛不同。详见第四部分HFSS简介。,数值求解引起的误差,(1-1),1、截断误差（truncation error）：在数值分析过程中，将无穷序列的和由有限求和代替所引入的误差。如一个泰勒级数的n阶近似。 2、离散误差(discretization error)：将连续空间及时间离散化所引入的误差。 3、舍入误差(round-off error)：在计算过程中，由于有限的数据位数（数值精度）所引入的误差。,思考：随着网格尺寸和时间步进的减小，舍入误差和离散误差的变化趋势？,随着离散间隔的减小，舍入误差会逐步变大，离散误差则会随之减

9、少。,电磁场计算机辅助分析,第一部分 有限差分时域法,第一章 引言,目录,第二章 FDTD算法介绍,第三章 吸收边界条件,第四章 激励源,有限差分时域法（Finite-Difference Time-Domain Method, FDTD）是计算电磁学仿真中所用的一种数值技术。最早由K. S. Yee 于1966年在其论文Numerical solution of initial boundary value problems involving Maxwells equations in isotropic mediaIEEE Trans. Antennas Propagat. Page(s

10、): 302-307, 1966, Volume: AP-14 中提出。,第一章 引言,FDTD的主要优点：1、由于采用时域方法，其运行一次仿真的解可以覆盖很宽的频带范围；2、差分的数学意义明确简单。基于差分的离散技术易于编程和计算。,FDTD的主要缺点：、存在数值色散误差，且误差大；2、收敛慢；3、为了确保运算稳定，网格尺寸和时间步进需要满足CFL条件。,FDTD的数学基础是麦克斯韦方程组。FDTD对电磁场E、H分量在空间和时间上采取交替抽样的离散方式，应用这种离散方式将含时间变量的麦克斯韦旋度方程转化为一组差分方程，并在时间轴上逐步推进地求解空间电磁场。,有关FDTD的基本简介以及基于FD

11、TD免费仿真软件的链接参考/wiki/Finite-difference_time-domain_method,随着电子计算机技术的发展，FDTD在电子学、光学、电磁学等领域得到了广泛的应用。FDTD在电磁领域的主要应用有：1、辐射天线的分析；2、微波器件和导行波结构的研究；3、散射和雷达截面计算；4、电磁波在复杂结构及多介质中的传播研究；5、电磁兼容问题的仿真分析。,设一函数f(x),其增量 称为f(x)的一阶前向差分。,函数的微分与自变量的微分之商等于该函数的导数。导数也叫微商。,差商,1.1 差分的概念,前向差分通常是微分在离散的函数中的等效

12、运算。,差分与微分的关系：,函数f(x)在xi点的导数可表示为,f(xi)可按如下方式进行近似,1.2 导数的差分近似,1.2.1 一阶导数,1）前向差分,(1-1a),这里我们将f(xi)简化表示为fi。,这三种不同的近似表达式引入的误差，可由泰勒级数来分析。式子f(xi+x)的泰勒级数展开为,将上式移项整理得,(1-2),(1-3),类似的,f(xi-x)的泰勒级数展开为,(1-5),因此,(1-6),可见前向及后向差分只有一阶精度。并且该精度大小取决于取样周期x。,(1-4),即,根据(1-2)和（1-4）式，可得中心差分和一阶导数间的关系为,f(x)的导数可以表示为,(1-7),(1-

13、8),可见中心差分的精度为二阶精度，其大小取决于取样周期x的平方。常规的FDTD算法采用有二阶精度的中心差分公式算法。,1.2.2 二阶导数,(1-10),根据(1-2)和（1-4）式,(1-9),将上式移项整理得,各向同性线性介质中的本构关系为,2.1 麦克斯韦方程,麦克斯韦旋度方程为,(2-1a),(2-1b),(2-1c),(2-1d),第二章 FDTD算法介绍,(2-1e),(2-1f),E为电场强度,单位V/m; D为电通量密度，单位C/m2; H为磁场强度，单位A/m; B为磁通密度，单位Wb/m2; J为电流密度，单位A/m2; Jm为磁流密度，单位V/m2.,直角坐标系中（2-

14、1a） 、（2-1b）可写为如下标量方程,（2-2a）,（2-2b）,（2-2c）,（2-2d）,（2-2e）,（2-2f）,2.2 麦克斯韦方程的FDTD方程,对于任何空间和时间的函数,我们可以写成,空间中任意一点的坐标我们可以离散化表示为,(2-3),(2-4),2.2.1 Yee元胞,Kane Yee1 将空间划分为若干个立方体,该立方体被称为”Yee Cell”. 电场分量处于元胞边线的中点，磁场分量处于元胞面的中点。,2.2.2 直角坐标中FDTD：三维情形,以(2-2d)式为例，,(2-5),其基于Yee元胞的有限差分方程如下,(2-6),上式中采用了如下平均值近似, 即,(2-5

15、)式整理得,(2-7),其中,(2-8),m为场点的坐标，对于(2-6)式而言m=(i+1/2,j,k)。,计算中通常采用Yee元胞为立方体，即,另外,这时(2-6)式可写为,(2-10),其中,(2-11),(2-9),(2-13),(2-12),类似的，(2-2e)(2-2f)离散后的形式为,(2-15),(2-14),磁场的相应FDTD公式为,(2-17),(2-16),其中,(2-18),通常采用迭代的方式来利用上述FDTD差分方程计算电磁场,若已知t1=t0=nt时刻空间各处E 的值,计算t2=t1+t/2时刻空间各处H 的值,计算t1=t2+t/2时刻空间各处E 的值,根据（2-1

16、4）、（2-15）、（2-16）,根据（2-9）、（2-12）、（2-13）,2.2.3 直角坐标中FDTD：二维情形,对于二维情况，设所有物理量与z坐标无关，即 。于是由（2-1a）及（2-1b）得,(2-19),(2-20),TE波,TM波,显然二维情况下电磁场可以划分为独立的两组，即Ex,Ey,Hz为一组，称为对于az的TE波；Hx,Hy,Ez为一组，称为对于az的TM波。二维情况下的Yee元胞如下图。,（a） TM波 ；（b）TE波,(2-23),(2-22),对于TE波，Hx=Hy=Ez=0, FDTD公式为,(2-21),(2-25),对于TM波，Ex=Ey=Hz=0, FDTD公

17、式为,(2-24),(2-26),利用迭代方式，对TE波和TM波的FDTD方程进行计算的程序流程如下,若已知t1=t0=nt时刻空间各处E及(n-1/2) t空间各处H 的值,计算t2=t1+t/2时刻空间各处H 的值,计算t1=t2+t/2时刻空间各处E 的值,若已知t1=t0=nt时刻空间各处H及(n-1/2) t空间各处E 的值,计算t2=t1+t/2时刻空间各处E 的值,计算t1=t2+t/2时刻空间各处H 的值,TE波,TM波,2.2.4 直角坐标中FDTD：一维情形,对于沿Z方向传播的TEM波，介质参数和场量与x, y无关，于是麦克斯韦标量方程为,(2-27),上面两式即为前面给出

18、的(2-2d)和(2-2b)。 一维情况E、H分量空间节点取样如下：,(2-28),（2-27）、（2-28）的FDTD离散为,(2-29),如果介质无耗，即=m=0, 则上面两式简化为,(2-30),(2-31),(2-32),利用麦克斯韦旋度方程导出的差分方程进行运算时, 时间变量步长t与空间变量步长x, y, z,之间必须满足一定条件, 即Courant-Friedrich-Levy条件 (简称为CFL或Courant 条件), 否则将出现数值不稳定性, 具体表现为随计算步数的增加, 计算场量无限制地增大.,三维情况下CFL条件为：,(2-33),2.3 数值稳定性要求,一般x=y=z=

19、, 则CFL条件变为,(2-34a),二维情况下CFL条件为：,(2-34b),一维情况下CFL条件为：,(2-34c),如果电磁波所在空间的媒质特性和频率有关, 则电磁波的传播速度也是频率的函数, 这种现象称为色散. 存在色散现象的媒质称为色散媒质.即,除了真空，任何其它实际媒质都是色散的。,2.4 数值色散,FDTD法只是Maxwell旋度方程的一种近似, 即使在非色散媒质中, 其数值解的波速与光速不同, 而且随频率变化. 由FDTD运算导致的色散现象称为数值色散.,由麦克斯韦方程可导出电磁场的任意直角分量满足齐次波动方程：,(2-35),将正弦平面波的表达式,(2-36),(2-37),

20、(2-38),代入(2-35)得理想的色散关系(解析解),即,可见相速vp只和,有关,对于非色散媒质(即,不是频率的函数), vp和频率无关.,将(2-35)进行二阶偏导的差分近似后, 如,(2-39),将 (2-36) 代入得,上式是FDTD数值色散关系的一般形式,表明 1. FDTD计算中的波的传播速度c与有关，这是离散后由于近似引起的色散 2. FDTD计算中的波的传播速度c与传播方向有关，这是离散后引起的各向异性,当x0, y0, z0, t0, (2-39)转化为(2-38). 这表明取样周期越小，数值色散就越小为使差分近似所带来的数值色散足够小，一般要求空间间隔满足,(2-40),

21、x, y, z/10,参考文献1中的计算例子：用FDTD计算二维空间内一正方形理想导体对TM波的散射情况。设入射电磁波只有Ez分量和Hy分量。并设其在时域上的波形为半个周期的正弦曲线，其宽度为。正方形导体的边长为。,2.5 数值运算例子,入射波可表示为,(2-41a),(2-41b),设Ez=1。并设,设：,TM波的差分方程为：,(2-42a),(2-42b),(2-42c),将(2-41)、(2-41b)代入(2-42a)(2-42c), 并考虑具体的边界条件，我们即可求得电磁波的数值运算结果。,）无障碍物时电磁波传播情况。,n为时间增量的数目。图中振荡波形及波形宽度变宽是由有限差分运算的误

22、差造成。,）有障碍物时，在j=30处电磁波传播情况,尽管在j=30处没有遇到障碍物，但是仍然受到障碍物反射电磁波的影响。反射波与入射波的相位相反，这是由于需要满足障碍物的边界条件。,3）有障碍物时，在j=50处电磁波传播情况,在j=50处，电磁波会遇到障碍物，并产生预期的反射波。该反射波遇到右边的边界时，又会被再次反射。,目前常用的吸收边界主要有Mur吸收边界条件和完全匹配层PML(Perfectly Matching Layer)吸收边界条件。,由于计算机容量的限制, FDTD只能处理有限区域的电磁问题对于开放域，需要在计算区域的截断边界处给出吸收边界条件来模拟电磁波连续的传播过计算空间并且

23、没有反射,第三章 吸收边界条件,吸收边界条件ABC(Absorbing Boundary Condition)是FDTD中的关键技术,二维空间波动方程为,其平面波的解为,(3-1),3.1 Engquist-Majda吸收边界条件（单向波吸收）,其中,(3-2),(3-3),考虑右图所示情况,设x=0平面为截断边界，在x0的右侧区域有入射波和反射波，,(3-4),f-为左行波，即入射波；f+为右行波，代表反射波。将(3-4)代入(3-1)，同时保留对x的导数，得：,定义微分算子,(3-5),(3-6),L可因式分解为,(3-7),L-称为左行波算子，L+为右行波算子。则(3-5)式可写为,要使

24、x=0处截断面处右行波、即反射波为0，需在截断面出设置,(3-8),(3-9),即,(3-10),根据(3-3)和(3-4)式，有,(3-11),(3-12),为了使(3-10)由频域转换到时域，将(3-11)代入(3-10)得,这就是Engquist Majda吸收边界条件。Engquist 和Majda证明，若满足(3-10)的边界条件，则以任意角度入射截断面的平面波都会被边界吸收。,类似地，对截断区域在讨论区域的右侧，若截断面位于x=b处，则相应的吸收边界条件为,(3-13),(3-10)和(3-12)式中算子的根号部分不适合直接进行数值计算。在时域有限差分法中，通过对根号部分取近似而实

25、现能够实际执行的吸收边界条件。即对根号部分进行Taylor展开，取其前一项为一阶近似，即,(3-10)的一阶近似为,(3-15),3.2 近似吸收边界条件,(3-14),3.2.1 一阶近似,若取根号部分Taylor展开的前两项，即,依据(3-16),，(3-10)可近似为,(3-17),(3-16),3.2.2 二阶近似（Mur吸收边界条件）,乘以jk，得,(3-18),利用(3-11)过渡到时域，可得,此为Mur建议的具有二阶近似，适用于二维问题的近似吸收边界条件。,(3-19),由于采取了近似，使得当入射角0时，吸收边界不能完全吸收入射波，从而产生残留的反射波。反射波随着增大而增大。,由

26、于采取了近似，采用Mur吸收边界条件进行FDTD运算时，只有当入射角较小时才能取得较为精确的结果。,为了克服这一缺点，Berenger(1994年)提出在截断边界处设置一特殊介质层（即完全匹配层，Perfectly Matched Layer），其特点为： PML为虚拟的有一定厚度（一般设为3-5个网格尺寸）的有耗介质，透射波会迅速衰减为零。 PML的波阻抗与相邻介质（自由空间）波阻抗完全匹配，这样入射波在任意角度下可无反射地进入PML。,3.2 完全匹配层PML(Perfectly Matched Layer),Berenger证明当介质的电导率和导磁率满足以下关系时,该介质为PML。,(3-20),从源所产生的电磁场结构来分类，主要分为近场源（如电压源、电流源）和远场源两大类。我们这里只讨论远场源。,对激励源的模拟是FDTD方法分析电磁问题的一个重要任务。,第四章 激励源,从源随时间变化看，主要有两类激励源：一类是随时间周期变化的时谐场源，一类是对时间呈脉冲函数形式的波源。,从空间分布来看，有面源、线源和点源等。,为了用FDTD方法来计算时谐场