時域有限差分法FDTD算法的基本原理及仿真

1、時域有限差分法(FDTD 算法)時域有限差分法是1966年K.S.Yee發(fā)表在AP上的一篇論文建立起來的， 后 被稱為Yee網(wǎng)格空間離散方式。這種方法通過將Maxwell旋度方程轉(zhuǎn)化為有限差 分式而直接在時域求解，通過建立時間離散的遞進序歹0,在相互交織的網(wǎng)格空 間中交替計算電場和磁場。FDTD算法的基本思想是把帶時間變量的Maxwell旋度方程轉(zhuǎn)化為差分形 式，模擬出電子脈沖和理想導體作用的時域響應(yīng)。需要考慮的三點是差分格式、 解的穩(wěn)定性、吸收邊界條件。有限差分通常采用的步驟是：采用一定的網(wǎng)格劃分 方式離散化場域；對場內(nèi)的偏微分方程及各種邊界條件進行差分離散化處理，建立差分格式，得到差分方程

2、組；結(jié)合選定的代數(shù)方程組的解法，編制程序，求邊 值問題的數(shù)值解。1. FDTDFDTD的基本原理FDTD方法由Maxwell旋度方程的微分形式出發(fā)，利用二階精度的中心差分 近似，直接將微分運算轉(zhuǎn)換為差分運算，這樣達到了在一定體積內(nèi)和一段時間上 對連續(xù)電磁場數(shù)據(jù)的抽樣壓縮。Maxwell方程的旋度方程組為：H m 1t(1)式可化為如下六個標量方程：；HzHy：Ex:V；:tHx：Hz _:Ey:zt:x；：tHy-Hxx=：Ez:x；：t-Ex； EzOEy.:Ez：Ey, ；:HX_- - = - mHx一y ：z:tM 一玉一土一。Hm y:z : x :t：Ey ：Ex I：H- -=r

3、*-(2)-：tz二mHz上面的六個偏微分方程是FDTDYee首先在空間上建立矩形差分網(wǎng)格，在時刻F(x,y,z,t) = F(i x, j y,k z,n t) = Fn(i,j,k)用中心差分取二階精度：對空間離散：算法的基礎(chǔ)。Mt時刻，F(xiàn)(x,y,z)可以寫成(3)-：F(x,y,z,t):x：F(x,y,z,t):yy n.yFn(i,j 1 2,k) - Fn(i,j -1 2,k) . o，-.：y2, fE L H =七ft在直角坐標系中，(DF(x,y,z,t) zL Fn(i,j,k 12)-Fn(i,j,k -12)Oz2zFn12(i，j，k)-Fn2(i，j，k)O廣t

4、Yee把空間任一網(wǎng)格上的E和H的六個分量，如下圖放置圖1 Yee氏網(wǎng)格及其電磁場分量分布在FDTD中，空間上連續(xù)分布的電磁場物理量離散的空間排布如圖所示。由圖可見，電場和磁場分量在空間交義放置，各分量的空間相對位置也適合丁Maxwell方程的差分計算，能夠恰當?shù)孛枋鲭姶艌龅膫鞑ヌ匦浴Ｍ瑫r，電場和磁 場在時間上交替抽樣，抽樣時間間隔相差半個時間步，使Maxwell旋度方程離散 以后構(gòu)成顯式差分方程，從而可以在時間上迭代求解，而不需要進行矩陣求逆運 算。因此，由給定相應(yīng)電磁問題的初始條件，F(xiàn)DTD就可以逐步推進地求得以后 各個時刻空間電磁場的分布。根據(jù)這一原則可以寫出六個差分方程：1；（i 1/2

5、，j，k）tE京幻1 .。 右，J，w i2（i 1/2，j，k）At_ 1（5）不1/2，j，k）.1 .i（my2 （i 1/2, j，k）H/2（i+1/2，j，k）-H/2（i+1/2，j -1/2，k）Hy*/2（i+1/2j，k-1/2）-H廠/2（i + 1/2，j，k + 1/2）yz其余的也如法可以寫出，每個網(wǎng)格點上的個場分兩的新值依賴丁該點在前一 時間步長時刻的值機該點周圍的臨近點上另一場量在早半個時間步長時的值。因此任一時刻可一次算出一個點，并行算法可計算出多個點。通過這些運算可以交 替算出電場磁場在各個時間步的值。根據(jù)上述FDTD差分方程組可得出計算電磁場的時域推進計算

