尹文祿 楊 虎 肖 科 柴舜連 毛鈞杰

（1.西南電子電信技術(shù)研究所，四川 成都610041；2.國防科技大學(xué)電子科學(xué)與工程學(xué)院，湖南 長沙410073）

1.引 言

有限元方法由于具有處理不均勻媒質(zhì)和復(fù)雜幾何體的能力，在電磁計(jì)算領(lǐng)域得到了廣泛的研究與應(yīng)用。如何構(gòu)造基函數(shù)及提高其收斂速度，是有限元發(fā)展的一個(gè)重要方向。從20世紀(jì)80年代末起，各種高階矢量元由于具有以較少的未知量獲得較高精度的優(yōu)點(diǎn)[1]，成為研究的熱點(diǎn)方向之一[2－3]。其中，基于重心坐標(biāo)的三角形、四面體單元能對(duì)任意復(fù)雜幾何體進(jìn)行精確、靈活的建模，且在未知量相同的情況下，計(jì)算精度遠(yuǎn)高于其它形式的基函數(shù)（如四邊形、六面體等）[4]，因此應(yīng)用最為廣泛。針對(duì) Nedelec隱式的構(gòu)造方法[5]，各種四面體矢量基函數(shù)層出不窮，可分為兩類：插值基[6]和疊層基[7]，文獻(xiàn)[1] 、[8] －[9] 對(duì)各種形式的高階四面體矢量元進(jìn)行了系統(tǒng)構(gòu)建、實(shí)現(xiàn)與性能比較，驗(yàn)證了滿足相同函數(shù)空間的插值（或疊層）矢量基函數(shù)具有相同計(jì)算精度。

有限元法分析電磁問題存在兩個(gè)近似：一是通過網(wǎng)格剖分對(duì)目標(biāo)模型進(jìn)行區(qū)域離散引起的模型誤差，二是通過選取基函數(shù)來近似場分布帶來的計(jì)算誤差，問題求解精度取決于這兩個(gè)近似誤差的疊加效果。對(duì)于后一個(gè)近似，可采用各種高階插值或疊層矢量基函數(shù)來提高場近似的精度；而對(duì)于前一個(gè)近似，當(dāng)采用直線單元分析曲線邊界時(shí)，對(duì)幾何模型的近似程度決定了求解精度。隨著高階單元剖分尺寸增加，線性規(guī)則單元對(duì)曲邊界不能很好地模擬。為了克服該缺點(diǎn)，各種曲線元被引入到高階矢量元中，通過協(xié)變映射，將曲線單元映射為直線單元進(jìn)行計(jì)算。目前，針對(duì)曲線矢量元實(shí)現(xiàn)的文獻(xiàn)主要集中在四邊形和六面體形式[10－12]，分析三角形和四面體的文獻(xiàn)較少[13－14]。針對(duì)這個(gè)問題，在文獻(xiàn)[14] 的基礎(chǔ)上，進(jìn)一步完善高階四面體曲線矢量元的實(shí)現(xiàn)過程，并討論了其實(shí)現(xiàn)流程及一些關(guān)鍵問題，最后通過數(shù)值實(shí)例驗(yàn)證了曲線矢量元極好的計(jì)算性能。

2.高階曲線矢量元的實(shí)現(xiàn)公式

如圖1所示，設(shè)真實(shí)的曲線不規(guī)則圖形所在的xyz坐標(biāo)系為全局坐標(biāo)系，而映射后標(biāo)準(zhǔn)規(guī)則圖形所在的uvw坐標(biāo)系為局部坐標(biāo)系，即原點(diǎn)為（0，0，0）的標(biāo)準(zhǔn)笛卡兒坐標(biāo)系。

則曲線矢量元的實(shí)現(xiàn)過程分為以下幾個(gè)步驟：

1）研究曲線坐標(biāo)變換的相關(guān)公式，得到局部坐標(biāo)系下兩種類型單元矩陣的顯式積分公式。

2）選擇合適階次的曲線映射函數(shù)，推導(dǎo)對(duì)應(yīng)的曲線映射Jacobian矩陣。

3）選擇合適的矢量基函數(shù)，推導(dǎo)該矢量元及對(duì)應(yīng)旋度公式在局部坐標(biāo)系下的表達(dá)形式。

