張 勇 林 皋 劉 俊 徐喜榮

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

引 言

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

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

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

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

1.B樣條和NURBS

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

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

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

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

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

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

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

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

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

NURBS曲線用向量形式表示為

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

對(duì)于NURBS曲面也有相似的形式，且為

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

2.1 靜電場(chǎng)問(wèn)題

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

2.2 求解域的等幾何離散

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

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

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

圖4 NURBS離散求解域

2.3 求解域的保形細(xì)化

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

2.4 電勢(shì)場(chǎng)的等參化

等參方法用于近似求解域的電勢(shì)場(chǎng)，求解域中的任一點(diǎn)的電勢(shì)具有與式（7）相似的表示形式。

式中:φi表示求解域控制點(diǎn)的電勢(shì)場(chǎng)變量。

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

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

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

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

Q=是等效控制點(diǎn)激勵(lì)，包括體電荷分布和第二類邊界條件的貢獻(xiàn)，表達(dá)式為

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

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

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

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

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

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

3.偏心柱面靜電場(chǎng)問(wèn)題算例分析

3.1 偏心圓柱面

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

上述問(wèn)題簡(jiǎn)化為平面問(wèn)題，如圖5所示，其電勢(shì)解析解為

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

圖5 偏心圓柱面示意圖

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

圖6 偏心圓柱面計(jì)算簡(jiǎn)圖

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

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

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

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

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

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

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

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

3.2 偏心橢圓柱面

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

圖11 偏心橢圓柱面計(jì)算簡(jiǎn)圖

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

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

4.結(jié) 論

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

［1］李 敬.級(jí)數(shù)法求解軸對(duì)稱的靜電場(chǎng)［J］.云南師范大學(xué)學(xué)報(bào)，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］馬西奎.最小二乘邊界配點(diǎn)法在電磁場(chǎng)邊值問(wèn)題數(shù)值分析中的應(yīng)用［J］.微波學(xué)報(bào)，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］盛劍霓.電磁場(chǎng)與波分析中半解析法的理論方法與應(yīng)用［M］.北京:科學(xué)出版社，2006:1-35.

［5］周 勤，茅乃豐.旋轉(zhuǎn)對(duì)稱靜電場(chǎng)的一種新的數(shù)值解法［J］.電波科學(xué)學(xué)報(bào)，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］金建銘.電磁場(chǎng)有限元方法［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òng)方程差分解法對(duì)波導(dǎo)螺釘調(diào)配器的分析［J］.電波科學(xué)學(xué)報(bào)，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)的空間靜電場(chǎng)分布［J］.重慶大學(xué)學(xué)報(bào)，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］梁志偉，趙國(guó)偉，徐 杰，等.柱形等離子體天線輻射特性的矩量法分析［J］.電波科學(xué)學(xué)報(bào)，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］張淮清.電磁場(chǎng)計(jì)算中的徑向基函數(shù)無(wú)網(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］陳云林，李 旭，蘭明建.偏心圓柱面靜電場(chǎng)電勢(shì)分布和電場(chǎng)強(qiáng)度分布的求解［J］.渝西學(xué)院學(xué)報(bào)，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）