李 建,劉宗信,王發(fā)年,李 斌,陳桂杰
(解放軍95972部隊(duì) 甘肅 酒泉 735018)
時(shí)域有限差分(FDTD)方法具有方法簡單、適用性強(qiáng)等優(yōu)點(diǎn),已被廣泛應(yīng)用。然而,由于常規(guī)FDTD算法的時(shí)間步長受穩(wěn)定性條件的限制,而時(shí)間步長的選取是由最小空間步長決定的[1],這樣小的步長會(huì)使迭代步數(shù)和模擬時(shí)間增加。1999年T.Namiki提出了一種新的交替方向隱格式時(shí)域有限差分算法(ADI-FDTD),該方法具有無條件穩(wěn)定的特性,因而時(shí)間步長不再受穩(wěn)定性條件的限制[2-3]。然而ADI-FDTD方法在一個(gè)時(shí)間步上需要兩次迭代,增加了每步的計(jì)算時(shí)間和計(jì)算機(jī)的存儲(chǔ)量,而且隨著時(shí)間步的增大,數(shù)值色散變大,降低了數(shù)值模擬的精度。
為了克服以上算法的缺點(diǎn),B.K.Huang等人提出了弱無條件穩(wěn)定算法[4]。該方法是有條件穩(wěn)定的,但其穩(wěn)定性條件比常規(guī)FDTD要寬松,而且該方法比ADI方法有更高的計(jì)算精度,計(jì)算時(shí)間上也比ADI方法有優(yōu)勢(shì)。CFS-PML吸收邊界以其簡單、直觀、高吸收性能(尤其在吸收凋落波方面)著稱,一經(jīng)提出就被廣泛使用,但是該邊界條件不便于實(shí)現(xiàn)。本文擬將CFS-PML技術(shù)引入到混合顯隱式弱無條件穩(wěn)定時(shí)域有限差分算法(HIE-FDTD)中,研究其具體差分格式、吸波性能及如何合理設(shè)置其復(fù)頻率參數(shù)。
借助輔助差分方程,可以將三維情形下,PML局域內(nèi)的麥克斯韋方程表示如下:
按照HIE方法[4]對(duì)式(1)~(6)進(jìn)行差分離散可以得到:
其中,△t為時(shí)間步長,△x,△y,△z是 x軸、y軸和 z軸方向的網(wǎng)格大小。使用標(biāo)準(zhǔn)Yee格式對(duì)公式(7)~(18)中的輔助變量進(jìn)行離散可以得到一系列差分方程,這里僅給出(7)的差分格式,其余各項(xiàng)按同樣的方法均可得到。
從公式(19)~(24)和(25)可以看出,對(duì)于電場場量 Ey和 Hy可以直接更新得到。將未知磁場量(24)式代入(19)式中,連立(15)、(16)、(21)、(22) 可得一關(guān)于電場分量 Ex在 n+1 時(shí)間步的一個(gè)隱式迭代表達(dá),如式(26)所示:
式(26)可以用Thomas方法解決,求解后可以通過顯式直接更新。同理可以求得和。從上述公式可以看出,三維HIEFDTD方法在整個(gè)循環(huán)迭代過程中,需要解算2個(gè)三對(duì)角矩陣,進(jìn)行4次顯式迭代和12個(gè)輔助變量的運(yùn)算。然而三維ADI-FDTD方法在整個(gè)循環(huán)迭代過程中,需要解算3個(gè)三對(duì)角矩陣,進(jìn)行3次顯式迭代和12個(gè)輔助變量的運(yùn)算。顯然,HIE-FDTD比ADI-FDTD有著更好的計(jì)算效率。
為了驗(yàn)證CFS-PML在三維HIE-FDTD算法中的精確性和有效性,采用圖1所示的幾何模型,采用均勻空間步長劃分△=△x=△y=△z=0.5 mm,網(wǎng)格為 66×66×66,每個(gè)方向上的CPML層均設(shè)置為8層用來截?cái)嘤?jì)算局域[5]。觀察點(diǎn)設(shè)置在A(57,33,33)點(diǎn)和B(57,57,33)。采用式(6)所示形式的高斯脈沖作為激烈源放在整個(gè)結(jié)構(gòu)的正中間,其中:f0=1 GHz ,t0=1.0×10-10,τ=1.0×10-10。
CPML層中的參數(shù)由式(7)確定,本文取m=4,α取固定值。
圖1 數(shù)值示例的幾何模型Fig.1 Geometry of the numerical example
為了便于比較,設(shè)立一266×266×266的參考空間,在觀察點(diǎn),關(guān)于時(shí)間函數(shù)的相對(duì)反射誤差及各參數(shù)的含義同式(31)。
圖2顯示了在 HIE-FDTD中使用Mur吸收邊界、在HIE-FDTD中使用CPML吸收邊界以及在傳統(tǒng)FDTD中使用CPML吸收邊界時(shí)相對(duì)反射誤差隨時(shí)間的變化曲線。從圖中顯示的計(jì)算結(jié)果可以得出本文所提出的CPML吸收邊界條件比一階Mur吸收邊界具有更好的吸波性能,在整個(gè)變化過程中,HIE-CPML比HIE-Mur的反射誤差平均要低30 dB。
下面,我們來觀察PML函數(shù)的基本構(gòu)成變量κmax,σmax以及α與其最大反射誤差的關(guān)系。圖3、圖4分別是α=0.05、α=0.003時(shí)κmax和σmax與反射誤差的等高圖。從圖可以看出,當(dāng)α=0.05、α=0.003 時(shí),最大誤差均為-72 dB,但當(dāng)取 α=0.05,可以在一個(gè)較大范圍內(nèi)合理選取κmax和σmax來實(shí)現(xiàn)最佳誤差,從而使得在選值時(shí)易于預(yù)測(cè)。顯而易見,選取κmax=14、σmax/σopt=1.0就可以實(shí)現(xiàn)低至-72 dB的最大相對(duì)誤差。
圖2 不同方法的相對(duì)反射誤差Fig.2 Relative reflection error for different methods
圖3 α=0.05、m=4時(shí)觀察點(diǎn)觀測(cè)到的關(guān)于 κmax和 σmax/σopt函數(shù)的最大反射誤差Fig.3 Maximum reflection error at observation point as a function of κmaxand σmax/σopt(α=0.05 and m=4)
圖4 α=0.003、m=4時(shí)觀察點(diǎn)觀測(cè)到的關(guān)于κmax和σmax/σopt函數(shù)的最大反射誤差Fig.4 Maximum reflection error at observation point as a function of κmaxand σmax/σopt(α=0.003 and m=4)
本文將坐標(biāo)伸縮完全匹配層CPML引入到3維弱無條件穩(wěn)定算法HIE-FDTD中研究其吸波性能。詳細(xì)推導(dǎo)了CPML在3維HIE-FDTD算法中的差分公式,建立了計(jì)算模型,并將其與幾種常用的吸收邊界條件的吸波性能進(jìn)行了綜合比較。數(shù)值結(jié)果表明:當(dāng)將本文所提方法的CPML層數(shù)設(shè)置為8時(shí),其最大反射誤差為-72 dB,遠(yuǎn)低于傳統(tǒng)FDTD方法的反射誤差。另外,當(dāng)匹配層參數(shù)設(shè)置為α=0.05,可以在一個(gè)較大范圍內(nèi)合理選取κmax和σmax來實(shí)現(xiàn)最佳誤差,從而使得在選值時(shí)易于預(yù)測(cè)反射情況。
[1]Taflove A,Hagness S C.Computational Electrodynamics:the Finite-DifferenceTime-domain Method[M].Boston,MA:Artech House,2000.
[2]Namiki T,Ito K.Numerical simulation using ADI-FDTD method to estimate shielding effectiveness of thin conductive enclosures[J].IEEE Trans.Microwave Theory Tech,2001,49(6):1060-1066.
[3]Yuan C H,Chen Z Z.Toward accurate time-domain simulation of highly conductive materials[J].IEEE Trans.Microwave Theory Tech.-s Digest,2002,(WE4E-4):1135-1138.
[4]Binke Huang,Gang Wang,Yansheng Jiang.A Hybrid implicitexplicit FDTD scheme with weakly conditional stability[J].Microwave and Optical Tech.Lett,2003,39(2):97-101.
[5]LIU Zong-xin,CHEN Yi-wang,SUN Xue-gang,et al.Implementation of CFS-PML for HIE-FDTD Method[J].IEEE Antennas and Wireless Propagation Letters,2012(11):381-384.
[6]Mao Y F,Chen B,Chen H L,et al.Unconditionally stable SFDTD algorithm for solving oblique incident wave on periodic structures[J].IEEE Microw.Wireless Comp.Lett.,2009,19(5):257-259.