Comparison of two discontinuous spectral element methods
-
摘要: 时域间断元方法是近年来电磁场计算领域的重要进展之一。将基函数的插值点和数值积分点重合的质量集中技术是降低该间断元方法质量矩阵存储开销和提高计算效率的重要手段。通过谐振腔、带通滤波器以及光子晶体内的电磁场等数值算例,在四边形网格上比较了传统的质量集中算法和近来提出的Weight-Adjust算法之间的差异。计算结果表明,尽管两种方法的存储量一样,但Weight-Adjust算法具有更高的精度。Abstract: This paper is concerned with the comparison of two spectral Discontinuous Galerkin Time-Domain (DGTD) methods for the solution of Maxwell's equations in two-dimensional space. The first scheme is based on the conventional mass-lumping technique, where the same set of points are chosen for both the interpolation base functions and the numerical integration of coefficients. The second scheme is a newly proposed approach, called the Weight-Adjusted discontinuous Galerkin (WADG) method. Several numerical examples are presented to evaluate the performance of the two methods. It is shown that although the two methods need the same storage capacity, the WADG method has higher accuracy.
-
Key words:
- discontinuous Galerkin method /
- Maxwell equations /
- mass-lumping /
- spectral method
-
间断伽辽金有限元(DG)方法是近年来计算电磁学领域发展最为迅速的方法之一[1-2]。它具有精度高[3]、不规则网格易处理[4]、时间步长选取灵活[5]、并行计算易实现[6]等优点,在雷达、通讯、微波、集成光学等领域有重要的应用[7-9]。为了提高间断元方法的计算精度,需要加密模型划分的单元格(h方法)或者提高基函数插值阶数(p方法)。一般来说, p方法要比h方法收敛更快。基函数阶数的增加必然会带来质量矩阵的增大,从而使存储开销急剧增长(增长速度与阶数的三次方成正比)。为了解决这一困难,人们提出了质量集中算法。其基本思想是,将基函数的插值点与控制方程弱形式中的数值积分点重合,从而达到增加质量矩阵中零元素的目的。对标量基函数而言,可以使质量矩阵对角化,而对矢量基函数,质量矩阵可化为3×3的块对角矩阵。显然,对两种类型的基函数,质量集中都极大地降低了存储和计算开销。文献[10]介绍了一种新的质量集中算法,并将其成功应用到波动问题上。关于该方法在电磁学领域中的应用,以及新算法与传统质量集中算法之间的差异性需要进一步探讨。本文通过数值算例对这两种算法进行了比较,重点考察了两种方法的精度和稳定特征,为质量集中算法在DG方法中的应用提供参考。
1. 间断有限元算法
1.1 二维问题的质量矩阵
本文仅讨论二维问题。TE波模式下(TM模式类似),无源区间的二维麦克斯韦方程组可以表示为
{ε∂Ex∂t=∂Hz∂yε∂Ey∂t=−∂Hz∂xμ∂Hz∂t=∂Ey∂x−∂Ex∂y (1) 式中: Ex和Ey分别是电场沿x轴和y轴的分量; Hz是磁场的分量; ε和μ分别是腔体内的电容率和磁导率。
将计算区域Ω划分为K个四边形网格单元,每个单元记作Ωk。单元Ωk内的电磁场量E和H可以近似为
Ekx(x,y,t)=N∑i=0ekx,i(t)Φi(x,y)Eky(x,y,t)=N∑i=0eky,i(t)Φi(x,y)Hkz(x,y,t)=N∑i=0hkz,i(t)Φi(x,y) (2) 式中: (ex, ik, ey, ik, hz, ik)为待求的电磁场量; Φ(x, y)是二维基矢量。将电磁场量的离散形式代入到方程组(1)中,用基矢量点乘公式两端,并在Ωk上积分,可以得到
εN∑i=0ϕki,j∂ekx,i∂t=N∑i=0ξki,jhkz,iεN∑i=0ϕki,j∂eky,i∂t=−N∑i=0δki,jhkz,iμN∑i=0ϕki,j∂hkz,i∂t=N∑i=0δki,jeky,i−N∑i=0ξki,jekx,i (3) 其中
(⋅,⋅)Ωk=∫Ωk⋅dsϕki,j=(Φki,Φkj)Ωkδki,j=(Φki,∂Φkj∂x)Ωk,ξki,j=(Φki,∂Φkj∂y)Ωk 两次使用分部积分可以进一步得到间断有限元法的强形式
N∑i=0(Mkij⋅∂ekx,i∂t−Sy,kij⋅hkz,i)=∫∂ΩkΦki(x,y)⋅ˆn×(Hk∗z−Hkz)dlN∑i=0(Mkij⋅∂eky,i∂t+Sx,kij⋅hkz,i)=∫∂ΩkΦki(x,y)⋅ˆn×(Hk∗z−Hkz)dlN∑i=0(Mkij⋅∂hkz,i∂t−Sx,kij⋅eky,i+Sy,kij⋅ekx,i)=∫∂ΩkΦki(x,y)[ˆn×(Ek∗x−Ekx)+ˆn×(Ek∗y−Eky)]dl (4) 式中: ˆn为∂Ωk外向法向量,(Hzk*, Exk*, Eyk*)是边界∂Ωk上的数值通量。
Mk表示该单元上的质量矩阵,Sx, k, Sy, k表示相应的刚度矩阵,具体表示为
Mkij=C∫ΩkΦki(x,y)Φkj(x,y)dsSx,kij=∫ΩkΦki(x,y)∂Φkj(x,y)∂x dsSy,kij=∫ΩkΦki(x,y)∂Φkj(x,y)∂y ds (5) 式中: C表示ε或者μ。
由于积分单元Ωk并不是一个标准的四边形,所以我们使用坐标变换,通过雅克比行列式|Jk|将其转换到[-1, 1]2的正方形单元ˆΩ上积分,质量矩阵的元素表示为
Mkij=∫ˆΩCΦkiΦkj|Jk|ds=CN′∑α=0ωβN′∑β=0ωαΦi(xα,yβ)Φj(xα,yβ)|Jk| (6) 式中: ωα, ωβ是对应积分点的权重。
通过观察Mijk的形式可以发现,对于选用p=N阶插值多项式的单元,需要一个(N+1)2×(N+1)2的方阵来存储。由于系数C和雅克比行列式|Jk|不同,每个单元的Mk也不尽相同。因此,对于拥有n个单元的网格,程序中总共需要n个(N+1)2×(N+1)2的二维数组来存储质量矩阵。
1.2 两种方法质量矩阵存储量比较
为了减少质量矩阵带来的内存开销,传统的质量集中间断有限元算法是选用相同的插值点和积分点来近似计算,则二维基矢量Φ (x, y)表示为
Φij(x,y)=ϕNi(x)ϕNj(y)=N∏k=0,k≠ix−xkxi−xk⋅N∏k=0,k≠jy−ykyj−yk (7) 当插值点与积分点重合时有
N∑α=0ωβN∑β=0ωαΦi(xα,yβ)Φj(xα,yβ){≠0,i=j=0,i≠j (8) 由此方法得到的质量矩阵Mk,是一个(N+1)2阶的对角矩阵,只需要用一个长度为(N+1)2的一维数组来存放其对角线上的元素,大大减少了存储空间的损耗。
文献[10]中介绍的Weight-Adjust算法与之类似,在其基础之上给出了另外一种ˆΩ单元上逆质量矩阵(Mk)-1的近似求法,可以表示为
(Mk)−1≈M−1M1/ε|Jk|M−1 (9) 其中
Mij=∫ˆΩΦiΦj ds (10) (M1/ε|Jk|)ij=∫ˆΩεΦiΦj|Jk|ds=N∑β=0ωβN∑α=0ωaΦi(xα,yβ)Φj(xα,yβ)ε|Jk| (11) 通过观察式(11)我们可以发现,M矩阵同样需要一个(N+1)2阶的方阵来存储,但是由于每个参考单元上的M矩阵均相同,拥有n个单元的网格可以共用同一个(N+1)2阶方阵,造成的开销可以忽略不计。当选用相同的插值点和积分点时,矩阵M1/ε|Jk|是一个(N+1)2阶的对角矩阵,只需要用一个长度为(N+1)2的一维数组来存放其对角线上的元素。
通过上面的分析,两种近似方法所需要的质量矩阵存储量相同,下文的数值算例部分则从精度以及计算效率等方面比较这两种方法。
2. 数值算例和讨论
2.1 谐振腔
首先考虑计算一个有解析解的方形谐振腔问题,来分析两种方法的精度以及计算效率的差异。设计一个[0, 1]2 m2的方形区域,划分为12个不规则的四边形单元。在区域内给定基模作为初始值后,两种方法均使用相同的时间步长(时间积分使用Runge-Kutta法[4])和6阶插值多项式进行模拟,并记录测量点(0.5, 0.85)处磁场分量Hz随着时间变化的值。最后分别将两种方法模拟得到的数据与参考值进行比较,得到绝对误差。
在图 1(a)中,给出了测量点处的参考时域波型以及两种近似方法下得到的实际波型,图 1(b)中进一步得到测量点处数值解与参考解比较的绝对误差。可以明显地看出,传统DG算法模拟结果的误差振幅大约为WADG算法的两倍。
另外,通过简单地调整程序的时间步长,我们发现在保证模拟结果收敛于解析解的前提下,传统的质量集中算法可以使用的最大时间步长是Weight-Adjust算法的4倍。
2.2 滤波器
下面我们将用上述两种方法去模拟计算一个带通滤波器的散射参数,具体的截面模型如图 2所示,各项参数均参考文献[11]。整个截面模型包括左右两端的PML吸收层,并且被划分为3488个长方形单元。
滤波器的散射系数可以表示为
S11=√pref pinc ,S21=√ptra pinc (12) 为了方便计算S参数所需的系数,我们添加了TF/SF界面,并在该界面上设置了TE10的模式波。通过电磁场的能流密度公式,我们求得TF/SF面左侧面上的反射系数pref,右侧面上的发射系数pinc以及滤波器与右侧PML交界面上的透射系数ptra。
我们将最后模拟结果,即不同频率模式波下的散射系数表示在图 3中。两种方法在长方形网格下模拟得到的散射系数结果差别较小,并且与文献[11]中给出的模拟结果基本一致。
2.3 光子晶体
近年来,间断有限元法在光子晶体领域的应用比较广泛,为了进一步探讨两种算法处理更加复杂的问题时的性能,我们将其应用到模拟光子晶体功率分配器内的电磁场[12]。如图 4所示,整个截面被划分为3488个四边形单元,每个单元均选用4阶插值多项式。模型的截面尺寸如图 5所示,包括一层厚度为1 μm的PML吸收层以及其最外侧的PEC边界。
光子晶体中圆柱体的半径为0.18 μm,其材料的相对电容率和相对磁导率分别为ε(r)=11.56和μ(r)=1。相邻圆柱体中心的间距为1 μm,其中的填充介质相对介电常数和相对磁导率分别为ε(r)=1和μ(r)=1。在图中标注的位置我们给定了一个高斯脉冲E(r,t)=ˆzG(t),其中G(t)=e-t2/2τ2cos(ωt),频率ω=2.513 3×106c(c是光速),时延τ=1.591 5。
分别用两种方法进行模拟,记录下每点处的电磁场大小,并进行绘制。这里给出了Weight-Adjusted算法下,光子晶体功率分配器截面电磁场模拟结果,如图 5中所示,本文模拟得到的电磁场分布与文献[12]给出的参考分布相同。当然,传统的质量集中算法也能得到相似的模拟结果,这里便不再使用图片说明。
为了进一步说明两种方法在模拟精度上的差异,我们选取分配器通道内的一点,分别记录两组模拟中该点处的时域波型,并与参考数据进行比较(参考数据由高阶间断有限元法模拟得到)。从图 6给出的两种模拟方法各自的绝对误差中可以看出,Weight-Adjusted算法获得的数据精确度更高。
3. 结论
Weight-Adjusted算法是一种能够简化质量矩阵,从而节省内存开销的近似算法。本文首次尝试将其应用到处理四边形网格上的二维时域电磁场问题,并与该领域内一种已经成熟使用的质量集中算法在精度、内存、计算效率等方面进行了比较。谐振腔、带通滤波器以及光子晶体内电磁场的数值算例的结果表明,Weight-Adjusted算法的精度明显更高。
-
-
[1] Stephane D, Clement D, Stephane L, et al. Recent advances on a DGTD method for time-domain electromagnetics[J]. Photonics and Nanostructures Fundamentals and Applications, 2013, 11(4): 291-302. doi: 10.1016/j.photonics.2013.06.005 [2] Stephane D, Stephane L, Ludovic M. Temporal convergence analysis of a locally implicit discontinuous Galerkin time domain method for electromagnetic wave propagation in dispersive media[J]. Journal of Computational and Applied Mathematics, 2017, 316(15): 122-132. [3] Xiong Meng, Jennifer K R. Discontinuous Galerkin methods for nonlinear scalar hyperbolic conservation laws: divided difference estimates and accuracy enhancement[J]. Numerische Mathematik, 2017, 136(1): 27-73. [4] Zhu Jun, Zhong Xinghui, Shu Chiwang, et al. Runge-Kutta discontinuous Galerkin method with a simple and compact hermite WENO limiter on unstructured meshes[J]. Communications in Computational Physics, 2017, 21(3): 623-649. [5] Demirel A, Niegemann J, Busch K, et al. Efficient multiple time-stepping algorithms of higher order[J]. Journal of Computational Physics, 2015, 285(15): 133-148. [6] 夏轶栋, 伍贻兆, 吕宏强. 高阶间断有限元法的并行计算研究[J]. 空气动力学报, 2011, 29(5): 537-541. https://www.cnki.com.cn/Article/CJFDTOTAL-KQDX201105000.htmXia Yidong, Wu Yizhao, Lü Hongqiang. Parallel computation of a high-order discontinuous Galerkin method on unstructured grids. Acta Aerodynamica Sinica, 2011, 29(5): 537-541 https://www.cnki.com.cn/Article/CJFDTOTAL-KQDX201105000.htm [7] Liu Xiang, Kameni A, Serhir M, et al. 3D modeling with discontinuous Galerkin time domain method for ground penetration radar application[C]//GDR Ondes, Assemblee Generale Interferences d'Ondes, 2015. [8] Su Yan, Andrew D G, Jin Jianming. Modeling of plasma formation during high-power microwave breakdown in air using the discontinuous Galerkin time-domain method[J]. IEEE Journal on Multiscale and Multiphysics Computational Techniques, 2016, 27(1): 2-13. [9] Nikolai S, Claire S, Stephane L, et al. A DGTD method for the numerical modeling of the interaction of light with nanometer scale metallic structures taking into account non-local dispersion effects[J]. Journal of Computational Physics, 2016, 316(7): 396-415. [10] Chan J, Hewett R J, Warburton T. Weight-adjusted discontinuous Galerkin methods: wave propagation in heterogeneous media[OL]. arXiv, 2016. [11] Anderson W K, Wang Li, Kapadia S, et al. Petrov-Galerkin and discontinuous Galerkin methods for time-domain and frequency-domain electromagnetic simulations[J]. Journal of Computational Physics, 2011, 230(23): 8360-8385. [12] Liu Meilin, Sirenko K, Bagci H. An efficient discontinuous Galerkin finite element method for highly accurate solution of Maxwell equations[J]. IEEE Trans Antennas and Propagation, 2012, 60(8): 3992-3998. 期刊类型引用(3)
1. 段书超, 王刚华, 谢卫平, 阚明先, 李晶, 邹文康, 徐强, 但加坤. FOI-PERFECT程序对电磁驱动高能量密度系统的三维弛豫磁流体力学模拟. 强激光与粒子束. 2016(04): 146-152 . 本站查看
2. 章征伟, 魏懿, 孙奇志, 刘伟, 赵小明, 张朝辉, 王贵林, 郭帅, 谢卫平. 材料强度对电磁驱动固体套筒内爆过程的影响. 强激光与粒子束. 2016(04): 162-166 . 本站查看
3. 沈双晏, 金星. 磁流体动力学磁控进气道流场分布数值模拟. 强激光与粒子束. 2015(12): 214-222 . 本站查看
其他类型引用(0)
-