一種基于不連續(xù)伽遼金方法求解多區(qū)域目標散射問題的優(yōu)化預處理器

張慧雯 黃曉偉 吳比翼 盛新慶

（北京理工大學信息與電子學院應用電磁研究所，北京 100081）

引 言

在電磁計算應用中，由介質和金屬結合的多區(qū)域結構有著廣泛的應用價值.典型的介質金屬結合的多區(qū)域目標，例如天線-天線罩系統(tǒng)、微波天線陣列等，無法通過建立單一等效模型得到準確的電磁計算結果. 而一些曾經(jīng)使用單一區(qū)域模型來計算電磁特性的目標，也可以通過細化區(qū)域結構的方式得到更精確的結果，比如飛行器、導彈模型等. 因此，研究多區(qū)域目標的高效、高精度電磁計算方法是十分必要的.

為了提高多區(qū)域目標的電磁計算效率，我們使用面積分(surface integral equation, SIE)[1]矩量法(method of moment, MoM)[2]以更高效地求解分塊均勻目標的散射問題，因為其相較于體積分方程具有更少的未知數(shù). SIE 方程一般使用Rao-Wilton-Glisson(RWG)基函數(shù)進行離散[3]，但RWG 基函數(shù)的連續(xù)性使得其很難處理多區(qū)域目標的區(qū)域結合邊界[4]. 通常的解決方法是在邊界處增加區(qū)域連接模型(connectregion modeling, CRM)以保證邊界的連續(xù)性條件，但這需要共型區(qū)域剖分[5]. 然而當目標具有精細幾何結構時，共型剖分是十分繁瑣且耗費時間的[6-7].

為解決上述問題，區(qū)域分解技術(domain decomposition method, DDM)被前人提出，以在不同區(qū)域中使用非共型且相互獨立的幾何剖分[8]SIE 方程的區(qū)域分解可以分為兩大類：封閉式區(qū)域分解(體分解) (volume-based SIE-DDM)和非封閉式區(qū)域分解(面分解)(surface-based SIE-DDM)[9]. 封閉式區(qū)域分解以體積為基本區(qū)域單元，在子區(qū)域的連接處需要施加Robin傳輸條件(TCs)以保證子區(qū)域電流的連續(xù)性[10-11]，而這將帶來更多的未知數(shù). 非封閉式區(qū)域分解方法中，不連續(xù)伽遼金 (discontinuous Galerkin, DG)方法的研究最為深入[12-13]. 在不同面區(qū)域的邊界處，DG 方法使用半RWG 基函數(shù)進行剖分，同時將電流的連續(xù)性條件直接加在半RWG 的阻抗矩陣元素上. DG 方法無需增加未知數(shù)，因此該方法的非共型剖分可以更容易地運用在復雜目標中. 近年來，學者已將DG 方法成功運用于完美電導體[14-15]、均勻介質[16]和阻抗邊界表面[17-18]等目標的計算中，而對于多區(qū)域結構的目標，該方法還未有深入的研究.

對于電大目標的仿真，預處理技術是一種加速迭代收斂的有效方法[19-20]. 近場預處理器(near-field preconditioner, NFP)是多層快速多極子方法(multilevel fast multipole algorithm, MLFMA)[21]常用的預處理方式. 對于復雜的多區(qū)域目標，其近場矩陣元素非常多，增加了預處理矩陣求逆的運算時間. 距離稀疏預處理器(distance sparse preconditioner,DSP)使用距離限制減少了預處理器的矩陣元素，可以明顯減少計算時間和內存需求[22]. 進一步，我們提出了優(yōu)化距離稀疏預處理器(optimized DSP, ODSP)，不同積分算子的矩陣塊，其數(shù)值大小隨距離變化的趨勢不盡相同. 算子的預處理矩陣塊可以保留更少的矩陣元素而并不會大程度削弱其預處理的效果. 對于多區(qū)域目標，O-DSP 對不同的子區(qū)域矩陣塊使用不同的構建方法，使之合成為完整的預處理矩陣. 與DSP 相比，O-DSP 保持了相近的加速收斂效果，而且構建的預處理矩陣元素數(shù)目大幅降低.

