畢業(yè)設(shè)計(jì)（論文）時(shí)域有限差分法對(duì)平面TE波的MATLAB仿真

1、時(shí)域有限差分法對(duì)平面te波的matlab仿真摘 要 時(shí)域有限差分法是由有限差分法發(fā)展出來(lái)的數(shù)值計(jì)算方法。自1966年yee在其論文中首次提出時(shí)域有限差分以來(lái)，時(shí)域有限差分法在電磁研究領(lǐng)域得到了廣泛的應(yīng)用。主要有分析輻射條線、微波器件和導(dǎo)行波結(jié)構(gòu)的研究、散射和雷達(dá)截面計(jì)算、分析周期結(jié)構(gòu)、電子封裝和電磁兼容的分析、核電磁脈沖的傳播和散射以及在地面的反射及對(duì)電纜傳輸線的干擾、微光學(xué)元器件中光的傳播和衍射特性等等。由于電磁場(chǎng)是以場(chǎng)的形態(tài)存在的物質(zhì)，具有獨(dú)特的研究方法，采取重疊的研究方法是其重要的特點(diǎn)，即只有理論分析、測(cè)量、計(jì)算機(jī)模擬的結(jié)果相互佐證，才可以認(rèn)為是獲得了正確可信的結(jié)論。時(shí)域有限差分法就是實(shí)

2、現(xiàn)直接對(duì)電磁工程問(wèn)題進(jìn)行計(jì)算機(jī)模擬的基本方法。在近年的研究電磁問(wèn)題中，許多學(xué)者對(duì)時(shí)域脈沖源的傳播和響應(yīng)進(jìn)行了大量的研究，主要是描述物體在瞬態(tài)電磁源作用下的理論。另外，對(duì)于物體的電特性，理論上具有幾乎所有的頻率成分，但實(shí)際上，只有有限的頻帶內(nèi)的頻率成分在區(qū)主要作用。 文中主要談到了關(guān)于高斯制下完全匹配層的差分公式的問(wèn)題，通過(guò)matlab程序?qū)e波進(jìn)行了仿真，模擬了高斯制下完全匹配層中磁場(chǎng)分量瞬態(tài)分布。得到了相應(yīng)的磁場(chǎng)幅值效果圖。關(guān)鍵詞：時(shí)域有限差分 完全匹配層 matlab 磁場(chǎng)幅值效果圖目 錄摘 要1目 錄2第一章 緒 論31.1 課題背景與意義31.2 時(shí)域有限差分法的發(fā)展與應(yīng)用32.1

3、maxwell方程和yee氏算法62.2 fdtd的基本差分方程82.3 時(shí)域有限差分法相關(guān)技術(shù)102.3.1 數(shù)值穩(wěn)定性問(wèn)題102.3.2 數(shù)值色散112.3.3 離散網(wǎng)格的確定122.4 吸收邊界條件122.4.1 一階和二階近似吸收邊界條件132.4.2 二維棱邊及角頂點(diǎn)的處理162.4.3 完全匹配層182.5 fdtd計(jì)算所需時(shí)間步的估計(jì)22第三章 matlab的仿真的程序及模擬243.1 matlab程序及相應(yīng)說(shuō)明243.2 出圖及結(jié)果273.2.1程序部分273.2.2 所出的效果圖28第四章 結(jié) 論30參考文獻(xiàn)31第一章 緒 論1.1 課題背景與意義20世紀(jì)60年代以來(lái)，隨著計(jì)

4、算機(jī)技術(shù)的發(fā)展，一些電磁場(chǎng)的數(shù)值計(jì)算方法逐步發(fā)展起來(lái)，并得到廣泛應(yīng)用，其中主要有：屬于頻域技術(shù)的有限元法（fem）、矩量法(mm)和單矩法等；屬于時(shí)域技術(shù)方面的時(shí)域有限差分法(fdtd)、傳輸線矩陣法(tlm)和時(shí)域積分方程法等。此外，還有屬于高頻技術(shù)的幾何衍射理論(gtd)和衍射物理理論(pld)等。各種方法都具有自己的特點(diǎn)和局限性，在實(shí)際中經(jīng)常把它們相互配合而形成各種混合方法12。其中fdtd是一種已經(jīng)獲得廣泛應(yīng)用并且有很大發(fā)展前景的時(shí)域數(shù)值計(jì)算方法。時(shí)域有限差分（fdtd）方法于1966年由k.s.yee3 提出并迅速發(fā)展，且獲得廣泛應(yīng)用。k.s.yee用后來(lái)被稱作yee氏網(wǎng)格的空間離散

