共形HIE-FDTD 方法的穩(wěn)定性分析及改進

戴婧怡，楊 陽，李強兵，鮑卓如

(南京航空航天大學電子信息工程學院，江蘇南京210016)

0 引言

1966年由K S Yee［1］提出的傳統(tǒng)時域有限差分方法，已廣泛應用于各種電磁仿真問題。為了提高計算速度、改善仿真效果和節(jié)省計算機資源，近年來，國內(nèi)外學者提出了各種算法優(yōu)化策略。然而FDTD方法本身存在著2個缺陷:①階梯法模擬不能與網(wǎng)格共形的目標時，不可避免要進行階梯近似，從而產(chǎn)生誤差并可能會引起虛擬表面波;②柯朗－弗里德里希斯列維 (Courant－Friedrichs－Lewy，CFL)［2－3］穩(wěn)定條件限制了時間步長的選取，使得計算效率受到限制。

對于階梯近似的問題，Dey和Mittra［4］提出的金屬共形時域有限差分方法，可以很好地改善由于階梯近似產(chǎn)生的誤差。但是，該方法并非沒有弊端，在運用該方法時，需要滿足2個穩(wěn)定性條件，而且即便在滿足條件的前提下，也不能做到時間步長的選取完全達到CFL條件?？梢哉f該方法是以增加計算時間來換取較高的精度。對于CFL條件限制的改進，可以采用無條件穩(wěn)定或弱條件穩(wěn)定的方法，其中以交替方向隱式時域有限差分方法(alternating－direction implicit FDTD，ADI－FDTD)［5］和Crank－Nicolson時域有限差分方法(Crank－Nicolson FDTD，CNFDTD)［6－7］為主。1999年T Namiki提出的ADI－FDTD方法，將傳統(tǒng)的時域有限差分方法中電場顯式的迭代改成隱式迭代，這樣時間步長的選取就不再受到CFL的限制，但是該方法的計算精度會隨著時間步長的增加而降低。CN－FDTD方法在時間步長增大時離散誤差較小，精度較高，但該方法需要求解一個大型稀疏矩陣方程組，在求解微細結(jié)構時效率很低。為了提高計算精度和效率，基于隱格式思想的 HIE－FDTD 方法［8－9］被提出來。該方法在電場迭代時，任意2個方向均采用隱式，另外一個方向采用顯式，這樣時間步長的取值可以與顯式方向無關，即弱條件穩(wěn)定。HIE－FDTD方法比ADI－FDTD方法有更高的精度，特別適用于在某個方向有微細剖分的結(jié)構。本文提出的金屬共形HIE－FDTD方法采用x，y方向的隱式差分，著重分析在共形條件下減小HIE－FDTD時間步長的因子［10］，并給出修正方法。

1 Dey-Mittra的金屬共形方法應用在HIE-FDTD

HIE－FDTD方法在1999年被提出，是基于無條件穩(wěn)定思想的改進算法。算法采用電場分量在任意2個方向進行隱式差分，磁場分量差分格式不變的方法，可將CFL條件放寬，在第2節(jié)中進行推導［11］。

(1)—(2)式中:c為真空中光速;d為維數(shù)(對于二維問題，d=2;對于三維問題，d=3)。

本文分析的HIE－FDTD方法在x，y方向采用隱格式，z方向采用顯式的格式。下面給出x方向和z方向E，H的離散方程為

同理，可以寫出y方向電場與磁場的迭代公式。綜上，采用z方向顯格式x，y方向隱格式的HIE－FDTD方法的迭代公式，可以寫成矩陣(7)的形式

當用HIE－FDTD方法模擬不規(guī)則金屬邊緣時，需要模擬的結(jié)構邊緣不能與網(wǎng)格線重合，含有不規(guī)則邊緣的金屬板用階梯網(wǎng)格劃分邊緣的截面示意圖如圖1a所示，其中，深色部分為金屬，淺色部分為介質(zhì)，syz為共形網(wǎng)格中介質(zhì)部分的面積，ly，lz為共形網(wǎng)格中介質(zhì)部分的邊長。

圖1 含有不規(guī)則邊緣的金屬片用階梯網(wǎng)格劃分邊緣的截面示意圖Fig.1 Cross－sectional view of a PEC boundary and the distorted cells

用Dey－Mittra提出的金屬共形算法只需要對磁場表達式做如下修正。若含不規(guī)則邊緣的截面只在某一個平面，則只有一個方向的磁場需要修正，以x方向為例，磁場表達式由(5)修正為

