謝裕清，李　琳

(華北電力大學(xué) 新能源電力系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室，北京 102206)

?

計(jì)算多介質(zhì)電場(chǎng)節(jié)點(diǎn)場(chǎng)強(qiáng)的有限元-邊界元法

謝裕清，李琳

(華北電力大學(xué) 新能源電力系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室，北京 102206)

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

有限元-邊界元法；多介質(zhì)；介質(zhì)分界面；邊界場(chǎng)強(qiáng)

0　引　言

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

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

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

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

1　有限元-邊界元法計(jì)算節(jié)點(diǎn)場(chǎng)強(qiáng)流程

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

圖1　多介質(zhì)靜電場(chǎng)模型Fig.1　Multi-medium static electrical field model

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

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

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

2　有限元法計(jì)算全域節(jié)點(diǎn)電位

圖1所示的多介質(zhì)電場(chǎng)計(jì)算模型，標(biāo)量電位滿足的基本方程及邊界條件為

(1)

運(yùn)用變分原理，將靜電場(chǎng)邊值問題的微分方程式(1)轉(zhuǎn)化為泛函求極值問題。式(1)所對(duì)應(yīng)的等效極值泛函為[14]

(2)

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

(3)

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

(4)

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

3　邊界元法計(jì)算子域邊界場(chǎng)強(qiáng)及內(nèi)部節(jié)點(diǎn)場(chǎng)強(qiáng)

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

圖3　單介質(zhì)子域Fig.3　Single-medium subdomain

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

(5)

根據(jù)格林定理：

(6)

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

(7)

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

(8)

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

(9)

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

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

(10)

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

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

(11)

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

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

(12)

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

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

(13)

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

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

式中：ξ為空間坐標(biāo)x，y，z的一個(gè)方向，即Eξ,i為內(nèi)部節(jié)點(diǎn)i在ξ方向的電場(chǎng)強(qiáng)度的分量。

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

(15)

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

4　算例分析

4.1計(jì)算模型及面積加權(quán)平均算法簡(jiǎn)介

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

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

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

(16)

式中：r為場(chǎng)域任意點(diǎn)與坐標(biāo)原點(diǎn)的距離。

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

(17)

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

4.2場(chǎng)域邊界及介質(zhì)分界面節(jié)點(diǎn)場(chǎng)強(qiáng)精度比較

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

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

表1　場(chǎng)域邊界及介質(zhì)分界面處場(chǎng)強(qiáng)大小及誤差分析

4.3場(chǎng)域內(nèi)部節(jié)點(diǎn)場(chǎng)強(qiáng)精度分析

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

表2　沿半徑方向場(chǎng)強(qiáng)大小及誤差分析

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

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

圖5　沿半徑方向電場(chǎng)強(qiáng)度分布圖Fig.5　Electrical field intensity along cable radius

5　結(jié)　論

本文提出了一種有限元邊界元混合算法求解多介質(zhì)靜電場(chǎng)電場(chǎng)強(qiáng)度的計(jì)算方法。對(duì)一具有解析解的雙層介質(zhì)同軸電纜模型的電場(chǎng)進(jìn)行了計(jì)算分析，可以得出以下結(jié)論：

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

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

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

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

[1] 崔翔, 謝羲. 計(jì)算電磁場(chǎng)量 E 和 B 的互補(bǔ)變分方法 [J]. 中國電機(jī)工程學(xué)報(bào), 1988, 8(2): 22-32.

[2] 左鵬, 鄒軍, 袁建生. 三維有限元中提高由已知節(jié)點(diǎn)電位求場(chǎng)強(qiáng)計(jì)算精度的全域節(jié)點(diǎn)場(chǎng)強(qiáng)法 [J]. 中國電機(jī)工程學(xué)報(bào), 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é)點(diǎn)有限元的算法比較 [J]. 電工技術(shù)學(xué)報(bào), 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] 苑津莎, 張秀敏, 王志明. 棱邊有限元法在工程渦流問題中的應(yīng)用研究[J]. 華北電力大學(xué)學(xué)報(bào), 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] 王澤忠, 趙良, 劉之方. 二維開域靜電場(chǎng)有限元與邊界元迭代解法的研究[J]. 華北電力大學(xué)學(xué)報(bào), 2002, 29(z1): 36-40.

[12] 潘曉彤, 王澤忠, 方舟, 等. 直流輸電換流閥系統(tǒng)表面電場(chǎng)分析與均壓環(huán)設(shè)計(jì)[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ì)算的互補(bǔ)和互補(bǔ)—對(duì)偶能量法[J]. 西安交通大學(xué)學(xué)報(bào), 1986, 20(3): 91-101.

[15] 王澤忠. 簡(jiǎn)明電磁場(chǎng)數(shù)值計(jì)算 [M]. 北京：機(jī)械工業(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] 吳洪潭. 邊界元法在傳熱學(xué)中的應(yīng)用 [M]. 北京：國防工業(yè)出版社, 2008:14-17.

[18] 倪光正, 楊仕友, 邱捷. 工程電磁場(chǎng)數(shù)值計(jì)算 [M]. 北京：機(jī)械工業(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.

國家自然科學(xué)基金資助項(xiàng)目(51277064)；中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金資助項(xiàng)目(JB2015131).

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

TM151

A

1007-2691(2016)04-0021-06