戴婧怡，楊 陽(yáng)，李強(qiáng)兵，鮑卓如

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

0 引言

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

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

1 Dey-Mittra的金屬共形方法應(yīng)用在HIE-FDTD

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

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

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

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

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

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

用Dey－Mittra提出的金屬共形算法只需要對(duì)磁場(chǎng)表達(dá)式做如下修正。若含不規(guī)則邊緣的截面只在某一個(gè)平面，則只有一個(gè)方向的磁場(chǎng)需要修正，以x方向?yàn)槔?，磁?chǎng)表達(dá)式由(5)修正為

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

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

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

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

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

HIE－FDTD方法的時(shí)間步長(zhǎng)的選取只與隱格式方向的空間步長(zhǎng)有關(guān)，下面進(jìn)行簡(jiǎn)要的推導(dǎo)。

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

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

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

解得

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

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

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

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

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

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

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

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

則

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

3 數(shù)值證明

一個(gè)功分器結(jié)構(gòu)［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，相對(duì)介電常數(shù) εr=2.2。取dx=0.787 5 mm，dy=1.5 mm，dz=0.1 mm。圖3為端口2的時(shí)域波形，圖3上的曲線分別是應(yīng)用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，以及應(yīng)用本文提出的方法，取 dt=2.311 3 ps=7ΔtHIE－FDTD的端口時(shí)域波形，從圖3中可以得出結(jié)論:①減小時(shí)間步長(zhǎng)可有效地推遲發(fā)散;②本文提出的方法，較Dey－Mittra的方法，可以大大推遲發(fā)散，在時(shí)間步長(zhǎng)取相同的情況下，可將發(fā)散出現(xiàn)的時(shí)間步推遲3.18倍，(即，使時(shí)間步長(zhǎng)取 Dey－Mittra方法的1.75倍，發(fā)散可推遲1.25倍，時(shí)間步長(zhǎng)取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時(shí)端口2的S21和仿真結(jié)果相比較的對(duì)比圖。可以看出，本文提出的方法具有較高的精度。表1給出了本文提出的算法取dt=2.311 3 ps=7ΔtHIE－FDTD以及 Dey－Mittra 算法取 dt=1.320 7 ps=4ΔtHIE－FDTD時(shí)，CPU time 的比較，從表1 可以看出，時(shí)間節(jié)約878 s(約42%)。結(jié)果證明:本文提出算法，能有效推遲發(fā)散，并具有較高的精度和效率。

表1 Dey－Mittra 方法取 dt=1.320 7 ps=4ΔtHIE－FDTD時(shí)和本文提出方法取 dt=2.311 3 ps=7ΔtHIE－FDTD時(shí)的 CPU 用時(shí)比較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é)構(gòu)Fig.2 Geometry of the conventional TPD

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

4 結(jié)論

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

圖4 Dey－Mittra共形算法、本文提出算法和文獻(xiàn)結(jié)果的對(duì)比圖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)明.計(jì)算電磁學(xué)快速方法［M］.西安:西北工業(yè)大學(xué)出版社，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］蘭婧，楊陽(yáng)，戴婧怡.FDTD全波分析缺陷地結(jié)構(gòu)微帶線的高效算法實(shí)現(xiàn)［J］.重慶郵電大學(xué)學(xué)報(bào):自然科學(xué)版，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］葛德彪，閆玉波.時(shí)域有限差分方法［M］.2版.西安:西安電子科技大學(xué)出版社，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.