張 勇 林 皋 劉 俊 徐喜榮

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

引 言

由兩根圓柱導(dǎo)體所組成的傳輸線在實際工程中有著極為廣泛的應(yīng)用，由于特殊要求或者制造的誤差，常常呈現(xiàn)偏心狀態(tài)，是柱面或是橢圓形式的，其靜電場邊值問題成為靜電學(xué)中的典型問題。

計算電磁學(xué)是電磁學(xué)分析的一種重要工具，其中對靜電場求解是最為基礎(chǔ)和重要的環(huán)節(jié)。計算電磁學(xué)主要有解析法、半解析法和數(shù)值解法三大類。解析法雖然能夠得到精確的解析式，清晰地看出各物理量間的依賴關(guān)系，但其只適用于簡單的求解域、簡單的物理參數(shù)，不具有通用性。半解析法有級數(shù)法、配置點法、多極理論方法和新型等效源法［1-4］等方法。這些方法在計算電磁學(xué)中取得了較好的效果和應(yīng)用，但該類方法也只適用一些幾何形狀比較規(guī)則的靜電場問題。數(shù)值解法有有限差分法、有限元法、邊界元法、矩量法、無網(wǎng)格法［5-12］等。其中，有限元法是應(yīng)用最為廣泛的數(shù)值方法，是分析靜電場問題最強大、最通用方法之一，適合于含有復(fù)雜媒質(zhì)問題的數(shù)值分析。但該方法在處理含有場值奇異以及不同材料交接點等問題中遇到很大的挑戰(zhàn)，通常的處理方法是在這些點周圍處增加節(jié)點，但這勢必需要更多的前處理和計算工作時間，當(dāng)然也有很多改進的有限元方法，例如提高形函數(shù)或者使用超級單元等。各類數(shù)值方法都需要對求解域進行離散，而離散之后的計算模型可能與原模型不具有一致性，從而喪失精度。

等幾何分析方法（IGA）是2005年由 Hughes［13］率先提出的一種新型數(shù)值方法。該方法剔除了傳統(tǒng)數(shù)值方法的求解域的模型非一致性，從而實現(xiàn)了將問題的分析計算構(gòu)架于精確模型上，實現(xiàn)較高的計算精度與效率。目前，這種方法已被用于求解固體力學(xué)和流體力學(xué)問題［14-19］。本文首次將該方法應(yīng)用于靜電場問題的求解。從靜電場控制方程—拉普拉斯方程出發(fā)，結(jié)合加權(quán)余量法推導(dǎo)了靜電場問題的等幾何分析方程。

應(yīng)用本文方法求解偏心圓柱面靜電場問題，通過與傳統(tǒng)方法和解析解的比較，等幾何分析方法顯示出低自由度消耗、高效性和精確性。

1.B樣條和NURBS

1.1 B樣條和NURBS基函數(shù)

B樣條是由B樣條基函數(shù)Ni，p（i=0，1，…，n），控制點Pi（i=0，1，…，n）和節(jié)點向量Ξ=｛ξ0，ξ1，…，ξm｝定義的。B樣條基函數(shù)將節(jié)點向量所在的參數(shù)空間中的點映射到控制點所在的物理空間，從而B樣條參數(shù)空間的一個參數(shù)區(qū)間［ξi，ξi+1）與物理空間中的一個單元Vi相對應(yīng)。一維節(jié)點向量Ξ=｛ξ0，ξ1…，ξm｝是一組參數(shù)空間中的非遞減的坐標(biāo)序列，這里ξi是節(jié)點，i是節(jié)點標(biāo)號。

一維空間B樣條基函數(shù)Cox-de Boor公式［20］，遞歸定義為式（1），p為基函數(shù)的次數(shù)。

如果節(jié)點在參數(shù)空間中是等間距的則稱為均勻B樣條，否則，稱為非均勻B樣條。節(jié)點重數(shù)是指節(jié)點在節(jié)點向量中出現(xiàn)的次數(shù)，重數(shù)大于1的節(jié)點定義為多重節(jié)點。如果一個節(jié)點向量的第一個和最后一個節(jié)點是p+1重的，那么它是開放的。開放節(jié)點向量是計算機輔助設(shè)計（CAD）標(biāo)準(zhǔn)，本文也采用開放節(jié)點向量。B樣條通常沒有插值性，即控制點并不在B樣條曲線上，即使開放節(jié)點向量的B樣條也僅僅插值于端點。節(jié)點向量為｛0，0，0，1，2，3，4，4，4｝的基函數(shù)如圖1所示。

