孔才智+郝毫保
摘 要： 针对于大尺寸电磁辐射问题，将撕裂互连法应用于三维电大尺寸辐射问题的计算仿真。该算法是区域分解方法中的一种，求解区域划分成互不重叠的子区域，各子区域之间通过拉格朗日乘子将交界面的连续边界条件或者传输条件耦合。鉴于原来的传输条件存在不收敛或者收敛较慢情况，提出同时能在TE/TM凋落模式快速收敛的二阶传输边界条件作为子区域之间交界面的边界条件，并且引入虚拟激励流为交换信息。数值结果表明，该改进算法有效地提高了撕裂互连法在迭代求解中的收敛性。在求解三维电大辐射问题时，该算法仿真结果与有限元算法结果一致，表明撕裂互连算法是计算大尺寸电磁辐射问题的一种有效方法。
关键词： 撕裂互连法； 有限元法； 二阶传输边界条件； 电磁辐射
中图分类号： TN911?34； O441 文献标识码： A 文章编号： 1004?373X（2015）16?0004?05
Application of tearing and interconnecting method in large?scale electromagnetic radiation under second?order transmission conditions problem
KONG Caizhi， HAO Haobao
（South China Normal University， Guangzhou 510006， China）
Abstract： The tearing and interconnecting method is applied to the simulation of three?dimensional （3?D） large?scale electromagnetic radiation to solve the problems of large?scale electromagnetic radiation. As one of the domain decomposition methods， this method reduces the computational complexity by dividing the domain under solution into several non?overlapping subdomains. The continuous boundary conditions or transmission conditions of the sub?domain interfaces are coupled by using the Lagrange multipliers and an iterative process is used to realize the solution of the original global problem. For further increasing the convergence of this method， a second?order transmission condition that can transmit both TE and TM evanescent modes is developed A virtual excitation source is introduced as the exchange information. Numerical simulation results reveal that the modification has improved effectively convergence performance in the iteration solution of the tearing and interconnecting method. Regarding to the electromagnetic radiation problems， the results are well matched with those obtained by finite element method.
Keywords： tearing and interconnecting method； finite element method； second?order transmission condition； electromagnetic radiation
0 引 言
传统的有限元方法在处理电大尺寸辐射问题和复杂细节结构目标上出现未知数的数量很多，导致全局线性系统矩阵规模庞大，矩阵形态较差，迭代求解往往很慢甚至不收敛。区域分解法提出在并行计算上的优势，是有限元方法解决这些问题的重要方法之一。在众多区域分解算法当中，撕裂互连法被证明是非常有前景和高效的算法之一[1?3]。该算法将原求解区域划分成若干个不重叠的子区域，在子区域的交界面上使用拉格朗日乘子隐式强加场的连续性条件把子区域连接起来，通过消去法将全局未知量的矩阵变成子区域交界面未知量的系统矩阵，一旦获得交界面上的未知量，代回到相应的子区域，各子区域就可以单独求解。撕裂互连算法不仅具有区域分解法的“分而治之”的策略优点，而且减少全局未知量为交界面的未知量，使其系统矩阵应用Krylov子空间迭代求解[4?5]，避免了求解大矩阵的复杂问题。为了提高算法在高频电磁中的迭代收敛性，实现子区域交界面的一阶Robin类型传输条件（FOTC）替代连续性条件耦合，并将该算法应用到三维电大尺寸电磁问题上，模型用一阶吸收边界条件或者完全匹配层（PML）截断求解区域[6?8]。由于一阶传输条件（FOTC）只能用于传输模式，为了能够传输TE的凋落模式和提高迭代算法的收敛性，在子区域之间的交界面引入二阶传输条件（SOTC?TE）求解大规模多区域的电磁模型[9?10]。具体的计算实例表明，二阶传输条件具有更高的收敛效果。
为了达到提高算法精度和分析复杂结构的目标，网格剖分变密，导致分界面的单元中存在凋落模式，凋落模式会严重影响到传输条件的收敛模式，上述提到的二阶传输条件（SOTC?TE）仅能在TE模式下快速收敛。本文将推导出适用于撕裂互连法新型的SOTC?FULL二阶传输条件，使得该传输条件具备更明显的物理意义，该传输条件能同时在TE和TM凋落模式下快速收敛。将改进后的算法应用于求解三维电磁辐射问题中，采用一阶吸收边界条件截断求解区域。
1 撕裂互连原理与传输条件
一个有截断边界的有源电磁问题可以归结为在[Ω]空间上求解边值问题的矢量波动方程。不失一般性，将求解自由空间划分为互不重叠的子空间。每个子空间内建立有限元矢量波动方程以及限定的边界条件：
在[Ωs]子空间内满足：
[?×μ-1r?×E-k2oεrE=-jkozoJimp] （1）
在PEC边界满足的边界条件：
[n×E=0] （2）
在PMC边界满足的边界条件：
[n×?×E=0] （3）
一阶吸收边界条件（ABC）：
[n×（?×E）+jkin×（n×E）=Uinc] （4）
其中：[Ωs]表示第[s]子空间；[ko]和[zo]分别为自由空间的波数和波阻抗；[Jimp]为内部强加源；PEC为理想导体边界条件；PMC为理想磁场边界条件；[Uinc=][n×（?×Einc）+jkin×（n×Einc）]为入射波，[ki=koμrεr]。
当一个原求解区域划分为若干个子区域，只要保证子区域之间场的连续性，则由子区域与子区域之间的传输条件所构成的问题与原问题等效，这是区域分解法的基本思想。撕裂互连法是将拉格朗日乘子引入到传输条件中，最初的算法是用Dirichlet边界连续性条件表示传输条件，由于在该条件下存在无法收敛或者收敛很慢的情况，因此引用Robin一阶传输条件；进一步应用发现，当网格划分变密，频率变低时[11?12]，分界面单元中将存在凋落模式，该模式影响传输条件的收敛效果。为了改善收敛效果，本文引入非共形区域分解法（Non?conformal Domain Decomposition Method）[13?15]中能同时在TE和TM凋落模式下快速收敛的二阶传输条件（SOTC?FULL），推导出使该方法适用于撕裂互连算法的形式，表示如下：
[ni×（μ-1ri?×Ei）+αini×（ni×Ei）-βi?× [ni（?×Ei）n]+γi?τ?τ·[ni×（μ-1ri?×Ei）]=Λ] （5）
式中：[i]表示子区域之间的界面；[αi=jkoμsqrεsqr]，[μsqr]和[εsqr]分别表示子区域界面两边的相对介电常数和磁导率的均值；[βi]和[γi]由最小网格和基函数的阶数决定，[βi=-j（k+ktez）]，[γi=1（k2+kktmz）]，在这里取[kmax，ter=πhmin]，[kmax，tmr=π/（2×hmin）]，[hmin]是子区域界面的最小网格尺寸，若[βi=0]，[γi=0]，则式（5）为一阶传输条件FOTC，若[γi=0]，式（5）变为SOTC?TE传输条件，[Λ]为子区域界面引入的拉格朗日乘子未知量。
基于上述边界条件的基础，本文推导出一种新型的二阶传输条件，对于式（5）中的第4项[γi?τ?τ·[ni×（μ-1ri?×Ei）]]，因为[ni×（μ-1ri?×Ei）]只有切向分量，所以：
[?τ·[n×（μ-1r?×E）]=（?-?n）·[n×（μ-1r?×E）] =[?·[n×（μr-1?×E）]]]
根据矢量恒等式和在自由空间的无源二阶矢量波动方程得到：
[?·[n×（μ-1?×E）]=n·[?×（μ-1?×E）]=k2oεrn·E]
将式（5）中的传输条件展开分别表示为交界面2个不同子区域的传输条件，把分界面上的传输条件看成是一种激励，本文中以[gi]来表示，这一激励类似于面电流。因此式（5）可以分解为以下2种形式：
[ni1×（μ-1ri1?×Ei1）+αi1ni1×（ni1×Ei1）-βi1?× [ni1（?×Ei1）n]+γi1?τn·E=gi] （6）
[-ni2×（μ-1ri2?×Ei2）+αi2ni2×（ni2×Ei2）-βi2?× [ni2（?×Ei2）n]-γi2?τn·E=gi] （7）
式中：[i1]表示[Si∈Ω1]；[i2]表示[Si∈Ω2]。
对于子区域[Ωs]联立方程式（1）～式（4），式（6），式（7），将子区域矢量场用矢量棱边元展开，并运用伽辽金方法，获得子区域的矩阵方程：
[Ks{Εs}={fs}-i∈neibourΩssiNs·gids] （8）
式中：s表示子区域编号；[Εs]为待求的电场未知量；子区域有限元矩阵[Ks]和激励向量[fs]的具体形式如下：
[Ks=Ω[μ-1r（?×Ns）·（?×Ns）-k2oεrNs·Ns]dv+ jkoSABC（n×Ns）·（n×Ns）ds+αs（n×Ns）·（n×Ns）ds+ βs（?×Ns）n（?×Ns）nds+γs（n·Ns）（?τ·Ns）ds] （9）
[{fs}=-jkozoΩNs·Jimpdv-SABCNs·Uincds] （10）
每个子区域中的单元边未知量按照位置可以分为3种类型[E=[Ev，Ei，Ec]]，分别表示内部棱边，交界面和拐角处未知量，其中[Ev，Ei]可以归为一类[Er]。因此子区域的矩阵有限元系统可以表示为：
[Ks=KsrrKsrcKscrKscc， {Εs}=ΕsrΕsc] （11）
交界面[Si]上定义与交界面相关的矢量[ρi]以及布尔矩阵[Bi]。假定某一棱边在子区域[Ωi]和子区域[Ωj]的交界面[Si]上的编号[p]，在子区域内[Ωi]的编号为[q]，在子区域[Ωj]的编号为[t]，则[BTi（p，q）=1]，[BTj（p，t）=1]，因此：
[BTiρi=-BTjρi=siNs·gids] （12）
将所有的拐角（多个子区域重叠的边）处的棱边进行编号，并定义布尔矩阵[Bsc]为棱边对应在拐角上的编号与子区域内部编号的映射关系矩阵。假定在子区域内[Ωs]内某一拐角的编号为[m]，在全局拐角棱边编号为[n]，则[Bsc（m，n）=1]，用[Ec]表示全局拐角未知量：
[BscEc=Esc] （13）
联立式（8），式（11）～式（13），在子区域[Ωs]内矩阵形式可表示为：
[KsrrKsrcKscrKsccΕsrΕsc=fsrfsc--Biρi+Bjρi 0] （14）
不失一般性，将整个求解区域划分为4个子区域，如图1所示。
图1 区域分解图
在子区域之间交界面[S1]，子区域[Ω1]和[Ω2]之间的界面，同时显示强加Dirichlet连续性边界条件[B1E1+B2E2=0]。整合4个子区域系统矩阵方程式（14），得到全局系统矩阵[K]如下：
由式（16）可以看出，本文推导出一种新型的二阶传输条件以及引入的虚拟激励流后，获得的系统矩阵与原有撕裂互连算法中的系统矩阵一致，很好地保证了撕裂互连算法区别于其他区域分解法最大特点[5]，因此可以应用撕裂互连算法中的求解方法求解系统线性矩阵方程。应用高斯消去法，消去全局拐角未知量[Ec]，可以获得关于交界面上的对偶未知量[λ]的全局系统方程为：
[[Krr+KrcK-1ccKTrc]λ=dr-KrcK-1ccfc] （17）
[{Ec}=[Kcc]-1（{fc}+[Krc]T{λ}）] （18）
[{Esr}=[Ksrr]-1（{fsr}-[Ps]T{λ}-[Ksrc][Bsc]{Ec}）] （19）
式中[Krr]，[Krc]，[Kcc]，[Krc]，[dr]，[fc]的具体表达式以及以上公式的详细推导过程见文献[16]。式（17）称为全局界面系统方程，未知量[λ]只定义在区域之间的界面。经过验证，该方程求解通过Krylov子空间迭代方法如稳态双共轭梯度算法（BiGCSTAB）和广义最小余量算法（GMRES）求解更有效。获得交界面的未知量[λ]后，各个子域内部的电场值可以通过子域方程式（18），式（19）分别求出。原来三维有限元方程问题转化为二维界面问题，交界面上的未知量通常比较少，因此求解该方程的维度和计算量比原有限元系统方程减少很多，所以该方法能够求解电大电磁辐射和散射问题。
2 数值算例
算例1：为了验证撕裂互连法在不同的传输条件FOTC，SOTC?TE，SOTC?FULL下的收敛性和算法的正确性，仿真了自由空间波传播的算例，如图2所示。
图2 自由空间波传播求解区域
该模型内外都是空气，入射波[E=xe-jkz+ye-jkz]，入射波沿Z轴传播的与X轴成45°极化方向的平面波。求解空间用一阶吸收边界条件截断外边界，整体求解区域剖分为4个子区域，网格划分密度为[λ20]，入射波频率为10 GHz，收敛精度为[10-8]，采用完全相同的网格剖分和子区域划分方式，全局具有2 240个未知量。从图3迭代结果可以看出SOTC?FULL的收敛速度比FOTC，SOTC?TE传输条件的收敛效果好。
图3 自由空间波传播模型不同传输条件迭代结果对比图
算例2：因为撕裂互连法适用于大尺寸电磁分析特别是天线有限阵列，因为在多数典型的天线阵列中，很多天线单元是相同的。对天线阵列进行仿真，只要构造和计算几个子区域的有限元矩阵。整个求解区域划分为9个子区域，如图4所示。
图4 具有不同有限元矩阵的9个子区域图
由于应用外边界如吸收边界条件和完全匹配层，边和角的有限元矩阵不同于内部的有限元矩阵，为此，其他的有限元矩阵可由这9个矩阵确定下来。也就是说，这9个子区域的有限元矩阵[Ksrr]一旦构成和计算出来，这些矩阵可以在全局系统界面方程迭代过程和最后的各个子区域求解过程中替代其他子区域矩阵。
本文应用撕裂对接算法仿真和分析微带馈电的贴片天线，如图5所示。
图5 微带馈电贴片天线结构图
16 mm
图6给出了撕裂互连算法求解该模型的S11参数结果图，文章中的算法计算结果和ANSYS求解的结果吻合的非常好。
图6 贴片天线的S11参数算法仿真与ANSYS仿真结果对比图
3 结 语
本文推导和提出适用于撕裂互连算法的二阶传输条件（SOTC?FULL），并简要介绍了其公式推导过程。验证了该二阶传输条件（SOTC?FULL）比一阶传输条件（FOTC）有更好的收敛性，并且将其应用到撕裂互连算法求解电磁辐射问题上分析微带馈电贴片天线。仿真结果与ANSYS软件中的求解结果符合程度高，从而证明了撕裂互连算法是计算大尺寸电磁辐射问题的一种有效方法。
参考文献
[1] LI Y， JIN J M. A vector dual?primal finite element tearing and interconnecting method for solving 3?D large?scale electromagnetic problems [J]. IEEE Transactions on Antennas and Propagation， 2006， 54（10）： 3000?3009.
[2] LI Y J， JIN J M. A new dual?primal domain decomposition approach for finite element simulation of 3?D large?scale electromagnetic problems [J]. IEEE Transactions on Antennas and Propagation， 2007， 55（10）： 2803?2810.
[3] 杜磊.时域有限元电磁计算方法的研究[D].南京：南京理工大学，2010.
[4] ZHAO Kezhong. A domain decomposition method for solving electrically large electromagnetic problems [D]. Ohio State： The Ohio State University， 2007.
[5] JIN J M， RILEY D J. Finite element analysis of antennas and arrays[M]. John Wiley & Sons， 2009. 336?387.
[6] 贾会亮.频域有限元区域分解法的研究[D].南京：南京理工大学，2009.
[7] XUE M F， JIN J M. Nonconformal FETI?DP methods for large?scale electromagnetic simulation [J]. IEEE Transactions on Antennas and Propagation， 2012， 60（9）： 4291?4305.
[8] 高红伟，盛新庆.三维非均匀电大目标散射的有限元撕接区域分解法计算[J].北京理工大学学报，2012，32（6）：602?606.
[9] PENG Zhen， RAWAT Vineet， LEE Jin?Fa. One way domain decomposition method with second order transmission conditions for solving electromagnetic wave problesms [J]. Journal of Computational Physics， 2010， 229： 1181?1197.
[10] XUE M F， JIN J M. A hybrid conformal/nonconformal domain decomposition method for multi?region electromagnetic modeling[J]. IEEE Transactions on Antennas and Propagation， 2014， 62： 2009?2021.
[11] 马金.非均匀介质中低频近场的分析与探测应用[D].成都：电子科技大学，2013.
[12] YAO W， JIN J M， KREIN P T. Application of the LU recombination method to the FETI?DP method for solving low?frequency multiscale electromagnetic problems [J]. IEEE Transactions on Magnetics， 2013， 49（10）： 5346?5355.
[13] PENG Z， LEE J F. Non?conformal domain decomposition method with mixed true second order transmission condition for solving large finite antenna arrays [J]. IEEE Transactions on Antennas and Propagation， 2011， 59（5）： 1638?1651.
[14] PENG Z， LEE J F. A scalable nonoverlapping and nonconformal domain decomposition method for solving time?harmonic maxwell equations in R^3 [J]. SIAM Journal on Scientific Computing， 2012， 34（3）： A1266?A1295.
[15] MA J， NIE Z， SUN X. Efficient modeling of large?scale electromagnetic well?logging problems using an improved nonconformal FEM?DDM [J]. IEEE Transactions on Geoscience and Remote Sensing， 2014， 52（3）： 1825?1833.
[16] YAO W. Accurate， efficient， and stable domain decomposition methods for analysis of electromechanical problems [D]. Urbana： University of Illinois at Urbana?Champaign， 2014.
