一種高效率的RPIM法及其在電磁場中的應用

王立鵬，王欣彥，戰洪仁，張先珍，寇麗萍

(沈陽化工大學，遼寧沈陽110142)

0 引 言

近些年來發展的無單元法不需要生成網格，通過節點來構造插值函數，在國內外受到高度重視，成為熱門的數值計算方法。無單元法在分析裂紋擴展、特大變形、碰撞、金屬沖壓成型等很多問題具有優勢［1］。無單元法也可被用來解決電磁場數值計算網格限制的問題［2-7］。目前，已有多種無單元法，徑向基點插值型無單元法［8］是其中一種。徑向基函數(RBF，Radial Basis Fuction)是以點x到節點xi之間的距離為自變量構成的一類函數，也叫距離基函數。它具有求解效率高、運算簡單、各向同性和空間維數無關等優點，幾乎可以逼近所有的函數。因此，在多元逼近論中用徑向基函數進行插值已成為一種有效的工具。常用的徑向基函數有數種，其中的復合2次Multiquadrics(MQ)徑向基函數由于具有較高的收斂率［9］，應用較為廣泛。

然而，無單元法的計算量較大，其求解效率不及有限元法。通常，無單元法的求解采用迭代方法，包括高斯-賽德爾(Gauss-Seidel)迭代、松弛迭代、共軛梯度法(Conjugate Gradient，CG)、廣義最小余數法等。

目前，在迭代求解方法中，多重網格法是一種比較熱門的求解方法。簡單地說，多重網格法是采用不同尺度的網格，先在細網格上進行光滑迭代，消除高頻誤差，然后通過限制算子R將細網格上殘差轉換至粗網格，再通過粗網格來消除低頻誤差，這樣一層一層傳下去，然后殘差又通過延拓算子R一層一層返回到細網格上，對先前的細網格迭代解進行修正。這就是應用較為廣泛的V循環多重網格法的思想。

上世紀60年代，多重網格方法(MG)開始出現，起初并未受到重視。從20世紀70年代開始，多重網格法成為研究熱點，得到了快速的發展，現在已成為一種有效的數值方法，用來求解偏微分方程組離散化形成的大型代數方程組［10-18］。多重網格法的主要優點是它具有快速收斂性，收斂速度與網格尺寸無關。該方法可應用到電磁場的數值分析中［19］。

多重網格法是以有限差分法為架構來引入多重網格的校正與迭代技術，實質上是多重網格有限差分法。因而，同有限差分法一樣，多重網格法對復雜的定解條件、復雜的邊界幾何結構等方面問題模擬較為困難。

多重網格法求解技術同有限單元法相結合［20］，可提高求解效率，正被用來求解各種數值計算問題。目前，該方法在無單元法的應用還很有限，多重網格法已被用在無單元伽遼金法上來求解電磁場問題［21］，但這種高效的求解方法還未被應用到徑向基點插值型無單元法。本文將徑向基點插值型無單元法同多重網格法思想相結合，提出基于徑向基點插值型無單元的多重節點法，該方法在粗節點的構造方法及限制算子的確定等方面都與有限元多重網格法有所區別。通過算例分析，證明多重節點法可獲得較為理想的徑向基點插值型無單元法的求解效率。

1 徑向基點插值型無單元多重節點法(RPIM－MN)原理

1.1 多重節點法(MN)

多重網格法是基于網格基礎之上的，在有限單元法、有限差分法等以網格為基礎的離散方法中都有廣泛的應用。而無單元法本身指的是通過一組散布節點來構造場函數的近似表達式，可以消除網格，即使有背景網格也只是用來進行積分計算的。可見，以網格為基礎的多重網格法并不是為無單元法量身定做的。但是，可以借鑒高效求解的多重網格法的思想，將其引入到無單元法的求解中，這就是多重節點法。

為了簡單起見，只分析用粗細兩層無單元節點。在細節點的離散場域上，實施徑向基點插值型無單元法求解給定的線性方程組:

式中:KH為粗節點線性系統的系數矩陣;ΩH為粗節點求解域;rH為殘差。

該方法的計算步驟如下:

(1)前光滑迭代計算，對KHxH=bH進行細節點迭代，得到近似的解υH，并計算細節點近似解的殘差。

