時(shí)域有限差分法介紹

1、時(shí)域有限差分法(姚偉)介紹時(shí)域有限差分法(姚偉)介紹10/10時(shí)域有限差分法(姚偉)介紹伊犁師范學(xué)院碩士研究生期末核查科目：電磁波有限時(shí)域差分方法姓名：姚偉學(xué)號(hào)：學(xué)院：電子與信息工程學(xué)院專業(yè)：無(wú)線電物理時(shí)域有限差分法選題背景在多種可用的數(shù)值方法中，時(shí)域有限差分法(FDTD)是一種新近發(fā)展起來(lái)的可選方法。1966年，K.S.Yee首次提出電磁場(chǎng)數(shù)值計(jì)算的新方法時(shí)域有限差分法(FiniteDifference-TimeDomain，簡(jiǎn)稱FDTD)。經(jīng)歷了二十年的發(fā)展FDTD法才漸漸走向成熟。上世紀(jì)80年月后期以來(lái)FDTD法進(jìn)入了一個(gè)新的發(fā)展階段，即由成熟轉(zhuǎn)為被寬泛接受和應(yīng)用的階段。FDTD法是解決

2、復(fù)雜問(wèn)題的有效方法之一，是一種直接鑒于時(shí)域電磁場(chǎng)微分方程的數(shù)值算法，它直接在時(shí)域?qū)axwell旋度方程用二階精度的中心差分近似，進(jìn)而將時(shí)域微分方程的求解變換為差分方程的迭代求解。是電磁場(chǎng)和電磁波運(yùn)動(dòng)規(guī)律和運(yùn)動(dòng)過(guò)程的計(jì)算機(jī)模擬。原則上能夠求解隨意形式的電磁場(chǎng)和電磁波的技術(shù)和工程問(wèn)題，并且對(duì)計(jì)算機(jī)內(nèi)存容量要求較低、計(jì)算速度較快、尤其合用于并行算法?，F(xiàn)在FDTD法己被寬泛應(yīng)用于天線的解析與設(shè)計(jì)、目標(biāo)電磁散射、電磁兼容、微波電路和光路時(shí)域解析、生物電磁劑量學(xué)、瞬態(tài)電磁場(chǎng)研究等多個(gè)領(lǐng)域1。原理解析2.1FDTD的Yee元胞E,H場(chǎng)分量取樣節(jié)點(diǎn)在空間和時(shí)間上采取交替排布，利用電生磁，磁生電的原理DEBH

3、HtEttt圖1Yee模型如圖1所示，Yee單元有以下特點(diǎn)2：1）E與H分量在空間交錯(cuò)放置，相互垂直；每一坐標(biāo)平面上的E分量四周由H分量環(huán)繞，H分量的四周由E分量圍繞；場(chǎng)分量均與坐標(biāo)軸方向一致。2）每一個(gè)Yee元胞有8個(gè)節(jié)點(diǎn)，12條棱邊，6個(gè)面。棱邊上電場(chǎng)分量近似相等，用棱邊的中心節(jié)點(diǎn)表示，平面上的磁場(chǎng)分量近似相等，用面的中心節(jié)點(diǎn)表示。3）每一場(chǎng)分量自己相距一個(gè)空間步長(zhǎng)，E和H相距半個(gè)空間步長(zhǎng)4）每一場(chǎng)分量自己相距一個(gè)時(shí)間步長(zhǎng)，E和H相距半個(gè)時(shí)間步長(zhǎng)，電場(chǎng)取n時(shí)刻的值，磁場(chǎng)取n+0.5時(shí)刻的值；即：電場(chǎng)n時(shí)刻的值由n-1時(shí)刻的值獲得，磁場(chǎng)n+0.5時(shí)刻的值由n-0.5時(shí)刻的值獲得；電場(chǎng)n時(shí)刻的

4、旋度對(duì)應(yīng)n+0.5時(shí)刻的磁場(chǎng)值，磁場(chǎng)n+0.5時(shí)刻的旋度對(duì)應(yīng)(n+0.5)+0.5時(shí)刻的電場(chǎng)值，逐步外推。5）3個(gè)空間方向上的時(shí)間步長(zhǎng)相等，以保證平均介質(zhì)中場(chǎng)量的空間變量與時(shí)間變量完全對(duì)稱。應(yīng)用這種離散方式，將含時(shí)間變量的Maxwell方程轉(zhuǎn)變?yōu)橐唤M差分方程，并在時(shí)間軸上逐步推進(jìn)地求解空間電磁場(chǎng)。由電磁問(wèn)題的初值和邊界條件，就能夠逐步推進(jìn)地求解此后各時(shí)刻空間電磁場(chǎng)散布。2.2Maxwell方程FDTD的差分格式HDJ麥克斯韋第一、二方程t（1）BEJmt式中，J時(shí)電流密度，反應(yīng)電損耗，Jm是磁流密度，單位V/m2，反應(yīng)磁損耗。主要與上式對(duì)應(yīng)。各向同性介質(zhì)中的本構(gòu)關(guān)系：DEBHJEJmmH（2）

