等幾何分析方法求解靜電場非齊次邊值問題

張 勇 林 皋 胡志強 徐喜榮

(1.大連理工大學建設工程學部，遼寧 大連 116024； 2.大連理工大學電子信息與電氣工程學部，遼寧 大連 116024)

引 言

計算電磁學是電磁學分析的一種重要工具，其中對靜電場求解是最為基礎和重要的環節。計算電磁學的求解方法主要有解析法、半解析法和數值解法三大類。解析法和半解析法雖然能夠得到精確或半精確的解析式，但只適用于簡單的求解域、簡單的物理參數的問題，不具有通用性。數值解法雖然是在一定條件下的近似解法，但具有較強的通用性，最為常見的有有限差分法、有限元法、邊界元法、矩量法、無網格法[1-8]等。其中，有限元法是應用最為廣泛的數值方法。各類數值方法都需要對求解域進行二次離散或者建模，而離散之后的計算模型可能與原模型不具有一致性，從而喪失精度。

等幾何分析方法(IGA)是2005年由Hughes[10]率先提出的一種新型數值方法。該方法避免了傳統數值方法中求解模型與設計模型的非一致性，從而實現了將問題的分析計算構架于精確模型上，提高了計算精度與效率。目前這種方法已被用于求解固體力學和流體力學問題[11-16]。最近，張勇和林皋等將該方法應用于靜電場問題[17]，利用加權余量法推導了靜電場齊次邊值問題的等幾何分析方程，并實現了對自然邊界條件和齊次強制邊界條件的施加。另外，也將其應用于波導本征問題[18-19]的求解, 但在波導本征問題中，無論是橫向磁場(TM)波還是橫向電場(TE)波，邊界條件都是齊次的。雖然等幾何分析方法具有良好的計算特性和較高的計算效率，但是由于非均勻有理B樣條(NURBS)基函數不具有插值性[20]，使得其求解帶有非齊次強制邊界條件的邊值問題變得復雜和困難。

自然邊界條件和強制邊界條件的施加方式是不同的，自然邊界條件最終轉化為關于邊界積分，從而對齊次或者非齊次邊界并不敏感。但對于強制邊界條件，需要給定離散自由度值，這對于具有插值特性的數值方法而言，是極為容易的，比如在有限元方法中，只需通過結點坐標便可計算出邊界結點自由度值。但對于不具有插值特性的數值方法而言，這樣是無法做到的，比如在無網格方法，張淮清等[8]給出了基于拉格朗日乘子方法施加強制邊界條件，文獻[9]則給出了拉格朗日乘子、罰函數法、混合法等施加邊界條件，其中拉格朗日乘子最為簡潔。本文在文獻[17]基礎上將等幾何分析方法進一步拓展到靜電場非齊次邊值問題。建立帶有拉格朗日乘子的泛函，同時對拉格朗日乘子采用等幾何離散，最終得到關于控制點電勢及與非齊次邊界相關的控制點拉格朗日乘子的增廣代數方程。考慮到拉格朗日乘子僅賦存于非齊次強制邊界上，其增加的計算自由度相對于整體自由度而言較小，從而并不會影響等幾何分析方法的計算效率。

應用基于拉格朗日乘子的等幾何分析方法求解靜電場非齊次邊值問題，以頂部帶有正弦分布電勢的方形板以及內外壁均為非齊次邊值的圓環域為數值算例。計算結果表明本文方法能很好處理靜電場非齊次邊值問題，通過與傳統方法和解析解的比較，顯示出其低自由度消耗、高效和高精度的特點。

1 基于拉格朗日乘子的靜電場等幾何分析方法推導

1.1 靜電場非齊次邊值問題

靜電場問題是電磁學分析中的一個典型邊值問題(BVP)，描述靜電場(域內無電荷)特性的控制方程(φ為電勢)為

1.2 曲面的NURBS描述

NURBS是計算機輔助設計(CAD)中幾何形體的標準數學描述，是由NURBS基函數Ri,p,(i=0,1,…,n)，控制點Pi,(i=0,1,…,n)和節點向量Ξ={ξ0,ξ1…,ξm}三要素來定義的，其中，p為基函數次數，n為控制點數，m為節點數。NURBS基函數是通過B-樣條有理化得到的，如一維NURBS基函數可表示為

