CPML在3維HIE-FDTD 算法中的實現

李 建，劉宗信，王發年，李 斌，陳桂杰

（解放軍95972部隊 甘肅 酒泉 735018）

時域有限差分（FDTD）方法具有方法簡單、適用性強等優點，已被廣泛應用。然而，由于常規FDTD算法的時間步長受穩定性條件的限制，而時間步長的選取是由最小空間步長決定的[1]，這樣小的步長會使迭代步數和模擬時間增加。1999年T.Namiki提出了一種新的交替方向隱格式時域有限差分算法（ADI-FDTD），該方法具有無條件穩定的特性，因而時間步長不再受穩定性條件的限制[2-3]。然而ADI-FDTD方法在一個時間步上需要兩次迭代，增加了每步的計算時間和計算機的存儲量，而且隨著時間步的增大，數值色散變大，降低了數值模擬的精度。

為了克服以上算法的缺點，B.K.Huang等人提出了弱無條件穩定算法[4]。該方法是有條件穩定的，但其穩定性條件比常規FDTD要寬松，而且該方法比ADI方法有更高的計算精度，計算時間上也比ADI方法有優勢。CFS-PML吸收邊界以其簡單、直觀、高吸收性能（尤其在吸收凋落波方面）著稱，一經提出就被廣泛使用，但是該邊界條件不便于實現。本文擬將CFS-PML技術引入到混合顯隱式弱無條件穩定時域有限差分算法（HIE-FDTD）中，研究其具體差分格式、吸波性能及如何合理設置其復頻率參數。

1 公式推導

借助輔助差分方程，可以將三維情形下，PML局域內的麥克斯韋方程表示如下：

按照HIE方法[4]對式(1)～(6)進行差分離散可以得到：

其中，△t為時間步長，△x，△y，△z是 x軸、y軸和 z軸方向的網格大小。使用標準Yee格式對公式(7)～(18)中的輔助變量進行離散可以得到一系列差分方程，這里僅給出(7)的差分格式，其余各項按同樣的方法均可得到。

從公式(19)～(24)和(25)可以看出，對于電場場量 Ey和 Hy可以直接更新得到。將未知磁場量(24)式代入(19)式中，連立(15)、(16)、(21)、(22) 可得一關于電場分量 Ex在 n+1 時間步的一個隱式迭代表達，如式(26)所示：

式(26)可以用Thomas方法解決，求解后可以通過顯式直接更新。同理可以求得和。從上述公式可以看出，三維HIEFDTD方法在整個循環迭代過程中，需要解算2個三對角矩陣，進行4次顯式迭代和12個輔助變量的運算。然而三維ADI-FDTD方法在整個循環迭代過程中，需要解算3個三對角矩陣，進行3次顯式迭代和12個輔助變量的運算。顯然，HIE-FDTD比ADI-FDTD有著更好的計算效率。

2 數值結果

為了驗證CFS-PML在三維HIE-FDTD算法中的精確性和有效性，采用圖1所示的幾何模型，采用均勻空間步長劃分△=△x=△y=△z=0.5 mm，網格為 66×66×66，每個方向上的CPML層均設置為8層用來截斷計算局域[5]。觀察點設置在A(57,33,33)點和B(57,57,33)。采用式(6)所示形式的高斯脈沖作為激烈源放在整個結構的正中間，其中：f0=1 GHz ，t0=1.0×10-10,τ=1.0×10-10。

CPML層中的參數由式(7)確定，本文取m=4，α取固定值。

圖1 數值示例的幾何模型Fig.1 Geometry of the numerical example

為了便于比較，設立一266×266×266的參考空間，在觀察點，關于時間函數的相對反射誤差及各參數的含義同式(31)。

圖2顯示了在 HIE-FDTD中使用Mur吸收邊界、在HIE-FDTD中使用CPML吸收邊界以及在傳統FDTD中使用CPML吸收邊界時相對反射誤差隨時間的變化曲線。從圖中顯示的計算結果可以得出本文所提出的CPML吸收邊界條件比一階Mur吸收邊界具有更好的吸波性能，在整個變化過程中，HIE-CPML比HIE-Mur的反射誤差平均要低30 dB。

下面，我們來觀察PML函數的基本構成變量κmax，σmax以及α與其最大反射誤差的關系。圖3、圖4分別是α=0.05、α=0.003時κmax和σmax與反射誤差的等高圖。從圖可以看出，當α=0.05、α=0.003 時，最大誤差均為-72 dB，但當取 α=0.05，可以在一個較大范圍內合理選取κmax和σmax來實現最佳誤差，從而使得在選值時易于預測。顯而易見，選取κmax=14、σmax/σopt=1.0就可以實現低至-72 dB的最大相對誤差。

圖2 不同方法的相對反射誤差Fig.2 Relative reflection error for different methods

圖3 α=0.05、m=4時觀察點觀測到的關于 κmax和 σmax/σopt函數的最大反射誤差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時觀察點觀測到的關于κmax和σmax/σopt函數的最大反射誤差Fig.4 Maximum reflection error at observation point as a function of κmaxand σmax/σopt(α=0.003 and m=4)

3 結論

本文將坐標伸縮完全匹配層CPML引入到3維弱無條件穩定算法HIE-FDTD中研究其吸波性能。詳細推導了CPML在3維HIE-FDTD算法中的差分公式，建立了計算模型，并將其與幾種常用的吸收邊界條件的吸波性能進行了綜合比較。數值結果表明：當將本文所提方法的CPML層數設置為8時，其最大反射誤差為-72 dB，遠低于傳統FDTD方法的反射誤差。另外，當匹配層參數設置為α=0.05，可以在一個較大范圍內合理選取κmax和σmax來實現最佳誤差，從而使得在選值時易于預測反射情況。

[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.