5、其中m是磁阻率，計(jì)算磁損耗的。以E,H為變量，在直角坐標(biāo)中，展開(kāi)麥克斯韋第一、二方程，分別為HzHyExExyztHxHzEyEy（3）zxtHyHxEzEzxytEzEyHxmHxyztExEzHymHy（4）zxtEyExHzmHzxyt令fx,y,z,t代表E,H在直角坐標(biāo)中的任何一個(gè)分量，離散符號(hào)取為fx,y,z,tfix,jy,kz,ntfni,j,k（5）x,y,z,t對(duì)于時(shí)間和空間的一階偏導(dǎo)數(shù)取中心差分近似為fxfyixyjyfx1fynni21,j,kfi,j21,kfnni21,j,ki,j21,k（6）fzft1fni,j,k1fni,j,kzkzz2111fni,j,kf

6、n-i,j,k22tntt12能夠看出，每一節(jié)點(diǎn)上沿某一方向場(chǎng)分量的一階偏微分能夠用在該方向上相鄰兩點(diǎn)的一階中心差商來(lái)描繪，將式(1)用一階中心差商方程取代，整理后便獲得一階差分方程，它擁有二階精度3。Yee元胞如圖1所示，規(guī)定為1）剖分節(jié)點(diǎn)與場(chǎng)分量所在棱邊中點(diǎn)不同，場(chǎng)分量的地點(diǎn)，即E,H節(jié)點(diǎn)是Yee元胞節(jié)點(diǎn)的相對(duì)地點(diǎn)，不需要獨(dú)自編碼；2）當(dāng)空間存在媒質(zhì)分界面時(shí)，場(chǎng)量自動(dòng)知足場(chǎng)的連續(xù)性條件，E1tE2t,H1tH2t電磁分量的取樣方式不單吻合法拉第電磁感覺(jué)定律和安培環(huán)路定律的自然構(gòu)造，也吻合麥克斯韋方程的差分計(jì)算。其次，時(shí)間步長(zhǎng)能夠取為電磁波流傳一個(gè)空間步長(zhǎng)所需時(shí)間的一半，因此E與H在時(shí)間次序

7、上交替抽樣，時(shí)間間隔相差半個(gè)時(shí)間步長(zhǎng)。2.3一維問(wèn)題平均平面波（TEM波）是一維問(wèn)題，設(shè)電磁波沿z軸方向流傳，則Ez0,Hz0，0,0場(chǎng)量和介質(zhì)參數(shù)均與x，y無(wú)關(guān)，即xy，麥克斯韋方程為HyExEztx（7）ExHymHyzt和HxEyEzty（8）EyHxmHxzt旋轉(zhuǎn)坐標(biāo)軸后能夠只保存一組公式4，設(shè)保存（7）Yee元胞如圖2所示ExHyz圖2一維Yee元胞差分格式為n1n1n11n11kCBmHy2kHy2kExkCAmExz22（9）n1n112121nnHyCPmHyCQmkk2k2Exk1Exz（10）如果介質(zhì)無(wú)損耗，則0,m02.4二維問(wèn)題三維平時(shí)是散射問(wèn)題，二維是TE、TM波問(wèn)題

8、，一維是TEM波問(wèn)題。在二維場(chǎng)中，所有物理量與Z坐標(biāo)無(wú)關(guān)，既/z0。于是在TE和TM波的表達(dá)式分別為HzExExytTE波（Ez0）HzEyEy（11）xtEyExHzmHzxytEzHxmHxytTM波（Hz0）EzHymHy（12）xtHyHxEzEzxyt圖3分別給出了TM波和TE波的Yee元胞圖(i,j+1/2)HxEy(i+1/2,j+1/2)(i,j+1/2)HzyEx(i,j)(i+1/2,j)Ez圖3TM波的Yee元胞圖4TE波的Yee元胞對(duì)于TE波，只需令Ez0，在z上，Hx,Hy不隨z變化，m中去掉k即可獲得：111n111n111n1n2i2ii,j,jHz,jEx2,j

