2017畢業(yè)論文-時域有限差分法對平面te波的matlab仿真

1、<p>　　時域有限差分法對平面TE波的MATLAB仿真</p><p><b>　　摘 要</b></p><p>　　時域有限差分法是由有限差分法發(fā)展出來的數(shù)值計算方法。自1966年Yee在其論文中首次提出時域有限差分以來，時域有限差分法在電磁研究領(lǐng)域得到了廣泛的應(yīng)用。主要有分析輻射條線、微波器件和導(dǎo)行波結(jié)構(gòu)的研究、散射和雷達截面計算、分析周期結(jié)構(gòu)、電

2、子封裝和電磁兼容的分析、核電磁脈沖的傳播和散射以及在地面的反射及對電纜傳輸線的干擾、微光學(xué)元器件中光的傳播和衍射特性等等。</p><p>　　由于電磁場是以場的形態(tài)存在的物質(zhì)，具有獨特的研究方法，采取重疊的研究方法是其重要的特點，即只有理論分析、測量、計算機模擬的結(jié)果相互佐證，才可以認為是獲得了正確可信的結(jié)論。時域有限差分法就是實現(xiàn)直接對電磁工程問題進行計算機模擬的基本方法。在近年的研究電磁問題中，許多學(xué)者對時

3、域脈沖源的傳播和響應(yīng)進行了大量的研究，主要是描述物體在瞬態(tài)電磁源作用下的理論。另外，對于物體的電特性，理論上具有幾乎所有的頻率成分，但實際上，只有有限的頻帶內(nèi)的頻率成分在區(qū)主要作用。</p><p>　　文中主要談到了關(guān)于高斯制下完全匹配層的差分公式的問題，通過MATLAB程序?qū)E波進行了仿真，模擬了高斯制下完全匹配層中磁場分量瞬態(tài)分布。得到了相應(yīng)的磁場幅值效果圖。</p><p>　　

4、關(guān)鍵詞：時域有限差分 完全匹配層 MATLAB 磁場幅值效果圖</p><p><b>　　目 錄</b></p><p><b>　　摘 要1</b></p><p><b>　　目 錄2</b></p><p>　　第一章 緒 論3</p>

5、<p>　　1.1 課題背景與意義3</p><p>　　1.2 時域有限差分法的發(fā)展與應(yīng)用3</p><p>　　2.1 Maxwell方程和Yee氏算法6</p><p>　　2.2 FDTD的基本差分方程8</p><p>　　2.3 時域有限差分法相關(guān)技術(shù)10</p><p>　　2.3

6、.1 數(shù)值穩(wěn)定性問題10</p><p>　　2.3.2 數(shù)值色散11</p><p>　　2.3.3 離散網(wǎng)格的確定12</p><p>　　2.4 吸收邊界條件12</p><p>　　2.4.1 一階和二階近似吸收邊界條件13</p><p>　　2.4.2 二維棱邊及角頂點的處理16</p&g

7、t;<p>　　2.4.3 完全匹配層18</p><p>　　2.5 FDTD計算所需時間步的估計22</p><p>　　第三章 MATLAB的仿真的程序及模擬24</p><p>　　3.1 MATLAB程序及相應(yīng)說明24</p><p>　　3.2 出圖及結(jié)果27</p><p>　　

8、3.2.1程序部分27</p><p>　　3.2.2 所出的效果圖28</p><p>　　第四章 結(jié) 論30</p><p><b>　　參考文獻31</b></p><p><b>　　第一章 緒 論</b></p><p>　　1.1 課題背景與意義&

9、lt;/p><p>　　20世紀60年代以來，隨著計算機技術(shù)的發(fā)展，一些電磁場的數(shù)值計算方法逐步發(fā)展起來，并得到廣泛應(yīng)用，其中主要有：屬于頻域技術(shù)的有限元法（FEM）、矩量法(MM)和單矩法等；屬于時域技術(shù)方面的時域有限差分法(FDTD)、傳輸線矩陣法(TLM)和時域積分方程法等。此外，還有屬于高頻技術(shù)的幾何衍射理論(GTD)和衍射物理理論(PLD)等。各種方法都具有自己的特點和局限性，在實際中經(jīng)常把它們相互配合而形

10、成各種混合方法[1~2]。其中FDTD是一種已經(jīng)獲得廣泛應(yīng)用并且有很大發(fā)展前景的時域數(shù)值計算方法。時域有限差分（FDTD）方法于1966年由K.S.Yee[3] 提出并迅速發(fā)展，且獲得廣泛應(yīng)用。K.S.Yee用后來被稱作Yee氏網(wǎng)格的空間離散方式，把含時間變量的Maxwell旋度方程轉(zhuǎn)化為差分方程，并成功地模擬了電磁脈沖與理想導(dǎo)體作用的時域響應(yīng)。但是由于當(dāng)時理論的不成熟和計算機軟硬件條件的限制，該方法并未得到相應(yīng)的發(fā)展。20世紀80年代

11、中期以后，隨著上述兩個條件限制的逐步解除，F(xiàn)DTD便憑借其特有的優(yōu)勢得以迅速發(fā)展。它能方便、精確地預(yù)測實際工程中的大量復(fù)雜電磁問題，應(yīng)用范圍幾乎涉及所有電磁領(lǐng)域，成</p><p>　　另外，利用矩量法求解電磁場問題時，要用到并失Green函數(shù)。對于某些問題，可以找到其解析形式的并失Green函數(shù)；而對于復(fù)雜的問題，很難找到其解析形式的并失Green函數(shù)，這樣就使得問題無法解決。作為時域分析中的一個重要數(shù)值方法，

12、FDTD不存在這樣的問題。</p><p>　　1.2 時域有限差分法的發(fā)展與應(yīng)用</p><p>　　經(jīng)過四十多年的發(fā)展，F(xiàn)DTD已發(fā)展成為一種成熟的數(shù)值計算方法。在發(fā)展過程中，幾乎都是圍繞幾個重要問題展開的，即數(shù)值穩(wěn)定性、計算精度、數(shù)值色散、激勵源技術(shù)以及開域電磁問題的吸收邊界條件等。</p><p>　　數(shù)值穩(wěn)定和計算精度對任何一種數(shù)值計算方法都是至關(guān)重要的。