5、方式，把含時(shí)間變量的maxwell旋度方程轉(zhuǎn)化為差分方程，并成功地模擬了電磁脈沖與理想導(dǎo)體作用的時(shí)域響應(yīng)。但是由于當(dāng)時(shí)理論的不成熟和計(jì)算機(jī)軟硬件條件的限制，該方法并未得到相應(yīng)的發(fā)展。20世紀(jì)80年代中期以后，隨著上述兩個(gè)條件限制的逐步解除，fdtd便憑借其特有的優(yōu)勢(shì)得以迅速發(fā)展。它能方便、精確地預(yù)測(cè)實(shí)際工程中的大量復(fù)雜電磁問(wèn)題，應(yīng)用范圍幾乎涉及所有電磁領(lǐng)域，成為電磁工程界和理論界研究的一個(gè)熱點(diǎn)。目前，fdtd日趨成熟，并成為分析大部分實(shí)際電磁問(wèn)題的首選方法。另外，利用矩量法求解電磁場(chǎng)問(wèn)題時(shí)，要用到并失green函數(shù)。對(duì)于某些問(wèn)題，可以找到其解析形式的并失green函數(shù)；而對(duì)于復(fù)雜的問(wèn)題，很難找

6、到其解析形式的并失green函數(shù)，這樣就使得問(wèn)題無(wú)法解決。作為時(shí)域分析中的一個(gè)重要數(shù)值方法，fdtd不存在這樣的問(wèn)題。1.2 時(shí)域有限差分法的發(fā)展與應(yīng)用經(jīng)過(guò)四十多年的發(fā)展，fdtd已發(fā)展成為一種成熟的數(shù)值計(jì)算方法。在發(fā)展過(guò)程中，幾乎都是圍繞幾個(gè)重要問(wèn)題展開的，即數(shù)值穩(wěn)定性、計(jì)算精度、數(shù)值色散、激勵(lì)源技術(shù)以及開域電磁問(wèn)題的吸收邊界條件等。數(shù)值穩(wěn)定和計(jì)算精度對(duì)任何一種數(shù)值計(jì)算方法都是至關(guān)重要的。a.taylor和m.e.brodwin4利用本征值方法給出了直角坐標(biāo)系下fdtd的空間步長(zhǎng)與時(shí)間步長(zhǎng)之間的關(guān)系。x.min等5研究了存在邊界條件時(shí)fdtd的穩(wěn)定性問(wèn)題。對(duì)于數(shù)值色散，與實(shí)際的物理色散不同，

7、它是由電磁場(chǎng)量在空間和時(shí)間上的對(duì)波動(dòng)方程作差分近似處理造成的。這種色散引起的誤差造成在計(jì)算區(qū)域內(nèi)傳播的電磁波逐漸畸變67。k. l. shlager 等8比較了二維和三維空間中幾種正交網(wǎng)格算法的色散誤差。當(dāng)采用其他變形或非正交網(wǎng)格時(shí)，必須重新分析其數(shù)值穩(wěn)定性和色散特性911，p.monk 和 e.suli12分析了不均勻長(zhǎng)方體網(wǎng)格算法的穩(wěn)定性。激勵(lì)源的設(shè)計(jì)和引入也是fdtd的一個(gè)重要任務(wù)。目前，應(yīng)用最廣泛的激勵(lì)源引入技術(shù)是總場(chǎng)/散射場(chǎng)體系12。對(duì)于散射問(wèn)題，通常在fdtd計(jì)算空間中引入連接邊界，它將整個(gè)計(jì)算空間劃分為內(nèi)部的總場(chǎng)區(qū)和外部的散射場(chǎng)區(qū)，如圖1-1。利用huygens原理，可以在連接邊

