CPML在3維HIE-FDTD 算法中的實(shí)現(xiàn)

李 建，劉宗信，王發(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ù)。

1 公式推導(dǎo)

借助輔助差分方程，可以將三維情形下，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ì)于電場(chǎng)場(chǎng)量 Ey和 Hy可以直接更新得到。將未知磁場(chǎng)量(24)式代入(19)式中，連立(15)、(16)、(21)、(22) 可得一關(guān)于電場(chǎng)分量 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ì)算效率。

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

為了驗(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)

3 結(jié)論

本文將坐標(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.