4）選擇合適的積分法則，數(shù)值求解曲線矢量元的單元矩陣。

5）組合單元矩陣并強(qiáng)加邊界條件，得到整體矩陣并求解矩陣方程，對(duì)計(jì)算的數(shù)據(jù)后處理。

與直線形式的矢量元相比[1]，步驟①～④是曲線矢量元實(shí)現(xiàn)時(shí)所獨(dú)有的。

2.1 三維協(xié)變映射

定義逆協(xié)變映射為F∶（x，y，z）→（u，v，w），對(duì)應(yīng)的Jacobian矩陣J（F）為

則協(xié)變映射為F－1∶（u，v，w）→（x，y，z）。此時(shí)，

由JJ－1＝I，則det（J－1）＝1/det（J）。顯然，在uvw坐標(biāo)系計(jì)算偏導(dǎo)要比在xyz坐標(biāo)系更容易，因?yàn)樵谠撟鴺?biāo)系下曲線映射函數(shù)已知，可得到精確解析值；而在xyz坐標(biāo)系下，網(wǎng)格單元通常是非規(guī)則的，對(duì)應(yīng)積分非常繁瑣。

設(shè)r為xyz坐標(biāo)系從原點(diǎn)到任意點(diǎn)P（x，y，z）的矢量，由全微分的概念[15]，

則uvw坐標(biāo)系單位矢量定義為[16]

從而，互易單位矢量可定義為[16]

單位矢量和互易單位矢量滿足雙正交關(guān)系aiaj＝δij.

應(yīng)用曲線坐標(biāo)投影關(guān)系，任意矢量場W可用協(xié)變分量或逆協(xié)變分量表示[17]：

式中：Wu，v，w＝W·au，v，w為協(xié)變分量；Wu，v，w＝W·au，v，w為逆協(xié)變分量。顯然，協(xié)變分量表示矢量場W的切向分量，逆協(xié)變分量表示W(wǎng)的法向分量。

考慮到協(xié)變分量的定義，則協(xié)變分量之間的映射關(guān)系為

根據(jù)偏微分理論[18]，任意函數(shù)f的偏導(dǎo)分量之間的映射關(guān)系為

2.2 局部坐標(biāo)系下曲線單元矩陣積分公式

將兩種坐標(biāo)系中各分量的映射關(guān)系式（8）至式（10）代入xyz坐標(biāo)系下的旋度場可得

比較逆協(xié)變分量公式（7）與式（11），則可定義uvw坐標(biāo)系下▽×W的逆協(xié)變分量

同理，定義xyz坐標(biāo)系下▽×W的逆協(xié)變分量

進(jìn)一步推導(dǎo)可得兩者之間存在以下關(guān)系

顯然，與式（14）相比，文獻(xiàn)[17] P119式（42）、文獻(xiàn)[14] P441式（14）對(duì)應(yīng)的公式都是不對(duì)的，缺少這一項(xiàng)。

據(jù)此，矢量有限元方法中單元矩陣所包含的兩個(gè)積分公式可轉(zhuǎn)換為uvw坐標(biāo)系下的積分形式。

由三重積分的一般變量替換公式[18]

則單元矩陣

為了公式的簡單起見，式（16）將Vuvw縮寫為V.

對(duì)于早期大部分文獻(xiàn)而言，公式推導(dǎo)都是假定單元結(jié)點(diǎn)順序滿足右手系關(guān)系，此時(shí)，det（J）恒為正。為了滿足單元交界面處棱邊方向的一致性（均從較小點(diǎn)指向較大點(diǎn)），采用文獻(xiàn)[1] 提出的結(jié)點(diǎn)編碼方法，即規(guī)定結(jié)點(diǎn)中的4個(gè)頂點(diǎn)（i1，i2，i3，i4）輸出滿足結(jié)點(diǎn)遞增順序（設(shè)i1＜i2＜i3＜i4），顯然這將不滿足右手系關(guān)系，因此J的行列式要加絕對(duì)值。

2.3 曲線映射函數(shù)選取及Jacobian矩陣確定