(2)對殘差進行限制，rH←R(bυh-Khυh，由 Ωh向ΩH轉移。

(3)在ΩH上解殘差方程。設υH=0，解KHxH=rH，用υH作為近似的解。

粗節點系數矩陣:KH=RKhP

(4)插值，由向ΩH向Ωh轉移校正量，并修正Ωh上的解:

由R=σPT可確定P，其中σ為比例系數。

(5)后光滑迭代計算，ε0為一個任意小的正數，若 ε≥ε0，轉向步驟(1)。

1.2 徑向基點插值型無單元多重節點法(RPIMMN)

在粗節點的構造方法及限制算子的確定方法上，徑向基點插值型無單元多重節點法同有限元多重網格法有所區別。

(1)粗節點的構造算法

該法采用了聚集式的方法構造粗節點，也就是從初始細節點中，按一定地方式選出粗節點，不再另外布置新點。具體如下:

(a)為每個節點k確定一個“影響值”v(k)，其值等于節點k的影響域內節點數目。

(b)對于邊界節點，首先將邊界幾何結構突變點選作粗節點;然后按“影響值”從大到小選取N個節點為粗節點。N約為邊界節點總數的一半。

(c)對于內部節點，按“影響值”從大到小選取M個節點為粗節點。M約為節點總數的一半。

(d)對于均勻布點，可在細節點的基礎上間隔的選取粗節點。

(e)將新粗節點的影響域中所有未定點變為細節點。

圖1表示了均勻布點情況下，兩重節點法的粗、細節點布置及各節點的支撐域。圖中的圓形區域為各節點的影響域。

圖1 徑向基點插值型無單元多重節點法示意圖

(2)限制算子R的確定

對無單元法來說，限制算子R是指將細節點上殘差限制到粗節點上的映射矩陣。限制算子是根據粗、細節點的位置來確定的。由于沒有網格的束縛，徑向基點插值型無單元兩重節點法粗節點和細節點之間的映射關系較為簡單，可采用節點間的距離函數來構造。本文討論兩重節點(細節點層為k，粗節點層為k-1)的迭代情況。

限制算子矩陣元素:

式中:i為任意一個粗節點號;j為任意一個在粗節點i影響域內的無單元節點號(不包括i節點);m是粗節點i影響域內節點數;rij表示粗節點i與節點的距離。

按式(3)，限制算子R可以根據細節點和粗節點在場域中的相應位置來確定。

延拓算子P可根據P=RT來確定。

上述算法也可以類推，形成徑向基點插值型無單元多重節點法。

2 計算實例

算例一:一長直接地金屬槽，a=0.1 m，槽內電位函數 φ(x，y)滿足 Laplace方程，φ0=10 V，邊界條件如圖2所示。取第一類邊界上罰因子為10 000;在計算中采用了二重節點，51×51個細節點，26×26個粗節點，計算后再增加節點數，比較多重節點法(MN)與共軛梯度法(ICCG)的計算結果，如圖3所示。可以看出，同 ICCG法相比，MN法的計算時間要少得多。

算例二:一臺長為55 mm、外徑為53 mm的永磁直流電動機，該電機為瓦片形磁極結構，4對極8槽，釹鐵硼永磁體磁化方向為平行磁化。其磁鋼矯頑力為7.81×105A/m，所建立的場域模型如圖4所示。

圖4 永磁直流電動機空載磁場

對電機的空載電磁場進行分析，取一個極距為求解區域，AB邊、CD邊滿足第二類齊次邊界條件:

用無單元法離散電機電磁場域，分別采用MN法和ICCG法對其進行求解。MN法采用二重節點法，無單元細節點數為59 731，無單元粗節點數為29460。

通過計算，得到10-4精度要求下的電機磁力線分布，如圖4所示。其中，二重節點法計算時間為187 s，ICCG法計算時間為294 s。

由于多重節點法的迭代計算是分散在各層中進行的，且粗節點數目較少，所以計算迭代次數時，粗節點的迭代次數必須折合到細節點，方便與單層ICCG法的迭代次數進行對比。

將W定義為總迭代數，其表達式:

式中:M為多重節點的層數;Nk為在第k層中的迭代次數。

由于在多層節點中，W相當于把各層的迭代計算次數都折算到最細層，然后疊加。由于忽略了限制算子和延拓算子的工作量，總迭代數可方便、粗略地估計多重節點法計算量。

比較多重節點法總迭代數與ICCG法的總迭代數，得到總迭代數與收斂容差之間的關系曲線，如圖5所示。容易看出，在滿足相同的精度條件下，多重節點法的總迭代數要比ICCG法的小很多。也就是說，徑向基點插值型無單元多重節點法較ICCG法有更高的計算效率。

圖5 總迭代數和收斂容差關系曲線

3 結 論

本文對于用多重節點法求解RPIM無單元法離散的電磁場域進行研究，總結如下:

(1)通過引入多重網格法的思想，構造了適用于徑向基點插值型無單元法的多重節點法，并實現了用多重節點法求解徑向基點插值型無單元法問題，從而在快速收斂的同時又兼有無單元法消除網格限制的特點。

(2)對同一問題來說，在相同的精度要求下，在同一自由度下，RPIM-MN法要較ICCG法具有更高的計算效率。

［1］ 張雄，宋康祖，陸明萬.無網格法研究進展及其應用［J］.計算力學學報，2003，20(6):730-742.

［2］ Maréchal Y.Some meshless methods for electromagnetic field computations［J］.IEEE Trans.Mag.，1998，34(5):3351-3354.

［3］ Cingoski V，Miyamoto N，Yamashita H.Element-free Galerkin method for electromagnetic field Computations［J］.IEEE Trans.Mag.，1998，34(5):3236-3239.

［4］ 劉素貞，楊慶新，陳海燕，等.無單元法在電磁場數值計算中的應用研究［J］.電工技術學報.2001，16(2):30-33.

［5］ Lu Huifen，Jiang Hao.Application of meshless local Petrov-Galerkin methods in electromagnetics［C］//Electrical Machines and Systems.2003，2:798-801.

［6］ Ho S L，Yang S Y，Wong H C，et al.A meshless collocation method based on radial basis functions and wavelets［J］.IEEE Trans.Magn.，2004，40(2):1021-1024.

［7］ Li Gang，Aluru N R.Linear，nonlinear and mixed-regime analysis of electrostatic MEMS［J］.Sensor and Actuators(A)，2001，91(3):278-291.

［8］ Wang J G，Liu G R.A Point Interpolation Meshless Method Based on Radial Basis Functions［J］.Int.J.Numer.Methods Engrg.，2002，54:1623-1648.

［9］ Schaback R.Error estimates and condition numbers for radial basis function interpolation［J］.In Proceedings of Adv.Comput.Math.，1995:251-264.

［10］ Hackbusch W，Trottenberg U.Multigrid methods III［C］//Proc.3rd Int.Conf.on Multigrid methods.Int.Series of Num Math 98.Birkhauser，Basel 1991.

［11］ Hackbusch W，Trottenberg U.Multigrid methods［C］//Lecture Notes in Mathematics 960.Springer，Berlin，1986.

［12］ Iserles A.A first course in the numerical analysis of differential equations［M］.Cambridge:Cambridge Univ.Press，1996.

［13］ Hackbusch W.Convergence of multi-grid iterations applied to difference equation［J］.Math.Comp.，1980，34:425-440.

［14］ Hackbusch W.On the convergence of multi-grid iterations［J］.Beitrage Nuimer.Muth.，1981，9:213-239.

［15］ Hackbusch W，Wittum G.Multigrid methods V［C］//Proc.5th European Multigrid Conf.in Stuttgart Germany，Lecture Notes in Computational Science and Engineering 3.Springer，Berlin，1998.

［16］ Hackbusch W.Multi-grid methods and applications［M］.New York:Springer-Verlag，1985.

［17］ Hackbusch W，Trottenberg U.Multigrid methods II［C］//Lecture Notes in Mathematics 1228.Springer，Berlin，1982.

［18］ Hackbusch W.Robust multigrid methods［M］.Friedr Vieweg 8 Sohn，1988.

［19］ 周浩，倪光正.多重網格法及其在靜態電磁場數值分析中的應用［J］.中國電機工程學報，1990，10(5):20-26.

［20］ Braess D，Verfurth R.Multigrid methods for nonconforming finite element methods［J］.IAM J.Numer.Anal.，1990，27:979-986.

［21］ 王立鵬，王欣彥，唐任遠.EFG-MG法及其在電磁場數值計算中的應用［J］.電工技術學報.2010，25(1):1-5.