計(jì)算電磁學(xué)之FDTD算法的MATLAB語言實(shí)現(xiàn)

1、South China Normal University課程設(shè)計(jì)實(shí)驗(yàn)報(bào)告課程名稱： 計(jì)算電磁學(xué) 指導(dǎo)老師： 專業(yè)班級(jí)： 2014級(jí) 電路與系統(tǒng)姓 名： 學(xué) 號(hào)： FDTD算法的MATLAB語言實(shí)現(xiàn)摘要：時(shí)域有限差分(FDTD)算法是K.S.Yee于1966年提出的直接對(duì)麥克斯韋方程作差分處理,用來解決電磁脈沖在電磁介質(zhì)中傳播和反射問題的算法。其基本思想是:FDTD計(jì)算域空間節(jié)點(diǎn)采用Yee元胞的方法,同時(shí)電場(chǎng)和磁場(chǎng)節(jié)點(diǎn)空間與時(shí)間上都采用交錯(cuò)抽樣；把整個(gè)計(jì)算域劃分成包括散射體的總場(chǎng)區(qū)以及只有反射波的散射場(chǎng)區(qū),這兩個(gè)區(qū)域是以連接邊界相連接,最外邊是采用特殊的吸收邊界,同時(shí)在這兩個(gè)邊界之間有個(gè)輸出邊

2、界,用于近、遠(yuǎn)場(chǎng)轉(zhuǎn)換;在連接邊界上采用連接邊界條件加入入射波,從而使得入射波限制在總場(chǎng)區(qū)域；在吸收邊界上采用吸收邊界條件,盡量消除反射波在吸收邊界上的非物理性反射波。本文主要結(jié)合FDTD算法邊界條件特點(diǎn),在特定的參數(shù)設(shè)置下，用MATLAB語言進(jìn)行編程，在二維自由空間TEz網(wǎng)格中，實(shí)現(xiàn)脈沖平面波。關(guān)鍵詞：FDTD；MATLAB；算法1 緒論1.1 課程設(shè)計(jì)背景與意義 20世紀(jì)60年代以來，隨著計(jì)算機(jī)技術(shù)的發(fā)展，一些電磁場(chǎng)的數(shù)值計(jì)算方法逐步發(fā)展起來，并得到廣泛應(yīng)用，其中主要有：屬于頻域技術(shù)的有限元法（FEM）、矩量法(MM)和單矩法等；屬于時(shí)域技術(shù)方面的時(shí)域有限差分法(FDTD)、傳輸線矩陣法(T

3、LM)和時(shí)域積分方程法等。其中FDTD是一種已經(jīng)獲得廣泛應(yīng)用并且有很大發(fā)展前景的時(shí)域數(shù)值計(jì)算方法。時(shí)域有限差分（FDTD）方法于1966年由K.S.Yee提出并迅速發(fā)展，且獲得廣泛應(yīng)用。K.S.Yee用后來被稱作Yee氏網(wǎng)格的空間離散方式，把含時(shí)間變量的Maxwell旋度方程轉(zhuǎn)化為差分方程，并成功地模擬了電磁脈沖與理想導(dǎo)體作用的時(shí)域響應(yīng)。但是由于當(dāng)時(shí)理論的不成熟和計(jì)算機(jī)軟硬件條件的限制，該方法并未得到相應(yīng)的發(fā)展。20世紀(jì)80年代中期以后，隨著上述兩個(gè)條件限制的逐步解除，F(xiàn)DTD便憑借其特有的優(yōu)勢(shì)得以迅速發(fā)展。它能方便、精確地預(yù)測(cè)實(shí)際工程中的大量復(fù)雜電磁問題，應(yīng)用范圍幾乎涉及所有電磁領(lǐng)域，成為電