為了將xyz坐標(biāo)系中的曲線單元映射為uvw坐標(biāo)系下的規(guī)則直線單元，必須要定義相互之間的曲線映射法則。采用二階Lagrangian多項(xiàng)式函數(shù)作為曲線映射函數(shù)，實(shí)現(xiàn)對(duì)曲邊界的二階近似。

對(duì)于n階四面體結(jié)點(diǎn)單元，曲線映射函數(shù)定義為[4]

因此，二階曲線映射函數(shù)的具體表達(dá)式如下[4]

為了獲得正確的曲線映射關(guān)系，要使單元結(jié)點(diǎn)輸出格式滿足正確的對(duì)應(yīng)關(guān)系，如圖1所示。

則xyz坐標(biāo)系下的任意結(jié)點(diǎn)坐標(biāo)采用uvw坐標(biāo)系可表示為

式中 （xi，yi，zi）為xyz坐標(biāo)系下的結(jié)點(diǎn)坐標(biāo)。根據(jù)Jacobian矩陣定義式（1），則

因此，只要計(jì)算出曲線映射函數(shù)對(duì)u、v、w的偏導(dǎo)數(shù)，則J矩陣就很容易計(jì)算得到。求出J后，由于J矩陣僅為3×3階，因此可采用直接方法求解單元矩陣計(jì)算時(shí)所需的J－1、J－1T、det（J）。由于J矩陣與數(shù)值積分點(diǎn)的位置相關(guān)，因此在計(jì)算單元矩陣時(shí)必須在每一個(gè)積分點(diǎn)處都要進(jìn)行計(jì)算。

2.4 矢量基函數(shù)及相關(guān)旋度的局部坐標(biāo)表示

由文獻(xiàn)[1] [9] ，只要滿足相同的Nedelec函數(shù)空間，則基函數(shù)可獲得相同的計(jì)算精度。因此，選取二階非歸一化Webb99疊層基[7]作為單元內(nèi)場近似基函數(shù)。則第i條棱邊上的邊元基函數(shù)可表示為

第i個(gè)面元基函數(shù)可表示為

式中：Li為重心坐標(biāo)。顯然，采用基函數(shù)分類技術(shù)，二階矢量基函數(shù)可分為四類：e1、e2、f1、f2。為了保證四面體單元之間切向場連續(xù)，采用文獻(xiàn)[1] [19] 定義的二階四面體矢量元的四類分量取法，即在每個(gè)面內(nèi)（i1，i2，i3）（i1＜i2＜i3），棱邊e1方向從較小點(diǎn)指向較大點(diǎn)，e2方向與之相反，f1在棱邊i1－i2上，取點(diǎn)順序?yàn)椋╥1，i2，i3）；f2在棱邊i1－i3上，取點(diǎn)順序?yàn)椋╥3，i1，i2）。

如圖1所示，規(guī)則四面體棱邊長度為1時(shí)，u、v、w分別表示重心坐標(biāo)Li2、Li3、Li4，即

由梯度在曲線坐標(biāo)系中的定義[16]，

則式（5）定義的互易單位矢量可改寫為

以e1的第一個(gè)分量（即對(duì)應(yīng)棱邊為i1－i2）為例，代入式（25）、式（27）可得

從式（28）可以看出，基函數(shù)的協(xié)變分量[Wu Wv Ww]1＝[1－v－wuu] 。由式（12），基函數(shù)旋度的逆協(xié)變分量 [Vu Vv Vw]1＝ [0－2 2] 。

采用相同的方法，對(duì)Webb99疊層基的其它基函數(shù)進(jìn)行局部坐標(biāo)表示，可得到基函數(shù)的協(xié)變分量、基函數(shù)旋度的逆協(xié)變分量（詳細(xì)結(jié)果可參考文獻(xiàn)[20] ，在此不再贅述）。需要指出的是，與文獻(xiàn)[14] 推導(dǎo)結(jié)果不同的原因在于，文中采用的基函數(shù)對(duì)面元形式進(jìn)行了歸一化，同時(shí)f1、f2面元分量分別按（i1，i2，i3）、（i3，i1，i2）結(jié)點(diǎn)順序進(jìn)行選取，從而確保了面元選擇的唯一性，便于程序?qū)崿F(xiàn)與單元矩陣組合。

