一種開域靜磁場雙標量位混合有限元邊界元法研究

王占輝 ，高俊吉

（1. 海軍后勤部后勤辦，北京 100841；2. 海軍工程大學，武漢 430033）

0 引言

電磁場開域問題是電氣工程與電磁散射和輻射領域中常遇到的一類問題。有限元法目前在電磁場分析計算中居主導地位，但從原理而言，它不適于解決開域問題。邊界元法可以解決開域問題，但處理非線性問題及多種介質問題比較困難。如果在鐵磁區采用有限元進行離散，而在鐵磁邊界上采用邊界元法進行離散，然后把兩種方法在邊界上有機結合起來，即可完成對無界域的求解。這就是有限元－邊界元混合方法。該方法綜合了有限元法與邊界元法的優點，成為目前解決復雜開域問題最有效且最受人們重視的方法之一。

應用混合有限元邊界元法來求解開域問題的基本思想[1]就是將整個求解區域（即無限空間）劃分為有限元區域與邊界元區域，有限元區域含源，含多種線性或非線性介質；邊界元區域則為外部的無限自由空間。在有限元區域內采用有限單元剖分。對于兩區域的相交邊界采用邊界單元剖分。最后通過兩區域相交邊界上的交界條件，將有限元方程與邊界元方程結合起來，得到“混合方程”。求解混合方程確定兩區域相交邊界上的信息后，即可確定邊界元區域即外部無限自由空間中任一點的磁位值或磁場值。

本文對全標量位與部分標量位相結合的雙標量位混合有限元邊界元法進行了研究及仿真計算，計算結果表明了該方法的正確性及高精度。

1 雙標量位的定義[2]

圖1所示為開域靜磁場問題，區域A為鐵磁區域（區域內不存在傳導電流），B為電流區域（磁導率μ=μ0為常數），C為導體及空氣區。區域A和區域C兩種介質的交界面為Γjk。

圖1 開域靜磁場示意圖

1.1 全標量位

無電流的鐵磁區域A為無旋場，其場強可以用一個標量位的負梯度表示，則存在全標量磁位設為ψ，使H=-ψ▽ ，其中ψ稱為全標量位，它可以用來表征磁場中總的場強[3]。

利用全標量位，則磁感應強度可以表示為

當磁場區域存在鐵磁介質時，上式為非線性拉普拉斯方程，它描述了靜磁場中不存在電流區域的情況。

1.2 部分標量位

一般情況下可以認為空間任何一點的實際磁場是由電流源部分所產生的磁場和物質被磁化所產生的磁場兩部分的疊加，即（式中為空間任意點的磁場強度；為由宏觀傳導電流在真空介質中產生的磁場強度；為鐵磁介質被磁化后分子電流磁矩產生的磁場強度）。

對于電流J產生的磁場強度有；由于，因此，磁化強度的旋度為零，即

則為一無旋場，可以用一個標量位φ的負梯度來表示，即；標量位φ部分地描述了磁場的性質，它的負梯度表征了磁化部分的場強，因此稱為部分標量位。

采用部分標量位方法，顯然在計算中帶來不少方便，由于每一個節點只有一個未知量，大大降低了對內存的要求。

在有電流的區域B中，采用部分標量位設為φ，其方程為

2 有限元公式

將鐵磁區域A進行剖分，解函數ψ用基函數Ni和節點函數值ψi展開，即（其中，n為求解區域中的節點總數）

3 邊界元公式

對于有電流的區域B中的磁場微分方程(2)，其二維問題的格林函數為；三維問題的格林函數為；r是源點到場點的距離。

則通過加權余量法可以得出區域B的邊界積分方程為[4]

二維情況t=θ2π，θ為點i所張平面角；

三維情況t=Ω/4π，Ω為點i所張立體角。

將交界面為Γjk剖分成有限數量的單元。將解函數φ（?φ/n?）用基函數和節點函數值iφ展開，即（其中，k為單元的插值節點數）。這里插值基函數的選取可采用線性單元，也可用常單元。