4、磁工程界和理論界研究的一個(gè)熱點(diǎn)。目前，F(xiàn)DTD日趨成熟，并成為分析大部分實(shí)際電磁問題的首選方法。1.2 時(shí)域有限差分法的發(fā)展與應(yīng)用經(jīng)過四十多年的發(fā)展，F(xiàn)DTD已發(fā)展成為一種成熟的數(shù)值計(jì)算方法。在發(fā)展過程中，幾乎都是圍繞幾個(gè)重要問題展開的，即數(shù)值穩(wěn)定性、計(jì)算精度、數(shù)值色散、激勵(lì)源技術(shù)以及開域電磁問題的吸收邊界條件等。數(shù)值穩(wěn)定和計(jì)算精度對(duì)任何一種數(shù)值計(jì)算方法都是至關(guān)重要的。其中激勵(lì)源的設(shè)計(jì)和引入也是FDTD的一個(gè)重要任務(wù)。目前，應(yīng)用最廣泛的激勵(lì)源引入技術(shù)是總場(chǎng)/散射場(chǎng)體系。對(duì)于散射問題，通常在FDTD計(jì)算空間中引入連接邊界，它將整個(gè)計(jì)算空間劃分為內(nèi)部的總場(chǎng)區(qū)和外部的散射場(chǎng)區(qū)，如圖1.1。利用Huy

5、gens原理，可以在連接邊界處引入入射場(chǎng)，使入射場(chǎng)的加入變得簡(jiǎn)單易行。圖1.1總場(chǎng)/散射場(chǎng)區(qū)2 FDTD法基本原理 2.1 Maxwell方程和Yee氏算法根據(jù)電磁場(chǎng)基本方程組的微分形式，若在無源空間，其空間中的煤質(zhì)是各向同性、線性和均勻的，即煤質(zhì)的參數(shù)不隨時(shí)間變化且各向同性，則Maxwell旋度方程可寫成： 式中，E是電場(chǎng)強(qiáng)度，單位為伏/米（V/m）；H是磁場(chǎng)強(qiáng)度，單位為安/米（A/m）；表示介質(zhì)介電系數(shù)，單位為法拉/米（F/m）；表示磁導(dǎo)系數(shù)，單位為亨利/米（H/m）；表示介質(zhì)電導(dǎo)率，單位為西門子/米（S/m）；m表示導(dǎo)磁率，單位為歐姆/米（/m）。 K.S.Yee在1966年建立了如圖2

6、.1所示的空間網(wǎng)格，這就是著名的Yee氏元胞網(wǎng)格。圖2.1 Yee氏網(wǎng)格及其電磁場(chǎng)分量分布 電場(chǎng)和磁場(chǎng)被交叉放置,電場(chǎng)分量位于網(wǎng)格單元每條棱的中心,磁場(chǎng)分量位于網(wǎng)格單元每個(gè)面的中心,每個(gè)磁場(chǎng)(電場(chǎng))分量都有4個(gè)電場(chǎng)(磁場(chǎng))分量環(huán)繞。這樣不僅保證了介質(zhì)分界面上切向場(chǎng)分量的連續(xù)性條件得到自然滿足,而且還允許旋度方程在空間上進(jìn)行中心差分運(yùn)算,同時(shí)也滿足了法拉第電磁感應(yīng)定律和安培環(huán)路積分定律,也可以很恰當(dāng)?shù)啬M電磁波的實(shí)際傳播過程。2.2 時(shí)域有限差分法相關(guān)技術(shù)2.2.1數(shù)值穩(wěn)定性問題FDTD方程是一種顯式差分方程，在執(zhí)行時(shí)，存在一個(gè)重要的問題：即算法的穩(wěn)定性問題。這種不穩(wěn)定性表現(xiàn)為在解顯式方程時(shí)，隨

7、著時(shí)間步數(shù)的繼續(xù)增加，計(jì)算結(jié)果也將無限制地增加。Taflove于1975年對(duì)Yee氏差分格式的穩(wěn)定性進(jìn)行了討論，并導(dǎo)出了對(duì)時(shí)間步長(zhǎng)的限制條件。數(shù)值解是否穩(wěn)定主要取決于時(shí)間步長(zhǎng)t與空間步長(zhǎng)x、y、z的關(guān)系。對(duì)于非均勻媒質(zhì)構(gòu)成的計(jì)算空間選用如下的穩(wěn)定性條件：上式是空間和時(shí)間離散之間應(yīng)當(dāng)滿足的關(guān)系，又稱為Courant穩(wěn)定性條件。若采用均勻立方體網(wǎng)格：x=y=z=sD其中，v為計(jì)算空間中的電磁波的最大速度。2.2.2數(shù)值色散 FDTD方程組是對(duì)Maxwell旋度方程進(jìn)行差分近似，在進(jìn)行數(shù)值計(jì)算時(shí)，將會(huì)在計(jì)算網(wǎng)格中引起數(shù)字波模的色散，即在FDTD網(wǎng)格中，電磁波的相速與頻率有關(guān)，電磁波的相速度隨波長(zhǎng)、傳

