比例邊界有限元法分析靜電場(chǎng)問(wèn)題

姜宇馳

（常熟理工學(xué)院 物理與電子工程學(xué)院，江蘇 常熟 2 15500）

1 引言

在許多電工設(shè)備、絕緣結(jié)構(gòu)以及微波領(lǐng)域的TEM傳輸線特性分析中，常常需要確定靜電場(chǎng)分布.對(duì)于靜電場(chǎng)問(wèn)題求解常見(jiàn)的數(shù)值方法有：有限元[1]、有限差分法[2]和邊界元法[3]等.有限元和有限差分法都屬于全域算法，需要對(duì)整個(gè)求解域離散，數(shù)據(jù)準(zhǔn)備和求解的工作量較大；另外，在靜電場(chǎng)分析中常常遇到開(kāi)域問(wèn)題，更增加了問(wèn)題的處理難度.如文獻(xiàn)[4]采用相似剖分有限元法，利用遞推關(guān)系引入轉(zhuǎn)移矩陣來(lái)構(gòu)造包含整個(gè)無(wú)限空間的有限元節(jié)點(diǎn)信息的等效剛度陣處理開(kāi)域問(wèn)題.文獻(xiàn)[5]討論了如何設(shè)置截?cái)噙吔鐥l件來(lái)獲得高精度的有限元解.邊界元法，只需對(duì)求解域邊界上劃分單元，克服了有限元法和有限差分的全域離散的不足，但存在尋找基本解和奇異積分計(jì)算的問(wèn)題.由Song和Wolf[6]在上世紀(jì)90年代率先提出和發(fā)展起來(lái)的比例邊界有限元法，不僅集成了有限元法和邊界元法的優(yōu)點(diǎn)，而且具有獨(dú)特的優(yōu)勢(shì)：首先，和邊界元法一樣，此法僅離散計(jì)算域的邊界，空間上降低一個(gè)維度，并且它不需要基本解；其次，與比例中心相連的邊界和材料界面不需要離散，這有助于降低計(jì)算成本；最后，解在徑向是解析的，方便處理開(kāi)域問(wèn)題.比例邊界有限元法的這些顯著優(yōu)點(diǎn)，已被最近報(bào)道的若干研究證實(shí).如文獻(xiàn)[7]使用比例邊界有限元法計(jì)算了雙材料界面裂縫瞬態(tài)動(dòng)應(yīng)力強(qiáng)度因子.文獻(xiàn)[8]利用比例邊界有限元法模擬了自由面滲流問(wèn)題.但是，比例邊界有限元法在靜電場(chǎng)分析中的應(yīng)用還很少見(jiàn)報(bào)道.本文根據(jù)在實(shí)際靜電場(chǎng)問(wèn)題求解中，經(jīng)常遇到在某直線邊界上給定均一電位的情況，如果比例中心位于這樣的邊界上，該邊界可以不進(jìn)行離散，推導(dǎo)了比例中心位于這類(lèi)邊界上的比例邊界有限元處理公式，并將其用到有限域和無(wú)限開(kāi)域邊界的處理中.通過(guò)與解析解的比較，說(shuō)明了比例邊界有限元法求解靜電場(chǎng)問(wèn)題不僅具有滿意的精度，方便處理開(kāi)域問(wèn)題，而且能極大減少數(shù)據(jù)準(zhǔn)備工作量.

2 靜電場(chǎng)問(wèn)題的數(shù)學(xué)模型及比例邊界有限元方程

2.1 靜電場(chǎng)問(wèn)題數(shù)學(xué)模型

為了討論方便，考慮二維平面靜電場(chǎng)問(wèn)題，它滿足的拉普拉斯方程和邊界條件為

式中?為梯度算子，φ為電位，n為邊界外法線方向，S1為給定電位邊界，也稱(chēng)為第一類(lèi)邊界，S2為第二類(lèi)邊界.

2.2 靜電場(chǎng)比例邊界有限元方程

圖1 比例邊界坐標(biāo)定義