(3)

式中:ωi是對應于Ni,p(ξ)的權值；Ni,p(ξ)為B-樣條基函數，可由Cox-de Boor公式[20]遞歸定義為式(4)：

(4)

NURBS基函數將節點向量所在的參數空間中的點映射到控制點所在的物理空間，即NURBS參數空間的一個參數區間[ξi,ξi+1)與物理空間中的一個單元Vi相對應。

NURBS具有與B樣條一樣的性質，如非負性、單位分解性、緊支性、靈活連續性、線性無關性和非插值性等。但與B樣條相比，它可以更準確地描述圓錐曲線和曲面。當權值相等時NURBS基函數則退化為B樣條基函數。

高維NURBS基函數可由低維NURBS基函數張量積得到，如式(5)所示的二維NURBS基函數形式，p、q分別為兩個方向基函數的次數。

(5)

一般地，NURBS曲面的向量形式可寫為

(6)

式中，Pi,j是NURBS曲面控制點。圖1給出了典型的NURBS曲面，圖1(a)為空間NURBS曲面，圖1(b)為平面內NURBS曲面，后者適用于二維靜電場問題。

(a) 二次NURBS曲面

(b) 二次NURBS平面內曲面圖1 NURBS曲面(控制點標記為·)

1.3 求解域的等幾何離散與細分

在等幾何分析方法中，求解域的離散模型不是二次建模得到的，而是繼承于求解域的設計模型，即 CAD系統中的NURBS描述，如上節所述，對于二維問題，兩個張量積方向的節點向量可表示為Ξ={ξ0,ξ1，…,ξku-1},Η={η0,η1，…,ηkv-1}.參數域內的節點區間[ξh,ξh+1)×[ηk,ηk+1)映射成求解域上的一個單元Vh,k。求解域NURBS單元Vh,k的任一點坐標表示為

=NBh,k(ξ,η)xh,k

(7)

從CAD系統導入的初始模型一般并不能滿足解的精度需求，需要對求解域進行細分。而在等幾何分析方法中，求解域的細分是非通信的、保形的，即細分的過程無需再與原來的設計模型進行通信，并且細分前后的模型是幾何一致的，這是等幾何分析的重要特征，也是其與傳統有限元方法的重要區別。保形細分算法有h-細化方法，p-細化方法，k-細化方法等[10]。

1.4 電勢場的等參化

等參方法用于近似求解域的電勢場，求解域中任一點的電勢采用和式(7)相似的表示形式為：

(8)

式中，φi表示求解域控制點的電勢場變量。

1.5 等幾何分析方法離散方程

采用虛功原理[21]將控制微分方程式(1)及非齊次邊界條件式(2)表達成帶有拉格朗日乘子的泛函

(9)

對泛函Π(φ)取一次變分得

=0

(10)

式(10)可進一步離散化為

=0

(11)

對引入的拉格朗日乘子采用NURBS等幾何離散，即

=NB(ξ,η)λ

(12)

將式(8)和(12)代入式(11)，注意到δφ和δλ的任意性，得到等幾何分析方法離散方程的形式為

(13)

(14)

(15)

(16)

2 靜電場非齊次邊值問題的算例分析

2.1 方形槽

通過分離變量法，可得該問題的電勢及電場強度解析解如式(17)所示：

(17)

數值計算時，上述問題簡化為平面問題，如圖2所示，其中紅色點為等幾何分析方法的初始網格控制點。

圖2 方形槽非齊次邊值問題及IGA初始網格示意圖

等幾何分析方法的初始網格只有一個等幾何單元，不能滿足計算的精度要求，利用h-細化方法，可以得到較細的網格，即網格A(含有36個控制點自由度和6個拉格朗日乘子自由度)和網格B(含有100個控制點自由度和10個拉格朗日乘子自由度)。表1為本文方法與文獻[8]列出的有限差分法(FDM)、徑向基無網格法(RBF)計算電勢最大值誤差、相對均方根誤差對比。本文方法中場值計算點為65×65個，在方形槽內均勻分布。相對均方根誤差ERMS和最大誤差Emax計算公式為：

