謝裕清，李　琳

(華北電力大學 新能源電力系統(tǒng)國家重點實驗室，北京 102206)

?

計算多介質電場節(jié)點場強的有限元-邊界元法

謝裕清，李琳

(華北電力大學 新能源電力系統(tǒng)國家重點實驗室，北京 102206)

為了解決節(jié)點有限元法在多介質電場計算中場域邊界及介質分界面處場強計算精度不高的問題，提出了一種有限元-邊界元混合算法。該算法首先應用節(jié)點有限元法計算多介質場域中的節(jié)點電位，然后將場域劃分成多個單一介質的子域，子域邊界處的電位通過有限元法計算的節(jié)點電位映射插值得到。接下來對各個子域分別采用邊界元法計算邊界處的法向場強，子域內部任意位置的場強則通過邊界積分方程計算得到。數(shù)值計算結果表明該方法在場域邊界以及分界面處節(jié)點場強的計算精度比傳統(tǒng)的面積加權平均算法高，場域內部節(jié)點場強的計算誤差也較低。所提方法不僅適合于靜電場的計算問題，也適合于其它場強在介質分界面處不連續(xù)的矢量場計算問題。

有限元-邊界元法；多介質；介質分界面；邊界場強

0　引　言

在電磁場數(shù)值計算領域中，通過節(jié)點有限元方法求解電位函數(shù)的理論及技術已經非常成熟，電位的計算精度也較高。然而，怎樣通過節(jié)點位函數(shù)獲得較高精度的節(jié)點場強一直是學者較為關注的問題。文獻[1]從互補變分原理出發(fā)得到了求解全域節(jié)點場強的有限元離散公式。文獻[2]基于面矢基函數(shù)，應用加權余量法對電位的梯度函數(shù)進行有限元離散，通過散度定理對電位梯度進行轉換，避免了直接對電位進行數(shù)值微分，提高了計算精度。文獻[3-5]基于超收斂理論求得單元內精度最高點的場強，然后利用該點的場強通過插值理論得到網(wǎng)格的節(jié)點場強。上述方法對于單一介質模型能夠獲得較高的計算精度，但是對于多介質模型在分界面處的節(jié)點場強的計算效果較差。

在實際工程問題中，靜電場場域的邊界以及介質交界面處的法向場強是研究介質絕緣問題的重要參數(shù)。為了解決節(jié)點有限元法不能很好地計算介質分界面處的場強的問題，棱邊有限元法得到迅速發(fā)展。該方法直接選用場矢量作為求解變量，采用棱單元插值保證場強切向方向連續(xù)，而不強制法向方向連續(xù)，可以避免通過位函數(shù)求解場量而引起的微分誤差等問題[6-8]。該方法計算精度較高，但是與節(jié)點有限元相比，棱邊有限元法計算代價較高。

對于場域邊界及內部的場強，邊界元法可以得到比較高精度的計算結果。該方法將靜電場的邊值問題轉換成邊界積分方程問題，利用有限元離散技術對場域邊界進行離散,結合邊界積分方程將靜電場邊值問題轉換成代數(shù)方程問題求解。該方法僅僅需要對邊界單元進行網(wǎng)格剖分，其求解問題的空間維數(shù)比有限元低，計算精度高，易處理開域問題等特點[9-13]。然而，邊界元法對于多介質模型處理起來也較為困難。

結合有限元以及邊界元法的計算優(yōu)勢，本文提出求解多介質靜電場節(jié)點場強的有限元-邊界元混合算法。首先以節(jié)點有限元法求解場域節(jié)點電位，再以邊界元法計算從場域提取的單一介質子域的節(jié)點場強。由于子域是由單一介質組成，場域分界面同時也是兩個不同子域的交界面，因此可以避免由于介質分界面法向場強不連續(xù)而需要在分界面處額外布置節(jié)點的困難，從而可以提高介質分界面處的電場強度的計算精度。此外，場域邊界處的節(jié)點場強也比傳統(tǒng)的節(jié)點有限元法的計算精度高。

