周國華，劉勝道，肖昌漢，劉大明

(海軍工程大學(xué)電氣工程學(xué)院，湖北武漢430033)

鋼鐵結(jié)構(gòu)艦船周圍存在的磁異常信號是水中磁性兵器攻擊和空中磁性探測的重要信號源。為了提高艦船生命力，各國海軍都在致力于研究艦船磁性防護(hù)方法及措施，其中重點(diǎn)研究內(nèi)容之一就是艦船受地磁場磁化作用產(chǎn)生的感應(yīng)磁場數(shù)值建模及其防護(hù)問題［1-4］。

鑒于磁場積分法只需剖分源區(qū)、易于處理開域問題等優(yōu)點(diǎn)，常被應(yīng)用于艦船感應(yīng)磁場計(jì)算領(lǐng)域［5-7］。艦船三維復(fù)雜結(jié)構(gòu)使得艦船感應(yīng)磁場數(shù)值計(jì)算中需較多的剖分單元才能獲得相對穩(wěn)定精確的計(jì)算結(jié)果。然而由于磁場積分法系數(shù)矩陣為非對稱稠密陣，基于磁場積分法的船磁計(jì)算時(shí)間隨艦船剖分單元數(shù)目的增加而快速增長，從而使得過去艦船感應(yīng)磁場建模大多停留于小規(guī)模計(jì)算［2-4］，影響了船磁計(jì)算效率及磁場積分法在船磁計(jì)算領(lǐng)域的推廣應(yīng)用。20世紀(jì)末期，Rokhlin和Greengard創(chuàng)造性地提出了加速矩陣與向量乘積運(yùn)算的快速多極子方法(fast multipole method，F(xiàn)MM)［8-9］。隨著快速多極技術(shù)的不斷發(fā)展，求解大規(guī)模問題的快速計(jì)算優(yōu)勢逐漸被進(jìn)一步開發(fā)，目前已被廣泛應(yīng)用于電磁散射、彈性位勢、聲學(xué)技術(shù)等領(lǐng)域［10-12］。

為提高船磁積分方程計(jì)算效率，本文將快速多極子技術(shù)與磁場積分法相耦合，提出了基于快速多極子和積分方程的艦船感應(yīng)磁場快速預(yù)報(bào)方法，設(shè)計(jì)了球體感應(yīng)磁場解析算例和地磁作用下簡化艇模主艙段感應(yīng)磁場實(shí)驗(yàn)。

1 磁場積分法

如圖1所示，在外磁場B0作用下的鐵磁物體(如鐵質(zhì)艦船)在空間點(diǎn)P產(chǎn)生的磁感應(yīng)強(qiáng)度可表示為［13］

式中:▽P為對場點(diǎn)坐標(biāo)的梯度算子，▽Q為對源點(diǎn)坐標(biāo)的梯度算子，v為鐵磁物體所占體積，M為由外磁場磁化引起的鐵磁物體內(nèi)部磁化強(qiáng)度，P為場點(diǎn)，Q為源點(diǎn)。

圖1 靜磁計(jì)算問題示意圖Fig.1 The sketch map of magnetostatic computation

對于均勻磁化體，式(1)中體積分可簡化為面積分

式中:s為均勻磁化體表面積，n為均勻磁化體表面外法線方向。

為得到鐵磁物體內(nèi)部的磁化強(qiáng)度，通常將鐵磁物體離散為若干足夠小的均勻磁化體。以每個(gè)離散單元中心為計(jì)算場點(diǎn)，可建立以單元內(nèi)部磁感應(yīng)強(qiáng)度為未知量的代數(shù)方程組

式中:j=1，2，…，N，N為鐵磁物體離散單元數(shù)，μri為第i個(gè)單元內(nèi)部的相對磁導(dǎo)率。采用矩陣形式可表示為

其中:C為系數(shù)矩陣并由式(3)決定，X為鐵磁物體各離散單元內(nèi)部磁感應(yīng)強(qiáng)度組成的向量，Y為鐵磁物體各離散單元中心處外加磁場組成的向量。求解方程組(4)即可得到鐵磁物體內(nèi)部磁感應(yīng)強(qiáng)度分布，再根據(jù)單元內(nèi)部場量關(guān)系M=(μr－1)B/(μ0μr)及式(1)即可求得空間任意點(diǎn)P的磁感應(yīng)強(qiáng)度值。上述即為用于艦船感應(yīng)磁場計(jì)算的基于磁感應(yīng)強(qiáng)度B的磁場積分法的基本原理。

