fdtd講座1125

1、FDTD時域有限差分法一.時域有限差分法概述二.時域有限差分法的特點三.非分裂完美吸收邊界FDTD方法 的研究一一.時域有限差分法概述時域有限差分法概述1.FDTD法的發(fā)展背景法的發(fā)展背景 FDTD方法最初由K.S.Yee提出，它是一種電磁場數值計算，通過將時間和空間離散化來直接求解麥克斯韋方程組的方法。對麥克斯韋方程的偏微分形式進行中心差分。 盡管Yee早在1966年發(fā)表了他的論文，但是這種方法沒有立刻得到認可，這不僅因為當時計算機性能（CPU速度和內存等）的限制，而且缺少一種有效的高精度的吸收邊界條件（ABC）來實現自由空間的無反射截斷。每一個磁場分量由四個電場分量環(huán)繞，同樣每一個電場分量

2、由四個磁場分量環(huán)繞。電場和磁場在時間順序上交替抽樣，抽樣時間間隔彼此相差半個時間步不需要進行矩陣求逆運算，由給定相應電磁問題的初始值，就可以逐步推進地求得以后各個時刻空間電磁場的分布 這種電磁場分量的空間取樣方式不僅符合法拉弟感應定律和安培環(huán)路定律的自然結構，而且這種電磁場各分量的空間相對位置也適合于麥克斯韋方程的差分計算，能夠恰當地描述電磁場的傳播特性。 作為一種電磁場的數值計算方法，時域有限差分法具有一些非常突出的特點，也是它的優(yōu)點。 正是由于這些，使得越來越多的人對它產生了濃厚的興趣，并得到越來越廣泛的應用。 二二.時域有限差分法的特點時域有限差分法的特點 (1)直接時域計算直接時域計算

3、 時域有限差分法直接把含時間變量的 Maxwell旋度方程在Yee氏網格空間中轉換為差分方程。 在這種差分格式中每個網格點上的電場(或磁場)分量僅與它相鄰的磁場(或電場)分量及上一時間步該點的場值有關。 在每一時間步計算網格空間各點的電場和磁場分量，隨著時間步的推進，即能直接模擬電磁波的傳播及其與物體的相互作用過程。 這一特點使它能直接給出非常豐富的電磁場問題的時域信息，給復雜的物理過程描繪出清晰的物理圖像如果需要頻域信息，則只需對時域信息進行傅立葉變換，為獲得寬頻帶的信息，只需在寬頻譜的脈沖（高斯脈沖）激勵下進行一次計算 (2)廣泛的適用性廣泛的適用性 由于時域有限差分法的直接出發(fā)點是概括電

4、磁場普通規(guī)律的Maxwell方程，這就預示著這一方法應具有最廣泛的適用性。 從具體的算法看，在時域有限差分法的差分格式中被模擬空間電磁性質的參量是按空間網格給出的，因此，只需設定相應空間點以適當的參數，就可模擬各種復雜的電磁結構。 媒質的非均勻性、各向異性、色散特性和非線性等均能很容易地進行精確模擬。 由于在網格空間中電場和磁場分量是被交叉放置的，而且計算中用差分代替了微商，使得介質交界面上的邊界條件能自然得到滿足，這就為模擬復雜的結構提供了極大的方便。 (3)節(jié)約存儲空間和計算時間節(jié)約存儲空間和計算時間 在時域有限差分法中每個網格電場和磁場的六個分量及其上一時間步的值是必須存儲的。此外還有描

5、述各網格電磁性質的參數以及吸收邊界條件和連接條件的有關參量，它們一般是空間網格總數N的數倍。 所以，時域有限差分法所需要的存儲空間直接由所需的網格空間決定，與網格總數N成正比。 在計算時，每個網格的電磁場都按同樣的差分格式計算，所以，就所需的主要計算時間而言，也是與網格總數N成正比 相比之下，若離散單元也是N，則矩量法所需的存儲空間與3N的2次方成正比，而所需的CPU時間則與(3N)的2次方至(3N)的3次方成正比 (4)適合并行計算適合并行計算 很多復雜的電磁場問題不能計算往往不是沒有可選用的方法，而是由于計算條件的限制。 當代電子計算機的發(fā)展方向是運用并行處理技術，以進一步提高計算速度。

