求解三維電磁場問題的改進的自適應區域分解算法

李 鵬， 孟令琴

(1.上海大學特種光纖與光接入網省部共建重點實驗室，上海200072;2.東南大學毫米波國家重點實驗室，南京210096)

很多電磁問題用理論分析來計算十分復雜，甚至不可能得到結果，因此需要通過數值計算得到問題的結果.1966年，Yee［1］首次提出一種新的電磁場數值計算方法——時域有限差分(finite difference time domain，FDTD)算法.該算法已廣泛應用于電磁散射、電磁兼容、天線、波導與諧振腔系統、微波集成電路等問題的研究中.FDTD算法是以差分原理為基礎，用離散點上的數值解來逼近連續場區域的真實解.隨著研究的深入，各種提高FDTD算法的精度和速度的技術不斷涌現.Xu等［2-4］提出了一種新的區域分解算法，該算法首先把計算區域分割成幾個子區域，在每個子區域里面采用適合本區域的共形網格，然后分別計算各個子區域，最后在各個子區域邊界上采取一種數據傳遞方案，綜合可得整個區域的解.文獻［5］將拓撲網絡應用到FDTD算法中，把一個大問題區域分解成為幾個相互獨立的子區域，分別計算各個區域的散射參數，然后通過矩陣方程計算出總的散射參數，從而計算大區域上的電磁耦合和傳輸問題.文獻［6］提出并驗證了一種空間區域劃分FDTD算法.該算法把一個大的計算區域分解為幾個小的區域，并重復利用電腦有限的內存，從而可以解決大的電尺寸問題.然而這種方法會占用很大的硬盤空間和很長的計算時間.為了計算稀疏填充的大型多物體，文獻［7-8］提出了一種多區域的FDTD(multiple region/FDTD，MR/FDTD)算法，達到了減少內存需要量和縮短求解時間的目的.文獻［9-10］提出并驗證了一種自適應區域分解(adaptive domain decomposition-FDTD，ADD-FDTD)算法.該算法將三維電磁場問題的求解區域分解為若干獨立的子區域，在各個子區域之間通過連接邊界條件進行自適應檢測，以判斷波是否能傳播到子區域，以及子區域是否需要進行運算.該算法可以顯著縮短計算時間.

本研究針對文獻［9-10］中的自適應區域分解方法，提出了一種改進的自適應區域分解算法.在判斷波是否傳播到檢測面的時候，不采用某點的電場值和預先設定的閾值比較，而是通過對檢測面上的電場積分得到相應的電壓，然后選取其中最大的一個電壓值和預先設定的閾值比較.

隨著半導體材料科學和通信電子技術的快速發展，20世紀60年代出現了微波集成電路(hybrid microwave integrated circuit，HMIC)，之后又出現了單片微波集成電路(monolithic microwave integrated circuit，MMIC)、三維微波集成電路，最后出現了微波多芯片組件(microwave multi-chip module，MMCM).由于現代通信的需要，微波電路工作頻率不斷提高，這使得MMIC的互連和封裝越來越困難.互連是多芯片組件(multi-chip module，MCM)中的最關鍵技術之一，研究互連的特性問題可以提高集成電路集成度、速度和可靠性.對高速電路互連結構的分析過程中，采用頻域分析方法計算所要求的寬頻上的特性是十分耗時和繁瑣的，而時域方法可以通過一次計算得到寬頻上的頻率響應.

為了驗證所提出的算法在高速互連結構中應用的正確性，本研究選取了3個在高速電路互連結構中的常見結構，分別采用HFSS、原始FDTD算法、ADD-FDTD算法和IADD-FDTD算法進行了分析，證明了本研究提出的IADD-FDTD算法的正確性.

1 一種改進的ADD-FDTD算法

原始的ADD-FDTD算法中，各子區域之間在進行自適應區域檢測時，通過選擇相鄰子區域邊界面上的一點的電場值和預先設定的閾值進行比較.如微帶結構中，取微帶導帶下方的中心線上一點的電場值作為檢測值，取電場激勵源最大幅值的0.001%作為閾值.但該算法有2個比較明顯的缺點.

(1)在某些結構中，如圖1所示的多條平行微帶線結構中，有N條平行的微帶線，原始的ADDFDTD算法只是簡單地選擇檢測面上其中一條微帶線導帶下的電場Ei作為檢測值.如果Ei在檢測面上的所有電場值中并不是最大的，就可能出現波其實已經傳播到此子區域，但是當前的計算并沒有包含此子區域，從而造成計算結果的錯誤.