采用上述磁場積分法求解外磁場作用下的鐵磁物體磁化場問題時(shí)，形成的系數(shù)矩陣為非對稱稠密陣，當(dāng)剖分單元N較大時(shí)，迭代求解代數(shù)方程組(4)將消耗大量的時(shí)間，大大影響了計(jì)算效率。本文提出了基于快速多極子和磁場積分法的艦船感應(yīng)磁場快速預(yù)報(bào)方法，以改善艦船感應(yīng)磁場計(jì)算效率。

2 快速多極子加速磁場積分算法

快速多極子技術(shù)的核心思想是將積分方程核函數(shù)展開成2個(gè)分別僅與源點(diǎn)和場點(diǎn)有關(guān)的函數(shù)，再通過聚合、轉(zhuǎn)移和發(fā)散過程來減少系數(shù)矩陣每次迭代的計(jì)算量，從而提高計(jì)算速度。快速多極子加速積分方程計(jì)算過程是建立在上述磁場積分法和線性代數(shù)方程組迭代法求解基礎(chǔ)上的。如圖2所示(為直觀清晰，采用四叉樹結(jié)構(gòu)繪圖)，該方法先將鐵磁區(qū)進(jìn)行離散，然后用一個(gè)立方體將所有離散單元罩住，并將其等分為8個(gè)子立方體，依次類推直至最小立方體包含的離散單元數(shù)目至多Nmax個(gè)為止，這樣就得到了罩住所有離散單元的一批立方體，按照相應(yīng)規(guī)則即可建立對應(yīng)八叉樹。

對于某個(gè)參考立方體boxj而言，其他立方體又被區(qū)分為近距立方體和遠(yuǎn)距立方體。每個(gè)立方體內(nèi)及近距立方體內(nèi)離散單元間的相互作用仍采用磁場積分法直接進(jìn)行計(jì)算，而大量遠(yuǎn)距立方體內(nèi)離散單元間的作用通過聚合、轉(zhuǎn)移、發(fā)散的快速多極算法來實(shí)現(xiàn)。因而，這就避免了顯式地產(chǎn)生積分方程系數(shù)矩陣，且加速了迭代法求解代數(shù)方程中矩陣向量乘積運(yùn)算。

圖2 用四叉樹表示的立方體關(guān)系示意圖Fig.2 The relationship ofcubicsdipicted with quadtree

2.1 積分核函數(shù)多極展開

為適應(yīng)三維復(fù)雜艦船的剖分，采用非規(guī)則六面體單元。當(dāng)剖分單元數(shù)目足夠多時(shí)，可將每個(gè)離散單元磁場強(qiáng)度矢量B視為常量，被積函數(shù)中各物理量之間關(guān)系如圖3所示。

圖3 被積函數(shù)中各物理量之間關(guān)系圖Fig.3 The relationship of each symbol in integral function

為便于加速計(jì)算，根據(jù)快速多極子理論，需將式(2)被積函數(shù)中Green函數(shù)的基本解用雙諧函數(shù)進(jìn)行展開［14］:

式中:rP為場點(diǎn)P的矢徑，rQ為源點(diǎn)Q的矢徑，n為多極子展開階數(shù)，取值影響著計(jì)算精度和計(jì)算速度; Rn，m和Sn，m分別為球面調(diào)和函數(shù)表示Sn，m的共軛復(fù)數(shù)，不同應(yīng)用領(lǐng)域其定義略有不同，此處定義為

式中:(ρ，α，β)和(r，θ，φ)分別為rQ和rP的球坐標(biāo)，為連帶勒讓德函數(shù)，其定義為

其中，Pn為勒讓德函數(shù)，且連帶勒讓德函數(shù)滿足:

將式(2)采用式(5)展開，經(jīng)推導(dǎo)不難得到

其中

式(13)稱為以O(shè)點(diǎn)為中心的源點(diǎn)矩或多極子展開公式(multipole expansion，ME)。

2.2 源點(diǎn)矩到本地展開系數(shù)的轉(zhuǎn)變

