靜電場(chǎng)分析的比例邊界有限元法

劉 俊， 林 皋＊， 王 復(fù) 明， 李 建 波

（1.大連理工大學(xué) 建設(shè)工程學(xué)部 水利工程學(xué)院，遼寧 大連 116024；2.大連理工大學(xué) 海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024；3.鄭州大學(xué) 水利與環(huán)境學(xué)院，河南 鄭州 450002）

0 引 言

靜電場(chǎng)求解主要有解析法、半解析法和數(shù)值解法三大類.解析法的優(yōu)點(diǎn)是能夠得到精確的解析式，直觀地看出各物理量間的變化關(guān)系，但只適用于簡(jiǎn)單的幾何邊界、簡(jiǎn)單的物理參數(shù).20世紀(jì)80年代初期，人們對(duì)半解析解進(jìn)行了研究，其主要方法有級(jí)數(shù)法、多極點(diǎn)法、廣義多極技術(shù)、多極理論方法和新型等效源法[1～4]等.這些方法在電磁計(jì)算方面取得了較好的效果和較廣泛的應(yīng)用，但該類方法有時(shí)也只適用一些幾何形狀比較規(guī)則的靜電場(chǎng)問(wèn)題.數(shù)值計(jì)算方法是隨著計(jì)算機(jī)的問(wèn)世應(yīng)運(yùn)而生的.該類方法中有有限差分法、有限元法、邊界元法、矩量法、無(wú)網(wǎng)格法[5～9]等.其中有限差分法網(wǎng)格剖分容易，編制程序方便，但對(duì)不規(guī)則的復(fù)雜邊界，網(wǎng)格剖分缺少靈活性.有限元法是分析靜電場(chǎng)問(wèn)題適應(yīng)性最強(qiáng)、最通用的方法之一，因?yàn)樵摲椒ㄟm合于含有復(fù)雜媒質(zhì)（包括非線性媒質(zhì)以及各向異性媒質(zhì)等）問(wèn)題的數(shù)值分析.但該方法在處理含有場(chǎng)值奇異點(diǎn)以及不同材料交接點(diǎn)等問(wèn)題中遇到很大的挑戰(zhàn)，通常的處理方法是在這些點(diǎn)周圍增加節(jié)點(diǎn)，但這勢(shì)必需要更多的前處理和計(jì)算工作時(shí)間，當(dāng)然也有很多改進(jìn)的有限元法，例如提高形函數(shù)階段或者使用超級(jí)單元等.與此同時(shí)，邊界元法也可以很好地解決該類問(wèn)題.它降低了待求問(wèn)題的維數(shù)，簡(jiǎn)化了數(shù)據(jù)的處理，直接求解無(wú)界邊值問(wèn)題精度較高.但該方法的缺點(diǎn)是有時(shí)找不到對(duì)應(yīng)的基本解.無(wú)網(wǎng)格法只需節(jié)點(diǎn)，不需單元，適合處理復(fù)雜邊界問(wèn)題，計(jì)算精度高、收斂速度快.但針對(duì)某些問(wèn)題，其基函數(shù)的選取及節(jié)點(diǎn)布置對(duì)計(jì)算精度的影響等問(wèn)題還有待于進(jìn)一步研究.從以上分析可以看出，各種解法有各自的優(yōu)缺點(diǎn)，采用一種方法，有時(shí)得不到滿意的結(jié)果，為此，近年來(lái)出現(xiàn)了將不同方法相結(jié)合的方法.例如，有限元－邊界元法、有限元與無(wú)單元耦合法、多極理論－邊界元耦合法、邊界元法與無(wú)網(wǎng)格局部Petrov－Galerkin法的耦合法[4、10～12]等.

比例邊界有限元法（scaled boundary finite element method，SBFEM）是20世紀(jì)90年代由Wolf和Song等[13]率先提出和發(fā)展起來(lái)的一種半解析的數(shù)值方法.該方法綜合了有限元法和邊界元法的優(yōu)點(diǎn)，只需用有限元離散部分邊界，從而實(shí)現(xiàn)了將問(wèn)題降低一維的目的，而在沒(méi)有離散的坐標(biāo)方向利用解析方法求解.與傳統(tǒng)有限元法比較，它顯著降低了計(jì)算工作量.相對(duì)于邊界元法，它不需要基本解，也不存在積分的奇異性問(wèn)題，具有較高的計(jì)算精度，對(duì)于各向異性，以及物理性質(zhì)沿某一方向發(fā)生特定變化的問(wèn)題處理也比較方便.目前這種方法已被用于求解有限域、無(wú)限域的彈性靜力問(wèn)題、動(dòng)力學(xué)問(wèn)題、水庫(kù)壩體與水體的相互作用問(wèn)題、斷裂力學(xué)、繞流場(chǎng)和勢(shì)流場(chǎng)問(wèn)題、液體晃蕩問(wèn)題、波浪與結(jié)構(gòu)相互作用問(wèn)題、聲學(xué)問(wèn)題、滲流問(wèn)題[14]等.本文首次將該方法應(yīng)用于靜電場(chǎng)問(wèn)題的求解.從靜電場(chǎng)控制方程——拉普拉斯方程出發(fā)，結(jié)合加權(quán)余量法推導(dǎo)靜電場(chǎng)分析的比例邊界有限元方程，引入特征值問(wèn)題對(duì)該方程進(jìn)行求解并得出電位和電場(chǎng)的計(jì)算公式；以二維靜電場(chǎng)為例通過(guò)數(shù)值算例驗(yàn)證比例邊界有限元法求解靜電場(chǎng)問(wèn)題的高效性和精確性.

