時域有限差分法姚偉介紹

1、伊犁師范學院碩士研究生期末考核科目：電磁波有限時域差分方法姓名：姚住學號：1076411203009學院：由子與信息工程學院專業(yè)：無線電物理時域有限差分法1選題背景在多種可用的數值方法中，時域有限差分法 （FDTD）是一種新近發(fā)展起來的可選方法。 1966年，K.S.Yee首次提出電磁場數值計算的新方法一時域有限差分法（FiNte Difference- Time Domain ,簡稱FDTD）經歷了二十年的發(fā)展 FDTDt才逐漸走向成熟。 上世紀80年代后期以來FDTDt進入了一個新的發(fā)展階段，即由成熟轉為被廣泛接受和應 用的階段。FDTDt是解決復雜問題的有效方法之一，是一種直接基于時域電

2、磁場微分方程的數值算法，它直接在時域將 Maxwell旋度方程用二階精度的中心差分近似，從而將時域 微分方程的求解轉換為差分方程的迭代求解。是電磁場和電磁波運動規(guī)律和運動過程的計 算機模擬。原則上可以求解任意形式的電磁場和電磁波的技術和工程問題，并且對計算機 內存容量要求較低、計算速度較快、尤其適用于并行算法。現在FDTD法己被廣泛應用于天線的分析與設計、目標電磁散射、電磁兼容、微波電路和光路時域分析、生物電磁劑量 學、瞬態(tài)電磁場研究等多個領域1。2原理分析2.1 FDTD 的 Yee元胞E,H場分量取樣節(jié)點在空間和時間上采取交替排布，利用電生磁，磁生電的原理.D 壬=S.:t ft'

3、、E =汨 ,:H=r1.:t 2 t圖1 Yee模型如圖1所示，Yee單元有以下特點2:1） E與H分量在空間交叉放置，相互垂直；每一坐標平面上的E分量四周由H分量環(huán)繞，H分量的四周由E分量環(huán)繞；場分量均與坐標軸方向一致。2）每一個Yee元胞有8個節(jié)點，12條棱邊，6個面。棱邊上電場分量近似相等，用 棱邊的中心節(jié)點表示，平面上的磁場分量近似相等，用面的中心節(jié)點表示。3）每一場分量自身相距一個空間步長，E和H相距半個空間步長4）每一場分量自身相距一個時間步長，E和H相距半個時間步長，電場取n時刻的值， 磁場取n+0.5時刻的值；即：電場n時刻的值由n-1時刻的值得到，磁場n+0.5時刻的值 由

4、n-0.5時刻的值得到；電場n時刻的旋度對應n+0.5時刻的磁場值，磁場n+0.5時刻的旋度對應(n+0.5)+0.5時刻的電場值，逐步外推。5) 3個空間方向上的時間步長相等，以保證均勻介質中場量的空間變量與時間變量完全對稱。應用這種離散方式，將含時 間變量的Maxwell方程轉化為一組差分方程，并在時間軸上逐步推進地求解空間電磁場。 由電磁問題的初值和邊界條件，就可以逐步推進地求解以后各時刻空間電磁場分布。2.2 Maxwell方程FDTD勺差分格式VxH ="+J麥克斯韋第一、二方程 t(1)' E = - - J mT式中，J時電流密度，反映電損耗，Jm是磁流密度，單