8、界處引入入射場(chǎng)，使入射場(chǎng)的加入變得簡(jiǎn)單易行。圖1-1開域電磁問(wèn)題中，為了在有限的計(jì)算空間內(nèi)模擬無(wú)限空間中的電磁問(wèn)題，必須在計(jì)算空間的截?cái)噙吔缣幵O(shè)置吸收邊界條件。吸收邊界條件從開始簡(jiǎn)單的插值邊界，已經(jīng)發(fā)展了多種吸收邊界條件。在早期得到廣泛應(yīng)用的是g.mur13的一階和二階吸收邊界條件，它是基于b.engquist和a.majda14的單向波方程而提出的差分格式，在fdtd仿真區(qū)域外邊界具有0.5%到5%的反射系數(shù)。目前應(yīng)用最廣泛的是j.p.berenger15-17的分裂式完全匹配層，以及z.s.sacks等18和s.d.gedney20的各向異性介質(zhì)的完全匹配層，它們可使fdtd模擬的最大動(dòng)態(tài)

9、范圍達(dá)到80db。另一方面，為了更好的擬合研究對(duì)象的形狀，克服臺(tái)階逼近帶來(lái)的誤差，d.e.merewether19提出了柱坐標(biāo)系下的網(wǎng)格剖分方法，r.holland20提出了球坐標(biāo)系下的網(wǎng)格剖分方法，p.monk和e.suli12提出了變網(wǎng)格步長(zhǎng)方法，s.s.zivanovic等21和p.thoma等22提出了亞網(wǎng)格技術(shù)（即在一般區(qū)域采用粗網(wǎng)格，在電磁場(chǎng)快變區(qū)域采用精細(xì)網(wǎng)格）。利用這些技術(shù)，可以更精確地模擬各種復(fù)雜的結(jié)構(gòu)，適應(yīng)各種復(fù)雜的介質(zhì)，提高了復(fù)雜介質(zhì)中數(shù)值計(jì)算的精度。時(shí)域模擬一般獲得的是近場(chǎng)電磁信息，為了得到諸如天線方向圖或散射體雷達(dá)散射截面之類的遠(yuǎn)場(chǎng)信息，必須獲得計(jì)算區(qū)域以外的頻域場(chǎng)或

10、瞬態(tài)場(chǎng)。多位學(xué)者在這方面做了許多工作，發(fā)展了一種高效的時(shí)域近遠(yuǎn)場(chǎng)變換方法23-26。借助這種方法，可以實(shí)現(xiàn)由計(jì)算區(qū)域內(nèi)近場(chǎng)數(shù)據(jù)到計(jì)算區(qū)域外遠(yuǎn)場(chǎng)數(shù)據(jù)的外推。目前，粗糙面散射的fdtd，傳遞函數(shù)在fdtd中的應(yīng)用，周期介質(zhì)、各向異性介質(zhì)、色散介質(zhì)和含有集中元件的fdtd，以及網(wǎng)絡(luò)并行fdtd技術(shù)等方面也取得了很大進(jìn)展。fdtd在迅速發(fā)展的同時(shí)，也獲得了非常廣泛的應(yīng)用。目前，它幾乎被應(yīng)用到了電磁場(chǎng)工程中的各個(gè)方面，例如：電磁散射、生物電磁計(jì)量學(xué)、輻射天線的分析、微波器件和導(dǎo)行波結(jié)構(gòu)的研究、散射和雷達(dá)截面的計(jì)算、周期結(jié)構(gòu)的分析、電子封裝和電磁兼容的分析、核電磁脈沖傳播和散射的分析、以及微光學(xué)元器件中光

11、的傳播和衍射特性的分析等。隨著新技術(shù)的不斷提出，其應(yīng)用范圍和成效正在迅速地?cái)U(kuò)大和提高。第二章 時(shí)域有限差分法的基本原理maxwell方程是描述宏觀電磁現(xiàn)象的一組基本方程。這組方程即可以寫成微分形式，又可以寫成積分形式。fdtd方法由maxwell旋度方程的微分形式出發(fā)，利用二階精度的中心差分近似，直接將微分運(yùn)算轉(zhuǎn)換為差分運(yùn)算，這樣達(dá)到了在一定體積內(nèi)和一段時(shí)間上對(duì)連續(xù)電磁場(chǎng)數(shù)據(jù)的抽樣壓縮。2.1 maxwell方程和yee氏算法 根據(jù)27中電磁場(chǎng)基本方程組的微分形式，若在無(wú)源空間，其空間中的媒質(zhì)是各向同性、線性和均勻的，即媒質(zhì)的參數(shù)不隨時(shí)間變化且各向同性，則maxwell旋度方程可寫成： （2-