(18)

Emax=max(abs(φi_exact-φi_calc))

(19)

式中：φi_exact為解析解；φi_calc為數值解。

表1數據表明：網格A的結果和有限差分法的精度相當，網格 B的結果和無網格法接近，這表明本文采用拉格朗日乘子處理等幾何分析方法中的非齊次邊值問題是有效的。對于電場強度，由于文獻[8]中沒有給出計算結果和誤差，故直接給出本文方法計算的電場強度相對均方根誤差和最大誤差，見表2.與電勢數值解的精度相比，電場強度數值解的精度稍低，這是源于電場強度是電勢的梯度場相關，也就是它的數值精度要比電勢低。但網格A到網格B的誤差減小的程度表明，本文方法的數值解有較快收斂速度。表1和表2中的Emax和ERMS就是最大誤差和相對均方根誤差，計算式分別為式(18)和(19)。

表1 電勢計算值的最大值誤差和均方根誤差列表

表2 電場強度計算值的最大值誤差和均方根誤差列表

上述比較初步確定了本文方法的有效性，為了進一步比較其收斂性和收斂速度，采用五種尺寸的單元劃分，對于每一種尺寸的網格，分別對應一個等幾何單元模型和一個傳統的有限元單元模型。計算自由度數量和單元數量在表3中列出，其中，對于等幾何分析方法，其計算自由度數為N+M，N為控制點自由度數，M為非齊次邊界附件自由度數，即拉格朗日乘子離散自由度數。在圖3中用雙對數坐標表示自由度數和單元數的關系，從圖3可以看出等幾何分析方法的單位單元的自由度消耗要比有限元少，并且隨著網格的加密，這種差距在增大。對于非齊次邊值問題，需要附加自由度，但其所占的比重并不大。

表3 各種網格下的單元數量和自由度數量

圖3 計算自由度與單元數量圖

(a) 電勢相對均方根誤差

(b) x方向電場強度相對均方根誤差

(c) y方向電場強度相對均方根誤差圖4 域內電勢和電場強度均方根誤差收斂圖

對每個網格得出65×65個計算點上電勢和電場強度值，并根據公式(18)和(19)計算其與解析解(17)的相對均方根誤差。圖4給出了這些誤差隨網格細化的收斂，圖4(a)、(b)、(c)分別代表電勢、x方向和y方向電場強度的相對均方根誤差收斂圖。從整體上看，電勢場的誤差收斂較快(在網格5中為10-7數量級)，電場強度收斂較慢(在網格5中為10-4數量級)；從電勢或者電場強度收斂斜率看，等幾何分析方法的收斂斜率均比有限元大，也就是說等幾何分析方法的收斂速度快。

最后給出整個域上的等幾何分析方法在網格5中的電勢數值解分布圖，如圖5(看1063頁)所示，(a)與(b)分別是數值解及其和解析解的相對誤差分布圖，圖5(b)中的相對誤差在10-6數量級，表明其數值解精度高。另外從誤差在域內的分布可以看出，非齊次邊界上解的精度要比域內解低，這表明對于非齊次邊值，其上拉格朗日乘子的離散數目不能過少。

從這個簡單的靜電場非齊次邊值問題可以看出，等幾何分析方法的優勢是固有的等幾何性，即計算模型與設計模型一致，以及細化非通信和保形特點。但是由于其基函數的非插值性，對于非齊次邊界需要額外的附加自由度，所幸的是邊界自由度數是O(N)，而域內控制點的自由度數是O(N2)，其中N為一個方向上自由度的度量，即附加自由度占計算自由度的比重較小，特別是對于細網格劃分情況。此外，雖然增加了計算自由度，但和傳統方法相比，等幾何分析方法在計算效率上具有單位單元的計算自由度消耗小、計算精度高、收斂速度快的特點。

2.2 同心圓柱環面

電位和電場強度的解析解可由極坐標下的分離變量法求得

k(Crk+Dr-k)sin(kθ)eθ

(20)

