基于WENO重構的時域無網格算法及其應用研究

高煜堃，陳紅全

(南京航空航天大學航空宇航學院 南京 210016)

基于WENO重構的時域無網格算法及其應用研究

高煜堃，陳紅全

(南京航空航天大學航空宇航學院 南京 210016)

研究了用于求解電磁散射問題的WENO重構時域無網格算法?；跓o網格點云結構，引入沿點云中心點和衛星點連線方向的局部一維坐標，并結合虛擬點的設置，構成在點云中實施三階WENO重構所要求的模板，以便利用WENO重構計算中心點與衛星點連線中點處物理量的左右狀態值，供通量運算。設置的虛擬點上的物理量則利用點云中已有的最近點插值系數直接插值確定。用發展的算法對典型的一維和二維電磁散射問題進行數值模擬，并與理論解和傳統的基于線性重構的計算結果進行比較，驗證了WENO重構獲得的數值解比線性重構更接近理論解。給出了多體干擾電磁散射場算例，展示出用發展的算法處理復雜多體干擾情形的效果。

點云; 電磁散射場; 時域無網格算法; 多體干擾; WENO重構

無網格方法是繼網格方法之后出現的一種新型數值方法。計算區域的離散只涉及布點填充，既可以采用網格點，也可以打破傳統的網格方法[1]中網格拓撲的約束而根據需要直接布點，因此，該方法在對復雜外形的處理上更加靈活[2]。近年來，無網格方法已在計算電磁學(computational electromagnetics, CEM)領域得到了應用和發展，用于求解電磁場問題的代表性方法有局部弱式Petrov-Galerkin無網格法和徑向基點插值無網格法等[3-5]。這些方法通常會涉及形函數(或基函數)的選取。形函數的選擇不僅會影響矩陣的質量[6]，還會影響到邊界條件的施加，有時邊界條件需要采用特殊的方法強制給定[7]。在計算流體力學(computational fluid dynamics, CFD)領域，近年來也出現了與CEM中做法不同的無網格方法。無網格點上的空間導數可基于局部點云結構通過極小曲面逼近得到?？臻g離散涉及的通量運算常采用近似黎曼解處理。該類無網格方法已成功應用于求解Euler方程，能模擬出復雜外形的繞流[2]。

鑒于所研究的麥克斯韋方程和流體力學中的Euler方程具有雙曲型等相同的數學特性，本文意在借鑒上述CFD中的無網格算法，發展基于加權本質無振蕩(weighted essentially non-oscillatory, WENO)重構的時域無網格算法，用于求解麥克斯韋方程。CFD中的WENO重構是在本質無振蕩(essentially non-oscillatory，ENO)重構方法[8]的基礎上發展而來的高階重構方法，大多用于結構網格[9]。本文基于無網格點云和局部一維坐標，給出了WENO重構實施所要求模板的構造方法。通量運算涉及的無網格點云中心點與衛星點連線中點處物理量的左右狀態值，改用三階WENO重構計算，以替代通常的線性重構?？臻g離散涉及的通量運算采用基于點云結構的近似黎曼解方法確定，時間離散則按照四步Runge-Kutta法推進求解。然后結合算例分析，驗證了本文基于WENO重構的數值解比通常的線性重構結果更接近理論解。最后，本文研究了雙NACA0012翼型及飛機投彈模擬的多體干擾電磁散射特性，在一定程度上展示了算法處理多體干擾等復雜情形的能力。

1 時域麥克斯韋方程

對于橫磁(transverse magnetic, TM)波，在直角坐標系中守恒型無量綱時域麥克斯韋方程可寫為：

2 時域麥克斯韋方程的無網格離散求解

2.1 計算區域布點離散及點云生成

無網格算法計算區域的離散只涉及布點填充，既可以方便地選取結構網格或非結構網格點，也可以打破傳統的網格拓撲約束，按要求直接布點。計算區域布點填充后，在點的局部需構成無網格算法實施所要求的點云結構。以點云ci為例，如圖1所示，中心點i自然納入點云ci中，再合理選取一定數目的衛星點。本文以中心點i為圓心，Ri為半徑畫圓，將圓內除中心點i以外的點(1、2、3、4、5和6)視為衛星點一并納入點云ci中。其中圓半徑Ri與點云當地點與點之間的距離相關，通?？赏ㄟ^擴大或縮小ir來選取合適數目的衛星點。