3.高階曲線矢量元實(shí)現(xiàn)流程及關(guān)鍵問題

圖2給出了采用固定階曲線矢量元的程序?qū)崿F(xiàn)流程。具體實(shí)現(xiàn)時(shí)，先通過合適的判別準(zhǔn)則，得到曲線單元編號(hào)的數(shù)組；然后在計(jì)算單元矩陣時(shí)分別采用合適的積分公式進(jìn)行處理。當(dāng)單元為曲線單元時(shí)，采用高斯數(shù)值積分進(jìn)行計(jì)算；當(dāng)單元為直線單元時(shí)，則采用高斯解析公式進(jìn)行計(jì)算。與直線型矢量元的實(shí)現(xiàn)過程相比[20]，除了要增加圖2中虛線框所示的內(nèi)容外，其它工作完全相同。

圖2 二階四面體曲線矢量元的實(shí)現(xiàn)流程

在混合階矢量元的實(shí)現(xiàn)過程中，需要解決以下幾個(gè)關(guān)鍵問題：

3.1 曲線矢量單元的判別

選擇合適的判別準(zhǔn)則，對(duì)每個(gè)剖分單元進(jìn)行判斷，若滿足曲線單元條件（單元至少含有一條曲線即可），則取為曲線單元，否則為直線單元。

這是因?yàn)?，采用高階數(shù)值積分法則進(jìn)行曲線單元積分時(shí)，極大增加了計(jì)算量。因此，為了降低計(jì)算成本，只在含曲邊的區(qū)域（如曲線邊界、不同媒質(zhì)之間的曲交界面處等）上采用曲線單元，而在其它位置采用直線單元。由于曲線矢量元滿足了相鄰單元切向場分量連續(xù)（該條件即使在曲邊上也是滿足的），因此在不同單元形式的交界面處進(jìn)行正常組合即可得到整體矩陣。

3.2 曲線單元矩陣的數(shù)值積分實(shí)現(xiàn)

在計(jì)算單元矩陣時(shí)，對(duì)每一個(gè)單元編號(hào)，先判斷屬于曲線元還是直線元，直線單元采用二階矢量元高斯解析積分公式[8]進(jìn)行計(jì)算。對(duì)于曲線單元而言，由于曲線映射函數(shù)Ni、基函數(shù)Wi都是u、v、w的函數(shù)，這導(dǎo)致了單元矩陣中的待積分系數(shù)（如J、J－1、V、W 等）都是函數(shù)矩陣，使得很難求出積分的解析表達(dá)式，從而只能通過數(shù)值積分來求解曲線單元矩陣。在有限元法中，由于被積函數(shù)復(fù)雜，一般采用高斯數(shù)值積分，即可用較少的積分點(diǎn)達(dá)到較高的計(jì)算精度，從而節(jié)省計(jì)算時(shí)間。

基于四面體單元的三維高斯數(shù)值積分公式為[4]

式中：F（L1i，L2i，L3i，L4i）表示被積函數(shù)；（L1i，L2i，L3i，L4i）表示積分點(diǎn)；Wi為加權(quán)因子；ns為積分點(diǎn)數(shù)；V為積分單元的體積。該積分公式一個(gè)重要性質(zhì)在于[21]，如果被積函數(shù)多項(xiàng)式的階數(shù)≤（1，2，3，4，5），那么，選擇ns＝（1，4，5，11，14）則可精確求出積分。

具體而言，曲線單元矩陣高斯數(shù)值積分的實(shí)現(xiàn)流程如圖3所示。

圖3 曲線單元矩陣高斯數(shù)值積分的實(shí)現(xiàn)流程

當(dāng)采用局部坐標(biāo)分量表示單元矩陣積分點(diǎn)時(shí)，則

顯然，積分點(diǎn)與L1i無關(guān)，因?yàn)長1i＝1－u－vw.當(dāng)規(guī)則直四面體的棱邊長度為1時(shí)，高斯積分對(duì)應(yīng)體積V＝1/6。

