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

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

（北京理工大學(xué)信息與電子學(xué)院應(yīng)用電磁研究所，北京 100081）

引 言

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

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

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

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

1 理論與公式

1.1 多區(qū)域目標(biāo)的電磁散射計(jì)算

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

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

1.2 DG 方法

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

對(duì)于復(fù)雜目標(biāo)，使用MLFMA 來(lái)加速計(jì)算[23]. 對(duì)于DG 建立的離散方程同樣適用. 唯一的區(qū)別在于邊界的半RWG 需要進(jìn)行一個(gè)三角形的積分計(jì)算.

1.3 預(yù)處理方法

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

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

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

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

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

2 數(shù)值結(jié)果

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

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

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

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

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

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

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

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

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

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

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

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

圖8 介質(zhì)球內(nèi)嵌不同個(gè)數(shù)的導(dǎo)體塊模型Fig. 8 Different numbers of conducting spheres embedded in one dielectric sphere

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

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

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

頻率為 40 GHz 的平面波沿-Z方向入射，使用CTF 方程計(jì)算，使用O-DSP 作為預(yù)處理器， VV極化方向的雙棧 RCS 如圖11 所示，驗(yàn)證了方程的計(jì)算正確性. 迭代收斂速度如圖12 所示， DSP 和O-DSP 的迭代次數(shù)分別是 149 和153 ，略有差異，兩者明顯快于NP 的收斂速度. DSP 的矩陣元素?cái)?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 使用不同預(yù)處理器求解射頻芯片模型的迭代收斂情況Fig. 12 Convergence histories for the radio frequency chip model under different preconditioners

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

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

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

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

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

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

表4 DSP 和O-DSP 預(yù)處理器計(jì)算四旋翼無(wú)人機(jī)迭代次數(shù)和矩陣元素?cái)?shù)量對(duì)比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 結(jié) 論

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