時(shí)域有限差分法的Matlab仿真

張通　孫晶

【摘 要】文章介紹了時(shí)域有限差分法的基本原理，利用matlab仿真，實(shí)現(xiàn)了用時(shí)域有限差分程序來計(jì)算二維問題空間中的電場分布。

【關(guān)鍵詞】時(shí)域有限差分法；Matlab；電場分布

Simulation of Finite Difference Time Domain Method Using Matlab

ZHANG Tong SUN Jing

（College of Physics，Mechanical and Electrical Engineering，Jishou University， JiShou Hunan 416000，China）

【Abstract】The basic principle of finite difference time domain is introduced in this paper.With two-dimensional finite difference time domain program to calculate the electric field distribution of the problem space is implemented using Matlab.

【Key words】FDTD；Matlab；Electric field distribution

0 引言

時(shí)域有限差分（Finite Difference Time Domain，F(xiàn)DTD）法是K.S.Yee在1966年給出的利用有限差分式把麥克斯韋（Maxwell）旋度方程替換為一組差分方程[1]，并提供所解問題中電磁特性物理意義的算法，可直接在時(shí)域中求解。

Matlab是一種功能強(qiáng)大、高效的高級(jí)技術(shù)計(jì)算語言和交互式環(huán)境[2]，在科學(xué)和工程領(lǐng)域中贏得了極為廣泛的應(yīng)用，將其用于FDTD法的數(shù)值計(jì)算及仿真，不僅可以簡化程序設(shè)計(jì)、操作方便，另外運(yùn)算結(jié)果也更簡潔。

因此，本文將結(jié)合Matlab強(qiáng)大的數(shù)組運(yùn)算和繪圖功能，通過對(duì)FDTD法編程來模擬出二維問題空間中電場分布。

1 FDTD法的基本原理

FDTD算法將問題空間離散為電場和磁場分量在其位置上交叉放置的空間網(wǎng)格點(diǎn)，并以中心差分的方式近似Maxwell方程中關(guān)于空間和時(shí)間的導(dǎo)數(shù)，通過時(shí)間向前推進(jìn)的差分方程模擬出電磁場在時(shí)域的進(jìn)程?？臻g網(wǎng)格中，電場分量位于Yee元胞網(wǎng)格單元每條棱的中心，磁場分量位于網(wǎng)格單元每個(gè)面的中心[3]，如圖1所示。

1.1 Maxwell方程的差分形式

Maxwell旋度方程為：

？犖×H=+J；？犖×E=--Jm（1）

已知本構(gòu)關(guān)系表達(dá)式為：

D=？著E；B=？滋H；J=？啄E；Jm=？啄mH

在直角坐標(biāo)系中，根據(jù)本構(gòu)關(guān)系把（1）式寫為：

（2）

下面我們求解（2）式的中心差分，令f（x，y，z，t）表示E或H某一分量，離散形式寫為：

f（x，y，z，t）=f（i？駐x，j？駐y，k？駐z，n？駐t）=fn（i，j，k）（3）

在二維問題空間中，假定任意的電磁場分量只與x，y坐標(biāo)有關(guān)，與z坐標(biāo)無關(guān)，即？墜/？墜z=0，以TE波為例，Hx=Hy=Ez=0，由（2）式可得

（4）

用中心差分式來近似（4）式中的導(dǎo)數(shù)，根據(jù)場分量的位置，并采??？駐x=？駐y=？駐z=？啄離散方式，得到了關(guān)于TE波的FDTD公式為：

Ex=CA·Ex+CB'·Hz-Hz（5）

Ey=CA·Ey+CB'·Hz-Hz（6）

（7）

式中系數(shù)CA，CB'，CQ'的定義為：

CA=1-/1+；CB'=/1+

CP=1-/1+；CQ'=/1+

為了統(tǒng)一TE波、TM波兩者方程的離散形式，分別將（5）、（6）、（7）式中的空間位置標(biāo)號(hào)移動(dòng)1/2，時(shí)間移動(dòng)？駐t/2，以上式子分別寫為：

Ex=CA·Ex+CB'·Hz-Hz（8）

Ey=CA·Ey+CB'·Hz-Hz（9）

（10）

利用TE波與TM波之間的對(duì)偶關(guān)系，寫出通用于計(jì)算求解二維問題空間中TE波與TM波的FDTD程序。