1　有限元-邊界元法計算節(jié)點場強流程

圖1為多介質靜電場模型示意圖，其包含n種介質，場域邊界Γ1滿足第一類邊界條件，邊界Γ2滿足第二類邊界條件。在靜電場數(shù)值計算中，電位在場域內處處連續(xù)，場域中的電位可以通過節(jié)點有限元法得到精度較高的數(shù)值解。在介質交界面處，電場強度僅僅切向連續(xù)，切向電場強度可以通過電位微分方程計算得到，而法向場強在交界面處不連續(xù)，采用電位微分法的方式求得的場強誤差較大。

圖1　多介質靜電場模型Fig.1　Multi-medium static electrical field model

本文結合有限元法與邊界元法的特點，提出有限元邊界元混合算法計算場域中的節(jié)點場強，其計算流程圖如圖2所示。首先對于整體靜電場模型采用有限元網(wǎng)格離散，建立基于電位的節(jié)點有限元方程，得到整體場域的節(jié)點電位。然后，將場域劃分成多個由單一介質構成的子域，提取子域的幾何信息。使用邊界單元離散子域邊界，子域的邊界節(jié)點電位則通過有限元法求得的場域節(jié)點電位映射插值得到。最后，采用邊界元法建立關于子域法向場強的邊界元方程，求得子域的邊界法向場強。由于在介質分界面上電場的切向分量連續(xù)，仍然可以通過節(jié)點電位的微分進行計算并得到較高精度的解。子域內部任意位置的場強則可以結合子域邊界上的節(jié)點電位以及節(jié)點法向場強通過邊界積分方程計算得到。

圖2　有限元-邊界元計算節(jié)點場強流程圖Fig.2　Flowchart of finite element and boundary element mixed method for solving node electrical field intensity

實際計算過程中，子域的選取可以根據(jù)實際計算需要選擇，僅需要滿足單一介質條件即可。如果僅關注整體場域中某一小部分區(qū)域的電場強度情況，可以圍繞這部分區(qū)域構造子域，此時該子域的計算規(guī)模比整體場域的計算規(guī)模要小得多，節(jié)約了計算資源。

2　有限元法計算全域節(jié)點電位

圖1所示的多介質電場計算模型，標量電位滿足的基本方程及邊界條件為

(1)

運用變分原理，將靜電場邊值問題的微分方程式(1)轉化為泛函求極值問題。式(1)所對應的等效極值泛函為[14]

(2)

上式泛函取得極值的解也即偏微分方程(1)式的解，第一類邊界條件為變分約束條件。應用有限元單元剖分場域，選取合適的插值函數(shù)，構造電位φ的近似函數(shù)為

(3)

式中：n為場域離散網(wǎng)格節(jié)點總數(shù)；Nj為節(jié)點電位φj對應的插值基函數(shù)。將式(3)帶入式(2)，并根據(jù)泛函極值條件?J(φ)/?φi=0(i=1,2,…,n)，可以得到以節(jié)點電位φ1,φ2, …,φn為變量的有限元代數(shù)方程組：

(4)

第一類邊界條件采用強加邊界條件法帶入上式的有限元方程組式[15]，應用高斯消去法或者迭代法求解該方程組可得到場域各個離散網(wǎng)格點處的節(jié)點電位。

3　邊界元法計算子域邊界場強及內部節(jié)點場強

如圖3所示的單介質子域是全場域的部分求解區(qū)域，該子域的介電常數(shù)εi為一定值，子域內部的電荷密度為ρi，子域邊界Γi上的節(jié)點電位φi可從全域節(jié)點有限元法所計算的子域邊界上的節(jié)點電位映射得到。

圖3　單介質子域Fig.3　Single-medium subdomain

該子域的電位控制方程及邊界條件如下所示：