1 理論與公式

1.1 多區(qū)域目標的電磁散射計算

考慮圖1 所示的交界面模型，這里 Ωi表示不同的體積區(qū)域，可以是均勻介質或金屬. 金屬 δi=0，均勻介質 δi=1. 任意兩個空間區(qū)域的交界面用S ij表示，默認法相方向n?ij以區(qū)分面的外表面和內表面. 若交界面的兩個空間區(qū)域至少一個是金屬區(qū)域，則只需要在表面上定義電流，用 Δij=0表示；當交界面的兩個空間區(qū)域均為介質時，還需要定義磁流，用Δij=1表示.

圖1 任意交界表面與其兩側體積區(qū)域示意圖Fig. 1 An arbitrary surface with its internal and external regions

1.2 DG 方法

在式(1)、式(2)所表示的原方程中所有邊界上建立方程(10)，并乘以 ?jαijαii0/ki加到原方程中，用同樣的基函數(shù)進行離散，以消去式(7)L算子元素中的第五項基函數(shù)是半RWG 的計算.

對于復雜目標，使用MLFMA 來加速計算[23]. 對于DG 建立的離散方程同樣適用. 唯一的區(qū)別在于邊界的半RWG 需要進行一個三角形的積分計算.

1.3 預處理方法

多區(qū)域目標SIE 方程迭代計算的收斂速度較慢，因此需要進行合適的預處理. 預處理矩陣的定義如下：

圖2 多區(qū)域交界處示意圖Fig. 2 The junction of multiple regions

圖3 介質球體不同算子矩陣塊內元素模值隨距離約束的變化Fig. 3 The variation of modulus value with distance constraint of a dielectric sphere

數(shù)值實驗發(fā)現(xiàn)，在K 算子矩陣塊中使用式(13)的預處理方式和原DSP 的預處理方式在迭代速率上相近，但是預處理矩陣元素進一步減少. 結合之前L算子矩陣塊的預處理方法，簡稱這種方法為 ODSP.

矩陣塊中的未知數(shù)為N，則近場矩陣元素個數(shù)為O(N2). 式(12) 中距離判斷與介質波長有關，因此其矩陣元素數(shù)量為O(N)，由于有很多元素仍被約束在距離判斷內，O(N)的實際系數(shù)很大. 式(13) O-DSP 中矩陣每一行只會保留4 個或2 個元素，因此整體元素數(shù)將小于4N，最終整體的預處理矩陣元素數(shù)將會明顯減少，且介質表面未知數(shù)越多，減少效果越明顯.

2 數(shù)值結果

利用廣義最小殘差 (generalized minimum residual,GMRES)求解器來進行仿真，Krylov 子空間的大小為100，迭代殘差為1 0?3. 預處理求逆使用大規(guī)模并行求解器(multifrontal massively parallel solver, MUMPS)，計算平臺為2.2 GHz CPUs 和 1 TB 內存的 Intel Xeon E7-4850.

如圖4 所示，半徑為0.5 m 的介質球，相對介電常數(shù)為5，被分為4 個子區(qū)域，剖分尺寸為0.005 m，為介質波長的11/100，頻率為300 MHz 的平面波沿-Z方向入射. 利用DSP 和O-DSP 對 VV 極化方向的雙棧雷達散射截面(radar cross section, RCS)(記為VV-RCS)進行對比計算，計算結果如圖5 所示. 可以看出，DSP 和O-DSP 數(shù)值計算結果均很好地吻合了理論計算的Mie 求解，不會對求解結果的精度產(chǎn)生影響.

圖4 介質球的區(qū)域分解形式Fig. 4 Domain partitioning for the dielectric sphere

圖5 介質球的雙棧VV-RCSFig. 5 The bistatic VV-RCS of the dielectric sphere