當(dāng)|OPO|＞|POP|，有球面調(diào)和函數(shù)關(guān)系式［15］:

代入以O(shè)為中心的展開式(12)，經(jīng)推導(dǎo)可得:

其中

式(15)稱為將以O(shè)為中心的源點(diǎn)矩轉(zhuǎn)變?yōu)橐訮O為中心的展開系數(shù)的傳遞公式(moment to local translation，M2L)。

2.3 本地展開系數(shù)轉(zhuǎn)移

基于式(14)可將本地展開系數(shù)發(fā)散至場點(diǎn)P (local expansion，LE)，實(shí)現(xiàn)相應(yīng)單元之間貢獻(xiàn)量的計(jì)算，簡寫成:

在利用快速多極子加速磁場積分法計(jì)算中，源區(qū)單元對場區(qū)單元的聚合、轉(zhuǎn)移和發(fā)散等過程如圖4所示。

圖4 聚合、轉(zhuǎn)移、發(fā)散等過程示意圖Fig.4 The sketch map of ME，M2L and LE

2.4 算法基本步驟

采用快速多極子加速積分方程求解三維靜磁問題具體步驟如下:

1)網(wǎng)格剖分。對計(jì)算對象進(jìn)行離散剖分，以獲取單元信息和節(jié)點(diǎn)信息。目前可用來進(jìn)行網(wǎng)格剖分的軟件較多，為提高效率，本文不獨(dú)立編寫剖分程序而直接選用TrueGrid軟件來進(jìn)行剖分。鑒于非規(guī)則六面體單元具有剖分通用性強(qiáng)等優(yōu)點(diǎn)，因而將其作為計(jì)算對象網(wǎng)格模型的基本單元。

2)建立八叉樹。將整個(gè)計(jì)算對象用一個(gè)足夠大的立方體罩住;將該立方體分為8個(gè)子立方體，每個(gè)立方體再分為8個(gè)更小的子立方體，依次類推直到所需層數(shù)，該層滿足每個(gè)子立方體包含的單元數(shù)不超過給定的數(shù)目Nmax。按照要求將該層立方體編號并建立描述立方體相互關(guān)系的近距組、遠(yuǎn)距組信息數(shù)組［14］，如第i個(gè)立方體內(nèi)部單元數(shù)numt(i)、第i個(gè)立方體父立方體的編號ifath(i)、第i個(gè)立方體在本層中的編碼itree(i)、第i個(gè)立方體開始的單元個(gè)數(shù)編號loct(i)、單元原始編號在八叉樹結(jié)構(gòu)中的編碼ielem(i)等。

3)采用Gmres迭代法求解方程組(4)。迭代過程中涉及大量C3N×3NX3N矩陣向量乘積運(yùn)算，但采用快速多極子與積分方程相耦合可大大加速矩陣向量乘積運(yùn)算，具體步驟如下:

①初始化迭代參數(shù)、快速多極子參數(shù);

②將各單元作用聚合到每個(gè)相應(yīng)立方體中心，即按式(13)計(jì)算源點(diǎn)矩;

③按照原始積分式(3)直接計(jì)算自身立方體內(nèi)及近距立方體內(nèi)單元之間的相互作用;

④采用快速多極算法計(jì)算遠(yuǎn)距立方體內(nèi)單元之間的相互作用，以boxj為例，根據(jù)式(15)將其任一遠(yuǎn)距立方體boxi的源點(diǎn)矩轉(zhuǎn)移到自身立方體boxj中心，根據(jù)式(16)將以boxj為中心的展開系數(shù)發(fā)散至自身立方體內(nèi)每一個(gè)單元中心，通過這一過程即實(shí)現(xiàn)了boxi內(nèi)所有單元到boxj內(nèi)所有單元貢獻(xiàn)系數(shù)的計(jì)算;

⑤將上述2部分貢獻(xiàn)量相加即可實(shí)現(xiàn)式(4)中的矩陣矢量相乘，可以寫為

式中:w1為包含在葉子立方體內(nèi)自身及其近距立方體內(nèi)的源點(diǎn)單元，w2為包含在遠(yuǎn)距立方體內(nèi)源點(diǎn)單元。等式右邊第1項(xiàng)表示w1中的單元對第j個(gè)單元貢獻(xiàn)量，采用磁場積分法直接計(jì)算;等式右邊第2項(xiàng)表示w2中的單元對第j個(gè)單元貢獻(xiàn)量，采用快速多極算法加速計(jì)算。

