麦克斯韦方程时域无网格算法及其应用研究

高煜堃 陈红全 蒲赛虎

（南京航空航天大学航空宇航学院，江苏 南京210016）

引 言

求解麦克斯韦方程的时域数值方法主要有时域有限体积法（Finite-Volume Time-Domain，FVTD）［1］、时域有限差分法（Finite-Difference Time-Domain，FDTD）［2］、时域有限元法（Finite-Element Time-Domain，FETD）［3］以及时域多分辨率分析法（Multi-Resolution Time-Domain，MRTD）［4］等，大多基于网格进行计算.其中FVTD方法的提出是借鉴了计算流体力学（Computational Fluid Dynamics，CFD）中经典的有限体积做法，该方法计算区域的空间离散基于网格填充.按照所基于的网格，FVTD方法大体上可分为结构网格算法［1］和非结构网格算法［5］.结构网格算法网格划分简单，数据结构清晰，但一般只适合处理相对简单的外形，对于复杂的几何外形一般需要用到块结构等特殊处理技术.非结构网格算法适合对复杂外形进行有效的描述和离散，但网格生成相对复杂和耗时.

可以注意到，上述算法从基于结构网格发展到基于非结构网格，其主要特征是打破了网格结构化的约束.也有学者考虑进一步打破网格化约束，尝试利用无网格技术进行电磁场数值模拟的研究，发展出不同的电磁场无网格求解方法［6-11］，其中无单元Galerkin方法［6］和利用径向基函数的无网格法［7-9］最有代表性.无单元Galerkin方法是以打破有限元单元约束为特征，其形函数的构造主要采用滑动最小二乘法，已在静电箱体［6］等静态电磁场问题的求解中取得了成功，目前，已从全局弱式无网格法发展为局部弱式Petrov-Galerkin无网格法［10］，并且，求解问题也从静态电磁场发展用于求解时变电磁场［11］.利用径向基函数的无网格法常见的有径向基点插值无网格法［7-8］，该法是为避免采用多项式基点插值所引起的奇异性问题而提出的［12］，其显式无网格法已成功用于求解三维电磁场问题［7］，目前，算法已从显式发展到无条件稳定的隐式无网格法，突破了显式方法时间步长的限制，提高了计算效率［8］.但也必须注意到，上述两类无网格方法与近似解涉及的形函数（或基函数）的选取相关，因求解一般会涉及到矩阵的求逆运算，形函数的选择会影响到矩阵的质量［9］，也会影响到边界条件的施加，有时边界条件需要采用特殊的方法强制给定［12］.

我们注意到，在CFD领域，近年来也出现了一种新的无网格算法.空间导数的近似及通量运算等方面与上述两类无网格方法不同，无网格点上的空间导数可基于局部点云结构通过极小曲面逼近计算得到.空间离散涉及的通量运算采用CFD中的近似黎曼解处理.由于无网格算法计算域的离散只涉及布点填充，既可以采用网格点来布点，也可以打破传统的网格算法在拓扑结构上的约束，而根据需要直接布点，因此，在对复杂外形的处理上更加灵活［13］.目前，许多研究者将无网格方法用于求解Euler方程，并成功应用于复杂外形的流场计算［13-14］.

鉴于所研究的Maxwell方程和流体力学中的Euler方程具有双曲型等相同的数学特性，本文借鉴CFD中的无网格方法来发展用于求解麦克斯韦方程的时域无网格算法.为此，针对本文无网格算法涉及的计算区域布点问题，提出了一种基于直角网格点的直接布点方法.采用该方法，可以对计算区域进行快速布点填充，同时得到计算所需的点云结构.在当地点云上，引入二次极小曲面逼近确定计算点上的空间导数.空间离散涉及的通量运算借鉴近似黎曼解方法构造确定，给出了具体的基于点云结构的通量运算处理方法.时间离散则采用四步Runge-Kutta法推进求解.最后，结合算例，给出了圆柱和NACA0012翼型典型算例的验证与分析，并展示了算法处理多体及多部件电磁干扰问题的能力.

1 时域麦克斯韦方程

对于横磁波（Transverse Magnetic Wave，TM）在直角坐标系中，守恒型无量纲时域麦克斯韦方程可写为

2 计算区域布点离散及点云生成

如前所述，无网格算法计算区域的空间离散只涉及布点填充.一般来说，既可以方便地选取结构网格或非结构网格点，也可以打破传统的网格方法在拓扑结构上的约束，按照一定的要求直接布点.鉴于无网格算法布点灵活的特点，这里给出一种基于直角网格点的直接布点离散方法.