(5)

根據(jù)格林定理：

(6)

式中:G為自有空間中的格林函數(shù)，滿足2G=-δ，式中δ=δ(P,Q)為狄拉克函數(shù)；P，Q分別為場域中的場點及源點。在二維自由空間中的格林函數(shù)[16]為

(7)

式中：r為從源點到場點之間的距離。在三維自由空間中，格林函數(shù)為[16]

(8)

令ρi/εi=f，將子域電位控制方程式(3)帶入格林公式式(4)中，結合格林函數(shù)的基本性質可以得到子域電位的邊界積分方程為[17]

(9)

式中：在光滑邊界上ci=1/2，場域內部ci=1，場域外部ci=0。φi為空間任意位置的電位，而式(9)右邊的電位φ則表示子域邊界上的電位。

對于子域邊界上，式(9)中的φi以及φ均已給定，而?φ/?n為待求量。根據(jù)場強與電位的梯度關系，子域邊界的法向場強為En=-?φ/?n，將其帶入邊界積分方程式(9)，此時將形成待求量為子域邊界法向場強的邊界積分方程。應用邊界單元將子域邊界離散成一系列的邊界元單元，使用單元插值函數(shù)構造邊界上電位φ以及邊界法向場強En的近似值為

(10)

式中：n為子域邊界節(jié)點總數(shù)；φj，En,j分別為節(jié)點電位以及節(jié)點表面法向場強。

將式(10)帶入式(9)結合子域的邊界電位通過配點的方式可以得到求解子域邊界場強的離散代數(shù)方程組。

(11)

式中：e為離散邊界單元編號；ne為離散邊界單元的總數(shù)量；Gi為場點固定于第i個邊界節(jié)點時的格林函數(shù)；δij為克羅內克函數(shù)(當場點i與源點j重合時其值為1，其余為0)，區(qū)域D為子域中含有電荷密度的區(qū)域。

將式(11)寫成矩陣形式為

(12)

該式即為求解電位邊界法向場強的矩陣方程。式中En為邊界節(jié)點法向場強組成的向量，φ為邊界節(jié)點電位組成的向量，A和B為與式(11)對應的系數(shù)矩陣。通過求解矩陣方程式(12)即可得到子域邊界處的節(jié)點法向場強。對于子域邊界處的切向節(jié)點場強可以通過邊界電位的微分得到。子域內部任意位置的電位與節(jié)點場強可通過邊界積分方程計算得到。

子域內部任意位置的節(jié)點電位可以通過全域有限元計算的節(jié)點電位映射插值得到，也可以通過式(9)的邊界積分方程計算得到。子域內部求解電位的邊界積分方程為

(13)

式中：φi為子域內部待求的節(jié)點電位。

根據(jù)電場強度與電位的梯度關系E=-φ，結合式(13)，可以得到求解子域內部節(jié)點場強分量的邊界積分方程為

式中：ξ為空間坐標x，y，z的一個方向，即Eξ,i為內部節(jié)點i在ξ方向的電場強度的分量。

比較式(14)及式(15)可以看到采用邊界元法計算的場域電位與場強具有相同的計算精度。使用邊界單元離散子域邊界，將邊界電位及法向場強的插值公式(10)式帶入(14)式可以得到計算子域內點的電場強度分量Eξ,i的計算方程式為

(15)

從上式可以看出計算子域內部任意位置的電場強度僅需要對子域的邊界以及帶有電荷源的區(qū)域進行網(wǎng)格剖分，并不需要對全部場域進行離散。計算子域某個內點的電場強度僅需要將該內點的空間坐標帶入式(15)即可。與節(jié)點有限元法計算節(jié)點場強的方法相比，該方法靈活方便，而且計算效率更快。

4　算例分析

4.1計算模型及面積加權平均算法簡介