圖1 多條平行微帶線結構Fig.1 Parallel multi-microstrips structure

(2)原始ADD-FDTD算法中，如果檢測面與不連續結構距離過大，則整個計算區域只能分成很少的幾個子區域，不能充分發揮ADD-FDTD算法通過劃分子區域來提高計算效率的作用.如果檢測面距離不連續結構很近，則檢測面上有的點的電場值可能存在不穩定的情況，這就會出現虛假的較大的電場值.在檢測時如果選擇檢測面上的最大電場值的點作為檢測點，則容易出現其實波并沒有傳播到檢測面上，但是由于某點不穩定而造成了較大的電場值，程序就會認為波已經傳播到了下一個子區域，從而將下一個子區域加入迭代，這樣就會削弱ADDFDTD算法運算時間短、運算效率高的優點.

針對原始ADD-FDTD算法中的不足，本研究提出一種改進的ADD-FDTD算法，其流程圖如圖2所示.從圖中可以看出，IADD-FDTD算法的步驟如下：首先，對整個區域初始化賦0值，根據實際計算的物理結構將計算區域分為Z0個區域，令Z的初始值等于1;然后，開始時間循環，對Z個子區域內的電場值和電壓值進行更新，其余Z0－Z個子區域的值不變.對于子區域之間邊界上的電場進行積分得到電壓，根據這些電壓值判斷電磁波是否傳播到了下一個子區域.如果沒傳播到，則直接開始下一個時間循環，如果傳播到了下一個子區域，則Z值加1以后再進入下一個時間循環.這樣一直循環下去，直到時間步進結束.

圖2 IADD-FDTD算法流程圖Fig.2 Flow chart of IADD-FDTD method

本算法在判斷波是否傳播到檢測面時，并不采用原始ADD-FDTD算法中的方法，即將某點的電場值和預先設定的閾值Eon比較，而是通過對檢測面上的電場積分得到相應的電壓，然后選取其中最大的一個電壓值和預先設定的閾值電壓進行比較.圖3是一個沿著z軸放置的微帶結構，ABCD面代表其中的一個檢測面.設結構沿x方向劃分了i+1個網格，那么在ABCD面上沿虛線積分可以得到i個電壓值V1，V2，…，Vi，選擇其中最大的值Vx和閾值Von相比較，如果Vx＞Von，則認為波已經傳播到了檢測面;如果Vx＜Von，則認為波還未傳播到檢測面.閾值Von的選取是令Ey等于電場源的最大值的0.001%，然后乘以dy再乘以y方向上的網格數目J得到的.

圖3 檢測面示意圖Fig.3 Schematic drawing of test surface

首先，由于對檢測面上的1，2，…，i的所有Ey分別進行積分得到V1，V2，…，Vi，然后取其中的最大值和Von進行比較，這樣在計算有多條平行微帶線結構時，就不會出現選擇的點不是邊界面上的最大場值點的問題，從而確保了不會出現波已經傳播到下一子區域而計算中卻沒有包含下一子區域的情況.

其次，由于采用了電壓檢測，即使某點電場由于不穩定而造成了較大的電場值，該電場值在整個積分中的影響也很小.也就不會出現程序由于錯誤的最大電場值，而認為波已經傳播到了下一個子區域，從而將下一個子區域加入迭代的情況.這表明本研究提出的算法具備一定的抗干擾性.另外，由于檢測面和不連續結構不需要相隔太遠，因此可以對計算區域進行更為密集的子區域劃分.

2 數值計算

下面舉例對高速電路中常見的微帶結構進行分析，包括三維微帶結構和多層微帶電路互連結構(單通孔結構和雙通孔結構)，以證明本研究提出的IADD-FDTD算法的正確性和有效性.

2.1 微帶線結構的分析

在高速微波電路中，微帶線是應用最廣泛的微波信號傳輸線.本算例采用FR4材料的PCB板，尺寸為20 mm×20 mm，介電常數εr=4.5，介質高h= 0.7 mm，微帶線寬W=1.2 mm.具體結構如圖4所示.

圖4 微帶互連線Fig.4 Microstrip interconnection line

為了將數值色散減少到可以接受的水平，空間網格步長需要滿足［11］：