⑥根據(jù)Gmres迭代算法更新每個(gè)單元內(nèi)部磁感應(yīng)強(qiáng)度，判別誤差要求，若不滿足迭代結(jié)束條件，返回②重新迭代計(jì)算。

2.5 編程流程及數(shù)值實(shí)現(xiàn)的注意點(diǎn)

為提高代碼執(zhí)行效率，基于快速多極子和磁場積分法的船磁計(jì)算軟件源代碼采用Fortran90編寫，在Intel Visual Fortran編譯環(huán)境下調(diào)試運(yùn)行。計(jì)算硬件選用8核3.4GHz CPU和4GB RAM臺式計(jì)算機(jī)，操作系統(tǒng)為32位Windows XP。算法實(shí)現(xiàn)流程如圖5所示。

圖5 快速多極子加速積分方程計(jì)算流程圖Fig.5 The flow chart of fast multipole method and integral equation method

3 計(jì)算實(shí)例

3.1 均勻外場作用下鐵質(zhì)球體感應(yīng)磁場解析算例

采用均勻外磁場作用下鐵質(zhì)球體產(chǎn)生的感應(yīng)磁場，驗(yàn)證所提出算法及程序的正確性。

如圖6所示，半徑為10 m、相對磁導(dǎo)率為150的鐵質(zhì)球體放置于均勻外磁場B=［34 500 0 0］中。鐵質(zhì)球體受外磁場磁化在其周圍產(chǎn)生的感應(yīng)磁場可解析求得［16］，將其作為用于數(shù)值模擬比較的理論值。61個(gè)計(jì)算場點(diǎn)沿x軸方向關(guān)于原點(diǎn)對稱均勻分布在球體下方、間距0.5 m、距球體中心高度15 m。

圖6 均勻地磁場作用下鐵質(zhì)球體剖分Fig.6 The mesh of the iron sphere submerged in an applied magnetic field

為考核基于快速多極子和磁場積分法的船磁快速算法的計(jì)算效率及所編程序正確性，分別將該鐵質(zhì)球體剖分為729、1 331、1 728、2 197、3 375個(gè)非規(guī)則六面體單元。在上述5種剖分條件下，分別采用直接磁場積分法(DIEM)［16］和采用快速多極子加速積分方程快速算法(本文方法)來計(jì)算場點(diǎn)處感應(yīng)磁場。本文方法算法中單個(gè)立方體包含的最大單元數(shù)Nmax取，其中EleNum為單元個(gè)數(shù)。迭代參數(shù)設(shè)置為:迭代結(jié)束相對殘差取10－8，最大迭代次數(shù)50。

圖7給出了采用2種數(shù)值模擬方法所求得感應(yīng)磁場與理論值的對比曲線，其中Theory、DIEM、本文方法分別表示由解析方法、直接磁場積分法和采用快速多極子加速積分方程求得的感應(yīng)磁場，Bx、By和Bz分別表示球體感應(yīng)磁場的X分量、Y分量和Z分量。由圖7結(jié)果表明，與感應(yīng)磁場理論值相比，本文方法和DIEM計(jì)算的最大相對誤差分別為0.63%和0.49%，均可實(shí)現(xiàn)高精度計(jì)算。圖中局部放大處表明本文方法的計(jì)算誤差略大于DIEM，主要是由于遠(yuǎn)區(qū)單元之間作用近似計(jì)算所致。圖8給出了2種數(shù)值模擬方法求解場點(diǎn)感應(yīng)磁場計(jì)算時(shí)間比較曲線。

圖7 球體感應(yīng)磁場計(jì)算值與理論值比較圖Fig.7 The comparison of the calculated and the theory induced magnetic field of the iron sphere

圖8 2種數(shù)值模擬方法計(jì)算時(shí)間比較曲線Fig.8 The comparison of the computation time of the two numerical methods

由圖8可以看出，本文方法計(jì)算效率遠(yuǎn)高于DIEM，隨著剖分單元個(gè)數(shù)的增加，本文方法方法計(jì)算效率優(yōu)勢越來越明顯。上述計(jì)算結(jié)果表明了采用快速多極子加速積分方程快速算法所編程序正確性。

