求解三维麦克斯韦方程的时域无网格算法研究

高煜堃 陈红全

求解三维麦克斯韦方程的时域无网格算法研究

高煜堃 陈红全

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

发展了用于求解三维麦克斯韦方程的时域无网格算法.算法基于生成的无网格点云，通过泰勒级数展开结合加权最小二乘逼近计算点云中心点上的空间导数，并构造近似黎曼解处理空间离散涉及的通量运算；空间离散后的半离散方程则采用四步Runge－Kutta格式推进求解.结合求解三维麦克斯韦方程，给出了时域无网格算法的具体实施过程，并基于发展的算法，成功地模拟出金属球、立方体及进气道模型等三维散射目标的电磁散射场，获得的雷达散射截面能与理论解、矩量法或精确控制法等结果吻合.

时域无网格算法；麦克斯韦方程；点云；加权最小二乘；三维

引 言

无网格方法是在网格方法之后出现并发展的，其计算域的离散只涉及布点，既可以方便地选用已有的结构或非结构网格的格点来填充计算域，也可以打破传统网格方法［1－2］在拓扑结构上的约束，而根据需要直接进行布点离散，因此，无网格方法在处理复杂外形时更加灵活［3］.鉴于无网格方法具有布点灵活的特点，近年来，该法已被引入到计算电磁学（Computational Electromagnetics，CEM）领域，用于求解电磁场问题，其中代表性的方法有径向基点插值无网格法和局部弱式Petrov－Galerkin无网格法等［4－6］.这些方法在实施过程中通常会涉及到形函数（或基函数）的选取.形函数的选择不仅会影响到矩阵的质量［7］，还会影响到边界条件的施加，有时边界条件需采用特殊方法强制给定［8］.我们注意到在计算流体力学（Computational Fluid Dynamics，CFD）领域，近年来也出现了一种与CEM中做法不同的无网格方法.无网格点处的空间导数可基于方法实施中构造的局部点云结构通过极小曲面逼近得到，而空间离散涉及的通量运算则常采用近似黎曼解处理.目前，这种无网格方法已被成功应用于求解欧拉方程［9］，能够模拟出复杂外形的绕流［3］.

鉴于所研究的麦克斯韦方程和流体力学中的欧拉方程具有双曲型等相同的数学特性，可借鉴上述CFD中的无网格方法，发展出求解麦克斯韦方程的时域无网格算法.作者根据这一思想，已成功地发展出二维时域无网格算法［10］.本文意在进一步将二维算法拓展到处理面向实用的三维问题，发展出求解三维麦克斯韦方程的时域无网格算法.为此，基于无网格三维点云结构，改用适应性更好的函数加权最小二乘逼近计算离散点上的空间导数，空间离散涉及的通量运算则沿用文献［10］的做法，采用近似黎曼解方法确定，通量运算涉及的无网格点云中心点与卫星点连线中点处物理量的左右状态值由线性函数逼近给出，时间离散则按照四步Runge－Kutta法推进求解，形成三维时域无网格算法的具体实施方法.接着，本文先用金属球及立方体等典型三维散射体对发展的算法进行散射场及目标特性的考核运算，再给出进气道模型模拟的复杂散射体算例，以展示算法处理三维散射目标的能力.

1 三维时域无网格算法

1.1控制方程

在直角坐标系中，三维守恒型无量纲时域麦克斯韦方程可写为

式中：

Ez分别为电场强度E在x、y和z方向上的分量，Hx、Hy和Hz分别为磁场强度H在x、y和z方向上的分量，ε为介电常数，μ为磁导率，σm为磁阻率，σ为电导率.

1.2布点及点云生成简介

图1 点云Ci示意图

如前述，无网格方法计算区域的离散只涉及布点.一般来说，既可以方便地选取已有的结构网格或非结构网格的格点，也可以打破传统网格方法在拓扑结构上的约束，而根据需要直接布点.计算区域布点后，在点的局部需构成时域无网格算法实施所要求的点云结构.为了描述方便，先以二维情形为例，对无网格点云Ci的构造过程进行介绍.如图1所示，先将中心点i纳入Ci中，接下来再合理选取若干卫星点.本文通过以中心点i为圆心、ri为半径画圆，将圆内包含的所有点（点1、2、3、4、5和6，点i除外）视为卫星点一起纳入Ci.Ci中卫星点的数目可通过ri来控制.对于本文的三维情形，点云Ci的构造方法可描述为：以中心点i为球心、ri为球半径画球面，将球面内包含的所有点（点i除外）视为卫星点同中心点i一起纳入点云Ci中.