為進(jìn)行比例邊界有限元分析，考察圖1所示的內(nèi)外邊界相對(duì)于比例中心O具有相似關(guān)系的控制域.從比例中心出發(fā)，定義徑向比例坐標(biāo)ξ，在比例中心點(diǎn)O，ξ=0，在內(nèi)邊界：ξ=ξ0，而在外邊界曲線S上坐 標(biāo) ξ=ξ1. 整 個(gè) 求 解 域 定 義 為 ξ0≤ξ≤ξ1，s0≤s≤s1，其中s表示環(huán)向坐標(biāo)，如果s0，s1位置重合，則屬于封閉域問(wèn)題，否則，屬于非封閉域問(wèn)題；如果 ξ0=0，ξ1=1表示有限區(qū)域問(wèn)題，如果 ξ0=1，ξ1=∞表示開(kāi)域問(wèn)題.為討論方便，比例中心O點(diǎn)設(shè)為笛卡兒坐標(biāo)系原點(diǎn).將邊界曲線S劃分成若干個(gè)有限單元，曲線S上的點(diǎn)位置可以利用插值函數(shù)[ ]N(s)近似表示為

其中，{x}，{y}為邊界曲線上單元結(jié)點(diǎn)坐標(biāo).與邊界上的點(diǎn)相對(duì)應(yīng)的在比例坐標(biāo)為ξ的曲線S（ξ）上的點(diǎn)坐標(biāo)可以表示為

利用式（4）和（5），計(jì)算域內(nèi)任意一點(diǎn)的位置可由比例坐標(biāo)系中的坐標(biāo)分量ξ和s確定.

笛卡兒坐標(biāo)系下的梯度算子在比例坐標(biāo)系下可變換為

其中

雅克比行列式定義為

利用加權(quán)余量法方程（1）可以寫(xiě)成

其中W為任意權(quán)函數(shù).

引進(jìn)有限元方法中的插值函數(shù)[N ( s)]，勢(shì)函數(shù)φ(ξ , s)可表示為

式中{a ( ξ) } 為結(jié)點(diǎn)勢(shì)向量.應(yīng)用伽遼金法，權(quán)函數(shù)可采用勢(shì)函數(shù)同樣的插值函數(shù)

比例坐標(biāo)系中，dΩ=J ξdξds，且ξ為常數(shù)時(shí)有dS=τξds，其中τξ是比例坐標(biāo)與總體坐標(biāo)間的比例因子；在 s為常數(shù)時(shí)有ds=τsdξ.將（6）、（10）和（11）代入（9）式可得[9]

方程（14）即為所求的比例邊界有限元方程，（12）（13）為內(nèi)、外邊界條件.其中

而

與比例坐標(biāo)ξ無(wú)關(guān)，只取決于s.

這里[E0]、[E1]、[E2]盡管是對(duì)曲線S（ξ）上的單元給出的，由于比例坐標(biāo)的性質(zhì)對(duì)于任意ξ曲線上的單元，它們恒等于所對(duì)應(yīng)邊界單元上的值.因此在對(duì)區(qū)域離散時(shí)僅離散邊界單元即可.而{Fs(ξ)}表示側(cè)面邊界的貢獻(xiàn)，當(dāng)域封閉即不存在側(cè)面邊界或側(cè)面為絕緣邊界時(shí)該項(xiàng)為0.

3 方程求解

3.1 比例中心位于非指定均一電位邊界上

從式（12）和（13）可以看出，對(duì)任意比例坐標(biāo) ξ可以將穿過(guò)曲面S（ξ）的外法向通量{Q ( ξ)}寫(xiě)成下面的形式：

其中“±“分別表示有限域（ξ坐標(biāo)方向與表面外法向一致）和開(kāi)域問(wèn)題（ξ坐標(biāo)方向與表面外法向相反）.式（14）是二階Euler—Cauchy方程，這里只考慮側(cè)面為絕緣邊界，即{Fs(ξ)}=0的情況，利用（20）式，將其降為一階，構(gòu)造如下方程求解[9]

其中