3.2 簡化艇模主艙段感應(yīng)磁場實(shí)驗(yàn)算例

為方便，采用地磁場作用下簡化艇模主艙段產(chǎn)生的感應(yīng)磁場，驗(yàn)證所提出算法的工程實(shí)用性。

如圖9所示，相對磁導(dǎo)率155的簡化艇模主艙段，置于均勻地磁場中，地磁場水平分量和垂直分量分別為34 500 nT和34 800 nT。計(jì)算場點(diǎn)沿z軸方向關(guān)于原點(diǎn)對稱均勻分布于簡化艇模主艙段下方，縱向間距100 mm、橫向間距200 mm、垂向距離358 mm。

圖9 地磁場作用下簡化艇模主艙段剖分Fig.9 The mesh of the main cabin of the simple submarine mockup submerged in an applied magnetic field

為消除簡化艇模主艙段剩磁對考核數(shù)值模擬誤差的影響，利用測磁精度為1 nT的多通道測磁系統(tǒng)在地磁北向和南向分別測量得到了簡化艇模主艙段產(chǎn)生的磁感應(yīng)強(qiáng)度，并根據(jù)所測得的數(shù)據(jù)分解出簡化艇模主艙段受地磁場水平磁化作用而在場點(diǎn)產(chǎn)生的磁感應(yīng)強(qiáng)度值［2］，將其作為數(shù)值模擬的比對值。

將簡化艇模主艙段剖分為8 160個(gè)非規(guī)則六面體單元，單個(gè)立方體包含的最大單元數(shù) Nmax取迭代參數(shù)設(shè)置與計(jì)算實(shí)例1相同。

采用本文方法來計(jì)算簡化艇模主艙段感應(yīng)磁場共耗時(shí)214 s。圖10給出了簡化艇模主艙段受地磁場水平分量磁化作用而產(chǎn)生的感應(yīng)磁場測量值與計(jì)算值對比曲線。

圖10 簡化艇模主艙段受地磁場水平分量磁化作用而產(chǎn)生的感應(yīng)磁場測量值與計(jì)算值對比曲線圖Fig.10 The comparison of the calculated and the measured induced magnetic field of the simple submarine mockup magneitzed by horizontal geomagnetic field

由圖10可以看出，測量值與計(jì)算值吻合很好，最大相對誤差為4.14%，可以滿足工程需求。然而，在相同條件下采用直接磁場積分法進(jìn)行計(jì)算時(shí)，耗時(shí)超過1 h。相比較而言，本文方法可大大提高船磁計(jì)算效率。

4 結(jié)束語

本文提出了基于磁場積分法和快速多極子的艦船感應(yīng)磁場快速預(yù)報(bào)方法，并編制了相應(yīng)程序包，為解決艦船感應(yīng)磁場快速計(jì)算提供了一條技術(shù)途徑。解析球體和簡化艇模主艙段感應(yīng)磁場計(jì)算實(shí)例表明，該方法計(jì)算精度高、速度快，具有很高的工程實(shí)用價(jià)值。與直接磁場積分法相比，在相同計(jì)算條件下，當(dāng)離散單元數(shù)目達(dá)到3 000時(shí)，即可減少計(jì)算時(shí)間近5倍，隨著離散單元數(shù)目的增加，該方法加速效果更加顯著。值得指出的是，文中采用的是單層快速多極子來加速磁場積分船磁計(jì)算過程，若采用多層快速多極子、自適應(yīng)多層快速多極子等技術(shù)效果將更加明顯，這也是將來進(jìn)一步研究的內(nèi)容。

［1］HOLMES J J.Modeling a ship＇s ferromagnetic signatures［M］.［S.l.］:Morgan and Claypool Publishers，2007:1-4.

［2］郭成豹，劉大明，肖昌漢，等.艦載活動(dòng)設(shè)備磁場建模分析研究［J］.兵工學(xué)報(bào)，2009，30(2):196-200.

GUO Chengbao，LIU Daming，XIAO Changhan，et al.Magnetic field modelling of shipborne moving equipment［J］.Acta Armamentarii，2009，30(2):196-200.