6、(5)計算程序的通用性計算程序的通用性 由于Maxwell方程是時域有限差分法計算任何問題的數學模型，因而它的基本差分方程對廣泛的問題是不變的。 因此一個基礎的時域有限差分法計算程序，對廣泛的電磁場問題具有通用性，對不同的問題或不同的計算對象只需修改有關部分，而大部分是共同的。 頻率升高? (6)簡單、直觀、容易掌握簡單、直觀、容易掌握 由于時域有限差分法直接從Maxwell方程出發(fā)，不需要任何導出方程，這樣就避免了使用更多的數學工具，使得它成為所有電磁場的計算方法中最簡單的一種。 其次，由于它能直接在時域中模擬電磁波的傳播及其與物體作用的物理過程，所以它又是非常直觀的一種方法。 由于它既簡單

7、又直觀，掌握它就不是件很困難的事情，只要有電磁場電磁場的基本理論知識，不需要數學上的很多準備，就可以學習運用這一方法解決很復雜的電磁場問題。這樣，這一方法很容易得到推廣，并在很廣泛的領域發(fā)揮作用。1.麥克斯韋方程組和麥克斯韋方程組和Yee算法算法在無源、線性、各向同性、非色散介質中，麥克斯韋方程組表示為： （1.1a） （1.1b）*HEHt EHEt 三三.非分裂完美吸收邊界FDTD方法的研究 在三維直角坐標系中，方程1.1a和1.1b可表示成如下微分形式： （1.2a）1yxzHEEtxz（1.2b）1yxzEEHtyx（1.2c）1yxzEHEtzy1yxzxHEHEtyz（1.2d）1

