陳 衛(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é)果，就需要進行細網(wǎng)格劃分，從而造成很大的計算負擔。針對第1個問題，亞網(wǎng)格技術(shù)［2－3］、材質(zhì)平均技術(shù)［4－5］等已用于傳統(tǒng)的FDTD方法；對于第2個問題，文獻［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］，其收斂率可達4階，但后時不穩(wěn)定現(xiàn)象和繁瑣的操作，使其無法推廣?；旌蟻喚W(wǎng)格（Hybrid－Subgridding）［9］技術(shù)的思想是在目標邊界附近采用FDTD（2，2）方法，網(wǎng)格是細網(wǎng)格；遠離邊界處采用高階FDTD（2，4）方法，網(wǎng)格是粗網(wǎng)格。由于差分精度、網(wǎng)格尺寸及穩(wěn)定度等因素的不同，其需要采用時間和空間的雙插值方案來減少粗細網(wǎng)格之間的非物理反射，這必然帶來計算精度的降低。

在前人研究的基礎(chǔ)上，本文提出了一種新的介質(zhì)體曲面問題的FDTD（2，4）建模方法，通過對原網(wǎng)格內(nèi)外環(huán)路所圍區(qū)域進行二次剖分，分別得到內(nèi)外環(huán)路的平均電磁參數(shù)，再通過加權(quán)平均得到等效電參數(shù)，避免了積分運算，簡化了建模過程。FDTD（2，4）運算仍在原網(wǎng)格上進行，所以沒有增加計算量。相關(guān)數(shù)值結(jié)果表明，在相同計算量的條件下，本文方法具有更高的精度。

1 高階子網(wǎng)格材質(zhì)平均FDTD（2，4）方法

高階差分本質(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)完成。

2 數(shù)值結(jié)果

算例1 無耗介質(zhì)球的頻域遠場分析。

球半徑為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)格）。以細網(wǎng)格劃分（Δδ＝0.5cm）的傳統(tǒng)FDTD（2，2）解為參考，有耗介質(zhì)球Ex場沿z方向的頻域近場值如圖3所示，從圖3可以看出，采用平均材質(zhì)模型后，其分界面處的電場值和參考解吻合良好。而階梯近似模型得到的曲線，有明顯的“上翹”現(xiàn)象。

圖3 有耗介質(zhì)球Ex場沿z方向的頻域近場值

3 結(jié)束語

本文給出了一種基于子網(wǎng)格材質(zhì)平均的高階FDTD（2，4）方法，在原網(wǎng)格基礎(chǔ)上劃分了更多的子網(wǎng)格，更加逼近介質(zhì)體曲面，提高了精度。該方法需在迭代前進行求等效電參數(shù)的預(yù)處理，迭代仍在網(wǎng)格上進行，除了需要存儲子網(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.