1 靜電場(chǎng)的比例邊界有限元方程推導(dǎo)

描述靜電場(chǎng)（域內(nèi)無(wú)電荷）特性的控制方程為拉普拉斯方程（ 為電位）：

存在兩類邊界條件：一類是在邊界上給定電位值，一類是給定電場(chǎng)值，即

對(duì)于以上靜電場(chǎng)控制方程和邊界條件問(wèn)題的比例邊界有限元方程推導(dǎo)，首先需建立笛卡兒坐標(biāo)系統(tǒng)和比例邊界有限元坐標(biāo)系統(tǒng)轉(zhuǎn)換關(guān)系.考察圖1所示的有限區(qū)域和無(wú)限區(qū)域，在域內(nèi)選擇一點(diǎn)O作為比例邊界有限元的相似中心（要求整個(gè)計(jì)算域必須在相似中心可視化的范圍之內(nèi)，也就是邊界各點(diǎn)與相似中心均可用直線進(jìn)行連接），建立以點(diǎn)O為中心的ξ和s坐標(biāo)體系.其中有限區(qū)域0＝ξ0≤ξ≤ξ1＝1、無(wú)限區(qū)域1＝ξ0≤ξ≤ξ1＝∞.計(jì)算區(qū)域內(nèi)的點(diǎn)在比例坐標(biāo)系統(tǒng)下可以表示為

圖1 比例邊界有限元法計(jì)算示意圖Fig.1 Sketch for scaled boundary finite elementmethod

兩坐標(biāo)系統(tǒng)可通過(guò)雅克比矩陣聯(lián)系起來(lái)：

其中

在比例坐標(biāo)下梯度的算子可以表示為

其中

對(duì)方程（1）、（2）、（3）利用加權(quán)余量法和格林第一公式可得（w為權(quán)函數(shù)）

根據(jù)等參變換概念，電位函數(shù)可采用相同的插值函數(shù)N（s）進(jìn)行離散：

方程（11）中權(quán)函數(shù)w也用同樣的形函數(shù)表達(dá)：

將式（12）、（13）代入式（11）可得

其中

比例坐標(biāo)下的空間微元面積

比例坐標(biāo)下的邊界微元線長(zhǎng)度

其中τξ為關(guān)于ξ的函數(shù).

對(duì)式（14）中含w（ξ，s），ξ項(xiàng)做分部積分，并將式（17）、（18）代入可得

方程（19）的系數(shù)為

考慮方程（19）中w（ξ）任意性可得

方程（23）為比例邊界有限元的基本方程.方程（24）、（25）為需要滿足的計(jì)算域內(nèi)、外邊界條件.對(duì)于有限區(qū)域（0＝ξ0≤ξ≤ξ1＝1），內(nèi)邊界條件式（24）為相似中心一點(diǎn)，外邊界條件式（25）為離散邊界；對(duì)于無(wú)限區(qū)域（1＝ξ0≤ξ≤ξ1＝∞），內(nèi)邊界條件式（24）為離散邊界，外邊界條件式（25）為無(wú)窮遠(yuǎn)處外邊界條件.

2 比例邊界有限元方程求解

式（23）為二階Euler－Cauchy齊次方程，為了便于降階求解，從式（24）、（25）特性看出，引入Q（ξ）作為 （ξ）的對(duì)偶變量：

由此可得具有兩倍未知數(shù)的變量：

可將式（23）化為狀態(tài)方程：

其中

由于Z陣為Hamilton陣，可以通過(guò)求解Z的特征值問(wèn)題來(lái)得到式（30）的解（式（31）特征值中成對(duì)出現(xiàn)，互為負(fù)號(hào)）：

由此方程（30）的解為

將方程（27）代入方程（31）可得

