王占輝 ，高俊吉

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

0 引言

電磁場(chǎng)開(kāi)域問(wèn)題是電氣工程與電磁散射和輻射領(lǐng)域中常遇到的一類(lèi)問(wèn)題。有限元法目前在電磁場(chǎng)分析計(jì)算中居主導(dǎo)地位，但從原理而言，它不適于解決開(kāi)域問(wèn)題。邊界元法可以解決開(kāi)域問(wèn)題，但處理非線性問(wèn)題及多種介質(zhì)問(wèn)題比較困難。如果在鐵磁區(qū)采用有限元進(jìn)行離散，而在鐵磁邊界上采用邊界元法進(jìn)行離散，然后把兩種方法在邊界上有機(jī)結(jié)合起來(lái)，即可完成對(duì)無(wú)界域的求解。這就是有限元－邊界元混合方法。該方法綜合了有限元法與邊界元法的優(yōu)點(diǎn)，成為目前解決復(fù)雜開(kāi)域問(wèn)題最有效且最受人們重視的方法之一。

應(yīng)用混合有限元邊界元法來(lái)求解開(kāi)域問(wèn)題的基本思想[1]就是將整個(gè)求解區(qū)域（即無(wú)限空間）劃分為有限元區(qū)域與邊界元區(qū)域，有限元區(qū)域含源，含多種線性或非線性介質(zhì)；邊界元區(qū)域則為外部的無(wú)限自由空間。在有限元區(qū)域內(nèi)采用有限單元剖分。對(duì)于兩區(qū)域的相交邊界采用邊界單元剖分。最后通過(guò)兩區(qū)域相交邊界上的交界條件，將有限元方程與邊界元方程結(jié)合起來(lái)，得到“混合方程”。求解混合方程確定兩區(qū)域相交邊界上的信息后，即可確定邊界元區(qū)域即外部無(wú)限自由空間中任一點(diǎn)的磁位值或磁場(chǎng)值。

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

1 雙標(biāo)量位的定義[2]

圖1所示為開(kāi)域靜磁場(chǎng)問(wèn)題，區(qū)域A為鐵磁區(qū)域（區(qū)域內(nèi)不存在傳導(dǎo)電流），B為電流區(qū)域（磁導(dǎo)率μ=μ0為常數(shù)），C為導(dǎo)體及空氣區(qū)。區(qū)域A和區(qū)域C兩種介質(zhì)的交界面為Γjk。

圖1 開(kāi)域靜磁場(chǎng)示意圖

1.1 全標(biāo)量位

無(wú)電流的鐵磁區(qū)域A為無(wú)旋場(chǎng)，其場(chǎng)強(qiáng)可以用一個(gè)標(biāo)量位的負(fù)梯度表示，則存在全標(biāo)量磁位設(shè)為ψ，使H=-ψ▽ ，其中ψ稱(chēng)為全標(biāo)量位，它可以用來(lái)表征磁場(chǎng)中總的場(chǎng)強(qiáng)[3]。

利用全標(biāo)量位，則磁感應(yīng)強(qiáng)度可以表示為

當(dāng)磁場(chǎng)區(qū)域存在鐵磁介質(zhì)時(shí)，上式為非線性拉普拉斯方程，它描述了靜磁場(chǎng)中不存在電流區(qū)域的情況。

1.2 部分標(biāo)量位

一般情況下可以認(rèn)為空間任何一點(diǎn)的實(shí)際磁場(chǎng)是由電流源部分所產(chǎn)生的磁場(chǎng)和物質(zhì)被磁化所產(chǎn)生的磁場(chǎng)兩部分的疊加，即（式中為空間任意點(diǎn)的磁場(chǎng)強(qiáng)度；為由宏觀傳導(dǎo)電流在真空介質(zhì)中產(chǎn)生的磁場(chǎng)強(qiáng)度；為鐵磁介質(zhì)被磁化后分子電流磁矩產(chǎn)生的磁場(chǎng)強(qiáng)度）。

對(duì)于電流J產(chǎn)生的磁場(chǎng)強(qiáng)度有；由于，因此，磁化強(qiáng)度的旋度為零，即

則為一無(wú)旋場(chǎng)，可以用一個(gè)標(biāo)量位φ的負(fù)梯度來(lái)表示，即；標(biāo)量位φ部分地描述了磁場(chǎng)的性質(zhì)，它的負(fù)梯度表征了磁化部分的場(chǎng)強(qiáng)，因此稱(chēng)為部分標(biāo)量位。

采用部分標(biāo)量位方法，顯然在計(jì)算中帶來(lái)不少方便，由于每一個(gè)節(jié)點(diǎn)只有一個(gè)未知量，大大降低了對(duì)內(nèi)存的要求。

