MATLAB編程在FDTD算法中的應(yīng)用

摘要：文章介紹了時(shí)域有限差分（FDTD）法的基本原理，推導(dǎo)了二維TM波的FDTD表達(dá)式。給出了MATLAB語(yǔ)言編程的步驟和應(yīng)注意的問題，并結(jié)合算例闡述了基于MATLAB仿真的基本方法，最后得出用MATLAB語(yǔ)言對(duì)FDTD算法仿真的幾點(diǎn)結(jié)論。

關(guān)鍵詞：MATLAB；FDTD；編程

時(shí)域有限差分（FDTD）法是在21世紀(jì)60年代由K.S.Yee首先提出并用于求解電磁場(chǎng)散射問題，其主要思路是在空間軸和時(shí)間軸上對(duì)場(chǎng)量進(jìn)行離散，并用中心差分代替偏微分，這就將麥克斯韋方程組轉(zhuǎn)化為了差分方程，通過在時(shí)間軸和空間軸上采取蛙跳法（leapfrog）逐步推進(jìn)地求解，最終求得在一定邊值與初值條件下的空間場(chǎng)解。隨著計(jì)算機(jī)技術(shù)的發(fā)展，F(xiàn)DTD的應(yīng)用越來(lái)越多，對(duì)于FDTD算法的編程求解，最常用的程序語(yǔ)言有VC和FORTRUN，而MATLAB作為一種可視化效果好的軟件，在FDTD計(jì)算中可視化程度較高，并具有能顯示動(dòng)態(tài)場(chǎng)效果的特點(diǎn)。文章采用FDTD法對(duì)高壓開關(guān)柜內(nèi)超高頻電磁波的輻射和傳播特性進(jìn)行仿真，仿真中將對(duì)二維TM電磁波進(jìn)行FDTD表達(dá)式的推導(dǎo)，并結(jié)合FDTD算法邊界條件的特點(diǎn)，用MATLAB語(yǔ)言進(jìn)行編程。

1二維TM電磁波FDTD算法

1.1算式推導(dǎo)

在自由空間中，對(duì)于二維TM電磁波，，MAXWELL的兩個(gè)旋度方程可以分解為：

＝-0 （1）

＝0（2）

０=-（3）

構(gòu)造二維TM波FDTD cell如圖1所示：

按照FDTD元胞對(duì)以上三式中的偏導(dǎo)用中心差商代替，分別可得：

=-0 （4）

=0 （5）

０

=-（6）

由于方程中出現(xiàn)了半網(wǎng)格和半時(shí)間步，為了便于編程，可以將上面差分式改為如下FDTD算式：

Hx（i，j，k+1）=Hx（i，j，k）-[EZ（i，j+1，n）-EZ（i，j，n）]

（7）

Hy（i，j，n+1）=Hy（i，j，n）+[EZ（i+1，j，n）-EZ（i，j，n）]

（8）

Ez（i，j，n+1）

=Ez（i，j，n）+[-

] （9）

1.2數(shù)值穩(wěn)定性的條件

FDTD方法是以一組有限差分方程來(lái)替代MAXWELL旋度方程來(lái)進(jìn)行數(shù)值計(jì)算的方法，在執(zhí)行形如FDTD算法時(shí)，隨著時(shí)間步長(zhǎng)的增長(zhǎng)，如何保證該算法的穩(wěn)定性是一個(gè)重要問題。數(shù)值解是否穩(wěn)定主要取決于時(shí)間步長(zhǎng)?駐t與空間步長(zhǎng)?駐x、?駐y、?駐z之間的關(guān)系。按照Courant穩(wěn)定性條件，F(xiàn)DTD算法中空間和時(shí)間間隔應(yīng)滿足的關(guān)系為：

C?駐t≤[++]-1/2（10）

式（10）中C為真空中的光速。對(duì)于TM波的二維情形，令式中?駐z=∞，網(wǎng)格單元尺寸通常可取為?駐x＝?駐y，此時(shí)，C?駐t≤?駐x，一般選C?駐t=?駐x/2。當(dāng)?駐x、?駐y，不相等時(shí)，應(yīng)取二者中的較小值。當(dāng)沿兩個(gè)軸向的網(wǎng)格元可變時(shí)，即?駐x、?駐y分別是i，j的函數(shù)，則應(yīng)取每個(gè)軸向上的最小值，然后再先二者中的較小值。時(shí)間步長(zhǎng)應(yīng)滿足下式：

?駐t=（11）

才能保證該算法在較長(zhǎng)時(shí)間步長(zhǎng)上運(yùn)行的穩(wěn)定性。