12、1a） （2-1b）式中，是電場(chǎng)強(qiáng)度，單位為伏/米（v/m）；是磁場(chǎng)強(qiáng)度，單位為安/米（a/m）；表示介質(zhì)介電系數(shù)，單位為法拉/米（f/m）； 表示磁導(dǎo)系數(shù)，單位為亨利/米（h/m）；表示介質(zhì)電導(dǎo)率，單位為西門子/米（s/m）；表示導(dǎo)磁率，單位為歐姆/米（）。在直角坐標(biāo)系中，（2-1）式可化為如下六個(gè)標(biāo)量方程： （2-2） （2-3）這六個(gè)偏微分方程是fdtd算法的基礎(chǔ)。 k.s.yee3在1966年建立了如圖2-1所示的空間網(wǎng)格，這就是著名的yee氏元胞網(wǎng)格。圖2-1 yee氏網(wǎng)格及其電磁場(chǎng)分量分布并引入如下的差分近似方法對(duì)（2-2）、(2-3)式中的六個(gè)偏微分方程進(jìn)行了差分離散。令代表或在

13、直角坐標(biāo)系中某一分量，在時(shí)間和空間域中的離散可記為 (2-4)式中，、和分別是長(zhǎng)方體網(wǎng)格沿x、y、z方向的空間步長(zhǎng)，是時(shí)間步長(zhǎng)，i、j、k分別是沿x、y、z方向的網(wǎng)格編號(hào)，n是時(shí)間步數(shù)。對(duì)關(guān)于時(shí)間和空間的一階偏導(dǎo)數(shù)取中心差分近似，具有二階精度，即 (2-5a) (2-5b) (2-5c) (2-5d) 在fdtd中，空間上連續(xù)分布的電磁場(chǎng)物理量離散的空間排布如圖2-1所示。由圖可見，電場(chǎng)和磁場(chǎng)分量在空間交叉放置，使得在每個(gè)坐標(biāo)平面上每個(gè)電場(chǎng)分量被磁場(chǎng)環(huán)繞，每個(gè)磁場(chǎng)分量也被電場(chǎng)環(huán)繞。這種電磁場(chǎng)的空間結(jié)構(gòu)與電磁感應(yīng)和電磁波傳播的規(guī)律相符，在每一個(gè)網(wǎng)格單元都能滿足法拉第感應(yīng)定律和安培環(huán)流定律。各分量

14、的空間相對(duì)位置也適合于maxwell方程的差分計(jì)算，能夠恰當(dāng)?shù)孛枋鲭姶艌?chǎng)的傳播特性。同時(shí)，電場(chǎng)和磁場(chǎng)在時(shí)間上交替抽樣，抽樣時(shí)間間隔相差半個(gè)時(shí)間步，使maxwell旋度方程離散以后構(gòu)成顯式差分方程，從而可以在時(shí)間上迭代求解，而不需要進(jìn)行矩陣求逆運(yùn)算。因此，由給定相應(yīng)電磁問(wèn)題的初始條件，fdtd就可以逐步推進(jìn)地求得以后各個(gè)時(shí)刻空間電磁場(chǎng)的分布。2.2 fdtd的基本差分方程根據(jù)上述原則，可將（2-2）、（2-3）式離散為如下的差分方程形式： （2-6a） (2-6b) (2-6c) (2-6d) (2-6e) (2-6f) 式中， (2-7a)， (2-7b)(2-6)式就是fdtd的基本差分方程