圖1 無網格點云ci示意圖

2.2 空間導數逼近

當點云生成后，在點云Ci上假定中心點i附近的函數值分布滿足函數f=f(x,y)。那么，沿用文獻[10]的方法，f的二次逼近函數可寫為：

式中，fi=f(xi,yi)；h=x?xi；l=y?yi；ai(i=1,2,,5)為待定系數，表示在中心點i處相應的各階偏導數。在點云Ci上，記衛星點k(k=1,2,,M)處的函數值為fk，則定義函數逼近的總體誤差為：

基于總體誤差G極小，逼近函數可整理為[10]：

式中，系數αk、βk、γk、ηk和ζk僅與離散點的坐標相關，可在迭代計算前一次計算得到。函數f在中心點i處的空間導數ai(i=1,2,,5)可表示為：

同樣地，也可以用中心點與衛星點連線中點ik處的函數值fik進行表示：

式中，整理得到的系數αik、βik、γik、ηik和ζik也僅與離散點的坐標相關，可在迭代計算前一次計算得到。本文將根據式(6)處理下面空間離散涉及的通量運算。

2.3 基于無網格點云的WENO重構

為了提高時域無網格算法的計算精度，本文借鑒CFD中三階WENO重構的有限差分方法[11]構造在點云中實施WENO重構所要求的模板。為此，基于點云結構，在點云ci上以中心點i和衛星點k為例，建立局部一維坐標，如圖2所示。圖中點ii是點i關于點k的對稱點(虛擬點)，點kk是點k關于點i的對稱點。記點i和點k處的函數值分別為fi和fk，虛擬點ii和kk處的函數值分別記為fii和fkk，則點i和點k的中點ik處的函數值fLik和fRik可按照下式WENO重構進行計算：

式中，加權因子ξ1=1/3；ξ2=2/3。于是式(7)可進一步整理為：

圖2 基于無網格點云的WENO重構模板示意圖

對于虛擬點上的函數值，最直接的做法是針對虛擬點，逐一選點并構造插值模板，再按照插值方法確定。但是，要逐一搜索插值點，且在計算插值系數的過程中還會涉及矩陣求逆，都需要花費較多的計算時間。為了減小不必要的計算量，本文利用已知的離虛擬點最近的點云結構及其空間導數信息，通過采用形如式(4)的二次函數逼近來插值計算出虛擬點上的函數值。

2.4 通量運算

在ci上，應用式(6)近似中心點i處的通量為：

再采用近似黎曼解確定數值通量Qik，即：

2.5 時間離散

采用無網格方法對麥克斯韋方程進行空間離散后，在點云Ci上，可以得到其半離散形式為：

式中，Ri為點云中心點，即計算點上的殘差。本文利用式(13)按照4步Runge-Kutta格式[10]進行時間推進求解。迭代計算時，如式(9)～式(12)所示，殘差計算涉及的各空間離散項則利用當前計算值計算確定。對于計算涉及的邊界條件，本文在物面應用良導體邊界條件，而在計算域的截斷位置則采用完全匹配層邊界條件，具體公式及參數取值詳見文獻[10]。

3 算例與分析

本文按上述方法研制了對應的計算程序，現結合算例進行電磁散射模擬。算例涉及的計算域中x、y坐標均為以入射波波長λ為特征長度、無量綱化后的空間位置。如圖3所示，TM平面波沿k方向照射目標，并定義k軸與x軸之間的夾角為入射角?，那么歸一化后的入射波分量可表示為：

式中，k=xcos?+ysin ?；t為無量綱參數。

圖3 TM平面波

3.1 一維電磁波傳播數值模擬

本文先選用具有解析解的一維電磁波傳播問題對發展的WENO重構算法進行考核。假定電磁波沿x方向傳播，那么式(1)的麥克斯韋方程可簡化為：

對應的解析解為：

在t=10時刻提取相應的電場信息進行研究。定義計算誤差為：

式中，Ez,cal為計算值。

圖4 計算誤差與布點密度的關系

圖4給出了計算誤差與布點密度的關系。整體上看，不論基于WENO重構還是線性重構，error隨 dom的增加而減??；dom一定時，本文數值解的error顯著小于線性重構的結果，計算得到的電場幅值和相位均與理論解更加接近。