式中，Δhmax表示Δx，Δy和Δz的最大值，λmin為計算頻率范圍內波長的最小值.綜合考慮所計算結構的尺寸，這里取Δx=Δy=Δz=0.1 mm.

空間和時間離散間隔應該滿足Courant穩定性條件［11］：

式中，vmax為工作模式的最大相速.

根據以上條件，并考慮到計算頻率和合理的迭代次數，選擇Δt=0.18×10－12s.

在導帶下面的切面上均勻設置激勵源，即

式中，t0=6T，T=25Δt.

微帶地板和導帶均設置成理想金屬導體(perfect electric conductor，PEC)，在介質板上方設置10層空氣層，四周包裹單軸各向異性介質完全匹配層(uniaxial perfectly matched layer，UPML)，總網格數為220×220×27.首先，在Ansoft HFSS軟件中建模仿真微帶線結構;然后，分別用原始FDTD，ADDFDTD和IADD-FDTD 3種算法通過Matlab編程進行仿真.采用ADD-FDTD算法仿真時，將微帶結構沿z軸分成5個子區域，而采用IADD-FDTD算法仿真時，將微帶結構沿z軸分成8個子區域.圖5給出了采用4種方法仿真得到的散射參數曲線.

圖5 微帶線的頻率響應Fig.5 Frequency response of microstrip

由圖5可知，采用3種算法計算的結果完全相同，另外，3種算法和HFSS仿真得到的結果也比較吻合.ADD-FDTD算法比原始FDTD算法節省約6.80%的計算時間，而本研究提出的IADD-FDTD算法比原始FDTD算法節省約10.98%的計算時間.

2.2 多層印刷電路板結構中微帶單通孔互連結構的分析

隨著電路小型化的需求，印刷電路板(printed circuit board，PCB)更多地采用多層結構，這樣就可以將傳輸線放在不同的層，在不同走線之間用接地板進行屏蔽以減少相互間的干擾.這些不同層的傳輸線之間通過通孔連接，由于電路的結構越來越復雜、傳輸速率越來越快、工作頻率越來越高，使得由于互連和封裝所引起的寄生效應在實際電路設計過程中變得越來越不容忽視，其中比較突出的問題就是多芯片結構的垂直通孔互連的分析.

本算例的具體結構如圖6所示，其中 L= 25.2 mm，W=13.9 mm，T=3.2 mm，Wm=3.3 mm，D1=0.7 mm，D2=3.9 mm，D3=3.9 mm，介質的介電常數為3.4.

圖6 單通孔互連結構Fig.6 Interconnection structure of single through-hole

在進行空間網格劃分時，x方向和y方向的空間步長均為0.1 mm，z方向的空間步長為0.2 mm;時間步長Δt=0.1×10－12s.激勵源和2.1節中的算例相同，將激勵源加在10Δx處微帶線下方的平面上.波出射面選擇在距離右邊界10Δx的平面上.

在計算區域周圍采用一階Mur吸收邊界條件，總劃分網格為252×139×32.由于原始ADD-FDTD算法采用電場檢測，要求檢測邊界與不連續結構的距離較遠，所以只能劃分成4個區域.而 IADDFDTD算法由于采用電壓檢測，不需要檢測邊界和不連續結構相隔太遠，所以在程序中可以劃分6個區域.圖7給出了采用4種算法仿真得到的0～20 GHz的散射參數曲線.

由圖7可知，原始FDTD，ADD-FDTD和IADDFDTD 3種算法的計算結果完全相同，另外，3種算法和HFSS仿真得到的結果也比較吻合，從而驗證了本研究提出的算法的正確性.ADD-FDTD算法比傳統FDTD算法節省約5.88%的計算時間，本研究提出的IADD-FDTD算法比傳統FDTD算法節省約10.61%的計算時間.

圖7 單通孔的頻率響應Fig.7 Frequency response of single through-hole

2.3 多層PCB結構中雙通孔互連結構的分析

在多層PCB結構中往往采用多個通孔對微帶進行互連，多個通孔之間會產生干擾，本算例以2個通孔為例進行分析.由于FDTD算法中的基本計算單元是矩形網格，為了計算簡便，本算例用方形柱來代替實際應用中采用的圓柱.采用的板材為FR4材料，其介電常數εr=4.5，PCB板尺寸為30 mm× 17.6 mm，此結構共有2層，每層高為0.8 mm.矩形方柱的尺寸為0.8 mm×0.8 mm×1.6 mm，地層上的矩形反焊盤尺寸為1.6 mm×1.6 mm，頂層和底層上的所有微帶線長度均為 10.8 mm，寬度為1.6 mm，具體結構和尺寸如圖8所示.

