时域有限元法在三维电磁辐射和散射问题中的应用

时域有限元法在三维电磁辐射和散射问题中的应用

时域有限元法近年来，通信系统和雷达系统的工作频率范围有所增加，这要求天线和相关微波设备具有频带特性。超宽带移动通信设备基本工作在微波波段,内部几何结构和材料分布十分复杂。很多先进设计还需要外加电子材料,如人工电磁带隙材料、左手材料等,以期进一步提高其性能。如此一来,对这一类具有复杂结构和非均匀介质的超宽带天线、超宽带雷达的瞬态辐射、目标的瞬态散射问题进行建模与仿真提出了新要求。十多年来,有大量的文献采用不同的时域分析计算方法,例如时域有限差分法(FDTD)、时域积分法(TDIE)和时域有限元法(TDFEM)等,对各种不同的结构进行了理论分析。FDTD法因其理论简单、方法通用而成为流行的数值模拟工具。然而,它不太适用于复杂几何结构和介电材料分布。TDIE法在分析无界均匀介质中的金属或均匀介质目标的散射和辐射时具有优势,原因有两方面:一是它把未知量限制在目标表面积上;二是,它通过格林函数自然地满足辐射条件。然而,当数值模拟的区域包含复杂非均匀介质时,TDIE法就会面临非均匀介质体内带来的大型矩阵皆为满阵而非稀疏矩阵,从而需要消耗大量的计算资源的困难。相反地,在众多的方法中,TDFEM因其继承了FEM特别适用于复杂的几何结构和介电特性分布的优点及在非均匀介质体内产生的也是稀疏矩阵的长处,因而在解决宽带特性的电磁场问题上弥补了FDTD和TDIE法的某些不足而被人们渐渐关注。时域有限元法起源于上世纪80年代,最初是由A.C.Cangellaris等人提出的微分形式的点匹配时域有限元法。该方法结合了FDTD显式积分法的简明性和传统有限元的灵活性,但是由于采用的是节点元,在强加边界条件时存在困难;又因同时采用电场和磁场,变量增加了一倍;此外还会出现源于电场和磁场点在网格上的交替出现而产生的类似FDTD法的越级(leap-frog)问题。为了避免或解决这些问题,1990年代提出无条件稳定的隐式时域有限元法[11,12,13,14,15,16],该方法基于电场(或磁场)的矢量有限元空间基函数和时间基函数展开的二阶矢量波方程,其优点是简单而稳定,缺点是每个时间步必须求解一个矩阵方程;对于开域空间问题,采用传统的吸收边界条件,有限元计算空间十分庞大,这个问题更加严重,以致无法进行。由于这些原因,这类时域有限元法发展缓慢。直到2000年代初,由于完全匹配层(PML)和正交基函数概念的提出和应用,时域有限元法才得到较大的发展。目前时域有限元法在国外,特别是以J.M.Jin的科研组在PML的应用、正交基函数的引进及其在二维和三维电磁辐射问题的应用方面,作了若干开创性的工作;在国内除外,仅见二维时域有限元法关于稳定性分析的报道。以下就时域有限元法应用于三维电磁辐射和散射问题的原理、发展中遇到的挑战和解决方案、尚待解决的问题和发展趋势作简要介绍,以期推动它在我国的应用和发展。全文共分5个部分,本节是引言;第二部分给出基本原理和公式;第三部分论述时域有限元法在应用于三维电磁问题时遇到的主要问题以及相应的解决方法;第四部分讨论时域有限元法在该领域的应用和发展的前景;最后一节为结论。1理论发展和公式1.1单元状态方程的匹配和点匹配时域有限元法起源于上世纪80年代,由A.C.Cangellaris等人首次提出,它将FDTD法作有限元插值扩展得到了节点时域有限元算法,称为点匹配FEM。很明显,这种方法的提出和发展的源动力是希望结合FDTD显式积分法的简明性和传统有限元空间离散过程的灵活性,其原理和公式如下:与FDTD方法相似,点匹配时域有限元法计算区域被2个互补的网格(分别是电场和磁场)离散化,网格的构造规则是使每个磁场网格的一个单元包含另一个单元的一个节点,反之亦然。按照有限元法的过程,电场E和磁场H在计算区域上可以被近似表达为:E(r,t)=Μ∑i=1ϕi(r)Ei(t),Η(r,t)=Ν∑j=1ϕj(r)Ηj(t)‚(1)这里M和N分别是电场和磁场节点的数目,Ei和Hj分别是电场和磁场在节点上的值,ϕi,ϕj是已知的有限元插值函数。利用(1)式可以表达在单元内部的任意一点的场变量值。然后,将Maxwell旋度方程化为下面的状态方程组:μdΗjdt=-L∑l=1∇φl(rj)×El(t)j=1,2,⋯Ν‚σEi+εdEidt=L∑l=1∇φl(ri)×Ηl(t)i=1,2,⋯Μ‚(2)其中φl(r)是有关单元节点l的单元形状函数,L是单元的节点数目。∇φl(rj)表示在r=rj处∇φl(r)的值。将(2)式中的时间微分用差分代替,就得到微分方程形式的点匹配TDFEM计算公式。值得注意的是此方法主要缺陷在于:每个场的三个分量都位于相同的节点上,这就导致了在媒质的发生电磁特性突变的边界表面上强加合适的边条件变得十分复杂。1.2whit充换电子束考虑到点匹配有限元法的缺陷,1990年之后有人提出另一途径得到的TDFEM,此类TDFEM是离散化二阶矢量波方程(也即旋度-旋度方程),从Maxwell’s方程组中消去一个场变量而得到的[11,12,13,14,15,16],其中典型的是StephenD.Gedney等人提出的一种无条件稳定的隐式FETD方法,其原理和公式简述如下:时变电场必须满足依赖时间的非齐次波动方程:∇×μ-1r∇×E+μ0σ∂E∂t+εrc20∂2E∂t2=-μ0∂J∂t(3)式中,c0是自由空间光速;σ是电导率;εr和μr是相对介电常数和磁导率。J是有限域Ω的表面∂Ω上的电流密度。所求问题可以简化成把∂Ω界定为是一个理想导电面(PEC)或理想导磁面(PMC)的电磁场边值问题。作(4)式与测试函数T的内积,应用Green第一恒等式可以导出(3)式的弱形式∭Ω[μ-1r(∇×Τ)(∇×E)+μ0σΤ∂E∂t+εrc20Τ∂2E∂t2]dΩ=-∭Ωμ0Τ∂J∂tdΩ(4)将矢量场用Whitney边棱元基函数作权重函数Wi展开,展开系数ej是时间连续函数,同时,测试函数也取为Whitney基。代回(4)式,并由此导出二阶微分方程[Τε]1c20d2edt2+[Τσ]η0c0dedt+[S]e=-f(5)这里[Tε]、[Tσ]和[S]是不依赖时间的矩阵:[Τε,σ]i,j=∭Ω{εrσ}Wi⋅WjdΩ(6)[S]i,j=∭Ω1μr∇×Wi⋅∇×WjdΩ(7)fi=∭Ωμ0Wi⋅∂J∂tdΩ(8)将(5)式中的时间微分用Newmark-β差分公式展开,可以得到积分形式的隐式TDFEM迭代公式,据此,由第n-1和第n个时间步的场值,可计算出第n+1个时间步的场值。2积分形式的隐式时域有限元法上一节中阐述的第一类直接求解Maxwell方程组的显式时域有限元法,优点是不需要在每个时间步上求解一个矩阵方程;但是它一般同时在电场和磁场上工作,变量增加了一倍,而且需要不同种类的基函数来展开它们;另外,由于电场和磁场点在网格上的交替出现,会导致类似FDTD法的越级(leap-frog)问题。由于这些原因,这一方法未得到推广应用,人们转而寻求另外的解决方案——积分形式的隐式时域有限元法。积分形式的隐式时域有限元法的优点是简单而且稳定,其缺点是只适用于理想导电面或理想导磁面;对于无界空间的开域问题,采用传统的吸收边界条件,会使计算区域变得十分庞大;同时,它在每个时间步上需要求解一个空间矩阵方程。由于上述原因,加之当时计算机速度和内存的限制,时域有限元法长期不能广泛流行。这种情况直到完全匹配层(PML)的出现和应用才得到改变。下面分两节阐述时域有限元法是如何解决上述问题的,首先阐述无界自由空间的截断问题——完全匹配层法,然后阐述庞大计算量问题的两个解决方案——正交基函数和子域分解并行计算。2.1时域有限元法微波辐射和散射问题的重要特点是开域,因此在用时域有限元法来分析和计算电磁场时,所面临的一个实际问题就是如何确定网格划分的边界和如何处理场在这些边界上的值。不像FDTD法在吸收边界条件(ABC)上有很好的发展,很长一段时间里,很少有对时域有限元法的精确ABCs进行研究的,也仅仅实现了一阶和二阶ABCs。直到把完全匹配层(perfectlymatchedlayer:PML)的概念引入到时域有限元法,这一问题才得到很好的解决。可以说,Berenger在1994年提出的完全匹配层把这项研究推向了新阶段。PML摆脱了以往吸收边界条件研究中拘泥于实际可实现的物理吸收边界的框框,提出了一种基于各向异性介质的数学边界。这种吸收边界的性能比以往提出的边界优越,而且为进一步研究能与不同介质相匹配的完全匹配层提供了新的思路。1995年,ZacharySacks等提出了一种易于实现和调控的PML。这种各向异性介质的性能参数[ε]和[μ]为复对角张量,为了使PML层的内表面成为一个完全的无反射界面,对角张量应具有以下形式:[μ]=μ0[Λ][ε]=ε0[Λ](9)Λ(r,ω)=diag{γyγzγx,γzγxγy,γxγyγz},γξ=1+σξ(ξ)jωε0ξ=x,y,z(10)这样电场在PML中满足∇×μ-1(Λ-1⋅∇×E)-ω2εΛ⋅E=0(11)(11)式看似简单,但是由于Λ是张量,而且是频率ω的函数,在进行时域计算时,必须利用Laplace或Fourier变换将(11)式转化为时域下的方程,数学实现十分繁琐。因此,虽然PML是在十多年前提出的,并已成为FDTD通常用的截断方法,并已应用于频域有限元法,但其实现于时域有限元法在2002年之后才有报道。为了得到(11)式对应的时域离散矩阵方程,我们采用Laplace变换及其反变换。同时,将电场根据空间基函数Nj展开,并利用测试基函数Ni得到时域版本的半离散的矩阵方程:[Τ]d2dt2{e}+[B]ddt{e}+[S]{e}-{h}+{g}=0(12)Τij=∭VεΝi⋅ΝjdV(13)Bij=∭Vεε0Νi⋅J⋅ΝjdV(14)Sij=∭V1μ(∇×Νi)⋅(∇×Νj)dV+∭Vεε20Νi⋅Κ1⋅ΝjdV(15)hi=∭Vεε30Νi⋅∑j(Κ2⋅uj⋅Νj)dV(16)gi=∭V1με0(∇×Νi)⋅∑j[(L1⋅u+j+L2⋅u++j)⋅(∇×Νj)]dV(17)其中J,K1,K2,L1,L2都是对角矩阵,矩阵元是σξ(ξ=x,y,z)的函数;并且uj,u+j,u++j是含有卷积的对角矩阵diage-σξε0t*ej(t),e-σξε0t*ej(t),e-σξε0t*ej(t)(18)其中σ的下标ξ代表{x,y,z}的轮换对称表达,在uj中分别对应于{x,y,z};uj+中分别对应于{y,z,x};而在uj++中分别对应于{z,x,y}。在利用PML和TDFEM求解电磁辐射和散射问题时,首先在辐射体或散射体的外侧设置封闭的PML层,在PML层内部,求解由(12)式所代表的TDFEM矩阵方程,在PML层内表面所包围包括辐射体和散射体内的空间,则求解由(5)所代表的TDFEM矩阵方程。由于PML层很薄(约n分之一波长)而且可以离辐射体或散射体很近,因此与ABCs相比极大的压缩了计算空间。2.2计算域和大型天线阵列尽管TDFEM应用了PML算法后大大缩小了计算区域,然而,对于模拟电大尺寸天线和大型天线阵列,计算域和未知场变量个数十分庞大的问题依然存在;另外在每个时间步上都要求解一个阶数极高的矩阵方程,这大大增加了计算的复杂度。为解决计算内存和时间问题人们进一步采取了各种办法,比较典型的有引入正交基函数和子域分解技术:2.2.1第i号边的面基函数有学者提出构造正交基函数的办法来解决在每个时间步上都要求解一个矩阵方程问题,通过正交基函数可以使基于Newmark-β的时域有限元法的隐式计算公式中的系数矩阵成为对角阵或高度稀疏的矩阵,从而迅速缩小每个时间步的计算时间。下面简要介绍它的构造原理:考虑一个矢量波动方程,用Whitney基函数对电场和磁场进行逼近。对于二维问题,原始的线性Whitney边基函数定义为Wi=ξj∇ξk-ξk∇ξj(19)这里的ξk是与节点k有关的线性节点基函数,且第i号边连接节点j和k,这些函数自然满足单元边界的切向场的连续性,对法向分量的不连续性没有影响。对于三维问题,面基函数一般定义为:Ui=2(ξj∇ξk×∇ξm+ξk∇ξm×∇ξj+ξm∇ξj×∇ξk)(20)其中,编号为i的面由节点j,k和m连接而成,ξj,ξk,ξm表示体积坐标。Ui的法向分量在单元边界是连续的,并且它的通量仅在穿过编号为i的面时为1,其他为0。显然这两组边基函数都不是正交的。然后对已有的基函数进行扩充,可构造出新的一组正交基函数。首先我们要定义新的内积形式:〈u,v〉=∑i=13αiu(mi)⋅v(mi)(21)这里mi代表第i号边的中点,ml代表编号为i的面的中心点。对于二维的正交基函数,要再定义3个辅助的基函数:Bi=ξjξkn^ii=1,2,3(22)ξi是第i号边的单位法向矢量。之所以构造这3个辅助基函数是因为它们具有以下性质:1)沿剖分的棱边的切向分量为零;2)Bi仅在第i号边上有非零法向分量;3)根据(21)式的内积定义,这些基函数是两两正交的。那么我们可以定义新的一组基函数Z:Ζi=Wi-∑j=13〈Wi,Bj〉〈Bi,Bj〉Bj.(23)由它的定义可知,新的边基函数Zi是两两正交的。同样的思路,文献中构造了一组三维正交完备的基函数。2.2.2子域和局部时间步另一种有效解决TDFEM庞大计算域的方法是子域分解技术,即将原始的计算域分成几个较小的子域并采用并行计算。每个子域问题包含了缩小规模的计算域,可以用稀疏矩阵求解器直接求解,并且与原来的单一区域问题相比,它减少了整体的计算复杂性和计算时间。文献中提出了一种新型的子域分解方案。以一个一般的计算域V为例,以金属表面S1和阻抗表面S2为边界。通过任意的人工边界SD将V分裂成V1和V2两个子域。子域上的电场和磁场值用Ei和Hi表示,i=1,2。在子域上二阶矢量波方程可以写成如下的形式:∇×(μr-1⋅∇×Ei)+1c02εr∂2Ei∂t2+μ0σ∂Ei∂t=-μ0∂Ji∂t(24)∇×(εr-1⋅∇×Ei)+1c02μr∂2Ηi∂t2+μεrσ∂Ηi∂t=∇×(εr-1Ji)(25)S1上的边界条件是n^×E=0,n^×(∇×Η)=0(26)S2上的边界条件是n^×(μ-1∇×E)+γe∂∂t(n^×n^×E)=0n^×(ε∇×Η)+γh∂∂t(n^×n^×Η)=0(27)如果取γe=ε/μ和γh=μ/ε就是ABCs,对(24)(25)式以矢量基函数Ni测试得到弱形式的波方程,强加边条件并离散化,就可以得到矩阵形式的方程并求解。该方法是基于每个子域的双场二阶矢量波方程,对同一空间网格却不同的时间域上求解电场和磁场。每个子域问题可以小到足以靠一个稀疏矩阵求解器直接求解。这样,通过对每个子域的预先分解,并在时间步进之前将信息存储在内存。在每个时间步,子域的问题靠高效利用局部的预先分解矩阵解决。它的独特之处在于并不需要解全局的矩阵方程,缺陷是由于该方法将原来计算域分解成了几个单独的子域,而各个子域问题由于离散误差,将不再与原有的问题精确匹配,导致数值色散误差的增加。3正演中间层存在的问题如上节所述,时域有限元法在PML的应用、正交基函数的引进、子域分解技术等方面取得了一些进展,但尚有若干问题需要解决,还存在很大的发展空间,例如:3.1混合方法3.1.1求解各区域的边界积分文献提出了一种新型混合时域有限元-边界积分方法(TDFE-BI)分析三维电磁散射问题。它使用了一种人工边界将无限的解区域划分成内部和外部区域,分别用时域有限元(TDFEs)和边界积分(BIs)求解,两者通过在人工边界上的场连续性连接起来。它的优点是可以使用一种快速算法,即多层平面波时域(PWTD)方法来求解边界积分。在考虑一个电大尺寸目标时,引用这个方法大大降低了计算费用;缺点是每个时间步上需要求解一个矩阵方程。有限元-边界积分混合方法(FEM-BI)在频域有着较好的应用,在时域则刚开始起步,尚待发展。3.1.2tdfem-floquet定理混合算法周期性结构在各种电磁散射和辐射问题的应用有着悠久的历史。特别是新型的、人工构造的周期结构——电磁带隙材料和左手材料(其介电常数和磁导率均为负值)的提出,掀起了一股研究周期性结构的热潮。为了分析这类周期结构的电磁散射和辐射特性,文献、提出了在二维和三维问题上实现TDFEM与FloquetABC混合的方法:当一个周期性结构包含许多个周期单元晶格时,一般认为对该结构在周期重复的方向上无限延展的近似分析是合理的。TDFEM-Floquet定理混合算法通过引入一个转换场变量使得在单元晶格的四个面上的周期性边条件很容易实现,能从对一个周期单元晶格的分析得到无限周期结构的特性,该方法的优点是将有限元计算区域限制在一个周期的范围,大大缩小了计算量。缺点是在每个时间步上也需要解一个矩阵方程,而且转换场变量的引入对时域稳定性的影响还没有数学上严格的证明。3.1.3新型混合方法PSTD法借助平面波谱展开理论,利用快速傅里叶变换(FFT)及其逆变换代替麦克斯韦方程中的空间差商,该法最初应用于FDTD方法,在FDTD方法中PSTD法在处理复杂的曲线边界的时候采用阶梯近似,有较大的精度误差。时域有限元法(TDFEM)使用于处理复杂的曲线边界,但是,它需要更多的内存。为了收益这两种方法的好处,文献提出了一种新型时域拟谱法(PSTD)和时域有限元法(TDFEM)混合方法。相比TDFEM-FDTD方法,TDFEM-PSTD大大缓和了空间离散上的条件限制。3.2层级单元的求解除了混合算法的研究,我们也可以通过自适应FEM网格剖分技术改进TDFEM的精度。通常的改善有限元分析精度的方法是利用高阶基函数,这些基函数可分为两类:插值基函数和层级基函数。层级单元可以通过两种方式来提高精度:一种方式是在误差估计大的区域中提高基函数的阶(p-细化);另一种方式是保持基函数的阶数不变,来细分单元(h-细化)。hp自适应有限元法(hp-FEM),顾名思义,是允许h和p都能在网格剖分上独立地改变。这就使得用更少的未知量获得想要的精度成为可能。然而,一个主要的难题是如何决定h和p的分布,以及如何确定“细化的方向”。另外网格的局部加密从原理上来说是很简单的事情,但是在三维的网格