2基于MATLAB的編程方法

2.1MATLAB編程方法

用MATLAB語(yǔ)言對(duì)本算例FDTD算法進(jìn)行求解的步驟可歸納為：

①正確建立二維TM波FDTD算法的數(shù)學(xué)模型。其中包括正確地劃分網(wǎng)格、推導(dǎo)出正確的FDTD算式等。本例中的算式已經(jīng)給出，下面介紹網(wǎng)格的劃分：根據(jù)算例的尺寸，可以將區(qū)域劃分為100×100個(gè)網(wǎng)格。其中空間步長(zhǎng)取為?駐x=?駐y=0.005，時(shí)間步長(zhǎng)取為?駐t=?駐x=2C，其中，C為自由空間的光速。?駐x、?駐y、?駐t在選取過程中必須滿足Courant穩(wěn)定性條件。雖然時(shí)間和空間步長(zhǎng)越小，F(xiàn)DTD運(yùn)算精度越高，但如果網(wǎng)格太細(xì)會(huì)使計(jì)算機(jī)運(yùn)算時(shí)間過長(zhǎng)，同時(shí)，在網(wǎng)格分辨率（即波長(zhǎng)與空間步長(zhǎng)之比）小到一定程度時(shí)，數(shù)值誤差就不再變化。

②確定邊界條件和設(shè)置激勵(lì)源。FDTD方法可用于求解電磁學(xué)中的散射和輻射等問題，為了在有限空間內(nèi)模擬無(wú)限大空間的問題，必須考慮吸收邊界條件。學(xué)用的吸收邊界有Mur 和PML吸收邊界。在一般情況下，對(duì)角區(qū)要求不高時(shí)，可用Mur作為邊界條件，而對(duì)計(jì)算精度要求較高時(shí)，可以采用PML吸收條件。對(duì)于本算例中，求解的是二維邊界條件，可以在程序中直接在四個(gè)邊界上將電場(chǎng)值賦0，同時(shí)將Ex、Ev、Ez賦0（TM波）即可。在FDTD算法中激勵(lì)源的設(shè)置非常重要，常用的激勵(lì)源有正弦激勵(lì)和高斯激勵(lì)，本算例中使用高斯脈沖，其表達(dá)式為：

Ei（t）=exp（-0.5×） （12）

在高斯脈沖中，根據(jù)經(jīng)驗(yàn)選擇合適的t0和t值。取t0=20 nS，t=5 s。

③繪制編寫MATLAB程序的流程圖，并編寫程序。MATLAB具有較強(qiáng)的可視化效果，本算例中，通過所編程序?qū)崿F(xiàn)三維彩圖的輸出，程序流程圖如圖2所示。

2.2運(yùn)行結(jié)果

圖3、4、5、6分別是運(yùn)算時(shí)間步為N=180、N=210、N=250、N=500時(shí)場(chǎng)的顯示圖。其中，N=180步時(shí)，波未到達(dá)邊界；N=210步時(shí)，電磁波剛好到達(dá)邊界；N=250步時(shí)，電磁波在邊界上反射較短時(shí)間后的情況；N=500步時(shí)，電磁波在邊界上反射較長(zhǎng)時(shí)間后的情況。

3結(jié)語(yǔ)

文章結(jié)合算例分析了使用MATLAB對(duì)FDTD算法進(jìn)行編程來(lái)求解電磁場(chǎng)問題的方法。FDTD在電磁場(chǎng)數(shù)值分析過程中，具有概念清晰、不用對(duì)大型矩陣求逆。因此編程比較簡(jiǎn)單，結(jié)合MATLAB強(qiáng)大的數(shù)據(jù)處理功能和圖形處理功能，能夠使FDTD算法的編程更加簡(jiǎn)明，編程效率更高，結(jié)果顯示更直觀。

參考文獻(xiàn)：

[1] 葛德彪，閆玉波.電磁波時(shí)域有限差分法[M].西安:西安電子 科技大學(xué)出版社，2002.

[2] 郭春波.時(shí)域有限差分法的MATLAB仿真[J].現(xiàn)代電子技 術(shù)，2005，8(199):45-46.

[3] A#8226;Taflove.Computational Electrodynamics-The finite-Difference Time-Domain Method[M].Artech Hourse，1995.

[4] [美]DeitelHM著.薛萬(wàn)明譯.C程序設(shè)計(jì)教程[M].北京:機(jī)械 工業(yè)出版社，2000.

[5] 石澄賢，鄭明芳.MATLAB語(yǔ)言及其在電子信息工程中的應(yīng)用[M].北京:清華大學(xué)出版社，2004.