本文研究表明，高斯積分點(diǎn)的數(shù)目與曲線映射函數(shù)的階次無關(guān)，僅與基函數(shù)階次有關(guān)。當(dāng)采用n階矢量基函數(shù)時(shí)，可獲得2n階計(jì)算精度，則對(duì)應(yīng)高斯積分只要滿足2n階精度即可。如棱邊元可獲得2階計(jì)算精度，因此采用2階4點(diǎn)Gellert積分[21]即可獲得精確結(jié)果；而二階矢量元可獲得4階精度，因此采用4階11點(diǎn)Gellert積分即可。當(dāng)然，采用更高階高斯積分同樣可獲得相同的計(jì)算精度，但無疑增加了計(jì)算代價(jià)。如采用6階24點(diǎn)或8階45點(diǎn)Keast積分[22]也可獲得與4階11點(diǎn)Gellert積分相同精度。需要指出的是，所有Keast積分與體積無關(guān)，因此被積函數(shù)不需乘以體積1/6；同時(shí)，8階45點(diǎn)Keast積分的首個(gè)積分點(diǎn)權(quán)值應(yīng)為－0.039327（即加一個(gè)負(fù)號(hào)），文獻(xiàn)[23] 糾正了該錯(cuò)誤。

4.數(shù)值結(jié)果

以下給出幾個(gè)算例的計(jì)算結(jié)果，分別從計(jì)算精度、矩陣條件數(shù)等角度驗(yàn)證曲線矢量元較之直線單元具有更好的計(jì)算性能。其中，直線單元分別采用一階非歸一化棱邊元、二階非歸一化Webb99疊層基；曲線單元在采用二階曲線映射函數(shù)的基礎(chǔ)上，基函數(shù)分別采用一階非歸一化棱邊元、二階非歸一化Webb99疊層基。

4.1 驗(yàn)證結(jié)果

采用與文獻(xiàn)[1] 相同的數(shù)值實(shí)驗(yàn)方法分析r＝1 m空氣填充球腔的諧振本征值。在不同網(wǎng)格尺寸h的情況下，計(jì)算前八個(gè)非零本征值與解析解之間的歸一化均方誤差ε.以整體矩陣T為例，計(jì)算不同矢量元對(duì)應(yīng)的矩陣條件數(shù)。其中，含邊界棱邊的單元采用曲線元，其余單元采用直線元。

圖4繪出了歸一化均方誤差ε與網(wǎng)格尺寸h之間的關(guān)系?？梢钥闯觯治銮€模型的精度主要受限于對(duì)模型曲邊界的近似精度。在相同尺寸的前提下，若采用直線元分析曲線區(qū)域，高階元只能達(dá)到低階元的精度；若采用曲線元分析曲線區(qū)域，二階元的計(jì)算精度遠(yuǎn)高于棱邊元。值得一提的是，曲線二階矢量元的收斂速度為O（h－2.07），而直線二階矢量元分析直線模型時(shí)可達(dá)到 O（h－4）[1]。這是因?yàn)椋A曲線映射函數(shù)雖比直線元更能提高近似曲線邊界的精度，但不能準(zhǔn)確模擬球邊界，必須要提高曲線映射函數(shù)的階次或采用Bezier四面體[24]等特殊曲線映射單元。

圖4 球腔相對(duì)誤差與平均棱邊尺寸的關(guān)系

圖5繪出了歸一化均方誤差ε與未知量數(shù)目之間的關(guān)系?？梢钥闯?，當(dāng)分析曲線目標(biāo)時(shí)，在相同未知量數(shù)目的前提下，采用高階曲線矢量元精度遠(yuǎn)優(yōu)于采用直線形式的高、低階元或曲線單元形式的低階元；而高階直線元雖然在網(wǎng)格尺寸相同的情況下獲得與低階直線元相比擬的計(jì)算精度，但未知量數(shù)目急劇增加。這充分驗(yàn)證了，只能先采用曲線元對(duì)曲線模型進(jìn)行高精度建模，高階矢量基函數(shù)才能發(fā)揮以較少未知量獲得較高精度的優(yōu)點(diǎn)。

圖5 球腔相對(duì)誤差與未知量數(shù)目的關(guān)系

圖6給出了在不同未知量數(shù)目情況下，采用直線、曲線矢量單元得到的整體矩陣T條件數(shù)。可以看出，不論采用曲線元還是直線元形式，只要基函數(shù)相同，則矩陣條件數(shù)幾乎完全相同 （圖中相同類型基函數(shù)的條件數(shù)曲線完全重合）。這意味著，在極大提高了計(jì)算精度的同時(shí)，曲線元形式并不影響計(jì)算效率。

