楊利霞 梁 慶 于萍萍 王 剛
(1.江蘇大學(xué) 通信工程系,江蘇 鎮(zhèn)江 212013;2.東南大學(xué) 毫米波國家重點(diǎn)實(shí)驗(yàn)室,江蘇 南京 210096)
近年來,隨著計(jì)算機(jī)性能的不斷提高和數(shù)值理論的不斷發(fā)展,計(jì)算電磁學(xué)取得了很大的發(fā)展,目前已經(jīng)形成了多種電磁學(xué)計(jì)算方法,例如矩量法(MoM)[1]、有限元法(FEM)[2]、時(shí)域有限差分法(FDTD)[3-4]等等。并且,隨著電磁波應(yīng)用的廣泛和計(jì)算機(jī)技術(shù)的發(fā)展,各種方法的研究也更加深入。在FDTD計(jì)算中,為了在數(shù)值上模擬一個(gè)開放式系統(tǒng),必須用吸收邊界條件來截?cái)嘤?jì)算區(qū)域。經(jīng)過多年的研究與發(fā)展,吸收邊界的方法也在不斷改進(jìn),從1969年的Taylor首次提出用吸收邊界來吸收外向行波,當(dāng)時(shí)采用的吸收邊界為簡單插值方法,到1981年Mur提出了在計(jì)算區(qū)域截?cái)噙吔缣幍囊浑A和二階吸收邊界條件,我們把這種方法稱為Mur邊界[5],這是一種比較有效的吸收邊界,1994年由Berenger提出將麥克斯韋方程擴(kuò)展為場分裂形式的完全匹配層(PML)[6-7],并由Sacks和Gedney等人又提出了各向異性介質(zhì)的PML,即UPML[8-10]。隨著研究的深入更多的方法被學(xué)者們提出,如卷積完全匹配層吸收邊界(CPML)[11],駐行波吸收邊界條件
(STWBC)[12]等。
2003年由Cummer提出了近似完全匹配層(Nearly Perfectly Matched Layer,NPML)[13]吸收邊界,后來在2004年經(jīng)Hu和Cummer證實(shí)NPML是一種完全匹配層吸收邊界[14],但是他們只給出了簡單的二維算例分析,本文在此基礎(chǔ)上仿真分析了三維情況下的NPML吸收邊界。這種方法在PML吸收邊界的基礎(chǔ)上通過拉伸坐標(biāo)系進(jìn)行坐標(biāo)變換來實(shí)現(xiàn)[15],在截?cái)喑R?guī)線性電磁介質(zhì)材料時(shí)與PML具有相同的偏微分方程,該方法避免了PML的場分裂形式,編程復(fù)雜度有了較大程度的降低,因而在計(jì)算處理過程中具有非常明顯的優(yōu)勢,尤其是編程較為簡單,因而更容易實(shí)現(xiàn)。
麥克斯韋旋度方程為
(1)
(2)
在直角坐標(biāo)系中麥克斯韋旋度方程寫為
(3)
(4)
(5)
(6)
(7)
(8)
在直角坐標(biāo)系FDTD計(jì)算中,NPML吸收邊界可以分為棱邊區(qū)、面區(qū)和角頂區(qū)三種情況。首先以x,y方向交叉的棱邊為例分析。根據(jù)文獻(xiàn)[13]依據(jù)如下坐標(biāo)變換規(guī)則
(9)
分別用?x′和?y′代替?x和?y,依據(jù)式(9)的處理方式,在時(shí)域,棱邊區(qū)的電磁場的Maxwell支配方程可以寫為
(10)
(11)
(12)
(13)
(14)
(15)
以式(10)為例,令Sx=1+σx/ jωε0,Sy=1+σy/ jωε0,則有
故該式可化為
(16)
通過文獻(xiàn)[14],可知
因此,上式又可以寫為
(17)
(18)
式(10)可化為
(19)
同理式(11)~(15)可化為
(20)
(21)
(22)
(23)
(24)
下面以式(19)為例進(jìn)行中心差分離散,有
(25)
同理,對式(20)~(24)離散可得到
(26)
(27)
(28)
(29)
(30)
式中:
m代表觀察點(diǎn)(x,y,z)處的一組整數(shù)或半整數(shù)。而引入變量的離散方式也可以按照中心差分方式進(jìn)行離散如下
(31)
(32)
綜上,可以得到三維情況下NPML吸收邊界的FDTD迭代計(jì)算的推進(jìn)過程如下
NPML其他區(qū)域的FDTD公式可參照以上形式進(jìn)行推導(dǎo),角頂區(qū)的坐標(biāo)變換規(guī)則是將?x,?y,?z同時(shí)進(jìn)行變換?x′,?y′,?z′,棱邊區(qū)則變換其中兩個(gè)方向,而面區(qū)則只變換一個(gè)方向,根據(jù)此種變換方式可以將NPML吸收邊界的遞推公式分別推導(dǎo)出,這里就不再贅述。
通過以上分析,可以看出,在編程實(shí)現(xiàn)NPML吸收邊界的時(shí)候,可以用其角頂區(qū)的遞推公式作為一般公式,通過對σM(M=x,y,z)的討論來區(qū)分不同的區(qū)域,從而避免了對吸收邊界內(nèi)不同方向、不同棱邊以及角點(diǎn)不同公式的區(qū)分,降低了程序上的復(fù)雜度,方便編程實(shí)現(xiàn),如圖1所示三維NPML吸收邊界參數(shù)分布。而對于整體的FDTD計(jì)算區(qū)域來講,同樣可根據(jù)NPML的特點(diǎn),在模擬計(jì)算區(qū)域和吸收邊界內(nèi)利用NPML角點(diǎn)處的方程而統(tǒng)一進(jìn)行編程,在模擬計(jì)算區(qū),即σM(M=x,y,z)為0時(shí),設(shè)置F=F′(F′代表E′和H′,F代表E和H),這樣即使吸收邊界內(nèi)的方程與模擬計(jì)算區(qū)域內(nèi)的迭代式不同,也可以將兩者的計(jì)算用統(tǒng)一的公式實(shí)現(xiàn)編程。從而回避了NPML區(qū)與模擬區(qū)域邊界數(shù)據(jù)處理問題,可以簡化編程,不易出錯(cuò)。
圖1 三維NPML吸收邊界參數(shù)分布
考慮三維空間中的電偶極子輻射。設(shè)垂直電偶極子位于計(jì)算區(qū)域中心Ez(0, 0, 1/2)處,F(xiàn)DTD計(jì)算區(qū)域?yàn)?4×24×24個(gè)元胞,四周被吸收層包圍。元胞尺寸為5 cm×5 cm×5 cm,即δ=5 cm。按照δ=2cΔt,時(shí)間間隔為Δt=83.333 ps??疾炀嚯x輻射源10δ處觀察點(diǎn)Q點(diǎn)的電場Ez(10, 10, 1/2)。吸收層內(nèi)側(cè)距離觀察點(diǎn)Q兩個(gè)元胞。計(jì)算中輻射源采用高斯脈沖
(33)
計(jì)算結(jié)果如圖2所示。圖中給出了NPML、UPML、 傳統(tǒng)PML結(jié)果與瞬態(tài)偶極子的解析解進(jìn)行對比,由圖可見曲線符合很好,驗(yàn)證了NPML算法具有良好的吸收效果。
圖2 Q點(diǎn)處電偶極子的輻射場
金屬球半徑為1 m,δ=0.05 m, Δt=δ/2c,目標(biāo)區(qū)域?yàn)?0×40×40δ3,入射波采用高斯脈沖波。圖3為球的后向遠(yuǎn)區(qū)散射電場隨時(shí)間的變化, 金屬球的單站RCS隨頻率的變化,如圖4所示。作為比較,圖中還給出了Mie級數(shù)解,以及UPML方法得帶的RCS圖。在0~350 MHz范圍內(nèi)三者符合得很好。
計(jì)算模型為方形平板,邊長為29 cm,厚度為1 cm。FDTD網(wǎng)格δ=1 cm,Δt=δ/2c,入射波為高斯脈沖,入射波沿z方向,電場E分量沿+x方向。圖5給出其后向遠(yuǎn)區(qū)散射電場隨時(shí)間的變化,圖6為將計(jì)算結(jié)果經(jīng)傅立葉變換后得到的單站RCS隨頻率的變化。為了對照,圖中還給出了矩量法的結(jié)果,以及UPML方法的結(jié)果,可見三者符合得很好。
圖3 時(shí)域響應(yīng)
圖4 頻域響應(yīng)
圖5 時(shí)域響應(yīng)
圖6 頻域響應(yīng)
通過以上兩個(gè)三維算例,可以看出NPML吸收邊界能夠很好地達(dá)到吸收要求,程序編制實(shí)現(xiàn)比較簡單。通過采用NPML吸收邊界,UPML吸收邊界及相關(guān)解析數(shù)值解比較,可以看出在低頻部分,采用NPML,UPML的兩個(gè)算例都能與解析解很好地吻合,但在高頻部分,NPML方法得到的結(jié)果與MOM方法的結(jié)果吻合得比UPML方法的結(jié)果好。該方法在截?cái)鄰?fù)雜的色散介質(zhì)中具有較大的優(yōu)勢,是我們下一步的研究方向。
[1] 李世智. 電磁輻射與散射問題的矩量法[M]. 北京: 電子工業(yè)出版社, 1985.
[2] 吳 霞, 周樂柱. 時(shí)域有限元法在計(jì)算電磁問題上的應(yīng)用及發(fā)展[J]. 電波科學(xué)學(xué)報(bào), 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] 葛德彪, 閆玉波. 電磁波時(shí)域有限差分方法[M]. 西安: 西安電子科技大學(xué)出版社, 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] 譚懷英, 梁甸農(nóng), 劉克成. FDTD的一種新型 吸收邊界-STWBC[J].電波科學(xué)學(xué)報(bào), 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.