時域有限差分法的Matlab仿真

張通　孫晶

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

【關鍵詞】時域有限差分法；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 引言

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

Matlab是一種功能強大、高效的高級技術計算語言和交互式環境[2]，在科學和工程領域中贏得了極為廣泛的應用，將其用于FDTD法的數值計算及仿真，不僅可以簡化程序設計、操作方便，另外運算結果也更簡潔。

因此，本文將結合Matlab強大的數組運算和繪圖功能，通過對FDTD法編程來模擬出二維問題空間中電場分布。

1 FDTD法的基本原理

FDTD算法將問題空間離散為電場和磁場分量在其位置上交叉放置的空間網格點，并以中心差分的方式近似Maxwell方程中關于空間和時間的導數，通過時間向前推進的差分方程模擬出電磁場在時域的進程。空間網格中，電場分量位于Yee元胞網格單元每條棱的中心，磁場分量位于網格單元每個面的中心[3]，如圖1所示。

1.1 Maxwell方程的差分形式

Maxwell旋度方程為：

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

已知本構關系表達式為：

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

在直角坐標系中，根據本構關系把（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坐標有關，與z坐標無關，即？墜/？墜z=0，以TE波為例，Hx=Hy=Ez=0，由（2）式可得

（4）

用中心差分式來近似（4）式中的導數，根據場分量的位置，并采取？駐x=？駐y=？駐z=？啄離散方式，得到了關于TE波的FDTD公式為：

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

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

（7）

式中系數CA，CB'，CQ'的定義為：

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

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

為了統一TE波、TM波兩者方程的離散形式，分別將（5）、（6）、（7）式中的空間位置標號移動1/2，時間移動？駐t/2，以上式子分別寫為：

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

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

（10）

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

1.2 數值色散及穩定性條件

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

為保證數值計算穩定性，根據Cournant穩定條件，算法中的時間步長應滿足：c？駐t≤

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

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

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

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

（1）定義幾何體。設置圓柱的位置范圍及材料類型，包括中心坐標、半徑及電介質。

（2）定義激勵源。設置激勵源波形及參量，包括初始值矢量坐標、波形類型及量級。

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

（4）顯示取樣參數，計算并獲取在節點上的電場。

（5）獲取取樣電場給定頻率的頻域響應，規定此程序中場量在6000個時間步以后獲取。

（6）顯示頻率邊界輸出，計算得到的響應曲線。

使用Matlab運行仿真后，仿真結果如下圖所示：

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

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