本文采用一個場強具有解析解的雙層介質同軸電纜模型，比較本文所提出的有限元―邊界元混合算法與傳統(tǒng)的面積加權平均算法求解節(jié)點場強的計算精度。同軸電纜的模型示意圖如圖4所示。

圖4　雙層介質同軸電纜模型Fig.4　Two-layer medium coaxial-cable model

圖中所示的同軸電纜模型內半徑R1=1 mm，介質分界面處的半徑R2=2 mm，外半徑R3=4 mm。介電常數(shù)ε1=ε0，ε2=5ε0。所加外電壓為U0=100 V。該計算模型軸向場強的解析解為

(16)

式中：r為場域任意點與坐標原點的距離。

傳統(tǒng)的求取場域節(jié)點場強的方法是對單元位函數(shù)直接進行數(shù)值微分，該方法不僅精度低而且由于一個節(jié)點與多個單元相連接使得連續(xù)介質中計算的節(jié)點場強也是多值的，這與實際物理規(guī)律相悖。面積加權平均算法計算節(jié)點場強的基本思路是對位函數(shù)進行數(shù)值微分取得單元重心處的場強，而節(jié)點場強則是與該節(jié)點相關聯(lián)的單元重心處的場強的面積加權平均，其具體計算式為[18]

(17)

式中：Ek,i表示第i個節(jié)點k(k可為空間坐標方向x，y，z任一方向)方向的場強分量；Ne表示與節(jié)點i相關聯(lián)的單元個數(shù)；Ek,eij表示與節(jié)點i相關聯(lián)的第j個單元eij重心處k方向的場強分量；Seij為單元eij的面積。

4.2場域邊界及介質分界面節(jié)點場強精度比較

首先，對圖4所示的同軸電纜模型進行有限元網(wǎng)格劃分求取節(jié)點電位以及運用面積加權平均法求取節(jié)點場強。網(wǎng)格劃分單元采用四邊形一次單元，場域網(wǎng)格數(shù)為7 200，節(jié)點數(shù)為7 320。接著，將場域按照介質類型劃分成兩個子域，然后對子域進行邊界元計算子域的邊界法向場強及子域內部節(jié)點場強。使用一次線性單元對子域邊界進行邊界網(wǎng)格劃分，兩個子域的網(wǎng)格數(shù)及節(jié)點數(shù)均為120。

在本模型的計算中，場強的大小僅與節(jié)點處的半徑有關。應用有限元-邊界元法(FBM)與面積加權平均法(AWA)計算模型邊界及分界面處場強的大小以及誤差分析結果如表1所示。表中位置r=1.00 mm處表示場域內表面邊界處，r=2.00(-)表示分界面處偏向內表面一側的位置，r=2.00(+)表示分界面處偏向外表面一側的位置，r=4.00 mm處表示場域外表面邊界處。從表中的計算結果可以看出，本文所提出的有限元―邊界元混合算法計算場域邊界及分界面處的場強大小誤差均小于1%，計算精度較高。而面積加權平均算法在場域邊界處的計算精度要比本文所提出的算法計算精度低，而在介質分界面處的計算結果與解析解相差很大，其主要原因是由于介質分界面處的法向場強不連續(xù)所導致的。

表1　場域邊界及介質分界面處場強大小及誤差分析

4.3場域內部節(jié)點場強精度分析

應用有限元-邊界元法與面積加權平均法計算場域內部區(qū)域的節(jié)點場強，選取場域一條沿半徑方向上的節(jié)點場強進行分析，選取的節(jié)點位置為距原點半徑分別為r=1.20，1.50，1.80，2.50，3.00，3.50 mm的位置。選取節(jié)點所計算的場強大小及誤差分析如表2所示。

表2　沿半徑方向場強大小及誤差分析

