② 西安电子科技大学信息感知技术协同创新中心 西安 710071
② Collaborative Innovation Center of Information Sensing and Understanding, Xidian University, Xi’an 710071, China
在目标特性的精确预估中,矩量法(Method of Moment,MoM)、有限元法(Finite Element Method,FEM)和时域有限差分方法(Finite Difference Time Domain,FDTD)均是应用广泛的计算方法。MoM主要长于金属目标散射问题的处理,FEM和FDTD方法则更适宜于复杂媒质的计算。MoM和FEM需要求解大型线性方程组,应用于电大尺寸目标时求解效率低。FDTD采用显式迭代的方式,易于并行计算,并且计算复杂度低,但其在对目标拟合上存在台阶误差,即使采取亚网格或共形网格也很难精确拟合。近年来,离散伽辽金法(Discontinuous Galerkin method)算法受到了人们的重视[1]。该算法在流体力学及有限体积法(Finite volume method)中早有应用[2],近年来被引入有限元计算当中来形成了时域离散伽辽金算法(Discontinuous Galerkin Time Domain,DGTD)[3, 4]。DGTD结合了FDTD显式迭代和FEM网格精确拟合目标几何性形状的特点,是一种有广阔应用前景的计算方法。近十余年,国际上的同行对该方法已经进行了多方面的研究,但在国内将基于矢量基函数的DGTD方法用于目标散射特性的分析还未见相关报道。本文首先简要介绍了基于棱边矢量基函数的3维DGTD方法的基本思想。然后,给出了在该方法实现中的平面波引入及总场边界条件。最后,给出目标单站散射特性的算例表明本文算法的正确性。
2 基于棱边基函数的3维DGTDDGTD可以看作由时域有限元方法(Finite Element Time Domain,FETD)改变边界条件的处理方式得到。FETD是从支配方程和边界条件出发将计算区域划分为多个单元后导出矩阵方程并求解的方法,J. M. Jin,J. F. Lee等人在FETD方面做了很多工作[5, 6, 7]。一般来说,FETD方程推导有两种途径,即变分法与Galerkin法。变分法寻找一个适当泛函对应于支配方程和边界条件,Galerkin法通过加权余量来寻找支配方程和边界条件的相应弱解形式。两种方法都会得到一个大型矩阵方程。在FETD中的若干基本概念如基函数、单元积分、局域、全域、Galerkin加权等均可在DGTD中继承使用,而两种有限元算法的局域边界处理方式有差异。
下面从Maxwell方程组出发描述DGTD的基本思想,Maxwell旋度方程为:
ε∂E∂t−∇×H+σE=−Jμ∂E∂t+∇×E+σmH=−Jm} | (1) |
本文中采用四面体作为DGTD算法的基本离散单元(如图 1所示),采用Galerkin加权法[6],对式(1)在四面体单元内积分,可得
![]() |
图 1 四面体单元示意图 Fig.1 Schematic diagram of tetrahedron |
∫v⋅[ε∂E∂t−∇×H+σE+J]dΩ=0∫v⋅[μ∂E∂t+∇×E+σmH+Jm]dΩ=0} | (2) |
∫v⋅(∇×E)dΩ=∫∇⋅(E×v)dΩ+∫E⋅(∇×v)dΩ=∫∧n⋅(E×v)ds+∫E⋅(∇×v)dΩ=∫v⋅(∧n×E)ds+∫E⋅(∇×v)dΩ∫v⋅(∇×H)dΩ=∫∇⋅(H×v)dΩ+∫H⋅(∇×v)dΩ=∫∧n⋅(H×v)ds+∫H⋅(∇×v)dΩ=∫v⋅(∧n×H)ds+∫H⋅(∇×v)dΩ} | (3) |
在流体力学中为求解双曲方程,边界条件采用Numerical Flux[2],后又被引入了时域有限体积法(FVTD)[8],又被引入DGTD作核心思想[4]。在边界上定义E*与H*为:
∧n×E∗=∧n×E+KeE[∧n×(E+−E)]+veH[∧n×(∧n×(H+−H))]∧n×H∗=∧n×H+KeH[∧n×(H+−H)]−veE[∧n×(∧n×(E+−E))]} | (4) |
keH,veH,veE的取值如表 1[4]。
![]() |
表 1 Numerical Flux系数表达式 Tab. 1 Coefficient of Numerical Flux |
表 1中Ye=√εe/μe=1/Ze,e为当前单元,e+为相邻单元。
由于E*,H*仅定义在边界上,可得
∫v⋅(∇×E∗)dΩ=∫v⋅(∇×E)dΩ+∫{keE[∧n×(E+−E)]+veH[∧n×(∧n×(H+−H))]}ds∫v⋅(∇×H∗)dΩ=∫v⋅(∇×H)dΩ+∫{keH[∧n×(H+−H)]+veE[∧n×(∧n×(E+−E))]}ds} | (5) |
将式(5)代入式(2),
ε∫v⋅∂E∂tdΩ−∫v⋅(∇×H)dΩ−∫{keH[∧n×(H+−H)]ds+∫veE[∧n×(∧n×(E+−E))]}ds+σ∫v⋅EdΩ+∫v⋅JdΩ=0μ∫v⋅∂H∂tdΩ+∫v⋅(∇×E)dΩ+∫{keE[∧n×(E+−E)]ds+∫veH[∧n×(∧n×(H+−H))]}ds+σm∫v⋅HdΩ+∫v⋅JmdΩ=0} | (6) |
经整理后得矩阵形式的偏微分方程组,此方程组针对单个四面体,实际计算中在计算域内循环所有四面体迭代求解。
ε[Me]∂∂t{Ee}−[Se]{He}−{Fhe}+σ[Me]{Ee}+{Je}=0μ[Me]∂∂t{He}+[Se]{Ee}+{Fee}+σm[Me]{He}+{Jem}=0} | (7) |
本文采用Whitney-I型矢量基函数[9],Me,Se,Fhe,Fee,Je,Jm矩阵元素表达式如下:
Meij=∭ΩeNei⋅NejdΩ,Seij=∭ΩeNei⋅(∇×Nej)dΩFhei=∑[keH]⋅{He+j−Hej}−∑[νeE]⋅{Ee+j−Eej}=∑keH[Mfeij]⋅{He+j−Hej}−∑νeE[Mgeij]⋅{Ee+j−Eej}Feei=∑[keH]⋅{Ee+j−Eej}+∑[νeH]⋅{He+j−Hej}=∑keE[Mfeij]⋅{Ee+j−Eej}−∑νeH[Mgeij]⋅{He+j−Hej}Jei=∭ΩeNei⋅⋅JeidΩ,Jemi=∭ΩeNei⋅⋅JemdΩ} | (8) |
式(8)中的Nei,Nej为基函数,
3 DGTD中的平面波引入及总场边界条件用DGTD计算散射问题时可将计算区域划分为总场区和散射场区,如图 2所示。采用等效原理,类似FDTD中的总场边界条件,在总场边界上设置入射波的切向分量即可将入射波引入到总场区[10]。
![]() |
图 2 总场区与散射场区的划分 Fig.2 The total field region and the scattered field region |
为便于说明总场边界条件的加入,图 3中以三角形代替四面体,规定散射场区单元的场均属于散射场,总场区单元的场均属于总场。在总场区内部及散射场区内部,DGTD计算方式如前文所述,需要特殊处理的是边界处相邻的单元计算式,具体来讲是在Numerical Flux表达式中增加入射波的计算部分。
![]() |
图 3 DGTD总场散射场边界 Fig.3 TF/SF boundary for DGTD |
(1) 当四面体处于总场边界的总场区一侧时,应在Numerical Flux上加上入射波值,此时Eej,Hej属于总场,Ee+j,He+j属于散射场,Ei,Hi为入射波。
Fhej=∑keH[Mfeij]⋅{He+j−Hej+Hi}−∑veE[Mgeij]⋅{Ee+j−Eej+Ei}Feej=∑keE[Mfeij]⋅{Ee+j−Eej+Ei}+∑veH[Mgeij]⋅{He+j−Hej+Hi}} | (9) |
(2) 当四面体处于总场边界的散射场区一侧时,应在Numerical Flux上减去入射波值,此时3Eej,Hej属于散射场,Ee+j,He+j属于总场,Ei,Hi为入射波。
Fhej=∑keH[Mfeij]⋅{He+j−Hej−Hi}−∑veE[Mgeij]⋅{Ee+j−Eej−Ei}Feej=∑keE[Mfeij]⋅{Ee+j−Eej−Ei}+∑veH[Mgeij]⋅{He+j−Hej−Hi}} | (10) |
采取上面的方法即可将平面波引入总场区域。
4 算例为说明本文算法的正确性,用DGTD方法计算半径为1 m的金属球的后向RCS并与解析解结果相对比。计算区域被离散成为319308个四面体(计算域剖分截面图如图 4所示),空间离散尺度为0.15 m,采用1阶Mur吸收边界。图 5是金属球后向单站RCS的计算结果(圆圈),作为比较,图中还给出了Mie级数计算结果(实线),由图可见两者的结果吻合较好。本算例DGTD消耗内存约280 MB,计算耗时约20 min,计算平台为Intel Core i5 3470;采用FETD计算所需内存消耗为2 G,计算时间在6 h以上。显然,本文方法有更高的计算效率。
![]() |
图 4 计算域剖分截面图 Fig.4 Sectional view of the computational domain |
![]() |
图 5 金属球单站RCS Fig.5 The monostatic RCS of PEC sphere |
设介质球参数为e=4e0,m=m0,s=0.00015 S/m,其余参数同上。图 6是本文方法(圆圈)与Mie级数计算结果(实线)的比较,由图可见,两者相吻合。
![]() |
图 6 介质球单站RCS Fig.6 The monostatic RCS of medium sphere |
为说明本文算法对于计算较为复杂目标的适用性,下面计算金属弹头的后向RCS。弹头模型如图 7所示(底面半径为1 m,高为2 m)。电磁波沿z轴负方向,电场沿x方向极化入射。图 8为弹头的后向RCS,实线为DGTD算法的结果,斜十字为MoM算法扫频结果。由图可见,两种算法的结果符合较好。
![]() |
图 7 弹头模型 Fig.7 Model of PEC bullet |
![]() |
图 8 金属弹头单站RCS Fig.8 The monostatic RCS of PEC bullet |
在散射问题的DGTD的分析中吸收边界条件、时间离散方式、近远场外推和四面体相邻单元的快速判断等方面都需要精细考虑,限于篇幅,本文没有讨论这些内容。目前,由于吸收边界采用1阶Mur吸收边界,本文算例中DGTD的计算精度有待提高。采用UPML吸收边界、高阶基函数及精度更高的时间离散方案将是我们下一步的工作。可以预期,作为一种新兴算法,DGTD将在目标特性的精确预估,特别是含复杂媒质的低RCS目标特性的精确预估方面扮演重要角色。
[1] | Ji X, Lu T, Cai W, et al.. Discontinuous Galerkin Time Domain (DGTD) methods for the study of 2-D waveguidecoupled microring resonators[J]. Journal of Lightwave Technology, 2005, 23(11): 3864-3874.(![]() |
[2] | Shu C W. A brief survey on discontinuous Galerkin methods in computational fluid dynamics[J]. Advances in Mechanics, 2013, 43(6): 541-553.(![]() |
[3] | Gedney S D, Kramer T, Luo C, et al.. The Discontinuous Galerkin Finite Element Time Domainmethod(DGFETD)[C]. IEEE International Symposium on Electromagnetic Compatibility, 2008: 1-4.(![]() |
[4] | Alvarez J. A Discontinuous Galerkin Finite Element Method for the Time-Domain Solution of Maxwell Equations[D]. [Ph.D. dissertation], University of Granada (Spain), 2013: 31-39.(![]() |
[5] | Lee Jin-Fa, Lee Robert, and Cangellaris Andreas. Timedomain finite-element methods[J]. IEEE Transactions on Antennas and Propagation, 1997, 45(3): 430-442.(![]() |
[6] | Jin Jian-ming.The Finite Element Method in Electromagnetic[M]. New York: John Wiley & Sons, 2002: 22-23.(![]() |
[7] | Riley D J, Jin Jian-ming, Lou Z, et al.. Total-and scatteredfield decomposition technique for the finite-element timedomain method[J]. IEEE Transactions on Antennas and Propagation, 2006, 54(1): 35-41.(![]() |
[8] | Shankar V, Mohammadian A H, and Hall W F. A timedomain, finite-volume treatment for the Maxwell equations[J]. Electromagnetics, 1990, 10(1/2): 127-145.(![]() |
[9] | Bossavit A. Whitney forms: a class of finite elements for three-dimensional computations in electromagnetism[J]. IEE Proceedings A (Physical Science, Measurement and Instrumentation, Management and Education, Reviews), 1988, 135(8): 493-500.(![]() |
[10] | 葛德彪, 魏兵. 电磁波时域计算方法(下册)[M]. 西安电子科技 大学出版社, 2014: 188-191. Ge D B and Wei B. Time Domain Computational Method for Electromagnetic Wave (Volume II)[M]. Xi'an: Xidian University Press, 2014: 188-191.(![]() |