1、時域有限差分法對平面te波的matlab仿真摘 要 時域有限差分法是由有限差分法發展出來的數值計算方法。自1966年yee在其論文中首次提出時域有限差分以來，時域有限差分法在電磁研究領域得到了廣泛的應用。主要有分析輻射條線、微波器件和導行波結構的研究、散射和雷達截面計算、分析周期結構、電子封裝和電磁兼容的分析、核電磁脈沖的傳播和散射以及在地面的反射及對電纜傳輸線的干擾、微光學元器件中光的傳播和衍射特性等等。由于電磁場是以場的形態存在的物質，具有獨特的研究方法，采取重疊的研究方法是其重要的特點，即只有理論分析、測量、計算機模擬的結果相互佐證，才可以認為是獲得了正確可信的結論。時域有限差分法就是實

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

20、考慮fdtd的數值色散要求 （2-13）通常。上式是根據已知所關心頻率上限情況下來確定fdtd網格尺寸的；反之，若給定，則fdtd計算結果可用的上限頻率也隨之確定。3入射波的要求。入射波的上限截止頻率應包含所關心頻率范圍，即。2.4 吸收邊界條件 由時域有限差分法的基本原理可知，在利用時域有限差分法研究電磁場時，需在全部問題空間建立yee氏網格空間，并存儲每個單元網格上任一時間步的六個場分量用于下一時間步的計算。而在對于輻射、散射這類開放系統的實際研究中，不可能有無限大的存儲空間。因此，必須在某處將網格空間截斷，且在截斷邊界網格點處運用特殊的場分量計算方法，使得向邊界面行進的波在邊界處保持“外

21、向行進”特性、無明顯的反射現象，并且不會使內部空間的場產生畸變，從而用有限網格空間模擬電磁波在無界空間中傳播的情況。具有這種功能的邊界條件稱之為吸收邊界條件，或輻射邊界條件，或網格截斷條件2931，如圖2-3所示。圖 2-3 附加截斷邊界使計算區域變為有限域從fdtd的基本差分方程組可以看出，在截斷邊界面上切向場分量的計算需要利用計算空間以外的電磁場分量，因此fdtd基本差分方程對這些截斷邊界面上的場分量失效。如何處理截斷邊界上的場分量，使之與需要考慮的無限空間有盡量小的差異，是fdtd中必須很好解決的一個重要問題。實際上，這是要求在誤差可容忍的范圍內，計算空間中的外向波能夠順利通過截斷邊界面

22、而不引起波的明顯反射，使有限計算空間的數值模擬與實際情況趨于一致，對外向波而言，就像在無限大空間中傳播一樣。所以，需要一種截斷邊界網格處的特殊計算方法，它不僅要保證邊界場計算的必要精度，而且還要大大消除非物理因素引起的波反射，使得用有限的網格空間就能模擬電磁波在無限空間中的傳播。但是如果處理不當，截斷邊界面可能造成較大反射，構成數值模擬誤差的一部分，甚至可能造成算法不穩定。加于截斷邊界場分量符合上述要求的算法就稱為吸收邊界條件（absorbing boundary conditions）。2.4.1 一階和二階近似吸收邊界條件 在截斷邊界附近通常沒有激勵源。考慮齊次波動方程 (2-14)式中，

23、表示直角坐標系下任意電磁場分量。 b.engquist和a.majda15利用偏微分算子對式(2-13)作因式分解，并分別取其taylor級數展開式中的第一項和前兩項近似，導出了適合直角坐標系下fdtd吸收邊界條件的單向波動方程，這就是engquist-majda吸收邊界條件。設三維長方體fdtd區域0xa，0yb，0z=np+1&i=n1-np di=0; elseif i=n1-np+1 di=np+i-n1-0.5; end %di是采樣點橫向距pml內邊界的距離 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是采樣點縱向距pml內邊界的距離 sigma_my=sigma_max*(dj/np)m; if sigma_mx+sigma_my=0 %真空區 if j=100&i=100 t=30; %可選擇高斯脈沖 term=(tt-t); pulse=exp(-4*pi*term2/202); pulse=sin(2*pi*tt/40); %可選正弦時諧源 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區 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是采樣點橫向距pml內邊界的距離 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是采樣點縱向距pml內邊界的距離 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); %對靠近邊界的中央磁場點采樣 bzxx(tt,2)=bz(50,90);3.2 出圖及結果3.2.1程序部分 figure(1); %可視化處理 clf; mesh(bz); %磁場的幅值 axis(0 n 0 n -0.5 0.5); xlabe

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

30、max的設置上面，在高斯制下長度單位變為cm，所以同樣的空間步長在高斯制下數量級就擴大了100倍，而sigma_max的計算式里面分母含有空間步長，所以如果不加以修正，sigma_max在高斯制下就變成了國際制下的1%，太大的導電率和太小的導電率都不能實現較好的吸收效果，下階段應將其做適當修正并且寫粒子仿真時用高斯制maxwell方程可以參考角點反射大，多設置幾層pml試一試。參考文獻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