表2中顯示的計算結果表明對于子域內部節(jié)點場強的計算，傳統(tǒng)的面積加權平均法與有限元-邊界元計算誤差均比較小，但是傳統(tǒng)的面積加權平均法計算誤差更小。實質上，這兩種算法子域內部的計算誤差可比性不大，原因在于面積加權平均法計算節(jié)點場強的精度與子域內部網(wǎng)格的剖分的密度有很大的關系。而有限元-邊界元法計算子域內部的節(jié)點場強不需要對子域進行網(wǎng)格剖分，其計算精度與子域表面的節(jié)點法向場強和電位的計算精度以及數(shù)值積分的計算誤差有關。

有限元-邊界元法與面積加權平均法計算的節(jié)點場強沿半徑方向的計算結果曲線圖如圖5所示。從這兩種算法求解的電纜沿著半徑方向場強大小的計算結果可以看出面積加權算法計算的場強大小在分界面處與解析解相差很大，計算值是介質分界面兩邊單元場強的面積加權平均，無法反映場強在介質分界面處不連續(xù)的特點。與此相比，有限元-邊界元法求解場強的結算結果在場域邊界、內部以及分界面處均與解析解一致，具有較高的計算精度。

圖5　沿半徑方向電場強度分布圖Fig.5　Electrical field intensity along cable radius

5　結　論

本文提出了一種有限元邊界元混合算法求解多介質靜電場電場強度的計算方法。對一具有解析解的雙層介質同軸電纜模型的電場進行了計算分析，可以得出以下結論：

(1)有限元邊界元混合算法計算多介質節(jié)點場強可以很好地解決節(jié)點場強在分界面處法向場強不連續(xù)所引起的計算困難。算法在介質分界面及場域邊界處的節(jié)點場強計算精度要比傳統(tǒng)的面積加權平均算法高。

(2)應用邊界元法計算子域內部場強，不需要再對子域內部進行網(wǎng)格劃分。僅需要輸入所計算點的空間坐標然后運用邊界積分方程即可計算得到，比有限元法需要通過網(wǎng)格節(jié)點場強進行插值計算方便。

(3)本方法子域的劃分僅需要滿足單一介質的閉合區(qū)域條件即可。實際計算過程中可以僅對部分需要計算場強的區(qū)域提取子域信息，并不一定需要以介質分界面作為子域的邊界，也可在單一介質內部劃分子域，從而提高計算效率。

(4)本文僅對二維模型進行了算法驗證，對于三維模型也可以利用本文所提出的算法進行計算。此外，對于其它類似的問題，例如多介質磁場問題也可以通過本文提出的算法進行計算。

[1] 崔翔, 謝羲. 計算電磁場量 E 和 B 的互補變分方法 [J]. 中國電機工程學報, 1988, 8(2): 22-32.

[2] 左鵬, 鄒軍, 袁建生. 三維有限元中提高由已知節(jié)點電位求場強計算精度的全域節(jié)點場強法 [J]. 中國電機工程學報, 2015, 35(5): 1243-1249.

[3] CAO W, ZHANG Z, ZOU Q. Superconvergence of discontinuous Galerkin methods for linear hyperbolic equations [J]. SIAM Journal on Numerical Analysis, 2014, 52(5): 2555-2573.

[4] SHI D Y, LI M H. Superconvergence analysis of the stable conforming rectangular mixed finite elements for the linear elasticity problem[J]. Journal of Computational Mathematics, 2014, 32(2): 205-214.

[5] ZHANG Z, NAGA A. A new finite element gradient recovery method: Superconvergence property [J]. SIAM Journal on Scientific Computing, 2005, 26(4): 1192-213.

[6] 張秀敏, 苑津莎, 徐永生, 等. 基于工程損耗模型的棱邊有限元與節(jié)點有限元的算法比較 [J]. 電工技術學報, 2003, 18(3): 41-47.

[7] MUR G. Edge elements, their advantages and their disadvantages [J]. IEEE Transactions on Magnetics, 1994, 30(5): 3552-3557.

[8] 苑津莎, 張秀敏, 王志明. 棱邊有限元法在工程渦流問題中的應用研究[J]. 華北電力大學學報, 2004, 31(3): 1-7.