Dey－Mittra的金屬共形方法需要滿足2個條件以保持算法的基本穩(wěn)定性［12］。

1)共形網(wǎng)格面積s(i，j，k)應大于規(guī)則網(wǎng)格面積δ2的5%。

2)共形網(wǎng)格中的最長相對邊長l'(i，j，k)與該網(wǎng)格相對面積 s'(i，j，k)的比率小于12。

當滿足上述條件，Δt也需要適當減小，以保證穩(wěn)定。因此，將該共形方法應用在HIE－FDTD方法中，也會導致時間步長的減小，或發(fā)散的提前。下面就導致發(fā)散提前的不穩(wěn)定因子進行分析，并給出改進方法，進行數(shù)值驗證。

2 金屬共形HIE-FDTD方法的穩(wěn)定性分析

HIE－FDTD方法的時間步長的選取只與隱格式方向的空間步長有關，下面進行簡要的推導。

參考HIE－FDTD方法的時域基本推進公式的矩陣形式(7)，考慮平面波本征模的解以及有限差分近似，即

(9)式中:f(i，j，k)=e0xyz;ku為波矢量;ζ為增長因子;u=x，y，z;V=E，H 帶入矩陣(7)式得到新的矩陣形式為

(11)式中，ru=(cΔt)2W2u。

解得

為了保證數(shù)值結(jié)果的穩(wěn)定性，要求增長因子|ζi|≤1，即

經(jīng)過推導可以得出:HIE－FDTD方法的數(shù)值穩(wěn)定性條件為

將變形后的修正磁場表達式(15)帶入矩陣(7)得到

為使該矩陣方程有非零解，其系數(shù)行列式必須為零，因此，我們得到與增長因子ζ有關的多項式為

K是由同一個網(wǎng)格中的介質(zhì)與金屬的分界產(chǎn)生的，K的大小取決于小邊長與小面積的比值。當網(wǎng)格為規(guī)則網(wǎng)格時，K≡1，(16)式即還原為常規(guī)HIEFDTD方程，即時間步長的選取只需滿足(14)式。

Dey－Mittra的方法中，在圖1b的情況下，有

在圖1c的情況下，有即，在Dey－Mittra方法中，需要分情況討論且不能保證穩(wěn)定性。要使共形方法盡量不破壞原迭代方程的穩(wěn)定性，則需要因子K趨近于1。

根據(jù)K的計算公式，要使得K≡1，可令

則

容易求得:syz(Dey－Mittra)≤sequ。更容易滿足穩(wěn)定條件1)。根據(jù)(20)式，可以看出本文提出的方法較Dey－Mittra的方法可以保持HIE－FDTD固有的穩(wěn)定條件。

3 數(shù)值證明

一個功分器結(jié)構［8］如圖2所示。器件尺寸分別為 W50=4.5 mm，L50=10 mm，d=4.3 mm，hl=4.73 mm，Wl=28.4 mm，ll=11.92 mm，介質(zhì)板厚度為1.575 mm，相對介電常數(shù) εr=2.2。取dx=0.787 5 mm，dy=1.5 mm，dz=0.1 mm。圖3為端口2的時域波形，圖3上的曲線分別是應用Dey－Mittra 的方法分別取 dt=2.311 3 ps=7ΔtHIE－FDTD，dt=1.650 9 ps=5ΔtHIE－FDTD，dt=1.320 7 ps=4ΔtHIE－FDTD，以及應用本文提出的方法，取 dt=2.311 3 ps=7ΔtHIE－FDTD的端口時域波形，從圖3中可以得出結(jié)論:①減小時間步長可有效地推遲發(fā)散;②本文提出的方法，較Dey－Mittra的方法，可以大大推遲發(fā)散，在時間步長取相同的情況下，可將發(fā)散出現(xiàn)的時間步推遲3.18倍，(即，使時間步長取 Dey－Mittra方法的1.75倍，發(fā)散可推遲1.25倍，時間步長取Dey－Mittra方法的1.4倍，發(fā)散可推遲1.4倍)。圖4給出了本文提出的算法在取 dt=2.311 3 ps=7ΔtHIE－FDTD以及 Dey－Mittra 算法取 dt=1.320 7 ps=4ΔtHIE－FDTD時端口2的S21和仿真結(jié)果相比較的對比圖?？梢钥闯?，本文提出的方法具有較高的精度。表1給出了本文提出的算法取dt=2.311 3 ps=7ΔtHIE－FDTD以及 Dey－Mittra 算法取 dt=1.320 7 ps=4ΔtHIE－FDTD時，CPU time 的比較，從表1 可以看出，時間節(jié)約878 s(約42%)。結(jié)果證明:本文提出算法，能有效推遲發(fā)散，并具有較高的精度和效率。