式中系數A、B、C、D可根據邊界條件計算得出。在計算中，取Ra=0.5 m，Rb=1.0 m，φa1=1.75 V，φa2=0.5 V，φb1=-0.5 V，φb2=-0.5 V，k=2.

圖6 圓柱環面非齊次邊值問題及初始網格圖

等幾何分析的初始模型如圖6所示，紅色點表示等幾何模型的控制點，其不在域內體現了控制點的非插值性。對于圓柱環面非齊次邊值問題，其內壁和外壁均為非齊次邊界條件，比算例1的齊次和非齊次混合邊值條件要簡單，但本算例模型的邊界是更為復雜的曲線形式。采用六種尺寸的單元劃分，對于每一種尺寸的網格，分別對應一個等幾何單元模型和一個傳統的有限元單元模型，各種網格下的單元數量和自由度數量如表4所列，對于網格6，有限元模型的自由度數已經超過10 000，而等幾何分析方法只有不到5 000.

表4 各種網格下的單元數量和自由度數量

對每個網格計算出260(環向)×65(徑向)個計算點上電勢和電場強度值，并根據式(18)和(19)計算其與解析解(20)的相對均方根誤差。在圖7(a)、(b)、(c)中分別顯示了電勢和電場強度的相對均方根誤差隨計算自由度的收斂性，從中可以看出在網格較粗(網格1)的情況下，由于非齊次邊界上附加自由度占的比重較大，使得等幾何分析方法比有限元要消耗更多的自由度數，但從圖中誤差-自由度的斜率可以看出，隨著網格加密，等幾何分析方法的收斂較快。需指出是圖7為雙對數坐標軸，雖然兩條線的距離不是很大，但體現的差別為倍數關系。

圖8(看1064頁)顯示了等幾何分析方法在網格6中的電勢數值解分布，(a)與(b)分別是數值解及其和解析解的相對誤差分布圖，圖8(b)中的相對誤差在10-5數量級，與算例1相比，其精度有所下降，這說明求解域的形狀及非齊次邊界的數量對結果的精度有一定的影響。

(a) 電勢相對均方根誤差

(b) x方向電場強度相對均方根誤差

(c) y方向電場強度相對均方根誤差圖7 域內電勢和電場強度均方根誤差收斂圖

3 結 論

本文將等幾何分析方法擴展到靜電場非齊次邊值問題，在虛功方程中引入拉格朗日乘子，推導出靜電場非齊次邊值問題的等幾何分析方法求解式。修正后的等幾何分析方法雖然增加了計算自由度(非齊次邊界上的離散拉格朗日乘子)，但與有限元方法相比，該方法仍然具有計算自由度少、收斂速度快、計算精度和計算效率高的特點，并且非通訊保形細化方法使得模型的細分可以自動進行，更為簡潔與方便。應用此方法求解了方形槽與同心圓環的靜電場非齊次邊值問題，結果表明該方法能夠有效的求解非齊次邊值問題，且具有上述的優點，可進一步在計算電磁學中推廣應用。

[1] 周 勤,茅乃豐.旋轉對稱靜電場的一種新的數值解法[J].電波科學學報,1990,5(4):26-34.

ZHOU Qin,MAO Naifeng.The electrostatic field problems with representations of spline functions and their quasi-charge density analysis method[J].Chinese Journal of Radio Science,1990,5(4):26-34.(in Chinese)

[2] SONG B,FU J.Modified indirect boundary element technique and its application to electromagnetic potential problems[J].IEE proceeding-H of Microwaves Antennas and Propagation,1992,139(3):292-296.

[3] 金建銘.電磁場有限元方法[M].西安電子科技大學出版社,1998:51-70.

[4] 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.

[5] 甄蜀春,曹 蕾,張繼龍.波動方程差分解法對波導螺釘調配器的分析[J].電波科學學報,2005,20(1):125-127.

ZHEN Shuchun,CAO Lei,ZHANG Jilong.Analysis of screw tuner based on wave equation finite difference method[J].Chinese Journal of Radio Science,2005,20(1):125-127.( in Chinese)

