二維DGTD方法中UPML吸收邊界的實現

李林茜，魏 兵，楊 謙，葛德彪

(1. 西安電子科技大學 物理與光電工程學院，陜西 西安 710071; 2. 西安電子科技大學 信息感知協同創新中心，陜西 西安 710071)

?

二維DGTD方法中UPML吸收邊界的實現

李林茜1,2，魏 兵1,2，楊 謙1,2，葛德彪1,2

(1. 西安電子科技大學 物理與光電工程學院，陜西 西安 710071; 2. 西安電子科技大學 信息感知協同創新中心，陜西 西安 710071)

針對二維情形單軸各向異性完全匹配層吸收邊界條件，研究了橫磁波情形時域離散伽略金算法單軸各向異性完全匹配層吸收邊界的理論基礎和實現方案．借鑒時域有限差分方法——時域離散伽略金算法中吸收邊界阻抗匹配、各向異性介質參數設置和匹配矩陣等思想，結合時域離散伽略金算法空間離散網格的非結構特性和離散單元之間場量傳遞的特點，給出了在單軸各向異性完全匹配層中電磁場量時域迭代公式，進一步離散成為離散時域迭代計算式．由于時域離散伽略金算法網格的非結構特性，一階SM吸收邊界條件一般對入射電磁波的反射率在 -24 dB 左右．仿真算例說明，給出的時域離散伽略金算法單軸各向異性完全匹配層吸收邊界對電磁波的雙重衰減達 130 dB，表明單軸各向異性完全匹配層具有良好的吸收效果．關鍵詞：時域離散伽略金方法;單軸各向異性完全匹配層;橫磁波

截斷邊界條件是許多電磁場數值方法精度保證的關鍵．在過去幾十年，學著們提出許多種吸收邊界條件，并成功應用于時域有限差分(Finite Difference Time Domain， FDTD)等方法中[1-2]．文獻[3]中提出了完全匹配層(Perfectly Matched Layer， PML)吸收邊界，該邊界是在計算域最外面增加一層非物理吸收層，該吸收層與相鄰區域阻抗匹配，進入PML層的電磁波因為層內介質參數漸變而迅速衰減．后來有人提出的坐標伸縮完全匹配層(Coordinate stretched Perfectly Matched Layer， CPML)以及各向異性完全匹配層(Uniaxial anisotropy Perfectly Matched Layer， UPML)，均具有良好效果．

時域離散伽略金(Discontinuous Galerkin Time Domain， DGTD)方法[4-5]是基于時域體積元發展起來的一種兼備時域有限元(Finite Element Time Domian，FETD)方法網格剖分的靈活性和FDTD顯式迭代特點的算法．文獻[6]中將PML技術應用于DGTD，其后，文獻[7]將PML技術應用于共形的三維DGTD計算．文獻[8]提出一種DGTD提高PML穩定性的方案．文獻[9-10]中討論了高階結點DGTD方法中的PML問題．上述文獻主要討論三維問題，然而在很多的情形，例如金屬搭接縫中的電波傳播問題、層狀介質中的波傳播問題、地面上方長傳輸線的電磁場分布、地下長線纜的電磁場分布問題和光纖里的光傳輸問題等都可以簡化成二維問題來處理．考慮到二維方案在解決許多問題中方便快捷的特性，系統地研究DGTD算法的二維UPML吸收邊界具有重要的意義．

文中從各向異性介質中的Maxwell方程組和UPML基本理論出發，得到UPML介質層中的支配方程和輔助方程．考慮到DGTD空間離散網格的非結構特性，結合離散單元之間場量傳遞的特點，將方程加權積分得到單元矩陣方程．進一步，考慮相鄰單元之間電磁場矢量切向分量突變的通量邊界條件，通過基函數展開和矩陣單元計算，得到具有通量形式的時域公式，進而得到橫磁(Transverse Magnetic，TM)波情形UPML中的時域迭代計算式．數值結果說明，文中UPML吸收邊界對電磁波的衰減達 130 dB，吸收效果良好，滿足二維情形電磁場精確仿真的需求．

1 二維UPML基本理論

UPML是指在計算域外部增加一層非物理區域，該區域中為各向異性介質．實際操作中適當選擇各向異性介質的本構參數，PML層波阻抗與相鄰區域匹配，使入射到分界面上的電磁波無反射地進入該區域．在PML層中，電場和磁場滿足的含輔助變量的支配方程為[11]