式中：c1、c2為積分常數(shù).對(duì)于有限域問(wèn)題，ξ＝0處的 取有限值，故c2＝0；對(duì)于無(wú)限域問(wèn)題，ξ＝∞處的 取有限值，故c1＝0.有限域的積分常數(shù)c1和無(wú)限域的積分常數(shù)c2均可由邊界條件確定.積分常數(shù)確定后，可由式（32）、（12）確定域內(nèi)任意點(diǎn)的電位和由E＝－確定域內(nèi)任意點(diǎn)電場(chǎng)強(qiáng)度.

3 算例計(jì)算分析

為了驗(yàn)證該方法數(shù)值模擬精度和高效性，本文首先選擇了一個(gè)具有解析解和其他數(shù)值解的經(jīng)典算例以便于對(duì)比.

算例1 一長(zhǎng)直金屬槽，側(cè)壁與底面的電位為0，而頂面蓋電位 （x，b）＝U0sin（πx/a），需求出域內(nèi)的電位和電場(chǎng).比例邊界有限元計(jì)算示意圖如圖2所示，在仿真中，取a＝3.0 m，b＝1.0 m，U0＝1.0 V，相似中心放在域中心.

圖2 比例邊界有限元法算例1計(jì)算示意圖Fig.2 Sketch for the first example using scaled boundary finite element method

表1為本文方法與文獻(xiàn)[9]列出的有限差分法（FDM）、徑向基無(wú)網(wǎng)格法（RBF）計(jì)算電位最大誤差、相對(duì)均方根誤差對(duì)比（電位單位為V）.相對(duì)均方根誤差計(jì)算公式 （i＿exact為解析解，i＿calc為數(shù)值解）為

表1數(shù)據(jù)表明，比例邊界有限元法精度非常高.在電場(chǎng)強(qiáng)度（單位為V/m）分析中，由于文獻(xiàn)[9]沒(méi)有給出相似節(jié)點(diǎn)分布下的計(jì)算結(jié)果，本文單獨(dú)給出了本方法計(jì)算電場(chǎng)強(qiáng)度的相對(duì)均方根誤差.邊界為32節(jié)點(diǎn)的Ex相對(duì)均方根誤差為1.8720%；邊界為64節(jié)點(diǎn)的Ex相對(duì)均方根誤差為0.5217%.邊界為32節(jié)點(diǎn)的Ey相對(duì)均方根誤差為2.0310%；邊界為64節(jié)點(diǎn)的Ey相對(duì)均方根誤差為0.7617%.同樣可以看出本方法精度也比較高.電位解析解、SBFEM數(shù)值解（邊界64節(jié)點(diǎn)）電位等值線分布如圖3所示.從圖中可以看出：比例邊界有限元法和解析解非常吻合.與此同時(shí)，由于文獻(xiàn)[9]沒(méi)給出其他兩種計(jì)算方法的計(jì)算耗時(shí)，為此，本文采用Intel Core Q8200（2.33 GHz）計(jì)算機(jī)對(duì)本算例的SBFEM 64節(jié)點(diǎn)模型進(jìn)行計(jì)算時(shí)間測(cè)試，其消耗的時(shí)間不到1 s，由此看出計(jì)算效率比較高，而且比例邊界有限元法前期準(zhǔn)備（即單元?jiǎng)澐郑┮采?，主要的?jì)算時(shí)間在特征值問(wèn)題求解上.

表1 不同方法的電位計(jì)算最大誤差和相對(duì)均方根誤差對(duì)比Tab.1 Maximum and mean square error′s comparison between different methods for potential calculation

圖3 算例1解析解和SBFEM數(shù)值解電位等值線Fig.3 Potential isolines of analytical and SBFEM solutions（Example 1）

比例邊界有限元法對(duì)處理復(fù)雜邊界條件問(wèn)題計(jì)算精度也很高，為此本文選擇如下算例.

算例2 有一二維靜電場(chǎng)邊值問(wèn)題，邊界由一段x＝a＝1.0 m的直線和x＝y(tǒng)2的拋物線圍成，其定解問(wèn)題及比例邊界有限元示意圖如圖4所示.網(wǎng)格劃分時(shí)，為了盡可能準(zhǔn)確地描述拋物線邊界形狀，需要采用較多邊界節(jié)點(diǎn)進(jìn)行離散，為此選取在拋物線段上以x＝0.05 m為間距，在直線段上以y＝0.05 m為間距，得到80個(gè)邊界節(jié)點(diǎn).表2為本文方法與文獻(xiàn)[15]列出的半解析方法——多極理論、邊界元法計(jì)算所得電位、電場(chǎng)強(qiáng)度對(duì)比結(jié)果.可以看出比例邊界有限元法計(jì)算復(fù)雜邊界也具有很高的精度.