可以利用矩陣[Z]的特征值和特征向量求（21）式的解，我們求解如下特征值問(wèn)題

其中特征值[λ1]的實(shí)部為負(fù)值，[λ2]的實(shí)部為正值.（21）式的通解為[9]

其中{c1}、{c2}為常數(shù)，對(duì)于有限域問(wèn)題 ξ=0位置的φ取有限值，故{c2}＝0，{Q ( ξ) } 前面取“＋”號(hào)，而對(duì)于開(kāi)域問(wèn)題，{c1}＝0，{Q ( ξ) } 前面取“－”號(hào).由（25）（26）式，我們可以建立起勢(shì){a ( ξ) } 和法向通量{Q ( ξ) } 之間的關(guān)系：

對(duì)于有限域問(wèn)題：

對(duì)于開(kāi)域問(wèn)題：

這樣我們就可以通過(guò)（27）式獲得邊界電位值，通過(guò)（25）式確定{c1}或{c2}，然后確定計(jì)算域內(nèi)任意點(diǎn)電位值.

3.2 比例中心位于指定均一電位直線邊界

當(dāng)比例中心選擇在指定均一電位直線邊界時(shí)，由于比例邊界有限元在徑向解是解析的，故可以不離散這一邊界，減少數(shù)據(jù)準(zhǔn)備和計(jì)算工作量，但需對(duì)該邊界對(duì)其他單元的影響做特殊處理.為了簡(jiǎn)化，這里也僅考慮絕緣邊界，將方程（14）改寫(xiě)為

其中，下標(biāo)c表示指定均一電位側(cè)邊界所對(duì)應(yīng)的邊界節(jié)點(diǎn)，u是其余邊界節(jié)點(diǎn).由于ac(ξ)=ac，即整個(gè)側(cè)面指定電位為定值與ξ無(wú)關(guān)，則（30）可以改寫(xiě)為

同樣，（20）式也可寫(xiě)成

由于ac(ξ) =ac，則（32）式可變?yōu)?/p>

將（33）式對(duì) ξ求導(dǎo)，然后乘以ξ得

將（34）式代入（31）式，得

將（33）式變形為

將（36）式代入（35）式，得

聯(lián)合（36）、（37）式，寫(xiě)成如下形式：

其中

類(lèi)似地可以利用矩陣[Z]的特征值和特征向量求（38）式的解，我們求解如下特征值問(wèn)題

其中特征值[λ1]的實(shí)部為負(fù)值，[λ2]的實(shí)部為正值.（38）式的解為

其中{c1}、{c2}為常數(shù)，對(duì)于有限域問(wèn)題 ξ=0位勢(shì)φ取有限值，故{c2}＝0，{Q ( ξ) } 前面取“＋”號(hào)，而對(duì)于開(kāi)域問(wèn)題，{c1}＝0，{Q ( ξ) } 前面取“－”號(hào).而

{c1}、{c2}由邊界確定后，域內(nèi)任意坐標(biāo)ξ處的電位和通量可根據(jù)（42）（43）式確定.從而，我們可以建立起電位{a ( ξ) } 和通量{Q ( ξ) } 之間的關(guān)系

對(duì)于有限域問(wèn)題，

對(duì)于開(kāi)域問(wèn)題，

這樣可以通過(guò)（45）式獲得邊界電位值，通過(guò)（42）式確定{c1}或{c2}，然后確定計(jì)算域內(nèi)任意點(diǎn)電位值.

4 算例分析

4.1 比例中心選擇在非指定均一電位邊界

為了驗(yàn)證所推公式的正確性和算法的精度，首先選擇一個(gè)具有解析結(jié)果的經(jīng)典算例.如圖2所示，在半徑為R的圓柱形截面邊界處，指定φ=U0cosθ，求圓柱內(nèi)外靜電場(chǎng)分布.此問(wèn)題有解析解[10]