其中，μ1和ε1為主要計算域中的介電系數和磁導率；a、b和c均為對角張量(具體表示為 a= diag(axx，ayy，azz)，其余類似)；輔助變量Ph和Pe與電場E和磁場H之間滿足的微分方程為

其中，d和κ同樣為對角陣．UPML實現過程中，需要獲得式(1)和式(2)的分量式并離散．

對于TM波，關注的電磁場分量為Ez、Hx和Hy．因此，式(1)可轉化為

其中，輔助函數Ph和Pe的分量式為

將式(4)代入式(3)，得TM波情形支配方程式(1)所滿足的分量方程，即

方程式(2)的分量式為

圖1 UPML區劃分示意圖

分量方程式(5)和式(6)在具體實現數值計算時，需要改寫為離散的時域遞推公式．離散方案依賴于空間離散單元和基函數的選取等．在FDTD方法中常采用六面體離散計算域，并采用點匹配法，而DGTD與FETD類似常采用三維四面體或二維三角形離散單元和伽遼金法．由于采用非結構性網格離散，DGTD的UPML實現更為復雜．

2 二維情形UPML離散DGTD公式

UPML區域劃分如圖1所示．在包含目標的計算域外部，吸收邊界被劃分成為4個棱邊區和4個平面區．

不同于FETD中單元交界面處采用電磁場連續的強制邊界條件，DGTD采用數值通量的方式實現網格之間的場量交互．TM波情形在邊界上交互時數值通量定義為[5]

在式(7)通量定義的基礎上，通過加權積分[12]和基函數展開，可得分量方程式(5)對應的矩陣方程，即

其中，上標m為單元編號; Mm為單元質量矩陣；Sx和Sy為剛性矩陣；Fe，z、Fh，x和Fh，y分別指電場e和磁場h的通量變化量函數．

將式(8)和式(9)在t=nΔt處離散，式(10)在t=(n+1/2) Δt處離散，可得

其中，αgx、βgx、αgy、βgy、αgz和βgz為常系數．

輔助方程式(6)無需加權積分，直接離散可得

其中，ch、ci、cj和ck為常系數．

式(11)～(16)就是二維TM波情形DGTD算法中UPML區域的遞推計算公式．

3 數值算例

設二維計算域為8 m×8 m的矩形，區域離散為 10 902 個結點，共 21 474 個三角形單元，UPML層厚為 1 m．線源在計算域中心處．

UPML層中電導率σ為非均勻分布，其漸變公式為

σ=σmaxd-d0(Δd)mm,

圖2 線電流源在自由空間中的輻射電場

其中，d為UPML區中單元的重心坐標位置；d0為UPML內側界面位置；Δd為UPML層厚度；m為整數，本算例中 m=4．

計算時，時間離散間隔Δt=0.16×10-10s，當輻射電磁波的頻率 f=0.3 GHz 時，文中方法的計算結果與利用Hankel函數得到的解析解的比較如圖2所示．圖2中三角形和實線分別表示文中算法和解析解的結果，由圖可見，兩種算法的計算結果一致，這說明了文中吸收邊界算法的有效性．

為更好地說明電磁波進入UPML層內的衰減情況，下面給出在輻射線源和計算域最外端連線上場值的分布情形．圖3(a)是線元輻射的場值快照，線源位于A點處．線段AB上的電場如圖3(b)所示，由圖可見，電磁波在剛進入UPML層的階段衰減不大(這是由于界面附近參數接近于內部計算域，便于電磁波無反射地進入吸收層)．隨著電磁波深入UPML層，衰減量逐步加大(吸收層參數設置漸變所致)．吸收層界面上場值約為 -25 dB，而到吸收層外側已下降到 -90 dB 以下．場值衰減為入射場的 1/500，考慮到返回計算域的電場值還需經過相同的路徑，這一衰減的效果將加倍．

圖3 線段AB處的電場Ez

4 結 束 語

文中從二維UPML的基本理論出發，根據空間離散網格的非結構特性，結合DGTD中離散單元之間場量傳遞的特點，在TM波情形下推導了DGTD各向異性完全匹配層吸收邊界中的迭代計算式，實現了波在UPML層中的衰減．仿真算例說明，文中給出的DGTD UPML吸收邊界對電磁波的雙重衰減達 130 dB，吸收效果良好，滿足二維情形電磁場精確仿真的需求．

[1] 王林年, 褚慶昕, 梁昌洪. 交替方向隱式時域有限差分法中的Berenger理想匹配層[J]. 西安電子科技大學學報, 2006, 33(3): 458-461.