圖1 二次NURBS基函數(shù)

非均勻有理B樣條 （NURBS）是通過B樣條有理化實現(xiàn)的，關(guān)系式為

式中:ωi是Ni，p（ξ）相應(yīng)的權(quán)值，所以，NURBS具有與B樣條一樣的性質(zhì)，如非負性、單位分解性、緊支性、靈活連續(xù)性和線性無關(guān)性等。與B樣條相比，它可以更準(zhǔn)確地建立圓錐曲線和曲面。若權(quán)值相等，則NURBS基函數(shù)退化為B樣條基函數(shù)。

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

1.2 基于NURBS的幾何形體表示方法

NURBS曲線用向量形式表示為

式中:Pi是NURBS曲線控制點。圖2、圖3給出了一典型的NURBS曲線，采用的基函數(shù)如圖1所示。

對于NURBS曲面也有相似的形式，且為

2.基于NURBS的靜電場等幾何分析方法推導(dǎo)

2.1 靜電場問題

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

2.2 求解域的等幾何離散

求解域采用NURBS離散，節(jié)點向量表示為Ξ=｛ξ0，ξ1，…，ξku-1｝，Η=｛η0，η1，…，ηkv-1｝。參數(shù)域內(nèi)的節(jié)點區(qū)間［ξh，ξh+1）×［ηk，ηk+1）映射成求解域上的一個單元Vh，k，相應(yīng)的非零NURBS基函數(shù)索引為（h1，k1），其中，h1=p-h，p-h+1，…，p，k1=qk，q-k+1，…，q，p、q分別為兩個方向基函數(shù)次數(shù)。根據(jù)NURBS基函數(shù)的緊支性，單元Vh，k的所有非零基函數(shù)個數(shù)為（p+1）（q+1），與該單元的控制點數(shù)量相同。求解域NURBS單元Vh，k的任一點坐標(biāo)表示為

式中:X（ξ，η）h，k=｛X（ξ，η）h，k，Y（ξ，η）h，k）｝表示求解域中 NURBS單元Vh，k任一點坐標(biāo)；=（，）表示相應(yīng)與該單元的控制點坐標(biāo)；xh，k表示該單元所有控制點坐標(biāo)列向量；是等幾何分析基函數(shù)，可表示B樣條基函數(shù)或NURBS基函數(shù)的張量積形式；NBh，k（ξ，η）是相應(yīng)的矩陣形式。

為簡化，以后的符號中省掉單元上標(biāo)h，k.圖4給出了一典型被NURBS離散的求解域，控制點標(biāo)記為·，黑色虛線表示控制網(wǎng)，綠色曲線表示NURBS單元的邊界，控制點的非插值性如圖4所示。

圖4 NURBS離散求解域

2.3 求解域的保形細化

求解域的保形細分是等幾何分析的重要特征，對于初始的模型，可以通過保形細分得到與原模型相一致的細化模型。保形細分算法有h-細化方法，p-細化方法，k-細化方法等［13］。

2.4 電勢場的等參化

等參方法用于近似求解域的電勢場，求解域中的任一點的電勢具有與式（7）相似的表示形式。

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

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

控制微分方程式（6a）是采用Galerkin加權(quán)余量法或者虛功原理［21］求得的。推導(dǎo)過程從略，等幾何分析方法離散方程可表示為

關(guān)于控制點電勢的方程，K=表 示NURBS單元系數(shù)矩陣集成的總體系數(shù)矩陣為

公式（10）的積分區(qū)間是 NURBS單元Vh，k的參數(shù)區(qū)間［ξh，ξh+1）×［ηk，ηk+1），不是非標(biāo)準(zhǔn)的高斯積分區(qū)間，計算積分時需要轉(zhuǎn)化，對于不同的單元這個積分區(qū)間是不同的。

Q=是等效控制點激勵，包括體電荷分布和第二類邊界條件的貢獻，表達式為

2.6 邊界條件施加及離散方程求解

