靜電場分析的比例邊界有限元法

劉 俊， 林 皋＊， 王 復 明， 李 建 波

（1.大連理工大學 建設工程學部 水利工程學院，遼寧 大連 116024；2.大連理工大學 海岸和近海工程國家重點實驗室，遼寧 大連 116024；3.鄭州大學 水利與環境學院，河南 鄭州 450002）

0 引 言

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

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

1 靜電場的比例邊界有限元方程推導

描述靜電場（域內無電荷）特性的控制方程為拉普拉斯方程（ 為電位）：

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

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

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

兩坐標系統可通過雅克比矩陣聯系起來：

其中

在比例坐標下梯度的算子可以表示為

其中

對方程（1）、（2）、（3）利用加權余量法和格林第一公式可得（w為權函數）

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

方程（11）中權函數w也用同樣的形函數表達：

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

其中

比例坐標下的空間微元面積

比例坐標下的邊界微元線長度

其中τξ為關于ξ的函數.

對式（14）中含w（ξ，s），ξ項做分部積分，并將式（17）、（18）代入可得

方程（19）的系數為

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

方程（23）為比例邊界有限元的基本方程.方程（24）、（25）為需要滿足的計算域內、外邊界條件.對于有限區域（0＝ξ0≤ξ≤ξ1＝1），內邊界條件式（24）為相似中心一點，外邊界條件式（25）為離散邊界；對于無限區域（1＝ξ0≤ξ≤ξ1＝∞），內邊界條件式（24）為離散邊界，外邊界條件式（25）為無窮遠處外邊界條件.

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

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

由此可得具有兩倍未知數的變量：

可將式（23）化為狀態方程：

其中

由于Z陣為Hamilton陣，可以通過求解Z的特征值問題來得到式（30）的解（式（31）特征值中成對出現，互為負號）：

由此方程（30）的解為

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

式中：c1、c2為積分常數.對于有限域問題，ξ＝0處的 取有限值，故c2＝0；對于無限域問題，ξ＝∞處的 取有限值，故c1＝0.有限域的積分常數c1和無限域的積分常數c2均可由邊界條件確定.積分常數確定后，可由式（32）、（12）確定域內任意點的電位和由E＝－確定域內任意點電場強度.

3 算例計算分析

為了驗證該方法數值模擬精度和高效性，本文首先選擇了一個具有解析解和其他數值解的經典算例以便于對比.

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

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

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

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

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

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

比例邊界有限元法對處理復雜邊界條件問題計算精度也很高，為此本文選擇如下算例.

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

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

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

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

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

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

4 結 論

本文提出了靜電場分析的比例邊界有限元法，推導了相應的控制方程，并利用特征值問題進行了求解.由該方法數值算例的結果與解析解、半解析解以及其他數值方法結果的對比發現該方法具有很高的精度和計算效率.

[1]李 敬.級數法求解軸對稱的靜電場[J].云南師范大學學報，1997，17（3）：43－45

[2]馬西奎.最小二乘邊界配點法在電磁場邊值問題數值分析中的應用[J].微波學報，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]盛劍霓.電磁場與波分析中半解析法的理論方法與應用[M].北京：科學出版社，2006：1－35

[5]SONG B，FU 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]金建銘.電磁場有限元方法[M].西安：西安電子科技大學出版社，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]樊德森.靜電場邊值問題的矩量法解[J].計算物理，

1989，6（1）：1－8

[9]張淮清.電磁場計算中的徑向基函數無網格法研究[D].重慶：重慶大學，2008：25－36

[10]王江忠，趙 良，劉之方.二維開域靜電場有限元與邊界元迭代解法的研究[J].華北電力大學學報，2002，29（增刊）：36－40

[11]王志華.有限元與無單元耦合法在電磁場數值計算中的應用研究[D].石家莊：河北工業大學，2006：1－25

[12]李茂軍.基于邊界元法與無網格局部Petrov－Galerkin法的耦合法和區域分解法[D].重慶：重慶大學，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].電子科技大學學報，1997，26（6）：599－604