圖6 矩陣條件數(shù)與未知量數(shù)目的關(guān)系

4.2 具體應(yīng)用

為了進(jìn)一步說明曲線矢量元在提高計(jì)算精度的能力，采用曲線矢量元分析圖7所示的圓柱腔，在曲邊界、內(nèi)部填充媒質(zhì)的曲交界面處采用曲線元，其余部分采用直線元。

如圖8所示，采用各種直線、曲線形式的基函數(shù)對(duì)該腔體進(jìn)行分析，得到不同未知量數(shù)目情況下的主模諧振頻率的計(jì)算結(jié)果。其中，測量結(jié)果為6.943 GHz[25].

圖9繪出了對(duì)Lebaric介質(zhì)填充圓柱腔計(jì)算的主模諧振頻率歸一化均方誤差ε與未知量數(shù)目之間的關(guān)系?？梢钥闯觯唠A曲線矢量元由于同時(shí)具有對(duì)曲邊界模擬程度好、對(duì)場近似程度高的特點(diǎn)，因此獲得了比其它形式矢量元更高的計(jì)算精度。

圖9 Lebaric圓柱腔主模相對(duì)誤差與未知量數(shù)目的關(guān)系

5.結(jié) 論

系統(tǒng)研究了曲線矢量元的實(shí)現(xiàn)過程，詳細(xì)分析了實(shí)現(xiàn)過程中的關(guān)鍵問題及相關(guān)細(xì)節(jié)。具體而言，以二階曲線映射函數(shù)結(jié)合二階非歸一化Webb99疊層基為例，系統(tǒng)而顯式地研究了曲線矢量元的實(shí)現(xiàn)流程和關(guān)鍵問題 （如協(xié)變映射的概念、曲線單元矩陣的公式推導(dǎo)、曲線映射函數(shù)選取與Jacobian矩陣求解、矢量基函數(shù)及其旋度在局部坐標(biāo)系中的展開、曲線單元矩陣的高斯數(shù)值積分技術(shù)等）。然后，通過分析r＝1 m球腔諧振本征值問題，從計(jì)算精度、條件數(shù)等方面對(duì)各種直線或曲線形式的高、低階矢量元進(jìn)行比較，得到了一些有益的結(jié)論。最后，將其用于分析復(fù)雜腔體問題，獲得了極好的計(jì)算結(jié)果。上述分析方法和結(jié)論可直接推廣到任意形式、更高階曲線矢量元的研究與應(yīng)用。

[1] 尹文祿，鄧 聰，柴舜連，等.基于Nedelec條件的高階矢量元構(gòu)造與實(shí)現(xiàn) [J] .微波學(xué)報(bào)，2009，25（3）：7－12.YIN Wenlu，DENG Cong，CHAI Shunlian，et al.Construction and implementation of higher order vector elements based on Nedelec constraints [J] .Journal of Microwaves，2009，25 （3）：7－12.（in Chinese）

[2] NOTAROS B M.Higher order frequency－domain computational electromagnetics [J] .IEEE Transactions on Antennas and Propagation，2008，56 （8）：2251－2276.

[3] 尹文祿，葉良豐，鄧 聰，等.基于高階四面體矢量元的大規(guī)模本征值求解 [J] .微波學(xué)報(bào)，2010，26（1）：12－18.YIN Wenlu，YE Liangfeng，DENG Cong，et al.Efficient solution of large－scale eigenvalue problems based on higher order tetrahedral vector elements[J] .Journal of Microwaves，2010，26 （1）：12－18.（in Chinese）

[4] JIN J M.The Finite Element Method in Electromagnetics[M] .2nd ed.New York：John Wiley and Sons，2002.

[5] NEDELEC J C.Mixed finite elements in R3[J] .Numerische Mathematik，1980，35 （9）：315－341.

[6] GRAGLIA R D，WILTON D R，PETERSON A F.Higher order interpolatory vector bases for computational electromagnetics [J] .IEEE Transactions on Antennas and Propagation，1997，45 （3）：329－342.

