基于快速多極子的艦船感應磁場快速預報方法

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

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

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

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

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

1 磁場積分法

如圖1所示，在外磁場B0作用下的鐵磁物體(如鐵質艦船)在空間點P產生的磁感應強度可表示為［13］

式中:▽P為對場點坐標的梯度算子，▽Q為對源點坐標的梯度算子，v為鐵磁物體所占體積，M為由外磁場磁化引起的鐵磁物體內部磁化強度，P為場點，Q為源點。

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

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

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

為得到鐵磁物體內部的磁化強度，通常將鐵磁物體離散為若干足夠小的均勻磁化體。以每個離散單元中心為計算場點，可建立以單元內部磁感應強度為未知量的代數方程組

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

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

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

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

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

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

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

2.1 積分核函數多極展開

為適應三維復雜艦船的剖分，采用非規則六面體單元。當剖分單元數目足夠多時，可將每個離散單元磁場強度矢量B視為常量，被積函數中各物理量之間關系如圖3所示。

圖3 被積函數中各物理量之間關系圖Fig.3 The relationship of each symbol in integral function

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

式中:rP為場點P的矢徑，rQ為源點Q的矢徑，n為多極子展開階數，取值影響著計算精度和計算速度; Rn，m和Sn，m分別為球面調和函數表示Sn，m的共軛復數，不同應用領域其定義略有不同，此處定義為

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

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

將式(2)采用式(5)展開，經推導不難得到

其中

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

2.2 源點矩到本地展開系數的轉變

當|OPO|＞|POP|，有球面調和函數關系式［15］:

代入以O為中心的展開式(12)，經推導可得:

其中

式(15)稱為將以O為中心的源點矩轉變為以PO為中心的展開系數的傳遞公式(moment to local translation，M2L)。

2.3 本地展開系數轉移

基于式(14)可將本地展開系數發散至場點P (local expansion，LE)，實現相應單元之間貢獻量的計算，簡寫成:

在利用快速多極子加速磁場積分法計算中，源區單元對場區單元的聚合、轉移和發散等過程如圖4所示。

圖4 聚合、轉移、發散等過程示意圖Fig.4 The sketch map of ME，M2L and LE

2.4 算法基本步驟

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

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

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

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

①初始化迭代參數、快速多極子參數;

②將各單元作用聚合到每個相應立方體中心，即按式(13)計算源點矩;

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

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

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

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

⑥根據Gmres迭代算法更新每個單元內部磁感應強度，判別誤差要求，若不滿足迭代結束條件，返回②重新迭代計算。

2.5 編程流程及數值實現的注意點

為提高代碼執行效率，基于快速多極子和磁場積分法的船磁計算軟件源代碼采用Fortran90編寫，在Intel Visual Fortran編譯環境下調試運行。計算硬件選用8核3.4GHz CPU和4GB RAM臺式計算機，操作系統為32位Windows XP。算法實現流程如圖5所示。

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

3 計算實例

3.1 均勻外場作用下鐵質球體感應磁場解析算例

采用均勻外磁場作用下鐵質球體產生的感應磁場，驗證所提出算法及程序的正確性。

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

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

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

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

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

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

由圖8可以看出，本文方法計算效率遠高于DIEM，隨著剖分單元個數的增加，本文方法方法計算效率優勢越來越明顯。上述計算結果表明了采用快速多極子加速積分方程快速算法所編程序正確性。

3.2 簡化艇模主艙段感應磁場實驗算例

為方便，采用地磁場作用下簡化艇模主艙段產生的感應磁場，驗證所提出算法的工程實用性。

如圖9所示，相對磁導率155的簡化艇模主艙段，置于均勻地磁場中，地磁場水平分量和垂直分量分別為34 500 nT和34 800 nT。計算場點沿z軸方向關于原點對稱均勻分布于簡化艇模主艙段下方，縱向間距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

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

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

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

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

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

4 結束語

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

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

［2］郭成豹，劉大明，肖昌漢，等.艦載活動設備磁場建模分析研究［J］.兵工學報，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］翁行泰，曹梅芬.潛艇感應磁場的三維有限元法計算研究［J］.上海交通大學學報，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］陳杰，魯習文.艦船感應磁場預測的一種新方法［J］.物理學報，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］.海軍工程大學學報，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］周國華，肖昌漢，閆輝，等.一種弱磁作用下鐵磁物體感應磁場的計算方法［J］.哈爾濱工程大學學報，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］聶在平，胡俊，姚海英，等.用于復雜目標三維矢量散射分析的快速多極子方法［J］.電子學報，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］寧德志，滕斌，勾瑩.快速多極子展開技術在高階邊界元方法中的實現［J］.計算力學學報，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］吳海軍，蔣偉康，魯文波.三維聲學多層快速多極子邊界元及其應用［J］.物理學報，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］.北京:機械工業出版社，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］.電工技術學報，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］盛新慶.計算電磁學要論［M］.北京:科學出版社，2004:29-31.