WANG Linnian, CHU Qingxin, LIANG Changhong. Berenger’s Perfectly Matched Layer for the Alternating Direction Implicit Finite-difference Time-domain Method[J]. Journal of Xidian University, 2006, 33(3): 458-461.

[2]姜彥南,葛德彪,楊利霞.衛星模型散射FDTD計算的共形邊界研究[J]. 西安電子科技大學學報, 2007, 34(4): 587-589.

JIANG Yannan, GE Debiao, YANG Lixia. Conformal Boundary Scheme in FDTD Computation of Satellite Model Scattering[J]. Journal of Xidian University, 2007, 34(4): 587-589.

[3]BERENGER J P. A Perfectly Matched Layer for the Absorption of Electromagnetic Waves[J]. Journal of Computational Physics, 1994, 114(2): 185-200.

[4]CHEN J, LIU Q H. Discontinuous Galerkin Time-domain Methods for Multiscale Electromagnetic Simulations: a Review[J]. Proceedings of the IEEE, 2013, 101(2): 242-254.

[5]ALVAREZ J. A Discontinuous Galerkinfinite Element Method for the Time-domain Solution of Maxwell Equations[D]. Granada: University of Granada, 2013: 31-39.

[6]XIAO T, LIU Q H. Three-dimensional Unstructured-grid Discontinuous Galerkin Method for Maxwell’s Equations with Well-posed Perfectly Matched Layer[J]. Microwave and Optical Technology Letters, 2005, 46(5): 459-463.

[7]STYLIANOS D, LEE J F. Interior Penalty Discontinuous Galerkin Finite ElementMethod for the Time-dependent First Order Maxwell’s Equations[J]. IEEE Transactions on Antennas and Propagation, 2010, 58(12): 4085-4090.

[8]PENG D, CHEN L, YIN W L, et al. An Efficient DGTD Implementation of the Uniaxial Perfectly Matched Layer[C]//Proceedings of the 2013 International Symposium on Antennas & Propagation: 2. Piscataway: IEEE, 2013: 1272-1275.

[9]JI X, CAI W, ZHANG P. Reflection/Transmission Characteristics of a Discontinuous Galerkin Method for Maxwell’s Equations in Dispersive Inhomogeneous Media[J]. Journal of Computational Mathematics-international Edition, 2008, 26(3): 347-364.

[10]NIEGEMANN J, K?NIG M, STANNIGEL K, et al. Higher-order Time-domain Methods for the Analysis of Nano-photonic Systems [J]. Photonics and Nanostructures-Fundamentals and Applications, 2009, 7(1): 2-11.

[11]GEDNEY S D, KRAMER T, LUO C, et al. The Discontinuous Galerkin Finite Element Time Domain Method (DGFETD)[C]//IEEE International Symposium on Electromagnetic Compatibility. Piscataway: IEEE, 2008: 1-4.

[12]JIN J M. The Finite Element Method in Electromagnetics [M]. 3rd Edition. New York: John Wiley &Sons, 2014: 37.

(編輯：齊淑娟)

Implementation of UPML absorbing boundary conditions in DGTD method

LILinqian1,2，WEIBing1,2，YANGQian1,2，GEDebiao1,2

(1. School of Physics and Optoelectronic Engineering, Xidian Univ., Xi’an 710071, China; 2. Collaborative Innovation Center of Information Sensing and Understanding, Xidian Univ., Xi’an 710071, China)

For the DGTD-UPML absorb boundary condition,the basic theory of UPML in the TM case is studied and implementationis presented. By lessons from the ideas of impedance matching, anisotropic parameters set and matched matrix in the FDTD method, based on the characteristics of the unstructured meshing and features of flux transfer among the elements, the formulas for wave propagation in the UPML are presented, and then the iterative formulas are derived. Compared with the generalreflectivity of -24 dB in the first step SMABC due to the discretization of the unstructured grid in the DGTD, the double attenuation of 130 dB demonstrates the good absorption of the 2D UPML layer for the electromagnetic wave.Key Words: discontinuous Galerkin time domain(DGTD); uniaxial anisotropy perfectly matched layer(UPML); TM wave

2015-10-05

時間：2016-04-01

國家自然科學基金資助項目(61231003，61571348，61401344)

李林茜(1985-)，男，西安電子科技大學博士研究生，E-mail: 395106835@qq.com．

魏 兵(1970-)，男，教授，E-mail: bwei@xidian.edu.cn．

http://www.cnki.net/kcms/detail/61.1076.tn.20160401.1622.030.html

10.3969/j.issn.1001-2400.2016.06.015

O441.4

A

1001-2400(2016)06-0086-05