三維新型非分裂場完全匹配層吸收邊界條件

楊利霞 梁 慶 于萍萍 王 剛

(1.江蘇大學 通信工程系，江蘇 鎮江 212013；2.東南大學 毫米波國家重點實驗室，江蘇 南京 210096)

1. 引 言

近年來，隨著計算機性能的不斷提高和數值理論的不斷發展，計算電磁學取得了很大的發展，目前已經形成了多種電磁學計算方法，例如矩量法(MoM)[1]、有限元法(FEM)[2]、時域有限差分法(FDTD)[3-4]等等。并且，隨著電磁波應用的廣泛和計算機技術的發展，各種方法的研究也更加深入。在FDTD計算中，為了在數值上模擬一個開放式系統，必須用吸收邊界條件來截斷計算區域。經過多年的研究與發展，吸收邊界的方法也在不斷改進，從1969年的Taylor首次提出用吸收邊界來吸收外向行波，當時采用的吸收邊界為簡單插值方法，到1981年Mur提出了在計算區域截斷邊界處的一階和二階吸收邊界條件，我們把這種方法稱為Mur邊界[5]，這是一種比較有效的吸收邊界，1994年由Berenger提出將麥克斯韋方程擴展為場分裂形式的完全匹配層(PML)[6-7]，并由Sacks和Gedney等人又提出了各向異性介質的PML，即UPML[8-10]。隨著研究的深入更多的方法被學者們提出，如卷積完全匹配層吸收邊界(CPML)[11],駐行波吸收邊界條件

(STWBC)[12]等。

2003年由Cummer提出了近似完全匹配層(Nearly Perfectly Matched Layer，NPML)[13]吸收邊界，后來在2004年經Hu和Cummer證實NPML是一種完全匹配層吸收邊界[14]，但是他們只給出了簡單的二維算例分析，本文在此基礎上仿真分析了三維情況下的NPML吸收邊界。這種方法在PML吸收邊界的基礎上通過拉伸坐標系進行坐標變換來實現[15]，在截斷常規線性電磁介質材料時與PML具有相同的偏微分方程，該方法避免了PML的場分裂形式，編程復雜度有了較大程度的降低，因而在計算處理過程中具有非常明顯的優勢，尤其是編程較為簡單，因而更容易實現。

2. 理論分析及公式推導

2.1 Maxwell支配方程

麥克斯韋旋度方程為

(1)

(2)

在直角坐標系中麥克斯韋旋度方程寫為

(3)

(4)

(5)

(6)

(7)

(8)

2.2 NPML吸收邊界FDTD離散式的推導

在直角坐標系FDTD計算中，NPML吸收邊界可以分為棱邊區、面區和角頂區三種情況。首先以x,y方向交叉的棱邊為例分析。根據文獻[13]依據如下坐標變換規則

(9)

分別用?x′和?y′代替?x和?y，依據式(9)的處理方式，在時域，棱邊區的電磁場的Maxwell支配方程可以寫為

(10)

(11)

(12)

(13)

(14)

(15)

以式(10)為例，令Sx=1+σx/ jωε0，Sy=1+σy/ jωε0，則有

故該式可化為

(16)

通過文獻[14]，可知

因此，上式又可以寫為

(17)

(18)

式(10)可化為

(19)

同理式(11)～(15)可化為

(20)

(21)

(22)

(23)

(24)

下面以式(19)為例進行中心差分離散，有

(25)

同理，對式(20)～(24)離散可得到

(26)

(27)

(28)

(29)

(30)

式中：

m代表觀察點(x,y,z)處的一組整數或半整數。而引入變量的離散方式也可以按照中心差分方式進行離散如下

(31)

(32)

綜上，可以得到三維情況下NPML吸收邊界的FDTD迭代計算的推進過程如下

NPML其他區域的FDTD公式可參照以上形式進行推導，角頂區的坐標變換規則是將?x，?y，?z同時進行變換?x′，?y′，?z′，棱邊區則變換其中兩個方向，而面區則只變換一個方向，根據此種變換方式可以將NPML吸收邊界的遞推公式分別推導出，這里就不再贅述。

2.3 編程思想

通過以上分析，可以看出，在編程實現NPML吸收邊界的時候，可以用其角頂區的遞推公式作為一般公式，通過對σM(M=x,y,z)的討論來區分不同的區域，從而避免了對吸收邊界內不同方向、不同棱邊以及角點不同公式的區分，降低了程序上的復雜度，方便編程實現，如圖1所示三維NPML吸收邊界參數分布。而對于整體的FDTD計算區域來講，同樣可根據NPML的特點，在模擬計算區域和吸收邊界內利用NPML角點處的方程而統一進行編程，在模擬計算區，即σM(M=x,y,z)為0時，設置F=F′(F′代表E′和H′,F代表E和H)，這樣即使吸收邊界內的方程與模擬計算區域內的迭代式不同，也可以將兩者的計算用統一的公式實現編程。從而回避了NPML區與模擬區域邊界數據處理問題，可以簡化編程，不易出錯。

圖1 三維NPML吸收邊界參數分布

3. 算例驗證及分析

3.1 與傳統吸收邊界的對比