13、A.Taylor和M.E.Brodwin[4]利用本征值方法給出了直角坐標系下FDTD的空間步長與時間步長之間的關(guān)系。X.Min等[5]研究了存在邊界條件時FDTD的穩(wěn)定性問題。對于數(shù)值色散，與實際的物理色散不同，它是由電磁場量在空間和時間上的對波動方程作差分近似處理造成的。這種色散引起的誤差造成在計算區(qū)域內(nèi)傳播的電磁波逐漸畸變[6~7]。K. L. Shlager 等[8]比較了二維和三維空間中幾種正交網(wǎng)格算法的色散誤差。當(dāng)采用其他變

14、形或非正交網(wǎng)格時，必須重新分析其數(shù)值穩(wěn)定性和色散特性[9~11]，P.Monk 和 E.Suli[12]分析了不均勻長方體網(wǎng)格算法的穩(wěn)定性。</p><p>　　激勵源的設(shè)計和引入也是FDTD的一個重要任務(wù)。目前，應(yīng)用最廣泛的激勵源引入技術(shù)是總場/散射場體系[12]。對于散射問題，通常在FDTD計算空間中引入連接邊界，它將整個計算空間劃分為內(nèi)部的總場區(qū)和外部的散射場區(qū)，如圖1-1。利用Huygens原理，可以在連

15、接邊界處引入入射場，使入射場的加入變得簡單易行。</p><p><b>　　圖1-1</b></p><p>　　開域電磁問題中，為了在有限的計算空間內(nèi)模擬無限空間中的電磁問題，必須在計算空間的截斷邊界處設(shè)置吸收邊界條件。吸收邊界條件從開始簡單的插值邊界，已經(jīng)發(fā)展了多種吸收邊界條件。在早期得到廣泛應(yīng)用的是G.Mur[13]的一階和二階吸收邊界條件，它是基于B.Eng

16、quist和A.Majda[14]的單向波方程而提出的差分格式，在FDTD仿真區(qū)域外邊界具有0.5%到5%的反射系數(shù)。目前應(yīng)用最廣泛的是J.P.Berenger[15-17]的分裂式完全匹配層，以及Z.S.Sacks等[18]和S.D.Gedney[20]的各向異性介質(zhì)的完全匹配層，它們可使FDTD模擬的最大動態(tài)范圍達到80dB。</p><p>　　另一方面，為了更好的擬合研究對象的形狀，克服臺階逼近帶來的誤差

17、，D.E.Merewether[19]提出了柱坐標系下的網(wǎng)格剖分方法，R.Holland[20]提出了球坐標系下的網(wǎng)格剖分方法，P.Monk和E.Suli[12]提出了變網(wǎng)格步長方法，S.S.Zivanovic等[21]和P.Thoma等[22]提出了亞網(wǎng)格技術(shù)（即在一般區(qū)域采用粗網(wǎng)格，在電磁場快變區(qū)域采用精細網(wǎng)格）。利用這些技術(shù)，可以更精確地模擬各種復(fù)雜的結(jié)構(gòu)，適應(yīng)各種復(fù)雜的介質(zhì)，提高了復(fù)雜介質(zhì)中數(shù)值計算的精度。</p>

18、<p>　　時域模擬一般獲得的是近場電磁信息，為了得到諸如天線方向圖或散射體雷達散射截面之類的遠場信息，必須獲得計算區(qū)域以外的頻域場或瞬態(tài)場。多位學(xué)者在這方面做了許多工作，發(fā)展了一種高效的時域近遠場變換方法[23-26]。借助這種方法，可以實現(xiàn)由計算區(qū)域內(nèi)近場數(shù)據(jù)到計算區(qū)域外遠場數(shù)據(jù)的外推。目前，粗糙面散射的FDTD，傳遞函數(shù)在FDTD中的應(yīng)用，周期介質(zhì)、各向異性介質(zhì)、色散介質(zhì)和含有集中元件的FDTD，以及網(wǎng)絡(luò)并行FDTD技

19、術(shù)等方面也取得了很大進展。</p><p>　　FDTD在迅速發(fā)展的同時，也獲得了非常廣泛的應(yīng)用。目前，它幾乎被應(yīng)用到了電磁場工程中的各個方面，例如：電磁散射、生物電磁計量學(xué)、輻射天線的分析、微波器件和導(dǎo)行波結(jié)構(gòu)的研究、散射和雷達截面的計算、周期結(jié)構(gòu)的分析、電子封裝和電磁兼容的分析、核電磁脈沖傳播和散射的分析、以及微光學(xué)元器件中光的傳播和衍射特性的分析等。隨著新技術(shù)的不斷提出，其應(yīng)用范圍和成效正在迅速地擴大和提高

20、。</p><p>　　第二章 時域有限差分法的基本原理</p><p>　　Maxwell方程是描述宏觀電磁現(xiàn)象的一組基本方程。這組方程即可以寫成微分形式，又可以寫成積分形式。FDTD方法由Maxwell旋度方程的微分形式出發(fā)，利用二階精度的中心差分近似，直接將微分運算轉(zhuǎn)換為差分運算，這樣達到了在一定體積內(nèi)和一段時間上對連續(xù)電磁場數(shù)據(jù)的抽樣壓縮。</p><p>

21、;　　2.1 Maxwell方程和Yee氏算法</p><p>　　根據(jù)[27]中電磁場基本方程組的微分形式，若在無源空間，其空間中的媒質(zhì)是各向同性、線性和均勻的，即媒質(zhì)的參數(shù)不隨時間變化且各向同性，則Maxwell旋度方程可寫成：</p><p><b>　?。?-1a）</b></p><p><b>　　（2-1b）</b

22、></p><p>　　式中，是電場強度，單位為伏/米（V/m）；是磁場強度，單位為安/米（A/m）；表示介質(zhì)介電系數(shù)，單位為法拉/米（F/m）； 表示磁導(dǎo)系數(shù)，單位為亨利/米（H/m）；表示介質(zhì)電導(dǎo)率，單位為西門子/米（S/m）；表示導(dǎo)磁率，單位為歐姆/米（）。</p><p>　　在直角坐標系中，（2-1）式可化為如下六個標量方程：</p><p><