本文中φ的基函數采用一次線性單元，?φ/?n的基函數采用常單元。方程離散為

4 混合方程

4.1 交界面條件的處理

由電磁學[3]的知識可知，在兩種不同介質的交界面上，磁感應強度的切向分量是連續的，如果交界面上不存在面電流，則磁場強度的切向分量也是連續的，即

其中，下標j，k表示兩種不同的介質；n表示法向分量；t表示切向分量。

由于在鐵磁區域A中，任一點磁場強度H=-▽ψ，而在電流的區域B中H=-▽φ+Hs；（其中Hs為電流產生的磁場），則在兩區域的交界面Γjk上，根據方程(7)，可以得到

從而可以導出磁位交界面條件。

這里的全標量位ψ的零點選在坐標原點上，部分標量位φ的零點選在無窮遠處。假設地磁場為均勻場，設地磁場向量為則在邊界上任一點A(xi,yi,zi)處，坐標原點到該點的向量，則地磁場在A點的標量磁位可以直接求出，為

由于邊界上任一點的全標量位與部分標量位之差即為地磁場在該點產生的標量磁位，則直接可以得到邊界點A處的全標量位與部分標量位的關系:

綜上，交界面的條件為

4.2 形成有限元－邊界元混合方程

對交界面條件進行處理后，聯立有限元邊界元方程，結合交界面條件(10)可以得到有限元－邊界元混合方程如下：

解之即可得到邊界Γjk上離散點處的φ、?φ/?n；以及區域A中的ψ。利用求得的邊界Γjk上的φ、?φ/?n即可求得空間中任一點的磁場。

5 仿真算例

對地磁場中各向同性的鐵質均勻無限長空心圓柱體進行了仿真計算。

尺寸－外徑：2 m；內徑：1 m。地磁場縱向分量27.85（A/m）；橫向、垂向分量為零；圓柱體相對磁導率：100。

取空心圓柱體的任意一個與中軸垂直的截面，將其剖分成48個三角形單元，36個節點。

環外計算點－縱向：-4~4 m；橫向：5 m；

環內計算點－縱向：-0.5~0.5 m；橫向：0.6 m。

計算結果見圖2～圖3。

計算誤差－鐵區誤差：1.03%；內部誤差：7.02%；外部誤差：1.23%。（最大值誤差[5]）

從仿真結果可以看出，計算結果比較準確，說明該算法對開域靜磁場問題有較好的適用性。

圖2 無限長空心圓柱體鐵區中的磁位值

6 結束語

本文對開域靜磁場的雙標量位混合有限元邊界元方法進行了研究。在該方法中，在鐵磁區域采用全標量磁位，而在電流區和自由空間采用部分標量磁位。通過對交界面條件的處理，導出全標量磁位與部分標量磁位在自由空間和鐵磁區交界面的耦合邊界條件。進而推導出雙標量位混合有限元與邊界元耦合算法。仿真計算說明本方法對開域靜磁場問題有較好的適用性。

混合有限元邊界元方法的主要缺點是所得到的代數方程組的系數矩陣不再對稱且不定，系數矩陣的條件數易變得較大，病態較嚴重，這使得求解比較困難。這就需要采用預條件雙共軛梯度法等求解方法進行處理。這是提高該方法計算效率的有效途徑之一。

[1] S.J.Salon. The hybrid finite element-boundary element method in electromagnetics, IEEE Transac- tions on Magnetics.Vo1.Mag.21, No.5, September, 1986,1829～ 1834.

[2] 顏威利, 楊慶新, 汪友華等.電氣工程電磁場數值分析.北京：機械工業出版社,2005：110～115.

[3] 馮慈璋. 電磁場.北京：高等教育出版社,1983：134-138.

[4] 倪光正,楊仕友,錢秀英等.工程電磁場數值計算.北京：機械工業出版社, 2004：237～242.

[5] 周耀忠,宋武昌,唐申生.潛艇磁場外推的數學模型研究. 海軍工程大學學報, 2003, 15(4)：31～35.