[7] WEBB J P.Hierarchal vector basis functions of arbitrary order for triangular and tetrahedral finite elements [J] .IEEE Transactions on Antennas and Propagation，1999，47 （8）：1244－1253.

[8] 尹文祿，鄧 聰，楊 虎，等.高階四面體矢量元的實(shí)現(xiàn)與性能比較 [J] .微波學(xué)報(bào)，2010，26 （3）：15－20.YIN Wenlu，DENG Cong，YANG Hu，et al.Implementation and performance comparison for higher order tetrahedral vector elements [J] .Journal of Microwaves，2010，26 （3）：15－20.（in Chinese）

[9] YIN W L，DONG J，CHEN L，et al.Systematic construction and performance comparison for higher order hierarchical tetrahedral vector elements[C] //International Conference on Wireless Communications and Signal Processing.Suzhou China，2010：1－4.

[10] CROWLEY C W，SILVESTER P P，HURWITZ H.Covariant projection elements for 3D vector field problems[J] .IEEE Transactions on Magnetics，1988，24 （1）：397－400.

[11] WANG J S，IDA N.Curvilinear and higher order＇edge＇finite elements in electromagnetic field computation [J] .IEEE Transactions on Magnetics，1993，29 （2）：1491－1494.

[12] PALMA M S，SARKER T K，CASTILLO L E G，et al.Iterative and Self Adaptive Finite Elements in Electromagnetic Modelling [M] .Boston：Artech House，1998.

[13] MARAIS N，DAVIDSON D B.Numerical evaluation of hierarchical vector finite elements on curvilinear domains in 2－D [J] .IEEE Transactions on Antennas and Propagation，2006，54 （2）：734－738.

[14] SWARTZ J P，DAVIDSON D B.Curvilinear vector finite elements using a set of hierarchical basis functions [J] .IEEE Transactions on Antennas and Propagation，2007，55 （2）：440－446.

[15] SILVESTER P P，F(xiàn)ERRARI R L.Finite Elements for Electrical Engineers [M] .3rd ed.Cambridge：Cambridge University Press，1996.

[16] STRATTON J A.Electromagnetic Theory [M] .New York：McGraw－Hill，1941.

[17] ITOH T，PELOSI G，SILVESTER P P.Finite Element Software for Microwave Engineering [M] .New York：John Wiley ＆Sons，1996.

[18] 汪 浩.高等數(shù)學(xué) （下冊(cè)）[M] .長沙：國防科技大學(xué)出版社，1988.

[19] 尹文祿，鄧 聰，趙 菲，等.高階矢量有限元方法實(shí)現(xiàn)及關(guān)鍵問題 [J] .電波科學(xué)學(xué)報(bào)，2009，24（2）：349－353.YIN Wenlu，DENG Cong，ZHAO Fei，et al.Implementation for higher order vector FEM and some key issues [J] .Chinese Journal of Radio Science，2009，24 （2）：349－353.（in Chinese）

[20] 尹文祿.高階矢量有限元方法在電磁領(lǐng)域中的研究及應(yīng)用 [D] .長沙：國防科技大學(xué)，2010.YIN Wenlu.Study and Application on Higher Order Vector FEM in Electromagnetic Field [D] .Changsha，P.R.China：National University of Defense Technology，2010.

[21] GELLERT M，HARBORD R.Moderate degree cubature formulas for 3－D tetrahedral finite－element approximations [J] .Communications in Applied Numerical Methods.1991，7 （6）：487－495.

[22] KEAST P.Moderate－degree tetrahedral quadrature formulas [J] .Computer Methods in Applied Mechanics and Engineering，1986，55 （3）：339－348.

[23] SAVAGE J S，PETERSON A F.Quadrature rules for numerical integration over triangles and tetrahedra[J] .IEEE Antennas and Propagation Magazine，1996，38 （3）：100－102.

[24] MARTINI E，SELLERI S.Innovative class of curvilinear tetrahedral elements [J] .IEE Electronics Letters，2001，37 （9）：557－558.

[25] LEBARIC J E，KAJFEZ D.Analysis of dielectric resonator cavities using the finite integration technique [J] .IEEE Transactions on Microwave Theory and Techniques，1989，37 （11）：1740－1748.