9、CAmExi2CBmHz2222y式中：mi21,j（13）n11n11n111n111i,ji,jCBmHz2i,jHz2i,jEy2CAmEy2x2222式中：mi,j1（14）2n111n111Hz2iCPmHz2i,j2,j222Eyni1,j1Eyni,j1Exni21,j1Exni21,jCQm22xy式中，m1,j1（15）i22對(duì)波，只需令Hz0，在z上，Ex,Ey不隨z變化，m中去掉k，即可獲得：TMn11n111nnHx2i,jCPmHxi,jCQmi,j222yEzi,j1Ez1式中，mi,j（16）2n1i21,jn1i21,jCQm1Ezni1,jEzni,jHy2C

10、PmHy2x式中，mi21,j（17）Ezn1i,jCAmEzni,jn11n11n11n11Hy2i,jHy2i,jHx2i,jHx2i,j2222CBmxy式中：mi,j（18）為了編寫(xiě)統(tǒng)一的TE和TM波二維FDTD程序，可將描繪TE波差分公式(13)(15)中相應(yīng)的標(biāo)號(hào)整體移動(dòng)1/2，即坐標(biāo)（x,y）分別沿x和y軸方向移動(dòng)半個(gè)網(wǎng)格，并將離散時(shí)間也移動(dòng)半個(gè)時(shí)間步長(zhǎng)，式(13)(15)能夠從頭寫(xiě)為n1n111nnEx1CAmExi,ji,j2i,j222CBmHzi,j1Hzy式中：m1（19）i,j2n1i21,jn1i21,jCBm1Hzni1,jHzni,jEy2CAmEy2x式中：m

11、i21,jHzn1i,jCPmHzni,jn11n11n11n1Ey2i,jEy2i,jEx2i,jEx2i,j222CQmxy（20）12式中，mi,j（21）能夠看出，TE波的FDTD公式(19)(21)與TM波的FDTD公式(16)(18)形式相同，給編程帶來(lái)極大方便。注意TE波和TM波之間的對(duì)偶關(guān)系5，即EHmmHE這樣就能夠編寫(xiě)統(tǒng)一的計(jì)算程序了。2.5三維問(wèn)題（直角坐標(biāo)系）2.5.1電場(chǎng)時(shí)間推進(jìn)差分格式節(jié)點(diǎn)i,j,k的3個(gè)電場(chǎng)分量分別用i1,j,k、12i,j2,k、i,j,k12地點(diǎn)上的Ex,Ey,Ez表示，以式（3）中第一個(gè)公式為例：HzHyExEyzt在ttn1i,j,k的離散

12、公式為：2時(shí)間步，對(duì)節(jié)點(diǎn)i1,j,kExn1i21,j,kExni21,j,k2t11Exn1i21,j,kExni21,j,kn21,j21,kn21,j21,ki1Hz2iHz2i2,j,k2y11n11n11Hy2i,j,kHy2i,j,k2222zn121,j,k上式中的第二項(xiàng)用平均值來(lái)替代Ex2i是因?yàn)殡x散方程中電場(chǎng)的時(shí)間取樣是整數(shù)n，磁場(chǎng)的時(shí)間取樣是n+1/2，所以只能取n及n+1時(shí)電場(chǎng)的平均值。實(shí)際也證明這個(gè)平均值使FDTD算法擁有數(shù)值穩(wěn)定性。整理后，將En1i1,j,k作為未知數(shù)，其余作為迭代計(jì)算的已知數(shù)x2Exn1i21,j,kCAmExni21,j,kCBmn111n111

13、n111n111Hz2i,j,kHz2i,j,kHy2i,j,kHy2i,j,k22222222yzmi1,j,k（22）2mm1mtt22mCAmmmmt1t22mtCBm1mmmmt1t22m同理，式（3）中其余兩個(gè)公式的離散形式為Eyn1i,j21,kCAmEyni,j21,kCBmn111n111n111n11122i,j2i2iHxi,j,kHx,kHz,j,kHz,j,k22222222zxmi,j21,k（23）n1i,j,k1n1CBmEz2CAmEzi,j,k2n111n111n111n11122i2i,j2i,jHyiHyHx,kHx,k2,j,k22,j,k22222xy