23、;b>　?。?-2）</b></p><p><b>　?。?-3）</b></p><p>　　這六個偏微分方程是FDTD算法的基礎(chǔ)。</p><p>　　K.S.Yee[3]在1966年建立了如圖2-1所示的空間網(wǎng)格，這就是著名的Yee氏元胞網(wǎng)格。</p><p>　　圖2-1 Yee氏網(wǎng)格及其電磁

24、場分量分布</p><p>　　并引入如下的差分近似方法對（2-2）、(2-3)式中的六個偏微分方程進行了差分離散。令代表或在直角坐標系中某一分量，在時間和空間域中的離散可記為</p><p><b>　　(2-4)</b></p><p>　　式中，、和分別是長方體網(wǎng)格沿x、y、z方向的空間步長，是時間步長，i、j、</p>&

25、lt;p>　　k分別是沿x、y、z方向的網(wǎng)格編號，n是時間步數(shù)。對關(guān)于時間和空間的一階偏導(dǎo)數(shù)取中心差分近似，具有二階精度，即</p><p><b>　　(2-5a) </b></p><p><b>　　(2-5b)</b></p><p><b>　　(2-5c)</b></p&g

26、t;<p><b>　　(2-5d)</b></p><p>　　在FDTD中，空間上連續(xù)分布的電磁場物理量離散的空間排布如圖2-1所示。由圖可見，電場和磁場分量在空間交叉放置，使得在每個坐標平面上每個電場分量被磁場環(huán)繞，每個磁場分量也被電場環(huán)繞。這種電磁場的空間結(jié)構(gòu)與電磁感應(yīng)和電磁波傳播的規(guī)律相符，在每一個網(wǎng)格單元都能滿足法拉第感應(yīng)定律和安培環(huán)流定律。各分量的空間相對位置也適

27、合于Maxwell方程的差分計算，能夠恰當(dāng)?shù)孛枋鲭姶艌龅膫鞑ヌ匦?。同時，電場和磁場在時間上交替抽樣，抽樣時間間隔相差半個時間步，使Maxwell旋度方程離散以后構(gòu)成顯式差分方程，從而可以在時間上迭代求解，而不需要進行矩陣求逆運算。因此，由給定相應(yīng)電磁問題的初始條件，F(xiàn)DTD就可以逐步推進地求得以后各個時刻空間電磁場的分布。</p><p>　　2.2 FDTD的基本差分方程</p><p>

28、;　　根據(jù)上述原則，可將（2-2）、（2-3）式離散為如下的差分方程形式：</p><p><b>　?。?-6a）</b></p><p><b>　　(2-6b)</b></p><p><b>　　(2-6c)</b></p><p><b>　　(2-6d)&

29、lt;/b></p><p><b>　　(2-6e)</b></p><p><b>　　(2-6f)</b></p><p><b>　　式中</b></p><p><b>　　， (2-7a)</b></p><p>

30、<b>　　， (2-7b)</b></p><p>　　(2-6)式就是FDTD的基本差分方程組。從式中可以看出，方程組中含有半個空間步和半個時間步，為了便于編程，可將（2-6）式改寫成如下形式[28]：</p><p><b>　　(2-8a)</b></p><p><b>　　(2-8b)</b&g

31、t;</p><p><b>　　(2-8c)</b></p><p><b>　　(2-8d)</b></p><p><b>　　(2-8e)</b></p><p><b>　　(2-8f)</b></p><p>　　根據(jù)上

32、述FDTD差分方程組可得出計算電磁場的時域推進計算方法，如圖2-2所示。</p><p>　　圖2-2 FDTD在時域的交叉半步逐步推進計算</p><p>　　式(2-8a)~(2-8c)的等號左邊的電場值是第n次循環(huán)的電場值，等號右邊的電場值是第n-1次循環(huán)存儲在內(nèi)存中的電場值，磁場值是本次循環(huán)計算得到的磁場值；式(2-8d)~(2-8f) 等號左邊的磁場值是第n次循環(huán)的磁場值，等號

33、右邊的磁場值和電場值都是第n-1次循環(huán)存儲在內(nèi)存中的場值。這樣，就解決了半個時間步在程序中無法表示的問題，而且也沒有破壞電磁場在時間上逐步推進的邏輯關(guān)系。</p><p>　　2.3 時域有限差分法相關(guān)技術(shù)</p><p>　　2.3.1 數(shù)值穩(wěn)定性問題</p><p>　　上述FDTD方程是一種顯式差分方程，在執(zhí)行時，存在一個重要的問題：即算法的穩(wěn)定性問題。這種不

34、穩(wěn)定性表現(xiàn)為在解顯式方程時，隨著時間步數(shù)的繼續(xù)增加，計算結(jié)果也將無限制地增加。Taflove等[4]于1975年對Yee氏差分格式的穩(wěn)定性進行了討論，并導(dǎo)出了對時間步長的限制條件。數(shù)值解是否穩(wěn)定主要取決于時間步長與空間步長、、的關(guān)系。對于非均勻媒質(zhì)構(gòu)成的計算空間選用如下的穩(wěn)定性條件：</p><p><b>　　(2-9)</b></p><p>　　(2-9)式是空

35、間和時間離散之間應(yīng)當(dāng)滿足的關(guān)系，又稱為Courant穩(wěn)定性條件。若采用均勻立方體網(wǎng)格：</p><p><b>　　(2-10)</b></p><p>　　其中，為計算空間中的電磁波的最大速度。</p><p>　　2.3.2 數(shù)值色散</p><p>　　FDTD方程組是對Maxwell旋度方程進行差分近似，在進行數(shù)

