陳 衛(wèi), 饒 元, 朱 軍, 傅雷揚
(安徽農(nóng)業(yè)大學(xué) 信息與計算機學(xué)院,安徽 合肥 230036)
時域有限差分方法[1](Finite Difference Time Doman,簡稱FDTD)作為一種經(jīng)典的時域電磁仿真方法,直接求解依賴時間的麥克斯韋旋度方程,利用二階精度的中心差分近似旋度方程中的時間及空間微分算符,且自動滿足無散條件和場量切向連續(xù)性條件,已成為計算與電磁波相關(guān)的各類問題的非常有效的方法[2]。時域有限差分法主要有兩大缺點:① 由于采用規(guī)則的結(jié)構(gòu)網(wǎng)格和階梯近似方法模擬介質(zhì)曲面和材質(zhì)不連續(xù)性物體,難以獲得滿意的結(jié)果;② 由于數(shù)值色散、各向異性的影響,會產(chǎn)生明顯的積累誤差,造成了計算結(jié)果的失真、歪曲。因此,為了獲得精確的結(jié)果,就需要進(jìn)行細(xì)網(wǎng)格劃分,從而造成很大的計算負(fù)擔(dān)。針對第1個問題,亞網(wǎng)格技術(shù)[2-3]、材質(zhì)平均技術(shù)[4-5]等已用于傳統(tǒng)的FDTD方法;對于第2個問題,文獻(xiàn)[6]提出了FDTD(2,4)方法,其粗網(wǎng)格劃分,大大節(jié)約了計算機資源。但是,在高階差分情況下精確模擬介質(zhì)曲面是一個值得研究的問題。為了解決這一問題,一些學(xué)者做了大量的工作,4階交錯差分的一側(cè)導(dǎo)數(shù)(One-Sided Difference)技術(shù)[7]和導(dǎo)數(shù)匹配技術(shù)(Derivative Matching)[8],其收斂率可達(dá)4階,但后時不穩(wěn)定現(xiàn)象和繁瑣的操作,使其無法推廣?;旌蟻喚W(wǎng)格(Hybrid-Subgridding)[9]技術(shù)的思想是在目標(biāo)邊界附近采用FDTD(2,2)方法,網(wǎng)格是細(xì)網(wǎng)格;遠(yuǎn)離邊界處采用高階FDTD(2,4)方法,網(wǎng)格是粗網(wǎng)格。由于差分精度、網(wǎng)格尺寸及穩(wěn)定度等因素的不同,其需要采用時間和空間的雙插值方案來減少粗細(xì)網(wǎng)格之間的非物理反射,這必然帶來計算精度的降低。
在前人研究的基礎(chǔ)上,本文提出了一種新的介質(zhì)體曲面問題的FDTD(2,4)建模方法,通過對原網(wǎng)格內(nèi)外環(huán)路所圍區(qū)域進(jìn)行二次剖分,分別得到內(nèi)外環(huán)路的平均電磁參數(shù),再通過加權(quán)平均得到等效電參數(shù),避免了積分運算,簡化了建模過程。FDTD(2,4)運算仍在原網(wǎng)格上進(jìn)行,所以沒有增加計算量。相關(guān)數(shù)值結(jié)果表明,在相同計算量的條件下,本文方法具有更高的精度。
高階差分本質(zhì)上符合廣義的安培定理和廣義的法拉第電磁感應(yīng)定理,只不過將與空間2階差分等價的單環(huán)路積分轉(zhuǎn)化成與空間4階差分等價的雙環(huán)路積分。
其中,εr1和εr2分別為區(qū)域S1和區(qū)域S2的局部相對介電常數(shù);σ1和σ2分別為區(qū)域S1和區(qū)域S2的局部電導(dǎo)率;S1為由內(nèi)環(huán)路L1圍繞的深灰色區(qū)域;S2為由外環(huán)路L2圍繞的淺灰色區(qū)域。
圖1 廣義安培環(huán)路
在對材質(zhì)不連續(xù)性處理前,本文作如下假設(shè):
(1)在內(nèi)環(huán)路或外環(huán)路的棱邊中點的磁場值等于沿著該條棱邊磁場的平均值。
(2)區(qū)域S1或區(qū)域S2中心的電場值等于電場在整個區(qū)域的平均值。
(3)區(qū)域S1或區(qū)域S2中心的介電常數(shù)與電導(dǎo)率分別等于區(qū)域內(nèi)的平均值。
首先,將時間連續(xù)的空間積分方程(1)、(2)轉(zhuǎn)化為時間連續(xù)的空間差分方程,即
為了保持和FDTD(2,4)形式上的一致,根據(jù)麥克斯韋方程的積分形式和微分形式的等價性,運算(3)式和(4)式可得:
對上述積分離散求和得:
其中,n1、n2分別為內(nèi)、外環(huán)路劃分子網(wǎng)格的數(shù)量,通常取n2=3n1。
最后,將(7)式在時間上離散化,便可得到^Ex場的迭代公式為:
顯然,(11)式形式上和一般的FDTD(2,4)保持一致,平均電磁參數(shù)計算可在整個電磁仿真的預(yù)處理環(huán)節(jié)執(zhí)行,同主程序相比,能在很短的時間內(nèi)完成。
算例1 無耗介質(zhì)球的頻域遠(yuǎn)場分析。
球半徑為0.5m,相對介電常數(shù)為4。入射波沿z方向傳播,x方向極化,頻率300MHz。一致空間步長設(shè)置為:Δδ=0.05m,穩(wěn)定度常數(shù)CFLδ=0.4。階梯近似 FDTD(2,4)和高階子網(wǎng)格材質(zhì)平均 AM-FDTD(2,4)算法的E面雙站RCS計算結(jié)果如圖2所示。
從圖2可以看出,在一般的網(wǎng)格剖分條件下,基于高階子網(wǎng)格材質(zhì)平均技術(shù)的AM-FDTD(2,4)算法總體精度明顯得到了提高。
算例2 有耗介質(zhì)球的頻域近場分析。
電場z方向傳播,x方向極化。介質(zhì)球直徑為20cm,相對介電常數(shù)為30,電導(dǎo)率為0.3。采用一致的空間步長,參數(shù)設(shè)置為:Δδ=1cm,CFLδ=0.4,f=200MHz。本文采用離散時間傅里葉變換,計算x=Δδ/2,y=0處Ex場沿z軸的頻域波形。
圖2 無耗介質(zhì)球雙站RCS
采用FDTD(2,4)算法,空氣介質(zhì)分界面的處理分別采用階梯近似和子網(wǎng)格平均材質(zhì)模型(3×3子網(wǎng)格)。以細(xì)網(wǎng)格劃分(Δδ=0.5cm)的傳統(tǒng)FDTD(2,2)解為參考,有耗介質(zhì)球Ex場沿z方向的頻域近場值如圖3所示,從圖3可以看出,采用平均材質(zhì)模型后,其分界面處的電場值和參考解吻合良好。而階梯近似模型得到的曲線,有明顯的“上翹”現(xiàn)象。
圖3 有耗介質(zhì)球Ex場沿z方向的頻域近場值
本文給出了一種基于子網(wǎng)格材質(zhì)平均的高階FDTD(2,4)方法,在原網(wǎng)格基礎(chǔ)上劃分了更多的子網(wǎng)格,更加逼近介質(zhì)體曲面,提高了精度。該方法需在迭代前進(jìn)行求等效電參數(shù)的預(yù)處理,迭代仍在網(wǎng)格上進(jìn)行,除了需要存儲子網(wǎng)格電參數(shù)外,不涉及函數(shù)積分等復(fù)雜的運算過程,處理相對簡單,未增加計算量。相關(guān)算例表明該方法有效,簡化了建模過程,顯著提高了計算的精度。
[1]Yee K S.Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media[J].IEEE Transactions on Antennas and Propagation,1966,14(5):302-307.
[2]王 志,江谷傳.高階FDTD方法在三維散射問題中的應(yīng)用[J].合肥工業(yè)大學(xué)學(xué)報:自然科學(xué)版,2009,32(11):1776-1779.
[3]Okoniewski M,Okoniewska E,Stuchly M A.Three-dimensional subgridding algorithm for FDTD[J].IEEE Transactions on Antennas and Propagation,1997,45(3):422-429.
[3]Chevalier M W,Luebbers R J,Cable V P.FDTD local grid with materials transverse[J].IEEE Transactions Antennas and Propagat,1997,45(3):411-421.
[4]Sullivan D M.Electromagnetic simulation using the FDTD method[M].New York:Wiley-IEEE Press,2000:34-147.
[5]張世芳,褚慶昕.介質(zhì)體曲面的時域有限差分子網(wǎng)格建模[J].西安電子科技大學(xué)學(xué)報,2003,30(5):623-624.
[6]Fang J.Time domain finite difference computation for Maxwell’s equations[D].Berkeley:Univ California,1989.
[7]Yefet A,Petropoulos P G.A staggered fourth-order accurate explicit finite difference scheme for the time-domain Maxwell's equations[J].Journal of Computational Physics,2001,168(2):286-315.
[8]Zhao S,Wei G W.High-order FDTD methods via derivative matching for Maxwell’s equations with material interfaces[J].Journal of Computational Physics,2004,200(1):60-103.
[9]Georgakopoulos S V,Birtcher C R,Balanis C A,et al.HIRF penetration and PED coupling analysis for scaled fuselage models using a hybrid subgrid FDTD(2,2)/FDTD(2,4)method[J].IEEE Transactions on Electromagnetic Compatibility,2003,45(2):293-305.