1.3空间导数的计算

基于上述点云Ci，若中心点i附近的函数值分布满足函数f＝f（x，y，z），那么f可用点i处的函数值fi＝f（xi，yi，zi）通过泰勒级数展开逼近

式中：h＝x－xi；l＝y－yi；m＝z－zi；ai（i＝1，2，3）为函数f在点i处的偏导数.

采用线性逼近，式（2）中函数f的近似值¯f可写为

逼近时，卫星点k（k＝1，…，M）处的函数值为当前值，可认为是已知的，记作fk，那么这M个卫星点函数逼近的总体误差通过加权可表示为

式中，ωk为权函数.记这M个卫星点k到中心点i的距离分别为r1，r2，…，rM，并将其中距离最大的rmax定义为该点云的参考半径.那么，权函数ωk可定义为ωk＝1／r－2k，其中r－k＝rk／rmax.

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

那么逼近函数可整理为

式中，αk、βk和γk仅与点云Ci中卫星点及中心点的坐标相关，在时间推进迭代计算前可以一次求出.由式（6）可以确定函数f在点i处的空间导数为

同样地，函数f在点i处的空间导数a1、a2和a3也可以通过中心点与卫星点连线中点ik处的函数值fik进行逼近表示：

式（8）中的系数αik、βik和γik也仅与点云Ci中离散点的坐标相关，在时间推进迭代计算前也可以一次求出.借鉴文献［9］的做法，通量运算是在ik处进行的，因此，本文将用式（8）来处理空间离散涉及的通量运算.

1.4通量运算

基于点云Ci，应用式（8），则中心点i处的通量相关项可近似为

鉴于∑（αikF1i＋βikF2i＋γikF3i）为已知，这里只介绍∑（αikF1ik＋βikF2ik＋γikF3ik）的处理方法.在每一个卫星点k与中心点i连线的中点ik处，定义一个数值通量

采用近似黎曼解确定Qik，即

式中：rik为从点i指向点k的矢量；U表示上述矢量中的任一分量；▽Ui和▽Uk为计算点上物理量的梯度，由式（7）计算得到.

1.5时间推进与边界条件

空间离散完成后，在点云Ci上，可得到控制方程的半离散形式为

式中，Ri为计算点上的残差.对于式（15），采用四步Runge－Kutta法进行时间推进求解［9］.

数值求解时还涉及到边界条件.本文物面处应用良导体边界条件［11］，在计算域的截断处则采用完全匹配层（Perfectly Matched Layer，PML）边界条件，相关公式及参数取值详见文献［12］.

2 算例与分析

本文基于上述算法已研制出对应的计算程序，这里结合算例对算法进行考核与分析.以下算例所涉及的计算域中的x、y和z坐标均以波长λ为特征长度进行无量纲化.如图2（a）所示，沿用文献［13］的做法，假定平面入射波沿k轴方向照射目标，并定义k轴与z轴之间的夹角为θ，将k轴在xoy平面内的投影定义为k′轴，并将k′轴与x轴之间的夹角定义为φ；以k轴为坐标轴建立球坐标系，并记各坐标轴方向的单位矢量分别为er、eφ和eθ.如图2（b）所示，Ei和Hi分别为入射电场矢量和入射磁场矢量，α为Ei与eθ的夹角（称为极化角）.那么，归一化的入射波电场分量可写为：

式中：E0＝cos［2π（k－t）］，k＝xsinθcosφ＋ysin θsinφ＋zcosθ.

入射波磁场分量计算公式可类似给出，具体参见文献［13］.

图2 线极化平面入射波

2.1球的散射

先选用有级数解［14］供比较的金属球散射算例对算法进行考核.计算时金属球半径取为λ，计算域大小设置为8λ×8λ×8λ，总点数为463 056，平面入射波沿x轴方向传播，对应θ＝90°，φ＝0°，α＝0°.计算得到的金属球的散射场（z＝0截面）及其对应的双站雷达散射截面（Radar Cross Section，RCS）分布分别见图3和图4.从图4中不难看出，本文结果的峰值及其在整个双站角范围内的分布都与级数解吻合.