3.2 二維圓柱散射數值模擬

本文選用有級數解[12]可供比較的金屬圓柱算例進行電磁散射模擬。先采用如圖5a所示的均勻布點方式，布點密度dom20=。數值模擬時，圓柱半徑為2.0λ，計算域設置為12λ×12λ，總點數為60 698，TM波入射角?=0°。計算得到的圓柱散射場和雙站雷達散射截面(radar cross section, RCS)分布分別如圖6和圖7所示，圖7給出了級數解比較。從圖7不難看出，本文基于WENO重構計算得到的雙站RCS和傳統基于線性重構的計算結果，其峰值及其在整個雙站角范圍內的分布整體上都能與級數解一致；但是在各極值點處，本文結果與級數解吻合更好，如圖7b所示。

圖5 圓柱計算點云示意圖

圖6 圓柱的散射場等值線分布

圖7 圓柱的雙站RCS分布及局部放大圖

圖8 不同布點密度對應的圓柱雙站RCS分布及局部放大圖

本例均勻布點方式離散計算區域總點數達60 698，計算量較大?？紤]到本文在物面處引入入射波及提取時域信息，因此物面附近的布點密度應予以保證，而對于截斷邊界附近，布點密度適當放寬。改用如圖5b所示的非均勻布點方式重新計算圓柱散射場，其他條件設置同圖5a，計算得到的雙站RCS分布如圖8所示。從圖中不難看出，當dom10=時，計算得到的雙站RCS分布與均勻布點的結果基本吻合；當dom8=時，除個別極值點處數值有差異外，計算值大體上也能與均勻布點情形一致；當dom5=時，計算值與均勻布點的結果存在明顯差異，如圖8b所示。當dom10=時，總點數由均勻布點的60 698下降到16 380，計算量顯著減小，而計算結果幾乎不受影響。

3.3 雙NACA0012翼型散射場數值模擬

選用雙NACA0012翼型算例[13]驗證本文算法處理多體干擾問題的能力。沿用文獻[13]的做法，取NACA0012翼型弦長為4λ，雙翼型弦線之間的距離為λ。計算域取12λ×12λ，電磁波入射角分別取?=0°,45°，對應計算得到的雙NACA0012翼型的散射場等值線分布如圖9所示，基本與文獻[13]的精確控制法結果一致。從圖9可以看到，在雙翼型之間的區域，散射波疊加作用顯著，電磁干擾較強；而且隨著入射角的改變，干擾作用的效果也存在差異。

圖9 雙NACA0012翼型散射場等值線分布

3.4 飛機投彈多體散射場數值模擬

為了進一步驗證本文算法處理多體干擾問題的能力，給出飛機投彈時出現的復雜多體干擾電磁散射場算例。飛機模型源自文獻[14]，機體長取為8λ，通過對飛機下表面修形來模擬彈艙，炸彈用橢圓體表示。計算域為20λ×20λ，電磁波入射角取?=135°，對應電磁波從飛機前下方照射，計算得到的散射場等值線分布如圖10所示。從圖中可以看出飛機投彈與不投彈時的散射場多體干擾情況。不投彈時，只在飛機前部的下方位置附近，由于機頭與機身的干擾，局部存在散射場的疊加效應外，其他地方的散射場大體上呈現出均勻的環繞分布，強散射方向出現在雙站角135°方向和225°方向；當飛機彈艙打開，進行投彈時，在飛機下方，除了飛機自身各個部件的干擾以外，炸彈與飛機之間也存在強烈的電磁散射干擾，散射場方向發生改變，最強散射僅出現在雙站角135°方向。

圖10 飛機投彈與不投彈時的散射場分布

4 結 束 語

本文研究和發展了求解電磁散射問題的基于WENO重構的時域無網格算法。算例表明：本文算法相較于通常的線性重構算法，計算得到的數值解能更接近理論解。由于該算法是基于無網格點云發展的，計算區域的離散只涉及布點，適合處理多體干擾等復雜情形。

[1] 蘇敏, 陳剛, 童創明. FVTD在電磁散射問題中的應用[J].航天電子對抗, 2008, 24(3): 56-58. SU Min, CHEN Gang, TONG Chuang-ming. Application of FVTD method in CEM problems[J]. Aerospace Electronic Warfare, 2008, 24(3): 56-58.