1.2 數(shù)值色散及穩(wěn)定性條件

FDTD方法為場的行為提供了一種解，連續(xù)函數(shù)的導(dǎo)數(shù)有限差分近似給解引入了誤差，我們把用FDTD數(shù)值方法得到的相速與實(shí)際的相速之間的差別稱為數(shù)值色散。為保證結(jié)果準(zhǔn)確性，空間網(wǎng)格大小應(yīng)滿足？姿min≥10？駐，？駐=min（？駐x，？駐y，？駐z），？姿min表示媒質(zhì)空間中最小波長值，減小網(wǎng)格大小雖然會(huì)減小數(shù)值色散，但在計(jì)算中將會(huì)占用更多的內(nèi)存。

為保證數(shù)值計(jì)算穩(wěn)定性，根據(jù)Cournant穩(wěn)定條件，算法中的時(shí)間步長應(yīng)滿足：c？駐t≤

其中c=1/為介質(zhì)中的光速，一般選??？駐t=δ/（2c）。

2 問題空間中電場分布的模擬

以圖2所示二維問題為例，幾何圖形中包含半徑為0.2m、介電常數(shù)為4的圓柱，激勵(lì)信號(hào)頻率為1GHz。此問題空間由邊長為5mm的正方形網(wǎng)格構(gòu)成，端接為8mm的PML，圓柱與PML邊界問題的空氣隙在xn、yn和yp方向上為30個(gè)網(wǎng)格，在xp方向上為80個(gè)網(wǎng)格，其激勵(lì)源為一正弦波形的外加電流密度。

其中，J1為激勵(lì)線源，E1為電場取樣點(diǎn)。此例中，定義了兩種輸出類型，包括瞬態(tài)電場分布和某頻率下的電場分布。通過Matlab編程來實(shí)現(xiàn)問題的定義與模擬，其主要編程如下：

（1）定義幾何體。設(shè)置圓柱的位置范圍及材料類型，包括中心坐標(biāo)、半徑及電介質(zhì)。

（2）定義激勵(lì)源。設(shè)置激勵(lì)源波形及參量，包括初始值矢量坐標(biāo)、波形類型及量級(jí)。

（3）定義并初始化輸出參量。設(shè)置頻率邊界，定義并初始化取樣電場、磁場及瞬態(tài)電場、某頻率下的電場，取更新頻率為10個(gè)時(shí)間步。

（4）顯示取樣參數(shù)，計(jì)算并獲取在節(jié)點(diǎn)上的電場。

（5）獲取取樣電場給定頻率的頻域響應(yīng)，規(guī)定此程序中場量在6000個(gè)時(shí)間步以后獲取。

（6）顯示頻率邊界輸出，計(jì)算得到的響應(yīng)曲線。

使用Matlab運(yùn)行仿真后，仿真結(jié)果如下圖所示：

假定正弦激勵(lì)的時(shí)域響應(yīng)在6000時(shí)間步達(dá)到穩(wěn)定后獲得穩(wěn)態(tài)場的幅度，所以用給定的程序能夠獲得場頻域響應(yīng)的幅度，而且場的幅度響應(yīng)僅屬于此單一的激勵(lì)頻率。由圖5發(fā)現(xiàn)，某些頻率的取樣電場幅度和相位響應(yīng)并不能同時(shí)求解，只能顯示1GHz單一頻率幅度。為解決此問題，我們將采用更有效地離散傅里葉變換（DFT）[4]。

如前面所討論的，我們將求解多頻率下電場幅度，如果激勵(lì)信號(hào)是具有一定頻譜寬度的波形，那么使用DFT就可以得到多個(gè)頻率的結(jié)果。首先，對(duì)于激勵(lì)波形，其頻譜應(yīng)包含所想要的頻率分量。其次，為計(jì)算多個(gè)頻率的場分布，必須實(shí)施實(shí)時(shí)DFT。使用DFT技術(shù)同時(shí)求解出多頻率取樣電場的幅度和相位響應(yīng)如圖8所示。另外，如圖9所示求解出的在1GHz頻率下的電場幅度分布與圖6結(jié)果相同，同時(shí)也求解出在2GHz頻率下的電場幅度分布，如圖10所示，驗(yàn)證了此方法的可行性。