表1 Dey－Mittra 方法取 dt=1.320 7 ps=4ΔtHIE－FDTD時和本文提出方法取 dt=2.311 3 ps=7ΔtHIE－FDTD時的 CPU 用時比較Tab.1 CPU time of Dey－Mittra when dt=1.320 7 ps=4ΔtHIE－FDTD and proposed method when dt=2.311 3 ps=7ΔtHIE－FDTD

圖2 功分器結(jié)構Fig.2 Geometry of the conventional TPD

圖3 端口2的時域波形圖Fig.3 Wave shape of time domain electric field

4 結(jié)論

本文將Dey－Mittra的金屬共形方法應用在HIEFDTD，并分析了該金屬共形HIE－FDTD方法不能取到CFL條件的原因，提出了一個可以增加算法穩(wěn)定性的面積計算方法。數(shù)值實驗結(jié)果證明:提出的方法易于實現(xiàn)，當2種方法取相同時間步長，本文提出方法較Dey－Mittra的算法發(fā)散出現(xiàn)時刻推遲3.18倍，在相同精度情況下，計算時間可以節(jié)約42%，是一種高效的計算方法。

圖4 Dey－Mittra共形算法、本文提出算法和文獻結(jié)果的對比圖Fig.4 Comparation of S21with the Dey－Mittra’s method，the proposed method and the reference result

［1］YEE Kane S.Numerical solution of initial boundary value problems involving Maxwell＇s equations in isotropic media［J］.IEEE Transaction on Antennas and Propagation，1966，14(3):302－307.

［2］TAFLOVE A，HAQNESS S C.Computational Electromagnetics:The finite－difference time－domain method ［M］.Boston，MA:Artech House，2000.

［3］童創(chuàng)明.計算電磁學快速方法［M］.西安:西北工業(yè)大學出版社，2010.TONG Chuangming.Fast method for computational electromagnetics［M］.Xi＇an:Northwest Polytechnic University Press，2010.

［4］DEY Supriyo ，MITTRA Raj.A locally conformal finite difference time－domain(FDTD)algorithm for modeling three－dimensional perfectly conducting objects［J］.IEEE Microwave and Guided Wave Letters，1997，7(9):273－275.

［5］NAMIKI Takefumi.A new FDTD algorithm based on alternating－direction implicit method ［J］.IEEE Transactions on Microwave theory and technology，1999，47(10):2003－2007.

［6］YANG Y，CHEN R S，WANG D X，et al.Unconditionally stable Crank－Nicolson FDTD method for simulation of three－dimensional microwave circuits ［J］.IET Microw.Antennas Propag，2007，1(4):937－942.

［7］YANG Y，F(xiàn)AN Z H，DING D Z，et al.Extend two－step preconditioning technique for the Crank－Nicolson finitedifference time－domain method to analyze the 3D microwave circuits［J］.Int J RF Microw C E，2009，19(4):460－469.

［8］van VIEN，CHAUDHURI Sujeet K.A hybrid implicit－ex－plicit FDTD scheme for nonlinear optical waveguide modeling［J］.IEEE Transactions on Microwave Theory and Techniques，1999，47(5):540－545.

［9］HUANG B K，WANG G，JIANG Y S，et al.A Hybrid Implicit－Explicit FDTD scheme with weakly conditional stability［J］.Microwave and Optical Technology Letters，2003，39(2):97－101.

［10］DAI Jian，CHEN Zhizhang，SU Donglin，et al.Stability Analysis and Improvement of the Conformal ADI－FDTD Methods［J］.IEEE Transactions on Antennas and Propagation，2011，59(6):2248－2258.

［11］蘭婧，楊陽，戴婧怡.FDTD全波分析缺陷地結(jié)構微帶線的高效算法實現(xiàn)［J］.重慶郵電大學學報:自然科學版，2012，24(4):426－430，466.LAN Jing，YANG yang，DAI Jingyi.Highly－efficient FDTD method of full－wave analysis on the microstrip defected ground structure［J］.Journal of Chongqing University of Posts and Telecommunications:Natural Science Edition，2012，24(4):426－430，466.

［12］葛德彪，閆玉波.時域有限差分方法［M］.2版.西安:西安電子科技大學出版社，2002:184－185.GE Debiao，YAN Yubo.Finite－Difference Time－Domain Method for Electromagnetic Waves［M］.2nd ed.Xi＇an:Xi dian Press，2006:184－185.