图3 金属球的散射电场分布（z＝0截面）

图4 金属球的双站RCS分布（θ＝90°）

为了定量分析计算结果的误差，这里定义计算的算术平均误差EAM和平均相对百分误差EARP为：

式中，σi，ser和σi，cal分别为双站RCS的级数解和计算值.本例计算的EAM为0.095dB，EARP也控制在了1%以内.

必须指出，本例属于时谐场数值模拟，计算最终得到的是稳定的周期解.图5给出了周期间隔的残值收敛历程，图中横坐标为迭代周期数，纵坐标为平均残值的对数值.从图中不难看出，当迭代计算达到48个周期时（相当于入射波在自由空间传播了48个波长距离），平均残值已经由初始的10－1量级下降到10－6量级，这与文献［13］认为的一般入射波传播距离至少达到计算域线性尺度的6倍才能收敛到时谐场的稳态（对应本例约为48个波长距离）是一致的.

图5 残值收敛历程

2.2立方体的散射

接着选用文献［15］的立方体散射算例来进一步考核算法.为了便于同文献［15］的结果进行比较，本例立方体的边长取为λ，置于6λ×6λ×6λ计算域中进行数值模拟，空间布点如图6所示，图中可看出立方体表面布点及z＝－0.5截面附近点的空间分布，本算例总点数为113 225，入射波的设置同上.

图7给出了计算得到的立方体的散射场分布（z＝0截面），其远场外推得到的双站RCS分布见图8.从图8中不难看出，本文计算得到的双站RCS的峰值及其在整个双站角范围内的分布与文献［15］中矩量法的结果一致.

图6 立方体的表面布点及z＝－0.5截面附近点的空间分布

图7 立方体的散射电场分布（z＝0截面）

图8 立方体的双站RCS分布（θ＝90°）

2.3进气道模型的散射

这里给出进气道模型模拟的复杂散射体算例，以展示算法处理三维散射目标的能力.本例进气道模型（见图9）源自文献［16］，长4λ，数值模拟时计算域取为8λ×8λ×8λ，总点数为736 587，线极化平面入射波设置为θ＝45°，φ＝0°，α＝0°.图10给出了计算得到的进气道模型在y＝0截面的散射场分布.从图中可以看到进气道内部散射场的干扰情况，强散射方向出现在θ＝45°和θ＝135°.从图10可以看到，在y＝0截面，进气道模型从几何上已被截成弦长为4λ的双NACA0012翼型.为了展示二维与三维的差异，本文也计算给出了该二维双NACA0012翼型的散射场供比较.图11给出了进气道的双站RCS分布.从图中可以看到，本文结果与文献［16］结果大体一致，最强散射都出现在θ＝45°方向，约为26dB.进气道三维实体的散射特性与对应二维截面的差异一定程度上可从图11或从图12与图10的比较中看出.

图9 进气道模型及其表面布点

图10 进气道模型的散射场分布（y＝0截面）

图11 进气道模型的双站RCS分布（φ＝0°）

图12 4λ弦长的双NACA0012翼型散射电场分布

3 结 论

本文结合三维无网格点云结构及加权最小二乘函数逼近法，发展出求解三维麦克斯韦方程的时域无网格算法，并成功地用于金属球和立方体等三维典型散射体的散射场计算.算例表明发展的算法计算结果能与级数解、文献矩量法或精确控制法等结果一致.给出的进气道模型模拟的复杂散射体算例，一定程度上展示出发展的算法处理三维实际散射目标的能力.

［1］ 苏 敏，陈 刚，童创明.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）

［2］ 邓 聪，尹文禄，柴舜连，等.两种适用于时域有限体积方法的截断边界［J］.电波科学学报，2009，24（3）：537－540＋545.

DENG Cong，YIN Wenlu，CHAI Shunlian，et al.Two boundary conditions for the truncation of cell－centered FVTD algorithm on unstructured lattices［J］.Chinese Journal of Radio Science，2009，24（3）：537－540＋545.（in Chinese）

［3］ 陈红全.求解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）

［4］ 赵美玲，聂玉峰，林世明.电磁场计算中的MLPG－FE耦合新方法［J］.电波科学学报，2006，21（6）：959－964.