L算子的距離約束不變，為介質波長的1 /2，令K算子矩陣塊內元素受不同的距離約束，對比其預處理內存和迭代次數(shù)，結果如圖6 所示. 可以看到，當距離閾值在介質波長的1/10～1/2 時，迭代次數(shù)幾乎不變，但元素個數(shù)隨距離的減小大大減少，說明去掉含奇異點殘留項的元素不會對迭代產(chǎn)生很大的影響.當距離閾值為介質波長的3/50 時，處于是否包含奇異點殘留元素的臨界，迭代次數(shù)發(fā)生了較為明顯的改變，因此有保留這些元素的必要性.

圖6 K 算子矩陣迭代次數(shù)和元素數(shù)量隨距離約束的變化Fig. 6 The iteration numbers and element numbers of K operator matrix block, as a function of distance constraint

表1 對幾種不同預處理器的數(shù)值效果進行了對比. 對比的稀疏策略(sparse strategies, ST)有如下四種：

表1 不同預處理器的DG 方法計算數(shù)值效果對比Tab. 1 Comparison of the n umerical performance of the DG solution with different preconditioner

由表1 可知，和NP 相比，所有的預處理器都增加了收斂速度，其中NFP，DSP，O-DSP 效果最好. 這三種收斂速度較好的預處理器中，O-DSP 矩陣元素數(shù)量最少. 和DS-2BDP 相比O-DSP 收斂速度極快，說明K算子矩陣塊中有奇異點殘留項的元素對于矩陣的收斂十分重要，即使數(shù)量并不多. 同樣使用ST3 作用在L算子矩陣塊，迭代速率明顯下降，說明L算子矩陣塊需要保留較多的矩陣元素.

圖7 為圖4 中的介質球子區(qū)域數(shù)量分別為4、8、16 和32 的分解情況. 不同分解區(qū)域數(shù)對DSP 和O-DSP 的迭代次數(shù)和元素比率的影響如表2 所示.區(qū)域數(shù)越多意味著邊界半RWG 越多，K算子矩陣塊會保留更少的元素. 從表2 中可以看出，兩個預處理器的迭代速率幾乎是一致的，說明子區(qū)域的數(shù)量并不會影響O-DSP 的效率.

圖7 介質球的不同區(qū)域分解Fig. 7 Domain partition schemes for a dielectric sphere

表2 不同分解區(qū)域數(shù)對DSP 和O-DSP 的迭代次數(shù)和元素比率的影響Tab. 2 Comparison of the iteration numbers and element ratios of the DG solution with DSP and O-DSP for different decomposed spheres

如圖8 所示，區(qū)域數(shù)量為4 的介質球內嵌不同個數(shù)半徑為0.15 m、剖分尺寸均為0.005 m 的金屬球.頻率為300 MHz 的平面波沿-Z方向入射，介質球內嵌不同個數(shù)的導體塊模型的雙棧VV-RCS 結果如圖9所示. 從圖9 可以看出，使用DG 和O-DSP 的預處理器計算的數(shù)值結果很好地吻合了商業(yè)軟件FEKO 使用共型剖分計算的結果. 導體塊數(shù)量分別為2，4 和8 時三種預處理器迭代次數(shù)和元素個數(shù)對比如表3所示. 對比DSP，O-DSP 的迭代數(shù)有輕微的增加，而矩陣元素數(shù)量明顯減少. 由此說明，O-DSP 仍然適用于多區(qū)域的目標.

圖8 介質球內嵌不同個數(shù)的導體塊模型Fig. 8 Different numbers of conducting spheres embedded in one dielectric sphere

圖9 介質球內嵌不同個數(shù)導體塊模型的雙棧VV-RCSFig. 9 The bistatic VV-RCS for the model of multiple conducting spheres embedded in one dielectric sphere