采用直角网格点离散计算区域具有布点速度快、质量好等特点，因此大部分计算区域尽可能用直角网格点进行布点离散；但是完全采用直角网格点在物面处精度会受到影响，为此，本文在物面附近局部采用沿物面法向直接布点的做法（见图1）.这样做不仅能够弥补直角网格点对物体外形适应性差的弱点，而且这种布点方法与物面边界和远场边界的正交性均可控，便于边界条件的施加.例如，在数值计算中，通常在物面处引入入射波分量，并根据物面边界条件计算出对应的散射波分量，然后根据控制方程模拟散射波的传播，因此物面处正交性的可控在一定程度上可以降低因边界条件的处理而引起的计算误差.

图1 无网格算法布点填充示意图

计算区域布点填充后，在点的局部需构成无网格算法实施所要求的点云结构.以点云Ci为例（见图2），中心点Vi自然纳入点云Ci中，下面需要合理选取一定数目的卫星点.本文的做法是以中心点Vi为圆心，ri为半径画圆，将圆内除Vi以外的点（V1、V2、V3、V4、V5和V6）视为卫星点一并纳入点云Ci中.其中圆半径ri与点云当地点与点之间的距离相关，通常可通过扩大或缩小ri来选取合适数目的卫星点.

图2 无网格点云Ci示意图

3 时域麦克斯韦方程的无网格离散求解

3.1 空间导数的确定

基于上述点云结构，在点云Ci上，假定中心点Vi附近的函数值分布满足函数f＝f（x，y），则f可用Vi处的值fi＝f（xi，yi）通过泰勒级数展开逼近：

式中：h＝x-xi；l＝y-yi；ai（i＝1，2，…，5）即为函数f在中心点Vi处相应的各阶偏导数.

采用时域无网格方法，卫星点Vk（k＝1，…，M）处的函数值为已知，记为fk，则这M个卫星点函数逼近的总体误差可表示为

基于总体误差G极小，可以得到空间导数满足

逼近函数可整理为

式中αk和βk仅与离散点坐标相关，在时间推进计算前可以一次求出.由式（6）可以确定函数f在Vi处的空间导数为

空间导数确定以后，即可代入通量计算式进行通量运算.有关本文时域无网格算法的基于点云结构的通量运算处理方法将在下一节介绍.

3.2 通量运算

在点云Ci上，应用式（7），中心点Vi处的通量相关项可近似写为

观察式（8）中等号的右端项，由于中心点Vi处的函数值为已知，系数αik和βik根据式（5）已经计算得到，因此求和项∑（αikF1i＋βikF2i）可以直接计算.下面介绍求和项∑（αikF1ik＋βikF2ik）的具体处理方法.在中心点Vi与每一个卫星点Vk连线的中点处，定义一个数值通量为

采用近似黎曼解确定中点处的数值通量Qik，即：

rik是从点i指向点k的矢量，注意这里U表示上述矢量中的任一分量，▽Ui和▽Uk为计算点上物理量的梯度，由式（7）计算给出.

3.3 时间离散

采用无网格方法对麦克斯韦方程进行空间离散以后，在点云Ci上，可以得到其半离散形式为

式中，Ri为计算点上的残差.对于式（14），按照四步Runge-Kutta方法进行时间推进求解［1］.

数值求解时还需补充计算涉及的边界条件，即物面边界条件和截断边界条件.本文物面应用良导体边界条件［1］，即总电场的切向量和总磁场的法向量在良导体表面为零；计算域的截断位置则采用完全匹配层（Perfectly Matched Layer，PML）边界条件，该条件中参数的取值详见文献［5］.

4 算例与分析

本文已按上述方法研制了对应的计算程序，并结合算例对方法进行了分析与考核.算例所涉及的计算域中的x和y坐标均为以入射波波长λ为特征长度、无量纲化后的空间位置.如图3所示，TM平面波沿k方向照射目标，定义k轴与x轴之间的夹角为入射角φ，那么归一化后的入射波分量可表示为

式中k＝xcosφ＋ysinφ.

4.1 典型算例的考核计算

这里先选用有解析级数解［15］供比较的金属圆柱算例进行算法考核.计算时采用直接布点方法进行计算域的离散（如图1所示），总点数为20 911，圆柱半径设为λ，远场截断边界取在6λ处，TM波入射角取为φ＝0°.计算得到的圆柱双站雷达散射截面（Radar Cross Section，RCS）分布见图4.从图中不难看出，计算区域的离散不论是基于直接布点，还是基于网格点（总点数为18 741），计算得到的双站RCS的峰值及其在整个双站角范围内的分布都能与级数解吻合.