[9] SIMPSON R N, BORDAS S P A, TREVELYAN J, et al. A two-dimensional isogeometric boundary element method for elastostatic analysis[J]. Computer Methods in Applied Mechanics and Engineering, 2012, 209: 87-100.

[10] LIU Y J, MUKHERJEE S, NISHIMURA N, et al. Recent advances and emerging applications of the boundary element method [J]. Applied Mechanics Reviews, 2011, 64(3): 1-38.

[11] 王澤忠, 趙良, 劉之方. 二維開域靜電場有限元與邊界元迭代解法的研究[J]. 華北電力大學學報, 2002, 29(z1): 36-40.

[12] 潘曉彤, 王澤忠, 方舟, 等. 直流輸電換流閥系統(tǒng)表面電場分析與均壓環(huán)設計[J]. 現(xiàn)代電力, 2014, 31(1): 68-73.

[13] REN Z, KALSCHEUER T, GREENHALGH S, et al. A hybrid boundary element-finite element approach to modeling plane wave 3D electromagnetic induction responses in the Earth[J]. Journal of Computational Physics, 2014, 258: 705-717.

[14] 倪成, 馮慈璋, 倪光正. 電磁能量和參數(shù)計算的互補和互補—對偶能量法[J]. 西安交通大學學報, 1986, 20(3): 91-101.

[15] 王澤忠. 簡明電磁場數(shù)值計算 [M]. 北京：機械工業(yè)出版社, 2011:107-116.

[16] AIELLO G, ALFONZETTI S, BORZG, et al. Comparing FEM-BEM and FEM-DBCI for open-boundary electrostatic field problems [J]. The European Physical Journal Applied Physics, 2007, 39(2): 143-148.

[17] 吳洪潭. 邊界元法在傳熱學中的應用 [M]. 北京：國防工業(yè)出版社, 2008:14-17.

[18] 倪光正, 楊仕友, 邱捷. 工程電磁場數(shù)值計算 [M]. 北京：機械工業(yè)出版社, 2010:140-148.

Finite Element Method-boundary Element Method for Calculation of Nodal Electric Field Intensity in Multi-medium Electric Field

XIE Yuqing, LI Lin

(State Key Laboratory of Alternate Electrical Power System with Renewable Energy Sources, North China Electric Power University, Beijing 102206, China)

A finite element-boundary element hybrid method is proposed in this paper to solve the problem of inaccuracy when calculating the electric field intensity of field boundary and dielectric interface in the multi-medium electric field with nodal finite element method. Finite element method is firstly used to calculate the node electric potential of the field domain. Then, the whole domain is divided into several single-medium subdomains in which the boundary electric potential can be calculated by the mapping interpolation of node electric potential from the result of finite element method. Next, the boundary element method is used to calculate the normal field intensity at the boundaries in subdomains, and the interior electric field intensity will be calculated by boundary integral equation. The results show that such calculation method yields more accurate results when calculating the nodal electric field intensity of the field boundary and the interface than the traditional area-weighted average method, and the calculation error of the node intensity in the interior field is also of low level. The presented method could not only be used in the static electric field, but also can be used to solve calculation problems in other vector field with discontinuous field intensity in the dielectric interface.Key words：finite element method-boundary element method; multi-medium; dielectric interface; boundary electric field intensity

10.3969/j.ISSN.1007-2691.2016.04.04

2015-10-11.

國家自然科學基金資助項目(51277064)；中央高校基本科研業(yè)務費專項資金資助項目(JB2015131).

謝裕清(1991-)，男，博士研究生，主要研究方為電磁場數(shù)值計算及多物理場耦合計算；李琳(1962-)，男，教授，博士生導師，主要研究方向為電磁場理論及其應用和電力系統(tǒng)電磁兼容。

TM151

A

1007-2691(2016)04-0021-06