15、組。從式中可以看出，方程組中含有半個(gè)空間步和半個(gè)時(shí)間步，為了便于編程，可將（2-6）式改寫成如下形式28：(2-8a)(2-8b)(2-8c)(2-8d)(2-8e)(2-8f)根據(jù)上述fdtd差分方程組可得出計(jì)算電磁場(chǎng)的時(shí)域推進(jìn)計(jì)算方法，如圖2-2所示。已知t1=t0= 0時(shí)刻空間各處的電磁場(chǎng)初始值計(jì)算t2=t1+ 時(shí)刻空間各處的磁場(chǎng)值計(jì)算t1=t2+ 時(shí)刻空間各處的電場(chǎng)值循環(huán)n次 圖2-2 fdtd在時(shí)域的交叉半步逐步推進(jìn)計(jì)算式(2-8a)(2-8c)的等號(hào)左邊的電場(chǎng)值是第n次循環(huán)的電場(chǎng)值，等號(hào)右邊的電場(chǎng)值是第n-1次循環(huán)存儲(chǔ)在內(nèi)存中的電場(chǎng)值，磁場(chǎng)值是本次循環(huán)計(jì)算得到的磁場(chǎng)值；式(2-8d

16、)(2-8f) 等號(hào)左邊的磁場(chǎng)值是第n次循環(huán)的磁場(chǎng)值，等號(hào)右邊的磁場(chǎng)值和電場(chǎng)值都是第n-1次循環(huán)存儲(chǔ)在內(nèi)存中的場(chǎng)值。這樣，就解決了半個(gè)時(shí)間步在程序中無(wú)法表示的問(wèn)題，而且也沒(méi)有破壞電磁場(chǎng)在時(shí)間上逐步推進(jìn)的邏輯關(guān)系。2.3 時(shí)域有限差分法相關(guān)技術(shù)2.3.1 數(shù)值穩(wěn)定性問(wèn)題上述fdtd方程是一種顯式差分方程，在執(zhí)行時(shí)，存在一個(gè)重要的問(wèn)題：即算法的穩(wěn)定性問(wèn)題。這種不穩(wěn)定性表現(xiàn)為在解顯式方程時(shí)，隨著時(shí)間步數(shù)的繼續(xù)增加，計(jì)算結(jié)果也將無(wú)限制地增加。taflove等4于1975年對(duì)yee氏差分格式的穩(wěn)定性進(jìn)行了討論，并導(dǎo)出了對(duì)時(shí)間步長(zhǎng)的限制條件。數(shù)值解是否穩(wěn)定主要取決于時(shí)間步長(zhǎng)與空間步長(zhǎng)、的關(guān)系。對(duì)于非均勻媒

17、質(zhì)構(gòu)成的計(jì)算空間選用如下的穩(wěn)定性條件： (2-9)(2-9)式是空間和時(shí)間離散之間應(yīng)當(dāng)滿足的關(guān)系，又稱為courant穩(wěn)定性條件。若采用均勻立方體網(wǎng)格： (2-10)其中，為計(jì)算空間中的電磁波的最大速度。2.3.2 數(shù)值色散 fdtd方程組是對(duì)maxwell旋度方程進(jìn)行差分近似，在進(jìn)行數(shù)值計(jì)算時(shí)，將會(huì)在計(jì)算網(wǎng)格中引起數(shù)字波模的色散，即在fdtd網(wǎng)格中，電磁波的相速與頻率有關(guān)，電磁波的相速度隨波長(zhǎng)、傳播方向及變量離散化的情況不同而改變。這種關(guān)系由非物理因素引起，且色散將導(dǎo)致非物理因素引起的脈沖波形畸變、人為的各向異性和虛假折射等現(xiàn)象。顯然，色散與空間、時(shí)間的離散間隔有關(guān)，如下式所示： (2-11

18、)式(2-11)是三維情況下在fdtd方法中的單色平面波數(shù)值色散關(guān)系的一般形式，它表明fdtd計(jì)算中波的傳播速度與傳播方向有關(guān)。式中、分別是波矢量沿、方向的分量，是角頻率，是被模擬的均勻介質(zhì)中的光速。與數(shù)值色散關(guān)系相對(duì)應(yīng)，在無(wú)耗介質(zhì)中的單色平面波，色散解析關(guān)系是： (2-12)由式(2-11)可知，當(dāng)式(2-11)中的、均趨于零時(shí)，它就趨于式(2-12)。也就是說(shuō)數(shù)值色散是由于用近似差分替代連續(xù)微分而引起的，而且在理論上可以減小到任意程度，只要此時(shí)時(shí)間步長(zhǎng)和空間步長(zhǎng)都足夠小，但這將大大增加所需的計(jì)算機(jī)存儲(chǔ)空間和計(jì)算時(shí)間，并使累積誤差增加。因此，在實(shí)際計(jì)算中要根據(jù)問(wèn)題的性質(zhì)和計(jì)算機(jī)的軟硬件條件來(lái)