[2] 陳紅全. 求解Euler方程的隱式無網格算法[J]. 計算物理, 2003, 20(1): 9-13. CHEN Hong-quan. Implicit gridless method for Euler equations[J]. Chinese Journal of Computational Physics, 2003, 20(1): 9-13.

[3] ZHAO Mei-ling, NIE Yu-feng. A study of boundary conditions in the meshless local Petrov-Galerkin (MLPG) method for electromagnetic field computations[J]. CMES, 2008, 37(2): 97-112.

[4] YU Yi-qiang, CHEN Zhi-zhang. A 3D radial point interpolation method for meshless time-domain modeling[J]. IEEE Transactions on Microwave Theory and Techniques, 2009, 57(8): 2015-2020.

[5] LAI Sheng-jian, WANG Bing-zhong, DUAN Yong. Meshless radial basis function method for transient electromagnetic computations[J]. IEEE Transactions on Magnetics, 2008, 44(10): 2288-2295.

[6] 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.

[7] LIU G R, GU Y T. 無網格法理論及程序設計[M]. 王建明,周學軍, 譯. 濟南: 山東大學出版社, 2008. LIU G R, GU Y T. An introduction to meshfree methods and their programming[M]. Translated by WANG Jian-ming, ZHOU Xue-jun. Jinan: Shandong University Press, 2008.

[8] LIU Xu-dong, OSHER S, CHAN T. Weighted essentially non-oscillatory schemes[J]. Journal of Computational Physics, 1994, 115(1): 200-212.

[9] GEROLYMOS G A, SENECHAL D, VALLET I. Very-high-order weno schemes[J]. Journal of Computational Physics, 2009, 228: 8481-8524.

[10] 堃高煜, 陳紅全. 基于非結構網格格點FVTD算法的電磁散射模擬[J]. 南京航空航天大學學報, 2013, 45(3): 415-423. GAO Yu-kun, CHEN Hong-quan. Electromagnetic scattering simulation based on cell-vertex unstructured-grid FVTD algorithm[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2013, 45(3): 415-423.

[11] JIANG Guang-shan, SHU Chi-wang. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics, 1996, 126: 202-228.

[12] HALLEROD T, 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.

[13] CHEN Hong-quan, SUI Hong-tao. New time-explicit asymptotic method for the solution of an exact controllability problem of scattering waves[J]. Chinese Journal of Aeronautics, 2000, 13(2): 80-85.

[14] BESPALOV A. Application of fictitious domain method to the solution of the helmholtz equation in unbounded domain[R]. INRIA 1797. Le Chesnay: INRIA, 1992.

編輯稅 紅

Study of a WENO Reconstruction Based Meshless
Time-Domain Algorithm and Its Applications

GAO Yu-kun and CHEN Hong-quan
(College of Aerospace Engineering, Nanjing University of Aeronautics & Astronautics Nanjing 210016)

A meshless time-domain algorithm based on weighted essentially non-oscillatory (WENO) reconstruction for solving the electromagnetic scattering problems is studied. The stencil required for implementing third-order WENO reconstruction in the gridless cloud is obtained by setting a local coordinate in the direction of each satellite point and introducing a virtual point so as to use WENO reconstruction to approximate the physical quantities at the midpoint between the central and satellite points of the gridless cloud. Additionally, physical quantities at the virtual point are determined by direct interpolation based on available interpolation coefficients of the nearest point. Then, approximate Riemann solver is introduced in dealing with the computation of the flux related to the governing equations, and an explicit four-stage Runge-Kutta scheme is employed in time-marching. After that, based on the developed algorithm, typical 1-D and 2-D cases for solving the electromagnetic fields are simulated. The simulation results verify that the numerical results calculated by using WENO reconstruction are closer to the theoretical solutions than that based on linear function reconstruction. The paper ends with the presentation of the electromagnetic scattering fields with multi-bodies, which show the developed algorithm has the ability to accommodate complex geometries.

cloud of points; electromagnetic scattering field; meshless time-domain method; multi-bodies; WENO reconstruction

O441.4

A doi:10.3969/j.issn.1001-0548.2015.04.013

2013 ? 11 ? 27；

2014 ? 12 ? 12

國家自然科學基金(11172134)；江蘇省高校優勢學科建設工程資助項目

高煜堃(1984 ? )，男，博士生，主要從事計算電磁學方面的研究.