NURBS單元的控制點不具有插值性，邊界條件的施加相對與傳統(tǒng)的有限元有較大的區(qū)別。對于第一類邊界條件，利用NURBS基函數(shù)的非負性將與邊界相關(guān)的單元的控制點分成兩組，在Γφ上恒為零的和不恒為零的。

恒為零的基函數(shù)對應(yīng)的控制點對電勢沒有貢獻，即式（12）只剩下非負基函數(shù)對應(yīng)的項，再利用NURBS基函數(shù)的單位分解性，可得

從而將求解域上的邊界約束條件轉(zhuǎn)化到控制點上。

第二類邊界條件貢獻，即式（11）的第二項，利用邊界積分，將邊界微元轉(zhuǎn)化到參數(shù)域中進行。

對于不同的邊界面，邊界微元形式和積分區(qū)間也不盡相同，例如，對于U參數(shù)ξh邊界面，形如式（15）.

3.偏心柱面靜電場問題算例分析

3.1 偏心圓柱面

偏心圓柱面靜電場問題［22］，即一圓柱面（半徑為R2，軸線為L2）位于另一個圓柱面（半徑為R1，軸線為L1）的內(nèi)部，兩圓柱面上的電勢分別為V2和V1，但兩圓柱面的軸線L1與L2不重合，且軸線間距為L，求解兩圓柱面夾層之間的電勢分布問題。

上述問題簡化為平面問題，如圖5所示，其電勢解析解為

式中:C0、A0、x1、x2為計算系數(shù)，與偏心圓柱面的幾何尺寸R1、R2、L及電場參數(shù)V2和V1有關(guān)。

圖5 偏心圓柱面示意圖

在仿真中，取R1=1.0m，R2=0.5m，L=0.2m，V1=0V，V2=1V.利用求解域與邊界條件的對稱性，只取上半部分進行分析，如圖6所示。

圖6 偏心圓柱面計算簡圖

為了比較解的收斂性和收斂速度，采用五種尺寸的單元劃分，對于每一種尺寸的網(wǎng)格，分別對應(yīng)一個等幾何單元模型和一個傳統(tǒng)的有限元單元模型，計算自由度數(shù)量和單元數(shù)量在表1中列出，并在圖7中顯示。隨著單元數(shù)量的增加，傳統(tǒng)有限元單元模型的計算自由度迅速增加，而等幾何單元模型的計算自由度增加的速率明顯緩慢，這源于等幾何分析單元之間可以共享更多的自由度。在圖6所示的典型截面上電場強度的分布如圖8所示。圖8（a）為有限元計算的結(jié)果，在單元過度處出現(xiàn)震蕩，即使很細的網(wǎng)格也未能幸免，而圖8（b）所示的等幾何分析計算結(jié)果，和解析解的誤差甚小，且不出現(xiàn)震蕩。

表1 各種網(wǎng)格下的單元數(shù)量和自由度數(shù)量

圖7 計算自由度與單元數(shù)量圖

整個求解域的相對解析解的誤差計算公式為式（17），φi-exact為解析解，φi-calc為數(shù)值解。

電勢的相對誤差隨網(wǎng)格尺寸變化如圖9所示，隨著單元尺寸的減少，傳統(tǒng)有限元Lagrange單元和NRUBS等幾何單元的計算結(jié)果的相對誤差都減少，而在每個網(wǎng)格尺寸下后者比前者小得多，可見其精度高，收斂速度快。

圖9 電勢的相對誤差隨網(wǎng)格尺寸變化圖

整個域上的電勢（單位:V）分布圖如圖10（見212頁）所示，圖10（a）與（b）分別是IGA的數(shù)值解及其和解析解的相對誤差分布圖，誤差很小，只有10-7數(shù)量級，可見其求解精度高。

從偏心圓柱面靜電場問題可以看出，一方面等幾何分析方法天然地具有等幾何性，即計算模型與設(shè)計模型一致，且還擁有細化保形性；另一方面等幾何分析方法在計算方面還具有單位單元的計算自由度消耗小，計算精度高，收斂速度快的特點。

3.2 偏心橢圓柱面

偏心橢圓柱面與偏心圓柱面類似，只是帶電體是橢圓柱面，如圖11所示。

圖11 偏心橢圓柱面計算簡圖

在仿真中，取長半軸a1=1.0m，a2=0.5m，b1／a1=b2／a2=0.5，L=0.2m，V1=0V，V2=1V.