8、播方向及變量離散化的情況不同而改變。這種關(guān)系由非物理因素引起，且色散將導(dǎo)致非物理因素引起的脈沖波形畸變、人為的各向異性和虛假折射等現(xiàn)象。為獲得理想的色散關(guān)系，問題空間分割應(yīng)按照小于正常網(wǎng)格的原則進(jìn)行。一般選取的最大空間步長(zhǎng)為max=min/20，min為所研究范圍內(nèi)電磁波的最小波長(zhǎng)。由上分析說明，數(shù)值色散在用FDTD法分析電磁場(chǎng)傳播中的影響是不可能避免的，但我們可以盡可能的減小數(shù)值色散的影響。2.2.3離散網(wǎng)格的確定 無論是簡(jiǎn)單目標(biāo)還是復(fù)雜目標(biāo)，在進(jìn)行FDTD離散時(shí)網(wǎng)格尺寸的確定，除了受計(jì)算資源的限制不可能取得很小外，還需要考慮以下幾個(gè)因素：1目標(biāo)離散精確度的要求。網(wǎng)格應(yīng)當(dāng)足夠小以便能精確模擬

9、目標(biāo)幾何形狀和電磁參數(shù)。2FDTD方法本身的要求。主要是考慮色散誤差的影響。設(shè)網(wǎng)格為立方體x=y=z=，所關(guān)心頻段的頻率上限為f max，對(duì)應(yīng)波長(zhǎng)為f minl，則考慮FDTD的數(shù)值色散要求min /N通常N10。上式是根據(jù)已知所關(guān)心頻率上限情況下來確定FDTD網(wǎng)格尺寸d的；反之，若給定，則FDTD計(jì)算結(jié)果可用的上限頻率也隨之確定。2.3 吸收邊界條件由時(shí)域有限差分法的基本原理可知，在利用時(shí)域有限差分法研究電磁場(chǎng)時(shí)，需在全部問題空間建立Yee氏網(wǎng)格空間，并存儲(chǔ)每個(gè)單元網(wǎng)格上任一時(shí)間步的六個(gè)場(chǎng)分量用于下一時(shí)間步的計(jì)算。而在對(duì)于輻射、散射這類開放系統(tǒng)的實(shí)際研究中，不可能有無限大的存儲(chǔ)空間。因此，必

10、須在某處將網(wǎng)格空間截?cái)啵以诮財(cái)噙吔缇W(wǎng)格點(diǎn)處運(yùn)用特殊的場(chǎng)分量計(jì)算方法，使得向邊界面行進(jìn)的波在邊界處保持“外向行進(jìn)”特性、無明顯的反射現(xiàn)象，并且不會(huì)使內(nèi)部空間的場(chǎng)產(chǎn)生畸變，從而用有限網(wǎng)格空間模擬電磁波在無界空間中傳播的情況。具有這種功能的邊界條件稱之為吸收邊界條件，或輻射邊界條件，或網(wǎng)格截?cái)鄺l件，如圖2.2所示圖2.2附加截?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)分量，使之與需要考慮的無限空間有盡量小的差異，是FDTD中必須很好解決的

11、一個(gè)重要問題。實(shí)際上，這是要求在誤差可容忍的范圍內(nèi)，計(jì)算空間中的外向波能夠順利通過截?cái)噙吔缑娑灰鸩ǖ拿黠@反射，使有限計(jì)算空間的數(shù)值模擬與實(shí)際情況趨于一致，對(duì)外向波而言，就像在無限大空間中傳播一樣。所以，需要一種截?cái)噙吔缇W(wǎng)格處的特殊計(jì)算方法，它不僅要保證邊界場(chǎng)計(jì)算的必要精度，而且還要大大消除非物理因素引起的波反射，使得用有限的網(wǎng)格空間就能模擬電磁波在無限空間中的傳播。但是如果處理不當(dāng)，截?cái)噙吔缑婵赡茉斐奢^大反射，構(gòu)成數(shù)值模擬誤差的一部分，甚至可能造成算法不穩(wěn)定。加于截?cái)噙吔鐖?chǎng)分量符合上述要求的算法就稱為吸收邊界條件（AbsorbingBoundaryConditions）。2.4 FDTD計(jì)