14、1mi,j,k（24）2以上三式是電場(chǎng)的時(shí)間推進(jìn)計(jì)算公式。2.5.2磁場(chǎng)時(shí)間推進(jìn)差分格式節(jié)點(diǎn)i,j,k的3個(gè)磁場(chǎng)分量分別用i,j1，k1、i1,j,k122221i2,j2,k地點(diǎn)上的Hx,Hy,Hz表示，同樣，議論式（4）中第一個(gè)公式，設(shè)察看點(diǎn)x,y,z為Hx的節(jié)點(diǎn)，即在時(shí)刻tnt，對(duì)節(jié)點(diǎn)i,j21,k21的離散公式為：n111n111Hx2i,jCPmHx2i,j2,k22,k2n1ni,j,k1n1,k1ni,j1,kCQmEzi,j1,k2Ez2Eyi,j2Ey2yzmi,j1,k122mmm1t2CPmmmm1t2（25）mmtmmmtmtCQm1mmmtmt212m同理，式（4）中

15、其余兩個(gè)公式的離散形式為n111Hy2i,j,k22CQmn121,j,kCPmHy2iExni21,j,k1Exniz121,j,kEzni1,j,k1Ezni,j,k1222xmi1,j,k1（26）22n121,j21n121,j21,kHz2i,kCPmHz2iEyni1,j21,kEyni,j21,kExni21,j1,kExni21,j,kCQmxymi1,j1,k（27）22以上三式是磁場(chǎng)的時(shí)間推進(jìn)計(jì)算公式。時(shí)域推進(jìn)計(jì)算框圖（交錯(cuò)半步逐步推進(jìn)）若已知t1t0nt時(shí)空間各節(jié)點(diǎn)處的電場(chǎng)值（賦初值）計(jì)算t2t1t25）（27）時(shí)空間各節(jié)點(diǎn)處的磁場(chǎng)值，式（2計(jì)算t1t2t時(shí)空間各節(jié)點(diǎn)處的

16、電場(chǎng)值，式（22）（24）2在編程中，為了使電場(chǎng)和磁場(chǎng)有相同的數(shù)量級(jí)（為減小誤差），可對(duì)H或E進(jìn)行“歸一化”辦理，即：用E/Z0取代E，式中Z0是自由空間波阻抗。HZ0H取代H，用E6。計(jì)算結(jié)果再分別除以和乘以Z0即可。能夠看出，這種離散方法電場(chǎng)和磁場(chǎng)在時(shí)間次序上交替抽樣，抽樣間隔相差半個(gè)時(shí)間步長(zhǎng)，使麥克斯韋方程離散后成為顯示差分方程，進(jìn)而能夠在時(shí)間上迭代求解，不需矩陣求逆。給定初值后，能夠逐步推進(jìn)，求得此后各個(gè)時(shí)刻點(diǎn)的空間電磁散布。這是FDTD法的最大特點(diǎn)。2.6解的穩(wěn)定性在FDTD中，時(shí)間增量t和空間增量x、y、z之間不是相互獨(dú)立的，它們的取值必須知足一定的關(guān)系，以防備數(shù)值結(jié)果的不穩(wěn)定，表

17、現(xiàn)為隨著時(shí)間步數(shù)的增加，計(jì)算結(jié)果發(fā)散。造成解不穩(wěn)定的因素有多種：誤差因素：計(jì)算機(jī)在計(jì)算過(guò)程中，原始數(shù)據(jù)可能有誤差，如系數(shù)陣成立過(guò)程中產(chǎn)生的誤差，而每次運(yùn)算由于只能保存有限位數(shù)而又產(chǎn)生誤差，誤差的積累有可能吞沒(méi)真實(shí)解，使計(jì)算結(jié)果不可靠，即不穩(wěn)定；計(jì)算方法不合適；t、h離散間隔不當(dāng)?shù)?。為了確定數(shù)值解穩(wěn)定的條件，有很多推導(dǎo)方式，結(jié)論相同2.6.1時(shí)間步長(zhǎng)穩(wěn)定性要求一般情況下T（28）t2.6.2時(shí)間步長(zhǎng)與空間步長(zhǎng)的關(guān)系三維t11（29）222111xyz在非平均地區(qū)，v取最大值。真空中v=c(光速)。假如正方體Yee元胞，xyz，那么t1（30）3假如正方形Yee元胞（二維），xy，那么t1（）231假如線段平分Yee元胞（一維），那么t1（32）2.7數(shù)值色散當(dāng)波流傳的速度是頻次的函數(shù)，即速度與頻次相關(guān)時(shí)，稱其波為色散波。色散的原因有多種：由于媒質(zhì)是金