求解三維麥克斯韋方程的時域無網格算法研究

高煜堃 陳紅全

求解三維麥克斯韋方程的時域無網格算法研究

高煜堃 陳紅全

（南京航空航天大學航空宇航學院，江蘇南京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