19、選擇合適的時(shí)間步長(zhǎng)和空間步長(zhǎng)。為獲得理想的色散關(guān)系，問(wèn)題空間分割應(yīng)按照小于正常網(wǎng)格的原則進(jìn)行。一般選取的最大空間步長(zhǎng)為，為所研究范圍內(nèi)電磁波的最小波長(zhǎng)。由上分析說(shuō)明，數(shù)值色散在用fdtd法分析電磁場(chǎng)傳播中的影響是不可能避免的，但我們可以盡可能的減小數(shù)值色散的影響。2.3.3 離散網(wǎng)格的確定 無(wú)論是簡(jiǎn)單目標(biāo)還是復(fù)雜目標(biāo)，在進(jìn)行fdtd離散時(shí)網(wǎng)格尺寸的確定，除了受計(jì)算資源的限制不可能取得很小外，還需要考慮以下幾個(gè)因素：1目標(biāo)離散精確度的要求。網(wǎng)格應(yīng)當(dāng)足夠小以便能精確模擬目標(biāo)幾何形狀和電磁參數(shù)。2fdtd方法本身的要求。主要是考慮色散誤差的影響。設(shè)網(wǎng)格為立方體，所關(guān)心頻段的頻率上限為，對(duì)應(yīng)波長(zhǎng)為，則

20、考慮fdtd的數(shù)值色散要求 （2-13）通常。上式是根據(jù)已知所關(guān)心頻率上限情況下來(lái)確定fdtd網(wǎng)格尺寸的；反之，若給定，則fdtd計(jì)算結(jié)果可用的上限頻率也隨之確定。3入射波的要求。入射波的上限截止頻率應(yīng)包含所關(guān)心頻率范圍，即。2.4 吸收邊界條件 由時(shí)域有限差分法的基本原理可知，在利用時(shí)域有限差分法研究電磁場(chǎng)時(shí)，需在全部問(wèn)題空間建立yee氏網(wǎng)格空間，并存儲(chǔ)每個(gè)單元網(wǎng)格上任一時(shí)間步的六個(gè)場(chǎng)分量用于下一時(shí)間步的計(jì)算。而在對(duì)于輻射、散射這類開放系統(tǒng)的實(shí)際研究中，不可能有無(wú)限大的存儲(chǔ)空間。因此，必須在某處將網(wǎng)格空間截?cái)啵以诮財(cái)噙吔缇W(wǎng)格點(diǎn)處運(yùn)用特殊的場(chǎng)分量計(jì)算方法，使得向邊界面行進(jìn)的波在邊界處保持“外

21、向行進(jìn)”特性、無(wú)明顯的反射現(xiàn)象，并且不會(huì)使內(nèi)部空間的場(chǎng)產(chǎn)生畸變，從而用有限網(wǎng)格空間模擬電磁波在無(wú)界空間中傳播的情況。具有這種功能的邊界條件稱之為吸收邊界條件，或輻射邊界條件，或網(wǎng)格截?cái)鄺l件2931，如圖2-3所示。圖 2-3 附加截?cái)噙吔缡褂?jì)算區(qū)域變?yōu)橛邢抻驈膄dtd的基本差分方程組可以看出，在截?cái)噙吔缑嫔锨邢驁?chǎng)分量的計(jì)算需要利用計(jì)算空間以外的電磁場(chǎng)分量，因此fdtd基本差分方程對(duì)這些截?cái)噙吔缑嫔系膱?chǎng)分量失效。如何處理截?cái)噙吔缟系膱?chǎng)分量，使之與需要考慮的無(wú)限空間有盡量小的差異，是fdtd中必須很好解決的一個(gè)重要問(wèn)題。實(shí)際上，這是要求在誤差可容忍的范圍內(nèi)，計(jì)算空間中的外向波能夠順利通過(guò)截?cái)噙吔缑?/p>

22、而不引起波的明顯反射，使有限計(jì)算空間的數(shù)值模擬與實(shí)際情況趨于一致，對(duì)外向波而言，就像在無(wú)限大空間中傳播一樣。所以，需要一種截?cái)噙吔缇W(wǎng)格處的特殊計(jì)算方法，它不僅要保證邊界場(chǎng)計(jì)算的必要精度，而且還要大大消除非物理因素引起的波反射，使得用有限的網(wǎng)格空間就能模擬電磁波在無(wú)限空間中的傳播。但是如果處理不當(dāng)，截?cái)噙吔缑婵赡茉斐奢^大反射，構(gòu)成數(shù)值模擬誤差的一部分，甚至可能造成算法不穩(wěn)定。加于截?cái)噙吔鐖?chǎng)分量符合上述要求的算法就稱為吸收邊界條件（absorbing boundary conditions）。2.4.1 一階和二階近似吸收邊界條件 在截?cái)噙吔绺浇ǔ](méi)有激勵(lì)源?？紤]齊次波動(dòng)方程 (2-14)式中，

