
《fdtd基本原理》由会员分享,可在线阅读,更多相关《fdtd基本原理(43页珍藏版)》请在文档大全上搜索。
1、 FDTD基本原理直角坐标系中的时域有限差分方法直角坐标系中的时域有限差分方法 简单地说,将麦克斯韦方程组中的旋度方程分解为6个(电磁场各三个)标量方程以后,运用差分近似来代替各微分运算,即可得到麦克斯韦方程的时域有限差分计算格式。 首先回顾麦克斯韦方程组: 0 BDJtDHJtBEem1. 一一维情形维情形(1) FDTD基本原理其旋度方程可分别写为:),(),(),(),(),(),(tRJtREtRBttRJtRHtRDtme ( 2)相应的本构关系为: ),(),(),(),(tRHtRBtREtRD( 3)则可以得到: ),(),(),(),(),(),(tRJtREtRHttRJt
2、RHtREtme( 4) FDTD基本原理特殊地,真空中的本构关系为: ),(),(),(),(00tRHutRBtREtRD),(),(),(),(),(),(00tRJtREtRHttRJtRHtREtme则可知:( 6)( 5)这种电磁场之间的耦合关系可以用图1表示为: FDTD基本原理 源法拉第定律安培-麦克斯韦电路定律),(),(),(tRJtREtRHtm),(),(),(tRJtRHtREte场 图1 电磁场耦合变化关系示意图源法拉第定律安培-麦克斯韦电路定律),(),(),(tRJtREtRBtm),(),(),(tRJutRBtREte场 FDTD基本原理 为简单起见,我们首
3、先考虑1维情况下麦克斯韦旋度方程以及其FDTD差分计算。由真空中的本构关系可得到: ),(),(),(),(),(),(),(),(),(),(00tRJtREtRBttRJtRHtRDttRHtRBtREtRDme( 7)可得到1维麦克斯韦方程组中的旋度方程为:),(1),(1),(),(1),(1),( ),(),( ),(),(0000tzJtzEztzHttzJtzHztzEtetzHtRHetzEtREmyxyexyxyyxx(8b)(8a)设: FDTD基本原理 对上式中的空间变量以及时间变量分别利用具有二阶精度的中心差分近似:)()2()2()()()2()2()(22tOttt
4、fttfdtxdfzOzzzfzzfdzzdf(9b)(9a) 则由(8a)式,将磁场分量与电场分量在不同位置处离散,如图2所示。 )21(znyH)(znxE2z2zzn21zn21zn)21(znyH)21(znyH)(znxE)1(znxE2z2zzn1zn21zn(a) 电场与其周围磁场 (b) 磁场与其周围电场 FDTD基本原理图 2 电磁场空间离散位置示意图 )2(znxE)2/3(znyH)2/3(znyH)2/1(znyH)2/1(znyH)1(znxE)(znxE)1(znxE2z2z2z2z2z2z2z2zn2/ 3zn1zn2/ 1znzn1zn2/ 1zn2/ 3zn(
5、c)电磁场分量在空间上的分布 =1,., (10a)znzEzx:znznzHzy)2/1(:znzN=1,., (10b)zN由( 8)- ( 9)可得: FDTD基本原理 (14 )(1)()(11)(2101021tJtEtEztHtzzzznmynxnxny)(1)()(11)(021210tJtHtHztEtzzzzneynynynx(11)(12)(13)()() 1(1),(2)21(zOnEnEztzEzztxztxznzxz)()21()21(1),(2zOnHnHztzHzztyztyznzyz FDTD基本原理 下面我们接着讨论(12)与(14)对时间变量的微分运算。如图
6、3所示,若我们令磁场的采样时刻为: tntEtx:tntNtntHty)2/1(:tntN=1,., (15a)=1,., (15b)2zn1zn2/1znzn1zn2/1zn2/3zn),2(tznnxE) 2/ 1, 2/ 3(tznnyH) 2/ 1, 2/ 3(tznnyH) 2/ 1, 2/ 1(tznnyH)2/ 1, 2/ 1(tznnyH), 1(tznnxE),(tznnxE), 1(tznnxE2z2z2z2z2z2z2zz2/3zn图 3 电磁场时间离散示意图 FDTD基本原理 由(9b)可得: )()(2),()1,()(tOtEEtEttztzznnxnnxnx)()
7、(2)21,21()21,21()21(tOtHHtHttztzznnynnyny( 17)( 18) 由(12)与(17)可得(8a)的差分形式:)(1)()(11021210),()1,(tJtHtHztEEzzztztzneynynynnxnnx( 19)由(14)与(18)可得(8b)的差分形式: FDTD基本原理 )(1)()(1121010)21,21()21,21(tJtEtEztHHzzztztznmynxnxnnynny ( 20)由(19)整理可得:)21,(0)21,21()21,21(0),()1,(tztztztztznnexnnynnynnxnnxJtHHztEE(
8、 21) ),21(0),(), 1(0)21,21()21,21(tztztztztznnmynnxnnxnnynnyJtEEztHH( 22) 至此,(21)与(22)构成了1维麦克斯韦方程组的显式差分运算方程。 由(20)整理可得: FDTD基本原理 这种电磁场分量的空间取样方式不仅符合法拉第感应定律和安培环路定律的自然结构,而且这种电磁场量的空间相对位置也适合于麦克斯韦方程的差分计算,能够恰当的描述电磁场的传播特性。此外,电磁场在时间顺序上交替抽样,抽样时间间隔彼此相差半个时间步,使得麦克斯韦旋度方程离散以后构成了显式差分方程,从而可以在时间上迭代求解。这种电磁场量随时间交替变化可用图
9、4形象表示。 因此,由给定相应电磁问题的初始值,FDTD方法就可以逐步推进地求得以后各个时刻空间电磁场的分布。我们仍以1维麦克斯韦旋度方程的求解为例,其电磁场分量的空间与时间离散如图5所示。 FDTD基本原理tntttnxE1tnxE21tnyH21tnyHtnxE21tnyH图 4 电磁场量交替迭代的蛙跳算法 FDTD基本原理(离散)显式一维FDTD在时空交错网格中的蛙跳算法 )2/1(tn)1(tnxEyH)2/1(zn)2/1(zn)1(zn)(zn)2/3(zn)1(zn)2/3(znxE)2(zn)(zn)1(zn)1(zn)2(zn)(tn)2/1(tnyH)2/3(zn)2/1(
10、zn)2/1(zn)2/3(zn图 5 蛙跳算法示意图 FDTD基本原理 对(21,22)式引入归一化量,则可得到归一化场量的离散方程式: )21,()21,21()21,21(),() 1,(ntnzexntnzyntnzyntnzxntnzxJtHHtEE(23) ),21(),(), 1()21,21()21,21(ntnzmyntnzxntnzxntnzyntnzyJtEEtHH, , ;tttrefrefrefrefcxttcxtrefref , , , , ;zxzrefcccrefrefref0ref , , ;xrefxEEEyrefyHHHrefrefrefrefrefref
11、refrefrefrefrefrefrefZEEEcEH其中, FDTD基本原理exerefexJJJrefrefreferefEtJmymrefmyJJJrefrefrefrefrefrefmrefctEHtJ则可得到归一化场量的离散方程式: (24)21,()21,21()21,21(),()1,(ntnzexntnzyntnzyntnzxntnzxJtHHtEE (25) ),21(),(), 1()21,21()21,21(ntnzmyntnzxntnzxntnzyntnzyJtEEtHH其具体的计算流程图6所示。 FDTD基本原理图 6 一维FDTD算法流程图 开始电磁场分量赋初值结
12、束计算模拟区域内的磁场分量:21ttnnttNn YesNo21ttnn)() 1()21()21(2121znxznxznyznynEnEztcnHnHtttt源激励边界条件计算模拟区域内的电场分量:)21()21()()(21211znyznyznxznxnHnHztcnEnEtttt FDTD基本原理2. 三维情形三维情形 Yee算法算法 下面我们讨论三维FDTD 离散方程的推导。由麦克斯韦旋度方程 ( 26),(),(),(),(),(),(tRJtRHtRDttRJtREtRBtem直角坐标系下旋度算子运算可以写为zxyyzxxyzzyxzyxeytRExtREextREztREez
13、tREytREtREtREtREzyxeeetRE),(),(),(),(),(),(),(),(),(),( 27) FDTD基本原理( 28)zxyyzxxyzzyxzyxeytRHxtRHextRHztRHeztRHytRHtRHtRHtRHzyxeeetRH),(),(),(),(),(),(),(),(),(),(把后面两项考虑进去,则),(),(),(tRJtREtRBtm)(),(),(),(),(),(),(),(),(),(zmzymyxmxzxyyzxxyzzzyyxxeJeJeJeytRExtREextREztREeztREytREetRBetRBetRBt( 30) (
14、 29) FDTD基本原理 ( 31),(),(),(tRJtRHtRDte( 32)(),(),(),(),(),(),(),(),(),(zezyeyxexzxyyzxxyzzzyyxxeJeJeJeytRHxtRHextRHztRHeztRHytRHetRDetRDetRDt 若将矢量旋度方程中各方向分量独立写出,我们可得六个标量方程: ( 33a),(),(),(),(tRJztREytREtRBtmxyzx ( 33b),(),(),(),(tRJxtREztREtRBtmyzxy FDTD基本原理 ( 33c),(),(),(),(tRJytRExtREtRBtexxyz ( 33
15、d),(),(),(),(tRJztRHytRHtRDtexyzx ( 33e),(),(),(),(tRJxtRHztRHtRDteyzxy ( 33f),(),(),(),(tRJytRHxtRHtRDtezxyz根据本构关系式: FDTD基本原理 ),(),(tRHtRBxx),(),(tREtRDxx ( 34),(),(tRHtRByy),(),(tREtRDyy ),(),(tRHtRBzz),(),(tREtRDzz ( 35b),(),(),(),(tRJxtREztREtRHtmyzxy ( 35c),(),(),(),(tRJytRExtREtRHtexxyz ( 35d)
16、,(),(),(),(tRJztRHytRHtREtexyzx ( 35e),(),(),(),(tRJxtRHztRHtREteyzxy FDTD基本原理 ( 35f),(),(),(),(tRJytRHxtRHtREtezxyz 式(35)中六个耦合的偏微分方程形成了FDTD方法求解电磁波与物体相互作用的算法基础。 为了将上述公式在空间和时间上离散化,K.S.Yee首先将空间按立方体分割,因此这种分割计算空间的网格通常称为Yee网格。这是实现时域有限差分的关键。电磁场的六个分量在空间的取样点分别放在立方体的边沿和表面中心点上,电磁场通过电场和磁场的耦合传播,如图 7所示。 FDTD基本原理
17、xyzxExExEyHzHxHyEyEyEzEzEzEzExEyEyHxHzH图 7 基本Yee网格以及其中的电磁场配置 FDTD基本原理 FDTD基本原理 从图 7中可以看出,各个电磁分量配置在Yee网格的特殊位置上:电场分量位于网格棱边中心并且平行于棱边,每个电场分量环绕有四个磁场分量;磁场分量位于网格面中心并且垂直于这个面,每个磁场分量环绕有四个电场分量。在空间取样上,电场和磁场分量在任何方向上始终相差半个网格步长;在时间取样上,磁场分量与电场分量相互错开半个时间步。这种场量配置不仅允许旋度方程作中心差分近似,也满足在网格上执行Faraday定律和Ampere定律的自然几何结构,因而能恰
18、当地模拟电磁波传播,而且可以自然满足媒质边界面上连续性条件。时域有限差分实际上就是在空间和时间上离散取样电磁场。 假设网格元顶点坐标 可记为:),(zyx ( 36),(),(zkyjxikji FDTD基本原理 其中 分别表示在 坐标方向网格步长, 为整数。zyx,zyx,kji, 在时间上,取n时刻的时间步为 , 为时间步长。电场分量 在时刻取样,而磁场分量在与电场相差半个时间步长处取样,即磁场的取样点为 。根据时间和空间网格划分的规律,任意一个空间和时间的函数可表示为tnttn)21(21 ( 37),(),(zkyjxiFkjiFtnn于是可以得到,磁场和电场的取样值分别为 )21,2
19、1,(21kjiHnx, )21,21(21kjiHny, , , , ),21,21(21kjiHnz),21(kjiEnx),21,(kjiEny ,考虑到时间上 和 有半个时间步的变化,按照)21,(kjiEnzEHtnttntn FDTD基本原理Yee网格上的电磁场量配置,需要利用一阶导数的二阶精度中心差分近似式 ( 38)()2/()2/()(2xOxxxfxxfxxf 其中 是步长间隔。Yee对时间导数和空间导数采用中心差分格式近似,这样程序简单,而且具有二阶精度。 将(35)标量方程中的电磁场时间和空间导数利用(38),得到各个电磁场分量的时域有限差分方程。 x以第一式为例:,
20、( 39),(),(tRHtRHtxx),(),(),(),(tRJztREytREtRHmxyzx ( 40)(),()(tHtRHmxx)(),()(tJtRJmmxmx记 FDTD基本原理 ( 41) )()()(),(2)()(yOytEtEytREbzfzz)()()(),(2)()(zOztEtEztREuydyy)()()()()()()()()()()()(tJztEtEytEtEtHmmxuydybzfzmx差分旋度算子的一部分( 42)( 43) 其它分量可以类似得到,由此即可得到 微分方程的差分离散形式。对于无源情形,则各个分量具体为: Maxwell FDTD基本原理z
21、y)(dyE)(mxH)(dyE)( fzE)(uyE)(bzE)( fzE)(bzE)(uyE)(mxHx 图 8 电磁场量形成的交链图 FDTD基本原理 ( 44a)zkjiHkjiHykjiHkjiHmCBkjiEmCAkjiEnynynznznxnx)21,21()21,21(),21,21(),21,21()(),21()(),21(212121211 ( 44b)xkjiHkjiHzkjiHkjiHmCBkjiEmCAkjiEnznznxnxnyny),21,21(),21,21()21,21,()21,21,()(),21,()(),21,(212121211 FDTD基本原理
22、( 44c)ykjiHkjiHxkjiHkjiHmCBkjiEmCAkjiEnxnxnynynznz)21,21,()21,21,()21,21()21,21()()21,()()21,(212121211 ( 44d)ykjiEkjiEzkjiEkjiEmCQkjiHmCPkjiHnznznynynxnx)21,()21, 1,(),21,() 1,21,()()21,21,()()21,21,(2121 FDTD基本原理 ( 44e) zkjiEkjiExkjiEkjiEmCQkjiHmCPkjiHnxnxnznznyny),21() 1,21()21,()21, 1()()21,21()
23、()21,21(2121 ( 44f)xkjiEkjiEykjiEkjiEmCQkjiHmCPkjiHnynynxnxnznz),21,(),21, 1(),21(), 1,21()(),21,21()(),21,21(2121 FDTD基本原理)(),(),(),(mCQmCPmCBmCA 在上组公式中,系数 分别表示各场量位置处的媒质参数, 为各场量所对应的数组数组下标下标。m)(2)(1)(2)(12)()(2)()()(mtmmtmmtmmtmmCA)(2)(1)(2)()(1)(mtmmtmtmmCB ( 45)(2)(1)(2)(12)()(2)()()(mtmmtmmtmmtmm
24、CPmmmm)(2)(1)(2)()(1)(mtmmtmtmmCQmm FDTD基本原理 上述时域有限差分方程表明,任何时刻的电磁场取决任何时刻的电磁场取决于上一时间步的电磁场、与此电磁场正交的面上前半个时于上一时间步的电磁场、与此电磁场正交的面上前半个时间步相邻的磁电场以及媒质参数。间步相邻的磁电场以及媒质参数。由于采用了中心差分近似,时域有限差分方程在空间和时间上具有二阶精度。 作为时域方法,时域有限差分法把所有研究的电磁问题作为初值问题,初始时刻模拟区域内的电磁场为零,在源激励下,以蛙跳的方式迭代时域有限差分方程,在时间上逐步向前推进电场和磁场。随着时间的发展,在有限计算区域内,时间和空
25、间上离散取样电磁场量,数值模拟电磁波传播以及与媒质间的相互作用,近似实际连续的电磁波,获得整个计算区域内时域电磁信息。 FDTD基本原理 为了方便编程计算,用计算机语言表示上述六个差分方程,得到: ( 46a)zkjiHkjiHykjiHkjiHmCBkjiEmCAkjiEnynynznznxnx 1, 1,)(,)(,1 ( 46b)xkjiHkjiHzkjiHkjiHmCBkjiEmCAkjiEnznznxnxnyny, 1, 1,)(,)(,1 ( 46c)ykjiHkjiHxkjiHkjiHmCBkjiEmCAkjiEnxnxnynynznz), 1,(,), 1(,)(,)(,1 F
26、DTD基本原理 ( 46d)ykjiEkjiEzkjiEkjiEmCQkjiHmCPkjiHnznznynynxnx, 1, 1,)(,)(,1 ( 46e)zkjiEkjiExkjiEkjiEmCQkjiHmCPkjiHnxnxnznznyny, 1, 1)(,)(,1 ( 46f)xkjiEkjiEykjiEkjiEmCQkjiHmCPkjiHnynynxnxnznz, 1, 1,)(,)(,1 图9形象的给出了“计算空间的Yee网格离散化电磁场分量空间分布编程实现时各分量表示”这一实现过程 FDTD基本原理有助于我们清晰的理解程序上FDTD的实现方法。 ) 1,(kjixE) 1,(kj
27、iyE), 1,(kjizE),(kjixH),(kjiyH),(kjizHCube(i,j,k),(kjixE),(kjiyE),(kjizE), 1,(kjixE), 1(kjizE), 1(kjiyEzyxCube(i,j,k)(i,j,k+1)(i,j+1,k+1)(i+1,j,k+1)(i+1,j,k)(i,j+1,k)(i+1,j+1,k)(i,j,k)21,(kjizE)21,21,(kjixH)21,21(kjiyH),21,21(kjizH),21,(kjiyE),21(kjixExzyzyx图9三维空间内Yee网格的划分 依据这些三维空间内差分方程组可得出计算电磁场的时域推
28、进计算方法,其流程图可以参考图6给出的流程图。 FDTD基本原理3. 媒质不连续性处理媒质不连续性处理 介质交界面的边界条件要求电场E和磁场H在穿过边界时连续。假设计算目标非磁性材料,即 ,因此交界面总位于电场分量交界面。相应的电场迭代方程也应该做些修正,以考虑相邻网格不同的介电常数和电导率。0),(zyx 接下来我们推导有四种不同介质网格包围的电场Ex分量的迭代公式。假设环线C围绕 分量。环线C包含四个路径, ,每一个路径落入一个网格内,令S1处于网格Cube(i,j,k)内,S2处于网格Cube(i,j-1,k)内,S3处于网格Cube(i,j-1,k-1)内,S4处于网格Cube(i,j
29、,k-1)内,如10图所示:nkjixE),(4321,SSSS FDTD基本原理,22,11,33,44kjiEx,Path CCube(i,j-1,k)Cube(i,j-1,k-1)Cube(i,j,k)Cube(i,j,k-1)Interface planeszyx,22,11,33,441A2A3A4A1S2S3S4SzykjiHny,21kjiHnz,211,21kjiHnykjiHnz, 1,21kjiEnx,图 10 媒质不连续性处理 沿着路径C应用积分形式安培麦克斯韦电路定应用积分形式安培麦克斯韦电路定律律 ,可得: EHtE1 FDTD基本原理 ( 47a)sdEtEsdEt
30、EsdEtEsdEtEl dHl dHl dHl dHl dHAAAAssssC43214321)()()()(44332211利用时间导数的中心差分近似,这个方程可以近似为:( 47b)4,2,4,2,4,2,4,2, 1,1,141413131212111121212121zytkjiEkjiEkjiEkjiEzytkjiEkjiEkjiEkjiEzytkjiEkjiEkjiEkjiEzytkjiEkjiEkjiEkjiEzkjiHkjiHykjiHkjiHnxnxnxnxnxnxnxnxnxnxnxnxnxnxnxnxnznznyny FDTD基本原理重新组织上式各项,可以得到由四种不同媒质包围的电场Ex电场分量的一般迭代公式:( 47c) zkjiHkjiHykjiHkjiHttkjiEttkjiEnynynznznxnx1, 1,21,2121,212121211其中: 和。4432144321 类似可得Ey与Ez分量。概括而言,我们可以通过介电常数和电导率的平均值来修改电磁场量的迭代公式。(47c)可以自动保证切向电场连续,因为穿过交界面时磁导率没有变化,这也保证了磁场的连续性。 FDTD基本原理 不连续的磁导率可以按照类似于介电常数和电导率按照平均的方法处理,但必须指出:目前的空间网格结构不能容许介电常数与磁阻率同时不连续(1994 -Hockanson)。