图3 TM平面波

图4 圆柱的双站雷达散射截面分布

接下来给出NACA0012金属翼型算例.该算例源自文献［16］，常用于电磁散射算法的验证，其翼型弦长取为10λ.计算时将NACA0012翼型置于20λ×20λ的计算域中间，计算域离散采用直接布点和取用非结构网格点两种布点方法，对应的点数分别为33 829和36 889.该算例10λ弦长对应高频段，较半径为λ的圆柱算例需更密的点云，为了显示清楚，这里只给出以稀点云分布的示意图（如图5所示），TM波入射角取为φ＝90°.

图6给出了计算得到的NACA0012翼型双站RCS分布.从图中不难看出，采用本文方法计算得到的双站RCS的峰值及其在整个双站角范围内的分布与文献［16］三阶结构网格FVTD的计算结果一致.

4.2 主体带外挂的多体干扰电磁散射计算

图5 NACA0012翼型的点云分布示意图

图6 NACA0012翼型的双站雷达散射截面分布

为了考察本文发展的时域无网格算法求解多体干扰电磁散射问题的能力，对飞行器主体带外挂的情形进行数值模拟.沿用文献［17］的做法，选用NACA0012翼型模拟飞行器主体，选用缩小的翼型模拟外挂.图7为主体带外挂的计算点云分布的局部放大图，计算域为12λ×16λ，总点数为25 520.

为了与文献［17］结果进行比较，大翼型弦长取为4λ，小翼型弦长取为2λ，TM波入射角取为φ＝0°.计算得到的散射场等值线分布见图8.从图中可以看出：外挂对主体的电磁散射特性影响显著，主体和外挂之间的电磁散射干扰清晰可见，大体上与文献［17］的结果一致.

图9为计算得到的双站RCS分布图，从图中可以看到本文计算得到的双站RCS的峰值及其在整个双站角范围内的分布与文献［17］精确控制法结果接近.

4.3 飞机模型的电磁散射场数值模拟

图7 计算点云分布的局部放大图

图8 主体带外挂的散射场分布

图9 主体带外挂的双站雷达散射截面分布

该算例考察本文方法处理多部件干扰电磁散射问题的能力.多部件飞机模型的几何外形源自文献［18］，机体长取为8λ.将飞机模型置于20λ×20λ的计算域中，计算点云分布的局部放大图如图10（a）所示，总点数为25 889，TM波入射角分别取为φ＝0°、φ＝90°和φ＝180°，计算得到的散射场等值线分布分别对应图10（b）、（c）和（d）.从图中可以看到飞机的各个部件之间存在电磁散射干扰：当φ＝0°时，对应电磁波从飞机尾部照射机体，电磁散射干扰作用最强的区域发生在垂尾和机身之间，强散射方向出现在双站角0°和180°附近（见图10（b））；当φ＝90°时，对应电磁波从下方照射飞机，机头与机身之间存在明显的电磁散射干扰作用，在双站角90°方向和270°附近散射较强（见图10（c））；当φ＝180°时，对应电磁波从机头照射飞机，机头与机身、机身与垂尾之间均存在散射波的叠加作用，最强散射只出现在双站角180°方向（见图10（d））.该飞机模型在不同方向电磁波照射时的上述散射特征，大体上与文献［18］的虚拟区域法结果一致（图中没有画出），同时，在一定程度上也反映出本文发展的时域无网格方法适合处理多部件干扰等复杂外形.

图10 飞机模型的点云分布和散射场等值线分布图

5 结 论

借鉴CFD中求解Euler方程的无网格方法，成功地发展出求解麦克斯韦方程的时域无网格算法.该算法可快速布点离散计算区域，基于点云结构借鉴近似黎曼解方法处理通量运算，其空间导数由二次极小曲面逼近确定.算例验证分析表明：不论基于网格点，还是直接布点，采用本文方法得到的计算结果都能与级数解或文献结果吻合；得到的飞行器主体带外挂以及飞机模型的电磁散射场，一定程度上展示出本文所发展的算法具有处理多体及多部件干扰等复杂问题的能力.