23、表示直角坐標(biāo)系下任意電磁場(chǎng)分量。 b.engquist和a.majda15利用偏微分算子對(duì)式(2-13)作因式分解，并分別取其taylor級(jí)數(shù)展開式中的第一項(xiàng)和前兩項(xiàng)近似，導(dǎo)出了適合直角坐標(biāo)系下fdtd吸收邊界條件的單向波動(dòng)方程，這就是engquist-majda吸收邊界條件。設(shè)三維長(zhǎng)方體fdtd區(qū)域0xa，0yb，0z=np+1&i=n1-np di=0; elseif i=n1-np+1 di=np+i-n1-0.5; end %di是采樣點(diǎn)橫向距pml內(nèi)邊界的距離 sigma_mx=sigma_max*(di/np)m; for j=1:n1 if j=np+1&j=n1-np dj=0

24、; elseif j=n1-np+1 dj=np+j-n1-0.5; end %dj是采樣點(diǎn)縱向距pml內(nèi)邊界的距離 sigma_my=sigma_max*(dj/np)m; if sigma_mx+sigma_my=0 %真空區(qū) if j=100&i=100 t=30; %可選擇高斯脈沖 term=(tt-t); pulse=exp(-4*pi*term2/202); pulse=sin(2*pi*tt/40); %可選正弦時(shí)諧源 c_miu=c*del_t/del_s; eterm1=c_miu*(ex(i,j+1)-ex(i,j); eterm2=c_miu*(ey(i+1,j)-ey(

25、i,j); bz(i,j)=bz(i,j)+eterm1-eterm2+pulse;%加入脈沖源 else c_miu=c*del_t/del_s; eterm1=c_miu*(ex(i,j+1)-ex(i,j); eterm2=c_miu*(ey(i+1,j)-ey(i,j); bz(i,j)=bz(i,j)+eterm1-eterm2; end else %pml區(qū) cpm=(1-2*c*sigma_mx*del_t)/(1+2*c*sigma_mx*del_t); cqm=c*del_t/(1+2*c*sigma_mx*del_t)/del_s; bzx(i,j)=cpm*bzx(i,j

26、)-cqm*(ey(i+1,j)-ey(i,j); cpm=(1-2*c*sigma_my*del_t)/(1+2*c*sigma_my*del_t); cqm=c*del_t/(1+2*c*sigma_my*del_t)/del_s; bzy(i,j)=cpm*bzy(i,j)+cqm*(ex(i,j+1)-ex(i,j); bz(i,j)=bzx(i,j)+bzy(i,j); end end end for i=2:n1 if i=np+1&i=n1-np di=0; elseif i=n1-np+1 di=np+i-n1-1; end %di是采樣點(diǎn)橫向距pml內(nèi)邊界的距離 sigma_

27、ex=sigma_max*(di/np)m; for j=1:n1 cam=(1-2*c*sigma_ex*del_t)/(1+2*c*sigma_ex*del_t); cbm=c*del_t/(1+2*c*sigma_ex*del_t)/del_s; ey(i,j)=cam*ey(i,j)-cbm*(bz(i,j)-bz(i-1,j); end end for i=1:n1 for j=2:n1 if j=np+1&j=n1-np dj=0; elseif j=n1-np+1 dj=np+j-n1-1; end %dj是采樣點(diǎn)縱向距pml內(nèi)邊界的距離 sigma_ey=sigma_max*(

28、dj/np)m; cam=(1-2*c*sigma_ey*del_t)/(1+2*c*sigma_ey*del_t); cbm=c*del_t/(1+2*c*sigma_ey*del_t)/del_s; ex(i,j)=cam*ex(i,j)+cbm*(bz(i,j)-bz(i,j-1); end end bzxx(tt,1)=bz(90,50); %對(duì)靠近邊界的中央磁場(chǎng)點(diǎn)采樣 bzxx(tt,2)=bz(50,90);3.2 出圖及結(jié)果3.2.1程序部分 figure(1); %可視化處理 clf; mesh(bz); %磁場(chǎng)的幅值 axis(0 n 0 n -0.5 0.5); xlabe

29、l(i) ylabel(j) drawnow;endfigure(2);plot(bzxx);figure(3);title(磁場(chǎng)幅值分布圖);surface(bz);shading interp;axis square;toc3.2.2 所出的效果圖正弦時(shí)諧場(chǎng)源輻射效果圖高斯脈沖輻射效果圖第四章 結(jié) 論通過(guò)對(duì)時(shí)域差分法的學(xué)習(xí)而進(jìn)行的關(guān)于二維fdtd的te波仿真，鞏固了迭代式的形式與相應(yīng)的循環(huán)編程技巧。但也存在一定問(wèn)題，由于在吸收上我的編程思路不是很規(guī)范，所以在吸收效果上與理論上的完全匹配層吸收效果還有差距，但是已經(jīng)基本達(dá)到了我的預(yù)期目的?？赡苁俏业姆謪^(qū)方法不對(duì)的緣故，原因就是出在sigma_

30、max的設(shè)置上面，在高斯制下長(zhǎng)度單位變?yōu)閏m，所以同樣的空間步長(zhǎng)在高斯制下數(shù)量級(jí)就擴(kuò)大了100倍，而sigma_max的計(jì)算式里面分母含有空間步長(zhǎng)，所以如果不加以修正，sigma_max在高斯制下就變成了國(guó)際制下的1%，太大的導(dǎo)電率和太小的導(dǎo)電率都不能實(shí)現(xiàn)較好的吸收效果，下階段應(yīng)將其做適當(dāng)修正并且寫粒子仿真時(shí)用高斯制maxwell方程可以參考角點(diǎn)反射大，多設(shè)置幾層pml試一試。參考文獻(xiàn)1j.m.jin, m.zunoubi, k.c.donepudi, w.c.chew, frequency-domain and time-domain finite-element solution of m

31、axwells equations using spectral lanczos decomposition method. comput. methods appl. mech. engrg. 1999,169:2792962k.dongsoo,h.b.lee. hybird full-wave analysis of via-hole grounds using finit-difference and finite-element time-domain methods. ieee transactions on microwave theory and techniques. 1997

32、,45(12):221722223 kane s. yee numerical solution of initial boundary value problems involving maxwells equations in isotropic mediaj.ieee trans antennas propagation, 1966,ap-14(3):3023074 a.taylor and m.e.brodwin. numerical solution of steady-state em scattering problems using the time-maxwells equa

33、tions. j ieee trans.microwave theroy tech.,aug. 1975，mtt-23:6236305 x.min,w.sun and k.m.chen. stability analysis of the finite difference time domain method applied to unbounded electromagnetic problem j. ieee antennas and propagation society international symposium,dallas:may 1990,(4):164016436a.ta

34、flove. review of the formulation and applications of the finite-difference time-domain method for numerical modeling of electromagnetic wave interactions with arbitrary structures j. wave motion,june 1988,10(6): 5475827 i.s.kim and w.j.r.hoefer. numerical dispersion characteristics and stability fac

35、tor for the td-fd method j.electronics letters, july 1990, 26(7): 4854878 k.l.shlager,j.g.msloney, s.l.ray,and etc. relative accuracy of several finite-difference time-domain methods in two and three dimensions j. ieee transactions on antennas and propagation,dec.1993,ap-41(12)：173217379 s.l.ray. nu

36、merical dispersion and stability characteristics of finite-difference time-domain methods on nonorthogonal meshes j. ieee transactions on antennas and propagation,feb.1993,ap-41(2)：23323510 h.shi and j.l.drewriak. dispersive comparison for dsi- and tenor-based nonorthogonal fdtd j. ieee microwave an

37、d guided wave letters, may 1996,6(5):19319511 f.xiao and h.yabe. numerical dispersion relation for fdtd method in general curvilinear coordinates j. ieee microwave and guided wave letters, feb.1997,7(2):485012 r.holland and j.williamstotal field versus scattered field finite difference codes:a compa

38、rative assessmentj.ieee transactions on nuclear science，dec.1983，ns-30(6)：4583458813g.murabsorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic field equationsjieee transactions on electromagnetic compatibility，nov.1981，emc-23(4)：37738214b.engquist a

39、nd a.majda. absorbing boundary conditions for the numerical simulation of waves j.mathematics of the computation，july 1977，31(139)：62965115j.p.berengera perfectly matched layer for the absorption of electromagnetic wavesjjournal of computational physics，oct.1994，114(2)：18520016j.p.berengerperfectly

40、matched layer for the fdtd solution of wave-structure interaction problemjieee transactions on antennas propagation，jan.1996，ap-44(1)：11011717j.p.berengerthree-dimensional perfectly matched layer for the absorption of electromagnetic wavesjjournal of computational physics，sept.1996，127(2)：36337918z.s.s