高階曲線矢量有限元方法實現及關鍵問題

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

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

1.引 言

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

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

2.高階曲線矢量元的實現公式

如圖1所示，設真實的曲線不規則圖形所在的xyz坐標系為全局坐標系，而映射后標準規則圖形所在的uvw坐標系為局部坐標系，即原點為（0，0，0）的標準笛卡兒坐標系。

則曲線矢量元的實現過程分為以下幾個步驟：

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

2）選擇合適階次的曲線映射函數，推導對應的曲線映射Jacobian矩陣。

3）選擇合適的矢量基函數，推導該矢量元及對應旋度公式在局部坐標系下的表達形式。

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

5）組合單元矩陣并強加邊界條件，得到整體矩陣并求解矩陣方程，對計算的數據后處理。

與直線形式的矢量元相比[1]，步驟①～④是曲線矢量元實現時所獨有的。

2.1 三維協變映射

定義逆協變映射為F∶（x，y，z）→（u，v，w），對應的Jacobian矩陣J（F）為

則協變映射為F－1∶（u，v，w）→（x，y，z）。此時，

由JJ－1＝I，則det（J－1）＝1/det（J）。顯然，在uvw坐標系計算偏導要比在xyz坐標系更容易，因為在該坐標系下曲線映射函數已知，可得到精確解析值；而在xyz坐標系下，網格單元通常是非規則的，對應積分非常繁瑣。

設r為xyz坐標系從原點到任意點P（x，y，z）的矢量，由全微分的概念[15]，

則uvw坐標系單位矢量定義為[16]

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

單位矢量和互易單位矢量滿足雙正交關系aiaj＝δij.

應用曲線坐標投影關系，任意矢量場W可用協變分量或逆協變分量表示[17]：

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

考慮到協變分量的定義，則協變分量之間的映射關系為

根據偏微分理論[18]，任意函數f的偏導分量之間的映射關系為

2.2 局部坐標系下曲線單元矩陣積分公式

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

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

同理，定義xyz坐標系下▽×W的逆協變分量

進一步推導可得兩者之間存在以下關系

顯然，與式（14）相比，文獻[17] P119式（42）、文獻[14] P441式（14）對應的公式都是不對的，缺少這一項。

據此，矢量有限元方法中單元矩陣所包含的兩個積分公式可轉換為uvw坐標系下的積分形式。

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

則單元矩陣

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

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

2.3 曲線映射函數選取及Jacobian矩陣確定

為了將xyz坐標系中的曲線單元映射為uvw坐標系下的規則直線單元，必須要定義相互之間的曲線映射法則。采用二階Lagrangian多項式函數作為曲線映射函數，實現對曲邊界的二階近似。

對于n階四面體結點單元，曲線映射函數定義為[4]

因此，二階曲線映射函數的具體表達式如下[4]

為了獲得正確的曲線映射關系，要使單元結點輸出格式滿足正確的對應關系，如圖1所示。

則xyz坐標系下的任意結點坐標采用uvw坐標系可表示為

式中 （xi，yi，zi）為xyz坐標系下的結點坐標。根據Jacobian矩陣定義式（1），則

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

2.4 矢量基函數及相關旋度的局部坐標表示

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

第i個面元基函數可表示為

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

如圖1所示，規則四面體棱邊長度為1時，u、v、w分別表示重心坐標Li2、Li3、Li4，即

由梯度在曲線坐標系中的定義[16]，

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

以e1的第一個分量（即對應棱邊為i1－i2）為例，代入式（25）、式（27）可得

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

采用相同的方法，對Webb99疊層基的其它基函數進行局部坐標表示，可得到基函數的協變分量、基函數旋度的逆協變分量（詳細結果可參考文獻[20] ，在此不再贅述）。需要指出的是，與文獻[14] 推導結果不同的原因在于，文中采用的基函數對面元形式進行了歸一化，同時f1、f2面元分量分別按（i1，i2，i3）、（i3，i1，i2）結點順序進行選取，從而確保了面元選擇的唯一性，便于程序實現與單元矩陣組合。