36、值計算時，將會在計算網(wǎng)格中引起數(shù)字波模的色散，即在FDTD網(wǎng)格中，電磁波的相速與頻率有關(guān)，電磁波的相速度隨波長、傳播方向及變量離散化的情況不同而改變。這種關(guān)系由非物理因素引起，且色散將導(dǎo)致非物理因素引起的脈沖波形畸變、人為的各向異性和虛假折射等現(xiàn)象。顯然，色散與空間、時間的離散間隔有關(guān)，如下式所示：</p><p><b>　　(2-11)</b></p><p>　

37、　式(2-11)是三維情況下在FDTD方法中的單色平面波數(shù)值色散關(guān)系的一般形式，它表明FDTD計算中波的傳播速度與傳播方向有關(guān)。式中、、分別是波矢量沿、、方向的分量，是角頻率，是被模擬的均勻介質(zhì)中的光速。與數(shù)值色散關(guān)系相對應(yīng)，在無耗介質(zhì)中的單色平面波，色散解析關(guān)系是：</p><p><b>　　(2-12)</b></p><p>　　由式(2-11)可知，當(dāng)式(2

38、-11)中的、、、均趨于零時，它就趨于式(2-12)。也就是說數(shù)值色散是由于用近似差分替代連續(xù)微分而引起的，而且在理論上可以減小到任意程度，只要此時時間步長和空間步長都足夠小，但這將大大增加所需的計算機存儲空間和計算時間，并使累積誤差增加。因此，在實際計算中要根據(jù)問題的性質(zhì)和計算機的軟硬件條件來選擇合適的時間步長和空間步長。為獲得理想的色散關(guān)系，問題空間分割應(yīng)按照小于正常網(wǎng)格的原則進行。一般選取的最大空間步長為，為所研究范圍內(nèi)電磁波的最

39、小波長。由上分析說明，數(shù)值色散在用FDTD法分析電磁場傳播中的影響是不可能避免的，但我們可以盡可能的減小數(shù)值色散的影響。</p><p>　　2.3.3 離散網(wǎng)格的確定</p><p>　　無論是簡單目標還是復(fù)雜目標，在進行FDTD離散時網(wǎng)格尺寸的確定，除了受計算資源的限制不可能取得很小外，還需要考慮以下幾個因素：</p><p>　　1．目標離散精確度的要求。網(wǎng)格

40、應(yīng)當(dāng)足夠小以便能精確模擬目標幾何形狀和電磁參數(shù)。</p><p>　　2．FDTD方法本身的要求。主要是考慮色散誤差的影響。設(shè)網(wǎng)格為立方體，所關(guān)心頻段的頻率上限為，對應(yīng)波長為，則考慮FDTD的數(shù)值色散要求</p><p><b>　?。?-13）</b></p><p>　　通常。上式是根據(jù)已知所關(guān)心頻率上限情況下來確定FDTD網(wǎng)格尺寸的；反之

41、，若給定，則FDTD計算結(jié)果可用的上限頻率也隨之確定。</p><p>　　3．入射波的要求。入射波的上限截止頻率應(yīng)包含所關(guān)心頻率范圍，即。</p><p>　　2.4 吸收邊界條件</p><p>　　由時域有限差分法的基本原理可知，在利用時域有限差分法研究電磁場時，需在全部問題空間建立Yee氏網(wǎng)格空間，并存儲每個單元網(wǎng)格上任一時間步的六個場分量用于下一時間步的計

42、算。而在對于輻射、散射這類開放系統(tǒng)的實際研究中，不可能有無限大的存儲空間。因此，必須在某處將網(wǎng)格空間截斷，且在截斷邊界網(wǎng)格點處運用特殊的場分量計算方法，使得向邊界面行進的波在邊界處保持“外向行進”特性、無明顯的反射現(xiàn)象，并且不會使內(nèi)部空間的場產(chǎn)生畸變，從而用有限網(wǎng)格空間模擬電磁波在無界空間中傳播的情況。具有這種功能的邊界條件稱之為吸收邊界條件，或輻射邊界條件，或網(wǎng)格截斷條件[29~31]，如圖2-3所示。</p><

43、p>　　圖 2-3 附加截斷邊界使計算區(qū)域變?yōu)橛邢抻?lt;/p><p>　　從FDTD的基本差分方程組可以看出，在截斷邊界面上切向場分量的計算需要利用計算空間以外的電磁場分量，因此FDTD基本差分方程對這些截斷邊界面上的場分量失效。如何處理截斷邊界上的場分量，使之與需要考慮的無限空間有盡量小的差異，是FDTD中必須很好解決的一個重要問題。實際上，這是要求在誤差可容忍的范圍內(nèi)，計算空間中的外向波能夠順利通過截斷

44、邊界面而不引起波的明顯反射，使有限計算空間的數(shù)值模擬與實際情況趨于一致，對外向波而言，就像在無限大空間中傳播一樣。所以，需要一種截斷邊界網(wǎng)格處的特殊計算方法，它不僅要保證邊界場計算的必要精度，而且還要大大消除非物理因素引起的波反射，使得用有限的網(wǎng)格空間就能模擬電磁波在無限空間中的傳播。但是如果處理不當(dāng)，截斷邊界面可能造成較大反射，構(gòu)成數(shù)值模擬誤差的一部分，甚至可能造成算法不穩(wěn)定。</p><p>　　加于截斷邊界

45、場分量符合上述要求的算法就稱為吸收邊界條件（Absorbing Boundary Conditions）。</p><p>　　2.4.1 一階和二階近似吸收邊界條件</p><p>　　在截斷邊界附近通常沒有激勵源。考慮齊次波動方程</p><p><b>　　(2-14)</b></p><p>　　式中，表示直角坐

46、標系下任意電磁場分量。</p><p>　　B.Engquist和A.Majda[15]利用偏微分算子對式(2-13)作因式分解，并分別取其Taylor級數(shù)展開式中的第一項和前兩項近似，導(dǎo)出了適合直角坐標系下FDTD吸收邊界條件的單向波動方程，這就是Engquist-Majda吸收邊界條件。設(shè)三維長方體FDTD區(qū)域0<x<a，0<y<b，0<z<d，這時有六個截斷邊界，其一階和