在有電流的區(qū)域B中，采用部分標(biāo)量位設(shè)為φ，其方程為

2 有限元公式

將鐵磁區(qū)域A進(jìn)行剖分，解函數(shù)ψ用基函數(shù)Ni和節(jié)點(diǎn)函數(shù)值ψi展開(kāi)，即（其中，n為求解區(qū)域中的節(jié)點(diǎn)總數(shù)）

3 邊界元公式

對(duì)于有電流的區(qū)域B中的磁場(chǎng)微分方程(2)，其二維問(wèn)題的格林函數(shù)為；三維問(wèn)題的格林函數(shù)為；r是源點(diǎn)到場(chǎng)點(diǎn)的距離。

則通過(guò)加權(quán)余量法可以得出區(qū)域B的邊界積分方程為[4]

二維情況t=θ2π，θ為點(diǎn)i所張平面角；

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

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

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

4 混合方程

4.1 交界面條件的處理

由電磁學(xué)[3]的知識(shí)可知，在兩種不同介質(zhì)的交界面上，磁感應(yīng)強(qiáng)度的切向分量是連續(xù)的，如果交界面上不存在面電流，則磁場(chǎng)強(qiáng)度的切向分量也是連續(xù)的，即

其中，下標(biāo)j，k表示兩種不同的介質(zhì)；n表示法向分量；t表示切向分量。

由于在鐵磁區(qū)域A中，任一點(diǎn)磁場(chǎng)強(qiáng)度H=-▽?duì)?，而在電流的區(qū)域B中H=-▽?duì)?Hs；（其中Hs為電流產(chǎn)生的磁場(chǎng)），則在兩區(qū)域的交界面Γjk上，根據(jù)方程(7)，可以得到

從而可以導(dǎo)出磁位交界面條件。

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

由于邊界上任一點(diǎn)的全標(biāo)量位與部分標(biāo)量位之差即為地磁場(chǎng)在該點(diǎn)產(chǎn)生的標(biāo)量磁位，則直接可以得到邊界點(diǎn)A處的全標(biāo)量位與部分標(biāo)量位的關(guān)系:

綜上，交界面的條件為

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

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

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

5 仿真算例

對(duì)地磁場(chǎng)中各向同性的鐵質(zhì)均勻無(wú)限長(zhǎng)空心圓柱體進(jìn)行了仿真計(jì)算。

尺寸－外徑：2 m；內(nèi)徑：1 m。地磁場(chǎng)縱向分量27.85（A/m）；橫向、垂向分量為零；圓柱體相對(duì)磁導(dǎo)率：100。

取空心圓柱體的任意一個(gè)與中軸垂直的截面，將其剖分成48個(gè)三角形單元，36個(gè)節(jié)點(diǎn)。

環(huán)外計(jì)算點(diǎn)－縱向：-4~4 m；橫向：5 m；

環(huán)內(nèi)計(jì)算點(diǎn)－縱向：-0.5~0.5 m；橫向：0.6 m。

計(jì)算結(jié)果見(jiàn)圖2～圖3。

計(jì)算誤差－鐵區(qū)誤差：1.03%；內(nèi)部誤差：7.02%；外部誤差：1.23%。（最大值誤差[5]）

從仿真結(jié)果可以看出，計(jì)算結(jié)果比較準(zhǔn)確，說(shuō)明該算法對(duì)開(kāi)域靜磁場(chǎng)問(wèn)題有較好的適用性。

圖2 無(wú)限長(zhǎng)空心圓柱體鐵區(qū)中的磁位值

6 結(jié)束語(yǔ)

本文對(duì)開(kāi)域靜磁場(chǎng)的雙標(biāo)量位混合有限元邊界元方法進(jìn)行了研究。在該方法中，在鐵磁區(qū)域采用全標(biāo)量磁位，而在電流區(qū)和自由空間采用部分標(biāo)量磁位。通過(guò)對(duì)交界面條件的處理，導(dǎo)出全標(biāo)量磁位與部分標(biāo)量磁位在自由空間和鐵磁區(qū)交界面的耦合邊界條件。進(jìn)而推導(dǎo)出雙標(biāo)量位混合有限元與邊界元耦合算法。仿真計(jì)算說(shuō)明本方法對(duì)開(kāi)域靜磁場(chǎng)問(wèn)題有較好的適用性。

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

[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] 顏威利, 楊慶新, 汪友華等.電氣工程電磁場(chǎng)數(shù)值分析.北京：機(jī)械工業(yè)出版社,2005：110～115.

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

[4] 倪光正,楊仕友,錢(qián)秀英等.工程電磁場(chǎng)數(shù)值計(jì)算.北京：機(jī)械工業(yè)出版社, 2004：237～242.

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