表3 導體塊數(shù)量不同情況下三種預處理器迭代次數(shù)和元素個數(shù)對比Tab. 3 Comparison of the iteration numbers and element numbers by 3 different preconditioners with different numbers of conducting spheres

圖10 所示為一個簡化的射頻芯片模型，包含了所有射頻芯片內部的基本結構. 介質長方體模型尺寸為20 mm×20 mm×5 mm ，介電常數(shù)為 2.5，內嵌間距為 1 mm 的 5 片金屬片，其中第一層和第五層用 4個金屬圓柱相接作為整體的支撐，第一層和第二層之間由過孔相連，第三層和第四層之間是介電常數(shù)為4.2 的介質體. 剖分尺寸為細金屬柱0.3 mm，其余部分0.5 mm. 一共分為 17 個 DG 區(qū)域，總未知數(shù)為 85 970 個.

頻率為 40 GHz 的平面波沿-Z方向入射，使用CTF 方程計算，使用O-DSP 作為預處理器， VV極化方向的雙棧 RCS 如圖11 所示，驗證了方程的計算正確性. 迭代收斂速度如圖12 所示， DSP 和O-DSP 的迭代次數(shù)分別是 149 和153 ，略有差異，兩者明顯快于NP 的收斂速度. DSP 的矩陣元素數(shù)量是 78 061 756，O-DSP 為 46 922 730，數(shù)量明顯減少.

圖11 射頻芯片模型的雙棧VV-RCSFig. 11 The bistatic VV-RCS of the radio frequency chip model

圖12 使用不同預處理器求解射頻芯片模型的迭代收斂情況Fig. 12 Convergence histories for the radio frequency chip model under different preconditioners

圖13 所示為四旋翼無人機模型的區(qū)域分解示意圖，四個旋翼和支架是金屬的，機身是相對介電常數(shù)為4 的介質體. 飛行器長寬約為2 000 mm 和830 mm.12 個子區(qū)域的總未知數(shù)為415 223 . 介質處的剖分為6 mm，一般導體的剖分為8 mm，旋翼處由于太細剖分尺寸為6 mm.

圖13 四旋翼無人機模型的區(qū)域分解示意圖Fig. 13 Illustration of subdomain partition scheme of a fourrotor aircraft model

頻率為 4 GHz 的平面波沿-Z方向入射， VV極化方向的雙棧 RCS 如圖14 所示. 可以看出，O-DSP的 DG 方法和FEKO 符合良好.

圖14 四旋翼無人機模型的雙棧VV-RCSFig. 14 The bistatic VV-RCS of the four-rotor aircraft model

圖15 為利用NP，DSP 和O-DSP 求解四旋翼無人機模型的迭代收斂情況. 可以看出，DSP 和O-DSP均快速收斂，二者速度相近. 表4 為 DSP 和O-DSP預處理器計算四旋翼無人機迭代次數(shù)和矩陣元素數(shù)量對比，其中O-DSP 矩陣元素數(shù)量明顯減少.

圖15 使用不同預處理器求解四旋翼無人機模型的迭代收斂情況Fig. 15 Convergence histories for the four-rotor aircraft model under different preconditioners

表4 DSP 和O-DSP 預處理器計算四旋翼無人機迭代次數(shù)和矩陣元素數(shù)量對比Tab. 4 Comparison of the iteration numbers and element numbers of the DG solution with DSP and O-DSP for the fourrotor aircraft model

3 結 論

本文使用DG 方法來求解多區(qū)域目標問題，并使用O-DSP 來提高仿真效率. 優(yōu)化的預處理矩陣根據(jù)積分微分算子的不同數(shù)值特性，分塊進行更恰當?shù)念A處理構建. 數(shù)值算例證明了預處理器的靈活性和準確性. O-DSP 相比DSP 擁有相近的迭代次數(shù)和更稀疏的結構特性，由于K算子的預處理矩陣塊使用了一種更精準的判斷方式保留了恰當?shù)脑?，因此ODSP 是一種更加有效的計算多區(qū)域目標散射的預處理方式.