［1］许 勇，乐嘉陵.基于CFD的电磁散射数值模拟［J］.空气动力学学报，2004，22（2）：185-189.XU Yong，LE Jialing.CFD-based numerical simulation of electromagnetic scattering［J］.Acta Aerodynamica Sinica，2004，22（2）：185-189.（in Chinese）

［2］胡成林，窦文斌.新型亚毫米波扫描天线的FDTD分析［J］.电波科学学报，2010，25（4）：725-727＋814.HU Chenglin，DOU Wenbin.FDTD analysis of a novel sub-millimeter scanning antenna structure［J］.Chinese Journal of Radio Science，2010，25（4）：725-727＋814.（in Chinese）

［3］吴 霞，周乐柱.时域有限元法在计算电磁问题上的应用及发展［J］.电波科学学报，2008，23（6）：1209-1216.WU Xia，ZHOU Lezhu.Application and development of time-domain finite element method on EM analysis［J］.Chinese Journal of Radio Science，2008，23（6）：1209-1216.（in Chinese）

［4］李炎红，杨 峰，聂在平，等.基于MRTD的微带天线的辐射分析［J］.电波科学学报，2010，25（4）：674-678.LI Yanhong，YANG Feng，NIE Zaiping，et al.Microstrip antenna radiation based on MRTD method［J］.Acta Electronica Sinica，2010，25（4）：674-678.（in Chinese）

［5］SHI Yan，LIANG Changhong.The finite-volume timedomain algorithm using least square method in solving Maxwell’s equations［J］.Journal of Computational Physics，2007，226（2）：1444-1457.

［6］PARREIRA G F，SILVA E J，FONSECA A R，et al.The element-free Galerkin method in three-dimensional electromagnetic problems［J］.IEEE Transactions on Magnetics，2006，42（4）：711-714.

［7］YU Yiqiang，CHEN Zhizhang.A 3-D radial point interpolation method for meshless time-domain modeling［J］.IEEE Transactions on Microwave Theory and Techniques，2009，57（8）：2015-2020.

［8］CHEN Xiaojie，CHEN Zhizhang.YU Yiqiang，et al.An unconditionally stable radial point interpolation meshless method with laguerre polynomials［J］.IEEE Transactions on Antennas and Propagation，2011，59（10）：3756-3763.

［9］MIRZAVAND R，ABDIPOUR A，MORADI G，et al.Full-wave semiconductor devices simulation using meshless and finite-difference time-domain approaches［J］.IET Microwaves，Antennas ＆ Propagation，2011，5（6）：685-691.

［10］ZHAO Meiling，NIE Yufeng.A study of boundary conditions in the meshless local Petrov-Galerkin（MLPG）method for electromagnetic field computations［J］.CMES，2008，37（2）：97-112.

［11］SOARES D.Numerical modelling of electromagnetic wave propagation by meshless local Petrov-Galerkin formulations［J］.CMES，2009，50（2）：97-114.

［12］LIU G R，GU Y T.无网格法理论及程序设计［M］.王建明，周学军，译.济南：山东大学出版社，2008.

［13］陈红全.求解Euler方程的隐式无网格算法［J］.计算物理，2003，20（1）：9-13.CHEN Hongquan.Implicit gridless method for Euler equations［J］.Chinese Journal of Computational Physics，2003，20（1）：9-13.（in Chinese）

［14］MA Zhihua，CHEN Hongquan，WU Xiaojun.A gridless-finite volume hybrid algorithm for Euler equations［J］.Chinese Journal of Aeronautics，November 2006，19（4）：286-294.

［15］HALLEROD T and RYLANDER T.Electric and magnetic losses modeled by a stable hybrid with explicit-implicit time-stepping for Maxwell’s equations［J］.Journal of Computational Physics，2008，227：4499-4511.

［16］苏 敏，陈 刚，童创明.FVTD在电磁散射问题中的应用［J］.航天电子对抗，2008，24（3）：56-58.SU Min，CHEN Gang，TONG Chuangming.Application of FVTD method in CEM problems［J］.Aerospace Electronic Warfare，2008，24（3）：56-58.（in Chinese）

［17］陈红全，黄明恪.用精确控制法计算复杂区域的二维电磁散射场［J］.计算物理，2000，17（4）：414-420.CHEN Hongquan，HUANG Mingke.Computations of scattering waves in 2-D complicated fields by using exact controllability approach［J］.Chinese Journal of Computational Physics，2000，17（4）：414-420.（in Chinese）

［18］BESPALOV A.Application of fictitious domain method to the solution of the helmholtz equation in unbounded domain［R］.Le Chesnay：INRIA，1992.