6、方法，如圖2已知& f =以0時刻空間各處的電磁場初始對時間離散:cF(x, y,z,t)(4)所示。2.數(shù)值穩(wěn)定性條件時間步長At,空間步長Ax ,留,Az必須滿足一定的關(guān)系，否則就使得數(shù)值表現(xiàn)不穩(wěn)定，表現(xiàn)為：隨著計算步數(shù)的增加，計算場量的數(shù)值會無限的增大，這種增大不是由丁誤差積累造成的，而是由丁電磁波的傳播關(guān)系被破壞造成的。所以At, Ax , Ay , Az必須滿足一定的關(guān)系以保證穩(wěn)定性。Taflove等在1975年對Yee氏差分格式的穩(wěn)定性進行了討論，并導出了對時間步長的限制條件。數(shù) 值解是否穩(wěn)定主要取決丁時間步長At與空間步長Ax、Ay、Az的關(guān)系。對丁非均勻媒質(zhì)構(gòu)成的計算空

7、間選用如下的穩(wěn)定性條件:s石木用均勻立方體網(wǎng)格：Ax = Ay = Az = As ,M v3x而一般?。篈t = C為光速。2c3.數(shù)值色散FDTD網(wǎng)格中，會導致數(shù)字波模在網(wǎng)格中發(fā)生改變，這種改變是由丁計算網(wǎng) 格本身引起的，而非物理因素，所以必須考慮。即在FDTD網(wǎng)格中，電磁波的 相速與頻率有關(guān)，電磁波的相速度隨波長、傳播方向及變量離散化的情況不同而 改變。色散將導致非物理因素引起的脈沖波形畸變、人為的各向異性和虛假折射 等現(xiàn)象。顯然，色散與空間、時間的離散間隔有關(guān)，如下式所示：圖2 FDTD在時域的交叉半步逐步推進計算.-:t 1v/1 * 3x)2(1y)2(1z)2(6)(7)當Ax

8、, Ay , Az不相等時，At =min（Ax，釵盆）2c(8)(9)sin2心2與數(shù)值色散關(guān)系相對應(yīng)，在無耗介質(zhì)中的單色平面波，色散解析關(guān)系是:仙/c )2 = k：+ k，k2由式(9)可知，當式(9)中的At、AX、Ay、AZ均趨丁零時，它就趨丁式(10)。也就是說數(shù)值色散是由丁用近似差分替代連續(xù)微分而引起的，而且在理論上可以減小到任意程度，只要此時時間步長和空間步長都足夠小。為獲得理 想的色散關(guān)系，問題空間分割應(yīng)按照小丁正常網(wǎng)格的原則進行。一般選取的最大空間步長為Amax = yin/20 ,f 為所研究范圍內(nèi)電磁波的最小波長。由上分析說明，數(shù)值色散在用FDTD法分析電磁場傳播中的影

9、響是不可能避免的，但我 們可以盡可能的減小數(shù)值色散的影響?，F(xiàn)在適當選取時間和空間步長，傳播方向，可以得到理想情況，如下所示：3-D方形網(wǎng)格：(數(shù)值穩(wěn)定的極限狀態(tài)，可得理想色散關(guān)系)取波沿對角線傳播kx= ky= kz= k/J3 , Ax = Ay = & = 8qt=臭欄(11)_._ L 8 g.2-D萬形網(wǎng)格：也是沿對角線傳播kx = ky =kz = k/J2 , At = -(12)1-D網(wǎng)格:4.吸收邊界條件在電磁場的輻射和散射問題中，邊界總是開放的，電磁場占據(jù)無限大空間， 而計算機內(nèi)存是有限的，所以只能模擬有限空間。即：時域有限差分網(wǎng)格將在某 處被截斷。這要求在網(wǎng)格截斷處

10、不能引起波的明顯反射，因而對向外傳播的波而言，就像在無限大的空間傳播一樣，一種行之有效的方法是在截斷處設(shè)置一種吸 收邊界條件。使傳播到截斷出的波被邊界吸收而不產(chǎn)生反射。下面只給出Engquist-Majda吸收邊界條件，采用Mur差分格式，其總體虛 假反射在1%5%之間。一維一階近似情形，x=0邊界：u(0)=un+件當u(1)-u(0) (14) c t x二維二階近似情形，x=0邊界：Wn1(0, j) = -W(1, j)王牛Wn1(1, j) (0, j)c t ：x2+2 xWn(0, j) +Wn(1, j) +(ct)x.(15)ct己x2(己y) (ct己x)Wn(0 j 1) 2Wn(0 j) Wn(0 j -1) Wn(1 j 1) 2Wn(1 j) Wn(1 j T),三維二階近似情形，x=0邊界：TVsin一=7T7Sin - TTsm(10)(13)Wn1(0,j,k) = W % j,k) W01(1, j,k)一時(0, j, k) c ：t以2 -+Wn(0,j,kWn(1,j,k) .( 16)c t：=x2( y) (c t：x)Wn(0,j 1,k2Wn(0,j,k) Wn(0,j -1,k) W(1,j 1,k)-2Wn(1, j,k) Wn(1, j -1,k) Wn(0,j,k 1)-2Wn(0,j,k)