[6] 汪朝暉,廖振方,陳德淑.有限元法分析尖板電極結構的空間靜電場分布[J].重慶大學學報,2010,33(5):41-47.

WANG Zhaohui,LIAO Zhenfang,CHEN Deshu.Analysis of spatial electric field with point-plate electrodes configuration using finite element method[J].Journal of Chongqing University,2010,33(5):41-47.(in Chinese)

[7] 梁志偉,趙國偉,徐 杰,等.柱形等離子體天線輻射特性的矩量法分析[J].電波科學學報,2008,23(4):749-753.

LIANG Zhiwei,ZHAO Guowei,XU Jie,et al.Analysis of plasma-column antenna using moment method[J].Chinese Journal of Radio Science,2008,23(4):749-753.(in Chinese)

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

ZHANG Huaiqing.Research on Radial Basis Function Meshless Method in Numercial Computation of Electromagentic Field[D].Chongqing： Chongqing University,2008.(in Chinese)

[9] FERNANDEZ-MENDEZ S,HUERTA A.Imposing essential boundary conditions in mesh-free methods[J].Comput Methods Appl Mech Engrg,2004,193(12/14):1257-1275.

[10] HUGHES T J R,COTTRELL J A and BAZILEVS Y.Isogeometric analysis CAD,finite elements,NURBS,exact geometry and mesh refinement[J].Comput Methods Appl Mech Engrg,2005,194(39/41):4135/4195.

[11] REALI A.An isogeometric analysis approach for the study of structural vibrations[J].Journal of Earthquake Engineering,2006,10(Special Issue 1):1-30.

[12] BAZILEVS Y,CALO V M,ZHANG Y,et al.Isogeometric fluid-structure interaction analysis with applications to arterial blood flow[J].Computational Mechanics,2006,38(4/5):310-322.

[13] COTTRELL J A,REALI A,BAZILEVS Y,et al.Isogeometric analysis of structural vibrations[J].Computer Methods in Applied Mechanics and Engineering,2006,195(41/43):5257-5296.

[14] COTTRELL J A,HUGHES T J R,REALI A.Studies of refinement and continuity in isogeometric structural analysis[J].Computer Methods in Applied Mechanics and Engineering,2007,196(41/44):4160-4183.

[15] ZHANG Y,BAZILEVS Y,GOSWAMI S,et al.Patient-specific vascular NURBS modeling for isogeometric analysis of blood flow[J].Computer Methods in Applied Mechanics and Engineering,2007,196(31/32):2943-2959.

[16] ZHANG Y,LIN G,HU Z Q.Isogeometric analysis based on scaled boundary finite element method[J].IOP Conference Series:Materials Science and Engineering,2010,10(1):012237.

[17] 張 勇,林 皋,劉 俊,等.等幾何分析方法及其應用于偏心柱面靜電場問題[J].電波科學學報,2012,27(1):177-183.

ZHANG Yong,LIN Gao,LIU Jun,et al.The isogeometric analysis and its application to eccentrically parallel cylinders electrostatic problem[J].Chinese Journal of Radio Science,2012,27(1):177-183.(in Chinese)

[18] 張 勇,林 皋,劉 俊,等.波導本征問題的等幾何分析方法[J].應用力學學報,2012,29(2):113-119.

ZHANG Yong,LIN Gao,LIU Jun,et al.The isogeometric analysis for eigen problem of waveguide[J].Chinese Journal of Applied Mechanics,2012,29(2):113-119.(in Chinese)

[19] 張 勇,林 皋,胡志強,等.基于等幾何分析方法求解任意截面波導本征問題[J].計算物理,2012,29(4):9-17.

ZHANG Yong,LIN Gao,HU Zhiqiang,et al.Eigenvalue Analysis of Waveguide with Arbitrary Cross-section by Isogeometric Analysis[J].Chinese Journal of Computational Physics,2012,29(4):9-17.(in Chinese)

[20] PIEGL L,TILLER W.The NURBS Book [M].Berlin:Springer Verlag,1997.

[21] HUGHES T J R.The Finite Element Method,Linear Static and Dynamic Finite Element Analysis[M].New York:Dover,2000.