5、位V/m反映磁損耗。主要 與上式對應。各向同性介質中的本構關系：(2)D = E B = H J = EJm = m H其中力是磁阻率，計算磁損耗的以E，H為變量，在直角坐標中，展開麥克斯韋第一、二方程，分別為-：Ey什-：Ey-Ez:t:t-HxEyEz(3)二 y-Ex:zft- mHx-Ez-：Hy：z'EyEx-yft=_ J以-:t- mH y-'mHz令f(x,y,z,t淤表E, H在直角坐標中的任何一個分量，離散符號取為(5)f x,y,z,ti；=f i x, j y,k z,n ：t )= fn i,j,kf(X, y,z,t 次于時間和空間的一階偏導數取中心

6、差分近似為 定;ifn(iT，j,k)_fn。：，"x=i4&X wf n(,j +,k)-f n(i,j -4，kjl(6),矽心 wf n(i , j ,k +； )- f n(i, j ,k 一： N友 z4& Az尸f1 n+n-1»- |f 2(i，j*f 2(i,j,k)、ct t 朝& |可以看出，每一節(jié)點上沿某一方向場分量的一階偏微分可以用在該方向上相鄰兩點的 一階中心差商來描述，將式（1）用一階中心差商方程取代，整理后便得到一階差分方程， 它具有二階精度3 oYee元胞如圖1所示，規(guī)定為1）剖分節(jié)點與場分量所在棱邊中點不同，場分量的

7、位置，即 E,H節(jié)點是Yee元胞節(jié) 點的相對位置，不需要單獨編碼；2）當空間存在媒質分界面時，場量自動滿足場的連續(xù)性條件，E1t = E2t，H1t = H2t電磁分量的取樣方式不僅符合法拉第電磁感應定律和安培環(huán)路定律的自然結構，也符 合麥克斯韋方程的差分計算。其次，時間步長可以取為電磁波傳播一個空間步長所需時間 的一半，因此E與H在時間順序上交替抽樣，時間間隔相差半個時間步長。2.3 一維問題均勻平面波（TEM波）是一維問題，設電磁波沿 z軸方向傳播，則Ez=°，Hz = ° ,cc=0，一 =0場量和介質參數均與x, y無關，即敵到,麥克斯韋方程為.:E:£；

8、：t汩Uft-+ ymH y；H y _;z：EX 二 一z-：Hx:z 王y:z,:E£ft.:H p：(8)-:t旋轉坐標軸后可以只保留一組公式4,設保留（7）Yee元胞如圖2所示差分格式為ExHy圖2 一維Yee元胞c1E；Rk)=CA(mE(k)CB(m)1 Hz1+1-h y2 k -1(9)11nn _Hy 2 k 1 : CP m Hy 2 k 2 -CQ m（i0）如果介質無損耗，則=0，m = 02.4二維問題三維通常是散射問題，二維是 TE、TM波問題，一維是TEM發(fā)問題。在二維場中，所有物理量與 Z坐標無關，既已/金=0。于是在TE和TM波的表達式分別為TE波(

9、Ez=0).:Hz：:Exz 二；一xEx.y .t.,Hz-Ey Eyx .t沱 y FEx -Hz改 ：:y：:t一 mHz(11)TM波(Hz =0)：Ez ：Hx一 mHx,t:H Vy - mHyFt)；：Ez Ezft:y-：Ez：x'Hy.x圖3分別給出了 TM波和TE波的Yee元胞圖圖3 TM波的Yee元胞=3(12)：:H；:yEy(i+1/2, j+1/2)Hz (i, j+1/2)Ex圖4 TE波的Yee元胞對于TE波，只要令巳二°,在Az上，Hx,Hy不隨z變化，m中去掉k即可得至IJ:111、n a 1n 1nnt -inn七 1ExT +2J )=