ZHAO Meiling，NIE Yufeng，LIN Shiming.A new coupled meshless local Petrov－Galerkin and finite element（MLPG－FE）method for electromagnetic field computations［J］.Chinese Journal of Radio Science，2006，21（6）：959－964.（in Chinese）

［5］ 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.

［6］ LAI Shengjian，WANG Bingzhong，DUAN Yong.Meshless radial basis function method for transient electromagnetic computations［J］.IEEE Transactions on Magnetics，2008，44（10）：2288－2295.

［7］ 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.

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

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

［10］ 高煜堃，陈红全，蒲赛虎.麦克斯韦方程时域无网格算法及其应用研究［J］.电波科学学报，2013，28（6）：1013－1020.

GAO Yukun，CHEN Hongquan，PU Saihu.Meshless time－domain method and its applications for solving Maxwell’s equations［J］.Chinese Journal of Radio Science，2013，28（6）：1013－1020.（in Chinese）

［11］ 高煜堃，陈红全.基于非结构网格格点FVTD算法的电磁散射模拟［J］.南京航空航天大学学报，2013，45（3）：415－423.

GAO Yukun，CHEN Hongquan.Electromagnetic scattering simulation based on cell－vertex unstructured－grid FVTD algorithm［J］.Journal of Nanjing University of Aeronautics ＆Astronautics，2013，45（3）：415－423.（in Chinese）

［12］ FAN Guoxin，LIU Qinghuo.A strongly well－posed PML in lossy media［J］.IEEE Antennas and Wireless Propagation Letters，2003，2：97－100.

［13］ 葛德彪，闫玉波.电磁波时域有限差分方法［M］.西安：西安电子科技大学出版社，2011.

［14］ 葛德彪，魏 兵.电磁波理论［M］.北京：科学出版社，2011.

［15］ 王 志，江谷传.高阶FDTD方法在三维散射问题中的应用［J］.合肥工业大学学报，2009，32（11）：1776－1779.

WANG Zhi，JIANG Guchuan.High－order FDTD for calculating three－dimensional scattering problems［J］.Journal of Hefei University of Technology，2009，32（11）：1776－1779.（in Chinese）

［16］ BRISTEAU M O，GLOWINSKI R，PÉRIAUX J.Controllability methods for the calculation of time－periodic solutions；application to scattering［J］.Journal of Computational Physics，1998，147（2）：265－292.

A meshless time－domain algorithm for solving the 3－D Maxwell’s equations

GAO Yukun CHEN Hongquan
（College of Aerospace Engineering，Nanjing University of Aeronautics ＆Astronautics，Nanjing Jiangsu 210016，China）

A meshless time－domain algorithm for a solution of the 3－D Maxwell’s equations is developed.The spatial derivatives related to the algorithm are approximated by using an expanded Taylor series and the weighted least square technique in each cloud of points，and then a particular approximate Riemann solver is constructed for computing the physical flux of the governing equations.After that，an explicit four－stage Runge－Kutta scheme is used to advance the Maxwell’s equations in time.Combined with solving the 3－D Maxwell’s equations，the implementations of the present algorithm are described in details.Based on the developed algorithm，numerical results for typical 3－D objects such as a metal sphere，a metal cube and an air－intake model are presented，which show that the obtained bistatic radar cross sections are in good agreement with the series solution or that of the reference.

meshless time－domain method；Maxwell’s equations；cloud of points；the weighted least square；three－dimension

O441.4

A

1005－0388（2015）02－0257－07

高煜堃（1984－），男，江苏人，南京航空航天大学博士研究生，主要研究方向为计算电磁学.

陈红全（1962－），男，浙江人，南京航空航天大学教授、博士生导师，主要研究方向为计算流体力学及计算电磁学.

高煜堃，陈红全.求解三维麦克斯韦方程的时域无网格算法研究［J］.电波科学学报，2015，30（2）：257－263.

10.13443／j.cjors.2014050701

GAO Yukun，CHEN Hongquan.Ameshless time－domain algorithm for solving the 3－D Maxwell’s equations［J］.Chinese Journal of Radio Science，2015，30（2）：257－263.（in Chinese）.doi：10.13443／j.cjors.2014050701

2014－05－07

国家自然科学基金（No.11172134）；江苏高校优势学科建设工程

联系人：陈红全E－mail：hqchenam＠nuaa.edu.cn