fdtd基本原理

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(c

5、)電磁場分量在空間上的分布 =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)對時間變量的微分運算。如圖3

6、所示，若我們令磁場的采樣時刻為： 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)()(2)21,21()2

7、1,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( 21) ),21(0

8、),(), 1(0)21,21()21,21(tztztztztznnmynnxnnxnnynnyJtEEztHH( 22) 至此，(21)與(22)構成了1維麥克斯韋方程組的顯式差分運算方程。 由(20)整理可得： FDTD基本原理 這種電磁場分量的空間取樣方式不僅符合法拉第感應定律和安培環(huán)路定律的自然結構，而且這種電磁場量的空間相對位置也適合于麥克斯韋方程的差分計算，能夠恰當?shù)拿枋鲭姶艌龅膫鞑ヌ匦?。此外，電磁場在時間順序上交替抽樣，抽樣時間間隔彼此相差半個時間步，使得麥克斯韋旋度方程離散以后構成了顯式差分方程，從而可以在時間上迭代求解。這種電磁場量隨時間交替變化可用圖4形象表示。 因此，由

9、給定相應電磁問題的初始值，F(xiàn)DTD方法就可以逐步推進地求得以后各個時刻空間電磁場的分布。我們仍以1維麥克斯韋旋度方程的求解為例，其電磁場分量的空間與時間離散如圖5所示。 FDTD基本原理tntttnxE1tnxE21tnyH21tnyHtnxE21tnyH圖 4 電磁場量交替迭代的蛙跳算法 FDTD基本原理（離散）顯式一維FDTD在時空交錯網(wǎng)格中的蛙跳算法 )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(zn)2/1(zn)2