這里假定U0=10，R=1.0，由于邊界條件的對(duì)稱(chēng)性，選取x軸上方半部分作為分析對(duì)象，沿環(huán)向均勻布置101個(gè)節(jié)點(diǎn)，將上半圓劃分100個(gè)單元，比例中心選在原點(diǎn)O處.對(duì)于圓柱內(nèi)采用有限域公式處理，而圓柱外場(chǎng)采用開(kāi)域比例邊界有限元公式處理，分別求得的圓柱內(nèi)外電位分布如圖3所示.

圖2 指定圓柱表面電位求內(nèi)外場(chǎng)問(wèn)題示意圖

圖3 圓柱內(nèi)外電位分布等值線圖

為了比較數(shù)值解和解析解的接近程度，選擇x軸上三個(gè)點(diǎn)的理論解和數(shù)值解及數(shù)值解相對(duì)誤差列于表1中.

從表1可以看出本文的比例邊界有限元法解的精度非常高，在（2.0，0.0）處相對(duì)誤差最大，僅為0.038%.

表1 指定圓柱表面電位問(wèn)題數(shù)值和理論解比較

4.2 比例中心選擇在指定均一電位直線邊界

4.2.1 有限域靜電場(chǎng)分析

為了考察所推導(dǎo)的比例中心選擇在指定均一直線邊界條件下比例邊界有限元公式的正確性和算法精度解析，考解慮[11一]個(gè)長(zhǎng)直接地金屬槽，其三壁電位為0，頂蓋電位為其截面如圖4所示，該問(wèn)題有

計(jì)算中取a=10，利用對(duì)稱(chēng)性，選擇圖4所示的點(diǎn)劃線的左半部分作為計(jì)算分析對(duì)象.比例中心選擇在坐標(biāo)原點(diǎn)，與x軸和y軸重合的兩個(gè)指定電位為零的邊界通過(guò)比例中心，在計(jì)算的時(shí)候不需要離散，只需離散其余邊界.邊界離散和采用本文算法計(jì)算所得的等位線分布如圖5所示.

為了考察算法的精度，選擇x＝5截面上幾個(gè)點(diǎn)的解析解和計(jì)算值進(jìn)行了比較，如表2所示.可以看出比例邊界元法具有很高的精度，表2中最大相對(duì)誤差僅0.45%.

4.2.2 開(kāi)域靜電場(chǎng)分析

為了考察比例邊界有限元法在開(kāi)域靜電場(chǎng)問(wèn)題分析的精度，我們考察如圖6所示的A、B兩塊半無(wú)限大平行板，它們?cè)谠c(diǎn)O處通過(guò)極薄的絕緣材料隔開(kāi)，左端A板電勢(shì)為10 V，右端B板電勢(shì)為0 V，求兩板上半空間的電位分布.該問(wèn)題的電位分布的理論解[12]為：

計(jì)算時(shí)，比例中心選擇在兩板的交點(diǎn)處，由于平板為指定均一電位的直線邊界，且通過(guò)比例中心，故不需要離散，選擇半徑為R＝10的半圓形進(jìn)行離散.沿該半圓形均勻布置51個(gè)節(jié)點(diǎn)，50個(gè)單元，整個(gè)空間被半圓分割成2個(gè)區(qū)域，一個(gè)是半圓內(nèi)部有限區(qū)域，另一個(gè)是半圓外部開(kāi)域.這兩個(gè)區(qū)域可以看成是兩個(gè)子結(jié)構(gòu)，對(duì)這兩個(gè)子結(jié)構(gòu)分別使用比例邊界有限元法，計(jì)算出矩陣[H]和右端項(xiàng)，然后對(duì)號(hào)累加獲得問(wèn)題的解.比例邊界有限元法求得的等位線分布如圖6所示.

圖4 長(zhǎng)直接地金屬槽邊界條件示意圖

圖5 金屬槽截面電位分布

表2 長(zhǎng)接地金屬槽電位問(wèn)題數(shù)值和理論解比較

比例邊界有限元算得的在半徑R=10的圓環(huán)上，在角度θ=π/5，2π/5，3π/5，4π/5上電位分別為2.0000，4.0000，6.0000和8.0000 V和解析解完全一致.可以看出比例邊界有限元法具有相當(dāng)高的精度，對(duì)開(kāi)域處理非常方便.