［3］翁行泰，曹梅芬.潛艇感應(yīng)磁場的三維有限元法計(jì)算研究［J］.上海交通大學(xué)學(xué)報(bào)，1994，28(5):69-76.

WENG Xingtai，CAO Meifen.A research on the calculation of the induced magnetic field of a submarine［J］.Journal of Shanghai Jiaotong University，1994，28(5):69-76.

［4］陳杰，魯習(xí)文.艦船感應(yīng)磁場預(yù)測的一種新方法［J］.物理學(xué)報(bào)，2010，59(1):239-245.

CHEN Jie，LU Xiwen.A new method for predicting the induced magnetic field of naval vessels［J］.Acta Physica Sinica，2010，59(1):239-245.

［5］CHADEBEC O，COULOMB J L，LECONTE V，et al.Modeling of static magnetic anomaly created by iron plates［J］.IEEE Transactions on Magnetics，2000，36(4):667-671.

［6］郭成豹，何明，周耀忠.積分方程法計(jì)算艦船感應(yīng)磁場［J］.海軍工程大學(xué)學(xué)報(bào)，2001，13(6):71-74.

GUO Chengbao，HE Ming，ZHOU Yaozhong.Calculation of induced magnetic fields of ships by integral equation method［J］.Journal of Naval University of Engineering，2001，13 (6):71-74.

［7］周國華，肖昌漢，閆輝，等.一種弱磁作用下鐵磁物體感應(yīng)磁場的計(jì)算方法［J］.哈爾濱工程大學(xué)學(xué)報(bào)，2009，30 (1):91-95.

ZHOU Guohua，XIAO Changhan，YAN Hui，et al.A method to calculate the induced magnetic field of ferromagnetic objects in a weak magnetic field［J］.Journal of Harbin Engineering University，2009，30(1):91-95.

［8］ROKHLIN V.Rapid solution of integral equations of classical potential theory［J］.Journal of Computational Physics，1985，60(2):187-207.

［9］GREENGARD L，ROKHLIN V.A fast algorithm for particle simulations［J］.Journal of Computational Physics，1987，73(2):325-348.

［10］聶在平，胡俊，姚海英，等.用于復(fù)雜目標(biāo)三維矢量散射分析的快速多極子方法［J］.電子學(xué)報(bào)，1999，27 (6):104-109.

NIE Zaiping，HU Jun，YAO Haiying，et al.The fast multipole methods for vector analysis of scattering from 3-dimensional objects with complex structure［J］.Acta Electronica Sinica，1999，27(6):104-109.

［11］寧德志，滕斌，勾瑩.快速多極子展開技術(shù)在高階邊界元方法中的實(shí)現(xiàn)［J］.計(jì)算力學(xué)學(xué)報(bào)，2005，22(6): 700-704.

NING Dezhi，TENG Bin，GOU Ying.Implementation of the fast multipole expansion technique in the higher order BEM［J］.Chinese Journal of Computational Mechanics，2005，22(6):700-704.

［12］吳海軍，蔣偉康，魯文波.三維聲學(xué)多層快速多極子邊界元及其應(yīng)用［J］.物理學(xué)報(bào)，2012，61(5):054301.

WU Haijun，JIANG Weikang，LU Wenbo.Multilevel fast multipole boundary element method for 3D acoustic problems and its applications［J］.Acta Physica Sinica，2012，61(5):054301.

［13］樊明武，顏威利.電磁場積分方程法［M］.北京:機(jī)械工業(yè)出版社，1988:11-19.

［14］LIU Yijun.Fast multipole boundary element method theory and applications in engineering［M］.New York:Cambridge University Press，2009:71-73.

［15］YOSHIDA K.Applications of fast multipole method to boundary integral equation method［M］.［S.l.］:Kyoto University Press，2001:106-108.

［16］周國華，肖昌漢，劉勝道，等.基于六面體單元表面磁場積分法求解三維靜磁場［J］.電工技術(shù)學(xué)報(bào)，2009，24(3):1-7.

ZHOU Guohua，XIAO Changhan，LIU Shengdao，et al.3D magnetostatic field computation with hexahedral surface integral equation method［J］.Transactions of China Electrotechnical Society，2009，24(3):1-7.

［17］盛新慶.計(jì)算電磁學(xué)要論［M］.北京:科學(xué)出版社，2004:29-31.