8、yxzyEHHEtzx（1.2e）1yxzzHHEEtxy（1.2d）1.1 1.1 一維一維MaxwellMaxwell方程的方程的FDTDFDTD算法（算法（1 1）一維自由空間Maxwell方程利用一階導數的二階中心差分近似，上面的方程變?yōu)?1 (1.1)yxHEtz 01 (1.2)yxHEtz (1.3) )()(1)()(212102/12/1zkHkHtkEkEnynynxnx )() 1(1)()(2/12/1021211zkEkEtkHkHnxnxnyny(1.4) n 時間步長k 空間步長1.1 1.1 一維一維MaxwellMaxwell方程的方程的FDTDFDTD算法（

9、算法（2 2）電場與磁場分量的空間-時間分布圖2/1nxE2k1kk1k2k2/1nxE2k1kk1k2k2/3k2/1k2/1k2/5k3/ 2k nxH0 xtn c 0023xxtcc 穩(wěn)定條件 時間步長t與空間步長x，y和z之間必須滿足一定條件，否則將出現數值不穩(wěn)定性1.1 1.1 一維一維MaxwellMaxwell方程的方程的FDTDFDTD算法（算法（3 3）001/c 1/21/211220()()( )( )1nnnnyyxxHkHkEkEktx 11/21/211220()()(1)( )1 nnnnyyxxHkHkEkEktx 1/21/211220( )( )()()n

10、nnnxxyytEkEkHkHkx11/21/211220()()(1)( ) nnnnyyxxtHkHkEkEkxEE00由于介電常數和磁導率不在一個數量級1/21/2112200( )( )()()nnnnxxyytEkEkHkHkx 11/21/2112200()()(1)( ) nnnnyyxxtHkHkEkEkx 1.1 1.1 一維一維MaxwellMaxwell方程的方程的FDTDFDTD算法（算法（3 3） 用計算機語言表示的FDTD公式式中，時間變量已隱含在迭代公式中, 只要給定了所有空間點上電和磁場的初值，就可以一步一步地求出任意時刻所有空間點上的電場和磁場值。 1*5 .

11、 0kExkExkHykHy 0.5* 1Ex kEx kHy kHy k(1.7) (1.8) 1/21/21122( )( )0.5 ()()nnnnxxyyEkEkHkHk11/21/21122()()0.5 (1)( ) nnnnyyxxHkHkEkEk一維邊界條件 102nxnxEE0023xxtcc 20 xtc距離穩(wěn)定性條件電磁波傳播距離設置的邊界條件一維電磁波在介質中傳播1/21/2112200( )( )()()nnnnxxyyrtEkEkHkHkx 11/21/2112200()()(1)( ) nnnnyyxxtHkHkEkEkx 1*5 . 0kExkExkHykHy

12、* 1Ex kEx kCb kHy kHy k 0.5/Cb kepsilon (1.20a)EHJt01 (1.20b)HEt JE在有耗媒質中001rrEHEt 00( )( )1( )yxxrrHtE tE ttz 000( )( )1( )yxxrrHtE tE ttz 00( )( )1yxHtE ttz 1/21/21/21/2000(1/ 2)(1/ 2)( )( )( )( )12nnnnnnyyxxxxrrHkHkEkEkEkEktx 我們對方程進行時域和空間的近似一維電磁波在有耗介質中傳播一維電磁波在有耗介質中傳播00112tx 1/21/2001/ 2( )1( )1(1

13、/ 2)(1/ 2)22nnnnxxyyrrrttEkEkHkHk 1/21/21/21/2000(1/ 2)(1/ 2)( )( )( )( )12nnnnnnyyxxxxrrHkHkEkEkEkEktx 1/21/21/21/20( )( )1( )( )(1/ 2)(1/ 2)22nnnnnnxxxxyyrrEkEktEkEkHkHk 1/21/2000(1)21/ 2( )( )(1/ 2)(1/ 2)(1)(1)22nnnnrxxyyrrrtEkEkHkHktt 因此我們得到計算機的計算公式 *(1 )ex kca kex kcb khy khy k 0.5*( 1)hy khy k

14、ex kex k一維電磁波在有耗介質中傳播一維電磁波在有耗介質中傳播1/21/2000(1)21/ 2( )( )(1/ 2)(1/ 2)(1)(1)22nnnnrxxyyrrrtEkEkHkHktt 因此我們得到計算機的計算公式 *(1 )ex kca kex kcb khy khy k 0.5*( 1)hy khy kex kex k*/(2*)eafdt sigmaepsz epsilon (1.)/(1.)ca keafeaf 0.5/(*(1.)cb kepsiloneaf一維電磁波在有耗介質中傳播一維電磁波在有耗介質中傳播對于任意媒質，麥克斯韋方程組的微分形式為 )()()(*0E

15、Dr1.2 1.2 一維有耗媒質的一維有耗媒質的FDTDFDTD算法算法(2.1a) (2.1b) (2.1c) （2.1b）寫成了頻域形式DHt 01HEtEE00DD001HtD001)()()(*EDrEtH100令這樣和我們上一節(jié)提到的公式是一樣的 ，只是D變成了E。1.2 1.2 一維有耗媒質的一維有耗媒質的FDTDFDTD算法算法(2.2a) (2.2b) (2.3a) (2.3b) (2.3c) *0()rrj 對于有耗媒質 ,我們可表示為下式0()()()rDEEj00( )( )( )trD tE tE tdt從傅立葉變換理論，我們知道頻域1/j是時域的積分，由于我們是對頻域

16、的采樣，上式變?yōu)?0nnnirtDEE1.2 1.2 一維有耗媒質的一維有耗媒質的FDTDFDTD算法算法(2.4) (2.5) 2.4代入2.3b(2.6) 0crj，我們分離電場,得到公式1000nnnnirttDEEE(2.7) 1000nninrtDEEt(2.8a) E n是當前的電壓值和當前的D值和以前E的和我們令00nnitIE(2.8b) 公式2.8可簡化為0nnnrDIEt10nnntIIE(2.8c) 1.2 1.2 一維有耗媒質的一維有耗媒質的FDTDFDTD算法算法現在,FDTD公式可以表示為:dxk=dxk+0.5*(hyk-1-hyk)exk=gaxk *(dxk-

17、ixk)ixk=ixk+gbxk*exkhyk=hyk+0.5*(exk-exk+1)其中,gaxk=1/(epsilon+(sigma*dt/epsz)gbxk=sigma*dt/epsz所有記錄媒質信息都包含在公式2.9b,2.9c中.對于自由空間,gax=1 和gbx=0；對于有耗媒質，gax和gbx由公式2.10計算。公式2.9a和2.9d，其中包含空間差分，并不隨媒質而改變。1.2 1.2 一維有耗媒質的一維有耗媒質的FDTDFDTD算法算法(2.9a) (2.9b) (2.9c) (2.9d) (2.10a) (2.10b) 0tIDErnnn10nnntIIE0( )()( )r

18、DEj我們可以避免在時間域上處理棘手的卷積積分, 我們可以利用Z變換,10( )( )( )( )rtD zE zz I zE z01/( )( )( )1rtD zE zE zz1010/( )( )( )( )1ttI zE zz I zE zz10( )( )( )rD zz I zE zt10nnnrDIEt10nnntIIE1.3 1.3 利用利用Z Z變換公式變換公式HtD001)()()(*EDrEtH100(2.3a) (2.3b) (2.3c) 1.3 1.3 二維二維MaxwellMaxwell方程的方程的FDTDFDTD算法算法 (1)(1)對于二維仿真，可以用TE和TM

19、波表示。對于TM波對于二維仿真，可以用TE和TM波表示。對于TM波，公式2.3可以簡化為下式)(100yHxHtDxyz)()()(*zrzEDyEtHzx001xEtHzy001，其中02 cxt1.3 1.3 二維二維MaxwellMaxwell方程的方程的FDTDFDTD算法算法 (1)(1)(3.1a) (3.1b) (3.1c) (3.1d) 我們可以得到公式： dz(i,j)=dz(i,j)+0.5*(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1);ez(i,j)=gaz(i,j)*(dz(i,j)-iz(i,j); iz (i,j)=iz(i,j)+gbz

20、(i,j)*ez(I,j)hx(i,j)= hx(i,j)+0.5*( ez(i,j)-ez(i,j+1);hy(i,j)= hy(i,j)+0.5*( ez(i+1,j)-ez(i,j);xH和zD與一維的簡單的有損介質表達方式是一樣的 1.3 1.3 二維二維MaxwellMaxwell方程的方程的FDTDFDTD算法算法 yH(3.2a) (3.2a) (3.2c) (3.2d) (3.2e) 基本理論：如果一束在媒質A中的波，接觸的媒質B，反射系數的大小通過兩種媒質的阻抗描述通常我們假設磁導率為常數，所以，當波照射從 到 時將會產生反射。但當 隨介電常數改變時， 將保持為常數， 為零，

21、沒有反射發(fā)生。但這并不能解決我們的問題，進入媒質的波會繼續(xù)傳播。這就需要一種有耗媒質，波在到達邊界之前，在這種媒質中很快衰減。我們可以通過設置磁導率和介電常數為復數，因為虛數部分代表著衰減。 BABA1.4 PML1.4 PML吸收邊界吸收邊界 （1 1）121.4 PML1.4 PML吸收邊界吸收邊界 （2 2）)(0yHxHcDjxyz)()()(*zrzEDyEcHjzx0 xEcHjzy0我們把公式3.1轉化為復數形式的麥克斯韋方程(3.3a) (3.3b) (3.3c) (3.3d) 對于3.3 式，我們引入有耗媒質，其介電常數和磁導率定義為 *Fz*Fx*Fy)()()(0*yHx

22、HcyxDjxyFzFzz)()()(*zrzEDyEcyxHjzFyFyx0*)()(xEcyxHjzFyFyy0*)()(（1） 只與電位移D有關，與E無關；（2）我們在PML中增加的有耗媒質與真實的介電 常數沒有關系。 1.4 PML1.4 PML吸收邊界吸收邊界 （3 3）(3.4) (3.5a) (3.5b) (3.5c) (3.5d) F(引入的有耗媒質只對吸收邊界起作用） 1.4.1 實現PML的兩個基本條件 (1)Sack提出的構造PML的兩個基本條件1、對于PML來說，媒質1的阻抗應該是一常數，1*0FxFx阻抗為1是歸一化 2、垂直邊界的方向，介電常數和相對磁導率的值必須是

23、相反的 *1FyFx*1FyFx0*mDFmFm0*mHFmFm我們假設介電常數和相對磁導率復數形式為：下面的公式滿足3.7式，其中1FmFm000DmHmD(3.6) (3.7a) (3.7b) (3.8a) (3.8b) (3.9a) (3.9b) 引入的參數1.4.1 實現PML的兩個基本條件 (2)我們把3.9代入3.61/ )(1/ )(100*0jxjxFxFx滿足第一個條件。如果 值隨著進入PML區(qū)而逐漸增大，那么電場和磁場將逐漸衰減。如果我們只在X方向應用PML，公式可以簡化只有X方向的)()(0*yHxHcxDjxyFzzxEcxHjzFyy0*)(*Fx*FxyEcxHjz

24、Fyx0*)(3.10a) (3.10b) (3.10c) 1.4.1 實現PML的兩個基本條件 (3)公式3.9代入3.10)()(1 (00yHxHcDjxjxyzDyEcHjxjzxD010)(1 (xEcHjxjzyD00)(1 (我們注意到公式3.11b和3.11cHx和公式Hy的磁導率相反，滿足第二條件。下面我們把公式3.11a左邊轉化為差分形式：2),(),()(),(),()(2/12/102/12/10jiDjiDitjiDjiDDitDnznzDnznzzDz2)(1 1),(2)(1 1),(02/102/1titjiDtitjiDDnzDnz(3.11a) (3.11b

25、) (3.12) 求Dz與Dz-1的關系我們令02/)(11)(2tiigiD002/)(12/)(1)(3titiigiDD同樣，我們可以處理3.11c，給出下式其中),2/1(5 .0)(2),()(3),(2/12/1jiHigijiDigijiDnynznz)2/1,()2/1,(),2/1(jiHjiHjiHnxnxny),(), 1(5 . 0) 2/ 1( 2),()( 3),(2/12/11jiEjiEifijiHifijiHnznznynz， 02/)2/1(11)2/1(2tiifiD002/)2/1(12/)2/1(1)2/1(3titiifiDD其中 i+1/2，是磁場

26、和電場相差1/2網格1.42 PML公式的推導 (1)(3.12a) (3.12b) (3.13) (3.14) (3.15a) (3.15b) 公式3.11b轉化為下式00( )1zDzxExEjHcyjy ，其中xecurlxjiEjiEyEnznzz_),()1,(2/12/1_)(_)2/1,()2/1,(0001TnDnxnxxecurltxxecurlctjiHjiH) 2/ 1,(2)(_) 2/ 1,() 2/ 1,(2/ 1001jiItxecurlxtcjiHjiHnHxDnxnx2/10 xtc02)()()( 1tiixnifiD在時域代表積分，j/1j在時域代表差分，

27、上式可化為) 2/ 1,()( 1_5 . 0) 2/ 1,() 2/ 1,(2/ 11jiIifiecurljiHjiHnHxnxnx1.42 PML公式的推導 (2)(3.16) (3.17) )()(1 (00yHxHcDjxjxyzDyEcHjxjzxD010)(1 (xEcHjxjzyD00)(1 (1/2( ,)(1 / 2,)( ,1 / 2)( ,1 / 2)nnnnzyxxDijHijHijHij),(), 1(5 . 0) 2/ 1( 2),()( 3),(2/12/11jiEjiEifijiHifijiHnznznynz) 2/ 1,()( 1_5 . 0) 2/ 1,(

28、) 2/ 1,(2/ 11jiIifiecurljiHjiHnHxnxnx最后引入匹配參數后的二維差分方程為了計算f 和g ，我們沒有必要去計算實際的電導率，引入輔助參數，3)_(*333. 0pmllengthixn pmllengthi_,.,2 , 1)()(1ixnifi)(11()(2ixnigi)(1)(1()(3ixnixnigi)(1 ifi)(2 igi)(3 igi這些輔助參數有一定的變化范圍0 0.333 1 0.75 1 0.5對于主要問題空間， =0， =1， =1。因此從主要問題空間到 PML空間是無縫轉換。)( 1 ifi)(2 igi)( 3igi1.42 PM

29、L公式的推導 (3)(3.18) (3.19a) (3.19b) (3.19c) PML區(qū)折疊部分12 ,3 , 2 , 3fjfj fj gj g j減小，增大 PML 邊界參數關系 以上我們應用PML在方向。顯然，在把y方向加到其中。)()(1)()(1 (000yHxHcDjyjxjxyzDDyEcHjyjxjzzDD0010)(1 ()(1 (xEcHjyjxjzyDD0100)(1)()(1 (和前面的處理過程相同，使用下面公式代替3.13), 2/1(5 . 0)( 2)( 2),()( 3)( 3),(2/12/1jiHjgjigijiDjgjigijiDnynznz)2/1,(

30、)2/1,(), 2/1(jiHjiHjiHnxnxny1.42 PML公式的推導 (4)(3.20a) (3.20b) (3.20c) (3.21) 1.42 PML公式的推導 (5)， 在Y方向的可給出如下式), 2/1()( 1_5 . 0)2/1(2), 2/1()2/1(3), 2/1(2/11jiIjfjecurlifijiHifijiHnHynynyecurljiIjiInHxnHx_)2/1,()2/1,(2/12/1),() 1,(_2/12/1jiEjiEecurlnznz在X方向的可給出如下式 ecurljiIjiInHxnHx_)2/1,()2/1,(2/12/1)2/

31、1,()( 3)2/1,(1jiHjfjjiHnxnx(3.22a) (3.22b) (3.22c) )2/1,()( 1_5 . 0)(22/1jiIifiecurljfjnHx(3.23a) (3.23b) (3.23c) 1.5 1.5 三維三維MaxwellMaxwell方程的方程的FDTDFDTD算法（算法（1 1） Yee 首先將空間按立方體分割，電磁場的六個分量在空間的取樣點分別放在立方體的邊沿和表面中心點上，電場與磁場分量在任何方向始終相差半個網格步長。Ey( i, j, k )xyzEzEyEyEzEzHxExExExHzHy下面是三維麥克斯韋方程 、 )(100zHyHtD

32、yzx)(100 xHzHtDzxy)(100yHxHtDxyz)(100yEzEtHzyx)(100zExEtHxzy)(100 xEyEtHyxz對于三維，和二維處理PML邊界完全相同，只是我們處理三個方向公式3.11變成下式)()(1)()(1)()(1 (01000yHxHcDjzjyjxjxyzzyx)()(1 ()(1)()(1 (0000yHxHjzcDjyjxjxyzzyx1.5 1.5 三維三維MaxwellMaxwell方程的方程的FDTDFDTD算法（算法（1 1）(5.1) (5.2) 比二維麥克斯韋方程 ，多出一項我們定義1.5 1.5 三維三維MaxwellMaxw

33、ell方程的方程的FDTDFDTD算法（算法（1 1）hcurljzchcurlcz_1)(_000)()(1 ()(1)()(1 (0000yHxHjzcDjyjxjxyzzyxhcurljIDz_1)(_()(1)()(1 (0000DzzzyxIzhcurlcDjyjxj公式右邊比二維多出一項(5.5) (5.4) (5.3) 1.6總場散射場方法 總場/散射場方法是基于Maxwell方程的線性特性和下列電磁場的分解 totincscattotincscatEEEHHH式中下標inc表示不存在任何材料時的入射波場。下標scat表示散射波場，最初它們是未知的，它們是入射波與材料相互作用產生的場。下標tot表示總場。根據Maxwell方程的線性特性，無論是入射場、散射場或總場都滿足Maxwell方程，所以FDTD法可以獨立地應用于入射場、散射場和總場。如圖4-3所示，將計算域分成兩個區(qū)：區(qū)域1和區(qū)域2。區(qū)域1中包含了所有散射體。在區(qū)域1中用FDTD法模擬總場，稱為總場區(qū)。區(qū)域2中為自由空間的一部分