12、算所需時(shí)間步的估計(jì) 為了使計(jì)算達(dá)到穩(wěn)定,通常計(jì)算所需要時(shí)間步按照電磁波往返穿越FDTD計(jì)算區(qū)對(duì)角線35次來估計(jì)。若FDTD計(jì)算區(qū)總元胞數(shù)為N，則對(duì)角線上元胞約為（三維)。按照Courant穩(wěn)定條件，設(shè)計(jì)算中心區(qū)t=/（2c），即穿越對(duì)角線一次需要時(shí)間步為?？傆?jì)算時(shí)間步約為 需() 。對(duì)于二維情況則約為 ?；蛘哒f，計(jì)算時(shí)間步大約等于FDTD計(jì)算區(qū)對(duì)角線上元胞數(shù)目的1220倍。實(shí)際上，計(jì)算所需時(shí)間步還與散射體具體形狀、結(jié)構(gòu)有關(guān)。圖2.3給出了應(yīng)用FDTD分析電磁場(chǎng)問題時(shí)的程序流程圖圖2.3FDTD程序流程圖3 課程設(shè)計(jì)要求及MATLAB語言實(shí)現(xiàn)3.1課程設(shè)計(jì)要求 采用總場(chǎng)/散射場(chǎng)技術(shù)，在二維自由空

13、間TEz網(wǎng)格中，實(shí)現(xiàn)脈沖平面波。 總場(chǎng)區(qū)域100*100個(gè)網(wǎng)格。 周邊散射場(chǎng)區(qū)域20個(gè)網(wǎng)格。 散射場(chǎng)區(qū)域外采用PML吸收邊界條件。 平面波0度角入射。使用高斯脈沖，中心頻率為1GHz,帶寬為1/e。 采用正方形Yee網(wǎng)格，離散步長(zhǎng)=1.5cm，時(shí)間離散步長(zhǎng)為=/2c.3.2 MATLAB程序% 此程序?qū)崿F(xiàn)自由空間的平面波 % 此程序采用總場(chǎng)/散射場(chǎng)技術(shù) % 此程序使用PML設(shè)置吸收邊界條件 %KE=128;%真空區(qū)域網(wǎng)格數(shù)IE=KE;JE=KE;T=0; NSTEP=300;%迭代次數(shù)e=2.71828;% 初始化 %dz=zeros(KE,KE);ez=zeros(KE,KE);hx=zer

14、os(KE,KE);hy=zeros(KE,KE);ihx=zeros(KE,KE);ihy=zeros(KE,KE);ga=ones(KE,KE);ez_inc=zeros(1,JE);hx_inc=zeros(1,JE);% PML 參數(shù) %gi2=ones(1,IE);gi3=gi2;fi1=zeros(1,IE);fi2=ones(1,IE);fi3=fi2;gj2=ones(1,JE);gj3=gj2;fj1=zeros(1,JE);fj2=ones(1,JE);fj3=fj2;npml=8;for i=1:8 %設(shè)置PML層中的參數(shù) xnum=npml-i; xd=npml; xx

15、n=xnum/xd; xn=0.33*xxn3; gi2(i)=1.0/(1.0+xn); gi2(IE-1-i)=1.0/(1.0+xn); gi3(i)=(1.0-xn)/(1.0+xn); gi3(IE-1-i)=(1.0-xn)/(1.0+xn); xxn=(xnum-0.5)/xd; xn=0.25*xxn3; fi1(i)=xn; fi1(IE-2-i)=xn; fi2(i)=1.0/(1.0+xn); fi2(IE-2-i)=1.0/(1.0+xn); fi3(i)=(1.0-xn)/(1.0+xn); fi3(IE-2-i)=(1.0-xn)/(1.0+xn); endfor