47、二階解</p><p>　　析吸收邊界條件的具體形式見表2-1，其中代表任一直角場分量。G. Mur[14]對表中的吸收邊界條件引入了一種簡單有效的差分數(shù)值算法，即對時間和空間的偏微分取二階中心差分近似，將單向波方程離散化，便形成了著名的G.Mur的一階和二階吸收邊界條件，其總體虛假反射在1%～5%。圖2-5給出了反射系數(shù)（反射波與入射波）與入射角的關(guān)系，由圖2-4可看出，僅當(dāng)入射角較小時其反射系數(shù)較小。<

48、/p><p>　　表2-1 三維長方體FDTD區(qū)域的一階和二階吸收邊界條件</p><p>　　圖2-4 近似吸收邊界條件作用后殘留的反射波與入射波之比</p><p>　　根據(jù)Yee氏元胞網(wǎng)格的特點，在FDTD截斷邊界面上只有電場切向分量和磁場法向分量。以界面為例，此界面僅有，，節(jié)點。由于FDTD中的計算式不涉及區(qū)域，即不涉及截斷邊界界面外的節(jié)點。所以，吸收邊界條件將

49、不考慮，而只考慮電場切向分量和。以為例，對一階和二階吸收邊界條件分別有</p><p><b>　　(2-15a)</b></p><p><b>　　(2-15b)</b></p><p>　　將式(2-15b)分別在距離邊界半個空間步長的輔助網(wǎng)格點處及時刻離散，此時的位置在，并對各項做差分近似，同時對輔助網(wǎng)格點處的場值

50、應(yīng)用線性插值，假設(shè)，可得到一階吸收邊界條件在三維情況的形式：</p><p><b>　　(2-16)</b></p><p>　　式中代表截斷邊界上切向場分量。</p><p>　　將式(2-16)所示二階吸收邊界推廣到長方體元胞各邊、和不相等情形。設(shè)邊界為，對于分量有</p><p><b>　　(2-17

51、)</b></p><p>　　同理，可以得到所有截斷邊界面上切向場分量的一階和二階差分式。需要注意的是，用二階近似條件計算界面上與棱邊相鄰的一列節(jié)點時會涉及棱邊上的場值。因此，要想避免用到棱邊上的場值，只需對截斷邊界面上切向場分量的計算按以下兩種情況區(qū)別對待：（1）截斷邊界面上與棱邊相鄰的一列場分量采用一階差分式；（2）截斷邊界面上其它場分量采用二階差分式。這樣，就完全不必考慮棱邊上的場分量，避免了

52、計算棱邊上場分量所帶來的誤差。實際計算表明，這樣做提高了吸收邊界條件的精度和FDTD計算的穩(wěn)定性。但是，就目前FDTD的發(fā)展來看，G.Mur的一階和二階吸收邊界條件已不能滿足目前高精度計算的要求。</p><p>　　2.4.2 二維棱邊及角頂點的處理</p><p>　　對于二維電磁場問題，在用FDTD計算邊界處的元胞時，將涉及到截斷邊界外側(cè)的節(jié)點。如對于二維TM情況，電磁場分量有，由圖

53、2-5可見，在用FDTD計算邊界處的TM元胞和時并不涉及截斷邊界以外或的節(jié)點。只有涉及截斷邊界外側(cè)的節(jié)點。因此，只需給出邊界處切向場分量的吸收邊界條件。同樣，對于TE波只需給出邊界處切向場分量的吸收邊界條件。表2-2為矩形截斷邊界面上節(jié)點的二階Mur吸收邊界條件。</p><p>　　圖2-5 二維TM左截斷邊界元胞</p><p>　　表2-2 矩形截斷邊界四邊上的二階Mur吸收邊界條件

54、</p><p>　　在二維矩形計算區(qū)域的角點，吸收邊界條件的離散式需特殊考慮。假設(shè)角點附近只有向外傳播的行波，且傳播方向沿角點處元胞的對角線，如圖2-6所示的矩形區(qū)域左下角點處的TM元胞為例，導(dǎo)出適用于角點的吸收邊界條件離散式。根據(jù)Courant穩(wěn)定條件式，有。在點與對角點之間取一點，對于沿對角線的外行波，有以下等式：</p><p><b>　?。?-18）</b>

55、;</p><p>　　由于傳播距離很小，上式中略去了振幅的衰減。在利用線性插值</p><p><b>　?。?-19）</b></p><p>　　由式（2-18）解出后代入（2-17）式得</p><p><b>　?。?-20）</b></p><p><b&g

56、t;　　或</b></p><p><b>　　（2-21）</b></p><p>　　對于矩形域的其它角點可作類似的處理。對于二維TE波情況，將式（2-19），（2-21）式中的換為即可。</p><p>　　圖2-6 矩形域四個角點</p><p>　　2.4.3 完全匹配層</p>&l

57、t;p>　　完全匹配層[16~18]（Perfectly Matched Layer, PML）是1994年由J.P.Berenger首先提出，并將其設(shè)置在FDTD計算區(qū)域截斷邊界處，用來吸收外向電磁波。Berenger假設(shè)將電磁場分量在PML介質(zhì)中分裂，并分別對各個分裂的場分量賦以不同的損耗。這就相當(dāng)于在FDTD區(qū)域截斷邊界外設(shè)置了一種特殊的非物理的吸收介質(zhì)層，該層介質(zhì)的波阻抗與相鄰介質(zhì)的波阻抗完全匹配，因而外向波將無反射地穿過

58、分界面進入PML。同時，由于PML為有耗介質(zhì)，而且不依賴于外向波的入射角和波阻抗，即使為有限厚度，外向波在其中也會迅速衰減。在實際計算中，PML是目前一種很常用的吸收邊界條件，有很好的吸收效果，其總的網(wǎng)格噪聲能量是使用普通吸收邊界條件時的1/107，可使FDTD模擬的最大動態(tài)范圍達到80dB。</p><p>　　在PML介質(zhì)中，其網(wǎng)格剖分方式與常規(guī)FDTD網(wǎng)格完全一致，每個場分量在網(wǎng)格中的位置也不變，只是都被分