3.高階曲線矢量元實現流程及關鍵問題

圖2給出了采用固定階曲線矢量元的程序實現流程。具體實現時，先通過合適的判別準則，得到曲線單元編號的數組；然后在計算單元矩陣時分別采用合適的積分公式進行處理。當單元為曲線單元時，采用高斯數值積分進行計算；當單元為直線單元時，則采用高斯解析公式進行計算。與直線型矢量元的實現過程相比[20]，除了要增加圖2中虛線框所示的內容外，其它工作完全相同。

圖2 二階四面體曲線矢量元的實現流程

在混合階矢量元的實現過程中，需要解決以下幾個關鍵問題：

3.1 曲線矢量單元的判別

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

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

3.2 曲線單元矩陣的數值積分實現

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

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

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

具體而言，曲線單元矩陣高斯數值積分的實現流程如圖3所示。

圖3 曲線單元矩陣高斯數值積分的實現流程

當采用局部坐標分量表示單元矩陣積分點時，則

顯然，積分點與L1i無關，因為L1i＝1－u－vw.當規則直四面體的棱邊長度為1時，高斯積分對應體積V＝1/6。

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

4.數值結果

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

4.1 驗證結果

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

圖4繪出了歸一化均方誤差ε與網格尺寸h之間的關系。可以看出，分析曲線模型的精度主要受限于對模型曲邊界的近似精度。在相同尺寸的前提下，若采用直線元分析曲線區域，高階元只能達到低階元的精度；若采用曲線元分析曲線區域，二階元的計算精度遠高于棱邊元。值得一提的是，曲線二階矢量元的收斂速度為O（h－2.07），而直線二階矢量元分析直線模型時可達到 O（h－4）[1]。這是因為，二階曲線映射函數雖比直線元更能提高近似曲線邊界的精度，但不能準確模擬球邊界，必須要提高曲線映射函數的階次或采用Bezier四面體[24]等特殊曲線映射單元。

圖4 球腔相對誤差與平均棱邊尺寸的關系

圖5繪出了歸一化均方誤差ε與未知量數目之間的關系。可以看出，當分析曲線目標時，在相同未知量數目的前提下，采用高階曲線矢量元精度遠優于采用直線形式的高、低階元或曲線單元形式的低階元；而高階直線元雖然在網格尺寸相同的情況下獲得與低階直線元相比擬的計算精度，但未知量數目急劇增加。這充分驗證了，只能先采用曲線元對曲線模型進行高精度建模，高階矢量基函數才能發揮以較少未知量獲得較高精度的優點。

圖5 球腔相對誤差與未知量數目的關系

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

圖6 矩陣條件數與未知量數目的關系

4.2 具體應用

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

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

圖9繪出了對Lebaric介質填充圓柱腔計算的主模諧振頻率歸一化均方誤差ε與未知量數目之間的關系。可以看出，高階曲線矢量元由于同時具有對曲邊界模擬程度好、對場近似程度高的特點，因此獲得了比其它形式矢量元更高的計算精度。

圖9 Lebaric圓柱腔主模相對誤差與未知量數目的關系

5.結 論

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

[1] 尹文祿，鄧 聰，柴舜連，等.基于Nedelec條件的高階矢量元構造與實現 [J] .微波學報，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] 尹文祿，葉良豐，鄧 聰，等.基于高階四面體矢量元的大規模本征值求解 [J] .微波學報，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] 尹文祿，鄧 聰，楊 虎，等.高階四面體矢量元的實現與性能比較 [J] .微波學報，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，FERRARI 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] 汪 浩.高等數學 （下冊）[M] .長沙：國防科技大學出版社，1988.

[19] 尹文祿，鄧 聰，趙 菲，等.高階矢量有限元方法實現及關鍵問題 [J] .電波科學學報，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] 尹文祿.高階矢量有限元方法在電磁領域中的研究及應用 [D] .長沙：國防科技大學，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.