16、j=1:8 xnum=npml-j; xd=npml; xxn=xnum/xd; xn=0.33*xxn3; gj2(j)=1.0/(1.0+xn); gj2(JE-1-j)=1.0/(1.0+xn); gj3(j)=(1.0-xn)/(1.0+xn); gj3(JE-1-j)=(1.0-xn)/(1.0+xn); xxn=(xnum-0.5)/xd; xn=0.25*xxn3; fj1(j)=xn; fj1(JE-2-j)=xn; fj2(j)=1.0/(1.0+xn); fj2(JE-2-j)=1.0/(1.0+xn); fj3(j)=(1-xn)/(1.0+xn); fj3(JE-2-j

17、)=(1.0-xn)/(1.0+xn); end% 總場(chǎng)/散射場(chǎng)邊界 %ia=28;ib=IE-ia-1;ja=28;jb=JE-ja-1;ez_inc_low_m1=0;ez_inc_low_m2=0;ez_inc_high_m1=0;ez_inc_high_m2=0;% 信號(hào)源 %t0=80;%脈沖高度spread=12;%脈沖寬度% 網(wǎng)格 %ddx=0.015; % 離散步長(zhǎng) %dt=ddx/(2*3e8); % 計(jì)算時(shí)間離散步長(zhǎng) %epsz=8.851*1e-12; % 自由空間的介電常數(shù) % 迭代求解電場(chǎng)與磁場(chǎng) % for n=1:NSTEP a=1; T=T+1; for j=2:

18、JE ez_inc(j)= ez_inc(j)+0.5*( hx_inc(j-1)- hx_inc(j); end % 入射波緩沖區(qū) % ez_inc(1) =ez_inc_low_m2; ez_inc_low_m2=ez_inc_low_m1; ez_inc_low_m1=ez_inc(2); ez_inc(JE-1) =ez_inc_high_m2; ez_inc_high_m2=ez_inc_high_m1; ez_inc_high_m1=ez_inc(JE-2); % 計(jì)算 dx 區(qū)域 % for i=2:KE for j=2:KE dz(i,j)=gi3(i)*gj3(j)*dz(i

19、,j)+gi2(i)*gj2(j)*0.5*(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1); end end % 輸入高斯脈沖信號(hào)源 % pulse=exp(-0.5*(t0-T)2/spread2); % ez_inc(3)=pulse; % 入射波 DZ 參數(shù) % for i=ia:ib dz(i,ja)= dz(i,ja)+0.5*hx_inc(ja-1); dz(i,jb)= dz(i,jb)-0.5*hx_inc(jb); end % 計(jì)算電場(chǎng)dx區(qū)域 % for i=1:KE for j=1:KE ez(i,j)=ga(i,j)*dz(i,j); end

20、end % 設(shè)置電場(chǎng) EZ 邊緣為0, 作為PML的參數(shù)% for j=1:JE-1 ez(1,j)=0; ez(IE,j)=0; end for i=1:IE-1 ez(i,1)=0; ez(i,JE)=0; end for j=1:JE-1 hx_inc(j)=hx_inc(j)+0.5*( ez_inc(j)-ez_inc(j+1); end % 計(jì)算磁場(chǎng)hx區(qū)域 % for i=1:KE for j=1:KE-1 curl_e=ez(i,j)-ez(i,j+1); ihx(i,j)=ihx(i,j)+fi1(i)*curl_e; hx(i,j)=fj3(j)*hx(i,j)+fj2(j

21、)*0.5* (curl_e+ihx(i,j); end end % 入射磁場(chǎng)Hx區(qū)域參數(shù) % for i=ia:ib hx(i,ja-1)= hx(i,ja-1)+0.5*ez_inc(ja); hx(i,jb)= hx(i,jb)-0.5*ez_inc(jb); end % 計(jì)算磁場(chǎng)hy區(qū)域% for i=1:KE-1 for j=1:KE curl_e=ez(i+1,j)-ez(i,j); ihy(i,j)=ihy(i,j)+fj1(j)*curl_e; hy(i,j)=fi3(i)*hy(i,j)+fi2(i)*0.5* (curl_e+ihy(i,j); end end % 入射磁場(chǎng)hy區(qū)域參數(shù) % for j=ja:jb hy(ia-1,j)= hy(ia-1,j)-0.5*ez_inc(j); hy(ib,j)= hy(ib,j)+0.5*ez_inc(j); end % 圖形展示 % m=moviein(500); x=1:KE; y=1:KE; X,Y=meshgrid(x,y); surf(X,Y,ez); axis(0 KE 0 KE -1 1 ); set(gcf,Color, white, Number, off, Name, sprintf(Simulation FDTD 2D, Iteration = %i, n); titl