59、裂為兩個子分量。這樣，式(2-1)可寫成</p><p><b>　　(2-22a)</b></p><p><b>　　(2-22b)</b></p><p><b>　　(2-22c)</b></p><p><b>　　(2-22d)</b><

60、/p><p><b>　　(2-22e)</b></p><p><b>　　(2-22f)</b></p><p><b>　　(2-22g)</b></p><p><b>　　(2-22h)</b></p><p><b&g

61、t;　　(2-22i)</b></p><p><b>　　(2-22j)</b></p><p><b>　　(2-17k)</b></p><p><b>　　(2-22l)</b></p><p>　　式中，，，，，為電導(dǎo)率和磁導(dǎo)率，描述了PML介質(zhì)的各向異性

62、。當(dāng)且時， PML介質(zhì)退化為普通有耗介質(zhì);當(dāng)且時，退化為自由空間中的Maxwell方程。同時，要滿足PML介質(zhì)的重要基本條件即阻抗匹配條件，如下式所示：</p><p><b>　　(2-23)</b></p><p>　　式中和分別為真空的介電常數(shù)和真空的磁導(dǎo)率。</p><p>　　對式(2-18)做差分處理，可得到PML介質(zhì)中的FDTD差

63、分方程表達式：</p><p><b>　　(2-23a)</b></p><p><b>　　(2-23b)</b></p><p><b>　　(2-23c)</b></p><p><b>　　(2-23d)</b></p><p

64、><b>　　(2-23e)</b></p><p><b>　　(2-23f)</b></p><p><b>　　(2-23g)</b></p><p><b>　　(2-23h)</b></p><p><b>　　(2-23i)&l

65、t;/b></p><p><b>　　(2-23j)</b></p><p><b>　　(2-23k)</b></p><p><b>　　(2-23l)</b></p><p><b>　　式中</b></p><p>

66、　　， (2-24)</p><p><b>　　， (2-25)</b></p><p><b>　　其中，；</b></p><p>　　在實際計算中，在FDTD計算中PML的設(shè)置不可能延伸到半無限空間，只能是有限厚度，因此PML的外側(cè)邊界需要特殊處理，通常情況下采用理想導(dǎo)體邊界截斷。這樣，透入PML中的外向波到

67、達理想導(dǎo)體邊界處會反射回來，重新進入計算區(qū)域，PML的反射系數(shù)不再等于零。PML介質(zhì)內(nèi)沿方向的電導(dǎo)率分布通常采取以下函數(shù)形式</p><p><b>　?。?-26）</b></p><p>　　式中，為PML層的厚度，為PML層靠近FDTD分界面的距離，是電導(dǎo)率分布階數(shù)，表示PML中電導(dǎo)率變化的程度，取整數(shù)。式（2-26）說明了實際計算中PML介質(zhì)層厚度必須取若干個

68、空間步長，使電導(dǎo)率從FDTD-PML分界面的0漸變到PML最外面的，避免電導(dǎo)率躍變太大，盡量消除數(shù)值反射。這樣，如果相對于FDTD-PML分界面定義的外向波入射角為，則PML內(nèi)側(cè)表面反射系數(shù)為</p><p><b>　?。?-27）</b></p><p>　　使用PML吸收邊界條件時，首先要選定三個參數(shù)：PML層數(shù)，電導(dǎo)率分布階數(shù)，PML表面反射系數(shù)。大量的數(shù)值實

69、驗表明：（1）當(dāng)層數(shù)N固定時，減小，即增加PML的衰減，可以使局部及總體誤差都單調(diào)地減小。然而，當(dāng)減小到一定程度后，這種現(xiàn)象不再出現(xiàn)，原因是存在由空間網(wǎng)格引起的固有誤差。（2）增加PML層數(shù)，可以使局部及總體誤差都單調(diào)地減小，但是N過大又會使計算量劇增，需折衷考慮吸收效果和計算量。（3）電導(dǎo)率分布階數(shù)的改變不會影響計算量，卻會影響PML的吸收效果，一般情況下，越大，吸收效果越好。因此，在實際計算中參數(shù)的選擇隨所處理問題的不同而不同，需綜

70、合考慮，并由數(shù)值實驗尋找最佳值。</p><p>　　2.5 FDTD計算所需時間步的估計</p><p>　　為了使計算達到穩(wěn)定,通常計算所需要時間步按照電磁波往返穿越FDTD計算區(qū)對角線3～5次來估計。若FDTD計算區(qū)總元胞數(shù)為，則對角線上元胞約為 （三維)。按照Courant穩(wěn)定條件，設(shè)計算中心區(qū)，即穿越對角線一次需要時間步為?？傆嬎銜r間步約需步。對于二維情況則約為?；蛘哒f，計算時間

71、步大約等于FDTD計算區(qū)對角線上元胞數(shù)目的12～20倍。實際上，計算所需時間步還與散射體具體形狀、結(jié)構(gòu)有關(guān)。</p><p>　　圖2-9給出了應(yīng)用FDTD分析電磁場問題時的程序流程圖</p><p>　　圖2-9 FDTD 程序流程圖</p><p>　　第三章 MATLAB的仿真的程序及模擬</p><p>　　3.1 MATLAB程序

72、及相應(yīng)說明</p><p>　　% 二維FDTD TE波仿真</p><p>　　clear all;</p><p><b>　　% 定義常數(shù)</b></p><p>　　pi=3.1415;</p><p>　　c=3.0e10; %高斯制下光速</p

73、><p>　　f=1.0e15; %頻率</p><p>　　lambda=c/f; %波長</p><p>　　nmax=400; %時間步數(shù)</p><p>　　del_s=lambda/20; %每最小波長20個采樣點<

74、;/p><p>　　del_t=0.5*del_s/c; %迭代時間步長</p><p>　　n=182; %真空區(qū)域網(wǎng)格數(shù)</p><p>　　np=9; %pml層數(shù)</p><p>　　N1=n+2*np;

75、%總網(wǎng)格數(shù)</p><p>　　N=N1+1; %采樣點數(shù)</p><p>　　M=4; %導(dǎo)電率漸變指數(shù)</p><p>　　sigma_max=(M+1)/1.50/pi/del_s; %最大導(dǎo)電率</p><p>　　% TE波的分量初始化</p>