圖8 雙通孔互連結構Fig.8 Interconnection structure of double through-holes

對該結構進行空間網格劃分時，x，y，z方向上的空間步長均為0.2 mm;時間步長 Δt=0.385× 10－12s.激勵源和2.1節中算例相同，將激勵源加在10Δy處微帶線下方的平面上.

在z方向上，器件上下包裹10層空氣層，四周各設置10層UPML吸收邊界條件.總的網格數為170×108×48.分別采用HFSS，FDTD，ADD-FDTD和IADD-FDTD算法對該結構進行仿真.采用ADDFDTD算法時將待解問題分成4個區域，采用IADDFDTD算法時將待解問題分成6個區域.圖9給出了4種算法仿真得到的0～20 GHz的散射參數曲線.

圖9 雙通孔結構的頻率響應Fig.9 Frequency response of double through-holes

由圖9可知，FDTD，ADD-FDTD和IADD-FDTD 3種算法的計算結果完全相同，另外，3種算法和HFSS仿真得到的結果也比較吻合，從而驗證了本研究提出的算法的正確性.ADD-FDTD算法比傳統FDTD算法節省約7.12%的計算時間，本研究提出的IADD-FDTD算法比傳統 FDTD算法節省約12.67%的計算時間.

3 結束語

相對于傳統FDTD和ADD-FDTD算法，本研究提出的IADD-FDTD算法消除了不穩定點對程序的影響，進一步提高了運算速度.對高速互連電路中常見的3種結構，分別用原始FDTD，ADD-FDTD和IADD-FDTD算法進行了Matlab編程仿真.3種算法的計算結果完全一致，并且和HFSS仿真得到的結果也非常吻合.

由于本研究提出的算法可以減少檢測面上不穩定點造成的影響，從而可以進行更密集的分區，這也使得本研究提出的算法更能節省計算時間.這個優點在計算復雜結構和電的大尺寸問題時表現得更為明顯，可以節省更多的計算時間.

另外，由于本研究提出的算法采用了分區域檢測的方法，因此可以在不同的子區域內分別建立不同的計算模型，提高了解決復雜問題的精度.這對于計算區域較大、集成度較高的高速電路也是非常適合的.

［1］ YEEK S.Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media［J］.IEEE Transactions on Antennas Propagation，1966，14(5)：302-307.

［2］ XUF，HONGW.Domain decomposition FDTD algorithm for the analysis of a new type of E-plane sectorial horn with aperture field distribution optimization［J］.IEEE Transactions on Antenna and Propagation，2004，52 (2)：426-434.

［3］ XUF，HONGW，ZHUX W.A new time domain near field to farfield transformation forFDTD in two dimension［C］∥ IEEE MTT-S International Microwave Symposium Digest.2002：2057-2060.

［4］ XUF，HONGW.Accurate analysis of microstrip patch antennas using FDTD method［C］∥ International Conference on Computational Electromagnectics and Its Applications Proceedings.1999：214-217.

［5］ WANGJ Y，GAOB Q.A hybrid method on calculating EM coupling in large problem space［C］∥International Symposium on Electromagnetic Compatibility Proceedings.1997：442-445.

［6］ XUEZ H，GAOB Q，LIUR X，et al.Space-divided FDTD method［C］∥ International Symposium on Electromagnetic Compatibility Proceedings.1999：568-571.

［7］ JOHNSONJ M，YAHYAR S.Multiple region FDTD (MR/FDTD)and its application to microwave analysis and modeling［C］∥IEEE MTT-S Symposium Digest.1996：1475-1478.

［8］ JOHNSONJ M，YAHYAR S.MR/FDTD：a multipleregion finite difference time domain method［J］.Microwave and Optical Technology Letters，1997，14 (2)：101-105.

［9］ 張華，洪偉，郝張成，等.一種采用邊界檢測的自適應區域分解時域有限差分方法［J］.應用科學學報，2005，23(5)：493-496.

［10］ 張華，洪偉，郝張成.求解三維電磁問題的自適應區域分解FDTD方法［J］.電波科學學報，2006，21(3)：397-402.

［11］ 葛德彪，閆玉波.電磁波時域有限差分法［M］.西安：西安電子科技大學出版社，2002：25-31.