圖6 具有不同電位的兩塊半無(wú)限大金屬板所形成的空間電位

5 結(jié)論

本文推導(dǎo)了用于靜電場(chǎng)分析的比例邊界有限元法方程.尤其是對(duì)指定均一電位的直線邊界，推導(dǎo)了比例中心位于這類(lèi)邊界上的處理公式，并將其應(yīng)用到有限域和無(wú)限開(kāi)域邊界的處理中.通過(guò)與解析解的比較，說(shuō)明了比例邊界有限元法求解靜電場(chǎng)問(wèn)題不僅可得到滿意的精度，方便處理開(kāi)域問(wèn)題，而且能極大地減少數(shù)據(jù)準(zhǔn)備工作量，且在計(jì)算過(guò)程中發(fā)現(xiàn)：

（1）比例中心的選擇是比例邊界有限元分析的關(guān)鍵.當(dāng)比例邊界有限元的比例中心選擇在絕緣邊界或指定均一電位邊界上，這些邊界不需要離散，只需離散其余邊界.因此可以減少數(shù)據(jù)準(zhǔn)備工作量和計(jì)算時(shí)間.

（2）對(duì)于簡(jiǎn)單的靜電場(chǎng)問(wèn)題，比例邊界有限元法可以直接應(yīng)用.而對(duì)于材料特性分布比較復(fù)雜或無(wú)窮域靜電場(chǎng)問(wèn)題，單純的比例邊界有限元法難以處理，可采用比例邊界有限元－子結(jié)構(gòu)的方法進(jìn)行處理.

[1]崔翔.應(yīng)用有限元方法計(jì)算含有電位懸浮導(dǎo)體的電場(chǎng)分布[J].華北電力學(xué)院學(xué)報(bào)，1995，22（2）：1-7.

[2]汪琛，王保平，童林夙.計(jì)算二維靜電場(chǎng)的非正交有限差分算法[J].計(jì)算物理，1997，14（3）：305-310.

[3]何锃，呂浚潮，戴呈豪.新型快速多極邊界元法求解電荷任意分布的二維靜電場(chǎng)[J].計(jì)算物理，2007，24（4）：433-438.

[4]董興其，安同一.靜電問(wèn)題的相似剖分有限元方法[J].微波學(xué)報(bào)，1996，12（1）：15-21.

[5]馬西奎，韓社教.開(kāi)放域靜態(tài)電磁場(chǎng)問(wèn)題數(shù)值解的漸進(jìn)邊界條件技術(shù)研究[J].中國(guó)科學(xué)（E輯），2003，33（11）：1021-1027.

[6] Song Ch ，Wolf，J P. Consistent infinitesimal finite element cell method：three-dimensional vector wave equation [J]. International Journal for Numerical Methods in Engineering，1996，39：2189-2208.

[7]楊貞軍，D e e k sAJ.基于頻域比例邊界有限元法的雙材料界面裂縫瞬態(tài)動(dòng)應(yīng)力強(qiáng)度因子的計(jì)算[J].中國(guó)科學(xué)（G輯），2008，38（01）：77-8 8.

[8]李鳳志.滲流自由面分析的比例邊界有限元法[J].計(jì)算物理，2009，26（5）：665-670.

[9] Deeks A J，Cheng L. Potential flow around obstacles using the scaled boundary finite-element method [J]. Int J Numer Mech Fluids，2003，41：721-741.

[10]徐永斌，何國(guó)瑜，盧才成，等.工程電磁場(chǎng)基礎(chǔ)[M].北京：北京航空航天大學(xué)出版社，1992：9，165-166.

[11]倪光正，楊仕有，邱捷，等.工程電磁場(chǎng)數(shù)值計(jì)算[M].北京：機(jī)械工業(yè)出版社，2010：121-122.

[12]劉鵬程.電磁場(chǎng)解析方法[M].北京：電子工業(yè)出版社，1995：106-107.