76、;<p><b>　　tic;</b></p><p>　　figure(1);</p><p>　　axis([0 N 0 N -0.5 0.5]);</p><p>　　Ex=zeros(N1,N); %x方向為橫向，采樣點為網(wǎng)格的橫向邊，故行數(shù)+1</p><p>　　Ey=zer

77、os(N,N1); %y方向為縱向，采樣點為網(wǎng)格的縱向邊，故列數(shù)+1</p><p>　　Bz=zeros(N1,N1); %矩陣行為縱向網(wǎng)格數(shù)，矩陣列為橫向網(wǎng)格數(shù),循環(huán)中用j表示行數(shù)，i表示列數(shù)</p><p>　　Bzx=zeros(N1,N1);</p><p>　　Bzy=zeros(N1,N1);</p>

78、<p>　　Bzxx=zeros(nmax,2);</p><p>　　%進入電磁場迭代計算</p><p>　　for tt=1:nmax</p><p>　　for i=1:N1</p><p>　　if i>=np+1&&i<=N1-np</p><p><b>

79、　　di=0;</b></p><p>　　elseif i<=np</p><p>　　di=np-i+0.5;</p><p>　　elseif i>=N1-np+1</p><p>　　di=np+i-N1-0.5;</p><p>　　end %

80、di是采樣點橫向距PML內(nèi)邊界的距離</p><p>　　sigma_mx=sigma_max*(di/np)^M;</p><p>　　for j=1:N1</p><p>　　if j>=np+1&&j<=N1-np</p><p><b>　　dj=0;</b></p>&

81、lt;p>　　elseif j<=np</p><p>　　dj=np-j+0.5;</p><p>　　elseif j>=N1-np+1</p><p>　　dj=np+j-N1-0.5;</p><p>　　end %dj是采樣點縱向距PML內(nèi)邊界的距離</p>&

82、lt;p>　　sigma_my=sigma_max*(dj/np)^M;</p><p>　　if sigma_mx+sigma_my==0 %真空區(qū)</p><p>　　if j==100&&i==100</p><p>　　t=30; %可選擇高斯脈沖</p><p>　　

83、term=(tt-t);</p><p>　　pulse=exp(-4*pi*term^2/20^2);</p><p>　　pulse=sin(2*pi*tt/40); %可選正弦時諧源</p><p>　　c_miu=c*del_t/del_s;</p><p>　　Eterm1=c_miu*(Ex(i,j+1)-Ex(i,j));<

84、;/p><p>　　Eterm2=c_miu*(Ey(i+1,j)-Ey(i,j));</p><p>　　Bz(i,j)=Bz(i,j)+Eterm1-Eterm2+pulse;%加入脈沖源</p><p><b>　　else</b></p><p>　　c_miu=c*del_t/del_s;</p>&

85、lt;p>　　Eterm1=c_miu*(Ex(i,j+1)-Ex(i,j));</p><p>　　Eterm2=c_miu*(Ey(i+1,j)-Ey(i,j));</p><p>　　Bz(i,j)=Bz(i,j)+Eterm1-Eterm2;</p><p><b>　　end</b></p><p>　　

86、else %PML區(qū)</p><p>　　cpm=(1-2*c*sigma_mx*del_t)/(1+2*c*sigma_mx*del_t);</p><p>　　cqm=c*del_t/(1+2*c*sigma_mx*del_t)/del_s;</p><p>　　Bzx(i,j)=cpm*Bzx(i,j)-c

87、qm*(Ey(i+1,j)-Ey(i,j));</p><p>　　cpm=(1-2*c*sigma_my*del_t)/(1+2*c*sigma_my*del_t);</p><p>　　cqm=c*del_t/(1+2*c*sigma_my*del_t)/del_s;</p><p>　　Bzy(i,j)=cpm*Bzy(i,j)+cqm*(Ex(i,j+1)-

88、Ex(i,j));</p><p>　　Bz(i,j)=Bzx(i,j)+Bzy(i,j);</p><p><b>　　end</b></p><p><b>　　end</b></p><p><b>　　end</b></p><p>　　for

89、i=2:N1</p><p>　　if i>=np+1&&i<=N1-np</p><p><b>　　di=0;</b></p><p>　　elseif i<=np</p><p>　　di=np-i+1;</p><p>　　elseif i>=N1-

90、np+1</p><p>　　di=np+i-N1-1;</p><p>　　end %di是采樣點橫向距PML內(nèi)邊界的距離</p><p>　　sigma_ex=sigma_max*(di/np)^M;</p><p>　　for j=1:N1</p><p>　　ca

91、m=(1-2*c*sigma_ex*del_t)/(1+2*c*sigma_ex*del_t);</p><p>　　cbm=c*del_t/(1+2*c*sigma_ex*del_t)/del_s;</p><p>　　Ey(i,j)=cam*Ey(i,j)-cbm*(Bz(i,j)-Bz(i-1,j));</p><p><b>　　end</b

92、></p><p><b>　　end</b></p><p>　　for i=1:N1</p><p>　　for j=2:N1</p><p>　　if j>=np+1&&j<=N1-np</p><p><b>　　dj=0;</b>&

93、lt;/p><p>　　elseif j<=np</p><p>　　dj=np-j+1;</p><p>　　elseif j>=N1-np+1</p><p>　　dj=np+j-N1-1;</p><p>　　end %dj是采樣點縱向距PML內(nèi)邊界的距離<

94、/p><p>　　sigma_ey=sigma_max*(dj/np)^M;</p><p>　　cam=(1-2*c*sigma_ey*del_t)/(1+2*c*sigma_ey*del_t);</p><p>　　cbm=c*del_t/(1+2*c*sigma_ey*del_t)/del_s;</p><p>　　Ex(i,j)=cam*

95、Ex(i,j)+cbm*(Bz(i,j)-Bz(i,j-1));</p><p><b>　　end</b></p><p><b>　　end</b></p><p>　　Bzxx(tt,1)=Bz(90,50); %對靠近邊界的中央磁場點采樣</p><p&g