10、CA(mEx(i+2，j )+CB(m) Hz 2(i+2,j + 2 )-H z 2(i+2,j-2)八J式中：m = i +2, j(13)/njn i、!-n 書,1 ca/ u-nf. .上 1 o k 1- - 0 2 f , 1 , 1 0 2 (1 , 1 Ey q, j + 2 )=CA(m Ey i, j + 2 )-CB(m )- H z 2Q + 2 , j + 2 ) H z 2 U _ 2, j + 2 )八)式中：m =i, j +：(14)11Hz 2(i +2, j + 2 )=CP(mHz 7 +2，j + 2)E；(i +1,j +； )-E；(i,j +；

11、 ) EX(i+： , j +1)-EX(i+1 , j ) -CQ(m 1AxAy式中，m=i +2 ,j +2(15)對TM波，只要令Hz =0 ,在Az上，Ex,Ey不隨z變化，m中去掉k,即可得到:n Jn 1AhX 2 i,j ； =CPmH；W i,j .； -CQm - E； i,j 1 -E； i,j 1y式中,m=i, j +：(16)n/n 1H；消 +2，j )=CP(mH>(i +2，j)+CQ(m) 1 13+1/)E；(i,j )1x(17)式中，m =i jEL i,j =CAm E； i,j11111h；2G+.)h；2G：,j) H：2(i,j+1AHX

12、2(i,j(18)式中：m =i, j為了編寫統(tǒng)一的TE和TM波二維FDTLg序，可將描述TE波差分公式(13)(15)中相 應的標號整體移動1/2 ,即坐標(x,y )分別沿x和y軸方向移動半個網格，并將離散時間 也移動半個時間步長，式(13)(15)可以重新寫為n 1n1Ex 2 i,j J =CAmEx 2 i,j 弓 CB m H； i,j 1 -H； i,j 1(19)11n1nEy 2 i 2-,j 二CA m Ey 2 ijCBmini 1'j -Hl" 1(20)Hn1 i,j =CPmHn i,j-CQ(m)1n 5Ey 2 in 51n .j -Ey 2

13、i -7,jEx 2 i,jx1n2 1_ Ex i,j - 2y式中：m =i +；, j(21)式中，m=ij可以看出，TE波的FDT必式(19)(21)與TM波的FDT必式(16)(18)形式相同, 給編程帶來極大方便。注意TE波和TM波之間的對偶關系5 ,即這樣就可以編寫統(tǒng)一的計算程序了2.5三維問題(直角坐標系)2.5.1 電場時間推進差分格式(i +：, j,k )、(i, j +k X (i,j,k 十;位置上的 Ex,Ey,Ez節(jié)點(i,j,k的3個電場分量分別用 表示，以式(3)中第一個公式為例：必也:.y；z在"*"+："問步，對節(jié)點(i,j,

14、k )的離散公式為:ei 1 ,j,k -J i ； ,j,k -E'i ； ,j,k；,j,k - tn 工in 乂n i J jkEnijkH 2 i- j - J k-H 2 i- j - kx i 2 ,j,kExi 2 ,j,kH z i2 ,j 2 ,k， H z i2 ,j 2 ,k yn 1 iin =iiHy 2(i + 2 ,j,k %) -Hy 2 (i+ ,j,k-2A z上式中的第二項用平均值來替代 E： ”i +1, j,k)是因為離散方程中電場的時間取樣是 整數n,磁場的時間取樣是n+1/2,所以只能取n及n+1時電場的平均值。實際也證明這 個平均值使FD

15、TDJJ法具有數值穩(wěn)定性。整理后，將EXG+j,k肝為未知數，其余作為迭代計算的已知數E?卡。+； ,j,k J=CA(m£(i +； ,j,k 計CB(m). 1111n11n11n11n11. Hz xAym=i,j,k+；（24）以上三式是電場的時間推進計算公式。2.5.2磁場時間推進差分格式節(jié)點（i,j,k用勺3個磁場分量分別用（i, j+,k + ）、（i +；, j, k + ； ）、（i +1 , j十,k位置上的Hx,Hy,Hz表示，同樣，討論式（4）中第一個公式，設觀察點（x,y,z）為Hx的節(jié)點，即 在時刻t=nAt,對節(jié)點(i,j+1,k+2)的離散公式為：J1

16、 11211Hx 2 i, j 2,k 2 =CP m Hx 2 i, j 2,k 八七。+1,k+2 "E；(i,j,k + 2)E；(i,j +J,k + 1)-E”j+U-CQ(m )' (i +； ,j+； ,kHz 2 G +； ,j -1 ,k)_Hy 2 ( +；,j,k +；) Hy 2(i +； ,j,k yA zt一m =i ：, j,k；m m 1 m tCAm 一 ：t - 22； mCA m - ； m m - 1m ：t.：t 22 ； m.：tCB m =1: m；m m 一 1 . m ：tt 22 ； m同理，式(3)中其它兩個公式的離散形式

17、為Ey+(, j +2,k "CA(mE；,j +2,k 計 CB(m>1111H：4(i, j +；,k + ； )-H：*(i, j +；,k ； ) H：'(i + ； , j + 2 ,kH：%i 一； , j +；,k)zAxm =i,j ：,kE；1 i,j,k -； =CAmE；i,j,k ； CBm-1111-II n q2 / ,1 . I . 1 II nd2 /1 . I . 1 I_jn 與1 1 I 1 1 II n42 /- -1 |上1(22)(23)Hy (i +2"卜2)Hy « -2，j，k、2)Hx (i, j

18、+2 卜十2)-H? (i, j - 2，k+ 2)(25)1 Im =i,j 2，kmm m1m m tCP m .，t - 2_ _2mmmm1mm.；ttCQ m =1 ； m:m1 m rt 22 ； m同理，式(4)中其它兩個公式的離散形式為 11n o 11_n r 11Hy 2i 2,j,k 2 =CPmHy 2i 2，j,k2 E(i+；,j,k+1)-EXi +：,j,k) E；(i +1,j,k+6E；(i,j,k+：-CQ(m2222AzAxm=i+j,k+；(26)nJn 1H： 2i ；,j ；,k =CPmH：，i ；,j 晨k-E；(i+1,j +：kLE；G,j

19、 +：,k) E：(i +2，j +1,k)-EXi+2，j,k)1-CQ(m ) 2222AxAym = i +：，j +2*（27）以上三式是磁場的時間推進計算公式。時域推進計算框圖（交叉半步逐步推進）在編程中，為了使電場和磁場有相同的數量級（為減小誤差），可對H或E進行“歸一化”處理，即：用H =Z°H取代H ,用 = £/2。取代E,式中Zo是自由空間波阻抗。計算結果再分別除以和乘以Z0即可。可以看出，這種離散方法電場和磁場在時間順序上交替抽樣，抽樣間隔相差半個時間 步長，使麥克斯韋方程離散后成為顯示差分方程，從而可以在時間上迭代求解，不需矩陣 求逆。給定初值后，可以逐步推進，求得以后各個時刻點的空間電磁分布。這是FDTD法的最大特點。2.6解的穩(wěn)定性在FDTM,時間增量 組和空間增量Ax、Ay、Az之間不是相互獨立的，它們的取值必須 滿足一定的關系，以避免數值結果的不穩(wěn)定，表現為隨著時間步數的增加，計算結果發(fā)散 造成解不穩(wěn)定的因素有多種： 誤差因素：計算機在計算過程中，原始