10、/3(zn圖 5 蛙跳算法示意圖 FDTD基本原理 對(21,22)式引入歸一化量,則可得到歸一化場量的離散方程式： )21,()21,21()21,21(),() 1,(ntnzexntnzyntnzyntnzxntnzxJtHHtEE(23) ),21(),(), 1()21,21()21,21(ntnzmyntnzxntnzxntnzyntnzyJtEEtHH， ， ；tttrefrefrefrefcxttcxtrefref ， ， ， ， ；zxzrefcccrefrefref0ref ， ， ；xrefxEEEyrefyHHHrefrefrefrefrefrefrefrefrefre

11、frefrefrefZEEEcEH其中， FDTD基本原理exerefexJJJrefrefreferefEtJmymrefmyJJJrefrefrefrefrefrefmrefctEHtJ則可得到歸一化場量的離散方程式： (24)21,()21,21()21,21(),()1,(ntnzexntnzyntnzyntnzxntnzxJtHHtEE (25) ),21(),(), 1()21,21()21,21(ntnzmyntnzxntnzxntnzyntnzyJtEEtHH其具體的計算流程圖6所示。 FDTD基本原理圖 6 一維FDTD算法流程圖 開始電磁場分量賦初值結束計算模擬區(qū)域內的磁場

12、分量:21ttnnttNn YesNo21ttnn)() 1()21()21(2121znxznxznyznynEnEztcnHnHtttt源激勵邊界條件計算模擬區(qū)域內的電場分量:)21()21()()(21211znyznyznxznxnHnHztcnEnEtttt FDTD基本原理2. 三維情形三維情形 Yee算法算法 下面我們討論三維FDTD 離散方程的推導。由麥克斯韋旋度方程 ( 26),(),(),(),(),(),(tRJtRHtRDttRJtREtRBtem直角坐標系下旋度算子運算可以寫為zxyyzxxyzzyxzyxeytRExtREextREztREeztREytREtREt

13、REtREzyxeeetRE),(),(),(),(),(),(),(),(),(),( 27) FDTD基本原理( 28)zxyyzxxyzzyxzyxeytRHxtRHextRHztRHeztRHytRHtRHtRHtRHzyxeeetRH),(),(),(),(),(),(),(),(),(),(把后面兩項考慮進去，則),(),(),(tRJtREtRBtm)(),(),(),(),(),(),(),(),(),(zmzymyxmxzxyyzxxyzzzyyxxeJeJeJeytRExtREextREztREeztREytREetRBetRBetRBt( 30) ( 29) FDTD基本

14、原理 ( 31),(),(),(tRJtRHtRDte( 32)(),(),(),(),(),(),(),(),(),(zezyeyxexzxyyzxxyzzzyyxxeJeJeJeytRHxtRHextRHztRHeztRHytRHetRDetRDetRDt 若將矢量旋度方程中各方向分量獨立寫出，我們可得六個標量方程： ( 33a),(),(),(),(tRJztREytREtRBtmxyzx ( 33b),(),(),(),(tRJxtREztREtRBtmyzxy FDTD基本原理 ( 33c),(),(),(),(tRJytRExtREtRBtexxyz ( 33d),(),(),()

15、,(tRJztRHytRHtRDtexyzx ( 33e),(),(),(),(tRJxtRHztRHtRDteyzxy ( 33f),(),(),(),(tRJytRHxtRHtRDtezxyz根據(jù)本構關系式： 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首先將空間按立方體分割，因此這種分割計算空間的網(wǎng)格通常稱為Yee網(wǎng)格。這是實現(xiàn)時域有限差分的關鍵。電磁場的六個分量在空間的取樣點分別放在立方體的邊沿和表面中心點上，電磁場通過電場和磁場的耦合傳播，如圖 7所示。 FDTD基本原理xyzxExExEyH

17、zHxHyEyEyEzEzEzEzExEyEyHxHzH圖 7 基本Yee網(wǎng)格以及其中的電磁場配置 FDTD基本原理 FDTD基本原理 從圖 7中可以看出，各個電磁分量配置在Yee網(wǎng)格的特殊位置上：電場分量位于網(wǎng)格棱邊中心并且平行于棱邊，每個電場分量環(huán)繞有四個磁場分量；磁場分量位于網(wǎng)格面中心并且垂直于這個面，每個磁場分量環(huán)繞有四個電場分量。在空間取樣上，電場和磁場分量在任何方向上始終相差半個網(wǎng)格步長；在時間取樣上，磁場分量與電場分量相互錯開半個時間步。這種場量配置不僅允許旋度方程作中心差分近似，也滿足在網(wǎng)格上執(zhí)行Faraday定律和Ampere定律的自然幾何結構，因而能恰當?shù)啬M電磁波傳播，而

18、且可以自然滿足媒質邊界面上連續(xù)性條件。時域有限差分實際上就是在空間和時間上離散取樣電磁場。 假設網(wǎng)格元頂點坐標 可記為：),(zyx ( 36),(),(zkyjxikji FDTD基本原理 其中 分別表示在 坐標方向網(wǎng)格步長， 為整數(shù)。zyx,zyx,kji, 在時間上，取n時刻的時間步為 , 為時間步長。電場分量 在時刻取樣，而磁場分量在與電場相差半個時間步長處取樣，即磁場的取樣點為 。根據(jù)時間和空間網(wǎng)格劃分的規(guī)律，任意一個空間和時間的函數(shù)可表示為tnttn)21(21 ( 37),(),(zkyjxiFkjiFtnn于是可以得到，磁場和電場的取樣值分別為 )21,21,(21kjiHnx

19、， )21,21(21kjiHny， ， ， ， ),21,21(21kjiHnz),21(kjiEnx),21,(kjiEny ，考慮到時間上 和 有半個時間步的變化，按照)21,(kjiEnzEHtnttntn FDTD基本原理Yee網(wǎng)格上的電磁場量配置，需要利用一階導數(shù)的二階精度中心差分近似式 ( 38)()2/()2/()(2xOxxxfxxfxxf 其中 是步長間隔。Yee對時間導數(shù)和空間導數(shù)采用中心差分格式近似，這樣程序簡單，而且具有二階精度。 將(35)標量方程中的電磁場時間和空間導數(shù)利用(38)，得到各個電磁場分量的時域有限差分方程。 x以第一式為例：， ( 39),(),(t

20、RHtRHtxx),(),(),(),(tRJztREytREtRHmxyzx ( 40)(),()(tHtRHmxx)(),()(tJtRJmmxmx記 FDTD基本原理 ( 41) )()()(),(2)()(yOytEtEytREbzfzz)()()(),(2)()(zOztEtEztREuydyy)()()()()()()()()()()()(tJztEtEytEtEtHmmxuydybzfzmx差分旋度算子的一部分( 42)( 43) 其它分量可以類似得到，由此即可得到 微分方程的差分離散形式。對于無源情形，則各個分量具體為： Maxwell FDTD基本原理zy)(dyE)(mxH

21、)(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基本原理 ( 44c)ykjiH

22、kjiHxkjiHkjiHmCBkjiEmCAkjiEnxnxnynynznz)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()()21,21(212

23、1 ( 44f)xkjiEkjiEykjiEkjiEmCQkjiHmCPkjiHnynynxnxnznz),21,(),21, 1(),21(), 1,21()(),21,21()(),21,21(2121 FDTD基本原理)(),(),(),(mCQmCPmCBmCA 在上組公式中，系數(shù) 分別表示各場量位置處的媒質參數(shù)， 為各場量所對應的數(shù)組數(shù)組下標下標。m)(2)(1)(2)(12)()(2)()()(mtmmtmmtmmtmmCA)(2)(1)(2)()(1)(mtmmtmtmmCB ( 45)(2)(1)(2)(12)()(2)()()(mtmmtmmtmmtmmCPmmmm)(2)(

24、1)(2)()(1)(mtmmtmtmmCQmm FDTD基本原理 上述時域有限差分方程表明，任何時刻的電磁場取決任何時刻的電磁場取決于上一時間步的電磁場、與此電磁場正交的面上前半個時于上一時間步的電磁場、與此電磁場正交的面上前半個時間步相鄰的磁電場以及媒質參數(shù)。間步相鄰的磁電場以及媒質參數(shù)。由于采用了中心差分近似，時域有限差分方程在空間和時間上具有二階精度。 作為時域方法，時域有限差分法把所有研究的電磁問題作為初值問題，初始時刻模擬區(qū)域內的電磁場為零，在源激勵下，以蛙跳的方式迭代時域有限差分方程，在時間上逐步向前推進電場和磁場。隨著時間的發(fā)展，在有限計算區(qū)域內，時間和空間上離散取樣電磁場量，

25、數(shù)值模擬電磁波傳播以及與媒質間的相互作用，近似實際連續(xù)的電磁波，獲得整個計算區(qū)域內時域電磁信息。 FDTD基本原理 為了方便編程計算，用計算機語言表示上述六個差分方程，得到： ( 46a)zkjiHkjiHykjiHkjiHmCBkjiEmCAkjiEnynynznznxnx 1, 1,)(,)(,1 ( 46b)xkjiHkjiHzkjiHkjiHmCBkjiEmCAkjiEnznznxnxnyny, 1, 1,)(,)(,1 ( 46c)ykjiHkjiHxkjiHkjiHmCBkjiEmCAkjiEnxnxnynynznz), 1,(,), 1(,)(,)(,1 FDTD基本原理 ( 4

26、6d)ykjiEkjiEzkjiEkjiEmCQkjiHmCPkjiHnznznynynxnx, 1, 1,)(,)(,1 ( 46e)zkjiEkjiExkjiEkjiEmCQkjiHmCPkjiHnxnxnznznyny, 1, 1)(,)(,1 ( 46f)xkjiEkjiEykjiEkjiEmCQkjiHmCPkjiHnynynxnxnznz, 1, 1,)(,)(,1 圖9形象的給出了“計算空間的Yee網(wǎng)格離散化電磁場分量空間分布編程實現(xiàn)時各分量表示”這一實現(xiàn)過程 FDTD基本原理有助于我們清晰的理解程序上FDTD的實現(xiàn)方法。 ) 1,(kjixE) 1,(kjiyE), 1,(kj

27、izE),(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網(wǎng)格的劃分 依據(jù)這些三維空間內差分方程組可得出計算電磁場的時域推進計算方法，其流程圖可以參考圖6給出的流程圖。 FDTD基本原理3. 媒質不連續(xù)性處理媒質不連續(xù)性處理 介質交界面的邊界條件要求電場E和磁場H在穿過邊界時連續(xù)。假設計算目標非磁性材料，即 ，因此交界面總位于電場分量交界面。相應的電場迭代方程也應該做些修正，以考慮相鄰網(wǎng)格不同的介電常數(shù)和電導率。0),(zyx 接下來我們推導有四種不同介質網(wǎng)格包圍的電場Ex分量的迭代公式。假設環(huán)線C圍繞 分量。環(huán)線C包含四個路徑， ,每一個路徑落入一個網(wǎng)格內，令S1處于網(wǎng)格Cube(i,j,k)內，S2處于網(wǎng)格Cube(i,j-1,k)內，S3處于