96、t;　　Bzxx(tt,2)=Bz(50,90);</p><p><b>　　3.2 出圖及結(jié)果</b></p><p><b>　　3.2.1程序部分</b></p><p>　　figure(1); %可視化處理</p><p>&

97、lt;b>　　clf;</b></p><p>　　mesh(Bz); %磁場的幅值</p><p>　　axis([0 N 0 N -0.5 0.5]);</p><p>　　xlabel('i')</p><p>　　ylabel('j&

98、#39;)</p><p><b>　　drawnow;</b></p><p><b>　　end</b></p><p>　　figure(2);</p><p>　　plot(Bzxx);</p><p>　　figure(3);</p><p>

99、;　　title('磁場幅值分布圖');</p><p>　　surface(Bz);</p><p>　　shading interp;</p><p>　　axis square;</p><p><b>　　toc</b></p><p>　　3.2.2 所出的效果圖<

100、/p><p>　　正弦時諧場源輻射效果圖</p><p><b>　　高斯脈沖輻射效果圖</b></p><p><b>　　第四章 結(jié) 論</b></p><p>　　通過對時域差分法的學(xué)習(xí)而進行的關(guān)于二維FDTD的TE波仿真，鞏固了迭代式的形式與相應(yīng)的循環(huán)編程技巧。但也存在一定問題，由于在吸收上

101、我的編程思路不是很規(guī)范，所以在吸收效果上與理論上的完全匹配層吸收效果還有差距，但是已經(jīng)基本達到了我的預(yù)期目的?？赡苁俏业姆謪^(qū)方法不對的緣故，原因就是出在sigma_max的設(shè)置上面，在高斯制下長度單位變?yōu)閏m，所以同樣的空間步長在高斯制下數(shù)量級就擴大了100倍，而sigma_max的計算式里面分母含有空間步長，所以如果不加以修正，sigma_max在高斯制下就變成了國際制下的1%，太大的導(dǎo)電率和太小的導(dǎo)電率都不能實現(xiàn)較好的吸收效果，下階

102、段應(yīng)將其做適當(dāng)修正并且寫粒子仿真時用高斯制Maxwell方程可以參考角點反射大，多設(shè)置幾層PML試一試。</p><p><b>　　參考文獻</b></p><p>　　[1]J.M.Jin, M.Zunoubi, K.C.Donepudi, W.C.Chew, Frequency-Domain and Time-Domain Finite-Element Solu

103、tion of Maxwell’s Equations Using Spectral Lanczos Decomposition Method. Comput. Methods Appl. Mech. Engrg. 1999,169:279~296</p><p>　　[2]K.Dongsoo,H.B.Lee. Hybird Full-Wave Analysis of Via-hole Grounds Using

104、 Finit-Difference and Finite-Element Time-Domain Methods. IEEE Transactions on Microwave Theory and Techniques. 1997,45(12):2217~2222</p><p>　　[3] Kane S. Yee Numerical solution of initial boundary value pro

105、blems involving Maxwell’s equations in isotropic media[J].IEEE Trans Antennas Propagation, 1966,AP-14(3):302~307</p><p>　　[4] A.Taylor and M.E.Brodwin. Numerical solution of steady-state EM scattering probl

106、ems using the time-Maxwell’s equations. [J] IEEE Trans.Microwave Theroy Tech.,Aug. 1975，MTT-23:623~630</p><p>　　[5] X.Min,W.Sun and K.M.Chen. Stability analysis of the finite difference time domain method ap

107、plied to unbounded electromagnetic problem [J]. IEEE Antennas and Propagation Society International Symposium,Dallas:May 1990,(4):1640~1643</p><p>　　[6]A.Taflove. Review of the formulation and applications o

108、f the finite-difference time-domain method for numerical modeling of electromagnetic wave interactions with arbitrary structures [J]. Wave Motion,June 1988,10(6): 547~582</p><p>　　[7] I.S.Kim and W.J.R.Hoefe

109、r. Numerical dispersion characteristics and stability factor for the TD-FD method [J].Electronics Letters, July 1990, 26(7): 485~487</p><p>　　[8] K.L.shlager,J.G.Msloney, S.L.Ray,and etc.. Relative accuracy

110、of several finite-difference time-domain methods in two and three dimensions [J]. IEEE Transactions on Antennas and Propagation,Dec.1993,AP-41(12)：1732~1737</p><p>　　[9] S.L.Ray. Numerical dispersion and sta

111、bility characteristics of finite-difference time-domain methods on nonorthogonal meshes [J]. IEEE Transactions on Antennas and Propagation,Feb.1993,AP-41(2)：233~235</p><p>　　[10] H.Shi and J.L.Drewriak. Disp

112、ersive comparison for DSI- and tenor-based nonorthogonal FDTD [J]. IEEE Microwave and Guided Wave Letters, May 1996,6(5):193~195</p><p>　　[11] F.Xiao and H.Yabe. Numerical dispersion relation for FDTD method

113、 in general curvilinear coordinates [J]. IEEE Microwave and Guided Wave Letters, Feb.1997,7(2):48~50</p><p>　　[12] R.Holland and J.Williams．Total field versus scattered field finite difference codes:a compar

114、ative assessment[J].IEEE Transactions on Nuclear Science，Dec.1983，NS-30(6)：4583～4588</p><p>　　[13]G.Mur．Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagneti

115、c field equations[J]．IEEE Transactions on Electromagnetic Compatibility，Nov.1981，EMC-23(4)：377～382</p><p>　　[14]B.Engquist and A.Majda. Absorbing boundary conditions for the numerical simulation of waves [J

116、].Mathematics of the Computation，July 1977，31(139)：629～651</p><p>　　[15]J.P.Berenger．A perfectly matched layer for the absorption of electromagnetic waves[J]．Journal of Computational Physics，Oct.1994，114(2)：

117、185～200</p><p>　　[16]J.P.Berenger．Perfectly matched layer for the FDTD solution of wave-structure interaction problem[J]．IEEE Transactions on Antennas Propagation，Jan.1996，AP-44(1)：110～117</p><p&g