考慮三維空間中的電偶極子輻射。設垂直電偶極子位于計算區域中心Ez(0, 0, 1/2)處，FDTD計算區域為24×24×24個元胞，四周被吸收層包圍。元胞尺寸為5 cm×5 cm×5 cm，即δ=5 cm。按照δ=2cΔt，時間間隔為Δt=83.333 ps。考察距離輻射源10δ處觀察點Q點的電場Ez(10, 10, 1/2)。吸收層內側距離觀察點Q兩個元胞。計算中輻射源采用高斯脈沖

(33)

計算結果如圖2所示。圖中給出了NPML、UPML、 傳統PML結果與瞬態偶極子的解析解進行對比，由圖可見曲線符合很好，驗證了NPML算法具有良好的吸收效果。

圖2 Q點處電偶極子的輻射場

3.2 計算金屬球的RCS

金屬球半徑為1 m，δ=0.05 m, Δt=δ/2c，目標區域為40×40×40δ3，入射波采用高斯脈沖波。圖3為球的后向遠區散射電場隨時間的變化， 金屬球的單站RCS隨頻率的變化，如圖4所示。作為比較，圖中還給出了Mie級數解，以及UPML方法得帶的RCS圖。在0～350 MHz范圍內三者符合得很好。

3.3 計算金屬平板的RCS

計算模型為方形平板,邊長為29 cm,厚度為1 cm。FDTD網格δ=1 cm,Δt=δ/2c,入射波為高斯脈沖,入射波沿z方向，電場E分量沿+x方向。圖5給出其后向遠區散射電場隨時間的變化，圖6為將計算結果經傅立葉變換后得到的單站RCS隨頻率的變化。為了對照，圖中還給出了矩量法的結果,以及UPML方法的結果，可見三者符合得很好。

圖3 時域響應

圖4 頻域響應

圖5 時域響應

圖6 頻域響應

4. 結果分析

通過以上兩個三維算例，可以看出NPML吸收邊界能夠很好地達到吸收要求，程序編制實現比較簡單。通過采用NPML吸收邊界，UPML吸收邊界及相關解析數值解比較，可以看出在低頻部分，采用NPML，UPML的兩個算例都能與解析解很好地吻合，但在高頻部分，NPML方法得到的結果與MOM方法的結果吻合得比UPML方法的結果好。該方法在截斷復雜的色散介質中具有較大的優勢，是我們下一步的研究方向。

[1] 李世智. 電磁輻射與散射問題的矩量法[M]. 北京: 電子工業出版社, 1985.

[2] 吳 霞, 周樂柱. 時域有限元法在計算電磁問題上的應用及發展[J]. 電波科學學報, 2008, 23(6): 1209-1216.

WU Xia, ZHOU Lezhu. Application and development of time-domain finite element method on EM analysis[J]. Chinese Journal of Radio Science, 2008, 23(6): 1209-1216. (in Chinese)

[3] KUNZ K S, LUEBBERS R J. The Finite Difference Time Domain Method for Electromagnetics[M]. Boca Raton, FL: CRC Press, 1993.

[4] 葛德彪, 閆玉波. 電磁波時域有限差分方法[M]. 西安: 西安電子科技大學出版社, 2005.

[5] MUR G. Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations[J]. IEEE Trans. Electromagn. Compat., 1981, EMC-23(4): 377-382.

[6] B RENGER J P. A perfectly matched layer for the absorption of electromagnetic waves[J]. J. Comput. Phys., 1994, 114(2): 185-200.

[7]B RENGER J P. Perfectly matched layer for the FDTD solutions of wavestructure interaction problems[J]. IEEE Trans. Antennas Propagat., 1996, AP-44(1): 110-117.

[8] SACKS Z S, KINGSLAND D M, LEE R, et al. Aperfectly matched anisotropic absorber for use as an absorbing boundary condition[J]. IEEE Trans. Antennas Propagat., 1995, AP-43(12): 1460-1463.

[9] KATZ D C, THIELE E T, and TAFLOVE A.Validation and extension to three dimensions of the Bérenger absorbing boundary condition for FDTD meshes[J]. IEEE Microwave and Guided Wave Letters, 1994, 4(8): 268-270.

[10] GEDNEY S D. An anisotropic PML absorbing media for the FDTD simulation of fields in lossy and dispersive media[J]. Electromagn., 1996, 16 (4): 399-415.

[11] RODEN J A and GEDNEY S D. Convolutional PML (CPML): An efficient FDTD implementation of the CFS-PML for arbitrary media[J]. Microw. Opt. Technol. Lett., 2000, 27(5): 334-339.

[12] 譚懷英, 梁甸農, 劉克成. FDTD的一種新型 吸收邊界-STWBC[J].電波科學學報, 2001, 16(1): 412-413.

TAN Huaiying, LIANG Diannong, LIU Kecheng. Standing-traveling wave boundary condition (STWBC) for finite-difference time-domain mesh truncation[J]. Chinese Journal of Radio Science, 2001, 16(1): 412-413. (in Chinese)

[13] CUMMER S A. A simple, nearly perfectly matched layer for general electromagnetic media[J]. IEEE Micr-owave Wireless Components Lett., 2003, 13(3): 128-130.

[14] HU W Y and CUMMER S A. The nearly perfectly matched layer is a perfectly matched layer[J]. IEEE Microwave Wireless Components Lett., 2004, 3(1): 137-140.

[15] CHEW W C and WEEDON W H. A 3D perfectly matched medium from modified Maxwell’s equations with stretched coordinates[J]. IEEE Microwave Opt. Technol. Lett., 1994, 7(13): 599-604.