柱面形狀的改變使得解析法很難進行。在3.1節(jié)中的收斂性分析的基礎(chǔ)上，直接給出了等幾何分析的求解結(jié)果分布圖，如圖12（見212頁）所示。從圖中可以看出，電勢在偏心方向比較集中，且此處的電場強度較高。

4.結(jié) 論

將等幾何分析方法擴展到靜電場問題，推導(dǎo)了靜電場問題的等幾何分析方法求解格式。該方法具有細化保形性，計算自由度少，收斂速度快，計算精度和計算效率高的特點。并應(yīng)用此方法求解了偏心圓柱面和偏心橢圓柱面的靜電場問題，結(jié)果表明了該方法具有上述的優(yōu)點，適合進一步在計算電磁學(xué)中推廣。

［1］李 敬.級數(shù)法求解軸對稱的靜電場［J］.云南師范大學(xué)學(xué)報，1997，17（3）:43-45.LI Jing.Inprogression resolution statics field with symmetric axis［J］.Journal of Yunnan Normal University:Natural Sciences Edition，1997，17（3）:43-45.（in Chinese）

［2］馬西奎.最小二乘邊界配點法在電磁場邊值問題數(shù)值分析中的應(yīng)用［J］.微波學(xué)報，1994，37（2）:17-22.MA Xikui.Application of least square boundary point matching method to the numerical analysis of electromagnetic field problems［J］.Journal of Microwares，1994，37（2）:17-22.（in Chinese）

［3］BALLISTI R，HAFNER C.The multiple multipole method in electro-and magnetostatic problems［J］.IEEE Transactions on Magnetics，1983，19（6）:2367-2370.

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

［5］周 勤，茅乃豐.旋轉(zhuǎn)對稱靜電場的一種新的數(shù)值解法［J］.電波科學(xué)學(xué)報，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）

［6］SONG B，F(xiàn)U 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.

［7］金建銘.電磁場有限元方法［M］.西安:西安電子科技大學(xué)出版社，1998:51-70.

［8］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.

［9］甄蜀春，曹 蕾，張繼龍.波動方程差分解法對波導(dǎo)螺釘調(diào)配器的分析［J］.電波科學(xué)學(xué)報，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）

［10］汪朝暉，廖振方，陳德淑.有限元法分析尖板電極結(jié)構(gòu)的空間靜電場分布［J］.重慶大學(xué)學(xué)報，2010，33（5）:41-47.WANG Zhaohui，LIAO Zhenfang，CHEN Deshu.A-nalysis 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）

［11］梁志偉，趙國偉，徐 杰，等.柱形等離子體天線輻射特性的矩量法分析［J］.電波科學(xué)學(xué)報，2008，23（4）:749-753 LIANG Zhiwei，ZHAO Guowei，XU Jie，et al.A-nalysis of plasma-column antenna using moment method［J］.Chinese Journal of Radio Science，2008，23（4）:749-753.（in Chinese）

［12］張淮清.電磁場計算中的徑向基函數(shù)無網(wǎng)格法研究［D］.重慶:重慶大學(xué)，2008.ZHANG Huaiqing.Research on Radial Basis Function Meshless Method in Numercial Computation of Electromagentic Field［D］.Chongqing:Chongqing U-niversity，2008.（in Chinese）

［13］HUGHES T J R，COTTRELL J A，BAZILEVS Y.Isogeometric analysis CAD， finite elements，NURBS，exact geometry and mesh refinement［J］.Comput.Methods Appl.Mech.Engrg.，2005，194（39／40／41）:4135-4195.

［14］REALI A.An isogeometric analysis approach for the study of structural vibrations［J］.Journal of Earthquake Engineering，2006，10（Special Issue 1）:1-30.

［15］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.

［16］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／42／43）:5257-5296.

［17］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／42／43／44）:4160-4183.

［18］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.

［19］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.

［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.

［22］陳云林，李 旭，蘭明建.偏心圓柱面靜電場電勢分布和電場強度分布的求解［J］.渝西學(xué)院學(xué)報，2002，1（4）:28-31.CHEN Yunlin，LI Xu，LAN Mingjian.The solution for the electrical potential and electrical field intensity s distribution inside two parallel conductor cylinders whose axes are separated［J］.Journal of Western Chongqing University，2002，1（4）:28-31.（in Chinese）