算例3 本算例為求一無(wú)限域問(wèn)題.一無(wú)限長(zhǎng)、半徑R＝1.0 m圓筒被沿軸線切成兩半，上一半電位為U0＝1 V，下一半接地（電位為0），如圖5（a）所示.由于該問(wèn)題是關(guān)于y軸對(duì)稱問(wèn)題，選取右半部分進(jìn)行模擬.模擬中采用一有限子域和一無(wú)限子域，它們的相似中心放在同一點(diǎn)，離散邊界選取33個(gè)節(jié)點(diǎn)，如圖5（b）所示.筒內(nèi)和筒外電位解析解和SBFEM數(shù)值解等值線分布如圖6所示.從圖中可以看出：比例邊界有限元法數(shù)值解和解析解也非常吻合.

圖4 比例邊界有限元法算例2計(jì)算示意圖Fig.4 Sketch for the second example using scaled boundary finite element method

表2 不同方法的電位計(jì)算結(jié)果比較Tab.2 Potential calculation result comparison between different methods

圖5 模型和比例邊界有限元法計(jì)算示意圖Fig.5 Model and sketch for scaled boundary finite element method

圖6 算例3解析解和SBFEM數(shù)值解電位等值線Fig.6 Potential isolines of analytical and SBFEM solutions（Example 3）

4 結(jié) 論

本文提出了靜電場(chǎng)分析的比例邊界有限元法，推導(dǎo)了相應(yīng)的控制方程，并利用特征值問(wèn)題進(jìn)行了求解.由該方法數(shù)值算例的結(jié)果與解析解、半解析解以及其他數(shù)值方法結(jié)果的對(duì)比發(fā)現(xiàn)該方法具有很高的精度和計(jì)算效率.

[1]李 敬.級(jí)數(shù)法求解軸對(duì)稱的靜電場(chǎng)[J].云南師范大學(xué)學(xué)報(bào)，1997，17（3）：43－45

[2]馬西奎.最小二乘邊界配點(diǎn)法在電磁場(chǎng)邊值問(wèn)題數(shù)值分析中的應(yīng)用[J].微波學(xué)報(bào)，1994，37（2）：17－22

[3]BALLIST C H.The multiple multipole method in electro－and magnetostatic problems [J]. IEEE Transactions on Magnetics，1983，19（6）：2367－2370

[4]盛劍霓.電磁場(chǎng)與波分析中半解析法的理論方法與應(yīng)用[M].北京：科學(xué)出版社，2006：1－35

[5]SONG B，F(xiàn)U J.Modified indirect boundary element technique and its application to electromagnetic potential problems [J]. IEE Proceeding H：Microwaves， Antennas and Propagation， 1992，139（3）：292－296

[6]金建銘.電磁場(chǎng)有限元方法[M].西安：西安電子科技大學(xué)出版社，1998：51－70

[7]BAMJI S S， BULINSKI A T. Electric field

calculations with the boundary element method [J].

IEEE Transactions on Electrical Insulation，1993，

28（3）：420－424

[8]樊德森.靜電場(chǎng)邊值問(wèn)題的矩量法解[J].計(jì)算物理，

1989，6（1）：1－8

[9]張淮清.電磁場(chǎng)計(jì)算中的徑向基函數(shù)無(wú)網(wǎng)格法研究[D].重慶：重慶大學(xué)，2008：25－36

[10]王江忠，趙 良，劉之方.二維開(kāi)域靜電場(chǎng)有限元與邊界元迭代解法的研究[J].華北電力大學(xué)學(xué)報(bào)，2002，29（增刊）：36－40

[11]王志華.有限元與無(wú)單元耦合法在電磁場(chǎng)數(shù)值計(jì)算中的應(yīng)用研究[D].石家莊：河北工業(yè)大學(xué)，2006：1－25

[12]李茂軍.基于邊界元法與無(wú)網(wǎng)格局部Petrov－Galerkin法的耦合法和區(qū)域分解法[D].重慶：重慶大學(xué)，2009：1－36

[13]WOLF J P，SONG C M.Consistent infinitesimal finite element cell method：three－dimensional vector wave equation [J].International Journal for Numerical Methods in Engineering，1996，39（13）：2189－2208

[14]LIU J，LIN G，WANG F M，etal.The scaled boundary finite element method applied to electromagnetic field problems [C] // IOP Conference Series：Materials Science and Engineering，Syndey，WCCM/APCOM 2010.London：IOP，2010

[15]鄭勤紅.計(jì)算復(fù)雜場(chǎng)域靜電場(chǎng)問(wèn)題的多極理論[J].電子科技大學(xué)學(xué)報(bào)，1997，26（6）：599－604